新型圖形硬件支持下的動態規劃局部比對算法加速研究

張 林

(浙江中醫藥大學生命科學學院，浙江杭州310053)

生物序列分析是生物信息學研究的一個重要手段，它是了解基因組結構和功能的基本途徑［1］。而在生物序列分析中，生物序列比對是最常用和最經典的研究手段［2］。基于真核表達的剪接特點和蛋白質的模塊性質，目前，序列局部比對已成為生物序列比對中的主要分析方式。其通過將查詢序列與整個生物數據庫的所有序列進行局部比對，來完成一系列的生物學分析，探索生命起源和存在的內在規律。

而在后基因組時代，生物數據庫含有海量數據。因此，在保證比對準確的同時如何提高其效率，已成為序列局部比對的瓶頸。現在比較重要的是動態規劃局部比對算法［3］，準確性高是其最大優點，但存在著效率問題。為了增加動態規劃局部比對算法的效率，最早提出的解決方案是加入啟發式算法，例如現在常見的FASTA和BLAST算法［4］，但其準確性受到影響。現在的思路是采用專門的生物信息處理機器和專有硬件提高動態規劃局部比對算法的比對速度，但其缺點是投入太大，實現困難。

近年來，圖形硬件(計算機顯卡)發展迅速，主要具有以下優點。第一，浮點運算能力強，顯存與位寬大，非常適合高密度的運算，也非常適合流式處理［5］。第二，并行度很高，具有統一渲染構架，非常適合并行運算。第三，可編程，具備了通用計算能力［6－7］，在圖像處理［8］、計算幾何［9］、科學計算［10］、生物數據分析等領域已有多處應用［11－12］。第四，價格便宜，性價比非常高。

鑒于上述背景，本文擬以動態規劃局部比對算法為研究對象，將其實現，然后利用圖形硬件將其加速，并通過一系列實例進行加速性能評測，探索兼具準確性、高效性、通用性的生物序列局部比對方法。

1 實驗材料與方法

1.1 實驗材料

電腦配置如下:CPU:Intel酷睿i7 4960X(頻率4 GHz、15 MB 三級緩存);內存:DDR3(16 GB);硬盤:1 TB。

3款最新的計算機顯卡(圖形硬件)，品牌為NVIDIA，型號為 Geforce GTX 760、Geforce GTX 770、Geforce GTX 780Ti。安裝在上述配置的電腦上分別運行。

1.2 實驗方法

1.2.1 實現動態規劃局部比對算法

實現動態規劃局部比對算法，即尋找兩條序列局部最優匹配的算法，又稱Smith-Waterman算法。

1.2.2 利用圖形硬件加速動態規劃局部比對算法

將上述算法通過邏輯映射和物理映射移植到圖形硬件上，并研究在圖形硬件上研究對其加速。

1.2.3 加速后測評

查詢對象和數據庫選擇:在UniProtKB/Swiss-Prot蛋白質數據庫中選擇8條不同長度的蛋白質查詢序列，它們的accession編號分別為Q64372(長度73aa)、Q9CQW9(長 度 137aa)、Q60961(長 度233aa)、Q8WU90(長度 426aa)、Q8R2G6(長度949aa)、Q8CJ12(長度 1 009aa)、Q01815(長度2 139aa)、Q0F9K2(長度4 113aa)。將上述8條蛋白質序列比對搜索 UniProtKB/Swiss-Prot數據庫(2014_02版本，序列平均長度375aa，共有1 283 282個蛋白分子)。

實驗組:比對執行動態規劃局部比對算法(Smith-Waterman算法)，在上述三款圖形硬件(見實驗材料)上分別運行并加速。

對照組:作對照的是由CPU執行的SSEARCH程序。它采用的是FASTA啟發式算法，在動態規劃局部比對算法(Smith-Waterman算法)基礎上進行了啟發式算法加速，比對速度快，但準確性受到影響。對照組其他比對條件、電腦配置均與實驗組一致(見實驗材料)。

評價參數:時間性能(執行時間)、MCUPS(Million Cell Updates Per Second)性能值。其中MCUPS為“每秒百萬次格點更新”，在生物序列比對中，其代表了每秒能計算完成的得分矩陣的格點的數量，其值越高，說明比對的速度越快、效率越高。設查詢序列長度為L，比對的數據庫中有n個分子序列，其中第i個序列的長度為Li，查詢序列比對搜索數據庫所耗費的時間為t。則

2 實驗結果

2.1 實現動態規劃局部比對算法

2.1.1 動態規劃局部比對算法基本描述

設有兩條生物分子序列，分別為 Y=Y［1］…Y［p］(長度為 p)和序列 Z=Z［1］… Z［q］(長度為 q)進行比對。采用簡單打分:殘基匹配得2分，不匹配與空位均扣1分，利用c()函數來表示序列殘基間的得分。

最優比對得分矩陣(簡稱得分矩陣)的設置如下。設E為得分矩陣，E(i，j)表示序列Y的前綴Y［1］…Y［i－1］Y［i］和序列Z 的前綴Z［1］…Z［j－1］Z［j］之間的最佳比對得分(1＜=i＜=n，1＜=j＜=m)。E具備初始條件:E(i，0)=E(0，j)=0。

2.1.2 利用遞歸打分計算得分矩陣E

遞歸打分函數如下:

通過上述運算，獲得E。

2.1.3 利用回溯獲得最優比對序列

回溯算法實現如下:

2.2 利用圖形硬件加速動態規劃局部比對算法研究

2.2.1 算法障礙與瓶頸分析

動態規劃局部比對算法主要包含兩個步驟，一是遞歸打分計算得分矩陣E，二是回溯獲得最優比對序列。其中，步驟一的計算復雜度為O(p*q)，步驟二的計算復雜度為O［max(p，q)］。因此，步驟一是瓶頸，需利用圖形硬件進行加速，而步驟二對計算時間的影響幾乎為0，可以在繼續在CPU上執行。

2.2.2 遞歸打分計算得分矩陣的映射過程

經過分析發現，在遞歸打分計算得分矩陣E時，任意格點值的計算，只與其左邊、上方、斜上方的三個格點相關(見圖1)。因此，可將E的計算作為流式處理，其中，打分過程映射為流處理核心(Kernel)、計算得到的E數據映射為流數據(Stream Data)，循環進行(見圖2(a))，此為第一步——邏輯映射。此外，我們還要完成第二步，物理映射，將上述邏輯模型移植到圖形硬件上完成。圖形硬件具有良好的流處理能力，其片段著色器(Fragment Processors)可完成Kernel功能，紋理系統則可以紋理數據的形式保存Stream Data，循環操作完成E的計算(見圖2(b))。

圖1 遞歸打分中的數據相關性Fig.1 Data relativity of recursion marking

圖2 遞歸打分的映射過程Fig.2 The mapping process of recursion marking

另外，在真實情況下，生物序列之間的打分規則會比較復雜，形成打分矩陣。而對于蛋白質分子的比對，由于蛋白質的進化性質，還需引入不同的取代矩陣進行打分。鑒于此，可將打分矩陣也保存在紋理中，供片段著色器讀取。圖3給出了遞歸打分計算得分矩陣算法的最終映射過程。

圖3 遞歸打分算法的最終映射過程Fig.3 Final mapping process of recursion marking algorithm

2.2.3 利用逆對角線法對遞歸打分進行并行加速

通過上面分析發現，在遞歸打分計算得分矩陣E時，任意格點值的計算，只與其左邊、上方、斜上方的三個格點相關，因此每一條逆對角線上格點的計算互不干擾(見圖4(a))。所以，我們可對同一條逆對角線上的所有格點并行計算，而圖形硬件恰好具備良好的并行計算的能力，可完成此項操作(見圖4(b))。因為逆對角線的條數為p+q－1，從而將計算復雜度從原來的O(p*q)降低為O(p+q－1)。

圖4 利用逆對角線在圖形硬件上進行加速計算Fig.4 Accelerated computing by the inverse diagonal rule on the graphics hardware

綜上所述，映射后利用圖形硬件計算得分矩陣 并實現加速的具體算法如圖5:

圖5 在圖形硬件中遞歸計算得分矩陣Fig.5 Computing scoring matrix on the graphics hardware

2.3 動態規劃局部比對算法加速后性能測評

2.3.1 執行時間

表1中顯示的是:8條長度不同的序列分別由SSEARCH、Geforce GTX 760、Geforce GTX 770、Geforce GTX 780Ti執行數據庫比對搜索所耗費的時間。

表1 時間性能評測結果Table 1 Performance evaluation results in time

2.3.2 MCUPS 性能值評測

8條長度不同的序列分別由SSEARCH、Geforce GTX 760、Geforce GTX 770、Geforce GTX 780Ti執行數據庫比對搜索的MCUPS性能值見圖6。

圖6 MCUPS性能比較Fig.6 Performance comparison in MCUPS

3 討論

3.1 時間性能評測分析

第一，無論是利用SSEARCH還是圖形硬件進行序列比對，比對所需要時間都與序列長度成正比，序列長度越長，比對的時間也就越長，符合比對規律。

第二，NVIDIA Geforce GTX 760、NVIDIA Geforce GTX 770、NVIDIA Geforce GTX 780Ti中任意一款圖形硬件的比對時間，都要大大短于SSEARCH。

第三，在 NVIDIA GeforceGTX 760、NVIDIA Geforce GTX 770、NVIDIA Geforce GTX 780Ti三塊圖形硬件之間，比對所消耗的時間與硬件性能成反比，圖形硬件性能越強，耗時越短。其中，GTX 780Ti性能最強，耗時最短。

3.2 MCUPS值評測分析

第一，NVIDIA Geforce GTX 760、NVIDIA Geforce GTX 770、NVIDIA Geforce GTX 780Ti和 SSEARCH 相比，平均加速分別10.4 倍、12 倍、14.5 倍，比對速度都有飛躍性地提高。

第二，隨著查詢序列長度的增加，比對加速越發顯著，查詢序列長度越大、比對需要耗費的時間越多，圖形硬件的加速效率則越高。

第三，從760到780Ti，隨著性能不斷增強，比對速度提升顯著。以此次測試中性能最高的780Ti為例，對于不同長度的查詢序列，對比于SSEARCH，其平均加速越為14.5倍;當查詢序列長度增加到最高的4113aa時，比對加速達到了22.9倍。

3.3 綜合分析

從上述時間性能和MCUPS性能值均可知，利用圖形硬件進行加速的比對算法，其速度要遠超利用啟發式算法提速的比對算法，能極大地提升生物序列局部比對速度。此外，近年來圖形硬件每月都有新款推出，硬件性能不斷提升，通過圖形硬件進行比對加速還有更為廣闊的提升空間及更為光明的前景。第三，相比利用啟發式算法進行加速，圖形硬件加速算法嚴格執行動態規劃局部比對算法，比對準確性不會有任何地犧牲。最后，圖形硬件可在個人電腦上運行，價格非常便宜，性價比極高，相比于專門的生物信息處理機器和專有硬件來說，成本可以忽略不計，能大規模推廣和應用，通用性極強。

4 結論

綜上所述，將準確率最高的動態規劃局部比對算法予以實現，然后將其映射到圖形硬件(計算機顯卡)上實現算法加速，并通過實例評測其性能。結果表明該加速算法在保證比對準確性的同時，能顯著提高比對速度。從而為實現準確性、高效性、通用性并存的生物序列局部比對方法做出新的探索，為生物學、醫學、高性能計算的結合開拓新的方向。

References)

［1］ 陳銘，包家立.生物信息學［M］.北京:科學出版社，2013.CHEN Ming，BAO Jiali.Bioinformatics［M］.Beijing:Science Press，2013.

［2］ WIKIPEDIA S.Bioinformatics:proteomics，hidden markov model，biostatistics，proteome，sequence alignment，full genome sequencing［M］.USA:University-Press，2013.

［3］ SMITH T，WATERMAN M.Identification of common molecular subsequences［J］.Journal of Molecular Biology，1981，147(1):195 －197.

［4］ David J.Multiple sequence alignment methods［M］.New York:Springer，2014.

［5］ TANG M，TONG R F，NARAIN R，et al.GPU-based streaming algorithm for high-resolution cloth simulation［J］.Computer Graphics Forum，2013，32(7):21 －30.

［6］ MANOCHA D.General-purpose computations using graphics processors［J］.Computer，2005，38(8):85 －88.

［7］ 張寶華，王海水，許祿.DNA序列編碼及相似度計算［J］.高等學校化學學報，2006，27(12):2277－2280.ZHANG Baohua，WANG Haishui，XU lu.DNA sequence code and similarity computations［J］.Chemical Journal of Chinese Universities，2006，27(12):2277 －2280.

［8］ WEISS，R，SHRAGGE J.Solving 3D anisotropic elastic wave equations on parallel GPU devices［J］.Geophysics，2013，78(2):7 －15.

［9］ BAI H T，LI Y G，CHEN L Y，et al.Parallel optimization of geometric correction algorithm based on CPU-GPU hybrid architecture［J］.Applied Mechanics and Materials，2014，543(26):2804 －2808.

［10］ DANIE L，JESWIN G，JUSTIN H，et al.Stencil-Aware GPU optimization of iterative solvers［J］.SIAM Journal on Scientific Computing，2013，35(5):209 －228..

［11］ PAYNE B R，OWEN G S，WEBER I.Accelerating protein structure recovery using graphics processing units［A］.In Proceedings of International Conference on Computational Science(ICCS 05)［C］.Singapore:Springer-Verlag，2005.

［12］ MENG X D，JI Y Q.Modern computational techniques for the HMMER sequence analysis ［J］.Bioinformatics，2013，2013(03):13 －17.