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

CC BY-SA 4.0
2023-10-12T18:52:56.560Z
2023-10-12T19:07:20.660Z
2023-10-12T19:07:20.660Z
Alexlok
<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>

CC BY-SA 4.0
2023-10-12T18:56:56.207Z
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?
Ram RS
CC BY-SA 4.0
2023-10-12T18:59:01.810Z
Also by this point you should know better than to use `set` as a variable name.
Ram RS
<r><for-loop>
A function to produce M fragments of DNA with length of I
# 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
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.
# 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
Or same thing with a for loop (since you gave it as a tag):
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
And we can doublecheck the results:
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


---

