text: user3450: 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?
  
  user73: 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.
  
original_text: 3450: <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>
   73: <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>
  