From 2ef2c9e121d3973ad3cd4268e219b0b1bfdb7bd2 Mon Sep 17 00:00:00 2001
From: kiran <kiran@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Fri, 10 Apr 2009 04:17:53 +0000
Subject: [PATCH] Fixed an issue wherein the SQ field was only being pulled
 from the first read of the pileup, no matter what.  Fixed an issue wherein
 Andrew enumerates his bases as A:0, C:1, T:2, G:3, and Kiran's QualityUtils
 methods enumerate bases as A:0, C:1, G:2, T:3 (we should standardize this). 
 Fixed an issue wherein the remaining probability was being divided by 3
 rather than 2 when four-base probs are enabled.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@356 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../gatk/walkers/AlleleFrequencyWalker.java    | 18 ++++++++++++++----
 1 file changed, 14 insertions(+), 4 deletions(-)

diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java
index 003b944cd..91ddf973c 100755
--- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java
+++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java
@@ -150,11 +150,12 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
 
             // First, set the called base to the correct quality
             int callednum = nuc2num[read.getReadBases()[offset]]; // number of the called base (0,1,2,3)
+
             quals[i][callednum] = 1 - Math.pow(10.0, (double)read.getBaseQualities()[offset] / -10.0);
 
             // For other 3 bases, check if 2-base probs exist
             assert (reads.size() > 0);
-            Object SQ_field = reads.get(0).getAttribute("SQ");
+            Object SQ_field = reads.get(i).getAttribute("SQ");
             if (SQ_field == null || FORCE_1BASE_PROBS) {
                 // Set all nonref qual scores to their share of the remaining probality not "used" by the reference base's qual
                 double nonref_quals = (1.0 - quals[i][callednum]) / 3;
@@ -168,13 +169,22 @@ public class AlleleFrequencyWalker extends LocusWalker<AlleleFrequencyEstimate,
                 //System.out.printf("SAM record: %s\n", read.format());
 
                 byte hex_qual = hex_quals[offset];
-                int called2num = QualityUtils.compressedQualityToBaseIndex(hex_qual);
+                int baseIndex = QualityUtils.compressedQualityToBaseIndex(hex_qual);
+                // Todo: Andrew indexes his nucleotides as 0:A, 1:C, 2:T, 3:G.  Kiran indexes his nucleotides as 0:A, 1:C, 2:G, 3:T.  This must be reconciled.  Here, I just translate appropriately.
+                int called2num = 0;
+                switch (baseIndex) {
+                    case 0: called2num = 0; break;
+                    case 1: called2num = 1; break;
+                    case 2: called2num = 3; break;
+                    case 3: called2num = 2; break;
+                    default: assert(baseIndex >= 0 && baseIndex < 4); break;
+                }
+
                 double qual2 = QualityUtils.compressedQualityToProb(hex_qual);
                 //System.out.printf("2ND %x %d %f\n", hex_qual, called2num, qual2);
                 quals[i][called2num] = qual2;
 
-                // 
-                double nonref_quals = (1.0 - quals[i][callednum] - quals[i][called2num]) / 3;
+                double nonref_quals = (1.0 - quals[i][callednum] - quals[i][called2num]) / 2;
                 for (int b=0; b<4; b++)
                     if (b != callednum && b != called2num)
                         quals[i][b] = nonref_quals;