多目標遺傳算法的含假結RNA二級結構預測

顧 倜，蔡磊鑫，王 帥，呂 強

(1.蘇州大學 計算機科學與技術學院，蘇州 215006;2.江蘇省計算機信息處理技術重點實驗室(蘇州大學)，蘇州 215006)

多目標遺傳算法的含假結RNA二級結構預測

顧 倜1，蔡磊鑫1，王 帥1，呂 強2*

(1.蘇州大學 計算機科學與技術學院，蘇州 215006;2.江蘇省計算機信息處理技術重點實驗室(蘇州大學)，蘇州 215006)

假結是RNA中一種重要的結構，由于建模的困難導致它更難被預測。通過堿基之間的配對概率來預測含假結RNA二級結構的ProbKnot算法具有很高的精度，但該算法僅用了配對概率作為預測依據，導致陰性配對大量出現，因此精度中的特異性較低。實驗結合ProbKnot算法中堿基配對概率模型，通過使用多目標遺傳算法，從而提高預測含假結RNA二級結構的特異性，以此促進總體精度的提高。實驗過程中，首先計算出每個堿基成為單鏈的概率，作為新增的預測依據，然后使用遺傳算法對RNA二級結構進行交叉、變異和迭代，最后得到Pareto最優解，進一步得出最高的最大期望精度。實驗結果表明，在使用的RNA案例中，采用該方法比現有方法精度平均提高約4%。

RNA二級結構; 假結; 多目標優化; 遺傳算法; 最大期望精度

在最初的生物中心法則中，RNA被認為在表達遺傳信息時作為蛋白質翻譯，發揮了短暫的作用。后來發現RNA除了翻譯蛋白質，還具有多種其他功能，如調控基因表達，轉運RNA，催化肽鏈形成和指導蛋白質合成等。隨著研究的深入，RNA的形象已經由功能單一的簡單線性堿基序列演變成種類多樣、結構復雜、功能特異的生命核心物，同時它在中心法則中取得了與DNA和蛋白質同樣重要的地位。許多RNA序列具有準確定義的結構，想要理解這些RNA序列如何實現它們的功能，知道它們的結構是非常重要的[1]。

RNA中最基本的4種堿基為A(腺嘌呤)、U(尿嘧啶)、G(鳥嘌呤)、C(胞嘧啶)。RNA一級結構是堿基的有序序列。RNA二級結構是由堿基之間相互作用形成的3種規范堿基對AU，GC和GU以及未配對堿基的單鏈組成的結構[2]。RNA三級結構是原子的三維排列。因為RNA的結構一般是分層次的，所以二級結構可以在不用知道三級結構的情況下得到確定，同樣也可以為預測三級結構提供很大的幫助。RNA二級結構承擔著由RNA一級結構推測其三級結構的橋梁作用。研究RNA二級結構，可以為提高RNA及蛋白質結構預測的準確率提供幫助。因此，RNA結構研究的一個重要切入點是對其二級結構的研究。

鑒于二級結構的重要性，對二級結構的預測研究顯得尤為重要。預測RNA二級結構的方法主要有序列比對法和最小自由能法。其中, 序列比對法需要耗費大量人力，預測效率較低，因而最小自由能法是預測RNA二級結構時廣泛采用的一種方法。Zuker[3]提出的Mfold算法是早期基于最小自由能方法的二級結構預測算法，該算法有很高的正確率，但不能預測含假結的復雜結構。二級結構預測中，假結作為一種特殊的結構，是實現結構功能的重要因素。但是因為堿基相互重疊的特性，生物信息計算對假結的預測難以奏效。預測包含假結的RNA二級結構問題是NP難的[4]，但要分析RNA的真實結構，對假結的準確預測是必須解決的問題。

通過加入假結約束來預測假結的最小自由能算法是目前使用比較多的預測方法，其中Probknot算法[5-6]在幾種RNA含假結二級結構預測算法[7-10]中精準度較高。Probknot算法的配分函數(partition function)[11]利用機器學習法和熱力學模型[12]計算出每個堿基之間的配對概率，各個堿基通過配對概率相互配對，然后解除形成的結構中連續配對堿基數較少的莖區，以確保穩定性[13]，最終得出預測結果。本文運用多目標遺傳算法，通過提高特異性從而提高總體精度的方式來改進Probknot算法，并與Probknot算法以及其他算法的結果進行對比和總結。

1 打分函數與預測性能指標

1.1 打分函數MaxExpect

相比最小自由方法的最小自由能結構，MaxExpect利用RNA結構中各個堿基配對情況打分得到的最大期望精度結構比最小自由能結構具有更高的精度[14]。

MaxExpect的打分公式如下

(1)

式中:γ取1;Pbp(ij)為堿基i與j相互配對的概率;Pss(k)為堿基k是單鏈的概率。ProbKnot算法中Pbp(i,j)是已知的，本文通過Pbp(i,j)計算出Pss(k)為

(2)

1.2 預測性能指標

預測性能指標是用來評價RNA二級結構預測結果好壞的標準。

馬休茲相互作用系數(Matthews correlation coefficient, MCC)是通過比對天然結構與預測結構而計算出來的值，最小值為-1，表示預測結構與天然結構堿基配對相似度為0；最大值為1，表示預測結構與天然結構相似度為1。 具體衡量公式如下：

(3)

(4)

(5)

式中：TP(true positive)為正確預測堿基對的個數；FN(false negative)為真實結構中存在但沒有被正確預測出的堿基對個數；FP(false positive)為真實結構中不存在卻被錯誤預測到的堿基對個數；TN(true negative)為正確預測的不配對的堿基的個數；敏感性X(sensitivity)指真實結構中所有的堿基對中被正確預測到的百分比；特異性Y(specificity)指在所有預測到的堿基對中正確預測的百分比。

在評價RNA二級結構預測的結果時，如果天然結構中i與j配對，那么預測結構中除了i與j配對，i與j+1或j-1以及j與i-1或i+1配對都被認為是正確的配對，這是考慮到了多種因素，包括確定確切配對的困難性，以及這些預測的配對在質量上與正確配對的一致性[15]。

2 多目標遺傳算法

2.1 多目標優化問題

多目標優化問題的關鍵點在于如何使各個目標之間同時達到最優[16-17]。多目標優化問題的解不僅僅是一個全局最優解，而是一個解集，本文稱這個解集為Pareto最優解集[18]。在這個解集中彼此任意解都不相互支配，并且解集外的解都會被解集中至少一個解支配。

遺傳算法是一種可用于尋找最優解的搜索算法，它借鑒了生物界的進化規律，通過模擬自然進化過程以此來搜索最優解。由于遺傳算法運行一次能找到多目標優化問題的多個Pareto最優解，因而它是求解多目標優化問題的一個有效手段[19]。

NSGA-II是目前多目標優化領域中最優秀的多目標遺傳算法之一，它被國內外學者廣泛引用，本文將主要結合該算法的思想來進行實驗。

2.2 NSGA-II

NSGA-II(Non-dominated sorting genetic algorithm)是一種新型多目標優化遺傳算法，相對于 NSGA 的缺陷，NSGA-II 有如下改進：計算復雜性從O(mN3)降至O(mN2)，具備最優保留機制及無需確定一個共享參數的優點，從而進一步提高計算效率和算法的魯棒性。NSGA-II該算法求得的 Pareto 最優解分布均勻，收斂性和魯棒性好，具有良好的優化效果，是求解多目標優化問題的一種新思路[20-23]。

NSGA-II的遺傳算法中，通過適應度函數計算每個個體的適應度并對其進行排序，排序算法分為兩個步驟。

Step1Non-dominated sorting(非支配排序). 種群基于Pareto順序分為k個子集(Q1,…,Qk)，每個Qi中的個體都不會被Qj中的個體支配(i

Step2Crowding-distance sorting(擁擠距離排序).子群Qi用擁擠距離排序。在不擁擠局域內的解優先度更高。

2.2.1 非支配排序

快速非支配排序是基于非支配排序的改進算法見表1。

非支配排序。m個個體和種群中的其他個體進行支配關系比較，是否支配其他全部個體，復雜度為O(mN)；循環進行直到等級1中的非支配個體全部被搜索到，復雜度為O(mN2)；最壞的情況下，有N個等級，每個等級只存在一個解，復雜度為O(mN3)。

表1 快速非支配排序Table 1 Rapid non dominated sorting

2.2.2 擁擠距離排序

擁擠距離排序用于保持解的多樣性。每個個體的擁擠距離是通過計算與其相鄰的兩個個體在每個子目標函數上的距離差之和來求取,如圖1所示。

圖1 二目標函數的擁擠距離計算Fig.1 The crowding distance calculation of two objective functions

圖1所示個體i的擁擠距離為

Di=(fi+1,1-fi-1,1)+(fi-1,2-fi+1,2),

(6)

即圖1中虛線四邊形的長與寬之和。

排序時當兩個個體屬于不同等級的非支配解集時，等級較小的優先度高，當兩個個體屬于相同等級的非支配解集時，擁擠距離較大的優先度高,即

if(irankjdistance))then(i≥nj)

2.2.3 遺傳算法

NSGA-II中遺傳算法流程如圖2所示。

Step1創造一個初始父代種群P0，使用交叉和變異操作產生子代種群Q0。

Step2對P0和Q0組成的整體R0進行非支配排序，構造所有不同等級的非支配解集Z1,Z2,Z3,……。

Step3對分好等級的非支配解集進行擁擠距離排序，根據適應度高低得到前N個解，構成下一次迭代的父代種群P1。

Step4重復上述3個步驟，直到結果收斂。

圖2 NSGA-II步驟Fig.2 Steps of NSGA-II

3 數據與實驗流程

3.1 數據集

本文采用ASE以及CRW的RNA數據集作為實驗的研究對象，對ASE與CRW的RNA天然結構進行計算，提取出每個RNA結構的一級序列，通過ProbKnot算法的配分函數，得到包含堿基配對概率的pfs文件，將此文件作為本文輸入。經過計算，有效的數據集為877個。其中ASE案例450個，CRW案例427個(見表2)。ASE的序列長度大多數在200～500之間，CRW的序列長度大多數在1 000以上。

表2 實驗數據集表Table 2 The experimental dataset

3.2 實驗過程

基于多目標遺傳算法的含假結RNA二級結構預測方法，實驗過程如下。

Step1首先通過對單個RNA序列進行變異算法得到不同的RNA二級結構作為遺傳算法的初始種群，然后通過對初始種群內的個體進行變異交叉算法形成數量相同的新的RNA二級結構，再對每個RNA二級結構進行適應度計算排序，選擇適應度較高的個體組成新的種群。

Step2重復該步驟直到結果收斂，將種群中的最優個體作為結果輸出，得到這個RNA序列的二級結構預測結果。

Step3輸入其他每個RNA序列來完成上述兩個步驟，得到所有RNA序列的二級結構預測結果，作為實驗結果。其中，為了適應遺傳算法中的交叉與變異，算法對RNA二級結構采用如圖3所示的精英保留策略。圖3中，位于上方的結構(a)表示1號和4號堿基配對、2號和3號堿基為單鏈的RNA二級結構，它通過變異算法可能形成結構(b)與結構(c)。同樣的，位于下方的結構(a)與結構(b)通過交叉算法可能形成結構(c)與結構(d)。

圖3 RNA二級結構的交叉和變異Fig.3 Crossing and mutating of RNA secondary structure

3.3 實驗參數設置

實驗的參數設置會對實驗結果產生影響。本次實驗中有4個參數對實驗結果影響較大，它們是populations，iterations， minloop與α，其中：Populations是初始種群數，iterations是收斂后迭代次數。這兩個參數值設置的越大，遺傳算法搜索到最優解的可能性越高，但運行時間也會越長。根據RNA序列的長度大小，本次實驗populations的取值為1 000～2 000。由于RNA序列長度與種群數的區別，不同的結構使用遺傳算法的迭代次數相差較大，所以iterations表示收斂后的迭代次數而不是總迭代次數，它的取值為500～1 000。

Minloop是最小螺旋長度，α是用于比較配對概率與單鏈概率的一個數值。為了提高交叉變異算法的效率，算法對形成的結構進行了兩方面的約束：一是在每次交叉變異形成RNA二級結構后，算法會去除長度小于minloop的螺旋來保證RNA結構的穩定性，因為天然結構中螺旋長度小于3的結構較為少見，所以minloop取值為3；二是在交叉變異時，只進行概率大于α的配對與解除配對，這樣可以減少結構中配對概率與單鏈概率為0或很小的堿基對與單鏈，避免采用完全隨機交叉變異會消耗大量不必要資源的情況，α取值為0.001～0.1，這個值越大，算法收斂速度越快，但會降低結果的多樣性。

4 結果與分析

4.1 多目標遺傳算法與ProbKnot算法結果對比

在877個實驗結果中，有9個案例(其中ASE案例2個，CRW案例7個)的敏感性和特異性在兩種算法的計算下都為0，因此之后的數據統計中刪除了這9個案例。

圖4、5比較了實驗案例的敏感性與特異性，多目標遺傳算法的結果用黑色的點表示，Probknot算法的結果用灰色的點表示，點所在的位置越是接近右上，結果越好。圖4、5中表明,單從敏感性與特異性的結果上看，多目標遺傳算法的優勢并不明顯。

圖4 448個ASE案例敏感性與特異性值統計Fig.4 The statistic of sensitivity and specificity about 448 cases of ASE

MCC值是綜合敏感性與特異性的評價指標，如圖6、7所示。多目標遺傳算法的結果用黑色的點表示，Probknot算法的結果用灰色的點表示，點所在的位置越是接近上方，結果越好。從圖6、7中可以看出，大多數黑色的點在灰色點的上方，這意味著在大多數實驗案例上，多目標遺傳算法的結果優于Probknot算法。

圖5 420個CRW案例敏感性與特異性值統計Fig.5 The statistic of sensitivity and specificity about 420 cases of CRW

圖6 448個ASE案例MCC值統計Fig.6 The statistic of MCC about 448 cases of ASE

圖7 420個CRW案例MCC值統計Fig.7 The statistic of MCC about 420 cases of CRW

表3整理了ProbKnot算法和多目標遺傳算法的實驗數據的平均值。表3中數據表明，多目標遺傳算法在ASE和CRW數據集上有更高的特異性與MCC值。

表3 RNA二級結構預測方法實驗對比表Table 3 Comparison of methods for predicting RNA secondary structure

4.2多目標遺傳算法與其他RNA二級結構預測算法結果對比

本文選取了263個案例(其中ASE案例158個，CRW案例105個)，采用ProbKnot算法、多目標遺傳算法、MaxExpect算法、Fold算法[24]進行對比實驗。

表4整理了ProbKnot算法、多目標遺傳算法、MaxExpect算法、Fold算法的實驗數據的平均值。表4中數據表明,4種RNA二級結構預測方法中，多目標遺傳算法在ASE和CRW數據集上有更高的特異性與MCC值，在CRW數據集上的特異性也較高。

由于CRW數據集中RNA序列長度相差較大，第2次對比實驗選取了一些長度較短的案例，表3中CRW數據集的總體精度比表2要高。

表4 4種RNA二級結構預測方法實驗對比Table 4 Comparison of four methods for predicting RNA secondary structure

4.3 實驗小結

兩次對比實驗中，多目標遺傳算法在ASE和CRW數據集上的特異性與MCC的平均值更高，大多數案例在改進算法后MCC值得到提高，敏感性提升不明顯。因此本文可以得出結論，多目標遺傳算法的精度更高，得到的結果與天然結構更接近。

5 結 論

1)本文基于多目標遺傳算法，使用概率模型進行RNA二級結構打分，在多目標優化部分計算出堿基的單鏈概率，以此完善概率模型。在遺傳算法部分根據RNA二級結構的特征編寫適合的交叉變異算法，以此提高程序運行效率。

2)實驗結果表明，在ASE和CRW數據集上的對比試驗中，多目標遺傳算法的精度更高，相比ProbKnot算法，MCC值平均提高4%。

3)本文運用多目標優化算法對假結RNA二級結構進行了實驗預測，證明了多目標遺傳算法可以有效地增加含假結RNA二級結構預測的精度，但仍然具有發展和改進的空間。

References)

[1]WAN Yue, QU Kun, ZHANG Qiangfeng, et al. Landscape and variation of RNA secondary structure across the human transcriptome[J]. Nature, 2014, 505(7485): 706-709. DOI:10.1038/nature12946.

[2]PENALVA L O F. RNA secondary structure[M]. New York: Encyclopedia of Systems Biology, 2013: 1864-1864. DOI:10.1007/978-1-4419-9863-7_319.

[3]ZUKER M. Mfold web server for nucleic acid folding and hybridization prediction[J]. Nucleic acids research, 2003, 31(13): 3406-3415.DOI: 10.1093/nar/gkg595.

[4]RUAN Jianhua, STORMO G D, ZHANG Weixiong. An iterated loop matching approach to the prediction of RNA secondary structures with pseudoknots[J]. Bioinformatics, 2004, 20(1): 58-66. DOI: 10.1093/bioinformatics/btg373.

[5]MATHEWS D H, DISNEY M D, CHILDS J L, et al. Incorporating chemical modification constraints into a dynamic programming algorithm for prediction of RNA secondary structure[J]. Proceedings of the National Academy of Sciences of the United States of America, 2004, 101(19): 7287-7292. DOI: 10.1073/pnas.0401799101.

[6]BELLAOUSOV S, MATHEWS D H. ProbKnot: fast prediction of RNA secondary structure including pseudoknots[J]. Rna, 2010, 16(10): 1870-1880.DOI:10.1261/rna.2125310.

[7]BINDEWALD E, KLUTH T, SHAPIRO B A. CyloFold: secondary structure prediction including pseudoknots[J]. Nucleic Acids Research, 2010, 38(suppl_2): W368-W372.DOI:10.1093/nar/gkq432.

[8]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): i85-i93.DOI:10.1093/bioinformatics/btr215.

[9]JABBARI H, CONDON A. A fast and robust iterative algorithm for prediction of RNA pseudoknotted secondary structures[J]. BMC Bioinformatics, 2014, 15: 147. DOI: 10.1186/1471-2105-15-147.

[10]DAWSON W K, TAKAI T, ITO N, et al. A new entropy model for RNA: part III. Is the folding free energy landscape of RNA funnel shaped?[J]. Journal of Nucleic Acids Investigation, 2014, 5: 2652.DOI: 10.4081/jnai.2014.2652.

[11]MATHEWS D H. Using an RNA secondary structure partition function to determine confidence in base pairs predicted by free energy minimization[J]. Rna, 2004, 10(8): 1178-1190. DOI:10.1261/rna.7650904.

[12]SEETIN M G, MATHEWS D H. TurboKnot: rapid prediction of conserved RNA secondary structures including pseudoknots[J]. Bioinformatics, 2012, 28(6): 792-798. DOI: 10.1093/bioinformatics/bts044.

[13]YONEMOTO H, ASAI K, HAMADA M. A semi-supervised learning approach for RNA secondary structure prediction[J]. Computational Biology and Chemistry, 2015, 57: 72-79.DOI:10.1016/j.compbiolchem.2015.02.002.

[14]LU Z J, GLOOR J W, MATHEWS D H. Improved RNA secondary structure prediction by maximizing expected pair accuracy[J]. Rna, 2009, 15(10): 1805-1813. DOI:10.1261/rna.1643609.

[15]DEIGAN K E, LI T W, MATHEWS D H, et al. Accurate SHAPE-directed RNA structure determination[J]. Proceedings of the National Academy of Sciences, 2009, 106(1): 97-102. DOI:10.1073 pnas.0806929106.

[16]DEB K. Multi-objective optimization[M]. Springer US: Search Methodologies, 2014: 403-449.

[17]DEB K. Multi-objective optimization using evolutionary algorithms[M]. New York, NY: John Wiley & Sons, 2001: 3-34. DOI: 10.1007/978-0-85729-652-8_1.

[18]ARNOLD B C. Pareto distribution[M]. New York, NY: John Wiley & Sons, Ltd, 2015. DOI: 10.1002/9781118445112.stat01100.pub2.

[19]LI Hui, ZHANG Qingfu. Multiobjective optimization problems with complicated Pareto sets, MOEA/D and NSGA-II[J]. IEEE Transactions on Evolutionary Computation, 2009, 13(2): 284-302. DOI: 10.1109/TEVC.2008.925798.

[20]DEB K, PRATAP A, AGARWAL S, et al. A fast and elitist multiobjective genetic algorithm: NSGA-II[J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2): 182-197. DOI: 10.1109/4235.996017.

[21]DEB K, AGRAWAL S, PRATAP A, et al. A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II[C]//International Conference on Parallel Problem Solving From Nature. Berlin Heidelberg: Springer, 2000: 849-858. DOI: 10.1007/3-540-45356-3_83.

[22]DEB K, JAIN H. Handling many-objective problems using an improved NSGA-II procedure[C]//2012 IEEE Congress on Evolutionary Computation. Brisbane, QLD, Australia: IEEE, 2012: 1-8.DOI: 10.1109/CEC.2012.6256519

[23]SEADA H, DEB K. U-nsga-III: A unified evolutionary optimization procedure for single, multiple, and many objectives: Proof-of-principle results[C]//International Conference on Evolutionary Multi-Criterion Optimization. Cham: Springer, 2015: 34-49. DOI: 10.1007/978-3-319-15892-1_3.

[24]BELLAOUSOV S, REUTER J S, SEETIN M G, et al. RNAstructure: Web servers for RNA secondary structure prediction and analysis[J]. Nucleic acids research, 2013, 41:W471-W474.DOI:10.1093/nar/gkt290.

PredictionofRNAsecondarystructureincludingpseudoknotsbasedonmulti-objectivegeneticalgorithm

GU Ti1，CAI Leixin1，WANG Shuai1，Lü Qiang2*

(1.SchoolofComputerScience&Technology,SoochowUniversity,Suzhou215006,China;2.ProvincialKeyLaboratoryforComputerInformationProcessingTechnologyofJiangsu(SoochowUniversity),Suzhou215006,China)

Pseudoknot is a kind of important structure in RNA, which is more difficult to be predicted due to the difficulty in training the prediction model. ProbKnot algorithm has a high accuracy in predicting the secondary structure of pseudoknotted RNA based on base-pair probabilities. However, this algorithm only makes use of the base-pair probabilities as a predictive feature, which leads to a large number of negative pairs and lower specificity. The overall accuracy can be improved which follows the improvement of specificity by combining the probabilistic model of base pairing in ProbKnot algorithm and multi-objective genetic algorithm. In the experiments, we will first calculate the single-stranded probability as a new predictive feature, and then use genetic algorithm to cross, mutate and iterate the secondary structure of RNA. Finally, we can get the Pareto optimal solutions and the highest maximum expected accuracy. The experiment results showed that in the RNA cases, applying this method could achieve an average increase of 4% in the accuracy compared with the current existing methods.

RNA secondary structure; Pseudoknot; Multi-objective optimization; Genetic algorithm; Maximum expected accuracy

TP391

A

1672-5565(2017)03-142-07

10.3969/j.issn.1672-5565.201701006

2017-01-20；

2017-04-11.

國家自然科學基金(61170125).

顧倜，男，碩士研究生，研究方向：生物信息計算；E-mail：20145227032@stu.suda.edu.cn.

*通信作者：呂強，男，教授，研究方向：生物信息計算、元啟發搜索、并行計算；E-mail：qiang@suda.edu.cn.