Hi GATK Users,

Happy Thanksgiving!
Our staff will be observing the holiday and will be unavailable from 22nd to 25th November. This will cause a delay in reaching out to you and answering your questions immediately. Rest assured we will get back to it on Monday November 26th. We are grateful for your support and patience.
Have a great holiday everyone!!!

Regards
GATK Staff

MuTect2 beta --germline_resource for build b37

Jeff_GaitherJeff_Gaither Math Biosciences Institute, Columbus, OHMember

Hi - I'm looking to run MuTect2 beta using the --germline_resource option. However, I've consistently used the b37 genome build throughout my analaysis, while the suggested resource (gnomad) appears to only be available for the hg19 build. So I'm wondering whether I should
1. Go ahead and use the gnomad hg19 files, despite the fact that my whole analysis has used the b37 build?
2. Lift over my existing gnomad vcfs from hg19 to b37? (In this case, I'd need an hg19tob37 liftOver file - I can't find one anywhere).
3. Use another germline resource?

I wonder which option you would recommend? Many thanks for your time.

Best Answer

Answers

  • Jeff_GaitherJeff_Gaither Math Biosciences Institute, Columbus, OHMember
    edited July 2017

    Okay, I found an hg19 -> b37 liftover chain file at
    ftp://[email protected]/Liftover_Chain_Files

    However, when I tried to lift the gnomad vcf files from hg19 to b37 using Picard's LiftoverVcfs, every single site was rejected as unliftable. So this didn't work.

    Post edited by Jeff_Gaither on
  • Jeff_GaitherJeff_Gaither Math Biosciences Institute, Columbus, OHMember
    Accepted Answer

    Hey, I just found the needed file in the Resource bundle...
    ftp://[email protected]/bundle/beta/Mutect2/af-only-gnomad.raw.sites.b37.vcf.gz
    So, problem solved!

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Jeff_Gaither
    Hi,

    Glad you solved it on your own!

    -Sheila

  • Hi, could you give me af-only-gnomd.vcf for the hg19 build,I don't know how to find this vcf.
    Thanks~

  • KHL0702KHL0702 Member
    edited August 2017

    Hi, I would like to know where I can download the file af-only-gnomd.vcf for the hg38 build,the Resource bundle you given is only for b37,
    I hope to get some advice from you,thanks for your help!

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @luoshizhi @KHL0702
    Hi,

    Have a look at this thread. It is the same for hg38 too.

    -Sheila

  • Binfls6321DTUBinfls6321DTU Member
    edited October 2017

    @Sheila
    Hi,
    I'm new in using MuTect2 :p
    I'm a bit confused that the gnomAD file provided on gnomad.broadinstitute.org/downloads is about three times larger than the gnomAD file provided on ftp://[email protected]/bundle/beta/Mutect2/af-only-gnomad.raw.sites.b37.vcf.gz ?
    is there any difference between these two files ? could the difference affect the output of MuTect2??

    Thx in advance :)

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Binfls6321DTU
    Hi,

    Sorry for the confusion. That file is what we use in the Mutect2 tutorial. It only contains variants from one chromosome. You should use the entire gnomAD VCF if you are running on all chromosomes. I will make a note to make it more clear in the documentation for GATK4.

    Also, if you are new, have a look at the Mutect2 tutorials in the GATK presentations section. I would recommend running through the tutorial from July for a basic introduction to the workflow, then trying out the September version for a more in-depth look :smiley:

    -Sheila

  • Binfls6321DTUBinfls6321DTU Member
    edited October 2017

    @Sheila Thank you very much, :D I have tried the gnomAD file and succeeded. Note that the some AF values in gnomAD vcf file are only a dot, e.g. ".", and that will make the program crash.

    Another problem I get now is, it always reports the "segmentation fault" when I use the samtools-1.5/misc/seq_cache_populate.pl

    /var/spool/torque/mom_priv/jobs/18460427.risoe-r04-sn064.cm.cluster.SC: line 20: 13995 Segmentation fault /home/projects/cu_10098/data/ennbin/project/hg37_test/origin/copy/cromwell-executions/WES_normal_tumor_somatic_SNV_wf/aef23f13-da39-4c48-bc7f-1fe70c87d06a/call-ConvertToCram_tumor/inputs/home/projects/cu_10098/apps/src/samtools-1.3.1/samtools view -C -T /home/projects/cu_10098/data/ennbin/project/hg37_test/origin/copy/cromwell-executions/WES_normal_tumor_somatic_SNV_wf/aef23f13-da39-4c48-bc7f-1fe70c87d06a/call-ConvertToCram_tumor/inputs/home/databases/gatk-legacy-bundles/b37/human_g1k_v37_decoy.fasta /home/projects/cu_10098/data/ennbin/project/hg37_test/origin/copy/cromwell-executions/WES_normal_tumor_somatic_SNV_wf/aef23f13-da39-4c48-bc7f-1fe70c87d06a/call-ConvertToCram_tumor/inputs/home/projects/cu_10098/data/ennbin/project/hg37_test/origin/copy/cromwell-executions/WES_normal_tumor_somatic_SNV_wf/aef23f13-da39-4c48-bc7f-1fe70c87d06a/call-GatherBamFiles_tumor/execution/TCRBOA3-T.bam
    13996 Done | tee TCRBOA3-T.cram
    13997 Done | md5sum
    13998 Done | awk '{print $1}' > TCRBOA3-T.cram.md5

    below is the command

    /home/projects/cu_10098/data/ennbin/project/hg37_test/origin/copy/cromwell-executions/WES_normal_tumor_somatic_SNV_wf/aef23f13-da39-4c48-bc7f-1fe70c87d06a/call-ConvertToCram_tumor/inputs/home/projects/cu_10098/apps/src/samtools-1.3.1/samtools view -C -T /home/projects/cu_10098/data/ennbin/project/hg37_test/origin/copy/cromwell-executions/WES_normal_tumor_somatic_SNV_wf/aef23f13-da39-4c48-bc7f-1fe70c87d06a/call-ConvertToCram_tumor/inputs/home/databases/gatk-legacy-bundles/b37/human_g1k_v37_decoy.fasta /home/projects/cu_10098/data/ennbin/project/hg37_test/origin/copy/cromwell-executions/WES_normal_tumor_somatic_SNV_wf/aef23f13-da39-4c48-bc7f-1fe70c87d06a/call-ConvertToCram_tumor/inputs/home/projects/cu_10098/data/ennbin/project/hg37_test/origin/copy/cromwell-executions/WES_normal_tumor_somatic_SNV_wf/aef23f13-da39-4c48-bc7f-1fe70c87d06a/call-GatherBamFiles_tumor/execution/TCRBOA3-T.bam | \
    tee TCRBOA3-T.cram | \
    md5sum | awk '{print $1}' > TCRBOA3-T.cram.md5
    #Create REF_CACHE. Used when indexing a CRAM
    /home/projects/cu_10098/data/ennbin/project/hg37_test/origin/copy/cromwell-executions/WES_normal_tumor_somatic_SNV_wf/aef23f13-da39-4c48-bc7f-1fe70c87d06a/call-ConvertToCram_tumor/inputs/home/projects/cu_10098/apps/src/samtools-1.5/misc/seq_cache_populate.pl -root ./ref/cache /home/projects/cu_10098/data/ennbin/project/hg37_test/origin/copy/cromwell-executions/WES_normal_tumor_somatic_SNV_wf/aef23f13-da39-4c48-bc7f-1fe70c87d06a/call-ConvertToCram_tumor/inputs/home/databases/gatk-legacy-bundles/b37/human_g1k_v37_decoy.fasta
    export REF_PATH=:
    export REF_CACHE=./ref/cache/%2s/%2s/%s
    /home/projects/cu_10098/data/ennbin/project/hg37_test/origin/copy/cromwell-executions/WES_normal_tumor_somatic_SNV_wf/aef23f13-da39-4c48-bc7f-1fe70c87d06a/call-ConvertToCram_tumor/inputs/home/projects/cu_10098/apps/src/samtools-1.3.1/samtools index TCRBOA3-T.cram
    mv TCRBOA3-T.cram.crai TCRBOA3-T.crai

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @Binfls6321DTU
    Hi,

    Note that the some AF values in gnomAD vcf file are only a dot, e.g. ".", and that will make the program crash.

    What did you do to get around that? Can you post the exact error message?

    For Samtools, I would recommend asking on their maliling list, as we don't help with their tools. The only thing I can recommend is trying the very latest version of their tools :smile:

    -Sheila

  • Hi,
    I'm having the same problem with dots in the gnomad file when passing it through the germline_resource option. Here is the error message:

    java.lang.NumberFormatException: For input string: "."
    at sun.misc.FloatingDecimal.readJavaFormatString(FloatingDecimal.java:2043)
    at sun.misc.FloatingDecimal.parseDouble(FloatingDecimal.java:110)
    at java.lang.Double.parseDouble(Double.java:538)
    at java.lang.Double.valueOf(Double.java:502)
    at htsjdk.variant.variantcontext.CommonInfo.lambda$getAttributeAsDoubleList$2(CommonInfo.java:299)
    at java.util.stream.ReferencePipeline$3$1.accept(ReferencePipeline.java:193)
    at java.util.Collections$2.tryAdvance(Collections.java:4717)
    at java.util.Collections$2.forEachRemaining(Collections.java:4725)
    at java.util.stream.AbstractPipeline.copyInto(AbstractPipeline.java:481)
    at java.util.stream.AbstractPipeline.wrapAndCopyInto(AbstractPipeline.java:471)
    at java.util.stream.ReduceOps$ReduceOp.evaluateSequential(ReduceOps.java:708)
    at java.util.stream.AbstractPipeline.evaluate(AbstractPipeline.java:234)
    at java.util.stream.ReferencePipeline.collect(ReferencePipeline.java:499)
    at htsjdk.variant.variantcontext.CommonInfo.getAttributeAsList(CommonInfo.java:273)
    at htsjdk.variant.variantcontext.CommonInfo.getAttributeAsDoubleList(CommonInfo.java:293)
    at htsjdk.variant.variantcontext.VariantContext.getAttributeAsDoubleList(VariantContext.java:740)
    at org.broadinstitute.hellbender.tools.walkers.mutect.Mutect2Engine.isActive(Mutect2Engine.java:312)
    at org.broadinstitute.hellbender.engine.AssemblyRegionIterator.loadNextAssemblyRegion(AssemblyRegionIterator.java:158)
    at org.broadinstitute.hellbender.engine.AssemblyRegionIterator.next(AssemblyRegionIterator.java:129)
    at org.broadinstitute.hellbender.engine.AssemblyRegionIterator.next(AssemblyRegionIterator.java:32)
    at org.broadinstitute.hellbender.engine.AssemblyRegionWalker.processReadShard(AssemblyRegionWalker.java:250)
    at org.broadinstitute.hellbender.engine.AssemblyRegionWalker.traverse(AssemblyRegionWalker.java:231)
    at org.broadinstitute.hellbender.engine.GATKTool.doWork(GATKTool.java:838)
    at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:119)
    at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:176)
    at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:195)
    at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:131)
    at org.broadinstitute.hellbender.Main.mainEntry(Main.java:152)
    at org.broadinstitute.hellbender.Main.main(Main.java:233)

    Many thanks in advance for your help,

    Olivier

    Issue · Github
    by shlee

    Issue Number
    2654
    State
    closed
    Last Updated
    Closed By
    chandrans
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @olisand
    Hi Olivier,

    Can you check if this issue still occurs in GATK4 latest beta? Also, if the VCF validates with ValidateVariants and still errors in GATK4, can you submit a bug report?

    Thanks,
    Sheila

  • Hi Sheila,

    Thank you for your answer.
    It is still occurring with beta.6 and the VCF validates. I'll submit a bug report today.

    Olivier

  • shleeshlee CambridgeMember, Broadie, Moderator admin
    edited November 2017

    Hi Olivier (@olisand),

    I am following up on behalf of the team who are away for a workshop. Let me know when you submit your bug report.

    Just to be clear I understand, the gnomad file that is causing the error:

    java.lang.NumberFormatException: For input string: "."
    at sun.misc.FloatingDecimal.readJavaFormatString(FloatingDecimal.java:2043)

    is from ftp://[email protected]/bundle/beta/Mutect2/af-only-gnomad.raw.sites.b37.vcf.gz? Or does it come from a different source?

    Thanks.

  • Hi (@shlee),

    I got it from there: https://storage.googleapis.com/gnomad-public/release/2.0.2/vcf/exomes/gnomad.exomes.r2.0.2.sites.vcf.bgz
    and added 'chr' before chromosome numbers in it.
    I just tried with the one from the bundle (after adding 'chr' to each chromosome number) and it worked without error. So, I guess case is solved ?

    Thank you for your help,

    Olivier

  • shleeshlee CambridgeMember, Broadie, Moderator admin

    Hi @olisand,

    The URL you give isn't one of ours so I cannot comment. I'm glad you've solved your issue.

    The resources we guarantee work with our tools come from our Resource Bundle at https://software.broadinstitute.org/gatk/download/bundle and have already been formatted for either GRCh37 (sans chr) or hg19 (with chr) contig naming conventions. I recommend you use the files we provide for the convention that you need instead of converting between the two. The concern is there are some coordinate differences for the mitochondria. You can also use Picard LiftoverVcf with CHAIN files as well.

  • Hi @shlee,

    Ok, thanks again for your help.

  • kkshaxqdkkshaxqd chinaMember

    Hello @olisand, I download the gnomad.exomes.r2.0.2.sites.vcf.bgz ,but can not deal with the bgz file, so I want to know how to add chr in the gnomad.exomes.r2.0.2.sites.vcf.bgz, and let the mutect2 use the file....

    Thanks.

    JOSH ZHANG

  • shleeshlee CambridgeMember, Broadie, Moderator admin

    Hi @kkshaxqd,

    We provide a simplified gnomad resource file (for GRCh37) in the FTP Bundle, in beta>Mutect2.

  • kkshaxqdkkshaxqd chinaMember

    Hello, @shlee , thanks very much , en....That's right, I think it will work and resolve the problem. GRCh37/hg19.....

  • @olisand said:
    Hi (@shlee),

    I got it from there: https://storage.googleapis.com/gnomad-public/release/2.0.2/vcf/exomes/gnomad.exomes.r2.0.2.sites.vcf.bgz
    and added 'chr' before chromosome numbers in it.
    I just tried with the one from the bundle (after adding 'chr' to each chromosome number) and it worked without error. So, I guess case is solved ?

    Thank you for your help,

    Olivier

    Hi,Olisand!
    I have run into the same question as yours, and i wonder how did you make it that "add 'chr' before chromosome numbers in it"? and would you like to share the experience with me? thanks a lot!

  • huseehusee Member

    I can not open the site: ftp://[email protected]/Liftover_Chain_Files,where can I get a chain file to lift small_exac_common_3_b37.vcf.gz to hg19?

  • shleeshlee CambridgeMember, Broadie, Moderator admin
  • huseehusee Member

    Thanks ,@shlee ,I can open the ftp://[email protected]/Liftover_Chain_Files suddenly and download the b37tohg19.chain file.

  • shleeshlee CambridgeMember, Broadie, Moderator admin

    Great to hear you can now access the FTP folder @husee.

Sign In or Register to comment.