一種新的序列集保守信號(hào)識(shí)別算法

劉暢

摘 要：模體發(fā)現(xiàn)是計(jì)算機(jī)科學(xué)中的一個(gè)較為重要且具有一定挑戰(zhàn)的問(wèn)題，主要用于定位DNA序列集中的保守信號(hào)。首先，分析了已有的基于圖聚類的模體發(fā)現(xiàn)算法MCL-WMR，討論了它存在的兩個(gè)缺陷。其次，針對(duì)這兩個(gè)缺陷提出了MCL-WMR的改進(jìn)算法iMCL-WMR。實(shí)驗(yàn)結(jié)果表明，所提DNA模體發(fā)現(xiàn)算法的時(shí)間性能好于所比較的算法MCL-WMR和qPMS9，能夠在1個(gè)小時(shí)以內(nèi)處理數(shù)百條輸入序列，而且能夠應(yīng)對(duì)某些輸入序列不含模體實(shí)例的測(cè)試數(shù)據(jù)。

關(guān)鍵詞：保守信號(hào)；模體發(fā)現(xiàn)；圖聚類

中圖分類號(hào)：TP301.6 文獻(xiàn)標(biāo)識(shí)碼：A

Abstract：Motif discovery is an important and challenging issue in computer science，mainly used to locate conserved signals in a set of DNA sequences.At first，a graph clustering based motif discovery algorithm MCL-WMR is analyzed and its two drawbacks are discussed.Then，in order to overcome the two drawbacks，an improved algorithm of MCL-WMR，named as iMCL-WMR，is proposed.Experimental results show that，with a better time performance than the compared algorithms of MCL-WMR and qPMS9，the proposed algorithm can process hundreds of sequences within one hour and deal with any case that some input sequences do not contain motif instances.

Keywords：conserved signals；motif discovery；graph clustering

1 引言（Introduction）

序列集中保守的信號(hào)稱為模體。模體發(fā)現(xiàn)就是在一組給定的序列集合中找出未知的模體，是很多重要的數(shù)據(jù)挖掘任務(wù)的第一步[1-3]。本文主要關(guān)注DNA序列集合中的模體發(fā)現(xiàn)問(wèn)題[4]，可用于定位許多具有重要功能的基因區(qū)域。問(wèn)題形式化描述如下[5]：

給定字符表{A，C，G，T}上的n條長(zhǎng)為m的序列集合S={sl，s2，…，sn}，目標(biāo)是找出一個(gè)長(zhǎng)為l（0

此問(wèn)題是一個(gè)NP難問(wèn)題[6]，當(dāng)前已涌現(xiàn)了許多求解算法[7]。一類算法通過(guò)窮盡所有可能的l-mer來(lái)找出潛在的模體[8，9]。另一類算法將模體視為一個(gè)位置頻率矩陣（Position Frequency Matrix，PFM）[10]，一般對(duì)初始的PFM進(jìn)行迭代優(yōu)化來(lái)得到模體。還有一類算法將模體發(fā)現(xiàn)問(wèn)題抽象為圖聚類問(wèn)題，比如MCL-WMR[11]。

由于此問(wèn)題的NP難特性，目前它仍然是一個(gè)開(kāi)放性的難題；此外，當(dāng)前待處理的模體發(fā)現(xiàn)的數(shù)據(jù)規(guī)模在變大[12]，許多現(xiàn)有算法在時(shí)間性能上達(dá)不到實(shí)用要求。本文通過(guò)分析和改進(jìn)MCL-WMR，提出了一種新的模體發(fā)現(xiàn)算法，它的時(shí)間性能更為高效且能夠處理某些輸入序列不含模體實(shí)例的情況。

2 相關(guān)工作（Related works）

數(shù)據(jù)聚類是將一個(gè)給定的數(shù)據(jù)集中的數(shù)據(jù)項(xiàng)進(jìn)行分組，使得同一組中的數(shù)據(jù)項(xiàng)高度相似并顯著區(qū)分于其他組中的數(shù)據(jù)項(xiàng)。對(duì)于模體發(fā)現(xiàn)，由于同一轉(zhuǎn)錄因子的結(jié)合位點(diǎn)在某種程度上相互相似，采用合適的相似性度量的聚類算法可以將輸入的DNA序列集合中的結(jié)合位點(diǎn)聚集在一起。然而，由于模體在序列中往往是高度退化的，定義合適的相似性度量和聚類模型對(duì)于在沒(méi)有任何先驗(yàn)的情況下找到相關(guān)的結(jié)合位點(diǎn)是一個(gè)困難的任務(wù)。

模體的表示模型在設(shè)計(jì)模體發(fā)現(xiàn)算法中是至關(guān)重要的。一個(gè)最廣為使用的模體表示模型是位置頻率矩陣PFM。一個(gè)長(zhǎng)度為l的模體的PFM是一個(gè)矩陣Mpfm=[abi]4×l，其中b={A，C，G，T}，i=1，…，l，并且對(duì)于任意的i有∑abi=1。PFM表示一組對(duì)齊的等長(zhǎng)位點(diǎn)中一個(gè)核苷酸在某個(gè)位置上出現(xiàn)的概率，它也表示了蛋白質(zhì)與DNA進(jìn)行交互的強(qiáng)度。

MCL是一個(gè)知名的圖聚類算法[11]。它的基本原理是模擬隨機(jī)游走?？紤]一個(gè)圖中的聚類，聚類內(nèi)部將含有許多邊，而聚類之間將含有較少的邊；這意味著，從圖中的一個(gè)頂點(diǎn)出發(fā)隨機(jī)游走到有邊相連的另一個(gè)頂點(diǎn)時(shí)，更有可能停留在一個(gè)聚類內(nèi)部，而跨越到另一個(gè)聚類的概率會(huì)比較低。通過(guò)在圖中進(jìn)行隨機(jī)游走，可以發(fā)現(xiàn)流（flow）聚焦的地方，即聚類之所在。圖中的隨機(jī)游走使用馬爾可夫鏈進(jìn)行計(jì)算，即通過(guò)計(jì)算圖的概率矩陣的連續(xù)的乘方來(lái)模擬流。MCL每次迭代中主要包含擴(kuò)展和膨脹兩個(gè)步驟，其中膨脹參數(shù)r決定了聚類的粒度。

MCL-WMR[11]是一個(gè)基于圖聚類的模體發(fā)現(xiàn)算法，它的偽代碼如下所示。它將輸入的DNA序列集合轉(zhuǎn)為一個(gè)帶權(quán)圖，然后采用MCL對(duì)圖進(jìn)行聚類從而得出模體。序列中每個(gè)l-mer被表示為一個(gè)頂點(diǎn)，在圖的構(gòu)建中保證與模體實(shí)例相對(duì)應(yīng)的兩個(gè)頂點(diǎn)間要有邊相連接，這樣模體發(fā)現(xiàn)轉(zhuǎn)化為了在圖中尋找大小為n的團(tuán)。

Algorithm MCL-WMR

Input：l，d，S={s1，s2，…，sn}

Output：a motif

1：build a graph G according to the given data set S

2：select a reference sequence sr from S randomly

3：for each l-mer x in sr do

4：generate a subgraph G' induced by x

5：cluster G' using MCL

6：for each cluster do

7：if the weight of the current cluster is above threshold then

8：verify that if the cluster contains a motif

9：if find a motif，output it

建圖的具體過(guò)程為：（i）頂點(diǎn)集中的頂點(diǎn)vi，j表示第i條序列位置j處的l-mer，對(duì)于每個(gè)i，j=1，2，…，m-l+1，所以共有n（m-l+1）個(gè)頂點(diǎn)；（ii）對(duì)于每對(duì)頂點(diǎn)vi，j和vi'，j'（i≠i'），當(dāng)表示的兩個(gè)l-mer的海明距離小于等于2d時(shí)，vi，j和vi'，j'間有一條邊相連；（iii）當(dāng)有邊相連的兩個(gè)頂點(diǎn)的海明距離k滿足d

在MCL-WMR算法中，為了減少M(fèi)CL所處理的頂點(diǎn)數(shù)量，從輸入序列集合S中任取一條參考序列sr，然后由sr中的每個(gè)l-mer x誘導(dǎo)一個(gè)子圖G'，使得G'是由{S-sr}中與x滿足海明距離小于等于2d的所有l(wèi)-mer組成。然后分別對(duì)G'用MCL進(jìn)行聚類；由于任意兩個(gè)模體實(shí)例滿足海明距離小于等于2d，所以必然存在一個(gè)子圖G'，其中包含由n個(gè)模體實(shí)例組成的團(tuán)。

3 方法（Method）

3.1 算法思想

MCL-WMR存在如下的兩個(gè)缺陷。

首先，MCL-WMR假設(shè)每條輸入序列至少含有一個(gè)模體實(shí)例，因此不能應(yīng)對(duì)有部分序列不含有模體實(shí)例的情況。在真實(shí)的情況中，生物學(xué)家產(chǎn)生的數(shù)據(jù)或經(jīng)過(guò)初步處理的數(shù)據(jù)中并不能保證每條序列都含有模體的出現(xiàn)。此外，當(dāng)模體以高度退化的形式出現(xiàn)在某些序列中時(shí)，它并不能顯著地區(qū)分于背景中的子串。根據(jù)模體在序列中出現(xiàn)的情況，可以分為OOPS、ZOOPS和TCM三種模型，分別表示模體在每條序列中出現(xiàn)一次、模體在每條序列中出現(xiàn)零次或一次、模體在每條序列中出現(xiàn)零次或多次。顯然，MCL-WMR只能應(yīng)對(duì)OOPS模型的情況，而不能應(yīng)對(duì)更符合實(shí)際情況的ZOOPS和TCM模型的情況。

其次，MCL-WMR不能應(yīng)對(duì)規(guī)模較大的輸入數(shù)據(jù)。MCL-WMR在處理經(jīng)典的20條長(zhǎng)為600的序列集時(shí)，有較好的時(shí)間性能。近年來(lái)，隨著生物實(shí)驗(yàn)和測(cè)序技術(shù)的發(fā)展，染色質(zhì)免疫共深沉與高通量測(cè)序相結(jié)合的ChIP-seq技術(shù)[12]可以使人們?cè)诨蚪M水平上進(jìn)行模體發(fā)現(xiàn)，但處理的數(shù)據(jù)量包含數(shù)百條甚至更多條序列。當(dāng)數(shù)據(jù)規(guī)模變大時(shí)，即使MCL-WMR對(duì)多個(gè)子圖分別進(jìn)行處理，但每個(gè)子圖要處理的數(shù)據(jù)量依然很大，使得MCL聚類的時(shí)間效率變得很低。

為了解決這兩個(gè)問(wèn)題，本文對(duì)MCL-WMR進(jìn)行改進(jìn)，提出一個(gè)稱為iMCL-WMR的算法?；舅枷肴缦拢涸u(píng)估輸入序列中每個(gè)l-mer的信號(hào)強(qiáng)度，將信號(hào)較強(qiáng)的l-mer視為潛在的模體實(shí)例；由每個(gè)潛在的模體實(shí)例誘導(dǎo)一個(gè)子圖G'，使得G'中有邊相連的兩個(gè)l-mer間的海明距離小于等于d；對(duì)G'進(jìn)行聚類，對(duì)得到的各個(gè)聚類進(jìn)行求精，最后選擇得分最大的進(jìn)行輸出。

iMCL-WMR不受限于每條序列必須包含一個(gè)模體實(shí)例的OOPS模型，它能夠適應(yīng)ZOOPS和TCM模型的情況。此外，當(dāng)處理的數(shù)據(jù)規(guī)模較大時(shí)，存在一些較為相似的模體實(shí)例，即它們之間的海明距離小于等于d；因此，在由一個(gè)潛在的模體實(shí)例誘導(dǎo)一個(gè)子圖時(shí)，iMCL-WMR限定有邊相連的兩個(gè)l-mer間的海明距離小于等于d，這樣iMCL-WMR中的子圖的數(shù)據(jù)規(guī)模將小于MCL-WMR中的子圖的數(shù)據(jù)規(guī)模，從而提高M(jìn)CL聚類的時(shí)間效率。

3.2 劃分子序列集合

首先，評(píng)估各個(gè)l-mer的信號(hào)強(qiáng)度。令si，j表示第i條序列si位置j處的l-mer。一個(gè)任意的l-mer x與給定的序列si的距離dis（x，si），定義為x與si中的l-mer的可能的最小海明距離，如式（1）所示；其中，dH（x，si，j）表示兩個(gè)l-mer的海明距離，即兩個(gè)l-mer對(duì)齊后字符不相等的位置的個(gè)數(shù)。

（1）

令SD（x）表示輸入序列中的一個(gè)l-mer x的信號(hào)強(qiáng)度，定義為x與各個(gè)輸入序列si中最為相似的l-mer的總的相同的位置個(gè)數(shù)，如式（2）所示。

（2）

一個(gè)l-mer x的信號(hào)強(qiáng)度越大，它越有可能是一個(gè)潛在的模體實(shí)例。這可以從模體發(fā)現(xiàn)的問(wèn)題屬性來(lái)解釋：模體的各個(gè)實(shí)例是由同一模體變異而產(chǎn)生的，所以各個(gè)實(shí)例間是相互相似的；這樣，如果一個(gè)l-mer x是模體實(shí)例，那么部分序列中的模體實(shí)例將與其非常相似，使得x的信號(hào)強(qiáng)度顯著地大于背景中的或隨機(jī)的l-mer的信號(hào)強(qiáng)度。

計(jì)算出輸入序列中所有的l-mer si，j的信號(hào)強(qiáng)度SD（si，j）后，取其均值和標(biāo)準(zhǔn)差，分別記為μ和σ。令L表示輸入序列中所有的l-mer si，j的集合，即L={si，j：1≤i≤n，1≤j≤m}。將L中的l-mer近似成正態(tài)分布，如圖1所示，橫軸表示信號(hào)強(qiáng)度，縱軸表示對(duì)應(yīng)信號(hào)強(qiáng)度的l-mer的數(shù)量?；诖?，將L劃分成L1、L2和L3三部分。

（3）

（4）

（5）

根據(jù)序列中l(wèi)-mer的信號(hào)強(qiáng)度的定義，對(duì)所有子序列的集合的劃分解釋如下。L1是背景序列區(qū)，即它對(duì)應(yīng)很小的信號(hào)強(qiáng)度，包含模體實(shí)例的概率很低。L3是潛在的模體實(shí)例區(qū)，即它對(duì)應(yīng)很大的信號(hào)強(qiáng)度，包含模體實(shí)例的概率很高。L2是模體實(shí)例和背景序列的混合區(qū)，它既包含模體實(shí)例也包含背景序列，且難以進(jìn)行區(qū)分?；谶@個(gè)劃分，iMCL-WMR在構(gòu)建圖時(shí)：一方面，僅使用L2和L3中的元素構(gòu)建圖G，這樣可以減小數(shù)據(jù)規(guī)模；另一方面，僅使用L3中的元素誘導(dǎo)子圖G'，這樣可以控制子圖的數(shù)量。

3.3 聚類求精

為了保證良好的時(shí)間性能，iMCL-WMR執(zhí)行聚類的子圖中并不一定包含全部的模體實(shí)例，這樣得到的聚類C是一個(gè)包含部分模體實(shí)例的集合。因此，需要對(duì)C進(jìn)行求精，得到全部的模體實(shí)例及模體。求精算法如下。

Algorithm RefineCluster

Input：l，d，S={s1，s2，…，sn}，C

Output：a motif

1：generate a consensus c according to aligned elements in C

2：score←∞，curScore←0

3：while curScore>score do

4：score←curScore，C←Ф

5：for i←1 to n do

6.put an l-mer x in si with dH（c，x）=dis（c，si） to C

7：generate a new consensus c according to aligned elements in C

8：curScore ← IC（c）

9：output c and its score curScore

求精算法迭代地對(duì)當(dāng)前的模體進(jìn)行更新，直到得到信息量最高的模體。模體c的信息量記為IC（c），計(jì)算方法如下所示，其中fij表示C中l(wèi)-mer對(duì)齊后第j列的字符i的出現(xiàn)頻率，bi表示字符i的背景頻率。

（6）

3.4 整體算法

iMCL-WMR算法的偽代碼描述如下。

Algorithm iMCL-WMR

Input：l，d，S={s1，s2，…，sn}

Output：a motif

1：compute the signal degree for each l-mer in S

2：divide all l-mers in S to three sets L1，L2 and L3

3：build a graph G using the l-mers in L2 and L3

4：for each l-mer x in L3 do

5：generate a subgraph G' induced by x

6：cluster G' using MCL

7：for each obtained cluster C do

8：refine C using algorithm RefineCluster

9：obtain a motif c and its score curScore

10：rank the obtained motifs according to their scores

11：output the motif with maximum score

算法第1行偽代碼計(jì)算輸入序列中各個(gè)l-mer的信號(hào)強(qiáng)度。依據(jù)信號(hào)強(qiáng)度，第2行偽代碼將輸入序列中的l-mer劃分成L1、L2和L3三個(gè)集合。第3行偽代碼使用L2和L3中的l-mer構(gòu)建一個(gè)圖G，注意iMCL-WMR要求兩個(gè)l-mer間的海明距離k小于等于d時(shí)對(duì)應(yīng)的兩個(gè)頂點(diǎn)才有邊相連接；邊的權(quán)值設(shè)置為（d-k）（l-k），這比MCL-WMR的權(quán)值設(shè)置方法更為精細(xì)，使得海明距離越小時(shí)越以更大的概率視為一對(duì)模體實(shí)例。第4—9行偽代碼對(duì)L3中的各個(gè)l-mer誘導(dǎo)的子圖依次進(jìn)行聚類求精，得到多個(gè)模體。第10—11行偽代碼對(duì)得到的多個(gè)模體依照其得分進(jìn)行排序，最終輸出得分最高的模體。

4 實(shí)驗(yàn)結(jié)果（Experimental results）

首先，使用經(jīng)典的模擬數(shù)據(jù)集（隨機(jī)生成20條長(zhǎng)為600的DNA序列，然后在每條序列中植入一個(gè)（l，d）模體，使得模體與模體實(shí)例間最大差異d個(gè)位置），比較所提算法iMCL-WMR與現(xiàn)有算法的運(yùn)行時(shí)間和識(shí)別準(zhǔn)確率。除了與MCL-WMR進(jìn)行比較之外，又選取了qPMS9作為比較算法，它是最近提出的算法并且是時(shí)間性能最好的精確算法。所有算法都是用C++實(shí)現(xiàn)的，并且運(yùn)行于相同的平臺(tái)（3.0GHz的處理器和4GB的內(nèi)存）。選用的（l，d）為挑戰(zhàn)問(wèn)題實(shí)例，即（9，2）、（11，3）、（13，4）、（15，5）、（17，6）、（19，7）、（21，8）、（23，9）、（25，10）；對(duì)于每個(gè)（l，d）問(wèn)題實(shí)例，隨機(jī)生成了六組數(shù)據(jù)，實(shí)驗(yàn)結(jié)果取均值。識(shí)別準(zhǔn)確率用核苷酸水平的性能系數(shù)來(lái)衡量，它等于預(yù)測(cè)出的模體中的核苷酸與真實(shí)的模體中的核苷酸的交集的大小除以其交并集的大小。

表1給出了這些（l，d）問(wèn)題實(shí)例上的實(shí)驗(yàn)結(jié)果，其中s表示秒、m表示分鐘、h表示小時(shí)。整體來(lái)看，三個(gè)算法的識(shí)別準(zhǔn)確率是相當(dāng)?shù)?；qPMS9的識(shí)別準(zhǔn)確率略好于iMCL-WMR和MCL-WMR的識(shí)別準(zhǔn)確率，因?yàn)樗且粋€(gè)精確算法，窮盡地遍歷了所有的候選模體。對(duì)于時(shí)間性能，iMCL-WMR在整體上好于MCL-WMR和qPMS9；qPMS9在求解小的（l，d）問(wèn)題實(shí)例時(shí)具有優(yōu)勢(shì)，這是因?yàn)榍蠼庑〉模╨，d）問(wèn)題實(shí)例需要窮盡遍歷的候選模體的數(shù)量較少，但在求解大的（l，d）問(wèn)題實(shí)例時(shí)它的時(shí)間性能很差；iMCL-WMR的時(shí)間性能好于MCL-WMR的原因是對(duì)MCL聚類的子圖進(jìn)行了降維處理。endprint

其次，固定（l，d）為（15，5）問(wèn)題實(shí)例以及每條序列的長(zhǎng)度為600，變化序列數(shù)量n，同時(shí)只在80%的序列中植入模體實(shí)例。由于MCL-WMR不支持非OOPS模型的情況，這里只與qPMS9進(jìn)行了比較。實(shí)驗(yàn)結(jié)果如表2所示?？梢园l(fā)現(xiàn)，在求解非OOPS模型的輸入數(shù)據(jù)時(shí)，所提算法iMCL-WMR的時(shí)間性能顯著地好于qPMS9的時(shí)間性能。這是因?yàn)?，iMCL-WMR依據(jù)信號(hào)強(qiáng)度劃分輸入序列中的l-mer，再誘導(dǎo)子圖進(jìn)行聚類的方法可以很自然地處理非OOPS模型的情況；qPMS9是將輸入數(shù)據(jù)轉(zhuǎn)化成多個(gè)OOPS模型的數(shù)據(jù)再進(jìn)行窮舉，所以當(dāng)數(shù)據(jù)規(guī)模較大時(shí)，時(shí)間性能很差。

5 結(jié)論（Conclusion）

本文提出了一種新的模體發(fā)現(xiàn)算法iMCL-WMR。實(shí)驗(yàn)結(jié)果表明，一方面，iMCL-WMR的時(shí)間性能好于所比較的算法MCL-WMR和qPMS9，能夠在1個(gè)小時(shí)以內(nèi)處理數(shù)百條輸入序列；另一方面，iMCL-WMR能夠應(yīng)對(duì)某些輸入序列不含模體實(shí)例的測(cè)試數(shù)據(jù)?？傊?，iMCL-WMR在處理較大規(guī)模和更符合真實(shí)場(chǎng)景的DNA序列數(shù)據(jù)集時(shí)具有較好的實(shí)用性。

參考文獻(xiàn)（References）

[1] Wong KC，et al.A comparison study for DNA motif modeling on protein binding microarray[J].IEEE/ACM Transactions on Computational Biology and Bioinformatics，2016，13

（2）：261-271.

[2] Kranjc J，et al.ClowdFlows：Online workflows for distributed big data mining[J].Future Generation Computer Systems，2017，68：38-58.

[3] Liu B，et al.Efficient motif discovery for large-scale time series in healthcare[J].IEEE Transactions on Industrial Informatics. 2015，11（3）：583-590.

[4] Rajagopal N，et al.High-throughput mapping of regulatory DNA[J].Nature biotechnology，2016，34（2）：167.

[5] Pevzner PA，Sze SH.Combinatorial approaches to finding subtle signals in DNA sequences[C].In：Proceedings of the Eighth International Conference on Intelligent Systems for Molecular Biology，California，2000：269-278.

[6] Evans PA，Smith AD，Wareham HT.On the complexity of finding common approximate substrings[J].Theory Computer Science，2003，306：407-430.

[7] Zambelli F，Pesole G，Pavesi G.Motif discovery and transcription factor binding sites before and after the next generation sequencing era[J].Briefings in Bioinformatics，2013，14

（2）：225-237.

[8] Al-Okaily A，Huang CH.ET-motif：solving the exact（l，d）-planted motif problem using error tree structure[J].Journal of Computational Biology，2016，23（7）：615-623.

[9] Nicolae M，Rajasekaran S.qPMS9：an efficient algorithm for querum planted motif search[J].Scientific Reports，2015，5：7813.

[10] Machanick P，Bailey TL.MEME-ChIP： motif analysis of large DNA datasets[J]. Bioinformatics，2011，27（12）：1696-1697.

[11] Boucher C，Brown D，Church P.A graph clustering approach to weak motif recognition[C].In：Proceedings of the 7th International Workshop on Algorithms in Bioinformatics，Philadelphia，PA，USA，2007：149-160.

[12] 高山，等.下一代測(cè)序中ChIP-seq數(shù)據(jù)的處理與分析[J].遺傳， 2012，34（6）：773-783.

作者簡(jiǎn)介：

劉 暢（1974-），男，本科，工程師.研究領(lǐng)域：軟件開(kāi)發(fā).endprint