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

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

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

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

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

系統(tǒng)發(fā)育就是指生物譜系的分支演化歷史，或是指生命自起源后的整個(gè)遺傳進(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ù)來(lái)重現(xiàn)生命進(jìn)化史的科學(xué)探索。用偉大的進(jìn)化生物學(xué)家 Dobzhansky（1973）的名言“如果沒有進(jìn)化論，生物學(xué)的一切便毫無(wú)意義”來(lái)強(qiáng)調(diào)系統(tǒng)發(fā)育樹的重要性是恰如其分的。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

基于分子水平的系統(tǒng)發(fā)育推斷方法可以分為兩大類，即基于離散特征的方法和基于距離的方法。基于離散特征的系統(tǒng)發(fā)育樹重構(gòu)算法通過(guò)搜索各種可能的樹，從中選出最能夠解釋物種之間進(jìn)化關(guān)系的系統(tǒng)發(fā)育關(guān)系樹，這類方法利用統(tǒng)計(jì)技術(shù)定義一個(gè)最優(yōu)化標(biāo)準(zhǔn)，對(duì)樹的優(yōu)劣進(jìn)行評(píng)價(jià)，包括最大簡(jiǎ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)造一個(gè)距離矩陣來(lái)表示每?jī)蓚€(gè)物種之間的進(jìn)化距離，然后基于這個(gè)距離矩陣，采用聚類算法對(duì)研究的物種進(jìn)行分類。通過(guò)不斷的合并距離最小的兩個(gè)節(jié)點(diǎn)和構(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)組平均法比較簡(jiǎn)單，得出的系統(tǒng)發(fā)育樹不可加和，現(xiàn)在很少使用，常用鄰接法來(lái)構(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)化分支長(zhǎng)度總和最小的那棵樹，距離法通常不能找到精確的最小進(jìn)化樹，只能找到近似的最小進(jìn)化樹，但是它的計(jì)算速度非常快，而且準(zhǔn)確率較高，因此被廣泛應(yīng)用于系統(tǒng)發(fā)育分析（Zhang & Lai，2010）。當(dāng)可操作單元數(shù)量較多時(shí)，這種方法的計(jì)算量會(huì)大增，因此，又提出了啟發(fā)式搜索算法（Mucherino & Seref，2009）：從一個(gè)距離矩陣開始，采用一定的準(zhǔn)則，遞歸地合并矩陣中距離最短的節(jié)點(diǎn)，并重構(gòu)新的距離矩陣，直到只剩下最后一個(gè)分類單元為止。其中最常用的是鄰接法（Saitou & Nei，1986）。下面舉例說(shuō)明鄰接法重建系統(tǒng)發(fā)育樹的過(guò)程。假設(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個(gè)序列中每個(gè)序列都含有50個(gè)堿基，每?jī)蓚€(gè)序列之間的距離定義為失配堿基的個(gè)數(shù)（這里忽略刪除和插入事件）。則每次聚類可得出距離矩陣如表 2、3、4 所示。根據(jù)公式1

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

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

由表 5 可推斷出 5 條序列的系統(tǒng)發(fā)育樹拓?fù)鋱D和各分支長(zhǎng)度分別如圖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 估計(jì)各分支長(zhǎng)度Figure 4 Branch-Length Estimation

隨后的研究在鄰接法基礎(chǔ)上又提出了很多改進(jìn)算法：Studier & Keppler（1988）提出的改進(jìn)算法，引入了線性數(shù)組的概念，大幅降低了計(jì)算的時(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í)間。距離法速度快，適合于大型數(shù)據(jù)集和自舉分析,允許不同序列間有不同的分支長(zhǎng)度,允許多重替換，但當(dāng)序列差異很大時(shí)，轉(zhuǎn)換成距離矩陣會(huì)使序列信息減少，而且距離法只提供一棵可能的樹，并對(duì)模型的依賴比較強(qiáng)烈。

2.2 最大簡(jiǎn)約法

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

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

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

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

表 6 4條同源序列的比對(duì)Table 6 4 Homology sequences alignment

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖 7 三個(gè)簡(jiǎn)約樹對(duì)應(yīng)的一致樹Figure 7 Consensus tree form three MP trees

2.3 最大似然估計(jì)法

一般來(lái)講，如果模型合適，最大似然法的效果較好。最大似然法根據(jù)特定的 “替代模型”（substitution model）分析既定的一組序列數(shù)據(jù)，使所獲得的每一個(gè)拓?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ù) （如果有缺失則不計(jì)算），在運(yùn)算過(guò)程中僅考慮堿基取代而忽略缺失/插入，算法相對(duì)費(fèi)時(shí)。

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

圖8 4個(gè)DNA序列Figure 8 Four DNA sequences

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

因?yàn)橛?3 個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)可能的值是ATGC，所以有 43=64 個(gè)通路。

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

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

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

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

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

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

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

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

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

2.4 貝葉斯算法

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

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

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

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

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

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

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

2.5.1 PHYLIP

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

2.5.2 MEGA

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

2.5.3 MrBayes

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

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

表7 常用軟件及其特點(diǎn)Table 7 Frequently-used software and characteristic

3 總 結(jié)

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