汪方良，施慧彬

(南京航空航天大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，江蘇 南京 211100)

基于OpenCL的RNA二級(jí)結(jié)構(gòu)預(yù)測(cè)算法

汪方良，施慧彬

(南京航空航天大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，江蘇 南京 211100)

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

RNA二級(jí)結(jié)構(gòu)預(yù)測(cè)；假結(jié)；OpenCL；異構(gòu)計(jì)算

0 引 言

RNA在基因表達(dá)中起到了十分重要的作用，對(duì)于每種RNA的功能分析，解析其結(jié)構(gòu)特征是關(guān)鍵的一步。通過(guò)實(shí)驗(yàn)方法分析RNA結(jié)構(gòu)特征雖然精確，但成本較高，因此，通過(guò)計(jì)算方法預(yù)測(cè)RNA二級(jí)結(jié)構(gòu)一直是近年來(lái)計(jì)算分子生物學(xué)領(lǐng)域比較熱門的課題之一。

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

1 改進(jìn)的遺傳算法

1.1相關(guān)定義

首先定義遺傳算法中的基因、染色體、個(gè)體、種群與RNA二級(jí)結(jié)構(gòu)中相關(guān)概念的映射關(guān)系。

定義1：RNA二級(jí)結(jié)構(gòu)中的每個(gè)莖區(qū)(Stem)代表一個(gè)基因；多個(gè)莖區(qū)的集合組成染色體(Chromosome)；每個(gè)染色體代表一個(gè)個(gè)體(Individual)；個(gè)體的集合構(gòu)成種群(Population)。

對(duì)于莖區(qū)的概念與定義參見(jiàn)文獻(xiàn)[12]。改進(jìn)后的遺傳算法包括三種遺傳操作。

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

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

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

莖區(qū)池是在輸入RNA序列后根據(jù)堿基配對(duì)原則生成的所有候選莖區(qū)的集合，莖區(qū)池的構(gòu)建將在下文算法流程中詳細(xì)介紹。需要注意的是，由于每個(gè)個(gè)體的進(jìn)化都是獨(dú)立的，所以每個(gè)個(gè)體將獨(dú)立維護(hù)一個(gè)自有的莖區(qū)池。

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

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

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

(1)

其中，Estem表示個(gè)體中所有莖區(qū)能量值；Ehairpin表示所有發(fā)卡環(huán)能量值；Ebulge表示所有凸環(huán)能量值；Einternal表示所有內(nèi)環(huán)能量值；Emultibranch表示所有多分支環(huán)能量值。

前五個(gè)類型的子結(jié)構(gòu)能量值計(jì)算都按標(biāo)準(zhǔn)能量模型MT計(jì)算，但是MT模型中不包括假結(jié)的情況，因此對(duì)于假結(jié)的能量計(jì)算采用D&P09能量模型中假結(jié)的計(jì)算辦法。對(duì)于前五個(gè)子結(jié)構(gòu)的計(jì)算方法見(jiàn)文獻(xiàn)[14]。支持的假結(jié)類型定義如下：

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

RNA的假結(jié)結(jié)構(gòu)有很多種，最常見(jiàn)的有H-type型假結(jié)[15]。提出的算法支持兩種假結(jié)結(jié)構(gòu)的預(yù)測(cè)：H-type型假結(jié)和多分支環(huán)內(nèi)含有一個(gè)或多個(gè)H-type假結(jié)，如圖1所示。

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

圖1 兩類可預(yù)測(cè)假結(jié)結(jié)構(gòu)

從公開(kāi)RNA數(shù)據(jù)庫(kù)RNA STRAND中搜索到這兩種類型的真實(shí)RNA序列，見(jiàn)表1和表2。

表1 H-type假結(jié)序列

表2 多分支環(huán)包含假結(jié)序列

其中在Structure一行中，記錄的是該RNA的真實(shí)二級(jí)結(jié)構(gòu)，“.”表示未配對(duì)堿基，括號(hào)和字幕表示配對(duì)的堿基。將上述兩個(gè)RNA序列作為輸入，檢測(cè)算法的有效性。

1.2算法步驟

(1)構(gòu)建莖區(qū)池。

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

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

(2)初始化種群。

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

算法：初始化種群

輸入：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

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

(3)種群迭代進(jìn)化。

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

算法：種群迭代進(jìn)化算法

輸入：初始化種群

輸出：進(jìn)化完成的最優(yōu)種群

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

上述算法中，自由能即每個(gè)個(gè)體的適應(yīng)度，計(jì)算方法如定義5，三種遺傳操作如定義2～4。當(dāng)種群迭代進(jìn)化完成后，輸出最終的種群個(gè)體，按照自由能由小到大進(jìn)行排序，考慮到最優(yōu)解有可能出現(xiàn)在能量次小的個(gè)體中，此處輸出自由能最小的5個(gè)個(gè)體，并考量算法預(yù)測(cè)的精確度。

1.3預(yù)測(cè)結(jié)果評(píng)估

目前較為流行的評(píng)估辦法主要從兩個(gè)維度進(jìn)行測(cè)算：敏感性(Sensitivity)和陽(yáng)性預(yù)測(cè)率(Positive Predictive Value)[18]，其公式如下：

(2)

(3)

其中，TP表示正確預(yù)測(cè)的堿基對(duì)數(shù)；FP表示錯(cuò)誤預(yù)測(cè)的堿基對(duì)數(shù)；FN表示真實(shí)結(jié)構(gòu)中應(yīng)該存在但是預(yù)測(cè)結(jié)果中未能預(yù)測(cè)的堿基對(duì)數(shù)。

目前較為著名的開(kāi)源軟件及算法，比如RNAfold、CentroidFold、ILM、HotKnots、pknotsRG等平均的敏感性在0.5～0.75之間，陽(yáng)性預(yù)測(cè)率在0.65～0.75之間。

預(yù)測(cè)首先實(shí)現(xiàn)了一種改進(jìn)的遺傳算法，預(yù)測(cè)包含兩種假結(jié)的RNA二級(jí)結(jié)構(gòu)，但是在預(yù)測(cè)中，尤其是包含假結(jié)的預(yù)測(cè)，計(jì)算復(fù)雜度較高，耗時(shí)較長(zhǎng)，而遺傳算法本身具有一定的可并行性，因此運(yùn)用異構(gòu)并行計(jì)算對(duì)上述算法進(jìn)行改進(jìn)和提升。

2 基于OpenCL的并行加速

2.1OpenCL簡(jiǎn)介

OpenCL(Open Computing Language)即開(kāi)放計(jì)算語(yǔ)言，是非盈利技術(shù)聯(lián)盟Khronos Group管理的異構(gòu)編程框架。其主要作用是將異構(gòu)設(shè)備用于并行算法加速，具有較好的跨平臺(tái)性。與技術(shù)特征類似的CUDA相比，OpenCL支持更多種類的加速設(shè)備，比如FPGA，DSP，CPU，GPU(包括部分移動(dòng)端設(shè)備)。OpenCL能夠?qū)⒉煌軜?gòu)的設(shè)備整合到統(tǒng)一框架下，釋放加速設(shè)備的計(jì)算性能，為通用計(jì)算提供更多的計(jì)算資源，提高計(jì)算效率。

OpenCL架構(gòu)包含四種模型的定義：平臺(tái)模型、執(zhí)行模型、內(nèi)存模型、編程模型。平臺(tái)模型定義了異構(gòu)平臺(tái)、加速設(shè)備、主設(shè)備、計(jì)算單元之間的關(guān)系；執(zhí)行模型定義了內(nèi)核(kernel)如何在設(shè)備上執(zhí)行，以及如何配置執(zhí)行內(nèi)核所需的上下文等環(huán)境；內(nèi)存模型定義了一個(gè)抽象的內(nèi)存分層，可以保證開(kāi)發(fā)人員無(wú)需關(guān)心底層的內(nèi)存架構(gòu)；編程模型定義了并發(fā)模型與物理硬件之間的映射關(guān)系。關(guān)于OpenCL的詳細(xì)內(nèi)容參見(jiàn)文獻(xiàn)[19-21]。

2.2OpenCL執(zhí)行步驟

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

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

(3)為加速設(shè)備創(chuàng)建上下文；

(4)為加速設(shè)備創(chuàng)建命令隊(duì)列，以便宿主機(jī)向加速設(shè)備發(fā)送各類命令；

(5)創(chuàng)建程序(program)對(duì)象，編譯內(nèi)核代碼；

(6)創(chuàng)建內(nèi)核對(duì)象，關(guān)聯(lián)內(nèi)核代碼；

(7)設(shè)置內(nèi)核參數(shù)；

(8)執(zhí)行內(nèi)核函數(shù)；

(9)讀取執(zhí)行結(jié)果。

從上述流程可以發(fā)現(xiàn)，OpenCL程序最主要的部分在于kernel的設(shè)計(jì)與實(shí)現(xiàn)。針對(duì)需要加速的算法如何提取可并行部分，用kernel代替，在加速設(shè)備上實(shí)現(xiàn)并行運(yùn)行是算法改進(jìn)的關(guān)鍵。

2.3并行遺傳算法

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

圖2 并行遺傳算法

工作項(xiàng)(Work item)是OpenCL平臺(tái)模型中定義的最小計(jì)算單元，工作項(xiàng)最終會(huì)和加速設(shè)備上的計(jì)算核心相互映射。比如GPU平臺(tái)就有數(shù)千個(gè)計(jì)算核心，并行遺傳算法為大小為n的種群分配n個(gè)工作項(xiàng)，每個(gè)工作項(xiàng)獨(dú)立負(fù)責(zé)該個(gè)體的遺傳操作。這里的遺傳操作即是定義2～4中的過(guò)程。

但是三種遺傳操作中，交叉操作在并行情況下，有可能引起沖突：當(dāng)個(gè)體a與個(gè)體b同時(shí)選中個(gè)體c(a、b、c互不相等)進(jìn)行交叉操作時(shí)，即會(huì)引起沖突，導(dǎo)致畸形的后代。在提出算法中，畸形的后代體現(xiàn)為無(wú)效的莖區(qū)被保存在二級(jí)結(jié)構(gòu)中，這樣的結(jié)果并不是所需要的，所以在每輪迭代后，并行算法會(huì)進(jìn)行畸形后代的篩選，舍棄錯(cuò)誤無(wú)效的輸出。通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)，這樣確實(shí)會(huì)帶來(lái)部分性能的損失，但是從單輪測(cè)試數(shù)據(jù)來(lái)看，每一百個(gè)個(gè)體中會(huì)產(chǎn)生五個(gè)左右的畸形后代，損失在可接受范圍以內(nèi)。

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

3 實(shí)驗(yàn)及結(jié)果分析

3.1有效性測(cè)試

為了測(cè)試算法的通用性，實(shí)驗(yàn)數(shù)據(jù)不僅包括含假結(jié)的RNA序列，同時(shí)包括不含假結(jié)的RNA序列。實(shí)驗(yàn)測(cè)試源數(shù)據(jù)來(lái)自RNA STRAND公開(kāi)數(shù)據(jù)庫(kù)。測(cè)試長(zhǎng)度在43～125 nt之間。預(yù)測(cè)結(jié)果對(duì)比見(jiàn)圖3。圖中，(a)、(c)、(e)是來(lái)自數(shù)據(jù)庫(kù)的RNA序列真實(shí)二級(jí)結(jié)構(gòu)，(b)、(d)、(f)是所提算法預(yù)測(cè)結(jié)果通過(guò)RNA Movies軟件輸出的平面二級(jí)結(jié)構(gòu)。其中(a)是RFA_00730號(hào)RNA序列，(c)、(e)是表1、表2當(dāng)中定義的序列。

圖3 包含二級(jí)結(jié)構(gòu)預(yù)測(cè)結(jié)果對(duì)比

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

3.2并行算法加速測(cè)試

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

(2)軟件環(huán)境：Windows 8.1(64bit)，驅(qū)動(dòng)版本Catalyst 15.7，Visual Studio Ultimate 2012，AMD APP SDK 3.0。

以表2中的PDB_00447號(hào)序列為例，序列長(zhǎng)度為120 nt，種群大小設(shè)置為300，監(jiān)視周期為20，加速測(cè)試結(jié)果如圖4所示。

圖4 加速測(cè)試結(jié)果

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

4 結(jié)束語(yǔ)

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

[1] 張濤濤.基于比較序列分析的RNA二級(jí)結(jié)構(gòu)預(yù)測(cè)算法研究[D].哈爾濱：哈爾濱工業(yè)大學(xué),2007.

[2] 方小永.基于比較序列分析的RNA二級(jí)結(jié)構(gòu)預(yù)測(cè)與評(píng)估[D].長(zhǎng)沙：國(guó)防科學(xué)技術(shù)大學(xué),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二級(jí)結(jié)構(gòu)預(yù)測(cè)算法的研究[D].重慶：重慶大學(xué),2009.

[5] 邢 翀.RNA二級(jí)結(jié)構(gòu)預(yù)測(cè)算法的研究[D].長(zhǎng)春：吉林大學(xué),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] 彭 政.帶假結(jié)的RNA二級(jí)結(jié)構(gòu)預(yù)測(cè)算法研究[D].長(zhǎng)沙：湖南大學(xué),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二級(jí)結(jié)構(gòu)預(yù)測(cè)的研究[D].成都：電子科技大學(xué),2010.

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

[17] 高世樂(lè),丁克詮.含假結(jié)RNA二級(jí)結(jié)構(gòu)類的圖語(yǔ)法[J].計(jì)算機(jī)工程與應(yīng)用,2008,44(2):23-25.

[18] 劉振棟.包含假結(jié)的RNA結(jié)構(gòu)預(yù)測(cè)算法研究[D].濟(jì)南：山東大學(xué),2014.

[19] Gaster B.OpenCL異構(gòu)計(jì)算[M].張?jiān)迫?張先軼,龍國(guó)平,等,譯.第2版.北京：清華大學(xué)出版社,2012.

[20] 詹 云,趙新?tīng)N,譚同德.基于OpenCL的異構(gòu)系統(tǒng)并行編程[J].計(jì)算機(jī)工程與設(shè)計(jì),2012,33(11):4191-4195.

[21] 陳 鋼，吳百鋒.面向OpenCL模型的GPU性能優(yōu)化[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),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">網(wǎng)絡(luò)出版時(shí)間

時(shí)間：2017-07-11

國(guó)家“973”重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃項(xiàng)目(2014CB744900)

汪方良(1991-)，男，碩士研究生，研究方向?yàn)楫悩?gòu)計(jì)算、計(jì)算機(jī)體系結(jié)構(gòu)；施慧彬，博士，副教授，研究方向?yàn)橛?jì)算機(jī)體系結(jié)構(gòu)、可重構(gòu)計(jì)算。

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