基于GPU的遺傳退火多序列比對并行研究

韋樹烽，劉 羽，蔣財運

（1.桂林理工大學廣西礦冶與環(huán)境科學實驗中心，廣西桂林541004；2.桂林理工大學機械與控制工程學院，廣西桂林541004；3.桂林理工大學信息科學與工程學院，廣西桂林541004）

0 引 言

在生物信息學中，多序列比對（multiple sequences alignment）是一種基本的生物信息處理方法，在發(fā)現(xiàn)生物序列的結(jié)構(gòu)、功能和進化的信息方面具有重要意義［1］，它是一個尚未有效解決的NP完全的組合優(yōu)化問題。生物科學的不斷發(fā)展使得GenBank、EMBL、DDBJ等序列數(shù)據(jù)庫容量快速增長。截止2013年2月GenBank庫核酸序列已超過1.6億條，堿基對超過1.5千億對［2］。對于序列比對的結(jié)果評估具有計算速度和最佳敏感度兩個標準［3］，設計一個合理高效的序列比對算法是一個重要的研究方向。

為了提高比對效率，并行計算技術(shù)被應用到該領域。目前主要有BLAST、FASTA、Smith－Waterman和CLUSTALW等MSA的并行處理策略［4－7］。由于近年來GPU應用于并行計算領域上的優(yōu)勢，國內(nèi)外眾多測序比對軟件都在嘗試基于GPU的平臺進行研發(fā)。而本文主要工作就是基于CUDA模型研究遺傳退火多序列比對算法的并行加速模型，為生物信息研究分析提供支撐。

1 多序列比對算法

1.1 多序列比對算法數(shù)學描述

表1 多序列比對示例

1.2 遺傳模擬退火多序列比對算法概述與設計

在MSA算法中，主要有動態(tài)規(guī)劃法、漸進算法、迭代算法等。其中遺傳模擬退火多序列比對算法（SAGAMSA）是迭代算法中的典型全局比對算法，具有極高計算復雜性［8］。它通過精細迭代方式來實現(xiàn)多序列比對，利用模擬退火策略對遺傳操作進行改進使MSA可以獲得更好的效果，同時具有很好的魯棒性、對序列個數(shù)不敏感和不易陷入局部最小的特性。SAGA－MSA可以定義為一個8元組［9］

式中:C——個體的編碼方式，P——初始種群，F(xiàn)——個體的適應度評價函數(shù)，M——種群大小，Φ——選擇退火算子，?！徊嫱嘶鹚阕樱贰儺愅嘶鹚阕?，T——終止退火條件。下面簡要介紹本文幾個關(guān)鍵點:

（1）本文采用可變長的二維矩陣編碼方式，由1.1節(jié)知m條序列比對結(jié)果是一張m＊n矩陣，編碼矩陣為SEQ′ij（0＜i＜m，0＜j＜n），這種編碼方式具有簡單直觀、易于操作的優(yōu)點。

（2）初始化種群是SAGA－MSA搜索尋優(yōu)的出發(fā)點，本文采用隨機產(chǎn)生初始種群方法，對長度為Li（0＜i＜m）的SEQi（0＜i＜m）隨機選擇不大于Li的整數(shù)Xl（0≤l≤n－Li），在序列Xi位置插入“－”使得比對結(jié)果每條序列的長度都為n。

（3）SAGA－MSA的適應度函數(shù)采用SP目標函數(shù)，函數(shù)值越大則種群個體的質(zhì)量越好，進入下一代的概率就越高。目標函數(shù)是總評分減去空位仿射罰分［10］，總評分中，兩個空位符匹配評分為0，兩個非空位字符匹配時參照Blosum62評分矩陣［11］，空位與其它字符的比對評分值為－5?？瘴环律淞P分中，空位開放罰分為10，空位延伸罰分為0.2。

（4）選擇退火算子是在遺傳操作中加入了退火策略，采用Metropolis接受準則。為了保持種群多樣性，操作中不刪除適應度不夠好的種群個體。以概率pnew接受新一代個體，保存當代最佳個體

式中:f——適應度函數(shù)，t＞exp（t＊alpha），（alpha為溫度變化系數(shù)）為退火溫度控制參數(shù)。

在農(nóng)業(yè)機械化法律關(guān)系中，出現(xiàn)了責任主體和監(jiān)督主體不明晰的現(xiàn)象。目前在全國范圍內(nèi)，我國尚未建立一個完善的農(nóng)業(yè)機械安全管理法律制度。針對農(nóng)機管理部門對農(nóng)業(yè)機械的初次檢驗、年度檢驗和臨時檢驗；農(nóng)業(yè)機械的停駛登記、報廢登記；農(nóng)業(yè)機械駕駛操作人員的考核、安全教育、證照核發(fā)、變動登記、年審、農(nóng)機違章和責任事故的處理等一系列的問題在各個地區(qū)之間存在的差異不大，可以在全國范圍內(nèi)確立統(tǒng)一的標準，由政府主管的農(nóng)業(yè)部門作為主要監(jiān)管機構(gòu)，聯(lián)合交通運輸部門和公安部門，形成一個較為完善的農(nóng)業(yè)機械管理體系，對農(nóng)業(yè)機械化的安全管理問題明確一個相對清晰的責任劃分，方便農(nóng)業(yè)機械的安全管理。

（5）交叉操作是將老一代的序列進行重組，形成優(yōu)秀的序列遺傳到新一代。本文采用兩點交叉的遺傳操作，對老一代種群個體隨機從［0，Li－1］（0＜i＜m）中選取兩個整數(shù)rand1＞rand2作為交叉位置，以概率pnew對個體進行交叉操作產(chǎn)生新代。

1.3 SAGA－MSA算法串行測試與分析

根據(jù)1.2節(jié)的描述對SAGA－MSA實現(xiàn)串行算法如下:

對上述串行算法各階段耗時分析，測試數(shù)據(jù)來源于GenBank的序列文件gbbct1.fsa＿aa中隨機抽取的2048條序列，平均長度為256。設定初始參數(shù)種群數(shù)為5、迭代次數(shù)為5、初始溫度為10、溫度變化系數(shù)為0.5、終止溫度為1。測試結(jié)果見表2，算法的第三、四階段耗時比例較大，占用大部分的計算時間。

表2 SAGA－MSA串行算法各階段耗時百分比

2 GPU架構(gòu)和CUDA編程模型

自2007年NVIDIA推出通用并行計算模型CUDA以來，GPU（圖形處理器）能夠解決復雜計算問題而被應用到高性能并行計算領域。CUDA是用于GPU計算的開發(fā)環(huán)境，它是一個全新的軟硬件架構(gòu)，包含CUDA指令集架構(gòu)（ISA）和GPU內(nèi)部的并行計算引擎兩部分［12］。在CUDA的模型中可以將GPU視為一個并行數(shù)據(jù)計算設備，對計算進行分配和管理，具有并行效果明顯的特點。

在本文并行模型中，GPU專注于復雜的并行處理計算。由于CUDA內(nèi)部并行機制的多級性，算法采用任務分解方式，將需要并行處理的粗粒度任務劃分成相對獨立的模塊，不同模塊由不同的GPU線程塊處理，同時盡量減少通信。線程塊內(nèi)由多個處理線程進行細粒度的計算。因此如何準確劃分合理的求解粒度和解決模塊間通信問題是我們研究的重點。

3 SAGA－MSA算法并行化分析與實現(xiàn)

本文的研究目標是在保持序列比對結(jié)果在一定的敏感度的前提下，提高其比對效率。SAGA－MSA算法并行化的關(guān)鍵技術(shù)難點在設計合理的并行粒度劃分策略，解決并行線程的通信延遲問題、保持良好的負載平衡，以此保證MSA并行化實現(xiàn)的合理性和高效性。

3.1 算法并行化分析

表2的串行分析結(jié)果顯示，算法的stage 3和stage 4耗時最為明顯，當數(shù)據(jù)集很大時，stage 1、stage 2和stage 5幾乎可以忽略不計。因此本文主要針對stage 3和stage 4這兩階段的運算進行并行化。SAGA－MSA算法繼承了自然進化過程中固有的并行性，算法中的迭代不是從單個解開始計算，而是從問題解集（種群）開始。計算過程從一組解迭代到另一組解，能夠同時處理群體中多個個體，易于并行化。而本文擬采用SAGA－MSA的種群計算的粗粒度模型（coarse－grained model）［13］結(jié)合多序列比對中固有的細粒度模型（fine－grained model）來研究算法的并行化實現(xiàn)。

在SAGA－MSA中，粗粒度并行模型就是種群中若干個體并行計算，如圖1所示。每個個體分配對應的線程塊執(zhí)行并行運算，每個線程塊進行獨立的選擇、交叉和變異操作和適應度計算。定期地相互傳送適應度最好的個體，加快滿足終止條件的要求。這樣算法的并行化在不同的種群規(guī)模時都以較小的通信開銷獲得了更好的加速。

圖1 SAGA－MSA的粗粒度并行模型流程

由于MSA是雙序列比對的拓展，具有天然并行性［14］。根據(jù)兩條序列比對與另兩條序列比對的數(shù)據(jù)不相關(guān)特性，可是確定多序列比對的細粒度模型。既是從m條序列SEQ′ij（1≤i≤m－1，0＜j＜n）中選擇的SEQ′ij與SEQ′（i＋1）j成為一組比對，這樣比對組數(shù)是m＊（m－1）/2。分別對每一組分配不同的處理線程進行處理來實現(xiàn)并行化，如圖2所示。相對于串行的循環(huán)比對各組的最大優(yōu)點就是可以實現(xiàn)同時對多個組并行處理，理論上只要有足夠的處理核，可以實現(xiàn)很好的加速效果。

圖2 SAGA－MSA的細粒度并行策略

同時細粒度模型還有一層數(shù)據(jù)處理并行化，例如圖2的SEQ′4j和SEQ′5j（j＝1，2，…n）兩條序列比對時進行并行化，它們的第j對氨基酸比對，用兩個線程分別對兩條序列各自搜索實現(xiàn)并行處理。另外，在對m條序列生成比對矩陣時每條序列操作也都是并行操作的。

當然算法的并行化研究還需根據(jù)輸入數(shù)據(jù)本身的特點進行合理并行化。因此選擇輸入序列的種群規(guī)模p和序列條數(shù)m是處理核的整數(shù)倍時，理論上能達到較好的負載均衡。另外輸入序列的平均長度也影響著負載均衡，假如對輸入序列長度引入數(shù)學期望E，有E增大，則線程塊和線程處理等待同步操作時耗時增加，各并行處理線程能達到的負載均衡變差，序列比對效率有所下降。

3.2 算法并行化實現(xiàn)

算法的并行化代碼是C/C＋＋結(jié)合CUDA編程模型的開發(fā)語言。結(jié)合CUDA內(nèi)部并行機制與3.1節(jié)分析的并行方式設計了合理的三層并行算法。CUDA提供了HOST端調(diào)用DEVICE端的全局函數(shù)（＿＿global＿＿），調(diào)用格式例如調(diào)用函數(shù)change＜＜＜dimBlock1，dim Thread1，smem－Size1＞＞＞（），其中change函數(shù)是1.3第五點實現(xiàn)的種群交叉操作。CUDA中在每個線程塊的執(zhí)行寬度warp是32個線程，線程數(shù)的設置為32的倍數(shù)能達到最好效率。

計算適應度Fitness函數(shù)（如圖3所示）是并行實現(xiàn)的重點，函數(shù)的并行性能直接影響并行效率。這一階段分為兩個并行函數(shù)操作，第一個函數(shù)用CUDA線程塊內(nèi)的兩層并行級分別同時計算m＊（m－1）/2個序列比對組的適應度并保存到數(shù)組sum12。第二個是并行規(guī)約數(shù)組sum12［m＊（m－1）/2］，在規(guī)約運算中，大部分的warp都能滿載，可以有效提高效率。函數(shù)SafeTheBest的作用是找出最優(yōu)的適應度并保存最優(yōu)種群，根據(jù)概率pnew（見式（2））接受上一代種群。算法完成一代種群計算后，自動進入下一代的種群更新，直至滿足退火終止條件得到最優(yōu)適應度和種群后才退出，這相對于串行算法每計算一個種群個體都比較適應度和保存種群得到了很大的優(yōu)化。

圖3 Fitness函數(shù)并行實現(xiàn)

另外為了記錄算法的各階段和總程序的耗時情況，構(gòu)建了一個精確的時鐘函數(shù)，調(diào)用核函數(shù)后在HOST端用cuda Thread Synchronize（）來實現(xiàn)等待GPU時間同步，這樣可以準確地記錄程序的總耗時和各階段的運算時間。

4 測 試

在本節(jié)中，對SAGA－MSA串行算法和并行算法進行全方面的比較。測試環(huán)境以下:CPU是3.30 GHz的英特爾酷睿i3－2120，GPU為NVIDIA GeForce 405（512MB的RAM），操作系統(tǒng)是Windows XP（2GB內(nèi)存）。測試數(shù)據(jù)來源于GANBANK庫2013年2月版的gbbct1.fsa＿aa序列數(shù)據(jù)文件。從中隨機選擇序列并編制了20個序列數(shù)量從8到1024測試實例。其中每個測試實例有5種序列長度，平均長度為49、99、203、227、412，比對結(jié)果的序列長度為64、128、256、450、720。為了確保兩種算法的等效計算，使用了下列3組不同的起始輸入?yún)?shù):第一組種群大小為5、迭代次數(shù)為5、起始溫度為10、退火系數(shù)為0.5、終止溫度為1；第二組是第一組只增加迭代次數(shù)和起始溫度到10、20；第三組是在第二組基礎上只增加種群大小到20。另外為了保證結(jié)果不受系統(tǒng)不穩(wěn)定形成的誤差影響，每個實例運行5次取其平均值得到穩(wěn)定的結(jié)果。

首先對兩種算法的適應度（敏感度）進行比較。圖4是兩種算法在上述第一組輸入?yún)?shù)情況下得到的適應度比較圖，串行算法用實線表示，并行算法用虛線表示，不同的線上符號表示5種序列長度。可以觀察的到，串行與并行算法的適應度分別在上述的輸入序列長度各不相同時都能幾乎保持一致，隨著序列數(shù)量的增加，適應度都是以一定的趨勢在增加。應當強調(diào)的是，因為算法的隨機選擇機制，SAGA－MSA并不是一個等效的算法，它們的比對結(jié)果不完全相同。因此本文的并行算法能滿足與串行算法在相同敏感度的基礎上研究加速的前提。

圖4 串行算法與并行算法的適應度比較

圖5是兩種算法的總耗時的時間曲線圖，執(zhí)行時間以s為單位。兩個算法的執(zhí)行時間都是隨著序列規(guī)模增大（序列數(shù)量和序列長度）而增大，顯然串行算法的增長速度更快，例如序列長度為720時，它的平均輸入序列長度412，則插入的平均空位數(shù)是308?？梢杂^察到并行算法雖然在開始問題規(guī)模較小時耗時比串行的多，但是增加序列條數(shù)時耗時是平緩增加，串行耗時卻增長很快，當增加到704條序列以后并行算法的耗時開始比串行的時間耗時小。

圖5 串行算法與并行算法的總耗時比較

圖6是兩種算法耗時加速比曲線圖，可以觀察到不同序列規(guī)模有不同的加速比。比對結(jié)果長度如64、128、256的平均插入“－”空位數(shù)是15、29、53的加速效果明顯，加速比提升較快，表明在輸入序列數(shù)量不多的情況下也能有較好的加速效果，長度為128時的加速效果最明顯，但是由于受限于顯卡的性能，從序列數(shù)量為832后幾乎是平緩增長的，并在序列數(shù)量為1024時加速比能達到12.53倍。另外對于長度為450和720插入平均空位數(shù)是223和308的加速效果不太明顯，加速比只達到了2.19和1.57。影響它們的加速效果是由于輸入序列的最大長度和最小長度之間的相差很大，并且有多數(shù)的長序列和短序列交雜輸入，由3.1節(jié)的輸入數(shù)據(jù)分析可知這樣會造成在并行運算時不同的序列處理線程等待同步耗時太長。

圖6 串行算法與并行算法的加速比

最后，為了測試分析并行算法在增加種群數(shù)量和迭代次數(shù)對并行效率的影響，特別增加了第二、三組的參數(shù)測試，測試的加速比都與第一組的對比。在這兩個測試中，發(fā)現(xiàn)在圖7和圖8的加速效果都有些許增加，圖8比圖7更為明顯。并行算法在增加迭代次數(shù)和種群數(shù)量時都能保持良好的加速效果，且能更快達到加速極限值。如圖8中，序列長度為128時加速比平均增加了0.06，并在序列數(shù)量為384時能達到最大增加值0.19，因此初始參數(shù)的設定對加速效果也有不同程度的影響。

圖7 增加迭代數(shù)和初始溫度的加速比變化

圖8 增加種群數(shù)量的加速比變化

5 結(jié)束語

本文從多個層面進行了SAGA－MSA并行化的合理性和高效性處理，采用了合理的并行粒度劃分策略，并從線程同步、存儲訪問和負載均衡等方面對算法進行了優(yōu)化，用CUDA的共享內(nèi)存來緩解訪存操作的競爭影響，在處理插入空位和比對空位時排除了大部分的冗余同步操作。通過在運算開始就進行數(shù)據(jù)塊的初始化，減少了數(shù)據(jù)間交換引起的通信延遲問題。通過對序列插入空位位置進行排序，保證了良好的負載均衡。測試結(jié)果表明，實現(xiàn)的基于GPU的SAGA－MSA并行高效算法，相對于串行算法，計算速度有了明顯的提升。雖然受限于測試GPU硬件性能，還是實現(xiàn)了12.53倍的加速比，充分發(fā)揮了GPU在并行運算中的性能優(yōu)勢。

下一步，將考慮引入基于消息傳遞的MPI和基于共享內(nèi)存的OPENMP并行機制來實現(xiàn)多層次的混合并行模型，CPU多核和GPU眾核分別負責粗、細粒度的并行計算任務。粗細粒度相結(jié)合的并行方式可以更好地利用計算資源，從而進一步提高加速比。

［1］ZHOU Hong.Research on parallel algorithm of multiple DNA sequence alignment based on de Bruijn graph［D］.Tianjin:Tianjin University，2011（in Chinese）.［周紅.基于de Bruijn圖的DNA多序列比對并行算法研究［D］.天津:天津大學，2011.］

［2］GenBank.Growth of GenBank and WGS［DB/OL］.http://www.ncbi.nlm.nih.gov/genbank/statistics，2013.

［3］KANG Xiaojun.Gene sequence alignment algorithm research and implement in SNP［D］.Wuhan:Huazhong Agricultural University，2011（in Chinese）.［康曉軍.基因序列比對算法在SNP中的研究及應用［D］.武漢:華中農(nóng)業(yè)大學，2011］.

［4］Ali Khajeh－Saeed，Stephen Poole，J Blair Perot.Acceleration of the Smith－Waterman algorithm using single and multiple graphics processors［J］.Journal of Computational Physics，2010，229（11）:4247－4258.

［5］Lars Wienbrandt，Stefan Baumgart，Jost Bissel，et al.Massively parallel FPGA－based implementation of BLASTp with the two－h(huán)it method［C］//Proceedings of the International Conference on Computational Science，2011:1967－1976.

［6］Cheng Ling，Khaled Benkrid.Design and implementation of a CUDA－compatible GPU－based core for gapped BLAST algorithm［C］//Proceedings of the International Conference on Computational Science，2010:495－504.

［7］Baomin Xu，Jin Gao，Chunyan Li.An efficient algorithm for DNA fragment assembly in MapReduce［J］.Biochemical and Biophysical Research Communications，2012，426（3）:395－398.

［8］ZHANG Yi.Based on genetic annealing bioinformatics multiple sequence alignment algorithm［D］.Chengdu:University of Electronic Science and Technology，2006（in Chinese）.［張憶.基于遺傳退火的生物信息學多序列比對算法研究［D］.成都:電子科技大學，2006.］

［9］BIAN Xia，MI Liang.Development on genetic algorithm theory and its applications［J］.Application Research of Computers，2010，27（7）:2425－2429（in Chinese）.［邊霞，米良.遺傳算法理論及其應用研究進展［J］.計算機應用研究，2010，27（7）:2425－2429.］

［10］Chengpeng Bi.Deterministic local alignment methods improved by a simple genetic algorithm［J］.Neurocomputing，2010，73（13）:2394－2406.

［11］Jacek Blazewicz，Wojciech Frohmberg，Michal Kierzynka，et al.G－MSA－－A GPU－based，fast and accurate algorithm for multiple sequence alignment［J］.Journal of Parallel and Distributed Computing，2013，73（1）:32－41.

［12］BAI Hongtao.Research on high performance parallel algorithms based on GPU［D］.Jilin:Jilin University，2010（in Chinese）.［白洪濤.基于GPU的高性能并行算法研究［D］.吉林:吉林大學，2010.］

［13］TAN Caifeng，MA Anguo，Xing Zuocheng.Research on the parallel implementation of genetic algorithm on CUDA platform［J］.Computer Engineering and Science，2009，31（A1）:68－72（in Chinese）.［譚彩鳳，馬安國，邢座程.基于CUDA平臺的遺傳算法并行實現(xiàn)研究［J］.計算機工程與科學，2009，31（A1）:68－72.］

［14］YE Xiaochun，LIN Wei，F(xiàn)AN Dongrui，et al.Efficient parallelization and optimization of protein sequence comparison algorithm on many－core architecture［J］.Journal of Software，2010，21（12）:3094－3105（in Chinese）［葉笑春，林偉，范東睿，等.蛋白質(zhì)序列比對算法在眾核結(jié)構(gòu)上的并行優(yōu)化［J］.軟件學報，2010，21（12）:3094－3105.］