HaplotypeCaller does not call all variants

noeDnoeD Member
edited October 2015 in Ask the GATK team

Hello.
I am implementing my pipeline to Roche and I have a list of variants that are validated with Sanger.
In my pipeline I am using HaplotypeCaller to call variants, but there are variants that are not called. I have seen the bam file with IGV and I have found those variants in my reads.
For example:

  • I have SNP at this coordinate of BRCA2: 32906729 (A->C), the reads with this snp are 100% (150 reads in total)
  • I have SNP at this coordinate of BRCA1: 41215825 (G->A), the reads with this snp are 100% (260 reads in total)
    but I have not found they in my vcf.

Why is there this problem? Could you help me? Thank you in advance
Best regards

Tagged:

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
  • noeDnoeD Member

    Thank you very much for your response. I'm sorry but I have not understood if the recalibrated.bam file is the file after indelRealigner.
    Thank you in advance

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
  • noeDnoeD Member

    this is the screenshot of original bam for BRCA2: 32906729 (A->C)

  • noeDnoeD Member

    This is the screenshot of bamout.bam

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Can you please post the VCF output for the region surrounding this SNP? We'll need to see both the GVCF intermediate and the final VCF file, if you're doing joint calling, or just the VCF if you're not. And we'll also need to see the command line you ran.

  • noeDnoeD Member

    I am sorry but I have found errors in my bed file and I have resolved they. But there are a few number of variants that I have not found.
    For example:
    41219780 41219780 A->G, the reads with this snp are 100% (219)
    I have done the screenshot of bam file and bamout on IGV, and in the second are not reads in that position.
    This is the g.vcf ouput for the region surrounding the SNP:
    chr17 41215825 . C T, 8499.77 . DP=240;MLEAC=2,0;MLEAF=1.00,0.00;MQ=60.00 GT:AD:DP:GQ:PL:SB 1/1:0,240,0:240:99:8528,722,0,8528,722,8528:0,0,131,109
    chr17 41215826 . G . . END=41216036 GT:DP:GQ:MIN_DP:PL 0/0:247:99:246:0,120,1800
    chr17 41216037 . C . . END=41216046 GT:DP:GQ:MIN_DP:PL 0/0:247:0:103:0,0,0
    chr17 41219593 . T . . END=41219779 GT:DP:GQ:MIN_DP:PL 0/0:219:99:214:0,120,1800
    chr17 41219780 . T . . END=41219780 GT:DP:GQ:MIN_DP:PL 0/0:219:0:219:0,0,0
    chr17 41219781 . C . . END=41219803 GT:DP:GQ:MIN_DP:PL 0/0:219:99:219:0,120,1800
    chr17 41219804 . T . . END=41219804 GT:DP:GQ:MIN_DP:PL 0/0:219:0:219:0,0,0
    chr17 41219805 . G . . END=41219811 GT:DP:GQ:MIN_DP:PL 0/0:219:99:219:0,120,1800
    chr17 41219812 . T . . END=41219822 GT:DP:GQ:MIN_DP:PL 0/0:219:0:76:0,0,0
    chr17 41222907 . G . . END=41223093 GT:DP:GQ:MIN_DP:PL 0/0:195:99:187:0,120,1800
    chr17 41223094 . T C, 9943.77 . DP=319;MLEAC=2,0;MLEAF=1.00,0.00;MQ=60.00 GT:AD:DP:GQ:PL:SB 1/1:0,319,0:319:99:9972,959,0,9972,959,9972:0,0,176,143
    In the vcf I have not this region, I am reporting the nearest:
    chr13 32968808 . CT C 774.73 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.086;ClippingRankSum=-0.329;DP=156;FS=0.626;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=-0.078;QD=4.97;ReadPosRankSum=-4.378;SOR=0.765 GT:AD:DP:GQ:PL 0/1:93,63:156:99:812,0,1175
    chr17 41215825 . C T 8569.77 . AC=2;AF=1.00;AN=2;DP=242;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=28.81;SOR=0.910 GT:AD:DP:GQ:PL 1/1:0,242:242:99:8598,728,0
    chr17 41223094 . T C 9943.77 .

    Thank you very much for your help

  • noeDnoeD Member

    This is my command line:

    java -Xmx8g -jar GenomeAnalysisTK.jar -T RealignerTargetCreator -R hg19.fa -I 454RG8.bam -o 454RG8.IndelRealigner.intervals --known Mills_and_1000G_gold_standard.indels.hg19.sites.vcf -L coordinate.bed

    java -jar -Xmx8gGenomeAnalysisTK.jar -T IndelRealigner -R hg19.fa -I 454RG8.bam -targetIntervals 454RG8.IndelRealigner.intervals -o bwaAll8.realigned.bam -known Mills_and_1000G_gold_standard.indels.hg19.sites.vcf

    java -jar -Xmx8g GenomeAnalysisTK.jar -T BaseRecalibrator -Rhg19.fa -I bwaAll8.realigned.bam -knownSites dbsnp_138.hg19.vcf -knownSites Mills_and_1000G_gold_standard.indels.hg19.sites.vcf -o bwaAll8.recal.table -L coordinate.bed

    java -jar -Xmx8g GenomeAnalysisTK.jar -T PrintReads -R hg19.fa -I bwaAll8.realigned.bam -BQSR bwaAll8.recal.table -o bwaAll8.recal.bam

    java -jar -Xmx8g GenomeAnalysisTK.jar -T BaseRecalibrator -R hg19.fa -I bwaAll8.realigned.bam -knownSites dbsnp_138.hg19.vcf -knownSites Mills_and_1000G_gold_standard.indels.hg19.sites.vcf -BQSR bwaAll8.recal.table -o bwaAll8.after_recal.table -L coordinate.bed

    java -jar -Xmx8g GenomeAnalysisTK.jar -T HaplotypeCaller -R hg19.fa -I bwaAll8.recal.bam -o bwaAll8.g.vcf -ERC GVCF --doNotRunPhysicalPhasing -L coordinate.bed --variant_index_type LINEAR --variant_index_parameter 128000 -bamout bwaAll8.gvcf.bam

    java -jar -Xmx8g GenomeAnalysisTK.jar -T HaplotypeCaller -R hg19.fa -I bwaAll8.recal.bam -L coordinate.bed --genotyping_mode DISCOVERY -stand_emit_conf 10 -stand_call_conf 30 -o bwaAll8.vcf

  • noeDnoeD Member

    This is the screenshot of IGV

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Try forcing the bamout output for the region as described in point 4 of this document: https://www.broadinstitute.org/gatk/guide/article?id=5484

  • noeDnoeD Member

    I have forced the bamout output and the result is the same: I have not reads in the region.
    This is the command that I have used:
    java -jarGenomeAnalysisTK.jar -T HaplotypeCaller -R hg19.fa -I bwaAll8.recal.bam -o b2.vcf -L chr17:41219780 -bamout bamout8s.bam

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    No, there's a difference between emitting the bamout file (which only contains reads if there is a variant call) and forcing the program to output the reads even if there is no variant call there. Please read point 4 of the document I linked. You need to add some more arguments to your command line.

  • noeDnoeD Member

    thank you very very much for your help.
    Please find in the attached the screenshot of the new output.
    This is the command that I have used:
    java -jar GenomeAnalysisTK.jar -T HaplotypeCaller -R hg19.fa -I bwaAll8.recal.bam -L chr17:41219780 -o hc_forced2.vcf -bamout force_bamout8.bam -forceActive -disableOptimizations

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @noeD
    Hi,

    What does the VCF look like at the site after running the command? Can you also try running the same command with -L chr17:41219480-41220080. Usually we run with some padding around the site of interest.

    It looks like the reason this site is not being called is because the tool is confused between all three possible genotypes (all the PLs are 0). Can you confirm all the mapping qualities and base qualities are good at the site?

    Thanks,
    Sheila

  • noeDnoeD Member

    In the vcf I have not reported variant. The mapping qualities are 60 and base qualities are over 25.
    Thank you for your help

    bb.png 30.6K
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hmm, if you look at the sequence context it's pretty clear that there's a repeated region, so the reads could map in different places. The forced bamout is not informative because although one configuration is chosen that shows the variant, when the program goes back to calculate the genotype likelihoods, it cannot determine confidently which is the correct configuration. That is why it cannot output a variant call at this site, unfortunately. It is a limitation of the method in dealing with repeated sequence contexts. I think using longer reads would help break the tie. What kind of sequencing technology are you using?

  • noeDnoeD Member

    I am working with Roche 454. I have tried to use UnifiedGenotyper and I have found my variants. How is it possible? Anyway I have not INDEL in my vcf.
    The command that I have used is:
    java -jar GenomeAnalysisTK.jar -T UnifiedGenotyper -glm BOTH -R hg19.fa -I bwaAll8.recal.bam --dbsnp dbsnp_138.hg19.vcf -o resultMID8.raw.vcf -L coordinate.bed
    Thank you in advance

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hmm, if I remember correctly 454 data is pretty weak in areas with homopolymers; this could be causing your problem.

    You're not getting any indels with -glm BOTH? That's surprising...

  • noeDnoeD Member

    No, I have not indels, I have tried also -glm INDEL and the vcf is empty

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    That's very odd. My guess is that there's something wrong with the sequence data. Do you have other samples that you can test?

  • noeDnoeD Member

    @Geraldine_VdAuwera said:
    Hmm, if I remember correctly 454 data is pretty weak in areas with homopolymers; this could be causing your problem.

    You're not getting any indels with -glm BOTH? That's surprising...

    Anyway I have the same problem at this position 41219804 41219804 A/G and I have not a homopolymer...

  • noeDnoeD Member

    @Geraldine_VdAuwera said:
    That's very odd. My guess is that there's something wrong with the sequence data. Do you have other samples that you can test?

    I have already tried with other samples and I have the same empty vcf. So I do not understand why if I use HaplotyCaller I find they, but I do not find some variants, if use UnifiedGenotyper I find the variants that I have not found with HaplotyCaller but I don't find others variants. I'm going crazy :(

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @noeD
    Hi,

    Can you confirm you are using the latest version of GATK?

    Thanks,
    Sheila

  • noeDnoeD Member

    If I am not mistaken yes, but just to be on the safe side how can I verify the version? Thank you very much

  • tommycarstensentommycarstensen United KingdomMember ✭✭✭

    @noeD java GATK.jar -version and it's printed in all of your log files and the metadata lines of your output VCF files.

  • tommycarstensentommycarstensen United KingdomMember ✭✭✭

    @noeD Not sure if it will make any difference (actually doubt it), but can you try to set minPruning to 1 instead of the default 2 as per this thread?

Sign In or Register to comment.