以下內(nèi)容參考了博客https://www.cnblogs.com/bio-mary/p/12818888.html、前研究組師妹的總結(jié)資料和MCMCTree tutorials。
"MCMCTree performs Bayesian estimation of species divergence times using soft fossil constraints under various molecular clock models."
MCMCTree是PAML包的一部分，功能是在多種分子鐘模型下，利用較寬松的的化石約束，用貝葉斯方法估算物種分化時(shí)間。
安裝MCMCTree
在Linux系統(tǒng)下的安裝很簡(jiǎn)單，用conda安裝paml就可以：
conda install -c bioconda paml
準(zhǔn)備文件
"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."
運(yùn)行mcmctree需要準(zhǔn)備三個(gè)文件：比對(duì)序列文件、樹(shù)文件和控制文件。
在paml軟件包下有個(gè)example文件夾，里面有很多示例文件，用于mcmctree分析的示例文件在DatingSoftBound里面。序列文件名稱：mtCDNApri123，樹(shù)文件名稱：mtCDNApri.trees，控制文件的名稱mcmctree.ctl，該示例文件分析的是7個(gè)類人猿物種的線粒體蛋白編碼基因。
1 比對(duì)序列文件（核苷酸或蛋白質(zhì)）
我用的比對(duì)好的序列文件是.phy文件，在軟件Geneious里完成的文件格式轉(zhuǎn)化，先導(dǎo)入fasat文件，導(dǎo)出時(shí)彈出的對(duì)話框phylip alignment export中，選擇Relaxed phylip-Full length names followed by a single space，意思是導(dǎo)出的文件中每條序列的序列名和序列之間有一個(gè)空格。但是mcmctree要求序列名和序列之間至少兩個(gè)空格空格，如果序列少的化，可以手動(dòng)增加空格。我的序列比較多，用的方法是：使用geneious里batch rename功能統(tǒng)一給序列名后面加個(gè)序列名里沒(méi)有的字符，然后導(dǎo)出的phylip文件用編輯器打開(kāi)，用空格替換那個(gè)字符。
下圖是示例文件的展示，格式是txt。7表示的是物種樹(shù)，3331是比對(duì)序列的長(zhǎng)度。在示例文件中，比對(duì)序列被分成了3部分，對(duì)應(yīng)第一、第二和第三密碼子位點(diǎn)， 每一部分都如下所示。
2 有化石校準(zhǔn)信息的樹(shù)文件
需要注意的是樹(shù)文件必須是定根的、不帶有枝長(zhǎng)信息的二歧樹(shù)（rooted binary tree）
paml軟件包下面的example文件夾下的樹(shù)文件（.tree）如下：
7代表的是物種樹(shù)，1代表只有一顆樹(shù)。藍(lán)色框中顯示的是化石時(shí)間標(biāo)記，表示的是人（humna）和黑猩猩（chimpanzee，bonobo）的共同祖先出現(xiàn)的時(shí)間在0.06到0.08個(gè)100個(gè)百萬(wàn)年之間。
注：后來(lái)我一個(gè)師妹跟我說(shuō)這種化石時(shí)間標(biāo)記有程序bug，前面的>.06不能被程序讀取，推薦了另一種化石時(shí)間標(biāo)記方式：'B(.06,.08)'
（注意，這里化石時(shí)間標(biāo)記的單位是100個(gè)百萬(wàn)年，而不是1個(gè)百萬(wàn)年）
樹(shù)文件代表的樹(shù)的拓?fù)浣Y(jié)構(gòu)如下：
我們自己的樹(shù)文件常常是帶枝長(zhǎng)信息，想要去除枝長(zhǎng)信息的話可以用notepad等編輯器打開(kāi)樹(shù)文件，手動(dòng)刪除枝長(zhǎng)信息，還有一種簡(jiǎn)單的方法，在linux系統(tǒng)中運(yùn)行下面的命令行，tree.nwk是輸入的帶枝長(zhǎng)信息的樹(shù)文件，output_tree.nwk是輸出的不帶枝長(zhǎng)信息的樹(shù)文件
perl -ne '$_=~s/:[\d\.]+//g; print $_;' tree.nwk > output_tree.nwk
3 控制文件
控制文件常命名為mcmctree.ctl，含有對(duì)程序的指令。
各類型參數(shù)的含義如下：
1）輸入輸出參數(shù)
seed=-1 ：
#設(shè)置一個(gè)隨機(jī)數(shù)（正整數(shù)或負(fù)整數(shù)）來(lái)運(yùn)行程序，若設(shè)置為-1，表示使用系統(tǒng)當(dāng)前時(shí)間為隨機(jī)數(shù)，因此每次運(yùn)行mcmctree得到的結(jié)果會(huì)有所不同。
建議用seed=-1運(yùn)行程序至少兩次，檢驗(yàn)兩次運(yùn)行的結(jié)果是否相近。如果每次運(yùn)行的結(jié)果差別很大，可能是多種原因?qū)е碌模壕徛諗浚╯low convergence）、不良混合（poor mixing）、采樣不完全（insufficient samples taken）, 或程序錯(cuò)誤。
如果需要一個(gè)可重復(fù)的結(jié)果，可以設(shè)置一個(gè)特定的奇數(shù)或偶數(shù)
seqfile：比對(duì)序列文件名
treefile：樹(shù)文件文件名
mcmcfile：針對(duì)設(shè)定的參數(shù)運(yùn)行的MCMC的報(bào)告文件，可以用Tracer軟件查看。
outfile：一旦程序運(yùn)行完成，總結(jié)性的結(jié)果文件會(huì)寫(xiě)入該文件，該文件記錄了超度量樹(shù)和進(jìn)化速率等參數(shù)信息。
2）數(shù)據(jù)使用說(shuō)明參數(shù)
ndata：比對(duì)序列的分區(qū)數(shù)，在本例中，根據(jù)核苷酸在密碼子中的位置（第一位、第二位和第三位）分成了三個(gè)分區(qū)，因此ndata=3。
自己的核苷酸比對(duì)數(shù)據(jù)怎么根據(jù)密碼子第1、2、3位的位置進(jìn)行分區(qū)呢？推薦一個(gè)在線的工具Split codons：https://www.bioinformatics.org/sms2/split_codons.html
usedata：設(shè)置是否利用多序列比對(duì)的數(shù)據(jù)
usedata = 0：不適用多序列比對(duì)數(shù)據(jù)，不會(huì)進(jìn)行l(wèi)ikelihood計(jì)算，雖然能得到mcmc樹(shù)且計(jì)算速度飛快，但分歧時(shí)間結(jié)果有問(wèn)題
usedata = 1 使用多序列數(shù)據(jù)進(jìn)行l(wèi)ikelihood計(jì)算，正常進(jìn)行MCMC，是一般的使用參數(shù)
usedata = 2 進(jìn)行正常的approximation likelihood分析，此時(shí)不需要讀取多序列比對(duì)數(shù)據(jù)，直接讀取當(dāng)前目錄中的in.BV文件，該文件是使用usedata=3參數(shù)生成的out.BV文件重命名而來(lái)。
usedata = 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è)置可能不太好，尤其是蛋白序列，推薦自己修改該軟件自動(dòng)生成的baseml/codeml配置文件，然后手動(dòng)運(yùn)行baseml/codeml命令，再整合其結(jié)果文件為out.BV文件。
seqtype：比對(duì)序列的類型，=0代表核苷酸序列，=1代表密碼子序列，=2代表氨基酸序列
clock：選擇使用的分子鐘模型，
clock=1：global clock的方法，假設(shè)所有分支的進(jìn)化速率一致
clock=2：使用獨(dú)立速率模型（independent rates model），在該模型中，速率符合對(duì)數(shù)-正態(tài)分布（也就是說(shuō)速率的對(duì)數(shù)符合正態(tài)分布）
clock=3： 使用correlated rates的方法，和方法類似，但是log（r）的方差和時(shí)間（t）相關(guān)
RootAge：如果我們?cè)跇?shù)文件中沒(méi)有對(duì)根提供時(shí)間校準(zhǔn)，那么用參數(shù)對(duì)根提供一個(gè)時(shí)間校準(zhǔn)。在示例文件中使用了<1.0，意思是所有猿的最近共同祖先出現(xiàn)的時(shí)間不會(huì)大于100個(gè)百萬(wàn)年前。
cleandata:設(shè)置是否移除不明確的字符或gap后再進(jìn)行數(shù)據(jù)分析。
cleandata = 0不需要
cleandata = 1 需要
3）位點(diǎn)替換模型參數(shù)
model: 設(shè)置要使用的替換模型
model = 0: JC69,計(jì)算很快
model = 1: K80,
model = 2: F81,
model = 3: F84,
model = 4: HKY85
model = 7: GTR
alpha：核酸序列中不同位點(diǎn)，其進(jìn)化速率不一致，其變異速率服從Gamma分布，一般設(shè)置gamma分布的alpha值為0.5：alpha = 0.5，若該參數(shù)設(shè)置為0，則表示所有位點(diǎn)的進(jìn)化速率一致。
ncatG = 5 ：設(shè)置離散型GAMMA分布的categories值
BDparas= 1 1 0.1： 控制出生-死亡過(guò)程的參數(shù)，設(shè)置出生率、死亡率和取樣比例。
若輸入有根樹(shù)文件中的時(shí)間單位發(fā)生變化，則需相應(yīng)修改出生率和死亡率的值。例如時(shí)間單位由100Myr變換成1myr，則要設(shè)置成：.01.01 .01
kappa_gamma = 6 2 ：設(shè)置替換模型參數(shù)kappa（轉(zhuǎn)換/顛換比率）的GAMMA先驗(yàn)。
alpha_gamma = 1 1： 設(shè)置替換模型中g(shù)amma形狀參數(shù)（用于反應(yīng)位點(diǎn)上不同的速率）alpha的GAMMA先驗(yàn)。
注意：當(dāng)usedata = 2時(shí)，alpha、ncatG、alpha_gamma、kappa_gamma四個(gè)GAMMA參數(shù)無(wú)效，因?yàn)榇藭r(shí)不會(huì)用到多序列比對(duì)的數(shù)據(jù)
3）進(jìn)化速率參數(shù)
rgene_gamma = 2 20 1：設(shè)置序列中所所有位點(diǎn)平均（堿基/密碼子/氨基酸）替換率的Dirichlet-GAMMA分布參數(shù)：
alpha=2、beta=20、初始平均替換率為每100million年（取決于輸入有根樹(shù)文件中的時(shí)間單位）1個(gè)替換。若時(shí)間單位由100Myr變換為1Myr，則要設(shè)置成"2 2000 1"。總體上的平均進(jìn)化速率為：2 / 20 = 0.1 個(gè)替換 / 每100Myr，即每個(gè)位點(diǎn)每年的替換數(shù)為 1e-9。
sigma2_gamma = 1 10 1：設(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)化速率沒(méi)有差異，即方差為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：冒號(hào)前的值設(shè)置是否自動(dòng)進(jìn)行finetune，一般設(shè)置成1，然后程序自動(dòng)進(jìn)行優(yōu)化分析；冒號(hào)后面設(shè)置各個(gè)參數(shù)的步進(jìn)值：times, musigma2, rates, mixing, paras, FossilErr。由于有了自動(dòng)設(shè)置，該參數(shù)不像以前版本那么重要了，可能以后會(huì)取消該參數(shù)。
4）mcmc參數(shù)
print：設(shè)置打印mcmc的取樣信息
print = 0：不打印mcmc結(jié)果
print = 1：打印除了分支進(jìn)化速率的其他信息（各內(nèi)部節(jié)點(diǎn)分歧時(shí)間、平均進(jìn)化速率，sigma2值）
print = 2：打印所有信息
burnin = 2000 ：將前2000次迭代burnin后，再進(jìn)行取樣（即打印出打印出該次迭代計(jì)算的結(jié)果信息，各內(nèi)部節(jié)點(diǎn)分歧時(shí)間、平均進(jìn)化速率，sigma2值和各分支進(jìn)化速率等）
samplefreq = 10 ：每10次迭代取樣一次
nsample = 20000 ：當(dāng)取樣次數(shù)達(dá)到該次數(shù)時(shí)，則取樣結(jié)束（程序也運(yùn)行結(jié)束）
burnin= 2000，samplefreq= 10，nsample= 20000：總的意思是程序會(huì)去除掉（burnin)前2000次迭代的結(jié)果，然后每10次迭代取一次樣，直到取樣次數(shù)達(dá)到20000次，因此MCMC會(huì)運(yùn)算2000+10*20000次迭代。一般來(lái)說(shuō)，程序需要進(jìn)行10000-20000次取樣，才能獲得較好的數(shù)據(jù)總結(jié)。一般來(lái)說(shuō)，如果想通過(guò)增加MCMC長(zhǎng)度來(lái)提高收斂效果，一般是通過(guò)增加samplefreq，但保持nsample在一個(gè)比較合理的范圍。
別人的博客曾進(jìn)行過(guò)一個(gè)其個(gè)人的總結(jié)，照抄如下：
數(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
在Linux系統(tǒng)下運(yùn)行時(shí)，將序列文件、樹(shù)文件和控制文件放在同一個(gè)路徑下，在該路徑下運(yùn)行程序：
which mcmctree? #查找mcmctree的路徑
~/miniconda3/bin/mcmctree mcmctree.ctl? #~/miniconda3/bin/mcmctree是程序的絕對(duì)路徑，mcmctree.ctl是控制文件的名字
用approximate likelihood calculation進(jìn)行物種分化時(shí)間的估算
在對(duì)較大的數(shù)據(jù)進(jìn)行l(wèi)ikelihood function計(jì)算時(shí)，計(jì)算成本昂貴，耗時(shí)長(zhǎng)，下面介紹針對(duì)較大比對(duì)數(shù)據(jù)的approximate likelihood計(jì)算的方法。
總共分兩步，在第一步中，枝長(zhǎng)、gradient、Hessian由最大似然法估算出，在第二步中，利用gradient、Hessian進(jìn)行分化時(shí)間的估算，這一步用Taylor expansion的方法創(chuàng)造出近似于最大似然法的功能。
總的操作步驟如下：
新建一個(gè)文件夾Hessian，將樹(shù)文件、序列文件和控制文件拷貝過(guò)去，控制文件中的usedata = 3，然后運(yùn)行程序。程序運(yùn)行結(jié)束后生成文件中有一個(gè)out.BV文件
再建立一個(gè)新文件夾approx01，將上一步中的out.BV文件、樹(shù)文件、序列文件和控制文件拷貝過(guò)去，out.BV文件重命名為in.BV,控制文件中的usedata = 3改為usedata = 2。
