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

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

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

2.2 自適應(yīng)序列長度的模型框架

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

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

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

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

2.3 自適應(yīng)序列長度的LSTM訓(xùn)練算法

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

1)算法輸入:模型的輸入為RNA一級(jí)序列，經(jīng)過序列的擴(kuò)充及配分函數(shù)計(jì)算后，得到網(wǎng)絡(luò)輸入特征x，x為三維數(shù)組，x［i，j，k］表示 batch_size條序列中第 i條序列的第 j個(gè)堿基的第k個(gè)特征，N+5維向量x［i，j］為該堿基的輸入特征向量，N為補(bǔ)全后的序列長度，x［i，j］的第1個(gè)元素為該類堿基在第i條序列中出現(xiàn)的頻率，第2-5個(gè)元素為堿基的類型，即［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個(gè)元素為配分函數(shù)的輸出.

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

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

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

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

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

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

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

評(píng)價(jià)指標(biāo)為準(zhǔn)確率 acc［24］、敏感性 sen、特異性 ppv［23］和Matthews 相關(guān)系數(shù) mcc［25］.

3 實(shí)驗(yàn)與分析

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

本文使用的數(shù)據(jù)集來自權(quán)威數(shù)據(jù)集RNA STRAND［23］，從中選取五個(gè)子集:TMR、SPR、SRP、RFA、和 ASE.五個(gè)數(shù)據(jù)集共有2493條RNA序列，其中1158條序列包含假結(jié)，假結(jié)總數(shù)為3116，平均長度267.37，其中最長的序列有553個(gè)堿基.TMR、SPR、SRP、RFA和ASE五個(gè)數(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 定長模型與自適應(yīng)模型比較

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

網(wǎng)絡(luò)結(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］代替.模型的輸出類別為“配對(duì)”和“不配對(duì)”.

表2 自適應(yīng)LSTM與定長LSTM實(shí)驗(yàn)結(jié)果Table 2 Experiment results of adaptive LSTM and fixed LSTM

定長LSTM模型要求對(duì)數(shù)據(jù)進(jìn)行截?cái)嗵幚恚瑢⑿蛄薪財(cái)嘀翑?shù)據(jù)集中最短序列的長度，對(duì)于ASE，SPR和TMR數(shù)據(jù)集，分別截取長度為189，54，102的序列進(jìn)行訓(xùn)練的測試.

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

表2所為自適應(yīng)LSTM(adaptive)與定長LSTM(fixed)在五個(gè)數(shù)據(jù)集上的實(shí)驗(yàn)結(jié)果，圖2所示為mcc和acc指標(biāo)的散點(diǎn)圖，其中“·”表示自適應(yīng)LSTM，“+”表示定長LSTM.

圖2 mcc和acc散點(diǎn)圖Fig.2 mcc and acc scatter plot

從實(shí)驗(yàn)結(jié)果中可以看出，自適應(yīng)LSTM方法在5個(gè)數(shù)據(jù)集的mcc和acc指標(biāo)均高于定長LSTM方法.其原因在于LSTM方法對(duì)序列進(jìn)行截?cái)嗟念A(yù)處理，將數(shù)據(jù)集中所有序列截?cái)嘀料嗤L度，使得部分原本配對(duì)的堿基在截?cái)嗪笞兂刹慌鋵?duì)，從而丟失部分特征信息.

3.3 與其它四種經(jīng)典方法對(duì)比實(shí)驗(yàn)

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

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

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

4 總結(jié)

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

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

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