生物信息学序列拼接.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,基因组序列拼接,段倩倩,序列拼接,序列拼接任务即将测序生成的,reads,短片段拼接起来，恢复出原始的序列。该问题是序列分析的最基本任务，是基因组研究成功与失败的关键，拼接结果直接影响到序列标注，基因预测、基因组比较等后续任务。,基因组序列的拼接也是基因组研究必须解决的首要难题。其困难不仅来自它的海量数据,(,以人类基因组序列为例，从数量为,10,兆级的片断恢复出长度为亿级的原始序列,),，而且源于它含有高度重复的序列。,拼接问题的难点,DNA,测序数据有其固有的四个的特点，他们也正是解决实际的序列拼接问题的难点所在：,1,测序有误差,2,不完全覆盖性,3,序列所在链不确定,4,重复序列的干扰,1,测序有误差,由于测序技术的局限，难免会出现测序错误，尤其是在序列的末端，一般错误率可控制在,1,以下。所以对每个碱基一般有一个正确概率，以质量打分的形式给出。因此每个,r,i,都有个可信度。而,read,与,read,之间有不同程度的重叠，由此导致有的重叠可信度高，有的重叠可信度低。,2,不完全覆盖性,不是所有的碱基被测序的次数都等于平均测序覆盖度。极端的情况，可能会出现源基因组序列上部分区域未被测序的情况,(,这段区域称为,gap),。即，测序的,reads,集合不是原始基因组序列一个完整覆盖。此时需要借助于各种图谱如：基因组指纹图谱,(genome fingerprint map),，基因组级物理图谱,(genome-wide physical map),，细胞发生图谱,(cytogenetic maps),等协助对,reads,进行定位,3,序列所在链不确定,由于测序过程中无法确定特定片断属于,DNA,双链中的哪一条链上，所以我们在拼接过程中并不清楚使用的是,read,的正义链，还是其互补链。,4,重复序列的干扰,DNA,序列自身含有高度重复的子序列，它们一种表现为短序列的串级重复，比如：,(,GGAA),n,。或,A,m,T,n,等。另一种表现为大量相似序列,(,其拷贝数可达几十万),散布在基因组的各个地方。,Repeat,的存在，将导致,fragments,间,overlap,的不真实性，进而产生错拼的结果。因此在拼接过程中耍确定这些序列的形式及大小，才能保证以高概率恢复出其在原始真实序列中的位置,拼接算法评价,以上拼接问题的四个难点不仅极大的增加了解决实际拼接问题的难度，而且从某种程度上说无法完整地恢复出原始,DNA,序列来。即实际上仅能构建出若干个,contig,(,重建的,fragments,的一种排列形式，它覆盖基因组上一段连续区域,),这些,contig,将指导测序项目,finishing,阶段的实验方法最终构建,DNA,完整序列。,目前，国际上对拼接软件的公认评价标准包括两方面，即重建出的,contig,的数目和准确度。我们发展的基因组序列拼接新算法的目标是在确保准确性的前提下，构建尽量少的,contig,，以减少测序后期大量的人力和财力的投入。,基因组序列拼接算法研究现状,现在最常用的拼接程序使用的拼接算法可分成两类，一类是将拼接问题转化为在图中寻找的,Hamilton,路径的问题；另一类是将拼接问题在某种特殊情况下转化成寻求图中的,Euler,路径的问题。他们均有其成功的典型算法。,1.,转化为,Hamilton Path,问题,每个,DNA,片段,(read),相当于图中一个结点，如果两个片段之间存在着重叠,(overlap),关系，则在两个结点之间定义一条边，而沿着,DNA,原始序列从头到尾，则必然经过每个结点一次且仅一次，即是一条,Hamilton,路径。一条,contig,表示图中一条简单路，此类算法以,Phrap,，,TIGR Assembler,，,CAP3,，,GigAssemble,等为代表。,他们都是遵循“,overlap-layout-consensus”,的框架。首先，为了构建图。计算任意两个,read,间可能的比对情况。其次，通过去除歧义的或者不确信的边得到较为准确的图，并在其上寻找非交叉的简单路的集合，该集合对应于,contig,的集合。最终，通过对包含在一个简单路上的所有,read,进行多序列比对，为每一个,contig,构建一个一致性序列,(consensus sequence),。,2.,转化为,Euler Path,问题,EULER,是这类算法的代表。与传统方法沿着“,OverlapLayoutConsensus”,路线不同，它不计算各个,read,之间的,Overlap,，即没有,Overlap,步骤。,它的大致想法如下：,为了排除,read,中的错误，获得,Error-Free,的,read,，将所有的,read,切割成小片,n-mers,。,将每个,read,和,G,k,的近似进行比对，寻求,read,的最小改变能够使得,read,的所有,n-mers,包含在,G,k,的近似集合中。从而构建了高质量序列，而对于,Poor read,，直接抛弃，对,Chimeric,read(,两端在,n-mers,中但整体不在的,reads),进行特殊处理。,初始的想法是要实现去除,reads,中的测序错误的目的，如果知道原始序列,G,，那么直接使用测序获得的,read,和,G,进行比较即可。,但是实际上,G,并不可知，那么退而求其次，,G,的序列片断,G,k,亦可，事实上,G,k,亦不可知。所以将所有的,read,切割成小片,n-mers,，所有,Solid,的,n-mers,形成的集合称为,G,k,的近似。最后，构造,De,Bruijn,图。,现有算法的主要问题,虽然已经开发了以上的算法，基因组序列拼接问题尚未彻底解决，以上两类算法都存在着各自的缺陷。,对于第一类算法来说，实际上是在图中寻找一条使得评价函数值最优的,Hamilton,路径，这是一个,NP,完全问题。,一般都采用,greedy-merging,的算法近似求解。由于这种,step-by-step,的局部贪心算法，其明显的局部特性忽略了,reads,间“长距离”或者整体性的联系，从而导致了拼接错误，即拼接结果和真实的,DNA,原始序列不同。最近研究指出，在对已知序列的流行性感冒嗜血杆菌基因组的拼接过程中，无论是,Phrap,，,TIGR Assembler,，还是,CAP3,，都发生了拼接错误的现象。,对于第二类算法来说，它只能在特殊的情况下，才能将问题简化成寻求一条,Euler,路径，最终的结果是从多条候选的,Euler,超路中选择出来的。,EULER,算法依然存在拼接错误，且结果选择的过程没有理论依据。,EULER,软件在实际数据集的运行速度上和第一类算法相当。更重要的是，,EULER,采用的算法过于独立，很难利用其他辅助生物信息，导致其实用性和流行性大打折扣。,局部搜索,(Local Search),方法,胡杰,将局部搜索方法运用于一个具体的问题，需要对以下四项进行明确的定义。,1,将原同题表示成一个优化问题，即定义一个可行域以及在该可行域上的一个目标函数。,2,定义可行域中的邻域结构即说明满足什么条件的两个点相邻。,3,确定在邻域中的搜索方式。,4,局部极值点的处理。如果当前解点邻域中的所有点的目标函数值都比当前点大，那么这点就称为局部极值点。在一些问题中，局部极值点就是全局最优点。而对另一些问题而言，局部极值点已经满足求解实际问题的需求。,基于局部搜索的序列拼接算法框架,主要目标,是在,layout,阶段尽可能准确的前提下，获得更长的,contig,。具体的，使用局部搜索算法求得数据集上近似全局的最短超串；然后根据求得的最短公共超串对应的,fragment,的排列关系为基础获得“,consensus segment”,1.Overlap,定义,如果一个串的前缀是另一个串的后缀则认为这两个串之间存在,overlap,，并根据,over-lap,构建超串。对给定的串,f,和,g,，存在多个可能的,overlap,关系,比如说，若,f=ACTGGGAGCAGC,，,g=AGCAGCTTTTACT,，,那么他们之间至少存在两个,overlap,形式。,在我们的算法中，仅考虑两个串之间最大的,overlap,情况，并定义,overlap(f,，,g),表示,f,和,g,之间存在的多个,overlap,关系中最长的一个,overlap,所包含的字符个数。,在上面的例子中,overlap(f,，,g)=6,。如果,f,和,g,之间,overlap,区域长度小于,M(M,是一个足够小的正整数,),，则,overlap(f,，,g)=0,。,2.,优化目标定义,我们对,reads,集合,S,中的每个元素按照其左端在超串,t,中首次出现的位置进行排序，并沿超串,t,从左向右依次读入，并将其对应为序,S=,s,l,，,s,2,，,S,n,),。我们用,P(S),表示这个由字符串为元素构成的序。,在序,P(S),中，对于每一个连续的字符串元素对,s,i,，,s,i+1,，,都存在,overlap(s,i,，,S,i+1,),。,因此，字符串的一个排列等价于一个超串,t,及作用在其上的函数,overlap,，超串,t,的长度,length(t,)=t,=,s,s,Is,overlap(s,i,，,S,i+1,),为了求超串,t,具有最小的长度。因为在给定集合,s,中,s,s,Is,是确定的值，我们可以进一步把问题转化成寻找一个排列,P,使得,最大化。,3.,邻域定义,在使用局部搜索方法前，我们首先定义排列的邻域这一概念。我们称一个排列为问题的一个解。设,reads,集合,S,是一个具有,n,个元素的字符串集合，,P(S),是,reads,集合,S,所有解的集合。,定义操作,rshift(i,，,j),，即对一个解,PP(S),，将,P,第,i,位置的元素向右移至第,j,个位置,(1i,jn,),。,也就是说，如果我们应用这一操作于给定的解,P=,sl,s2,，,sn,),，则结果,P=,rshift(i,，,j)(P,),定义为；,类似的，定义操作,lshift(i,j,),将,P,第,i,位置的元素,向左,移至第,j,个位置,(1j,in,),，定义,P=,lshift(i,，,j)(P,),为；,这样，一个解,PP(s,),的邻域可以定义为,N(P)=,rshift(i,，,j)(P),|,1i,jn,),lshift(i,，,j)(P),|1,j0,时，,p,比,p,更好，则由,p,转移到,p,此外，对给定的解,P=s,1,，,s,2,S,n,内部的元素均包含前趋和后继。但是，头元,素仅有后继而尾元素仅有前驱。在算法中，为了,消除头尾元素差异,，我们将排列的头尾,元素连结起来并使用局部搜索方法寻找最优的“,Loop”,超串。接着在,overlap,最小的地方切分该环状超串，最终还原成一条线性超串。,5.,局部极值点的处理,当经过以,read,为单元的搜索后，可以获得一个当前邻域内的局部最优解,1,k,1,k,1,+1,k2,k,2,+1,，,，,k,3,k,m,。,它对应集合,s,上的一个,superstring,。该解对应的,reads,的排列中，任意相邻两个,reads,间的,overlop,关系薄弱，即,overlap(,i,，,i,+1,)M(,其中，,M,是一个足够小的正整数，,1i n),。,例如：,如果,overlap(r,kl,，,r,kI+1,)M,，,overlap(r,k,2,，,r,k,2,+1,)M,，,overlap(r,k,m,，,r,k,m,+1,)M1,，,overlap(r,1,，,r,3,)M1,，但,overlap(r,2,，,r,3,)M2,和,overlap(r,3,，,r,2,)50k),Contigs,平均长度,错误,contigs,数目,错误,contigs,总长,(,bp,),AE007872,542869,7.5,CAP3,541465,21,3,28763,0,0,PHRAP,540588,19,3,28452,0,0,B-LSA,541082,14,4,28648,0,0,AE000657,1551335,7.5,CAP3,1539089,77,7,19988,6,144539,PHRAP,1528509,71,7,21500,3,104164,B-LSA,1545953,50,11,30919,2,101102,INRUENSS,1830138,8,CAP3,1795752,68,12,26408,5,227941,PHRAP,1798645,63,14,28748,4,188705,B-LSA,1811169,44,12,41162,1,5219,AE000782,2178400,7.5,CAP3,2169693,103,11,21509,6,279380,PHRAP,2150650,94,11,22381,4,212060,B-LSA,2171203,75,13,28949,3,39270,AE007869,2841581,7.5,CAP3,2826991,126,14,22436,1,36274,PHRAP,2815557,116,15,24272,2,29589,B-LSA,2833684,99,20,31485,0,0,AL009126,4214814,7.5,CAP3,4165484,171,10,24358,6,138549,PHRAP,4137799,155,19,25695,5,168824,B-LSA,4190667,123,21,31484,9,122422,CAP3,、,PHRAP,、,B-LAS,的性能比较,结果一,contigs,的个数,contig,平均长度,B-LSA,408,33.661,PHRAP,533,25.483,CAP3,566,23.035,-,结果：,产生更少的,contig,，而且,contig,更长,有效的减少了测序项目后续工作的工作量,结果二,总,contigs,数,错误数目,错拼率,B-LSA,408,15,3.67,CAP3,745,52,6.97,PHRAP,566,24,4.24,-,结果：,包含更少的错误拼写现象,基于原型算法框架实现的,Basic LSA,系统在针对实际序列模拟的,shotgun reads,数据集上的拼接结果表明：,在不引入其他辅助信息的前提下，,Basic LSA,比,PHRAP(1999),及,CAP3,都有明显的改进这也证明了原型算法的有效性和可行性,基于原型算法的优化,原型算法的运行速度过慢，以至成为整个系统运行的,瓶颈,，从而导致该系统无法满足实验和可能的进一步应用的需要这也成了将原型系统实用化必须解决的问题,.,提速优化的目标是在不改变运行结果的前提下，提高系统运行速度,“邻域剪枝”,-(Neighborhood Pruning),尽可能通过排除那些目标函数值低于当前值的邻接点来缩小搜索空间这里的关键在于找到了一个关于能够增加目标函数值的变换的必要条件,After transposition,“领域剪枝”（,Neighborhood Pruning,）方法,Befor,transposition,在搜索过程中，每一个,transposition,将导致目标函数改变,overlap,，,Aoverlap,=Overlap4+Overlap5+Overlap6 Overlapl-Overlap2-Overlap3,若该,transposition,操作能改进当前解，必定有,Aoverlap,0,因此，也必然满足以下条件,(Overlap4-Overlapl-Overlap2)+(Overlap5+Overlap6)Overlap30,两种可能性：,1,Overlap4Overlapl+Overlap2,基于上述,Neighborhood Pruning,方法，我们优化了原型算法，由于局部搜索,(Local Search),是一种启发式的算法，其算法复杂度难以精确估计，随着真实数据特性的不同存在显著的差异,.,我们算法优化后的改进程度也会有所波动。我们选择了三个具有代表性的数据集，描述在不改变算法运行结果的前提下，提速优化前后的运行时间的比较,原始序列,序列长度,覆盖度,系统,Overlap,阶段用时,Layout,阶段用时,Conmmous,阶段用时,系统总用时,LSA layout,阶段优化后,Artificial,3M,8,PHRAP,300,4,790,1390,718,B-LSA,557,13506,434,14497,Im,-LSA,692,19,582,1293,BA000004,4.2M,8,PHRAP,1379,38,718,2135,191,B-LSA,1381,8234,664,6898,Im,-LSA,1899,43,517,1964,B4,2.76M,23,PHRAP,10730,757,2802,14289,89,B-LSA,10804,169385,2788,182975,Im,-LSA,10861,1903,2821,15585,LSA,提速优化前后性能比较,经过提速优化之后，,Layout,阶段不再是整个系统的时间瓶颈因此整个系统的运行速度与,Phrap,大致相当，从而朝着实现全基因组序列拼接系统的终极目标走出了关键一步,Thank You,！,