介绍
本文档描述了如何利用Sentieon®基因组学工具的分片能力将DNAseq®流程分布到多台 服务器上;将其他流程(如TNseq®)进行分布遵循相同原则,因为所有Sentieon®基因组学工具都具有相同的内置分布式处理能力。这种分布的目标是为了减少流程的总运行时间,以更快地生成结果;然而,这种分布也会带来一些额外的开销,使计算成本增加。 利用分布能力,流程的每个阶段被分成小任务;每个任务处理基因组的一部分,并可以在不同的服务器上并行运行。每个任务生成一个部分结果,需要按顺序合并为最终的单一输出;这种合并需要仔细进行,以确保考虑到边界并生成与没有分片运行的流程相同的结果。 分布的执行框架不在本文档的范围内,用户需要在保持正确的数据依赖关系的同时,分发数据/文件并启动正确的进程。
分片和分片化
我们将基因组分成许多连续且不重叠的部分,每个部分称为一个分片(shard)。每个分片可以定义为单个染色体contig的名称,或者是遵循contig:shard_start-shard_end约定的染色体的一部分。特殊的分片NO_COOR用于所有没有坐标的未映射读取。 Sentieon®二进制文件支持将分片分布到多个服务器,并且可以通过添加一个或多个带参数的分片选项在单个命令中处理多个分片。在单个命令行中使用多个选项时,这些分片需要按照参考染色体列表是连续的;例如,一个命令可以包含一个覆盖chr2结尾的分片和一个覆盖chr3开头的分片,但不能同时包含一个覆盖chr2结尾和chr22开头的分片。--shard SHARD
--shard
您可以参考附录中的脚本,根据与参考相关联的dict或输入bam文件中的bam文件表头,为基因组创建分片的示例文件。我们建议使用100M作为分片大小。generate_shards.sh
分布式处理的数据流和数据/文件依赖关系
按照推荐的工作流程,DNAseq®的流程将一对fastq文件通过以下阶段进行处理:BWA对齐生成sorted.bam,去重生成dedup.bam,BQSR生成recal.table,Haplotyper生成output.vcf.gz文件。如图1所示,说明了这样一个流程的数据流。
图1 DNAseq®流程的典型数据流程
为了将上述流程分布到多个服务器上,每个阶段被划分为处理分片数据的命令,这些命令需要来自特定分片以及相邻分片的输入;但是,有两个例外情况:Dedup阶段需要来自所有分片上的LocusCollector命令的所有score文件,而Haplotyper阶段需要一个完整的合并校正表。
图2 以4个分片进行分布式处理的DNAseq®流程的数据流程,不处理未映射的读取
以图2为例,说明了一个以4个分片进行分布的流水线的数据流程;这并不是一个典型的使用案例,因为使用推荐的分片大小为100M碱基会导致需要超过30个分片。在图2的示例中,各个阶段需要以下输入并生成以下输出:
- 分片的LocusCollector阶段(去重1)需要sorted.bam作为输入。该阶段生成一个文件。
i-th
part_deduped.bam$shard_i.score.gz
- 分片的Dedup阶段(去重2)需要sorted.bam输入,以及所有LocusCollector阶段的所有结果。该阶段生成一个文件。
i-th
part_deduped$shard_i.bam
- 分片的QualCal阶段需要文件,以及可用的文件和文件。该阶段生成一个文件。
i-th
part_deduped$shard_i.bam
part_deduped$shard_i+1.bam
part_deduped$shard_i-1.bam
part_recal_data$shard_i.table
- 在QualCal之后,所有文件需要合并为一个单一的校准表文件,用于变异调用。驱动程序和QualCal的选项可以用于执行边界感知合并。
part_recal_data$shard_i.table
--passthru --merge
- 分片的Haplotyper阶段需要文件,以及可用的文件和文件;此外,完整合并的校准表需要作为输入。该阶段生成一个文件。
i-th
part_deduped$shard_i.bam
part_deduped$shard_i+1.bam
part_deduped$shard_i-1.bam
part_output$shard_i.vcf.gz
- 在Haplotyper之后,所有文件需要合并为一个单一的输出VCF文件。驱动程序和Haplotyper的选项用于执行边界感知合并。
part_output$shard_i.vcf.gz
--passthru
--merge
- 如果在流程的任何阶段需要合并的输出bam文件,可以使用util二进制文件执行边界感知合并;需要在命令中添加选项,以便util merge不处理读取,而只是按块复制它们。
--mergemode=10
重要提示:在合并结果时,务必按照正确的顺序依次输入部分结果,否则结果将不正确。
分发BWA
上述说明中未包含有关使用BWA进行分发对齐的任何信息。为了分发BWA对齐,您可以使用Sentieon®工具中提供的工具为输入的FASTQ文件创建索引文件;然后,您可以使用fqidx命令的结果作为BWA mem的输入,在不同的服务器上处理FASTQ文件的特定部分;您需要确保在BWA命令中包含该选项,因为fqidx的输出包含交错的reads在单个输出中。这个流程的结果与在单次运行上执行的结果是相同的。fqidx -p
1.当无法创建FASTQ索引文件(版本201808.02+)时进行BWA分发 如果无法创建FASTQ索引文件,则可以使用带有分数选项的util fqidx工具,并同时使用该选项。此方法将输入的FASTQ文件分割成多个读取块的片段,并从每个片段中提取每个元素以在不同的服务器上进行处理,从起始位置到结束位置。-F m/n -K n m-th m 0 n-1
请注意,此功能仅建议在 文件存储速度足够快以支持每个fqidx进程同时读取输入FASTQ文件的IT环境中使用。这在云环境中很常见,或者在具有高带宽NFS系统的本地集群环境中使用。
大规模队列联合调用需要GVCFtyper做分布式处理
针对包含超过1000个样本的大规模联合调用,我们推荐在GVCFtyper中进行设置。虽然默认的基因分型模式在较小样本集中理论上更准确,但多项式模式在大样本集中同样准确,并且在大量样本中的扩展性更好。--genotype_mode multinomial
在运行大量样本(>3000)的联合队列调用时,建议使用类似上述描述的分发技术进行分发;然而,需要考虑几个挑战:
- 联合调用的总运行时间。
- 分布式GVCFtyper命令将在哪些机器上运行时的资源需求。
- 存储和访问大量GVCF输入文件的物流问题。
- 输出VCF文件的文件大小。
一般推荐使用Sentieon®的GVCFtyper的分片功能,将不同的基因组片段分别处理在不同的机器上。以下命令访问完整的GVCF输入列表,并假设这些输入文件位于共享位置,如NFS存储或可通过s3 http/https
协议进行访问的远程 对象存储位置。
重要提示: 在处理每个片段时,输入GVCF的列表需要保持一致的顺序,因为最终输出中的样本顺序将取决于输入顺序,并且合并需要所有部分文件具有相同的样本顺序。 使用--shard选项时,GVCFtyper生成的输出VCF文件不是有效的VCF文件,因此在下面描述的合并之前不应使用它们。 在处理完所有片段后,您需要运行GVCFtyper的合并命令来合并结果,确保中间VCF输入的顺序与参考基因组一致。输入文件需要在共享位置(例如NFS存储或可通过s3 http/https
协议访问的远程对象存储位置)中可用。
总体运行时间和资源要求
为了减少整体运行时间,建议使用足够的分片,允许将运行时细粒度分布为多个 服务器。分片可以按分片大小或预期创建复杂性。但是,单个分片的最终合并步骤结果无法分发,需要在单个服务器中运行; 此事实设置了用于分布,因为合并可以主导整个运行时。 GVCFtyper 算法是一种非常有效的算法,可能会 I/O 受 GVCF 存储位置性能的限制。因此,建议将 GVCF 输入存储在文件中提供 600 MBps 传输速率的系统。 对于极大量样本 (>10000) 联合队列调用,内存可能会成为一个问题,某些操作系统限制也可能发挥作用。在这种情况下,建议执行以下操作:
- 将操作系统打开文件限制设置为足够大的数量,以允许软件以打开足够的文件句柄。这是通过命令完成的。
ulimit -n NUM
- 将操作系统堆栈限制设置为足够大的数量,以允许软件为操作分配足够的内存,因为内存与样本数,并且可能太大而无法容纳在默认情况下分配的堆栈中。 这是通过命令完成的。
ulimit -s NUM
- 使用包含输入 GVCF 列表的文件,以防止命令过长 参数列出操作系统长度。您可以使用以下命令执行此操作:
- 按照 jemalloc 应用说明中所述使用 jemalloc 内存分配,并通过以下方式更新 vcf 缓存大小 将以下代码添加到脚本中:
- 将分片大小设置为较小的数字,即50M。此外,使用 GVCFtyper 命令中的驱动程序选项,以确保所有线程都得到充分利用。该命令将变为:
--traverse_param
大型输出VCF文件的挑战
当运行非常大的队列时,输出VCF文件将包含大量列,其中包含每个样本的基因型信息。这么多列可能使输出文件变得难以处理,例如在完整文件上运行VQSR可能不切实际,因此您可能需要考虑替代VCF存储输出的方法。 根据您用于存储输出的方法,将输出按样本组或基因组坐标分割可能会改善大型输出文件的问题;请随时与Sentieon支持团队联系,告知您如何存储联合调用的输出的具体要求。 (1)按基因组区域分割输出 为了使输出VCF文件变小,您可以在特定的基因组子区域(例如单个染色体)上执行片段合并。您可以通过仅合并子集的中间VCF文件来实现此目的。例如,您可以通过仅合并涵盖每个染色体的片段来创建每个染色体的VCF文件:如果片段1-4涵盖chr1,而片段5同时涵盖chr1和chr2,则以下代码将创建一个仅包含chr1变异体的VCF文件:
使用上述方法,您可以创建有效的VCF文件,其大小仅为完整VCF文件的一小部分。 为了在上述输出上运行VQSR,您需要提供一个包含完整基因组范围内所有变异记录的VCF文件,因为这些记录需要用于对VQSR模型进行适当的校准。然而,VQSR仅需要VCF数据的前8列,因此您无需将每个特定基因组子区域的所有VCF文件连接起来,可以通过提取和连接每个文件的前8列来创建一个包含必要信息的较小的VCF文件。使用下面的代码将创建所需的文件:VarCal VarCal
接下来,可以在包含VCF的前八列的合并VCF上运行算法,生成用于SNP和INDEL的分段和重新校准文件。然后,可以直接将VQSR应用于按基因组区域分割的VCF,使用算法进行操作:VarCal ApplyVarCal
(2)按样本组分割输出 另外,Sentieon® GVCFtyper合并工具提供了算法选项作为解决这一挑战的潜在方案。在合并步骤中使用算法选项可以生成有效的部分VCF文件,每个文件包含一部分样本,从而将完整的VCF文件分割成较小、更易处理的输出文件。使用算法选项的用法为:--split_by_sample--split_by_sample--split_by_sample
算法选项中使用的输入文件是一个以制表符分隔的文件,每行以特定样本组的输出文件开头,后跟以制表符分隔的相应样本组样本列表。您可以使用多行具有相同输出文件的方式,将多个行中的所有样本分组。下面的示例显示了两个样本组,其中样本1到5将输出到group1输出文件,而样本6到8将输出到group2输出文件:split.conf
--split_by_sample
上述GVCFtyper合并命令将生成以下输出:
- GVCFtyper_main.vcf.gz:一个包含所有VCF信息,包括INFO列在内的部分VCF文件。不包含样本信息。此文件对于运行VQSR非常有用,因为它包含了变异校准所需的所有必要信息。
- GVCFtyper_file_group1.vcf.gz:一个包含至INFO列、FORMAT列以及group1中所有样本的所有列的部分VCF文件。
- GVCFtyper_file_group2.vcf.gz:一个包含至INFO列、FORMAT列以及group2中所有样本的所有列的部分VCF文件。 然后,您可以使用bcftools合并部分VCF并选择您感兴趣的样本。您可以使用以下代码来实现:
在该脚本中,第三个参数是一个逗号分隔的感兴趣样本列表,而extract.sh脚本如下所示:
云环境的特殊考虑事项
在云环境中,通常没有能够容纳大量GVCF输入的NFS存储。对于联合基因分型,Sentieon® GVCFtyper允许托管在对象存储位置(如AWS S3)或通过HTTP访问的GVCF输入文件。然而,对于大型队列(100+),使用对象存储输入的GVCFtyper命令的内存需求可能会过高,因此不建议使用此方法访问输入文件。 在云环境中推荐的方法是根据节点处理的分片,将部分GVCF输入文件下载到计算节点上。可以使用bcftools进行GVCF输入的部分下载,但是需要在bcftools命令中添加--no-version选项,以确保不同分片的标头不会因差异而导致GVCFtyper拒绝合并它们。
为了获得额外的效率,可以使用“瀑布式”方法,如图3和图4所示。在这种方法中,计算节点将按顺序处理多个分片,同时运行下一个分片的部分GVCF输入文件的下载,并与当前分片的GVCFtyper处理并行进行,从而更有效地共享资源,因为这是一个CPU密集型和I/O密集型的过程。流程如下:
- 下载shard1的部分GVCF输入文件。
- 启动shard1的GVCFtyper,并与此同时并行下载shard2的部分GVCF输入文件。
- 在上一步完成之后,启动shard2的GVCFtyper,并与此同时使用bcftools下载shard3。
图3 在云环境中将GVCFtyper分发到多个服务器
图4 在云环境中将GVCFtyper分发到多个服务器,下载下一个分片的输入并与上一个分片的处理并行进行的瀑布式方法的详细说明
Shell示例
以下是用于分发命令的Shell示例,使用1G碱基的分片大小,仅选用于演示目的。推荐的分片大小是100M碱基。