diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReads.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReads.java index d2b9803ef..93d8319b7 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReads.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReads.java @@ -25,9 +25,6 @@ package org.broadinstitute.sting.gatk.walkers.compression.reducereads; -import net.sf.samtools.Cigar; -import net.sf.samtools.CigarElement; -import net.sf.samtools.CigarOperator; import net.sf.samtools.util.SequenceUtil; import org.broadinstitute.sting.commandline.Argument; import org.broadinstitute.sting.commandline.Hidden; @@ -183,7 +180,7 @@ public class ReduceReads extends ReadWalker, ReduceRea * A value of 0 turns downsampling off. */ @Argument(fullName = "downsample_coverage", shortName = "ds", doc = "", required = false) - protected int downsampleCoverage = 0; + protected int downsampleCoverage = 250; @Hidden @Argument(fullName = "", shortName = "dl", doc = "", required = false) @@ -535,81 +532,12 @@ public class ReduceReads extends ReadWalker, ReduceRea if (debugLevel == 1) System.out.println("BAM: " + read.getCigar() + " " + read.getAlignmentStart() + " " + read.getAlignmentEnd()); -// if (!DONT_USE_SOFTCLIPPED_BASES) -// reSoftClipBases(read); - if (!DONT_COMPRESS_READ_NAMES) compressReadName(read); out.addAlignment(read); } - private void reSoftClipBases(GATKSAMRecord read) { - Integer left = (Integer) read.getTemporaryAttribute("SL"); - Integer right = (Integer) read.getTemporaryAttribute("SR"); - if (left != null || right != null) { - Cigar newCigar = new Cigar(); - for (CigarElement element : read.getCigar().getCigarElements()) { - newCigar.add(new CigarElement(element.getLength(), element.getOperator())); - } - - if (left != null) { - newCigar = updateFirstSoftClipCigarElement(left, newCigar); - read.setAlignmentStart(read.getAlignmentStart() + left); - } - - if (right != null) { - Cigar invertedCigar = invertCigar(newCigar); - newCigar = invertCigar(updateFirstSoftClipCigarElement(right, invertedCigar)); - } - read.setCigar(newCigar); - } - } - - /** - * Facility routine to revert the first element of a Cigar string (skipping hard clips) into a soft-clip. - * To be used on both ends if provided a flipped Cigar - * - * @param softClipSize the length of the soft clipped element to add - * @param originalCigar the original Cigar string - * @return a new Cigar object with the soft clips added - */ - private Cigar updateFirstSoftClipCigarElement (int softClipSize, Cigar originalCigar) { - Cigar result = new Cigar(); - CigarElement leftElement = new CigarElement(softClipSize, CigarOperator.S); - boolean updated = false; - for (CigarElement element : originalCigar.getCigarElements()) { - if (!updated && element.getOperator() == CigarOperator.M) { - result.add(leftElement); - int newLength = element.getLength() - softClipSize; - if (newLength > 0) - result.add(new CigarElement(newLength, CigarOperator.M)); - updated = true; - } - else - result.add(element); - } - return result; - } - - /** - * Given a cigar string, returns the inverted cigar string. - * - * @param cigar the original cigar - * @return the inverted cigar - */ - private Cigar invertCigar(Cigar cigar) { - Stack stack = new Stack(); - for (CigarElement e : cigar.getCigarElements()) - stack.push(e); - Cigar inverted = new Cigar(); - while (!stack.empty()) { - inverted.add(stack.pop()); - } - return inverted; - } - - /** * Quality control procedure that checks if the consensus reads contains too many * mismatches with the reference. This should never happen and is a good trigger for diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SlidingWindow.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SlidingWindow.java index 60d7096a3..6a77b0c65 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SlidingWindow.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SlidingWindow.java @@ -499,7 +499,7 @@ public class SlidingWindow { result.addAll(addToSyntheticReads(0, start)); result.addAll(finalizeAndAdd(ConsensusType.BOTH)); - for (GATKSAMRecord read : result) { + for (GATKSAMRecord read : allReads) { readsInWindow.remove(read); // todo -- not optimal, but needs to be done so the next region doesn't try to remove the same reads from the header counts. } @@ -627,7 +627,7 @@ public class SlidingWindow { int locationIndex = startLocation < 0 ? 0 : readStart - startLocation; if (removeRead && locationIndex < 0) - throw new ReviewedStingException("read is behind the Sliding Window. read: " + read + " cigar: " + read.getCigarString() + " window: " + startLocation + "," + stopLocation); + throw new ReviewedStingException("read is behind the Sliding Window. read: " + read + " start " + read.getUnclippedStart() + "," + read.getUnclippedEnd() + " cigar: " + read.getCigarString() + " window: " + startLocation + "," + stopLocation); if (!removeRead) { // we only need to create new header elements if we are adding the read, not when we're removing it if (locationIndex < 0) { // Do we need to add extra elements before the start of the header? -- this may happen if the previous read was clipped and this alignment starts before the beginning of the window diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsIntegrationTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsIntegrationTest.java index 08f7ddd37..ee894b420 100755 --- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsIntegrationTest.java +++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReadsIntegrationTest.java @@ -21,28 +21,28 @@ public class ReduceReadsIntegrationTest extends WalkerTest { @Test(enabled = true) public void testDefaultCompression() { - RRTest("testDefaultCompression ", L, "323dd4deabd7767efa0f2c6e7fa4189f"); + RRTest("testDefaultCompression ", L, "72eb6db9d7a09a0cc25eaac1aafa97b7"); } @Test(enabled = true) public void testMultipleIntervals() { String intervals = "-L 20:10,100,000-10,100,500 -L 20:10,200,000-10,200,500 -L 20:10,300,000-10,300,500 -L 20:10,400,000-10,500,000 -L 20:10,500,050-10,500,060 -L 20:10,600,000-10,600,015 -L 20:10,700,000-10,700,110"; - RRTest("testMultipleIntervals ", intervals, "c437fb160547ff271f8eba30e5f3ff76"); + RRTest("testMultipleIntervals ", intervals, "104b1a1d9fa5394c6fea95cd32967b78"); } @Test(enabled = true) public void testHighCompression() { - RRTest("testHighCompression ", " -cs 10 -minvar 0.3 -mindel 0.3 " + L, "3a607bc3ebaf84e9dc44e005c5f8a047"); + RRTest("testHighCompression ", " -cs 10 -minvar 0.3 -mindel 0.3 " + L, "c55140cec60fa8c35161680289d74d47"); } @Test(enabled = true) public void testLowCompression() { - RRTest("testLowCompression ", " -cs 30 -minvar 0.01 -mindel 0.01 -minmap 5 -minqual 5 " + L, "afd39459c841b68a442abdd5ef5f8f27"); + RRTest("testLowCompression ", " -cs 30 -minvar 0.01 -mindel 0.01 -minmap 5 -minqual 5 " + L, "0f2e57b7f6de03cc4da1ffcc8cf8f1a7"); } @Test(enabled = true) public void testIndelCompression() { - RRTest("testIndelCompression ", " -cs 50 -L 20:10,100,500-10,100,600 ", "f7b9fa44c10bc4b2247813d2b8dc1973"); + RRTest("testIndelCompression ", " -cs 50 -L 20:10,100,500-10,100,600 ", "dda0c95f56f90e5f633c2437c2b21031"); } @Test(enabled = true)