陳 飛，郝福珍

（華北計(jì)算技術(shù)研究所，北京100083）

0 引 言

研究RNA二級結(jié)構(gòu)對于防治以RNA進(jìn)行基因傳遞的疾病有重要作用。長鏈RNA可能包含幾千到幾百萬不等數(shù)量的核苷酸，用X射線晶體衍射、核磁共振等方法測定其二級結(jié)構(gòu)，對時(shí)間和資金消耗都很大，當(dāng)前主流的做法是使用基于最小自由能的動(dòng)態(tài)規(guī)劃模型預(yù)測RNA二級結(jié)構(gòu)。在1981年Zuker和Stiegler提出動(dòng)態(tài)規(guī)劃模型求解該問題［1］后，一方面人們提出了啟發(fā)式算法和近似算法來優(yōu)化動(dòng)態(tài)規(guī)劃模型，減小了計(jì)算壓力；另一方面人們還提出了更多的分子熱力學(xué)細(xì)節(jié)以及相應(yīng)的算法，來更準(zhǔn)確地預(yù)測RNA二級結(jié)構(gòu)，但這樣也會(huì)相應(yīng)的增加整體算法的時(shí)間復(fù)雜度和空間復(fù)雜度。針對RNA鏈長度增加和預(yù)測準(zhǔn)確度增高帶來的雙重計(jì)算壓力的難題，并行算法是一種值得嘗試的解決辦法。

本課題選用的并行計(jì)算平臺(tái)是nVIDIA公司推出的基于CUDA編程模型的可編程GPU，型號為GTX550，它擁有384個(gè)處理器核心。相對于CPU，GPU擁有更多的可用計(jì)算核心，在浮點(diǎn)運(yùn)算和并行計(jì)算方面擁有更明顯的優(yōu)勢，同時(shí)GPU比小型機(jī)或計(jì)算機(jī)集群而言更容易獲取。CUDA編程模型提供的軟件環(huán)境使編寫程序、操作GPU就如同控制大規(guī)模并行機(jī)器一樣便捷，不需要像以前依賴圖形API接口來實(shí)現(xiàn)對GPU的訪問，使得在GPU上實(shí)現(xiàn)并行通用計(jì)算程序的難度大大降低。同時(shí)CUDA提供了內(nèi)存管理、高性能計(jì)算指令、設(shè)備訪問和執(zhí)行調(diào)度等函數(shù)，方便控制數(shù)據(jù)移動(dòng)及優(yōu)化程序。

文章中討論的并行算法包含了一系列依賴緊密的數(shù)據(jù)，包括5個(gè)大小為n2（n表示RNA序列長度）的數(shù)組以及一系列索引數(shù)組。在解決該問題中，合理利用共享存儲(chǔ)器（shared memory）、常量存儲(chǔ)器 （constant memory）等存儲(chǔ)數(shù)據(jù)，同時(shí)控制設(shè)備內(nèi)存中的數(shù)據(jù)進(jìn)行移動(dòng)，可以有效的加速算法的運(yùn)行。

1 RNA二級結(jié)構(gòu)預(yù)測算法及動(dòng)態(tài)規(guī)劃求解

RNA是由核苷酸組合成的聚合物，其包含A、C、G、U這4種堿基，這4種堿基由A－U、C－G和U－G這3種組合配對，并由氫鍵連接而成的結(jié)構(gòu)被稱為堿基對，配對全集 ｛（A－U），（U－A），（C－G），（G－C），（U－G），（G－U）｝。預(yù)測RNA二級結(jié)構(gòu)是指，取給定RNA堿基序列中的各個(gè)位點(diǎn)進(jìn)行配對，將得到的配對結(jié)果轉(zhuǎn)換成二級模型的過程。其中使用計(jì)算機(jī)用動(dòng)態(tài)規(guī)劃算法得到的結(jié)果不包括偽結(jié)。

1.1 RNA二級結(jié)構(gòu)預(yù)測

圖1 是由79個(gè)堿基組成的RNA的二級結(jié)構(gòu)，在這個(gè)結(jié)構(gòu)中包含發(fā)卡環(huán) （hairpin loop）、內(nèi)部環(huán) （internal loop）、多分支環(huán)（multiloop）、莖（stack）、 外部環(huán)（external loop）等5種子結(jié)構(gòu)。

圖1 RNA二級結(jié)構(gòu)示例

其中莖是由連續(xù)的堿基對堆積構(gòu)成的結(jié)構(gòu)，它在三維空間中的表現(xiàn)為A型雙螺旋，這種子結(jié)構(gòu)對整體的RNA二級結(jié)構(gòu)的穩(wěn)定性起促進(jìn)作用，即令二級結(jié)構(gòu)的自由能降低；除此之外的幾種子結(jié)構(gòu)，比如：莖中間出現(xiàn)一些無法配對的堿基形成的內(nèi)部環(huán)、相鄰連續(xù)序列因兩端互補(bǔ)形成發(fā)卡結(jié)構(gòu)的發(fā)卡環(huán)，會(huì)對整體RNA結(jié)構(gòu)的穩(wěn)定性造成破壞，即令二級結(jié)構(gòu)的自由能升高。

在最小自由能 （minimum free energy）模型中，一個(gè)RNA二級結(jié)構(gòu)的自由能，可以通過將組成它的所有獨(dú)立子結(jié)構(gòu)，如環(huán) （internal loop）、莖 （stack）等所對應(yīng)的自由能累加得到，其中每個(gè)子結(jié)構(gòu)的能量只依賴于環(huán)本身的結(jié)構(gòu)，而和整體結(jié)構(gòu)中其它部分無關(guān)［2］。每個(gè)RNA二級結(jié)構(gòu)都有自身對應(yīng)的自由能，擁有最小自由能的二級結(jié)構(gòu)最穩(wěn)定，被稱為一個(gè)RNA序列的最優(yōu)二級結(jié)構(gòu)。

1.2 動(dòng)態(tài)規(guī)劃求解算法

在1981年，Zuker和Stiegler提出了使用動(dòng)態(tài)規(guī)劃算法求解RNA二級結(jié)構(gòu)預(yù)測問題，并指出預(yù)測給定RNA序列的二級結(jié)構(gòu)問題是一個(gè)類似于Smith－Waterman local alignment算法的最優(yōu)化問題［3］，但更加復(fù)雜。動(dòng)態(tài)規(guī)劃算法作為計(jì)算RNA二級結(jié)構(gòu)的有效算法，其空間復(fù)雜度為O （n2），用V、VBI、VM、WM等至少4個(gè)數(shù)組保存運(yùn)算的中間結(jié)果。中間計(jì)算完成后同傳統(tǒng)動(dòng)態(tài)規(guī)劃一樣，通過回溯找到具有最小自由能的二級結(jié)構(gòu)。

給定長度為N的RNA堿基序列，設(shè)其自由能大小為W（N），i，j作為游標(biāo)，在滿足條件1＜i＜j＜N的情況下，對所有i，j的可行值進(jìn)行計(jì)算。

對于上述長度為N的RNA堿基序列s，其子序列s1…sj的最優(yōu)結(jié)構(gòu)的自由能W（j）為

在上述公式中，V（i，j）表示在給定RNA序列中從第i個(gè)核苷酸至第j個(gè)核苷酸對應(yīng)的子序列擁有的自由能。

如果第i個(gè)核苷酸和第j個(gè)核苷酸能夠配對成堿基對，則

eH（i，j）是序列si…sj形成發(fā)卡環(huán) （hairpin loop）結(jié)構(gòu)所包含的自由能，該發(fā)卡環(huán)序列中第i和第j兩個(gè)核苷酸組成堿基對，并作為子結(jié)構(gòu)整體對外的連接點(diǎn)。

eS（i，j）是由 （i，j）和 （i＋1，j－1）兩個(gè)堿基對組成的莖所包含的自由能。

VBI（i，j）是序列si…sj形成內(nèi)部環(huán) （internal loop）結(jié)構(gòu)所包含的自由能，這個(gè)內(nèi)部環(huán)以序列中第i和第j兩個(gè)核苷酸組成的堿基對作為子結(jié)構(gòu)整體對外的連接點(diǎn)。內(nèi)部環(huán)具體計(jì)算方式如式（3）所示

其中式 （3）還要滿足i′－i＋j－j′－2＞0，公式中eL（i，j）是包含外部堿基對 （i，j）和內(nèi)部堿基對 （i′，j′）的內(nèi)部環(huán)的自由能

VM（i，j）是序列si…sj構(gòu)成擁有最優(yōu)結(jié)構(gòu)的多分支環(huán)包含的自由能。其值大小同單鏈結(jié)構(gòu)的數(shù)量有線性關(guān)系，引入WM來計(jì)算VM，其關(guān)系式如下

用動(dòng)態(tài)規(guī)劃算法求解RNA二級結(jié)構(gòu)問題，其空間復(fù)雜度為O（n2），其中包含V、VBI、VM、WM等4個(gè)占用空間大小為（N＋1）＊（N＋1）／2的數(shù)組，數(shù)組大小與RNA序列長度有關(guān)，動(dòng)態(tài)規(guī)劃前期保存中間結(jié)果以求出的最小自由能 （minimum free energy，MFE），回溯階段作為得到最優(yōu)結(jié)構(gòu)的關(guān)鍵信息。在RNA二級結(jié)構(gòu)預(yù)測中回溯需要的時(shí)間遠(yuǎn)小于求最小自由能的時(shí)間。

使用動(dòng)態(tài)規(guī)劃算法計(jì)算最小自由能的時(shí)間復(fù)雜度為O（n4），預(yù)測一個(gè)給定序列的時(shí)間復(fù)雜度為：O（n2）＊子序列計(jì)算時(shí)間復(fù)雜度。

子序列的最小自由能的計(jì)算方法由上述式 （1）～式（5）給出，其中計(jì)算內(nèi)部環(huán)和多分支環(huán)的時(shí)間復(fù)雜度最大。當(dāng)計(jì)算子序列si…sj對應(yīng)內(nèi)部環(huán)的自由能大小時(shí)，設(shè)置游標(biāo)i′、j′，使它們滿足條件i＜i′＜j′＜j，測試所有i′、j′可行解中是否存在第i′和第j′個(gè)堿基可以形成內(nèi)部堿基對，如果可以形成內(nèi)部堿基對則計(jì)算對應(yīng)子結(jié)構(gòu)的自由。

2 相關(guān)工作

根據(jù)上述動(dòng)態(tài)規(guī)劃算法，預(yù)測RNA二級結(jié)構(gòu)算法的時(shí)間復(fù)雜度為O（n4），但當(dāng)子序列很長時(shí)，O（n4）的時(shí)間復(fù)雜度變得難以接受。目前有兩種比較好的優(yōu)化方案：一種由HoHacker［3］提出的啟發(fā)式算法的解決辦法是設(shè)置閥值，將i′和j′的間隔控制在一定大小，如閥值設(shè)置為30最為合適，獲得較好加速的同時(shí)又不會(huì)產(chǎn)生很大偏差，此時(shí)時(shí)間復(fù)雜度為O（k2n2）（k表示閥值的大?。?；另外Lyngs等［4］指出在高溫情況下啟發(fā)式算法將閥值設(shè)置為30太小，從而使結(jié)果產(chǎn)生較大偏差，他們認(rèn)為長度相同的內(nèi)部環(huán)擁有相等的自由能，可以將長度相同的內(nèi)部環(huán)的自由能計(jì)算后存儲(chǔ)起來以供后續(xù)調(diào)用，新提出的這種方法可以將時(shí)間復(fù)雜度降低到 O（n3）。

另外開發(fā)并行實(shí)現(xiàn)對計(jì)算規(guī)模較大的程序進(jìn)行加速也是一個(gè)趨勢，GTfold是Mathuriya等提出的基于OpenMP的并行預(yù)測RNA二級結(jié)構(gòu)的程序，它具有很好的可擴(kuò)展性，相對于串行程序也有著不錯(cuò)的加速比。GTfold是一個(gè)使用共享存儲(chǔ)的并行程序，實(shí)現(xiàn)該程序的穩(wěn)定運(yùn)行對計(jì)算機(jī)系統(tǒng)的要求并不高，但是要想獲得較好的加速比，需要將GTDOLD運(yùn)行在擁有較大CPU緩存的系統(tǒng)上，因?yàn)閷?shí)驗(yàn)表明，在CPU緩存較小的系統(tǒng)上運(yùn)行時(shí)，其結(jié)果較于串行實(shí)現(xiàn)加速并不明顯。

3 算法的并行，數(shù)據(jù)和計(jì)算順序的描述

動(dòng)態(tài)規(guī)劃解決RNA二級結(jié)構(gòu)問題的核心偽碼如下所示：

input：RNA Sequence of Length N

begin

for bin 0，...，N do

＃calling device－kernel parallel for schedule

for i in 1，...，N－b do

j＝i＋b；

calcVBI（i，j）；

calcVM（i，j）；

calcV（i，j）；

calcWM（i，j）；

end

calcW（b＋1）；

end

return W（N）；

end

動(dòng)態(tài)規(guī)劃算法的計(jì)算模式會(huì)出現(xiàn)明顯的波前 （wave front）特性［5］，即調(diào)用函數(shù)進(jìn)行計(jì)算時(shí)，后層循環(huán)需要使用之前循環(huán)計(jì)算得到的結(jié)果作為輸入數(shù)據(jù)，如圖2所示，想要得到陰影三角中A點(diǎn)的值，需要AB、AC邊上所有的點(diǎn)作為已知量。在試圖對動(dòng)態(tài)規(guī)劃算法進(jìn)行并行時(shí)，上述依賴關(guān)系會(huì)為數(shù)據(jù)和任務(wù)的劃分帶來很大的困難，并使得計(jì)算負(fù)載嚴(yán)重不均衡。

圖2 RNA二級結(jié)構(gòu)預(yù)測算法數(shù)據(jù)依賴關(guān)系

具有波前特性的算法在執(zhí)行期間，如圖3所示，一次執(zhí)行過程中，沿波峰即斜對角線方向的計(jì)算元素之間沒有依賴關(guān)系。根據(jù)圖2，從橫向和縱向分析元素依賴性，顯然右側(cè)元素計(jì)算依賴于左側(cè)元素的數(shù)值，上側(cè)元素計(jì)算依賴于下側(cè)元素的數(shù)值。從而得到結(jié)論：由于每個(gè)橫行或者縱列之間的元素有相互依賴關(guān)系，單純橫向或者縱向推進(jìn)計(jì)算只能按照串行方式運(yùn)行程序；而每一道波峰內(nèi)部的元素沒有依賴關(guān)系，沿波峰方向推進(jìn)算法運(yùn)行，可以有效實(shí)現(xiàn)并行［5］。

在對具有波前特性的算法進(jìn)行并行時(shí)，可依照對角線元素的數(shù)量分配計(jì)算核心，但隨著程序的推進(jìn)對角線長度會(huì)逐漸變短，對應(yīng)的循環(huán)需要的并行核心數(shù)量也會(huì)逐漸減少，這樣對集群的資源會(huì)造成極大的浪費(fèi)。為了解決這種浪費(fèi)，Stratton等提出一種動(dòng)態(tài)劃分的方式，根據(jù)運(yùn)行中當(dāng)前整體計(jì)算量的大小，動(dòng)態(tài)決定分配到每個(gè)計(jì)算節(jié)點(diǎn)上的計(jì)算量，均衡考慮單機(jī)運(yùn)算量和集群節(jié)點(diǎn)間的通信開銷［9］。

圖3 一種并行RNA二級結(jié)構(gòu)預(yù)測算法的計(jì)算順序

考慮上述問題，負(fù)載均衡、數(shù)據(jù)傳輸?shù)仁怯绊懖⑿兴惴ㄐ实闹匾蛩?，但是對于單一的GPU而言，每個(gè)GPU有很多的計(jì)算核心，計(jì)算核心之間擁有通信速率很高的設(shè)備共享內(nèi)存，這樣我們可以解決數(shù)據(jù)傳輸問題。另外針對計(jì)算資源浪費(fèi)問題，CUDA編程模型可以充分利用GPU的計(jì)算核心，將整體負(fù)載分布到各個(gè)核心單元上，雖然當(dāng)規(guī)模降低到一定程度時(shí)依舊會(huì)產(chǎn)生資源浪費(fèi)，但由于此時(shí)整體計(jì)算量已不大，這并不會(huì)對程序執(zhí)行效率帶來太大的影響，同時(shí)負(fù)載均衡也有效利用計(jì)算核心的并行能力。由于CUDA在負(fù)載均衡上的優(yōu)點(diǎn)，Stratton等提出了一個(gè)翻譯框架，這個(gè)框架通過將一個(gè)線程塊轉(zhuǎn)換成一個(gè)單線程中的循環(huán)，進(jìn)而把CUDA程序編譯成可在多核CPU上運(yùn)行的程序，并指出這種方式生成的代碼性能及伸縮性很好［8］。從這一角度，我們認(rèn)為使用GPU預(yù)測RNA二級結(jié)構(gòu)是一個(gè)不錯(cuò)的選擇。

4 在GPU上實(shí)現(xiàn)RNA二級結(jié)構(gòu)預(yù)測算法

CUDA模型包含：線程組層次結(jié)構(gòu)、共享存儲(chǔ)器、屏蔽同步這3個(gè)重要的抽象概念［6］。CUDA擴(kuò)展C語言應(yīng)用上述3個(gè)抽象，可以很好指導(dǎo)并行程序開發(fā)人員把問題分割，找到能夠獨(dú)立解決的粗粒度子問題，將程序分割成能被并行協(xié)同解決的一個(gè)個(gè)小部分。這樣使開發(fā)人員，很方便的將傳統(tǒng)C編寫出的單線程層次結(jié)構(gòu)程序，轉(zhuǎn)為并行程序。

另外CUDA編寫的程序?qū)μ幚砥鞯暮诵臄?shù)量透明，一個(gè)CUDA程序可以在擁有任何數(shù)量流處理器的GPU設(shè)備上運(yùn)行。運(yùn)行時(shí)系統(tǒng)根據(jù)當(dāng)前硬件處理器的數(shù)量，由計(jì)算任務(wù)分配 （CWD）單元?jiǎng)討B(tài)進(jìn)行負(fù)載均衡。

約定在一臺(tái)計(jì)算機(jī)中，處理器和內(nèi)存等被稱為主機(jī)端（Host），搭載在計(jì)算機(jī)上的GPU被稱為設(shè)備端 （Device）。在基于CUDA模型的GPU上實(shí)現(xiàn)預(yù)測RAN二級結(jié)構(gòu)的并行算法時(shí)，算法中復(fù)雜指令優(yōu)化、獨(dú)立子問題的選擇以及運(yùn)行時(shí)的數(shù)據(jù)交換是影響程序執(zhí)行效率的幾個(gè)關(guān)鍵因素。

RNA二級結(jié)構(gòu)預(yù)測算法僅需要使用整型數(shù)據(jù)進(jìn)行計(jì)算，并不需要使用單精度或雙精度運(yùn)算，在這一點(diǎn)上無法利用GPU遠(yuǎn)超CPU的浮點(diǎn)性能。在獨(dú)立子問題選擇上，按照沿波峰方向?qū)栴}進(jìn)行。最后數(shù)據(jù)的交換方面主要分為兩個(gè)子問題，主機(jī)端和設(shè)備端之間的數(shù)據(jù)傳輸以及設(shè)備端計(jì)算數(shù)據(jù)的組織，而這也成為影響算法運(yùn)行效率的關(guān)鍵因素。

在一個(gè)支持CUDA模型的GPU中，除包含一定數(shù)量的流處理器外，還包含對應(yīng)的獨(dú)立設(shè)備端存儲(chǔ)單元。表1以GPU中的處理器為主視角，描述了各類存儲(chǔ)器的一些特性。其中GPU中的流處理器在主機(jī)端的內(nèi)存直接存儲(chǔ)數(shù)據(jù)延遲最大，其次是GPU內(nèi)部的全局存儲(chǔ)器，另外GPU內(nèi)部的有很多流多處理器，每個(gè)流多處理器配有一定數(shù)量自己獨(dú)立的寄存器、共享存儲(chǔ)器，這些存儲(chǔ)器的存取延時(shí)最小。當(dāng)主機(jī)端調(diào)用核函數(shù) （kernel）令設(shè)備端即GPU執(zhí)行程序時(shí)，每個(gè)線程塊 （block）都對應(yīng)分配到一個(gè)流多處理器 （SM）上，當(dāng)一個(gè)線程塊結(jié)束運(yùn)行，計(jì)算任務(wù)分配（CWD）單元為流多處理器分配一個(gè)新的線程塊。流多處理器中的共享存儲(chǔ)器可以視為大小固定，存取時(shí)間只需要一個(gè)時(shí)鐘的高速存儲(chǔ)器，且線程塊中的多個(gè)線程共享使用該存儲(chǔ)器，對于計(jì)算性能在1.2以上的CUDA來說，其共享存儲(chǔ)器大小至少為16KB。盡可能將數(shù)據(jù)放到更快的處理器上能夠有效提高處理器的計(jì)算效率。

表1 與GPU設(shè)備連接的各種存儲(chǔ)器特性比較

在實(shí)現(xiàn)中，GPU和CPU之間的數(shù)據(jù)傳輸考慮利用固定內(nèi)存 （pinned memory），即在主機(jī)端直接調(diào)用CUDA的cudaHostAlloc（）函數(shù)分配存儲(chǔ)空間，這樣分配出的內(nèi)存數(shù)據(jù)在程序運(yùn)行過程中，將一直保留在主機(jī)端內(nèi)存中，不會(huì)再因?yàn)閮?nèi)存不足被交換到對應(yīng)的虛擬內(nèi)存，權(quán)衡整機(jī)運(yùn)行效率和該類數(shù)據(jù)占用內(nèi)存空間大小，在本實(shí)驗(yàn)中能獲得一些性能提升。

GPU中流多處理器使用一種被稱為單指令、多線程（SIMT）的架構(gòu)，用來管理并行運(yùn)行的數(shù)百、上千個(gè)線程。在流多處理器SIMT單元中，并行線程以32個(gè)作為一組，形成warp，統(tǒng)一進(jìn)行創(chuàng)建、管理、調(diào)度和執(zhí)行。在分析程序運(yùn)行時(shí)數(shù)據(jù)的使用情況時(shí)，我們發(fā)現(xiàn)一些子函數(shù)中，如：計(jì)算內(nèi)部環(huán)能量時(shí)調(diào)用的函數(shù)eL（），經(jīng)過并行處理后，第j個(gè)線程在運(yùn)算時(shí)會(huì)需要使用數(shù)組RNA［N］中的RNA［j－2］、RNA［j－1］、RNA［j］、RNA［j＋1］、RNA［j＋2］等元素，也就是說第j個(gè)進(jìn)程和第j＋1個(gè)進(jìn)程之間至少會(huì)對RNA數(shù)組中的4個(gè)元素進(jìn)行重復(fù)的訪問。線程組中的多個(gè)線程在調(diào)用數(shù)據(jù)情況如圖4所示，并不會(huì)造成bank conflict，可以考慮把這些數(shù)據(jù)放到共享存儲(chǔ)器中。相對于從全局存儲(chǔ)器直接提取數(shù)據(jù)交給處理器，速度會(huì)有很大提升。另外，常量內(nèi)存也是所有流多處理器共享的存儲(chǔ)器，其中值由主機(jī)端直接操作而不允許在設(shè)備端更改，相對于全局存儲(chǔ)器其速度略高。在GPU上實(shí)現(xiàn)的RNA二級結(jié)構(gòu)預(yù)測程序中，包含一些大小固定、占用空間有限的索引數(shù)組，如包含256個(gè)元素用來計(jì)算stack loop能量的數(shù)組stack，保存發(fā)卡環(huán)不穩(wěn)定能量的數(shù)組hairpin、保存內(nèi)部環(huán)對應(yīng)的不穩(wěn)定能量的數(shù)組inter等。因?yàn)樯鲜鲇脕砑铀俪绦蛩惴ǖ臄?shù)組，在程序執(zhí)行期間，數(shù)組值并不會(huì)改變，放在常量內(nèi)存中 （對任意計(jì)算能力的CUDA，常量內(nèi)存大小為64KB）可以提高一些性能。

圖4 共享存儲(chǔ)器中數(shù)據(jù)讀取無沖突

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

RNA二級結(jié)構(gòu)預(yù)測算法可以分為兩個(gè)方面進(jìn)行評測：預(yù)測結(jié)果準(zhǔn)確度和算法的執(zhí)行時(shí)間。

當(dāng)前RNA二級結(jié)構(gòu)預(yù)測主要采用敏感性 （sensitivity）、陽性預(yù)測率 （positive predictive value，PPV）和馬休茲相互作用系數(shù)［7］（Matthews correlation coefficient）這3種方式來衡量預(yù)測準(zhǔn)確率。這里使用陽性預(yù)測率PPV來對實(shí)驗(yàn)產(chǎn)生的結(jié)果數(shù)據(jù)進(jìn)行評價(jià)。定義TP（true positive）：正確預(yù)測的堿基對個(gè)數(shù)；FP（false positive）沒有正確預(yù)測的堿基對個(gè)數(shù)，陽性預(yù)測率PPV＝TP／（TP＋FP）［10］。

算法性能分析，選擇同Mathuriya等提出的基于共享存儲(chǔ)的多線程并行實(shí)現(xiàn)GTFOLD進(jìn)行比較，其加速比相對與最好的串行算法約有19.8倍。

5.1 準(zhǔn)確度分析

準(zhǔn)確度分析選用Gutell database［7］中的數(shù)據(jù)作為基準(zhǔn)，表2是一組沉降系數(shù)為16S和23S（S是大分子物質(zhì)在超速離心沉降中的一個(gè)物理學(xué)單位，可反映分子質(zhì)量的大小［11］）的RNA序列，在這3種算法下得出的自由能的結(jié)果集。

表2 右側(cè)對應(yīng)描述了上述序列按照3種方法預(yù)測的自由能，以及以Gutell database為基準(zhǔn)得到的陽性預(yù)測率。

表2 GPUfold、UNAfold、RNAfold對一組RNA預(yù)測最小自由能和陽性預(yù)測率

使用PPV進(jìn)行準(zhǔn)確度測量，預(yù)測準(zhǔn)確率只和被正確預(yù)測的堿基對個(gè)數(shù)有關(guān)，從表格可以看到，3個(gè)程序?qū)NA二級結(jié)構(gòu)預(yù)測結(jié)果的準(zhǔn)確度相差不大。但是從整體看，使用計(jì)算機(jī)進(jìn)行RNA二級結(jié)構(gòu)預(yù)測的準(zhǔn)確度還是比較低，在這些程序中實(shí)現(xiàn)更多的分子熱力學(xué)細(xì)節(jié)能夠提高準(zhǔn)確度，同時(shí)實(shí)現(xiàn)這些細(xì)節(jié)對這3個(gè)程序的時(shí)間復(fù)雜度會(huì)帶來同等幅度增加，另外如果實(shí)現(xiàn)的細(xì)節(jié)包含對假結(jié)的處理，那么時(shí)間復(fù)雜度的增長幅度會(huì)很大。

5.2 性 能

實(shí)驗(yàn)選用的計(jì)算機(jī)包含一個(gè)主頻為3.3GHz的CPU，這個(gè)CPU有兩個(gè)大小為256KB的二級緩存以及一個(gè)大小為3MB的三級緩存，整機(jī)的主內(nèi)存大小為4GB，選用的并行GPU設(shè)備是GTX 550，它包含192個(gè)流處理器，設(shè)備內(nèi)部帶寬98.5GB／s，操作系統(tǒng)使用Fedora 15，內(nèi)核版本為linux 2.6，安裝的CUDA SDK （software development kit）版本是4.2。

圖5 中是使用在GPU上實(shí)現(xiàn)的程序GPUfold同GT－fold程序進(jìn)行對比測試，測試對象是一組長度為500－5000的RNA序列。從結(jié)果可以看到，隨著RNA序列長度增加，GPU上實(shí)現(xiàn)的程序GPUfold略優(yōu)于基于多核系統(tǒng)的實(shí)現(xiàn)GTfold。

圖5 GPUfold同GTfold效率對比

最后我們還對GPU核函數(shù)調(diào)度參數(shù)進(jìn)行變化，分析線程塊 （numBlocks）數(shù)量、線程塊中線程 （threadsPer－Block）的數(shù)量對執(zhí)行效率的影響。在這個(gè)試驗(yàn)中我們選用包含5000個(gè)堿基的RNA序列作為實(shí)驗(yàn)對象，圖6中的4條曲線選取的線程塊數(shù)量分別為8、4、2、1，我們發(fā)現(xiàn)隨著線程塊數(shù)量增加，程序能以更快速度得到結(jié)果，但是當(dāng)線程塊數(shù)量達(dá)到GPU設(shè)備中流多處理器數(shù)量兩倍之后，即使再增加線程塊數(shù)量，計(jì)算花費(fèi)時(shí)間將不會(huì)有明顯降低。當(dāng)核函數(shù)中線程塊數(shù)量固定，隨著線程數(shù)量增加，運(yùn)行速率也會(huì)有明顯提高，通過最終試驗(yàn)結(jié)果，我們發(fā)現(xiàn)取256作為線程數(shù)的數(shù)值能夠取得較好的效果。

圖6 線程塊數(shù)量及線程數(shù)對運(yùn)行時(shí)間的影響

6 結(jié)束語

本文提出了一種在基于CUDA模型GPU設(shè)備上并行計(jì)算預(yù)測RNA二級結(jié)構(gòu)的實(shí)現(xiàn)。在這個(gè)實(shí)現(xiàn)中，我們充分利用了GPU上共享存儲(chǔ)器、常量內(nèi)存等數(shù)據(jù)存儲(chǔ)設(shè)備，相較串行算法MFOLD等該實(shí)現(xiàn)運(yùn)行速度更快，同時(shí)與GT－fold并行算法在小型機(jī)上運(yùn)行有相同的準(zhǔn)確度。最終我們得到一個(gè)在GPU上實(shí)現(xiàn)的RNA二級結(jié)構(gòu)預(yù)測程序GPU－fold，該并行程序的運(yùn)行效率只和運(yùn)行設(shè)備中計(jì)算核心數(shù)量及其性能有關(guān)，執(zhí)行時(shí)間比較穩(wěn)定，沒有如同GTfold對緩存大小依賴。

［1］Markham NR，Zuker M.UNAFold：Software for nucleic acid folding and hybridization ［J］.Methods Mol Biol，2008：453：3－31.

［2］Mathews D H.Revolutions in RNA secondary structure prediction［J］.Journal of Molecular Biology，2006，359 （3）：526－532.

［3］Mathuriya A，Bader D A，Heitsch C E，et al.GTfold：a scalable multicore code for RNA secondary structure prediction［C］／／Hawaii，USA：ACM Symposium on Applied Computing，2009：981－988.

［4］Mathews DH，Turner DH，Zuker M.RNA secondary structure prediction ［M］.Current Protocols in Nucleic Acid Chemistry，2007：2944－2960.

［5］Du Z，Yin Z，Bader D.A tile－based parallel viterbi algorithm for biological sequence alignment on GPU with CUDA ［C］／／International Workshop on High Performance Computational Biology，2010：1－8.

［6］Nickolls J，Buck I，Skadron K，et al.Scalable parallel programming with CUDA ［J］.ACM Queue，2008，6 （2）：40－53.

［7］Cole J R，Wang Q.The ribosomal database project：Improved alignments and new tools for rRNA analysis ［J］.Nucleic Acids Res，2009，37 （suppl）：141－145.

［8］Stratton J A，Stone S S，Hwu W W.MCUDA：An efficient implementation of CUDA kernels on multicores［C］／／Canada：International Workshop on Languages and Compilers for Parallel Computing，2008：16－30.

［9］TAN Guangming，F(xiàn)ENG Shengzhong，SUN Ninghui.An optimized and efficiently parallelized dynamic programming for RNA secondary structure prediction ［J］Journal of Software，2006，17 （7）：1501－1509 （in Chinese）.［譚光明，馮圣中，孫寧輝.RNA二級結(jié)構(gòu)預(yù)測中動(dòng)態(tài)規(guī)劃的優(yōu)化和有效并行［J］.軟件學(xué)報(bào)，2006，17 （7）：1501－1509.］

［10］ZOU Quan，GUO Maozu，ZHANG Taotao.A review of RNA secondary structure prediction algorithms ［J］.Chinese Journal of Electronics，2008，36 （2）：331－336 （in Chinese）.［鄒權(quán)，郭茂祖，張濤濤.RNA二級結(jié)構(gòu)預(yù)測方法綜述 ［J］.電子學(xué)報(bào)，2008，36 （2）：331－336.］

［11］DONG Hao.The study on methods of RNA secondary structure prediction ［D］.Jilin：Jilin University，2011 （in Chinese）.［董浩.RNA二級結(jié)構(gòu)預(yù)測方法研究 ［D］.吉林：吉林大學(xué)，2011.］