From e6c262594fcb6c48aa4f1ced3f49827cfdb5543f Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@me.com>
Date: Tue, 5 Mar 2013 10:12:38 -0500
Subject: [PATCH] bwa-sw: ditch stdaln

---
 bwtsw2_aux.c | 38 +++++++++++++++++---------------------
 ksw.c        |  2 +-
 2 files changed, 18 insertions(+), 22 deletions(-)

diff --git a/bwtsw2_aux.c b/bwtsw2_aux.c
index 4455c5f..6527495 100644
--- a/bwtsw2_aux.c
+++ b/bwtsw2_aux.c
@@ -11,9 +11,9 @@
 #include "bwt_lite.h"
 #include "utils.h"
 #include "bwtsw2.h"
-#include "stdaln.h"
 #include "kstring.h"
 #include "bwa.h"
+#include "ksw.h"
 
 #include "kseq.h"
 KSEQ_DECLARE(gzFile)
@@ -94,13 +94,12 @@ bwtsw2_t *bsw2_dup_no_cigar(const bwtsw2_t *b)
 
 void bsw2_extend_left(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *_query, int lq, uint8_t *pac, bwtint_t l_pac, uint8_t *_mem)
 {
-	int i, matrix[25];
+	int i;
 	bwtint_t k;
 	uint8_t *target = 0, *query;
-	AlnParam par;
+	int8_t mat[25];
 
-	par.matrix = matrix;
-	__gen_ap(par, opt);
+	bwa_fill_scmat(opt->a, opt->b, mat);
 	query = calloc(lq, 1);
 	// sort according to the descending order of query end
 	ks_introsort(hit, b->n, b->hits);
@@ -111,8 +110,7 @@ void bsw2_extend_left(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *_query, int lq
 	for (i = 0; i < b->n; ++i) {
 		bsw2hit_t *p = b->hits + i;
 		int lt = ((p->beg + 1) / 2 * opt->a + opt->r) / opt->r + lq;
-		int score, j;
-		path_t path;
+		int score, j, qle, tle;
 		p->n_seeds = 1;
 		if (p->l || p->k == 0) continue;
 		for (j = score = 0; j < i; ++j) {
@@ -127,12 +125,12 @@ void bsw2_extend_left(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *_query, int lq
 		for (k = p->k - 1, j = 0; k > 0 && j < lt; --k) // FIXME: k=0 not considered!
 			target[j++] = pac[k>>2] >> (~k&3)*2 & 0x3;
 		lt = j;
-		score = aln_extend_core(target, lt, query + lq - p->beg, p->beg, &par, &path, 0, p->G, _mem);
+		score = ksw_extend(p->beg, &query[lq - p->beg], lt, target, 5, mat, opt->q, opt->r, opt->bw, -1, p->G, &qle, &tle, 0, 0, 0);
 		if (score > p->G) { // extensible
 			p->G = score;
-			p->len += path.i;
-			p->beg -= path.j;
-			p->k -= path.i;
+			p->k -= tle;
+			p->len += tle;
+			p->beg -= qle;
 		}
 	}
 	free(query); free(target);
@@ -140,29 +138,27 @@ void bsw2_extend_left(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *_query, int lq
 
 void bsw2_extend_rght(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *query, int lq, uint8_t *pac, bwtint_t l_pac, uint8_t *_mem)
 {
-	int i, matrix[25];
+	int i;
 	bwtint_t k;
 	uint8_t *target;
-	AlnParam par;
-	
-	par.matrix = matrix;
-	__gen_ap(par, opt);
+	int8_t mat[25];
+
+	bwa_fill_scmat(opt->a, opt->b, mat);
 	target = calloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1);
 	for (i = 0; i < b->n; ++i) {
 		bsw2hit_t *p = b->hits + i;
 		int lt = ((lq - p->beg + 1) / 2 * opt->a + opt->r) / opt->r + lq;
-		int j, score;
-		path_t path;
+		int j, score, qle, tle;
 		if (p->l) continue;
 		for (k = p->k, j = 0; k < p->k + lt && k < l_pac; ++k)
 			target[j++] = pac[k>>2] >> (~k&3)*2 & 0x3;
 		lt = j;
-		score = aln_extend_core(target, lt, query + p->beg, lq - p->beg, &par, &path, 0, 1, _mem);
+		score = ksw_extend(lq - p->beg, &query[p->beg], lt, target, 5, mat, opt->q, opt->r, opt->bw, -1, 1, &qle, &tle, 0, 0, 0);
 //		if (score < p->G) fprintf(stderr, "[bsw2_extend_hits] %d < %d\n", score, p->G);
 		if (score >= p->G) {
 			p->G = score;
-			p->len = path.i;
-			p->end = path.j + p->beg;
+			p->len = tle;
+			p->end = p->beg + qle;
 		}
 	}
 	free(target);
diff --git a/ksw.c b/ksw.c
index 5666a8f..e331390 100644
--- a/ksw.c
+++ b/ksw.c
@@ -426,7 +426,7 @@ int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target,
 			max_ie = gscore > h1? max_ie : i;
 			gscore = gscore > h1? gscore : h1;
 		}
-		if (m == 0 || max - m - abs((i - max_i) - (j - max_j)) * gape > zdrop) break; // drop to zero, or below Z-dropoff
+		if (m == 0 || (zdrop > 0 && max - m - abs((i - max_i) - (j - max_j)) * gape > zdrop)) break; // drop to zero, or below Z-dropoff
 		if (m > max) {
 			max = m, max_i = i, max_j = mj;
 			max_off = max_off > abs(mj - i)? max_off : abs(mj - i);