From b66c6dced17e891ee19843ed4814f71119841035 Mon Sep 17 00:00:00 2001
From: carneiro <carneiro@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Wed, 4 May 2011 19:26:55 +0000
Subject: [PATCH] - No longer prints out non confident calls (they were leading
 to tables that don't add up and confusing some Pacbio folk). - Added
 sensitivity and Specificity to the report. - With the changes in genotype
 likelihoods, the indel analysis only happens if the BAM file also has an
 extended event. Not great, but at least it's not broken.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5759 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../gatk/walkers/GenotypeAndValidateWalker.java   | 15 +++++++++++++--
 1 file changed, 13 insertions(+), 2 deletions(-)

diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/GenotypeAndValidateWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/GenotypeAndValidateWalker.java
index b2279c522..9e768d166 100755
--- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/GenotypeAndValidateWalker.java
+++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/GenotypeAndValidateWalker.java
@@ -156,6 +156,7 @@ public class GenotypeAndValidateWalker extends RodWalker<GenotypeAndValidateWalk
         if (emitConf >= 0) uac.STANDARD_CONFIDENCE_FOR_EMITTING = emitConf;
         if (callConf >= 0) uac.STANDARD_CONFIDENCE_FOR_CALLING = callConf;
 
+        uac.GLmodel = GenotypeLikelihoodsCalculationModel.Model.SNP;
         snpEngine = new UnifiedGenotyperEngine(getToolkit(), uac);
 
         // Adding the INDEL calling arguments for UG
@@ -205,6 +206,9 @@ public class GenotypeAndValidateWalker extends RodWalker<GenotypeAndValidateWalk
             return counter;
         }
 
+
+        boolean writeVariant = true;
+
         if (bamIsTruth) {
             if (call.confidentlyCalled) {
                 // If truth is a confident REF call
@@ -224,6 +228,7 @@ public class GenotypeAndValidateWalker extends RodWalker<GenotypeAndValidateWalk
             }
             else {
                 counter.nNotConfidentCalls = 1L;
+                writeVariant = false;
             }
         }
         else {
@@ -246,10 +251,11 @@ public class GenotypeAndValidateWalker extends RodWalker<GenotypeAndValidateWalk
             }
             else {
                 counter.nNotConfidentCalls = 1L;
+                writeVariant = false;
             }
         }
 
-        if (vcfWriter != null) {
+        if (vcfWriter != null && writeVariant) {
             if (!vcComp.hasAttribute("callStatus")) {
                 MutableVariantContext mvc = new MutableVariantContext(vcComp);
                 mvc.putAttribute("callStatus", call.isCalledAlt(callConf) ? "ALT" : "REF" );
@@ -284,6 +290,8 @@ public class GenotypeAndValidateWalker extends RodWalker<GenotypeAndValidateWalk
     public void onTraversalDone( CountedData reduceSum ) {
         double ppv = 100 * ((double) reduceSum.nAltCalledAlt /( reduceSum.nAltCalledAlt + reduceSum.nRefCalledAlt));
         double npv = 100 * ((double) reduceSum.nRefCalledRef /( reduceSum.nRefCalledRef + reduceSum.nAltCalledRef));
+        double sensitivity = 100 * ((double) reduceSum.nAltCalledAlt /( reduceSum.nAltCalledAlt + reduceSum.nAltCalledRef));
+        double specificity = (reduceSum.nRefCalledRef + reduceSum.nRefCalledAlt > 0) ? 100 * ((double) reduceSum.nRefCalledRef /( reduceSum.nRefCalledRef + reduceSum.nRefCalledAlt)) : 100;
         logger.info(String.format("Resulting Truth Table Output\n\n" +
                                   "---------------------------------------------------\n" +
                                   "\t\t|\tALT\t|\tREF\t\n"  +
@@ -294,8 +302,11 @@ public class GenotypeAndValidateWalker extends RodWalker<GenotypeAndValidateWalk
                                   "positive predictive value: %f%%\n" +
                                   "negative predictive value: %f%%\n" +
                                   "---------------------------------------------------\n" +
+                                  "sensitivity: %f%%\n" +
+                                  "specificity: %f%%\n" +
+                                  "---------------------------------------------------\n" +
                                   "not confident: %d\n" +
                                   "not covered: %d\n" +
-                                  "---------------------------------------------------\n", reduceSum.nAltCalledAlt, reduceSum.nRefCalledAlt, reduceSum.nAltCalledRef, reduceSum.nRefCalledRef, ppv, npv, reduceSum.nNotConfidentCalls, reduceSum.nUncovered));
+                                  "---------------------------------------------------\n", reduceSum.nAltCalledAlt, reduceSum.nRefCalledAlt, reduceSum.nAltCalledRef, reduceSum.nRefCalledRef, ppv, npv, sensitivity, specificity, reduceSum.nNotConfidentCalls, reduceSum.nUncovered));
     }
 }