Fix reported bug where no-call genotypes with no reads get genotype (#…

…5667) * Fix reported bug where no-call genotypes with no reads get genotype posterior probabilities and calls
broadinstitute · Feb 14, 2019 · a7ecde4 · a7ecde4
1 parent 6e28a92
commit a7ecde4
Show file tree

Hide file tree

Showing 5 changed files with 36 additions and 2 deletions.
diff --git a/...org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java b/...org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java
@@ -187,7 +187,21 @@ public static double[] parsePosteriorsIntoProbSpace(Genotype genotype) {
 
     // return the double[] of likelihoods if available, otherwise null
     private static double[] getLikelihoodsVector(Genotype genotype) {
-        return genotype.hasLikelihoods() ? genotype.getLikelihoods().getAsVector() : null;
+        return hasRealLikelihoods(genotype) ? genotype.getLikelihoods().getAsVector() : null;
+    }
+
+    /**
+     * return true if the genotype contains real data for likelihoods, not just zero placeholder
+     */
+    private static boolean hasRealLikelihoods(Genotype genotype) {
+        if (!genotype.hasLikelihoods()) {
+            return false;
+        }
+        if (genotype.hasDP() && genotype.getDP() == 0) {
+            //if there's no depth and the PLs are all zero, then there's really no data here; don't rely on no-call
+            return MathUtils.arrayMax(genotype.getPL()) > 0;
+        }
+        return true;
     }
 
     /**

diff --git a/src/main/java/org/broadinstitute/hellbender/utils/MathUtils.java b/src/main/java/org/broadinstitute/hellbender/utils/MathUtils.java
@@ -1221,6 +1221,11 @@ public static int maxElementIndex(final double[] array, final int endIndex) {
         return maxElementIndex(Utils.nonNull(array), 0, endIndex);
     }
 
+    public static int arrayMax(final int[] array) {
+        Utils.nonNull(array);
+        return array[maxElementIndex(array)];
+    }
+
     public static double arrayMax(final double[] array) {
         Utils.nonNull(array);
         return array[maxElementIndex(array)];

diff --git a/...oadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java b/...oadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java
@@ -358,7 +358,9 @@ public void testGVCFModeGenotypePosteriors() throws Exception {
 
         for (final VariantContext vc : results.getRight()) {
             final Genotype g = vc.getGenotype(0);
-            Assert.assertTrue(g.hasExtendedAttribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY));
+            if (g.hasDP() && g.getDP() > 0 && g.hasGQ() && g.getGQ() > 0) {
+                Assert.assertTrue(g.hasExtendedAttribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY));
+            }
             if (isGVCFReferenceBlock(vc) ) {
                 Assert.assertTrue(!vc.hasAttribute(GATKVCFConstants.GENOTYPE_PRIOR_KEY));
             }

diff --git a/...ute/hellbender/tools/walkers/variantutils/CalculateGenotypePosteriorsIntegrationTest.java b/...ute/hellbender/tools/walkers/variantutils/CalculateGenotypePosteriorsIntegrationTest.java
@@ -57,6 +57,7 @@ public void testNumRefWithPanel() throws IOException {
 
     @Test
     //use the first 20 variants to save time; they have a nice range of AC from 4 to over 4000
+    //three variants have GTs with zero depth and PLs=[1,0,0], which should get PPs
     public void testUsingDiscoveredAF() throws IOException {
         final IntegrationTestSpec spec = new IntegrationTestSpec(
                         " -O %s" +

diff --git a/...dinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtilsUnitTest.java b/...dinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtilsUnitTest.java
@@ -391,6 +391,18 @@ public void testDifferentRefAllelesAndNonRefCounts() {
         Assert.assertEquals(arraysEq(expectedPPs, _mleparse((List<Integer>)test1result.getGenotype(0).getAnyAttribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY))), "");
     }
 
+    @Test
+    public void testGenotypesWithNoDataDoNotChange() {
+        final GenotypeBuilder gb = new GenotypeBuilder("s1").
+                alleles(Arrays.asList(Allele.NO_CALL, Allele.NO_CALL)).DP(0).AD(new int[]{0,0}).GQ(0).PL(new int[]{0,0,0});
+        final VariantContext vc = makeVC("noCall", Arrays.asList(Aref, Allele.create("T")), gb.make());
+        final List<VariantContext> resources = new ArrayList<>();
+        resources.add(new VariantContextBuilder(makeVC("2", Arrays.asList(Aref,T))).attribute(GATKVCFConstants.MLE_ALLELE_COUNT_KEY,500).attribute(VCFConstants.ALLELE_NUMBER_KEY,1000).make());
+        final VariantContext result = PosteriorProbabilitiesUtils.calculatePosteriorProbs(vc, resources, 0, defaultOptions);
+        Assert.assertTrue(result.getGenotype("s1").isNoCall());
+        Assert.assertFalse(result.getGenotype("s1").hasExtendedAttribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY));
+    }
+
     private double[] pl2gl(final int[] pl) {
         final double[] gl = new double[pl.length];
         for ( int idx = 0; idx < gl.length; idx++ ) {