Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-08-04T11:09:14.447Z
Answers.Text: 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
Answers.UserDisplayName: Klemens Fröhlich
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-07-27T05:57:46.090Z
Answers.LastActivityDate: 2023-07-27T18:48:34.297Z
Answers.LastEditDate: 2023-07-27T18:48:34.297Z
Answers.UserDisplayName: jay.sf
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-07-25T16:33:16.813Z
Comments.Text: Say `"FLGERTCTRNERQP"` would be a candidate, it has three possible splits. Can you show all possible results for 1, 2, and 3 skipped splits?
Comments.UserDisplayName: jay.sf
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-07-26T08:47:29.533Z
Comments.Text: 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
Comments.UserDisplayName: Klemens Fröhlich
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-07-27T12:06:05.203Z
Comments.Text: @KlemensFröhlich the answer below looks cool, you've enough rep to vote - definitely recommend it.
Comments.UserDisplayName: M__
Tags: <r><proteins><sequence-analysis>
Title: missed cleavages combinations generation
ThreadText: # 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
  
  ```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
  
  ```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!
  
  username_1'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:
  
  FLGE**R**TCT**R**NE**R**QPFLGE**R**TCT**R**NE**R**QP  
  FLGE**1**TCT**2**NE**3**QPFLGE**4**TCT**5**NE**6**QP  
  FLGE**R**\_TCT**R**\_NE**R**\_QPFLGE**R**\_TCT**R**\_NE**R**\_QP (perfect digest)  
  FLGE**R**TCT**R**\_NE**R**QPFLGE**R**\_TCT**R**NE**R**\_QP (one missed cleavage and always the uneven R and Ks)  
  FLGE**R**\_TCT**R**NE**R**\_QPFLGE**R**TCT**R**\_NE**R**QP (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**: @username_0 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 `var`iance 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 username_0  
  
  ---
  
  