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论坛上深入讨论过，大体总结如下：

• CombineGVCFs对于一千以上样本量的joint genotyping是一个很耗时间与内存的过程， 不建议使用，而GenomicsDBImport则是一个很好的选择。(讨论1)
• 如果使用GenomicsDBImport，intervals建议一定要分割很小，即-L参数后的interval切分成更小（比如有提议推荐 1 interval per sample，即多少份样本多少份intervals)，同时处理多个样本可以用 --batchSize参数指定同时处理的样本数目；(讨论2；讨论3)。