多特征融合的植物長鏈非編碼RNA的預(yù)測

閆玲娟，陳穎麗，閆冬雪，范芷妤

(內(nèi)蒙古大學(xué) 物理科學(xué)與技術(shù)學(xué)院，呼和浩特 010021)

研究表明，編碼蛋白質(zhì)的基因只占整個基因組的一小部分，大部分都屬于非編碼區(qū)域，在非編碼區(qū)域中有大部分會經(jīng)過轉(zhuǎn)錄產(chǎn)生RNA，但是這些RNA不經(jīng)過翻譯的步驟，這類RNA稱為非編碼RNA(Non-coding RNA,ncRNA)[1]。非編碼RNA不僅數(shù)量龐大，種類也有很多[2]。其中長鏈非編碼RNA(Long non-coding RNA,lncRNA)屬于分子長度大于200個堿基的調(diào)控非編碼RNA[3]。最初lncRNA并不被重視，被認(rèn)為是轉(zhuǎn)錄過程中的“噪音”，但是隨著實驗技術(shù)和預(yù)測算法的改進，不僅鑒定和預(yù)測到越來越多的lncRNA，而且還發(fā)現(xiàn)lncRNA在動植物中都具有重要的生物學(xué)功能[4][5]。目前大量的lncRNA的研究還主要集中在人和動物等物種中，相對來說對于植物lncRNA的研究還比較落后[6]。隨著高通量測序技術(shù)和計算方法的發(fā)展，在過去的幾十年中，大量的lncRNA在不同的植物中被發(fā)現(xiàn)，例如擬南芥、水稻、玉米、小麥、黃瓜、番茄等[7]。通過生物信息學(xué)分析手段可以預(yù)測和鑒定越來越多的lncRNA[8]。2017年Mohan Singh等人開發(fā)的預(yù)測植物lncRNA的工具PLncPRO[9]，分別在干旱和鹽脅迫條件下，在水稻和鷹嘴豆中發(fā)現(xiàn)了3 714和3 457個高可信度的lncRNA。2018年常征等[10]也通過機器學(xué)習(xí)方法對植物lncRNA進行了預(yù)測，從PNRD數(shù)據(jù)庫中下載了2 565條lncRNA為正集數(shù)據(jù)，負(fù)集是從RefSeq數(shù)據(jù)庫下載的2 500條mRNA，采用5折交叉驗證的準(zhǔn)確率為89%，用了同樣的數(shù)據(jù)集在CPAT、CNCI、PLEK等軟件上進行分類預(yù)測的結(jié)果分別是85.7%、82.7%、71.4%。可見計算機預(yù)測的方法對植物lncRNA的預(yù)測起到了很好的作用，但是由于lncRNA特殊的序列屬性，使得lncRNA的鑒定工作仍然面臨著挑戰(zhàn)[11]，到目前來說，lncRNA的精確識別仍然是植物研究領(lǐng)域的主要問題之一[9]。本文通過新建植物lncRNA和mRNA數(shù)據(jù)集，提取lncRNA的序列及結(jié)構(gòu)特征，并將多特征融合，利用支持向量機算法對植物lncRNA進行了預(yù)測，取得了較好的預(yù)測效果。

1 材料和方法

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

正集數(shù)據(jù)來源于數(shù)據(jù)庫NONCODEv5(http://www.noncode.org/)中擬南芥的3 763條lncRNA，負(fù)集數(shù)據(jù)是從NCBI(https://www.ncbi.nlm.nih.gov/)中隨機下載的3 800條多種植物的mRNA序列。考慮到序列的相似性對預(yù)測結(jié)果會造成一定的影響，通過CD-Hit軟件[12]去除冗余，取相似性小于60%的序列，最后得到2 464條lncRNA，2 459條mRNA。

1.2 特征提取

1.2.1 lncRNA序列的k-mer核苷酸組分信息

從lncRNA序列中提取核苷酸組分特征[13]，如果將lncRNA序列表示為

S=A1，A2，A3，A4，...，AL-1，AL

(1)

其中Aj就是四種堿基(腺嘌呤A，胞嘧啶C，鳥嘌呤G，尿嘧啶U)中的一種。

對于一個長度為L的核苷酸序列，當(dāng)k=1時就代表四種堿基出現(xiàn)的頻數(shù)，則RNA序列就能表示成4維的向量，同樣當(dāng)k=4的時候，則有AAAA，AAAC，...，UUUU不同的組合，RNA序列就能表示成256維的特征向量：

X=[x1，x2，...，x256]

(2)

考慮到堿基的化學(xué)特性，將四種核苷酸進行了約化，約化分別是嘌呤嘧啶約化和強弱鍵約化，腺嘌呤和鳥嘌呤都是嘌呤，都用M來表示，胞嘧啶和尿嘧啶都是嘧啶，都用N來表示；又因為鳥嘌呤和胞嘧啶之間以3個氫鍵形成堿基配對，分子能量的穩(wěn)定性較高，所以將鳥嘌呤和胞嘧啶都用P來表示，而腺嘌呤與尿嘧啶是以兩個氫鍵形成堿基對，穩(wěn)定性相對較弱，將腺嘌呤和尿嘧啶都用Q表示。對RNA序列進行約化后重新提取k-mer信息，此時k取4的時候RNA序列就可以表示成24=16維的特征向量。

1.2.2 開放閱讀框

開放閱讀框(Open reading frame,ORF)的長度是常被用來區(qū)分lncRNA和mRNA的最基本的標(biāo)準(zhǔn)之一[14]，然而，專門用于預(yù)測ORF的生物信息學(xué)工具很少，sORF finder[15]根據(jù)編碼序列間的核苷酸組成偏見，并通過同義和非同義替換率評估的氨基酸水平上的潛在功能限制進行編碼sORFs的識別，但是它是將近十年前開發(fā)的。還有一些編碼潛能的評估工具，例如CPC、CPAT、CNCI、CPC2、LGC等都可以用來預(yù)測ORF[16]，尤其LGC是在基于ORF長度和GC含量之間的特征關(guān)系來評估編碼潛能的，它在計算分析從植物到哺乳動物等多種物種方面都具有廣泛的應(yīng)用潛力[17]，本文就是利用LGC來識別植物lncRNA的ORF。為了找到關(guān)于開放閱讀框的最佳的特征集，構(gòu)建了最長開放閱讀框的長度以及它的相對長度兩組特征，一個開放閱讀框的相對長度是由其長度除以相應(yīng)的轉(zhuǎn)錄本長度來定義的[18]。

1.2.3 二級結(jié)構(gòu)

RNA序列是由四種不同的堿基組成，RNA二級結(jié)構(gòu)也是由不同的結(jié)構(gòu)元素組成，像是莖區(qū)和環(huán)狀結(jié)構(gòu)等[19]。RNA的結(jié)構(gòu)在很大程度上決定了它的功能，識別RNA分子的二級結(jié)構(gòu)是了解其生物學(xué)功能的有效途徑[20]。結(jié)構(gòu)特征的提取是利用機器學(xué)習(xí)的方法預(yù)測lncRNA的重要步驟之一，但是目前還沒有合適的結(jié)構(gòu)特征提取工具[21]。本文使用了RNAfold軟件[22]預(yù)測了lncRNA的二級結(jié)構(gòu)，得到了二級結(jié)構(gòu)點括號表示形式，括號表示配對的堿基，形成莖結(jié)構(gòu)，點表示沒有配對的堿基，形成單鏈或環(huán)結(jié)構(gòu)，最后的數(shù)字表示釋放后的最小自由能[23]。將lncRNA的二級結(jié)構(gòu)頸環(huán)個數(shù)及歸一化的最小自由能等作為特征，歸一化的最小自由能是由釋放后的最小自由能除以對應(yīng)序列的長度來定義的[24]。

1.2.4 RNA的幾何柔性信息

PseKNC在核苷酸序列的分類上常被作為一種特征[25]，它用離散的模型或者向量表示核苷酸序列，而且通過其組成寡核苷酸的物理化學(xué)特性保留相當(dāng)多的序列順序信息，特別是全局或者局部的序列順序信息[26]。提取了核苷酸序列的遠(yuǎn)距離和近距離的信息后，能夠更加有助于對核苷酸序列進行分類預(yù)測。近幾年來，PseKNC被廣泛用于計算遺傳學(xué)和基因組學(xué)等各個領(lǐng)域，像是預(yù)測DNA甲基化位點，預(yù)測啟動子，預(yù)測基因組中核小體位置，鑒定microRNA前體等等[27]?？紤]RNA局部結(jié)構(gòu)特性，一般相鄰兩個堿基對的空間排列有六個參數(shù)，三個角度旋轉(zhuǎn)參數(shù)(Tilt，Roll，Twist)和三個局部平移參數(shù)(Rise，Slide，Shift)[28]，六個RNA結(jié)構(gòu)信息參數(shù)值詳細(xì)見表1。將RNA序列表示為4k+λ維的特征向量[26]：

表1 RNA二核苷酸的柔性結(jié)構(gòu)參數(shù)值[29]Table 1 Parameter values of flexible structure of RNA dinucleotide[29]

R=[d1d2…d4kd4k+1…d4k+λ]T

(3)

(4)

1.3 支持向量機算法(Support Vector Machine,SVM)

支持向量機的基本思想是尋找兩個類之間的最大邊界超平面，對于非線性數(shù)據(jù)，使用核函數(shù)將它映射到線性的高維空間中，然后在高維空間中擬合一個線性函數(shù)去解決非線性分類問題[30]。文中采用的是Chan和Lin小組開發(fā)的LIBSVM軟件包[31]。

1.4 預(yù)測性能評估

預(yù)測算法的性能是常用的敏感性(Sensitivity,Sn)、特異性(Specificity,Sp)、總體預(yù)測成功率(Acc)以及馬修相關(guān)系數(shù)(Mathew's Correlation Coefficient,MCC)，定義為：

(5)

(6)

(7)

MCC=

100%

(8)

其中，TP表示正集序列被預(yù)測正確的序列數(shù)，TN表示負(fù)集序列被預(yù)測正確的序列數(shù)，F(xiàn)N表示負(fù)集序列被預(yù)測錯誤的序列數(shù)，F(xiàn)P表示正集序列被預(yù)測錯誤的序列數(shù)。

2 結(jié)果分析

2.1 植物lncRNA的序列特征

分析了植物lncRNA序列特征，首先是計算2 464條植物lncRNA序列所含的單堿基的比例，與2 459條植物mRNA所含單堿基的比例相對比，見圖1。發(fā)現(xiàn)植物lncRNA富含堿基A和U，而mRNA富含堿基C和G。同樣計算了2 464條植物lncRNA序列所含的堿基二聯(lián)體的比例，與2 459條植物mRNA所含的堿基二聯(lián)體的比例相對比，見圖2。發(fā)現(xiàn)植物lncRNA的AA/AU/UA/UU二聯(lián)體的頻數(shù)也是明顯比mRNA高。

圖1 兩類序列中的四種核苷酸組分Fig.1 Four kinds of nucleotide components in two sequences

圖2 兩類序列中的二聯(lián)體組分Fig.2 Percentage of diomorphic component in two sequences

2.2 SVM算法識別植物lncRNA

根據(jù)計算分析植物lncRNA的序列特征，發(fā)現(xiàn)植物lncRNA也有一定的序列偏好特征。所以，本文提取lncRNA的k-mer序列特征作為SVM的輸入向量來識別植物lncRNA，基于Jackknife檢驗的不同特征的預(yù)測結(jié)果見圖3。k值的范圍是從1取到6，從圖3中可以看出，在k取4的時候總體預(yù)測成功率是最高的。以約化后的k-mer信息作為特征向量時，取得的總體預(yù)測成功率相對來說并不是很好，但是強弱鍵約化后的總體預(yù)測成功率明顯要比嘌呤嘧啶約化后的總體預(yù)測成功率高，可見強弱鍵約化在識別植物lncRNA時是比嘌呤嘧啶約化更好的一個序列特征。

圖3 不同k-mer組分的預(yù)測結(jié)果Fig.3 Prediction performance of different k-mer features

用基于ORF的兩個特征識別植物lncRNA時，以最長開放閱讀框的長度作為特征向量輸入到SVM算法中，總體預(yù)測成功率達到87.26%，而用其相對長度作為特征向量時，總體預(yù)測成功率達到88.26%，相對長度作為特征向量比直接以其長度作為特征向量時的總體預(yù)測成功率要高，相對長度通常用作長度的補充特征，在分類算法中有更好的表現(xiàn)[32]。

用RNAfold軟件預(yù)測出的lncRNA的二級結(jié)構(gòu)頸環(huán)個數(shù)及歸一化的最小自由能作為特征向量進行預(yù)測時，預(yù)測結(jié)果見圖4，其中單一特征中莖的個數(shù)預(yù)測效果相對來說是比較好的，總體預(yù)測成功率為72.17%，將這些單一特征融合后進行預(yù)測，最好的預(yù)測成功率達到78.36%。

圖4 不同二級結(jié)構(gòu)對預(yù)測結(jié)果的影響Fig.4 Influence of different secondary structures on prediction results

使用Pse-in-one[33]軟件對序列的核苷酸之間的空間位置信息進行提取時，有兩個參數(shù)ω和λ，ω是權(quán)重因子，取值范圍是0到1，為了找到最佳的ω和λ的值，計算了λ的步長為5時對應(yīng)的總體預(yù)測成功率，見圖5。從圖5中可以看出，λ步長為5且取值在1到30之間時，隨著ω的增加，總體預(yù)測成功率逐漸降低，在ω相同時，λ值越大，預(yù)測成功率是偏小的，總體來看，當(dāng)λ=5，ω=0.1時，得到的總體預(yù)測成功率最大為85.9%。

圖5 λ步長為5偽核苷酸特征分類準(zhǔn)確率Fig.5 Classification accuracy of pseudonucleotides with λ Steps 5

綜合這些序列和結(jié)構(gòu)特征信息的預(yù)測結(jié)果，結(jié)果比較好的特征有4-mer組分信息、最長開放閱讀框的長度和相對長度以及PseKNC在參數(shù)λ=5，ω=0.1時的結(jié)構(gòu)信息。將這些特征信息融合后進行預(yù)測，預(yù)測結(jié)果見表2。

表2 Jackknife檢驗下不同特征融合后的SVM預(yù)測結(jié)果Table 2 Prediction performance of SVM model fusing different features under Jackknife test

3 討 論

研究發(fā)現(xiàn)植物lncRNA也有一定的序列和結(jié)構(gòu)偏好特征。在提取k-mer特征信息時，隨著k的增加，特征向量的維數(shù)在以2k增加，由于特征向量的維數(shù)過大的時候會導(dǎo)致序列信息的冗余現(xiàn)象，因此將k的取值只取到6，結(jié)果發(fā)現(xiàn)在k=4的時候預(yù)測結(jié)果最好，在Jackknife檢驗下，總體預(yù)測成功率達到93.36%?？紤]到堿基的化學(xué)特性和植物lncRNA功能相關(guān)，將堿基約化后進行預(yù)測，分析兩種不同的堿基約化方式預(yù)測的結(jié)果發(fā)現(xiàn)，強弱鍵約化后預(yù)測取得的成功率比嘌呤嘧啶約化取得的成功率高，說明強/弱鍵約化(PQ約化)更能反應(yīng)植物lncRNA的序列信息，有利于植物lncRNA的識別。在RNA序列的幾何柔性信息中，用不同的參數(shù)預(yù)測結(jié)果表明，當(dāng)參數(shù)ω=0.1、λ=5的時候預(yù)測結(jié)果最好，總體預(yù)測成功率達到85.9%。

雖然這些特征信息的預(yù)測結(jié)果都比較好，但是單一的特征去預(yù)測總是有一定的局限性，所以為了能夠提取到更多的植物lncRNA序列中蘊藏的結(jié)構(gòu)和功能的信息，將這些特征信息進行了融合，用融合后的特征再去對植物lncRNA進行預(yù)測，發(fā)現(xiàn)結(jié)果有了一定的提高，但是融合信息太多也會造成信息冗余，所以只融合序列和結(jié)構(gòu)信息參數(shù)較優(yōu)的幾個特征，對植物lncRNA預(yù)測的總體成功率達到了96.14%，敏感性達到了96.51%，特異性達到了95.77%，馬修相關(guān)系數(shù)的值是0.92。說明最長開放閱讀框的相對長度以及4-mer組分信息的融合對植物lncRNA的預(yù)測很有效，而且發(fā)現(xiàn)最長開放閱讀框的相對長度和其它一些序列及結(jié)構(gòu)特征信息的組合對預(yù)測結(jié)果都能有一定程度的提高，由此可見，最長開放閱讀框的相對長度在植物lncRNA和mRNA的分類預(yù)測中是一個重要的特征信息。

隨著鑒定和預(yù)測出的lncRNA越來越多，植物lncRNA也開始越來越受到關(guān)注，雖然植物lncRNA的研究相對于動物lncRNA的研究還遠(yuǎn)遠(yuǎn)落后，但是動物lncRNA的研究策略為植物lncRNA的研究提供了很好的借鑒。下一步，為了更加準(zhǔn)確地識別植物lncRNA，可以再深入的挖掘一些更加全面的特征信息，將有效的特征信息融合后再去預(yù)測，或者合理的去融合一些比較優(yōu)越的算法提高預(yù)測結(jié)果。