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

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

(內蒙古大學 物理科學與技術學院，呼和浩特 010021)

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

1 材料和方法

1.1 數據集

正集數據來源于數據庫NONCODEv5(http://www.noncode.org/)中擬南芥的3 763條lncRNA，負集數據是從NCBI(https://www.ncbi.nlm.nih.gov/)中隨機下載的3 800條多種植物的mRNA序列。考慮到序列的相似性對預測結果會造成一定的影響，通過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的核苷酸序列，當k=1時就代表四種堿基出現的頻數，則RNA序列就能表示成4維的向量，同樣當k=4的時候，則有AAAA，AAAC，...，UUUU不同的組合，RNA序列就能表示成256維的特征向量：

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

(2)

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

1.2.2 開放閱讀框

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

1.2.3 二級結構

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

1.2.4 RNA的幾何柔性信息

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

表1 RNA二核苷酸的柔性結構參數值[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)

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

1.4 預測性能評估

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

(5)

(6)

(7)

MCC=

100%

(8)

其中，TP表示正集序列被預測正確的序列數，TN表示負集序列被預測正確的序列數，FN表示負集序列被預測錯誤的序列數，FP表示正集序列被預測錯誤的序列數。

2 結果分析

2.1 植物lncRNA的序列特征

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

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

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

2.2 SVM算法識別植物lncRNA

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

圖3 不同k-mer組分的預測結果Fig.3 Prediction performance of different k-mer features

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

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

圖4 不同二級結構對預測結果的影響Fig.4 Influence of different secondary structures on prediction results

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

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

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

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

3 討 論

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

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

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