From 3c5dc675ab2ff0e7b73a7bed7802bab3c2bfc37d Mon Sep 17 00:00:00 2001
From: ebanks <ebanks@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Wed, 6 Oct 2010 15:16:41 +0000
Subject: [PATCH] For Guillermo: only decide that something is a clear
 reference call if it is at least 10 times as likely as the next best genotype

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4441 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../genotyper/AlleleFrequencyCalculationModel.java        | 8 +++++---
 .../gatk/walkers/genotyper/UnifiedGenotyperEngine.java    | 2 +-
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/AlleleFrequencyCalculationModel.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/AlleleFrequencyCalculationModel.java
index d63dc5947..33b7e0940 100755
--- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/AlleleFrequencyCalculationModel.java
+++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/AlleleFrequencyCalculationModel.java
@@ -52,6 +52,8 @@ public abstract class AlleleFrequencyCalculationModel implements Cloneable {
         GRID_SEARCH
     }
 
+    private static final double LOG10_REFERENCE_CALL_EPSILON = 1.0;
+
     protected int N;
     protected AlleleFrequencyMatrix AFMatrix;
     protected Set<BiallelicGenotypeLikelihoods> refCalls;
@@ -179,8 +181,7 @@ public abstract class AlleleFrequencyCalculationModel implements Cloneable {
         AFMatrix.clear();
 
         for ( BiallelicGenotypeLikelihoods GL : GLs.values() ) {
-            // todo - gdebug workaround for crash
-            if (false) { //isClearRefCall(GL) ) {
+            if ( isClearRefCall(GL) ) {
                 refCalls.add(GL);
             } else {
                 AFMatrix.setLikelihoods(GL.getPosteriors(), GL.getSample());
@@ -193,7 +194,8 @@ public abstract class AlleleFrequencyCalculationModel implements Cloneable {
             return false;
 
         double[] likelihoods = GL.getLikelihoods();
-        return ( likelihoods[0] > likelihoods[1] && likelihoods[0] > likelihoods[2]);
+        double refLikelihoodMinusEpsilon = likelihoods[0] - LOG10_REFERENCE_CALL_EPSILON;
+        return ( refLikelihoodMinusEpsilon > likelihoods[1] && refLikelihoodMinusEpsilon > likelihoods[2]);
     }
 
     protected class CalculatedAlleleFrequency {
diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/UnifiedGenotyperEngine.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/UnifiedGenotyperEngine.java
index 9f155e285..4e1b3e5de 100755
--- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/UnifiedGenotyperEngine.java
+++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/genotyper/UnifiedGenotyperEngine.java
@@ -351,7 +351,7 @@ public class UnifiedGenotyperEngine {
             AFline.append(String.format("%.2f\t", ((float)i)/N));
             AFline.append(String.format("%.8f\t", log10AlleleFrequencyPriors[i]));
             if ( log10AlleleFrequencyPosteriors.get()[i] == AlleleFrequencyCalculationModel.VALUE_NOT_CALCULATED)
-                AFline.append("0.0\t");                
+                AFline.append("0.00000000\t");                
             else
                 AFline.append(String.format("%.8f\t", log10AlleleFrequencyPosteriors.get()[i]));
             AFline.append(String.format("%.8f\t", normalizedPosteriors[i]));