From e0ab4e4b306d2c8227f4986bd410fc9cdcb0311a Mon Sep 17 00:00:00 2001 From: Guillermo del Angel Date: Wed, 28 Mar 2012 21:01:31 -0400 Subject: [PATCH] Refactoring so that ConsensusAlleleCounter can use regular pileups and can operate correctly. This involved adding utility functions to ReadBackedPileup to count # of insertions/deletions right after current position. Added unit test for IndelGenotypeLikelihoods, esp. ConsensusAlleleCounter logic --- .../genotyper/ConsensusAlleleCounter.java | 56 ++++++++++++++----- .../pileup/AbstractReadBackedPileup.java | 20 +++++++ .../sting/utils/pileup/ReadBackedPileup.java | 14 +++++ 3 files changed, 75 insertions(+), 15 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ConsensusAlleleCounter.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ConsensusAlleleCounter.java index 4cf6586a6..51d3fb92b 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ConsensusAlleleCounter.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ConsensusAlleleCounter.java @@ -109,15 +109,10 @@ public class ConsensusAlleleCounter { delCount += indelPileup.getNumberOfDeletions(); } else { - // todo - this should be version to be used when extended events are removed - // todo - maybe we should create utility functions in ReadBackedPileup definition to do the equivalent thing? - for (PileupElement p: context.getBasePileup()) { - if (p.isBeforeDeletion()) - delCount++; - else if (p.isBeforeInsertion()) - insCount++; - } - } + final ReadBackedPileup indelPileup = context.getBasePileup(); + insCount += indelPileup.getNumberOfInsertionsAfterThisElement(); + delCount += indelPileup.getNumberOfDeletionsAfterThisElement(); + } } if (insCount < minIndelCountForGenotyping && delCount < minIndelCountForGenotyping) @@ -127,10 +122,20 @@ public class ConsensusAlleleCounter { // todo -- warning, can be duplicating expensive partition here AlignmentContext context = AlignmentContextUtils.stratify(sample.getValue(), contextType); - final ReadBackedExtendedEventPileup indelPileup = context.getExtendedEventPileup(); + final ReadBackedPileup indelPileup; - final int nIndelReads = indelPileup.getNumberOfInsertions() + indelPileup.getNumberOfDeletions(); - final int nReadsOverall = indelPileup.getNumberOfElements(); + final int nIndelReads, nReadsOverall; + + if (context.hasExtendedEventPileup()) { + indelPileup = context.getExtendedEventPileup(); + nIndelReads = ((ReadBackedExtendedEventPileup)indelPileup).getNumberOfInsertions() + indelPileup.getNumberOfDeletions(); + nReadsOverall = indelPileup.getNumberOfElements(); + } + else { + indelPileup = context.getBasePileup(); + nIndelReads = indelPileup.getNumberOfInsertionsAfterThisElement() + indelPileup.getNumberOfDeletionsAfterThisElement(); + nReadsOverall = indelPileup.getNumberOfElements(); + } if ( nIndelReads == 0 || (nIndelReads / (1.0 * nReadsOverall)) < minFractionInOneSample) { // if ( nIndelReads > 0 ) // logger.info("Skipping sample " + sample.getKey() + " with nIndelReads " + nIndelReads + " nReads " + nReadsOverall); @@ -139,7 +144,8 @@ public class ConsensusAlleleCounter { // logger.info("### Keeping sample " + sample.getKey() + " with nIndelReads " + nIndelReads + " nReads " + nReadsOverall); } - for (ExtendedEventPileupElement p : indelPileup.toExtendedIterable()) { + + for (PileupElement p : indelPileup) { final GATKSAMRecord read = ReadClipper.hardClipAdaptorSequence(p.getRead()); if (read == null) continue; @@ -154,7 +160,8 @@ public class ConsensusAlleleCounter { } */ String indelString = p.getEventBases(); - if (p.isInsertion()) { + + if (isInsertion(p)) { boolean foundKey = false; // copy of hashmap into temp arrayList ArrayList> cList = new ArrayList>(); @@ -222,7 +229,7 @@ public class ConsensusAlleleCounter { } } - else if (p.isDeletion()) { + else if (isDeletion(p)) { indelString = String.format("D%d",p.getEventLength()); int cnt = consensusIndelStrings.containsKey(indelString)? consensusIndelStrings.get(indelString):0; consensusIndelStrings.put(indelString,cnt+1); @@ -234,6 +241,25 @@ public class ConsensusAlleleCounter { return consensusIndelStrings; } + + // todo - helper routines to check for extended pileup elements, to remove when extended events are removed + private static final boolean isInsertion(final PileupElement p) { + if (p instanceof ExtendedEventPileupElement) + return ((ExtendedEventPileupElement) p).isInsertion(); + else + return p.isBeforeInsertion(); + + } + + private static boolean isDeletion(final PileupElement p) { + if (p instanceof ExtendedEventPileupElement) + return p.isDeletion(); + else + return p.isBeforeDeletion(); + + } + + private List consensusCountsToAlleles(final ReferenceContext ref, final Map consensusIndelStrings) { final GenomeLoc loc = ref.getLocus(); 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 c8f00778f..5a7e0f1c5 100644 --- a/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java +++ b/public/java/src/org/broadinstitute/sting/utils/pileup/AbstractReadBackedPileup.java @@ -873,6 +873,26 @@ public abstract class AbstractReadBackedPileup, HasGenomeLoca */ public int getNumberOfDeletions(); + /** + * Simple useful routine to count the number of deletion bases in at the next position this pileup + * + * @return + */ + public int getNumberOfDeletionsAfterThisElement(); + + /** + * Simple useful routine to count the number of insertions right after this pileup + * + * @return + */ + public int getNumberOfInsertionsAfterThisElement(); + public int getNumberOfMappingQualityZeroReads(); /**