Fix P-value calculation

src/main/java/com/antigenomics/mageri/core/output/VcfUtil.java

@@ -21,7 +21,7 @@
 import com.antigenomics.mageri.core.mutations.Mutation;
 
 public final class VcfUtil {
-    public static final int MAX_QUAL = 999, UNDEF_QUAL = -1;
+    public static final int MAX_QUAL = 100, UNDEF_QUAL = -1;
 
     private VcfUtil() {
     }

src/main/java/com/antigenomics/mageri/misc/AuxiliaryStats.java

@@ -38,33 +38,34 @@ public static double negativeBinomialPdf(int k, int r, double p) {
 
     public static double betaBinomialPdf(int k, int n, double alpha, double beta) {
         return Math.exp(
-                Beta.logBeta(k + alpha, n - k + beta) - Beta.logBeta(alpha, beta) +
-                Gamma.logGamma(n + 1) - Gamma.logGamma(k + 1) - Gamma.logGamma(n - k + 1)
+                (Beta.logBeta(k + alpha, n - k + beta) - Beta.logBeta(alpha, beta)) +
+                        (Gamma.logGamma(n + 1) - Gamma.logGamma(k + 1) - Gamma.logGamma(n - k + 1))
         );
     }
 
+    @Deprecated
     public static double betaBinomialCdf(int k, int n, double alpha, double beta) {
         double sum = 0;
 
         for (int i = 0; i < k; i++) {
             sum += betaBinomialPdf(i, n, alpha, beta);
         }
 
-        return Math.min(1.0, sum); // overflow possible ?
+        return Math.min(1.0, sum);
     }
 
     public static double betaBinomialPvalueFast(int k, int n, double alpha, double beta) {
-        return betaBinomialPvalueFast(k, n, alpha, beta, 1e-100);
+        return betaBinomialPvalueFast(k, n, alpha, beta, 1e-10);
     }
 
     public static double betaBinomialPvalueFast(int k, int n, double alpha, double beta, double pThreshold) {
-        double sum = 1 + 0.5 * betaBinomialPdf(k, n, alpha, beta);
+        double sum = 1.0 + 0.5 * betaBinomialPdf(k, n, alpha, beta);
 
         for (int i = 0; i < k; i++) {
             sum -= betaBinomialPdf(i, n, alpha, beta);
 
             if (sum <= pThreshold) {
-                break;
+                return pThreshold;
             }
         }
 

src/test/java/com/antigenomics/mageri/core/variant/model/ErrorModelStatTest.java

@@ -0,0 +1,35 @@
+/*
+* Copyright 2014-2016 Mikhail Shugay
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+
+package com.antigenomics.mageri.core.variant.model;
+
+import com.antigenomics.mageri.misc.AuxiliaryStats;
+import org.junit.Test;
+
+public class ErrorModelStatTest {
+    @Test
+    public void betaBinomTest() {
+        //System.out.println(calc(1447, 1447, 2.1314828, 24344.82));
+        System.out.println(calc(13, 2374, 2.1314828, 24344.82));
+        System.out.println(calc(500, 2374, 2.1314828, 24344.82));
+        //System.out.println(calc(1590, 2374, 2.1314828, 24344.82));
+    }
+
+    private static double calc(int n, int N, double a, double b) {
+        //return -10 * Math.log10(1.0 - AuxiliaryStats.betaBinomialCdf(n, N, a, b));
+        return -10 * Math.log10(AuxiliaryStats.betaBinomialPvalueFast(n, N, a, b));
+    }
+}