不依賴于剪接位點(diǎn)信號的高精度轉(zhuǎn)錄組序列比對算法①

張 勇, 徐 云

1(中國科學(xué)技術(shù)大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院, 合肥 230027)2(中國科學(xué)技術(shù)大學(xué) 安徽省高性能計算重點(diǎn)實(shí)驗(yàn)室, 合肥 230027)3(國防科學(xué)技術(shù)大學(xué) 高性能計算協(xié)同創(chuàng)新中心, 長沙 410073)

不依賴于剪接位點(diǎn)信號的高精度轉(zhuǎn)錄組序列比對算法①

張 勇1,2, 徐 云1,3

1(中國科學(xué)技術(shù)大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院, 合肥 230027)2(中國科學(xué)技術(shù)大學(xué) 安徽省高性能計算重點(diǎn)實(shí)驗(yàn)室, 合肥 230027)3(國防科學(xué)技術(shù)大學(xué) 高性能計算協(xié)同創(chuàng)新中心, 長沙 410073)

高通量轉(zhuǎn)錄組測序技術(shù)已經(jīng)發(fā)展成為分析不同細(xì)胞中選擇性剪接事件的最有效方法, 其測序數(shù)據(jù)處理的第一步是將數(shù)以百萬的測序片段準(zhǔn)確地比對到參考序列上, 稱之為轉(zhuǎn)錄組序列比對. 現(xiàn)有的比對工具基本上都是依賴于經(jīng)典的剪接位點(diǎn)信號, 一定程度上限制了轉(zhuǎn)錄組測序技術(shù)發(fā)現(xiàn)全新剪接位點(diǎn)的能力. 為此, 我們設(shè)計了一種不依賴于剪接位點(diǎn)信號的轉(zhuǎn)錄組序列比對方法RNAMap, 該方法按照重疊種子方式劃分測序片段, 使用帶有左右錨點(diǎn)的窗口掃描參考序列, 找出種子中含有的剪接位點(diǎn). 計算實(shí)驗(yàn)表明, RNAMap精確度高達(dá)95%, 召回率也明顯優(yōu)于其他算法.

選擇性剪接; 高通量轉(zhuǎn)錄組測序; 滑動窗口; 剪接位點(diǎn)

真核生物的基因是斷裂基因, 由內(nèi)含子序列和外顯子序列組成, 選擇性剪接是一種重要的轉(zhuǎn)錄后修飾過程, 在此期間, 前體RNA中的一個或多個內(nèi)含子片段被剪切除去, 然后剩余的外顯子拼接稱為成熟的mRNA, 如圖1. 選擇性剪接使得基因能夠產(chǎn)生多樣的轉(zhuǎn)錄本, 而且人類基因組中90%以上的多外顯子基因會發(fā)生選擇性剪接. 相關(guān)研究表明, RNA剪接發(fā)生異常與人類的許多疾病密切相關(guān)[1].

定性和定量研究轉(zhuǎn)錄組的傳統(tǒng)方法是構(gòu)建cDNA或表達(dá)序列標(biāo)簽(EST)文庫, 然后通過Sanger測序進(jìn)行后續(xù)分析. 但是, 因?yàn)镾anger測序技術(shù)的成本較高且通量較低, 所以這種方法十分昂貴和低效. 隨著下一代測序技術(shù)(next-generation sequencing, NGS)的迅猛發(fā)展, 高通量轉(zhuǎn)錄組測序(RNA-seq)在分析全基因組的剪接信息, 尤其是選擇性剪接事件方面展現(xiàn)出了極佳的性能, 并已經(jīng)發(fā)展成為研究剪接轉(zhuǎn)錄本的最有效的技術(shù)[2]. 目前, RNA-seq已能夠應(yīng)用于疾病的臨床診斷; 此外, 在基礎(chǔ)生物學(xué)研究中RNA-seq也有廣泛的應(yīng)用, 如分析不同的基因在不同階段的表達(dá)情況.

圖1 真核細(xì)胞基因結(jié)構(gòu)圖

RNA-seq分析軟件的一項(xiàng)重要功能便是重建剪接之前的mRNA在細(xì)胞中的形態(tài), 此外, 還應(yīng)該能夠評估每一種剪接異構(gòu)體的表達(dá)水平. 然而, 所有分析過程的第一步都是要將RNA-seq中得到的測序片段(reads)比對到基因組上的原始位置, 而這些短片段的長度從數(shù)十堿基到數(shù)百堿基不等, 數(shù)量有幾十萬甚至幾百萬和上千萬, 所以, 比對的過程是極其耗時的.

事實(shí)上, 如果測序片段完全來自于外顯子序列,那么常規(guī)的序列比對工具(BWA[3]、Bowtie[4]等)便可以應(yīng)對這種比對工作. 但是, 有大量的短片段是來自于兩個甚至多個外顯子序列, 在人類基因組中兩個外顯子序列一般間距20bp～500000bp, 這遠(yuǎn)遠(yuǎn)超過了常規(guī)序列比對工具處理的范圍. 因此, 研究的主要問題便是如何將跨越剪接位點(diǎn)的測序片段快速且準(zhǔn)確地比對到參考序列上.

為了解決上述問題, 早期的策略是根據(jù)已有的基因組注釋文件, 利用常規(guī)序列比對工具將測序片段定位到基因組上. 雖然這種方法可以定位大部分的測序片段, 但它的局限性也不容忽視. 畢竟, 即使是目前人們研究的最為深入的人類基因組, 它的注釋文件仍然是不完整的, 所以上述策略是無法識別未在注釋文件中出現(xiàn)的全新的剪接位點(diǎn), 而這也就使RNA-seq喪失了發(fā)現(xiàn)新剪接異構(gòu)體的能力.

事實(shí)上, 近年來也相繼出現(xiàn)了一些不依賴于基因組注釋文件的RNA-seq序列比對工具, 比如SpliceMap、MapSplice、TopHat[5]、CRAC[6]、 OLego[7]和HISAT[8]等. 其中, TopHat系列軟件是最具有代表性,也是目前使用最廣的比對工具. 它采用外顯子優(yōu)先的策略, 整個比對過程分為兩個階段. 第一階段, 利用Bowtie將測序片段定位到參考基因組上, 這樣, 含有剪接位點(diǎn)的測序片段就會被過濾出來; 然后通過MAQ中的組裝模塊將成功定位的短片段組裝起來.經(jīng)過這一階段, 供體位點(diǎn)和受體位點(diǎn)的側(cè)翼序列拼接起來組成潛在的剪接序列, 作為下一階段的參考序列.第二階段, 將在第一階段中未成功定位的測序片段比對到上述由外顯子拼接成的序列上. 然而, TopHat在拼接外顯子序列時僅僅考慮經(jīng)典的剪接位點(diǎn)(GT/C-AG), 雖然目前已知的具有經(jīng)典信號的剪接位點(diǎn)占了絕大多數(shù), 但是有研究表明非經(jīng)典剪接位點(diǎn)的比例很有可能被低估了[9,10]. 因此, TopHat存在的主要問題是會遺漏具有非經(jīng)典剪接信號的測序片段. 其余幾個RNA-seq序列比對軟件雖然分別采用了各自不同的比對策略, 但在默認(rèn)情況下也都是依賴了經(jīng)典的剪接信號, 所以也具有與TopHat類似的缺陷. 盡管個別軟件可以通過設(shè)定參數(shù)來窮盡所有類型的剪接位點(diǎn)信號, 但是算法的復(fù)雜度較高. 總之, 目前依賴于經(jīng)典剪接位點(diǎn)信號的比對算法已經(jīng)發(fā)展的較為完善, 但尚缺乏對不依賴于剪接信號的比對算法的研究.

為了克服以上不足, 我們使用帶有左右錨點(diǎn)的窗口掃描參考序列, 設(shè)計了命名為RNAMap的轉(zhuǎn)錄組序列比對工具.

1 方法

RNAMap的執(zhí)行過程分為兩個階段. 第一階段, RNAMap嘗試?yán)贸Ｒ?guī)的序列比對工具將原始數(shù)據(jù)集中的所有的測序片段定位到參考基因組上. 在這一階段中, 完全來自于一個外顯子序列的測序片段可以被直接比對到基因組上, 這樣沒有比對上的測序片段就可能含有剪接位點(diǎn). 第二階段, 對于這些未比對上的片段, 利用兩個表來尋找其中的剪接位點(diǎn). RNAMap的執(zhí)行流程如圖2所示. RNAMap將測序片段劃分成幾個重疊的種子, 每個種子可發(fā)現(xiàn)一個剪接位點(diǎn), 這樣我們的方法就能夠處理含有多個剪接位點(diǎn)的測序片段.

1.1 初始比對階段

RNAMap使用Bowtie來處理RNA-seq的reads,將它們比對到參考序列上. 如果存在基因組注釋文件(文件中記錄了原基因組中外顯子序列的位置), 那么可以此文件為基礎(chǔ)生成轉(zhuǎn)錄組序列(只包含外顯子序列), 并將其作為參考序列. 采用這種策略, 一方面可以提高序列比對的敏感性和準(zhǔn)確性; 另一方面也可以加速比對的過程. 如果無法獲得有效的注釋文件, 那么RNAMap會選擇基因組作為參考序列.

即使以轉(zhuǎn)錄組作為參考序列, 也仍然會有一些reads無法成功地定位到參考序列上, 可能因?yàn)檫@些reads中被錯誤測序的堿基數(shù)超出了RNAMap設(shè)定的閾值, 另一個重要的原因是產(chǎn)生這些短片段的轉(zhuǎn)錄本信息并沒有記錄在注釋文件中. 此外, 由于假基因的存在[11], 也會有一些短片段被錯誤的定位到參考序列上.

圖2 RNAMap流程圖

1.2 分段比對階段

一個剪接位點(diǎn)可以將一個read分成兩個片段(segments), 但事實(shí)上, 這些segments并不是完全隨機(jī)的分布在基因組上. 如果我們不考慮一些特殊的情況,比如基因融合, 剪接位點(diǎn)分割一個read產(chǎn)生的segments應(yīng)該被定位到同一個染色體上, 并且滿足一定的距離限制, 對于人類及其他哺乳動物, 一般為20bp～500000bp. 如果先分別獨(dú)立的定位這些segments, 然后再根據(jù)位置限制條件進(jìn)行過濾會, 那么這樣會增加搜索空間. 因此, 如果在比對segments的同時添加有一定的限制條件, 那么既可以減小搜索空間, 又可以精簡后續(xù)的篩選過程. 正是基于這樣的考慮, 分段比對階段分為以下三個步驟.

1.2.1 劃分種子

將測序得到的片段劃分成互相重疊的種子, 例如,將長度為100bp的測序片段劃分為三個長為50bp的種子, 它們在原測序片段上的區(qū)間分別為[1,50]、[26,75]和[51,100]. 轉(zhuǎn)錄組測序深度可以保證每一個剪接位點(diǎn)至少會被一個種子所覆蓋.

1.2.2 種子定位

上一步產(chǎn)生的種子可以分為兩類: 一類種子不含有剪接位點(diǎn); 另一類種子含有剪接位點(diǎn), 并且我們假設(shè)它們僅含有一個剪接位點(diǎn). RNAMap調(diào)用Bowtie來比對所有的種子, 第一類種子可以被成功地定位到參考序列上, 第二類種子被過濾出來, 然后建立兩個表進(jìn)行索引, 一個為靜態(tài)表, 另一個為動態(tài)表. 此外, 種子中的剪接位點(diǎn)既可能出現(xiàn)在種子的前半段, 也可能出現(xiàn)在后半段. 下面僅討論剪接位點(diǎn)出現(xiàn)在后半段的情況, 以種子的前半段序列作為左錨點(diǎn), 后1/4序列作為右錨點(diǎn); 至于另外一種情況, 可以用一種對稱的方法來實(shí)現(xiàn). 為了能夠處理含有誤配的情況, RNAMap采用了與PerM[12]類似的單周期空間種子的方法.

① 靜態(tài)表

以種子的前半段(左錨點(diǎn))作為鍵, 以種子的標(biāo)識號作為值, 建立靜態(tài)表. 所有種子的鍵-值對信息都需要加入到表中, 并且靜態(tài)表一旦建立, 在比對的過程中都將保持不變, 如圖3.

② 動態(tài)表

動態(tài)表用于存儲右錨點(diǎn)序列與左錨點(diǎn)比對信息組成的鍵-值對. 有兩個窗口沿著參考基因組進(jìn)行滑動,窗口S的長度與靜態(tài)表鍵的長度相等, 用來查詢靜態(tài)表; 窗口D的長度一般為種子長度的1/4, 用來查詢動態(tài)表.

圖3 靜態(tài)表示例

當(dāng)窗口S中的序列在靜態(tài)表中查詢到時, 表明該種子的前半段可以定位到此處, 之后繼續(xù)向后延伸比對, 直至遇到第一個無法匹配的位點(diǎn), 最后以該種子的后1/4片段為鍵, 以前面的定位信息(包括標(biāo)志號、定位起點(diǎn)、比對的堿基數(shù))為值, 插入動態(tài)表.

當(dāng)窗口D中的序列在動態(tài)表中查詢到時, 表明該種子的后1/4片段可以定位到此處, 之后繼續(xù)向前延伸比對, 直至遇到第一個無法匹配的位點(diǎn), 最后結(jié)合動態(tài)表中對應(yīng)的值中保存的定位信息, 就可以判定種子是否能夠分段比對到參考基因組上的兩個位置. 如果前后兩部分比對的位置超出了預(yù)設(shè)的距離范圍, 則需要將動態(tài)表中對應(yīng)的鍵-值信息刪除. 此外, 每當(dāng)掃描完一個染色體的序列, 也需要將動態(tài)表中的信息全部清空, 以保證種子的兩部分定位到同一個染色體上.

1.2.3 種子整合

將種子的定位結(jié)果組合成reads的定位信息. 在這一過程中需要檢查種子比對位置的一致性, 從而將符合要求的種子組合成完整的reads.

2 實(shí)驗(yàn)結(jié)果

我們分別在模擬數(shù)據(jù)集和真實(shí)數(shù)據(jù)集上測試RNAMap的性能, 并與其他主流軟件進(jìn)行對比. 為了保證實(shí)驗(yàn)的可靠性和公平性, 所有的軟件都在同一臺計算機(jī)上運(yùn)行, 其基本的配置為Intel(R) Core(TM) i7-4770K CPU, 24G RAM, 64-bit Ubuntu 14.04 OS.

2.1 模擬數(shù)據(jù)集測試

我們使用FluxSimulator[13]軟件, 以人類基因組GRCh38及其注釋文件為基礎(chǔ), 隨機(jī)模擬產(chǎn)生了1000000條長為100bp的測序片段(reads), 每條序列的來源信息保存在BED格式的文件中, 因此我們可以計算各軟件比對結(jié)果的召回率(Recall Rate)和精確度(Precision), 結(jié)果如表1.

表1 各軟件的模擬數(shù)據(jù)集比對結(jié)果統(tǒng)計

由表1可知, 雖然CRAC和HISAT可以將絕大多數(shù)的測序片段(分別為97.46%和96.36%)比對到基因組上, 但是精確度比較低, 尤其是CRAC的精確度不足70%. TopHat2、SpliceMap、OLego、MapSplice、CRAC的召回率都比較低, 雖然HISAT的召回率達(dá)到80.17%,但是它的精確度也僅僅高于SpliceMap和CRAC. 因?yàn)槲覀兊哪M數(shù)據(jù)集是隨機(jī)產(chǎn)生的, 所以其中有大量的reads含有非經(jīng)典的剪接位點(diǎn), 但上述軟件無法處理此類情況, 從而造成其比對的質(zhì)量相對較低.

RNAMap的精確度高達(dá)96.15%, 雖然稍低于OLego的96.22%, 但是RNAMap的召回率達(dá)到了85.60%, 明顯高于其它幾個軟件. 這是因?yàn)樵诒葘Φ恼麄€過程中, RNAMap并未受經(jīng)典剪接位點(diǎn)信號的限制, 因而可以更準(zhǔn)確地將各類reads比對到參考序列上.

2.1 真實(shí)數(shù)據(jù)集測試

我們在73685727條長為100bp的真實(shí)測序片段數(shù)據(jù)集(來源于K562細(xì)胞系, 是一種人類的白血病細(xì)胞, GEO序列號為GSM1838573)上比較各個軟件的性能,結(jié)果如圖4所示.

由圖4可知, RNAMap可以將68647397條測序片段(93.16%)定位到參考序列上, 明顯優(yōu)于TopHat2、SpliceMap、OLego和MapSplice, 僅次于CRAC和HISAT. K562細(xì)胞系是一種癌變的細(xì)胞, 其選擇性剪接事件也與正常細(xì)胞不同, 因此其測序得到的reads中會含有更多類型的剪接位點(diǎn). 雖然我們無法統(tǒng)計真實(shí)數(shù)據(jù)集中正確的匹配位置, 但是根據(jù)模擬數(shù)據(jù)集的結(jié)果, 我們知道RNAMap的精確度在95%以上, 因此可以推斷RNAMap的整體性能較佳.

圖4 各軟件的真實(shí)數(shù)據(jù)集比對結(jié)果統(tǒng)計

3 結(jié)語

本文針對高通量轉(zhuǎn)錄組測序的序列比對問題, 提出了一種使用帶錨點(diǎn)的滑動窗口掃描參考序列的比對方法, 通過在模擬數(shù)據(jù)集和真實(shí)數(shù)據(jù)集上對算法的性能進(jìn)行測試, RNAMap無論是在召回率, 還是在精確度和片段匹配率上都表現(xiàn)出較優(yōu)的性能. 雖然在測序片段為100bp的大小為1000000的數(shù)據(jù)集上, RNAMap比目前最快的軟件HISAT大約多耗時40%, 但是時間仍然在可接受的范圍內(nèi), 而且可以獲得更高的精確度.該方法不借助經(jīng)典的剪接位點(diǎn)信號, 因此可以充分發(fā)揮RNA-seq的優(yōu)勢, 識別基因組注釋文件中沒有記錄的全新的剪接位點(diǎn). 接下來的工作, 一是需要提高RNAMap的比對速度; 二是解決含有多個剪接位點(diǎn)的種子的比對問題, 這樣可以降低對測序深度和種子重疊度的要求.

1 Nagao K, Togawa N, Fujii K, et al. Detecting tissue-specific alternative splicing and disease-associated aberrant splicing of the PTCH gene with exon junction microarrays. Human Molecular Genetics, 2005, 14(22): 3379–3388.

2 Mcgettigan PA. Transcriptomics in the RNA-seq era. Current Opinion in Chemical Biology, 2013, 17(1): 4–11.

3 Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 2009, 25(14): 1754–1760.

4 Langmead B, Trapnell C, Pop M, et al. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome biology, 2009, 10(3): R25.

5 Garber M, Grabherr MG, Guttman M, et al. Computational methods for transcriptome annotation and quantification using RNA-seq. Nature Methods, 2011, 8(6): 469–477.

6 Philippe N, Salson M, Commes T, et al. CRAC: An integrated approach to the analysis of RNA-seq reads. Genome Biology, 2013, 14(3): R30.

7 Wu J, Anczukow O, Krainer AR, et al. OLego: Fast and sensitive mapping of spliced mRNA-seq reads using small seeds. Nucleic Acids Research, 2013, 41(10): 5149–5163.

8 Kim D, Langmead B, Salzberg SL. HISAT: A fast spliced aligner with low memory requirements. Nature Methods, 2015, 12(4): 357–360.

9 Filichkin SA, Priest HD, Givan SA, et al. Genome-wide mapping of alternative splicing in Arabidopsis thaliana. Genome Research, 2010, 20(1): 45–58.

10 Parada GE, Munita R, Cerda CA, et al. A comprehensive survey of non-canonical splice sites in the human transcriptome. Nucleic Acids Research, 2014, 42(16): 10564–10578.

11 Kalyana-Sundaram S, Kumar-Sinha C, Shankar S, et al. Expressed pseudogenes in the transcriptional landscape of human cancers. Cell, 2012, 149(7): 1622–1634.

12 Chen Y, Souaiaia T, Chen T. PerM: Efficient mapping of short sequencing reads with periodic full sensitive spaced seeds. Bioinformatics, 2009, 25(19): 2514–2521.

13 Griebel T, Zacher B, et al. Modelling and simulating generic RNA-Seq experiments with the flux simulator. Nucleic Acids Research, 2012, 40(20): 10073–10083.

Highly Precise Transcriptome Sequence Alignment Algorithm Independent From Splice Site Signals

ZHANG Yong1,2, XU Yun1,312
(School of Computer Science and Technology, University of Science and Technology of China, Hefei 230027, China) (Key Laboratory of High Performance Computing of Anhui Province, University of Science and Technology of China, Hefei 230027, China)3(Collaborative Innovation Center of High Performance Computing, National University of Defense Technology, Changsha 410073, China)

RNA-seq has become the most effective method of analyzing alternative splicing events in different types of cells. The first step of processing data of RNA-seq is to exactly align millions of sequencing fragments against the reference sequence, which is called transcriptome sequence alignment. The existing sequence alignment tools for RNA-seq almost rely on canonical splice site signals, which, to some extent, limits the ability to identify novel splice sites. Therefore, we design a method independent from splice site signals, named RNAMap. It divides the sequencing fragments according to overlapping seeds method and scans the reference sequence via sliding windows with left and right anchors. In this way, splice sites within seeds can be identified. The computational experiments indicate that RNAMap not only reaches a precision of over 95%, but also outperforms the existing softwares in recall rate.

alternative splicing; RNA-seq; sliding windows; splice sites

國家自然科學(xué)基金(60533020)

2016-03-17;收到修改稿時間:2016-04-11

10.15888/j.cnki.csa.005443