diff --git a/public/java/src/org/broadinstitute/sting/utils/recalibration/BaseRecalibration.java b/public/java/src/org/broadinstitute/sting/utils/recalibration/BaseRecalibration.java index d7598a826..a182e6eda 100644 --- a/public/java/src/org/broadinstitute/sting/utils/recalibration/BaseRecalibration.java +++ b/public/java/src/org/broadinstitute/sting/utils/recalibration/BaseRecalibration.java @@ -28,7 +28,6 @@ package org.broadinstitute.sting.utils.recalibration; import org.broadinstitute.sting.gatk.walkers.bqsr.*; import org.broadinstitute.sting.utils.MathUtils; import org.broadinstitute.sting.utils.QualityUtils; -import org.broadinstitute.sting.utils.exceptions.ReviewedStingException; import org.broadinstitute.sting.utils.sam.GATKSAMRecord; import java.io.File; @@ -42,15 +41,48 @@ import java.util.*; */ public class BaseRecalibration { - private final static String UNRECOGNIZED_REPORT_TABLE_EXCEPTION = "Unrecognized table. Did you add an extra required covariate? This is a hard check that needs propagate through the code"; private final static int MAXIMUM_RECALIBRATED_READ_LENGTH = 5000; private final ReadCovariates readCovariates; private final QuantizationInfo quantizationInfo; // histogram containing the map for qual quantization (calculated after recalibration is done) private final LinkedHashMap> keysAndTablesMap; // quick access reference to the read group table and its key manager + private final KeysAndTables keysAndTables; private final Covariate[] requestedCovariates; // list of all covariates to be used in this calculation + static class KeysAndTables { + + public enum Type { + READ_GROUP_TABLE(0), + QUALITY_SCORE_TABLE(1), + OPTIONAL_COVARIATE_TABLE(2); + + private final int index; + + private Type(int index) { + this.index = index; + } + } + + public final BQSRKeyManager[] managers = new BQSRKeyManager[Type.values().length]; + public final Map[] tables = new Map[Type.values().length]; + + public KeysAndTables(final Map> keysAndTablesMap) { + for (Map.Entry> mapEntry : keysAndTablesMap.entrySet()) { + Type type; + if (mapEntry.getKey().getNumRequiredCovariates() == 1) + type = Type.READ_GROUP_TABLE; + else if (mapEntry.getKey().getNumOptionalCovariates() == 0) + type = Type.QUALITY_SCORE_TABLE; + else + type = Type.OPTIONAL_COVARIATE_TABLE; + managers[type.index] = mapEntry.getKey(); + tables[type.index] = mapEntry.getValue(); + } + } + + } + /** * Constructor using a GATK Report file * @@ -61,6 +93,7 @@ public class BaseRecalibration { RecalibrationReport recalibrationReport = new RecalibrationReport(RECAL_FILE); keysAndTablesMap = recalibrationReport.getKeysAndTablesMap(); + keysAndTables = new KeysAndTables(keysAndTablesMap); requestedCovariates = recalibrationReport.getRequestedCovariates(); quantizationInfo = recalibrationReport.getQuantizationInfo(); if (quantizationLevels == 0) // quantizationLevels == 0 means no quantization, preserve the quality scores @@ -81,6 +114,7 @@ public class BaseRecalibration { protected BaseRecalibration(final QuantizationInfo quantizationInfo, final LinkedHashMap> keysAndTablesMap, final Covariate[] requestedCovariates) { this.quantizationInfo = quantizationInfo; this.keysAndTablesMap = keysAndTablesMap; + keysAndTables = new KeysAndTables(keysAndTablesMap); this.requestedCovariates = requestedCovariates; readCovariates = new ReadCovariates(MAXIMUM_RECALIBRATED_READ_LENGTH, requestedCovariates.length); } @@ -111,7 +145,7 @@ public class BaseRecalibration { } - + /** * Implements a serial recalibration of the reads using the combinational table. * First, we perform a positional recalibration, and then a subsequent dinuc correction. @@ -131,59 +165,60 @@ public class BaseRecalibration { */ protected byte performSequentialQualityCalculation(final long[] key, final EventType errorModel) { - final byte qualFromRead = (byte)(long)key[1]; + final double globalDeltaQ = calculateGlobalDeltaQ(keysAndTables.managers[KeysAndTables.Type.READ_GROUP_TABLE.index], keysAndTables.tables[KeysAndTables.Type.READ_GROUP_TABLE.index], key, errorModel); + final double deltaQReported = calculateDeltaQReported(keysAndTables.managers[KeysAndTables.Type.QUALITY_SCORE_TABLE.index], keysAndTables.tables[KeysAndTables.Type.QUALITY_SCORE_TABLE.index], key, errorModel, globalDeltaQ); + final double deltaQCovariates = calculateDeltaQCovariates(keysAndTables.managers[KeysAndTables.Type.OPTIONAL_COVARIATE_TABLE.index], keysAndTables.tables[KeysAndTables.Type.OPTIONAL_COVARIATE_TABLE.index], key, errorModel, globalDeltaQ, deltaQReported); - double globalDeltaQ = 0.0; - double deltaQReported = 0.0; - double deltaQCovariates = 0.0; - long masterKey; - - for (Map.Entry> mapEntry : keysAndTablesMap.entrySet()) { - final BQSRKeyManager keyManager = mapEntry.getKey(); - final Map table = mapEntry.getValue(); - - switch(keyManager.getNumRequiredCovariates()) { - case 1: // this is the ReadGroup table - masterKey = keyManager.createMasterKey(key, errorModel, -1); - final RecalDatum empiricalQualRG = table.get(masterKey); - if (empiricalQualRG != null) { - final double globalDeltaQEmpirical = empiricalQualRG.getEmpiricalQuality(); - final double aggregrateQReported = empiricalQualRG.getEstimatedQReported(); - globalDeltaQ = globalDeltaQEmpirical - aggregrateQReported; - } - break; - case 2: - final int numOptionalCovariates = keyManager.getNumOptionalCovariates(); - if (numOptionalCovariates == 0) { // this is the QualityScore table - masterKey = keyManager.createMasterKey(key, errorModel, -1); - final RecalDatum empiricalQualQS = table.get(masterKey); - if (empiricalQualQS != null) { - final double deltaQReportedEmpirical = empiricalQualQS.getEmpiricalQuality(); - deltaQReported = deltaQReportedEmpirical - qualFromRead - globalDeltaQ; - } - } - else { // this is the table with all the covariates - for (int i = 0; i < numOptionalCovariates; i++) { - masterKey = keyManager.createMasterKey(key, errorModel, i); - if (masterKey < 0) - continue; - final RecalDatum empiricalQualCO = table.get(masterKey); - if (empiricalQualCO != null) { - final double deltaQCovariateEmpirical = empiricalQualCO.getEmpiricalQuality(); - deltaQCovariates += (deltaQCovariateEmpirical - qualFromRead - (globalDeltaQ + deltaQReported)); - } - } - } - break; - default: - throw new ReviewedStingException(UNRECOGNIZED_REPORT_TABLE_EXCEPTION); - } - } - - double recalibratedQual = qualFromRead + globalDeltaQ + deltaQReported + deltaQCovariates; // calculate the recalibrated qual using the BQSR formula + final byte qualFromRead = (byte)key[1]; + double recalibratedQual = qualFromRead + globalDeltaQ + deltaQReported + deltaQCovariates; // calculate the recalibrated qual using the BQSR formula recalibratedQual = QualityUtils.boundQual(MathUtils.fastRound(recalibratedQual), QualityUtils.MAX_RECALIBRATED_Q_SCORE); // recalibrated quality is bound between 1 and MAX_QUAL return quantizationInfo.getQuantizedQuals().get((int) recalibratedQual); // return the quantized version of the recalibrated quality } + private double calculateGlobalDeltaQ(final BQSRKeyManager keyManager, final Map table, final long[] key, final EventType errorModel) { + double result = 0.0; + + final long masterKey = keyManager.createMasterKey(key, errorModel, -1); + final RecalDatum empiricalQualRG = table.get(masterKey); + if (empiricalQualRG != null) { + final double globalDeltaQEmpirical = empiricalQualRG.getEmpiricalQuality(); + final double aggregrateQReported = empiricalQualRG.getEstimatedQReported(); + result = globalDeltaQEmpirical - aggregrateQReported; + } + + return result; + } + + private double calculateDeltaQReported(final BQSRKeyManager keyManager, final Map table, final long[] key, final EventType errorModel, final double globalDeltaQ) { + double result = 0.0; + + final long masterKey = keyManager.createMasterKey(key, errorModel, -1); + final RecalDatum empiricalQualQS = table.get(masterKey); + if (empiricalQualQS != null) { + final double deltaQReportedEmpirical = empiricalQualQS.getEmpiricalQuality(); + final byte qualFromRead = (byte)key[1]; + result = deltaQReportedEmpirical - qualFromRead - globalDeltaQ; + } + + return result; + } + + private double calculateDeltaQCovariates(final BQSRKeyManager keyManager, final Map table, final long[] key, final EventType errorModel, final double globalDeltaQ, final double deltaQReported) { + double result = 0.0; + + final int numOptionalCovariates = keyManager.getNumOptionalCovariates(); + for (int i = 0; i < numOptionalCovariates; i++) { + final long masterKey = keyManager.createMasterKey(key, errorModel, i); + if (masterKey < 0) + continue; + final RecalDatum empiricalQualCO = table.get(masterKey); + if (empiricalQualCO != null) { + final double deltaQCovariateEmpirical = empiricalQualCO.getEmpiricalQuality(); + final byte qualFromRead = (byte)key[1]; + result += (deltaQCovariateEmpirical - qualFromRead - (globalDeltaQ + deltaQReported)); + } + } + return result; + } }