【比较基因组】从进化树到分化时间

。虽然不是学生物或者医学的，但是听了很多生命科学的报告。关于生命的起源、进化——从个体到宇宙，生命是如何从无到有，又如何实现了现在的千姿百态——是生命科学数百年来一直在研究的问题。其实总结来说，无非要抛出常见的三联问：

记得以前听一个老师的报告，问过大家一个问题：如果你碰到一个无所不知，文明远远高于我们的外星人，你可以问他一个问题，会问什么？我记得一个深刻的答案是：会问人类是不是人工智能？   好像有个科幻小说里面也有这方面的描述，问题例如：宇宙的目的是什么？意识的目的是什么？灵魂的目的是什么？等等等。

利用进化树将穿越数万年的物种进化浓缩于一图，更是生物学研究中不可或缺的必备工具，也是我们领域天天打交道的事情。

由于不同的基因或DNA片段的进化速率存在较大的差异，我们就可以利用这些基因或DNA片段来估算各个分类水平上有机体之间的进化关系。

分化时间是当前宏观进化分析的一个热点，以某一个或某几个特定类群的化石时间作为参照，通过基因序列间的分歧程度以及分子钟来估计分支间的分歧时间，同时计算系统发育树上其它节点的发生时间，从而推断相关类群的起源和不同类群的分歧时间。

 

根据分子钟理论，基因序列中密码子随着时间的推移以几乎恒定的比例相互替换，即具有恒定的演化速率，两个物种之间的遗传距离将与物种的分化时间成正比。我们通常采用单拷贝基因家族中的四重简并位点来估算分子钟（替换速率）以及物种间的分化时间，密码子中的四重简并位点由于第三碱基不改变所编码的氨基酸，属于中性进化，其中中性替换速率一般用每个位点每年的变异数来衡量。

 

利用贝叶斯统计或其他方法估计物种分歧时间的程序很多，如R8S、MCMCTREE、MULTIDIVTIME、BEAST等，通过不同的策略将化石时间信息整合到一个系统发育树中，从而计算得到Divergence time Tree。

 

下面就是一个带有分化时间的进化树（分支带有物种的分化时间，并添加了分化时间的置信度范围，蓝色数字为分化时间）：

=====构建系统发育树=====

前面讲过很多如何构建发育树的文章，比如基于orthologous的基因利用muscle/Gblocks/RAxml构建发育树。

muscle -in ../OG0012924.fa -out OG0012924.fa

Gblocks OG0012924.fa -t=p -b4=5 -b5=h -e=.gb

raxmlHPC-PTHREADS-SSE3  ­-f a ­-x 12345 ­-p 12345 ­-# 100 ­-m PROTGAMMALGX ­-s OG0012924.fa.gb ­-n ex -T 20 

下面我们测试用不同的几个工具尝试去计算分化时间：

=======两两物种的分化时间=======

通过网站：http://www.timetree.org/ ，可以查两两物种间的分化时间。有个单位是：1Mya  = 1 million years ago。

 

另外一个可以查询分化时间的网站：https://fossilcalibrations.org/

==================r8s=======================

准备输入文件：

前面是输入的传统方式的tree。

BLFORMAT

length     设置树的枝长单位。若枝长单位是位点替换率，则其值为 persize，则枝长 * 序列长度 = 替换数；若枝长单位是替换数，则该参数值为 total。默认参数是 total。

nsites     设置多序列比对的序列长度。

ultrametric 是否是超度量树，一般进化树选 no。默认参数是 no。

MRCA

该命令用来设置内部节点名。示例中设置了 tree 和 scer 的 most recent common ancestor 的节点名为 asc。

FIXAGE

该命令用来设置 MRCA 指定的节点的分歧时间，使用该指定的分歧时间作为校正，来预测其它节点的分歧时间。

r8s 需要至少有一个内部节点进行 fixage。

CONSTRAIN

该命令用来约束 MRCA 指定的节点的分歧时间，设置该节点允许的最大或最小分歧时间。

DIVTIME

该命令用来进行分歧时间和替换速率计算。总共有 4 种计算方法（LF | NPRS | PL）和 3 种数学算法（Powell | TN | QNEWT）。一般常用且较优，是使用 PL 和 TN。

但是使用 PL 方法，则需要设置参数 smoothing 的值。通过设置多个 smoothing 的值来进行一些计算，选择最优的值即可。一般情况下，该值位于 1~1000 能得到一个最佳（最低）的分值。

divtime method=PL crossv=yes cvstart=0cvinc=1 cvnum=18;

上述命令，是设置 smoothing 的值从 1, 10, 100, 1000 ... 1e17, 来计算，最后得到最佳的 smoothing值。

若使用 fixage 对 2 个节点的分歧时间进行了固定，则可以运行命令：

divtime method=PL crossv=yesfossilfixed=yes cvstart=0 cvinc=1 cvnum=18;

若使用 fixage 对 1 个节点进行分歧时间固定，同时使用 constrain 对 2 个节点进行了约束，则可以运行命令：

divtime method=PL crossv=yesfossilconstrained=yes cvstart=0 cvinc=1 cvnum=18;

得到最优的 smoothing 值后，使用 set 进行设置，然后运行 divtime 进行分歧时间和替换速率的计算。

SHOWAGE

使用该命令能得到各个节点的分歧时间和替换速率。

DESCRIBE

使用该命令得到进化树的图和文字结果。其 plot 参数如下：

cladogram   得到分支树的图，图上有各个节点的编号，和 showage 的结果结合观察。

phylogram   得到进化树的图，枝长表示替换数。

chronogram  得到超度量树的图，枝长表示时间。

ratogram    得到树图，枝长表示替换速率。

phylo_description     得到树的ASCII文字结果，枝长表示替换数。

chrono_description    得到树的ASCII文字结果，枝长表示时间。

rato_description      得到树的ASCII文字结果，枝长表示替换速率。

node_info   得到节点的信息表格

然后运行输入文件：/gpfs03/home/jingjing/software/r8s1.81/src/r8s -b -f r8s_ctl_file.txt > r8s_tmp.txt

=================MCMCtree==================

使用 PAML mcmctree（http://abacus.gene.ucl.ac.uk/software/paml.html） 估计分歧时间。

下载安装：

wget http://abacus.gene.ucl.ac.uk/software/paml4.9j.tgz

tar xf paml4.9j.tgz

cd paml4.9j

rm bin/*.exe

cd src

make -f Makefile

rm *.o

mv baseml basemlg codeml pamp mcmctree evolver yn00chi2  ../bin

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

序列比对文件（nucl或者pep都可以）

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

mcmctree.ctl；控制文件

(1) 序列比对文件

paml4.9j/examples/DatingSoftBound/mtCDNApri123.txt，是一个phlip格式，PAML要求输入的Phylip格式，其物种名和后面的序列之间至少间隔两个空格（是为了允许物种名的属名和种名之间有一个空格）。

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

 

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

 

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

示例文件如下：

第一行，7个物种，1颗树，中间用空格分开，并且不能有枝长。

第二行，其中包含有校准点信息。校准点信息一般指95%HPD（Highest Posterior Density）对应的置信区间，单位为100个百万年(MY)。尾部一定要有分号。

这棵树有两种化石校准，一种是人类黑猩猩的最近共同祖先：’>.06<.08’，另一种是大猩猩的最近共同祖先：’>.12<.16’。时间单位是100百万年(Myr)。因此，第一次校准将人类黑猩猩的共同祖先限制在6-8 myr之前。mcmctree中的校准界限是软的，并且很小的概率（默认为0.025）可能违反界限。注意，树的根部没有化石校准。MCMCTree需要在树的根上进行校准，当该校准不存在于树文件中时，必须在控制文件中指定它(variable RootAge)。

(3) 配置文件：

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:normalapproximation; 3:out.BV                         (in.BV) # 设置是否利用多序列比对的数据：

#0，表示不使用多序列比对数据，则不会进行likelihood计算，虽然能得到mcmc树且计算速度飞快，但是其分歧时间结果是有问题的；

#1，表示使用多序列比对数据进行likelihood计算，正常进行MCMC，是一般使用的参数; 

#2，进行正常的approximationlikelihood分析，此时不需要读取多序列比对数据，直接读取当前目录中的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 nofossil 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 (0or 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 # 当取样次数达到该次数时，则取样结束（程序也将运行结束）。

(4) 运行：

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信息，可以用于收敛分析。jnode是程序multipvtime使用的节点号，由JeffThorne编写。

倒数第二列值：比如第8节点；0.0436=0.1850-0.1414。

 

(5) 结果分析：

运行会得到以下4个文件

 

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

节点处的位分化时间，比如 chimpanzee 和bonobo分化时间为3.01 MYA。

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 包含由较多信息的结果文件。例如，各碱基频率、节点命名信息、各节点分歧时间、进化速率和进化树等。

所以，我们只需要把比对文件和生成的tree自己写个脚本转化成输入需求就行了。

本文使用 文章同步助手 同步