DNA序列中基于适应性后缀树的重复体识别算法.docx

DNA序列中基于适应性后缀树的重复体识别算法.docx

《DNA序列中基于适应性后缀树的重复体识别算法.docx》由会员分享，可在线阅读，更多相关《DNA序列中基于适应性后缀树的重复体识别算法.docx（13页珍藏版）》请在冰豆网上搜索。

DNA序列中基于适应性后缀树的重复体识别算法

DNA序列中基于适应性后缀树的重复体识别算法

霍红卫王小武

（西安电子科技大学计算机学院西安710071）

Email:

hwhuo@

摘要现有的在DNA序列中识别重复体的算法多数是基于比对的，对识别速度和吞吐量有很大的限制。

针对这个问题本文根据一个平衡重复体的长度和频率的定义，提出了一种基于Ukkonen后缀树的快速识别重复体的RepSeeker算法。

算法采用最低限制频率，最大程度地扩展了重复体的长度，同时为了进一步地提高RepSeeker算法的效率，对Ukkonen的后缀树构造算法进行了适应性改进，在构造时加入RepSeeker算法所需的结点信息并将叶子结点和分支结点加以区分，从而使得RepSeeker算法能通过直接读取结点信息来求得子串频率和子串位置。

这种改进较大地提高了RepSeeker算法的性能，而且空间开销不大。

实验中使用了NCBI中的9条典型DNA序列作为测试数据，并对后缀树改进前后的重复体识别算法做了比较分析。

结果表明RepSeeker在没有损失精度的情况下缩短了算法的运行时间。

实验结果与理论上的分析一致。

关键词重复体识别；适应性后缀树；Ukkonen算法；RepSeeker算法

中图分类号：

TP18

1.引言

基因组中含有许多重复元素。

例如，在人类基因组的约3.2109个碱基对中超过50%已被识别为各种重复元素[1,2]。

重复体识别对于分析新的基因组非常重要[3]，这是因为

（1）重复体以各种方式引导着基因组的进化过程；

（2）在进行同源查找之前需要对重复体进行掩模。

而在实际中有各种各样的重复体，大多数重复体的功能并未完全被理解和定义[4]。

当前研究表明某些重复体在基因表达和转录调控方面起着重要作用[5]。

与重复体结合的碱基可能导致基因重组，使基因组发生重大变化。

重复体种类繁多，它们包含着几个到数百的碱基对，有些可达上万碱基对。

一些人类的遗传疾病诸如脆性X染色体综合症、亨延顿氏症以及弗里德共济失调都与重复体长度的不规则性有关[6]。

一个重要的生物信息学问题是如何快速识别并有效的表示基因组中的重复体。

目前，解决重复体识别的方法大致有两类：

RepeatMasker[7]根据一个已经注释的重复体数据库对已知重复体进行查找。

这个数据库很大程度上依赖于同源序列的相似性。

这种方法不能用于处理新的基因组序列，是因为它不能为新测序的基因组构建所需的库信息。

而且，对于新的基因组，它的重复体库需要手工编撰，因为这种方法是面向特定基因组的。

对RepeatMasker数据库的从头识别仍然是生物信息学中的一个挑战问题[8]。

REPuter[9]是另一种方法，它摘取具有最大长度的所有重复体的相似对，并把重复体定义为一组具有最大长度的相似字符串对。

这两种方法都没有考虑重复体的出现次数。

一般而言，在真实生物序列中，重复体会出现多次。

例如，在人类基因组中Alu出现106次。

在复杂的基因组中，转座子一般出现几十万次。

因此，在识别方法中结合重复体的频率更合理。

具有生物意义的重复体的定义必须考虑重复体的长度和频率。

一些研究表明准确地定义重复体是不容易的。

一些方法只能找到短的重复序列或串联重复序列。

它们难以找出长且散布的重复序列。

最近的一些方法集中在识别重复体的边界上。

Price等人提出了RepeatScout算法[10]，该方法使用高频l-mer种子来查找重复体的边界，且用贪心法扩展每个种子，使之成为更长的同源序列。

Edgar和Myers研制了RILER软件包[11]，通过刻画重复体特征和局部比对来识别具有可靠边界的重复体。

意识到重复体出现频率的重要性，已有几种方法把长度和频率结合在重复体的定义中。

Zheng和Lonardi[12]给出了一种基于后缀树来查找DNA序列中重复体的方法，且时间复杂度为O（n2f）[13]。

Zheng和Lonardi算法的效率不高，是因为对于数十万碱基对长度的DNA序列，算法仍然难以有效工作。

尽管在试图定义和识别一个序列中重复体已有大量成果，重复体查找仍然是一个挑战性的问题。

文献[8]使用了局部序列比对策略和A-Bruijn图来解重复体查找问题。

然而，基于A-Bruijn图的方法的分析非常复杂和困难。

此外，这种方法需要对输入序列进行双序列局部比对。

其性能主要依赖于局部比对的结果。

Gu等人[14]提出了一种基于精确字统计法来估计大型真核基因组中重复结构出现频率的方法，该方法对比期望出现次数多的寡核苷酸进行分类。

本文中按照文献[12]中对重复体的定义，对Ukkonen后缀树构造算法做了适应性的改进，提出一种快速识别重复体的算法RepSeeker。

RepSeeker算法使用最低频率限制，并扩展重复体的长度使之达到最大。

在后缀树的构造过程中，对叶子节点进行编号，并把叶子节点的信息加入到分支节点中，叶子节点和分支节点所包含的信息不同，以使RepSeeker算法能够直接从节点中获得子串的频率和位置信息。

这种改进大大提高了RepSeeker算法的性能，而且空间开销不大。

RepSeeker算法使用来自NCBI中的9条序列进行了性能测试。

并对改进前后算法的性能作了比较。

实验结果表明，对其数据结构所做的改进大大降低了RepSeeker算法的运行时间，同时又保证了识别的精度。

实验结果与理论上的分析一致。

2.RepSeeker算法

2.1表示、约定及基本定义

表示一个有限非空字母表。

字符串和字母用斜体字体。

|S|表示字符串S的长度。

S[i]表示字符串S中的第i个字符，1i|S|。

S[i,j]表示S的子串S[i]S[i+1]..S[j],1ij|S|。

约定S[i,i]=S[i]。

L-串（或L-子串）是一个长为L的字符串（子串）。

Ai表示A在S中的第i个出现。

Ai表示A的第i个拷贝。

有时我们会交替使用这些表示，因为根据上下文就可以明确它们是表示子串还是表示子串所在位置。

如果字符串A是字符串S的一个子串，且在S中出现多次，则称A是S的一个重复体。

在分子生物学中，重复体就是一个碱基序列在染色体中多次出现的拷贝。

阈值是指一个重复体在一个序列中重复出现次数所指定的最小值。

滑动窗口（见图1）是某个定长的可视框。

2.2定义

定义1设A是S的一个重复体，（A1,A2,…,Am）是A在S中出现的一个有序表，Ai是A在S中的第i个出现，m是A在S中的出现次数，且m2。

设B是A的一个子串，（B1,B2,…,Bk）是B在S中出现的一个有序表，k是B在S中的出现次数。

如果B在A中从位置s开始，那么B=A[s,s-1+|B|]，0s|A|-|B|。

如果k=m且每个Bi在Ai中出现的偏移量相同，则称B是A的一个子重复，i=1,2,..,m。

定义2A是S的一个非平凡子串，当且仅当A是S的一个非空、真L-子串，0<|A|<|S|，L是重复体的最小长度。

定义3如果A是S的一个具有最大长度的非平凡子串，且A在S中至少出现Fm次，A的每个非平凡子串是A的一个子重复，其中Fm是指定的重复体最小出现频率，则称A是S的一个基本重复体。

性质1如果A是S的一个基本重复体，且出现频率为f，那么A的每个L-子串出现频率均为f。

证明略.结论可由定义直接而得。

■

上述性质使我们在计算重复体的出现频率时，可以删除大多数非候选的重复体。

在计算中，使用后缀树作为基本数据结构，来统计具有给定阈值的其L-子串的频率。

2.3算法描述

RepSeeker算法首先找出输入序列S中的所有基本重复体，然后输出重复体的一个有序表。

重复体表中的元素是一个数对，表示一个重复体在输入序列中的起始位置和结束位置。

因此，识别基本重复体的问题可以转换为寻找重复体的边界问题。

因而，RepSeeker算法检查输入序列中每个位置，并确定一个位置是否是重复体的一个边界。

使用穷尽算法查找基本重复体是不切实际的，因为在S中有O（n2）个子串，对于每个长为m的子串，需要检查O（m2）个子串。

因而，我们构造输入序列S的一棵后缀树，帮助进行频率统计。

按照性质1，可得：

如果A是S的一个基本重复体，那么A的所有L-子串出现频率相同，且至少为Fm。

于是，RepSeeker算法计算出所有L-子串的出现频率，并根据频率数组把频率相等且至少为Fm的放在一块中。

由于这只是一个必要条件，RepSeeker算法会检查这些块，并分裂那些包含非子重复的块。

进而，算法对所得重复体进行扩展，最终进行分类。

RepSeeker算法由以下5步组成。

第1步：

计算子串在滑动窗口中的频率

我们把重复体的最小长度L作为滑动窗口的宽度，并计算长度为W的子串在此窗口中出现的频率。

滑动窗口每次向右移动一个位置。

设频率数组为f，f[i]表示在位置i开始的L-子串的出现次数，即S[i,i+L-1]的频率。

图1中显示了输入序列S的频率数组f的值。

在图1中，f[0]=3表示在位置0开始的长为4的子串S[0,3]在S中出现3次，S[0,3]=“ABCD”；f[5]=3表示在位置5开始的长为4的子串S[5,8]在S中出现3次。

图1.频率数组f（L=4,Fm=3）

第2步：

求出频率相等的块

根据第1步中计算的频率数组，可以计算出出现次数至少为Fm的L-子串的起始位置和结束位置。

在RepSeeker算法中，使用l和r分别记录它的左右边界（即起始位置和结束位置）。

如果分别来自l数组和r数组的两个元素在输入序列S中的位置相同且位置为i，那么它们表示同频率的候选重复块，（l[i],r[i]）。

基于频率数组f，我们把序列划分成同频率的块。

特别是，对于任何位置i，如果f[i]≠f[i-1]，那么位置i是这个同频率块的起始位置；如果f[i]≠f[i+1]，说明同频率块的最后L-子串在位置i开始，也就是说，i+L-1是这块的结束位置。

每当得到一个起始位置或结束位置时，就把这个位置插入到l数组或r数组中。

完成l数组或r数组的计算之后，接下来是对l数组和r数组中的元素进行配对。

（l[i],r[i]）是我们所找到的第i个相同频率块。

例如，对于l[0]=0且r[0]=9，图1中的块S[0,9]=“ABCDEEBCAD”是一个频率相等块，它的长为4的7个子串出现频率一样。

对于图1中的示例，去掉频率小于Fm（=3）的子串后，最终频率相等块为S[0,9],S[12,21],S[24,28],S[29,33]和S[36,41]。

第3步：

子重复检查

对从第2步得到的块进行检查，设某一块Di=S[i1,i2]，块长度为length,出现频率为k，其块内所有长度为L的子串的频率为m，若k

从左到右依次求出块Di中所有长度为L的子串在S中的出现位置，然后对开始位置相邻的两个子串的出现位置序列进行比较，若存在对应次序上的位置不相邻，则分裂块Di。

例如，图1中，子串S[0,9]的发生频率为2，但在第2步中被识别为重复体，通过子重复检查，发现S[1,4]=“BCDE”的出现位置为{1,13,30}，S[2,5]=“CDEE”的出现位置为{2,14,36}，两子串的第3次出现位置不相邻，因此将2插入l数组，将5插入r数组。

对数组l和r的插入要求进行有序插入，即插入后数组仍保持有序。

所以S[0,9]可分裂成S[0,4]，S[2,7]，S[5,9]三个重复体。

第4步：

重复体扩展

为了尽可能得到更长的重复体，我们归并有重叠的重复体，如果归并后的重复体频率至少为Fm。

具有较高频率的重复体仍然被保留下来，较低频率的重复体被扩展。

令重复体A和B含有重叠块，重叠块称为C。

归并满足的条件如下：

merop（A,B）=

100%OP

且

frequ