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 951ab265c..47ce8165c 100644 --- a/public/java/src/org/broadinstitute/sting/utils/clipreads/ClippingOp.java +++ b/public/java/src/org/broadinstitute/sting/utils/clipreads/ClippingOp.java @@ -432,25 +432,37 @@ public class ClippingOp { } private int calculateAlignmentStartShift(Cigar oldCigar, Cigar newCigar) { - int shift = 0; + int newShift = 0; + int oldShift = 0; + int deletionShift = 0; - // Rewind to previous start (by counting everything that was already clipped in this read) - for (CigarElement cigarElement : oldCigar.getCigarElements()) { - if (!cigarElement.getOperator().consumesReferenceBases()) - shift -= cigarElement.getLength(); - else - break; - } - - // Advance to new start (by counting everything new that has been clipped ) for (CigarElement cigarElement : newCigar.getCigarElements()) { - if (!cigarElement.getOperator().consumesReferenceBases()) - shift += cigarElement.getLength(); + if (cigarElement.getOperator() == CigarOperator.HARD_CLIP || cigarElement.getOperator() == CigarOperator.SOFT_CLIP) + newShift += cigarElement.getLength(); else break; } - return shift; + for (CigarElement cigarElement : oldCigar.getCigarElements()) { + if (cigarElement.getOperator() == CigarOperator.HARD_CLIP || cigarElement.getOperator() == CigarOperator.SOFT_CLIP ) + oldShift += Math.min(cigarElement.getLength(), newShift - oldShift); + 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; } private int calculateHardClippingAlignmentShift(CigarElement cigarElement, int clippedLength) { @@ -462,8 +474,8 @@ public class ClippingOp { return -clippedLength; } - if (cigarElement.getOperator() == CigarOperator.DELETION) - return cigarElement.getLength(); +// if (cigarElement.getOperator() == CigarOperator.DELETION) +// return cigarElement.getLength(); return 0; } 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 9b9ce4fdf..d0e7f8371 100644 --- a/public/java/src/org/broadinstitute/sting/utils/clipreads/ReadClipper.java +++ b/public/java/src/org/broadinstitute/sting/utils/clipreads/ReadClipper.java @@ -71,21 +71,21 @@ public class ReadClipper { private SAMRecord hardClipByReferenceCoordinates(int refStart, int refStop) { int start = (refStart < 0) ? 0 : ReadUtils.getReadCoordinateForReferenceCoordinate(read, refStart); - int stop = (refStop < 0) ? read.getReadLength() - 1: ReadUtils.getReadCoordinateForReferenceCoordinate(read, refStop); + int stop = (refStop < 0) ? read.getReadLength() - 1 : ReadUtils.getReadCoordinateForReferenceCoordinate(read, refStop); - if (start < 0 || stop > read.getReadLength() - 1 + numDeletions(read)) + if (start < 0 || stop > read.getReadLength() - 1) throw new ReviewedStingException("Trying to clip before the start or after the end of a read"); - // TODO add check in the Hardclip function - if ( start > stop ) - stop = ReadUtils.getReadCoordinateForReferenceCoordinate(read, ReadUtils.getRefCoordSoftUnclippedEnd(read)); - + if ( start > stop ) { +// stop = ReadUtils.getReadCoordinateForReferenceCoordinate(read, ReadUtils.getRefCoordSoftUnclippedEnd(read)); + throw new ReviewedStingException("START > STOP -- this should never happen -- call Mauricio!"); + } //This tries to fix the bug where the deletion is counted a read base and as a result, the hardCLipper runs into //an endless loop when hard clipping the cigar string because the read coordinates are not covered by the read - stop -= numDeletions(read); - if ( start > stop ) - start -= numDeletions(read); +// stop -= numDeletions(read); +// if ( start > stop ) +// start -= numDeletions(read); //System.out.println("Clipping start/stop: " + start + "/" + stop); @@ -145,7 +145,7 @@ public class ReadClipper { cutLeft = readIndex + cigarElement.getLength() - 1; } } - else + else if (cigarElement.getOperator() != CigarOperator.HARD_CLIP) rightTail = true; if (cigarElement.getOperator().consumesReadBases()) diff --git a/public/java/test/org/broadinstitute/sting/utils/genotype/vcf/VCFWriterUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/genotype/vcf/VCFWriterUnitTest.java index a8e6593b1..35c6a4993 100644 --- a/public/java/test/org/broadinstitute/sting/utils/genotype/vcf/VCFWriterUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/genotype/vcf/VCFWriterUnitTest.java @@ -105,7 +105,6 @@ public class VCFWriterUnitTest extends BaseTest { public static VCFHeader createFakeHeader(Set metaData, Set additionalColumns) { metaData.add(new VCFHeaderLine(VCFHeaderVersion.VCF4_0.getFormatString(), VCFHeaderVersion.VCF4_0.getVersionString())); metaData.add(new VCFHeaderLine("two", "2")); - additionalColumns.add("FORMAT"); additionalColumns.add("extra1"); additionalColumns.add("extra2"); return new VCFHeader(metaData, additionalColumns); @@ -159,6 +158,6 @@ public class VCFWriterUnitTest extends BaseTest { Assert.assertTrue(additionalColumns.contains(key)); index++; } - Assert.assertEquals(index+1, additionalColumns.size() /* for the header field we don't see */); + Assert.assertEquals(index, additionalColumns.size()); } }