mcmctree进行分化时间的估算

以下内容参考了博客https://www.cnblogs.com/bio-mary/p/12818888.html、前研究组师妹的总结资料和MCMCTree tutorials。

"MCMCTree performs Bayesian estimation of species divergence times using soft fossil constraints under various molecular clock models."

MCMCTree是PAML包的一部分，功能是在多种分子钟模型下，利用较宽松的的化石约束，用贝叶斯方法估算物种分化时间。

安装MCMCTree

在Linux系统下的安装很简单，用conda安装paml就可以：

conda install -c bioconda paml

准备文件

"The program uses for input a sequence alighment (nucleotide or protein), a phylogenetic tree with fossil calibrations, and a control file (ususally called mcmctree.ctl) that contains the instructions for the program. "

"MCMCTree is part of the PAML package."

运行mcmctree需要准备三个文件：比对序列文件、树文件和控制文件。

在paml软件包下有个example文件夹，里面有很多示例文件，用于mcmctree分析的示例文件在DatingSoftBound里面。序列文件名称：mtCDNApri123，树文件名称：mtCDNApri.trees，控制文件的名称mcmctree.ctl，该示例文件分析的是7个类人猿物种的线粒体蛋白编码基因。

1 比对序列文件（核苷酸或蛋白质）

我用的比对好的序列文件是.phy文件，在软件Geneious里完成的文件格式转化，先导入fasat文件，导出时弹出的对话框phylip alignment export中，选择Relaxed phylip-Full length names followed by a single space，意思是导出的文件中每条序列的序列名和序列之间有一个空格。但是mcmctree要求序列名和序列之间至少两个空格空格，如果序列少的化，可以手动增加空格。我的序列比较多，用的方法是：使用geneious里batch rename功能统一给序列名后面加个序列名里没有的字符，然后导出的phylip文件用编辑器打开，用空格替换那个字符。

下图是示例文件的展示，格式是txt。7表示的是物种树，3331是比对序列的长度。在示例文件中，比对序列被分成了3部分，对应第一、第二和第三密码子位点， 每一部分都如下所示。

序列文件示例

2 有化石校准信息的树文件

需要注意的是树文件必须是定根的、不带有枝长信息的二歧树（rooted binary tree）

paml软件包下面的example文件夹下的树文件（.tree）如下：

树文件示例

7代表的是物种树，1代表只有一颗树。蓝色框中显示的是化石时间标记，表示的是人（humna）和黑猩猩（chimpanzee，bonobo）的共同祖先出现的时间在0.06到0.08个100个百万年之间。

（注意化石时间标记的单位是100个百万年）

树文件代表的树的拓扑结构如下：

示例树文件的拓扑结构

我们自己的树文件常常是带枝长信息，想要去除枝长信息的话可以用notepad等编辑器打开树文件，手动删除枝长信息，还有一种简单的方法，在linux系统中运行下面的命令行，tree.nwk是输入的带枝长信息的树文件，output_tree.nwk是输出的不带枝长信息的树文件

perl -ne '$_=~s/:[\d\.]+//g; print $_;' tree.nwk > output_tree.nwk

3 控制文件

控制文件常命名为mcmctree.ctl，含有对程序的指令。

控制文件示例

各类型参数的含义如下：

1）输入输出参数

seed=-1 ：

#设置一个随机数（正整数或负整数）来运行程序，若设置为-1，表示使用系统当前时间为随机数，因此每次运行mcmctree得到的结果会有所不同。

建议用seed=-1运行程序至少两次，检验两次运行的结果是否相近。如果每次运行的结果差别很大，可能是多种原因导致的：缓慢收敛（slow convergence）、不良混合（poor mixing）、采样不完全（insufficient samples taken）, 或程序错误。

如果需要一个可重复的结果，可以设置一个特定的奇数或偶数

seqfile：比对序列文件名

treefile：树文件文件名

mcmcfile：针对设定的参数运行的MCMC的报告文件，可以用Tracer软件查看。

outfile：一旦程序运行完成，总结性的结果文件会写入该文件，该文件记录了超度量树和进化速率等参数信息。

2）数据使用说明参数

ndata：比对序列的分区数，在本例中，根据核苷酸在密码子中的位置（第一位、第二位和第三位）分成了三个分区，因此ndata=3。

usedata：设置是否利用多序列比对的数据

usedata = 0：不适用多序列比对数据，不会进行likelihood计算，虽然能得到mcmc树且计算速度飞快，但分歧时间结果有问题

usedata = 1 使用多序列数据进行likelihood计算，正常进行MCMC，是一般的使用参数

usedata = 2 进行正常的approximation likelihood分析，此时不需要读取多序列比对数据，直接读取当前目录中的in.BV文件，该文件是使用usedata=3参数生成的out.BV文件重命名而来。

usedata = 3：程序利用多序列比对数据调用baseml/codeml命令对数据进行分析，生成out.BV文件。

由于mcmctree调用baseml/codeml进行计算的参数设置可能不太好，尤其是蛋白序列，推荐自己修改该软件自动生成的baseml/codeml配置文件，然后手动运行baseml/codeml命令，再整合其结果文件为out.BV文件。

seqtype：比对序列的类型，=0代表核苷酸序列，=1代表密码子序列，=2代表氨基酸序列

clock：选择使用的分子钟模型，

clock=1：global clock的方法，假设所有分支的进化速率一致

clock=2：使用独立速率模型（independent rates model），在该模型中，速率符合对数-正态分布（也就是说速率的对数符合正态分布）

clock=3： 使用correlated rates的方法，和方法类似，但是log（r）的方差和时间（t）相关

RootAge：如果我们在树文件中没有对根提供时间校准，那么用参数对根提供一个时间校准。在示例文件中使用了<1.0，意思是所有猿的最近共同祖先出现的时间不会大于100个百万年前。

cleandata:设置是否移除不明确的字符或gap后再进行数据分析。

cleandata = 0不需要

cleandata = 1 需要

3）位点替换模型参数

model: 设置要使用的替换模型

model = 0: JC69,计算很快

model = 1: K80,

model = 2: F81,

model = 3: F84,

model = 4: HKY85

model = 7: GTR

alpha：核酸序列中不同位点，其进化速率不一致，其变异速率服从Gamma分布，一般设置gamma分布的alpha值为0.5：alpha = 0.5，若该参数设置为0，则表示所有位点的进化速率一致。

ncatG = 5 ：设置离散型GAMMA分布的categories值

BDparas= 1 1 0.1： 控制出生-死亡过程的参数，设置出生率、死亡率和取样比例。

若输入有根树文件中的时间单位发生变化，则需相应修改出生率和死亡率的值。例如时间单位由100Myr变换成1myr，则要设置成：.01.01 .01

kappa_gamma = 6 2 ：设置替换模型参数kappa（转换/颠换比率）的GAMMA先验。

alpha_gamma = 1 1： 设置替换模型中gamma形状参数（用于反应位点上不同的速率）alpha的GAMMA先验。

注意：当usedata = 2时，alpha、ncatG、alpha_gamma、kappa_gamma四个GAMMA参数无效，因为此时不会用到多序列比对的数据

3）进化速率参数

rgene_gamma = 2 20 1：设置序列中所所有位点平均（碱基/密码子/氨基酸）替换率的Dirichlet-GAMMA分布参数：

alpha=2、beta=20、初始平均替换率为每100million年（取决于输入有根树文件中的时间单位）1个替换。若时间单位由100Myr变换为1Myr，则要设置成"2 2000 1"。总体上的平均进化速率为：2 / 20 = 0.1 个替换 / 每100Myr，即每个位点每年的替换数为 1e-9。

sigma2_gamma = 1 10 1：设置所有位点进化速率取对数后方差（sigma的平方）的Dirichlet-GAMMA分布参数：alpha=1、beta=10、初始方差值为1。

当clock参数值为1时，表示使用全局的进化速率，各分枝的进化速率没有差异，即方差为0，该参数无效；

当clock参数值为2时，若修改了时间单位，该参数值不需要改变；

当clock参数值为3时，若修改了时间单位，该参数值需要改变。

finetune = 1: .1 .1 .1 .1 .1 .1：冒号前的值设置是否自动进行finetune，一般设置成1，然后程序自动进行优化分析；冒号后面设置各个参数的步进值：times, musigma2, rates, mixing, paras, FossilErr。由于有了自动设置，该参数不像以前版本那么重要了，可能以后会取消该参数。

4）mcmc参数

print：设置打印mcmc的取样信息

print = 0：不打印mcmc结果

print = 1：打印除了分支进化速率的其他信息（各内部节点分歧时间、平均进化速率，sigma2值）

print = 2：打印所有信息

burnin = 2000 ：将前2000次迭代burnin后，再进行取样（即打印出打印出该次迭代计算的结果信息，各内部节点分歧时间、平均进化速率，sigma2值和各分支进化速率等）

samplefreq = 10 ：每10次迭代取样一次

nsample = 20000 ：当取样次数达到该次数时，则取样结束（程序也运行结束）

burnin= 2000，samplefreq= 10，nsample= 20000：总的意思是程序会去除掉（burnin)前2000次迭代的结果，然后每10次迭代取一次样，直到取样次数达到20000次，因此MCMC会运算2000+10*20000次迭代。一般来说，程序需要进行10000-20000次取样，才能获得较好的数据总结。一般来说，如果想通过增加MCMC长度来提高收敛效果，一般是通过增加samplefreq，但保持nsample在一个比较合理的范围。

别人的博客曾进行过一个其个人的总结，照抄如下：

数据简单时，进行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

在Linux系统下运行时，将序列文件、树文件和控制文件放在同一个路径下，在该路径下运行程序：

which mcmctree  #查找mcmctree的路径

~/miniconda3/bin/mcmctree mcmctree.ctl  #~/miniconda3/bin/mcmctree是程序的绝对路径，mcmctree.ctl是控制文件的名字

用approximate likelihood calculation进行物种分化时间的估算

在对较大的数据进行likelihood function计算时，计算成本昂贵，耗时长，下面介绍针对较大比对数据的approximate likelihood计算的方法。

总共分两步，在第一步中，枝长、gradient、Hessian由最大似然法估算出，在第二步中，利用gradient、Hessian进行分化时间的估算，这一步用Taylor expansion的方法创造出近似于最大似然法的功能。

总的操作步骤如下：

新建一个文件夹Hessian，将树文件、序列文件和控制文件拷贝过去，控制文件中的usedata = 1修改为usedata = 3，然后运行程序。

然后再建立一个新文件夹approx01，将上一步中的out.BV文件、树文件、序列文件和控制文件拷贝过去