X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bcftools%2Fprob1.c;h=37dc8ed32f5af45edc9af6fb8a2d2e306a275e97;hb=955605d6cf50290e32ff923da397021361ceb596;hp=1d0328f936c180ebb447e4006631fa2c018363b6;hpb=018b5036f1db2b928ee92931fb2ca9b81116be71;p=samtools.git

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index 1d0328f..37dc8ed 100644
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -97,6 +97,8 @@ int bcf_p1_read_prior(bcf_p1aux_t *ma, const char *fn)
 	for (k = 0; k <= ma->M; ++k) ma->phi[k] /= sum;
 	for (k = 0; k <= ma->M; ++k) fprintf(stderr, " %d:%.3lg", k, ma->phi[ma->M - k]);
 	fputc('\n', stderr);
+	for (sum = 0., k = 1; k <= ma->M; ++k) sum += k * ma->phi[ma->M - k];
+	fprintf(stderr, "[heterozygosity] %lf\n", (double)sum / ma->M);
 	return 0;
 }
 
@@ -212,7 +214,7 @@ int bcf_p1_call_gt(const bcf_p1aux_t *ma, double f0, int k)
 	}
 	max = 1. - max;
 	if (max < 1e-308) max = 1e-308;
-	q = (int)(-3.434 * log(max) + .499);
+	q = (int)(-4.343 * log(max) + .499);
 	if (q > 99) q = 99;
 	return q<<2|max_i;
 }
@@ -326,18 +328,28 @@ static void contrast(bcf_p1aux_t *ma, double pc[4]) // mc_cal_y() must be called
 		sum_alt += (long double)ma->phi[k] * ma->z[k];
 //	printf("* %lg, %lg *\n", (double)sum, (double)sum_alt); // DEBUG: sum should equal sum_alt
 	sum = sum_alt;
-	y = lgamma(2*n1 + 1) - lgamma(ma->M + 1);
-	for (k = 1, x = 0.; k <= 2 * n2; ++k)
-		x += ma->phi[k] * ma->z2[k] * exp(lgamma(ma->M - k + 1) - lgamma(2*n2 - k + 1) + y);
-	pc[0] = ma->z1[0] * x / sum;
+	// the variant is specific to group2
+//	printf("%lg %lg %lg %lg\n", ma->z[2*(n1+n2)]/exp(ma->t - (ma->t1 + ma->t2)), ma->z1[2*n1], ma->z2[2*n2], (double)sum);
 	y = lgamma(2*n2 + 1) - lgamma(ma->M + 1);
+	for (k = 0, x = 0.; k < 2 * n2; ++k)
+		x += ma->phi[2*n1+k] * ma->z2[k] * expl(lgamma(2*n1 + k + 1) - lgamma(k + 1) + y);
+	pc[1] = ma->z1[2*n1] * x / sum;
+	for (k = 1, x = 0.; k <= 2 * n2; ++k)
+		x += ma->phi[k] * ma->z2[k] * expl(lgamma(ma->M - k + 1) - lgamma(2*n2 - k + 1) + y);
+	pc[1] += ma->z1[0] * x / sum;
+	// the variant is specific to group1
+	y = lgamma(2*n1 + 1) - lgamma(ma->M + 1);
+	for (k = 0, x = 0.; k < 2 * n1; ++k)
+		x += ma->phi[2*n2+k] * ma->z1[k] * expl(lgamma(2*n2 + k + 1) - lgamma(k + 1) + y);
+	pc[0] = ma->z2[2*n2] * x / sum;
 	for (k = 1, x = 0.; k <= 2 * n1; ++k)
-		x += ma->phi[k] * ma->z1[k] * exp(lgamma(ma->M - k + 1) - lgamma(2*n1 - k + 1) + y);
-	pc[1] = ma->z2[0] * x / sum;
-	for (k = 0; k < 4; ++k) {
+		x += ma->phi[k] * ma->z1[k] * expl(lgamma(ma->M - k + 1) - lgamma(2*n1 - k + 1) + y);
+	pc[0] += ma->z2[0] * x / sum;
+	// rescale
+	for (k = 2; k < 4; ++k) {
 		y = 1. - pc[k];
 		if (y <= 0.) y = 1e-100;
-		pc[k] = (int)(-3.434 * log(y) + .499);
+		pc[k] = (int)(-4.343 * log(y) + .499);
 		if (pc[k] > 99.) pc[k] = 99.;
 	}
 }