diff --git a/java/src/org/broadinstitute/sting/alignment/bwa/BWAAligner.java b/java/src/org/broadinstitute/sting/alignment/bwa/BWAAligner.java index 115fb08f8..3bb27eec2 100644 --- a/java/src/org/broadinstitute/sting/alignment/bwa/BWAAligner.java +++ b/java/src/org/broadinstitute/sting/alignment/bwa/BWAAligner.java @@ -112,7 +112,7 @@ public class BWAAligner implements Aligner { List lowerBounds = alignment.negativeStrand ? reverseLowerBounds : forwardLowerBounds; // if z < D(i) then return {} - int mismatches = maxDiff - alignment.mismatches - alignment.gapOpens - alignment.gapExtensions; + int mismatches = maxDiff - alignment.getMismatches() - alignment.getGapOpens() - alignment.getGapExtensions(); if( alignment.position < lowerBounds.size()-1 && mismatches < lowerBounds.get(alignment.position+1).value ) continue; @@ -139,7 +139,7 @@ public class BWAAligner implements Aligner { successfulMatches.add(finalAlignment); bestScore = Math.min(finalAlignment.getScore(),bestScore); - bestDiff = Math.min(finalAlignment.mismatches+finalAlignment.gapOpens+finalAlignment.gapExtensions,bestDiff); + bestDiff = Math.min(finalAlignment.getMismatches()+finalAlignment.getGapOpens()+finalAlignment.getGapExtensions(),bestDiff); maxDiff = bestDiff + 1; } @@ -172,36 +172,36 @@ public class BWAAligner implements Aligner { } if( allowDifferences && - alignment.position+1 >= INDEL_END_SKIP-1+alignment.gapOpens+alignment.gapExtensions && - read.getReadLength()-1-(alignment.position+1) >= INDEL_END_SKIP+alignment.gapOpens+alignment.gapExtensions ) { + alignment.position+1 >= INDEL_END_SKIP-1+alignment.getGapOpens()+alignment.getGapExtensions() && + read.getReadLength()-1-(alignment.position+1) >= INDEL_END_SKIP+alignment.getGapOpens()+alignment.getGapExtensions() ) { if( alignment.getCurrentState() == AlignmentState.MATCH_MISMATCH ) { - if( alignment.gapOpens < MAXIMUM_GAP_OPENS ) { + if( alignment.getGapOpens() < MAXIMUM_GAP_OPENS ) { // Add a potential insertion extension. BWAAlignment insertionAlignment = createInsertionAlignment(alignment); - insertionAlignment.gapOpens++; + insertionAlignment.incrementGapOpens(); alignments.add(insertionAlignment); // Add a potential deletion by marking a deletion and augmenting the position. List deletionAlignments = createDeletionAlignments(bwt,alignment); for( BWAAlignment deletionAlignment: deletionAlignments ) - deletionAlignment.gapOpens++; + deletionAlignment.incrementGapOpens(); alignments.addAll(deletionAlignments); } } else if( alignment.getCurrentState() == AlignmentState.INSERTION ) { - if( alignment.gapExtensions < MAXIMUM_GAP_EXTENSIONS && mismatches > 0 ) { + if( alignment.getGapExtensions() < MAXIMUM_GAP_EXTENSIONS && mismatches > 0 ) { // Add a potential insertion extension. BWAAlignment insertionAlignment = createInsertionAlignment(alignment); - insertionAlignment.gapExtensions++; + insertionAlignment.incrementGapExtensions(); alignments.add(insertionAlignment); } } else if( alignment.getCurrentState() == AlignmentState.DELETION ) { - if( alignment.gapExtensions < MAXIMUM_GAP_EXTENSIONS && mismatches > 0 ) { + if( alignment.getGapExtensions() < MAXIMUM_GAP_EXTENSIONS && mismatches > 0 ) { // Add a potential deletion by marking a deletion and augmenting the position. List deletionAlignments = createDeletionAlignments(bwt,alignment); for( BWAAlignment deletionAlignment: deletionAlignments ) - deletionAlignment.gapExtensions++; + deletionAlignment.incrementGapExtensions(); alignments.addAll(deletionAlignments); } } @@ -225,7 +225,6 @@ public class BWAAligner implements Aligner { seed.position = -1; seed.loBound = 0; seed.hiBound = bwt.length(); - seed.mismatches = 0; return seed; } @@ -271,7 +270,7 @@ public class BWAAligner implements Aligner { newAlignment.position++; newAlignment.addState(AlignmentState.MATCH_MISMATCH); if( Bases.fromASCII(bases[newAlignment.position]) == null || base != Bases.fromASCII(bases[newAlignment.position]) ) - newAlignment.mismatches++; + newAlignment.incrementMismatches(); newAlignments.add(newAlignment); } diff --git a/java/src/org/broadinstitute/sting/alignment/bwa/BWAAlignment.java b/java/src/org/broadinstitute/sting/alignment/bwa/BWAAlignment.java index 6cc42ec68..0b49660f3 100644 --- a/java/src/org/broadinstitute/sting/alignment/bwa/BWAAlignment.java +++ b/java/src/org/broadinstitute/sting/alignment/bwa/BWAAlignment.java @@ -50,17 +50,17 @@ public class BWAAlignment implements Alignment, Cloneable { /** * Working variable. How many mismatches have been encountered at this point. */ - protected int mismatches; + private int mismatches; /** * Number of gap opens in alignment. */ - protected int gapOpens; + private int gapOpens; /** * Number of gap extensions in alignment. */ - protected int gapExtensions; + private int gapExtensions; /** * Working variable. The lower bound of the alignment within the BWT. @@ -72,6 +72,11 @@ public class BWAAlignment implements Alignment, Cloneable { */ protected int hiBound; + /** + * Cache the score. + */ + private int score; + /** * Gets the starting position for the given alignment. * @return Starting position. @@ -113,13 +118,35 @@ public class BWAAlignment implements Alignment, Cloneable { * @return BWA-style scores. 0 is best. */ public int getScore() { - return mismatches*aligner.MISMATCH_PENALTY + gapOpens*aligner.GAP_OPEN_PENALTY + gapExtensions*aligner.GAP_EXTENSION_PENALTY; + return score; } public int getMismatches() { return mismatches; } public int getGapOpens() { return gapOpens; } public int getGapExtensions() { return gapExtensions; } + public void incrementMismatches() { + this.mismatches++; + updateScore(); + } + + public void incrementGapOpens() { + this.gapOpens++; + updateScore(); + } + + public void incrementGapExtensions() { + this.gapExtensions++; + updateScore(); + } + + /** + * Updates the score based on new information about matches / mismatches. + */ + private void updateScore() { + score = mismatches*aligner.MISMATCH_PENALTY + gapOpens*aligner.GAP_OPEN_PENALTY + gapExtensions*aligner.GAP_EXTENSION_PENALTY; + } + /** * Create a new alignment with the given parent aligner. * @param aligner Aligner being used. @@ -164,12 +191,15 @@ public class BWAAlignment implements Alignment, Cloneable { public int compareTo(Alignment rhs) { BWAAlignment other = (BWAAlignment)rhs; - // If the scores are equal, use the score to disambiguate. - int scoreComparison = Integer.valueOf(getScore()).compareTo(other.getScore()); - if( scoreComparison != 0 ) - return scoreComparison; + // If the scores are different, disambiguate using the score. + if(score != other.score) + return score > other.score ? 1 : -1; - return -Long.valueOf(this.creationNumber).compareTo(other.creationNumber); + // Otherwise, use the order in which the elements were created. + if(this.creationNumber != other.creationNumber) + return this.creationNumber > other.creationNumber ? -1 : 1; + + return 0; } public String toString() { diff --git a/java/src/org/broadinstitute/sting/alignment/bwa/bwt/Counts.java b/java/src/org/broadinstitute/sting/alignment/bwa/bwt/Counts.java index ece0cfb9e..48a83ffa8 100644 --- a/java/src/org/broadinstitute/sting/alignment/bwa/bwt/Counts.java +++ b/java/src/org/broadinstitute/sting/alignment/bwa/bwt/Counts.java @@ -13,10 +13,15 @@ import java.util.Map; */ public class Counts implements Cloneable { /** - * Internal representation of counts, broken down by pack value. + * Internal representation of counts, broken down by ASCII value. */ private Map counts = new HashMap(); + /** + * Internal representation of cumulative counts, broken down by ASCII value. + */ + private Map cumulativeCounts = new HashMap(); + /** * Create an empty Counts object with values A=0,C=0,G=0,T=0. */ @@ -28,16 +33,22 @@ public class Counts implements Cloneable { * @param cumulative Whether the counts are cumulative, (count_G=numA+numC+numG,for example). */ public Counts( int[] data, boolean cumulative ) { - for( byte base: Bases.instance) - counts.put(base,data[Bases.toPack(base)]); - - // De-cumulatize data as necessary. if(cumulative) { - int previousCount = 0; - for( byte base: Bases.instance ) { - int count = counts.get(base); - counts.put(base,count-previousCount); - previousCount = count; + int priorCount = 0; + for(byte base: Bases.instance) { + int count = data[Bases.toPack(base)]; + counts.put(base,count-priorCount); + cumulativeCounts.put(base,priorCount); + priorCount = count; + } + } + else { + int priorCount = 0; + for(byte base: Bases.instance) { + int count = data[Bases.toPack(base)]; + counts.put(base,count); + cumulativeCounts.put(base,priorCount); + priorCount += count; } } } @@ -49,11 +60,18 @@ public class Counts implements Cloneable { */ public int[] toArray(boolean cumulative) { int[] countArray = new int[counts.size()]; - for(byte base: Bases.instance) - countArray[Bases.toPack(base)] = counts.get(base); if(cumulative) { - for( int i = 1; i < countArray.length; i++ ) - countArray[i] += countArray[i-1]; + int index = 0; + for(byte base: Bases.instance) { + if(index++ == 0) continue; + countArray[index] = getCumulative(base); + } + countArray[countArray.length-1] = getTotal(); + } + else { + int index = 0; + for(byte base: Bases.instance) + countArray[index] = get(base); } return countArray; } @@ -101,12 +119,7 @@ public class Counts implements Cloneable { * @return Number of bases of this type seen. */ public int getCumulative(byte base) { - int accum = 0; - for( byte current: Bases.allOf() ) { - if(base == current) break; - accum += counts.get(current); - } - return accum; + return cumulativeCounts.get(base); } /** @@ -115,8 +128,9 @@ public class Counts implements Cloneable { */ public int getTotal() { int accumulator = 0; - for( int count : counts.values() ) - accumulator += count; + for(byte base: Bases.instance) { + accumulator += get(base); + } return accumulator; } }