基因序列中的密码子随着时间的推移以几乎恒定的比例进行替换，具有恒定的演化速率，因此两个物种时间的遗传距离与物种分歧时间成正比。

本文就利用mcmctree估算物种分歧时间做一个介绍。mcmctree 属于PAML软件包，之前我们介绍过PAML中的codeml. paml计算 KaKs值。

关于mcmctree说明文档可查看http://abacus.gene.ucl.ac.uk/software/MCMCtree.Tutorials.pdf

关于PAML软件的安装可查看之前codml的笔记即可。

本文以软件的实例数据为例子进行介绍

1 示例文件练习

示例文件位于paml4.9j/examples/DatingSoftBound 下；

输入文件主要包括3个文件：

• 序列比对文件（nucl或者pep都可以）
• 带有化石校准的系统发育树
• mcmctree.ctl； 控制文件

（1）序列比对文件

mtCDNApri123.txt 内容如下所示，是一个phlip格式

该示例文件分为3个区块，密码子在3个位点的多序列比对结果。

此外，也可以输入密码子或氨基酸序列比对信息，则需要再配置文件中修改相应的参数来表明数据类型 。若使用氨基酸序列来进行分析，由于mcmctree命令不能选择较好的氨基酸替换模型进行分析，需要自己手动运行codeml进行分析后，在生成中间文件用于运行mcmctree，比较麻烦。因此，推荐按上述方法使用密码子三位点的碱基序列作为输入信息进行PAML分析。

(2) 带有校准点的有根树

示例文件如下

 7  1

((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);

第一行，7个物种，1颗树，中间用空格分开。
第二行，其中包含有校准点信息。校准点信息一般指95%HPD（Highest Posterior Density）对应的置信区间，单位为。尾部一定要有分号。

（3）配置文件mcmctree.ctl

seed = -1  # 表示使用系统当前时间为随机数
       seqfile = mtCDNApri123.txt # 多序列比对文件
      treefile = mtCDNApri.trees # 带有校准信息有根树
      mcmcfile = mcmc.txt #输出mcmc信息文本文件，可以用Tracer软件查看
       outfile = out.txt # 输出文件，记录了超度量树和进化速率等参数信息

         ndata = 3 # 输入的多序列比对的数据个数，这里密码子3个位置的数据；如果有一个，则设置为1
       seqtype = 0    * 0: nucleotides; 1:codons; 2:AAs #数据类型
       usedata = 1    * 0: no data; 1:seq like; 2:normal approximation; 3:out.BV (in.BV) # 设置是否利用多序列比对的数据：\
#0，表示不使用多序列比对数据，则不会进行likelihood计算，虽然能得到mcmc树且计算速度飞快，但是其分歧时间结果是有问题的；\
#1，表示使用多序列比对数据进行likelihood计算，正常进行MCMC，是一般使用的参数; \
#2，进行正常的approximation likelihood分析，此时不需要读取多序列比对数据，直接读取当前目录中的in.BV文件。该文件是使用usedata = 3参数生成的out.BV文件重命名而来的。\
#此外，由于程序BUG，当设置usedata = 2时，一定要在改行参数后加 *，否则程序报错 Error: file name empty.. \
#3，程序利用多序列比对数据调用baseml/codeml命令对数据进行分析，生成out.BV文件。由于mcmctree调用baseml/codeml进行计算的参数设置可能不太好（特别时对蛋白序列进行计算时），\
#推荐自己修该软件自动生成的baseml/codeml配置文件，然后再手动运行baseml/codeml命令，再整合其结果文件为out.BV文件。
         clock = 2    * 1: global clock; 2: independent rates; 3: correlated rates # (1) 全局分子钟,适用于近缘物种（两物种分化时间小于20个million分化时间或者序列差异小于5%即近缘物种）\
# (2) 树枝上进化速率服从独立同分布的对数正态分布（推荐使用）
# (3) 树枝上的进化速率自相关。
       RootAge = '<1.0'  * safe constraint on root age, used if no fossil for root. # 设置root节点的分歧时间，一般设置一个最大值。

         model = 0    * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
# 碱基替换模型：
  #（0）无参数（适用于近缘物种）
  #（1）参数为转换颠换比Kappa
  #（2）参数为T，C，A，G的频率
 #（3）最复杂的进化模型
 #（4）参数为T，C，A，G的频率+kappa（适用于远缘物种）
         alpha = 0.5    * alpha for gamma rates at sites 
# 核酸序列中不同位点，其进化速率不一致，其变异速率服从GAMMA分布。一般设置GAMMA分布的alpha值为0.5。\
# 若该参数值设置为0，则表示所有位点的进化速率一致。此外，当userdata = 2时，alpha、ncatG、alpha_gamma、kappa_gamma这4个GAMMA参数无效。因为userdata = 2时，不会利用到多序列比对的数据。
         ncatG = 5    * No. categories in discrete gamma # 设置离散型GAMMA分布的categories值。

     cleandata = 0    * remove sites with ambiguity data (1:yes, 0:no)?

       BDparas = 1 1 0.1  * birth, death, sampling # 设置出生率、死亡率和取样比例。若输入有根树文件中的时间单位发生改变，则需要相应修改出生率和死亡率的值。
例如，时间单位由100Myr变换为1Myr，则要设置成".01 .01 0.1"。
   kappa_gamma = 6 2      * gamma prior for kappa  #设置kappa（转换/颠换比率）的GAMMA分布参数
   alpha_gamma = 1 1      * gamma prior for alpha #设置GAMMA形状参数alpha的GAMMA分布参数

   rgene_gamma = 2 20 1   * gammaDir prior for rate for genes #设置序列中所所有位点平均[碱基/密码子/氨基酸]替换率的Dirichlet-GAMMA分布参数：
  sigma2_gamma = 1 10 1   * gammaDir prior for sigma^2     (for clock=2 or 3) #设置所有位点进化速率取对数后方差（sigma的平方）的Dirichlet-GAMMA分布参数：alpha=1、beta=10、初始方差值为1。
# 当clock参数值为1时，表示使用全局的进化速率，各分枝的进化速率没有差异，即方差为0，该参数无效；
# 当clock参数值为2时，若修改了时间单位，该参数值不需要改变；
# 当clock参数值为3时，若修改了时间单位，该参数值需要改变。

      finetune = 1: .1 .1 .1 .1 .1 .1 * auto (0 or 1): times, musigma2, rates, mixing, paras, FossilErr #冒号前的值设置是否自动进行finetune，一般设置成1，然后程序自动进行优化分析

         print = 1   * 0: no mcmc sample; 1: everything except branch rates 2: everything # 设置打印mcmc的取样信息
#0，不打印mcmc结果；
#1，打印除了分支进化速率的其它信息（各内部节点分歧时间、平均进化速率、sigma2值）；
#2，打印所有信息。 
        burnin = 2000 #将前2000次迭代burnin后，再进行取样（即打印出该次迭代计算的结果信息，各内部节点分歧时间、平均进化速率、sigma2值和各分支进化速率等）
      sampfreq = 10 #每10次迭代则取样一次。
       nsample = 20000 # 当取样次数达到该次数时，则取样结束（程序也将运行结束）。

运行

~/paml4.9j/src/mcmctree mcmctree.ctl

序在运行过程当中，屏幕也会生成一些信息。

查看每列中的值（接受比例、时间和速率）很重要。它们应该在MCMC运行期间保持稳定。如果验收比例变化太大（特别是在开始时），这意味着预烧时间不够长，因此应增加控制文件中的预烧变量

在最后，给出各个内部节点的分歧时间(t)、平均变异速率(mu)、变异速率方差(sigma2)和r_left的Posterior信息：
均值(mean)、95%双侧置信区间(95% Equal-tail CI)和95% HPD置信区间(95% HPD CI)等信息。
此外，第二列给出了各个内部节点的Posterior mean time信息，可以用于。jnode是程序multipvtime使用的节点号，由JeffThorne编写。
倒数第二列值： 比如第8节点； 0.0436=0.1850-0.1414

结果分析

运行会得到以下4个文件

• FigTree.tre 生成含有分歧时间的超度量树文件
  适用于Andrew Rambaut的figtree程序（http://tree.bio.ed.ac.uk/software/figtree/）

节点处的位分化时间，比如 chimpanzee 和bonobo分化时间为3.02 MYA。该数值与out .tee文件最后面的数值对应，也可查看95%置信区间

• mcmc.txt MCMC取样信息，包含各内部节点分歧时间、平均进化速率、sigma2值等信息，可以在Tracer软件中打开。通过查看各参数的ESS值，若ESS值大于200，则从一定程度上表示MCMC过程能收敛，结果可靠。

在我们的示例中，该文件有50列(如果print=2)或14列(如果print=1)和20002行。第一列是样本的生成(迭代)数，接下来的48列(print=2)或12列(print=1)分别对应于所分析的48个或12个参数；最后一列是似然值(第50列，如果print=2；或第14列，如果print=1)。行数与MCMC期间采集的样本数相对应。“mcmc.txt”文件适用于Tracer程序分析（http://beast.bio.ed.ac.uk/tracer）

• out.txt 包含由较多信息的结果文件。例如，各碱基频率、节点命名信息、各节点分歧时间、进化速率和进化树等

打开此文件。文件开头有很多垃圾，你通常可以忽略。向下滚动该文件，直到您看到3个系统发育树(如果print=1；如果print=2，您将看到6个)
第一棵树只包含节点标签；第二棵树包含以时间单位表示的分支长度；第三棵树包含以时间单位表示的分支长度，以及节点年龄的可信度区间（记住只有节点有）。最后三棵树（仅当您将“print”设置为2时写入输出文件）具有替换率，而不是分支长度，每个树表示每个位点（在我们的示例中，三个密码子位置中的每一个）。在树下面，是后验均值，与屏幕输出最后的一致。

• SeedUsed 文件
  文件包含用于初始化MCMC的随机种子。如果将文件中的值（在我的例子中是1847523971）复制到控制文件的seed变量中，则可以再次运行分析并获得相同的结果

注意事项

（1） 如何检测MCMC的结果是否达到收敛状态？
收敛的意思，即经过很多次迭代后，得到的MCMC树的各枝长趋于一个定值，变动非常小。
于是，最直接的检测方法是：分别使用不同的Seed值进行mcmctree ()，进行两次或多次分析，然后比较两个结果树是否一致，实际就是比较树文件中各内部节点的Height值（分歧时间 / Posterior time）。(可查看out.txt文件后面分化时间（第二列）)
可通过以下3种方法进行：

• 结果文件mcmc.txt中每行记录了一个取样的结果，包含各内部节点的Posterior time值，即进化树的Height值，对该mcmc.txt文件进行分析，得到每个内部节点的Posterior mean time值。然后作图 （官方说明文档第6页图示），横坐标表示第一次运行的各内部节点的分歧时间均值、纵坐标表示第二次运行的各内部节点的分歧时间均值。该散点图要表现出非常明显的对角线，才认为达到收敛。这时可以考虑计算其相关系数来判断是否符合线性。
• 直接比较两次结果树文件中的各枝长。计算各枝长总的偏差百分比，当偏差百分比低于0.1%，则认为两次结果非常吻合，差异低于0.1%，认为达到收敛。(如上述所述)
• 此外，还可以使用Tracer分析mcmc.txt文件，检测其ESS值，一般认为该值高于200，则可能达到收敛。该方法可用于辅助检测。最后，若不收敛，则需要提高burnin、nsample值，重新运行程序。

以下图片是根据两次运行结果的分化时间进行作图：

（2）如何设置burnin、sampfreq和nsample值？

• 一般推荐设置nsample值不小于20k(官方示例数据中设置的值)，只有当该值较大时，求得的均值才有意义。
• sampfreq表示取样间隔，一般设置为10，100或1000。nsample 和sampfreq 的乘积表示有效迭代的次数，这些迭代越准确，次数越大，则结果越好，越趋于收敛。当然，次数越大，越消耗时间。若数据量很小，可以考虑设置sampfreq为1000；若数据量大，每次迭代耗时多，则考虑设置sampfre为10
• 一般最开始的迭代结果很差，需要去除(burnin)其结果，则设置burnin值。要设置足够大的burnin值，保证在程序运行时当迭代比例达到0% (即burnin结束) 时，其参数的接受比例值较好，在30%左右且随迭代次数增多时变动幅度不大。推荐设置burnin的迭代次数为有效次数的10~40%。PAML软件时先burnin后再记录有效迭代的参数值；其它软件如BEAST则记录所有的参数值后，最后汇总时burnin掉一定比例的数据。
• 总体上，其实就是两个参数：burnin掉的迭代次数和用于计算结果的有效迭代次数。由于需要根据有效迭代数据来进行均值计算，若记录每次迭代的参数，则生成的mcmc.txt文件行数过多，文件太大，汇总时也消耗计算。于是设置没隔一定迭代轮数取样一次，这样生成的mcmc.txt文件不会太大，对最后的均值计算影响不大。

我个人推荐有效迭代次数 1~10M 次，burnin 0.2~4M次。按时间消耗从少到多的参数值：
数据简单时，进行0.5M迭代次数，burnin比例40%：{ burnin 200k; sampfreq 10; nsample 50k }
数据中等时，进行1M迭代次数，burnin比例40%：{ burnin 400k; sampfreq 10; nsample 100k }
数据复杂时，进行5M迭代次数，burnin比例20%：{ burnin 1000k; sampfreq 10; nsample 500k }
数据巨大时，进行10M迭代次数，burnin比例20%：{ burnin 2000k; sampfreq 100; nsample 100k }

使用相应参数进行分析后，若不满意其收敛程度，则选用更多迭代次数的参数。

• 带有化石校准的系统发育树

利用物种间单拷贝基因进行构建系统发育树，如果要选择核苷酸/氨基酸替换模型，则可以使用IQ-TREE,而后使用RAxML进行构建进化树。

• 两物种间的化石分化时间
  通过time tree网站,可进行物种化石分化时间查询，比如
  查询，人和黑猩猩
  
  
  查询以后会得到下面结果
  

可以看到其为： 5.1-11.8 MYA, 则将其添加到进化树中就是>.051<0.118。因为其单位为100百万年

2 使用approximate likelihood方法更快速地计算分歧时间

按照以上方法进行MCMC计算非常消耗计算。为了加快计算速度，可以将分析分成两个步骤：
（1）首先，使用Maximum Likelihood方法计算枝长和Hessian信息；
（2）再使用MCMC方法计算分歧时间。

具体步骤：
第一步：返回上一次运行的文件夹，构建Hessian的文件夹,拷贝树，比对文件，配置文件

mkdir Hessian
cp ../mcmctree.ctl ../mtCDNApri123.txt ../mtCDNApri.trees ./

打开配置文件，将usedate 设置为3，并运行

usedata = 3    * 0: no data; 1:seq like; 2:normal approximation; 3:out.BV (in.BV)
/soft/paml4.9j/src/mcmctree mcmctree.ctl

第二步：新建文件夹，拷贝树，配置文件，比对文件以及out.BV到新文件夹
上一步得到out.BV中含有各个枝的长度，gradient和Hessian的值，将out.BV更改为in.BV. 再次更改配置文件usedata=2

mkdir approx01
cp  ../mcmctree.ctl ../mtCDNApri123.txt ../mtCDNApri.trees ../out.BV ./
mv out.BV in.BV
usedata = 2    * 0: no data; 1:seq like; 2:normal approximation; 3:out.BV (in.BV)
/soft/paml4.9j/src/mcmctree mcmctree.ctl

本次mcmctree 根据上一步的gradient 和Hessian，进行估计分化时间。此外，还可以再次建立新文件夹approx02, 再次运行第二步，两次结果进行比较。

⚠️； 如果clock=1，则不能用approximate 方法。

参考

• MCMCtree——最详尽的物种进化历史分歧时间估算教程
• 使用PAML进行分歧时间计算
• 进化树专题（九）| 从DNA序列到分化时间——如何构建带分歧时间的进化树？