丁勝楠，吳 鳴，徐 云

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

0 引言

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

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

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

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

1 相關(guān)工作

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

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

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

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

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

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

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

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

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

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

1.2.1定長(zhǎng)種子過濾方法

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

1.2.2變長(zhǎng)種子過濾方法

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

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

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

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

2.2 方法描述

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

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

圖2 滑動(dòng)窗口示意圖

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

表1 核苷酸成分表

2.2.2索引表上檢索符合篩選條件的滑動(dòng)窗口

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

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

(1)

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

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

(2)

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

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

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

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

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

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

3 實(shí)驗(yàn)與分析

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

3.1 過濾效果

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

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

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

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

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

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

表4 線蟲基因組上計(jì)算時(shí)間對(duì)比

表5 人類基因組計(jì)算時(shí)間對(duì)比

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

4 結(jié)論

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

參考文獻(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.