基于位圖索引和B+樹的BLAST改進(jìn)算法

1.中山大學(xué) 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，廣州 510275

2.北京師范大學(xué)珠海分校 信息技術(shù)學(xué)院，廣東 珠海 519085

1.中山大學(xué) 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，廣州 510275

2.北京師范大學(xué)珠海分校 信息技術(shù)學(xué)院，廣東 珠海 519085

1 引言

BLAST是Basic Local Alignment Search Tool的簡寫。是采用啟發(fā)式算法的一個(gè)生物序列比對的軟件工具包，最初由Altschul S在1990年發(fā)展起來[1-2]。由于其在序列比對之前進(jìn)行一次定位，并對定位片段使用得分矩陣進(jìn)行打分的方式找出高分片段HSP（High Score Pair），因此雖然BLAST軟件采用的算法喪失了一些敏感度，可是它也達(dá)到了計(jì)算速度的極大提高，從而BLAST逐漸成為生物信息學(xué)領(lǐng)域的最重要的一種軟件，使用BLAST算法進(jìn)行序列搜索與比對也成為了分子生物學(xué)的必備。

雖然BLAST算法速度非常快，可是面對爆炸式增長的生物數(shù)據(jù)庫就顯得力不從心。現(xiàn)代生物數(shù)據(jù)庫大小增長速度十分快速，NCBI的GenBank核酸堿基數(shù)目大概每14個(gè)月就翻一倍。而即使依照摩爾定律，計(jì)算機(jī)的性能也要18個(gè)月才能更新一倍。因此，如何能夠更好地使用序列比對算法將是一個(gè)重要的問題。

為此，本文設(shè)想使用位圖索引組織BLAST程序使用的生物數(shù)據(jù)庫，并依靠B+樹進(jìn)行二次索引，使得BLAST在進(jìn)行序列搜索的時(shí)候，能夠更快地找到HSP。

2 問題描述

BLAST算法的過程[3]分為3個(gè)步驟。分別是播種（Seeding），擴(kuò)展（Extension），以及評估（Evaluation）。

BLAST算法假定序列比對中查詢序列和數(shù)據(jù)庫序列中的重要序列有共同之處，因此可以依靠兩者之間的共同點(diǎn)找出應(yīng)該搜索的區(qū)域。根據(jù)這個(gè)條件，BLAST程序從查詢序列中提取出所有長度為一固定值W的子串，并定義這些長度固定的子串為單詞（Word），然后BLAST程序利用這些單詞建立一個(gè)單詞表，再根據(jù)這張表，對數(shù)據(jù)庫序列進(jìn)行一次遍歷，找出所有與單詞表中單詞相匹配的片段。如DNA序列S=GCACCTACTA，當(dāng)W=8時(shí)序列S可以生成三個(gè)單詞：GCACCTAC，CACCTACT，ACCTACTA。BLAST程序首先找出搜索數(shù)據(jù)庫中有哪個(gè)位置有這三個(gè)單詞，然后標(biāo)記這些單詞。這些搜索數(shù)據(jù)庫中找到的共同單詞被稱為命中（Word Hits），如圖1。之后BLAST只考慮在這些命中區(qū)域附近搜索。這些命中就像種子一般，這就是BLAST算法的第一個(gè)步驟：播種。后面的計(jì)算都從這里發(fā)展，依靠這樣BLAST能夠節(jié)省大量搜索不必要區(qū)域的時(shí)間。

圖1 命中

對于蛋白序列，命中的范圍除了包括準(zhǔn)確的匹配，還包括相似的單詞。因?yàn)橹豢紤]準(zhǔn)確的匹配可能會(huì)忽略一些重要的匹配。例如搜索數(shù)據(jù)庫中有段序列T=MGKGD，明顯MGR，GRG和RGD都和序列T不匹配。可是在生物學(xué)中兩個(gè)序列是十分相似的。為了避免出現(xiàn)這種情況。BLAST定義單詞除了包括準(zhǔn)確的匹配之外，還包括依據(jù)得分矩陣生成的鄰域（Neighborhood）。每個(gè)鄰域都有自己的得分，BLAST在搜索前會(huì)設(shè)置一個(gè)值T，并搜索時(shí)自動(dòng)生成所有得分大于T的單詞。由于蛋白序列搜索與核酸序列搜索有很大相似性，以下主要討論DNA序列搜索的BLAST算法。

在實(shí)際的操作過程中，BLAST程序會(huì)先利用查詢序列根據(jù)用戶設(shè)定的W值和T值生成一個(gè)單詞表。然后順序查詢數(shù)據(jù)庫序列，對于每一個(gè)長度為W的子序列，BLAST都在單詞表中進(jìn)行查找，如果找到，則說明這是一個(gè)命中，否則跳過。假設(shè)單詞長度為W，序列生成單詞數(shù)目為v。序列數(shù)據(jù)庫共有n條序列，序列平均長度為l，假定沒有長度小于w的序列，則每個(gè)序列能夠產(chǎn)生的片段數(shù)為l-w+1。總共產(chǎn)生片段數(shù)為n(l-w+1)。對于每一個(gè)片段，BLAST都要從單詞庫中查找一次，在對單詞表排序的情況下，采用二分法查找則使用的平均時(shí)間復(fù)雜度為O(lnv)。因此整個(gè)過程時(shí)間復(fù)雜度為：

無論比對的序列為蛋白序列還是DNA序列，其W值通常為3～12左右。而數(shù)據(jù)庫中序列的長度通常都是超過100。因此可以忽略上式的w值。可記時(shí)間復(fù)雜度為：

實(shí)際上，n與l的乘積代表序列數(shù)據(jù)庫的總大小，用C代替，上式可化簡為：

3 基于位圖索引和B+樹的BLAST改進(jìn)算法

由式（1）知上述過程耗時(shí)取決于數(shù)據(jù)庫大小C與生成單詞數(shù)目v。其搜索速度與數(shù)據(jù)庫大小成線性關(guān)系。而現(xiàn)代生物數(shù)據(jù)庫發(fā)展的速度將會(huì)對計(jì)算資源提出巨大的要求。因此逐字順序檢索數(shù)據(jù)庫序列是不現(xiàn)實(shí)的。而對于一個(gè)DNA查詢序列，其生成的單詞數(shù)理論最多有4w個(gè)，當(dāng)W=11時(shí)約為4×106個(gè)。如果能夠找出一種辦法，能夠加快單詞從數(shù)據(jù)庫中找到匹配的時(shí)間，則效果可能會(huì)好很多。

如果在搜索數(shù)據(jù)庫中對所有可能的單詞建立其所在數(shù)據(jù)庫位置的表，并對這個(gè)表對單詞字段建立位圖索引[4-5]，則可以極大地加快單詞匹配的速度。為了加快單詞的檢索，可以對單詞字段再建立一個(gè)表，采用B+樹[6-7]存儲(chǔ)。如此可以大幅降低單詞匹配的耗時(shí)。

改進(jìn)的算法主要在于對序列數(shù)據(jù)庫存儲(chǔ)方式的改進(jìn)。首先把數(shù)據(jù)庫序列進(jìn)行編號(hào)，如果序列長度過大，可能還要分段。在此，規(guī)定長度大于L=250的數(shù)據(jù)庫序列都進(jìn)行分段，使長度大小不超過L。設(shè)編號(hào)總數(shù)為N，即該數(shù)據(jù)庫有N段。然后對于每一個(gè)可能的單詞，都根據(jù)數(shù)據(jù)庫生成一個(gè)長度為N的位向量。位向量的生成方法是查找數(shù)據(jù)庫的每一段，如果第i條有該單詞的命中，則對應(yīng)位向量第i個(gè)元素為1，否則記0。

以表1的DNA數(shù)據(jù)庫作為例子。這個(gè)數(shù)據(jù)庫只包含3條序列，而且每條長度為60，因此不需要對序列進(jìn)行分段。

首先對數(shù)據(jù)庫按序列進(jìn)行編號(hào)。這里片段總數(shù)N=3。然后第1個(gè)DNA段第1個(gè)長度為W的片段TACTTGCCAAG，其在第1個(gè)DNA段有1個(gè)命中，因此位向量第1個(gè)元素為1。對于第2個(gè)DNA段，片段TACTTGCCAAG并沒有命中，因此位向量第2個(gè)元素為0。依此類推，片段TACTTGCCAAG的位向量為100。然后對于這個(gè)DNA段的第2個(gè)長度為W的片段ACTTGCCAAGT，因?yàn)樗诘?個(gè)DNA段有命中，因此這個(gè)單詞位向量第3個(gè)元素置1。可得出這個(gè)片段的位向量為101。依此類推，遍歷整個(gè)數(shù)據(jù)庫之后，可以得出一張“單詞-位向量”表。

為了更快地進(jìn)行查找。可以為這張單詞-位向量表進(jìn)行進(jìn)一步的索引。例如以“單詞”為鍵，建立B+樹，樹的葉子節(jié)點(diǎn)指向這個(gè)單詞所對應(yīng)的位向量。假如建立的B+樹的一個(gè)塊有255個(gè)指針，那么3層255階的B+樹能夠有2553≈1.66×107個(gè)指向記錄的指針。當(dāng)BLAST數(shù)據(jù)庫建立的索引表單詞長度W=11時(shí)，上限單詞數(shù)目為411≈4.2×106，小于上述B+樹可提供指針數(shù)目，因此只需要建立一個(gè)三層B+樹就能夠進(jìn)行索引了。

表1 DNA序列數(shù)據(jù)庫

當(dāng)BLAST根據(jù)查詢序列生成單詞表時(shí)，給出任意一個(gè)長度為W的單詞，可以根據(jù)B+樹找到單詞對應(yīng)的位向量，然后利用位向量知道數(shù)據(jù)庫第幾條DNA段有命中。從而快速地找到命中的序列。

圖2 三階B+樹

4 算法分析

假定B+樹每一個(gè)塊都是充滿的（事實(shí)上，因?yàn)楫?dāng)W=11時(shí)單詞數(shù)目的上限小于葉節(jié)點(diǎn)指針數(shù)，因此這個(gè)B+樹是不可能充滿的），那么給出一個(gè)單詞要找到其在B+樹對應(yīng)的鍵時(shí)，程序需要從根節(jié)點(diǎn)開始逐層往下比較，每層都順序比對所有的鍵，平均比對次數(shù)為255/2次，一直到葉節(jié)點(diǎn)所在的第三層，共需比對255/2×3=382.5次。

找到對應(yīng)位向量之后，程序只需讀入位向量一次，然后根據(jù)位向量查找命中的序列段。每在位向量中發(fā)現(xiàn)一個(gè)1，則說明有一條序列內(nèi)部有命中。然后遍歷這一段序列，找到所有的命中。設(shè)數(shù)據(jù)庫有N個(gè)序列段，每個(gè)序列段長度最大值為L，假定數(shù)據(jù)庫每個(gè)序列段的長度都達(dá)到最大，從而對于每一個(gè)有命中的序列段，程序都需要比對(L-W+1)次以找出命中的具體位置。

又設(shè)查詢序列能夠生成v個(gè)單詞，每個(gè)單詞在數(shù)據(jù)庫中有h個(gè)序列段是有命中的，那么平均每個(gè)單詞查找命中需要比對的次數(shù)為：

這里的h可以用任意長度為W的單詞在數(shù)據(jù)庫有命中的序列段個(gè)數(shù)的期望近似代替，即：h?H=N·P。其中P是任意長度為W的單詞在數(shù)據(jù)庫的任意一個(gè)序列段有命中的概率近似估計(jì)（因?yàn)椴煌瑔卧~在一個(gè)隨機(jī)序列出現(xiàn)的概率不同，因此這里計(jì)算概率時(shí)不考慮重復(fù)情況只給出概率上限）：

因此平均每個(gè)單詞需進(jìn)行比對的次數(shù)可估計(jì)為：

而查詢序列生成v個(gè)單詞，兩步共需比對次數(shù)為：

因?yàn)閿?shù)據(jù)庫的序列條數(shù)通常十分巨大，因此可以省略382.5，同時(shí)對于較大的L值，W值通常比較小。因此T可以近似估計(jì)為：

這里N與L的乘積實(shí)際上就是數(shù)據(jù)庫的規(guī)模，即C=N×L，則上式可化簡為：

把式（3）與式（1）比較，可以得出T與T0的關(guān)系：

取W=11，L=250，則 β與v有關(guān)，可以繪出 β(v)的圖像如圖3。

由圖3可知隨著查詢序列長度的增長原算法與改進(jìn)算法耗時(shí)比在降低，即對于較長的查詢序列，改進(jìn)算法的效果將會(huì)變差，不過對于DNA序列，查詢序列的長度一般小于1 000，根據(jù)圖像可得 v＜1 000時(shí) β(v)＞100，即原算法耗時(shí)T至少為改進(jìn)算法耗時(shí)T0的100倍，換句話說，執(zhí)行任何一次搜索的時(shí)候，算法查詢命中所耗費(fèi)時(shí)間將變?yōu)樵瓉淼?/100。

除此之外，只讀入有命中的序列段而不是全部讀入整個(gè)序列數(shù)據(jù)庫將會(huì)節(jié)省大量的時(shí)間。另外，為了更快地進(jìn)行檢索，還可以采用把整個(gè)B+樹索引載入內(nèi)存，這樣只需要一次的磁盤I/O就能找到單詞所對應(yīng)的位向量。

5 速度提升的代價(jià)

事實(shí)上，為序列數(shù)據(jù)庫建立位圖索引將會(huì)需要大量空間用于存儲(chǔ)單詞-位向量表，即使采用壓縮的算法[8]存儲(chǔ)位向量，根據(jù)對幾個(gè)生物數(shù)據(jù)庫的建表實(shí)現(xiàn)的結(jié)果，建立單詞-位向量表時(shí)最好預(yù)留30倍原文件大小的空間。也即是說，對于一個(gè)大小為1 GB的序列數(shù)據(jù)庫，至少需要30 GB的硬盤空間用于存放改進(jìn)的數(shù)據(jù)庫。而這正是速度提升的代價(jià)所在。

不過鑒于現(xiàn)今PB級的存儲(chǔ)設(shè)備已相當(dāng)普及，相對于為了加快搜索速度而增添計(jì)算資源，采購大容量存儲(chǔ)設(shè)備反而更加劃算，因此這個(gè)代價(jià)是完全可以承受的。

另一個(gè)速度提升的代價(jià)在于對序列數(shù)據(jù)庫的預(yù)處理時(shí)間，易知如果需要為一個(gè)序列數(shù)據(jù)庫建立位圖索引，將會(huì)需要大量的時(shí)間用于建立單詞-位向量表。不過在實(shí)際應(yīng)用中，對于一個(gè)序列數(shù)據(jù)庫來說，建表將會(huì)是一次性的。因此一旦建表完成，之后的任何一次序列搜索將都會(huì)獲得速度的提升。而B+樹的特性又能夠使得對數(shù)據(jù)庫的日常維護(hù)變得十分容易，例如，當(dāng)要修改某條序列時(shí)，只需要考慮受到影響的單詞的位向量即可。

6 結(jié)束語

本文分析了BLAST算法其在播種階段的計(jì)算步驟，并提出采用位圖索引和B+樹重新組織數(shù)據(jù)庫的方法提高計(jì)算速度。采用位圖索引的思想對序列數(shù)據(jù)庫進(jìn)行處理，生成單詞-位向量表，在使用B+樹進(jìn)行二次索引，從而改變BLAST在播種階段的處理方式。理論證明，新的算法能夠使BLAST查找命中耗費(fèi)的時(shí)間至少縮小為原算法的百分之一。

但還有一些問題尚待進(jìn)一步研究，例如改進(jìn)序列數(shù)據(jù)庫的存儲(chǔ)效率以及改進(jìn)算法實(shí)現(xiàn)兩個(gè)問題，這將是下一步研究的重心。

[1]Altschul S，Gish W.Basic local alignment search tool[J].Journal of Molecular Biology，1990，215：403-410.

[2]Altschul S，Madden T，Sch?ffer A.Gapped BLAST and PSIBLAST：a new generation of protein database search programs[J].Nucleic Acids Research，1997，25（17）：3389-3402.

[3]Korf I，Yandell M，Bedell J.BLAST[M].America：O’Reilly Media，2003.

[4]Bayer R，McCreight E M.Organization and maintenance of large ordered indexes[J].Acta Informatica，1972，1：173-189.

[5]Comar D.The ubiquitous B-tree[J].Computing Surveys，1979，11（2）：121-137.

[6]O’Neil P.Model 204 architecture and performance[C]//Proceedings of the 2nd International Workshop on High Performance Transaction Systems，Asilomar，CA，USA，1987：40-59.

[7]O’Neil P，Quass D.Improved query performance with variant indexes[C]//SIGMOD’97.New York：ACM Press，1997：38-49.

[8]Wu K，Otoo E，Shoshani A.Optimizing bitmap indices with efficient compression[J].ACM Transactions on Database Systems，2006，31（1）：1-38.

基于位圖索引和B+樹的BLAST改進(jìn)算法

黃志洪1，呂 威2，黃 俊1

HUANG Zhihong1,LV Wei2,HUANG Jun1

1.School of Mathematics&Computational Science,Sun Yat-sen University,Guangzhou 510275,China
2.School of Information Technology,Beijing Normal University Zhuhai Campus,Zhuhai,Guangdong 519085,China

It concentrates on the time consuming procedure that goes through the database in the first step of BLAST.In order to speed up the program,it introduces a new approach which using bit map index and B+tree.The developed method builds up a word-bit_vector table according to the database,and reorganizes it with B+tree.It proves theoretically that it decreases the word searching time of BLAST substantially.

sequence alignment;BLAST algorithm;bitmap index;B+tree

針對BLAST算法在查找命中的過程中需要遍歷數(shù)據(jù)庫造成計(jì)算資源消耗的問題，提出了基于位圖索引和B+樹的數(shù)據(jù)存儲(chǔ)方式以加快數(shù)據(jù)的檢索。改進(jìn)算法利用位圖索引的原理建立數(shù)據(jù)庫的單詞-位向量表，并對這個(gè)表使用B+樹再次進(jìn)行索引，最終達(dá)到加快BLAST程序的運(yùn)算速度。對于DNA序列這個(gè)方法能夠使BLAST查找命中耗費(fèi)的時(shí)間得到極大的減少。

序列比對；BLAST算法；位圖索引；B+樹

A

TP391

10.3778/j.issn.1002-8331.1110-0392

HUANG Zhihong,LV Wei,HUANG Jun.Improved BLAST algorithm based on bitmap indexes and B+tree.Computer Engineering and Applications,2013,49（11）：118-120.

國家自然科學(xué)基金（No.11126039）。

黃志洪（1967—），男，講師，主要研究領(lǐng)域?yàn)閿?shù)據(jù)庫、數(shù)據(jù)分析。E-mail：luwei00@126.com

2011-10-20

2012-02-06

1002-8331（2013）11-0118-03

CNKI出版日期：2012-07-19 http://www.cnki.net/kcms/detail/11.2127.TP.20120719.1511.006.html