雖然不是學(xué)生物或者醫(yī)學(xué)的，但是聽了很多生命科學(xué)的報告撑刺。關(guān)于生命的起源鹉胖、進(jìn)化——從個體到宇宙，生命是如何從無到有，又如何實現(xiàn)了現(xiàn)在的千姿百態(tài)——是生命科學(xué)數(shù)百年來一直在研究的問題甫菠。其實總結(jié)來說挠铲，無非要拋出常見的三聯(lián)問：
記得以前聽一個老師的報告，問過大家一個問題：如果你碰到一個無所不知寂诱，文明遠(yuǎn)遠(yuǎn)高于我們的外星人拂苹，你可以問他一個問題，會問什么痰洒？我記得一個深刻的答案是：會問人類是不是人工智能瓢棒？? ?好像有個科幻小說里面也有這方面的描述，問題例如：宇宙的目的是什么丘喻？意識的目的是什么脯宿？靈魂的目的是什么？等等等泉粉。
利用進(jìn)化樹將穿越數(shù)萬年的物種進(jìn)化濃縮于一圖连霉，更是生物學(xué)研究中不可或缺的必備工具，也是我們領(lǐng)域天天打交道的事情嗡靡。
由于不同的基因或DNA片段的進(jìn)化速率存在較大的差異窘面，我們就可以利用這些基因或DNA片段來估算各個分類水平上有機(jī)體之間的進(jìn)化關(guān)系。
分化時間是當(dāng)前宏觀進(jìn)化分析的一個熱點叽躯，以某一個或某幾個特定類群的化石時間作為參照，通過基因序列間的分歧程度以及分子鐘來估計分支間的分歧時間肌括，同時計算系統(tǒng)發(fā)育樹上其它節(jié)點的發(fā)生時間点骑，從而推斷相關(guān)類群的起源和不同類群的分歧時間。
根據(jù)分子鐘理論谍夭，基因序列中密碼子隨著時間的推移以幾乎恒定的比例相互替換黑滴，即具有恒定的演化速率，兩個物種之間的遺傳距離將與物種的分化時間成正比紧索。我們通常采用單拷貝基因家族中的四重簡并位點來估算分子鐘（替換速率）以及物種間的分化時間袁辈，密碼子中的四重簡并位點由于第三堿基不改變所編碼的氨基酸，屬于中性進(jìn)化珠漂，其中中性替換速率一般用每個位點每年的變異數(shù)來衡量晚缩。
利用貝葉斯統(tǒng)計或其他方法估計物種分歧時間的程序很多，如R8S媳危、MCMCTREE荞彼、MULTIDIVTIME、BEAST等待笑，通過不同的策略將化石時間信息整合到一個系統(tǒng)發(fā)育樹中鸣皂，從而計算得到Divergence time Tree。
下面就是一個帶有分化時間的進(jìn)化樹（分支帶有物種的分化時間，并添加了分化時間的置信度范圍寞缝，藍(lán)色數(shù)字為分化時間）：
=====構(gòu)建系統(tǒng)發(fā)育樹=====
前面講過很多如何構(gòu)建發(fā)育樹的文章癌压，比如基于orthologous的基因利用muscle/Gblocks/RAxml構(gòu)建發(fā)育樹。
muscle -in ../OG0012924.fa -out OG0012924.fa
Gblocks OG0012924.fa -t=p -b4=5 -b5=h -e=.gb
raxmlHPC-PTHREADS-SSE3? --f a --x 12345 --p 12345 --# 100 --m PROTGAMMALGX --s OG0012924.fa.gb --n ex -T 20?
下面我們測試用不同的幾個工具嘗試去計算分化時間：
=======兩兩物種的分化時間=======
通過網(wǎng)站：http://www.timetree.org/ 荆陆，可以查兩兩物種間的分化時間滩届。有個單位是：1Mya ?= 1 million years ago。
另外一個可以查詢分化時間的網(wǎng)站：https://fossilcalibrations.org/
==================r8s=======================
準(zhǔn)備輸入文件：
前面是輸入的傳統(tǒng)方式的tree慎宾。
BLFORMAT
length ????設(shè)置樹的枝長單位丐吓。若枝長單位是位點替換率，則其值為 persize趟据，則枝長 * 序列長度 = 替換數(shù)券犁；若枝長單位是替換數(shù)，則該參數(shù)值為 total汹碱。默認(rèn)參數(shù)是 total粘衬。
nsites ????設(shè)置多序列比對的序列長度。
ultrametric 是否是超度量樹咳促，一般進(jìn)化樹選 no稚新。默認(rèn)參數(shù)是 no。
MRCA
該命令用來設(shè)置內(nèi)部節(jié)點名跪腹。示例中設(shè)置了 tree 和 scer 的 most recent common ancestor 的節(jié)點名為 asc褂删。
FIXAGE
該命令用來設(shè)置 MRCA 指定的節(jié)點的分歧時間，使用該指定的分歧時間作為校正冲茸，來預(yù)測其它節(jié)點的分歧時間屯阀。
r8s 需要至少有一個內(nèi)部節(jié)點進(jìn)行 fixage。
CONSTRAIN
該命令用來約束 MRCA 指定的節(jié)點的分歧時間轴术，設(shè)置該節(jié)點允許的最大或最小分歧時間难衰。
DIVTIME
該命令用來進(jìn)行分歧時間和替換速率計算《涸裕總共有 4 種計算方法（LF | NPRS | PL）和 3 種數(shù)學(xué)算法（Powell | TN | QNEWT）盖袭。一般常用且較優(yōu)，是使用 PL 和 TN彼宠。
但是使用 PL 方法鳄虱，則需要設(shè)置參數(shù) smoothing 的值。通過設(shè)置多個 smoothing 的值來進(jìn)行一些計算凭峡，選擇最優(yōu)的值即可醇蝴。一般情況下，該值位于 1~1000 能得到一個最佳（最低）的分值想罕。
divtime method=PL crossv=yes cvstart=0cvinc=1 cvnum=18;
上述命令悠栓，是設(shè)置 smoothing 的值從 1, 10, 100, 1000 ... 1e17, 來計算霉涨，最后得到最佳的 smoothing值。
若使用 fixage 對 2 個節(jié)點的分歧時間進(jìn)行了固定惭适，則可以運(yùn)行命令：
divtime method=PL crossv=yesfossilfixed=yes cvstart=0 cvinc=1 cvnum=18;
若使用 fixage 對 1 個節(jié)點進(jìn)行分歧時間固定笙瑟，同時使用 constrain 對 2 個節(jié)點進(jìn)行了約束，則可以運(yùn)行命令：
divtime method=PL crossv=yesfossilconstrained=yes cvstart=0 cvinc=1 cvnum=18;
得到最優(yōu)的 smoothing 值后癞志，使用 set 進(jìn)行設(shè)置往枷，然后運(yùn)行 divtime 進(jìn)行分歧時間和替換速率的計算。
SHOWAGE
使用該命令能得到各個節(jié)點的分歧時間和替換速率凄杯。
DESCRIBE
使用該命令得到進(jìn)化樹的圖和文字結(jié)果错洁。其 plot 參數(shù)如下：
cladogram ??得到分支樹的圖，圖上有各個節(jié)點的編號戒突，和 showage 的結(jié)果結(jié)合觀察屯碴。
phylogram ??得到進(jìn)化樹的圖，枝長表示替換數(shù)膊存。
chronogram ?得到超度量樹的圖导而，枝長表示時間。
ratogram ???得到樹圖隔崎，枝長表示替換速率今艺。
phylo_description ??? 得到樹的ASCII文字結(jié)果，枝長表示替換數(shù)爵卒。
chrono_description ?? 得到樹的ASCII文字結(jié)果虚缎，枝長表示時間。
rato_description ???? 得到樹的ASCII文字結(jié)果钓株，枝長表示替換速率遥巴。
node_info???得到節(jié)點的信息表格
然后運(yùn)行輸入文件：/gpfs03/home/jingjing/software/r8s1.81/src/r8s -b -f r8s_ctl_file.txt > r8s_tmp.txt
=================MCMCtree==================
使用 PAML mcmctree（http://abacus.gene.ucl.ac.uk/software/paml.html） 估計分歧時間。
下載安裝：
wget http://abacus.gene.ucl.ac.uk/software/paml4.9j.tgz
tar xf paml4.9j.tgz
cd?paml4.9j
rm bin/*.exe
cd src
make -f Makefile
rm *.o
mv baseml basemlg codeml pamp mcmctree evolver yn00chi2? ../bin
輸入文件主要包括3個文件：
序列比對文件（nucl或者pep都可以）
帶有化石校準(zhǔn)的系統(tǒng)發(fā)育樹
mcmctree.ctl享幽；控制文件
(1) 序列比對文件
paml4.9j/examples/DatingSoftBound/mtCDNApri123.txt，是一個phlip格式拾弃，PAML要求輸入的Phylip格式值桩，其物種名和后面的序列之間至少間隔兩個空格（是為了允許物種名的屬名和種名之間有一個空格）。
該示例文件分為3個區(qū)塊豪椿，密碼子在3個位點的多序列比對結(jié)果奔坟。
此外，也可以輸入密碼子或氨基酸序列比對信息搭盾，則需要再配置文件中修改相應(yīng)的參數(shù)來表明數(shù)據(jù)類型咳秉。若使用氨基酸序列來進(jìn)行分析，由于mcmctree命令不能選擇較好的氨基酸替換模型進(jìn)行分析鸯隅，需要自己手動運(yùn)行codeml進(jìn)行分析后澜建，在生成中間文件用于運(yùn)行mcmctree向挖，比較麻煩。因此炕舵，推薦按上述方法使用密碼子三位點的堿基序列作為輸入信息進(jìn)行PAML分析何之。
(2) 帶有校準(zhǔn)點的有根樹
示例文件如下：
第一行，7個物種咽筋，1顆樹溶推，中間用空格分開，并且不能有枝長奸攻。
第二行蒜危，其中包含有校準(zhǔn)點信息。校準(zhǔn)點信息一般指95%HPD（Highest Posterior Density）對應(yīng)的置信區(qū)間睹耐，單位為100個百萬年(MY)唤殴。尾部一定要有分號滞谢。
這棵樹有兩種化石校準(zhǔn)，一種是人類黑猩猩的最近共同祖先：’>.06<.08’，另一種是大猩猩的最近共同祖先：’>.12<.16’孕似。時間單位是100百萬年(Myr)。因此匠抗，第一次校準(zhǔn)將人類黑猩猩的共同祖先限制在6-8 myr之前澳叉。mcmctree中的校準(zhǔn)界限是軟的，并且很小的概率（默認(rèn)為0.025）可能違反界限窄绒。注意贝次，樹的根部沒有化石校準(zhǔn)。MCMCTree需要在樹的根上進(jìn)行校準(zhǔn)彰导，當(dāng)該校準(zhǔn)不存在于樹文件中時蛔翅，必須在控制文件中指定它(variable RootAge)。
(3) 配置文件：
seed = -1?# 表示使用系統(tǒng)當(dāng)前時間為隨機(jī)數(shù)
seqfile = mtCDNApri123.txt # 多序列比對文件
treefile = mtCDNApri.trees # 帶有校準(zhǔn)信息有根樹
mcmcfile = mcmc.txt #輸出mcmc信息文本文件位谋，可以用Tracer軟件查看
outfile = out.txt # 輸出文件山析，記錄了超度量樹和進(jìn)化速率等參數(shù)信息
ndata = 3 # 輸入的多序列比對的數(shù)據(jù)個數(shù)，這里密碼子3個位置的數(shù)據(jù)掏父；如? ? ? ? ? ? ? ? ? ? ? 果有一個笋轨，則設(shè)置為1
seqtype = 0??? * 0: nucleotides; 1:codons; 2:AAs #數(shù)據(jù)類型
usedata = 1 ?? * 0: no data; 1:seq like; 2:normalapproximation; 3:out.BV ?? ? ? ? ? ? ? ? ? ? ? ?(in.BV) # 設(shè)置是否利用多序列比對的數(shù)據(jù)：
#0，表示不使用多序列比對數(shù)據(jù)赊淑，則不會進(jìn)行l(wèi)ikelihood計算爵政，雖然能得到mcmc樹且計算速度飛快，但是其分歧時間結(jié)果是有問題的陶缺；
#1钾挟，表示使用多序列比對數(shù)據(jù)進(jìn)行l(wèi)ikelihood計算，正常進(jìn)行MCMC饱岸，是一般使用的參數(shù);?
#2掺出，進(jìn)行正常的approximationlikelihood分析徽千，此時不需要讀取多序列比對數(shù)據(jù)，直接讀取當(dāng)前目錄中的in.BV文件蛛砰。該文件是使用usedata = 3參數(shù)生成的out.BV文件重命名而來的罐栈。
#此外，由于程序BUG泥畅，當(dāng)設(shè)置usedata =2時荠诬，一定要在改行參數(shù)后加 *，否則程序報錯 Error: file name empty..?
#3位仁，程序利用多序列比對數(shù)據(jù)調(diào)用baseml/codeml命令對數(shù)據(jù)進(jìn)行分析柑贞，生成out.BV文件。由于mcmctree調(diào)用baseml/codeml進(jìn)行計算的參數(shù)設(shè)置可能不太好（特別時對蛋白序列進(jìn)行計算時）聂抢，
#推薦自己修該軟件自動生成的baseml/codeml配置文件钧嘶，然后再手動運(yùn)行baseml/codeml命令，再整合其結(jié)果文件為out.BV文件琳疏。
clock = 2???* 1: global clock; 2: independent rates; 3: correlated rates # (1) 全局分子鐘,適用于近緣物種（兩物種分化時間小于20個million分化時間或者序列差異小于5%即近緣物種）
# (2) 樹枝上進(jìn)化速率服從獨立同分布的對數(shù)正態(tài)分布（推薦使用）
# (3) 樹枝上的進(jìn)化速率自相關(guān)有决。
RootAge = '<1.0'? * safe constraint on root age, used if nofossil for root. # 設(shè)置root節(jié)點的分歧時間，一般設(shè)置一個最大值空盼。
model = 0???* 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
# 堿基替換模型：
???#（0）無參數(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
# 核酸序列中不同位點苟呐，其進(jìn)化速率不一致，其變異速率服從GAMMA分布俐筋。一般設(shè)置GAMMA分布的alpha值為0.5牵素。\
# 若該參數(shù)值設(shè)置為0，則表示所有位點的進(jìn)化速率一致澄者。此外笆呆，當(dāng)userdata = 2時，alpha单起、ncatG、alpha_gamma嘀倒、kappa_gamma這4個GAMMA參數(shù)無效局冰。因為userdata = 2時，不會利用到多序列比對的數(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è)置出生率产雹、死亡率和取樣比例蔓挖。若輸入有根樹文件中的時間單位發(fā)生改變拷获，則需要相應(yīng)修改出生率和死亡率的值邪财。
例如，時間單位由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è)置序列中所所有位點平均[堿基/密碼子/氨基酸]替換率的Dirichlet-GAMMA分布參數(shù)：
sigma2_gamma = 1 10 1 ?* gammaDir prior for sigma^2 ???(for clock=2 or 3) #設(shè)置所有位點進(jìn)化速率取對數(shù)后方差（sigma的平方）的Dirichlet-GAMMA分布參數(shù)：alpha=1胶滋、beta=10、初始方差值為1镀迂。
# 當(dāng)clock參數(shù)值為1時丁溅，表示使用全局的進(jìn)化速率，各分枝的進(jìn)化速率沒有差異探遵，即方差為0窟赏，該參數(shù)無效；
# 當(dāng)clock參數(shù)值為2時箱季，若修改了時間單位涯穷，該參數(shù)值不需要改變；
# 當(dāng)clock參數(shù)值為3時藏雏，若修改了時間單位拷况，該參數(shù)值需要改變。
finetune = 1: .1 .1 .1 .1 .1 .1 * auto (0or 1): times, musigma2, rates, mixing, paras, FossilErr #冒號前的值設(shè)置是否自動進(jìn)行finetune掘殴，一般設(shè)置成1赚瘦，然后程序自動進(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é)點分歧時間、平均進(jìn)化速率病瞳、sigma2值）揽咕；
#2，打印所有信息套菜。
burnin = 2000 #將前2000次迭代burnin后亲善，再進(jìn)行取樣（即打印出該次迭代計算的結(jié)果信息，各內(nèi)部節(jié)點分歧時間逗柴、平均進(jìn)化速率蛹头、sigma2值和各分支進(jìn)化速率等）
sampfreq = 10 #每10次迭代則取樣一次。
nsample = 20000 # 當(dāng)取樣次數(shù)達(dá)到該次數(shù)時，則取樣結(jié)束（程序也將運(yùn)行結(jié)束）掘而。
(4) 運(yùn)行：
mcmctree mcmctree.ctl
在運(yùn)行過程當(dāng)中，屏幕會生成一些信息于购。
其中：
第一列：取樣的百分比進(jìn)度袍睡。
第2~6列：參數(shù)的接受比例。一般肋僧，其值應(yīng)該在30%左右斑胜。20~40%是很好的結(jié)果， 15~70%是可以接受的范圍嫌吠。若這些值在開始時變動較大止潘，則表示burnin數(shù)設(shè)置太小。
第7~x列：各內(nèi)部節(jié)點的平均分歧時間辫诅，第7列則是root節(jié)點的平均分歧時間凭戴。若有y個物種，則總共有y-1個內(nèi)部節(jié)點炕矮，從第7列開始后的y-1列對應(yīng)這些內(nèi)部節(jié)點么夫。
倒數(shù)第3~x列：r_left值。若輸入3各多序列比對結(jié)果肤视，則有3列档痪。x列的前一列是中劃線 - 。
倒數(shù)第1~2列：likelihood值和時間消耗邢滑。
注意：查看每列中的值（接受比例腐螟、時間和速率）很重要。它們應(yīng)該在MCMC運(yùn)行期間保持穩(wěn)定困后。如果驗收比例變化太大（特別是在開始時）乐纸，這意味著預(yù)燒時間不夠長，因此應(yīng)增加控制文件中的預(yù)燒變量操灿。
在輸出最后锯仪，給出各個內(nèi)部節(jié)點的分歧時間(t)、平均變異速率(mu)趾盐、變異速率方差(sigma2)和r_left的Posterior信息：
均值(mean)庶喜、95%雙側(cè)置信區(qū)間(95% Equal-tail CI)和95% HPD置信區(qū)間(95% HPD CI)等信息。
此外救鲤，第二列給出了各個內(nèi)部節(jié)點的Posterior mean time信息久窟，可以用于收斂分析。jnode是程序multipvtime使用的節(jié)點號本缠，由JeffThorne編寫斥扛。
倒數(shù)第二列值：比如第8節(jié)點；0.0436=0.1850-0.1414。
(5) 結(jié)果分析：
運(yùn)行會得到以下4個文件
FigTree.tre 生成含有分歧時間的超度量樹文件稀颁。適用于Andrew Rambaut的figtree程序（http://tree.bio.ed.ac.uk/software/figtree/）
節(jié)點處的位分化時間芬失，比如 chimpanzee 和bonobo分化時間為3.01 MYA。
mcmc.txt MCMC取樣信息匾灶，包含各內(nèi)部節(jié)點分歧時間棱烂、平均進(jìn)化速率、sigma2值等信息阶女，可以在Tracer軟件中打開颊糜。通過查看各參數(shù)的ESS值，若ESS值大于200秃踩，則從一定程度上表示MCMC過程能收斂衬鱼，結(jié)果可靠。
在我們的示例中憔杨，該文件有50列(如果print=2)或14列(如果print=1)和20002行鸟赫。第一列是樣本的生成(迭代)數(shù)，接下來的48列(print=2)或12列(print=1)分別對應(yīng)于所分析的48個或12個參數(shù)消别；最后一列是似然值(第50列惯疙，如果print=2；或第14列妖啥，如果print=1)霉颠。行數(shù)與MCMC期間采集的樣本數(shù)相對應(yīng)【Ｊ“mcmc.txt”文件適用于Tracer程序分析（http://beast.bio.ed.ac.uk/tracer）
out.txt 包含由較多信息的結(jié)果文件蒿偎。例如，各堿基頻率怀读、節(jié)點命名信息诉位、各節(jié)點分歧時間、進(jìn)化速率和進(jìn)化樹等菜枷。
所以苍糠，我們只需要把比對文件和生成的tree自己寫個腳本轉(zhuǎn)化成輸入需求就行了。
本文使用 文章同步助手 同步
