劉佳俊，胡大裟，蔣玉明

（四川大學(xué)計(jì)算機(jī)學(xué)院，成都610065）

0 引言

基因測(cè)序與變異鑒定技術(shù)在基因研究領(lǐng)域中發(fā)揮著非常重要的作用。隨著下一代測(cè)序技術(shù)[1]誕生，測(cè)序成本下降，運(yùn)用重測(cè)序研究群體遺傳變異成為一種可行方法。然而隨著測(cè)序深度和樣本數(shù)目劇增，面對(duì)海量的基因序列，鑒定并尋找單核苷酸多態(tài)性（Single Nucleotide Polymorphism，SNP）[2-3]位點(diǎn)需要大量的計(jì)算資源和時(shí)間，基因數(shù)據(jù)的分析速度成為研究的瓶頸之一。在保證鑒定的高精度前提下，如何提高基因鑒定的速度，降低時(shí)間成本，已成為基因鑒定中的關(guān)鍵技術(shù)之一。

近年來(lái)，已經(jīng)有許多學(xué)者將基因測(cè)序技術(shù)與計(jì)算機(jī)軟件研究相結(jié)合，如：BWA[4]軟件在短序列的分析應(yīng)用中能達(dá)到較高的準(zhǔn)確率。Heng Li[5]使用無(wú)衍生物最小化（Minimization Without Derivatives，MWD）等方法[6]估算等位基因頻率，并提出了一種突變識(shí)別（variant calling）的統(tǒng)計(jì)學(xué)框架。Aaron McKenna 等[7]提出了一個(gè)使用分布式和共享內(nèi)存并行化的結(jié)構(gòu)化編程框架—基因組分析工具包（The Genome Analysis Toolkit，GATK）。Dries Decap 等[8]則根據(jù)MapReduce[9]編程模型在多節(jié)點(diǎn)環(huán)境下實(shí)現(xiàn)并行分布式內(nèi)存計(jì)算。M.C.Schatz[10]使用隨機(jī)微衛(wèi)星擴(kuò)增多態(tài)DNA（RMAP）建模，提出一種通過(guò)Hadoop 實(shí)現(xiàn)的高度敏感并行種子擴(kuò)增讀取映射算法。

面對(duì)日益增長(zhǎng)的基因數(shù)據(jù)，多線程的基因分析工具出現(xiàn)了計(jì)算耗時(shí)過(guò)長(zhǎng)的問(wèn)題。而在實(shí)際應(yīng)用中，往往又需要進(jìn)行多個(gè)樣本的SNP 鑒定，以SNP 鑒定中主要使用的工具SAMtools 為例，對(duì)一個(gè)約為4GB 的測(cè)序樣本，鑒定單個(gè)樣本的計(jì)算耗時(shí)約為三至四天；而當(dāng)樣本數(shù)量增加至10 個(gè)時(shí)，計(jì)算時(shí)間隨之增至約一個(gè)月；當(dāng)樣本數(shù)量增加至200 個(gè)，計(jì)算時(shí)間至少需要兩個(gè)月。

本文設(shè)計(jì)了一個(gè)基于SPMD（Single Program，Multiple Data）[11]模式的并行化模型，在高性能計(jì)算平臺(tái)上結(jié)合傳統(tǒng)鑒定工具的優(yōu)點(diǎn)，根據(jù)MapReduce 框架將數(shù)據(jù)劃分和映射，把計(jì)算任務(wù)拆分并利用多節(jié)點(diǎn)進(jìn)行計(jì)算，并實(shí)現(xiàn)多樣本的并行計(jì)算，充分利用計(jì)算平臺(tái)多節(jié)點(diǎn)多CPU 資源，解決傳統(tǒng)工具在高性能計(jì)算平臺(tái)上性能較低的問(wèn)題，提高遺傳變異鑒定的計(jì)算速度以及加速變異檢測(cè)過(guò)程。

1 多節(jié)點(diǎn)并行化模型

MapReduce 是用于處理大量數(shù)據(jù)的一種支持并行化的分布式模型，MapReduce 主要有Map（映射）和Reduce（歸約）兩部分組成。Map 函數(shù)根據(jù)輸入基因數(shù)據(jù)，劃分出若干個(gè)參考序列——樣本序列組合，Reduce函數(shù)收集包含相同序列區(qū)間并對(duì)相應(yīng)位點(diǎn)進(jìn)行單核苷酸多態(tài)性鑒定。

SNP 鑒定主要是通過(guò)將參考序列和對(duì)應(yīng)的若干個(gè)樣本序列進(jìn)行組合，對(duì)每個(gè)基因位點(diǎn)的匹配程度進(jìn)行統(tǒng)計(jì)。假設(shè)F 為某一染色體參考序列，Si為樣本代號(hào)i在該染色體上的對(duì)應(yīng)序列，則{Si,Sj,…,Sn}是n 個(gè)樣本在該染色體上的樣本序列集合。將F 和{Si,Sj,…,Sn}的參考序列——樣本序列進(jìn)行組合；假設(shè)參考序列上的某一基因位點(diǎn)為GF，GS是樣本序列上與GF對(duì)應(yīng)位置的位點(diǎn)，將包含GF位點(diǎn)信息的read 序列與GS進(jìn)行匹配等計(jì)算統(tǒng)計(jì)，評(píng)估GS處發(fā)生SNP 的可能性。

1.1 多樣本多節(jié)點(diǎn)的并行計(jì)算模型

Map 函數(shù)將參考序列和樣本序列分別進(jìn)行劃分，得到若干個(gè)Fk和對(duì)應(yīng)的{Si-k,Sj-k,…,Sn-k}的參考序列——樣本序列組合，然后將具有相同序列區(qū)間的多個(gè)樣本數(shù)據(jù)進(jìn)行合并，再由Reduce 函數(shù)收集和歸約并計(jì)算這些參考序列——樣本序列組合的基于Phred 的質(zhì)量得分Q（Phred-like Quality Scores），經(jīng)序列對(duì)齊等相關(guān)處理，最后過(guò)濾篩選，輸出SNP 候選位點(diǎn)信息[12]。各個(gè)節(jié)點(diǎn)上的計(jì)算流程如圖1 所示。

圖1 各計(jì)算節(jié)點(diǎn)上的計(jì)算流程

（1）參考數(shù)據(jù)：參考的基因組數(shù)據(jù)，包含模板生物序列的DNA 或者蛋白質(zhì)序列信息。

（2）樣本數(shù)據(jù)：實(shí)際樣本測(cè)序得到的基因數(shù)據(jù)，包含了多個(gè)read 比對(duì)信息，用于基因鑒定與分析。

（3）參考數(shù)據(jù)劃分：以計(jì)算節(jié)點(diǎn)數(shù)量和參考數(shù)據(jù)大小為參考，將參考序列分割為合理大小用于并行計(jì)算。

（4）樣本數(shù)據(jù)劃分：以參考分割程序的劃分地址為參考，將樣本序列分割為合理大小用于并行計(jì)算。

（5）變異鑒定：將劃分后的參考數(shù)據(jù)和樣本數(shù)據(jù)映射匹配，對(duì)序列中的每一個(gè)位點(diǎn)進(jìn)行SNP 概率估計(jì)。

（6）多節(jié)點(diǎn)鑒定結(jié)果合并：將多個(gè)計(jì)算節(jié)點(diǎn)生成的結(jié)果文件按照劃分地址的順序進(jìn)行排列與組合。

1.2 任務(wù)分配方案

動(dòng)態(tài)負(fù)載優(yōu)先（Dynamic Load Priority，DLP）方案首次任務(wù)分配時(shí)隨機(jī)選擇若干節(jié)點(diǎn)，通過(guò)控制主機(jī)進(jìn)行主動(dòng)探測(cè)，獲得當(dāng)前子節(jié)點(diǎn)的負(fù)載情況和計(jì)算進(jìn)度，再對(duì)子節(jié)點(diǎn)的負(fù)載和資源量進(jìn)行優(yōu)先級(jí)判定。在計(jì)算平臺(tái)上有若干計(jì)算節(jié)點(diǎn)N={N1,N2,…,Nn}，其中節(jié)點(diǎn)Ni的空余資源為Rf，假設(shè)節(jié)點(diǎn)Ni上已有k 個(gè)任務(wù)正在計(jì)算，其中任務(wù)Tj占用的資源為Rj，對(duì)應(yīng)的計(jì)算進(jìn)度百分比為Sj；當(dāng)前進(jìn)行任務(wù)分配的計(jì)算資源請(qǐng)求為R，則節(jié)點(diǎn)Ni的分配優(yōu)先級(jí)Pi為：

當(dāng)控制主機(jī)進(jìn)行再分配時(shí)會(huì)結(jié)合當(dāng)前各節(jié)點(diǎn)資源的優(yōu)先級(jí)來(lái)指派任務(wù)，選擇優(yōu)先級(jí)別高的節(jié)點(diǎn)進(jìn)行再分配。

1.3 Map函數(shù)和Reduce函數(shù)

根據(jù)SNP 鑒定中數(shù)據(jù)文件的種類[13]，進(jìn)行數(shù)據(jù)劃分的Map 處理可以分為參考序列和樣本序列兩種：參考序列的Map 處理首先需要按照染色體種類進(jìn)行拆分，然后通過(guò)設(shè)置并行數(shù)量，對(duì)拆分后的單染色體參考序列，進(jìn)行平均分割操作以實(shí)現(xiàn)各節(jié)點(diǎn)負(fù)載均衡；樣本序列的Map 處理會(huì)根據(jù)參考序列Map 的分割結(jié)果，提取并劃分對(duì)應(yīng)區(qū)間的樣本序列。

Map 函數(shù)的輸入為（key,value）鍵值對(duì)，其中key 為參考序列每個(gè)分割區(qū)間的起始地址，value 為對(duì)應(yīng)區(qū)間的參考序列。Map 函數(shù)將參考基因組根據(jù)區(qū)間長(zhǎng)度或者測(cè)序深度進(jìn)行平均劃分，保證每段數(shù)據(jù)量相近。Map函數(shù)完成數(shù)據(jù)劃分后，生成對(duì)應(yīng)的（key,value）文件存儲(chǔ)在磁盤上，進(jìn)入下一步計(jì)算操作。

參考序列的Map 函數(shù)描述如下：

函數(shù)1 參考序列的Map 函數(shù)。

輸入：參考序列數(shù)據(jù)。

輸出：參考序列鍵值對(duì)。

Step1 控制節(jié)點(diǎn)將輸入的參考序列文件按照其中數(shù)據(jù)的染色體名稱進(jìn)行分離，得到多個(gè)單染色體參考序列文件。

Step2 控制節(jié)點(diǎn)獲取指定的劃分?jǐn)?shù)量，確定進(jìn)行劃分的參考序列的最小閾值。對(duì)上述每個(gè)參考序列進(jìn)行Step3 到Step4 的劃分操作。

Step3 獲取數(shù)據(jù)文件的大小，選擇根據(jù)劃分?jǐn)?shù)量進(jìn)行平均劃分或者根據(jù)測(cè)序深度進(jìn)行劃分。

Step4 建立起始地址——對(duì)應(yīng)劃分區(qū)間數(shù)據(jù)的鍵值對(duì)，輸出至硬盤上并結(jié)束。

參考序列的Map 函數(shù)運(yùn)行在控制節(jié)點(diǎn)上，時(shí)間開銷取決于輸入的參考序列大小，參考序列主要由字符組成，遍歷的速度較快?？刂乒?jié)點(diǎn)統(tǒng)計(jì)輸出的起始地址——?jiǎng)澐謹(jǐn)?shù)據(jù)鍵值對(duì)數(shù)量，根據(jù)DLP 方案進(jìn)行任務(wù)分配。在每個(gè)計(jì)算節(jié)點(diǎn)均需上執(zhí)行樣本序列——Map函數(shù)以及相關(guān)處理。

在參考序列的Map 函數(shù)完成參考序列的分割后，每個(gè)計(jì)算任務(wù)并行進(jìn)行每個(gè)任務(wù)的排序和提取等操作。在提取對(duì)應(yīng)區(qū)間內(nèi)的序列信息時(shí)，可能會(huì)導(dǎo)致在區(qū)間兩端的信息出現(xiàn)缺失，從而影響鑒定結(jié)果。為了減小和避免這種情況，在提取序列信息時(shí)額外增加提取區(qū)間兩端的范圍，可確保序列信息的完整。

樣本序列的Map 函數(shù)描述如下：

函數(shù)2 樣本序列的Map 函數(shù)。

輸入：樣本序列數(shù)據(jù)以及起始地址——?jiǎng)澐謹(jǐn)?shù)據(jù)鍵值對(duì)。

輸出：樣本序列鍵值對(duì)。

Step1 計(jì)算節(jié)點(diǎn)把樣本序列數(shù)據(jù)根據(jù)參考序列中的染色體進(jìn)行分離。

Step2 對(duì)每個(gè)染色體的樣本序列進(jìn)行區(qū)間統(tǒng)計(jì)，判斷當(dāng)前地址區(qū)間內(nèi)是否有相關(guān)的讀長(zhǎng)read 信息[14]，如果區(qū)間內(nèi)包含相關(guān)信息則根據(jù)參考序列中的地址區(qū)間進(jìn)行劃分，劃分時(shí)額外增加的區(qū)間長(zhǎng)度為l，通過(guò)根據(jù)各個(gè)read 的長(zhǎng)度和出現(xiàn)的頻次的加權(quán)計(jì)算，如公式（2）；如果沒有則輸出空白文件。

式（2）中，Li為編號(hào)為i 的read 序列的長(zhǎng)度，Ni為編號(hào)為i 的read 序列出現(xiàn)的次數(shù)。

Step3 輸出樣本序列鍵值對(duì)至硬盤上，完成參考序列-樣本序列組的劃分映射并結(jié)束。

在經(jīng)過(guò)Map 函數(shù)劃分后，計(jì)算節(jié)點(diǎn)對(duì)輸出的參考序列——樣本序列組合進(jìn)行區(qū)間地址的校驗(yàn)，校驗(yàn)通過(guò)后使用Reduce 函數(shù)歸約并進(jìn)行SNP 鑒定計(jì)算。若校驗(yàn)未通過(guò)重新進(jìn)行Step2。

Reduce 函數(shù)描述如下：

函數(shù)3 Reduce 函數(shù)。

輸入：參考序列——樣本序列組合。

輸出：SNP 鑒定結(jié)果文件。

Step1 計(jì)算節(jié)點(diǎn)多線程讀入多個(gè)起始區(qū)間的參考序列，并收集所有樣本中具有相同區(qū)間的樣本序列。其中每個(gè)線程讀入一個(gè)計(jì)算單元數(shù)據(jù)，即一個(gè)參考序列和與其對(duì)應(yīng)的若干個(gè)樣本序列。

Step2 對(duì)每一組鑒定數(shù)據(jù)單元，計(jì)算序列中的堿基位點(diǎn)的等位基因的先驗(yàn)概率，根據(jù)概率選出最有可能變異的位點(diǎn)。

Step3 計(jì)算節(jié)點(diǎn)將SNP 結(jié)果輸出至硬盤并結(jié)束。

控制節(jié)點(diǎn)會(huì)持續(xù)詢問(wèn)所有計(jì)算節(jié)點(diǎn)直到完成SNP鑒定，然后控制節(jié)點(diǎn)把所有結(jié)果根據(jù)劃分的地址順序進(jìn)行合并。

2 實(shí)驗(yàn)測(cè)試與結(jié)果分析

實(shí)驗(yàn)測(cè)試有2 臺(tái)曙光A840R-G，4×AMD 6172 12核CPU，512G 內(nèi)存，共4 個(gè)計(jì)算節(jié)點(diǎn)，操作系統(tǒng)為L(zhǎng)inux version 2.6.32-431.el6.x86_64；Red Hat Enterprise Linux Server release 6.5，作業(yè)管理系統(tǒng)為L(zhǎng)SF。

實(shí)驗(yàn)數(shù)據(jù)以某地某品種（系）煙草材料為作為測(cè)試樣本，參考序列為隨機(jī)選擇的其第17 號(hào)染色體數(shù)據(jù)，樣本序列則為該品種（系）煙草的128 個(gè)樣本測(cè)序數(shù)據(jù)，每個(gè)樣本的文件大小約為4GB。實(shí)驗(yàn)主要分析樣本數(shù)和并行數(shù)對(duì)計(jì)算時(shí)間的影響，以SAMtools 工具在單節(jié)點(diǎn)的計(jì)算時(shí)間為對(duì)比。第一部分實(shí)驗(yàn)，固定計(jì)算的樣本規(guī)模，同時(shí)改變參與計(jì)算的節(jié)點(diǎn)數(shù)量；第二部分實(shí)驗(yàn)，固定每個(gè)節(jié)點(diǎn)所計(jì)算的樣本數(shù)量，通過(guò)增加節(jié)點(diǎn)數(shù)量，分析計(jì)算效率的變化。

采用本文提出的并行優(yōu)化模型在多節(jié)點(diǎn)上設(shè)置并行數(shù)依次為2 組、3 組、5 組、10 組、15 組、20 組、30 組、50 組、100 組，其中，每組并行數(shù)再分別計(jì)算樣本數(shù)1個(gè)、2 個(gè)、4 個(gè)、8 個(gè)、16 個(gè)、32 個(gè)、64 個(gè)、128 個(gè)，每個(gè)樣本的大小約4GB。通過(guò)上述計(jì)算，以運(yùn)算時(shí)間和加速比[15]兩個(gè)指標(biāo)來(lái)評(píng)價(jià)本文模型，統(tǒng)計(jì)運(yùn)算時(shí)間如圖2 所示，加速比如圖3 所示。

圖2 樣本數(shù)和并行數(shù)對(duì)計(jì)算效率的影響

由圖2 可看出：在并行數(shù)保持不變的情況下，樣本數(shù)與計(jì)算時(shí)間成正比，即隨著樣本數(shù)的增多，排序等操作所需遍歷的文件總大小也在增大，平均每個(gè)基因位點(diǎn)的堆疊深度在增加，呈現(xiàn)出計(jì)算時(shí)間成倍增加；而在樣本數(shù)保持不變的情況下，在并行數(shù)小于20 時(shí)，增加并行數(shù)可顯著降低計(jì)算時(shí)間，即通過(guò)實(shí)現(xiàn)并行化，增加并行計(jì)算數(shù)量，計(jì)算時(shí)間會(huì)呈現(xiàn)明顯下降的趨勢(shì)。

當(dāng)樣本數(shù)小于16 個(gè)時(shí)，多組數(shù)據(jù)相互對(duì)比，計(jì)算時(shí)間差距并不明顯，使用并行化加速效果并不明顯，在實(shí)際應(yīng)用中可以考慮使用傳統(tǒng)計(jì)算方法計(jì)算，避免不必要的計(jì)算誤差；但當(dāng)樣本數(shù)為16 個(gè)至128 個(gè)時(shí)，采用并行加速模型可明顯增加提高計(jì)算效率。

由圖3 可看出隨著并行數(shù)的增加，加速比在并行數(shù)20 以內(nèi)呈現(xiàn)明顯上升趨勢(shì)，說(shuō)明增加并行數(shù)可有效地減小每個(gè)任務(wù)的計(jì)算時(shí)間，從而減少總體計(jì)算時(shí)間；加速比在并行數(shù)量30 到50 這個(gè)區(qū)間出現(xiàn)緩慢上升甚至持平的趨勢(shì)，這是因?yàn)楫?dāng)每個(gè)節(jié)點(diǎn)的計(jì)算時(shí)間縮小到一定程度后，SNP 鑒定的時(shí)間主要由相對(duì)固定耗時(shí)的I/O 操作效率所確定，并且已經(jīng)相當(dāng)接近，所以導(dǎo)致加速比不再明顯增加；當(dāng)并行數(shù)在50 以上后，加速比將保持恒定。由圖2 和圖3 綜合也說(shuō)明并行數(shù)量并不是并非簡(jiǎn)單的越大越好，而是與樣本數(shù)有直接關(guān)系。

圖3 樣本數(shù)和并行數(shù)對(duì)多節(jié)點(diǎn)計(jì)算加速比的影響

3 結(jié)語(yǔ)

針對(duì)基因變異鑒定耗時(shí)較大的問(wèn)題，通過(guò)MapReduce 框架模型將計(jì)算任務(wù)劃分，在多節(jié)點(diǎn)上實(shí)現(xiàn)并行加速計(jì)算，提出了一種針對(duì)多樣本的多節(jié)點(diǎn)并行化優(yōu)化算法，在海量基因數(shù)據(jù)的計(jì)算與分析中具有廣泛適用性。

通過(guò)SNP 鑒定的計(jì)算實(shí)驗(yàn)發(fā)現(xiàn)，樣本數(shù)平均每增加一倍，計(jì)算時(shí)間也平均增加約一倍；通過(guò)增加并行數(shù)加速計(jì)算，可使計(jì)算時(shí)間明顯下降：并行數(shù)平均每增加一倍，計(jì)算時(shí)間平均減少約32%，加速比平均上升約46%。驗(yàn)證了該優(yōu)化算法是有效的。同時(shí)，并行計(jì)算的加速效果受樣本數(shù)和并行數(shù)影響，在小數(shù)據(jù)量或多余并行下無(wú)法顯示其優(yōu)越性，實(shí)際計(jì)算中應(yīng)根據(jù)樣本數(shù)和硬件條件，合理規(guī)劃才能獲得更有效的加速效果。