Release minimap2-2.0 (r275)

fix self-comparison in index parameter override check

Makefile

@@ -2,7 +2,7 @@ CC=			gcc
 CFLAGS=		-g -Wall -O2 -Wc++-compat
 CPPFLAGS=	-DHAVE_KALLOC
 INCLUDES=	-I.
-OBJS=		kthread.o kalloc.o ksw2_extz2_sse.o ksw2_extd2_sse.o misc.o bseq.o \
+OBJS=		kthread.o kalloc.o ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_ll_sse.o misc.o bseq.o \
 			sketch.o sdust.o index.o chain.o align.o hit.o map.o format.o
 PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite
@@ -51,6 +51,7 @@ index.o: kthread.h bseq.h minimap.h mmpriv.h kvec.h kalloc.h khash.h
 kalloc.o: kalloc.h
 ksw2_extd2_sse.o: ksw2.h kalloc.h
 ksw2_extz2_sse.o: ksw2.h kalloc.h
+ksw2_ll_sse.o: ksw2.h kalloc.h
 main.o: bseq.h minimap.h mmpriv.h
 map.o: kthread.h kvec.h kalloc.h sdust.h mmpriv.h minimap.h bseq.h
 misc.o: minimap.h ksort.h

NEWS.md

@@ -0,0 +1,47 @@
+Release 2.0-r275 (8 August 2017)
+--------------------------------
+
+This release is identical to version 2.0rc1, except the version number. It is
+described and evaluated in the following technical report:
+
+ * Li, H. (2017). Minimap2: fast pairwise alignment for long DNA sequences.
+   [arXiv:1708.01492v1](https://arxiv.org/abs/1708.01492v1).
+
+(2.0: 8 August 2017, r275)
+
+
+
+Release 2.0rc1-r232 (30 July 2017)
+----------------------------------
+
+This release improves the accuracy of long-read alignment and added several
+minor features.
+
+ * Improved mapping quality estimate for short alignments containing few seed
+   hits.
+
+ * Fixed a minor bug that affects the chaining accuracy towards the ends of a
+   chain. Changed the gap cost for chaining to reduce false seeding.
+
+ * Skip potentially wrong seeding and apply dynamic programming more frequently.
+   This slightly increases run time, but greatly reduces false long gaps.
+
+ * Perform local alignment at Z-drop break point to recover potential inversion
+   alignment. Output the SA tag in the SAM format. Added scripts to evaluate
+   mapping accuracy for reads simulated with pbsim.
+
+This release completes features intended for v2.0. No major features will be
+added to the master branch before the final v2.0.
+
+(2.0rc1: 30 July 2017, r232)
+
+
+
+Release r191 (19 July 2017)
+---------------------------
+
+This is the first public release of minimap2, an aligner for long reads and
+assemblies. This release has a few issues and is generally not recommended for
+production uses.
+
+(19 July 2017, r191)

README.md

@@ -31,6 +31,10 @@ and [longISLND][longislnd]), better chaining and the ability to produce CIGAR
 with fast extension alignment (see also [libgaba][gaba] and [ksw2][ksw2]) and
 piece-wise affine gap cost.
 
+If you use minimap2 in your work, please consider to cite:
+
+> Li, H. (2017). Minimap2: fast pairwise alignment for long DNA sequences. [arXiv:1708.01492](https://arxiv.org/abs/1708.01492).
+
 ## Installation
 
 For modern x86-64 CPUs, just type `make` in the source code directory. This
@@ -40,23 +44,64 @@ will run a little slower. At present, minimap2 does not work with non-x86 CPUs
 or ancient CPUs that do not support SSE2. SSE2 is critical to the performance
 of minimap2.
 
+## Algorithm Overview
+
+In the following, minimap2 command line options have a dash ahead and are
+highlighted in bold.
+
+1. Read **-I** [=*4G*] reference bases, extract (**-k**,**-w**)-minimizers and
+   index them in a hash table.
+
+2. Read **-K** [=*200M*] query bases. For each query sequence, do step 3
+   through 7:
+
+3. For each (**-k**,**-w**)-minimizer on the query, check against the reference
+   index. If a reference minimizer is not among the top **-f** [=*2e-4*] most
+   frequent, collect its the occurrences in the reference, which are called
+   *seeds*.
+
+4. Sort seeds by position in the reference. Chain them with dynamic
+   programming. Each chain represents a potential mapping. For read
+   overlapping, report all chains and then go to step 8. For reference mapping,
+   do step 5 through 7:
+
+5. Let *P* be the set of primary mappings, which is an empty set initially. For
+   each chain from the best to the worst according to their chaining scores: if
+   on the query, the chain overlaps with a chain in *P* by **--mask-level**
+   [=*0.5*] or higher fraction of the shorter chain, mark the chain as
+   *secondary* to the chain in *P*; otherwise, add the chain to *P*.
+
+6. Retain all primary mappings. Also retain up to **-N** [=*5*] top secondary
+   mappings if their chaining scores are higher than **-p** [=*0.8*] of their
+   corresponding primary mappings.
+
+7. If alignment is requested, filter out an internal seed if it potentially
+   leads to both a long insertion and a long deletion. Extend from the
+   left-most seed. Perform global alignments between internal seeds.  Split the
+   chain if the accumulative score along the global alignment drops by **-z**
+   [=*400*], disregarding long gaps. Extend from the right-most seed.  Output
+   chains and their alignments.
+
+8. If there are more query sequences in the input, go to step 2 until no more
+   queries are left.
+
+9. If there are more reference sequences, reopen the query file from the start
+   and go to step 1; otherwise stop.
+
 ## Limitations
 
-* At the alignment phase, minimap2 performs global alignments between minimizer
-  hits. If the positions of these minimizer hits are incorrect, the final
-  alignment may be suboptimal or unnecessarily fragmented.
-
-* Minimap2 may produce poor alignments that may need post-filtering. We are
-  still exploring a reliable and consistent way to report good alignments.
+* Minimap2 may produce suboptimal alignments through long low-complexity
+  regions where seed positions may be suboptimal. This should not be a big
+  concern because even the optimal alignment may be wrong in such regions.
 
 * Minimap2 does not work well with Illumina short reads as of now.
 
 * Minimap2 requires SSE2 instructions to compile. It is possible to add
   non-SSE2 support, but it would make minimap2 slower by several times.
 
-In general, minimap2 is a young project with most code written since June,
-2017. It may have bugs and room for improvements. Bug reports and suggestions
-are warmly welcomed.
+In general, minimap2 is a young project with most code written since June, 2017.
+It may have bugs and room for improvements. Bug reports and suggestions are
+warmly welcomed.
 
 
 

align.c

@@ -126,6 +126,12 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) //
 
 static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const int8_t *mat, int w, int flag, ksw_extz_t *ez)
 {
+	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
+		int i;
+		fprintf(stderr, "===> q=(%d,%d), e=(%d,%d), bw=%d, flag=%d, zdrop=%d <===\n", opt->q, opt->q2, opt->e, opt->e2, w, flag, opt->zdrop);
+		for (i = 0; i < tlen; ++i) fputc("ACGTN"[tseq[i]], stderr); fputc('\n', stderr);
+		for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr); fputc('\n', stderr);
+	}
 	if (opt->q == opt->q2 && opt->e == opt->e2)
 		ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, opt->zdrop, flag, ez);
 	else
@@ -158,27 +164,85 @@ static inline void mm_adjust_minier(const mm_idx_t *mi, uint8_t *const qseq0[2],
 	}
 }
 
-static void mm_filter_bad_seeds(int n, mm128_t *a, int max_space, int max_diff)
+static void mm_filter_bad_seeds(void *km, int as1, int cnt1, mm128_t *a, int min_gap, int diff_thres, int max_ext_len, int max_ext_cnt)
 {
-	int i;
-	if (n < 3) return;
-	for (i = 0; i < n - 2; ++i) {
-		int32_t q[3], t[3], gap01, gap12, gap02;
-		t[0] = (int32_t)a[i].x, q[0] = (int32_t)a[i].y;
-		t[1] = (int32_t)a[i+1].x, q[1] = (int32_t)a[i+1].y;
-		t[2] = (int32_t)a[i+2].x, q[2] = (int32_t)a[i+2].y;
-		if (t[2] - t[0] > max_space || q[2] - q[0] > max_space) continue;
-		gap01 = (t[1] - t[0]) - (q[1] - q[0]), gap01 = gap01 > 0? gap01 : -gap01;
-		gap12 = (t[2] - t[1]) - (q[2] - q[1]), gap12 = gap12 > 0? gap12 : -gap12;
-		gap02 = (t[2] - t[0]) - (q[2] - q[0]), gap02 = gap02 > 0? gap02 : -gap02;
-		if (gap01 + gap12 - gap02 > max_diff)
-			a[++i].y |= 1ULL << 41;
+	int max_st, max_en, n, i, k, max, *K;
+	for (i = 1, n = 0; i < cnt1; ++i) { // count the number of gaps longer than min_gap
+		int gap = ((int32_t)a[as1 + i].y - a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - a[as1 + i - 1].x);
+		if (gap < -min_gap || gap > min_gap) ++n;
+	}
+	if (n <= 1) return;
+	K = (int*)kmalloc(km, n * sizeof(int));
+	for (i = 1, n = 0; i < cnt1; ++i) { // store the positions of long gaps
+		int gap = ((int32_t)a[as1 + i].y - a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - a[as1 + i - 1].x);
+		if (gap < -min_gap || gap > min_gap)
+			K[n++] = i;
+	}
+	max = 0, max_st = max_en = -1;
+	for (k = 0;; ++k) { // traverse long gaps
+		int gap, l, n_ins = 0, n_del = 0, qs, rs, max_diff = 0, max_diff_l = -1;
+		if (k == n || k >= max_en) {
+			if (max_en > 0)
+				for (i = K[max_st]; i < K[max_en]; ++i)
+					a[as1 + i].y |= MM_SEED_IGNORE;
+			max = 0, max_st = max_en = -1;
+			if (k == n) break;
+		}
+		i = K[k];
+		gap = ((int32_t)a[as1 + i].y - a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - a[as1 + i - 1].x);
+		if (gap > 0) n_ins += gap;
+		else n_del += -gap;
+		qs = (int32_t)a[as1 + i - 1].y;
+		rs = (int32_t)a[as1 + i - 1].x;
+		for (l = k + 1; l < n && l <= k + max_ext_cnt; ++l) {
+			int j = K[l], diff;
+			if ((int32_t)a[as1 + j].y - qs > max_ext_len || (int32_t)a[as1 + j].x - rs > max_ext_len) break;
+			gap = ((int32_t)a[as1 + j].y - (int32_t)a[as1 + j - 1].y) - (a[as1 + j].x - a[as1 + j - 1].x);
+			if (gap > 0) n_ins += gap;
+			else n_del += -gap;
+			diff = n_ins + n_del - abs(n_ins - n_del);
+			if (max_diff < diff)
+				max_diff = diff, max_diff_l = l;
+		}
+		if (max_diff > diff_thres && max_diff > max)
+			max = max_diff, max_st = k, max_en = max_diff_l;
+	}
+	kfree(km, K);
+}
+
+static void mm_fix_bad_ends(const mm_reg1_t *r, const mm128_t *a, int bw, int32_t *as, int32_t *cnt)
+{
+	int32_t i, l;
+	*as = r->as, *cnt = r->cnt;
+	if (r->cnt < 3) return;
+	l = a[r->as].y >> 32 & 0xff;
+	for (i = r->as + 1; i < r->as + r->cnt - 1; ++i) {
+		int32_t lq, lr, min, max;
+		lr = (int32_t)a[i].x - (int32_t)a[i-1].x;
+		lq = (int32_t)a[i].y - (int32_t)a[i-1].y;
+		min = lr < lq? lr : lq;
+		max = lr > lq? lr : lq;
+		if (max - min > l >> 1) *as = i;
+		l += min;
+		if (l >= bw << 1) break;
+	}
+	*cnt = r->as + r->cnt - *as;
+	l = a[r->as + r->cnt - 1].y >> 32 & 0xff;
+	for (i = r->as + r->cnt - 2; i > *as; --i) {
+		int32_t lq, lr, min, max;
+		lr = (int32_t)a[i+1].x - (int32_t)a[i].x;
+		lq = (int32_t)a[i+1].y - (int32_t)a[i].y;
+		min = lr < lq? lr : lq;
+		max = lr > lq? lr : lq;
+		if (max - min > l >> 1) *cnt = i + 1 - *as;
+		l += min;
+		if (l >= bw) break;
 	}
 }
 
 static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], mm_reg1_t *r, mm_reg1_t *r2, mm128_t *a, ksw_extz_t *ez)
 {
-	int32_t rid = a[r->as].x<<1>>33, rev = a[r->as].x>>63;
+	int32_t rid = a[r->as].x<<1>>33, rev = a[r->as].x>>63, as1, cnt1;
 	uint8_t *tseq, *qseq;
 	int32_t i, l, bw, dropped = 0, rs0, re0, qs0, qe0;
 	int32_t rs, re, qs, qe;
@@ -190,13 +254,14 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	bw = (int)(opt->bw * 1.5 + 1.);
 
 	r2->cnt = 0;
-	mm_adjust_minier(mi, qseq0, &a[r->as], &rs, &qs);
-	mm_adjust_minier(mi, qseq0, &a[r->as + r->cnt - 1], &re, &qe);
-	mm_filter_bad_seeds(r->cnt, &a[r->as], opt->max_gap, 100);
+	mm_fix_bad_ends(r, a, opt->bw, &as1, &cnt1);
+	mm_filter_bad_seeds(km, as1, cnt1, a, 10, 40, opt->max_gap>>1, 10);
+	mm_adjust_minier(mi, qseq0, &a[as1], &rs, &qs);
+	mm_adjust_minier(mi, qseq0, &a[as1 + cnt1 - 1], &re, &qe);
 
 	// compute rs0 and qs0
-	if (r->split && r->as > 0) {
-		mm_adjust_minier(mi, qseq0, &a[r->as-1], &rs0, &qs0);
+	if (r->split && as1 > 0) {
+		mm_adjust_minier(mi, qseq0, &a[as1-1], &rs0, &qs0);
 	} else {
 		if (qs > 0 && rs > 0) { // actually this is always true
 			l = qs < opt->max_gap? qs : opt->max_gap;
@@ -238,13 +303,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	re1 = rs, qe1 = qs;
 	assert(qs1 >= 0 && rs1 >= 0);
 
-	for (i = 1; i < r->cnt; ++i) { // gap filling
-		if (a[r->as+i].y>>41&1) continue;
-		mm_adjust_minier(mi, qseq0, &a[r->as + i], &re, &qe);
+	for (i = 1; i < cnt1; ++i) { // gap filling
+		if ((a[as1+i].y & (MM_SEED_IGNORE|MM_SEED_TANDEM)) && i != cnt1 - 1) continue;
+		mm_adjust_minier(mi, qseq0, &a[as1 + i], &re, &qe);
 		re1 = re, qe1 = qe;
-		if (i == r->cnt - 1 || (a[r->as+i].y>>40&1) || qe - qs >= opt->min_ksw_len || re - rs >= opt->min_ksw_len) {
+		if (i == cnt1 - 1 || (a[as1+i].y&MM_SEED_LONG_JOIN) || (qe - qs >= opt->min_ksw_len && re - rs >= opt->min_ksw_len)) {
 			int bw1 = bw;
-			if (a[r->as+i].y>>40&1)
+			if (a[as1+i].y & MM_SEED_LONG_JOIN)
 				bw1 = qe - qs > re - rs? qe - qs : re - rs;
 			qseq = &qseq0[rev][qs];
 			mm_idx_getseq(mi, rid, rs, re, tseq);
@@ -256,13 +321,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			if (ez->zdropped) { // truncated by Z-drop; TODO: sometimes Z-drop kicks in because the next seed placement is wrong. This can be fixed in principle.
 				int j;
 				for (j = i - 1; j >= 0; --j)
-					if ((int32_t)a[r->as + j].x < re + ez->max_t)
+					if ((int32_t)a[as1 + j].x < re + ez->max_t)
 						break;
 				dropped = 1;
 				r->p->dp_score += ez->max;
 				re1 = rs + (ez->max_t + 1);
 				qe1 = qs + (ez->max_q + 1);
-				if (r->cnt - (j + 1) >= opt->min_cnt)
+				if (cnt1 - (j + 1) >= opt->min_cnt)
 					mm_split_reg(r, r2, j + 1, qlen, a);
 				break;
 			} else r->p->dp_score += ez->score;
@@ -296,10 +361,68 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	kfree(km, tseq);
 }
 
+static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez)
+{
+	int tl, ql, score, ret = 0, q_off, t_off;
+	uint8_t *tseq, *qseq;
+	int8_t mat[25];
+	void *qp;
+
+	memset(r_inv, 0, sizeof(mm_reg1_t));
+	if (!(r1->split&1) || !(r2->split&2)) return 0;
+	if (r1->id != r1->parent && r1->parent != MM_PARENT_TMP_PRI) return 0;
+	if (r2->id != r2->parent && r2->parent != MM_PARENT_TMP_PRI) return 0;
+	if (r1->rid != r2->rid || r1->rev != r2->rev) return 0;
+	ql = r2->qs - r1->qe;
+	tl = r2->rs - r1->re;
+	if (ql < opt->min_chain_score || ql > opt->max_gap) return 0;
+	if (tl < opt->min_chain_score || tl > opt->max_gap) return 0;
+
+	ksw_gen_simple_mat(5, mat, opt->a, opt->b);
+	tseq = (uint8_t*)kmalloc(km, tl);
+	mm_idx_getseq(mi, r1->rid, r1->re, r2->rs, tseq);
+	qseq = &qseq0[!r1->rev][qlen - r2->qs];
+
+	mm_seq_rev(ql, qseq);
+	mm_seq_rev(tl, tseq);
+	qp = ksw_ll_qinit(km, 2, ql, qseq, 5, mat);
+	score = ksw_ll_i16(qp, tl, tseq, opt->q, opt->e, &q_off, &t_off);
+	kfree(km, qp);
+	mm_seq_rev(ql, qseq);
+	mm_seq_rev(tl, tseq);
+	if (score < opt->min_dp_max) goto end_align1_inv;
+	q_off = ql - (q_off + 1), t_off = tl - (t_off + 1);
+	mm_align_pair(km, opt, ql - q_off, qseq + q_off, tl - t_off, tseq + t_off, mat, (int)(opt->bw * 1.5), KSW_EZ_EXTZ_ONLY, ez);
+	if (ez->n_cigar == 0) goto end_align1_inv; // should never be here
+	mm_append_cigar(r_inv, ez->n_cigar, ez->cigar);
+	r_inv->p->dp_score = ez->max;
+	mm_update_extra(r_inv->p, qseq + q_off, tseq + t_off, mat, opt->q, opt->e);
+	r_inv->id = -1;
+	r_inv->parent = MM_PARENT_UNSET;
+	r_inv->inv = 1;
+	r_inv->rev = !r1->rev;
+	r_inv->qs = r1->qe + q_off, r_inv->qe = r_inv->qs + ez->max_q + 1;
+	r_inv->rs = r1->re + t_off, r_inv->re = r_inv->rs + ez->max_t + 1;
+	ret = 1;
+end_align1_inv:
+	kfree(km, tseq);
+	return ret;
+}
+
+static inline mm_reg1_t *mm_insert_reg(const mm_reg1_t *r, int i, int *n_regs, mm_reg1_t *regs)
+{
+	regs = (mm_reg1_t*)realloc(regs, (*n_regs + 1) * sizeof(mm_reg1_t));
+	if (i + 1 != *n_regs)
+		memmove(&regs[i + 2], &regs[i + 1], sizeof(mm_reg1_t) * (*n_regs - i - 1));
+	regs[i + 1] = *r;
+	++*n_regs;
+	return regs;
+}
+
 mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a)
 {
 	extern unsigned char seq_nt4_table[256];
-	int32_t i, r, n_regs = *n_regs_;
+	int32_t i, n_regs = *n_regs_;
 	uint8_t *qseq0[2];
 	ksw_extz_t ez;
 
@@ -313,15 +436,13 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 
 	// align through seed hits
 	memset(&ez, 0, sizeof(ksw_extz_t));
-	for (r = 0; r < n_regs; ++r) {
+	for (i = 0; i < n_regs; ++i) {
 		mm_reg1_t r2;
-		mm_align1(km, opt, mi, qlen, qseq0, &regs[r], &r2, a, &ez);
-		if (r2.cnt > 0) {
-			regs = (mm_reg1_t*)realloc(regs, (n_regs + 1) * sizeof(mm_reg1_t)); // this should be very rare
-			if (r + 1 != n_regs)
-				memmove(&regs[r + 2], &regs[r + 1], sizeof(mm_reg1_t) * (n_regs - r - 1));
-			regs[r + 1] = r2;
-			++n_regs;
+		mm_align1(km, opt, mi, qlen, qseq0, &regs[i], &r2, a, &ez);
+		if (r2.cnt > 0) regs = mm_insert_reg(&r2, i, &n_regs, regs);
+		if (i > 0 && mm_align1_inv(km, opt, mi, qlen, qseq0, &regs[i-1], &regs[i], &r2, &ez)) {
+			regs = mm_insert_reg(&r2, i, &n_regs, regs);
+			++i; // skip the inserted INV alignment
 		}
 	}
 	*n_regs_ = n_regs;

chain.c

@@ -7,9 +7,9 @@
 
 static const char LogTable256[256] = {
 #define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
-    -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
-    LT(4), LT(5), LT(5), LT(6), LT(6), LT(6), LT(6),
-    LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7)
+	-1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
+	LT(4), LT(5), LT(5), LT(6), LT(6), LT(6), LT(6),
+	LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7)
 };
 
 static inline int ilog2_32(uint32_t v)
@@ -21,22 +21,27 @@ static inline int ilog2_32(uint32_t v)
 
 int mm_chain_dp(int max_dist, int bw, int max_skip, int min_cnt, int min_sc, int64_t n, mm128_t *a, uint64_t **_u, void *km)
 { // TODO: make sure this works when n has more than 32 bits
-	int32_t st = 0, j, k, *f, *p, *t, *v, n_u, n_v;
-	int64_t i;
-	uint64_t *u, *u2;
+	int32_t st = 0, k, *f, *p, *t, *v, n_u, n_v;
+	int64_t i, j;
+	uint64_t *u, *u2, sum_qspan = 0;
+	float avg_qspan;
 	mm128_t *b, *w;
 
 	if (_u) *_u = 0;
 	f = (int32_t*)kmalloc(km, n * 4);
 	p = (int32_t*)kmalloc(km, n * 4);
 	t = (int32_t*)kmalloc(km, n * 4);
+	v = (int32_t*)kmalloc(km, n * 4);
 	memset(t, 0, n * 4);
 
+	for (i = 0; i < n; ++i) sum_qspan += a[i].y>>32&0xff;
+	avg_qspan = (float)sum_qspan / n;
+
 	// fill the score and backtrack arrays
 	for (i = 0; i < n; ++i) {
 		uint64_t ri = a[i].x;
 		int32_t qi = (int32_t)a[i].y, q_span = a[i].y>>32&0xff; // NB: only 8 bits of span is used!!!
-		int32_t max_f = -INT32_MAX, max_j = -1, n_skip = 0, min_d;
+		int32_t max_f = q_span, max_j = -1, n_skip = 0, min_d, max_f_past = -INT32_MAX;
 		while (st < i && ri - a[st].x > max_dist) ++st;
 		for (j = i - 1; j >= st; --j) {
 			int64_t dr = ri - a[j].x;
@@ -44,10 +49,10 @@ int mm_chain_dp(int max_dist, int bw, int max_skip, int min_cnt, int min_sc, int
 			if (dr == 0 || dq <= 0 || dq > max_dist) continue;
 			dd = dr > dq? dr - dq : dq - dr;
 			if (dd > bw) continue;
+			max_f_past = max_f_past > f[j]? max_f_past : f[j];
 			min_d = dq < dr? dq : dr;
 			sc = min_d > q_span? q_span : dq < dr? dq : dr;
-			sc -= dd? ilog2_32(dd) * 2 : 0;
-			if (min_d > q_span) sc -= ilog2_32(min_d) / 2;
+			sc -= (int)(dd * .01 * avg_qspan) + (ilog2_32(dd)>>1);
 			sc += f[j];
 			if (sc > max_f) {
 				max_f = sc, max_j = j;
@@ -58,8 +63,7 @@ int mm_chain_dp(int max_dist, int bw, int max_skip, int min_cnt, int min_sc, int
 			}
 			if (p[j] >= 0) t[p[j]] = i;
 		}
-		if (max_j >= 0) f[i] = max_f, p[i] = max_j;
-		else f[i] = q_span, p[i] = -1;
+		f[i] = max_f, p[i] = max_j, v[i] = max_f_past; // v[] keeps the max score in the previous chain
 	}
 
 	// find the ending positions of chains
@@ -67,16 +71,21 @@ int mm_chain_dp(int max_dist, int bw, int max_skip, int min_cnt, int min_sc, int
 	for (i = 0; i < n; ++i)
 		if (p[i] >= 0) t[p[i]] = 1;
 	for (i = n_u = 0; i < n; ++i)
-		if (t[i] == 0 && f[i] >= min_sc)
+		if (t[i] == 0 && v[i] >= min_sc)
 			++n_u;
 	if (n_u == 0) {
-		kfree(km, f); kfree(km, p); kfree(km, t);
+		kfree(km, f); kfree(km, p); kfree(km, t); kfree(km, v);
 		return 0;
 	}
 	u = (uint64_t*)kmalloc(km, n_u * 8);
-	for (i = n_u = 0; i < n; ++i)
-		if (t[i] == 0 && f[i] >= min_sc)
-			u[n_u++] = (uint64_t)f[i] << 32 | i;
+	for (i = n_u = 0; i < n; ++i) {
+		if (t[i] == 0 && v[i] >= min_sc) {
+			j = i;
+			while (j >= 0 && f[j] < v[j]) j = p[j]; // find the point that maximizes f[]
+			if (j < 0) j = i; // TODO: this should really be assert(j>=0)
+			u[n_u++] = (uint64_t)f[j] << 32 | j;
+		}
+	}
 	radix_sort_64(u, u + n_u);
 	for (i = 0; i < n_u>>1; ++i) { // reverse, s.t. the highest scoring chain is the first
 		uint64_t t = u[i];
@@ -85,7 +94,6 @@ int mm_chain_dp(int max_dist, int bw, int max_skip, int min_cnt, int min_sc, int
 
 	// backtrack
 	memset(t, 0, n * 4);
-	v = (int32_t*)kmalloc(km, n * 4);
 	for (i = n_v = k = 0; i < n_u; ++i) { // starting from the highest score
 		int32_t n_v0 = n_v, k0 = k;
 		j = (int32_t)u[i];

format.c

@@ -56,7 +56,8 @@ static void mm_sprintf_lite(kstring_t *s, const char *fmt, ...)
 
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 {
-	mm_sprintf_lite(s, "\tcm:i:%d\ts1:i:%d", r->cnt, r->score);
+	int type = r->inv? 'I' : r->id == r->parent? 'P' : 'S';
+	mm_sprintf_lite(s, "\ttp:A:%c\tcm:i:%d\ts1:i:%d", type, r->cnt, r->score);
 	if (r->parent == r->id) mm_sprintf_lite(s, "\ts2:i:%d", r->subsc);
 	if (r->split) mm_sprintf_lite(s, "\tzd:i:%d", r->split);
 	if (r->p) mm_sprintf_lite(s, "\tNM:i:%d\tms:i:%d\tAS:i:%d\tnn:i:%d", r->p->n_diff, r->p->dp_max, r->p->dp_score, r->p->n_ambi);
@@ -105,7 +106,7 @@ static void sam_write_sq(kstring_t *s, char *seq, int l, int rev, int comp)
 	} else str_copy(s, seq, seq + l);
 }
 
-void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r)
+void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int n_regs, const mm_reg1_t *regs)
 {
 	int flag = 0;
 	s->l = 0;
@@ -120,7 +121,7 @@ void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m
 		if (r->parent != r->id) flag |= 0x100;
 		else if (!r->sam_pri) flag |= 0x800;
 		mm_sprintf_lite(s, "%s\t%d\t%s\t%d\t%d\t", t->name, flag, mi->seq[r->rid].name, r->rs+1, r->mapq);
-		if (r->p) { // TODO: using hard clippings
+		if (r->p) { // actually this should always be true for SAM output
 			uint32_t k, clip_len = r->rev? t->l_seq - r->qe : r->qs;
 			int clip_char = (flag&0x800)? 'H' : 'S';
 			if (clip_len) mm_sprintf_lite(s, "%d%c", clip_len, clip_char);
@@ -144,6 +145,31 @@ void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m
 			else mm_sprintf_lite(s, "*");
 		}
 		write_tags(s, r);
+		if (r->parent == r->id && r->p && n_regs > 1 && regs && r >= regs && r - regs < n_regs) { // supplementary aln may exist
+			int i, n_sa = 0; // n_sa: number of SA fields
+			for (i = 0; i < n_regs; ++i)
+				if (i != r - regs && regs[i].parent == regs[i].id && regs[i].p)
+					++n_sa;
+			if (n_sa > 0) {
+				mm_sprintf_lite(s, "\tSA:Z:");
+				for (i = 0; i < n_regs; ++i) {
+					const mm_reg1_t *q = &regs[i];
+					int l_M, l_I = 0, l_D = 0, clip5 = 0, clip3 = 0;
+					if (r == q || q->parent != q->id || q->p == 0) continue;
+					if (q->qe - q->qs < q->re - q->rs) l_M = q->qe - q->qs, l_D = (q->re - q->rs) - l_M;
+					else l_M = q->re - q->rs, l_I = (q->qe - q->qs) - l_M;
+					clip5 = q->rev? t->l_seq - q->qe : q->qs;
+					clip3 = q->rev? q->qs : t->l_seq - q->qe;
+					mm_sprintf_lite(s, "%s,%d,%c,", mi->seq[q->rid].name, q->rs+1, "+-"[q->rev]);
+					if (clip5) mm_sprintf_lite(s, "%dS", clip5);
+					if (l_M) mm_sprintf_lite(s, "%dM", l_M);
+					if (l_I) mm_sprintf_lite(s, "%dI", l_I);
+					if (l_D) mm_sprintf_lite(s, "%dD", l_D);
+					if (clip3) mm_sprintf_lite(s, "%dS", clip3);
+					mm_sprintf_lite(s, ",%d,%d;", q->mapq, q->p->n_diff);
+				}
+			}
+		}
 	}
 	s->s[s->l] = 0; // we always have room for an extra byte (see str_enlarge)
 }

hit.c

@@ -125,11 +125,11 @@ void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r)
 	aux = (uint64_t*)kmalloc(km, n * 8);
 	t = (mm_reg1_t*)kmalloc(km, n * sizeof(mm_reg1_t));
 	for (i = n_aux = 0; i < n; ++i) {
-		if (r[i].cnt > 0) { // squeeze out elements with cnt==0 (soft deleted)
+		if (r[i].inv || r[i].cnt > 0) { // squeeze out elements with cnt==0 (soft deleted)
 			assert(r[i].p);
 			aux[n_aux++] = (uint64_t)r[i].p->dp_max << 32 | i;
 		} else if (r[i].p) {
-			kfree(km, r[i].p);
+			free(r[i].p);
 			r[i].p = 0;
 		}
 	}
@@ -177,13 +177,14 @@ void mm_sync_regs(void *km, int n_regs, mm_reg1_t *regs) // keep mm_reg1_t::{id,
 	mm_set_sam_pri(n_regs, regs);
 }
 
-void mm_select_sub(void *km, float mask_level, float pri_ratio, int best_n, int *n_, mm_reg1_t *r)
+void mm_select_sub(void *km, float mask_level, float pri_ratio, int min_diff, int best_n, int *n_, mm_reg1_t *r)
 {
 	if (pri_ratio > 0.0f && *n_ > 0) {
 		int i, k, n = *n_, n_2nd = 0;
 		for (i = k = 0; i < n; ++i)
 			if (r[i].parent == i) r[k++] = r[i];
-			else if (r[i].score >= r[r[i].parent].score * pri_ratio && n_2nd++ < best_n) r[k++] = r[i];
+			else if ((r[i].score >= r[r[i].parent].score * pri_ratio || r[i].score + min_diff >= r[r[i].parent].score) && n_2nd++ < best_n)
+				r[k++] = r[i];
 			else if (r[i].p) free(r[i].p);
 		if (k != n) mm_sync_regs(km, k, r); // removing hits requires sync()
 		*n_ = k;
@@ -196,7 +197,7 @@ void mm_filter_regs(void *km, const mm_mapopt_t *opt, int *n_regs, mm_reg1_t *re
 	for (i = k = 0; i < *n_regs; ++i) {
 		mm_reg1_t *r = &regs[i];
 		int flt = 0;
-		if (r->cnt < opt->min_cnt) flt = 1;
+		if (!r->inv && r->cnt < opt->min_cnt) flt = 1;
 		if (r->p) {
 			if (r->p->blen - r->p->n_ambi - r->p->n_diff < opt->min_chain_score) flt = 1;
 			else if (r->p->dp_max < opt->min_dp_max) flt = 1;
@@ -265,7 +266,7 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
 		if (r1->re - r1->rs < max_gap>>1 || r1->qe - r1->qs < max_gap>>1) continue;
 
 		// all conditions satisfied; join
-		a[r1->as].y |= 1ULL<<40;
+		a[r1->as].y |= MM_SEED_LONG_JOIN;
 		r0->cnt += r1->cnt, r0->score += r1->score;
 		mm_reg_set_coor(r0, qlen, a);
 		r1->cnt = 0;
@@ -287,17 +288,22 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
 	}
 }
 
-void mm_set_mapq(int n_regs, mm_reg1_t *regs)
+void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc)
 {
+	static const float q_coef = 30.0f;
 	int i;
 	for (i = 0; i < n_regs; ++i) {
 		mm_reg1_t *r = &regs[i];
-		if (r->parent == r->id) {
-			int mapq;
+		if (r->inv) {
+			r->mapq = 0;
+		} else if (r->parent == r->id) {
+			int mapq, subsc;
+			float pen_cm = r->cnt >= 10? 1.0f : 0.1f * r->cnt;
+			subsc = r->subsc > min_chain_sc? r->subsc : min_chain_sc;
 			if (r->p && r->p->dp_max2 > 0 && r->p->dp_max > 0) {
 				float identity = (float)(r->p->blen - r->p->n_diff - r->p->n_ambi) / (r->p->blen - r->p->n_ambi);
-				mapq = (int)(identity * 30.0 * (1. - (float)r->p->dp_max2 * r->subsc / r->p->dp_max / r->score) * logf(r->score));
-			} else mapq = (int)(30.0 * (1. - (float)r->subsc / r->score) * logf(r->score));
+				mapq = (int)(identity * pen_cm * q_coef * (1. - (float)r->p->dp_max2 * subsc / r->p->dp_max / r->score) * logf(r->score));
+			} else mapq = (int)(pen_cm * q_coef * (1. - (float)subsc / r->score) * logf(r->score));
 			mapq = mapq > 0? mapq : 0;
 			r->mapq = mapq < 60? mapq : 60;
 		} else r->mapq = 0;

index.c

@@ -25,6 +25,7 @@ mm_idx_t *mm_idx_init(int w, int k, int b, int is_hpc)
 	mi = (mm_idx_t*)calloc(1, sizeof(mm_idx_t));
 	mi->w = w, mi->k = k, mi->b = b, mi->is_hpc = is_hpc;
 	mi->B = (mm_idx_bucket_t*)calloc(1<<b, sizeof(mm_idx_bucket_t));
+	if (!(mm_dbg_flag & 1)) mi->km = km_init();
 	return mi;
 }
 
@@ -37,9 +38,12 @@ void mm_idx_destroy(mm_idx_t *mi)
 		free(mi->B[i].a.a);
 		kh_destroy(idx, (idxhash_t*)mi->B[i].h);
 	}
-	for (i = 0; i < mi->n_seq; ++i)
-		free(mi->seq[i].name);
-	free(mi->seq); free(mi->B); free(mi->S); free(mi);
+	if (!mi->km) {
+		for (i = 0; i < mi->n_seq; ++i)
+			free(mi->seq[i].name);
+		free(mi->seq);
+	} else km_destroy(mi->km);
+	free(mi->B); free(mi->S); free(mi);
 }
 
 const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n)
@@ -228,7 +232,7 @@ static void *worker_pipeline(void *shared, int step, void *in)
 			old_m = p->mi->n_seq, m = p->mi->n_seq + s->n_seq;
 			kroundup32(m); kroundup32(old_m);
 			if (old_m != m)
-				p->mi->seq = (mm_idx_seq_t*)realloc(p->mi->seq, m * sizeof(mm_idx_seq_t));
+				p->mi->seq = (mm_idx_seq_t*)krealloc(p->mi->km, p->mi->seq, m * sizeof(mm_idx_seq_t));
 			// make room for p->mi->S
 			for (i = 0, sum_len = 0; i < s->n_seq; ++i) sum_len += s->seq[i].l_seq;
 			old_max_len = (p->sum_len + 7) / 8;
@@ -244,7 +248,8 @@ static void *worker_pipeline(void *shared, int step, void *in)
 				uint32_t j;
 				if (p->keep_name) {
 					assert(strlen(s->seq[i].name) <= 254); // a long query name breaks BAM
-					seq->name = strdup(s->seq[i].name);
+					seq->name = (char*)kmalloc(p->mi->km, strlen(s->seq[i].name) + 1);
+					strcpy(seq->name, s->seq[i].name);
 				} else seq->name = 0;
 				seq->len = s->seq[i].l_seq;
 				seq->offset = p->sum_len;
@@ -364,12 +369,12 @@ mm_idx_t *mm_idx_load(FILE *fp)
 	if (fread(x, 4, 5, fp) != 5) return 0;
 	mi = mm_idx_init(x[0], x[1], x[2], x[4]);
 	mi->n_seq = x[3];
-	mi->seq = (mm_idx_seq_t*)calloc(mi->n_seq, sizeof(mm_idx_seq_t));
+	mi->seq = (mm_idx_seq_t*)kcalloc(mi->km, mi->n_seq, sizeof(mm_idx_seq_t));
 	for (i = 0; i < mi->n_seq; ++i) {
 		uint8_t l;
 		mm_idx_seq_t *s = &mi->seq[i];
 		fread(&l, 1, 1, fp);
-		s->name = (char*)malloc(l + 1);
+		s->name = (char*)kmalloc(mi->km, l + 1);
 		fread(s->name, 1, l, fp);
 		s->name[l] = 0;
 		fread(&s->len, 4, 1, fp);

ksw2.h

@@ -39,8 +39,8 @@ typedef struct {
  * @param mat       m*m scoring mattrix in one-dimension array
  * @param gapo      gap open penalty; a gap of length l cost "-(gapo+l*gape)"
  * @param gape      gap extension penalty
- * @param w         band width
- * @param zdrop     off-diagonal drop-off to stop extension (positive)
+ * @param w         band width (<0 to disable)
+ * @param zdrop     off-diagonal drop-off to stop extension (positive; <0 to disable)
  * @param flag      flag (see KSW_EZ_* macros)
  * @param ez        (out) scores and cigar
  */
@@ -53,6 +53,8 @@ void ksw_extd(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t
 void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
 				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int flag, ksw_extz_t *ez);
 
+void ksw_extf2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t mch, int8_t mis, int8_t e, int w, int xdrop, ksw_extz_t *ez);
+
 /**
  * Global alignment
  *
@@ -67,6 +69,9 @@ int ksw_gg(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *ta
 int ksw_gg2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
 int ksw_gg2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_);
 
+void *ksw_ll_qinit(void *km, int size, int qlen, const uint8_t *query, int m, const int8_t *mat);
+int ksw_ll_i16(void *q, int tlen, const uint8_t *target, int gapo, int gape, int *qe, int *te);
+
 #ifdef __cplusplus
 }
 #endif
@@ -101,16 +106,26 @@ static inline uint32_t *ksw_push_cigar(void *km, int *n_cigar, int *m_cigar, uin
 //   bit 0-2: which type gets the max - 0 for H, 1 for E, 2 for F, 3 for \tilde{E} and 4 for \tilde{F}
 //   bit 3/0x08: 1 if a continuation on the E state (bit 5/0x20 for a continuation on \tilde{E})
 //   bit 4/0x10: 1 if a continuation on the F state (bit 6/0x40 for a continuation on \tilde{F})
-static inline void ksw_backtrack(void *km, int is_rot, int is_rev, const uint8_t *p, const int *off, int n_col, int i0, int j0, int *m_cigar_, int *n_cigar_, uint32_t **cigar_)
+static inline void ksw_backtrack(void *km, int is_rot, int is_rev, const uint8_t *p, const int *off, const int *off_end, int n_col, int i0, int j0, int *m_cigar_, int *n_cigar_, uint32_t **cigar_)
 { // p[] - lower 3 bits: which type gets the max; bit
 	int n_cigar = 0, m_cigar = *m_cigar_, i = i0, j = j0, r, state = 0;
 	uint32_t *cigar = *cigar_, tmp;
 	while (i >= 0 && j >= 0) { // at the beginning of the loop, _state_ tells us which state to check
-		if (is_rot) r = i + j, tmp = p[r * n_col + i - off[r]];
-		else tmp = p[i * n_col + j - off[i]];
+		int force_state = -1;
+		if (is_rot) {
+			r = i + j;
+			if (i < off[r]) force_state = 2;
+			if (off_end && i > off_end[r]) force_state = 1;
+			tmp = force_state < 0? p[r * n_col + i - off[r]] : 0;
+		} else {
+			if (j < off[i]) force_state = 2;
+			if (off_end && j > off_end[i]) force_state = 1;
+			tmp = force_state < 0? p[i * n_col + j - off[i]] : 0;
+		}
 		if (state == 0) state = tmp & 7; // if requesting the H state, find state one maximizes it.
 		else if (!(tmp >> (state + 2) & 1)) state = 0; // if requesting other states, _state_ stays the same if it is a continuation; otherwise, set to H
 		if (state == 0) state = tmp & 7; // TODO: probably this line can be merged into the "else if" line right above; not 100% sure
+		if (force_state >= 0) state = force_state;
 		if (state == 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 0, 1), --i, --j; // match
 		else if (state == 1 || state == 3) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i; // deletion
 		else cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, 1), --j; // insertion

ksw2_extd2_sse.c

@@ -42,7 +42,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	a2= _mm_sub_epi8(a2, tmp); \
 	b2= _mm_sub_epi8(b2, tmp);
 
-	int r, t, qe = q + e, n_col_, *off = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc, long_thres, long_diff;
+	int r, t, qe = q + e, n_col_, *off = 0, *off_end = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc, long_thres, long_diff;
 	int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX);
 	int32_t *H = 0, H0 = 0, last_H0_t = 0;
 	uint8_t *qr, *sf, *mem, *mem2 = 0;
@@ -96,7 +96,8 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	if (with_cigar) {
 		mem2 = (uint8_t*)kmalloc(km, ((qlen + tlen - 1) * n_col_ + 1) * 16);
 		p = (__m128i*)(((size_t)mem2 + 15) >> 4 << 4);
-		off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int));
+		off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2);
+		off_end = off + qlen + tlen - 1;
 	}
 
 	for (t = 0; t < qlen; ++t) qr[t] = query[qlen - 1 - t];
@@ -105,7 +106,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) {
 		int st = 0, en = tlen - 1, st0, en0, st_, en_;
 		int8_t x1, x21, v1;
-		uint8_t *qrr = qr + (qlen - 1 - r), p_en0 = 0;
+		uint8_t *qrr = qr + (qlen - 1 - r);
 		int8_t *u8 = (int8_t*)u, *v8 = (int8_t*)v, *x8 = (int8_t*)x, *x28 = (int8_t*)x2;
 		__m128i x1_, x21_, v1_;
 		// find the boundaries
@@ -199,18 +200,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 			}
 		} else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment
 			__m128i *pr = p + r * n_col_ - st_;
-			off[r] = st;
-			if (en0 < r && en0 < tlen - 1) { // to avoid backtracking out of the band; this assumes a fixed band
-				int8_t a, a2, z = ((uint8_t*)s)[en0];
-				a = x8[en0-1] + v8[en0-1];
-				p_en0 = a > z? 1 : 0;
-				z = a > z? a : z;
-				p_en0 |= a - (z - q) > 0? 1<<4 : 0;
-				a2 = x28[en0-1] + v8[en0-1];
-				p_en0 = a2 > z? 3 : p_en0;
-				z = a2 > z? a2 : z;
-				p_en0 |= a2 - (z - q2) > 0? 1<<6 : 0;
-			}
+			off[r] = st, off_end[r] = en;
 			for (t = st_; t <= en_; ++t) {
 				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp;
 				__dp_code_block1;
@@ -257,18 +247,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 			}
 		} else { // gap right-alignment
 			__m128i *pr = p + r * n_col_ - st_;
-			off[r] = st;
-			if (en0 < r && en0 < tlen - 1) { // to avoid backtracking out of the band; this assumes a fixed band
-				int8_t a, a2, z = ((uint8_t*)s)[en0];
-				a = x8[en0-1] + v8[en0-1];
-				p_en0 = a >= z? 1 : 0;
-				z = a >= z? a : z;
-				p_en0 |= a - (z - q) >= 0? 1<<4 : 0;
-				a2 = x28[en0-1] + v8[en0-1];
-				p_en0 = a2 >= z? 3 : p_en0;
-				z = a2 >= z? a2 : z;
-				p_en0 |= a2 - (z - q2) >= 0? 1<<6 : 0;
-			}
+			off[r] = st, off_end[r] = en;
 			for (t = st_; t <= en_; ++t) {
 				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp;
 				__dp_code_block1;
@@ -314,7 +293,6 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				_mm_store_si128(&pr[t], d);
 			}
 		}
-		if (with_cigar && en0 < r && en0 < tlen - 1) ((uint8_t*)(p + r * n_col_))[en0 - st] = p_en0;
 		if (!approx_max) { // find the exact max with a 32-bit score array
 			int32_t max_H, max_t;
 			// compute H[], max_H and max_t
@@ -384,9 +362,9 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	if (with_cigar) { // backtrack
 		int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR);
 		if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY))
-			ksw_backtrack(km, 1, rev_cigar, (uint8_t*)p, off, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+			ksw_backtrack(km, 1, rev_cigar, (uint8_t*)p, off, off_end, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
 		else if (ez->max_t >= 0 && ez->max_q >= 0)
-			ksw_backtrack(km, 1, rev_cigar, (uint8_t*)p, off, n_col_*16, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+			ksw_backtrack(km, 1, rev_cigar, (uint8_t*)p, off, off_end, n_col_*16, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
 		kfree(km, mem2); kfree(km, off);
 	}
 }

ksw2_extz2_sse.c

@@ -33,7 +33,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	a = _mm_sub_epi8(a, z); \
 	b = _mm_sub_epi8(b, z);
 
-	int r, t, qe = q + e, n_col_, *off = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc;
+	int r, t, qe = q + e, n_col_, *off = 0, *off_end = 0, tlen_, qlen_, last_st, last_en, wl, wr, max_sc, min_sc;
 	int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX);
 	int32_t *H = 0, H0 = 0, last_H0_t = 0;
 	uint8_t *qr, *sf, *mem, *mem2 = 0;
@@ -76,7 +76,8 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	if (with_cigar) {
 		mem2 = (uint8_t*)kmalloc(km, ((qlen + tlen - 1) * n_col_ + 1) * 16);
 		p = (__m128i*)(((size_t)mem2 + 15) >> 4 << 4);
-		off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int));
+		off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2);
+		off_end = off + qlen + tlen - 1;
 	}
 
 	for (t = 0; t < qlen; ++t) qr[t] = query[qlen - 1 - t];
@@ -85,7 +86,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	for (r = 0, last_st = last_en = -1; r < qlen + tlen - 1; ++r) {
 		int st = 0, en = tlen - 1, st0, en0, st_, en_;
 		int8_t x1, v1;
-		uint8_t *qrr = qr + (qlen - 1 - r), *u8 = (uint8_t*)u, *v8 = (uint8_t*)v, *x8 = (uint8_t*)x, p_en0 = 0;
+		uint8_t *qrr = qr + (qlen - 1 - r), *u8 = (uint8_t*)u, *v8 = (uint8_t*)v;
 		__m128i x1_, v1_;
 		// find the boundaries
 		if (st < r - qlen + 1) st = r - qlen + 1;
@@ -152,14 +153,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 			}
 		} else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment
 			__m128i *pr = p + r * n_col_ - st_;
-			off[r] = st;
-			if (en0 < r && en0 < tlen - 1) { // to avoid backtracking out of the band; this assumes a fixed band
-				int8_t a, z = ((uint8_t*)s)[en0] + 2 * qe;
-				a = x8[en0-1] + v8[en0-1];
-				p_en0 = a > z? 1 : 0;
-				z = a > z? a : z;
-				p_en0 |= a - (z - q) > 0? 0x08 : 0;
-			}
+			off[r] = st, off_end[r] = en;
 			for (t = st_; t <= en_; ++t) {
 				__m128i d, z, a, b, xt1, vt1, ut, tmp;
 				__dp_code_block1;
@@ -185,14 +179,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 			}
 		} else { // gap right-alignment
 			__m128i *pr = p + r * n_col_ - st_;
-			off[r] = st;
-			if (en0 < r && en0 < tlen - 1) {
-				int8_t a, z = ((uint8_t*)s)[en0] + 2 * qe;
-				a = x8[en0-1] + v8[en0-1];
-				p_en0 = a >= z? 1 : 0;
-				z = a >= z? a : z;
-				p_en0 |= a - (z - q) >= 0? 0x08 : 0;
-			}
+			off[r] = st, off_end[r] = en;
 			for (t = st_; t <= en_; ++t) {
 				__m128i d, z, a, b, xt1, vt1, ut, tmp;
 				__dp_code_block1;
@@ -217,7 +204,6 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				_mm_store_si128(&pr[t], d);
 			}
 		}
-		if (with_cigar && en0 < r && en0 < tlen - 1) ((uint8_t*)(p + r * n_col_))[en0 - st] = p_en0;
 		if (!approx_max) { // find the exact max with a 32-bit score array
 			int32_t max_H, max_t;
 			// compute H[], max_H and max_t
@@ -289,9 +275,9 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	if (with_cigar) { // backtrack
 		int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR);
 		if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY))
-			ksw_backtrack(km, 1, rev_cigar, (uint8_t*)p, off, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+			ksw_backtrack(km, 1, rev_cigar, (uint8_t*)p, off, off_end, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
 		else if (ez->max_t >= 0 && ez->max_q >= 0)
-			ksw_backtrack(km, 1, rev_cigar, (uint8_t*)p, off, n_col_*16, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+			ksw_backtrack(km, 1, rev_cigar, (uint8_t*)p, off, off_end, n_col_*16, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
 		kfree(km, mem2); kfree(km, off);
 	}
 }

ksw2_ll_sse.c

@@ -0,0 +1,147 @@
+#include <stdlib.h>
+#include <stdint.h>
+#include <string.h>
+#include <emmintrin.h>
+#include "ksw2.h"
+
+#ifdef __GNUC__
+#define LIKELY(x) __builtin_expect((x),1)
+#define UNLIKELY(x) __builtin_expect((x),0)
+#else
+#define LIKELY(x) (x)
+#define UNLIKELY(x) (x)
+#endif
+
+typedef struct {
+	int qlen, slen;
+	uint8_t shift, mdiff, max, size;
+	__m128i *qp, *H0, *H1, *E, *Hmax;
+} kswq_t;
+
+/**
+ * Initialize the query data structure
+ *
+ * @param size   Number of bytes used to store a score; valid valures are 1 or 2
+ * @param qlen   Length of the query sequence
+ * @param query  Query sequence
+ * @param m      Size of the alphabet
+ * @param mat    Scoring matrix in a one-dimension array
+ *
+ * @return       Query data structure
+ */
+void *ksw_ll_qinit(void *km, int size, int qlen, const uint8_t *query, int m, const int8_t *mat)
+{
+	kswq_t *q;
+	int slen, a, tmp, p;
+
+	size = size > 1? 2 : 1;
+	p = 8 * (3 - size); // # values per __m128i
+	slen = (qlen + p - 1) / p; // segmented length
+	q = (kswq_t*)kmalloc(km, sizeof(kswq_t) + 256 + 16 * slen * (m + 4)); // a single block of memory
+	q->qp = (__m128i*)(((size_t)q + sizeof(kswq_t) + 15) >> 4 << 4); // align memory
+	q->H0 = q->qp + slen * m;
+	q->H1 = q->H0 + slen;
+	q->E  = q->H1 + slen;
+	q->Hmax = q->E + slen;
+	q->slen = slen; q->qlen = qlen; q->size = size;
+	// compute shift
+	tmp = m * m;
+	for (a = 0, q->shift = 127, q->mdiff = 0; a < tmp; ++a) { // find the minimum and maximum score
+		if (mat[a] < (int8_t)q->shift) q->shift = mat[a];
+		if (mat[a] > (int8_t)q->mdiff) q->mdiff = mat[a];
+	}
+	q->max = q->mdiff;
+	q->shift = 256 - q->shift; // NB: q->shift is uint8_t
+	q->mdiff += q->shift; // this is the difference between the min and max scores
+	// An example: p=8, qlen=19, slen=3 and segmentation:
+	//  {{0,3,6,9,12,15,18,-1},{1,4,7,10,13,16,-1,-1},{2,5,8,11,14,17,-1,-1}}
+	if (size == 1) {
+		int8_t *t = (int8_t*)q->qp;
+		for (a = 0; a < m; ++a) {
+			int i, k, nlen = slen * p;
+			const int8_t *ma = mat + a * m;
+			for (i = 0; i < slen; ++i)
+				for (k = i; k < nlen; k += slen) // p iterations
+					*t++ = (k >= qlen? 0 : ma[query[k]]) + q->shift;
+		}
+	} else {
+		int16_t *t = (int16_t*)q->qp;
+		for (a = 0; a < m; ++a) {
+			int i, k, nlen = slen * p;
+			const int8_t *ma = mat + a * m;
+			for (i = 0; i < slen; ++i)
+				for (k = i; k < nlen; k += slen) // p iterations
+					*t++ = (k >= qlen? 0 : ma[query[k]]);
+		}
+	}
+	return q;
+}
+
+int ksw_ll_i16(void *q_, int tlen, const uint8_t *target, int _gapo, int _gape, int *qe, int *te)
+{
+	kswq_t *q = (kswq_t*)q_;
+	int slen, i, gmax = 0, qlen8;
+	__m128i zero, gapoe, gape, *H0, *H1, *E, *Hmax;
+	uint16_t *H8;
+
+#define __max_8(ret, xx) do { \
+		(xx) = _mm_max_epi16((xx), _mm_srli_si128((xx), 8)); \
+		(xx) = _mm_max_epi16((xx), _mm_srli_si128((xx), 4)); \
+		(xx) = _mm_max_epi16((xx), _mm_srli_si128((xx), 2)); \
+		(ret) = _mm_extract_epi16((xx), 0); \
+	} while (0)
+
+	// initialization
+	*qe = *te = -1;
+	zero = _mm_set1_epi32(0);
+	gapoe = _mm_set1_epi16(_gapo + _gape);
+	gape = _mm_set1_epi16(_gape);
+	H0 = q->H0; H1 = q->H1; E = q->E; Hmax = q->Hmax;
+	slen = q->slen, qlen8 = slen * 8;
+	memset(E,    0, slen * sizeof(__m128i));
+	memset(H0,   0, slen * sizeof(__m128i));
+	memset(Hmax, 0, slen * sizeof(__m128i));
+	// the core loop
+	for (i = 0; i < tlen; ++i) {
+		int j, k, imax;
+		__m128i e, h, f = zero, max = zero, *S = q->qp + target[i] * slen; // s is the 1st score vector
+		h = _mm_load_si128(H0 + slen - 1); // h={2,5,8,11,14,17,-1,-1} in the above example
+		h = _mm_slli_si128(h, 2);
+		for (j = 0; LIKELY(j < slen); ++j) {
+			h = _mm_adds_epi16(h, *S++);
+			e = _mm_load_si128(E + j);
+			h = _mm_max_epi16(h, e);
+			h = _mm_max_epi16(h, f);
+			max = _mm_max_epi16(max, h);
+			_mm_store_si128(H1 + j, h);
+			h = _mm_subs_epu16(h, gapoe);
+			e = _mm_subs_epu16(e, gape);
+			e = _mm_max_epi16(e, h);
+			_mm_store_si128(E + j, e);
+			f = _mm_subs_epu16(f, gape);
+			f = _mm_max_epi16(f, h);
+			h = _mm_load_si128(H0 + j);
+		}
+		for (k = 0; LIKELY(k < 16); ++k) {
+			f = _mm_slli_si128(f, 2);
+			for (j = 0; LIKELY(j < slen); ++j) {
+				h = _mm_load_si128(H1 + j);
+				h = _mm_max_epi16(h, f);
+				_mm_store_si128(H1 + j, h);
+				h = _mm_subs_epu16(h, gapoe);
+				f = _mm_subs_epu16(f, gape);
+				if(UNLIKELY(!_mm_movemask_epi8(_mm_cmpgt_epi16(f, h)))) goto end_loop_i16;
+			}
+		}
+end_loop_i16:
+		__max_8(imax, max);
+		if (imax >= gmax) {
+			gmax = imax; *te = i;
+			memcpy(Hmax, H1, slen * sizeof(__m128i));
+		}
+		S = H1; H1 = H0; H0 = S;
+	}
+	for (i = 0, H8 = (uint16_t*)Hmax; i < qlen8; ++i)
+		if ((int)H8[i] == gmax) *qe = i / 8 + i % 8 * slen;
+	return gmax;
+}

main.c

@@ -8,7 +8,7 @@
 #include "minimap.h"
 #include "mmpriv.h"
 
-#define MM_VERSION "2.0-r191-dirty"
+#define MM_VERSION "2.0-r275"
 
 void liftrlimit()
 {
@@ -30,6 +30,7 @@ static struct option long_options[] = {
 	{ "print-seed",     no_argument,       0, 0 },
 	{ "max-chain-skip", required_argument, 0, 0 },
 	{ "min-dp-len",     required_argument, 0, 0 },
+	{ "print-aln-seq",  no_argument,       0, 0 },
 	{ "version",        no_argument,       0, 'V' },
 	{ "min-count",      required_argument, 0, 'n' },
 	{ "min-chain-score",required_argument, 0, 'm' },
@@ -85,6 +86,7 @@ int main(int argc, char *argv[])
 		else if (c == 0 && long_idx == 6) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED; // --print-seed
 		else if (c == 0 && long_idx == 7) opt.max_chain_skip = atoi(optarg); // --max-chain-skip
 		else if (c == 0 && long_idx == 8) opt.min_ksw_len = atoi(optarg); // --min-dp-len
+		else if (c == 0 && long_idx == 9) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ; // --print-aln-seq
 		else if (c == 'V') {
 			puts(MM_VERSION);
 			return 0;
@@ -192,7 +194,7 @@ int main(int argc, char *argv[])
 		mm_idx_t *mi = 0;
 		if (fpr) {
 			mi = mm_idx_load(fpr);
-			if (idx_par_set && mm_verbose >= 2 && (mi->k != k || mi->w != w || mi->is_hpc != mi->is_hpc))
+			if (idx_par_set && mm_verbose >= 2 && (mi->k != k || mi->w != w || mi->is_hpc != is_hpc))
 				fprintf(stderr, "[W::%s::%.3f*%.2f] Indexing parameters on the command line (-k/-w/-H) overridden by parameters in the prebuilt index.\n",
 						__func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0));
 		} else if (!mm_bseq_eof(fp)) {

map.c

@@ -31,7 +31,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 
 	opt->a = 2, opt->b = 4, opt->q = 4, opt->e = 2, opt->q2 = 24, opt->e2 = 1;
 	opt->zdrop = 400;
-	opt->min_dp_max = opt->min_chain_score;
+	opt->min_dp_max = opt->min_chain_score * opt->a;
 	opt->min_ksw_len = 200;
 }
 
@@ -209,8 +209,10 @@ mm_reg1_t *mm_map_frag(const mm_mapopt_t *opt, const mm_idx_t *mi, mm_tbuf_t *b,
 		mm128_t *p = &b->mini.a[i + m_st];
 		mm_match_t *q = &m[i];
 		const uint64_t *r = q->x.cr;
-		int k, q_span = p->x & 0xff;
+		int k, q_span = p->x & 0xff, is_tandem = 0;
 		if (q->n >= opt->mid_occ) continue;
+		if (i > 0 && p->x>>8 == b->mini.a[m_st + i - 1].x>>8) is_tandem = 1;
+		if (i < n - 1 && p->x>>8 == b->mini.a[m_st + i + 1].x>>8) is_tandem = 1;
 		for (k = 0; k < q->n; ++k) {
 			const char *tname = mi->seq[r[k]>>32].name;
 			int32_t rpos = (uint32_t)r[k] >> 1;
@@ -227,6 +229,7 @@ mm_reg1_t *mm_map_frag(const mm_mapopt_t *opt, const mm_idx_t *mi, mm_tbuf_t *b,
 				p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | (uint32_t)r[k]>>1;
 				p->y = (uint64_t)q_span << 32 | (qlen - ((q->qpos>>1) + 1 - q_span) - 1);
 			}
+			if (is_tandem) p->y |= MM_SEED_TANDEM;
 		}
 	}
 	n_a = j;
@@ -237,25 +240,33 @@ mm_reg1_t *mm_map_frag(const mm_mapopt_t *opt, const mm_idx_t *mi, mm_tbuf_t *b,
 
 	if (mm_dbg_flag & MM_DBG_PRINT_SEED)
 		for (i = 0; i < n_a; ++i)
-			fprintf(stderr, "SD\t%s\t%d\t%c\t%d\t%d\n", mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff));
+			fprintf(stderr, "SD\t%s\t%d\t%c\t%d\t%d\t%d\n", mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff),
+					i == 0? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));
 
 	n_u = mm_chain_dp(opt->max_gap, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, n_a, a, &u, b->km);
 	regs = mm_gen_regs(b->km, qlen, n_u, u, a);
 	*n_regs = n_u;
+
+	if (mm_dbg_flag & MM_DBG_PRINT_SEED)
+		for (j = 0; j < n_u; ++j)
+			for (i = regs[j].as; i < regs[j].as + regs[j].cnt; ++i)
+				fprintf(stderr, "CN\t%d\t%s\t%d\t%c\t%d\t%d\t%d\n", j, mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff),
+						i == regs[j].as? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));
+
 	if (!(opt->flag & MM_F_AVA)) { // don't choose primary mapping(s) for read overlap
 		mm_set_parent(b->km, opt->mask_level, *n_regs, regs);
-		mm_select_sub(b->km, opt->mask_level, opt->pri_ratio, opt->best_n, n_regs, regs);
+		mm_select_sub(b->km, opt->mask_level, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 		mm_join_long(b->km, opt, qlen, n_regs, regs, a); // TODO: this can be applied to all-vs-all in principle
 	}
 	if (opt->flag & MM_F_CIGAR) {
 		regs = mm_align_skeleton(b->km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs()
 		if (!(opt->flag & MM_F_AVA)) {
 			mm_set_parent(b->km, opt->mask_level, *n_regs, regs);
-			mm_select_sub(b->km, opt->mask_level, opt->pri_ratio, opt->best_n, n_regs, regs);
+			mm_select_sub(b->km, opt->mask_level, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 			mm_set_sam_pri(*n_regs, regs);
 		}
 	}
-	mm_set_mapq(*n_regs, regs);
+	mm_set_mapq(*n_regs, regs, opt->min_chain_score);
 
 	// free
 	kfree(b->km, a);
@@ -335,15 +346,15 @@ static void *worker_pipeline(void *shared, int step, void *in)
 			mm_bseq1_t *t = &s->seq[i];
 			for (j = 0; j < s->n_reg[i]; ++j) {
 				mm_reg1_t *r = &s->reg[i][j];
-				if (p->opt->flag & MM_F_OUT_SAM) mm_write_sam(&p->str, mi, t, r);
+				if (p->opt->flag & MM_F_OUT_SAM) mm_write_sam(&p->str, mi, t, r, s->n_reg[i], s->reg[i]);
 				else mm_write_paf(&p->str, mi, t, r);
 				puts(p->str.s);
-				free(r->p);
 			}
 			if (s->n_reg[i] == 0 && (p->opt->flag & MM_F_OUT_SAM)) {
-				mm_write_sam(&p->str, 0, t, 0);
+				mm_write_sam(&p->str, 0, t, 0, 0, 0);
 				puts(p->str.s);
 			}
+			for (j = 0; j < s->n_reg[i]; ++j) free(s->reg[i][j].p);
 			free(s->reg[i]);
 			free(s->seq[i].seq); free(s->seq[i].name);
 			if (s->seq[i].qual) free(s->seq[i].qual);

minimap.h

@@ -46,6 +46,7 @@ typedef struct {
 	mm_idx_seq_t *seq;  // sequence name, length and offset
 	uint32_t *S;        // 4-bit packed sequence
 	mm_idx_bucket_t *B; // index
+	void *km;
 } mm_idx_t;
 
 typedef struct {
@@ -60,7 +61,7 @@ typedef struct {
 typedef struct {
 	int32_t id;
 	uint32_t cnt:31, rev:1;
-	uint32_t rid:31, rep:1;
+	uint32_t rid:31, inv:1;
 	int32_t score;
 	int32_t qs, qe, rs, re;
 	int32_t parent, subsc;

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "19 July 2017" "minimap2-2.0-r190-dirty" "Bioinformatics tools"
+.TH minimap2 1 "8 August 2017" "minimap2-2.0-r275" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -137,10 +137,8 @@ number of minimizers [3]
 .BI -m \ INT
 Discard chains with chaining score
 .RI < INT
-[40]. Chaining score equals the approximate number of matching bases (exact if
-not using
-.BR -H )
-minus base-2 logarithm gap penalty. It is computed with dynamic programming.
+[40]. Chaining score equals the approximate number of matching bases minus a
+concave gap penalty. It is computed with dynamic programming.
 .TP
 .B -X
 Perform all-vs-all mapping. In this mode, if the query sequence name is
@@ -249,35 +247,49 @@ are:
 .RS
 .TP 8
 .B map-pb
-PacBio/Oxford Nanopore read to reference mapping (-Hk19)
+PacBio/Oxford Nanopore read to reference mapping
+.RB ( -Hk19 )
 .TP
 .B map10k
 The same as
 .B map-pb
-(-Hk19)
+.RB ( -Hk19 )
 .TP
 .B map-ont
-Slightly more sensitive for Oxford Nanopore to reference mapping (-k15). For
-PacBio reads, HPC minimizers consistently leads to faster performance and more
-sensitive results in comparison to normal minimizers. For Oxford Nanopore data,
-normal minimizers are better, though not much. The effectiveness of HPC is
-determined by the sequencing error mode.
+Slightly more sensitive for Oxford Nanopore to reference mapping
+.RB ( -k15 ).
+For PacBio reads, HPC minimizers consistently leads to faster performance and
+more sensitive results in comparison to normal minimizers. For Oxford Nanopore
+data, normal minimizers are better, though not much. The effectiveness of HPC
+is determined by the sequencing error mode.
 .TP
 .B asm5
-Long assembly to reference mapping (-k19 -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200).
+Long assembly to reference mapping
+.RB ( -k19
+.B -w19 -A1 -B19 -O39,81 -E3,1 -s200
+.BR -z200 ).
 Typically, the alignment will not extend to regions with 5% or higher sequence
 divergence. Only use this preset if the average divergence is far below 5%.
 .TP
 .B asm10
-Long assembly to reference mapping (-k19 -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200). Up
-to 10% sequence divergence.
+Long assembly to reference mapping
+.RB ( -k19
+.B -w19 -A1 -B9 -O16,41 -E2,1 -s200
+.BR -z200 ).
+Up to 10% sequence divergence.
 .TP 8
 .B ava-pb
-PacBio all-vs-all overlap mapping (-Hk19 -w5 -Xp0 -m100 -K500m -g10000 --max-chain-skip 25)
+PacBio all-vs-all overlap mapping
+.RB ( -Hk19
+.B -w5 -Xp0 -m100 -K500m -g10000 --max-chain-skip
+.BR 25 ).
 .TP 8
 .B ava-ont
-Oxford Nanopore all-vs-all overlap mapping (-k15 -w5 -Xp0 -m100 -K500m -g10000
---max-chain-skip 25). Similarly, the major difference from
+Oxford Nanopore all-vs-all overlap mapping
+.RB ( -k15
+.B -w5 -Xp0 -m100 -K500m -g10000 --max-chain-skip
+.BR 25 ).
+Similarly, the major difference from
 .B ava-pb
 is that this preset is not using HPC minimizers.
 .RE
@@ -331,6 +343,7 @@ cb | cb | cb
 r | c | l .
 Tag	Type	Description
 _
+tp	A	Type of aln: P/primary, S/secondary and I/inversion
 cm	i	Number of minimizers on the chain
 s1	i	Chaining score
 s2	i	Chaining score of the best secondary chain
@@ -344,13 +357,9 @@ cg	Z	CIGAR string (only in PAF)
 .SH LIMITATIONS
 .TP 2
 *
-At the alignment phase, minimap2 performs global alignments between minimizer
-hits. If the positions of these minimizer hits are incorrect, the final
-alignment may be suboptimal or unnecessarily fragmented.
-.TP
-*
-Minimap2 may produce poor alignments that may need post-filtering. We are still
-exploring a reliable and consistent way to report good alignments.
+Minimap2 may produce suboptimal alignments through long low-complexity regions
+where seed positions may be suboptimal. This should not be a big concern
+because even the optimal alignment may be wrong in such regions.
 .TP
 *
 Minimap2 does not work well with Illumina short reads as of now.

misc/mapstat.js

@@ -144,7 +144,7 @@ while (file.readline(buf) >= 0) {
 		}
 		if (n_cigar > 65535) ++n_cigar_64k;
 		if (ql + sclip != aqlen)
-			warn("WARNING: aligned query length is inconsistent with CIGAR at line " + lineno);
+			warn("WARNING: aligned query length is inconsistent with CIGAR at line " + lineno + " (" + (ql+sclip) + " != " + aqlen + ")");
 		if (atlen != null && atlen != tl)
 			warn("WARNING: aligned reference length is inconsistent with CIGAR at line " + lineno);
 		if (is_sam) {

misc/pbsim-eval.js

@@ -0,0 +1,143 @@
+var getopt = function(args, ostr) {
+	var oli; // option letter list index
+	if (typeof(getopt.place) == 'undefined')
+		getopt.ind = 0, getopt.arg = null, getopt.place = -1;
+	if (getopt.place == -1) { // update scanning pointer
+		if (getopt.ind >= args.length || args[getopt.ind].charAt(getopt.place = 0) != '-') {
+			getopt.place = -1;
+			return null;
+		}
+		if (getopt.place + 1 < args[getopt.ind].length && args[getopt.ind].charAt(++getopt.place) == '-') { // found "--"
+			++getopt.ind;
+			getopt.place = -1;
+			return null;
+		}
+	}
+	var optopt = args[getopt.ind].charAt(getopt.place++); // character checked for validity
+	if (optopt == ':' || (oli = ostr.indexOf(optopt)) < 0) {
+		if (optopt == '-') return null; //  if the user didn't specify '-' as an option, assume it means null.
+		if (getopt.place < 0) ++getopt.ind;
+		return '?';
+	}
+	if (oli+1 >= ostr.length || ostr.charAt(++oli) != ':') { // don't need argument
+		getopt.arg = null;
+		if (getopt.place < 0 || getopt.place >= args[getopt.ind].length) ++getopt.ind, getopt.place = -1;
+	} else { // need an argument
+		if (getopt.place >= 0 && getopt.place < args[getopt.ind].length)
+			getopt.arg = args[getopt.ind].substr(getopt.place);
+		else if (args.length <= ++getopt.ind) { // no arg
+			getopt.place = -1;
+			if (ostr.length > 0 && ostr.charAt(0) == ':') return ':';
+			return '?';
+		} else getopt.arg = args[getopt.ind]; // white space
+		getopt.place = -1;
+		++getopt.ind;
+	}
+	return optopt;
+}
+
+var c, max_mapq = 60, mode = 0, err_out_q = 256, print_err = false, ovlp_ratio = 0.333;
+while ((c = getopt(arguments, "Q:r:m:")) != null) {
+	if (c == 'Q') err_out_q = parseInt(getopt.arg), print_err = true;
+	else if (c == 'r') ovlp_ratio = parseFloat(getopt.arg);
+	else if (c == 'm') mode = parseInt(getopt.arg);
+}
+
+var file = arguments.length == getopt.ind? new File() : new File(arguments[getopt.ind]);
+var buf = new Bytes();
+
+var tot = [], err = [];
+for (var q = 0; q <= max_mapq; ++q)
+	tot[q] = err[q] = 0;
+
+function is_correct(s, b)
+{
+	if (s[0] != b[0] || s[3] != b[3]) return false;
+	var o, l;
+	if (s[1] < b[1]) {
+		if (s[2] <= b[1]) return false;
+		o = (s[2] < b[2]? s[2] : b[2]) - b[1];
+		l = (s[2] > b[2]? s[2] : b[2]) - s[1];
+	} else {
+		if (b[2] <= s[1]) return false;
+		o = (s[2] < b[2]? s[2] : b[2]) - s[1];
+		l = (s[2] > b[2]? s[2] : b[2]) - b[1];
+	}
+	return o/l > ovlp_ratio? true : false;
+}
+
+function count_err(qname, a, tot, err, mode)
+{
+	var s = qname.split("!");
+	if (a.length == 0) return;
+	if (s.length < 5 || (s[4] != '+' && s[4] != '-'))
+		throw Error("Failed to parse pbsim2fa read names '" + qname + "'");
+	s[2] = parseInt(s[2]);
+	s[3] = parseInt(s[3]);
+	s.shift(); // skip pbsim orginal read name
+	if (mode == 0 || mode == 1) { // longest only or first only
+		var max_i = 0;
+		if (mode == 0) {
+			var max = 0;
+			for (var i = 0; i < a.length; ++i)
+				if (a[i][2] - a[i][1] > max)
+					max = a[i][2] - a[i][1], max_i = i;
+		}
+		var mapq = a[max_i][4];
+		++tot[mapq];
+		if (!is_correct(s, a[max_i])) {
+			if (mapq >= err_out_q)
+				print('E', qname, a[max_i].join("\t"));
+			++err[mapq];
+		}
+	} else if (mode == 2) { // all primary mode
+		var max_err_mapq = -1, max_mapq = 0, max_err_i = -1;
+		for (var i = 0; i < a.length; ++i) {
+			max_mapq = max_mapq > a[i][4]? max_mapq : a[i][4];
+			if (!is_correct(s, a[i]))
+				if (a[i][4] > max_err_mapq)
+					max_err_mapq = a[i][4], max_err_i = i;
+		}
+		if (max_err_mapq >= 0) {
+			++tot[max_err_mapq], ++err[max_err_mapq];
+			if (max_err_mapq >= err_out_q)
+				print('E', qname, a[max_err_i].join("\t"));
+		} else ++tot[max_mapq];
+	}
+}
+
+var lineno = 0, last = null, a = [];
+while (file.readline(buf) >= 0) {
+	var line = buf.toString();
+	++lineno;
+	if (line[0] != '@') {
+		var t = line.split("\t");
+		if (last != t[0]) {
+			if (last != null) count_err(last, a, tot, err, mode);
+			a = [], last = t[0];
+		}
+		if (t[4] == '+' || t[4] == '-') { // PAF
+			if (/\ts1:i:\d+/.test(line) && !/\ts2:i:\d+/.test(line)) // secondary alignment in minimap2 PAF
+				continue;
+			var mapq = parseInt(t[11]);
+			if (mapq > max_mapq) mapq = max_mapq;
+			a.push([t[5], parseInt(t[7]), parseInt(t[8]), t[4], mapq]);
+		}
+	}
+}
+if (last != null) count_err(last, a, tot, err, mode);
+
+buf.destroy();
+file.close();
+
+var sum_tot = 0, sum_err = 0, q_out = -1, sum_tot2 = 0, sum_err2 = 0;
+for (var q = max_mapq; q >= 0; --q) {
+	if (tot[q] == 0) continue;
+	if (q_out < 0 || err[q] > 0) {
+		if (q_out >= 0) print('Q', q_out, sum_tot, sum_err, (sum_err2/sum_tot2).toFixed(9));
+		sum_tot = sum_err = 0, q_out = q;
+	}
+	sum_tot += tot[q], sum_err += err[q];
+	sum_tot2 += tot[q], sum_err2 += err[q];
+}
+print('Q', q_out, sum_tot, sum_err, (sum_err2/sum_tot2).toFixed(9));

misc/pbsim2fa.js

@@ -0,0 +1,81 @@
+Bytes.prototype.reverse = function()
+{
+	for (var i = 0; i < this.length>>1; ++i) {
+		var tmp = this[i];
+		this[i] = this[this.length - i - 1];
+		this[this.length - i - 1] = tmp;
+	}
+}
+
+// reverse complement a DNA string
+Bytes.prototype.revcomp = function()
+{
+	if (Bytes.rctab == null) {
+		var s1 = 'WSATUGCYRKMBDHVNwsatugcyrkmbdhvn';
+		var s2 = 'WSTAACGRYMKVHDBNwstaacgrymkvhdbn';
+		Bytes.rctab = [];
+		for (var i = 0; i < 256; ++i) Bytes.rctab[i] = 0;
+		for (var i = 0; i < s1.length; ++i)
+			Bytes.rctab[s1.charCodeAt(i)] = s2.charCodeAt(i);
+	}
+	for (var i = 0; i < this.length>>1; ++i) {
+		var tmp = this[this.length - i - 1];
+		this[this.length - i - 1] = Bytes.rctab[this[i]];
+		this[i] = Bytes.rctab[tmp];
+	}
+	if (this.length&1)
+		this[this.length>>1] = Bytes.rctab[this[this.length>>1]];
+}
+
+if (arguments.length < 2) {
+	print("Usage: k8 pbsim2paf.js <chr.list> <pbsim1.maf> [[pbsim2.maf] ...]");
+	exit(1);
+}
+
+var file, buf = new Bytes(), buf2 = new Bytes();
+file = new File(arguments[0]);
+var chr_list = [];
+while (file.readline(buf) >= 0) {
+	var t = buf.toString().split(/\s+/);
+	chr_list.push(t[0]);
+}
+file.close();
+
+for (var k = 1; k < arguments.length; ++k) {
+	var fn = arguments[k];
+	file = new File(fn);
+	var state = 0, reg;
+	while (file.readline(buf) >= 0) {
+		var line = buf.toString();
+		if (state == 0 && line.charAt(0) == 'a') {
+			state = 1;
+		} else if (state == 1 && line.charAt(0) == 's') {
+			var t = line.split(/\s+/);
+			var st = parseInt(t[2]);
+			reg = [st, st + parseInt(t[3])];
+			state = 2;
+		} else if (state == 2 && line.charAt(0) == 's') {
+			var m, t = line.split(/\s+/);
+			if ((m = /S(\d+)_\d+/.exec(t[1])) == null) throw Error("Failed to parse the read name");
+			var chr_id = parseInt(m[1]) - 1;
+			if (chr_id >= chr_list.length) throw Error("Index outside the chr list");
+			var name = [t[1], chr_list[chr_id], reg[0], reg[1], t[4]].join("!");
+			var seq = t[6].replace(/\-/g, "");
+			if (seq.length != parseInt(t[5])) throw Error("Inconsistent read length");
+			if (seq.indexOf("NN") < 0) {
+				if (t[4] == '-') {
+					buf2.set(seq, 0);
+					buf2.length = seq.length;
+					buf2.revcomp();
+					seq = buf2.toString();
+				}
+				print(">" + name);
+				print(seq);
+			}
+			state = 0;
+		}
+	}
+	file.close();
+}
+buf.destroy();
+buf2.destroy();

mmpriv.h

@@ -8,9 +8,14 @@
 #define MM_PARENT_UNSET   (-1)
 #define MM_PARENT_TMP_PRI (-2)
 
-#define MM_DBG_NO_KALLOC   0x1
-#define MM_DBG_PRINT_QNAME 0x2
-#define MM_DBG_PRINT_SEED  0x4
+#define MM_DBG_NO_KALLOC     0x1
+#define MM_DBG_PRINT_QNAME   0x2
+#define MM_DBG_PRINT_SEED    0x4
+#define MM_DBG_PRINT_ALN_SEQ 0x8
+
+#define MM_SEED_LONG_JOIN  (1ULL<<40)
+#define MM_SEED_IGNORE     (1ULL<<41)
+#define MM_SEED_TANDEM     (1ULL<<42)
 
 #ifndef kroundup32
 #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
@@ -36,7 +41,7 @@ void radix_sort_64(uint64_t *beg, uint64_t *end);
 uint32_t ks_ksmall_uint32_t(size_t n, uint32_t arr[], size_t kk);
 
 void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r);
-void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r);
+void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int n_regs, const mm_reg1_t *regs);
 int mm_chain_dp(int max_dist, int bw, int max_skip, int min_cnt, int min_sc, int64_t n, mm128_t *a, uint64_t **_u, void *km);
 mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
 
@@ -45,11 +50,11 @@ void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a);
 void mm_sync_regs(void *km, int n_regs, mm_reg1_t *regs);
 int mm_set_sam_pri(int n, mm_reg1_t *r);
 void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r);
-void mm_select_sub(void *km, float mask_level, float pri_ratio, int best_n, int *n_, mm_reg1_t *r);
+void mm_select_sub(void *km, float mask_level, float pri_ratio, int min_diff, int best_n, int *n_, mm_reg1_t *r);
 void mm_filter_regs(void *km, const mm_mapopt_t *opt, int *n_regs, mm_reg1_t *regs);
 void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs, mm128_t *a);
 void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r);
-void mm_set_mapq(int n_regs, mm_reg1_t *regs);
+void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc);
 
 #ifdef __cplusplus
 }