Intermediate commit on simplifying AFCalcResult

-- Renamed old class AFCalcResultTracker.  This object is now allocated by the AFCalc itself, since it is heavy-weight and was badly optimized in the UG with a thread-local variable. Now, since there's already a AFCalc thread-local there, we get that optimization for free.
-- Removed the interface to provide the AFCalcResultTracker to getlog10PNonRef.
-- Wrote new, clean but unused AFCalcResult object that will soon replace the tracker as the external interface to the AFCalc model results, leaving the tracker as an internal tracker structure.  This will allow me to (1) finally test things exhaustively, as the contracts on this class are clear (2) finalize the IndependentAllelesDiploidExactAFCalc class as it can work with a meaningfully defined result across each object
This commit is contained in:
Mark DePristo 2012-10-08 11:04:57 -04:00
parent c82aa01e0e
commit 06687bfaf6
13 changed files with 560 additions and 341 deletions

View File

@ -61,7 +61,7 @@ public class ExactAFCalculationPerformanceTest {
final VariantContext vc = testBuilder.makeACTest(ACs, 0, nonTypePL); final VariantContext vc = testBuilder.makeACTest(ACs, 0, nonTypePL);
timer.start(); timer.start();
final AFCalcResult result = calc.getLog10PNonRef(vc, priors); final AFCalcResultTracker resultTracker = calc.getLog10PNonRef(vc, priors);
final long runtime = timer.getElapsedTimeNano(); final long runtime = timer.getElapsedTimeNano();
int otherAC = 0; int otherAC = 0;
@ -72,7 +72,7 @@ public class ExactAFCalculationPerformanceTest {
} }
final List<Object> columns = new LinkedList<Object>(coreValues); final List<Object> columns = new LinkedList<Object>(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); report.addRowList(columns);
} }
} }
@ -127,11 +127,11 @@ public class ExactAFCalculationPerformanceTest {
vcb.genotypes(genotypes); vcb.genotypes(genotypes);
timer.start(); timer.start();
final AFCalcResult result = calc.getLog10PNonRef(vcb.make(), priors); final AFCalcResultTracker resultTracker = calc.getLog10PNonRef(vcb.make(), priors);
final long runtime = timer.getElapsedTimeNano(); final long runtime = timer.getElapsedTimeNano();
final List<Object> columns = new LinkedList<Object>(coreValues); final List<Object> columns = new LinkedList<Object>(coreValues);
columns.addAll(Arrays.asList(runtime, result.getnEvaluations(), nonTypePL, position)); columns.addAll(Arrays.asList(runtime, resultTracker.getnEvaluations(), nonTypePL, position));
report.addRowList(columns); report.addRowList(columns);
} }
} }
@ -157,11 +157,11 @@ public class ExactAFCalculationPerformanceTest {
final VariantContext vc = testBuilder.makeACTest(ac, nNonInformative, nonTypePL); final VariantContext vc = testBuilder.makeACTest(ac, nNonInformative, nonTypePL);
timer.start(); timer.start();
final AFCalcResult result = calc.getLog10PNonRef(vc, priors); final AFCalcResultTracker resultTracker = calc.getLog10PNonRef(vc, priors);
final long runtime = timer.getElapsedTimeNano(); final long runtime = timer.getElapsedTimeNano();
final List<Object> columns = new LinkedList<Object>(coreValues); final List<Object> columns = new LinkedList<Object>(coreValues);
columns.addAll(Arrays.asList(runtime, result.getnEvaluations(), nonTypePL, nNonInformative)); columns.addAll(Arrays.asList(runtime, resultTracker.getnEvaluations(), nonTypePL, nNonInformative));
report.addRowList(columns); report.addRowList(columns);
} }
} }
@ -219,9 +219,9 @@ public class ExactAFCalculationPerformanceTest {
final VariantContext vc = testBuilder.makeACTest(new int[]{ac}, 0, 100); final VariantContext vc = testBuilder.makeACTest(new int[]{ac}, 0, 100);
final SimpleTimer timer = new SimpleTimer().start(); 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(); final long runtime = timer.getElapsedTimeNano();
logger.info("result " + result.getNormalizedPosteriorOfAFGTZero()); logger.info("result " + resultTracker.getNormalizedPosteriorOfAFGTZero());
logger.info("runtime " + runtime); logger.info("runtime " + runtime);
} }

View File

@ -78,8 +78,8 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
@Override @Override
public void computeLog10PNonRef(final VariantContext vc, public void computeLog10PNonRef(final VariantContext vc,
final double[] log10AlleleFrequencyPriors, final double[] log10AlleleFrequencyPriors,
final AFCalcResult result) { final AFCalcResultTracker resultTracker) {
combineSinglePools(vc.getGenotypes(), vc.getNAlleles(), ploidy, log10AlleleFrequencyPriors, result); combineSinglePools(vc.getGenotypes(), vc.getNAlleles(), ploidy, log10AlleleFrequencyPriors, resultTracker);
} }
@ -180,13 +180,13 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
* @param numAlleles Number of alternate alleles * @param numAlleles Number of alternate alleles
* @param ploidyPerPool Number of samples per pool * @param ploidyPerPool Number of samples per pool
* @param log10AlleleFrequencyPriors Frequency priors * @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, protected static void combineSinglePools(final GenotypesContext GLs,
final int numAlleles, final int numAlleles,
final int ploidyPerPool, final int ploidyPerPool,
final double[] log10AlleleFrequencyPriors, final double[] log10AlleleFrequencyPriors,
final AFCalcResult result) { final AFCalcResultTracker resultTracker) {
final ArrayList<double[]> genotypeLikelihoods = getGLs(GLs); final ArrayList<double[]> genotypeLikelihoods = getGLs(GLs);
@ -203,9 +203,9 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
combinedPoolLikelihoods.add(set); combinedPoolLikelihoods.add(set);
for (int p=1; p<genotypeLikelihoods.size(); p++) { for (int p=1; p<genotypeLikelihoods.size(); p++) {
result.reset(); // TODO -- why is this here? It makes it hard to track the n evaluation resultTracker.reset(); // TODO -- why is this here? It makes it hard to track the n evaluation
combinedPoolLikelihoods = fastCombineMultiallelicPool(combinedPoolLikelihoods, genotypeLikelihoods.get(p), combinedPloidy, ploidyPerPool, combinedPoolLikelihoods = fastCombineMultiallelicPool(combinedPoolLikelihoods, genotypeLikelihoods.get(p), combinedPloidy, ploidyPerPool,
numAlleles, log10AlleleFrequencyPriors, result); numAlleles, log10AlleleFrequencyPriors, resultTracker);
combinedPloidy = ploidyPerPool + combinedPloidy; // total number of chromosomes in combinedLikelihoods combinedPloidy = ploidyPerPool + combinedPloidy; // total number of chromosomes in combinedLikelihoods
} }
@ -213,7 +213,7 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
public static CombinedPoolLikelihoods fastCombineMultiallelicPool(final CombinedPoolLikelihoods originalPool, double[] newGL, int originalPloidy, int newGLPloidy, int numAlleles, public static CombinedPoolLikelihoods fastCombineMultiallelicPool(final CombinedPoolLikelihoods originalPool, double[] newGL, int originalPloidy, int newGLPloidy, int numAlleles,
final double[] log10AlleleFrequencyPriors, final double[] log10AlleleFrequencyPriors,
final AFCalcResult result) { final AFCalcResultTracker resultTracker) {
@ -235,10 +235,10 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
// keep processing while we have AC conformations that need to be calculated // keep processing while we have AC conformations that need to be calculated
StateTracker stateTracker = new StateTracker(); StateTracker stateTracker = new StateTracker();
while ( !ACqueue.isEmpty() ) { while ( !ACqueue.isEmpty() ) {
result.incNEvaluations(); resultTracker.incNEvaluations();
// compute log10Likelihoods // compute log10Likelihoods
final ExactACset ACset = ACqueue.remove(); final ExactACset ACset = ACqueue.remove();
final double log10LofKs = calculateACConformationAndUpdateQueue(ACset, newPool, originalPool, newGL, log10AlleleFrequencyPriors, originalPloidy, newGLPloidy, result, stateTracker, ACqueue, indexesToACset); final double log10LofKs = calculateACConformationAndUpdateQueue(ACset, newPool, originalPool, newGL, log10AlleleFrequencyPriors, originalPloidy, newGLPloidy, resultTracker, stateTracker, ACqueue, indexesToACset);
// adjust max likelihood seen if needed // adjust max likelihood seen if needed
if ( log10LofKs > stateTracker.getMaxLog10L()) if ( log10LofKs > stateTracker.getMaxLog10L())
@ -263,7 +263,7 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
* @param log10AlleleFrequencyPriors Prior object * @param log10AlleleFrequencyPriors Prior object
* @param originalPloidy Total ploidy of original combined pool * @param originalPloidy Total ploidy of original combined pool
* @param newGLPloidy Ploidy of GL vector * @param newGLPloidy Ploidy of GL vector
* @param result AFResult object * @param resultTracker AFResult object
* @param stateTracker max likelihood observed so far * @param stateTracker max likelihood observed so far
* @param ACqueue Queue of conformations to compute * @param ACqueue Queue of conformations to compute
* @param indexesToACset AC indices of objects in queue * @param indexesToACset AC indices of objects in queue
@ -276,7 +276,7 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
final double[] log10AlleleFrequencyPriors, final double[] log10AlleleFrequencyPriors,
final int originalPloidy, final int originalPloidy,
final int newGLPloidy, final int newGLPloidy,
final AFCalcResult result, final AFCalcResultTracker resultTracker,
final StateTracker stateTracker, final StateTracker stateTracker,
final LinkedList<ExactACset> ACqueue, final LinkedList<ExactACset> ACqueue,
final HashMap<ExactACcounts, ExactACset> indexesToACset) { final HashMap<ExactACcounts, ExactACset> indexesToACset) {
@ -284,7 +284,7 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
// compute likeihood in "set" of new set based on original likelihoods // compute likeihood in "set" of new set based on original likelihoods
final int numAlleles = set.getACcounts().getCounts().length; final int numAlleles = set.getACcounts().getCounts().length;
final int newPloidy = set.getACsum(); 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 // add to new pool
@ -339,11 +339,11 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
* @param ploidy2 Ploidy of second pool * @param ploidy2 Ploidy of second pool
* @param numAlleles Number of alleles * @param numAlleles Number of alleles
* @param log10AlleleFrequencyPriors Array of biallelic priors * @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, public static void combineMultiallelicPoolNaively(CombinedPoolLikelihoods originalPool, double[] yy, int ploidy1, int ploidy2, int numAlleles,
final double[] log10AlleleFrequencyPriors, final double[] log10AlleleFrequencyPriors,
final AFCalcResult result) { final AFCalcResultTracker resultTracker) {
/* /*
final int dim1 = GenotypeLikelihoods.numLikelihoods(numAlleles, ploidy1); final int dim1 = GenotypeLikelihoods.numLikelihoods(numAlleles, ploidy1);
final int dim2 = GenotypeLikelihoods.numLikelihoods(numAlleles, ploidy2); final int dim2 = GenotypeLikelihoods.numLikelihoods(numAlleles, ploidy2);
@ -397,7 +397,7 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
* @param numAlleles Number of alleles (including ref) * @param numAlleles Number of alleles (including ref)
* @param ploidy1 Ploidy of original pool (combined) * @param ploidy1 Ploidy of original pool (combined)
* @param ploidy2 Ploidy of new pool * @param ploidy2 Ploidy of new pool
* @param result AFResult object * @param resultTracker AFResult object
* @return log-likehood of requested conformation * @return log-likehood of requested conformation
*/ */
private static double computeLofK(final ExactACset set, private static double computeLofK(final ExactACset set,
@ -405,7 +405,7 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
final double[] secondGL, final double[] secondGL,
final double[] log10AlleleFrequencyPriors, final double[] log10AlleleFrequencyPriors,
final int numAlleles, final int ploidy1, final int ploidy2, final int numAlleles, final int ploidy1, final int ploidy2,
final AFCalcResult result) { final AFCalcResultTracker resultTracker) {
final int newPloidy = ploidy1 + ploidy2; final int newPloidy = ploidy1 + ploidy2;
@ -423,8 +423,8 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
final double log10Lof0 = firstGLs.getGLOfACZero() + secondGL[HOM_REF_INDEX]; final double log10Lof0 = firstGLs.getGLOfACZero() + secondGL[HOM_REF_INDEX];
set.getLog10Likelihoods()[0] = log10Lof0; set.getLog10Likelihoods()[0] = log10Lof0;
result.setLog10LikelihoodOfAFzero(log10Lof0); resultTracker.setLog10LikelihoodOfAFzero(log10Lof0);
result.setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]); resultTracker.setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]);
return log10Lof0; return log10Lof0;
} else { } else {
@ -467,14 +467,14 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
// update the MLE if necessary // update the MLE if necessary
final int altCounts[] = Arrays.copyOfRange(set.getACcounts().getCounts(),1, set.getACcounts().getCounts().length); 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 // apply the priors over each alternate allele
for (final int ACcount : altCounts ) { for (final int ACcount : altCounts ) {
if ( ACcount > 0 ) if ( ACcount > 0 )
log10LofK += log10AlleleFrequencyPriors[ACcount]; log10LofK += log10AlleleFrequencyPriors[ACcount];
} }
result.updateMAPifNeeded(log10LofK, altCounts); resultTracker.updateMAPifNeeded(log10LofK, altCounts);
return log10LofK; return log10LofK;
} }
@ -506,12 +506,12 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
* @param ploidy1 Ploidy of first pool (# of chromosomes in it) * @param ploidy1 Ploidy of first pool (# of chromosomes in it)
* @param ploidy2 Ploidy of second pool * @param ploidy2 Ploidy of second pool
* @param log10AlleleFrequencyPriors Array of biallelic priors * @param log10AlleleFrequencyPriors Array of biallelic priors
* @param result Af calculation result object * @param resultTracker Af calculation result object
* @return Combined likelihood vector * @return Combined likelihood vector
*/ */
public static ProbabilityVector combineBiallelicPoolsNaively(final ProbabilityVector originalPool, final double[] newPLVector, public static ProbabilityVector combineBiallelicPoolsNaively(final ProbabilityVector originalPool, final double[] newPLVector,
final int ploidy1, final int ploidy2, final double[] log10AlleleFrequencyPriors, final int ploidy1, final int ploidy2, final double[] log10AlleleFrequencyPriors,
final AFCalcResult result) { final AFCalcResultTracker resultTracker) {
final int newPloidy = ploidy1 + ploidy2; final int newPloidy = ploidy1 + ploidy2;
@ -536,8 +536,8 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
final double log10Lof0 = x[0]+y[0]; final double log10Lof0 = x[0]+y[0];
result.setLog10LikelihoodOfAFzero(log10Lof0); resultTracker.setLog10LikelihoodOfAFzero(log10Lof0);
result.setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]); resultTracker.setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]);
double maxElement = log10Lof0; double maxElement = log10Lof0;
int maxElementIdx = 0; int maxElementIdx = 0;
@ -579,8 +579,8 @@ public class GeneralPloidyExactAFCalc extends ExactAFCalc {
} }
alleleCounts[0] = k; alleleCounts[0] = k;
result.updateMLEifNeeded(combinedLikelihoods[k],alleleCounts); resultTracker.updateMLEifNeeded(combinedLikelihoods[k],alleleCounts);
result.updateMAPifNeeded(combinedLikelihoods[k] + log10AlleleFrequencyPriors[k],alleleCounts); resultTracker.updateMAPifNeeded(combinedLikelihoods[k] + log10AlleleFrequencyPriors[k],alleleCounts);
} }

View File

@ -76,11 +76,11 @@ public class ExactAFCalculationModelUnitTest extends BaseTest {
} }
} }
public AFCalcResult execute() { public AFCalcResultTracker execute() {
return getCalc().getLog10PNonRef(getVC(), getPriors()); return getCalc().getLog10PNonRef(getVC(), getPriors());
} }
public AFCalcResult executeRef() { public AFCalcResultTracker executeRef() {
final ExactAFCalc ref = new ReferenceDiploidExactAFCalc(getCalc().nSamples, getCalc().getMaxAltAlleles()); final ExactAFCalc ref = new ReferenceDiploidExactAFCalc(getCalc().nSamples, getCalc().getMaxAltAlleles());
return ref.getLog10PNonRef(getVC(), getPriors()); return ref.getLog10PNonRef(getVC(), getPriors());
} }
@ -209,8 +209,8 @@ public class ExactAFCalculationModelUnitTest extends BaseTest {
@Test(enabled = true, dataProvider = "GLsWithNonInformative", dependsOnMethods = "testGLs") @Test(enabled = true, dataProvider = "GLsWithNonInformative", dependsOnMethods = "testGLs")
public void testGLsWithNonInformative(GetGLsTest onlyInformative, GetGLsTest withNonInformative) { public void testGLsWithNonInformative(GetGLsTest onlyInformative, GetGLsTest withNonInformative) {
final AFCalcResult expected = onlyInformative.execute(); final AFCalcResultTracker expected = onlyInformative.execute();
final AFCalcResult actual = withNonInformative.execute(); final AFCalcResultTracker actual = withNonInformative.execute();
testResultSimple(withNonInformative); testResultSimple(withNonInformative);
@ -225,22 +225,22 @@ public class ExactAFCalculationModelUnitTest extends BaseTest {
} }
private void testResultSimple(final GetGLsTest cfg) { private void testResultSimple(final GetGLsTest cfg) {
final AFCalcResult refResult = cfg.executeRef(); final AFCalcResultTracker refResultTracker = cfg.executeRef();
final AFCalcResult result = cfg.execute(); 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(); // final int minNumberOfEvaluations = cfg.getVC().getCalledChrCount();
// Assert.assertTrue(result.getnEvaluations() >= minNumberOfEvaluations, // Assert.assertTrue(result.getnEvaluations() >= minNumberOfEvaluations,
// "Number of evaluations " + result.getnEvaluations() + " must be at least " + minNumberOfEvaluations); // "Number of evaluations " + result.getnEvaluations() + " must be at least " + minNumberOfEvaluations);
Assert.assertNotNull(result.getAllelesUsedInGenotyping()); Assert.assertNotNull(resultTracker.getAllelesUsedInGenotyping());
Assert.assertTrue(cfg.getAlleles().containsAll(result.getAllelesUsedInGenotyping()), "Result object has alleles not in our initial allele list"); 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++ ) { for ( int altAlleleI = 0; altAlleleI < cfg.numAltAlleles; altAlleleI++ ) {
int expectedAlleleCount = cfg.getExpectedAltAC(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); final Allele allele = cfg.getAlleles().get(altAlleleI+1);
Assert.assertEquals(calcAC_MLE, expectedAlleleCount, "MLE AC not equal to expected AC for allele " + allele); 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, private void compareToRefResult(final AFCalcResultTracker refResultTracker,
final AFCalcResult result) { final AFCalcResultTracker resultTracker) {
final double TOLERANCE = 1; final double TOLERANCE = 1;
// MAP may not be equal // MAP may not be equal
// Assert.assertEquals(result.getAlleleCountsOfMAP(), refResult.getAlleleCountsOfMAP()); // Assert.assertEquals(result.getAlleleCountsOfMAP(), refResult.getAlleleCountsOfMAP());
Assert.assertEquals(result.getAlleleCountsOfMLE(), refResult.getAlleleCountsOfMLE()); Assert.assertEquals(resultTracker.getAlleleCountsOfMLE(), refResultTracker.getAlleleCountsOfMLE());
Assert.assertEquals(result.getAllelesUsedInGenotyping(), refResult.getAllelesUsedInGenotyping()); Assert.assertEquals(resultTracker.getAllelesUsedInGenotyping(), refResultTracker.getAllelesUsedInGenotyping());
Assert.assertEquals(result.getLog10LikelihoodOfAFzero(), refResult.getLog10LikelihoodOfAFzero(), TOLERANCE); Assert.assertEquals(resultTracker.getLog10LikelihoodOfAFzero(), refResultTracker.getLog10LikelihoodOfAFzero(), TOLERANCE);
// Assert.assertEquals(result.getLog10MAP(), refResult.getLog10MAP(), TOLERANCE); // Assert.assertEquals(result.getLog10MAP(), refResult.getLog10MAP(), TOLERANCE);
// Assert.assertEquals(result.getLog10MLE(), refResult.getLog10MLE(), TOLERANCE); // Assert.assertEquals(result.getLog10MLE(), refResult.getLog10MLE(), TOLERANCE);
// Assert.assertEquals(result.getLog10PosteriorOfAFzero(), refResult.getLog10PosteriorOfAFzero(), TOLERANCE); // Assert.assertEquals(result.getLog10PosteriorOfAFzero(), refResult.getLog10PosteriorOfAFzero(), TOLERANCE);
// Assert.assertEquals(result.getLog10PosteriorsMatrixSumWithoutAFzero(), refResult.getLog10PosteriorsMatrixSumWithoutAFzero(), TOLERANCE); // Assert.assertEquals(result.getLog10PosteriorsMatrixSumWithoutAFzero(), refResult.getLog10PosteriorsMatrixSumWithoutAFzero(), TOLERANCE);
Assert.assertEquals(result.getNormalizedPosteriorOfAFGTZero(), refResult.getNormalizedPosteriorOfAFGTZero(), 0.5); Assert.assertEquals(resultTracker.getNormalizedPosteriorOfAFGTZero(), refResultTracker.getNormalizedPosteriorOfAFGTZero(), 0.5);
Assert.assertEquals(result.getNormalizedPosteriorOfAFzero(), refResult.getNormalizedPosteriorOfAFzero(), 0.5); Assert.assertEquals(resultTracker.getNormalizedPosteriorOfAFzero(), refResultTracker.getNormalizedPosteriorOfAFzero(), 0.5);
} }
@Test(enabled = true, dataProvider = "Models") @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); 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"); 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); 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); 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"); 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(resultTracker.getAlleleCountsOfMAP()[0], 1);
Assert.assertEquals(result.getAlleleCountsOfMAP()[1], 1); Assert.assertEquals(resultTracker.getAlleleCountsOfMAP()[1], 1);
} }
// -------------------------------------------------------------------------------- // --------------------------------------------------------------------------------
@ -400,9 +400,9 @@ public class ExactAFCalculationModelUnitTest extends BaseTest {
final VariantContextBuilder vcb = new VariantContextBuilder(vcRoot); final VariantContextBuilder vcb = new VariantContextBuilder(vcRoot);
vcb.genotypes(genotypes); 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"); "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 refPrior = 1 - QualityUtils.qualToErrorProb(log10NonRefPrior);
final double[] priors = MathUtils.toLog10(new double[]{refPrior, (1-refPrior) / 2, (1-refPrior) / 2}); 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); GetGLsTest cfg = new GetGLsTest(model, 1, Arrays.asList(AB), priors, "pNonRef" + log10NonRefPrior);
final AFCalcResult result = cfg.execute(); final AFCalcResultTracker resultTracker = cfg.execute();
final int actualAC = result.getAlleleCountsOfMAP()[0]; final int actualAC = resultTracker.getAlleleCountsOfMAP()[0];
final double pRefWithPrior = AB.getLikelihoods().getAsVector()[0] + priors[0]; final double pRefWithPrior = AB.getLikelihoods().getAsVector()[0] + priors[0];
final double pHetWithPrior = AB.getLikelihoods().getAsVector()[1] + priors[1]; final double pHetWithPrior = AB.getLikelihoods().getAsVector()[1] + priors[1];
final boolean expectNonRef = pRefWithPrior <= pHetWithPrior; final boolean expectNonRef = pRefWithPrior <= pHetWithPrior;
if ( expectNonRef ) if ( expectNonRef )
Assert.assertTrue(result.getNormalizedPosteriorOfAFGTZero() > 0.5); Assert.assertTrue(resultTracker.getNormalizedPosteriorOfAFGTZero() > 0.5);
else else
Assert.assertTrue(result.getNormalizedPosteriorOfAFGTZero() < 0.5); Assert.assertTrue(resultTracker.getNormalizedPosteriorOfAFGTZero() < 0.5);
final int expectedAC = expectNonRef ? 1 : 0; final int expectedAC = expectNonRef ? 1 : 0;
Assert.assertEquals(actualAC, expectedAC, Assert.assertEquals(actualAC, expectedAC,
@ -468,8 +468,8 @@ public class ExactAFCalculationModelUnitTest extends BaseTest {
final double nonRefPrior = (1-refPrior) / 2; final double nonRefPrior = (1-refPrior) / 2;
final double[] priors = MathUtils.toLog10(new double[]{refPrior, nonRefPrior, nonRefPrior, nonRefPrior, nonRefPrior, nonRefPrior}); 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); GetGLsTest cfg = new GetGLsTest(model, 2, Arrays.asList(AB, AC), priors, "pNonRef" + log10NonRefPrior);
final AFCalcResult result = cfg.execute(); final AFCalcResultTracker resultTracker = cfg.execute();
final int actualAC_AB = result.getAlleleCountsOfMAP()[0]; final int actualAC_AB = resultTracker.getAlleleCountsOfMAP()[0];
final double pRefABWithPrior = AB.getLikelihoods().getAsVector()[0] + priors[0]; final double pRefABWithPrior = AB.getLikelihoods().getAsVector()[0] + priors[0];
final double pHetABWithPrior = AB.getLikelihoods().getAsVector()[1] + priors[1]; 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 nonRefPriorSecondAllele = Math.pow(nonRefPrior, 2);
final double refPriorSecondAllele = 1 - nonRefPriorSecondAllele; 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 pRefACWithPrior = AB.getLikelihoods().getAsVector()[0] + Math.log10(refPriorSecondAllele);
final double pHetACWithPrior = AC.getLikelihoods().getAsVector()[3] + Math.log10(nonRefPriorSecondAllele); final double pHetACWithPrior = AC.getLikelihoods().getAsVector()[3] + Math.log10(nonRefPriorSecondAllele);
final int expectedAC_AC = pRefACWithPrior <= pHetACWithPrior ? 1 : 0; final int expectedAC_AC = pRefACWithPrior <= pHetACWithPrior ? 1 : 0;

View File

@ -138,15 +138,15 @@ public class GeneralPloidyAFCalculationModelUnitTest extends BaseTest {
@Test(dataProvider = "getGLs") @Test(dataProvider = "getGLs")
public void testGLs(GetGLsTest cfg) { 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()); final int len = GeneralPloidyGenotypeLikelihoods.getNumLikelihoodElements(1 + cfg.numAltAlleles, cfg.ploidy * cfg.GLs.size());
double[] priors = new double[len]; // flat priors 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; int nameIndex = 1;
for ( int allele = 0; allele < cfg.numAltAlleles; allele++, nameIndex+=2 ) { for ( int allele = 0; allele < cfg.numAltAlleles; allele++, nameIndex+=2 ) {
int expectedAlleleCount = Integer.valueOf(cfg.name.substring(nameIndex, nameIndex+1)); 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); // System.out.format( "%s Expected:%d Calc:%d\n",cfg.toString(),expectedAlleleCount, calculatedAlleleCount);
Assert.assertEquals(calculatedAlleleCount, expectedAlleleCount); Assert.assertEquals(calculatedAlleleCount, expectedAlleleCount);

View File

@ -35,7 +35,7 @@ import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.annotator.VariantAnnotatorEngine; 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.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.*;
import org.broadinstitute.sting.utils.baq.BAQ; import org.broadinstitute.sting.utils.baq.BAQ;
import org.broadinstitute.sting.utils.classloader.PluginManager; import org.broadinstitute.sting.utils.classloader.PluginManager;
@ -82,9 +82,6 @@ public class UnifiedGenotyperEngine {
// the model used for calculating p(non-ref) // the model used for calculating p(non-ref)
private ThreadLocal<AFCalc> afcm = new ThreadLocal<AFCalc>(); private ThreadLocal<AFCalc> afcm = new ThreadLocal<AFCalc>();
// the allele frequency likelihoods and posteriors (allocated once as an optimization)
private ThreadLocal<AFCalcResult> alleleFrequencyCalculationResult = new ThreadLocal<AFCalcResult>();
// because the allele frequency priors are constant for a given i, we cache the results to avoid having to recompute everything // 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[] log10AlleleFrequencyPriorsSNPs;
private final double[] log10AlleleFrequencyPriorsIndels; private final double[] log10AlleleFrequencyPriorsIndels;
@ -355,9 +352,7 @@ public class UnifiedGenotyperEngine {
// initialize the data for this thread if that hasn't been done yet // initialize the data for this thread if that hasn't been done yet
if ( afcm.get() == null ) { if ( afcm.get() == null ) {
afcm.set(getAlleleFrequencyCalculationObject(N, logger, verboseWriter, UAC)); 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 // estimate our confidence in a reference call and return
if ( vc.getNSamples() == 0 ) { if ( vc.getNSamples() == 0 ) {
@ -368,7 +363,7 @@ public class UnifiedGenotyperEngine {
generateEmptyContext(tracker, refContext, stratifiedContexts, rawContext)); 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? // is the most likely frequency conformation AC=0 for all alternate alleles?
boolean bestGuessIsRef = true; boolean bestGuessIsRef = true;
@ -474,7 +469,7 @@ public class UnifiedGenotyperEngine {
// the forward lod // the forward lod
VariantContext vcForward = calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.FORWARD, allAllelesToUse, false, model, perReadAlleleLikelihoodMap); 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[] normalizedLog10Posteriors = MathUtils.normalizeFromLog10(AFresult.log10AlleleFrequencyPosteriors, true);
double forwardLog10PofNull = AFresult.getLog10PosteriorOfAFzero(); double forwardLog10PofNull = AFresult.getLog10PosteriorOfAFzero();
double forwardLog10PofF = AFresult.getLog10PosteriorsMatrixSumWithoutAFzero(); double forwardLog10PofF = AFresult.getLog10PosteriorsMatrixSumWithoutAFzero();
@ -482,7 +477,7 @@ public class UnifiedGenotyperEngine {
// the reverse lod // the reverse lod
VariantContext vcReverse = calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.REVERSE, allAllelesToUse, false, model, perReadAlleleLikelihoodMap); 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); //normalizedLog10Posteriors = MathUtils.normalizeFromLog10(AFresult.log10AlleleFrequencyPosteriors, true);
double reverseLog10PofNull = AFresult.getLog10PosteriorOfAFzero(); double reverseLog10PofNull = AFresult.getLog10PosteriorOfAFzero();
double reverseLog10PofF = AFresult.getLog10PosteriorsMatrixSumWithoutAFzero(); double reverseLog10PofF = AFresult.getLog10PosteriorsMatrixSumWithoutAFzero();
@ -622,8 +617,6 @@ public class UnifiedGenotyperEngine {
AFline.append(i + "/" + N + "\t"); AFline.append(i + "/" + N + "\t");
AFline.append(String.format("%.2f\t", ((float)i)/N)); AFline.append(String.format("%.2f\t", ((float)i)/N));
AFline.append(String.format("%.8f\t", getAlleleFrequencyPriors(model)[i])); 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()); verboseWriter.println(AFline.toString());
} }

View File

@ -73,6 +73,7 @@ public abstract class AFCalc implements Cloneable {
private SimpleTimer callTimer = new SimpleTimer(); private SimpleTimer callTimer = new SimpleTimer();
private PrintStream callReport = null; private PrintStream callReport = null;
private final AFCalcResultTracker resultTracker;
protected AFCalc(final UnifiedArgumentCollection UAC, final int nSamples, final Logger logger, final PrintStream verboseWriter) { 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); 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; this.verboseWriter = verboseWriter;
if ( exactCallsLog != null ) if ( exactCallsLog != null )
initializeOutputFile(exactCallsLog); initializeOutputFile(exactCallsLog);
} this.resultTracker = new AFCalcResultTracker(Math.max(maxAltAlleles, maxAltAllelesForIndels));
/**
* @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()));
} }
/** /**
@ -111,30 +103,27 @@ public abstract class AFCalc implements Cloneable {
* *
* @param vc the VariantContext holding the alleles and sample information * @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 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) * @return result (for programming convenience)
*/ */
public final AFCalcResult getLog10PNonRef(final VariantContext vc, public AFCalcResultTracker getLog10PNonRef(final VariantContext vc, final double[] log10AlleleFrequencyPriors) {
final double[] log10AlleleFrequencyPriors,
final AFCalcResult result) {
if ( vc == null ) throw new IllegalArgumentException("VariantContext cannot be null"); if ( vc == null ) throw new IllegalArgumentException("VariantContext cannot be null");
if ( log10AlleleFrequencyPriors == null ) throw new IllegalArgumentException("priors vector 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 // reset the result, so we can store our new result there
result.reset(); resultTracker.reset();
final VariantContext vcWorking = reduceScope(vc); final VariantContext vcWorking = reduceScope(vc);
callTimer.start(); callTimer.start();
computeLog10PNonRef(vcWorking, log10AlleleFrequencyPriors, result); computeLog10PNonRef(vcWorking, log10AlleleFrequencyPriors, resultTracker);
final long nanoTime = callTimer.getElapsedTimeNano(); final long nanoTime = callTimer.getElapsedTimeNano();
if ( callReport != null ) if ( callReport != null )
printCallInfo(vcWorking, log10AlleleFrequencyPriors, nanoTime, result.getLog10PosteriorOfAFzero()); printCallInfo(vcWorking, log10AlleleFrequencyPriors, nanoTime, resultTracker.getLog10PosteriorOfAFzero());
result.setAllelesUsedInGenotyping(vcWorking.getAlleles()); resultTracker.setAllelesUsedInGenotyping(vcWorking.getAlleles());
return result; return resultTracker;
} }
// --------------------------------------------------------------------------- // ---------------------------------------------------------------------------
@ -163,12 +152,12 @@ public abstract class AFCalc implements Cloneable {
* *
* @param vc variant context with alleles and genotype likelihoods * @param vc variant context with alleles and genotype likelihoods
* @param log10AlleleFrequencyPriors priors * @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 // TODO -- add consistent requires among args
public abstract void computeLog10PNonRef(final VariantContext vc, public abstract void computeLog10PNonRef(final VariantContext vc,
final double[] log10AlleleFrequencyPriors, final double[] log10AlleleFrequencyPriors,
final AFCalcResult result); final AFCalcResultTracker resultTracker);
/** /**
* Must be overridden by concrete subclasses * Must be overridden by concrete subclasses

View File

@ -26,38 +26,36 @@
package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc; package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc;
import com.google.java.contract.Ensures; import com.google.java.contract.Ensures;
import com.google.java.contract.Requires;
import org.broadinstitute.sting.utils.MathUtils; import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.variantcontext.Allele; import org.broadinstitute.sting.utils.variantcontext.Allele;
import java.util.Arrays; import java.util.Arrays;
import java.util.List; import java.util.List;
/** /**
* Created by IntelliJ IDEA. * Describes the results of the AFCalc
* User: ebanks
* Date: Dec 14, 2011
* *
* 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 { 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 private final static int AF0 = 0;
protected double log10MLE; private final static int AF1p = 1;
protected double log10MAP; 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[] 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; int nEvaluations = 0;
@ -68,36 +66,28 @@ public class AFCalcResult {
/** /**
* Create a results object capability of storing results for calls with up to maxAltAlleles * Create a results object capability of storing results for calls with up to maxAltAlleles
*
* @param maxAltAlleles an integer >= 1
*/ */
public AFCalcResult(final int maxAltAlleles) { public AFCalcResult(final int[] alleleCountsOfMLE,
if ( maxAltAlleles < 1 ) throw new IllegalArgumentException("maxAltAlleles must be >= 0, saw " + maxAltAlleles); final int nEvaluations,
final List<Allele> 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]; this.alleleCountsOfMLE = alleleCountsOfMLE;
alleleCountsOfMAP = new int[maxAltAlleles]; this.nEvaluations = nEvaluations;
this.allelesUsedInGenotyping = allelesUsedInGenotyping;
reset(); this.log10LikelihoodsOfAC = Arrays.copyOf(log10LikelihoodsOfAC, LOG_10_ARRAY_SIZES);
} this.log10PriorsOfAC = Arrays.copyOf(log10PriorsOfAC, LOG_10_ARRAY_SIZES);
this.log10PosteriorsOfAC = computePosteriors(log10LikelihoodsOfAC, log10PriorsOfAC);
/**
* 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;
} }
/** /**
@ -115,18 +105,6 @@ public class AFCalcResult {
return alleleCountsOfMLE; 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 * Returns the number of cycles used to evaluate the pNonRef for this AF calculation
* *
@ -136,36 +114,6 @@ public class AFCalcResult {
return nEvaluations; 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. * 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 * Get the log10 normalized -- across all ACs -- posterior probability of AC == 0
* @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 * @return
*/ */
protected void reset() { @Ensures({"goodLog10Value(result)"})
log10MLE = log10MAP = log10LikelihoodOfAFzero = log10PosteriorOfAFzero = AFCalc.VALUE_NOT_CALCULATED; public double getLog10PosteriorOfAFEq0() {
for ( int i = 0; i < alleleCountsOfMLE.length; i++ ) { return log10PosteriorsOfAC[AF0];
alleleCountsOfMLE[i] = 0;
alleleCountsOfMAP[i] = 0;
}
currentPosteriorsCacheIndex = 0;
log10PosteriorMatrixSum = null;
allelesUsedInGenotyping = null;
nEvaluations = 0;
} }
/** /**
* 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() { @Ensures({"goodLog10Value(result)"})
nEvaluations++; public double getLog10PosteriorOfAFGT0() {
return log10PosteriorsOfAC[AF1p];
} }
protected void updateMLEifNeeded(final double log10LofK, final int[] alleleCountsForK) { /**
if ( log10LofK > log10MLE ) { * Get the log10 unnormalized -- across all ACs -- likelihood of AC == 0
log10MLE = log10LofK; *
for ( int i = 0; i < alleleCountsForK.length; i++ ) * @return
alleleCountsOfMLE[i] = alleleCountsForK[i]; */
@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) { if ( shouldSumToOne || MathUtils.compareDoubles(MathUtils.sumLog10(vector), 0.0, 1e-2) != 0 )
addToPosteriorsCache(log10LofK); return false;
if ( log10LofK > log10MAP ) { return true; // everything is good
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<Allele> 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) { private static boolean goodLog10Value(final double result) {
return result <= 0.0 || Double.isInfinite(result) || Double.isNaN(result); return result <= 0.0 && ! Double.isInfinite(result) && ! Double.isNaN(result);
}
protected void setAClimits(int[] AClimits) {
this.AClimits = AClimits;
} }
} }

View File

@ -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<Allele> 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<Allele> 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<Allele> 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;
}
}

View File

@ -19,9 +19,9 @@ public class ConstrainedDiploidExactAFCalc extends DiploidExactAFCalc {
super(UAC, N, logger, verboseWriter); 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); final int[] maxACsToConsider = computeMaxACs(vc);
result.setAClimits(maxACsToConsider); resultTracker.setAClimits(maxACsToConsider);
return new StateTracker(maxACsToConsider); return new StateTracker(maxACsToConsider);
} }

View File

@ -42,12 +42,12 @@ public abstract class DiploidExactAFCalc extends ExactAFCalc {
super(UAC, N, logger, verboseWriter); 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 @Override
public void computeLog10PNonRef(final VariantContext vc, public void computeLog10PNonRef(final VariantContext vc,
final double[] log10AlleleFrequencyPriors, final double[] log10AlleleFrequencyPriors,
final AFCalcResult result) { final AFCalcResultTracker resultTracker) {
final int numAlternateAlleles = vc.getNAlleles() - 1; final int numAlternateAlleles = vc.getNAlleles() - 1;
final ArrayList<double[]> genotypeLikelihoods = getGLs(vc.getGenotypes()); final ArrayList<double[]> genotypeLikelihoods = getGLs(vc.getGenotypes());
final int numSamples = genotypeLikelihoods.size()-1; final int numSamples = genotypeLikelihoods.size()-1;
@ -66,16 +66,16 @@ public abstract class DiploidExactAFCalc extends ExactAFCalc {
indexesToACset.put(zeroSet.getACcounts(), zeroSet); indexesToACset.put(zeroSet.getACcounts(), zeroSet);
// keep processing while we have AC conformations that need to be calculated // 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() ) { while ( !ACqueue.isEmpty() ) {
result.incNEvaluations(); // keep track of the number of evaluations resultTracker.incNEvaluations(); // keep track of the number of evaluations
// compute log10Likelihoods // compute log10Likelihoods
final ExactACset set = ACqueue.remove(); final ExactACset set = ACqueue.remove();
if ( stateTracker.withinMaxACs(set.getACcounts()) ) { 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 // adjust max likelihood seen if needed
stateTracker.update(log10LofKs, set.getACcounts()); stateTracker.update(log10LofKs, set.getACcounts());
@ -161,13 +161,13 @@ public abstract class DiploidExactAFCalc extends ExactAFCalc {
final LinkedList<ExactACset> ACqueue, final LinkedList<ExactACset> ACqueue,
final HashMap<ExactACcounts, ExactACset> indexesToACset, final HashMap<ExactACcounts, ExactACset> indexesToACset,
final double[] log10AlleleFrequencyPriors, final double[] log10AlleleFrequencyPriors,
final AFCalcResult result) { final AFCalcResultTracker resultTracker) {
//if ( DEBUG ) //if ( DEBUG )
// System.out.printf(" *** computing LofK for set=%s%n", set.ACcounts); // System.out.printf(" *** computing LofK for set=%s%n", set.ACcounts);
// compute the log10Likelihoods // compute the log10Likelihoods
computeLofK(set, genotypeLikelihoods, log10AlleleFrequencyPriors, result); computeLofK(set, genotypeLikelihoods, log10AlleleFrequencyPriors, resultTracker);
final double log10LofK = set.getLog10Likelihoods()[set.getLog10Likelihoods().length-1]; 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, private void computeLofK(final ExactACset set,
final ArrayList<double[]> genotypeLikelihoods, final ArrayList<double[]> genotypeLikelihoods,
final double[] log10AlleleFrequencyPriors, final double[] log10AlleleFrequencyPriors,
final AFCalcResult result) { final AFCalcResultTracker resultTracker) {
set.getLog10Likelihoods()[0] = 0.0; // the zero case set.getLog10Likelihoods()[0] = 0.0; // the zero case
final int totalK = set.getACsum(); 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]; set.getLog10Likelihoods()[j] = set.getLog10Likelihoods()[j-1] + genotypeLikelihoods.get(j)[HOM_REF_INDEX];
final double log10Lof0 = set.getLog10Likelihoods()[set.getLog10Likelihoods().length-1]; final double log10Lof0 = set.getLog10Likelihoods()[set.getLog10Likelihoods().length-1];
result.setLog10LikelihoodOfAFzero(log10Lof0); resultTracker.setLog10LikelihoodOfAFzero(log10Lof0);
result.setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]); resultTracker.setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]);
return; return;
} }
@ -284,14 +284,14 @@ public abstract class DiploidExactAFCalc extends ExactAFCalc {
double log10LofK = set.getLog10Likelihoods()[set.getLog10Likelihoods().length-1]; double log10LofK = set.getLog10Likelihoods()[set.getLog10Likelihoods().length-1];
// update the MLE if necessary // update the MLE if necessary
result.updateMLEifNeeded(log10LofK, set.getACcounts().getCounts()); resultTracker.updateMLEifNeeded(log10LofK, set.getACcounts().getCounts());
// apply the priors over each alternate allele // apply the priors over each alternate allele
for ( final int ACcount : set.getACcounts().getCounts() ) { for ( final int ACcount : set.getACcounts().getCounts() ) {
if ( ACcount > 0 ) if ( ACcount > 0 )
log10LofK += log10AlleleFrequencyPriors[ACcount]; log10LofK += log10AlleleFrequencyPriors[ACcount];
} }
result.updateMAPifNeeded(log10LofK, set.getACcounts().getCounts()); resultTracker.updateMAPifNeeded(log10LofK, set.getACcounts().getCounts());
} }
private void pushData(final ExactACset targetSet, private void pushData(final ExactACset targetSet,

View File

@ -52,31 +52,31 @@ public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc {
} }
@Override @Override
protected StateTracker makeMaxLikelihood(VariantContext vc, AFCalcResult result) { protected StateTracker makeMaxLikelihood(VariantContext vc, AFCalcResultTracker resultTracker) {
return refModel.makeMaxLikelihood(vc, result); return refModel.makeMaxLikelihood(vc, resultTracker);
} }
@Override @Override
public void computeLog10PNonRef(final VariantContext vc, public void computeLog10PNonRef(final VariantContext vc,
final double[] log10AlleleFrequencyPriors, final double[] log10AlleleFrequencyPriors,
final AFCalcResult result) { final AFCalcResultTracker resultTracker) {
final List<AFCalcResult> independentResults = computeLog10PNonRefForEachAllele(vc, log10AlleleFrequencyPriors); final List<AFCalcResultTracker> independentResultTrackers = computeLog10PNonRefForEachAllele(vc, log10AlleleFrequencyPriors);
combineIndependentPNonRefs(vc, independentResults, log10AlleleFrequencyPriors, result); combineIndependentPNonRefs(vc, independentResultTrackers, log10AlleleFrequencyPriors, resultTracker);
} }
protected List<AFCalcResult> computeLog10PNonRefForEachAllele(final VariantContext vc, protected List<AFCalcResultTracker> computeLog10PNonRefForEachAllele(final VariantContext vc,
final double[] log10AlleleFrequencyPriors) { final double[] log10AlleleFrequencyPriors) {
final int nAltAlleles = vc.getNAlleles() - 1; final int nAltAlleles = vc.getNAlleles() - 1;
final List<AFCalcResult> results = new ArrayList<AFCalcResult>(nAltAlleles); final List<AFCalcResultTracker> resultTrackers = new ArrayList<AFCalcResultTracker>(nAltAlleles);
for ( int altI = 0; altI < nAltAlleles; altI++ ) { for ( int altI = 0; altI < nAltAlleles; altI++ ) {
final List<Allele> biallelic = Arrays.asList(vc.getReference(), vc.getAlternateAllele(altI)); final List<Allele> biallelic = Arrays.asList(vc.getReference(), vc.getAlternateAllele(altI));
final VariantContext subvc = biallelicCombinedGLs(vc, biallelic, altI + 1); final VariantContext subvc = biallelicCombinedGLs(vc, biallelic, altI + 1);
final AFCalcResult result = refModel.getLog10PNonRef(subvc, log10AlleleFrequencyPriors); final AFCalcResultTracker resultTracker = refModel.getLog10PNonRef(subvc, log10AlleleFrequencyPriors);
results.add(result); resultTrackers.add(resultTracker);
} }
return results; return resultTrackers;
} }
protected VariantContext biallelicCombinedGLs(final VariantContext rootVC, final List<Allele> biallelic, final int allele2) { protected VariantContext biallelicCombinedGLs(final VariantContext rootVC, final List<Allele> 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 * Takes each independent result and merges it into the final result object
* *
* @param independentPNonRefs the pNonRef result for each allele independently * @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, protected void combineIndependentPNonRefs(final VariantContext vc,
final List<AFCalcResult> independentPNonRefs, final List<AFCalcResultTracker> independentPNonRefs,
final double[] log10AlleleFrequencyPriors, final double[] log10AlleleFrequencyPriors,
final AFCalcResult result) { final AFCalcResultTracker resultTracker) {
final int nChrom = vc.getNSamples() * 2; 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 // both the likelihood and the posterior of AF=0 are the same for all alleles
// TODO -- check and ensure this is true // TODO -- check and ensure this is true
result.setLog10LikelihoodOfAFzero(independentPNonRefs.get(0).getLog10LikelihoodOfAFzero()); resultTracker.setLog10LikelihoodOfAFzero(independentPNonRefs.get(0).getLog10LikelihoodOfAFzero());
result.setLog10PosteriorOfAFzero(independentPNonRefs.get(0).getLog10PosteriorOfAFzero()); resultTracker.setLog10PosteriorOfAFzero(independentPNonRefs.get(0).getLog10PosteriorOfAFzero());
result.log10PosteriorMatrixSum = 0.0; resultTracker.log10PosteriorMatrixSum = 0.0;
int altI = 0; int altI = 0;
for ( final AFCalcResult independentPNonRef : independentPNonRefs ) { for ( final AFCalcResultTracker independentPNonRef : independentPNonRefs ) {
result.log10MLE += independentPNonRef.getLog10MLE(); resultTracker.log10MLE += independentPNonRef.getLog10MLE();
// TODO -- technically double counting some posterior mass // TODO -- technically double counting some posterior mass
result.log10MAP += independentPNonRef.getLog10MAP(); resultTracker.log10MAP += independentPNonRef.getLog10MAP();
// TODO -- technically double counting some posterior mass // TODO -- technically double counting some posterior mass
result.log10PosteriorMatrixSum += independentPNonRef.getLog10PosteriorsMatrixSumWithoutAFzero(); resultTracker.log10PosteriorMatrixSum += independentPNonRef.getLog10PosteriorsMatrixSumWithoutAFzero();
result.getAlleleCountsOfMAP()[altI] = independentPNonRef.getAlleleCountsOfMAP()[0]; resultTracker.getAlleleCountsOfMAP()[altI] = independentPNonRef.getAlleleCountsOfMAP()[0];
result.getAlleleCountsOfMLE()[altI] = independentPNonRef.getAlleleCountsOfMLE()[0]; resultTracker.getAlleleCountsOfMLE()[altI] = independentPNonRef.getAlleleCountsOfMLE()[0];
result.nEvaluations += independentPNonRef.nEvaluations; resultTracker.nEvaluations += independentPNonRef.nEvaluations;
altI++; altI++;
} }
} }

View File

@ -15,7 +15,7 @@ public class ReferenceDiploidExactAFCalc extends DiploidExactAFCalc {
super(UAC, N, logger, verboseWriter); 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(); return new StateTracker();
} }
} }

View File

@ -23,7 +23,7 @@
*/ */
package org.broadinstitute.sting.gatk.walkers.validation.validationsiteselector; 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.DiploidExactAFCalc;
import org.broadinstitute.sting.gatk.walkers.genotyper.afcalc.ReferenceDiploidExactAFCalc; import org.broadinstitute.sting.gatk.walkers.genotyper.afcalc.ReferenceDiploidExactAFCalc;
import org.broadinstitute.sting.utils.variantcontext.VariantContext; import org.broadinstitute.sting.utils.variantcontext.VariantContext;
@ -54,10 +54,10 @@ public class GLBasedSampleSelector extends SampleSelector {
flatPriors = new double[1+2*samples.size()]; flatPriors = new double[1+2*samples.size()];
AFCalculator = new ReferenceDiploidExactAFCalc(samples.size(), 4); AFCalculator = new ReferenceDiploidExactAFCalc(samples.size(), 4);
} }
AFCalcResult result = new AFCalcResult(vc.getAlternateAlleles().size()); AFCalcResultTracker resultTracker = new AFCalcResultTracker(vc.getAlternateAlleles().size());
AFCalculator.computeLog10PNonRef(subContext, flatPriors, result); AFCalculator.computeLog10PNonRef(subContext, flatPriors, resultTracker);
// do we want to let this qual go up or down? // do we want to let this qual go up or down?
if ( result.getLog10PosteriorOfAFzero() < referenceLikelihood ) { if ( resultTracker.getLog10PosteriorOfAFzero() < referenceLikelihood ) {
return true; return true;
} }