基于AVX指令集BWT算法在DNA序列比對中應用

孫紅敏，杜博雅，鄭萍，李東野，曹延杰，侯星辰

（1.東北農(nóng)業(yè)大學電氣與信息學院，哈爾濱 150030；2.武漢理工大學計算機科學與技術學院，武漢 430070）

基于AVX指令集BWT算法在DNA序列比對中應用

孫紅敏1，杜博雅1，鄭萍1，李東野2，曹延杰1，侯星辰1

（1.東北農(nóng)業(yè)大學電氣與信息學院，哈爾濱 150030；2.武漢理工大學計算機科學與技術學院，武漢 430070）

新一代高通量測序技術發(fā)展產(chǎn)生大規(guī)模DNA序列片段，快速準確地將短序列比對到參考基因組成為生物信息學重要研究課題之一。針對BWT索引技術序列比對算法研究，提出基于Intel微架構AVX指令集優(yōu)化BWT算法，通過改進計算方式實現(xiàn)算法并優(yōu)化。結果表明，應用AVX指令集可減少CPU訪存次數(shù)，降低算法時間復雜度，提高序列比對效率，為基因數(shù)據(jù)分析提供更高效快速序列比對方法，加快對全基因組序列處理。

序列比對；AVX指令集；BWT算法；并行優(yōu)化

隨著新一代測序平臺（NGS）廣泛使用，測序通量高、時間和成本下降使DNA序列數(shù)據(jù)顯著增加，產(chǎn)生TB級測序片段。快速準確地將大規(guī)模DNA序列比對到參考基因組上成為生物信息學亟待解決難題[1]。將短序列與參考基因組序列比對并在參考基因組上定位，根據(jù)比對結果分析預測序列間相似性，可為后續(xù)開展表達量分析、預測SNP位點、選擇性剪接分析、疾病相關性、藥物研發(fā)、鑒定基因組中功能基因等研究奠定基礎[2]。DNA序列比對現(xiàn)已成為生物信息學重要內容之一[3]。

在常用比對算法中，根據(jù)索引結構不同，采用空位種子片段索引算法Maq時間效率較低，序列中插入缺失（Indel）處理困難[4]；BFAST序列比對算法速度慢[5]?；诠１肀葘λ惴≒ASS[6]和SOAP[7]查找速度快，但內存消耗大；而采用Burrows-Wheeler轉換基于后綴、前綴索引算法Bowtie[8]、BWA[9]、BWASW等即使采用BWT技術，可解決內存消耗問題，但訪存次數(shù)多、計算量大，在數(shù)據(jù)規(guī)模較大時查找時間長。因此如何減少時間消耗成為優(yōu)化改進BWT方法研究核心。本文針對BWT索引技術耗時問題，優(yōu)化多次遞歸調用函數(shù)部分，引入AVX指令技術方法解決運算時耗問題，提高BWT算法序列比對速率，實現(xiàn)大量基因數(shù)據(jù)快速處理。

1 BWT算法

BWT算法通過對字符串循環(huán)移位得到字符矩陣排序和變換，將基因組序列壓縮并建立索引，通過查找和回溯定位字符串DNA序列比對。傳統(tǒng)BWT構造和查找具體過程如下：假定參數(shù)T為字符串‘ATGGACT’，定義$字典序小于其他字符。在輪轉矩陣（Burrows-Wheeler matrix）中，最后一列第i個出現(xiàn)字符x，與第一列中第i個出現(xiàn)x是同一個，這是BWT算法LF Mapping性質，通過最后一列BWT（T）還原字符串T。BWT循環(huán)排序過程如圖1所示。

圖1 BWT循環(huán)排序過程Fig.1Rotation sorting process of BWT

原串‘ATGGACT$’經(jīng)循環(huán)排序得到最后一列BWT（T）=TG$AGTCA。生成后綴數(shù)組（Suffix Array）SA[r]={7，4，0，5，3，2，6，1}。BWT中定義兩個參數(shù)C（c）和Occ（c，r），c代表堿基字符‘ATGC’，r代表矩陣行。C（c）表示字典序小于字符c所有字符個數(shù)，Occ（c，r）函數(shù)表示在BWT（T）中第r行之前出現(xiàn)字符c個數(shù)。Occ（r，c）固定不變，可將其保存在內存中供查找。LF（r）Matching算法過程：

在LF mapping基礎上構建索引，BWT字符串匹配。如圖1所示第一列BWT（T）前添加數(shù)組SA[r]可以記錄第r行在參考基因組中位置。利用BWT向后搜索（Backward）算法模擬前綴樹自頂向下遍歷實現(xiàn)字符串精確匹配。前綴樹每條邊表示一個字符，從樹根結點到葉子節(jié)點構成參數(shù)T前綴。結點數(shù)字對應以此結點到根結點之間字符作后綴在右邊排序后綴數(shù)組中范圍，定義兩個公式計算以參數(shù)T作為前綴后綴數(shù)組范圍：

字符串匹配即為查找SA區(qū)間坐標位置，即求解后綴數(shù)組值。字符串‘ATGGACT’前綴樹與反向文本‘TCAGGTA’后綴樹相同，圖2所示為反向文本T后綴樹，其中符號^表示T起始位置。

圖2 字符串‘TCAGGTA’后綴樹Fig.2Suffix tree of string‘TCAGGTA’

圖3 [-R（W），-R（W）]字符串精確匹配步驟Fig.3[-R（W）,-R（W）]String matching procedure

綜上可知，LF（r）Matching算法和精確匹配Backward算法反復調用Occ（r，c），而Occ（r，c）函數(shù)計算密集、訪問頻繁、數(shù)據(jù)量大、計算過程時間長。應用AVX指令集改進該部分運算方式，對bwt_occ函數(shù)內部加法運算并改寫，加快算法運算速度，提升算法效率。

2 AVX指令并優(yōu)化BWT過程

全新指令擴展集AVX（Intel advanced vector extensions）是基于單指令多數(shù)據(jù)集（SIMD）技術思想設計實現(xiàn)[10]，即一個CPU指令執(zhí)行周期內同一指令下同時處理內存中多道數(shù)據(jù)，實現(xiàn)多數(shù)據(jù)處理，提高數(shù)據(jù)處理效率。AVX指令支持256位矢量運算，可以一次處理8個單精度浮點型數(shù)據(jù)，使浮點數(shù)運算速度提高近8倍。相比于傳統(tǒng)指令，AVX指令節(jié)省訪問時間，提高數(shù)據(jù)處理速度。

2.1 函數(shù)耗時統(tǒng)計

通過熱點分析工具perf對程序分析，通過BWT算法中序列比對過程測試各函數(shù)耗時比例結果見圖4。

圖4 Occ函數(shù)中各模塊耗時情況Fig.4Time consuming of each module in Occ function

由圖4可知，bwt_occ函數(shù)運算時間占比42.78%，耗時最多。為提高算法執(zhí)行效率，采用結合AVX指令集求和計算思想，即每次計算8個數(shù)據(jù)，對bwt_occ函數(shù)運算部分優(yōu)化。

2.2 AVX優(yōu)化設計流程

256位AVX寄存器對應_m256數(shù)據(jù)類型，即256位單精度浮點型數(shù)據(jù)，該寄存器可一次并處理8個單精度浮點數(shù)或4個雙精度浮點數(shù)。在實際編碼過程中，occ每隔32行存儲一次，將此Occ（c，r）值記為一個“checkpoint”，這樣存儲結構能減少內存開銷。計算occ時需通過遍歷BWT（T）字符串得到結果，直至得到BWT矩陣行數(shù)為32倍數(shù)時取出相應checkpoint值。在求解checkpoint過程中，需要判斷并計算每次遍歷位置值。將這些數(shù)值存儲到數(shù)組avx_result中，當遍歷到checkpoint后，使用AVX指令集計算數(shù)據(jù)，可充分發(fā)揮AVX指令集高速性，實現(xiàn)一個CPU指令執(zhí)行周期內同一指令下同時處理多道數(shù)據(jù)。多數(shù)據(jù)處理，減少occ計算次數(shù)，減少程序運行時間。

AVX指令優(yōu)化分四部分：①變量定義與初始化。②AVX批量處理。③合并，將__m256上多個數(shù)值合并到求和變量。④處理剩余數(shù)據(jù)，對剩余數(shù)據(jù)采用基本逐項相加。算法優(yōu)化流程如圖5所示。

代碼優(yōu)化具體實現(xiàn)如下：

①定義一個浮點型數(shù)組存放整型數(shù)據(jù)_occ_aux值，對這些值求和運算。

②使用AVX數(shù)據(jù)類型存放浮點型數(shù)組地址。AVX指令集在通過load指令加載數(shù)據(jù)時，需指針加減操作，通過Load指令加載pp指針，將指向存放_occ_aux值數(shù)組指針pp加載到y(tǒng)fsLoad中。AVX指令為_mm256_load_ps（）；表示對ymm寄存器中內容單精度浮點數(shù)加法運算。

圖5 AVX改進算法流程Fig.5Flowchart of improved algorithm by AVX

③通過add指令求和：使用AVX指令_mm256_ setzero_ps（）；將yfsSum數(shù)值置零。其中yfsLoad中存放occ值，yfsSum中存放yfsLoad中值和。AVX求和運算指令為_mm256_add_ps（）。

④將每一步驟計算結果存放在浮點型變量中，通過強制類型轉換將該浮點型變量值賦值回原始int型變量。

AVX過程共用到三個函數(shù)：

①_mm256_setzero_ps：將__m256中每一個單精度浮點數(shù)均賦值為零。偽代碼：for（i=0；i＜8；++ i）C[i]=0.0f；

②_mm256_load_ps：從內存中對齊加載8個單精度浮點數(shù)到__m256變量。偽代碼：for（i=0；i＜8；++i）C[i]=_A[i]；

③_mm256_add_ps：相加指令，對2個__m256變量8個單精度浮點數(shù)垂直相加。偽代碼：for（i=0；i＜8；++i）C[i]=A[i]+B[i]。

將原有bwt_int數(shù)據(jù)類型轉換為AVX數(shù)據(jù)類型_m256并行運算。加法運算通過AVX指令集load加載指令和add求和指令實現(xiàn)并行運算優(yōu)化。計算bwt_occAVX指令集實現(xiàn)具體操作為：

yfsSum=_mm256_setzero_ps（）;

yfsLoad=_mm256_load_ps（pp）;

yfsSum=_mm256_add_ps（yfsSum，yfsLoad）;

AVX指令加法運算工作模式如圖6所示。

圖6 單精度浮點型AVX加法工作模式Fig.6Single-precision floating-point addition AVX mode

通過整合上述參數(shù)，就可得到BWT算法計算occ值AVX編程，算法具體實現(xiàn)如下：

Begin

設置avx塊寬，賦值broadwidth

計算塊數(shù)，賦值cntBlock

計算剩余數(shù)量，賦值cntRem

計算所需內存空間，賦值aux_result

for p到end

把__occ_aux4（bwt，*p）保存到aux_result

end for

for i=0到cntBlock

使用avx指令計算

end for

將avx指令處理后數(shù)據(jù)合并

for i=0到cntRem

處理剩余數(shù)據(jù)

end for

應用AVX指令BWT算法，減少主循環(huán)次數(shù)和occ計算次數(shù)，提高occ函數(shù)運算速率，使數(shù)據(jù)處理速度顯著提升。理論上能使算法復雜度降低至0（n/8），其中n為遍歷次數(shù)。但在實際代碼實現(xiàn)中，由于每次均需初始化AVX指令并保存遍歷數(shù)據(jù)，因此實際算法時間復雜度降低至0（n/2）。改進后算法性能穩(wěn)定，在使用不同數(shù)據(jù)規(guī)模測試時，不會使時間復雜度隨數(shù)據(jù)規(guī)模變化而改變，時間復雜度穩(wěn)定趨近于0（n/2）。

3 結果與分析

測試平臺為：Intel（R）Core（TM）i7-3770 CPU@ 3.40GHZ，4核，內存8 G，Ubuntu13.04 OS。試驗使用BWA（Burrows-Wheeler aligner）軟件對BWT算法優(yōu)化對比驗證。使用BWA需要兩種輸入文件，參考基因組數(shù)據(jù)Reference genome data（fasta格式為. fa、.fasta、.fna）和短序列數(shù)據(jù)Short reads data（fastaq格式為fastaq、.fq）。

為確定方法有效性，采用真實數(shù)據(jù)集測試。試驗中采用Illumina/Solexa測序平臺大豆基因數(shù)據(jù)——ensemble數(shù)據(jù)庫參考基因組序列文件ftp.ensemblgenomes.org，該物種測序序列為Glycine max.dna. genome.Fa.gz，短序列數(shù)據(jù)為NCBI數(shù)據(jù)庫下載大豆高通量數(shù)據(jù)sra格式文件。通過sratoolkit工具轉換為序列比對可識別單末端測序序列數(shù)據(jù)（Single-end）fastq文件。序列比對到基因組后通常采用SAM（SequenceAlignment/Map）格式存儲。

試驗中分別截取1 G fastq序列比對文件中長度為200 bp 200～1 000 M不同大小序列片段比對試驗。BWT算法與基于AVX改進算法時間對比如表1所示。

表1 BWT算法與AVX-BWT算法時間對比Table 1Time comparison between BWT algorithm and AVX-BWT

在同一序列規(guī)模情況下，BWT與AVX代碼運行時間與序列大小間關系如圖7所示。

圖7 BWT算法優(yōu)化運算效率Fig.7BWT algorithm optimization algorithm efficiency

如圖7所示，當fq序列大小為200 Mb時，BWT源碼運行時間為4 246.14 s，AVX優(yōu)化后運行時間縮減為2 533.39 s；當fq數(shù)據(jù)大小為1 000 Mb時，比對時間由21 331.92 s縮減至12 604.30 s。設BWT時間曲線斜率為k1，AVX時間曲線斜率為k2，兩條曲線分別為y1=k1x+b1，y2=k2x+b2。BWT方法曲線函數(shù)為y1=4284.9x+264.42，AVX方法曲線函數(shù)為y2= 2520.2x-85.12。

如圖7可知，隨序列規(guī)模增加，兩種方法運行時間均呈線性增加，逐漸趨近于0（n/2），該試驗結果與給出時間復雜度理論分析結果接近，即兩種方法時間復雜度與序列長度呈線性相關，AVX優(yōu)化后時間復雜度趨近于0（n/2）。得出AVX指令優(yōu)化后，耗時減少一半，序列比對時間明顯減少。

4 討論

對于BWT算法改進，采用AVX指令技術實現(xiàn)單線程內加速，而Intel采用Cilk工具對BWT算法存在負載失衡從多線程角度算法加速，但存在硬件條件限制，是否適于處理數(shù)據(jù)需深入研究。二階BWT索引方法[11]與傳統(tǒng)BWT方法中逐位索引查找不同，從算法結構角度優(yōu)化，改進后BWT方法按雙位索引查找，序列比對速度提高10%，而使用AVX指令從計算方式改進，優(yōu)化提高40%～50%，速度提升優(yōu)于二階索引方法。

分別對C代碼BWT源程序和AVX指令優(yōu)化代碼測試比較，驗證優(yōu)化后函數(shù)功能正確性。序列比對過程中AVX指令改寫后程序將occ值輸出時，同樣需要將數(shù)值類型由整型轉換為浮點型輸出，否則無法將序列比對到SA后綴數(shù)組位置坐標輸出為SAM格式。將輸出文件與源碼輸出文件通過文本比對工具比較，2 000 000行數(shù)據(jù)存在約200行差異序列，即結果存在0.01%誤差。

出現(xiàn)誤差原因是AVX處理浮點型數(shù)據(jù)速度快，因此在優(yōu)化過程中改寫循環(huán)部分使用float型加法運算，而BWT源碼中使用int64數(shù)據(jù)類型，整型與浮點型直接作加法操作會出現(xiàn)錯誤。因此在AVX加法操作結束后需將浮點型累加并轉換為整型后與BWT中int64相加。float在計算機中二進制存儲方式遵從IEEE規(guī)范，存儲中分符號位（Sign）、指數(shù)位（Exponent）、尾數(shù)部分（Mantissa）。由于DNA序列中不存在負數(shù)情況，因此只需截取指數(shù)位和尾數(shù)位。將內存存儲float二進制格式轉化為整型，在數(shù)據(jù)類型轉換時當數(shù)值超過float可存儲32位時出現(xiàn)溢出現(xiàn)象，導致結果差異顯著。對于大數(shù)據(jù)溢出現(xiàn)象而產(chǎn)生誤差可采取分解存儲解決方案，但移位和拆分操作耗費時間，使算法整體執(zhí)行時間消耗未降低?？紤]較大數(shù)據(jù)出現(xiàn)概率較小，故使用存在0.01%誤差float數(shù)據(jù)類型AVX代碼優(yōu)化。

改寫B(tài)WT中耗時比重較大、計算密集、訪問頻繁的bwt_occ函數(shù)AVX指令，在計算occ時將AVX指令應用在循環(huán)迭代中優(yōu)化加法運算，減少函數(shù)計算次數(shù)，減少訪問保存次數(shù)，使一個線程中一條指令并行做8次運算，提升算法執(zhí)行效率。利用AVX指令集優(yōu)化occ函數(shù)指令，可實現(xiàn)數(shù)據(jù)同步處理。在相同數(shù)據(jù)條件下，對比BWT源碼，AVX優(yōu)化后編碼速率在原算法基礎上平均提高2倍，AVX技術使BWT算法執(zhí)行效率得到顯著提升。

5 結論

本文在DNA序列比對算法BWT基礎上，提出一種應用AVX指令技術改變數(shù)據(jù)運算方式方法。在相同數(shù)據(jù)結構下，應用AVX指令能夠減少序列比對索引算法中函數(shù)內部循環(huán)次數(shù)和計算次數(shù)，使算法時間性能顯著提升，提高序列比對過程查找效率。本文針對BWT算法中字符串查找部分指令優(yōu)化，運算效率提升有限，后續(xù)研究可探討構建索引運算部分優(yōu)化及Occ函數(shù)其他模塊AVX指令改寫，進一步提高DNA序列比對速度。應用算法并行化將成為序列分析重要發(fā)展方向，AVX指令集可應用于數(shù)據(jù)壓縮及數(shù)據(jù)分類等其他基因序列處理和研究領域，提高基因序列數(shù)據(jù)挖掘及全基因組序列處理效率。

[1]Kahn S D.On the future of genomic data[J].Science,2011,331 (6018):728-729．

[2]李麗文,朱延明,李杰.RNAi技術在植物功能基因組學中的研究進展[J].東北農(nóng)業(yè)大學學報,2007,38(1):119-124.

[3]國宏哲,王亞東.基因組Mapping系統(tǒng)索引構建原理[J].智能計算機與應用,2012,47(4):3.

[4]Li H,Ruan J,Durbin R.Mapping short DNA sequencing reads and calling variants using mapping quality scores[J].Genome Res, 2008,18(11):1851-1858

[5]Homer N,Merriman B,Nelson S F.BFAST:an alignment tool for large scale genome resequencing[J].PLoS One,2009,4(11):7767. [6]Campagna D,Albiero A,Bilardi A,et al.PASS:a program to align short sequences[J].Bio-informatics,2009,25(7):967-968.

[7]Li R Q,Li Y R,Kristiansen K,et al.SOAP:Short Oligonucleotide Alignmen Program[J].Bio-informatics,2008,24(5):713-714．

[8]Langmead B,Trapnell C,Pop M,et al.Ultrafast and memoryefficient alignment of short DNA sequences to the human genome [J].Genome Biol,2009,10(3):1-10.

[9]Li H,Durbin R.Fast and accurate short read alignment with Burrows-Wheeler transform[J].Bio-informatics,2009,25(14):1754-1760.

[10]朱林,馮燕.基于單指令多數(shù)據(jù)技術的H.264編碼優(yōu)化[J].計算機應用,2005,25(12):2799-2802.

[11]趙雅男,徐云,程昊宇.序列比對算法中的BW變換索引技術研

BWT algorithm based on AVX instructions in the application of theDNA sequence alignment/

SUN Hongmin1,DU Boya1,ZHENG Ping1,LI Dongye2,CAO Yanjie1,HOU Xingchen1(1.School of Electrical and Information,Northeast Agricultural University, Harbin 150030,China;2.School of Computer Science and Technology,Wuhan University of Technology,Wuhan 430070,China)

The development of new generation of high throughput sequencing technology had produced large-scale DNA sequence fragments,how to quickly and accurately match to the reference genome had become one of the important research topics in bioinformatics.According to the sequence alignment algorithm of BWT index,new BWT algorithm based on Intel micro architecture was proposed to optimize the AVX instruction set,improved computation method to achieve parallel optimization algorithm.The results showed that the application of AVX instruction set CPU could reduce the number of memory accesses and reduced the time complexity of the algorithm,improved the efficiency of gene sequence alignment,data analysis to provide more efficient fast sequence alignment methods,to speed up the processing and analysis of the whole genome sequence.

sequence alignment;AVX instruction;BWT index;parallel optimization

TP399

A

1005-9369（2016）11-0093-07

時間2016-12-1 14:48:43[URL]http://www.cnki.net/kcms/detail/23.1391.S.20161201.1448.022.html

孫紅敏,杜博雅,鄭萍,等.基于AVX指令集BWT算法在DNA序列比對中應用[J].東北農(nóng)業(yè)大學學報,2016,47(11)：93-99.

Sun Hongmin,Du Boya,Zheng Ping,et al.BWT algorithm based on AVX instructions in the application of the DNA sequence alignment[J].Journal of Northeast Agricultural University,2016,47(11)：93-99.(in Chinese with English abstract)

2016-08-14

國家“863計劃”項目（2013AA10230304）

孫紅敏（1971-），女，教授，碩士生導師，研究方向為生物信息技術。E-mail:sunhongmin111@126.com