X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=pysam.git;a=blobdiff_plain;f=samtools%2Fbcftools%2Fem.c.pysam.c;h=914dc0587b8bc83486ea7b7232347a1bf87b4905;hp=94b6d7532ca012c698986b7f668e198331b211a5;hb=e1756c41e7a1d7cc01fb95e42df9dd04da2d2991;hpb=ca46ef4ba4a883c57cea62d5bf1bc021f1185109

diff --git a/samtools/bcftools/em.c.pysam.c b/samtools/bcftools/em.c.pysam.c
index 94b6d75..914dc05 100644
--- a/samtools/bcftools/em.c.pysam.c
+++ b/samtools/bcftools/em.c.pysam.c
@@ -170,7 +170,7 @@ static double lk_ratio_test(int n, int n1, const double *pdg, double f3[3][3])
 // x[5..6]: group1 freq, group2 freq
 // x[7]: 1-degree P-value
 // x[8]: 2-degree P-value
-int bcf_em1(const bcf1_t *b, int n1, int flag, double x[10])
+int bcf_em1(const bcf1_t *b, int n1, int flag, double x[9])
 {
 	double *pdg;
 	int i, n, n2;
@@ -182,7 +182,7 @@ int bcf_em1(const bcf1_t *b, int n1, int flag, double x[10])
 	n = b->n_smpl; n2 = n - n1;
 	pdg = get_pdg3(b);
 	if (pdg == 0) return -1;
-	for (i = 0; i < 10; ++i) x[i] = -1.; // set to negative
+	for (i = 0; i < 9; ++i) x[i] = -1.;
 	{
 		if ((x[0] = est_freq(n, pdg)) < 0.) {
 			free(pdg);
@@ -190,7 +190,7 @@ int bcf_em1(const bcf1_t *b, int n1, int flag, double x[10])
 		}
 		x[0] = freqml(x[0], 0, n, pdg);
 	}
-	if (flag & (0xf<<1|3<<8)) { // estimate the genotype frequency and test HWE
+	if (flag & (0xf<<1|1<<8)) { // estimate the genotype frequency and test HWE
 		double *g = x + 1, f3[3], r;
 		f3[0] = g[0] = (1 - x[0]) * (1 - x[0]);
 		f3[1] = g[1] = 2 * x[0] * (1 - x[0]);
@@ -209,24 +209,20 @@ int bcf_em1(const bcf1_t *b, int n1, int flag, double x[10])
 		x[6] = freqml(x[0], n1, n, pdg);
 	}
 	if ((flag & 1<<7) && n1 > 0 && n1 < n) { // 1-degree P-value
-		double f[3], f3[3][3], tmp;
+		double f[3], f3[3][3];
 		f[0] = x[0]; f[1] = x[5]; f[2] = x[6];
 		for (i = 0; i < 3; ++i)
 			f3[i][0] = (1-f[i])*(1-f[i]), f3[i][1] = 2*f[i]*(1-f[i]), f3[i][2] = f[i]*f[i];
-		tmp = log(lk_ratio_test(n, n1, pdg, f3));
-		if (tmp < 0) tmp = 0;
-		x[7] = kf_gammaq(.5, tmp);
+		x[7] = kf_gammaq(.5, log(lk_ratio_test(n, n1, pdg, f3)));
 	}
-	if ((flag & 3<<8) && n1 > 0 && n1 < n) { // 2-degree P-value
-		double g[3][3], tmp;
+	if ((flag & 1<<8) && n1 > 0 && n1 < n) { // 2-degree P-value
+		double g[3][3];
 		for (i = 0; i < 3; ++i) memcpy(g[i], x + 1, 3 * sizeof(double));
 		for (i = 0; i < ITER_MAX; ++i)
 			if (g3_iter(g[1], pdg, 0, n1) < EPS) break;
 		for (i = 0; i < ITER_MAX; ++i)
 			if (g3_iter(g[2], pdg, n1, n) < EPS) break;
-		tmp = log(lk_ratio_test(n, n1, pdg, g));
-		if (tmp < 0) tmp = 0;
-		x[8] = kf_gammaq(1., tmp);
+		x[8] = kf_gammaq(1., log(lk_ratio_test(n, n1, pdg, g)));
 	}
 	// free
 	free(pdg);