X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=samtools.git;a=blobdiff_plain;f=bcftools%2Fprob1.c;fp=bcftools%2Fprob1.c;h=a380484310ccd7dc0ee72f922ba9f97aef6be80a;hp=fc9cb29911c1abc746e63b9001a23e853aaa58f5;hb=e2bd290b8643f0a728c282a8b825392307750210;hpb=9f4bebab2e0917c676ae739b2d05cb22ad6c4ed5

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index fc9cb29..a380484 100644
--- 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);