diff --git a/java/src/org/broadinstitute/sting/gatk/refdata/rodVariants.java b/java/src/org/broadinstitute/sting/gatk/refdata/rodVariants.java index 2d458b6af..9240e1d85 100755 --- a/java/src/org/broadinstitute/sting/gatk/refdata/rodVariants.java +++ b/java/src/org/broadinstitute/sting/gatk/refdata/rodVariants.java @@ -5,6 +5,8 @@ import org.broadinstitute.sting.utils.GenomeLoc; import java.io.IOException; public class rodVariants extends BasicReferenceOrderedDatum { + private enum Genotype { AA, AC, AG, AT, CC, CG, CT, GG, GT, TT } + private GenomeLoc loc; private char refBase = 'N'; private int depth; @@ -75,4 +77,40 @@ public class rodVariants extends BasicReferenceOrderedDatum { public float getLodBtnb() { return lodBtnb; } public float[] getGenotypeLikelihoods() { return genotypeLikelihoods; } + + public void adjustLikelihoods(double[] likelihoods) { + for (int likelihoodIndex = 0; likelihoodIndex < likelihoods.length; likelihoodIndex++) { + genotypeLikelihoods[likelihoodIndex] += likelihoods[likelihoodIndex]; + } + + String bestGenotype = "NN"; + double bestLikelihood = Double.NEGATIVE_INFINITY; + double nextBestLikelihood = Double.NEGATIVE_INFINITY; + double refLikelihood = Double.NEGATIVE_INFINITY; + + for (int likelihoodIndex = 0; likelihoodIndex < likelihoods.length; likelihoodIndex++) { + if (genotypeLikelihoods[likelihoodIndex] > bestLikelihood) { + bestLikelihood = genotypeLikelihoods[likelihoodIndex]; + + bestGenotype = Genotype.values()[likelihoodIndex].toString(); + } + } + + for (int likelihoodIndex = 0; likelihoodIndex < likelihoods.length; likelihoodIndex++) { + if (genotypeLikelihoods[likelihoodIndex] > nextBestLikelihood && genotypeLikelihoods[likelihoodIndex] < bestLikelihood) { + nextBestLikelihood = genotypeLikelihoods[likelihoodIndex]; + } + } + + for (int likelihoodIndex = 0; likelihoodIndex < likelihoods.length; likelihoodIndex++) { + if (refBase == Genotype.values()[likelihoodIndex].toString().charAt(0) && + refBase == Genotype.values()[likelihoodIndex].toString().charAt(1)) { + refLikelihood = genotypeLikelihoods[likelihoodIndex]; + } + } + + this.bestGenotype = bestGenotype; + this.lodBtr = (float) (bestLikelihood - refLikelihood); + this.lodBtnb = (float) (bestLikelihood - nextBestLikelihood); + } } diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VariantFiltrationWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VariantFiltrationWalker.java index b289fc2ea..a03673950 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VariantFiltrationWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/variants/VariantFiltrationWalker.java @@ -47,7 +47,7 @@ public class VariantFiltrationWalker extends LocusWalker { if (feature.equalsIgnoreCase("binomial")) { ivf = new IVFBinomialStrand(); } else { throw new StingException(String.format("Cannot understand feature '%s'\n", feature)); } - variant.multiplyLikelihoods(ivf.compute(ref, context)); + variant.adjustLikelihoods(ivf.compute(ref, context)); vwriter.println(variant); }