UG calling a SNP as heterozygous when it (most likely) isn't

Jlowry

The title kind of explains the situation, but basically I've got a SNP that shows up in IGV that I would call homozygous that the Unified Genotyper has labeled as heterozygous. The total read depth is 35, 32 of which were called as a SNP (A-->T), 2 were called the reference base (A), and one read contained a G. I went through your article describing why a SNP visible in IGV might not get called, and none of those five questions explained this situation. I didn't alter the --hets option at all either. Any help you might be able to offer would be greatly appreciated.

Geraldine_VdAuwera

OK, that makes sense. You need to update to upgrade to a more recent version.

What's happening is that you probably have some low-grade contamination in your sequence data (those non-ref reads). The UG can't just dismiss them, so it calls the site as heterozygous, but if you look at the genotype quality (GQ) I bet you'll find that it's really low -- which means that the UG knows it's a crappy call.

In version 2.2 we added a filter that downsamples per allele, which has the net effect of filtering out low-grade contamination (see the release highlights for 2.2 for more details). This feature is turned on by default, so if you run again with the latest version, those contaminating reads should get downsampled and the UG will be able to make the hom var call that you want to see.

In future, keep in mind that you should look at the GQ before getting too focused on the actual genotype.

Good luck!

Geraldine_VdAuwera

Hi there, can you tell me what version you are using and what is your command line?

Jlowry

I'm using v2.1-11-g13c0244, and here's my command line usage (I'm using it via Galaxy wrapper, not sure if that makes a difference):

java -XX:DefaultMaxRAMFraction=1 -XX:+UseParallelGC -jar /galaxy/galaxy-dist/tool-data/shared/jars/gatk/GenomeAnalysisTK.jar -T UnifiedGenotyper --num_threads 4 --out /galaxy/galaxy-dist/database/files/001/dataset_1046.dat --metrics_file /galaxy/galaxy-dist/database/files/001/dataset_1047.dat -et STANDARD --genotype_likelihoods_model BOTH --standard_min_confidence_threshold_for_calling 20.0 --standard_min_confidence_threshold_for_emitting 20.0 --pedigreeValidationType STRICT --read_filter MappingQuality --min_mapping_quality_score 10 --interval_set_rule UNION --downsampling_type NONE --baq OFF --baqGapOpenPenalty 40.0 --defaultBaseQualities -1 --validation_strictness STRICT --interval_merging ALL --p_nonref_model EXACT --heterozygosity 0.001 --pcr_error_rate 0.0001 --genotyping_mode DISCOVERY --output_mode EMIT_VARIANTS_ONLY --min_base_quality_score 15 --max_deletion_fraction 0.05 --max_alternate_alleles 5 --min_indel_count_for_genotyping 5 --indel_heterozygosity 0.000125 --indelGapContinuationPenalty 10 --indelGapOpenPenalty 45 --indelHaplotypeSize 80 -I /tmp/tmp-gatk-e2zCxM/gatk_input_0.bam -R /tmp/tmp-gatk-e2zCxM/gatk_input.fasta

Jlowry

Thanks for the help, I've made the upgrade and am waiting to see if this resolves the issue. Regardless, it seems like I'd want to check the GQ for possible het var calls, is there a score below which you would say, "Hmm, that seems questionable"?

Geraldine_VdAuwera

It's all relative to the quality of your dataset. One way to empirically derive a threshold is to plot the distribution of GQs for your callset and look in closer detail at a few subsets of calls along the distribution.

Jlowry

That makes sense. Thanks for the suggestion, I'll give it a shot.