有效的Common Motif 識別算法

木妮娜·玉素甫++古麗娜·玉素甫

摘要：模體發(fā)現(xiàn)在揭示基因組水平上的基因表達調(diào)控規(guī)律以及在蛋白質(zhì)序列中定位保守結(jié)構(gòu)域中起著重要作用。本文提出一種在生物序列中識別Common Motif（公共模體）的算法。算法采用基于后綴數(shù)組或QSA數(shù)組的重復(fù)模式識別算法挖掘串中最大重復(fù)模式作為基元，對基元進行過濾與剪枝后，根據(jù)約束條件對優(yōu)化后基元進行計算與處理從而得到公共模體。算法與基于后綴樹或Trie樹的同類算法相比在時間和空間效率上都得到了提高。

關(guān)鍵詞： 模體發(fā)現(xiàn)；重復(fù)模式；約束條件；生物計算；后綴數(shù)組

中圖分類號：TP391 文獻標識碼：A 文章編號：1009-3044（2016）11-0164-05

Abstract： Motif finding plays an important role on revealing the regulation of gene expression in the genomic level and targeting the conserved domains in the protein sequence. This paper presents an algorithm for finding Common Motif in biological sequences. The algorithm uses the repeat detection algorithms which based on suffix array or QSA array to mining the maximal repeats as primitives. After filtering and pruning， optimized primitives are calculated and processed according to constraints to obtain the common motif. The algorithm is more time and space efficient than the algorithms based on suffix tree or Trie.

Key words： Motif finding； repeats； constraints； bioinformatics； suffix array

1 概述

揭示基因組水平上的基因表達調(diào)控規(guī)律是現(xiàn)代生物信息學(xué)面臨的重大挑戰(zhàn)之一，模體發(fā)現(xiàn)（Motif Finding）問題被認為是這一挑戰(zhàn)中的重要任務(wù)[1]，其出發(fā)點是找出序列中或不同序列間的相似性片段，從而歸結(jié)出序列片段中蘊涵的特征模式，進而推斷出該特征模式與已知的結(jié)構(gòu)和功能之間的內(nèi)在聯(lián)系。這種具有保守特征的序列模式，即相似性高的序列片段，通常稱為模體（Motif）[2]。Motif不僅在核酸序列中存在，在蛋白質(zhì)序列中也存在這種特征模式的序列片段。模體發(fā)現(xiàn)方法在核酸序列分析中發(fā)揮著重要作用[3]，如啟動子識別、DNA結(jié)合蛋白質(zhì)的靶基因和靶位的識別等。對于蛋白質(zhì)序列而言，它們通常表明蛋白質(zhì)序列中特異性結(jié)合位點，如核酸酶和轉(zhuǎn)錄因子（TF）。模體發(fā)現(xiàn)可以利用蛋白質(zhì)序列或結(jié)構(gòu)中的某些特征模式識別相關(guān)蛋白質(zhì)的性質(zhì)。另外，模體發(fā)現(xiàn)還可以用于某些相似性較低的相關(guān)序列檢測。因此，分析和識別Motif及了解它們的功能對于理解和解釋整個基因組行為的意義重大。

現(xiàn)有文獻中對生物序列中Motif有不同的定義，比如序列模體（Sequence motifs）、結(jié)構(gòu)（化）模體（Structured motifs， Structural motifs）、網(wǎng)絡(luò)模體（Network motifs）、公共模體/復(fù)合模體（Common motifs）等分類。由于模體的特征差別很大，不同的方法適用于不同的問題，而且測試序列的選擇也十分困難，因而對于眾多的模體識別方法與算法并沒有統(tǒng)一的評價標準。文獻[4]對模體識別算法進行了總結(jié)，可以從中看出研究者已經(jīng)在此領(lǐng)域取得了很多成果，隨著人們對序列的生物學(xué)意義的不斷了解，新的方法和算法也將會不斷產(chǎn)生。

一般地，根據(jù)算法搜索策略與設(shè)計中所用的組合方法的不同，查找Motif的算法粗略地分為兩類，即第一類是概率序列模型的方法，又稱為統(tǒng)計方法，典型的有基于期望最大化（Expectation Maximization） [5]，采樣（Sampling） [6]，隨機投影（Random Projections） [7]等技術(shù)的方法。第二類是基于串（詞/模式）的方法，又稱為組合方法。其中包括簡單字串列舉法YMF（Yeast Motif Finder）[8]，不匹配（前綴）樹法MITRA[9]以及WINNOWER算法[10]等。

Marsan和Sagot[11]提出了兩種基于后綴樹的算法用于查找DNA序列中保守結(jié)構(gòu)模體，考慮到了結(jié)合位點的結(jié)構(gòu)，允許一定數(shù)量的不匹配堿基位。在文獻[12，13，14]中，作者提出了挖掘帶有間隙的公共模體的算法，其中文獻[12]利用自動機構(gòu)造算法，并基于模體在所有序列中都存在的設(shè)定，而[13，14]中算法則是基于后綴樹。

盡管已有線性時間和空間的創(chuàng)建后綴樹算法，但是對于長度為n的文本，建立后綴樹一般需要20n～40n，建樹過程中頻繁的內(nèi)存分配操作也導(dǎo)致速度很慢。在信息過濾和計算生物學(xué)等領(lǐng)域中，由于一條數(shù)據(jù)信息的長度可能達到幾個G，對于這種超長的數(shù)據(jù)信息，如果構(gòu)造后綴樹結(jié)構(gòu)，則所需的存儲空間是其應(yīng)用的一個瓶頸。如文獻[14]中算法還需建立兩次后綴樹，即對正向與反向字串各建立一次。文獻[13]中算法建立后綴樹后，經(jīng)搜索獲取深度為k的子樹以標示基元。DNA或蛋白質(zhì)序列是自相似度較高序列，序列的自相似度越高，就有越多的模體，因而運行時的速度及輸出的規(guī)模大小也都是模體發(fā)現(xiàn)算法需要考慮的問題。

后綴數(shù)組是一種解決字符串問題的有力工具，相比于后綴樹，它更易于實現(xiàn)且占用內(nèi)存更少，因而使用后綴數(shù)組比使用其他索引結(jié)構(gòu)（后綴樹、字典樹等）更為高效。在本文中我們提出一種公共模體識別的有效算法，使用基于后綴數(shù)組或QSA數(shù)組的重復(fù)模式識別算法挖掘串中最大重復(fù)模式作為基元；對基元進行過濾與剪枝后，根據(jù)約束條件對優(yōu)化后基元進行計算與處理從而得到公共模體。通過理論分析，算法在時間和空間效率上都比同類算法得到了提高。

2 基本定義及挖掘Common Motif的模型

給定一個有限字符集Σ，Σ上任一條序列S可看作是有限個字符順次排列形成的字符串，|S|稱為S的長度。本文中Σ表示生物序列字符集，若為DNA序列，則Σ={A，C，G，T}；若為RNA序列Σ={A，C，G，U}；若為蛋白質(zhì)序列，則Σ為20個簡單氨基酸分子。S[i]表示S的第i個元素，其中1≤i≤|S|。S[i，j]表示S中一個子串，其起始位置為i，結(jié)束位置為j，其中1≤i≤j≤|S|，并且|S[i，j]|=j–i+1。S從i開始的后綴表示為Suffix（i），即Suffix（i）=S[i，|S|]。

我們考慮具有間隙的公共模體的挖掘問題，定義如下。

定義1（Common Motif，公共模體）給定串集合{S1， S2， ... ， Sr}及正整數(shù)k，m，d，其中d≥1，k≥1，m>1，查找具有間隙的Common Motif問題就是查找具有最大長度的非空子串B1， B2 ， ... ，Bm，并且每個Bi（1≤i≤m）滿足以下條件：

1. [B1=B2=...=Bm=k；]

2. [B1*d1B2*d2....*dm-1Bm]出現(xiàn)在Si中（i = 1..r）。

在定義1中，如果d為零，則[B1*d1B2*d2....*dm-1Bm]成為Si中精確重復(fù)模式，因此設(shè)d≥1；同時設(shè)k≥1，m>1。我們稱[B1*d1B2*d2....*dm-1Bm]為帶有間隙的模體，由于出現(xiàn)在多個串Si中，又稱為公共模體。子串組B1， B2 ， ... ，Bm，又稱為“鏈”（Chain），每個Bi稱為帶有間隙的模體的“塊”（Block），每個Bi的長度是固定的，且長度相等。這里的d1，d2，…，dm-1如果不同，則為可變間隙，本文我們只討論帶有固定間隙的公共模體。

3 Common Motif識別算法

3.1 重復(fù)模式的形式化描述

有許多方法可以檢測一條序列中重復(fù)的序列模式，這些方法經(jīng)過改進后可用于從一組序列中提取相同或者相近的短序列，作為候選的Motif。因而查找具有間隙的公共模體問題，類似于查找序列中重復(fù)模式，其中各個重復(fù)模式子串之間帶有連續(xù)的通配符，稱為間隙。一個Motif被認為是最大的，如果它既不能向左擴展也不能向右擴展。如果是公共模體，根據(jù)定義1，非空子串B1， B2 ， ... ，Bm具有最大長度，即組成模體的各個子串是最大的。在本文算法中將B1， B2 ， ... ，Bm稱為基元。

文獻[15]給出了最大重復(fù)模式完整定義，我們現(xiàn)結(jié)合本文算法對其中幾個定義作一些格式上的修改并簡單描述如下：

定義1：設(shè)u為串S的一個子串，則重復(fù)模式（Repeats）是一個多元組： R = （ p； u； i1， i2，…， ie ），其中p≥1， e≥2；1≤i1≤i2≤…≤ie ≤n；u= S[i1… i1 + p-1] = S[i2… i2 + p-1] = … = S[ie … ie + p-1]；稱p為R的周期（長度）；e為R的指數(shù)；u為R的生成元，也稱為S的重復(fù)子串，i1，i2，…，ie為重復(fù)子串起始位置。

定義2：不可左擴展（NLE）：至少存在一對s，t （1≤s

定義3：不可右擴展（NRE）：至少存在一對s，t （1≤s

定義4：如果一個Repeat既是NLE又是NRE，則稱其NE，或最大（Maximal Repeats）。

3.2 Common Motif識別算法步驟

Common Motif識別算法包括以下6個步驟：

步驟1：

1）對于r個不同的串集合{S1， S2， . . . ， Sr}，我們首先將其合并成一個串，各串之間用字符#1，#2，…，#r-1分隔，即S=S1#1S2#2…#r-1Sr，其中#j（1≤j≤r-1）不屬于任何串，且規(guī)定其值小于Σ中任意字符，并且滿足條件#1<#2<…<#r-1。設(shè)|Si|=ni（1≤i≤r），則|S|=n，其中n為所有ni平均值。

2）利用已有算法（或作相應(yīng)修改），識別S中所有最大重復(fù)模式。

根據(jù)定義，組成Common Motif中的“塊”即基元是最大重復(fù)模式，因而Common Motif識別的關(guān)鍵問題就是查找Maximal Repeats。最大重復(fù)模式的識別也是公共模體識別過程中最耗時的一步?，F(xiàn)有的很多方法選用基于后綴樹的算法來查找重復(fù)子序列，而根據(jù)第1節(jié)中分析，后綴數(shù)組相比于后綴樹，更易于實現(xiàn)且占用內(nèi)存更少。因此在此步驟，我們利用文獻[15]或[16]中算法。在文獻[15]中，我們提出了一種基于QSA數(shù)組計算所有帶有約束條件的NE重復(fù)模式（亦即最大重復(fù)模式）的算法RPT。RPT算法的優(yōu)點是可通過約束條件最小周期pmin和最大間距gmax，篩選符合條件的長度為n的串中NE重復(fù)模式，算法最多使用略大于5n字節(jié)的內(nèi)存來存儲序列S、數(shù)組QSA及位矢量PROC，因而空間效率很高，缺點是最差時間復(fù)雜度略高。在文獻[16]中，我們提出識別長度為n生物序列中完全重復(fù)模式的有效算法CRFinder，算法基于后綴數(shù)組SA[17]和最長公共前綴數(shù)組LCP[18]，因而空間復(fù)雜度是9n。算法經(jīng)過改寫即可輸出最大重復(fù)模式。如圖1所示，例如對蛋白質(zhì)序列S=ATGCAATGCCVGGCATTGCATV，算法CRFinder首先計算S的SA與LCP值，然后根據(jù)LCP的上升與下降規(guī)律，通過入棧與出棧操作，識別S中重復(fù)模式。則改寫CRFinder需進一步作：1）在入棧與出棧操作時，只輸出NRE重復(fù)模式；2）對輸出后的重復(fù)模式做NLE檢測。

可以使用線性時間構(gòu)造算法計算后綴數(shù)組SA（如使用算法[17]）和最長公共前綴數(shù)組LCP（如使用算法[18]），因而對于串S及用戶給定的正整數(shù)pmin與fmin，我們可以通過改寫的CRFinder算法識別長度p≥pmin的所有完整最大重復(fù)模式。算法執(zhí)行時間與字符集無關(guān)，是串長度的線性時間，即為O（n）。算法輸出格式為（p；u；i，j，?），其中p和u與定義1中參數(shù)相同；?為滿足條件?≥?min的模式頻率，i…j是SA中一個區(qū)間。這種輸出表示法的缺點是重復(fù)子串的出現(xiàn)位置不是順序的，需要先利用SA數(shù)組將（p；u；i，j，?）轉(zhuǎn)換為（p；u；SA[i]，SA[i+1]，…SA[j]，?）的格式，并采用排序法對所有最大重復(fù)模式子串的出現(xiàn)位置（即SA[i]，SA[i+1]，…SA[j]）進行排序，因此會增加算法的時間復(fù)雜度。所以我們在步驟1（3）中設(shè)計檢測算法Check，既避免對重復(fù)子串進行排序，同時也可檢測其是否滿足限制條件。

3）判斷步驟1（2）中識別出的重復(fù)模式是否滿足限制條件：至少在每個串中出現(xiàn)mmin次且至少在q個串中出現(xiàn)，如表1中算法所示。

數(shù)組count=count[1..r]對位于各個Si中的重復(fù)子串進行計數(shù)。通過使用這兩個數(shù)組，可用O（（j-i）log2r）時間復(fù)雜度的算法計算出R中重復(fù)子串是否滿足滿足限制條件：至少在每個串中出現(xiàn)mmin次且至少在q個串中出現(xiàn)。值得注意的是，算法首先檢測重復(fù)模式出現(xiàn)頻率?（即j?i+1），如果小于mminq，則沒有必要進行進一步檢測。函數(shù)BinarySearch返回索引t，表明重復(fù)子串位置SA[h]出現(xiàn)在串St中。設(shè)最大重復(fù)模式個數(shù)為O（α），對每個重復(fù)模式，數(shù)組count需清零r次，因而時間復(fù)雜度為O（αr）。對S中每個重復(fù)子串，二分查找最壞時間為O（log2r），設(shè)S中識別出的重復(fù)子串總數(shù)為φ個，則算法Check的時間復(fù)雜度為O（αr+φlog2r）。算法Check還可進一步簡化，如設(shè)置一個初始值為零的鏈表L，將二分查找返回值t加入進去；在最后一條語句output（R）后執(zhí)行另一個循環(huán)，將L中t值清零，并執(zhí)行count[t]=0；則算法時間復(fù)雜度可降為O（φlog2r）。

4）對S中所有的重復(fù)模式提取出第二字段即子串部分，作為公共模體的基元，并進行編號；編號可以依據(jù)輸出順序，也可以隨機編號。

步驟2：建立一個新的串集合T={T1，T2，…，Tr}，其中Ti （1≤i≤r）的值為步驟1（4）中取得的值，即Ti[j]=Si中重復(fù)模式在j位置上的值。

步驟3：在Ti （1≤i≤r）中做m-1次長度為d+k的跳越搜索，查找所有長度為m的“塊”的鏈；即

[Ci[j]=Ti[j]Ti[j+d+k]Ti[j+2（d+k）]...Ti[j+（m-1）（d+k）]，] [?j∈1..n-m（d+k）+d]

注意，只要有一項Ti的值為零，則Ci的值為零。

從而使得motif具有如下長度：

[B1*dB2*d...*dBm=m（d+k）-d]

步驟4：設(shè)Ri為鏈Ci[j]的集合，其中[?j∈1..n-m（d+k）+d]，并對所有的i求出Ri。

步驟5：對集合Ri（為鏈Ci[j]集合的集合）建立Motif樹（即另一種形式的廣義后綴樹），深度為m，葉子節(jié)點為Ci[j]所在的Ri。

步驟6：含有公共Ri的葉子節(jié)點的路徑，即為公共模體。

3.3 算法舉例

下面根據(jù)條件k=2， m=3， d=1，以圖2中S1，S2，S3為例，根據(jù)步驟1-6識別S1，S2，S3中公共模體。

3.4 算法的性能分析

現(xiàn)在分析算法的復(fù)雜度。第一步驟已在前面分析，為O（rn+αr+φlog2r），經(jīng)過改進可降為O（rn+φlog2r）；步驟2的時間復(fù)雜度為O（rn）；在步驟3中，創(chuàng)建了長度為m（d+k）-d的n-m（d+k）+d個鏈，因此步驟3的最差時間復(fù)雜度為O（rn），假設(shè)m，k和d為常數(shù)。步驟5中建立Motif樹最差時間復(fù)雜度O（rn），步驟6也是線性時間。因而總的時間復(fù)雜度為O（rn+φlog2r）。

3.5 與現(xiàn)有算法的比較與改進

本文研究工作與文獻[13]的研究工作（問題1的算法）最為接近，其區(qū)別在于：

1）由于采用了基于后綴數(shù)組或QSA數(shù)組的重復(fù)模式識別算法挖掘串中最大重復(fù)模式作為基元，空間效率得到很大提高，與建立后綴樹需20rn～40rn相比，采用基于QSA數(shù)組的算法空間復(fù)雜度為6rn，基于后綴數(shù)組算法空間復(fù)雜度為9rn。同時也避免了在廣義后綴樹上搜索得到k-子樹的過程。

2）由于采用步驟1（4），可判斷基元是否在所有序列中出現(xiàn)，一方面通過剪枝，在后續(xù)步驟中（步驟4和5）避免了大量的匹配操作，可有效限制搜索空間，極大地提高了實際運行時空效率。另一方面，也可回答下列問題，即“給定的一組序列，可能的motif僅在部分序列中出現(xiàn)，如何解決？”，此時可取q

3）算法可以很方便地進行擴展，對不同間隙的Common Motif進行識別；由于篇幅所限，本文沒有討論此問題。

4 結(jié)束語

本文提出了一種公共模體識別的有效算法，算法與基于后綴樹或Trie樹的同類算法相比在時間和空間效率上都得到了提高。

參考文獻：

[1] Kellis， M.， Patterson， N.， Endrizzi， M.， Birren， B. & Lander， E. S. Sequencing and comparison of yeast species to identify genes and regulatory elements[J]. Nature， 2003， 423：241-254.

[2] Patrik D'haeseleer1. What are DNA sequence motifs？[J]. Nature Biotechnology， 2006， 24：423-425.

[3] GuhaThakurta D. Computational identification of transcriptional regulatory elements in DNA sequence[J]. Nucleic Acids Res. 2006， 34（12）： 3585-3598.

[4] Modan K Das， Ho-Kwok Dai. A survey of DNA motif finding algorithms[J]. BMC Bioinformatics， 2007，8（Suppl 7）：S21.

[5] T. L. Bailey， C. Elkan. Unsupervised Learning of Multiple Motifs in Biopolymers using EM[J]. Machine Learning， 1995， 21（1-2）：51-80.

[6] W. Thompson， E. C. Rouchka， C. E. Lawrence. Gibbs Recursive Sampler： Finding Transcription Factor Binding Sites[J]. Nucleic Acids Research， 2003， 31（13）：3580-3585.

[7] J. Buhler， M. Tompa. Finding Motifs Using Random Projections[J]. Journal Computational Biology， 2002， 9（2）：225-242.

[8] S. Sinha， M. Tompa. YMF： A Program for Discovery of Novel Transcription Factor Binding Sites by Statistical Overrepresentation[J]. Nucleic Acids Research， 2003， 31（13）： 3586-3588.

[9] E. Eskin， P. A. Pevzner. Finding Composite Regulatory Patterns in DNA Sequences[J]. Bioinformatics， 2002， 18 Suppl 1：S354-363.

[10] Pevzner P， Sze S. Combinatorial approaches to finding subtle signals in DNA sequences[C]. In Proceedings of the Eighth International Conference on Intelligent Systems for Molecular Biology. Menlo Park， California： AAAI Press， 2000： 269-278.

[11] L. Marsan， M.-F. Sagot. Algorithms for Extracting Structured Motifs Using a Suffix Tree with Application to Promoter and Regulatory Site Consensus Identification[J]. Journal of Computational Biology， 2000， 7（3-4）： 345-360.

[12] Pavlos Antoniou， Jan Holub， Costas S. Iliopoulos， Bo ivoj Melichar， Pierre Peterlongo. Finding Common Motifs with Gaps Using Finite Automata[J]. Lecture Notes in Computer Science， 2006， 4094：69-77.

[13] C. S. Iliopoulos， J. Mchugh， P. Peterlongo， N. Pisanti， W. Rytter， M.-F. Sagot. A first approach to finding common motifs with gaps[J]. International Journal of Foundations of Computer Science， 2005， 16（6）：1145-1154.

[14] P. Antoniou， M. Crochemore， C. Iliopoulos， P. Peterlongo. Application of suffix trees for the acquisition of common motifs with gaps in a set of strings[C]. Proceedings of the International Conference on Language and Automata Theory and Applications， 2007， 88-97.

[15] 木妮娜·玉素甫，古麗娜·玉素甫，張海軍.基于QSA數(shù)組計算序列中所有NE重復(fù)模式的算法[J]. 計算機科學(xué)， 2014， 41（3）： 249-252-262.

[16] Munina Yusufu， Gulina Yusufu. Efficient Algorithm for Extracting Complete Repeats from Biological Sequences[J]. International Journal of Computer Applications， 2015， 128（16）：33-37.

[17] Juha Karkkainen， Peter Sanders. Simple linear work suffix array construction[C]， Proc. of 30th ICALP， LNCS 2719， 2003， 943-955.

[18] Kasai， G. Lee， H. Arimura， S. Arikawa， K. Park. Linear-time longest-common-prefix computation in suffix arrays and its applications[C]， Proc. of 12th CPM， LNCS 2089， 2001， 181-192.