基于CUDA平臺的數據庫序列比對算法加速

摘要：在生物信息學中，數據庫序列比對是極為常用的操作，Smith-Waterman算法是最流行的序列比對算法，精確度高，但是計算復雜度高，在進行大量的序列比對非常耗時。另外，生物技術的發展使得已知的序列數據庫變得越來越龐大，這導致進行數據庫序列比對所消耗的時間也越來越長，因而有必要加速數據庫序列比對算法。NVIDIA提出了CUDA編程架構，相比之前的GPGPU具有更好的可編程性，用戶可以更輕松地發掘出GPU強大的計算能力。在CUDA平臺上實現了Smith-Waterman的數據庫序列比對算法的并行加速，速度優于已有的基于GPU的實現，超過了基于啟發式算法的BLAST算法執行速度。

關鍵詞：序列比對； Smith-Waterman算法； CUDA； GPU計算

中圖分類號：TP37 文獻標識碼：A文章編號：2095-2163（2013）02-0044-06

0引言

數據庫序列比對是生物信息學研究領域中基本性的熱點應用之一。研究人員經常需要面對將新發現的基因序列與已知的基因數據庫進行搜索比對的課題狀況。而通過基因序列比對，或者可以找出兩個序列間的相似區域，判斷該基因序列的新穎性，或者將兩類生物體進行全基因組（full genome）序列比對，以判斷兩者間的同源性等。在已有的眾多序列比對算法中，Smith-Waterman算法[1]最為著名。該算法基于Needleman-Wunsch算法[2]的改進，可用于局部序列比對，能夠得到精確的序列比對結果。而Smith-Waterman算法在經過Gotoh[3]的進一步完善后，算法復雜度達到了O（M×N）。其中，M和N分別為兩個比對序列的長度。但是，Smith-Waterman算法仍然表現了復雜度大，計算耗時的不足，而與此同時，隨著各類生物序列數據庫規模的逐年增大，其計算速度已經越來越難以滿足實際應用要求。

近年來，出現了很多加速方法，按其原理可將其組織為兩類：一是在算法層次進行改進，在降低算法復雜性后，提高速度；二是通過硬件加速平臺，提升算法速度。

對于算法改進，BLAST[4]和FASTA[5]是算法改進后的典型代表。這些算法都基于啟發式算法，可查找出兩個比對序列中同時存在的某個特定長度的字或者種子，對這些字或種子所在的區域實現擴展，進行相似性比對。只有當種子在兩個序列中都有命中的區域時，才進行拓展比對，因此大大縮短了運行時間，但是卻無法得到與Smith-Waterman算法相同的精確結果。

對于通過特殊硬件平臺加速，P-NAC[6]、BISP[7]以及SAMBA[8]等利用FPGA硬件實現核心算法，速度快，但是價格昂貴，缺乏靈活性，不易于擴展。在多處理器上并行Smith-Waterman算法也是常用的方式[9-12]，但是商用多核處理器價位較高，相對低價的GPU逐漸成為時新的加速平臺。Liu[13]利用GPGPU實現了基于Smith-Waterman的數據庫序列比對，并相比于SSEARCH獲得了多達5倍的性能提升，而Manavski[14]則第一次利用NVIDIA的計算統一架構（CUDA）平臺付諸實現，相較Liu的研究結論又獲得18倍的加速。之后，CUDA編程環境的GPU即成為研究人員聚焦關注的平臺之一[15-17]，采用CUDA編程模型[18]，可方便地完成劃分和并行任務。本文在基于CUDA平臺的NVIDIA GeForce 8800GT上加速Smith-Waterman算法，獲得了約3 000MCUPS（Million Cell Updates Per Second）的比對速度，其性能超過了基于啟發式算法的BLAST在通用CPU上的速度。

本文組織如下：第1節和第2節分別簡單介紹基于Smith-Waterman的序列比對算法，以及NVIDIA 8800GT顯卡的結構特征和CUDA編程環境，第3節提出了基于CUDA編程環境的算法實現和優化，并對性能數據進行了分析，第4節概述本文的相關工作，最后一節總結全文工作。

1序列比對算法

數據庫序列比對是將查詢序列與數據庫中各個序列通過逐一運行序列比對算法，在完成比對后計算得到各次比對的結果分值，如果分值超過了指定閾值，就說明數據庫中的該條待檢序列與查詢序列匹配成功，可以進入后續的其他操作流程。

序列比對算法的輸入為兩個序列文件，其輸出為序列間的相似區域，下面是一個蛋白質序列比對的例子，具體細節如圖1所示。第2期袁競杰，等：基于CUDA平臺的數據庫序列比對算法加速智能計算機與應用第3卷 圖1 蛋白質序列比對示例

Fig.1 Example of protein sequence comparison

由圖1可見，字符串代表由20個不同氨基酸組成的蛋白質序列，“|” 代表序列間這兩個字符完全匹配，“：”代表序列間這兩個字符近似匹配，可以根據指定的置換矩陣（substitution matrix）來判斷哪些字符屬于近似匹配，“-“表示插入了一個空位。

Smith-Waterman算法是進行序列比對的有效方法，通過利用動態規劃方法來計算最優的局部序列匹配，而且多是采用Gotoh改進后的算法。算法具體過程如下：

由式（1）可以看出，矩陣中各個點的計算存在著嚴格的依賴關系，每個點的計算都取決于同一行左側一列的點、同一列上一行的點以及左側一列、上面一行的左上方點的數據值。因此，在實際計算時，可以按照行、列或者是反對角線的順序依次計算。算法復雜度為O（|A||B|）。

2GPU體系結構和CUDA編程環境

根據圖像應用要求，GPU的體系結構能高效、可靠地支持數量龐大的細粒度的線程調度，而這一點在傳統的通用并行計算領域也存在著廣闊的應用空間。CUDA是NVIDIA于2006年推出的計算架構，憑此可發掘GPU強大的計算能力，給密集型應用計算帶來了新的解決契機。

如圖2所示，GeForce 8800 GT芯片由14個流多處理器（SM）組成，每個SM均包含有8個小的處理器核（processor），這8個小核共享一個指令單元，可通過SIMD的方式執行指令，運行頻率為1.5GHz。同時，每個小核又包括一個32-bit的單精度浮點乘/加計算部件。此外，每個SM上還有2個特殊功能單元（SFU），用來完成復雜的浮點運算，如開方、求取正/余弦等。由此，GPU的整體運算能力已可達到378GFLOPS。

CUDA是對標準C/C++的擴展，并提供了一些額外的數據類型、關鍵字以及函數調用方式等。在CUDA模式下，GPU就是一個專門用于加速的協處理器，程序中的串行部分由主機CPU來完成，而并行部分則包裝成一個kernel函數，由主機調用GPU執行。CUDA上的線程組織方式分為grid-block-thread三級[19]。最頂層是grid，執行同一kernel函數的全部thread構成一個grid，Grid可包含多個block，每個block又包含多個thread。通常，一個SM上最多含有768個thread，而對block，則是最多容納8個。執行時，位于同一block內的thread組織成warp形式，每個warp可包含32個連續thread。warp是GPU調度的基本單位，當SM上的小核執行到長延時指令時，系統將快速切換到另一個可以待命執行的warp，直到SM上沒有可執行的warp時，則暫停，以此減少等待耗時。表1GeForce 8800GT 存儲系統

Tab.1 GeForce 8800GT memory system

存儲類型位置大小命中時延是否只讀作用域Register片上8 192/SM1 拍否Thread內部Shared memory片上16KB/SM幾拍否Block內部Global memory片外，無片上cache512M200多拍否全局Constant memory片外，有片上cache總共64KB幾拍是全局Texture memory片外，有片上cache8KB/SM幾拍是全局3算法實現及性能評估

對于訪存性能提升策略，CUDA與通用CPU存在很大的不同。CPU主要是通過設置片內cache來加速訪存，而CUDA遇到長延時訪存指令時，將快速切換到其他線程執行，切換過程開銷很小，因此，需要創建足夠多的線程用來實現切換調度，從而充分利用訪存操作的長延時。

在進行數據庫比對時，查詢序列和數據庫序列的比對操作相互完全獨立，直接的方法就是使每個線程都要完成一個查詢序列與一個數據庫序列的比對過程，這樣，m個查詢序列和n個數據庫序列就可以創建m×n個線程。實際數據庫中，n值常常較大，數量級在5以上，如此數目的線程數完全滿足CUDA的調度需求。本文實驗使用Swiss-Prot[19]蛋白質數據庫（數據庫共含有366 226個蛋白質，132 054 191個氨基酸），在數據庫中提取查詢序列，Smith-Waterman算法所使用的置換矩陣為BLOSUM50，首次空位插入開銷為10，后續空位插入開銷為2。

3.1訪存操作的優化

CUDA上訪存操作對性能的影響主要表現在兩個方面，一是訪存延時大，需要超過200個時鐘周期；二是對訪存帶寬的要求高，GPU的計算能力強達378GFLOPS，但是訪存帶寬卻不到100GB/S，如果計算訪存比低，將導致帶寬形成短板。因此，除了需要創建大量調度線程來匹配訪存延時外，還應該盡量減少程序中的訪存操作。這樣做目的主要有三個：一是降低對訪存帶寬的需求；二可降低對線程數的需求，避免出現活動的線程數不足引發的等待；三能保持程序的連貫性，減少線程執行中的切換。

在本文的實現中，借鑒了文獻[11]和[16] 中采用的query-sequence profile方法，即根據查詢序列和置換矩陣計算得到的每種字符對應的匹配分值，再按照查詢序列的字符排列完成連續存放好。如圖4所示，列序列是一個蛋白質查詢序列示例，行字符代表氨基酸的種類（A～Z表示）。每一列記錄了查詢序列對應某種蛋白質的比對分值，矩陣中所有數據按列連續存放。當查詢序列和數據庫序列的某一個字符進行比對計算時，該字符對應的所有置換分值均可從內存中連續讀到，無需再查找置換矩陣。query-sequence profile方法可以減少（1）、（3）種訪存操作的總數量，同時，將隨機訪問轉換成連續訪問能夠更好地利用cache或者數據讀寫合并，來實現加速訪存。

profile信息表大小取決于查詢序列和序列中字符種類的乘積。大多情況下，shared memory、constant memory以及texture memory都可以滿足需要，但不同的數據存放方式決定了其性能也將有所不同。在CUDA模式下，可以顯式指定數據的存放位置。圖5給出了3種對應不同查詢序列長度時的存儲結果，可以看到，texture memory在這些長度下的速度都是最快的。下面將分析討論最佳的存儲方式。

使用shared memory時，要避免bank沖突，沖突發生時訪問速度將會下降。而本文使用的查詢序列長度卻會使得warp內的并行thread訪問同一bank，導致所有訪問串行發生，勢必影響了數據讀取速度。此外，受到shared memory容量的限制，查詢序列不可能太長；而且寫shared memory時，所有線程都需要增加一個同步操作，同步操作會對單個block內的最佳thread數帶來影響。

而對constant memory的訪問，warp內的線程訪問同一地址時，速度最快，否則，訪問時間將隨著訪問的地址數呈線性增長。在本文中，形式各異的線程都在和數據庫中的不同序列進行比對，因其訪問不同種類字符的概率很高，決定了每個線程的訪問地址都不一致。隨著查詢序列長度的增加，沖突越發明顯，性能也逐漸變差。圖5（b）顯示，在查詢序列長度達到127后，其比對速度即開始低于shared memory的性能，使用constant memory的擴展性也不好。

第四種訪存操作針對計算時臨時生成的兩列H值和E值（假設按列進行計算，按行計算則為H值和F值），需要寫HE_out和讀HE_in。由表1可知，支持讀寫操作的只有global memory、local memory或者shared memory。計算過程中，每個線程均需要保存長度等于查詢序列的兩列數據， 首先，shared memory空間限度無法滿足該要求。若使用global memory，需要預先做好分配，當kernel中的線程數目較多時，同樣會出現空間無法應對的問題，雖然可通過減少單個kernel線程數并增加kernel數來完成相同的任務，但線程數的減少卻容易造成調度困難，從而影響整體效率提升。另外，多個kernel函數還將引入創建和退出的額外開銷。實驗統計表明，只有在查詢序列較短時，global memory利用向量訪存操作性能高于local memory。因此，文中采用local memory來保存H值和E值。但是local memory沒有cache，速度很慢，因此要盡量減少訪問次數。實際上，計算出來的分值都不大，用16位整數完全可以表示，將這種位數較少的讀寫操作進行合并，可以明顯減少訪存次數。目前，CUDA尚不支持對local memory向量讀寫，因此可將2個16位的訪存進行合并，使用一條32位的讀寫指令即可完成。合并后的性能比SW_prof提高了66.3%，具體可參見圖6中的SW_HEpk。

與對HE值的合并類似，第二種訪存操作也可以合并。數據庫序列都是8bit字符，存在global memory中，CUDA提供對其128位向量的取數操作，可一次提取16個字符到寄存器，這樣就將原有的|B|次訪存降為|B|/16次。結果性能如圖6中的CUDA.dbpk項所示，相比SW_HEpk又獲得了18.8%的提升。

3.2其它優化

在CUDA模式下，一個warp內的所有thread采用SIMD的形式，執行相同的指令。為了充分利用計算資源，同一warp內thread的計算量應該盡量接近，使負載達到均衡。在執行算法時，計算量與兩個序列的長度乘積成正比，且不同thread的查詢序列相同，因此就需要保證相鄰thread得到的數據庫序列長度基本接近。為此，將對所有數據庫序列進行預先降階重排，降階排列可以使得比對序列長、計算任務重的thread所在的block優先執行，留下的則是計算任務較輕的block，這樣，GPU的多個SM就可以在更加接近的時間內全部完成，從而實現良好的負載均衡。其性能如圖6中的SW_sort所示，這是在SW_init的基礎上進行降排序得到的結果，排序后獲得了1.82倍的速度提升。圖6中的SW_sort右側的各項結果都是在已排序的前提下而得到的。

在其他資源確定不變情況下，需要盡可能地提高GPU的占用率（occupancy），也就是要提高每個SM上正在執行的thread數（以warp為單位）。不同數據庫序列和查詢序列之間的比較完全獨立，block內的thread并不存在同步的模式，所以可將所有thread組織在一個block內，使得單個SM上活動的thread數盡可能多，由此提高占用率。

最后，根據編譯器生成的PTX文件，對kernel函數的內存循環代碼進行手工的優化，重點在于減少寄存器的使用以及消除冗余的操作。此時，發現代碼中存在大量的數據轉換指令，可通過調整變量的類型，而將其消除。優化過后，每個線程使用的寄存器由24個減少到21個，而每個SM上的可活動線程數由320增加到了384個，由此獲得了大約4%的性能提升，如圖6中的SW_tune所示。圖6中的其它結果卻是在每個block包含320個thread的情況下得到的。

經過上述優化后，最終的計算訪存比就達到了18：1。在這種計算訪存比的意義作用下，可假設訪問device memory的開銷是300個時鐘周期，那么每個SM上只需要6個warp，即192個thread就可以匹配訪存延時，而目前卻有384個thread，已經完全能夠滿足為匹配訪存延遲所需進行的規定調度，因此，在訪存上已經不會存在明顯的性能瓶頸。

3.3性能對比

在實驗最后，將本文得到的結果和已有的加速方法進行了對比，對比結果如表2所示。表中，SW_tune為本文的結果實現，其配置為：使用了單塊NVIDIA GeForce 8800GT型號的GPU，主機為Intel Core 2 6300的雙核CPU，運行頻率1.86GHz。Nvcc版本為release 1.1， V0.2.1221，且程序使用了-O3優化。SW_SAM為文獻[16]中基于單塊GeForce 8800GTX的CUDA實現。SW_Liu來源于文獻[17]，是基于GeForce 7800GTX，使用OpenGL接口的實現。Ssearch和Blast的結果數據均來自文獻[16]，基于3 GHz的Intel Pentium 4處理器平臺。表 2 與已有的其它實現方法的性能對比（單位：MCUPS）

Tab.2 Performance comparison with existed methods（Units：MCUPS）Seq nameLengthSW_tuneSW_SAMSW_LiuSsearchBlastO29181633 0691 8491971191 488P036301272 9951 8893171191 948P537652553 0151 8114281212 027Q8ZGB43612 9171 810486114 1 936P582295112 9241 7955331232 691由表2可以看到，本文的實驗實現在性能上均優于其它幾類的實現成果。SW_Liu是基于GPU的Smith-Waterman實現，本文在對應不同查詢序列的長度下實現速度是前者的5.5～15.6倍。SW_SAM也是基于CUDA的實現，但其在降低訪存開銷、減少寄存器使用以及block/thread的組織等方面并未得到充分優化，與該方法相比，本文依然獲得了約60%的性能提升。 SW_SAM實驗中的GPU包含16個SM，文中的GPU只包含14個SM，如果在相同的平臺背景下，本文方法的優勢將更加明顯。另外，即使相比于BLAST[6]這種快速的基于啟發式的算法，本文也穩定地獲得8.7%～106%的速度優勢。

4結束語

本文在基于NVIDIA GeForce 8800GT平臺的CUDA模式下實現了Smith-Waterman算法的加速，使用了一系列方法優化程序的訪存操作、線程劃分和組織等。結論顯示，算法在得到精確結果的同時，也保證了較快的速度，達到了約3 000MCUPS，現已超過了基于啟發式的BLAST算法在主流通用CPU上的性能。

研究工作表明，CUDA平臺對于生物信息學中的序列比對算法獲得了良好的加速效果，因此，未來還可以對其它類型的序列比對算法實現加速。同時，也要繼續將CUDA的加速效果和在當前主流使用的多核/眾核平臺上的實驗結果進行對比研究，以分析不同平臺的優、劣特性，為今后計算平臺的選擇提供全面、有益的參考依據。

參考文獻：

[1]SMITH T F， WATERMAN. SM. Identification of common molecular subsequences[J]. Jounal of Molecular Biology， 1981， 147（1）： 195-197.

[2]NEEDLEMAN S B， WUNSCH C D. A general method applicable to the search for similarities in the amino acid sequence of two proteins[J]. Jounal of Molecular Biology， 1970（48）： 443-453.

[3]GOTOH. An improved algorithm for matching biological sequences[J]. Jounal of Molecular Biology， 1982， 162：705-708.

[4]ALTSCHUL S， MADDEN T， SCHAFFER A， et al. Gapped BLAST and PSI-BLAST ： A new generation of protein database search programs[J]. Nucleic Acids Research， 1997， 25（17）： 3389-3402.

[5]PEARSON W R， LIPMAN D J. Improved tools for biological sequence comparison[J]. Proceedings of the National Academy of Sciences， USA，1988， 85（8）： 2444-2448.

[6]LOPRESTI D. P-NAC. A systolic array for comparing nucleic acid sequences[J]. Computer， 1987， 20（7）： 81-88.

[7]CHOW E， HUNKAPILLER T， PETERSON J， et al. Biological information signal processor[C]. //Proc. of the International Conference on Application Specific Array Processor. Barcelona， Spain：IEEE， 1991：144-160.

[8]GUERDOUX P， LAVENIER D. Samba： hardware accelerator for biological sequence comparison[J]. CABIOS， 1997， 13（6）： 609～615.

[9]ROGES T， SEEBERG E. Six-fold speed-up of Smith-Waterman sequence database searches using parallel processing on common microprocessors[J]. Bioinformatics， 2000， 16：699-706.

[10]FARRAR M. Striped Smith-Waterman speeds database searches six times over other SIMD implementations [J].Bioinformatics， 2007， 23（2）： 156-161.

[11]MARTIN W S， DELCUVILLO J B， USECHE F J， et al. A multithreaded parallel implementation of a dynamic programming algorithm for sequence comparison[C]// Pacific Symposium on Biocomputing. Hawaii， 2001.

[12]YAP T K， FRIEDER O， MARTINO L R. Parallel computation in biological sequence analysis[J]. IEEE Transactions on Parallel and Distributed Systems， 1998， 9（3）： 283-294.

[13]LIU W， SCHMIDT B， VOSS G， et al. Bio-sequence database scanning on a GPU[C] //Proceeding of the 20th IEEE International Parallel Distributed Processing Symposium：2006（IPDPS2006）（HICOMB Workshop）. Rhode Island， Greece：IEEE，2006.

[14]MANAVSKI S A， VALLE G. CUDA compatible GPU cards as efficient hardware accelerators for Smith-Waterman sequence alignment[J]. BMC Bioinformatics， 2008，9（Suppl 2）：S10.

[15]ZHENG F； XU X B； YANG Y H， et al. Accelerating Biological Sequence Alignment Algorithm on GPU with CUDA[C] //Proc. of ICCIS.Chengdu， China： IEEE ， 2011：18 – 21.

[16]LEE S T， LIN C Y， HUNG C L， et al. Using frequency distance filteration for reducing database search workload on GPU-based cloud service[C]//Proc.of CloudCom. Taipei： IEEE，2012：735 – 740.

[17]HAINS D， CASHERO Z， OTTENBERG M， et al. Improving CUDASW++，a Parallelization of Smith-Waterman for CUDA Enabled Devices[C] //Proc. of IPDPS Workshop.Shanghai：IEEE，2011： 490-501.

[18]NVIDIA CUDA. http：//developer.nvidia.com/object/cuda.html

[19]ExPASy， http：//www.expasy.ch/sprot/