林 敏, 鐘一文, 林 娟

(福建農(nóng)林大學(xué)計算機與信息學(xué)院，福建 福州 350002)

?

基于GPU并行的生物序列局部比對算法

林 敏, 鐘一文, 林 娟

(福建農(nóng)林大學(xué)計算機與信息學(xué)院，福建 福州 350002)

對Smith-Waterman算法的計算公式進行了改進以適應(yīng)GPU并行的特點,并提出新的基于BLOCK分塊的并行前綴掃描法；通過UP-DOWN步驟、BLOCK間調(diào)整、Eij微調(diào)等步驟在O(logn)時間內(nèi)計算出行中每一個元素的前綴最大值；最后將回溯過程置于GPU端，避免了CPU與GPU間內(nèi)存的拷貝.與傳統(tǒng)的Smith-Waterman算法相比,該算法在低端的GPU平臺性能提升90倍；與同樣基于GPU的SWAT算法相比，性能也有較大的提升.

Smith-Waterman算法; 并行前綴掃描; 通用圖形處理器; 序列比對

生物序列比對需在2個序列中查找相似序列，以便對基因或蛋白質(zhì)所具備的功能做出分析[1].常見的比對算法有BLAST[2]、FASTA, Smith-Waterman[3]等算法.BLAST算法和FASTA算法運行效率較高，但其敏感度均不如Smith-Waterman算法.Smith-Waterman算法是一種局部比對算法，它具有匹配敏感度高的特點，但其運算效率為O(m*n).而在生物序列比對中，為了尋求更高的準(zhǔn)確度和敏感性，常常又需要選用Smith-Waterman算法進行更精確的匹配.Zhao et al[4]選用Smith-Waterman算法進行生物序列比對.

Smith-Waterman算法被提出以來，研究者提出了各種改進的方案以提高其執(zhí)行效率.Wozniak[5]提出了采用Wave Front方式計算動態(tài)規(guī)劃矩陣(DP).這種方式是利用對角線元素可以并行計算的特性，但對角線的并行化程度受制于對角線的長度，只有當(dāng)連接2個角的對角線時才能達到最大的并行，所以大部分情況下大量線程處于空閑狀態(tài).并且，根據(jù)GPU全局內(nèi)存的合并訪問原理，這種以對角線方式訪問矩陣元素的方法，將被拆分成多次訪問，不適合GPU的并行.

Farrar[6]提出了一種基于單指令多數(shù)據(jù)流(single instruction multiple data, SIMD)的加速方式，但是最多只有6倍的加速.Akoglu[7]提出了第1個基于統(tǒng)一計算設(shè)備架構(gòu)(compute unified device architecture, CUDA)的加速方案，最快有20倍的加速，但在一般情況下只能達到8-10倍的加速.KHAJEH-SAEED[8]提出基于單GPU或GPU的并行加速方式，通過改造Smith-Waterman算法的計算方法，實現(xiàn)了并行計算F值；但其在計算E值時采用了并行前綴掃描[9-10]的策略，這種掃描算法并未考慮GPU存儲器的特點，并行的效率未大幅提高.馬海晨等[11]和Liu et al[12]提出了粗粒度的并行，將1條序列同時與多個序列組成的數(shù)據(jù)庫序列進行比對，這種算法的本質(zhì)還是雙序列比對，有待于通過對雙序列比對算法的優(yōu)化進一步提高其性能.因此，本文結(jié)合GPU并行的特點，對Smith-Waterman的并行步驟做了改進，提出了新的GPU加速的Smith-Waterman方案.

1 GPU并行的Smith-Waterman算法

1.1 Smith-Waterman計算公式的改進

Smith-Waterman算法是基于局部聯(lián)配的動態(tài)規(guī)劃算法，DP中每一元素的得分由其對角線、列和行方向的得分決定，建立的步驟如式(1)所示：

(1)

圖1 構(gòu)造DPFig.1 Creating DP

式中，Hij表示DP中第i行第j列的得分，Sij表示第i行第j列堿基的比對得分，Gs表示開放空位得分，Ge表示擴展空位得分.根據(jù)比對的第1個步驟建立DP，建立DP的步驟就是序列比對的過程，DP中每一元素由3個方向的值決定.第2個步驟通過DP尋找最長相似子序列，先在DP中查找最大值，由最大值元素位置開始回溯，直到得分為0為止，其回溯所經(jīng)歷的路徑就是其比對結(jié)果.比對過程中需要加入空位罰分、錯配罰分、空位擴展罰分，以使匹配結(jié)果更接近生物特性.圖1為構(gòu)造DP矩陣的過程；圖2表示回溯的過程.

圖2 回溯的過程

式(1)中的計算步驟不適合GPU的并行，Khajeh-saeed et al[8]對式(1)進行了重新設(shè)計，使GPU的線程可以并行計算DP矩陣的每列或每行，其計算步驟如下：

通過式(2)、式(3)先求每一列的最大值，采用式(4)計算行最大值,最后通過式(5)計算Hij.

Fij=Max(Fi-1,j,Hi-1,j-Gs)-Ge

(2)

Uij=Max(Hi-1,j-1+Sij,Fij,0)

(3)

Eij=Max1≤k

(4)

Hij=Max(Eij-Gs,Uij)

(5)

Fij和Eij分別表示列和行最大值，Hij為DP矩陣，Sij為比對得分，Gs和Ge分別表示空位罰分和空位擴展罰分.本文中的算法是根據(jù)GPU全局內(nèi)存合并訪問的特點，采用了列并行的方式.同時為了采用GPU并行前綴掃描算法計算Eij，將式(3)、式(4)重新設(shè)計成式(6)、式(7)，無需再保存Uij變量，需要的存儲空間更少.

Eij=Max(Hi-1,j-1+Sij,Fij,0)-(n-j)Ge

(6)

(7)

1.2GPU并行Smith-Waterman算法的實現(xiàn)

1.2.1GPU并行計算架構(gòu)平臺 本文的GPU并行算法是基于CUDA架構(gòu)的，而CUDA是一種通用并行計算架構(gòu)，它包含了指令集架構(gòu)以及GPU內(nèi)部的并行計算引擎，該架構(gòu)使GPU能夠解決復(fù)雜的計算問題.CUDA采用單指令多線程(SIMT)的執(zhí)行模型，它同時能產(chǎn)生大量的線程，通過大量線程的并行計算來隱藏存儲器訪問延遲，而不必使用較大的數(shù)據(jù)緩存[11].

運行在CUDA中的程序稱為Kernel,每個Kernel相當(dāng)于1個計算網(wǎng)格Grid，每個Grid由一維、二維或三維線程塊(BLOCK)組成,而每個BLOCK被組織成一維、二維或三維的線程網(wǎng)格.BLOCK內(nèi)的所有線程能共同使用共享存儲器，且其中的所有線程同步執(zhí)行，而BLOCK間的線程共享全局存儲器不能直接同步，只能通過重新啟動Kernel以達到同步的效果.

1.2.2 算法的主要結(jié)構(gòu) 算法的核心部分是并行計算DP矩陣中每行的F值和E值，然后最終計算出H值.為了計算F值，將DP矩陣以一維數(shù)組方式存儲在GPU端，以較長序列作為矩陣的行，以較短序列作為矩陣的列.采用列并行的方式計算的F值剛好滿足合并訪問(圖3)的要求，并且每次無需記錄所有行的F值，僅保留上一次計算的F值即可.而E值的計算則較為復(fù)雜，本文提出了一種基于BLOCK的并行前綴掃描算法，可快速計算出每行的E值，整個算法的執(zhí)行步驟如算法1所描述.

圖3 合并內(nèi)存訪問(實線部分為可以1次合并訪問，虛線部分則要分2次完成)

算法1：主程序

begin

載入序列，將序列二值化壓縮存儲;

在GPU端為Hij,Fij,Eij分配存儲空間;

計算blocknum=ceil(較長序列長度/256);

for i=1 to較短序列長度{

cacl_Fij<<>>;//并行計算方法Fij

并行前綴掃描算法計算Eij;}

回溯算法尋找最長相似子序列;

輸出比對結(jié)果;

end

算法1中每次計算Fij和Eij時將重新啟動新的內(nèi)核，因為在GPU中要實現(xiàn)所有線程在某個代碼處同步，只能通過在該位置結(jié)束內(nèi)核的運行，并啟動新的內(nèi)核執(zhí)行下一步代碼.

1.3 基于BLOCK分塊的并行前綴掃描算法

算法1中E值的計算是一種典型的并行掃描算法，這種算法僅需要O(logn)次計算就能找出所有元素的前綴最大值.但隨著掃描的進行，同一個線程束內(nèi)的線程所訪問的元素將無法保持在同一個128byte段內(nèi)，即原來的并行訪問將變成n次串行訪問.并且為了保證線程間數(shù)據(jù)的同步，掃描的每一步都要啟動一次內(nèi)核，這也將影響掃描算法的效率.

本文提出了一種基于BLOCK分塊的并行前綴掃描算法，將每行的Eij以256為單位劃分成1個BLOCK，不足256的加0補足到256，計算時將相應(yīng)BLOCK的Eij元素加載到Share Memory中，分別經(jīng)過UP-DOWN步驟、BLOCK間調(diào)整、Eij微調(diào)3個步驟最終計算出Eij.算法中以256作為劃分單位的原因:一是每個BLOCK內(nèi)256個線程是比較適中的值，同時又解決了合并訪問的問題；二是減少了計算的復(fù)雜度，在塊內(nèi)可以使用移位指令來加速算法的執(zhí)行.Eij的計算分別經(jīng)歷了3個步驟,這3個步驟的具體執(zhí)行過程分別描述如下.

1.3.1 UP-DOWN步驟 采用UP-DOWN步驟計算出每個BLOCK內(nèi)Eij元素的前綴最大值.圖4、5以8個Eij元素劃分為1個BLOCK，分析了其算法過程，圖中虛線部分為直接賦值，實線部分為取兩者的較大值.

圖4 UP步驟

圖5 DOWN步驟

算法2：UP步驟

begin

tx=塊內(nèi)線程號;

將Eij從全局內(nèi)存加載到共享內(nèi)存

shareEij中;

塊內(nèi)線程并行 ford=0 to

{

if(tx<(256>>d+1))

{

e1=tx*(1<<(d+1))+(1<

e2=(tx+1)*(1<<(d+1))-1;

if(shareEij[e2]

shareEij[e2]=shareEij[e1];

}

__syncthreads();//塊內(nèi)同步}

End

算法3：DOWN步驟

begin

將每個BLOCK內(nèi)最后一個元素的Eij值存儲到

BLOCKEij變量中，以便進行塊間調(diào)整;

shareEij[255]=極小值;

塊內(nèi)線程并行：ford=7 to 0 {

if(tx<(1<<7-d)) {

e1=tx*(1<<(d+1))+(1<

e2=(tx+1)*(1<<(d+1))-1;

temp=shareEij[e1];

shareEij[e1]=shareEij[e2];

if(shareEij[e2]

shareEij[e2]=temp; }

__syncthreads();}

將共享內(nèi)存shareEij寫回Eij;

end

圖6 每個分塊的前綴最大值調(diào)整Fig.6 Adjustment of prefix max value of each block

1.3.2 BLOCK間調(diào)整 采用UP-DOWN步驟得到的前綴最大值，是基于每個BLOCK的，為了得到全局的前綴最大值，還必須從第2個BLOCK開始對Eij值進行調(diào)整.即將第i個BLOCK的每個Eij與第i-1個BLOCK的最后一個元素進行比較，取其較大值作為新的Eij值，其偽代碼見算法4，運行過程如圖6所示.

算法4：調(diào)整BLOCK間Eij

begin

256個線程并行執(zhí)行：fori=2 to最后一個BLOCK{

Eij[i*256+tx]=max(blockEij[i-1],

Eij[i*256+tx],Eij[i*256-1]);

__syncthreads();

}

end

1.3.3Eij微調(diào) 經(jīng)過前述步驟計算得到的Eij多計算了n-i個空位得分Ge,所以將Eij加上(n-i)*Ge得到最終的Eij值，具體過程如算法5所描述.

算法5：Eij微調(diào)

begin

for 并行：{

計算線程號tid

Eij[tid]=Eij[tid]+(n-tid-1)*Ge-Gs;

score[i*(n+1)+tid]

=Max(Hij[tid],Eij[tid])；

Eij[tid]=0;

}

end

1.4 Smith-Waterman回溯算法的改進

回溯算法需要先掃描整個DP矩陣以尋找最大值，算法中得到的DP矩陣存儲在GPU端，直接在GPU端并行查找最大值，避免了GPU與CPU之間的內(nèi)存?zhèn)鬏?同時在回溯過程中并行查找最大值，然后啟動一線程從最大值開始回溯.

2 結(jié)果與分析

2.1 各個算法的運行時間

硬件環(huán)境如下：i3-2310處理器，8G內(nèi)存，GPU分別為GT520M和GTX650TI BOOST.在相同的試驗環(huán)境下，Khajeh-saeed et al[8]提出的GPU并行算法(SWAT)的運行時間與本文改進的并行算法相比較，選取的序列長度分別為200、103、104、105、106.從表1可以看出，本文提出的算法的性能大大優(yōu)于SW算法，運算時間只有SWAT算法的十分之一，而在GTX650TI BOOST上運行的效率則進一步提高.

表1 各算法的運行時間

2.2 各算法CUPS性能的比較

在序列比對中，一般采用CUPS衡量算法性能，CUPS越高表示算法的性能越好，計算公式為(M×N)/T，M和N分別代表2個序列的長度，T表示比對時間.圖7比較了傳統(tǒng)的Smith-Waterman算法、SWAT算法及本文提出的算法在GT520及GTX650TI BOOST上的CUPS.從圖7可以看出，在GTX650平臺上本文的算法最優(yōu)，且在GT520平臺上也展示了較好的性能，均優(yōu)于SWAT算法在GT520和GTX650平臺上的結(jié)果；而SW算法性能最差.

圖7 各算法的CUPS與M×N的關(guān)系

2.3 各種算法的加速比分析

為了進一步比較算法的性能，比較了SWAT算法與本文算法的加速比.加速比公式為TACC=TCPU/TGPU，TCPU代表在CPU平臺上的運行時間，TGPU表示在GPU上的運行時間.圖8展示了序列長度和加速比的關(guān)系，可以看到本文算法的加速比隨著序列的增長迅速增長，在GT520平臺上達到90倍左右后逐漸達到飽和.當(dāng)序列長度達到一定時，加速比逐漸達到緩和，而在GTX650平臺上則加速比隨著序列的增長而增加，最多達到250倍左右；而SWAT算法在GT520和GTX650平臺上分別最多只能達到10倍和50倍.

圖8 各算法的序列長度與加速比的關(guān)系

3 小結(jié)

本文在GPU-CPU異構(gòu)體上提出了GPU并行的Smith-Waterman算法，通過改進算法求解公式，采用按列并行的方式，并且提出基于BLOCK的并行前綴掃描算法.該并行前綴掃描算法經(jīng)過3個步驟得到E值.這3個步驟分別是UP-DOWN步驟、BLOCK間調(diào)整、Eij微調(diào).算法在GT520和GTX650TI BOOST上進行運行，序列長度分別為200、103、104、105、106；并在此環(huán)境下同時比較了傳統(tǒng)的Smith-Waterman算法、SWAT算法和本文算法.

結(jié)果表明：本文算法在保證原有算法敏感度的基礎(chǔ)上，在低端顯卡GT520M上最快達到90倍的加速比；而在GTX650TI BOOST中端顯卡上最快有262倍的加速比，而同樣基于GPU并行的SWAT算法[8]最多只有100倍的加速比，大大優(yōu)于傳統(tǒng)的Smith-Waterman算法.從各算法的CUPS性能指標(biāo)上看，本文算法的CUPS增長最快，其次是SWAT算法；并且隨著序列長度的增加，本文算法的CUPS增長性能隨之加快，而SWAT算法增長緩慢.從各算法的具體執(zhí)行時間來看，采用本文算法比對較短序列的性能，與傳統(tǒng)的Smith-Waterman算法比較接近，略優(yōu)于SWAT算法；但在比對較長序列時，本文算法的運行時間大大低于其他2種算法.

[1] 陳華,蔡鐵城,張沖,等.花生葉全長cDNA文庫的構(gòu)建和部分表達序列標(biāo)簽分析[J].福建農(nóng)林大學(xué)學(xué)報:自然科學(xué)版,2014,43(6):596-601.

[2] ALTSCHUL S F, GISH W, MILLER W, et al. Basic local alignment search tool[J]. Journal of Molecular Biology, 1990,215(3):403-410.

[3] SMITH T F, WATERMAN M S. Identification of common molecular subsequences[J]. Journal of Molecular Biology, 1981,147(1):195-197.

[4] ZHAO X, LI H, BAO T. Analysis on the interaction between post-spliced introns and corresponding protein coding sequences in ribosomal protein genes[J]. Journal of Theoretical Biology, 2013,328:33-42.

[5] WOZNIAK A. Using video-oriented instructions to speed up sequence comparison[J]. Computer Applications in the Biosciences: CABIOS, 1997,13(2):145-150.

[6] FARRAR M. Striped Smith-Waterman speeds database searches six times over other SIMD implementations[J]. Bioinformatics, 2007,23(2):156-161.

[7] AKOGLU A, STRIEMER G M. Scalable and highly parallel implementation of Smith-Waterman on graphics processing unit using CUDA[J]. Cluster Computing, 2009,12(3):341-352.

[8] KHAJEH-SAEED A, POOLE S, BLAIR PEROT J. Acceleration of the Smith-Waterman algorithm using single and multiple graphics processors[J]. Journal of Computational Physics, 2010,229(11):4247-4258.

[9] HARRIS M, SENGUPTA S, OWENS J D. Parallel prefix sum (scan) with CUDA[J]. GPU Gems, 2007,3(39):851-876.

[10] ZHANG Y, OWENS J D, SENGUPTA S, et al. Scan primitives for GPU computing [J]. Graphics Hardware, 2007(3):97-106.

[11] 馬海晨,韋剛,吳百峰.基于GPGPU的生物序列快速比對[J].計算機工程,2012,38(4):241-244.

[12] LIU W, SCHMIDT B, VOSS G, et al. Streaming algorithms for biological sequence alignment on GPUs[J]. IEEE Transactions on Parallel & Distributed Systems, 2007, 18(9):1270-1281.

(責(zé)任編輯:葉濟蓉)

Parallel Smith-Waterman algorithm based on GPU

LIN Min, ZHONG Yi-wen, LIN Juan

(College of Computer and Information, Fujian Agriculture and Forestry University, Fuzhou, Fujian 350002, China)

The formulae of Smith-Waterman algorithm was improved to make it adapt to the parallel characteristics of GPU, and a novel strategy of parallel prefix scan was proposed. The prefix maximum of each element in the row within O (logn) time was caculated through UP-DOWN steps and the adjustment between blocks andEijfine tuning. At last, the backtrack procudure ran at GPU side to avoid the memory copies between GPU and CPU. Compared with traditional Smith-Waterman algorithm, this algorithm performance increased 90 times on the lower-end GPU platform, and also had larger ascension in comparison with SWAT algorithm.

Smith-Waterman algorithm; parallel prefix scan; graphics processing unit; sequences alignment

2014-08-28

2015-02-15

福建省自然科學(xué)基金資助項目(2013J01216、2014J01219).

林敏(1981-),男,講師,碩士.研究方向:智能信息處理、并行計算及生物信息學(xué).Email:jsjxy_lm@126.com.通訊作者鐘一文(1968-),男, 教授,博士.研究方向:智能計算.Email:yiwzhong@fafu.edu.cn.

TP391

A

1671-5470(2015)04-0442-07

10.13323/j.cnki.j.fafu(nat.sci.).2015.04.019