From 038cbcf80e69c8a43140c4fc9c16df1a842a4b6e Mon Sep 17 00:00:00 2001 From: kiran Date: Thu, 23 Jul 2009 20:59:52 +0000 Subject: [PATCH] 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 --- .../gatk/walkers/variants/IVFSecondaryBases.java | 5 ++++- .../sting/playground/utils/GenotypeLikelihoods.java | 9 ++++++++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/IVFSecondaryBases.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/IVFSecondaryBases.java index 283745bc9..be8c03340 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/IVFSecondaryBases.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/IVFSecondaryBases.java @@ -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; diff --git a/java/src/org/broadinstitute/sting/playground/utils/GenotypeLikelihoods.java b/java/src/org/broadinstitute/sting/playground/utils/GenotypeLikelihoods.java index 0ab42def3..ae183ffc8 100644 --- a/java/src/org/broadinstitute/sting/playground/utils/GenotypeLikelihoods.java +++ b/java/src/org/broadinstitute/sting/playground/utils/GenotypeLikelihoods.java @@ -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(); }