CN111583997B - 杂合变异下校正第三代测序数据中测序错误的混合方法 - Google Patents

杂合变异下校正第三代测序数据中测序错误的混合方法 Download PDF

Info

Publication number
CN111583997B
CN111583997B CN202010373513.2A CN202010373513A CN111583997B CN 111583997 B CN111583997 B CN 111583997B CN 202010373513 A CN202010373513 A CN 202010373513A CN 111583997 B CN111583997 B CN 111583997B
Authority
CN
China
Prior art keywords
generation sequencing
long
corrected
heterozygosity
site
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
Application number
CN202010373513.2A
Other languages
English (en)
Other versions
CN111583997A (zh
Inventor
王嘉寅
刘佳琦
赖欣
萧笑
张选平
朱晓燕
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN202010373513.2A priority Critical patent/CN111583997B/zh
Publication of CN111583997A publication Critical patent/CN111583997A/zh
Application granted granted Critical
Publication of CN111583997B publication Critical patent/CN111583997B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • 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
    • G16B20/30—Detection of binding sites or motifs
    • 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
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Medical Informatics (AREA)
  • Biophysics (AREA)
  • Theoretical Computer Science (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Chemical & Material Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Molecular Biology (AREA)
  • Genetics & Genomics (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)

Abstract

本发明公开了一种杂合变异下校正第三代测序数据中测序错误的混合方法，输入数据为第二代测序数据和第三代测序数据，利用已有的比对软件和组装软件对输入数据进行处理，基于贝叶斯分类器原理对基因位点的杂合性进行判断，结合杂合判断的结果对第三代测序数据中的读段进行校正，解决了现有校正算法在处理杂合变异时的低准确度和无效性的问题。本发明在校正测序错误时考虑了杂合变异，设计了一系列概率模型对杂合性进行判断和分类，再针对不同的杂合性分类采用不同的校正策略，解决了已有校正方法遇到杂合变异时出现校正错误的问题。

Description

杂合变异下校正第三代测序数据中测序错误的混合方法
技术领域
本发明属于第三代测序技术领域，具体涉及一种杂合变异下校正第三代测序 数据中测序错误的混合方法。
背景技术
基因组测序技术，尤其是单分子长读测序技术，也称为第三代测序(英文名 称：Third Generation Sequence，英文缩写：TGS)，彻底改变了基因组学的研究。 TGS技术不仅延续了下一代测序(英文名称：Next Generation Sequence,英文缩写： NGS)技术的高通量等优点，而且产生的读段长度更长，可达到10kbp。因此， TGS技术为许多领域带来了巨大推动力，例如结构变异的检测，甲基化的鉴定以 及疾病的诊断。虽然TGS在读段长度和其他许多方面处于领先地位，但由于测 序技术的限制，其测序错误率远高于NGS，这也给进一步的研究带来困难。目前 来看，虽然TGS的测序错误率随着技术的进步而逐渐降低，但在过去十年中积 累的PB级测序数据也不能丢弃，因此，对TGS数据的测序错误校正仍是本领域 研究的热点与难点。
目前，国内外的科研人员在第三代测序数据校正算法的研究上已取得了一些 成果，尽管这些算法在整体校正准确性上的表现尚可，但是对杂合变异的校正效 果较差。现有的方法一般分为两类：自校正(英文名称：self-correction)和混合 校正(英文名称：hybrid-correction)。自校正的核心思想是通过在第三代测序长 读段(英文名称：LongRead，英文缩写：LR)数据之间建立多重比对并计算局 部比对，进而找到LR之间的一致序列。这种方法实际上是基于多重比对和局部 比对来估计杂合变异，但其对杂合变异的校正性能通常被LR的覆盖度所限制。 当前，由于测序成本的原因，已有的TGS数据集的覆盖度往往很低，极大限制 了自校正的应用，也使得对杂合变异的准确校正更难实现。例如，当LR的覆盖 度小于2时，从数学期望的角度考虑，几乎不可能将杂合变异与测序错误或纯合 变异区分开。
由于自校正存在上述问题，混合校正在实践中更受欢迎。混合校正的基本思 想是：给定LR和第二代测序短读段(英文名称：Short Read，英文缩写：SR)， 将SR比对到LR上，然后针对该比对结果进行投票，获得最多投票的等位基因 就是最终的校正结果。可以看出，混合校正的核心是投票，一些最新研究也是基 于投票过程进行改进，而当前的混合校正算法不能解决杂合变异校正的原因也在 于算法本身的结构问题。图1展示了错误校正的示例，混合校正算法错误处理杂 合变异示意图。L和S分别表示第三代测序和下一代测序数据集，Li表示L中的第 i个长读段，Si表示S中的第i个短读段，
Figure BDA0002479252190000021
表示Li的第j个碱基，
Figure BDA0002479252190000022
表示Si的第j个 碱基。可以看出，由于仅通过碱基个数比进行投票校正，本不应该被校正的杂合变异
Figure BDA0002479252190000023
被错误地从等位基因A校正为等位基因C。此外，即使在SR上考虑杂合 变异，由于每个长读段被独立处理，也不能实现杂合变异和噪声的区分。另一方 面，算法的目的是校正LR，因此为了控制成本，SR的覆盖率低，这也不利于杂 合变异的识别。除此之外，每校正一个长读段都需要将所有SR重新比对一遍， 进而导致了校正的低效率。上述这些都使得现有方法校正杂合变异时的表现很 差。
然而，在许多情况下，杂合变异比纯合变异更为普遍，并且杂合性在疾病的 基因型-表型分析和遗传学研究中起着重要作用。区分杂合变异与测序错误是正 确处理杂合变异的关键和难点，现有方法中简单的投票过程无法处理此类复杂状 况。
发明内容
本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种杂合 变异下校正第三代测序数据中测序错误的混合方法，解决由于现有校正算法的单 一投票机制等算法结构问题所带来的杂合变异误校正，突破测序数据的读段长度 以及覆盖度对校正算法准确度尤其是对杂合变异的校正准确度的负面影响，实现 低覆盖度测序数据杂合变异校正的精确度提高。
本发明采用以下技术方案：
杂合变异下校正第三代测序数据中测序错误的混合方法，包括以下步骤：
S1、输入第二代测序短读段数据S和第三代测序长读段数据L，得到重叠群 并按顺序首尾相接得到一个基因序列；
S2、把第二代测序短读段数据S和第三代测序长读段数据L与伪参考序列 Ref比对，分别获得第二代测序短读段数据S和第三代测序长读段数据L的比对 成功长读段数据集Lm和比对失败长读段数据集Lu，并储存；
S3、分别为第二代测序短读段数据S和第三代测序长读段数据L建立两组概 率模型，计算纯合性和杂合性的后验概率，并以后验概率的值高的一方作为判断 结果；
S4、根据步骤S3得出的位点Ri杂合判断结果，校正位点Ri下对应的比对成 功长读段数据集Lm中的所有位点，遍历伪参考序列Ref上的所有位点，重复以上 校正策略，直到比对成功长读段数据集Lm中的所有读段都被校正；
S5、使用步骤S3得出的第二代测序短读段数据S的杂合性判断结果和对比 对失败长读段数据集Lu进行校正，将校正后的结果混合，得到L的完整校正结 果。
具体的，步骤S2中，比对成功长读段数据集Lm为比对一致比例大于等于 90％的读段，比对失败长读段数据集Lu为剩余的读段。
具体的，具体步骤如下：
S301、统计位点Ri下第二代测序短读段数据S和第三代测序长读段数据L的 比对结果中的碱基具体类型和数量和出位点Ri下的长读段深度RDi和短读段深度 rdi，其中，第三代测序长读段数据L的比对结果中的碱基类型按出现频次从大到 小排序，用Xq表示，对应出现频次用|Xq|表示，第二代测序短读段数据S的比对 结果中的碱基类型和数量分别用xq和|xq|表示，q＝1,2,3,4,5；
S302、明确第三代测序长读段数据L比对结果中的碱基类型Xq和第二代测序 短读段数据S比对结果中碱基类型xq中各碱基的条件概率，根据基因组学碱基分 布模型，碱基分布模型服从二项分布；
S303、根据先验概率和条件概率基于贝叶斯概率模型计算当前比对结果下位 点Ri下第二代测序短读段数据S和第三代测序长读段数据L比对结果的杂合与纯 合的概率；
S304、根据步骤S303计算出的位点Ri下第二代测序短读段数据S和第三代 测序长读段数据L比对结果的杂合与纯合的概率，取概率值大的作为杂合性判断 结果。
进一步的，步骤S302中，对于Xq中的各碱基，假设服从概率为P1的二项分 布，即Xq～Bin(Dq,P1)，则Xq发生|Xq|次的概率为：
Figure BDA0002479252190000041
Figure BDA0002479252190000051
其中，P1是第三代测序技术的测序误差先验概率，Dq为位点Ri下来自长读段 的碱基的比对深度；
对于xq中的各碱基，假设服从概率为P2的二项分布，即xq～Bin(dq,P2)，则xq发 生|xq|次的概率为：
Figure BDA0002479252190000052
Figure BDA0002479252190000053
其中，P2是第二代测序技术的测序误差先验概率，dq为位点Ri下来自短读段 的碱基的比对深度。
进一步的，步骤S303中，对应位点Ri下的比对结果碱基种类为：
kind＝1时，第二代测序短读段数据S和第三代测序长读段数据L均判定为纯 合；
kind＝2时，对于第三代测序长读段数据L，即|X1|+|X2|＝RDi，对于第二代测 序短读段数据S，即|x1|+|x2|＝rdi，根据贝叶斯概率模型分别计算纯合与杂合的后 验概率；
kind＝3时，对于第三代测序长读段数据L，即|X1|+|X2|+|X3|＝RDi，根据贝叶 斯概率模型分别计算纯合与杂合的后验概率P(c is homozygosity|{X1,X2,X3})和 P(c isheterozygosity|{X1,X2,X3})；对于第二代测序短读段数据S，即|x1|+|x2|+|x3|＝rdi， 根据贝叶斯概率模型分别计算纯合与杂合的后验概率P(c is homozygosity|{x1,x2,x3}) 和P(c is heterozygosity|{x1,x2,x3})；
kind＝4时，对于第三代测序长读段数据L，即|X1|+|X2|+|X3|+|X4|＝RDi，根据 贝叶斯概率模型分别计算纯合与杂合的后验概率P(c is homozygosity|{X1,X2,X3,X4}) 和P(c is heterozygosity|{X1,X2,X3,X4})；对于第二代测序短读段数据S，即 |x1|+|x2|+|x3|+|x4|＝rdi，根据贝叶斯概率模型分别计算纯合与杂合的后验概率 P(c is homozygosity|{x1,x2,x3,x4})和P(c is heterozygosity|{x1,x2,x3,x4})；
kind＝5时，对于第三代测序长读段数据L，即|X1|+|X2|+|X3|+|X4|+|X5|＝RDi， 根据贝叶斯概率模型分别计算纯合与杂合的后验概率 P(c is homozygosity|{X1,X2,X3,X4,X5})和P(c is heterozygosity|{X1,X2,X3,X4,X5})；对于第 二代测序短读段数据S，即|x1|+|x2|+|x3|+|x4|+|x5|＝rdi，根据贝叶斯概率模型分别计 算纯合与杂合的后验概率P(c ishomozygosity|{x1,x2,x3,x4,x5})和 P(c is heterozygosity|{x1,x2,x3,x4,x5})。
具体的，步骤S4中，当第二代测序短读段数据S和第三代测序长读段数据L 均判断为纯合，则位点Ri的判断结果为纯合，对应待校正位点按照第二代测序短 读段数据S的比对结果校正，即，将待校正位点校正为第二代测序短读段数据S 在位点Ri下的比对结果中出现频次最高的碱基；
当第二代测序短读段数据S和第三代测序长读段数据L均判断为杂合，则位 点Ri的判断结果为杂合，对应待校正位点按照第二代测序短读段数据S的比对结 果校正，即，将待校正位点校正为第二代测序短读段数据S在位点Ri下的比对结 果中出现频次前二高的碱基之一
当第二代测序短读段数据S为杂合，第三代测序长读段数据L判断为纯合， 则位点Ri的判断结果为纯合，对应待校正位点按照第三代测序长读段数据L的比 对结果校正，即，将待校正位点校正为第三代测序长读段数据L在位点Ri下的比 对结果中出现频次最高的碱基；
当第二代测序短读段数据S为纯合，第三代测序长读段数据L判断为杂合， 则位点Ri的判断结果为纯合，对应待校正位点按照第二代测序短读段数据S的比 对结果校正，即，将待校正位点校正为第二代测序短读段数据S在位点Ri下的比 对结果中出现频次最高的碱基。
具体的，步骤S5具体为：
S501、任意提取一个对比对失败长读段数据集Lu中的待校正长读段read1， 用比对工具把第二代测序短读段数据S比对到read1上；
S502、提取read1的第一个位点r1以及该位点下第二代测序短读段数据S的 比对结果，根据步骤S3的杂合性判断步骤对该位点进行杂合性判断；
S503、根据步骤S502的杂合判断结果对位点r1进行校正；
S504、按顺序遍历read1上的所有位点，重复步骤S502和S503，直到read1 上的所有位点都被校正；
S505、遍历提取对比对失败长读段数据集Lu中的所有待校正读段，重复步 骤S501～S504，直到对比对失败长读段数据集Lu中的所有读段都被校正。
进一步的，步骤S503中，若判断位点r1为杂合，则将位点r1校正为S在该位 点下的比对结果中出现频次前二高的碱基之一；若判断位点r1为纯合，则将位点r1校正为S在该位点下的比对结果中出现频次最高的碱基。
与现有技术相比，本发明至少具有以下有益效果：
本发明一种考虑了杂合变异的基于综合概率模型的混合校正方法，解决了由 于现有校正算法的单一投票机制等算法结构问题所带来的杂合变异误校正，客服 了测序数据的读段长度以及覆盖度对校正算法准确度尤其是对杂合变异的校正 准确度的负面影响，实现低覆盖度测序数据杂合变异的精确校正。
进一步的，对第三代测序长读段数据集L进行组装，并按顺序首尾连接输出 的重叠群contigs，得到一个伪参考序列Ref，可以选择该领域中现有的任何基因 组组装工具完成；构建Ref的过程可以初步获得L中各个读段的分歧，是后续计 算的统计性指标。
进一步的，将第二代测序短读段数据集S和L同时比对到Ref上，针对L的 比对结果，根据比对一致比例划分为比对成功长读段数据集Lm和比对失败长读 段数据集Lu。用Ref配合读段比对可以避免重复比对，显著改善了现有方法由于 多次重复比对导致计算效率低的问题，提高本发明的计算速度；同时，定位比对 到同一位点的长读段和短读段，用于后续分析。
进一步的，在同一Ref位点i下，分别针对Lm和S比对结果中的碱基占比 建立贝叶斯概率模型，并根据贝叶斯概率的大小进行该位点的杂合性判断。相较 于现有校正方法，该在校正操作前加入概率模型进行位点的杂合性判断，基本解 决了现有方法在杂合变异校正上的低准确度问题，并且还不损失其他正常位点的 校正准确率。
进一步的，针对Ref位点i的杂合判断结果，对比对到该位点下的Lm中所 有长读段的对应位点进行校正，校正策略依据杂合判断结果的不同而有区别，这 也反映了本发明与现有校正算法的不同之处在于对待杂合变异和纯合位点采取 了不同的校正策略。顺序遍历Ref的每一个位点并重复操作直到Lm中的所有长 读段完成校正。
进一步的，对于Lu，提取集合中的长读段，将S比对到其上，建立贝叶斯概 率模型进行杂合性判断，进一步根据结果进行校正。遍历Lu中每一个长读段并 重复操作直到Lu中的所有长读段完成校正。
综上所述，本发明在校正测序错误时考虑了杂合变异，设计了一系列概率模 型对杂合性进行判断和分类，再针对不同的杂合性分类采用不同的校正策略，解 决了已有校正方法遇到杂合变异时出现校正错误的问题。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
图1为混合校正算法错误处理杂合变异示意图；
图2为本发明流程图；
图3为Lk、St、
Figure BDA0002479252190000091
Ri、Ref的关系与分布展示图。
具体实施方式
本发明提供了一种杂合变异下校正第三代测序数据中测序错误的混合方法 QIHC(QI Heterozygosity Correction)，输入数据为第二代测序数据(以下简称为 S)和第三代测序数据(以下简称为L)，利用已有的比对软件和组装软件对输入 数据进行处理，基于贝叶斯分类器原理对基因位点的杂合性进行判断，结合杂合 判断的结果对L中的读段进行校正，解决了现有校正算法在处理杂合变异时的低 准确度和无效性的问题。
本发明基于以下在学术界具有普遍共识的假设：
1.根据现有测序技术标准，一个位点的测序结果可能有5种，包括碱基：A、 T、C、G和1个未知状态N。
2.根据基因组学的碱基分布模型，碱基分布的概率模型服从二项分布。
请参阅图2，本发明一种杂合变异下校正第三代测序数据中测序错误的混合 方法，包括以下步骤：
S1、制备伪参考序列
本发明的输入数据为第二代测序短读段数据(以下简称为S)和第三代测序 长读段数据(以下简称为L)，无需提前准备参考序列。首先制备一个伪参考序列 (以下简称为Ref)。具体地，使用任意一种常用的基因组组装软件，例如Canu， 将L输入组装软件进行组装，得到重叠群(英文名称：contigs)并按顺序首尾相 接，得到一个长度约为10000bp(其中bp是基因组长度单位)的基因序列。
S2、获得读段比对文件
使用任意一种常用的基因组比对软件，例如Blasr，把输入数据S和L比对 到伪参考序列Ref上，分别获得S和L的比对结果SAM文件(英文名称： SequenceAlignment/Map，英文缩写：SAM)并储存。对于L的比对结果，根据 比对一致比例划分为比对成功长读段数据集(以下简称为Lm)和比对失败长读 段数据集(以下简称为Lu)。具体地，比对一致比例大于等于90％的读段划分到 Lm，否则划分到Lu。
S3、建立概率模型判断杂合性
不同于已有方法，本发明在校正步骤前加入杂合性判断，以便于根据被判断 位点的杂合性进行校正。基于贝叶斯分类器原理，分别为S和L建立两组概率模 型，计算纯合性和杂合性的后验概率，并以概率值较高的一方作为判断结果。
首先，为了方便描述和理解，同时也为了规范数据格式，做出以下定义和说 明：
Lk：L中的第k个长读段；
St：S中的第t个短读段；
Figure BDA0002479252190000101
Lk上的第m个碱基；
Figure BDA0002479252190000102
St上的第n个碱基；
Ri：Ref上的第i个位点；
请参阅图3，直观展示了Lk、St、
Figure BDA0002479252190000114
Ri和Ref的关系与分布。
为了便于理解，我们先针对位点Ri进行描述，其余位点同理。具体包含以下 步骤：
S301、统计位点Ri的对应比对结果。这一步要统计出位点Ri下的L和S的比 对结果中的碱基具体类型和数量。其中，L的比对结果中的碱基类型按出现频次 从大到小排序，用Xq表示，对应出现频次用|Xq|表示，q＝1,2,3,4,5；S的比对结果 中的碱基类型和数量分别用xq和|xq|表示，q＝1,2,3,4,5。另外，还要统计出位点Ri下的长读段深度和短读段深度，分别用RDi和rdi表示。
S302、明确碱基分布的概率模型。这一步要明确Xq和xq中各碱基的条件概率， 根据基因组学碱基分布模型，碱基分布模型服从二项分布。
具体地，对于Xq中的各碱基，假设服从概率为P1的二项分布，即 Xq～Bin(Dq,P1)，则Xq发生|Xq|次的概率为
Figure BDA0002479252190000111
Figure BDA0002479252190000112
其中，P1是第三代测序技术的测序误差先验概率。
具体地，对于xq中的各碱基，假设服从概率为P2的二项分布，即xq～Bin(dq,P2)， 则xq发生|xq|次的概率为
Figure BDA0002479252190000113
其中，P2是第 二代测序技术的测序误差先验概率。
S303、计算L和S的后验概率。这一步骤根据之前的先验概率和条件概率基 于贝叶斯概率模型计算当前比对结果下位点Ri分别是杂合与纯合的后验概率。原 理公式如下：
Figure BDA0002479252190000121
其中，c表示纯合(英文名称：homozygosity)或杂合(英文名称：heterozygosity)。需要注意的是，并不是每个位点下的比对结果中出现的碱基种 类都能多达全部类型，对于没有出现的碱基类型，直接在公式(1)中去掉对应部分 即可。
具体地，在实际计算中，将情况分为5种，分别对应位点Ri下的比对结果碱 基种类为1、2、3、4和5的情况，为了描述简便，用kind表示比对结果中的碱 基种类：
1、kind＝1时，无论L还是S，均判定为纯合。
2、kind＝2时，对于L，即|X1|+|X2|＝RDi，根据公式(1)分别计算纯合与杂合的 后验概率，即
Figure RE-GDA0002531370130000111
Figure RE-GDA0002531370130000112
其中，
Figure BDA0002479252190000124
p1为第三代测序技术的测序误差先验概率，可由具体 的测序技术得到，p3为纯合的先验概率，p4为杂合的先验概率，p5＝0.5，表示当 前碱基来自父母的概率各是0.5。
对于S，即|x1|+|x2|＝rdi，根据公式(1)分别计算纯合与杂合的后验概率，即
Figure BDA0002479252190000131
Figure BDA0002479252190000132
其中，
Figure BDA0002479252190000133
p2为第二代测序技术的测序误差先验概率，可由具体 的测序技术得到。
3、kind＝3时，对于L，即|X1|+|X2|+|X3|＝RDi，根据公式(1)分别计算纯合与杂 合的后验概率P(c is homozygosity|{X1,X2,X3})和P(c is heterozygosity|{X1,X2,X3})。
对于S，即|x1|+|x2|+|x3|＝rdi，根据公式(1)分别计算纯合与杂合的后验概率 P(c is homozygosity|{x1,x2,x3})和P(c is heterozygosity|{x1,x2,x3})。
4、kind＝4时，对于L，即|X1|+|X2|+|X3|+|X4|＝RDi，根据公式(1)分别计算纯合 与杂合的后验概率P(c is homozygosity|{X1,X2,X3,X4})和 P(c is heterozygosity|{X1,X2,X3,X4})。
对于S，即|x1|+|x2|+|x3|+|x4|＝rdi，根据公式(1)分别计算纯合与杂合的后验概率 P(c is homozygosity|{x1,x2,x3,x4})和P(c is heterozygosity|{x1,x2,x3,x4})。
5、kind＝5时，对于L，即|X1|+|X2|+|X3|+|X4|+|X5|＝RDi，根据公式(1)分别计 算纯合与杂合的后验概率P(c is homozygosity|{X1,X2,X3,X4,X5})和 P(c isheterozygosity|{X1,X2,X3,X4,X5})。
对于S，即|x1|+|x2|+|x3|+|x4|+|x5|＝rdi，根据公式(1)分别计算纯合与杂合的后验 概率P(c is homozygosity|{x1,x2,x3,x4,x5})和P(c is heterozygosity|{x1,x2,x3,x4,x5})。
S304、根据S303步骤计算出的位点Ri下的L的比对结果的杂合与纯合的概 率，取概率值较大者作为杂合性判断结果，同理，S的结果也做同样处理。
S4、校正Lm
这一步骤校正Lm中的长读段，根据步骤S3得出的位点Ri杂合判断结果，校 正位点Ri下对应的Lm中的所有位点。
具体地，在校正过程中，由L和S给出的杂合性判断结果有四种可能组合， 本发明针对这四种组合分别给出不同校正策略：
L和S均判断为纯合，则位点Ri的判断结果为纯合，对应待校正位点按照S 的比对结果校正，即，将待校正位点校正为S在位点Ri下的比对结果中出现频次 最高的碱基；
L和S均判断为杂合，则位点Ri的判断结果为杂合，对应待校正位点按照S 的比对结果校正，即，将待校正位点校正为S在位点Ri下的比对结果中出现频次 前二高的碱基之一；
L判断为纯合，S判断为杂合，则位点Ri的判断结果为纯合，对应待校正位 点按照L的比对结果校正，即，将待校正位点校正为L在位点Ri下的比对结果中 出现频次最高的碱基；
L判断为杂合，S判断为纯合，则位点Ri的判断结果为纯合，对应待校正位 点按照S的比对结果校正，即，将待校正位点校正为S在位点Ri下的比对结果中 出现频次最高的碱基。
遍历Ref上的所有位点，重复以上校正策略，直到Lm中的所有读段都被校 正。
S5、校正Lu
这一步骤校正Lu中的长读段。Lu中的长读段由于未比对成功，因此认为它 们之间没有一致性等关系，本发明只用S的杂合性判断结果和比对结果对Lu进 行校正。
S501、任意提取一个Lu中的待校正长读段read1，用比对工具把S比对到read1 上；
S502、提取read1的第一个位点r1以及该位点下S的比对结果，根据步骤S3 的杂合性判断步骤对该位点进行杂合性判断；
S503、根据S502的杂合判断结果对位点r1进行校正；
具体地，若判断位点r1为杂合，则将位点r1校正为S在该位点下的比对结果 中出现频次前二高的碱基之一；若判断位点r1为纯合，则将位点r1校正为S在该 位点下的比对结果中出现频次最高的碱基；
S504、按顺序遍历read1上的所有位点，重复S502和S503，直到read1上 的所有位点都被校正；
S505、遍历提取Lu中的所有待校正读段，重复S501-S504，直到Lu中的所 有读段都被校正。
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实 施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所 描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中 的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因 此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的 本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本 领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属 于本发明保护的范围。
为了验证本发明的有效性，我们对本发明与已有校正算法Canu和Jabba在 杂合变异处的校正效果进行比较。在模拟实验中，为了有效集中地展示校正算法 对杂合变异的校正效果，在此仅展示杂合变异的校正结果。对于每个杂合变异， 我们关注校正后的杂合性变化。具体地，判断一个位点是否仍为杂合变异的标准 如下：将校正后的长读段比对到参考序列(hg19)上，观察杂合变异下相应的比 对碱基的分布，如果该分布满足杂合性分布，则该位置保持杂合性；否则，该位 点失去其杂合性。评价指标如下：真阳性(TP)位点是在校正后仍保持杂合性的位 点、假阴性(FN)位点是在校正后没有保持杂合性的位点、精确度(accuracy＝1-error rate)。
为了进行对比实验，本发明将第三代测序数据L的覆盖度设置为3×，5×， 10×，12×和15×，读段长度设置为10000bps；第二代测序数据S的覆盖度设置为 5×，10×，15×，20×和50×，读段长度设置为100bps。实验中使用的第三代测序 数据集L包含500个杂合变异位点。仿真数据生成工具选用PBSIM，比对工具 选用Blasr，参数设置为-header，-m5。针对L的不同覆盖度，分别用QIHC、Canu 和Jabba对其进行校正。结果如表1所示：
表1：QIHC、Canu和Jabba的校正精确度结果比较
Figure BDA0002479252190000161
Figure BDA0002479252190000171
可见，三种方法的结果相比较，Jabba的精确度明显低于QIHC和Canu，这 也说明早期的校正方法在校正过程中完全没有考虑杂合变异。对于Canu和QIHC 的结果，当L的覆盖度为3×时，QIHC的精确度比Canu高15个百分点。随着覆 盖度的升高，QIHC的表现也一直比Canu出色，直到覆盖度达到12×后，Canu 的精确度才略有反超，但差距并不大。这也说明，本发明更适合第三代测序数据 的低覆盖度特征。
在校正后的杂合性质量方面，本发明也进行了一系列的实验来验证其有效 性。由于上一个实验已经验证Jabba在杂合变异的校正效果远不及QIHC和Canu， 因此这里只比较QIHC和Canu的校正结果。将L的覆盖度设置为15×。结果如 表2所示：
表2：QIHC和Canu的杂合性质量结果比较
Figure BDA0002479252190000172
杂合性质量是指校正后的杂合变异在保持了杂合性的同时，其校正碱基是否 落在正确的范围之内。例如，由等位基因A和C组成的杂合变异，在校正后， 比对到该位点的碱基仍应是A和C占比更大；否则，尽管该位点在校正后仍保 持杂合性，但其杂合性质量非常低。为了更清楚地分析杂合性质量，我们对其进 行量化。具体来说，对于一个A-C杂合变异，我们将比对到该位点的碱基A和C 的占比与碱基T和G的占比进行比较。如果前者大于后者，即差值为正，则杂 合性质量高(表示为positive)，否则杂合性质量低(表示为negative)。更具体地， 对具有高杂合性质量的位点进行细分，将其划分为优秀(表示为excellence)和良好(表示为good)，将差值在0与0.3之间定义为良好，在0.3与1之间定义为 优秀。
由表2结果可见，相较于Canu，QIHC在杂合性质量方面更占优势，这也是 本实验选择15×覆盖度进一步分析杂合性质量的原因，即，当QIHC在15×覆盖 度的校正精确度结果不占优势时，它仍然在杂合性质量方面明显领先于Canu。
综上所述，本发明一种杂合变异下校正第三代测序数据中测序错误的混合方 法。
以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡 是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发 明权利要求书的保护范围之内。

Claims (7)

1.杂合变异下校正第三代测序数据中测序错误的混合方法，其特征在于，包括以下步骤：
S1、输入第二代测序短读段数据S和第三代测序长读段数据L，得到重叠群并按顺序首尾相接得到一个基因序列；
S2、把第二代测序短读段数据S和第三代测序长读段数据L与伪参考序列Ref比对，分别获得第二代测序短读段数据S和第三代测序长读段数据L的比对成功长读段数据集Lm和比对失败长读段数据集Lu，并储存；
S3、分别为第二代测序短读段数据S和第三代测序长读段数据L建立两组概率模型，计算纯合性和杂合性的后验概率，并以后验概率的值高的一方作为判断结果，具体步骤如下：
S301、统计位点Ri下第二代测序短读段数据S和第三代测序长读段数据L的比对结果中的碱基具体类型和数量和出位点Ri下的长读段深度RDi和短读段深度rdi，其中，第三代测序长读段数据L的比对结果中的碱基类型按出现频次从大到小排序，用Xq表示，对应出现频次用|Xq|表示，第二代测序短读段数据S的比对结果中的碱基类型和数量分别用xq和|xq|表示，q＝1,2,3,4,5；
S302、明确第三代测序长读段数据L比对结果中的碱基类型Xq和第二代测序短读段数据S比对结果中碱基类型xq中各碱基的条件概率，根据基因组学碱基分布模型，碱基分布模型服从二项分布；
S303、根据先验概率和条件概率基于贝叶斯概率模型计算当前比对结果下位点Ri下第二代测序短读段数据S和第三代测序长读段数据L比对结果的杂合与纯合的概率；
S304、根据步骤S303计算出的位点Ri下第二代测序短读段数据S和第三代测序长读段数据L比对结果的杂合与纯合的概率，取概率值大的作为杂合性判断结果；
S4、根据步骤S3得出的位点Ri杂合判断结果，校正位点Ri下对应的比对成功长读段数据集Lm中的所有位点，遍历伪参考序列Ref上的所有位点，重复以上校正策略，直到比对成功长读段数据集Lm中的所有读段都被校正；
S5、使用步骤S3得出的第二代测序短读段数据S的杂合性判断结果和对比对失败长读段数据集Lu进行校正，将校正后的结果混合，得到L的完整校正结果。
2.根据权利要求1所述的杂合变异下校正第三代测序数据中测序错误的混合方法，其特征在于，步骤S2中，比对成功长读段数据集Lm为比对一致比例大于等于90％的读段，比对失败长读段数据集Lu为剩余的读段。
3.根据权利要求1所述的杂合变异下校正第三代测序数据中测序错误的混合方法，其特征在于，步骤S302中，对于Xq中的各碱基，假设服从概率为P1的二项分布，即Xq～Bin(Dq,P1)，则Xq发生|Xq|次的概率为：
Figure FDA0003457332510000021
Figure FDA0003457332510000022
其中，P1是第三代测序技术的测序误差先验概率，Dq为位点Ri下来自长读段的碱基的比对深度；
对于xq中的各碱基，假设服从概率为P2的二项分布，即xq～Bin(dq,P2)，则xq发生|xq|次的概率为：
Figure FDA0003457332510000023
Figure FDA0003457332510000024
其中，P2是第二代测序技术的测序误差先验概率，dq为位点Ri下来自短读段的碱基的比对深度。
4.根据权利要求1所述的杂合变异下校正第三代测序数据中测序错误的混合方法，其特征在于，步骤S303中，用kind表示比对结果中的碱基种类，对应位点Ri下的比对结果碱基种类为：
kind＝1时，第二代测序短读段数据S和第三代测序长读段数据L均判定为纯合；
kind＝2时，对于第三代测序长读段数据L，即|X1|+|X2|＝RDi，对于第二代测序短读段数据S，即|x1|+|x2|＝rdi，根据贝叶斯概率模型分别计算纯合与杂合的后验概率；
kind＝3时，对于第三代测序长读段数据L，即|X1|+|X2|+|X3|＝RDi，根据贝叶斯概率模型分别计算纯合与杂合的后验概率P(c is homozygosity|{X1,X2,X3})和P(c isheterozygosity|{X1,X2,X3})；对于第二代测序短读段数据S，即|x1|+|x2|+|x3|＝rdi，根据贝叶斯概率模型分别计算纯合与杂合的后验概率P(c is homozygosity|{x1,x2,x3})和P(cis heterozygosity|{x1,x2,x3})；
kind＝4时，对于第三代测序长读段数据L，即|X1|+|X2|+|X3|+|X4|＝RDi，根据贝叶斯概率模型分别计算纯合与杂合的后验概率P(c is homozygosity|{X1,X2,X3,X4})和P(c isheterozygosity|{X1,X2,X3,X4})；对于第二代测序短读段数据S，即|x1|+|x2|+|x3|+|x4|＝rdi，根据贝叶斯概率模型分别计算纯合与杂合的后验概率P(c is homozygosity|{x1,x2,x3,x4})和P(c is heterozygosity|{x1,x2,x3,x4})；
kind＝5时，对于第三代测序长读段数据L，即|X1|+|X2|+|X3|+|X4|+|X5|＝RDi，根据贝叶斯概率模型分别计算纯合与杂合的后验概率P(c is homozygosity|{X1,X2,X3,X4,X5})和P(c is heterozygosity|{X1,X2,X3,X4,X5})；对于第二代测序短读段数据S，即|x1|+|x2|+|x3|+|x4|+|x5|＝rdi，根据贝叶斯概率模型分别计算纯合与杂合的后验概率P(c ishomozygosity|{x1,x2,x3,x4,x5})和P(c is heterozygosity|{x1,x2,x3,x4,x5})。
5.根据权利要求1所述的杂合变异下校正第三代测序数据中测序错误的混合方法，其特征在于，步骤S4中，当第二代测序短读段数据S和第三代测序长读段数据L均判断为纯合，则位点Ri的判断结果为纯合，对应待校正位点按照第二代测序短读段数据S的比对结果校正，即，将待校正位点校正为第二代测序短读段数据S在位点Ri下的比对结果中出现频次最高的碱基；
当第二代测序短读段数据S和第三代测序长读段数据L均判断为杂合，则位点Ri的判断结果为杂合，对应待校正位点按照第二代测序短读段数据S的比对结果校正，即，将待校正位点校正为第二代测序短读段数据S在位点Ri下的比对结果中出现频次前二高的碱基之一
当第二代测序短读段数据S为杂合，第三代测序长读段数据L判断为纯合，则位点Ri的判断结果为纯合，对应待校正位点按照第三代测序长读段数据L的比对结果校正，即，将待校正位点校正为第三代测序长读段数据L在位点Ri下的比对结果中出现频次最高的碱基；
当第二代测序短读段数据S为纯合，第三代测序长读段数据L判断为杂合，则位点Ri的判断结果为纯合，对应待校正位点按照第二代测序短读段数据S的比对结果校正，即，将待校正位点校正为第二代测序短读段数据S在位点Ri下的比对结果中出现频次最高的碱基。
6.根据权利要求1所述的杂合变异下校正第三代测序数据中测序错误的混合方法，其特征在于，步骤S5具体为：
S501、任意提取一个对比对失败长读段数据集Lu中的待校正长读段read1，用比对工具把第二代测序短读段数据S比对到read1上；
S502、提取read1的第一个位点r1以及该位点下第二代测序短读段数据S的比对结果，根据步骤S3的杂合性判断步骤对该位点进行杂合性判断；
S503、根据步骤S502的杂合判断结果对位点r1进行校正；
S504、按顺序遍历read1上的所有位点，重复步骤S502和S503，直到read1上的所有位点都被校正；
S505、遍历提取对比对失败长读段数据集Lu中的所有待校正读段，重复步骤S501～S504，直到对比对失败长读段数据集Lu中的所有读段都被校正。
7.根据权利要求6所述的杂合变异下校正第三代测序数据中测序错误的混合方法，其特征在于，步骤S503中，若判断位点r1为杂合，则将位点r1校正为S在该位点下的比对结果中出现频次前二高的碱基之一；若判断位点r1为纯合，则将位点r1校正为S在该位点下的比对结果中出现频次最高的碱基。
CN202010373513.2A 2020-05-06 2020-05-06 杂合变异下校正第三代测序数据中测序错误的混合方法 Active CN111583997B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010373513.2A CN111583997B (zh) 2020-05-06 2020-05-06 杂合变异下校正第三代测序数据中测序错误的混合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010373513.2A CN111583997B (zh) 2020-05-06 2020-05-06 杂合变异下校正第三代测序数据中测序错误的混合方法

Publications (2)

Publication Number Publication Date
CN111583997A CN111583997A (zh) 2020-08-25
CN111583997B true CN111583997B (zh) 2022-03-01

Family

ID=72126225

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010373513.2A Active CN111583997B (zh) 2020-05-06 2020-05-06 杂合变异下校正第三代测序数据中测序错误的混合方法

Country Status (1)

Country Link
CN (1) CN111583997B (zh)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013028739A1 (en) * 2011-08-25 2013-02-28 Complete Genomics Phasing of heterozygous loci to determine genomic haplotypes
CN108460245A (zh) * 2017-02-21 2018-08-28 深圳华大基因科技服务有限公司 使用三代序列优化二代组装结果的方法和装置
CN108460246A (zh) * 2018-03-08 2018-08-28 北京希望组生物科技有限公司 一种基于三代测序平台的hla基因分型方法
CN108629156A (zh) * 2017-03-21 2018-10-09 深圳华大基因科技服务有限公司 三代测序数据纠错的方法、装置和计算机可读存储介质
CN110010193A (zh) * 2019-05-06 2019-07-12 西安交通大学 一种基于混合策略的复杂结构变异检测方法
CN110621785A (zh) * 2017-06-20 2019-12-27 深圳华大生命科学研究院 基于三代捕获测序对二倍体基因组单倍体分型的方法和装置

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2013028739A1 (en) * 2011-08-25 2013-02-28 Complete Genomics Phasing of heterozygous loci to determine genomic haplotypes
CN108460245A (zh) * 2017-02-21 2018-08-28 深圳华大基因科技服务有限公司 使用三代序列优化二代组装结果的方法和装置
CN108629156A (zh) * 2017-03-21 2018-10-09 深圳华大基因科技服务有限公司 三代测序数据纠错的方法、装置和计算机可读存储介质
CN110621785A (zh) * 2017-06-20 2019-12-27 深圳华大生命科学研究院 基于三代捕获测序对二倍体基因组单倍体分型的方法和装置
CN108460246A (zh) * 2018-03-08 2018-08-28 北京希望组生物科技有限公司 一种基于三代测序平台的hla基因分型方法
CN110010193A (zh) * 2019-05-06 2019-07-12 西安交通大学 一种基于混合策略的复杂结构变异检测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Draft Sequencing of the Heterozygous Diploid Genome of Satsuma (Citrus unshiu Marc.) Using a Hybrid Assembly Approach;Tokurou Shimizu1 等;《frontiers in Genetics》;20171205;全文 *
Jabba: hybrid error correction for long sequencing reads;Giles Miclotte 等;《Algorithms for Molecular Biology》;20160503;全文 *
复杂基因组测序技术研究进展;高胜寒 等;《遗传》;20181130;第40卷(第11期);全文 *

Also Published As

Publication number Publication date
CN111583997A (zh) 2020-08-25

Similar Documents

Publication Publication Date Title
US8725422B2 (en) Methods for estimating genome-wide copy number variations
EP3298523B1 (en) Methods and systems for copy number variant detection
US10777301B2 (en) Hierarchical genome assembly method using single long insert library
NZ759659A (en) Deep learning-based variant classifier
US20220130488A1 (en) Methods for detecting copy-number variations in next-generation sequencing
US20160117444A1 (en) Methods for determining absolute genome-wide copy number variations of complex tumors
CN110832597A (zh) 基于深度神经网络的变体分类器
WO2010059235A2 (en) Algorithms for sequence determination
US20130138358A1 (en) Algorithms for sequence determination
CN107133493B (zh) 基因组序列的组装方法、结构变异探测方法和相应的系统
CN107480470B (zh) 基于贝叶斯与泊松分布检验的已知变异检出方法和装置
WO2013166517A1 (en) Methods for determining absolute genome-wide copy number variations of complex tumors
CN103114150A (zh) 基于酶切建库测序与贝叶斯统计的单核苷酸多态性位点鉴定的方法
WO2010051320A2 (en) Methods for assembling panels of cancer cell lines for use in testing the efficacy of one or more pharmaceutical compositions
CN111583998B (zh) 一种考虑拷贝数变异因素的基因组结构变异分型方法
US20200350037A1 (en) System, method and computer accessible-medium for multiplexing base calling and/or alignment
CN111583997B (zh) 杂合变异下校正第三代测序数据中测序错误的混合方法
WO2019132010A1 (ja) 塩基配列における塩基種を推定する方法、装置及びプログラム
WO2015004016A1 (en) Transcript determination method
CN109935274B (zh) 一种基于k-mer分布特征的长读数重叠区检测方法
CN104424398A (zh) 碱基序列对准系统及方法
US20240127905A1 (en) Integrating variant calls from multiple sequencing pipelines utilizing a machine learning architecture
US20230207050A1 (en) Machine learning model for recalibrating nucleotide base calls corresponding to target variants
WO2024073519A1 (en) Machine-learning model for refining structural variant calls
AU2022301321A1 (en) Machine-learning model for generating confidence classifications for genomic coordinates

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