If the result from the secondary-base test is 0.0, replace the result with a minimum likelihood such that the log-likelihood doesn't underflow.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1303 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
kiran 2009-07-23 20:59:52 +00:00
parent 093550a3f2
commit 038cbcf80e
2 changed files with 12 additions and 2 deletions

View File

@ -75,7 +75,10 @@ public class IVFSecondaryBases implements IndependentVariantFeature {
double offPrior = MathUtils.binomialProbability(offIsGenotypic, offTotal, p2off[genotypeIndex]);
double onPrior = MathUtils.binomialProbability(onIsGenotypic, onTotal, p2on[genotypeIndex]);
likelihoods[genotypeIndex] = Math.log10(offPrior) + Math.log10(onPrior);
double logOffPrior = MathUtils.compareDoubles(offPrior, 0.0, 1e-10) == 0 ? Math.log10(Double.MIN_VALUE) : Math.log10(offPrior);
double logOnPrior = MathUtils.compareDoubles(onPrior, 0.0, 1e-10) == 0 ? Math.log10(Double.MIN_VALUE) : Math.log10(onPrior);
likelihoods[genotypeIndex] += logOffPrior + logOnPrior;
}
return likelihoods;

View File

@ -173,6 +173,10 @@ public class GenotypeLikelihoods {
double likelihood = calculateAlleleLikelihood(ref, read, genotypes[i], qual);
likelihoods[i] += likelihood;
coverage += 1;
// if (i == 8) {
// System.out.println(genotypes[8] + " " + likelihood + " " + likelihoods[8] + " " + coverage);
// }
}
}
@ -325,7 +329,10 @@ public class GenotypeLikelihoods {
double offPrior = MathUtils.binomialProbability(offIsGenotypic, offTotal, offNextBestBasePriors.get(genotypes[genotypeIndex]));
double onPrior = MathUtils.binomialProbability(onIsGenotypic, onTotal, onNextBestBasePriors.get(genotypes[genotypeIndex]));
likelihoods[genotypeIndex] += Math.log10(offPrior) + Math.log10(onPrior);
double logOffPrior = MathUtils.compareDoubles(offPrior, 0.0, 1e-10) == 0 ? Math.log10(Double.MIN_VALUE) : Math.log10(offPrior);
double logOnPrior = MathUtils.compareDoubles(onPrior, 0.0, 1e-10) == 0 ? Math.log10(Double.MIN_VALUE) : Math.log10(onPrior);
likelihoods[genotypeIndex] += logOffPrior + logOnPrior;
}
this.sort();
}