Revert AF calculation bug and yet another reblocking fix #7670

ldgauthier · 2022-02-11T21:15:18Z

The warp pipeline tests caught some cases that we apparently didn't have in our integration tests, but now we do!

ldgauthier · 2022-02-11T21:48:27Z

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<Allele> 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), false);
+                    AlleleSubsettingUtils.subsetAlleles(vc.getGenotypes(), defaultPloidy, vc.getAlleles(), allelesToKeep, gpc, GenotypeAssignmentMethod.SET_TO_NO_CALL, vc.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), true);


These go with reverting the boolean variable to force PLs to be output

ldgauthier · 2022-02-11T21:49:05Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtils.java

@@ -94,7 +94,7 @@ public static GenotypesContext subsetAlleles(final GenotypesContext originalGs,
                newLog10GQ = -0.1*g.getGQ();
            }

-            final boolean useNewLikelihoods = !suppressUninformativePLs ||
+            final boolean useNewLikelihoods = forceEmitUninformativePLs ||


The old conditional wasn't right -- I trusted the IntelliJ refactoring too much. This is also a more logical name.

ldgauthier · 2022-02-11T21:49:53Z

.../org/broadinstitute/hellbender/tools/walkers/genotyper/afcalc/AlleleFrequencyCalculator.java

@@ -77,7 +77,7 @@ public static AlleleFrequencyCalculator makeCalculator(final DragstrParams drags
        final double[] log10Likelihoods;
        if (g.hasLikelihoods()) {
            log10Likelihoods = g.getLikelihoods().getAsVector();
-        } else if (g.isHomRef()) {
+        } else if ( g.isHomRef() || g.isNoCall()) {


Taking out the no-call shouldn't have happened in 4.2.5.0, but it did because of poor test coverage. Added new multi-allelic tests to GGVCFs.

droazen

@ldgauthier Back to you with my comments

droazen · 2022-02-14T17:17:30Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtils.java

@@ -55,7 +55,7 @@ public static GenotypesContext subsetAlleles(final GenotypesContext originalGs,
                                                 final GenotypePriorCalculator gpc,
                                                 final GenotypeAssignmentMethod assignmentMethod,
                                                 final int depth,
-                                                 final boolean suppressUninformativePLs) {
+                                                 final boolean forceEmitUninformativePLs) {


Did you check all usages of this method to make sure that the booleans have been flipped when necessary, since this arg now has the opposite meaning?

I took this out instead. The reason it was still there was because I didn't want to change the SelectVariants integration test expected outputs, but that's not a good reason. I want to prefer explicitness and consistency over space savings or whatever the motivation was for dropping PLs sometimes.

droazen · 2022-02-14T17:23:35Z

.../org/broadinstitute/hellbender/tools/walkers/genotyper/afcalc/AlleleFrequencyCalculator.java

@@ -77,7 +77,7 @@ public static AlleleFrequencyCalculator makeCalculator(final DragstrParams drags
        final double[] log10Likelihoods;
        if (g.hasLikelihoods()) {
            log10Likelihoods = g.getLikelihoods().getAsVector();
-        } else if (g.isHomRef()) {
+        } else if ( g.isHomRef() || g.isNoCall()) {


Add a brief comment here documenting the reason why hom refs and no calls are handled in the same case here, to discourage our future selves.

This is just too dangerous. One of the integration tests was doing the wrong thing, so who knows what QUAL came out. Instead I still prohibit no-calls here but call genotypes in the internal subsetting so we should always have calls when we have the right information.

droazen · 2022-02-14T17:29:11Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/SelectVariants.java

@@ -1033,7 +1033,7 @@ private VariantContext subsetRecord(final VariantContext vc, final boolean prese
            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);
+                            vc.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), true);


If the force emit arg is critical here too, add a comment to that effect

droazen · 2022-02-14T17:31:52Z

src/main/java/org/broadinstitute/hellbender/utils/GenotypeUtils.java

@@ -143,7 +143,7 @@ public static GenotypeCounts computeDiploidGenotypeCounts(final VariantContext v
     * @return  true if we have enough info for AF calculation
     */
    public static boolean genotypeIsUsableForAFCalculation(Genotype g) {
-        return g.hasLikelihoods() || (g.isHomRef() && g.hasGQ() && 2 == g.getPloidy());
+        return g.hasLikelihoods() || g.hasGQ() || g.getAlleles().stream().anyMatch(a -> a.isCalled() && a.isNonReference() && !a.isSymbolic());


So this condition got reverted to its previous state too -- can you add a brief comment here explaining the logic for posterity?

droazen · 2022-02-14T17:32:40Z

src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVariantContextUtils.java

@@ -1687,7 +1687,7 @@ public static void splitASFilters(VariantContext vc, List<Map<String, Object>> a
                        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), true));
+                builder.genotypes(AlleleSubsettingUtils.subsetAlleles(vc.getGenotypes(),2,vc.getAlleles(), alleles, null, genotypeAssignmentMethodUsed,vc.getAttributeAsInt("DP",0), false));


force emit is false here -- is that correct?

droazen · 2022-02-15T18:38:12Z

src/test/java/org/broadinstitute/hellbender/tools/walkers/GenotypeGVCFsIntegrationTest.java

+                //23 highly multi-allelic sites across 54 1000G exomes to test allele subsetting and QUAL calculation
+                {getTestFile("multiallelicQualRegression.vcf "),
+                        getTestFile("multiallelicQualRegression.expected.vcf"),
+                        NO_EXTRA_ARGS, hg38Reference}


You've confirmed that this new test case fails with 4.2.5.0?

Yes -- it fails like the Hindenburg.

droazen · 2022-02-15T18:42:10Z

...ava/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCFIntegrationTest.java

+                        " -do-qual-approx --floor-blocks -GQB 20 -GQB 30 -GQB 40 " +
+                        " --" + StandardArgumentDefinitions.ADD_OUTPUT_VCF_COMMANDLINE + " false",
+                Arrays.asList(getToolTestDataDir() + "prodWesOutput.g.vcf"));
+        spec3.executeTest("testWithExactComparison", this);


Perhaps this test method should use a DataProvider to feed it the three command lines, instead of embedding the three cases within one method like this?

ldgauthier · 2022-02-16T17:55:55Z

@droazen there might be some new test failures (which I will address if/when they arise), but this is ready for re-review.

droazen

@ldgauthier Back to you -- just needs a few additional comments and then I think we should be good, provided the WARP and GATK tests pass

droazen · 2022-02-16T19:51:40Z

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<Allele> 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), false);
+                    AlleleSubsettingUtils.subsetAlleles(vc.getGenotypes(), defaultPloidy, vc.getAlleles(), allelesToKeep, gpc, GenotypeAssignmentMethod.BEST_MATCH_TO_ORIGINAL);


Add a comment documenting why we want to use GenotypeAssignmentMethod.BEST_MATCH_TO_ORIGINAL here

droazen · 2022-02-16T20:15:14Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/GenotypingEngine.java

@@ -391,11 +391,12 @@ boolean isVcCoveredByDeletion(final VariantContext vc) {
     * {@link GenotypeLikelihoods#MAX_DIPLOID_ALT_ALLELES_THAT_CAN_BE_GENOTYPED}.
     */
    protected final boolean cannotBeGenotyped(final VariantContext vc) {
-        // protect against too many alternate alleles that we can't even run AF on:
+        //


Fill in or delete empty comment

droazen · 2022-02-16T20:21:45Z

.../org/broadinstitute/hellbender/tools/walkers/genotyper/afcalc/AlleleFrequencyCalculator.java

-        } else if (g.isHomRef()) {
+        } else if ( g.isHomRef()) {
+            //no-call with no PLs seems risky, but there are a few places in the QUAL/AF code where we subset alleles,
+            // but then leave the genotypes as no-calls


This comment seems misplaced/confusing, since no-calls with no PLs will now get handled by the else block below, not this block, and will always produce the Genotype does not contain likelihoods necessary to calculate posteriors exception, right?

else { throw new IllegalStateException("Genotype " + g + " does not contain likelihoods necessary to calculate posteriors."); }

Recommend adding a comment to the else block below explaining clearly why we want to always throw an IllegalStateException here for no calls with no PLs.

droazen · 2022-02-16T20:23:58Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCF.java

@@ -568,8 +566,7 @@ VariantContext cleanUpHighQualityVariant(final VariantContext variant) {
        if(allelesNeedSubsetting && !keepAllAlts) {
            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), true);
+                    newAlleleSetUntrimmed, null, GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN);


Add a comment documenting why we're using GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN here instead of the GenotypeAssignmentMethod.BEST_MATCH_TO_ORIGINAL we're using everywhere else.

droazen · 2022-02-16T20:25:35Z

src/main/java/org/broadinstitute/hellbender/utils/GenotypeUtils.java

@@ -139,6 +139,7 @@ public static GenotypeCounts computeDiploidGenotypeCounts(final VariantContext v
    /**
     * Do we have (or can we infer) likelihoods necessary for allele frequency calculation?
     * Some reblocked and/or DRAGEN GVCFs omit likelihoods for ref blocks, but we can estimate them
+     * If GenomicsDB max alt threshold is too low, variants may also be missimg PLs -- reject those


missimg -> missing

…re's no non-ref Reenable EM in multi-allelic QUAL calculation

because that can cause problems

droazen

@ldgauthier Back to you with a couple of final comments

droazen · 2022-03-08T16:03:39Z

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<Allele> 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), false);
+                    AlleleSubsettingUtils.subsetAlleles(vc.getGenotypes(), defaultPloidy, vc.getAlleles(), allelesToKeep, gpc, GenotypeAssignmentMethod.BEST_MATCH_TO_ORIGINAL);


Confirming that this should be hardcoded to BEST_MATCH_TO_ORIGINAL rather than respecting configuration.genotypeArgs.genotypeAssignmentMethod?

Yes, with no PLs no-calls are just going to cause problems, so keep any 0/0 genotypes as such. I added a comment to that extent.

src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/GenotypingEngine.java

droazen · 2022-03-08T16:07:21Z

.../org/broadinstitute/hellbender/tools/walkers/genotyper/afcalc/AlleleFrequencyCalculator.java

-        } else if (g.isHomRef()) {
+        } else if ( g.isHomRef()) {
+            //no-call with no PLs seems risky, but there are a few places in the QUAL/AF code where we subset alleles,
+            // but then leave the genotypes as no-calls


This comment still seems confusing/misplaced. It talks about handling no-calls in the case that now handles hom-ref (no-calls now fall through to the else block below).

droazen · 2022-03-08T16:11:24Z

src/main/java/org/broadinstitute/hellbender/utils/GenotypeUtils.java

@@ -139,6 +139,7 @@ public static GenotypeCounts computeDiploidGenotypeCounts(final VariantContext v
    /**
     * Do we have (or can we infer) likelihoods necessary for allele frequency calculation?
     * Some reblocked and/or DRAGEN GVCFs omit likelihoods for ref blocks, but we can estimate them
+     * If GenomicsDB max alt threshold is too low, variants may also be missimg PLs -- reject those


missimg -> missing

Is this comment still relevant?

Fixed typo. Yes, still relevant. I clarified that here variants means non-reference genotypes, where we can't estimate PLs.

droazen · 2022-03-08T16:22:29Z

src/test/java/org/broadinstitute/hellbender/tools/walkers/GenotypeGVCFsIntegrationTest.java

+    @Test
+    public void testGDBMaxAltsGreaterThanGGVCFsMaxAlts() throws IOException {
+
+    }


Fill in empty test method

I realized this case was redundant with testMaxAltsToCombineInGenomicsDB. Deleted and added comment to testMaxAlts

droazen

👍 Merge when tests pass @ldgauthier

ldgauthier requested a review from droazen February 11, 2022 21:16

ldgauthier force-pushed the ldg_yetAnotherReblockingFix branch 2 times, most recently from 6f61408 to 483ec41 Compare February 11, 2022 21:47

ldgauthier commented Feb 11, 2022

View reviewed changes

droazen requested changes Feb 15, 2022

View reviewed changes

droazen self-assigned this Feb 15, 2022

droazen requested changes Feb 16, 2022

View reviewed changes

ldgauthier added 3 commits February 17, 2022 10:23

Fix automatic zeroing of array end/non-ref data because sometimes the…

5562afe

…re's no non-ref Reenable EM in multi-allelic QUAL calculation

Previous version had bad AN values

7339c05

Be more clear about max alts arg; don't let no-calls get inferred

a245499

because that can cause problems

ldgauthier force-pushed the ldg_yetAnotherReblockingFix branch from 08150df to a245499 Compare February 17, 2022 15:24

This was referenced Feb 22, 2022

NPE When using GenomicsDB and GenotypeGVCFs with --max-alternates #7687

Closed

Add companion index to MultiVariantWalker input variant arg. #7689

Merged

GenotypeGVCFs error in 4.2.4.1 java.lang.IllegalStateException #7639

Closed

ldgauthier added 5 commits March 2, 2022 13:18

Hook up new max alt GDB arg to Gnarly

baab476

Allow Gnarly's max alt count to be changed via GenotypeCalculateArgs

04cd52c

Oops, broke GGVCFs

d107744

Add regression test for GGVCFs max alt arg mixup

b8456e5

Every. Single. Time.

b51cb54

droazen requested changes Mar 8, 2022

View reviewed changes

ldgauthier added 2 commits March 8, 2022 11:53

Address review comments; just cleanup

dc9ccd5

Actually delete empty test

8f71a0d

droazen approved these changes Mar 8, 2022

View reviewed changes

ldgauthier merged commit 2d2cca2 into master Mar 8, 2022

ldgauthier deleted the ldg_yetAnotherReblockingFix branch March 8, 2022 19:07

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Revert AF calculation bug and yet another reblocking fix #7670

Revert AF calculation bug and yet another reblocking fix #7670

ldgauthier commented Feb 11, 2022

ldgauthier Feb 11, 2022

ldgauthier Feb 11, 2022

ldgauthier Feb 11, 2022

droazen left a comment

droazen Feb 14, 2022

ldgauthier Feb 16, 2022

droazen Feb 14, 2022

ldgauthier Feb 16, 2022

droazen Feb 14, 2022

droazen Feb 14, 2022

droazen Feb 14, 2022

droazen Feb 15, 2022

ldgauthier Feb 16, 2022

droazen Feb 15, 2022

ldgauthier Feb 16, 2022

ldgauthier commented Feb 16, 2022

droazen left a comment

droazen Feb 16, 2022

droazen Feb 16, 2022

droazen Feb 16, 2022

droazen Feb 16, 2022

droazen Feb 16, 2022

droazen left a comment

droazen Mar 8, 2022

ldgauthier Mar 8, 2022

droazen Mar 8, 2022

ldgauthier Mar 8, 2022

droazen Mar 8, 2022

ldgauthier Mar 8, 2022

droazen Mar 8, 2022

ldgauthier Mar 8, 2022

droazen left a comment

Revert AF calculation bug and yet another reblocking fix #7670

Revert AF calculation bug and yet another reblocking fix #7670

Conversation

ldgauthier commented Feb 11, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

droazen left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ldgauthier commented Feb 16, 2022

droazen left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

droazen left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

droazen left a comment

Choose a reason for hiding this comment