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

Fix for mismatching indel quals erro: need to adjust for softclips just like we do for bases and normal quals.

This commit is contained in:
Eric Banks 2013-03-04 23:23:18 -05:00
parent 0697594778
commit b715218bfe
1 changed files with 16 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
protected/java/src/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java
View File

@ -343,8 +343,10 @@ public class PairHMMIndelErrorModel {
} }
} }
else { else {
final byte[] readBases = Arrays.copyOfRange(unclippedReadBases,numStartSoftClippedBases, unclippedReadBases.length-numEndSoftClippedBases); final int endOfCopy = unclippedReadBases.length - numEndSoftClippedBases;
final byte[] readQuals = Arrays.copyOfRange(unclippedReadQuals,numStartSoftClippedBases, unclippedReadBases.length-numEndSoftClippedBases); final byte[] readBases = Arrays.copyOfRange(unclippedReadBases, numStartSoftClippedBases, endOfCopy);
final byte[] readQuals = Arrays.copyOfRange(unclippedReadQuals, numStartSoftClippedBases, endOfCopy);
int j=0; int j=0;
byte[] previousHaplotypeSeen = null; byte[] previousHaplotypeSeen = null;
@ -356,6 +358,16 @@ public class PairHMMIndelErrorModel {
getContextHomopolymerLength(readBases,hrunProfile); getContextHomopolymerLength(readBases,hrunProfile);
fillGapProbabilities(hrunProfile, contextLogGapOpenProbabilities, contextLogGapContinuationProbabilities); fillGapProbabilities(hrunProfile, contextLogGapOpenProbabilities, contextLogGapContinuationProbabilities);
// get the base insertion and deletion qualities to use
final byte[] baseInsertionQualities, baseDeletionQualities;
if ( read.hasBaseIndelQualities() ) {
baseInsertionQualities = Arrays.copyOfRange(read.getBaseInsertionQualities(), numStartSoftClippedBases, endOfCopy);
baseDeletionQualities = Arrays.copyOfRange(read.getBaseDeletionQualities(), numStartSoftClippedBases, endOfCopy);
} else {
baseInsertionQualities = contextLogGapOpenProbabilities;
baseDeletionQualities = contextLogGapOpenProbabilities;
}
boolean firstHap = true; boolean firstHap = true;
for (Allele a: haplotypeMap.keySet()) { for (Allele a: haplotypeMap.keySet()) {
@ -385,7 +397,7 @@ public class PairHMMIndelErrorModel {
if (previousHaplotypeSeen == null) { if (previousHaplotypeSeen == null) {
//no need to reallocate arrays for each new haplotype, as length won't change //no need to reallocate arrays for each new haplotype, as length won't change
pairHMM.initialize(X_METRIC_LENGTH, Y_METRIC_LENGTH); pairHMM.initialize(Y_METRIC_LENGTH, X_METRIC_LENGTH);
} }
int startIndexInHaplotype = 0; int startIndexInHaplotype = 0;
@ -394,8 +406,7 @@ public class PairHMMIndelErrorModel {
previousHaplotypeSeen = haplotypeBases.clone(); previousHaplotypeSeen = haplotypeBases.clone();
readLikelihood = pairHMM.computeReadLikelihoodGivenHaplotypeLog10(haplotypeBases, readBases, readQuals, readLikelihood = pairHMM.computeReadLikelihoodGivenHaplotypeLog10(haplotypeBases, readBases, readQuals,
(read.hasBaseIndelQualities() ? read.getBaseInsertionQualities() : contextLogGapOpenProbabilities), baseInsertionQualities, baseDeletionQualities,
(read.hasBaseIndelQualities() ? read.getBaseDeletionQualities() : contextLogGapOpenProbabilities),
contextLogGapContinuationProbabilities, startIndexInHaplotype, firstHap); contextLogGapContinuationProbabilities, startIndexInHaplotype, firstHap);