BEAST2推断分化时间

分子钟：一个特定的生物大分子（蛋白质或DNA）在所有的演化谱系中具有恒定的演化速率。其中，演化谱系是指不同类群。Beast2软件基于多序列比对后的结果，按照mcmc（马尔可夫蒙特卡洛方法）构建贝叶斯进化树，而ML（最大似然法）和贝叶斯进化树对氨基酸/核苷酸替代模型的选择非常敏感，故在进行进化树或分化时间构建之前，需对氨基酸/核苷酸替代模型进行选择。
分析案例：

该属物种间存在明显的网状进化信号 | Edelman et al., 2019

在序列数据的贝叶斯分析中，先验扮演重要角色。 如果未正确指定先验，则可能会导致运行花费很长时间来收敛，根本无法收敛，或者会导致推断的树和模型参数有偏差。本例中，我们探讨如何选择先验条件，以及如何使用2009年在美国传播的流感病毒中的H3N2甲型流感病毒数据进行分子时钟校准。分子时钟模型旨在估算数据的替代率，重要的是要了解在什么情况下使用哪种模型以及何时进行分子校准。

1. 所需软件

（1）BEAST2 - Bayesian Evolutionary Analysis Sampling Trees: 是一个免费软件包，可使用MCMC进行贝叶斯分子序列的进化分析，使用有根树进行推理。
（2）BEAUti2 - Bayesian Evolutionary Analysis Utility: 用于生成BEAST2 XML配置文件的图形用户界面工具。BEAST2和BEAUti2都是Java程序，这两个程序在Windows和Linux上将具有相同的布局和功能。 BEAUti2作为BEAST2软件包的一部分提供，无需单独安装。
（3）TreeAnnotator: 用于汇总树的后验样本，以生成最大的进化枝可信度树。 它也可以用于总结和可视化其他树参数（例如节点高度）的后验估计。为BEAST2软件包的一部分，无需单独安装。
（4）Tracer
（5）FigTree

2. 使用BEAUti创建分析文件

以Wang et.al,2015中L148_r树的taxa为基础，去除重复的taxa，最终得到122个taxa。DNA序列数据L122.nex, 起始树为ML法推断出的besttree。
打开BEAUti，将DNA序列数据和起始树导入程序，按以下步骤生成跑BEAST所需的xml文件。
（1）定义tMRCA,共4个。

Outgroup（2taxa,2种派模蛛,非强制单系）
Linyphiidae（120taxa，所有传统分类中定义的皿蛛科物种，非强制单系）。
Linyphiidae（117taxa，except stemonyphantes亚科的3个taxa, 强制单系，用于设定化石矫正点）
Orsonwelles（2taxa, include Orsonwelles属两种，强制单系）

（2）核苷酸分子替换模型。Site Model：GTR。事先使用jModelTest等软件对该最适模型及其参数进行判定。此部分可参考 http://blog.sciencenet.cn/blog-460481-791370.html
（3）分子钟模型Clock Models。 将Molecular Clock Model选为“Relaxed Clock:Uncorrelated Lognormal”, estimate不选择。
（4）设置Priors。选择Tree Prior为Speciation:Birth-Death Model。
（5）设置Prior。在tmrca（Linyphiinae）中选择Lognormal，参数为0，2，125（参考Dimitrov，2012）。 在tmrca（Osonwelles）中选择Lognormal，参数为0，0.001，2.7。
（6）设置起始树。在Beast2中，无法从BEAUti的GUI界面设置起始树。唯一的做法是生成XML文件后，手工修改代码，将默认使用随机树做起始树的代码改为如下格式
<init id="StartingTree.t:FPV" initial="@Tree.t:FPV" spec="beast.util.TreeParser" IsLabelledNewick="true">
<input name="newick">yournewick; </input>
</init>
在引入自定义的Starting tree后，tMRCA的设定便有了限制，tMRCA的taxa分组不能和Starting tree拓扑结构冲突，否则在BEAST的时候会报错。
（7）Operators。如果你仅仅希望用Beast2推算分化时间，树拓扑要保持与ML分析的最优树完全一致，在这部分你要将XML文件代码中<OPERATORS>部分Subtree slide，Narrow exchange，Wide exchange， Wilson Balding四行代码全部删掉，以保证在整个mcmc过程中树拓扑保持一致，只是支长有变化。
（8）设置MCMC。链长设为40000000，Echo设为1000，Log_every设为1000（一般来讲，tracer中显示所有参数的ESS值均大于200表示BEAST分析结果可信，需要根据此参数确定最适合链长）。
（9）生成XML文档。

3. 运行BEAST

在BEAST中打开xml文件，运行后生成后缀为trees和log的文件。在/beast/bin目录下找到beast程序，运行BEAST。弹出的对话框打开编辑好的xml文件，就会自动运行beast，结果trees和log文件存储在bin目录下。

4. 运行Tracer分析结果。

用Tracer打开log文件分析，Tracer显示了BEAST后得到的各种参数，包括likelihood，meanRate(平均进化速率)，coefficientOfVariation(每一分支分子进化的变异系数，表示分子变化的速率在不同分支上是否等速)，root和每一个tMRCA分化时间的估计值等等。同时选中root.height和tMRCA，可以通过箱线图和边界密度图显示预设的化石校正点估计值是否一致。

5. 运行TreeAnnotator得到合议树。

用TreeAnnotator打开trees文件，设定参数得到合议树。TreeAnnotator中Output File文件默认是让选择，其实直接取一个输出的文件名点ok即可。需要注意的是，由于jvm有内存限制，对于超过10000棵树需要和议时，不能直接运行treeannotator，而是进入beast的lib目录，运行以下代码启动treeannotator。

java -Xmx2048m -Xms2048m -classpath beast.jar dr.app.tools.TreeAnnotator

合议树第一个参数设为4000（40000000代，每1000代生成1棵树，共40000棵。取10%为4000）.第二个参数改为0

6. 用FIGTREE打开tree文件

为了使结果更为直观，可以在Figtree中显示tMRCA中node的值及95%HPD.
参考：
https://www.jianshu.com/p/ddea8e486121
https://blog.csdn.net/zhanyongjia_cnu/article/details/50618131
https://github.com/Taming-the-BEAST/Prior-selection
https://www.sohu.com/a/209191217_761120