Simulating DNA sequence evolution in R

Question

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?

Answer 1

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.