<p>What are you expecting these lines to do?</p>

<pre><code>num.transi &lt;- rbinom(n = 1, size = length.seqs, prob = transi.rate)
...
num.transv &lt;- rbinom(n = 1, size = length.seqs, prob = transv.rate)
</code></pre>

<p>Using the variables that are specified in your code, they are extremely zero-skewed. This is not too surprising given the mutation rate and sequence length.</p>

<p>I expect that increasing the sequence length, or increasing the transition rate, should result in more haplotypes being generated.</p>

<p>After running this through myself, I see that increasing the mutation rate does increase the number of haplotypes, but the skew is retained. This is because each mutation creates a new haplotype, following something like an infinite allele model.</p>

<p>I guess this is not particularly surprising. It's unlikely that two random mutations from a particular haplotype will result in the same haplotype, and therefore any mutations (from <em>any</em> existing haplotype) will generally create something new.</p>

CC BY-SA 4.0
2018-10-03T22:33:59.497Z
These lines compute the number of transitions and transversions across the sequence alignment of a given length according to a binomial distribution. You're correct that more haplotypes will likely be generated when length.seqs and/or transi.rate are increased. However, I'm looking to ensure that there are more sequences per haplotype (so that the distribution is not so skewed).
compbiostats
CC BY-SA 4.0
2018-10-04T02:18:57.143Z
You're complaining about a skew towards the primary haplotype, i.e. that there are not enough mutations. My recommended fix is to either increase the sequence length or increase the transition rate.
gringer
CC BY-SA 4.0
2018-10-04T02:31:32.343Z
Been playing with the code... still get the same effect when mutation rate and/or sequence length is increased. There should be a better way - I'll ponder it and experiment some more.
compbiostats
CC BY-SA 4.0
2018-10-04T02:36:03.557Z
My specific application is in DNA barcoding, which itself relies on short (~500 - ~650 bp) DNA fragments. Increasing length.seqs is not really an option - unless, should it work, I then trim the alignment evenly on both 5' and 3' ends to the required length. I'll ponder this.
compbiostats
CC BY-SA 4.0
2018-10-03T20:54:34.027Z
2018-10-05T05:10:26.430Z
2018-10-05T05:10:26.430Z
gringer
<p>I hope someone can lend their thoughts on the below code to generate DNA sequences under the Kimura-2-Parameter model of DNA substitution.</p>

<p>The issue is that each time the code is run and the haplotype distribution is examined, there always is a very skewed distribution. </p>

<p>What I would like is for a different distribution to be generated each time the code is run. that is, there should be multiple specimens sharing the same haplotype (instead of most sharing the most common haplotype). </p>

<p>Here is my code:</p>

<pre><code>library(pegas)
set.seed(17)


sim.seqs &lt;- TRUE
length.seqs &lt;- 500
num.seqs &lt;- 100 # number of DNA sequences
subst.model &lt;- "K80" # nucleotide substitution model
transi.rate &lt;- 1e-4 # transition rate
transv.rate &lt;- transi.rate / 2 # transversion rate


if (sim.seqs == TRUE) {

  nucl &lt;- as.DNAbin(c('a','c','g','t'))
  res &lt;- sample(nucl, size = length.seqs, replace = TRUE, prob = rep(0.25, 4))

 if (subst.model == "K80") {
    transi.set &lt;- list('a' = as.DNAbin('g'),
                       'c' = as.DNAbin('t'),
                       'g' = as.DNAbin('a'),
                       't' = as.DNAbin('c'))

    transv.set &lt;- list('a' = as.DNAbin(c('c', 't')),
                       'c' = as.DNAbin(c('a', 'g')),
                       'g' = as.DNAbin(c('c', 't')),
                       't' = as.DNAbin(c('a', 'g')))

    transi &lt;- function(res) {
      unlist(transi.set[as.character(res)])
    }

    transv &lt;- function(res) {
      sapply(transv.set[as.character(res)], sample, 1)
    }

    duplicate.seq &lt;- function(res) {
      num.transi &lt;- rbinom(n = 1, size = length.seqs, prob = transi.rate) # total number of transitions

      if (num.transi &gt; 0) {
        idx &lt;- sample(length.seqs, size = num.transi, replace = FALSE)
        res[idx] &lt;- transi(res[idx])
      }

      num.transv &lt;- rbinom(n = 1, size = length.seqs, prob = transv.rate) # total number of transversions

      if (num.transv &gt; 0) {
        idx &lt;- sample(length.seqs, size = num.transv, replace = FALSE)
        res[idx] &lt;- transv(res[idx])
      }

      res
      }

    }

  res &lt;- matrix(replicate(num.seqs, duplicate.seq(res)), byrow = TRUE, nrow = num.seqs)

  class(res) &lt;- "DNAbin"

  # write.dna(res, file = "seqs.fas", format = "fasta")

  h &lt;- sort(haplotype(res), decreasing = TRUE, what = "frequencies")
  rownames(h) &lt;- 1:nrow(h)

}
</code></pre>

<h2>Output</h2>

<pre><code>h

Haplotypes extracted from: res

Number of haplotypes: 5
     Sequence length: 500

Haplotype labels and frequencies:

 1  2  3  4  5
96  1  1  1  1
</code></pre>

<p>If the code is run multiple times (without the seed), the same pattern emerges: h is always skewed toward the most dominant haplotype.</p>

<p>For example, I want to be able to run the code (without seed) and get something like</p>

<pre><code>h

Haplotypes extracted from: res

Number of haplotypes: 5
     Sequence length: 500

Haplotype labels and frequencies:

  1    2    3    4   5
  35  25   20   15   5
</code></pre>

<p>i think I would have to specify a distribution for the haplotypes, but the number of haplotypes generated by the mutation process is not known beforehand.</p>

<p>Any thoughts?</p>

CC BY-SA 4.0
2018-10-04T07:24:16.123Z
Similar to what gringer said: Perhaps I understood wrong your model, but why do you only expect one transition and transversion? Shouldn't that be `rbinom( n = 500, ...`? So, in a sequence of 500 you expect to know how many substitutions are in all the sequence, not if there is a mutation on the sequence.
llrs
CC BY-SA 4.0
2018-10-04T14:47:33.010Z
Hi Jarrett, the empty double line spacing unfortunately makes your code harder to read than necessary. It would help if you used empty lines only to logically separate blocks, not between all lines.
Konrad Rudolph
CC BY-SA 4.0
2018-10-04T14:53:21.083Z
@KonradRudolph Nice catch! I have fixed the line spacing.
compbiostats
CC BY-SA 4.0
2018-10-04T19:38:18.753Z
@Llopis The number of substitutions per sequence will not be known in advance. This is determined with rbinom(...).
compbiostats
<r><phylogenetics>
Simulating DNA sequence evolution in R
# Simulating DNA sequence evolution in R

(Asked by: username_0 on 2018-10-03)

I hope someone can lend their thoughts on the below code to generate DNA sequences under the Kimura-2-Parameter model of DNA substitution.
The issue is that each time the code is run and the haplotype distribution is examined, there always is a very skewed distribution. 
What I would like is for a different distribution to be generated each time the code is run. that is, there should be multiple specimens sharing the same haplotype (instead of most sharing the most common haplotype). 
Here is my code:
library(pegas)
set.seed(17)


sim.seqs <- TRUE
length.seqs <- 500
num.seqs <- 100 # number of DNA sequences
subst.model <- "K80" # nucleotide substitution model
transi.rate <- 1e-4 # transition rate
transv.rate <- transi.rate / 2 # transversion rate


if (sim.seqs == TRUE) {

  nucl <- as.DNAbin(c('a','c','g','t'))
  res <- sample(nucl, size = length.seqs, replace = TRUE, prob = rep(0.25, 4))

 if (subst.model == "K80") {
    transi.set <- list('a' = as.DNAbin('g'),
                       'c' = as.DNAbin('t'),
                       'g' = as.DNAbin('a'),
                       't' = as.DNAbin('c'))

    transv.set <- list('a' = as.DNAbin(c('c', 't')),
                       'c' = as.DNAbin(c('a', 'g')),
                       'g' = as.DNAbin(c('c', 't')),
                       't' = as.DNAbin(c('a', 'g')))

    transi <- function(res) {
      unlist(transi.set[as.character(res)])
    }

    transv <- function(res) {
      sapply(transv.set[as.character(res)], sample, 1)
    }

    duplicate.seq <- function(res) {
      num.transi <- rbinom(n = 1, size = length.seqs, prob = transi.rate) # total number of transitions

      if (num.transi > 0) {
        idx <- sample(length.seqs, size = num.transi, replace = FALSE)
        res[idx] <- transi(res[idx])
      }

      num.transv <- rbinom(n = 1, size = length.seqs, prob = transv.rate) # total number of transversions

      if (num.transv > 0) {
        idx <- sample(length.seqs, size = num.transv, replace = FALSE)
        res[idx] <- transv(res[idx])
      }

      res
      }

    }

  res <- matrix(replicate(num.seqs, duplicate.seq(res)), byrow = TRUE, nrow = num.seqs)

  class(res) <- "DNAbin"

  # write.dna(res, file = "seqs.fas", format = "fasta")

  h <- sort(haplotype(res), decreasing = TRUE, what = "frequencies")
  rownames(h) <- 1:nrow(h)

}

Output
h

Haplotypes extracted from: res

Number of haplotypes: 5
     Sequence length: 500

Haplotype labels and frequencies:

 1  2  3  4  5
96  1  1  1  1

If the code is run multiple times (without the seed), the same pattern emerges: h is always skewed toward the most dominant haplotype.
For example, I want to be able to run the code (without seed) and get something like
h

Haplotypes extracted from: res

Number of haplotypes: 5
     Sequence length: 500

Haplotype labels and frequencies:

  1    2    3    4   5
  35  25   20   15   5

i think I would have to specify a distribution for the haplotypes, but the number of haplotypes generated by the mutation process is not known beforehand.
Any thoughts?

**Tags:** r, phylogenetics
> **Comments:**
> **username_1**: Similar to what gringer said: Perhaps I understood wrong your model, but why do you only expect one transition and transversion? Shouldn't that be `rbinom( n = 500, ...`? So, in a sequence of 500 you expect to know how many substitutions are in all the sequence, not if there is a mutation on the sequence.  
> **username_2**: Hi Jarrett, the empty double line spacing unfortunately makes your code harder to read than necessary. It would help if you used empty lines only to logically separate blocks, not between all lines.  
> **username_0**: @KonradRudolph Nice catch! I have fixed the line spacing.  
> **username_0**: @Llopis The number of substitutions per sequence will not be known in advance. This is determined with rbinom(...).  

---

## 1 answers

### Answer by username_3 on 2018-10-03. Score: 1 | Accepted answer

What are you expecting these lines to do?
num.transi <- rbinom(n = 1, size = length.seqs, prob = transi.rate)
...
num.transv <- rbinom(n = 1, size = length.seqs, prob = transv.rate)

Using the variables that are specified in your code, they are extremely zero-skewed. This is not too surprising given the mutation rate and sequence length.
I expect that increasing the sequence length, or increasing the transition rate, should result in more haplotypes being generated.
After running this through myself, I see that increasing the mutation rate does increase the number of haplotypes, but the skew is retained. This is because each mutation creates a new haplotype, following something like an infinite allele model.
I guess this is not particularly surprising. It's unlikely that two random mutations from a particular haplotype will result in the same haplotype, and therefore any mutations (from any existing haplotype) will generally create something new.

> **Comments:**
> **username_0**: These lines compute the number of transitions and transversions across the sequence alignment of a given length according to a binomial distribution. You're correct that more haplotypes will likely be generated when length.seqs and/or transi.rate are increased. However, I'm looking to ensure that there are more sequences per haplotype (so that the distribution is not so skewed).  
> **username_3**: You're complaining about a skew towards the primary haplotype, i.e. that there are not enough mutations. My recommended fix is to either increase the sequence length or increase the transition rate.  
> **username_0**: Been playing with the code... still get the same effect when mutation rate and/or sequence length is increased. There should be a better way - I'll ponder it and experiment some more.  
> **username_0**: My specific application is in DNA barcoding, which itself relies on short (~500 - ~650 bp) DNA fragments. Increasing length.seqs is not really an option - unless, should it work, I then trim the alignment evenly on both 5' and 3' ends to the required length. I'll ponder this.  

---

