这是我们"大规模基因组分析"系列文章的第二篇，在第一篇文章中，我们探讨了一个简单的问题：当对大量基因组进行测序时，如何提供实时聚合。我们通过使用Delta-Lake和使用sparksql构建的流媒体管道解决了这个问题。在这个博客中，我们关注更高级的联合基因分型过程，包括将来自多个个体的变异调用合并到一个群体的单一视图中。这是基因组学中最常见和最复杂的问题之一。在Databricks，我们利用分布式计算、存储和云基础设施的创新，并将其应用于基因组学，以帮助解决阻碍组织执行联合基因分型能力的问题，"N+1"问题，以及扩展到群体水平的挑战。我们针对基因组学的统一分析平台提供了一个优化的管道，只需点击一下就可以扩展到大规模的集群和数千个样本。在这个博客中，我们探讨如何将这些创新应用于联合基因分型。在我们深入研究联合基因分型之前，首先让我们讨论一下人们为什么要进行大规模测序。大多数人都熟悉23andMe或祖先dna产生的遗传数据。这些测试使用基因分型阵列，它读取基因组中固定数量的变异，通常约有1000000个众所周知的变异，这些变异通常发生在正常人群中。通过测序，我们可以得到一个个体所有变异的无偏图片，不管它们是我们以前在健康人身上见过很多次的变异，还是我们以前从未见过的导致或预防疾病的变体。图1展示了这两种方法之间的区别。 图1：这张图说明了基因型阵列（左）和测序（中）以及联合基因分型（右）产生的变异数据之间的差异。基因分型阵列被限制为"读取"固定数量的已知变体，但保证每个变体的每个样本都有一个基因型。在测序中，我们能够发现非常罕见的变异，它们只存在于一个人身上，但是确定一个新的变异是否真的在这个人身上是独一无二的，还是仅仅用当前的技术很难检测到，这是一个非常重要的问题。虽然测序提供了更高的分辨率，但当我们试图检查一个基因变体对许多患者的影响时，我们遇到了一个问题。由于基因分型阵列测量的是所有样本中相同的变体，所以对许多个体进行观察是一个直截了当的命题：所有变异都是在所有个体中测量的。在处理测序数据时，我们有一个更为棘手的命题：如果我们在患者1中看到了一个变体，但在患者2中没有看到该变体，这说明了什么？病人2没有这种变异的等位基因吗？或者，当我们对病人2进行测序时，是否发生了导致测序器无法读取我们感兴趣的变体的错误？联合基因分型以三种不同的方式解决这个问题：结合来自多个样本的证据，我们可以拯救不符合严格统计显著性的变异，以便在单个样本中准确检测出来随着你在人类基因组中每个位点的预测准确性的提高，你能够更好地模拟测序错误和过滤虚假的变异联合基因分型为所有样本提供了一个通用的变体表示，简化了询问个体X中的变体是否也存在于个体Y中的问题联合基因分型在规模上准确识别遗传变异联合基因分型的工作原理是在计算每个个体基因型的可能性时，将所有个体的数据集中在一起。这为我们提供了一个统一的表示每个个体中每个变体有多少个拷贝，这是研究遗传变异与疾病之间联系的关键踏脚石。当我们计算这些新的可能性时，我们也能够计算出在总体中出现的给定变量的先验概率分布，我们可以用它来消除边界变量调用的歧义。对于一个更具体的例子，表1显示了通过GATK变体调用管道调用的样本HG002的indel（插入/删除）和单核苷酸变体（SNV）的精确度和召回率统计数据，与高置信区域中的瓶内基因组（GIAB）高置信度变体调用进行了比较。表1:HG002的变量调用精度，作为单个样本处理召回精密度索引96.25%98.32%SNV公司99.72%99.40%相比之下，表2显示了当我们联合调用HG002和两个亲属（HG003和HG004）的变体时，索引精确度和召回率以及SNV召回率的提高。这个错误率减半是非常重要的，尤其是在临床应用中。表2:HG002与HG003和HG004联合基因分型后的变体调用准确率召回精密度索引98.21%98.98%SNV公司99.78%99.34%最初，联合基因分型是直接从所有个体的原始测序数据中进行的，但是随着研究规模的扩大，这种方法变得不切实际。现代方法从基因组变异调用文件（gVCF）开始，gVCF是一个制表符分隔的文件，包含在单个样本中看到的所有变体，以及在没有发现变异的每个位置上有关测序数据质量的信息。虽然基于gVCF的方法所涉及的数据比查看序列要少，但是一个中等规模的项目仍然可以拥有数十兆字节的gVCF数据。这种基于gVCF的方法消除了返回原始读取的需要，但仍然需要在联合调用变量时重新处理所有N+1个gVCFs，并将新样本添加到队列中。我们的方法使用Delta Lakes来实现队列中N+1个样本的增量平方，同时使用apachesparktm并行化regenotyping。在群体水平上调用变体的挑战尽管联合变异调用很重要，但生物信息学团队经常推迟这一步，因为GATK4周围现有的基础设施使得这些工作负载难以运行，甚至更难以扩展。GATK4的联合基因分型算法的默认实现是单线程的，这种实现的扩展依赖于使用工作流语言和WDL和Cromwell等运行程序手动并行化联合基因分型内核。虽然GATK4支持基于Spark的单倍型调用者，但它不支持运行使用Spark并行化的geneogegvcfs。此外，为了可伸缩性，GATK4最佳实践联合基因分型工作流依赖于在genomisdb中存储数据。不幸的是，GenomicsDB对aws3或azureblob存储等云存储系统的支持有限，研究表明，在一些大型数据集上，新的GenomicsDB工作流比旧的CombineGVCFs/geneogegvcfs工作流慢。我们的解决方案基因组学联合基因分型管道的统一分析平台（Azure | AWS）为这些共同需求提供了解决方案。图2显示了联合基因分型流水线的计算架构。这个管道以笔记本（Azure | AWS）的形式提供，可以作为Databricks作业调用，联合变量调用管道运行起来很简单：用户只需提供他们的输入文件和输出目录。当管道运行时，它首先通过VCF阅读器Azure | AWS）将输入的gVCF数据附加到Delta Lake。Delta-Lake提供了廉价的增量更新，这使得在现有队列中添加N+1个样本非常便宜。当管道运行时，它使用sparksql来存储变量调用。然后，联合变量调用算法在每个bin上并行运行，并随变量的数量线性伸缩。图2:Databricks联合基因分型流水线的计算流程。在第1阶段，gVCF数据被摄取到一个三角洲湖的柱状存储器中，该存储由基因组库划分。随着新样本的到来，这个三角洲湖表可以逐步更新。在第二阶段，我们从Delta-Lake表中加载变量和参考模型，并直接从GATK4的基因型egvcfs工具中运行核心重分型算法。默认情况下，最终平方的基因型矩阵保存到Delta Lake，但也可以写成VCF。并行连接变量调用步骤通过sparksql实现，使用与DNASeq和TNSeq管道类似的架构。具体地说，我们将来自gVCF文件的所有输入基因型/参考模型放入参考基因组的相邻区域。在每个bin中，我们根据参考位置和样本ID对数据进行排序，然后直接从GATK4的GenotypeGVCFs工具中调用联合基因分型算法，并通过排序的迭代器为基因组bin调用。然后我们将这些数据保存到Delta表中，也可以选择作为VCF文件。有关技术实现的更多细节，请参阅我们的Spark Summit 2019 talk。标杆管理为了对我们的方法进行基准测试，我们使用来自1000个基因组项目的低覆盖率WGS数据进行规模测试，并使用来自瓶内基因组联盟的数据进行准确度基准测试。为了生成输入的gVCF文件，我们使用DNASeq管道对齐并调用变量。图3展示了我们的方法在数据集和集群大小上都是有效的可伸缩的。有了这个架构，我们能够在79小时内，在13台c5.9X大型机器上，联合调用来自1000个基因组项目的2504个样本全基因组测序数据的变体。到目前为止，我们已经与客户合作，在超过3000个全基因组样本的项目中扩展这条管道。图3：通过将输入数据保持在10个样本不变并增加执行器的数量来评估强缩放（左）。弱标度（右）通过将簇大小固定为13个i3.8x大工作线程并增加处理的gVCFs数量来评估。在默认设置下运行，该管道与HG002、HG003和HG004三元组上的GATK4联合变量调用管道高度一致。表3描述了当我们的管道与"地面t"比较时，变异和基因型水平的一致性