基于長種子的二代測序序列找全比對算法①

吳 邪,劉 歡,徐 云

1(中國科學技術大學 計算機科學與技術學院,合肥 230026)

2(安徽省高性能計算重點實驗室,合肥 230026)

第二代測序技術(NGS)[1]的出現使生物學家能夠以較低的成本快速確定DNA 分子的核苷酸排列.如今,NGS 平臺每天均會產生海量的長度為75 至300 個堿基對(bp)的測序片段(reads,也稱之為讀段)需要進行序列比對分析,這是下游生物學應用的基礎步驟.序列比對分析結果可以解釋人類種群之間的基因多樣性,對研究物種進化過程和引發遺傳病的基因變異有著深遠意義[2].

一般地,序列比對算法可以分為找最佳比對算法(best-mapper)和找全比對算法(all-mapper)兩類,前者用于找到每條讀段在參考基因組上匹配概率最高的一個或幾個位置,而后者則用于找到滿足誤差的全部匹配位置.這兩類算法應用場景不同,如全轉錄組測序分析以及DNA 和蛋白質之間的作用關系分析,利用找最佳比對算法找出讀段的一個或幾個最佳匹配位置已能達到要求.而對于結構變異探測以及拷貝數變異(CNV)分析這類需要找出所有基因組重復區域信息的應用,找全比對算法是必不可少的.由于找全比對算法需要枚舉所有可能的匹配位置,故而其計算成本更高,更迫切需要研究高性能的算法.目前主流的找全比對算法存在驗證階段時間開銷過大的問題,提高驗證階段的時間效率,成為了測序序列找全比對算法中的重要優化環節.

主流的測序序列找全比對算法普遍采用種子擴展(seed-and-extend)方法,該方法分為過濾和驗證兩個階段.在過濾階段,首先從讀段中選取多個長度固定的堿基序列片段作為種子(seeds),然后在參考基因組索引中檢索種子的出現位置并過濾出合適的候選位置; 在驗證階段,需要利用動態規劃算法驗證候選位置是否是讀段匹配參考基因組的實際位置.該方法過濾階段使用的參考基因組索引主要包括哈希索引和FM 索引[3]兩類.哈希索引的優勢是種子定位速度快,而問題是空間開銷大,像人類基因組的哈希索引超過11 GB,并且種子越長其空間開銷會指數級地增長.FM 索引優勢是空間開銷要小得多,人類基因組的索引僅需要大約3 GB,其缺陷是定位操作時間相對前者要大得多,并且種子越長時間開銷越大.受這些因素的限制,現有找全比對算法通常選取較短的種子來處理讀段,一般為11 至14 bp.由于短種子在參考基因組上的出現頻率高,導致過濾和驗證的位置數量過多,而驗證階段時間開銷占比更大,是提升比對算法時間性能的關鍵.

為此,本文提出了一種基于長種子的二代測序序列找全比對算法.該算法中的參考基因組索引是結合普通哈希索引與布隆過濾器(Bloom filter,BF)[4]構建的面向長種子的哈希索引.該索引方法通過模運算縮減哈?？臻g,避免了哈希索引在長種子情況下空間開銷過大,并應用布隆過濾器識別同一存儲位置上的不同種子.對于31 億堿基位點的人類基因組,本文的索引方法可以在約8 GB 空間上索引長度達30 bp 的種子,我們的找全比對算法相比于已有最好的找全比對算法時間減少26%.

1 相關研究

1.1 基因組索引技術

(1)經典索引方法

參考基因組索引是大部分序列比對算法中的基本數據結構,索引的空間開銷和檢索效率都是影響序列比對算法性能的重要因素.依據索引結構不同,可將基因組索引分為后綴數組、FM 索引和哈希索引3 類.

后綴數組的思路是對給定文本串S的所有后綴按照字典序進行排序,并將所有后綴在S中的起始位置存在數組SA中.對于一個待檢索的模式串P,所有以P為前綴的文本串在S上的位置就存儲在后綴數組的一個連續區間之中.雖然后綴數組支持定位非固定長度的種子,但是其空間開銷相比于哈希索引沒有優勢并且檢索速度要慢得多.因此,目前只有少數序列比對算法采用后綴數組來索引參考基因組,這類算法包括Vmatch[5]和Segemehl[6]等.

Ferragina 和Manzini 將后綴數組與BWT 變換(Burrows-Wheeler transform)結合,提出了FM 索引[3].該索引方法對后綴數組中的特定位置進行采樣,并利用BWT 變換來求解非采樣位置.FM 索引的定位速度相比于后綴數組進一步下降,但通過僅存儲采樣位置使所需存儲空間大大減少,是目前空間消耗最小的參考基因組索引.由于FM 索引定位種子速度較慢,目前常用于找最佳的序列比對算法,如BWA[7]、SOAP2[8]、Bowtie2[9]等BLAST[10]算法首次提出用哈希表構建基因組索引.哈希索引保存了參考基因組上每條q-gram(長度為q的連續子序列)的位置信息,將q-gram 作為哈希表的key,并把此q-gram 在基因組上出現的位置保存為value.由于哈希索引僅需一次哈希函數運算即可定位種子,非常適合需要頻繁定位種子的找全比對問題,因此大部分找全比對算法都采用哈希索引方法來構建參考基因組索引,包括mrFAST[11]、Hobbes2[12]、FEM[13]、BitMapper[14]等.哈希索引結構如圖1 所示.

圖1 哈希索引結構

(2)索引技術新發展

基因組索引技術的改進是提升序列比對算法性能的重要一環.目前針對FM 索引的改進主要集中在種子定位速度方面.Chacón 等人提出了n-step FM 索引[15],該索引方法利用二維的BWT 結構使得FM 索引中的向后搜索算法(backward search algorithm)能一次向后搜索n個字符,較大程度地提升了定位速度.Fernandes等人將采樣后的最長公共前綴數組應用到FM 索引方法中[16],提升了定位速度并節省了空間.Cheng 等人提出了FMtree 種子定位算法[17],將FM 索引定位操作的搜索空間組織成多叉樹,使多個采樣位置能被同時計算出,極大地提升了定位速度.

對于哈希索引,其種子定位操作的時間復雜度已經是常數級.因此,改進策略主要集中于減少空間消耗.稀疏化的哈希索引[13,18]是降低哈希索引空間開銷的有效方法.文獻[18]中按照特定規則采樣出部分q-gram,并將它們的位置保存在哈希表中,后續再計算出未采樣q-gram 的位置.該方法以較小的時間復雜度代價大大降低了哈希索引的空間消耗.文獻[13]同樣是基于采樣q-gram 的思路,與前者不同的是,該方法在選取種子時采用分組選種策略盡量選擇已被采樣的qgram 作為種子,避免了計算未采樣q-gram 的位置的時間消耗.

1.2 主流的過濾階段方法

早期的找全比對算法,如BLAST 算法,采用連續種子方法,將讀段的每一個q-gram 當作種子在索引中檢索,之后擴展并合并這些種子的位置.這種方法需要定位大量種子,時間開銷大,尤其對于大型基因組,其種子出現頻率較高,導致需要驗證的候選位置過多,極大增加了后續驗證階段的計算量.另外,連續種子方法也無法正確匹配包含插入刪除錯誤的讀段.

MAQ[19]和SOAP[20]等方法采用了填充種子方法(spaced seed),其利用多個填充種子模板來覆蓋一條讀段.每個模板由特定長度的0 和1 組成,忽略掉0 對應的堿基字符.這類方法一般會設計多個0 位位置不同的種子模板,可以保證當讀段有少量置換錯誤時至少有一個模板可以命中.填充種子方法時間開銷比連續種子方法小,但同樣無法正確匹配包含插入刪除錯誤的讀段.

目前性能最好的幾種找全比對算法普遍采用基于鴿巢原理的過濾方法.這類方法不僅降低了過濾階段的時間開銷,還能正確匹配包含插入刪除錯誤的讀段.mrFAST 是采用這類方法的經典算法之一,其假定每條讀段最多允許e個誤配,并在讀段上不重疊地劃分出e+2 個種子,由鴿巢原理可知,至少存在2 個種子不包含誤配.這表示對于一個讀段與參考基因組匹配的位置,該位置至少是兩個或兩個以上的種子的出現位置.mrFAST 將只出現一次的種子位置過濾掉,從而減少候選位置數量.Hobbes2 的過濾階段方法目前效果最好,其在mrFAST 的基礎上采用動態規劃算法選取種子出現頻次和最小的組合,進一步減少了候選位置數量.

2 問題定義及相關數據結構

2.1 找全序列比對問題定義

測序序列找全比對問題本質上是求解海量短字符串在長文本串中的全部匹配位置.從數學的角度上可以定義為四元組MS A=(Σ,S,A,F).對于DNA 序列,Σ={A,C,G,T}分別表示4 種堿基不同的脫氧核糖核酸;S={S1,S2,···,Sn}表示測序平臺產生的讀段集合,其中,Si=(bi1,bi2,···,biLi),bik∈Σ,Li是讀段Si的長度;A=(Aik)N×M,其中(M≥max(L1,L2,···,LN)),Aik表示讀段Si的第k個堿基經過驗證后的比對結果,F是打分函數,F(A)表示讀段最終的比對結果.找全比對的目的是找到F(A)分值滿足一定閾值的全部位置.

2.2 布隆過濾器

布隆過濾器能夠查詢一個數據是否包含于此過濾器對應的數據集合,它由k個哈希函數和1 個二進制位向量組成,位向量的初始值全為0.插入新數據時,k個哈希函數會將該數據映射到位向量上的k個不同位置,并將相應位置的值置1.查詢時,用相同的k個哈希函數進行映射,檢查映射的位置上的值是否全為1 即可確定集合中是否包含該元素.布隆過濾器的結構如圖2 所示.

圖2 布隆過濾器結構

布隆過濾器是一種時間和空間效率都比較高的數據結構,其時間復雜度僅與哈希函數個數有關,空間復雜度為 ,其中m是位向量的長度.當兩個數據元素的k個哈希函數映射的位置均相同時,布隆過濾器會產生假陽性結果,將原本不存在于集合中的數據報告為存在.為降低假陽性率,在實際使用中需權衡k、m和數據集合大小n之間的關系,一般先考慮好m的大小,然后即可求得假陽性概率最小時的哈希函數個數k,根據已有研究工作[21],其計算式可表示為:

假陽性概率p與k、m、n相關,根據文獻[21],其計算式如式(2):

3 比對算法設計

本文首先利用Jellyfish 軟件[22]統計種子的長度與出現頻率的關系,進而確定合適的種子長度.隨后設計并實現了面向長種子的哈希索引,該索引結構的優勢在于支持長種子的同時定位速度快,能夠得到較少的候選位置.然后固定間隔選取q-gram 作為種子,并在長種子哈希索引中定位.最后利用目前效率最高的位并行算法來驗證候選位置.該算法流程可分為4 個步驟: (1)構建索引; (2)選取種子; (3)定位及過濾; (3)驗證候選位置.

3.1 構建索引

在構建長種子索引之前,首先需要確定種子的長度.由Jellyfish 軟件得到的統計結果可知,對于包含31 億個字符的人類基因組hg19,30 bp 長度的種子中出現頻率在10 次以內的占比達90%,如圖3 所示.

圖3 種子出現頻率與其長度變化關系曲線

由圖3 可知,種子長度達到30 bp 時,出現頻率小于10 的種子占比趨于穩定.另一方面,二代測序平臺產生的讀段長度較短,過長的種子容易包含測序錯誤.因此,本文將種子長度設定為30 bp.

長種子哈希索引主要包括一張哈希表以及多個布隆過濾器.其中,每個布隆過濾器對應參考基因組的一個區域.采用這種結構的根本目的是縮減哈希索引存儲長種子的空間開銷.以種子長度取30 bp 為例,普通哈希索引中需要存儲430個不同key 的信息.而在實現索引結構時每個key 對應一個4 B 大小的指針,僅指針信息的空間開銷就達到430×4 B=222TB.這樣的存儲代價過大,顯然無法投入實際應用.為此,本文通過模運算將哈?？臻g限制為固定值M,極大減小了長種子情況下哈希索引的空間開銷,以M取225為例(實際實現時,為減少哈希碰撞,M應取質數),存儲全部指針信息僅需225×4 B=128 MB.索引構建過程如圖4 所示,首先用長度為30 的滑動窗口在參考基因組上逐位選取q-gram.將每個q-gram 轉換成64 位無符號整型數值后,再模M轉換為對應的key 值,然后將q-gram 在基因組上的位置記錄在哈希表中該key 對應的表項中.

圖4 索引構建過程

利用模運算限制哈?？臻g大小的做法會引起新的問題: 若兩個不同的q-gram 模M后的key 值相同,則它們在基因組上的位置會記錄在哈希表中同一個key 對應的value 中,導致定位種子時產生大量錯誤的候選位置.本文利用布隆過濾器來篩選掉錯誤候選位置,首先將整個參考基因組均勻劃分為多個不重疊區域,為每個區域建立一個對應的布隆過濾器,具體過程如圖5 所示.

圖5 索引構建過程

每個布隆過濾器存儲的元素為其對應的基因組區域內所有q-gram.其作用是驗證某個q-gram 是否屬于該布隆過濾器對應的基因組區域.對于某個待定位種子,檢索哈希表得到其位置集合后,需利用布隆過濾器驗證每個位置所在的基因組區域是否包含該種子.這樣可判斷此位置是否為待定位種子的實際出現位置,從而篩選掉錯誤的候選位置.

3.2 選取種子

由于二代測序產生的讀段長度短,基于鴿巢原理的種子選取方法不再適用于基于長種子的算法.本文采用滑動窗口方法: 設定步長l-step,并用長度30 的滑動窗口在讀段上每次右移l-step 位選取種子,具體過程如圖6 所示.

圖6 種子的生成

采用滑動窗口方法是為了盡可能選出不包含測序錯誤的種子,能有效避免定位得到的位置集合中不包含讀段與基因組實際匹配位置的情況.參數l-step 會影響比對結果,過大會使部分讀段無法匹配到參考基因組上或漏掉部分匹配位置; 過小會提升比對精度,但需要定位的種子數量會成倍增多,導致定位過程的時間開銷增大.本文經過多次實驗比較后發現,將l-step 設定為10 時,我們的比對算法能達到與現有找全比對算法同等的精度.

3.3 定位及過濾

將選取的種子轉換為key,在哈希表中檢索.對于檢索出的每個位置,利用該位置所在基因組區域的布隆過濾器驗證此位置是否為對應種子實際的出現位置,驗證通過則保留下來作為候選位置.種子定位及過濾過程如圖7 所示.

圖7 種子定位及過濾

3.4 驗證候選位置

候選位置是讀段與參考基因組最有可能匹配上的位置,但其中仍包含部分非實際匹配的位置,需要在驗證階段將這些錯誤位置篩選掉.本文采用“帶狀”Myers算法[23]的一種高效版本[14]來校驗每個候選位置,此版本的算法通過位并行提高了計算速度,并利用CPU 上的128 位寄存器和SSE 指令集來加速驗證過程,實現了同一時間驗證多個候選位置的操作,極大降低了驗證階段的時間開銷.

4 實驗分析

本文所有實驗均在同一機器的相同配置環境下運行.實驗平臺為: Intel Xeon Gold 5120 CPU (14 核心28 線程,2.20 GHz),503 GB DDR4 RAM,Ubuntu 18.04 64 位操作系統.

實驗所用數據均下載自NCBI 數據庫.其中參考基因組為千人基因組項目中的Genome Reference Consortium GRCh37 參考序列,版本號為 hg19.測序數據集說明如下:

(1)R1: 千人基因組項目中第NA17871 號測序序列SRR007347,包含30 萬條長度100 bp 的讀段.

(2)R2: 千人基因組項目中第NA11840 號測序序列ERR012100,包含10 萬條長度100 bp 的讀段.

4.1 候選位置數量對比實驗

為驗證本文所采用的長種子方法能得到更少量的候選位置,選擇mrFAST 和Hobbes2 來進行讀段平均候選位置數量的對比實驗,所用數據集為R1.其中Hobbes2 的過濾方法是目前效果最好的方法.實驗結果如表1 所示.

表1 各算法的讀段平均候選位置數量結果

由表1 可知,本文算法得到的讀段平均候選位置數比mrFAST 少約87%,比Hobbes2 少約69%.采用長種子定位的方法在候選位置數量上有明顯的優勢.

4.2 算法性能對比實驗

在找全比對算法的整體流程中,驗證階段的時間占比最大.以目前速度最快的BitMapper 為例,其驗證階段的時間占比在70%以上.因此,相比于定位長種子的額外時間消耗,本文方法在驗證階段帶來的時間性能提升更為顯著.為驗證此觀點,選擇目前速度最快的幾種找全比對算法與本文算法做對比,實驗結果如表2 所示.

表2 各算法性能測試結果

本文利用匹配上的讀段比例和召回率來衡量找全比對算法的精度,其中,匹配上的讀段比例表示成功匹配到參考基因組上的讀段占比; 召回率表示讀段的正確匹配位置被算法報告出來的比例.由表2 可知,對于匹配上的讀段比例,mrFAST 的結果最好,其他4 個算法相差不大.對于召回率,FEM 和BitMapper 的結果較好,本文算法則略優于Hobbes2 和mrFAST.在內存開銷上,mrFAST 表現最好,FEM 次之,本文算法優于Hobbes2 和BitMapper.在時間性能方面,本文算法比目前速度最快的BitMapper 快了約26%.

5 結論與展望

本文提出一種二代測序序列找全比對算法,在算法中采用哈希索引與布隆過濾器結合的方法實現了長種子的定位.利用長種子的低頻特性可得到更少量的候選位置,從而提升算法的運行速度.實驗結果表明,本文算法相比于其他主流的找全比對算法,候選位置數量更少,在維持相同水平精度的同時速度更快.由于定位長種子帶來的額外時間開銷,我們的算法在過濾階段比其他算法需要更長的運行時間,如何優化布隆過濾器篩選位置的策略從而提高長種子的定位速度將是下一步的工作重點.