CN108491687B - 基于contig质量评估分类及图优化的scaffolding方法 - Google Patents
基于contig质量评估分类及图优化的scaffolding方法 Download PDFInfo
- Publication number
- CN108491687B CN108491687B CN201810242418.1A CN201810242418A CN108491687B CN 108491687 B CN108491687 B CN 108491687B CN 201810242418 A CN201810242418 A CN 201810242418A CN 108491687 B CN108491687 B CN 108491687B
- Authority
- CN
- China
- Prior art keywords
- edge
- contig
- scaffold
- nodes
- edges
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Medical Informatics (AREA)
- General Health & Medical Sciences (AREA)
- Theoretical Computer Science (AREA)
- Biophysics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Bioethics (AREA)
- Data Mining & Analysis (AREA)
- Chemical & Material Sciences (AREA)
- Molecular Biology (AREA)
- Genetics & Genomics (AREA)
- Artificial Intelligence (AREA)
- Analytical Chemistry (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Databases & Information Systems (AREA)
- Epidemiology (AREA)
- Evolutionary Computation (AREA)
- Public Health (AREA)
- Software Systems (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于contig质量评估分类及图优化的scaffolding方法,采用序列比对信息以及contig的GC含量信息对contig集合进行质量评估并分类,再将每个contig作为一个节点,根据双端读数比对到contig上的数量期望值以及实际值之间的差异判断是否在两个节点之间构建边,并计算边的权值,构建加权的scaffold图。最后通过为节点分配方向以及剪切节点来消除scaffold图中的方向冲突,通过为节点分配顺序来消除scaffold图中的顺序冲突。本发明简单易用,在四组真实测序数据上表现出良好的拼接结果,较其他序列拼接方法具有更高的F‑score值。
Description
技术领域
本发明属于生物信息学领域,涉及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,设Mc是两端读数都能比对到该contig上的双端读数的数量,Md是只有一端读数能比对到该contig上的读数数目,定义Pc作为判断测序深度较低的区域的阈值,通过公式(3)进行计算:
Pc=Mc/(MC+Md) (3)
2.3)用Pg代表全基因组的GC含量,对于某一contig,设PGC为该contig上的GC含量,计算Pg和PGC;
2.4)对于某一contig,若同时满足fc(j)<fc*,rc(j)<rc*,Pc<Min_rate,Pg-0.1≤PGC≤Pg+0.1,则将其标记为错误类型的节点,同时将该contig上的位置j标记为错误位置;若同时满足fc(j)<fc*,rc(j)<rc*,Pc>Min_rate,PGC<Pg-0.1或PGC>Pg+0.1,则将该contig标记为不确定类型的节点,同时将该contig上的位置j标记为潜在错误位置;将其他所有的contig标记为正确类型的节点;其中Min_rate表示最低比率,为经验参数。
进一步地,所述步骤4包括以下步骤:
4.1)通过为节点分配方向来检测冲突边以及剪切节点;具体步骤为:
4.11)用Oi来表示节点Ci的方向,当取值0表示Ci为正向,当取值1表示Ci为反向;ηij是一个松弛变量,用来表示节点Ci和Cj之间的边eij是否为冲突边,当ηij为1,则该边不是冲突边,否则是冲突边;用w来表示每轮迭代的边的权值阈值,w的取值范围是[wmin,wmax],wmax是初始scaffold图中所有边的权值中的最大值,wmin是初始scaffold图中所有边的权值中的最小值;
对于某条边eij,如果这条边不在Ω1中,且Oi≠Oj,那么需要满足约束条件公式(4);如果这条边不在Ω1中,且Oi=Oj,那么需要满足约束条件公式(5);如果这条边在Ω1中,且Oi≠Oj,那么需要满足约束条件公式(4);如果这条边在Ω1中,且Oi=Oj,那么需要满足约束条件公式(5);
ηij≤Oi+Oj≤2-ηij (4)
ηij-1≤Oi-Oj≤1-ηij (5)
其中,ηij的取值通过求解以下优化函数得到:
其中,wij表示边eij的权重;
对于Ω1中的边,如果ηij=0,那么该边有可能是一条冲突边,接下来则检查该边的两个节点,此时有三种情况:
①若只有一个节点被标记为不确定类型的节点,则按照标记的潜在错误位置剪切该不确定节点(即在该contig标记的潜在错误的位置,将contig断开成两条新的contig),然后返回步骤2;
②若两个节点都被标记为不确定类型的节点,那么分别计算与两个节点相连的边的权值总和,选择权值总和较小的一个节点进行剪切,然后返回步骤2;
③若两个节点都被标记为错误类型的节点或正确类型的节点,或一个被标记为正确类型的节点另一个被标记为错误类型的节点,说明该边确实存在冲突,不受节点影响,则从scaffold图中删除该边,不对节点进行任何处理;
4.13)令w=w-0.1,若w<wmin,则结束迭代,将scaffold图替换为子图Gs,否则返回步骤4.12);
4.2)通过为节点分配顺序来检测冲突边;具体步骤为:
4.21)用Xi和Xj分别表示分配给Ci和Cj的起始位置,且Xi和Xj均为[0,C]内的整数,且Xi≥Xj,C代表所有节点长度之和的两倍;φij是一个松弛变量,用来表示节点Ci和Cj之间的边eij是否为冲突边,当φij为1,则该边不是冲突边,否则是冲突边;用w来表示每轮迭代的边的权值阈值,w的取值范围是[wmin,wmax],wmax是初始scaffold图中所有边的权值中的最大值,wmin是初始scaffold图中所有边的权值中的最小值;
4.23)建立以下约束条件和优化函数;
如果边eij不在Ω2集合中,则边满足约束条件公式(7),如果边eij在Ω2集合中,则边满足约束条件公式(8);
0<Xj-Xi-Leni≤μis+3σis (8)
优化函数如公式(9)所示;
通过求解优化函数得到φij的取值,当φij取值为1时,将该边保留在集合Ω2中,当φij取值为0时,将该边从集合Ω2中删除;
4.24)令w=w-0.1,若w<wmin,则结束迭代,此时若边不存在于Ω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)进行计算:
对于测序深度较低的区域,定义Mc以及Md来对contig质量进行评估,对于某一contig,Mc是两端读数都能比对到该contig上的双端读数的数量,Md是只有一端读数能比对到该contig上的读数数目,定义Pc作为判断测序深度较低的区域的阈值,可以通过公式(3)进行计算:
Pc=Mc/(MC+Md) (3)
为了减少由于GC含量(在DNA的4种碱基中,鸟嘌呤和胞嘧啶所占的比率称为GC含量)偏差造成的覆盖度不平衡,用Pg代表全基因组的GC含量,PGC代表该contig上的GC含量(即统计该区域内鸟嘌呤和胞嘧啶所占的比率),综合以上提及的质量评价指标,Min_rate表示最低比率(根据经验或实验结果人为设置),可以将contig分为以下三类:正确、错误、不确定。当同时满足fc(j)<fc*,rc(j)<rc*,Pc<Min_rate,Pg-0.1≤PGC≤Pg+0.1,则将该contig标记为错误(False)类型的节点,同时将contig上的位置j标记为错误位置并进行存储;当同时满足fc(j)<fc*,rc(j)<rc*,Pc>Min_rate,PGC<Pg-0.1或PGC>Pg+0.1,则将该contig标记为不确定(Uncertain)类型的节点,同时将该contig上的位置j标记为潜在错误位置并进行存储;将剩下的所有contig标记为正确(True)类型的节点。
三、构建加权的scaffold图
本步骤利用双端读数的比对信息来决定是否在两个节点,即两个contig之间添加边,以及如何给边赋予恰当的权值。对于每两个contig,比对上的双端读数有四种情况,即两条读数比对到两条正向contig,两条读数比对到两条反向contig,两条读数比对到两条方向相对的contig,两条读数比对到两条方向相反的contig。对于前两种情况,两个contig在同一方向上;对于后两种情况,两个contig在相反方向上。本步骤通过统计这四种情况出现的数量,选取数量最多的一种来进行计算,而对于数量较少的其他情况则舍弃。
对于两个contigCi和Cj,用Gt表示由能够比对到Ci和Cj上的第t对双端读数计算得到的Ci和Cj的之间的空白间距大小,用R1和R2表示能够比对到Ci和Cj上的第t对双端读数,并且R1比对到Ci端,R2比对到Cj端,用p1表示R1比对到Ci上的起始位置,p2表示R2比对到Cj的起始位置,r表示读数的长度(本方法采用的是原始双端读数文件,双端读数中两个读数的长度是一致的,统一用r表示),则Gt可以通过公式(4)进行计算。
Gt=μis-(leni-p1)-(p2+r) (4)
其中,leni表示第i个contig,即Ci的长度。
由所有能够比对到Ci和Cj上的双端读数计算Ci和Cj之间的空白间距Gij,可以通过公式(5)进行估算,
其中,n表示所有能够比对到Ci和Cj上的双端读数的数目;
对于双端读数Rt和Rs,如果Rt比对到Ci的位置是pt,那么Rt的配偶读数Rs能够比对到Cj的概率为插入长度落在区间[Dt+Gij+r,Dt+Gij+lenj]内的概率,其中Dt=leni-pt,假设插入长度服从正态分布,那么该概率可以用公式(6)进行计算,f(x)表示正态分布的密度函数。
通常,插入长度的长度小于μis+3σis,所以可以将计算范围缩小至[max(μis+3σis-Gij-r,0),leni],所以双端读数比对到Ci和Cj的期望值可以用公式(7)进行计算,s表示区间[max(μis+3σis-Gij-r,0),leni]内比对上的读数的数目,e代表测序错误率(由数据预处理阶段得到),Pt由公式(6)计算所得。
接下来,计算比对期望值与真实值之间的比率ρij,如果Ci与Cj邻接的可能性越大,那么比对期望值与真实值越接近,ρij的值越小,ρij可以通过公式(8)计算。ρji的计算与ρij相似,不同在于将R1比对到Cj端,R2比对到Ci端。此处可以由用户设定一个判定阈值,默认值为0.2,当ρij和ρji的平均值超过设定的阈值,则不添加边,否则在两个节点之间添加一条边。
当对所有节点进行处理后,得到初始scaffold图。初始scaffold图中边的权值即为ρij和ρji的平均值。
四、对scaffold图进行优化
(1)通过为节点分配方向来检测冲突边以及剪切节点。在构建好的初始scaffold图中,如果一条边连接着一个contig的5’端和另一个contig的3’端,那么定义这条边为同向边,反之,定义这条边为反向边。当scaffold图中的其中一个节点确定了方向,根据边的同向或反向信息,就可以推测出所有节点的方向。然而,在scaffold图中往往存在着一些会导致方向冲突的边。
本步骤构建了一个整数线性规划模型LP(OB,CS),采用迭代策略来检测冲突边,OB是优化函数,CS是约束集合,用Oi来表示节点Ci的方向,当取值0表示Ci为正向,当取值1表示Ci为反向,ηij是一个松弛变量,用来表示节点Ci和Cj之间的边eij是否为冲突边,当ηij为1,则该边不是冲突边,否则是冲突边;用w来表示每轮迭代的边的权值阈值,w的取值范围是[wmin,wmax],wmax是初始scaffold图中所有边的权值中的最大值,wmin是初始scaffold图中所有边的权值中的最小值;
对于某条边eij,如果这条边不在Ω1中,且Oi≠Oj,那么需要满足约束条件公式(9);如果这条边不在Ω1中,且Oi=Oj,那么需要满足约束条件公式(10);如果这条边在Ω1中,且Oi≠Oj,那么需要满足约束条件公式(9);如果这条边在Ω1中,且Oi=Oj,那么需要满足约束条件公式(10);在找到满足要求的边的同时,求解优化目标函数,即求使得值最大的ηij取值,其中,wij表示边eij的权重。
ηij≤Oi+Oj≤2-ηij (9)
ηij-1≤Oi-Oj≤1-ηij (10)
优化函数如公式(11)所示。
对于Ω1中的边,如果ηij=0,那么该边有可能是一条冲突边,接下来则检查该边的两个节点,此时有三种情况:
①若只有一个节点被标记为不确定类型的节点,则按照标记的潜在错误位置剪切该不确定节点(即在该contig标记的潜在错误的位置,将contig断开成两条新的contig),然后返回步骤2(重新进行节点质量评估及分类,并重新构建scaffold图后,再次检测是否存在冲突边);
②若两个节点都被标记为不确定类型的节点,那么分别计算与两个节点相连的边的权值总和,选择权值总和较小的一个节点进行剪切,然后返回步骤2(重新进行节点质量评估及分类,并重新构建scaffold图后,再次检测是否存在冲突边);
③若两个节点都被标记为错误类型的节点或正确类型的节点,或一个被标记为正确类型的节点另一个被标记为错误类型的节点,说明该边确实存在冲突,不受节点影响,则删除该边,不对节点进行任何处理。
每次迭代后w都减去0.1,并进行下一次迭代,直到w<wmin,结束迭代;所有迭代结束后,将scaffold图替换为子图Gs,以节省计算机内存资源。
(2)通过为节点分配顺序来检测冲突边。在scaffold图中,边不仅约束着两个节点的方向,还约束着两个节点的顺序关系,本步骤通过为每个节点分配起始位置来确定节点间的顺序关系,并且使得节点间的空白距离(gap)尽可能小。两个contig之间的距离可以通过分配的起始位置计算得到,当计算得到的距离与边约束的距离相差太大,则认为该边有可能是冲突边,并删除该边。Xi代表Ci的起始位置,范围位于[0,C]之间,Xi定义为整数,C代表所有节点长度之和的两倍。用w来表示每轮迭代的边的权值阈值,将w初始化为0.9,每次迭代中,只考虑权值大于w的边,并且每次迭代后w值都减去0.1。Ω2表示初始设置的边集合,每次迭代中,都将权值大于w的边放入集合中,Φij是松弛变量(取值为0或1),Xi和Xj分别表示分配给Ci和Cj的起始位置,且Xi≥Xj。如果边eij不在Ω2集合中,则边满足约束条件公式(12),如果边eij在Ω2集合中,则边满足约束条件公式(13)。
0<Xj-Xi-Leni≤μis+3σis (13)
优化函数如公式(14)所示;通过优化函数求解出Φij的取值,当Φij取值为1时,将该边保留在集合中,当Φij取值为0时,将该边从集合中删除;
多次迭代直到w<wmin,结束迭代,此时若边不存在于Ω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 (5)
1.一种基于contig质量评估分类及图优化的scaffolding方法,其特征在于,包括以下步骤:
步骤1、数据预处理:首先将双端读数比对到已有的contig集合上,得到比对结果;然后对比对结果进行过滤;
步骤2、节点质量评估及分类:利用步骤1中保留的比对信息以及contig中的GC含量特征来对contig的质量进行评估,并根据质量评估结果将contig分为正确,错误,不确定三大类;
步骤3、构建加权的scaffold图:将每个contig作为一个节点,根据双端读数比对到contig上的数量期望值以及实际值之间的差异判断是否在两个节点之间构建边,并计算边的权值;
步骤4、对scaffold图进行优化:在比对信息的基础上结合节点的分类与边的权值对scaffold图进行优化;
所述步骤4包括以下步骤:
4.1)通过为节点分配方向来检测冲突边以及剪切节点;具体步骤为:
4.11)用Oi来表示节点Ci的方向,当取值0表示Ci为正向,当取值1表示Ci为反向;ηij是一个松弛变量,用来表示节点Ci和Cj之间的边eij是否为冲突边,当ηij为1,则该边不是冲突边,否则是冲突边;用w来表示每轮迭代的边的权值阈值,w的取值范围是[wmin,wmax],wmax是初始scaffold图中所有边的权值中的最大值,wmin是初始scaffold图中所有边的权值中的最小值;
对于某条边eij,如果这条边不在Ω1中,且Oi≠Oj,那么需要满足约束条件公式(4);如果这条边不在Ω1中,且Oi=Oj,那么需要满足约束条件公式(5);如果这条边在Ω1中,且Oi≠Oj,那么需要满足约束条件公式(4);如果这条边在Ω1中,且Oi=Oj,那么需要满足约束条件公式(5);
ηij≤Oi+Oj≤2-ηij (4)
ηij-1≤Oi-Oj≤1-ηij (5)
其中,ηij的取值通过求解以下优化函数得到:
其中,wij表示边eij的权重;
对于Ω1中的边,如果ηij=0,那么该边有可能是一条冲突边,接下来则检查该边的两个节点,此时有三种情况:
①若只有一个节点被标记为不确定类型的节点,则按照标记的潜在错误位置剪切该不确定节点,即在该contig标记的潜在错误的位置,将contig断开成两条新的contig,然后返回步骤2;
②若两个节点都被标记为不确定类型的节点,那么分别计算与两个节点相连的边的权值总和,选择权值总和较小的一个节点进行剪切,然后返回步骤2;
③若两个节点都被标记为错误类型的节点或正确类型的节点,或一个被标记为正确类型的节点另一个被标记为错误类型的节点,说明该边确实存在冲突,不受节点影响,则从scaffold图中删除该边,不对节点进行任何处理;
4.13)令w=w-0.1,若w<wmin,则结束迭代,将scaffold图替换为子图Gs,否则返回步骤4.12);
4.2)通过为节点分配顺序来检测冲突边;具体步骤为:
4.21)用Xi和Xj分别表示分配给Ci和Cj的起始位置,且Xi和Xj均为[0,C]内的整数,且Xi≥Xj,C代表所有节点长度之和的两倍;Φij是一个松弛变量,用来表示节点Ci和Cj之间的边eij是否为冲突边,当Φij为1,则该边不是冲突边,否则是冲突边;用w来表示每轮迭代的边的权值阈值,w的取值范围是[wmin,wmax],wmax是初始scaffold图中所有边的权值中的最大值,wmin是初始scaffold图中所有边的权值中的最小值;
4.23)建立以下约束条件和优化函数;
如果边eij不在Ω2集合中,则边满足约束条件公式(7),如果边eij在Ω2集合中,则边满足约束条件公式(8);
0<Xj-Xi-Leni≤μis+3σis (8)
优化函数如公式(9)所示;
通过求解优化函数得到Φij的取值,当φij取值为1时,将该边保留在集合Ω2中,当Φij取值为0时,将该边从集合Ω2中删除;
4.24)令w=w-0.1,若w<wmin,则结束迭代,此时若边不存在于Ω2集合中,那么该边视为冲突边,从scaffold图中删除该边,删除冲突边后的scaffold图即为优化后的scaffold图,否则返回步骤4.22);
步骤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,设Mc是两端读数都能比对到该contig上的双端读数的数量,Md是只有一端读数能比对到该contig上的读数数目,定义Pc作为判断测序深度较低的区域的阈值,通过公式(3)进行计算:
Pc=Mc/(MC+Md) (3)
2.3)用Pg代表全基因组的GC含量,对于某一contig,设PGC为该contig上的GC含量,计算Pg和PGC;
2.4)对于某一contig,若同时满足fc(j)<fc*,rc(j)<rc*,Pc<Min_rate,Pg-0.1≤PGC≤Pg+0.1,则将其标记为错误类型的节点,同时将该contig上的位置j标记为错误位置;若同时满足fc(j)<fc*,rc(j)<rc*,Pc>Min_rate,PGC<Pg-0.1或PGC>Pg+0.1,则将该contig标记为不确定类型的节点,同时将该contig上的位置j标记为潜在错误位置;将其他所有的contig标记为正确类型的节点;其中Min_rate表示最低比率,为经验参数。
3.根据权利要求1~2中任一项所述的基于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为经验参数。
4.根据权利要求3所述的基于contig质量评估分类及图优化的scaffolding方法,其特征在于,所述步骤1.3)中,k=2。
5.根据权利要求1所述的基于contig质量评估分类及图优化的scaffolding方法,其特征在于,所述步骤4.12)和步骤4.22)中,设定w的初始值为0.9。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810242418.1A CN108491687B (zh) | 2018-03-22 | 2018-03-22 | 基于contig质量评估分类及图优化的scaffolding方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810242418.1A CN108491687B (zh) | 2018-03-22 | 2018-03-22 | 基于contig质量评估分类及图优化的scaffolding方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108491687A CN108491687A (zh) | 2018-09-04 |
CN108491687B true CN108491687B (zh) | 2021-07-13 |
Family
ID=63319459
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810242418.1A Active CN108491687B (zh) | 2018-03-22 | 2018-03-22 | 基于contig质量评估分类及图优化的scaffolding方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108491687B (zh) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110066862B (zh) * | 2019-05-22 | 2021-02-12 | 中南大学 | 一种基于高通量测序读数的重复dna序列识别方法 |
CN110544510B (zh) * | 2019-05-31 | 2023-03-24 | 中南大学 | 基于邻接代数模型及质量等级评估的contig集成方法 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104200133A (zh) * | 2014-09-19 | 2014-12-10 | 中南大学 | 一种基于读数和距离分布的基因组De novo序列拼接方法 |
CN106022003A (zh) * | 2016-05-17 | 2016-10-12 | 杭州和壹基因科技有限公司 | 一种基于三代PacBio测序数据的scaffold构建方法 |
CN106355000A (zh) * | 2016-08-25 | 2017-01-25 | 中南大学 | 基于双端读数insert size统计特征的scaffolding方法 |
WO2017214461A1 (en) * | 2016-06-08 | 2017-12-14 | The Broad Institute, Inc. | Linear genome assembly from three dimensional genome structure |
-
2018
- 2018-03-22 CN CN201810242418.1A patent/CN108491687B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104200133A (zh) * | 2014-09-19 | 2014-12-10 | 中南大学 | 一种基于读数和距离分布的基因组De novo序列拼接方法 |
CN106022003A (zh) * | 2016-05-17 | 2016-10-12 | 杭州和壹基因科技有限公司 | 一种基于三代PacBio测序数据的scaffold构建方法 |
WO2017214461A1 (en) * | 2016-06-08 | 2017-12-14 | The Broad Institute, Inc. | Linear genome assembly from three dimensional genome structure |
CN106355000A (zh) * | 2016-08-25 | 2017-01-25 | 中南大学 | 基于双端读数insert size统计特征的scaffolding方法 |
Non-Patent Citations (3)
Title |
---|
MEC: Misassembly error correction in contigs using a combination of paired-end reads and GC-contents;Binbin Wu et al.;《2017 IEEE International Conference on Bioinformatics and Biomedicine (BIBM)》;20171218;第216-219页 * |
基于2b-RAD技术的辅助基因组组装和标记分型研究;窦锦壮;《中国博士学位论文全文数据库 基础科学辑》;20161015(第10期);A006-17 * |
高效的分布式大规模基因组序列组装;徐魁;《中国优秀硕士学位论文全文数据库 基础科学辑》;20170215(第02期);A006-792 * |
Also Published As
Publication number | Publication date |
---|---|
CN108491687A (zh) | 2018-09-04 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Di Genova et al. | Efficient hybrid de novo assembly of human genomes with WENGAN | |
US20130131995A1 (en) | System And Method to Correct Out of Phase Errors In DNA Sequencing Data by Use of a Recursive Algorithm | |
WO2018218788A1 (zh) | 一种基于全局种子打分优选的三代测序序列比对方法 | |
CN108491687B (zh) | 基于contig质量评估分类及图优化的scaffolding方法 | |
CN108660200B (zh) | 一种检测短串联重复序列扩张的方法 | |
CN112365922B (zh) | 用于检测msi的微卫星位点、其筛选方法及应用 | |
US20150178446A1 (en) | Iterative clustering of sequence reads for error correction | |
CN111627501A (zh) | 检测msi的微卫星位点、其筛选方法及应用 | |
KR20190082854A (ko) | 데이터 판독 재정렬을 시퀀싱하는 방법 | |
Goltsman et al. | Meraculous-2D: Haplotype-sensitive assembly of highly heterozygous genomes | |
CN106355000B (zh) | 基于双端读数insert size统计特征的scaffolding方法 | |
US11250931B2 (en) | Systems and methods for detecting recombination | |
WO2023124779A1 (zh) | 基于三代测序数据检测点突变的分析方法和装置 | |
EP3663890B1 (en) | Alignment method, device and system | |
CN112669902B (zh) | 检测基因组结构变异的方法、计算设备和存储介质 | |
CN114896557A (zh) | 融合接近度和度的网络节点重要性评估方法及其应用 | |
CN110544510B (zh) | 基于邻接代数模型及质量等级评估的contig集成方法 | |
CN108830047A (zh) | 一种基于长读数和contig分类的scaffolding方法 | |
CN109637586B (zh) | 测序深度的矫正方法及装置 | |
CN115602244B (zh) | 一种基于序列比对骨架的基因组变异检测方法 | |
CN112599251B (zh) | 疾病筛查模型的构建方法、疾病筛查模型及筛查装置 | |
CN117497056B (zh) | 一种无对照hrd检测方法、系统及装置 | |
CN109584959B (zh) | 测序深度的矫正方法及装置 | |
CN114694755B (zh) | 基因组组装方法、装置、设备及存储介质 | |
CN109637585B (zh) | 测序深度的矫正方法及装置 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
TR01 | Transfer of patent right | ||
TR01 | Transfer of patent right |
Effective date of registration: 20221107 Address after: 518000 a1002 Tian'an innovation and Technology Plaza, chegong temple, Shatou street, Futian District, Shenzhen City, Guangdong Province Patentee after: SHENZHEN ZAOZHIDAO TECHNOLOGY CO.,LTD. Address before: Yuelu District City, Hunan province 410083 Changsha Lushan Road No. 932 Patentee before: CENTRAL SOUTH University |