環形RNA全長組裝與定量工具的研究進展

謝宏文，吳 靜,2，宋曉峰 *

(1.南京航空航天大學 自動化學院，南京 211106；2.南京醫科大學 生物醫學工程與信息學院，南京 211166)

環形RNA是一種特殊的內源性非編碼RNA，其由前體RNA(Pre-mRNA)經過反向剪接形成，呈封閉環狀結構，不具有5’末端帽子和3’末端poly(A)尾巴。研究發現環形RNA參與許多重要的生物學過程，如作為miRNA海綿體，參與基因調控，編碼蛋白等。環形RNA的反向剪接位點是鑒別和定量環形RNA的關鍵，而環形RNA可以通過可變剪接產生不同的環形轉錄本，這些轉錄本的全長序列信息以及精確定量對于進一步研究環形RNA功能具有重要作用。生物信息學方法因其能夠高效便捷的處理高通量RNA-seq數據，被普遍用來鑒別和分析環形RNA。本文介紹了真核生物環形RNA的產生及功能，對環形RNA檢測以及全長組裝和定量方面的研究工具進行了綜述。

1 環形RNA的產生及其功能

1.1 環形RNA的產生

Sanger等在20世紀70 年代發現某些高等植物中存在可致病的單鏈環狀類病毒，這是人類首次發現環形RNA[1]。但在過去的幾十年中，人們把環形RNA看作剪接副產物，環形RNA的研究進展十分緩慢。近年來，隨著第二代高通量測序技術的出現及環形RNA分子純化方法的運用和發展，人們可以重新認識和研究環形RNA。

環形RNA不具有5’末端帽子和3’末端poly(A)尾巴，是以反向剪接的方式形成的一種共價環狀結構。環形RNA主要來源于蛋白編碼基因，大多由其中的外顯子衍生而來并積累在細胞質中，少部分為包含外顯子和內含子序列的環形RNA，以及僅來自蛋白編碼基因內含子區域的環形RNA[2]。環形RNA的產生機制可分為內含子環化和外顯子環化。內含子環化發生在外顯子剪接過程中，該過程會產生內含子套索結構的中間產物，這些套索在剪接完成后仍然存在，形成內含子環形RNA(見圖1a)，如ci-ankrd52由內含子套索結構形成，可以與RNA合成酶Ⅱ結合促進其轉錄作用[3]；外顯子環化又可分為三種，即內含子配對驅動環化、套索驅動環化和RNA結合蛋白驅動環化[4]。內含子配對驅動環化將環化外顯子的側翼內含子互補配對，使環化外顯子的剪接供體和剪接受體位置更接近，形成環狀結構，然后將環狀結構中的內含子切除后形成環形RNA(見圖1b)，如秀麗隱桿線蟲側翼內含子上反向互補重復序列可以使轉錄本形成發夾結構，促進外顯子環化[5]；套索驅動環化發生在外顯子跳躍過程中，該過程可以使下游5’剪接供體與上游3’剪接受體結合，形成包含外顯子的套索結構，再通過套索內的剪接切除內含子，形成環形RNA(見圖1c)，如Barrett對酵母mrps16基因進行質粒表達，通過刪除剪接位點證明了包含外顯子的套索結構對于酵母細胞中環形RNA的生成是必要的[6]；RNA結合蛋白可以結合到環化外顯子兩側的內含子的結合位點上，拉近兩端的內含子促進環化，形成環形RNA(見圖1d)，如可變剪接調控因子quaking可以與環化外顯子側翼內含子上的quaking結合位點結合，促進兩側內含子相互靠近，連接成環[7]。

圖1 環形RNA環化模型

1.2 環形RNA的功能

由于環形RNA獨特的剪接方式及其環狀結構，曾經被認為是剪接副產物，不具有生物功能。近年來的研究表明環形RNA廣泛存在、結構穩定并且具有組織特異性，在生物體生長發育過程中發揮重要作用[8]。以下將從環形RNA與miRNA相互作用、調控基因表達、與RNA結合蛋白相互作用、翻譯蛋白質四個方面介紹環形RNA的功能。

1.2.1 環形RNA與miRNA相互作用

近期研究表明，環形RNA可以作為miRNA分子海綿，通過競爭性的結合miRNA，降低miRNA對其靶分子的抑制作用，進而調控基因表達水平(見圖2 a)。Hansen等人研究發現CDR1as(Antisense to the cerebellar degeneration-related protein 1 transcript)含有miR-7的超過70個保守結合位點，可以調節miR-7靶基因的表達[9]。Memczak等人在CDR1as敲除的人類細胞中觀察到miR-7靶點的下調，認為CDR1as是作為miRNA海綿來減弱miRNA介導的反應[10]。You等人經過研究認為大多數circRNA并不充當miRNA海綿，并不比其同源mRNA具有更多的miRNA結合位點[11]。

1.2.2 環形RNA與RNA結合蛋白相互作用

研究表明，RNA結合蛋白(RBPs)如Argonaute、RNA聚合酶II和MBL可與環形RNA結合。一些環形RNA可以儲存、分類或定位RBP，并且可能像它們調節miRNA活性一樣，通過作為競爭元素來調節RBP的功能(見圖2a)[12]。環形RNA獨特的三級結構在RNA或RBP復合物的組裝中起重要作用，其可以作為腳手架(Scaffolding)為蛋白質與RNA、蛋白質與DNA、蛋白質與蛋白質之間的相互作用提供平臺[13]。

1.2.3 環形RNA調控基因表達

環形RNA可以認為是一種選擇性剪接產物，因此環形RNA可能在選擇性剪接水平上起到調節基因表達的作用。環形RNA可與U1 snRNP相互作用，增強其親本基因的表達(見圖2 b)[14]。Burd等人通過研究推測cANRIL(Circular ANRIL)具有轉錄調控功能，cANRIL的形成降低了ANRIL轉錄本的表達，從而對編碼基因INK4/ARF的轉錄進行調控[15]。Yang等人發現環形RNA circ-ITCH可以作為海綿體結合miR-17和miR-224，上調靶基因p21和PTEN的表達[16]。

1.2.4 環形RNA翻譯蛋白質

環形RNA缺乏帽依賴性翻譯的必需元件，例如5'帽子和poly(A)尾，而有些環形RNA也可以翻譯產生蛋白(見圖2 c)。環形RNA的非帽依賴性翻譯可以通過內部核糖體進入位點(IRES)或在5'非翻譯區(UTR)中加入m6A修飾后發生[17]。Legnini等人證明circ-ZNF609包含開放閱讀框(ORF)，可以通過非帽依賴性翻譯編碼蛋白控制肌細胞的增殖[18]。Yang等人通過Northern blotting和質譜檢測等技術發現circ-FBXW7能夠編碼蛋白FBXW7-185 aa，可以控制癌細胞周期和增殖[19]。

2 環形RNA檢測工具

生物信息學方法因其能夠高效便捷的處理高通量RNA-seq數據，被普遍用來鑒別和分析環形RNA。為了研究和探索環形RNA的普遍性質和多樣功能，研究者們開發了各種計算工具從RNA測序數據中檢測環形RNA。這些工具根據檢測方法可以為兩類：基于注釋檢測和不基于注釋檢測。

2.1 基于注釋的環形RNA檢測工具

基于注釋的環形RNA檢測工具根據鑒定策略可分為兩類，第一類是基于“偽參考序列”的檢測策略，第二類是基于“比對片段”的檢測策略。

“偽參考序列”檢測策略需要根據基因注釋信息來構建環形RNA的偽序列，通過偽序列檢測環形RNA?；谶@種策略的檢測工具有KNIFE和PTESFinder等[20-21]。其中KNIFE應用較為廣泛，該工具首先從基因組注釋信息中構建所有可能的無序“外顯子-外顯子”連接序列，在Bowtie2的幫助下，將讀取的數據分別映射到基因組、rRNA序列、線性“外顯子-外顯子”連接序列和無序“外顯子-外顯子”連接序列[22]。真實的BSJ讀段必須覆蓋指定的核苷酸數量且不能比對到基因組和rRNA序列以及規范的線性“外顯子-外顯子”連接序列。KNIFE增加了從頭分析模塊檢測來自未注釋剪接位點的環形RNA，但Zeng等人認為這種從頭檢測方法是基于統計學進行推斷，不能提供準確的斷點[23]。

“比對片段”的檢測策略的大致思路便是將測序讀段與參考基因組比對，跨越BSJ(Back-spliced junction)的讀段被分裂成片段并以相反的順序與參考基因組對齊，根據這些讀段去識別反向剪接位點。基于該策略的工具有CIRCexplorer、UROBORUS和DCC等[24-26]。其中CIRCexplorer使用較為廣泛，該工具首先使用TopHat將讀段比對到參考基因組，然后提取未映射的讀段與TopHat Fusion的結果進行比對，若讀段映射到染色體上的順序相反且映射位置與已知基因注釋中的外顯子剪接邊界一致，便得到環形RNA的反向剪接位點。

相較于不基于注釋的檢測方法，基于注釋的檢測方法能夠更可靠的檢測反向剪接位點，但無法應用于缺乏基因組注釋的物種上。

2.2 不基于注釋的環形RNA檢測工具

不基于注釋的環形RNA檢測工具有circRNA finder、find_circ和CIRI等[10,27-28]。其中CIRI使用較為廣泛，該工具使用BWA-MEM將讀段映射到參考基因組[29]。與上述工具需要從未映射上的讀段中提取錨序列來檢測反向剪接不同，BWA-MEM可以自動確定環形RNA派生讀段的斷點。CIRI對BWA-MEM比對完成的結果進行兩次掃描計算，過濾掉假陽性的反向剪接位點，精度和靈敏度優于其他檢測工具，且運算時間耗費不大，實現了更好的平衡性能。相較于基于注釋的檢測方法，不基于注釋的檢測方法應用范圍更廣，能夠檢測新的候選環狀RNA，但需要提高檢測和過濾過程中的靈敏度和準確性。

3 環形RNA全長組裝與定量工具

研究表明環形RNA在同一個反向剪接位點內部可通過可變剪接形成多個不同的轉錄本(見圖3)[19]。若要深入研究環形RNA的功能特性，必須獲取環形RNA轉錄本內部全長序列信息以及對不同環形RNA內部可變剪接產物進行精確定量。研究者們開發了多種計算工具用于環形RNA內部結構的探索、全長序列的組裝及表達量的分析，目前已知的工具(見表1)。

表1 環形RNA下游分析計算工具

圖3 環形RNA內部可變剪接示意圖

3.1 環形RNA內部結構探索及全長序列組裝工具

為了更好的揭示環形RNA的調控機制，需要深入研究環形RNA內部復雜的選擇性剪接結構，構建環形RNA全長序列。目前識別環形RNA內部可變剪接結構的工具主要有CIRI-AS、FUCHS和CircSplice[30-32]；針對于環形RNA全長序列組裝的工具主要有CIRCexplorer2、CIRI-full、CircAST、isoCirc和CIRI-long[33-37]。

Gao等人開發了CIRI-AS工具研究環形RNA內部的剪接結構[30]。該工具首先利用CIRI找到環形RNA反向剪接位點以及跨越反向剪接位點的讀段，然后利用雙端測序信息，分析跨越反向剪接讀段的另一端讀段是否支持環形RNA內部的前向剪接位點，根據支持前向剪接位點讀段信息構建前向剪接圖，列出可能發生可變剪接的外顯子和所有可能的路徑，以檢測環形RNA內部的可變剪接事件。通過計算與實驗驗證相結合的手段，揭示了外顯子跳躍，5’或3’可變剪接位點以及內含子保留這四種可變剪接事件在環形RNA中普遍存在。CIRI-AS證實了環形RNA內部存在可變剪接事件，但并沒有給出環形轉錄本內部的組成。

Metge等人開發了工具FUCHS對環狀RNA內可變剪切等信息進行分析和解讀[31]。FUCHS基于長讀段測序(大于150 bp),在運行前需要用戶手動輸入比對和環形RNA位點檢測信息。FUCHS從輸入的比對信息中提取出嵌合比對的讀段信息，識別出同一反向剪接位點內的前向可變剪接事件。類似于CIRI-AS，FUCHS利用雙端測序信息來剪接驗證環形RNA內部的可變剪接，篩選掉只有一段跨越反向剪接位點的雙端測序讀段實現對假陽性環形RNA的過濾。

Feng等人開發了CircSplice工具用于識別環形RNA內部的可變剪接事件[32]。該工具通過STAR進行比對，通過GT-AG和CT-AC兩種剪切位點過濾并識別環形RNA反向剪接位點，保留兩端都落在反向剪接位點之間的讀段，根據一端比對到反向剪接位點，另一端比對到反向剪接位點內部的讀段對外顯子跳躍、內含子保留、5’選擇性剪接，3’選擇性剪接四種可變剪接類型進行判斷，計算支持可變剪接事件的讀段數量，最后根據基因組坐標融合不同的可變剪接事件，并比較不同樣本間的可變剪接差異。相較于上文提到的CIRI-AS，該方法支持不同樣本環形RNA可變剪接差異的比較。

CIRCexplorer2是環形RNA檢測軟件CIRCexplorer升級版[33]。CIRCexplorer2將Tophat未比對上但映射到Tophat-fusion的讀段重新比對到已知和從頭組裝的注釋上，可以檢測來自注釋或新外顯子邊界的反向剪接位點；增加了檢測環形RNA中的可變剪接事件功能用于確定環形RNA全長，該工具通過比較分析ploy(A)+和poly(A)-的數據集識別環形RNA內部可變剪接事件，重構環形RNA轉錄本序列。雖然CIRCexplorer2能夠通過這種方法給出環形RNA的全長轉錄本，其使用的工具卻是線性轉錄本的組裝工具Cufflinks，會得到錯誤的環形RNA序列，給計算帶來偏差。

Zheng等人開發了工具CIRI-full，提出了一種環形RNA轉錄本組裝的新方法[34]。由于二代測序數據讀長較短，限制了研究者們只能定位環形RNA的反向剪接位點，難以獲得環形RNA內部的全長結構信息。Zheng等人于是增加測序讀長到250 bp，在較長的讀長下，若某一雙端測序讀段來自與環形RNA，兩端讀段會存在反向重疊區(Reverse overlap,RO)，該特征不僅可以識別環形RNA的反向剪接位點，而且可以判斷該序列是否覆蓋整個環形RNA，從而獲取環形RNA的全長序列。

2020年，Wu等人開發了組裝環形RNA全長序列的工具CircAST[35]。該工具首先利用Tophat的比對結果找到跨越外顯子發生剪接的讀段信息，后根據用戶通過環形RNA識別工具(如UROBORUS、CIRCexplorer2、CIRI2等)檢測到的反向剪接位點信息以及基因組注釋文件中的外顯子邊界信息，對每一個基因位點構建反向剪接事件的有向無環圖(Directed acyclic graph，DAG)，圖中的起點和終點對應環形RNA中的兩個發生反向剪接的外顯子。相較于傳統的線性轉錄本拼接算法，CircAST若檢測到某一基因位點存在超過一個反向剪接事件，則構建多個剪接圖，因此可以進行準確拼接和組裝。CircAST利用擴展最小路徑覆蓋算法(Extended minimum path cover，EMPC)結合測序讀段信息計算出最優拼接方式，實現環形RNA的全長序列的拼接組裝。

Xin等人提出了識別環形RNA全長序列的工具isoCirc[36]。由于短讀段測序不能通過實驗確定環形RNA全長序列的組成，該工具將數據進行線性RNA消化，并進行滾環擴增，提高低豐度環形RNA的比例，后利用納米孔長讀段測序技術得到包含環形RNA全長序列的長度段測序數據。isoCirc基于上述長讀段測序數據，識別出其中重復共有的片段，并將兩個相同的重復共有片段連接起來。將連接片段比對到參考基因組以識別環形RNA反向剪接位點和環形RNA內部前向剪接位點信息。通過多重策略(如比對質量，BSJ/FSJ 精確度等)對比對片段進行過濾，之后便得到環形RNA反向剪接位點及全長序列信息。Xin等人使用isoCirc對12種人類組織及人類HEK293細胞系進行測試，一共檢測到107 147個環形RNA全長序列，其中包含40 628個長度大于500 nt的環形RNA亞型。isoCirc工具利用納米孔測序技術，彌補了短讀段測序的缺陷，提供了一種研究環形RNA全長序列的新策略。

Zhang等人同樣基于納米孔測序技術，開發了CIRI-long工具[37]。在環形RNA建庫過程中，CIRI-long與isoCirc稍有不同，在對數據進行線性RNA消化和環形RNA滾環擴增之后，針對于長度更長cDNA序列加入了片段長度篩選流程。經過納米孔測序后得到長讀段測序數據，利用k-mer匹配方法得到長度段中的重復片段，并使用偏序比對修正測序錯誤，然后將重復片段比對到參考基因組，結合剪接信號信息，識別反向剪接位點和環形RNA內部正向剪接信息，實現了環形RNA的識別及全長序列的重建。Zhang等人評估了該方法在模擬數據上的效果，并與Illumina測序和實時定量RT-PCR進行了比較，驗證了該方法的準確性。

3.2 環形RNA定量工具

為了深入研究環形RNA的調控機制以及不同環形RNA轉錄本之間的表達差異，需要對環形轉錄本進行精確定量。目前對于環形RNA進行定量的工具主要有sailfish-cir、CIRCexplorer3、CIRIquant、CIRI-full和CircAST[38-40]。

Li等人開發了工具sailfish-cir對環形RNA的表達量進行計算[38]。先前的研究主要根據檢測到的反向剪接讀段數量來量化環形RNA的表達，這種檢測方法由于反向剪接讀段較少，精度較低。因此Li等人將檢測到的環形RNA在反向剪接位點處切開，將3’端的序列追加到5’端，構成偽線性轉錄本，并加入到真實線性轉錄本集合中構成混合轉錄本集合。后利用對線性轉錄本定量的工具sailfish對混合轉錄本進行定量分析。sailfish-cir基于期望最大化模型，通過迭代將讀段分配到不同轉錄本上，并估計這些轉錄本的表達量。與基于反向剪接讀段計數的直接定量法相比，sailfish-cir與qRT-PCR定量結果有更強的相關性。

Ma等人對CIRCexplorer分析流程進行了升級，開發了環形RNA與線性RNA定量比較的工具CIRCexplorer3-CLEAR[39]。該方法提出了一種定量轉錄本表達量的參數FPB(Fragments per billion mapped bases)，該參數相較于FPKM能夠不受測序讀段長度的影響。該工具首先利用HISAT2將測序數據比對到參考基因組，利用StringTie對線性轉錄本進行組裝，計算跨越前向剪接位點的測序片段得到線性轉錄本的表達FPBlinear；然后利用TopHat-Fusion處理未比對上的測序片段，通過計算跨越反向剪接位點片段得到環形轉錄本的表達量FPBcirc，通過環形轉錄本表達量與對應線性轉錄本表達量的比值CIRCscore(FPBcirc/FPBliner)衡量環形RNA的相對表達量。相較于其他基于跨越反向剪接位點讀段進行定量的方法，該工具提出了新的定量參數FPB，可以利用CIRCscore對環形RNA進行相對定量，具有更廣泛的適應性。

Zhang等人開發了鑒定和定量環形RNA的算法CIRIquant[40]。CIRIquant首先利用CIRI等工具從測序數據中識別環形RNA的反向剪接位點，然后基于環形RNA反向剪接位點構造環形轉錄本偽參考序列，重新映射讀段到偽參考序列。CIRIquant運用這種方法可以更準確的識別環形RNA的反向剪接序列，根據反向剪接讀段和前向剪接讀段的比例對環形RNA進行定量?？紤]到環形RNA建庫時RNase R處理存在偏差，CIRIquant結合未經RNase R處理過的數據利用高斯混合模型(Gaussian mixture model,GMM)對RNase R處理數據進行校正，使定量更為準確。

除了這些工具外，上文提到的工具CIRI-full和CircAST不僅可以組裝環形轉錄本全長序列，還可以對環形轉錄本進行表達量估計。CIRI-full在構建環形RNA的全長后，采用蒙特卡洛方法模擬不同環形RNA剪接產物讀段在全長序列上的分布，通過梯度下降法求得最優的表達量組合，實現了對不同環形轉錄本相對豐度的預測評估，并在模擬和真實數據上驗證了該方法的準確性。但該方法受到測序長度的限制，需要更高的測序成本和更先進的測序技術來獲得較長讀長的測序數據；CircAST在實現環形RNA全長序列組裝之后，通過期望最大化(Expectation maximization，EM)算法來估計每個組裝的環形RNA轉錄本的豐度。其在似然函數的構造等方面充分注意到環形RNA的結構特點，在絕對定量和相對定量中均表現出了良好的性能。

4 總結與展望

隨著深度測序技術和分子純化方法的發展，人們對于環形RNA的認識正在逐漸發展。在對環形RNA的研究中，計算方法因其在高通量RNA-seq數據分析中的便利性以及在分析環形RNA表達方面的優勢而起著重要的作用。由于環形RNA自身的結構特點和其特有的剪接事件，我們不能簡單地用線性RNA的相關工具來解決環形RNA的相關計算，需要開發更多針對于環形RNA的計算工具。為了更好的揭示環形RNA的特殊功能，需要對環形RNA內部全長序列進行準確重建。由于二代測序的讀段較短，基于二代測序數據重建環形RNA全長序列難度較大；基于納米孔長讀段測序技術構建環形RNA全長的工具將會是未來研究開發的方向。在環形RNA定量方面，現有的工具或基于線性RNA開發而來，或對于輸入數據存在特定要求。而在全轉錄組測序數據上由于環形RNA與線性RNA重疊部分的干擾，精準定量環形RNA轉錄本仍然是一個挑戰。開發基于全轉錄組測序數據環形RNA特定的計算工具是目前環形RNA研究中迫切需要解決的問題，這對于進一步研究環形RNA的生物學功能至關重要。隨著計算算法、測序技術、基因組學和生物信息學的發展，相信會有更多環形RNA相關計算工具的出現，促進對于環形RNA機制和功能的深入研究。