生物信息学 序列拼接

基因组序列拼接段倩倩序列拼接序列拼接任务即将测序生成的reads短片段拼接起来，恢复出原始的序列。该问题是序列分析的最基本任务，是基因组研究成功与失败的关键，拼接结果直接影响到序列标注，基因预测、基因组比较等后续任务。基因组序列的拼接也是基因组研究必须解决的首要难题。其困难不仅来自它的海量数据(以人类基因组序列为例，从数量为10兆级的片断恢复出长度为亿级的原始序列)，而且源于它含有高度重复的序列。拼接问题的难点DNA测序数据有其固有的四个的特点，他们也正是解决实际的序列拼接问题的难点所在：1．测序有误差2．不完全覆盖性3．序列所在链不确定4．重复序列的干扰1．测序有误差由于测序技术的局限，难免会出现测序错误，尤其是在序列的末端，一般错误率可控制在1％以下。所以对每个碱基一般有一个正确概率，以质量打分的形式给出。因此每个ri都有个可信度。而read与read之间有不同程度的重叠，由此导致有的重叠可信度高，有的重叠可信度低。2．不完全覆盖性不是所有的碱基被测序的次数都等于平均测序覆盖度。极端的情况，可能会出现源基因组序列上部分区域未被测序的情况(这段区域称为gap)。即，测序的reads集合不是原始基因组序列一个完整覆盖。此时需要借助于各种图谱如：基因组指纹图谱(genomefingerprintmap)，基因组级物理图谱(genome-widephysicalmap)，细胞发生图谱(cytogeneticmaps)等协助对reads进行定位．3．序列所在链不确定由于测序过程中无法确定特定片断属于DNA双链中的哪一条链上，所以我们在拼接过程中并不清楚使用的是read的正义链，还是其互补链。4．重复序列的干扰DNA序列自身含有高度重复的子序列，它们一种表现为短序列的串级重复，比如：(GGAA)n。或AmTn等。另一种表现为大量相似序列(其拷贝数可达几十万)散布在基因组的各个地方。Repeat的存在，将导致fragments间overlap的不真实性，进而产生错拼的结果。因此在拼接过程中耍确定这些序列的形式及大小，才能保证以高概率恢复出其在原始真实序列中的位置．拼接算法评价以上拼接问题的四个难点不仅极大的增加了解决实际拼接问题的难度，而且从某种程度上说无法完整地恢复出原始DNA序列来。即实际上仅能构建出若干个contig(重建的fragments的一种排列形式，它覆盖基因组上一段连续区域)这些contig将指导测序项目finishing阶段的实验方法最终构建DNA完整序列。目前，国际上对拼接软件的公认评价标准包括两方面，即重建出的contig的数目和准确度。我们发展的基因组序列拼接新算法的目标是在确保准确性的前提下，构建尽量少的contig，以减少测序后期大量的人力和财力的投入。基因组序列拼接算法研究现状现在最常用的拼接程序使用的拼接算法可分成两类，一类是将拼接问题转化为在图中寻找的Hamilton路径的问题；另一类是将拼接问题在某种特殊情况下转化成寻求图中的Euler路径的问题。他们均有其成功的典型算法。1.转化为HamiltonPath问题每个DNA片段(read)相当于图中一个结点，如果两个片段之间存在着重叠(overlap)关系，则在两个结点之间定义一条边，而沿着DNA原始序列从头到尾，则必然经过每个结点一次且仅一次，即是一条Hamilton路径。一条contig表示图中一条简单路，此类算法以Phrap，TIGRAssembler，CAP3，GigAssemble等为代表。他们都是遵循“overlap-layout-consensus”的框架。首先，为了构建图。计算任意两个read间可能的比对情况。其次，通过去除歧义的或者不确信的边得到较为准确的图，并在其上寻找非交叉的简单路的集合，该集合对应于contig的集合。最终，通过对包含在一个简单路上的所有read进行多序列比对，为每一个contig构建一个一致性序列(consensussequence)。2.转化为EulerPath问题EULER是这类算法的代表。与传统方法沿着“Overlap—Layout—Consensus”路线不同，它不计算各个read之间的Overlap，即没有Overlap步骤。它的大致想法如下：为了排除read中的错误，获得Error-Free的read，将所有的read切割成小片n-mers。将每个read和Gk的近似进行比对，寻求read的最小改变能够使得read的所有n-mers包含在Gk的近似集合中。从而构建了高质量序列，而对于Poorread，直接抛弃，对Chimericread(两端在n-mers中但整体不在的reads)进行特殊处理。初始的想法是要实现去除reads中的测序错误的目的，如果知道原始序列G，那么直接使用测序获得的read和G进行比较即可。但是实际上G并不可知，那么退而求其次，G的序列片断Gk亦可，事实上Gk亦不可知。所以将所有的read切割成小片n-mers，所有Solid的n-mers形成的集合称为Gk的近似。最后，构造DeBruijn图。现有算法的主要问题虽然已经开发了以上的算法，基因组序列拼接问题尚未彻底解决，以上两类算法都存在着各自的缺陷。对于第一类算法来说，实际上是在图中寻找一条使得评价函数值最优的Hamilton路径，这是一个NP完全问题。一般都采用greedy-merging的算法近似求解。由于这种step-by-step的局部贪心算法，其明显的局部特性忽略了reads间“长距离”或者整体性的联系，从而导致了拼接错误，即拼接结果和真实的DNA原始序列不同。最近研究指出，在对已知序列的流行性感冒嗜血杆菌基因组的拼接过程中，无论是Phrap，TIGRAssembler，还是CAP3，都发生了拼接错误的现象。对于第二类算法来说，它只能在特殊的情况下，才能将问题简化成寻求一条Euler路径，最终的结果是从多条候选的Euler超路中选择出来的。EULER算法依然存在拼接错误，且结果选择的过程没有理论依据。EULER软件在实际数据集的运行速度上和第一类算法相当。更重要的是，EULER采用的算法过于独立，很难利用其他辅助生物信息，导致其实用性和流行性大打折扣。局部搜索(LocalSearch)方法胡杰将局部搜索方法运用于一个具体的问题，需要对以下四项进行明确的定义。1．将原同题表示成一个优化问题，即定义一个可行域以及在该可行域上的一个目标函数。2．定义可行域中的邻域结构，即说明满足什么条件的两个点相邻。3．确定在邻域中的搜索方式。4．局部极值点的处理。如果当前解点邻域中的所有点的目标函数值都比当前点大，那么这点就称为局部极值点。在一些问题中，局部极值点就是全局最优点。而对另一些问题而言，局部极值点已经满足求解实际问题的需求。基于局部搜索的序列拼接算法框架主要目标是在layout阶段尽可能准确的前提下，获得更长的contig。具体的，使用局部搜索算法求得数据集上近似全局的最短超串；然后根据求得的最短公共超串对应的fragment的排列关系为基础获得“consensussegment”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={sl，s2…．，Sn)。我们用P(S)表示这个由字符串为元素构成的序。在序P(S)中，对于每一个连续的字符串元素对si，si+1，都存在overlap(si，Si+1)。因此，字符串的一个排列等价于一个超串t及作用在其上的函数overlap，超串t的长度length(t)=│t│=∑s∈sIs|-overlap(si，Si+1),为了求超串t具有最小的长度。因为在给定集合s中∑s∈sIs│是确定的值，我们可以进一步把问题转化成寻找一个排列P使得最大化。ni13.邻域定义在使用局部搜索方法前，我们首先定义排列的邻域这一概念。我们称一个排列为问题的一个解。设reads集合S是一个具有n个元素的字符串集合，P(S)是reads集合S所有解的集合。定义操作rshift(i，j)，即对一个解P∈P(S)，将P第i位置的元素向右移至第j个位置(1≤ij≤n)。也就是说，如果我们应用这一操作于给定的解P={sl，s2…．，sn)，则结果P’=rshift(i，j)(P)定义为；类似的，定义操作lshift(i,j)将P第i位置的元素向左移至第j个位置(1≤ji≤n)，定义P’=lshift(i，j)(P)为；这样，一个解P∈P(s)的邻域可以定义为N(P)={rshift(i，j)(P)|1≤ij≤n)｝∪{lshift(i，j)(P)|1≤ji≤n｝．4．局部搜索算法通过计算从可行解P到其邻居P’∈N(P)的过程中的∆overlap，并采用“first-fit的搜索策略，我们在当前解的邻域范围内搜索第一个更好的解，并转移到该解，因此可以通过迭代的方法改善解的质量，并实现潜在的纠错功能。例如，在γshift(i,j)操作下将p移至p’，∆overlap=overlap(si-1,si+1)+overlap(sj,sj+1)+overlap(si,sj+1)-overlap(si-1,si)-overlap(si,si+1)-overlap(sj,sj+1)。当∆overlap0时，p’比p更好，则由p转移到p’此外，对给定的解P={s1’，s2’…．Sn’｝,内部的元素均包含前趋和后继。但是，头元素仅有后继而尾元素仅有前驱。在算法中，为了消除头尾元素差异，我们将排列的头尾元素连结起来并使用局部搜索方法寻找最优的“Loop”超串。接着在overlap最小的地方切分该环状超串，最终还原成一条线性超串。5.局部极值点的处理当经过以read为单元的搜索后，可以获得一个当前邻域内的局部最优解{γ1…,γk1,γk1+1,…,γk2,γk2+1，…，γk3γkm｝。它对应集合s上的一个superstring。该解对应的reads的排列中，任意相邻两个reads间的overlop关系薄弱，即overlap(γi，γi+1)M(其中，M是一个足够小的正整数，1≤i≤n)。例如：如果overlap(rkl，rkI+1)M，overlap(rk2，rk2+1)M,…，overlap(rkm，rkm+1)M，令fragment集合F=｛f1,f2,f3….,fm}，fi即为在序rki-1,…,rki下对应的一个子超串。则在此薄弱处将superstring分割若干“sub-superstrings.它们构成新的fragment集合F。此时，我们将F作为当前最短超串求解的子串集合，可以获得一个新的是最短公共超串(SCS)问题的实例。我们继续应用前述的Overlap定义和邻域的定义。最终得到reads集合上的一个局部最优解P’。6．获得“consensussegments”在处理对“contiguoussegments”的识别过程中，有关键的两个步骤：(1)在求得的超串中识别“weakjoins”．(2)在“weakjoins”位置，获得contiguoussegments．7.Repeat的识别策略为了鉴别潜在的falsejoins，我们必须分析所有的包含不一致的前缀或者后缀的fragment。具体的fragm