From 94e324b8442bd3180b1e4d016a7c4d68c33aa343 Mon Sep 17 00:00:00 2001
From: jmaguire <jmaguire@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Sun, 17 May 2009 13:58:11 +0000
Subject: [PATCH] Write N for the alt allele when we're hom-ref. Stop EM loop
 when we've converged (likelihood[t-1] == likelihood[t]).

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@737 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../playground/gatk/walkers/PoolCaller.java   | 21 +++++++++++++++----
 1 file changed, 17 insertions(+), 4 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 1f014a665..7c2a2c90b 100644
--- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCaller.java
+++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/PoolCaller.java
@@ -114,6 +114,7 @@ public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
         double[] trajectory = new double[MAX_ITERATIONS + 1]; trajectory[0] = EM_alt_freq;
         double[] likelihood_trajectory = new double[MAX_ITERATIONS + 1]; likelihood_trajectory[0] = 0.0;
         boolean is_a_snp = false;
+		char alt = 'N';
         for (int iterations = 0; iterations < MAX_ITERATIONS; iterations++)
         { 
 	        // 6. Re-call from shallow coverage using the estimated frequency as a prior, 
@@ -131,7 +132,9 @@ public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
 	                calls[i] = callers.get(i).map(tracker, ref, contexts[i]);
 	                String genotype = calls[i].genotype();
 
-                    likelihood += calls[i].posterior();
+					if (calls[i].alt != 'N') { alt = calls[i].alt; }
+
+                    likelihood += calls[i].pBest;
 	
                     if (! FRACTIONAL_COUNTS)
                     {
@@ -156,9 +159,18 @@ public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
             }
 	        EM_alt_freq = EM_sum / EM_N;
             trajectory[iterations+1] = EM_alt_freq;
-            likelihood_trajectory[iterations+1] = likelihood/(double)EM_N;
+            likelihood_trajectory[iterations+1] = likelihood;
 
-            //System.out.printf("DBGTRAJ %f %f %f %f\n", EM_sum, EM_N, trajectory[iterations], trajectory[iterations+1]);
+			if (likelihood_trajectory[iterations] == likelihood_trajectory[iterations+1]) { break; }
+
+            //System.out.printf("DBGTRAJ %s %f %f %f %f %f %f\n", 
+			//						context.getLocation(),
+			//						EM_sum, 
+			//						EM_N, 
+			//						trajectory[iterations], 
+			//						trajectory[iterations+1], 
+			//						likelihood_trajectory[iterations],
+			//						likelihood_trajectory[iterations+1]);
         }
 
         // 7. Output some statistics.
@@ -171,7 +183,8 @@ public class PoolCaller extends LocusWalker<AlleleFrequencyEstimate, String>
 			//System.out.printf("DBG %f %f %c %s\n", calls[i].pBest, calls[i].pRef, ref, calls[i].bases);
 		}
 		double discovery_lod = discovery_posterior - discovery_null;
-		discovery_output_file.printf("%s %f %f %f %f\n", context.getLocation(), EM_alt_freq, discovery_posterior, discovery_null, discovery_lod);
+		if (discovery_lod == 0) { alt = 'N'; }
+		discovery_output_file.printf("%s %c %c %f %f %f %f\n", context.getLocation(), ref, alt, EM_alt_freq, discovery_posterior, discovery_null, discovery_lod);
 
         for (int i = 0; i < sample_names.size(); i++)
         {