diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtils.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtils.java index af901b2dff2..e633cac19ed 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtils.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtils.java @@ -46,7 +46,7 @@ private AlleleSubsettingUtils() {} // prevent instantiation * @param allelesToKeep the subset of alleles to use with the new Genotypes * @param assignmentMethod assignment strategy for the (subsetted) PLs * @param depth the original variant DP or 0 if there was no DP - * @param emitEmptyPLs force the output of a PL array even if there is no data + * @param suppressUninformativePLs force the output of a PL array even if there is no data * @return a new non-null GenotypesContext */ public static GenotypesContext subsetAlleles(final GenotypesContext originalGs, final int defaultPloidy, @@ -55,7 +55,7 @@ public static GenotypesContext subsetAlleles(final GenotypesContext originalGs, final GenotypePriorCalculator gpc, final GenotypeAssignmentMethod assignmentMethod, final int depth, - final boolean emitEmptyPLs) { + final boolean suppressUninformativePLs) { Utils.nonNull(originalGs, "original GenotypesContext must not be null"); Utils.nonNull(allelesToKeep, "allelesToKeep is null"); Utils.nonEmpty(allelesToKeep, "must keep at least one allele"); @@ -90,7 +90,7 @@ public static GenotypesContext subsetAlleles(final GenotypesContext originalGs, newLog10GQ = -0.1*g.getGQ(); } - final boolean useNewLikelihoods = emitEmptyPLs || + final boolean useNewLikelihoods = !suppressUninformativePLs || (newLikelihoods != null && (depth != 0 || GATKVariantContextUtils.isInformative(newLikelihoods))); final GenotypeBuilder gb = new GenotypeBuilder(g); diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/GenotypingEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/GenotypingEngine.java index af0ee663543..0af78f89bb1 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/GenotypingEngine.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/GenotypingEngine.java @@ -133,7 +133,7 @@ public VariantContext calculateGenotypes(final VariantContext vc, final Genotype if (maxAltAlleles < vc.getAlternateAlleles().size()) { final List allelesToKeep = AlleleSubsettingUtils.calculateMostLikelyAlleles(vc, defaultPloidy, maxAltAlleles); final GenotypesContext reducedGenotypes = allelesToKeep.size() == 1 ? GATKVariantContextUtils.subsetToRefOnly(vc, defaultPloidy) : - AlleleSubsettingUtils.subsetAlleles(vc.getGenotypes(), defaultPloidy, vc.getAlleles(), allelesToKeep, gpc, GenotypeAssignmentMethod.SET_TO_NO_CALL, vc.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), true); + AlleleSubsettingUtils.subsetAlleles(vc.getGenotypes(), defaultPloidy, vc.getAlleles(), allelesToKeep, gpc, GenotypeAssignmentMethod.SET_TO_NO_CALL, vc.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), false); reducedVC = new VariantContextBuilder(vc).alleles(allelesToKeep).genotypes(reducedGenotypes).make(); } @@ -181,7 +181,7 @@ && noAllelesOrFirstAlleleIsNotNonRef(outputAlternativeAlleles.alleles) && givenA // create the genotypes //TODO: omit subsetting if output alleles is not a proper subset of vc.getAlleles final GenotypesContext genotypes = outputAlleles.size() == 1 ? GATKVariantContextUtils.subsetToRefOnly(vc, defaultPloidy) : - AlleleSubsettingUtils.subsetAlleles(vc.getGenotypes(), defaultPloidy, vc.getAlleles(), outputAlleles, gpc, configuration.genotypeArgs.genotypeAssignmentMethod, vc.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), true); + AlleleSubsettingUtils.subsetAlleles(vc.getGenotypes(), defaultPloidy, vc.getAlleles(), outputAlleles, gpc, configuration.genotypeArgs.genotypeAssignmentMethod, vc.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), false); if (configuration.genotypeArgs.usePosteriorProbabilitiesToCalculateQual && hasPosteriors(genotypes)) { final double log10NoVariantPosterior = phredNoVariantPosteriorProbability(outputAlleles, genotypes) * -.1; diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCF.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCF.java index 77a2aff4a50..3321e19caee 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCF.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCF.java @@ -507,9 +507,14 @@ protected GenotypeBuilder changeCallToHomRefVersusNonRef(final VariantContext lo } else { final List bestAlleles = AlleleSubsettingUtils.calculateMostLikelyAlleles(lowQualVariant, genotype.getPloidy(), 1); final Allele bestAlt = bestAlleles.stream().filter(a -> !a.isReference()).findFirst().orElse(Allele.NON_REF_ALLELE); //allow span dels + //we care about the best alt even though it's getting removed because NON_REF should get the best likelihoods + //it shouldn't matter that we're passing in different alt alleles since the GenotypesContext only knows + // the called alleles and this is a reference genotype that will stay hom-ref final GenotypesContext context = AlleleSubsettingUtils.subsetAlleles(lowQualVariant.getGenotypes(), genotype.getPloidy(), lowQualVariant.getAlleles(), Arrays.asList(inputRefAllele, bestAlt), - null, GenotypeAssignmentMethod.BEST_MATCH_TO_ORIGINAL, lowQualVariant.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), false); //BEST_MATCH to avoid no-calling low qual genotypes + null, GenotypeAssignmentMethod.BEST_MATCH_TO_ORIGINAL, //BEST_MATCH to avoid no-calling low qual genotypes + lowQualVariant.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), + false); //emitEmptyPLs = true to make sure we always subset final Genotype subsetG = context.get(0); gb = new GenotypeBuilder(subsetG).noAttributes(); //remove attributes because hom ref blocks shouldn't have posteriors //subsetting may strip GQ and PLs for low qual genotypes @@ -564,7 +569,7 @@ VariantContext cleanUpHighQualityVariant(final VariantContext variant) { newAlleleSetUntrimmed.removeAll(allelesToDrop); final GenotypesContext gc = AlleleSubsettingUtils.subsetAlleles(variant.getGenotypes(), genotype.getPloidy(), variant.getAlleles(), newAlleleSetUntrimmed, null, GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN, - variant.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), false); + variant.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), true); if (gc.get(0).isHomRef() || !gc.get(0).hasGQ() || gc.get(0).getAlleles().contains(Allele.NO_CALL)) { //could be low quality or no-call after subsetting if (dropLowQuals) { return null; @@ -775,7 +780,7 @@ private static void addQualAnnotations(final Map destination, fi //TODO: this isn't going to work for DRAGEN's genotype posteriors final GenotypesContext gc = AlleleSubsettingUtils.subsetAlleles(updatedAllelesVC.getGenotypes(), updatedAllelesGenotype.getPloidy(), updatedAllelesVC.getAlleles(), Arrays.asList(updatedAllelesVC.getReference(), alt), null, - GenotypeAssignmentMethod.BEST_MATCH_TO_ORIGINAL, 0, false); + GenotypeAssignmentMethod.BEST_MATCH_TO_ORIGINAL, 0, true); //assignment method doens't really matter as long as we don't zero out PLs; don't need depth to get PLs for quals final Genotype subsettedGenotype = gc.get(0); diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariants.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariants.java index a4468a791ef..42367fe6619 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariants.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariants.java @@ -5,7 +5,6 @@ import htsjdk.variant.variantcontext.Genotype; import htsjdk.variant.variantcontext.GenotypeBuilder; import htsjdk.variant.variantcontext.GenotypesContext; -import htsjdk.variant.variantcontext.GenotypeLikelihoods; import htsjdk.variant.variantcontext.VariantContext; import htsjdk.variant.variantcontext.VariantContextBuilder; import htsjdk.variant.variantcontext.VariantContextUtils; @@ -34,9 +33,7 @@ import org.broadinstitute.hellbender.tools.walkers.genotyper.AlleleSubsettingUtils; import org.broadinstitute.hellbender.tools.walkers.genotyper.GenotypeAssignmentMethod; import org.broadinstitute.hellbender.utils.samples.MendelianViolation; -import org.broadinstitute.hellbender.utils.samples.PedigreeValidationType; import org.broadinstitute.hellbender.utils.samples.SampleDB; -import org.broadinstitute.hellbender.utils.samples.SampleDBBuilder; import org.broadinstitute.hellbender.utils.Utils; import org.broadinstitute.hellbender.utils.variant.GATKVariantContextUtils; import org.broadinstitute.hellbender.utils.variant.*; @@ -1033,7 +1030,7 @@ private VariantContext subsetRecord(final VariantContext vc, final boolean prese // fix the PL and AD values if sub has fewer alleles than original vc final GenotypesContext subGenotypesWithOldAlleles = sub.getGenotypes(); //we need sub for the right samples, but PLs still go with old alleles newGC = sub.getNAlleles() == vc.getNAlleles() ? subGenotypesWithOldAlleles : - AlleleSubsettingUtils.subsetAlleles(subGenotypesWithOldAlleles, 0, vc.getAlleles(), sub.getAlleles(), null, GenotypeAssignmentMethod.DO_NOT_ASSIGN_GENOTYPES, vc.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), false); + AlleleSubsettingUtils.subsetAlleles(subGenotypesWithOldAlleles, 0, vc.getAlleles(), sub.getAlleles(), null, GenotypeAssignmentMethod.DO_NOT_ASSIGN_GENOTYPES, vc.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), true); } else { newGC = sub.getGenotypes(); } diff --git a/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVariantContextUtils.java b/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVariantContextUtils.java index 04c49743810..5f4e706917a 100644 --- a/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVariantContextUtils.java +++ b/src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVariantContextUtils.java @@ -3,7 +3,6 @@ import com.google.common.annotations.VisibleForTesting; import htsjdk.samtools.SAMSequenceDictionary; import htsjdk.samtools.util.Locatable; -import htsjdk.tribble.TribbleException; import htsjdk.utils.ValidationUtils; import htsjdk.variant.variantcontext.*; import htsjdk.variant.variantcontext.writer.Options; @@ -1688,7 +1687,7 @@ public static List splitVariantContextToBiallelics(final Variant genotypeAssignmentMethodUsed != GenotypeAssignmentMethod.SET_TO_NO_CALL) AlleleSubsettingUtils.addInfoFieldAnnotations(vc, builder, keepOriginalChrCounts); - builder.genotypes(AlleleSubsettingUtils.subsetAlleles(vc.getGenotypes(),2,vc.getAlleles(), alleles, null, genotypeAssignmentMethodUsed,vc.getAttributeAsInt("DP",0), false)); + builder.genotypes(AlleleSubsettingUtils.subsetAlleles(vc.getGenotypes(),2,vc.getAlleles(), alleles, null, genotypeAssignmentMethodUsed,vc.getAttributeAsInt("DP",0), true)); final VariantContext trimmed = trimAlleles(builder.make(), trimLeft, true); biallelics.add(trimmed); } diff --git a/src/test/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtilsUnitTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtilsUnitTest.java index 2b53fff7d04..e2ef9d0b4e4 100644 --- a/src/test/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtilsUnitTest.java +++ b/src/test/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtilsUnitTest.java @@ -143,7 +143,7 @@ public void testUpdatePLsAndADData(final VariantContext originalVC, AlleleSubsettingUtils.subsetAlleles(oldGs, 0, originalVC.getAlleles(), selectedVCwithGTs.getAlleles(), null, GenotypeAssignmentMethod.DO_NOT_ASSIGN_GENOTYPES, - originalVC.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), false); + originalVC.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), true); Assert.assertEquals(actual.size(), expectedGenotypes.size()); for ( final Genotype expected : expectedGenotypes ) { @@ -327,7 +327,7 @@ public void testUninformativePLsAreKeptWhenDepthIsNotZero(){ final List alleles = Arrays.asList(Aref); final Genotype uniformativePL = new GenotypeBuilder("sample", alleles).PL(new int[] {0}).make(); final GenotypesContext result = AlleleSubsettingUtils.subsetAlleles(GenotypesContext.create(uniformativePL), 2, - alleles, alleles, null, GenotypeAssignmentMethod.DO_NOT_ASSIGN_GENOTYPES, 10, false); + alleles, alleles, null, GenotypeAssignmentMethod.DO_NOT_ASSIGN_GENOTYPES, 10, true); final Genotype genotype = result.get(0); Assert.assertTrue(genotype.hasPL()); Assert.assertEquals(genotype.getPL(), new int[]{0}); @@ -412,7 +412,7 @@ public void testAlleleReorderingBehavior() { final GenotypesContext newGs = AlleleSubsettingUtils.subsetAlleles(GenotypesContext.create(g5), 2, threeAlleles, threeAllelesSorted, null, - GenotypeAssignmentMethod.DO_NOT_ASSIGN_GENOTYPES, 10, false); + GenotypeAssignmentMethod.DO_NOT_ASSIGN_GENOTYPES, 10, true); Assert.assertEquals(newGs.get(0).getPL(), new int[] {50, 20, 0, 40, 10, 30}); } diff --git a/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCFUnitTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCFUnitTest.java index c083a620891..c6044669f3f 100644 --- a/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCFUnitTest.java +++ b/src/test/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCFUnitTest.java @@ -110,6 +110,12 @@ public void testLowQualVariantToGQ0HomRef() { Assert.assertTrue(!modified.filtersWereApplied()); Assert.assertEquals(modified.getLog10PError(), VariantContext.NO_LOG10_PERROR); + final Genotype longPls = VariantContextTestUtils.makeG("sample1", Allele.NO_CALL, Allele.NO_CALL, 0,0,0,0,0,0); + final VariantContext lotsOfZeroPls = makeNoDepthVC("lowQualVar", Arrays.asList(LONG_REF, DELETION, Allele.NON_REF_ALLELE), LONG_REF.length(), longPls); + final VariantContext properlySubset = reblocker.lowQualVariantToGQ0HomRef(lotsOfZeroPls).make(); + Assert.assertEquals(properlySubset.getGenotype(0).getPL().length, 3); + Assert.assertEquals(properlySubset.getGenotype(0).getPL(), new int[]{0,0,0}); + //No-calls were throwing NPEs. Now they're not. final Genotype g2 = VariantContextTestUtils.makeG("sample1", Allele.NO_CALL,Allele.NO_CALL); final VariantContext noData = makeDeletionVC("noData", Arrays.asList(LONG_REF, DELETION, Allele.NON_REF_ALLELE), LONG_REF.length(), g2); @@ -425,6 +431,13 @@ private VariantContext makeDeletionVC(final String source, final List al .genotypes(Arrays.asList(genotypes)).unfiltered().log10PError(-3.0).attribute(VCFConstants.DEPTH_KEY, EXAMPLE_DP).make(); } + private VariantContext makeNoDepthVC(final String source, final List alleles, final int refLength, final Genotype... genotypes) { + final int start = DEFAULT_START; + final int stop = start+refLength-1; + return new VariantContextBuilder(source, "1", start, stop, alleles) + .genotypes(Arrays.asList(genotypes)).unfiltered().log10PError(-3.0).make(); + } + private Genotype addAD(final Genotype g, final int... ads) { return new GenotypeBuilder(g).AD(ads).make(); } diff --git a/src/testUtils/java/org/broadinstitute/hellbender/testutils/VariantContextTestUtils.java b/src/testUtils/java/org/broadinstitute/hellbender/testutils/VariantContextTestUtils.java index a651c21419e..febf245b2ba 100644 --- a/src/testUtils/java/org/broadinstitute/hellbender/testutils/VariantContextTestUtils.java +++ b/src/testUtils/java/org/broadinstitute/hellbender/testutils/VariantContextTestUtils.java @@ -186,7 +186,7 @@ public static VariantContext sortAlleles(final VariantContext vc, final VCFHeade result.alleles(sortedAlleles); GenotypesContext newGT = AlleleSubsettingUtils.subsetAlleles(vc.getGenotypes(),2,vc.getAlleles(),sortedAlleles, null, - GenotypeAssignmentMethod.SET_TO_NO_CALL, vc.getAttributeAsInt(VCFConstants.DEPTH_KEY,0), false); + GenotypeAssignmentMethod.SET_TO_NO_CALL, vc.getAttributeAsInt(VCFConstants.DEPTH_KEY,0), true); // Asserting that the new genotypes were calculated properly in case AlleleSubsettingUtils behavior changes if (newGT.getSampleNames().size() != vc.getGenotypes().size()) throw new IllegalStateException("Sorting this variant context resulted in a different number of genotype alleles, check that AlleleSubsettingUtils still supports reordering:" + vc.toString());