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

林 敏, 鐘一文, 林 娟

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

?

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

林 敏, 鐘一文, 林 娟

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

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

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

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

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

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

1 GPU并行的Smith-Waterman算法

1.1 Smith-Waterman計(jì)算公式的改進(jìn)

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

(1)

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

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

圖2 回溯的過(guò)程

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

通過(guò)式(2)、式(3)先求每一列的最大值，采用式(4)計(jì)算行最大值,最后通過(guò)式(5)計(jì)算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為比對(duì)得分，Gs和Ge分別表示空位罰分和空位擴(kuò)展罰分.本文中的算法是根據(jù)GPU全局內(nèi)存合并訪問(wèn)的特點(diǎn)，采用了列并行的方式.同時(shí)為了采用GPU并行前綴掃描算法計(jì)算Eij，將式(3)、式(4)重新設(shè)計(jì)成式(6)、式(7)，無(wú)需再保存Uij變量，需要的存儲(chǔ)空間更少.

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

(6)

(7)

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

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

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

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

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

算法1：主程序

begin

載入序列，將序列二值化壓縮存儲(chǔ);

在GPU端為Hij,Fij,Eij分配存儲(chǔ)空間;

計(jì)算blocknum=ceil(較長(zhǎng)序列長(zhǎng)度/256);

for i=1 to較短序列長(zhǎng)度{

cacl_Fij<<>>;//并行計(jì)算方法Fij

并行前綴掃描算法計(jì)算Eij;}

回溯算法尋找最長(zhǎng)相似子序列;

輸出比對(duì)結(jié)果;

end

算法1中每次計(jì)算Fij和Eij時(shí)將重新啟動(dòng)新的內(nèi)核，因?yàn)樵贕PU中要實(shí)現(xiàn)所有線程在某個(gè)代碼處同步，只能通過(guò)在該位置結(jié)束內(nèi)核的運(yùn)行，并啟動(dòng)新的內(nèi)核執(zhí)行下一步代碼.

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

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

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

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

圖4 UP步驟

圖5 DOWN步驟

算法2：UP步驟

begin

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

將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

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

BLOCKEij變量中，以便進(jìn)行塊間調(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寫(xiě)回Eij;

end

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

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

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

begin

256個(gè)線程并行執(zhí)行：fori=2 to最后一個(gè)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)過(guò)前述步驟計(jì)算得到的Eij多計(jì)算了n-i個(gè)空位得分Ge,所以將Eij加上(n-i)*Ge得到最終的Eij值，具體過(guò)程如算法5所描述.

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

begin

for 并行：{

計(jì)算線程號(hào)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回溯算法的改進(jìn)

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

2 結(jié)果與分析

2.1 各個(gè)算法的運(yùn)行時(shí)間

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

表1 各算法的運(yùn)行時(shí)間

2.2 各算法CUPS性能的比較

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

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

2.3 各種算法的加速比分析

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

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

3 小結(jié)

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

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

[1] 陳華,蔡鐵城,張沖,等.花生葉全長(zhǎng)cDNA文庫(kù)的構(gòu)建和部分表達(dá)序列標(biāo)簽分析[J].福建農(nóng)林大學(xué)學(xué)報(bào):自然科學(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的生物序列快速比對(duì)[J].計(jì)算機(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é)任編輯:葉濟(jì)蓉)

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é)基金資助項(xiàng)目(2013J01216、2014J01219).

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

TP391

A

1671-5470(2015)04-0442-07

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