Imported Upstream version 0.1.17

[samtools.git] / bcftools / prob1.c

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index fc9cb29911c1abc746e63b9001a23e853aaa58f5..a380484310ccd7dc0ee72f922ba9f97aef6be80a 100644 (file)
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -497,6 +497,13 @@ int bcf_p1_cal(const bcf1_t *b, int do_contrast, bcf_p1aux_t *ma, bcf_p1rst_t *r
        for (k = 0, sum = 0.; k < ma->M; ++k)
                sum += ma->afs1[k];
        rst->p_var = (double)sum;
+       { // compute the allele count
+               double max = -1;
+               rst->ac = -1;
+               for (k = 0; k <= ma->M; ++k)
+                       if (max < ma->z[k]) max = ma->z[k], rst->ac = k;
+               rst->ac = ma->M - rst->ac;
+       }
        // calculate f_flat and f_em
        for (k = 0, sum = 0.; k <= ma->M; ++k)
                sum += (long double)ma->z[k];
@@ -509,16 +516,27 @@ int bcf_p1_cal(const bcf1_t *b, int do_contrast, bcf_p1aux_t *ma, bcf_p1rst_t *r
        { // estimate equal-tail credible interval (95% level)
                int l, h;
                double p;
-               for (i = 0, p = 0.; i < ma->M; ++i)
+               for (i = 0, p = 0.; i <= ma->M; ++i)
                        if (p + ma->afs1[i] > 0.025) break;
                        else p += ma->afs1[i];
                l = i;
-               for (i = ma->M-1, p = 0.; i >= 0; --i)
+               for (i = ma->M, p = 0.; i >= 0; --i)
                        if (p + ma->afs1[i] > 0.025) break;
                        else p += ma->afs1[i];
                h = i;
                rst->cil = (double)(ma->M - h) / ma->M; rst->cih = (double)(ma->M - l) / ma->M;
        }
+       if (ma->n1 > 0) { // compute LRT
+               double max0, max1, max2;
+               for (k = 0, max0 = -1; k <= ma->M; ++k)
+                       if (max0 < ma->z[k]) max0 = ma->z[k];
+               for (k = 0, max1 = -1; k <= ma->n1 * 2; ++k)
+                       if (max1 < ma->z1[k]) max1 = ma->z1[k];
+               for (k = 0, max2 = -1; k <= ma->M - ma->n1 * 2; ++k)
+                       if (max2 < ma->z2[k]) max2 = ma->z2[k];
+               rst->lrt = log(max1 * max2 / max0);
+               rst->lrt = rst->lrt < 0? 1 : kf_gammaq(.5, rst->lrt);
+       } else rst->lrt = -1.0;
        rst->cmp[0] = rst->cmp[1] = rst->cmp[2] = rst->p_chi2 = -1.0;
        if (do_contrast && rst->p_var > 0.5) // skip contrast2() if the locus is a strong non-variant
                rst->p_chi2 = contrast2(ma, rst->cmp);