From 9decd20f46e3d96420b4c8729fba6523786df104 Mon Sep 17 00:00:00 2001
From: depristo <depristo@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Fri, 29 Jan 2010 15:36:12 +0000
Subject: [PATCH] Fix to priors to allow lower het values for mouse guys; no
 intergration test changes

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@2734 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../genotyper/DiploidGenotypePriors.java      | 23 +++++++++++--------
 .../walkers/genotyper/UnifiedGenotyper.java   |  2 +-
 .../papergenotyper/GATKPaperGenotyper.java    |  2 +-
 3 files changed, 15 insertions(+), 12 deletions(-)

diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypePriors.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypePriors.java
index 0973341f4..5fe2a8390 100755
--- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypePriors.java
+++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypePriors.java
@@ -17,9 +17,11 @@ public class DiploidGenotypePriors {
     public static final double YRI_HETEROZYGOSITY = 1.0 / 850;
 
     /**
-     * Default value of the prob of seeing a true B/C het when the reference is A
+     * Default value of the prob of seeing a reference error.  Used to calculation the
+     * chance of seeing a true B/C het when the reference is A, which we assume is the product
+     * of the ref error rate and the het. Default value is Q60
      */
-    public static final double PROB_OF_TRISTATE_GENOTYPE = 1e-5;
+    public static final double PROB_OF_REFERENCE_ERROR = 1e-6;  // the reference is
 
     private final static double[] flatPriors = new double[DiploidGenotype.values().length];
 
@@ -52,7 +54,7 @@ public class DiploidGenotypePriors {
             double[] vals = heterozygosity2DiploidProbabilities(heterozygosity);
             priors = getReferenceIndependentPriors(ref, vals[0], vals[1], vals[2]);
         } else {
-            priors = getReferencePolarizedPriors(ref, heterozygosity, PROB_OF_TRISTATE_GENOTYPE);
+            priors = getReferencePolarizedPriors(ref, heterozygosity, PROB_OF_REFERENCE_ERROR);
         }
     }
 
@@ -198,16 +200,17 @@ public class DiploidGenotypePriors {
      *
      * @param ref
      * @param heterozyosity
-     * @param pTriStateGenotype
+     * @param pRefError
      */
-    public static double[] getReferencePolarizedPriors(char ref, double heterozyosity, double pTriStateGenotype ) {
-        if ( ! MathUtils.isBounded(pTriStateGenotype, 0.0, 1.0) ) {
-            throw new RuntimeException(String.format("p Tristate genotype is out of bounds %f", pTriStateGenotype));
+    public static double[] getReferencePolarizedPriors(char ref, double heterozyosity, double pRefError ) {
+        if ( ! MathUtils.isBounded(pRefError, 0.0, 0.01) ) {
+            throw new RuntimeException(String.format("BUG: p Reference error is out of bounds (0.0 - 0.01) is allow range %f", pRefError));
         }
 
-        if ( pTriStateGenotype >= heterozyosity ) {
-            throw new RuntimeException(String.format("p Tristate genotype %f is greater than the heterozygosity %f", pTriStateGenotype, heterozyosity));
-        }
+        double pTriStateGenotype = heterozyosity * PROB_OF_REFERENCE_ERROR;
+//        if ( pTriStateGenotype >= heterozyosity ) {
+//            throw new RuntimeException(String.format("p Tristate genotype %f is greater than the heterozygosity %f", pTriStateGenotype, heterozyosity));
+//        }
 
         double pHomRef = heterozygosity2HomRefProbability(heterozyosity);
         double pHet    = heterozygosity2HetProbability(heterozyosity);
diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java
index 1dbc28f3f..92c853f47 100755
--- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java
+++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java
@@ -297,7 +297,7 @@ public class UnifiedGenotyper extends LocusWalker<VariantCallContext, UnifiedGen
         if ( stratifiedContexts == null )
             return null;
 
-        DiploidGenotypePriors priors = new DiploidGenotypePriors(ref, UAC.heterozygosity, DiploidGenotypePriors.PROB_OF_TRISTATE_GENOTYPE);
+        DiploidGenotypePriors priors = new DiploidGenotypePriors(ref, UAC.heterozygosity, DiploidGenotypePriors.PROB_OF_REFERENCE_ERROR);
         VariantCallContext call = gcm.get().callLocus(tracker, ref, rawContext.getLocation(), stratifiedContexts, priors);
 
         // annotate the call, if possible
diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/papergenotyper/GATKPaperGenotyper.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/papergenotyper/GATKPaperGenotyper.java
index 2fd3ba669..ac40f4616 100644
--- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/papergenotyper/GATKPaperGenotyper.java
+++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/papergenotyper/GATKPaperGenotyper.java
@@ -45,7 +45,7 @@ public class GATKPaperGenotyper extends LocusWalker<SimpleCall, Integer> impleme
         ReadBackedPileup pileup = context.getPileup();
         double likelihoods[] = DiploidGenotypePriors.getReferencePolarizedPriors(ref.getBase(),
                                                                                  DiploidGenotypePriors.HUMAN_HETEROZYGOSITY,
-                                                                                 DiploidGenotypePriors.PROB_OF_TRISTATE_GENOTYPE);
+                                                                                 0.01);
         // get the bases and qualities from the pileup
         byte bases[] = pileup.getBases();
         byte quals[] = pileup.getQuals();