From 7a921c908cf6d9796d1989f1e5cfc3f0f59db497 Mon Sep 17 00:00:00 2001
From: kiran <kiran@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Thu, 18 Jun 2009 07:26:37 +0000
Subject: [PATCH] Can now adjust the genotype likelihoods of a variant returned
 from the rod. This automatically causes the lodBtr, lodBtnb, and genotype to
 be recomputed.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1041 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../sting/gatk/refdata/rodVariants.java       | 38 +++++++++++++++++++
 .../variants/VariantFiltrationWalker.java     |  2 +-
 2 files changed, 39 insertions(+), 1 deletion(-)

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<Integer, Integer> {
             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);
         }