随着DNA测序和下游数据处理成本的不断降低，群体测序研究开始应用于之前未曾处理过的分析规模。队列层面变异目录是祖源研究、罕见变异分析和基因型/表型关联发现的关键资源。因此，检出集必须高度准确，但数据处理和分析挑战仍然存在。
目前有多种工具可以汇总来自大量队列的检出集，其中一个示例是GATK及其Joint Genotyping工作流程1，该品牌多年来一直是行业标杆。但是，它包含几个只有GATK检出的变异集才需要的步骤，不适用于其他变异检出程序，例如变异评分重新校准（VQSR）。此外，还有DNAnexus开发的GLnexus2，此软件与其云端平台紧密结合。至于如何在其他云端平台或HPC上进行大规模部署，目前尚不清楚。许多现有工具面临的另一个挑战是如何进行高效的迭代分析，以便将新样本及其变异检出添加到现有的检出集中。
DRAGEN gVCF Genotyper是因美纳用于在群体规模上对小生殖系变异进行汇总和基因分型的解决方案。gVCF Genotyper的输入是一组由DRAGEN生殖系变异检出程序编写的gVCF文件3。输出结果是一个多样本VCF文件（msVCF），其中包含队列中发现的所有变异的基因型。输出结果还将包含队列指标，可用于质量控制和变异筛选。
DRAGEN gVCF Genotyper采用迭代工作流程，可将新样本添加到现有队列中。该工作流程可让用户高效地将新批次样本与现有批次样本合并，而无需重复处理。
用户可以将输出的msVCF输入下游工具，如Hail (https://hail.is) 和PLINK (https://www.cog-genomics.org/plink), 进行数据探索或进行关联研究。
本文概述了gVCF Genotyper的工作流程，展示如何进行大规模部署，并使用千人基因组计划（1KGP）队列演示了DRAGEN gVCF Genotyper的性能。
我们在多个内部项目中使用了gVCF Genotyper，例如使用DRAGEN 4.0重新分析了1KGP队列 (https://www.internationalgenome.org) 和多达500,000个样本的更大队列。外部用户包括Genomics England4,5和其他组织。
这款软件是因美纳新一代测序数据分析解决方案DRAGEN的一部分。在DRAGEN 4.0中，使用分布式处理模式时，gVCF Genotyper可扩展至数十万样本。执行gVCF Genotyper的推荐（也是最具可扩展性的）平台是Illumina Connected Analytics（ICA），但我们也支持在DRAGEN服务器和HPC部署上进行本地分析。
大规模基因分型
DRAGEN gVCF Genotyper依赖于gVCF输入格式，其中既包含类似VCF的变异信息，也包含特定位置（纯合参考位置）不存在变异的置信度。
基因分型工作流程由几个步骤组成：首先，从每个输入文件中提取变异等位基因，存储变异基因型和其他指标。下一步，将收集某个位置上所有等位基因的变异位点。我们会使用过滤器，剔除低质量的等位基因和未被检出基因型支持的等位基因，包括使用该位置纯合参照样本的信息。
随后，对变异等位基因进行归一化处理。我们从输入中的等位基因标签到位点中的等位基因顺序计算等位基因图谱，下面给出了一个示例（表1）。我们还更新了每个样本的基因型和QC指标，以匹配归一化的等位基因集。
可扩展的VCF表示形式
gVCF Genotyper将输出一个msVCF，其中包含为队列中所有样本计算的指标。
随着队列规模的增加，msVCF可能会变得很大——在某些情况下需要的内存会超过VCF解析器可以分配的内存量，从而对内存和运行时间造成影响。这是由VCF条目引起的，这些条目存储了等位基因的每个可能组合的值，例如FORMAT/PL，因此随着等位基因数量的增加而呈现平方增长。在基因组复杂度较低的区域，每个位点的等位基因数量可能很大，它还会随着样本数量的增加而急剧增加。
由于上述原因，我们用FORMAT/LPL标签代替了FORMAT/PL，FORMAT/LPL标签只存储样本中出现的等位基因的值，而不是队列中在该位点发现的所有等位基因的值。我们还为每个样本添加了一个新的FORMAT/LAA字段，该字段将列出相对于全局等位基因列表在该样本中出现的替代等位基因的单碱基标签。
这种方法被称为本地等位基因，也被BCFtools（https://samtools.github.io/bcftools/bcftools.html）和Hail（https://hail.is/docs/0.2/experimental/vcf_combiner.html#local-alleles）等开源软件所采用。
msVCF还包含每个等位基因的等位基因计数和近交系数。等位基因计数用于计算等位基因频率，除其他外，等位基因频率是决定一个变异是良性还是可能致病的有用指标。
近交系数（IC）衡量的是变异位点上的过高杂合性。此系数还可用于检测和过滤read定位杂峰：6基因组中杂合性过高的位置可能是read定位不正确的区域，这些区域看起来是杂合的，因为一些read实际上属于基因组的不同部分。
近交系数的计算公式为：1-（观察到的杂合基因型数量）/（预期杂合基因型数量）。如果这个数量是负数，则意味着我们观察到的杂合子数量比根据哈迪-温伯格假设所预期的数量要多（https://en.wikipedia.org/wiki/Hardy%E2%80%93Weinberg_principle）。
一个良好的经验法则是过滤近交系数为负数的变异位点6。图1展示了为什么这是一个有用的方法：负的近交系数可以捕捉到没有被其他过滤器（如检出率（NS_GT））过滤掉的高杂合性位点。
下图是DRAGEN gVCF Genotyper编写的msVCF文件，该文件使用了之前解释的“本地等位基因”形式（图2）。
GATK等开源软件编写的是全套基因型似然值（FORMAT/PL），而DRAGEN gVCF Genotyper和Hail等其他工具编写的则是仅针对样本中实际存在的等位基因（“本地等位基因”）的精简集。
与GATK最佳实践的比较
GATK（基因组分析工具包）最佳实践建议采取几个步骤来分析来自大型队列的变异检出，包括重新校准质量分值、输入数据库和完善基因型（也称为联合基因分型）。
联合基因分型是指一类利用队列信息提高基因分型准确性的算法，GATK团队是这一方法的应用先驱。
提高基因分型的准确性很重要，但我们已经证明7，DRAGEN变异检出并不需要GATK式的联合基因分型算法，因为它并不能提高准确性。原因是GATK算法会尝试去除变异杂峰，但这些杂峰已经在DRAGEN的上游进行了过滤。
由于GATK联合基因分型算法也是一项计算昂贵的操作，我们建议用户只运行DRAGEN gVCF Genotyper，而不对DRAGEN变异检出进行GATK式联合基因分型。
VQSR也是如此，对于使用新版DRAGEN的变异检出集，不再推荐使用VQSR（https://community.illumina.com/s/question/0D53l00006lvBDBCA2/why-vqsr-was-removed-from-latest-versions-of-dragen）。
单服务器分析与分布式处理
如果有多个批次的样本，gVCF Genotyper工作流程可采用分步模式运行；如果只有一个批次的样本，则可采用单服务器分析模式运行。
用于gVCF Genotyper的FASTA参考必须与用于构建DRAGEN哈希表（用于比对与变异检出）的参考相同。
单服务器分析模式采用一组gVCF输入文件，并将其汇总为单一msVCF，这对于针对有限数量的样本快速生成包含指标和基因型的msVCF非常有用。在此模式下，无法进行迭代分析。
DRAGEN 4.0用户指南中提供了单服务器分析工作流程的详细信息8。
如图3所示的分步模式包括几个步骤，用于在多个节点之间分配计算。工作流程可扩展到数十万个样本。
步骤1（gVCF汇总）：将一批gVCF文件汇总为队列文件和Census文件。队列文件是一种压缩数据格式，用于存储来自多个样本的gVCF数据，类似于多样本gVCF。Census文件存储了队列中每个样本的所有变异和纯合参考基因型区块的汇总统计数据。如果在处理的样本较多，用户应将样本分成多个批次，每个批次的样本量相近（例如1000个样本），然后分别对每个批次运行步骤1。
步骤2（Census汇总）：生成所有批次Census文件后，用户可将它们汇总为一个全局Census文件。这一步骤可扩展为数千个批次。
步骤3（生成msVCF）：生成全局Census文件后，用户可将其与每批队列和Census文件结合使用，为每批样本生成msVCF文件。
为了便于在多个计算节点上并行处理，对于上述每个步骤，用户都可以选择将基因组分割成大小相等的片段（使用用户指南中描述的命令行选项“-shard”），并在不同的计算节点上使用一个迭代gVCF Genotyper的实例处理每个片段。基因组片段会沿着染色体边界断开，从而对片段数量施加了一个下限。
此工作流程会生成（样本批次数）*（基因组片段数）msVCF输出文件。如果需要合并所有批次的msVCF，可以运行一个附加步骤，将所有批次的msVCF文件合并为一个包含所有样本的msVCF。
有新一批样本时，用户只需对该批样本执行步骤1，然后将新一批样本的Census文件与之前所有批次的全局Census文件合并，并写入新的全局Census文件。每次更新全局Census文件，发现新的变异位点和/或更新现有变异位点的变异统计信息时，都可以生成包含新信息的新分批msVCF文件。
下图4显示了将新一批样本添加到现有队列中所需的步骤。对于新批次中的gVCF文件，只需执行步骤1即可建立新的每批次队列和Census文件。之后，工作流程将执行步骤2以建立新的全局Census文件，该文件需要更新步骤3中每个批次的msVCF文件。
云端和本地可用性
DRAGEN gVCF Genotyper是DRAGEN 4.0的一部分，可以像所有其他组件一样在DRAGEN服务器上运行。但是，这只适用于小规模队列：在单服务器分析模式下，处理1000个gVCF文件大约需要24小时。
运行迭代gVCF Genotyper的工作流程在ICA中实施，是我们向拥有大规模队列（＞10,000个样本）的用户推荐的选项。ICA工作流程已在100,000多个生殖系样本中进行了测试。
gVCF Genotyper可在ICA平台上使用。对于ICA v2，可使用命令行界面（CLI）11。
对于从gVCF到msVCF的单服务器分析工作流程，在ICA中运行gVCF Genotyper的成本约为0.3 iCredits/样本。每增加一个新样本，迭代分析的成本约为0.3 iCredits，每更新一个已处理样本的msVCF，成本约为0.06 iCredits。这一估算基于在500个内核上运行的DRAGEN 4.0和在AWS区域us-east-1中运行的ICA。
在大多数情况下，该平台允许用户在成本与运行时间之间进行权衡——例如，使用较少的ICA节点可以降低成本，但会增加获得结果的时间。
我们还提供软件模式的gVCF Genotyper二进制版本，可在商用硬件和HPC上运行。
gVCF Genotyper的运行时间取决于磁盘I/O的速度。输入文件的大小也是一个重要因素。
因此，在规划gVCF基因分型项目时，请考虑使用DRAGEN生殖系变异检出程序中的紧凑型gVCF输出模式（命令行选项--vc-compact-gvcf）。这使得gVCF文件的体积缩小了40%-50%，并大大提高了文件解析软件（如gVCF Genotyper）的运行速度，存储成本也相应降低。这种紧凑型gVCF模式省略了gVCF中的一些指标，这些指标仅用于从头检测谱系中的变异，这意味着紧凑型gVCF不能与DRAGEN谱系检出程序联合使用。
使用案例：ICA千人基因组计划分析
七年前，1KGP发布了一套由低覆盖度测序组成的开放式数据集（来自26个群体的2504个无关个体的数据）。这是初次为编制人类基因变异目录而进行的大规模测序工作。
纽约基因组中心的研究人员与其他科研伙伴合作，将这一资源扩展到许多亲子家系。扩展集包含3202个样本，其中有602个家系样本，可公开获取9。
纽约基因组中心使用GATK检出该队列中的变异。因美纳使用DRAGEN 4.0重新分析了样本，并公开了这些样本。此数据发布将在单独的博客文章中介绍。
我们使用DRAGEN gVCF Genotyper对所有3202个样本进行了汇总和基因分型。分析成本为990 iCredits，相当于每个样本0.31 iCredits。此分析获得一个msVCF，其中包含队列中所有变异的基因型和频率。分析中还收集了有关缺失基因型和偏离哈代-温伯格假设的统计数据，用近交系数表示（图1）。
我们通过两个示例展示基因分型msVCF的多功能性和准确性：
隐性家系的孟德尔不一致基因型
我们计算了属于队列10（样本NA20891、NA20882和NA20900）的隐性家系的孟德尔错误率。计算违背家族关系的基因型检出是评估准确性的有用指标，因为它们并不局限于基因组的高置信度区域。
下表3显示了在家系中至少一个成员具有变异的情况下，孟德尔错误数量占位点总数量的比例。表中比较了GATK检出集与DRAGEN 4.0和3.5.7生成的检出集。
1KGP中假定的细胞系杂峰分析
纽约基因组中心对整个队列中只出现一次的变异（单例，等位基因数[AC]=1）进行了分析9。他们的假设是，这些变异中有很多是细胞系的杂峰。他们还比较了队列中有亲属样本中单例的数量与无亲属样本中单例的数量。子女中单例数量最少，其次是父母，然后是队列中没有任何亲属的样本。
DRAGEN gVCF Genotyper可即时计算许多变异指标，其中包括等位基因计数。因此，无需借助其他工具即可轻松提取和分析单例变异。上图5显示了每个样本的单例变异分类。该图显示了样本如何根据其亲缘关系状态进行分组；子女往往拥有最少的单例，其次是父母和无亲缘关系的样本。
参考文献
1. Derek Caetano-Anolles (2022) “Calling variants on cohorts of samples using the HaplotypeCaller in GVCF mode” https://gatk.broadinstitute.org/hc/en-us/articles/360035890411-Calling-variants-on-cohorts-of-samples-using-the-HaplotypeCaller-in-GVCF-mode
2. Michsel F Lin (2018) “GLnexus: joint variant calling for large cohort sequencing “ https://www.biorxiv.org/content/10.1101/343970v1
3. Derek Caetano-Anolles (2022) gVCF – Genomic Variant Call Format https://gatk.broadinstitute.org/hc/en-us/articles/360035531812-GVCF-Genomic-Variant-Call-Format
4. Kousathanas et al (2022) “Whole-genome sequencing reveals host factors underlying critical COVID-19” https://www.nature.com/articles/s41586-022-04576-6 Nature volume 607, pages 97–103
5. Genomics England press release on covid 19 patient study https://www.genomicsengland.co.uk/news/study-insights-severe-covid-19
6. Derek Caetano-Anolles (2022) “Inbreeding coefficient” https://gatk.broadinstitute.org/hc/en-us/articles/360035531992-Inbreeding-Coefficient
7. Population genetics data processing with the DRAGENTM Bio-IT Platform https://jp.support.illumina.com/content/dam/illumina/gcs/assembled-assets/marketing-literature/dragen-popgen-tech-note-m-gl-00561/dragen-popgen-tech-note-m-gl-00561.pdf
8. Instructions for using the Illumina DRAGEN Bio-IT Platform https://support-docs.illumina.com/SW/DRAGEN_v40/Content/SW/DRAGEN/gVCFGenotyper.htm
9. Marta Byrska-Bishop et al (2022) “High-coverage whole-genome sequencing of the expanded 1000 Genomes Project cohort including 602 trios” Cell, Volume 185, https://www.cell.com/cell/fulltext/S0092-8674(22)00991-6#%20
10. Roslin N, Li W, Paterson AD, Strug LJ. Quality control analysis of the 1000 Genomes Project Omni2.5 genotypes. bioRxiv. 2016;078600–078600.
11.PopGen CLI在ICA中以捆绑包形式提供，其名称取决于ICA订阅区域。目前我们在以下3个地区有授权捆绑包：
a.美国东部（use1）授权包名称为“popgen-cli-release-use1”
b.伦敦（euw2）授权包名称为“popgen-cli-release-euw2”
c.新加坡（apse1）授权包名称为“popgen-cli-release-apse1”