This is less correct than my previous change but it's what UGv1 does and now is not the right time to start mucking with things.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@4506 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
ebanks 2010-10-15 18:56:45 +00:00
parent 471c18054f
commit f78ff08e2b
4 changed files with 17 additions and 11 deletions

View File

@ -81,12 +81,14 @@ public abstract class AlleleFrequencyCalculationModel implements Cloneable {
* @param GLs genotype likelihoods
* @param log10AlleleFrequencyPriors priors
* @param log10AlleleFrequencyPosteriors array (pre-allocated) to store results
* @param minFrequencyToCalculate the minimum frequency which needs to be calculated
*/
public abstract void getLog10PNonRef(RefMetaDataTracker tracker,
ReferenceContext ref,
Map<String, BiallelicGenotypeLikelihoods> GLs,
double[] log10AlleleFrequencyPriors,
double[] log10AlleleFrequencyPosteriors);
double[] log10AlleleFrequencyPosteriors,
int minFrequencyToCalculate);
/**
* Can be overridden by concrete subclasses

View File

@ -29,7 +29,6 @@ import org.apache.log4j.Logger;
import org.broad.tribble.util.variantcontext.Allele;
import org.broad.tribble.util.variantcontext.Genotype;
import org.broad.tribble.util.variantcontext.GenotypeLikelihoods;
import org.broad.tribble.util.variantcontext.VariantContext;
import org.broad.tribble.vcf.VCFConstants;
import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
import org.broadinstitute.sting.utils.*;
@ -54,7 +53,8 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel {
ReferenceContext ref,
Map<String, BiallelicGenotypeLikelihoods> GLs,
double[] log10AlleleFrequencyPriors,
double[] log10AlleleFrequencyPosteriors) {
double[] log10AlleleFrequencyPosteriors,
int minFrequencyToCalculate) {
// Math requires linear math to make efficient updates.
double[] alleleFrequencyPriors = MathUtils.normalizeFromLog10(log10AlleleFrequencyPriors);
@ -233,7 +233,7 @@ public class ExactAFCalculationModel extends AlleleFrequencyCalculationModel {
HashMap<String, Object> attributes = new HashMap<String, Object>();
ArrayList<Allele> myAlleles = new ArrayList<Allele>();
attributes.put(VCFConstants.DEPTH_KEY, contexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).size());
double qual = 0.0;
double qual;
double[] posteriors = GLs.get(sample).getPosteriors();
if (bestGTguess == 0) {

View File

@ -48,7 +48,8 @@ public class GridSearchAFEstimation extends AlleleFrequencyCalculationModel {
ReferenceContext ref,
Map<String, BiallelicGenotypeLikelihoods> GLs,
double[] log10AlleleFrequencyPriors,
double[] log10AlleleFrequencyPosteriors) {
double[] log10AlleleFrequencyPosteriors,
int minFrequencyToCalculate) {
initializeAFMatrix(GLs);
@ -68,12 +69,17 @@ public class GridSearchAFEstimation extends AlleleFrequencyCalculationModel {
// an optimization to speed up the calculation: if we are beyond the local maximum such
// that subsequent likelihoods won't factor into the confidence score, just quit
if ( maxLikelihoodSeen - log10AlleleFrequencyPosteriors[i] > LOG10_OPTIMIZATION_EPSILON )
if ( i >= minFrequencyToCalculate && maxLikelihoodSeen - log10AlleleFrequencyPosteriors[i] > LOG10_OPTIMIZATION_EPSILON )
return;
if ( log10AlleleFrequencyPosteriors[i] > maxLikelihoodSeen )
maxLikelihoodSeen = log10AlleleFrequencyPosteriors[i];
}
// TODO: we really need to get rid of this and the minFrequencyToCalculate argument
// it's possible that we need to calculate higher frequencies
for (int i = maxAlleleFrequencyToTest; i <= minFrequencyToCalculate; i++)
log10AlleleFrequencyPosteriors[i] = log10AlleleFrequencyPosteriors[maxAlleleFrequencyToTest];
}
/**

View File

@ -157,7 +157,7 @@ public class UnifiedGenotyperEngine {
// 'zero' out the AFs (so that we don't have to worry if not all samples have reads at this position)
clearAFarray(log10AlleleFrequencyPosteriors.get());
afcm.get().getLog10PNonRef(tracker, refContext, GLs, log10AlleleFrequencyPriors, log10AlleleFrequencyPosteriors.get());
afcm.get().getLog10PNonRef(tracker, refContext, GLs, log10AlleleFrequencyPriors, log10AlleleFrequencyPosteriors.get(), 0);
// find the most likely frequency
int bestAFguess = MathUtils.maxElementIndex(log10AlleleFrequencyPosteriors.get());
@ -234,9 +234,8 @@ public class UnifiedGenotyperEngine {
GLs.clear();
glcm.get().getLikelihoods(tracker, refContext, stratifiedContexts, StratifiedAlignmentContext.StratifiedContextType.FORWARD, genotypePriors, GLs);
clearAFarray(log10AlleleFrequencyPosteriors.get());
afcm.get().getLog10PNonRef(tracker, refContext, GLs, log10AlleleFrequencyPriors, log10AlleleFrequencyPosteriors.get());
afcm.get().getLog10PNonRef(tracker, refContext, GLs, log10AlleleFrequencyPriors, log10AlleleFrequencyPosteriors.get(), bestAFguess);
double forwardLog10PofNull = log10AlleleFrequencyPosteriors.get()[0];
bestAFguess = MathUtils.maxElementIndex(log10AlleleFrequencyPosteriors.get());
double forwardLog10PofF = log10AlleleFrequencyPosteriors.get()[bestAFguess];
//System.out.println("forwardLog10PofNull=" + forwardLog10PofNull + ", forwardLog10PofF=" + forwardLog10PofF);
@ -244,9 +243,8 @@ public class UnifiedGenotyperEngine {
GLs.clear();
glcm.get().getLikelihoods(tracker, refContext, stratifiedContexts, StratifiedAlignmentContext.StratifiedContextType.REVERSE, genotypePriors, GLs);
clearAFarray(log10AlleleFrequencyPosteriors.get());
afcm.get().getLog10PNonRef(tracker, refContext, GLs, log10AlleleFrequencyPriors, log10AlleleFrequencyPosteriors.get());
afcm.get().getLog10PNonRef(tracker, refContext, GLs, log10AlleleFrequencyPriors, log10AlleleFrequencyPosteriors.get(), bestAFguess);
double reverseLog10PofNull = log10AlleleFrequencyPosteriors.get()[0];
bestAFguess = MathUtils.maxElementIndex(log10AlleleFrequencyPosteriors.get());
double reverseLog10PofF = log10AlleleFrequencyPosteriors.get()[bestAFguess];
//System.out.println("reverseLog10PofNull=" + reverseLog10PofNull + ", reverseLog10PofF=" + reverseLog10PofF);