多特征融合的lncRNA識(shí)別與其功能預(yù)測(cè)

常征，孟軍，施云生，莫馮然

（大連理工大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，遼寧 大連 116023）

近年來(lái)，非編碼RNA(non-coding RNA,ncRNA)識(shí)別的相關(guān)研究已成為人們關(guān)注的熱點(diǎn)。一直以來(lái)，轉(zhuǎn)錄本被大家普遍認(rèn)為只起到翻譯蛋白質(zhì)的作用，但隨著人類(lèi)基因組注釋工作的不斷推進(jìn)，研究結(jié)果表明只有大約1%～2%的基因參與了編碼蛋白的工作[1]，而以往被大家忽略的非編碼序列也在整個(gè)生命活動(dòng)中扮演著至關(guān)重要的角色。這些非編碼序列中，有一種長(zhǎng)度大于200 nt、無(wú)法編碼蛋白質(zhì)的轉(zhuǎn)錄本尤其受到關(guān)注，被稱(chēng)為長(zhǎng)鏈非編碼RNA(long non-coding RNA,lncRNA)[2]。近年來(lái)發(fā)現(xiàn)lncRNA具有調(diào)節(jié)生物體生命活動(dòng)的重要作用[3-4]，而各種傳統(tǒng)的實(shí)驗(yàn)方法，一方面需要花費(fèi)大量時(shí)間和高額費(fèi)用，另一方面，因?yàn)閘ncRNA的低表達(dá)和低保守性等原因，在識(shí)別lncRNA方面受到不同程度的影響。研究人員對(duì)人和動(dòng)物進(jìn)行了大量的實(shí)驗(yàn)，并且出現(xiàn)了具有良好魯棒性的lncRNA識(shí)別軟件。

RNAseq和全基因組陣列分析顯示，植物體內(nèi)也存在大量的lncRNA，它們?cè)谥参锏拈_(kāi)花、雄性不育、營(yíng)養(yǎng)代謝、生物和非生物脅迫等生物過(guò)程中起著調(diào)節(jié)因子的作用[5]。與哺乳動(dòng)物相比，植物ncRNA的研究起步比較晚，且多數(shù)集中在短鏈非編碼RNA上，這為植物lncRNA識(shí)別與分析帶來(lái)了困難。研究植物lncRNA將幫助生命學(xué)科的工作者進(jìn)一步揭示植物內(nèi)部生命活動(dòng)，因此深入研究植物lncRNA并預(yù)測(cè)其功能具有非常重要的意義。

目前，在計(jì)算預(yù)測(cè)lncRNA方面，許多研究工作都利用機(jī)器學(xué)習(xí)算法建立預(yù)測(cè)模型，通過(guò)輸入各類(lèi)序列特征、結(jié)構(gòu)特征，構(gòu)建識(shí)別lncRNA的分類(lèi)器模型。研究表明，對(duì)于lncRNA識(shí)別，通過(guò)提取開(kāi)放閱讀框、密碼子頻率偏好性、與已知蛋白質(zhì)相似度等特征作為輸入，對(duì)線性回歸、支持向量機(jī)以及其他模型進(jìn)行訓(xùn)練得到的分類(lèi)器具有良好的分類(lèi)效果[6]。近年來(lái)衍生出的預(yù)測(cè)軟件多采用以上特征。其中，CPC[7]和 CPAT[8]都是通過(guò)序列特征來(lái)區(qū)分編碼和非編碼RNA；CNCI[9]能夠?qū)⒂?xùn)練好的分類(lèi)器運(yùn)用到近親物種的lncRNA識(shí)別；PLEK[10]可以從高通量測(cè)序的轉(zhuǎn)錄本中識(shí)別lncRNA。 然而，大多數(shù)軟件只在動(dòng)物數(shù)據(jù)集上得到良好的驗(yàn)證，專(zhuān)門(mén)為植物lncRNA識(shí)別設(shè)計(jì)的軟件目前還比較稀缺。

隨著基因組學(xué)研究的不斷深入，產(chǎn)生了大量未被標(biāo)注的基因序列。由于生物實(shí)驗(yàn)方法驗(yàn)證基因功能的代價(jià)十分昂貴，如何通過(guò)計(jì)算機(jī)方法對(duì)基因序列功能進(jìn)行大規(guī)模預(yù)測(cè)成了近年來(lái)生物信息學(xué)的研究熱點(diǎn)之一[11]。

為了進(jìn)一步提高植物lncRNA預(yù)測(cè)的準(zhǔn)確性，基于機(jī)器學(xué)習(xí)分類(lèi)算法，通過(guò)對(duì)下載的高可信度數(shù)據(jù)提取開(kāi)放閱讀框、k核苷酸頻率以及二級(jí)結(jié)構(gòu)特征等多特征融合[12]作為輸入特征，訓(xùn)練樸素貝葉斯、支持向量機(jī)和梯度提升決策樹(shù)3種分類(lèi)模型，并采用加權(quán)投票的多分類(lèi)器集成方法，集成分類(lèi)結(jié)果以得到更好的分類(lèi)性能。利用標(biāo)注測(cè)序數(shù)據(jù)測(cè)試集對(duì)模型進(jìn)行驗(yàn)證、分析并選擇性能最好的將其作為最終分類(lèi)器。提出的方法通過(guò)五折交叉驗(yàn)證，得到了較好的性能。在功能預(yù)測(cè)方面，根據(jù)lncRNA-microRNA、microRNA-mRNA相互作用關(guān)系，建立調(diào)控網(wǎng)絡(luò)，利用相關(guān)聯(lián)的RNA預(yù)測(cè)lncRNA的功能。

1 基于多特征融合的lncRNA預(yù)測(cè)

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

擬南芥的生物學(xué)實(shí)驗(yàn)數(shù)據(jù)和基因注釋信息相對(duì)比較豐富，常被廣泛用于植物脅迫響應(yīng)的研究中[13]。本文使用的正集數(shù)據(jù)為PNRD[14](http://structuralbiology.cau.edu.cn/PNRD/) 2 565條具有高可信度的擬南芥lncRNA序列。負(fù)集數(shù)據(jù)是從Ref-Seq數(shù)據(jù)庫(kù)下載的48 148條mRNA序列。為了保證正負(fù)樣本均衡，從負(fù)集原始數(shù)據(jù)中隨機(jī)采樣出2 500條mRNA作為最終訓(xùn)練集，如表1所示。

表 1 數(shù)據(jù)集信息Table 1 Dataset information

1.2 開(kāi)放閱讀框

在分子遺傳學(xué)中，開(kāi)放閱讀框(open reading frame，ORF)是閱讀框的一部分，具有潛在的翻譯能力[15]。研究表明，mRNA的ORF覆蓋率明顯高于lncRNA，且mRNA具有更多的完整性開(kāi)放閱讀框[16]。首先利用TransDecoder軟件包計(jì)算得到每個(gè)序列的開(kāi)放閱讀框信息，然后對(duì)其分別提取完整性、ORF覆蓋率以及歸一化的ORF值3種特征。將完整性定義為一個(gè)布爾變量bool，0代表不存在完整性O(shè)RF，1表示存在完整性O(shè)RF。覆蓋率Coverage等于所有的合法ORF(本文只考慮正鏈情況下)的長(zhǎng)度與RNA序列長(zhǎng)度之比，歸一化ORF Normalized_ORF是序列中ORF個(gè)數(shù)n與RNA序列長(zhǎng)度L的比值，分別定義為

式中Li代表序列中第i個(gè)ORF的長(zhǎng)度。

整合3種特征得到特征向量：

1.3 二級(jí)結(jié)構(gòu)

二級(jí)結(jié)構(gòu)(secondary structures，SS)是單條序列通過(guò)堿基配對(duì)自身形成莖區(qū)和環(huán)球，與RNA的功能息息相關(guān)，可以作為識(shí)別lncRNA的重要依據(jù)。目前，預(yù)測(cè)二級(jí)結(jié)構(gòu)的計(jì)算方法分為兩類(lèi)：基于熱力學(xué)和基于系統(tǒng)發(fā)生學(xué)。前者認(rèn)為：生物體在形成高維結(jié)構(gòu)的時(shí)候，將使自身達(dá)到穩(wěn)態(tài)結(jié)構(gòu)，因此釋放的能量應(yīng)更多。研究表明，二級(jí)結(jié)構(gòu)越是穩(wěn)定(釋放的自由能越多)其潛在的編碼能力越強(qiáng)。同時(shí)，二級(jí)結(jié)構(gòu)的穩(wěn)定性與RNA序列中配對(duì)堿基個(gè)數(shù)以及核苷酸C和G的含量有一定的關(guān)系。本文使用ViennaRNA[17]工具包對(duì)序列形成二級(jí)結(jié)構(gòu)釋放的自由能進(jìn)行計(jì)算并得到二級(jí)結(jié)構(gòu)的點(diǎn)括號(hào)表示形式，然后從中提取出配對(duì)堿基的個(gè)數(shù)以及C和G堿基的含量。歸一化最小自由能由如下：

式中：MFE是釋放的自由能，L是RNA序列的長(zhǎng)度。

整合上述3個(gè)特征得到如下特征向量：

式中np為配對(duì)堿基的個(gè)數(shù)。

1.4 k-mers

密碼子是遺傳物質(zhì)編碼的信息由活細(xì)胞轉(zhuǎn)化為蛋白質(zhì)的一套規(guī)則，蛋白質(zhì)為保證其某些生物功能，在自然選擇下會(huì)表現(xiàn)出對(duì)某些密碼子的偏好性。因此mRNA在密碼子方面表現(xiàn)出一定的保守性，而不編碼蛋白的lncRNA其保守性較差。所以可以使用密碼子頻率當(dāng)作識(shí)別lncRNA的一個(gè)特征。然而，因?yàn)闊o(wú)法準(zhǔn)確定位mRNA編碼區(qū)域，且lncRNA有多個(gè)編碼區(qū)域，直接計(jì)算密碼子頻率存在一定困難。為解決以上問(wèn)題，使用一個(gè)近似的解決方法：k-mers特征計(jì)算。

一個(gè)k-mer具有k個(gè)核苷酸，每個(gè)核苷酸可以是 A、C、G 或 T，k取值為 1、2和 3，則有 4+16+64=84種模式：4個(gè) 1-mer，16個(gè) 2-mer，64個(gè)3-mer。使用一個(gè)長(zhǎng)度為k的滑動(dòng)窗口來(lái)匹配上述k-mer。滑動(dòng)窗口沿RNA序列以步長(zhǎng)為1核苷酸進(jìn)行滑動(dòng)匹配，使用ci表示匹配到的次數(shù)()，并且為每個(gè)k-mer分配一個(gè)系數(shù)wk，從而使得每類(lèi)頻率對(duì)預(yù)測(cè)效果的影響一樣，具體如下：

式中：sk是總的匹配次數(shù)，L為RNA序列長(zhǎng)度，fi為每種k-mer出現(xiàn)的頻率，則得到二級(jí)結(jié)構(gòu)特征向量：

本文選擇融合上述3類(lèi)特征組成含90維的特征集作為最終的特征向量：

1.5 構(gòu)建分類(lèi)模型

樸素貝葉斯方法(Naive Bayes)是基于貝葉斯定理的監(jiān)督學(xué)習(xí)算法，即簡(jiǎn)單地假設(shè)每對(duì)特征之間相互獨(dú)立。相比于其他更復(fù)雜的方法，樸素貝葉斯學(xué)習(xí)器和分類(lèi)器非常快，并且有助于解決高維數(shù)據(jù)問(wèn)題。支持向量機(jī)是一種基于統(tǒng)計(jì)學(xué)習(xí)的分類(lèi)方法[18]，其模型參數(shù)確定會(huì)對(duì)應(yīng)到一個(gè)凸最優(yōu)化問(wèn)題，因此可以保證得到最優(yōu)解。目前流行的CPC、CNCI等軟件都使用SVM作為分類(lèi)器。梯度提升決策樹(shù)(gradient boost decision tree，GBDT)是對(duì)于任意可微損失函數(shù)的提升算法的泛化，它具有強(qiáng)大的預(yù)測(cè)能力以及在輸出空間中對(duì)異常點(diǎn)的魯棒性。

結(jié)合訓(xùn)練集高維度、非連續(xù)等特征，以及模型自身的魯棒性，本文選擇基于高斯分布的樸素貝葉斯模型、支持向量機(jī)以及梯度提升決策樹(shù)3個(gè)模型進(jìn)行訓(xùn)練。然后使用網(wǎng)格搜索法分別調(diào)整3個(gè)分類(lèi)器的超參數(shù)。并且采用加權(quán)投票分法來(lái)融合上述3個(gè)分類(lèi)模型的輸出得到最終的預(yù)測(cè)結(jié)果。

1.6 性能評(píng)價(jià)指標(biāo)

本文選擇使用準(zhǔn)確率(Accuracy，ACC)、精確率 (Precison，P)、召回率 (Recall，R)、F1值 (F1_score)來(lái)評(píng)估訓(xùn)練出的分類(lèi)模型。定義如下：

式中：TP指將正類(lèi)預(yù)測(cè)為正類(lèi)數(shù)，F(xiàn)N指將正類(lèi)預(yù)測(cè)為負(fù)類(lèi)數(shù)，F(xiàn)P指將負(fù)類(lèi)預(yù)測(cè)為正類(lèi)數(shù)，TN指將負(fù)類(lèi)預(yù)測(cè)為負(fù)類(lèi)數(shù)。

2 LncRNA功能預(yù)測(cè)

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

用于構(gòu)建互作網(wǎng)絡(luò)的microRNA序列是從miRBase[19](http://www.mirbase.org/index.shtml)下載的427條成熟擬南芥microRNA序列。lncRNA以及mRNA則選用上述下載的具有高可信度的2 565條lncRNA與2 500條mRNA。

2.2 靶向預(yù)測(cè)提取互作對(duì)

研究證明，脅迫作用下，植物的性狀將發(fā)生改變，而這個(gè)過(guò)程是由多個(gè)基因相互作用形成的，lncRNA也參加其中。作為競(jìng)爭(zhēng)性?xún)?nèi)源RNA或者microRNA內(nèi)源性模擬靶標(biāo)的lncRNA，可以有效抑制microRNA的功能，從而間接作用mRNA影響生物形狀及蛋白表達(dá)[20]。首先使用RNAhybrid[21]預(yù)測(cè)lncRNA-microRNA相互作用對(duì)。本文綜合考慮microRNA靶標(biāo)結(jié)合的特征，設(shè)置RNAhybrid參數(shù)：最小自由能?25 kcal/mol，保證種子區(qū)域2～8位完全配對(duì)，p-value值小于0.05。然后，在雜交區(qū)內(nèi)結(jié)合RNA結(jié)構(gòu)信息[22]進(jìn)行篩選：

1) microRNA序列5’端開(kāi)始的10～12位必須有突起點(diǎn)；

2) 中間位置突起點(diǎn)只允許包含lncRNA序列2～4個(gè)核苷酸；

3) 除了中間位置的突起外，microRNA上的錯(cuò)配和G:U配對(duì)總數(shù)少于4并且連續(xù)錯(cuò)配小于2。

利用psRNAtarget靶向預(yù)測(cè)工具預(yù)測(cè)出microRNA-mRNA相互作用對(duì)，并且篩選出有l(wèi)ncRNA靶點(diǎn)的microRNA靶向mRNA的數(shù)據(jù)。

2.3 構(gòu)建調(diào)控網(wǎng)絡(luò)與功能預(yù)測(cè)

融合兩類(lèi)相互作用對(duì)，基于Cytoscape[23]工具包構(gòu)建初級(jí)的lncRNA-microRNA-mRNA互作網(wǎng)絡(luò)，然后對(duì)該網(wǎng)絡(luò)進(jìn)行模塊分析，利用GO[24]數(shù)據(jù)庫(kù)中的術(shù)語(yǔ)了解基因特性。這些GO術(shù)語(yǔ)被劃分為3類(lèi)：細(xì)胞成分(cellular component)，分子功能(molecular function)和生物過(guò)程(biological process)，因此可以基于GO術(shù)語(yǔ)對(duì)各個(gè)模塊進(jìn)行注釋預(yù)測(cè)lncRNA的功能。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 標(biāo)注數(shù)據(jù)測(cè)試集驗(yàn)證結(jié)果

交叉驗(yàn)證(cross validation，CV)是一種模型驗(yàn)證技術(shù)，把給定的數(shù)據(jù)進(jìn)行切分，將切分的數(shù)據(jù)集組合為訓(xùn)練集與測(cè)試集，用于驗(yàn)證模型的泛化能力，有效降低模型的過(guò)擬合問(wèn)題。交叉驗(yàn)證方法可分為簡(jiǎn)單交叉驗(yàn)證法、K折交叉驗(yàn)證以及留一交叉驗(yàn)證。其中，應(yīng)用最多的是K折交叉驗(yàn)證。標(biāo)注數(shù)據(jù)即從公共數(shù)據(jù)庫(kù)采集到的被驗(yàn)證與標(biāo)注的序列數(shù)據(jù)。為了減少計(jì)算消耗的時(shí)間并評(píng)估分類(lèi)模型的泛化性能，本文直接對(duì)分類(lèi)模型在數(shù)據(jù)集上進(jìn)行5折交叉驗(yàn)證，得到的性能效果作為標(biāo)注數(shù)據(jù)測(cè)試集的測(cè)試結(jié)果。

為了驗(yàn)證本文提出方法的有效性，使用目前比較流行的CPAT、CNCI與PLEK軟件在本文采用的數(shù)據(jù)集上進(jìn)行分類(lèi)預(yù)測(cè)，將得到的結(jié)果進(jìn)行比較。CPAT使用邏輯回歸模型；CNCI通過(guò)分析序列的內(nèi)在組成來(lái)區(qū)分蛋白編碼和非蛋白轉(zhuǎn)錄本，使用ATN分?jǐn)?shù)矩陣以及序列結(jié)構(gòu)兩類(lèi)特征；PLEK使用k-mer和滑動(dòng)窗口來(lái)分析轉(zhuǎn)錄本，選取k-mers頻率作為其特征。后兩款軟件都使用支持向量機(jī)作為其分類(lèi)器，結(jié)果如表2所示。可以看出，本文提出的基于多特征融合的集成方法在精確率上超過(guò)90%，優(yōu)于CPAT、CNCI與PLEK；召回率分別比CPAT、CNCI、PLEK高出6.8%、7.4%和8.8%；F1得分也優(yōu)于另外三者。這些結(jié)果表明本文提出的方法可以有效地預(yù)測(cè)植物lncRNA。

表 2 基于不同方法的分類(lèi)結(jié)果比較Table 2 Classification results comparison based on different methods

為了進(jìn)一步驗(yàn)證本文基于多特征融合所構(gòu)建的分類(lèi)模型的有效性，分別給出單獨(dú)使用開(kāi)放閱讀框、二級(jí)結(jié)構(gòu)、k-mers作為特征訓(xùn)練分類(lèi)器得出的預(yù)測(cè)結(jié)果，交叉驗(yàn)證的結(jié)果如表3所示。

表 3 基于不同特征的分類(lèi)結(jié)果比較Table 3 Classification results comparison based on different features

可以看出，本文提出的方法整體的準(zhǔn)確率為89.0%，比單獨(dú)使用開(kāi)放閱讀框、二級(jí)結(jié)構(gòu)、kmers要分別高出6.2%、17%、7.6%，這表明提出的方法對(duì)于識(shí)別lncRNA相較于使用單一類(lèi)特征是有效的。并且可以看出，使用ORF得到的預(yù)測(cè)結(jié)果要優(yōu)于其他兩類(lèi)，這意味著ORF在識(shí)別lncRNA上具有更好的區(qū)分度。

3.2 網(wǎng)絡(luò)構(gòu)建與功能預(yù)測(cè)

經(jīng)過(guò)兩個(gè)靶向預(yù)測(cè)軟件包的預(yù)測(cè)并且對(duì)預(yù)測(cè)結(jié)果按上述規(guī)則進(jìn)行篩選后得到數(shù)據(jù)如表4～表5所示。

表 4 篩選后的microRNA-lncRNA靶點(diǎn)數(shù)據(jù)Table 4 Filtered microRNA-lncRNA target data

表 5 篩選后的microRNA-mRNA靶點(diǎn)數(shù)據(jù)Table 5 Filtered microRNA-mRNA target data

融合以上兩類(lèi)數(shù)據(jù)構(gòu)建的初級(jí)調(diào)控網(wǎng)絡(luò)如圖1所示。

調(diào)控網(wǎng)絡(luò)中的每個(gè)模塊以microRNA為中心，形成microRNA同時(shí)與lncRNA、mRNA相互作用的調(diào)控子網(wǎng)絡(luò)。調(diào)控子網(wǎng)絡(luò)根據(jù)RNA作用數(shù)目和類(lèi)型的不同可以分為：1)單microRNA作用網(wǎng)絡(luò)，即單個(gè)microRNA作為結(jié)點(diǎn)與lncRNA、mRNA相互作用，但與網(wǎng)絡(luò)中其他microRNA結(jié)點(diǎn)沒(méi)有聯(lián)系，如圖2所示；2)多microRNA相互作用網(wǎng)絡(luò)，不同的microRNA通過(guò)靶向同一個(gè)mRNA、lncRNA形成相互作用的模塊，如圖3。

圖1 擬南芥初級(jí)調(diào)控網(wǎng)絡(luò)(三角形代表microRNA，矩形代表lncRNA，圓形代表mRNA)Fig. 1 Primary regulatory network of Arabidopsis( Triangles represent microRNAs, rectangles represent lncRNAs, and circles represent mRNAs)

圖2 單microRNA作用網(wǎng)絡(luò)Fig. 2 Single microRNA interaction network

圖3 多microRNA作用網(wǎng)絡(luò)Fig. 3 Multiple microRNA interaction network

在構(gòu)建調(diào)控網(wǎng)絡(luò)并進(jìn)行模塊分析后，使用GO術(shù)語(yǔ)檢查模塊中的mRNA的功能注釋?zhuān)?duì)和mRNA相關(guān)的lncRNA可能參與的生物調(diào)控過(guò)程進(jìn)行預(yù)測(cè)，部分結(jié)果如表6。可以看到根據(jù)相關(guān)聯(lián)的RNA，本文預(yù)測(cè)的lncRNA所具有的生物調(diào)控功能。例如NONATHT002539參與到氮化合物代謝、分解代謝以及生物合成過(guò)程；NONATHT 000372促進(jìn)蛋白質(zhì)磷酸化；NONATHT002765和NONATHT002470、NONATHT002469都會(huì)影響細(xì)胞轉(zhuǎn)化的過(guò)程等。

表 6 lncRNA功能預(yù)測(cè)Table 6 lncRNA function prediction

4 結(jié)束語(yǔ)

本文基于植物RNA序列，提取開(kāi)放閱讀框、二級(jí)結(jié)構(gòu)和k-mers 3類(lèi)特征，并將它們?nèi)诤铣梢粋€(gè)90維的特征向量作為輸入，訓(xùn)練樸素貝葉斯、支持向量機(jī)、梯度提升決策樹(shù)3種機(jī)器學(xué)習(xí)模型，并采用加權(quán)投票分法來(lái)集成分類(lèi)結(jié)果。通過(guò)與現(xiàn)有的識(shí)別軟件CNCI和PLEK相比，本文提出方法取得了較好的性能，可以有效地識(shí)別預(yù)測(cè)植物lncRNA。基于內(nèi)源性競(jìng)爭(zhēng)規(guī)則，篩選lncRNA-microRNA、microRNA-mRNA相互作用數(shù)據(jù)，并整合兩類(lèi)數(shù)據(jù)構(gòu)建調(diào)控網(wǎng)絡(luò)，基于互作網(wǎng)絡(luò)利用GO術(shù)語(yǔ)對(duì)各個(gè)模塊的mRNA注釋?zhuān)M(jìn)而通過(guò)mRNA預(yù)測(cè)lncRNA功能。未來(lái)將結(jié)合深度學(xué)習(xí)技術(shù)，進(jìn)一步改善預(yù)測(cè)的準(zhǔn)確率。