崔英博，黃春，唐滔，楊燦群，廖湘科，彭紹亮

1. 國防科技大學計算機學院，湖南 長沙 410073；2. 湖南大學信息科學與工程學院，湖南 長沙 410082；3. 國家超級計算長沙中心，湖南 長沙 410082

1 引言

近年來，高通量測序技術(shù)的迅猛發(fā)展給生命科學帶來了巨大變革，測序通量不斷提高，測序成本不斷下降，如單個人類基因組的測序成本逐年降低，基因組數(shù)據(jù)規(guī)模不斷增長。基因組數(shù)據(jù)大約每7個月就會增加一倍[1]?，F(xiàn)有服務(wù)器以及序列分析算法已經(jīng)無法及時有效地處理如此大規(guī)模的數(shù)據(jù)，需要借助并行計算機的強大算力以及并行算法的有效支撐來實現(xiàn)對基因組大數(shù)據(jù)的有效處理[2]。

測序儀產(chǎn)生讀段（reads）后，首先進行質(zhì)量控制，去除測序質(zhì)量較差的數(shù)據(jù)；然后進行序列比對，將讀段回帖到基因組，找到讀段最可能的起源位置，并輸出比對文件；再基于比對文件檢測基因組的變異情況，主要包括單核苷酸多態(tài)性（single nucleotide polymorphism，SNP）、插入刪除變異（indel）、結(jié)構(gòu)變異（structure variation，SV）和拷貝數(shù)變異（copy number variation，CNV）等；最后根據(jù)需要，進行特定的功能分析。序列比對和變異檢測是基因組數(shù)據(jù)分析的基礎(chǔ)環(huán)節(jié)，是后續(xù)功能性分析的基礎(chǔ)，也是數(shù)據(jù)分析流程中最耗時的步驟[3]。

為有效處理高通量測序技術(shù)帶來的海量基因組大數(shù)據(jù)，本文選取基因組變異檢測中常用的序列比對和SNP檢測兩個步驟，對相關(guān)算法進行了改進，并利用OpenMP（open multi-processing）和消息傳遞接口（message passing interface，MPI）等并行技術(shù)實現(xiàn)了多級并行。在不同數(shù)據(jù)集和并行規(guī)模下的測試中，核心算法的加速比達到9倍以上，大規(guī)模測試中算法的并行效率保持在60%以上，在保證精度的前提下獲得了良好的并行性能和可擴展性，有效提高了基因組大數(shù)據(jù)變異檢測的能力。

2 相關(guān)工作

2.1 序列比對算法

序列比對算法以讀段和參考基因組為輸入，將讀段比對到基因組，找到讀段最可能的起源位置?；蚪M規(guī)模一般較大，且高通量測序產(chǎn)生的讀段數(shù)量龐大，現(xiàn)有的比對工具大部分是借助索引結(jié)構(gòu)來處理大規(guī)模數(shù)據(jù)的。根據(jù)索引結(jié)構(gòu)的不同，可以將序列比對算法分為兩類：基于哈希表的比對算法和基于后綴樹的比對算法[4]。

最早的基于哈希索引的比對工具是BLAST[5]，后續(xù)出現(xiàn)的基于哈希索引的比對工具有SOAP[6]、SeqMap[7]、RMAP[8]、BFAST[9]、CloudBurst[10]、PerM[11]和GNUMAP[12]等。

后綴樹是一種存儲著一個字符串所有后綴的數(shù)據(jù)結(jié)構(gòu)，在后綴樹中，所有相同的后綴都會合并成一條路徑。因此，基于后綴樹的序列比對可以合并對相同后綴的比較，提高效率?；诤缶Y樹的比對算法查詢時間為O(q)，其中q為子串的長度，時間復雜度與母串無關(guān)，具有較高的查詢效率。但是，后綴樹的空間復雜度較高，內(nèi)存需求很大，無法在常用的服務(wù)器上實現(xiàn)索引構(gòu)建，因此基于后綴樹的比對算法在早期應(yīng)用較少。

FM索引[13]的出現(xiàn)大大降低了后綴樹的空間開銷，使得基于后綴樹的比對算法得到了廣泛應(yīng)用。FM索引可以被理解為壓縮的后綴樹，其能夠在壓縮狀態(tài)下實現(xiàn)字符串搜索，因此極大地降低了索引的空間開銷，如人類基因組的FM索引大小為3 GB左右，大多數(shù)服務(wù)器甚至臺式機能滿足需求。基于FM索引的比對算法具有時間和空間的雙重優(yōu)勢，迅速在二代序列比對領(lǐng)域流行起來。常用的基于FM索引的比對工具有BWA[14]、Bowtie[15]和SOAP2[16]等，本文的工作基于BWA開展。

2.2 SNP檢測算法

SNP是基因組中單個堿基位點的變異。序列比對完成后，基因組中大多數(shù)位點會對應(yīng)多條讀段，SNP檢測通過逐個位點比較參考基因組和輸入序列來確定SNP位點。每個位點的SNP檢測過程完全相同，且相互獨立。SNP檢測算法一般會綜合位點的堿基類型、測序質(zhì)量、在基因組中的位置等信息，計算當前位點是SNP位點的概率。基因組規(guī)模一般較大，而內(nèi)存空間有限，無法一次性完成全基因組的SNP檢測，因此SNP檢測算法一般會將基因組分為若干窗口，對窗口進行逐個檢測。

SNP檢測算法將比對完成的讀段、參考序列作為輸入，有時也將已知SNP的數(shù)據(jù)庫信息作為輸入[17]。研究人員開發(fā)了一些基于Web的SNP檢測工具，只要通過瀏覽器上傳相關(guān)數(shù)據(jù)，即可完成分析，包括HaploSNPer[18]和SNiPlay[19]等。但是基于Web的分析工具在數(shù)據(jù)安全性方面存在一定的風險，因此，更常用的工具是能夠獨立運行的SNP檢測工具，包括QualitySNP[20]、SAMtools[21]、SOAPsnp[17]、GATK[22]、Isaac[23]、SNVer[24]、VarScan[25]和VarScan2[26]等，本文的工作基于SOAPsnp開展。

SNP分析是一個相對耗時的過程，因此，也存在不少對SNP算法的優(yōu)化，如Crossbow[27]利用Hadoop和云計算加速SNP檢測，Rainbow[28]針對大數(shù)據(jù)集對Crossbow做了進一步優(yōu)化，GSNP[29]通過GPU來加速SNP的計算過程。

2.3 OpenMP和MPI技術(shù)

如圖1、圖2所示，從硬件實現(xiàn)的角度，可以將并行計算機系統(tǒng)分為共享內(nèi)存系統(tǒng)和分布式內(nèi)存系統(tǒng)兩種。在共享內(nèi)存系統(tǒng)中，多個計算核心或CPU共享內(nèi)存，可以提高內(nèi)存數(shù)據(jù)的利用率；在分布式內(nèi)存系統(tǒng)中，每個節(jié)點擁有自己獨立的CPU和內(nèi)存，節(jié)點間通過互聯(lián)網(wǎng)絡(luò)連接，CPU只能訪問自己節(jié)點上的內(nèi)存空間，無法訪問其他節(jié)點上的內(nèi)存空間。

針對兩種不同的并行系統(tǒng)，存在兩種不同的并行編程模型。面向共享內(nèi)存系統(tǒng)的并行編程模型以多線程為主，主要包括OpenMP和pthread，其中pthread需要程序員控制多線程的發(fā)起、數(shù)據(jù)分配、同步等操作，較為煩瑣。而OpenMP僅需要在程序中可并行的模塊前添加指導語句即可，由編譯器和運行時庫負責對數(shù)據(jù)進行劃分、發(fā)起多線程并行、管理私有和共享數(shù)據(jù)、進行線程同步等。相比之下，OpenMP更易于實現(xiàn)，而且具有良好的可移植性[30-31]，因此本文主要基于OpenMP實現(xiàn)多線程并行。

面向分布式內(nèi)存系統(tǒng)的并行編程模型以MPI為主，支持多進程并行，進程間通過傳遞消息的方式實現(xiàn)跨節(jié)點通信。MPI并不是一種新的編程語言，而是一個通信庫，支持C、C++和Fortran等語言，可以看作對這些語言的擴展。MPI除了支持點對點通信外，還支持廣播和規(guī)約等聚合通信，非常便利，是分布式系統(tǒng)并行編程的事實標準。

3 相關(guān)算法分析

3.1 序列比對算法分析

序列比對算法一般包括讀入讀段、比對和寫出結(jié)果3個步驟。為提高比對效率，一般會事先對參考基因組建立索引。對于給定物種，參考基因組的索引生成后可以被重復使用，對整體比對時間的影響并不大，因此，本文主要對比對過程進行并行優(yōu)化。

BWA是使用非常廣泛的、基于Burrows-Wheeler變換（Burrows-Wheeler transform，BWT）的短序列比對工具之一[32]。BWT是一種數(shù)據(jù)壓縮算法，其主要思想是對一個給定文本的所有輪轉(zhuǎn)進行排序，返回的排序的最后一列被稱為BWT字符串，該列為BWT的結(jié)果。文本經(jīng)過BWT后會將很多相同的字符排序到一起，因此更加容易進行壓縮。FM索引是一種基于BWT的索引結(jié)構(gòu)，支持在壓縮狀態(tài)下對文本進行字符串查詢。在BWA中，精確匹配采用后向搜索方法實現(xiàn)，即對于一個字符串，從后向前逐個字符地進行比對，后向搜索實際上完成了對基因組構(gòu)成的后綴樹的自頂向下遍歷。

3.2 SNP檢測算法分析

圖3為SOAPsnp檢測SNP的流程。SOAPsnp的輸入包括比對的讀段和參考基因組，有時也包括已知SNP的數(shù)據(jù)庫信息，SOAPsnp的輸出為保守基因型信息。SOAPsnp主要包括7個模塊：cal_p_mat、read_site、counting、likelihood、posterior、output和recycle，核心數(shù)據(jù)結(jié)構(gòu)包括p_matrix、base_info和type_likely。

cal_p_mat模塊掃描讀段比對文件，計算滿足四元組合(堿基質(zhì)量分數(shù), 讀段中的偏移, 參考基因組堿基類型, 讀段中的堿基類型)的堿基的數(shù)量，并生成校準矩陣p_matrix，用于調(diào)整概率計算中的測序質(zhì)量分數(shù)。其他6個模塊通過多輪處理，完成對全基因組的SNP檢測，其中每輪會處理一個窗口：首先，read_site模塊從輸入文件中加載固定數(shù)量的位點（一個窗口寬度）；然后，counting模塊收集堿基信息，并保存到base_info內(nèi)，在每個窗口中，逐個位點依次進行SNP檢測；likelihood模塊以base_info和p_matrix為輸入，計算每個位點每種基因型為SNP的概率，并輸出概率值type_likely；posterior模塊根據(jù)貝葉斯模型計算每種基因型的后驗概率，這里將基因型概率以及讀段和參考基因組之間的SNP估計值作為先驗概率，后驗概率計算完成后，當前位點的檢測結(jié)果被輸出到文件，然后從likelihood模塊開始，處理當前窗口的下一位點；recycle模塊負責窗口切換，即處理跨窗口的讀段，并為下一窗口初始化緩沖區(qū)。

p_matrix是一個浮點類型的四維矩陣，規(guī)模為256×256×4×4，大小為8 MB。矩陣的4個維度分別對應(yīng)堿基質(zhì)量分數(shù)、讀段中的偏移位置、參考基因組堿基類型以及讀段中的堿基類型。p_matrix矩陣被用于在似然計算中校正測序的質(zhì)量分數(shù)。base_info是一個規(guī)模為4×64×256×2的矩陣，存儲了堿基類型、測序質(zhì)量分數(shù)、讀段中的偏移坐標和讀段的方向等信息。base_info是一個計數(shù)矩陣，矩陣中的一個位置對應(yīng)一種四元組合，矩陣中記錄的數(shù)值就是滿足相應(yīng)四元組合的堿基的數(shù)量。

4 并行優(yōu)化方法

4.1 序列比對算法的并行優(yōu)化

利用BWA進行序列比對主要包括生成參考基因組索引和序列比對兩個步驟。其中，對于給定的基因組，參考基因組的索引僅需在比對前生成一次，可以被重復使用，因此，研究生成索引過程的并行化對于序列比對來說意義不大，其起到的加速效果十分有限。本文主要研究的內(nèi)容是將序列比對步驟并行化。

序列比對步驟具體包括讀入?yún)⒖蓟蚪M索引、讀入短讀序列、序列比對和輸出比對結(jié)果4個子步驟，下面分別介紹各個步驟的并行化方法。

（1）讀入?yún)⒖蓟蚪M索引

在進行比對之前，首先必須獲得短讀序列和參考基因組。BWA雖然實現(xiàn)了多線程比對，但是讀入?yún)⒖夹蛄幸廊皇菃尉€程實現(xiàn)，而且BWA讀入索引的方法相對簡單。在進行多線程比對時，每個線程會從主線程獲得一個指向索引數(shù)據(jù)結(jié)構(gòu)的指針。這種方法雖然簡單，但是行之有效，因為僅需幾秒鐘就可以將人類基因組索引讀入內(nèi)存中。

多線程可以共享內(nèi)存空間，但是對于跨節(jié)點的MPI程序來說，各個進程的內(nèi)存空間相互獨立，如果仍使用BWA原本的方法，那么需要每個進程分別讀一次索引文件，效率低下，而且，當進程規(guī)模較大時，會對文件系統(tǒng)造成巨大的壓力，嚴重影響程序性能。

基于以上考慮，本文使用主從模式讀取索引。首先指定一個進程為主進程，由主進程將參考基因組索引讀入內(nèi)存，然后再用MPI函數(shù)將索引廣播到所有其他的進程，這樣每個進程都有了索引的一個副本。

這種方法的好處是只需要一個進程訪問索引文件即可，缺點是在主進程讀入索引時，其他進程都處于空閑狀態(tài)。不過主進程讀取索引文件的速度很快，而且廣播可以在lgn（n為進程數(shù)量）時間內(nèi)完成，因此這是一種有效的方法。當大規(guī)模系統(tǒng)的互聯(lián)帶寬較高時，這種方法的優(yōu)勢更加明顯。

（2）讀入短讀序列

讀入索引后，需要讀入進行比對的短讀序列。BWA支持FASTQ格式的讀段文件。受限于內(nèi)存空間，當文件中讀段數(shù)量較大時，BWA會分批進行迭代處理，直到所有序列比對完畢。類似于讀入?yún)⒖蓟蚪M索引，BWA在讀入短讀序列時也采用單線程方式實現(xiàn)。讀入序列時，BWA會一直讀取，直到識別到FASTQ格式的行開頭信息（>、@或+字符）。

本文讀入短讀序列采用類似主從模式的方式，首先由一個指定的進程掃描短讀序列的FASTQ文件，獲取文件中的序列數(shù)量、文件行數(shù)等信息，然后根據(jù)進程數(shù)量對讀段進行劃分。如p號進程處理序列文件的前l(fā)條序列，其中N為文件行數(shù)，n為進程數(shù)量，p從0開始計數(shù)。FASTQ文件格式規(guī)定，每4行代表一條序列，因此，可以很容易地從序列標號計算出相應(yīng)的文件行號。主進程每次跳過l條序列時，會將讀取文件的位置發(fā)給相應(yīng)的進程，這樣每個進程都獲得了自己讀取文件的位置，也就完成了序列分發(fā)。這個過程實際上相當于為短讀序列文件建立了一個索引，只是這個索引并未被存儲下來。

（3）序列比對

經(jīng)過前兩個步驟，各進程已經(jīng)有了一份參考基因組索引副本以及分配給本進程的讀段，可以進行序列比對。由于序列之間不存在依賴性，各進程間的序列比對是相互獨立的，在每個進程內(nèi)部，BWA本身實現(xiàn)了多線程比對，進一步提高了比對的并發(fā)度和速度。

（4）輸出比對結(jié)果

并行程序的多進程間無法共享文件指針、變量和內(nèi)存，因此只能每個進程輸出單獨的文件。各短讀序列之間沒有相關(guān)性，當所有序列的比對都完成后，可以通過簡單的shell命令將所有結(jié)果文件直接合并成一個總的文件。

另一種方案是采用主從模式，由一個進程負責寫文件，其他進程將比對結(jié)果發(fā)給寫文件進程，但是這樣會對寫文件進程造成很大的壓力，因為在各進程并行比對且進程規(guī)模較大時，所有進程都使用一個I/O通道，非常擁擠，這會成為系統(tǒng)的瓶頸，影響整個軟件的性能。因此，本文采用每個進程單獨輸出文件的方式。

4.2 SNP檢測算法優(yōu)化

（1）程序熱點分析

首先，為了找到SOAPsnp程序的性能瓶頸，筆者分別測試了7個模塊的運行時間，發(fā)現(xiàn)likelihood和recycle兩個模塊耗時較長，分別占總時間的59%和30%，排在第三位的是output模塊，占比為8.3%。

然后筆者利用性能分析工具進行測試，測試結(jié)果顯示，likelihood模塊花費了大量時間在訪問主存，特別是對base_info的訪問。base_info用于存儲比對的堿基，likelihood模塊需要遍歷base_info來檢測SNP。recycle模塊負責初始化下一窗口，處理跨窗口的讀段，也包括復制base_info信息。

（2）四維矩陣壓縮降維優(yōu)化

在進行SNP檢測時，各位點是相互獨立的，在每個位點的計算中，base_info矩陣的每個元素都會被訪問一次，將產(chǎn)生131 072次訪存（4×64×256×2）。也就是說，整個人類基因組的32億個位點將產(chǎn)生4.19×1014次訪存。

通過分析base_info，筆者發(fā)現(xiàn)這是一個高度稀疏的矩陣。base_info的每個元素代表對四元組合(堿基, 測序質(zhì)量分數(shù),讀段中的偏移, 讀段方向)的計數(shù)，初始值為0。在堿基中，每發(fā)現(xiàn)一個滿足特定四元組合的堿基，就會將相應(yīng)元素增加1?？紤]到每個位點比對的堿基數(shù)量主要受測序深度的影響，而測序深度通常低于100層，base_info有131 072個元素，那么base_info矩陣的最大非零比例為100/131 072，也就是0.076%。從這些分析可以看出，base_info矩陣99.9%以上的元素為0，對于計算而言毫無意義，只有不到0.1%的訪存是有效的，但是算法仍會遍歷整個矩陣，這造成了大量的無用計算。

基于base_info的高度稀疏性，為提高訪存效率，筆者直接使用一維數(shù)組base_array替代四維矩陣base_info。由于堿基類型、讀段方向、測序質(zhì)量分數(shù)和讀段中的偏移這些信息的最大值都是確定的，直接通過位操作將這些信息封裝在一個32位的字中。在counting模塊中，直接將比對的堿基信息依次存儲到base_array中，從而避免對base_info矩陣的多次訪問，提高程序的空間局部性和cache命中率。

采用一維數(shù)組替換四維矩陣后，likelihood計算模塊由原來的5層循環(huán)變?yōu)?層循環(huán)（見表1），大大提高了計算效率。另外，放棄使用base_info矩陣后，recycle模塊復制的數(shù)據(jù)也減少了很多，不到原來的0.1%。

（3）基于快表的去冗計算

likelihood模塊最耗時的部分是更新type_likely的數(shù)據(jù)結(jié)構(gòu)。在算法中，每次更新type_likely需要分別計算10種基因型（見表2）的概率值。每個堿基位點需要計算一次type_likely，對于人類基因組而言，需要訪問type_likely數(shù)十億次。每次更新type_likely時，需要訪問p_matrix矩陣兩次，取出相應(yīng)數(shù)值，再進行一個對數(shù)操作來計算概率，這些操作的時鐘周期都較長，且type_likely更新的次數(shù)很高。

表1 優(yōu)化前后的堿基信息遍歷算法

為提高type_likely的計算速度，筆者采用了一種用空間換時間的策略。首先，筆者觀察到type_likely每次更新只涉及10種基因型，數(shù)量固定，而訪問的p_matrix矩陣值和規(guī)模也是固定的。因此，對于p_matrix矩陣的每個元素，可以提前計算出10種基因型對應(yīng)的10個值，并將其保存到內(nèi)存中，當type_likely需要訪問p_matrix的某個元素時，直接從預(yù)先算好的矩陣（預(yù)計算矩陣）中讀取。p_matrx矩陣的大小為8 MB，對應(yīng)10種基因型，則預(yù)計算矩陣的大小為80 MB。雖然這樣犧牲了一定的內(nèi)存空間，但是可以將原來的兩次訪存及一次對數(shù)操作變?yōu)橐淮卧L存操作，降低了算法的復雜度。

（4）I/O壓縮策略

對SOAPsnp進行算法改進后，計算速度顯著提高。但是隨之而來的一個問題是，優(yōu)化前占程序運行總時間8.3%的輸出模塊占優(yōu)化后程序總時間的50%以上，成為新的瓶頸。為提高程序的I/O性能，筆者進一步分析了SOAPsnp程序。輸入文件主要為比對完成的讀段，是其他軟件的輸出，因此，暫時不考慮優(yōu)化。SOAPsnp的輸出為一個17列的文本文件，如圖4所示。

通過觀察SOAPsnp的輸出文件發(fā)現(xiàn)，每一列都存在很多重復的值，可以針對每列采取不同的壓縮方式。其中第1列為參考序列名稱，圖4中的chr22是染色體名稱，這一列的值完全相同，因此只需要存儲一次。第2列為堿基在參考序列上的偏移，每行遞增，因此也只需要存儲第一個值。對于存儲堿基值的第3、4、6列，由于SNP的占比非常低，因此這3列的絕大部分值是相同的，可以合并存儲，只需要額外存儲那些是SNP的位點。與基因型相關(guān)的列都是高度稀疏的，只需要存儲非零值即可。第11、12、13、15、17列為默認值，只有存在SNP時該值才發(fā)生變化。因此，這幾列也只需要存儲那些不同的值和一個默認值即可。其他列則根據(jù)特性，可分別采用游程編碼或字典等方式進行壓縮。

4.3 SNP檢測算法的并行優(yōu)化

各個位點的SNP檢測是相互獨立的，具有天然的并行性，筆者采用OpenMP和MPI技術(shù)實現(xiàn)了多級并行SNP檢測。在窗口內(nèi)部，各位點之間采用基于OpenMP的多線程并行，每個線程計算一個位點，這樣可以利用多線程共享數(shù)據(jù)的特性，降低內(nèi)存壓力。此外，為提高SNP檢測速度，筆者還利用MPI技術(shù)實現(xiàn)了多窗口間的并行，不同的進程運行在不同的計算節(jié)點上，分別處理不同的窗口。對于跨窗口的讀段，在開始每個窗口的計算前，相鄰窗口會進行兩輪通信，第一次通信是確認每個窗口的位置，確定沒有重疊和遺漏；第二次通信是交換跨窗口的讀段的信息，相當于圖3中的recycle模塊，對窗口進行初始化，這樣就可以保證讀段信息的完整性和SNP檢測的正確性。

表2 基因型的表示方法

4.4 變異檢測分析流水線I/O整合優(yōu)化

在基因組變異檢測分析流水線中，序列比對和SNP檢測是相鄰的兩個步驟，序列比對軟件首先將比對結(jié)果輸出到文件，由于SNP檢測需要基于排序好的讀段進行，因此，還需要對輸出的比對結(jié)果進行排序，并行輸出排序好的比對文件，SNP檢測以排好序的比對文件為輸入進行分析。以上流程需要3次大規(guī)模讀操作和3次大規(guī)模寫操作，而基因組規(guī)模一般較大，且高通量測序通常需要測多層，以保證數(shù)據(jù)的魯棒性，因此，比對、排序和SNP檢測對I/O造成了巨大壓力，特別是當程序不斷優(yōu)化，計算時間不斷縮短，而數(shù)據(jù)規(guī)模不斷增長時，I/O成為程序的主要瓶頸。

針對這一問題，本文提出了一種I/O整合策略。當序列比對完成后，將比對結(jié)果暫時保存在內(nèi)存中，對結(jié)果進行排序，然后直接基于排序后的結(jié)果進行SNP檢測，這樣可以避免兩輪I/O操作。另外，程序也會正常輸出比對結(jié)果，用于后續(xù)其他分析。這樣一來，只需要讀入一次參考基因組和序列文件，即可完成比對和SNP檢測兩個步驟，大大降低了I/O壓力，提高了數(shù)據(jù)處理的速度。

5 測試分析

5.1 測試環(huán)境與數(shù)據(jù)

本文的測試主要在“天河二號”超級計算機上進行，節(jié)點配置見表3。

實驗采用的物種為人類，參考基因組版本為hg38，讀段數(shù)據(jù)包括真實數(shù)據(jù)和模擬數(shù)據(jù)，真實數(shù)據(jù)由Illumina測序儀產(chǎn)生，模擬數(shù)據(jù)由wgsim生成，長度均為100堿基對（base pair，BP）。

5.2 序列比對測試分析

（1）強可擴展性

本節(jié)通過對500萬條序列進行比對來測試程序的強可擴展性，進程數(shù)分別為1個、2個、4個、8個、16個、32個、64個。測試數(shù)據(jù)見表4。

表3 “天河二號”超級計算機節(jié)點配置

如圖5所示，若保持問題規(guī)模不變，加入更多的處理器，程序處理時間接近于線性減少，這證明程序具有良好的可擴展性。雖然每個進程都要讀入一份參考索引，但是采取進程廣播的形式來分發(fā)參考索引并不會讓程序總體處理時間有明顯增加，而總體的程序處理時間大大減少，運行時間接近于理想情況的1/n，表現(xiàn)出良好的強可擴展性。

（2）弱可擴展性

在弱可擴展性測試中，令每個進程處理50萬條序列，進程數(shù)分別為1個、2個、4個、8個、16個、32個、64個，數(shù)據(jù)規(guī)模隨進程數(shù)的增加同比擴大。測試數(shù)據(jù)見表5。

如圖6所示，若在進程數(shù)增加的同時同比擴大數(shù)據(jù)規(guī)模，每個進程平均處理50萬條序列，程序運行時間保持穩(wěn)定。從單進程到多進程，處理數(shù)據(jù)的時間會有明顯增加，這主要是因為多個進程在輸入、輸出數(shù)據(jù)時需要更多時間，并且進程之間需要通信，這導致了程序執(zhí)行時間的增加。在程序并行執(zhí)行時，平均每個進程運行50萬條序列，運行時間波動不大，基本保持穩(wěn)定。這說明，在增加數(shù)據(jù)規(guī)模的同時增加進程數(shù)可以使處理機在相同的時間內(nèi)處理更多的數(shù)據(jù)，程序的弱可擴展性良好。

5.3 SNP檢測測試分析

（1）算法優(yōu)化效果測試

本文首先測試了各種算法優(yōu)化策略的加速效果，以SOAPsnp單線程程序為基準，逐一添加優(yōu)化策略，并測試運行時間，如圖7所示。在圖7中，橫軸的SOAPsnp表示原始算法，o1表示base_info壓縮優(yōu)化，o1+o2表示在o1的基礎(chǔ)上添加快表優(yōu)化，o1+o2+o3表示在o1+o2的基礎(chǔ)上進行了輸出壓縮優(yōu)化。

從圖7可以看到，效果最顯著的是base_info壓縮優(yōu)化，因為壓縮后不僅大大減少了計算量，還降低了空間開銷，一次優(yōu)化降低了likelihood和recycle兩個模塊的運行時間。快表優(yōu)化是對likelihood計算核心的優(yōu)化，以時間換空間，起到了很好的效果。輸出壓縮優(yōu)化則在計算模塊加速充分的情況下，進一步降低了程序運行時間，提高了處理速度。

（2）并行程序可擴展性

本文利用OpenMP+MPI的并行技術(shù)實現(xiàn)了多級并行。圖8展示的是程序的強可擴展性。以128個節(jié)點的性能為基準，程序運行時間近似呈線性下降，表明程序具有良好的可擴展性。圖9展示的是程序弱可擴展性的并行效率，測試數(shù)據(jù)集的大小隨節(jié)點規(guī)模同比擴大，最大數(shù)據(jù)規(guī)模達到了542 GB，從圖9可以看出，程序并行效率隨著計算規(guī)模的擴大呈下降趨勢，雖然筆者對程序的通信和輸出進行了優(yōu)化，但是當數(shù)據(jù)量和進程規(guī)模較大時，掃描讀段文件、進程通信和結(jié)果輸出所占的比例會不斷提高，計算占比會相對降低，導致程序的并行效率降低，但是總體并行效率仍保持在60%以上。

6 結(jié)束語

本文選取基因組大數(shù)據(jù)變異檢測中最耗時的兩個環(huán)節(jié)——序列比對和SNP檢測，首先從算法角度進行了優(yōu)化和改進，然后利用OpenMP、MPI等并行技術(shù)實現(xiàn)了多級并行處理。測試表明，并行優(yōu)化后的算法具有良好的并行性能和可擴展性，為快速處理高通量測序技術(shù)帶來的海量基因組大數(shù)據(jù)提供了有力的支持。

表4 同一讀段長度下不同進程數(shù)的程序運行時間

表5 不同進程數(shù)和讀段長度下程序運行時間