From 7084ecdeb67b96dd4015c7e78a3d65fddecd7d4f Mon Sep 17 00:00:00 2001
From: jmaguire <jmaguire@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Wed, 13 May 2009 15:50:48 +0000
Subject: [PATCH] a few changes; checked in to allow debugging.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@688 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../playground/gatk/walkers/PoolCaller.java   | 16 +++----
 .../gatk/walkers/SingleSampleGenotyper.java   | 42 +++++++++----------
 2 files changed, 28 insertions(+), 30 deletions(-)

diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCaller.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCaller.java
index 6c4721b69..7789a0343 100644
--- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCaller.java
+++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCaller.java
@@ -19,6 +19,9 @@ import java.util.List;
 import java.util.Random;
 import java.io.*;
 
+// Draft iterative pooled caller
+// j.maguire 4-27-2009
+
 public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String> 
 {
     List<SingleSampleGenotyper> callers = null;
@@ -28,11 +31,11 @@ public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
     @Argument(required=false, shortName="max_iterations", doc="Maximum number of iterations for EM") public int MAX_ITERATIONS = 10;
     @Argument(fullName="lodThreshold", shortName="lod", required=false, doc="lod threshold for outputting individual genotypes")       public Double lodThreshold = 2.0;
     @Argument(fullName="discovery_output", shortName="discovery_output", required=true, doc="file to write SNP discovery output to")       public String DISCOVERY_OUTPUT;
+    @Argument(fullName="individual_output_prefix", shortName="individual_output_prefix", required=true, doc="prefix to write individual SNP calls to") public String INDIVIDUAL_OUTPUT_PREFIX;
+
     
     private Random random;
-
     private SAMFileHeader header;
-
 	private PrintStream discovery_output_file;
 
     public void initialize() 
@@ -69,27 +72,26 @@ public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
             System.out.println("SAMPLE: " + sample_name);
 
             SingleSampleGenotyper caller = new SingleSampleGenotyper();
-            caller.metricsFileName = "/dev/null";
+			caller.callsFileName = INDIVIDUAL_OUTPUT_PREFIX + "." + sample_name + ".calls";
+			caller.metricsFileName = INDIVIDUAL_OUTPUT_PREFIX + "." + sample_name + ".metrics";
             caller.lodThreshold = lodThreshold;
             caller.fourBaseMode = false;
             caller.printMetrics = false;
-            caller.initialize();
+            //caller.initialize();
             callers.add(caller);
         } 
     }
 
     public AlleleFrequencyEstimate map(RefMetaDataTracker tracker, char ref, LocusContext context) 
     {
-		if (ref == 'N') { return null; }
+
 
         // 1. seperate each context.
         LocusContext[] contexts = new LocusContext[sample_names.size()];
         for (int i = 0; i < sample_names.size(); i++)
         {
             contexts[i] = filterLocusContext(context, sample_names.get(i), 0);
-            //System.out.printf("DEPTH %s %d\n", sample_names.get(i), contexts[i].getReads().size());
         }
-        //System.out.printf("DEPTH %s %d\n", "TOTAL", context.getReads().size());
 
         // EM Loop:
 	    AlleleFrequencyEstimate[] calls = null;
diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/SingleSampleGenotyper.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/SingleSampleGenotyper.java
index 687d8d40b..2da90e3b1 100644
--- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/SingleSampleGenotyper.java
+++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/SingleSampleGenotyper.java
@@ -15,6 +15,7 @@ import org.broadinstitute.sting.utils.QualityUtils;
 
 import net.sf.samtools.SAMRecord;
 
+import java.io.*;
 import java.util.*;
 
 // Draft single sample genotyper
@@ -22,6 +23,7 @@ import java.util.*;
 
 public class SingleSampleGenotyper extends LocusWalker<AlleleFrequencyEstimate, String>  
 {
+    @Argument(fullName="calls",        shortName="calls",        doc="File to write SNP calls to",      required=true) public String callsFileName;
     @Argument(fullName="metrics",      shortName="met",          doc="metrics",      required=false) public String metricsFileName = "/dev/null";
     @Argument(fullName="metInterval",  shortName="mi",           doc="metInterval",  required=false)     public Integer metricsInterval = 50000;
     @Argument(fullName="printMetrics", shortName="printMetrics", doc="printMetrics", required=false)      public Boolean printMetrics = true;
@@ -34,38 +36,32 @@ public class SingleSampleGenotyper extends LocusWalker<AlleleFrequencyEstimate,
     @Argument(fullName="qHomNonRef",   shortName="qHomNonRef",   doc="qHomNonRef",   required=false)      public Double qHomNonRef = 0.97;
 
     public AlleleMetrics metrics;
+	public PrintStream   calls_file;
 
 	public String sample_name;
     
     public boolean filter(RefMetaDataTracker tracker, char ref, LocusContext context) { return true; }
     public boolean requiresReads() { return true; }
-    public void initialize() { metrics = new AlleleMetrics(metricsFileName, lodThreshold); sample_name = null; }
+    public void initialize() 
+	{ 
+		try 
+		{
+			metrics = new AlleleMetrics(metricsFileName, lodThreshold); 
+			sample_name = null; 
+			calls_file = new PrintStream(callsFileName);
+		}
+		catch (Exception e)
+		{
+			e.printStackTrace(); 
+			System.exit(-1);
+		}
+	}
 
     public AlleleFrequencyEstimate map(RefMetaDataTracker tracker, char ref, LocusContext context) {
         String rodString = getRodString(tracker);
 
 		if (ref == 'N') { return null; }
 
-        /*
-        AlleleFrequencyEstimate freq;
-        if (fourBaseMode) {
-            if (retest) {
-                // Compute single quality score genotype likelihoods
-                freq = getOneProbAlleleFrequency(ref, context, rodString);
-
-                if (freq.qhat > 0.1) {
-                    // If we found a putative variant, recompute four-base prob genotype likelihoods
-                    freq = getFourProbAlleleFrequency(ref, context, rodString, tracker);
-                }
-            } else {
-                freq = getFourProbAlleleFrequency(ref, context, rodString, tracker);
-            }
-        } else  {
-            // Compute single quality score genotype likelihoods
-            freq = getOneProbAlleleFrequency(ref, context, rodString);
-        }
-        */
-		
 		if (context.getReads().size() != 0)
 		{
 			SAMRecord read = context.getReads().get(0);
@@ -367,7 +363,7 @@ public class SingleSampleGenotyper extends LocusWalker<AlleleFrequencyEstimate,
             // No longer hom-ref, so output a ref line.
             double lod = confident_ref_interval_LOD_sum / confident_ref_interval_length;
 
-            out.format("%s\tCALLER\tREFERENCE\t%d\t%d\t%f\t.\t.\tLENGTH %d\n",
+            calls_file.format("%s\tCALLER\tREFERENCE\t%d\t%d\t%f\t.\t.\tLENGTH %d\n",
                             confident_ref_interval_contig,
                             confident_ref_interval_start,
                             last_position_considered,
@@ -399,7 +395,7 @@ public class SingleSampleGenotyper extends LocusWalker<AlleleFrequencyEstimate,
         last_position_considered = current_offset;
         
         if (alleleFreq.lodVsRef >= 5) {
-            out.print(alleleFreq.asGFFString());
+            calls_file.print(alleleFreq.asGFFString());
 
             /*
             String gtype = genotypeTypeString(alleleFreq.qstar, alleleFreq.N);