missed cleavages combinations generation

Question

I am looking for an R based solution to the following problem:

Trypsin does not cleave perfectly in proteomics experiments. So starting from a protein sequence FASTA I want to generate all possible peptides which COULD be generated if trypsin was allowed n number of "missed cleavages"

In informatics terms I would describe this as skipping a split during a strsplit() function

Example to clarify:

We start out with the protein sequence: GRGKA which is usually split by strsplit() (digested by trypsin) after every K and R

"GRGKA" %>% strsplit(., "(?<=K|R)", perl = TRUE)

[[1]]
[1] "GR" "GK" "A" 

now imagine that trypsin did not cleave perfectly, for example missing always one K or R

"GRGK" "GKA"

would be generated in this case. Or if trypsin digested this protein really badly and would end up missing two cleavages we would end up with the sequence:

"GRGKA"

I would love to have a generalizable solution that allows to generate all combinations of sequences that arise with n number of missed cleavages.

If you need any more information please feel free to ask, sorry if I missed something!

jay.sf's requested example:

"FLGERTCTRNERQP"  

perfect digest:

"FLGER" "TCTR" "NER" "QP"   

one missed cleavage all combinations:

"FLGERTCTR" "TCTRNER" "NERQP"

two missed cleavages all combinations:

"FLGERTCTRNER" "QP" "FLGER" "TCTRNERQP" 

three missed cleavages (there are a total of three K and Rs so trypsin would miss all)

"FLGERTCTRNERQP"

For longer sequences I think this becomes a frameshift like problem: If we count the K and Rs in the sequence:

FLGERTCTRNERQPFLGERTCTRNERQP
FLGE1TCT2NE3QPFLGE4TCT5NE6QP
FLGER_TCTR_NER_QPFLGER_TCTR_NER_QP (perfect digest)
FLGERTCTR_NERQPFLGER_TCTRNER_QP (one missed cleavage and always the uneven R and Ks)
FLGER_TCTRNER_QPFLGERTCTR_NERQP (one missed cleavage and always the even R and Ks)

So I think for 1 missed cleavage we have two frames over the whole sequence and for two missed cleavages we would have three frames and so on....

Answer 1

Essentially we want to combine 1, 2, ..., n - 1 splits, and if the variance of split lengths is greater than one we have a second solution. Here is a function that, after the strsplit, first creates indices cb to later paste the respective splits together. I'm not quite sure what the maximum number of splits is, but this might be a start.

mccg <- function(x) {
  ## missed cleavages combinations generation ++++++++++++++++++++++++++++++++++
  gs <- \(n, k, length=n*k, right=FALSE) {
    ## helper, similar to `gl()` but w/o factorization
    if (right) {
      cbr <- gs(n, k, n, right=FALSE)
      rev(max(cbr) - cbr + 1L)
    } else {
      rep_len(rep.int(seq_len(n), rep.int(k, n)), length)
    }
  }
  ## split string after [K, R]
  spl_x <- el(strsplit(x, '(?<=K|R)', perl=TRUE))
  len <- length(spl_x)
  res <- if (len <= 1L) {
    ## if [K, R] at the end or missing
    x
  } else {
    ## nested list of combinations
    cb <- lapply(1:(len - 1L) |> {\(.) setNames(., . - 1L)}(), \(i, .len=len) {
      cb <- list(gs(.len, i, .len))
      if (i != 1L && var(tabulate(el(cb)))) {
        ## second solution, starting from right
        c(cb, list(gs(.len, i, .len, right=TRUE)))
      } else {
        cb
      }
    }) |> c(setNames(list(rep_len(1L, len)), len - 1L))
    ## paste according to sequences
    rapply(cb, \(x) vapply(unname(split(spl_x, x)), paste, collapse='', 
                           character(1)), 
           how='list')
  }
  return(res)
}

Usage

Two splits:

mccg(x='GRGKA')
# $`0`
# [1] "GR" "GK" "A" 
# 
# $`1`
# $`1`[[1]]
# [1] "GRGK" "A"   
# 
# $`1`[[2]]
# [1] "GR"  "GKA"
# 
# 
# $`2`
# [1] "GRGKA"

Three splits:

mccg(x='FLGERTCTRNERQP')
# $`0`
# [1] "FLGER" "TCTR"  "NER"   "QP"   
# 
# $`1`
# [1] "FLGERTCTR" "NERQP"    
# 
# $`2`
# $`2`[[1]]
# [1] "FLGERTCTRNER" "QP"          
# 
# $`2`[[2]]
# [1] "FLGER"     "TCTRNERQP"
# 
# 
# $`3`
# [1] "FLGERTCTRNERQP"

More splits:

mccg(x='FLGERTCTRNERQPFLGERTCTRNERQP')
# $`0`
# [1] "FLGER"   "TCTR"    "NER"     "QPFLGER" "TCTR"    "NER"     "QP"     
# 
# $`1`
# $`1`[[1]]
# [1] "FLGERTCTR"  "NERQPFLGER" "TCTRNER"    "QP"        
# 
# $`1`[[2]]
# [1] "FLGER"       "TCTRNER"     "QPFLGERTCTR" "NERQP"      
# 
# 
# $`2`
# $`2`[[1]]
# [1] "FLGERTCTRNER"   "QPFLGERTCTRNER" "QP"            
# 
# $`2`[[2]]
# [1] "FLGER"          "TCTRNERQPFLGER" "TCTRNERQP"     
# 
# 
# $`3`
# $`3`[[1]]
# [1] "FLGERTCTRNERQPFLGER" "TCTRNERQP"          
# 
# $`3`[[2]]
# [1] "FLGERTCTRNER"     "QPFLGERTCTRNERQP"
# 
# 
# $`4`
# $`4`[[1]]
# [1] "FLGERTCTRNERQPFLGERTCTR" "NERQP"                  
# 
# $`4`[[2]]
# [1] "FLGERTCTR"           "NERQPFLGERTCTRNERQP"
# 
# 
# $`5`
# $`5`[[1]]
# [1] "FLGERTCTRNERQPFLGERTCTRNER" "QP"                        
# 
# $`5`[[2]]
# [1] "FLGER"                   "TCTRNERQPFLGERTCTRNERQP"
# 
# 
# $`6`
# [1] "FLGERTCTRNERQPFLGERTCTRNERQP"

Simulation

Here are some toy data,

set.seed(42)
n <- 1e3
min_chr <- 5; max_chr <- 50
strings <- replicate(n, {
  repeat {
    smp <- sample(setdiff(LETTERS, c("B", "Z", "J", "O", "U", "X")),
                  sample(min_chr:max_chr, 1L), replace=TRUE)
    if (any(grepl('K|R', smp))) break
  }
  paste(smp, collapse='')
}) |> setNames(nm=_)

over which the execution of the function can be tested in a loop.

r <- lapply(strings, \(x) tryCatch(mccg(x), error=\(e) NA))

Answer 2

1) gsubfn is like gsub except the replacement can be a proto object containing a fun proto method which is a function whose first argument receives the proto object and which inherits a count (and other) built-ins. The match is input to the function as the second argument and is replaced by the output of the function. Here we replace each R or K with itself and if it is more than the k-th occurrence then with a semicolon as well. Then split by semicolon.

library(gsubfn)

# inputs
x <- c("GRGKA", "GRGKAK")  # vector of inputs
k <- 1 # no of splits to skip

p <- proto(fun = function(self, x) paste0(x, if (count > k) ";"))
strsplit(gsubfn("[RK]", p, x), ";")

giving

[[1]]
[1] "GRGK" "A"   

[[2]]
[1] "GRGK" "AK"  

2) Insert a semicolon after each R or K and then remove the first k. strsplit what is left. This uses base R only.

xx <- gsub("([RK])", "\\1;", x)
for(i in seq_len(k)) xx <- sub(";", "", xx)
strsplit(xx, ";")
# same output as above