GATK4 多个样本GenotypeGVCFs前用 CombineGVCFs还是GenomicsDBImport

我们知道,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论坛上深入讨论过,大体总结如下:

  • CombineGVCFs对于一千以上样本量的joint genotyping是一个很耗时间与内存的过程, 不建议使用,而GenomicsDBImport则是一个很好的选择。(讨论1)
  • 如果使用GenomicsDBImport,intervals建议一定要分割很小,即-L参数后的interval切分成更小(比如有提议推荐 1 interval per sample,即多少份样本多少份intervals),同时处理多个样本可以用 --batchSize参数指定同时处理的样本数目;(讨论2讨论3)。
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 194,242评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,769评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,484评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,133评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,007评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,080评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,496评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,190评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,464评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,549评论 2 309
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,330评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,205评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,567评论 3 298
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,889评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,160评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,475评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,650评论 2 335

推荐阅读更多精彩内容