DNAscope 白皮书: 基于机器学习的高精度胚系变异检测流程

2023-08-25
广东
本文字数：4604 字
阅读完需：约 15 分钟

DNAscope 模块，是 Sentieon 软件的一个精准高效的胚系变异检测模块。其在 GATK 基础上优化了核心算法，在继承 GATK 成熟且完整的 BAM 预处理流程的同时，引入机器学习基因分型模型。相比于 GATK 金标准而言，在大幅降低计算成本的情况下，DNAscope 流程能够大幅度提升 SNP 和 Indel 的检测准确度和稳定性。

GATK HaplotypeCaller 由于其准确性高已成为胚系短变异检测的行业标准，在公开发表的第三方基准测试中都取得了最佳性能。然而，包括 HaplotypeCaller 在内的基于短读长序列技术的变异检测流程所检测出来的变异位点与高置信度真集并不能完美匹配，尤其是在基因组的复杂区域，如单碱基多聚体和其他重复位点。这些复杂区域很多都是与临床相关的，随着 NGS 在临床分析中的使用越发频繁，提高这些区域变异识别的准确度将变得越来越重要。

GATK 在复杂区域的限制

GATK 在 reads 比对排列复杂区域无法进行局部重头组装；
GATK 在基因组的特定区域无法进行组装，产生了特定的盲点；
GATK 可能会丢失单倍型之前的联系，无法生成最准确的单倍型序列。

此前，大多数变异检测工具使用的是构建统计模型进行变异位点基因分型，并通过设置手动调节过滤参数来消除变异识别的假阳性。虽然这些过滤器在大多数情况下运行良好，但基于机器学习的假阳性位点过滤及基因分型则可通过学习变异特征之间更复杂的关系来提高准确性，因此相较于传统的显式统计模型，机器学习模型在改进变异过滤及基因分型中的应用越来越广泛，并逐渐成为主流的解决方案。

DNAscope 将基于单倍型拼接的变异检测与机器学习模型结合，从而实现更高的准确性。DNAscope 在 GATK HaplotypeCaller 类似的逻辑体系结构基础上，对活性区域检测和局部组装等模块进行了优化，尤其是在高复杂度区域，可有效提高灵敏度和鲁棒性。当应用机器学习模型时，DNAscope 会输出带有额外注释信息的候选变异列表，然后将这些被注释的候选变异传递到机器学习模型中进行基因分型，从而提高了变异检测和基因分型的准确性。

DNAscope 方法概述

DNAscope 遵循与 GATK HaplotypeCaller 类似的算法流。首先将可能存在遗传变异的位点确定为活跃区域，然后使用 de Bruijn 图对跨活动区域比对的序列 reads 进行局部组装，并通过 PairHMM 计算 read-haplotype 的可能性，从而输出带有额外注释信息的候选变异列表，并传递到基于机器学习的基因分型模型中进行处理，以确定最终正确的变异基因型。

DNAscope 评测

为了评估 DNAscope 的变异检测准确度，我们将 3 个测序深度约为 36x 的 GIAB 样本（HG002，HG003 和 HG004）分别进行降采样得到 5 个覆盖深度 15x-36x 的数据集作为测试集，并采用 Sentieon 的 DNAscope 和 DNAseq（匹配 GATK 结果）标准流程进行测试。DNAscope 的性能评测参考的是美国国家标准与技术研究所 (NIST) 于 2021 年下半年发布更新的 GIAB high-confidence calls v4.2.1，该版本的真集涵盖了对二代测序来说较为困难的区域，所以准确率等指标不能和旧版真集的评测数据进行直接对比。此外，这里我们还省略了会导致 F1-scores 下降的变异位点质量值重新校正(VQSR)的步骤。

评测结果显示，对于 SNP 和 INDEL，DNAscope 在所有测试集中的灵敏度（recalls）和特异性（precisions）都显著高于 DNAseq。例如在 30x 的 HG002 数据集中，DNAscope 的 SNP 和 INDEL 的 F1-score 分别高达 99.57% 和 99.46%。此外，相较于 DNAseq，DNAscope 的错误变异检测结果降低了 1 倍以上（87607 到 34634）；同时在多个数据集的测试结果中表明，DNAscope 没有出现模型过拟合的情况。另外值得一提的是，模型文件和测序平台可以特异性适配，这在已发表的华大测序仪模型的性能评测文章中有印证，说明了 DNAscope 可以适配新的测序技术。

DNAscope VS DNAseq 灵敏度曲线

图 a-b 分别为 30x HG002、HG003 和 HG004 样本的 SNPs 及 Indels 评测结果。DNAscope 结果准确度总体高于 DNAseq。

对全基因组进行多层次的分区分析能够更好地体现分析工具的准确性和稳健度，特别是在低可比对区域等的复杂区域，这些区域的变异检测准确性评估能够更加全面地展示分析工具的检测性能。因此，我们采用了全球基因组学与健康联盟（Global Alliance for Genomics and Health，GA4GH）的标准对基因复杂区域进行多层次分区分析，以评估 DNAscope 和 DNAseq 的准确度。下图显示 DNAscope 的综合分析性能较 DNAseq 更为出色，尤其在复杂区域的分析精确度较 DNAseq 有较大幅度的提升。

30x HG002 的 GA4GH 分区评估

图 a-b 分别为 30x HG002、HG003 和 HG004 样本的 SNPs 及 Indels 评测结果。DNAscope 结果准确度总体高于 DNAseq。

我们同样使用 HG002、HG003 和 HG004 三个样本分别进行降采样得到的 5 个不同测序深度的数据集来评估两个分析流程变异检测精确度随测序深度的变化波动情况。测试结果符合预期，两个流程在低深度的精确度相较于高深度都有一定程度的下降。而 DNAscope 由于对算法架构进行了整体优化，同时引入了机器学习模型，以至于 DNAscope 下降的幅度相对较低，值得一提的是，DNAscope 在 20x 的准确度表现已经优于与 DNAseq 流程在 36x 的表现。

HG002、HG003、HG004 多测序深度分析

图 a-b 分别显示 DNAscope 和 DNAseq 分析结果 SNP 及 Indel 的 F1 评分。

多倍体重测序在农业基因组研究中十分普遍，但 DNAscope 的机器学习模型是基于二倍体样本建立的，所以当前版本的 DNAscope 并不支持非二倍体样本的分析。对此，我们尝试用 DNAscope 内置的贝叶斯基因分型模型（Bayesian genotyping model）对 15x 测序深度样本进行精确度评估，来模拟较难进行高深度测序的大型基因组序列分析。正如预期，DNAscope 的贝叶斯基因分型模型整体的精确度比机器学习模型低，但在 INDEL 的表现上仍优于 DNAseq，在 SNP 方面则与 DNAseq 表现相当。

15X HG002、HG003 和 HG004 样本分析

DNAscope_Bayesian 以及 DNAseq 对 15x 测序深度的 HG002、HG003 和 HG004 样本进行分析的精确度对比，图中显示 SNP 及 Indel 结果的 F1 评分。

最后我们也在标准云计算环境下对 DNAscope 的运行速度进行了基准测试。我们在 4 个 AWS C6i 实例上分别对 30x 的 HG002 测序数据进行分析。如下图所示，在计算规模为 96 线程或以上时，DNAscope 完成一个 30x 的测序数据分析用时低于 1 小时，分析速度与 DNAseq 相当且比开源的 BWA/GATA 分析流程快了将近 5 倍。此外，测试结果显示，DNAscope 分析速度几乎与使用线程数呈线性关系，由此说明 DNAscope 的可扩展性优异，可对运算资源进行有效利用。

DNAscope 并行计算扩展性评测

DNAScope 在多个 AWSC6i 实例上的运行耗时。图中显示每个实例类型不同运行阶段的时长。

DNAscope 评测结论

总体而言，DNAscope 在不同的样本和测序深度下，精准度高于 DNAseq。GA4GH 基因区域多层次分区分析的结果表明，DNAscope 在大多数复杂区域都能维持卓越的准确度表现。此外，DNAscope 在融合了 GATK HaplotypeCaller 完整数学模型的基础上，引入机器学习技术，在实现高准确性的同时，保持优异的计算效率。

在 30x 的测序深度下，DNAscope 在 SNP 与 Indel 变异检测中都达到了 99.5%左右的准确度(F1-score)，相比于 GATK 总体错误个数降低了两倍多，展现了显著的准确率提升。
在基因组分区分析中，发现 DNAscope 在比对(mapping)困难的高复杂区域，例如 MHC 区域中，相比 GATK 有着更加显著的准确度优势。
当测序深度降低时，DNAscope 也可以保持相对较高的准确率，例如 20x 的 DNAscope 分析准确率已经超越了 36x 的 GATK 结果。
DNAscope 速度比 GATK 快 5-10 倍，在 96 线程节点上可以 1 小时完成 30x 全基因组数据分析。
支持非人物种分析，在使用内置贝叶斯模型分析 15x 深度数据时，准确率高于 GATK 流程。

DNAscope 方法学

候选变异检测中的方法改进

Sentieon 的 Haplotyper （DNAseq）复现了 GATK HaplotypeCaller 数学模型，在做到结果高度匹配的同时还提高了算法的运行效率。Haplotyper 和 DNAscope 都使用与 HaplotypeCaller 相同的算法架构，包括活性区域检测、局部组装、基因型计算；同时通过去除降采样，实现软件的方法改进，资源管理优化等技术提高鲁棒性。此外， DNAscope 还包含其他的算法改进及下文所介绍的变异位点注释。

变异位点注释

DNAscope 新增了 GATK 流程中所没有的额外的变异位点注释，这些额外添加的注释与原有标准注释将共同在基因分型过程中被输入到机器学习模型中。新增注释如下：

Entropy：定义为在局部组装过程中，所有被识别的单倍型的香农熵（Shannon entropy）。高熵（high entropy）是可能存在序列比对位置错误或比对结果错误的指标。

UDP：所有输入样本未经过滤的测序深度。该指标类似于深度“DP”，但是“DP”统计中会删除比对质量低或比对错误可能性高的序列。

NBQ：变异位点本身及其两边各 5bp 区域的的所有位点的平均碱基质量。

ReadPosEndDist：该变异位点到所有覆盖序列两端的平均距离。

基于模型的基因分型

DNAscope 基因分型模型的设计能够将真实的胚系变异从系统噪声中分辨出来。通过机器学习，识别文库构建，测序和序列比对过程中出现的错误范式（error patterns），从检测困难区域提取有用信息。为了在分析准确性与计算效率方面超过深度学习方法，我们选择梯度提升机（Gradient Boosting Machines，GBMs）作为框架来学习结构化基因数据的错误范式。GBMs 通过连续建立决策树，以序列集成的方式训练弱基学习器（base learner），逐步减少残差。在此过程中，模型通过学习率来控制每个基学习器的逐步改进过程，一旦在验证过程中观察到过拟合就可以尽早中止训练过程。

我们使用来自不同数据源的 GIAB HG001 和 HG005（不含 20 号染色体）的标准品样本训练模型，其中 90%的变异位点用作训练集，剩余的 10%作为验证集，以优化模型的复杂度，在性能与泛用性方面寻求平衡，同时防止过拟合。模型过滤步骤的输入是高灵敏度模式下运行的 DNAscope 生成的 VCF 文件，此模式下只需要少量的序列支持就能够检测出候选的变异位点。随后，假定样本在所有位点上都是二倍体条件下，模型利用 VCF 输出文件的位点注释作为特征将候选位点分类为 4 种基因型：“0/0”、“0/1”、“1/1”、“1/2”。基因型“0/0”表示报告的候选变异位点实际和参考基因组一致，并不是真实的变异位点。从最终变异位点列表 VCF 中可以看出，DNAscope 基因分型模型大大降低了错误发现率（false discovery rate），在赋予样本基因型的同时保持较高的灵敏度。此外，模型还在最终 VCF 文件中添加了“ML_PROB”注释，即模型计算出的此位点为假阳性的可能性。

测试基准

测试所使用的 FASTQ 数据均来自于公开数据，在“数据”部分有详细说明。

所有的 FASTQ 数据均使用标准化数据处理流程进行处理，包括使用经 Sentieon 优化的 BWA-MEM v202112.01 进行序列比对；使用 v202112.01 版本的 Sentieon 模块标记重复；使用 SAMtools v1.3.1（适用于某些样本）进行降采样；使用 DNAscope v202112.01、DNAscope Illumina 模型 v1.0 和 Haplotyper 202112.01 进行变异检测；使用 hap.py v0.3.10（ RTGtools vcfeval v3.8.2 作为引擎），比对 NIST 真集 v4.2.1 用于进行变异位点准确性评估。

文件和工具

分析期间生成的中间文件和输出文件可以在 gs://dnascope-benchmarking 中找到。

使用 Google Colaboratory 上的交互式和协作式 Jupyter 笔记本分析 Hap.py 输出和其他输出文件。

创作场景