模體發現問題中OOPS模型的EM算法

黃影

摘 要 生物序列模體發現是從生物序列集合找出具有特定功能子序列的問題，它是生物信息學中最重要問題之一。OOPS模型是一種簡單常用的模體發現的序列模型，它表示每條序列中有一個模體出現。本文提出了一種適用于OOPS序列模型的EM算法，用于查找每條序列中所出現的模體。實驗結果表示，所提出的算法的時間性能和識別準確率均好于經典的Project算法。

關鍵詞 模體發現 生物序列 EM算法

中圖分類號：TP181 文獻標識碼：A DOI：10.16400/j.cnki.kjdkx.2015.08.010

EM Algorithm of OOPS Model in Motif Discovery Problem

HUANG Ying

（School of Information Engineering， Xi'an University， Xi'an， Shaanxi 710065）

Abstract Biological sequence motif discovery from biological sequence set to find out having a problem specific functional sub-sequence， it is one of the most important issues in bioinformatics. OOPS model series model is a simple common motif finding， which represents each sequence has a motif appears. This paper presents a suitable OOPS series model EM algorithm， for finding that appears in each sequence motif. The results indicated that the time performance of the algorithm and the recognition accuracy rate of the proposed algorithm is better than the classical Project.

Key words motif discovery; biological sequence; EM algorithm

1背景

EM（expectation-maximization）算法是Dempster， Laird和Rubin（DLR）三個人在1977年正式提出的。主要是用于在不完全數據的情況下計算最大似然估計。在EM算法正式提出以來，人們對EM算法的性質有更加深入的研究;并且在此基礎上，提出了很多改進的算法，在數理統計，數據挖掘，機器學習以及模式識別等領域有廣泛的應用。

生物序列模體（motif）發現問題是從一個從相關生物序列集合找出具有某種特定功能的子序列的問題，它是生物信息學中最重要問題之一，也是最基本問題之一。針對motif在每條序列中出現或不出現或者出現多次，有三種對應的模型，即：若每條序列中有且僅有一個motif出現，則被稱為OOPS模型，這是最簡單的模型;若每條序列中至多有一個motif出現，則被稱為ZOOPS模型;若每條序列允許多個motif出現或者沒有motif出現，則被稱為TCM模型，這是最復雜的模型。

本文實現基于最簡單的OOPS模型實現EM算法，來查找每條序列中所出現的motif。

2模體發現的EM算法

在最大化似然的問題中，給定被稱為數據參數似然值的函數（ ; ），或者就稱為似然函數：

（∣ ） = （∣ ） = （ ; ）

我們的目標就是找到使最大的 值，即

* = ? （ ; ）

EM算法的核心思想就是根據已有的數據來遞歸估計似然函數。

所有的EM算法都由兩個主要的步驟組成：

E-step

M-step

在E-step中，對于未知的潛在變量，使用當前參數的估計值和當前的觀察。在M-step中，重新計算出來一個新的參數的估計值。在每次迭代中，似然值都會增加，因此該過程總會達到一個漸近最大值;但EM可能陷入一個局部最優。通過交替使用這兩個步驟，EM算法逐步改進模型的參數，使參數和訓練樣本的似然概率逐漸增大，最后終止于一個極大點。

直觀地理解EM算法，它也可被看作為一個逐次逼近算法：事先并不知道模型的參數，可以隨機的選擇一套參數或者事先粗略地給定某個初始參數，確定出對應于這組參數的最可能的狀態，計算每個訓練樣本的可能結果的概率，在當前的狀態下再由樣本對參數修正，重新估計參數，并在新的參數下重新確定模型的狀態，這樣，通過多次的迭代，循環直至某個收斂條件滿足為止，就可以使得模型的參數逐漸逼近真實參數。

EM算法的主要目的是提供一個簡單的迭代算法計算后驗密度函數，它的最大優點是簡單和穩定，但容易陷入局部最優。

如果給定一組已比對好的序列集，可以直接構造一個表征motif的profile矩陣，即PWM;圖1是一個例子。

但對于未比對的序列，我們怎么構造profile矩陣？因為在通常情況下，我們并不知道motif是什么，這時我們就可以使用EM算法。

使用概率矩陣 = 表示motif概率，這里表示列的字符所出現的概率。用概率 = 來表示背景（即，序列中motif之外的部分）概率，用表示背景中的字符的出現概率。并定義矩陣，其元素表示motif出現在序列的位置的概率。例如，在圖2中，給定長度為6的3條DNA序列，這里motif的長度W=3，其對應的、和Z：

圖1 給定一組已比對好的序列集，直接構造表征motif的profile矩陣

圖2 給定長度為6的3條DNA序列，motif的長度W=3，對應的、和

整個EM算法的框架如下：

EM algorithm

Given： length parameter W， training set of sequences

set initial values for P

do ? ?re-estimate Z from P ? ? ? ? ? （E-step）

re-estimate P from Z ? ? ? ? ? （M-step）

until change in P <

Return： P， Z

在E-step中， 由計算的估計即motif出現的最可能位置;在M-step中，由重新計算的估計即一個新的motif的PWM。

假定給定條長度為的序列： = （， ，…，）， = 1， …，，我們找出長度為的motif。則該序列在位置j出現motif的概率為：

在E-step：由計算

在E-step的第次迭代中，由下式計算獲得：

= （ = 1∣ ， ， ）

=

假定motif等可能的出現在序列的任意位置，則上式可化簡為：

=

在M-step：由計算

對于motif，對任意的，有

=

這里 = ?。

對于背景，上式的可用 = 計算得到。這里表示序列中堿基的總個數，是為避免計算中出現零概率而引入的擾動。

3實驗結果

我們在隨機生成的模擬數據上對我們的程序進行了測試。隨機生成t條滿足均勻分布的DNA序列，在每條序列中植入（l， ）motif，即motif長度為l，發生個變異。對所實現算法的正確性，我們使用性能系數進行評估：

performace coefficient =

這里表示所植入的motif在序列中的位點集， 表示由算法所找到的motif在序列中的位點集。

我們將我們的實現同project方法進行了比較，結果如表1所示：

表1 EM和project的pc值

表1中每項表示在十個樣本上的平均的性能系數，其中每組包含20條數據，每條序列長為600。每輪測試進行10次迭代。對于（12， 3）、（15， 4）、 （18， 5）、（21， 6）、（24， 7），我們的實現的性能系數分別為0.6474、0.8740、0.9810、0.9723、1，同相應的project的方法相差不大，而它的性能系數分別為0.6576、0.8740、0.9905、0.9723、1。

同時，我們改變motif的長度，對兩個算法的執行時間進行比較（單位：s），如表2所示：

表2 EM和project的運行時間

從表2中我們可以看出，我們的方法隨著所找motif長度的增加，時間增長比較慢;而相應的project方法所用的時間增長很快。

參考文獻

[1] Arthur Dempster， Nan Laird， and Donald Rubin. Maximum likelihood from incomplete data via the EM algorithm.Journal of the Royal Statistical Society，Series B，1977.39（1）：1-38.