From 5d39cd5db8f52bf0f29a786fbb8f7ea0b2ccb105 Mon Sep 17 00:00:00 2001
From: rpoplin <rpoplin@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Thu, 8 Jul 2010 22:44:10 +0000
Subject: [PATCH] Added --fdr_filter_level to ApplyVariantCuts so that you can
 create beautiful tranche plots and also decide which tranche level to filter
 at. The previous version always filtered at the smallest tranche. The tranche
 filter names are appropriately added to the VCF header.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@3747 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../variantrecalibration/ApplyVariantCuts.java | 18 +++++++++++-------
 .../VariantRecalibrator.java                   | 16 +++++++---------
 2 files changed, 18 insertions(+), 16 deletions(-)

diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/ApplyVariantCuts.java b/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/ApplyVariantCuts.java
index f70512f10..b58e5ce18 100755
--- a/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/ApplyVariantCuts.java
+++ b/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/ApplyVariantCuts.java
@@ -65,6 +65,9 @@ public class ApplyVariantCuts extends RodWalker<Integer, Integer> {
     private String TRANCHE_FILENAME = "optimizer.dat.tranches";
     @Argument(fullName="outputVCFFile", shortName="outputVCF", doc="The output filtered VCF file", required=true)
     private String OUTPUT_FILENAME = "optimizer.vcf";
+    @Argument(fullName="fdr_filter_level", shortName="fdr_filter_level", doc="The FDR level at which to start filtering.", required=true)
+    private double FDR_FILTER_LEVEL = 0.0;
+
 
     /////////////////////////////
     // Private Member Variables
@@ -87,8 +90,10 @@ public class ApplyVariantCuts extends RodWalker<Integer, Integer> {
             for( final String line : new XReadLines(new File( TRANCHE_FILENAME )) ) {
                 if( !firstLine ) {
                     final String[] vals = line.split(",");
-                    qCuts.add(Double.parseDouble(vals[2]));
-                    filterName.add(vals[4]);
+                    if(Double.parseDouble(vals[0]) >= FDR_FILTER_LEVEL) {
+                        qCuts.add(Double.parseDouble(vals[2]));
+                        filterName.add(vals[4]);
+                    }
                 }
                 firstLine = false;
             }
@@ -119,11 +124,10 @@ public class ApplyVariantCuts extends RodWalker<Integer, Integer> {
             }
         }
 
-        // BUGBUG : why doesn't this work?  VCFFilterHeaderLine is protected?
-        //for( int iii = 1; iii < filterName.size(); iii++ ) {
-        //    hInfo.add(new VCFFilterHeaderLine(filterName.get(iii)), String.format("FDR tranche level at qual " + qCuts.get(iii)));
-        //}
-        //hInfo.add(new VCFFilterHeaderLine(filterName.get(0)+"+"), String.format("FDR tranche level at qual > " + qCuts.get(filterName.size()-1)));
+        for( int iii = 1; iii < filterName.size(); iii++ ) {
+            hInfo.add(new VCFFilterHeaderLine(filterName.get(iii), String.format("FDR tranche level at qual " + qCuts.get(iii))));
+        }
+        hInfo.add(new VCFFilterHeaderLine(filterName.get(0)+"+", String.format("FDR tranche level at qual > " + qCuts.get(filterName.size()-1))));
 
         final VCFHeader vcfHeader = new VCFHeader(hInfo, samples);
         vcfWriter.writeHeader(vcfHeader);
diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantRecalibrator.java b/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantRecalibrator.java
index ec777d751..652c0005c 100755
--- a/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantRecalibrator.java
+++ b/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantRecalibrator.java
@@ -155,15 +155,13 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
 
         // Set up default values for the FDR tranches if necessary
         if( FDR_TRANCHES == null ) {
-            FDR_TRANCHES = new Double[3];
-            FDR_TRANCHES[0] = 12.5;
-            FDR_TRANCHES[1] = 20.0;
-            FDR_TRANCHES[2] = 30.0;
-            //FDR_TRANCHES[0] = 1.0;
-            //FDR_TRANCHES[1] = 5.0;
-            //FDR_TRANCHES[2] = 10.0;
-            //FDR_TRANCHES[3] = 20.0;
-            //FDR_TRANCHES[4] = 30.0;
+            FDR_TRANCHES = new Double[6];
+            FDR_TRANCHES[0] = 0.1;
+            FDR_TRANCHES[1] = 1.0;
+            FDR_TRANCHES[2] = 5.0;
+            FDR_TRANCHES[3] = 12.5;
+            FDR_TRANCHES[4] = 25.0;
+            FDR_TRANCHES[5] = 35.0;
         }
     }