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