Simulating DNA sequences in R with a given value of $\theta = 4N_{e}\mu$

Question

This might not be the most appropriate site to be asking such a question, but perhaps someone has a solution.

My question is: is there an R package or function for simulation of DNA sequences of a given basepair length generated according to a prespecified value of the population mutation rate $\theta = 4N_{e}\mu$, where $N_{e}$ is the effective population size and $\mu$ is the per-generation mutation rate

I have real sequence data from GenBank, but I would also like to simulate some random DNA sequences according to a coalescent process.

Essentially the function would output the generated aligned DNA sequences as a FASTA file.

I have been rather unsuccessful in tracking down a suitable package, as most only output generated phylogenies.

Answer 1

There are a number of tools but I don't know any that comes with an R library. They are usually all called from the command line. I don't either of any that produce FASTA file. SimBit (my own software) can produce vcf files which can easily be converted into FASTA with PGDspider.

However, if your simulations are so easy, you might not need these individual-based softwares. Just write the code yourself. It might only be take a few tens of lines. Here is a simple example

freq = 0.1 # initial frequency
N = 1000   # Constant population size
mu = 1e-7  # mutation rate

for (generation in 1:nbGenerations)
{
   # Drift
   freq = rbinom(1,N,freq) / N

   # Mutations
   oneWayMutations   = rbinom(1,N*freq,mu) / N
   otherWayMutations = rbinom(1,N*(1-freq),mu) / N
   freq = freq - oneWayMutations + otherWayMutations
}

If you need something more complex (such as several loci), then you should probably use one of the existing individual-based software. See the post Sequence evolution simulation tool