We've moved!
This site is now read-only. You can find our new documentation site and support forum for posting questions here.
Be sure to read our welcome blog!

Select variants says VCF header is missing

I'm running select variants, to select one sample from a multi sample vcf with the following command
java -Xmx2g -jar $gatk_dir/GenomeAnalysisTK.jar -R $reference_dir -T SelectVariants --variant $multi_vcf -sn OKJ042_HC2012 -o "OKJ042_unfiltered_vcf"

I get the following error (pasted below) saying a key isn't defined in the VCFHeader, but when I examine the the vcf header, I see the key is defined in the filter ID (see below). The same file will process if I add the commands --excludeFiltered and --excludeNonVariants

All processing has been done within GaTK, except one step to pull out indels and snps into two files which was done in vcftools.

Any advice on what is wrong with the header would be greatly appreciated.

ERROR A GATK RUNTIME ERROR has occurred (version 2.4-7-g5e89f01):
ERROR
ERROR Please visit the wiki to see if this is a known problem
ERROR If not, please post the error, with stack trace, to the GATK forum
ERROR Visit our website and forum for extensive documentation and answers to
ERROR commonly asked questions http://www.broadinstitute.org/gatk
ERROR
ERROR MESSAGE: Key RMSMapppingQuality<40 found in VariantContext field FILTER at Y_unplaced:13251 but this key isn't defined in the VCFHeader. We require all VCFs to have complete VCF headers by default.

First 30 lines of the multi_vcf

fileformat=VCFv4.1

FILTER=<ID="ReadPos < -8.0",Description="ReadPosRankSum < -8.0">

FILTER=<ID=FisherStrand>60,Description="FS > 60.0">

FILTER=<ID=HapScore>13.0,Description="HaplotypeScore > 13.0">

FILTER=<ID=Low<10Coverage,Description="DP < 10">

FILTER=<ID=LowQual,Description="Low quality">

FILTER=<ID=MQRankSum<-12.5,Description="MQRankSum < -12.5">

FILTER=<ID=Q<30,Description="QUAL < 30">

FILTER=<ID=QualityByDepth<2.0,Description="QD < 2.0">

FILTER=<ID=RMSMapppingQuality<40,Description="MQ < 40.0">

FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">

FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">

FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">

FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">

INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">

INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">

INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">

INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities">

INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">

INFO=<ID=DS,Number=0,Type=Flag,Description="Were any of the samples downsampled?">

INFO=<ID=Dels,Number=1,Type=Float,Description="Fraction of Reads Containing Spanning Deletions">

INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">

INFO=<ID=HaplotypeScore,Number=1,Type=Float,Description="Consistency of the site with at most two segregating haplotypes">

INFO=<ID=InbreedingCoeff,Number=1,Type=Float,Description="Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hardy-Weinberg expectation">

INFO=<ID=MLEAC,Number=A,Type=Integer,Description="Maximum likelihood expectation (MLE) for the allele counts (not necessarily the same as the AC), for each ALT allele, in the same order as listed">

INFO=<ID=MLEAF,Number=A,Type=Float,Description="Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT allele, in the same order as listed">

INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">

INFO=<ID=MQ0,Number=1,Type=Integer,Description="Total Mapping Quality Zero Reads">

INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">

INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">

INFO=<ID=RPA,Number=.,Type=Integer,Description="Number of times tandem repeat unit is repeated, for each allele (including reference)">

INFO=<ID=RU,Number=1,Type=String,Description="Tandem repeat unit (bases)">

INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">

INFO=<ID=STR,Number=0,Type=Flag,Description="Variant is a short tandem repeat">

Best Answer

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Will add in a check for that, thanks for pointing it out.

  • prepagamprepagam Member

    Thanks. I removed the < > from within the filterName, and all is ok now.

  • baescbaesc Member

    Hi Dear GATK Team,
    I've got a similar problem, only I haven't got angle brackets.
    Here's my code:
    RIGHT_NOW=$(date +"%x %r %Z") echo "-----------------------------------------------------------------" echo "FUNCTION: BeagleOutputToVCF" ${RIGHT_NOW} echo "-----------------------------------------------------------------" $JAVA -Xmx2g -jar $ProgDir$GATK/GenomeAnalysisTK.jar \ -L ${region} \ -T BeagleOutputToVCF \ -R ${Dir}UMD31.fasta \ -V ./${Reg}/${Reg}.${bull}.flt.vcf \ -beagleR2:BEAGLE ./${Reg}/${Reg}.${bull}.beagle.in.r2 \ -beaglePhased:BEAGLE ./${Reg}/${Reg}.${bull}.beagle.in.phased \ -beagleProbs:BEAGLE ./${Reg}/${Reg}.${bull}.beagle.in.gprobs \ -o ./${Reg}/${Reg}.${bull}.FINAL.vcf done

    ...here's my error:

    `##### ERROR ------------------------------------------------------------------------------------------

    ERROR stack trace

    java.lang.IllegalStateException: Key AC found in VariantContext field INFO at Chr1:4000371 but this key isn't defined in the VCFHeader. We require all VCFs to have complete VCF headers by default.
    at org.broadinstitute.variant.variantcontext.writer.VCFWriter.fieldIsMissingFromHeaderError(VCFWriter.java:601)
    at org.broadinstitute.variant.variantcontext.writer.VCFWriter.add(VCFWriter.java:266)
    at org.broadinstitute.sting.gatk.io.storage.VariantContextWriterStorage.add(VariantContextWriterStorage.java:175)
    at org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub.add(VariantContextWriterStub.java:242)
    at org.broadinstitute.sting.gatk.walkers.beagle.BeagleOutputToVCF.map(BeagleOutputToVCF.java:359)
    at org.broadinstitute.sting.gatk.walkers.beagle.BeagleOutputToVCF.map(BeagleOutputToVCF.java:78)
    at org.broadinstitute.sting.gatk.traversals.TraverseLociNano$TraverseLociMap.apply(TraverseLociNano.java:267)
    at org.broadinstitute.sting.gatk.traversals.TraverseLociNano$TraverseLociMap.apply(TraverseLociNano.java:255)
    at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.executeSingleThreaded(NanoScheduler.java:274)
    at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.execute(NanoScheduler.java:245)
    at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:144)
    at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:92)
    at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:48)
    at org.broadinstitute.sting.gatk.executive.LinearMicroScheduler.execute(LinearMicroScheduler.java:99)
    at org.broadinstitute.sting.gatk.GenomeAnalysisEngine.execute(GenomeAnalysisEngine.java:311)
    at org.broadinstitute.sting.gatk.CommandLineExecutable.execute(CommandLineExecutable.java:113)
    at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:245)
    at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:152)
    at org.broadinstitute.sting.gatk.CommandLineGATK.main(CommandLineGATK.java:91)

    ERROR ------------------------------------------------------------------------------------------
    ERROR A GATK RUNTIME ERROR has occurred (version 2.6-3-gdee51c4):
    ERROR
    ERROR Please check the documentation guide to see if this is a known problem
    ERROR If not, please post the error, with stack trace, to the GATK forum
    ERROR Visit our website and forum for extensive documentation and answers to
    ERROR commonly asked questions http://www.broadinstitute.org/gatk
    ERROR
    ERROR MESSAGE: Key AC found in VariantContext field INFO at Chr1:4000371 but this key isn't defined in the VCFHeader. We require all VCFs to have complete VCF headers by default.
    ERROR ------------------------------------------------------------------------------------------`

    ... and here's my header:

    `sh-3.2$ more SAMTOOLS_BTA1_4000k-4100k_SING.55069.flt.vcf

    fileformat=VCFv4.1

    samtoolsVersion=0.1.18 (r982:295)

    INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">

    INFO=<ID=DP4,Number=4,Type=Integer,Description="# high-quality ref-forward bases, ref-reverse, alt-forward and alt-reverse bases">

    INFO=<ID=MQ,Number=1,Type=Integer,Description="Root-mean-square mapping quality of covering reads">

    INFO=<ID=FQ,Number=1,Type=Float,Description="Phred probability of all samples being the same">

    INFO=<ID=AF1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele frequency (assuming HWE)">

    INFO=<ID=AC1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele count (no HWE assumption)">

    INFO=<ID=G3,Number=3,Type=Float,Description="ML estimate of genotype frequencies">

    INFO=<ID=HWE,Number=1,Type=Float,Description="Chi^2 based HWE test P-value based on G3">

    INFO=<ID=CLR,Number=1,Type=Integer,Description="Log ratio of genotype likelihoods with and without the constraint">

    INFO=<ID=UGT,Number=1,Type=String,Description="The most probable unconstrained genotype configuration in the trio">

    INFO=<ID=CGT,Number=1,Type=String,Description="The most probable constrained genotype configuration in the trio">

    INFO=<ID=PV4,Number=4,Type=Float,Description="P-values for strand bias, baseQ bias, mapQ bias and tail distance bias">

    INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">

    INFO=<ID=PC2,Number=2,Type=Integer,Description="Phred probability of the nonRef allele frequency in group1 samples being larger (,smaller) than in group2.">

    INFO=<ID=PCHI2,Number=1,Type=Float,Description="Posterior weighted chi^2 P-value for testing the association between group1 and group2 samples.">

    INFO=<ID=QCHI2,Number=1,Type=Integer,Description="Phred scaled PCHI2.">

    INFO=<ID=PR,Number=1,Type=Integer,Description="# permutations yielding a smaller PCHI2.">

    INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias">

    FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

    FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">

    FORMAT=<ID=GL,Number=3,Type=Float,Description="Likelihoods for RR,RA,AA genotypes (R=ref,A=alt)">

    FORMAT=<ID=DP,Number=1,Type=Integer,Description="# high-quality bases">

    FORMAT=<ID=SP,Number=1,Type=Integer,Description="Phred-scaled strand bias P-value">

    FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">

    CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 55069

    Chr1 4000371 . C T 130 . DP=11;VDB=0.0380;AF1=1;AC1=2;DP4=0,0,4,5;MQ=26;FQ=-54 GT:PL:DP:SP:GQ 1/1:163,27,0:9:0:51
    Chr1 4000375 . A C 140 . DP=12;VDB=0.0398;AF1=1;AC1=2;DP4=0,0,5,5;MQ=26;FQ=-57 GT:PL:DP:SP:GQ 1/1:173,30,0:10:0:57
    `

    So, the GATK is telling me that the Key "AC" found in the INFO part of my VCF isn't included in my Header, but the problem seems to be that a part of the Key (namely, the "1") is not being read.... what am I doing wrong? I used exaclty the same .bam file in the ProduceBeagleInput walker... is that my mistake?

    Thanks for any help. Cheers, Chris

  • pdexheimerpdexheimer Member ✭✭✭✭

    What happens if you just add an AC into your header? I think what's happening here is that BeagleOutputToVCF copies your existing header and adds a couple of new fields (but not AC). Then it generates an AC INFO field for you, and triggers a validator at some point that throws your exception. I'm specifically looking at the initialize() method and line 344 of BeagleOutputToVCF.java (VariantContextUtils.calculateChromosomeCounts(builder, false);)

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    If that is indeed what's happening it's a bug; we should be generating fully valid VCFs by our own standards.

  • baescbaesc Member

    Thanks for the message @pdexheimer. I've tried your suggestion... after including "AC INFO" field, I had to also include an "AF INFO" field... and now I'm getting the same error message, except this time it looks like this:
    `##### ERROR ------------------------------------------------------------------------------------------

    ERROR stack trace

    java.lang.IllegalStateException: Key AN found in VariantContext field INFO at Chr1:4000371 but this key isn't defined in the VCFHeader. We require all VCFs to have complete VCF headers by default.
    at org.broadinstitute.variant.variantcontext.writer.VCFWriter.fieldIsMissingFromHeaderError(VCFWriter.java:601)
    at org.broadinstitute.variant.variantcontext.writer.VCFWriter.add(VCFWriter.java:266)
    at org.broadinstitute.sting.gatk.io.storage.VariantContextWriterStorage.add(VariantContextWriterStorage.java:175)
    at org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub.add(VariantContextWriterStub.java:242)
    at org.broadinstitute.sting.gatk.walkers.beagle.BeagleOutputToVCF.map(BeagleOutputToVCF.java:359)
    at org.broadinstitute.sting.gatk.walkers.beagle.BeagleOutputToVCF.map(BeagleOutputToVCF.java:78)
    at org.broadinstitute.sting.gatk.traversals.TraverseLociNano$TraverseLociMap.apply(TraverseLociNano.java:267)
    at org.broadinstitute.sting.gatk.traversals.TraverseLociNano$TraverseLociMap.apply(TraverseLociNano.java:255)
    at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.executeSingleThreaded(NanoScheduler.java:274)
    at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.execute(NanoScheduler.java:245)
    at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:144)
    at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:92)
    at org.broadinstitute.sting.gatk.traversals.TraverseLociNano.traverse(TraverseLociNano.java:48)
    at org.broadinstitute.sting.gatk.executive.LinearMicroScheduler.execute(LinearMicroScheduler.java:99)
    at org.broadinstitute.sting.gatk.GenomeAnalysisEngine.execute(GenomeAnalysisEngine.java:311)
    at org.broadinstitute.sting.gatk.CommandLineExecutable.execute(CommandLineExecutable.java:113)
    at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:245)
    at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:152)
    at org.broadinstitute.sting.gatk.CommandLineGATK.main(CommandLineGATK.java:91)

    ERROR ------------------------------------------------------------------------------------------
    ERROR A GATK RUNTIME ERROR has occurred (version 2.6-3-gdee51c4):
    ERROR
    ERROR Please check the documentation guide to see if this is a known problem
    ERROR If not, please post the error, with stack trace, to the GATK forum
    ERROR Visit our website and forum for extensive documentation and answers to
    ERROR commonly asked questions http://www.broadinstitute.org/gatk
    ERROR
    ERROR MESSAGE: Key AN found in VariantContext field INFO at Chr1:4000371 but this key isn't defined in the VCFHeader. We require all VCFs to have complete VCF headers by default.
    ERROR ------------------------------------------------------------------------------------------

    `
    ...whereby I see no "AN" Key anywhere in the vcf file....

    Thanks in advance for your help!
    Cheers,
    Chris

  • baescbaesc Member

    ... ok thats weird.
    I relented and added an "AN" in the header and now it seems to swallow it....

    Thanks, @pdexheimer!

    Cheers,
    Chris

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hmm, can you post the line you got for Chr1:4000371 ?

  • pdexheimerpdexheimer Member ✭✭✭✭

    That makes sense - AC, AF, and AN are all pretty reasonable things for something called "calculateChromosomeCounts" to output. Glad that worked!

  • baescbaesc Member

    Hi @Geraldine_VdAuwera, here's the line for position 4000371:
    Chr1 4000371 . C T 130 . AC=2;AC1=2;AF=1.00;AF1=1;AN=2;DP=11;DP4=0,0,4,5;FQ=-54;MQ=26;NumGenotypesChanged=0;VDB=0.0380 GT:DP:GQ:OG:PL:SP 1|1:9:60:.:163,27,0:0
    The affected part of the header now reads:
    `##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">

    INFO=<ID=AC1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele count (no HWE assumption)">

    INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">

    INFO=<ID=AF1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele frequency (assuming HWE)">

    INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">

    `
    Thanks for your help!
    Cheers,
    Chris

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Ah, there's your AN... Just checking that the program is not being unreasonable.

    I'll look into why the header is not coming out correctly in the first place.

  • baescbaesc Member

    Hi @Geraldine_VdAwera, that would be great.
    Thanks for your help.
    Cheers,
    Chris

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Of course. We're not seeing this bad behavior locally though, would it be possible for you to upload a test file so we can debug? Instructions are here: http://www.broadinstitute.org/gatk/guide/article?id=1894

  • baescbaesc Member

    Hi @Geraldine_VdAuwera, just a quick question before delving into the details: The .vcf File being read into BeagleOutputToVCF originates from SAMTOOLS and contains an AC1 and AF1 instead of a normal AC and AF. Furthermore, AN is missing completely in the SAMTOOLS vcf. Could this be the reason for the (bit of) confusion?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    I don't think so because it doesn't explain why the GATK would emit a VCF with the correct fields without documenting those fields in the header.

  • baescbaesc Member

    Ok @Geraldine_VdAuwera, I've uploaded a folder (Comment_7947) to the ftp site. Hope I did this right! Cheers,
    Chris

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Thanks Chris! We'll look into this and let you know.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hey @baesc, I finally found some time to look at your bug report but I can't find any folder named "Comment_7947" anywhere in our incoming directory. Are you sure the upload completed successfully? Can you try zipping it first and uploading again?

  • baescbaesc Member

    Hi @Geraldine_VdAuwera, sorry for the delayed response - just had a bit of an extended weekend... I had a look at ftp://ftp.broadinstitute.org/incoming/ and found my folder "Comment_7947" uploaded on 29.07.2013... it contains a .zip file with a little example.... Should I have saved this somewhere else? Cheers, Chris

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    No worries, glad you're enjoying summer!

    I see the problem -- you uploaded the files to the general Broad FTP instead of our team's GATK-specific FTP. I'm not set up to access the general Broad FTP. Could you please redo the upload, this time with the address/credentials stated here: http://www.broadinstitute.org/gatk/guide/article?id=1215

    Thanks and sorry for the inconvenience!

  • baescbaesc Member

    Hi @ Geraldine_VdAuwera, my bad ;-). Ok, I'm now uploading my .zip file to the correct ftp address, only now its called "zipfile.zip". Cheers, Chris

  • baescbaesc Member

    Hi @Geraldine_VdAuwera, could you find the zipfile? Its called "zipfile.zip" and I hope I put in the right directory. Cheers,
    Chris

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @baesc,

    Sorry this comes so late, but yes I did find your file. I just had a look at it though and it looks like you didn't include the input VCF file. But in any case, after re-reading the thread I realized that it's the VCF file produced by Samtools, not the one produced by GATK, that had the bad header, which is then rejected by GATK. Is that correct?

  • baescbaesc Member

    Hi @Geraldine_VdAuwera, thats correct. I first read the samtools .vcf file into GATK, run beagle, and then try to convert back to a .vcf file using GATK again...
    Cheers,
    Chris

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Okay then, I must have misread the issue earlier. If it's Samtools' fault there's nothing we can do here. Sorry for all the bother!

Sign In or Register to comment.