From 7d63e76245ba668cf625a327972eb33fc5e95140 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sun, 16 Mar 2014 23:25:04 -0400 Subject: [PATCH] r444: more debugging output in CIGAR generation Also found a potential issue which should not affect accuracy but may hurt speed. Will investigate later. --- bwa.c | 8 ++++++-- bwamem.c | 12 ++++++++---- main.c | 2 +- 3 files changed, 15 insertions(+), 7 deletions(-) diff --git a/bwa.c b/bwa.c index 140d57e..0a02905 100644 --- a/bwa.c +++ b/bwa.c @@ -106,6 +106,7 @@ uint32_t *bwa_gen_cigar(const int8_t mat[25], int q, int r, int w_, int64_t l_pa tmp = rseq[i], rseq[i] = rseq[rlen - 1 - i], rseq[rlen - 1 - i] = tmp; } if (l_query == re - rb && w_ == 0) { // no gap; no need to do DP + // FIXME: due to an issue in mem_reg2aln(), we never come to this block. This does not affect accuracy, but it hurts performance. cigar = malloc(4); cigar[0] = l_query<<4 | 0; *n_cigar = 1; @@ -113,8 +114,6 @@ uint32_t *bwa_gen_cigar(const int8_t mat[25], int q, int r, int w_, int64_t l_pa *score += mat[rseq[i]*5 + query[i]]; } else { int w, max_gap, min_w; - //printf("[Q] "); for (i = 0; i < l_query; ++i) putchar("ACGTN"[(int)query[i]]); putchar('\n'); - //printf("[R] "); for (i = 0; i < re - rb; ++i) putchar("ACGTN"[(int)rseq[i]]); putchar('\n'); // set the band-width max_gap = (int)((double)(((l_query+1)>>1) * mat[0] - q) / r + 1.); max_gap = max_gap > 1? max_gap : 1; @@ -123,6 +122,11 @@ uint32_t *bwa_gen_cigar(const int8_t mat[25], int q, int r, int w_, int64_t l_pa min_w = abs(rlen - l_query) + 3; w = w > min_w? w : min_w; // NW alignment + if (bwa_verbose >= 4) { + printf("* Global bandwidth: %d\n", w); + printf("* Global ref: "); for (i = 0; i < rlen; ++i) putchar("ACGTN"[(int)rseq[i]]); putchar('\n'); + printf("* Global query: "); for (i = 0; i < l_query; ++i) putchar("ACGTN"[(int)query[i]]); putchar('\n'); + } *score = ksw_global(l_query, query, rlen, rseq, 5, mat, q, r, w, n_cigar, &cigar); } {// compute NM and MD diff --git a/bwamem.c b/bwamem.c index c9b32b0..0ca116f 100644 --- a/bwamem.c +++ b/bwamem.c @@ -948,14 +948,14 @@ mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t * exit(1); } w2 = infer_bw(qe - qb, re - rb, ar->truesc, opt->a, opt->q, opt->r); - if (bwa_verbose >= 4) printf("Band width: infer=%d, opt=%d, alnreg=%d\n", w2, opt->w, ar->w); + if (bwa_verbose >= 4) printf("* Band width: inferred=%d, cmd_opt=%d, alnreg=%d\n", w2, opt->w, ar->w); if (w2 > opt->w) w2 = w2 < ar->w? w2 : ar->w; - else w2 = opt->w; + else w2 = opt->w; // TODO: check if we need this line. Need to test on many reads. i = 0; a.cigar = 0; do { free(a.cigar); a.cigar = bwa_gen_cigar(opt->mat, opt->q, opt->r, w2, bns->l_pac, pac, qe - qb, (uint8_t*)&query[qb], rb, re, &score, &a.n_cigar, &NM); - if (bwa_verbose >= 4) printf("Final alignment: w2=%d, global_sc=%d, local_sc=%d\n", w2, score, ar->truesc); + if (bwa_verbose >= 4) printf("* Final alignment: w2=%d, global_sc=%d, local_sc=%d\n", w2, score, ar->truesc); if (score == last_sc) break; // it is possible that global alignment and local alignment give different scores last_sc = score; w2 <<= 1; @@ -1010,11 +1010,13 @@ 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)) { + if (bwa_verbose >= 4) printf("=====> Processing read '%s' <=====\n", w->seqs[i].name); w->regs[i] = mem_align1_core(w->opt, w->bwt, w->bns, w->pac, w->seqs[i].l_seq, w->seqs[i].seq); } else { + if (bwa_verbose >= 4) printf("=====> Processing read '%s'/1 <=====\n", w->seqs[i<<1|0].name); w->regs[i<<1|0] = mem_align1_core(w->opt, w->bwt, w->bns, w->pac, w->seqs[i<<1|0].l_seq, w->seqs[i<<1|0].seq); + if (bwa_verbose >= 4) printf("=====> Processing read '%s'/2 <=====\n", w->seqs[i<<1|1].name); w->regs[i<<1|1] = mem_align1_core(w->opt, w->bwt, w->bns, w->pac, w->seqs[i<<1|1].l_seq, w->seqs[i<<1|1].seq); } } @@ -1024,10 +1026,12 @@ static void worker2(void *data, int i, int tid) extern int mem_sam_pe(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, const mem_pestat_t pes[4], uint64_t id, bseq1_t s[2], mem_alnreg_v a[2]); worker_t *w = (worker_t*)data; if (!(w->opt->flag&MEM_F_PE)) { + if (bwa_verbose >= 4) printf("=====> Finalizing read '%s' <=====\n", w->seqs[i].name); mem_mark_primary_se(w->opt, w->regs[i].n, w->regs[i].a, w->n_processed + i); mem_reg2sam_se(w->opt, w->bns, w->pac, &w->seqs[i], &w->regs[i], 0, 0); free(w->regs[i].a); } else { + if (bwa_verbose >= 4) printf("=====> Finalizing read pair '%s' <=====\n", w->seqs[i<<1|0].name); mem_sam_pe(w->opt, w->bns, w->pac, w->pes, (w->n_processed>>1) + i, &w->seqs[i<<1], &w->regs[i<<1]); free(w->regs[i<<1|0].a); free(w->regs[i<<1|1].a); } diff --git a/main.c b/main.c index 2fd0138..0c355fc 100644 --- a/main.c +++ b/main.c @@ -4,7 +4,7 @@ #include "utils.h" #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.7.7-r443" +#define PACKAGE_VERSION "0.7.7-r444" #endif int bwa_fa2pac(int argc, char *argv[]);