From 1e3cadbfc29714af36cb9fe9475ebaad9eb5a546 Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@me.com>
Date: Mon, 18 Mar 2013 20:49:32 -0400
Subject: [PATCH] r368: bugfix - wrong CIGAR when bridging 3 contigs

In this case, bwa_fix_xref() will return insane coordinates. The old version
did not check the return status and write wrong CIGAR. This bug only happen to
very short assembly contigs.
---
 bwa.c    | 2 +-
 bwamem.c | 9 ++++++++-
 main.c   | 2 +-
 3 files changed, 10 insertions(+), 3 deletions(-)

diff --git a/bwa.c b/bwa.c
index 08d96b8..e86a87e 100644
--- a/bwa.c
+++ b/bwa.c
@@ -144,7 +144,7 @@ int bwa_fix_xref(const int8_t mat[25], int q, int r, int w, const bntseq_t *bns,
 {
 	int ib, ie, is_rev;
 	int64_t fb, fe, mid = -1;
-	if (*rb < bns->l_pac && *re > bns->l_pac) { // cross the for-rev boundary
+	if (*rb < bns->l_pac && *re > bns->l_pac) { // cross the for-rev boundary; actually with BWA-MEM, we should never come to here
 		*qb = *qe = *rb = *re = -1;
 		return -1; // unable to fix
 	} else {
diff --git a/bwamem.c b/bwamem.c
index 82dbe45..1af7f87 100644
--- a/bwamem.c
+++ b/bwamem.c
@@ -772,6 +772,10 @@ void mem_reg2sam_se(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pa
 			if (p->secondary >= 0 && p->score < a->a[p->secondary].score * .5) continue;
 			q = kv_pushp(mem_aln_t, aa);
 			*q = mem_reg2aln(opt, bns, pac, s->l_seq, s->seq, p);
+			if (q->rid < 0) { // unfixable cross-reference alignment
+				--aa.n;
+				continue;
+			}
 			q->flag |= extra_flag | (p->secondary >= 0? 0x100 : 0); // flag secondary
 			if (p->secondary >= 0) q->sub = -1; // don't output sub-optimal score
 			if ((opt->flag&MEM_F_NO_MULTI) && k && p->secondary < 0) q->flag |= 0x10000;
@@ -848,7 +852,10 @@ mem_aln_t mem_reg2aln(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *
 		query[i] = query_[i] < 5? query_[i] : nst_nt4_table[(int)query_[i]];
 	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);
+	if (bwa_fix_xref(opt->mat, opt->q, opt->r, opt->w, bns, pac, (uint8_t*)query, &qb, &qe, &rb, &re) < 0) { // unfixable cross-reference alignment
+		a.rid = -1; a.pos = -1; a.flag |= 0x4;
+		return a;
+	}
 	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);
diff --git a/main.c b/main.c
index d9e75e2..0a79581 100644
--- a/main.c
+++ b/main.c
@@ -3,7 +3,7 @@
 #include "utils.h"
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.7.3a-r367"
+#define PACKAGE_VERSION "0.7.3-r368-beta"
 #endif
 
 int bwa_fa2pac(int argc, char *argv[]);