張 勇, 徐 云

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

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

張 勇1,2, 徐 云1,3

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

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

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

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

定性和定量研究轉(zhuǎn)錄組的傳統(tǒng)方法是構(gòu)建cDNA或表達(dá)序列標(biāo)簽(EST)文庫(kù), 然后通過(guò)Sanger測(cè)序進(jìn)行后續(xù)分析. 但是, 因?yàn)镾anger測(cè)序技術(shù)的成本較高且通量較低, 所以這種方法十分昂貴和低效. 隨著下一代測(cè)序技術(shù)(next-generation sequencing, NGS)的迅猛發(fā)展, 高通量轉(zhuǎn)錄組測(cè)序(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)該能夠評(píng)估每一種剪接異構(gòu)體的表達(dá)水平. 然而, 所有分析過(guò)程的第一步都是要將RNA-seq中得到的測(cè)序片段(reads)比對(duì)到基因組上的原始位置, 而這些短片段的長(zhǎng)度從數(shù)十堿基到數(shù)百堿基不等, 數(shù)量有幾十萬(wàn)甚至幾百萬(wàn)和上千萬(wàn), 所以, 比對(duì)的過(guò)程是極其耗時(shí)的.

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

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

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

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

1 方法

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

1.1 初始比對(duì)階段

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

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

圖2 RNAMap流程圖

1.2 分段比對(duì)階段

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

1.2.1 劃分種子

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

1.2.2 種子定位

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

① 靜態(tài)表

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

② 動(dòng)態(tài)表

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

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

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

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

1.2.3 種子整合

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

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

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

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

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

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

由表1可知, 雖然CRAC和HISAT可以將絕大多數(shù)的測(cè)序片段(分別為97.46%和96.36%)比對(duì)到基因組上, 但是精確度比較低, 尤其是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), 但上述軟件無(wú)法處理此類情況, 從而造成其比對(duì)的質(zhì)量相對(duì)較低.

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

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

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

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

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

3 結(jié)語(yǔ)

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

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

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

2016-03-17;收到修改稿時(shí)間:2016-04-11

10.15888/j.cnki.csa.005443