比较基因组学基因组学遗传学

mcmctree估算物种分歧时间

2021-01-10  本文已影响0人  斩毛毛

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

本文就利用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个文件:

(1)序列比对文件

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



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

此外,也可以输入密码子或氨基酸序列比对信息,则需要再配置文件中修改相应的参数来表明数据类型 \color{red}{ndata}。若使用氨基酸序列来进行分析,由于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)对应的置信区间,单位为\color{red}{100个百万年(MY)}。尾部一定要有分号。

(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

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


第一列:取样的百分比进度。
第2~6列: 参数的接受比例。一般,其值应该在30%左右。20~40%是很好的结果, 15~70%是可以接受的范围。若这些值在开始时变动较大,则表示burnin数设置太小。
第7~x列:各内部节点的平均分歧时间,第7列则是root节点的平均分歧时间。若有y个物种,则总共有y-1个内部节点,从第7列开始后的y-1列对应这些内部节点。
倒数第3~x列:r_left值。若输入3各多序列比对结果,则有3列。x列的前一列是中划线 - 。
倒数第1~2列:likelihood值和时间消耗。

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

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

结果分析

运行会得到以下4个文件

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

在我们的示例中,该文件有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

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

注意事项

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

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


(2)如何设置burnin、sampfreq和nsample值?

我个人推荐有效迭代次数 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进行构建进化树。

可以看到其为: 5.1-11.8 MYA, 则将其添加到进化树中就是>.51<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 方法。

参考

上一篇 下一篇

猜你喜欢

热点阅读