一種自適應序列長度的RNA二級結(jié)構(gòu)深度預測方法

吳宏杰，湯 燁，陸衛(wèi)忠，崔志明，2，付保川，GAO Zhen，3(蘇州科技大學電子與信息工程學院，江蘇蘇州25009)

2(蘇州大學智能信息處理及應用研究所，江蘇蘇州215006)

3(School of Engineering Technology McMaster University，Hamilton Ontario Canada 45011)

E-mail:luwz@mail.usts.edu.cn

1 引言

RNA(核糖核苷酸)是遺傳信息的關鍵載體，其二維結(jié)構(gòu)與基因的成熟、調(diào)節(jié)和功能密切相關［1］.研究RNA功能與結(jié)構(gòu)的關系，確定RNA折疊的形式和頻率對于揭示RNA分子在生命過程中的作用具有重要意義［2，3］RNA分子的二級結(jié)構(gòu)描述了堿基之間的相互作用關系.確定RNA的二級結(jié)構(gòu)不僅是RNA分子從一級結(jié)構(gòu)到三級結(jié)構(gòu)的必經(jīng)之路，可以有效縮小RNA三級結(jié)構(gòu)的解空間，降低搜索復雜度;而且二級結(jié)構(gòu)也是RNA功能推斷的重要依據(jù).

從計算機角度，將長度為n的RNA堿基序列視為字符序列.序列中的每個堿基用i(0≤i≤n－1)作為下標，二級結(jié)構(gòu)中的兩個相互作用的堿基稱為堿基對(i，j)［5］.我們把 RNA的堿基序列稱為一級結(jié)構(gòu)，把序列中的堿基對信息稱為二級結(jié)構(gòu).從RNA一級結(jié)構(gòu)推斷其相應二級結(jié)構(gòu)的過程稱為二級結(jié)構(gòu)預測.確定RNA的二級結(jié)構(gòu)后，該RNA在二維平面中的實際排列圖形也就隨之確定.

假結(jié)是指序列中存在交叉的堿基對(i，j)和(k，l)且i＜k＜j＜l，它是RNA二級結(jié)構(gòu)的一種子結(jié)構(gòu)，在生命科學鄰域中具有重要意義［4］.不帶假結(jié)的RNA二級結(jié)構(gòu)預測問題，已被 Zuker［5］和 Mathews［6，7］的 動 態(tài) 規(guī) 劃 算 法 及 其 改 進 算 法mfold［8］，GTfold［9］等充分研究，其時間復雜度是 O(n3)，空間復雜度是O(n2).帶假結(jié)的RNA二級結(jié)構(gòu)預測問題是一個NP完全的優(yōu)化問題［10］.依據(jù)組合數(shù)學方法，對于長度為n的RNA序列可以折疊出1.8n種可能的RNA二級結(jié)構(gòu)，隨著序列長度n的增長，其計算復雜度將快速增長［4］.

常見的RNA二級結(jié)構(gòu)預測模型可以分為熱力學模型、同源對比模型和統(tǒng)計學習模型三類.在熱力學模型中，假設RNA分子服從熱力學規(guī)律，自由能越低的RNA結(jié)構(gòu)越穩(wěn)定.預測方法從所有可能的候選結(jié)構(gòu)中搜索出自由能最小的結(jié)構(gòu)，作為預測結(jié)果.同源對比模型通過從同源序列中尋找共突變的堿基對進行結(jié)構(gòu)預測，需要找到同源序列作為輸入［11］.統(tǒng)計學習模型將RNA二級結(jié)構(gòu)預測問題轉(zhuǎn)化為序列中堿基的配對結(jié)果的分類問題，學習已知的RNA二級結(jié)構(gòu)的規(guī)律，構(gòu)建學習模型，其預測準確率一般都能超過前兩類方法.本文主要討論統(tǒng)計學習模型.

ProbKnot［12］是RNA structure軟件中基于堿基配對概率的統(tǒng)計學習預測方法，通過計算堿基之間的配對概率進行預測，能夠預測出所有拓撲結(jié)構(gòu)的假結(jié).Cylofold［13］是一種不受假結(jié)復雜程度限制的RNA二級結(jié)構(gòu)預測方法，它是以莖區(qū)為單位的計算型方法.CentroidFold［14］是一種基于統(tǒng)計決策理論的后驗解碼方法，用于預測RNA二級結(jié)構(gòu)，屬于統(tǒng)計學習模型，但只能預測長度不大于400的序列且預測結(jié)果不包含假結(jié).IPknot是一種基于整數(shù)規(guī)劃的快速RNA二級結(jié)構(gòu)預測方法［15］.

上述方法及常用統(tǒng)計方法取得了一定的效果，但仍然存在三個難點:第一，RNA堿基配對信息是生物序列長期不斷進化形成的結(jié)果，進化過程具有生物多樣性，無法用簡單的公式或規(guī)律進行描述.導致大多數(shù)的淺層學習方法難以挖掘序列中隱含的深層序列信息;第二，這些學習方法的經(jīng)典模型只能接受固定大小的特征向量，而在面向序列的預測問題中，特征向量往往源自于序列信息，特征向量的大小隨著序列長度的變化而變化.大多數(shù)方法將序列截斷至相同長度后進行訓練，但是這樣不僅丟失了截斷部分的信息，破壞了序列的整體性，而且部分原本配對的堿基由于與之配對的堿基被截斷，變?yōu)椴慌鋵A基，嚴重影響訓練集的準確性;第三，在序列的真實結(jié)構(gòu)中，配對與不配對堿基分布不均衡導致訓練數(shù)據(jù)集合不平衡.

長短期記憶 LSTM［17］(Long Short-Term Memory)是一種典型的深度神經(jīng)網(wǎng)絡模型，其遞歸神經(jīng)網(wǎng)絡模型結(jié)構(gòu)，能有效挖掘數(shù)據(jù)中深層隱含特征［17，18］.近幾年，在自然語言處理、圖像識別、語音識別領域都有了出色較好的表現(xiàn).本文嘗試用LSTM來建模RNA二級結(jié)構(gòu)模型.更進一步，經(jīng)典的深度神經(jīng)網(wǎng)絡同樣受限于固定長度的特征輸入，難以充分發(fā)揮深層挖掘能力的優(yōu)勢.由此，本文在經(jīng)典長短期記憶LSTM網(wǎng)絡的基礎上，建立了一種自適應長度的模型框架，引入了Mask向量，建立了相應的訓練算法，從而避免了長序列被截斷問題，能夠自適應RNA序列長度.同時引入了加權(quán)向量，動態(tài)調(diào)節(jié)序列中各個堿基損失函數(shù)占序列總損失函數(shù)的比重，緩解了樣本分布不均衡的現(xiàn)象.在訓練特征方面，除了常用特征外，新引入了配分函數(shù)的輸出矩陣作為訓練特征.Mathews［19］的研究表明，配分函數(shù)的輸出矩陣與堿基配對概率顯著相關.本方法的預測結(jié)果與四種經(jīng)典方法ProbKnot，Cylofold，CentroidFold和IPknot比較，acc和 mcc指標分別提升了17.9%和14.7%以上.本方法在五個經(jīng)典數(shù)據(jù)集上進行了驗證，并提供了源碼、數(shù)據(jù)集與結(jié)果.

2 基于自適應序列長度模型的RNA二級結(jié)構(gòu)預測方法

2.1 序列長度固定的LSTM網(wǎng)絡

基于LSTM的RNA二級結(jié)構(gòu)預測模型的結(jié)構(gòu)如圖1右圖所示(圖中省略了各個功能層之間用于數(shù)據(jù)形狀轉(zhuǎn)化的reshape層)，此模型能預測出假結(jié)且允許堿基缺失，模型由輸入層、LSTM層、全連接層和輸出層組成.LSTM層用于提取RNA序列特征，全連接層進行特征的進一步提取和分類.模型的輸入特征向量xI代表序列中第i個堿基的特征向量.輸入層將n個輸入特征向量映射為n個d維向量后分別輸入LSTM層中的前向LSTM單元lstm_f和反向LSTM單元lstm_b，將兩者的輸出拼接后得到n個2d維向量，再輸入由全連接層和輸出層組成的BP神經(jīng)網(wǎng)絡，進行分類.全連接層使用tanh作為激活函數(shù)，同時加入了dropout層以提高模型對測試數(shù)據(jù)的擬合程度.輸出層使用softmax作為激活函數(shù)，將全連接層的輸出轉(zhuǎn)化為概率分布向量y［i，j］，用于表示序列中每個堿基屬于各個輸出類別的概率.

模型使用分類問題中常用的交叉熵作為損失函數(shù)，對每一條輸入模型的序列，它的總損失函數(shù)等于所有堿基預測結(jié)果損失函數(shù)的均值，其計算公式如下:

其中n為序列長度，m為類別個數(shù)，y和y'均為二維數(shù)組，分別表示模型預測結(jié)果和樣本實際標簽，y［i，j］表示第i個堿基預測屬于第j個類別的概率.得到序列的總損失函數(shù)后，使用批量梯度下降算法，更新網(wǎng)絡模型中的權(quán)值和偏置參數(shù)，使總損失函數(shù)最小.在此LSTM模型的基礎上，為了能適應序列長度變化，本文增加了自適應模塊，通過Mask向量來控制序列長度.

2.2 自適應序列長度的模型框架

本文設計的自適應序列長度的LSTM模型框架如圖1所示，右側(cè)為LSTM模塊，左側(cè)為自適應模塊，包括加權(quán)向量和mask向量的計算.

本模型的自適應模塊以RNA一級序列為輸入，對于第i條序列Si首先提取出序列中堿基類型等特征，生成序列的特征矩陣Xi，然后進行預處理，將Xi擴充至指定大小，同時為每條序列計算出它的Mask向量Mi和權(quán)值向量Wi，Mask向量用于控制序列長度和輸出類別個數(shù)，權(quán)值向量用于控制序列中各個堿基損失函數(shù)占序列總損失函數(shù)的比重，再將序列特征輸入訓練好的LSTM模型，得到序列中每個堿基屬于各個類別的概率，最后對模型的輸出進行處理，便可得到配對結(jié)果，即RNA序列的二級結(jié)構(gòu).

在計算第k條序列中第i個堿基的交叉熵時，加入了系數(shù)Mk(i)，Mk(i)也會以乘積的形式出現(xiàn)在計算梯度公式中，當該堿基為原序列中的堿基(即i≤nk)時，Mk(i)=1，梯度值不變，網(wǎng)絡權(quán)值更新與傳統(tǒng)方法相同，當該堿基不是原序列中的堿基(即nk＜i≤N)時，Mk(i)=0，梯度值為0，網(wǎng)絡權(quán)值不會更新. 因此，模型對擴展堿基的無效預測結(jié)果不會影響模型的更新.

圖1 序列長度自適應的模型框架圖Fig.1 Adaptive sequence length model framework

2.3 自適應序列長度的LSTM訓練算法

在RNA二級結(jié)構(gòu)預測的多分類問題中，不僅要預測每個堿基是否配對，還要詳細預測出與之配對的堿基序號，需要模型能夠提取序列中更深層次的特征信息，因此將圖1模型中的全連接層的隱藏層數(shù)增加到二層，與雙向LSTM層相連的隱藏層神經(jīng)元個數(shù)為2×d_fc，與輸出層相連的隱藏層神經(jīng)元個數(shù)為d_fc.同時在第i堿基的輸入特征中加入堿基頻率和配分函數(shù).配分函數(shù)可由RNA Structure軟件中的pfunction模塊計算.

1)算法輸入:模型的輸入為RNA一級序列，經(jīng)過序列的擴充及配分函數(shù)計算后，得到網(wǎng)絡輸入特征x，x為三維數(shù)組，x［i，j，k］表示 batch_size條序列中第 i條序列的第 j個堿基的第k個特征，N+5維向量x［i，j］為該堿基的輸入特征向量，N為補全后的序列長度，x［i，j］的第1個元素為該類堿基在第i條序列中出現(xiàn)的頻率，第2-5個元素為堿基的類型，即［0，0，0，0］，［0，0，0，1］，［0，0，1，0］，［0，1，0，0］和［1，0，0，0］五類，分別表示缺失堿基，A，G，C和U，后N個元素為配分函數(shù)的輸出.

2)算法輸出:模型的輸出 Y 為二維數(shù)組，y［i，j］表示batch_size條序列中第i條序列的第j個堿基預測類別，其值為0表示不配對，不為0時表示第j個堿基與第y［i，j］個堿基配對.

3)訓練算法:在多分類問題中使用了動態(tài)加權(quán)方法的自適應LSTM模型訓練算法如下:

其中T為迭代次數(shù)上限，N為最大序列長度，n(i)為第i條序列的長度，x為輸入序列特征，y和y'分別預測結(jié)果和實際標簽，M和W分別為Mask向量和權(quán)值向量，運算符“·”向量中對應位置元素的乘積運算，get_data()函數(shù)用于從數(shù)據(jù)集抽取數(shù)據(jù)，preprocess(x，n［i］，N)用于將原始序列補全至長度為 N，并生成 Mask向量和權(quán)值向量，lstm_f()和

2.4 網(wǎng)絡模型參數(shù)

實驗以python作為編程語言，在TensorFlow框架下進行模型的搭建，訓練和測試.采用五折交叉驗證，訓練集與測試集比例為1:4.模型中各個功能層的輸入張量和輸出張量形狀如表1所示.

表1中batch_size為批量梯度下降時的每批樣本包含的序列個數(shù);d_in為序列中每個堿基的輸入特征長度，d_step為LSTM中的循環(huán)次數(shù)，即最大序列長度;d_lstm為LSTM中隱藏層神經(jīng)元個數(shù);d_fc為全連接層神經(jīng)元個數(shù);d_out為輸出類別個數(shù).

表1 各層張量形狀表Table 1 Tensor shape of layers

評價指標為準確率 acc［24］、敏感性 sen、特異性 ppv［23］和Matthews 相關系數(shù) mcc［25］.

3 實驗與分析

3.1 數(shù)據(jù)集

本文使用的數(shù)據(jù)集來自權(quán)威數(shù)據(jù)集RNA STRAND［23］，從中選取五個子集:TMR、SPR、SRP、RFA、和 ASE.五個數(shù)據(jù)集共有2493條RNA序列，其中1158條序列包含假結(jié)，假結(jié)總數(shù)為3116，平均長度267.37，其中最長的序列有553個堿基.TMR、SPR、SRP、RFA和ASE五個數(shù)據(jù)集中包含的序列總數(shù)分別為 721、622、383、313和454，序列長度范圍分別為102～463、54～93、66～533、40～553和189～486，包含假結(jié)的序列總數(shù)分別為 713、0、0、29 和 416.

3.2 定長模型與自適應模型比較

通過比較定長LSTM模型和自適應LSTM模型在RNA二級結(jié)構(gòu)預測的二分類問題上的預測效果，即僅預測堿基是否配對，以說明自適應LSTM模型的有效性.

網(wǎng)絡結(jié)構(gòu)采用圖1右圖中的模型結(jié)構(gòu).模型的輸入特征為堿基類型，用［0，0，0，1］，［0，0，1，0］，［0，1，0，0］和［1，0，0，0］分別代表 A，G，C 和 U 四種堿基，缺失堿基用［0，0，0，0］代替.模型的輸出類別為“配對”和“不配對”.

表2 自適應LSTM與定長LSTM實驗結(jié)果Table 2 Experiment results of adaptive LSTM and fixed LSTM

定長LSTM模型要求對數(shù)據(jù)進行截斷處理，將序列截斷至數(shù)據(jù)集中最短序列的長度，對于ASE，SPR和TMR數(shù)據(jù)集，分別截取長度為189，54，102的序列進行訓練的測試.

兩種模型中LSTM隱藏層神經(jīng)元個數(shù)d_lstm設置為16，全連接層神經(jīng)元個數(shù)d_fc設置為8，學習率為0.01，batch_size為10，迭代次數(shù)上限為600倍訓練集樣本數(shù).

表2所為自適應LSTM(adaptive)與定長LSTM(fixed)在五個數(shù)據(jù)集上的實驗結(jié)果，圖2所示為mcc和acc指標的散點圖，其中“·”表示自適應LSTM，“+”表示定長LSTM.

圖2 mcc和acc散點圖Fig.2 mcc and acc scatter plot

從實驗結(jié)果中可以看出，自適應LSTM方法在5個數(shù)據(jù)集的mcc和acc指標均高于定長LSTM方法.其原因在于LSTM方法對序列進行截斷的預處理，將數(shù)據(jù)集中所有序列截斷至相同長度，使得部分原本配對的堿基在截斷后變成不配對，從而丟失部分特征信息.

3.3 與其它四種經(jīng)典方法對比實驗

本實驗比較自適應LSTM(adaptive)與ProbKnot算法、Cylofold算法、CentroidFold算法和IPknot算法在TMR，SPR，SRP，RFA，ASE五個的數(shù)據(jù)集上的預測結(jié)果.

模型的輸入輸出與2.3節(jié)相同，實驗使用五折交叉驗證，d_steps為 560，d_lstm 為 128，學習率為 0.01，batch_size 為10，迭代次數(shù)上限為2000倍訓練集樣本數(shù).訓練完成后分別統(tǒng)計模型在五個數(shù)據(jù)集上的各項評價指標.

實驗結(jié)果如下頁表3所示，從表3中可以看出，自適應LSTM算法的acc指標在5個數(shù)據(jù)集上均高于其他算法，acc和mcc的均值分別比其他方法高13.6%和14.8%.在序列平均長度最大，含假結(jié)序列最多的TMR數(shù)據(jù)集上優(yōu)勢最明顯，分別提高了18.7%和38.9%.

4 總結(jié)

本文在遞歸神經(jīng)網(wǎng)絡的基礎上，提出了一種自適應序列長度的RNA二級結(jié)構(gòu)預測深度學習方法，引入了Mask向量，構(gòu)建了新的損失函數(shù)，形成了相應的訓練算法.本方法較好的解決了經(jīng)典模型不支持可變序列數(shù)據(jù)與序列長度可變的兩難問題.同時加入了動態(tài)加權(quán)算法，緩解了樣本分布不均衡的問題.實驗結(jié)果表明，該方法能有效提高RNA二級結(jié)構(gòu)預測的精度，在五個測試集上，比其他四種典型方法提高了13.6%和14.8%;在序列平均長度最大，假結(jié)最多的TMR數(shù)據(jù)集上優(yōu)勢最明顯，分別提高了18.7%和38.9%.

本文進一步的研究方向包括:第一，此方法的預測結(jié)果中可能存在一些不符合配對原則的配對結(jié)果，使用合理的算法對神經(jīng)網(wǎng)絡的輸出結(jié)果進行篩選處理.第二，改進網(wǎng)絡結(jié)構(gòu)，以進一步提高模型的預測精度.

表3 與其它四種經(jīng)典方法對比實驗結(jié)果Table 3 Experiment results of four classic methods