1 引言

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

2 隱馬爾科夫模型法

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

2.1 馬爾科夫鏈

考慮只取有限個或可數(shù)個狀態(tài)的隨機過程{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}成立，則稱此隨機過程為離散狀態(tài)馬爾科夫鏈。簡單的說，就是系統(tǒng)未來的狀態(tài)僅依賴于當(dāng)前狀態(tài)。一個馬爾科夫鏈的概率分布完全由它的初始分布P(X0)與轉(zhuǎn)移矩陣P=(pij)決定。

2.2 HMM基本原理

隱馬爾科夫模型HMM是由馬爾科夫鏈發(fā)展擴充而來的一種隨機模型。HMM可以被理解為一個雙重隨機過程，一個是不可觀察的（隱含的）狀態(tài)變化序列，另一個是由該不可觀察的狀態(tài)產(chǎn)生的可觀察符號序列。隱馬爾科夫模型形式描述如下：一個HMM模型是一個三元組M=(A，S，Q)，其中A是字母表，S是有限狀態(tài)集合，每個狀態(tài)可以釋放字母表中的字符。Q為概率集合，包括兩個部分：一是狀態(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的一個相繼狀態(tài)序列，X=（x1，x2，…，xL）是一個字符序列，按下述方式定義狀態(tài)轉(zhuǎn)換概率和字符釋放概率：

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

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

對于給定的路徑Π，可以按下面的公式計算出產(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模型時，用方框表示各個狀態(tài)，方框之間的連線表示狀態(tài)轉(zhuǎn)換。對于每個狀態(tài)，詳細(xì)地描述各個字符的釋放概率，而對于狀態(tài)之間的轉(zhuǎn)換，也給出相應(yīng)轉(zhuǎn)換動作發(fā)生的概率，即狀態(tài)轉(zhuǎn)換概率。表示DNA序列的HMM如圖1所示。

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

圖1 隱馬爾科夫模型

3 并行算法

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

給定一個字符序列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）對于每個i=0，1，…，L-1及每個l∈S，按下式進行遞歸計算

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

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

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

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

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

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

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

4 結(jié)束語

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