From e6091454768385d722644774042e822658d20713 Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@live.co.uk>
Date: Wed, 10 Nov 2010 03:05:46 +0000
Subject: [PATCH] bcftools: compute equal-tail (Bayesian) credible interval

---
 bcftools/call1.c |  3 ++-
 bcftools/prob1.c | 13 +++++++++++++
 bcftools/prob1.h |  1 +
 3 files changed, 16 insertions(+), 1 deletion(-)

diff --git a/bcftools/call1.c b/bcftools/call1.c
index 2b28452..5a5c1cd 100644
--- 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 04a0e15..176a0fc 100644
--- 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 7158fe2..6d9d28e 100644
--- 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;
-- 
2.39.2