From 9346acde1b2c930e13948b7cb43076f2ebc63dfc Mon Sep 17 00:00:00 2001 From: Heng Li Date: Fri, 15 Mar 2013 21:26:37 -0400 Subject: [PATCH] Release bwa-0.7.3a-r367 In 0.7.3, the wrong CIGAR bug was only fixed in one scenario, but not fixed in another corner case. --- NEWS | 10 ++++++++++ bwa.1 | 2 +- bwamem.c | 24 +++++++++++++++++------- bwamem.h | 3 ++- main.c | 2 +- 5 files changed, 31 insertions(+), 10 deletions(-) diff --git a/NEWS b/NEWS index 788ae404..0cf0591b 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,13 @@ +Release 0.7.3a (15 March, 2013) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +In 0.7.3, the wrong CIGAR bug was only fixed in one scenario, but not fixed +in another corner case. + +(0.7.3a: 15 March 2013, r367) + + + Release 0.7.3 (15 March, 2013) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/bwa.1 b/bwa.1 index a02aebe5..1edbf12a 100644 --- a/bwa.1 +++ b/bwa.1 @@ -1,4 +1,4 @@ -.TH bwa 1 "15 March 2013" "bwa-0.7.3" "Bioinformatics tools" +.TH bwa 1 "15 March 2013" "bwa-0.7.3a" "Bioinformatics tools" .SH NAME .PP bwa - Burrows-Wheeler Alignment Tool diff --git a/bwamem.c b/bwamem.c index da9f2bed..82dbe454 100644 --- a/bwamem.c +++ b/bwamem.c @@ -542,7 +542,7 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int a = kv_pushp(mem_alnreg_t, *av); memset(a, 0, sizeof(mem_alnreg_t)); a->w = aw[0] = aw[1] = opt->w; - a->score = -1; + a->score = a->truesc = -1; if (s->qbeg) { // left extension uint8_t *rs, *qs; @@ -560,10 +560,15 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int 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 - if (gscore <= 0 || gscore <= a->score - opt->pen_clip) a->qb = s->qbeg - qle, a->rb = s->rbeg - tle; // local hits - else a->qb = 0, a->rb = s->rbeg - gtle; // reach the end + if (gscore <= 0 || gscore <= a->score - opt->pen_clip) { // local extension + a->qb = s->qbeg - qle, a->rb = s->rbeg - tle; + a->truesc = a->score; + } else { // to-end extension + a->qb = 0, a->rb = s->rbeg - gtle; + a->truesc = gscore; + } free(qs); free(rs); - } else a->score = s->len * opt->a, a->qb = 0, a->rb = s->rbeg; + } else a->score = a->truesc = s->len * opt->a, a->qb = 0, a->rb = s->rbeg; if (s->qbeg + s->len != l_query) { // right extension int qle, tle, qe, re, gtle, gscore, sc0 = a->score; @@ -578,8 +583,13 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int if (a->score == prev || max_off[1] < (aw[1]>>1) + (aw[1]>>2)) break; } // similar to the above - if (gscore <= 0 || gscore <= a->score - opt->pen_clip) a->qe = qe + qle, a->re = rmax[0] + re + tle; - else a->qe = l_query, a->re = rmax[0] + re + gtle; + if (gscore <= 0 || gscore <= a->score - opt->pen_clip) { // local extension + a->qe = qe + qle, a->re = rmax[0] + re + tle; + a->truesc += a->score - sc0; + } else { // to-end extension + a->qe = l_query, a->re = rmax[0] + re + gtle; + 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); } @@ -839,7 +849,7 @@ mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t * a.mapq = ar->secondary < 0? mem_approx_mapq_se(opt, ar) : 0; if (ar->secondary >= 0) a.flag |= 0x20000; bwa_fix_xref(opt->mat, opt->q, opt->r, opt->w, bns, pac, (uint8_t*)query, &qb, &qe, &rb, &re); - w2 = infer_bw(qe - qb, re - rb, ar->score, opt->a, opt->q, opt->r); + w2 = infer_bw(qe - qb, re - rb, ar->truesc, opt->a, opt->q, opt->r); w2 = w2 < opt->w? w2 : opt->w; 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); a.NM = NM; diff --git a/bwamem.h b/bwamem.h index ae201e6a..b7a4e795 100644 --- a/bwamem.h +++ b/bwamem.h @@ -43,7 +43,8 @@ typedef struct { typedef struct { int64_t rb, re; // [rb,re): reference sequence in the alignment int qb, qe; // [qb,qe): query sequence in the alignment - int score; // best SW score + int score; // best local SW score + int truesc; // actual score corresponding to the aligned region; possibly smaller than $score int sub; // 2nd best SW score int csub; // SW score of a tandem hit int sub_n; // approximate number of suboptimal hits diff --git a/main.c b/main.c index a03c49f0..d9e75e2a 100644 --- a/main.c +++ b/main.c @@ -3,7 +3,7 @@ #include "utils.h" #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.7.3-r366" +#define PACKAGE_VERSION "0.7.3a-r367" #endif int bwa_fa2pac(int argc, char *argv[]);