張麗娜，榮昌鶴，何 遠(yuǎn)，關(guān) 瓊，何 彬，朱興文，劉佳妮，陳紅菊3,,*

1. 大理學(xué)院 數(shù)學(xué)與計算機學(xué)院，云南 大理 671003

2. 云南林業(yè)職業(yè)技術(shù)學(xué)院，云南 昆明 650224

3. 紅河學(xué)院，云南 蒙自 661100

4. 中國科學(xué)院昆明動物研究所 遺傳資源與進(jìn)化國家重點實驗室，計算生物學(xué)與醫(yī)學(xué)生態(tài)學(xué)研究組，云南 昆明 650223

系統(tǒng)發(fā)育就是指生物譜系的分支演化歷史，或是指生命自起源后的整個遺傳進(jìn)化史 （Avise，2006），系統(tǒng)發(fā)育樹是描述物種間或操作分類單元間（operation taxonomic units，OTUs）系統(tǒng)發(fā)育關(guān)系的圖論模型。操作分類單元可以是現(xiàn)存物種、基因、基因組或者是任何其他可操作單元。系統(tǒng)發(fā)育樹的構(gòu)建就是從現(xiàn)存物種和古生物學(xué)記錄存留的證據(jù)來重現(xiàn)生命進(jìn)化史的科學(xué)探索。用偉大的進(jìn)化生物學(xué)家 Dobzhansky（1973）的名言“如果沒有進(jìn)化論，生物學(xué)的一切便毫無意義”來強調(diào)系統(tǒng)發(fā)育樹的重要性是恰如其分的。

由于技術(shù)限制，最初分類學(xué)家只能依靠生物的形態(tài)特征來推斷物種間的親緣關(guān)系。但表型特征存在一定的局限性，由于趨同進(jìn)化現(xiàn)象，有時候親緣關(guān)系很遠(yuǎn)的生物體也表現(xiàn)出很大的相似性，如鯨和蝙蝠，雖然形態(tài)差異很大，但都具有發(fā)達(dá)的高頻回聲定位能力。同時，許多生物個體可能由于體型較小，數(shù)量多而導(dǎo)致對其表型特征的研究較困難，如各類微生物。另外許多生物體間的共同特征少之又少，很難發(fā)現(xiàn)何種表型特征能用來研究比較。隨著分子生物學(xué)研究的不斷發(fā)展和檢測核苷酸序列和各種氨基酸序列技術(shù)的成熟，使得從小分子層面上構(gòu)建系統(tǒng)發(fā)育樹成為可能。近年來測序技術(shù)的迅猛發(fā)展，使得測序成本降低，涌現(xiàn)的海量核酸序列、氨基酸序列數(shù)據(jù)也被收集于如 GenBank、EMBL和DDBJ 等大型數(shù)據(jù)庫中，促使人們可從更大范圍上建立物種間的遺傳進(jìn)化關(guān)系。分子水平的進(jìn)化研究具有傳統(tǒng)方法不可比擬的優(yōu)勢，可從核酸和氨基酸序列差異程度來精確判斷物種進(jìn)化的時期和速度，確定親緣關(guān)系極遠(yuǎn)的生物體間的進(jìn)化關(guān)系，同時能對體型較小的微生物間的進(jìn)化關(guān)系進(jìn)行深入研究。

目前許多系統(tǒng)發(fā)育樹構(gòu)建算法都是從解決最優(yōu)化問題出發(fā)，如最大簡約法、最大似然法等，但是這些方法受物種數(shù)量嚴(yán)格限制，當(dāng)物種數(shù)量較多時，構(gòu)建系統(tǒng)發(fā)育樹是一個典型的 NP-complete 難題 （Foulds & Graham，1982）。這意味著在多項式時間內(nèi)不能被計算機求解，只能被非確定機求解；不能得到絕對數(shù)值解，只能通過比較相對解來確定最合適的答案。然而慶幸的是人們后來發(fā)明了改進(jìn)算法：啟發(fā)式搜索算法，通過分割數(shù)據(jù)集 （操作單元）變成小的子集，再對小的子集使用最優(yōu)化算法 （最大似然或最大簡約算法等）求出每個子集對應(yīng)的最優(yōu)樹，然后合并每個子集得到的最優(yōu)樹，最終形成整個數(shù)據(jù)集的最優(yōu)樹。

隨著生物信息學(xué)的發(fā)展，使用計算機技術(shù)處理系統(tǒng)發(fā)育樹成為不可或缺的理論，構(gòu)建系統(tǒng)發(fā)育樹的軟件包的相繼出現(xiàn)，并得到了廣泛的應(yīng)用。對構(gòu)建進(jìn)化樹程序包的算法、運用限制條件及其優(yōu)缺點的了解，有助于我們選用合適的建樹方法和分析軟件，更進(jìn)一步說，為我們對現(xiàn)有方法的改進(jìn)和編寫性能更完善的軟件提供思想源泉和幫助。

1 構(gòu)建系統(tǒng)發(fā)育樹的一般過程

不同的領(lǐng)域?qū)溆胁煌亩x，下面簡單列舉了部分樹的定義及生物信息中與系統(tǒng)發(fā)育樹相關(guān)的基本術(shù)語。

樹（圖論中定義）：連通的無環(huán)圖稱為樹。度為 1 的叫葉子節(jié)點，度大于等于 1 的為根節(jié)點，節(jié)點間的連線叫樹枝。

樹（數(shù)據(jù)結(jié)構(gòu)中定義）：由一個集合以及在該集合上定義的一種非線性結(jié)構(gòu)關(guān)系。

樹（生物信息中定義）：表示物種之間的進(jìn)化關(guān)系的樹狀圖譜。由樹枝和節(jié)點組成。節(jié)點分為內(nèi)部節(jié)點和外部節(jié)點，內(nèi)部節(jié)點代表的是進(jìn)化事件發(fā)生的位置或進(jìn)化過程中的共同祖先，外部節(jié)點又叫葉子節(jié)點，代表的是不同物種或是可操作單元。樹枝是連接各節(jié)點的邊，樹枝長度代表的是生物進(jìn)化時間或進(jìn)化距離。葉子節(jié)點的度為 1，內(nèi)部節(jié)點的度至少為 3。如圖 1a 所示，節(jié)點 A?D 為葉子節(jié)點，節(jié)點 1、2 為內(nèi)部節(jié)點，節(jié)點 0 為根節(jié)點。根據(jù)拓?fù)浣Y(jié)構(gòu)的不同系統(tǒng)發(fā)育樹可以分為有根樹和無根樹。有根樹（圖 1a）有一個根節(jié)點，代表所有其他節(jié)點的共同祖先，從根節(jié)點只有唯一路徑經(jīng)進(jìn)化到達(dá)其他任何節(jié)點；無根樹（圖 1b）只表明了節(jié)點之間的關(guān)系，沒有進(jìn)化方向，但是通過引入外群或外部參考物種可以在無根樹中指派根節(jié)點（Gregory，2008）。

圖1 系統(tǒng)發(fā)育樹Figure 1 Phylogenetic treea：有根樹；b：無根樹。 a: Rooted tree; b: Unrooted tree.

構(gòu)建系統(tǒng)發(fā)育樹包括選擇同源序列、序列比對、計算推斷進(jìn)化樹、評估進(jìn)化樹四個步驟。具體流程如圖 2 所示。

圖2 構(gòu)建系統(tǒng)發(fā)育樹流程圖Figure 2 Phylogenetic tree flowchart

構(gòu)建系統(tǒng)發(fā)育樹的第一步是選擇同源序列作為計算數(shù)據(jù)。這一步實際上包含兩個過程：一是收集序列數(shù)據(jù)，二是確定數(shù)據(jù)的同源性。序列數(shù)據(jù)可以通過實驗或通過公共數(shù)據(jù)庫下載獲得。目前公共數(shù)據(jù)庫主要有美國的 GenBank、歐洲的 EMBL、日本的 DDBJ 等。

序列比對提供一種衡量核酸或蛋白質(zhì)序列之間相關(guān)性的度量方法。將兩條或多條序列寫成兩行或多行，使盡可能多的相同字符出現(xiàn)在同一列中，將不同序列中的每一位點進(jìn)行逐一比對，構(gòu)建一個打分矩陣來表示序列間的相似性或同源性。DNA序列在進(jìn)化中由于替換、插入/刪除、突變事件使其發(fā)生改變，所以在比對中，錯配與突變相應(yīng)，而空位與插入或缺失對應(yīng)。最常用的比對工具有Blast（Altschul et al，1990）、Clustal（Larkin et al，2007）、 Muscle（Edgar，2004）和FASTA（Lipman & Pearson，1985）等。

計算推斷系統(tǒng)發(fā)育樹的主要任務(wù)是求出最優(yōu)樹的拓?fù)浣Y(jié)構(gòu)和估計分支長度。這部分算法及常用軟件在后面詳細(xì)介紹。

評估的目的是對已經(jīng)得出的系統(tǒng)發(fā)育樹的置信度進(jìn)行評估，常用的方法有自舉檢驗法 （bootstrap methods）（Felsenstein，1985；Penny & Hendy，1985）及刀切法 （jackknife methods）（Shao & Tu，1996）。自舉檢驗法是從原始序列中隨機選取堿基組成和原始序列相同長度的新序列，這樣在每個序列中有些堿基被重復(fù)選擇，而有些堿基未被選擇，按這樣的方法取出和原始數(shù)據(jù)序列數(shù)相同的新序列組成新的組。將所有的新序列組用某種算法生成多個新的進(jìn)化樹。將生成的許多進(jìn)化樹進(jìn)行比較，把所有新的樹中相同拓?fù)浣Y(jié)構(gòu)最多的樹認(rèn)為是最真實的樹，樹中分支位置的數(shù)值表示該種結(jié)構(gòu)占所有樹中的百分比值，該值小于 75通常都認(rèn)為是置信度較低的分支。刀切法是對原始數(shù)據(jù)進(jìn)行“不放回式”隨機抽取，從數(shù)據(jù)集里去除一部分序列數(shù)據(jù)或每次去掉一個分類群對象，然后對剩下的數(shù)據(jù)進(jìn)行系統(tǒng)發(fā)育分析。刀切法產(chǎn)生的數(shù)據(jù)小于原始數(shù)據(jù)，（delete-half-jackknifing）（Felsenstein，1985；Wu，1986）。兩類檢測方法的差別在于，前者是對全部數(shù)據(jù)進(jìn)行“重置式”隨機抽取，數(shù)據(jù)抽到的概率是相等的，且建立的和原始數(shù)據(jù)大小相等，而后者是“不放回式”抽取，產(chǎn)生的數(shù)據(jù)小于原始數(shù)據(jù)。

2 構(gòu)建系統(tǒng)發(fā)育樹常用算法原理

基于分子水平的系統(tǒng)發(fā)育推斷方法可以分為兩大類，即基于離散特征的方法和基于距離的方法。基于離散特征的系統(tǒng)發(fā)育樹重構(gòu)算法通過搜索各種可能的樹，從中選出最能夠解釋物種之間進(jìn)化關(guān)系的系統(tǒng)發(fā)育關(guān)系樹，這類方法利用統(tǒng)計技術(shù)定義一個最優(yōu)化標(biāo)準(zhǔn)，對樹的優(yōu)劣進(jìn)行評價，包括最大簡約法（maximum parsimony methods）（Mount，2008）、最大似然法 （maximum likelihood methods）（Myung，2003）和貝葉斯法（Bayesian methods）（Holder & Lewis，2003）。距離法的理論基礎(chǔ)是最小進(jìn)化原理（minimum evolution，ME）（Saitou & Nei，1986），這類方法首先構(gòu)造一個距離矩陣來表示每兩個物種之間的進(jìn)化距離，然后基于這個距離矩陣，采用聚類算法對研究的物種進(jìn)行分類。通過不斷的合并距離最小的兩個節(jié)點和構(gòu)建新的距離矩陣，最終得出進(jìn)化樹。距離法包括非加權(quán)組平均（unweighted pair-group method with arithmetic mean，UPGMA）、鄰接法 （neighbor- joining，NJ）、距離變換法（transformed distance method）和鄰接關(guān)系法（neighbors relation method）等（Takezaki，1998）。非加權(quán)組平均法比較簡單，得出的系統(tǒng)發(fā)育樹不可加和，現(xiàn)在很少使用，常用鄰接法來構(gòu)建系統(tǒng)發(fā)育樹。表1列出了常用構(gòu)建系統(tǒng)發(fā)育樹的算法及支持軟件。

2.1 鄰接法

Kidd & Sgaramelh-Zonta（1971）最早提出基于距離數(shù)據(jù)的系統(tǒng)發(fā)育樹重構(gòu)算法，從所有可能的進(jìn)化樹中選擇進(jìn)化分支長度總和最小的那棵樹，距離法通常不能找到精確的最小進(jìn)化樹，只能找到近似的最小進(jìn)化樹，但是它的計算速度非?？?，而且準(zhǔn)確率較高，因此被廣泛應(yīng)用于系統(tǒng)發(fā)育分析（Zhang & Lai，2010）。當(dāng)可操作單元數(shù)量較多時，這種方法的計算量會大增，因此，又提出了啟發(fā)式搜索算法（Mucherino & Seref，2009）：從一個距離矩陣開始，采用一定的準(zhǔn)則，遞歸地合并矩陣中距離最短的節(jié)點，并重構(gòu)新的距離矩陣，直到只剩下最后一個分類單元為止。其中最常用的是鄰接法（Saitou & Nei，1986）。下面舉例說明鄰接法重建系統(tǒng)發(fā)育樹的過程。假設(shè)有以下 5 組同源序列：

S1: GTGCTGCACGGCTCAGTATAGCATTTA CCCTTCCATCTTCAGATCCTGAA

S2: ACGCTGCACGGCTCAGTGCGGTGCTTA CCCTCCCATCTTCAGATCCTGAA

S3: GTGCTGCACGGCTCGGCGCAGCATTTAC CCTCCCATCTTCAGATCCTATC

S4: GTATCACACGACTCAGCGCAGCATTTGC CCTCCCGTCTTCAGATCCTAAA

S5: GTATCACATAGCTCAGCGCAGCATTTG CCCTCCCGTCTTCAGATCTAAAA

表1 系統(tǒng)發(fā)育樹常用算法及支持軟件Table 1 Frequently-used algorithms and software for phylogeny reconstruction(http://evolution.genetics.washington.edu/ phylip/software.html)

以上5個序列中每個序列都含有50個堿基，每兩個序列之間的距離定義為失配堿基的個數(shù)（這里忽略刪除和插入事件）。則每次聚類可得出距離矩陣如表 2、3、4 所示。根據(jù)公式1

求出Q值。公式中 n 為物種個數(shù)或序列個數(shù)，在 n個序列組成的所有可能的無根樹中找出 Q 值最小的兩個序列組成鄰近關(guān)系，重新構(gòu)建距離矩陣，根據(jù)新的距離矩陣再找最小的 Q 值組成一組，反復(fù)上面的過程直到所有的序列都找到了自己的鄰居（Studier & Keppler，1988）。根據(jù)表2、3、4求出所有的 Q 值，見表5。

表 2 序列間距離矩陣Table 2 Pairwise distance matrix

由表 5 可推斷出 5 條序列的系統(tǒng)發(fā)育樹拓?fù)鋱D和各分支長度分別如圖3和圖4所示：

表3 第一次聚類得到的距離矩陣Table 3 Distance matrix after the first clustering

表 4 第二次聚類得到的距離矩陣Table 4 Distance matrix after the second clusterin

表 5 Studier J 和 Keppler K 方法得到的 Q 值表Table 5 Q value from Studier J and Keppler K

圖 3 NJ算法得到的系統(tǒng)發(fā)育樹拓?fù)鋱DFigure 3 Topology of Phylogenetic Tree with NJ Approach

圖4 估計各分支長度Figure 4 Branch-Length Estimation

隨后的研究在鄰接法基礎(chǔ)上又提出了很多改進(jìn)算法：Studier & Keppler（1988）提出的改進(jìn)算法，引入了線性數(shù)組的概念，大幅降低了計算的時間復(fù)雜度（Chen et al，2006）；Bruno et al（2000）提出了加權(quán)鄰接法（weighted neighbor-joining）算法、Gascuel（1997）提出了 BIONJ 算法、Desper &Gascuel（2012）提出的 FASTME 算法和Criscuolo& Gascuel（2008）提出了快速鄰接法算法，均縮短了建立系統(tǒng)發(fā)育樹的時間。距離法速度快，適合于大型數(shù)據(jù)集和自舉分析,允許不同序列間有不同的分支長度,允許多重替換，但當(dāng)序列差異很大時，轉(zhuǎn)換成距離矩陣會使序列信息減少，而且距離法只提供一棵可能的樹，并對模型的依賴比較強烈。

2.2 最大簡約法

最大簡約法是基于奧卡姆剃刀原則 （Occam's razor）而發(fā)展起來的一種進(jìn)化樹重構(gòu)的方法，即突變越少的進(jìn)化關(guān)系就越有可能是物種之間的真實的進(jìn)化關(guān)系，系統(tǒng)發(fā)生突變越少得到的系統(tǒng)發(fā)生結(jié)論就越可信（Sober，1988）。最大簡約法首先是由Camin & Sokal （1965）提出來的，經(jīng)過 H ein （1990，1993）的研究發(fā)展使得用最大簡約法來建立進(jìn)化樹得到極大的發(fā)展及應(yīng)用。

最大簡約法采用5個假設(shè) （Felsenstein，1978，1979，1981a,b）：（1）序列中的每個位點獨立進(jìn)化；（2）不同世系（lineage）獨立進(jìn)化；（3）序列上的位點（堿基或氨基酸）的替換概率小于該分枝系統(tǒng)發(fā)生時間的長度；（4）系統(tǒng)發(fā)生的不同分支改變有不同，但高變化率的分支和低變化率的分支間的變化大小不會相差很大；（5）位點間變化不會相差太大。一個位點的刪除和插入各算一個變化，當(dāng)然連續(xù)的刪除 N 個位點，應(yīng)該算作獨立的 N 個事件。

用簡約法推斷系統(tǒng)發(fā)生關(guān)系，首先判斷信息位點。信息位點是那些產(chǎn)生突變能把其中的一棵樹同其他樹區(qū)別開來的位點。如果一個位點是信息位點，那么該位點至少有兩種以上的核苷酸，并且每種核苷酸至少出現(xiàn)兩次（見表6）。簡約法中只考慮信息位點而不考慮非信息位點。

其次確定每棵樹的替換數(shù)目 （Fitch，1971）。這里以3棵樹為例來說明構(gòu)建過程，如圖5。要確定每棵樹的替換數(shù)目，就要從 5 個已知的外部節(jié)點上的核苷酸推斷出 4 個內(nèi)部節(jié)點上最可能的核苷酸。尋找內(nèi)部節(jié)點的算法如下：如果一個內(nèi)部節(jié)點的兩個直接后代節(jié)點上的核苷酸的交集為非空，g那么這個節(jié)點的最可能的候選核苷酸就是這個交集；否則為它的兩個后代節(jié)點上核苷酸的并集。當(dāng)一個并集成為一個節(jié)點的核苷酸集時，通向該節(jié)點的分支的某個位置必定發(fā)生一個核苷酸替換。故而并集中核苷酸的數(shù)目也是生成外部節(jié)點上的核苷酸的最小替換數(shù)，外部節(jié)點從它們的共同祖先出發(fā)，通過這些替換，形成當(dāng)前的核苷酸狀態(tài)。找好內(nèi)部節(jié)點后，即可計算該內(nèi)部節(jié)點后代的替換數(shù)。計算信息位點的替換數(shù)，是通過計算外部節(jié)點上不同核苷酸的數(shù)目減去 1 即可得到。考慮所有可能的樹，分別對每棵樹中的變化打分，統(tǒng)計每個位點的核苷酸最小替換數(shù)目，所有信息位點替換數(shù)的總和最小的樹即為最簡約樹。

表 6 4條同源序列的比對Table 6 4 Homology sequences alignment

圖 5 3 棵有根樹及內(nèi)部節(jié)點Figure 5 Three rooted trees and internal nodes

隨著序列數(shù)量的增加，可能的樹的拓?fù)浣Y(jié)構(gòu)呈現(xiàn)爆炸性增加 （如10 個物種，存在34 459 425 棵可能的無根樹[，遍歷這些可能的樹的拓?fù)浣Y(jié)構(gòu)，計算出最小替換數(shù)而找到最簡約樹，無疑計算量是相當(dāng)龐大的。對序列數(shù)據(jù)集較多的建樹，一般選用分支約束算法 （branchand-bound algorithm）（Land & Doig，1960）和啟發(fā)式算法 （heuristic algorithm）（Mucherino & Seref，2009）進(jìn)行樹的拓?fù)浣Y(jié)構(gòu)查找。

分支約束算法查找的樹，首先是從只有兩個物種組成的樹開始 （如果是無根樹，從3個物種的樹開始）；其次程序試著在合適的位置增加下一個物種，并對增加物種后的樹進(jìn)行替換數(shù)目的評價，迭代直到將所有的物種都加到樹上。它是一個深度優(yōu)先搜尋的過程 （depth-first search）（Even & Even，2011）。首先把第三個物種加在第一個可能的位置，這時第四個物種加在它的第一個可能的位置，再次是第五個物種，依次遍歷直到樹的第一個可能的樹產(chǎn)生。對樹的步數(shù)進(jìn)行衡量。改變物種的位置，直到遍歷所有的位置。四棵樹的深度優(yōu)先搜尋的過程如下：

首先建立兩個物種的樹：（A,B）

把C加到第一個可能的位置：（（A,B）,C）

把D加到第一個可能的位置：（（（A,D）,B）,C）

把D加到第二個可能的位置：（（A,（B,D））,C）

把D加到第三個可能的位置：（（（A,B）,D）,C）

把D加到第四個可能的位置：（（A,B）,（C,D））

把D加到第五個可能的位置：（（（A,B）,C）,D）

把C加到第二個可能的位置：（（A,C）,B）

把D加到第一個可能的位置：（（（A,D）,C）,B）

把D加到第二個可能的位置：（（A,（C,D））,B）

把D加到第三個可能的位置：（（（A,C）,D）,B）

把D加到第四個可能的位置：（（A,C）,（B,D））

把D加到第五個位置：（（（A,C）,B）,D）

把C加到第三個可能的位置：（A,（B,C））

把D加到第一個可能的位置：（（A,D）,（B,C））

把D加到第二個可能的位置：（A,（（B,D）,C））

把D加到第三個可能的位置：（A,（B,（C,D）））

把D加到第四個可能的位置：（A,（（B,C）,D））

把D加到第五個可能的位置：（（A,（B,C））,D）

如上所示，深度優(yōu)先搜尋也只不過是另外一種一次產(chǎn)生所有可能的樹的算法。即使物種數(shù)量中等，生成的可能樹的數(shù)量也是非常大的。當(dāng)然這種情況實際中是不會發(fā)生的，因為樹會以一個特定的順序生成，一些可能樹的拓?fù)浣Y(jié)構(gòu)是不會產(chǎn)生的。分支約束算法也是由這些深度優(yōu)先搜索步驟組成，只不過有一點改變，在樹的構(gòu)建過程中，部分樹如 （A,（B,C））的步數(shù)也被衡量。增加物種，預(yù)測會增加的步數(shù)，取增加步數(shù)的位置為增加的物種所在位置。分支約束算法會計算增加物種后不變的位點數(shù)量和變化的位點數(shù)量。因而如果 A、B 和C 及根有20個可變的位點，且如果樹（（A,C）,B）要求 24 步，當(dāng) D 增加有 8 個可變位點，那么，無論 D 加到哪個位置，最終的樹都不會少于 32 步。如果發(fā)現(xiàn)樹（（A,B）,（C,D））僅僅只有 30 步，那么我們就可以確定（（A,C）,B）上沒有位置可以讓 D 加上。分支約束算法會保留一個最簡約樹列表，這樣就可以砍掉一部分，從而避免一些可能的特定的樹的分支生成。因而分支約束算法能讓我們不必生成所有可能的樹而又能得到最簡約的樹，從而減少計算時間。

啟發(fā)式搜索算法通過子樹分支交換 （branch swapping），把分支嫁接到此步分析中找到的最好的那棵樹的其他位置，而產(chǎn)生一棵拓?fù)浣Y(jié)構(gòu)和初始樹相似的樹 （見圖6）。

圖 6 啟發(fā)式搜索剪除與嫁接Figure 6 Pruning and grafting of heuristic search

對于有7條序列的啟發(fā)式搜索在第一輪會產(chǎn)生上百棵新樹，計算突變數(shù)總和，其中比初始樹突變數(shù)更少的新樹被保留并在第二輪分析中被剪除和嫁接。重復(fù)這個過程，直到無法再產(chǎn)生比前一輪總突變數(shù)更少的樹，則此樹為最簡約樹。啟發(fā)式搜索能大大減少查找的可能樹的數(shù)量，從而解決對大量數(shù)據(jù)搜索樹的數(shù)量過大的問題。

最大簡約法可能會產(chǎn)生多棵簡約樹，此時通常選取一棵能概括這些簡約樹的一致樹（consensus tree）作為代表（Taylor et al，2011）。這種做法是將所有樹中都一致的分支點作為二叉分支點，不一致的分支點變?yōu)檫B接多個分支的內(nèi)部節(jié)點（如圖 7）。

圖 7 三個簡約樹對應(yīng)的一致樹Figure 7 Consensus tree form three MP trees

2.3 最大似然估計法

一般來講，如果模型合適，最大似然法的效果較好。最大似然法根據(jù)特定的 “替代模型”（substitution model）分析既定的一組序列數(shù)據(jù)，使所獲得的每一個拓?fù)浣Y(jié)構(gòu)的似然值最大。選出最大似然值最大的拓?fù)浣Y(jié)構(gòu)作為最優(yōu)系統(tǒng)樹。其分析的核心在于替代模型，常用的有 Jukes-Cantor 模型（Jukes & Cantor，1969），Kimura 雙參數(shù)模型（Kimura，1980）等。算法要求所有分類單元有完整的 DNA序列數(shù)據(jù) （如果有缺失則不計算），在運算過程中僅考慮堿基取代而忽略缺失/插入，算法相對費時。

在最大似然算法中，考慮拓?fù)浣Y(jié)構(gòu)和枝長兩個參數(shù)，并對似然率求最大值來估計枝長。算法基于統(tǒng)計特性，有很好的數(shù)學(xué)理論支持。在進(jìn)化速率可變的假設(shè)下，最大簡約法略差于轉(zhuǎn)換距離法和鄰接法的結(jié)果，最大似然法的結(jié)果最優(yōu) （Zhong et al，2001）。也就是說極大似然算法允許各分支進(jìn)化速率不同。極大似然算法原理如下：似然函數(shù)：給定進(jìn)化模型 M，模型的 K 個參數(shù)，進(jìn)化樹拓?fù)浣Y(jié)構(gòu)，枝長，當(dāng)前序列出現(xiàn)的可能性：如何取這些參數(shù)，使得該序列出現(xiàn)的可能性最大，

圖8 4個DNA序列Figure 8 Four DNA sequences

圖9 4個DNA序列可能的拓?fù)浣Y(jié)構(gòu)Figure 9 All possible trees come from four DNA sequences

因為有 3 個節(jié)點，每個節(jié)點可能的值是ATGC，所以有 43=64 個通路。

L（第6列）= SUM L （所有可能的進(jìn)化路徑）= L（路徑1）+ L（路徑2）+ L（路徑3）+ … + L（路徑64）

圖 10中節(jié)點 1、2、3、4 為葉子節(jié)點，5、6 為內(nèi)部節(jié)點，0 為根節(jié)點，vi為枝長，是進(jìn)化樹的參數(shù)，參數(shù)的值由似然函數(shù)通過觀察到的序列來估計。節(jié)點 K 的似然函數(shù)：

其中g(shù)x0表示節(jié)點 0 為核苷酸x0時的先驗概率，常常等于核苷酸在整個序列中的相對頻率，它可以用 ML 法來估計。Pij（v）為給定位點在時間0時的核苷酸 i 到時間 t 變?yōu)楹塑账?j 的概率，i、j 指 A、T、G、C 的任一種，在ML算法中允許各分支的替代速率 r 不同，用 vi=riti來表示第 i 個分支的預(yù)期替代數(shù)。計算 Pij（v）需要使用特定的替換模型。Felsenstein（Felsenstein，1981a）使用了等輸入模型。在此模型中 Pii（v）和 Pij（v）為：

當(dāng) gi=1/4，v=4rt 時，上述模型演變?yōu)?Jukes-Cantor模型。針對不同類型的數(shù)據(jù)選擇合適的模型可以增加準(zhǔn)確度。以上過程分析了有根樹的算法，如果使用一個可逆模型，即不論向前還是向后進(jìn)化核苷酸的替代過程不變。用數(shù)學(xué)表述為：

這樣節(jié)點 5 和 6 之間的核苷酸替代數(shù)（v5+v6）恒定而與根節(jié)點 0 的位置無關(guān)。計算 Lk時，指定圖 10 的 v5+v6為 v5，并假設(shè)進(jìn)化開始于該樹的某一點，為方便起見，假定從節(jié)點 5 開始，大大簡化了樹的復(fù)雜度，具體如圖 11 所示。這樣1 就可以簡化為：（5）

圖10 TTAG可能的進(jìn)化通路圖Figure 10 The evolutionary pathway of TTAG

圖11 有根樹轉(zhuǎn)為無根樹Figure 11 Rooted tree into a unrooted tree

到此我們只考慮了一個核苷酸位點，在整個建樹過程中我們必須考慮包括不變位點在內(nèi)的所有核苷酸位點。整個序列的似然率 L 是對所有位點的 Lk求積，整個樹的似然率對數(shù)為：

通過改變參數(shù) Vi，使 lnL 最大化，計算方法可以使用 Newton 方法或其他數(shù)值計算方法實現(xiàn)。最后選出似然值最大的拓?fù)浣Y(jié)構(gòu)作為最優(yōu)系統(tǒng)樹。

2.4 貝葉斯算法

基于統(tǒng)計學(xué)規(guī)律運作的算法還有貝葉斯算法，與極大似然估計算法不同的是，后者指定樹的結(jié)構(gòu)和進(jìn)化模型，計算序列組成的概率，從而推斷出對應(yīng)的進(jìn)化樹。前者正好相反，是由給定的序列組成，計算進(jìn)化樹和進(jìn)化模型的概率。

其中，P(T,θ) 為給定的樹T和參數(shù)θ的先驗概率/邊緣概率，是不考慮序列時的概率。 P (T,θ|D)為給定序列下的后驗概率， P (D |T,θ)為給定的樹T和參數(shù)θ的似然值，分母P(D)是一正則化常數(shù)。該定理表明后驗信息可由前驗信息和堿

基序列信息所得（Yang & Rannala，2012）。具體原理如圖 12 所示。

圖12 貝葉斯算法進(jìn)化樹原理圖Figure 12 Schematic of phylogenetic tree fromBayesian algorithm

開始不知道樹的概率，先假設(shè)每棵樹的可能性都是相等的，將 DNA 序列信息和進(jìn)化模型代入貝葉斯公式計算每棵樹的可能性，取概率最大者為最后的進(jìn)化樹。圖12的拓?fù)渲校╔,（Y,W））的進(jìn)化樹概率最大，所以為最后的進(jìn)化樹。每個系統(tǒng)樹的拓?fù)浣Y(jié)構(gòu)分布在不同區(qū)間；每棵樹的位置受到拓?fù)浣Y(jié)構(gòu)及枝長的影響 （Sanmartín et al，2008）。對系統(tǒng)發(fā)生問題，難以得到各概率的解析解，現(xiàn)有的解決辦法主要是 MCMC （Markov chain monte carlo sampling）方法。將進(jìn)化樹（拓?fù)浣Y(jié)構(gòu)與進(jìn)化模型參數(shù)）轉(zhuǎn)換為馬爾科夫鏈，待馬爾科夫鏈?zhǔn)諗坑诤篁灨怕史植技纯伞?/p>

2.5 系統(tǒng)發(fā)育樹重建常用的軟件包介紹

目前有很多軟件包可以進(jìn)行系統(tǒng)發(fā)生樹推斷及可靠性檢驗，還有像Unifrac和ITOL（interactive tree of life）等在線畫樹和分析樹的工具。網(wǎng)站http://evolution.genetics.washington.edu/phylip/softw are.html 列出了150多種相關(guān)軟件包，并可以對軟件進(jìn)行按類別查詢，如按軟件的運行系統(tǒng)、使用的算法等進(jìn)行查詢，對軟件進(jìn)行簡單介紹同時提供了下載的鏈接。具體使用時可按需求用不同的軟件，這里簡單介紹3種最常用的軟件。

2.5.1 PHYLIP

PHYLIP （phylogeny inference package）是由美國華盛頓大學(xué) Felsenstein 用 C 語言編寫的系統(tǒng)發(fā)生推斷軟件包,它提供免費的源代碼，支持Windows 和 Linux 等多種系統(tǒng)。在 3.69 版本中，由 35 個子程序組成，可以實現(xiàn)最大似然法、最大簡約法和距離法建樹。最大似然法有兩類程序：帶生物鐘的建樹子程序（dnamlk、promlk），可對進(jìn)化似然距離進(jìn)行估計；不帶生物鐘建樹程序 （dnaml、proml）。最大簡約法也有帶分子鐘建樹子程序（dnapennys），可以對進(jìn)化距離進(jìn)行估計；和不帶生物鐘的建樹子程序（dnapars、protpars）。距離法建樹由 dnadist、prodist、fitch、kitsch、neighbor 等子程序組成，dnadist 和 prodist 可實現(xiàn) F84、Kimura、Jukes-Cantor、LogDet 模型計算距離矩陣，fitch 子程序可實現(xiàn)不帶分子鐘的 Fitch-Margoliash法畫樹，而 neighbor 子程序帶有鄰接法和非加權(quán)組平均法兩種畫樹方法。每種建樹方法都帶有各自許多不同的選項供研究人員根據(jù)自己研究的目的進(jìn)行選擇優(yōu)化。軟件包帶有畫樹的子程序：可以畫三角形有根樹及矩形有根樹（drawgram），也可以畫無根樹（drawtree）。子程序 seqboot 使用自舉檢驗法或刀切法對構(gòu)建的樹進(jìn)行標(biāo)準(zhǔn)誤估計及可靠性檢驗，提供分析報告。此程序包還可以實現(xiàn)一致樹的構(gòu)建（consensus），以及樹的重構(gòu)（retree）等等。唯一不方便的是該程序包基于命令行形式，操作界面不夠友好。

2.5.2 MEGA

MEGA（molecular evolutionary genetics analysis）是由美國賓夕法尼亞州立大學(xué)Masatoshi Nei 等編寫的進(jìn)行分子進(jìn)化遺傳分析的軟件包。目前最新的版本為 5.0。它能對核酸序列及氨基酸序列進(jìn)行系統(tǒng)發(fā)生分析。在建樹方法上,提供了距離法中的非加權(quán)組平均和鄰接法及 MP 法，5.0 版本還提供了最大似然法算法，對構(gòu)建的樹可進(jìn)行自舉檢驗及標(biāo)準(zhǔn)誤估計的可靠性檢驗，并提供分析報告。該軟件不僅可以對本地序列文件進(jìn)行分析，而且可Web在線搜索分析，可以分析 NCBI 數(shù)據(jù)庫中的序列文件來重建進(jìn)化樹。該軟件可畫出矩形、三角形、圓形等多種形狀的系統(tǒng)發(fā)育樹。

2.5.3 MrBayes

MrBayes （Bayesian inference of phylogeny）是由John Huelsenbeck 等編寫，使用馬爾可夫鏈方法來估計參數(shù)模型的后驗概率分布。該軟件采用命令行形式，支持 Windows 和 UNIX 等多種系統(tǒng)，能夠處理核苷酸、氨基酸、限制性酶切位點和形態(tài)數(shù)據(jù)等多種數(shù)據(jù)，同時集成了多物種溯祖算法，支持正向、負(fù)向和總線形拓?fù)浣Y(jié)構(gòu)，支持 BEAGLE 數(shù)據(jù)庫，在使用兼容的硬件（NVIDIA圖形卡）條件下可以提高運行速度。表 7 列出了常用的建樹軟件及其特點。

當(dāng)序列間的分歧度不高，且序列足夠長時，鄰接法、最大簡約法和最大似然法得到的進(jìn)化樹往往具有相似的拓?fù)浣Y(jié)構(gòu) （Saitou & Imanishi，1989）。當(dāng)序列之間的分歧度比較高，將DNA序列轉(zhuǎn)為距離矩陣時往往會丟失一些信息 （Penny，1982）。而距離法的性能依賴于距離矩陣的質(zhì)量，因此，距離法只能當(dāng)序列滿足某些條件時才會有較高的準(zhǔn)確性。簡約法不依賴任何進(jìn)化模型，但進(jìn)化樹的簡約計分完全決定于重建祖先序列中的最小突變數(shù)，而突變是否按照事先約定的核苷酸最少替代途徑進(jìn)行是不得而知的。再者，所有分支的突變數(shù)不可能相同。由于沒有考慮核苷酸的突變過程，使得長分支末端的序列由于趨同進(jìn)化而顯示較好的相似性，導(dǎo)致對“長枝吸引” （Holder & Lewis，2003）的敏感。因此，當(dāng)序列分歧度較高時，最大簡約法極可能得出錯誤的拓?fù)浣Y(jié)構(gòu)。最大簡約法只適用于序列相似性較高的序列建立進(jìn)化樹，其次，最大簡約法在序列數(shù)據(jù)量較大的時候，建立進(jìn)化樹相當(dāng)耗時（是個 NP-complete 問題）（Foulds & Graham,1982）。最大似然法是一種建立在進(jìn)化模型上的統(tǒng)計方法，具有統(tǒng)計一致性、健壯性，能夠在一個統(tǒng)計框架內(nèi)比較不同的樹以及充分利用原始數(shù)據(jù)等優(yōu)點 （Bryant & Galtier，2005）。但它與鄰接法一樣需要選擇模型，一般選擇 Kimura-2 參數(shù)模型。但對于不同模型會得出不同的結(jié)果，算法相對比較耗時，適用于序列不是很多的情況。貝葉斯法因為后驗概率不僅涉及所有的樹， 而且對于每一棵樹還整合了枝長和替代模型參數(shù)值的所有可能組合，所以不可能采用常規(guī)的分析方法解決。所幸的是，一系列數(shù)值方法可用于近似地獲取后驗概率，其中最有用的就是馬爾可夫鏈·蒙特卡羅算法。其基本思想是建立馬爾可夫鏈，以替代模型參數(shù)作為狀態(tài)空間，其靜態(tài)分布就是參數(shù)的后驗概率分布。通過計算機模擬和抽樣技術(shù)獲得分支格局的后驗概率。同以往的最大似然法相比，貝葉斯推論的優(yōu)越性在于：能夠以很高的計算速度處理大型數(shù)據(jù)集，同時還使用后驗概率衡量樹的置信度 （Huelsenbeck &Ronquist，2001）。

表7 常用軟件及其特點Table 7 Frequently-used software and characteristic

3 總 結(jié)

近年人們在構(gòu)建系統(tǒng)發(fā)育樹方面已經(jīng)取得了很大進(jìn)展，構(gòu)建系統(tǒng)發(fā)育樹的算法和軟件也一直在不斷完善。通過對生物系統(tǒng)發(fā)生分析重建進(jìn)化樹，使人們更進(jìn)一步了解了生命進(jìn)化的歷史。但構(gòu)建系統(tǒng)發(fā)育樹是一個復(fù)雜的任務(wù)，目前還沒有哪一種算法能完全揭露真實的歷史進(jìn)化關(guān)系，同時現(xiàn)有算法存在的一些問題：如“長枝吸引”、非勻速分子鐘的修正、搜索算法的優(yōu)化等還有待進(jìn)一步解決。系統(tǒng)發(fā)育分析的算法和軟件的不斷完善將會在更多的研究領(lǐng)域得到應(yīng)用，如在基因多序列比對、分子鐘以及統(tǒng)計系統(tǒng)地理學(xué)分析的過程中（Yang &Rannala，2012）。