diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ExactAFCalculationPerformanceTest.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ExactAFCalculationPerformanceTest.java index 7a8a2389a..628b4f880 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ExactAFCalculationPerformanceTest.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ExactAFCalculationPerformanceTest.java @@ -61,7 +61,7 @@ public class ExactAFCalculationPerformanceTest { final VariantContext vc = testBuilder.makeACTest(ACs, 0, nonTypePL); timer.start(); - final AFCalcResult result = calc.getLog10PNonRef(vc, priors); + final AFCalcResultTracker resultTracker = calc.getLog10PNonRef(vc, priors); final long runtime = timer.getElapsedTimeNano(); int otherAC = 0; @@ -72,7 +72,7 @@ public class ExactAFCalculationPerformanceTest { } final List columns = new LinkedList(coreValues); - columns.addAll(Arrays.asList(runtime, result.getnEvaluations(), nonTypePL, ACs[0], nAltSeg, otherAC)); + columns.addAll(Arrays.asList(runtime, resultTracker.getnEvaluations(), nonTypePL, ACs[0], nAltSeg, otherAC)); report.addRowList(columns); } } @@ -127,11 +127,11 @@ public class ExactAFCalculationPerformanceTest { vcb.genotypes(genotypes); timer.start(); - final AFCalcResult result = calc.getLog10PNonRef(vcb.make(), priors); + final AFCalcResultTracker resultTracker = calc.getLog10PNonRef(vcb.make(), priors); final long runtime = timer.getElapsedTimeNano(); final List columns = new LinkedList(coreValues); - columns.addAll(Arrays.asList(runtime, result.getnEvaluations(), nonTypePL, position)); + columns.addAll(Arrays.asList(runtime, resultTracker.getnEvaluations(), nonTypePL, position)); report.addRowList(columns); } } @@ -157,11 +157,11 @@ public class ExactAFCalculationPerformanceTest { final VariantContext vc = testBuilder.makeACTest(ac, nNonInformative, nonTypePL); timer.start(); - final AFCalcResult result = calc.getLog10PNonRef(vc, priors); + final AFCalcResultTracker resultTracker = calc.getLog10PNonRef(vc, priors); final long runtime = timer.getElapsedTimeNano(); final List columns = new LinkedList(coreValues); - columns.addAll(Arrays.asList(runtime, result.getnEvaluations(), nonTypePL, nNonInformative)); + columns.addAll(Arrays.asList(runtime, resultTracker.getnEvaluations(), nonTypePL, nNonInformative)); report.addRowList(columns); } } @@ -219,9 +219,9 @@ public class ExactAFCalculationPerformanceTest { final VariantContext vc = testBuilder.makeACTest(new int[]{ac}, 0, 100); final SimpleTimer timer = new SimpleTimer().start(); - final AFCalcResult result = testBuilder.makeModel().getLog10PNonRef(vc, testBuilder.makePriors()); + final AFCalcResultTracker resultTracker = testBuilder.makeModel().getLog10PNonRef(vc, testBuilder.makePriors()); final long runtime = timer.getElapsedTimeNano(); - logger.info("result " + result.getNormalizedPosteriorOfAFGTZero()); + logger.info("result " + resultTracker.getNormalizedPosteriorOfAFGTZero()); logger.info("runtime " + runtime); } diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/GeneralPloidyExactAFCalc.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/GeneralPloidyExactAFCalc.java index 77dff98c6..73c393c68 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/GeneralPloidyExactAFCalc.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/GeneralPloidyExactAFCalc.java @@ -78,8 +78,8 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { @Override public void computeLog10PNonRef(final VariantContext vc, final double[] log10AlleleFrequencyPriors, - final AFCalcResult result) { - combineSinglePools(vc.getGenotypes(), vc.getNAlleles(), ploidy, log10AlleleFrequencyPriors, result); + final AFCalcResultTracker resultTracker) { + combineSinglePools(vc.getGenotypes(), vc.getNAlleles(), ploidy, log10AlleleFrequencyPriors, resultTracker); } @@ -180,13 +180,13 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { * @param numAlleles Number of alternate alleles * @param ploidyPerPool Number of samples per pool * @param log10AlleleFrequencyPriors Frequency priors - * @param result object to fill with output values + * @param resultTracker object to fill with output values */ protected static void combineSinglePools(final GenotypesContext GLs, final int numAlleles, final int ploidyPerPool, final double[] log10AlleleFrequencyPriors, - final AFCalcResult result) { + final AFCalcResultTracker resultTracker) { final ArrayList genotypeLikelihoods = getGLs(GLs); @@ -203,9 +203,9 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { combinedPoolLikelihoods.add(set); for (int p=1; p stateTracker.getMaxLog10L()) @@ -263,7 +263,7 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { * @param log10AlleleFrequencyPriors Prior object * @param originalPloidy Total ploidy of original combined pool * @param newGLPloidy Ploidy of GL vector - * @param result AFResult object + * @param resultTracker AFResult object * @param stateTracker max likelihood observed so far * @param ACqueue Queue of conformations to compute * @param indexesToACset AC indices of objects in queue @@ -276,7 +276,7 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { final double[] log10AlleleFrequencyPriors, final int originalPloidy, final int newGLPloidy, - final AFCalcResult result, + final AFCalcResultTracker resultTracker, final StateTracker stateTracker, final LinkedList ACqueue, final HashMap indexesToACset) { @@ -284,7 +284,7 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { // compute likeihood in "set" of new set based on original likelihoods final int numAlleles = set.getACcounts().getCounts().length; final int newPloidy = set.getACsum(); - final double log10LofK = computeLofK(set, originalPool, newGL, log10AlleleFrequencyPriors, numAlleles, originalPloidy, newGLPloidy, result); + final double log10LofK = computeLofK(set, originalPool, newGL, log10AlleleFrequencyPriors, numAlleles, originalPloidy, newGLPloidy, resultTracker); // add to new pool @@ -339,11 +339,11 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { * @param ploidy2 Ploidy of second pool * @param numAlleles Number of alleles * @param log10AlleleFrequencyPriors Array of biallelic priors - * @param result Af calculation result object + * @param resultTracker Af calculation result object */ public static void combineMultiallelicPoolNaively(CombinedPoolLikelihoods originalPool, double[] yy, int ploidy1, int ploidy2, int numAlleles, final double[] log10AlleleFrequencyPriors, - final AFCalcResult result) { + final AFCalcResultTracker resultTracker) { /* final int dim1 = GenotypeLikelihoods.numLikelihoods(numAlleles, ploidy1); final int dim2 = GenotypeLikelihoods.numLikelihoods(numAlleles, ploidy2); @@ -397,7 +397,7 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { * @param numAlleles Number of alleles (including ref) * @param ploidy1 Ploidy of original pool (combined) * @param ploidy2 Ploidy of new pool - * @param result AFResult object + * @param resultTracker AFResult object * @return log-likehood of requested conformation */ private static double computeLofK(final ExactACset set, @@ -405,7 +405,7 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { final double[] secondGL, final double[] log10AlleleFrequencyPriors, final int numAlleles, final int ploidy1, final int ploidy2, - final AFCalcResult result) { + final AFCalcResultTracker resultTracker) { final int newPloidy = ploidy1 + ploidy2; @@ -423,8 +423,8 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { final double log10Lof0 = firstGLs.getGLOfACZero() + secondGL[HOM_REF_INDEX]; set.getLog10Likelihoods()[0] = log10Lof0; - result.setLog10LikelihoodOfAFzero(log10Lof0); - result.setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]); + resultTracker.setLog10LikelihoodOfAFzero(log10Lof0); + resultTracker.setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]); return log10Lof0; } else { @@ -467,14 +467,14 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { // update the MLE if necessary final int altCounts[] = Arrays.copyOfRange(set.getACcounts().getCounts(),1, set.getACcounts().getCounts().length); - result.updateMLEifNeeded(log10LofK, altCounts); + resultTracker.updateMLEifNeeded(log10LofK, altCounts); // apply the priors over each alternate allele for (final int ACcount : altCounts ) { if ( ACcount > 0 ) log10LofK += log10AlleleFrequencyPriors[ACcount]; } - result.updateMAPifNeeded(log10LofK, altCounts); + resultTracker.updateMAPifNeeded(log10LofK, altCounts); return log10LofK; } @@ -506,12 +506,12 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { * @param ploidy1 Ploidy of first pool (# of chromosomes in it) * @param ploidy2 Ploidy of second pool * @param log10AlleleFrequencyPriors Array of biallelic priors - * @param result Af calculation result object + * @param resultTracker Af calculation result object * @return Combined likelihood vector */ public static ProbabilityVector combineBiallelicPoolsNaively(final ProbabilityVector originalPool, final double[] newPLVector, final int ploidy1, final int ploidy2, final double[] log10AlleleFrequencyPriors, - final AFCalcResult result) { + final AFCalcResultTracker resultTracker) { final int newPloidy = ploidy1 + ploidy2; @@ -536,8 +536,8 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { final double log10Lof0 = x[0]+y[0]; - result.setLog10LikelihoodOfAFzero(log10Lof0); - result.setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]); + resultTracker.setLog10LikelihoodOfAFzero(log10Lof0); + resultTracker.setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]); double maxElement = log10Lof0; int maxElementIdx = 0; @@ -579,8 +579,8 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc { } alleleCounts[0] = k; - result.updateMLEifNeeded(combinedLikelihoods[k],alleleCounts); - result.updateMAPifNeeded(combinedLikelihoods[k] + log10AlleleFrequencyPriors[k],alleleCounts); + resultTracker.updateMLEifNeeded(combinedLikelihoods[k],alleleCounts); + resultTracker.updateMAPifNeeded(combinedLikelihoods[k] + log10AlleleFrequencyPriors[k],alleleCounts); } diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ExactAFCalculationModelUnitTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ExactAFCalculationModelUnitTest.java index ebab8d7e2..6402ca6c5 100644 --- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ExactAFCalculationModelUnitTest.java +++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ExactAFCalculationModelUnitTest.java @@ -76,11 +76,11 @@ public class ExactAFCalculationModelUnitTest extends BaseTest { } } - public AFCalcResult execute() { + public AFCalcResultTracker execute() { return getCalc().getLog10PNonRef(getVC(), getPriors()); } - public AFCalcResult executeRef() { + public AFCalcResultTracker executeRef() { final ExactAFCalc ref = new ReferenceDiploidExactAFCalc(getCalc().nSamples, getCalc().getMaxAltAlleles()); return ref.getLog10PNonRef(getVC(), getPriors()); } @@ -209,8 +209,8 @@ public class ExactAFCalculationModelUnitTest extends BaseTest { @Test(enabled = true, dataProvider = "GLsWithNonInformative", dependsOnMethods = "testGLs") public void testGLsWithNonInformative(GetGLsTest onlyInformative, GetGLsTest withNonInformative) { - final AFCalcResult expected = onlyInformative.execute(); - final AFCalcResult actual = withNonInformative.execute(); + final AFCalcResultTracker expected = onlyInformative.execute(); + final AFCalcResultTracker actual = withNonInformative.execute(); testResultSimple(withNonInformative); @@ -225,22 +225,22 @@ public class ExactAFCalculationModelUnitTest extends BaseTest { } private void testResultSimple(final GetGLsTest cfg) { - final AFCalcResult refResult = cfg.executeRef(); - final AFCalcResult result = cfg.execute(); + final AFCalcResultTracker refResultTracker = cfg.executeRef(); + final AFCalcResultTracker resultTracker = cfg.execute(); - compareToRefResult(refResult, result); + compareToRefResult(refResultTracker, resultTracker); - Assert.assertEquals(result.getNormalizedPosteriorOfAFzero() + result.getNormalizedPosteriorOfAFGTZero(), 1.0, 1e-4); + Assert.assertEquals(resultTracker.getNormalizedPosteriorOfAFzero() + resultTracker.getNormalizedPosteriorOfAFGTZero(), 1.0, 1e-4); // final int minNumberOfEvaluations = cfg.getVC().getCalledChrCount(); // Assert.assertTrue(result.getnEvaluations() >= minNumberOfEvaluations, // "Number of evaluations " + result.getnEvaluations() + " must be at least " + minNumberOfEvaluations); - Assert.assertNotNull(result.getAllelesUsedInGenotyping()); - Assert.assertTrue(cfg.getAlleles().containsAll(result.getAllelesUsedInGenotyping()), "Result object has alleles not in our initial allele list"); + Assert.assertNotNull(resultTracker.getAllelesUsedInGenotyping()); + Assert.assertTrue(cfg.getAlleles().containsAll(resultTracker.getAllelesUsedInGenotyping()), "Result object has alleles not in our initial allele list"); for ( int altAlleleI = 0; altAlleleI < cfg.numAltAlleles; altAlleleI++ ) { int expectedAlleleCount = cfg.getExpectedAltAC(altAlleleI); - int calcAC_MLE = result.getAlleleCountsOfMLE()[altAlleleI]; + int calcAC_MLE = resultTracker.getAlleleCountsOfMLE()[altAlleleI]; final Allele allele = cfg.getAlleles().get(altAlleleI+1); Assert.assertEquals(calcAC_MLE, expectedAlleleCount, "MLE AC not equal to expected AC for allele " + allele); @@ -257,20 +257,20 @@ public class ExactAFCalculationModelUnitTest extends BaseTest { // } } - private void compareToRefResult(final AFCalcResult refResult, - final AFCalcResult result) { + private void compareToRefResult(final AFCalcResultTracker refResultTracker, + final AFCalcResultTracker resultTracker) { final double TOLERANCE = 1; // MAP may not be equal // Assert.assertEquals(result.getAlleleCountsOfMAP(), refResult.getAlleleCountsOfMAP()); - Assert.assertEquals(result.getAlleleCountsOfMLE(), refResult.getAlleleCountsOfMLE()); - Assert.assertEquals(result.getAllelesUsedInGenotyping(), refResult.getAllelesUsedInGenotyping()); - Assert.assertEquals(result.getLog10LikelihoodOfAFzero(), refResult.getLog10LikelihoodOfAFzero(), TOLERANCE); + Assert.assertEquals(resultTracker.getAlleleCountsOfMLE(), refResultTracker.getAlleleCountsOfMLE()); + Assert.assertEquals(resultTracker.getAllelesUsedInGenotyping(), refResultTracker.getAllelesUsedInGenotyping()); + Assert.assertEquals(resultTracker.getLog10LikelihoodOfAFzero(), refResultTracker.getLog10LikelihoodOfAFzero(), TOLERANCE); // Assert.assertEquals(result.getLog10MAP(), refResult.getLog10MAP(), TOLERANCE); // Assert.assertEquals(result.getLog10MLE(), refResult.getLog10MLE(), TOLERANCE); // Assert.assertEquals(result.getLog10PosteriorOfAFzero(), refResult.getLog10PosteriorOfAFzero(), TOLERANCE); // Assert.assertEquals(result.getLog10PosteriorsMatrixSumWithoutAFzero(), refResult.getLog10PosteriorsMatrixSumWithoutAFzero(), TOLERANCE); - Assert.assertEquals(result.getNormalizedPosteriorOfAFGTZero(), refResult.getNormalizedPosteriorOfAFGTZero(), 0.5); - Assert.assertEquals(result.getNormalizedPosteriorOfAFzero(), refResult.getNormalizedPosteriorOfAFzero(), 0.5); + Assert.assertEquals(resultTracker.getNormalizedPosteriorOfAFGTZero(), refResultTracker.getNormalizedPosteriorOfAFGTZero(), 0.5); + Assert.assertEquals(resultTracker.getNormalizedPosteriorOfAFzero(), refResultTracker.getNormalizedPosteriorOfAFzero(), 0.5); } @Test(enabled = true, dataProvider = "Models") @@ -278,9 +278,9 @@ public class ExactAFCalculationModelUnitTest extends BaseTest { final Genotype BB = makePL(Arrays.asList(C, C), 20000000, 20000000, 0); GetGLsTest cfg = new GetGLsTest(calc, 1, Arrays.asList(BB, BB, BB), FLAT_3SAMPLE_PRIORS, "flat"); - final AFCalcResult result = cfg.execute(); + final AFCalcResultTracker resultTracker = cfg.execute(); - int calculatedAlleleCount = result.getAlleleCountsOfMAP()[0]; + int calculatedAlleleCount = resultTracker.getAlleleCountsOfMAP()[0]; Assert.assertEquals(calculatedAlleleCount, 6); } @@ -290,10 +290,10 @@ public class ExactAFCalculationModelUnitTest extends BaseTest { final Genotype AC = makePL(Arrays.asList(A, G), 100, 100, 100, 0, 100, 100); GetGLsTest cfg = new GetGLsTest(calc, 2, Arrays.asList(AB, AC), FLAT_3SAMPLE_PRIORS, "flat"); - final AFCalcResult result = cfg.execute(); + final AFCalcResultTracker resultTracker = cfg.execute(); - Assert.assertEquals(result.getAlleleCountsOfMAP()[0], 1); - Assert.assertEquals(result.getAlleleCountsOfMAP()[1], 1); + Assert.assertEquals(resultTracker.getAlleleCountsOfMAP()[0], 1); + Assert.assertEquals(resultTracker.getAlleleCountsOfMAP()[1], 1); } // -------------------------------------------------------------------------------- @@ -400,9 +400,9 @@ public class ExactAFCalculationModelUnitTest extends BaseTest { final VariantContextBuilder vcb = new VariantContextBuilder(vcRoot); vcb.genotypes(genotypes); - final AFCalcResult result = testBuilder.makeModel().getLog10PNonRef(vcb.make(), testBuilder.makePriors()); + final AFCalcResultTracker resultTracker = testBuilder.makeModel().getLog10PNonRef(vcb.make(), testBuilder.makePriors()); - Assert.assertEquals(result.getNormalizedPosteriorOfAFGTZero(), expectedPNonRef, tolerance, + Assert.assertEquals(resultTracker.getNormalizedPosteriorOfAFGTZero(), expectedPNonRef, tolerance, "Actual pNonRef not within tolerance " + tolerance + " of expected"); } @@ -432,17 +432,17 @@ public class ExactAFCalculationModelUnitTest extends BaseTest { final double refPrior = 1 - QualityUtils.qualToErrorProb(log10NonRefPrior); final double[] priors = MathUtils.toLog10(new double[]{refPrior, (1-refPrior) / 2, (1-refPrior) / 2}); GetGLsTest cfg = new GetGLsTest(model, 1, Arrays.asList(AB), priors, "pNonRef" + log10NonRefPrior); - final AFCalcResult result = cfg.execute(); - final int actualAC = result.getAlleleCountsOfMAP()[0]; + final AFCalcResultTracker resultTracker = cfg.execute(); + final int actualAC = resultTracker.getAlleleCountsOfMAP()[0]; final double pRefWithPrior = AB.getLikelihoods().getAsVector()[0] + priors[0]; final double pHetWithPrior = AB.getLikelihoods().getAsVector()[1] + priors[1]; final boolean expectNonRef = pRefWithPrior <= pHetWithPrior; if ( expectNonRef ) - Assert.assertTrue(result.getNormalizedPosteriorOfAFGTZero() > 0.5); + Assert.assertTrue(resultTracker.getNormalizedPosteriorOfAFGTZero() > 0.5); else - Assert.assertTrue(result.getNormalizedPosteriorOfAFGTZero() < 0.5); + Assert.assertTrue(resultTracker.getNormalizedPosteriorOfAFGTZero() < 0.5); final int expectedAC = expectNonRef ? 1 : 0; Assert.assertEquals(actualAC, expectedAC, @@ -468,8 +468,8 @@ public class ExactAFCalculationModelUnitTest extends BaseTest { final double nonRefPrior = (1-refPrior) / 2; final double[] priors = MathUtils.toLog10(new double[]{refPrior, nonRefPrior, nonRefPrior, nonRefPrior, nonRefPrior, nonRefPrior}); GetGLsTest cfg = new GetGLsTest(model, 2, Arrays.asList(AB, AC), priors, "pNonRef" + log10NonRefPrior); - final AFCalcResult result = cfg.execute(); - final int actualAC_AB = result.getAlleleCountsOfMAP()[0]; + final AFCalcResultTracker resultTracker = cfg.execute(); + final int actualAC_AB = resultTracker.getAlleleCountsOfMAP()[0]; final double pRefABWithPrior = AB.getLikelihoods().getAsVector()[0] + priors[0]; final double pHetABWithPrior = AB.getLikelihoods().getAsVector()[1] + priors[1]; @@ -480,7 +480,7 @@ public class ExactAFCalculationModelUnitTest extends BaseTest { final double nonRefPriorSecondAllele = Math.pow(nonRefPrior, 2); final double refPriorSecondAllele = 1 - nonRefPriorSecondAllele; - final int actualAC_AC = result.getAlleleCountsOfMAP()[1]; + final int actualAC_AC = resultTracker.getAlleleCountsOfMAP()[1]; final double pRefACWithPrior = AB.getLikelihoods().getAsVector()[0] + Math.log10(refPriorSecondAllele); final double pHetACWithPrior = AC.getLikelihoods().getAsVector()[3] + Math.log10(nonRefPriorSecondAllele); final int expectedAC_AC = pRefACWithPrior <= pHetACWithPrior ? 1 : 0; diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/GeneralPloidyAFCalculationModelUnitTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/GeneralPloidyAFCalculationModelUnitTest.java index 7381349ca..48f282901 100644 --- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/GeneralPloidyAFCalculationModelUnitTest.java +++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/GeneralPloidyAFCalculationModelUnitTest.java @@ -138,15 +138,15 @@ public class GeneralPloidyAFCalculationModelUnitTest extends BaseTest { @Test(dataProvider = "getGLs") public void testGLs(GetGLsTest cfg) { - final AFCalcResult result = new AFCalcResult(cfg.numAltAlleles); + final AFCalcResultTracker resultTracker = new AFCalcResultTracker(cfg.numAltAlleles); final int len = GeneralPloidyGenotypeLikelihoods.getNumLikelihoodElements(1 + cfg.numAltAlleles, cfg.ploidy * cfg.GLs.size()); double[] priors = new double[len]; // flat priors - GeneralPloidyExactAFCalc.combineSinglePools(cfg.GLs, 1 + cfg.numAltAlleles, cfg.ploidy, priors, result); + GeneralPloidyExactAFCalc.combineSinglePools(cfg.GLs, 1 + cfg.numAltAlleles, cfg.ploidy, priors, resultTracker); int nameIndex = 1; for ( int allele = 0; allele < cfg.numAltAlleles; allele++, nameIndex+=2 ) { int expectedAlleleCount = Integer.valueOf(cfg.name.substring(nameIndex, nameIndex+1)); - int calculatedAlleleCount = result.getAlleleCountsOfMAP()[allele]; + int calculatedAlleleCount = resultTracker.getAlleleCountsOfMAP()[allele]; // System.out.format( "%s Expected:%d Calc:%d\n",cfg.toString(),expectedAlleleCount, calculatedAlleleCount); Assert.assertEquals(calculatedAlleleCount, expectedAlleleCount); 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 cbe50b951..92e1c31f0 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 @@ -35,7 +35,7 @@ import org.broadinstitute.sting.gatk.contexts.ReferenceContext; import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; import org.broadinstitute.sting.gatk.walkers.annotator.VariantAnnotatorEngine; import org.broadinstitute.sting.gatk.walkers.genotyper.afcalc.AFCalc; -import org.broadinstitute.sting.gatk.walkers.genotyper.afcalc.AFCalcResult; +import org.broadinstitute.sting.gatk.walkers.genotyper.afcalc.AFCalcResultTracker; import org.broadinstitute.sting.utils.*; import org.broadinstitute.sting.utils.baq.BAQ; import org.broadinstitute.sting.utils.classloader.PluginManager; @@ -82,9 +82,6 @@ public class UnifiedGenotyperEngine { // the model used for calculating p(non-ref) private ThreadLocal afcm = new ThreadLocal(); - // the allele frequency likelihoods and posteriors (allocated once as an optimization) - private ThreadLocal alleleFrequencyCalculationResult = new ThreadLocal(); - // because the allele frequency priors are constant for a given i, we cache the results to avoid having to recompute everything private final double[] log10AlleleFrequencyPriorsSNPs; private final double[] log10AlleleFrequencyPriorsIndels; @@ -355,9 +352,7 @@ public class UnifiedGenotyperEngine { // initialize the data for this thread if that hasn't been done yet if ( afcm.get() == null ) { afcm.set(getAlleleFrequencyCalculationObject(N, logger, verboseWriter, UAC)); - alleleFrequencyCalculationResult.set(new AFCalcResult(UAC.MAX_ALTERNATE_ALLELES)); } - AFCalcResult AFresult = alleleFrequencyCalculationResult.get(); // estimate our confidence in a reference call and return if ( vc.getNSamples() == 0 ) { @@ -368,7 +363,7 @@ public class UnifiedGenotyperEngine { generateEmptyContext(tracker, refContext, stratifiedContexts, rawContext)); } - afcm.get().getLog10PNonRef(vc, getAlleleFrequencyPriors(model), AFresult); + AFCalcResultTracker AFresult = afcm.get().getLog10PNonRef(vc, getAlleleFrequencyPriors(model)); // is the most likely frequency conformation AC=0 for all alternate alleles? boolean bestGuessIsRef = true; @@ -474,7 +469,7 @@ public class UnifiedGenotyperEngine { // the forward lod VariantContext vcForward = calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.FORWARD, allAllelesToUse, false, model, perReadAlleleLikelihoodMap); - afcm.get().getLog10PNonRef(vcForward, getAlleleFrequencyPriors(model), AFresult); + AFresult = afcm.get().getLog10PNonRef(vcForward, getAlleleFrequencyPriors(model)); //double[] normalizedLog10Posteriors = MathUtils.normalizeFromLog10(AFresult.log10AlleleFrequencyPosteriors, true); double forwardLog10PofNull = AFresult.getLog10PosteriorOfAFzero(); double forwardLog10PofF = AFresult.getLog10PosteriorsMatrixSumWithoutAFzero(); @@ -482,7 +477,7 @@ public class UnifiedGenotyperEngine { // the reverse lod VariantContext vcReverse = calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.REVERSE, allAllelesToUse, false, model, perReadAlleleLikelihoodMap); - afcm.get().getLog10PNonRef(vcReverse, getAlleleFrequencyPriors(model), AFresult); + AFresult = afcm.get().getLog10PNonRef(vcReverse, getAlleleFrequencyPriors(model)); //normalizedLog10Posteriors = MathUtils.normalizeFromLog10(AFresult.log10AlleleFrequencyPosteriors, true); double reverseLog10PofNull = AFresult.getLog10PosteriorOfAFzero(); double reverseLog10PofF = AFresult.getLog10PosteriorsMatrixSumWithoutAFzero(); @@ -622,8 +617,6 @@ public class UnifiedGenotyperEngine { AFline.append(i + "/" + N + "\t"); AFline.append(String.format("%.2f\t", ((float)i)/N)); AFline.append(String.format("%.8f\t", getAlleleFrequencyPriors(model)[i])); - AFline.append(String.format("%.8f\t", alleleFrequencyCalculationResult.get().getLog10MLE())); - AFline.append(String.format("%.8f\t", alleleFrequencyCalculationResult.get().getLog10MAP())); verboseWriter.println(AFline.toString()); } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalc.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalc.java index 6ba73e59f..8245726b1 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalc.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalc.java @@ -73,6 +73,7 @@ public abstract class AFCalc implements Cloneable { private SimpleTimer callTimer = new SimpleTimer(); private PrintStream callReport = null; + private final AFCalcResultTracker resultTracker; protected AFCalc(final UnifiedArgumentCollection UAC, final int nSamples, final Logger logger, final PrintStream verboseWriter) { this(nSamples, UAC.MAX_ALTERNATE_ALLELES, UAC.MAX_ALTERNATE_ALLELES_FOR_INDELS, UAC.exactCallsLog, logger, verboseWriter); @@ -94,16 +95,7 @@ public abstract class AFCalc implements Cloneable { this.verboseWriter = verboseWriter; if ( exactCallsLog != null ) initializeOutputFile(exactCallsLog); - } - - /** - * @see #getLog10PNonRef(org.broadinstitute.sting.utils.variantcontext.VariantContext, double[], AFCalcResult) - * - * Allocates a new results object. Useful for testing but slow in practice. - */ - public final AFCalcResult getLog10PNonRef(final VariantContext vc, - final double[] log10AlleleFrequencyPriors) { - return getLog10PNonRef(vc, log10AlleleFrequencyPriors, new AFCalcResult(getMaxAltAlleles())); + this.resultTracker = new AFCalcResultTracker(Math.max(maxAltAlleles, maxAltAllelesForIndels)); } /** @@ -111,30 +103,27 @@ public abstract class AFCalc implements Cloneable { * * @param vc the VariantContext holding the alleles and sample information * @param log10AlleleFrequencyPriors a prior vector nSamples x 2 in length indicating the Pr(AF = i) - * @param result a pre-allocated (for efficiency) object to hold the result of the calculation * @return result (for programming convenience) */ - public final AFCalcResult getLog10PNonRef(final VariantContext vc, - final double[] log10AlleleFrequencyPriors, - final AFCalcResult result) { + public AFCalcResultTracker getLog10PNonRef(final VariantContext vc, final double[] log10AlleleFrequencyPriors) { if ( vc == null ) throw new IllegalArgumentException("VariantContext cannot be null"); if ( log10AlleleFrequencyPriors == null ) throw new IllegalArgumentException("priors vector cannot be null"); - if ( result == null ) throw new IllegalArgumentException("Results object cannot be null"); + if ( resultTracker == null ) throw new IllegalArgumentException("Results object cannot be null"); // reset the result, so we can store our new result there - result.reset(); + resultTracker.reset(); final VariantContext vcWorking = reduceScope(vc); callTimer.start(); - computeLog10PNonRef(vcWorking, log10AlleleFrequencyPriors, result); + computeLog10PNonRef(vcWorking, log10AlleleFrequencyPriors, resultTracker); final long nanoTime = callTimer.getElapsedTimeNano(); if ( callReport != null ) - printCallInfo(vcWorking, log10AlleleFrequencyPriors, nanoTime, result.getLog10PosteriorOfAFzero()); + printCallInfo(vcWorking, log10AlleleFrequencyPriors, nanoTime, resultTracker.getLog10PosteriorOfAFzero()); - result.setAllelesUsedInGenotyping(vcWorking.getAlleles()); - return result; + resultTracker.setAllelesUsedInGenotyping(vcWorking.getAlleles()); + return resultTracker; } // --------------------------------------------------------------------------- @@ -163,12 +152,12 @@ public abstract class AFCalc implements Cloneable { * * @param vc variant context with alleles and genotype likelihoods * @param log10AlleleFrequencyPriors priors - * @param result (pre-allocated) object to store results + * @param resultTracker (pre-allocated) object to store results */ // TODO -- add consistent requires among args public abstract void computeLog10PNonRef(final VariantContext vc, final double[] log10AlleleFrequencyPriors, - final AFCalcResult result); + final AFCalcResultTracker resultTracker); /** * Must be overridden by concrete subclasses diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResult.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResult.java index 5a8cab80b..e80dbc3d7 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResult.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResult.java @@ -26,38 +26,36 @@ package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc; import com.google.java.contract.Ensures; +import com.google.java.contract.Requires; import org.broadinstitute.sting.utils.MathUtils; +import org.broadinstitute.sting.utils.Utils; import org.broadinstitute.sting.utils.variantcontext.Allele; import java.util.Arrays; import java.util.List; /** - * Created by IntelliJ IDEA. - * User: ebanks - * Date: Dec 14, 2011 + * Describes the results of the AFCalc * - * Useful helper class to communicate the results of the allele frequency calculation + * Only the bare essentials are represented here, as all AFCalc models must return meaningful results for + * all of these fields. * - * TODO -- WHAT IS THE CONTRACT ON MAP AC AND P NON REF? + * Note that all of the values -- i.e. priors -- are checked now that they are meaningful, which means + * that users of this code can rely on the values coming out of these functions. */ public class AFCalcResult { - // These variables are intended to contain the MLE and MAP (and their corresponding allele counts) of the site over all alternate alleles - protected double log10MLE; - protected double log10MAP; + private final static int AF0 = 0; + private final static int AF1p = 1; + private final static int LOG_10_ARRAY_SIZES = 2; + + private final double[] log10LikelihoodsOfAC; + private final double[] log10PriorsOfAC; + private final double[] log10PosteriorsOfAC; + + /** + * The AC values for all ALT alleles at the MLE + */ private final int[] alleleCountsOfMLE; - private final int[] alleleCountsOfMAP; - - // The posteriors seen, not including that of AF=0 - private static final int POSTERIORS_CACHE_SIZE = 5000; - private final double[] log10PosteriorMatrixValues = new double[POSTERIORS_CACHE_SIZE]; - private int currentPosteriorsCacheIndex = 0; - protected Double log10PosteriorMatrixSum = null; - - // These variables are intended to contain the likelihood/posterior probability for the site's being monomorphic (i.e. AF=0 for all alternate alleles) - private double log10LikelihoodOfAFzero; - private double log10PosteriorOfAFzero; - private int[] AClimits; int nEvaluations = 0; @@ -68,36 +66,28 @@ public class AFCalcResult { /** * Create a results object capability of storing results for calls with up to maxAltAlleles - * - * @param maxAltAlleles an integer >= 1 */ - public AFCalcResult(final int maxAltAlleles) { - if ( maxAltAlleles < 1 ) throw new IllegalArgumentException("maxAltAlleles must be >= 0, saw " + maxAltAlleles); + public AFCalcResult(final int[] alleleCountsOfMLE, + final int nEvaluations, + final List allelesUsedInGenotyping, + final double[] log10LikelihoodsOfAC, + final double[] log10PriorsOfAC) { + if ( allelesUsedInGenotyping == null || allelesUsedInGenotyping.size() < 1 ) throw new IllegalArgumentException("allelesUsedInGenotyping must be non-null list of at least 1 value " + allelesUsedInGenotyping); + if ( alleleCountsOfMLE == null ) throw new IllegalArgumentException("alleleCountsOfMLE cannot be null"); + if ( alleleCountsOfMLE.length != allelesUsedInGenotyping.size() ) throw new IllegalArgumentException("alleleCountsOfMLE.length " + alleleCountsOfMLE.length + " != allelesUsedInGenotyping.size() " + allelesUsedInGenotyping.size()); + if ( nEvaluations < 0 ) throw new IllegalArgumentException("nEvaluations must be >= 0 but saw " + nEvaluations); + if ( log10LikelihoodsOfAC.length != 2 ) throw new IllegalArgumentException("log10LikelihoodsOfAC must have length equal 2"); + if ( log10PriorsOfAC.length != 2 ) throw new IllegalArgumentException("log10PriorsOfAC must have length equal 2"); + if ( ! goodLog10ProbVector(log10LikelihoodsOfAC, LOG_10_ARRAY_SIZES, false) ) throw new IllegalArgumentException("log10LikelihoodsOfAC are bad " + Utils.join(",", log10LikelihoodsOfAC)); + if ( ! goodLog10ProbVector(log10PriorsOfAC, LOG_10_ARRAY_SIZES, true) ) throw new IllegalArgumentException("log10priors are bad " + Utils.join(",", log10PriorsOfAC)); - alleleCountsOfMLE = new int[maxAltAlleles]; - alleleCountsOfMAP = new int[maxAltAlleles]; + this.alleleCountsOfMLE = alleleCountsOfMLE; + this.nEvaluations = nEvaluations; + this.allelesUsedInGenotyping = allelesUsedInGenotyping; - reset(); - } - - /** - * Get the log10 value of the probability mass at the MLE - * - * @return a log10 prob - */ - @Ensures("goodLog10Value(result)") - public double getLog10MLE() { - return log10MLE; - } - - /** - * Get the log10 value of the probability mass at the max. a posterior (MAP) - * - * @return a log10 prob - */ - @Ensures("goodLog10Value(result)") - public double getLog10MAP() { - return log10MAP; + this.log10LikelihoodsOfAC = Arrays.copyOf(log10LikelihoodsOfAC, LOG_10_ARRAY_SIZES); + this.log10PriorsOfAC = Arrays.copyOf(log10PriorsOfAC, LOG_10_ARRAY_SIZES); + this.log10PosteriorsOfAC = computePosteriors(log10LikelihoodsOfAC, log10PriorsOfAC); } /** @@ -115,18 +105,6 @@ public class AFCalcResult { return alleleCountsOfMLE; } - /** - * Returns a vector with maxAltAlleles values containing AC values at the MAP - * - * @see #getAlleleCountsOfMLE() for the encoding of results in this vector - * - * @return a non-null vector of ints - */ - @Ensures("result != null") - public int[] getAlleleCountsOfMAP() { - return alleleCountsOfMAP; - } - /** * Returns the number of cycles used to evaluate the pNonRef for this AF calculation * @@ -136,36 +114,6 @@ public class AFCalcResult { return nEvaluations; } - /** - * TODO -- eric what is this supposed to return? my unit tests don't do what I think they should - * - * @return - */ - public double getLog10PosteriorsMatrixSumWithoutAFzero() { - if ( log10PosteriorMatrixSum == null ) { - log10PosteriorMatrixSum = MathUtils.log10sumLog10(log10PosteriorMatrixValues, 0, currentPosteriorsCacheIndex); - } - return log10PosteriorMatrixSum; - } - - /** - * TODO -- eric what is this supposed to return? my unit tests don't do what I think they should - * - * @return - */ - public double getLog10LikelihoodOfAFzero() { - return log10LikelihoodOfAFzero; - } - - /** - * TODO -- eric what is this supposed to return? my unit tests don't do what I think they should - * - * @return - */ - public double getLog10PosteriorOfAFzero() { - return log10PosteriorOfAFzero; - } - /** * Get the list of alleles actually used in genotyping. * @@ -183,126 +131,107 @@ public class AFCalcResult { } /** - * Get the normalized -- across all AFs -- of AC == 0, NOT LOG10 - * @return - */ - // TODO -- this ensure cannot be enabled right now because the log10 inputs can be infinity, etc. - // TODO -- we should own these values in a more meaningful way and return good values in the case - // TODO -- where this happens, or instead thrown an error and have a function to say "was this calculation successful -// @Ensures({"result >= 0.0", "result <= 1.0"}) - public double getNormalizedPosteriorOfAFzero() { - return getNormalizedPosteriors()[0]; - } - - /** - * Get the normalized -- across all AFs -- of AC > 0, NOT LOG10 - * @return - */ - // TODO -- this ensure cannot be enabled right now because the log10 inputs can be infinity, etc. - // TODO -- we should own these values in a more meaningful way and return good values in the case - // TODO -- where this happens, or instead thrown an error and have a function to say "was this calculation successful - //@Ensures({"result >= 0.0", "result <= 1.0"}) - public double getNormalizedPosteriorOfAFGTZero() { - return getNormalizedPosteriors()[1]; - } - - private double[] getNormalizedPosteriors() { - final double[] posteriors = new double[]{ getLog10PosteriorOfAFzero(), getLog10PosteriorsMatrixSumWithoutAFzero() }; - return MathUtils.normalizeFromLog10(posteriors); - } - - public int[] getAClimits() { - return AClimits; - } - - // -------------------------------------------------------------------------------- - // - // Protected mutational methods only for use within the calculation models themselves - // - // -------------------------------------------------------------------------------- - - /** - * Reset the data in this results object, so that it can be used in a subsequent AF calculation + * Get the log10 normalized -- across all ACs -- posterior probability of AC == 0 * - * Resetting of the data is done by the calculation model itself, so shouldn't be done by callers any longer + * @return */ - protected void reset() { - log10MLE = log10MAP = log10LikelihoodOfAFzero = log10PosteriorOfAFzero = AFCalc.VALUE_NOT_CALCULATED; - for ( int i = 0; i < alleleCountsOfMLE.length; i++ ) { - alleleCountsOfMLE[i] = 0; - alleleCountsOfMAP[i] = 0; - } - currentPosteriorsCacheIndex = 0; - log10PosteriorMatrixSum = null; - allelesUsedInGenotyping = null; - nEvaluations = 0; + @Ensures({"goodLog10Value(result)"}) + public double getLog10PosteriorOfAFEq0() { + return log10PosteriorsOfAC[AF0]; } /** - * Tell this result we used one more evaluation cycle + * Get the log10 normalized -- across all ACs -- posterior probability of AC > 0 + * + * @return */ - protected void incNEvaluations() { - nEvaluations++; + @Ensures({"goodLog10Value(result)"}) + public double getLog10PosteriorOfAFGT0() { + return log10PosteriorsOfAC[AF1p]; } - protected void updateMLEifNeeded(final double log10LofK, final int[] alleleCountsForK) { - if ( log10LofK > log10MLE ) { - log10MLE = log10LofK; - for ( int i = 0; i < alleleCountsForK.length; i++ ) - alleleCountsOfMLE[i] = alleleCountsForK[i]; + /** + * Get the log10 unnormalized -- across all ACs -- likelihood of AC == 0 + * + * @return + */ + @Ensures({"goodLog10Value(result)"}) + public double getLog10LikelihoodOfAFEq0() { + return log10LikelihoodsOfAC[AF0]; + } + + /** + * Get the log10 unnormalized -- across all ACs -- likelihood of AC > 0 + * + * @return + */ + @Ensures({"goodLog10Value(result)"}) + public double getLog10LikelihoodOfAFGT0() { + return log10LikelihoodsOfAC[AF1p]; + } + + /** + * Get the log10 unnormalized -- across all ACs -- prior probability of AC == 0 + * + * @return + */ + @Ensures({"goodLog10Value(result)"}) + public double getLog10PriorOfAFEq0() { + return log10PriorsOfAC[AF0]; + } + + /** + * Get the log10 unnormalized -- across all ACs -- prior probability of AC > 0 + * + * @return + */ + @Ensures({"goodLog10Value(result)"}) + public double getLog10PriorOfAFGT0() { + return log10PriorsOfAC[AF1p]; + } + + /** + * Returns the log10 normalized posteriors given the log10 likelihoods and priors + * + * @param log10LikelihoodsOfAC + * @param log10PriorsOfAC + * + * @return freshly allocated log10 normalized posteriors vector + */ + @Requires("log10LikelihoodsOfAC.length == log10PriorsOfAC.length") + @Ensures("goodLog10ProbVector(result, LOG_10_ARRAY_SIZES, true)") + private static double[] computePosteriors(final double[] log10LikelihoodsOfAC, final double[] log10PriorsOfAC) { + final double[] log10UnnormalizedPosteriors = new double[log10LikelihoodsOfAC.length]; + for ( int i = 0; i < log10LikelihoodsOfAC.length; i++ ) + log10UnnormalizedPosteriors[i] = log10LikelihoodsOfAC[i] + log10PriorsOfAC[i]; + + return MathUtils.normalizeFromLog10(log10UnnormalizedPosteriors, true); + } + + /** + * Check that the log10 prob vector vector is well formed + * + * @param vector + * @param expectedSize + * @param shouldSumToOne + * + * @return true if vector is well-formed, false otherwise + */ + private static boolean goodLog10ProbVector(final double[] vector, final int expectedSize, final boolean shouldSumToOne) { + if ( vector.length != expectedSize ) return false; + + for ( final double pr : vector ) { + if ( pr > 0 ) return false; // log10 prob. vector should be < 0 + if ( Double.isInfinite(pr) || Double.isNaN(pr) ) return false; } - } - protected void updateMAPifNeeded(final double log10LofK, final int[] alleleCountsForK) { - addToPosteriorsCache(log10LofK); + if ( shouldSumToOne || MathUtils.compareDoubles(MathUtils.sumLog10(vector), 0.0, 1e-2) != 0 ) + return false; - if ( log10LofK > log10MAP ) { - log10MAP = log10LofK; - for ( int i = 0; i < alleleCountsForK.length; i++ ) - alleleCountsOfMAP[i] = alleleCountsForK[i]; - } - } - - private void addToPosteriorsCache(final double log10LofK) { - // add to the cache - log10PosteriorMatrixValues[currentPosteriorsCacheIndex++] = log10LofK; - - // if we've filled up the cache, then condense by summing up all of the values and placing the sum back into the first cell - if ( currentPosteriorsCacheIndex == POSTERIORS_CACHE_SIZE ) { - final double temporarySum = MathUtils.log10sumLog10(log10PosteriorMatrixValues, 0, currentPosteriorsCacheIndex); - log10PosteriorMatrixValues[0] = temporarySum; - currentPosteriorsCacheIndex = 1; - } - } - - protected void setLog10LikelihoodOfAFzero(final double log10LikelihoodOfAFzero) { - this.log10LikelihoodOfAFzero = log10LikelihoodOfAFzero; - if ( log10LikelihoodOfAFzero > log10MLE ) { - log10MLE = log10LikelihoodOfAFzero; - Arrays.fill(alleleCountsOfMLE, 0); - } - } - - protected void setLog10PosteriorOfAFzero(final double log10PosteriorOfAFzero) { - this.log10PosteriorOfAFzero = log10PosteriorOfAFzero; - if ( log10PosteriorOfAFzero > log10MAP ) { - log10MAP = log10PosteriorOfAFzero; - Arrays.fill(alleleCountsOfMAP, 0); - } - } - - protected void setAllelesUsedInGenotyping(List allelesUsedInGenotyping) { - if ( allelesUsedInGenotyping == null || allelesUsedInGenotyping.isEmpty() ) - throw new IllegalArgumentException("allelesUsedInGenotyping cannot be null or empty"); - - this.allelesUsedInGenotyping = allelesUsedInGenotyping; + return true; // everything is good } private static boolean goodLog10Value(final double result) { - return result <= 0.0 || Double.isInfinite(result) || Double.isNaN(result); - } - - protected void setAClimits(int[] AClimits) { - this.AClimits = AClimits; + return result <= 0.0 && ! Double.isInfinite(result) && ! Double.isNaN(result); } } \ No newline at end of file diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResultTracker.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResultTracker.java new file mode 100644 index 000000000..97e69be92 --- /dev/null +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResultTracker.java @@ -0,0 +1,308 @@ +/* + * Copyright (c) 2010. + * + * Permission is hereby granted, free of charge, to any person + * obtaining a copy of this software and associated documentation + * files (the "Software"), to deal in the Software without + * restriction, including without limitation the rights to use, + * copy, modify, merge, publish, distribute, sublicense, and/or sell + * copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following + * conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES + * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT + * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, + * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR + * THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc; + +import com.google.java.contract.Ensures; +import org.broadinstitute.sting.utils.MathUtils; +import org.broadinstitute.sting.utils.variantcontext.Allele; + +import java.util.Arrays; +import java.util.List; + +/** + * Created by IntelliJ IDEA. + * User: ebanks + * Date: Dec 14, 2011 + * + * Useful helper class to communicate the results of the allele frequency calculation + * + * TODO -- WHAT IS THE CONTRACT ON MAP AC AND P NON REF? + */ +public class AFCalcResultTracker { + // These variables are intended to contain the MLE and MAP (and their corresponding allele counts) of the site over all alternate alleles + protected double log10MLE; + protected double log10MAP; + private final int[] alleleCountsOfMLE; + private final int[] alleleCountsOfMAP; + + // The posteriors seen, not including that of AF=0 + private static final int POSTERIORS_CACHE_SIZE = 5000; + private final double[] log10PosteriorMatrixValues = new double[POSTERIORS_CACHE_SIZE]; + private int currentPosteriorsCacheIndex = 0; + protected Double log10PosteriorMatrixSum = null; + + // These variables are intended to contain the likelihood/posterior probability for the site's being monomorphic (i.e. AF=0 for all alternate alleles) + private double log10LikelihoodOfAFzero; + private double log10PosteriorOfAFzero; + private int[] AClimits; + + int nEvaluations = 0; + + /** + * The list of alleles actually used in computing the AF + */ + private List allelesUsedInGenotyping = null; + + /** + * Create a results object capability of storing results for calls with up to maxAltAlleles + * + * @param maxAltAlleles an integer >= 1 + */ + public AFCalcResultTracker(final int maxAltAlleles) { + if ( maxAltAlleles < 1 ) throw new IllegalArgumentException("maxAltAlleles must be >= 0, saw " + maxAltAlleles); + + alleleCountsOfMLE = new int[maxAltAlleles]; + alleleCountsOfMAP = new int[maxAltAlleles]; + + reset(); + } + + /** + * Get the log10 value of the probability mass at the MLE + * + * @return a log10 prob + */ + @Ensures("goodLog10Value(result)") + public double getLog10MLE() { + return log10MLE; + } + + /** + * Get the log10 value of the probability mass at the max. a posterior (MAP) + * + * @return a log10 prob + */ + @Ensures("goodLog10Value(result)") + public double getLog10MAP() { + return log10MAP; + } + + /** + * Returns a vector with maxAltAlleles values containing AC values at the MLE + * + * The values of the ACs for this call are stored in the getAllelesUsedInGenotyping order, + * starting from index 0 (i.e., the first alt allele is at 0). The vector is always + * maxAltAlleles in length, and so only the first getAllelesUsedInGenotyping.size() - 1 values + * are meaningful. + * + * @return a vector with allele counts, not all of which may be meaningful + */ + @Ensures("result != null") + public int[] getAlleleCountsOfMLE() { + return alleleCountsOfMLE; + } + + /** + * Returns a vector with maxAltAlleles values containing AC values at the MAP + * + * @see #getAlleleCountsOfMLE() for the encoding of results in this vector + * + * @return a non-null vector of ints + */ + @Ensures("result != null") + public int[] getAlleleCountsOfMAP() { + return alleleCountsOfMAP; + } + + /** + * Returns the number of cycles used to evaluate the pNonRef for this AF calculation + * + * @return the number of evaluations required to produce the answer for this AF calculation + */ + public int getnEvaluations() { + return nEvaluations; + } + + /** + * TODO -- eric what is this supposed to return? my unit tests don't do what I think they should + * + * @return + */ + public double getLog10PosteriorsMatrixSumWithoutAFzero() { + if ( log10PosteriorMatrixSum == null ) { + log10PosteriorMatrixSum = MathUtils.log10sumLog10(log10PosteriorMatrixValues, 0, currentPosteriorsCacheIndex); + } + return log10PosteriorMatrixSum; + } + + /** + * TODO -- eric what is this supposed to return? my unit tests don't do what I think they should + * + * @return + */ + public double getLog10LikelihoodOfAFzero() { + return log10LikelihoodOfAFzero; + } + + /** + * TODO -- eric what is this supposed to return? my unit tests don't do what I think they should + * + * @return + */ + public double getLog10PosteriorOfAFzero() { + return log10PosteriorOfAFzero; + } + + /** + * Get the list of alleles actually used in genotyping. + * + * Due to computational / implementation constraints this may be smaller than + * the actual list of alleles requested + * + * @return a non-empty list of alleles used during genotyping + */ + @Ensures({"result != null", "! result.isEmpty()"}) + public List getAllelesUsedInGenotyping() { + if ( allelesUsedInGenotyping == null ) + throw new IllegalStateException("allelesUsedInGenotyping requested but not yet set"); + + return allelesUsedInGenotyping; + } + + /** + * Get the normalized -- across all AFs -- of AC == 0, NOT LOG10 + * @return + */ + // TODO -- this ensure cannot be enabled right now because the log10 inputs can be infinity, etc. + // TODO -- we should own these values in a more meaningful way and return good values in the case + // TODO -- where this happens, or instead thrown an error and have a function to say "was this calculation successful +// @Ensures({"result >= 0.0", "result <= 1.0"}) + public double getNormalizedPosteriorOfAFzero() { + return getNormalizedPosteriors()[0]; + } + + /** + * Get the normalized -- across all AFs -- of AC > 0, NOT LOG10 + * @return + */ + // TODO -- this ensure cannot be enabled right now because the log10 inputs can be infinity, etc. + // TODO -- we should own these values in a more meaningful way and return good values in the case + // TODO -- where this happens, or instead thrown an error and have a function to say "was this calculation successful + //@Ensures({"result >= 0.0", "result <= 1.0"}) + public double getNormalizedPosteriorOfAFGTZero() { + return getNormalizedPosteriors()[1]; + } + + private double[] getNormalizedPosteriors() { + final double[] posteriors = new double[]{ getLog10PosteriorOfAFzero(), getLog10PosteriorsMatrixSumWithoutAFzero() }; + return MathUtils.normalizeFromLog10(posteriors); + } + + public int[] getAClimits() { + return AClimits; + } + + // -------------------------------------------------------------------------------- + // + // Protected mutational methods only for use within the calculation models themselves + // + // -------------------------------------------------------------------------------- + + /** + * Reset the data in this results object, so that it can be used in a subsequent AF calculation + * + * Resetting of the data is done by the calculation model itself, so shouldn't be done by callers any longer + */ + protected void reset() { + log10MLE = log10MAP = log10LikelihoodOfAFzero = log10PosteriorOfAFzero = AFCalc.VALUE_NOT_CALCULATED; + for ( int i = 0; i < alleleCountsOfMLE.length; i++ ) { + alleleCountsOfMLE[i] = 0; + alleleCountsOfMAP[i] = 0; + } + currentPosteriorsCacheIndex = 0; + log10PosteriorMatrixSum = null; + allelesUsedInGenotyping = null; + nEvaluations = 0; + } + + /** + * Tell this result we used one more evaluation cycle + */ + protected void incNEvaluations() { + nEvaluations++; + } + + protected void updateMLEifNeeded(final double log10LofK, final int[] alleleCountsForK) { + if ( log10LofK > log10MLE ) { + log10MLE = log10LofK; + for ( int i = 0; i < alleleCountsForK.length; i++ ) + alleleCountsOfMLE[i] = alleleCountsForK[i]; + } + } + + protected void updateMAPifNeeded(final double log10LofK, final int[] alleleCountsForK) { + addToPosteriorsCache(log10LofK); + + if ( log10LofK > log10MAP ) { + log10MAP = log10LofK; + for ( int i = 0; i < alleleCountsForK.length; i++ ) + alleleCountsOfMAP[i] = alleleCountsForK[i]; + } + } + + private void addToPosteriorsCache(final double log10LofK) { + // add to the cache + log10PosteriorMatrixValues[currentPosteriorsCacheIndex++] = log10LofK; + + // if we've filled up the cache, then condense by summing up all of the values and placing the sum back into the first cell + if ( currentPosteriorsCacheIndex == POSTERIORS_CACHE_SIZE ) { + final double temporarySum = MathUtils.log10sumLog10(log10PosteriorMatrixValues, 0, currentPosteriorsCacheIndex); + log10PosteriorMatrixValues[0] = temporarySum; + currentPosteriorsCacheIndex = 1; + } + } + + protected void setLog10LikelihoodOfAFzero(final double log10LikelihoodOfAFzero) { + this.log10LikelihoodOfAFzero = log10LikelihoodOfAFzero; + if ( log10LikelihoodOfAFzero > log10MLE ) { + log10MLE = log10LikelihoodOfAFzero; + Arrays.fill(alleleCountsOfMLE, 0); + } + } + + protected void setLog10PosteriorOfAFzero(final double log10PosteriorOfAFzero) { + this.log10PosteriorOfAFzero = log10PosteriorOfAFzero; + if ( log10PosteriorOfAFzero > log10MAP ) { + log10MAP = log10PosteriorOfAFzero; + Arrays.fill(alleleCountsOfMAP, 0); + } + } + + protected void setAllelesUsedInGenotyping(List allelesUsedInGenotyping) { + if ( allelesUsedInGenotyping == null || allelesUsedInGenotyping.isEmpty() ) + throw new IllegalArgumentException("allelesUsedInGenotyping cannot be null or empty"); + + this.allelesUsedInGenotyping = allelesUsedInGenotyping; + } + + private static boolean goodLog10Value(final double result) { + return result <= 0.0 || Double.isInfinite(result) || Double.isNaN(result); + } + + protected void setAClimits(int[] AClimits) { + this.AClimits = AClimits; + } +} \ No newline at end of file diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ConstrainedDiploidExactAFCalc.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ConstrainedDiploidExactAFCalc.java index 3257be97b..1b021aa77 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ConstrainedDiploidExactAFCalc.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ConstrainedDiploidExactAFCalc.java @@ -19,9 +19,9 @@ public class ConstrainedDiploidExactAFCalc extends DiploidExactAFCalc { super(UAC, N, logger, verboseWriter); } - protected StateTracker makeMaxLikelihood(final VariantContext vc, final AFCalcResult result) { + protected StateTracker makeMaxLikelihood(final VariantContext vc, final AFCalcResultTracker resultTracker) { final int[] maxACsToConsider = computeMaxACs(vc); - result.setAClimits(maxACsToConsider); + resultTracker.setAClimits(maxACsToConsider); return new StateTracker(maxACsToConsider); } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/DiploidExactAFCalc.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/DiploidExactAFCalc.java index 48e4e8359..0dac2653d 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/DiploidExactAFCalc.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/DiploidExactAFCalc.java @@ -42,12 +42,12 @@ public abstract class DiploidExactAFCalc extends ExactAFCalc { super(UAC, N, logger, verboseWriter); } - protected abstract StateTracker makeMaxLikelihood(final VariantContext vc, final AFCalcResult result); + protected abstract StateTracker makeMaxLikelihood(final VariantContext vc, final AFCalcResultTracker resultTracker); @Override public void computeLog10PNonRef(final VariantContext vc, final double[] log10AlleleFrequencyPriors, - final AFCalcResult result) { + final AFCalcResultTracker resultTracker) { final int numAlternateAlleles = vc.getNAlleles() - 1; final ArrayList genotypeLikelihoods = getGLs(vc.getGenotypes()); final int numSamples = genotypeLikelihoods.size()-1; @@ -66,16 +66,16 @@ public abstract class DiploidExactAFCalc extends ExactAFCalc { indexesToACset.put(zeroSet.getACcounts(), zeroSet); // keep processing while we have AC conformations that need to be calculated - final StateTracker stateTracker = makeMaxLikelihood(vc, result); + final StateTracker stateTracker = makeMaxLikelihood(vc, resultTracker); while ( !ACqueue.isEmpty() ) { - result.incNEvaluations(); // keep track of the number of evaluations + resultTracker.incNEvaluations(); // keep track of the number of evaluations // compute log10Likelihoods final ExactACset set = ACqueue.remove(); if ( stateTracker.withinMaxACs(set.getACcounts()) ) { - final double log10LofKs = calculateAlleleCountConformation(set, genotypeLikelihoods, stateTracker, numChr, ACqueue, indexesToACset, log10AlleleFrequencyPriors, result); + final double log10LofKs = calculateAlleleCountConformation(set, genotypeLikelihoods, stateTracker, numChr, ACqueue, indexesToACset, log10AlleleFrequencyPriors, resultTracker); // adjust max likelihood seen if needed stateTracker.update(log10LofKs, set.getACcounts()); @@ -161,13 +161,13 @@ public abstract class DiploidExactAFCalc extends ExactAFCalc { final LinkedList ACqueue, final HashMap indexesToACset, final double[] log10AlleleFrequencyPriors, - final AFCalcResult result) { + final AFCalcResultTracker resultTracker) { //if ( DEBUG ) // System.out.printf(" *** computing LofK for set=%s%n", set.ACcounts); // compute the log10Likelihoods - computeLofK(set, genotypeLikelihoods, log10AlleleFrequencyPriors, result); + computeLofK(set, genotypeLikelihoods, log10AlleleFrequencyPriors, resultTracker); final double log10LofK = set.getLog10Likelihoods()[set.getLog10Likelihoods().length-1]; @@ -250,7 +250,7 @@ public abstract class DiploidExactAFCalc extends ExactAFCalc { private void computeLofK(final ExactACset set, final ArrayList genotypeLikelihoods, final double[] log10AlleleFrequencyPriors, - final AFCalcResult result) { + final AFCalcResultTracker resultTracker) { set.getLog10Likelihoods()[0] = 0.0; // the zero case final int totalK = set.getACsum(); @@ -261,8 +261,8 @@ public abstract class DiploidExactAFCalc extends ExactAFCalc { set.getLog10Likelihoods()[j] = set.getLog10Likelihoods()[j-1] + genotypeLikelihoods.get(j)[HOM_REF_INDEX]; final double log10Lof0 = set.getLog10Likelihoods()[set.getLog10Likelihoods().length-1]; - result.setLog10LikelihoodOfAFzero(log10Lof0); - result.setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]); + resultTracker.setLog10LikelihoodOfAFzero(log10Lof0); + resultTracker.setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]); return; } @@ -284,14 +284,14 @@ public abstract class DiploidExactAFCalc extends ExactAFCalc { double log10LofK = set.getLog10Likelihoods()[set.getLog10Likelihoods().length-1]; // update the MLE if necessary - result.updateMLEifNeeded(log10LofK, set.getACcounts().getCounts()); + resultTracker.updateMLEifNeeded(log10LofK, set.getACcounts().getCounts()); // apply the priors over each alternate allele for ( final int ACcount : set.getACcounts().getCounts() ) { if ( ACcount > 0 ) log10LofK += log10AlleleFrequencyPriors[ACcount]; } - result.updateMAPifNeeded(log10LofK, set.getACcounts().getCounts()); + resultTracker.updateMAPifNeeded(log10LofK, set.getACcounts().getCounts()); } private void pushData(final ExactACset targetSet, diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/IndependentAllelesDiploidExactAFCalc.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/IndependentAllelesDiploidExactAFCalc.java index 56ef1ed3b..b74923086 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/IndependentAllelesDiploidExactAFCalc.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/IndependentAllelesDiploidExactAFCalc.java @@ -52,31 +52,31 @@ public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc { } @Override - protected StateTracker makeMaxLikelihood(VariantContext vc, AFCalcResult result) { - return refModel.makeMaxLikelihood(vc, result); + protected StateTracker makeMaxLikelihood(VariantContext vc, AFCalcResultTracker resultTracker) { + return refModel.makeMaxLikelihood(vc, resultTracker); } @Override public void computeLog10PNonRef(final VariantContext vc, final double[] log10AlleleFrequencyPriors, - final AFCalcResult result) { - final List independentResults = computeLog10PNonRefForEachAllele(vc, log10AlleleFrequencyPriors); - combineIndependentPNonRefs(vc, independentResults, log10AlleleFrequencyPriors, result); + final AFCalcResultTracker resultTracker) { + final List independentResultTrackers = computeLog10PNonRefForEachAllele(vc, log10AlleleFrequencyPriors); + combineIndependentPNonRefs(vc, independentResultTrackers, log10AlleleFrequencyPriors, resultTracker); } - protected List computeLog10PNonRefForEachAllele(final VariantContext vc, + protected List computeLog10PNonRefForEachAllele(final VariantContext vc, final double[] log10AlleleFrequencyPriors) { final int nAltAlleles = vc.getNAlleles() - 1; - final List results = new ArrayList(nAltAlleles); + final List resultTrackers = new ArrayList(nAltAlleles); for ( int altI = 0; altI < nAltAlleles; altI++ ) { final List biallelic = Arrays.asList(vc.getReference(), vc.getAlternateAllele(altI)); final VariantContext subvc = biallelicCombinedGLs(vc, biallelic, altI + 1); - final AFCalcResult result = refModel.getLog10PNonRef(subvc, log10AlleleFrequencyPriors); - results.add(result); + final AFCalcResultTracker resultTracker = refModel.getLog10PNonRef(subvc, log10AlleleFrequencyPriors); + resultTrackers.add(resultTracker); } - return results; + return resultTrackers; } protected VariantContext biallelicCombinedGLs(final VariantContext rootVC, final List biallelic, final int allele2) { @@ -138,36 +138,36 @@ public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc { * Takes each independent result and merges it into the final result object * * @param independentPNonRefs the pNonRef result for each allele independently - * @param result the destination for the combined result + * @param resultTracker the destination for the combined result */ protected void combineIndependentPNonRefs(final VariantContext vc, - final List independentPNonRefs, + final List independentPNonRefs, final double[] log10AlleleFrequencyPriors, - final AFCalcResult result) { + final AFCalcResultTracker resultTracker) { final int nChrom = vc.getNSamples() * 2; - result.reset(); + resultTracker.reset(); // both the likelihood and the posterior of AF=0 are the same for all alleles // TODO -- check and ensure this is true - result.setLog10LikelihoodOfAFzero(independentPNonRefs.get(0).getLog10LikelihoodOfAFzero()); - result.setLog10PosteriorOfAFzero(independentPNonRefs.get(0).getLog10PosteriorOfAFzero()); - result.log10PosteriorMatrixSum = 0.0; + resultTracker.setLog10LikelihoodOfAFzero(independentPNonRefs.get(0).getLog10LikelihoodOfAFzero()); + resultTracker.setLog10PosteriorOfAFzero(independentPNonRefs.get(0).getLog10PosteriorOfAFzero()); + resultTracker.log10PosteriorMatrixSum = 0.0; int altI = 0; - for ( final AFCalcResult independentPNonRef : independentPNonRefs ) { - result.log10MLE += independentPNonRef.getLog10MLE(); + for ( final AFCalcResultTracker independentPNonRef : independentPNonRefs ) { + resultTracker.log10MLE += independentPNonRef.getLog10MLE(); // TODO -- technically double counting some posterior mass - result.log10MAP += independentPNonRef.getLog10MAP(); + resultTracker.log10MAP += independentPNonRef.getLog10MAP(); // TODO -- technically double counting some posterior mass - result.log10PosteriorMatrixSum += independentPNonRef.getLog10PosteriorsMatrixSumWithoutAFzero(); + resultTracker.log10PosteriorMatrixSum += independentPNonRef.getLog10PosteriorsMatrixSumWithoutAFzero(); - result.getAlleleCountsOfMAP()[altI] = independentPNonRef.getAlleleCountsOfMAP()[0]; - result.getAlleleCountsOfMLE()[altI] = independentPNonRef.getAlleleCountsOfMLE()[0]; + resultTracker.getAlleleCountsOfMAP()[altI] = independentPNonRef.getAlleleCountsOfMAP()[0]; + resultTracker.getAlleleCountsOfMLE()[altI] = independentPNonRef.getAlleleCountsOfMLE()[0]; - result.nEvaluations += independentPNonRef.nEvaluations; + resultTracker.nEvaluations += independentPNonRef.nEvaluations; altI++; } } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ReferenceDiploidExactAFCalc.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ReferenceDiploidExactAFCalc.java index 7ae710e73..9aa93061f 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ReferenceDiploidExactAFCalc.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ReferenceDiploidExactAFCalc.java @@ -15,7 +15,7 @@ public class ReferenceDiploidExactAFCalc extends DiploidExactAFCalc { super(UAC, N, logger, verboseWriter); } - protected StateTracker makeMaxLikelihood(final VariantContext vc, final AFCalcResult result) { + protected StateTracker makeMaxLikelihood(final VariantContext vc, final AFCalcResultTracker resultTracker) { return new StateTracker(); } } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/GLBasedSampleSelector.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/GLBasedSampleSelector.java index 11b4ca3cc..006c303dc 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/GLBasedSampleSelector.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/GLBasedSampleSelector.java @@ -23,7 +23,7 @@ */ package org.broadinstitute.sting.gatk.walkers.validation.validationsiteselector; -import org.broadinstitute.sting.gatk.walkers.genotyper.afcalc.AFCalcResult; +import org.broadinstitute.sting.gatk.walkers.genotyper.afcalc.AFCalcResultTracker; import org.broadinstitute.sting.gatk.walkers.genotyper.afcalc.DiploidExactAFCalc; import org.broadinstitute.sting.gatk.walkers.genotyper.afcalc.ReferenceDiploidExactAFCalc; import org.broadinstitute.sting.utils.variantcontext.VariantContext; @@ -54,10 +54,10 @@ public class GLBasedSampleSelector extends SampleSelector { flatPriors = new double[1+2*samples.size()]; AFCalculator = new ReferenceDiploidExactAFCalc(samples.size(), 4); } - AFCalcResult result = new AFCalcResult(vc.getAlternateAlleles().size()); - AFCalculator.computeLog10PNonRef(subContext, flatPriors, result); + AFCalcResultTracker resultTracker = new AFCalcResultTracker(vc.getAlternateAlleles().size()); + AFCalculator.computeLog10PNonRef(subContext, flatPriors, resultTracker); // do we want to let this qual go up or down? - if ( result.getLog10PosteriorOfAFzero() < referenceLikelihood ) { + if ( resultTracker.getLog10PosteriorOfAFzero() < referenceLikelihood ) { return true; }