Beast系统发育关系和分化时间的估算(上)
以下内容参考了同研究组师弟整理的笔记、和网络资源：Beast v1.8简易手册-修订版 · 语雀 (yuque.com)，根据自己的实际操作整理。
Beast软件是采用贝叶斯演化分析的原理，用来估算系统发育关系和分化时间的软件，使用前提是需要有化石证据或其他的先验信息进行时间校准。
整体流程用到的软件安装都很简单，就是一般的windows系统下安装软件的正常操作，不再赘述。
Beats官网 https://www.beast2.org/
需要的软件包括：
Beast：该软件包包括Beats程序、BEAUTi、TreeAnnotator和其他的应用程序
Tracer：该程序是用来探究Beast结果的，对连续型参数的分布进行量化和可视化，提供可以诊断的信息。
Figtree：用来呈现系统发育树，尤其是通过Beast获取的树。
下面用我自己的一个实例展示一下整个流程。
数据的介绍
我用到的数据是线粒体基因组中的13个蛋白质编码基因和2个rRNA基因（12S和16S）。
因为是多基因，又有蛋白编码基因（编码核苷酸的密码子对应的三个碱基的演化速率存在差异），所以我选择对数据进行分区估算最优核苷酸替换模型。在进行多序列比对时，对两个rRNA基因分别采用mafft alignment，而对13个蛋白质编码基因分别translation align，保证不破坏编码同一核苷酸的三个碱基。之后，将15个比对后的多重序列串联在一起（整个过程在Geneious里操作）。
通过查阅文献，我找到了这里用到的类群的化石证据（包括年代信息），以及前人研究发现的几个类群的分化时间。
除了查阅文献，可以去这两个网站寻找化石证据和分化时间信息：
化石查询网站：https://fossilcalibrations.org/
分化时间查询网站：http://timetree.org/
第零步 估算分区模型
这一步不是必要的一步，所以叫第零步。这一步里我用PartitionFinder这个软件估算分区和每个分区最优的核苷酸替换模型。具体操作可以参考我的这篇文章：
https://www.jianshu.com/p/1a0214ef79a7
第一步 用BEAUTi建立BEAST XML 输入文件
1 载入比对序列的.nex文件（load NEXUS format alignment）
这一步的目的是将我们的序列信息和模型参数信息整合在一起，建立一个输入BEAST进行运算的文件。
首先要把序列文件导出为.nex文件（Geneious里操作）。
打开.nex文件，将分区的信息加入进去，这些信息在PartitionFinder输出的结果文件analysis的best_schem.txt里：
best_scheme.txt文件中展示分区和每个分区的最有替换模型的信息
best_scheme.txt文件中可以直接加入进.nex文件的分区信息
我们打开导出的序列数据的.nex文件，将上述分区信息加入进去，整体的效果如下：
打开BEAUTi，上方工具栏：File-Import Data，导入.nex文件：
可以看到，是按照给出的分区信息被拆分成12个部分。
2 定义校准节点（defining the calibration node）
主要参数的设置选项
我们可以看到主窗口上有很多个选键，我们从左向右依次设置。
点击Taxa，这里用来设置分类单元子集，之后你就可以为设置的分类单元的共同最近祖先（MRCA,most recent common ancester）添加分化时间校准信息。
点击左下角的+，添加分类单元，双击taxon set一列的行名，可以修改名称。刚打开窗口，所有分类群都处于左边的excluded taxa一栏，选择属于该分类单元的分类群，点击右箭头，把它们导入右栏included taxa。
Taxa我们看到在taxon set右边有一个Mono？复选框，如果你确定你设置的分类单元是一个单系类群，可以在复选框里打勾。在之后及进行的MCMC分析中，这一类群就会保持为单系类群。
Mono？
3 设置核苷酸替换模型（setting the substitution model）
因为我这里选择了对数据进行分区，因此需要对之前partitionfinder运算得到的12个subset分别进行模型的设置。点击Partitions，全选后，点击窗口上方Unlink subst. Models，
点击Sites，左面substitution model下面是所有的subset，一个个分别选中后，根据partitionfinder的结果文件内容，在右面设置模型参数：
这里要修改的参数是Substitution Model和Site Heterogeneity Model。很简单，如果最优模型是GTR+I+G+X，Substitution Model选择GTR，Site Heterogeneity Model选择Gamma+Invariant Sites。因为我之前已经对编码氨基酸的三个碱基进行了分区，这里就不用考虑partition into codon positions了，选择off。
4 设置分子钟模型(setting the clock model)
点击Clocks，clock type下拉菜单中选择uncorrelated relaxed clock （注：这里的选择根据个人数据情况）。
5 设置树先验(tree prior)
树先验用于设置种群或物种规模随时间变化的模型。在tree prior的下拉菜单中选择Yule Process。这是一个简单的物种形成模型，更适用于分析来自不同物种的序列。
6 设置先验(Priors)
Priors用于设置模型参数的先验。我们需要基于已知的化石证据或参考文献信息，为一些分化时间设置一个先验分布，点击prior一栏的Using Tree Prior，在弹出的对话框中，设置Prior Distribution为Normal。在mean中填55，stdev中填0.5，意思是最近共同祖先距今55百万年，标准差为0.5个百万年。这里需要根据自己的先验信息设置。
对于一些我们设置的分类群子集，我们可能没有先验信息，就不用设置，可以根据已有的信息进行估算。
7 设置MCMC选项(setting the MCMC options)
Length of chain: this is the number of steps the MCMC will make in the chain before finishing. How long this should be depends on the size of the data set, the complexity of the model and the quality of answer required.
length of chain 指的是MCMC分析运行的总代数，可以通过增加代数提高ESS（effective sample size）值，使参数收敛。默认是10,000,000。
一个粗略的计算方法是总代数=3000*序列条数的平方。这里有20条序列，就是1,200,000代
Log parameters every：for the log file, the value should be set relative to the total length of the chain. Sampling too often will result in very large files with little extra benefit in terms of the precision of the analysis. Sample too infrequently and the log file will not contain much information about the distribution of the parameters.
log parameters every 指的是采样的频率，导致太多无用或意义不大的信息产生，采样不够频繁又不能包含足够多的参数分布信息。一般要根据length of chain的数值来设置，要保证样本数量（样本数量=length of chain/log parameter every）至少有10,000个。
8 生成BEAST XML文件(generating the BEAST XML file)
点击右下角的Generate BEAST File，生成一个.XML文件，然后通过BEAST运算该文件。
打开beast，点击BEAST XML File选项的Choose File，选择生成的XML文件，点击Run。
