From 6866a41914c808e9b1d63edbf4860349143cbd14 Mon Sep 17 00:00:00 2001 From: Guillermo del Angel Date: Thu, 23 Feb 2012 09:45:47 -0500 Subject: [PATCH] Added functionality in pileups to not only determine whether there's an insertion or deletion following the current position, but to also get the indel length and involved bases - definitely needed for extended event removal, and needed for pool caller indel functionality. --- .../gatk/iterators/LocusIteratorByState.java | 15 ++++++--- .../pileup/AbstractReadBackedPileup.java | 1 + .../pileup/ExtendedEventPileupElement.java | 2 +- .../sting/utils/pileup/PileupElement.java | 31 ++++++++++++++++++- .../utils/pileup/ReadBackedPileupImpl.java | 10 ++++-- 5 files changed, 50 insertions(+), 9 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/gatk/iterators/LocusIteratorByState.java b/public/java/src/org/broadinstitute/sting/gatk/iterators/LocusIteratorByState.java index 6edae3816..b8dd03317 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/iterators/LocusIteratorByState.java +++ b/public/java/src/org/broadinstitute/sting/gatk/iterators/LocusIteratorByState.java @@ -175,8 +175,8 @@ public class LocusIteratorByState extends LocusIterator { return String.format("%s ro=%d go=%d co=%d cec=%d %s", read.getReadName(), readOffset, genomeOffset, cigarOffset, cigarElementCounter, curElement); } - public CigarOperator peekForwardOnGenome() { - return ( cigarElementCounter + 1 > curElement.getLength() && cigarOffset + 1 < nCigarElements ? cigar.getCigarElement(cigarOffset + 1) : curElement ).getOperator(); + public CigarElement peekForwardOnGenome() { + return ( cigarElementCounter + 1 > curElement.getLength() && cigarOffset + 1 < nCigarElements ? cigar.getCigarElement(cigarOffset + 1) : curElement ); } public CigarOperator stepForwardOnGenome() { @@ -462,15 +462,19 @@ public class LocusIteratorByState extends LocusIterator { final SAMRecordState state = iterator.next(); // state object with the read/offset information final GATKSAMRecord read = (GATKSAMRecord) state.getRead(); // the actual read final CigarOperator op = state.getCurrentCigarOperator(); // current cigar operator - final CigarOperator nextOp = state.peekForwardOnGenome(); // next cigar operator + final CigarElement nextElement = state.peekForwardOnGenome(); // next cigar element + final CigarOperator nextOp = nextElement.getOperator(); final int readOffset = state.getReadOffset(); // the base offset on this read + byte[] insertedBases = Arrays.copyOfRange(read.getReadBases(), readOffset + 1, readOffset + 1 + nextElement.getLength()); + int nextElementLength = nextElement.getLength(); if (op == CigarOperator.N) // N's are never added to any pileup continue; if (op == CigarOperator.D) { if (readInfo.includeReadsWithDeletionAtLoci()) { // only add deletions to the pileup if we are authorized to do so - pile.add(new PileupElement(read, readOffset, true, nextOp == CigarOperator.D, nextOp == CigarOperator.I, nextOp == CigarOperator.S || (state.getGenomeOffset() == 0 && read.getSoftStart() != read.getAlignmentStart()))); + pile.add(new PileupElement(read, readOffset, true, nextOp == CigarOperator.D, nextOp == CigarOperator.I, nextOp == CigarOperator.S || (state.getGenomeOffset() == 0 && read.getSoftStart() != read.getAlignmentStart()), + null,nextOp == CigarOperator.D? nextElementLength:-1)); size++; nDeletions++; if (read.getMappingQuality() == 0) @@ -479,7 +483,8 @@ public class LocusIteratorByState extends LocusIterator { } else { if (!filterBaseInRead(read, location.getStart())) { - pile.add(new PileupElement(read, readOffset, false, nextOp == CigarOperator.D, nextOp == CigarOperator.I, nextOp == CigarOperator.S || (state.getGenomeOffset() == 0 && read.getSoftStart() != read.getAlignmentStart()))); + pile.add(new PileupElement(read, readOffset, false, nextOp == CigarOperator.D, nextOp == CigarOperator.I, nextOp == CigarOperator.S || (state.getGenomeOffset() == 0 && read.getSoftStart() != read.getAlignmentStart()), + new String(insertedBases),nextElementLength)); size++; if (read.getMappingQuality() == 0) nMQ0Reads++; diff --git a/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java b/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java index 70ad70f43..7c2a67aba 100644 --- a/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java +++ b/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java @@ -205,6 +205,7 @@ public abstract class AbstractReadBackedPileup createNewPileup(GenomeLoc loc, PileupElementTracker pileupElementTracker); protected abstract PE createNewPileupElement(GATKSAMRecord read, int offset, boolean isDeletion, boolean isBeforeDeletion, boolean isBeforeInsertion, boolean isNextToSoftClip); + protected abstract PE createNewPileupElement(GATKSAMRecord read, int offset, boolean isDeletion, boolean isBeforeDeletion, boolean isBeforeInsertion, boolean isNextToSoftClip, String nextEventBases, int nextEventLength ); // -------------------------------------------------------- // diff --git a/public/java/src/org/broadinstitute/sting/utils/pileup/ExtendedEventPileupElement.java b/public/java/src/org/broadinstitute/sting/utils/pileup/ExtendedEventPileupElement.java index 506442d03..8df0aa0b8 100644 --- a/public/java/src/org/broadinstitute/sting/utils/pileup/ExtendedEventPileupElement.java +++ b/public/java/src/org/broadinstitute/sting/utils/pileup/ExtendedEventPileupElement.java @@ -48,7 +48,7 @@ public class ExtendedEventPileupElement extends PileupElement { public ExtendedEventPileupElement(GATKSAMRecord read, int offset, int eventLength, String eventBases, Type type) { - super(read, offset, type == Type.DELETION, false, false, false); // extended events are slated for removal + super(read, offset, type == Type.DELETION, false, false, false,null,-1); // extended events are slated for removal this.read = read; this.offset = offset; this.eventLength = eventLength; diff --git a/public/java/src/org/broadinstitute/sting/utils/pileup/PileupElement.java b/public/java/src/org/broadinstitute/sting/utils/pileup/PileupElement.java index 9df22700e..022cadbbe 100755 --- a/public/java/src/org/broadinstitute/sting/utils/pileup/PileupElement.java +++ b/public/java/src/org/broadinstitute/sting/utils/pileup/PileupElement.java @@ -27,6 +27,10 @@ public class PileupElement implements Comparable { protected final boolean isBeforeDeletion; protected final boolean isBeforeInsertion; protected final boolean isNextToSoftClip; + protected final int eventLength; + protected final String eventBases; // if it is a deletion, we do not have information about the actual deleted bases + // in the read itself, so we fill the string with D's; for insertions we keep actual inserted bases + /** * Creates a new pileup element. @@ -37,12 +41,15 @@ public class PileupElement implements Comparable { * @param isBeforeDeletion whether or not this base is before a deletion * @param isBeforeInsertion whether or not this base is before an insertion * @param isNextToSoftClip whether or not this base is next to a soft clipped base + * @param nextEventBases bases in event in case element comes before insertion or deletion + * @param nextEventLength length of next event in case it's insertion or deletion */ @Requires({ "read != null", "offset >= -1", "offset <= read.getReadLength()"}) - public PileupElement(final GATKSAMRecord read, final int offset, final boolean isDeletion, final boolean isBeforeDeletion, final boolean isBeforeInsertion, final boolean isNextToSoftClip) { + public PileupElement(final GATKSAMRecord read, final int offset, final boolean isDeletion, final boolean isBeforeDeletion, final boolean isBeforeInsertion, final boolean isNextToSoftClip, + final String nextEventBases, final int nextEventLength) { if (offset < 0 && isDeletion) throw new ReviewedStingException("Pileup Element cannot create a deletion with a negative offset"); @@ -52,6 +59,14 @@ public class PileupElement implements Comparable { this.isBeforeDeletion = isBeforeDeletion; this.isBeforeInsertion = isBeforeInsertion; this.isNextToSoftClip = isNextToSoftClip; + if (isBeforeInsertion) + eventBases = nextEventBases; + else + eventBases = null; // ignore argument in any other case + if (isBeforeDeletion || isBeforeInsertion) + eventLength = nextEventLength; + else + eventLength = -1; } public boolean isDeletion() { @@ -104,6 +119,20 @@ public class PileupElement implements Comparable { return getBaseDeletionQual(offset); } + /** + * Returns length of the event (number of inserted or deleted bases + */ + public int getEventLength() { + return eventLength; + } + + /** + * Returns actual sequence of inserted bases, or a null if the event is a deletion or if there is no event in the associated read. + */ + public String getEventBases() { + return eventBases; + } + public int getMappingQual() { return read.getMappingQuality(); } diff --git a/public/java/src/org/broadinstitute/sting/utils/pileup/ReadBackedPileupImpl.java b/public/java/src/org/broadinstitute/sting/utils/pileup/ReadBackedPileupImpl.java index 7a6ebef21..759d64b2f 100644 --- a/public/java/src/org/broadinstitute/sting/utils/pileup/ReadBackedPileupImpl.java +++ b/public/java/src/org/broadinstitute/sting/utils/pileup/ReadBackedPileupImpl.java @@ -71,7 +71,13 @@ public class ReadBackedPileupImpl extends AbstractReadBackedPileup