From 59df32977609f0fa8ed874482954356741411b39 Mon Sep 17 00:00:00 2001 From: Mark DePristo Date: Tue, 5 Feb 2013 10:11:05 -0500 Subject: [PATCH] Fast path for biallelic variants in IndependentAllelesDiploidExactAFCalc -- If the VariantContext is a bi-allelic variant already, don't split up the VC (it doesn't do anything) and then combine it back together. This saves us a lot of work on average -- Be more protective of calls to AFCalc with a VariantContext that might only have ref allele, throwing an exception --- .../gatk/walkers/genotyper/afcalc/AFCalc.java | 11 +++--- .../IndependentAllelesDiploidExactAFCalc.java | 34 ++++++++++++++----- 2 files changed, 32 insertions(+), 13 deletions(-) diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalc.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalc.java index 5334847ae..642a2b32f 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalc.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalc.java @@ -113,12 +113,14 @@ public abstract class AFCalc implements Cloneable { /** * Compute the probability of the alleles segregating given the genotype likelihoods of the samples in vc * - * @param vc the VariantContext holding the alleles and sample information + * @param vc the VariantContext holding the alleles and sample information. The VariantContext + * must have at least 1 alternative allele * @param log10AlleleFrequencyPriors a prior vector nSamples x 2 in length indicating the Pr(AF = i) * @return result (for programming convenience) */ public AFCalcResult getLog10PNonRef(final VariantContext vc, final double[] log10AlleleFrequencyPriors) { if ( vc == null ) throw new IllegalArgumentException("VariantContext cannot be null"); + if ( vc.getNAlleles() == 1 ) throw new IllegalArgumentException("VariantContext has only a single reference allele, but getLog10PNonRef requires at least one at all " + vc); if ( log10AlleleFrequencyPriors == null ) throw new IllegalArgumentException("priors vector cannot be null"); if ( stateTracker == null ) throw new IllegalArgumentException("Results object cannot be null"); @@ -170,18 +172,19 @@ public abstract class AFCalc implements Cloneable { * @param vc the initial VC provided by the caller to this AFcalculation * @return a potentially simpler VC that's more tractable to genotype */ - @Requires("vc != null") + @Requires({"vc != null", "vc.getNAlleles() > 1"}) @Ensures("result != null") protected abstract VariantContext reduceScope(final VariantContext vc); /** * Actually carry out the log10PNonRef calculation on vc, storing results in results * - * @param vc variant context with alleles and genotype likelihoods + * @param vc variant context with alleles and genotype likelihoods, + * must have at least one alt allele * @param log10AlleleFrequencyPriors priors * @return a AFCalcResult object describing the results of this calculation */ - @Requires({"vc != null", "log10AlleleFrequencyPriors != null"}) + @Requires({"vc != null", "log10AlleleFrequencyPriors != null", "vc.getNAlleles() > 1"}) protected abstract AFCalcResult computeLog10PNonRef(final VariantContext vc, final double[] log10AlleleFrequencyPriors); diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/IndependentAllelesDiploidExactAFCalc.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/IndependentAllelesDiploidExactAFCalc.java index 56d50ceba..af5c79230 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/IndependentAllelesDiploidExactAFCalc.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/IndependentAllelesDiploidExactAFCalc.java @@ -156,16 +156,25 @@ import java.util.*; public AFCalcResult computeLog10PNonRef(final VariantContext vc, final double[] log10AlleleFrequencyPriors) { final List independentResultTrackers = computeAlleleIndependentExact(vc, log10AlleleFrequencyPriors); - final List withMultiAllelicPriors = applyMultiAllelicPriors(independentResultTrackers); - return combineIndependentPNonRefs(vc, withMultiAllelicPriors); - } + if ( independentResultTrackers.size() == 0 ) + throw new IllegalStateException("Independent alleles model returned an empty list of results at VC " + vc); + + if ( independentResultTrackers.size() == 1 ) { + // fast path for the very common bi-allelic use case + return independentResultTrackers.get(0); + } else { + // we are a multi-allelic, so we need to actually combine the results + final List withMultiAllelicPriors = applyMultiAllelicPriors(independentResultTrackers); + return combineIndependentPNonRefs(vc, withMultiAllelicPriors); + } + } /** * Compute the conditional exact AFCalcResult for each allele in vc independently, returning * the result of each, in order of the alt alleles in VC * - * @param vc the VariantContext we want to analyze + * @param vc the VariantContext we want to analyze, with at least 1 alt allele * @param log10AlleleFrequencyPriors the priors * @return a list of the AFCalcResults for each bi-allelic sub context of vc */ @@ -208,13 +217,20 @@ import java.util.*; @Ensures("result.size() == vc.getNAlleles() - 1") protected final List makeAlleleConditionalContexts(final VariantContext vc) { final int nAltAlleles = vc.getNAlleles() - 1; - final List vcs = new LinkedList(); - for ( int altI = 0; altI < nAltAlleles; altI++ ) { - vcs.add(biallelicCombinedGLs(vc, altI + 1)); + if ( nAltAlleles == 1 ) { + // fast path for bi-allelic case. + return Collections.singletonList(vc); + } else { + // go through the work of ripping up the VC into its biallelic components + final List vcs = new LinkedList(); + + for ( int altI = 0; altI < nAltAlleles; altI++ ) { + vcs.add(biallelicCombinedGLs(vc, altI + 1)); + } + + return vcs; } - - return vcs; } /**