From 183e4cb41c6136b9bff3dc83fcfa9f023e1d7629 Mon Sep 17 00:00:00 2001 From: Petr Danecek Date: Fri, 25 May 2012 14:16:25 +0100 Subject: [PATCH] The -E switch to recalculate extended BAQ on the fly, ignoring existing BQs --- bam_md.c | 7 ++++++- bam_plcmd.c | 10 +++++----- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/bam_md.c b/bam_md.c index d42aa8f..ce40a12 100644 --- a/bam_md.c +++ b/bam_md.c @@ -188,7 +188,7 @@ int bam_cap_mapQ(bam1_t *b, char *ref, int thres) int bam_prob_realn_core(bam1_t *b, const char *ref, int flag) { - int k, i, bw, x, y, yb, ye, xb, xe, apply_baq = flag&1, extend_baq = flag>>1&1; + int k, i, bw, x, y, yb, ye, xb, xe, apply_baq = flag&1, extend_baq = flag>>1&1, redo_baq = flag&4; uint32_t *cigar = bam1_cigar(b); bam1_core_t *c = &b->core; kpa_par_t conf = kpa_par_def; @@ -197,6 +197,11 @@ int bam_prob_realn_core(bam1_t *b, const char *ref, int flag) // test if BQ or ZQ is present if ((bq = bam_aux_get(b, "BQ")) != 0) ++bq; if ((zq = bam_aux_get(b, "ZQ")) != 0 && *zq == 'Z') ++zq; + if (bq && redo_baq) + { + bam_aux_del(b, bq-1); + bq = 0; + } if (bq && zq) { // remove the ZQ tag bam_aux_del(b, zq-1); zq = 0; diff --git a/bam_plcmd.c b/bam_plcmd.c index be2c87d..07f0a4f 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -67,7 +67,7 @@ static inline void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, cons #define MPLP_NO_ORPHAN 0x40 #define MPLP_REALN 0x80 #define MPLP_NO_INDEL 0x400 -#define MPLP_EXT_BAQ 0x800 +#define MPLP_REDO_BAQ 0x800 #define MPLP_ILLUMINA13 0x1000 #define MPLP_IGNORE_RG 0x2000 #define MPLP_PRINT_POS 0x4000 @@ -133,7 +133,7 @@ static int mplp_func(void *data, bam1_t *b) } has_ref = (ma->ref && ma->ref_id == b->core.tid)? 1 : 0; skip = 0; - if (has_ref && (ma->conf->flag&MPLP_REALN)) bam_prob_realn_core(b, ma->ref, (ma->conf->flag & MPLP_EXT_BAQ)? 3 : 1); + if (has_ref && (ma->conf->flag&MPLP_REALN)) bam_prob_realn_core(b, ma->ref, (ma->conf->flag & MPLP_REDO_BAQ)? 7 : 3); if (has_ref && ma->conf->capQ_thres > 10) { int q = bam_cap_mapQ(b, ma->ref, ma->conf->capQ_thres); if (q < 0) skip = 1; @@ -448,7 +448,7 @@ int bam_mpileup(int argc, char *argv[]) mplp.max_depth = 250; mplp.max_indel_depth = 250; mplp.openQ = 40; mplp.extQ = 20; mplp.tandemQ = 100; mplp.min_frac = 0.002; mplp.min_support = 1; - mplp.flag = MPLP_NO_ORPHAN | MPLP_REALN | MPLP_EXT_BAQ; + mplp.flag = MPLP_NO_ORPHAN | MPLP_REALN; while ((c = getopt(argc, argv, "Agf:r:l:M:q:Q:uaRC:BDSd:L:b:P:o:e:h:Im:F:EG:6OsV")) >= 0) { switch (c) { case 'f': @@ -467,7 +467,7 @@ int bam_mpileup(int argc, char *argv[]) case 'S': mplp.fmt_flag |= B2B_FMT_SP; break; case 'V': mplp.fmt_flag |= B2B_FMT_DV; break; case 'I': mplp.flag |= MPLP_NO_INDEL; break; - case 'E': mplp.flag |= MPLP_EXT_BAQ; break; + case 'E': mplp.flag |= MPLP_REDO_BAQ; break; case '6': mplp.flag |= MPLP_ILLUMINA13; break; case 'R': mplp.flag |= MPLP_IGNORE_RG; break; case 's': mplp.flag |= MPLP_PRINT_MAPQ; break; @@ -508,7 +508,7 @@ int bam_mpileup(int argc, char *argv[]) fprintf(stderr, " -b FILE list of input BAM files [null]\n"); fprintf(stderr, " -C INT parameter for adjusting mapQ; 0 to disable [0]\n"); fprintf(stderr, " -d INT max per-BAM depth to avoid excessive memory usage [%d]\n", mplp.max_depth); -// fprintf(stderr, " -E extended BAQ for higher sensitivity but lower specificity\n"); + fprintf(stderr, " -E recalculate extended BAQ on the fly thus ignoring existing BQs\n"); fprintf(stderr, " -f FILE faidx indexed reference sequence file [null]\n"); fprintf(stderr, " -G FILE exclude read groups listed in FILE [null]\n"); fprintf(stderr, " -l FILE list of positions (chr pos) or regions (BED) [null]\n"); -- 2.39.2