<p>Mutation is not a random process because there are many mutational biases, e.g. the <a href="https://en.wikipedia.org/wiki/Transition_(genetics)" rel="nofollow noreferrer">transition vs. transversion ratio</a>. Ideally use a proper evolutionary model. Maybe try a tool that incorporates substitution and indel models. I'm not up-to-date on the best packages, but e.g.</p>
<ul>
<li>INDELible <a href="http://abacus.gene.ucl.ac.uk/software/indelible/" rel="nofollow noreferrer">http://abacus.gene.ucl.ac.uk/software/indelible/</a></li>
<li>Evolver <a href="https://www.drive5.com/evolver/" rel="nofollow noreferrer">https://www.drive5.com/evolver/</a> (not to be confused with the PAML package of the same name)</li>
</ul>

CC BY-SA 4.0
2020-08-10T10:58:04.283Z
2020-08-10T10:58:04.283Z
Chris_Rands
<p>I decided that this was too involved for a comment, so I'm making it an answer.</p>
<p>It is still quite difficult to answer the question given the absence of methods explanation; but I think that there's some misunderstanding of phylogenetic modeling going on here, and I can at least address that. I'm going to try to break it down into its components, starting with the basic statistical idea of simulation and then going back to phylogenies.</p>
<p><strong>Statistical reasoning</strong></p>
<p>Let's say that I flip 10 coins (2 outcomes: head (H) and tails (T)), and I observe the sequence <code>HTHTTTHHTT</code>. I then fit a maximum likelihood model, and find that the best model for the data has</p>
<p><span class="math-container">$Pr(H) = 0.4$</span>,</p>
<p><span class="math-container">$Pr(T) = 0.6$</span>,</p>
<p>and that these two outcomes are mutually exhaustive:</p>
<p><span class="math-container">$Pr(H) + Pr(T) = 1.0$</span></p>
<p>I can then use a binomial model to re-simulate data based on these parameters. I will observe various sequences, for example:</p>
<p><code>HTTTTTTHHT</code>
<code>TTTTHHTHHH</code>
<code>HTHHTHHTHH</code>...</p>
<p>However, based purely on the parameters of the model, I don't expect every new simulation to recapitulate the exact same <code>HTHTTTHHTT</code> sequence, because I have reduced its information to some rate parameters for the simulation.</p>
<p>Each of those new sequences will yield model parameters (<span class="math-container">$Pr(H), Pr(T)$</span>) that are similar to the training data, but they will not necessarily be particularly similar to the training data in sequence.</p>
<p><strong>Phylogenetic modeling</strong></p>
<p>Phylogenetic modeling works in the same way, such that you reduce the observations at the tips of the tree into some model with some parameters. Most of the time you are estimating e.g. a rate matrix of parameters for transition between states for a set of residues (in simplest case for DNA 4 nucleotide states and for proteins, 20 amino acid states usually). If you are working with gene sequences you are probably using the 4 nucleotide states or some derivation thereof.</p>
<p>Simulations usually work by starting with some fixed ancestral state at the root (which can take various states), and then evolving forward from there along each branch of the tree using the rate matrix in a <a href="https://en.wikipedia.org/wiki/Continuous-time_Markov_chain" rel="nofollow noreferrer">continuous-time Markov chain</a> to &quot;<a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3509945/" rel="nofollow noreferrer">replay the tape of life</a>&quot;, by assuming that the model you've fit somewhat reflects truth.</p>
<p>While opinions differ somewhat on what to expect from &quot;replaying the tape of life&quot; (see that linked review for more), basically no one expects to get the same results as the input data, or even 90% similar to the input data. That would indicate that evolution is highly deterministic, which just isn't reliably true. Of course there are situations where such determinism happens to be the case, but they are relatively rare, and I don't think that sequence evolution is considered to be very deterministic at all.</p>
<p>It's a little different in that the tree defines a correlation structure of the tips, which constrains the simulation results to be somewhat more similar to the training data than we expect from coin-flipping. But there's still no reason to expect particularly similar results on the tips.</p>
<p>But you can also model discrete characters on a tree, which can be even simpler. For example, I used the <a href="https://rdrr.io/cran/ape/" rel="nofollow noreferrer">R APE package</a> to simulate the following character states (A or B for each tip) on a simulated tree of 10 tips with the following commands:</p>
<pre><code>tr = rtree(10)
c_tr = rTraitDisc(tr)
</code></pre>
<p><a href="https://i.sstatic.net/sfQVk.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/sfQVk.png" alt="enter image description here" /></a></p>
<p>I can then build a model for this character and this tree using the following:</p>
<pre><code>rates = ace(c_tr, tr, type=&quot;discrete&quot;)
</code></pre>
<p>and use it to simulate new sets of states:</p>
<pre><code># A (index 1) is estimated to be ancestral state with 99% confidence
new_sim = rTraitDisc(tr, rate=rates$rates, states=c(&quot;A&quot;, &quot;B&quot;), root.value=1)
plot(tr)
tiplabels(new_sim, adj = c(1,0))
</code></pre>
<p><a href="https://i.sstatic.net/ei9Rk.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ei9Rk.png" alt="enter image description here" /></a></p>
<p>Look! Not a single observation of state B! Let's try it again:</p>
<pre><code>new_sim2 = rTraitDisc(tr, rate=rates$rates, states=c(&quot;A&quot;, &quot;B&quot;), root.value=1)
plot(tr)
tiplabels(new_sim2, adj = c(1,0))
</code></pre>
<p><a href="https://i.sstatic.net/CYM3y.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/CYM3y.png" alt="enter image description here" /></a></p>
<p>Ok, this is more similar, but there's still another state B on the tips.</p>
<p><strong>Conclusion</strong></p>
<p>Evolution is not very deterministic. The correlation structure of the tree makes the outcomes less random than independent trials as in coin-flipping, but it is not at all expected that simulation results from phylogenetic models will match the training observations well. The best that can be said is that each simulated dataset can hopefully be <em>reduced to the same model parameters</em> as the input data, but that's a very different statement from yielding the exact same outcomes.</p>
<p>If what you are looking for from simulation is to yield the exact same sequence states at the tips, then you are unlikely to get that from phylogenetic simulation. You might occasionally get pretty close, but that is not at all the expectation.</p>

CC BY-SA 4.0
2020-08-16T21:34:26.843Z
2020-08-16T21:34:26.843Z
Maximilian Press
<p>I'm working on understanding the evolution of a gene in a phylogenetic tree. I know the rates of evolution of each organism (from the tree).</p>
<p>I am taking a random DNA sequence with my gene, evolving it along the tree with the said rates and comparing this to the actual gene in these various organisms.
I am 'evolving' the gene if randomly at a probability determined by the rate at each step in the tree.</p>
<p>The gene from my simulation is different from what is actually in these organisms. How do I modify my model, so that my simulation gives the same results?</p>

CC BY-SA 4.0
2020-08-09T21:18:02.753Z
If mutation is a random process, why would you expect to reproduce the same result at the end of a stochastic simulation? Most simulation methods will not give you the exact same results as the training data; the "truth" is only one instantiation of many possibilities. Can you say anything about the methods you are using?
Maximilian Press
CC BY-SA 4.0
2020-08-10T08:33:43.173Z
@MaximilianPress I agree, there will not be a 100% fit. Even 90% similarity is good enough.
Using a phylogenetic tree, every time there is a branch, I copy the genome and evolve it independently based on branch length. i.e. if a random number is less than branch length, the nucleotide changes (for each position)
Anthony Guterres
CC BY-SA 4.0
2020-08-11T10:31:48.620Z
What is the purpose of the sequence simulation?
M__
CC BY-SA 4.0
2020-08-17T17:55:30.337Z
@AnthonyGuterres see my new answer, it was going to be too hard to fit into a comment. Again, hard to know how to modify your model without knowing anything about your model.
Maximilian Press
<phylogenetics><simulated-data>
How to modify DNA evolution model to fit actual data?
# How to modify DNA evolution model to fit actual data?

(Asked by: username_0 on 2020-08-09)

I'm working on understanding the evolution of a gene in a phylogenetic tree. I know the rates of evolution of each organism (from the tree).

I am taking a random DNA sequence with my gene, evolving it along the tree with the said rates and comparing this to the actual gene in these various organisms.
I am 'evolving' the gene if randomly at a probability determined by the rate at each step in the tree.

The gene from my simulation is different from what is actually in these organisms. How do I modify my model, so that my simulation gives the same results?

**Tags:** phylogenetics, simulated-data
> **Comments:**
> **username_1**: If mutation is a random process, why would you expect to reproduce the same result at the end of a stochastic simulation? Most simulation methods will not give you the exact same results as the training data; the "truth" is only one instantiation of many possibilities. Can you say anything about the methods you are using?  
> **username_0**: @MaximilianPress I agree, there will not be a 100% fit. Even 90% similarity is good enough. Using a phylogenetic tree, every time there is a branch, I copy the genome and evolve it independently based on branch length. i.e. if a random number is less than branch length, the nucleotide changes (for each position)  
> **username_2**: What is the purpose of the sequence simulation?  
> **username_1**: @AnthonyGuterres see my new answer, it was going to be too hard to fit into a comment. Again, hard to know how to modify your model without knowing anything about your model.  

---

## 2 answers

### Answer by username_3 on 2020-08-10. Score: 1

Mutation is not a random process because there are many mutational biases, e.g. the [transition vs. transversion ratio](https://en.wikipedia.org/wiki/Transition_(genetics)). Ideally use a proper evolutionary model. Maybe try a tool that incorporates substitution and indel models. I'm not up-to-date on the best packages, but e.g.

* INDELible <http://abacus.gene.ucl.ac.uk/software/indelible/>
* Evolver <https://www.drive5.com/evolver/> (not to be confused with the PAML package of the same name)


---

### Answer by username_1 on 2020-08-16. Score: 1

I decided that this was too involved for a comment, so I'm making it an answer.
It is still quite difficult to answer the question given the absence of methods explanation; but I think that there's some misunderstanding of phylogenetic modeling going on here, and I can at least address that. I'm going to try to break it down into its components, starting with the basic statistical idea of simulation and then going back to phylogenies.
Statistical reasoning
Let's say that I flip 10 coins (2 outcomes: head (H) and tails (T)), and I observe the sequence HTHTTTHHTT. I then fit a maximum likelihood model, and find that the best model for the data has
$Pr(H) = 0.4$,
$Pr(T) = 0.6$,
and that these two outcomes are mutually exhaustive:
$Pr(H) + Pr(T) = 1.0$
I can then use a binomial model to re-simulate data based on these parameters. I will observe various sequences, for example:
HTTTTTTHHT
TTTTHHTHHH
HTHHTHHTHH...
However, based purely on the parameters of the model, I don't expect every new simulation to recapitulate the exact same HTHTTTHHTT sequence, because I have reduced its information to some rate parameters for the simulation.
Each of those new sequences will yield model parameters ($Pr(H), Pr(T)$) that are similar to the training data, but they will not necessarily be particularly similar to the training data in sequence.
Phylogenetic modeling
Phylogenetic modeling works in the same way, such that you reduce the observations at the tips of the tree into some model with some parameters. Most of the time you are estimating e.g. a rate matrix of parameters for transition between states for a set of residues (in simplest case for DNA 4 nucleotide states and for proteins, 20 amino acid states usually). If you are working with gene sequences you are probably using the 4 nucleotide states or some derivation thereof.
Simulations usually work by starting with some fixed ancestral state at the root (which can take various states), and then evolving forward from there along each branch of the tree using the rate matrix in a continuous-time Markov chain to "replay the tape of life", by assuming that the model you've fit somewhat reflects truth.
While opinions differ somewhat on what to expect from "replaying the tape of life" (see that linked review for more), basically no one expects to get the same results as the input data, or even 90% similar to the input data. That would indicate that evolution is highly deterministic, which just isn't reliably true. Of course there are situations where such determinism happens to be the case, but they are relatively rare, and I don't think that sequence evolution is considered to be very deterministic at all.
It's a little different in that the tree defines a correlation structure of the tips, which constrains the simulation results to be somewhat more similar to the training data than we expect from coin-flipping. But there's still no reason to expect particularly similar results on the tips.
But you can also model discrete characters on a tree, which can be even simpler. For example, I used the R APE package to simulate the following character states (A or B for each tip) on a simulated tree of 10 tips with the following commands:
tr = rtree(10)
c_tr = rTraitDisc(tr)


I can then build a model for this character and this tree using the following:
rates = ace(c_tr, tr, type="discrete")

and use it to simulate new sets of states:
# A (index 1) is estimated to be ancestral state with 99% confidence
new_sim = rTraitDisc(tr, rate=rates$rates, states=c("A", "B"), root.value=1)
plot(tr)
tiplabels(new_sim, adj = c(1,0))


Look! Not a single observation of state B! Let's try it again:
new_sim2 = rTraitDisc(tr, rate=rates$rates, states=c("A", "B"), root.value=1)
plot(tr)
tiplabels(new_sim2, adj = c(1,0))


Ok, this is more similar, but there's still another state B on the tips.
Conclusion
Evolution is not very deterministic. The correlation structure of the tree makes the outcomes less random than independent trials as in coin-flipping, but it is not at all expected that simulation results from phylogenetic models will match the training observations well. The best that can be said is that each simulated dataset can hopefully be reduced to the same model parameters as the input data, but that's a very different statement from yielding the exact same outcomes.
If what you are looking for from simulation is to yield the exact same sequence states at the tips, then you are unlikely to get that from phylogenetic simulation. You might occasionally get pretty close, but that is not at all the expectation.


---

