基于禁忌遺傳算法的RNA二級結構預測

劉勇國，朱 嬋，晏 華

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

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

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

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

1 RNA二級結構

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

圖1 RNA二級結構示意圖

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

2 TGARNA算法

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

圖2 TGARNA算法流程

2.1 種群初始化

圖3 個體編碼方式

2.2 相容性檢測

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

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

2.3 適應度計算

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

2.4 選擇操作

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

2.5 交叉操作

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

2.6 變異操作

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

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

TGARNA算法步驟如下。

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

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

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

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

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

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

3 實驗結果

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

3.1 算法參數(shù)選擇

圖4 不同代溝值Gg的算法進化結果

表1 不同代溝值Gg的結果比較

圖5 不同交叉概率pc的算法進化結果

表2 不同交叉概率pc的結果比較

圖6 不同禁忌表長度因子α的算法進化結果

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

3.2 實驗結果討論

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

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

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

圖7 Giardia virus序列二級結構的自由能計算

表4 Giardia virus序列二級結構的預測結果

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

圖8 Y1 scRNA的二級結構的自由能計算

表5 Y1 scRNA序列二級結構的預測結果

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

圖9 U17 snoRNA序列二級結構的自由能計算

表6 U17 snoRNA二級結構的預測結果

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

圖10 PSTVd RNA的二級結構自由能計算

表7 PSTVd RNA二級結構的預測結果

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

4 結 論

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

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

[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] 鄒權, 郭茂祖, 張濤濤. RNA二級結構預測方法綜述[J].電子學報, 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] 任清華, 莫忠息, 陶玉敏. 預測RNA二級結構的一種遺傳模擬退火算法[J]. 武漢大學學報(理學版), 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二級結構的研究進展[J]. 生物信息學, 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] 胡桂武, 彭宏. 利用混沌差分進化算法預測RNA二級結構[J]. 計算機科學, 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二級結構預測算法[J]. 計算機工程與應用, 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.