Using BEAST 2 for Bayesian Tip Dating

研究背景

生物类群的历史分异时间一直是生物学研究的基础和热点问题。但由于逝去历史的诸多不确定性，如何准确推断分异时间长期以来也面临各种挑战。伴随分子数据的涌现、化石信息的积累、以及演化模型 (如分子钟) 的发展，贝叶斯定年法近些年得到了广泛关注和应用。该方法将置换模型、演化速率、类群分异过程等纳入一个整体的统计学分析框架中，将化石信息纳入先验模型，从而估算类群的绝对历史分异时间。
        常用的贝叶斯定年方法可以分为两类，一类是节点定年 (node dating)，另一类是支端定年 (tip dating)，其主要区别在于如何利用化石的信息。节点定年法把化石信息转化为概率分布来校准系统发生树上的部分内部节点，继而估计其它内部节点的分异时间 (e.g., Drummond et al., 2006)。这一类方法通常应用于仅包含现生类群的分析。支端定年法则直接利用化石信息，把化石和现生类群同时作为系统发生树的支端进行分析，从而估计树中内部节点的分异时间 (Pyron, 2011; Ronquist et al., 2012)。
        支端定年相比节点定年在理论方面存在诸多优势。比如：首先，支端定年可以充分利用化石和现生类群的形态数据，从而可以开展全证据支端定年 (total-evidence tip dating) 分析。由于一般无法获得化石的分子数据，因此化石的系统发生位置由形态数据决定。其次，支端定年可以潜在利用与所研究类群相关的所有化石记录，无需像节点定年仅挑选与某节点有关的最古老的化石记录。第三，支端定年只需根据化石自身年代信息设置时间校准先验，无需对某节点的历史分异时间特别挑选统计分布以描述其时间先验。特别伴随描述类群分异及化石采样过程的石化生灭过程 (fossilized birth-death process, FBD) 模型的提出 (Stadler, 2010)，支端定年得到了广泛应用。具体分析还可根据实际情况针对各历史时间段设置不同分异和采样速率，以及根据现生类群取样情况设置多样化采样策略等 (Gavryushkina et al., 2014; Zhang et al., 2016)。多项研究提示：贝叶斯支端定年在推断类群历史分异时间方面具有良好的准确性和精确性 (如Gavryushkina et al., 2014; Zhang et al., 2016; Luo et al., 2020)。
        借助BEAST 2 软件包 (Bouckaert et al., 2019) 及Tracer (Rambaut et al., 2018)、FigTree (http://tree.bio.ed.ac.uk/software/figtree/) 等软件，本文将展示贝叶斯支端定年的基本操作步骤和要点信息。

仪器设备

BEAST 2、Tracer、Figtree在Windows、Linux、macOS等主流计算机操作系统下均可运行，界面基本一致。由于篇幅所限，本文仅展示macOS Big Sur (11.0.1) 系统下的操作界面和流程。

软件版本信息及下载地址

BEAST是根据马尔可夫链蒙特卡罗 (MCMC) 算法开展贝叶斯进化分析的开源免费软件。其虽然主要用于推断生物类群的历史分异时间，但也可以用于构建系统发生树、重建类群祖先性状、估计种群大小、实施模型选择等等 (Bouckaert et al., 2019)。和BEAST 1不同，BEAST 2采用崭新的架构编写，突出模块化，可以通过加载模块或软件包从而拥有多种功能。BEAST 2软件包中除主要实施MCMC分析的BEAST程序外，还拥有BEAUti、LogComiber、TreeAnnotator等多个程序可以在前期参数设置、后期数据整理等方面提供帮助。
        BEAST 2软件包可以从官方网站https://www.beast2.org下载最新版本。其安装步骤请参照下载文件夹中的README.txt文件。本文采用撰写时的最新版本BEAST v.2.6.3。BEAST 2的正常运行依赖于Java v8或更高版本。因此，需首先安装Java以保证BEAST 2的正常运行。Tracer和FigTree在本文用于将BEAST分析结果可视化，从而易于对结果进行判定和解析。其分别可以从https://www.beast2.org/tracer-2/和http://tree.bio.ed.ac.uk/software/figtree/下载最新版本。本文采用Tracer v1.7.1和FigTree v1.4.4。

实验步骤

1. 实验数据
   出于演示目的，本文采用计算机模拟的数据，其源于Luo et al. (2020)，包含50个现生类群 (或物种) 和7个化石类群 (或物种)。分子数据 (基因序列) 部分包含5个分区 (partition)，分别为1_1st，1_2nd，1_3rd，1_4th，1_5th，化石物种无分子数据，由"?"替代。形态数据部分由"0"或"1"代表化石和现生类群的离散性状特征 (图1，附件信息https://github.com/ArongLuo/Protocol-)。
   
   
   图1. BEAUti加载数据及初步设置后所示
   
   
2. BEAUti参数设置
   2.1

   在BEAST 2安装目录 (如/Applications/BEAST 2.6.3/) 中双击BEAUti图标，打开BEAUti。
   2.2

   基本的全证据支端定年分析依赖模块SA (sampled ancestor) 和MM (morphological model)，所以需要首先加载安装。在File中单击Manage Packages，在软件包列表中查找SA和MM并进行安装 (图2)。
   2.3

   在BEAUti中通过Import Alignment和Add Morphological Data分别打开分子和形态数据。假设所有的数据源于相同的类群演化过程，因此对6个数据分区 (分子数据5个分区和形态数据) Link Trees，并命名为tree；假设形态和分子数据具不同的进化速率，则只对分子数据的5个分区Link Clock Models，并分别对分子和形态部分命名为mol和1_morpho_v。由于描述各分区的碱基置换模型 (或特征变换模型) 不同或参数相对独立，则对各分区Unlink Site Models (图1)。
   
   
   图2. BEAST 2 软件包列表及所需安装的SA和MM
   
   
   2.4

   点击BEAUti主选项Tip Dates，选中Use tip dates。在准备该数据文件时，已经将序列名称"_"后设置为类群的历史时间，如t9_0代表现生物种t9，et50_84代表存在于84百万年左右的化石物种et50。在时间单位为年的假设下，选择Before the present，并实际以百万年为真实时间单位。根据序列名称特点，在此通过Auto-configure增加各物种的时间信息 (图3)。值得注意的是，类群的时间也可以手动填写或修改，后续对速率先验的设置需与真实的时间单位保持一致。
   
   
   图3. 在BEAUti中设置各物种的历史时间信息
   
   
   2.5

   点击BEAUti主选项Site Model开始位点变换模型设置。在此，分子数据各分区均设为HKY + G，形态数据采用MK模型 (图4)。
   2.6

   点击BEAUti主选项Clock Model开始演化速率设置。对形态和分子数据分别采用宽松钟模型，可具体设为Relaxed Clock Log Normal。该模型假设演化速率服从独立的对数正态分布。
   2.7

   点击BEAUti主选项Priors开始先验设置。多数参数可使用其默认先验，但1)需特别选择FBD模型作为支端定年的树先验，并可通过Condition On Rho Sampling使树先验受限于已知的现生类群取样百分比Rho，本数据分析中Rho为1 (图5)；2) 需考虑FBD模型所涉参数的先验设置，如无明确信息，可对diversificationRate、samplingProportion、turnover分别采用信息较弱的统计分布如exponential (10)、beta (1,1)、beta (1,1)；3) FBD模型一般基于类群根部时间或起源时间，在此对类群起源时间设置先验Uniform (94,300)，其中94是所有7个化石物种历史时间的最大值。4) 在已选宽松分子钟模型基础上设置速率先验，将ucldMean设置为uniform (10e-6, 1)，ucldStdev保留默认。值得注意的是，由于各类群年龄在2.4中是以百万年为真实单位，速率在此以每百万年为单位。
   
   
   图4. 在BEAUti中对形态和分子数据分别设置碱基置换模型HKY + G和特征变换模型Mk
   
   
   2.8

   点击BEAUti主选项MCMC开始MCMC参数设置，可包括代数、记录频率、文件名称等。在此设Chain Length为100百万代，树文件和文本文件的Log Every是5000。上述各项设置完成后，即可保存.xml文件。
   
   
3. BEAST数据运行
   在BEAST 2安装目录中双击BEAST图标，打开BEAST后在弹出页面通过Choose File加载之前保存成功的.xml文件，点击Run。如无意外，BEAST即开始MCMC分析 (图6)。为了保证结果的可靠性，一般需要对同一个数据至少独立运行两次以保证后验概率得到收敛。
   
   
   图5. 在BEAUti中设置FBD模型为树先验
   
   
   图6. BEAST成功运行界面
   
   
4. 结果分析
   BEAST一个独立运行结果一般包含三个文件：.log文件，.trees文件，.state文件。结果分析中一般仅需考虑.log文件和.trees文件。在此，r1和r2两个文件夹 (两次独立运行) 分别有1_beast.log和1_beast.trees。
   4.1

   在Tracer安装目录 (一般为/Applications/) 双击Tracer图标，通过File下拉菜单Import Trace File分别加载r1和r2中的1_beast.log文件。加载后，即可查看posterior、prior等参数的各种统计值如Mean，Median，95% HPD等。特别的，可通过查看参数的ESS数值来获悉各参数在MCMC分析中的有效独立取样大小；理想情况下，ESS值一般大于200。在逐一查看.log文件时，也可以同时选中两个.log文件比较两者的分析结果，也可以点击Combined检查两者的结果是否收敛。通过Tracer分析，如果判定两次独立分析结果比较理想 (图7)，即可以开展后续分析。
   4.2

   在BEAST 2安装目录中双击LogCombiner图标打开LogCombiner。将File type设定为Tree Files，点击"+"加载r1和r2中两个.trees文件并可将Burnin percentage设为25 (与Tracer默认Burnin percentage为10不同)。对Output File设置路径和命名后，点击run即可将两个.trees文件合并为一个文件。
   
   
   图7. Tracer提示BEAST两次独立运行情况良好
   
   
   4.3

   在BEAST 2安装目录中双击TreeAnnotator图标打开TreeAnnotator。如果需要获得Maximum clade credibility tree，一般将Node heights设置为Median heights，在对Input Tree File和Output File分别设置后，点击Run即可开始分析。由于在树文件合并时已经对Burnin有所设置，在此只需保留Burnin percentage为0。
   4.4

   在FigTree安装目录双击FigTree图标，通过File/Open加载4.3产生的Maximum clade credibility tree文件。通过FigTree各选项可将类群的系统发生关系和历史分异时间进行展示。如图8可较清晰呈现估算的类群根部时间95%置信区间、化石类群的系统发生位置等等。该步骤当然也可以用其它软件如DensiTree进行实现。
   
   
   图8. FigTree软件显示类群的系统发生关系和历史分异时间
   
   

小结与建议

本文向读者展示了使用BEAST 2软件包进行贝叶斯全证据支端定年分析的基本流程。值得注意的是，贝叶斯支端定年分析也可在化石和现生类群形态特征数据缺失的情况下实施 (Heath et al., 2014)。该情况下，类似节点定年，一般需要对化石的系统发生位置设置拓扑限制；后验概率则会在兼顾拓扑限制同时，根据FBD先验统计化石系统发生位置的各种可能。具体操作可参考：https://taming-the-beast.org/tutorials/FBD-tutorial/。这种 (非全证据) 的贝叶斯支端定年特别适用于分子数据不断涌现但形态特征数据相对稀缺的系统基因组学时代。
        另需注意的是，虽然支端定年存在诸多理论优势，但由于其涉及形态特征和化石数据，所以在实际应用中也可能面临不少挑战。比如，化石的形态特征数据会存在不完整或片碎化等特点，容易对化石的系统发生位置推断造成偏倚；用于描述形态特征变换的Mk/Mkv模型过于简单，往往不能充分反映形态特征的复杂特点；形态特征是否遵循形态钟仍有待研究，等等。所以，对于贝叶斯支端定年的实际分析结果，特别是依据形态数据推断的化石系统发育位置，仍需结合其它证据 (如古生物学证据) 等综合考虑。

竞争性利益声明

本文作者无利益纷争。

致谢

作者首先感谢bio-protocol为本文提供了发表平台，同时也特别感谢同行专家提出的宝贵修改意见。

参考文献

1. Bouckaert, R., Vaughan, T. G., Barido-Sottani, J., Duchene, S., Fourment, M., Gavryushkina, A., Heled, J., Jones, G., Kuhnert, D., De Maio, N., Matschiner, M., Mendes, F. K., Muller, N. F., Ogilvie, H. A., du Plessis, L., Popinga, A., Rambaut, A., Rasmussen, D., Siveroni, I., Suchard, M. A., Wu, C. H., Xie, D., Zhang, C., Stadler, T. and Drummond, A. J. (2019). BEAST 2.5: an advanced software platform for Bayesian evolutionary analysis. PLOS Comput Biol 15: e1006650.
2. Drummond, A. J., Ho, S. Y. W., Phillips, M. J. and Rambaut, A. (2006). Relaxed phylogenetics and dating with confidence. PLOS Biol 4: e88.
3. Gavryushkina, A., Welch, D., Stadler, T. and Drummond, A. J. (2014). Bayesian inference of sampled ancestor trees for epidemiology and fossil calibration. PLOS Comput Biol 10: e1003919.
   
4. Heath, T. A., Huelsenbeck, J. P. and Stadler, T. (2014). The fossilized birth-death process for coherent calibration of divergence-time estimates. Proc Natl Acad Sci USA 111: E2957-E2966.