申请/专利权人：中南大学

申请日：2016-03-17

公开（公告）日：2018-03-06

公开（公告）号：CN105787295B

主分类号：G06F19/20(2011.01)I

分类号：G06F19/20(2011.01)I

优先权：

专利状态码：有效-授权

法律状态：2018.03.06#授权;2016.08.17#实质审查的生效;2016.07.20#公开

摘要：本发明公开了一种基于双端读数insert size分布的contig错误连接区域识别方法，包括以下步骤：1输入contigs集合和双端读数文库，使用序列比对工具将双端读数文库的双端读数比对到contigs集合上，得到比对结果；2根据比对结果，得到双端支持稀疏的区域；将这些区域作为错误连接的候选区域；3并通过双端读数的分布检验对候选区域进行延伸，最终通过区域长度判定候选区域是否是错误连接位置；4确定错误连接区域的边界。本发明方法具有较高的准确度，通过错误位点切割能够明显减少contig中的拼接错误，有效地提高了contig的质量。

主权项：一种基于双端读数insert size分布的contig错误连接区域识别方法，其特征在于，包括以下步骤：步骤一：输入contigs集合和双端读数文库，使用序列比对工具将双端读数文库的双端读数比对到contigs集合上，得到比对结果；contigs表示多个基因组连续片段，contig表示一个基因组连续片段；对每一对双端读数RL，RR，记它们比对到contig上的位置为PL，PR，若PL或PR未比对到任何位置，记位置为NA；根据比对结果，在每条contig的每个位点p维持3个集合，分别记为SMp、SLp和SRp：SMp＝{PR‑PL|RL,RR mapped accordantly,andPR+PL2＝p}，表示以p为中点，一致地比对到该contig上的所有双端支持中，每对双端读数之间的距离，即insert size；SLp＝{PR‑PL|RL,RR mapped accordantly,and PL＝p}，表示以p为左顶点，一致地比对到该contig上的所有双端支持中，每对双端读数之间的距离，即insert size；SRp＝{PR‑PL|RL,RR mapped accordantly,and PR＝p}，表示以p为右顶点，一致地比对到该contig上的所有双端支持中，每对双端读数之间的距离，即insert size；并记DLp＝|{PL,PR|PL＝p andPR＝NA or RR mapped in other contig}|，表示左读数比对到p位置，右读数没有比对到这条contig的读数数量；DRp＝|{PL,PR|PR＝p andPL＝NA or RL mapped in other contig}|，表示右读数比对到p位置，左读数没有比对到这条contig的读数数量；步骤二：根据比对结果，得到双端支持稀疏的区域；将这些区域作为错误连接的候选区域；如果一对双端读数跨过某个contig区域，则称该对读数是对这个区域的一个支持；定义集合检查集合Z中的两两区间，若两个区间的距离小于较小区间的15，则将两个区间所包含的范围进行合并，得到更新后的集合Z’；取集合Z’中区间长度大于指定长度的元素，组成候选区域的集合C；步骤三：根据比对结果，计算每一个候选区域比对的不一致率，并使用假设检验计算候选区域附近的双端读数insert size服从已知正态分布Nμ,σ2的概率，根据计算出的不一致率和假设检验结果决定是否对各候选区域r向左或向右延伸μ2，μ为已知的正态分布均值；根据各候选区域的最终长度判定候选区域是否是错误连接位置，如果不是，则将其从候选区域集合C中去除；步骤四、根据比对结果确定错误连接区域的边界；所述步骤三具体包括以下步骤：3.1对C中的每个候选区域r＝[a,b]，在候选区域左边，定义集合：SL＝{x|x∈SLa‑μ2∪SLa‑μ2+1∪..∪SLa‑2∪SLa‑1}，表示所有左端落在候选区域外面并且离候选区域距离小于μ2的双端读数中，每对双端读数之间的距离；计算不一致比对率其中，|SL|表示所有左端落在候选区域外面并且离候选区域距离小于μ2的双端读数数量；若u小于指定阈值，则放弃延伸左边；否则，使用Kolmogorov‑Smirnov检验计算SL服从Nμ,σ2的概率；若概率小于指定概率阈值，则认为SL不服从Nμ,σ2，支持了这个区域是错误连接位置，将候选区域r向左延伸μ2；3.2对C中的每个候选区域r＝[a,b]，在候选区域右边，定义集合：SR＝{x|x∈SRb+1∪SRb+2∪..∪SRb+μ2‑1∪SRb+μ2}，表示所有右端落在候选区域外面并且离候选区域距离小于μ2的双端读数中，每对双端读数之间的距离；计算不一致比对率其中，|SR|表示所有右端落在候选区域外面并且离候选区域距离小于μ2的双端读数数量；若u小于指定阈值，则放弃延伸右边；否则，使用Kolmogorov‑Smirnov检验计算SR服从Nμ,σ2的概率；若概率小于指定概率阈值，则认为SR不服从Nμ,σ2，支持了这个区域是错误连接位置，将候选区域r向右延伸μ2；根据各候选区域的最终长度判定候选区域是否是错误连接位置，如果候选区域的最终长度小于μ，说明该候选区域不是错误连接位置，则将其从候选区域集合C中去除；所述步骤四具体包括以下步骤：4.1在候选区域r＝[a,b]的左端，取μ2的长度，定义集合BL＝{x|x∈SRa∪SRa+1∪..∪SRa+μ2‑1∪SRa+μ2}，执行下列步骤：i.检查BL是否服从Nμ,σ2，若服从，则跳过步骤ii，进入步骤iii；否则，初始化平移步长为μ4，将BL向左移动一步，同时将步长减半，进入步骤ii；ii.检查BL是否服从Nμ,σ2，若服从，则将BL向右平移一步，并将步长减半；否则，将BL向左平移一步，并将步长减半；重复ii直到步长减小到步长指定值，进入步骤iii；iii.在r中移除被BL包含的部分；4.2在候选区域r＝[a,b]的右端，取μ2的长度，定义集合BR＝{x|x∈SLb∪SLb‑1∪..∪SLb‑μ2+1∪SLb‑μ2}，执行下列步骤：i.检查BR是否服从Nμ,σ2，若服从，则跳过步骤ii，进入步骤iii；否则，初始化平移步长为μ4，将BR向右移动一步，同时将步长减半，进入步骤ii；ii.检查BR是否服从Nμ,σ2，若服从，则将BR向左平移一步，并将步长减半；否则，将BR向右平移一步，并将步长减半；重复ii直到步长减小到步长指定值，进入步骤iii；iii.在r中移除被BR包含的部分；候选区域r中剩余的部分即为确定的错误连接区域。

全文数据：基于双端读数insertsize分布的contig错误连接区域识别方法技术领域[0001]本发明属于生物信息学领域，涉及基于双端读数insertsize分布的contig错误连接区域识别方法。背景技术[0002]获得一个物种的基因组是许多生命科学领域研究的基本前提，由于受物种基因组的复杂性以及测序技术未解决的测序长度制约，当前，对测序出的短序列进行组装是必不可少的任务。第三代测序技术由于测序成本高且错误率高，目前主流的测序还是采用以illumina基因测序仪为代表的第二代技术，或是二三代混合测序再拼接。第二代测序技术的特点是通量大，读长短，测序错误率低。iIIumina基因测序仪可以产生pairedend或matepair双端测序读数，解决一部分由读长过短带来的问题。Pairedend读数和matepair读数的每一对由两条读数组成，左读数的左端点到右读数的右端点的距离称为insertsize插入片段长度）。已有研究表明，一个文库的insertsize服从正态分布。采用第二代读数的组装算法可以分为两类:基于overlap-layout-consensus0LC，重叠-排列-生成一致序列）的方法和基于deBruijn图的方法。OLC方法第一步先计算读数与读数之间的所有交叠，接着根据交叠关系构造一个图，最后从图中获得长序列。deBruijn图采用一种违反直觉的方法，首先将读数切成更短的k-mer，再利用k-mer的交叠构造deBruijn图，接着从deBruijn图中得到拼接序列。[0003]OLC方法在测序倍数低、读长长的小型基因组上能很好地工作，但随着倍数和基因组规模增加，OLC方法需要的计算时间和内存太大，很难在普通机器上运转。而deBruijn图方法则可以在测序倍数高、读长readlength短的大型基因组上很好地工作，而随着测序成本的降低，数据的通量一直在提高，并且随着生物实验的需要，越来越复杂的物种也需要进行基因组测序。因此，在当前环境下，基于deBruijn图的方法被更为广泛地采用。[0004]基于deBruijn图方法，DanielR.Zerbino等人于2008发布了组装工具Velvet，通过双端读数，解决了一些长度小于insertsize双端读数中每对读数之间的间距）的repeat重复）问题，提高了拼接质量。而后，JaredT.Simpson等人发布了ABySS，提出了一种并行化方法，使复杂物种的大规模测序数据的组装成为可能。SanteGnerre等人于2011在ALLPATHS，ALLPATHS2的基础上发布了ALLPATHS-LG，针对哺乳动物基因组规模大，重复区复杂的特点，结合overlapping文库、jumpping文库以及PacBio数据，在高倍的哺乳动物数据集上取得了很好的结果。2011年，PramilaNuwanthaAriyaratne等人发布了PE-Assembler，提出了一种使用双端读数进行扩展的方法，取代传统的图遍历方法，在控制内存耗用的同时，提高了准确性，并使并行化变得简单。2014年，JunweiLuo等人发布了组装算法EPGA，通过双端读数的分布，提高了序列扩展时的准确性，提高了拼接质量。[0005]尽管可用的拼接工具越来越多，序列组装依然存在几个挑战：一个是物种自身基因组的复杂性，包括两个因素：重复区和多倍体的杂合度。对于重复区问题，一旦重复区长度大于双端数据的insertsize，则无法可靠地对重复区附近的序列进行链接，对于杂合度问题，杂合度对基因组组装的影响主要体现在不能合并同源染色体，杂合度高的区域，会把多条同源染色体都组装出来，同时也增加了deBruijn图的复杂性，使得获取高质量contig更加困难。另一个是测序不均问题，这是第二代测序技术的测序过程中PCR阶段引入的问题。这种不均会导致一些区域的覆盖非常低，以至于这些少量读数很可能被当成错误丢弃；而在另一些区域非常高，导致拼接过程中很可能产生误判。由于这些问题的存在，目前许多复杂物种的序列组装结果存在拼接碎片化，长度过短，拼接中带有错误等问题。[0006]此外，由于生物本身的复杂性和测序数据的不同规模，一个拼接工具往往在某些物种、某些数据上拼接结果很好，在另一些物种的数据上拼接结果就很差。为此，一些致力于解决这个问题的contig集成工具也陆续发布出来。包括GAM，MAIA，minimus2,CISA等，这些contig整合工具的初衷是集各个拼接工具输出的contig之所长，得到更好的contig。理想情况下，通过集成确实可以有效提高contig的质量，而现实中，由于拼接工具不能保证拼出的结果完全正确无误，一旦错误累积，将导致contig质量严重下降，因此，如何避免contig中的错误积累到整合结果中是不得不考虑的问题。发明内容[0007]本发明提出了一种新的基于双端读数insertsize分布的contig异常区域识别方法，通过将双端读数比对到contigs【基因组连续片段】上，找到双端支持稀疏的区域作为种子区域候选错误区域），并通过双端读数insertsize的分布检验对种子区域进行延伸，最终通过区域长度确定异常位点。本发明对contig中的错误识别具有较高的准确度，通过错误位点切割能够明显减少contig中的拼接错误，有效地提高了contig的质量。[0008]本发明的技术方案为：[0009]—种基于双端读数insertsize分布的contig错误连接区域识别方法，包括以下步骤：[0010]步骤一:输入contigs集合和双端读数文库，使用序列比对工具将双端读数文库的双端读数比对到contigs集合上，得到比对结果；[0011]对每一对双端读数Rl，Rr，记它们比对到contig上的位置为PL，Pr，若Pl或比对到任何位置，记位置为NA;[0012]根据比对结果，在每条contig的每个位点P维持3个集合，分别记为SMp、SLt^PSRp:[0013]SMp={PR-PLIRl，Rrmappedaccordantly，andPR+PL2=p}，表不以P为中点，一致地比对到该contig上的所有双端支持中，每对双端读数之间的距离，即insertsize;[0014]SLp={PR_PLIRl,Rrmappedaccordantly,andPL=p}，表示以p为左顶点，一致地比对到该contig上的所有双端支持中，每对双端读数之间的距离，即insertsize;[0015]SRp={PR-PLIRl,Rrmappedaccordantly,andPR=p}，表不以p为右顶点，一致地比对到该contig上的所有双端支持中，每对双端读数之间的距离，即insertsize;[0016]并记DLp=I{PL，PRIPL=pandPr=NAorRrmappedinothercontig}I，表示左读数比对到P位置，右读数没有比对到这条contig的读数数量；[0017]DRp=I{PL，PRIPR=pandPl=NAorRlmappedinothercontig}I，表不右读数比对到P位置，左读数没有比对到这条contig的读数数量；[0018]步骤二:根据比对结果，得到双端支持稀疏的区域;将这些区域作为错误连接的候选区域；[0019]如果一对双端读数跨过某个contig区域，则称该对读数是对这个区域的一个支持；[0020]定义集启[0021]检查集合Z中的两两区间，若两个区间的距离小于较小区间的15,则将两个区间所包含的范围进行合并，得到更新后的集合Z’；[0022]取集合Z’中区间长度大于指定长度的元素，组成候选区域的集合C;C中的每个区间都是双端支持相对其他位置薄弱的区域，因此有更大的可能性是错误连接的位置。[0023]步骤三:根据比对结果，计算每一个候选区域比对的不一致率，并使用假设检验计算候选区域附近的双端读数insertsize服从已知正态分布的概率，根据计算出的不一致率和假设检验结果决定是否对各候选区域r向左或向右延伸μ2;根据各候选区间的最终长度判定候选区域是否是错误连接位置，如果不是，则将其从候选区域集合C中去除；[0024]步骤四、根据比对结果确定错误连接区域的边界，边界确定后，即确定了所有错误连接区域;将所有错误连接区域剪除。[0025]所述步骤三包括以下步骤：[0026]3.1如图2a，对C中的每个区间r=[a,b]，在区间左边，定义集合：[0027]Sl={x|xeSLa—u2USLa—u2+1U··USLa-2USLa—，表示所有左端落在候选区间外面并且离候选区间距离小于μ2的双端读数中，每对双端读数之间的距离；[0028]计算不一致比对率[0029]其中，IS11表示所有左端落在候选区间外面并且离候选区间距离小于μ2的双端读数数量；[0030]若u小于指定阈值，则放弃延伸左边；否则，使用Kolmogorov-Smirnov检验计算Sl服从Nμ，〇2的概率p-value;若概率小于指定概率阈值，则认为S1不服从Nμ，〇2，支持了这个区域是错误连接位置，将区间r向左延伸μ2，μ为已知的正态分布均值；[0031]3.2同样，如图2⑹，对C中的每个区间r=[a,b]，在区间右边，定义集合：[0032]Sr={x|xeSRb+iUSRb+2U..USRbW2-iUSRbW2}，表示所有右端落在候选区间外面并且离候选区间距离小于μ2的双端读数中，每对双端读数之间的距离；[0033]计算不一致比对率[0034]其中，ISRI表示所有右端落在候选区间外面并且离候选区间距离小于μ2的双端读数数量；[0035]若u小于指定阈值，则放弃延伸右边；否则，使用Kolmogorov-Smirnov检验计算31?服从Nμ，〇2的概率p-value;若概率小于指定概率阈值，则认SSr不服从Nμ，〇2，支持了这个区域是错误连接位置，将区间r向右延伸μ2，μ为已知的正态分布均值；[0036]根据各候选区间的最终长度判定候选区域是否是错误连接位置，如果候选区间的最终长度小于μ，说明该候选区域不是错误连接位置，则将其从候选区域集合C中去除;基于假设，一个错误连接至少将导致长度μ的范围双端分布出现问题。[0037]所述步骤四包括以下步骤：[0038]为了确定错误连接区域的边界，依然使用双端读数分布；[0039]4.1如图3所示，在区间r=[a,b]的左端，取μ2的长度，定义集合Bl={xIXeSRaUSRa+iU..USRa+u2—IUSRa+u2}，执行下列步骤：[0040]i.检查Bl是否服从Νμ，σ2，若服从，则跳过步骤ii，进入步骤iii;否则，初始化平移步长为μ4，将Bl向左移动一步，同时将步长减半，进入步骤ii;[0041]ii.检查Bl是否服从Nμ，σ2，若服从，则将Bl向右平移一步，并将步长减半；否则，将Bl向左平移一步，并将步长减半;重复ii直到步长减小到步长指定值，进入步骤iii;[0042]iii.在r中移除被Bl包含的部分；[0043]4.2对于区间r的右边，也使用同样的方法确定边界。[0044]在区间r=[a，b]的右端，取μ2的长度，定义集合Br={xIxeSLbUSLb—IU..USLb—μ2+ιUSLb—μ2}，执行下列步骤：[0045]i.检查Br是否服从Nμ，σ2，若服从，则跳过步骤ii，进入步骤iii;否则，初始化平移步长为μ4，将Br向右移动一步，同时将步长减半，进入步骤ii;[0046]ii.检查Br是否服从Nμ，σ2，若服从，则将Br向左平移一步，并将步长减半；否则，将Br向右平移一步，并将步长减半;重复ii直到步长减小到步长指定值，进入步骤iii;[0047]iii.在r中移除被Br包含的部分。[0048]区间r中剩余的部分即为确定的错误连接区域；[0049]所述步骤一中的序列比对工具采用现有序列比对工具bowtie2。[0050]所述步骤一中的双端读数文库为jumpinglibrary，即跳查文库，双端距离较大的文库。[0051]所述步骤二中，取集合Z’中区间长度大于μ10的元素，组成候选区域的集合C，其中μ为已知的正态分布均值。[0052]所述步骤3.1和3.2中，指定阈值取全局不一致比对率％的4倍;全局不一致比对率Ug的计算方法为：[0053]Ug=1-NaN[0054]其中，N表示全部能比对到contigs上的读数数量单端能比对到contigs上的双端读数数量）；Na表示能一致比对到contigs上的读数数量双端能比对到contigs上的双端读数数量），即不仅读数本身比对到contig上，读数的对端读数mate也比对到同一条contig上，并且比对的方向和距离符合双端读数的方向和insertsize范围。[0055]所述步骤3.1和3.2中，指定概率阈值设为0.01。[0056]所述步骤四中，步长指定值取1。[0057]本发明：[0058]a提出了一种新的基于合并双端支持稀疏区间的候选区域筛选方法，如果一对双端读数跨过某个contig区域，则称该对读数是对这个区域的一个支持，在该方法中以双端读数的中点为代表保存这个支持。根据比对结果，在每个contig位置上保存以该位置为中点的双端支持数量，将连续支持数为零的位点合并成区间，再根据区间与区间的距离大小将邻近的区间进行合并，得到候选区域。该方法有效提高了错误识别的效率，并提高了识别准确度。[0059]b提出了一种新的基于双端支持分布的错误位置判定方法，对每一个候选区间，在候选区间的左边，收集所有左端落在候选区间外面并且离区间距离小于μ2的双端读数，这些双端读数有单端比对上的，也有双端比对上的，计算二者的比值为不一致比对率;再利用KoImogorov-Smirnov检验计算双端比对上的读数服从已知分布的p-value概率）。若不一致比对率大于阈值且p-value概率大于阈值，则将候选区间往左延伸μ2。在候选区域的右边做同样的计算。根据候选区间的最终长度对其进行取舍，得到待切除区间。[0060]提出了一种新的基于双端支持分布的边界确定方法，对每一个待切除区间，在区间的左边选取μ2长度的子区间，收集所有右端落在该子区间的双端支持，折半地左右平移该区间，直到区间内的双端支持服从已知分布，将这部分序列在切除区间去除，在待切除区间的右边做同样的计算，最终确定剪切边界。[0061]有益效果：[0062]本发明的方法不仅通过错误位点切割有效地减少c〇ntig中的拼接错误，对基于contig集成的序列拼接方法，本发明也能有效地避免这些错误积累到集成过程中，在高倍数据集上有效减少了contig连接错误的数量，为下游处理提供了更纯净的contig，有助于正确的contig进行合并，通过实验验证，有效提高了最终拼接的长度N50并减少了错误，提高序列拼接的质量和准确率。附图说明[0063]图l:Contig连接错误实例；[0064]图2:双端读数分布检验，图2a为区间左边检验，图2⑹为区间左边检验；[0065]图3:错误剪切边界确定；具体实施方式[0066]—、Contig错误分析[0067]无论是基于deBruijn图的de-novo拼接方法，还是基于OLC的拼接方法，从deBruijn图和overlap图中获得的长序列（contig通常很难保证完全正确，尤其是在基因组复杂、repeat重复）富集的物种中，错误更是难以避免。这些错误以contig的错误连接为主，即两段来自基因组不同位置的序列，被错误地拼在了一起，这种错误通常会在后续处理过程中保留下来，更严重的是，由于后续步骤对这种错误敏感，因此，contig中的错误很可能被放大。[0068]例如图1所示，使用某一拼接软件得到a，，c，d，e5条3〇111:188，这些contigs在基因组上的顺序为a，bl，c，d，b2，e，即，是错误连接的contig。[0069]这将导致几个问题：Icontig中的错误被继承到最终拼接的结果中，形成更长的带错误的序列，降低了拼接结果的可靠性;2错误的连接阻断了正确的连接，如contigbl和c本应是紧邻的，由于bl错误地和b2接在了一起，使得bl无法连接到c，阻碍了更长序列的形成，降低了拼接结果的质量;3Contig中的错误还可能在scaffolding超长序列片段组装阶段导致二次错误，例如，contigd的后续contig应为b2,但b2已经被contigbl粘住，这可能导致contigd选择另一个有少量双端支持的错误的后续contig，扩大了序列拼接错误。[0070]因此，减少contig中的错误连接具有重要意义。本文提出了一种基于双端读数分布的contig错误连接识别和纠正方法，有效地减少contig中的错误连接数，提高了contig的质量。经contig集成工具和scaffolding工具验证，有效提高了最终拼接出的序列长度并较少了最终结果的错误数。[0071]二、双端读数mapping[0072]使用现有序列比对工具bowtie2将jumpinglibrary跳查文库，即双端距离较大的文库）的双端读数比对到contig集合上，对每一对读数RlSRr，记它们比对到contig上的位置为P、PR，若P1SPr未比对到任何位置，记位置为NA。根据比对结果，在每条contig的每个位点P维持3个集合：[0073]SMp={PR-PLIRl，Rrmappedaccordantly,andPR+PL2=p}，表不以P为中点，一致地比对到该contig上的所有双端支持，这里只需要记录双端的距离，下同；[0074]SLp={PR-PLIRl，Rrmappedaccordantly,andPL=p}，表示以P为左顶点，一致地比对到该contig上的所有双端支持；[0075]SRp={PR-PLIRl，Rrmappedaccordantly,andPR=p}，表不以P为右顶点，一致地比对到该contig上的所有双端支持；[0076]并记DLp=I{PL，PRIPL=pandPr=NAorRrmappedinothercontig}I，表示左读数比对到P位置，右读数没有比对到这条contig的双端读数数量；[0077]DRp=I{PL，PRIPR=pandPl=NAorRlmappedinothercontig}I，表不右读数比对到P位置，左读数没有比对到这条contig的双端读数数量。[0078]三、确定候选区域[0079]定义集合·.检查集合中的两两区间，若两个区间的距离小于较小区间的15,则将两个区间所包含的范围进行合并。取区间长度大于指定长度预设值μ1〇,其中μ为已知的正态分布均值的元素，组成候选区域的集合CX中的每个区间都是双端支持相对其他位置薄弱的区域，因此有更大的可能性是错误连接的位置。[0080]如图2a，对C中的每个区间r=[a,b]，在区间左边，定义集合：[0081]Sl={xIXeSLa-u2USLa-u2+lU..USLa-2USLa-l}[0082]计算不一致比对率。若u小于指定阈值，则放弃延伸左边，其中阈值按比对结果的全局不一致率计算，预设值取全局不一致率的4倍。否则，使用KolmogoroV-5111；[1'110¥检验计算51服从1'1^，02的口-¥31116。若口-¥31116小于指定阈值取0.01，则认为51不服从Nμ，02，支持了这个区域是错误连接位置，将区间r向左延伸μ2，μ为已知的正态分布均值。[0083]同样，如图2⑹，在区间右边，定义集合：[0084]Sr={xIXeSRb+iUSRb+2U··USRb+u2-1USRb+u2}[0085]计算不一致比对率u。并同样根据u和假设检验结果决定是否对区间r向右延伸μ2〇[0086]最后去除区间长度小于μ的区间，基于假设，一个错误连接至少将导致长度μ的范围双端分布出现问题。[0087]四、确定剪切边界[0088]为了确定剪切的边界，依然使用双端读数分布。如图3所示，在区间r=[a，b]的左端，取μ2的长度，定义集合妒={x|xeSRaUSRa+1U..USRaWHUSRaW2I，执行下列步骤：[0089]i.检查Bl是否服从Νμ，σ2，若服从，则跳过步骤ii，进入步骤iii，否则，初始化平移步长为μ4，将Bl向左移动一步，同时将步长减半。[0090]ii.检查Bl是否服从Nμ，σ2，若服从，则将Bl向右平移一步，并将步长减半；否则，将向左平移一步，并将步长减半。重复ii直到步长减小到指定值【指定值取1】。[0091]iii.在r中移除被Bl包含的部分。[0092]对于区间r的右边，也使用同样的方法确定边界。[0093]边界确定后，即确定了所有异常区域。移除所有异常区域，将错误连接的contigs剪开，获得最终结果。[0094]五、实验验证[0095]为了验证本方法的有效性，我们使用PEAssemblIer、Velvet、SOAPdenovo和Abyss四个组装工具对三个物种的iIlumina测序数据进行组装产生contigs。然后用jumppinglibrary数据对产生的contigs分别进行改错，staph，ecoli，spome的jumppinglibrary倍数分别为70倍，61倍，179倍。然后，使用本发明的改错方法对其进行改错。为了检验改错的准确性，我们将改错过程中剪切下来的区域收集起来，并使用权威的序列拼接评价工具quast检查被剪除的区域是否是错误序列，并做了统计。[0096]对于每一组contigs，我们使用改错工具进行改错，统计改错的次数，同时，每次改错后我们都使用quast检查错误是否减少，即改错是否正确，结果如表1所示。[0097]表1·错误识别准确率表中每一项表不（errorconfirmedbyquasterrorcalledW]~从表I中可以看出，在3个数据集的12组测试实验中，有4组改正了contig的错误，_改错准确率均是100%。我们还可以看到，Spome数据集由于基因组复杂，更容易产生组装错误，并且由于jumppingIibrary倍数高，改错效果最好。[0Ί00]为了进一步验证contigs改错的现实意义，我们将上述改错前后的contig按每个物种进行集成。集成工具使用了CISA，使用quast评估整合后的contig。[0101]表2.改错前后的contig整合结果，使用整合工具CISA[0104]从表2中可以看到，staph数据集和spome数据集上，改错之后在保持拼接长度的前提下，减少了最终集成结果的拼接错误。在基因组最为复杂的spome物种上，不仅拼接错误减少了58,拼接长度和覆盖度也有了提高。从这组实验结果可以看出，本文提出的方法对拼接容易出错的复杂物种具有较好的适应性，能有效改善contig整合结果。[0105]Contigs改错不仅对contig集成具有重要意义，对scaffolding也会产生影响。我们选择上述实验中改错工具产生作用的4组contig用scaffold工具SSPACE进行scaffolding测试，使用quast统计了改错前后的scaffolding结果，如表3所示。[0106]表3.改错前后的scaffolding结果，使用scaffold工具SSPACE[0108]从表3中可以看出，在4组实验中，最终组装结果错误数都得到降低。不仅如此，第3组和第4组实验中，由于错误的剔除，组装过程中干扰因素减少，导致拼接长度也得到提高。

权利要求：I.一种基于双端读数insertsize分布的contig错误连接区域识别方法，其特征在于，包括以下步骤：步骤一:输入contigs集合和双端读数文库，使用序列比对工具将双端读数文库的双端读数比对到contigs集合上，得到比对结果;contigs表示多个基因组连续片段，contig表示一个基因组连续片段；对每一对双端读数R^Rr，记它们比对到contig上的位置为P、Pr，若PlSPr未比对到任何位置，记位置为NA;根据比对结果，在每条contig的每个位点p维持3个集合，分别记为SMP、SLp和SRp:,表示以P为中点，一致地比对到该contig上的所有双端支持中，每对双端读数之间的距离，S卩insertsize;，表示以P为左顶点，一致地比对到该contig上的所有双端支持中，每对双端读数之间的距离，S卩insertsize;，表示以P为右顶点，一致地比对到该contig上的所有双端支持中，每对双端读数之间的距离，S卩insertsize;并记，表示左读数比对到P位置，右读数没有比对到这条contig的读数数量；，表示右读数比对到P位置，左读数没有比对到这条contig的读数数量；步骤二:根据比对结果，得到双端支持稀疏的区域;将这些区域作为错误连接的候选区域；如果一对双端读数跨过某个contig区域，则称该对读数是对这个区域的一个支持；定义集忐检查集合Z中的两两区间，若两个区间的距离小于较小区间的15,则将两个区间所包含的范围进行合并，得到更新后的集合Z’；取集合Z’中区间长度大于指定长度的元素，组成候选区域的集合C;步骤三:根据比对结果，计算每一个候选区域比对的不一致率，并使用假设检验计算候选区域附近的双端读数insertsize服从已知正态分布Nμ，σ2的概率，根据计算出的不一致率和假设检验结果决定是否对各候选区域r向左或向右延伸μ2，μ为已知的正态分布均值;根据各候选区域的最终长度判定候选区域是否是错误连接位置，如果不是，则将其从候选区域集合C中去除；步骤四、根据比对结果确定错误连接区域的边界；所述步骤三具体包括以下步骤：3.1对C中的每个候选区域r=[a,b]，在候选区域左边，定义集合：，表示所有左端落在候选区域外面并且离候选区域距离小于μ2的双端读数中，每对双端读数之间的距离；计算不一致比对率其中，IS11表示所有左端落在候选区域外面并且离候选区域距离小于μ2的双端读数数量；若U小于指定阈值，则放弃延伸左边；否则，使用Kolmogorov-Smirnov检验计算31服从Nμ，〇2的概率;若概率小于指定概率阈值，则认为S1不服从Nμ，〇2，支持了这个区域是错误连接位置，将候选区域r向左延伸μ2;3.2对C中的每个候选区域r=[a,b]，在候选区域右边，定义集合：，表示所有右端落在候选区域外面并且离候选区域距离小于μ2的双端读数中，每对双端读数之间的距离；计算不一致比对；其中，ISrI表示所有右端落在候选区域外面并且离候选区域距离小于μ2的双端读数数量；若u小于指定阈值，则放弃延伸右边；否则，使用Kolmogorov-Smirnov检验计算31?服从Nμ，〇2的概率;若概率小于指定概率阈值，则认SSr不服从Nμ，〇2，支持了这个区域是错误连接位置，将候选区域r向右延伸μ2;根据各候选区域的最终长度判定候选区域是否是错误连接位置，如果候选区域的最终长度小于y，说明该候选区域不是错误连接位置，则将其从候选区域集合C中去除；所述步骤四具体包括以下步骤：4.1在候选区域r=[a，b]的左端，取μ2的长度，定义集合Bl={xIxeSRaUSRa+iU..USRa+u2-lUSRa+u2}，执行下列步骤：i.检查是否服从Nμ，〇2，若服从，则跳过步骤ii，进入步骤iii;否则，初始化平移步长为μ4，将Bl向左移动一步，同时将步长减半，进入步骤ii;ii.检查Bl是否服从Nμ，σ2，若服从，则将Bl向右平移一步，并将步长减半;否则，将Bl向左平移一步，并将步长减半;重复ii直到步长减小到步长指定值，进入步骤iii;iii.在r中移除被Bl包含的部分；4.2在候选区域『=|^，13]的右端，取42的长度，定义集合1^={14£51^1^1^-：^.^SLb—μ2+lUSLb—μ2}，执行下列步骤：i.检查Br是否服从Nμ，〇2，若服从，则跳过步骤ii，进入步骤iii;否则，初始化平移步长为μ4，将Br向右移动一步，同时将步长减半，进入步骤ii;ii.检查Br是否服从Nμ，σ2，若服从，则将Br向左平移一步，并将步长减半;否则，将Br向右平移一步，并将步长减半;重复ii直到步长减小到步长指定值，进入步骤iii;iii.在r中移除被Br包含的部分；候选区域r中剩余的部分即为确定的错误连接区域。2.根据权利要求1所述的基于双端读数insertsize分布的contig错误连接区域识别方法，其特征在于，所述步骤一中的序列比对工具采用现有序列比对工具bowtie2。3.根据权利要求2所述的基于双端读数insertsize分布的contig错误连接区域识别方法，其特征在于，所述步骤一中的双端读数文库为jumpingIibrary，即跳查文库。4.根据权利要求3所述的基于双端读数insertsize分布的contig错误连接区域识别方法，其特征在于，所述步骤二中，取集合Z’中区间长度大于μ10的元素，组成候选区域的集合C，其中μ为已知的正态分布均值。5.根据权利要求4所述的基于双端读数insertsize分布的contig错误连接区域识别方法，其特征在于，所述步骤3.1和3.2中，指定阈值取全局不一致比对率也的4倍;全局不一致比对率Ug的计算方法为：其中，N表示全部能比对到contigs上的读数数量;Na表示能一致比对到contigs上的读数数量，即不仅读数本身比对到contig上，读数的对端读数也比对到同一条contig上，并且比对的方向和距离符合双端读数的方向和insertsize范围。6.根据权利要求5所述的基于双端读数insertsize分布的contig错误连接区域识别方法，其特征在于，所述步骤3.1和3.2中，指定概率阈值设为0.01。7.根据权利要求6所述的基于双端读数insertsize分布的contig错误连接区域识别方法，其特征在于，所述步骤四中，步长指定值取1。

百度查询： 中南大学 基于双端读数insert size分布的contig错误连接区域识别方法

免责声明
1、本报告根据公开、合法渠道获得相关数据和信息，力求客观、公正，但并不保证数据的最终完整性和准确性。
2、报告中的分析和结论仅反映本公司于发布本报告当日的职业理解，仅供参考使用，不能作为本公司承担任何法律责任的依据或者凭证。

阅读全文 双屏查看 官方信息 专利公告 收藏专利 下载PDF 下载WORD