From 72dfb0c99e2f66c8e16c03ab935ca67592f2fa03 Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@me.com>
Date: Sun, 25 Jun 2017 10:22:13 -0400
Subject: [PATCH] fixed a bug in ksw2

---
 align.c          | 2 ++
 ksw2_extz2_sse.c | 2 +-
 main.c           | 4 +++-
 3 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/align.c b/align.c
index 4d55871..7a0f298 100644
--- a/align.c
+++ b/align.c
@@ -180,6 +180,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			qseq = &qseq0[rev][qs];
 			mm_idx_getseq(mi, rid, rs, re, tseq);
 			#if 0
+			int k;
 			fprintf(stderr, "===> [%d] %d-%d %c (%s:%d-%d) <===\n", i, qs, qe, "+-"[rev], mi->seq[rid].name, rs, re);
 			for (k = 0; k < re - rs; ++k) fputc("ACGTN"[tseq[k]], stderr); fputc('\n', stderr);
 			for (k = 0; k < qe - qs; ++k) fputc("ACGTN"[qseq[k]], stderr); fputc('\n', stderr);
@@ -189,6 +190,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			mm_update_extra(r->p, qseq, tseq, ez->n_cigar, ez->cigar);
 			if (ez->score == KSW_NEG_INF) { // truncated by Z-drop
 				int j;
+				abort();
 				for (j = i - 1; j >= 0; --j)
 					if ((int32_t)a[r->as + j].x < re + ez->max_t)
 						break;
diff --git a/ksw2_extz2_sse.c b/ksw2_extz2_sse.c
index f291331..1bdf962 100644
--- a/ksw2_extz2_sse.c
+++ b/ksw2_extz2_sse.c
@@ -53,7 +53,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	m1_     = _mm_set1_epi8(m - 1); // wildcard
 
 	ez->max_q = ez->max_t = ez->mqe_t = ez->mte_q = -1;
-	ez->max = 0, ez->mqe = ez->mte = KSW_NEG_INF;
+	ez->max = 0, ez->score = ez->mqe = ez->mte = KSW_NEG_INF;
 	ez->n_cigar = 0;
 
 	wl = wr = w;
diff --git a/main.c b/main.c
index 4b9f766..9ef7119 100644
--- a/main.c
+++ b/main.c
@@ -54,7 +54,7 @@ int main(int argc, char *argv[])
 	mm_realtime0 = realtime();
 	mm_mapopt_init(&opt);
 
-	while ((c = getopt(argc, argv, "w:k:B:b:t:r:f:Vv:Ng:I:d:ST:s:Dx:Hp:m:")) >= 0) {
+	while ((c = getopt(argc, argv, "w:k:B:b:t:r:f:Vv:Ng:I:d:ST:s:Dx:Hp:m:z:")) >= 0) {
 		if (c == 'w') w = atoi(optarg);
 		else if (c == 'k') k = atoi(optarg);
 		else if (c == 'b') b = atoi(optarg);
@@ -72,6 +72,7 @@ int main(int argc, char *argv[])
 		else if (c == 'S') opt.flag |= MM_F_AVA | MM_F_NO_SELF;
 		else if (c == 'T') opt.sdust_thres = atoi(optarg);
 		else if (c == 's') opt.min_score = atoi(optarg);
+		else if (c == 'z') opt.zdrop = atoi(optarg);
 		else if (c == 'V') {
 			puts(MM_VERSION);
 			return 0;
@@ -113,6 +114,7 @@ int main(int argc, char *argv[])
 //		fprintf(stderr, "    -D         skip self mappings but keep dual mappings\n"); // too confusing to expose to end users
 		fprintf(stderr, "    -S         skip self and dual mappings\n");
 		fprintf(stderr, "    -p FLOAT   threshold to output a mapping [%g]\n", opt.pri_ratio);
+		fprintf(stderr, "    -z INT     Z-drop score [%d]\n", opt.zdrop);
 		fprintf(stderr, "    -x STR     preset (recommended to be applied before other options) []\n");
 		fprintf(stderr, "               ava10k: -Sw5 -L100 -m0 (PacBio/ONT all-vs-all read mapping)\n");
 		fprintf(stderr, "  Input/Output:\n");