From 0ee687e49d4bf952e8594d5c51a19fd1ed9ba8dc Mon Sep 17 00:00:00 2001
From: ebanks <ebanks@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Tue, 22 Mar 2011 02:42:28 +0000
Subject: [PATCH] For Mauricio: now, even in GENOTYPE_GIVEN_ALLELES mode, the
 VariantCallContext (which now inherits directly from VC) will report
 reference calls as confidently called if they pass the threshold even if the
 QUAL of the record itself is low because we were forced to have an ALT
 allele.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5485 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../walkers/genotyper/UnifiedGenotyper.java   |  6 +++---
 .../genotyper/UnifiedGenotyperEngine.java     | 15 ++++++++-----
 .../walkers/genotyper/VariantCallContext.java | 21 +++++++++++--------
 3 files changed, 25 insertions(+), 17 deletions(-)

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 e2caf87c8..096aa3af1 100755
--- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java
+++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java
@@ -202,14 +202,14 @@ public class UnifiedGenotyper extends LocusWalker<VariantCallContext, UnifiedGen
         // the base was confidently callable
         sum.nBasesCalledConfidently += value.confidentlyCalled ? 1 : 0;
 
-        // can't make a confident variant call here
-        if ( value.vc == null )
+        // can't make a call here
+        if ( !value.shouldEmit )
             return sum;
 
         try {
             // we are actually making a call
             sum.nCallsMade++;
-            writer.add(value.vc, value.refBase);
+            writer.add(value, value.refBase);
         } catch (IllegalArgumentException e) {
             throw new IllegalArgumentException(e.getMessage() + "; this is often caused by using the --assume_single_sample_reads argument with the wrong sample name");
         }
diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java
index f54509c4a..5f9974123 100755
--- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java
+++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java
@@ -297,7 +297,7 @@ public class UnifiedGenotyperEngine {
         // estimate our confidence in a reference call and return
         if ( vc.getNSamples() == 0 )
             return (UAC.OutputMode != OUTPUT_MODE.EMIT_ALL_SITES ?
-                    estimateReferenceConfidence(stratifiedContexts, genotypePriors.getHeterozygosity(), false, 1.0) :
+                    estimateReferenceConfidence(vc, stratifiedContexts, genotypePriors.getHeterozygosity(), false, 1.0) :
                     new VariantCallContext(generateEmptyContext(tracker, refContext, stratifiedContexts, rawContext), refContext.getBase(), false));
 
         // 'zero' out the AFs (so that we don't have to worry if not all samples have reads at this position)
@@ -336,7 +336,7 @@ public class UnifiedGenotyperEngine {
         if ( UAC.OutputMode != OUTPUT_MODE.EMIT_ALL_SITES && !passesEmitThreshold(phredScaledConfidence, bestAFguess) ) {
             // technically, at this point our confidence in a reference call isn't accurately estimated
             //  because it didn't take into account samples with no data, so let's get a better estimate
-            return estimateReferenceConfidence(stratifiedContexts, genotypePriors.getHeterozygosity(), true, 1.0 - PofF);
+            return estimateReferenceConfidence(vc, stratifiedContexts, genotypePriors.getHeterozygosity(), true, 1.0 - PofF);
         }
 
         // create the genotypes
@@ -423,7 +423,7 @@ public class UnifiedGenotyperEngine {
             vcCall = variantContexts.iterator().next(); // we know the collection will always have exactly 1 element.
         }
 
-        VariantCallContext call = new VariantCallContext(vcCall, passesCallThreshold(phredScaledConfidence));
+        VariantCallContext call = new VariantCallContext(vcCall, confidentlyCalled(phredScaledConfidence, PofF));
         call.setRefBase(refContext.getBase());
         return call;
     }
@@ -491,7 +491,7 @@ public class UnifiedGenotyperEngine {
             AFs[i] = AlleleFrequencyCalculationModel.VALUE_NOT_CALCULATED;
     }
 
-    private VariantCallContext estimateReferenceConfidence(Map<String, StratifiedAlignmentContext> contexts, double theta, boolean ignoreCoveredSamples, double initialPofRef) {
+    private VariantCallContext estimateReferenceConfidence(VariantContext vc, Map<String, StratifiedAlignmentContext> contexts, double theta, boolean ignoreCoveredSamples, double initialPofRef) {
         if ( contexts == null )
             return null;
 
@@ -518,7 +518,7 @@ public class UnifiedGenotyperEngine {
             P_of_ref *= 1.0 - (theta / 2.0) * MathUtils.binomialProbability(0, depth, 0.5);
         }
 
-        return new VariantCallContext(QualityUtils.phredScaleErrorRate(1.0 - P_of_ref) >= UAC.STANDARD_CONFIDENCE_FOR_CALLING);
+        return new VariantCallContext(vc, QualityUtils.phredScaleErrorRate(1.0 - P_of_ref) >= UAC.STANDARD_CONFIDENCE_FOR_CALLING, false);
     }
 
     protected void printVerboseData(String pos, VariantContext vc, double PofF, double phredScaledConfidence, double[] normalizedPosteriors) {
@@ -651,6 +651,11 @@ public class UnifiedGenotyperEngine {
         return conf >= UAC.STANDARD_CONFIDENCE_FOR_CALLING;
     }
 
+    protected boolean confidentlyCalled(double conf, double PofF) {
+        return conf >= UAC.STANDARD_CONFIDENCE_FOR_CALLING ||
+                (UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES && QualityUtils.phredScaleErrorRate(PofF) >= UAC.STANDARD_CONFIDENCE_FOR_CALLING);
+    }
+
     protected void computeAlleleFrequencyPriors(int N) {
         // calculate the allele frequency priors for 1-N
         double sum = 0.0;
diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/VariantCallContext.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/VariantCallContext.java
index f280b75b7..16dac6fbe 100755
--- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/VariantCallContext.java
+++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/VariantCallContext.java
@@ -35,29 +35,32 @@ import org.broad.tribble.util.variantcontext.VariantContext;
  *
  * Useful helper class to communicate the results of calculateGenotype to framework
  */
-public class VariantCallContext {
-    public VariantContext vc = null;
+public class VariantCallContext extends VariantContext {
     public byte refBase;
 
     // Was the site called confidently, either reference or variant?
     public boolean confidentlyCalled = false;
 
+    // Should this site be emitted?
+    public boolean shouldEmit = true;
+
     VariantCallContext(VariantContext vc, boolean confidentlyCalledP) {
-        this.vc = vc;
+        super(vc);
         this.confidentlyCalled = confidentlyCalledP;
     }
 
+    VariantCallContext(VariantContext vc, boolean confidentlyCalledP, boolean shouldEmit) {
+        super(vc);
+        this.confidentlyCalled = confidentlyCalledP;
+        this.shouldEmit = shouldEmit;
+    }
+
     VariantCallContext(VariantContext vc, byte ref, boolean confidentlyCalledP) {
-        this.vc = vc;
+        super(vc);
         this.refBase = ref;
         this.confidentlyCalled = confidentlyCalledP;
     }
 
-    // blank variant context => we're a ref site
-    VariantCallContext(boolean confidentlyCalledP) {
-        this.confidentlyCalled = confidentlyCalledP;
-    }
-
     public void setRefBase(byte ref) {
         this.refBase = ref;
     }