基于MapReduce 的多樣本基因鑒定并行優化模型

劉佳俊，胡大裟，蔣玉明

（四川大學計算機學院，成都610065）

0 引言

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

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

面對日益增長的基因數據，多線程的基因分析工具出現了計算耗時過長的問題。而在實際應用中，往往又需要進行多個樣本的SNP 鑒定，以SNP 鑒定中主要使用的工具SAMtools 為例，對一個約為4GB 的測序樣本，鑒定單個樣本的計算耗時約為三至四天；而當樣本數量增加至10 個時，計算時間隨之增至約一個月；當樣本數量增加至200 個，計算時間至少需要兩個月。

本文設計了一個基于SPMD（Single Program，Multiple Data）[11]模式的并行化模型，在高性能計算平臺上結合傳統鑒定工具的優點，根據MapReduce 框架將數據劃分和映射，把計算任務拆分并利用多節點進行計算，并實現多樣本的并行計算，充分利用計算平臺多節點多CPU 資源，解決傳統工具在高性能計算平臺上性能較低的問題，提高遺傳變異鑒定的計算速度以及加速變異檢測過程。

1 多節點并行化模型

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

SNP 鑒定主要是通過將參考序列和對應的若干個樣本序列進行組合，對每個基因位點的匹配程度進行統計。假設F 為某一染色體參考序列，Si為樣本代號i在該染色體上的對應序列，則{Si,Sj,…,Sn}是n 個樣本在該染色體上的樣本序列集合。將F 和{Si,Sj,…,Sn}的參考序列——樣本序列進行組合；假設參考序列上的某一基因位點為GF，GS是樣本序列上與GF對應位置的位點，將包含GF位點信息的read 序列與GS進行匹配等計算統計，評估GS處發生SNP 的可能性。

1.1 多樣本多節點的并行計算模型

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

圖1 各計算節點上的計算流程

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

（2）樣本數據：實際樣本測序得到的基因數據，包含了多個read 比對信息，用于基因鑒定與分析。

（3）參考數據劃分：以計算節點數量和參考數據大小為參考，將參考序列分割為合理大小用于并行計算。

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

（5）變異鑒定：將劃分后的參考數據和樣本數據映射匹配，對序列中的每一個位點進行SNP 概率估計。

（6）多節點鑒定結果合并：將多個計算節點生成的結果文件按照劃分地址的順序進行排列與組合。

1.2 任務分配方案

動態負載優先（Dynamic Load Priority，DLP）方案首次任務分配時隨機選擇若干節點，通過控制主機進行主動探測，獲得當前子節點的負載情況和計算進度，再對子節點的負載和資源量進行優先級判定。在計算平臺上有若干計算節點N={N1,N2,…,Nn}，其中節點Ni的空余資源為Rf，假設節點Ni上已有k 個任務正在計算，其中任務Tj占用的資源為Rj，對應的計算進度百分比為Sj；當前進行任務分配的計算資源請求為R，則節點Ni的分配優先級Pi為：

當控制主機進行再分配時會結合當前各節點資源的優先級來指派任務，選擇優先級別高的節點進行再分配。

1.3 Map函數和Reduce函數

根據SNP 鑒定中數據文件的種類[13]，進行數據劃分的Map 處理可以分為參考序列和樣本序列兩種：參考序列的Map 處理首先需要按照染色體種類進行拆分，然后通過設置并行數量，對拆分后的單染色體參考序列，進行平均分割操作以實現各節點負載均衡；樣本序列的Map 處理會根據參考序列Map 的分割結果，提取并劃分對應區間的樣本序列。

Map 函數的輸入為（key,value）鍵值對，其中key 為參考序列每個分割區間的起始地址，value 為對應區間的參考序列。Map 函數將參考基因組根據區間長度或者測序深度進行平均劃分，保證每段數據量相近。Map函數完成數據劃分后，生成對應的（key,value）文件存儲在磁盤上，進入下一步計算操作。

參考序列的Map 函數描述如下：

函數1 參考序列的Map 函數。

輸入：參考序列數據。

輸出：參考序列鍵值對。

Step1 控制節點將輸入的參考序列文件按照其中數據的染色體名稱進行分離，得到多個單染色體參考序列文件。

Step2 控制節點獲取指定的劃分數量，確定進行劃分的參考序列的最小閾值。對上述每個參考序列進行Step3 到Step4 的劃分操作。

Step3 獲取數據文件的大小，選擇根據劃分數量進行平均劃分或者根據測序深度進行劃分。

Step4 建立起始地址——對應劃分區間數據的鍵值對，輸出至硬盤上并結束。

參考序列的Map 函數運行在控制節點上，時間開銷取決于輸入的參考序列大小，參考序列主要由字符組成，遍歷的速度較快。控制節點統計輸出的起始地址——劃分數據鍵值對數量，根據DLP 方案進行任務分配。在每個計算節點均需上執行樣本序列——Map函數以及相關處理。

在參考序列的Map 函數完成參考序列的分割后，每個計算任務并行進行每個任務的排序和提取等操作。在提取對應區間內的序列信息時，可能會導致在區間兩端的信息出現缺失，從而影響鑒定結果。為了減小和避免這種情況，在提取序列信息時額外增加提取區間兩端的范圍，可確保序列信息的完整。

樣本序列的Map 函數描述如下：

函數2 樣本序列的Map 函數。

輸入：樣本序列數據以及起始地址——劃分數據鍵值對。

輸出：樣本序列鍵值對。

Step1 計算節點把樣本序列數據根據參考序列中的染色體進行分離。

Step2 對每個染色體的樣本序列進行區間統計，判斷當前地址區間內是否有相關的讀長read 信息[14]，如果區間內包含相關信息則根據參考序列中的地址區間進行劃分，劃分時額外增加的區間長度為l，通過根據各個read 的長度和出現的頻次的加權計算，如公式（2）；如果沒有則輸出空白文件。

式（2）中，Li為編號為i 的read 序列的長度，Ni為編號為i 的read 序列出現的次數。

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

在經過Map 函數劃分后，計算節點對輸出的參考序列——樣本序列組合進行區間地址的校驗，校驗通過后使用Reduce 函數歸約并進行SNP 鑒定計算。若校驗未通過重新進行Step2。

Reduce 函數描述如下：

函數3 Reduce 函數。

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

輸出：SNP 鑒定結果文件。

Step1 計算節點多線程讀入多個起始區間的參考序列，并收集所有樣本中具有相同區間的樣本序列。其中每個線程讀入一個計算單元數據，即一個參考序列和與其對應的若干個樣本序列。

Step2 對每一組鑒定數據單元，計算序列中的堿基位點的等位基因的先驗概率，根據概率選出最有可能變異的位點。

Step3 計算節點將SNP 結果輸出至硬盤并結束。

控制節點會持續詢問所有計算節點直到完成SNP鑒定，然后控制節點把所有結果根據劃分的地址順序進行合并。

2 實驗測試與結果分析

實驗測試有2 臺曙光A840R-G，4×AMD 6172 12核CPU，512G 內存，共4 個計算節點，操作系統為Linux version 2.6.32-431.el6.x86_64；Red Hat Enterprise Linux Server release 6.5，作業管理系統為LSF。

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

采用本文提出的并行優化模型在多節點上設置并行數依次為2 組、3 組、5 組、10 組、15 組、20 組、30 組、50 組、100 組，其中，每組并行數再分別計算樣本數1個、2 個、4 個、8 個、16 個、32 個、64 個、128 個，每個樣本的大小約4GB。通過上述計算，以運算時間和加速比[15]兩個指標來評價本文模型，統計運算時間如圖2 所示，加速比如圖3 所示。

圖2 樣本數和并行數對計算效率的影響

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

當樣本數小于16 個時，多組數據相互對比，計算時間差距并不明顯，使用并行化加速效果并不明顯，在實際應用中可以考慮使用傳統計算方法計算，避免不必要的計算誤差；但當樣本數為16 個至128 個時，采用并行加速模型可明顯增加提高計算效率。

由圖3 可看出隨著并行數的增加，加速比在并行數20 以內呈現明顯上升趨勢，說明增加并行數可有效地減小每個任務的計算時間，從而減少總體計算時間；加速比在并行數量30 到50 這個區間出現緩慢上升甚至持平的趨勢，這是因為當每個節點的計算時間縮小到一定程度后，SNP 鑒定的時間主要由相對固定耗時的I/O 操作效率所確定，并且已經相當接近，所以導致加速比不再明顯增加；當并行數在50 以上后，加速比將保持恒定。由圖2 和圖3 綜合也說明并行數量并不是并非簡單的越大越好，而是與樣本數有直接關系。

圖3 樣本數和并行數對多節點計算加速比的影響

3 結語

針對基因變異鑒定耗時較大的問題，通過MapReduce 框架模型將計算任務劃分，在多節點上實現并行加速計算，提出了一種針對多樣本的多節點并行化優化算法，在海量基因數據的計算與分析中具有廣泛適用性。

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