基因序列中的密碼子隨著時(shí)間的推移以幾乎恒定的比例進(jìn)行替換迹卢，具有恒定的演化速率裸影，因此兩個(gè)物種時(shí)間的遺傳距離與物種分歧時(shí)間成正比。

本文就利用mcmctree估算物種分歧時(shí)間做一個(gè)介紹。mcmctree 屬于PAML軟件包柠掂，之前我們介紹過(guò)PAML中的codeml. paml計(jì)算 KaKs值。

關(guān)于mcmctree說(shuō)明文檔可查看http://abacus.gene.ucl.ac.uk/software/MCMCtree.Tutorials.pdf

關(guān)于PAML軟件的安裝可查看之前codml的筆記即可依沮。

本文以軟件的實(shí)例數(shù)據(jù)為例子進(jìn)行介紹

1 示例文件練習(xí)

示例文件位于paml4.9j/examples/DatingSoftBound 下涯贞；

輸入文件主要包括3個(gè)文件：

• 序列比對(duì)文件（nucl或者pep都可以）
• 帶有化石校準(zhǔn)的系統(tǒng)發(fā)育樹
• mcmctree.ctl枪狂； 控制文件

（1）序列比對(duì)文件

mtCDNApri123.txt 內(nèi)容如下所示，是一個(gè)phlip格式

該示例文件分為3個(gè)區(qū)塊宋渔，密碼子在3個(gè)位點(diǎn)的多序列比對(duì)結(jié)果州疾。

此外，也可以輸入密碼子或氨基酸序列比對(duì)信息皇拣，則需要再配置文件中修改相應(yīng)的參數(shù)來(lái)表明數(shù)據(jù)類型 严蓖。若使用氨基酸序列來(lái)進(jìn)行分析，由于mcmctree命令不能選擇較好的氨基酸替換模型進(jìn)行分析氧急，需要自己手動(dòng)運(yùn)行codeml進(jìn)行分析后颗胡，在生成中間文件用于運(yùn)行mcmctree，比較麻煩吩坝。因此杭措，推薦按上述方法使用密碼子三位點(diǎn)的堿基序列作為輸入信息進(jìn)行PAML分析。

(2) 帶有校準(zhǔn)點(diǎn)的有根樹

示例文件如下

 7  1

((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);

第一行钾恢，7個(gè)物種手素，1顆樹，中間用空格分開瘩蚪。
第二行泉懦，其中包含有校準(zhǔn)點(diǎn)信息。校準(zhǔn)點(diǎn)信息一般指95%HPD（Highest Posterior Density）對(duì)應(yīng)的置信區(qū)間疹瘦，單位為崩哩。尾部一定要有分號(hào)。

（3）配置文件mcmctree.ctl

seed = -1  # 表示使用系統(tǒng)當(dāng)前時(shí)間為隨機(jī)數(shù)
       seqfile = mtCDNApri123.txt # 多序列比對(duì)文件
      treefile = mtCDNApri.trees # 帶有校準(zhǔn)信息有根樹
      mcmcfile = mcmc.txt #輸出mcmc信息文本文件言沐，可以用Tracer軟件查看
       outfile = out.txt # 輸出文件邓嘹，記錄了超度量樹和進(jìn)化速率等參數(shù)信息

         ndata = 3 # 輸入的多序列比對(duì)的數(shù)據(jù)個(gè)數(shù)，這里密碼子3個(gè)位置的數(shù)據(jù)险胰；如果有一個(gè)汹押，則設(shè)置為1
       seqtype = 0    * 0: nucleotides; 1:codons; 2:AAs #數(shù)據(jù)類型
       usedata = 1    * 0: no data; 1:seq like; 2:normal approximation; 3:out.BV (in.BV) # 設(shè)置是否利用多序列比對(duì)的數(shù)據(jù)：\
#0，表示不使用多序列比對(duì)數(shù)據(jù)起便，則不會(huì)進(jìn)行l(wèi)ikelihood計(jì)算棚贾，雖然能得到mcmc樹且計(jì)算速度飛快，但是其分歧時(shí)間結(jié)果是有問(wèn)題的榆综；\
#1妙痹，表示使用多序列比對(duì)數(shù)據(jù)進(jìn)行l(wèi)ikelihood計(jì)算，正常進(jìn)行MCMC鼻疮，是一般使用的參數(shù); \
#2怯伊，進(jìn)行正常的approximation likelihood分析，此時(shí)不需要讀取多序列比對(duì)數(shù)據(jù)判沟，直接讀取當(dāng)前目錄中的in.BV文件耿芹。該文件是使用usedata = 3參數(shù)生成的out.BV文件重命名而來(lái)的崭篡。\
#此外，由于程序BUG猩系，當(dāng)設(shè)置usedata = 2時(shí)，一定要在改行參數(shù)后加 *中燥，否則程序報(bào)錯(cuò) Error: file name empty.. \
#3寇甸，程序利用多序列比對(duì)數(shù)據(jù)調(diào)用baseml/codeml命令對(duì)數(shù)據(jù)進(jìn)行分析，生成out.BV文件疗涉。由于mcmctree調(diào)用baseml/codeml進(jìn)行計(jì)算的參數(shù)設(shè)置可能不太好（特別時(shí)對(duì)蛋白序列進(jìn)行計(jì)算時(shí)）拿霉，\
#推薦自己修該軟件自動(dòng)生成的baseml/codeml配置文件，然后再手動(dòng)運(yùn)行baseml/codeml命令咱扣，再整合其結(jié)果文件為out.BV文件绽淘。
         clock = 2    * 1: global clock; 2: independent rates; 3: correlated rates # (1) 全局分子鐘,適用于近緣物種（兩物種分化時(shí)間小于20個(gè)million分化時(shí)間或者序列差異小于5%即近緣物種）\
# (2) 樹枝上進(jìn)化速率服從獨(dú)立同分布的對(duì)數(shù)正態(tài)分布（推薦使用）
# (3) 樹枝上的進(jìn)化速率自相關(guān)。
       RootAge = '<1.0'  * safe constraint on root age, used if no fossil for root. # 設(shè)置root節(jié)點(diǎn)的分歧時(shí)間闹伪，一般設(shè)置一個(gè)最大值沪铭。

         model = 0    * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
# 堿基替換模型：
  #（0）無(wú)參數(shù)（適用于近緣物種）
  #（1）參數(shù)為轉(zhuǎn)換顛換比Kappa
  #（2）參數(shù)為T，C偏瓤，A杀怠，G的頻率
 #（3）最復(fù)雜的進(jìn)化模型
 #（4）參數(shù)為T，C厅克，A赔退，G的頻率+kappa（適用于遠(yuǎn)緣物種）
         alpha = 0.5    * alpha for gamma rates at sites 
# 核酸序列中不同位點(diǎn)，其進(jìn)化速率不一致证舟，其變異速率服從GAMMA分布硕旗。一般設(shè)置GAMMA分布的alpha值為0.5。\
# 若該參數(shù)值設(shè)置為0女责，則表示所有位點(diǎn)的進(jìn)化速率一致漆枚。此外，當(dāng)userdata = 2時(shí)抵知，alpha浪读、ncatG、alpha_gamma辛藻、kappa_gamma這4個(gè)GAMMA參數(shù)無(wú)效碘橘。因?yàn)閡serdata = 2時(shí)，不會(huì)利用到多序列比對(duì)的數(shù)據(jù)吱肌。
         ncatG = 5    * No. categories in discrete gamma # 設(shè)置離散型GAMMA分布的categories值痘拆。

     cleandata = 0    * remove sites with ambiguity data (1:yes, 0:no)?

       BDparas = 1 1 0.1  * birth, death, sampling # 設(shè)置出生率、死亡率和取樣比例氮墨。若輸入有根樹文件中的時(shí)間單位發(fā)生改變纺蛆，則需要相應(yīng)修改出生率和死亡率的值吐葵。
例如，時(shí)間單位由100Myr變換為1Myr桥氏，則要設(shè)置成".01 .01 0.1"温峭。
   kappa_gamma = 6 2      * gamma prior for kappa  #設(shè)置kappa（轉(zhuǎn)換/顛換比率）的GAMMA分布參數(shù)
   alpha_gamma = 1 1      * gamma prior for alpha #設(shè)置GAMMA形狀參數(shù)alpha的GAMMA分布參數(shù)

   rgene_gamma = 2 20 1   * gammaDir prior for rate for genes #設(shè)置序列中所所有位點(diǎn)平均[堿基/密碼子/氨基酸]替換率的Dirichlet-GAMMA分布參數(shù)：
  sigma2_gamma = 1 10 1   * gammaDir prior for sigma^2     (for clock=2 or 3) #設(shè)置所有位點(diǎn)進(jìn)化速率取對(duì)數(shù)后方差（sigma的平方）的Dirichlet-GAMMA分布參數(shù)：alpha=1、beta=10字支、初始方差值為1凤藏。
# 當(dāng)clock參數(shù)值為1時(shí)，表示使用全局的進(jìn)化速率堕伪，各分枝的進(jìn)化速率沒有差異揖庄，即方差為0，該參數(shù)無(wú)效欠雌；
# 當(dāng)clock參數(shù)值為2時(shí)蹄梢，若修改了時(shí)間單位，該參數(shù)值不需要改變富俄；
# 當(dāng)clock參數(shù)值為3時(shí)禁炒，若修改了時(shí)間單位，該參數(shù)值需要改變霍比。

      finetune = 1: .1 .1 .1 .1 .1 .1 * auto (0 or 1): times, musigma2, rates, mixing, paras, FossilErr #冒號(hào)前的值設(shè)置是否自動(dòng)進(jìn)行finetune齐苛，一般設(shè)置成1，然后程序自動(dòng)進(jìn)行優(yōu)化分析

         print = 1   * 0: no mcmc sample; 1: everything except branch rates 2: everything # 設(shè)置打印mcmc的取樣信息
#0桂塞，不打印mcmc結(jié)果凹蜂；
#1，打印除了分支進(jìn)化速率的其它信息（各內(nèi)部節(jié)點(diǎn)分歧時(shí)間阁危、平均進(jìn)化速率玛痊、sigma2值）；
#2狂打，打印所有信息擂煞。 
        burnin = 2000 #將前2000次迭代burnin后，再進(jìn)行取樣（即打印出該次迭代計(jì)算的結(jié)果信息趴乡，各內(nèi)部節(jié)點(diǎn)分歧時(shí)間对省、平均進(jìn)化速率、sigma2值和各分支進(jìn)化速率等）
      sampfreq = 10 #每10次迭代則取樣一次晾捏。
       nsample = 20000 # 當(dāng)取樣次數(shù)達(dá)到該次數(shù)時(shí)蒿涎，則取樣結(jié)束（程序也將運(yùn)行結(jié)束）。

運(yùn)行

~/paml4.9j/src/mcmctree mcmctree.ctl

序在運(yùn)行過(guò)程當(dāng)中惦辛，屏幕也會(huì)生成一些信息劳秋。

查看每列中的值（接受比例、時(shí)間和速率）很重要乱灵。它們應(yīng)該在MCMC運(yùn)行期間保持穩(wěn)定塑崖。如果驗(yàn)收比例變化太大（特別是在開始時(shí)），這意味著預(yù)燒時(shí)間不夠長(zhǎng)痛倚，因此應(yīng)增加控制文件中的預(yù)燒變量

在最后规婆，給出各個(gè)內(nèi)部節(jié)點(diǎn)的分歧時(shí)間(t)、平均變異速率(mu)蝉稳、變異速率方差(sigma2)和r_left的Posterior信息：
均值(mean)抒蚜、95%雙側(cè)置信區(qū)間(95% Equal-tail CI)和95% HPD置信區(qū)間(95% HPD CI)等信息。
此外耘戚，第二列給出了各個(gè)內(nèi)部節(jié)點(diǎn)的Posterior mean time信息嗡髓，可以用于。jnode是程序multipvtime使用的節(jié)點(diǎn)號(hào)收津，由JeffThorne編寫饿这。
倒數(shù)第二列值： 比如第8節(jié)點(diǎn)； 0.0436=0.1850-0.1414

結(jié)果分析

運(yùn)行會(huì)得到以下4個(gè)文件

• FigTree.tre 生成含有分歧時(shí)間的超度量樹文件
  適用于Andrew Rambaut的figtree程序（http://tree.bio.ed.ac.uk/software/figtree/）

節(jié)點(diǎn)處的位分化時(shí)間撞秋，比如 chimpanzee 和bonobo分化時(shí)間為3.02 MYA长捧。該數(shù)值與out .tee文件最后面的數(shù)值對(duì)應(yīng)，也可查看95%置信區(qū)間

• mcmc.txt MCMC取樣信息吻贿，包含各內(nèi)部節(jié)點(diǎn)分歧時(shí)間唆姐、平均進(jìn)化速率、sigma2值等信息，可以在Tracer軟件中打開奉芦。通過(guò)查看各參數(shù)的ESS值赵抢，若ESS值大于200，則從一定程度上表示MCMC過(guò)程能收斂声功，結(jié)果可靠烦却。

在我們的示例中，該文件有50列(如果print=2)或14列(如果print=1)和20002行先巴。第一列是樣本的生成(迭代)數(shù)其爵，接下來(lái)的48列(print=2)或12列(print=1)分別對(duì)應(yīng)于所分析的48個(gè)或12個(gè)參數(shù)；最后一列是似然值(第50列伸蚯，如果print=2摩渺；或第14列，如果print=1)剂邮。行數(shù)與MCMC期間采集的樣本數(shù)相對(duì)應(yīng)摇幻。“mcmc.txt”文件適用于Tracer程序分析（http://beast.bio.ed.ac.uk/tracer）

• out.txt 包含由較多信息的結(jié)果文件挥萌。例如绰姻，各堿基頻率、節(jié)點(diǎn)命名信息引瀑、各節(jié)點(diǎn)分歧時(shí)間狂芋、進(jìn)化速率和進(jìn)化樹等

打開此文件。文件開頭有很多垃圾憨栽，你通持姆可以忽略。向下滾動(dòng)該文件屑柔，直到您看到3個(gè)系統(tǒng)發(fā)育樹(如果print=1屡萤；如果print=2，您將看到6個(gè))
第一棵樹只包含節(jié)點(diǎn)標(biāo)簽锯蛀；第二棵樹包含以時(shí)間單位表示的分支長(zhǎng)度灭衷；第三棵樹包含以時(shí)間單位表示的分支長(zhǎng)度，以及節(jié)點(diǎn)年齡的可信度區(qū)間（記住只有節(jié)點(diǎn)有）旁涤。最后三棵樹（僅當(dāng)您將“print”設(shè)置為2時(shí)寫入輸出文件）具有替換率翔曲，而不是分支長(zhǎng)度，每個(gè)樹表示每個(gè)位點(diǎn)（在我們的示例中劈愚，三個(gè)密碼子位置中的每一個(gè)）瞳遍。在樹下面，是后驗(yàn)均值菌羽，與屏幕輸出最后的一致掠械。

• SeedUsed 文件
  文件包含用于初始化MCMC的隨機(jī)種子。如果將文件中的值（在我的例子中是1847523971）復(fù)制到控制文件的seed變量中，則可以再次運(yùn)行分析并獲得相同的結(jié)果

注意事項(xiàng)

（1） 如何檢測(cè)MCMC的結(jié)果是否達(dá)到收斂狀態(tài)猾蒂？
收斂的意思均唉，即經(jīng)過(guò)很多次迭代后，得到的MCMC樹的各枝長(zhǎng)趨于一個(gè)定值肚菠，變動(dòng)非常小舔箭。
于是，最直接的檢測(cè)方法是：分別使用不同的Seed值進(jìn)行mcmctree ()层扶，進(jìn)行兩次或多次分析，然后比較兩個(gè)結(jié)果樹是否一致烙荷，實(shí)際就是比較樹文件中各內(nèi)部節(jié)點(diǎn)的Height值（分歧時(shí)間 / Posterior time）镜会。(可查看out.txt文件后面分化時(shí)間（第二列）)
可通過(guò)以下3種方法進(jìn)行：

• 結(jié)果文件mcmc.txt中每行記錄了一個(gè)取樣的結(jié)果，包含各內(nèi)部節(jié)點(diǎn)的Posterior time值终抽，即進(jìn)化樹的Height值戳表，對(duì)該mcmc.txt文件進(jìn)行分析，得到每個(gè)內(nèi)部節(jié)點(diǎn)的Posterior mean time值拿诸。然后作圖 （官方說(shuō)明文檔第6頁(yè)圖示）扒袖，橫坐標(biāo)表示第一次運(yùn)行的各內(nèi)部節(jié)點(diǎn)的分歧時(shí)間均值塞茅、縱坐標(biāo)表示第二次運(yùn)行的各內(nèi)部節(jié)點(diǎn)的分歧時(shí)間均值亩码。該散點(diǎn)圖要表現(xiàn)出非常明顯的對(duì)角線，才認(rèn)為達(dá)到收斂野瘦。這時(shí)可以考慮計(jì)算其相關(guān)系數(shù)來(lái)判斷是否符合線性描沟。
• 直接比較兩次結(jié)果樹文件中的各枝長(zhǎng)。計(jì)算各枝長(zhǎng)總的偏差百分比鞭光，當(dāng)偏差百分比低于0.1%吏廉，則認(rèn)為兩次結(jié)果非常吻合，差異低于0.1%惰许，認(rèn)為達(dá)到收斂席覆。(如上述所述)
• 此外，還可以使用Tracer分析mcmc.txt文件汹买，檢測(cè)其ESS值佩伤，一般認(rèn)為該值高于200，則可能達(dá)到收斂晦毙。該方法可用于輔助檢測(cè)生巡。最后，若不收斂见妒，則需要提高burnin孤荣、nsample值，重新運(yùn)行程序。

以下圖片是根據(jù)兩次運(yùn)行結(jié)果的分化時(shí)間進(jìn)行作圖：

（2）如何設(shè)置burnin盐股、sampfreq和nsample值钱豁？

• 一般推薦設(shè)置nsample值不小于20k(官方示例數(shù)據(jù)中設(shè)置的值)，只有當(dāng)該值較大時(shí)疯汁，求得的均值才有意義寥院。
• sampfreq表示取樣間隔，一般設(shè)置為10涛目，100或1000秸谢。nsample 和sampfreq 的乘積表示有效迭代的次數(shù)，這些迭代越準(zhǔn)確霹肝，次數(shù)越大估蹄，則結(jié)果越好，越趨于收斂沫换。當(dāng)然臭蚁，次數(shù)越大，越消耗時(shí)間讯赏。若數(shù)據(jù)量很小垮兑，可以考慮設(shè)置sampfreq為1000；若數(shù)據(jù)量大漱挎，每次迭代耗時(shí)多系枪，則考慮設(shè)置sampfre為10
• 一般最開始的迭代結(jié)果很差，需要去除(burnin)其結(jié)果磕谅，則設(shè)置burnin值私爷。要設(shè)置足夠大的burnin值，保證在程序運(yùn)行時(shí)當(dāng)?shù)壤_(dá)到0% (即burnin結(jié)束) 時(shí)膊夹，其參數(shù)的接受比例值較好衬浑，在30%左右且隨迭代次數(shù)增多時(shí)變動(dòng)幅度不大。推薦設(shè)置burnin的迭代次數(shù)為有效次數(shù)的10~40%放刨。PAML軟件時(shí)先burnin后再記錄有效迭代的參數(shù)值工秩；其它軟件如BEAST則記錄所有的參數(shù)值后，最后匯總時(shí)burnin掉一定比例的數(shù)據(jù)进统。
• 總體上助币，其實(shí)就是兩個(gè)參數(shù)：burnin掉的迭代次數(shù)和用于計(jì)算結(jié)果的有效迭代次數(shù)。由于需要根據(jù)有效迭代數(shù)據(jù)來(lái)進(jìn)行均值計(jì)算麻昼，若記錄每次迭代的參數(shù)奠支，則生成的mcmc.txt文件行數(shù)過(guò)多，文件太大抚芦，匯總時(shí)也消耗計(jì)算倍谜。于是設(shè)置沒隔一定迭代輪數(shù)取樣一次迈螟，這樣生成的mcmc.txt文件不會(huì)太大，對(duì)最后的均值計(jì)算影響不大尔崔。

我個(gè)人推薦有效迭代次數(shù) 1~10M 次答毫，burnin 0.2~4M次。按時(shí)間消耗從少到多的參數(shù)值：
數(shù)據(jù)簡(jiǎn)單時(shí)季春，進(jìn)行0.5M迭代次數(shù)洗搂，burnin比例40%：{ burnin 200k; sampfreq 10; nsample 50k }
數(shù)據(jù)中等時(shí)，進(jìn)行1M迭代次數(shù)载弄，burnin比例40%：{ burnin 400k; sampfreq 10; nsample 100k }
數(shù)據(jù)復(fù)雜時(shí)耘拇，進(jìn)行5M迭代次數(shù)，burnin比例20%：{ burnin 1000k; sampfreq 10; nsample 500k }
數(shù)據(jù)巨大時(shí)宇攻，進(jìn)行10M迭代次數(shù)惫叛，burnin比例20%：{ burnin 2000k; sampfreq 100; nsample 100k }

使用相應(yīng)參數(shù)進(jìn)行分析后，若不滿意其收斂程度逞刷，則選用更多迭代次數(shù)的參數(shù)嘉涌。

• 帶有化石校準(zhǔn)的系統(tǒng)發(fā)育樹

利用物種間單拷貝基因進(jìn)行構(gòu)建系統(tǒng)發(fā)育樹，如果要選擇核苷酸/氨基酸替換模型夸浅，則可以使用IQ-TREE,而后使用RAxML進(jìn)行構(gòu)建進(jìn)化樹仑最。

• 兩物種間的化石分化時(shí)間
  通過(guò)time tree網(wǎng)站,可進(jìn)行物種化石分化時(shí)間查詢，比如
  查詢帆喇，人和黑猩猩
  
  
  查詢以后會(huì)得到下面結(jié)果
  

可以看到其為： 5.1-11.8 MYA, 則將其添加到進(jìn)化樹中就是>.051<0.118警医。因?yàn)槠鋯挝粸?00百萬(wàn)年

2 使用approximate likelihood方法更快速地計(jì)算分歧時(shí)間

按照以上方法進(jìn)行MCMC計(jì)算非常消耗計(jì)算休讳。為了加快計(jì)算速度，可以將分析分成兩個(gè)步驟：
（1）首先源梭，使用Maximum Likelihood方法計(jì)算枝長(zhǎng)和Hessian信息氓侧；
（2）再使用MCMC方法計(jì)算分歧時(shí)間。

具體步驟：
第一步：返回上一次運(yùn)行的文件夾贤重，構(gòu)建Hessian的文件夾,拷貝樹，比對(duì)文件，配置文件

mkdir Hessian
cp ../mcmctree.ctl ../mtCDNApri123.txt ../mtCDNApri.trees ./

打開配置文件路星，將usedate 設(shè)置為3，并運(yùn)行

usedata = 3    * 0: no data; 1:seq like; 2:normal approximation; 3:out.BV (in.BV)
/soft/paml4.9j/src/mcmctree mcmctree.ctl

第二步：新建文件夾诱桂，拷貝樹洋丐，配置文件，比對(duì)文件以及out.BV到新文件夾
上一步得到out.BV中含有各個(gè)枝的長(zhǎng)度挥等，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 根據(jù)上一步的gradient 和Hessian，進(jìn)行估計(jì)分化時(shí)間肝劲。此外迁客，還可以再次建立新文件夾approx02, 再次運(yùn)行第二步郭宝，兩次結(jié)果進(jìn)行比較。

??掷漱； 如果clock=1粘室，則不能用approximate 方法。

參考

• MCMCtree——最詳盡的物種進(jìn)化歷史分歧時(shí)間估算教程
• 使用PAML進(jìn)行分歧時(shí)間計(jì)算
• 進(jìn)化樹專題（九）| 從DNA序列到分化時(shí)間——如何構(gòu)建帶分歧時(shí)間的進(jìn)化樹卜范？