GATK4 多个样本GenotypeGVCFs前用 Combin
我们知道,GATK 4
多个样本joint genotyping用模块GenotypeGVCFs
, 目前GenotypeGVCFs
只支持以下三种形式的输入文件:
1)a single single-sample GVCF
2)a single multi-sample GVCF created by CombineGVCFs
3)a GenomicsDB workspace created by GenomicsDBImport.
即:单个样本的GVCF文件;由CombineGVCFs
模块将多个样本的GVCF文件生成在一起的文件;由GenomicsDB
模块将多个样本GVCF处理生成一起的工作空间。当然这里的GVCF 文件是由HaplotypeCaller
模块的-ERC GVCF
或者 -ERC BP_RESOLUTION
参数产生, 如果是其他工具生成的GVCF可能会因为缺少某些GenotypeGVCFs
需要的重要信息导致出错。
因此对于多个样本joint genotyping,在用GenotypeGVCFs前需要将多个样本的gvcf文件用CombineGVCFs
方式或者GenomicsDBImport
方式合并成一个文件,前者是一个总的gvcf文件,后者是一个GenomicsDB工作目录。
按照GATK的官方说明,CombineGVCFs
效率比较低,需要许多内存,推荐使用 GenomicsDBImport
。不过需注意,GenomicsDBImport
是需要分开interval计算的,如分染色体计算,其每次只能处理一个interval。鉴于GATK极力推荐GenomicsDBImport
,我们以染色体chr10
为例测试CombineGVCFs
和GenomicsDBImport
对一个trio家系的外显子数据效果,这两个模块的命令分别如下:
CombineGVCFs:
java -Xmx4g -jar gatk-package-4.1.2.0-local.jar CombineGVCFs -R GRCh38.fa -L chr10.bed --variant father_chr10.g.vcf.gz --variant mother_chr10.g.vcf.gz --variant child_chr10.g.vcf.gz -O family_chr10.g.vcf.gz
GenomicsDBImport:
java -Xmx4g -Xms4g -jar gatk-package-4.1.2.0-local.jar GenomicsDBImport -R GRCh38.fa -L chr10.bed --variant father_chr10.g.vcf.gz --variant mother_chr10.g.vcf.gz --variant child_chr10.g.vcf.gz --genomicsdb-workspace-path chr10_database.db
# -L是必须参数,上文已提到;
# --genomicsdb-workspace-path后面接的必须是一个不存在的空目录;
分别运行上述两个命令,但是问题来了...,发现CombineGVCFs
不到5分钟即可完成,而GenomicsDBImport
已经超过24个小时还在运行中,而且在其genomicsdb-workspace-path产生的chr10_database.db目录下生成超过185,752个子文件,并还在继续生成中!!!这只是对一个染色体的运行,产生如此多的子文件对于有文件数目limit的存储是一个很大的限制,更重要的是运行时间如此之久。这是为什么呢?
首先检查了上述命令是没有问题的,Google粗略搜索了下没找到原因,直接去GATK 论坛上询问,给出的解释如下,比较简单直接:
GenomicsDBImport is used for samples in the order of thousands. For <1000 samples it is better to use CombineGVCFs.
也就是说 GenomicsDBImport
更适用于1000个样本以上的joint genotyping!好吧,这点在GATK的官方使用文档中并没有说明。带着这个问题的疑虑,我又搜索了下发现其实先前已有很多人问过相同的问题并在GATK论坛上深入讨论过,大体总结如下: