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

張 勇, 徐 云

1(中國科學技術大學 計算機科學與技術學院, 合肥 230027)2(中國科學技術大學 安徽省高性能計算重點實驗室, 合肥 230027)3(國防科學技術大學 高性能計算協同創新中心, 長沙 410073)

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

張 勇1,2, 徐 云1,3

1(中國科學技術大學 計算機科學與技術學院, 合肥 230027)2(中國科學技術大學 安徽省高性能計算重點實驗室, 合肥 230027)3(國防科學技術大學 高性能計算協同創新中心, 長沙 410073)

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

選擇性剪接; 高通量轉錄組測序; 滑動窗口; 剪接位點

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

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

圖1 真核細胞基因結構圖

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

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

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

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

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

1 方法

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

1.1 初始比對階段

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

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

圖2 RNAMap流程圖

1.2 分段比對階段

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

1.2.1 劃分種子

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

1.2.2 種子定位

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

① 靜態表

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

② 動態表

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

圖3 靜態表示例

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

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

1.2.3 種子整合

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

2 實驗結果

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

2.1 模擬數據集測試

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

表1 各軟件的模擬數據集比對結果統計

由表1可知, 雖然CRAC和HISAT可以將絕大多數的測序片段(分別為97.46%和96.36%)比對到基因組上, 但是精確度比較低, 尤其是CRAC的精確度不足70%. TopHat2、SpliceMap、OLego、MapSplice、CRAC的召回率都比較低, 雖然HISAT的召回率達到80.17%,但是它的精確度也僅僅高于SpliceMap和CRAC. 因為我們的模擬數據集是隨機產生的, 所以其中有大量的reads含有非經典的剪接位點, 但上述軟件無法處理此類情況, 從而造成其比對的質量相對較低.

RNAMap的精確度高達96.15%, 雖然稍低于OLego的96.22%, 但是RNAMap的召回率達到了85.60%, 明顯高于其它幾個軟件. 這是因為在比對的整個過程中, RNAMap并未受經典剪接位點信號的限制, 因而可以更準確地將各類reads比對到參考序列上.

2.1 真實數據集測試

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

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

圖4 各軟件的真實數據集比對結果統計

3 結語

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

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

國家自然科學基金(60533020)

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

10.15888/j.cnki.csa.005443