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

Officially removing all code associated with extended events. Note that I still have a longer term project on my plate to refactor the ReadBackedPileup, but that's a much larger effort.

This commit is contained in:
Eric Banks 2012-06-15 15:55:03 -04:00
parent 783b7f6899
commit 677babf546
35 changed files with 183 additions and 1108 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
public/java/src/org/broadinstitute/sting/gatk/ReadProperties.java
View File

@ -51,11 +51,6 @@ public class ReadProperties {
return includeReadsWithDeletionAtLoci;
}
@Deprecated
public boolean generateExtendedEvents() {
return false;
}
/**
* Gets a list of the files acting as sources of reads.
* @return A list of files storing reads data.

28
public/java/src/org/broadinstitute/sting/gatk/contexts/AlignmentContext.java
View File

@ -28,7 +28,6 @@ package org.broadinstitute.sting.gatk.contexts;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.HasGenomeLocation;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.broadinstitute.sting.utils.pileup.ReadBackedExtendedEventPileup;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
@ -89,36 +88,9 @@ public class AlignmentContext implements HasGenomeLocation {
* @return
*/
public ReadBackedPileup getBasePileup() {
if(!hasBasePileup())
throw new ReviewedStingException("No base pileup is available. Please check for a base pileup with hasBasePileup() before attempting to retrieve a pileup.");
return basePileup;
}
/** Returns extended event (indel) pileup over the current genomic location. May return null if this context keeps
* only base pileup.
* @return
*/
@Deprecated
public ReadBackedExtendedEventPileup getExtendedEventPileup() {
if(!hasExtendedEventPileup())
throw new ReviewedStingException("No extended event pileup is present.");
return (ReadBackedExtendedEventPileup)basePileup;
}
/**
* Returns true if this alignment context keeps base pileup over the current genomic location.
* TODO: Syntax of AlignmentContext uses hasBasePileup() / hasExtendedEventPileup() as an enumeration mechanism. Change this to a more sensible interface.
* @return
*/
public boolean hasBasePileup() { return !(basePileup instanceof ReadBackedExtendedEventPileup); }
/** Returns true if this alignment context keeps extended event (indel) pileup over the current genomic location.
*
* @return
*/
@Deprecated
public boolean hasExtendedEventPileup() { return basePileup instanceof ReadBackedExtendedEventPileup; }
/**
* Returns true if any reads have been filtered out of the pileup due to excess DoC.
* @return True if reads have been filtered out. False otherwise.

44
public/java/src/org/broadinstitute/sting/gatk/contexts/AlignmentContextUtils.java
View File

@ -116,19 +116,15 @@ public class AlignmentContextUtils {
*
**/
public static Map<SAMReadGroupRecord, AlignmentContext> splitContextByReadGroup(AlignmentContext context, Collection<SAMReadGroupRecord> readGroups) {
if ( ! context.hasBasePileup() ) {
return Collections.emptyMap();
} else {
HashMap<SAMReadGroupRecord, AlignmentContext> contexts = new HashMap<SAMReadGroupRecord, AlignmentContext>();
HashMap<SAMReadGroupRecord, AlignmentContext> contexts = new HashMap<SAMReadGroupRecord, AlignmentContext>();
for (SAMReadGroupRecord rg : readGroups) {
ReadBackedPileup rgPileup = context.getBasePileup().getPileupForReadGroup(rg.getReadGroupId());
if ( rgPileup != null ) // there we some reads for RG
contexts.put(rg, new AlignmentContext(context.getLocation(), rgPileup));
}
return contexts;
for (SAMReadGroupRecord rg : readGroups) {
ReadBackedPileup rgPileup = context.getBasePileup().getPileupForReadGroup(rg.getReadGroupId());
if ( rgPileup != null ) // there we some reads for RG
contexts.put(rg, new AlignmentContext(context.getLocation(), rgPileup));
}
return contexts;
}
public static Map<String, AlignmentContext> splitContextBySampleName(ReadBackedPileup pileup) {
@ -139,32 +135,16 @@ public class AlignmentContextUtils {
public static AlignmentContext joinContexts(Collection<AlignmentContext> contexts) {
// validation
GenomeLoc loc = contexts.iterator().next().getLocation();
boolean isExtended = contexts.iterator().next().basePileup instanceof ReadBackedExtendedEventPileup;
for(AlignmentContext context: contexts) {
if(!loc.equals(context.getLocation()))
throw new ReviewedStingException("Illegal attempt to join contexts from different genomic locations");
if(isExtended != (context.basePileup instanceof ReadBackedExtendedEventPileup))
throw new ReviewedStingException("Illegal attempt to join simple and extended contexts");
}
AlignmentContext jointContext;
if(isExtended) {
List<ExtendedEventPileupElement> pe = new ArrayList<ExtendedEventPileupElement>();
for(AlignmentContext context: contexts) {
for(PileupElement pileupElement: context.basePileup)
pe.add((ExtendedEventPileupElement)pileupElement);
}
jointContext = new AlignmentContext(loc, new ReadBackedExtendedEventPileupImpl(loc,pe));
List<PileupElement> pe = new ArrayList<PileupElement>();
for(AlignmentContext context: contexts) {
for(PileupElement pileupElement: context.basePileup)
pe.add(pileupElement);
}
else {
List<PileupElement> pe = new ArrayList<PileupElement>();
for(AlignmentContext context: contexts) {
for(PileupElement pileupElement: context.basePileup)
pe.add(pileupElement);
}
jointContext = new AlignmentContext(loc, new ReadBackedPileupImpl(loc,pe));
}
return jointContext;
return new AlignmentContext(loc, new ReadBackedPileupImpl(loc,pe));
}
}

273
public/java/src/org/broadinstitute/sting/gatk/iterators/LocusIteratorByState.java
View File

@ -40,9 +40,7 @@ import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.ReservoirDownsampler;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.pileup.ExtendedEventPileupElement;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.pileup.ReadBackedExtendedEventPileupImpl;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileupImpl;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
import org.broadinstitute.sting.utils.sam.ReadUtils;
@ -63,7 +61,6 @@ public class LocusIteratorByState extends LocusIterator {
// member fields
//
// -----------------------------------------------------------------------------------------------------------------
private boolean hasExtendedEvents = false; // will be set to true if at least one read had an indel right before the current position
/**
* Used to create new GenomeLocs.
@ -92,26 +89,10 @@ public class LocusIteratorByState extends LocusIterator {
// stepForwardOnGenome(). The next call to stepForwardOnGenome() will clear that memory (as we remember only extended
// events immediately preceding the current reference base).
boolean generateExtendedEvents = true; // should we generate an additional, special pile for indels between the ref bases?
// the only purpose of this flag is to shield away a few additional lines of code
// when extended piles are not needed, it may not be even worth it...
byte[] insertedBases = null; // remember full inserted sequence if we are generating piles of extended events (indels)
int eventLength = -1; // will be set to the length of insertion/deletion if we are generating piles of extended events
byte eventDelayedFlag = 0; // will be set to non-0 if there was an event (indel) right before the
// current base on the ref. We use a counter-like variable here since clearing the indel event is
// delayed by one base, so we need to remember how long ago we have seen the actual event
int eventStart = -1; // where on the read the extended event starts (i.e. the last position on the read prior to the
// event, or -1 if alignment starts with an insertion); this one is easy to recompute on the fly,
// we cache it here mainly for convenience
public SAMRecordState(SAMRecord read, boolean extended) {
public SAMRecordState(SAMRecord read) {
this.read = read;
cigar = read.getCigar();
nCigarElements = cigar.numCigarElements();
generateExtendedEvents = extended;
//System.out.printf("Creating a SAMRecordState: %s%n", this);
}
@ -150,27 +131,6 @@ public class LocusIteratorByState extends LocusIterator {
return curElement.getOperator();
}
/**
* Returns true if we just stepped over insertion/into a deletion prior to the last return from stepForwardOnGenome.
*
* @return
*/
public boolean hadIndel() {
return (eventLength > 0);
}
public int getEventLength() {
return eventLength;
}
public byte[] getEventBases() {
return insertedBases;
}
public int getReadEventStartOffset() {
return eventStart;
}
public String toString() {
return String.format("%s ro=%d go=%d co=%d cec=%d %s", read.getReadName(), readOffset, genomeOffset, cigarOffset, cigarElementCounter, curElement);
}
@ -208,19 +168,6 @@ public class LocusIteratorByState extends LocusIterator {
genomeOffset++; // extended events need that. Logically, it's legal to advance the genomic offset here:
// we do step forward on the ref, and by returning null we also indicate that we are past the read end.
if (generateExtendedEvents && eventDelayedFlag > 0) {
// if we had an indel right before the read ended (i.e. insertion was the last cigar element),
// we keep it until next reference base; then we discard it and this will allow the LocusIterator to
// finally discard this read
eventDelayedFlag--;
if (eventDelayedFlag == 0) {
eventLength = -1; // reset event when we are past it
insertedBases = null;
eventStart = -1;
}
}
return null;
}
}
@ -232,17 +179,6 @@ public class LocusIteratorByState extends LocusIterator {
cigarElementCounter = curElement.getLength();
break;
case I: // insertion w.r.t. the reference
if (generateExtendedEvents) {
// we see insertions only once, when we step right onto them; the position on the read is scrolled
// past the insertion right after that
if (eventDelayedFlag > 1)
throw new UserException.MalformedBAM(read, String.format("Adjacent I/D events in read %s -- cigar: %s. This is an indication of a malformed file, but the SAM spec allows reads with adjacent insertion/deletion. If you are sure you want to use this read, re-run your analysis with the extra option: -rf BadCigar", read.getReadName(), read.getCigarString()));
insertedBases = Arrays.copyOfRange(read.getReadBases(), readOffset + 1, readOffset + 1 + curElement.getLength());
eventLength = curElement.getLength();
eventStart = readOffset;
eventDelayedFlag = 2; // insertion causes re-entry into stepForwardOnGenome, so we set the delay to 2
// System.out.println("Inserted "+(new String (insertedBases)) +" after "+readOffset);
} // continue onto the 'S' case !
case S: // soft clip
cigarElementCounter = curElement.getLength();
readOffset += curElement.getLength();
@ -250,19 +186,6 @@ public class LocusIteratorByState extends LocusIterator {
case D: // deletion w.r.t. the reference
if (readOffset < 0) // we don't want reads starting with deletion, this is a malformed cigar string
throw new UserException.MalformedBAM(read, "Read starting with deletion. Cigar: " + read.getCigarString() + ". This is an indication of a malformed file, but the SAM spec allows reads starting in deletion. If you are sure you want to use this read, re-run your analysis with the extra option: -rf BadCigar");
if (generateExtendedEvents) {
if (cigarElementCounter == 1) {
// generate an extended event only if we just stepped into the deletion (i.e. don't
// generate the event at every deleted position on the ref, that's what cigarElementCounter==1 is for!)
if (eventDelayedFlag > 1)
throw new UserException.MalformedBAM(read, String.format("Adjacent I/D events in read %s -- cigar: %s. This is an indication of a malformed file, but the SAM spec allows reads with adjacent insertion/deletion. If you are sure you want to use this read, re-run your analysis with the extra option: -rf BadCigar", read.getReadName(), read.getCigarString()));
eventLength = curElement.getLength();
eventDelayedFlag = 2; // deletion on the ref causes an immediate return, so we have to delay by 1 only
eventStart = readOffset;
insertedBases = null;
// System.out.println("Deleted "+eventLength +" bases after "+readOffset);
}
}
// should be the same as N case
genomeOffset++;
done = true;
@ -280,21 +203,6 @@ public class LocusIteratorByState extends LocusIterator {
throw new IllegalStateException("Case statement didn't deal with cigar op: " + curElement.getOperator());
}
if (generateExtendedEvents) {
if (eventDelayedFlag > 0 && done) {
// if we did make a successful step on the ref, decrement delayed flag. If, upon the decrementing the,
// the flag is 1, we are standing on the reference base right after the indel (so we have to keep it).
// Otherwise, we are away from the previous indel and have to clear our memories...
eventDelayedFlag--; // when we notice an indel, we set delayed flag to 2, so now
// if eventDelayedFlag == 1, an indel occured right before the current base
if (eventDelayedFlag == 0) {
eventLength = -1; // reset event when we are past it
insertedBases = null;
eventStart = -1;
}
}
}
return done ? curElement.getOperator() : stepForwardOnGenome();
}
}
@ -374,147 +282,69 @@ public class LocusIteratorByState extends LocusIterator {
// this call will set hasExtendedEvents to true if it picks up a read with indel right before the current position on the ref:
readStates.collectPendingReads();
int size = 0;
int nDeletions = 0;
int nInsertions = 0;
int nMQ0Reads = 0;
final GenomeLoc location = getLocation();
final Map<String, ReadBackedPileupImpl> fullPileup = new HashMap<String, ReadBackedPileupImpl>();
boolean hasBeenSampled = false;
for (final String sample : samples) {
final Iterator<SAMRecordState> iterator = readStates.iterator(sample);
final List<PileupElement> pile = new ArrayList<PileupElement>(readStates.size(sample));
hasBeenSampled |= location.getStart() <= readStates.getDownsamplingExtent(sample);
int size = 0; // number of elements in this sample's pileup
int nDeletions = 0; // number of deletions in this sample's pileup
int nMQ0Reads = 0; // number of MQ0 reads in this sample's pileup (warning: current implementation includes N bases that are MQ0)
// if extended events are requested, and if previous traversal step brought us over an indel in
// at least one read, we emit extended pileup (making sure that it is associated with the previous base,
// i.e. the one right *before* the indel) and do NOT shift the current position on the ref.
// In this case, the subsequent call to next() will emit the normal pileup at the current base
// and shift the position.
if (readInfo.generateExtendedEvents() && hasExtendedEvents) {
Map<String, ReadBackedExtendedEventPileupImpl> fullExtendedEventPileup = new HashMap<String, ReadBackedExtendedEventPileupImpl>();
while (iterator.hasNext()) {
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 CigarElement nextElement = state.peekForwardOnGenome(); // next cigar element
final CigarElement lastElement = state.peekBackwardOnGenome(); // last cigar element
final CigarOperator nextOp = nextElement.getOperator(); // next cigar operator
final CigarOperator lastOp = lastElement.getOperator(); // last cigar operator
final int readOffset = state.getReadOffset(); // the base offset on this read
// get current location on the reference and decrement it by 1: the indels we just stepped over
// are associated with the *previous* reference base
GenomeLoc loc = genomeLocParser.incPos(getLocation(), -1);
final boolean isBeforeDeletion = nextOp == CigarOperator.DELETION;
final boolean isAfterDeletion = lastOp == CigarOperator.DELETION;
final boolean isBeforeInsertion = nextOp == CigarOperator.INSERTION;
final boolean isAfterInsertion = lastOp == CigarOperator.INSERTION;
final boolean isNextToSoftClip = nextOp == CigarOperator.S || (state.getGenomeOffset() == 0 && read.getSoftStart() != read.getAlignmentStart());
boolean hasBeenSampled = false;
for (final String sample : samples) {
Iterator<SAMRecordState> iterator = readStates.iterator(sample);
List<ExtendedEventPileupElement> indelPile = new ArrayList<ExtendedEventPileupElement>(readStates.size(sample));
hasBeenSampled |= loc.getStart() <= readStates.getDownsamplingExtent(sample);
int nextElementLength = nextElement.getLength();
size = 0;
nDeletions = 0;
nInsertions = 0;
nMQ0Reads = 0;
int maxDeletionLength = 0;
if (op == CigarOperator.N) // N's are never added to any pileup
continue;
while (iterator.hasNext()) {
final SAMRecordState state = iterator.next();
final GATKSAMRecord read = (GATKSAMRecord) state.getRead(); // the actual read
final CigarOperator op = state.getCurrentCigarOperator(); // current cigar operator
final int readOffset = state.getReadOffset(); // the base offset on this read
final int eventStartOffset = state.getReadEventStartOffset(); // this will be -1 if base is not a deletion, or if base is the first deletion in the event. Otherwise, it will give the last base before the deletion began.
final int eventLength = state.getEventLength();
if (op == CigarOperator.N) // N's are never added to any pileup
continue;
if (state.hadIndel()) { // this read has an indel associated with the previous position on the ref
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, isBeforeDeletion, isAfterDeletion, isBeforeInsertion, isAfterInsertion, isNextToSoftClip, null, nextOp == CigarOperator.D ? nextElementLength : -1));
size++;
ExtendedEventPileupElement pileupElement;
if (state.getEventBases() == null) { // Deletion event
nDeletions++;
maxDeletionLength = Math.max(maxDeletionLength, state.getEventLength());
pileupElement = new ExtendedEventPileupElement(read, eventStartOffset, eventLength);
}
else { // Insertion event
nInsertions++;
pileupElement = new ExtendedEventPileupElement(read, eventStartOffset, eventLength, state.getEventBases());
}
if (read.getMappingQuality() == 0)
nMQ0Reads++;
indelPile.add(pileupElement);
}
// this read has no indel so add it to the pileup as a NOEVENT:
// a deletion that didn't start here (therefore, not an extended event)
// we add (mis)matches as no events.
else if (op != CigarOperator.D || readInfo.includeReadsWithDeletionAtLoci()) {
size++;
indelPile.add(new ExtendedEventPileupElement((GATKSAMRecord) state.getRead(), readOffset));
nDeletions++;
if (read.getMappingQuality() == 0)
nMQ0Reads++;
}
}
else {
if (!filterBaseInRead(read, location.getStart())) {
String insertedBaseString = null;
if (nextOp == CigarOperator.I)
insertedBaseString = new String(Arrays.copyOfRange(read.getReadBases(), readOffset + 1, readOffset + 1 + nextElement.getLength()));
if (indelPile.size() != 0)
fullExtendedEventPileup.put(sample, new ReadBackedExtendedEventPileupImpl(loc, indelPile, size, maxDeletionLength, nInsertions, nDeletions, nMQ0Reads));
}
hasExtendedEvents = false; // we are done with extended events prior to current ref base
nextAlignmentContext = new AlignmentContext(loc, new ReadBackedExtendedEventPileupImpl(loc, fullExtendedEventPileup), hasBeenSampled);
}
else { // this is a regular event pileup (not extended)
final GenomeLoc location = getLocation();
final Map<String, ReadBackedPileupImpl> fullPileup = new HashMap<String, ReadBackedPileupImpl>();
boolean hasBeenSampled = false;
for (final String sample : samples) {
final Iterator<SAMRecordState> iterator = readStates.iterator(sample);
final List<PileupElement> pile = new ArrayList<PileupElement>(readStates.size(sample));
hasBeenSampled |= location.getStart() <= readStates.getDownsamplingExtent(sample);
size = 0; // number of elements in this sample's pileup
nDeletions = 0; // number of deletions in this sample's pileup
nMQ0Reads = 0; // number of MQ0 reads in this sample's pileup (warning: current implementation includes N bases that are MQ0)
while (iterator.hasNext()) {
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 CigarElement nextElement = state.peekForwardOnGenome(); // next cigar element
final CigarElement lastElement = state.peekBackwardOnGenome(); // last cigar element
final CigarOperator nextOp = nextElement.getOperator(); // next cigar operator
final CigarOperator lastOp = lastElement.getOperator(); // last cigar operator
final int readOffset = state.getReadOffset(); // the base offset on this read
final boolean isBeforeDeletion = nextOp == CigarOperator.DELETION;
final boolean isAfterDeletion = lastOp == CigarOperator.DELETION;
final boolean isBeforeInsertion = nextOp == CigarOperator.INSERTION;
final boolean isAfterInsertion = lastOp == CigarOperator.INSERTION;
final boolean isNextToSoftClip = nextOp == CigarOperator.S || (state.getGenomeOffset() == 0 && read.getSoftStart() != read.getAlignmentStart());
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, isBeforeDeletion, isAfterDeletion, isBeforeInsertion, isAfterInsertion, isNextToSoftClip, null, nextOp == CigarOperator.D ? nextElementLength : -1));
size++;
nDeletions++;
if (read.getMappingQuality() == 0)
nMQ0Reads++;
}
}
else {
if (!filterBaseInRead(read, location.getStart())) {
String insertedBaseString = null;
if (nextOp == CigarOperator.I)
insertedBaseString = new String(Arrays.copyOfRange(read.getReadBases(), readOffset + 1, readOffset + 1 + nextElement.getLength()));
pile.add(new PileupElement(read, readOffset, false, isBeforeDeletion, isAfterDeletion, isBeforeInsertion, isAfterInsertion, isNextToSoftClip, insertedBaseString, nextElementLength));
size++;
if (read.getMappingQuality() == 0)
nMQ0Reads++;
}
pile.add(new PileupElement(read, readOffset, false, isBeforeDeletion, isAfterDeletion, isBeforeInsertion, isAfterInsertion, isNextToSoftClip, insertedBaseString, nextElementLength));
size++;
if (read.getMappingQuality() == 0)
nMQ0Reads++;
}
}
if (pile.size() != 0) // if this pileup added at least one base, add it to the full pileup
fullPileup.put(sample, new ReadBackedPileupImpl(location, pile, size, nDeletions, nMQ0Reads));
}
updateReadStates(); // critical - must be called after we get the current state offsets and location
if (!fullPileup.isEmpty()) // if we got reads with non-D/N over the current position, we are done
nextAlignmentContext = new AlignmentContext(location, new ReadBackedPileupImpl(location, fullPileup), hasBeenSampled);
if (pile.size() != 0) // if this pileup added at least one base, add it to the full pileup
fullPileup.put(sample, new ReadBackedPileupImpl(location, pile, size, nDeletions, nMQ0Reads));
}
updateReadStates(); // critical - must be called after we get the current state offsets and location
if (!fullPileup.isEmpty()) // if we got reads with non-D/N over the current position, we are done
nextAlignmentContext = new AlignmentContext(location, new ReadBackedPileupImpl(location, fullPileup), hasBeenSampled);
}
}
@ -546,9 +376,7 @@ public class LocusIteratorByState extends LocusIterator {
while (it.hasNext()) {
SAMRecordState state = it.next();
CigarOperator op = state.stepForwardOnGenome();
if (state.hadIndel() && readInfo.generateExtendedEvents())
hasExtendedEvents = true;
else if (op == null) {
if (op == null) {
// we discard the read only when we are past its end AND indel at the end of the read (if any) was
// already processed. Keeping the read state that retunred null upon stepForwardOnGenome() is safe
// as the next call to stepForwardOnGenome() will return null again AND will clear hadIndel() flag.
@ -757,12 +585,9 @@ public class LocusIteratorByState extends LocusIterator {
int readCount = 0;
for (SAMRecord read : reads) {
if (readCount < maxReads) {
SAMRecordState state = new SAMRecordState(read, readInfo.generateExtendedEvents());
SAMRecordState state = new SAMRecordState(read);
state.stepForwardOnGenome();
newReadStates.add(state);
// TODO: What if we downsample the extended events away?
if (state.hadIndel())
hasExtendedEvents = true;
readCount++;
}
}

13
public/java/src/org/broadinstitute/sting/gatk/traversals/TraverseLoci.java
View File

@ -53,19 +53,6 @@ public class TraverseLoci<M,T> extends TraversalEngine<M,T,LocusWalker<M,T>,Locu
dataProvider.getShard().getReadMetrics().incrementNumIterations();
if ( locus.hasExtendedEventPileup() ) {
// if the alignment context we received holds an "extended" pileup (i.e. pileup of insertions/deletions
// associated with the current site), we need to update the location. The updated location still starts
// at the current genomic position, but it has to span the length of the longest deletion (if any).
location = engine.getGenomeLocParser().setStop(location,location.getStop()+locus.getExtendedEventPileup().getMaxDeletionLength());
// it is possible that the new expanded location spans the current shard boundary; the next method ensures
// that when it is the case, the reference sequence held by the ReferenceView will be reloaded so that
// the view has all the bases we are gonna need. If the location fits within the current view bounds,
// the next call will not do anything to the view:
referenceView.expandBoundsToAccomodateLoc(location);
}
// create reference context. Note that if we have a pileup of "extended events", the context will
// hold the (longest) stretch of deleted reference bases (if deletions are present in the pileup).
ReferenceContext refContext = referenceView.getReferenceContext(location);

14
public/java/src/org/broadinstitute/sting/gatk/walkers/PileupWalker.java
View File

@ -34,9 +34,7 @@ import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.collections.Pair;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.pileup.ReadBackedExtendedEventPileup;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import java.io.PrintStream;
@ -79,13 +77,11 @@ public class PileupWalker extends LocusWalker<Integer, Integer> implements TreeR
String rods = getReferenceOrderedData( tracker );
if ( context.hasBasePileup() ) {
ReadBackedPileup basePileup = context.getBasePileup();
out.printf("%s %s", basePileup.getPileupString((char)ref.getBase()), rods);
if ( SHOW_VERBOSE )
out.printf(" %s", createVerboseOutput(basePileup));
out.println();
}
ReadBackedPileup basePileup = context.getBasePileup();
out.printf("%s %s", basePileup.getPileupString((char)ref.getBase()), rods);
if ( SHOW_VERBOSE )
out.printf(" %s", createVerboseOutput(basePileup));
out.println();
return 1;
}

4
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/AlleleBalance.java
View File

@ -30,11 +30,9 @@ import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.AnnotatorCompatibleWalker;
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.InfoFieldAnnotation;
import org.broadinstitute.sting.utils.BaseUtils;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.codecs.vcf.VCFHeaderLineType;
import org.broadinstitute.sting.utils.codecs.vcf.VCFInfoHeaderLine;
import org.broadinstitute.sting.utils.pileup.ReadBackedExtendedEventPileup;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.variantcontext.Genotype;
import org.broadinstitute.sting.utils.variantcontext.GenotypesContext;
@ -72,7 +70,7 @@ public class AlleleBalance extends InfoFieldAnnotation {
// we care only about het calls
AlignmentContext context = stratifiedContexts.get(genotype.getSampleName());
if ( context == null || !context.hasBasePileup() )
if ( context == null )
continue;
final ReadBackedPileup pileup = context.getBasePileup();

3
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/AlleleBalanceBySample.java
View File

@ -49,9 +49,6 @@ public class AlleleBalanceBySample extends GenotypeAnnotation implements Experim
if ( altAlleles.size() == 0 )
return null;
if ( !stratifiedContext.hasBasePileup() )
return null;
final String bases = new String(stratifiedContext.getBasePileup().getBases());
if ( bases.length() == 0 )
return null;

2
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/BaseCounts.java
View File

@ -59,8 +59,6 @@ public class BaseCounts extends InfoFieldAnnotation {
int[] counts = new int[4];
for ( Map.Entry<String, AlignmentContext> sample : stratifiedContexts.entrySet() ) {
if ( !sample.getValue().hasBasePileup() )
continue;
for (byte base : sample.getValue().getBasePileup().getBases() ) {
int index = BaseUtils.simpleBaseToBaseIndex(base);
if ( index != -1 )

2
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthOfCoverage.java
View File

@ -44,7 +44,7 @@ public class DepthOfCoverage extends InfoFieldAnnotation implements StandardAnno
int depth = 0;
for ( Map.Entry<String, AlignmentContext> sample : stratifiedContexts.entrySet() )
depth += sample.getValue().hasBasePileup() ? sample.getValue().getBasePileup().depthOfCoverage() : 0;
depth += sample.getValue().getBasePileup().depthOfCoverage();
Map<String, Object> map = new HashMap<String, Object>();
map.put(getKeyNames().get(0), String.format("%d", depth));
return map;

7
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthPerAlleleBySample.java
View File

@ -1,6 +1,5 @@
package org.broadinstitute.sting.gatk.walkers.annotator;
import org.apache.commons.lang.ArrayUtils;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
@ -58,9 +57,6 @@ public class DepthPerAlleleBySample extends GenotypeAnnotation implements Standa
private void annotateSNP(AlignmentContext stratifiedContext, VariantContext vc, GenotypeBuilder gb) {
if ( ! stratifiedContext.hasBasePileup() )
return;
HashMap<Byte, Integer> alleleCounts = new HashMap<Byte, Integer>();
for ( Allele allele : vc.getAlleles() )
alleleCounts.put(allele.getBases()[0], 0);
@ -81,9 +77,6 @@ public class DepthPerAlleleBySample extends GenotypeAnnotation implements Standa
}
private void annotateIndel(AlignmentContext stratifiedContext, VariantContext vc, GenotypeBuilder gb) {
if ( ! stratifiedContext.hasBasePileup() )
return;
ReadBackedPileup pileup = stratifiedContext.getBasePileup();
if ( pileup == null )
return;

2
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/FisherStrand.java
View File

@ -296,7 +296,7 @@ public class FisherStrand extends InfoFieldAnnotation implements StandardAnnotat
for ( String sample : stratifiedContexts.keySet() ) {
final AlignmentContext context = stratifiedContexts.get(sample);
if ( context == null || !context.hasBasePileup() )
if ( context == null )
continue;
final ReadBackedPileup pileup = context.getBasePileup();

5
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/HaplotypeScore.java
View File

@ -74,9 +74,6 @@ public class HaplotypeScore extends InfoFieldAnnotation implements StandardAnnot
final int locus = ref.getLocus().getStart() + (ref.getLocus().getStop() - ref.getLocus().getStart()) / 2;
if ( !context.hasBasePileup() )
return null;
final ReadBackedPileup pileup = context.getBasePileup();
// Compute all haplotypes consistent with the current read pileup
@ -86,7 +83,7 @@ public class HaplotypeScore extends InfoFieldAnnotation implements StandardAnnot
if (haplotypes != null) {
for (final Genotype genotype : vc.getGenotypes()) {
final AlignmentContext thisContext = stratifiedContexts.get(genotype.getSampleName());
if (thisContext != null && thisContext.hasBasePileup()) {
if (thisContext != null) {
final ReadBackedPileup thisPileup = thisContext.getBasePileup();
if (vc.isSNP())
scoreRA.add(scoreReadsAgainstHaplotypes(haplotypes, thisPileup, contextSize, locus)); // Taking the simple average of all sample's score since the score can be negative and the RMS doesn't make sense

3
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/LowMQ.java
View File

@ -31,9 +31,6 @@ public class LowMQ extends InfoFieldAnnotation {
double total = 0;
for ( Map.Entry<String, AlignmentContext> sample : stratifiedContexts.entrySet() )
{
if ( !sample.getValue().hasBasePileup() )
continue;
for ( PileupElement p : sample.getValue().getBasePileup() )
{
if ( p.getMappingQual() == 0 ) { mq0 += 1; }

10
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/MappingQualityZero.java
View File

@ -31,12 +31,10 @@ public class MappingQualityZero extends InfoFieldAnnotation implements StandardA
int mq0 = 0;
for ( Map.Entry<String, AlignmentContext> sample : stratifiedContexts.entrySet() ) {
final AlignmentContext context = sample.getValue();
if ( context.hasBasePileup() ) {
final ReadBackedPileup pileup = context.getBasePileup();
for (PileupElement p : pileup ) {
if ( p.getMappingQual() == 0 )
mq0++;
}
final ReadBackedPileup pileup = context.getBasePileup();
for (PileupElement p : pileup ) {
if ( p.getMappingQual() == 0 )
mq0++;
}
}
Map<String, Object> map = new HashMap<String, Object>();

12
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/MappingQualityZeroBySample.java
View File

@ -40,9 +40,7 @@ import org.broadinstitute.sting.utils.variantcontext.GenotypeBuilder;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import java.util.Arrays;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
/**
* Count for each sample of mapping quality zero reads
@ -55,12 +53,10 @@ public class MappingQualityZeroBySample extends GenotypeAnnotation {
return;
int mq0 = 0;
if ( context.hasBasePileup() ) {
final ReadBackedPileup pileup = context.getBasePileup();
for (PileupElement p : pileup ) {
if ( p.getMappingQual() == 0 )
mq0++;
}
final ReadBackedPileup pileup = context.getBasePileup();
for (PileupElement p : pileup ) {
if ( p.getMappingQual() == 0 )
mq0++;
}
gb.attribute(getKeyNames().get(0), mq0);

10
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/MappingQualityZeroFraction.java
View File

@ -31,12 +31,10 @@ public class MappingQualityZeroFraction extends InfoFieldAnnotation implements E
for ( Map.Entry<String, AlignmentContext> sample : stratifiedContexts.entrySet() ) {
AlignmentContext context = sample.getValue();
depth += context.size();
if ( context.hasBasePileup() ) {
final ReadBackedPileup pileup = context.getBasePileup();
for (PileupElement p : pileup ) {
if ( p.getMappingQual() == 0 )
mq0++;
}
final ReadBackedPileup pileup = context.getBasePileup();
for (PileupElement p : pileup ) {
if ( p.getMappingQual() == 0 )
mq0++;
}
}
if (depth > 0) {

16
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/NBaseCount.java
View File

@ -28,15 +28,13 @@ public class NBaseCount extends InfoFieldAnnotation {
int countRegularBaseSolid = 0;
for( final AlignmentContext context : stratifiedContexts.values() ) {
if ( context.hasBasePileup() ) { // must be called as getBasePileup may throw error when pileup has no bases
for( final PileupElement p : context.getBasePileup()) {
final String platform = p.getRead().getReadGroup().getPlatform();
if( platform != null && platform.toUpperCase().contains("SOLID") ) {
if( BaseUtils.isNBase( p.getBase() ) ) {
countNBaseSolid++;
} else if( BaseUtils.isRegularBase( p.getBase() ) ) {
countRegularBaseSolid++;
}
for( final PileupElement p : context.getBasePileup()) {
final String platform = p.getRead().getReadGroup().getPlatform();
if( platform != null && platform.toUpperCase().contains("SOLID") ) {
if( BaseUtils.isNBase( p.getBase() ) ) {
countNBaseSolid++;
} else if( BaseUtils.isRegularBase( p.getBase() ) ) {
countRegularBaseSolid++;
}
}
}

2
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/QualByDepth.java
View File

@ -48,7 +48,7 @@ public class QualByDepth extends InfoFieldAnnotation implements StandardAnnotati
if ( context == null )
continue;
depth += context.hasBasePileup() ? context.getBasePileup().depthOfCoverage() : 0;
depth += context.getBasePileup().depthOfCoverage();
}
if ( depth == 0 )

10
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/RMSMappingQuality.java
View File

@ -42,12 +42,10 @@ public class RMSMappingQuality extends InfoFieldAnnotation implements StandardAn
for ( Map.Entry<String, AlignmentContext> sample : stratifiedContexts.entrySet() ) {
AlignmentContext context = sample.getValue();
if ( context.hasBasePileup() ) {
final ReadBackedPileup pileup = context.getBasePileup();
for (PileupElement p : pileup ) {
if ( p.getMappingQual() != QualityUtils.MAPPING_QUALITY_UNAVAILABLE )
qualities[index++] = p.getMappingQual();
}
final ReadBackedPileup pileup = context.getBasePileup();
for (PileupElement p : pileup ) {
if ( p.getMappingQual() != QualityUtils.MAPPING_QUALITY_UNAVAILABLE )
qualities[index++] = p.getMappingQual();
}
}

3
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/RankSumTest.java
View File

@ -63,9 +63,6 @@ public abstract class RankSumTest extends InfoFieldAnnotation implements Standar
continue;
}
if (!context.hasBasePileup())
continue;
final ReadBackedPileup pileup = context.getBasePileup();
if (pileup == null)
continue;

8
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/SpanningDeletions.java
View File

@ -35,11 +35,9 @@ public class SpanningDeletions extends InfoFieldAnnotation implements StandardAn
int depth = 0;
for ( Map.Entry<String, AlignmentContext> sample : stratifiedContexts.entrySet() ) {
AlignmentContext context = sample.getValue();
if ( context.hasBasePileup() ) {
final ReadBackedPileup pileup = context.getBasePileup();
deletions += pileup.getNumberOfDeletions();
depth += pileup.getNumberOfElements();
}
final ReadBackedPileup pileup = context.getBasePileup();
deletions += pileup.getNumberOfDeletions();
depth += pileup.getNumberOfElements();
}
Map<String, Object> map = new HashMap<String, Object>();
map.put(getKeyNames().get(0), String.format("%.2f", depth == 0 ? 0.0 : (double)deletions/(double)depth));

22
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/TechnologyComposition.java
View File

@ -39,18 +39,16 @@ public class TechnologyComposition extends InfoFieldAnnotation implements Experi
for ( Map.Entry<String, AlignmentContext> sample : stratifiedContexts.entrySet() ) {
AlignmentContext context = sample.getValue();
if ( context.hasBasePileup() ) {
final ReadBackedPileup pileup = context.getBasePileup();
for ( PileupElement p : pileup ) {
if(ReadUtils.is454Read(p.getRead()))
reads454++;
else if (ReadUtils.isSOLiDRead(p.getRead()))
readsSolid++;
else if (ReadUtils.isIlluminaRead(p.getRead()))
readsIllumina++;
else
readsOther++;
}
final ReadBackedPileup pileup = context.getBasePileup();
for ( PileupElement p : pileup ) {
if(ReadUtils.is454Read(p.getRead()))
reads454++;
else if (ReadUtils.isSOLiDRead(p.getRead()))
readsSolid++;
else if (ReadUtils.isIlluminaRead(p.getRead()))
readsIllumina++;
else
readsOther++;
}
}

10
public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/VariantAnnotator.java
View File

@ -305,12 +305,10 @@ public class VariantAnnotator extends RodWalker<Integer, Integer> implements Ann
// if the reference base is not ambiguous, we can annotate
Map<String, AlignmentContext> stratifiedContexts;
if ( BaseUtils.simpleBaseToBaseIndex(ref.getBase()) != -1 ) {
if ( context.hasBasePileup() ) {
stratifiedContexts = AlignmentContextUtils.splitContextBySampleName(context.getBasePileup());
annotatedVCs = new ArrayList<VariantContext>(VCs.size());
for ( VariantContext vc : VCs )
annotatedVCs.add(engine.annotateContext(tracker, ref, stratifiedContexts, vc));
}
stratifiedContexts = AlignmentContextUtils.splitContextBySampleName(context.getBasePileup());
annotatedVCs = new ArrayList<VariantContext>(VCs.size());
for ( VariantContext vc : VCs )
annotatedVCs.add(engine.annotateContext(tracker, ref, stratifiedContexts, vc));
}
for ( VariantContext annotatedVC : annotatedVCs )

11
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ConsensusAlleleCounter.java
View File

@ -98,11 +98,9 @@ public class ConsensusAlleleCounter {
for ( Map.Entry<String, AlignmentContext> sample : contexts.entrySet() ) {
final AlignmentContext context = AlignmentContextUtils.stratify(sample.getValue(), contextType);
if ( context.hasBasePileup() ) {
final ReadBackedPileup indelPileup = context.getBasePileup();
insCount += indelPileup.getNumberOfInsertionsAfterThisElement();
delCount += indelPileup.getNumberOfDeletionsAfterThisElement();
}
final ReadBackedPileup indelPileup = context.getBasePileup();
insCount += indelPileup.getNumberOfInsertionsAfterThisElement();
delCount += indelPileup.getNumberOfDeletionsAfterThisElement();
}
if ( insCount < minIndelCountForGenotyping && delCount < minIndelCountForGenotyping )
@ -112,9 +110,6 @@ public class ConsensusAlleleCounter {
// todo -- warning, can be duplicating expensive partition here
AlignmentContext context = AlignmentContextUtils.stratify(sample.getValue(), contextType);
if ( !context.hasBasePileup() )
continue;
final ReadBackedPileup indelPileup = context.getBasePileup();
final int nIndelReads = indelPileup.getNumberOfInsertionsAfterThisElement() + indelPileup.getNumberOfDeletionsAfterThisElement();

28
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java
View File

@ -30,13 +30,11 @@ import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.AlignmentContextUtils;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.Requires;
import org.broadinstitute.sting.gatk.walkers.indels.PairHMMIndelErrorModel;
import org.broadinstitute.sting.utils.BaseUtils;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.Haplotype;
import org.broadinstitute.sting.utils.codecs.vcf.VCFConstants;
import org.broadinstitute.sting.utils.collections.Pair;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
@ -138,21 +136,19 @@ public class IndelGenotypeLikelihoodsCalculationModel extends GenotypeLikelihood
for (Map.Entry<String, AlignmentContext> sample : contexts.entrySet()) {
AlignmentContext context = AlignmentContextUtils.stratify(sample.getValue(), contextType);
if (context.hasBasePileup()) {
final ReadBackedPileup pileup = context.getBasePileup();
if (pileup != null) {
final GenotypeBuilder b = new GenotypeBuilder(sample.getKey());
final double[] genotypeLikelihoods = pairModel.computeDiploidReadHaplotypeLikelihoods(pileup, haplotypeMap, ref, eventLength, getIndelLikelihoodMap());
b.PL(genotypeLikelihoods);
b.DP(getFilteredDepth(pileup));
genotypes.add(b.make());
final ReadBackedPileup pileup = context.getBasePileup();
if (pileup != null) {
final GenotypeBuilder b = new GenotypeBuilder(sample.getKey());
final double[] genotypeLikelihoods = pairModel.computeDiploidReadHaplotypeLikelihoods(pileup, haplotypeMap, ref, eventLength, getIndelLikelihoodMap());
b.PL(genotypeLikelihoods);
b.DP(getFilteredDepth(pileup));
genotypes.add(b.make());
if (DEBUG) {
System.out.format("Sample:%s Alleles:%s GL:", sample.getKey(), alleleList.toString());
for (int k = 0; k < genotypeLikelihoods.length; k++)
System.out.format("%1.4f ", genotypeLikelihoods[k]);
System.out.println();
}
if (DEBUG) {
System.out.format("Sample:%s Alleles:%s GL:", sample.getKey(), alleleList.toString());
for (int k = 0; k < genotypeLikelihoods.length; k++)
System.out.format("%1.4f ", genotypeLikelihoods[k]);
System.out.println();
}
}
}

64
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java
View File

@ -252,7 +252,7 @@ public class UnifiedGenotyperEngine {
vc = new VariantContextBuilder("UG_call", ref.getLocus().getContig(), ref.getLocus().getStart(), ref.getLocus().getStart(), alleles).make();
}
if ( annotationEngine != null && rawContext.hasBasePileup() ) {
if ( annotationEngine != null ) {
// Note: we want to use the *unfiltered* and *unBAQed* context for the annotations
final ReadBackedPileup pileup = rawContext.getBasePileup();
stratifiedContexts = AlignmentContextUtils.splitContextBySampleName(pileup);
@ -422,7 +422,7 @@ public class UnifiedGenotyperEngine {
if ( myAlleles.size() != vc.getAlleles().size() && !limitedContext ) // TODO - this function doesn't work with mixed records or records that started as mixed and then became non-mixed
vcCall = VariantContextUtils.reverseTrimAlleles(vcCall);
if ( annotationEngine != null && !limitedContext && rawContext.hasBasePileup() ) {
if ( annotationEngine != null && !limitedContext ) {
// Note: we want to use the *unfiltered* and *unBAQed* context for the annotations
final ReadBackedPileup pileup = rawContext.getBasePileup();
stratifiedContexts = AlignmentContextUtils.splitContextBySampleName(pileup);
@ -441,7 +441,7 @@ public class UnifiedGenotyperEngine {
private Map<String, AlignmentContext> getFilteredAndStratifiedContexts(UnifiedArgumentCollection UAC, ReferenceContext refContext, AlignmentContext rawContext, final GenotypeLikelihoodsCalculationModel.Model model) {
if ( !BaseUtils.isRegularBase(refContext.getBase()) || !rawContext.hasBasePileup() )
if ( !BaseUtils.isRegularBase(refContext.getBase()) )
return null;
Map<String, AlignmentContext> stratifiedContexts = null;
@ -507,9 +507,7 @@ public class UnifiedGenotyperEngine {
int depth = 0;
if ( isCovered ) {
AlignmentContext context = contexts.get(sample);
if ( context.hasBasePileup() )
depth = context.getBasePileup().depthOfCoverage();
depth = contexts.get(sample).getBasePileup().depthOfCoverage();
}
P_of_ref *= 1.0 - (theta / 2.0) * getRefBinomialProb(depth);
@ -571,37 +569,35 @@ public class UnifiedGenotyperEngine {
final List<GenotypeLikelihoodsCalculationModel.Model> models = new ArrayList<GenotypeLikelihoodsCalculationModel.Model>(2);
if ( rawContext.hasBasePileup() ) {
// if we're genotyping given alleles and we have a requested SNP at this position, do SNP
if ( UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES ) {
final VariantContext vcInput = UnifiedGenotyperEngine.getVCFromAllelesRod(tracker, refContext, rawContext.getLocation(), false, logger, UAC.alleles);
if ( vcInput == null )
return models;
// if we're genotyping given alleles and we have a requested SNP at this position, do SNP
if ( UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES ) {
final VariantContext vcInput = UnifiedGenotyperEngine.getVCFromAllelesRod(tracker, refContext, rawContext.getLocation(), false, logger, UAC.alleles);
if ( vcInput == null )
return models;
if ( vcInput.isSNP() ) {
// ignore SNPs if the user chose INDEL mode only
if ( UAC.GLmodel == GenotypeLikelihoodsCalculationModel.Model.BOTH )
models.add(GenotypeLikelihoodsCalculationModel.Model.SNP);
else if ( UAC.GLmodel.name().toUpperCase().contains("SNP") )
models.add(UAC.GLmodel);
}
else if ( vcInput.isIndel() || vcInput.isMixed() ) {
// ignore INDELs if the user chose SNP mode only
if ( UAC.GLmodel == GenotypeLikelihoodsCalculationModel.Model.BOTH )
models.add(GenotypeLikelihoodsCalculationModel.Model.INDEL);
else if (UAC.GLmodel.name().toUpperCase().contains("INDEL"))
models.add(UAC.GLmodel);
}
// No support for other types yet
if ( vcInput.isSNP() ) {
// ignore SNPs if the user chose INDEL mode only
if ( UAC.GLmodel == GenotypeLikelihoodsCalculationModel.Model.BOTH )
models.add(GenotypeLikelihoodsCalculationModel.Model.SNP);
else if ( UAC.GLmodel.name().toUpperCase().contains("SNP") )
models.add(UAC.GLmodel);
}
else if ( vcInput.isIndel() || vcInput.isMixed() ) {
// ignore INDELs if the user chose SNP mode only
if ( UAC.GLmodel == GenotypeLikelihoodsCalculationModel.Model.BOTH )
models.add(GenotypeLikelihoodsCalculationModel.Model.INDEL);
else if (UAC.GLmodel.name().toUpperCase().contains("INDEL"))
models.add(UAC.GLmodel);
}
// No support for other types yet
}
else {
if ( UAC.GLmodel == GenotypeLikelihoodsCalculationModel.Model.BOTH ) {
models.add(GenotypeLikelihoodsCalculationModel.Model.SNP);
models.add(GenotypeLikelihoodsCalculationModel.Model.INDEL);
}
else {
if ( UAC.GLmodel == GenotypeLikelihoodsCalculationModel.Model.BOTH ) {
models.add(GenotypeLikelihoodsCalculationModel.Model.SNP);
models.add(GenotypeLikelihoodsCalculationModel.Model.INDEL);
}
else {
models.add(UAC.GLmodel);
}
models.add(UAC.GLmodel);
}
}

48
public/java/src/org/broadinstitute/sting/gatk/walkers/indels/RealignerTargetCreator.java
View File

@ -185,38 +185,36 @@ public class RealignerTargetCreator extends RodWalker<RealignerTargetCreator.Eve
}
// look at the normal context to get deletions and positions with high entropy
if ( context.hasBasePileup() ) {
final ReadBackedPileup pileup = context.getBasePileup();
final ReadBackedPileup pileup = context.getBasePileup();
int mismatchQualities = 0, totalQualities = 0;
final byte refBase = ref.getBase();
for ( PileupElement p : pileup ) {
int mismatchQualities = 0, totalQualities = 0;
final byte refBase = ref.getBase();
for ( PileupElement p : pileup ) {
// check the ends of the reads to see how far they extend
furthestStopPos = Math.max(furthestStopPos, p.getRead().getAlignmentEnd());
// check the ends of the reads to see how far they extend
furthestStopPos = Math.max(furthestStopPos, p.getRead().getAlignmentEnd());
// is it a deletion or insertion?
if ( p.isDeletion() || p.isBeforeInsertion() ) {
hasIndel = true;
if ( p.isBeforeInsertion() )
hasInsertion = true;
}
// look for mismatches
else if ( lookForMismatchEntropy ) {
if ( p.getBase() != refBase )
mismatchQualities += p.getQual();
totalQualities += p.getQual();
}
// is it a deletion or insertion?
if ( p.isDeletion() || p.isBeforeInsertion() ) {
hasIndel = true;
if ( p.isBeforeInsertion() )
hasInsertion = true;
}
// make sure we're supposed to look for high entropy
if ( lookForMismatchEntropy &&
pileup.getNumberOfElements() >= minReadsAtLocus &&
(double)mismatchQualities / (double)totalQualities >= mismatchThreshold )
hasPointEvent = true;
// look for mismatches
else if ( lookForMismatchEntropy ) {
if ( p.getBase() != refBase )
mismatchQualities += p.getQual();
totalQualities += p.getQual();
}
}
// make sure we're supposed to look for high entropy
if ( lookForMismatchEntropy &&
pileup.getNumberOfElements() >= minReadsAtLocus &&
(double)mismatchQualities / (double)totalQualities >= mismatchThreshold )
hasPointEvent = true;
// return null if no event occurred
if ( !hasIndel && !hasPointEvent )
return null;

10
public/java/src/org/broadinstitute/sting/gatk/walkers/phasing/ReadBackedPhasingWalker.java
View File

@ -269,10 +269,7 @@ public class ReadBackedPhasingWalker extends RodWalker<PhasingStatsAndOutput, Ph
logger.debug("Unprocessed variant = " + VariantContextUtils.getLocation(getToolkit().getGenomeLocParser(), vc));
}
int numReads = 0;
if (context.hasBasePileup()) {
numReads = context.getBasePileup().getNumberOfElements();
}
int numReads = context.getBasePileup().getNumberOfElements();
PhasingStats addInPhaseStats = new PhasingStats(numReads, 1);
phaseStats.addIn(addInPhaseStats);
}
@ -1107,10 +1104,7 @@ public class ReadBackedPhasingWalker extends RodWalker<PhasingStatsAndOutput, Ph
this.sampleReadBases = new HashMap<String, ReadBasesAtPosition>();
if (alignment != null) {
ReadBackedPileup pileup = null;
if (alignment.hasBasePileup()) {
pileup = alignment.getBasePileup();
}
ReadBackedPileup pileup = alignment.getBasePileup();
if (pileup != null) {
// filter the read-base pileup based on min base and mapping qualities:
pileup = pileup.getBaseAndMappingFilteredPileup(MIN_BASE_QUALITY_SCORE, MIN_MAPPING_QUALITY_SCORE);

2
public/java/src/org/broadinstitute/sting/gatk/walkers/validation/GenotypeAndValidateWalker.java
View File

@ -365,7 +365,7 @@ public class GenotypeAndValidateWalker extends RodWalker<GenotypeAndValidateWalk
return counter;
// Do not operate on variants that are not covered to the optional minimum depth
if (!context.hasReads() || !context.hasBasePileup() || (minDepth > 0 && context.getBasePileup().getBases().length < minDepth)) {
if (!context.hasReads() || (minDepth > 0 && context.getBasePileup().getBases().length < minDepth)) {
counter.nUncovered = 1L;
if (vcComp.getAttribute("GV").equals("T"))
counter.nAltNotCalled = 1L;

5
public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java
View File

@ -448,10 +448,7 @@ public abstract class AbstractReadBackedPileup<RBP extends AbstractReadBackedPil
UnifiedPileupElementTracker<PE> filteredTracker = new UnifiedPileupElementTracker<PE>();
for (PE p : pileupElementTracker) {
if (p.getRead().getMappingQuality() >= minMapQ &&
(p.isDeletion() ||
((p instanceof ExtendedEventPileupElement) && ((ExtendedEventPileupElement) p).getType() == ExtendedEventPileupElement.Type.NOEVENT) ||
p.getQual() >= minBaseQ)) {
if (p.getRead().getMappingQuality() >= minMapQ && (p.isDeletion() || p.getQual() >= minBaseQ)) {
filteredTracker.add(p);
}
}

154
public/java/src/org/broadinstitute/sting/utils/pileup/ExtendedEventPileupElement.java
View File

@ -1,154 +0,0 @@
package org.broadinstitute.sting.utils.pileup;
import net.sf.samtools.SAMRecord;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
import java.util.Arrays;
/**
* In the "standard" locus traversal mode,
* the traversal is performed striclty over the reference bases. Thus, only pileups of bases (and hence local events
* such as point mutations) are "seen" at every invocation of the walker's map() function at every (genomic) locus. Deletions
* are seen on the base-by-base basis (i.e. the pileup does keep the information about the current reference base being deleted
* in some reads), but the information about the extended event (deletion length, string of all deleted bases) is not kept.
* The insertions that may be present in some reads are not seen at all in such strict reference traversal mode.
* <p/>
* By convention, any extended event (indel) is mapped onto the reference at the last base prior to the event (i.e.
* last base before the insertion or deletion). If the special "extended" traversal mode is turned on and there is
* an indel in at least one read that maps onto the reference position Z, the walker's map function will be called twice:
* first call will be performed in a "standard" mode, with a pileup of bases over the position Z, and then the additional
* call will be made at the same position with a pileup of event/noevent calls, where events are extended and contain
* full information about insertions/deletions. Then the next, "standard", call to map() will be performed at the next
* (covered) reference position. Note that if the extended event at Z was a deletion, the "standard" base pileup at
* Z+1 and following bases may still contain deleted bases. However the fully extended event call will be performed
* only once, at the position where the indel maps (starts).
* <p/>
* This class wraps an "extended" event (indel) so that in can be added to a pileup of events at a given location.
* <p/>
* Created by IntelliJ IDEA.
* User: asivache
* Date: Dec 21, 2009
* Time: 2:57:55 PM
* To change this template use File | Settings | File Templates.
*/
// Extended events are slated for removal
public class ExtendedEventPileupElement extends PileupElement {
public enum Type {
NOEVENT, DELETION, INSERTION
}
private Type type = null;
private int eventLength = -1;
private String eventBases = null; // 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
private SAMRecord read;
private int offset; // position in the read immediately BEFORE the event
// This is broken! offset is always zero because these member variables are shadowed by base class
public ExtendedEventPileupElement(GATKSAMRecord read, int offset, int eventLength, String eventBases, Type type) {
super(read, offset, type == Type.DELETION, false, false, false, false, false, null, -1); // extended events are slated for removal
this.read = read;
this.offset = offset;
this.eventLength = eventLength;
this.eventBases = eventBases;
this.type = type;
}
/**
* Quick constructor for insertions.
*
* @param read the read, in which the indel is observed
* @param offset position in the read immediately before the indel (can be -1 if read starts with an insertion)
* @param length length of the indel (number of inserted or deleted bases); length <=0 indicates that the read has no indel (NOEVENT)
* @param eventBases inserted bases. null indicates that the event is a deletion; ignored if length<=0 (noevent)
*/
public ExtendedEventPileupElement(GATKSAMRecord read, int offset, int length, byte[] eventBases) {
this(read, offset, length, new String(eventBases).toUpperCase(), Type.INSERTION);
}
/**
* Quick constructor for non indels (matches)
*
* @param read the read
* @param offset where in the read the match is
*/
public ExtendedEventPileupElement(GATKSAMRecord read, int offset) {
this(read, offset, -1, null, Type.NOEVENT);
}
/**
* Quick constructor for deletions
*
* @param read the read
* @param offset the last base before the deletion starts (left aligned deletion)
* @param length length of this deletion
*/
public ExtendedEventPileupElement(GATKSAMRecord read, int offset, int length) {
this(read, offset, length, null, Type.DELETION);
}
public boolean isDeletion() {
return type == Type.DELETION;
}
public boolean isInsertion() {
return type == Type.INSERTION;
}
public boolean isIndel() {
return isDeletion() || isInsertion();
}
public Type getType() {
return type;
}
// The offset can be negative with insertions at the start of the read, but a valid base does exist at this position with
// a valid base quality. The following code attempts to compensate for that.'
@Override
public byte getBase() {
return getBase(offset >= 0 ? offset : offset + eventLength);
}
@Override
public int getBaseIndex() {
return getBaseIndex(offset >= 0 ? offset : offset + eventLength);
}
@Override
public byte getQual() {
return getQual(offset >= 0 ? offset : offset + eventLength);
}
/**
* 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;
}
@Override
public String toString() {
char c = '.';
String fillStr = null;
if (isDeletion()) {
c = '-';
char[] filler = new char[eventLength];
Arrays.fill(filler, 'D');
fillStr = new String(filler);
} else if (isInsertion()) c = '+';
return String.format("%s @ %d = %c%s MQ%d", getRead().getReadName(), getOffset(), c, isIndel() ?
(isInsertion() ? eventBases : fillStr) : "", getMappingQual());
}
}

208
public/java/src/org/broadinstitute/sting/utils/pileup/ReadBackedExtendedEventPileup.java
View File

@ -1,208 +0,0 @@
/*
* Copyright (c) 2010, The Broad Institute
*
* 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.sting.utils.pileup;
import net.sf.samtools.SAMRecord;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.collections.Pair;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
import java.util.Collection;
import java.util.List;
/**
* A clean interface for working with extended event pileups.
*
* @author mhanna
* @version 0.1
*/
public interface ReadBackedExtendedEventPileup extends ReadBackedPileup {
/**
* Returns a new ReadBackedPileup that is free of deletion spanning reads in this pileup. Note that this
* does not copy the data, so both ReadBackedPileups should not be changed. Doesn't make an unnecessary copy
* of the pileup (just returns this) if there are no deletions in the pileup.
*
* @return
*/
public ReadBackedExtendedEventPileup getPileupWithoutDeletions();
/**
* Returns a new ReadBackedPileup where only one read from an overlapping read
* pair is retained. If the two reads in question disagree to their basecall,
* neither read is retained. If they agree on the base, the read with the higher
* quality observation is retained
*
* @return the newly filtered pileup
*/
public ReadBackedExtendedEventPileup getOverlappingFragmentFilteredPileup();
/**
* Returns a new ReadBackedPileup that is free of mapping quality zero reads in this pileup. Note that this
* does not copy the data, so both ReadBackedPileups should not be changed. Doesn't make an unnecessary copy
* of the pileup (just returns this) if there are no MQ0 reads in the pileup.
*
* @return
*/
public ReadBackedExtendedEventPileup getPileupWithoutMappingQualityZeroReads();
/**
* Gets the pileup consisting of only reads on the positive strand.
* @return A read-backed pileup consisting only of reads on the positive strand.
*/
public ReadBackedExtendedEventPileup getPositiveStrandPileup();
/**
* Gets the pileup consisting of only reads on the negative strand.
* @return A read-backed pileup consisting only of reads on the negative strand.
*/
public ReadBackedExtendedEventPileup getNegativeStrandPileup();
/**
* Gets a pileup consisting of all those elements passed by a given filter.
* @param filter Filter to use when testing for elements.
* @return a pileup without the given filtered elements.
*/
public ReadBackedExtendedEventPileup getFilteredPileup(PileupElementFilter filter);
/** Returns subset of this pileup that contains only bases with quality >= minBaseQ, coming from
* reads with mapping qualities >= minMapQ. This method allocates and returns a new instance of ReadBackedPileup.
* @param minBaseQ
* @param minMapQ
* @return
*/
public ReadBackedExtendedEventPileup getBaseAndMappingFilteredPileup( int minBaseQ, int minMapQ );
/** Returns subset of this pileup that contains only bases with quality >= minBaseQ.
* This method allocates and returns a new instance of ReadBackedPileup.
* @param minBaseQ
* @return
*/
public ReadBackedExtendedEventPileup getBaseFilteredPileup( int minBaseQ );
/** Returns subset of this pileup that contains only bases coming from reads with mapping quality >= minMapQ.
* This method allocates and returns a new instance of ReadBackedPileup.
* @param minMapQ
* @return
*/
public ReadBackedExtendedEventPileup getMappingFilteredPileup( int minMapQ );
/**
* Returns a pileup randomly downsampled to the desiredCoverage.
*
* @param desiredCoverage
* @return
*/
public ReadBackedExtendedEventPileup getDownsampledPileup(int desiredCoverage);
/**
* Gets a list of all the samples stored in this pileup.
* @return List of samples in this pileup.
*/
public Collection<String> getSamples();
public Iterable<ExtendedEventPileupElement> toExtendedIterable();
/**
* Returns the number of deletion events in this pileup
*
* @return
*/
public int getNumberOfDeletions();
/**
* Returns the number of insertion events in this pileup
*
* @return
*/
public int getNumberOfInsertions();
/** Returns the length of the longest deletion observed at the site this
* pileup is associated with (NOTE: by convention, both insertions and deletions
* are associated with genomic location immediately before the actual event). If
* there are no deletions at the site, returns 0.
* @return
*/
public int getMaxDeletionLength();
/**
* Returns the number of mapping quality zero reads in this pileup.
* @return
*/
public int getNumberOfMappingQualityZeroReads();
/**
* @return the number of elements in this pileup
*/
public int getNumberOfElements();
/**
* @return the location of this pileup
*/
public GenomeLoc getLocation();
/**
* Returns a list of the reads in this pileup. Note this call costs O(n) and allocates fresh lists each time
* @return
*/
public List<GATKSAMRecord> getReads();
/**
* Returns a list of the offsets in this pileup. Note this call costs O(n) and allocates fresh lists each time
* @return
*/
public List<Integer> getOffsets();
/**
* Returns an array of the events in this pileup ('I', 'D', or '.'). Note this call costs O(n) and allocates fresh array each time
* @return
*/
public byte[] getEvents();
/** A shortcut for getEventStringsWithCounts(null);
*
* @return
*/
public List<Pair<String,Integer>> getEventStringsWithCounts();
/** Returns String representation of all distinct extended events (indels) at the site along with
* observation counts (numbers of reads) for each distinct event. If refBases is null, a simple string representation for
* deletions will be generated as "<length>D" (i.e. "5D"); if the reference bases are provided, the actual
* deleted sequence will be used in the string representation (e.g. "-AAC").
* @param refBases reference bases, starting with the current locus (i.e. the one immediately before the indel), and
* extending far enough to accomodate the longest deletion (i.e. size of refBases must be at least 1+<length of longest deletion>)
* @return list of distinct events; first element of a pair is a string representation of the event, second element
* gives the number of reads, in which that event was observed
*/
public List<Pair<String,Integer>> getEventStringsWithCounts(byte[] refBases);
public String getShortPileupString();
/**
* Get an array of the mapping qualities
* @return
*/
public byte[] getMappingQuals();
}

248
public/java/src/org/broadinstitute/sting/utils/pileup/ReadBackedExtendedEventPileupImpl.java
View File

@ -1,248 +0,0 @@
/*
* Copyright (c) 2010, The Broad Institute
*
* 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.sting.utils.pileup;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.collections.Pair;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
import java.util.*;
public class ReadBackedExtendedEventPileupImpl extends AbstractReadBackedPileup<ReadBackedExtendedEventPileupImpl, ExtendedEventPileupElement> implements ReadBackedExtendedEventPileup {
private int nInsertions;
private int maxDeletionLength; // cached value of the length of the longest deletion observed at the site
public ReadBackedExtendedEventPileupImpl(GenomeLoc loc, List<ExtendedEventPileupElement> pileupElements) {
super(loc, pileupElements);
}
public ReadBackedExtendedEventPileupImpl(GenomeLoc loc, PileupElementTracker<ExtendedEventPileupElement> tracker) {
super(loc, tracker);
}
/**
* Optimization of above constructor where all of the cached data is provided
*
* @param loc
* @param pileup
*/
public ReadBackedExtendedEventPileupImpl(GenomeLoc loc, List<ExtendedEventPileupElement> pileup, int size,
int maxDeletionLength, int nInsertions, int nDeletions, int nMQ0Reads) {
super(loc, pileup, size, nDeletions, nMQ0Reads);
this.maxDeletionLength = maxDeletionLength;
this.nInsertions = nInsertions;
}
// this is the good new one
public ReadBackedExtendedEventPileupImpl(GenomeLoc loc, Map<String, ? extends ReadBackedExtendedEventPileupImpl> pileupElementsBySample) {
super(loc, pileupElementsBySample);
}
/**
* Calculate cached sizes, nDeletion, and base counts for the pileup. This calculation is done upfront,
* so you pay the cost at the start, but it's more efficient to do this rather than pay the cost of calling
* sizes, nDeletion, etc. over and over potentially.
*/
@Override
protected void calculateCachedData() {
super.calculateCachedData();
nInsertions = 0;
nMQ0Reads = 0;
for (ExtendedEventPileupElement p : this.toExtendedIterable()) {
if (p.isDeletion()) {
maxDeletionLength = Math.max(maxDeletionLength, p.getEventLength());
} else {
if (p.isInsertion()) nInsertions++;
}
}
}
@Override
protected void addPileupToCumulativeStats(AbstractReadBackedPileup<ReadBackedExtendedEventPileupImpl, ExtendedEventPileupElement> pileup) {
super.addPileupToCumulativeStats(pileup);
ReadBackedExtendedEventPileup extendedEventPileup = ((ReadBackedExtendedEventPileup) pileup);
this.nInsertions += extendedEventPileup.getNumberOfInsertions();
this.maxDeletionLength += extendedEventPileup.getMaxDeletionLength();
}
@Override
protected ReadBackedExtendedEventPileupImpl createNewPileup(GenomeLoc loc, PileupElementTracker<ExtendedEventPileupElement> tracker) {
return new ReadBackedExtendedEventPileupImpl(loc, tracker);
}
@Override
protected ExtendedEventPileupElement createNewPileupElement(final GATKSAMRecord read, final int offset, final boolean isDeletion, final boolean isBeforeDeletion, final boolean isAfterDeletion, final boolean isBeforeInsertion, final boolean isAfterInsertion, final boolean isNextToSoftClip) {
throw new UnsupportedOperationException("Not enough information provided to create a new pileup element");
}
@Override
protected ExtendedEventPileupElement createNewPileupElement(final GATKSAMRecord read, final int offset, final boolean isDeletion, final boolean isBeforeDeletion, final boolean isAfterDeletion, final boolean isBeforeInsertion, final boolean isAfterInsertion, final boolean isNextToSoftClip, final String nextEventBases, final int nextEventLength ) {
throw new UnsupportedOperationException("Not enough information provided to create a new pileup element");
}
/**
* Returns the number of insertion events in this pileup
*
* @return
*/
@Override
public int getNumberOfInsertions() {
return nInsertions;
}
/**
* Returns the length of the longest deletion observed at the site this
* pileup is associated with (NOTE: by convention, both insertions and deletions
* are associated with genomic location immediately before the actual event). If
* there are no deletions at the site, returns 0.
*
* @return
*/
@Override
public int getMaxDeletionLength() {
return maxDeletionLength;
}
public Iterable<ExtendedEventPileupElement> toExtendedIterable() {
return new Iterable<ExtendedEventPileupElement>() {
public Iterator<ExtendedEventPileupElement> iterator() {
return pileupElementTracker.iterator();
}
};
}
/**
* Returns an array of the events in this pileup ('I', 'D', or '.'). Note this call costs O(n) and allocates fresh array each time
*
* @return
*/
@Override
public byte[] getEvents() {
byte[] v = new byte[getNumberOfElements()];
int i = 0;
for (ExtendedEventPileupElement e : this.toExtendedIterable()) {
switch (e.getType()) {
case INSERTION:
v[i] = 'I';
break;
case DELETION:
v[i] = 'D';
break;
case NOEVENT:
v[i] = '.';
break;
default:
throw new ReviewedStingException("Unknown event type encountered: " + e.getType());
}
i++;
}
return v;
}
/**
* A shortcut for getEventStringsWithCounts(null);
*
* @return
*/
@Override
public List<Pair<String, Integer>> getEventStringsWithCounts() {
return getEventStringsWithCounts(null);
}
@Override
public String getShortPileupString() {
// In the pileup format, each extended event line has genomic position (chromosome name and offset),
// reference "base" (always set to "E" for E(xtended)), then 'I','D', or '.' symbol for each read representing
// insertion, deletion or no-event, respectively.
return String.format("%s %s E %s",
getLocation().getContig(), getLocation().getStart(), // chromosome name and coordinate
new String(getEvents()));
}
/**
* Returns String representation of all distinct extended events (indels) at the site along with
* observation counts (numbers of reads) for each distinct event. If refBases is null, a simple string representation for
* deletions will be generated as "<length>D" (i.e. "5D"); if the reference bases are provided, the actual
* deleted sequence will be used in the string representation (e.g. "-AAC").
*
* @param refBases reference bases, starting with the current locus (i.e. the one immediately before the indel), and
* extending far enough to accomodate the longest deletion (i.e. size of refBases must be at least 1+<length of longest deletion>)
* @return list of distinct events; first element of a pair is a string representation of the event, second element
* gives the number of reads, in which that event was observed
*/
@Override
public List<Pair<String, Integer>> getEventStringsWithCounts(byte[] refBases) {
Map<String, Integer> events = new HashMap<String, Integer>();
for (ExtendedEventPileupElement e : this.toExtendedIterable()) {
Integer cnt;
String indel = null;
switch (e.getType()) {
case INSERTION:
indel = "+" + e.getEventBases();
break;
case DELETION:
indel = getDeletionString(e.getEventLength(), refBases);
break;
case NOEVENT:
continue;
default:
throw new ReviewedStingException("Unknown event type encountered: " + e.getType());
}
cnt = events.get(indel);
if (cnt == null) events.put(indel, 1);
else events.put(indel, cnt.intValue() + 1);
}
List<Pair<String, Integer>> eventList = new ArrayList<Pair<String, Integer>>(events.size());
for (Map.Entry<String, Integer> m : events.entrySet()) {
eventList.add(new Pair<String, Integer>(m.getKey(), m.getValue()));
}
return eventList;
}
/**
* Builds string representation of the deletion event. If refBases is null, the representation will be
* "<length>D" (e.g. "5D"); if the reference bases are provided, a verbose representation (e.g. "-AAC")
* will be generated. NOTE: refBases must start with the base prior to the actual deletion (i.e. deleted
* base(s) are refBase[1], refBase[2], ...), and the length of the passed array must be sufficient to accomodate the
* deletion length (i.e. size of refBase must be at least length+1).
*
* @param length
* @param refBases
* @return
*/
private String getDeletionString(int length, byte[] refBases) {
if (refBases == null) {
return Integer.toString(length) + "D"; // if we do not have reference bases, we can only report something like "5D"
} else {
return "-" + new String(refBases, 1, length).toUpperCase();
}
}
}

5
public/java/test/org/broadinstitute/sting/gatk/iterators/LocusIteratorByStateUnitTest.java
View File

@ -17,8 +17,6 @@ import org.broadinstitute.sting.gatk.datasources.reads.SAMReaderID;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.baq.BAQ;
import org.broadinstitute.sting.utils.pileup.ReadBackedExtendedEventPileup;
import org.broadinstitute.sting.utils.classloader.JVMUtils;
import org.broadinstitute.sting.utils.sam.ArtificialSAMUtils;
import org.testng.annotations.BeforeClass;
import org.testng.annotations.Test;
@ -74,9 +72,6 @@ public class LocusIteratorByStateUnitTest extends BaseTest {
boolean foundIndel = false;
while (li.hasNext()) {
AlignmentContext context = li.next();
if(!context.hasBasePileup())
continue;
ReadBackedPileup pileup = context.getBasePileup().getBaseFilteredPileup(10);
for (PileupElement p : pileup) {
if (p.isBeforeInsertion()) {