X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=pysam.git;a=blobdiff_plain;f=samtools%2Fbcftools%2Fprob1.c.pysam.c;h=5bc52a564e85bb44a09cc4dedb65783b1762bc6b;hp=ac331848f624e466091552f6291ea46e0d239062;hb=ca46ef4ba4a883c57cea62d5bf1bc021f1185109;hpb=d02fe5283ed7a93a2f76a5d6dc6e37b40c11b9b1

diff --git a/samtools/bcftools/prob1.c.pysam.c b/samtools/bcftools/prob1.c.pysam.c
index ac33184..5bc52a5 100644
--- a/samtools/bcftools/prob1.c.pysam.c
+++ b/samtools/bcftools/prob1.c.pysam.c
@@ -499,6 +499,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];
@@ -511,16 +518,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);