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

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

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

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

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

1 引言

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

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

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

2 問題描述

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

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

圖1 命中

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 三階B+樹

4 算法分析

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

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

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

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

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

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

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

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

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

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

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

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

5 速度提升的代價

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

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

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

6 結(jié)束語

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

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

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

序列比對；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)域為數(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