使用PAML進行分歧時間計算
PAML軟件中的mcmctree命令可以使用Bayesian方法估算物種分歧時間。
對程序輸入帶有校准點信息的有根樹、多序列比對結果，即可得到進化樹各內部節點95%置信區間分歧時間信息。
PAML軟件也能利用帶有枝長信息的樹，快速進行分歧時間計算。
1. 輸入文件（1/3）：帶有校准點的有根樹 input.trees 輸入一個帶有校准點的有根樹文件，示例如下：
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格式） input.txt
7 3331 human AGACTGTTAGCAGCGCCGGGCAACTCCCTACTCAA... chimpanzee AGACTGTTGGCAACGTCGGGCAACTCCCCGCCCAA... bonobo AGACTGTTGGCAACGCCGGGCAACTCCCCGCCCAA... gorilla AGATTGTTAGCAACGTCGGGTAACCCCCCACTCAA... orangutan AGGCTACTAACAGCGCCGGACGATTCCTCGCCTAA... sumatran AGACTACTAACAGCGCCGGGCGATTCCTCACCCAA... gibbon AGACTATTGACAACGTCGGGCAACTCTCTACTCAA... 7 3331 human AAATTCCTTCCCTTGTCCCTTTTTTCCTTTCATTA... chimpanzee AAATTCCTCCCCTTGTCCCTTTTTTCCTTTCATTA... bonobo AAATTCCTCCCCTTGTCCCTTTTTTCCTTTCATTA... gorilla AAATTCCTTCCCTTGTCCCTTTTTTCCTTTCATTA... orangutan AAATTCCTCCCCTTGTCCCTTTTTTCCTTTCATTA... sumatran AAGTTCCTTCCCTTGTCCCTTTTTTCCTTTCATTA... gibbon AAATTCCTCCCCTTGTCCCTCTTTTCCTTTCATTA... 7 3331 human CATGCTACTCCACACACCAAGCTATCTAGCCTCCC... chimpanzee CATACTACTCCACACACCAAACTACCTAGCCTCCC... bonobo CATGCTACTCCACACACCAAGCTACCTAGTCTCCC... gorilla CATACTACTCCACACACCAAATCATCTAGCCTCCC... orangutan CATACCACTCCACACCCTATACCATCCAACTTCCC... sumatran CATATCACTCCAAACCCCAAACCATCCAGCCTCCC... gibbon CATACTACTCCATACACCAAATTATCCAACTCCCC...
PAML要求輸入的Phylip格式，其物種名和后面的序列之間至少間隔兩個空格（是為了允許物種名的屬名和種名之間有一個空格）。
當然，也可以輸入一個多序列比對結果進行分析。在軟件說明文檔中，是密碼子三個位點分別的多序列比對結果，程序可以分別計算密碼子各個位點的變異速率。
此外，也可以輸入密碼子或氨基酸序列比對信息，則需要再配置文件中修改相應的參數來表明數據類型。若使用氨基酸序列來進行分析，由於mcmctree命令不能選擇較好的氨基酸替換模型進行分析，需要自己手動運行codeml進行分析后，在生成中間文件用於運行mcmctree，比較麻煩。因此，推薦按上述方法使用密碼子三位點的鹼基序列作為輸入信息進行PAML分析。
3. 輸入文件（3/3）：mcmctree命令的配置文件 mcmctree.ctl
軟件難點在於配置文件的理解，示例如下：
*輸入輸出參數：
seed = -1 *設置一個隨機數來運行程序。若設置為-1，表示使用系統當前時間為隨機數。
seqfile = input.txt *設置輸入的多序列比對文件路徑
treefile = input.trees *設置輸入的帶校准點信息有根樹的文件路徑
mcmcfile = mcmc.txt *設置輸出的mcmc信息文本文件，可以用Tracer軟件查看
outfile = out.txt *設置輸出文件路徑，該文件中記錄了超度量樹和進化速率等參數信息。
*數據使用說明參數：
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 *當外部節點由取樣時間時使用該參數進行設置，同時該參數也設置了時間單位。具體數據示例請見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結果文件；可以在此處手動修改tmp0001.ctl配置文件，例如，修改其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命令進行分歧時間計算
程序運行命令：
mcmctree mcmctree.ctl
結果文件解釋
程序在運行過程中，會在屏幕生生成一些信息。
比較耗時間的步驟主要在於取樣的百分比進度：
第一列：取樣的百分比進度。
第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信息，可以用於收斂分析。
在當前目錄下，生成幾個主要結果：
FigTree.tre 生成含有分歧時間的超度量樹文件
mcmc.txt MCMC取樣信息，包含各內部節點分歧時間、平均進化速率、sigma2值等信息，可以在Tracer軟件中打開。
通過查看各參數的ESS值，若ESS值大於200，則從一定程度上表示MCMC過程能收斂，結果可靠。
out.txt 包含由較多信息的結果文件。例如，各鹼基頻率、節點命名信息、各節點分歧時間、進化速率和進化樹等。
使用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命令，要求額外多輸入一個名為 FixedDsClock23.txt 的文件。該文件中存放了相應的帶有枝長信息的系統發育樹。該文件中系統發育樹的拓撲結構要和input.trees一致。推薦使用baseml命令輸入input.trees和多序列比對結果計算枝長。
FixedDsClock23.txt 文件內容示例：
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命令進行分歧時間計算：
infinitesites mcmctree.ctl使用infinitesites進行分歧時間計算時，程序要求輸入多序列比對文件。雖然程序讀取了序列信息，但在計算時會忽略其序列信息。
其它注意事項
(1) 如何檢測MCMC的結果是否達到收斂狀態？
收斂的意思，即經過很多次迭代后，得到的MCMC樹的各枝長趨於一個定值，變動非常小。
於是，最直接的檢測方法是：分別使用不同的Seed值進行mcmctree或infinitesites進行兩次或多次分析，然后比較兩個結果樹是否一致，實際就是比較樹文件中各內部節點的Height值（分歧時間 / Posterior time）。
比較方法可以有多種:
(a) 結果文件mcmc.txt中每行記錄了一個取樣的結果，包含各內部節點的Posterior time值，即進化樹的Height值，對該mcmc.txt文件進行分析，得到每個內部節點的Posterior mean time值。然后作圖 （官方說明文檔第6頁圖示），橫坐標表示第一次運行的各內部節點的分歧時間均值、縱坐標表示第二次運行的各內部節點的分歧時間均值。該散點圖要表現出非常明顯的對角線，才認為達到收斂。這時可以考慮計算其相關系數來判斷是否符合線性。
(b) 我覺得第一種官方文檔給的方法比較麻煩且是否符合對角線沒有明確的定義，於是就直接比較兩次結果樹文件中的各枝長。計算各枝長總的偏差百分比，當偏差百分比低於0.1%，則認為兩次結果非常吻合，差異低於0.1%，認為達到收斂。
(c) 此外，還可以使用Tracer分析mcmc.txt文件，檢測其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掉的迭代次數和用於計算結果的有效迭代次數。由於需要根據有效迭代數據來進行均值計算，若記錄每次迭代的參數，則生成的mcmc.txt文件行數過多，文件太大，匯總時也消耗計算。於是設置沒隔一定迭代輪數取樣一次，這樣生成的mcmc.txt文件不會太大，對最后的均值計算影響不大。
(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
那如何查兩兩物種間的化石分化時間？
通過網站http://www.timetree.org/ ，該網站根據多篇文獻支持提供兩兩物種間的分化時間，並給出置信度范圍，單位是Mya(million years ago)。
另外一個可以查詢分化時間的網站是https://fossilcalibrations.org/ ，可以互相參考。
在已經構建好的樹中添加化石時間，主要形式是 '>0.22<6.75'，'>'后接分化時間的最小值，'<'后接最大值，如下實例：
當樣本比較多時，要考慮不同距離的節點，只有一個節點估算的可能不准確，可以提供三到五個的節點（物種）的分化時間。
03分歧時間估計
這里主要介紹3種程序：
1. 使用 r8s（https://sourceforge.net/projects/r8s/）估計分歧時間，速度很快，一般時間<10min
准備配置文件r8s.ctl：
運行程序：
2. 使用multidivtime（https://brcwebportal.cos.ncsu.edu/thorne/multidivtime.html） 估計分歧時間，時間~5hours
-> 先利用PAML baseml估計分支的核苷酸頻率、轉換/顛換率等參數
再利用paml2modelinf將baseml輸出文件轉換為estbranches輸入文件
利用estbranche對樹分支長度進行最大似然率估計
輸入文件是前一步的輸出結果
-> 運行multidivtime進行分歧時間估計
利用前一步的輸出，准備配置文件divtime.ctl，運行程序
3. 使用 PAML mcmctree（http://abacus.gene.ucl.ac.uk/software/paml.html） 估計分歧時間，強烈推薦，時間~20hours
准備配置文件 mcmctree.ctl
運行程序：
計算完成后，就獲得了帶分歧時間的進化樹了，后續可以對進化樹進行一系列美化操作。
來源
http://www.chenlianfu.com/?p=2974
https://www.biomart.cn/65423/news/2908793.htm