From 489c15b99d29907f1c864cac87f434a5e7294cbf Mon Sep 17 00:00:00 2001 From: Mauricio Carneiro Date: Mon, 15 Aug 2011 01:42:34 -0400 Subject: [PATCH] Fixed indexing issue in coordinate conversion When a read had been previously soft clipped, the UnclippedEnd could not be used directly as Reference Coordinate for clipping , because the read does not go that far. --- .../sting/utils/clipreads/ReadClipper.java | 11 +++++++---- .../org/broadinstitute/sting/utils/sam/ReadUtils.java | 6 +++--- 2 files changed, 10 insertions(+), 7 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 1307f1b6f..00405a98c 100644 --- a/public/java/src/org/broadinstitute/sting/utils/clipreads/ReadClipper.java +++ b/public/java/src/org/broadinstitute/sting/utils/clipreads/ReadClipper.java @@ -3,6 +3,7 @@ package org.broadinstitute.sting.utils.clipreads; import net.sf.samtools.Cigar; import net.sf.samtools.CigarElement; import net.sf.samtools.SAMRecord; +import org.broad.tribble.util.PositionalStream; import org.broadinstitute.sting.utils.sam.ReadUtils; import org.jets3t.service.multi.ThreadedStorageService; @@ -48,8 +49,10 @@ public class ReadClipper { } public SAMRecord hardClipByReferenceCoordinates(int refStart, int refStop) { - int start = ReadUtils.getReadCoordinateForReferenceCoordinate(read, refStart); - int stop = ReadUtils.getReadCoordinateForReferenceCoordinate(read, refStop); + 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); this.addOp(new ClippingOp(start, stop)); return clipRead(ClippingRepresentation.HARDCLIP_BASES); } @@ -60,9 +63,9 @@ public class ReadClipper { } public SAMRecord hardClipBothEndsByReferenceCoordinates(int left, int right) { - this.read = hardClipByReferenceCoordinates(read.getUnclippedStart(), left); + this.read = hardClipByReferenceCoordinates(-1, left); this.ops = null; // reset the operations - return hardClipByReferenceCoordinates(right, read.getUnclippedEnd()); + return hardClipByReferenceCoordinates(right, -1); } /** diff --git a/public/java/src/org/broadinstitute/sting/utils/sam/ReadUtils.java b/public/java/src/org/broadinstitute/sting/utils/sam/ReadUtils.java index 3d803804a..5d39a01a5 100644 --- a/public/java/src/org/broadinstitute/sting/utils/sam/ReadUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/sam/ReadUtils.java @@ -650,13 +650,13 @@ public class ReadUtils { return ReadAndIntervalOverlap.RIGHT_OVERLAP; } - @Requires({"refCoord >= read.getUnclippedStart()", "refCoord <= read.getUnclippedEnd()"}) + @Requires({"refCoord >= read.getAlignmentStart()", "refCoord <= read.getAlignmentEnd()"}) @Ensures({"result >= 0", "result < read.getReadLength()"}) public static int getReadCoordinateForReferenceCoordinate(SAMRecord read, int refCoord) { int readBases = 0; int refBases = 0; - int goal = refCoord - read.getUnclippedStart(); // read coords are 0-based! - boolean goalReached = false; + int goal = refCoord - read.getAlignmentStart(); // read coords are 0-based! + boolean goalReached = refBases == goal; Iterator cigarElementIterator = read.getCigar().getCigarElements().iterator(); while (!goalReached && cigarElementIterator.hasNext()) {