Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-03T22:33:59.497Z
Answers.Text: 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).
Answers.UserDisplayName: compbiostats
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-04T02:18:57.143Z
Answers.Text: 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.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-04T02:31:32.343Z
Answers.Text: 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.
Answers.UserDisplayName: compbiostats
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-04T02:36:03.557Z
Answers.Text: 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.
Answers.UserDisplayName: compbiostats
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-10-03T20:54:34.027Z
Answers.LastActivityDate: 2018-10-05T05:10:26.430Z
Answers.LastEditDate: 2018-10-05T05:10:26.430Z
Answers.UserDisplayName: gringer
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-10-04T07:24:16.123Z
Comments.Text: 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.
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-10-04T14:47:33.010Z
Comments.Text: 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.
Comments.UserDisplayName: Konrad Rudolph
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-10-04T14:53:21.083Z
Comments.Text: @KonradRudolph Nice catch! I have fixed the line spacing.
Comments.UserDisplayName: compbiostats
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-10-04T19:38:18.753Z
Comments.Text: @Llopis The number of substitutions per sequence will not be known in advance. This is determined with rbinom(...).
Comments.UserDisplayName: compbiostats
Tags: <r><phylogenetics>
Title: Simulating DNA sequence evolution in R
ThreadText: # 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 username_3 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**: @username_2 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.  
  
  ---
  
  