PAML软件中的mcmctree命令估算物种分歧时间 PAML软件中的mcmctree命令可以使⽤Bayesian⽅法估算物种分歧时间。
对程序输⼊带有校准点信息的有根树、多序列⽐对结果，即可得到进化树各内部节点95%置信区间分歧时间信息。南京海关
PAML软件也能利⽤带有枝长信息的树，快速进⾏分歧时间计算。
1. 输⼊⽂件（1/3）：带有校准点的有根树 s 输⼊⼀个带有校准点的有根树⽂件，⽰例如下：
氧化铟
7 1
((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);
⽂件内容分两⾏：第⼀⾏表述树中有7个物种，共计1个树，两个数值之间⽤空分割；
第⼆⾏则是Newick格式树信息，其中包含有校准点信息。校准点信息⼀般指95%HPD（Highest Posterior Density）对应的置信区间；
校准点单位是100MYA（软件说明⽂档中使⽤该单位，也推荐使⽤该单位，若使⽤其它单位，后续配置⽂件中的相关参数也需要对应修改）。此外，Newick格式的树尾部⼀定要有分号，没有的话程序可能不能正常运⾏。 2. 输⼊⽂件（2/3）：密码⼦在3个位点的多序列⽐对结果（Phylip格式）
7 3331
human
chimpanzee
bonobo
gorilla
orangutan
sumatran
gibbon
7 3331
human
chimpanzee
bonobo
gorilla
orangutan
sumatran
gibbon
7 3331
human
chimpanzee
bonobo
gorilla
orangutan
sumatran
gibbon 贵州民族大学教务管理系统
PAML要求输⼊的Phylip格式，其物种名和后⾯的序列之间⾄少间隔两个空格（是为了允许物种名的属名和种名之间有⼀个空格）。
当然，也可以输⼊⼀个多序列⽐对结果进⾏分析。在软件说明⽂档中，是密码⼦三个位点分别的多序列⽐对结果，程序可以分别计算密码⼦各个位点的变异速率。
此外，也可以输⼊密码⼦或氨基酸序列⽐对信息，则需要再配置⽂件中修改相应的参数来表明数据类型。若使⽤氨基酸序列来进⾏分析，由于mcmctree命令不能选择较好的氨基酸替换模型进⾏分析，需要⾃⼰⼿动运⾏codeml进⾏分析后，在⽣成中间⽂件⽤于运⾏mcmctree，⽐较⿇烦。因此，推荐按上述⽅法使⽤密码⼦三位点的碱基序列作为输⼊信息进⾏PAML分析。
3. 输⼊⽂件（3/3）：mcmctree命令的配置⽂件 l
软件难点在于配置⽂件的理解，⽰例如下：
*输⼊输出参数：
seed = -1 *设置⼀个随机数来运⾏程序。若设置为-1，表⽰使⽤系统当前时间为随机数。
seqfile = *设置输⼊的多序列⽐对⽂件路径
treefile = s *设置输⼊的带校准点信息有根树的⽂件路径
mcmcfile = *设置输出的mcmc信息⽂本⽂件，可以⽤Tracer软件查看
outfile = *设置输出⽂件路径，该⽂件中记录了超度量树和进化速率等参数信息。
*数据使⽤说明参数：
ndata = 3 *设置输⼊的多序列⽐对的数据个数。这⾥指密码⼦3个位置的数据。
seqtype = 0 *设置多序列⽐对的数据类型：
0，表⽰核酸数据；
1，表⽰密码⼦⽐对数据；
2，表⽰氨基酸数据。
根据不同的数据类型，程序调⽤相应的baseml或codeml进⾏相应的参数计算。
usedata = 1 *设置是否利⽤多序列⽐对的数据：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⽂件。
cleandata = 0 *设置是否移除不明确字符（N、？、W、R和Y等）或含以后gap的列后再进⾏数据分析：
0，不需要，但在序列两两⽐较的时候，还是会去除后进⾏⽐较；
1，需要。
clock = 2 *设置分⼦钟⽅法：
1，global clock⽅法，表⽰所有分⽀进化速率⼀致；
2，independent rates⽅法，各分⽀的进化速率独⽴且进化速率的对数log(r)符合正态分布;
3，correlated rates⽅法，和⽅法2类似，但是log(r)的⽅差和时间t相关。
* TipDate = 1 100 *当外部节点由取样时间时使⽤该参数进⾏设置，同时该参数也设置了时间单位。
2013小企业会计准则具体数据⽰例请见examples/TipData⽂件夹。
RootAge = '<1.0' *设置root节点的分歧时间，⼀般设置⼀个最⼤值。
*位点替换模型参数：
model = 4 *设置碱基替换模型：0，JC69；1，K80；2，F81；3，F84；4，HKY85。
当设置usedata = 1时，不能使⽤其它碱基替换模型。
若需要选择其它碱基替换模型，则考虑使⽤usedata = 3参数，就可以设置model参数值为其它的碱基替换模型。
此时，程序会将数据进⾏分割，例如，ndata = 3，则程序分别对3个数据进⾏分析，⽣成baseml的输⼊⽂件，
并调⽤baseml进⾏分析（若数据量较⼤，分析较慢，推荐按ctrl + c 强⾏终⽌，程序则终⽌baseml的运⾏，进⾏下⼀个数据的输⼊⽂件⽣成并运⾏下个数据的baseml命令，
再按 ctrl + c 强⾏终⽌，这样可以让程序⽣成3各数据的输⼊⽂件和baseml配置⽂件，然后⼿动并⾏化运⾏，加快时间）；
每个分析结果中会得到rst2结果⽂件；可以在此处⼿动修改l配置⽂件，例如，修改其model参数值，推荐修改outfile参数值，从⽽可以⼿动并⾏化运⾏baseml命令；
然后将所有数据的rst2结果使⽤cat命令合并成⽂件in.BV，⽤于下⼀步设置参数usedata = 2，进⾏MCMC分析。
alpha = 0.5 *核酸序列中不同位点，其进化速率不⼀致，其变异速率服从GAMMA分布。⼀般设置GAMMA分布的alpha值为0.5。
若该参数值设置为0，则表⽰所有位点的进化速率⼀致。
伽玛分布（Gamma Distribution）是统计学的⼀种连续概率函数，是概率统计中⼀种⾮常重要的分布。“指数分布”和“χ2分布”都是伽马分布的特例。 Gamma分布中的参数α称为形状参数（shape parameter）Gamma分布的特殊形式，当形状参数α=1时，伽马分布就是参数为γ的指数分布，X~Exp（γ）当α=n/2，β=1/2时，伽马分布就是⾃由度为n的卡⽅分布，X^2(n)
此外，当userdata = 2时，alpha、ncatG、alpha_gamma、kappa_gamma这4个GAMMA参数⽆效。
因为userdata = 2时，不会利⽤到多序列⽐对的数据。
ncatG = 5 *设置离散型GAMMA分布的categories值。
洛桑多吉
BDparas = 1 1 0.1 *设置出⽣率、死亡率和取样⽐例。若输⼊有根树⽂件中的时间单位发⽣改变，则需要相应修改出⽣率和死亡率的值。
例如，时间单位由100Myr变换为1Myr，则要设置成".01 .01 0.1"。
kappa_gamma = 6 2 *设置kappa（转换/颠换⽐率）的GAMMA分布参数。
alpha_gamma = 1 1 *设置GAMMA形状参数alpha的GAMMA分布参数。
*进化速率参数：
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。
由于有了⾃动设置，该参数不像以前版本那么重要了，可能以后会取消该参数。
*MCMC参数：
print = 1 *设置打印mcmc的取样信息：
0，不打印mcmc结果；
1，打印除了分⽀进化速率的其它信息（各内部节点分歧时间、平均进化速率、sigma2值）；
2，打印所有信息。
burnin = 2000 *将前2000次迭代burnin后，再进⾏取样（即打印出该次迭代计算的结果信息，各内部节点分歧时间、平均进化速率、sigma2值和各分⽀进化速率等）。
sampfreq = 10 *每10次迭代则取样⼀次。
nsample = 20000 *当取样次数达到该次数时，则取样结束（程序也将运⾏结束）。
运⾏mcmctree命令进⾏分歧时间计算
程序运⾏命令：
l
结果⽂件解释
程序在运⾏过程中，会在屏幕⽣⽣成⼀些信息。
⽐较耗时间的步骤主要在于取样的百分⽐进度：
第⼀列：取样的百分⽐进度。
第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值和时间消耗。
屏幕信息最后，给出各个内部节点的分歧时间(t)、平均变异速率(mu)、变异速率⽅差(sigma2)和r_left的Posterior信息：
均值(mean)、95%双侧置信区间(95% Equal-tail CI)和95% HPD置信区间(95% HPD CI)等信息。
此外，倒数第⼆列给出了各个内部节点的Posterior mean time信息，可以⽤于收敛分析。
在当前⽬录下，⽣成⼏个主要结果：
< ⽣成含有分歧时间的超度量树⽂件正交试验设计
< MCMC取样信息，包含各内部节点分歧时间、平均进化速率、sigma2值等信息，可以在Tracer软件中打开。
通过查看各参数的ESS值，若ESS值⼤于200，则从⼀定程度上表⽰MCMC过程能收敛，结果可靠。
< 包含由较多信息的结果⽂件。例如，各碱基频率、节点命名信息、各节点分歧时间、进化速率和进化树等。
使⽤approximate likelihood⽅法更快速地计算分歧时间
按照以上⽅法进⾏MCMC计算⾮常消耗计算。
为了加快计算速度，可以将分析分成两个步骤：
（1）⾸先，使⽤Maximum Likelihood⽅法计算枝长和Hessian信息；
（2）再使⽤MCMC⽅法计算分歧时间。
运⾏步骤也相应分成两步：
（1）设置配置⽂件中参数 usedata = 3，然后和上⾯同样运⾏mcmctree命令，软件会对各个多序列⽐对结果调⽤baseml/codeml命令进⾏分歧时间计算，⽣成结果⽂件rst2，将各个多序列⽐对结果⽂件内容直接cat（linux系统的常⽤命令）合并成out.BV⽂件；
（2）然后将out.BV⽂件重命名为in.BV，并将配置⽂件中参数 usedata = 2，再次运⾏同样的mcmctree命令，程序会运⾏MCMC分析并得到结果。使⽤该⽅法，其MCMC速度快很多。
若输⼊数据是protein序列，则mcmctree调⽤codeml命令进⾏分析时的model参数选择是不正确的，需要⾃⼰⼿动修改codeml命令的参数配置⽂件后再运⾏codeml命令，并将其结果⽂件rst2的内容合并到out.BV⽂件中。
若想使⽤更复杂的模型，如GTR进⾏mcmctree分析，则按同样⽅法⼿动修改配置⽂件并允许baseml命令，最后将rst2的内容整合到out.BV⽂件中。此外，可以将分别来⾃核酸(⾮编码RNA)和蛋⽩的rst2内容同时整合到out.BV中进⾏计算。
使⽤infinitesites输⼊含有枝长信息的系统进化树快速进⾏分歧时间计算
infinitesites程序进⾏分歧时间计算时，其假设前提为多序列⽐对的序列长度为⽆限长。
相⽐于正常的mcmctree命令，要求额外多输⼊⼀个名为 的⽂件。该⽂件中存放了相应的带有枝长信息的系统发育树。该⽂件中系统发育树的拓扑结构要和s⼀致。推荐使⽤baseml命令输⼊s和多序列⽐对结果计算枝长。
< ⽂件内容⽰例：
7
((((human: 0.029043, (chimpanzee: 0.014557, bonobo: 0.010908): 0.016729): 0.015344, gorilla: 0.033888): 0.033816, (orangutan: 0.026872, sumatran: 0.022437): 0.069648): 0.073309, gibbon: 0.024637); ((((human: 0.012463, (chimpanzee: 0.002782, bonobo: 0.003835): 0.003331): 0.004490, gorilla: 0.014278): 0.006308, (orangutan: 0.010818, sumatran: 0.008845): 0.030551): 0.004363, gibbon: 0.029246); ((((human: 0.270862, (chimpanzee: 0.066698, bonobo: 0.056883): 0.124104): 0.139082, gorilla: 0.310797): 0.391342, (orangutan: 0.152555, sumatran: 0.114176): 0.696518): 0.017607, gibbon: 1.394718)；运⾏infinitesites命令进⾏分歧时间计算：
l
使⽤infinitesites进⾏分歧时间计算时，程序要求输⼊多序列⽐对⽂件。虽然程序读取了序列信息，但在计算时会忽略其序列信息。
其它注意事项
(1) 如何检测MCMC的结果是否达到收敛状态？
收敛的意思，即经过很多次迭代后，得到的MCMC树的各枝长趋于⼀个定值，变动⾮常⼩。
于是，最直接的检测⽅法是：分别使⽤不同的Seed值进⾏mcmctree或infinitesites进⾏两次或多次分析，然后⽐较两个结果树是否⼀致，实际就是⽐较树⽂件中各内部节点的Height值（分歧时间 / Posterior time）。
⽐较⽅法可以有多种:
(a) 结果⽂件中每⾏记录了⼀个取样的结果，包含各内部节点的Posterior time值，即进化树的Height值，对该⽂件进⾏分析，得到每个内部节点的Posterior mean time值。然后作图（官⽅说明⽂档第6页图⽰），横坐标表⽰第⼀次运⾏的各内部节点的分歧时间均值、纵坐标表⽰第⼆次运⾏的各内部节点的分歧时间均值。该散点图要表现出⾮常明显的对⾓线，才认为达到收敛。这时可以考虑计算其相关系数来判断是否符合线性。
(b) 我觉得第⼀种官⽅⽂档给的⽅法⽐较⿇烦且是否符合对⾓线没有明确的定义，于是就直接⽐较两次结果树⽂件中的各枝长。计算各枝长总的偏差百分⽐，当偏差百分⽐低于0.1%，则认为两次结果⾮常
吻合，差异低于0.1%，认为达到收敛。
(c) 此外，还可以使⽤Tracer分析⽂件，检测其ESS值，⼀般认为该值⾼于200，则可能达到收敛。该⽅法可⽤于辅助检测。最后，若不收敛，则需要提⾼burnin、nsample值，重新运⾏程序。
(2) 如何设置burnin、sampfreq和nsample值？
(a) ⼀般推荐设置nsample值不⼩于20k(官⽅⽰例数据中设置的值)，只有当该值较⼤时，求得的均值才有意义。
(b) sampfreq表⽰取样间隔，⼀般设置为10，100或1000。nsample 和sampfreq 的乘积表⽰有效迭代的次数，这些迭代越准确，次数越⼤，则结果越好，越趋于收敛。当然，次数越⼤，越消耗时间。若数据量很⼩，可以考虑设置sampfreq为1000；若数据量⼤，每次迭代耗时多，则考虑设置sampfre为10。
(c) ⼀般最开始的迭代结果很差，需要去除(burnin)其结果，则设置burnin值。要设置⾜够⼤的burnin值，保证在程序运⾏时当迭代⽐例达到0% (即burnin结束) 时，其参数的接受⽐例值较好，在30%左右且随迭代次数增多时变动幅度不⼤。推荐设置burnin的迭代次数为有效次数的10~40%。PAML软件时先burnin后再记录有效迭代的参数值；其它软件如BEAST则记录所有的参数值后，最后汇总时burnin掉⼀定⽐例的数据。
(d) 总体上，其实就是两个参数：burnin掉的迭代次数和⽤于计算结果的有效迭代次数。由于需要根据有效迭代数据来进⾏均值计算，若记录每次迭代的参数，则⽣成的⽂件⾏数过多，⽂件太⼤，汇总时也消耗计算。于是设置没隔⼀定迭代轮数取样⼀次，这样⽣成的⽂件不会太⼤，对最后的均值计算影响不⼤。(e) 我个⼈推荐有效迭代次数 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 }
使⽤相应参数进⾏分析后，若不满意其收敛程度，则选⽤更多迭代次数的参数。
从DNA序列到分化时间——如何构建带分歧时间的进化树
分⽀带有物种的分化时间，有些还可添加分化时间的置信度范围
分化时间是当前宏观进化分析的⼀个热点，以某⼀个或某⼏个特定类的化⽯时间作为参照，通过基因序列间的分歧程度以及分⼦钟来估计分⽀间的分歧时间，同时计算系统发育树上其它节点的发⽣时间，从⽽推断相关类的起源和不同类的分歧时间。
根据分⼦钟理论，基因序列中密码⼦随着时间的推移以⼏乎恒定的⽐例相互替换，即具有恒定的演化速率，两个物种之间的遗传距离将与物种的分化时间成正⽐。我们通常采⽤单拷贝基因家族中的四重简并位点来估算分⼦钟（替换速率）以及物种间的分化时间，密码⼦中的四重简并位点由于第三碱基不改变所编码的氨基酸，属于中性进化，其中中性替换速率⼀般⽤每个位点每年的变异数来衡量。
利⽤贝叶斯统计或其他⽅法估计物种分歧时间的程序很多，如R8S、MCMCTREE、MULTIDIVTIME、BEAST等，通过不同的策略将化⽯时间信息整合到⼀个系统发育树中，从⽽计算得到Divergence time Tree。
下⾯我们对构建Divergence time Tree的⽅法和步骤做⼀个简述：
利⽤单拷贝基因的CDS或PEP序列构建系统发育树，推荐PhyML或Mrbayes 。因为ML树和贝叶斯进化树对核苷酸 / 氨基酸替代模型的选择⾮常敏感，故在进⾏进化树或分化时间构建之前，需对核苷酸 / 氨基酸替代模型进⾏选择
构建获得 newick格式的进化树，如下例⼦：
02获得fossil calibration time
那如何查两两物种间的化⽯分化时间？
在已经构建好的树中添加化⽯时间，主要形式是 '>0.22<6.75'，'>'后接分化时间的最⼩值，'<'后接最⼤值，如下实例：
当样本⽐较多时，要考虑不同距离的节点，只有⼀个节点估算的可能不准确，可以提供三到五个的节点（物种）的分化时间。03分歧时间估计
这⾥主要介绍3种程序：
准备配置⽂件l：
运⾏程序：
-> 先利⽤PAML baseml估计分⽀的核苷酸频率、转换/颠换率等参数
再利⽤paml2modelinf将baseml输出⽂件转换为estbranches输⼊⽂件
利⽤estbranche对树分⽀长度进⾏最⼤似然率估计
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议