丁勝楠，吳 鳴，徐 云

(1.中國科學(xué)技術(shù)大學(xué) 軟件學(xué)院，安徽 合肥 230027；2.中國科學(xué)技術(shù)大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院，安徽 合肥 230027；3.安徽省高性能計算重點實驗室，安徽 合肥 230027)

0 引言

DNA測序已經(jīng)成為生物醫(yī)學(xué)研究中不可或缺的工具，可直接應(yīng)用于個性化醫(yī)療和基因變異導(dǎo)致的疾病研究。目前，二代測序技術(shù)，比較常見的有：羅氏454測序、Illumina等。而現(xiàn)在最常用的NGS技術(shù)就是Illumina測序技術(shù)，它可以確保在幾十個小時內(nèi)產(chǎn)生幾百GB甚至上TB的測序數(shù)據(jù)，數(shù)據(jù)量特別巨大，完全能夠滿足高通量測序的通量要求。隨著二代測序平臺的不斷發(fā)展，較第一代測序技術(shù)測序通量迅速增加，如目前常用的Illumina測序技術(shù)，它能在幾個小時內(nèi)產(chǎn)生數(shù)GB長度為100 bp(位點)的片段測序數(shù)據(jù)，將這些大量片段測序數(shù)據(jù)reads比對到參考基因組上，成為了測序序列比對算法面臨的嚴(yán)峻挑戰(zhàn)。

測序平臺測出的短DNA序列片段(又稱為reads)數(shù)量極高，達(dá)到數(shù)十億個。在序列比對過程中，每個read在基因組中可能有一個或多個相似的片段，因此均需要在參考基因組中進(jìn)行比對。同時，由于測序平臺的局限性，每個測得的read可能會出現(xiàn)一些錯誤，以及基因組產(chǎn)生變異，如DNA序列中堿基可能會發(fā)生相應(yīng)的插入、刪除和替換(通常來說，read的這一類變化率小于等于5%)。因此，允許一定數(shù)量誤配的字符串匹配是一個比較困難的問題，并會耗費(fèi)大量的時間。以上問題導(dǎo)致基因序列比對在基因數(shù)據(jù)分析中占據(jù)了大量的時間。

測序序列比對分為找最好和找全比對兩種應(yīng)用問題。對于給定的一個編輯距離(插入、刪除、替換)或海明距離的閾值k，找出每條測序序列(read，讀數(shù))在單個或多個基因組上滿足閾值k的位置。如果需要找出所有位置，這類問題稱為找全問題，對應(yīng)的比對算法為找全算法(all-mapper)；如果只需要找出分值(編輯距離或海明距離與位點質(zhì)量的綜合函數(shù))最高或較高的一些位置，這類問題稱為找最好問題，對應(yīng)的比對算法為找最好算法(best-mapper)。找全測序序列比對是其中的重要問題，已有的mrFAST[1]、RazerS3[2]、GEM[3]、Bitmapper[4]、Bowtie2[5]、Masai[6]、Hobbes[7]等就是找全比對算法的代表。這些找全測序序列比對算法大多采用稱為“種子-擴(kuò)展(seed-and-extend)”的方法。該方法是先進(jìn)行種子過濾，即將reads劃分為若干個不重疊的長度為11 bp左右的種子，應(yīng)用鴿巢原理確定出候選位置，對于人類基因組這樣的候選位置往往會突現(xiàn)數(shù)千之多，嚴(yán)重增加了測序序列比對算法后續(xù)驗證階段的處理時間。因此，對于找全問題的測序序列比對算法，研究一種有效的過濾方法來減少候選位置來實現(xiàn)比對算法效率的提升，具有重要的研究意義和應(yīng)用價值。

傳統(tǒng)的種子過濾方法可以稱為細(xì)粒度過濾法，這種方法產(chǎn)生的候選位置較多，有必要研究一種粗粒度過濾方法(比如區(qū)域過濾)來進(jìn)一步地減少候選位置。文獻(xiàn)[7]中提到的過濾算法，雖然上下邊界限制條件比較好，但過濾是對每個候選位置的操作，而驗證是對最終結(jié)果的操作，所以過濾需要盡可能快。這里采用了計數(shù)索引來加速過濾操作。Hobbes也采用了類似的方法，但其過濾未采用索引技術(shù)，增加了過濾的總時間。而本文的方法在較好的過濾效果和高效的計算中取得了一個平衡，綜合來說，本文的區(qū)域過濾算法時間和性能較好。

1 相關(guān)工作

1.1 種子擴(kuò)展方法(seed-and-extend)

種子擴(kuò)展方法的第一步是將read分割為幾個相對較短的，非重疊的seed。在查找過程中即使read中包含錯誤，但只要存在一個seed正確，則read的正確位置仍然可以被找到。而無錯誤的種子可以通過鴿巢原理進(jìn)行選取，對于一個read，若允許最大誤配數(shù)為k，那么由鴿巢原理得知，只需將其分為k+1個無重疊的seed，必然存在一個seed為正確的。對于一個不包含錯誤的seed來說，通過索引查找其在基因組中匹配的位置，也就得到了read的位置。

第二步，借助索引對seed進(jìn)行快速地查找。在序列比對的過程中，seed的位置通常存儲在對基因組預(yù)處理過的索引中(常用的索引有Hash索引和基于Burrows-Wheeler變換(BWT[8])的FM-index[9])，用來加快查找的速度。

最后一步，將查找到的位置返回基因組中做擴(kuò)展驗證。目前常用驗證方法是動態(tài)規(guī)劃算法(如Smith-Waterman[10]等)。通過擴(kuò)展驗證，檢查read在基因組中匹配的質(zhì)量，得出最終的位置，可以用于后續(xù)的生物學(xué)意義分析。

種子擴(kuò)展方法流程如圖1所示。

圖1 種子擴(kuò)展方法流程

從圖1可以看出，將read切分成若干seed后，會產(chǎn)生很多重復(fù)以及錯誤的候選位置，而這些重復(fù)或者錯誤的候選位置也是需要過濾掉的種子。而參考基因組(待比對基因組)某一片段上堿基的個數(shù)是固定已知的，需要有一種方法將read中對應(yīng)堿基含有誤配的上下界找出來(且保證所找到的候選位置是無遺漏的)，并以此做限定條件，去匹配含有編輯距離長度的基因片段里堿基的個數(shù)，進(jìn)而篩選掉不在這個范圍內(nèi)的候選位置，來縮短比對和驗證時間。

對于一個給定的read，假設(shè)允許的誤配數(shù)為k，由鴿巢原理可知，將其分為k+1個seed，之后利用索引查找seed在基因組中的位置，鴿巢過濾后，再通過動態(tài)規(guī)劃的驗證擴(kuò)展在基因組中進(jìn)行精確比對，求出其最終匹配位置。

1.2 現(xiàn)有各種種子過濾方法

目前各種種子選取方法按種子長度可分為定長和變長，按種子數(shù)量可分為k+1、k+2等，其中k為編輯距離。不同的過濾方法往往帶來非常不同的過濾效果，以下介紹近些年來種子過濾方法的研究和進(jìn)展。

1.2.1定長種子過濾方法

FastHash提出，對于一個長度為100 bp的read，若seed的長度為12，則最多可選取8個種子。當(dāng)k=5時，由鴿巢原理可知選取6個種子至少有1個種子一定是準(zhǔn)確匹配的。先選取8個種子，先去除頻次最大的2個種子，這樣剩余的6個種子的候選位置就大大減少。由于種子的位置不同，可以組合出各種不同的k+1種子選取方案，Hobbes提出了動態(tài)規(guī)劃選取k+1個定長種子，使其各種子的頻次和為最小。Hobbes2在Hobbes的基礎(chǔ)上，提出動態(tài)規(guī)劃選取k+2頻次最小的定長種子。由鴿巢原理易知，至少2個種子的位置是正確的，且這兩個種子應(yīng)該保持在reads上的準(zhǔn)確相對位置，因此可對于k+2的種子進(jìn)行相對位置再過濾，這樣大大減少了候選位置。

1.2.2變長種子過濾方法

基于定長種子的過濾算法盡管時間復(fù)雜度低，但過濾后得到的位置數(shù)仍很多，而基于變長種子的過濾算法雖然時間消耗較定長過濾大，但仍在可接受范圍內(nèi)，而且可以找到最佳候選位置?；谧冮L種子的過濾算法如下：OSS提出了通過動態(tài)規(guī)劃，選取k+1個頻次最小的變長種子進(jìn)行擴(kuò)展驗證；GEM提出了一種基于FM-index的變長種子的選取方法，先設(shè)定一個閾值t，從read的初始位置開始選取種子，若種子頻率高于t，則擴(kuò)展種子長度，否則，進(jìn)行下一個種子的查找，直至選取k+1個種子。相比較變長種子過濾而言，由于定長種子速度快，種子位置多，對應(yīng)的過濾效果更顯著，因此本文基于定長種子過濾方法進(jìn)行實驗。

2 區(qū)域過濾方法設(shè)計

2.1 區(qū)域過濾方法的思想

為確保read在參考基因組的所有候選位置不漏選，制定了一個上下界來過濾，并對該上下界進(jìn)行正確性證明，同時對該上下界進(jìn)行優(yōu)化選取。在此基礎(chǔ)上設(shè)計出一種帶有間隔且與read長度相當(dāng)?shù)幕瑒哟翱?，建立一種有效的滑動窗口方法。由于每個read包含的每種核苷酸(A/G/C/T)個數(shù)是固定的，經(jīng)過k編輯距離后的read中每種核苷酸個數(shù)至多變化為k，最終目標(biāo)是在這些滑動窗口中找出滿足核苷酸個數(shù)變化的所有窗口。

2.2 方法描述

2.2.1構(gòu)建滑動窗口核苷酸組成索引表

以線蟲基因組為例，將其分成1 000份，每份有100 000長度的區(qū)段，設(shè)計方案中就以這100 000長度的區(qū)段為基因組作為示例，并將若干條長度為100 bp的read以編輯距離k=5比對到該區(qū)段上?；瑒哟翱谠O(shè)計如圖2所示。為方便說明，以長度為110的滑動窗口為例。選取長度為110的滑動窗口，以間隔k=5進(jìn)行滑動，這樣在該區(qū)段上可以產(chǎn)生19 979個滑動窗口。

圖2 滑動窗口示意圖

按基因組建立如表1所示的A成分索引表，用于快速定位符合要求的滑動窗口。該索引表由每個滑動窗口的各個核苷酸成分個數(shù)組成。該索引表第1列為A成分個數(shù)，第2列為符合A成分個數(shù)的滑動窗口起始位置。

表1 核苷酸成分表

2.2.2索引表上檢索符合篩選條件的滑動窗口

現(xiàn)假設(shè)ρr(x)為x在readr中出現(xiàn)的個數(shù)，x表示核苷酸A/C/G/T。那么在k編輯距離的情況下，可以得出如下過濾公式：

ρr(x)-3k≤ρw(x)≤ρr(x)+3k，x∈{A,C,G,T}

(1)

公式(1)是一個過濾公式組，每個核苷酸都應(yīng)該將該公式作為一個相對獨(dú)立的情形進(jìn)行約束。

公式(1)是比較直觀的，但過濾條件比較寬松?，F(xiàn)在再考慮另一個相對緊湊些的過濾公式(2)，過濾公式(2)是將4個核苷酸進(jìn)行整體考慮。由于總堿基數(shù)目是固定的，每個堿基修改、刪除、插入都有相對于另一個堿基的增減。比如：將A替換為C，這就意味著A減少1個而C增加1個。這樣對于k次編輯距離，就得出過濾公式(2)：

(2)

公式(2)的邊界對公式(1)做了進(jìn)一步的限定，在保證正確性的前提下對候選位置從整體進(jìn)行了限定，使得過濾得到的區(qū)域更窄，進(jìn)一步提升過濾率。

綜上，正確候選位置需要滿足上述條件，從而可以過濾錯誤的候選位置。

2.2.3與k+1過濾方法結(jié)合獲取候選種子位置

將符合要求的區(qū)域(窗口)與seed位置進(jìn)行重合檢查，將在符合要求的區(qū)域(窗口)之外的seed位置過濾掉，從而減少種子候選集的大小。

2.2.4將區(qū)域過濾和BitMapper算法融合

將粗粒度和細(xì)粒度過濾后得到的種子候選位置應(yīng)用到BitMapper算法中，去檢測對應(yīng)融合后的BitMapper算法時間性能的提升。具體做法為，對BitMapper比對算法中的過濾部分代碼進(jìn)行修改，在線蟲基因組和人類基因組數(shù)據(jù)上進(jìn)行計算實驗，以檢驗算法的時間效率和精確度。

3 實驗與分析

在線蟲基因組和人類基因組上進(jìn)行計算實驗，檢測區(qū)域過濾方法的過濾效果及其計算時間的提升。線蟲基因組長度為1×108bp，相對人類基因組長度3×108bp而言規(guī)模要小些。兩個不同規(guī)模的基因組能較為全面地反映區(qū)域過濾方法的效果。

3.1 過濾效果

過濾效果包括窗口和種子的過濾效果，通過過濾來減少在驗證和比對過程所消耗的時間，作為性能是否提升的一個判別準(zhǔn)則。在線蟲基因組和人類基因組上的過濾效果分別如表2和表3所示，這里過濾后窗口(種子)比例對應(yīng)的公式為：過濾后窗口(種子)數(shù)/過濾前窗口(種子)數(shù)。

表2 線蟲基因組過濾結(jié)果

表3 人類基因組過濾結(jié)果

從表2可以看出，區(qū)域過濾方法在線蟲基因組上是有效的，并且種子過濾效果要比窗口過濾的效果略好。從表3可以看出，區(qū)域過濾方法在人類基因組上也同樣有效，相比于線蟲基因組上效果要差些，主要原因是線蟲基因組的結(jié)構(gòu)性比人類基因組的結(jié)構(gòu)性強(qiáng)些。

3.2 基于區(qū)域過濾的BitMapper算法時間性能

將區(qū)域過濾方法融合到原BitMapper算法中，對原BitMapper算法和融入?yún)^(qū)域過濾的BitMapper算法進(jìn)行計算時間的比較，比較結(jié)果如表4和表5所示。

表4 線蟲基因組上計算時間對比

表5 人類基因組計算時間對比

盡管區(qū)域過濾方法比單一k+1過濾方法要消耗稍多些的過濾時間，但候選種子位置的顯著減少，會大幅減少比對算法的驗證時間，總體使得比對算法時間明顯減少。表4表明線蟲基因組計算時間提升了36%左右，表5表明人類基因組計算時間提升了30%左右。

4 結(jié)論

本文針對二代測序序列比對算法，提出了一種基于區(qū)域的粗粒度過濾方法。將該方法應(yīng)用于線蟲和人類基因組上，使得找全序列比對BitMapper算法有明顯的時間性能提升。下一步的工作將考慮如何設(shè)計更精巧的索引,從而加快過濾速度，進(jìn)一步提升比對算法的時間性能。

參考文獻(xiàn)

[1] MARSCHALL T, MARZ M, ABEEL T, et al. Computational pan-genomics: status, promises and challenges[J]. bioRxiv, 2016: 043430.

[2] 1000 Genomes Project Consortium. An integrated map of genetic variation from 1,092 human genomes[J]. Nature, 2012, 491(7422): 56-65.

[3] BALL M P, THAKURIA J V, ZARANEK A W, et al. A public resource facilitating clinical use of genomes[J]. Proceedings of the National Academy of Sciences, 2012, 109(30): 11920-11927.

[4] LEVY S, SUTTON G, NG P C, et al. The diploid genome sequence of an individual human[J]. PLoS Biol, 2007, 5(10): e254.

[5] DANECEK P, AUTON A, ABECASIS G, et al. The variant call format and VCFtools[J]. Bioinformatics, 2011, 27(15): 2156-2158.

[6] MARCO-SOLA S, SAMMETH M, GUIGR, et al. The GEM mapper: fast, accurate and versatile alignment by filtration[J]. Nature Methods, 2012, 9(12): 1185-1188.

[7] LI H, DURBIN R. Fast and accurate short read alignment with Burrows-Wheeler transform[J]. Bioinformatics, 2009, 25(14): 1754-1760.

[8] RAHN R, WEESE D, REINERT K. Journaled string tree-a scalable data structure for analyzing thousands of similar genomes on your laptop[J]. Bioinformatics, 2014, 30(24): 3499-3505.

[9] HUANG L, POPIC V, BATZOGLOU S. Short read alignment with populations of genomes[J]. Bioinformatics, 2013, 29(13): i361-i370.

[10] DANEK A, DEOROWICZ S, GRABOWSKI S. Indexes of large genome collections on a PC[J]. PloS One, 2014, 9(10): e109384.