Imported Upstream version 0.5

[pysam.git] / samtools / bcftools / prob1.c.pysam.c

diff --git a/samtools/bcftools/prob1.c.pysam.c b/samtools/bcftools/prob1.c.pysam.c
index 5bc52a564e85bb44a09cc4dedb65783b1762bc6b..ac331848f624e466091552f6291ea46e0d239062 100644 (file)
--- a/samtools/bcftools/prob1.c.pysam.c
+++ b/samtools/bcftools/prob1.c.pysam.c
@@ -499,13 +499,6 @@ 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];
@@ -518,27 +511,16 @@ 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, p = 0.; i >= 0; --i)
+               for (i = ma->M-1, 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);