Interpretation of GQ, PL and GT for sites with low-confidence calls

[andrewparkermorgan]

DeepVariant developers and experienced users,

I am a recent adopter of DeepVariant -- all my prior experience is with "traditional" non-CNN-based callers and I am having trouble developing an intuition for what sample-level features should make me suspicious of a genotype call. For context, I am working with WGS data for a reasonable large cohort (~1600) of wild house mice, sequenced on several different platforms (mostly Illumina NextSeq or HiSeq) at a range of depths (~2x - 50x). We are interested in hard genotype calls, of course, but for sites/samples where a high-confidence call is not available, we hope to aggregate information from the PLs to estimate the site frequency spectrum.

There are a couple of example cases that I find confusing.

• GT = 0/0 but no reads (AD = 0/0). GQ is low and PL = 0,0,0 as expected, but why is this not reported as missing? The PL communicates the lack of information but if we are hoping to filter on missingness, it seems like we are going to overestimate the call rate, possibly by a lot.
• GT = ./. with decent number of reads (ie AD = 7,0) and PL = 8,0,158 suggesting there is enough information to make a call, but GQ = 0. Why is this marked as missing? How should we interpret the PL?

Below is a snippet from the VCF at a site which includes both of the above cases, as well as some nonmissing calls with reasonable support.

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  LEJ7226 LEJ7269 LEJ7346 LEJ7402 LEJ7436 LEJ7444 GR2     GR4     GR6     MVZ-Mamm-192074
19      4000046 19_4000046_G_A  G       A       50      PASS    AF=0.00973618;AQ=50;F_MISSING=0.00871731;DP=27219;AD_REF=26865;AD_ALT=346;MIN_GQ_ALT=26;NS=1592;AN=12;MAF=0.00973618;AC=5;AC_Het=27;AC_Hom=4;AC_Hemi=0;HWE=0.0086604;ExcHet=0.99968      GT:DP:AD:GQ:PL:RNC      ./.:5:5,0:0:14,0,104:II 0/0:0:0,0:1:0,0,0:..    ./.:5:5,0:0:14,0,104:II 0/0:0:0,0:1:0,0,0:..     ./.:6:6,0:0:11,0,131:II ./.:7:7,0:0:8,0,158:II  0/1:4:2,2:10:10,0,18:.. 0/1:4:2,2:15:16,0,20:.. 0/1:9:5,4:22:24,0,24:.. 1/1:11:0,11:32:50,44,0:..

Thanks,
Andrew Morgan

[akolesnikov]

Hi Andrew,

The PL values are calculated directly from the posterior probabilities output by the CNN. However, these values (and the corresponding GT) can be updated during the postprocess_variants stage. For example, when multi-allelic sites are reconciled or certain candidate alleles are filtered out.

For your second example: the genotype failed the quality threshold. The site was initially classified as a RefCall, but since --cnn_homref_call_min_gq defaults to 20, any RefCall with GQ < 20 is converted to ./. while retaining the original PL array.

To give you a definitive answer on the first example, could you share the raw, unmerged DeepVariant output before downstream processing? Alternatively, if you can provide a small BAM snippet reproducing the behavior, we'd be happy to inspect it locally.

[andrewparkermorgan]

Thanks for the explanation, but I'm still confused about whether the PL can be trusted for these sites which are affected by the postprocessing step. For example, in the snippet in my original post, sample LEJ7226 gets a missing call and has 7 reads for the reference and 0 for the alternate allele, but the PL array (14,0,104) implies that the most likely genotype would have been heterozygous. My group is planning some analyses that use just the PLs in order to account for uncertainty of discrete genotype calls in lower-coverage samples, so it is really important for us to understand in what cases the we should ignore or zero-out the PL.

I have an embedded IGV viewer with the relevant site and supporting reads here.

Looking at the reads, it appears that the overturned reference calls are cases when most reads carry the reference base, one read carries a non-reference base (with low base quality) but that base is different than what ultimately shows up as the ALT allele in the merged VCF. I went back to the unmerged gVCFs and can see that the no-call originates before merging. So I guess I understand what is happening here -- DeepVariant does not believe the site is truly a variant but also sees evidence for a non-reference base and therefore reports a missing call. We could get rid of these by adjusting the cnn_homref_call_min_gq threshold.

Incidentally, it looks like postprocess_variants updates the genotype call but does not alter the PL array (I think the relevant point in the code is here. GLnexus then seems to just carry those forward, unchanged. That creates a situation where the PL is describing likelihoods for diploid genotypes form a different set of alleles than is ultimately reported in the merged VCF. I find that confusing, and I haven't seen clear documentation of that behavior elsewhere (unless I have just missed it.)

My question about why DeepVariant would report a homozygous reference call for a site with zero read support remains, though.

[kishwarshafin]

@andrewparkermorgan

Is it possible to look at the record before GLNexus to see what exactly DeepVariant omits in this case? I agree this is a bit confusing.

When you say "has 7 reads for the reference and 0 for the alternate allele", are you referring to IGV? Because internally DeepVariant does a de-bruijn graph-based realignment which can alter the alignment. If you have the record from DeepVariant's output the AF would be helpful to look at.