基于块排序索引的生物序列局部比对查询技术参考模板

1、基于块排序索引的生物序列局部比对查询技术*Block Sorting Index-based Techniques for Local Alignment Searches on Biological Sequences李永光 王 镝 王国仁 马宜菲 (东北大学信息科学与工程学院计算机系统研究所 沈阳 110004)Abstract A common query against large protein and gene sequence data sets is to locate targets that are similar to an input query sequence. T

2、he current set popular search tools, such as BLAST, employs heuristics to improve the speed of such searches. However, such heuristics can sometimes miss targets, which in many cases is undesirable. The alternative to BLAST is to use an accurate algorithm, Such as Smith-Waterman(S-W) algorithm. Howe

3、ver, these accurate algorithms are computationally very expensive. Recently, a new technique, OASIS, has been proposed to improve the efficiency and accuracy by employing dynamical programming during traversing suffix tree and its speed is comparable to BLAST. But its main drawback is too much memor

4、y consuming. We propose an efficient and accurate algorithm for locally aligning genome sequences. We construct a block sorting index structure for the large sequence. The index structure is less than the suffix tree index and can be fit for large data size. Experimental results show that our algori

5、thm has better performance than OASIS.Keywords sequence; block sorting index; accuracy1 / 141. 引言随着人类基因组计划的不断发展，在结构基因学和功能基因学的研究过程中，生物序列的相似性分析成为一种有效的手段1。通常情况下，两条DNA长序列，可能只在很小的区域内(密码区)存在关系；不同家族的蛋白质往往具有功能和结构上的相同的一些区域。因而，研究序列局部相似性比研究全序列相似性往往更有意义。序列局部比对是一种关于片段相似性的定性描述1。通常的方法是通过动态规划2进行精确查找两个序列的最优局部比对，但因其代

6、价太大，对超长生物序列直接应用这种技术是不可行的。为了提高查询速度，启发式算法当前被广泛应用，这些算法可以分为三类：基于哈希表的算法、基于频率空间过滤的算法、基于后缀树索引的算法。基于哈希表的典型算法有FASTA5和BLAST6。FASTA的核心思想是对数据库序列中的所有模式进行哈希索引。在查询时基于哈希索引可以快速检索出可能的匹配模式。然后根据这些模式的扩展得到更长的序列。BLAST算法是建立在严格的统计学的基础之上，它主要用于发现具有较高的相似性的局部比对。在大多数情况下，根据局部比对参数会产生若干个HSP(High-score Sequence Pairs)。然后把所有匹配分值大于阈值的

7、HSP作为结果。这种基于启发式算法的过滤原理尽管在一定程度上尽量减少丢失局部比对的结果，但精确率却不可能达到100%，一些符合条件的匹配序列可能不在查询结果中。MRS4是基于频率空间过滤的代表方法，它采用滑动窗口和小波变换技术，设计了一种高效的生物序列搜索方法。根据频率距离是编辑距离的上界的过滤原理，可以节省大约50%的编辑距离计算。这种技术比较适合长的DNA序列，但不适合蛋白质序列的比对分析。另外一些方法是基于后缀树技术的, 该类方法一般是通过对一条序列构建后缀树，然后利用后缀树快速查找到匹配的模式。对于后缀树中的某个分支，如果完全匹配的模式的个数超过给定阈值，那么就对该分支使用BLAST作

8、进一步的查询。这类算法对阈值比较低的查询处理效率太低。OASIS3通过在遍历后缀树的同时运用动态规划算法，提供了一种精确的序列局部比对的搜索方法。其实验结果表明，它比动态规划算法快很多，而且可以和BLAST相提并论，但存储后缀树需要很大的存储空间，其构造过程也是复杂的。本文首先提出了一种基于块排序索引算法来计算序列比对，我们引入基于这种索引结构的片断向量的概念，通过片断向量的扩展完成动态规划算法，从而完成序列比对查找过程。算法采用A*8搜索策略，进行过滤，当计算查询序列部分匹配时，搜索算法同时计算匹配剩余查询所得分值的上界，根据阈值，如果上界分值小于阈值，那么该片段向量就被过滤掉，从而减少搜索

9、空间的计算。在搜索的每一步，算法首先扩展可能产生最优局部比对的片断向量，这就能保证返回的结果按照比对的分值由达到小排序。实验结果表明，我们的算法在性能和过滤能力上都优于OASIS。本文的主要贡献在于通过运用较小的索引结构，我们提出一种精确有效的局部比对搜索算法，因而这种算法为在大的生物数据集提供了一种可行有效的精确局部比对查找策略。本文余下部分的组织结构如下：第2部分给出了本文的背景知识；第3部分具体介绍了核心数据结构和算法；第4部分给出了测试结果和性能分析；第5部分总结全文。2. 背景知识动态规划算法能精确计算生物序列全局和局部比对，在这一部分，我们首先介绍生物序列的局部比对的概念，然后讨论

10、动态规划算法的相关工作，最后介绍块排序索引结构。2. 1 序列局部比对给定两条字符序列Q=q1q2qm, T=t1t2tn,这两条序列的局部比对是指通过某种方式是“排列”Q和T的子序列。图1给出一个局部比对例子，同时给出三种类型的局部比对操作。l 替换：用相同或者不同字符进行替换（如图1中的1和2标记）。l 删除：允许在目标序列中忽略一个字符（如图1中的3标记）。l 插入：允许在查询序列中忽略一个字符（如图1中的4标记）。31124Target：C A B I N 1Query： D R A I N 图1. 局部比对的样例序列是通过三种操作所得到的分值之和来确定比对得分。每种操作可以形象描述为

11、替换操作“”，其中插入操作表示为“ ”,删除操作表示为“ ”，操作的罚分可以标记为S（）。可以通过替换矩阵S存储这些分值，所以“”的操作代价可以简单表示为S,。表1给出一个通用的替换矩阵的样例。称为单位编辑距离矩阵（因为完全匹配分值为1，否则为-1）。ACGTA1-1-1-1-1C-11-1-1-1G-1-11-1-1T-1-1-11-1-1-1-1-1-表1：单位编辑距离矩阵2. 2 Smith-Waterman(S-W)算法Smith-Waterman算法通过动态规划2的方法计算两条序列局部比对的最大的可能分值，计算过程需要O(mn)空间和时间的复杂度。算法形成一个m×n的矩阵G

12、，Gi,j保存查询序列和目标序列分别结束在qi和tj的最大比对分值。每一个Gi,j是通过下面的公式(1)计算出来的： （1）2. 3块排序索引结构对于长度为n的序列S（我们的应用中，S以非序列中出现的字符作为结束字符，我们选用结束字符$，并规定$的字典序小于序列中的所有其它字符），进行块排序时对序列S循环移动0n-1位，得到一个n行列表，将各行按升序排序，就得到了矩阵M，如图2所示，左侧是对序列AGTACGCCTAG$的块排序结果, 这里我们不详细说明构造的过程，细节可看见7等文献。图2：序列AGTACGCCTAG的块排序索引结构下面我们描述基于块排序的索引结构。为了描述的方便我们将矩阵M的第

13、一列称为F列，最后一列成为L列。Fi和Lj分别表示F列中的第i个字符和L列中的第j个字符。另外，我们引入函数count (F,i)在F列的前i个字符中，Fi出现的次数。类似的，我们也定义count (L,j)。现在，我们给出在块排序结构中一个字符的后继字符的表达。对于一个字符Fi，如果Fi=Lj并且count（F，i）=count（L，j），则Fj是Fi在给定序列中的后继字符。块排序结构由F列和T列构成，Ti是一个整数，它表示Fi的后继字符在F列中的位置，即FTi是Fi在原序列中的后继字符。例如，图2给出序列AGTACCGCCTAG$的块排序结果。F3的后继字符是FT3，即F5。块排序索引就有

14、如下的保序性：性质1对于F中的任意字符Fi和Fj，如果Fi=Fj并且i<j，则Ti<Tj。证明：令si是矩阵M中的第i行代表的序列，则Fi是si的第一个字符。对于两个序列si和sj，它们的先后关系是由两序列中的第一个不相同的字符决定的。如果Fi=Fj，则si和sj的顺序与sTi和sTj的顺序是相同的。因此，若i<j，则Ti<Tj。例如，在图2中 F1=F2=F3=A，则（T1=3）< (T2=5) < (T3=7)。在这边文章中，我们对DNA序列建立索引，因此序列的字符集为A，C，G，T，$。令s是一长度为n的DNA序列，l1，l2，l3，l4分别代表A，C

15、，G，T在F列中的首个字符的位置。则索引中的F列具有如下性质： 因此，根据块排序索引的性质，我们仅需要存储5个整数和T列作为索引就可以完成任意子序列的查找过程，因而基于该索引，我们能完成局部比对的查找过程。3. 基于块排序索引（Block Sorting）局部比对查找的算法根据块排序索引结构特点，本部分首先给出有关定义，然后根据问题的定义，给出算法的具体实现过程，最后给出算法的削减策略。3. 1 定义定义1 (片断向量) 根据索引结构中的T列，我们会得到一组向量，其中每一个向量保存某个字符在索引中的位置信息，该向量称为片断向量。例如，在图2中 10，11 代表字符T在索引中的第10和11位置上

16、，1 ，2代表在索引第一次扩展后，字符T的后继字符在索引中的位置为1和2。定义2（根向量）根据块排序索引，我们首先会得到每个字符在索引中的范围，根据这个范围我们会得到每个字符初始状态的位置向量，该向量成为根向量例如，在图2中，1，2，3 为字符A的根向量，4，5，6 为字符C的根向量。定义3（后继向量和前驱向量）根据索引结构中的T列，片断向量frag1的后继产生的向量称为frag1的后继向量，frag1称为前驱向量。例如，在图2中，1，2 为 10，11 的后继向量，3 和5 是1，2的后继向量，同时 10，11 为 1，2 的前驱向量，1，2 为 3 和 5 的前驱向量。定义4（向量扩展）假

17、设frag1为一片断向量，我们得到frag1所有后继向量的过程称为向量frag1的扩展。例如，在图2中，我们开始得到字符T的根向量 10，11，然后我们根据索引结构中的T列信息就会得到其后继向量 1，2，依次类推，我们通过扩展就能得到所有以字符T开始的所有子序列。图3：向量10，11扩展过程定义5（扩展路径）从根向量开始，片段向量扩展过程会产生一些路径，其中每一条路径对应于目标序列中的某一子序列，该子序列的长度即为路径的长度，该路径称为扩展路径。例如，在图3中，10158和11270称为索引扩展路径，分别代表子序列TACG和TAC$。3.2算法描述算法核心为在局部比对的搜索过程中执行优者为先A

18、*算法。在块排序索引结构扩展过程中完成动态规划算法，描述如算法1所示，其中输入参数如下：l BS：目标序列的块排序（Block Sorting）索引结构；l Q：查询序列，q1q2qm；l S：任意的一个替换矩阵；l minScore：比对要求的最小值；这些参数同时也是初始化函数的输入参数，所有分值大于minScore的比对按分值由大到小返回。算法1：局部比对搜索过程（）Search Process ( )输入：目标序列的块排序索引BS；查询序列Q；打分矩阵S；阈值minScore输出：目标寻列中与查询序列相似的子序列；算法描述：1: H,PQ= Initialize(BS,Q,S,minSc

19、ore);2: while ( !Empty(PQ ) ) do3: SearchVector PQ.Pop();4: if ( SearchVector.tag = Viable ) then5: for ( successor Successors(SearchVector) do6: newBlock Expand(SearchVector,successor,H,Q,S,minScore);7: if ( newBlock.tag = Accepted newBlock.tag=Viable);8: PQ.Push(newBlock);9: end if10: end for11: e

20、lse if ( SearchVector.tag=Accepted) then12: Return subsequence containing Path ( SearchVector.vec);13: end if14: end while在算法中，每个SearchVector与块排序索引扩展过程中的片断向量一一对应，如图4表示图2索引结构一步扩展后，得到片段向量1,2对应的SearchVector。SearchVector表示部分目标序列与查询序列的比对，即查询序列和目标序列结束在对应于其片断向量的路径对应的子序列的比对。每个SearchVector存储查询序列从开始到每个字符比对得分，

21、同时计算计算匹配查询序列剩余部分的最大的可能比对得分。这些得分的总和用来在优先级队列PQ中组织这些SearchVector。算法总是选择优先级队列中的首个元素然后扩展它的后继向量。因为算法总是扩展在队列首部的SearchVector，所以可以近似看成优者为先的A*算法。算法中SearchVector为核心数据结构，现简要说明一下它包含以下数据域：Vec：扩展过程中的片断向量，其扩展路径与目标序列中子序列一一对应。Z：为某一向量 Z0,Z1.,Zm，其中Zi 记录对应与从根向量扩展到片段向量Vec的扩展路径所对应序列与查询序列Q结束在qi的子序列之间的最好局部比对得分。如果当前的比对被削减（削减

22、策略将会后面讨论），则Zi被设为-。向量Z，粗略的讲，是和表2中S-W算法矩阵中的某一列对应的。maxScore：沿当前路径所得到的最大比对得分。这个分值是是对应于当前路径的序列与查询序列的任何子序列的最好的比对得分。maxf：通过进一步扩展该Vec所能得到的最大的可能比对得分。tag：标记当前SearchVector的状态，可能的取值为Accepted、Viable、Unviable，如果查询序列与当前路径的序列比对最好得分已经得到，并且达到阈值的要求，则该SearchVector被标记为Accepted，当该SearchVector在优先队列PQ的首部时，我们就会得到该查询的最好的局部比对

23、结果。如果tag标记为Viable，则说明在当前路径上还可能得到更好的比对得分。如果当前路径无论怎么扩展所得到的最好局部比对得分都比都达不到阈值的要求，则tag标记为Unviable，对应的片断向量也就被削减掉。标记为Accepted和Viable的SearchVector被放入优先级队列。优先队列是以maxf的值进行组织的，所以首先扩展的SearchVector必须保证队比队列中的其它SearchVector可能产生更高比对得分。通过优先队列的排序，算法扩展队列首部的SearchVector的后继向量，如算法3所示。扩展过程会增加队列中的元素，搜索结束的条件是队列为空（具体看算法1）。图4：

24、对应于片断向量1，2的SearchVector结构3.2.1初始化过程算法开始计算一个启发式向量H，同时根据索引结构中的根向量初始化优先队列。具体步骤如算法2所示。H中的每个元素hi代表qi+1qi+2qm与任意目标最大的可能比对得分（假设Q=q1q2qm）。因此，计算H是比较简单的。假设插入和删除操作的罚分都为非正值，hm可设为0，因为查询序列左部分时空字符串，然后我们就能逐步计算H中其它元素的值：hi-1=hi + 替换qi-1的最大得分。每个RootVector对应于索引结构的根向量，其路径长度为1。因为比对分值保存在Z向量中，如果分值设为-则这部分比对可以过滤掉。被过滤掉的意思即当H中

25、的hi比给定的阈值小，这意味继续计算也不可能得到一个更有意义的序列局部比对。算法2：初始化（H , PQ）Initialize( )输入：目标序列的块排序索引SA；查询序列Q；打分矩阵S；阈值minScore；输出：启发式向量H；优先队列PQ；算法描述：1: hi hi+1 + max S(qi );2: for ( every RootVector ) do3. compute Zi , maxScore, maxf, tag;3: if ( RootVector.tag=Accepted RootVector.tag= Viable ) then4: PQ.Push (RootVector

26、);5: end if6: end for7: return H,PQ;3.2.2扩展过程算法3给出了完整的块排序索引扩展算法。输入当前的片断向量和其后继片断向量，算法开始时初始化Successor（为一SearchVector结构）的Vec和maxScore。其中矩阵M的各个元素根据S-W算法得到的，稍微不同的是，每个元素不会设为0。在S-W算法中,0代表重新进行一个新的比对。在我们的算法中不需要这种设定，因为从索引中的根向量开始，所有目标序列的子序列都会找到，继续设为0意味着在别的地方重复以前的工作。因为与查询所有的比对都会在计算M时考虑到，并且根据块排序的重要性质，每条子序列都对应于索引

27、扩展中的每一条路径，所以所有从查询序列开始的比对都会考虑，因而不会丢失任何比对的结果。算法3：基于块排序过滤的比对扩展算法 Expand ()输入：当前的片断向量；后继片断向量；输出：可能产生结果的后继向量算法描述：1: Successor.Vec the first successor SearchVector;2: Successor.maxScore SearchVector.maxScore;3: for ( i 1.m ) do4: Mi,0 the current SearchVector.Zi; 5: Mi,1 max( Mi-1,0 + S( qi c), Mi-1,1 + S

28、( qi -), Mi,0 + S( - c)6: if ( Mi,1 > Successor.maxScore ) then7: Successor.maxScore Mi,1;8: Successor.maxf max (Mi,1 + hi );9: end if10: end for11: if (Successor.maxScore >= Successor.maxf) && (Successor.maxScore> minScore ) ) then 12: Successor. tag Accepted;13: return Successor;

29、14: else if (Successor.maxScore < minScore) 15: Successor. tag Unviable;16: return Successor;17: end if18: if ( IsEnd (Successor) then19: if (minScore Successor.maxScore) then 20: Successor.tag Accepted;21: Successor.maxf Successor.maxScore;22: else 23: Successor.tag Viable;24: Successor.Z Mi,1;2

30、5: end if26: return Successor;3.3削减策略计算mi,j后，算法进行比对计算削减，在不能继续扩展的比对位置上设置比对得分-。这个过程允许算法丢弃一些SearchVector，因为这些SearchVector或者不能继续扩展，或者被其它扩展路径包含。在下列三种情况下执行削减策略：1. 比对得分为非正数（mi,j 0）：为了避免重复工作，所有比对得分为负数的片断向量就被消减掉。考虑一部分查询qaqa+1qb和一部分目标序列tctc+1td进行局部比对，如果开始位置位qatc，结束位置为qbtd的最大比对得分为负数，那么我们就能保证任何关于qaqb和tctd的比对得分会

31、比qb+1和td+1比对得分小，而第二个比对会在另外的扩展路径得到。所以没有必要保存第一个比对，更没必要扩展第一个比对结果。2. 存在的比对已经是最优的（mi，j+hiblock.maxScore）：基于良好的启发性，在当前的扩展路径上，对比对进行扩展不会得到比当前结果更好的结果，即当前片断向量的后继向量不会产生更好的比对。3. 阈值限制（mi，j+hi minScore）：当前比对即使继续扩展也不会产生一个能大于阈值的比对得分。在这些条件的限制下，算法就能在某条扩展路径上停止继续扩展，通过启发式向量H，我们能够确定作进一步扩展的最大可能比对得分的上界，即maxf的数值。如果当前的扩展路径的比

32、对得分不小于这个商界（maxf maxScore），算法就停止继续扩展，然后根据是否达到minScore的限制，返回一个Accepted或者Unviable的结果。当比对得分还能进一步提高（maxf > maxScore）,则要看有无必要在该路径上作进一步扩展：如果maxf小于minScore，则返回标记为“Unviable”的SearchVector，否则标记为“Viable”。 4. 性能测试和分析为了测试基于块排序过滤器的局部比对搜索算法的性能，我们在同样的平台上实现了我们的算法和基于后缀树的算法（OASIS），所有的测试都是在一台Intel Pentium 4 *2.0GHz C

33、PU，2G内存，1T硬盘IBM x Series 345服务器上进行的。我们使用比较流行的人类染色体chr22 和chr18的片断和Drosophila melanogeneses作为测试数据，使用当前比较流行的单位编辑距离替代矩阵进行比对分值计算。系统使用的操作系统为Windows 2003，使用Microsoft Visual Studio.2003编译器编译。4.1阈值minScore对算法性能的影响OASIS 是基于后缀树技术进行精确局部比对查询的算法，在实验中，我们算法与OASIS性能比较，我们通过minScore 参数控制查询的敏感度。公式（2）给出期望得到的比对个数E与比对得分S

34、的关系： （2） 其中，m是查询序列的长度，n为数据库的大小，K和是通过BLAST计算的常数，所以算法中的minScore可以用公式（3）进行计算： （3） 图 5 给出根据minScore的影响下，两种算法的比较。对于查询序列，我们用Drosophila melanogenetic chromosome NT_033778 的子序列（TAAACTTGGCCTGCTT）,用Drosophila melanogenetic chr4 全序列作为目标序列。查询时间在图中是取以2为底数的对数数值，查询序列的长度为16，minScore的范围412，因为我们算法的启发性比较好，所以从图中不难看出，mi

35、nScore对我们的算法的影响比较小。图5：阈值minScore对查询性能的影响4.2 查询长度和数据库数据大小对算法性能的影响我们根据查询序列的长度和目标序列的大小进行性能比较。图6说明了不同查询长度对查询性能的影响，图7给出不同目标序列大小对查询性能的影响。我们使用8M的目标序列，目标序列是人类的基因序列chr22片段，查询序列为随机截取chr18的子序列。因为在我们实验机器的内存限制下，OASIS无图6：查询长度对查询性能的影响法为这些数据集大于10M的目标序列构建后缀树索引，而Block Sorting (BS)我们可以为17M的数据构造索引结构。所以在图6中我们的目标序列大小为8M的

36、数据，图7中OASIS算法我们最大的目标序列为9M的数据，BS算法则测试到17M的数据。扩展对于OASIS算法对于返回的Accepted 结果，如果为未非叶子节点，则必须通过遍历后缀树后，然后根据得到的叶子节点才能得到子序列在目标序列中的位置，而基于Block Sorting（BS）算法中，返回Accepted结果的Vec中就已经保存了子序列的位置信息。所以在一定程度上提高了速度。图7： 目标长度对查询性能的影响4.3过滤效率的性能比较我们同时计算所要求扩展的列，即目标序列中需要进行动态规划计算的部分，从而反映出两种算法的过滤效率。图8给出两种算法的过滤能力比较，对于Block Sorting（BS）算法每次扩展一步后就会加以判断当前可能的最优扩展的片段向量，然后对此片段向量进行继续扩展，而对于OASIS，根据后缀树结构，需要根据节点的对应的子序列的长度进行多步扩展后才能对当前的所有可扩展后缀树节点做出最优的判断，而算法所有扩展步数之和即为动