文章

Reconstruction of full-length circular RNAs enables isoform-level quantification

作者單位

中國(guó)科學(xué)院北京生命科學(xué)研究院計(jì)算基因組學(xué)實(shí)驗(yàn)室趙方慶團(tuán)隊(duì)（第一作者：博士研究生鄭毅&助理研究員冀培豐）

發(fā)表期刊

Genome Medicine

發(fā)表時(shí)間

2019年01月19日

引言

“君子生非異也，善假于物也?！边@句話很適合該文的方法CIRI-full，它巧妙地應(yīng)用了i)環(huán)狀RNA的本性——環(huán)狀且較短（大部分約300~500bp）; ii)環(huán)化建庫(kù)以及更長(zhǎng)雙末端測(cè)序（PE250, PE300）——環(huán)狀RNA序列全覆蓋。

名詞解釋

FRJ：forward splice junction，前向剪接位點(diǎn)，比如 mRNA 是pre-mRNA 通過(guò)剪接體復(fù)合物通過(guò)前向剪切去除內(nèi)含子將外顯子連在一起的產(chǎn)物。

BRJ：back splice junction，反向剪接位點(diǎn)，比如環(huán)狀RNA通過(guò)剪接體復(fù)合物將供體3’末端外顯子與受體外顯子5’末端連接在一起的產(chǎn)物。

RO：reverse overlap，反向交疊，如圖灰色陰影所示，表示雙末端測(cè)序（paired-end sequencing） Read1 與 Read2 的5’末端或3’末端存在反向一致序列。

Cirexon：circRNA’s exon，構(gòu)成環(huán)狀RNA的外顯子

Isoform：來(lái)自同一基因組區(qū)間的環(huán)狀RNA轉(zhuǎn)錄本

FSG：forward splice graph，前向剪接圖

文章亮點(diǎn)

?  提出了識(shí)別環(huán)狀RNA全長(zhǎng)序列的方法CIRI-full

?  相比基于BSJ的方法，RO方法對(duì)低表達(dá)豐度的環(huán)狀RNA 更敏感

?  利用FSG的方法精確識(shí)別和量化環(huán)狀RNA的轉(zhuǎn)錄本（isoform）

?  6個(gè)物種（包括人，恒河猴，小鼠，大鼠，兔，雞）大腦全長(zhǎng)環(huán)狀RNA

?  相比BSJ，CIRI-full獲得的轉(zhuǎn)錄本（isoform）能過(guò)濾假陽(yáng)性的差異表達(dá)環(huán)狀RNA

數(shù)據(jù)資源

測(cè)序數(shù)據(jù)：PRJNA475651

軟件：https:///projects/ciri

環(huán)狀RNA研究的攔路虎——無(wú)法準(zhǔn)確捕獲全長(zhǎng)序列

環(huán)狀RNA 功能的重要性毋庸置疑，然而對(duì)其功能的探究仍舉步維艱，主要原因除了人們?nèi)跃窒抻谘芯烤€性RNA分子的思維模式，另一個(gè)主要原因就是無(wú)法大規(guī)模獲得環(huán)狀RNA的全長(zhǎng)序列并準(zhǔn)確定量。我們都知道，從序列同源以及物種進(jìn)化地角度能夠有效地探究一個(gè)遺傳分子可能存在地功能；另外，一個(gè)遺傳分子表達(dá)的準(zhǔn)確定量是反應(yīng)其在生物體內(nèi)發(fā)揮功能的有效途徑。然而，大規(guī)模從 RNA-seq 中識(shí)別全長(zhǎng) circRNA 的方法仍待開(kāi)發(fā)。這篇文章就提出了一個(gè)新的方法CIRI-full從RNA-seq中探究環(huán)狀RNA轉(zhuǎn)錄本的全長(zhǎng)序列以及表達(dá)定量。

實(shí)驗(yàn)基礎(chǔ)

當(dāng)插入片段長(zhǎng)度（即測(cè)序文庫(kù)片段大小）大于環(huán)狀RNA序列長(zhǎng)度，那么就可以觀察到雙末端 reads 的RO事件。

建庫(kù)：環(huán)化建庫(kù)，最后插入片段大小為300-800bp

測(cè)序平臺(tái)：Illumina HiSeq 2500 platform

Read長(zhǎng)度：250 x 2, 300 x 2

方法概述

該方法總共分為四部分，包括：

?  RO read 的檢測(cè)與檢驗(yàn)（圖A-D）

檢測(cè)雙末端 reads（paired-end reads）Read1與Read2的5’末端是否存在反向交疊RO（圖A-B），如果存在則作為候選RO-merged read；隨后分割比對(duì)RO-merged read 到參考基因組，最長(zhǎng)的作為位置標(biāo)簽——錨定位點(diǎn)anchor，區(qū)分異常以及無(wú)法比對(duì)的區(qū)域（圖C）；校正剪切信號(hào)GT/AG確定環(huán)狀RNA邊界（圖D）。

?  BSJ 以及 cirexon 檢測(cè)

利用 CIRI2 檢測(cè)反向剪切事件BSJ，并通過(guò)CIRI-AS識(shí)別BSJ中的可變剪接事件（single-splice event）。如果構(gòu)成BSJ的兩個(gè)reads只落在BSJ內(nèi)部的外顯子即cirexons，那么環(huán)狀RNA的全長(zhǎng)序列就能由cirexons線性重構(gòu)。

?  組裝RO與BSJ的reads（圖E）

i） RO-merged read存在3’-RO或兩個(gè)末端落在同一個(gè)cirexon上

ii） 每個(gè)BSJ的所有雙末端reads都落在cirexons

以上兩種情況可以直接確定環(huán)狀RNA全長(zhǎng)，而剩余的reads，CIRI-full結(jié)合兩者優(yōu)點(diǎn)——BSJ確定邊界，RO-reads確定內(nèi)部cirexons——進(jìn)一步組裝全長(zhǎng)。

?  轉(zhuǎn)錄本isoform重構(gòu)與定量（圖F）

CIRI-full用BSJ cirexons以及RO-merged reads都構(gòu)建了一個(gè)前向剪接圖FSG，利用深度優(yōu)先搜索算法（adapted depth-first search）將FSG分解成通道（path）——候選轉(zhuǎn)錄本（isoform），然后用蒙特卡洛仿真（Monte Carlo simulation）方法模擬每個(gè)通道（轉(zhuǎn)錄本）的表達(dá)豐度，最后用近乎窮舉的方式尋找每一條通道相對(duì)表達(dá)豐度（設(shè)置了每個(gè)環(huán)狀RNA最多10個(gè)轉(zhuǎn)錄本以提高效率）。

方法詳述

1. RO read 的檢測(cè)與檢驗(yàn)

5’-RO識(shí)別策略 對(duì)于每對(duì) Read1 與 Read2，Read1 與 Read2 其中一個(gè)的 5’ 末端的前 10 bp 被分成步長(zhǎng)為 1 bp、窗口為 8 bp 的三份（如示例圖，假設(shè)為 Read1 的 5’ 端）。

這些子序列被用來(lái)作為種子搜索 Read2 的 5’ 末端，一旦所有種子匹配堿基數(shù) >= 7bp，那么從 Read2 的 5’ 最后一個(gè)堿基到 Read1 匹配到的那個(gè)位點(diǎn)的序列被提取出來(lái)，然后再比對(duì)到 Read1 上。如果滿足以下條件，則這對(duì) reads 被認(rèn)為 5’-RO：reads pair上的長(zhǎng)度 >= 13 bp，堿基一致性 >= 95%。隨后，Read1 與 Read2 根據(jù)交疊部分合并為一條長(zhǎng) read 作為候選 RO-merged read以待進(jìn)一步驗(yàn)證。3’-RO與此類(lèi)似。

過(guò)濾假陽(yáng)性 候選RO-merged read 通過(guò)BWA-MEM比對(duì)確定在參考基因組上的位置，然后將比對(duì)長(zhǎng)度進(jìn)行排序，其中最長(zhǎng)且比對(duì)得分 > 15 的序列被作為anchor（錨定位點(diǎn)）來(lái)確定 reads 的位置。如果anchor 的兩頭在基因組上的范圍 < 100kbp，那么將計(jì)算ro-merged="" read="" 比對(duì)上基因組的長(zhǎng)度；長(zhǎng)度如若超過(guò)了="" read 的一半，將用于再次做局部比對(duì)（local="" realignment）來(lái)確定ro-merged="" read在基因組上的精確位置——尋找bsj（如果沒(méi)有將用anchor），并用動(dòng)態(tài)規(guī)劃算法計(jì)算無(wú)法比對(duì)以及異常比對(duì)的片段。之后將過(guò)濾假陽(yáng)性ro-merged="" read，原則包括：read不包含bsj以及read兩端等長(zhǎng)子序列比對(duì)位置沒(méi)有落在anchor="">

校正GT/AG 由于BWA-MEM無(wú)法獲得GT/AG剪切信號(hào)以及精確的剪切邊界，另外RO-merged reads可能包含套索結(jié)構(gòu)（lariat structure）。因此，需要對(duì)每條候選RO- merged read檢測(cè)和校正GT/AG位點(diǎn)，如果不存在將丟棄。另外，對(duì)每條read剪接位點(diǎn)（junction site）上下游5-bp比對(duì)質(zhì)量進(jìn)行檢測(cè)，如果存在gap或錯(cuò)配，也將丟棄不用。

2. BSJ 以及 cirexon 檢測(cè)

BSJ 用CIRI2檢測(cè)，而B(niǎo)SJ中的可變剪切事件（single-splice events）用CIRI-AS來(lái)推測(cè)。BSJ的檢測(cè)方法在所有識(shí)別軟件中都比較類(lèi)似，簡(jiǎn)單描述就是截取read的5’與3’末端的子序列比對(duì)到參考基因組上，如果比對(duì)方向相反將作為候選circRNA。詳細(xì)請(qǐng)?jiān)L問(wèn)CIRI 以及CIRI2原文。每個(gè)BSJ中的cirexons通過(guò)CIRI-AS捕獲的剪切事件推測(cè)。

3. 組裝RO與BSJ的reads

基于所識(shí)別到的RO-merged reads 以及BSJ邊界信息，利用以下原則構(gòu)建全長(zhǎng)環(huán)狀RNA。

?  如果RO-merged read存在3’-RO或兩個(gè)末端落在同一個(gè)cirexon上，那么該readd被認(rèn)定為一個(gè)全長(zhǎng)環(huán)狀RNA；否則等待進(jìn)一步組裝。

?  如果每個(gè)BSJ的所有雙末端reads都落在cirexons上，那么該BSJ的全長(zhǎng)環(huán)狀RNA由cirexons線性重構(gòu)；否則BSJ以及cirexons等待進(jìn)一步組裝。

?  RO與BSJ實(shí)際上存在互補(bǔ)信息——BSJ獲得環(huán)狀RNA精確的邊界，而RO識(shí)別內(nèi)部結(jié)構(gòu)。

CIRI-full結(jié)合兩者優(yōu)點(diǎn)，將等待進(jìn)一步組裝的RO-merged reads以及BSJ reads根據(jù)BSJ進(jìn)行排序聚類(lèi)，如果兩種read被發(fā)現(xiàn)落在同一個(gè)BSJ內(nèi)，那么這些reads將被用來(lái)重構(gòu)全長(zhǎng)環(huán)狀RNA。另外，RO-merged reads被用來(lái)確定BSJ未識(shí)別的額外cirexons，如果所有reads都落在cirexons上，那么全長(zhǎng)將通過(guò) cirexon 線性組裝；如若不然，BSJ 所識(shí)別的 cirexons 將標(biāo)記為部分重構(gòu)的 circRNA。

4. 轉(zhuǎn)錄本isoform重構(gòu)與定量

對(duì)每一個(gè)全長(zhǎng)環(huán)狀RNA，CIRI-full用BSJ cirexons以及RO-merged reads都構(gòu)建了一個(gè)前向剪接圖FSG，其中節(jié)點(diǎn)表示cirexons，邊表示兩個(gè)外顯子之間存在前向連接。

理論上，F(xiàn)SG包含了環(huán)狀RNA所有存在的轉(zhuǎn)錄本，然而值得注意的是，由于環(huán)狀RNA本身成環(huán)特性，所以FSG實(shí)際上是一個(gè)閉環(huán)。利用深度優(yōu)先搜索算法（adapted depth-first search）將FSG分解成通道（path），迭代地從每個(gè)節(jié)點(diǎn)開(kāi)始并在斷點(diǎn)或起始節(jié)點(diǎn)（cirexon）結(jié)束；其中短通道被合并成更長(zhǎng)的通道，而冗余的通道被過(guò)濾，為了避免過(guò)多可能的假陽(yáng)性，CIRI-full設(shè)置了固定的最大通道數(shù)（默認(rèn)10）。

詳細(xì)步驟：

1> FSG的邊包括4大類(lèi)：

?  BSJ

?  Phasing FSJ(forward splice junction)，表示只有一個(gè)環(huán)狀轉(zhuǎn)錄本

?  Co-occurred FSJs，表示剪切事件的數(shù)目與RO reads的數(shù)目相同

?  剩余的FSJ

其中，包含phasing FSJ和Co-occurred FSJs的通道優(yōu)先級(jí)最高，被稱(chēng)為phased isoforms。而對(duì)第四類(lèi)剩余的FSJs，根據(jù)節(jié)點(diǎn)的測(cè)序深度排序，只有top10被保留。所有保留的通道將作為候選轉(zhuǎn)錄本。

2> 為了確定每條通道的相對(duì)豐度，用蒙特卡洛仿真（Monte Carlo simulation）方法根據(jù) RNA-seq 插入片段的長(zhǎng)度分布（通過(guò) paired-end reads 在參考基因組上的距離推測(cè)）模擬每條通道的 BSJ-reads 分布。

3> 為了量化每條通道的相對(duì)豐度，CIRI-full 使用了窮舉法（approximate exhaustive search algorithm）：先為每條通道設(shè)置一個(gè)相對(duì)豐度（正整數(shù)），并且所有通道豐度和等于BSJ-reads總數(shù)?；诔跏枷鄬?duì)豐度以及模擬BSJ-reads，每條通道節(jié)點(diǎn)與邊的累積豐度被計(jì)算；根據(jù)mapped BSJ-reads計(jì)算節(jié)點(diǎn)與邊的真實(shí)豐度；隨后，計(jì)算模擬值與真實(shí)值的距離。隨后迭代地跟新模擬值，知道距離趨于收斂。

通過(guò)以上步驟，最終獲得每個(gè)環(huán)狀轉(zhuǎn)錄本的相對(duì)豐度。