[发明专利]基于双端读数insertsize分布的contig错误连接区域识别方法

摘要

本发明公开了一种基于双端读数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,and(PR+PL)/2＝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 and(PR＝NA or RR mapped in other contig)}|，表示左读数比对到p位置，右读数没有比对到这条contig的读数数量；DRp＝|{(PL,PR)|PR＝p and(PL＝NA or RL mapped in other contig)}|，表示右读数比对到p位置，左读数没有比对到这条contig的读数数量；步骤二：根据比对结果，得到双端支持稀疏的区域；将这些区域作为错误连接的候选区域；如果一对双端读数跨过某个contig区域，则称该对读数是对这个区域的一个支持；定义集合检查集合Z中的两两区间，若两个区间的距离小于较小区间的1/5，则将两个区间所包含的范围进行合并，得到更新后的集合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中剩余的部分即为确定的错误连接区域。