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= 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