Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.

GenomicsDBImport - Duplicate field name MQ0 found in vid map

sf21sf21 Member
edited October 2018 in Ask the GATK team

I am running into an issue running GenomicsDBImport (version v4.0.5.1) on a single gVCF. The gVCF was generated using gatk version, 3.5-0-g36282e4.

Here is the error,

...
15:22:35.273 DEBUG GenomeLocParser -  HLA-DRB1*16:02:01 (11005 bp)
15:22:35.626 INFO  IntervalArgumentCollection - Processing 64444167 bp from intervals
15:22:35.751 INFO  GenomicsDBImport - Done initializing engine
15:22:35.752 DEBUG IOUtils - Deleted /my_database_chr20
Created workspace /my_database_chr20
15:22:35.941 INFO  GenomicsDBImport - Vid Map JSON file will be written to /my_database_chr20/vidmap.json
15:22:35.941 INFO  GenomicsDBImport - Callset Map JSON file will be written to /my_database_chr20/callset.json
15:22:35.941 INFO  GenomicsDBImport - Complete VCF Header will be written to /my_database_chr20/vcfheader.vcf
15:22:35.941 INFO  GenomicsDBImport - Importing to array - /my_database_chr20/genomicsdb_array
15:22:35.955 INFO  ProgressMeter - Starting traversal
15:22:35.955 INFO  ProgressMeter -        Current Locus  Elapsed Minutes     Batches Processed   Batches/Minute
15:22:35.955 INFO  GenomicsDBImport - Starting batch input file preload
15:22:36.161 INFO  GenomicsDBImport - Finished batch preload
15:22:36.161 INFO  GenomicsDBImport - Importing batch 1 with 1 samples
Duplicate field name MQ0 found in vid map
terminate called after throwing an instance of 'ProtoBufBasedVidMapperException'
  what():  ProtoBufBasedVidMapperException : Duplicate fields exist in vid map

GenomicsDBImport command,

JAVA_OPTS="-Xmx${memory}g" ${gatk_exe} GenomicsDBImport \
  --genomicsdb-workspace-path ${run_dir}/my_database_chr20 \
  --intervals chr20:1-100000 \
  --verbosity DEBUG \
  --lenient \
  --reader-threads 5 \
  -ip 500 \
  --batch-size 1  \
  --sample-name-map ${data_dir}/samples.list \
  --overwrite-existing-genomicsdb-workspace \
  -R ${ref}

Haplotype caller command,

/Software/bin/java
 -XX:ParallelGCThreads=2
 -Djava.io.tmpdir=/scratch/68a3f0e84d564e729fc42934118dc642
 -Xmx8196M
 -jar /Software/GenomeAnalysisTK/GenomeAnalysisTK-3.5/GenomeAnalysisTK.jar
 -T HaplotypeCaller
 --genotyping_mode DISCOVERY
 -A AlleleBalanceBySample
 -A DepthPerAlleleBySample
 -A DepthPerSampleHC
 -A InbreedingCoeff
 -A MappingQualityZeroBySample
 -A StrandBiasBySample
 -A Coverage
 -A FisherStrand
 -A HaplotypeScore
 -A MappingQualityRankSumTest
 -A MappingQualityZero
 -A QualByDepth
 -A RMSMappingQuality
 -A ReadPosRankSumTest
 -A VariantType
 -l INFO
 --emitRefConfidence GVCF
 -rf BadCigar
 --variant_index_parameter 128000
 --variant_index_type LINEAR
 -R /Resources/GRCh38_full_analysis_set_plus_decoy_hla.fa
 -nct 1
 -I /data/NA12878.final.bam
 -o /data/NA12878.20.haplotypeCalls.raw.g.vcf
 -L chr20

Don't think this should be the issue, but FYI

$ bcftools view -h ../data/NA12878.haplotypeCalls.raw.g.vcf.gz | grep MQ0

##FORMAT=<ID=MQ0,Number=1,Type=Integer,Description="Number of Mapping Quality Zero Reads per sample">
##INFO=<ID=MQ0,Number=1,Type=Integer,Description="Total Mapping Quality Zero Reads">
Post edited by sf21 on

Answers

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @sf1

    1) I am looking into this issue.

    2) It is not useful to apply MappingQualityZeroBySample/ MappingQualityZero annotation with HaplotypeCaller because HC filters out all reads with MQ0 upfront, so the annotation will always return a value of 0. So maybe try it without those annotations for the time being and I will get back to you about this issue.

    Regards
    Bhanu

  • sf21sf21 Member

    Hi Bhanu,
    Thank you for looking into the issue. For this test case i can definitely re-run haplotype caller to regenerate the GVCF without those annotations and try importing. We have a lot of GVCFs (1000s) that already include those annotations and re-running haplotype caller would be a significant amount of work and so I was wondering if there is a way to either remove those annotations or have GenomicsDBImport ignore them? Thanks!

  • sf21sf21 Member

    I ran a couple of tests and it looks like the issue is connected to how the definition for MQ0 in INFO and FORMAT columns is interpreted from the header.
    Test 1: Re-ran haplotype caller without the MappingQualityZeroBySample and MappingQualityZero annotations and i was able to successfully import the resulting GVCF.
    Test 2: Removed the definition entry for MQ0 field in the FORMAT column from the header and was able to successfully import the GVCF.

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    @sf21

    Thats great. Looks like that should help with the issue for now while we look into this with our developers.

    Regards
    Bhanu

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @sf21

    This was an issue in gatk3 but in gatk4 it has been resolved by only using MappingQualityZero because with HaplotypeCaller because HC filters out all reads with MQ0 upfront, so the annotation will always return a value of 0.

    Regards
    Bhanu

  • sf21sf21 Member

    Got it. Thank you for clarifying.

  • lelimatlelimat San Francisco, CAMember

    Any updates on this? Is there an easier way to run "GenomicsDBImport" without having to rerun HaplotypeCaller or changing the header of the gVCF?

    Thanks!

Sign In or Register to comment.