Variant in VCF of multiple samples called by HaplotypeCaller absent in their respective BAM files

Arvand88Arvand88 Member
edited July 2017 in Ask the GATK team

Hello GATK team,

I followed GATK best practices and called variants with haplotypecaller in 6 exome samples. However, in 4 patients (total) I have a variant on Chr12 that is absent in the BAM file. the variant is a big deletion (bigger than 10 nucleotides) and ONLY in one read of one of the samples I can see the variant. I'm confused about what has happened here. Can you please explain how can a variant be called while not present in the BAM?

Best Answer

  • shleeshlee Cambridge admin
    Accepted Answer

    @Arvand88,

    No need to rerun your analysis. GVCF mode is only one valid approach. I assume you are running the latest GATKv3.7. Otherwise, please try your analysis again with v3.7.

    However, in 4 patients (total) I have a variant on Chr12 that is absent in the BAM file. the variant is a big deletion (bigger than 10 nucleotides) and ONLY in one read of one of the samples I can see the variant.

    When you say the above, are you looking at your raw BAM or the reassembled HaplotypeCaller BAMOUT? If the latter, then something is odd and we'll need to get a bug report from you. Instructions for bug reports are at https://software.broadinstitute.org/gatk/guide/article?id=1894.

Answers

  • shleeshlee CambridgeMember, Broadie, Moderator admin

    Hi @Arvand88,

    Can you please outline your workflow and provide the exact commands you use to generate GVCFs (if using -ERC GVCF mode) and the final VCF.

    If you haven't run VQSR, then you haven't followed the Best Practices. Calling with sensitivity and then filtering based on annotation profiles using VQSR completes are Best Practice recommendation. If you find that the variant is questionable, then see if VQSR removes it.

    Remember also that GATK's focus is on calling variant sites with sensitivity at the cohort level, whatever the allele and zygosity. This is certainly preferable to missing out on calling true positive sites. It is up to you to determine, via your VQSR filtering, what level of precision you need for your downstream application.

  • Arvand88Arvand88 Member
    edited July 2017

    Dear @shlee
    Thank you for getting back to me.
    I haven't used Haplotypecaller in GVCF mode. I used it to call variants on recalibrated BAM files of samples with the following command:

    GATK \ 
    -T HaplotypeCaller \ 
    -R ucsc.hg19.fasta \ 
    -I 1175recal.bam \ 
    --genotyping_mode DISCOVERY \ 
    -o *.vcf 
    

    Moreover, I applied hard filters because I was under the impression that once the number of samples is less than 30, VQSR cannot be applied accurately. I read this in Best Practices documentations.

    hard filters for SNPs:

    --filterExpression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || SOR > 3.0 || ReadPosRankSum < -8.0"
    

    hard filters for INDELs:

    --filterExpression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || ReadPosRankSum < -20.0"
    

    Afterward, I merged filtered SNPs and filtered INDELs and ran GATK VariantAnnotator:

    java -jar GenomeAnalysisTK.jar \ 
    -R reference.fasta \ 
    -T VariantAnnotator \ 
    -I input.bam \ 
    -o output.vcf \ 
    -A Coverage –A QualByDepth –A TandemRepeatAnnotator –A VariantType \ 
    -V input.vcf \ 
    -L input.vcf \ 
    --dbsnp dbsnp138.vcf
    

    Next, I applied the following custom filters by VCFfilter from VCFlib:

    vcffilter -f "DP > 9 & QUAL > 29 & FILTER = PASS" *gatkanno.vcf > *final.vcf
    

    Finallly, annotated all with ANNOVAR.

    So do you think I should go back and run HaplotypeCaller in gvcf mode? Could you point me to a document because I am not sure what I should do once I have the gVCFs

    Thank you

    Post edited by shlee on
  • shleeshlee CambridgeMember, Broadie, Moderator admin
    Accepted Answer

    @Arvand88,

    No need to rerun your analysis. GVCF mode is only one valid approach. I assume you are running the latest GATKv3.7. Otherwise, please try your analysis again with v3.7.

    However, in 4 patients (total) I have a variant on Chr12 that is absent in the BAM file. the variant is a big deletion (bigger than 10 nucleotides) and ONLY in one read of one of the samples I can see the variant.

    When you say the above, are you looking at your raw BAM or the reassembled HaplotypeCaller BAMOUT? If the latter, then something is odd and we'll need to get a bug report from you. Instructions for bug reports are at https://software.broadinstitute.org/gatk/guide/article?id=1894.

  • @shlee
    Yes I use the latest version (3.7)
    I checked the variants in the recalibrated BAM files (output from BQSR) and the deduped BAM files (output from PICARD) on IGV. All pairs for all samples showed the same result and there was no difference whatsoever.
    I did not include the BAMOUT option but now that you mentioned it and considering the outcome I think I should have. I will run the haplotypecaller again with the BAMOUT option and let you know what I find.
    I will also let you know what happens with the GVCF approach. I am trying everything because this strange variant was my best pathogenic one in this pedigree. Maybe VQSR can help me find something else.

  • Arvand88Arvand88 Member
    edited July 2017

    @shlee
    I checked the bamout files and they seem to be in order with the VCFs. I see the reads showing the deletion listed in VCF.
    As for GVCF approach, I did run GenotypeGVCFs on the 6 samples I have simultaneously, but I think I cannot trust the results due to low quantity of samples. Many good variants have been left out in comparison with HC normal mode.
    Thank you for mentioning bamout argument. It helped clear my doubts about some variants but raised questions about some others. I am opening a new discussion since the subject is irrelevant to the current title.

  • shleeshlee CambridgeMember, Broadie, Moderator admin

    Thanks for the followup @Arvand88. For good depth coverage samples, e.g. ~30x PCR-free, we expect the -ERC GVCF mode + GenotypeGVCFs to be as sensitive as the HC normal mode. I'll be interested to see your followup. Sheila's back to manning the forum so I'm only on call for one or few questions a day. I hope we can figure out what is going on with your data and you will be able to move ahead in your research.

Sign In or Register to comment.