<p>Essentially we want to combine <em>1, 2, ..., n - 1</em> splits, and if the <code>var</code>iance of split lengths is greater than one we have a second solution. Here is a function that,  after the <code>strsplit</code>, first creates indices <code>cb</code> to later <code>paste</code> the respective splits together. I'm not quite sure what the maximum number of splits is, but this might be a start.</p>
<pre><code>mccg &lt;- function(x) {
  ## missed cleavages combinations generation ++++++++++++++++++++++++++++++++++
  gs &lt;- \(n, k, length=n*k, right=FALSE) {
    ## helper, similar to `gl()` but w/o factorization
    if (right) {
      cbr &lt;- 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 &lt;- el(strsplit(x, '(?&lt;=K|R)', perl=TRUE))
  len &lt;- length(spl_x)
  res &lt;- if (len &lt;= 1L) {
    ## if [K, R] at the end or missing
    x
  } else {
    ## nested list of combinations
    cb &lt;- lapply(1:(len - 1L) |&gt; {\(.) setNames(., . - 1L)}(), \(i, .len=len) {
      cb &lt;- list(gs(.len, i, .len))
      if (i != 1L &amp;&amp; var(tabulate(el(cb)))) {
        ## second solution, starting from right
        c(cb, list(gs(.len, i, .len, right=TRUE)))
      } else {
        cb
      }
    }) |&gt; 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)
}
</code></pre>
<h2>Usage</h2>
<p><em>Two splits:</em></p>
<pre><code>mccg(x='GRGKA')
# <span class="math-container">$`0`
# [1] "GR" "GK" "A" 
# 
# $</span>`1`
# <span class="math-container">$`1`[[1]]
# [1] "GRGK" "A"   
# 
# $</span>`1`[[2]]
# [1] &quot;GR&quot;  &quot;GKA&quot;
# 
# 
# $`2`
# [1] &quot;GRGKA&quot;
</code></pre>
<p><em>Three splits:</em></p>
<pre><code>mccg(x='FLGERTCTRNERQP')
# <span class="math-container">$`0`
# [1] "FLGER" "TCTR"  "NER"   "QP"   
# 
# $</span>`1`
# [1] &quot;FLGERTCTR&quot; &quot;NERQP&quot;    
# 
# <span class="math-container">$`2`
# $</span>`2`[[1]]
# [1] &quot;FLGERTCTRNER&quot; &quot;QP&quot;          
# 
# <span class="math-container">$`2`[[2]]
# [1] "FLGER"     "TCTRNERQP"
# 
# 
# $</span>`3`
# [1] &quot;FLGERTCTRNERQP&quot;
</code></pre>
<p><em>More splits:</em></p>
<pre><code>mccg(x='FLGERTCTRNERQPFLGERTCTRNERQP')
# <span class="math-container">$`0`
# [1] "FLGER"   "TCTR"    "NER"     "QPFLGER" "TCTR"    "NER"     "QP"     
# 
# $</span>`1`
# <span class="math-container">$`1`[[1]]
# [1] "FLGERTCTR"  "NERQPFLGER" "TCTRNER"    "QP"        
# 
# $</span>`1`[[2]]
# [1] &quot;FLGER&quot;       &quot;TCTRNER&quot;     &quot;QPFLGERTCTR&quot; &quot;NERQP&quot;      
# 
# 
# <span class="math-container">$`2`
# $</span>`2`[[1]]
# [1] &quot;FLGERTCTRNER&quot;   &quot;QPFLGERTCTRNER&quot; &quot;QP&quot;            
# 
# <span class="math-container">$`2`[[2]]
# [1] "FLGER"          "TCTRNERQPFLGER" "TCTRNERQP"     
# 
# 
# $</span>`3`
# <span class="math-container">$`3`[[1]]
# [1] "FLGERTCTRNERQPFLGER" "TCTRNERQP"          
# 
# $</span>`3`[[2]]
# [1] &quot;FLGERTCTRNER&quot;     &quot;QPFLGERTCTRNERQP&quot;
# 
# 
# <span class="math-container">$`4`
# $</span>`4`[[1]]
# [1] &quot;FLGERTCTRNERQPFLGERTCTR&quot; &quot;NERQP&quot;                  
# 
# <span class="math-container">$`4`[[2]]
# [1] "FLGERTCTR"           "NERQPFLGERTCTRNERQP"
# 
# 
# $</span>`5`
# <span class="math-container">$`5`[[1]]
# [1] "FLGERTCTRNERQPFLGERTCTRNER" "QP"                        
# 
# $</span>`5`[[2]]
# [1] &quot;FLGER&quot;                   &quot;TCTRNERQPFLGERTCTRNERQP&quot;
# 
# 
# $`6`
# [1] &quot;FLGERTCTRNERQPFLGERTCTRNERQP&quot;
</code></pre>
<h2>Simulation</h2>
<p>Here are some toy data,</p>
<pre><code>set.seed(42)
n &lt;- 1e3
min_chr &lt;- 5; max_chr &lt;- 50
strings &lt;- replicate(n, {
  repeat {
    smp &lt;- sample(setdiff(LETTERS, c(&quot;B&quot;, &quot;Z&quot;, &quot;J&quot;, &quot;O&quot;, &quot;U&quot;, &quot;X&quot;)),
                  sample(min_chr:max_chr, 1L), replace=TRUE)
    if (any(grepl('K|R', smp))) break
  }
  paste(smp, collapse='')
}) |&gt; setNames(nm=_)
</code></pre>
<p>over which the execution of the function can be tested in a loop.</p>
<pre><code>r &lt;- lapply(strings, \(x) tryCatch(mccg(x), error=\(e) NA))
</code></pre>

CC BY-SA 4.0
2023-08-04T11:09:14.447Z
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
Klemens Fröhlich
CC BY-SA 4.0
2023-07-27T05:57:46.090Z
2023-07-27T18:48:34.297Z
2023-07-27T18:48:34.297Z
jay.sf
<p>I am looking for an R based solution to the following problem:</p>
<p>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 &quot;missed cleavages&quot;</p>
<p>In informatics terms I would describe this as skipping a split during a <code>strsplit()</code> function</p>
<p>Example to clarify:</p>
<p>We start out with the protein sequence: <code>GRGKA</code> which is usually split by <code>strsplit()</code> (digested by trypsin) after every K and R</p>
<pre class="lang-r prettyprint-override"><code>&quot;GRGKA&quot; %&gt;% strsplit(., &quot;(?&lt;=K|R)&quot;, perl = TRUE)

[[1]]
[1] &quot;GR&quot; &quot;GK&quot; &quot;A&quot; 
</code></pre>
<p>now imagine that trypsin did not cleave perfectly, for example missing always one K or R</p>
<pre class="lang-r prettyprint-override"><code>&quot;GRGK&quot; &quot;GKA&quot;
</code></pre>
<p>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:</p>
<pre><code>&quot;GRGKA&quot;
</code></pre>
<p>I would love to have a generalizable solution that allows to generate all combinations of sequences that arise with n number of missed cleavages.</p>
<p>If you need any more information please feel free to ask, sorry if I missed something!</p>
<p>jay.sf's requested example:</p>
<pre><code>&quot;FLGERTCTRNERQP&quot;  
</code></pre>
<p>perfect digest:</p>
<pre><code>&quot;FLGER&quot; &quot;TCTR&quot; &quot;NER&quot; &quot;QP&quot;   
</code></pre>
<p>one missed cleavage all combinations:</p>
<pre><code>&quot;FLGERTCTR&quot; &quot;TCTRNER&quot; &quot;NERQP&quot;
</code></pre>
<p>two missed cleavages all combinations:</p>
<pre><code>&quot;FLGERTCTRNER&quot; &quot;QP&quot; &quot;FLGER&quot; &quot;TCTRNERQP&quot; 
</code></pre>
<p>three missed cleavages (there are a total of three K and Rs so trypsin would miss all)</p>
<pre><code>&quot;FLGERTCTRNERQP&quot;
</code></pre>
<p>For longer sequences I think this becomes a frameshift like problem: If we count the K and Rs in the sequence:</p>
<p>FLGE<strong>R</strong>TCT<strong>R</strong>NE<strong>R</strong>QPFLGE<strong>R</strong>TCT<strong>R</strong>NE<strong>R</strong>QP<br />
FLGE<strong>1</strong>TCT<strong>2</strong>NE<strong>3</strong>QPFLGE<strong>4</strong>TCT<strong>5</strong>NE<strong>6</strong>QP<br />
FLGE<strong>R</strong>_TCT<strong>R</strong>_NE<strong>R</strong>_QPFLGE<strong>R</strong>_TCT<strong>R</strong>_NE<strong>R</strong>_QP (perfect digest)<br />
FLGE<strong>R</strong>TCT<strong>R</strong>_NE<strong>R</strong>QPFLGE<strong>R</strong>_TCT<strong>R</strong>NE<strong>R</strong>_QP (one missed cleavage and always the uneven R and Ks)<br />
FLGE<strong>R</strong>_TCT<strong>R</strong>NE<strong>R</strong>_QPFLGE<strong>R</strong>TCT<strong>R</strong>_NE<strong>R</strong>QP  (one missed cleavage and always the even R and Ks)</p>
<p>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....</p>

CC BY-SA 4.0
2023-07-25T16:33:16.813Z
Say `"FLGERTCTRNERQP"` would be a candidate, it has three possible splits. Can you show all possible results for 1, 2, and 3 skipped splits?
jay.sf
CC BY-SA 4.0
2023-07-26T08:47:29.533Z
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
Klemens Fröhlich
CC BY-SA 4.0
2023-07-27T12:06:05.203Z
@KlemensFröhlich the answer below looks cool, you've enough rep to vote - definitely recommend it.
M__
<r><proteins><sequence-analysis>
missed cleavages combinations generation
# missed cleavages combinations generation

(Asked by: username_0 on 2023-07-25)

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

**Tags:** r, proteins, sequence-analysis
> **Comments:**
> **username_1**: Say `"FLGERTCTRNERQP"` would be a candidate, it has three possible splits. Can you show all possible results for 1, 2, and 3 skipped splits?  
> **username_0**: 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  
> **username_2**: @KlemensFröhlich the answer below looks cool, you've enough rep to vote - definitely recommend it.  

---

## 1 answers

### Answer by username_1 on 2023-07-27. Score: 2

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

> **Comments:**
> **username_0**: 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  

---

