diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/indels/IndelRealigner.java b/java/src/org/broadinstitute/sting/gatk/walkers/indels/IndelRealigner.java index 8561711cb..51e39d768 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/indels/IndelRealigner.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/indels/IndelRealigner.java @@ -108,6 +108,8 @@ public class IndelRealigner extends ReadWalker { private FileWriter statsOutput = null; private FileWriter snpsOutput = null; +// private boolean DEBUG = true; + public void initialize() { if ( LOD_THRESHOLD < 0.0 ) @@ -409,11 +411,17 @@ public class IndelRealigner extends ReadWalker { final AlignedRead aRead = new AlignedRead(read); + // if ( DEBUG ) System.out.println("\nREAD "+read.getReadName()); + // first, move existing indels (for 1 indel reads only) to leftmost position within identical sequence int numBlocks = AlignmentUtils.getNumAlignmentBlocks(read); if ( numBlocks == 2 ) { +// if ( DEBUG ) System.out.println("Attempting to left-shift: read "+read.getReadName()+" alignment: "+read.getCigar().toString()+ +// " at "+read.getAlignmentStart()+"; ref chunk starts at "+leftmostIndex); + Cigar newCigar = indelRealignment(read.getCigar(), reference, read.getReadBases(), read.getAlignmentStart()-(int)leftmostIndex, 0); if ( aRead.setCigar(newCigar) ) { +// if ( DEBUG ) System.out.println("CIGAR CHANGED to "+newCigar.toString()) ; leftMovedIndels.add(aRead); } } @@ -431,7 +439,10 @@ public class IndelRealigner extends ReadWalker { if ( numBlocks == 2 ) { Consensus c = createAlternateConsensus(aRead.getAlignmentStart() - (int)leftmostIndex, aRead.getCigar(), reference, aRead.getRead().getReadBases()); if ( c == null ) {} //System.out.println("ERROR: Failed to create alt consensus for read "+aRead.getRead().getReadName()); - else altConsenses.add(c); + else { +// if ( DEBUG ) System.out.println("Adding consensus: read start "+aRead.getAlignmentStart()+" ref offset "+leftmostIndex+" read cigar "+aRead.getCigar()); + altConsenses.add(c); + } } else { // if ( debugOn ) System.out.println("Going to test..."); @@ -453,6 +464,7 @@ public class IndelRealigner extends ReadWalker { Consensus c = createAlternateConsensus(swConsensus.getAlignmentStart2wrt1(), swConsensus.getCigar(), reference, aRead.getRead().getReadBases()); if ( c != null ) { // if ( debugOn ) System.out.println("NEW consensus generated by SW: "+c.str ) ; + // if ( DEBUG ) System.out.println("Found SW consensus: ref starts at "+leftmostIndex+" alignment at "+swConsensus.getAlignmentStart2wrt1()+" cigar "+swConsensus.getCigar().toString()); altConsenses.add(c); } else { // if ( debugOn ) System.out.println("FAILED to create Alt consensus from SW"); @@ -480,6 +492,16 @@ public class IndelRealigner extends ReadWalker { while ( iter.hasNext() ) { Consensus consensus = iter.next(); +// if ( DEBUG ) { +// System.out.println("Checking consensus with alignment at "+consensus.positionOnReference+" cigar "+consensus.cigar); +// System.out.println(new String(consensus.str)); +// int z = 0; +// for ( ; z < consensus.positionOnReference; z++ ) System.out.print('.'); +// for ( z=0 ; z < consensus.cigar.getCigarElement(0).getLength() ; z++ ) System.out.print('.'); +// if ( consensus.cigar.getCigarElement(1).getOperator() == CigarOperator.I ) for ( z= 0; z < consensus.cigar.getCigarElement(1).getLength(); z++ ) System.out.print('I'); +// System.out.println(); +// } + // if ( debugOn ) System.out.println("Consensus: "+consensus.str); for ( int j = 0; j < altReads.size(); j++ ) { @@ -488,6 +510,17 @@ public class IndelRealigner extends ReadWalker { // the mismatch score is the min of its alignment vs. the reference and vs. the alternate int myScore = altAlignment.second; + +// if ( DEBUG ) { +// if ( toTest.getRead().getReadName().equals("302U3AAXX090430:1:35:429:1940#0") ) { +// System.out.println("READ: new score="+myScore+"; old score="+toTest.getMismatchScoreToReference() ); +// } +// if ( toTest.getRead().getReadName().equals("426RRAAXX090524:8:42:1730:53#0") ) { +// System.out.println("READ: new score="+myScore+"; old score="+toTest.getMismatchScoreToReference() +" start at "+toTest.getAlignmentStart()); +// } +// System.out.println(toTest.getRead().getReadName()+" old score="+toTest.getMismatchScoreToReference()+" new score="+myScore ); +// +// } if ( myScore >= toTest.getMismatchScoreToReference() ) myScore = toTest.getMismatchScoreToReference(); // keep track of reads that align better to the alternate consensus. @@ -505,12 +538,15 @@ public class IndelRealigner extends ReadWalker { break; } +// if ( DEBUG ) System.out.println("Consensus score: "+consensus.mismatchSum); + //logger.debug(consensus.str + " " + consensus.mismatchSum); if ( bestConsensus == null || bestConsensus.mismatchSum > consensus.mismatchSum) { // we do not need this alt consensus, release memory right away!! if ( bestConsensus != null ) bestConsensus.readIndexes.clear(); bestConsensus = consensus; + // if ( DEBUG ) System.out.println("Current consensus is better!"); //logger.debug(consensus.str + " " + consensus.mismatchSum); } else { // we do not need this alt consensus, release memory right away!! @@ -957,6 +993,7 @@ public class IndelRealigner extends ReadWalker { int indelIndexOnRead = readIndex+ce1.getLength(); // position of the indel on the READ (first insterted base, of first base after deletion) byte[] indelString = new byte[ce2.getLength()]; // inserted or deleted sequence + if ( ce2.getOperator() == CigarOperator.D ) System.arraycopy(refSeq, indelIndexOnRef, indelString, 0, ce2.getLength()); else if ( ce2.getOperator() == CigarOperator.I ) @@ -966,6 +1003,8 @@ public class IndelRealigner extends ReadWalker { // for now, we'll just punt the issue and not try to realign these return cigar; + // if ( DEBUG ) System.out.println("Indel sequence: "+new String(indelString)); + // now we have to check all WHOLE periods of the indel sequence: // for instance, if // REF: AGCTATATATAGCC @@ -985,11 +1024,13 @@ public class IndelRealigner extends ReadWalker { // position, we will move insertion left, to the position right after CA. This way, while moving the indel across the repeat // on the ref, we can theoretically move it across a non-repeat on the read if the latter has a mismtach. + // if ( DEBUG ) System.out.println("Starting with "+indelIndexOnRef+" on the ref"); while ( period < indel_length ) { // we will always get at least trivial period = indelStringLength period = BaseUtils.sequencePeriod(indelString, period+1); if ( indel_length % period != 0 ) continue; // if indel sequence length is not a multiple of the period, it's not gonna work +// if ( DEBUG ) System.out.println("Checking period of "+period); int newIndex = indelIndexOnRef; @@ -1007,9 +1048,14 @@ public class IndelRealigner extends ReadWalker { break; } } - if ( match ) + if ( match ) { +// if ( DEBUG ) System.out.println("Success! Moving left by "+period); newIndex -= period; // yes, they are the same, we can move indel farther left by at least period bases, go check if we can do more... - else break; // oops, no match, can not push indel farther left + } + else { +// if ( DEBUG ) System.out.println("Can not push further left with period "+period+"; adjusted position: "+newIndex); + break; // oops, no match, can not push indel farther left + } } final int newDifference = indelIndexOnRef - newIndex; @@ -1025,6 +1071,8 @@ public class IndelRealigner extends ReadWalker { if ( difference > 0 ) { + // if ( DEBUG ) System.out.println("Best offset found: "+difference); + // The following if() statement: this should've never happened, unless the alignment is really screwed up. // A real life example: // @@ -1181,12 +1229,18 @@ public class IndelRealigner extends ReadWalker { readIndexes = new ArrayList>(); } + @Override public boolean equals(Object o) { - return ( this == o || (o instanceof Consensus && this.str.equals(((Consensus)o).str)) ); + return ( this == o || (o instanceof Consensus && Arrays.equals(this.str,(((Consensus)o).str)) ) ); } public boolean equals(Consensus c) { - return ( this == c || this.str.equals(c.str) ); + return ( this == c || Arrays.equals(this.str,c.str) ) ; + } + + @Override + public int hashCode() { + return Arrays.hashCode(this.str); } }