時(shí)間樹

基因序列中的密碼子隨著時(shí)間的推移以幾乎恒定的比例進(jìn)行替換,具有恒定的演化速率,因此兩個(gè)物種的遺傳距離與物種分歧時(shí)間成正比榔至。序列比對(duì)時(shí),我們通常會(huì)利用單拷貝直系同源基因多序列比對(duì)文件用于構(gòu)建進(jìn)化樹欺劳。為了獲取更多的遺傳信息唧取,我們會(huì)將蛋白質(zhì)比對(duì)轉(zhuǎn)換成密碼子比對(duì),這時(shí)我們可以將這些密碼子比對(duì)文件串聯(lián)起來形成超級(jí)比對(duì)文件划提。而我們繪制時(shí)間樹枫弟,需要穩(wěn)定的核酸替換率,而三位密碼子所受到的選擇壓力各不相同鹏往,因此有必要使用 PartitionFinder2 將串聯(lián)而成的多序列比對(duì)文件進(jìn)行分區(qū)(可以分成三個(gè)部分)淡诗。將分區(qū)的比對(duì)文件可以分別計(jì)算進(jìn)化模型,RaxML伊履,MrBayes韩容,Beastiqtree 都支持分區(qū)建樹唐瀑。而推導(dǎo)時(shí)間樹的程序 mcmctree 和 BEAST 都支持分區(qū)信息群凶。

PROGRAM SUPPORT MULTITHREADING INPUT FILE(Fossil calibration information) TIME EASY INSTALLATION
Mega YES A tree without branching information,Multi-sequence comparison file About 1 hour YES
R8s NO A tree with branching information(root) Within 10 minutes NO
MCMCtree NO A tree without branching information(root),Multi-sequence comparison file More than four days NO
BEAST YES Multi-sequence comparison file More than four days YES

一.mcmctree

三個(gè)輸入文件(最好是核酸比對(duì)文件):

  • 1.序列比對(duì)文件
  • 2.帶有校準(zhǔn)點(diǎn)的有根樹
  • 3.配置文件mcmctree.ctl

1.序列比對(duì)文件

2.帶有校準(zhǔn)點(diǎn)的有根樹

Cl:置信區(qū)間(Confidence Interval)
MYA:百萬年以前(Million Years Ago)的縮寫
CI: (6.0 - 6.5 MYA)表示置信區(qū)間為從600萬年前到650萬年前的時(shí)間段。換句話說介褥,這個(gè)置信區(qū)間的長度是0.5百萬年座掘,也就是50萬年。


化石校準(zhǔn)點(diǎn)
 7  1

((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);

這棵帶有校準(zhǔn)點(diǎn)的有根樹中時(shí)間單位是:一億年 --100個(gè)百萬年
上圖這棵樹中:人類(human)和黑猩猩(chimpanzee)之間的分歧時(shí)間在600萬年-800萬年之間
要注意的點(diǎn)是:樹在根上沒有化石校正點(diǎn)

3.配置文件mcmctree.ctl

使用PAML進(jìn)行分歧時(shí)間計(jì)算 | 陳連福的生信博客 (chenlianfu.com)
估算系統(tǒng)樹分歧時(shí)間 —— paml.mcmctree,r8s | 生信技工 (yanzhongsino.github.io)
【比較基因組】從進(jìn)化樹到分化時(shí)間 - 簡書 (jianshu.com)

***輸入輸出參數(shù):**
        seed = -1           *設(shè)置一個(gè)隨機(jī)數(shù)來運(yùn)行程序柔滔。若設(shè)置為-1溢陪,表示使用系統(tǒng)當(dāng)前時(shí)間為隨機(jī)數(shù)。
     seqfile = input.txt    *設(shè)置輸入的多序列比對(duì)文件路徑
    treefile = input.trees  *設(shè)置輸入的帶校準(zhǔn)點(diǎn)信息有根樹的文件路徑 
    mcmcfile = mcmc.txt     *設(shè)置輸出的mcmc信息文本文件睛廊,可以用Tracer軟件查看
     outfile = out.txt      *設(shè)置輸出文件路徑形真,該文件中記錄了超度量樹和進(jìn)化速率等參數(shù)信息。

***數(shù)據(jù)使用說明參數(shù):**
       ndata = 3         *設(shè)置輸入的多序列比對(duì)的數(shù)據(jù)個(gè)數(shù)。這里指密碼子3個(gè)位置的數(shù)據(jù)咆霜。
     seqtype = 0         *設(shè)置多序列比對(duì)的數(shù)據(jù)類型:0邓馒,表示核酸數(shù)據(jù);1蛾坯,表示密碼子比對(duì)數(shù)據(jù)光酣;2,表示氨基酸數(shù)據(jù)脉课。根據(jù)不同的數(shù)據(jù)類型救军,程序調(diào)用相應(yīng)的baseml或codeml進(jìn)行相應(yīng)的參數(shù)計(jì)算。
     usedata = 1         *設(shè)置是否利用多序列比對(duì)的數(shù)據(jù):0倘零,表示不使用多序列比對(duì)數(shù)據(jù)唱遭,則不會(huì)進(jìn)行l(wèi)ikelihood計(jì)算,雖然能得到mcmc樹且計(jì)算速度飛快呈驶,但是其分歧時(shí)間結(jié)果是有問題的拷泽;1,表示使用多序列比對(duì)數(shù)據(jù)進(jìn)行l(wèi)ikelihood計(jì)算袖瞻,正常進(jìn)行MCMC司致,是一般使用的參數(shù); 2,進(jìn)行正常的approximation likelihood分析虏辫,此時(shí)不需要讀取多序列比對(duì)數(shù)據(jù)蚌吸,直接讀取當(dāng)前目錄中的in.BV文件锈拨。該文件是使用usedata = 3參數(shù)生成的out.BV文件重命名而來的砌庄。此外,由于程序BUG奕枢,當(dāng)設(shè)置usedata = 2時(shí)娄昆,一定要在改行參數(shù)后加 *,否則程序報(bào)錯(cuò) Error: file name empty.. 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è)置可能不太好(特別時(shí)對(duì)蛋白序列進(jìn)行計(jì)算時(shí))谷浅,推薦自己修該軟件自動(dòng)生成的baseml/codeml配置文件扒俯,然后再手動(dòng)運(yùn)行baseml/codeml命令,再整合其結(jié)果文件為out.BV文件一疯。
   cleandata = 0         *設(shè)置是否移除不明確字符(N撼玄、?墩邀、W掌猛、R和Y等)或含以后gap的列后再進(jìn)行數(shù)據(jù)分析:0,不需要眉睹,但在序列兩兩比較的時(shí)候荔茬,還是會(huì)去除后進(jìn)行比較废膘;1,需要慕蔚。 
       clock = 2         *設(shè)置分子種方法:1丐黄,global clock方法,表示所有分支進(jìn)化速率一致孔飒;2孵稽,independent rates方法,各分支的進(jìn)化速率獨(dú)立且進(jìn)化速率的對(duì)數(shù)log(r)符合正態(tài)分布; 3十偶,correlated rates方法菩鲜,和方法2類似,但是log(r)的方差和時(shí)間t相關(guān)惦积。
*    TipDate = 1 100     *當(dāng)外部節(jié)點(diǎn)由取樣時(shí)間時(shí)使用該參數(shù)進(jìn)行設(shè)置接校,同時(shí)該參數(shù)也設(shè)置了時(shí)間單位。具體數(shù)據(jù)示例請(qǐng)見examples/TipData文件夾狮崩。
     RootAge = '<1.0'    *設(shè)置root節(jié)點(diǎn)的分歧時(shí)間蛛勉,一般設(shè)置一個(gè)最大值。
 
***位點(diǎn)替換模型參數(shù):**
       model = 4         *設(shè)置堿基替換模型:0睦柴,JC69诽凌;1,K80坦敌;2侣诵,F(xiàn)81;3狱窘,F(xiàn)84杜顺;4,HKY85蘸炸。當(dāng)設(shè)置usedata = 1時(shí)躬络,不能使用其它堿基替換模型。若需要選擇其它堿基替換模型搭儒,則考慮使用usedata = 3參數(shù)穷当,就可以設(shè)置model參數(shù)值為其它的堿基替換模型。此時(shí)淹禾,程序會(huì)將數(shù)據(jù)進(jìn)行分割馁菜,例如,ndata = 3稀拐,則程序分別對(duì)3個(gè)數(shù)據(jù)進(jìn)行分析火邓,生成baseml的輸入文件,并調(diào)用baseml進(jìn)行分析(若數(shù)據(jù)量較大,分析較慢铲咨,推薦按ctrl + c 強(qiáng)行終止躲胳,程序則終止baseml的運(yùn)行,進(jìn)行下一個(gè)數(shù)據(jù)的輸入文件生成并運(yùn)行下個(gè)數(shù)據(jù)的baseml命令纤勒,再按 ctrl + c 強(qiáng)行終止坯苹,這樣可以讓程序生成3各數(shù)據(jù)的輸入文件和baseml配置文件,然后手動(dòng)并行化運(yùn)行摇天,加快時(shí)間)粹湃;每個(gè)分析結(jié)果中會(huì)得到rst2結(jié)果文件;可以在此處手動(dòng)修改tmp0001.ctl配置文件泉坐,例如为鳄,修改其model參數(shù)值,推薦修改outfile參數(shù)值腕让,從而可以手動(dòng)并行化運(yùn)行baseml命令孤钦;然后將所有數(shù)據(jù)的rst2結(jié)果使用cat命令合并成文件in.BV,用于下一步設(shè)置參數(shù)usedata = 2纯丸,進(jìn)行MCMC分析偏形。
       alpha = 0.5       *核酸序列中不同位點(diǎn),其進(jìn)化速率不一致觉鼻,其變異速率服從GAMMA分布俊扭。一般設(shè)置GAMMA分布的alpha值為0.5。若該參數(shù)值設(shè)置為0坠陈,則表示所有位點(diǎn)的進(jìn)化速率一致萨惑。此外,當(dāng)userdata = 2時(shí)畅姊,alpha咒钟、ncatG吹由、alpha_gamma若未、kappa_gamma這4個(gè)GAMMA參數(shù)無效。因?yàn)閡serdata = 2時(shí)倾鲫,不會(huì)利用到多序列比對(duì)的數(shù)據(jù)粗合。
       ncatG = 5         *設(shè)置離散型GAMMA分布的categories值。
     BDparas = 1 1 0.1   *設(shè)置出生率乌昔、死亡率和取樣比例隙疚。若輸入有根樹文件中的時(shí)間單位發(fā)生改變,則需要相應(yīng)修改出生率和死亡率的值磕道。例如供屉,時(shí)間單位由100Myr變換為1Myr,則要設(shè)置成".01 .01 0.1"。
 kappa_gamma = 6 2       *設(shè)置kappa(轉(zhuǎn)換/顛換比率)的GAMMA分布參數(shù)伶丐。
 alpha_gamma = 1 1       *設(shè)置GAMMA形狀參數(shù)alpha的GAMMA分布參數(shù)悼做。 

***進(jìn)化速率參數(shù):**
 rgene_gamma = 2 20 1    *設(shè)置序列中所所有位點(diǎn)平均[堿基/密碼子/氨基酸]替換率的Dirichlet-GAMMA分布參數(shù):alpha=2、beta=20哗魂、初始平均替換率為每100 million 年(取決于輸入有根樹文件中的時(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)化速率沒有差異峻呕,即方差為0利职,該參數(shù)無效;當(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ù)莫秆。

***MCMC參數(shù):**
       print = 1         *設(shè)置打印mcmc的取樣信息:0间雀,不打印mcmc結(jié)果;1镊屎,打印除了分支進(jìn)化速率的其它信息(各內(nèi)部節(jié)點(diǎn)分歧時(shí)間惹挟、平均進(jìn)化速率、sigma2值)缝驳;2连锯,打印所有信息归苍。 
      burnin = 2000      *將前2000次迭代burnin后,再進(jìn)行取樣(即打印出該次迭代計(jì)算的結(jié)果信息运怖,各內(nèi)部節(jié)點(diǎn)分歧時(shí)間霜医、平均進(jìn)化速率、sigma2值和各分支進(jìn)化速率等)驳规。
    sampfreq = 10        *每10次迭代則取樣一次肴敛。
     nsample = 20000     *當(dāng)取樣次數(shù)達(dá)到該次數(shù)時(shí),則取樣結(jié)束(程序也將運(yùn)行結(jié)束)吗购。至少要一萬個(gè)

二.MAGA

  • CLOCKS -> Compute Timetree -> RelTime-ML
  • LOAD SEQUENCE DATA : alignment data (fasta format)
  • LOAD TREE FILE : tree data (nwk format,no branch)
  • Select Taxa : add outgroup data(click OK)
  • Select Branch : (click a branch then click OK)
  • CALIBRATE NOTE -> (click clock icon ) -> SELECT CALIBRATE TYPE(Normal Distribution) -> (click OK)
    可以多線程運(yùn)行趴俘,但是由于電腦的內(nèi)存不足可能無法運(yùn)行(8條序列蚤吹,1101462 * 8個(gè)堿基,54GB內(nèi)存)

也可以下載運(yùn)行l(wèi)inux版本的maga。只不過下載的都是預(yù)編譯的版本令野,需要系統(tǒng)內(nèi)部特定版本的glibc 纳击,為此還得專門去編譯 glibc俏讹。

三.BEAST

BEAST1與BEAST2比較 - 高芳鑾的博文

  • 將多序列比對(duì)文件(fasta文件)轉(zhuǎn)換成nexus文件(可以使用mega的另存為
  • 使用BEAST下面的子程序 beauti (可視化頁面)將nex文件轉(zhuǎn)換成xml文件蝴罪。
# 導(dǎo)入beast 所需要的庫文件,否則不能正常運(yùn)行
git clone https://github.com/beagle-dev/beagle-lib.git
cd beagle-lib/
mkdir build
cd build/
cmake -DCMAKE_INSTALL_PREFIX=/home/bqxiao/program/beagle ..
make
make install
export LD_LIBRARY_PATH=/home/bqxiao/program/beagle/lib:$LD_LIBRARY_PATH

3.1 BEAST 1

BEAST 運(yùn)行的本身命令比較簡單:

nohup beast file.xml -threads N &
# beast(linux版) 本身是一個(gè)shell腳本埠偿,通過調(diào)用java來運(yùn)行beast.jar程序透罢,在beast腳本內(nèi)部,可以設(shè)置給java虛擬機(jī)的內(nèi)存大小冠蒋。當(dāng)選擇核算替換模型和線程數(shù)時(shí)羽圃,應(yīng)當(dāng)考慮堆內(nèi)存大小。

beauti 主要參數(shù):
patitions : 在這個(gè)頁面顯示導(dǎo)入的數(shù)據(jù)(xml抖剿,nex)的分區(qū)信息朽寞。在這里可以設(shè)置分區(qū)之間是否使用相同的模型,先驗(yàn)樹等信息
Taxa : 添加校準(zhǔn)節(jié)點(diǎn)斩郎,可以在序列數(shù)據(jù)中設(shè)置分類單元子集(相同子集拖入 Included Taxa 內(nèi))脑融,在后續(xù)可以添加校準(zhǔn)時(shí)間
Tip : 為每一條序列添加注釋(如:采集時(shí)間,采集地點(diǎn)等)
Sites : 進(jìn)化模型選擇
Clocks : 分子鐘模型設(shè)置
Tree : 種群或者物種規(guī)模隨時(shí)間的變化模型缩宜,主要考慮序列來自物種內(nèi)還是物種間肘迎。
Priors : 為前面的所有模型添加先驗(yàn)信息,在這里可以設(shè)置前面的 Taxa 校準(zhǔn)時(shí)間脓恕。
MCMC

  • Length of chain :MCMC分析運(yùn)行的總代數(shù)膜宋。可以通過以增加代數(shù)來提高ESS值(Effective Sample Size)使參數(shù)收斂炼幔。增加采樣頻率也可以提高ESS值,但是過于采樣頻繁也不會(huì)提高ESS值史简。
  • Log parameters every : 采樣頻率乃秀,即多少代生成一個(gè)樣本肛著。這個(gè)值的設(shè)置受Length of chain影響,最主要的一點(diǎn)需要保證最后的樣本最少得有10,000個(gè)跺讯。
  • Echo state to screen : 將狀態(tài)摘要輸出到BEAST控制臺(tái)窗口的頻率枢贿。
    最后點(diǎn)擊生成xml文件

3.2 BEAST 2

BEAST 2中beast的使用方法和BEAST 1 中的beast使用方法一致。而子程序 beauti 有些許偏差刀脏。

四.r8s

r8s 分析速度非尘旨裕快,對(duì)于一棵樹只需要幾秒鐘可以得出結(jié)果愈污。其中 pyr8s 是 r8s 的重寫版本耀态。其中 r8s 中的很多參數(shù)的(模型)算法只有部分在 pyr8s 得到的。

4.1 r8s v1.8.1

原本的r8s目前已經(jīng)不支持了暂雹,僅有 r8s download | SourceForge.net 這個(gè)網(wǎng)站還有源碼下載首装,但是依然有一些問題。
比如在編譯 fortran 一些文件時(shí)杭跪,有些代碼過于古老仙逻,特性已被刪除,不能通過編譯涧尿。
因此可以將make報(bào)錯(cuò)部分手動(dòng)編譯

# 原始make部分
gfortran -fopenmp -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /home/bqxiao/miniconda3/include -c -o tn.o tn.f
# 手動(dòng)編譯系奉,保留已刪除的特性: -std=legacy
# 其他標(biāo)準(zhǔn)還有 -std=f95 -std=f2003 -std=f2008
gfortran -std=legacy -fopenmp -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /home/bqxiao/miniconda3/include -c -o tn.o tn.f

另一部分編譯報(bào)錯(cuò)于全局變量重定義 -> 只需要在重定義部分前置 extern 修飾即可。
使用 :

r8s -b -f file.nex > file.out
# -b 使用命令行而不是交互式

配置文件如下:

#NEXUS

begin trees;
tree PAUP_9 = (Marchantia:0.033817,(Lycopodium:0.040281,((Equisetum:0.048533,(Osmunda:0.033640,Asplenium:0.036526):0.000425):0.011806,((((Cycas:0.009460,Zamia:0.018847):0.005021,Ginkgo:0.014702):1.687e-86,((Pinus:0.021500,(Podocarpac:0.015649,Taxus:0.021081):0.006473):0.002448,(Ephedra:0.029965,(Welwitsch:0.011298,Gnetum:0.014165):0.006883):0.016663):0.006309):0.010855,((Nymphaea:0.016835,(((((Saururus:0.019902,Chloranth:0.020151):1.687e-86,((Araceae:0.020003,(Palmae:0.006005,Oryza:0.031555):0.002933):0.007654,Acorus:0.038488):0.007844):1.777e-83,(Calycanth:0.013524,Lauraceae:0.035902):0.004656):1.687e-86,((Magnolia:0.015119,Drimys:0.010172):0.005117,(Ranunculus:0.029027,((Nelumbo:0.006180,Platanus:0.002347):0.003958,(Buxaceae:0.013294,((Pisum:0.035675,(Fagus:0.009848,Carya:0.008236):0.001459):0.001994,(Ericaceae:0.019136,Solanaceae:0.041396):0.002619):1.687e-86):0.004803):1.687e-86):0.006457):0.002918):0.007348,Austrobail:0.019265):1.687e-86):1.687e-86,Amborella:0.019263):0.003527):0.021625):0.012469):0.019372);
End;

begin rates;

blformat nsites=952 lengths=persite;

collapse;

mrca LP marchantia pisum;
mrca ANGIO amborella pisum;

mrca VP lycopodium pisum;
mrca V1 Equisetum pisum;
mrca FERN Osmunda Asplenium;
mrca SP Ginkgo pisum;
mrca GYMNO Ginkgo Gnetum;
mrca CYC Cycas Zamia;
mrca G2 Taxus Podocarpac;
mrca GNET Ephedra Gnetum;

mrca CL Calycanth Lauraceae;
mrca ChS Chloranth Saururus;
mrca MON Acorus Oryza;
mrca PO Palmae Oryza;
mrca MAG Magnolia Drimys;
mrca EUD Ranunculus Carya;
mrca NEL Nelumbo Platanus;
mrca BX Buxaceae Carya;
mrca ERIC Ericaceae Solanaceae;
mrca FC Fagus Carya;

fixage taxon=LP age=450;

constrain taxon=CL min_age=110;
constrain taxon=ChS min_age=125;
constrain taxon=MON min_age=120;
constrain taxon=PO  min_age=85;
constrain taxon=MAG min_age=125;
constrain taxon=NEL min_age=110;
constrain taxon=BX  min_age=105;
constrain taxon=ERIC min_age=91;
constrain taxon=FC  min_age=96;

constrain taxon=VP  min_age=420;
constrain taxon=V1  min_age=390;
constrain taxon=FERN  min_age=255;
constrain taxon=GYMNO  min_age=310;
constrain taxon=CYC  min_age=220;
constrain taxon=GNET  min_age=115;
constrain taxon=G2  min_age=210;

constrain taxon=EUD max_age=125;
constrain taxon=SP max_age=330;

unfixage taxon=lp;

set num_time_guesses=5;
unfixage taxon=platanus;
constrain taxon=platanus maxage=90;

divtime method=LF algorithm=tn;
showage;
describe plot=chronogram;
describe plot=tree_description;

end;

4.2 pyr8s

安裝非常簡單

pip install git+https://github.com/iTaxoTools/pyr8s.git

4.2.1 python 內(nèi)部使用

pyr8s 原文檔 使用方法有問題

import pyr8s.parse as p
a = p.parse('tests/legacy_1')
# 其中 parse.py 中根本沒有 parse 函數(shù)

通過查看 run.py 可以發(fā)現(xiàn) 姑廉,應(yīng)當(dāng)是以下代碼:

import pyr8s.parse as p
a=p.from_file("tests/legacy_1")
result=a.run()

4.2.2 shell 環(huán)境中使用

pyr8s file.nex > result.txt

配置文件如下:

#NEXUS


begin trees;
tree Apiaceae = [&U] (((Angelica_sinensis:0.0263287636,((Apium_graveolens:0.0460894249,Heracleum_sosnowskyi:0.0322190162):0.0024624336,(Coriandrum_sativum:0.0283896098,Peucedanum_praeruptorum_Dunn:0.0307326278):0.0032486520):0.0020207251):0.0171794647,Oenanthe_sinensis:0.0537308346):0.0110179388,Daucus_carota:0.0638858678,cacao:0.6170394091);
End;


begin rates;

[ 
    nsites : 核苷酸的個(gè)數(shù)
    persite : 最大似然樹
]
blformat nsites=1101462 lengths=persite;

collapse;


mrca DO Daucus_carota Oenanthe_sinensis;
mrca CH Coriandrum_sativum Heracleum_sosnowskyi;


constrain taxon=DO min_age=24 max_age=44;
constrain taxon=CH min_age=5 max_age=40;

set num_time_guesses=5;

divtime method=NP algorithm=PL;

[
    showage : 展示各個(gè)分支的節(jié)點(diǎn)的年齡
]
showage;

[
    describe plot=chronogram : 在終端畫出樹的形狀
]
describe plot=chronogram;

[
    describe plot=tree_description : 在終端導(dǎo)出這棵樹的nwk格式
]
describe plot=tree_description;

end;
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
  • 序言:七十年代末喜最,一起剝皮案震驚了整個(gè)濱河市,隨后出現(xiàn)的幾起案子庄蹋,更是在濱河造成了極大的恐慌瞬内,老刑警劉巖,帶你破解...
    沈念sama閱讀 218,941評(píng)論 6 508
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件限书,死亡現(xiàn)場(chǎng)離奇詭異虫蝶,居然都是意外死亡,警方通過查閱死者的電腦和手機(jī)倦西,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 93,397評(píng)論 3 395
  • 文/潘曉璐 我一進(jìn)店門能真,熙熙樓的掌柜王于貴愁眉苦臉地迎上來,“玉大人扰柠,你說我怎么就攤上這事粉铐。” “怎么了卤档?”我有些...
    開封第一講書人閱讀 165,345評(píng)論 0 356
  • 文/不壞的土叔 我叫張陵蝙泼,是天一觀的道長。 經(jīng)常有香客問我劝枣,道長汤踏,這世上最難降的妖魔是什么织鲸? 我笑而不...
    開封第一講書人閱讀 58,851評(píng)論 1 295
  • 正文 為了忘掉前任,我火速辦了婚禮溪胶,結(jié)果婚禮上搂擦,老公的妹妹穿的比我還像新娘。我一直安慰自己哗脖,他們只是感情好瀑踢,可當(dāng)我...
    茶點(diǎn)故事閱讀 67,868評(píng)論 6 392
  • 文/花漫 我一把揭開白布。 她就那樣靜靜地躺著才避,像睡著了一般橱夭。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上工扎,一...
    開封第一講書人閱讀 51,688評(píng)論 1 305
  • 那天徘钥,我揣著相機(jī)與錄音,去河邊找鬼肢娘。 笑死呈础,一個(gè)胖子當(dāng)著我的面吹牛,可吹牛的內(nèi)容都是我干的橱健。 我是一名探鬼主播而钞,決...
    沈念sama閱讀 40,414評(píng)論 3 418
  • 文/蒼蘭香墨 我猛地睜開眼,長吁一口氣:“原來是場(chǎng)噩夢(mèng)啊……” “哼拘荡!你這毒婦竟也來了臼节?” 一聲冷哼從身側(cè)響起,我...
    開封第一講書人閱讀 39,319評(píng)論 0 276
  • 序言:老撾萬榮一對(duì)情侶失蹤珊皿,失蹤者是張志新(化名)和其女友劉穎网缝,沒想到半個(gè)月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體蟋定,經(jīng)...
    沈念sama閱讀 45,775評(píng)論 1 315
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡粉臊,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 37,945評(píng)論 3 336
  • 正文 我和宋清朗相戀三年,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了驶兜。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片扼仲。...
    茶點(diǎn)故事閱讀 40,096評(píng)論 1 350
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡,死狀恐怖抄淑,靈堂內(nèi)的尸體忽然破棺而出屠凶,到底是詐尸還是另有隱情,我是刑警寧澤肆资,帶...
    沈念sama閱讀 35,789評(píng)論 5 346
  • 正文 年R本政府宣布矗愧,位于F島的核電站,受9級(jí)特大地震影響迅耘,放射性物質(zhì)發(fā)生泄漏贱枣。R本人自食惡果不足惜监署,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 41,437評(píng)論 3 331
  • 文/蒙蒙 一颤专、第九天 我趴在偏房一處隱蔽的房頂上張望纽哥。 院中可真熱鬧,春花似錦栖秕、人聲如沸春塌。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,993評(píng)論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽只壳。三九已至,卻和暖如春暑塑,著一層夾襖步出監(jiān)牢的瞬間吼句,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 33,107評(píng)論 1 271
  • 我被黑心中介騙來泰國打工事格, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留惕艳,地道東北人。 一個(gè)月前我還...
    沈念sama閱讀 48,308評(píng)論 3 372
  • 正文 我出身青樓驹愚,卻偏偏與公主長得像远搪,于是被迫代替她去往敵國和親。 傳聞我的和親對(duì)象是個(gè)殘疾皇子逢捺,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 45,037評(píng)論 2 355

推薦閱讀更多精彩內(nèi)容