I'm looking for a little bit of guidance.

My question is regarding the simulation of DNA sequences with a fix substitution rate.

The majority of the programs for simulating sequences use Continuous Markov Chain Models, with different instantaneous rate matrices that can be modified according to the free parameters of each model (JC69, HKY, GTR, etc..). For example, in HKY the free parameters are the transition/transversion (kappa) rate and nucleotide frequencies (pi).

In these programs (for example INDELible, Seq-Gen, Pyvolve), the substitution rate is expected to be represented by the branch length of the phylogenetic tree used to create the simulations.

INDELible paper.

"..Rate matrices are rescaled by INDELible such that the branch lengths represent the expected number of substitutions per site (or the average expected number of substitutions per site under a heterogeneous-sites model)." (P2 - Simulation of substitutions)

INDELible: a flexible simulator of biological sequence evolution

"...Each branch length is assumed to denote the mean number of nucleotide substitutions per site that will be simulated along that branch..." (P2 - Algorithm)

Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees

Let's say I want to simulate sequences with a mean substitution rate of 2.5 x 10^-7 per year. How would I prepare a tree that can accurately represent that substitution rate and be used in the simulation? In other words, how can the length be represented as a substitution rate accurately?

Any insight is appreciated!. Thanks in advance!


UPDATE based on comment

Simplest answer: just rescale your substitution rate from years to units of branch length. (Stop thinking in terms of years.)

original answer

The short answer is that you have to calibrate the branch lengths of the tree to yield the substitution rates per year you want. (Leave the topology alone.)

Tree branch lengths are usually agnostic to chronological time (e.g. years). Instead it is very focused on whatever the model used to construct the tree is. This can be simple substitution counts (e.g. Jukes-Cantor), or a sophisticated model that takes into account a lot of information about relative rates of different changes, possibly learned from the dataset itself (e.g. GTR).

You yourself will have to decide what the conversion factor is between your branch lengths (presumably derived from such a substitution model or a distance metric) and units of chronological time. In some trees, it may already be calibrated to units of chronological time. But those will be rare.

So what I'd do is to make sure that the branch lengths of the tree are of a length that will, on average, yield the number of substitutions that you expect based on how long (in years) a given branch on the tree is supposed to be. For example, if you have a tree of primates, then you want root to tip distances to be on the order of tens of millions of years, probably. If you have a tree of all bacteria, you want the root-to-tip distances to be on the order of billions of years.

Note please that this is identical to re-scaling the substitution matrix to substitutions per year. What is crucial is that you know what you mean by "substitutions per year" and how you are defining that in terms of the tree.

  • $\begingroup$ Please correct me if I'm wrong. Considering that branch length is agnostic to time, if i have a tree with certain "value" of length (e.g 0.04),Is the substitution model that defines the "time" itself? (years or millions of years)?. If I only want to simulate DNA sequences in a defined topology, it would be better to obtain the model parameters from a real dataset? (With IQTREE for example). Thanks again for your help $\endgroup$
    – Kako
    Nov 10 '21 at 20:03
  • $\begingroup$ @Kako Note that there is no scientific reason why chronological time has to be represented in a phylogeny or a substitution matrix. the substitution model is also usually thinking only in terms of branch length as the time scale (for phylogenetics at least). Neither branch lengths nor substitution models are usually understood in terms of years. So if you want your simulation to use years, then your value 0.04 should be years to match your substitution model. Or, you could just rescale your substitution rate to be subs/unit branch length (much simpler, what all those programs assume). $\endgroup$ Nov 10 '21 at 20:40

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.