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 93d8319b7..177050667 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 @@ -251,7 +251,7 @@ public class ReduceReads extends ReadWalker, ReduceRea LinkedList mappedReads; totalReads++; if (!debugRead.isEmpty() && read.getReadName().contains(debugRead)) - System.out.println("Found debug read!"); + System.out.println("Found debug read!"); if (debugLevel == 1) System.out.printf("\nOriginal: %s %s %d %d\n", read, read.getCigar(), read.getAlignmentStart(), read.getAlignmentEnd()); @@ -260,7 +260,14 @@ public class ReduceReads extends ReadWalker, ReduceRea // attribute hash so we can determine later if we need to write down the alignment shift to the reduced BAM file read.setTemporaryAttribute(GATKSAMRecord.REDUCED_READ_ORIGINAL_ALIGNMENT_START_SHIFT, read.getAlignmentStart()); read.setTemporaryAttribute(GATKSAMRecord.REDUCED_READ_ORIGINAL_ALIGNMENT_END_SHIFT, read.getAlignmentEnd()); - + + // Check if the read goes beyond the boundaries of the chromosome, and hard clip those boundaries. + int chromosomeLength = ref.getGenomeLocParser().getContigInfo(read.getReferenceName()).getSequenceLength(); + if (read.getSoftStart() < 0) + read = ReadClipper.hardClipByReadCoordinates(read, 0, -read.getSoftStart() - 1); + if (read.getSoftEnd() > chromosomeLength) + read = ReadClipper.hardClipByReadCoordinates(read, chromosomeLength - read.getSoftStart() + 1, read.getReadLength() - 1); + if (!DONT_SIMPLIFY_READS) read.simplify(); // Clear all unnecessary attributes if (!DONT_CLIP_ADAPTOR_SEQUENCES) 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 6a77b0c65..bdb9ef843 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 @@ -536,6 +536,10 @@ public class SlidingWindow { * @return a list of reads selected by the downsampler to cover the window to at least the desired coverage */ protected List downsampleVariantRegion(final List allReads) { + int nReads = allReads.size(); + if (nReads == 0) + return allReads; + double fraction = 100 / allReads.size(); if (fraction >= 1) return allReads; @@ -545,6 +549,7 @@ public class SlidingWindow { return downsampler.consumeDownsampledItems(); } + /** * Properly closes a Sliding Window, finalizing all consensus and variant * regions that still exist regardless of being able to fulfill the diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SyntheticRead.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SyntheticRead.java index 3b9fc0390..6134101d9 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SyntheticRead.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/SyntheticRead.java @@ -102,7 +102,7 @@ public class SyntheticRead { * @param base the base to add * @param count number of reads with this base */ - @Requires("count < Byte.MAX_VALUE") + @Requires("count <= Byte.MAX_VALUE") public void add(BaseIndex base, byte count, byte qual, byte insQual, byte delQual, double mappingQuality) { counts.add(count); bases.add(base); 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 ee894b420..5d16ba019 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 @@ -11,6 +11,8 @@ public class ReduceReadsIntegrationTest extends WalkerTest { final String DELETION_BAM = validationDataLocation + "filtered_deletion_for_reduce_reads.bam"; final String STASH_BAM = validationDataLocation + "ReduceReadsStashBug.bam"; final String STASH_L = " -L 14:73718184-73718284 -L 14:73718294-73718330 -L 14:73718360-73718556"; + final String DIVIDEBYZERO_BAM = validationDataLocation + "ReduceReadsDivideByZeroBug.bam"; + final String DIVIDEBYZERO_L = " -L " + validationDataLocation + "ReduceReadsDivideByZeroBug.intervals"; final String L = " -L 20:10,100,000-10,120,000 "; private void RRTest(String testName, String args, String md5) { @@ -64,5 +66,16 @@ public class ReduceReadsIntegrationTest extends WalkerTest { String base = String.format("-T ReduceReads %s -npt -R %s -I %s", STASH_L, REF, STASH_BAM) + " -o %s "; executeTest("testAddingReadAfterTailingTheStash", new WalkerTestSpec(base, Arrays.asList("886b43e1f26ff18425814dc7563931c6"))); } + + /** + * Divide by zero bug reported by GdA and users in the forum. Happens when the downsampler goes over a region where all reads get + * filtered out. + */ + @Test(enabled = true) + public void testDivideByZero() { + String base = String.format("-T ReduceReads %s -npt -R %s -I %s", DIVIDEBYZERO_L, REF, DIVIDEBYZERO_BAM) + " -o %s "; + executeTest("testDivideByZero", new WalkerTestSpec(base, Arrays.asList("137505c3efd1e9f8d9209dbdf8419ff9"))); + } + } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/diagnostics/targets/FindCoveredIntervals.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/diagnostics/targets/FindCoveredIntervals.java index 23e4d5ae0..373c8232e 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/diagnostics/targets/FindCoveredIntervals.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/diagnostics/targets/FindCoveredIntervals.java @@ -78,7 +78,7 @@ public class FindCoveredIntervals extends ActiveRegionWalker { public Long reduce(final GenomeLoc value, Long reduce) { if (value != null) { out.println(value.toString()); - return reduce++; + return ++reduce; } else return reduce; }