吳 邪,劉 歡,徐 云

1(中國科學(xué)技術(shù)大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院,合肥 230026)

2(安徽省高性能計算重點(diǎn)實(shí)驗(yàn)室,合肥 230026)

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

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

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

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

1 相關(guān)研究

1.1 基因組索引技術(shù)

(1)經(jīng)典索引方法

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

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

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

圖1 哈希索引結(jié)構(gòu)

(2)索引技術(shù)新發(fā)展

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

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

1.2 主流的過濾階段方法

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

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

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

2 問題定義及相關(guān)數(shù)據(jù)結(jié)構(gòu)

2.1 找全序列比對問題定義

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

2.2 布隆過濾器

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

圖2 布隆過濾器結(jié)構(gòu)

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

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

3 比對算法設(shè)計

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

3.1 構(gòu)建索引

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

圖3 種子出現(xiàn)頻率與其長度變化關(guān)系曲線

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

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

圖4 索引構(gòu)建過程

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

圖5 索引構(gòu)建過程

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

3.2 選取種子

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

圖6 種子的生成

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

3.3 定位及過濾

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

圖7 種子定位及過濾

3.4 驗(yàn)證候選位置

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

4 實(shí)驗(yàn)分析

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

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

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

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

4.1 候選位置數(shù)量對比實(shí)驗(yàn)

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

表1 各算法的讀段平均候選位置數(shù)量結(jié)果

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

4.2 算法性能對比實(shí)驗(yàn)

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

表2 各算法性能測試結(jié)果

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

5 結(jié)論與展望

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