diff --git a/java/src/org/broadinstitute/sting/gatk/contexts/StratifiedAlignmentContext.java b/java/src/org/broadinstitute/sting/gatk/contexts/StratifiedAlignmentContext.java index 5698b90fc..cc51a3dbf 100755 --- a/java/src/org/broadinstitute/sting/gatk/contexts/StratifiedAlignmentContext.java +++ b/java/src/org/broadinstitute/sting/gatk/contexts/StratifiedAlignmentContext.java @@ -42,63 +42,49 @@ import java.util.Map; */ public class StratifiedAlignmentContext { - public enum StratifiedContextType { OVERALL, FORWARD, REVERSE } + // Definitions: + // COMPLETE = full alignment context + // MQ0FREE = full context without MQ0 reads + // FORWARD = reads on forward strand (*no* MQ0 reads) + // REVERSE = reads on forward strand (*no* MQ0 reads) + // + public enum StratifiedContextType { COMPLETE, MQ0FREE, FORWARD, REVERSE } - private AlignmentContext overall = null; - private AlignmentContext forward = null; - private AlignmentContext reverse = null; private GenomeLoc loc; - - private ArrayList allReads = new ArrayList(); - private ArrayList forwardReads = new ArrayList(); - private ArrayList reverseReads = new ArrayList(); - - private ArrayList allOffsets = new ArrayList(); - private ArrayList forwardOffsets = new ArrayList(); - private ArrayList reverseOffsets = new ArrayList(); + private AlignmentContext[] contexts = new AlignmentContext[StratifiedContextType.values().length]; + private ArrayList[] reads = new ArrayList[StratifiedContextType.values().length]; + private ArrayList[] offsets = new ArrayList[StratifiedContextType.values().length]; public StratifiedAlignmentContext(GenomeLoc loc) { this.loc = loc; + for ( int i = 0; i < StratifiedContextType.values().length; i++) { + reads[i] = new ArrayList(); + offsets[i] = new ArrayList(); + } } public AlignmentContext getContext(StratifiedContextType context) { - switch ( context ) { - case OVERALL: return getOverallContext(); - case FORWARD: return getForwardContext(); - case REVERSE: return getReverseContext(); - } - return null; - } - - private AlignmentContext getOverallContext() { - if ( overall == null ) - overall = new AlignmentContext(loc, new ReadBackedPileup(loc, allReads, allOffsets)); - return overall; - } - - private AlignmentContext getForwardContext() { - if ( forward == null ) - forward = new AlignmentContext(loc, new ReadBackedPileup(loc, forwardReads, forwardOffsets)); - return forward; - } - - private AlignmentContext getReverseContext() { - if ( reverse == null ) - reverse = new AlignmentContext(loc, new ReadBackedPileup(loc, reverseReads, reverseOffsets)); - return reverse; + int index = context.ordinal(); + if ( contexts[index] == null ) + contexts[index] = new AlignmentContext(loc, new ReadBackedPileup(loc, reads[index], offsets[index])); + return contexts[index]; } public void add(SAMRecord read, int offset) { - if ( read.getReadNegativeStrandFlag() ) { - reverseReads.add(read); - reverseOffsets.add(offset); - } else { - forwardReads.add(read); - forwardOffsets.add(offset); + if ( read.getMappingQuality() > 0 ) { + reads[StratifiedContextType.MQ0FREE.ordinal()].add(read); + offsets[StratifiedContextType.MQ0FREE.ordinal()].add(offset); + if ( read.getReadNegativeStrandFlag() ) { + reads[StratifiedContextType.REVERSE.ordinal()].add(read); + offsets[StratifiedContextType.REVERSE.ordinal()].add(offset); + } else { + reads[StratifiedContextType.FORWARD.ordinal()].add(read); + offsets[StratifiedContextType.FORWARD.ordinal()].add(offset); + } } - allReads.add(read); - allOffsets.add(offset); + reads[StratifiedContextType.COMPLETE.ordinal()].add(read); + offsets[StratifiedContextType.COMPLETE.ordinal()].add(offset); } /** diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypeCalculationModel.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypeCalculationModel.java index ed8ce27a1..298a0fea6 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypeCalculationModel.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/DiploidGenotypeCalculationModel.java @@ -83,7 +83,7 @@ public class DiploidGenotypeCalculationModel extends JointEstimateGenotypeCalcul GenotypeCall call = GenotypeWriterFactory.createSupportedGenotypeCall(OUTPUT_FORMAT, ref, loc); if ( call instanceof ReadBacked ) { - ReadBackedPileup pileup = contexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.OVERALL).getPileup(); + ReadBackedPileup pileup = contexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.MQ0FREE).getPileup(); ((ReadBacked)call).setPileup(pileup); } if ( call instanceof SampleBacked ) { diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/EMGenotypeCalculationModel.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/EMGenotypeCalculationModel.java index 4b987a099..d513b1b20 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/EMGenotypeCalculationModel.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/EMGenotypeCalculationModel.java @@ -22,7 +22,7 @@ public abstract class EMGenotypeCalculationModel extends GenotypeCalculationMode public Pair> calculateGenotype(RefMetaDataTracker tracker, char ref, GenomeLoc loc, Map contexts, DiploidGenotypePriors priors) { // run the EM calculation - EMOutput overall = runEM(ref, contexts, priors, StratifiedAlignmentContext.StratifiedContextType.OVERALL); + EMOutput overall = runEM(ref, contexts, priors, StratifiedAlignmentContext.StratifiedContextType.MQ0FREE); double PofD = Math.pow(10, overall.getPofD()); double PofNull = Math.pow(10, overall.getPofNull()); @@ -93,11 +93,11 @@ public abstract class EMGenotypeCalculationModel extends GenotypeCalculationMode for ( String sample : GLs.keySet() ) { // create the call - AlignmentContext context = contexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.OVERALL); + AlignmentContext context = contexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.MQ0FREE); GenotypeCall call = GenotypeWriterFactory.createSupportedGenotypeCall(OUTPUT_FORMAT, ref, context.getLocation()); if ( call instanceof ReadBacked ) { - ReadBackedPileup pileup = contexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.OVERALL).getPileup(); + ReadBackedPileup pileup = contexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.MQ0FREE).getPileup(); ((ReadBacked)call).setPileup(pileup); } if ( call instanceof SampleBacked ) { diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeCalculationModel.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeCalculationModel.java index 9d594b755..902b79f72 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeCalculationModel.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeCalculationModel.java @@ -10,8 +10,6 @@ import org.apache.log4j.Logger; import java.io.*; import java.util.*; -import net.sf.samtools.SAMRecord; -import net.sf.samtools.SAMReadGroupRecord; /** * The model representing how we calculate a genotype given the priors and a pile diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/JointEstimateGenotypeCalculationModel.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/JointEstimateGenotypeCalculationModel.java index edd87adba..e252b4834 100644 --- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/JointEstimateGenotypeCalculationModel.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/JointEstimateGenotypeCalculationModel.java @@ -55,8 +55,8 @@ public abstract class JointEstimateGenotypeCalculationModel extends GenotypeCalc initializeAlleleFrequencies(frequencyEstimationPoints); - initialize(ref, contexts, StratifiedAlignmentContext.StratifiedContextType.OVERALL); - calculateAlleleFrequencyPosteriors(ref, frequencyEstimationPoints, contexts, StratifiedAlignmentContext.StratifiedContextType.OVERALL); + initialize(ref, contexts, StratifiedAlignmentContext.StratifiedContextType.MQ0FREE); + calculateAlleleFrequencyPosteriors(ref, frequencyEstimationPoints, contexts, StratifiedAlignmentContext.StratifiedContextType.MQ0FREE); calculatePofFs(ref, frequencyEstimationPoints); // print out stats if we have a writer @@ -74,7 +74,7 @@ public abstract class JointEstimateGenotypeCalculationModel extends GenotypeCalc int[] qualCounts = new int[4]; for ( String sample : contexts.keySet() ) { - AlignmentContext context = contexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.OVERALL); + AlignmentContext context = contexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.MQ0FREE); // calculate the sum of quality scores for each base ReadBackedPileup pileup = context.getPileup(); diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/PointEstimateGenotypeCalculationModel.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/PointEstimateGenotypeCalculationModel.java index 7470fa65c..3500468b0 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/PointEstimateGenotypeCalculationModel.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/PointEstimateGenotypeCalculationModel.java @@ -40,7 +40,7 @@ public class PointEstimateGenotypeCalculationModel extends EMGenotypeCalculation return null; // get the genotype likelihoods - Pair discoveryGL = getSingleSampleLikelihoods(sampleContext, priors, StratifiedAlignmentContext.StratifiedContextType.OVERALL); + Pair discoveryGL = getSingleSampleLikelihoods(sampleContext, priors, StratifiedAlignmentContext.StratifiedContextType.MQ0FREE); // find the index of the best genotype double[] posteriors = discoveryGL.second.getNormalizedPosteriors(); diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java index 7e9d5b11e..3f4a02949 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyper.java @@ -190,43 +190,40 @@ public class UnifiedGenotyper extends LocusWalker> map(RefMetaDataTracker tracker, ReferenceContext refContext, AlignmentContext fullContext) { + public Pair> map(RefMetaDataTracker tracker, ReferenceContext refContext, AlignmentContext context) { char ref = Character.toUpperCase(refContext.getBase()); if ( !BaseUtils.isRegularBase(ref) ) return null; - // remove mapping quality zero reads - AlignmentContext MQ0freeContext = filterAlignmentContext(fullContext); - // an optimization to speed things up when there is no coverage or when overly covered - if ( MQ0freeContext.getPileup().size() == 0 || - (UAC.MAX_READS_IN_PILEUP > 0 && MQ0freeContext.getPileup().size() > UAC.MAX_READS_IN_PILEUP) ) + if ( context.getPileup().size() == 0 || + (UAC.MAX_READS_IN_PILEUP > 0 && context.getPileup().size() > UAC.MAX_READS_IN_PILEUP) ) return null; // are there too many deletions in the pileup? - ReadBackedPileup pileup = MQ0freeContext.getPileup(); + ReadBackedPileup pileup = context.getPileup(); if ( isValidDeletionFraction(UAC.MAX_DELETION_FRACTION) && (double)pileup.getNumberOfDeletions() / (double)pileup.size() > UAC.MAX_DELETION_FRACTION ) return null; // stratify the AlignmentContext and cut by sample // Note that for testing purposes, we may want to throw multi-samples at pooled mode - Map stratifiedContexts = StratifiedAlignmentContext.splitContextBySample(MQ0freeContext, UAC.ASSUME_SINGLE_SAMPLE, (UAC.genotypeModel == GenotypeCalculationModel.Model.POOLED ? PooledCalculationModel.POOL_SAMPLE_NAME : null)); + Map stratifiedContexts = StratifiedAlignmentContext.splitContextBySample(context, UAC.ASSUME_SINGLE_SAMPLE, (UAC.genotypeModel == GenotypeCalculationModel.Model.POOLED ? PooledCalculationModel.POOL_SAMPLE_NAME : null)); if ( stratifiedContexts == null ) return null; DiploidGenotypePriors priors = new DiploidGenotypePriors(ref, UAC.heterozygosity, DiploidGenotypePriors.PROB_OF_TRISTATE_GENOTYPE); - Pair> call = gcm.calculateGenotype(tracker, ref, fullContext.getLocation(), stratifiedContexts, priors); + Pair> call = gcm.calculateGenotype(tracker, ref, context.getLocation(), stratifiedContexts, priors); // annotate the call, if possible if ( call != null && call.first != null && call.first instanceof ArbitraryFieldsBacked ) { Map annotations; if ( UAC.ALL_ANNOTATIONS ) - annotations = VariantAnnotator.getAllAnnotations(refContext, fullContext, call.first); + annotations = VariantAnnotator.getAllAnnotations(refContext, context, call.first); else - annotations = VariantAnnotator.getAnnotations(refContext, fullContext, call.first); + annotations = VariantAnnotator.getAnnotations(refContext, context, call.first); ((ArbitraryFieldsBacked)call.first).setFields(annotations); }