From 7ecc43e9a7d8e4246016ca7e3ad021317b963716 Mon Sep 17 00:00:00 2001
From: depristo <depristo@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Mon, 29 Jun 2009 18:57:42 +0000
Subject: [PATCH] Fixed subtle null ptr exception discovered by Kiran.  Now
 deals with the rare situation where you have only say Q28 bases at dbSNP
 sites, so you fail in the Table recalibration step with a null pointer error
 into the data structure indexed by quality score.  If you are Q score above
 those seen before you aren't modified in any way.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1118 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../TableRecalibrationWalker.java             | 40 +++++++++----------
 1 file changed, 20 insertions(+), 20 deletions(-)

diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/recalibration/TableRecalibrationWalker.java b/java/src/org/broadinstitute/sting/gatk/walkers/recalibration/TableRecalibrationWalker.java
index ca1aff4b6..f8847ba80 100755
--- a/java/src/org/broadinstitute/sting/gatk/walkers/recalibration/TableRecalibrationWalker.java
+++ b/java/src/org/broadinstitute/sting/gatk/walkers/recalibration/TableRecalibrationWalker.java
@@ -54,7 +54,7 @@ public class TableRecalibrationWalker extends ReadWalker<SAMRecord, SAMFileWrite
 
     private static Logger logger = Logger.getLogger(TableRecalibrationWalker.class);
 
-    private static String VERSION = "0.2.1";
+    private static String VERSION = "0.2.2";
 
     private final static boolean DEBUG = false;
 
@@ -266,11 +266,14 @@ public class TableRecalibrationWalker extends ReadWalker<SAMRecord, SAMFileWrite
      * @return
      */
     public byte[] recalibrateBasesAndQuals(final String readGroup, byte[] bases, byte[] quals) throws StingException {
+        if ( readGroup == null ) { throw new StingException("BUG: read group is null"); }
+        if ( bases == null ) { throw new StingException("BUG: bases array is null"); }
+        if ( bases.length == 0) { throw new StingException("BUG: bases array is size 0"); }
+        if ( quals == null ) { throw new StingException("BUG: quals array is null"); }
+
         byte[] recalQuals = new byte[quals.length];
         RecalMapping mapper = cache.get(readGroup);
-
-        //if ( mapper == null && DEBUG_ME )
-        //    return recalQuals;
+        if ( mapper == null ) { throw new StingException(String.format("BUG: couldn't find RecalMapping for readgroup %s", readGroup)); }
 
         recalQuals[0] = quals[0];               // can't change the first -- no dinuc
         for ( int cycle = 1; cycle < bases.length; cycle++ ) { // skip first and last base, qual already set because no dinuc
@@ -489,11 +492,19 @@ class SerialRecalMapping implements RecalMapping {
         int pos = manager.canonicalPos(cycle);
         int index = this.manager.getDinucIndex(prevBase, base);
 
-        double deltaQual = deltaQualMap[qual];
-        double deltaQPos = pos >= deltaQPosMap.length ? 0.0 : deltaQPosMap[pos][qual];  // == length indices no data for last base
-        double deltaQDinuc = index == -1 ? 0.0 : deltaQDinucMap[index][qual];           // -1 indices N in prev or current base
-        double newQual = qual + globalDeltaQ + deltaQual + deltaQPos + deltaQDinuc;
-        byte newQualByte = QualityUtils.boundQual((int)Math.round(newQual), QualityUtils.MAX_REASONABLE_Q_SCORE);
+        byte newQualByte = qual;
+        if ( qual < deltaQualMap.length ) {
+            // it's possible to see a qual not in the recal table -- i.e., we only see Q28 at a dbSNP site
+            double deltaQual = deltaQualMap[qual];
+            double deltaQPos = pos >= deltaQPosMap.length ? 0.0 : deltaQPosMap[pos][qual];  // == length indices no data for last base
+            double deltaQDinuc = index == -1 ? 0.0 : deltaQDinucMap[index][qual];           // -1 indices N in prev or current base
+            double newQual = qual + globalDeltaQ + deltaQual + deltaQPos + deltaQDinuc;
+            newQualByte = QualityUtils.boundQual((int)Math.round(newQual), QualityUtils.MAX_REASONABLE_Q_SCORE);
+
+            if ( newQualByte <= 0 && newQualByte >= QualityUtils.MAX_REASONABLE_Q_SCORE )
+                throw new RuntimeException(String.format("Illegal base quality score calculated: %s %c%c %d %d => %d + %.2f + %.2f + %.2f = %d",
+                        readGroup, prevBase, base, cycle, qual, qual, globalDeltaQ, deltaQPos, deltaQDinuc, newQualByte));
+        }
 
         //if ( printStateP(pos, RecalData.dinucIndex2bases(index), qual) )
         //    System.out.printf("%s %c%c %d %d => %d + %.2f + %.2f + %.2f + %.2f = %d%n",
@@ -501,17 +512,6 @@ class SerialRecalMapping implements RecalMapping {
         //           qual, globalDeltaQ, deltaQual, deltaQPos, deltaQDinuc,
         //            newQualByte);
 
-
-        //byte combiQualByte = combiMap.getNewQual(readGroup, prevBase, base, cycle, qual);
-        //System.out.printf("%s %c%c %d %d => %d + %.2f + %.2f + %.2f = %d vs %d (%d)%n",
-        //        readGroup, prevBase, base, cycle, qual,
-        //        qual, globalDeltaQ, deltaQPos, deltaQDinuc,
-        //        newQualByte, combiQualByte, newQualByte - combiQualByte);
-
-        if ( newQualByte <= 0 && newQualByte >= QualityUtils.MAX_REASONABLE_Q_SCORE )
-            throw new RuntimeException(String.format("Illegal base quality score calculated: %s %c%c %d %d => %d + %.2f + %.2f + %.2f = %d",
-                    readGroup, prevBase, base, cycle, qual, qual, globalDeltaQ, deltaQPos, deltaQDinuc, newQualByte));
-
         return newQualByte;
     }
 }
\ No newline at end of file