以下內容參考了博客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，含有對程序的指令。

控制文件示例

各類型參數(shù)的含義如下：

1）輸入輸出參數(shù)

seed=-1 ：

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

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

如果需要一個可重復的結果，可以設置一個特定的奇數(shù)或偶數(shù)

seqfile：比對序列文件名

treefile：樹文件文件名

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

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

2）數(shù)據使用說明參數(shù)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

cleandata = 0不需要

cleandata = 1 需要

3）位點替換模型參數(shù)

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，若該參數(shù)設置為0，則表示所有位點的進化速率一致。

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

BDparas= 1 1 0.1： 控制出生-死亡過程的參數(shù)，設置出生率、死亡率和取樣比例。

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

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

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

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

3）進化速率參數(shù)

rgene_gamma = 2 20 1：設置序列中所所有位點平均（堿基/密碼子/氨基酸）替換率的Dirichlet-GAMMA分布參數(shù)：

alpha=2、beta=20、初始平均替換率為每100million年（取決于輸入有根樹文件中的時間單位）1個替換。若時間單位由100Myr變換為1Myr，則要設置成"2 2000 1"?？傮w上的平均進化速率為：2 / 20 = 0.1 個替換 / 每100Myr，即每個位點每年的替換數(shù)為 1e-9。

sigma2_gamma = 1 10 1：設置所有位點進化速率取對數(shù)后方差（sigma的平方）的Dirichlet-GAMMA分布參數(shù)：alpha=1、beta=10、初始方差值為1。

當clock參數(shù)值為1時，表示使用全局的進化速率，各分枝的進化速率沒有差異，即方差為0，該參數(shù)無效；

當clock參數(shù)值為2時，若修改了時間單位，該參數(shù)值不需要改變；

當clock參數(shù)值為3時，若修改了時間單位，該參數(shù)值需要改變。

finetune = 1: .1 .1 .1 .1 .1 .1：冒號前的值設置是否自動進行finetune，一般設置成1，然后程序自動進行優(yōu)化分析；冒號后面設置各個參數(shù)的步進值：times, musigma2, rates, mixing, paras, FossilErr。由于有了自動設置，該參數(shù)不像以前版本那么重要了，可能以后會取消該參數(shù)。

4）mcmc參數(shù)

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

print = 0：不打印mcmc結果

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

print = 2：打印所有信息

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

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

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

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

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

數(shù)據簡單時，進行0.5M迭代次數(shù)，burnin比例40%：{ burnin 200k; sampfreq 10; nsample 50k }

數(shù)據中等時，進行1M迭代次數(shù)，burnin比例40%：{ burnin 400k; sampfreq 10; nsample 100k }

數(shù)據復雜時，進行5M迭代次數(shù)，burnin比例20%：{ burnin 1000k; sampfreq 10; nsample 500k }

數(shù)據巨大時，進行10M迭代次數(shù)，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進行物種分化時間的估算

在對較大的數(shù)據進行l(wèi)ikelihood function計算時，計算成本昂貴，耗時長，下面介紹針對較大比對數(shù)據的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。