3
$\begingroup$

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....

$\endgroup$
3
  • 1
    $\begingroup$ Say "FLGERTCTRNERQP" would be a candidate, it has three possible splits. Can you show all possible results for 1, 2, and 3 skipped splits? $\endgroup$
    – jay.sf
    Jul 25, 2023 at 16:33
  • $\begingroup$ 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 $\endgroup$ Jul 26, 2023 at 8:47
  • 1
    $\begingroup$ @KlemensFröhlich the answer below looks cool, you've enough rep to vote - definitely recommend it. $\endgroup$
    – M__
    Jul 27, 2023 at 12:06

2 Answers 2

2
$\begingroup$

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))
$\endgroup$
1
  • $\begingroup$ hi jay Thank you very much for looking into this! This already looks promising, but the main goal would not be to get a function that can go as deep as possible but rather that can be limited at any level that the user wants to choose. So basically it would be ideal to have a function that allows to set the missed cleavages to 0, 1, 2, etc and always receive all the combinations up to this point. Do you think your function could be modified to control the depth of allowed missed cleavages? Best Klemens $\endgroup$ Aug 4, 2023 at 11:09
0
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.