Release minimap2-2.8 (r672)

r671: cleanup command line options
r670: added PE support to mappy
2026-09-26 17:18:11 +08:00 · 2018-02-01 12:50:20 -05:00 · 2018-01-31 13:59:52 -05:00 · 2018-01-31 11:33:08 -05:00 · 2018-01-30 20:05:02 -05:00 · 2018-01-30 10:11:50 -05:00
40 changed files with 3721 additions and 454 deletions
@@ -10,10 +10,6 @@ matrix:
      python: "2.7"
      before_install: pip install cython
      script: python setup.py build_ext
-    - language: python
-      python: "3.3"
-      before_install: pip install cython
-      script: python setup.py build_ext
    - language: python
      python: "3.5"
      before_install: pip install cython
@@ -1,17 +1,24 @@
 CFLAGS=		-g -Wall -O2 -Wc++-compat
 CPPFLAGS=	-DHAVE_KALLOC
 INCLUDES=
-OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o index.o chain.o align.o hit.o map.o format.o pe.o ksw2_ll_sse.o
+OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o chain.o align.o hit.o map.o format.o pe.o esterr.o ksw2_ll_sse.o
 PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite
 LIBS=		-lm -lz -lpthread

+ifeq ($(arm_neon),)
 ifeq ($(sse2only),)
 	OBJS+=ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o
 else
 	OBJS+=ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_exts2_sse.o
 endif
+else
+    OBJS+=ksw2_extz2_neon.o ksw2_extd2_neon.o ksw2_exts2_neon.o
+    CFLAGS+=-D_FILE_OFFSET_BITS=64 -mfpu=neon -fsigned-char
+    INCLUDES+=-I sse2neon
+endif

+.PHONY:all extra clean depend
 .SUFFIXES:.c .o

 .c.o:
@@ -33,6 +40,8 @@ libminimap2.a:$(OBJS)
 sdust:sdust.c getopt.o kalloc.o kalloc.h kdq.h kvec.h kseq.h sdust.h
 		$(CC) -D_SDUST_MAIN $(CFLAGS) $< getopt.o kalloc.o -o $@ -lz

+# SSE-specific targets on x86/x86_64
+
 ksw2_extz2_sse41.o:ksw2_extz2_sse.c ksw2.h kalloc.h
 		$(CC) -c -msse4 $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@

@@ -54,8 +63,21 @@ ksw2_exts2_sse2.o:ksw2_exts2_sse.c ksw2.h kalloc.h
 ksw2_dispatch.o:ksw2_dispatch.c ksw2.h
 		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@

+# NEON-specific targets on ARM
+
+ksw2_extz2_neon.o:ksw2_extz2_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_SSE2_ONLY -D__SSE2__ $(INCLUDES) $< -o $@
+
+ksw2_extd2_neon.o:ksw2_extd2_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_SSE2_ONLY -D__SSE2__ $(INCLUDES) $< -o $@
+
+ksw2_exts2_neon.o:ksw2_exts2_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_SSE2_ONLY -D__SSE2__ $(INCLUDES) $< -o $@
+
+# other non-file targets
+
 clean:
-		rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist mappy.so mappy.c python/mappy.c mappy.egg*
+		rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist mappy*.so mappy.c python/mappy.c mappy.egg*

 depend:
 		(LC_ALL=C; export LC_ALL; makedepend -Y -- $(CFLAGS) $(CPPFLAGS) -- *.c)
@@ -65,6 +87,7 @@ depend:
 align.o: minimap.h mmpriv.h bseq.h ksw2.h kalloc.h
 bseq.o: bseq.h kvec.h kalloc.h kseq.h
 chain.o: minimap.h mmpriv.h bseq.h kalloc.h
+esterr.o: mmpriv.h minimap.h bseq.h
 example.o: minimap.h kseq.h
 format.o: kalloc.h mmpriv.h minimap.h bseq.h
 getopt.o: getopt.h
@@ -75,9 +98,12 @@ ksw2_extd2_sse.o: ksw2.h kalloc.h
 ksw2_exts2_sse.o: ksw2.h kalloc.h
 ksw2_extz2_sse.o: ksw2.h kalloc.h
 ksw2_ll_sse.o: ksw2.h kalloc.h
+kthread.o: kthread.h
 main.o: bseq.h minimap.h mmpriv.h getopt.h
 map.o: kthread.h kvec.h kalloc.h sdust.h mmpriv.h minimap.h bseq.h khash.h
-misc.o: minimap.h ksort.h
+map.o: ksort.h
+misc.o: mmpriv.h minimap.h bseq.h ksort.h
+options.o: mmpriv.h minimap.h bseq.h
 pe.o: mmpriv.h minimap.h bseq.h kvec.h kalloc.h ksort.h
 sdust.o: kalloc.h kdq.h kvec.h sdust.h
-sketch.o: kvec.h kalloc.h minimap.h
+sketch.o: kvec.h kalloc.h mmpriv.h minimap.h bseq.h
@@ -1,3 +1,89 @@
+Release 2.8-r672 (1 February 2018)
+----------------------------------
+
+Notable changes in this release include:
+
+ * Speed up short-read alignment by ~10%. The overall mapping accuracy stays
+   the same, but the output alignments are not always identical to v2.7 due to
+   unstable sorting employed during chaining. Long-read alignment is not
+   affected by this change as the speedup is short-read specific.
+
+ * Mappy now supports paired-end short-read alignment (#87). Please see
+   python/README.rst for details.
+
+ * Added option --for-only and --rev-only to perform alignment against the
+   forward or the reverse strand of the reference genome only (#91).
+
+ * Alleviated the issue with undesired diagonal alignment in the self mapping
+   mode (#10). Even if the output is not ideal, it should not interfere with
+   other alignments. Fully resolving the issue is intricate and may require
+   additional heuristic thresholds.
+
+ * Enhanced error checking against incorrect input (#92 and #96).
+
+For long query sequences, minimap2 should output identical alignments to v2.7.
+
+(2.8: 1 February 2018, r672)
+
+
+
+Release 2.7-r654 (9 January 2018)
+---------------------------------
+
+This release fixed a bug in the splice mode and added a few minor features:
+
+ * Fixed a bug that occasionally takes an intron as a long deletion in the
+   splice mode. This was caused by wrong backtracking at the last CIGAR
+   operator. The current fix eliminates the error, but it is not optimal in
+   that it often produces a wrong junction when the last operator is an intron.
+   A future version of minimap2 may improve upon this.
+
+ * Support high-end ARM CPUs that implement the NEON instruction set (#81).
+   This enables minimap2 to work on Raspberry Pi 3 and Odroid XU4.
+
+ * Added a C API to construct a minimizer index from a set of C strings (#80).
+
+ * Check scoring specified on the command line (#79). Due to the 8-bit limit,
+   excessively large score penalties fail minimap2.
+
+For genomic sequences, minimap2 should give identical alignments to v2.6.
+
+(2.7: 9 January 2018, r654)
+
+
+
+Release 2.6-r623 (12 December 2017)
+-----------------------------------
+
+This release adds several features and fixes two minor bugs:
+
+ * Optionally build an index without sequences. This helps to reduce the
+   peak memory for read overlapping and is automatically applied when
+   base-level alignment is not requested.
+
+ * Approximately estimate per-base sequence divergence (i.e. 1-identity)
+   without performing base-level alignment, using a MashMap-like method. The
+   estimate is written to a new dv:f tag.
+
+ * Reduced the number of tiny terminal exons in RNA-seq alignment. The current
+   setting is conservative. Increase --end-seed-pen to drop more such exons.
+
+ * Reduced the peak memory when aligning long query sequences.
+
+ * Fixed a bug that is caused by HPC minimizers longer than 256bp. This should
+   have no effect in practice, but it is recommended to rebuild HPC indices if
+   possible.
+
+ * Fixed a bug when identifying identical hits (#71). This should only affect
+   artifactual reference consisting of near identical sequences.
+
+For genomic sequences, minimap2 should give nearly identical alignments to
+v2.5, except the new dv:f tag.
+
+(2.6: 12 December 2017, r623)
+
+
+
 Release 2.5-r572 (11 November 2017)
 -----------------------------------

@@ -1,10 +1,7 @@
-[![Release](https://img.shields.io/badge/Release-v2.5-blue.svg?style=flat)](https://github.com/lh3/minimap2/releases)
-[![BioConda](https://img.shields.io/conda/vn/bioconda/minimap2.svg?style=flat)](https://anaconda.org/bioconda/minimap2)
+[![GitHub Downloads](https://img.shields.io/github/downloads/lh3/minimap2/total.svg?style=social&logo=github&label=Download)](https://github.com/lh3/minimap2/releases)
+[![BioConda Install](https://img.shields.io/conda/dn/bioconda/minimap2.svg?style=flag&label=BioConda%20install)](https://anaconda.org/bioconda/minimap2)
 [![PyPI](https://img.shields.io/pypi/v/mappy.svg?style=flat)](https://pypi.python.org/pypi/mappy)
-[![Python Version](https://img.shields.io/pypi/pyversions/mappy.svg?style=flat)](https://pypi.python.org/pypi/mappy)
-[![License](https://img.shields.io/badge/License-MIT-blue.svg?style=flat)](LICENSE.txt)
 [![Build Status](https://travis-ci.org/lh3/minimap2.svg?branch=master)](https://travis-ci.org/lh3/minimap2)
-[![Downloads](https://img.shields.io/github/downloads/lh3/minimap2/total.svg?style=flat)](https://github.com/lh3/minimap2/releases)
 ## <a name="started"></a>Getting Started
 ```sh
 git clone https://github.com/lh3/minimap2
@@ -68,18 +65,21 @@ Detailed evaluations are available from the [minimap2 preprint][preprint].

 ### <a name="install"></a>Installation

-Minimap2 only works on x86-64 CPUs. You can acquire precompiled binaries from
+Minimap2 is optimized for x86-64 CPUs. You can acquire precompiled binaries from
 the [release page][release] with:
 ```sh
-curl -L https://github.com/lh3/minimap2/releases/download/v2.5/minimap2-2.5_x64-linux.tar.bz2 \
+curl -L https://github.com/lh3/minimap2/releases/download/v2.8/minimap2-2.8_x64-linux.tar.bz2 \
  | tar -jxvf -
-./minimap2-2.5_x64-linux/minimap2
+./minimap2-2.8_x64-linux/minimap2
 ```
 If you want to compile from the source, you need to have a C compiler, GNU make
 and zlib development files installed. Then type `make` in the source code
 directory to compile. If you see compilation errors, try `make sse2only=1`
 to disable SSE4 code, which will make minimap2 slightly slower.

+Minimap2 also works with ARM CPUs supporting the NEON instruction sets. To
+compile, use `make arm_neon=1`.
+
 ### <a name="general"></a>General usage

 Without any options, minimap2 takes a reference database and a query sequence
@@ -154,6 +154,20 @@ slight improvement to accuracy. For noisy Nanopore Direct RNA-seq reads, it is
 recommended to use a smaller k-mer size for increased sensitivity to the first
 or the last exons.

+Minimap2 rates an alignment by the score of the max-scoring sub-segment,
+*excluding* introns, and marks the best alignment as primary in SAM. When a
+spliced gene also has unspliced pseudogenes, minimap2 does not intentionally
+prefer spliced alignment, though in practice it more often marks the spliced
+alignment as the primary. By default, minimap2 outputs up to five secondary
+alignments (i.e. likely pseudogenes in the context of RNA-seq mapping). This
+can be tuned with option **-N**.
+
+For long RNA-seq reads, minimap2 may produce chimeric alignments potentially
+caused by gene fusions/structural variations or by an intron longer than the
+max intron length **-G** (200k by default). For now, it is not recommended to
+apply an excessively large **-G** as this slows down minimap2 and sometimes
+leads to false alignments.
+
 It is worth noting that by default `-x splice` prefers GT[A/G]..[C/T]AG
 over GT[C/T]..[A/G]AG, and then over other splicing signals. Considering
 one additional base improves the junction accuracy for noisy reads, but
@@ -309,7 +323,7 @@ highlighted in bold. The description may help to tune minimap2 parameters.

 ### <a name="help"></a>Getting help

-Manpage [minimap2.1](minimap2.1) provides detailed description of minimap2
+Manpage [minimap2.1][manpage] provides detailed description of minimap2
 command line options and optional tags. If you encounter bugs or have further
 questions or requests, you can raise an issue at the [issue page][issue].
 There is not a specific mailing list for the time being.
@@ -341,12 +355,9 @@ mappy` or [from BioConda][mappyconda] via `conda install -c bioconda mappy`.
  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 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.
+* Minimap2 requires SSE2 instructions on x86 CPUs or NEON on ARM CPUs. It is
+  possible to add non-SIMD support, but it would make minimap2 slower by
+  several times.



@@ -363,3 +374,4 @@ warmly welcomed.
 [mappyconda]: https://anaconda.org/bioconda/mappy
 [issue]: https://github.com/lh3/minimap2/issues
 [k8]: https://github.com/attractivechaos/k8
+[manpage]: https://lh3.github.io/minimap2/minimap2.html
@@ -1,6 +1,7 @@
 #include <assert.h>
 #include <string.h>
 #include <stdlib.h>
+#include <math.h>
 #include "minimap.h"
 #include "mmpriv.h"
 #include "ksw2.h"
@@ -120,7 +121,7 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 	}
 }

-static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *qual, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e)
+static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e)
 {
 	uint32_t k, l, toff = 0, qoff = 0;
 	int32_t s = 0, max = 0, qshift, tshift;
@@ -209,6 +210,13 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 		ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, flag, ez);
 	else
 		ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez);
+	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
+		int i;
+		fprintf(stderr, "score=%d, cigar=", ez->score);
+		for (i = 0; i < ez->n_cigar; ++i)
+			fprintf(stderr, "%d%c", ez->cigar[i]>>4, "MIDN"[ez->cigar[i]&0xf]);
+		fprintf(stderr, "\n");
+	}
 }

 static inline int mm_get_hplen_back(const mm_idx_t *mi, uint32_t rid, uint32_t x)
@@ -222,7 +230,7 @@ static inline int mm_get_hplen_back(const mm_idx_t *mi, uint32_t rid, uint32_t x

 static inline void mm_adjust_minier(const mm_idx_t *mi, uint8_t *const qseq0[2], mm128_t *a, int32_t *r, int32_t *q)
 {
-	if (mi->is_hpc) {
+	if (mi->flag & MM_I_HPC) {
 		const uint8_t *qseq = qseq0[a->x>>63];
 		int i, c;
 		*q = (int32_t)a->y;
@@ -341,7 +349,49 @@ static void mm_max_stretch(const mm_mapopt_t *opt, const mm_reg1_t *r, const mm1
 	*as = max_i, *cnt = max_len;
 }

-static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], uint8_t *qual0[2], mm_reg1_t *r, mm_reg1_t *r2, int n_a, mm128_t *a, ksw_extz_t *ez, int splice_flag)
+static int mm_seed_ext_score(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, const int8_t mat[25], int qlen, uint8_t *qseq0[2], const mm128_t *a)
+{
+	uint8_t *qseq, *tseq;
+	int q_span = a->y>>32&0xff, qs, qe, rs, re, rid, score, q_off, t_off, ext_len = opt->anchor_ext_len;
+	void *qp;
+	rid = a->x<<1>>33;
+	re = (uint32_t)a->x + 1, rs = re - q_span;
+	qe = (uint32_t)a->y + 1, qs = qe - q_span;
+	rs = rs - ext_len > 0? rs - ext_len : 0;
+	qs = qs - ext_len > 0? qs - ext_len : 0;
+	re = re + ext_len < mi->seq[rid].len? re + ext_len : mi->seq[rid].len;
+	qe = qe + ext_len < qlen? qe + ext_len : qlen;
+	tseq = (uint8_t*)kmalloc(km, re - rs);
+	mm_idx_getseq(mi, rid, rs, re, tseq);
+	qseq = qseq0[a->x>>63] + qs;
+	qp = ksw_ll_qinit(km, 2, qe - qs, qseq, 5, mat);
+	score = ksw_ll_i16(qp, re - rs, tseq, opt->q, opt->e, &q_off, &t_off);
+	kfree(km, tseq);
+	kfree(km, qp);
+	return score;
+}
+
+static void mm_fix_bad_ends_splice(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, const mm_reg1_t *r, const int8_t mat[25], int qlen, uint8_t *qseq0[2], const mm128_t *a, int *as1, int *cnt1)
+{ // this assumes a very crude k-mer based mode; it is not necessary to use a good model just for filtering bounary exons
+	int score;
+	double log_gap;
+	*as1 = r->as, *cnt1 = r->cnt;
+	if (r->cnt < 3) return;
+	log_gap = log((int32_t)a[r->as + 1].x - (int32_t)a[r->as].x);
+	if ((a[r->as].y>>32&0xff) < log_gap + opt->anchor_ext_shift) {
+		score = mm_seed_ext_score(km, opt, mi, mat, qlen, qseq0, &a[r->as]);
+		if ((double)score / mat[0] < log_gap + opt->anchor_ext_shift) // a more exact format is "score < log_4(gap) + shift"
+			++(*as1), --(*cnt1);
+	}
+	log_gap = log((int32_t)a[r->as + r->cnt - 1].x - (int32_t)a[r->as + r->cnt - 2].x);
+	if ((a[r->as + r->cnt - 1].y>>32&0xff) < log_gap + opt->anchor_ext_shift) {
+		score = mm_seed_ext_score(km, opt, mi, mat, qlen, qseq0, &a[r->as + r->cnt - 1]);
+		if ((double)score / mat[0] < log_gap + opt->anchor_ext_shift)
+			--(*cnt1);
+	}
+}
+
+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, int n_a, mm128_t *a, ksw_extz_t *ez, int splice_flag)
 {
 	int is_sr = !!(opt->flag & MM_F_SR), is_splice = !!(opt->flag & MM_F_SPLICE);
 	int32_t rid = a[r->as].x<<1>>33, rev = a[r->as].x>>63, as1, cnt1;
@@ -351,23 +401,25 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	int32_t rs1, qs1, re1, qe1;
 	int8_t mat[25];

-	if (is_sr) assert(!mi->is_hpc); // HPC won't work with SR because with HPC we can't easily tell if there is a gap
+	if (is_sr) assert(!(mi->flag & MM_I_HPC)); // HPC won't work with SR because with HPC we can't easily tell if there is a gap

 	r2->cnt = 0;
 	if (r->cnt == 0) return;
 	ksw_gen_simple_mat(5, mat, opt->a, opt->b);
 	bw = (int)(opt->bw * 1.5 + 1.);

-	if (is_sr && !mi->is_hpc) {
+	if (is_sr && !(mi->flag & MM_I_HPC)) {
 		mm_max_stretch(opt, r, a, &as1, &cnt1);
 		rs = (int32_t)a[as1].x + 1 - (int32_t)(a[as1].y>>32&0xff);
 		qs = (int32_t)a[as1].y + 1 - (int32_t)(a[as1].y>>32&0xff);
 		re = (int32_t)a[as1+cnt1-1].x + 1;
 		qe = (int32_t)a[as1+cnt1-1].y + 1;
 	} else {
-		if (!is_splice)
+		if (is_splice) {
+			mm_fix_bad_ends_splice(km, opt, mi, r, mat, qlen, qseq0, a, &as1, &cnt1);
+		} else {
 			mm_fix_bad_ends(r, a, opt->bw, &as1, &cnt1);
-		else as1 = r->as, cnt1 = r->cnt;
+		}
 		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);
@@ -447,6 +499,14 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			re0 = re0 > re1? re0 : re1;
 		} else re0 = re, qe0 = qe;
 	}
+	if (a[r->as].y & MM_SEED_SELF) {
+		int max_ext = r->qs > r->rs? r->qs - r->rs : r->rs - r->qs;
+		if (r->rs - rs0 > max_ext) rs0 = r->rs - max_ext;
+		if (r->qs - qs0 > max_ext) qs0 = r->qs - max_ext;
+		max_ext = r->qe > r->re? r->qe - r->re : r->re - r->qe;
+		if (re0 - r->re > max_ext) re0 = r->re + max_ext;
+		if (qe0 - r->qe > max_ext) qe0 = r->qe + max_ext;
+	}

 	assert(re0 > rs0);
 	tseq = (uint8_t*)kmalloc(km, re0 - rs0);
@@ -470,7 +530,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int

 	for (i = is_sr? cnt1 - 1 : 1; i < cnt1; ++i) { // gap filling
 		if ((a[as1+i].y & (MM_SEED_IGNORE|MM_SEED_TANDEM)) && i != cnt1 - 1) continue;
-		if (is_sr && !mi->is_hpc) {
+		if (is_sr && !(mi->flag & MM_I_HPC)) {
 			re = (int32_t)a[as1 + i].x + 1;
 			qe = (int32_t)a[as1 + i].y + 1;
 		} else mm_adjust_minier(mi, qseq0, &a[as1 + i], &re, &qe);
@@ -533,7 +593,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	assert(re1 - rs1 <= re0 - rs0);
 	if (r->p) {
 		mm_idx_getseq(mi, rid, rs1, re1, tseq);
-		mm_update_extra(r, &qseq0[r->rev][qs1], qual0[r->rev]? &qual0[r->rev][qs1] : 0, tseq, mat, opt->q, opt->e);
+		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e);
 		if (rev && r->p->trans_strand)
 			r->p->trans_strand ^= 3; // flip to the read strand
 	}
@@ -541,10 +601,10 @@ 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], uint8_t *qual0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez)
+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, *qual;
+	uint8_t *tseq, *qseq;
 	int8_t mat[25];
 	void *qp;

@@ -562,7 +622,6 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	tseq = (uint8_t*)kmalloc(km, tl);
 	mm_idx_getseq(mi, r1->rid, r1->re, r2->rs, tseq);
 	qseq = &qseq0[!r1->rev][qlen - r2->qs];
-	qual = qual0[!r1->rev]? &qseq0[!r1->rev][qlen - r2->qs] : 0;

 	mm_seq_rev(ql, qseq);
 	mm_seq_rev(tl, tseq);
@@ -582,9 +641,10 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	r_inv->inv = 1;
 	r_inv->rev = !r1->rev;
 	r_inv->rid = r1->rid;
+	r_inv->div = -1.0f;
 	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;
-	mm_update_extra(r_inv, &qseq[q_off], qual? &qual[q_off] : 0, &tseq[t_off], mat, opt->q, opt->e);
+	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e);
 	ret = 1;
 end_align1_inv:
 	kfree(km, tseq);
@@ -601,11 +661,11 @@ static inline mm_reg1_t *mm_insert_reg(const mm_reg1_t *r, int i, int *n_regs, m
 	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, const char *qual, int *n_regs_, mm_reg1_t *regs, mm128_t *a)
+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, n_regs = *n_regs_, n_a;
-	uint8_t *qseq0[2], *qual0[2];
+	uint8_t *qseq0[2];
 	ksw_extz_t ez;

 	// encode the query sequence
@@ -615,12 +675,6 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 		qseq0[0][i] = seq_nt4_table[(uint8_t)qstr[i]];
 		qseq0[1][qlen - 1 - i] = qseq0[0][i] < 4? 3 - qseq0[0][i] : 4;
 	}
-	if (qual) {
-		qual0[0] = (uint8_t*)kmalloc(km, qlen * 2);
-		qual0[1] = qual0[0] + qlen;
-		for (i = 0; i < qlen; ++i)
-			qual0[0][i] = qual0[1][qlen - 1 - i] = qual[i] - 33;
-	} else qual0[0] = qual0[1] = 0;

 	// align through seed hits
 	n_a = mm_squeeze_a(km, n_regs, regs, a);
@@ -631,8 +685,8 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 			mm_reg1_t s[2], s2[2];
 			int which, trans_strand;
 			s[0] = s[1] = regs[i];
-			mm_align1(km, opt, mi, qlen, qseq0, qual0, &s[0], &s2[0], n_a, a, &ez, MM_F_SPLICE_FOR);
-			mm_align1(km, opt, mi, qlen, qseq0, qual0, &s[1], &s2[1], n_a, a, &ez, MM_F_SPLICE_REV);
+			mm_align1(km, opt, mi, qlen, qseq0, &s[0], &s2[0], n_a, a, &ez, MM_F_SPLICE_FOR);
+			mm_align1(km, opt, mi, qlen, qseq0, &s[1], &s2[1], n_a, a, &ez, MM_F_SPLICE_REV);
 			if (s[0].p->dp_score > s[1].p->dp_score) which = 0, trans_strand = 1;
 			else if (s[0].p->dp_score < s[1].p->dp_score) which = 1, trans_strand = 2;
 			else trans_strand = 3, which = (qlen + s[0].p->dp_score) & 1; // randomly choose a strand, effectively
@@ -645,13 +699,13 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 			}
 			regs[i].p->trans_strand = trans_strand;
 		} else { // one round of alignment
-			mm_align1(km, opt, mi, qlen, qseq0, qual0, &regs[i], &r2, n_a, a, &ez, opt->flag);
+			mm_align1(km, opt, mi, qlen, qseq0, &regs[i], &r2, n_a, a, &ez, opt->flag);
 			if (opt->flag&MM_F_SPLICE)
 				regs[i].p->trans_strand = opt->flag&MM_F_SPLICE_FOR? 1 : 2;
 		}
 		if (r2.cnt > 0) regs = mm_insert_reg(&r2, i, &n_regs, regs);
-		if (!(opt->flag&MM_F_SPLICE) && !(opt->flag&MM_F_SR) && i > 0) { // don't try inversion alignment for -xsplice or -xsr
-			if (mm_align1_inv(km, opt, mi, qlen, qseq0, qual0, &regs[i-1], &regs[i], &r2, &ez)) {
+		if (!(opt->flag&(MM_F_SPLICE|MM_F_SR)) && !(opt->flag&(MM_F_FOR_ONLY|MM_F_REV_ONLY)) && i > 0) { // don't try inversion alignment for -xsplice or -xsr, or --for-only/rev-only
+			if (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
 			}
@@ -659,7 +713,6 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 	}
 	*n_regs_ = n_regs;
 	kfree(km, qseq0[0]);
-	if (qual0[0]) kfree(km, qual0[0]);
 	kfree(km, ez.cigar);
 	mm_filter_regs(km, opt, n_regs_, regs);
 	mm_hit_sort_by_dp(km, n_regs_, regs);
@@ -2,6 +2,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <assert.h>
+#define __STDC_LIMIT_MACROS
 #include "bseq.h"
 #include "kvec.h"
 #include "kseq.h"
@@ -114,10 +115,15 @@ mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int
 	*n_ = 0;
 	if (n_fp < 1) return 0;
 	while (1) {
+		int n_read = 0;
 		for (i = 0; i < n_fp; ++i)
-			if (kseq_read(fp[i]->ks) < 0)
-				break;
-		if (i != n_fp) break; // some file reaches the end
+			if (kseq_read(fp[i]->ks) >= 0)
+				++n_read;
+		if (n_read < n_fp) {
+			if (n_read > 0)
+				fprintf(stderr, "[W::%s]\033[1;31m query files have different number of records; extra records skipped.\033[0m\n", __func__);
+			break; // some file reaches the end
+		}
 		if (a.m == 0) kv_resize(mm_bseq1_t, 0, a, 256);
 		for (i = 0; i < n_fp; ++i) {
 			mm_bseq1_t *s;
@@ -0,0 +1,64 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+#include "mmpriv.h"
+
+static inline int32_t get_for_qpos(int32_t qlen, const mm128_t *a)
+{
+	int32_t x = (int32_t)a->y;
+	int32_t q_span = a->y>>32 & 0xff;
+	if (a->x>>63)
+		x = qlen - 1 - (x + 1 - q_span); // revert the position to the forward strand of query
+	return x;
+}
+
+static int get_mini_idx(int qlen, const mm128_t *a, int32_t n, const uint64_t *mini_pos)
+{
+	int32_t x, L = 0, R = n - 1;
+	x = get_for_qpos(qlen, a);
+	while (L <= R) { // binary search
+		int32_t m = ((uint64_t)L + R) >> 1;
+		int32_t y = (int32_t)mini_pos[m];
+		if (y < x) L = m + 1;
+		else if (y > x) R = m - 1;
+		else return m;
+	}
+	return -1;
+}
+
+void mm_est_err(const mm_idx_t *mi, int qlen, int n_regs, mm_reg1_t *regs, const mm128_t *a, int32_t n, const uint64_t *mini_pos)
+{
+	int i;
+	uint64_t sum_k = 0;
+	float avg_k;
+
+	if (n == 0) return;
+	for (i = 0; i < n; ++i)
+		sum_k += mini_pos[i] >> 32 & 0xff;
+	avg_k = (float)sum_k / n;
+
+	for (i = 0; i < n_regs; ++i) {
+		mm_reg1_t *r = &regs[i];
+		int32_t st, en, j, k, n_match, n_tot, l_ref;
+		r->div = -1.0f;
+		if (r->cnt == 0) continue;
+		st = en = get_mini_idx(qlen, r->rev? &a[r->as + r->cnt - 1] : &a[r->as], n, mini_pos);
+		if (st < 0) {
+			if (mm_verbose >= 2)
+				fprintf(stderr, "[WARNING] logic inconsistency in mm_est_err(). Please contact the developer.\n");
+			continue;
+		}
+		l_ref = mi->seq[r->rid].len;
+		for (k = 1, j = st + 1, n_match = 1; j < n && k < r->cnt; ++j) {
+			int32_t x;
+			x = get_for_qpos(qlen, r->rev? &a[r->as + r->cnt - 1 - k] : &a[r->as + k]);
+			if (x == (int32_t)mini_pos[j])
+				++k, en = j, ++n_match;
+		}
+		n_tot = en - st + 1;
+		if (r->qs > avg_k && r->rs > avg_k) ++n_tot;
+		if (qlen - r->qs > avg_k && l_ref - r->re > avg_k) ++n_tot;
+		r->div = logf((float)n_tot / n_match) / avg_k;
+	}
+}
@@ -202,7 +202,9 @@ static void write_cs(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_

 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 {
-	int type = r->inv? 'I' : r->id == r->parent? 'P' : 'S';
+	int type;
+	if (r->id == r->parent) type = r->inv? 'I' : 'P';
+	else type = r->inv? 'i' : 'S';
 	if (r->p) {
 		mm_sprintf_lite(s, "\tNM:i:%d\tms:i:%d\tAS:i:%d\tnn:i:%d", r->blen - r->mlen + r->p->n_ambi, r->p->dp_max, r->p->dp_score, r->p->n_ambi);
 		if (r->p->trans_strand == 1 || r->p->trans_strand == 2)
@@ -210,6 +212,12 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 	}
 	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->div >= 0.0f && r->div <= 1.0f) {
+		char buf[8];
+		if (r->div == 0.0f) buf[0] = '0', buf[1] = 0;
+		else sprintf(buf, "%.4f", r->div);
+		mm_sprintf_lite(s, "\tdv:f:%s", buf);
+	}
 	if (r->split) mm_sprintf_lite(s, "\tzd:i:%d", r->split);
 }

@@ -80,6 +80,7 @@ mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u,
 		ri->hash = (uint32_t)z[i].x;
 		ri->cnt = (int32_t)z[i].y;
 		ri->as = z[i].y >> 32;
+		ri->div = -1.0f;
 		mm_reg_set_coor(ri, qlen, a);
 	}
 	kfree(km, z);
@@ -147,7 +148,7 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff
 				ri->parent = rp->parent;
 				rp->subsc = rp->subsc > ri->score? rp->subsc : ri->score;
 				if (ri->cnt >= rp->cnt) cnt_sub = 1;
-				if (rp->p && ri->p && (rp->rs != ri->rs || rp->re != ri->re || ol != min)) { // the last condition excludes identical hits after DP
+				if (rp->p && ri->p && (rp->rid != ri->rid || rp->rs != ri->rs || rp->re != ri->re || ol != min)) { // the last condition excludes identical hits after DP
 					rp->p->dp_max2 = rp->p->dp_max2 > ri->p->dp_max? rp->p->dp_max2 : ri->p->dp_max;
 					if (rp->p->dp_max - ri->p->dp_max <= sub_diff) cnt_sub = 1;
 				}
@@ -234,7 +235,7 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_,
 			if (p == i || r[i].inv) { // primary or inversion
 				r[k++] = r[i];
 			} else if ((r[i].score >= r[p].score * pri_ratio || r[i].score + min_diff >= r[p].score) && n_2nd < best_n) {
-				if (!(r[i].qs == r[p].qs && r[i].qe == r[p].qe && r[i].rs == r[p].rs && r[i].re == r[p].re)) // not identical hits
+				if (!(r[i].qs == r[p].qs && r[i].qe == r[p].qe && r[i].rid == r[p].rid && r[i].rs == r[p].rs && r[i].re == r[p].re)) // not identical hits
 					r[k++] = r[i], ++n_2nd;
 				else if (r[i].p) free(r[i].p);
 			} else if (r[i].p) free(r[i].p);
@@ -389,8 +390,10 @@ mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int
 	for (s = 0; s < n_segs; ++s) {
 		regs[s] = mm_gen_regs(km, hash, qlens[s], seg[s].n_u, seg[s].u, seg[s].a);
 		n_regs[s] = seg[s].n_u;
-		for (i = 0; i < n_regs[s]; ++i)
+		for (i = 0; i < n_regs[s]; ++i) {
 			regs[s][i].seg_split = 1;
+			regs[s][i].seg_id = s;
+		}
 	}
 	return seg;
 }
@@ -7,6 +7,7 @@
 #endif
 #include <fcntl.h>
 #include <stdio.h>
+#define __STDC_LIMIT_MACROS
 #include "kthread.h"
 #include "bseq.h"
 #include "minimap.h"
@@ -31,19 +32,19 @@ typedef struct mm_idx_bucket_s {
 void mm_idxopt_init(mm_idxopt_t *opt)
 {
 	memset(opt, 0, sizeof(mm_idxopt_t));
-	opt->k = 15, opt->w = 10, opt->is_hpc = 0;
+	opt->k = 15, opt->w = 10, opt->flag = 0;
 	opt->bucket_bits = 14;
 	opt->mini_batch_size = 50000000;
 	opt->batch_size = 4000000000ULL;
 }

-mm_idx_t *mm_idx_init(int w, int k, int b, int is_hpc)
+mm_idx_t *mm_idx_init(int w, int k, int b, int flag)
 {
 	mm_idx_t *mi;
 	if (k*2 < b) b = k * 2;
 	if (w < 1) w = 1;
 	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->w = w, mi->k = k, mi->b = b, mi->flag = flag;
 	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;
@@ -89,7 +90,7 @@ void mm_idx_stat(const mm_idx_t *mi)
 {
 	int i, n = 0, n1 = 0;
 	uint64_t sum = 0, len = 0;
-	fprintf(stderr, "[M::%s] kmer size: %d; skip: %d; is_HPC: %d; #seq: %d\n", __func__, mi->k, mi->w, mi->is_hpc, mi->n_seq);
+	fprintf(stderr, "[M::%s] kmer size: %d; skip: %d; is_hpc: %d; #seq: %d\n", __func__, mi->k, mi->w, mi->flag&MM_I_HPC, mi->n_seq);
 	for (i = 0; i < mi->n_seq; ++i)
 		len += mi->seq[i].len;
 	for (i = 0; i < 1<<mi->b; ++i)
@@ -214,7 +215,7 @@ static void mm_idx_post(mm_idx_t *mi, int n_threads)
 #include "bseq.h"

 typedef struct {
-	int mini_batch_size, keep_name;
+	int mini_batch_size;
 	uint64_t batch_size, sum_len;
 	mm_bseq_file_t *fp;
 	mm_idx_t *mi;
@@ -246,7 +247,6 @@ static void *worker_pipeline(void *shared, int step, void *in)
 		s->seq = mm_bseq_read(p->fp, p->mini_batch_size, 0, &s->n_seq); // read a mini-batch
 		if (s->seq) {
 			uint32_t old_m, m;
-			uint64_t sum_len, old_max_len, max_len;
 			assert((uint64_t)p->mi->n_seq + s->n_seq <= UINT32_MAX); // to prevent integer overflow
 			// make room for p->mi->seq
 			old_m = p->mi->n_seq, m = p->mi->n_seq + s->n_seq;
@@ -254,29 +254,34 @@ static void *worker_pipeline(void *shared, int step, void *in)
 			if (old_m != m)
 				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;
-			max_len = (p->sum_len + sum_len + 7) / 8;
-			kroundup64(old_max_len); kroundup64(max_len);
-			if (old_max_len != max_len) {
-				p->mi->S = (uint32_t*)realloc(p->mi->S, max_len * 4);
-				memset(&p->mi->S[old_max_len], 0, 4 * (max_len - old_max_len));
+			if (!(p->mi->flag & MM_I_NO_SEQ)) {
+				uint64_t sum_len, old_max_len, max_len;
+				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;
+				max_len = (p->sum_len + sum_len + 7) / 8;
+				kroundup64(old_max_len); kroundup64(max_len);
+				if (old_max_len != max_len) {
+					p->mi->S = (uint32_t*)realloc(p->mi->S, max_len * 4);
+					memset(&p->mi->S[old_max_len], 0, 4 * (max_len - old_max_len));
+				}
 			}
 			// populate p->mi->seq
 			for (i = 0; i < s->n_seq; ++i) {
 				mm_idx_seq_t *seq = &p->mi->seq[p->mi->n_seq];
 				uint32_t j;
-				if (p->keep_name) {
+				if (!(p->mi->flag & MM_I_NO_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;
 				// copy the sequence
-				for (j = 0; j < seq->len; ++j) { // TODO: this is not the fastest way, but let's first see if speed matters here
-					uint64_t o = p->sum_len + j;
-					int c = seq_nt4_table[(uint8_t)s->seq[i].seq[j]];
-					mm_seq4_set(p->mi->S, o, c);
+				if (!(p->mi->flag & MM_I_NO_SEQ)) {
+					for (j = 0; j < seq->len; ++j) { // TODO: this is not the fastest way, but let's first see if speed matters here
+						uint64_t o = p->sum_len + j;
+						int c = seq_nt4_table[(uint8_t)s->seq[i].seq[j]];
+						mm_seq4_set(p->mi->S, o, c);
+					}
 				}
 				// update p->sum_len and p->mi->n_seq
 				p->sum_len += seq->len;
@@ -288,7 +293,10 @@ static void *worker_pipeline(void *shared, int step, void *in)
        step_t *s = (step_t*)in;
 		for (i = 0; i < s->n_seq; ++i) {
 			mm_bseq1_t *t = &s->seq[i];
-			mm_sketch(0, t->seq, t->l_seq, p->mi->w, p->mi->k, t->rid, p->mi->is_hpc, &s->a);
+			if (t->l_seq > 0)
+				mm_sketch(0, t->seq, t->l_seq, p->mi->w, p->mi->k, t->rid, p->mi->flag&MM_I_HPC, &s->a);
+			else if (mm_verbose >= 2)
+				fprintf(stderr, "[WARNING] the length database sequence '%s' is 0\n", t->name);
 			free(t->seq); free(t->name);
 		}
 		free(s->seq); s->seq = 0;
@@ -301,16 +309,15 @@ static void *worker_pipeline(void *shared, int step, void *in)
    return 0;
 }

-mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int is_hpc, int mini_batch_size, int n_threads, uint64_t batch_size, int keep_name)
+mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int flag, int mini_batch_size, int n_threads, uint64_t batch_size)
 {
 	pipeline_t pl;
 	if (fp == 0 || mm_bseq_eof(fp)) return 0;
 	memset(&pl, 0, sizeof(pipeline_t));
 	pl.mini_batch_size = mini_batch_size < batch_size? mini_batch_size : batch_size;
-	pl.keep_name = keep_name;
 	pl.batch_size = batch_size;
 	pl.fp = fp;
-	pl.mi = mm_idx_init(w, k, b, is_hpc);
+	pl.mi = mm_idx_init(w, k, b, flag);

 	kt_pipeline(n_threads < 3? n_threads : 3, worker_pipeline, &pl, 3);
 	if (mm_verbose >= 3)
@@ -323,17 +330,60 @@ mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int is_hpc, int mi
 	return pl.mi;
 }

-mm_idx_t *mm_idx_build(const char *fn, int w, int k, int is_hpc, int n_threads) // a simpler interface
+mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads) // a simpler interface; deprecated
 {
 	mm_bseq_file_t *fp;
 	mm_idx_t *mi;
 	fp = mm_bseq_open(fn);
 	if (fp == 0) return 0;
-	mi = mm_idx_gen(fp, w, k, 14, is_hpc, 1<<18, n_threads, UINT64_MAX, 1);
+	mi = mm_idx_gen(fp, w, k, 14, flag, 1<<18, n_threads, UINT64_MAX);
 	mm_bseq_close(fp);
 	return mi;
 }

+mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const char **seq, const char **name)
+{
+	uint64_t sum_len = 0;
+	mm128_v a = {0,0,0};
+	mm_idx_t *mi;
+	int i, flag = 0;
+	if (n <= 0) return 0;
+	for (i = 0; i < n; ++i) // get the total length
+		sum_len += strlen(seq[i]);
+	if (is_hpc) flag |= MM_I_HPC;
+	if (name == 0) flag |= MM_I_NO_NAME;
+	if (bucket_bits < 0) bucket_bits = 14;
+	mi = mm_idx_init(w, k, bucket_bits, flag);
+	mi->n_seq = n;
+	mi->seq = (mm_idx_seq_t*)kcalloc(mi->km, n, sizeof(mm_idx_seq_t)); // ->seq is allocated from km
+	mi->S = (uint32_t*)calloc((sum_len + 7) / 8, 4);
+	for (i = 0, sum_len = 0; i < n; ++i) {
+		const char *s = seq[i];
+		mm_idx_seq_t *p = &mi->seq[i];
+		uint32_t j;
+		if (name && name[i]) {
+			p->name = (char*)kmalloc(mi->km, strlen(name[i]) + 1);
+			strcpy(p->name, name[i]);
+		}
+		p->offset = sum_len;
+		p->len = strlen(s);
+		for (j = 0; j < p->len; ++j) {
+			int c = seq_nt4_table[(uint8_t)s[j]];
+			uint64_t o = sum_len + j;
+			mm_seq4_set(mi->S, o, c);
+		}
+		sum_len += p->len;
+		if (p->len > 0) {
+			a.n = 0;
+			mm_sketch(0, s, p->len, w, k, i, is_hpc, &a);
+			mm_idx_add(mi, a.n, a.a);
+		}
+	}
+	free(a.a);
+	mm_idx_post(mi, 1);
+	return mi;
+}
+
 /*************
 * index I/O *
 *************/
@@ -344,7 +394,7 @@ void mm_idx_dump(FILE *fp, const mm_idx_t *mi)
 	uint32_t x[5];
 	int i;

-	x[0] = mi->w, x[1] = mi->k, x[2] = mi->b, x[3] = mi->n_seq, x[4] = mi->is_hpc;
+	x[0] = mi->w, x[1] = mi->k, x[2] = mi->b, x[3] = mi->n_seq, x[4] = mi->flag;
 	fwrite(MM_IDX_MAGIC, 1, 4, fp);
 	fwrite(x, 4, 5, fp);
 	for (i = 0; i < mi->n_seq; ++i) {
@@ -371,7 +421,8 @@ void mm_idx_dump(FILE *fp, const mm_idx_t *mi)
 			fwrite(x, 8, 2, fp);
 		}
 	}
-	fwrite(mi->S, 4, (sum_len + 7) / 8, fp);
+	if (!(mi->flag & MM_I_NO_SEQ))
+		fwrite(mi->S, 4, (sum_len + 7) / 8, fp);
 	fflush(fp);
 }

@@ -421,8 +472,10 @@ mm_idx_t *mm_idx_load(FILE *fp)
 			kh_val(h, k) = x[1];
 		}
 	}
-	mi->S = (uint32_t*)malloc((sum_len + 7) / 8 * 4);
-	fread(mi->S, 4, (sum_len + 7) / 8, fp);
+	if (!(mi->flag & MM_I_NO_SEQ)) {
+		mi->S = (uint32_t*)malloc((sum_len + 7) / 8 * 4);
+		fread(mi->S, 4, (sum_len + 7) / 8, fp);
+	}
 	return mi;
 }

@@ -476,10 +529,10 @@ mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads)
 	mm_idx_t *mi;
 	if (r->is_idx) {
 		mi = mm_idx_load(r->fp.idx);
-		if (mi && mm_verbose >= 2 && (mi->k != r->opt.k || mi->w != r->opt.w || mi->is_hpc != r->opt.is_hpc))
+		if (mi && mm_verbose >= 2 && (mi->k != r->opt.k || mi->w != r->opt.w || (mi->flag&MM_I_HPC) != (r->opt.flag&MM_I_HPC)))
 			fprintf(stderr, "[WARNING]\033[1;31m Indexing parameters (-k, -w or -H) overridden by parameters used in the prebuilt index.\033[0m\n");
 	} else
-		mi = mm_idx_gen(r->fp.seq, r->opt.w, r->opt.k, r->opt.bucket_bits, r->opt.is_hpc, r->opt.mini_batch_size, n_threads, r->opt.batch_size, 1);
+		mi = mm_idx_gen(r->fp.seq, r->opt.w, r->opt.k, r->opt.bucket_bits, r->opt.flag, r->opt.mini_batch_size, n_threads, r->opt.batch_size);
 	if (mi) {
 		if (r->fp_out) mm_idx_dump(r->fp_out, mi);
 		++r->n_parts;
@@ -189,6 +189,10 @@ void km_stat(const void *_km, km_stat_t *s)
 			panic("[km_stat] The end of a free block enters another free block.");
 		if (p->ptr == km->loop_head) break;
 	}
-	for (p = km->core_head; p != NULL; p = p->ptr)
-		++s->n_cores, s->capacity += p->size * sizeof(header_t);
+	for (p = km->core_head; p != NULL; p = p->ptr) {
+		size_t size = p->size * sizeof(header_t);
+		++s->n_cores;
+		s->capacity += size;
+		s->largest = s->largest > size? s->largest : size;
+	}
 }
@@ -8,7 +8,7 @@ extern "C" {
 #endif

 typedef struct {
-	size_t capacity, available, n_blocks, n_cores;
+	size_t capacity, available, n_blocks, n_cores, largest;
 } km_stat_t;

 void *kmalloc(void *km, size_t size);
@@ -39,7 +39,24 @@ typedef struct {

 #define KSORT_SWAP(type_t, a, b) { register type_t t=(a); (a)=(b); (b)=t; }

-#define KSORT_INIT(name, type_t, __sort_lt)								\
+#define KSORT_INIT(name, type_t, __sort_lt) \
+	void ks_heapdown_##name(size_t i, size_t n, type_t l[]) \
+	{ \
+		size_t k = i; \
+		type_t tmp = l[i]; \
+		while ((k = (k << 1) + 1) < n) { \
+			if (k != n - 1 && __sort_lt(l[k], l[k+1])) ++k; \
+			if (__sort_lt(l[k], tmp)) break; \
+			l[i] = l[k]; i = k; \
+		} \
+		l[i] = tmp; \
+	} \
+	void ks_heapmake_##name(size_t lsize, type_t l[]) \
+	{ \
+		size_t i; \
+		for (i = (lsize >> 1) - 1; i != (size_t)(-1); --i) \
+			ks_heapdown_##name(i, lsize, l); \
+	} \
 	type_t ks_ksmall_##name(size_t n, type_t arr[], size_t kk)			\
 	{																	\
 		type_t *low, *high, *k, *ll, *hh, *mid;							\
@@ -116,7 +116,7 @@ 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, int with_N, const uint8_t *p, const int *off, const int *off_end, int n_col, int i0, int j0,
+static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int min_intron_len, 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;
@@ -138,11 +138,11 @@ static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int with_N, c
 		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 && !with_N)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i; // deletion
-		else if (state == 3 && with_N) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 3, 1), --i; // intron
+		else if (state == 1 || (state == 3 && min_intron_len <= 0)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i; // deletion
+		else if (state == 3 && min_intron_len > 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 3, 1), --i; // intron
 		else cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, 1), --j; // insertion
 	}
-	if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, i + 1); // first deletion
+	if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, min_intron_len > 0 && i >= min_intron_len? 3 : 2, i + 1); // first deletion
 	if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, j + 1); // first insertion
 	if (!is_rev)
 		for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR
@@ -367,9 +367,9 @@ void ksw_exts2_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, 1, (uint8_t*)p, off, off_end, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+			ksw_backtrack(km, 1, rev_cigar, long_thres, (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, 1, (uint8_t*)p, off, off_end, n_col_*16, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+			ksw_backtrack(km, 1, rev_cigar, long_thres, (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);
 	}
 }
@@ -2,6 +2,7 @@
 #include <stdlib.h>
 #include <limits.h>
 #include <stdint.h>
+#include "kthread.h"

 #if (defined(WIN32) || defined(_WIN32)) && defined(_MSC_VER)
 #define __sync_fetch_and_add(ptr, addend)     _InterlockedExchangeAdd((void*)ptr, addend)
@@ -6,7 +6,7 @@
 #include "mmpriv.h"
 #include "getopt.h"

-#define MM_VERSION "2.5-r572"
+#define MM_VERSION "2.8-r672"

 #ifdef __linux__
 #include <sys/resource.h>
@@ -28,7 +28,7 @@ static struct option long_options[] = {
 	{ "seed",           required_argument, 0, 0 },
 	{ "no-kalloc",      no_argument,       0, 0 },
 	{ "print-qname",    no_argument,       0, 0 },
-	{ "no-self",        no_argument,       0, 0 },
+	{ "no-self",        no_argument,       0, 'D' },
 	{ "print-seeds",    no_argument,       0, 0 },
 	{ "max-chain-skip", required_argument, 0, 0 },
 	{ "min-dp-len",     required_argument, 0, 0 },
@@ -37,12 +37,19 @@ static struct option long_options[] = {
 	{ "cost-non-gt-ag", required_argument, 0, 'C' },
 	{ "no-long-join",   no_argument,       0, 0 },
 	{ "sr",             no_argument,       0, 0 },
-	{ "frag",           optional_argument, 0, 0 },
-	{ "secondary",      optional_argument, 0, 0 },
+	{ "frag",           required_argument, 0, 0 },
+	{ "secondary",      required_argument, 0, 0 },
 	{ "cs",             optional_argument, 0, 0 },
 	{ "end-bonus",      required_argument, 0, 0 },
 	{ "no-pairing",     no_argument,       0, 0 },
-	{ "splice-flank",   optional_argument, 0, 0 },
+	{ "splice-flank",   required_argument, 0, 0 },
+	{ "idx-no-seq",     no_argument,       0, 0 },
+	{ "end-seed-pen",   required_argument, 0, 0 },   // 21
+	{ "for-only",       no_argument,       0, 0 },   // 22
+	{ "rev-only",       no_argument,       0, 0 },   // 23
+	{ "heap-sort",      required_argument, 0, 0 },   // 24
+	{ "all-chain",      no_argument,       0, 'P' },
+	{ "dual",           required_argument, 0, 0 },   // 26
 	{ "help",           no_argument,       0, 'h' },
 	{ "max-intron-len", required_argument, 0, 'G' },
 	{ "version",        no_argument,       0, 'V' },
@@ -65,9 +72,22 @@ static inline int64_t mm_parse_num(const char *str)
 	return (int64_t)(x + .499);
 }

+static inline void yes_or_no(mm_mapopt_t *opt, int flag, int long_idx, const char *arg, int yes_to_set)
+{
+	if (yes_to_set) {
+		if (strcmp(arg, "yes") == 0 || strcmp(arg, "y") == 0) opt->flag |= flag;
+		else if (strcmp(arg, "no") == 0 || strcmp(arg, "n") == 0) opt->flag &= ~flag;
+		else fprintf(stderr, "[WARNING]\033[1;31m option '--%s' only accepts 'yes' or 'no'.\033[0m\n", long_options[long_idx].name);
+	} else {
+		if (strcmp(arg, "yes") == 0 || strcmp(arg, "y") == 0) opt->flag &= ~flag;
+		else if (strcmp(arg, "no") == 0 || strcmp(arg, "n") == 0) opt->flag |= flag;
+		else fprintf(stderr, "[WARNING]\033[1;31m option '--%s' only accepts 'yes' or 'no'.\033[0m\n", long_options[long_idx].name);
+	}
+}
+
 int main(int argc, char *argv[])
 {
-	const char *opt_str = "2aSw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:";
+	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:";
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
 	int i, c, n_threads = 3, long_idx;
@@ -94,7 +114,7 @@ int main(int argc, char *argv[])
 	while ((c = getopt_long(argc, argv, opt_str, long_options, &long_idx)) >= 0) {
 		if (c == 'w') ipt.w = atoi(optarg);
 		else if (c == 'k') ipt.k = atoi(optarg);
-		else if (c == 'H') ipt.is_hpc = 1;
+		else if (c == 'H') ipt.flag |= MM_I_HPC;
 		else if (c == 'd') fnw = optarg; // the above are indexing related options, except -I
 		else if (c == 'r') opt.bw = (int)mm_parse_num(optarg);
 		else if (c == 't') n_threads = atoi(optarg);
@@ -106,7 +126,9 @@ int main(int argc, char *argv[])
 		else if (c == 'p') opt.pri_ratio = atof(optarg);
 		else if (c == 'M') opt.mask_level = atof(optarg);
 		else if (c == 'c') opt.flag |= MM_F_OUT_CG | MM_F_CIGAR;
-		else if (c == 'X') opt.flag |= MM_F_AVA | MM_F_NO_SELF;
+		else if (c == 'D') opt.flag |= MM_F_NO_DIAG;
+		else if (c == 'P') opt.flag |= MM_F_ALL_CHAINS;
+		else if (c == 'X') opt.flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN; // -D -P --no-long-join --dual=no
 		else if (c == 'a') opt.flag |= MM_F_OUT_SAM | MM_F_CIGAR;
 		else if (c == 'Q') opt.flag |= MM_F_NO_QUAL;
 		else if (c == 'Y') opt.flag |= MM_F_SOFTCLIP;
@@ -128,24 +150,23 @@ int main(int argc, char *argv[])
 		else if (c == 0 && long_idx == 2) opt.seed = atoi(optarg); // --seed
 		else if (c == 0 && long_idx == 3) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
 		else if (c == 0 && long_idx == 4) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
-		else if (c == 0 && long_idx == 5) opt.flag |= MM_F_NO_SELF; // --no-self
 		else if (c == 0 && long_idx == 6) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --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 == 0 && long_idx == 9) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
 		else if (c == 0 && long_idx ==10) opt.flag |= MM_F_SPLICE; // --splice
 		else if (c == 0 && long_idx ==12) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
 		else if (c == 0 && long_idx ==13) opt.flag |= MM_F_SR; // --sr
 		else if (c == 0 && long_idx ==17) opt.end_bonus = atoi(optarg); // --end-bonus
 		else if (c == 0 && long_idx ==18) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
+		else if (c == 0 && long_idx ==20) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
+		else if (c == 0 && long_idx ==21) opt.anchor_ext_shift = atoi(optarg); // --end-seed-pen
+		else if (c == 0 && long_idx ==22) opt.flag |= MM_F_FOR_ONLY; // --for-only
+		else if (c == 0 && long_idx ==23) opt.flag |= MM_F_REV_ONLY; // --rev-only
 		else if (c == 0 && long_idx == 14) { // --frag
-			if (optarg == 0 || strcmp(optarg, "yes") == 0 || strcmp(optarg, "y") == 0)
-				opt.flag |= MM_F_FRAG_MODE;
-			else opt.flag &= ~MM_F_FRAG_MODE;
+			yes_or_no(&opt, MM_F_FRAG_MODE, long_idx, optarg, 1);
 		} else if (c == 0 && long_idx == 15) { // --secondary
-			if (optarg == 0 || strcmp(optarg, "yes") == 0 || strcmp(optarg, "y") == 0)
-				opt.flag &= ~MM_F_NO_PRINT_2ND;
-			else opt.flag |= MM_F_NO_PRINT_2ND;
+			yes_or_no(&opt, MM_F_NO_PRINT_2ND, long_idx, optarg, 0);
 		} else if (c == 0 && long_idx == 16) { // --cs
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR;
 			if (optarg == 0 || strcmp(optarg, "short") == 0) {
@@ -158,9 +179,11 @@ int main(int argc, char *argv[])
 				fprintf(stderr, "[WARNING]\033[1;31m --cs only takes 'short' or 'long'. Invalid values are assumed to be 'short'.\033[0m\n");
 			}
 		} else if (c == 0 && long_idx == 19) { // --splice-flank
-			if (optarg == 0 || strcmp(optarg, "yes") == 0 || strcmp(optarg, "y") == 0)
-				opt.flag |= MM_F_SPLICE_FLANK;
-			else opt.flag &= ~MM_F_SPLICE_FLANK;
+			yes_or_no(&opt, MM_F_SPLICE_FLANK, long_idx, optarg, 1);
+		} else if (c == 0 && long_idx == 24) { // --heap-sort
+			yes_or_no(&opt, MM_F_HEAP_SORT, long_idx, optarg, 1);
+		} else if (c == 0 && long_idx == 26) { // --dual
+			yes_or_no(&opt, MM_F_NO_DUAL, long_idx, optarg, 0);
 		} else if (c == 'S') {
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR | MM_F_OUT_CS_LONG;
 			if (mm_verbose >= 2)
@@ -196,6 +219,10 @@ int main(int argc, char *argv[])
 		fprintf(stderr, "[ERROR]\033[1;31m --splice and --frag should not be specified at the same time.\033[0m\n");
 		return 1;
 	}
+	if (!fnw && !(opt.flag&MM_F_CIGAR))
+		ipt.flag |= MM_I_NO_SEQ;
+	if (mm_check_opt(&ipt, &opt) < 0)
+		return 1;

 	if (argc == optind || fp_help == stdout) {
 		fprintf(fp_help, "Usage: minimap2 [options] <target.fa>|<target.idx> [query.fa] [...]\n");
@@ -244,14 +271,18 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "                 map-ont: -k15 (Oxford Nanopore vs reference mapping)\n");
 		fprintf(fp_help, "                 asm5: -k19 -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200 (asm to ref mapping; break at 5%% div.)\n");
 		fprintf(fp_help, "                 asm10: -k19 -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200 (asm to ref mapping; break at 10%% div.)\n");
-		fprintf(fp_help, "                 ava-pb: -Hk19 -w5 -Xp0 -m100 -g10000 --max-chain-skip 25 (PacBio read overlap)\n");
-		fprintf(fp_help, "                 ava-ont: -k15 -w5 -Xp0 -m100 -g10000 --max-chain-skip 25 (ONT read overlap)\n");
+		fprintf(fp_help, "                 ava-pb: -Hk19 -Xw5 -m100 -g10000 --max-chain-skip 25 (PacBio read overlap)\n");
+		fprintf(fp_help, "                 ava-ont: -k15 -Xw5 -m100 -g10000 --max-chain-skip 25 (ONT read overlap)\n");
 		fprintf(fp_help, "                 splice: long-read spliced alignment (see minimap2.1 for details)\n");
 		fprintf(fp_help, "                 sr: short single-end reads without splicing (see minimap2.1 for details)\n");
 		fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of command-line options.\n");
 		return fp_help == stdout? 0 : 1;
 	}

+	if ((opt.flag & MM_F_SR) && argc - optind > 3) {
+		fprintf(stderr, "[ERROR] incorrect input: in the sr mode, please specify no more than two query files.\n");
+		return 1;
+	}
 	idx_rdr = mm_idx_reader_open(argv[optind], &ipt, fnw);
 	if (idx_rdr == 0) {
 		fprintf(stderr, "[ERROR] failed to open file '%s'\n", argv[optind]);
@@ -259,11 +290,18 @@ int main(int argc, char *argv[])
 	}
 	if (!idx_rdr->is_idx && fnw == 0 && argc - optind < 2) {
 		fprintf(stderr, "[ERROR] missing input: please specify a query file to map or option -d to keep the index\n");
+		mm_idx_reader_close(idx_rdr);
 		return 1;
 	}
 	if (opt.best_n == 0 && (opt.flag&MM_F_CIGAR) && mm_verbose >= 2)
 		fprintf(stderr, "[WARNING]\033[1;31m `-N 0' reduces alignment accuracy. Please use --secondary=no to suppress secondary alignments.\033[0m\n");
 	while ((mi = mm_idx_reader_read(idx_rdr, n_threads)) != 0) {
+		if ((opt.flag & MM_F_CIGAR) && (mi->flag & MM_I_NO_SEQ)) {
+			fprintf(stderr, "[ERROR] the prebuilt index doesn't contain sequences.\n");
+			mm_idx_destroy(mi);
+			mm_idx_reader_close(idx_rdr);
+			return 1;
+		}
 		if ((opt.flag & MM_F_OUT_SAM) && idx_rdr->n_parts == 1) {
 			if (mm_idx_reader_eof(idx_rdr)) {
 				mm_write_sam_hdr(mi, rg, MM_VERSION, argc, argv);
@@ -9,121 +9,7 @@
 #include "bseq.h"
 #include "khash.h"

-void mm_mapopt_init(mm_mapopt_t *opt)
-{
-	memset(opt, 0, sizeof(mm_mapopt_t));
-	opt->seed = 11;
-	opt->mid_occ_frac = 2e-4f;
-	opt->sdust_thres = 0; // no SDUST masking
-
-	opt->min_cnt = 3;
-	opt->min_chain_score = 40;
-	opt->bw = 500;
-	opt->max_gap = 5000;
-	opt->max_gap_ref = -1;
-	opt->max_chain_skip = 25;
-
-	opt->mask_level = 0.5f;
-	opt->pri_ratio = 0.8f;
-	opt->best_n = 5;
-
-	opt->max_join_long = 20000;
-	opt->max_join_short = 2000;
-	opt->min_join_flank_sc = 1000;
-
-	opt->a = 2, opt->b = 4, opt->q = 4, opt->e = 2, opt->q2 = 24, opt->e2 = 1;
-	opt->zdrop = 400;
-	opt->end_bonus = -1;
-	opt->min_dp_max = opt->min_chain_score * opt->a;
-	opt->min_ksw_len = 200;
-	opt->mini_batch_size = 500000000;
-
-	opt->pe_ori = 0; // FF
-	opt->pe_bonus = 33;
-}
-
-void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)
-{
-	if ((opt->flag & MM_F_SPLICE_FOR) && (opt->flag & MM_F_SPLICE_REV))
-		opt->flag |= MM_F_SPLICE;
-	if (opt->mid_occ <= 0)
-		opt->mid_occ = mm_idx_cal_max_occ(mi, opt->mid_occ_frac);
-	if (mm_verbose >= 3)
-		fprintf(stderr, "[M::%s::%.3f*%.2f] mid_occ = %d\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), opt->mid_occ);
-}
-
-void mm_mapopt_max_intron_len(mm_mapopt_t *opt, int max_intron_len)
-{
-	if ((opt->flag & MM_F_SPLICE) && max_intron_len > 0)
-		opt->max_gap_ref = opt->bw = max_intron_len;
-}
-
-int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
-{
-	if (preset == 0) {
-		mm_idxopt_init(io);
-		mm_mapopt_init(mo);
-	} else if (strcmp(preset, "ava-ont") == 0) {
-		io->is_hpc = 0, io->k = 15, io->w = 5;
-		mo->flag |= MM_F_AVA | MM_F_NO_SELF;
-		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
-	} else if (strcmp(preset, "ava-pb") == 0) {
-		io->is_hpc = 1, io->k = 19, io->w = 5;
-		mo->flag |= MM_F_AVA | MM_F_NO_SELF;
-		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
-	} else if (strcmp(preset, "map10k") == 0 || strcmp(preset, "map-pb") == 0) {
-		io->is_hpc = 1, io->k = 19;
-	} else if (strcmp(preset, "map-ont") == 0) {
-		io->is_hpc = 0, io->k = 15;
-	} else if (strcmp(preset, "asm5") == 0) {
-		io->is_hpc = 0, io->k = 19, io->w = 19;
-		mo->a = 1, mo->b = 19, mo->q = 39, mo->q2 = 81, mo->e = 3, mo->e2 = 1, mo->zdrop = 200;
-		mo->min_dp_max = 200;
-		mo->best_n = 50;
-	} else if (strcmp(preset, "asm10") == 0) {
-		io->is_hpc = 0, io->k = 19, io->w = 19;
-		mo->a = 1, mo->b = 9, mo->q = 16, mo->q2 = 41, mo->e = 2, mo->e2 = 1, mo->zdrop = 200;
-		mo->min_dp_max = 200;
-		mo->best_n = 50;
-	} else if (strcmp(preset, "short") == 0 || strcmp(preset, "sr") == 0) {
-		io->is_hpc = 0, io->k = 21, io->w = 11;
-		mo->flag |= MM_F_SR | MM_F_FRAG_MODE | MM_F_NO_PRINT_2ND | MM_F_2_IO_THREADS;
-		mo->pe_ori = 0<<1|1; // FR
-		mo->a = 2, mo->b = 8, mo->q = 12, mo->e = 2, mo->q2 = 24, mo->e2 = 1;
-		mo->zdrop = 100;
-		mo->end_bonus = 10;
-		mo->max_frag_len = 800;
-		mo->max_gap = 100;
-		mo->bw = 100;
-		mo->pri_ratio = 0.5f;
-		mo->min_cnt = 2;
-		mo->min_chain_score = 25;
-		mo->min_dp_max = 40;
-		mo->best_n = 20;
-		mo->mid_occ = 1000;
-		mo->max_occ = 5000;
-		mo->mini_batch_size = 50000000;
-	} else if (strcmp(preset, "splice") == 0 || strcmp(preset, "cdna") == 0) {
-		io->is_hpc = 0, io->k = 15, io->w = 5;
-		mo->flag |= MM_F_SPLICE | MM_F_SPLICE_FOR | MM_F_SPLICE_REV | MM_F_SPLICE_FLANK;
-		mo->max_gap = 2000, mo->max_gap_ref = mo->bw = 200000;
-		mo->a = 1, mo->b = 2, mo->q = 2, mo->e = 1, mo->q2 = 32, mo->e2 = 0;
-		mo->noncan = 9;
-		mo->zdrop = 200;
-	} else return -1;
-	return 0;
-}
-
-typedef struct {
-	uint32_t n;
-	uint32_t qpos;
-	uint32_t seg_id;
-	const uint64_t *cr;
-} mm_match_t;
-
 struct mm_tbuf_s {
-	sdust_buf_t *sdb;
-	mm128_v mini;
 	void *km;
 };

@@ -132,24 +18,23 @@ mm_tbuf_t *mm_tbuf_init(void)
 	mm_tbuf_t *b;
 	b = (mm_tbuf_t*)calloc(1, sizeof(mm_tbuf_t));
 	if (!(mm_dbg_flag & 1)) b->km = km_init();
-	b->sdb = sdust_buf_init(b->km);
 	return b;
 }

 void mm_tbuf_destroy(mm_tbuf_t *b)
 {
 	if (b == 0) return;
-	kfree(b->km, b->mini.a);
-	sdust_buf_destroy(b->sdb);
 	km_destroy(b->km);
 	free(b);
 }

-static int mm_dust_minier(int n, mm128_t *a, int l_seq, const char *seq, int sdust_thres, sdust_buf_t *sdb)
+static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *seq, int sdust_thres)
 {
 	int n_dreg, j, k, u = 0;
 	const uint64_t *dreg;
-	if (sdust_thres <= 0 || sdb == 0) return n;
+	sdust_buf_t *sdb;
+	if (sdust_thres <= 0) return n;
+	sdb = sdust_buf_init(km);
 	dreg = sdust_core((const uint8_t*)seq, l_seq, sdust_thres, 64, &n_dreg, sdb);
 	for (j = k = 0; j < n; ++j) { // squeeze out minimizers that significantly overlap with LCRs
 		int32_t qpos = (uint32_t)a[j].y>>1, span = a[j].x&0xff;
@@ -163,103 +48,209 @@ static int mm_dust_minier(int n, mm128_t *a, int l_seq, const char *seq, int sdu
 				l += ee - ss;
 			}
 			if (l <= span>>1) a[k++] = a[j]; // keep the minimizer if less than half of it falls in masked region
-		}
+		} else a[k++] = a[j];
 	}
+	sdust_buf_destroy(sdb);
 	return k; // the new size
 }

-static void collect_minimizers(const mm_mapopt_t *opt, const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, mm_tbuf_t *b)
+static void collect_minimizers(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, mm128_v *mv)
 {
 	int i, j, n, sum = 0;
-	b->mini.n = 0;
+	mv->n = 0;
 	for (i = n = 0; i < n_segs; ++i) {
-		mm_sketch(b->km, seqs[i], qlens[i], mi->w, mi->k, i, mi->is_hpc, &b->mini);
-		for (j = n; j < b->mini.n; ++j)
-			b->mini.a[j].y += sum << 1;
+		mm_sketch(km, seqs[i], qlens[i], mi->w, mi->k, i, mi->flag&MM_I_HPC, mv);
+		for (j = n; j < mv->n; ++j)
+			mv->a[j].y += sum << 1;
 		if (opt->sdust_thres > 0) // mask low-complexity minimizers
-			b->mini.n = n + mm_dust_minier(b->mini.n - n, b->mini.a + n, qlens[i], seqs[i], opt->sdust_thres, b->sdb);
-		sum += qlens[i], n = b->mini.n;
+			mv->n = n + mm_dust_minier(km, mv->n - n, mv->a + n, qlens[i], seqs[i], opt->sdust_thres);
+		sum += qlens[i], n = mv->n;
 	}
 }

-static mm128_t *collect_seed_hits(const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, int qlen, int64_t *n_a, int *rep_len, mm_tbuf_t *b)
-{
-	int rep_st = 0, rep_en = 0, i;
-	mm_match_t *m;
-	mm128_t *a;
+#include "ksort.h"
+#define heap_lt(a, b) ((a).x > (b).x)
+KSORT_INIT(heap, mm128_t, heap_lt)

-	m = (mm_match_t*)kmalloc(b->km, b->mini.n * sizeof(mm_match_t));
-	for (i = 0; i < b->mini.n; ++i) {
+typedef struct {
+	uint32_t n;
+	uint32_t q_pos, q_span;
+	uint32_t seg_id:31, is_tandem:1;
+	const uint64_t *cr;
+} mm_match_t;
+
+static mm_match_t *collect_matches(void *km, int *_n_m, int max_occ, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos)
+{
+	int i, rep_st = 0, rep_en = 0, n_m;
+	mm_match_t *m;
+	*n_mini_pos = 0;
+	*mini_pos = (uint64_t*)kmalloc(km, mv->n * sizeof(uint64_t));
+	m = (mm_match_t*)kmalloc(km, mv->n * sizeof(mm_match_t));
+	for (i = n_m = 0, *rep_len = 0, *n_a = 0; i < mv->n; ++i) {
+		const uint64_t *cr;
+		mm128_t *p = &mv->a[i];
+		uint32_t q_pos = (uint32_t)p->y, q_span = p->x & 0xff;
 		int t;
-		mm128_t *p = &b->mini.a[i];
-		m[i].qpos = (uint32_t)p->y;
-		m[i].cr = mm_idx_get(mi, p->x>>8, &t);
-		m[i].n = t;
-		m[i].seg_id = p->y >> 32;
-	}
-	for (i = 0, *n_a = 0; i < b->mini.n; ++i) // find the length of a[]
-		if (m[i].n < max_occ) *n_a += m[i].n;
-	a = (mm128_t*)kmalloc(b->km, *n_a * sizeof(mm128_t));
-	for (i = *rep_len = 0, *n_a = 0; i < b->mini.n; ++i) {
-		mm128_t *p = &b->mini.a[i];
-		mm_match_t *q = &m[i];
-		const uint64_t *r = q->cr;
-		int k, q_span = p->x & 0xff, is_tandem = 0;
-		if (q->n >= max_occ) {
-			int en = (q->qpos>>1) + 1, st = en - q_span;
+		cr = mm_idx_get(mi, p->x>>8, &t);
+		if (t >= max_occ) {
+			int en = (q_pos >> 1) + 1, st = en - q_span;
 			if (st > rep_en) {
 				*rep_len += rep_en - rep_st;
 				rep_st = st, rep_en = en;
 			} else rep_en = en;
-			continue;
-		}
-		if (i > 0 && p->x>>8 == b->mini.a[i - 1].x>>8) is_tandem = 1;
-		if (i < b->mini.n - 1 && p->x>>8 == b->mini.a[i + 1].x>>8) is_tandem = 1;
-		for (k = 0; k < q->n; ++k) {
-			int32_t rpos = (uint32_t)r[k] >> 1;
-			mm128_t *p;
-			if (qname && (opt->flag&(MM_F_NO_SELF|MM_F_AVA))) {
-				const char *tname = mi->seq[r[k]>>32].name;
-				int cmp;
-				cmp = strcmp(qname, tname);
-				if ((opt->flag&MM_F_NO_SELF) && cmp == 0 && rpos == (q->qpos>>1)) // avoid the diagonal
-					continue;
-				if ((opt->flag&MM_F_AVA) && cmp > 0) // all-vs-all mode: map once
-					continue;
-			}
-			p = &a[(*n_a)++];
-			if ((r[k]&1) == (q->qpos&1)) { // forward strand
-				p->x = (r[k]&0xffffffff00000000ULL) | rpos;
-				p->y = (uint64_t)q_span << 32 | q->qpos >> 1;
-			} else { // reverse strand
-				p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | rpos;
-				p->y = (uint64_t)q_span << 32 | (qlen - ((q->qpos>>1) + 1 - q_span) - 1);
-			}
-			p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
-			if (is_tandem) p->y |= MM_SEED_TANDEM;
+		} else {
+			mm_match_t *q = &m[n_m++];
+			q->q_pos = q_pos, q->q_span = q_span, q->cr = cr, q->n = t, q->seg_id = p->y >> 32;
+			q->is_tandem = 0;
+			if (i > 0 && p->x>>8 == mv->a[i - 1].x>>8) q->is_tandem = 1;
+			if (i < mv->n - 1 && p->x>>8 == mv->a[i + 1].x>>8) q->is_tandem = 1;
+			*n_a += q->n;
+			(*mini_pos)[(*n_mini_pos)++] = (uint64_t)q_span<<32 | q_pos>>1;
 		}
 	}
 	*rep_len += rep_en - rep_st;
-	kfree(b->km, m);
+	*_n_m = n_m;
+	return m;
+}
+
+static inline int skip_seed(int flag, uint64_t r, const mm_match_t *q, const char *qname, int qlen, const mm_idx_t *mi, int *is_self)
+{
+	*is_self = 0;
+	if (qname && (flag & (MM_F_NO_DIAG|MM_F_NO_DUAL))) {
+		const mm_idx_seq_t *s = &mi->seq[r>>32];
+		int cmp;
+		cmp = strcmp(qname, s->name);
+		if ((flag&MM_F_NO_DIAG) && cmp == 0 && s->len == qlen) {
+			if ((uint32_t)r>>1 == (q->q_pos>>1)) return 1; // avoid the diagnonal anchors
+			if ((r&1) == (q->q_pos&1)) *is_self = 1; // this flag is used to avoid spurious extension on self chain
+		}
+		if ((flag&MM_F_NO_DUAL) && cmp > 0) // all-vs-all mode: map once
+			return 1;
+	}
+	if (flag & (MM_F_FOR_ONLY|MM_F_REV_ONLY)) {
+		if ((r&1) == (q->q_pos&1)) { // forward strand
+			if (flag & MM_F_REV_ONLY) return 1;
+		} else {
+			if (flag & MM_F_FOR_ONLY) return 1;
+		}
+	}
+	return 0;
+}
+
+static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, const mm128_v *mv, int qlen, int64_t *n_a, int *rep_len,
+								  int *n_mini_pos, uint64_t **mini_pos)
+{
+	int i, n_m, heap_size = 0;
+	int64_t j, n_for = 0, n_rev = 0;
+	mm_match_t *m;
+	mm128_t *a, *heap;
+
+	m = collect_matches(km, &n_m, max_occ, mi, mv, n_a, rep_len, n_mini_pos, mini_pos);
+
+	heap = (mm128_t*)kmalloc(km, n_m * sizeof(mm128_t));
+	a = (mm128_t*)kmalloc(km, *n_a * sizeof(mm128_t));
+
+	for (i = 0, heap_size = 0; i < n_m; ++i) {
+		if (m[i].n > 0) {
+			heap[heap_size].x = m[i].cr[0];
+			heap[heap_size].y = (uint64_t)i<<32;
+			++heap_size;
+		}
+	}
+	ks_heapmake_heap(heap_size, heap);
+	while (heap_size > 0) {
+		mm_match_t *q = &m[heap->y>>32];
+		mm128_t *p;
+		uint64_t r = heap->x;
+		int32_t is_self, rpos = (uint32_t)r >> 1;
+		if (skip_seed(opt->flag, r, q, qname, qlen, mi, &is_self)) continue;
+		if ((r&1) == (q->q_pos&1)) { // forward strand
+			p = &a[n_for++];
+			p->x = (r&0xffffffff00000000ULL) | rpos;
+			p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
+		} else { // reverse strand
+			p = &a[(*n_a) - (++n_rev)];
+			p->x = 1ULL<<63 | (r&0xffffffff00000000ULL) | rpos;
+			p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
+		}
+		p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
+		if (q->is_tandem) p->y |= MM_SEED_TANDEM;
+		if (is_self) p->y |= MM_SEED_SELF;
+		// update the heap
+		if ((uint32_t)heap->y < q->n - 1) {
+			++heap[0].y;
+			heap[0].x = m[heap[0].y>>32].cr[(uint32_t)heap[0].y];
+		} else {
+			heap[0] = heap[heap_size - 1];
+			--heap_size;
+		}
+		ks_heapdown_heap(0, heap_size, heap);
+	}
+	kfree(km, m);
+	kfree(km, heap);
+
+	// reverse anchors on the reverse strand, as they are in the descending order
+	for (j = 0; j < n_rev>>1; ++j) {
+		mm128_t t = a[(*n_a) - 1 - j];
+		a[(*n_a) - 1 - j] = a[(*n_a) - (n_rev - j)];
+		a[(*n_a) - (n_rev - j)] = t;
+	}
+	if (*n_a > n_for + n_rev) {
+		memmove(a + n_for, a + (*n_a) - n_rev, n_rev * sizeof(mm128_t));
+		*n_a = n_for + n_rev;
+	}
+	return a;
+}
+
+static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, const mm128_v *mv, int qlen, int64_t *n_a, int *rep_len,
+								  int *n_mini_pos, uint64_t **mini_pos)
+{
+	int i, k, n_m;
+	mm_match_t *m;
+	mm128_t *a;
+	m = collect_matches(km, &n_m, max_occ, mi, mv, n_a, rep_len, n_mini_pos, mini_pos);
+	a = (mm128_t*)kmalloc(km, *n_a * sizeof(mm128_t));
+	for (i = 0, *n_a = 0; i < n_m; ++i) {
+		mm_match_t *q = &m[i];
+		const uint64_t *r = q->cr;
+		for (k = 0; k < q->n; ++k) {
+			int32_t is_self, rpos = (uint32_t)r[k] >> 1;
+			mm128_t *p;
+			if (skip_seed(opt->flag, r[k], q, qname, qlen, mi, &is_self)) continue;
+			p = &a[(*n_a)++];
+			if ((r[k]&1) == (q->q_pos&1)) { // forward strand
+				p->x = (r[k]&0xffffffff00000000ULL) | rpos;
+				p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
+			} else { // reverse strand
+				p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | rpos;
+				p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
+			}
+			p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
+			if (q->is_tandem) p->y |= MM_SEED_TANDEM;
+			if (is_self) p->y |= MM_SEED_SELF;
+		}
+	}
+	kfree(km, m);
+	radix_sort_128x(a, a + (*n_a));
 	return a;
 }

 static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_idx_t *mi, void *km, int qlen, int n_segs, const int *qlens, int *n_regs, mm_reg1_t *regs, mm128_t *a)
 {
-	if (!(opt->flag & MM_F_AVA)) { // don't choose primary mapping(s) for read overlap
+	if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
 		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b);
 		if (n_segs <= 1) mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 		else mm_select_sub_multi(km, opt->pri_ratio, 0.2f, 0.7f, max_chain_gap_ref, mi->k*2, opt->best_n, n_segs, qlens, n_regs, regs);
-		if (!(opt->flag & MM_F_SPLICE) && !(opt->flag & MM_F_SR) && !(opt->flag & MM_F_NO_LJOIN))
+		if (!(opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN))) // long join not working well without primary chains
 			mm_join_long(km, opt, qlen, n_regs, regs, a);
 	}
 }

-static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *km, int qlen, const char *seq, const char *qual, int *n_regs, mm_reg1_t *regs, mm128_t *a)
+static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *km, int qlen, const char *seq, int *n_regs, mm_reg1_t *regs, mm128_t *a)
 {
 	if (!(opt->flag & MM_F_CIGAR)) return regs;
-	regs = mm_align_skeleton(km, opt, mi, qlen, seq, qual, n_regs, regs, a); // this calls mm_filter_regs()
-	if (!(opt->flag & MM_F_AVA)) {
+	regs = mm_align_skeleton(km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs()
+	if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
 		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b);
 		mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 		mm_set_sam_pri(*n_regs, regs);
@@ -267,15 +258,17 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k
 	return regs;
 }

-void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, const char **quals, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
+void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
 {
-	int i, j, rep_len, qlen_sum, n_regs0;
+	int i, j, rep_len, qlen_sum, n_regs0, n_mini_pos;
 	int max_chain_gap_qry, max_chain_gap_ref, is_splice = !!(opt->flag & MM_F_SPLICE), is_sr = !!(opt->flag & MM_F_SR);
 	uint32_t hash;
 	int64_t n_a;
-	uint64_t *u;
+	uint64_t *u, *mini_pos;
 	mm128_t *a;
+	mm128_v mv = {0,0,0};
 	mm_reg1_t *regs0;
+	km_stat_t kmst;

 	for (i = 0, qlen_sum = 0; i < n_segs; ++i)
 		qlen_sum += qlens[i], n_regs[i] = 0, regs[i] = 0;
@@ -286,9 +279,9 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	hash ^= __ac_Wang_hash(qlen_sum) + __ac_Wang_hash(opt->seed);
 	hash  = __ac_Wang_hash(hash);

-	collect_minimizers(opt, mi, n_segs, qlens, seqs, b);
-	a = collect_seed_hits(opt, opt->mid_occ, mi, qname, qlen_sum, &n_a, &rep_len, b);
-	radix_sort_128x(a, a + n_a);
+	collect_minimizers(b->km, opt, mi, n_segs, qlens, seqs, &mv);
+	if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
+	else a = collect_seed_hits(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);

 	if (mm_dbg_flag & MM_DBG_PRINT_SEED) {
 		fprintf(stderr, "RS\t%d\n", rep_len);
@@ -327,8 +320,9 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		if (rechain) { // redo chaining with a higher max_occ threshold
 			kfree(b->km, a);
 			kfree(b->km, u);
-			a = collect_seed_hits(opt, opt->max_occ, mi, qname, qlen_sum, &n_a, &rep_len, b);
-			radix_sort_128x(a, a + n_a);
+			kfree(b->km, mini_pos);
+			if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
+			else a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 			a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 		}
 	}
@@ -342,9 +336,10 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 						i == regs0[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));

 	chain_post(opt, max_chain_gap_ref, mi, b->km, qlen_sum, n_segs, qlens, &n_regs0, regs0, a);
+	if (!is_sr) mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos);

 	if (n_segs == 1) { // uni-segment
-		regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], quals? quals[0] : 0, &n_regs0, regs0, a);
+		regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], &n_regs0, regs0, a);
 		mm_set_mapq(n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr);
 		n_regs[0] = n_regs0, regs[0] = regs0;
 	} else { // multi-segment
@@ -353,7 +348,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		free(regs0);
 		for (i = 0; i < n_segs; ++i) {
 			mm_set_parent(b->km, opt->mask_level, n_regs[i], regs[i], opt->a * 2 + opt->b); // update mm_reg1_t::parent
-			regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], quals? quals[i] : 0, &n_regs[i], regs[i], seg[i].a);
+			regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a);
 			mm_set_mapq(n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr);
 		}
 		mm_seg_free(b->km, n_segs, seg);
@@ -361,14 +356,27 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 			mm_pair(b->km, max_chain_gap_ref, opt->pe_bonus, opt->a * 2 + opt->b, opt->a, qlens, n_regs, regs); // pairing
 	}

+	kfree(b->km, mv.a);
 	kfree(b->km, a);
 	kfree(b->km, u);
+	kfree(b->km, mini_pos);
+
+	if (b->km) {
+		km_stat(b->km, &kmst);
+		if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
+			fprintf(stderr, "QM\t%s\t%d\tcap=%ld,nCore=%ld,largest=%ld\n", qname, qlen_sum, kmst.capacity, kmst.n_cores, kmst.largest);
+		assert(kmst.n_blocks == kmst.n_cores); // otherwise, there is a memory leak
+		if (kmst.largest > 1U<<28) {
+			km_destroy(b->km);
+			b->km = km_init();
+		}
+	}
 }

 mm_reg1_t *mm_map(const mm_idx_t *mi, int qlen, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
 {
 	mm_reg1_t *regs;
-	mm_map_frag(mi, 1, &qlen, &seq, 0, n_regs, &regs, b, opt, qname);
+	mm_map_frag(mi, 1, &qlen, &seq, n_regs, &regs, b, opt, qname);
 	return regs;
 }

@@ -396,26 +404,23 @@ typedef struct {
 static void worker_for(void *_data, long i, int tid) // kt_for() callback
 {
    step_t *s = (step_t*)_data;
-	int *qlens, j, off = s->seg_off[i], pe_ori = s->p->opt->pe_ori, is_sr = !!(s->p->opt->flag & MM_F_SR);
-	const char **qseqs, **quals = 0;
+	int qlens[MM_MAX_SEG], j, off = s->seg_off[i], pe_ori = s->p->opt->pe_ori;
+	const char *qseqs[MM_MAX_SEG];
 	mm_tbuf_t *b = s->buf[tid];
+	assert(s->n_seg[i] <= MM_MAX_SEG);
 	if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
-		fprintf(stderr, "QR\t%s\t%d\n", s->seq[off].name, tid);
-	qlens = (int*)kmalloc(b->km, s->n_seg[i] * sizeof(int));
-	qseqs = (const char**)kmalloc(b->km, s->n_seg[i] * sizeof(const char**));
-	quals = (const char**)kmalloc(b->km, s->n_seg[i] * sizeof(const char**));
+		fprintf(stderr, "QR\t%s\t%d\t%d\n", s->seq[off].name, tid, s->seq[off].l_seq);
 	for (j = 0; j < s->n_seg[i]; ++j) {
 		if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1))))
 			mm_revcomp_bseq(&s->seq[off + j]);
 		qlens[j] = s->seq[off + j].l_seq;
 		qseqs[j] = s->seq[off + j].seq;
-		quals[j] = is_sr? s->seq[off + j].qual : 0;
 	}
 	if (s->p->opt->flag & MM_F_INDEPEND_SEG) {
 		for (j = 0; j < s->n_seg[i]; ++j)
-			mm_map_frag(s->p->mi, 1, &qlens[j], &qseqs[j], &quals[j], &s->n_reg[off+j], &s->reg[off+j], b, s->p->opt, s->seq[off+j].name);
+			mm_map_frag(s->p->mi, 1, &qlens[j], &qseqs[j], &s->n_reg[off+j], &s->reg[off+j], b, s->p->opt, s->seq[off+j].name);
 	} else {
-		mm_map_frag(s->p->mi, s->n_seg[i], qlens, qseqs, quals, &s->n_reg[off], &s->reg[off], b, s->p->opt, s->seq[off].name);
+		mm_map_frag(s->p->mi, s->n_seg[i], qlens, qseqs, &s->n_reg[off], &s->reg[off], b, s->p->opt, s->seq[off].name);
 	}
 	for (j = 0; j < s->n_seg[i]; ++j) // flip the query strand and coordinate to the original read strand
 		if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1)))) {
@@ -429,9 +434,6 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
 				r->rev = !r->rev;
 			}
 		}
-	kfree(b->km, qlens);
-	kfree(b->km, qseqs);
-	kfree(b->km, quals);
 }

 static void *worker_pipeline(void *shared, int step, void *in)
@@ -5,26 +5,34 @@
 #include <stdio.h>
 #include <sys/types.h>

-#define MM_F_NO_SELF     0x001
-#define MM_F_AVA         0x002
-#define MM_F_CIGAR       0x004
-#define MM_F_OUT_SAM     0x008
-#define MM_F_NO_QUAL     0x010
-#define MM_F_OUT_CG      0x020
-#define MM_F_OUT_CS      0x040
-#define MM_F_SPLICE      0x080 // splice mode
-#define MM_F_SPLICE_FOR  0x100 // match GT-AG
-#define MM_F_SPLICE_REV  0x200 // match CT-AC, the reverse complement of GT-AG
-#define MM_F_NO_LJOIN    0x400
-#define MM_F_OUT_CS_LONG 0x800
-#define MM_F_SR          0x1000
-#define MM_F_FRAG_MODE   0x2000
+#define MM_F_NO_DIAG       0x001 // no exact diagonal hit
+#define MM_F_NO_DUAL       0x002 // skip pairs where query name is lexicographically larger than target name
+#define MM_F_CIGAR         0x004
+#define MM_F_OUT_SAM       0x008
+#define MM_F_NO_QUAL       0x010
+#define MM_F_OUT_CG        0x020
+#define MM_F_OUT_CS        0x040
+#define MM_F_SPLICE        0x080 // splice mode
+#define MM_F_SPLICE_FOR    0x100 // match GT-AG
+#define MM_F_SPLICE_REV    0x200 // match CT-AC, the reverse complement of GT-AG
+#define MM_F_NO_LJOIN      0x400
+#define MM_F_OUT_CS_LONG   0x800
+#define MM_F_SR            0x1000
+#define MM_F_FRAG_MODE     0x2000
 #define MM_F_NO_PRINT_2ND  0x4000
 #define MM_F_2_IO_THREADS  0x8000
 #define MM_F_LONG_CIGAR    0x10000
 #define MM_F_INDEPEND_SEG  0x20000
 #define MM_F_SPLICE_FLANK  0x40000
 #define MM_F_SOFTCLIP      0x80000
+#define MM_F_FOR_ONLY      0x100000
+#define MM_F_REV_ONLY      0x200000
+#define MM_F_HEAP_SORT     0x400000
+#define MM_F_ALL_CHAINS    0x800000
+
+#define MM_I_HPC          0x1
+#define MM_I_NO_SEQ       0x2
+#define MM_I_NO_NAME      0x4

 #define MM_IDX_MAGIC   "MMI\2"

@@ -46,7 +54,7 @@ typedef struct {
 } mm_idx_seq_t;

 typedef struct {
-	int32_t b, w, k, is_hpc;
+	int32_t b, w, k, flag;
 	uint32_t n_seq;            // number of reference sequences
 	mm_idx_seq_t *seq;         // sequence name, length and offset
 	uint32_t *S;               // 4-bit packed sequence
@@ -64,23 +72,25 @@ typedef struct {
 } mm_extra_t;

 typedef struct {
-	int32_t id;                     // ID for internal uses (see also parent below)
-	uint32_t cnt:28, rev:1, seg_split:1, sam_pri:1, proper_frag:1; // number of minimizers; if on the reverse strand
-	uint32_t rid:31, inv:1;         // reference index; if this is an alignment from inversion rescue
-	int32_t score;                  // DP alignment score
-	int32_t qs, qe, rs, re;         // query start and end; reference start and end
-	int32_t parent, subsc;          // parent==id if primary; best alternate mapping score
-	int32_t as;                     // offset in the a[] array (for internal uses only)
-	int32_t mlen, blen;             // seeded exact match length; seeded alignment block length
-	uint32_t mapq:8, split:2, n_sub:22; // mapQ; split pattern; number of suboptimal mappings
-	uint32_t pe_thru:1, score0:31;
+	int32_t id;             // ID for internal uses (see also parent below)
+	int32_t cnt;            // number of minimizers; if on the reverse strand
+	int32_t rid;            // reference index; if this is an alignment from inversion rescue
+	int32_t score;          // DP alignment score
+	int32_t qs, qe, rs, re; // query start and end; reference start and end
+	int32_t parent, subsc;  // parent==id if primary; best alternate mapping score
+	int32_t as;             // offset in the a[] array (for internal uses only)
+	int32_t mlen, blen;     // seeded exact match length; seeded alignment block length
+	int32_t n_sub;          // number of suboptimal mappings
+	int32_t score0;         // initial chaining score (before chain merging/spliting)
+	uint32_t mapq:8, split:2, rev:1, inv:1, sam_pri:1, proper_frag:1, pe_thru:1, seg_split:1, seg_id:8, dummy:8;
 	uint32_t hash;
+	float div;
 	mm_extra_t *p;
 } mm_reg1_t;

 // indexing and mapping options
 typedef struct {
-	short k, w, is_hpc, bucket_bits;
+	short k, w, flag, bucket_bits;
 	int mini_batch_size;
 	uint64_t batch_size;
 } mm_idxopt_t;
@@ -110,6 +120,7 @@ typedef struct {
 	int end_bonus;
 	int min_dp_max;  // drop an alignment if the score of the max scoring segment is below this threshold
 	int min_ksw_len;
+	int anchor_ext_len, anchor_ext_shift;

 	int pe_ori, pe_bonus;

@@ -148,6 +159,7 @@ extern double mm_realtime0; // wall-clock timer
 * @return 0 if success; -1 if _present_ unknown
 */
 int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo);
+int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo);

 /**
 * Update mm_mapopt_t::mid_occ via mm_mapopt_t::mid_occ_frac
@@ -200,6 +212,21 @@ void mm_idx_reader_close(mm_idx_reader_t *r);

 int mm_idx_reader_eof(const mm_idx_reader_t *r);

+/**
+ * Create an index from strings in memory
+ *
+ * @param w            minimizer window size
+ * @param k            minimizer k-mer size
+ * @param is_hpc       use HPC k-mer if true
+ * @param bucket_bits  number of bits for the first level of the hash table
+ * @param n            number of sequences
+ * @param seq          sequences in A/C/G/T
+ * @param name         sequence names; could be NULL
+ *
+ * @return minimap2 index
+ */
+mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const char **seq, const char **name);
+
 /**
 * Print index statistics to stderr
 *
@@ -253,6 +280,8 @@ void mm_tbuf_destroy(mm_tbuf_t *b);
 */
 mm_reg1_t *mm_map(const mm_idx_t *mi, int l_seq, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *name);

+void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname);
+
 /**
 * Align a fasta/fastq file and print alignments to stdout
 *
@@ -269,7 +298,7 @@ int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_

 // deprecated APIs for backward compatibility
 void mm_mapopt_init(mm_mapopt_t *opt);
-mm_idx_t *mm_idx_build(const char *fn, int w, int k, int is_hpc, int n_threads);
+mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads);

 #ifdef __cplusplus
 }
@@ -1,4 +1,4 @@
-.TH minimap2 1 "11 November 2017" "minimap2-2.5 (r572)" "Bioinformatics tools"
+.TH minimap2 1 "1 February 2018" "minimap2-2.8 (r672)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -99,6 +99,14 @@ multiple times to map it against each batch of target sequences.
 may be ending with k/K/m/M/g/G. NB: mapping quality is incorrect given a
 multi-part index.
 .TP
+.B --idx-no-seq
+Don't store target sequences in the index. It saves disk space and memory but
+the index generated with this option will not work with
+.B -a
+or
+.BR -c .
+When base-level alignment is not requested, this option is automatically applied.
+.TP
 .BI -d \ FILE
 Save the minimizer index of
 .I target.fa
@@ -121,7 +129,7 @@ Ignore top
 fraction of most frequent minimizers [0.0002]
 .TP
 .BI -g \ INT
-Stop chain enlongation if there are no minimizers in
+Stop chain enlongation if there are no minimizers within
 .IR INT -bp
 [10000].
 .TP
@@ -140,11 +148,28 @@ Discard chains with chaining score
 [40]. Chaining score equals the approximate number of matching bases minus a
 concave gap penalty. It is computed with dynamic programming.
 .TP
+.B -D
+If query sequence name/length are identical to the target name/length, ignore
+diagonal anchors. This option also reduces DP-based extension along the
+diagonal.
+.TP
+.B -P
+Retain all chains and don't attempt to set primary chains. Options
+.B -p
+and
+.B -N
+have no effect when this option is in use.
+.TP
+.BR --dual = yes | no
+During chaining, whether to skip pairs wherein the query name is
+lexicographically greater than the target name [yes]
+.TP
 .B -X
-Perform all-vs-all mapping. In this mode, if the query sequence name is
-lexicographically larger than the target sequence name, the hits between them
-will be suppressed; if the query sequence name is the same as the target name,
-diagonal minimizer hits will also be suppressed.
+Equivalent to
+.RB ' -DP
+.BR --dual = no
+.BR --no-long-join '.
+Primarily used for all-vs-all read overlapping.
 .TP
 .BI -p \ FLOAT
 Minimal secondary-to-primary score ratio to output secondary mappings [0.8].
@@ -154,6 +179,9 @@ the chain with a lower score is secondary to the chain with a higher score.
 If the ratio of the scores is below
 .IR FLOAT ,
 the secondary chain will not be outputted or extended with DP alignment later.
+This option has no effect when
+.B -X
+is applied.
 .TP
 .BI -N \ INT
 Output at most
@@ -171,9 +199,14 @@ Increasing this option slows down spliced alignment. [200k]
 .TP
 .BI -F \ NUM
 Maximum fragment length (aka insert size; effective with
-.BR -xsr / --frag)
+.BR -xsr / --frag = yes )
 [800]
 .TP
+.BI -M \ FLOAT
+Mark as secondary a chain that overlaps with a better chain by
+.I FLOAT
+or more of the shorter chain [0.5]
+.TP
 .BI --max-chain-skip \ INT
 A heuristics that stops chaining early [50]. Minimap2 uses dynamic programming
 for chaining. The time complexity is quadratic in the number of seeds. This
@@ -196,8 +229,20 @@ applies a second round of chaining with a higher minimizer occurrence threshold
 if no good chain is found. In addition, minimap2 attempts to patch gaps between
 seeds with ungapped alignment.
 .TP
-.BR --frag [= no | yes ]
+.BR --frag = no | yes
 Whether to enable the fragment mode [no]
+.TP
+.B --for-only
+Only map to the forward strand of the reference sequences. For paired-end
+reads in the forward-reverse orientation, the first read is mapped to forward
+strand of the reference and the second read to the reverse stand.
+.TP
+.B --rev-only
+Only map to the reverse complement strand of the reference sequences.
+.TP
+.BR --heap-sort = no | yes
+If yes, sort anchors with heap merge, instead of radix sort. Heap merge is
+faster for short reads, but slower for long reads. [no]
 .SS Alignment options
 .TP 10
 .BI -A \ INT
@@ -247,14 +292,13 @@ no attempt to match GT-AG [n]
 .BI --end-bonus \ INT
 Score bonus when alignment extends to the end of the query sequence [0].
 .TP
-.BR --splice-flank [= yes | no ]
+.BR --splice-flank = yes | no
 Assume the next base to a
 .B GT
 donor site tends to be A/G (91% in human and 92% in mouse) and the preceding
 base to a
 .B AG
-acceptor tends to be C/T [yes with
-.BR --splice ].
+acceptor tends to be C/T [no].
 This trend is evolutionarily conservative, all the way to S. cerevisiae
 (PMID:18688272). Specifying this option generally leads to higher junction
 accuracy by several percents, so it is applied by default with
@@ -264,6 +308,18 @@ However, the SIRV control does not honor this trend
 on SIRV data, please add
 .B --splice-flank=no
 to the command line.
+.TP
+.BI --end-seed-pen \ INT
+Drop a terminal anchor if
+.IR s <log( g )+ INT ,
+where
+.I s
+is the local alignment score around the anchor and
+.I g
+the length of the terminal gap in the chain. This option is only effective
+with
+.BR --splice .
+It helps to avoid tiny terminal exons. [6]
 .SS Input/output options
 .TP 10
 .B -a
@@ -332,7 +388,7 @@ K/M/G/k/m/g suffix is accepted. A large
 helps load balancing in the multi-threading mode, at the cost of increased
 memory.
 .TP
-.BR --secondary [= yes | no ]
+.BR --secondary = yes | no
 Whether to output secondary alignments [yes]
 .TP
 .B --version
@@ -379,13 +435,13 @@ Up to 10% sequence divergence.
 .B ava-pb
 PacBio all-vs-all overlap mapping
 .RB ( -Hk19
-.B -w5 -Xp0 -m100 -g10000 --max-chain-skip
+.B -Xw5 -m100 -g10000 --max-chain-skip
 .BR 25 ).
 .TP
 .B ava-ont
 Oxford Nanopore all-vs-all overlap mapping
 .RB ( -k15
-.B -w5 -Xp0 -m100 -g10000 --max-chain-skip
+.B -Xw5 -m100 -g10000 --max-chain-skip
 .BR 25 ).
 Similarly, the major difference from
 .B ava-pb
@@ -407,8 +463,8 @@ tag ignores introns to demote hits to pseudogenes.
 .B sr
 Short single-end reads without splicing
 .RB ( -k21
-.B -w11 --sr --frag -A2 -B8 -O12,32 -E2,1 -r50 -p.5 -N20 -f1000,5000 -n2 -m20
-.B -s40 -g200 -2K50m
+.B -w11 --sr --frag=yes -A2 -B8 -O12,32 -E2,1 -r50 -p.5 -N20 -f1000,5000 -n2 -m20
+.B -s40 -g200 -2K50m --heap-sort=yes
 .BR --secondary=no ).
 .RE
 .SS Miscellaneous options
@@ -461,7 +517,7 @@ cb | cb | cb
 r | c | l .
 Tag	Type	Description
 _
-tp	A	Type of aln: P/primary, S/secondary and I/inversion
+tp	A	Type of aln: P/primary, S/secondary and I,i/inversion
 cm	i	Number of minimizers on the chain
 s1	i	Chaining score
 s2	i	Chaining score of the best secondary chain
@@ -502,8 +558,8 @@ 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 requires SSE2 instructions to compile. It is possible to add
-non-SSE2 support, but it would make minimap2 slower by several times.
+Minimap2 requires SSE2 or NEON instructions to compile. It is possible to add
+non-SSE2/NEON support, but it would make minimap2 slower by several times.
 .SH SEE ALSO
 .PP
 miniasm(1), minimap(1), bwa(1).
@@ -1,4 +1,4 @@
-#include "minimap.h"
+#include "mmpriv.h"

 int mm_verbose = 1;
 int mm_dbg_flag = 0;
@@ -0,0 +1,258 @@
+/*******************************
+ * Command line option parsing *
+ *******************************/
+
+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;
+}
+
+/***********************
+ * Interval operations *
+ ***********************/
+
+Interval = {};
+
+Interval.sort = function(a)
+{
+	if (typeof a[0] == 'number')
+		a.sort(function(x, y) { return x - y });
+	else a.sort(function(x, y) { return x[0] != y[0]? x[0] - y[0] : x[1] - y[1] });
+}
+
+Interval.merge = function(a, sorted)
+{
+	if (typeof sorted == 'undefined') sorted = true;
+	if (!sorted) Interval.sort(a);
+	var k = 0;
+	for (var i = 1; i < a.length; ++i) {
+		if (a[k][1] >= a[i][0])
+			a[k][1] = a[k][1] > a[i][1]? a[k][1] : a[i][1];
+		else a[++k] = a[i].slice(0);
+	}
+	a.length = k + 1;
+}
+
+Interval.dedup = function(a, sorted)
+{
+	if (typeof sorted == 'undefined') sorted = true;
+	if (!sorted) Interval.sort(a);
+	var k = 0;
+	for (var i = 1; i < a.length; ++i)
+		if (a[k][0] != a[i][0] || a[k][1] != a[i][1])
+			a[++k] = a[i].slice(0);
+	a.length = k + 1;
+}
+
+Interval.index_end = function(a, sorted)
+{
+	if (a.length == 0) return;
+	if (typeof sorted == 'undefined') sorted = true;
+	if (!sorted) Interval.sort(a);
+	a[0].push(0);
+	var k = 0, k_en = a[0][1];
+	for (var i = 1; i < a.length; ++i) {
+		if (k_en <= a[i][0]) {
+			for (++k; k < i; ++k)
+				if (a[k][1] > a[i][0])
+					break;
+			k_en = a[k][1];
+		}
+		a[i].push(k);
+	}
+}
+
+Interval.find_intv = function(a, x)
+{
+	var left = -1, right = a.length;
+	if (typeof a[0] == 'number') {
+		while (right - left > 1) {
+			var mid = left + ((right - left) >> 1);
+			if (a[mid] > x) right = mid;
+			else if (a[mid] < x) left = mid;
+			else return mid;
+		}
+	} else {
+		while (right - left > 1) {
+			var mid = left + ((right - left) >> 1);
+			if (a[mid][0] > x) right = mid;
+			else if (a[mid][0] < x) left = mid;
+			else return mid;
+		}
+	}
+	return left;
+}
+
+Interval.find_ovlp = function(a, st, en)
+{
+	if (a.length == 0 || st >= en) return [];
+	var l = Interval.find_intv(a, st);
+	var k = l < 0? 0 : a[l][a[l].length - 1];
+	var b = [];
+	for (var i = k; i < a.length; ++i) {
+		if (a[i][0] >= en) break;
+		else if (st < a[i][1])
+			b.push(a[i]);
+	}
+	return b;
+}
+
+/*****************
+ * Main function *
+ *****************/
+
+function read_bed(fn, to_merge, to_dedup)
+{
+	var file = new File(fn);
+	var buf = new Bytes();
+	var h = {};
+	while (file.readline(buf) >= 0) {
+		var t = buf.toString().split("\t");
+		if (h[t[0]] == null)
+			h[t[0]] = [];
+		var bst = parseInt(t[1]);
+		var ben = parseInt(t[2]);
+		if (t.length >= 12 && /^\d+$/.test(t[9])) {
+			t[9] = parseInt(t[9]);
+			var sz = t[10].split(",");
+			var st = t[11].split(",");
+			for (var i = 0; i < t[9]; ++i) {
+				st[i] = parseInt(st[i]);
+				sz[i] = parseInt(sz[i]);
+				h[t[0]].push([bst + st[i], bst + st[i] + sz[i], 0, 0, 0]);
+			}
+		} else {
+			h[t[0]].push([bst, ben, 0, 0, 0]);
+		}
+	}
+	buf.destroy();
+	file.close();
+	for (var chr in h) {
+		if (to_merge) Interval.merge(h[chr], false);
+		else if (to_dedup) Interval.dedup(h[chr], false);
+		else Interval.sort(h[chr]);
+		Interval.index_end(h[chr]);
+	}
+	return h;
+}
+
+function main(args)
+{
+	var c, print_len = false, to_merge = true, to_dedup = false, fn_excl = null;
+	while ((c = getopt(args, "pde:")) != null) {
+		if (c == 'p') print_len = true;
+		else if (c == 'd') to_dedup = true, to_merge = false;
+		else if (c == 'e') fn_excl = getopt.arg;
+	}
+
+	if (args.length - getopt.ind < 2) {
+		print("Usage: k8 cnt-feat.js [options] <target.bed> <feature.bed>");
+		print("Options:");
+		print("  -e FILE    exclude features overlapping regions in BED FILE []");
+		print("  -p         print number of covered bases for each feature");
+		exit(1);
+	}
+
+	var excl = fn_excl != null? read_bed(fn_excl, true, false) : null;
+	var target = read_bed(args[getopt.ind], to_merge, to_dedup);
+
+	var file, buf = new Bytes();
+	var tot_len = 0, hit_len = 0;
+	file = args[getopt.ind+1] != '-'? new File(args[getopt.ind+1]) : new File();
+	while (file.readline(buf) >= 0) {
+		var t = buf.toString().split("\t");
+		var a = [];
+		var bst = parseInt(t[1]);
+		var ben = parseInt(t[2]);
+		if (t.length >= 12 && /^\d+$/.test(t[9])) { // BED12
+			t[9] = parseInt(t[9]);
+			var sz = t[10].split(",");
+			var st = t[11].split(",");
+			for (var i = 0; i < t[9]; ++i) {
+				st[i] = parseInt(st[i]);
+				sz[i] = parseInt(sz[i]);
+				a.push([bst + st[i], bst + st[i] + sz[i], false]);
+			}
+		} else a.push([bst, ben, false]); // 3-column BED
+		var feat_len = 0;
+		for (var i = 0; i < a.length; ++i) {
+			if (excl != null && excl[t[0]] != null) {
+				var oe = Interval.find_ovlp(excl[t[0]], a[i][0], a[i][1]);
+				if (oe.length > 0)
+					continue;
+			}
+			a[i][2] = true;
+			feat_len += a[i][1] - a[i][0];
+		}
+		tot_len += feat_len;
+		if (target[t[0]] == null) continue;
+		var b = [];
+		for (var i = 0; i < a.length; ++i) {
+			if (!a[i][2]) continue;
+			var o = Interval.find_ovlp(target[t[0]], a[i][0], a[i][1]);
+			for (var j = 0; j < o.length; ++j) {
+				var max_st = o[j][0] > a[i][0]? o[j][0] : a[i][0];
+				var min_en = o[j][1] < a[i][1]? o[j][1] : a[i][1];
+				b.push([max_st, min_en]);
+				o[j][2] += min_en - max_st;
+				++o[j][3];
+				if (max_st == o[j][0] && min_en == o[j][1])
+					++o[j][4];
+			}
+		}
+		// find the length covered
+		var feat_hit_len = 0;
+		if (b.length > 0) {
+			b.sort(function(a,b) {return a[0]-b[0]});
+			var st = b[0][0], en = b[0][1];
+			for (var i = 1; i < b.length; ++i) {
+				if (b[i][0] <= en) en = en > b[i][1]? en : b[i][1];
+				else feat_hit_len += en - st, st = b[i][0], en = b[i][1];
+			}
+			feat_hit_len += en - st;
+		}
+		hit_len += feat_hit_len;
+		if (print_len) print('F', t.slice(0, 4).join("\t"), feat_len, feat_hit_len);
+	}
+	file.close();
+
+	buf.destroy();
+
+	warn("# feature bases: " + tot_len);
+	warn("# feature bases overlapping targets: " + hit_len + ' (' + (100.0 * hit_len / tot_len).toFixed(2) + '%)');
+}
+
+main(arguments);
@@ -0,0 +1,150 @@
+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, fn_ucsc_fai = null, is_short = false;
+while ((c = getopt(arguments, "u:s")) != null) {
+	if (c == 'u') fn_ucsc_fai = getopt.arg;
+	else if (c == 's') is_short = true;
+}
+
+if (getopt.ind == arguments.length) {
+	print("Usage: k8 gff2bed.js [-u ucsc-genome.fa.fai] <in.gff>");
+	exit(1);
+}
+
+var ens2ucsc = {};
+if (fn_ucsc_fai != null) {
+	var buf = new Bytes();
+	var file = new File(fn_ucsc_fai);
+	while (file.readline(buf) >= 0) {
+		var t = buf.toString().split("\t");
+		var s = t[0];
+		if (/_(random|alt|decoy)$/.test(s)) {
+			s = s.replace(/_(random|alt|decoy)$/, '');
+			s = s.replace(/^chr\S+_/, '');
+		} else {
+			s = s.replace(/^chrUn_/, '');
+		}
+		s = s.replace(/v(\d+)/, ".$1");
+		if (s != t[0]) ens2ucsc[s] = t[0];
+	}
+	file.close();
+	buf.destroy();
+}
+
+var colors = {
+		'protein_coding':'0,128,255',
+		'lincRNA':'0,192,0',
+		'snRNA':'0,192,0',
+		'miRNA':'0,192,0',
+		'misc_RNA':'0,192,0'
+	};
+
+function print_bed12(exons, cds_st, cds_en, is_short)
+{
+	if (exons.length == 0) return;
+	var name = is_short? exons[0][7] + "|" + exons[0][5] : exons[0].slice(4, 7).join("|");
+	var a = exons.sort(function(a,b) {return a[1]-b[1]});
+	var sizes = [], starts = [], st, en;
+	st = a[0][1];
+	en = a[a.length - 1][2];
+	if (cds_st == 1<<30) cds_st = st;
+	if (cds_en == 0) cds_en = en;
+	if (cds_st < st || cds_en > en)
+		throw Error("inconsistent thick start or end for transcript " + a[0][4]);
+	for (var i = 0; i < a.length; ++i) {
+		sizes.push(a[i][2] - a[i][1]);
+		starts.push(a[i][1] - st);
+	}
+	var color = colors[a[0][5]];
+	if (color == null) color = '196,196,196';
+	print(a[0][0], st, en, name, 1000, a[0][3], cds_st, cds_en, color, a.length, sizes.join(",") + ",", starts.join(",") + ",");
+}
+
+var re_gtf = /(transcript_id|transcript_type|transcript_biotype|gene_name|transcript_name) "([^"]+)";/g;
+var re_gff3 = /(transcript_id|transcript_type|transcript_biotype|gene_name|transcript_name)=([^;]+)/g;
+var buf = new Bytes();
+var file = new File(arguments[getopt.ind]);
+
+var exons = [], cds_st = 1<<30, cds_en = 0, last_id = null;
+while (file.readline(buf) >= 0) {
+	var t = buf.toString().split("\t");
+	if (t[0].charAt(0) == '#') continue;
+	if (t[2] != "CDS" && t[2] != "exon") continue;
+	t[3] = parseInt(t[3]) - 1;
+	t[4] = parseInt(t[4]);
+	var id = null, type = "", gname = "N/A", biotype = "", m, tname = "N/A";
+	while ((m = re_gtf.exec(t[8])) != null) {
+		if (m[1] == "transcript_id") id = m[2];
+		else if (m[1] == "transcript_type") type = m[2];
+		else if (m[1] == "transcript_biotype") biotype = m[2];
+		else if (m[1] == "gene_name") name = m[2];
+		else if (m[1] == "transcript_name") tname = m[2];
+	}
+	while ((m = re_gff3.exec(t[8])) != null) {
+		if (m[1] == "transcript_id") id = m[2];
+		else if (m[1] == "transcript_type") type = m[2];
+		else if (m[1] == "transcript_biotype") biotype = m[2];
+		else if (m[1] == "gene_name") name = m[2];
+		else if (m[1] == "transcript_name") tname = m[2];
+	}
+	if (type == "" && biotype != "") type = biotype;
+	if (id == null) throw Error("No transcript_id");
+	if (id != last_id) {
+		print_bed12(exons, cds_st, cds_en, is_short);
+		exons = [], cds_st = 1<<30, cds_en = 0;
+		last_id = id;
+	}
+	if (t[2] == "CDS") {
+		cds_st = cds_st < t[3]? cds_st : t[3];
+		cds_en = cds_en > t[4]? cds_en : t[4];
+	} else if (t[2] == "exon") {
+		if (fn_ucsc_fai != null) {
+			if (ens2ucsc[t[0]] != null)
+				t[0] = ens2ucsc[t[0]];
+			else if (/^[A-Z]+\d+\.\d+$/.test(t[0]))
+				t[0] = t[0].replace(/([A-Z]+\d+)\.(\d+)/, "chrUn_$1v$2");
+		}
+		exons.push([t[0], t[3], t[4], t[6], id, type, name, tname]);
+	}
+}
+if (last_id != null)
+	print_bed12(exons, cds_st, cds_en, is_short);
+
+file.close();
+buf.destroy();
@@ -161,7 +161,8 @@ for (var tid in tr) {
 	if (anno[t[0]] == null) anno[t[0]] = [];
 	var s = t[4];
 	for (var i = 0; i < s.length - 1; ++i) {
-		if (s[i][1] >= s[i+1][0]) throw Error("ERROR: wrong annotation!");
+		if (s[i][1] >= s[i+1][0])
+			warn("WARNING: incorrect annotation for transcript "+tid+" ("+s[i][1]+" >= "+s[i+1][0]+")")
 		anno[t[0]].push([s[i][1], s[i+1][0]]);
 	}
 }
@@ -0,0 +1,105 @@
+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, min_ovlp = 2000, min_frac = 0.95, min_mapq = 10;
+while ((c = getopt(arguments, "q:l:f:")) != null) {
+	if (c == 'q') min_mapq = parseInt(getopt.arg);
+	else if (c == 'l') min_ovlp = parseInt(getopt.arg);
+	else if (c == 'f') min_frac = parseFloat(getopt.arg);
+}
+if (arguments.length - getopt.ind < 2) {
+	print("Usage: sort -k6,6 -k8,8n to-ref.paf | k8 ov-eval.js [options] - <ovlp.paf>");
+	print("Options:");
+	print("  -l INT     min overlap length [2000]");
+	print("  -q INT     min mapping quality [10]");
+	print("  -f FLOAT   min fraction of mapped length [0.95]");
+	exit(1);
+}
+
+var buf = new Bytes();
+var file = arguments[getopt.ind] == '-'? new File() : new File(arguments[getopt.ind]);
+var a = [], h = {};
+while (file.readline(buf) >= 0) {
+	var t = buf.toString().split("\t");
+	var is_pri = false;
+	if (parseInt(t[11]) < min_mapq) continue;
+	for (var i = 12; i < t.length; ++i)
+		if (t[i] == 'tp:A:P')
+			is_pri = true;
+	if (!is_pri) continue;
+	for (var i = 1; i <= 3; ++i)
+		t[i] = parseInt(t[i]);
+	for (var i = 6; i <= 8; ++i)
+		t[i] = parseInt(t[i]);
+	if (t[3] - t[2] < min_ovlp || t[8] - t[7] < min_ovlp || (t[3] - t[2]) / t[1] < min_frac)
+		continue;
+	var ctg = t[5], st = t[7], en = t[8];
+	while (a.length > 0) {
+		if (a[0][0] == ctg && a[0][2] > st)
+			break;
+		else a.shift();
+	}
+	for (var j = 0; j < a.length; ++j) {
+		if (a[j][3] == t[0]) continue;
+		var len = (en > a[j][2]? a[j][2] : en) - st;
+		if (len >= min_ovlp) {
+			var key = a[j][3] < t[0]? a[j][3] + "\t" + t[0] : t[0] + "\t" + a[j][3];
+			h[key] = len;
+		}
+	}
+	a.push([ctg, st, en, t[0]]);
+}
+file.close();
+
+file = new File(arguments[getopt.ind + 1]);
+while (file.readline(buf) >= 0) {
+	var t = buf.toString().split("\t");
+	var key = t[0] < t[5]? t[0] + "\t" + t[5] : t[5] + "\t" + t[0];
+	if (h[key] > 0) h[key] = -h[key];
+}
+file.close();
+buf.destroy();
+
+var n_ovlp = 0, n_missing = 0;
+for (var key in h) {
+	++n_ovlp;
+	if (h[key] > 0) ++n_missing;
+}
+print(n_ovlp + " overlaps inferred from the reference mapping");
+print(n_missing + " missed by the read overlapper");
+print((100 * (1 - n_missing / n_ovlp)).toFixed(2) + "% sensitivity");
@@ -36,20 +36,21 @@ var getopt = function(args, ostr) {
 	return optopt;
 }

-var c, maf_out = false, line_len = 80;
-while ((c = getopt(arguments, "ml:")) != null) {
-	if (c == 'm') maf_out = true;
-	else if (c == 'l') line_len = parseInt(getopt.arg); // TODO: not implemented yet
+var c, line_len = 80, fmt = "aln";
+while ((c = getopt(arguments, "f:l:")) != null) {
+	if (c == 'f') {
+		fmt = getopt.arg;
+		if (fmt != "aln" && fmt != "lastz-cigar" && fmt != "maf")
+			throw Error("format must be one of aln, lastz-cigar and maf");
+	} else if (c == 'l') line_len = parseInt(getopt.arg);
 }
 if (line_len == 0) line_len = 0x7fffffff;

 if (getopt.ind == arguments.length) {
-	print("Usage: k8 paf2aln.js [options] <with-cs.paf>");
+	print("Usage: k8 paf2aln.js [options] <in.paf>");
 	print("Options:");
-	print("  -m        MAF output (BLAST-like output by default)");
+	print("  -f STR    output format: aln (BLAST-like), maf or lastz-cigar [aln]");
 	print("  -l INT    line length in BLAST-like output [80]");
-	print("");
-	print("Note: this script only works when minimap2 is run with option '-S'");
 	exit(1);
 }

@@ -66,7 +67,7 @@ function padding_str(x, len, right)
 function update_aln(s_ref, s_qry, s_mid, type, seq, slen)
 {
 	var l = type == '*'? 1 : seq.length;
-	if (type == '=') {
+	if (type == '=' || type == ':') {
 		s_ref.set(seq);
 		s_qry.set(seq);
 		s_mid.set(Array(l+1).join("|"));
@@ -96,49 +97,73 @@ function print_aln(rs, qs, strand, slen, elen, s_ref, s_qry, s_mid)
 	var st, en;
 	if (strand == '+') st = qs + slen[1] + 1, en = qs + elen[1];
 	else st = qs - slen[1], en = qs - elen[1] + 1;
-	print(["Qry" + strand + ":", padding_str(st,               10, false), s_qry.toString(), padding_str(en          , 10, true)].join(" "));
+	print(["Qry" + strand + ":", padding_str(st, 10, false), s_qry.toString(), padding_str(en, 10, true)].join(" "));
 }

-var s_ref = new Bytes(), s_qry = new Bytes(), s_mid = new Bytes();
-var re = /([=\-\+\*])([A-Za-z]+)/g;
+var s_ref = new Bytes(), s_qry = new Bytes(), s_mid = new Bytes(); // these are used to show padded alignment
+var re_cs = /([:=\-\+\*])(\d+|[A-Za-z]+)/g;
+var re_cg = /(\d+)([MIDNSH])/g;

 var buf = new Bytes();
-var file = new File(arguments[getopt.ind]);
-if (maf_out) print("##maf version=1\n");
+var file = arguments[getopt.ind] == "-"? new File() : new File(arguments[getopt.ind]);
+var lineno = 0;
+if (fmt == "maf") print("##maf version=1\n");
 while (file.readline(buf) >= 0) {
 	var m, line = buf.toString();
 	var t = line.split("\t", 12);
-	if ((m = /\tcs:Z:(\S+)/.exec(line)) == null) continue;
-	var cs = m[1];
+	++lineno;
 	s_ref.length = s_qry.length = s_mid.length = 0;
 	var slen = [0, 0], elen = [0, 0];
-	if (maf_out) {
-		while ((m = re.exec(cs)) != null)
-			update_aln(s_ref, s_qry, s_mid, m[1], m[2], elen);
-		if (maf_out) {
-			var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : 0;
-			var len = t[0].length > t[5].length? t[0].length : t[5].length;
-			print("a " + score);
-			print(["s", padding_str(t[5], len, true), padding_str(t[7], 10, false), padding_str(parseInt(t[8]) - parseInt(t[7]), 10, false),
-				   "+", padding_str(t[6], 10, false), s_ref.toString()].join(" "));
-			var qs, qe, ql = parseInt(t[1]);
-			if (t[4] == '+') {
-				qs = parseInt(t[2]);
-				qe = parseInt(t[3]);
-			} else {
-				qs = ql - parseInt(t[3]);
-				qe = ql - parseInt(t[2]);
-			}
-			print(["s", padding_str(t[0], len, true), padding_str(qs, 10, false), padding_str(qe - qs, 10, false),
-				   t[4], padding_str(ql, 10, false), s_qry.toString()].join(" "));
-			print("");
+	if (fmt == "lastz-cigar") { // LASTZ-cigar output
+		var cg = (m = /\tcg:Z:(\S+)/.exec(line)) != null? m[1] : null;
+		if (cg == null) {
+			warn("WARNING: converting to LASTZ-cigar format requires the 'cg' tag, which is absent on line " + lineno);
+			continue;
 		}
-	} else {
-		line = line.replace(/\tc[sg]:Z:\S+/g, "");
+		var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? m[1] : 0;
+		var out = ['cigar:', t[0], t[2], t[3], t[4], t[5], t[7], t[8], '+', score];
+		while ((m = re_cg.exec(cg)) != null)
+			out.push(m[2], m[1]);
+		print(out.join(" "));
+	} else if (fmt == "maf") { // MAF output
+		var cs = (m = /\tcs:Z:(\S+)/.exec(line)) != null? m[1] : null;
+		if (cs == null) {
+			warn("WARNING: converting to MAF requires the 'cs' tag, which is absent on line " + lineno);
+			continue;
+		}
+		while ((m = re_cs.exec(cs)) != null) {
+			if (m[1] == ':')
+				throw Error("converting to MAF only works with 'minimap2 --cs=long'");
+			update_aln(s_ref, s_qry, s_mid, m[1], m[2], elen);
+		}
+		var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : 0;
+		var len = t[0].length > t[5].length? t[0].length : t[5].length;
+		print("a " + score);
+		print(["s", padding_str(t[5], len, true), padding_str(t[7], 10, false), padding_str(parseInt(t[8]) - parseInt(t[7]), 10, false),
+				"+", padding_str(t[6], 10, false), s_ref.toString()].join(" "));
+		var qs, qe, ql = parseInt(t[1]);
+		if (t[4] == '+') {
+			qs = parseInt(t[2]);
+			qe = parseInt(t[3]);
+		} else {
+			qs = ql - parseInt(t[3]);
+			qe = ql - parseInt(t[2]);
+		}
+		print(["s", padding_str(t[0], len, true), padding_str(qs, 10, false), padding_str(qe - qs, 10, false),
+				t[4], padding_str(ql, 10, false), s_qry.toString()].join(" "));
+		print("");
+	} else { // BLAST-like output
+		var cs = (m = /\tcs:Z:(\S+)/.exec(line)) != null? m[1] : null;
+		if (cs == null) {
+			warn("WARNING: converting to BLAST-like alignment requires the 'cs' tag, which is absent on line " + lineno);
+			continue;
+		}
+		line = line.replace(/\tc[sg]:Z:\S+/g, ""); // get rid of cs or cg tags
 		print('>' + line);
 		var rs = parseInt(t[7]), qs = t[4] == '+'? parseInt(t[2]) : parseInt(t[3]);
 		var n_blocks = 0;
-		while ((m = re.exec(cs)) != null) {
+		while ((m = re_cs.exec(cs)) != null) {
+			if (m[1] == ':') m[2] = Array(parseInt(m[2]) + 1).join("=");
 			var start = 0, rest = m[1] == '*'? 1 : m[2].length;
 			while (rest > 0) {
 				var l_proc;
@@ -0,0 +1,188 @@
+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 re_cs = /([:=*+-])(\d+|[A-Za-z]+)/g;
+var c, min_cov_len = 10000, min_var_len = 50000, gap_thres = 50, min_mapq = 5;
+while ((c = getopt(arguments, "l:L:g:q:")) != null) {
+	if (c == 'l') min_cov_len = parseInt(getopt.arg);
+	else if (c == 'L') min_var_len = parseInt(getopt.arg);
+	else if (c == 'g') gap_thres = parseInt(getopt.arg);
+	else if (c == 'q') min_mapq = parseInt(getopt.arg);
+}
+
+if (arguments.length == getopt.ind) {
+	print("Usage: k8 paf2diff.js [options] <with-cs.paf>");
+	print("Options:");
+	print("  -l INT    min alignment length to compute coverage ["+min_cov_len+"]");
+	print("  -L INT    min alignment length to call variants ["+min_var_len+"]");
+	print("  -q INT    min mapping quality ["+min_mapq+"]");
+	print("  -g INT    short/long gap threshold (for statistics only) ["+gap_thres+"]");
+	exit(1);
+}
+
+var file = arguments[getopt.ind] == '-'? new File() : new File(arguments[getopt.ind]);
+var buf = new Bytes();
+var tot_len = 0, n_sub = [0, 0, 0], n_ins = [0, 0, 0, 0], n_del = [0, 0, 0, 0];
+
+function count_var(o)
+{
+	if (o[3] > 1) return;
+	if (o[5] == '-' && o[6] == '-') return;
+	if (o[5] == '-') { // insertion
+		var l = o[6].length;
+		if (l == 1) ++n_ins[0];
+		else if (l == 2) ++n_ins[1];
+		else if (l < gap_thres) ++n_ins[2];
+		else ++n_ins[3];
+	} else if (o[6] == '-') { // deletion
+		var l = o[5].length;
+		if (l == 1) ++n_del[0];
+		else if (l == 2) ++n_del[1];
+		else if (l < gap_thres) ++n_del[2];
+		else ++n_del[3];
+	} else {
+		++n_sub[0];
+		var s = o[5] + o[6];
+		if (s == 'ag' || s == 'ga' || s == 'ct' || s == 'tc')
+			++n_sub[1];
+		else ++n_sub[2];
+	}
+}
+
+var a = [], out = [];
+var c1_ctg = null, c1_start = 0, c1_end = 0, c1_counted = false, c1_len = 0;
+while (file.readline(buf) >= 0) {
+	var line = buf.toString();
+	if (!/\ts2:i:/.test(line)) continue; // skip secondary alignments
+	var m, t = line.split("\t", 12);
+	for (var i = 6; i <= 11; ++i)
+		t[i] = parseInt(t[i]);
+	if (t[10] < min_cov_len || t[11] < min_mapq) continue;
+	var ctg = t[5], x = t[7], end = t[8];
+	// compute regions covered by 1 contig
+	if (ctg != c1_ctg || x >= c1_end) {
+		if (c1_counted && c1_end > c1_start) {
+			c1_len += c1_end - c1_start;
+			print('R', c1_ctg, c1_start, c1_end);
+		}
+		c1_ctg = ctg, c1_start = x, c1_end = end;
+		c1_counted = (t[10] >= min_var_len);
+	} else if (end > c1_end) { // overlap
+		if (c1_counted && x > c1_start) {
+			c1_len += x - c1_start;
+			print('R', c1_ctg, c1_start, x);
+		}
+		c1_start = c1_end, c1_end = end;
+		c1_counted = (t[10] >= min_var_len);
+	} else { // contained
+		if (c1_counted && x > c1_start) {
+			c1_len += x - c1_start;
+			print('R', c1_ctg, c1_start, x);
+		}
+		c1_start = end;
+	}
+	// output variants ahead of this alignment
+	while (out.length) {
+		if (out[0][0] != ctg || out[0][2] <= x) {
+			count_var(out[0]);
+			print('V', out[0].join("\t"));
+			out.shift();
+		} else break;
+	}
+	// update coverage
+	for (var i = 0; i < out.length; ++i)
+		if (out[i][1] >= x && out[i][2] <= end)
+			++out[i][3];
+	// drop alignments that don't overlap with the current one
+	var k = 0;
+	for (var i = 0; i < a.length; ++i)
+		if (a[0][0] == ctg && a[0][2] > x)
+			a[k++] = a[i];
+	a.length = k;
+	// core loop
+	if (t[10] >= min_var_len) {
+		if ((m = /\tcs:Z:(\S+)/.exec(line)) == null) continue; // no cs tag
+		var cs = m[1];
+		var blen = 0, n_diff = 0;
+		tot_len += t[10];
+		while ((m = re_cs.exec(cs)) != null) {
+			var cov = 1;
+			if (m[1] == '*' || m[1] == '+' || m[1] == '-')
+				for (var i = 0; i < a.length; ++i)
+					if (a[0][2] > x) ++cov;
+			if (m[1] == '=' || m[1] == ':') {
+				var l = m[1] == '='? m[2].length : parseInt(m[2]);
+				x += l, blen += l;
+			} else if (m[1] == '*') {
+				out.push([t[5], x, x+1, cov, t[11], m[2].charAt(0), m[2].charAt(1)]);
+				++x, ++blen, ++n_diff;
+			} else if (m[1] == '+') {
+				out.push([t[5], x, x, cov, t[11], '-', m[2]]);
+				++blen, ++n_diff;
+			} else if (m[1] == '-') {
+				out.push([t[5], x, x + m[2].length, cov, t[11], m[2], '-']);
+				x += m[2].length, ++blen, ++n_diff;
+			}
+		}
+	}
+	a.push([t[5], t[7], t[8]]);
+}
+if (c1_counted && c1_end > c1_start) {
+	c1_len += c1_end - c1_start;
+	print('R', c1_ctg, c1_start, c1_end);
+}
+while (out.length) {
+	count_var(out[0]);
+	print('V', out[0].join("\t"));
+	out.shift();
+}
+
+//warn(tot_len + " alignment columns considered in calling");
+warn(c1_len + " reference bases covered by exactly one contig");
+warn(n_sub[0] + " substitutions; ts/tv = " + (n_sub[1]/n_sub[2]).toFixed(3));
+warn(n_del[0] + " 1bp deletions");
+warn(n_ins[0] + " 1bp insertions");
+warn(n_del[1] + " 2bp deletions");
+warn(n_ins[1] + " 2bp insertions");
+warn(n_del[2] + " [3,"+gap_thres+") deletions");
+warn(n_ins[2] + " [3,"+gap_thres+") insertions");
+warn(n_del[3] + " >="+gap_thres+" deletions");
+warn(n_ins[3] + " >="+gap_thres+" insertions");
+
+buf.destroy();
+file.close();
@@ -0,0 +1,148 @@
+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 colors = ["0,128,255", "255,0,0", "0,192,0"];
+
+function print_lines(a, fmt) {
+	if (a.length == 0) return;
+	if (fmt == "bed") {
+		var n_pri = 0;
+		for (var i = 0; i < a.length; ++i)
+			if (a[i][8] == 0) ++n_pri;
+		if (n_pri > 1) {
+			for (var i = 0; i < a.length; ++i)
+				if (a[i][8] == 0) a[i][8] = 1;
+		} else if (n_pri == 0) {
+			warn("Warning: " + a[0][3] + " doesn't have a primary alignment");
+		}
+		for (var i = 0; i < a.length; ++i) {
+			a[i][8] = colors[a[i][8]];
+			print(a[i].join("\t"));
+		}
+	}
+	a.length = 0;
+}
+
+function main(args) {
+	var re = /(\d+)([MIDNSH])/g;
+	var c, fmt = "bed", fn_name_conv = null;
+	while ((c = getopt(args, "f:n:")) != null) {
+		if (c == 'f') fmt = getopt.arg;
+		else if (c == 'n') fn_name_conv = getopt.arg;
+	}
+	if (getopt.ind == args.length) {
+		warn("Usage: k8 splice2bed.js <in.paf>");
+		exit(1);
+	}
+
+	var conv = null;
+	if (fn_name_conv != null) {
+		conv = new Map();
+		var file = new File(fn_name_conv);
+		var buf = new Bytes();
+		while (file.readline(buf) >= 0) {
+			var t = buf.toString().split("\t");
+			conv.put(t[0], t[1]);
+		}
+		buf.destroy();
+		file.close();
+	}
+
+	var file = new File(args[getopt.ind]);
+	var buf = new Bytes();
+	var a = [];
+	while (file.readline(buf) >= 0) {
+		var line = buf.toString();
+		if (line.charAt(0) == '@') continue; // skip SAM header lines
+		var t = line.split("\t");
+		var is_pri = false, cigar = null, a1;
+		var qname = conv != null? conv.get(t[0]) : null;
+		if (qname != null) t[0] = qname;
+		if (t.length >= 10 && t[4] != '+' && t[4] != '-' && /^\d+/.test(t[1])) { // SAM
+			var flag = parseInt(t[1]);
+			if (flag&1) t[0] += '/' + (flag>>6&3);
+		}
+		if (a.length && a[0][3] != t[0]) {
+			print_lines(a, fmt);
+			a = [];
+		}
+		if (t.length >= 12 && (t[4] == '+' || t[4] == '-')) { // PAF
+			for (var i = 12; i < t.length; ++i) {
+				if (t[i].substr(0, 5) == 'cg:Z:') {
+					cigar = t[i].substr(5);
+				} else if (t[i].substr(0, 5) == 's2:i:') {
+					is_pri = true;
+				}
+			}
+			a1 = [t[5], t[7], t[8], t[0], Math.floor(t[9]/t[10]*1000), t[4]];
+		} else if (t.length >= 10) { // SAM
+			var flag = parseInt(t[1]);
+			if ((flag&4) || a[2] == '*') continue;
+			cigar = t[5];
+			is_pri = (flag&0x100)? false : true;
+			a1 = [t[2], parseInt(t[3])-1, null, t[0], 1000, (flag&16)? '-' : '+'];
+		} else {
+			throw Error("unrecognized input format");
+		}
+		if (cigar == null) throw Error("missing CIGAR");
+		var m, x0 = 0, x = 0, bs = [], bl = [];
+		while ((m = re.exec(cigar)) != null) {
+			if (m[2] == 'M' || m[2] == 'D') {
+				x += parseInt(m[1]);
+			} else if (m[2] == 'N') {
+				bs.push(x0);
+				bl.push(x - x0);
+				x += parseInt(m[1]);
+				x0 = x;
+			}
+		}
+		bs.push(x0);
+		bl.push(x - x0);
+		// write the BED12 line
+		if (a1[2] == null) a1[2] = a1[1] + x;
+		a1.push(a1[1], a1[2]); // thick start/end is the same as start/end
+		a1.push(is_pri? 0 : 2, bs.length, bl.join(",")+",", bs.join(",")+",");
+		a.push(a1);
+	}
+	print_lines(a, fmt);
+	buf.destroy();
+	file.close();
+	if (conv != null) conv.destroy();
+}
+
+main(arguments);
@@ -16,6 +16,7 @@
 #define MM_SEED_LONG_JOIN  (1ULL<<40)
 #define MM_SEED_IGNORE     (1ULL<<41)
 #define MM_SEED_TANDEM     (1ULL<<42)
+#define MM_SEED_SELF       (1ULL<<43)

 #define MM_SEED_SEG_SHIFT  48
 #define MM_SEED_SEG_MASK   (0xffULL<<(MM_SEED_SEG_SHIFT))
@@ -64,7 +65,7 @@ const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);
 int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f);
 mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, 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, const char *qual, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
+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);

 mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a);
 void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a);
@@ -79,6 +80,8 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_re
 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, int min_chain_sc, int match_sc, int rep_len, int is_sr);

+void mm_est_err(const mm_idx_t *mi, int qlen, int n_regs, mm_reg1_t *regs, const mm128_t *a, int32_t n, const uint64_t *mini_pos);
+
 mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int n_regs0, const mm_reg1_t *regs0, int *n_regs, mm_reg1_t **regs, const mm128_t *a);
 void mm_seg_free(void *km, int n_segs, mm_seg_t *segs);
 void mm_pair(void *km, int max_gap_ref, int dp_bonus, int sub_diff, int match_sc, const int *qlens, int *n_regs, mm_reg1_t **regs);
@@ -0,0 +1,140 @@
+#include <stdio.h>
+#include "mmpriv.h"
+
+void mm_mapopt_init(mm_mapopt_t *opt)
+{
+	memset(opt, 0, sizeof(mm_mapopt_t));
+	opt->seed = 11;
+	opt->mid_occ_frac = 2e-4f;
+	opt->sdust_thres = 0; // no SDUST masking
+
+	opt->min_cnt = 3;
+	opt->min_chain_score = 40;
+	opt->bw = 500;
+	opt->max_gap = 5000;
+	opt->max_gap_ref = -1;
+	opt->max_chain_skip = 25;
+
+	opt->mask_level = 0.5f;
+	opt->pri_ratio = 0.8f;
+	opt->best_n = 5;
+
+	opt->max_join_long = 20000;
+	opt->max_join_short = 2000;
+	opt->min_join_flank_sc = 1000;
+
+	opt->a = 2, opt->b = 4, opt->q = 4, opt->e = 2, opt->q2 = 24, opt->e2 = 1;
+	opt->zdrop = 400;
+	opt->end_bonus = -1;
+	opt->min_dp_max = opt->min_chain_score * opt->a;
+	opt->min_ksw_len = 200;
+	opt->anchor_ext_len = 20, opt->anchor_ext_shift = 6;
+	opt->mini_batch_size = 500000000;
+
+	opt->pe_ori = 0; // FF
+	opt->pe_bonus = 33;
+}
+
+void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)
+{
+	if ((opt->flag & MM_F_SPLICE_FOR) && (opt->flag & MM_F_SPLICE_REV))
+		opt->flag |= MM_F_SPLICE;
+	if (opt->mid_occ <= 0)
+		opt->mid_occ = mm_idx_cal_max_occ(mi, opt->mid_occ_frac);
+	if (mm_verbose >= 3)
+		fprintf(stderr, "[M::%s::%.3f*%.2f] mid_occ = %d\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), opt->mid_occ);
+}
+
+void mm_mapopt_max_intron_len(mm_mapopt_t *opt, int max_intron_len)
+{
+	if ((opt->flag & MM_F_SPLICE) && max_intron_len > 0)
+		opt->max_gap_ref = opt->bw = max_intron_len;
+}
+
+int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
+{
+	if (preset == 0) {
+		mm_idxopt_init(io);
+		mm_mapopt_init(mo);
+	} else if (strcmp(preset, "ava-ont") == 0) {
+		io->flag = 0, io->k = 15, io->w = 5;
+		mo->flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN;
+		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
+	} else if (strcmp(preset, "ava-pb") == 0) {
+		io->flag |= MM_I_HPC, io->k = 19, io->w = 5;
+		mo->flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN;
+		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
+	} else if (strcmp(preset, "map10k") == 0 || strcmp(preset, "map-pb") == 0) {
+		io->flag |= MM_I_HPC, io->k = 19;
+	} else if (strcmp(preset, "map-ont") == 0) {
+		io->flag = 0, io->k = 15;
+	} else if (strcmp(preset, "asm5") == 0) {
+		io->flag = 0, io->k = 19, io->w = 19;
+		mo->a = 1, mo->b = 19, mo->q = 39, mo->q2 = 81, mo->e = 3, mo->e2 = 1, mo->zdrop = 200;
+		mo->min_dp_max = 200;
+		mo->best_n = 50;
+	} else if (strcmp(preset, "asm10") == 0) {
+		io->flag = 0, io->k = 19, io->w = 19;
+		mo->a = 1, mo->b = 9, mo->q = 16, mo->q2 = 41, mo->e = 2, mo->e2 = 1, mo->zdrop = 200;
+		mo->min_dp_max = 200;
+		mo->best_n = 50;
+	} else if (strcmp(preset, "short") == 0 || strcmp(preset, "sr") == 0) {
+		io->flag = 0, io->k = 21, io->w = 11;
+		mo->flag |= MM_F_SR | MM_F_FRAG_MODE | MM_F_NO_PRINT_2ND | MM_F_2_IO_THREADS | MM_F_HEAP_SORT;
+		mo->pe_ori = 0<<1|1; // FR
+		mo->a = 2, mo->b = 8, mo->q = 12, mo->e = 2, mo->q2 = 24, mo->e2 = 1;
+		mo->zdrop = 100;
+		mo->end_bonus = 10;
+		mo->max_frag_len = 800;
+		mo->max_gap = 100;
+		mo->bw = 100;
+		mo->pri_ratio = 0.5f;
+		mo->min_cnt = 2;
+		mo->min_chain_score = 25;
+		mo->min_dp_max = 40;
+		mo->best_n = 20;
+		mo->mid_occ = 1000;
+		mo->max_occ = 5000;
+		mo->mini_batch_size = 50000000;
+	} else if (strcmp(preset, "splice") == 0 || strcmp(preset, "cdna") == 0) {
+		io->flag = 0, io->k = 15, io->w = 5;
+		mo->flag |= MM_F_SPLICE | MM_F_SPLICE_FOR | MM_F_SPLICE_REV | MM_F_SPLICE_FLANK;
+		mo->max_gap = 2000, mo->max_gap_ref = mo->bw = 200000;
+		mo->a = 1, mo->b = 2, mo->q = 2, mo->e = 1, mo->q2 = 32, mo->e2 = 0;
+		mo->noncan = 9;
+		mo->zdrop = 200;
+	} else return -1;
+	return 0;
+}
+
+int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
+{
+	if (mo->best_n < 0) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m -N must be no less than 0\033[0m\n");
+		return -4;
+	}
+	if (mo->best_n == 0 && mm_verbose >= 2)
+		fprintf(stderr, "[WARNING]\033[1;31m '-N 0' reduces mapping accuracy. Please use '--secondary=no' instead.\033[0m\n");
+	if (mo->pri_ratio < 0.0f || mo->pri_ratio > 1.0f) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m -p must be within 0 and 1 (including 0 and 1)\033[0m\n");
+		return -4;
+	}
+	if ((mo->flag & MM_F_FOR_ONLY) && (mo->flag & MM_F_REV_ONLY)) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m --for-only and --rev-only can't be applied at the same time\033[0m\n");
+		return -3;
+	}
+	if ((mo->q != mo->q2 || mo->e != mo->e2) && !(mo->e > mo->e2 && mo->q + mo->e < mo->q2 + mo->e2)) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m dual gap penalties violating E1>E2 and O1+E1<O2+E2\033[0m\n");
+		return -2;
+	}
+	if ((mo->q + mo->e) + (mo->q2 + mo->e2) > 127) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m scoring system violating ({-O}+{-E})+({-O2}+{-E2}) <= 127\033[0m\n");
+		return -1;
+	}
+	return 0;
+}
@@ -81,10 +81,13 @@ This constructor accepts the following arguments:

 .. code:: python

-	mappy.Aligner.map(seq)
+	mappy.Aligner.map(seq, seq2=None)

 This method aligns :code:`seq` against the index. It is a generator, *yielding*
-a series of :code:`mappy.Alignment` objects.
+a series of :code:`mappy.Alignment` objects. If :code:`seq2` is present, mappy
+performs paired-end alignment, assuming the two ends are in the FR orientation.
+Alignments of the two ends can be distinguished by the :code:`read_num` field
+(see below).

 Class mappy.Alignment
 ~~~~~~~~~~~~~~~~~~~~~
@@ -118,6 +121,9 @@ properties:
 * **is_primary**: if the alignment is primary (typically the best and the first
  to generate)

+* **read_num**: read number that the alignment corresponds to; 1 for the first
+  read and 2 for the second read
+
 * **cigar_str**: CIGAR string

 * **cigar**: CIGAR returned as an array of shape :code:`(n_cigar,2)`. The two
@@ -15,6 +15,7 @@ typedef struct {
 	int32_t blen, mlen, NM, ctg_len;
 	uint8_t mapq, is_primary;
 	int8_t strand, trans_strand;
+	int32_t seg_id;
 	int32_t n_cigar32;
 	uint32_t *cigar32;
 } mm_hitpy_t;
@@ -32,6 +33,7 @@ static inline void mm_reg2hitpy(const mm_idx_t *mi, mm_reg1_t *r, mm_hitpy_t *h)
 	h->NM = r->blen - r->mlen + r->p->n_ambi;
 	h->trans_strand = r->p->trans_strand == 1? 1 : r->p->trans_strand == 2? -1 : 0;
 	h->is_primary = (r->id == r->parent);
+	h->seg_id = r->seg_id;
 	h->n_cigar32 = r->p->n_cigar;
 	h->cigar32 = r->p->cigar;
 }
@@ -68,4 +70,35 @@ static inline void mm_reset_timer(void)
 	mm_realtime0 = realtime();
 }

+extern unsigned char seq_comp_table[256];
+static inline mm_reg1_t *mm_map_aux(const mm_idx_t *mi, const char *seq1, const char *seq2, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt)
+{
+	if (seq2 == 0) {
+		return mm_map(mi, strlen(seq1), seq1, n_regs, b, opt, NULL);
+	} else {
+		int _n_regs[2];
+		mm_reg1_t *regs[2];
+		char *seq[2];
+		int i, len[2];
+		len[0] = strlen(seq1);
+		len[1] = strlen(seq2);
+		seq[0] = (char*)seq1;
+		seq[1] = strdup(seq2);
+		for (i = 0; i < len[1]>>1; ++i) {
+			int t = seq[1][len[1] - i - 1];
+			seq[1][len[1] - i - 1] = seq_comp_table[(uint8_t)seq[1][i]];
+			seq[1][i] = seq_comp_table[t];
+		}
+		if (len[1]&1) seq[1][len[1]>>1] = seq_comp_table[(uint8_t)seq[1][len[1]>>1]];
+		mm_map_frag(mi, 2, len, (const char**)seq, _n_regs, regs, b, opt, NULL);
+		for (i = 0; i < _n_regs[1]; ++i)
+			regs[1][i].rev = !regs[1][i].rev;
+		*n_regs = _n_regs[0] + _n_regs[1];
+		regs[0] = (mm_reg1_t*)realloc(regs[0], sizeof(mm_reg1_t) * (*n_regs));
+		memcpy(&regs[0][_n_regs[0]], regs[1], _n_regs[1] * sizeof(mm_reg1_t));
+		free(regs[1]);
+		return regs[0];
+	}
+}
+
 #endif
@@ -5,7 +5,7 @@ cdef extern from "minimap.h":
 	# Options
 	#
 	ctypedef struct mm_idxopt_t:
-		short k, w, is_hpc, bucket_bits
+		short k, w, flag, bucket_bits
 		int mini_batch_size
 		uint64_t batch_size

@@ -30,6 +30,7 @@ cdef extern from "minimap.h":
 		int end_bonus
 		int min_dp_max
 		int min_ksw_len
+		int anchor_ext_len, anchor_ext_shift
 		int pe_ori, pe_bonus
 		float mid_occ_frac
 		int32_t mid_occ
@@ -51,7 +52,7 @@ cdef extern from "minimap.h":
 		pass

 	ctypedef struct mm_idx_t:
-		int32_t b, w, k, is_hpc
+		int32_t b, w, k, flag
 		uint32_t n_seq
 		mm_idx_seq_t *seq
 		uint32_t *S
@@ -78,7 +79,6 @@ cdef extern from "minimap.h":

 	mm_tbuf_t *mm_tbuf_init()
 	void mm_tbuf_destroy(mm_tbuf_t *b)
-	mm_reg1_t *mm_map(const mm_idx_t *mi, int l_seq, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *name)

 #
 # Helper header (because it is hard to expose mm_reg1_t with Cython)
@@ -91,11 +91,13 @@ cdef extern from "cmappy.h":
 		int32_t blen, mlen, NM, ctg_len
 		uint8_t mapq, is_primary
 		int8_t strand, trans_strand
+		int32_t seg_id
 		int32_t n_cigar32
 		uint32_t *cigar32

 	void mm_reg2hitpy(const mm_idx_t *mi, mm_reg1_t *r, mm_hitpy_t *h)
 	void mm_free_reg1(mm_reg1_t *r)
+	mm_reg1_t *mm_map_aux(const mm_idx_t *mi, const char *seq1, const char *seq2, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt)

 	ctypedef struct kstring_t:
 		unsigned l, m
@@ -10,9 +10,10 @@ cdef class Alignment:
 	cdef int _NM, _mlen, _blen
 	cdef int8_t _strand, _trans_strand
 	cdef uint8_t _mapq, _is_primary
+	cdef int _seg_id
 	cdef _ctg, _cigar # these are python objects

-	def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand):
+	def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand, seg_id):
 		self._ctg = ctg if isinstance(ctg, str) else ctg.decode()
 		self._ctg_len, self._r_st, self._r_en = cl, cs, ce
 		self._strand, self._q_st, self._q_en = strand, qs, qe
@@ -21,6 +22,7 @@ cdef class Alignment:
 		self._cigar = cigar
 		self._is_primary = is_primary
 		self._trans_strand = trans_strand
+		self._seg_id = seg_id

 	@property
 	def ctg(self): return self._ctg
@@ -64,6 +66,9 @@ cdef class Alignment:
 	@property
 	def cigar(self): return self._cigar

+	@property
+	def read_num(self): return self._seg_id + 1
+
 	@property
 	def cigar_str(self):
 		return "".join(map(lambda x: str(x[0]) + 'MIDNSH'[x[1]], self._cigar))
@@ -125,7 +130,7 @@ cdef class Aligner:
 	def __bool__(self):
 		return (self._idx != NULL)

-	def map(self, seq, buf=None):
+	def map(self, seq, seq2=None, buf=None):
 		cdef cmappy.mm_reg1_t *regs
 		cdef cmappy.mm_hitpy_t h
 		cdef ThreadBuffer b
@@ -134,7 +139,8 @@ cdef class Aligner:
 		if self._idx is NULL: return None
 		if buf is None: b = ThreadBuffer()
 		else: b = buf
-		regs = cmappy.mm_map(self._idx, len(seq), str.encode(seq), &n_regs, b._b, &self.map_opt, NULL)
+		if seq2 is None: regs = cmappy.mm_map_aux(self._idx, str.encode(seq), NULL, &n_regs, b._b, &self.map_opt)
+		else: regs = cmappy.mm_map_aux(self._idx, str.encode(seq), str.encode(seq2), &n_regs, b._b, &self.map_opt)

 		for i in range(n_regs):
 			cmappy.mm_reg2hitpy(self._idx, &regs[i], &h)
@@ -142,7 +148,7 @@ cdef class Aligner:
 			for k in range(h.n_cigar32):
 				c = h.cigar32[k]
 				cigar.append([c>>4, c&0xf])
-			yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.mlen, h.blen, h.NM, h.trans_strand)
+			yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.mlen, h.blen, h.NM, h.trans_strand, h.seg_id)
 			cmappy.mm_free_reg1(&regs[i])
 		free(regs)

@@ -56,6 +56,7 @@ sdust_buf_t *sdust_buf_init(void *km)
 	buf = (sdust_buf_t*)kcalloc(km, 1, sizeof(sdust_buf_t));
 	buf->km = km;
 	buf->w = kdq_init(int, buf->km);
+	kdq_resize(int, buf->w, 8);
 	return buf;
 }

@@ -23,7 +23,7 @@ def readme():

 setup(
 	name = 'mappy',
-	version = '2.5',
+	version = '2.8',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),
@@ -33,13 +33,13 @@ setup(
 	keywords = 'sequence-alignment',
 	scripts = ['python/minimap2.py'],
    ext_modules = [Extension('mappy',
-		sources = [module_src, 'align.c', 'bseq.c', 'chain.c', 'format.c', 'hit.c', 'index.c', 'pe.c',
+		sources = [module_src, 'align.c', 'bseq.c', 'chain.c', 'format.c', 'hit.c', 'index.c', 'pe.c', 'options.c',
 				   'ksw2_extd2_sse.c', 'ksw2_exts2_sse.c', 'ksw2_extz2_sse.c', 'ksw2_ll_sse.c',
-				   'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c'],
+				   'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c', 'esterr.c'],
 		depends = ['minimap.h', 'bseq.h', 'kalloc.h', 'kdq.h', 'khash.h', 'kseq.h', 'ksort.h',
 				   'ksw2.h', 'kthread.h', 'kvec.h', 'mmpriv.h', 'sdust.h',
 				   'python/cmappy.h', 'python/cmappy.pxd'],
-		extra_compile_args = ['-msse4'], # WARNING: ancient x86_64 CPUs don't have SSE4
+		extra_compile_args = ['-DHAVE_KALLOC', '-msse4'], # WARNING: ancient x86_64 CPUs don't have SSE4
 		include_dirs = ['.'],
 		libraries = ['z', 'm', 'pthread'])],
 	classifiers = [
@@ -2,8 +2,9 @@
 #include <stdlib.h>
 #include <assert.h>
 #include <string.h>
+#define __STDC_LIMIT_MACROS
 #include "kvec.h"
-#include "minimap.h"
+#include "mmpriv.h"

 unsigned char seq_nt4_table[256] = {
 	0, 1, 2, 3,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
@@ -113,22 +114,22 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i
 			}
 		} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
 		buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
-		if (l == w + k - 1) { // special case for the first window - because identical k-mers are not stored yet
+		if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
 			for (j = buf_pos + 1; j < w; ++j)
 				if (min.x == buf[j].x && buf[j].y != min.y) kv_push(mm128_t, km, *p, buf[j]);
 			for (j = 0; j < buf_pos; ++j)
 				if (min.x == buf[j].x && buf[j].y != min.y) kv_push(mm128_t, km, *p, buf[j]);
 		}
 		if (info.x <= min.x) { // a new minimum; then write the old min
-			if (l >= w + k) kv_push(mm128_t, km, *p, min);
+			if (l >= w + k && min.x != UINT64_MAX) kv_push(mm128_t, km, *p, min);
 			min = info, min_pos = buf_pos;
 		} else if (buf_pos == min_pos) { // old min has moved outside the window
-			if (l >= w + k - 1) kv_push(mm128_t, km, *p, min);
+			if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(mm128_t, km, *p, min);
 			for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers
 				if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer
 			for (j = 0; j <= buf_pos; ++j)
 				if (min.x >= buf[j].x) min = buf[j], min_pos = j;
-			if (l >= w + k - 1) { // write identical k-mers
+			if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
 				for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
 					if (min.x == buf[j].x && min.y != buf[j].y) kv_push(mm128_t, km, *p, buf[j]);
 				for (j = 0; j <= buf_pos; ++j)
@@ -268,18 +268,10 @@
 	Title = {Mason -- a read simulator for second generation sequencing data},
 	Year = {2010}}

-@article{Langmead:2012fk,
-	Author = {Langmead, Ben and Salzberg, Steven L},
-	Journal = {Nat Methods},
-	Pages = {357-9},
-	Title = {Fast gapped-read alignment with Bowtie 2},
-	Volume = {9},
-	Year = {2012}}
-
@article{Zaharia:2011aa,
 	Author = {Zaharia, Matei and others},
 	Journal = {arXiv:1111:5572},
-	Title = {Faster and More Accurate Sequence Alignment with SNAP},
+	Title = {Faster and More Accurate Sequence Alignment with {SNAP}},
 	Year = {2011}}

@article{Irimia:2008aa,
@@ -305,3 +297,19 @@
 	Title = {Versatile and open software for comparing large genomes},
 	Volume = {5},
 	Year = {2004}}
+
+@article {Li223297,
+	author = {Li, Heng and others},
+	title = {New synthetic-diploid benchmark for accurate variant calling evaluation},
+	year = {2017},
+	note = {doi:10.1101/223297},
+	journal = {bioRxiv}
+}
+
+@article{Berlin:2015xy,
+	Author = {Berlin, Konstantin and others},
+	Journal = {Nat Biotechnol},
+	Pages = {623-30},
+	Title = {Assembling large genomes with single-molecule sequencing and locality-sensitive hashing},
+	Volume = {33},
+	Year = {2015}}
@@ -13,7 +13,6 @@

 \usepackage{natbib}
 \bibliographystyle{apalike}
-\usepackage{hyperref}

 \DeclareMathOperator*{\argmax}{argmax}

@@ -138,7 +137,7 @@ $h=50$; even if the heuristic fails, the optimal chain is often close.

 \subsubsection{Backtracking}
 Let $P(i)$ be the index of the best predecessor of anchor $i$. It equals 0 if
-$f(i)=w_i$ or $\argmax_j\{f(j)+\eta(j,i)-\gamma(j,i)\}$ otherwise. For each
+$f(i)=w_i$ or $\argmax_j\{f(j)+\alpha(j,i)-\beta(j,i)\}$ otherwise. For each
 anchor $i$ in the descending order of $f(i)$, we apply $P(\cdot)$ repeatedly to
 find its predecessor and mark each visited $i$ as `used', until $P(i)=0$ or we
 reach an already `used' $i$. This way we find all chains with no anchors used
@@ -157,6 +156,53 @@ add the chain to $Q$. In the end, $Q$ contains all the primary chains. We did
 not choose a more sophisticated data structure (e.g. range tree or k-d tree)
 because this step is not the performance bottleneck.

+\subsubsection{Estimating per-base sequence divergence}
+Suppose a query sequence harbors $n$ seeds of length $k$, $m$ of which are
+present in a chain. We want to estimate the sequence divergence $\epsilon$
+between the query and the reference sequences in the chain. This is useful
+when base-level alignment is too expensive to perform.
+
+If we model substitutions with a homogeneous Poisson process along the query
+sequence, the probablity of seeing $k$ consecutive bases without substitutions
+is $e^{-k\epsilon}$. On the assumption that all $k$-mers are independent of
+each other, the likelihood function of $\epsilon$ is
+\[
+\mathcal{L}(\epsilon|n,m,k)=e^{-m\cdot k\epsilon}(1-e^{-k\epsilon})^{n-m}
+\]
+The maximum likelihood estimate of $\epsilon$ is
+\[
+\hat{\epsilon}=\frac{1}{k}\log\frac{n}{m}
+\]
+In reality, sequencing errors are sometimes clustered and $k$-mers are not
+independent of each other, especially when we take minimizers as seeds. These
+violate the assumptions in the derivation above. As a result, $\hat{\epsilon}$
+is only approximate and can be biased. It also ignores long deletions from the
+reference sequence. In practice, fortunately, $\hat{\epsilon}$ is often close
+to and strongly correlated with the sequence divergence estimated from
+base-level alignments. On the several datasets used in
+Section~\ref{sec:long-genomic}, the Spearman correlation coefficient is around
+$0.9$.
+
+\subsubsection{Indexing with homopolymer compressed $k$-mers}
+SmartDenovo
+(\href{https://github.com/ruanjue/smartdenovo}{https://github.com/ruanjue/smartdenovo};
+J Ruan, personal communication) indexes reads with homopolymer-compressed (HPC)
+$k$-mers and finds the strategy improves overlap sensitivity for SMRT reads.
+Minimap2 adopts the same heuristic.
+
+The HPC string of a string $s$, denoted by ${\rm HPC}(s)$, is constructed by
+contracting homopolymers in $s$ to a single base.  An HPC $k$-mer of $s$ is a
+$k$-long substring of ${\rm HPC}(s)$. For example, suppose $s={\tt GGATTTTCCA}$,
+${\rm HPC}(s)={\tt GATCA}$ and the first HPC 4-mer is ${\tt GATC}$.
+
+To demonstrate the effectiveness of HPC $k$-mers, we performed read overlapping
+for the example {\it E. coli} SMRT reads from PBcR~\citep{Berlin:2015xy}, using
+different types of $k$-mers. With normal 15bp minimizers per 5bp window,
+minimap2 finds 90.9\% of $\ge$2kb overlaps inferred from the read-to-reference
+alignment. With HPC 19-mers, minimap2 finds 97.4\% of overlaps. It achieves this
+higher sensitivity by indexing 1/3 fewer minimizers, which further helps
+performance. HPC-based indexing reduces the sensitivity for ONT reads, though.
+
 \subsection{Aligning genomic DNA}\label{sec:genomic}

 \subsubsection{Alignment with 2-piece affine gap cost}
@@ -397,7 +443,7 @@ consistent paired-end alignments.

 \section{Results}

-\subsection{Aligning long genomic reads}
+\subsection{Aligning long genomic reads}\label{sec:long-genomic}

 \begin{figure}[!tb]
 \centering
@@ -500,14 +546,18 @@ more junctions with a higher percentage being exactly or approximately correct.
 Minimap2 is over 40 times faster than GMAP and SpAln. While STAR is close to
 minimap2 in speed, it does not work well with noisy reads.

-We have also evaluated spliced aligners on public Iso-Seq data (human Alzheimer
-brain from \href{http://bit.ly/isoseqpub}{http://bit.ly/isoseqpub}). The
-observation is similar: minimap2 is faster at higher junction accuracy.
-On a private Nanopore Direct RNA data set with $\sim$17\% sequencing error rate
-(N. Loman, personal communication), minimap2 aligned 96\,467 introns
-from 37\,068 mapped reads with 95.4\% of them consistent with human gene
-annotations. In comparison, only 74.8\% of GMAP introns found in known gene
-annotations.
+We have also evaluated spliced aligners on a human Nanopore Direct RNA-seq
+dataset (\href{http://bit.ly/na12878ont}{http://bit.ly/na12878ont}). Minimap2
+aligned 10 million reads in $<$1 wall-clock hour using 16 CPU cores. 94.2\% of
+aligned splice junctions consistent with gene annotations. In comparison,
+GMAP under option `-k 14 -n 0 --min-intronlength 30 --cross-species' is 160
+times slower; 68.7\% of GMAP junctions are found in known gene annotations. The
+percentage increases to 84.1\% if an aligned junction within 10bp from an
+annotated junction is considered to be correct. On a public Iso-Seq dataset
+(human Alzheimer brain from
+\href{http://bit.ly/isoseqpub}{http://bit.ly/isoseqpub}), minimap2 is also
+faster at higher junction accuracy in comparison to other aligners in
+Table~\ref{tab:intron}.

 We noted that GMAP and SpAln have not been optimized for noisy reads. We are
 showing the best setting we have experimented, but their developers should be
@@ -551,8 +601,7 @@ with GATK HaplotypeCaller v3.5~\citep{Depristo:2011vn}. This run was sequenced
 from experimentally mixed CHM1 and CHM13 cell lines. Both of them are homozygous
 across the whole genome and have been \emph{de novo} assembled with SMRT reads
 to high quality. This allowed us to construct an independent truth variant
-data set
-(\href{https://github.com/lh3/CHM-eval}{https://github.com/lh3/CHM-eval}) for
+dataset~\citep{Li223297} for
 ERR1341796. In this evaluation, minimap2 has higher SNP false negative rate
 (FNR; 2.5\% of minimap2 vs 2.2\% of BWA-MEM), but fewer false positive SNPs per
 million bases (FPPM; 3.0 vs 3.9), lower 2--50bp INDEL FNR (7.3\% vs 7.5\%) and
@@ -597,7 +646,7 @@ uniqueness and reduce unsuccessful extensions. Minimap2 indexes reference
 k-mers with a hash table instead. Such fixed-length seeds are inferior to
 variable-length seeds in theory, but can be computed much more efficiently in
 practice. When a query sequence has multiple seed hits, we can afford to skip
-some highly repetitive seeds without affecting the final accuracy. This further
+highly repetitive seeds without affecting the final accuracy. This further
 alleviates the concern with the uniqueness of seeds. Hash table is the ideal
 data structure for mapping long query sequences.

@@ -605,8 +654,8 @@ data structure for mapping long query sequences.
 We owe a debt of gratitude to H. Suzuki and M. Kasahara for releasing their
 masterpiece and insightful notes before formal publication. We thank M.
 Schatz, P. Rescheneder and F. Sedlazeck for pointing out the limitation of
-BWA-MEM. We are also grateful to early minimap2 testers who have greatly helped
-to suggest features and to fix various issues.
+BWA-MEM. We are also grateful to minimap2 users who have greatly helped to
+suggest features and to fix various issues.

 \bibliography{minimap2}
Author	SHA1	Message	Date
Heng Li	53ce317e59	Release minimap2-2.8 (r672)	2018-02-01 12:50:20 -05:00
Heng Li	da6947cfa3	r671: cleanup command line options	2018-01-31 13:59:52 -05:00
Heng Li	46d6349af4	r670: added PE support to mappy and minor code cleanup	2018-01-31 11:33:08 -05:00
Heng Li	12a5a5fa3c	r669: improved self chain extension (#10 ) This has not fully resolved #10, only alleviated the issue.	2018-01-30 20:05:02 -05:00
Heng Li	ad18fa490d	fixed typos	2018-01-30 10:11:50 -05:00
Heng Li	43bfa6199d	r667: warn if one query file has fewer records #92	2018-01-28 17:36:21 -05:00
Heng Li	72b9b0e3b6	r666: report if >=3 query files in SR mode #92	2018-01-28 17:15:57 -05:00
Heng Li	6205fa6f21	document --heap-sort	2018-01-26 15:15:40 -05:00
Heng Li	d676a5314b	r664: use --heat-sort for sr by default	2018-01-26 12:25:42 -05:00
Heng Li	dfc78b39d3	refactor the old sorting	2018-01-26 09:37:48 -05:00
Heng Li	7b57c9a619	heap sort working on MT	2018-01-26 09:21:45 -05:00
Heng Li	123bc1d91d	put option operations in another file	2018-01-26 08:38:37 -05:00
Heng Li	dd18307e66	code backup	2018-01-25 21:52:49 -05:00
Heng Li	543fa12e68	r659: for C++ compatibility	2018-01-19 10:40:18 -05:00
Heng Li	af1a871270	r658: gives a warning if -N0 is used	2018-01-19 08:33:20 -05:00
Heng Li	2b71181a37	r657: check -p (#96 ) Well, in principle, every option should be checked. Will do when someone raise issues...	2018-01-19 01:03:38 -05:00
Heng Li	0454e6be91	explain -M in the manpage	2018-01-18 11:47:11 -05:00
Heng Li	33f8157961	r655: options to map to one strand of the ref #91	2018-01-16 10:34:30 -05:00
Heng Li	eecc06086f	Released minimap2-2.7 (r654)	2018-01-09 13:16:00 -05:00
Heng Li	dfea113f28	r653: the last change may write "N" wrongly	2018-01-08 11:33:53 -05:00
Heng Li	1842d7f5b5	allow to exclude regions	2018-01-07 22:35:56 -05:00
Heng Li	f5cfd439ee	r651: incorrectly treat introns as deletions This happened when the last operation during backtracking is an intron.	2018-01-07 19:42:50 -05:00
Heng Li	248b43cc47	work with targets in BED12	2018-01-07 19:39:25 -05:00
Heng Li	bf72969ab1	added a bed counter I know there are tools for this purpose, but they don't quite meet my exact need.	2018-01-07 15:19:18 -05:00
Heng Li	7b5a601d48	support paired-end reads this gives unnecessary warnings, which will be fixed later.	2018-01-07 13:53:26 -05:00
Heng Li	405d531100	don't test python v3.3	2018-01-05 22:17:00 -05:00
Heng Li	e9607fcd9b	allow to convert read names ONT read names are just too long and too hard to compress	2018-01-05 22:05:15 -05:00
Heng Li	a465a920ec	bug in block starts; added color and short-name	2018-01-05 21:45:56 -05:00
Heng Li	b20839be77	more robust ID conversion	2018-01-05 17:55:24 -05:00
Heng Li	209beb9955	convert EnsEMBL to UCSC name (optional)	2018-01-05 17:41:18 -05:00
Heng Li	cfe87f50c1	convert GTF/GFF3 to BED12	2018-01-05 17:10:43 -05:00
Heng Li	680b971bb0	deleted duplicated entries	2018-01-01 19:52:21 -05:00
Heng Li	6b0d3c1fa8	added PHONY and comments to makefile	2017-12-30 20:54:11 -05:00
Heng Li	dc9e3dcf4a	r639: changed -O/-E validation	2017-12-30 20:39:29 -05:00
Heng Li	cc75c12905	r638: disabled scoring checking I haven't figured out the exact bounds...	2017-12-30 07:50:40 -05:00
Heng Li	f159e1c2d3	new section on HPC k-mers	2017-12-24 19:11:23 -05:00
Heng Li	3a375d3436	renamed paf2ovlp to ov-eval	2017-12-24 18:04:00 -05:00
Heng Li	626f10e0d0	evaluate sensitivity in the same script	2017-12-24 17:55:42 -05:00
Heng Li	b997578078	find reads overlaps based on reference alignment	2017-12-24 17:20:27 -05:00
Heng Li	ce8a48d715	fixed two minor typos in references	2017-12-24 13:14:47 -05:00
Heng Li	99879e9e75	a new section on estimating sequence divergence	2017-12-24 12:57:38 -05:00
Heng Li	c969d1a1ce	updated direct RNA-seq results; cite syndip and a few minor changes	2017-12-24 11:06:17 -05:00
Heng Li	fcac296c4a	clarify ARM-NEON support in README	2017-12-18 23:01:44 -05:00
Heng Li	e420b17496	r629: API to construct index from strings	2017-12-18 22:29:46 -05:00
Heng Li	23a846c594	Merge pull request #81 from hasindu2008/master minimap2 on ARM processors	2017-12-16 09:40:41 -05:00
Hasindu Gamaarachchi	8995e2e078	added support for arm neon	2017-12-15 17:42:19 +11:00
Heng Li	ab345e600b	r626: function to check incorrect scoring system	2017-12-13 12:23:43 -05:00
Heng Li	d003a00d71	r625: HPC sketch still has one minor issue	2017-12-13 09:40:42 -05:00
Heng Li	ae85dcde76	added esterr.o dependencies	2017-12-13 09:01:48 -05:00
Heng Li	eb819c29e8	Release minimap2-2.6 (r623)	2017-12-12 11:09:59 -05:00
Heng Li	fb630de40a	r622: fixed bug in sdust due to recent refactor	2017-12-11 15:32:28 -05:00
Heng Li	43960a8ca7	r621: --print-qname also shows kalloc status	2017-12-11 12:30:08 -05:00
Heng Li	f6608fe99c	r620: revamped thread-local memory management * Don't preallocate sdust_buf or minizer list. kalloc should be fast enough - benchmarks needed to confirm. * Fixed a memory leak caused by divergence estimate (post v2.5) * Reset the kalloc buffer after mapping a long query. This reduces peak memory when large chunks of memory are allocated, at the cost of performance, though.	2017-12-11 12:11:10 -05:00
Heng Li	98a6e52c06	r618: heuristics to avoid tiny terminal exons	2017-12-11 00:57:55 -05:00
Heng Li	824712a4ee	r617: removed some unused code	2017-12-10 17:54:50 -05:00
Heng Li	0e42628ef6	r611: document --idx-no-seq; better inv aln	2017-12-08 13:16:18 -05:00
Heng Li	98a999fe44	r611: added pseudocount when est divergence	2017-12-08 12:57:57 -05:00
Heng Li	fec7bd713f	r610: warning if db sequence is 0-lengthed (#69 )	2017-12-07 21:05:39 -05:00
Heng Li	2f693e8ca4	r609: bugfix - SDUST masking not working	2017-12-07 11:45:38 -05:00
Heng Li	704ff9f4c6	r607: estimate sequence divergence Currently using the simplest method. There may be a more accurate estimate.	2017-12-06 16:14:39 -05:00
Heng Li	68c63f2d68	r606: fixed a sketch bug for long 256bp k-mer sketch() writes {-1,-1} to the output array.	2017-12-06 16:13:29 -05:00
Heng Li	76206f574f	support both SAM and PAF as input	2017-12-02 21:57:19 -05:00
Heng Li	e575f884e1	convert spliced PAF to BED	2017-12-02 21:23:41 -05:00
Heng Li	571161d5a2	changed badges	2017-12-01 16:16:27 -05:00
Heng Li	07d41efc2b	explain secondary/supplementary aln for RNA-seq	2017-11-30 23:02:20 -05:00
Heng Li	984f7846c0	r601: bugfix - a similar issue to r600 This bug unsets the alignment score of suboptimal alignments.	2017-11-30 11:51:34 -05:00
Heng Li	af1d6afba9	r600: bugfix - missing secondary alignments (#71 ) This should very rarely happen to typical data, but has a higher chance in artifactual data.	2017-11-30 11:34:10 -05:00
Heng Li	cbdb6c069f	when there are incorrect anno, warn but not abort	2017-11-24 12:48:49 -05:00
Heng Li	35b6d9f7d5	direct manpage to HTML	2017-11-24 11:05:20 -05:00
Heng Li	39a9666246	convert PAF to LAST's cigar output also added the support of BLAST-like output for "--cs=short".	2017-11-18 20:30:44 -05:00
Heng Li	662d05dc02	fixed incorrect var coordinate	2017-11-12 20:25:33 -05:00
Heng Li	379457c18b	use mapq threshold	2017-11-12 19:01:34 -05:00
Heng Li	03169d590b	print cov-1 regions (not BED output any more)	2017-11-12 18:45:49 -05:00
Heng Li	8c8d446820	find regions covered by one contig	2017-11-12 18:41:01 -05:00
Heng Li	0ddb064f17	dnadiff-like script; improvement coming	2017-11-12 15:07:29 -05:00
Heng Li	131cfc6938	r574: build index without sequences	2017-11-11 21:38:38 -05:00
Heng Li	2f463b1db0	r573: prepare to generalize index	2017-11-11 19:54:06 -05:00