Answers.Body: <p>There are different ways to approach the GC bias, the easiest is probably to decide beforehand how many &quot;A/T&quot; and how many &quot;G/C&quot; we need in each sequence, then choose them randomly.</p>
  <pre class="lang-r prettyprint-override"><code># set global parameters
  M &lt;- 7
  I &lt;- 15
  
  # for each sequence, choose a proportion of GC
  prop_gc &lt;- runif(M,
                   min = .7,
                   max = 1)
  
  # randomly choose bases, ensure we have a number that gives us at least `prop_gc`
  bases &lt;- lapply(prop_gc,
                  \(.prop_gc){
                    c(
                      sample(c(&quot;A&quot;,&quot;T&quot;),
                             size = floor(I*(1 - .prop_gc)), replace = TRUE),
                      sample(c(&quot;G&quot;,&quot;C&quot;),
                             size = ceiling(I*.prop_gc), replace = TRUE)
                    )
                  })
  
  # now we can simply take our bases and combine them (in a random order)
  sequences &lt;- sapply(bases,
                      \(.nt) sample(.nt) |&gt; paste(collapse = &quot;&quot;))
  
  sequences
  #&gt; [1] &quot;CGGCTTCCGCGCGGGT&quot; &quot;CTACCCACGGCGCCGG&quot; &quot;GGGCGGTGTATGTGCC&quot; &quot;GGCCCGGGTGCCCAGC&quot;
  #&gt; [5] &quot;GCGGGCGCGGGCGGGT&quot; &quot;CACCCCCTCCCGGGAC&quot; &quot;GTCGAAGCGCGCCAGC&quot;
  </code></pre>
  <p><sup>Created on 2023-10-12 with <a href="https://reprex.tidyverse.org" rel="nofollow noreferrer">reprex v2.0.2</a></sup></p>
  <p>Or same thing with a <code>for</code> loop (since you gave it as a tag):</p>
  <pre class="lang-r prettyprint-override"><code>M &lt;- 7
  I &lt;- 15
  
  sequences &lt;- character(length = M)
  for(i in seq_len(M)){
    prop_gc &lt;- runif(1,
                     min = .7,
                     max = 1)
    
    bases &lt;- c(
      sample(c(&quot;A&quot;,&quot;T&quot;),
             size = floor(I * (1 - prop_gc)), replace = TRUE),
      sample(c(&quot;G&quot;,&quot;C&quot;),
             size = ceiling(I * prop_gc), replace = TRUE)
    )
    
    sequences[[i]] &lt;- sample(bases) |&gt; paste(collapse = &quot;&quot;)
    
  }
  sequences
  #&gt; [1] &quot;CCCCGGTTGCGTTTCC&quot; &quot;CCTCGCTGCGCGTGGC&quot; &quot;CAGTCTCCCGCGGGGC&quot; &quot;GGGCCCGCCGAAAGGC&quot;
  #&gt; [5] &quot;GCGCGGGGCCGGTCAG&quot; &quot;CGGGCGCCCTGGACAC&quot; &quot;ACTCTGCGGGGGAGTC&quot;
  </code></pre>
  <p><sup>Created on 2023-10-12 with <a href="https://reprex.tidyverse.org" rel="nofollow noreferrer">reprex v2.0.2</a></sup></p>
  <p>And we can doublecheck the results:</p>
  <pre class="lang-r prettyprint-override"><code>M &lt;- 7000
  I &lt;- 500
  
  [...]
  
  # check
  gc_prop &lt;- function(x){
    toupper(x) |&gt;
      strsplit(&quot;&quot;) |&gt;
      sapply(\(.char) mean(.char == &quot;C&quot; | .char == &quot;G&quot;))
  }
  
  sapply(sequences, nchar) |&gt; table()
  #&gt; 
  #&gt;  500 
  #&gt; 7000
  
  min(gc_prop(sequences))
  #&gt; [1] 0.702
  </code></pre>
  <p><sup>Created on 2023-10-12 with <a href="https://reprex.tidyverse.org" rel="nofollow noreferrer">reprex v2.0.2</a></sup></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-10-12T18:52:56.560Z
Answers.LastActivityDate: 2023-10-12T19:07:20.660Z
Answers.LastEditDate: 2023-10-12T19:07:20.660Z
Answers.UserDisplayName: Alexlok
Body: <p>I want a script in R, which creates <code>M</code> fragments of DNA with a length of <code>I</code> but with <strong>GC content</strong> of &gt;70%</p>
  <p>I know this gives something but not exactly what I am looking for</p>
  <pre class="lang-r prettyprint-override"><code>require(Biostrings)
  
  set = NULL
  for (i in 1:4)
    set = c(set, paste(sample(DNA_ALPHABET[1:4], 3, replace=T), collapse=&quot;&quot;))
  rep(set,2)
  </code></pre>
  <p>Any idea please?</p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-12T18:56:56.207Z
Comments.Text: Why R though? If the need is to simulate sequences, why do you need to do it R - why not look for the best way to do it using any language/tool?
Comments.UserDisplayName: Ram RS
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-10-12T18:59:01.810Z
Comments.Text: Also by this point you should know better than to use `set` as a variable name.
Comments.UserDisplayName: Ram RS
Tags: <r><for-loop>
Title: A function to produce M fragments of DNA with length of I
ThreadText: # A function to produce M fragments of DNA with length of I
  
  (Asked by: username_0 on 2023-10-12)
  
  I want a script in R, which creates `M` fragments of DNA with a length of `I` but with **GC content** of >70%
  
  I know this gives something but not exactly what I am looking for
  
  ```r
  require(Biostrings)
  
  set = NULL
  for (i in 1:4)
    set = c(set, paste(sample(DNA_ALPHABET[1:4], 3, replace=T), collapse=""))
  rep(set,2)
  
  ```
  
  Any idea please?
  
  **Tags:** r, for-loop
  
  > **Comments:**
  > **username_1**: Why R though? If the need is to simulate sequences, why do you need to do it R - why not look for the best way to do it using any language/tool?  
  > **username_1**: Also by this point you should know better than to use `set` as a variable name.  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2023-10-12. Score: 3 | Accepted answer
  
  There are different ways to approach the GC bias, the easiest is probably to decide beforehand how many "A/T" and how many "G/C" we need in each sequence, then choose them randomly.
  
  ```r
  # set global parameters
  M <- 7
  I <- 15
  
  # for each sequence, choose a proportion of GC
  prop_gc <- runif(M,
                   min = .7,
                   max = 1)
  
  # randomly choose bases, ensure we have a number that gives us at least `prop_gc`
  bases <- lapply(prop_gc,
                  \(.prop_gc){
                    c(
                      sample(c("A","T"),
                             size = floor(I*(1 - .prop_gc)), replace = TRUE),
                      sample(c("G","C"),
                             size = ceiling(I*.prop_gc), replace = TRUE)
                    )
                  })
  
  # now we can simply take our bases and combine them (in a random order)
  sequences <- sapply(bases,
                      \(.nt) sample(.nt) |> paste(collapse = ""))
  
  sequences
  #> [1] "CGGCTTCCGCGCGGGT" "CTACCCACGGCGCCGG" "GGGCGGTGTATGTGCC" "GGCCCGGGTGCCCAGC"
  #> [5] "GCGGGCGCGGGCGGGT" "CACCCCCTCCCGGGAC" "GTCGAAGCGCGCCAGC"
  
  ```
  
  Created on 2023-10-12 with [reprex v2.0.2](https://reprex.tidyverse.org)
  
  Or same thing with a `for` loop (since you gave it as a tag):
  
  ```r
  M <- 7
  I <- 15
  
  sequences <- character(length = M)
  for(i in seq_len(M)){
    prop_gc <- runif(1,
                     min = .7,
                     max = 1)
    
    bases <- c(
      sample(c("A","T"),
             size = floor(I * (1 - prop_gc)), replace = TRUE),
      sample(c("G","C"),
             size = ceiling(I * prop_gc), replace = TRUE)
    )
    
    sequences[[i]] <- sample(bases) |> paste(collapse = "")
    
  }
  sequences
  #> [1] "CCCCGGTTGCGTTTCC" "CCTCGCTGCGCGTGGC" "CAGTCTCCCGCGGGGC" "GGGCCCGCCGAAAGGC"
  #> [5] "GCGCGGGGCCGGTCAG" "CGGGCGCCCTGGACAC" "ACTCTGCGGGGGAGTC"
  
  ```
  
  Created on 2023-10-12 with [reprex v2.0.2](https://reprex.tidyverse.org)
  
  And we can doublecheck the results:
  
  ```r
  M <- 7000
  I <- 500
  
  [...]
  
  # check
  gc_prop <- function(x){
    toupper(x) |>
      strsplit("") |>
      sapply(\(.char) mean(.char == "C" | .char == "G"))
  }
  
  sapply(sequences, nchar) |> table()
  #> 
  #>  500 
  #> 7000
  
  min(gc_prop(sequences))
  #> [1] 0.702
  
  ```
  
  Created on 2023-10-12 with [reprex v2.0.2](https://reprex.tidyverse.org)
  
  
  ---
  
  