diff --git a/public/gatk-utils/src/main/java/org/broadinstitute/gatk/utils/clipping/ClippingOp.java b/public/gatk-utils/src/main/java/org/broadinstitute/gatk/utils/clipping/ClippingOp.java index 26147e2e9..1873bc1cb 100644 --- a/public/gatk-utils/src/main/java/org/broadinstitute/gatk/utils/clipping/ClippingOp.java +++ b/public/gatk-utils/src/main/java/org/broadinstitute/gatk/utils/clipping/ClippingOp.java @@ -26,6 +26,7 @@ package org.broadinstitute.gatk.utils.clipping; import com.google.java.contract.Requires; +import com.google.common.annotations.VisibleForTesting; import htsjdk.samtools.Cigar; import htsjdk.samtools.CigarElement; import htsjdk.samtools.CigarOperator; @@ -196,35 +197,43 @@ public class ClippingOp { } /** - * Given a cigar string, get the number of bases hard or soft clipped at the start + * Given two cigar strings corresponding to read before and after soft-clipping, returns an integer + * corresponding to the number of reference bases that the new string must be offset by in order to have the + * correct start according to the reference. + * + * @param clippedCigar the new cigar string after clipping + * @param oldCigar the cigar string of the read before it was soft clipped */ - private int getNewAlignmentStartOffset(final Cigar __cigar, final Cigar __oldCigar) { - int num = 0; - for (CigarElement e : __cigar.getCigarElements()) { + @VisibleForTesting + static int getNewAlignmentStartOffset(final Cigar clippedCigar, final Cigar oldCigar) { + int readBasesBeforeReference = 0; // The number of read bases consumed on the new cigar before reference bases are consumed + + int basesBeforeReferenceOld = 0; // The number of read bases consumed on the old cigar before reference bases are consumed + int curReadCounter = 0; // A measure of the reference offset between the oldCigar and the clippedCigar + + for (final CigarElement e : clippedCigar.getCigarElements()) { if (!e.getOperator().consumesReferenceBases()) { if (e.getOperator().consumesReadBases()) { - num += e.getLength(); + readBasesBeforeReference += e.getLength(); } } else { - break; + if (!e.getOperator().consumesReadBases()) { + basesBeforeReferenceOld -= e.getLength(); // Accounting for any D or N cigar operators at the front of the string + } else { + break; + } } } - int oldNum = 0; - int curReadCounter = 0; - for (CigarElement e : __oldCigar.getCigarElements()) { + for (final CigarElement e : oldCigar.getCigarElements()) { int curRefLength = e.getLength(); - int curReadLength = e.getLength(); - if (!e.getOperator().consumesReadBases()) { - curReadLength = 0; - } + int curReadLength = e.getOperator().consumesReadBases() ? e.getLength() : 0; - boolean truncated = false; - if (curReadCounter + curReadLength > num) { - curReadLength = num - curReadCounter; - curRefLength = num - curReadCounter; - truncated = true; + final boolean truncated = curReadCounter + curReadLength > readBasesBeforeReference; + if (truncated) { + curReadLength = readBasesBeforeReference - curReadCounter; + curRefLength = readBasesBeforeReference - curReadCounter; } if (!e.getOperator().consumesReferenceBases()) { @@ -232,14 +241,14 @@ public class ClippingOp { } curReadCounter += curReadLength; - oldNum += curRefLength; + basesBeforeReferenceOld += curRefLength; - if (curReadCounter > num || truncated) { + if (curReadCounter > readBasesBeforeReference || truncated) { break; } } - return oldNum; + return Math.abs(basesBeforeReferenceOld); // if oldNum is negative it means some of the preceding N/Ds were trimmed but not all so we take absolute value } /** diff --git a/public/gatk-utils/src/test/java/org/broadinstitute/gatk/utils/clipping/ClippingOpUnitTest.java b/public/gatk-utils/src/test/java/org/broadinstitute/gatk/utils/clipping/ClippingOpUnitTest.java new file mode 100644 index 000000000..f0610429d --- /dev/null +++ b/public/gatk-utils/src/test/java/org/broadinstitute/gatk/utils/clipping/ClippingOpUnitTest.java @@ -0,0 +1,81 @@ +/* +* Copyright 2012-2016 Broad Institute, Inc. +* +* Permission is hereby granted, free of charge, to any person +* obtaining a copy of this software and associated documentation +* files (the "Software"), to deal in the Software without +* restriction, including without limitation the rights to use, +* copy, modify, merge, publish, distribute, sublicense, and/or sell +* copies of the Software, and to permit persons to whom the +* Software is furnished to do so, subject to the following +* conditions: +* +* The above copyright notice and this permission notice shall be +* included in all copies or substantial portions of the Software. +* +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR +* THE USE OR OTHER DEALINGS IN THE SOFTWARE. +*/ + +package org.broadinstitute.gatk.utils.clipping; + +import htsjdk.samtools.Cigar; +import htsjdk.samtools.TextCigarCodec; +import org.broadinstitute.gatk.utils.BaseTest; +import org.testng.Assert; +import org.testng.annotations.DataProvider; +import org.testng.annotations.Test; + +public final class ClippingOpUnitTest extends BaseTest { + + @Test (dataProvider = "SoftClipedReadsNewStart") + public void testGetNewAlignmentStartOffset(final String preClip, final String postClip, final int expectedResult) { + Cigar cpreClip = TextCigarCodec.decode(preClip); + Cigar cpostClip = TextCigarCodec.decode(postClip); + Assert.assertEquals(ClippingOp.getNewAlignmentStartOffset(cpostClip, cpreClip), expectedResult, + "getNewAlignmentStartOffset returned "+ClippingOp.getNewAlignmentStartOffset(cpostClip, cpreClip)+ + " when "+expectedResult+" was expected for "+preClip.toString()+" which was clipped to "+postClip.toString()); + } + + // provider fields: cigar string before clip, cigar string after clip, expected result for getNewAlignmentStartOffset() method + @DataProvider(name = "SoftClipedReadsNewStart") + public Object[][] makeRevertSoftClipsBeforeContig() { + return new Object[][] { + {"70M", "10S60M", 10}, + {"70M", "60M10S", 0}, + + {"30M10N30M", "30S10N30M", 30}, + {"30M10N30M", "30S5N30M", 35}, + {"30M10N30M", "30M10N30S", 0}, + {"30M10N30M", "30S30M", 40}, + {"30M10N30M", "30M30S", 0}, + {"30M10N30M", "15S15M10N30M", 15}, + {"30M10N30M", "30M10N15M15S", 0}, + // Testing multiple sequential reference but not sequence consuming base types + {"10N10D40M", "40M", 20}, + {"10N10D40M", "20S20M", 40}, + {"10N10D40M", "5N10D40M", 5}, + {"10N10D40M", "5D40M", 15}, + {"10N10D40M", "10N10D20M20S", 0}, + + {"10S10N20M10N10M", "10S20M10N10M", 10}, + {"10S10N20M10N10M", "10S10N20M10S", 0}, + + {"10S10I20M10I10M", "20S20M10I10M", 0}, + {"10S10I20M10I10M", "10S10I20M20S", 0}, + {"10S10I20M10I10M", "15S5I20M10I10M", 0}, + + {"10H60M", "10H10S50M", 10}, + {"10H60M", "10H50M10S", 0}, + {"10H10S50M", "10H20S40M", 10}, + {"10X60M", "20S50M", 20}, + {"10I40N20M","10S20M",40} + }; + } +}