From d742c10178d8981ce062cb1f993149cfe9876613 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sat, 24 Jul 2010 18:07:27 +0000 Subject: [PATCH] minor tweaks --- bam_mcns.c | 14 ++++++++++++-- bam_mcns.h | 1 + bam_plcmd.c | 46 ++++++++++++++++++++++++++++++++-------------- bamtk.c | 2 +- 4 files changed, 46 insertions(+), 17 deletions(-) diff --git a/bam_mcns.c b/bam_mcns.c index d98d480..d337e5a 100644 --- 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 a86fc1d..5c7c420 100644 --- 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 ee49cd3..0cc2d6b 100644 --- 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 1de3920..ae972d5 100644 --- 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[]); -- 2.39.2