環狀RNA計算預測方法的研究進展

譚超俊，顧萬君，謝雪英

(東南大學 生物科學與醫學工程學院，南京 210096)

環狀RNA是一類在前體mRNA(pre-mRNA)剪接過程中形成的新型內源RNA，與傳統線性RNA不同的是其通過反向剪接使5’端和3’端共價連接，形成閉合環狀結構(見圖1)，因而不受核酸外切酶介導的降解的影響，比線性RNA更穩定[1-2]。多年來，環狀RNA被認為是前體mRNA剪接過程中的副產品[3]，但越來越多的研究表明環狀RNA在生物體中通過多種途徑發揮著重要的作用。環狀RNA作為競爭性內源RNA，可以調控基因的轉錄和表達[4-5]；一些環狀RNA包含多個microRNA結合位點，可以充當microRNA海綿[6]；環狀RNA可以與蛋白相互作用，從而參與多種生物過程的調控[5,7]；另外，環狀RNA還能夠翻譯蛋白[8]。

圖1 環狀RNA的形成Fig.1 Formation of circRNA

為了更好地研究環狀RNA，準確有效地鑒定環狀RNA是至關重要的。多種因素可以促進環狀RNA的形成，如側翼區域的互補序列[9]、反向重復序列[10]、ALU和串聯重復序列[2]以及SNP密度[11]等。這些因素與RNA分子的進化保守性和結構特征被認為是鑒定環狀RNA的重要特征。目前，常見的環狀RNA識別工具是通過識別高通量測序(RNA-seq)數據中的反向剪接位點來鑒定環狀RNA，如find_circ[12]、CIRI[13]、circRNA_finder[14]、MapSplice[15]和CIRCexplorer[16]等。然而，已有相關研究[17-18]表明這些工具普遍存在較高的假陽性率和假陰性率，不同工具對同一個測序數據的檢出重合率非常低，這是因為首先這些工具是基于高通量測序數據的，因此對表達豐度非常敏感，但大部分環狀RNA通常是低表達的，在測序覆蓋率較低情況下難以捕獲，其次它們只利用了反向剪接位點信息來識別環狀RNA，忽略了環化過程中其他因素的影響。

近年來，機器學習方法越來越多地應用于生物信息學研究。一些研究[19-20]分析環狀RNA形成過程中的影響因素，通過訓練傳統的機器學習算法(支持向量機、隨機森林和多核學習等)來鑒別環狀RNA和長鏈非編碼RNA(Long non-coding RNA，lncRNA)，取得了較高的識別正確率。但是這些方法需要先進行特征分析，而且這些選取的特征不能全面充分地表征反向剪接過程。深度學習算法能夠處理大規模數據并自動提取有效特征，可以彌補傳統機器學習模型的不足。本文將從傳統機器學習方法和深度學習方法這兩個方面來介紹基于序列計算預測環狀RNA的8種工具，并比較分析它們在測試數據集上的識別結果。

1 基于傳統機器學習方法的工具

1.1 PredcircRNA

PredcircRNA[19](https://github.com/xypan1232/predcircrna)采用了基于多種特征訓練的多核學習框架模型，用來區分環狀RNA和其他長鏈非編碼RNA(其運行流程見圖2(a))。首先，從轉錄本中提取不同的特征—圖特征、保守性分數、序列組成、ALU和串聯(Tandem)重復序列、SNP密度和開放閱讀框(Open Reading Frame，ORF)等。圖特征用節點表示核苷酸，用邊表示核苷酸之間的鍵關系，可以表示RNA分子的序列和結構。保守性分數是根據UCSC中下載的phyloP(Phylogenetic p-values)保守分數自定義的。序列組成特征包含三核苷酸特征、GC含量、序列長度、GT、AG、GTAG和AGGT的頻率等。另外計算每個轉錄本的Alu重復序列數目，利用Tandem Repeats Finder檢測串聯重復序列并計算其頻率，txCdsPredict獲得每個轉錄本的開放閱讀框并提取其長度和比例。為了融合多種特征，PredcircRNA使用多核學習方法，能有效地對環狀RNA和其他lncRNA進行分類。PredcircRNA分析的這些特征(如保守性分數和GT/AG序列)對于鑒別環狀RNA和lncRNA有重要作用，而且不同類型的特征相較于單一特征能相互補充從而提高模型能力。PredcircRNA正樣本來自circBase數據庫中的14 084條環狀RNA，負樣本來自GENCODE數據庫中19 722條其他類型的lncRNA，隨機選取10 000條環狀RNA和相同數量的其他lncRNA進行模型訓練，剩余的數據作為獨立的測試數據集，通過5折交叉驗證，準確率達到0.862。

1.2 H-ELM

H-ELM[20]利用有特征選擇功能的H-ELM(hierarchical extreme learning machine，層次極限學習器)算法來提取特征，進一步識別環狀RNA和lncRNA(其流程原理見圖2(b))。該方法沿用了PredcircRNA中定義的特征，使用mRMR(minimum redundancy maximum relevance，最小冗余最大相關)方法對這些特征進行分析。利用獲得的特征列表、IFS(incremental feature selection，增量特征選擇)方法和H-ELM算法，建立最優分類模型。相較于PredcircRNA，H-ELM用mRMR方法對特征進行了分析選擇，利用了IFS方法和H-ELM算法來建立模型。H-ELM使用PredcircRNA的數據，雖然在同樣的數據集上通過十折交叉驗證，H-ELM模型的準確率為0.789，低于PredcircRNA，但它通過特征分析，發現進化保守性、序列特征、特異性序列和結構是區分環狀RNA和其他lncRNA的重要因素。

1.3 predict_circ

predict_circ[21]通過選取剪接位點側翼上下游內含子的長度、A-to-I密度、ALU重復序列和RNA結合蛋白(RBP)作用位點等100個與RNA成環相關的序列特征，建立機器學習模型來識別環狀RNA(其原理見圖2(c))，并比較了隨機森林和支持向量機的分類效果。結果表明，選取的序列特征能有效地鑒別RNA能否成環，同時不同序列特征對模型的分類預測能力的貢獻也不同。predict_circ共選取了3組數據集和1組獨立測試集，正樣本來自5種環狀RNA預測工具的檢測結果交集、circBase中的人類環狀RNA以及收集文獻中的經過PCR驗證的環狀RNA，負樣本來自UCSC人類基因數據hg19版本的編碼蛋白序列(去除與環狀RNA重合的轉錄本)，隨機抽取與每組正樣本集數量相近的序列作為對應的負樣本集。將每組數據的正負數據集中隨機抽取2/3作為訓練集，1/3作測試集，5折交叉驗證后，每組數據的分類準確率均達到0.85以上。

圖2 基于傳統統計學習方法的工具流程圖Fig.2 Flowchart of tools based on statistical learning

1.4 CirRNAPL

CirRNAPL[22](http://server.malab.cn/CirRNAPL/)使用了基于粒子群優化(particle swarm optimization，PSO)算法的極限學習機(extreme learning machine，ELM)模型，能準確識別circRNA(其流程見圖2(d))。首先，提取序列數據的四個特征：核糖核酸組成、序列自相關特征、偽核糖核酸組成和預測結構組成。核糖核酸組成包括k-mer(k參數為2和3)、錯配和子序列；序列自相關特征主要表征序列中核苷酸間的相關性；預測結構組成主要表示序列的結構特征。然后使用ELM識別環狀RNA，并且通過利用PSO算法優化ELM的參數，提高其泛化能力，達到更高的識別精度。CirRNAPL選取了三組數據對模型進行了訓練。第一組和第二組的正樣本均使用PredCircRNA方法中使用的14 084條環狀RNA數據，第一組和第二組的負樣本分別為GENCODE v19版本的9 533條編碼蛋白基因(protein-coding genes，PCGs)和1 973條lncRNAs；第三組數據的正樣本是H1hsec干細胞中表達的2 082條circRNA，反之，負樣本為H1hsec干細胞中未表達的相同數量的circRNA。將三組數據分為了訓練集和測試集，并進行了十折交叉驗證，最終在三組數據上的準確率分別為0.815，0.802和0.782。

2 基于深度學習方法的工具

2.1 DeepCirCode

DeepCirCode[23](https://github.com/BioDataLearning/DeepCirCode)是第一個采用深度學習模型來預測mRNA是否能反向剪接形成環狀RNA的分類工具(其流程見圖3(a))。該方法采用卷積神經網絡自動地從序列中學習相關特征——序列基序(Sequence motif)。因為有研究表明某些RBP等作用因子能通過特定的結合位點(序列基序)來促進RNA環化。DeepCirCode將候選反向剪接位點側翼的內含子和外顯子序列轉換成二進制向量作為網絡輸入，通過識別側翼序列中是否存在能促進環化的序列基序來預測環狀RNA。通過分析DeepCirCode檢測出的序列基序，發現其中一些確實與已知的RNA剪接、轉錄或翻譯的基序相匹配。此外，通過對小鼠和果蠅數據進行測試，發現一些人類序列基序在小鼠和果蠅的序列中也存在，這說明這些基序在進化過程中存在保守性，很可能在環狀RNA的生物發生過程中起著重要作用。DeepCirCode將circBase和circRNADb[24]兩個公共數據庫的環狀RNA作為正樣本數據集，按條件篩選出共7 964條人類外顯子環狀RNA，負樣本是從GENCODE的人類參考基因組注釋信息中隨機選取相應的剪接位點，10折交叉驗證的準確率為0.852 4，AUC為0.905 8。

2.2 circDeep

circDeep[25](https://github.com/UofLBioinformatics/circDeep)采用端到端(End-to-End)的深度學習框架來區別環狀RNA和lncRNA(流程原理見圖3(b))。circDeep引入了三種描述符(descriptor)：RCM(Reverse Complement Matching，反向互補匹配)描述符、ACNN-BLSTM序列描述符和保守性描述符。RCM描述符目的是選取促進環化過程的潛在側翼序列。ACNN-BLSTM序列描述符結合了ACNN(Asymmetric Convolution Neural Network，非對稱卷積神經網絡)和BLSTM(Bidirectional Long Short-Term Memory Network，雙向長短期記憶網絡)，能夠從每個序列中提取局部模式和遠程作用(Long-range dependencies)。保守性描述符包含物種間特殊序列的保守信息和保守基序特征。為了融合三種不同的描述符，它使用了一種從不同方面的信息構建非線性表示的深度學習架構。其正樣本來自circRNADb數據庫的31 939條人類環狀RNA，負樣本來自GENCODE的19 683條其他類型lncRNA。將每個數據集劃分為訓練集(75%)、驗證集(10%)和測試集(15%)，模型訓練后，測試集上的結果準確率達到0.941 7。

2.3 CRC(Contextual Regression for CircRNA)

CRC[26](https://github.com/chl556/Contextual_Regression_for_CircRNA)基于環狀RNA反向剪接位點的側翼區域特征—CpG島(Where long noncoding RNAs meet DNA methylation.)、RBP結合位點、簡單重復序列、A-to-I RNA編輯位點和序列，通過上下文回歸(contextual regression)模型來預測環狀RNA的形成(流程原理見圖3(c))，接著還通過特征提取(feature extraction technique)和PCA獲得10個特征主成分，在此基礎上運用K均值聚類，將環狀RNA分成7種亞型，這些亞型分別對應于已有的環狀RNA生物發生機制。因此作者認為人類環狀RNA具有多種不同的生物發生機制，可以分成多個不同的亞型。此外，CRC還發現環狀RNA生物發生與側翼區域CpG島之間新的關聯以及鑒定了相關的RNA結合蛋白。CRC從circNet數據庫中收集55 689個人類環狀RNA反向剪接位點作為正樣本，在hg19人類基因組上隨機選擇等量的位點作為負樣本，然后將數據隨機劃分為訓練集和測試集(比例為7∶3)，通過十次訓練，達到平均準確率為0.726和AUC值為0.801。

2.4 JEDI

JEDI(Junction Encoder with Deep Interaction)[27](https://github.com/hallogameboy/JEDI)用深度學習方法對剪接位點及其深層相互作用建立模型，直接從基因或轉錄本序列中預測環狀RNA。JEDI對序列里每個外顯子和內含子連接位點進行基于深度雙向循環神經網絡的編碼，然后用交叉注意層(Cross-attention layer)對反向剪接位點的深層相互作用建模(見圖3(d))。JEDI不僅能夠預測環狀RNA，而且能夠解釋剪接位點間的關系，從而發現基因內的反向剪接。另外對小鼠環狀RNA的研究結果表明，JEDI預測人類環狀RNA的模型也適用于小鼠環狀RNA數據。JEDI選取了三組數據。第一組：正樣本來自circRNADb的31 939條人類環狀RNA，負樣本來自GENCODE參考注釋的19 683條其他lncRNA，進行5次交叉驗證后所得準確率達到0.989 9；第二組：正樣本來自circRNADb的每條環狀RNA對應的7 777條基因，對于負樣本，去除了正樣本中選擇的基因，得到7 000條基因，經過5折交叉驗證，模型達到準確率0.964 6；第三組：正樣本來自circBase的1 522條小鼠環狀RNA，負樣本來自GENCODE參考注釋的1 522條其他lncRNA，模型訓練后，準確率為0.886 8。

圖3 基于深度學習方法的工具流程圖Fig.3 Flowchart of tools based on deep learning

3 環狀RNA預測工具的性能比較

為了在同一標準下比較分析以上算法的功能，我們利用公共數據庫circRNADb中人類環狀RNA數據對以上工具進行測試。其中，H-ELM未提供代碼下載鏈接，PredcircRNA工具所需比對數據庫目前不提供下載，CRC需要對數據計算一組特征值，但相應的計算代碼并未提供，CirRNAPL雖然提供了網站服務，但是網站上不方便處理大量的數據，也沒有提供工具的下載，所以本文最后只進行了基于深度學習算法的三種工具(DeepCirCode，circDeep和JEDI)的比較。circRNADb公共數據庫中共有32 194條人類環狀RNA，去除長度短于200 nt的環狀RNA，并過濾掉剪接位點側翼內含子序列和兩端外顯子序列短于50 nt的環狀RNA，最終獲得13 264條序列作為正樣本數據集。提取GENCODE[28]v19版本人類參考基因組注釋的其他類型的lncRNA，剔除與circBase[29]和circRNADb中環狀RNA序列相重疊的序列，得到8 125條lncRNA作為負樣本數據集。分別將正負樣本集中75%的數據作為訓練集，25%的數據作為測試數據，對以上三種模型進行訓練和測試，并采用以下指標來評估模型在測試集上的性能：準確度(Acc)、靈敏度(Sn)、特異性(Sp)和馬修斯相關系數(MCC)，分別定義如下：

(1)

(2)

(3)

(4)

結果見表1和圖4，這表明這三種基于深度學習算法的工具對于識別環狀RNA都有較好的效果，尤其是JEDI，在測試集上的識別正確率達到了97.89%。這三個工具中，circDeep運行時間最長，因為該算法需要耗費大量時間提取特征。

表1 模型分類性能打分Table 1 Performance in model classification %

圖4 測試分類結果ROC曲線Fig.4 ROC curves of prediction results of testing dataset

4 總結和展望

環狀RNA通過非經典方式進行反向剪接而成，通常認為剪接位點側翼區域的反向互補序列和RBP結合位點等是促進內含子區域配對從而介導反向剪接形成環狀RNA分子。本文主要介紹了8種基于序列預測環狀RNA的工具。這8種工具均基于RNA序列來挖掘其內在特征，利用不同的機器學習算法來識別環狀RNA，其優缺點見表2所示。PredcircRNA、H-ELM以及predict_circ發展了不同的策略來提取特征，并使用了傳統的統計學習算法(PredcircRNA的多核學習，H-ELM的層次極端學習機，predict_circ的支持向量機和隨機森林及CirRNAPL的基于粒子群優化算法的極限學習機)來構建分類器。這一類方法需要預先進行選擇和計算特征，而且提取的特征是專門用于描述序列某方面的性質，因此需要一定的先驗知識為基礎。DeepCirCode、circDeep、CRC和JEDI使用深度學習算法可以自動地從原始序列中學習復雜模式。DeepCirCode使用卷積神經網絡對反向剪接位點的側翼序列進行學習，circDeep使用卷積神經網絡和雙向長短時記憶網絡對序列進行編碼，CRC對反向剪接位點的側翼區域特征建立基于卷積神經網絡的上下文回歸模型，JEDI使用深度雙向循環神經網絡編碼序列并通過交叉注意層構建反向剪接的深層相互作用模型。卷積神經網絡能夠獲得重要的序列局部模式來進行預測，但是無法識別每個剪接位點的位置信息。circDeep通過應用循環神經網絡學習序列信息，彌補了卷積神經網絡的不足，但是忽略了一些基本的信息(如剪接位點)。CRC雖然能通過深度學習方式識別環狀發生過程，但輸入特征中包含了基于統計的信息。JEDI相較于前幾種深度學習預測工具，只對序列剪接位點周圍的側翼區域進行建模，不需要其他特征信息，充分挖掘了序列的剪接位點信息及其深度相互作用信息，可以自動發現反向剪接的位點，而無需任何注釋，還能夠很好地保留形成環狀RNA的剪接位點信息和其他重要信息，因此在模型評估的各衡量指標中都取得了最好的表現。

表2 預測環狀RNA 工具優缺點Table 2 Advantages and disadvantages of circRNA prediction tools

無論是基于傳統機器學習方法還是深度學習方法，以上模型都是從序列中挖掘局部信息，但受限于知識和方法的不足，已經利用的序列信息(序列基序，ALU序列和剪接位點等)還是不足以完全地解釋RNA成環機制。序列的上下游調控信息、遠程調控信息、RNA與蛋白質互作信息和RNA結構等是現有工具未挖掘到的一些信息。如何更全面地挖掘信息并有效地表征，是環狀RNA識別工具開發的一個可能方向。

本文通過比較分析現有工具各自特征提取的側重點和方法的優劣，目的在于幫助大家在研究過程中選擇合適的工具，也希望能對開發出更好的預測環狀RNA的算法和工具有所啟發，從而推進對環狀RNA形成機制的研究和功能的探索。