李金光 唐友 李丹

基金項目：吉林省科技發(fā)展計劃項目；項目名稱：基于數(shù)據(jù)挖掘技術(shù)的全基因組選擇方法研發(fā)及云計算平臺體系構(gòu)建；項目編號：YDZJ202201ZYTS692。

作者簡介：李金光（1999— ），男，河北廊坊人，碩士研究生；研究方向：先進控制技術(shù)與控制系統(tǒng)集成。

*通信作者：李丹（1980— ），男，吉林吉林人，副教授，博士研究生；研究方向：電氣工程。

摘要：近年來，多核技術(shù)在基因測序方面有著廣泛的應(yīng)用。因此，基因測序的并行化成為目前國內(nèi)熱點。在基因測序領(lǐng)域，隨著高通量測序技術(shù)的迅猛發(fā)展和應(yīng)用越加廣泛，其無時無刻不在產(chǎn)生大量的數(shù)據(jù)。對此，串行程序進行單個程序運行已經(jīng)無法滿足社會對基因測序的需求。為了同時分析數(shù)千個樣本，測序需要快速、靈活和內(nèi)存高效地實現(xiàn)。文章認為，對一個名為ANGSD的程序進行OpenMP的程序并行的軟件優(yōu)化設(shè)計和可行性研究非常必要。

關(guān)鍵詞：OpenMP；下一代測序；多線程

中圖分類號：TP39 文獻標(biāo)志碼：A

0 引言

測序技術(shù)經(jīng)過多年發(fā)展擁有了3個階段成果，分別是第一代測序技術(shù)、下一代測序技術(shù)（NGS）和最新的第三代單分子測序技術(shù)。下一代測序技術(shù)（Next Generation Sequencing，NGS）又被叫作高通量測序技術(shù)，此技術(shù)一次測序可以檢測幾十萬到幾百萬條序列。由于此技術(shù)具有高通量、強穩(wěn)定、高準(zhǔn)確度等優(yōu)點，因此在商業(yè)上被廣泛應(yīng)用［1］。

現(xiàn)有的通用多樣本NGS分析程序的例子有單線程SAMtools （C）和多線程GATK（Java）以及ANGSD。這3個程序之間有許多不同之處，但ANGSD的關(guān)鍵優(yōu)勢在于：（1）允許與原始測序數(shù)據(jù)直接相關(guān)的多種輸入數(shù)據(jù)類型（文本堆積、二進制基因型可能性文件、VCF文件）；（2）允許用戶選擇多種方法中間分析（例如計算GL的不同方法）；（3）包括一組在任何其他軟件中都沒有實現(xiàn)的下游分析的實現(xiàn)。ANGSD中很多方法都是基于GLs。ANGSD支持4種不同的模型來計算GLs：（1）重新校準(zhǔn)的SOAPsnp模型；（2）原GATK模型；（3）SAMtools 1.16 +修改Maq模型；（4）特定類型的錯誤模型。這些GL模型中的排序錯誤率要么是固定的（由qscores獲得），要么是從數(shù)據(jù)中估計出來的。4個實現(xiàn)的GL模型假設(shè)為二倍體樣本。工作流程分為2個步驟：（1）ANGSD生成特定的分析輸入數(shù)據(jù)。（2）根據(jù)ANGSD輸出，使用二級關(guān)聯(lián)程序進行下游分析。對于簡單的測試，例如ABBA-BABA/D-statistic，次要程序可以是簡單的Rscript；對于計算密集型的方法，次要程序可以是多線程的c/c++程序［2］。OpenMP是多核系統(tǒng)最流行的并行編程之一。因此，本文就OpenMP并行在ANGSD軟件的應(yīng)用可行性做展開論述。

1 研究意義

20世紀70年代，Walter Gilbert和Frederick Sanger發(fā)明了人類史上第一臺測序儀，并通過測序儀探測到第一個基因組序列噬菌體X174，全長擁有5 375個堿基。

在Sanger基因測序技術(shù)的催動下，人們對生命的本質(zhì)產(chǎn)生極大的興趣。在此動力下，人們對生命科學(xué)的研究逐漸步入了基因組學(xué)的時代。40多年來，人們通過不斷努力已將測序技術(shù)發(fā)展到了可觀的境界，已經(jīng)由第一代發(fā)展到了第三代測序技術(shù)。

由Sanger發(fā)明的最初的測序方法被人們稱為第一代測序技術(shù)，雖然是第一代技術(shù)，但至今仍在被廣泛應(yīng)用。同時，它也有著無法掩蓋的缺陷。它進行一次測序只能得到一條長度為700～1 000個堿基的序列，這個缺陷成為在世界舞臺上更進一步的掣肘，無法滿足現(xiàn)代科學(xué)發(fā)展所需要的龐大生物基因序列。

高通量測序（High-Throughput Sequencing，HTS）彌補了第一代基因測序的缺陷，運行一次便可以同時獲得幾十萬至幾百萬條核酸分子的序列，因為在Sanger基因測序技術(shù)的基礎(chǔ)上進行了突破，所以高通量測序也被稱為新一代測序 （Next Generation Sequencing，NGS）或第二代測序。

第二代測序技術(shù)雖然對各個領(lǐng)域產(chǎn)生了極大的影響，但是也有缺點。第二代測序技術(shù)獲得單條序列長度很短，如果想得到更為準(zhǔn)確的基因序列信息，就不得不依賴于較高的測序覆蓋度和更準(zhǔn)確的序列拼接技術(shù)。由于技術(shù)問題，最終得到的結(jié)果中會存在一定量的錯誤信息，誤導(dǎo)研究人員的判斷。

為了防止結(jié)果中出現(xiàn)錯誤信息擾亂人的認知，科研人員在保證高通量測序的基礎(chǔ)上，發(fā)明了第三代測序技術(shù)，也稱為單分子測序技術(shù)（Single-molecule sequencing technology）。該技術(shù)在上一代測序技術(shù)上做出了突破，解決了單條序列長度過短的問題，如今能夠直接得到長度在數(shù)萬個堿基的核酸序列信息。

目前，基因測序技術(shù)在包括胚胎的植入前遺傳學(xué)診斷研究、法醫(yī)學(xué)的少量DNA測序、表觀遺傳學(xué)和考古學(xué)、物種進化演替過程等眾多領(lǐng)域得到廣泛應(yīng)用。

就當(dāng)前的市場勢態(tài)而言，第二代短讀長序列檢測技術(shù)在世界測序技術(shù)市場上一直占據(jù)著絕對的優(yōu)勢地位，而第三代測序技術(shù)也在近年來實驗中迅速發(fā)展。

若要進行ANGSD軟件的簡化提速，就需要對ANGSD軟件運行的流程進行分析：可以通過程序運行來確定ANGSD讀取數(shù)據(jù)、運算以及輸出數(shù)據(jù)的時間占比；還可以通過Pthread和OpenMP的兩種并行語言的對比進行分析，以此來判斷ANGSD軟件是否擁有簡化的可能性。

2 ANGSD數(shù)據(jù)處理流程

ANGSD是一個多線程程序。該程序可以計算各種匯總統(tǒng)計，利用下一代測序數(shù)據(jù)中的完整信息，通過直接處理原始測序數(shù)據(jù)或使用基因型可能性，進行關(guān)聯(lián)映射和群體遺傳分析。該程序不僅支持多種輸入格式，包括BAM和beagle基因型概率文件，而且允許用戶在現(xiàn)有方法的組合之間進行選擇，并可以執(zhí)行其他地方?jīng)]有實現(xiàn)的分析。此篇文章便對輸入為BAM基因型概率的文件進行研究。

本文通過利用gdb調(diào)試工具，來判斷程序函數(shù)的調(diào)用關(guān)系，得出圖1為ANGSD軟件的函數(shù)調(diào)用流程。命令如下：

gdb ./angsd

set args-GL 1-doGlf 2-b bam2.filelist-doMajorMinor 1 -doMaf 1

b angsd.cpp：main

b shared.cpp：main_analysis

set print pretty on

r

以上已經(jīng)得出ANGSD軟件運行處理的流程，接下來為了更直觀地表現(xiàn)出哪一部分在軟件中可以進行簡化更改，本文對程序選取不同規(guī)格的Bam文件進行測序，觀測在虛擬機以及服務(wù)器的運行速度。

從表1可以發(fā)現(xiàn)不同容量Bam文件分別在reads運行時間占CPU總時間的5%，5%，7.5%，7.8%，所以可以得出結(jié)論讀取Bam文件reads記錄時間約為總時間的7%。

從表2可以發(fā)現(xiàn)不同容量Bam文件分別在輸出beagle基因型概率文件時間占CPU總時間的22%，20%，11.4%，10.9%，所以可以得出結(jié)論輸出beagle基因型概率文件約為總時間的10%。

從表3可以發(fā)現(xiàn)不同容量Bam文件分別在計算基因型概率時間占CPU總時間的7%，10%，16%，19%，所以可以得出結(jié)論：計算基因型概率時間約為總時間的20%。

通過上述數(shù)據(jù)可以發(fā)現(xiàn)ANGSD軟件擁有進一步提速的可能性，接下來本文將對其軟件加速的方案進行探討分析。

3 兩種并行方式簡介

3.1 OpenMP并行

并行計算機可簡單分成共享內(nèi)存和分布式存儲器。共享內(nèi)存是指幾個內(nèi)核之間共用的一條存儲器。而對只注重計算、只需要對線程間關(guān)系做出最基本控制的應(yīng)用而言，OpenMP十分適用。OpenMP屬于半自動多核并行，它介于手工并行和全自動并行之間。它通過兩方面實現(xiàn)并行化：（1）依靠編譯器來實現(xiàn)多核并行化；（2）需要程序員指定并行的代碼段或分配的粒度。

通常情況下，人工介入的部分只需要加上注釋或是標(biāo)注，相比手工并行，半自動并行更加通俗易懂［3］ 。

OpenMP是一個應(yīng)用程序接口（API），由一組主要的計算機硬件和軟件供應(yīng)商聯(lián)合定義。OpenMP為共享內(nèi)存并行應(yīng)用程序的開發(fā)人員提供了一個可移植、可伸縮的模型。該API在多種體系結(jié)構(gòu)上支持C/C++和Fortran。

OpenMP還給出了對并行算法的高層抽象說明，更適用于在多核CPU計算機上的并行程序?qū)崿F(xiàn)。它在多核心CPU的機器中采用OpenMP程序方式實現(xiàn)了并行計算，編譯時通過在程序中附加的pragma命令，可以自動地對程序并行處理，使OpenMP減少了對并行程序?qū)崿F(xiàn)的困難度和復(fù)雜性［4］。

3.2 Pthreads并行

POSIX線程（POSIX Threads，常被縮寫為Pthr-eads）是POSIX的線程標(biāo)準(zhǔn)，定義了創(chuàng)建和操縱線程的一套API。Pthreads作為手工多核并行，程序員需要在程序中手動調(diào)用API來實現(xiàn)并行化。

學(xué)術(shù)界將符合POSIX 線程標(biāo)準(zhǔn)的庫稱為Pthreads，一般應(yīng)用于Unix-likePOSIX 系統(tǒng)，例如Linux系統(tǒng)。但是Windows上的實現(xiàn)也存在，例如直接使用Windows API實現(xiàn)的第三方庫pthreads-w32；而利用Windows的SFU/SUA子系統(tǒng)，則可以使用微軟提供的一部分原生POSIX API 。

Pthreads API中大致共有100個函數(shù)調(diào)用，全都以“pthread_”開頭，并可以分為4類：

（1）線程管理，例如創(chuàng)建線程、等待（join）線程、查詢線程狀態(tài)等。

（2）互斥鎖（Mutex）：創(chuàng)建、摧毀、鎖定、解鎖、設(shè)置屬性等操作。

（3）條件變量（Condition Variable）：創(chuàng)建、摧毀、等待、通知、設(shè)置與查詢屬性等操作。

（4）使用了互斥鎖的線程間的同步管理。

不但如此，Pthreads還擁有一套貼近底層的API，可以用來實現(xiàn)對線程的精準(zhǔn)控制，例如創(chuàng)建/取消/同步某一個線程，設(shè)置各個線程的屬性，操作互斥鎖和條件變量，進行線程間同步等功能。

3.3 OpenMP和Pthreads共用模型

Pthreads和OpenMP都是多線程編程模型，這兩種并行語言通常采用如圖2所示的fork-join執(zhí)行。fork遞歸地將任務(wù)分解為較小的獨立子任務(wù)，直到它們足夠簡單以便異步執(zhí)行。join將所有子任務(wù)的結(jié)果遞歸地連接成單個結(jié)果，或者在返回void的任務(wù)的情況下，程序只是等待每個子任務(wù)執(zhí)行完畢。

4 ANGSD軟件兩種并行方式區(qū)別

ANGSD作為高通量測序的一個軟件在基因測序領(lǐng)域有著不可忽視的作用。該程序使用的是Pthreads并行，而本文進行對比的是OpenMP并行，雖然都是針對共享內(nèi)存而編寫的API，但是這兩者存在著本質(zhì)上的區(qū)別。Pthreads要求程序員直接表明對一個線程的執(zhí)行情況；而OpenMP則相對而言要求寬松很多，它允許程序員只需簡短地說明哪條程序需要并行執(zhí)行，屆時將由編譯器和所執(zhí)行時的操作系統(tǒng)來決定具體的線程去執(zhí)行某個任務(wù)。這種差異也證明了為什么共享內(nèi)存程序同時具有兩種完全不同的API。不過在Pthreads所有線程的行為中的細節(jié)都必須由程序員確定。與此相反，OpenMP最主要的優(yōu)勢之一，便是它的設(shè)計使編寫者能夠逐步地并行現(xiàn)有的串聯(lián)型程式，而不會再從零開始寫并行程序。這可以極大地簡化編程煩瑣的程序，因此把Pthreads更換為OpenMP并行進行簡化是可行的［5］ 。

5 結(jié)語

本文探討的內(nèi)容對ANGSD軟件加速的可能性進行了初步的分析，ANGSD軟件執(zhí)行的流程圖以及Bam文件在虛擬機以及服務(wù)器的運行也證明了其軟件是有加速的可行性的。對比Pthreads和OpenMP的優(yōu)缺點也可以說明OpenMP相對Pthreads而言具有更為簡化的編程方式，因此接下來的工作就可以對ANGSD軟件進行程序的更改，提高其運行的速度。

參考文獻

［1］黃芝準(zhǔn).組學(xué)大數(shù)據(jù)環(huán)境下的基因信息并行處理與分析方法研究［D］.合肥：中國科學(xué)技術(shù)大學(xué)，2017.

［2］DURVASULA A，HOFFMAN P J，KENT T V，et al. Angsd-wrapper：utilities for analysing next-generation sequencing data.［J］. Molecular ecology resources，2016（6）：1449-1450.

［3］戴晨，陳鵬，楊冬蕾，等.面向多核的并行編程和優(yōu)化研究［J］.計算機應(yīng)用與軟件，2013（12）：198-202，279.

［4］產(chǎn)院東.基于多核和眾核平臺的并行DNA序列比對算法［D］.濟南：山東大學(xué)，2019.

［5］段皞一.MPI、OpenMP、Taichi并行編程語言探究［J］.電子元器件與信息技術(shù)，2022（4）：123-134.

（編輯 王永超）

Abstract： In recent years， multi-core technology has been widely used in gene sequencing， so the parallelization of gene sequencing has become a hot topic in China. In the field of gene sequencing， with the rapid development and application of high-throughput sequencing technology， a large amount of data is generated all the time， and the single program operation of this serial program cannot meet the needs of society for gene sequencing. In order to analyze thousands of samples simultaneously， we need fast， flexible and memory-efficient implementations. OpenMP program parallelism for a program called ANGSD is very necessary.

Key words： OpenMP; next-generation sequencing; multithreading