[发明专利]一种基因组变异检测方法及检测装置

说明书

本发明公开了一种基因组变异检测方法及检测装置，涉及生物信息学研究领域，以解决现有基因组变异检测效率低下、且检测精度不高的问题。该方法包括：确定基因组的潜在变异区域，对落在潜在变异区域内的测序序列进行局部组装创建迭代DeBruijn图，根据迭代DeBruijn图得到单倍体；将测序序列与所有单倍体进行比对，将测序序列对应的最佳单倍体与潜在变异区域对应的参考序列进行比对后，得到一组比对结果，采用预设变换规则对该组比对结果中的两个比对结果进行变换，得到测序序列与潜在变异区域对应的参考序列的校正对齐结果，根据测序序列与潜在变异区域对应的参考序列的校正对齐结果，确定基因组的变异类型。

技术领域

本发明实施例涉及生物信息研究领域，尤其涉及一种基因组变异检测方法及检测装置。

背景技术

随着基因组测序成本的持续下降，高通量测序仪(如：Illumina Hiseq系列测序仪)产生的基因组测序序列(read)呈现爆照式增长，尤其是人类基因组测序序列的积累速度更为明显，如何从大量基因组测序序列得到高质量的基因组变异结果，成为一项富有挑战性的工作。

目前，人们通常利用短序列比对软件(如：Burrows Wheeler Aligner，BWA)将每条测序序列比对到参考序列(reference sequence)上，得到每条测序序列与参考序列的双序列比对结果(包括测序序列相对于参考序列详细的匹配、错配、插入和删除等信息)，然后根据所有测序序列与参考序列的双序列比对结果，得到基因组变异结果。但由于BWA等短序列比对软件是把每个read单独比对到参考序列上，并不考虑测序序列之间整体上是否对齐，很容易把原本属于同一种变异类型的测序序列错误地比对成包含不同变异类型的测序序列，造成误判。

为解决该问题，现有技术人员通常采用基因分析工具(英文：The GenomeAnalysis Toolkit，GATK)中的子模块HaplotypeCaller(简称：HC)来进行基因组变异检测，以消除测序序列间比对不齐造成的影响。其中，HC检测主要包括下述4个步骤：1)确定基因组中的潜在变异区域(ActiveRegion)；2)对每个潜在变异区域执行局部组装(localassemby)，根据组装后的序列以及kmer长度(Size)得到DeBruijn图(DeBruijn graph)，并对DeBruijn graph进行处理得到至少一个单倍体(haplotype)；3)利用PairHMM计算每个read与每个haplotype的最佳比对、以及每个read与每个haplotype对应的条件概率P(read|haplotype)；4)对每个“ActiveRegion”中的每个潜在变异位置，基于二倍体假设，计算出贝叶斯后验概率Pr(D|G)最大的一对haplotypes，将该对haplotypes作为该潜在变异位置的变异结果，其中，G表示基因型(Genotype)，D表示覆盖潜在变异区域中的每个变异位置的测序序列。

虽然HC检测能有效避免测序序列间不对齐造成的影响，但仍存在下述几种问题：1)由于在上述local assembly阶段使用了多个不同的kmer size，且对每个不同的kmersize分别建立一张DeBruijn graph并独立地遍历生成haplotypes，所以会导致haplotypes数目偏多，进而导致haplotypes跟参考序列进行Smith-Waterman(简称SW)比对的次数偏多，影响了运行时间。2)局部组装精度不够高；3)由于是基于二倍体假设，所以只支持二倍体变异，不支持多倍体变异，无法很好地支持多倍体变异检测；4)即便在二倍体情况下，计算出的Pr(D|G)也不够精准，例如：在二倍体假设下，覆盖每个变异位置的测序序列，应该是一部分来自于H1，剩下一部分来自于H2，且每个测序序列间相互独立，所以Pr(D|G)的最大值应该是来自于覆盖每个变异位置的测序序列的最佳分割：一部分来自H1，剩下一部分来自H2，且连乘的积最大，然而，GATK只是通过下述计算公式，通过计算均值的连乘来近似得到Pr(D|G)：

所以，由上可知现有HC基因组变异检测方法运行时间较长，效率低下，且检测精度不高。