基于貪心EM的模體預(yù)測(cè)算法?

張 斐

1 引言

當(dāng)人類(lèi)基因組研究進(jìn)入系統(tǒng)測(cè)序階段時(shí)，大量已測(cè)定的但未知功能或未經(jīng)注釋的DNA序列，急需可靠的自動(dòng)的基因組序列注釋方法和技術(shù)進(jìn)行處理。

人們利用相關(guān)測(cè)序后DNA序列建立了大量轉(zhuǎn)錄因子（TF）數(shù)據(jù)庫(kù)，為轉(zhuǎn)錄調(diào)控的進(jìn)一步研究提供了豐富的資源，它是用計(jì)算機(jī)方法構(gòu)建轉(zhuǎn)錄預(yù)測(cè)模型的基礎(chǔ)。大量的數(shù)據(jù)促進(jìn)了識(shí)別Motif算法和相關(guān)預(yù)測(cè)程序包的發(fā)展。目前，對(duì)識(shí)別Motif的算法研究是一個(gè)熱點(diǎn)，它們大部分都基于一致性序列和矩陣模型。

在應(yīng)用軟件預(yù)測(cè)Motif時(shí)，可以運(yùn)用的方法，有用不同的預(yù)測(cè)軟件比較后取長(zhǎng)補(bǔ)短，參考各軟件的預(yù)測(cè)結(jié)果，從而提出新的預(yù)測(cè)方案，最終得到最佳的預(yù)測(cè)結(jié)果［1］。采用不同的算法來(lái)比較預(yù)測(cè)，通過(guò)預(yù)測(cè)結(jié)果對(duì)比來(lái)選擇合適的算法。還可以在已有算法的基礎(chǔ)上來(lái)進(jìn)行改進(jìn)，在原有序列的基礎(chǔ)上進(jìn)行比較預(yù)測(cè)，結(jié)果更具有說(shuō)服力。本文就是在貪心EM算法的基礎(chǔ)上，通過(guò)篩選優(yōu)化參數(shù)修定序列的預(yù)測(cè)長(zhǎng)度范圍，從而改進(jìn)算法的整體預(yù)測(cè)正確率。

2 貪心EM算法的應(yīng)用基礎(chǔ)

貪心EM算法是實(shí)現(xiàn)機(jī)器學(xué)習(xí)（GMM）綜合學(xué)習(xí)的一種有效方法，它以遞增的方式線性組合多個(gè)“弱”模型（單高斯分量），進(jìn)而得到一個(gè)具有很高似然度的“強(qiáng)”模型（高斯混合模型），使得最后所得到的混合模型可以有效地逼近任意復(fù)雜的概率分布［2］。

2.1 混合Motif模型

首先定義一個(gè)有限字符集 Σ={α1,…,αΩ}，其中 Ω= ||Σ ，對(duì)于任一序列 S=a1a2…aL（ L≥1且ai=Σ），在由序列組成的字符串集Σ中，從位置1到位置L= ||S表示字符串長(zhǎng)度。字符序列aiai+1…ai+w-1表示S中的xi開(kāi)始的長(zhǎng)度為W 的字符串。這里有n=L-W+1表示從序列S中找到的字符串的長(zhǎng)度［3］。

序列S的概率為

由條件概率公式，可得：

定義一個(gè) yi=aiai+1…ai+w+2表示長(zhǎng)度為W-1的字符串。取對(duì)數(shù)后得到似然函數(shù)L：

假設(shè)有N條不同的序列s={s1,…,sN}，其長(zhǎng)度分別為L(zhǎng)1,…,LN，為了識(shí)別長(zhǎng)度為W 的Motif，序列S中的字符串長(zhǎng)度應(yīng)為W［4］。由于每條長(zhǎng)度為W原始序列SS有ms=Ls-W+1種的可能的組合，因而需要建立一個(gè)由n個(gè)字符串n=∑Ns=1ms組

為了解決模型參數(shù)的初始化問(wèn)題，出現(xiàn)了許多算法［5］。其中MEME就是在一次EM迭代后，用動(dòng)態(tài)算法來(lái)估計(jì)基于似然函數(shù)模型的可能起始點(diǎn)。而貪心EM算法是在局部搜索中找到了模型的合適的參數(shù)后再進(jìn)行全局搜索［5］，因此本文引入了EM算法的貪心學(xué)習(xí)過(guò)程。

2.2 貪心混合學(xué)習(xí)

假定一個(gè)θg+1(xi;θg+1)，將其增加到g劃分的混合模型 f(xi;ψg)中。 g+1部分在Motif模型的PWMg+1中含有參數(shù)集θg+1。平衡方程如下：

ψg+1表示一個(gè)新的參數(shù)集，它是由原來(lái)混合模型的參數(shù)集ψg、權(quán)重α和參數(shù)集θg+1組成。對(duì)ψg+1取似然得到：成的數(shù)據(jù)集X={x1,…,xn}。

對(duì)任意一個(gè)字符串xi建立一個(gè)混合Motif模型如下：

重復(fù)執(zhí)行上面過(guò)程得到一個(gè)g值，最后得到接近混合分布的對(duì)數(shù)似然度，上式有兩部分，第一個(gè)分量是原有混合分布 f(xi;ψg)，另一個(gè)分量是φg+1(xi;θg+1)。問(wèn)題就轉(zhuǎn)為運(yùn)用適當(dāng)?shù)姆椒▉?lái)優(yōu)化參數(shù) a 和 θg+1，使得 ψg在最大［6］。

考慮到以上因素，本文將局部和全局搜索結(jié)合起來(lái)，再執(zhí)行EM算法局部搜L(ψg+1)索來(lái)尋找a和θg+1的最大似然，這里通過(guò)學(xué)習(xí)過(guò)程可以出求混合權(quán)a和新加項(xiàng)分量的概率。

引入全局搜索，用泰勒（Taylor）相似a=a0來(lái)代替似然函數(shù)，我們可以通過(guò)從結(jié)果來(lái)查找最優(yōu)θg+1值，在a0=0.5處展開(kāi)L(ψg+1)的二階泰勒公式，并最大化a的二次函數(shù)：

2.3 改進(jìn)初始化模型參數(shù)

用搜索候選字符集的方法直接覆蓋含全部字符串的數(shù)據(jù)集X=，（ τ=1,…,n），這 里xτ=aτ1…aτW

［7］。這樣每一個(gè)字符串 xτ就與 PWM的θτ建立了聯(lián)系。

搜索候選字符串n增加了時(shí)間復(fù)雜度(o(n2))，由于每個(gè)字符串在候選參數(shù)模型中都需要計(jì)算，所以o(n2)是必須的［8］。為了減少時(shí)間復(fù)雜度，我們引入了基于Kd-樹(shù)結(jié)構(gòu)，通過(guò)Kd-樹(shù)對(duì)EM初始化參數(shù)進(jìn)行優(yōu)化，從而建立了基于Kd-樹(shù)的貪心EM預(yù)測(cè)算法（Production Algorithm of Kd-tree Greedy-EM，PKGE算法）。Kd-樹(shù)是通過(guò)定義一個(gè)遞歸的二進(jìn)制的K維數(shù)據(jù)集，它的根節(jié)點(diǎn)包括所有的數(shù)據(jù)，每一層通過(guò)檢測(cè)不同的屬性（關(guān)鍵字）值以決定選擇分枝的方向，從而加快最近臨近查詢(xún)。選擇一棵樹(shù)在水平方向上一個(gè)適當(dāng)?shù)钠矫嬖诖怪钡姆较蚓涂梢员憩F(xiàn)出主要數(shù)據(jù)的變化。從包括全部的字符串集X的根節(jié)點(diǎn)開(kāi)始，每一步都用迭代，在每個(gè)節(jié)點(diǎn)上進(jìn)行分隔。當(dāng)計(jì)算出字符串的k位置上每一個(gè)字符的相對(duì)概率后就包含在這個(gè)節(jié)點(diǎn)上，在位置q上表現(xiàn)的最大熵值［9］：

其中分隔是在搜索的字符位置q上，根據(jù)標(biāo)識(shí)為“奇數(shù)”或“偶數(shù)”。生成了兩個(gè)子集，在位置q上的字符串含有“偶數(shù)”和“奇數(shù)”字符，在圖1中第二個(gè)字符位置出現(xiàn)最大熵。以上的迭代生成了含有多個(gè)節(jié)點(diǎn)的樹(shù)，當(dāng)滿(mǎn)足其包含的字符串?dāng)?shù)量小于T時(shí)，節(jié)點(diǎn)就成為葉子節(jié)點(diǎn)［10］。樹(shù)上每個(gè)節(jié)點(diǎn)都是一個(gè)已經(jīng)初始化的字符串，每一個(gè)這樣的字符串集合通過(guò)它的質(zhì)心被特征化，這些葉子節(jié)點(diǎn)全部是由質(zhì)心和相應(yīng)的PWM組成。

圖1 Kd-樹(shù)結(jié)構(gòu)

3 基于貪心EM算法的改進(jìn)預(yù)測(cè)算法（PKGE算法）

設(shè)定 N條序列S={SS}，長(zhǎng)度LS= ||SS,(s=1,…,N)，假設(shè) Σ={α1,…,αΩ}，長(zhǎng)度從 Ω= ||Σ字符集中取值。算法流程如圖2所示。

初始化：

1）設(shè)定長(zhǎng)度W，在集合S中生成一個(gè)字符串X={xi}(i=1,…,n)，滿(mǎn)足n=∑Ns=1(Ls-W+1)。

2）用Kd-樹(shù)來(lái)處理數(shù)據(jù)集X，通過(guò)質(zhì)心C找到最終的一致性序列，利用全局搜索來(lái)定義C的相關(guān)初始參數(shù)值 θτ，τ=1，…，C［11］。

3）初始化設(shè)定 ?τ=0,?τ=1,…,C ，計(jì)算矩陣數(shù)量 ξτ,i。

4）初始化設(shè)定g=1，設(shè)置背景的參數(shù)等于字符的概率 αl∈ Σ ，如，，這里 fl表示字符αl在序列S中的出現(xiàn)頻率。

參數(shù)迭代部分：

1）EM的收斂條件

2）如果 g≥2，從Motif的候選序列集合 xi（當(dāng)zig＞0.9）中找到與之相鄰的 Ni={xj}(j=i-K,…,i+k)g≥2，并設(shè)置每個(gè)葉子節(jié)點(diǎn)?τ=1，其中τ可包含任意個(gè)xj。

3）增加 g+1部分，搜索 θτ（ τ=1，…，C且?τ=0）并令 θ?g+1等于 θτ，用 ξτ,i來(lái)代替 φ(xi;θτ)，用它的最大似然函數(shù)。再利用公式θτ=[pτl,k]，中的 θ?g+1的值來(lái)計(jì)算權(quán)a?。

4）用 a?和 θ?g+1的初始值執(zhí)行EM，直到執(zhí)行到收斂條件 ψg+1為止［13］。

5）如果 L(ψg+1)＞L(ψg)，則利用 g+1劃分的混合模型并重新轉(zhuǎn)向1步，否則終止。直到求出Motif的最大數(shù)量，或者

算法描述了在Motif預(yù)測(cè)中的一個(gè)混合Motif模型。由于EM算法不能減小似然，所以只有當(dāng)滿(mǎn)足 L(ψg+1)＞L(ψg)時(shí)才利用局部EM迭代來(lái)解決問(wèn)題，從而保證了在數(shù)據(jù)集中似然函數(shù)的單調(diào)增加。算法終止的條件包含g的取值和向量集?的最大值。在這里 ?τ=1，?τ=1,…,C［14］，候選序列部分的參數(shù)空間被徹底搜索，因此其它Motif在序列大量存在的可能性是很低的。

圖2 PKGE算法流程圖

4 算法測(cè)試結(jié)果

從PRINTS數(shù)據(jù)庫(kù)下載部分蛋白質(zhì)家族序列，隨機(jī)產(chǎn)生8條長(zhǎng)度在150bp～300bp之間［8，12］的蛋白質(zhì)序列，通過(guò)多次比對(duì)設(shè)定Kd-樹(shù)劃分值為6，葉子節(jié)點(diǎn)字符串的長(zhǎng)度設(shè)定為5，Motifs長(zhǎng)度為10，由于選定序列源于生物數(shù)據(jù)庫(kù)，雖經(jīng)過(guò)人工拼接，但仍有部分片段保留，因而程序執(zhí)行后仍能找到表1所示數(shù)量的Motifs。

表1 預(yù)測(cè)的PROSITE蛋白質(zhì)家族［15］

PKGE算法能顯著降低Motif的重疊率，而且還能保證Motif候選序列參照一致的規(guī)則進(jìn)行處理，此外通過(guò)迭代還減少了全局搜索的時(shí)間復(fù)雜度，與MEME算法相比，相鄰的Motif字符串從數(shù)據(jù)集中剔除了，重疊的字符串被排除在候選Motif集合之外，這就相對(duì)MEME算法有了改進(jìn)。PKGE算法主要是針對(duì)MEME的局限性來(lái)設(shè)計(jì)的，例如每一次輸入一個(gè)新的Motif序列時(shí)需要假定這個(gè)序列是真實(shí)存在的，從而限制了MEME模型中兩部分的聯(lián)系，PKGE算法通過(guò)增加混合密度估計(jì)改進(jìn)了這個(gè)問(wèn)題。

表2 算法性能比較

5 結(jié)語(yǔ)

算法實(shí)現(xiàn)了預(yù)測(cè)性能的一定改進(jìn)，但其特異性的絕對(duì)水平依然很低。對(duì)眾多預(yù)測(cè)算法的評(píng)估研究表明，它們通常僅在特定的小規(guī)模數(shù)據(jù)集上表現(xiàn)出良好的性能，而對(duì)更為廣泛和通用的數(shù)據(jù)集，算法產(chǎn)生了大量假陽(yáng)性結(jié)果，這應(yīng)該是導(dǎo)致此現(xiàn)象的直接原因，此外，限于當(dāng)前的認(rèn)識(shí)水平，假陽(yáng)性結(jié)果中有些可能確實(shí)是尚未發(fā)現(xiàn)的結(jié)合位點(diǎn)，其功能需要后續(xù)的實(shí)驗(yàn)來(lái)驗(yàn)證。

［1］張斐，譚軍，謝競(jìng)博.基于不同算法的motif預(yù)測(cè)比較分析與優(yōu)化［J］.計(jì)算機(jī)工程，2009，35（22）：94-96.

ZHANG Fei，TAN Jun，XIE Jingbo.Comparison，Analysis and Optimization of Motif Finding Based on Different Algorithms［J］.Computer Engineering，2009，35（22）：94-96.

［2］王維彬，鐘潤(rùn)添.一種基于貪心EM算法學(xué)習(xí)GMM的聚類(lèi)算法［J］.計(jì)算機(jī)仿真，2007，24（2）：65-68.

WANG Weibin，ZHONG Runtian.A Clustering Algorithm Based on Greedy EM Algorithm Learning GMM［J］.Computer Simulation，2007，24（2）：65-68.

［3］J J Verbeek，N Vlassis，B Krose.Efficient greedy learning of Gaussian mixture［J］.Neural Computation，2003，2（15）：469-485.

［4］S K Palaniswamy，S James，H Sun，et al.AGRIS and AtRegNet.A platform to link cis-regulatory elements and transcription factors into regulatory networks［J］.Plant Physiol，2006（140）：818-829.

［5］MTompa，NLi，TLBailey，et al.assessing computational tools for the discovery of transcription factor binding sites［J］.Nat Biotechnol，2005，23：137-144.

［6］T LBailey.Discovering Motifs in DNA and protein sequences：theappoximate common substring problem［D］.San Diego：University of California，1995，12-13.

［7］T LBailey，CElkan.Unsupervised learning of multipleMotifs in biopolymers using expectation maximization［J］.Machine Learning，1995，21：51-83.

［8］L BTimothy，WNadya，M Chris，et al.MEME：discovering and analyzing DNA and protein sequence Motifs［J］.Nucleic Acids Research，2006，34（Web Server issue）：369-373.

［9］A P Dempster，N M Laird，D B Rubin.Maximum likelihood from incomplete data via the EM algorithm［C］//J.Roy.Statist.Soc，B，1977，39：1-38.

［10］吳昕，羅靜初，李伍舉.基因調(diào)控元件的計(jì)算機(jī)識(shí)別和基因調(diào)控網(wǎng)絡(luò)構(gòu)建［J/OL］.http：//www.cbi.pku.edu.cn/chinese/documents/papers/WuX.pdf.

WU Xin，LUO Jingchu，LI Wuju.Computer recognition of regulatory elements and gene regulatory network［J/OL］.http：//www.cbi.pku.edu.cn/chinese/documents/papers/WuX.pdf.

［11］韓華，劉婉璐，吳翎燕.基于模體的復(fù)雜網(wǎng)絡(luò)測(cè)度量研究［J］.物理學(xué)報(bào)，2013（16）：1-9.

HAN Hua，LIU Wanlu，WU Lingyan.The measurement of complex network based on motif［J］.Acta Physica Sinica，2013（16）：1-9.

［12］TObayashi，K Kinoshita，K Nakai，et al.ATTED-II：a database of co-expressed genes and cis elements for identifying co-regulated gene groups in Arabidopsis［J］.Nucleic Acids Res，January 12，2007，35（suppl_1）：863-869.

［13］覃桂敏，高琳，呼加璐.生物網(wǎng)絡(luò)模體發(fā)現(xiàn)算法研究綜述［J］.電子學(xué)報(bào)，2009（10）：2258-2265.

QIN Guimin，GAO Lin，HU Jialu.A Review on Algorithms for Network Motif Discovery in Biological Networks［J］.Acta Electronica Sinica，2009（10） ：2258-2265.

［14］霍紅衛(wèi)，郭丹丹，于強(qiáng)，等.（l，d）-模體識(shí)別問(wèn)題的遺傳優(yōu)化算法［J］.計(jì)算機(jī)學(xué)報(bào)，2012，35（7）：1429-1439.

HUO Hongwei，GUO Dandan，YU Qiang，et al.Genetic Optimization for（l，d）-Motif Discovery［J］.Chinese Journal of Computers，2012，35（7）：1429-1439.

［15］宋濤.基于譜隱馬爾可夫模型的蛋白質(zhì)序列模體識(shí)別方法研究［D］.大連：大連理工大學(xué)，2015：60-61.

SONG Tao.Research on Methods for Protein Sequence Motif Discovery Based on Profile Hidden Markov Model［D］.Dalian：Dalian University of Technology，2015：60-61.