I am simulating some DNA sequences (containing characters A, C, G, and T) in R through specifying

i) the number of sequences to generate (num.seqs)

ii) the length of the generated sequences (length.seqs)

ii) the DNA mutation rate (mu.rate)

The mutation rate is distributed as a binomial (can be approximated by a Poisson distribution) in the following way

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

If num.muts is nonzero, I then do the following

idx <- sample(length.seqs, size = num.muts, replace = FALSE) # randomly place mutations along the DNA sequences

The end result is a distribution:

Number of unique DNA sequences (h) | 1 2 3 4 ....
Number of non-unique sequences     | a b c d ... (a + b + c + d + ... = num.seqs; a, b, c, d... > 0)

To illustrate, suppose

num.seqs <- 100
length.seqs <- 500
mu.rate <- 1e-4

which may give the distribution

 1 2 3 4
97 1 1 1

What I'm really after is a more 'realistic' distribution:

 1  2  3 4
53 26 13 8

which is clearly a multinomial.

Caveat: the total number of unique DNA sequences (h) is not known a priori. This is clearly a problem for the multinomial, as the prob argument to that function is a vector of length h.

rmultinom(n, size, prob)

My question is: With the code I currently have, how can I obtain a distribution akin to the second one shown above using the rmultinom() function?

Note: Increasing either num.seqs, length.seqs and/or mu.rate always results in the first distribution shown above (or a similarly positive skewed one).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.