基于蟻群算法的雙序列比對及其實現(xiàn)

序列比對是生物信息學的基礎(chǔ)，本文首先分析了基于動態(tài)規(guī)劃的經(jīng)典雙序列比對方法，然后根據(jù)蟻群算法求解旅行商問題的思路，將旅行商問題與雙序列比對問題簡單比較了它們的異同點，并詳述了基于蟻群算法的雙序列比對步驟，根據(jù)雙序列比對過程的具體特點，對基本蟻群算法各個參量以及更新公式進行了相應的修改，成功實現(xiàn)了基于蟻群算法的雙序列比對過程。

【關(guān)鍵詞】蟻群算法 雙序列比對 更新公式

序列比對是生物信息學的基礎(chǔ)，通過在基因比對中獲得大量的序列信息，可以推斷基因的結(jié)構(gòu)、功能和進化關(guān)系。生物序列比對主要應用在生物信息學上，例如DNA雙序列比對其相似性能夠反映出兩條DNA親緣關(guān)系，實現(xiàn)生物基因識別技術(shù)；通過研究病變的生物蛋白序列和正常的序列，能夠有效地研究病變機理甚至采取一定的預防措施。在計算機領(lǐng)域，一般首先是將序列元素用不同的字符表示，一整條序列就可以用字符串來表示，因此，雙序列比對問題就是比較兩條字符串的相似性。

目前對于雙序列比對方法主要是基于動態(tài)規(guī)劃算法有經(jīng)典的Needleman-Wunsch算法和Smith-Waterman算法，利用一定的打分規(guī)則和填充分值矩陣等步驟，可以進行序列的精確比對，主要包括比對和回溯兩個過程，需要耗費大量的時間，對于較長序列，需要耗費大量存儲空間來存儲分值矩陣，因此，經(jīng)典序列比對算法對機器性能有一定的要求。對于此種方法的不足，主要有兩方面的改進與優(yōu)化，一是算法上的優(yōu)化，二是硬件上的優(yōu)化。萬文利用CUDA開發(fā)平臺實現(xiàn)了基于CPU與GPU異構(gòu)系統(tǒng)來加速Smith-Waterman算法，并且對算法進行了深入分析，設(shè)計CPU+GPU的協(xié)同并行方案，對應用程序性能與系統(tǒng)性能有顯著的提升效果。夏飛利用FPGA器件結(jié)合通用處理器實現(xiàn)了算法加速器，并且設(shè)計了細粒度并行算法，成功解決了FPGA邏輯單元和存儲資源在進行長序列比對時資源受限問題。

另一方面，已經(jīng)衍生了很多優(yōu)化算法做近似比對，比如遺傳算法，模擬退火算法等現(xiàn)代智能算法，它們不能保證每次都能找到最優(yōu)解，但是采用這些算法收斂速度快，效率高。本文詳細分析了蟻群算法求解旅行商問題，并且與雙序列比對問題做了具體的比較，得出它們的異同點，并將蟻群算法運用到求解雙序列比對問題中。

1 經(jīng)典比對算法及其實現(xiàn)

經(jīng)典的Needleman-Wunsch算法和Smith-Waterman算法對于雙序列比對分為兩個過程：比對和回溯，比對結(jié)果可以得到兩條序列的相似程度，回溯過程則可以得到序列中元素的配對結(jié)果，但回溯過程不是必須的。

1.1 雙序列比對過程

假設(shè)對于給定的兩條DNA序列，其中一條序列字符串表示為，另一條序列字符串表示為，兩條序列中各元素sj和lj分別表示DNA四種堿基{A，G，C，T}。對于這兩條序列的比對過程分三種編輯操作：替代、刪除與插入，序列比對可以通過這三種編輯操作得到最大的相似性，即最優(yōu)比對結(jié)果。

通常，為表示這個最優(yōu)比對結(jié)果會采用一定的打分規(guī)則來表示比對過程的進行以及結(jié)果分析，設(shè)打分函數(shù)為：

如果序列元素sj與lj匹配相同，而且sj和lj都不是空位，則當前元素比對結(jié)果為+2分；如果sj與lj匹配不相同，而且sj和lj都不是空位，則比對結(jié)果為-1分，如果sj與lj其中一個為空位-，則比對結(jié)果為-2分。因此，通過對兩條序列中的元素逐個比對，最終會得到一個總分值來表示雙序列的匹配程度。經(jīng)典的NW算法和SW算法就是利用打分規(guī)則和分值矩陣，填充完分值矩陣，然后是回溯過程，從分值右下角開始往上搜索局部最大分值最終可得到最優(yōu)比對結(jié)果。

2 基于蟻群系統(tǒng)的雙序列比對

2.1 問題比較與算法設(shè)計

蟻群算法是一種現(xiàn)代智能仿生算法，它通過模擬蟻群在覓食過程中尋找最短路徑的方法來求解最優(yōu)化問題，最早是應用在TSP問題上，取得了很好的效果，又來又與多種優(yōu)化算法相結(jié)合，開始得到了廣泛應用。

用蟻群算法求解序列比對問題是在該領(lǐng)域的一種嘗試，在利用蟻群算法求解TSP問題的基礎(chǔ)上，將其算法思想用于求解序列比對問題上。雙序列比對問題可以歸結(jié)于：在分值矩陣中，求解出一條分值最大且路徑最短的問題。如下，列出了利用蟻群算法求解旅行商問題和求解雙序列比對問題上的主要異同點。

（1）在TSP問題中，螞蟻的初始位置隨機放置，而且在尋找下一個節(jié)點時也是隨機選擇，而在序列比對中，螞蟻的初始位置位于分值矩陣的左上角，其搜索移動方向也只有三個：右邊，下邊和右下角；

（2）在TSP問題中，啟發(fā)信息是城市間路徑的長度，將每輪迭代中所有螞蟻走的路徑進行比較，選取其中一只螞蟻行走的最短路徑作為局部最優(yōu)結(jié)果，而在序列比對中，將螞蟻走的每一步所產(chǎn)生的字符匹配得分作為啟發(fā)信息，當所有螞蟻都走到分值矩陣的右下角時，比較每只螞蟻所走路線的最終匹配得分，選取得分最高的作為局部最優(yōu)解，并記錄其行走路線；

（3）每輪搜索的終止方式不同：在TSP問題中，城市的數(shù)目一定，根據(jù)搜索禁忌表使得每只螞蟻在每輪搜索中有且僅有一次機會就可以遍歷完所有城市，而在序列比對中，由于每只螞蟻行走的路線不同，因此，最終到達終點（分值矩陣的右下角）時存在著先后到達順序，在此問題中，可以設(shè)置終點標志位，在未到達終點前，標志位為0，到達了終點，將標志位置為1；

（4）移動方向的選取原則類似（計算轉(zhuǎn)移概率并依據(jù)輪盤賭法選擇方向），但移動方向的概率公式不同。

2.2 更新公式

因此，根據(jù)基本的蟻群算法，在雙序列比對問題中需要修改的參量以及更新公式如下：

2.3 算法實現(xiàn)

利用蟻群算法求解雙序列比對問題詳細步驟如下：

步驟1：規(guī)則設(shè)定。輸入兩條待比對序列，并設(shè)定其打分規(guī)則，再根據(jù)分值構(gòu)建對應的比對字符匹配矩陣，將分值矩陣中的正數(shù)在字符匹配矩陣中相應的位置設(shè)為原值，負數(shù)取絕對值的倒數(shù)，則字符相同的位置為2，不相同的為1，有空位的為0.5；

步驟2：初始化。初始化m只人工螞蟻，并設(shè)置各參數(shù)：迭代次數(shù)、方向權(quán)值、初始化信息素矩陣等，從分值矩陣的左上角出發(fā)向矩陣右下角按一定的規(guī)則開始搜索；

步驟3：有規(guī)則搜索。螞蟻在前進搜索過程中，其搜索移動的方向只能三個：向下、向右和向右下方（以螞蟻當前位置為參考點）。首先設(shè)定一個固定值q0，螞蟻在選擇前進方向時，先產(chǎn)生一個隨機數(shù)，如果這個隨機數(shù)小于q0，則螞蟻選擇右下方位置移動，如果這個隨機數(shù)大于q0，則采用輪盤賭法來進行選擇移動方向。

步驟4：移動一步之后，判斷螞蟻是否走到邊界。螞蟻如果走到最右邊，則之后的移動方向只能選擇向下；如果走到最下邊，則之后的移動方向只能選擇向右，因此，考慮到此種邊緣情況，下一次迭代時，步驟3就不再需要計算概率。

步驟5：當螞蟻移動到了分值矩陣的右下角，設(shè)定其結(jié)束標志位，結(jié)束當前一輪的搜索，同時記錄螞蟻的行走路線。

步驟6：當m只螞蟻都搜索完之后，比較它們的行走路徑和得分，與歷史局部最優(yōu)解比較，更新當前局部最優(yōu)解，并保存這一輪的最高分及其對應的路徑（局部最優(yōu)解），然后更新全局最優(yōu)解，同時更新信息素矩陣。

步驟7：判斷當前迭代次數(shù)是否達到預先設(shè)定的次數(shù)，或者是，全局最優(yōu)解連續(xù)T次未更新（更新為相同的值），如果條件滿足，則算法結(jié)束，否則，進行下一次迭代，跳轉(zhuǎn)到步驟2。

最終可根據(jù)全局最優(yōu)解對應的螞蟻行走路線可推出最佳序列比對結(jié)果，如果螞蟻是向下方移動，則在序列L中對應的位置插入-，如果螞蟻是向右方移動，則在序列S中對應的位置插入-，如果螞蟻向右下角方向移動，則說明當前匹配的序列元素相同。

由此可以看出，利用蟻群算法求解雙序列比對問題之后，不再需要回溯過程，而且在比對過程中，每只螞蟻的行走路線不必太過依賴于前一步所走的位置，僅僅只是根據(jù)概率公式進行一個方向的選擇而已，從一定程度上體現(xiàn)了蟻群算法的搜索隨機性。

3 實驗結(jié)果與分析

C語言編程算法流程，并利用蟻群算法進行DNA雙序列比對，采用GCC編譯器，運行環(huán)境：Intel（R） i3-4160 CPU@3.60GHz。蟻群算法的各參數(shù)設(shè)置如下：迭代次數(shù)NcMax =100，蟻群規(guī)模M =20，信息素重要Alpha = 5.0， 字符匹配重要程度Beta = 3.0，空位數(shù)目重要程度 Gamma = 2.0，轉(zhuǎn)移概率設(shè)定值q0= 0.2，信息素揮發(fā)程度Rou = 0.05，信息素增加強度系數(shù)Q = 0.1， 信息素增量轉(zhuǎn)化參數(shù)Qc = 50，方向權(quán)重參考值H = 5，比對結(jié)果如下：

在GenBank數(shù)據(jù)庫中選取兩條DNA序列進行比對，驗證算法的可行性，其編號分別是DQ482171和DQ482174，序列長度為382，圖1顯示了比對過程中，局部最優(yōu)解與迭代次數(shù)的動態(tài)變化關(guān)系。

由以上結(jié)果分析可知，隨著迭代次數(shù)的增加，最優(yōu)比對分值結(jié)果不斷增加，說明螞蟻在搜索過程中，信息素的更新起到了一定的作用，對螞蟻行走的路線不斷優(yōu)化。然后再進行重復實驗10次，表1顯示了對于不同長度的幾組序列利用蟻群算法比對結(jié)果：最優(yōu)結(jié)果、最差結(jié)果、10次實驗的算術(shù)平均值，以及程序運行十次的累計耗時。

其中編號1的雙序列基因序列號分別是DQ482171和DQ482174；編號2的雙序列基因序列號分別是AB023276和AB023278；編號3的雙序列基因序列號分別是AF039225和AF039226，其中AF039225長度為613bp，而AF039226長度為612bp，以觀察不同長度序列利用蟻群算法比對的結(jié)果；編號4的雙序列對應的基因序列號分別是HIV1U39234和HIV1U39238，其長度較長，以觀察在長序列中算法的適用性；編號5的雙序列是對編號2序列進行了一定的修改，降低序列的相似度，以觀察在低相似度中算法的適用性。

在進行長序列比對時，可以看到蟻群算法并不能保證求得最優(yōu)解，這是因為，隨著序列的增長，螞蟻在前期搜索過程中如果偏離了最優(yōu)路線，導致后期搜索很難收斂，甚至遠離最優(yōu)解。由此可以看出，長序列不如短序列的比對效果好。在進行不同長度序列比對時，蟻群算法仍具有很好的收斂性，并盡可能的求得最優(yōu)解。在進行低相似度序列比對時，蟻群算法尋得路徑也可能會導致并不是最優(yōu)解，這是因為，在搜索時，如果螞蟻盡可能的保證最優(yōu)解而逐漸遠離實際的最優(yōu)解，不同字符匹配的數(shù)目越來越多，卻增加了螞蟻搜索的隨機性，使得這一影響更為明顯，搜索不到最優(yōu)解的可能性也會越大。

由以上結(jié)果可以看出，利用蟻群算法求解雙序列比對問題是可實現(xiàn)的，可以看到，在利用蟻群算法求解雙序列比對問題時，還是具有一定的隨機性，并且隨著迭代的次數(shù)的增加，算法開始收斂直到最優(yōu)結(jié)果，至于最優(yōu)結(jié)果的最終表達取決于字符匹配過程中對規(guī)則的設(shè)定。

4 結(jié)束語

文中對比了TSP問題和雙序列比對問題，并且成功將蟻群算法實現(xiàn)應用在序列比對之中，對于更新公式的替換，在很多參考文獻中稱之為改進的蟻群算法，在這方面還可以有很多值得研究的地方，而對于蟻群算法的參數(shù)尋優(yōu)卻仍是一個很有意義的研究方向，另外針對蟻群算法易陷入局部最優(yōu)結(jié)果的弊端，此處也可以利用其它的算法來優(yōu)化。如果對這幾個方面有更大的改進，蟻群算法應該會使得雙序列比對結(jié)果更加精確。

（通訊作者：汪楠）

參考文獻

[1]Smith T.F.，Waterman M.S.Identification of Common Molecular Subsequences[J].Journal of Molecular Biology，1981，147（01）.

[2]萬文.生物序列分析算法的CPU+GPU異構(gòu)并行優(yōu)化關(guān)鍵技術(shù)研究[D].國防科學技術(shù)大學，2012（03）.

[3]楊春燕，鐘誠.CPU和GPU協(xié)同并行加速多生物序列比對[J].小型微型計算機系統(tǒng)，2016（12）.

[4]夏飛.生物序列分析算法硬件加速器關(guān)鍵技術(shù)研究[D].國防科學技術(shù)大學，2011（03）.

[5]李方潔.基于智能算法的DNA序列比對研究[D].山東師范大學，2011（06）.

[6]李娟.生物序列相似性搜索算法研究與實現(xiàn)[D].華南理工大學，2013.

作者簡介

李宏周（1980-），男，桂林電子科技大學碩士生導師。研究方向為電路與系統(tǒng)。

汪楠（1993-），男，湖北省黃岡市人。桂林電子科技大學研究生。主要研究方向為嵌入式系統(tǒng)設(shè)計。

作者單位

1.桂林電子科技大學 計算機與信息安全學院 廣西壯族自治區(qū)桂林市 541004

2.桂林電子科技大學 電子工程與自動化學院 廣西壯族自治區(qū)桂林市 541004

3.桂林電子科技大學 生命與環(huán)境科學學院 廣西壯族自治區(qū)桂林市 541004