劉勇國，朱 嬋，晏 華

(1. 電子科技大學(xué)計算機(jī)科學(xué)與工程學(xué)院 成都 611731; 2. 蘇州大學(xué)江蘇省計算機(jī)信息處理技術(shù)重點實驗室 江蘇 蘇州 215006;3. 四川建筑職業(yè)技術(shù)學(xué)院圖書館 四川 德陽 618000)

生物體根據(jù)脫氧核糖核酸(deoxy ribonucleic acid，DNA)的遺傳信息合成蛋白質(zhì)并完成各項生命活動過程，如生物信號的識別和傳遞，營養(yǎng)物質(zhì)的運輸?shù)?。在該過程中，蛋白質(zhì)合成的直接模板是核糖核酸(ribonucleic acid，RNA)而非DNA。隨著對RNA分子結(jié)構(gòu)探索的不斷深入，研究發(fā)現(xiàn)某些簡單生物體如煙草花葉病毒的RNA已成為遺傳信息的直接載體，負(fù)責(zé)存儲遺傳信息和蛋白質(zhì)合成。目前，研究RNA分子結(jié)構(gòu)成為認(rèn)知生物遺傳信息和蛋白質(zhì)合成的重要途徑。RNA分子由一級、二級和三級結(jié)構(gòu)[1]組成。一級結(jié)構(gòu)為多核苷酸鏈的核苷酸序列；二級結(jié)構(gòu)為多核苷酸鏈回折形成的配對雙螺旋區(qū)和不對稱區(qū)域；三級結(jié)構(gòu)為二級結(jié)構(gòu)單元相互作用形成的三維空間結(jié)構(gòu)，確定二級結(jié)構(gòu)單元的空間定位取向和蛋白質(zhì)合成。RNA一級結(jié)構(gòu)預(yù)測起步最早，方法較完善，主要通過生物化學(xué)方法直接預(yù)測。二級結(jié)構(gòu)預(yù)測主要采用物理化學(xué)和生物信息學(xué)方法。三級結(jié)構(gòu)預(yù)測已被證明為NP-hard問題[2]。由于二級結(jié)構(gòu)介于一級結(jié)構(gòu)和三級結(jié)構(gòu)之間，存儲了較多高級結(jié)構(gòu)信息，因此研究RNA二級結(jié)構(gòu)成為預(yù)測RNA分子結(jié)構(gòu)的重要途徑。

RNA二級結(jié)構(gòu)預(yù)測包括基于物理化學(xué)的方法和基于生物信息學(xué)的方法[3]兩類。物理化學(xué)的方法通過X射線結(jié)晶衍射和核磁共振確定RNA分子結(jié)構(gòu)，雖然測量精度高，但對軟硬件要求較高，而且RNA分子降解較快，難以結(jié)晶，造成預(yù)測困難，所以，該類方法適用于堿基數(shù)量小于100的RNA二級結(jié)構(gòu)預(yù)測問題。生物信息學(xué)方法包含系統(tǒng)發(fā)育比較技術(shù)和自由能最小技術(shù)[2]兩類預(yù)測技術(shù)。系統(tǒng)發(fā)育比較技術(shù)對比待測序列與已知同源序列，通過同源序列二級結(jié)構(gòu)的形成方式預(yù)測待測序列的二級結(jié)構(gòu)[1]。自由能最小技術(shù)源于熱動力學(xué)的能量耗散原理，通過設(shè)置RNA二級結(jié)構(gòu)中莖環(huán)結(jié)構(gòu)的熱動力學(xué)參數(shù)，模擬實際RNA分子熱運動，建立二級結(jié)構(gòu)的能量模型，計算RNA序列折疊后的自由能，確定具備最小自由能的二級結(jié)構(gòu)為預(yù)測結(jié)果，實現(xiàn)對RNA二級結(jié)構(gòu)的預(yù)測[2]，預(yù)測方法包括動態(tài)規(guī)劃算法[4]、退火模擬算法[5]、遺傳算法[6]、上下文無關(guān)文法預(yù)測算法[7]、協(xié)變信息預(yù)測算法[3]等。由于遺傳算法具有隱含并行性、自適應(yīng)性等特點，已應(yīng)用于RNA二級結(jié)構(gòu)預(yù)測問題。研究發(fā)現(xiàn)，遺傳算法的種群進(jìn)化過程缺乏記憶性，搜索過程易出現(xiàn)種群早熟現(xiàn)象，導(dǎo)致預(yù)測精度降低[4,8-9]。文獻(xiàn)[10]設(shè)計并行遺傳算法模擬RNA序列形成二級結(jié)構(gòu)時的折疊情況，推測RNA序列的二級結(jié)構(gòu)。文獻(xiàn)[11]采用十進(jìn)制編碼設(shè)計遺傳算法預(yù)測RNA二級結(jié)構(gòu)，以改善預(yù)測效率和精度[8,11-12]。文獻(xiàn)[5]組合模擬退火與遺傳算法，以增強(qiáng)預(yù)測算法的局部搜索能力，避免種群早熟現(xiàn)象。文獻(xiàn)[13]基于混沌差分進(jìn)化算法預(yù)測RNA二級結(jié)構(gòu)，利用混沌映射的偽隨機(jī)性和遍歷性提高算法的全局搜索能力。文獻(xiàn)[14]提出基于免疫粒子群集成的RNA二級結(jié)構(gòu)預(yù)測方法，通過并行群落優(yōu)化方式實現(xiàn)協(xié)同演化。

本文提出基于禁忌遺傳算法的RNA二級結(jié)構(gòu)預(yù)測算法TGARNA。TGARNA算法給出莖區(qū)相容性檢測過程，保留最長莖區(qū)構(gòu)造莖區(qū)相容個體；將禁忌思想引入遺傳操作，記錄個體家族特征和進(jìn)化過程，防止近親繁殖，保持種群多樣性。仿真實驗表明，TGARNA算法能夠獲取最小自由能并預(yù)測RNA二級結(jié)構(gòu)。

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

除少數(shù)病毒RNA為雙鏈結(jié)構(gòu)外，RNA分子通常為多核苷酸單鏈結(jié)構(gòu)，包含腺嘌呤(adenine，A)、鳥嘌呤(guanine，G)、胞嘧啶(cytosine，C)和尿嘧啶(uralic，U)4類含氮堿基。堿基與戊糖結(jié)合形成核苷，戊糖側(cè)鏈與磷酸分子結(jié)合形成核苷酸，核苷酸分子以3’和5’磷酸二酯鍵連接形成核苷酸序列[1]。RNA單鏈自身回折呈現(xiàn)堿基配對和單鏈交替出現(xiàn)的莖環(huán)結(jié)構(gòu)，形成RNA二級結(jié)構(gòu)，如圖1所示。

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

RNA二級結(jié)構(gòu)中，堿基以A-U、G-C或G-U互補(bǔ)配對形成的連續(xù)雙螺旋區(qū)域稱為莖區(qū)，不構(gòu)成互補(bǔ)配對的單鏈結(jié)構(gòu)稱為環(huán)區(qū)，環(huán)區(qū)類型包括發(fā)夾環(huán)、凸環(huán)、內(nèi)環(huán)和內(nèi)分支環(huán)等。莖區(qū)由氫鍵相連的配對堿基疊加形成，氫鍵數(shù)量越多堿基對自由能越小，越有利于增強(qiáng)二級結(jié)構(gòu)穩(wěn)定性。環(huán)區(qū)數(shù)量增加導(dǎo)致自由能增長，造成二級結(jié)構(gòu)穩(wěn)定性下降[2]。RNA二級結(jié)構(gòu)數(shù)目隨序列長度增加呈指數(shù)級增長，因此有效預(yù)測RNA二級結(jié)構(gòu)成為RNA分子結(jié)構(gòu)研究急待解決的問題[2]。堿基是RNA二級結(jié)構(gòu)的基本單元，連續(xù)配對堿基集形成莖區(qū)，未配對堿基構(gòu)成環(huán)區(qū)，相容莖區(qū)集唯一確定RNA二級結(jié)構(gòu)[2]。RNA二級結(jié)構(gòu)的相關(guān)描述如下[5]。

2 TGARNA算法

TGARNA算法過程基于遺傳算法，通過相容性檢測構(gòu)造合法種群，引入禁忌思想防止近親繁殖，提高種群多樣性[15]，算法流程如圖2所示。

圖2 TGARNA算法流程

2.1 種群初始化

圖3 個體編碼方式

2.2 相容性檢測

由于不相容莖區(qū)組合將生成非法RNA二級結(jié)構(gòu)，本文設(shè)計相容性檢測環(huán)節(jié)建立合法個體。常用相容性檢測方法根據(jù)莖區(qū)排序位置確定莖區(qū)是否保留，未考慮其能量值，易造成低能量莖區(qū)漏選[5]。鑒于RNA二級結(jié)構(gòu)中莖區(qū)數(shù)目越多，莖區(qū)長度越長，自由能越小，故應(yīng)盡量保留長相容莖區(qū)。本文相容性檢測步驟如下：

通過相容性檢測避免非法個體，保留低能量莖區(qū)，改善種群性能。

2.3 適應(yīng)度計算

本文采用文獻(xiàn)[16]提出的最小自由能模型評價RNA二級結(jié)構(gòu)穩(wěn)定度，預(yù)測RNA二級結(jié)構(gòu)。RNA二級結(jié)構(gòu)自由能計算為：

2.4 選擇操作

本文采用比例選擇方式，個體適應(yīng)度越大，其選擇概率越大。個體Xi的選擇概率為：

2.5 交叉操作

交叉操作采用融合禁忌思想的單點交叉方式，針對個體Xi的莖區(qū)組合Xsi進(jìn)行。TGARNA算法將家族特征和禁忌思想引入交叉操作，防止近親個體繁殖后代，保持種群多樣性；同時，通過期望準(zhǔn)則允許禁忌個體產(chǎn)生高適應(yīng)度后代，增加選擇壓力，提高算法收斂速度。交叉操作步驟如下。

2.6 變異操作

本文采用動態(tài)變異方式產(chǎn)生后代個體，針對莖區(qū)組合更新變異后個體的家族特征和禁忌表。變異操作步驟如下。

2.7 TGARNA算法的實現(xiàn)

TGARNA算法步驟如下。

1) 給定RNA序列，建立莖區(qū)池并生成初始種群，設(shè)g=1；

2) 執(zhí)行種群相容性檢測；

3) 執(zhí)行選擇操作；

4) 執(zhí)行交叉操作；

5) 執(zhí)行變異操作；

6) 將父代種群中Nd個最優(yōu)個體更新后代種群中適應(yīng)度最低的Nd個個體；

3 實驗結(jié)果

仿真實驗平臺采用Windows XP SP3，Matlab7.1語言，Intel Core 2 Duo 2.20 GHz CPU，2G物理內(nèi)存，仿真算法每次實驗運行20次。首先討論算法參數(shù)的選擇過程，然后分析TGARNA算法的預(yù)測結(jié)果。

3.1 算法參數(shù)選擇

圖4 不同代溝值Gg的算法進(jìn)化結(jié)果

表1 不同代溝值Gg的結(jié)果比較

圖5 不同交叉概率pc的算法進(jìn)化結(jié)果

表2 不同交叉概率pc的結(jié)果比較

圖6 不同禁忌表長度因子α的算法進(jìn)化結(jié)果

表3 不同禁忌表長度因子α的結(jié)果比較

3.2 實驗結(jié)果討論

比較基于遺傳算法的RNA二級結(jié)構(gòu)預(yù)測方法(genetic algorithm based prediction of rna secondary structure，GARNA)、融入相容性檢測的基于遺傳算法的RNA二級結(jié)構(gòu)預(yù)測方法GARNA*和TGARNA算法，GARNA算法采用比例選擇、單點交叉和單點變異實現(xiàn)遺傳操作；GARNA*算法在GARNA算法的基礎(chǔ)上改進(jìn)莖區(qū)相容性檢測；TGARNA算法在GARNA*的基礎(chǔ)上融入禁忌思想。算法實驗種群規(guī)模均為100?？紤]到RNA序列長度差異，針對Giardia virus和Y1 scRNA序列，進(jìn)化代數(shù)設(shè)為100；針對U17 snoRNA和PSTVd RNA序列，進(jìn)化代數(shù)設(shè)為1 000。通過調(diào)整進(jìn)化代數(shù)適應(yīng)不同長度序列RNA分子的二級結(jié)構(gòu)預(yù)測。

仿真實驗采用Giardia virus、Y1 scRNA、U17 snoRNA和PSTVd RNA共4個RNA測試序列。Giardia virus序列源于賈第蟲病毒體，由77個堿基構(gòu)成，包含4個莖區(qū)、2個發(fā)夾環(huán)和1個多分支環(huán)。Y1 scRNA序列源于小家鼠的線形染色體組RNA，由111個堿基構(gòu)成，包含2個莖區(qū)、1個發(fā)夾環(huán)和1個內(nèi)環(huán)[17]。U17 snoRNA序列源于沼澤側(cè)頸龜?shù)霓D(zhuǎn)錄RNA，由240個堿基構(gòu)成，包含9個莖區(qū)、3個發(fā)夾環(huán)、3個凸環(huán)、5個內(nèi)環(huán)和1個多分支環(huán)。PSTVd RNA序列源于馬鈴薯紡錘塊莖類病毒體的線形染色體組RNA，由359個堿基構(gòu)成，包含25個莖區(qū)、1個發(fā)卡環(huán)和23個內(nèi)環(huán)[17]。

Giardia virus RNA二級結(jié)構(gòu)自由能的計算過程如圖7所示，該序列二級結(jié)構(gòu)的預(yù)測結(jié)果如表4所示。預(yù)測指標(biāo)包括平均最小自由能、最小自由能均方差、算法首次獲得最小自由能的平均運行時間和正確預(yù)測莖區(qū)數(shù)與實際二級結(jié)構(gòu)中的莖區(qū)數(shù)比?？梢姡珿ARNA算法的自由能下降緩慢，GARNA*與TGARNA算法的自由能下降迅速。GARNA*算法的運行時間下降到 0.6 s時，預(yù)測精度提高到56.52%，TGARNA算法的預(yù)測精度進(jìn)一步提高到60.12%，正確預(yù)測3個莖區(qū)、1個發(fā)夾環(huán)和1個多分支環(huán)。

圖7 Giardia virus序列二級結(jié)構(gòu)的自由能計算

表4 Giardia virus序列二級結(jié)構(gòu)的預(yù)測結(jié)果

Y1 scRNA二級結(jié)構(gòu)自由能計算如圖8所示。該序列二級結(jié)構(gòu)的預(yù)測結(jié)果如表5所示。GARNA算法收斂速度較慢，GARNA*和TGARNA算法能較快獲得最小自由能。GARNA*算法的運行時間下降到2.1 s時，預(yù)測精度達(dá)到100%，正確預(yù)測Y1 scRNA序列的二級結(jié)構(gòu)。TGARNA算法的運行時間進(jìn)一步下降到0.9 s。

圖8 Y1 scRNA的二級結(jié)構(gòu)的自由能計算

表5 Y1 scRNA序列二級結(jié)構(gòu)的預(yù)測結(jié)果

U17 snoRNA二級結(jié)構(gòu)自由能的計算過程如圖9所示。該序列二級結(jié)構(gòu)的預(yù)測結(jié)果如表6所示。GARNA*算法性能改善明顯，獲取最小自由能的速度較TGARNA算法快，GARNA算法預(yù)測效率最低。TGARNA算法較GARNA*算法的預(yù)測精度高，正確預(yù)測5個莖區(qū)、2個發(fā)夾環(huán)和1個凸環(huán)。

圖9 U17 snoRNA序列二級結(jié)構(gòu)的自由能計算

表6 U17 snoRNA二級結(jié)構(gòu)的預(yù)測結(jié)果

PSTVd RNA二級結(jié)構(gòu)自由能計算如圖10所示。該序列二級結(jié)構(gòu)的預(yù)測結(jié)果如表7所示。GARNA算法性能最差，TGARNA算法較GARNA*算法的最小自由能下降。TGARNA算法預(yù)測效率最高并達(dá)到最高預(yù)測精度，正確預(yù)測21個莖區(qū)、1個發(fā)夾環(huán)和16個內(nèi)環(huán)。

圖10 PSTVd RNA的二級結(jié)構(gòu)自由能計算

表7 PSTVd RNA二級結(jié)構(gòu)的預(yù)測結(jié)果

實驗表明，改進(jìn)莖區(qū)相容性檢測能夠改善種群性能，引入禁忌思想能夠防止近親個體繁殖，預(yù)測算法的收斂速度和預(yù)測效率得到改善。

4 結(jié) 論

RNA分子結(jié)構(gòu)是生物信息學(xué)領(lǐng)域的重要研究對象，通過RNA結(jié)構(gòu)研究病毒的遺傳信息，探索RNA參與合成蛋白質(zhì)及轉(zhuǎn)錄遺傳信息的過程。本文給出基于禁忌遺傳算法的RNA二級結(jié)構(gòu)預(yù)測方法TGARNA，通過設(shè)計莖區(qū)相容性檢測，保留長莖區(qū)以降低個體自由能，將家族特征和禁忌思想融入個體編碼，防止近親繁殖，保持種群多樣性。仿真實驗表明，除U17 snoRNA序列外，TGARNA算法能夠?qū)崿F(xiàn)較高的預(yù)測效率和精度。

由于RNA二級結(jié)構(gòu)預(yù)測的最小自由能模型忽略了某些結(jié)構(gòu)(如假結(jié))，造成預(yù)測算法對U17 snoRNA序列的預(yù)測精度較低。假結(jié)結(jié)構(gòu)形成復(fù)雜，文獻(xiàn)[1]證明使用最小自由能模型預(yù)測包含任意假結(jié)的RNA二級結(jié)構(gòu)為NP完全問題[1]。后續(xù)研究工作中，將考慮基于TGARNA算法將個體設(shè)計為由單位結(jié)構(gòu)組成，利用TGARNA算法計算相容單元結(jié)構(gòu)集合形成初始RNA二級結(jié)構(gòu)，通過環(huán)區(qū)配對確定被測序列中的假結(jié)結(jié)構(gòu)，預(yù)測含假結(jié)的RNA二級結(jié)構(gòu)。

[1] PAR S K, BANDYOPADHYAY S, RAY S S. Evolutionary computation in bioinformatics: a review[J]. IEEE Transactions on Systems, Man, and Cybernetics, Part C:Applications and Reviews, 2006, 36(5): 601-615.

[2] COHEN J. Bioinformatics - an introduction for computer scientists[J]. ACM Computer Surveys, 2004, 36(2):122-158.

[3] 鄒權(quán), 郭茂祖, 張濤濤. RNA二級結(jié)構(gòu)預(yù)測方法綜述[J].電子學(xué)報, 2008, 36(2): 331-337.ZOU Quan, GUO Mao-zu, ZHANG Tao-tao. A review of RNA secondary structure prediction algorithms[J]. Acta Electronica Sinica, 2008, 36(2): 331-337.

[4] RIVAS E, EDDY S R. A dynamic programming algorithm for RNA structure prediction including pseudoknots[J].Journal of Molecular Biology, 1999, 285(5): 2053-2068.

[5] 任清華, 莫忠息, 陶玉敏. 預(yù)測RNA二級結(jié)構(gòu)的一種遺傳模擬退火算法[J]. 武漢大學(xué)學(xué)報(理學(xué)版), 2004, 50(1):23-28.REN Qing-hua, MO Zhong-xi, TAO Yu-min. A geneticsimulated-annealing algorithm for predicting RNA secondary structure[J]. Journal of Wuhan University(Natural Science Edition), 2004, 50(1): 23-28.

[6] VAN BATENBURG F H D, GULTYAEV A P, PLEIJ C W A.An APL programmed genetic algorithm for the prediction of RNA secondary structure[J]. Journal of Theoretical Biology,1995, 174(3): 269-280.

[7] 唐四薪, 劉艷波, 尹軍. 文法推斷RNA二級結(jié)構(gòu)的研究進(jìn)展[J]. 生物信息學(xué), 2008, 6(4): 190-192.TANG Si-xin, LIU Yan-bo, YIN Jun. Research advances of grammatical inference of RNA secondary structure[J]. China Journal of Bioinformatics, 2008, 6(4): 190-192.

[8] WIESE K C, DESCHENES A A, HENDRIKS A G.RnaPredict—an evolutionary algorithm for RNA secondary structure prediction[J]. IEEE/ACM Transactions on Computational Biology and Bioinformatics, 2008, 5(1):25-41.

[9] LIU Y G, CHEN K F, LIAO X F, ZHANG W. A genetic clustering method for intrusion detection[J]. Pattern Recognition, 2004, 37(5): 927-942.

[10] SHAPIRO B A, WU J C, BENGALI D, POTTS M J. The massively parallel genetic algorithm for RNA folding MIMD implementation and population variation[J].Bioinformatics, 2001, 17(2): 137-148.

[11] WIESE K C, GLEN E. A permutation-based genetic algorithm for the RNA folding problem: a critical look at selection strategies, crossover operators and representation issues[J]. Biosystems, 2003, 72(1-2): 29-41.

[12] WIESE K C, GLEN E. A permutation based genetic algorithm for RNA secondary structure prediction[C]//In:Proceedings of the 2nd International Conference on Hybrid Intelligent Systems. Chile: [s.n.], 2002, 173-182.

[13] 胡桂武, 彭宏. 利用混沌差分進(jìn)化算法預(yù)測RNA二級結(jié)構(gòu)[J]. 計算機(jī)科學(xué), 2007, 34(9): 163-166.HU Gui-wu, PENG Hong. An algorithm-base chaos differential evolution for predicting RNA secondary structure[J]. Computer Science, 2007, 34(9): 163-166.

[14] 胡桂武, 彭宏. 基于免疫粒子群集成的RNA二級結(jié)構(gòu)預(yù)測算法[J]. 計算機(jī)工程與應(yīng)用, 2007, 43(3): 26-29.HU Gui-wu, PENG Hong. Algorithm based on immune PSO ensemble for predicting RNA secondary structure[J].Computer Engineering and Applications, 2007, 43(3):26-29.

[15] TING K C, LI H T, LEE H N. On the harmonious mating strategy throug tabu search[J]. Information Sciences, 2003,156(3-4): 189-214.

[16] MATHEWS D H, SABINA J, ZUKER M, TURNER D H.Expanded sequence dependence of thermodynamic parameters improves prediction of RNA secondary structure[J]. Journal of Molecular Biology, 1999, 288(5):911-940.