Ｍｏｔｉｆ識別算法簡介及軟件性能研究

摘 要：Motif在轉錄和后轉錄水平的基因表達調控中起著重要的作用。目前，識別Motif的算法和相應的軟件已有不少，但是卻鮮有對各種算法及軟件性能共同評測的研究和報告。介紹了算法的分類以及三種常見的Motif識別算法Wordup，MM和Gibbs采樣，并對AlignACE，MEME，MotifSampler，Weeder等13種Motif尋找軟件進行性能比較分析。通過生物學意義的研究和性能比較結果可以得出：由于唯有Weeder算法考慮了Motif保守核心位置，因而它在各種軟件中識別效果較好；大部分算法只考慮簡單而且短的Motif，所以各種軟件對酵母菌這種單細胞生物的Motif識別性能比多細胞生物要高。

關鍵詞：Motif； Wordup； MM； Gibbs采樣

中圖法分類號：TP301．6 文獻標識碼：A 文章編號：1001-3695(2006)10-0066-04

Introduction of Algorithms and Performance Research of

Softwares for Motif Discovery

ZHU Ji1，2， YANG Hua1，2， NIU Beifang1，2， LANG Xianyu1，2， LU Zhonghua1， CHI Xuebin1

(1.Supercomputing Center， Computer Network Information Center， Chinese Academy of Sciences， Beijing 100080， China; 2.Graduate School， Chinese Academy of Sciences， Beijing 100049， China)

Abstract:Motif plays a key role in the geneexpression regulating on both transcriptional and posttranscriptional levels. Nowadays there are several algorithms and softwares on detecting Motif， but， however， there is few papers on comparing the performance of these algorithms and softwares. This paper comes up with this background to introduce the classification of the algorithms in general and three common algorithms: Wordup， MM， Gibbs samplingin details. And a performance comparison is made on the thirteen softwares for Motif detecting such as AlignACE， MEME， MotifSampler， Weeder， etc. Based on the biological research and the performance report， this paper ends with a conclusion that Weeder is the most effective one of these softwares， for it is the only algorithm that takes account of the conserved core positions of Motifs; Most algorithms only consider simple and short Motifs， so their Motif detecting performance on monadic yeast is significantly higher than on metazoans.

Key words:Motif; Wordup; MM(Mixture Model); Gibbs Sampling

基因非編碼區的一個主要研究方向是對Motif的研究。因為轉錄和后轉錄水平，其基因的表達在很大程度上受到一些Motif的控制。它們本質上是一些比較短的DNA序列，這些序列一般均處在受調控基因的上游區域，轉錄因子可識別這些Motif并與之結合，從而調節DNA的代謝和轉錄；或者由RNA結合蛋白識別并與之結合，從而影響RNA的修飾、定位、翻譯和降解。因此，分析和識別Motif及了解它們的功能對于理解和解釋整個基因組行為的意義重大。

Motif的分析主要涉及三類問題：①在給定基因組序列中尋找已知的Motif；②在一系列共表達或者共調控基因的上游區域中發現未知的Motif；③尋找由一個已知轉錄因子調控的未知基因。本文主要討論第二類問題，即在一系列共表達基因的啟動子區域中探測新的Motif，通過分析和提取DNA序列特征來識別Motif。

一般原核生物的Motif特征比較明顯，容易識別。但是真核生物的Motif相對復雜，其Motif長度和空間分布變化較大，出現沒有固定的位置，相同蛋白質因子作用的結合位點也存在差異，這給識別Motif帶來了很大的困難。因此，要設計一個能識別所有Motif的方法幾乎是不可能的，而針對不同的生物和不同特點的Motif，出現了很多算法和軟件。

1 Motif尋找算法

1．1 算法的分類

根據算法搜索策略的不同，研究Motif的計算方法主要分為兩大類：①確定性的方法，即基于字串的方法，也稱單詞數數法，它包括簡單字串列舉法（YMF）、模式驅動列舉法（Wordup[2]，Oligo/Dyadanalysis，QuickScore）、樣本驅動列舉法（MOPAC）、輪廓列舉法（Consensus）、后綴樹法（Weeder）、不匹配（前綴）樹法（Mitra）等。單詞數數法是基于上游序列中的低聚核甘酸的頻率分析，將一個單詞出現的次數與其期望次數進行比較來衡量超代表性，然后將幾個相似的單詞組合起來形成一個Motif。②似然說的方法。它又稱為概率序列模型的方法，包括EM算法[1]（MEME，Improbizer）、Gibbs采樣算法[3，4]（AlignACE，ANNSpec，GLAM，MotifSampler，SeSiMCMC）等。它是一種近似算法，這類算法首先對Motif的信息進行某種近似描述，然后通過不斷迭代的過程對Motif信息進行調整優化，直至滿足迭代終止條件。

1．2 三種常見的算法

有許多方法可以檢測一條序列中重復的序列模式，這些方法經過改進后可用于從一組序列中提取相同或者相近的短序列，作為候選的Motif。這些方法所依據的主要原理是，Motif比隨機出現的序列片段具有更高的出現頻率。

1．2．1 Wordup算法

Wordup算法是在DNA序列中找出具有顯著統計特性的單詞。它是基于一級馬爾可夫鏈的分析方法，其允許對一組已知在功能上有聯系的（如啟動子區域、內含子等）且沒有進行過比對的序列進行分析，選出共同的具有特定生物學功能的單詞或者子序列。

假設有一條長度為n的序列S，S＝s1s2…sn（si=A，C，G，T；i=1，2，…，n），它含有n－w＋1個長度為w的寡聚核苷酸單詞，表示為tj＝sjsj＋1…sj＋W－1，j＝1，2…n－w＋1，寡聚核苷酸實際上是短的DNA片段，長度為w的寡聚核苷酸單詞的個數為4W個?，F在考慮有N條核酸序列，S1，S2，…，SN，令Li表示各序列的長度。假設寡聚核苷酸的分布符合泊松分布，用πi(tk)表示長度為w的第k種寡聚核苷酸單詞在序列Si中至少出現一次的概率（k＝1，2，…，4W），則

πi(tk)=1-e-μik（1）

其中μik為在第i條序列中長度為w的第k種單詞期望出現的平均次數，它可表示為

μik=qi(tk)×(Li-w+1)（2）

式(2)中qi(tk)是寡聚核苷酸單詞tk在序列i中的出現概率。令f(ux，y)和f(uz)分別表示tk中兩個相鄰位x，y二核苷酸和單核苷酸z在每一條序列中的出現概率。該方法在計算寡聚核苷酸單詞出現頻率時考慮了核苷酸間的相互影響，因為在DNA中有一個相當普遍的現象，就是在核苷酸使用上具有相當強的二核苷酸偏性。這是一個典型的一階馬爾可夫鏈模型，可以得出

qi(tk)=f(μ12)f(μ23)…f(μw-1，w)/f(μ2)…f(μw-1)（3）

而tk在不同序列中期望出現的次數為

∏(tk)=∑iπi(tk)（4）

假設一個序列集合中有N條DNA序列，長度為w的第k種單詞在不同序列中實際出現的次數由下式給出：

P(tk)=∑ipi(tk)（5）

其中，如果序列Si中存在tk ，則pi(tk)=1；否則pi(tk)=0。

根據基本的χ2統計可以計算出第k種單詞tk的統計學顯著性：

χ2(k)=(P(tk)-∏(tk))2/∏(tk)（6）

對于給定的χ2值及已知的分布函數P(χ2)，可以計算出超過這個值的長度為w的單詞個數：

N(χ2)=[1-P(χ2)]×4w（7）

1．2．2 MM算法

MM（Mixture Model）算法是EM（Expectation Maximization）算法的一種改進。該算法主要解決的問題是在一系列不知其Motif位置信息和特征矩陣（位置頻率矩陣）的共調控序列中，如果存在共同的Motif，確定Motif的位置和對應的特征矩陣。其基本思想在于Motif具有保守性，且有對應的特征矩陣，在不斷迭代的過程中只有當兩者適應時最大似然函數值才能達到最大。

MM算法的基本步驟是：先對序列集建立二元有限混合模型，再運用最大似然估計法對模型的參數值進行估計。

設序列集為S={S1，S2，…， SN}，取所有Sq中一切長度為w的子序列構成新的序列集X={X1，…，Xn}。由于Xi可能是Motif序列，也可能是非Motif序列（背景序列），因此X可以看成是一個隨機向量。當Xi（i=1，2，…，n）是Motif時，可用每個位置上堿基出現的概率來表征，即字母表A={a1，a2，…，an}中每個字母as出現在序列中的位置j的概率為fjs（j=1，2，…，W，s=1，2，…，L），再定義fj=（fj1，fj2，…，fjL）；當Xi為背景序列時可以用概率分布表示為f0=（f01，f02，…，f0L）。二元有限混合模型是建立在序列集合X的基礎上的，它可描述為以概率λ1選取模式序列，或以概率λ2=1-λ1選取背景序列構成的原序列。因此簡單地說，二元即表示二組，即Motif序列組和背景序列組。模型中待估計的參數為λ=(λ1，λ2)和θ=(θ1，θ2)，其中θ=(f1，f2，…，fw)，θ=f0。

定義建立指示變量：

Z=(Z1，Z2，…，Zn)

Zi=(Zi1，Zi2)

Zik=1 Xi來自第k組

0 其他（k=1，2)

當Xi為Motif時，Zi1=1，Zi2=0；而當Xi為背景序列時，Zi1=0，Zi2=1。顯然Z也為隨機變量，Zik=1說明Xi具有分布p(Xi|θk)。

模型的似然函數可表示為變量X和Z的聯合概率分布：

L(θ，λ|X，Z)=p(Xi，Zi|θ，λ)（8）

log L(θ，λ|X，Z)=log p(Xi，Zi|θ，λ)=

∑ni=1∑2k=1(Zik|X，θ，λ)log(p(Xi|θk)λk)（9）

其中p(Xi|θ1)=∏Wj=1∏Ls=1fI(s，Xij)js，p(Xi|θ0)=∏Wj=1∏Ls=1f I(s，Xij)0s。

Xij表示序列Xi的第j個位置的字母，I(s，x)是一個指示函數，其定義為

I(s，x)=1 x=as

0 其他

MM算法采用期望值最大來求解最大似然估計。首先假定Motif序列在原序列中的初始位置，計算特征矩陣；其次用這個特征矩陣重新計算初始位置；最后這兩個過程不斷迭代直至收斂或達到最大迭代次數。

（1）初始化模型參數，記為θ(0)和λ(0)，分別取0~1之間的隨機值，令迭代次數t=0。

(2)根據模型參數θ(i)和λ(i)，求Z(i)。首先計算出式（9）的期望值：

E[log L(θ，λ|X，Z)]=∑ni=1∑2k=1Z(i)iklog P(Xi|θk)+

∑ni=1 ∑2k=1Z(i)iklog λk（10）

其中Z(i)ik=E[Zik|X，θ(i)，λ(i)]=p(Xi|θ(i)k)λ(i)k∑2k=1p(Xi|θ(i)k)λ(i)k（11）

(3)根據Z(i)估計模型的參數θ(i+1)和λ(i+1)，使式（10）的值最大。

因為式（10）“+”號前后兩部分分別只與一個參數有關，因此可以對它們分別求最大值，對于右邊即求

arg maxλ∑ni=1∑2k=1Z(i)iklog λk

可得λ(i+1)k=∑ni=1Z(i)ikn。

對于式（10）“+”號左邊有

θ(i+1)k=arg maxθk∑ni=1Z(i)iklog p(Xi|θk)=

∑Wj=1∑Ls=1log f(i)js∑ni=1Z(i)ikI(s，Xij)(12)

定義c(i)0s=∑ni=1∑Wj=1Z(i)i2I(s，Xij)，c(i)js=∑ni=1EiZ(i)i1I(s，Xij)。其中，c(i)0s即字母as出現在背景序列的期望次數，c(i)js即字母as出現在Motif中位置j的期望次數；Ei是一個衰減系數，在0~1間變化，表示序列Xi是否是序列模式，若是，Ei不斷減少至0，否則Ei為1。這使得本算法可以從共調控序列中盡可能找到所有的Motif。

轉變后的式（12）為

θ(i+1)=(θ(i+1)1，θ(i+1)2)=(f(i+1)1，f(i+1)2，…，f(i+1)W，f(i+1)0)=arg maxθ∑Wj=0∑Ls=1c(i)jslog f(i)js(13)

利用拉格朗日算法可以解出

f(i+1)js=c(i)js∑Ls=1c(i)js（14）

(4)令t=t+1，重復(2)、(3)，直至θ前后兩次的差值小于用戶定義的閾值或達到最大迭代次數。

在經過一個完整的迭代過程后，可以根據θ的值得到序列模式的一致性序列（每一位取出現頻率最高的堿基）；根據Zi1的大小確定模式序列及其起始位置。

1．2．3 Gibbs采樣算法

Gibbs 采樣算法是一種特殊的馬爾可夫鏈蒙特卡羅方法（MCMC），該算法最早是用于蛋白質序列中的序列模式識別。后來將Gibbs采樣整合進貝葉斯模型并應用于多重序列比較，獲得了較好的結果。采用Gibbs采樣算法的軟件有MotifSampler，AlignACE，BioProspector，Gibbs Motif Sampler等。Gibbs 采樣算法識別Motif的基本原理是，通過隨機采樣不斷更新Motif模型和在各條序列中的出現位置以優化目標函數，當滿足一定的迭代終止條件時就得到了最終的候選Motif。

這里假設每類Motif在每條序列中出現且僅出現一次?；镜腉ibbs采樣具體的算法過程如下：

(1)初始化。包括Motif模型和背景模型的建立。采用位置頻率矩陣（PSFM）來表示Motif模型。首先隨機生成n個Motif在各條序列中的起始位點，記為SS={ssi}，i=1，2…，n。然后根據起始位點和定義的長度W得到候選Motif，并根據這些候選Motif建立位置頻率矩陣PSFM。背景模型通常采用獨立性模型，記為B={qA，qG，qC，qT}，表示各堿基在背景序列中的出現頻率。

(2)更新。從輸入序列集中順序選取一條序列Si，序列長度為Li(i=1，2…，n)，從SS中刪除屬于該序列的起始位點，根據變化后的SS重新計算PSFM。然后分別根據Motif模型和背景模型計算Si中所有可能的候選Motif，即Si[j，j+W-1]（j=1，2…，Li-W+1）的得分Scoree(i，j)和Scoreb(i，j)，也就是計算序列Si中從第j位到第j+L-1位的序列片段是Motif的可能性及是背景序列的可能性。這里：

Scoree(i，j)=Scoree(Si[j，j+W-1]|M)=

∏W-1t=0PSFM(t，Si[t+j，t+j])（15）

Scoreb(i，j)=Scoreb(Si[j，j+W-1]|B)=∏W-1t=0qSi[t+j，t+j]（16）

其中PSFM(i，b)是指Motif位置i上字母b的出現概率。

(3)采樣。計算兩種得分的比值ri，j=Scoree(i，j)Scoreb(i，j)，并對每一個j選取新的候選Motif，即以較大的概率選取比值較高的候選Motif，將其起始位點加入到SS中。根據新的起始位點集SS和Motif長度得到候選Motif并計算得分F：

F=∑Wi=1∑Lk=1piklogpikqak（17）

其中pik=PSFM(i，ak)， ak為字母表A={a1，a2，…，aL}中的第k個字母）。若得分F大于前一次得分，則轉到(2)，繼續迭代；否則重復(3)，直至重復次數大于一個預設定值；若所有序列均處理完，則轉到(4)。

(4)終止。若得分連續多次沒有改進或達到最大迭代次數，則終止程序。

1．3 四種算法的性能評價

Wordup算法屬于確定性的方法，這類方法非常簡單；但缺點是具有非常復雜的計算復雜度，只適合搜索短的Motif。MM算法和Gibbs采樣算法屬于似然說的方法，這類算法具有較低的計算復雜度，適合在大空間中搜索解；其缺點是不能保證得到問題的最優解。確定性算法的另一個缺點是嚴格要求Motif的連續性，即要求這個Motif中每個位置均是嚴格保守的，不能有替換發生。解決這一問題的方法是可以識別兩端保守，中間有固定長度的不保守序列Motif。該方法的思想是分別計算兩端保守區域的統計顯著性，然后計算給定中間不保守序列長度的條件下這兩端保守區域的聯合概率，從而計算其統計顯著性。但這種新方法也有明顯的缺點，就是搜索到的Motif內部不能有變化，兩端的保守區域不能有堿基替換。

Wordup算法中如果Motif出現位置與具有顯著統計意義的單詞重疊，那么其χ2值會比期望值要高得多。解決這個問題可以用下面的迭代過程對每一條序列進行處理。在第一次迭代中，對所有長度為W且χ2大于一定值的寡聚核苷酸單詞按χ2值從大到小進行排序。假設共有M個這樣的單詞，將后面的M－1個單詞與第一個（具有最大χ2值）單詞進行比較，那些與第一個單詞在位置上有重疊的單詞將被去掉。后面的M-1個單詞的χ2值將被重新計算，并且根據新的χ2值進行排序。同樣的過程重復M-1次，直到所有M個具有顯著統計意義的單詞不再相互重疊為止。

MM算法是EM算法的改進。它解決了EM嚴重依賴初始化條件，容易陷入局部極值而得不到最優解的問題。該算法的每個子序列均作一步EM迭代，保留似然最大的那個Motif模型作進一步EM迭代直到收斂，記錄對應Motif的位置，然后再反復這個過程直至得到優化結果。

Gibbs采樣算法有兩個缺點：①它從序列Si中隨機選取一個起始位置ssi，這個位置使得進度很慢，且不能確定最后能達到一個好的結果；②該算法計算每個起始位點ssi的F值，但是大部分位置不能達到很好的收斂效果，所以這些F值是沒有意義的。用螞蟻群聚最佳化算法可以找到每條序列更好的候選位置，這樣可以大大提高Gibbs采樣算法的效率和質量。

綜上所述，相對于Wordup算法，MM和Gibbs采樣算法確定性不如前者，但運行速度大大提高，而且很多實際應用均證明了這類得到的近似解基本上能滿足解決問題的需要。

2 13種軟件的性能比較

Motif尋找的軟件很多，它們尋找Motif的性能也不一樣，在此對13個軟件（AlignACE，ANNSpec，Consensus，GLAM，Improbizer，MEME，Mitra，MotifSampler，Oligo/Dyadanalysis，QuickScore，SeSiMCMC，Weeder，YMF）進行了比較[6，7]。它們處理的序列集共有56個，分別來自酵母菌、老鼠、人類和果蠅，序列集的來源有三種：①真實的有確認Motif生物體基因序列；②用HMM產生的隨機序列，其中插入了用權矩陣產生的Motif；③用一個簡單的多項模型產生的序列，其中插入了用權矩陣產生的Motif。這里用tp表示已知的與預測的位置中重疊的長度；fp表示沒有與已知核苷重疊的預測到的核苷數；tn表示不是已知的也沒有預測到的核苷數；fn表示已知的沒有預測到的核苷數。定義靈敏度Sens=tp/(tp+fn)，該值給出了已知的結合位點核苷中預測出的比例；位置預測值ppv=tp/(tp+fp)，該值給出了預測出的位點中已知位點所占的比例；特征值Spec=tn/(tn+fp)，該值給出了未預測出的沒有結合位點的核苷數。另外定義確定的核苷位置與預測的核苷位置之間的相關系數為

nCC=tp·tn-fn·fp(tp+fn)(tn+fp)(tp+fp)(tn+fn)

其值從-1（兩者正好反相關）~+1（兩者完全相關），如果預測的核苷位置在序列中隨機出現，那么nCC的值為0。

具備了這些參數以及相應的實驗數據，可以作出如圖1、圖2所示的圖形。

這里MEME使用了兩個版本，分別獨立運行。可以看到，盡管分別自由地選擇不同的參數和一條序列中的Motif個數（一個或沒有），但MEME和MEME3的結果非常一致。從圖1和圖2可以看出，Weeder的性能比其他的軟件要高，Weeder成功的一個因素是它以警告模式運行，在這種模式下只有最強的Motif才會被報告。當然也有例外，圖2中SeSiMCMC在處理果蠅數據集上效果最好；MEME3和YMF在老鼠數據集上效果更好。待處理的序列集中的Motif的長度可能比較長，而這些算法一般只考慮簡單了比較短的Motif，所以圖中單細胞生物酵母菌的nCC值比其余的三種動物都明顯高，而且最復雜的多細胞生物果蠅的nCC值都很低，甚至低于0。

從上面的分析可以看出，很多不同的搜索算法可以得到非常相似的結果。在Motif的長度不超過20 bp的情況下，即使應用簡單列舉法的軟件YMF都能處理復雜的樣本。另外，一個好的得分函數得到的結果會好得多，如Weeder。因為一般不知道真實的長度是多少，所以該得分函數應該包含比較不同長度Motif的可能性。另外由于無法預先知道序列集的情況，故每一個參數均應該從序列中考慮。

3 結束語

Motif尋找的算法軟件很多，本文在介紹這些算法的基礎上，對尋找Motif的13種軟件進行了性能比較。由于考慮Motif保守核心位置的原因，Weeder在各種軟件中顯示出了比較好的性能。隨著各種技術的發展和人們對分子生物學認識的深入，出現了越來越多的其他方法來識別Motif，如采用比較基因組學來發現在進化過程中保守的結合位點；考慮Motif之間的協同作用而設計的Motif模塊識別方法等。但是，目前還沒有一種軟件能尋找出所有正確的Motif，所以單獨使用一種方法難免出現偏差，將兩種不同的軟件結合起來尋找Motif的性能將會比單獨使用一種軟件高得多。另外一般生物基因組序列均比較長，程序運行需要很長時間，設計并行模式的軟件能大大提高軟件的運行速度。

參考文獻：

[1]Bailey TL， Elkan C. Fitting a Mixture Model by Expectation Maximization to Discover Motifs in Biopolymers[C]. Proc.of the 2nd International Conf. Int. Sys. Mol. Biol.， 1994.28-36.

[2]Pesole G， Prunella N， Liuni S， et al. Wordup: An Efficient Algorithm for Discovering Statistically Significant Patterns in DNA Sequences[J]. Nucleic Acids Res.， 1992，20(11):28712875.

[3]Thijs G， Marchal K， Lescot M， et al. A Gibbs Sampling Method to Detect Overrepresented Motifs in Upstream Regions of Coexpressed Genes[J]. Journal of Computational Biology(Special Issue Recomb 2001)， 2001，9(2):447-464.

[4]Lawrence CE， Altschul SF， Bogouski MS， et al. Detecting Subtle Sequence Signals: A Gibbs Sampling Strategy for Multiple Alignment[J]. Science， 1993，262:208-214.

[5]Pavesi G， Mereghetti P， Mauri G， et al. Weeder Web: Discovery of Transcription Factor Binding Sites in a Set of Sequences from Coregulated Genes[J]. Nucleic Acids Res.， 2004，32:199-203.

[6]Tompa M， Li N， Bailey TL， et al. Assessing Computational Tools for the Discovery of Transcription Factor Binding Sites[J]. Nature Biotechnology， 2005，23:137144.

[7]Maximilian Haubler. Motif Discovery in Promoter Sequences[EB/OL]. http://www.stud.unipotsdam.de/%7Ehau ssler/master/masterthesis.pdf.

作者簡介：

朱驥（1981-），男，江蘇海安人，碩士，主要研究方向為并行計算在生物信息領域中的應用；楊華，女，碩士，主要研究方向為并行計算在生物信息學領域的應用；牛北方，男，碩士，主要研究方向為并行計算在生物信息領域中的應用；郎顯宇，女，博士，主要研究方向為并行計算在生物信息領域中的應用；陸忠華，女，博士后，主要研究方向為生物數學、并行算法；遲學斌，男，中心主任，研究員，博導，博士，主要研究方向為高性能計算方法與軟件。

注:本文中所涉及到的圖表、注解、公式等內容請以PDF格式閱讀原文