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!

Why does HaplotypeCaller call heterozygotes when all reads are identical?

Hello, I've been exploring the vcf output from HaplotypeCaller-in-gVCF-mode, and I noticed that for some SNPs (37/~250K), HaplotypeCaller has called a heterozygote when according to the VCF, all the reads support either the reference or the alternative allele (90% of the time, it is when all reads support the alternative allele).

There are for fairly low-coverage sites (all < 10), although from exploring the data, it looks like "GQ" values are quite high (20+) and that they are all loci that have been phased by HaplotypeCaller.

Is HaplotypeCaller doing something smart here with the phasing to "recover" information that isn't in the raw sequences? I can't really figure out what that might be, but I'm also don't really know what is going on under the hood with physical phasing.

In case it matters, I am working with mosquito exome data, 150bp PE illumina reads, GATK 4.1.2.0, java 1.8.0_181. I've been following the Broad best practices pretty closely, although I haven't managed to figure out how to bootstrap base recalibration yet.

An example line from the VCF, with the individual of concern bolded (note that there are some bulks in here, and I think they might be exhibiting the same behavior, but so far I'm only looking a diploid single individuals):

NW_021837111.1 53637 . T G 107.95 . AC=1;AF=9.434e-03;AN=106;DP=221;ExcessHet=3.0103;FS=0.000;InbreedingCoeff=-0.1312;MLEAC=1;MLEAF=9.434e-03;MQ=60.00;QD=32.06;SOR=2.833GT:AD:DP:GQ:PGT:PID:PL:PS 0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0:55,0:55:8:.:.:0,8,17,25,35,45,55,67,79,91,105,120,137,155,176,199,226,257,296,346,417,537,1800 0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0:51,0:51:11:.:.:0,11,23,36,50,65,82,100,120,144,170,202,241,291,361,482,1800 0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0:36,0:36:7:.:.:0,7,14,22,30,38,47,57,67,78,89,102,116,132,149,169,192,219,252,294,354,456,1530 0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0:45,0:45:7:.:.:0,7,15,23,32,41,50,60,70,82,94,107,120,136,152,170,191,214,241,272,311,361,432,552,1800 0|1:0,3:3:25:0|1:53637_T_G:123,0,25:53637 ./.:0,0:0:.:.:.:0,0,0 ./.:0,0:0:.:.:.:0,0,0 ./.:0,0:0:.:.:.:0,0,0 ./.:0,0:0:.:.:.:0,0,0 ./.:0,0:0:.:.:.:0,0,0 ./.:0,0:0:.:.:.:0,0,0 ./.:0,0:0:.:.:.:0,0,0 0/0:1,0:1:3:.:.:0,3,22 0/0:1,0:1:3:.:.:0,3,36 0/0:7,0:7:18:.:.:0,18,270 0/0:2,0:2:6:.:.:0,6,75 0/0:2,0:2:6:.:.:0,6,77 0/0:1,0:1:3:.:.:0,3,40 0/0:5,0:5:15:.:.:0,15,205 0/0:3,0:3:9:.:.:0,9,112 ./.:0,0:0:.:.:.:0,0,0 0/0:5,0:5:15:.:.:0,15,213 ./.:0,0:0:.:.:.:0,0,0 0/0:4,0:4:12:.:.:0,12,124

Answers

  • Tiffany_at_BroadTiffany_at_Broad Cambridge, MAMember, Administrator, Broadie, Moderator admin

    Hi @jhb Thanks for the question. I will get back to you soon!

  • Tiffany_at_BroadTiffany_at_Broad Cambridge, MAMember, Administrator, Broadie, Moderator admin

    Hi @jhb would you be able to post a screenshot of your bamout at the site you are questioning?

  • jhbjhb Member

    Thanks, Tiffany, I'll try to get one. I've discovered that I've been having some problems upstream of this step as well, so maybe fixing those will take care of this too? Wishful thinking, I know....

  • SkyWarriorSkyWarrior TurkeyMember ✭✭✭

    I noticed that the site you show here is clearly not a biallelic site or you are working with pooled samples an ploidy is not set properly. Former is more likely but it seems there is more in that region that it cannot be explained by a biallelic snp thats why the allelic representation is off. Can you elaborate on your experiment design?

  • jhbjhb Member

    It's a mix of diploid individuals and bulks. The first four "individuals" listed are pools of some number of individuals (~10, but different for each), and the remainder are diploid individuals.

  • Tiffany_at_BroadTiffany_at_Broad Cambridge, MAMember, Administrator, Broadie, Moderator admin

    Hi @jhb when you say "bulks" do you mean the same thing as pools? Ok, let us know if you get the same results after you make some fixes upstream. A bamout will be helpful either way.

  • jhbjhb Member

    OK, the upstream stuff is now (apparently) in order.

    This phenomenon still persists at about 5% of sites, but based on the bamouts, I think that I have figured things out, but I would appreciate your comments.

    I looked at a few of these loci, and it seems to me that what is happening is that HaplotypeCaller does its local realignment of the original bowtie/whatever alignment. This realigns what was a homozygous site into a heterozygous one. E.g.:

    or

    Then, when you get a vcf file from GenotypeGVCFs a couple steps downstream, it has forgotten about HaplotypeCaller's realignment, and reports stats about how many reads support various genotypes, etc, based on the original bowtie alignment.

    This would explain the discrepancy. Followup questions:

    1) Am I right that GenotypeGVCFs provides stats based on the original alignment and not HaplotypeCaller's local alignment?
    2) If #1 is yes, is there a way to change this to get stats based on HaplotypeCaller's alignment?
    3) If #1 is yes, does this have any consequences for SNP-filtering? I.e., is it meaningful to do something like:

    VariantFiltration --genotypeFilterExpression "DP < 10"
    

    If VariantFiltration is filtering based on the information encoded in the GenotypeGVCFs-produced VCFs, then I would have my doubts about this, but I'm also not really sure about what else you could do instead.

Sign In or Register to comment.