一種自適應序列長度的RNA二級結構深度預測方法

吳宏杰，湯 燁，陸衛忠，崔志明，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(核糖核苷酸)是遺傳信息的關鍵載體，其二維結構與基因的成熟、調節和功能密切相關［1］.研究RNA功能與結構的關系，確定RNA折疊的形式和頻率對于揭示RNA分子在生命過程中的作用具有重要意義［2，3］RNA分子的二級結構描述了堿基之間的相互作用關系.確定RNA的二級結構不僅是RNA分子從一級結構到三級結構的必經之路，可以有效縮小RNA三級結構的解空間，降低搜索復雜度;而且二級結構也是RNA功能推斷的重要依據.

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

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

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

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

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

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

2 基于自適應序列長度模型的RNA二級結構預測方法

2.1 序列長度固定的LSTM網絡

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

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

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

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

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

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

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

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

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

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

1)算法輸入:模型的輸入為RNA一級序列，經過序列的擴充及配分函數計算后，得到網絡輸入特征x，x為三維數組，x［i，j，k］表示 batch_size條序列中第 i條序列的第 j個堿基的第k個特征，N+5維向量x［i，j］為該堿基的輸入特征向量，N為補全后的序列長度，x［i，j］的第1個元素為該類堿基在第i條序列中出現的頻率，第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個元素為配分函數的輸出.

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

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

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

2.4 網絡模型參數

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

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

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

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

3 實驗與分析

3.1 數據集

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

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

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

網絡結構采用圖1右圖中的模型結構.模型的輸入特征為堿基類型，用［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實驗結果Table 2 Experiment results of adaptive LSTM and fixed LSTM

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

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

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

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

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

3.3 與其它四種經典方法對比實驗

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

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

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

4 總結

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

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

表3 與其它四種經典方法對比實驗結果Table 3 Experiment results of four classic methods