CUDA加速的DNA-蛋白質(zhì)匹配及其優(yōu)化

陳春雷,慕德俊,張慧翔,胡 偉

(西北工業(yè)大學(xué) 自動(dòng)化學(xué)院,陜西 西安 710072)

結(jié)構(gòu)生物信息學(xué)已成為當(dāng)前計(jì)算機(jī)科學(xué)領(lǐng)域的一個(gè)研究熱點(diǎn)。其主要研究目標(biāo)是獲取生物系統(tǒng)的高分辨率結(jié)構(gòu)信息，從而精確地推理系統(tǒng)的功能、預(yù)測(cè)某些修改或擾動(dòng)對(duì)系統(tǒng)造成的影響。與生物信息學(xué)的其他領(lǐng)域相比，結(jié)構(gòu)生物信息學(xué)面臨眾多新的挑戰(zhàn)。其中之一是多數(shù)結(jié)構(gòu)生物信息學(xué)問(wèn)題的搜索空間是連續(xù)的；換言之，生物系統(tǒng)的微觀結(jié)構(gòu)一般通過(guò)原子在空間坐標(biāo)系中的位置來(lái)表示，而坐標(biāo)通常是連續(xù)變量,因此搜索空間是無(wú)窮的。雖然某些近似方法可以在一定程度上應(yīng)對(duì)結(jié)構(gòu)生物信息學(xué)問(wèn)題的這種連續(xù)特性，但是從搜索空間中求得最終解仍然需要進(jìn)行大量的運(yùn)算[1]。

DNA-蛋白質(zhì)匹配是生物信息學(xué)的一個(gè)重要問(wèn)題。傳統(tǒng)的序列分析方法往往會(huì)導(dǎo)致較高的多重檢驗(yàn)錯(cuò)誤率 (false-positive rate)[2]。多數(shù)基于結(jié)構(gòu)生物信息學(xué)的DNA-蛋白質(zhì)匹配方法把獲得轉(zhuǎn)錄因子-DNA組合體的信息作為一個(gè)必備條件，而轉(zhuǎn)錄因子-DNA組合體的信息通常需要通過(guò)實(shí)驗(yàn)獲得，這種實(shí)驗(yàn)通常需要耗費(fèi)大量的時(shí)間。參考文獻(xiàn)[2]提出的基于退火算法的匹配算法避開(kāi)了這個(gè)問(wèn)題,但是使用該方法進(jìn)行DNA-蛋白質(zhì)匹配仍需要較大的運(yùn)算量。參考文獻(xiàn)[3]基于參考集索引技術(shù)提出了一種快速的序列分析算法，并將其應(yīng)用于DNA-蛋白質(zhì)匹配。參考文獻(xiàn)[4]采用多CPU的方式設(shè)計(jì)實(shí)現(xiàn)了并行化的DNA-蛋白質(zhì)序列分析方法，并取得了較高的加速比。但上述工作均未涉及基于結(jié)構(gòu)生物信息學(xué)的匹配方法的加速問(wèn)題。

應(yīng)對(duì)大運(yùn)算量問(wèn)題的一個(gè)常用方法是并行計(jì)算。NVIDIA公司統(tǒng)一計(jì)算設(shè)備架構(gòu)CUDA(Compute Unified Device Architecture)是一種全新的并行計(jì)算架構(gòu)。在CUDA的支持下，通用計(jì)算圖形處理單元GPGPU(General Purpose Graphic Processing Unit)可以作為一種通用的效率、硬件加速器，發(fā)揮出強(qiáng)大的運(yùn)算能力[5]。

本文針對(duì)參考文獻(xiàn)[2]的DNA-蛋白質(zhì)匹配方法，從計(jì)算的角度對(duì)其進(jìn)行分析，設(shè)計(jì)實(shí)現(xiàn)并優(yōu)化了基于CUDA的加速方法，通過(guò)實(shí)驗(yàn)驗(yàn)證了其加速性能。

1 相關(guān)知識(shí)介紹

1.1 基于退火算法的DNA-蛋白質(zhì)匹配算法

給定一條DNA鏈和一條蛋白質(zhì)鏈,兩者之間的相對(duì)位置存在很多種可能。假定兩者均為剛體，DNA-蛋白質(zhì)組合體的物理能量主要受兩者相對(duì)位置的影響，組合體的物理能量越低，其穩(wěn)定性越強(qiáng)。該方法試圖尋找使得組合體最穩(wěn)定的相對(duì)位置，這一結(jié)果對(duì)生物信息學(xué)的研究具有重要意義。

從計(jì)算的角度看，該方法的基礎(chǔ)是退火算法。為了提高搜索精度，該匹配方法通常需要隨機(jī)產(chǎn)生大量初始“種子”，為每個(gè)“種子”啟動(dòng)一個(gè)基于退火算法的搜索過(guò)程，在整個(gè)解空間中尋找最穩(wěn)定的相對(duì)位置。這種匹配方法的流程如圖1所示。

圖1 基于退火算法的DNA-蛋白質(zhì)匹配方法流程圖

圖1中的初始值為T0，溫度T及溫度衰減周期interval為正整數(shù),溫度衰減系數(shù)為a。每當(dāng)算法所經(jīng)歷的平移和旋轉(zhuǎn)次數(shù)達(dá)到interval的非零整數(shù)倍時(shí)，T衰減為原來(lái)的a倍(0＜a＜1)，Tmin為溫度的最小值。

DNA-蛋白質(zhì)組合體的物理能量E1、E2分別代表第n次(1≤n≤Smax，Smax為平移和旋轉(zhuǎn)次數(shù)的最大值)平移和旋轉(zhuǎn)前后組合體的物理能量。Emin代表通過(guò)算法得到的組合體物理能量的最小值。如果E2＜E1，則接受第n次平移和旋轉(zhuǎn)的結(jié)果；否則,計(jì)算指數(shù)式exp(-(E2-E1)/T)的值，并使用C++庫(kù)函數(shù)drand48()產(chǎn)生一個(gè)在區(qū)間[0,1]上均勻分布的隨機(jī)數(shù)。如果指數(shù)式的值小于隨機(jī)數(shù)的值，則接受第n次平移旋轉(zhuǎn)的結(jié)果，反之不接受。

step為當(dāng)前進(jìn)行到了第幾次平移和旋轉(zhuǎn)。

1.2 CUDA編程模型及線程調(diào)度方式的特點(diǎn)

CUDA程序通常包含CPU串行代碼和GPU并行代碼，后者被稱作CUDA程序的內(nèi)核。在執(zhí)行內(nèi)核時(shí)，CUDA會(huì)產(chǎn)生大量并行工作的線程，以單指令多數(shù)據(jù)SIMD(Single Instruction Multiple Data)方式完成整個(gè)內(nèi)核的計(jì)算任務(wù)。CUDA采用網(wǎng)格(grid)和線程塊(block)來(lái)組織線程[6]。一個(gè)block的線程又被劃分為一個(gè)或多個(gè)線程組(warp)。warp是CUDA程序最基本的調(diào)度單位，屬于同一個(gè)warp的各個(gè)線程會(huì)從CUDA程序代碼段中相同的程序地址同時(shí)開(kāi)始執(zhí)行，但是各線程具有獨(dú)立的當(dāng)前指令指針和寄存器狀態(tài)。因此，各線程可能會(huì)有不同的執(zhí)行路徑，例如執(zhí)行不同的分支選擇結(jié)構(gòu)代碼或循環(huán)結(jié)構(gòu)代碼[7]。warp執(zhí)行特點(diǎn)如圖2所示。

圖2 warp的執(zhí)行特點(diǎn)示意圖

假設(shè)一個(gè)warp中共有4個(gè)線程：線程1～線程4，它們的執(zhí)行時(shí)間分別是t1～t4，并且t3＜t1＜t2＜t4。因?yàn)镃UDA的基本調(diào)度單位是一個(gè)完整的warp而非單個(gè)線程，所以整個(gè)warp的執(zhí)行時(shí)間取決于執(zhí)行時(shí)間最長(zhǎng)的線程t4。其他3個(gè)線程在執(zhí)行完畢后必須等待還沒(méi)有執(zhí)行完畢的線程，因此有一段時(shí)間處于空閑狀態(tài)，相應(yīng)的計(jì)算資源也就被閑置了。例如，線程3閑置的計(jì)算資源最多，其空閑時(shí)長(zhǎng)為(t4-t3)。造成計(jì)算資源閑置的原因是同一warp中各個(gè)線程的執(zhí)行路徑不同；而CUDA采用的是SIMD并行方式，執(zhí)行路徑的不同是由每個(gè)線程所處理的數(shù)據(jù)差異造成的。因此，依照一定規(guī)則對(duì)數(shù)據(jù)進(jìn)行重新排列是減少這種計(jì)算資源閑置狀況的常用方法[8]。重排規(guī)則通常視具體應(yīng)用而定。

2 使用CUDA加速DNA-蛋白質(zhì)匹配

2.1 從計(jì)算角度分析匹配算法

匹配方法的流程已由圖1給出，除參數(shù)初始化外，該方法可分為三部分：(1)平移、旋轉(zhuǎn)DNA和蛋白質(zhì)；(2)計(jì)算DNA-蛋白質(zhì)組合體的能量；(3)后續(xù)處理。這3部分在普通CPU平臺(tái)上所耗時(shí)間依次為∶73.624 s、1138.945 s、110.809 s(僅使用一個(gè)隨機(jī)初始種子，退火算法參數(shù)的預(yù)設(shè)值及軟硬件平臺(tái)參數(shù)指標(biāo)見(jiàn)本文第3部分),實(shí)驗(yàn)使用的DNA、蛋白質(zhì)結(jié)構(gòu)數(shù)據(jù)來(lái)自參考文獻(xiàn)[9]編號(hào)為2GLI的文件。其他的DNA、蛋白質(zhì)結(jié)構(gòu)數(shù)據(jù)的實(shí)驗(yàn)結(jié)果也顯示了計(jì)算能量所耗費(fèi)的時(shí)間遠(yuǎn)超過(guò)其他兩個(gè)部分。計(jì)算能量時(shí)，DNA-蛋白質(zhì)組合體所處的三維空間被劃分成若干個(gè)棱長(zhǎng)為埃米(10-10m)級(jí)的晶格。一條DNA鏈由若干個(gè)DNA原子構(gòu)成，一條蛋白質(zhì)鏈由若干個(gè)蛋白質(zhì)原子構(gòu)成；以2GLI組合體為例,共有855個(gè)DNA原子和1 270個(gè)蛋白質(zhì)原子。每次平移和旋轉(zhuǎn)前后，每個(gè)原子所屬的晶格都可能發(fā)生改變，每個(gè)晶格所包含的原子個(gè)數(shù)也可能發(fā)生改變。Di(i=1,2,…,ND)，Pj(j=1,2,…,NP)分別代表組合體中的DNA原子數(shù)量和蛋白質(zhì)原子數(shù)量。DNA原子i(i=1,2,…,ND)在三維空間中所處的晶格和相鄰的晶格共有27個(gè)，統(tǒng)稱為原子i的臨近晶格(Neighboring Lattice)，以pi,js(x,1,2,…,ni)代表這27個(gè)晶格中的所有蛋白質(zhì)原子，這些原子即原子i的相鄰蛋白質(zhì)原子。以di,jx代表DNA原子i與蛋白質(zhì)原子jx之間的歐式距離，兩個(gè)原子之間的能量是di,jx的函數(shù)：E(di,jx)。則組合體的總能量為：

根據(jù)以上分析可知：(1)計(jì)算能量是匹配方法中最耗時(shí)部分；(2)E(di,jx)(i=1,2,…,ND;x=1,2,…,ni)的計(jì)算是相互獨(dú)立的，無(wú)依賴關(guān)系，適于并行。根據(jù)阿姆達(dá)爾定律[10]，選擇計(jì)算能量部分進(jìn)行并行化。

2.2 基于CUDA的并行匹配算法

以下偽代碼給出了串行方式計(jì)算能量的步驟：

算法1：串行方式計(jì)算能量

(1)E←0

(2)for i=1 to ND

(3) for all neighboring lattice of atomi∶la=1 to 27

(4) for all protein atoms in lattice la

(5) compute E(di,jx)and E←E+E(di,jx)

(6) end for

(7) end for

(8)end for

令CUDA block為一維結(jié)構(gòu)，一個(gè)block中包含B個(gè)線程，一個(gè)初始種子的搜索過(guò)程僅由一個(gè)block完成?；谝陨洗写a，給出CUDA內(nèi)核的偽代碼：

算法2：并行方式計(jì)算能量

(1)tid←threadIdx.x

(2)Ethis_thread←0,stride←「ND/B?

(3)for s=0 to stride

(4)i←s.B+tid

(5)if i＜ND

(6) for all neighboring lattice of atom i∶la=1 to 27

(7) for all protein atoms in lattice la

(8) compute and E(di,jx)and Ethis_thread←Ethis_thread+E(di,jx)

(9) end for

(10) end for

(11)end for

算法2實(shí)現(xiàn)了算法1第2行循環(huán)操作的并行化。每個(gè)CUDA線程對(duì)應(yīng)一個(gè)DNA原子，依次進(jìn)入該原子的27個(gè)相鄰晶格，找到各個(gè)相鄰蛋白質(zhì)原子,計(jì)算E(di,jx)，并累加到能量的部分和：Ethis_thread。在CUDA內(nèi)核返回后，通過(guò)累加各線程的能量部分和可以得到組合體的總能量。另外，如果B＜ND，則需要把ND個(gè)DNA原子平均分成「ND/B?塊（最后一塊可能不足B個(gè)原子），然后，由整個(gè)block的線程依次處理各塊，算法2第3行的循環(huán)結(jié)構(gòu)即為了達(dá)到這個(gè)目的。

2.3 優(yōu)化并行算法

考慮到本文1.2節(jié)所述warp的執(zhí)行特點(diǎn)，上述并行方式存在一定的不足。算法2中第7行的循環(huán)次數(shù)取決于當(dāng)前晶格中蛋白質(zhì)原子的個(gè)數(shù)，同一warp中各線程的循環(huán)次數(shù)可能會(huì)不同，如果差異過(guò)大，會(huì)造成嚴(yán)重的計(jì)算資源閑置；對(duì)所有線程而言，算法2第3行、第6行循環(huán)結(jié)構(gòu)的循環(huán)次數(shù)是確定的。

假設(shè)DNA存儲(chǔ)在一個(gè)數(shù)組中，且該數(shù)組位于GPU主存(global memory)中。為了盡量減少計(jì)算資源閑置，重排DNA原子的原有次序(DNA鏈中的次序),處于同一個(gè)晶格中的DNA原子被存儲(chǔ)在數(shù)組的某一段連續(xù)區(qū)域內(nèi)。采用這種優(yōu)化措施是基于以下原因：(1)由式(1)可知，總能量由E(di,jx)做算術(shù)加法得到，與E(di,jx)的計(jì)算次序無(wú)關(guān)；(2)為了提高GPU主存的帶寬利用率，同一warp中的線程通常從主存中的臨近區(qū)域讀取數(shù)據(jù)，即內(nèi)存訪問(wèn)對(duì)齊 (coalescing)[6]；(3)DNA的雙螺旋結(jié)構(gòu)是非線性的，DNA在鏈中相鄰的原子未必處于同一晶格中;(4)同一晶格中的DNA原子的臨近蛋白質(zhì)原子數(shù)量相同，重排可以減少因線程執(zhí)行路徑差異造成的計(jì)算資源閑置。重排的效果如圖3所示。

圖3 DNA原子存儲(chǔ)次序重排示意圖

3 實(shí)驗(yàn)

3.1 實(shí)驗(yàn)平臺(tái)

硬 件 平 臺(tái)：Intel?CoreTM2 E7500 CPU，2.93 GHz，NVIDIA GTX 660圖形處理器(單個(gè)warp包含32個(gè)線程)，CPU主存4 GB。

軟件平臺(tái)：Linux操作系統(tǒng)內(nèi)核版本2.6.18-194.el5，gcc版本4.1.2，CUDA版本5.0。

3.2 實(shí)驗(yàn)參數(shù)設(shè)置與結(jié)果

DNA-蛋白質(zhì)結(jié)構(gòu)數(shù)據(jù)編號(hào)為2GLI。CPU版程序以單線程串行方式執(zhí)行；GPU版程序block size固定為512；各block從不同的初始種子出發(fā)，分別基于退火算法進(jìn)行搜索，每個(gè)block對(duì)應(yīng)一個(gè)初始種子。退火算法參數(shù)：初始溫度T0=120℃,最低溫度Tmin=0.001℃,溫度衰減周期interval=800,溫度衰減系數(shù)a=0.99,最大平移旋轉(zhuǎn)次數(shù)Smax=106。

實(shí)驗(yàn)結(jié)果如表1所示。與單線程CPU程序相比，未經(jīng)優(yōu)化的GPU程序?qū)@得最高可達(dá)8倍左右的加速比；而經(jīng)過(guò)重排優(yōu)化后，加速比在此基礎(chǔ)上進(jìn)一步顯著提升，最高可達(dá)15倍左右。

表1 實(shí)驗(yàn)結(jié)果

本文針對(duì)一種典型的DNA-蛋白質(zhì)匹配算法，設(shè)計(jì)實(shí)現(xiàn)了基于CUDA的并行化方法，從線程調(diào)度的角度對(duì)該方法進(jìn)行優(yōu)化，并通過(guò)實(shí)驗(yàn)驗(yàn)證了加速性能。實(shí)際應(yīng)用往往需要為一個(gè)組合體產(chǎn)生大量的初始種子(數(shù)百個(gè))，并為每個(gè)種子開(kāi)啟一個(gè)基于退火算法的搜索過(guò)程；其目的是達(dá)到較高的匹配精度。實(shí)驗(yàn)顯示，使用單個(gè)GPGPU加速，在單個(gè)block包含的線程數(shù)不變的前提下，隨著初始種子數(shù)量的增加，加速比逐漸趨于穩(wěn)定。例如，當(dāng)初始種子個(gè)數(shù)超過(guò)40后，加速比基本穩(wěn)定在15倍左右。其原因在于單個(gè)GPGPU的計(jì)算能力存在上限,當(dāng)種子足夠多時(shí)，其計(jì)算能力已得到較充分利用，無(wú)法繼續(xù)提高加速比。為了滿足實(shí)際應(yīng)用的需求，下一步的工作將考慮使用基于GPGPU的集群來(lái)加速匹配算法，以進(jìn)一步提高加速比。

[1]BOURNE P E，WEISSIG H.Structural bioinformatics[M].Hoboken∶Wiley-Liss Inc,2003.

[2]Liu Zhijie，Guo Juntao，Li Ting,et al.Structure-based prediction of transcription factor binding sites using a protein-DNA docking approach[J].Proteins,2008,72(4)∶1114-1124.

[3]戴東波，熊赟，朱揚(yáng)勇.基于參考集索引的高效序列相似性查找算法[J].軟件學(xué)報(bào)，2011(4)∶718-731.

[4]Zhang Zhang,Xiao Jingfa,Wu Jiayan,et al.ParaAT∶A parallel tool for constructing multiple protein-coding DNA alignments[J].Biochemical Biophysical Research Communications,2012,419(4)∶779-781.

[5]徐新海，楊學(xué)軍，林宇斐，等.一種面向CPU-GPU異構(gòu)系統(tǒng)的容錯(cuò)方法[J].軟件學(xué)報(bào)，2011，22(10)：2538-2552.

[6]NVIDIA Cooperation.CUDA programming guide version 5.0[EB/OL].[2013-05-15].http://docs.nvidia.com/cuda/cudac-programming-guide/.

[7]TIANYI D H，TAREK S A.Reducing branch divergence in GPU programs[A].In Proceedings of the Fourth Workshop on General Purpose Processing on Graphics Processing Units[C].London∶ACM.2012.

[8]IMEN C,AHCENE B,NOUREDINE M.Reducing thread divergence in GPU-based b&b applied to the flow-shop problem[A].In Proceedings of the 9th International Conference on Parallel[C].Berlin∶Springer-Verlag.2011.

[9]Rutgers and UCSD.Protein Data Bank[DB/OL].[2009-02-24].http://www.rcsb.org/pdb/explore/explore.do?structureId=2gli.

[10]GENE A.Validity of the single processor approach to achieving large-scale computing capabilities[A].In Proceedings of the April 18-20(AFIPS′67)[C].New York∶ACM.1967.