From 084a27854cc562b0ec2c649851a405fe894e860a Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@live.co.uk>
Date: Wed, 28 Jul 2010 15:59:19 +0000
Subject: [PATCH]  * samtools-0.1.8-8 (r644)  * mpileup becomes a little stable
 again  * the method is slightly different, but is more theoretically correct 
 * snp calling is O(n^2) instead of O(n^3)

---
 bam_mcns.c | 29 ++++-------------------------
 bamtk.c    |  2 +-
 2 files changed, 5 insertions(+), 26 deletions(-)

diff --git a/bam_mcns.c b/bam_mcns.c
index f2b284a..893d2ef 100644
--- a/bam_mcns.c
+++ b/bam_mcns.c
@@ -14,31 +14,12 @@ struct __mc_aux_t {
 	int n, M;
 	int ref, alt, alt2;
 	double *q2p, *pdg; // pdg -> P(D|g)
-	double *phi, *alpha, *CMk; // CMk=\binom{M}{k}
+	double *phi, *CMk; // CMk=\binom{M}{k}
 	double *z, *zswap; // aux for afs
 	double *afs, *afs1; // afs: accumulative AFS; afs1: site posterior distribution
 	int *qsum, *bcnt;
 };
 
-static void precal_alpha(mc_aux_t *ma) // \alpha[k]=\binom{M}{k}\sum_l\phi_l(l/M)^k(1-l/M)^{M-k}
-{
-	int k, l;
-	memset(ma->alpha, 0, sizeof(double) * (ma->M + 1));
-	for (l = 0; l <= ma->M; ++l)
-		ma->alpha[0] += ma->phi[l] * pow(1. - (double)l / ma->M, ma->M);
-	for (k = 1; k < ma->M; ++k) {
-		for (l = 1; l < ma->M; ++l) { // for k!=0 and k!=ma->M, l=0 and l=ma->M leads to zero
-			double x = exp((lgamma(ma->M + 1) - lgamma(k + 1) - lgamma(ma->M - k + 1))
-						   + k * log((double)l / ma->M)
-						   + (ma->M - k) * log(1. - (double)l / ma->M));
-			ma->alpha[k] += x * ma->phi[l];
-		}
-	}
-	for (l = 0; l <= ma->M; ++l)
-		ma->alpha[ma->M] += ma->phi[l] * pow((double)l / ma->M, ma->M);
-	fflush(stdout);
-}
-
 void mc_init_prior(mc_aux_t *ma, int type, double theta)
 {
 	int i;
@@ -54,7 +35,6 @@ void mc_init_prior(mc_aux_t *ma, int type, double theta)
 			sum += (ma->phi[i] = theta / (2 * ma->n - i));
 		ma->phi[2 * ma->n] = 1. - sum;
 	}
-	precal_alpha(ma);
 }
 
 mc_aux_t *mc_init(int n) // FIXME: assuming diploid
@@ -68,7 +48,6 @@ mc_aux_t *mc_init(int n) // FIXME: assuming diploid
 	ma->bcnt = calloc(4 * ma->n, sizeof(int));
 	ma->pdg = calloc(3 * ma->n, sizeof(double));
 	ma->phi = calloc(ma->M + 1, sizeof(double));
-	ma->alpha = calloc(ma->M + 1, sizeof(double));
 	ma->CMk = calloc(ma->M + 1, sizeof(double));
 	ma->z = calloc(2 * ma->n + 1, sizeof(double));
 	ma->zswap = calloc(2 * ma->n + 1, sizeof(double));
@@ -87,7 +66,7 @@ void mc_destroy(mc_aux_t *ma)
 	if (ma) {
 		free(ma->qsum); free(ma->bcnt);
 		free(ma->q2p); free(ma->pdg);
-		free(ma->phi); free(ma->alpha); free(ma->CMk);
+		free(ma->phi); free(ma->CMk);
 		free(ma->z); free(ma->zswap);
 		free(ma->afs); free(ma->afs1);
 		free(ma);
@@ -247,9 +226,9 @@ static double mc_add_afs(mc_aux_t *ma)
 	memset(ma->afs1, 0, sizeof(double) * (ma->M + 1));
 	mc_cal_z(ma);
 	for (k = 0, sum = 0.; k <= ma->M; ++k)
-		sum += (long double)ma->alpha[k] * ma->z[k] / ma->CMk[k];
+		sum += (long double)ma->phi[k] * ma->z[k] / ma->CMk[k];
 	for (k = 0; k <= ma->M; ++k) {
-		ma->afs1[k] = ma->alpha[k] * ma->z[k] / ma->CMk[k] / sum;
+		ma->afs1[k] = ma->phi[k] * ma->z[k] / ma->CMk[k] / sum;
 		if (isnan(ma->afs1[k]) || isinf(ma->afs1[k])) return -1.;
 	}
 	for (k = 0, sum = 0.; k <= ma->M; ++k) {
diff --git a/bamtk.c b/bamtk.c
index c70b41d..7b64b8f 100644
--- a/bamtk.c
+++ b/bamtk.c
@@ -9,7 +9,7 @@
 #endif
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.1.8-7 (r639)"
+#define PACKAGE_VERSION "0.1.8-8 (r644)"
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);
-- 
2.39.5