Treat missing truth samples as no-call; remove no-call from AF calcul…

…ations
broadinstitute · Feb 1, 2023 · 69ec8ad · 69ec8ad
1 parent 52017bb
commit 69ec8ad
Show file tree

Hide file tree

Showing 7 changed files with 21 additions and 6 deletions.
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/sv/SVAlleleCounter.java b/src/main/java/org/broadinstitute/hellbender/tools/sv/SVAlleleCounter.java
@@ -50,6 +50,7 @@ public int getNumber() {
      */
     private static Map<Allele, Long> computeCounts(final Collection<Genotype> genotypes) {
         return genotypes.stream().map(Genotype::getAlleles).flatMap(Collection::stream)
+                .filter(a -> !(a == null || a.equals(Allele.NO_CALL)))
                 .collect(Collectors.groupingBy(a -> a, Collectors.counting()));
     }
 

diff --git a/src/main/java/org/broadinstitute/hellbender/tools/sv/concordance/SVConcordanceAnnotator.java b/src/main/java/org/broadinstitute/hellbender/tools/sv/concordance/SVConcordanceAnnotator.java
@@ -164,7 +164,9 @@ public Boolean cnvGenotypesMatch(final String sample, final SVCallRecord eval, f
     @VisibleForTesting
     protected GenotypeConcordanceStates.TruthState getTruthState(final Genotype g) {
         // Treat non-existent genotype as hom-ref
-        if (g == null || g.isHomRef()) {
+        if (g == null) {
+            return GenotypeConcordanceStates.TruthState.NO_CALL;
+        } else if (g.isHomRef()) {
             return GenotypeConcordanceStates.TruthState.HOM_REF;
         } else if (g.isHet()) {
             return GenotypeConcordanceStates.TruthState.HET_REF_VAR1;

diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVConcordance.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVConcordance.java
@@ -45,7 +45,7 @@
  * after meeting minimum overlap criteria. Evaluation variants that are sucessfully matched are annotated with
  * genotype concordance metrics, including allele frequency of the truth variant. See output header for descriptions
  * of the specific fields. Note that genotypes of samples that are present in the evaluation VCF but not the truth
- * VCF are assumed to be homozygous-reference. For multi-allelic CNVs, only a copy state concordance metric is
+ * VCF are assumed to be no-calls for concordance. For multi-allelic CNVs, only a copy state concordance metric is
  * annotated.
  *
  * <h3>Inputs</h3>

diff --git a/src/test/java/org/broadinstitute/hellbender/tools/sv/SVAlleleCounterTest.java b/src/test/java/org/broadinstitute/hellbender/tools/sv/SVAlleleCounterTest.java
@@ -99,6 +99,17 @@ public Object[][] testCounterData() {
                         new double[]{1.},
                         2
                 },
+                // No-call not counted toward AN
+                {
+                        new Allele[]{Allele.SV_SIMPLE_INS},
+                        new Allele[][]{
+                                {Allele.SV_SIMPLE_INS},
+                                {Allele.NO_CALL}
+                        },
+                        new int[]{1},
+                        new double[]{1.},
+                        1
+                },
                 // Multi-allelic
                 {
                         new Allele[]{Allele.SV_SIMPLE_DEL, Allele.SV_SIMPLE_DUP},

diff --git a/src/test/java/org/broadinstitute/hellbender/tools/sv/concordance/ClosestSVFinderTest.java b/src/test/java/org/broadinstitute/hellbender/tools/sv/concordance/ClosestSVFinderTest.java
@@ -92,7 +92,8 @@ public void testMetrics() {
         final SVCallRecord outEval1 = out1.get(0);
         Assert.assertEquals(outEval1.getId(), eval1.getId());
         assertConcordanceMembers(outEval1, truth1.getId());
-        Assert.assertEquals(outEval1.getAttributes().get(GATKSVVCFConstants.GENOTYPE_CONCORDANCE_INFO), 1.);
+        // Would be 1 if missing samples were treated as hom-ref
+        Assert.assertEquals(outEval1.getAttributes().get(GATKSVVCFConstants.GENOTYPE_CONCORDANCE_INFO), 0.5);
 
         final SVCallRecord eval2 = SVTestUtils.makeRecord(
                 "eval2",

diff --git a/...t/java/org/broadinstitute/hellbender/tools/sv/concordance/SVConcordanceAnnotatorTest.java b/...t/java/org/broadinstitute/hellbender/tools/sv/concordance/SVConcordanceAnnotatorTest.java
@@ -283,7 +283,7 @@ public void testGetStateFromGenotype(final Allele[] genotypeAlleles,
     public void testGetTruthNullGenotypeState() {
         final SVConcordanceAnnotator collapser = new SVConcordanceAnnotator(false);
         final GenotypeConcordanceStates.TruthState actualTruthState = collapser.getTruthState(null);
-        Assert.assertEquals(actualTruthState, GenotypeConcordanceStates.TruthState.HOM_REF);
+        Assert.assertEquals(actualTruthState, GenotypeConcordanceStates.TruthState.NO_CALL);
     }
 
     @Test(expectedExceptions = IllegalArgumentException.class)

diff --git a/src/test/java/org/broadinstitute/hellbender/tools/walkers/sv/SVAnnotateEngineUnitTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/sv/SVAnnotateEngineUnitTest.java
@@ -427,8 +427,8 @@ public Object[][] getSVTypesAndSegmentsTestData() {
                         null},
                 { createVariantContext("chr22", 10510000, 10694100, "chr22", null, "N",
                         "<DEL:ME:LINE1>", 184100, null, null, null),
-                        SVAnnotateEngine.StructuralVariantAnnotationType.DEL,
-                        createListOfSVSegments(SVAnnotateEngine.StructuralVariantAnnotationType.DEL,
+                        GATKSVVCFConstants.StructuralVariantAnnotationType.DEL,
+                        createListOfSVSegments(GATKSVVCFConstants.StructuralVariantAnnotationType.DEL,
                                 new SimpleInterval[]{ new SimpleInterval("chr22", 10510000, 10694100)}),
                         null},
                 { createVariantContext("chr22", 10524000, 10710000, "chr22", null, "N",