From ea3dc2f00300131baa7303c92326a7ab099d0288 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 29 Jan 2014 10:51:02 -0500 Subject: [PATCH] r430: fix a bug producing incorrect alignment Ksw uses two rounds of SSE2-SW to find the boundaries of an alignment. If the second round gives a different score from the first round, it will fail. The fix checks if this happens, though I have not dig into an example to understand why this may happen in the first place. --- bwamem.c | 4 +++- bwamem_pair.c | 2 +- main.c | 2 +- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/bwamem.c b/bwamem.c index 9b05528..27b9c4c 100644 --- a/bwamem.c +++ b/bwamem.c @@ -910,7 +910,7 @@ mem_alnreg_v mem_align1(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t * mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, int l_query, const char *query_, const mem_alnreg_t *ar) { mem_aln_t a; - int i, w2, qb, qe, NM, score, is_rev; + int i, w2, qb, qe, NM, score, is_rev, last_sc = -(1<<30); int64_t pos, rb, re; uint8_t *query; @@ -939,6 +939,8 @@ mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t * 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 (score == last_sc) break; // it is possible that global alignment and local alignment give different scores + last_sc = score; w2 <<= 1; } while (++i < 3 && score < ar->truesc - opt->a); a.NM = NM; diff --git a/bwamem_pair.c b/bwamem_pair.c index 1729b8c..f1aa73a 100644 --- a/bwamem_pair.c +++ b/bwamem_pair.c @@ -147,7 +147,7 @@ int mem_matesw(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, const me int tmp, xtra = KSW_XSUBO | KSW_XSTART | (l_ms * opt->a < 250? KSW_XBYTE : 0) | opt->min_seed_len; aln = ksw_align(l_ms, seq, len, ref, 5, opt->mat, opt->q, opt->r, xtra, 0); memset(&b, 0, sizeof(mem_alnreg_t)); - if (aln.score >= opt->min_seed_len) { + if (aln.score >= opt->min_seed_len && aln.qb >= 0) { // something goes wrong if aln.qb < 0 b.qb = is_rev? l_ms - (aln.qe + 1) : aln.qb; b.qe = is_rev? l_ms - aln.qb : aln.qe + 1; b.rb = is_rev? (l_pac<<1) - (rb + aln.te + 1) : rb + aln.tb; diff --git a/main.c b/main.c index 69bf765..20cbf83 100644 --- a/main.c +++ b/main.c @@ -4,7 +4,7 @@ #include "utils.h" #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.7.5a-r428" +#define PACKAGE_VERSION "0.7.5a-r430" #endif int bwa_fa2pac(int argc, char *argv[]);