基于OpenCL的RNA二級結構預測算法

汪方良，施慧彬

(南京航空航天大學 計算機科學與技術學院，江蘇 南京 211100)

基于OpenCL的RNA二級結構預測算法

汪方良，施慧彬

(南京航空航天大學 計算機科學與技術學院，江蘇 南京 211100)

包含假結的RNA二級結構預測在計算分子生物學中一直是一個重要的研究領域，而預測包含任意類型假結結構已被證明為NP完全問題。為了解決此類問題，在CPU平臺上實現了一種改進的遺傳算法。該算法可預測包含兩類假結結構的RNA序列，敏感性可達到0.775，陽性預測率可達到0.822 5。針對基于遺傳算法帶假結的RNA二級結構預測低效的問題，提出了基于OpenCL的異構并行加速算法。該算法在分析串行算法并行性的基礎上，在種群迭代進化階段進行異構加速，并基于GPU設備和OpenCL編程框架改進算法過程。為驗證所提算法的可行性和有效性，基于相同的測試集進行了實驗測試。測試結果表明，相對于串行算法，改進后的異構并行加速算法平均可實現2.72倍的速度提升，有效降低了RNA二級結構預測的耗時，提高了算法模擬預測效率。

RNA二級結構預測；假結；OpenCL；異構計算

0 引 言

RNA在基因表達中起到了十分重要的作用，對于每種RNA的功能分析，解析其結構特征是關鍵的一步。通過實驗方法分析RNA結構特征雖然精確，但成本較高，因此，通過計算方法預測RNA二級結構一直是近年來計算分子生物學領域比較熱門的課題之一。

目前預測RNA結構的算法大致可分為兩類：基于多序列比對的預測算法和單序列預測算法。多序列比對算法利用了同源RNA序列具有的相近遺傳信息與結構特性，預測精度較高，但是需要較多的先驗信息[1-2]。而單序列預測算法只需要輸入目標序列的一級結構，就可以預測出對應的二級結構且能保證一定的精確度。其中廣泛應用的算法思想是基于最小自由能(Minimum Free Energy,MFE)的二級結構預測。1999年，Zuker等提出用動態規劃算法(Dynamic Programming)計算各類環區自由能之和，從而預測一個單鏈RNA的二級結構[3]。在此基礎上，國內外研究人員進行了大量探索。何靜媛將人工魚群算法運用到RNA二級結構預測中[4]；邢翀運用模擬退火算法[5]提高了結構預測的效率；Sato K等在MFE模型的基礎上運用分割函數(Partition Function)進一步產生最大預期精度，進行二級結構預測[6]。但上述算法都沒有包含RNA結構中一類重要的子結構—假結。包含任意假結的RNA二級結構預測已被證實為NP-complete問題[7]，目前大部分算法會對假結的類型進行限制，啟發式算法在假結預測中效果較好。Ruan J等提出一種環匹配算法，不僅可以以單鏈作為輸入，同時支持輸入多條序列進行比對預測，從而提高精度[8]。類似的算法還有HotKnots[9]和FlexStem[10]，都是通過迭代添加子結構完成預測。Tong K K等提出了一種改進的遺傳算法[11]。基于改進遺傳算法中的三種遺傳操作，實現了包含兩類假結的RNA二級結構預測，并基于OpenCL實現了算法的并行加速，有效提高了結構預測的效率。

1 改進的遺傳算法

1.1相關定義

首先定義遺傳算法中的基因、染色體、個體、種群與RNA二級結構中相關概念的映射關系。

定義1：RNA二級結構中的每個莖區(Stem)代表一個基因；多個莖區的集合組成染色體(Chromosome)；每個染色體代表一個個體(Individual)；個體的集合構成種群(Population)。

對于莖區的概念與定義參見文獻[12]。改進后的遺傳算法包括三種遺傳操作。

定義2 交叉(Crossover)：在兩個個體i,j之間，從i中隨機選取一個莖區p，莖區p必須對個體j是有效莖區；從個體j中隨機選取一個莖區q，莖區q對于個體i是有效莖區，如能在i,j中找到這樣一對莖區p,q，則進行莖區交換。

假設有兩個個體{3,2,6,1,4}與{4,5,8,0,2}，個體{3,2,6,1,4}代表該個體由5個莖區組成，從左至右代表莖區在構成個體時被添加的先后順序。個體{3,2,6,1,4}與個體{4,5,8,0,2}進行交叉操作，交換了stem 1與stem 0得到個體{3,2,6,4,0}與個體{4,5,8,2,1}。一個莖區對于一個個體是有效的，稱為有效莖區，在不包含假結預測的算法中，莖區不存在交叉與重疊則是有效的，而提出的算法包含了假結結構預測，允許莖區交叉(具體交叉的類型在假結類型中詳細介紹)，允許莖區重疊兩個堿基對。根據Mathew&Tuner的標準能量模型中的定義，莖區的堿基對越多，長度越長，則自由能越小，那么在MFE的算法思想下，自然會選擇最長且沒有重疊的莖區構成個體，但是目前已經有部分研究證明，最優結果往往不完全出現在自由能最小的個體中，而是一些自由能次小的個體[13]。因此采用的算法中允許兩個堿基對的重疊存在，而當重疊發生后，后加入的莖區舍棄重疊區的堿基對，保留原有個體內莖區結構不變。

定義3 替換(Replacement)：從個體中隨機選取一個莖區p，從莖區池中隨機選取莖區q，如果莖區q對于該個體是有效莖區，則從該個體中剔除莖區p并加入莖區池，將莖區q從莖區池取出添加到個體中。

莖區池是在輸入RNA序列后根據堿基配對原則生成的所有候選莖區的集合，莖區池的構建將在下文算法流程中詳細介紹。需要注意的是，由于每個個體的進化都是獨立的，所以每個個體將獨立維護一個自有的莖區池。

定義4 添加(Addition)：從莖區池中隨機選取一個莖區p，若莖區p對于該個體是有效的，則從莖區池取出添加到個體中。

定義5 適應度(Fitness)：以個體的自由能值作為適應度，自由能越小，適應度越高。自由能的計算公式如下：

E=Estem+Ehairpin+Ebulge+Einternal+Emultibranch+Epknots

(1)

其中，Estem表示個體中所有莖區能量值；Ehairpin表示所有發卡環能量值；Ebulge表示所有凸環能量值；Einternal表示所有內環能量值；Emultibranch表示所有多分支環能量值。

前五個類型的子結構能量值計算都按標準能量模型MT計算，但是MT模型中不包括假結的情況，因此對于假結的能量計算采用D&P09能量模型中假結的計算辦法。對于前五個子結構的計算方法見文獻[14]。支持的假結類型定義如下：

定義6 假結(Pseudoknots)：一個RNA的二級結構就是一組堿基對的集合，對于所有的堿基對，如果存在兩對堿基對(i,j),(k,l)，其中i

RNA的假結結構有很多種，最常見的有H-type型假結[15]。提出的算法支持兩種假結結構的預測：H-type型假結和多分支環內含有一個或多個H-type假結，如圖1所示。

圖中箭頭左邊為RNA序列起始端即5’端，右邊帶箭頭的末端表示3’端[16-17]，兩端之間的圓點表示堿基，弧線表示兩端的堿基配對形成堿基對，相鄰的配對堿基構成Stem。如圖1(a)所示，H-type假結由兩個莖區交叉形成，圖1(b)中多分支環內包含一個H-type假結，提出的算法也支持包含多個H-type假結結構的預測。

圖1 兩類可預測假結結構

從公開RNA數據庫RNA STRAND中搜索到這兩種類型的真實RNA序列，見表1和表2。

表1 H-type假結序列

表2 多分支環包含假結序列

其中在Structure一行中，記錄的是該RNA的真實二級結構，“.”表示未配對堿基，括號和字幕表示配對的堿基。將上述兩個RNA序列作為輸入，檢測算法的有效性。

1.2算法步驟

(1)構建莖區池。

構建莖區池(StemList)的通常做法是根據堿基配對原則，找出所有可能的莖區，但實際上，很多長度較小的莖區并不會出現在真實的二級結構中。此外，窮舉所有可能的莖區構成莖區池會增大算法搜索莖區的范圍，增加算法的計算復雜度。文中算法采用最大莖區構建莖區池，所謂最大莖區就是指莖區池中只保留長度最大的莖區。比如有N對連續堿基對的莖區稱為長度為N的莖區，則該情況包含N*(N-1)/2個子莖區，若只保留最大莖區(長度為N的莖區)，則可以大大縮小莖區池的容量，在文中算法采用的數據集中，運用最大莖區構建的莖區池是通常辦法的1/3。

輸入長度為n的RNA一級序列，先構造一個n*n的矩陣，矩陣的i行j列即代表第i個堿基與第j個堿基。若這對堿基滿足Watson-Crick堿基配對原則(AU，GC)或者Wobble堿基配對原則(GU)，則矩陣中的(i,j)元素置1，否則置0。填充完矩陣所有元素后，以矩陣的下三角子矩陣為范圍，沿著次對角線的方向，逐一尋找連續的值為1的元素序列，每個元素序列即是最大莖區，所有的最大莖區構成莖區池。

(2)初始化種群。

構建完莖區池后，將莖區池中的莖區按照自由能從小到大排序，初始化初代個體算法如下：

算法：初始化種群

輸入：StemList

輸出：初代種群

Set all stems being unused

for each Individual in population

for each Stem in StemList

if Stem is used before then

continue

end if

if Stem is valid for Individual

then

Create temporary individualtempid

tempid=Individual.addstem(Stem)

if tempid.energy

then

Individual=tempid

Stem is used

end if

end if

end for

end for

為了保持種群個體的多樣性，在構造初始種群時，莖區池中的莖區每個只使用一次，這樣就保證了初始種群中每個個體的基因都不一樣。種群個體之間的差異性保證了種群進化的有效和健康。其次，在初始化初代種群時，每個個體都會維護一個自有的莖區池，比如自有莖區池當中有1～10，10個最大莖區，Individuali包含{2,3,5,7,9}五個莖區，那么其自有的莖區池則包含剩下的{1,4,6,8,10}五個莖區，自有莖區池在遺傳操作中使用的方式已在定義2～4中介紹了。

(3)種群迭代進化。

完成種群的初始化后，以初始化種群為輸入，進行迭代進化。對于種群中的所有個體，都會隨機發生Crossover、Replacement、Addition三種遺傳操作中的一個，即三種遺傳操作是等概率(1/3)發生。其次，在迭代進化過程中，常見的遺傳算法會設置迭代次數，文中算法設置了監視周期，監視周期之后，會監視種群內所有個體，如果個體的自由能不再發生變化，則結束迭代。監視周期設置為20輪迭代。迭代次數超過監視周期，并且種群內所有個體自由能不再降低，即是種群迭代進化停止的條件。算法步驟如下：

算法：種群迭代進化算法

輸入：初始化種群

輸出：進化完成的最優種群

While halting condition do not meet do

for each Individualiin population

Prob=rand()%100

Create 2 temporary individualt1,t2

if Prob<=crossover rate then

Randomly pick an individuala,wherea!=i

t1=Individuali,t2=Individuala

Crossover(t1,t2)

Ift2.energy

Individuala=t2

end if

else if Prob>crossover rate AND Prob<=replacement rate then

t1=Individuali

Replacement(t1)

else if Prob>replacement rate AND Prob<=addition rate then

t1=Individuali

Addition(t1)

end if

ift1.energy

Individuali=t1

end if

end for

end while

上述算法中，自由能即每個個體的適應度，計算方法如定義5，三種遺傳操作如定義2～4。當種群迭代進化完成后，輸出最終的種群個體，按照自由能由小到大進行排序，考慮到最優解有可能出現在能量次小的個體中，此處輸出自由能最小的5個個體，并考量算法預測的精確度。

1.3預測結果評估

目前較為流行的評估辦法主要從兩個維度進行測算：敏感性(Sensitivity)和陽性預測率(Positive Predictive Value)[18]，其公式如下：

(2)

(3)

其中，TP表示正確預測的堿基對數；FP表示錯誤預測的堿基對數；FN表示真實結構中應該存在但是預測結果中未能預測的堿基對數。

目前較為著名的開源軟件及算法，比如RNAfold、CentroidFold、ILM、HotKnots、pknotsRG等平均的敏感性在0.5～0.75之間，陽性預測率在0.65～0.75之間。

預測首先實現了一種改進的遺傳算法，預測包含兩種假結的RNA二級結構，但是在預測中，尤其是包含假結的預測，計算復雜度較高，耗時較長，而遺傳算法本身具有一定的可并行性，因此運用異構并行計算對上述算法進行改進和提升。

2 基于OpenCL的并行加速

2.1OpenCL簡介

OpenCL(Open Computing Language)即開放計算語言，是非盈利技術聯盟Khronos Group管理的異構編程框架。其主要作用是將異構設備用于并行算法加速，具有較好的跨平臺性。與技術特征類似的CUDA相比，OpenCL支持更多種類的加速設備，比如FPGA，DSP，CPU，GPU(包括部分移動端設備)。OpenCL能夠將不同架構的設備整合到統一框架下，釋放加速設備的計算性能，為通用計算提供更多的計算資源，提高計算效率。

OpenCL架構包含四種模型的定義：平臺模型、執行模型、內存模型、編程模型。平臺模型定義了異構平臺、加速設備、主設備、計算單元之間的關系；執行模型定義了內核(kernel)如何在設備上執行，以及如何配置執行內核所需的上下文等環境；內存模型定義了一個抽象的內存分層，可以保證開發人員無需關心底層的內存架構；編程模型定義了并發模型與物理硬件之間的映射關系。關于OpenCL的詳細內容參見文獻[19-21]。

2.2OpenCL執行步驟

(1)獲取可用的平臺信息；

(2)獲取可用的設備列表(GPU/CPU)；

(3)為加速設備創建上下文；

(4)為加速設備創建命令隊列，以便宿主機向加速設備發送各類命令；

(5)創建程序(program)對象，編譯內核代碼；

(6)創建內核對象，關聯內核代碼；

(7)設置內核參數；

(8)執行內核函數；

(9)讀取執行結果。

從上述流程可以發現，OpenCL程序最主要的部分在于kernel的設計與實現。針對需要加速的算法如何提取可并行部分，用kernel代替，在加速設備上實現并行運行是算法改進的關鍵。

2.3并行遺傳算法

對于一個算法的并行化改進，可以有兩種思路：任務并行、數據并行。任務并行是指，算法的兩個任務階段的輸入輸出不存在依賴關系，可以獨立執行，那么就可以安排兩個任務并行計算；數據并行是指，前后兩次計算的指令相同，只是輸入的數據不同。在上述介紹的改進遺傳算法中，可以發現種群中的每個個體的進化實際上是獨立自主的，并不一定需要等待其他個體進化完成，所以所有個體進化可以看作任務并行進行改進。其次，每個個體的遺傳操作(交叉、替換、添加)都是三種之一，只是輸入的數據不同，因此可以看作數據并行進行改進。算法流程如圖2所示。

圖2 并行遺傳算法

工作項(Work item)是OpenCL平臺模型中定義的最小計算單元，工作項最終會和加速設備上的計算核心相互映射。比如GPU平臺就有數千個計算核心，并行遺傳算法為大小為n的種群分配n個工作項，每個工作項獨立負責該個體的遺傳操作。這里的遺傳操作即是定義2～4中的過程。

但是三種遺傳操作中，交叉操作在并行情況下，有可能引起沖突：當個體a與個體b同時選中個體c(a、b、c互不相等)進行交叉操作時，即會引起沖突，導致畸形的后代。在提出算法中，畸形的后代體現為無效的莖區被保存在二級結構中，這樣的結果并不是所需要的，所以在每輪迭代后，并行算法會進行畸形后代的篩選，舍棄錯誤無效的輸出。通過實驗發現，這樣確實會帶來部分性能的損失，但是從單輪測試數據來看，每一百個個體中會產生五個左右的畸形后代，損失在可接受范圍以內。

對于健康的后代根據自由能的大小進行判斷，如果自由能小于父代，則用子代代替父代個體進入下一輪遺傳操作；如果子代自由能比父代大，則舍棄，仍由父代進入下一輪遺傳操作，直至迭代超出監視周期(20輪)，且所有個體的自由能都不再減小，則結束遺傳操作，進化完成。輸出自由能最小的五個個體。

3 實驗及結果分析

3.1有效性測試

為了測試算法的通用性，實驗數據不僅包括含假結的RNA序列，同時包括不含假結的RNA序列。實驗測試源數據來自RNA STRAND公開數據庫。測試長度在43～125 nt之間。預測結果對比見圖3。圖中，(a)、(c)、(e)是來自數據庫的RNA序列真實二級結構，(b)、(d)、(f)是所提算法預測結果通過RNA Movies軟件輸出的平面二級結構。其中(a)是RFA_00730號RNA序列，(c)、(e)是表1、表2當中定義的序列。

圖3 包含二級結構預測結果對比

從圖中可以看出，不含假結及單個H-type假結序列的預測準確率較高，而包含多個H-type假結組成多分支環的序列相較前兩種序列，準確度有所下降。從RNA的真實結構分析，多假結的結構復雜度更高，并且(e)的真實結構用目前的能量模型進行度量，并不是自由能最小的結構，也就是說多假結結構序列的最優解往往不在能量最小的情況，唯有進一步優化能量模型，才能提高多假結結構的預測準確率。綜合測試數據集，所提算法平均陽性預測率達到0.822 5，平均敏感性達到0.775。說明該算法有效，預測結果有一定的參考價值。

3.2并行算法加速測試

(1)硬件環境：AMD A10-7400P 2.5-3.4 GHz，集成顯卡AMD Radeon R6 Graphic，獨立顯卡 AMD Radeon R9 M280X，可用內存6.94 GB。

(2)軟件環境：Windows 8.1(64bit)，驅動版本Catalyst 15.7，Visual Studio Ultimate 2012，AMD APP SDK 3.0。

以表2中的PDB_00447號序列為例，序列長度為120 nt，種群大小設置為300，監視周期為20，加速測試結果如圖4所示。

圖4 加速測試結果

如圖4所示，該序列在140輪迭代后完成種群進化，以5輪迭代為間隔，對比了CPU(串行算法)和GPU(并行算法)的耗時。雖然在并行算法中，需要對因交叉沖突產生的畸形后代進行過濾，損失部分性能，但是從圖中可以看出，并行算法相對串行算法仍有較大優勢。從數據來看，平均加速比達到2.72x。

4 結束語

為了解決包含假結在內的RNA二級結構預測問題，提出并實現了一種改進的遺傳算法。該算法能夠預測包含兩種假結結構在內的RNA二級結構。為了提升算法的預測效率，在保證預測結果有效的情況下，基于OpenCL對算法進行并行化改進與加速，平均獲得2.72倍的速度提升，有效提高了算法的計算效率。但對于多假結結構的預測仍有提升空間，下一步工作是優化能量模型，提高多假結結構的預測精度。另外，并行算法每輪迭代之間有繁瑣的內存傳遞過程，目前部分加速設備已經支持OpenCL 2.0中最新的共享虛擬內存(Shared Virtual Memory)，運用SVM技術優化訪存過程，提升加速效果也是后續研究工作之一。

[1] 張濤濤.基于比較序列分析的RNA二級結構預測算法研究[D].哈爾濱：哈爾濱工業大學,2007.

[2] 方小永.基于比較序列分析的RNA二級結構預測與評估[D].長沙：國防科學技術大學,2007.

[3] Mathews D H,Sabina J,Zuker M,et al.Expanded sequence dependence of thermodynamic parameters improves prediction of RNA secondary structure[J].Journal of Molecular Biology，1999,288(5):911-940.

[4] 何靜媛.RNA二級結構預測算法的研究[D].重慶：重慶大學,2009.

[5] 邢 翀.RNA二級結構預測算法的研究[D].長春：吉林大學,2012.

[6] Sato K,Kato Y,Hamada M,et al.IPknot:fast and accurate prediction of RNA secondary structures with pseudoknots using integer programming[J].Bioinformatics，2011,27(13):85-93.

[7] Lyngso R B,Pedersen C N.RNA pseudoknot prediction in energy-based models[J].Journal of Computational Biology，2000,7(3-4):409-427.

[8] Ruan J,Stormo G D,Zhang W.An iterated loop matching approach to the prediction of RNA secondary structures with pseudoknots[J].Bioinformatics,2004,20(1):58-66.

[9] Ren J,Rastegari B,Condon A,et al.HotKnots:heuristic prediction of RNA secondary structures including pseudoknots[J].RNA，2005,11(10):1494-1504.

[10] Chen X,He S M,Bu D,et al.FlexStem:improving predictions of RNA secondary structures with pseudoknots by reducing the search space[J].Bioinformatics，2008,24(18):1994-2001.

[11] Tong K K,Cheung K Y,Lee K H,et al.GAknot:RNA secondary structures prediction with pseudoknots using genetic algorithm[C]//IEEE symposium on computational intelligence in bioinformatics and computational biology.[s.l.]:IEEE,2013:136-142.

[12] 彭 政.帶假結的RNA二級結構預測算法研究[D].長沙：湖南大學,2008.

[13] Staple D W,Butcher S E.Pseudoknots:RNA structures with diverse functions[J].Plos Biology,2005,3(6):213.

[14] Andronescu M S,Pop C,Condon A E.Improved free energy parameters for RNA pseudoknotted secondary structure prediction[J].RNA,2009,16(1):26-42.

[15] 胥 杰.基于混沌模擬退火的RNA二級結構預測的研究[D].成都：電子科技大學,2010.

[16] 劉元寧，張 浩，李 誌，等.RNA假結結構分析[J].吉林大學學報:工學版,2009,39(S1):265-269.

[17] 高世樂,丁克詮.含假結RNA二級結構類的圖語法[J].計算機工程與應用,2008,44(2):23-25.

[18] 劉振棟.包含假結的RNA結構預測算法研究[D].濟南：山東大學,2014.

[19] Gaster B.OpenCL異構計算[M].張云泉,張先軼,龍國平,等,譯.第2版.北京：清華大學出版社,2012.

[20] 詹 云,趙新燦,譚同德.基于OpenCL的異構系統并行編程[J].計算機工程與設計,2012,33(11):4191-4195.

[21] 陳 鋼，吳百鋒.面向OpenCL模型的GPU性能優化[J].計算機輔助設計與圖形學學報,2011,23(4):571-581.

Secondary Structure Prediction of RNA Based on OpenCL

WANG Fang-liang，SHI Hui-bin

(College of Computer Science and Technology,Nanjing University of Aeronautics and Astronautics,Nanjing 211100,China)

Predicting RNA secondary structure is an important field in computational molecular biology especially including pseudoknots.However,predicting RNA secondary structure with all kinds of pseudoknots has been proven to be an NP-complete problem.To solve it,an improved genetic algorithm is proposed in CPU platform,which can predict two kinds of pseudoknots.Its sensitivity can reach 0.775 and its positive predictive value can reach 0.822 5.The prediction of RNA secondary structure with pseudoknots based on genetic algorithm is inefficient.To solve it,an accelerated algorithm based on OpenCL is presented,which accelerates the period of individual evolution according to the analysis of parallelizability of serial prediction algorithm.Then the algorithm established with GPU based on OpenCL is promoted.The contrast experiments with the same test set have been conducted compared with other algorithms.The experimental results show that the improved heterogeneous parallel algorithm has acquired 2.72 times faster average operation rate than others,reducing the computing time effectively and improving the efficiency of prediction.

RNA secondary structure;pseudoknots;OpenCL;heterogeneous computing

2016-10-31

：2017-02-17 < class="emphasis_bold">網絡出版時間

時間：2017-07-11

國家“973”重點基礎研究發展計劃項目(2014CB744900)

汪方良(1991-)，男，碩士研究生，研究方向為異構計算、計算機體系結構；施慧彬，博士，副教授，研究方向為計算機體系結構、可重構計算。

http://kns.cnki.net/kcms/detail/61.1450.TP.20170711.1456.078.html

TP311

:A

:1673-629X(2017)09-0001-06

10.3969/j.issn.1673-629X.2017.09.001