生物信息学 序列拼接

《生物信息学 序列拼接》由会员分享，可在线阅读，更多相关《生物信息学 序列拼接（56页珍藏版）》请在金锄头文库上搜索。

1、基因组序列拼接段倩倩序列拼接序列拼接任务即将测序生成的reads短 片段拼接起来，恢复出原始的序列。该问 题是序列分析的最基本任务，是基因组研 究成功与失败的关键，拼接结果直接影响 到序列标注，基因预测、基因组比较等后 续任务。基因组序列的拼接也是基因组研究必须 解决的首要难题。其困难不仅来自它的海 量数据(以人类基因组序列为例，从数量为 10兆级的片断恢复出长度为亿级的原始序 列)，而且源于它含有高度重复的序列。拼接问题的难点DNA测序数据有其固有的四个的特点， 他们也正是解决实际的序列拼接问题的难点 所在：1测序有误差2不完全覆盖性3序列所在链不确定4重复序列的干扰1测序有误差由于测序技术

2、的局限，难免会出现测序 错误，尤其是在序列的末端，一般错误率 可控制在1以下。所以对每个碱基一般有 一个正确概率，以质量打分的形式给出。 因此每个ri都有个可信度。而read与read之 间有不同程度的重叠，由此导致有的重叠 可信度高，有的重叠可信度低。2不完全覆盖性不是所有的碱基被测序的次数都等于 平均测序覆盖度。极端的情况，可能会出 现源基因组序列上部分区域未被测序的情 况(这段区域称为gap)。即，测序的reads 集合不是原始基因组序列一个完整覆盖。 此时需要借助于各种图谱如：基因组指纹 图谱(genome fingerprint map)， 基因组级 物理图谱(genome-wide

3、 physical map)，细 胞发生图谱(cytogenetic maps)等协助对 reads进行定位3序列所在链不确定由于测序过程中无法确定特定片断属于DNA 双链中的哪一条链上，所以我们在拼接过程中并 不清楚使用的是read的正义链，还是其互补链。4重复序列的干扰DNA序列自身含有高度重复的子序列，它们 一种表现为短序列的串级重复，比如：(GGAA)n 。或AmTn等。另一种表现为大量相似序列(其拷贝 数可达几十万)散布在基因组的各个地方。 Repeat的存在，将导致fragments间overlap的不 真实性，进而产生错拼的结果。因此在拼接过程 中耍确定这些序列的形式及大小，才能

4、保证以高 概率恢复出其在原始真实序列中的位置拼接算法评价以上拼接问题的四个难点不仅极大的增 加了解决实际拼接问题的难度，而且从某种 程度上说无法完整地恢复出原始DNA序列来 。即实际上仅能构建出若干个contig(重建的 fragments的一种排列形式，它覆盖基因组 上一段连续区域)这些contig将指导测序项目 finishing阶段的实验方法最终构建DNA完整 序列。目前，国际上对拼接软件的公认评价 标准包括两方面，即重建出的contig的数目 和准确度。我们发展的基因组序列拼接新 算法的目标是在确保准确性的前提下，构 建尽量少的contig，以减少测序后期大量的 人力和财力的投入。基因

5、组序列拼接算法研究现状现在最常用的拼接程序使用的拼接算 法可分成两类，一类是将拼接问题转化为 在图中寻找的Hamilton路径的问题；另一 类是将拼接问题在某种特殊情况下转化成 寻求图中的Euler路径的问题。他们均有其 成功的典型算法。1.转化为Hamilton Path问题每个DNA片段(read)相当于图中一个结 点，如果两个片段之间存在着重叠(overlap) 关系，则在两个结点之间定义一条边，而沿 着DNA原始序列从头到尾，则必然经过每个 结点一次且仅一次，即是一条Hamilton路径 。一条contig表示图中一条简单路，此类算 法以Phrap，TIGR Assembler，CAP

6、3， GigAssemble等为代表。他们都是遵循“overlap-layout- consensus”的框架。首先，为了构建图。计 算任意两个read间可能的比对情况。其次， 通过去除歧义的或者不确信的边得到较为准 确的图，并在其上寻找非交叉的简单路的集 合，该集合对应于contig的集合。最终，通过 对包含在一个简单路上的所有read进行多序 列比对，为每一个contig构建一个一致性序列 (consensus sequence)。2.转化为Euler Path问题EULER是这类算法的代表。与传统方 法沿着“OverlapLayoutConsensus”路 线不同，它不计算各个read之

7、间的Overlap ，即没有Overlap步骤。它的大致想法如下：为了排除read中的错误，获得Error- Free的read，将所有的read切割成小片n- mers。将每个read和Gk的近似进行比对，寻 求read的最小改变能够使得read的所有n- mers包含在Gk的近似集合中。从而构建了 高质量序列，而对于Poor read，直接抛弃 ，对Chimeric read(两端在n-mers中但整体 不在的reads)进行特殊处理。初始的想法是要实现去除reads中的 测序错误的目的，如果知道原始序列G， 那么直接使用测序获得的read和G进行比 较即可。但是实际上G并不可知，那么退而

8、求 其次， G的序列片断Gk亦可，事实上Gk亦 不可知。所以将所有的read切割成小片n- mers，所有Solid的n-mers形成的集合称 为Gk的近似。最后，构造De Bruijn图。现有算法的主要问题虽然已经开发了以上的算法，基因组 序列拼接问题尚未彻底解决，以上两类算 法都存在着各自的缺陷。对于第一类算法来说，实际上是在图中寻找 一条使得评价函数值最优的Hamilton路径，这是 一个NP完全问题。一般都采用greedy-merging的算法近似求解 。由于这种step-by-step的局部贪心算法，其明 显的局部特性忽略了reads间“长距离”或者整体性 的联系，从而导致了拼接错误

9、，即拼接结果和真 实的DNA原始序列不同。最近研究指出，在对已 知序列的流行性感冒嗜血杆菌基因组的拼接过程 中，无论是Phrap，TIGR Assembler，还是 CAP3，都发生了拼接错误的现象。对于第二类算法来说，它只能在特殊 的情况下，才能将问题简化成寻求一条 Euler路径，最终的结果是从多条候选的 Euler超路中选择出来的。EULER算法依然 存在拼接错误，且结果选择的过程没有理 论依据。EULER软件在实际数据集的运行 速度上和第一类算法相当。更重要的是， EULER采用的算法过于独立，很难利用其 他辅助生物信息，导致其实用性和流行性 大打折扣。局部搜索(Local Searc

10、h)方法胡杰将局部搜索方法运用于一个具体的问题，需要对以下四项进行明确的定义。1将原同题表示成一个优化问题，即定义一个可行域以及在该可行域上的一个目标函数。2定义可行域中的邻域结构，即说明满足什么条件的两个点相邻。3确定在邻域中的搜索方式。4局部极值点的处理。如果当前解点邻域中 的所有点的目标函数值都比当前点大，那 么这点就称为局部极值点。在一些问题中 ，局部极值点就是全局最优点。而对另一 些问题而言，局部极值点已经满足求解实 际问题的需求。基于局部搜索的序列拼接算法框架主要目标是在layout阶段尽可能准确的 前提下，获得更长的contig。具体的，使用局部搜索算法求得数据集上近似全局的最

11、短超串；然后根据求得的最短公共超串对 应的fragment的排列关系为基础获得“consensus segment”1.Overlap定义如果一个串的前缀是另一个串的后缀则 认为这两个串之间存在overlap，并根据 over-lap构建超串。对给定的串f和g，存在 多个可能的overlap关系比如说，若f=ACTGGGAGCAGC， g=AGCAGCTTTTACT，那么他们之间至少存在两个overlap形式。在我们的算法中，仅考虑两个串之间最 大的overlap情况，并定义overlap(f，g)表 示f和g之间存在的多个overlap关系中最长 的一个overlap所包含的字符个数。在上面

12、的例子中overlap(f，g)=6。如果 f和g之间overlap区域长度小于M(M是一个 足够小的正整数)，则overlap(f，g)=0。2.优化目标定义我们对reads集合S中的每个元素按照 其左端在超串t中首次出现的位置进行排序 ，并沿超串t从左向右依次读入，并将其对 应为序S=sl，s2，Sn)。我们用P(S)表 示这个由字符串为元素构成的序。在序 P(S)中，对于每一个连续的字符串元素对 si， si+1，都存在overlap(si，Si+1)。因此，字符串的一个排列等价于一个超串t及作用在其上的函数overlap，超串t的长度length(t)=t=ssIs|- overlap

13、(si，Si+1),为了求超串t具有最小的长度。因为在给定集合s中ssIs是确定的值，我们可以进一步把问题转化成寻找一个排列P使得最大化。3.邻域定义在使用局部搜索方法前，我们首先定义 排列的邻域这一概念。我们称一个排列为 问题的一个解。设reads集合S是一个具有n 个元素的字符串集合，P(S)是reads集合S 所有解的集合。定义操作rshift(i，j) ，即 对一个解PP(S)，将P第i位置的元素向右 移至第j个位置(1i0时，p比p更好，则由p转移 到p此外，对给定的解P=s1，s2Sn,内部的元素均包含前趋和后继。但是，头元素仅有后继而尾元素仅有前驱。在算法中，为了消除头尾元素差异

14、，我们将排列的头尾元素连结起来并使用局部搜索方法寻找最优 的“Loop”超串。接着在overlap最小的地方切分该环状超串，最终还原成一条线性超串。5.局部极值点的处理当经过以read为单元的搜索后，可以获 得一个当前邻域内的局部最优解 1,k1,k1+1, ,k2,k2+1，k3km。 它对应集合s上的一个superstring。该解对 应的reads的排列中，任意相邻两个reads 间的overlop关系薄弱，即overlap(i， i+1)M1 ， overlap(r1，r3)M1，但overlap(r2，r3)50k)Contigs 平均长长度错误错误 contigs数 目错误错误 c

15、ontigs 总长总长 (bp)AE0078 725428697.5CAP35414652132876300 PHRAP5405881932845200 B-LSA5410821442864800 AE0006 5715513357.5CAP31539089777199886144539 PHRAP1528509717215003104164 B-LSA15459535011309192101102 INRUEN SS18301388CAP317957526812264085227941 PHRAP17986456314287484188705 B-LSA1811169441241162152

16、19 AE0007 8221784007.5CAP3216969310311215096279380 PHRAP21506509411223814212060 B-LSA2171203751328949339270 AE0078 6928415817.5CAP328269911261422436136274 PHRAP28155571161524272229589 B-LSA283368499203148500 AL0091 2642148147.5CAP3416548417110243586138549 PHRAP413779915519256955168824 B-LSA419066712321314849122422CAP3、PHRAP、B-LAS的性能比较结果一contigs的个数contig平均长度 B-LSA40833.661PHRAP53325.483CAP356623.035-结果：