# Simulating DNA sequence evolution in R

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?

• 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 Oct 4 '18 at 7:24
• 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 Oct 4 '18 at 14:47
• @KonradRudolph Nice catch! I have fixed the line spacing. – compbiostats Oct 4 '18 at 14:53
• @Llopis The number of substitutions per sequence will not be known in advance. This is determined with rbinom(...). – compbiostats Oct 4 '18 at 19:38

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.

• 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 Oct 3 '18 at 22:33
• 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 Oct 4 '18 at 2:18
• 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 Oct 4 '18 at 2:31
• 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 Oct 4 '18 at 2:36