From 993ecb85da7e5ae15f811c551312138d53c0144c Mon Sep 17 00:00:00 2001 From: Mauricio Carneiro Date: Mon, 15 Aug 2011 15:22:54 -0400 Subject: [PATCH] Added Hard Clipping Tail Ends Added functionality to hard clip the low quality tail ends of reads (lowQual <= 2) --- .../sting/utils/clipreads/ReadClipper.java | 29 +++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/utils/clipreads/ReadClipper.java b/public/java/src/org/broadinstitute/sting/utils/clipreads/ReadClipper.java index 85063ad14..62b25ddf3 100644 --- a/public/java/src/org/broadinstitute/sting/utils/clipreads/ReadClipper.java +++ b/public/java/src/org/broadinstitute/sting/utils/clipreads/ReadClipper.java @@ -15,6 +15,7 @@ import java.util.List; */ public class ReadClipper { SAMRecord read; + boolean wasClipped; List ops = null; /** @@ -24,6 +25,7 @@ public class ReadClipper { */ public ReadClipper(final SAMRecord read) { this.read = read; + this.wasClipped = false; } /** @@ -41,7 +43,7 @@ public class ReadClipper { } public boolean wasClipped() { - return ops != null; + return wasClipped; } public SAMRecord getRead() { @@ -52,7 +54,7 @@ public class ReadClipper { int start = (refStart < 0) ? 0 : ReadUtils.getReadCoordinateForReferenceCoordinate(read, refStart); int stop = (refStop < 0) ? read.getReadLength() - 1 : ReadUtils.getReadCoordinateForReferenceCoordinate(read, refStop); - System.out.println("DEBUG -- clipping start/stop: " + start + "/" + stop); + System.out.println("Clipping start/stop: " + start + "/" + stop); this.addOp(new ClippingOp(start, stop)); return clipRead(ClippingRepresentation.HARDCLIP_BASES); } @@ -68,6 +70,28 @@ public class ReadClipper { return hardClipByReferenceCoordinates(right, -1); } + public SAMRecord hardClipLowQualEnds(byte lowQual) { + byte [] quals = read.getBaseQualities(); + int leftClipIndex = 0; + int rightClipIndex = read.getReadLength() - 1; + + // check how far we can clip both sides + while (quals[rightClipIndex] <= lowQual) rightClipIndex--; + while (quals[leftClipIndex] <= lowQual) leftClipIndex++; + + // if the entire read should be clipped, then return an empty read. (--todo: maybe null is better? testing this for now) + if (leftClipIndex > rightClipIndex) + return (new SAMRecord(read.getHeader())); + + if (rightClipIndex < read.getReadLength() - 1) { + this.addOp(new ClippingOp(rightClipIndex + 1, read.getReadLength() - 1)); + } + if (leftClipIndex > 0 ) { + this.addOp(new ClippingOp(0, leftClipIndex - 1)); + } + return this.clipRead(ClippingRepresentation.HARDCLIP_BASES); + } + /** * Return a new read corresponding to this.read that's been clipped according to ops, if any are present. * @@ -83,6 +107,7 @@ public class ReadClipper { for (ClippingOp op : getOps()) { clippedRead = op.apply(algorithm, clippedRead); } + wasClipped = true; return clippedRead; } catch (CloneNotSupportedException e) { throw new RuntimeException(e); // this should never happen