From 3fc1d94619e06935e518c05fee198678ca7da482 Mon Sep 17 00:00:00 2001 From: Laura Gauthier Date: Tue, 12 Feb 2019 15:18:49 -0500 Subject: [PATCH 1/3] Fix reported bug where no-call genotypes with no reads get genotype posterior probabilities and calls --- .../PosteriorProbabilitiesUtils.java | 21 ++++++++++++++++++- ...lateGenotypePosteriorsIntegrationTest.java | 1 + .../PosteriorProbabilitiesUtilsUnitTest.java | 12 +++++++++++ 3 files changed, 33 insertions(+), 1 deletion(-) diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java index ac2a287581c..3d1ca5a7703 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java @@ -187,7 +187,26 @@ 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) { + int[] pls = genotype.getPL(); + for (int i = 0; i < pls.length; i++) { + if (pls[i] > 0) { + return true; + } + } + return false; //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 true; } /** diff --git a/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/CalculateGenotypePosteriorsIntegrationTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/CalculateGenotypePosteriorsIntegrationTest.java index 25420825ecd..7fb6508318e 100644 --- a/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/CalculateGenotypePosteriorsIntegrationTest.java +++ b/src/test/java/org/broadinstitute/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/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtilsUnitTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtilsUnitTest.java index 5c43604b682..7b473e620d2 100644 --- a/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtilsUnitTest.java +++ b/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtilsUnitTest.java @@ -391,6 +391,18 @@ public void testDifferentRefAllelesAndNonRefCounts() { Assert.assertEquals(arraysEq(expectedPPs, _mleparse((List)test1result.getGenotype(0).getAnyAttribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY))), ""); } + @Test + public void testGenotypesWithNoDataDoNotChange() { + final GenotypeBuilder gb = new GenotypeBuilder("s1"); + gb.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(Allele.create("A",true), Allele.create("T")), gb.make()); + final List 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++ ) { From 23f40de8c6d3809bcfbd5835eb932402aea956d1 Mon Sep 17 00:00:00 2001 From: Laura Gauthier Date: Wed, 13 Feb 2019 10:38:03 -0500 Subject: [PATCH 2/3] Missed a HC posteriors test for a file where there was missing data --- .../haplotypecaller/HaplotypeCallerIntegrationTest.java | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/test/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java index 4f226c41913..c3ec39862ad 100644 --- a/src/test/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java +++ b/src/test/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java @@ -357,7 +357,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)); } From a51b63b932da2155ab2bc64e6e85c5f4d381c6f4 Mon Sep 17 00:00:00 2001 From: Laura Gauthier Date: Thu, 14 Feb 2019 10:33:26 -0500 Subject: [PATCH 3/3] Address review comments --- .../variantutils/PosteriorProbabilitiesUtils.java | 9 ++------- .../org/broadinstitute/hellbender/utils/MathUtils.java | 5 +++++ .../PosteriorProbabilitiesUtilsUnitTest.java | 6 +++--- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java index 3d1ca5a7703..4f50b907498 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java @@ -198,13 +198,8 @@ private static boolean hasRealLikelihoods(Genotype genotype) { return false; } if (genotype.hasDP() && genotype.getDP() == 0) { - int[] pls = genotype.getPL(); - for (int i = 0; i < pls.length; i++) { - if (pls[i] > 0) { - return true; - } - } - return false; //if there's no depth and the PLs are all zero, then there's really no data here; don't rely on no-call + //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 index 4bd8230ddc8..d1906b9e563 100644 --- 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/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtilsUnitTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtilsUnitTest.java index 7b473e620d2..894ef9e64ff 100644 --- a/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtilsUnitTest.java +++ b/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtilsUnitTest.java @@ -393,9 +393,9 @@ public void testDifferentRefAllelesAndNonRefCounts() { @Test public void testGenotypesWithNoDataDoNotChange() { - final GenotypeBuilder gb = new GenotypeBuilder("s1"); - gb.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(Allele.create("A",true), Allele.create("T")), gb.make()); + 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 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);