劉暢

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

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

中圖分類號：TP301.6 文獻(xiàn)標(biāo)識碼：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）

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

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

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

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

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

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

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

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

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

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

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

3 方法（Method）

3.1 算法思想

MCL-WMR存在如下的兩個缺陷。

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

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

為了解決這兩個問題，本文對MCL-WMR進(jìn)行改進(jìn)，提出一個稱為iMCL-WMR的算法。基本思想如下：評估輸入序列中每個l-mer的信號強(qiáng)度，將信號較強(qiáng)的l-mer視為潛在的模體實例；由每個潛在的模體實例誘導(dǎo)一個子圖G'，使得G'中有邊相連的兩個l-mer間的海明距離小于等于d；對G'進(jìn)行聚類，對得到的各個聚類進(jìn)行求精，最后選擇得分最大的進(jìn)行輸出。

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

3.2 劃分子序列集合

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

（1）

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

（2）

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

計算出輸入序列中所有的l-mer si，j的信號強(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所示，橫軸表示信號強(qiáng)度，縱軸表示對應(yīng)信號強(qiáng)度的l-mer的數(shù)量?；诖?，將L劃分成L1、L2和L3三部分。

（3）

（4）

（5）

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

3.3 聚類求精

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

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

求精算法迭代地對當(dāng)前的模體進(jìn)行更新，直到得到信息量最高的模體。模體c的信息量記為IC（c），計算方法如下所示，其中fij表示C中l(wèi)-mer對齊后第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行偽代碼計算輸入序列中各個l-mer的信號強(qiáng)度。依據(jù)信號強(qiáng)度，第2行偽代碼將輸入序列中的l-mer劃分成L1、L2和L3三個集合。第3行偽代碼使用L2和L3中的l-mer構(gòu)建一個圖G，注意iMCL-WMR要求兩個l-mer間的海明距離k小于等于d時對應(yīng)的兩個頂點才有邊相連接；邊的權(quán)值設(shè)置為（d-k）（l-k），這比MCL-WMR的權(quán)值設(shè)置方法更為精細(xì)，使得海明距離越小時越以更大的概率視為一對模體實例。第4—9行偽代碼對L3中的各個l-mer誘導(dǎo)的子圖依次進(jìn)行聚類求精，得到多個模體。第10—11行偽代碼對得到的多個模體依照其得分進(jìn)行排序，最終輸出得分最高的模體。

4 實驗結(jié)果（Experimental results）

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

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

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

5 結(jié)論（Conclusion）

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

參考文獻(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] 高山，等.下一代測序中ChIP-seq數(shù)據(jù)的處理與分析[J].遺傳， 2012，34（6）：773-783.

作者簡介：

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