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

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

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

1 材料和方法

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

正集數(shù)據(jù)來(lái)源于數(shù)據(jù)庫(kù)NONCODEv5(http://www.noncode.org/)中擬南芥的3 763條lncRNA，負(fù)集數(shù)據(jù)是從NCBI(https://www.ncbi.nlm.nih.gov/)中隨機(jī)下載的3 800條多種植物的mRNA序列?？紤]到序列的相似性對(duì)預(yù)測(cè)結(jié)果會(huì)造成一定的影響，通過(guò)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)中的一種。

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

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

(2)

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

1.2.2 開放閱讀框

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

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

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

1.2.4 RNA的幾何柔性信息

PseKNC在核苷酸序列的分類上常被作為一種特征[25]，它用離散的模型或者向量表示核苷酸序列，而且通過(guò)其組成寡核苷酸的物理化學(xué)特性保留相當(dāng)多的序列順序信息，特別是全局或者局部的序列順序信息[26]。提取了核苷酸序列的遠(yuǎn)距離和近距離的信息后，能夠更加有助于對(duì)核苷酸序列進(jìn)行分類預(yù)測(cè)。近幾年來(lái)，PseKNC被廣泛用于計(jì)算遺傳學(xué)和基因組學(xué)等各個(gè)領(lǐng)域，像是預(yù)測(cè)DNA甲基化位點(diǎn)，預(yù)測(cè)啟動(dòng)子，預(yù)測(cè)基因組中核小體位置，鑒定microRNA前體等等[27]。考慮RNA局部結(jié)構(gòu)特性，一般相鄰兩個(gè)堿基對(duì)的空間排列有六個(gè)參數(shù)，三個(gè)角度旋轉(zhuǎn)參數(shù)(Tilt，Roll，Twist)和三個(gè)局部平移參數(shù)(Rise，Slide，Shift)[28]，六個(gè)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 支持向量機(jī)算法(Support Vector Machine,SVM)

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

1.4 預(yù)測(cè)性能評(píng)估

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

(5)

(6)

(7)

MCC=

100%

(8)

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

2 結(jié)果分析

2.1 植物lncRNA的序列特征

分析了植物lncRNA序列特征，首先是計(jì)算2 464條植物lncRNA序列所含的單堿基的比例，與2 459條植物mRNA所含單堿基的比例相對(duì)比，見圖1。發(fā)現(xiàn)植物lncRNA富含堿基A和U，而mRNA富含堿基C和G。同樣計(jì)算了2 464條植物lncRNA序列所含的堿基二聯(lián)體的比例，與2 459條植物mRNA所含的堿基二聯(lián)體的比例相對(duì)比，見圖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算法識(shí)別植物lncRNA

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

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

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

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

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

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

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

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

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

3 討 論

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

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

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