吳 邪,劉 歡,徐 云

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

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

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

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

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

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

1 相關(guān)研究

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

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

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

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

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

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

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

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

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

1.2 主流的過(guò)濾階段方法

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

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

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

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

2.1 找全序列比對(duì)問(wèn)題定義

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

2.2 布隆過(guò)濾器

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

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

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

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

3 比對(duì)算法設(shè)計(jì)

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

3.1 構(gòu)建索引

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

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

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

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

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

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

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

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

3.2 選取種子

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

圖6 種子的生成

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

3.3 定位及過(guò)濾

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

圖7 種子定位及過(guò)濾

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

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

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

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

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

(1)R1: 千人基因組項(xiàng)目中第NA17871 號(hào)測(cè)序序列SRR007347,包含30 萬(wàn)條長(zhǎng)度100 bp 的讀段.

(2)R2: 千人基因組項(xiàng)目中第NA11840 號(hào)測(cè)序序列ERR012100,包含10 萬(wàn)條長(zhǎng)度100 bp 的讀段.

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

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

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

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

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

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

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

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

5 結(jié)論與展望

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