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

Restored optimization in Pair HMM only to compute HMM matrices starting in index where haplotypes start to diverge - saves about 15-20% of runtime which is what we lost by disabling banding in latest version, so runtime should be now about the same as what it was before refactoring. Output is bit-true to previous commit

This commit is contained in:
Guillermo del Angel 2012-04-19 19:39:43 -04:00
parent 3fa9089085
commit c44c7b9a97
1 changed files with 8 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
public/java/src/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java
View File

@ -347,7 +347,6 @@ public class PairHMMIndelErrorModel {
// initialize path metric and traceback memories for likelihood computation // initialize path metric and traceback memories for likelihood computation
double[][] matchMetricArray = null, XMetricArray = null, YMetricArray = null; double[][] matchMetricArray = null, XMetricArray = null, YMetricArray = null;
byte[] previousHaplotypeSeen = null; byte[] previousHaplotypeSeen = null;
int startIndexInHaplotype = 0;
final byte[] contextLogGapOpenProbabilities = new byte[readBases.length]; final byte[] contextLogGapOpenProbabilities = new byte[readBases.length];
final byte[] contextLogGapContinuationProbabilities = new byte[readBases.length]; final byte[] contextLogGapContinuationProbabilities = new byte[readBases.length];
@ -376,12 +375,7 @@ public class PairHMMIndelErrorModel {
indStart, indStop, ref.getWindow().getStart(), ref.getWindow().getStop(), startLocationInRefForHaplotypes, stopLocationInRefForHaplotypes, read.getReadLength(), read.getCigar().toString()); indStart, indStop, ref.getWindow().getStart(), ref.getWindow().getStop(), startLocationInRefForHaplotypes, stopLocationInRefForHaplotypes, read.getReadLength(), read.getCigar().toString());
if (indStart < 0 || indStop >= haplotype.getBases().length || indStart > indStop) {
// read spanned more than allowed reference context: we currently can't deal with this
throw new ReviewedStingException("BUG! bad read clipping");
// readLikelihood =0;
} else
{
final byte[] haplotypeBases = Arrays.copyOfRange(haplotype.getBases(), final byte[] haplotypeBases = Arrays.copyOfRange(haplotype.getBases(),
(int)indStart, (int)indStop); (int)indStart, (int)indStop);
@ -394,28 +388,26 @@ public class PairHMMIndelErrorModel {
XMetricArray = new double[X_METRIC_LENGTH][Y_METRIC_LENGTH]; XMetricArray = new double[X_METRIC_LENGTH][Y_METRIC_LENGTH];
YMetricArray = new double[X_METRIC_LENGTH][Y_METRIC_LENGTH]; YMetricArray = new double[X_METRIC_LENGTH][Y_METRIC_LENGTH];
PairHMM.initializeArrays(matchMetricArray, XMetricArray, YMetricArray, X_METRIC_LENGTH);
} }
PairHMM.initializeArrays(matchMetricArray, XMetricArray, YMetricArray, X_METRIC_LENGTH); int startIndexInHaplotype = 0;
if (previousHaplotypeSeen != null)
startIndexInHaplotype = computeFirstDifferingPosition(haplotypeBases, previousHaplotypeSeen);
previousHaplotypeSeen = haplotypeBases.clone();
readLikelihood = pairHMM.computeReadLikelihoodGivenHaplotype(haplotypeBases, readBases, readQuals, readLikelihood = pairHMM.computeReadLikelihoodGivenHaplotype(haplotypeBases, readBases, readQuals,
contextLogGapOpenProbabilities, contextLogGapOpenProbabilities, contextLogGapContinuationProbabilities, contextLogGapOpenProbabilities, contextLogGapOpenProbabilities, contextLogGapContinuationProbabilities,
startIndexInHaplotype, matchMetricArray, XMetricArray, YMetricArray); startIndexInHaplotype, matchMetricArray, XMetricArray, YMetricArray);
previousHaplotypeSeen = haplotypeBases.clone();
/* double l2 = computeReadLikehoodGivenHaplotype(haplotypeBases, readBases, readQuals, contextLogGapOpenProbabilities,
contextLogGapContinuationProbabilities, bandedLikelihoods);
*/
if (DEBUG) { if (DEBUG) {
System.out.println("H:"+new String(haplotypeBases)); System.out.println("H:"+new String(haplotypeBases));
System.out.println("R:"+new String(readBases)); System.out.println("R:"+new String(readBases));
System.out.format("L:%4.2f\n",readLikelihood); System.out.format("L:%4.2f\n",readLikelihood);
// System.out.format("Lorig:%4.2f\n",r2);
System.out.format("StPos:%d\n", startIndexInHaplotype); System.out.format("StPos:%d\n", startIndexInHaplotype);
} }
}
readEl.put(a,readLikelihood); readEl.put(a,readLikelihood);
readLikelihoods[readIdx][j++] = readLikelihood; readLikelihoods[readIdx][j++] = readLikelihood;
} }