bcftools: compute equal-tail (Bayesian) credible interval

diff --git a/bcftools/call1.c b/bcftools/call1.c
index 2b284527185c0e12747ca1dfb5387aa0bd615c36..5a5c1cd1aba3e1c7e15a00638ee468186d052a52 100644 (file)
--- a/bcftools/call1.c
+++ b/bcftools/call1.c
@@ -161,7 +161,8 @@ static int update_bcf1(int n_smpl, bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p
        kputs(b->alt, &s); kputc('\0', &s); kputc('\0', &s);
        kputs(b->info, &s);
        if (b->info[0]) kputc(';', &s);
-       ksprintf(&s, "AF1=%.3lf;AFE=%.3lf", 1.-pr->f_em, 1.-pr->f_exp);
+//     ksprintf(&s, "AF1=%.4lg;AFE=%.4lg;CI95=%.4lg,%.4lg", 1.-pr->f_em, 1.-pr->f_exp, pr->cil, pr->cih);
+       ksprintf(&s, "AF1=%.4lg;CI95=%.4lg,%.4lg", 1.-pr->f_em, pr->cil, pr->cih);
        ksprintf(&s, ";DP4=%d,%d,%d,%d;MQ=%d", a.d[0], a.d[1], a.d[2], a.d[3], a.mq);
        if (a.is_tested) {
                if (pr->pc[0] >= 0.) ksprintf(&s, ";PC4=%lg,%lg,%lg,%lg", pr->pc[0], pr->pc[1], pr->pc[2], pr->pc[3]);

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index 04a0e1555a16832d9237c1a47b9168931717997a..176a0fc96d4e0b12f14172bc7c0d68297e01a50a 100644 (file)
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -391,6 +391,19 @@ int bcf_p1_cal(bcf1_t *b, bcf_p1aux_t *ma, bcf_p1rst_t *rst)
                        flast = rst->f_em;
                }
        }
+       { // estimate equal-tail credible interval (95% level)
+               int l, h;
+               double p;
+               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)
+                       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;
+       }
        rst->g[0] = rst->g[1] = rst->g[2] = -1.;
        contrast(ma, rst->pc);
        return 0;

diff --git a/bcftools/prob1.h b/bcftools/prob1.h
index 7158fe29f89bee73dda2f2bdc9426af038a363c6..6d9d28e3cee6c7d0a66efbd5314912ac511cd384 100644 (file)
--- a/bcftools/prob1.h
+++ b/bcftools/prob1.h
@@ -9,6 +9,7 @@ typedef struct __bcf_p1aux_t bcf_p1aux_t;
 typedef struct {
        int rank0;
        double f_em, f_exp, f_flat, p_ref;
+       double cil, cih;
        double pc[4];
        double g[3];
 } bcf_p1rst_t;