From 8929bd1c25af96b0b4456a33652bd81756b6623c Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sun, 16 Mar 2014 15:18:58 -0400 Subject: [PATCH] r443: more verbose debugging information --- bwamem.c | 35 +++++++++++++++++++++++++---------- main.c | 2 +- 2 files changed, 26 insertions(+), 11 deletions(-) diff --git a/bwamem.c b/bwamem.c index 19ca561..c9b32b0 100644 --- a/bwamem.c +++ b/bwamem.c @@ -215,7 +215,7 @@ static void mem_insert_seed(const mem_opt_t *opt, int64_t l_pac, kbtree_t(chn) * s.rbeg = tmp.pos = bwt_sa(itr->bwt, p->x[0] + k); // this is the base coordinate in the forward-reverse reference s.qbeg = p->info>>32; s.len = slen; - if (bwa_verbose >= 5) printf("SEED l=%d,qb=%d,rb=%ld\n", s.len, s.qbeg, (long)s.rbeg); + if (bwa_verbose >= 5) printf("* Found SEED: length=%d,query_beg=%d,ref_beg=%ld\n", s.len, s.qbeg, (long)s.rbeg); if (s.rbeg < l_pac && l_pac < s.rbeg + s.len) continue; // bridging forward-reverse boundary; skip if (kb_size(tree)) { kb_intervalp(chn, tree, &tmp, &lower, &upper); // find the closest chain @@ -257,14 +257,14 @@ void mem_print_chain(const bntseq_t *bns, mem_chain_v *chn) int i, j; for (i = 0; i < chn->n; ++i) { mem_chain_t *p = &chn->a[i]; - err_printf("CHAIN(%d) n=%d w=%d", i, p->n, mem_chain_weight(p)); + err_printf("* Found CHAIN(%d): n=%d; weight=%d", i, p->n, mem_chain_weight(p)); for (j = 0; j < p->n; ++j) { bwtint_t pos; int is_rev, ref_id; pos = bns_depos(bns, p->seeds[j].rbeg, &is_rev); if (is_rev) pos -= p->seeds[j].len - 1; bns_cnt_ambi(bns, pos, p->seeds[j].len, &ref_id); - err_printf("\t%d,%d,%ld(%s:%c%ld)", p->seeds[j].len, p->seeds[j].qbeg, (long)p->seeds[j].rbeg, bns->anns[ref_id].name, "+-"[is_rev], (long)(pos - bns->anns[ref_id].offset) + 1); + err_printf("\t%d;%d,%ld(%s:%c%ld)", p->seeds[j].len, p->seeds[j].qbeg, (long)p->seeds[j].rbeg, bns->anns[ref_id].name, "+-"[is_rev], (long)(pos - bns->anns[ref_id].offset) + 1); } err_putchar('\n'); } @@ -518,7 +518,7 @@ int mem_chain2aln_short(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, a.score = x.score; a.csub = x.score2; kv_push(mem_alnreg_t, *av, a); - if (bwa_verbose >= 4) printf("chain2aln(short): [%d,%d) <=> [%ld,%ld)\n", a.qb, a.qe, (long)a.rb, (long)a.re); + if (bwa_verbose >= 4) printf("** Added alignment region via mem_chain2aln_short(): [%d,%d) <=> [%ld,%ld)\n", a.qb, a.qe, (long)a.rb, (long)a.re); return 0; } @@ -584,7 +584,9 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int w = max_gap < opt->w? max_gap : opt->w; if (qd - rd < w && rd - qd < w) break; } - if (i < av->n) { // the seed is (almost) contained in an existing alignment + if (i < av->n) { // the seed is (almost) contained in an existing alignment; further testing is needed to confirm it is not leading to a different aln + if (bwa_verbose >= 4) + printf("** Seed(%d) [%ld;%ld,%ld] is almost contained in an existing alignment. Confirming whether extension is needed...\n", k, (long)s->len, (long)s->qbeg, (long)s->rbeg); for (i = k + 1; i < c->n; ++i) { // check overlapping seeds in the same chain const mem_seed_t *t; if (srt[i] == 0) continue; @@ -597,6 +599,8 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int srt[k] = 0; // mark that seed extension has not been performed continue; } + if (bwa_verbose >= 4) + printf("** Seed(%d) might lead to a different alignment even though it is contained. Extension will be performed.\n", k); } a = kv_pushp(mem_alnreg_t, *av); @@ -604,7 +608,7 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int a->w = aw[0] = aw[1] = opt->w; a->score = a->truesc = -1; - if (bwa_verbose >= 4) err_printf("Extending from seed [%ld,%ld,%ld]\n", (long)s->len, (long)s->qbeg, (long)s->rbeg); + if (bwa_verbose >= 4) err_printf("** ---> Extending from seed(%d) [%ld;%ld,%ld] <---\n", k, (long)s->len, (long)s->qbeg, (long)s->rbeg); if (s->qbeg) { // left extension uint8_t *rs, *qs; int qle, tle, gtle, gscore; @@ -616,8 +620,13 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int for (i = 0; i < MAX_BAND_TRY; ++i) { int prev = a->score; aw[0] = opt->w << i; + if (bwa_verbose >= 4) { + int j; + printf("*** Left ref: "); for (j = 0; j < tmp; ++j) putchar("ACGTN"[(int)rs[j]]); putchar('\n'); + printf("*** Left query: "); for (j = 0; j < s->qbeg; ++j) putchar("ACGTN"[(int)qs[j]]); putchar('\n'); + } a->score = ksw_extend(s->qbeg, qs, tmp, rs, 5, opt->mat, opt->q, opt->r, aw[0], opt->pen_clip5, opt->zdrop, s->len * opt->a, &qle, &tle, >le, &gscore, &max_off[0]); - if (bwa_verbose >= 4) { printf("L\t%d < %d; w=%d; max_off=%d\n", prev, a->score, aw[0], max_off[0]); fflush(stdout); } + if (bwa_verbose >= 4) { printf("*** Left extension: prev_score=%d; score=%d; bandwidth=%d; max_off_diagonal_dist=%d\n", prev, a->score, aw[0], max_off[0]); fflush(stdout); } if (a->score == prev || max_off[0] < (aw[0]>>1) + (aw[0]>>2)) break; } // check whether we prefer to reach the end of the query @@ -639,8 +648,13 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int for (i = 0; i < MAX_BAND_TRY; ++i) { int prev = a->score; aw[1] = opt->w << i; + if (bwa_verbose >= 4) { + int j; + printf("*** Right ref: "); for (j = 0; j < rmax[1] - rmax[0] - re; ++j) putchar("ACGTN"[(int)rseq[re+j]]); putchar('\n'); + printf("*** Right query: "); for (j = 0; j < l_query - qe; ++j) putchar("ACGTN"[(int)query[qe+j]]); putchar('\n'); + } a->score = ksw_extend(l_query - qe, query + qe, rmax[1] - rmax[0] - re, rseq + re, 5, opt->mat, opt->q, opt->r, aw[1], opt->pen_clip3, opt->zdrop, sc0, &qle, &tle, >le, &gscore, &max_off[1]); - if (bwa_verbose >= 4) { printf("R\t%d < %d; w=%d; max_off=%d\n", prev, a->score, aw[1], max_off[1]); fflush(stdout); } + if (bwa_verbose >= 4) { printf("*** Right extension: prev_score=%d; score=%d; bandwidth=%d; max_off_diagonal_dist=%d\n", prev, a->score, aw[1], max_off[1]); fflush(stdout); } if (a->score == prev || max_off[1] < (aw[1]>>1) + (aw[1]>>2)) break; } // similar to the above @@ -652,7 +666,7 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int a->truesc += gscore - sc0; } } else a->qe = l_query, a->re = s->rbeg + s->len; - if (bwa_verbose >= 4) { printf("[%d]\taw={%d,%d}\tscore=%d\t[%d,%d) <=> [%ld,%ld)\n", k, aw[0], aw[1], a->score, a->qb, a->qe, (long)a->rb, (long)a->re); fflush(stdout); } + if (bwa_verbose >= 4) printf("*** Added alignment region: [%d,%d) <=> [%ld,%ld); score=%d; {left,right}_bandwidth={%d,%d}\n", a->qb, a->qe, (long)a->rb, (long)a->re, a->score, aw[0], aw[1]); // compute seedcov for (i = 0, a->seedcov = 0; i < c->n; ++i) { @@ -887,7 +901,7 @@ mem_alnreg_v mem_align1_core(const mem_opt_t *opt, const bwt_t *bwt, const bntse for (i = 0; i < chn.n; ++i) { mem_chain_t *p = &chn.a[i]; int ret; - if (bwa_verbose >= 4) err_printf("===> Processing chain(%d) <===\n", i); + if (bwa_verbose >= 4) err_printf("* ---> Processing chain(%d) <---\n", i); ret = mem_chain2aln_short(opt, bns->l_pac, pac, l_seq, (uint8_t*)seq, p, ®s); if (ret > 0) mem_chain2aln(opt, bns->l_pac, pac, l_seq, (uint8_t*)seq, p, ®s); free(chn.a[i].seeds); @@ -996,6 +1010,7 @@ typedef struct { static void worker1(void *data, int i, int tid) { worker_t *w = (worker_t*)data; + if (bwa_verbose >= 4) printf("=====> Processing read '%s' <=====\n", w->seqs[i].name); if (!(w->opt->flag&MEM_F_PE)) { w->regs[i] = mem_align1_core(w->opt, w->bwt, w->bns, w->pac, w->seqs[i].l_seq, w->seqs[i].seq); } else { diff --git a/main.c b/main.c index a8df9c0..2fd0138 100644 --- a/main.c +++ b/main.c @@ -4,7 +4,7 @@ #include "utils.h" #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.7.7-r441" +#define PACKAGE_VERSION "0.7.7-r443" #endif int bwa_fa2pac(int argc, char *argv[]);