1
$\begingroup$

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).

$\endgroup$

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.