From b4749757f81f099dd1086f867de02fb1ebd50f5c Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Mon, 12 Mar 2012 01:07:07 -0400 Subject: [PATCH] Fixes for SLOD: 1) didn't work properly for multi-allelics (randomly chose an allele, possibly one that wasn't genotyped in the full context); 2) in cases when there were more alt alleles than the max allowed and the user is calculating SB, we would recompute the best alt alleles(s); 3) for some reason, we were recomputing the LOD for the full context when we'd already done that. Given that this passes integration tests on my end, this should be the last commit before the release. --- .../GenotypeLikelihoodsCalculationModel.java | 33 ++++++++++--------- ...elGenotypeLikelihoodsCalculationModel.java | 15 +++++---- ...NPGenotypeLikelihoodsCalculationModel.java | 6 ++-- .../genotyper/UnifiedGenotyperEngine.java | 14 ++++---- 4 files changed, 34 insertions(+), 34 deletions(-) diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeLikelihoodsCalculationModel.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeLikelihoodsCalculationModel.java index ace780dd0f..fb2428258a 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeLikelihoodsCalculationModel.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/GenotypeLikelihoodsCalculationModel.java @@ -38,6 +38,7 @@ import org.broadinstitute.sting.utils.variantcontext.Allele; import org.broadinstitute.sting.utils.variantcontext.VariantContext; +import java.util.List; import java.util.Map; @@ -76,24 +77,24 @@ protected GenotypeLikelihoodsCalculationModel(UnifiedArgumentCollection UAC, Log /** * Can be overridden by concrete subclasses * - * @param tracker rod data - * @param ref reference context - * @param contexts stratified alignment contexts - * @param contextType stratified context type - * @param priors priors to use for GLs - * @param alternateAlleleToUse the alternate allele to use, null if not set - * @param useBAQedPileup should we use the BAQed pileup or the raw one? - * @param locParser Genome Loc Parser + * @param tracker rod data + * @param ref reference context + * @param contexts stratified alignment contexts + * @param contextType stratified context type + * @param priors priors to use for GLs + * @param alternateAllelesToUse the alternate allele to use, null if not set + * @param useBAQedPileup should we use the BAQed pileup or the raw one? + * @param locParser Genome Loc Parser * @return variant context where genotypes are no-called but with GLs */ - public abstract VariantContext getLikelihoods(RefMetaDataTracker tracker, - ReferenceContext ref, - Map contexts, - AlignmentContextUtils.ReadOrientation contextType, - GenotypePriors priors, - Allele alternateAlleleToUse, - boolean useBAQedPileup, - GenomeLocParser locParser); + public abstract VariantContext getLikelihoods(final RefMetaDataTracker tracker, + final ReferenceContext ref, + final Map contexts, + final AlignmentContextUtils.ReadOrientation contextType, + final GenotypePriors priors, + final List alternateAllelesToUse, + final boolean useBAQedPileup, + final GenomeLocParser locParser); protected int getFilteredDepth(ReadBackedPileup pileup) { diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java index 7ee7b0752c..1b73ef1d70 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java @@ -284,13 +284,14 @@ else if (p.isDeletion()) { private final static EnumSet allowableTypes = EnumSet.of(VariantContext.Type.INDEL, VariantContext.Type.MIXED); - public VariantContext getLikelihoods(RefMetaDataTracker tracker, - ReferenceContext ref, - Map contexts, - AlignmentContextUtils.ReadOrientation contextType, - GenotypePriors priors, - Allele alternateAlleleToUse, - boolean useBAQedPileup, GenomeLocParser locParser) { + public VariantContext getLikelihoods(final RefMetaDataTracker tracker, + final ReferenceContext ref, + final Map contexts, + final AlignmentContextUtils.ReadOrientation contextType, + final GenotypePriors priors, + final List alternateAllelesToUse, + final boolean useBAQedPileup, + final GenomeLocParser locParser) { if (tracker == null) return null; diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SNPGenotypeLikelihoodsCalculationModel.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SNPGenotypeLikelihoodsCalculationModel.java index 154612d255..dd21681f04 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SNPGenotypeLikelihoodsCalculationModel.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SNPGenotypeLikelihoodsCalculationModel.java @@ -64,7 +64,7 @@ public VariantContext getLikelihoods(final RefMetaDataTracker tracker, final Map contexts, final AlignmentContextUtils.ReadOrientation contextType, final GenotypePriors priors, - final Allele alternateAlleleToUse, + final List alternateAllelesToUse, final boolean useBAQedPileup, final GenomeLocParser locParser) { @@ -95,8 +95,8 @@ public VariantContext getLikelihoods(final RefMetaDataTracker tracker, } // find the alternate allele(s) that we should be using - if ( alternateAlleleToUse != null ) { - alleles.add(alternateAlleleToUse); + if ( alternateAllelesToUse != null ) { + alleles.addAll(alternateAllelesToUse); } else if ( useAlleleFromVCF ) { final VariantContext vc = UnifiedGenotyperEngine.getVCFromAllelesRod(tracker, ref, ref.getLocus(), true, logger, UAC.alleles); diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java index a60cc64f7b..05a977add7 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java @@ -237,14 +237,14 @@ public VariantCallContext calculateGenotypes(RefMetaDataTracker tracker, Referen // --------------------------------------------------------------------------------------------------------- // private method called by both UnifiedGenotyper and UGCalcLikelihoods entry points into the engine - private VariantContext calculateLikelihoods(RefMetaDataTracker tracker, ReferenceContext refContext, Map stratifiedContexts, AlignmentContextUtils.ReadOrientation type, Allele alternateAlleleToUse, boolean useBAQedPileup, final GenotypeLikelihoodsCalculationModel.Model model) { + private VariantContext calculateLikelihoods(RefMetaDataTracker tracker, ReferenceContext refContext, Map stratifiedContexts, AlignmentContextUtils.ReadOrientation type, List alternateAllelesToUse, boolean useBAQedPileup, final GenotypeLikelihoodsCalculationModel.Model model) { // initialize the data for this thread if that hasn't been done yet if ( glcm.get() == null ) { glcm.set(getGenotypeLikelihoodsCalculationObject(logger, UAC)); } - return glcm.get().get(model).getLikelihoods(tracker, refContext, stratifiedContexts, type, getGenotypePriors(model), alternateAlleleToUse, useBAQedPileup && BAQEnabledOnCMDLine, genomeLocParser); + return glcm.get().get(model).getLikelihoods(tracker, refContext, stratifiedContexts, type, getGenotypePriors(model), alternateAllelesToUse, useBAQedPileup && BAQEnabledOnCMDLine, genomeLocParser); } private VariantCallContext generateEmptyContext(RefMetaDataTracker tracker, ReferenceContext ref, Map stratifiedContexts, AlignmentContext rawContext) { @@ -398,16 +398,14 @@ else if ( UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_M //final boolean DEBUG_SLOD = false; // the overall lod - VariantContext vcOverall = calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.COMPLETE, vc.getAlternateAllele(0), false, model); - clearAFarray(AFresult.log10AlleleFrequencyLikelihoods); - clearAFarray(AFresult.log10AlleleFrequencyPosteriors); - afcm.get().getLog10PNonRef(vcOverall, getAlleleFrequencyPriors(model), AFresult); //double overallLog10PofNull = AFresult.log10AlleleFrequencyPosteriors[0]; double overallLog10PofF = MathUtils.log10sumLog10(AFresult.log10AlleleFrequencyPosteriors[0], 0); //if ( DEBUG_SLOD ) System.out.println("overallLog10PofF=" + overallLog10PofF); + List alternateAllelesToUse = builder.make().getAlternateAlleles(); + // the forward lod - VariantContext vcForward = calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.FORWARD, vc.getAlternateAllele(0), false, model); + VariantContext vcForward = calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.FORWARD, alternateAllelesToUse, false, model); clearAFarray(AFresult.log10AlleleFrequencyLikelihoods); clearAFarray(AFresult.log10AlleleFrequencyPosteriors); afcm.get().getLog10PNonRef(vcForward, getAlleleFrequencyPriors(model), AFresult); @@ -417,7 +415,7 @@ else if ( UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_M //if ( DEBUG_SLOD ) System.out.println("forwardLog10PofNull=" + forwardLog10PofNull + ", forwardLog10PofF=" + forwardLog10PofF); // the reverse lod - VariantContext vcReverse = calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.REVERSE, vc.getAlternateAllele(0), false, model); + VariantContext vcReverse = calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.REVERSE, alternateAllelesToUse, false, model); clearAFarray(AFresult.log10AlleleFrequencyLikelihoods); clearAFarray(AFresult.log10AlleleFrequencyPosteriors); afcm.get().getLog10PNonRef(vcReverse, getAlleleFrequencyPriors(model), AFresult);