系統(tǒng)發(fā)育樹構建簡明教程

1. 前言

構建一棵系統(tǒng)發(fā)育樹是研究系統(tǒng)學和進化的基礎。然而，發(fā)育樹的構建是對分類單元進化歷史的推測，因此對發(fā)育樹的可靠性檢驗也是重中之重。一棵發(fā)育樹的獲得大致分為數據輸入、數據處理 、算法計算 、樹的獲取 、可靠性檢驗、樹的可視化 和樹的注釋等等。

數據輸入一般分為兩類，序列數據如DNA、AA、RNA，性狀數據如形態(tài)數據、各種生態(tài)學數據。數據處理包括序列比對、飽和檢驗、堿基替換模型比較等。算法計算涉及到各種構建發(fā)育樹的算法，如UPGMA、NJ、MP、ML、BI等，以及多基因數據的串聯(lián)方法和溯祖理論之間的比較。樹的獲取指得是如何從多棵樹中總結出一棵最優(yōu)樹。可靠性檢驗顧名思義，指對發(fā)育樹拓撲結構的檢驗，包括自展支持、后驗概率等，以及對不同發(fā)育樹之間的比較。樹的可視化是使用工具展示發(fā)育樹以及使其更加美觀的方法。樹的注釋包羅萬象，如分化時間校準、祖先序列推測和各種生物地理學的各種注釋，加深了樹的深刻程度，和各種具體問題聯(lián)系起來，使樹的思想成為進化研究的基礎方法。

Yang, Z. and Rannala, B., 2012. Molecular phylogenetics: principles and practice. Nature reviews. Genetics, 13(5), p.303.

2. 數據準備

2.1 序列比對

序列比對是系統(tǒng)發(fā)育樹構建的基礎，旨在找到理論上的同源位點。存在多種多樣的比對算法，可以由不同軟件來實現。常用的如Muscle、MAFFT等。此教程以Mega中的Muscle算法比對DNA非編碼序列為例。

Wong, K.M., Suchard, M.A. and Huelsenbeck, J.P., 2008. Alignment uncertainty and genomic analysis. Science, 319(5862), pp.473-476.

  1. 打開Mega， 調出Align模塊> edit/build Alignment > creat a new alignment > 選擇DNA，出現Alignment explorer。所有的比對編輯操作均在這里面實現。界面如下：


  2. 通過edit導入新序列，edit> edit sequences form file


此處的file需要特定的格式，如常見的fasta（.fas）格式。格式如下：

>seq_1_name
ATCT
>seq_2_name
ACCT
>seq_3_name
ACCG
......

fasta格式比較簡單，> 后跟隨序列名稱，不支持空格等大部分特殊字符，僅支持部分字符如_ . 等，序列另起一行。導入文件前，需要將測序得到的序列以此序列保存，擴展名 .txt， .fas均可。

  1. 序列導入以后即可比對，比對前界面如下：



    然后 Muscle > Align DNA (如為編碼序列，可截取完ORF以后，選擇Align Codons)，選擇默認設置，點擊compute即可。比對參數界面：



    比對后則如下：

注：當序列明顯無法比對整齊時，可檢查是否為反向互補序列。

比對后的序列兩段常常不整齊，此比對導致的gap（即-）會部分影響發(fā)育樹結果，兩段大部分截取刪除后，少部分gap可用‘‘?’’填充整齊。

  1. 按需保存和導出比對結果。

比對結果的保存優(yōu)先使用fasta格式，而一般發(fā)育樹構建軟件通常使用nexus格式作為輸入文件。Mega支持這兩種格式的輸出。


export of alignment.png

2.2 堿基替換模型選擇

堿基替換模型是對堿基突變的量化描述，一般為一個4*4的矩陣。

A G C T
A 0 a b c
G a 0 d e
C b d 0 f
T c e f 0

4個堿基的突變情況，共產生6個不同的速率，這是一個最一般的模型，稱為 General time reverse model，即GTR，該模型具有普適性，符合絕大部分數據集。該模型假設A -- T和T -- A具有相同的速率，即突變沒有時間方向性，是可逆的。此等假設并未考慮其是否符合生物學的本質，僅僅是計算可行性上的需要。

該一般模型的種種簡化，會產生各種模型，如F81、JC69、HKY等。一個重要的簡化是把6個速率分為兩類，即轉換和顛換兩種不同的速率也就是HKY模型。一般認為轉換比顛換更為容易。

指的一提的是，一條序列上的全部位點并不具有完全相同的替換矩陣，即不同位置的堿基有著不同的突變速率，對于這種異質的突變速率通常通過Gamma分布(G)來描述，Gamma分布是一個靈活而強大的分布，具有較好的計算性。此外可能仍有部分堿基很少有突變事件發(fā)生，這些堿基會共同由一個不變位點比例（I）來描述。

（此段存疑）一般認為替換模型估算的突變數量僅和枝長相關，不影響拓撲結構。不同模型對拓撲結構的影響并不大，然而在計算時間節(jié)點之類對枝長敏感的算法時模型的選擇和預設就更為重要。

2.2.1 PAUP + MrModelTest模型選擇法

  1. 輸入文件需要nexus (*.nex)格式，該格式可以由Mega比對完成后生成。首先用PAUP軟件讀取nexus格式文件。打開PAUP > 選擇nexus文件。


  2. 運行模型計算批處理文件。PAUP中 File > open > mrblock.dat 運行自動進行，生成mrmodel.scores文件，該文件通常和mrblock.dat,modelGUI.exe一同位于MrMTGUI文件夾下。



    Mrblock.dat會對24個模型分別計算其得分，得分記錄于mrmodel.scores中。

  3. 最優(yōu)模型評價。打開ModelGUI.exe > 點擊右側select file... > 選擇mrmodel.scores文件 > 點擊右側MrModeltest!!!。界面如下：



    不同模型間的比較評價遵循特定的標準，如hLRTs，AIC，AICc，BIC，BF等等。MrModeltest僅支持前兩種評價標準，通常選用AIC標準。結果重點解讀最優(yōu)模型及其參數，自動生成的ML，BI的命令模塊。
    最優(yōu)模型：

-----------------------------------------------------
*                                                  *
*        AKAIKE INFORMATION CRITERION (AIC)        *
*                                                  *
-----------------------------------------------------
 Model selected: HKY+G
   -lnL =   1489.4283
    K =     5
    AIC =   2988.8567
   Base frequencies: 
     freqA =     0.1933
     freqC =     0.3464
     freqG =     0.3042
     freqT =     0.1561
   Substitution model: 
    Ti/tv ratio =    2.5316
   Among-site rate variation
     Proportion of invariable sites = 0
     Gamma distribution shape parameter = 0.2838

包括各種參數如-ln*L*， AIC值，堿基頻率，堿基替換速率，Gamma   shape參數以及不變位點比例等。

貝葉斯命令模塊：

BEGIN MRBAYES;
    Lset  nst=2  rates=gamma;
    Prset statefreqpr=dirichlet(1,1,1,1);
END;

如構建貝葉斯發(fā)育樹，此段命令直接貼入nexus文件末尾即可。

ML命令模塊（此處僅針對Garli軟件）：

BEGIN PAUP;
    Lset  Base=(0.1933 0.3464 0.3042)  Nst=2  TRatio=2.5316  Rates=gamma  Shape=0.2838  Pinvar=0;
END;

同理，此命令也需貼入nexus文件末尾。

2.2.2 jModeltest法選擇模型

MrModeltest僅在24個模型中選擇，Modeltest也不過是48個，這可能并不能選擇到最優(yōu)模型。Jmodeltest支持多達88個甚至1624個模型比較。jModeltest具有友好的圖形界面，操作簡單，但極耗硬件資源，對于較大的數據集幾乎不能完成計算。簡要操作如下：

  1. 打開jModeltest，File > Load DNA alignment, 輸入數據，此處使用fasta格式。


  2. 計算模型得分。Analysis > compute likelihood scores > compute likelihood（選擇默認參數，如有需要可自行修改），等待計算完成。



  3. 模型評價。jModeltest支持多種方法評價模型，一般來說BIC優(yōu)于AIC，然而AIC更為常用。jModeltest 的結果可以保存到HTML格式使用瀏覽器查看。



    結果可以在result > show results table 中查看，紅色標記的即為最優(yōu)模型。


2.2.3 形態(tài)數據的模型選擇

除了較為常用的DNA序列的堿基突變模型選擇，氨基酸數據的突變模型復雜的多，此處并不涉及?？紤]到越來越多的形態(tài)數據的構樹，形態(tài)數據的模型選擇仍然需要考慮，然而以上兩種方法均不支持。此處介紹ModelFinder來評估其模型選擇。

待續(xù)...

3. 系統(tǒng)發(fā)育樹構建

系統(tǒng)發(fā)育樹構建方法通常分為兩類，基于距離的方法和基于性狀的方法。

基于距離的方法是系統(tǒng)學早期發(fā)展起來的，將序列轉化為距離矩陣然后根據距離矩陣構建聚類樹，優(yōu)點是速度極快，缺點很多，模型考慮簡單，不適合遠緣序列，不適合復雜序列，理論上不總是可以得到一個最優(yōu)樹。

基于性狀的方法是系統(tǒng)發(fā)育的主流，不轉化為距離矩陣，避免了數據的丟失，直接基于堿基序列計算。常見的包括最大簡約法、貝葉斯法和最大似然法。

最大簡約法不基于任何假設，不進行模型描述，認為具有最少突變步驟的發(fā)育樹是最優(yōu)樹，計算強度較小，缺點同樣是不適合遠緣序列，無法考慮到復雜的突變事件。

最大似然法的基礎是統(tǒng)計學的最大似然估計，把拓撲結構和枝長均視為參數，使觀測數據（即堿基序列）有最大的似然值的參數為最優(yōu)參數，即最優(yōu)樹。缺點是計算強度較大，可能會得到次優(yōu)樹。

貝葉斯方法則剛剛相反，基于觀測數據，得分最高的拓撲結構被認為是最優(yōu)樹。蒙特卡洛（MC）和馬爾科夫鏈（MC）的引入使得貝葉斯方法的得到極大的發(fā)展。貝葉斯方法具有較快的運算速度，多個鏈同時運行也可較大限度的避免局部最優(yōu)化，因此被認為是最好的發(fā)育樹構建方法。

3.1 最大簡約法

最大簡約法一般通過[PAUP]( http://paup.phylosolutions.com/ )來實現，操作較為簡單。

  1. 準備nexus文件，可通過Mega導出。nexus文件末尾需添加命令模塊(block)。

    begin paup;
log file=****.log;
set autoclose =yes;
set criterion = parsimony;
set root = outgroup;
set storebrlens = yes;
set increase=prompt;
outgroup *** ***;
bootstrap nreps=1000 search=heuristic conlevel=50 / addseq=random nreps=10 swap=tbr hold=1;
savetrees from=1 to=1 file=****.tre format= altnex brlens = yes
savebootp = NodeLabels MaxDecimals=0;
pscores /TL=yes CI=yes RI=yes RC=Yes HI=Yes;
end;

    此命令需要修改三處命令，1，指定log文件名，建議與nexus文件一致；2，指定外類群，多個外類群通過空格隔開；3， 指定tre文件名，同樣建議和nexus文件一致。

  2. PAUP直接打開保存好的nexus文件，就可直接運算。當bootstrap replicates達到預置的重復數時，程序結束。生成log文件和tre文件。


    paup compute.png

  3. 在此，我們涉及到兩個概念，nexus文件和外類群，我們對它稍做解釋。

nexus文件剛才選模型時已有所接觸，此處再次解釋。nexus和fasta文件一樣是系統(tǒng)發(fā)育處理中常見的格式，但遠比fasta復雜，可以記錄序列文件和樹文件。一個典型的nexus文件如下：

#NEXUS

begin taxa;
 dimensions ntax=12;
 taxlabels
     Lemur_catta 
     Homo_sapiens
     Pan
     Gorilla
     Pongo
     Hylobates
     Macaca_fuscata
     M._mulatta
     M._fascicularis 
     M._sylvanus
     Saimiri_sciureus
     Tarsius_syrichta
 ;
 end;

begin characters;
 dimensions nchar=898;
 format missing=? gap=- matchchar=. interleave=yes datatype=dna;
 options gapmode=missing;
 matrix

Lemur_catta      AAGCTTCATAGGAGCA
Homo_sapiens     AAGCTTCACCGGCGCA
Pan              AAGCTTCACCGGCGCA
Gorilla          AAGCTTCACCGGCGCA
Pongo            AAGCTTCACCGGCGCA
Hylobates        AAGCTTTACAGGTGCA
Macaca_fuscata   AAGCTTTTCCGGCGCA
M._mulatta       AAGCTTTTCTGGCGCA
M._fascicularis  AAGCTTCTCCGGCGCA
M._sylvanus      AAGCTTCTCCGGTGCA
Saimiri_sciureus AAGCTTCACCGGCGCA
Tarsius_syrichta AAGTTTCATTGGAGCC

Lemur_catta      CTGTCTAGCCAACTCT
Homo_sapiens     CTGCCTAGCAAACTCA
Pan              CTGCCTAGCAAACTCA
Gorilla          CTGCCTAGCAAACTCA
Pongo            CTGCCTAGCAAACTCA
Hylobates        CTGCCTTGCAAACTCA
Macaca_fuscata   CTGCCTAGCCAATTCA
M._mulatta       CTGCCTAGCCAATTCA
M._fascicularis  CTGCTTGGCCAATTCA
M._sylvanus      CTGCTTGGCCAACTCA
Saimiri_sciureus CTGCCTAGCAAACTCA
Tarsius_syrichta TTGCCTAGCAAATACA

begin assumptions;
 charset coding = 2-457 660-896;
 charset noncoding = 1 458-659 897-898;
 charset 1stpos = 2-457\3 660-896\3;
 charset 2ndpos = 3-457\3 661-896\3;
 charset 3rdpos = 4-457\3 662-.\3;
 
 exset coding = noncoding;
 exset noncoding = coding;
 
 usertype 2_1 = 4        [weights transversions 2 times transitions]
     a  c  g  t
 [a] .  2  1  2
 [c] 2  .  2  1
 [g] 1  2  .  2
 [t] 2  1  2  .
     ;
 usertype 3_1 = 4        [weights transversions 3 times transitions]
     a  c  g  t
 [a] .  3  1  3
 [c] 3  .  3  1
 [g] 1  3  .  3
 [t] 3  1  3  .
     ;
     
 taxset hominoids = Homo_sapiens Pan Gorilla Pongo Hylobates;
end;

begin paup;
 constraints ch = ((Homo_sapiens,Pan));
 constraints chg = ((Homo_sapiens,Pan,Gorilla));
end;

通常nexus包括文件說明頭：

#NEXUS

begin taxa;

矩陣維度說明：

dimensions ntax=12;
 taxlabels
     Lemur_catta 
     Homo_sapiens
     ......
 ;
 end;

begin characters;
 dimensions nchar=898;

字符串說明：(包括是否數據分段，此處為分段)

format missing=? gap=- matchchar=. interleave=yes datatype=dna;
options gapmode=missing;

序列矩陣：

Lemur_catta      AAGCTTCATAGGAGCA
Homo_sapiens     AAGCTTCACCGGCGCA
......

Lemur_catta      CTGTCTAGCCAACTCT
Homo_sapiens     CTGCCTAGCAAACTCA
......

除了以上基本模塊，還有一些其他的block塊。例如數據分段(按基因、按密碼子第幾位等)：

begin assumptions;
 charset coding = 2-457 660-896;
 charset noncoding = 1 458-659 897-898;
 charset 1stpos = 2-457\3 660-896\3;
 charset 2ndpos = 3-457\3 661-896\3;
 charset 3rdpos = 4-457\3 662-.\3;
 
 exset coding = noncoding;
 exset noncoding = coding;

算法命令模塊：(例如上文提到的PAUP的MP命令block，下面列出的替換模型模塊以及單系約束命令，當然這些block并非所有的程序均支持。)

usertype 2_1 = 4     [weights transversions 2 times transitions]
     a  c  g  t
 [a] .  2  1  2
 [c] 2  .  2  1
 [g] 1  2  .  2
 [t] 2  1  2  .
     ;
 usertype 3_1 = 4        [weights transversions 3 times transitions]
     a  c  g  t
 [a] .  3  1  3
 [c] 3  .  3  1
 [g] 1  3  .  3
 [t] 3  1  3  .
     ;
     
 taxset hominoids = Homo_sapiens Pan Gorilla Pongo Hylobates;
end;

begin paup;
 constraints ch = ((Homo_sapiens,Pan));
 constraints chg = ((Homo_sapiens,Pan,Gorilla));
end;

此外nexus也可以保存樹，如常見的.tre擴展名的文件實際上就是nexus的格式，例如：

#NEXUS 
Begin trees;  [Treefile saved Thu Oct 05 11:34:24 2017]
[!
>Data file = C:\Users\zz\Desktop\8.nex
>Bootstrap method with heuristic search:
>  Number of bootstrap replicates = 1000
>  Starting seed = 1945405146
>  Optimality criterion = parsimony
>    Character-status summary:
>      Of 646 total characters:
>        All characters are of type 'unord'
>        All characters have equal weight
>        547 characters are constant
>        69 variable characters are parsimony-uninformative
>        Number of parsimony-informative characters = 30
>    Gaps are treated as "missing"
>  Starting tree(s) obtained via stepwise addition
>  Addition sequence: random
>    Number of replicates = 10
>    Starting seed = 1945405146
>  Number of trees held at each step during stepwise addition = 1
>  Branch-swapping algorithm: tree-bisection-reconnection (TBR)
>  Steepest descent option not in effect
>  Initial 'MaxTrees' setting = 100
>  Branches collapsed (creating polytomies) if maximum branch length is zero
>  'MulTrees' option in effect
>  Topological constraints not enforced
>  Trees are unrooted
>
>   1000 bootstrap replicates completed
>   Time used = 0.75 sec
>   ]
>   tree PAUP_1 = [&U] (Ficus_maxima_ITS,((((2009198_Ficus_pumila_ITS_AF7,2010058_Ficus_sarmentosa_ITS_AH3)65,(2009274_Ficus_sarmentosa_ITS_AV9,GBOWS0029_Ficus_pubigera_ITS_AL1)62)89,2010066_Ficus_pumila_ITS_AH4)93,(2009415_Ficus_sagittata_ITS_AX5,2014059_Ficus_hederacea_ITS_AP7)99)67,2009411_Ficus_laevis_ITS_AX4)0;
>   End;

另外一個概念是外類群，指定外類群的目的是為了置根，置根有多種方法，如中點置根法、分子種賦根法，當然也還有一些較新的方法。置根之后的發(fā)育樹才有方向，才可以看到祖裔關系。這里我們關注的外類群的挑選規(guī)則：外類群應該是所有內類群的姐妹群，關系越近越好。親緣關系較遠的外類群容易與內類群形成長枝吸引。

Tria, F.D.K., Landan, G. and Dagan, T., 2017. Phylogenetic rooting using minimal ancestor deviation. Nature Ecology & Evolution, 1, pp.s41559-017.

3.2 最大似然法

最大似然法最初是為了解決簡約法的長枝吸引而引入的一種系統(tǒng)發(fā)育重建方法，其理論基礎是發(fā)展非常成熟的最大似然估計方法。一開始，在PhyML和Garli等軟件中實現最大似然法較為耗時，其后以RXaML、IQ-tree為代表的執(zhí)行快速自展的算法，極大地提高了運行速度，基本上是最快的系統(tǒng)發(fā)育方法。以RXaML-master為例介紹最大似然樹的構建：

  1. 運行RXaML-master。其沒有用戶界面，需要在cmd（命令提示符）運行?？截惓绦蛑罜盤，打開執(zhí)行文件（*.exe）所在文件夾。


  2. 拷貝比對好的fasta文件至該文件夾下。


  3. 打開cmd（命令提示符）程序。輸入cd fasta文件所在路徑，然后回車。cd此處意為change directory。

  1. 輸入程序參數，回車后即開始運行，結果在同一文件夾下。此處先示例，然后再介紹常用參數。




  2. 結果會生成多個.tre文件，選擇bipartion.tre即可在各種可視化工具中調整美化。
  3. 參數簡介如下 
    raxmlHPC 
-s 9.fas 
-m GTRGAMMAI 
-N autoMR 
-x 23155001098 
-p 23100551098 
-T 4 
-o Ficus_maxima_ITS 
-f a 
-n result
  • rxamlHPC 程序頭調用rxaml程序。
  • -s 輸入文件參數，支持phylip、fasta格式。
  • -m 模型參數，不同數據源支持不同的模型，常見為DNA的GTR系列模型，如GTRGAMMA，GTRCAT，GTRGAMMAI等等。不支持非GTR模型。具體參見軟件說明書。
  • -N 自展重復參數，通常使用1000，也可以讓軟件決定多少為合適，參數為antoMR。
  • -x, -p 隨機數種子，可隨意輸入。
  • -T 處理器線程數，僅在HPC版本及配置多線程環(huán)境后方才生效。
  • -o 外類群參數，多個外類群間用， 隔開。
  • -f 執(zhí)行分析參數，a 此處為系統(tǒng)發(fā)育常用的快速自展分析。
  • -n 輸出結果擴展名參數，可任意指定。

3.3 貝葉斯方法

Mrbayes是實現貝葉斯算法的主要軟件。

3.3.1 輸入文件準備

Mega所導出的nexus的格式和貝葉斯所支持的nex格式略有區(qū)別。Mega所導出nexus前文已列出，下面列出Mrbayes的nex文件頭：

begin data;
    dimensions ntax=12 nchar=898;
    format datatype=dna interleave=no gap=-;
    matrix

差異主要包括三點：

  • 修改begin taxa 為begin data

  • 修改datatype=nucleotide至data=dna （僅針對DNA數據）

  • 刪除taxalabels模塊

此外還包括文件末尾的Mrbayes模塊。該模塊靈活度很大，可簡單也可復雜，最簡單的情況僅為一個堿基突變模型指定模塊（模型選擇部分已提及），最復雜的情況可以包含從log文件、模型指定、運行參數指定、樹總結參數指定等等。

一個相對詳細的Mrbayes模塊如下：

# Mrbayes block 模塊，將此模塊按照序列數據要求修改后，附于nex文件末尾即可。
# []內為每行參數的注釋內容，可保留，軟件默認不讀取[]內的內容。

begin mrbayes;

    [writting to log file]
    log start filename=log.log; 

    [constraint some taxa as monophyly]
    constraint outgroups = 36 37 157 158;

    [designate outgroups as taxa index]
    outgroup 36;
    outgroup 37;
    outgroup 157;
    outgroup 158; 

    [designate data partitions as start and end points]
    CHARSET its = 1-721;
    CHARSET ets = 722-1218;
    CHARSET g3pdh = 1219-1985;
    partition gene=3: its,ets,g3pdh;
    set partition=gene;

    [designate model parameters as partitions]
    lset applyto=(1) nst=6 rates=invgamma;
    lset applyto=(2) nst=6 rates=gamma;
    lset applyto=(3) nst=6 rates=invgamma;
    Prset applyto=(all) statefreqpr=dirichlet(1,1,1,1); 

    [designate model parameter when data not partition]
    [lset nst=6 rates=invgamma;
    Prset statefreqpr=dirichlet(1,1,1,1);]

    [designate mcmc chain parameters] 
    mcmcp ngen=15000000 samplefreq=1000 nruns=4 checkpoint=yes checkfreq = 10000 stoprule=yes stopval=0.01;

    [execute mcmc]
    mcmc;
    
    [summary tree and postiror prabability after mcmc runs]
    sumt relburnin=yes burninfrac=0.25;
    sump relburnin=yes burninfrac=0.25;

end;

各命令解釋：

  • log start filename：指定運行過程中的log文件名；

  • constraint outgroups：約束單系類群，后接類群序號，可指定多個單系群，另起一行繼續(xù)寫即可，但當想約束多個外類群為單系時，軟件不生效；

  • outgroup：指定外類群，后接類群序號，可指定多個類群；

  • lset： 模型指定參數，包括若干子參數：

    • nst：堿基突變速率個數，在1，2，6之間選擇；

      • 1時，對應JC、F81等模型；

      • 2時，對應K80、HKY等模型；

      • 6時，對應TrN、K81、GTR等模型；

    • rates：突變速率模型等等；

      • +G時，為gamma

      • +I時，為pinv

      • +I+G時，為invgamma

  • Prset：突變速率先驗分布類型，Mrbayes中使用狄利克雷分布作為其先驗；

Dir(\vec{p}|\vec{α})= {Γ(\sum*{K=1}^Kα*k) \over Π*{k=1}^KΓ(α*k)} {\displaystyle \prod*{k=1}^K}p_k^{α*k-1}= {1 \over Δ(\vec {α})}\displaystyle \prod*{k=1}^Kp_k^{α*k-1}

  • mcmc：蒙特卡洛-馬爾科夫鏈參數，包括若干重要子參數：

    • ngen：鏈的代數；

    • samplefreq：取樣頻率，太高則最終取樣的樹太少，不足以估算后續(xù)的樹枝長和后驗概率參數，太小則計算壓力大且不容易收斂，一般1000-10000之間；

    • nruns：同時獨立運行的分析個數，默認為2，一般選擇4；

    • checkpoint：檢查點，設置為yes時軟件會定時記錄當前運行狀態(tài)，中斷恢復運算時會從此狀態(tài)恢復，Mrbayes 3.2之后此參數默認為yes

    • checkfreq：檢查頻率，默認2000，可設置為10000；

    • stoprule：終止規(guī)則，設置為yes時，收斂值小于設定值時會自動停止，盡管總代數可能還沒有達到設置的總代數（ngen）

    • stopval：設置的收斂值閾值，和stoprule=yes配合使用；

    • 此時設置各參數后會立刻運行，如果想僅設置參數而不立刻運行，可使用mcmcp參數代替mcmc；

  • sumt：summary tree的縮寫，總結樹，包括若干子參數:

    • burnin：老化值，mcmc鏈在收斂前存在一段變動較為明顯的過程，這個數據對發(fā)育樹枝長和后驗概率的估算影響很大，需要舍去，一般為總數的1/10-1/4。此處針對的時抽樣后樹的個數，如總代數1000000，抽樣頻率1000，安裝1/4的比例來去掉老化值的話，則應輸入250。為了簡便，直接指定老化比例，使用參數relburnin=0.25指定舍去前25%的代數，免去計算抽樣的過程。

    • contype：一致樹的類型，halfcompat為50%多數一致樹（默認），allcompat為嚴格一致二叉樹；

    • conformat：一致樹的格式，Figtree格式為默認，包含多個參數，simple較為簡單，僅包括枝長和后驗概率，為大多數軟件所支持。注意：如iTol等可視化軟件僅支持simple模式。

  • sump：summary probability的縮寫，估算后驗概率，參數同sumt。

  • CHARSET：對數據進行分段，每段可以單獨指定堿基突變模型，一個分段模塊如下：

CHARSET its = 1-721;
CHARSET ets = 722-1218;
CHARSET g3pdh = 1219-1985;
partition gene=3: its,ets,g3pdh;
set partition=gene;</pre>

配合數據分段的模型分段指定的模塊如下：

lset applyto=(1) nst=6 rates=invgamma;
lset applyto=(2) nst=6 rates=gamma;
lset applyto=(3) nst=6 rates=invgamma;
Prset applyto=(all) statefreqpr=dirichlet(1,1,1,1); </pre>

3.3.2 軟件操作

將準備好的nexus文件放入到Mrbayes程序文件夾下即可準備運行。

Mrbayes的運行操作分為兩種，一種打開Mrbayes后手動逐行輸入命令行，如外類群，鏈長等參數，此時不需要準備復雜的Mrbayes block。另一個即為前面的提到的準備一個詳盡的Mrbayes block，預先根據數據提前指定好其參數，此時重復分析時就不用每次重復手動輸入命令行。此外，第二種方法還有一個優(yōu)勢，即支持中斷續(xù)跑，如果在分析過程中，程序意外中斷，僅需在mcmcp 命令中加入 append=yes，如下：

[designate mcmc chain parameters] 
mcmcp ngen=15000000 samplefreq=1000 nruns=4 checkpoint=yes checkfreq = 10000 stoprule=yes stopval=0.01 append=yes;

Mrbayes免安裝，沒有圖形化界面，所有操作均通過輸入操作命令行進行，主界面如下：


如果選擇在nexus文件中已經準備好了詳盡的Mrbayes Block，此時就可以直接輸入exe name.nex，回車后即可等待運行結束，在當前文件夾下可見.tre格式的樹文件。運行過程界面如下：

3.3.3 Mrbayes的收斂問題

Mrbayes由于使用MCMC算法，與上述兩個系統(tǒng)發(fā)育軟件有所不同，會存在一個收斂問題。Mrbayes軟件使用Average standard deviation of split frequencies來判斷收斂，默認每5000代計算一次該參數，如果該參數小于0.01，軟件運行到預設代數則會自動停止，如果沒有就會提示是否需要繼續(xù)增加代數！

但是由于數據集的差別，即使增加到很大代數后，該參數仍然大于0.01，此時有以下幾種方法可以嘗試解決：

  1. 繼續(xù)提高代數，直至時間成本難以承受，或者代數增加到30000000-50000000代（約數，在該代數下，大部分數據集都應該完成收斂；

  2. Average standard deviation of split frequencies只是一種判斷收斂的方法，對于MCMC算法的輸出結果，ESS（ Effective Sample Size ）也是一個參見的判讀收斂的方式，該方式可以在Tracer中實現。當代數已經足夠大的時候，使用Tracer軟件打開.p文件即可，如果所有的ESS都大于200時，即可認為運算已經收斂；

  3. 提高抽樣頻率：可在1000-10000之間選擇，提高抽樣頻率可以加快收斂速度。同時最好加長運算鏈長，以保證最終有足夠多的樹用來總結樹的拓撲結構和后驗概率；

  4. 如果在Tracer中檢查仍然沒有收斂的趨勢，就需要對運行參數進行調節(jié)了?？梢哉{節(jié)的參數如下：

    • temp參數：熱鏈的溫度，也就是熱鏈冷鏈交換的頻率，默認為0.1，可以稍微調高一些；
    • Ngammacat：默認為4，可改為16，提交Gamma分布的離散化估計精細度，會明顯提高運算強度；

  5. 指定一顆起始樹：一般使用UPGMA樹作為起始樹，一顆合適的起始樹可以避免由不恰當起始樹造成的局部最優(yōu)和收斂緩慢。起始樹的格式如下：

    [begin trees;
tree usertree=(A:0.01,B:0.01,((C:0.01,D:0.01):0.01,E:0.01):0.01);
end;]

4. 發(fā)育樹的可視化

發(fā)育樹的可視化較為復雜，分為多個層面，從僅僅把發(fā)育樹本身展示出來，到對發(fā)育樹本身進行顏色標識強調等，再到加上各種各樣的注釋數據，如物種信息、分布地信息、形態(tài)特征以及基因結構多種多樣。

對于發(fā)育樹本身的展示和強調可以在Figtree、Mega等常用軟件中進行。

而對于添加注釋信息則iTol、treeio等工具較為常用。

4.1 Figtree

Figtree是一個圖形化的發(fā)育樹可視化軟件，使用較為簡單。界面如下，可自行研究：

4.2 iTol

iTol是一個發(fā)育樹可視化的在線工具，可以方便對發(fā)育樹進行注釋，可實現的效果如下：

具體教程待續(xù)……

5. 發(fā)育樹的檢驗

系統(tǒng)發(fā)育是一種歷史過程，任何基于分子數據集得到的發(fā)育樹都是對真實系統(tǒng)發(fā)生的推測。一個合理地假設，總是應該接收來自各種證據的檢驗。

系統(tǒng)發(fā)育假設檢驗（phylogenetic hypothesis testing）是用統(tǒng)計學方法檢驗兩個或多個不同發(fā)育樹的差異是否有統(tǒng)計學上的顯著性。系統(tǒng)發(fā)育檢驗需要數據集、模型、兩棵以上的發(fā)育樹。已有有大量的檢驗方法，主要包括頻率檢驗或者貝葉斯檢驗。一般來說，檢驗方法包括Approximately unbiased test，Approximate Bayesian posterior probability test，bootstrap probability test，Kishino-Hasegawa test，weighted Kishino-Hasegawa test，Shimodaira-Hasegawa test和weighted Shimodaira-Hasegawa test等。常用的為Approximately unbiased test (AU)和Kishino-Hasegawa test (KH)。

黃原. (2012). 分子系統(tǒng)發(fā)生學. 科學出版社. pp 381-393.

多個軟件都可以用于執(zhí)行這種檢驗，如 PAUP，TREE-PUZZLE等。此處，我們介紹consel 01j.

Shimodaira, H., & Hasegawa, M. (2001). CONSEL: for assessing the confidence of phylogenetic tree selection. Bioinformatics, 17(12), 1246-1247.

具體教程請點擊這里！

最后編輯于 ：
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成，瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明：文章內容（如有圖片或視頻亦包括在內）由作者上傳并發(fā)布，文章內容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關閱讀更多精彩內容

友情鏈接更多精彩內容