* 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)

diff --git a/bam_mcns.c b/bam_mcns.c
index f2b284a719074e8a7451668a49b6dd4e3d4331b9..893d2efdcd7a26f6202e0d235f2c263df1cb1772 100644 (file)
--- 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 c70b41d34c793afe23b57b63d9f33ba892f01617..7b64b8f3d0c79fb928facf6cc73cd69c94726de6 100644 (file)
--- 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[]);