基于CUDA對RNA二級結構預測的并行研究

陳 飛，郝福珍

（華北計算技術研究所，北京100083）

0 引 言

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

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

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

1 RNA二級結構預測算法及動態規劃求解

RNA是由核苷酸組合成的聚合物，其包含A、C、G、U這4種堿基，這4種堿基由A－U、C－G和U－G這3種組合配對，并由氫鍵連接而成的結構被稱為堿基對，配對全集 ｛（A－U），（U－A），（C－G），（G－C），（U－G），（G－U）｝。預測RNA二級結構是指，取給定RNA堿基序列中的各個位點進行配對，將得到的配對結果轉換成二級模型的過程。其中使用計算機用動態規劃算法得到的結果不包括偽結。

1.1 RNA二級結構預測

圖1 是由79個堿基組成的RNA的二級結構，在這個結構中包含發卡環 （hairpin loop）、內部環 （internal loop）、多分支環（multiloop）、莖（stack）、 外部環（external loop）等5種子結構。

圖1 RNA二級結構示例

其中莖是由連續的堿基對堆積構成的結構，它在三維空間中的表現為A型雙螺旋，這種子結構對整體的RNA二級結構的穩定性起促進作用，即令二級結構的自由能降低；除此之外的幾種子結構，比如：莖中間出現一些無法配對的堿基形成的內部環、相鄰連續序列因兩端互補形成發卡結構的發卡環，會對整體RNA結構的穩定性造成破壞，即令二級結構的自由能升高。

在最小自由能 （minimum free energy）模型中，一個RNA二級結構的自由能，可以通過將組成它的所有獨立子結構，如環 （internal loop）、莖 （stack）等所對應的自由能累加得到，其中每個子結構的能量只依賴于環本身的結構，而和整體結構中其它部分無關［2］。每個RNA二級結構都有自身對應的自由能，擁有最小自由能的二級結構最穩定，被稱為一個RNA序列的最優二級結構。

1.2 動態規劃求解算法

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

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

對于上述長度為N的RNA堿基序列s，其子序列s1…sj的最優結構的自由能W（j）為

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

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

eH（i，j）是序列si…sj形成發卡環 （hairpin loop）結構所包含的自由能，該發卡環序列中第i和第j兩個核苷酸組成堿基對，并作為子結構整體對外的連接點。

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

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

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

VM（i，j）是序列si…sj構成擁有最優結構的多分支環包含的自由能。其值大小同單鏈結構的數量有線性關系，引入WM來計算VM，其關系式如下

用動態規劃算法求解RNA二級結構問題，其空間復雜度為O（n2），其中包含V、VBI、VM、WM等4個占用空間大小為（N＋1）＊（N＋1）／2的數組，數組大小與RNA序列長度有關，動態規劃前期保存中間結果以求出的最小自由能 （minimum free energy，MFE），回溯階段作為得到最優結構的關鍵信息。在RNA二級結構預測中回溯需要的時間遠小于求最小自由能的時間。

使用動態規劃算法計算最小自由能的時間復雜度為O（n4），預測一個給定序列的時間復雜度為：O（n2）＊子序列計算時間復雜度。

子序列的最小自由能的計算方法由上述式 （1）～式（5）給出，其中計算內部環和多分支環的時間復雜度最大。當計算子序列si…sj對應內部環的自由能大小時，設置游標i′、j′，使它們滿足條件i＜i′＜j′＜j，測試所有i′、j′可行解中是否存在第i′和第j′個堿基可以形成內部堿基對，如果可以形成內部堿基對則計算對應子結構的自由。

2 相關工作

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

另外開發并行實現對計算規模較大的程序進行加速也是一個趨勢，GTfold是Mathuriya等提出的基于OpenMP的并行預測RNA二級結構的程序，它具有很好的可擴展性，相對于串行程序也有著不錯的加速比。GTfold是一個使用共享存儲的并行程序，實現該程序的穩定運行對計算機系統的要求并不高，但是要想獲得較好的加速比，需要將GTDOLD運行在擁有較大CPU緩存的系統上，因為實驗表明，在CPU緩存較小的系統上運行時，其結果較于串行實現加速并不明顯。

3 算法的并行，數據和計算順序的描述

動態規劃解決RNA二級結構問題的核心偽碼如下所示：

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

動態規劃算法的計算模式會出現明顯的波前 （wave front）特性［5］，即調用函數進行計算時，后層循環需要使用之前循環計算得到的結果作為輸入數據，如圖2所示，想要得到陰影三角中A點的值，需要AB、AC邊上所有的點作為已知量。在試圖對動態規劃算法進行并行時，上述依賴關系會為數據和任務的劃分帶來很大的困難，并使得計算負載嚴重不均衡。

圖2 RNA二級結構預測算法數據依賴關系

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

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

圖3 一種并行RNA二級結構預測算法的計算順序

考慮上述問題，負載均衡、數據傳輸等是影響并行算法效率的重要因素，但是對于單一的GPU而言，每個GPU有很多的計算核心，計算核心之間擁有通信速率很高的設備共享內存，這樣我們可以解決數據傳輸問題。另外針對計算資源浪費問題，CUDA編程模型可以充分利用GPU的計算核心，將整體負載分布到各個核心單元上，雖然當規模降低到一定程度時依舊會產生資源浪費，但由于此時整體計算量已不大，這并不會對程序執行效率帶來太大的影響，同時負載均衡也有效利用計算核心的并行能力。由于CUDA在負載均衡上的優點，Stratton等提出了一個翻譯框架，這個框架通過將一個線程塊轉換成一個單線程中的循環，進而把CUDA程序編譯成可在多核CPU上運行的程序，并指出這種方式生成的代碼性能及伸縮性很好［8］。從這一角度，我們認為使用GPU預測RNA二級結構是一個不錯的選擇。

4 在GPU上實現RNA二級結構預測算法

CUDA模型包含：線程組層次結構、共享存儲器、屏蔽同步這3個重要的抽象概念［6］。CUDA擴展C語言應用上述3個抽象，可以很好指導并行程序開發人員把問題分割，找到能夠獨立解決的粗粒度子問題，將程序分割成能被并行協同解決的一個個小部分。這樣使開發人員，很方便的將傳統C編寫出的單線程層次結構程序，轉為并行程序。

另外CUDA編寫的程序對處理器的核心數量透明，一個CUDA程序可以在擁有任何數量流處理器的GPU設備上運行。運行時系統根據當前硬件處理器的數量，由計算任務分配 （CWD）單元動態進行負載均衡。

約定在一臺計算機中，處理器和內存等被稱為主機端（Host），搭載在計算機上的GPU被稱為設備端 （Device）。在基于CUDA模型的GPU上實現預測RAN二級結構的并行算法時，算法中復雜指令優化、獨立子問題的選擇以及運行時的數據交換是影響程序執行效率的幾個關鍵因素。

RNA二級結構預測算法僅需要使用整型數據進行計算，并不需要使用單精度或雙精度運算，在這一點上無法利用GPU遠超CPU的浮點性能。在獨立子問題選擇上，按照沿波峰方向對問題進行。最后數據的交換方面主要分為兩個子問題，主機端和設備端之間的數據傳輸以及設備端計算數據的組織，而這也成為影響算法運行效率的關鍵因素。

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

表1 與GPU設備連接的各種存儲器特性比較

在實現中，GPU和CPU之間的數據傳輸考慮利用固定內存 （pinned memory），即在主機端直接調用CUDA的cudaHostAlloc（）函數分配存儲空間，這樣分配出的內存數據在程序運行過程中，將一直保留在主機端內存中，不會再因為內存不足被交換到對應的虛擬內存，權衡整機運行效率和該類數據占用內存空間大小，在本實驗中能獲得一些性能提升。

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

圖4 共享存儲器中數據讀取無沖突

5 實驗結果

RNA二級結構預測算法可以分為兩個方面進行評測：預測結果準確度和算法的執行時間。

當前RNA二級結構預測主要采用敏感性 （sensitivity）、陽性預測率 （positive predictive value，PPV）和馬休茲相互作用系數［7］（Matthews correlation coefficient）這3種方式來衡量預測準確率。這里使用陽性預測率PPV來對實驗產生的結果數據進行評價。定義TP（true positive）：正確預測的堿基對個數；FP（false positive）沒有正確預測的堿基對個數，陽性預測率PPV＝TP／（TP＋FP）［10］。

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

5.1 準確度分析

準確度分析選用Gutell database［7］中的數據作為基準，表2是一組沉降系數為16S和23S（S是大分子物質在超速離心沉降中的一個物理學單位，可反映分子質量的大?。?1］）的RNA序列，在這3種算法下得出的自由能的結果集。

表2 右側對應描述了上述序列按照3種方法預測的自由能，以及以Gutell database為基準得到的陽性預測率。

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

使用PPV進行準確度測量，預測準確率只和被正確預測的堿基對個數有關，從表格可以看到，3個程序對RNA二級結構預測結果的準確度相差不大。但是從整體看，使用計算機進行RNA二級結構預測的準確度還是比較低，在這些程序中實現更多的分子熱力學細節能夠提高準確度，同時實現這些細節對這3個程序的時間復雜度會帶來同等幅度增加，另外如果實現的細節包含對假結的處理，那么時間復雜度的增長幅度會很大。

5.2 性 能

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

圖5 中是使用在GPU上實現的程序GPUfold同GT－fold程序進行對比測試，測試對象是一組長度為500－5000的RNA序列。從結果可以看到，隨著RNA序列長度增加，GPU上實現的程序GPUfold略優于基于多核系統的實現GTfold。

圖5 GPUfold同GTfold效率對比

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

圖6 線程塊數量及線程數對運行時間的影響

6 結束語

本文提出了一種在基于CUDA模型GPU設備上并行計算預測RNA二級結構的實現。在這個實現中，我們充分利用了GPU上共享存儲器、常量內存等數據存儲設備，相較串行算法MFOLD等該實現運行速度更快，同時與GT－fold并行算法在小型機上運行有相同的準確度。最終我們得到一個在GPU上實現的RNA二級結構預測程序GPU－fold，該并行程序的運行效率只和運行設備中計算核心數量及其性能有關，執行時間比較穩定，沒有如同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，FENG 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二級結構預測中動態規劃的優化和有效并行［J］.軟件學報，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）.［鄒權，郭茂祖，張濤濤.RNA二級結構預測方法綜述 ［J］.電子學報，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二級結構預測方法研究 ［D］.吉林：吉林大學，2011.］