生信分析工具重测序群体比较

GATK4 多个样本GenotypeGVCFs前用 Combin

2019-09-30  本文已影响0人  danria

我们知道,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为例测试CombineGVCFsGenomicsDBImport对一个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论坛上深入讨论过,大体总结如下:

上一篇下一篇

猜你喜欢

热点阅读