From cc23b0b8a9b184eee6908a2fcbfdf63318be11b1 Mon Sep 17 00:00:00 2001 From: Matt Hanna Date: Fri, 23 Sep 2011 14:52:31 -0400 Subject: [PATCH 4/6] Fix for recent change modelling unmapped shards: don't invoke optimization to combine mapped and unmapped shards. --- .../sting/gatk/datasources/reads/LowMemoryIntervalSharder.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/LowMemoryIntervalSharder.java b/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/LowMemoryIntervalSharder.java index ba6321121..bf5f33dc3 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/LowMemoryIntervalSharder.java +++ b/public/java/src/org/broadinstitute/sting/gatk/datasources/reads/LowMemoryIntervalSharder.java @@ -59,7 +59,7 @@ public class LowMemoryIntervalSharder implements Iterator { */ public FilePointer next() { FilePointer current = wrappedIterator.next(); - while(wrappedIterator.hasNext() && current.minus(wrappedIterator.peek()) == 0) + while(wrappedIterator.hasNext() && current.isRegionUnmapped == wrappedIterator.peek().isRegionUnmapped && current.minus(wrappedIterator.peek()) == 0) current = current.combine(parser,wrappedIterator.next()); return current; } From 9ea40f2e416727564849b749c0746fb688f6f1be Mon Sep 17 00:00:00 2001 From: Mauricio Carneiro Date: Fri, 23 Sep 2011 16:37:08 -0400 Subject: [PATCH 5/6] Deletions/Insertions in hard clip and bug fixes * Deletions now count as hard clipped bases in order to recover the original alignment start of a clipped read. * Insertions do not count as hard clipped bases for the same reason. * This created a bug in the previous cigar cleaning function. Fixed. --- .../sting/utils/clipreads/ClippingOp.java | 59 ++++++++++--------- 1 file changed, 32 insertions(+), 27 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/utils/clipreads/ClippingOp.java b/public/java/src/org/broadinstitute/sting/utils/clipreads/ClippingOp.java index 47ce8165c..81d00d9d7 100644 --- a/public/java/src/org/broadinstitute/sting/utils/clipreads/ClippingOp.java +++ b/public/java/src/org/broadinstitute/sting/utils/clipreads/ClippingOp.java @@ -398,7 +398,9 @@ public class ClippingOp { for (int i = 1; i <= 2; i++) { int shift = 0; + int totalHardClip = 0; boolean readHasStarted = false; + boolean addedHardClips = false; while(!cigarStack.empty()) { CigarElement cigarElement = cigarStack.pop(); @@ -408,14 +410,33 @@ public class ClippingOp { cigarElement.getOperator() != CigarOperator.DELETION && cigarElement.getOperator() != CigarOperator.HARD_CLIP) readHasStarted = true; + + else if ( !readHasStarted && cigarElement.getOperator() == CigarOperator.HARD_CLIP) + totalHardClip += cigarElement.getLength(); + else if ( !readHasStarted && cigarElement.getOperator() == CigarOperator.INSERTION) shift += cigarElement.getLength(); - if (readHasStarted || cigarElement.getOperator() == CigarOperator.HARD_CLIP) { - if (i==1) + else if ( !readHasStarted && cigarElement.getOperator() == CigarOperator.DELETION) + totalHardClip += cigarElement.getLength(); + + if (readHasStarted) { + if (i==1) { + if (!addedHardClips) { + if (totalHardClip > 0) + inverseCigarStack.push(new CigarElement(totalHardClip, CigarOperator.HARD_CLIP)); + addedHardClips = true; + } inverseCigarStack.push(cigarElement); - else + } + else { + if (!addedHardClips) { + if (totalHardClip > 0) + cleanCigar.add(new CigarElement(totalHardClip, CigarOperator.HARD_CLIP)); + addedHardClips = true; + } cleanCigar.add(cigarElement); + } } } // first pass (i=1) is from end to start of the cigar elements @@ -434,7 +455,6 @@ public class ClippingOp { private int calculateAlignmentStartShift(Cigar oldCigar, Cigar newCigar) { int newShift = 0; int oldShift = 0; - int deletionShift = 0; for (CigarElement cigarElement : newCigar.getCigarElements()) { if (cigarElement.getOperator() == CigarOperator.HARD_CLIP || cigarElement.getOperator() == CigarOperator.SOFT_CLIP) @@ -449,34 +469,19 @@ public class ClippingOp { else break; } - - int basesClipped = 0; - for (CigarElement cigarElement : oldCigar.getCigarElements()) { - if (basesClipped > newShift) // are we beyond the clipped region? - break; - - else if (cigarElement.getOperator() == CigarOperator.DELETION) // if this is a deletion, we have to adjust the starting shift - deletionShift += cigarElement.getLength(); - - else - basesClipped += cigarElement.getLength(); - } - - return newShift - oldShift + deletionShift; + return newShift - oldShift; } private int calculateHardClippingAlignmentShift(CigarElement cigarElement, int clippedLength) { - if (cigarElement.getOperator() == CigarOperator.INSERTION) { - int cigarElementLength = cigarElement.getLength(); - if (clippedLength >= cigarElementLength) - return -cigarElement.getLength(); - else - return -clippedLength; - } + // Insertions should be discounted from the total hard clip count + if (cigarElement.getOperator() == CigarOperator.INSERTION) + return -clippedLength; -// if (cigarElement.getOperator() == CigarOperator.DELETION) -// return cigarElement.getLength(); + // Deletions should be added to the total hard clip count + else if (cigarElement.getOperator() == CigarOperator.DELETION) + return cigarElement.getLength(); + // There is no shift if we are not clipping an indel return 0; }