劉清雪, 王 亮
(吉林建筑科技學(xué)院 計算機科學(xué)與工程學(xué)院, 吉林 長春 130111)
系統(tǒng)發(fā)生是指生物形成或進化的歷史[1]。系統(tǒng)發(fā)生分析是在所有可以系統(tǒng)發(fā)生過程中尋找最能真實反映歷史的發(fā)生樹,研究物種之間的進化關(guān)系。系統(tǒng)發(fā)生樹的構(gòu)建是一個NP問題,現(xiàn)在最常用的構(gòu)建進化樹的方法有距離法、最大簡約法和最大似然法。最大簡約法是基于最優(yōu)原則的建樹方法之一,其前期基礎(chǔ)是“多序列比對”,研究結(jié)果以樹的形式進行描述,并根據(jù)有無共同的祖先將其分為有根樹和無根樹。
系統(tǒng)發(fā)生樹也稱進化樹,一般以二叉樹的形式表示,它是由一系列代表當(dāng)前物種的葉子節(jié)點和代表推斷的祖先或進化事件發(fā)生的位置內(nèi)部節(jié)點組成。根據(jù)系統(tǒng)發(fā)生樹能否推斷出共同祖先和進化方向分為有根樹和無根樹。樹的分支式樣無論有根還是無根均被稱為拓撲結(jié)構(gòu)。系統(tǒng)發(fā)生樹可能的個數(shù)隨序列個數(shù)急劇增加。由N個物種建立系統(tǒng)發(fā)生樹,則可能的有根樹(NR)和無根樹(NU)的數(shù)目為:
(1)
(2)
二叉樹的遍歷是二叉樹最基本的操作,一次完整的遍歷可以產(chǎn)生樹中所有結(jié)點的一個線性序列。設(shè)L,V,R分別表示在位于某個結(jié)點時移到左子女、訪問結(jié)點數(shù)據(jù)和移到右子女的操作,按從左到右的順序遍歷,則存在三種遍歷方式,即LVR、VLR、LRV。根據(jù)V與L和R的相對位置,分別稱這三種為中序、前序和后序遍歷。如果用二叉樹表示表達式,則對其進行中序、前序和后序遍歷產(chǎn)生的線性序列與表達式的中綴、前綴和后綴自然對應(yīng)。
簡約法的理論基礎(chǔ)是解釋一個過程,最好的理論是所需假設(shè)數(shù)目最少的一個。最大簡約法通過分析待比較序列的信息位點,確定相應(yīng)的系統(tǒng)發(fā)生樹,該樹能用最小的動作產(chǎn)生序列的差異。樹的代價計算就是沿著各個分支累加最小核苷酸替換總數(shù)。具體實現(xiàn)算法為當(dāng)內(nèi)部節(jié)點的兩個直接后代節(jié)點的核苷酸相交為非空時,該節(jié)點最可能的候選核苷酸集就是這個交集;否則為它的兩個后代節(jié)點上核苷酸的并集。當(dāng)一個并集成為一個節(jié)點的核苷酸集時,通向該節(jié)點分支的某個位置上必定發(fā)生一個核苷酸替換,在統(tǒng)計過程中,只考慮信息位點的替換即可。一般來說,最大簡約法適用于以下情況:待檢序列的堿基數(shù)目大,但差異?。粔A基變異率近似相等;沒有過多的顛換、轉(zhuǎn)換傾向[2-3]。
遺傳算法是模擬生命進化的隨機搜索算法,是一種解決復(fù)雜優(yōu)化問題的通用框架。該算法首先按規(guī)則產(chǎn)生一定規(guī)模的初始群體,然后使其中的個體按一定的概率進行選擇、交叉與變異,并按照特定的適應(yīng)度評價函數(shù)進行判斷,選擇復(fù)制優(yōu)秀個體,組成新的一代,從多個可行解開始,按一定的遺傳方式進行迭代,并趨于最優(yōu)解,但其具有隨機漫游性、早熟性等特點,會使局部最優(yōu)。
模擬退火算法模擬固體退火過程,將固體加熱至高溫時,其內(nèi)能增大,內(nèi)部粒子變?yōu)闊o序狀,在其慢慢冷卻直至常溫的過程中,內(nèi)能逐漸減少,并在基態(tài)時達到最小。由初始解i和初始溫度t0開始,對當(dāng)前解重復(fù)“產(chǎn)生新解→計算內(nèi)能模擬函數(shù)差(fi-fj)→接受或舍棄”的迭代,并逐步衰減溫度t值,算法終止時,當(dāng)前解即為所得近似最優(yōu)解,其以一定的概率P接受惡化解的特性將有效規(guī)避遺傳算法的早熟性。
(3)
2.1.1 編碼方式
對于有n個物種的有根二叉樹,其中內(nèi)部結(jié)點用“+”表示,葉子結(jié)點用數(shù)字n表示,連續(xù)的2個數(shù)字和1個“+”表示一個最小的分枝。
2.1.2 適應(yīng)度函數(shù)
以簡約法中整個樹的信息位點突變數(shù)量為適應(yīng)度函數(shù),樹長為適應(yīng)值,并以歷史最大簡約樹為最終的最大簡約樹。
2.1.3 種群初始化
根據(jù)輸入的n個物種及n-1個內(nèi)部結(jié)點,隨機排列產(chǎn)生后綴編碼序列。
根據(jù)后綴法的編碼方式及特點,設(shè)計了選擇退火操作、復(fù)制操作、葉子交換退火操作、葉子與分枝交換退火操作、分枝與分枝交換退火操作。
2.2.1 選擇退火操作
采用隨機遍歷抽樣法隨機選擇兩個初始個體Pi和Pj,計算其對應(yīng)二叉樹的樹長fi和fj,如fi≤fj,則選擇個體Pi遺傳到下一代;否則,以一定的概率P接受Pj,以此保存物種的多樣性。
2.2.2 復(fù)制操作
將種群中序列按后序遍歷的方法轉(zhuǎn)換成對應(yīng)的二進制樹,再通過適應(yīng)度函數(shù)計算樹長,并根據(jù)適應(yīng)度值進行排序。
2.2.3 交換退火操作[11]
對群體中的個體,以一定的概率將指定兩位點的序列進行交換。同時要確保經(jīng)過變異后產(chǎn)生的解依然滿足二叉樹的定義,不會產(chǎn)生不可行解,因此設(shè)計以下三種交換子。
2.2.3.1 葉子交換退火操作
在親代的后綴編碼序列中隨機選擇兩個葉子位置,以2,7為例進行交換,8個葉子二叉樹及其后綴編碼如圖1所示。
圖1 8個葉子二叉樹及其后綴編碼
交換后形成新子代樹結(jié)構(gòu)及編碼序列如圖2所示。
圖2 葉子交換后二叉樹及其后綴編碼
此種交換是單位點的交換,樹的拓撲結(jié)構(gòu)沒有變化,進化速度較慢,概率值較高。計算交換前后兩個體適應(yīng)度值fi和fj,適應(yīng)度值低的絕對保留,但也以一定的概率P接受惡化的適應(yīng)度值高的解。
2.2.3.2 葉子與分枝交換退火操作
在親代的樹型結(jié)構(gòu)中隨機選擇一個葉子位置和一個分枝位置,以2和(7,8)為例,將其交換。交換產(chǎn)生的子樹結(jié)構(gòu)及編碼序列如圖3所示。
圖3 葉子與分枝交換后二叉樹
圖4 分枝與分枝交換后二叉樹
此種交換后樹的拓撲結(jié)構(gòu)發(fā)生變化,進化速度最快,為了控制結(jié)果的魯棒性,概率值應(yīng)較低。計算交換前后兩個體適應(yīng)度值fi和fj,選優(yōu)遺傳到下一代的同時,同樣以一定的概率P接受惡化解。
將新產(chǎn)生的序列和原有序列共同作為下一次迭代的初始群體。并將當(dāng)前群的最佳個體與歷史最佳個體進行比較,以低者作為歷史最佳個體。
群體迭代次數(shù)達到預(yù)定值或進化N代以后,最大簡約樹的樹長不再變化,這時,便可將最后一代進化中最佳個體作為最大簡約樹。
算法流程如圖5所示。
文中序列資料來源于https://treebase.org/treebase-web/home.html[12]的TreeBASE數(shù)據(jù)源。以算法的空間、時長和最大簡約樹的樹長為衡量標準,Phylip是一個綜合的系統(tǒng)發(fā)生分析軟件包,其中Dnapars是采用最大簡約法對原始數(shù)據(jù)和重采樣數(shù)據(jù)進行系統(tǒng)發(fā)生樹構(gòu)建,文中將改進算法與Dnapars軟件相比較。實驗環(huán)境:CPU,intel(R)Core(TM) i3-2100@3.10 GHz,RAM,8 GB,WIN32位操作系統(tǒng)。鑒于各算子對早熟現(xiàn)象的影響程度分析,經(jīng)過反復(fù)調(diào)整,算法參數(shù)如下:群體規(guī)模為500,最大迭代次數(shù)為500次,選擇算子為隨機選擇,葉子交換k1=0.003,葉子與分枝交換k2=0.002,分枝與分枝交換k3=0.001[13-14]。
圖5 改進算法流程
用簡約法構(gòu)建系統(tǒng)發(fā)生樹時,運行時間和樹長一直是兩個非常重要的指標,且隨著物種的增多和序列長度的增加,這兩個指標顯示尤為重要。樹長越小,樹越準確,時間越小,算法更優(yōu)。兩種算法實驗結(jié)果對比見表1。
表1 兩種算法實驗結(jié)果對比表
從表1中數(shù)據(jù)明顯可以得出,改進算法在這兩個指標上都有明顯提升,且隨物種和序列長度的增加,這種提升越明顯。因為在編碼過程中,后綴法編碼采用數(shù)組的方式進行存貯,使遺傳操作在線性時間內(nèi)完成,退火算子的引入又保持了物種的多樣性,跳出了局部最優(yōu)解,本算法能在時間和空間兩個維度將進化樹的構(gòu)建算法得以改進。