minor tweaks

diff --git a/bam_mcns.c b/bam_mcns.c
index d98d48030560a42fd65c00bda4f1633bc90ac169..d337e5a290aa76a1f1faa2c8a687db4639962a16 100644 (file)
--- a/bam_mcns.c
+++ b/bam_mcns.c
@@ -264,8 +264,8 @@ static void mc_add_afs(mc_aux_t *ma, double PD, double *f_map, double *p_map)
                        if (ma->afs1[k] > max) max = ma->afs1[k], max_k = k;
                        e += k * ma->afs1[k];
                }
-               *f_map = .5 * max_k / ma->n; *p_map = max;
-               printf(" * %.3lg:%.3lg:%.3lg:%.3lg * ", sum, 1.-.5*max_k/ma->n, max, 1.-.5*e/ma->n);
+               *f_map = .5 * max_k / ma->n; *p_map = max; // e should equal mc_rst_t::f_exp
+//             printf(" * %.3lg:%.3lg:%.3lg:%.3lg * ", sum, 1.-.5*max_k/ma->n, max, 1.-.5*e/ma->n);
        }
 }
 
@@ -297,3 +297,13 @@ int mc_cal(int ref, int *n, const bam_pileup1_t **plp, mc_aux_t *ma, mc_rst_t *r
                mc_add_afs(ma, rst->PD, &rst->f_map, &rst->p_map);
        return tot;
 }
+
+void mc_dump_afs(mc_aux_t *ma)
+{
+       int k;
+       fprintf(stderr, "[afs]");
+       for (k = 0; k <= ma->M; ++k)
+               fprintf(stderr, " %d:%.3lf", k, ma->afs[ma->M - k]);
+       fprintf(stderr, "\n");
+       memset(ma->afs, 0, sizeof(double) * (ma->M + 1));
+}

diff --git a/bam_mcns.h b/bam_mcns.h
index a86fc1ddbd62e575ffeca6d2d59917e0b7e2f9a6..5c7c42026def5133b2468e882785d756dfaee750 100644 (file)
--- a/bam_mcns.h
+++ b/bam_mcns.h
@@ -28,6 +28,7 @@ extern "C" {
        void mc_destroy(mc_aux_t *ma);
        int mc_cal(int ref, int *n, const bam_pileup1_t **plp, mc_aux_t *ma, mc_rst_t *rst, int level);
        int mc_call_gt(const mc_aux_t *ma, double f0, int k);
+       void mc_dump_afs(mc_aux_t *ma);
 
 #ifdef __cplusplus
 }

diff --git a/bam_plcmd.c b/bam_plcmd.c
index ee49cd3cff69afb2c0872978293fd3eaf0e13d9c..0cc2d6b6710dbedc53f89803971e8170233396eb 100644 (file)
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -449,8 +449,15 @@ int bam_pileup(int argc, char *argv[])
  * mpileup *
  ***********/
 
+#define MPLP_VCF   0x1
+#define MPLP_VAR   0x2
+#define MPLP_AFS   0x4
+#define MPLP_AFALL 0x8
+
+#define MPLP_AFS_BLOCK 0x10000
+
 typedef struct {
-       int vcf, max_mq, min_mq, var_only, prior_type, afs;
+       int max_mq, min_mq, prior_type, flag;
        double theta;
        char *reg, *fn_pos;
        faidx_t *fai;
@@ -481,6 +488,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
        const bam_pileup1_t **plp;
        bam_mplp_t iter;
        bam_header_t *h = 0;
+       uint64_t N = 0;
        char *ref;
        khash_t(64) *hash = 0;
        // allocate
@@ -518,7 +526,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
        }
        if (conf->fn_pos) hash = load_pos(conf->fn_pos, h);
        // write the VCF header
-       if (conf->vcf) {
+       if (conf->flag & MPLP_VCF) {
                kstring_t s;
                s.l = s.m = 0; s.s = 0;
                puts("##fileformat=VCFv4.0");
@@ -538,7 +546,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                free(s.s);
        }
        // mpileup
-       if (conf->vcf) {
+       if (conf->flag & MPLP_VCF) {
                ma = mc_init(n);
                mc_init_prior(ma, conf->prior_type, conf->theta);
        }
@@ -556,11 +564,11 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                        if (conf->fai) ref = fai_fetch(conf->fai, h->target_name[tid], &ref_len);
                        ref_tid = tid;
                }
-               if (conf->vcf) {
+               if (conf->flag & MPLP_VCF) {
                        mc_rst_t r;
                        int j, _ref0, depth, rms_q, _ref0b, is_var = 0, qref = 0, level = 2, tot;
                        uint64_t sqr_sum;
-                       if (conf->afs) level = 3;
+                       if (conf->flag & MPLP_AFS) level = 3;
                        _ref0 = _ref0b = (ref && pos < ref_len)? ref[pos] : 'N';
                        _ref0 = bam_nt16_nt4_table[bam_nt16_table[_ref0]];
                        tot = mc_cal(_ref0, n_plp, plp, ma, &r, level);
@@ -572,7 +580,8 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                                qref = (int)(-3.434 * log(q) + .499);
                                if (qref > 99) qref = 99;
                        }
-                       if (conf->var_only && !is_var) continue;
+                       if ((conf->flag & MPLP_VAR) && !is_var) continue;
+                       ++N; // number of processed lines
                        printf("%s\t%d\t.\t%c\t", h->target_name[tid], pos + 1, _ref0b);
                        if (is_var) {
                                if (_ref0 == r.ref) putchar("ACGTN"[r.alt]);
@@ -591,9 +600,13 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                        }
                        rms_q = (int)(sqrt((double)sqr_sum / depth) + .499);
                        printf("DP=%d;MQ=%d", depth, rms_q);
-                       printf(";AF=%.3lg", 1. - r.f_em);
-                       if (conf->afs)
-                               printf(";AF0=%.3lg;AFN=%.3lg;AFE=%.3lg", 1-r.f_naive, 1-r.f_nielsen, 1-r.f_exp);
+                       if (tot) {
+                               printf(";AF=%.3lf", 1. - r.f_em);
+                               if (conf->flag & MPLP_AFALL) {
+                                       printf(";AF0=%.3lf;AFN=%.3lf;AFE=%.3lf", 1-r.f_naive, 1-r.f_nielsen, 1-r.f_exp);
+                                       if (conf->flag & MPLP_AFS) printf(";AFB=%.3lf", 1-r.f_map);
+                               }
+                       }
                        printf("\tGT:GQ:DP");
                        if (tot) {
                                for (i = 0; i < n; ++i) {
@@ -602,6 +615,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                                }
                        } else for (i = 0; i < n; ++i) printf("\t./.:0:0");
                        putchar('\n');
+                       if (N % MPLP_AFS_BLOCK == 0) mc_dump_afs(ma);
                } else {
                        printf("%s\t%d\t%c", h->target_name[tid], pos + 1, (ref && pos < ref_len)? ref[pos] : 'N');
                        for (i = 0; i < n; ++i) {
@@ -623,6 +637,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                        putchar('\n');
                }
        }
+       if ((conf->flag&MPLP_VCF) && (conf->flag&MPLP_AFS)) mc_dump_afs(ma);
        if (hash) { // free the hash table
                khint_t k;
                for (k = kh_begin(hash); k < kh_end(hash); ++k)
@@ -649,7 +664,7 @@ int bam_mpileup(int argc, char *argv[])
        mplp.max_mq = 60;
        mplp.prior_type = MC_PTYPE_FULL;
        mplp.theta = 1e-3;
-       while ((c = getopt(argc, argv, "f:r:l:VvM:q:t:2F")) >= 0) {
+       while ((c = getopt(argc, argv, "vVcF2Sf:r:l:VM:q:t:")) >= 0) {
                switch (c) {
                case 't': mplp.theta = atof(optarg); break;
                case '2': mplp.prior_type = MC_PTYPE_COND2; break;
@@ -659,11 +674,13 @@ int bam_mpileup(int argc, char *argv[])
                        break;
                case 'r': mplp.reg = strdup(optarg); break;
                case 'l': mplp.fn_pos = strdup(optarg); break;
-               case 'v': mplp.var_only = 1; break;
-               case 'V': mplp.vcf = 1; break;
+               case 'V':
+               case 'c': mplp.flag |= MPLP_VCF; break;
+               case 'S': mplp.flag |= MPLP_AFS; break;
+               case 'F': mplp.flag |= MPLP_AFALL; break;
+               case 'v': mplp.flag |= MPLP_VAR; break;
                case 'M': mplp.max_mq = atoi(optarg); break;
                case 'q': mplp.min_mq = atoi(optarg); break;
-               case 'F': mplp.afs = 1; break;
                }
        }
        if (argc == 1) {
@@ -675,8 +692,9 @@ int bam_mpileup(int argc, char *argv[])
                fprintf(stderr, "         -M INT      cap mapping quality at INT [%d]\n", mplp.max_mq);
                fprintf(stderr, "         -q INT      filter out alignment with MQ smaller than INT [%d]\n", mplp.min_mq);
                fprintf(stderr, "         -t FLOAT    scaled mutation rate [%lg]\n", mplp.theta);
-               fprintf(stderr, "         -V          generate VCF output (SNP calling)\n");
+               fprintf(stderr, "         -c          generate VCF output (consensus calling)\n");
                fprintf(stderr, "         -v          show variant sites only\n");
+               fprintf(stderr, "         -S          calculate AFS (slow, to stderr)\n");
                fprintf(stderr, "         -2          conditional prior\n");
                fprintf(stderr, "\n");
                fprintf(stderr, "Notes: Assuming error independency and diploid individuals.\n\n");

diff --git a/bamtk.c b/bamtk.c
index 1de39203b25a4ee526588358c7ce19cd6916b650..ae972d54598ecee1fdd7102e36712192d6ae6178 100644 (file)
--- a/bamtk.c
+++ b/bamtk.c
@@ -9,7 +9,7 @@
 #endif
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.1.8-4 (r631)"
+#define PACKAGE_VERSION "0.1.8-5 (r636)"
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);