zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Refactoring mostly for readability (and small performance improvement)

This commit is contained in:
Eric Banks 2012-06-15 10:36:41 -04:00
parent c54e84e739
commit 0c218e4822
1 changed files with 87 additions and 52 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

135
public/java/src/org/broadinstitute/sting/utils/recalibration/BaseRecalibration.java
View File

@ -28,7 +28,6 @@ package org.broadinstitute.sting.utils.recalibration;
import org.broadinstitute.sting.gatk.walkers.bqsr.*; import org.broadinstitute.sting.gatk.walkers.bqsr.*;
import org.broadinstitute.sting.utils.MathUtils; import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.QualityUtils; import org.broadinstitute.sting.utils.QualityUtils;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord; import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
import java.io.File; import java.io.File;
@ -42,15 +41,48 @@ import java.util.*;
*/ */
public class BaseRecalibration { 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 static int MAXIMUM_RECALIBRATED_READ_LENGTH = 5000;
private final ReadCovariates readCovariates; private final ReadCovariates readCovariates;
private final QuantizationInfo quantizationInfo; // histogram containing the map for qual quantization (calculated after recalibration is done) private final QuantizationInfo quantizationInfo; // histogram containing the map for qual quantization (calculated after recalibration is done)
private final LinkedHashMap<BQSRKeyManager, Map<Long, RecalDatum>> keysAndTablesMap; // quick access reference to the read group table and its key manager private final LinkedHashMap<BQSRKeyManager, Map<Long, RecalDatum>> 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 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<Long, RecalDatum>[] tables = new Map[Type.values().length];
public KeysAndTables(final Map<BQSRKeyManager, Map<Long, RecalDatum>> keysAndTablesMap) {
for (Map.Entry<BQSRKeyManager, Map<Long, RecalDatum>> 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 * Constructor using a GATK Report file
* *
@ -61,6 +93,7 @@ public class BaseRecalibration {
RecalibrationReport recalibrationReport = new RecalibrationReport(RECAL_FILE); RecalibrationReport recalibrationReport = new RecalibrationReport(RECAL_FILE);
keysAndTablesMap = recalibrationReport.getKeysAndTablesMap(); keysAndTablesMap = recalibrationReport.getKeysAndTablesMap();
keysAndTables = new KeysAndTables(keysAndTablesMap);
requestedCovariates = recalibrationReport.getRequestedCovariates(); requestedCovariates = recalibrationReport.getRequestedCovariates();
quantizationInfo = recalibrationReport.getQuantizationInfo(); quantizationInfo = recalibrationReport.getQuantizationInfo();
if (quantizationLevels == 0) // quantizationLevels == 0 means no quantization, preserve the quality scores 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<BQSRKeyManager, Map<Long, RecalDatum>> keysAndTablesMap, final Covariate[] requestedCovariates) { protected BaseRecalibration(final QuantizationInfo quantizationInfo, final LinkedHashMap<BQSRKeyManager, Map<Long, RecalDatum>> keysAndTablesMap, final Covariate[] requestedCovariates) {
this.quantizationInfo = quantizationInfo; this.quantizationInfo = quantizationInfo;
this.keysAndTablesMap = keysAndTablesMap; this.keysAndTablesMap = keysAndTablesMap;
keysAndTables = new KeysAndTables(keysAndTablesMap);
this.requestedCovariates = requestedCovariates; this.requestedCovariates = requestedCovariates;
readCovariates = new ReadCovariates(MAXIMUM_RECALIBRATED_READ_LENGTH, requestedCovariates.length); readCovariates = new ReadCovariates(MAXIMUM_RECALIBRATED_READ_LENGTH, requestedCovariates.length);
} }
@ -131,59 +165,60 @@ public class BaseRecalibration {
*/ */
protected byte performSequentialQualityCalculation(final long[] key, final EventType errorModel) { 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);
double globalDeltaQ = 0.0; 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 deltaQReported = 0.0;
double deltaQCovariates = 0.0;
long masterKey;
for (Map.Entry<BQSRKeyManager, Map<Long, RecalDatum>> mapEntry : keysAndTablesMap.entrySet()) {
final BQSRKeyManager keyManager = mapEntry.getKey();
final Map<Long, RecalDatum> 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);
}
}
final byte qualFromRead = (byte)key[1];
double recalibratedQual = qualFromRead + globalDeltaQ + deltaQReported + deltaQCovariates; // calculate the recalibrated qual using the BQSR formula 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 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 return quantizationInfo.getQuantizedQuals().get((int) recalibratedQual); // return the quantized version of the recalibrated quality
} }
private double calculateGlobalDeltaQ(final BQSRKeyManager keyManager, final Map<Long, RecalDatum> 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<Long, RecalDatum> 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<Long, RecalDatum> 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;
}
} }