虽然不是学生物或者医学的，但是听了很多生命科学的报告。关于生命的起源、进化——从个体到宇宙，生命是如何从无到有，又如何实现了现在的千姿百态——是生命科学数百年来一直在研究的问题。其实总结来说，无非要抛出常见的三联问：
记得以前听一个老师的报告，问过大家一个问题：如果你碰到一个无所不知，文明远远高于我们的外星人，你可以问他一个问题，会问什么？我记得一个深刻的答案是：会问人类是不是人工智能？ 好像有个科幻小说里面也有这方面的描述，问题例如：宇宙的目的是什么？意识的目的是什么？灵魂的目的是什么？等等等。
利用进化树将穿越数万年的物种进化浓缩于一图，更是生物学研究中不可或缺的必备工具，也是我们领域天天打交道的事情。
由于不同的基因或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自己写个脚本转化成输入需求就行了。
本文使用 文章同步助手 同步
