CN110299185B - 一种基于新一代测序数据的插入变异检测方法及系统 - Google Patents

Abstract

本发明属于基因组测序技术领域，公开了一种基于新一代测序数据的插入变异检测方法；在确定变异发生位点时，有插入变异发生的区域一定会产生分裂读段，针对新序列插入、序列串联倍增、序列散在倍增等插入变异类型及缺失变异、倒置变异的分裂读段分布不同的特性；在确定插入变异发生种类及位点之后，通过利用部分匹配、完全匹配、以及未匹配的读段信息来构造一条虚拟参考序列，与原始参考序列比较得到插入序列的相关信息；利用拷贝数状态信息获得变异基因型。本发明能够解决插入变异位点判定不准确的问题；能够解决SR方法检测插入变异造成遗漏的问题；能够解决现有技术遇到重复序列可能会检测出错的问题。

Description

一种基于新一代测序数据的插入变异检测方法及系统

技术领域

本发明属于基因组测序技术领域，尤其涉及一种基于新一代测序数据的插入变异检测方法。

背景技术

目前，最接近的现有技术：基于新一代测序技术的split read分析方法。新一代测序是一种DNA测序技术，在测序过程中，将完整的样本DNA序列打碎，从中筛选出满足特定长度(通常为数百bp)的片段，在每个片段的一端或两端各读取一段长度为数十至数百bp的序列。读取出的序列长度通常远远小于被测样本DNA序列的长度，但是新一代测序技术可以同时读取大量这样的短序列，使得全部短序列的总长度达到样本DNA长度的数倍至数十倍，使获得样本DNA序列成为可能。插入变异是基因组中的一种重要的变异现象，是人类基因组结构变异的一种形式，并且与人类的疾病发生密切相关。

目前主要有4种检测基因组上插入变异的策略，分别为：(1)Read pair(也称为Pair-end Mapping，简称PEM，双端映射)；(2)split read(简称SR，分裂读段)；(3)ReadDepth(简称RD，读段覆盖深度)；(4)de novo Assembly(简称AS，序列从头组装)。

PEM分析方法：Pair-End(PE)测序的两条read(通常称为Read1和Read2)来自于同一个序列片段，因此，Read1和Read2之间存在着客观的物理关联，而Read1与Read2之间的距离，称为插入片段长度(insert size)。对于PEM方法，插入片段长度的分布是进行变异检测的一个关键信息，将样本Pair-end读段比对到参考序列上，其插入片段长度一般符合正态分布。因此，若某一对读段插入片段长度有异常，组成Read1和Read2的这个序列片段和参考基因组相比就存在着对应序列上变异的可能，比如，若某一对读段的映射长度小于样本平均插入长度，则在其对应的序列上可能存在插入变异。但是，PEM方法所拥有的缺陷在于，对于插入变异序列的检测，受限于Read1与Read2之间的客观物理关联，其检测长度无法超过插入片段的长度，另外，对于小于插入长度的插入变异，其检测精度也受限于插入片段长度的标准差。

SR分析方法：splitread是一类特殊的read，其出现通常是由基因组中的结构变异造成的，这类read在映射中不再保持连续序列的形式，而是包含了一定长度的不匹配部分，因此具有较高的映射难度。SR分析方法首先提取具有以下特点的pair-end读段，一条可以正常比对到参考序列上，另外一条不能比对，这里的不能比对指的是在这一条read上，只有其中一部分可以匹配到参考序列的某个位置，而另外一部分不能正确匹配，或者匹配的位置与前一部分的匹配位置不连续。提取出这些特点的pair-end读段后，利用正常比对的读段位置和插入长度确定一个查找范围，在这个范围内寻找未比对上的读段与参考序列的最佳匹配，通过最佳匹配点把未匹配的读段分割成多段，确定插入变异发生的位置。SR分析方法的缺陷在于单纯利用splitread信息只能初步确定可能发生变异的位点，现有方法大多利用疑似变异发生位点进行序列局部组装，利用split read信息局部组装contig，但是对于拥有重复片段的插入变异区域，即变异来源于基因内部片段的区域，组装contig时会遇到组装终止异常问题，造成组装无法终止或提前终止，导致插入片段的错误检测。

RD分析方法：主要利用覆盖度的信息进行分析，是目前解决基因组拷贝数变异检测(Copy number variantion，简称CNV)的主要方法，该方法通过检测样本在参考基因组上read的深度分布情况来检测序列重复。该方法缺陷在于只能确定较大片段(>1000bp)的重复序列变异，且其只能检测出变异发生的大致区域，无法使变异边界值检测达到bp级别。

AS分析方法：de novo Assembly的方法理论上能够提供对于较长插入片段的最好检测，但是目前的局限却也发生在Assembly本身，其对于人体基因组的全局拼装的时间复杂度呈现指数级增长，对于全局拼装的大规模应用是一大阻碍，且基因组上所存在的重复性序列会严重影响组装的质量，也在很大程度上阻碍了利用组装的方法在基因组变异检测方面的应用。

近年来，已经开发了许多用于检测插入变异的计算方法，除了使用不同的模型之外，通常提取不同的比对信息来推断变异发生的位置和内容，使用SR信息分析的方法主要有Pindel、PRISM、IMSindel、Pamir等。

Pindel使用基于字符串匹配的模式增长算法来从长度为36bp的短读段中识别变异发生的位点，通过分析测序序列对的映射结果，提出一种带锚的测序序列分裂比对方法。其以测序序列对中能够正常映射到参考基因组上的序列末端为锚点，在距离锚点的一定区域内完成分裂比对，其缺点是只能检测长度小于20bp的插入变异。

PRISM利用不一致的配对末端簇来缩小要搜索的变异发生区间并且指导分裂读段来进行映射，这种方法可以降低计算复杂度并且准确度有所提升，但是其仅限于检测小的插入变异。

IMSindel主要利用split read以及未匹配的读段信息来检测插入变异，其核心思想是利用这些读段来构建共识片段，然后利用共识片段构建共识序列，将共识序列与参考序列对比来检测插入变异。同样，该方法对于较大片段的插入变异也无法检测，当读段长度为100bp时，其所能检测的最大插入变异长度不超过250bp。

Pamir方法用于检测任意长度插入变异的插入位点及插入内容。其核心思想是首先将孤立读段组装成孤立contigs，根据OEA读段对的比对位置对其进行聚类，将OEA读段重新比对到孤立contigs以将OEA簇匹配到孤立contigs中，从而定位insertion位点及其内容。Pamir方法主要还是依赖于本地组装来完成插入变异的检测，不能排除重复序列带来的检测误差。同样利用本地组装的方法还有SvABA、MindTheGap等方法。

综上所述，现有技术存在的问题是：

(1)现有技术大多只检测单一种类的插入变异，不符合癌症样本中插入变异类型多样的状况，极大限制了癌症诊断及靶向药物选定的能力。

(2)现有技术大多对大片度插入变异的检测能力不足，侧重于小片段插入及缺失变异的检测，使得方法对于插入变异的检测不够充分，难以全面提炼癌症样本的DNA变异信息。

(3)现有技术的部分方法使用局部从头组装算法检测大片段插入变异，但易受重复序列区域影响造成组装出错，以致得到错误的变异检测结果。

解决上述技术问题的难度：

(1)因为新一代测序技术获取DNA片段的客观物理问题，其所获取的read读段通常为100-250bp，所以利用这些read信息来对基因插入变异进行检测时，对于大片段(50-1000bp)插入变异的检测不能通过简单的read比对来获得，使其检测模型的建立变得尤为困难，所以大片段插入变异的检测是一大挑战。

(2)基因中重复序列较为常见，因为新一代测序技术的短读段特性，对于重复序列的比对可能出现较大偏差，使得后续检测结果出错，所以解决重复区域对于插入变异检测结果的影响是另一大挑战。

解决上述技术问题的意义：

(1)构建新的模型解决重复序列及短读段对于插入变异检测的影响，为DNA插入变异检测提供新的思路和手段，对于变异检测方法研究有重要理论意义。

(2)对于大片段多类型插入变异的准确检测，可以为探究插入变异与癌症发生相关性提供依据，对于从生物角度上研究癌症的发生、发展有重要意义。

(3)高效准确的插入变异检测算法可为相关癌症诊断及靶向药物选定提供高质量的数据基础，为研究癌症机理提供数据分析平台。

发明内容

针对现有技术存在的问题，本发明提供了一种基于新一代测序数据的插入变异检测方法。

本发明是这样实现的，一种基于新一代测序数据的插入变异检测方法，所述基于新一代测序数据的插入变异检测方法包括以下步骤：

第一步，在确定变异发生位点时，有插入变异发生的区域一定会产生分裂读段，针对新序列插入、序列串联倍增、序列散在倍增等插入变异类型及缺失变异、倒置变异的分裂读段分布不同的特性；

第二步，在确定插入变异发生种类及位点之后，通过利用部分匹配、完全匹配、以及未匹配的读段信息构造一条虚拟参考序列，与原始参考序列比较得到插入序列的相关信息；

第三步，利用拷贝数状态信息获得变异基因型。

进一步，所述第一步具体包括：

(1)对于新序列插入，变异位点所覆盖的分裂读段为单一对应混合读段，该变异点下既有前半部分匹配的读段，又有后半部分匹配的读段，且读段在参考序列上均为单一对应；

(2)对于序列串联倍增，只有其倍增序列连接部分会产生分裂读段，分裂读段会在一定区域内进行双映射，同一读段在参考序列上映射到一定区域内的两个位点，且在映射的位点上存在的分裂读段匹配情况一致；

(3)对于序列散在倍增，在其倍增序列左右位点处会产生分裂读段，分裂读段会进行多重映射，同一读段在参考序列上映射到两个以上位点，且其中一个位点是原始序列的开始或终止处，对于映射到原始序列开始或结束处的分裂读段出现原始序列部分匹配。

进一步，所述第一步通过分裂读段特性对变异信号进行筛选分类后，定义每个变异点必须至少覆盖N_SR条reads确定为一个插入变异发生点，N_SR的定义由下式；

其中，R_e表示每个碱基的测序错误率，L_R表示read的长度，N_R表示总的reads数，Q_i表示第i条read的比对质量，

表示下取整；R_e设置为0.01。

进一步，所述第二步检测较长的插入变异时，从每个变异点开始拓展，每次选择左半部分匹配得分最高的分裂读段进行动态拼接后重新比对得到新的分裂读段，当没有分裂读段拼接时，则停止迭代。

进一步，所述第三步利用每个位点的覆盖度信息获得拷贝数状态信息，所有的插入区域其拷贝数状态都为正常，则样本的read覆盖度服从高斯分布N(μ,σ²)，其中μ为样本读段深度均值，样本均值和方差可由固定大小滑窗获得每一滑窗内位点覆盖度均值后计算所得；将零分布转换成一个标准的正态分布后，通过求得插入区域p值，通过对p值的分析可得变异基因型相关信息；

r_i′＝(r_i-μ)/σ；

综上所述，本发明的优点及积极效果为：本发明建立了基于NGS数据特性的变异信息分类方案，为准确靶定插入变异位点及类型提供基础。建立了基于动态模型的插入序列迭代对比检测模型，检测变异类型及插入序列。

本发明能够解决插入变异位点判定不准确的问题；能够解决SR方法检测插入变异造成遗漏的问题；能够解决现有技术遇到重复序列可能会检测出错的问题。

本发明通过建立拷贝数状态数学模型确定变异的基因型相关信息；本发明能有效提取插入变异的特征情况，对插入变异类型进行区分，达到对多类型的插入变异均有较好的检测效果；本发明使用迭代拼接方法，在新一代测序技术所得读段长度偏小的限制下，检测大片段的插入变异及避免了重复区域对插入变异检测的干扰；本发明通过拷贝数状态数学模型，对所得变异检测结果进行基因分型，进一步确定了插入变异发生情况。

对21号染色体中加入101个插入变异后，生成六种不同覆盖度下的仿真数据，本发明与其他方法的检测结果如下表：

由表中数据可以发现，本发明对于仿真数据插入变异的检测在各个覆盖度下准确率与敏感度方面均明显优于其他三种方法。

附图说明

图1是本发明实施例提供的基于新一代测序数据的插入变异检测方法流程图。

图2是本发明实施例提供的基于新一代测序数据的插入变异检测方法实现流程图。

图3是本发明实施例提供的新序列插入示意图。

图4是本发明实施例提供的序列串联倍增示意图。

图5是本发明实施例提供的序列散在倍增示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

针对现有技术大不符合癌症样本中插入变异类型多样的状况，极大限制了癌症诊断及靶向药物选定的能力；对大片度插入变异的检测能力不足；得到错误的变异检测结果的问题。本发明使用双端读端split read与insert size信息准确靶定插入变异发生位点及类型。本发明在技术上使用插入序列迭代拼接方法对插入变异片段进行检测提取。

下面结合附图对本发明的应用原理作详细的描述。

如图1所示，本发明实施例提供的基于新一代测序数据的插入变异检测方法包括以下步骤：

S101：在确定变异发生位点时，有插入变异发生的区域一定会产生分裂读段，针对新序列插入、序列串联倍增、序列散在倍增等插入变异类型及缺失变异、倒置变异等其他变异类型的分裂读段分布不同的特性，设计了变异信息筛选分类方案；

S102：在确定插入变异发生种类及位点之后，通过利用部分匹配、完全匹配、以及未匹配的读段信息来构造一条虚拟参考序列，与原始参考序列比较得到插入序列的相关信息；

S103：利用拷贝数状态信息获得变异基因型。

下面结合附图对本发明的应用原理作进一步的描述。

如图2所示，本发明实施例提供的基于新一代测序数据的插入变异检测方法具体包括以下步骤：

(1)插入变异位点信息分类方案；

确定插入变异发生的位点是检测插入变异非常重要的一个步骤，在确定变异发生位点时，分裂读段(splitread)可以提供非常重要的信息，有插入变异发生的区域一定会产生分裂读段，针对新序列插入、序列串联倍增、序列散在倍增等插入变异类型及缺失变异、倒置变异等其他变异类型的分裂读段分布不同的特性，设计了变异信息筛选分类方案。

如图3所示，对于新序列插入，其变异位点所覆盖的分裂读段为单一对应混合读段(R2、R3、R6、R8)，即该变异点下既有前半部分匹配的读段(R2、R3)，又有后半部分匹配的读段(R6、R8)，且这些读段在参考序列上均为单一对应。

如图4所示，对于序列串联倍增，只有其倍增序列连接部分会产生分裂读段，即图中R2、R3读段为分裂读段，其分裂读段会在一定区域内进行双映射，即同一读段可以在参考序列上映射到一定区域内的两个位点，且在映射的位点上存在的分裂读段匹配情况一致，如图中所示，参考序列左侧位点R2、R3均表现为后半部分匹配，在图中参考序列右侧位点R2、R3均表现为前半部分匹配。

如图5所示，对于序列散在倍增，在其倍增序列左右位点处会产生分裂读段(R2、R3、R6、R7)，其分裂读段会进行多重映射，即同一读段可以在参考序列上映射到两个以上位点，且其中一个位点是原始序列的开始或终止处，对于映射到原始序列开始或结束处的分裂读段出现原始序列部分匹配，如图5所示，对于原始序列开始处，R2、R3均表现为后半部分匹配，对于原始序列结束处，R6、R7均表现为前半部分匹配。其他映射位点为倍增序列散在插入处，在其位点上，两种匹配情况的分裂读段均会出现，即图中插入位点R2、R3表现为前半部分匹配，R6、R7表现为后半部分匹配。

通过分裂读段特性对变异信号进行筛选分类后，为了保证插入变异点检测的可靠性，在考虑测序错误和比对质量的影响下，定义每个变异点必须至少覆盖N_SR条reads才能确定其为一个插入变异发生点，N_SR的定义由式(1)所示；

其中，R_e表示每个碱基的测序错误率，L_R表示read的长度，N_R表示总的reads数，Q_i表示第i条read的比对质量，

表示下取整。R_e这里可以设置为0.01，由于测序平台不同，所以测序错误率也不尽相同，这个值一定程度可以反映测序技术的发展，最后加1是为了保证有相对足够的分裂读段来支持变异发生点。

(2)虚拟序列构建；

在确定插入变异发生种类及位点之后，通过利用部分匹配、完全匹配、以及未匹配的读段信息来构造一条虚拟参考序列，与原始参考序列比较得到插入序列的相关信息。检测较长的插入变异时，从每个变异点开始拓展，每次选择左半部分匹配得分最高的分裂读段进行动态拼接后重新比对得到新的分裂读段，当没有分裂读段可以用来拼接时，则停止迭代。为了提高变异检测的准确性，在开始阶段将参考序列进行处理，使得迭代构造虚拟序列后，可以直接从虚拟序列中提取出插入变异位点及插入序列。

(3)利用拷贝数状态信息获得变异基因型；

利用每个位点的覆盖度信息可以获得拷贝数状态信息，即假定所有的插入区域其拷贝数状态都为正常，则样本的read覆盖度服从高斯分布N(μ,σ²)，其中μ为样本读段深度均值，样本均值和方差可由固定大小滑窗获得每一滑窗内位点覆盖度均值后计算所得。利用式(2)，将零分布转换成一个标准的正态分布后，通过式(3)可求得插入区域p值，通过对p值的分析可得变异基因型相关信息。

r_i′＝(r_i-μ)/σ (2)

将本发明应用于NA18507生物数据，该数据有用生物方法得到的检测结果，相对于基于新一代测序的算法检测，生物方法得到的结果较为准确，以其为标准，将本发明检测结果与其他算法检测方法结果进行对比，如下表所示：

方法	检测数	确认数	确认比
				MindTheGap	12398	160	1.29％
IMSindel	8430	2554	30.30％
				本发明	74908	12499	16.69％

上表对各个方法的插入变异检测数、以生物方法检测变异为标准的确认数及确认比进行了统计，从中可以发现本发明的确认数最为理想，确认比相对较优，综合来看，本发明的检测性能最好。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于新一代测序数据的插入变异检测方法，其特征在于，所述基于新一代测序数据的插入变异检测方法包括以下步骤：

第一步，在确定变异发生位点时，有插入变异发生的区域一定会产生分裂读段，针对新序列插入、序列串联倍增、序列散在倍增插入变异类型及缺失变异、倒置变异的分裂读段分布不同的特性；

第二步，在确定插入变异发生种类及位点之后，通过利用部分匹配、完全匹配、以及未匹配的读段信息来构造一条虚拟参考序列，与原始参考序列比较得到插入序列的相关信息；

第三步，利用拷贝数状态信息获得变异基因型；

所述第一步通过分裂读段特性对变异信号进行筛选分类后，定义每个变异点必须至少覆盖N_SR条reads确定为一个插入变异发生点，N_SR的定义由下式；

其中，R_e表示每个碱基的测序错误率，L_R表示read的长度，N_R表示总的reads数，Q_i表示第i条read的比对质量，

表示下取整；R_e设置为0.01。

2.如权利要求1所述的基于新一代测序数据的插入变异检测方法，其特征在于，所述第一步具体包括：

(1)对于新序列插入，变异位点所覆盖的分裂读段为单一对应混合读段，该变异点下既有前半部分匹配的读段，又有后半部分匹配的读段，且读段在参考序列上均为单一对应；

(2)对于序列串联倍增，只有其倍增序列连接部分会产生分裂读段，分裂读段会在一定区域内进行双映射，同一读段在参考序列上映射到一定区域内的两个位点，且在映射的位点上存在的分裂读段匹配情况一致；

(3)对于序列散在倍增，在其倍增序列左右位点处会产生分裂读段，分裂读段会进行多重映射，同一读段在参考序列上映射到两个以上位点，且其中一个位点是原始序列的开始或终止处，对于映射到原始序列开始或结束处的分裂读段出现原始序列部分匹配。

3.如权利要求1所述的基于新一代测序数据的插入变异检测方法，其特征在于，所述第二步检测较长的插入变异时，从每个变异点开始拓展，每次选择左半部分匹配得分最高的分裂读段进行动态拼接后重新比对得到新的分裂读段，当没有分裂读段拼接时，则停止迭代。

4.如权利要求1所述的基于新一代测序数据的插入变异检测方法，其特征在于，所述第三步利用每个位点的覆盖度信息获得拷贝数状态信息，所有的插入区域其拷贝数状态都为正常，则样本的read覆盖度服从高斯分布N(μ,σ²)，其中μ为样本读段深度均值，样本均值和方差可由固定大小滑窗获得每一滑窗内位点覆盖度均值后计算所得；将零分布转换成一个标准的正态分布后，通过求得插入区域p值，通过对p值的分析可得变异基因型信息。

5.一种利用权利要求1-4任意一项所述基于新一代测序数据的插入变异检测方法的检测系统。

Images