Bayesian Tip Dating using MrBayes

研究背景

推断生物类群的系统发生关系和分异时间是分支系统学分析的基础。如何合理利用化石的形态和地质年代数据来完成此类推断一直是一个棘手的问题。近些年开发的贝叶斯支端定年法 (Bayesian tip dating) 为该问题提供了一个可行的解决方案。支端定年法把化石形态和年代数据整合在一次完整的计算过程中，能够尽可能充分地利用数据信息，同时考虑了树的拓扑结构、分异时间、演化速率以及化石年代的不确定性。该方法通过统计模型来描述特征演化、类群生灭、化石采样等过程，并借助相对成熟的贝叶斯统计框架和计算方法来进行参数估计和模型选择。
       本文着重分析离散的形态特征数据，这类数据是支端定年法的数据核心。除了形态特征之外，支端定年法法还可以结合现生类群的分子序列数据。这类整合的数据分析也被称为全证据定年 (total-evidence tip dating) (Ronquist et al., 2012; Gavryushkina et al., 2014; Zhang et al., 2016)。使用MrBayes进行全证据定年分析的流程可以参考另一篇教程 (Zhang, 2019)。
       下面按照MrBayes命令的顺序，首先设置描述特征状态演替的Mk模型 (Lewis, 2001)，然后设置描述时间树的石化生灭过程 (fossilized birth-death process) 模型 (Stadler, 2010)，接着设置描述特征演化速率的宽松形态钟 (relaxed clock) 模型，最后设置估计参数后验分布的马氏链蒙特卡罗 (Markov chain Monte Carlo，MCMC) 算法。除了下文提到的命令外，MrBayes还有强大的帮助系统便于查询相关命令的解释，即在MrBayes里使用help命令。本文所用数据为中生代鸟类的形态特征矩阵 (Zhang and Wang, 2019)。

仪器设备

1. 个人电脑
2. 操作系统：Windows、macOS或Linux。
3. MrBayes采用命令行执行，没有图形界面。
4. Windows中使用命令提示符 (command prompt)，macOS或Linux使用终端 (terminal) (图1)。
   

软件信息

MrBayes目前最新的版本为3.2.7，可以从如下GitHub链接下载：https://github.com/NBISweden/MrBayes/releases。下载后找到可执行文件，放在任意所需的目录下即可。此外，软件也可以从源代码自行编译，INSTALL文件中提供了具体方法，不在此处赘述。


图1. Windows命令提示符和macOS终端的界面

操作步骤

1. 准备数据和MrBayes命令文件
   数据和命令采用NEXUS格式 (Maddison et al., 1997)。这是一种纯文本格式，可以使用文本编辑器打开。本文所用的中生代鸟类的形态矩阵完整数据见附件 (Zhang and Wang, 2019)。该矩阵包含68个物种，280个特征 (图2)。
   
   
   图2. NEXUS数据部分截图
   
   MrBayes命令可以写在数据文件的结尾，也可以写在单独的文件中。本文假设命令写在文件run.nex里，独立于数据文件birds.nex。这样的好处是可以设置多个命令文件对同一数据完成不同的分析。下面逐行解释MrBayes命令文件。
   
2. MrBayes命令区块以Begin MrBayes; 开始，以End; 结尾。每条命令以分号结束。为了便于区分，本文中命令使用等款字体。
   2.1

   首先读取数据文件 (中括号内容为注释，软件不执行)。这里假设数据和命令 文件在同一文件夹下。
   execute birds.nex; [read in data matrix]
   
   2.2

   接着设置特征演替模型为Mkv (Lewis, 2001)，并且特征之间速率异质性服从伽马分布 (Yang, 1994)。Mkv模型是离散形态特征目前唯一可以使用的模型。特征默认都是无序的 (unordered，即从当前状态可直接变为任意其它状态)， 对于有序的 (ordered) 特征 (即从当前状态只能直接变为相邻的状态)，需要 根据其在矩阵的列号单独设置如下：
   [substitution model]
       ctype ordered: 1 3 8 28 31 43 51 56 64 67 69 70 72 74 92 107 117 159 168 176 183 193 205 213 214 216 219 222 229 233 234 249 261 265 268 270;
       lset coding = variable rates = gamma; [Mkv+Gamma]
   
   2.3

   数据中的地质年代以百万年前 (Ma) 为单位。根据化石的地质年代的上下界， 设定为均匀分布 (uniform)。如果化石年代是确定的，则设为固定的时间 (fixed)。
   [tip dates]
        calibrate
        Dromaeosauridae = uniform(66.0, 167.7)
        Archaeopteryx = uniform(145.0, 152.1)
        Jeholornis = uniform(110.6, 125.0)
        Chongmingia = uniform(110.6, 120.3)
        Sapeornis = uniform(110.6, 125.0)
        Confuciusornis_sanctus = uniform(110.6, 125.0)
        Changchengornis = uniform(120.3, 125.0)
        Eoconfuciusornis = uniform(129.0, 130.7)
        Confuciusornis_dui = uniform(120.3, 125.0)
        Yangavis = uniform(120.3, 125.0)
        Boluochia = uniform(110.6, 120.0)
        Concornis = uniform(120.3, 125.0)
        Elsornis = uniform(70.6, 84.9)
        Eoalulavis = uniform(120.3, 125.0)
        Cathayornis = uniform(110.6, 120.3)
        Eocathayornis = uniform(110.6, 120.3)
        Eoenantiornis = uniform(120.3, 125.0)
        Gobipteryx = uniform(72.1, 83.6)
        Longipteryx = uniform(110.6, 120.3)
        Longirostravis = uniform(120.3, 125.0)
        Neuquenornis = uniform(79.5, 83.5)
        Pengornis = uniform(110.6, 120.3)
        Eopengornis = uniform(129.0, 130.7)
        Protopteryx = uniform(129.0, 130.7)
        Rapaxavis = uniform(110.6, 120.3)
        Shanweiniao = uniform(120.3, 125.0)
        Vescornis = uniform(121.6, 122.5)
        Qiliania = uniform(112.6, 125.4)
        Dunhuangia = uniform(112.6, 125.4)
        Piscivorenantiornis = uniform(110.6, 120.3)
        Linyiornis = uniform(110.6, 120.3)
        Sulcavis = uniform(110.6, 120.3)
        Bohaiornis = uniform(110.6, 120.3)
        Longusunguis = uniform(110.6, 120.3)
        Shenqiornis = uniform(121.6, 122.5)
        Zhouornis = uniform(110.6, 120.3)
        Parabohaiornis = uniform(110.6, 120.3)
        Fortunguavis = uniform(110.6, 120.3)
        Pterygornis = uniform(110.6, 120.3)
        Cruralispennia = uniform(129.0, 130.7)
        Monoenantiornis = uniform(120.3, 125.0)
        Archaeorhynchus = uniform(110.6, 125.0)
        Schizooura = uniform(110.6, 120.3)
        Bellulornis = uniform(110.6, 120.3)
        Vorona = uniform(66.0, 72.1)
        Jianchangornis = uniform(110.6, 120.3)
        Songlingornis = uniform(110.6, 120.3)
        Longicrusavis = uniform(120.3, 125.0)
        Apsaravis = uniform(72.1, 83.6)
        Hongshanornis = uniform(120.3, 125.0)
        Archaeornithura = uniform(129.0, 130.7)
        Parahongshanornis = uniform(110.6, 120.3)
        Tianyuornis = uniform(120.3, 125.0)
        Yanornis = uniform(110.6, 120.3)
        Patagopteryx = uniform(79.5, 83.5)
        Yixianornis = uniform(120.3, 125.0)
        Piscivoravis = uniform(110.6, 120.3)
        Iteravis = uniform(110.6, 120.3)
        Gansus = uniform(113.0, 120.3)
        Ichthyornis = uniform(70.6, 94.3)
        Hesperornis = uniform(66.0, 84.9)
        Parahesperornis = uniform(70.6, 84.9)
        Enaliornis = uniform(100.5, 113.0)
        Baptornis_advenus = uniform(85.8, 89.3)
        Baptornis_varneri = uniform(70.6, 84.9)
        Vegavis = uniform(66.0, 68.0)
        ;
        prset nodeagepr = calibrated;
   
   2.4

   石化生灭过程 (fossilized birth-death process) 作为时间树 (timetree) 的先验，是用于描述类群分异 (speciation or birth)、灭绝 (extinction or death)、 化石采样 (fossil sampling) 和现生类群采样 (extant sampling) 的过程。对于现生类群 (鸡和鸭)，这里使用了多样化采样 (diversified sampling) (Zhang et al., 2016)，它比随机采样 (random sampling) 更符合这个数据。对于化石采样速率和生灭速率，可以假设是恒定的 (=diversity后面直接分号结束)，也可以是随时间变化的 (=diversity后面按化石采样、分异和灭绝为顺序分别设置变化时间点个数和对应的时间)。例如，这里设置了化石采样速率按时间分为四段，分异速率分为两段，灭绝速率恒定。现生类群的采样比例固定为0.0002，即样本中鸡和鸭两种占所有现生鸟物种的比例。
          MrBayes对生灭速率和化石采样速率重新参数化为：a净成种速率 (分异速率减去灭绝速率)，b相对灭绝速率 (灭绝速率除以分异速率)，c相对化石采样速率 (化石采样速率除以灭绝速率和化石采样速率之和)。净成种速率取值从0到无穷，后两者取值都是从0到1，所以对净成种速率使用了均值为0.01的指数分布先验 (范围从0到无穷)，对相对灭绝速率和相对化石采样速率使用了贝塔分布先验 (参数都为1；此时该分布为从0到1的均匀分布)。
           最后，石化生命过程以根部 (最近共同祖先) 时间为起始条件，也需设置先验分布。这里参考数据中最古老的化石 (驰龙)，设置为最小值153百万年前，均值169百万年前的偏移指数分布 (范围从153到无穷，且随时间增加概率密度降低)。除了指数分布外，MrBayes还有均匀分布，伽马分布，正态分布，对数正态分布可根据实际情况来选择。
   [fossilized birth-death model]
   prset brlenspr = clock:fossilization;
   [diversified extant sampling, followed by 3 fossil sampling rate shifts and 1 speciation rate shift at 66 Ma, extinction rate is assumed constant]
   prset samplestrat = diversity 3: 145 100 66, 1: 66;
   prset sampleprob = 0.0002; [two out of about 10000 extant species]
   prset speciationpr = exp(100); [net diversification rate, d]
   prset extinctionpr = beta(1, 1); [turnover rate, v]
   prset fossilizationpr = beta(1, 1); [relative fossil-sampling rate, s]
   prset treeagepr = offsetexp(153, 169); [root age]
   
   2.5

   宽松钟模型描述了形态特征在树枝上演化速率的变化模式。这里采用了独立伽马分布 (即白噪音white noise) 模型 (Lepage et al., 2007)，它比自相关的对数正态分布模型更适合这个数据 (Zhang and Wang, 2019)。基准速率 (base clock rate) 设置了一个弱信息先验 (伽马分布)，均值为0.02标准差为 0.014， 即每百万年平均每50个特征中有一个发生变化但方差也较大。除了伽马分布外，基准速率还有指数分布、正态分布和对数正态分布备选。不同的分布在概率密度函数的形状上有所区别，但此处作为基准速率的先验 (只要是弱信息的先验) 来说一般影响不大，因为化石的形态和年代数据提供了较多的信息。更多细节可参考 (Zhang and Wang, 2019)。
   [relaxed clock model]
        prset clockratepr = gamma(2, 100); [base clock rate]
        prset clockvarpr = igr; [independent gamma rates]
   
   2.6

   在定年分析中，构建的树为有根树，树根的位置默认由时钟模型确定。但是由于演化速率在树上可能变化很剧烈从而导致单单依靠宽松钟无法准确定树根。这时，需要使用附加的拓扑结构限制。下例通过限制鸟类为内群，把第一个类群 (驰龙) 排除为外群。接着还定义了几个大的分类单元。最后使用了其中四个限制来约束树的拓扑结构，以增加树的解析程度。
   [topology constraints]
        constraint Aves = 2-.; [ingroup]
        constraint Pygostylia = 4-.;
        constraint Ornithothoraces = 11-68;
        constraint Enantiornithines = 11-41;
        constraint Ornithuromorpha = 42-68;
        prset topologypr = constraint(Aves, Pygostylia, Enantiornithines, Ornithuromorpha);
   
   2.7

   最后来设置MCMC运行的参数。MrBayes默认同时运行两个独立运算，每个运算使用四条MCMC链 (Metropolis-coupled MCMC)，其中一条为冷链，用于收集参数样本，其余为热链，用于帮助MCMC跨越后验分布的多峰。这里设置了MCMC链长为5千万代，每2千代取样一次，每5万代往屏幕输出一条信息，每50万代计算一次两个运算中支系频率的平均标准差 (average standard deviation of split frequencies，ASDSF)。ASDSF体现了两次运算中树的拓扑结构的相似程度，如果两次运算都收敛，理论上ASDSF应该接近于0。一般把ASDSF小于0.01作为MCMC运算收敛的指标之一。另外一个重要指标是参数的有效样本大小 (ESS)，一般至少需要ESS都大于100。ESS 连同参数的后验分布估计 (包括均值，中位数和95%可信区间) 由sump得到，而sumt用于总结50%多数合意树 (majority-rule consensus tree)。它们默认都舍去MCMC样本的前25%作为burn-in。
   [MCMC settings]
      mcmcp nrun = 2 nchain = 4 ngen = 50000000 samplefr = 2000 printfr = 50000 diagnfr = 500000;
      mcmc filename = birds;
   
      [summarize tree and parameters, default burnin fraction 0.25]
      sumt;
      sump;
   
3. 准备好数据和命令文件后，就可以开始计算了。这里假设Windows中MrBayes执行文件的名称为mb.exe，macOS或Linux中的名称为mb，数据文件名称为birds.nex，命令文件名称为run.nex，它们都在同一目录下。
          在命令行中，先用cd命令进入该文件夹，然后运行MrBayes软件开始运算 (图1) 。正常情况下，应该即可看到MrBayes的输出 (图3)。等待计算结束得到结果。 程序在运行过程中会输出剩余时间，这个数据使用个人电脑大概需要一到两天的时间。使用MPI并行计算可以在几小时完成 (并行版本的MrBayes需要使用 macOS或Linux系统从源代码编译)。
          计算结束后，会生成一个以 .tre结尾的文件，这就是从后验样本总结的50% 多数合意树。该文件可用FigTree打开浏览 (http://tree.bio.ed.ac.uk/software/figtr ee/)。结果可以参考 Zhang and Wang (2019)。
   
   
   图3. MrBayes起始运算的屏幕输出
   
   

小结与建议

本文演示了如何利用MrBayes进行贝叶斯支端定年分析。大致流程主要包括下载软件，准备数据，准备命令文件、运行软件和查看结果。
        由于MCMC算法是随机的计算机方法，因此每次计算的结果多少会有不同。理想情况下，不管哪次运算都是对真实后验分布的估计，因此，不同独立运算的结果应该相差不大。这可以通过比较两次运算的参数样本 (如使用Tracer软件打开 .p结尾的文 件)，同时参考ASDSF和ESS值。
        如果发现两次运算的结果差异很大，很可能MCMC算法还没有收敛 (convergence)，或者不同运算卡在了后验分布的不同峰值附近导致混合 (mixing) 很差。前者一般可以通过增加MCMC链长 (ngen) 来改善。后者则相对困难，有时通过增加热链个数 (nchain) 并调整热链温度 (temperature) 能改善，但大多数情况下需要调整模型和先验的设置，选择更合理的模型往往能改善后验分布的形状，从而使MCMC能够更快收敛并很好混合。

致谢

感谢bio-protocol 的编辑工作以及同行专家的宝贵修改意见。

参考文献

1. 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.
2. Lepage, T., Bryant, D., Philippe, H. and Lartillot, N. (2007). A general comparison of relaxed molecular clock models. Mol Biol Evol 24: 2669-2680.
3. Lewis, P. O. (2001). A likelihood approach to estimating phylogeny from discrete morphological character data. Syst Biol 50: 913-925.
4. Maddison, D. R., Swofford, D. L. and Maddison, W. P. (1997). NEXUS: an extensible file format for systematic information. Syst Biol 46: 590-621.
5. Ronquist, F., Klopfstein, S., Vilhelmsen, L., Schulmeister, S., Murray, D. L. and Rasnitsyn, A. P. (2012). A total-evidence approach to dating with fossils, applied to the early radiation of the Hymenoptera. Syst Biol 61: 973-999.
6. Stadler, T. (2010). Sampling-through-time in birth-death trees. J Theor Biol 267: 396-404.
7. Yang, Z. (1994). Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. J Mol Evol 39: 306-314.
8. Zhang, C. (2019). Molecular Clock Dating using MrBayes. Vertebrata PalAsiatica 57: 241-252.
9. Zhang, C. and Wang, M. (2019). Bayesian tip dating reveals heterogeneous morphological clocks in Mesozoic birds. Roy Soc Open Sci 6: 182062.
10. Zhang, C., Stadler, T., Klopfstein, S., Heath, T. and Ronquist, F. (2016). Total-evidence dating under the fossilized birth-death process. Syst Biol65: 228-249.