劉清雪， 王 亮

(吉林建筑科技學院 計算機科學與工程學院， 吉林 長春 130111)

0 引 言

系統(tǒng)發(fā)生是指生物形成或進化的歷史[1]。系統(tǒng)發(fā)生分析是在所有可以系統(tǒng)發(fā)生過程中尋找最能真實反映歷史的發(fā)生樹，研究物種之間的進化關系。系統(tǒng)發(fā)生樹的構建是一個NP問題，現(xiàn)在最常用的構建進化樹的方法有距離法、最大簡約法和最大似然法。最大簡約法是基于最優(yōu)原則的建樹方法之一，其前期基礎是“多序列比對”，研究結果以樹的形式進行描述，并根據(jù)有無共同的祖先將其分為有根樹和無根樹。

1 基礎理論與算法分析

系統(tǒng)發(fā)生樹也稱進化樹，一般以二叉樹的形式表示，它是由一系列代表當前物種的葉子節(jié)點和代表推斷的祖先或進化事件發(fā)生的位置內(nèi)部節(jié)點組成。根據(jù)系統(tǒng)發(fā)生樹能否推斷出共同祖先和進化方向分為有根樹和無根樹。樹的分支式樣無論有根還是無根均被稱為拓撲結構。系統(tǒng)發(fā)生樹可能的個數(shù)隨序列個數(shù)急劇增加。由N個物種建立系統(tǒng)發(fā)生樹，則可能的有根樹(NR)和無根樹(NU)的數(shù)目為：

(1)

(2)

1.1 后綴法

二叉樹的遍歷是二叉樹最基本的操作，一次完整的遍歷可以產(chǎn)生樹中所有結點的一個線性序列。設L,V,R分別表示在位于某個結點時移到左子女、訪問結點數(shù)據(jù)和移到右子女的操作，按從左到右的順序遍歷，則存在三種遍歷方式，即LVR、VLR、LRV。根據(jù)V與L和R的相對位置，分別稱這三種為中序、前序和后序遍歷。如果用二叉樹表示表達式，則對其進行中序、前序和后序遍歷產(chǎn)生的線性序列與表達式的中綴、前綴和后綴自然對應。

1.2 簡約法

簡約法的理論基礎是解釋一個過程，最好的理論是所需假設數(shù)目最少的一個。最大簡約法通過分析待比較序列的信息位點，確定相應的系統(tǒng)發(fā)生樹，該樹能用最小的動作產(chǎn)生序列的差異。樹的代價計算就是沿著各個分支累加最小核苷酸替換總數(shù)。具體實現(xiàn)算法為當內(nèi)部節(jié)點的兩個直接后代節(jié)點的核苷酸相交為非空時，該節(jié)點最可能的候選核苷酸集就是這個交集；否則為它的兩個后代節(jié)點上核苷酸的并集。當一個并集成為一個節(jié)點的核苷酸集時，通向該節(jié)點分支的某個位置上必定發(fā)生一個核苷酸替換，在統(tǒng)計過程中，只考慮信息位點的替換即可。一般來說，最大簡約法適用于以下情況：待檢序列的堿基數(shù)目大，但差異小；堿基變異率近似相等；沒有過多的顛換、轉換傾向[2-3]。

1.3 遺傳算法[4-5]

遺傳算法是模擬生命進化的隨機搜索算法，是一種解決復雜優(yōu)化問題的通用框架。該算法首先按規(guī)則產(chǎn)生一定規(guī)模的初始群體，然后使其中的個體按一定的概率進行選擇、交叉與變異，并按照特定的適應度評價函數(shù)進行判斷，選擇復制優(yōu)秀個體，組成新的一代，從多個可行解開始，按一定的遺傳方式進行迭代，并趨于最優(yōu)解，但其具有隨機漫游性、早熟性等特點，會使局部最優(yōu)。

1.4 模擬退火算法[6-9]

模擬退火算法模擬固體退火過程，將固體加熱至高溫時，其內(nèi)能增大，內(nèi)部粒子變?yōu)闊o序狀，在其慢慢冷卻直至常溫的過程中，內(nèi)能逐漸減少，并在基態(tài)時達到最小。由初始解i和初始溫度t0開始，對當前解重復“產(chǎn)生新解→計算內(nèi)能模擬函數(shù)差(fi-fj)→接受或舍棄”的迭代，并逐步衰減溫度t值，算法終止時，當前解即為所得近似最優(yōu)解，其以一定的概率P接受惡化解的特性將有效規(guī)避遺傳算法的早熟性。

(3)

2 改進算法----基于后綴法編碼的遺傳退火簡約法

2.1 編碼方式、適應度函數(shù)和種群初始化[10]

2.1.1 編碼方式

對于有n個物種的有根二叉樹，其中內(nèi)部結點用“+”表示，葉子結點用數(shù)字n表示，連續(xù)的2個數(shù)字和1個“+”表示一個最小的分枝。

2.1.2 適應度函數(shù)

以簡約法中整個樹的信息位點突變數(shù)量為適應度函數(shù)，樹長為適應值，并以歷史最大簡約樹為最終的最大簡約樹。

2.1.3 種群初始化

根據(jù)輸入的n個物種及n-1個內(nèi)部結點，隨機排列產(chǎn)生后綴編碼序列。

2.2 遺傳退火算子設計

根據(jù)后綴法的編碼方式及特點，設計了選擇退火操作、復制操作、葉子交換退火操作、葉子與分枝交換退火操作、分枝與分枝交換退火操作。

2.2.1 選擇退火操作

采用隨機遍歷抽樣法隨機選擇兩個初始個體Pi和Pj，計算其對應二叉樹的樹長fi和fj,如fi≤fj，則選擇個體Pi遺傳到下一代；否則，以一定的概率P接受Pj，以此保存物種的多樣性。

2.2.2 復制操作

將種群中序列按后序遍歷的方法轉換成對應的二進制樹，再通過適應度函數(shù)計算樹長，并根據(jù)適應度值進行排序。

2.2.3 交換退火操作[11]

對群體中的個體，以一定的概率將指定兩位點的序列進行交換。同時要確保經(jīng)過變異后產(chǎn)生的解依然滿足二叉樹的定義，不會產(chǎn)生不可行解，因此設計以下三種交換子。

2.2.3.1 葉子交換退火操作

在親代的后綴編碼序列中隨機選擇兩個葉子位置，以2，7為例進行交換，8個葉子二叉樹及其后綴編碼如圖1所示。

圖1 8個葉子二叉樹及其后綴編碼

交換后形成新子代樹結構及編碼序列如圖2所示。

圖2 葉子交換后二叉樹及其后綴編碼

此種交換是單位點的交換，樹的拓撲結構沒有變化，進化速度較慢，概率值較高。計算交換前后兩個體適應度值fi和fj,適應度值低的絕對保留，但也以一定的概率P接受惡化的適應度值高的解。

2.2.3.2 葉子與分枝交換退火操作

在親代的樹型結構中隨機選擇一個葉子位置和一個分枝位置，以2和(7，8)為例，將其交換。交換產(chǎn)生的子樹結構及編碼序列如圖3所示。

圖3 葉子與分枝交換后二叉樹

圖4 分枝與分枝交換后二叉樹

此種交換后樹的拓撲結構發(fā)生變化，進化速度最快，為了控制結果的魯棒性，概率值應較低。計算交換前后兩個體適應度值fi和fj,選優(yōu)遺傳到下一代的同時，同樣以一定的概率P接受惡化解。

2.3 更新初始群體

將新產(chǎn)生的序列和原有序列共同作為下一次迭代的初始群體。并將當前群的最佳個體與歷史最佳個體進行比較，以低者作為歷史最佳個體。

2.4 結束條件

群體迭代次數(shù)達到預定值或進化N代以后,最大簡約樹的樹長不再變化，這時，便可將最后一代進化中最佳個體作為最大簡約樹。

2.5 算法流程

算法流程如圖5所示。

3 模擬實驗與結果分析

文中序列資料來源于https://treebase.org/treebase-web/home.html[12]的TreeBASE數(shù)據(jù)源。以算法的空間、時長和最大簡約樹的樹長為衡量標準,Phylip是一個綜合的系統(tǒng)發(fā)生分析軟件包，其中Dnapars是采用最大簡約法對原始數(shù)據(jù)和重采樣數(shù)據(jù)進行系統(tǒng)發(fā)生樹構建，文中將改進算法與Dnapars軟件相比較。實驗環(huán)境：CPU,intel(R)Core(TM) i3-2100@3.10 GHz,RAM,8 GB,WIN32位操作系統(tǒng)。鑒于各算子對早熟現(xiàn)象的影響程度分析，經(jīng)過反復調(diào)整，算法參數(shù)如下：群體規(guī)模為500，最大迭代次數(shù)為500次，選擇算子為隨機選擇，葉子交換k1=0.003，葉子與分枝交換k2=0.002，分枝與分枝交換k3=0.001[13-14]。

圖5 改進算法流程

用簡約法構建系統(tǒng)發(fā)生樹時，運行時間和樹長一直是兩個非常重要的指標，且隨著物種的增多和序列長度的增加，這兩個指標顯示尤為重要。樹長越小，樹越準確，時間越小，算法更優(yōu)。兩種算法實驗結果對比見表1。

表1 兩種算法實驗結果對比表

從表1中數(shù)據(jù)明顯可以得出，改進算法在這兩個指標上都有明顯提升，且隨物種和序列長度的增加，這種提升越明顯。因為在編碼過程中，后綴法編碼采用數(shù)組的方式進行存貯，使遺傳操作在線性時間內(nèi)完成，退火算子的引入又保持了物種的多樣性，跳出了局部最優(yōu)解，本算法能在時間和空間兩個維度將進化樹的構建算法得以改進。