Why does --sample-ploidy/-ploidy option not work in HaplotypeCaller or GenotypeGVCFs?
GATK HaplotypeCaller and GenotypeGVCFs both continue to call heterozygous SNPs even when either haploid flag (--sample-ploidy/-ploidy 1) is selected.
Can you please provide
a) GATK version used =
4.1.4.1
b) Exact GATK commands used
./gatk-4.1.4.1/gatk --java-options "-Xmx8g" HaplotypeCaller \
-R $ref \
-ploidy 1 \
-I $indir/accepted_hits_sortedRG.bam \
--emit-ref-confidence GVCF \
-O $indir/alignedreads.g.vcf \
-L chr1:2100000-3400000 -L chr2:500000-1500000 -L chr2:2500000-4500000
./gatk-4.1.4.1/gatk --java-options "-Xmx35g" GenotypeGVCFs \
-R $ref \
-V $indir/alignedreads.g.vcf \
--read-filter NotSecondaryAlignmentReadFilter \
--output $indir/${pref}_genotyped.vcf \
-ploidy 1
c) The entire error log if applicable.
no error log - haploid haplotyping/genotyping just doesn't work. both programs return "heterozygous" (i.e. false) SNPs. e.g. lines below:
chr2 3417270 . C A 710.01 PASS . GT:AD:DP:GQ:PL 1:30,46:76:99:720,0
chr2 3417271 . G T 710.01 PASS . GT:AD:DP:GQ:PL 1:30,46:76:99:720,0
chr2 3417272 . C A 710.01 PASS . GT:AD:DP:GQ:PL 1:30,46:76:99:720,0
chr2 3417277 . C T 665.01 PASS . GT:AD:DP:GQ:PL 1:30,46:76:99:675,0
chr2 3417281 . C G 710.01 PASS . GT:AD:DP:GQ:PL 1:30,46:76:99:720,0
chr2 3417282 . A C 710.01 PASS . GT:AD:DP:GQ:PL 1:30,46:76:99:720,0
chr2 3417315 . C T 579.01 PASS . GT:AD:DP:GQ:PL 1:31,46:77:99:589,0
chr2 3417387 . C T 663.01 PASS . GT:AD:DP:GQ:PL 1:26,41:67:99:673,0
chr2 3417396 . C A 696.01 PASS . GT:AD:DP:GQ:PL 1:26,41:67:99:706,0
chr2 3417405 . A G 696.01 PASS . GT:AD:DP:GQ:PL 1:26,42:68:99:706,0
chr2 3417418 . G T 395.01 PASS . GT:AD:DP:GQ:PL 1:24,34:58:99:405,0
chr2 3417468 . G A 665.01 PASS . GT:AD:DP:GQ:PL 1:20,36:56:99:675,0
chr2 3417471 . G A 710.01 PASS . GT:AD:DP:GQ:PL 1:19,36:55:99:720,0
chr2 3417474 . C T 755.01 PASS . GT:AD:DP:GQ:PL 1:19,36:55:99:765,0
chr2 3417513 . C T 850.01 PASS . GT:AD:DP:GQ:PL 1:17,35:52:99:860,0
chr2 3417521 . G C 928.01 PASS . GT:AD:DP:GQ:PL 1:16,38:54:99:938,0
chr2 3417542 . G A 980.01 PASS . GT:AD:DP:GQ:PL 1:17,39:56:99:990,0
-
Hi,
This site might not be haploid to begin with hence the AD values. HaplotypeCaller is doing the right thing though and you can tell by the GT filed that it is providing Haploid genotype calls.
If you question is based on results from a pileup based tool showing a different result then we can look into why that might be.
-
I'm working with a haploid organism. Therefore, it can’t possibly have multiple alleles at all of the sites in the pasted data snippet, as is shown in the AD field. Instead, what the data tell us is that the region in question is duplicated in the test genome but not in the reference. BUT IN TERMS OF (HAPLOTYPING/GENOTYPING, THE TEST GENOME IS WILD-TYPE (I.E. NO SNPS) AT EACH OF THE REFERENCE GENOME SITES REPORTED. Therefore, these are all false SNP calls. GATK should know this because it’s been told it’s calling in a HAPLOID organism. This is an ERROR in the program. You can’t say an organism only has the ALT allele (which is what the GT=1 field implies), when you’ve detected two alleles. In fact, based on it having returned these false calls, as well as many other bizarre genotyping decisions, I can tell that GATK is making invalid assumptions, and arbitrary decisions, in the calculations it makes when genotyping in triplicated/quadruplicated/etc. chromosome regions in haploids. This in turn tells me the genotyper doesn’t properly consider haploid genome organization and dynamics and their impacts on the haplotyping/genotyping process.
-
Hi Mark Farman
Sorry for the delay. We are looking into this and will get back to you shortly.
-
HI Mark Farman
I checked with our dev team and this is what they said:
HaplotypeCaller is expecting messy reads that may not always align correctly. In that situation, even in a haploid organism, it's possible to see reads supporting two or even more alleles at a given site: one of the alleles is real, and the others are driven by mapping artifacts or sequencing errors. In this case it's seeing a lot of support for the reference allele and even more support for the alt allele, so it calls alt. It always lists the other possible genotypes it considered, but it isn't saying that the sample is het, it's saying there's some evidence it's ref, and some evidence it's alt, this is how much for each. This isn't fundamentally a problem of ploidy, but that something messy is happening at those sites.
It sounds like you are expecting duplications, which can lead to messy alignments that can be mistaken for SNPs. HaplotypeCaller should be able to handle small duplications in high-complexity regions, but if those conditions don't hold, it can make errors. Not knowing the nature of these duplications, they could also have resulted in a paralogous sequence variant, which would look "het". And of course there could be a bug. Showing images of pile ups (or sharing a link to a bam trimmed down to these regions) would help us figure out what's going on in this case.
Please sign in to leave a comment.
4 comments