CN108491687A - 基于contig质量评估分类及图优化的scaffolding方法 - Google Patents

Abstract

本发明公开了一种基于contig质量评估分类及图优化的scaffolding方法，采用序列比对信息以及contig的GC含量信息对contig集合进行质量评估并分类，再将每个contig作为一个节点，根据双端读数比对到contig上的数量期望值以及实际值之间的差异判断是否在两个节点之间构建边，并计算边的权值，构建加权的scaffold图。最后通过为节点分配方向以及剪切节点来消除scaffold图中的方向冲突，通过为节点分配顺序来消除scaffold图中的顺序冲突。本发明简单易用，在四组真实测序数据上表现出良好的拼接结果，较其他序列拼接方法具有更高的F‑score值。

Description

基于contig质量评估分类及图优化的scaffolding方法

技术领域

本发明属于生物信息学领域，涉及contig质量评估分类以及scaffold图优化的scaffolding方法。

背景技术

从头序列组装(De Novo Sequence Assembly)是基因组学领域一项重要的研究方向，同时也是基因组学下游分析的一项重要基础。基因组学对基因组的组成、组内各基因精确结构、表达调控以及相互关系等方面进行了研究，序列拼接作为研究的基础条件，其准确性对整个基因组学的研究起着关键作用。由于基因组DNA序列结构比较复杂，特别是重复区(即一段DNA片段多次出现在基因组的不同位置)问题，测序错误问题(即读数中包含一定的错误碱基)，以及读数长度问题等限制了序列组装方法的应用。

序列组装包括以下三大步骤：(1)contig构建阶段：一条contig就是一条DNA序列片段，是根据读数之间的重叠关系对种子序列进行左右扩展得到的较长的序列。目前已经提出了很多contig构建的方法，一种是基于读数重叠图的方法，另一种是基于De Bruijin图的方法。(2)scaffolding阶段：对于第一阶段产生的大量contig，本阶段确定这些contig的方向及顺序关系，从而产生长度更长的序列片段scaffolds，scaffolds之间的空白区域用“N”来填充。(3)gap填充阶段：该阶段确定scaffolds中gap区域的序列，进而减少scaffolds中未知区域的长度。

由于序列工具产生的contig可能分布在基因组序列的任意区域，并且由于DNA是双链结构，这些contig可能处于双链上的任意一条链上，如果两个contig处在同一条链上，那么这两个contig就是同向的。Scaffolding方法用来确定contig之间的方向以及顺序关系，将它们组装成一些更长的序列片段(scaffold)。Scaffolding基于双端读数以及contig集合的支持，是序列拼接过程中十分重要的阶段，scaffolding能够使序列组装的结果更连续更完整，有助于后续基因识别，基因组比对，结构变异检测等研究，是序列组装研究中的热点之一。由于第二代测序技术比较成熟，并且具有正确率高、成本低和通量高的优势，所以在国内外得到了广泛的应用。虽然第二代测序技术产生的读数比较短，但是测序得到的双端读数的插入长度可以达到数千碱基，能够克服重复去带来的问题。所以采用双端读数来推断contig之间的方向和顺序关系是scaffolding方法研究的热点。

现有的scaffolding方法通常可以分为两大类：

(1)基于图的scaffolding方法。其基本思路是：首先将双端文库比对到contig集合上，将contig视为顶点，然后根据双端读数的比对情况，对存在双端支持的contig之间加一条有向边，由比对结果可推导出边的方向，以此构建有向图。由于构建的原始有向图较为复杂，并且包含许多可靠性很差的边，因此需要将图进行化简。化简图的方法一般包括移除矛盾边和含有重复区的contig节点，以及子图分割等，它们的目的都是在保证高质量边不被删除的情形下降低图的复杂性。在化简图之后，再从图中抽取路径，每条路径代表一条scaffold。基于图的scaffolding算法主要有GRASS，MIP，SOPRA，Bambus2，SCARPA，Opera，SGA，ABySS，ScaffoldMatch等。

(2)基于贪心策略的scaffolding方法。其基本思路是：首先输入contig集合，选取最长的contig作为种子，再利用contig不断扩充原始序列，最终得到scaffolds集合。在扩充延伸的过程中，该方法会碰到测序不均以及多条候选延伸路径或多个候选延伸节点的问题，如何处理这种问题将成为整个算法输出结果质量好坏的关键。基于贪心策略的scaffolding算法只考虑相邻的候选节点，每次选取得分最高的候选节点作为延伸方向。基于贪心策略的scaffolding算法主要有PE、Bambus、SSPACE、GigAssembler，Huson的方法，ISEA等。

目前，基于图的scaffolding算法准确度较高，但都面临着以下三个问题：(1)如何解决由于测序不均衡导致的边添加错误。现有的scaffolding方法将比对上的双端读数的数目作为边的权值，当出现重复区就会导致重复区段的边权值较大，此时无法区分该段是重复区域还是可信度较高的区域。当出现测序深度较低就会导致该区段的边权值较小，此时无法区分该段是测序深度低的区域还是可信度较低的区域。(2)如何充分利用高权值边对其他边的影响作用。现有方法同时考虑所有边来对图进行优化，忽视了高权值边对于其他边的影响作用，权值较高的边往往具有较高的可信度，利用这些高权值边来判断其他边的方向及顺序更为准确(3)如何有效减少由于重复区、测序不均衡和测序错误导致的contig集合中的错误。contig集合的拼接质量直接关系到scaffolding环节的准确度，然而目前的scaffolding方法都是直接将contig集合作为输入，通过删减边的形式来优化最后的scaffold图，然而scaffold图中的很多冲突边是由于contig的错误拼接导致的，单纯的删减边并不能从根本上解决冲突，反而会使有用的比对信息减少。

发明内容

本发明所解决的技术问题是，针对现有技术的不足，公开了一种基于contig质量评估分类及图优化的scaffolding方法(SCOP)，本发明简单易用，在不同的真实测序数据上得到了较好的scaffolds结果，较其他scaffolding方法具有更高的F-score值，并且在CJ，IJ等其他评价指标上也同样取得了较好的结果。

本发明的技术方案为：

一种基于contig质量评估分类及图优化的scaffolding方法，包括以下步骤：

步骤1、数据预处理：首先将双端读数比对到已有的contig集合上，得到比对结果；然后对比对结果进行过滤：

1.1)对于双端读数中的每个读数，只保留比对得分值最高的比对位置信息；

1.2)判断比对到同一个contig上的双端读数的插入长度是否在[μ_is-3σ_is,μ_is+3σ_is]范围内，若超过该范围则删除该条比对信息；其中μ_is表示比对到同一个contig上的所有双端读数的插入长度的均值，σ_is表示比对到同一个contig上的所有双端读数的插入长度的方差；

1.3)基于读数比对信息计算contig上每个位置的读数覆盖度，以及该contig上所有位置读数覆盖度的平均值μ和标准差σ；若该contig上某个位置的读数覆盖度与μ的差值大于2σ，则删除该条比对信息；

步骤2、节点质量评估及分类：利用步骤1中保留的比对信息以及contig中的GC含量特征来对contig的质量进行评估，并根据质量评估结果将contig分为正确(True)，错误(False)，不确定(Uncertain)三大类；

步骤3、构建加权的scaffold图：将每个contig作为一个节点，根据双端读数比对到contig上的数量期望值以及实际值之间的差异判断是否在两个节点之间构建边，并且根据一系列统计学方法计算边的权值；将所有的边按照权值从大到小进行排序，以便后续进行优化；

步骤4、对scaffold图进行优化：在比对信息的基础上结合节点的分类与边的权值对scaffold图进行优化；

本发明提出了一种新的优化scaffold图的方法，该方法分为节点方向的优化以及边顺序的优化。在节点方向优化的环节，通过给每个节点预分配一个方向，找出造成方向冲突的边，若该边连接的节点中存在不确定类型的节点，则分割该节点，再重新构造scaffold图，以检测是否存在方向冲突。本方法是一个迭代的过程，每次选择一个阈值，若边的权值大于该阈值则加入检测分析，每次迭代结束该阈值减0.1。

步骤5、从图中提取scaffolds：采用广度遍历的方法从优化后的scaffold图中尽可能提取较长的scaffolds，作为最后的输出。

进一步地，所述步骤2具体包括以下步骤：

2.1)对于某一contig，设j是contig中的某一位置，j∈[1,L]，L是该contig的长度，用fc(j)表示两端读数都能比对到该contig上的位置j的双端读数的数目，rc(j)表示任意一端读数能比对到该contig上的j位置的双端读数的数目；定义最小覆盖比率fc^*和rc^*作为判断该contig上正确或错误位置的阈值，通过公式(1)、(2)进行计算：

2.2)对于某一contig，设M_c是两端读数都能比对到该contig上的双端读数的数量，M_d是只有一端读数能比对到该contig上的读数数目，定义P_c作为判断测序深度较低的区域的阈值，通过公式(3)进行计算：

P_c＝M_c/(M_C+M_d) (3)

2.3)用P_g代表全基因组的GC含量，对于某一contig，设P_GC为该contig上的GC含量，计算P_g和P_GC；

2.4)对于某一contig，若同时满足fc(j)＜fc^*，rc(j)＜rc^*，P_c＜Min_rate，P_g-0.1≤P_GC≤P_g+0.1，则将其标记为错误类型的节点，同时将该contig上的位置j标记为错误位置；若同时满足fc(j)＜fc^*,rc(j)＜rc^*，P_c＞Min_rate，P_GC＜P_g-0.1或P_GC＞P_g+0.1，则将该contig标记为不确定类型的节点，同时将该contig上的位置j标记为潜在错误位置；将其他所有的contig标记为正确类型的节点；其中Min_rate表示最低比率，为经验参数。

进一步地，所述步骤4包括以下步骤：

4.1)通过为节点分配方向来检测冲突边以及剪切节点；具体步骤为：

4.11)用O_i来表示节点C_i的方向，当取值0表示C_i为正向，当取值1表示C_i为反向；η_ij是一个松弛变量，用来表示节点C_i和C_j之间的边e_ij是否为冲突边，当η_ij为1，则该边不是冲突边，否则是冲突边；用w来表示每轮迭代的边的权值阈值，w的取值范围是[w_min,w_max]，w_max是初始scaffold图中所有边的权值中的最大值，w_min是初始scaffold图中所有边的权值中的最小值；

4.12)构建一个边的集合Ω1，在w的取值范围内设定w的初始值；将scaffold图中权值大于w的边放入集合Ω1中，并把符合要求的边以及与边连接的节点加入到G_s子图中；边的加入要求如下：

对于某条边e_ij，如果这条边不在Ω1中，且O_i≠O_j，那么需要满足约束条件公式(4)；如果这条边不在Ω1中，且O_i＝O_j，那么需要满足约束条件公式(5)；如果这条边在Ω1中，且O_i≠O_j，那么需要满足约束条件公式(4)；如果这条边在Ω1中，且O_i＝O_j，那么需要满足约束条件公式(5)；

在找到满足要求的边的同时，求解优化目标函数，即求使得值最大的η_ij取值，其中，

η_ij≤O_i+O_j≤2-η_ij (4)

η_ij-1≤O_i-O_j≤1-η_ij (5)

其中，η_ij的取值通过求解以下优化函数得到：

其中，w_ij表示边e_ij的权重；

对于Ω1中的边，如果η_ij＝0，那么该边有可能是一条冲突边，接下来则检查该边的两个节点，此时有三种情况：

①若只有一个节点被标记为不确定类型的节点，则按照标记的潜在错误位置剪切该不确定节点(即在该contig标记的潜在错误的位置，将contig断开成两条新的contig)，然后返回步骤2；

②若两个节点都被标记为不确定类型的节点，那么分别计算与两个节点相连的边的权值总和，选择权值总和较小的一个节点进行剪切，然后返回步骤2；

③若两个节点都被标记为错误类型的节点或正确类型的节点，或一个被标记为正确类型的节点另一个被标记为错误类型的节点，说明该边确实存在冲突，不受节点影响，则从scaffold图中删除该边，不对节点进行任何处理；

4.13)令w＝w-0.1，若w＜w_min，则结束迭代，将scaffold图替换为子图G_s，否则返回步骤4.12)；

4.2)通过为节点分配顺序来检测冲突边；具体步骤为：

4.21)用X_i和X_j分别表示分配给C_i和C_j的起始位置，且X_i和X_j均为[0,C]内的整数，且X_i≥X_j，C代表所有节点长度之和的两倍；φ_ij是一个松弛变量，用来表示节点C_i和C_j之间的边e_ij是否为冲突边，当φ_ij为1，则该边不是冲突边，否则是冲突边；用w来表示每轮迭代的边的权值阈值，w的取值范围是[w_min,w_max]，w_max是初始scaffold图中所有边的权值中的最大值，w_min是初始scaffold图中所有边的权值中的最小值；

4.22)构建一个边的集合Ω2，在w的取值范围内设定w的初始值；将scaffold图中权值大于w的边放入集合Ω2中；

4.23)建立以下约束条件和优化函数；

如果边e_ij不在Ω2集合中，则边满足约束条件公式(7)，如果边e_ij在Ω2集合中，则边满足约束条件公式(8)；

0＜X_j-X_i-Len_i≤μ_is+3σ_is (8)

优化函数如公式(9)所示；

通过求解优化函数得到φ_ij的取值，当φ_ij取值为1时，将该边保留在集合Ω2中，当φ_ij取值为0时，将该边从集合Ω2中删除；

4.24)令w＝w-0.1，若w＜w_min，则结束迭代，此时若边不存在于Ω2集合中，那么该边视为冲突边，从scaffold图中删除该边，删除冲突边后的scaffold图即为优化后的scaffold图，否则返回步骤4.22)。

进一步地，所述步骤1中对比对结果进行过滤包括以下步骤：

1.1)对于双端读数中的每个读数，只保留比对得分值最高的比对位置信息；

1.2)判断比对到同一个contig上的双端读数的插入长度是否在[μ_is-3σ_is,μ_is+3σ_is]范围内，若超过该范围则删除该条比对信息；其中μ_is表示比对到同一个contig上的所有双端读数的插入长度的均值，σ_is表示比对到同一个contig上的所有双端读数的插入长度的方差；

1.3)基于读数比对信息计算contig上每个位置的读数覆盖度，以及该contig上所有位置读数覆盖度的平均值μ和标准差σ；若该contig上某个位置的读数覆盖度与μ的差值大于k倍σ，则删除该条比对信息；其中k为经验参数。

进一步地，所述步骤1.3)中，k＝2。

进一步地，所述步骤4.12)和步骤4.22)中，设定w的初始值为0.9。

有益效果：

本发明的方法不仅结合了双端读数的比对信息和contig集合的GC含量来对节点进行分类，并且通过节点与边的结合的方式来对scaffold图进行优化，该方法能够有效解决普遍存在的由于复杂重复区、测序深度不均衡以及测序错误造成的contig错误拼接的问题。

附图说明

图1：本发明流程图

图2：双端读数比对信息图

图3：节点切割示例图

图4：节点方向分配及切割流程示例

具体实施方式

一、数据预处理

读入fastq格式的双端读数文库，以及fasta格式的contig数据集合，将双端读数比对到contig集合上，获得.sam格式的比对结果，在该结果文件中，一条读数通常会有多个比对位置，这是由于测序错误或者序列中的重复区段导致的。结合比对得分值统计、双端读数的插入长度(insert size)、contigs读数覆盖度对序列比对信息进行评估并过滤掉质量较低的部分。

对于双端读数中的每个读数，只保留比对得分值最高的(最优的)比对位置信息，删除其它比对位置信息(删除非最优的比对位置信息)；

由于比对到同一个contig上的双端读数有很多对，每一对双端读数之间有一个距离，即插入长度(insert size)，这些双端读数的插入长度通常符合正态分布N(μ_is,σ_is)，即比对到同一个contig上的双端读数的插入长度应落在[μ_is-3σ_is,μ_is+3σ_is]范围之间，若不满足条件，则删除该条比对信息。然后基于读数比对信息计算contig上每个位置的读数覆盖度，以及该contig上所有位置读数覆盖度的平均值μ和标准差σ；若该contig上某个位置的读数覆盖度与μ的差值大于2σ，则同样删除该条比对信息。为了降低由于测序错误带来的影响，采用移除了低质量读数比对信息后剩下的读数数目除以比对结果中所有能够比对到同一个contig上的读数数目，得到序列的错误率e，测序错误率e将在后续步骤发挥作用。

二、节点质量评估及分类

在本步骤中，采用比对信息以及contig集合的GC含量对contig进行质量评估和分类。

对于某一contig，设j是contig中的某一位置，L是该contig的长度，那么j的范围在[1,L]之间，fc(j)表示两端读数都能比对到该contig上的位置j的双端读数的数目，rc(j)表示任意一端读数能比对到该contig上的j位置的双端读数的数目；定义最小覆盖比率fc^*和rc^*作为判断该contig上正确或错误位置的阈值，可以通过公式(1)、(2)进行计算：

对于测序深度较低的区域，定义M_c以及M_d来对contig质量进行评估，对于某一contig，M_c是两端读数都能比对到该contig上的双端读数的数量，M_d是只有一端读数能比对到该contig上的读数数目，定义P_c作为判断测序深度较低的区域的阈值，可以通过公式(3)进行计算：

P_c＝M_c/(M_C+M_d) (3)

为了减少由于GC含量(在DNA的4种碱基中，鸟嘌呤和胞嘧啶所占的比率称为GC含量)偏差造成的覆盖度不平衡，用P_g代表全基因组的GC含量，P_GC代表该contig上的GC含量(即统计该区域内鸟嘌呤和胞嘧啶所占的比率)，综合以上提及的质量评价指标，Min_rate表示最低比率(根据经验或实验结果人为设置)，可以将contig分为以下三类：正确、错误、不确定。当同时满足fc(j)＜fc^*，rc(j)＜rc^*，P_c＜Min_rate，P_g-0.1≤P_GC≤P_g+0.1，则将该contig标记为错误(False)类型的节点，同时将contig上的位置j标记为错误位置并进行存储；当同时满足fc(j)＜fc^*,rc(j)＜rc^*，P_c＞Min_rate，P_GC＜P_g-0.1或P_GC＞P_g+0.1，则将该contig标记为不确定(Uncertain)类型的节点，同时将该contig上的位置j标记为潜在错误位置并进行存储；将剩下的所有contig标记为正确(True)类型的节点。

三、构建加权的scaffold图

本步骤利用双端读数的比对信息来决定是否在两个节点，即两个contig之间添加边，以及如何给边赋予恰当的权值。对于每两个contig，比对上的双端读数有四种情况，即两条读数比对到两条正向contig，两条读数比对到两条反向contig，两条读数比对到两条方向相对的contig，两条读数比对到两条方向相反的contig。对于前两种情况，两个contig在同一方向上；对于后两种情况，两个contig在相反方向上。本步骤通过统计这四种情况出现的数量，选取数量最多的一种来进行计算，而对于数量较少的其他情况则舍弃。

对于两个contigC_i和C_j，用G_t表示由能够比对到C_i和C_j上的第t对双端读数计算得到的C_i和C_j的之间的空白间距大小，用R₁和R₂表示能够比对到C_i和C_j上的第t对双端读数，并且R₁比对到C_i端，R₂比对到C_j端，用p₁表示R₁比对到C_i上的起始位置，p₂表示R₂比对到C_j的起始位置，r表示读数的长度(本方法采用的是原始双端读数文件，双端读数中两个读数的长度是一致的，统一用r表示)，则G_t可以通过公式(4)进行计算。

G_t＝μ_is-(len_i-p₁)-(p₂+r) (4)

其中，len_i表示第i个contig，即C_i的长度。

由所有能够比对到C_i和C_j上的双端读数计算C_i和C_j之间的空白间距G_ij，可以通过公式(5)进行估算，

其中，n表示所有能够比对到C_i和C_j上的双端读数的数目；

对于双端读数R_t和R_s，如果R_t比对到C_i的位置是p_t，那么R_t的配偶读数R_s能够比对到C_j的概率为插入长度落在区间[D_t+G_ij+r,D_t+G_ij+len_j]内的概率，其中D_t＝len_i-p_t，假设插入长度服从正态分布，那么该概率可以用公式(6)进行计算，f(x)表示正态分布的密度函数。

通常，插入长度的长度小于μ_is+3σ_is，所以可以将计算范围缩小至[max(μ_is+3σ_is-G_ij-r,0),len_i]，所以双端读数比对到C_i和C_j的期望值可以用公式(7)进行计算，s表示区间[max(μ_is+3σ_is-G_ij-r,0),len_i]内比对上的读数的数目，e代表测序错误率(由数据预处理阶段得到)，P_t由公式(6)计算所得。

接下来，计算比对期望值与真实值之间的比率ρ_ij，如果C_i与C_j邻接的可能性越大，那么比对期望值与真实值越接近，ρ_ij的值越小，ρ_ij可以通过公式(8)计算。ρ_ji的计算与ρ_ij相似，不同在于将R₁比对到C_j端，R₂比对到C_i端。此处可以由用户设定一个判定阈值，默认值为0.2，当ρ_ij和ρ_ji的平均值超过设定的阈值，则不添加边，否则在两个节点之间添加一条边。

当对所有节点进行处理后，得到初始scaffold图。初始scaffold图中边的权值即为ρ_ij和ρ_ji的平均值。

四、对scaffold图进行优化

(1)通过为节点分配方向来检测冲突边以及剪切节点。在构建好的初始scaffold图中，如果一条边连接着一个contig的5’端和另一个contig的3’端，那么定义这条边为同向边，反之，定义这条边为反向边。当scaffold图中的其中一个节点确定了方向，根据边的同向或反向信息，就可以推测出所有节点的方向。然而，在scaffold图中往往存在着一些会导致方向冲突的边。

本步骤构建了一个整数线性规划模型LP(OB,CS)，采用迭代策略来检测冲突边，OB是优化函数，CS是约束集合，用O_i来表示节点C_i的方向，当取值0表示C_i为正向，当取值1表示C_i为反向，η_ij是一个松弛变量，用来表示节点C_i和C_j之间的边e_ij是否为冲突边，当η_ij为1，则该边不是冲突边，否则是冲突边；用w来表示每轮迭代的边的权值阈值，w的取值范围是[w_min,w_max]，w_max是初始scaffold图中所有边的权值中的最大值，w_min是初始scaffold图中所有边的权值中的最小值；

初始状态下，构建一个边的集合Ω1，设置w的初始值为0.9，每次迭代都将权值大于w的边放入集合中，并把符合要求的边以及与边连接的节点加入到G_s子图中；边的加入要求如下：

对于某条边e_ij，如果这条边不在Ω1中，且O_i≠O_j，那么需要满足约束条件公式(9)；如果这条边不在Ω1中，且O_i＝O_j，那么需要满足约束条件公式(10)；如果这条边在Ω1中，且O_i≠O_j，那么需要满足约束条件公式(9)；如果这条边在Ω1中，且O_i＝O_j，那么需要满足约束条件公式(10)；在找到满足要求的边的同时，求解优化目标函数，即求使得值最大的η_ij取值，其中，w_ij表示边e_ij的权重。

η_ij≤O_i+O_j≤2-η_ij (9)

η_ij-1≤O_i-O_j≤1-η_ij (10)

优化函数如公式(11)所示。

对于Ω1中的边，如果η_ij＝0，那么该边有可能是一条冲突边，接下来则检查该边的两个节点，此时有三种情况：

①若只有一个节点被标记为不确定类型的节点，则按照标记的潜在错误位置剪切该不确定节点(即在该contig标记的潜在错误的位置，将contig断开成两条新的contig)，然后返回步骤2(重新进行节点质量评估及分类，并重新构建scaffold图后，再次检测是否存在冲突边)；

②若两个节点都被标记为不确定类型的节点，那么分别计算与两个节点相连的边的权值总和，选择权值总和较小的一个节点进行剪切，然后返回步骤2(重新进行节点质量评估及分类，并重新构建scaffold图后，再次检测是否存在冲突边)；

③若两个节点都被标记为错误类型的节点或正确类型的节点，或一个被标记为正确类型的节点另一个被标记为错误类型的节点，说明该边确实存在冲突，不受节点影响，则删除该边，不对节点进行任何处理。

每次迭代后w都减去0.1，并进行下一次迭代，直到w＜w_min，结束迭代；所有迭代结束后，将scaffold图替换为子图G_s，以节省计算机内存资源。

(2)通过为节点分配顺序来检测冲突边。在scaffold图中，边不仅约束着两个节点的方向，还约束着两个节点的顺序关系，本步骤通过为每个节点分配起始位置来确定节点间的顺序关系，并且使得节点间的空白距离(gap)尽可能小。两个contig之间的距离可以通过分配的起始位置计算得到，当计算得到的距离与边约束的距离相差太大，则认为该边有可能是冲突边，并删除该边。X_i代表C_i的起始位置，范围位于[0,C]之间，X_i定义为整数，C代表所有节点长度之和的两倍。用w来表示每轮迭代的边的权值阈值，将w初始化为0.9，每次迭代中，只考虑权值大于w的边，并且每次迭代后w值都减去0.1。Ω2表示初始设置的边集合，每次迭代中，都将权值大于w的边放入集合中，Φ_ij是松弛变量(取值为0或1)，X_i和X_j分别表示分配给C_i和C_j的起始位置，且X_i≥X_j。如果边e_ij不在Ω2集合中，则边满足约束条件公式(12)，如果边e_ij在Ω2集合中，则边满足约束条件公式(13)。

0＜X_j-X_i-Len_i≤μ_is+3σ_is (13)

优化函数如公式(14)所示；通过优化函数求解出Φ_ij的取值，当Φ_ij取值为1时，将该边保留在集合中，当Φ_ij取值为0时，将该边从集合中删除；

多次迭代直到w＜w_min，结束迭代，此时若边不存在于Ω2集合中，那么该边视为冲突边，从scaffold图中删除该边，删除冲突边后的scaffold图即为优化后的scaffold图。

五、从图中提取scaffolds

本步骤采用广度优先遍历算法提取scaffolds，将长度大于μ+3σ的节点定义为长节点，首先提取长节点以及与节点相关联的边作为简单路径，然后将短节点加入到简单路径中。如果两个长节点之间有多个短节点与他们相连，则根据距离进行排序，并且依次将短节点插入到长节点之间。若输入数据集中包含多个双端读数文库，则依次将前一个文库得到的scaffold集合输出作为下一个集合的contig输入。

六、实验分析

为了验证本方法在真实数据上的效果，我们在四组真实数据上进行验证。这四组真实数据是由Illumina技术测序得到，包括金黄色酿脓葡萄球菌(S.aureus)，红假单胞菌(R.sphaeroides)，人类14号染色体(Human 14)和恶性疟原虫(P.falciparum)。前两个物种仅包含1个文库，后两个物种分别包含两个文库。这四个数据集的详细特征信息见表1。

表1数据集详细特征信息

原始的双端读数文件可以从GAGE下载，contig数据集由Velvet工具产生。为了更加全面的对本方法进行评价，实验过程中采用了由Hunt提出的方法来对实验结果进行评价。

表2.S.aureus和R.sphaeroides的评价结果

表3.P.falciparum和Human 14的评价结果

从表2，表3可以看出，本方法与其他12种方法相比，在这四个数据集上能够获取最高的F-score值，本方法在CJ(正确连接数量)这一指标上具有明显优势，除了数据集Human14，本方法在其他三个方法上都获得了最高的CJ值，表明经过本方法的拼接，能够得到更多的正确连接的scaffolds，并且scaffolds的质量相对其他工具更高。

Claims (6)

1.一种基于contig质量评估分类及图优化的scaffolding方法，其特征在于，包括以下步骤：

步骤1、数据预处理：首先将双端读数比对到已有的contig集合上，得到比对结果；然后对比对结果进行过滤；

步骤2、节点质量评估及分类：利用步骤1中保留的比对信息以及contig中的GC含量特征来对contig的质量进行评估，并根据质量评估结果将contig分为正确，错误，不确定三大类；

步骤3、构建加权的scaffold图：将每个contig作为一个节点，根据双端读数比对到contig上的数量期望值以及实际值之间的差异判断是否在两个节点之间构建边，并计算边的权值；

步骤4、对scaffold图进行优化：在比对信息的基础上结合节点的分类与边的权值对scaffold图进行优化；

步骤5、从图中提取scaffolds：采用广度遍历的方法从优化后的scaffold图中尽可能提取较长的scaffolds，作为最后的输出。

2.根据权利要求1所述的基于contig质量评估分类及图优化的scaffolding方法，其特征在于，所述步骤2具体包括以下步骤：

2.1)对于某一contig，设j是contig中的某一位置，j∈[1,L]，L是该contig的长度，用fc(j)表示两端读数都能比对到该contig上的位置j的双端读数的数目，rc(j)表示任意一端读数能比对到该contig上的j位置的双端读数的数目；定义最小覆盖比率fc^*和rc^*作为判断该contig上正确或错误位置的阈值，通过公式(1)、(2)进行计算：

2.2)对于某一contig，设M_c是两端读数都能比对到该contig上的双端读数的数量，M_d是只有一端读数能比对到该contig上的读数数目，定义P_c作为判断测序深度较低的区域的阈值，通过公式(3)进行计算：

P_c＝M_c/(M_C+M_d) (3)

2.3)用P_g代表全基因组的GC含量，对于某一contig，设P_GC为该contig上的GC含量，计算P_g和P_GC；

2.4)对于某一contig，若同时满足fc(j)＜fc^*，rc(j)＜rc^*，P_c＜Min_rate，P_g-0.1≤P_GC≤P_g+0.1，则将其标记为错误类型的节点，同时将该contig上的位置j标记为错误位置；若同时满足fc(j)＜fc^*,rc(j)＜rc^*，P_c＞Min_rate，P_GC＜P_g-0.1或P_GC＞P_g+0.1，则将该contig标记为不确定类型的节点，同时将该contig上的位置j标记为潜在错误位置；将其他所有的contig标记为正确类型的节点；其中Min_rate表示最低比率，为经验参数。

3.根据权利要求1所述的基于contig质量评估分类及图优化的scaffolding方法，其特征在于，所述步骤4包括以下步骤：

4.1)通过为节点分配方向来检测冲突边以及剪切节点；具体步骤为：

4.11)用O_i来表示节点C_i的方向，当取值0表示C_i为正向，当取值1表示C_i为反向；η_ij是一个松弛变量，用来表示节点C_i和C_j之间的边e_ij是否为冲突边，当η_ij为1，则该边不是冲突边，否则是冲突边；用w来表示每轮迭代的边的权值阈值，w的取值范围是[w_min,w_max]，w_max是初始scaffold图中所有边的权值中的最大值，w_min是初始scaffold图中所有边的权值中的最小值；

4.12)构建一个边的集合Ω1，在w的取值范围内设定w的初始值；将scaffold图中权值大于w的边放入集合Ω1中，并把符合要求的边以及与边连接的节点加入到G_s子图中；边的加入要求如下：

对于某条边e_ij，如果这条边不在Ω1中，且O_i≠O_j，那么需要满足约束条件公式(4)；如果这条边不在Ω1中，且O_i＝O_j，那么需要满足约束条件公式(5)；如果这条边在Ω1中，且O_i≠O_j，那么需要满足约束条件公式(4)；如果这条边在Ω1中，且O_i＝O_j，那么需要满足约束条件公式(5)；

在找到满足要求的边的同时，求解优化目标函数，即求使得值最大的η_ij取值，其中，

η_ij≤O_i+O_j≤2-η_ij (4)

η_ij-1≤O_i-O_j≤1-η_ij (5)

其中，η_ij的取值通过求解以下优化函数得到：

其中，w_ij表示边e_ij的权重；

对于Ω1中的边，如果η_ij＝0，那么该边有可能是一条冲突边，接下来则检查该边的两个节点，此时有三种情况：

①若只有一个节点被标记为不确定类型的节点，则按照标记的潜在错误位置剪切该不确定节点(即在该contig标记的潜在错误的位置，将contig断开成两条新的contig)，然后返回步骤2；

②若两个节点都被标记为不确定类型的节点，那么分别计算与两个节点相连的边的权值总和，选择权值总和较小的一个节点进行剪切，然后返回步骤2；

③若两个节点都被标记为错误类型的节点或正确类型的节点，或一个被标记为正确类型的节点另一个被标记为错误类型的节点，说明该边确实存在冲突，不受节点影响，则从scaffold图中删除该边，不对节点进行任何处理；

4.13)令w＝w-0.1，若w＜w_min，则结束迭代，将scaffold图替换为子图G_s，否则返回步骤4.12)；

4.2)通过为节点分配顺序来检测冲突边；具体步骤为：

4.21)用X_i和X_j分别表示分配给C_i和C_j的起始位置，且X_i和X_j均为[0,C]内的整数，且X_i≥X_j，C代表所有节点长度之和的两倍；Φ_ij是一个松弛变量，用来表示节点C_i和C_j之间的边e_ij是否为冲突边，当Φ_ij为1，则该边不是冲突边，否则是冲突边；用w来表示每轮迭代的边的权值阈值，w的取值范围是[w_min,w_max]，w_max是初始scaffold图中所有边的权值中的最大值，w_min是初始scaffold图中所有边的权值中的最小值；

4.22)构建一个边的集合Ω2，在w的取值范围内设定w的初始值；将scaffold图中权值大于w的边放入集合Ω2中；

4.23)建立以下约束条件和优化函数；

如果边e_ij不在Ω2集合中，则边满足约束条件公式(7)，如果边e_ij在Ω2集合中，则边满足约束条件公式(8)；

0＜X_j-X_i-Len_i≤μ_is+3σ_is (8)

优化函数如公式(9)所示；

通过求解优化函数得到Φ_ij的取值，当φ_ij取值为1时，将该边保留在集合Ω2中，当φ_ij取值为0时，将该边从集合Ω2中删除；

4.24)令w＝w-0.1，若w＜w_min，则结束迭代，此时若边不存在于Ω2集合中，那么该边视为冲突边，从scaffold图中删除该边，删除冲突边后的scaffold图即为优化后的scaffold图，否则返回步骤4.22)。

4.根据权利要求1～3中任一项所述的基于contig质量评估分类及图优化的scaffolding方法，其特征在于，所述步骤1中对比对结果进行过滤包括以下步骤：

1.1)对于双端读数中的每个读数，只保留比对得分值最高的比对位置信息；

1.2)判断比对到同一个contig上的双端读数的插入长度是否在[μ_is-3σ_is,μ_is+3σ_is]范围内，若超过该范围则删除该条比对信息；其中μ_is表示比对到同一个contig上的所有双端读数的插入长度的均值，σ_is表示比对到同一个contig上的所有双端读数的插入长度的方差；

1.3)基于读数比对信息计算contig上每个位置的读数覆盖度，以及该contig上所有位置读数覆盖度的平均值μ和标准差σ；若该contig上某个位置的读数覆盖度与μ的差值大于k倍σ，则删除该条比对信息；其中k为经验参数。

5.根据权利要求1～4中任一项所述的基于contig质量评估分类及图优化的scaffolding方法，其特征在于，所述步骤1.3)中，k＝2。

6.根据权利要求5所述的基于contig质量评估分类及图优化的scaffolding方法，其特征在于，所述步骤4.12)和步骤4.22)中，设定w的初始值为0.9。