Released minimap2-2.7 (r654)

r653: the last change may write "N" wrongly
allow to exclude regions
2026-09-26 04:18:12 +08:00 · 2018-01-09 13:16:00 -05:00 · 2018-01-08 11:33:53 -05:00 · 2018-01-07 22:35:56 -05:00 · 2018-01-07 19:42:50 -05:00 · 2018-01-07 19:39:25 -05:00
32 changed files with 3205 additions and 225 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 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,6 +63,19 @@ 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*
@@ -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
@@ -1,3 +1,60 @@
 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)
+[![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](https://img.shields.io/conda/vn/bioconda/minimap2.svg?style=flat)](https://anaconda.org/bioconda/minimap2)
+[![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.7/minimap2-2.7_x64-linux.tar.bz2 \
  | tar -jxvf -
-./minimap2-2.5_x64-linux/minimap2
+./minimap2-2.7_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,8 +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
+* Minimap2 requires SSE2 instructions on x86 CPUs or NEON on ARM CPUs. It is
-  non-SSE2 support, but it would make minimap2 slower by several times.
+  possible to add non-SIMD 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
@@ -363,3 +378,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);
@@ -470,7 +522,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 +585,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 +593,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 +614,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 +633,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 +653,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 +667,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 +677,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, &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[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 +691,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 (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 +705,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);
@@ -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);
@@ -31,19 +31,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 +89,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 +214,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 +246,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 +253,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;
+			if (!(p->mi->flag & MM_I_NO_SEQ)) {
-			old_max_len = (p->sum_len + 7) / 8;
+				uint64_t sum_len, old_max_len, max_len;
-			max_len = (p->sum_len + sum_len + 7) / 8;
+				for (i = 0, sum_len = 0; i < s->n_seq; ++i) sum_len += s->seq[i].l_seq;
-			kroundup64(old_max_len); kroundup64(max_len);
+				old_max_len = (p->sum_len + 7) / 8;
-			if (old_max_len != max_len) {
+				max_len = (p->sum_len + sum_len + 7) / 8;
-				p->mi->S = (uint32_t*)realloc(p->mi->S, max_len * 4);
+				kroundup64(old_max_len); kroundup64(max_len);
-				memset(&p->mi->S[old_max_len], 0, 4 * (max_len - old_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
+				if (!(p->mi->flag & MM_I_NO_SEQ)) {
-					uint64_t o = p->sum_len + j;
+					for (j = 0; j < seq->len; ++j) { // TODO: this is not the fastest way, but let's first see if speed matters here
-					int c = seq_nt4_table[(uint8_t)s->seq[i].seq[j]];
+						uint64_t o = p->sum_len + j;
-					mm_seq4_set(p->mi->S, o, c);
+						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 +292,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 +308,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 +329,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 +393,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 +420,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 +471,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);
+	if (!(mi->flag & MM_I_NO_SEQ)) {
-	fread(mi->S, 4, (sum_len + 7) / 8, fp);
+		mi->S = (uint32_t*)malloc((sum_len + 7) / 8 * 4);
 		fread(mi->S, 4, (sum_len + 7) / 8, fp);
 	}
 	return mi;
 }
@@ -476,10 +528,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)
+	for (p = km->core_head; p != NULL; p = p->ptr) {
-		++s->n_cores, s->capacity += p->size * sizeof(header_t);
+		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);
@@ -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 == 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 && with_N) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 3, 1), --i; // intron
+		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);
 	}
 }
@@ -6,7 +6,7 @@
 #include "mmpriv.h"
 #include "getopt.h"
-#define MM_VERSION "2.5-r572"
+#define MM_VERSION "2.7-r654"
 #ifdef __linux__
 #include <sys/resource.h>
@@ -43,6 +43,8 @@ static struct option long_options[] = {
 	{ "end-bonus",      required_argument, 0, 0 },
 	{ "no-pairing",     no_argument,       0, 0 },
 	{ "splice-flank",   optional_argument, 0, 0 },
 	{ "idx-no-seq",     no_argument,       0, 0 },
 	{ "end-seed-pen",   required_argument, 0, 0 },   // 21
 	{ "help",           no_argument,       0, 'h' },
 	{ "max-intron-len", required_argument, 0, 'G' },
 	{ "version",        no_argument,       0, 'V' },
@@ -94,7 +96,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);
@@ -132,12 +134,14 @@ int main(int argc, char *argv[])
 		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 == 14) { // --frag
 			if (optarg == 0 || strcmp(optarg, "yes") == 0 || strcmp(optarg, "y") == 0)
 				opt.flag |= MM_F_FRAG_MODE;
@@ -196,6 +200,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");
@@ -259,11 +267,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);
@@ -36,6 +36,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	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
@@ -64,29 +65,29 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		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;
+		io->flag = 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;
+		io->flag |= MM_I_HPC, 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;
+		io->flag |= MM_I_HPC, io->k = 19;
 	} else if (strcmp(preset, "map-ont") == 0) {
-		io->is_hpc = 0, io->k = 15;
+		io->flag = 0, io->k = 15;
 	} else if (strcmp(preset, "asm5") == 0) {
-		io->is_hpc = 0, io->k = 19, io->w = 19;
+		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->is_hpc = 0, io->k = 19, io->w = 19;
+		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->is_hpc = 0, io->k = 21, io->w = 11;
+		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;
 		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;
@@ -104,7 +105,7 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		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;
+		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;
@@ -114,6 +115,21 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 	return 0;
 }
 int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
 {
 	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;
 }
 typedef struct {
 	uint32_t n;
 	uint32_t qpos;
@@ -122,8 +138,6 @@ typedef struct {
 } mm_match_t;
 struct mm_tbuf_s {
 	sdust_buf_t *sdb;
 	mm128_v mini;
 	void *km;
 };
@@ -132,24 +146,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,45 +176,49 @@ 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);
+		mm_sketch(km, seqs[i], qlens[i], mi->w, mi->k, i, mi->flag&MM_I_HPC, mv);
-		for (j = n; j < b->mini.n; ++j)
+		for (j = n; j < mv->n; ++j)
-			b->mini.a[j].y += sum << 1;
+			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);
+			mv->n = n + mm_dust_minier(km, mv->n - n, mv->a + n, qlens[i], seqs[i], opt->sdust_thres);
-		sum += qlens[i], n = b->mini.n;
+		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)
+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 rep_st = 0, rep_en = 0, i;
 	mm_match_t *m;
 	mm128_t *a;
-	m = (mm_match_t*)kmalloc(b->km, b->mini.n * sizeof(mm_match_t));
+	*n_mini_pos = 0;
-	for (i = 0; i < b->mini.n; ++i) {
+	*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 = 0; i < mv->n; ++i) {
 		int t;
-		mm128_t *p = &b->mini.a[i];
+		mm128_t *p = &mv->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[]
+	for (i = 0, *n_a = 0; i < mv->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));
+	a = (mm128_t*)kmalloc(km, *n_a * sizeof(mm128_t));
-	for (i = *rep_len = 0, *n_a = 0; i < b->mini.n; ++i) {
+	for (i = *rep_len = 0, *n_a = 0; i < mv->n; ++i) {
-		mm128_t *p = &b->mini.a[i];
+		mm128_t *p = &mv->a[i];
 		mm_match_t *q = &m[i];
 		const uint64_t *r = q->cr;
 		int k, q_span = p->x & 0xff, is_tandem = 0;
@@ -213,8 +230,9 @@ static mm128_t *collect_seed_hits(const mm_mapopt_t *opt, int max_occ, const mm_
 			} else rep_en = en;
 			continue;
 		}
-		if (i > 0 && p->x>>8 == b->mini.a[i - 1].x>>8) is_tandem = 1;
+		(*mini_pos)[(*n_mini_pos)++] = (uint64_t)q_span<<32 | q->qpos>>1;
-		if (i < b->mini.n - 1 && p->x>>8 == b->mini.a[i + 1].x>>8) is_tandem = 1;
+		if (i > 0 && p->x>>8 == mv->a[i - 1].x>>8) is_tandem = 1;
 		if (i < mv->n - 1 && p->x>>8 == mv->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;
@@ -240,7 +258,7 @@ static mm128_t *collect_seed_hits(const mm_mapopt_t *opt, int max_occ, const mm_
 		}
 	}
 	*rep_len += rep_en - rep_st;
-	kfree(b->km, m);
+	kfree(km, m);
 	return a;
 }
@@ -258,7 +276,7 @@ static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_i
 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)
 {
 	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()
+	regs = mm_align_skeleton(km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs()
 	if (!(opt->flag & MM_F_AVA)) {
 		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);
@@ -269,13 +287,15 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k
 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)
 {
-	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,8 +306,8 @@ 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);
+	collect_minimizers(b->km, opt, mi, n_segs, qlens, seqs, &mv);
-	a = collect_seed_hits(opt, opt->mid_occ, mi, qname, qlen_sum, &n_a, &rep_len, b);
+	a = collect_seed_hits(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 	radix_sort_128x(a, a + n_a);
 	if (mm_dbg_flag & MM_DBG_PRINT_SEED) {
@@ -327,7 +347,8 @@ 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);
+			kfree(b->km, mini_pos);
 			a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 			radix_sort_128x(a, a + n_a);
 			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,6 +363,7 @@ 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);
@@ -361,8 +383,21 @@ 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)
@@ -396,14 +431,13 @@ 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);
+	int qlens[MM_MAX_SEG], 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;
+	const char *qseqs[MM_MAX_SEG], *quals[MM_MAX_SEG];
 	mm_tbuf_t *b = s->buf[tid];
 	assert(s->n_seg[i] <= MM_MAX_SEG);
 	memset(quals, 0, sizeof(char*) * MM_MAX_SEG);
 	if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
-		fprintf(stderr, "QR\t%s\t%d\n", s->seq[off].name, tid);
+		fprintf(stderr, "QR\t%s\t%d\t%d\n", s->seq[off].name, tid, s->seq[off].l_seq);
 	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**));
 	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]);
@@ -429,9 +463,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)
@@ -26,6 +26,10 @@
 #define MM_F_SPLICE_FLANK  0x40000
 #define MM_F_SOFTCLIP      0x80000
 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2
 #define MM_I_NO_NAME      0x4
 #define MM_IDX_MAGIC   "MMI\2"
 #define MM_MAX_SEG       255
@@ -46,7 +50,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 +68,25 @@ typedef struct {
 } mm_extra_t;
 typedef struct {
-	int32_t id;                     // ID for internal uses (see also parent below)
+	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
+	int32_t cnt;            // 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 rid;            // reference index; if this is an alignment from inversion rescue
-	int32_t score;                  // DP alignment score
+	int32_t score;          // DP alignment score
-	int32_t qs, qe, rs, re;         // query start and end; reference start and end
+	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 parent, subsc;  // parent==id if primary; best alternate mapping score
-	int32_t as;                     // offset in the a[] array (for internal uses only)
+	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 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
+	int32_t n_sub;          // number of suboptimal mappings
-	uint32_t pe_thru:1, score0:31;
+	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, dummy:16;
 	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 +116,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 +155,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 +208,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
 *
@@ -269,7 +292,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 "9 January 2018" "minimap2-2.7 (r654)" "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
@@ -264,6 +272,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
@@ -461,7 +481,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
@@ -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;
+var c, line_len = 80, fmt = "aln";
-while ((c = getopt(arguments, "ml:")) != null) {
+while ((c = getopt(arguments, "f:l:")) != null) {
-	if (c == 'm') maf_out = true;
+	if (c == 'f') {
-	else if (c == 'l') line_len = parseInt(getopt.arg); // TODO: not implemented yet
+		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 s_ref = new Bytes(), s_qry = new Bytes(), s_mid = new Bytes(); // these are used to show padded alignment
-var re = /([=\-\+\*])([A-Za-z]+)/g;
+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]);
+var file = arguments[getopt.ind] == "-"? new File() : new File(arguments[getopt.ind]);
-if (maf_out) print("##maf version=1\n");
+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;
+	++lineno;
 	var cs = m[1];
 	s_ref.length = s_qry.length = s_mid.length = 0;
 	var slen = [0, 0], elen = [0, 0];
-	if (maf_out) {
+	if (fmt == "lastz-cigar") { // LASTZ-cigar output
-		while ((m = re.exec(cs)) != null)
+		var cg = (m = /\tcg:Z:(\S+)/.exec(line)) != null? m[1] : null;
-			update_aln(s_ref, s_qry, s_mid, m[1], m[2], elen);
+		if (cg == null) {
-		if (maf_out) {
+			warn("WARNING: converting to LASTZ-cigar format requires the 'cg' tag, which is absent on line " + lineno);
-			var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : 0;
+			continue;
 			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 {
+		var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? m[1] : 0;
-		line = line.replace(/\tc[sg]:Z:\S+/g, "");
+		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(optarg.arg);
 	else if (c == 'g') gap_thres = parseInt(optarg.arg);
 	else if (c == 'q') min_mapq = parseInt(optarg.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 = 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);
@@ -64,7 +64,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 +79,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);
@@ -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
@@ -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.7',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),
@@ -35,7 +35,7 @@ setup(
    ext_modules = [Extension('mappy',
 		sources = [module_src, 'align.c', 'bseq.c', 'chain.c', 'format.c', 'hit.c', 'index.c', 'pe.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'],
@@ -113,22 +113,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
+We have also evaluated spliced aligners on a human Nanopore Direct RNA-seq
-brain from \href{http://bit.ly/isoseqpub}{http://bit.ly/isoseqpub}). The
+dataset (\href{http://bit.ly/na12878ont}{http://bit.ly/na12878ont}). Minimap2
-observation is similar: minimap2 is faster at higher junction accuracy.
+aligned 10 million reads in $<$1 wall-clock hour using 16 CPU cores. 94.2\% of
-On a private Nanopore Direct RNA data set with $\sim$17\% sequencing error rate
+aligned splice junctions consistent with gene annotations. In comparison,
-(N. Loman, personal communication), minimap2 aligned 96\,467 introns
+GMAP under option `-k 14 -n 0 --min-intronlength 30 --cross-species' is 160
-from 37\,068 mapped reads with 95.4\% of them consistent with human gene
+times slower; 68.7\% of GMAP junctions are found in known gene annotations. The
-annotations. In comparison, only 74.8\% of GMAP introns found in known gene
+percentage increases to 84.1\% if an aligned junction within 10bp from an
-annotations.
+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
+dataset~\citep{Li223297} for
 (\href{https://github.com/lh3/CHM-eval}{https://github.com/lh3/CHM-eval}) 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
+BWA-MEM. We are also grateful to minimap2 users who have greatly helped to
-to suggest features and to fix various issues.
+suggest features and to fix various issues.
 \bibliography{minimap2}
Author	SHA1	Message	Date
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