并行計(jì)算在生物序列比對中的應(yīng)用

郝靜 劉雅坤

摘 要：生物信息學(xué)作為一門綜合運(yùn)用分子生物學(xué)、數(shù)學(xué)和計(jì)算機(jī)等學(xué)科的理論和方法的交叉學(xué)科為闡明和理解海量DNA序列數(shù)據(jù)所包含的生物意義提供了可能。而海量的數(shù)據(jù)信息，使得生物序列比對計(jì)算需要耗費(fèi)大量的時(shí)間，本文中針對此問題，實(shí)現(xiàn)Windows API，OpenMP的并行算法，比較各自的加速比和最優(yōu)比對序列。

關(guān)鍵詞：生物序列比對 Windows API OpenMP

一、問題背景與提出

生物信息學(xué)最首要的任務(wù)就是從大量的生物信息數(shù)據(jù)中，提取有生物價(jià)值的信息，而序列比對是當(dāng)前最重要，最基本的方法，是指兩個(gè)或多個(gè)序列按字母比較， 盡可能確切地反映它們之間的相似和相異性，闡明序列之間的同源關(guān)系。在序列分析中， 將未知序列同已知序列進(jìn)行相似性比較是一種強(qiáng)有力的研究手段， 從序列的片段測定， 拼接， 基因的表達(dá)分析， 到 RNA 和蛋白質(zhì)的結(jié)構(gòu)功能預(yù)測， 物種親緣樹的構(gòu)建都需要進(jìn)行生物分子序列的相似性比較。

一般來說， 評價(jià)生物序列比對算法的標(biāo)準(zhǔn)有兩個(gè)： 一為算法的運(yùn)算速度， 二為獲得最佳比對結(jié)果的敏感性或準(zhǔn)確性。然而，隨著測序技術(shù)的發(fā)展，生物序列分析要求新的軟件方法和計(jì)算平臺，在分析大規(guī)模序列數(shù)據(jù)和解決復(fù)雜問題上，并行計(jì)算發(fā)揮著巨大的作用。

二、模型建立

序列比對是運(yùn)用某種特定的數(shù)學(xué)模型或算法， 找出兩個(gè)或多個(gè)序列之間的最大匹配堿基或殘基數(shù)， 比對算法的結(jié)果在很大程度上反映了序列之間的相似性程度以及它們的生物學(xué)特征。序列比對根據(jù)同時(shí)進(jìn)行比對的序列數(shù)目多少可分為雙序列比對和多序列比對，從比對范圍考慮也可分為全局比對和局部比對。

2.1 兩兩序列比對

兩兩序列比對，就是把兩條未知的序列進(jìn)行排列，通過字母的匹配，刪除和插入操作，使得兩條序列達(dá)到同樣長度，在操作的過程中，盡可能保持相同的字母對應(yīng)在同一個(gè)位置。

通常用得分矩陣描述序列兩兩比對， 兩條序列分別作為矩陣的兩維， 矩陣點(diǎn)記錄兩個(gè)維上對應(yīng)的兩個(gè)DNA序列的相似性分?jǐn)?shù)， 分?jǐn)?shù)越高則說明兩個(gè)DNA序列越相似。

2.2 序列比對的基本定義

在生物分子信息處理過程中， 將生物分子序列抽象為字符串， 其中的字符取自特定的字母表。例如，DNA 序列由四種核苷酸組成， 用“A”， “T”， “C”， “G”代表四種堿基， 其復(fù)雜度為 4。

生物序列比對可以看作字符串的比對，用如下的定義來描述序列的比對與相似性。

空位的引入是為了補(bǔ)償由于變異而產(chǎn)生的插入或缺失，每引入一個(gè)空位，比對的分值都會有所扣除：

其中， 表示空位的總罰分； 表示初始化一個(gè)空位的罰分； 表示空位延伸一個(gè)間隔的罰分。

Def4 全局最優(yōu)比對和局部最優(yōu)比對

對于兩條序列S和T的全局比對可以用序列 和 來表示，其中，

l） 和 分別是在S和T中加入一些空位而得到的序列；

2） ， 為序列 的長度。

將序列 和 相應(yīng)的位置進(jìn)行一一比較，其分值Score可用如下的公式來表示：

序列S和T的全局最優(yōu)比對：指在S和T的所有比對中相似性分值Socer最大時(shí)所對應(yīng)的比對。

（2）序列S和T的局部最優(yōu)比對：用 來表示，其中 分別為S和T的子序列，且 的Score分值是S和T中所有子序列比對分值的最大值。

三、局部最優(yōu)序列比對串行算法描述

3.1 基本思想

本文針對兩兩序列最優(yōu)局部比對的Smith-Waterman方法，利用初始條件和迭代關(guān)系得到一個(gè)得分矩陣，選取其中得分最高的子序列的末端，然后利用動態(tài)規(guī)劃的方法回溯尋找，得到局部最優(yōu)比對序列。

對于兩個(gè)長度分別為n和m 的DNA序列：

構(gòu)造得分矩陣 ，用來存放所有可能的對比結(jié)果。

初始條件：D（i，0）=D（0，i）=0

遞歸關(guān)系：

其中， ， 分別為在序列S和T中添加一個(gè)長度為x，y 的空位的罰分。

在得分矩陣D中找出 ： ，則 表示序列S和T的局部最優(yōu)比對分值。從 開始，按照從右下到

左上的方向，對矩陣D進(jìn)行回溯就可以求得局部最優(yōu)比對序列。

3.2 算法描述

為簡化問題，本文中假設(shè)序列S和T中分別只可加入一個(gè)空位，且空位間隔為1，則設(shè) 。

從而遞歸關(guān)系簡化如下：

四、某問題的并行算法描述

4.1 基于Windows API的多核并行算法的設(shè)計(jì)

Windows API多核并行算法利用WINAPI定義線程函數(shù)，然后主函數(shù)里創(chuàng)建線程，將線程函數(shù)導(dǎo)入創(chuàng)建的線程中運(yùn)行。具體步驟如下：

（1）主函數(shù)中利用CreateThread函數(shù)創(chuàng)建線程；

（2）定義線程函數(shù)，兩個(gè)線程分別調(diào)用，并根據(jù)線程號均分任務(wù)；

（3）為避免數(shù)據(jù)競爭，利用臨界區(qū)比較中間得分矩陣的最大值，得到最終結(jié)果；

（4）利用回溯在（3）得到的得分矩陣中找到局部最優(yōu)序列比對。

4.2 基于Open MP的多核并行算法的設(shè)計(jì)

OpenMP是一種面向共享內(nèi)存以及分布式共享內(nèi)存的多處理器多線程并行編程語言，其執(zhí)行模式采用Fork-Join的形式，當(dāng)程序執(zhí)行時(shí)，主線程生成（Fork）一組線程，隨著程序的執(zhí)行，在并行結(jié)束后，這組線程匯合（Join）成主線程。具體步驟如下：

（1）利用threadprivate指令將全局變量thread_xl復(fù)制到各個(gè)線程中；

（2）使用#pragma omp parallel for 將已知序列S均分到兩個(gè)線程中，實(shí)現(xiàn)循環(huán)并行化，分別得到兩個(gè)得分矩陣；

（3）利用#pragma omp parallel實(shí)現(xiàn)并行，求各線程得分矩陣中的最大元素及其位置；

（4）利用原子操作對兩個(gè)得分矩陣中的最大值進(jìn)行比較，選其較大者；

利用回溯在（4）得到的得分矩陣中找到局部最優(yōu)序列比對。

五、算法實(shí)現(xiàn)

5.1 Windows API

1.任務(wù)分配：

HANDLE hThread[numthread]；

int myid[numthread]；

InitializeCriticalSection（&cs）；

for（int i=0；i

{myid[i]=i；hThread[i]=CreateThread（NULL，0，API_getD，&myid[i]，0，NULL）；}

WaitForMultipleObjects（numthread，hThread，TRUE，INFINITE）；

DeleteCriticalSection（&cs）；

2.線程函數(shù)：

DWORD WINAPI API_getD（LPVOID arg） {

int threadID=*（（int*）arg）；

XL thread_xl={{0}，{0}，{0}，{0}}；

int p=1，n=5，m=6；

if（threadID==1）

{p=5；n=9；}

getD（p，n，m，thread_xl）；

D_Max（n，m，thread_xl）；

EnterCriticalSection（&cs）；

if（thread_xl.max>xl_2.max）

xl_2=thread_xl；

LeaveCriticalSection（&cs）；

return 0；}

5.2 OpenMP

1.并行計(jì)算得分矩陣

#pragma omp parallel for

for（int i=1；i<9；i++）{

for（int j=1；j<6；j++）{

thread_xl.score[i][j]=Score（s[i]，t[j]）；

int a=thread_xl.D[i-1][j-1]+thread_xl.score[i][j]；

int b=thread_xl.D[i-1][j]-W；

int c=thread_xl.D[i][j-1]-W；

thread_xl.D[i][j]=Max（a，b，c）；}}

2.計(jì)算最終得分矩陣中最大值以及位置

#pragma omp parallel

int myid=omp_get_thread_num（）；

D_Max（9，6，thread_xl）；

#pragma omp critical

if（thread_xl.max>xl_1.max）{xl_1=thread_xl；}

數(shù)值實(shí)驗(yàn)和結(jié)論

為了比較串行與多種并行算法的效率，以串長分別為9和6的序列為例，進(jìn)行數(shù)值試驗(yàn)，得到結(jié)果如下。

以串行計(jì)算時(shí)間為基礎(chǔ)，采用API，OpenMP算法使用兩個(gè)核計(jì)算的加速比分別為：0.430，0.667；加速比不高，這是因?yàn)獒槍Ρ纠械腟序列和T序列，序列長度短，在尋找最優(yōu)比對的過程中，計(jì)算量不大，所以串行也可以很快的完成計(jì)算；而并行在創(chuàng)建線程時(shí)消耗了大量時(shí)間，因此出現(xiàn)加速比小于1的情況。

得到的最優(yōu)序列比對分別為：串行與OpenMP：tcggc，t-ggc；API：ggc。結(jié)果的不同是由于這幾種并行計(jì)算算法的差異導(dǎo)致的。在具體分配任務(wù)時(shí)，API與MPI真正的將S序列分成兩段，得到兩個(gè)得分矩陣選擇含有最大得分的矩陣進(jìn)行回溯；而在OpenMP中，雖然將S序列分成兩段，但是得到的兩個(gè)得分矩陣是存放在一個(gè)大矩陣中進(jìn)行回溯。

從以上結(jié)果的分析中，可以看出，OpenMP算法的并行設(shè)計(jì)更適合此類問題。

為了比較不同迭代次數(shù)下多種計(jì)算方式的計(jì)算效率，在程序?qū)崿F(xiàn)中，增加隨機(jī)賦值已知序列與未知序列部分，增加其序列長度，得到計(jì)算結(jié)果如上所示。總體上，三種并行方式的加速比隨著迭代次數(shù)的增加是提高的，但受到計(jì)算量的限制，加速比仍未體現(xiàn)出并行計(jì)算在大計(jì)算量運(yùn)算中的優(yōu)勢，但是較上例已經(jīng)有所提高。

參考文獻(xiàn)

[1] 陳華.高性能并行計(jì)算.中國石油大學(xué)（華東）[M].2016

[2] 張法.基于Smith_Waterman算法的并行分而治之生物序列比對算法[J].中國科學(xué)：技術(shù)科學(xué).2004.34（2）：190-199

作者簡介：

郝靜（1994—），女，漢族，山東煙臺人，中國石油大學(xué)（華東）理學(xué)院，2013級本科生，數(shù)學(xué)與應(yīng)用數(shù)學(xué)專業(yè)。

劉雅坤（1995—），女，漢族，山東青島人，中國石油大學(xué)（華東）理學(xué)院，2013級本科生，數(shù)學(xué)與應(yīng)用數(shù)學(xué)專業(yè)。