From 16417bbf34e918d51780be490fe43b52a7fb5fce Mon Sep 17 00:00:00 2001
From: Valentin Ruano Rubio <valentin@broadinstitute.org>
Date: Tue, 11 Oct 2016 11:54:17 -0400
Subject: [PATCH] Fixes NaN issue in new Qual calculator

Fixes issue #1491
---
 .../genotyper/afcalc/AlleleFrequencyCalculator.java      | 9 ++++++++-
 .../variantutils/GenotypeGVCFsIntegrationTest.java       | 9 +++++++++
 2 files changed, 17 insertions(+), 1 deletion(-)

diff --git a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/afcalc/AlleleFrequencyCalculator.java b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/afcalc/AlleleFrequencyCalculator.java
index febdebb9c..abfddab82 100644
--- a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/afcalc/AlleleFrequencyCalculator.java
+++ b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/afcalc/AlleleFrequencyCalculator.java
@@ -159,8 +159,15 @@ public final class AlleleFrequencyCalculator extends AFCalculator {
                         pOfNonZeroAltAlleles[alleleIndex] += genotypePosterior);
             }
 
+            // Make sure that we handle appropriately pOfNonZeroAltAlleles that are close to 1; values just over 1.0 due to
+            // rounding error would result in NaN.
+            //  As every allele is present in at least one genotype, the p-non-zero-count for
+            // any allele is bound above by 1.0 - minimum genotype posterior because at least one genotype
+            // does not contain this allele.
+            final double maximumPNonZeroCount = 1.0 - MathUtils.arrayMin(genotypePosteriors);
+
             for (int allele = 0; allele < numAlleles; allele++) {
-                log10POfZeroCountsByAllele[allele] += Math.log10(1 - pOfNonZeroAltAlleles[allele]);
+                log10POfZeroCountsByAllele[allele] += Math.log10(1.0 - Math.min(maximumPNonZeroCount, pOfNonZeroAltAlleles[allele]));
             }
         }
 
diff --git a/protected/gatk-tools-protected/src/test/java/org/broadinstitute/gatk/tools/walkers/variantutils/GenotypeGVCFsIntegrationTest.java b/protected/gatk-tools-protected/src/test/java/org/broadinstitute/gatk/tools/walkers/variantutils/GenotypeGVCFsIntegrationTest.java
index 98b806f40..47e813204 100644
--- a/protected/gatk-tools-protected/src/test/java/org/broadinstitute/gatk/tools/walkers/variantutils/GenotypeGVCFsIntegrationTest.java
+++ b/protected/gatk-tools-protected/src/test/java/org/broadinstitute/gatk/tools/walkers/variantutils/GenotypeGVCFsIntegrationTest.java
@@ -681,4 +681,13 @@ public class GenotypeGVCFsIntegrationTest extends WalkerTest {
         spec.disableShadowBCF();
         executeTest("testGenotypingSpanningDeletionAcrossLines", spec);
     }
+
+    @Test
+    public void testNewQualNaNBugFix() {
+        final WalkerTestSpec spec = new WalkerTestSpec(
+            baseTestString(" -newQual -V " + privateTestDir + "input-newqual-nan-bug-fix.vcf", b37KGReferenceWithDecoy),
+                Collections.singletonList("503f4193c22fbcc451bd1c425b8b6bf8"));
+        spec.disableShadowBCF();
+        executeTest("testNewQualNaNBugFix", spec);
+    }
 }
\ No newline at end of file