zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #1535 from broadinstitute/rhl_clipping_op_align_offset_1447 

Fixed ClippingOP.getNewAlignmentStartOffset() to handle N cases properly
This commit is contained in:
Ron Levine 2017-03-06 17:53:12 -05:00 committed by GitHub
parent faeac6e5a4 1e88d849e3
commit b3df32047f
2 changed files with 111 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
public/gatk-utils/src/main/java/org/broadinstitute/gatk/utils/clipping/ClippingOp.java
View File

@ -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
}
/**

81
public/gatk-utils/src/test/java/org/broadinstitute/gatk/utils/clipping/ClippingOpUnitTest.java 100644
View File

@ -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}
};
}
}