以下內容參考了博客https://www.cnblogs.com/bio-mary/p/12818888.html竹勉、前研究組師妹的總結資料和MCMCTree tutorials名船。

"MCMCTree performs Bayesian estimation of species divergence times using soft fossil constraints under various molecular clock models."

MCMCTree是PAML包的一部分，功能是在多種分子鐘模型下栈雳，利用較寬松的的化石約束薪贫，用貝葉斯方法估算物種分化時間。

安裝MCMCTree

在Linux系統(tǒng)下的安裝很簡單，用conda安裝paml就可以：

conda install -c bioconda paml

準備文件

"The program uses for input a sequence alighment (nucleotide or protein), a phylogenetic tree with fossil calibrations, and a control file (ususally called mcmctree.ctl) that contains the instructions for the program. "

"MCMCTree is part of the PAML package."

運行mcmctree需要準備三個文件：比對序列文件趣避、樹文件和控制文件。

在paml軟件包下有個example文件夾扶欣，里面有很多示例文件鹅巍，用于mcmctree分析的示例文件在DatingSoftBound里面千扶。序列文件名稱：mtCDNApri123，樹文件名稱：mtCDNApri.trees骆捧，控制文件的名稱mcmctree.ctl澎羞，該示例文件分析的是7個類人猿物種的線粒體蛋白編碼基因。

1 比對序列文件（核苷酸或蛋白質）

我用的比對好的序列文件是.phy文件敛苇，在軟件Geneious里完成的文件格式轉化妆绞，先導入fasat文件，導出時彈出的對話框phylip alignment export中枫攀，選擇Relaxed phylip-Full length names followed by a single space括饶，意思是導出的文件中每條序列的序列名和序列之間有一個空格。但是mcmctree要求序列名和序列之間至少兩個空格空格来涨，如果序列少的化图焰，可以手動增加空格。我的序列比較多蹦掐，用的方法是：使用geneious里batch rename功能統(tǒng)一給序列名后面加個序列名里沒有的字符技羔，然后導出的phylip文件用編輯器打開，用空格替換那個字符卧抗。

下圖是示例文件的展示藤滥，格式是txt。7表示的是物種樹社裆，3331是比對序列的長度拙绊。在示例文件中，比對序列被分成了3部分泳秀，對應第一标沪、第二和第三密碼子位點， 每一部分都如下所示晶默。

序列文件示例

2 有化石校準信息的樹文件

需要注意的是樹文件必須是定根的谨娜、不帶有枝長信息的二歧樹（rooted binary tree）

paml軟件包下面的example文件夾下的樹文件（.tree）如下：

樹文件示例

7代表的是物種樹，1代表只有一顆樹磺陡。藍色框中顯示的是化石時間標記趴梢，表示的是人（humna）和黑猩猩（chimpanzee，bonobo）的共同祖先出現(xiàn)的時間在0.06到0.08個100個百萬年之間币他。

注：后來我一個師妹跟我說這種化石時間標記有程序bug坞靶，前面的>.06不能被程序讀取，推薦了另一種化石時間標記方式：'B(.06,.08)'

（注意蝴悉，這里化石時間標記的單位是100個百萬年彰阴，而不是1個百萬年）

樹文件代表的樹的拓撲結構如下：

示例樹文件的拓撲結構

我們自己的樹文件常常是帶枝長信息，想要去除枝長信息的話可以用notepad等編輯器打開樹文件拍冠，手動刪除枝長信息尿这，還有一種簡單的方法簇抵，在linux系統(tǒng)中運行下面的命令行，tree.nwk是輸入的帶枝長信息的樹文件射众，output_tree.nwk是輸出的不帶枝長信息的樹文件

perl -ne '$_=~s/:[\d\.]+//g; print $_;' tree.nwk > output_tree.nwk

3 控制文件

控制文件常命名為mcmctree.ctl碟摆，含有對程序的指令。

控制文件示例

各類型參數的含義如下：

1）輸入輸出參數

seed=-1 ：

#設置一個隨機數（正整數或負整數）來運行程序叨橱，若設置為-1典蜕，表示使用系統(tǒng)當前時間為隨機數，因此每次運行mcmctree得到的結果會有所不同罗洗。

建議用seed=-1運行程序至少兩次愉舔，檢驗兩次運行的結果是否相近。如果每次運行的結果差別很大伙菜，可能是多種原因導致的：緩慢收斂（slow convergence）轩缤、不良混合（poor mixing）、采樣不完全（insufficient samples taken）, 或程序錯誤贩绕。

如果需要一個可重復的結果典奉，可以設置一個特定的奇數或偶數

seqfile：比對序列文件名

treefile：樹文件文件名

mcmcfile：針對設定的參數運行的MCMC的報告文件，可以用Tracer軟件查看丧叽。

outfile：一旦程序運行完成，總結性的結果文件會寫入該文件公你，該文件記錄了超度量樹和進化速率等參數信息踊淳。

2）數據使用說明參數

ndata：比對序列的分區(qū)數，在本例中陕靠，根據核苷酸在密碼子中的位置（第一位迂尝、第二位和第三位）分成了三個分區(qū)，因此ndata=3剪芥。

自己的核苷酸比對數據怎么根據密碼子第1垄开、2、3位的位置進行分區(qū)呢税肪？推薦一個在線的工具Split codons：https://www.bioinformatics.org/sms2/split_codons.html

usedata：設置是否利用多序列比對的數據

usedata = 0：不適用多序列比對數據溉躲，不會進行l(wèi)ikelihood計算，雖然能得到mcmc樹且計算速度飛快益兄，但分歧時間結果有問題

usedata = 1 使用多序列數據進行l(wèi)ikelihood計算锻梳，正常進行MCMC，是一般的使用參數

usedata = 2 進行正常的approximation likelihood分析净捅，此時不需要讀取多序列比對數據疑枯，直接讀取當前目錄中的in.BV文件，該文件是使用usedata=3參數生成的out.BV文件重命名而來蛔六。

usedata = 3：程序利用多序列比對數據調用baseml/codeml命令對數據進行分析荆永，生成out.BV文件废亭。

由于mcmctree調用baseml/codeml進行計算的參數設置可能不太好，尤其是蛋白序列具钥，推薦自己修改該軟件自動生成的baseml/codeml配置文件豆村，然后手動運行baseml/codeml命令，再整合其結果文件為out.BV文件氓拼。

seqtype：比對序列的類型你画，=0代表核苷酸序列，=1代表密碼子序列桃漾，=2代表氨基酸序列

clock：選擇使用的分子鐘模型坏匪，

clock=1：global clock的方法，假設所有分支的進化速率一致

clock=2：使用獨立速率模型（independent rates model）撬统，在該模型中适滓，速率符合對數-正態(tài)分布（也就是說速率的對數符合正態(tài)分布）

clock=3： 使用correlated rates的方法，和方法類似恋追，但是log（r）的方差和時間（t）相關

RootAge：如果我們在樹文件中沒有對根提供時間校準凭迹，那么用參數對根提供一個時間校準。在示例文件中使用了<1.0苦囱，意思是所有猿的最近共同祖先出現(xiàn)的時間不會大于100個百萬年前嗅绸。

cleandata:設置是否移除不明確的字符或gap后再進行數據分析。

cleandata = 0不需要

cleandata = 1 需要

3）位點替換模型參數

model: 設置要使用的替換模型

model = 0: JC69,計算很快

model = 1: K80,

model = 2: F81,

model = 3: F84,

model = 4: HKY85

model = 7: GTR

alpha：核酸序列中不同位點撕彤，其進化速率不一致鱼鸠，其變異速率服從Gamma分布，一般設置gamma分布的alpha值為0.5：alpha = 0.5羹铅，若該參數設置為0蚀狰，則表示所有位點的進化速率一致。

ncatG = 5 ：設置離散型GAMMA分布的categories值

BDparas= 1 1 0.1： 控制出生-死亡過程的參數职员，設置出生率麻蹋、死亡率和取樣比例。

若輸入有根樹文件中的時間單位發(fā)生變化焊切，則需相應修改出生率和死亡率的值扮授。例如時間單位由100Myr變換成1myr，則要設置成：.01.01 .01

kappa_gamma = 6 2 ：設置替換模型參數kappa（轉換/顛換比率）的GAMMA先驗专肪。

alpha_gamma = 1 1： 設置替換模型中gamma形狀參數（用于反應位點上不同的速率）alpha的GAMMA先驗糙箍。

注意：當usedata = 2時，alpha牵祟、ncatG深夯、alpha_gamma、kappa_gamma四個GAMMA參數無效，因為此時不會用到多序列比對的數據

3）進化速率參數

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，然后程序自動進行優(yōu)化分析霎冯；冒號后面設置各個參數的步進值：times, musigma2, rates, mixing, paras, FossilErr。由于有了自動設置钞瀑，該參數不像以前版本那么重要了沈撞，可能以后會取消該參數。

4）mcmc參數

print：設置打印mcmc的取樣信息

print = 0：不打印mcmc結果

print = 1：打印除了分支進化速率的其他信息（各內部節(jié)點分歧時間雕什、平均進化速率缠俺，sigma2值）

print = 2：打印所有信息

burnin = 2000 ：將前2000次迭代burnin后，再進行取樣（即打印出打印出該次迭代計算的結果信息贷岸，各內部節(jié)點分歧時間壹士、平均進化速率，sigma2值和各分支進化速率等）

samplefreq = 10 ：每10次迭代取樣一次

nsample = 20000 ：當取樣次數達到該次數時偿警，則取樣結束（程序也運行結束）

burnin= 2000躏救，samplefreq= 10，nsample= 20000：總的意思是程序會去除掉（burnin)前2000次迭代的結果，然后每10次迭代取一次樣盒使，直到取樣次數達到20000次崩掘，因此MCMC會運算2000+10*20000次迭代。一般來說少办，程序需要進行10000-20000次取樣苞慢，才能獲得較好的數據總結。一般來說英妓，如果想通過增加MCMC長度來提高收斂效果挽放，一般是通過增加samplefreq，但保持nsample在一個比較合理的范圍蔓纠。

別人的博客曾進行過一個其個人的總結辑畦，照抄如下：

數據簡單時，進行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

在Linux系統(tǒng)下運行時乖杠，將序列文件、樹文件和控制文件放在同一個路徑下澄成，在該路徑下運行程序：

which mcmctree? #查找mcmctree的路徑

~/miniconda3/bin/mcmctree mcmctree.ctl? #~/miniconda3/bin/mcmctree是程序的絕對路徑胧洒，mcmctree.ctl是控制文件的名字

用approximate likelihood calculation進行物種分化時間的估算

在對較大的數據進行l(wèi)ikelihood function計算時，計算成本昂貴墨状，耗時長卫漫，下面介紹針對較大比對數據的approximate likelihood計算的方法。

總共分兩步肾砂，在第一步中列赎，枝長、gradient镐确、Hessian由最大似然法估算出包吝，在第二步中，利用gradient源葫、Hessian進行分化時間的估算诗越，這一步用Taylor expansion的方法創(chuàng)造出近似于最大似然法的功能。

總的操作步驟如下：

新建一個文件夾Hessian息堂，將樹文件嚷狞、序列文件和控制文件拷貝過去，控制文件中的usedata = 3，然后運行程序感耙。程序運行結束后生成文件中有一個out.BV文件

再建立一個新文件夾approx01褂乍，將上一步中的out.BV文件、樹文件即硼、序列文件和控制文件拷貝過去逃片，out.BV文件重命名為in.BV,控制文件中的usedata = 3改為usedata = 2。