fix no data hom refs #8715

ldgauthier · 2024-03-05T19:15:57Z

No description provided.

gatk-bot · 2024-03-05T19:49:44Z

Github actions tests reported job failures from actions build 8161782218
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
integration	17.0.6+10	8161782218.11	logs
integration	17.0.6+10	8161782218.0	logs

gatk-bot · 2024-03-05T23:16:23Z

Github actions tests reported job failures from actions build 8164035876
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
integration	17.0.6+10	8164035876.11	logs
integration	17.0.6+10	8164035876.0	logs

ldgauthier · 2024-03-06T14:51:45Z

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

@@ -306,7 +306,7 @@ public static void makeGenotypeCall(final int ploidy,
                                        final GenotypeAssignmentMethod assignmentMethod,
                                        final double[] genotypeLikelihoods,
                                        final List<Allele> allelesToUse,
-                                        final List<Allele> originalGT,
+                                        final Genotype originalGT,


I had to change this method signature so that we could use the GQ info to distinguish hom-ref genotypes that we want to no-call from hom-refs that stay hom-refs because with reblocking there are no PLs to use.

add to method doc

ldgauthier · 2024-03-06T14:52:46Z

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

@@ -588,7 +588,7 @@ public void testCombineReblockedGVCFs() {
        Assert.assertFalse(g0.hasPL());
        Assert.assertFalse(g0.hasGQ());
        Assert.assertFalse(g0.hasExtendedAttribute(VCFConstants.DEPTH_KEY));
-        Assert.assertTrue(g1.isHomRef());
+        Assert.assertTrue(g1.isNoCall());


In a lot of these tests I was explicitly testing the behavior of the GQ0 hom-ref, but that's exactly what we're trying to change here

ldgauthier · 2024-03-06T14:53:15Z

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

        Assert.assertTrue(vc1.hasAttribute(GATKVCFConstants.INBREEDING_COEFFICIENT_KEY));
-        Assert.assertEquals(vc1.getAttributeAsInt(VCFConstants.ALLELE_NUMBER_KEY, 0), 362);
+        Assert.assertEquals(vc1.getAttributeAsInt(VCFConstants.ALLELE_NUMBER_KEY, 0), 300);


This changes because of the new no-calls

ldgauthier · 2024-03-06T14:53:46Z

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

        final VariantContext vc1 = outputVCs.get(1);
-        Assert.assertTrue(vc1.getAttributeAsDouble(GATKVCFConstants.EXCESS_HET_KEY, 1000.0) < 10.0); //will be ~72 if homRefs aren't counted
+        Assert.assertTrue(vc1.getAttributeAsDouble(GATKVCFConstants.EXCESS_HET_KEY, 0.0) > 50.0); //will be ~72 if homRefs aren't counted


Some hom-refs in previous test result go to no call, so InbreedingCoeff and ExcessHet also change

ldgauthier · 2024-03-06T14:55:29Z

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

@@ -976,4 +982,39 @@ public void dbSNPError() {
        //make sure it doesn't throw an error
        runCommandLine(args);
    }
+
+    @Test
+    public void testNoReadsOutputAsNoCall() {


The input data for this test comes from a GVCF produced from an empty region of one of the NA12878 test bams and a GVCF that was edited to have a variant so we can force that position to be output. (I should probably put this test description in the code.)

jamesemery

Some very minor documentation/comment requests. Without getting very deep into the weeds of those integration test annotation changes (which look fine as per your comments) i think this looks good on the surface.

There are two more abstract questions I will pose:
Do we want to preserve the ability for anybody ever to output the intermidate outputs? We could create a new mode in GenotypeAssignmentMethod for that and change this default one to use the new behavior to satisfy twitter...

The second one is just to ask if these are the only two places (GnarlyGenotyper and CombineGVCFs) that need to change? Its only at the combine step we are throwing these out so the genomicsDB path shouldn't be a problem no? Do we want HC outputting ./. genotypes?

jamesemery · 2024-03-06T15:01:02Z

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

@@ -306,7 +306,7 @@ public static void makeGenotypeCall(final int ploidy,
                                        final GenotypeAssignmentMethod assignmentMethod,
                                        final double[] genotypeLikelihoods,
                                        final List<Allele> allelesToUse,
-                                        final List<Allele> originalGT,
+                                        final Genotype originalGT,


add to method doc

jamesemery · 2024-03-06T15:04:40Z

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

                if (assignmentMethod == GenotypeAssignmentMethod.PREFER_PLS) {
                    if (originalGT == null) {
                        throw new IllegalArgumentException("original GT cannot be null if assignmentMethod is PREFER_PLS");
+                    } else if (originalGT.hasGQ() && originalGT.getGQ() > 0){


i see. this makes sense and was certainly better than my hack attempt to fix this... so its now "prefer_PLs EXCEPT for nocalls" Can you add a note about this in GenotypeAssignmentMethod?

Is there a reason to preserve the intermediate behavior at all? We could split PREFER_PLS to also have a mode that treats the 0/0s like they used to be.

jamesemery · 2024-03-06T15:06:21Z

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

@@ -158,6 +158,7 @@ public Object[][] gvcfsToGenotype() {

                //23 highly multi-allelic sites across 54 1000G exomes to test allele subsetting and QUAL calculation
                //plus one 10-allele WGS variant that's all hom-ref with one GT that has unnormalized PLs from some sort of GenomicsDB corner case
+                //this VCF still has the haploid-looking GDB no-calls as in sample NA21137 at position chr1:3836468 -- allegedly GATK 4.2.5.0 from February 7, 2022, possibly due to --call-genotypes


jamesemery · 2024-03-06T15:20:45Z

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

@@ -804,6 +805,9 @@ public void testWithReblockedGVCF() {
        final List<VariantContext> actualVC = VariantContextTestUtils.getVariantContexts(output);
        Assert.assertFalse(actualVC.stream().anyMatch(vc -> vc.getGenotype(1).isHomRef() && vc.getGenotype(1).hasPL()));  //second sample has a bunch of 0/0s -- shouldn't have PLs

+        //this comes from a callset of NYGC 1000G samples plus syndip
+        //it seems likely that there's a variant that wasn't discovered in the graph because a bunch of samples are hom-ref with PL=[0,0,X]


a minor comment... but its possible to manufacture a sample with lots of assembly failures if you need to test this. You can set the kmer size to 1 and disallow stepped increases to get a sample that can assemble nothing.... Doesn't look necessary for you here though.

Ooh, that is good to know.

jamesemery · 2024-03-06T15:24:26Z

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

@@ -976,4 +982,39 @@ public void dbSNPError() {
        //make sure it doesn't throw an error
        runCommandLine(args);
    }
+
+    @Test
+    public void testNoReadsOutputAsNoCall() {


.../java/org/broadinstitute/hellbender/tools/walkers/gnarlyGenotyper/GnarlyGenotyperEngine.java

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

droazen · 2024-03-07T20:13:13Z

M2 WDL tests failed due to a transient 429 "too many requests" error -- re-running them.

ldgauthier added 2 commits March 5, 2024 11:53

Possible fix for no-call issue

e48051c

Add test with one sample variant, one sample empty

b8a7059

ldgauthier added 3 commits March 5, 2024 15:56

Fix InbreedingCoeff issue, but AN is wrong

7accbf9

Passing GQ fix seems to work; still need Gnarly test

771d38f

New Gnarly test shows good no-call behavior on GVCF with no read data

48647a7

Missed a CombineGVCFs test

eb35245

ldgauthier commented Mar 6, 2024

View reviewed changes

Add test for reblocked GVCF representing no reads

6c691ee

jamesemery approved these changes Mar 6, 2024

View reviewed changes

jamesemery reviewed Mar 6, 2024

View reviewed changes

src/test/java/org/broadinstitute/hellbender/tools/walkers/GenotypeGVCFsIntegrationTest.java Outdated Show resolved Hide resolved

ldgauthier added a commit that referenced this pull request Mar 6, 2024

Possibly more elegant fix than #8715

cdb8a9d

Address review comments

fec031b

ldgauthier merged commit b0463e4 into master Mar 7, 2024
20 checks passed

ldgauthier deleted the ldg_fix_no_data_hom_refs branch March 7, 2024 21:07

nalinigans mentioned this pull request Mar 23, 2024

Allow for GT to be a nocall if GQ and PL[0] are zero instead of homref GenomicsDB/GenomicsDB#332

Merged

nalinigans mentioned this pull request Apr 2, 2024

Allow for GT to be a nocall if GQ and PL[0] are zero instead of homref in GenomicsDB #8759

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

fix no data hom refs #8715

fix no data hom refs #8715

ldgauthier commented Mar 5, 2024

gatk-bot commented Mar 5, 2024 •

edited

Loading

gatk-bot commented Mar 5, 2024 •

edited

Loading

ldgauthier Mar 6, 2024

jamesemery Mar 6, 2024

ldgauthier Mar 6, 2024

ldgauthier Mar 6, 2024

ldgauthier Mar 6, 2024

ldgauthier Mar 6, 2024

jamesemery Mar 6, 2024

jamesemery left a comment

jamesemery Mar 6, 2024

jamesemery Mar 6, 2024

jamesemery Mar 6, 2024

jamesemery Mar 6, 2024

jamesemery Mar 6, 2024

ldgauthier Mar 6, 2024

jamesemery Mar 6, 2024

droazen commented Mar 7, 2024

fix no data hom refs #8715

fix no data hom refs #8715

Conversation

ldgauthier commented Mar 5, 2024

gatk-bot commented Mar 5, 2024 • edited Loading

gatk-bot commented Mar 5, 2024 • edited Loading

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

jamesemery 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

droazen commented Mar 7, 2024

gatk-bot commented Mar 5, 2024 •

edited

Loading

gatk-bot commented Mar 5, 2024 •

edited

Loading