论文部分内容阅读
摘要：模体发现在揭示基因组水平上的基因表达调控规律以及在蛋白质序列中定位保守结构域中起着重要作用。本文提出一种在生物序列中识别Common Motif（公共模体）的算法。算法采用基于后缀数组或QSA数组的重复模式识别算法挖掘串中最大重复模式作为基元，对基元进行过滤与剪枝后，根据约束条件对优化后基元进行计算与处理从而得到公共模体。算法与基于后缀树或Trie树的同类算法相比在时间和空间效率上都得到了提高。
关键词： 模体发现；重复模式；约束条件；生物计算；后缀数组
中图分类号：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 概述
揭示基因组水平上的基因表达调控规律是现代生物信息学面临的重大挑战之一，模体发现（Motif Finding）问题被认为是这一挑战中的重要任务[1]，其出发点是找出序列中或不同序列间的相似性片段，从而归结出序列片段中蕴涵的特征模式，进而推断出该特征模式与已知的结构和功能之间的内在联系。这种具有保守特征的序列模式，即相似性高的序列片段，通常称为模体（Motif）[2]。Motif不仅在核酸序列中存在，在蛋白质序列中也存在这种特征模式的序列片段。模体发现方法在核酸序列分析中发挥着重要作用[3]，如启动子识别、DNA结合蛋白质的靶基因和靶位的识别等。对于蛋白质序列而言，它们通常表明蛋白质序列中特异性结合位点，如核酸酶和转录因子（TF）。模体发现可以利用蛋白质序列或结构中的某些特征模式识别相关蛋白质的性质。另外，模体发现还可以用于某些相似性较低的相关序列检测。因此，分析和识别Motif及了解它们的功能对于理解和解释整个基因组行为的意义重大。
现有文献中对生物序列中Motif有不同的定义，比如序列模体（Sequence motifs）、结构（化）模体（Structured motifs， Structural motifs）、网络模体（Network motifs）、公共模体/复合模体（Common motifs）等分类。由于模体的特征差别很大，不同的方法适用于不同的问题，而且测试序列的选择也十分困难，因而对于众多的模体识别方法与算法并没有统一的评价标准。文献[4]对模体识别算法进行了总结，可以从中看出研究者已经在此领域取得了很多成果，随着人们对序列的生物学意义的不断了解，新的方法和算法也将会不断产生。
一般地，根据算法搜索策略与设计中所用的组合方法的不同，查找Motif的算法粗略地分为两类，即第一类是概率序列模型的方法，又称为统计方法，典型的有基于期望最大化（Expectation Maximization） [5]，采样（Sampling） [6]，随机投影（Random Projections） [7]等技术的方法。第二类是基于串（词/模式）的方法，又称为组合方法。其中包括简单字串列举法YMF（Yeast Motif Finder）[8]，不匹配（前缀）树法MITRA[9]以及WINNOWER算法[10]等。
Marsan和Sagot[11]提出了两种基于后缀树的算法用于查找DNA序列中保守结构模体，考虑到了结合位点的结构，允许一定数量的不匹配碱基位。在文献[12，13，14]中，作者提出了挖掘带有间隙的公共模体的算法，其中文献[12]利用自动机构造算法，并基于模体在所有序列中都存在的设定，而[13，14]中算法则是基于后缀树。
尽管已有线性时间和空间的创建后缀树算法，但是对于长度为n的文本，建立后缀树一般需要20n～40n，建树过程中频繁的内存分配操作也导致速度很慢。在信息过滤和计算生物学等领域中，由于一条数据信息的长度可能达到几个G，对于这种超长的数据信息，如果构造后缀树结构，则所需的存储空间是其应用的一个瓶颈。如文献[14]中算法还需建立两次后缀树，即对正向与反向字串各建立一次。文献[13]中算法建立后缀树后，经搜索获取深度为k的子树以标示基元。DNA或蛋白质序列是自相似度较高序列，序列的自相似度越高，就有越多的模体，因而运行时的速度及输出的规模大小也都是模体发现算法需要考虑的问题。
后缀数组是一种解决字符串问题的有力工具，相比于后缀树，它更易于实现且占用内存更少，因而使用后缀数组比使用其他索引结构（后缀树、字典树等）更为高效。在本文中我们提出一种公共模体识别的有效算法，使用基于后缀数组或QSA数组的重复模式识别算法挖掘串中最大重复模式作为基元；对基元进行过滤与剪枝后，根据约束条件对优化后基元进行计算与处理从而得到公共模体。通过理论分析，算法在时间和空间效率上都比同类算法得到了提高。 2 基本定义及挖掘Common Motif的模型
给定一个有限字符集Σ，Σ上任一条序列S可看作是有限个字符顺次排列形成的字符串，|S|称为S的长度。本文中Σ表示生物序列字符集，若为DNA序列，则Σ={A，C，G，T}；若为RNA序列Σ={A，C，G，U}；若为蛋白质序列，则Σ为20个简单氨基酸分子。S[i]表示S的第i个元素，其中1≤i≤|S|。S[i，j]表示S中一个子串，其起始位置为i，结束位置为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}及正整数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]出现在Si中（i = 1..r）。
在定义1中，如果d为零，则[B1*d1B2*d2....*dm-1Bm]成为Si中精确重复模式，因此设d≥1；同时设k≥1，m>1。我们称[B1*d1B2*d2....*dm-1Bm]为带有间隙的模体，由于出现在多个串Si中，又称为公共模体。子串组B1， B2 ， ... ，Bm，又称为“链”（Chain），每个Bi称为带有间隙的模体的“块”（Block），每个Bi的长度是固定的，且长度相等。这里的d1，d2，…，dm-1如果不同，则为可变间隙，本文我们只讨论带有固定间隙的公共模体。
3 Common Motif识别算法
3.1 重复模式的形式化描述
有许多方法可以检测一条序列中重复的序列模式，这些方法经过改进后可用于从一组序列中提取相同或者相近的短序列，作为候选的Motif。因而查找具有间隙的公共模体问题，类似于查找序列中重复模式，其中各个重复模式子串之间带有连续的通配符，称为间隙。一个Motif被认为是最大的，如果它既不能向左扩展也不能向右扩展。如果是公共模体，根据定义1，非空子串B1， B2 ， ... ，Bm具有最大长度，即组成模体的各个子串是最大的。在本文算法中将B1， B2 ， ... ，Bm称为基元。
文献[15]给出了最大重复模式完整定义，我们现结合本文算法对其中几个定义作一些格式上的修改并简单描述如下：
定义1：设u为串S的一个子串，则重复模式（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的指数；u为R的生成元，也称为S的重复子串，i1，i2，…，ie为重复子串起始位置。
定义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）不属于任何串，且规定其值小于Σ中任意字符，并且满足条件#1<#2<…<#r-1。设|Si|=ni（1≤i≤r），则|S|=n，其中n为所有ni平均值。
2）利用已有算法（或作相应修改），识别S中所有最大重复模式。
根据定义，组成Common Motif中的“块”即基元是最大重复模式，因而Common Motif识别的关键问题就是查找Maximal Repeats。最大重复模式的识别也是公共模体识别过程中最耗时的一步。现有的很多方法选用基于后缀树的算法来查找重复子序列，而根据第1节中分析，后缀数组相比于后缀树，更易于实现且占用内存更少。因此在此步骤，我们利用文献[15]或[16]中算法。在文献[15]中，我们提出了一种基于QSA数组计算所有带有约束条件的NE重复模式（亦即最大重复模式）的算法RPT。RPT算法的优点是可通过约束条件最小周期pmin和最大间距gmax，筛选符合条件的长度为n的串中NE重复模式，算法最多使用略大于5n字节的内存来存储序列S、数组QSA及位矢量PROC，因而空间效率很高，缺点是最差时间复杂度略高。在文献[16]中，我们提出识别长度为n生物序列中完全重复模式的有效算法CRFinder，算法基于后缀数组SA[17]和最长公共前缀数组LCP[18]，因而空间复杂度是9n。算法经过改写即可输出最大重复模式。如图1所示，例如对蛋白质序列S=ATGCAATGCCVGGCATTGCATV，算法CRFinder首先计算S的SA与LCP值，然后根据LCP的上升与下降规律，通过入栈与出栈操作，识别S中重复模式。则改写CRFinder需进一步作：1）在入栈与出栈操作时，只输出NRE重复模式；2）对输出后的重复模式做NLE检测。 可以使用线性时间构造算法计算后缀数组SA（如使用算法[17]）和最长公共前缀数组LCP（如使用算法[18]），因而对于串S及用户给定的正整数pmin与fmin，我们可以通过改写的CRFinder算法识别长度p≥pmin的所有完整最大重复模式。算法执行时间与字符集无关，是串长度的线性时间，即为O（n）。算法输出格式为（p；u；i，j，?），其中p和u与定义1中参数相同；?为满足条件?≥?min的模式频率，i…j是SA中一个区间。这种输出表示法的缺点是重复子串的出现位置不是顺序的，需要先利用SA数组将（p；u；i，j，?）转换为（p；u；SA[i]，SA[i 1]，…SA[j]，?）的格式，并采用排序法对所有最大重复模式子串的出现位置（即SA[i]，SA[i 1]，…SA[j]）进行排序，因此会增加算法的时间复杂度。所以我们在步骤1（3）中设计检测算法Check，既避免对重复子串进行排序，同时也可检测其是否满足限制条件。
3）判断步骤1（2）中识别出的重复模式是否满足限制条件：至少在每个串中出现mmin次且至少在q个串中出现，如表1中算法所示。
数组count=count[1..r]对位于各个Si中的重复子串进行计数。通过使用这两个数组，可用O（（j-i）log2r）时间复杂度的算法计算出R中重复子串是否满足满足限制条件：至少在每个串中出现mmin次且至少在q个串中出现。值得注意的是，算法首先检测重复模式出现频率?（即j?i 1），如果小于mminq，则没有必要进行进一步检测。函数BinarySearch返回索引t，表明重复子串位置SA[h]出现在串St中。设最大重复模式个数为O（α），对每个重复模式，数组count需清零r次，因而时间复杂度为O（αr）。对S中每个重复子串，二分查找最坏时间为O（log2r），设S中识别出的重复子串总数为φ个，则算法Check的时间复杂度为O（αr φlog2r）。算法Check还可进一步简化，如设置一个初始值为零的链表L，将二分查找返回值t加入进去；在最后一条语句output（R）后执行另一个循环，将L中t值清零，并执行count[t]=0；则算法时间复杂度可降为O（φlog2r）。
4）对S中所有的重复模式提取出第二字段即子串部分，作为公共模体的基元，并进行编号；编号可以依据输出顺序，也可以随机编号。
步骤2：建立一个新的串集合T={T1，T2，…，Tr}，其中Ti （1≤i≤r）的值为步骤1（4）中取得的值，即Ti[j]=Si中重复模式在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：设Ri为链Ci[j]的集合，其中[?j∈1..n-m（d k） d]，并对所有的i求出Ri。
步骤5：对集合Ri（为链Ci[j]集合的集合）建立Motif树（即另一种形式的广义后缀树），深度为m，叶子节点为Ci[j]所在的Ri。
步骤6：含有公共Ri的叶子节点的路径，即为公共模体。
3.3 算法举例
下面根据条件k=2， m=3， d=1，以图2中S1，S2，S3为例，根据步骤1-6识别S1，S2，S3中公共模体。
3.4 算法的性能分析
现在分析算法的复杂度。第一步骤已在前面分析，为O（rn αr φlog2r），经过改进可降为O（rn φlog2r）；步骤2的时间复杂度为O（rn）；在步骤3中，创建了长度为m（d k）-d的n-m（d k） d个链，因此步骤3的最差时间复杂度为O（rn），假设m，k和d为常数。步骤5中建立Motif树最差时间复杂度O（rn），步骤6也是线性时间。因而总的时间复杂度为O（rn φlog2r）。
3.5 与现有算法的比较与改进
本文研究工作与文献[13]的研究工作（问题1的算法）最为接近，其区别在于：
1）由于采用了基于后缀数组或QSA数组的重复模式识别算法挖掘串中最大重复模式作为基元，空间效率得到很大提高，与建立后缀树需20rn～40rn相比，采用基于QSA数组的算法空间复杂度为6rn，基于后缀数组算法空间复杂度为9rn。同时也避免了在广义后缀树上搜索得到k-子树的过程。
2）由于采用步骤1（4），可判断基元是否在所有序列中出现，一方面通过剪枝，在后续步骤中（步骤4和5）避免了大量的匹配操作，可有效限制搜索空间，极大地提高了实际运行时空效率。另一方面，也可回答下列问题，即“给定的一组序列，可能的motif仅在部分序列中出现，如何解决？”，此时可取q 3）算法可以很方便地进行扩展，对不同间隙的Common Motif进行识别；由于篇幅所限，本文没有讨论此问题。
4 结束语
本文提出了一种公共模体识别的有效算法，算法与基于后缀树或Trie树的同类算法相比在时间和空间效率上都得到了提高。
参考文献：
[1] Kellis， M.， Patterson， N.， Endrizzi， M.， Birren， B.
关键词： 模体发现；重复模式；约束条件；生物计算；后缀数组
中图分类号：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 概述
揭示基因组水平上的基因表达调控规律是现代生物信息学面临的重大挑战之一，模体发现（Motif Finding）问题被认为是这一挑战中的重要任务[1]，其出发点是找出序列中或不同序列间的相似性片段，从而归结出序列片段中蕴涵的特征模式，进而推断出该特征模式与已知的结构和功能之间的内在联系。这种具有保守特征的序列模式，即相似性高的序列片段，通常称为模体（Motif）[2]。Motif不仅在核酸序列中存在，在蛋白质序列中也存在这种特征模式的序列片段。模体发现方法在核酸序列分析中发挥着重要作用[3]，如启动子识别、DNA结合蛋白质的靶基因和靶位的识别等。对于蛋白质序列而言，它们通常表明蛋白质序列中特异性结合位点，如核酸酶和转录因子（TF）。模体发现可以利用蛋白质序列或结构中的某些特征模式识别相关蛋白质的性质。另外，模体发现还可以用于某些相似性较低的相关序列检测。因此，分析和识别Motif及了解它们的功能对于理解和解释整个基因组行为的意义重大。
现有文献中对生物序列中Motif有不同的定义，比如序列模体（Sequence motifs）、结构（化）模体（Structured motifs， Structural motifs）、网络模体（Network motifs）、公共模体/复合模体（Common motifs）等分类。由于模体的特征差别很大，不同的方法适用于不同的问题，而且测试序列的选择也十分困难，因而对于众多的模体识别方法与算法并没有统一的评价标准。文献[4]对模体识别算法进行了总结，可以从中看出研究者已经在此领域取得了很多成果，随着人们对序列的生物学意义的不断了解，新的方法和算法也将会不断产生。
一般地，根据算法搜索策略与设计中所用的组合方法的不同，查找Motif的算法粗略地分为两类，即第一类是概率序列模型的方法，又称为统计方法，典型的有基于期望最大化（Expectation Maximization） [5]，采样（Sampling） [6]，随机投影（Random Projections） [7]等技术的方法。第二类是基于串（词/模式）的方法，又称为组合方法。其中包括简单字串列举法YMF（Yeast Motif Finder）[8]，不匹配（前缀）树法MITRA[9]以及WINNOWER算法[10]等。
Marsan和Sagot[11]提出了两种基于后缀树的算法用于查找DNA序列中保守结构模体，考虑到了结合位点的结构，允许一定数量的不匹配碱基位。在文献[12，13，14]中，作者提出了挖掘带有间隙的公共模体的算法，其中文献[12]利用自动机构造算法，并基于模体在所有序列中都存在的设定，而[13，14]中算法则是基于后缀树。
尽管已有线性时间和空间的创建后缀树算法，但是对于长度为n的文本，建立后缀树一般需要20n～40n，建树过程中频繁的内存分配操作也导致速度很慢。在信息过滤和计算生物学等领域中，由于一条数据信息的长度可能达到几个G，对于这种超长的数据信息，如果构造后缀树结构，则所需的存储空间是其应用的一个瓶颈。如文献[14]中算法还需建立两次后缀树，即对正向与反向字串各建立一次。文献[13]中算法建立后缀树后，经搜索获取深度为k的子树以标示基元。DNA或蛋白质序列是自相似度较高序列，序列的自相似度越高，就有越多的模体，因而运行时的速度及输出的规模大小也都是模体发现算法需要考虑的问题。
后缀数组是一种解决字符串问题的有力工具，相比于后缀树，它更易于实现且占用内存更少，因而使用后缀数组比使用其他索引结构（后缀树、字典树等）更为高效。在本文中我们提出一种公共模体识别的有效算法，使用基于后缀数组或QSA数组的重复模式识别算法挖掘串中最大重复模式作为基元；对基元进行过滤与剪枝后，根据约束条件对优化后基元进行计算与处理从而得到公共模体。通过理论分析，算法在时间和空间效率上都比同类算法得到了提高。 2 基本定义及挖掘Common Motif的模型
给定一个有限字符集Σ，Σ上任一条序列S可看作是有限个字符顺次排列形成的字符串，|S|称为S的长度。本文中Σ表示生物序列字符集，若为DNA序列，则Σ={A，C，G，T}；若为RNA序列Σ={A，C，G，U}；若为蛋白质序列，则Σ为20个简单氨基酸分子。S[i]表示S的第i个元素，其中1≤i≤|S|。S[i，j]表示S中一个子串，其起始位置为i，结束位置为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}及正整数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]出现在Si中（i = 1..r）。
在定义1中，如果d为零，则[B1*d1B2*d2....*dm-1Bm]成为Si中精确重复模式，因此设d≥1；同时设k≥1，m>1。我们称[B1*d1B2*d2....*dm-1Bm]为带有间隙的模体，由于出现在多个串Si中，又称为公共模体。子串组B1， B2 ， ... ，Bm，又称为“链”（Chain），每个Bi称为带有间隙的模体的“块”（Block），每个Bi的长度是固定的，且长度相等。这里的d1，d2，…，dm-1如果不同，则为可变间隙，本文我们只讨论带有固定间隙的公共模体。
3 Common Motif识别算法
3.1 重复模式的形式化描述
有许多方法可以检测一条序列中重复的序列模式，这些方法经过改进后可用于从一组序列中提取相同或者相近的短序列，作为候选的Motif。因而查找具有间隙的公共模体问题，类似于查找序列中重复模式，其中各个重复模式子串之间带有连续的通配符，称为间隙。一个Motif被认为是最大的，如果它既不能向左扩展也不能向右扩展。如果是公共模体，根据定义1，非空子串B1， B2 ， ... ，Bm具有最大长度，即组成模体的各个子串是最大的。在本文算法中将B1， B2 ， ... ，Bm称为基元。
文献[15]给出了最大重复模式完整定义，我们现结合本文算法对其中几个定义作一些格式上的修改并简单描述如下：
定义1：设u为串S的一个子串，则重复模式（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的指数；u为R的生成元，也称为S的重复子串，i1，i2，…，ie为重复子串起始位置。
定义2：不可左扩展（NLE）：至少存在一对s，t （1≤s
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）不属于任何串，且规定其值小于Σ中任意字符，并且满足条件#1<#2<…<#r-1。设|Si|=ni（1≤i≤r），则|S|=n，其中n为所有ni平均值。
2）利用已有算法（或作相应修改），识别S中所有最大重复模式。
根据定义，组成Common Motif中的“块”即基元是最大重复模式，因而Common Motif识别的关键问题就是查找Maximal Repeats。最大重复模式的识别也是公共模体识别过程中最耗时的一步。现有的很多方法选用基于后缀树的算法来查找重复子序列，而根据第1节中分析，后缀数组相比于后缀树，更易于实现且占用内存更少。因此在此步骤，我们利用文献[15]或[16]中算法。在文献[15]中，我们提出了一种基于QSA数组计算所有带有约束条件的NE重复模式（亦即最大重复模式）的算法RPT。RPT算法的优点是可通过约束条件最小周期pmin和最大间距gmax，筛选符合条件的长度为n的串中NE重复模式，算法最多使用略大于5n字节的内存来存储序列S、数组QSA及位矢量PROC，因而空间效率很高，缺点是最差时间复杂度略高。在文献[16]中，我们提出识别长度为n生物序列中完全重复模式的有效算法CRFinder，算法基于后缀数组SA[17]和最长公共前缀数组LCP[18]，因而空间复杂度是9n。算法经过改写即可输出最大重复模式。如图1所示，例如对蛋白质序列S=ATGCAATGCCVGGCATTGCATV，算法CRFinder首先计算S的SA与LCP值，然后根据LCP的上升与下降规律，通过入栈与出栈操作，识别S中重复模式。则改写CRFinder需进一步作：1）在入栈与出栈操作时，只输出NRE重复模式；2）对输出后的重复模式做NLE检测。 可以使用线性时间构造算法计算后缀数组SA（如使用算法[17]）和最长公共前缀数组LCP（如使用算法[18]），因而对于串S及用户给定的正整数pmin与fmin，我们可以通过改写的CRFinder算法识别长度p≥pmin的所有完整最大重复模式。算法执行时间与字符集无关，是串长度的线性时间，即为O（n）。算法输出格式为（p；u；i，j，?），其中p和u与定义1中参数相同；?为满足条件?≥?min的模式频率，i…j是SA中一个区间。这种输出表示法的缺点是重复子串的出现位置不是顺序的，需要先利用SA数组将（p；u；i，j，?）转换为（p；u；SA[i]，SA[i 1]，…SA[j]，?）的格式，并采用排序法对所有最大重复模式子串的出现位置（即SA[i]，SA[i 1]，…SA[j]）进行排序，因此会增加算法的时间复杂度。所以我们在步骤1（3）中设计检测算法Check，既避免对重复子串进行排序，同时也可检测其是否满足限制条件。
3）判断步骤1（2）中识别出的重复模式是否满足限制条件：至少在每个串中出现mmin次且至少在q个串中出现，如表1中算法所示。
数组count=count[1..r]对位于各个Si中的重复子串进行计数。通过使用这两个数组，可用O（（j-i）log2r）时间复杂度的算法计算出R中重复子串是否满足满足限制条件：至少在每个串中出现mmin次且至少在q个串中出现。值得注意的是，算法首先检测重复模式出现频率?（即j?i 1），如果小于mminq，则没有必要进行进一步检测。函数BinarySearch返回索引t，表明重复子串位置SA[h]出现在串St中。设最大重复模式个数为O（α），对每个重复模式，数组count需清零r次，因而时间复杂度为O（αr）。对S中每个重复子串，二分查找最坏时间为O（log2r），设S中识别出的重复子串总数为φ个，则算法Check的时间复杂度为O（αr φlog2r）。算法Check还可进一步简化，如设置一个初始值为零的链表L，将二分查找返回值t加入进去；在最后一条语句output（R）后执行另一个循环，将L中t值清零，并执行count[t]=0；则算法时间复杂度可降为O（φlog2r）。
4）对S中所有的重复模式提取出第二字段即子串部分，作为公共模体的基元，并进行编号；编号可以依据输出顺序，也可以随机编号。
步骤2：建立一个新的串集合T={T1，T2，…，Tr}，其中Ti （1≤i≤r）的值为步骤1（4）中取得的值，即Ti[j]=Si中重复模式在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：设Ri为链Ci[j]的集合，其中[?j∈1..n-m（d k） d]，并对所有的i求出Ri。
步骤5：对集合Ri（为链Ci[j]集合的集合）建立Motif树（即另一种形式的广义后缀树），深度为m，叶子节点为Ci[j]所在的Ri。
步骤6：含有公共Ri的叶子节点的路径，即为公共模体。
3.3 算法举例
下面根据条件k=2， m=3， d=1，以图2中S1，S2，S3为例，根据步骤1-6识别S1，S2，S3中公共模体。
3.4 算法的性能分析
现在分析算法的复杂度。第一步骤已在前面分析，为O（rn αr φlog2r），经过改进可降为O（rn φlog2r）；步骤2的时间复杂度为O（rn）；在步骤3中，创建了长度为m（d k）-d的n-m（d k） d个链，因此步骤3的最差时间复杂度为O（rn），假设m，k和d为常数。步骤5中建立Motif树最差时间复杂度O（rn），步骤6也是线性时间。因而总的时间复杂度为O（rn φlog2r）。
3.5 与现有算法的比较与改进
本文研究工作与文献[13]的研究工作（问题1的算法）最为接近，其区别在于：
1）由于采用了基于后缀数组或QSA数组的重复模式识别算法挖掘串中最大重复模式作为基元，空间效率得到很大提高，与建立后缀树需20rn～40rn相比，采用基于QSA数组的算法空间复杂度为6rn，基于后缀数组算法空间复杂度为9rn。同时也避免了在广义后缀树上搜索得到k-子树的过程。
2）由于采用步骤1（4），可判断基元是否在所有序列中出现，一方面通过剪枝，在后续步骤中（步骤4和5）避免了大量的匹配操作，可有效限制搜索空间，极大地提高了实际运行时空效率。另一方面，也可回答下列问题，即“给定的一组序列，可能的motif仅在部分序列中出现，如何解决？”，此时可取q
4 结束语
本文提出了一种公共模体识别的有效算法，算法与基于后缀树或Trie树的同类算法相比在时间和空间效率上都得到了提高。
参考文献：
[1] Kellis， M.， Patterson， N.， Endrizzi， M.， Birren， B.