From a12ed404ce0b2a527748f776b7c08ba6c902427b Mon Sep 17 00:00:00 2001
From: kiran <kiran@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Wed, 10 Jun 2009 08:21:22 +0000
Subject: [PATCH] Changed method name from applyFourBaseDistributionPrior to
 applySecondBaseDistributionPrior.  'Cause that's how I roll.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@960 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../playground/utils/GenotypeLikelihoods.java | 31 ++++++++++++-------
 1 file changed, 20 insertions(+), 11 deletions(-)

diff --git a/java/src/org/broadinstitute/sting/playground/utils/GenotypeLikelihoods.java b/java/src/org/broadinstitute/sting/playground/utils/GenotypeLikelihoods.java
index b0c541acc..d46da5a51 100644
--- a/java/src/org/broadinstitute/sting/playground/utils/GenotypeLikelihoods.java
+++ b/java/src/org/broadinstitute/sting/playground/utils/GenotypeLikelihoods.java
@@ -197,7 +197,7 @@ public class GenotypeLikelihoods {
         this.sort();
     }
 
-    public void applyFourBaseDistributionPrior(String primaryBases, String secondaryBases) {
+    public void applySecondBaseDistributionPrior(String primaryBases, String secondaryBases) {
         for (int genotypeIndex = 0; genotypeIndex < genotypes.length; genotypeIndex++) {
             char firstAllele = genotypes[genotypeIndex].charAt(0);
             char secondAllele = genotypes[genotypeIndex].charAt(1);
@@ -210,24 +210,33 @@ public class GenotypeLikelihoods {
 
             for (int pileupIndex = 0; pileupIndex < primaryBases.length(); pileupIndex++) {
                 char primaryBase = primaryBases.charAt(pileupIndex);
-                char secondaryBase = secondaryBases.charAt(pileupIndex);
 
-                if (primaryBase != firstAllele && primaryBase != secondAllele) {
-                    if (secondaryBase == firstAllele || secondaryBase == secondAllele) {
-                        offIsGenotypic++;
+                if (secondaryBases != null) {
+                    char secondaryBase = secondaryBases.charAt(pileupIndex);
+
+                    if (primaryBase != firstAllele && primaryBase != secondAllele) {
+                        if (secondaryBase == firstAllele || secondaryBase == secondAllele) {
+                            offIsGenotypic++;
+                        }
+                        offTotal++;
+                    } else {
+                        if (secondaryBase == firstAllele || secondaryBase == secondAllele) {
+                            onIsGenotypic++;
+                        }
+                        onTotal++;
                     }
-                    offTotal++;
-                } else {
-                    if (secondaryBase == firstAllele || secondaryBase == secondAllele) {
-                        onIsGenotypic++;
-                    }
-                    onTotal++;
                 }
             }
 
             double offPrior = MathUtils.binomialProbability(offIsGenotypic, offTotal, offNextBestBasePriors.get(genotypes[genotypeIndex]));
             double onPrior = MathUtils.binomialProbability(onIsGenotypic, onTotal, onNextBestBasePriors.get(genotypes[genotypeIndex]));
 
+            //System.out.printf("%c%c [%d %d %f %e] [%d %d %f %e]\n",
+            //                  firstAllele, secondAllele,
+            //                  offIsGenotypic, offTotal, offNextBestBasePriors.get(genotypes[genotypeIndex]), offPrior,
+            //                  onIsGenotypic, onTotal, onNextBestBasePriors.get(genotypes[genotypeIndex]), onPrior);
+            //System.out.printf("%f\n", (new cern.jet.random.Binomial(onTotal, onNextBestBasePriors.get(genotypes[genotypeIndex]), cern.jet.random.engine.RandomEngine.makeDefault())).pdf(onIsGenotypic));
+
             likelihoods[genotypeIndex] += Math.log10(offPrior) + Math.log10(onPrior);
 
             //System.out.println(genotypes[genotypeIndex] + " " + offNextBestBasePriors.get(genotypes[genotypeIndex]) + " " + offIsGenotypic + " " + offTotal + " " + (((double) offIsGenotypic)/((double) offTotal)) + " " + offPrior);