基于機器學習的RNA編輯位點預測方法綜述

冷嘉承，吳凌云*

(1.中國科學院數(shù)學與系統(tǒng)科學研究院 應用數(shù)學研究所，管理、決策與信息系統(tǒng)重點實驗室，國家數(shù)學與交叉科學中心，北京 100190; 2.中國科學院大學 數(shù)學科學學院,北京100049)

RNA合成、加工、行使功能和降解是細胞生存的關鍵，并在許多不同的層面進行著調(diào)控[1]。RNA合成是基因表達的第一步，轉錄因子調(diào)控RNA聚合酶II(Pol II)與啟動子結合[2]，通過一套非常復雜的操作步驟將DNA轉錄為前體RNA[3]。前體RNA隨后被加工產(chǎn)生成熟mRNA、功能性tRNA和rRNA[4]。RNA的加工包括(1)加帽：將7-甲基鳥苷酸(m7G)添加到5’末端[5]；(2)聚腺苷酸化：在3’末端添加poly-A尾巴[6]；(3)剪接：去除內(nèi)含子之后拼接外顯子[7]；(4)RNA編輯：修改RNA分子序列并導致蛋白質(zhì)多樣性[4, 8]。而RNA編輯作為RNA加工中的一環(huán)，起著至關重要的作用。

RNA編輯有很多種改變RNA序列的機制，其中涉及了堿基的插入和缺失以及堿基替換，如胞苷(C)到尿苷(U)和腺苷(A)到肌苷(I)的脫氨基化，以及非模板化的核苷酸添加和插入[3]。到目前為止，人們已經(jīng)在真核生物、病毒、古細菌以及原核生物中發(fā)現(xiàn)了RNA編輯。所以，生物學家一直對RNA編輯保持著廣泛的關注和強烈的興趣[1-3, 8]。

在哺乳動物中已經(jīng)發(fā)現(xiàn)了兩種類型的RNA編輯。一種是由催化性類多肽載脂蛋白B mRNA編輯酶(APOBEC)催化的C至U編輯，其頻率相當?shù)?。另一種是A-to-I編輯，其中腺苷通過腺苷脫氨酶(ADAR)的作用脫氨基成肌苷，其頻率非常高。因為人類RNA編輯中絕大多數(shù)的編輯都屬于A-to-I編輯，故A-to-I編輯在生物細胞分子機制中尤為重要。研究表明，A-to-I RNA編輯同時也與腦功能、病毒感染和人類疾病有關。例如在人體中，RNA合成和加工中的錯誤可能引起神經(jīng)系統(tǒng)疾病，如三核苷酸擴張性疾病強直性肌營養(yǎng)不良和脆性X綜合征[8]。這些均顯現(xiàn)了RNA編輯研究的重要性與意義。

深度測序技術和生物信息學的發(fā)展使得人們可以在全局篩選A-to-I RNA編輯位點。但是目前為止，通過高通量測序數(shù)據(jù)準確識別RNA編輯事件仍然是一個巨大的挑戰(zhàn)?，F(xiàn)在的方法通常將短讀段映射到參考基因組或轉錄組，然后去除相同的讀段，過濾低質(zhì)量讀段,調(diào)用差異信息并且去除已知的單核苷酸多態(tài)性(SNP)[6, 7, 10, 15-17]。但是將大量的短讀段映射到參考基因組是非常耗時的，并且只有少數(shù)流水線和計算工具可公開用于處理RNA編輯[18]。最關鍵的是，由于成本問題，DNA測序與RNA測序一般不會一起進行，所以很難區(qū)分新的SNP位點和RNA編輯位點。實際上，大量的RNA編輯位點已被注釋為dbSNP中的SNP[19]。

針對這些現(xiàn)象，尤其是為了區(qū)分SNP位點和RNA編輯位點，陸續(xù)產(chǎn)生了許多能夠較為準確預測RNA編輯位點的方法，本文對一些常見方法進行了總結(見表1)，將在第三節(jié)對這些方法進行詳細的介紹。

表1 RNA編輯位點預測方法概覽Table 1 Overview of RNA editing site prediction methods

注：正樣本指發(fā)生RNA編輯的位點(包括所有類型的RNA編輯，以A-to-I編輯為主)，負樣本指發(fā)生堿基變化但沒發(fā)生RNA編輯的位點，其中SVM方法只針對A-to-I編輯。

1 RNA編輯數(shù)據(jù)庫

已有文獻在研究和評價RNA編輯位點預測方法時使用了不同的數(shù)據(jù)集，目前還沒有一個通用的、被廣泛接受的基準數(shù)據(jù)集。除了一些文獻中提供的獨立數(shù)據(jù)集，隨著越來越多的RNA編輯位點被發(fā)現(xiàn)，目前已經(jīng)出現(xiàn)了多個關于RNA編輯的數(shù)據(jù)庫。這些數(shù)據(jù)庫為研究基于計算的RNA編輯位點預測方法提供了很好的訓練數(shù)據(jù)和研究的基礎。

1.1 RADAR數(shù)據(jù)庫

Ramaswami和Li等人于2013年創(chuàng)建了RADAR數(shù)據(jù)庫[20]，包括了人類、蒼蠅、老鼠這三個物種(hg19/mm9/dm3)的A-to-I RNA編輯數(shù)據(jù)。該數(shù)據(jù)庫可以根據(jù)用戶的要求進行篩選，例如基因名稱、是否位于Alu區(qū)域、與其他物種之間的編輯保守性等，并且提供了一些常用數(shù)據(jù)的直接下載，如hg19中的全部RNA編輯位點等。在結果中點擊位置信息可以鏈接到UCSC基因組瀏覽器上瀏覽更詳細的信息，并且都附有數(shù)據(jù)的文獻來源。

1.2 DARNED數(shù)據(jù)庫

DARNED數(shù)據(jù)庫[21]是由Kiran和Baranov創(chuàng)建的。它包括多種數(shù)據(jù)來源：(1)生物信息學分析cDNA序列和基因組序列之間的差異;(2)SNP的分析數(shù)據(jù);(3)miRNA的分析數(shù)據(jù);(4)來自同一組織的RNA和DNA樣品的高通量測序結果[21]。最后將RNA編輯事件的位點映射到參考人類基因組。此數(shù)據(jù)庫不僅包括RADAR中的三個基因組還包括了hg18以及mm10。目前為止該數(shù)據(jù)庫支持三種RNA編輯位點的查詢方式：根據(jù)區(qū)域(染色體、位置、組織特異性等信息)查詢；根據(jù)基因名稱(如APOB)查詢；根據(jù)序列本身查詢。在此數(shù)據(jù)庫的查詢結果中，通過點擊位置信息不僅可以鏈接到UCSC上，還可以鏈接到ENSEMBL上，比較方便。對是否與SNP混淆，與哪個SNP混淆，基因名稱，基因區(qū)域都有詳細的標注，而且都可以通過點擊基因信息鏈接到NCBI上進行查詢，十分人性化。

2 RNA編輯位點預測方法

2.1 基于DNA序列數(shù)據(jù)的有監(jiān)督學習方法

Sun等人[9]提出了一個支持向量機(SVM)模型，基于DNA測序數(shù)據(jù)對A-to-I類型的RNA編輯位點進行預測。其核心思想是，RNA編輯位點主要是由該位點附近的序列決定的。該方法將序列與序列之間的關系進行相關性分析，得到相似矩陣，然后通過映射將其轉化到內(nèi)核空間得到核矩陣，最終利用SVM模型進行訓練，得到具有判別能力的RNA編輯位點分類器(見圖1)。

序列間的相關性是由字符串距離刻畫的：

(1)

其中a,b分別代表輸入的兩個字符串，DEdit代表的是編輯距離，即a,b之間互相轉換最少需要插入、刪除或者替換多少個字符，而DHamming代表的是漢明距離，與編輯距離不同，只允許替換。而L1,L2分別代表計算編輯距離和漢明距離所用的字符串長度(見圖1)。w是0到1之間取值的權重。D的值越大，代表兩個序列之間的相似性越低。該方法在LIBSVM[22]中使用字符串核函數(shù)將字符串數(shù)據(jù)轉換為向量空間。字符串核函數(shù)是對字符串類型數(shù)據(jù)進行操作的核函數(shù)，可寫為：

K(a,b)=exp(-gamma×D(a,b)2)

(2)

其中D是從等式(1)導出的組合距離。伽瑪參數(shù)定義單個訓練實例影響到達的距離,低值表示“遠”，高值表示“近”。

圖1 SVM方法流程圖[9]Fig.1 Flow chart of SVM method[9]

該論文還嘗試了單類SVM。與傳統(tǒng)的支持向量機相比，單類支持向量機嘗試學習一個決策邊界，實現(xiàn)樣本與原點之間的最大分離。他們根據(jù)Sch?lkopf等人的研究結果[23]，使用二進制值1作為編輯事件，-1作為非編輯事件。該文采用5折交叉驗證的方式得到，僅用正樣本作為訓練數(shù)據(jù)的單類SVM的性能要遠低于雙類SVM，單類SVM在參數(shù)nu=0.5時的精確度在果蠅、小鼠、人類數(shù)據(jù)集上分別是0.489, 0.495, 0.498，而雙類SVM模型卻分別達到了0.75, 0.85, 0.79左右。并且單類SVM方法通常無法權衡在正負樣本上性能的差異。該文還用基于人類樣本訓練的模型在一個Sanger測序集上進行了驗證試驗，其中79.3%(46/58)的位點被成功預測。

2.2 基于RNA-seq數(shù)據(jù)的無監(jiān)督學習方法

僅僅基于DNA序列的預測方法需要可靠的RNA編輯位點數(shù)據(jù)作為訓練樣本，而且這類方法對于訓練樣本的依賴性非常高?？紤]到不同物種和不同類型樣本中發(fā)生RNA編輯概率的差異，以及準確的RNA編輯位點數(shù)據(jù)的缺乏，目前更為常見的是從RNA-seq數(shù)據(jù)中識別出RNA編輯位點。按照使用的機器學習方法的不同，基于RNA-seq的RNA編輯位點識別方法又可以大致分為兩類：基于無監(jiān)督學習的過濾方法，和基于有監(jiān)督學習的機器學習方法。第一類方法通過比較RNA-seq數(shù)據(jù)和參考基因組的差異，獲得了潛在的RNA編輯位點，然后通過與SNP數(shù)據(jù)庫的比較以及其他指標進行過濾，去除假陽性位點，最終獲得較為可靠的RNA編輯位點。而第二種方法則通過RNA-seq數(shù)據(jù)和參考基因組提取出一些特征，建立機器學習模型，利用已知的RNA編輯位點數(shù)據(jù)進行有監(jiān)督訓練，獲得RNA編輯位點的預測模型。

Ramaswami等人[10]提出了整合多樣本RNA-seq數(shù)據(jù)識別RNA編輯位點的方法。為了利用大量可公開獲得的RNA-seq數(shù)據(jù)集來發(fā)現(xiàn)RNA編輯位點，該文章提出了兩種相關且互補的方法，以使用來自單個物種中的多個個體的RNA-seq數(shù)據(jù)準確鑒定RNA編輯位點。在第一種方法(見圖2(a))中，在每個RNA-seq樣品中分別將測序讀段映射到(非匹配的)基因組參考序列后，找出RNA的變化，并且將已知常見基因組SNP去除，以此將RNA編輯位點與其余稀有SNP區(qū)分開來。這主要是因為相同的編輯位點通常存在于不同的個體中，而罕見的SNP很可能不存在。

在第二種方法中(見圖2(b))，將不同樣本的RNA-seq的讀段匯總到一起進行比對，從而提高找出RNA變異的靈敏度。然后如第一種方法一樣，按接下來的步驟排除掉SNP，并且找出RNA編輯。由于罕見的SNP不可能出現(xiàn)在多個個體中，所以在匯聚后的比對中將以非常低的頻率存在。

圖2 Multi-Sampled Method的兩種方法[10]Fig.2 Two methods of Multi-Sampled Method[10]

而GIREMI方法[11]則是典型的僅需要單組RNA-seq的方法，GIREMI將RNA-seq讀段中單核苷酸變異(SNV)對的統(tǒng)計推斷模型與機器學習方法結合，以預測RNA編輯位點。GIREMI的輸入包括來自RNA-seq數(shù)據(jù)集的SNV(錯配)列表和公共數(shù)據(jù)庫(如dbSNP)中已知的SNP，輸出是預測的RNA編輯位點及其編輯水平。除了公開的SNP信息外，GIREMI僅僅使用感興趣的RNA-seq數(shù)據(jù)集進行所有分析，而不依賴于任何其他基因組或RNA-seq數(shù)據(jù)集，因此這種方法適用于更廣的范圍(見圖3)。

圖3 GIREMI流程圖[11]Fig.3 Flow chart of GIREMI[11]

GIREMI分為兩個步驟，首先進行SNV位點的互信息計算：

對于每個SNV，我們考慮所有可能的堿基A,T,C,G作為變量si的四種可能狀態(tài)。對于表示一對堿基(si,sj)的聯(lián)合變量，總共有16種可能狀態(tài)。各種狀態(tài)si，sj或(si,sj)的概率可以使用最大似然法進行估計得到?？紤]到所有可能的測序錯誤以及實際數(shù)據(jù)中的低測序深度，假設在實際數(shù)據(jù)中未觀察到的狀態(tài)的概率值為0.01。因此(si,sj)的互信息是:

(3)

其中N={A,T,C,G}，ni和nj分別表示si和sj的狀態(tài)。而SNPsi的信息值被定義為：

(4)

其中S代表帶有si的SNP對的集合，T代表集合S中SNP對的個數(shù)。

這樣，每個RNA-seq數(shù)據(jù)樣本均產(chǎn)生一個基于互信息的SNP信息值分布(I(si)(見圖4)。同樣，對這個樣本的每個SNV位點以同樣的方式(SNV對)求一個信息值，取95%的置信度，如果該位點的信息值落在SNP信息分布的置信區(qū)間外，則判定該位點是RNA編輯位點，否則為SNP位點。

其次，為了提高RNA編輯位點判定的精確度，GIREMI用第一步中識別出的RNA編輯位點作為正樣本來訓練廣義線性模型。該模型采用了兩個特征，一個是讀段中雜合SNP等位基因比率與SNV等位基因比率的差值的絕對值d，另一個則是基于序列本身特征的復合序列得分c(通過位置權重矩陣中+1和-1位置計算)。其回歸模型為：

(5)

其中，β0,β1,β2分別為需要學習的系數(shù)，g為邏輯回歸函數(shù)。

該文只對GM12878數(shù)據(jù)集進行了性能測試，并且達到了99.4%的準確度，但是其準確度的定義為100%-被預測為RNA編輯位點中SNP位點的百分比。

2.3 基于RNA-seq數(shù)據(jù)的有監(jiān)督學習方法

不同于GIREMI基于自主產(chǎn)生的正樣本來訓練廣義線性模型，RED-ML和RDDpred是基于已知樣本進行有監(jiān)督學習訓練的RNA編輯位點檢測模型(見圖5)。

圖4 SNP與RNA編輯位點的信息分布[11]Fig.4 Information distribution of SNP and RNA editing sites[11]

圖5 RED-ML流程圖[12]
Fig.5 Flow chart of RED-ML[12]

Xiong等人[12]建立了一個基于機器學習的RNA編輯位點預測工具RED-ML，并選擇了邏輯回歸(Logistic regression)進行模型的訓練。RED-ML使用的特征有三大類。第一類是基本讀段特征，包括候選位點的支持讀段數(shù)量和計算出的編輯頻率。第二類特征與可能的測序失誤和錯位有關，包括支持讀段的圖譜質(zhì)量、候選位點在定位讀段中的相對位置、鏈偏差的指示、候選位點是否落入簡單重復區(qū)域等。第三類是基于RNA編輯的已知屬性，如編輯類型(是否為A-to-I)，候選位點是否在Alu區(qū)域以及它的序列上下文。需要注意的是，與前兩類特征不同，第三類特征不能直接用來過濾非RNA編輯位點。然而，它們?nèi)匀豢梢蕴峁┯袃r值的信息，通過機器學習方法，將不同來源的信息結合起來做出分類決策。

作為一款RNA編輯檢測軟件工具，RED-ML的輸入是一個BAM文件，也可以利用相應的基因組差異信息。RED-ML將提取候選RNA編輯位點及其相應的特征，然后應用邏輯回歸分類器以相應的置信度檢測真正的RNA編輯位點。RED-ML可以僅基于人類RNA-seq數(shù)據(jù)執(zhí)行全基因組RNA編輯檢測，也可以利用匹配的DNA-Seq數(shù)據(jù)，并與其他常見的RNA-seq數(shù)據(jù)分析步驟很好地結合。

該文用從其他文章中已發(fā)表的數(shù)據(jù)自己篩選出正負樣本進行訓練，其中隨機選取了80%作為訓練集，20%作為測試集。在ROC曲線上的AUC達到了0.98，在PR曲線上的AUC達到了0.94。并且該文在CH24T、CH62T和HeLa樣本上做了RNA-seq驗證實驗，取閾值為0.5時成功驗證了90%的RNA編輯位點。

盡管RNA-seq數(shù)據(jù)可以用于RNA編輯位點檢測，但目前用RNA-seq進行RNA編輯位點檢測的算法也具有相當大的假陽性(False positive)風險，這是用RNA-seq檢測RNA編輯位點的最大挑戰(zhàn)之一。由于短讀段誤對齊而產(chǎn)生的假陽性，本質(zhì)上是由以下幾種因素導致：(1)基因組序列固有的重復片段；(2)模糊的剪切連接；(3)個體之間普遍的多態(tài)性；(4)測序讀段的短缺。RDDpred(RNA/DNA Differences prediction)[14]是一種基于隨機森林算法的RNA編輯位點預測方法，能夠大大減少樣本中的假陽性數(shù)據(jù)，從而提升RNA編輯位點的預測準確率。

RDDpred首先對輸入的測序數(shù)據(jù)進行初始比對，產(chǎn)生編輯位點候選者，然后從中選擇滿足特定條件的樣本作為訓練數(shù)據(jù)(見圖6)。該方法使用RNA編輯數(shù)據(jù)庫RADAR和DARNED中的RNA編輯位點作為正樣本。而負樣本則通過MES(mapping error set)方法來收集，這種方法可以在比對時計算基因組內(nèi)的容易導致錯誤的區(qū)域[13]。從RDD候選者中收集正/負樣本后，所有剩余的樣本被視為預測目標。然后RDDpred建立了一個包含15個特征的隨機森林預測器來預測RNA編輯位點。

圖6 RDDpred 整體流程圖[14]Fig.6 Overall flow chart of RDDpred[14]

RDDpred用來自Bahn和Peng的小組進行的獨立研究的兩個數(shù)據(jù)集進行了評估[13, 15]。在Bahn的研究中，RNA-seq產(chǎn)生了115 132 348個讀段，RDDpred檢測到6 856 440個初始RDD并預測了105 564個RNA編輯位點。在Peng的研究，RNA序列產(chǎn)生了583 640 030個讀段，RDDpred檢測到58 666 976個初始RDD并預測3 076 908個RNA編輯位點。雖然這兩項研究都使用人體組織，但它們產(chǎn)生了不同數(shù)量的RNA編輯位點(105 564與3 076 908)，這表明RNA編輯事件的表達模式在兩種環(huán)境下可能不同。而同時，該方法用自己的模型驗證了Bahn和Peng方法發(fā)現(xiàn)的編輯位點，分別成功驗證了95.32%(3 947/4 141)和90.37%(20 504/22 688)，并且都大幅度減少了錯誤的編輯位點，NPV分別達到了84.21%和75.86%。

3 討 論

本文介紹了五種RNA編輯位點預測方法。第一種方法僅僅需要DNA序列。基于DNA序列的方法主要利用了序列間的相似性進行預測，然而其缺點就是DNA序列中包含的信息是有限的，這使得其性能沒有達到較高的水準(精確度大概在0.7左右)。此外，這種方法無法用于研究RNA編輯在不同條件(例如疾病)、不同個體、不同組織中的差異。

后四類方法則基于RNA-seq測序數(shù)據(jù)和機器學習模型，它們的共同特點就是高度依賴于高通量測序的質(zhì)量與深度。與此同時，還有另外一個因素也限制了這些方法的效果，那就是RNA編輯水平，如果該位點處于一個比較低的編輯水平，那么預測難度將會大大提升。

RED-ML方法提出了大量可能與RNA編輯事件有關的特征，通過邏輯回歸模型進行整合。它的一個缺陷就是模型對訓練數(shù)據(jù)的依賴性很高，例如用人類數(shù)據(jù)得到的模型，在其他動物上的預測效果并不理想。而且這種方法對于比對工具有很強的依賴性(目前版本只優(yōu)化了BWA和TopHat2)，這導致用戶在選擇不同的比對軟件時會產(chǎn)生截然不同的結果阻礙研究的進行。

如何準確地區(qū)分SNP和RNA編輯，這是RNA編輯位點識別的一個核心問題。GIREMI方法通過計算不同SNV位點之間的互信息能夠更準確地區(qū)分SNP和RNA編輯。但是其代價就是計算量的增加，如果考慮到時間成本的話，可能會對效率有所影響。并且要求對測序讀段的長度不能太短，否則無法覆蓋兩個感興趣位點。

RDDpred方法通過RNA編輯數(shù)據(jù)庫和MSE對數(shù)據(jù)進行了篩選，然后用隨機森林模型進行訓練。該文章主要提出了組織特異性導致錯誤正樣本的問題，因為在RADAR和DARNED數(shù)據(jù)庫中有97%的編輯位點都是只存在于一個組織中，如果不清楚地將其篩選出來將會導致“預測危機”，因為這將從根本上(樣本上)導致訓練失敗，從而降低預測性能。如果將其剔除假陽性的方法引用到其他模型中，或許會產(chǎn)生更好的效果，值得讓人期待。

4 總 結

RNA編輯事件的識別對于理解轉錄后調(diào)控是非常重要的。本文首先介紹了RNA編輯的概念和意義，然后介紹了兩個現(xiàn)有的RNA編輯數(shù)據(jù)庫(DARNED、RADAR)。而隨著機器學習的發(fā)展，其在RNA編輯的相關研究中也起到了重要的作用，故本文對已有的RNA編輯位點預測方法進行了概述與討論，得到以下結論：(1)對于RNA編輯位點，我們更關注其在樣本中的表達；(2)對于RNA編輯來說，仍沒有一套像人類參考基因組一樣較為完備的標準；(3)雖然已經(jīng)有了一些基于機器學習的預測RNA編輯位點的方法，但是并沒揭示RNA編輯的本質(zhì)，即提取到判別RNA編輯位點的本質(zhì)特征。所以RNA編輯領域的研究還有很多亟待解決的問題和現(xiàn)象，希望以后能夠通過更深層次的模型去解釋RNA編輯，從而促進相關疾病的研究以及精準醫(yī)療的發(fā)展。