基于HMM的基因識別并行計(jì)算

1 引言

20世紀(jì)90年代以來，伴隨著各種基因組測序計(jì)劃的展開和分子結(jié)構(gòu)測定技術(shù)的突破，數(shù)以百計(jì)的生物學(xué)數(shù)據(jù)庫如雨后春筍般迅速出現(xiàn)和成長。如何利用這些不斷爆炸性增長的有關(guān)生物分子的原始數(shù)據(jù)，有效解決基因識別問題顯得越來越迫切。最初的基因分析方法是進(jìn)行簡單的核苷酸統(tǒng)計(jì)，而后加上剪切保守位點(diǎn)的檢測。以后采用了人工神經(jīng)網(wǎng)絡(luò)、隱馬爾科夫模型（HMM）[1，2]等先進(jìn)的信息處理和分析技術(shù)，提高基因識別的準(zhǔn)確率。但由于生物信息數(shù)據(jù)量巨大，傳統(tǒng)的串行算法往往無法處理或難以在滿意的時(shí)間內(nèi)得到結(jié)果。本文針對基因序列的識別，討論隱馬爾科夫模型分析算法的并行算法設(shè)計(jì)和并行效果分析。

2 隱馬爾科夫模型法

隱馬爾科夫模型[3]（Hidden Markov Models，HMM）是一種概率論模型，這種方法已經(jīng)成功應(yīng)用于多個(gè)領(lǐng)域，如語音識別、光學(xué)字符識別等。HMM在生物信息學(xué)領(lǐng)域中也有著重要的應(yīng)用，如序列分析、基因識別等。目前，基因識別的HMM方法也大致可以分為兩類，一類為按照內(nèi)容搜索的方法，通過核苷酸和三聯(lián)密碼子等在編碼區(qū)的分布規(guī)律來界定蛋白質(zhì)的編碼區(qū)；另一類為按照信號搜索的方法，通過編碼區(qū)周圍的信號界定蛋白質(zhì)編碼區(qū)。

2.1 馬爾科夫鏈

考慮只取有限個(gè)或可數(shù)個(gè)狀態(tài)的隨機(jī)過程{Xn，n=0，1，2，…}，假設(shè)對一切狀態(tài)i0，i1，…，in-1，i，j和一切n≥0，有P{Xn+1=j | Xn=i，Xn-1=in-1，…，X1=i1，X0=i0} = P{Xn+1=j | Xn=i}成立，則稱此隨機(jī)過程為離散狀態(tài)馬爾科夫鏈。簡單的說，就是系統(tǒng)未來的狀態(tài)僅依賴于當(dāng)前狀態(tài)。一個(gè)馬爾科夫鏈的概率分布完全由它的初始分布P(X0)與轉(zhuǎn)移矩陣P=(pij)決定。

2.2 HMM基本原理

隱馬爾科夫模型HMM是由馬爾科夫鏈發(fā)展擴(kuò)充而來的一種隨機(jī)模型。HMM可以被理解為一個(gè)雙重隨機(jī)過程，一個(gè)是不可觀察的（隱含的）狀態(tài)變化序列，另一個(gè)是由該不可觀察的狀態(tài)產(chǎn)生的可觀察符號序列。隱馬爾科夫模型形式描述如下：一個(gè)HMM模型是一個(gè)三元組M=(A，S，Q)，其中A是字母表，S是有限狀態(tài)集合，每個(gè)狀態(tài)可以釋放字母表中的字符。Q為概率集合，包括兩個(gè)部分：一是狀態(tài)轉(zhuǎn)換概率fkl，k，l∈S，表示從狀態(tài)k轉(zhuǎn)化到狀態(tài)l的概率；二是字符釋放概率，記為ek(b) (k∈S，b∈A)，表示在狀態(tài)k下釋放出字符b的概率。令路徑Π=(π1，π2，…，πL )是模型M的一個(gè)相繼狀態(tài)序列，X=（x1，x2，…，xL）是一個(gè)字符序列，按下述方式定義狀態(tài)轉(zhuǎn)換概率和字符釋放概率：

fkl = p(πi = l|πi-1 = k)

ek(b) = p(xi=b|πi= k)

對于給定的路徑Π，可以按下面的公式計(jì)算出產(chǎn)生序列X的概率：

P（X|Π）= fπ0，π1 eπi (xi)fπi，πi+1

這里，令π0為起始狀態(tài)，πi+1為終止?fàn)顟B(tài)。

在表示或分析HMM模型時(shí)，用方框表示各個(gè)狀態(tài)，方框之間的連線表示狀態(tài)轉(zhuǎn)換。對于每個(gè)狀態(tài)，詳細(xì)地描述各個(gè)字符的釋放概率，而對于狀態(tài)之間的轉(zhuǎn)換，也給出相應(yīng)轉(zhuǎn)換動(dòng)作發(fā)生的概率，即狀態(tài)轉(zhuǎn)換概率。表示DNA序列的HMM如圖1所示。

對生物序列而言，HMM的字符就是20個(gè)字母的氨基酸或4個(gè)字母的核苷酸。編碼蛋白質(zhì)的原始DNA序列，在生物的進(jìn)化過程中會受到自然環(huán)境和各種因素的影響，使翻譯出的蛋白質(zhì)序列[4]經(jīng)歷突變、遺失或引入外援序列等變化，最后按不同的進(jìn)化路徑分化，形成多種功能相近的蛋白質(zhì)。因此，可以把這些蛋白質(zhì)看作由一個(gè)基本蛋白質(zhì)序列經(jīng)過插入、刪除或替換了某些氨基酸殘基而形成。這個(gè)過程可以用HMM來表示。一個(gè)訓(xùn)練好的模型可以代表有共同特征的蛋白質(zhì)序列。HMM用于分析蛋白質(zhì)序列的原理是分析蛋白質(zhì)產(chǎn)生不同序列的概率，對于與模型相符合的序列，能以較大的概率產(chǎn)生。

圖1 隱馬爾科夫模型

3 并行算法

對于給定一個(gè)隱馬爾科夫模型M=(A，S，Q)和一個(gè)字符序列X（即基因序列），在M中尋找產(chǎn)生該序列的最優(yōu)路徑Π*，該路徑從起始狀態(tài)出發(fā)，結(jié)束于終止?fàn)顟B(tài)，在路徑中的每一個(gè)狀態(tài)都選擇釋放一個(gè)字符，使P(X|Π*)最大。這是基因識別中常用的一個(gè)方法，這里我們設(shè)計(jì)采用并行算法來求解HMM的最優(yōu)路徑問題。

給定一個(gè)字符序列X=（x1，x2，…，xL），以vk(i)代表序列前綴（x1，x2，…，xL）終止于k（k∈S，1≤i≤L）的最可能路徑的概率。求解過程如下：

（1）初始化 vbegin（0）=1

k≠begin vk(0) = 0

（2）對于每個(gè)i=0，1，…，L-1及每個(gè)l∈S，按下式進(jìn)行遞歸計(jì)算

vl(i+1) = el(xi+1)max{vk(i)fkl} k∈S

（3）最后，計(jì)算序列X終止于狀態(tài)“end”最可能的路徑概率，即P(X|Π*)的值

P(X|Π*) = max{vk(L)fk，end} k∈S

在實(shí)現(xiàn)中我們將隱馬爾科夫模型使用一顆狀態(tài)空間樹及一個(gè)字符釋放概率矩陣聯(lián)合表示。如圖2所示。

圖2 HMM的聯(lián)合表示

采用并行深度優(yōu)先搜索技術(shù)，在每一個(gè)前向分支處啟動(dòng)一個(gè)新的進(jìn)程，并行的計(jì)算多計(jì)算分支。在單CPU的情況下，算法的時(shí)間復(fù)雜度為O(L|S|2)，在具有N個(gè)計(jì)算節(jié)點(diǎn)的情況下，算法的時(shí)間復(fù)雜度為O(L|S|2/N)。

在理想的情況下，并行算法在理論上的加速比與計(jì)算節(jié)點(diǎn)數(shù)成正比。在大型基因結(jié)構(gòu)識別的問題域中，為實(shí)現(xiàn)并行計(jì)算而產(chǎn)生額外的啟動(dòng)、通信等時(shí)間與有效計(jì)算時(shí)間相比基本可以忽略，可近似達(dá)到理想加速比。

4 結(jié)束語

中國科學(xué)院院士張春霆指出生物信息學(xué)是生物學(xué)的核心和靈魂，數(shù)學(xué)與計(jì)算機(jī)技術(shù)則是它的基本工具。只有將并行計(jì)算研究和基因識別的理論研究有效聯(lián)系起來，在研究蛋白質(zhì)結(jié)構(gòu)預(yù)測與分析的方法基礎(chǔ)上，結(jié)合并行計(jì)算技術(shù)的特點(diǎn)，設(shè)計(jì)一系列的高效并行實(shí)現(xiàn)技術(shù)，實(shí)現(xiàn)高效、快速的基因識別，生物信息計(jì)算才能得到更快的發(fā)展。