Release minimap2-2.9 (r720)

r719: fixed bad memory access
r718: retrieve sequence from the index
2026-09-26 03:58:11 +08:00 · 2018-02-24 09:31:09 -05:00 · 2018-02-23 17:27:41 -05:00 · 2018-02-23 10:18:26 -05:00 · 2018-02-20 11:05:54 -05:00 · 2018-02-20 09:41:25 -05:00
40 changed files with 3212 additions and 2430 deletions
@@ -1,7 +1,7 @@
 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 esterr.o ksw2_ll_sse.o
+OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o chain.o align.o hit.o map.o format.o pe.o esterr.o ksw2_ll_sse.o
 PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite
 LIBS=		-lm -lz -lpthread
@@ -77,7 +77,7 @@ ksw2_exts2_neon.o:ksw2_exts2_sse.c ksw2.h kalloc.h
 # other non-file targets

 clean:
-		rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist mappy.so mappy.c python/mappy.c mappy.egg*
+		rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist mappy*.so mappy.c python/mappy.c mappy.egg*

 depend:
 		(LC_ALL=C; export LC_ALL; makedepend -Y -- $(CFLAGS) $(CPPFLAGS) -- *.c)
@@ -98,9 +98,12 @@ ksw2_extd2_sse.o: ksw2.h kalloc.h
 ksw2_exts2_sse.o: ksw2.h kalloc.h
 ksw2_extz2_sse.o: ksw2.h kalloc.h
 ksw2_ll_sse.o: ksw2.h kalloc.h
+kthread.o: kthread.h
 main.o: bseq.h minimap.h mmpriv.h getopt.h
 map.o: kthread.h kvec.h kalloc.h sdust.h mmpriv.h minimap.h bseq.h khash.h
-misc.o: minimap.h ksort.h
+map.o: ksort.h
+misc.o: mmpriv.h minimap.h bseq.h ksort.h
+options.o: mmpriv.h minimap.h bseq.h
 pe.o: mmpriv.h minimap.h bseq.h kvec.h kalloc.h ksort.h
 sdust.o: kalloc.h kdq.h kvec.h sdust.h
-sketch.o: kvec.h kalloc.h minimap.h
+sketch.o: kvec.h kalloc.h mmpriv.h minimap.h bseq.h
@@ -1,3 +1,55 @@
+Release 2.9-r720 (23 February 2018)
+-----------------------------------
+
+This release fixed multiple minor bugs.
+
+* Fixed two bugs that lead to incorrect inversion alignment. Also improved the
+  sensitivity to small inversions by using double Z-drop cutoff (#112).
+
+* Fixed an issue that may cause the end of a query sequence unmapped (#104).
+
+* Added a mappy API to retrieve sequences from the index (#126) and to reverse
+  complement DNA sequences. Fixed a bug where the `best_n` parameter did not
+  work (#117).
+
+* Avoided segmentation fault given incorrect FASTQ input (#111).
+
+* Combined all auxiliary javascripts to paftools.js. Fixed several bugs in
+  these scripts at the same time.
+
+(2.9: 24 February 2018, r720)
+
+
+
+Release 2.8-r672 (1 February 2018)
+----------------------------------
+
+Notable changes in this release include:
+
+ * Speed up short-read alignment by ~10%. The overall mapping accuracy stays
+   the same, but the output alignments are not always identical to v2.7 due to
+   unstable sorting employed during chaining. Long-read alignment is not
+   affected by this change as the speedup is short-read specific.
+
+ * Mappy now supports paired-end short-read alignment (#87). Please see
+   python/README.rst for details.
+
+ * Added option --for-only and --rev-only to perform alignment against the
+   forward or the reverse strand of the reference genome only (#91).
+
+ * Alleviated the issue with undesired diagonal alignment in the self mapping
+   mode (#10). Even if the output is not ideal, it should not interfere with
+   other alignments. Fully resolving the issue is intricate and may require
+   additional heuristic thresholds.
+
+ * Enhanced error checking against incorrect input (#92 and #96).
+
+For long query sequences, minimap2 should output identical alignments to v2.7.
+
+(2.8: 1 February 2018, r672)
+
+
+
 Release 2.7-r654 (9 January 2018)
 ---------------------------------

@@ -38,7 +38,7 @@ man ./minimap2.1
  - [Advanced features](#advanced)
    - [Working with >65535 CIGAR operations](#long-cigar)
    - [The cs optional tag](#cs)
-    - [Evaluation scripts](#eval)
+    - [Working with the PAF format](#paftools)
  - [Algorithm overview](#algo)
  - [Getting help](#help)
  - [Citing minimap2](#cite)
@@ -68,9 +68,9 @@ Detailed evaluations are available from the [minimap2 preprint][preprint].
 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.7/minimap2-2.7_x64-linux.tar.bz2 \
+curl -L https://github.com/lh3/minimap2/releases/download/v2.9/minimap2-2.9_x64-linux.tar.bz2 \
  | tar -jxvf -
-./minimap2-2.7_x64-linux/minimap2
+./minimap2-2.9_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
@@ -256,26 +256,13 @@ the alignment. The above example will become
 `=CGATCG-ata=AATAGAGTAG+gtc=GAAT*at=GCA`. The long form of `cs` encodes both
 reference and query sequences in one string.

-#### <a name="eval"></a>Evaluation scripts
+#### <a name="paftools"></a>Working with the PAF format

-Minimap2 comes with several (java)scripts for evaluating the accuracy of
-minimap2. These scripts require the [k8][k8] javascript shell to run.
-Recent minimap2 binary release tar-balls contain a copy of k8 executable, a
-single file. Here are a few examples on how to use these scripts:
-
-```sh
-# Generate reads from PBSIM alignment (truth encoded in read names)
-k8 misc/sim-pbsim.js ref.fa.fai pbsim-aln.maf > pbsim-reads.fq
-# Generate reads from mason2 alignment (not tested for simulated SVs)
-k8 misc/sim-mason2.js mason2-aln.sam > mason2-reads.fq
-# Evaluate mapping accuracy with ROC-like curve
-k8 misc/sim-eval.js my-aln.sam.gz > result.txt
-k8 misc/sim-eval.js my-aln.paf.gz > result.txt
-# Collect alignment statistics
-k8 misc/mapstat.js my-aln.sam > result.txt
-# Compare spliced junctions to existing gene annotations
-k8 misc/intron-eval.js anno.gtf my-spliced-aln.sam > result.txt
-```
+Minimap2 also comes with a (java)script [paftools.js](misc/paftools.js) that
+processes alignments in the PAF format. It calls variants from
+assembly-to-reference alignment, lifts over BED files based on alignment,
+converts between formats and provides utilities for various evaluations. For
+details, please see [misc/README.md](misc/README.md).

 ### <a name="algo"></a>Algorithm overview

@@ -359,10 +346,6 @@ mappy` or [from BioConda][mappyconda] via `conda install -c bioconda mappy`.
  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
-warmly welcomed.
-


 [paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
@@ -28,39 +28,63 @@ static inline void mm_seq_rev(uint32_t len, uint8_t *seq)
 		t = seq[i], seq[i] = seq[len - 1 - i], seq[len - 1 - i] = t;
 }

-static inline int test_zdrop_aux(int32_t score, int i, int j, int32_t *max, int *max_i, int *max_j, int e, int zdrop)
+static inline void update_max_zdrop(int32_t score, int i, int j, int32_t *max, int *max_i, int *max_j, int e, int *max_zdrop, int pos[2][2])
 {
 	if (score < *max) {
 		int li = i - *max_i;
 		int lj = j - *max_j;
 		int diff = li > lj? li - lj : lj - li;
-		if (*max - score > zdrop + diff * e)
-			return 1;
+		int z = *max - score - diff * e;
+		if (z > *max_zdrop) {
+			*max_zdrop = z;
+			pos[0][0] = *max_i, pos[0][1] = i + 1;
+			pos[1][0] = *max_j, pos[1][1] = j + 1;
+		}
 	} else *max = score, *max_i = i, *max_j = j;
-	return 0;
 }

-static int mm_check_zdrop(const uint8_t *qseq, const uint8_t *tseq, uint32_t n_cigar, uint32_t *cigar, const int8_t *mat, int8_t q, int8_t e, int zdrop)
+static int mm_test_zdrop(void *km, const mm_mapopt_t *opt, const uint8_t *qseq, const uint8_t *tseq, uint32_t n_cigar, uint32_t *cigar, const int8_t *mat)
 {
 	uint32_t k;
-	int32_t score = 0, max = 0, max_i = -1, max_j = -1, i = 0, j = 0;
-	for (k = 0; k < n_cigar; ++k) {
+	int32_t score = 0, max = INT32_MIN, max_i = -1, max_j = -1, i = 0, j = 0, max_zdrop = 0;
+	int pos[2][2] = {{-1, -1}, {-1, -1}}, q_len, t_len;
+
+	// find the score and the region where score drops most along diagonal
+	for (k = 0, score = 0; k < n_cigar; ++k) {
 		uint32_t l, op = cigar[k]&0xf, len = cigar[k]>>4;
 		if (op == 0) {
 			for (l = 0; l < len; ++l) {
 				score += mat[tseq[i + l] * 5 + qseq[j + l]];
-				if (test_zdrop_aux(score, i+l, j+l, &max, &max_i, &max_j, e, zdrop)) return 1;
+				update_max_zdrop(score, i+l, j+l, &max, &max_i, &max_j, opt->e, &max_zdrop, pos);
 			}
 			i += len, j += len;
-		} else if (op == 1) {
-			score -= q + e * len, j += len;
-			if (test_zdrop_aux(score, i, j, &max, &max_i, &max_j, e, zdrop)) return 1;
-		} else if (op == 2 || op == 3) {
-			score -= q + e * len, i += len;
-			if (test_zdrop_aux(score, i, j, &max, &max_i, &max_j, e, zdrop)) return 1;
+		} else if (op == 1 || op == 2 || op == 3) {
+			score -= opt->q + opt->e * len;
+			if (op == 1) j += len; // insertion
+			else i += len;         // deletion
+			update_max_zdrop(score, i, j, &max, &max_i, &max_j, opt->e, &max_zdrop, pos);
 		}
 	}
-	return 0;
+
+	// test if there is an inversion in the most dropped region
+	q_len = pos[1][1] - pos[1][0], t_len = pos[0][1] - pos[0][0];
+	if (!(opt->flag&(MM_F_SPLICE|MM_F_SR|MM_F_FOR_ONLY|MM_F_REV_ONLY)) && max_zdrop > opt->zdrop_inv && q_len < opt->max_gap && t_len < opt->max_gap) {
+		uint8_t *qseq2;
+		void *qp;
+		int q_off, t_off;
+		qseq2 = (uint8_t*)kmalloc(km, q_len);
+		for (i = 0; i < q_len; ++i) {
+			int c = qseq[pos[1][1] - i - 1];
+			qseq2[i] = c >= 4? 4 : 3 - c;
+		}
+		qp = ksw_ll_qinit(km, 2, q_len, qseq2, 5, mat);
+		score = ksw_ll_i16(qp, t_len, tseq + pos[0][0], opt->q, opt->e, &q_off, &t_off);
+		kfree(km, qseq2);
+		kfree(km, qp);
+		if (score >= opt->min_chain_score * opt->a && score >= opt->min_dp_max)
+			return 2; // there is a potential inversion
+	}
+	return max_zdrop > opt->zdrop? 1 : 0;
 }

 static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int *qshift, int *tshift)
@@ -193,9 +217,8 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) //
 	}
 }

-static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const int8_t *mat, int w, int end_bonus, int flag, ksw_extz_t *ez)
+static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
 {
-	int zdrop = opt->zdrop;
 	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
 		int i;
 		fprintf(stderr, "===> q=(%d,%d), e=(%d,%d), bw=%d, flag=%d, zdrop=%d <===\n", opt->q, opt->q2, opt->e, opt->e2, w, flag, opt->zdrop);
@@ -291,33 +314,39 @@ static void mm_filter_bad_seeds(void *km, int as1, int cnt1, mm128_t *a, int min
 	kfree(km, K);
 }

-static void mm_fix_bad_ends(const mm_reg1_t *r, const mm128_t *a, int bw, int32_t *as, int32_t *cnt)
+static void mm_fix_bad_ends(const mm_reg1_t *r, const mm128_t *a, int bw, int min_match, int32_t *as, int32_t *cnt)
 {
-	int32_t i, l;
+	int32_t i, l, m;
 	*as = r->as, *cnt = r->cnt;
 	if (r->cnt < 3) return;
-	l = a[r->as].y >> 32 & 0xff;
+	m = l = a[r->as].y >> 32 & 0xff;
 	for (i = r->as + 1; i < r->as + r->cnt - 1; ++i) {
 		int32_t lq, lr, min, max;
+		int32_t q_span = a[i].y >> 32 & 0xff;
+		if (a[i].y & MM_SEED_LONG_JOIN) break;
 		lr = (int32_t)a[i].x - (int32_t)a[i-1].x;
 		lq = (int32_t)a[i].y - (int32_t)a[i-1].y;
 		min = lr < lq? lr : lq;
 		max = lr > lq? lr : lq;
 		if (max - min > l >> 1) *as = i;
 		l += min;
-		if (l >= bw << 1) break;
+		m += min < q_span? min : q_span;
+		if (l >= bw << 1 || (m >= min_match && m >= bw) || m >= r->mlen >> 1) break;
 	}
 	*cnt = r->as + r->cnt - *as;
-	l = a[r->as + r->cnt - 1].y >> 32 & 0xff;
+	m = l = a[r->as + r->cnt - 1].y >> 32 & 0xff;
 	for (i = r->as + r->cnt - 2; i > *as; --i) {
 		int32_t lq, lr, min, max;
+		int32_t q_span = a[i+1].y >> 32 & 0xff;
+		if (a[i+1].y & MM_SEED_LONG_JOIN) break;
 		lr = (int32_t)a[i+1].x - (int32_t)a[i].x;
 		lq = (int32_t)a[i+1].y - (int32_t)a[i].y;
 		min = lr < lq? lr : lq;
 		max = lr > lq? lr : lq;
 		if (max - min > l >> 1) *cnt = i + 1 - *as;
 		l += min;
-		if (l >= bw) break;
+		m += min < q_span? min : q_span;
+		if (l >= bw << 1 || (m >= min_match && m >= bw) || m >= r->mlen >> 1) break;
 	}
 }

@@ -418,7 +447,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 		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);
+			mm_fix_bad_ends(r, a, opt->bw, opt->min_chain_score * 2, &as1, &cnt1);
 		}
 		mm_filter_bad_seeds(km, as1, cnt1, a, 10, 40, opt->max_gap>>1, 10);
 		mm_adjust_minier(mi, qseq0, &a[as1], &rs, &qs);
@@ -499,6 +528,14 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			re0 = re0 > re1? re0 : re1;
 		} else re0 = re, qe0 = qe;
 	}
+	if (a[r->as].y & MM_SEED_SELF) {
+		int max_ext = r->qs > r->rs? r->qs - r->rs : r->rs - r->qs;
+		if (r->rs - rs0 > max_ext) rs0 = r->rs - max_ext;
+		if (r->qs - qs0 > max_ext) qs0 = r->qs - max_ext;
+		max_ext = r->qe > r->re? r->qe - r->re : r->re - r->qe;
+		if (re0 - r->re > max_ext) re0 = r->re + max_ext;
+		if (qe0 - r->qe > max_ext) qe0 = r->qe + max_ext;
+	}

 	assert(re0 > rs0);
 	tseq = (uint8_t*)kmalloc(km, re0 - rs0);
@@ -508,7 +545,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 		mm_idx_getseq(mi, rid, rs0, rs, tseq);
 		mm_seq_rev(qs - qs0, qseq);
 		mm_seq_rev(rs - rs0, tseq);
-		mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, mat, bw, opt->end_bonus, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez);
+		mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, mat, bw, opt->end_bonus, r->split_inv? opt->zdrop_inv : opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez);
 		if (ez->n_cigar > 0) {
 			mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			r->p->dp_score += ez->max;
@@ -528,9 +565,10 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 		} else mm_adjust_minier(mi, qseq0, &a[as1 + i], &re, &qe);
 		re1 = re, qe1 = qe;
 		if (i == cnt1 - 1 || (a[as1+i].y&MM_SEED_LONG_JOIN) || (qe - qs >= opt->min_ksw_len && re - rs >= opt->min_ksw_len)) {
-			int j, bw1 = bw;
+			int j, bw1 = bw, zdrop_code;
 			if (a[as1+i].y & MM_SEED_LONG_JOIN)
 				bw1 = qe - qs > re - rs? qe - qs : re - rs;
+			// perform alignment
 			qseq = &qseq0[rev][qs];
 			mm_idx_getseq(mi, rid, rs, re, tseq);
 			if (is_sr) { // perform ungapped alignment
@@ -542,10 +580,12 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				}
 				ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, 0, qe - qs);
 			} else { // perform normal gapped alignment
-				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
+				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
 			}
-			if (mm_check_zdrop(qseq, tseq, ez->n_cigar, ez->cigar, mat, opt->q, opt->e, opt->zdrop))
-				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, extra_flag, ez); // second pass: lift approximate
+			// test Z-drop and inversion Z-drop
+			if ((zdrop_code = mm_test_zdrop(km, opt, qseq, tseq, ez->n_cigar, ez->cigar, mat)) != 0)
+				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, zdrop_code == 2? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate
+			// update CIGAR
 			if (ez->n_cigar > 0)
 				mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			if (ez->zdropped) { // truncated by Z-drop; TODO: sometimes Z-drop kicks in because the next seed placement is wrong. This can be fixed in principle.
@@ -557,8 +597,10 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				r->p->dp_score += ez->max;
 				re1 = rs + (ez->max_t + 1);
 				qe1 = qs + (ez->max_q + 1);
-				if (cnt1 - (j + 1) >= opt->min_cnt)
+				if (cnt1 - (j + 1) >= opt->min_cnt) {
 					mm_split_reg(r, r2, as1 + j + 1 - r->as, qlen, a);
+					if (zdrop_code == 2) r2->split_inv = 1;
+				}
 				break;
 			} else r->p->dp_score += ez->score;
 			rs = re, qs = qe;
@@ -568,7 +610,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	if (!dropped && qe < qe0 && re < re0) { // right extension
 		qseq = &qseq0[rev][qe];
 		mm_idx_getseq(mi, rid, re, re0, tseq);
-		mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, mat, bw, opt->end_bonus, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
+		mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
 		if (ez->n_cigar > 0) {
 			mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			r->p->dp_score += ez->max;
@@ -605,7 +647,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	if (r1->id != r1->parent && r1->parent != MM_PARENT_TMP_PRI) return 0;
 	if (r2->id != r2->parent && r2->parent != MM_PARENT_TMP_PRI) return 0;
 	if (r1->rid != r2->rid || r1->rev != r2->rev) return 0;
-	ql = r2->qs - r1->qe;
+	ql = r1->rev? r1->qs - r2->qe : r2->qs - r1->qe;
 	tl = r2->rs - r1->re;
 	if (ql < opt->min_chain_score || ql > opt->max_gap) return 0;
 	if (tl < opt->min_chain_score || tl > opt->max_gap) return 0;
@@ -613,7 +655,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	ksw_gen_simple_mat(5, mat, opt->a, opt->b);
 	tseq = (uint8_t*)kmalloc(km, tl);
 	mm_idx_getseq(mi, r1->rid, r1->re, r2->rs, tseq);
-	qseq = &qseq0[!r1->rev][qlen - r2->qs];
+	qseq = r1->rev? &qseq0[0][r2->qe] : &qseq0[1][qlen - r2->qs];

 	mm_seq_rev(ql, qseq);
 	mm_seq_rev(tl, tseq);
@@ -624,7 +666,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	mm_seq_rev(tl, tseq);
 	if (score < opt->min_dp_max) goto end_align1_inv;
 	q_off = ql - (q_off + 1), t_off = tl - (t_off + 1);
-	mm_align_pair(km, opt, ql - q_off, qseq + q_off, tl - t_off, tseq + t_off, mat, (int)(opt->bw * 1.5), -1, KSW_EZ_EXTZ_ONLY, ez);
+	mm_align_pair(km, opt, ql - q_off, qseq + q_off, tl - t_off, tseq + t_off, mat, (int)(opt->bw * 1.5), -1, opt->zdrop, KSW_EZ_EXTZ_ONLY, ez);
 	if (ez->n_cigar == 0) goto end_align1_inv; // should never be here
 	mm_append_cigar(r_inv, ez->n_cigar, ez->cigar);
 	r_inv->p->dp_score = ez->max;
@@ -634,8 +676,15 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	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;
+	if (r_inv->rev == 0) {
+		r_inv->qs = r2->qe + q_off;
+		r_inv->qe = r_inv->qs + ez->max_q + 1;
+	} else {
+		r_inv->qe = r2->qs - q_off;
+		r_inv->qs = r_inv->qe - (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], &tseq[t_off], mat, opt->q, opt->e);
 	ret = 1;
 end_align1_inv:
@@ -696,7 +745,7 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 				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 (i > 0 && regs[i].split_inv) {
 			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
@@ -706,7 +755,7 @@ 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]);
 	kfree(km, ez.cigar);
-	mm_filter_regs(km, opt, n_regs_, regs);
+	mm_filter_regs(km, opt, qlen, n_regs_, regs);
 	mm_hit_sort_by_dp(km, n_regs_, regs);
 	return regs;
 }
@@ -2,6 +2,7 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <assert.h>
+#define __STDC_LIMIT_MACROS
 #include "bseq.h"
 #include "kvec.h"
 #include "kseq.h"
@@ -53,15 +54,23 @@ void mm_bseq_close(mm_bseq_file_t *fp)
 	free(fp);
 }

+static inline char *kstrdup(const kstring_t *s)
+{
+	char *t;
+	t = (char*)malloc(s->l + 1);
+	memcpy(t, s->s, s->l + 1);
+	return t;
+}
+
 static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual)
 {
 	int i;
-	s->name = strdup(ks->name.s);
-	s->seq = strdup(ks->seq.s);
+	s->name = kstrdup(&ks->name);
+	s->seq = kstrdup(&ks->seq);
 	for (i = 0; i < ks->seq.l; ++i) // convert U to T
 		if (s->seq[i] == 'u' || s->seq[i] == 'U')
 			--s->seq[i];
-	s->qual = with_qual && ks->qual.l? strdup(ks->qual.s) : 0;
+	s->qual = with_qual && ks->qual.l? kstrdup(&ks->qual) : 0;
 	s->l_seq = ks->seq.l;
 }

@@ -114,10 +123,15 @@ mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int
 	*n_ = 0;
 	if (n_fp < 1) return 0;
 	while (1) {
+		int n_read = 0;
 		for (i = 0; i < n_fp; ++i)
-			if (kseq_read(fp[i]->ks) < 0)
-				break;
-		if (i != n_fp) break; // some file reaches the end
+			if (kseq_read(fp[i]->ks) >= 0)
+				++n_read;
+		if (n_read < n_fp) {
+			if (n_read > 0)
+				fprintf(stderr, "[W::%s]\033[1;31m query files have different number of records; extra records skipped.\033[0m\n", __func__);
+			break; // some file reaches the end
+		}
 		if (a.m == 0) kv_resize(mm_bseq1_t, 0, a, 256);
 		for (i = 0; i < n_fp; ++i) {
 			mm_bseq1_t *s;
@@ -94,6 +94,7 @@ void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a)
 	r2->id = -1;
 	r2->sam_pri = 0;
 	r2->p = 0;
+	r2->split_inv = 0;
 	r2->cnt = r->cnt - n;
 	r2->score = (int32_t)(r->score * ((float)r2->cnt / r->cnt) + .499);
 	r2->as = r->as + n;
@@ -245,16 +246,17 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_,
 	}
 }

-void mm_filter_regs(void *km, const mm_mapopt_t *opt, int *n_regs, mm_reg1_t *regs)
+void mm_filter_regs(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs)
 { // NB: after this call, mm_reg1_t::parent can be -1 if its parent filtered out
 	int i, k;
 	for (i = k = 0; i < *n_regs; ++i) {
 		mm_reg1_t *r = &regs[i];
 		int flt = 0;
 		if (!r->inv && !r->seg_split && r->cnt < opt->min_cnt) flt = 1;
-		if (r->p) {
+		if (r->p) { // these filters are only applied when base-alignment is available
 			if (r->mlen < opt->min_chain_score) flt = 1;
 			else if (r->p->dp_max < opt->min_dp_max) flt = 1;
+			else if (r->qs > qlen * opt->max_clip_ratio && qlen - r->qe > qlen * opt->max_clip_ratio) flt = 1;
 			if (flt) free(r->p);
 		}
 		if (!flt) {
@@ -337,7 +339,7 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
 					r->parent = regs[r->parent].parent;
 			}
 		}
-		mm_filter_regs(km, opt, n_regs_, regs);
+		mm_filter_regs(km, opt, qlen, n_regs_, regs);
 		mm_sync_regs(km, *n_regs_, regs);
 	}
 }
@@ -390,8 +392,10 @@ mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int
 	for (s = 0; s < n_segs; ++s) {
 		regs[s] = mm_gen_regs(km, hash, qlens[s], seg[s].n_u, seg[s].u, seg[s].a);
 		n_regs[s] = seg[s].n_u;
-		for (i = 0; i < n_regs[s]; ++i)
+		for (i = 0; i < n_regs[s]; ++i) {
 			regs[s][i].seg_split = 1;
+			regs[s][i].seg_id = s;
+		}
 	}
 	return seg;
 }
@@ -404,7 +408,33 @@ void mm_seg_free(void *km, int n_segs, mm_seg_t *segs)
 	kfree(km, segs);
 }

-void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr)
+static void mm_set_inv_mapq(void *km, int n_regs, mm_reg1_t *regs)
+{
+	int i, n_aux;
+	uint64_t *aux;
+	if (n_regs < 3) return;
+	for (i = 0; i < n_regs; ++i)
+		if (regs[i].inv) break;
+	if (i == n_regs) return; // no inversion hits
+
+	aux = (uint64_t*)kmalloc(km, n_regs * 8);
+	for (i = n_aux = 0; i < n_regs; ++i)
+		if (regs[i].parent == i || regs[i].parent < 0)
+			aux[n_aux++] = (uint64_t)regs[i].as << 32 | i;
+	radix_sort_64(aux, aux + n_aux);
+
+	for (i = 1; i < n_aux - 1; ++i) {
+		mm_reg1_t *inv = &regs[(int32_t)aux[i]];
+		if (inv->inv) {
+			mm_reg1_t *l = &regs[(int32_t)aux[i-1]];
+			mm_reg1_t *r = &regs[(int32_t)aux[i+1]];
+			inv->mapq = l->mapq < r->mapq? l->mapq : r->mapq;
+		}
+	}
+	kfree(km, aux);
+}
+
+void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr)
 {
 	static const float q_coef = 40.0f;
 	int64_t sum_sc = 0;
@@ -447,4 +477,5 @@ void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, in
 			if (r->p && r->p->dp_max > r->p->dp_max2 && r->mapq == 0) r->mapq = 1;
 		} else r->mapq = 0;
 	}
+	mm_set_inv_mapq(km, n_regs, regs);
 }
@@ -7,6 +7,7 @@
 #endif
 #include <fcntl.h>
 #include <stdio.h>
+#define __STDC_LIMIT_MACROS
 #include "kthread.h"
 #include "bseq.h"
 #include "minimap.h"
@@ -19,6 +20,8 @@
 KHASH_INIT(idx, uint64_t, uint64_t, 1, idx_hash, idx_eq)
 typedef khash_t(idx) idxhash_t;

+KHASH_MAP_INIT_STR(str, uint32_t)
+
 #define kroundup64(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, (x)|=(x)>>32, ++(x))

 typedef struct mm_idx_bucket_s {
@@ -28,15 +31,6 @@ typedef struct mm_idx_bucket_s {
 	void *h;     // hash table indexing _p_ and minimizers appearing once
 } mm_idx_bucket_t;

-void mm_idxopt_init(mm_idxopt_t *opt)
-{
-	memset(opt, 0, sizeof(mm_idxopt_t));
-	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 flag)
 {
 	mm_idx_t *mi;
@@ -53,6 +47,7 @@ void mm_idx_destroy(mm_idx_t *mi)
 {
 	int i;
 	if (mi == 0) return;
+	if (mi->h) kh_destroy(str, (khash_t(str)*)mi->h);
 	for (i = 0; i < 1<<mi->b; ++i) {
 		free(mi->B[i].p);
 		free(mi->B[i].a.a);
@@ -108,6 +103,34 @@ void mm_idx_stat(const mm_idx_t *mi)
 			__func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), n, 100.0*n1/n, (double)sum / n, (double)len / sum);
 }

+int mm_idx_index_name(mm_idx_t *mi)
+{
+	khash_t(str) *h;
+	uint32_t i;
+	int has_dup = 0, absent;
+	if (mi->h) return 0;
+	h = kh_init(str);
+	for (i = 0; i < mi->n_seq; ++i) {
+		khint_t k;
+		k = kh_put(str, h, mi->seq[i].name, &absent);
+		if (absent) kh_val(h, k) = i;
+		else has_dup = 1;
+	}
+	mi->h = h;
+	if (has_dup && mm_verbose >= 2)
+		fprintf(stderr, "[WARNING] some database sequences have identical sequence names\n");
+	return has_dup;
+}
+
+int mm_idx_name2id(const mm_idx_t *mi, const char *name)
+{
+	khash_t(str) *h = (khash_t(str)*)mi->h;
+	khint_t k;
+	if (h == 0) return -2;
+	k = kh_get(str, h, name);
+	return k == kh_end(h)? -1 : kh_val(h, k);
+}
+
 int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq)
 {
 	uint64_t i, st1, en1;
@@ -39,7 +39,24 @@ typedef struct {

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

-#define KSORT_INIT(name, type_t, __sort_lt)								\
+#define KSORT_INIT(name, type_t, __sort_lt) \
+	void ks_heapdown_##name(size_t i, size_t n, type_t l[]) \
+	{ \
+		size_t k = i; \
+		type_t tmp = l[i]; \
+		while ((k = (k << 1) + 1) < n) { \
+			if (k != n - 1 && __sort_lt(l[k], l[k+1])) ++k; \
+			if (__sort_lt(l[k], tmp)) break; \
+			l[i] = l[k]; i = k; \
+		} \
+		l[i] = tmp; \
+	} \
+	void ks_heapmake_##name(size_t lsize, type_t l[]) \
+	{ \
+		size_t i; \
+		for (i = (lsize >> 1) - 1; i != (size_t)(-1); --i) \
+			ks_heapdown_##name(i, lsize, l); \
+	} \
 	type_t ks_ksmall_##name(size_t n, type_t arr[], size_t kk)			\
 	{																	\
 		type_t *low, *high, *k, *ll, *hh, *mid;							\
@@ -2,6 +2,7 @@
 #include <stdlib.h>
 #include <limits.h>
 #include <stdint.h>
+#include "kthread.h"

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

-#define MM_VERSION "2.7-r654"
+#define MM_VERSION "2.9-r720"

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

+static inline void yes_or_no(mm_mapopt_t *opt, int flag, int long_idx, const char *arg, int yes_to_set)
+{
+	if (yes_to_set) {
+		if (strcmp(arg, "yes") == 0 || strcmp(arg, "y") == 0) opt->flag |= flag;
+		else if (strcmp(arg, "no") == 0 || strcmp(arg, "n") == 0) opt->flag &= ~flag;
+		else fprintf(stderr, "[WARNING]\033[1;31m option '--%s' only accepts 'yes' or 'no'.\033[0m\n", long_options[long_idx].name);
+	} else {
+		if (strcmp(arg, "yes") == 0 || strcmp(arg, "y") == 0) opt->flag &= ~flag;
+		else if (strcmp(arg, "no") == 0 || strcmp(arg, "n") == 0) opt->flag |= flag;
+		else fprintf(stderr, "[WARNING]\033[1;31m option '--%s' only accepts 'yes' or 'no'.\033[0m\n", long_options[long_idx].name);
+	}
+}
+
 int main(int argc, char *argv[])
 {
-	const char *opt_str = "2aSw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:";
+	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:";
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
 	int i, c, n_threads = 3, long_idx;
@@ -108,7 +127,9 @@ int main(int argc, char *argv[])
 		else if (c == 'p') opt.pri_ratio = atof(optarg);
 		else if (c == 'M') opt.mask_level = atof(optarg);
 		else if (c == 'c') opt.flag |= MM_F_OUT_CG | MM_F_CIGAR;
-		else if (c == 'X') opt.flag |= MM_F_AVA | MM_F_NO_SELF;
+		else if (c == 'D') opt.flag |= MM_F_NO_DIAG;
+		else if (c == 'P') opt.flag |= MM_F_ALL_CHAINS;
+		else if (c == 'X') opt.flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN; // -D -P --no-long-join --dual=no
 		else if (c == 'a') opt.flag |= MM_F_OUT_SAM | MM_F_CIGAR;
 		else if (c == 'Q') opt.flag |= MM_F_NO_QUAL;
 		else if (c == 'Y') opt.flag |= MM_F_SOFTCLIP;
@@ -118,7 +139,6 @@ int main(int argc, char *argv[])
 		else if (c == 'm') opt.min_chain_score = atoi(optarg);
 		else if (c == 'A') opt.a = atoi(optarg);
 		else if (c == 'B') opt.b = atoi(optarg);
-		else if (c == 'z') opt.zdrop = atoi(optarg);
 		else if (c == 's') opt.min_dp_max = atoi(optarg);
 		else if (c == 'C') opt.noncan = atoi(optarg);
 		else if (c == 'I') ipt.batch_size = mm_parse_num(optarg);
@@ -130,7 +150,6 @@ int main(int argc, char *argv[])
 		else if (c == 0 && long_idx == 2) opt.seed = atoi(optarg); // --seed
 		else if (c == 0 && long_idx == 3) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
 		else if (c == 0 && long_idx == 4) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
-		else if (c == 0 && long_idx == 5) opt.flag |= MM_F_NO_SELF; // --no-self
 		else if (c == 0 && long_idx == 6) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
 		else if (c == 0 && long_idx == 7) opt.max_chain_skip = atoi(optarg); // --max-chain-skip
 		else if (c == 0 && long_idx == 8) opt.min_ksw_len = atoi(optarg); // --min-dp-len
@@ -142,14 +161,13 @@ int main(int argc, char *argv[])
 		else if (c == 0 && long_idx ==18) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
 		else if (c == 0 && long_idx ==20) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
 		else if (c == 0 && long_idx ==21) opt.anchor_ext_shift = atoi(optarg); // --end-seed-pen
+		else if (c == 0 && long_idx ==22) opt.flag |= MM_F_FOR_ONLY; // --for-only
+		else if (c == 0 && long_idx ==23) opt.flag |= MM_F_REV_ONLY; // --rev-only
+		else if (c == 0 && long_idx ==27) opt.max_clip_ratio = atof(optarg); // --max-clip-ratio
 		else if (c == 0 && long_idx == 14) { // --frag
-			if (optarg == 0 || strcmp(optarg, "yes") == 0 || strcmp(optarg, "y") == 0)
-				opt.flag |= MM_F_FRAG_MODE;
-			else opt.flag &= ~MM_F_FRAG_MODE;
+			yes_or_no(&opt, MM_F_FRAG_MODE, long_idx, optarg, 1);
 		} else if (c == 0 && long_idx == 15) { // --secondary
-			if (optarg == 0 || strcmp(optarg, "yes") == 0 || strcmp(optarg, "y") == 0)
-				opt.flag &= ~MM_F_NO_PRINT_2ND;
-			else opt.flag |= MM_F_NO_PRINT_2ND;
+			yes_or_no(&opt, MM_F_NO_PRINT_2ND, long_idx, optarg, 0);
 		} else if (c == 0 && long_idx == 16) { // --cs
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR;
 			if (optarg == 0 || strcmp(optarg, "short") == 0) {
@@ -162,9 +180,11 @@ int main(int argc, char *argv[])
 				fprintf(stderr, "[WARNING]\033[1;31m --cs only takes 'short' or 'long'. Invalid values are assumed to be 'short'.\033[0m\n");
 			}
 		} else if (c == 0 && long_idx == 19) { // --splice-flank
-			if (optarg == 0 || strcmp(optarg, "yes") == 0 || strcmp(optarg, "y") == 0)
-				opt.flag |= MM_F_SPLICE_FLANK;
-			else opt.flag &= ~MM_F_SPLICE_FLANK;
+			yes_or_no(&opt, MM_F_SPLICE_FLANK, long_idx, optarg, 1);
+		} else if (c == 0 && long_idx == 24) { // --heap-sort
+			yes_or_no(&opt, MM_F_HEAP_SORT, long_idx, optarg, 1);
+		} else if (c == 0 && long_idx == 26) { // --dual
+			yes_or_no(&opt, MM_F_NO_DUAL, long_idx, optarg, 0);
 		} else if (c == 'S') {
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR | MM_F_OUT_CS_LONG;
 			if (mm_verbose >= 2)
@@ -188,6 +208,9 @@ int main(int argc, char *argv[])
 				fprintf(stderr, "[ERROR]\033[1;31m unrecognized cDNA direction\033[0m\n");
 				return 1;
 			}
+		} else if (c == 'z') {
+			opt.zdrop = opt.zdrop_inv = strtol(optarg, &s, 10);
+			if (*s == ',') opt.zdrop_inv = strtol(s + 1, &s, 10);
 		} else if (c == 'O') {
 			opt.q = opt.q2 = strtol(optarg, &s, 10);
 			if (*s == ',') opt.q2 = strtol(s + 1, &s, 10);
@@ -231,7 +254,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -B INT       mismatch penalty [%d]\n", opt.b);
 		fprintf(fp_help, "    -O INT[,INT] gap open penalty [%d,%d]\n", opt.q, opt.q2);
 		fprintf(fp_help, "    -E INT[,INT] gap extension penalty; a k-long gap costs min{O1+k*E1,O2+k*E2} [%d,%d]\n", opt.e, opt.e2);
-		fprintf(fp_help, "    -z INT       Z-drop score [%d]\n", opt.zdrop);
+		fprintf(fp_help, "    -z INT[,INT] Z-drop score and inversion Z-drop score [%d,%d]\n", opt.zdrop, opt.zdrop_inv);
 		fprintf(fp_help, "    -s INT       minimal peak DP alignment score [%d]\n", opt.min_dp_max);
 		fprintf(fp_help, "    -u CHAR      how to find GT-AG. f:transcript strand, b:both strands, n:don't match GT-AG [n]\n");
 		fprintf(fp_help, "  Input/Output:\n");
@@ -252,14 +275,18 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "                 map-ont: -k15 (Oxford Nanopore vs reference mapping)\n");
 		fprintf(fp_help, "                 asm5: -k19 -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200 (asm to ref mapping; break at 5%% div.)\n");
 		fprintf(fp_help, "                 asm10: -k19 -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200 (asm to ref mapping; break at 10%% div.)\n");
-		fprintf(fp_help, "                 ava-pb: -Hk19 -w5 -Xp0 -m100 -g10000 --max-chain-skip 25 (PacBio read overlap)\n");
-		fprintf(fp_help, "                 ava-ont: -k15 -w5 -Xp0 -m100 -g10000 --max-chain-skip 25 (ONT read overlap)\n");
+		fprintf(fp_help, "                 ava-pb: -Hk19 -Xw5 -m100 -g10000 --max-chain-skip 25 (PacBio read overlap)\n");
+		fprintf(fp_help, "                 ava-ont: -k15 -Xw5 -m100 -g10000 --max-chain-skip 25 (ONT read overlap)\n");
 		fprintf(fp_help, "                 splice: long-read spliced alignment (see minimap2.1 for details)\n");
 		fprintf(fp_help, "                 sr: short single-end reads without splicing (see minimap2.1 for details)\n");
 		fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of command-line options.\n");
 		return fp_help == stdout? 0 : 1;
 	}

+	if ((opt.flag & MM_F_SR) && argc - optind > 3) {
+		fprintf(stderr, "[ERROR] incorrect input: in the sr mode, please specify no more than two query files.\n");
+		return 1;
+	}
 	idx_rdr = mm_idx_reader_open(argv[optind], &ipt, fnw);
 	if (idx_rdr == 0) {
 		fprintf(stderr, "[ERROR] failed to open file '%s'\n", argv[optind]);
@@ -9,134 +9,6 @@
 #include "bseq.h"
 #include "khash.h"

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

-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;
+#include "ksort.h"
+#define heap_lt(a, b) ((a).x > (b).x)
+KSORT_INIT(heap, mm128_t, heap_lt)

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

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

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

-void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, const char **quals, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
+void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
 {
 	int i, j, rep_len, qlen_sum, n_regs0, 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);
@@ -307,8 +280,8 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	hash  = __ac_Wang_hash(hash);

 	collect_minimizers(b->km, opt, mi, n_segs, qlens, seqs, &mv);
-	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 (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
+	else a = collect_seed_hits(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);

 	if (mm_dbg_flag & MM_DBG_PRINT_SEED) {
 		fprintf(stderr, "RS\t%d\n", rep_len);
@@ -348,8 +321,8 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 			kfree(b->km, a);
 			kfree(b->km, u);
 			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);
+			if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
+			else a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 			a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 		}
 	}
@@ -366,8 +339,8 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	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);
-		mm_set_mapq(n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr);
+		regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], &n_regs0, regs0, a);
+		mm_set_mapq(b->km, n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr);
 		n_regs[0] = n_regs0, regs[0] = regs0;
 	} else { // multi-segment
 		mm_seg_t *seg;
@@ -375,8 +348,8 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		free(regs0);
 		for (i = 0; i < n_segs; ++i) {
 			mm_set_parent(b->km, opt->mask_level, n_regs[i], regs[i], opt->a * 2 + opt->b); // update mm_reg1_t::parent
-			regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], quals? quals[i] : 0, &n_regs[i], regs[i], seg[i].a);
-			mm_set_mapq(n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr);
+			regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a);
+			mm_set_mapq(b->km, n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr);
 		}
 		mm_seg_free(b->km, n_segs, seg);
 		if (n_segs == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR))
@@ -403,7 +376,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 mm_reg1_t *mm_map(const mm_idx_t *mi, int qlen, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
 {
 	mm_reg1_t *regs;
-	mm_map_frag(mi, 1, &qlen, &seq, 0, n_regs, &regs, b, opt, qname);
+	mm_map_frag(mi, 1, &qlen, &seq, n_regs, &regs, b, opt, qname);
 	return regs;
 }

@@ -431,11 +404,10 @@ typedef struct {
 static void worker_for(void *_data, long i, int tid) // kt_for() callback
 {
    step_t *s = (step_t*)_data;
-	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[MM_MAX_SEG], *quals[MM_MAX_SEG];
+	int qlens[MM_MAX_SEG], j, off = s->seg_off[i], pe_ori = s->p->opt->pe_ori;
+	const char *qseqs[MM_MAX_SEG];
 	mm_tbuf_t *b = s->buf[tid];
 	assert(s->n_seg[i] <= MM_MAX_SEG);
-	memset(quals, 0, sizeof(char*) * MM_MAX_SEG);
 	if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
 		fprintf(stderr, "QR\t%s\t%d\t%d\n", s->seq[off].name, tid, s->seq[off].l_seq);
 	for (j = 0; j < s->n_seg[i]; ++j) {
@@ -443,13 +415,12 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
 			mm_revcomp_bseq(&s->seq[off + j]);
 		qlens[j] = s->seq[off + j].l_seq;
 		qseqs[j] = s->seq[off + j].seq;
-		quals[j] = is_sr? s->seq[off + j].qual : 0;
 	}
 	if (s->p->opt->flag & MM_F_INDEPEND_SEG) {
 		for (j = 0; j < s->n_seg[i]; ++j)
-			mm_map_frag(s->p->mi, 1, &qlens[j], &qseqs[j], &quals[j], &s->n_reg[off+j], &s->reg[off+j], b, s->p->opt, s->seq[off+j].name);
+			mm_map_frag(s->p->mi, 1, &qlens[j], &qseqs[j], &s->n_reg[off+j], &s->reg[off+j], b, s->p->opt, s->seq[off+j].name);
 	} else {
-		mm_map_frag(s->p->mi, s->n_seg[i], qlens, qseqs, quals, &s->n_reg[off], &s->reg[off], b, s->p->opt, s->seq[off].name);
+		mm_map_frag(s->p->mi, s->n_seg[i], qlens, qseqs, &s->n_reg[off], &s->reg[off], b, s->p->opt, s->seq[off].name);
 	}
 	for (j = 0; j < s->n_seg[i]; ++j) // flip the query strand and coordinate to the original read strand
 		if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1)))) {
@@ -5,26 +5,30 @@
 #include <stdio.h>
 #include <sys/types.h>

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

 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2
@@ -55,7 +59,7 @@ typedef struct {
 	mm_idx_seq_t *seq;         // sequence name, length and offset
 	uint32_t *S;               // 4-bit packed sequence
 	struct mm_idx_bucket_s *B; // index (hidden)
-	void *km;
+	void *km, *h;
 } mm_idx_t;

 // minimap2 alignment
@@ -78,7 +82,7 @@ typedef struct {
 	int32_t mlen, blen;     // seeded exact match length; seeded alignment block length
 	int32_t n_sub;          // number of suboptimal mappings
 	int32_t score0;         // initial chaining score (before chain merging/spliting)
-	uint32_t mapq:8, split:2, rev:1, inv:1, sam_pri:1, proper_frag:1, pe_thru:1, seg_split:1, dummy:16;
+	uint32_t mapq:8, split:2, rev:1, inv:1, sam_pri:1, proper_frag:1, pe_thru:1, seg_split:1, seg_id:8, split_inv:1, dummy:7;
 	uint32_t hash;
 	float div;
 	mm_extra_t *p;
@@ -112,11 +116,12 @@ typedef struct {

 	int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties
 	int noncan;      // cost of non-canonical splicing sites
-	int zdrop;       // break alignment if alignment score drops too fast along the diagonal
+	int zdrop, zdrop_inv;   // break alignment if alignment score drops too fast along the diagonal
 	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;
+	float max_clip_ratio; // drop an alignment if BOTH ends are clipped above this ratio

 	int pe_ori, pe_bonus;

@@ -276,6 +281,8 @@ void mm_tbuf_destroy(mm_tbuf_t *b);
 */
 mm_reg1_t *mm_map(const mm_idx_t *mi, int l_seq, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *name);

+void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname);
+
 /**
 * Align a fasta/fastq file and print alignments to stdout
 *
@@ -290,6 +297,11 @@ int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int

 int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_mapopt_t *opt, int n_threads);

+// query sequence name and sequence in the minimap2 index
+int mm_idx_index_name(mm_idx_t *mi);
+int mm_idx_name2id(const mm_idx_t *mi, const char *name);
+int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
+
 // 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 flag, int n_threads);
@@ -1,4 +1,4 @@
-.TH minimap2 1 "9 January 2018" "minimap2-2.7 (r654)" "Bioinformatics tools"
+.TH minimap2 1 "24 February 2018" "minimap2-2.9 (r720)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -129,7 +129,7 @@ Ignore top
 fraction of most frequent minimizers [0.0002]
 .TP
 .BI -g \ INT
-Stop chain enlongation if there are no minimizers in
+Stop chain enlongation if there are no minimizers within
 .IR INT -bp
 [10000].
 .TP
@@ -148,11 +148,30 @@ Discard chains with chaining score
 [40]. Chaining score equals the approximate number of matching bases minus a
 concave gap penalty. It is computed with dynamic programming.
 .TP
+.B -D
+If query sequence name/length are identical to the target name/length, ignore
+diagonal anchors. This option also reduces DP-based extension along the
+diagonal.
+.TP
+.B -P
+Retain all chains and don't attempt to set primary chains. Options
+.B -p
+and
+.B -N
+have no effect when this option is in use.
+.TP
+.BR --dual = yes | no
+If
+.BR no ,
+skip query-target pairs wherein the query name is lexicographically greater
+than the target name [yes]
+.TP
 .B -X
-Perform all-vs-all mapping. In this mode, if the query sequence name is
-lexicographically larger than the target sequence name, the hits between them
-will be suppressed; if the query sequence name is the same as the target name,
-diagonal minimizer hits will also be suppressed.
+Equivalent to
+.RB ' -DP
+.BR --dual = no
+.BR --no-long-join '.
+Primarily used for all-vs-all read overlapping.
 .TP
 .BI -p \ FLOAT
 Minimal secondary-to-primary score ratio to output secondary mappings [0.8].
@@ -162,6 +181,9 @@ the chain with a lower score is secondary to the chain with a higher score.
 If the ratio of the scores is below
 .IR FLOAT ,
 the secondary chain will not be outputted or extended with DP alignment later.
+This option has no effect when
+.B -X
+is applied.
 .TP
 .BI -N \ INT
 Output at most
@@ -179,9 +201,14 @@ Increasing this option slows down spliced alignment. [200k]
 .TP
 .BI -F \ NUM
 Maximum fragment length (aka insert size; effective with
-.BR -xsr / --frag)
+.BR -xsr / --frag = yes )
 [800]
 .TP
+.BI -M \ FLOAT
+Mark as secondary a chain that overlaps with a better chain by
+.I FLOAT
+or more of the shorter chain [0.5]
+.TP
 .BI --max-chain-skip \ INT
 A heuristics that stops chaining early [50]. Minimap2 uses dynamic programming
 for chaining. The time complexity is quadratic in the number of seeds. This
@@ -204,8 +231,20 @@ applies a second round of chaining with a higher minimizer occurrence threshold
 if no good chain is found. In addition, minimap2 attempts to patch gaps between
 seeds with ungapped alignment.
 .TP
-.BR --frag [= no | yes ]
+.BR --frag = no | yes
 Whether to enable the fragment mode [no]
+.TP
+.B --for-only
+Only map to the forward strand of the reference sequences. For paired-end
+reads in the forward-reverse orientation, the first read is mapped to forward
+strand of the reference and the second read to the reverse stand.
+.TP
+.B --rev-only
+Only map to the reverse complement strand of the reference sequences.
+.TP
+.BR --heap-sort = no | yes
+If yes, sort anchors with heap merge, instead of radix sort. Heap merge is
+faster for short reads, but slower for long reads. [no]
 .SS Alignment options
 .TP 10
 .BI -A \ INT
@@ -232,11 +271,22 @@ Cost for a non-canonical GT-AG splicing (effective with
 .BR --splice )
 [0]
 .TP
-.BI -z \ INT
-Break an alignment if the running score drops too quickly along the diagonal of
-the DP matrix (diagonal X-drop, or Z-drop) [400]. Increasing the value improves
-the contiguity of the alignment at the cost of poor alignment in the middle
-(e.g. caused by a long inversion).
+.BI -z \ INT1[,INT2]
+Truncate an alignment if the running alignment score drops too quickly along
+the diagonal of the DP matrix (diagonal X-drop, or Z-drop) [400,200]. If the
+drop of score is above
+.IR INT2 ,
+minimap2 will reverse complement the query in the related region and align
+again to test small inversions. Minimap2 truncates alignment if there is an
+inversion or the drop of score is greater than
+.IR INT1 .
+Decrease
+.I INT2
+to find small inversions at the cost of performance and false positives.
+Increase
+.I INT1
+to improves the contiguity of alignment at the cost of poor alignment in the
+middle.
 .TP
 .BI -s \ INT
 Minimal peak DP alignment score to output [40]. The peak score is computed from
@@ -255,14 +305,13 @@ no attempt to match GT-AG [n]
 .BI --end-bonus \ INT
 Score bonus when alignment extends to the end of the query sequence [0].
 .TP
-.BR --splice-flank [= yes | no ]
+.BR --splice-flank = yes | no
 Assume the next base to a
 .B GT
 donor site tends to be A/G (91% in human and 92% in mouse) and the preceding
 base to a
 .B AG
-acceptor tends to be C/T [yes with
-.BR --splice ].
+acceptor tends to be C/T [no].
 This trend is evolutionarily conservative, all the way to S. cerevisiae
 (PMID:18688272). Specifying this option generally leads to higher junction
 accuracy by several percents, so it is applied by default with
@@ -352,7 +401,7 @@ K/M/G/k/m/g suffix is accepted. A large
 helps load balancing in the multi-threading mode, at the cost of increased
 memory.
 .TP
-.BR --secondary [= yes | no ]
+.BR --secondary = yes | no
 Whether to output secondary alignments [yes]
 .TP
 .B --version
@@ -399,13 +448,13 @@ Up to 10% sequence divergence.
 .B ava-pb
 PacBio all-vs-all overlap mapping
 .RB ( -Hk19
-.B -w5 -Xp0 -m100 -g10000 --max-chain-skip
+.B -Xw5 -m100 -g10000 --max-chain-skip
 .BR 25 ).
 .TP
 .B ava-ont
 Oxford Nanopore all-vs-all overlap mapping
 .RB ( -k15
-.B -w5 -Xp0 -m100 -g10000 --max-chain-skip
+.B -Xw5 -m100 -g10000 --max-chain-skip
 .BR 25 ).
 Similarly, the major difference from
 .B ava-pb
@@ -427,8 +476,8 @@ tag ignores introns to demote hits to pseudogenes.
 .B sr
 Short single-end reads without splicing
 .RB ( -k21
-.B -w11 --sr --frag -A2 -B8 -O12,32 -E2,1 -r50 -p.5 -N20 -f1000,5000 -n2 -m20
-.B -s40 -g200 -2K50m
+.B -w11 --sr --frag=yes -A2 -B8 -O12,32 -E2,1 -r50 -p.5 -N20 -f1000,5000 -n2 -m20
+.B -s40 -g200 -2K50m --heap-sort=yes
 .BR --secondary=no ).
 .RE
 .SS Miscellaneous options
@@ -522,8 +571,8 @@ where seed positions may be suboptimal. This should not be a big concern
 because even the optimal alignment may be wrong in such regions.
 .TP
 *
-Minimap2 requires SSE2 instructions to compile. It is possible to add
-non-SSE2 support, but it would make minimap2 slower by several times.
+Minimap2 requires SSE2 or NEON instructions to compile. It is possible to add
+non-SSE2/NEON support, but it would make minimap2 slower by several times.
 .SH SEE ALSO
 .PP
 miniasm(1), minimap(1), bwa(1).
@@ -1,4 +1,4 @@
-#include "minimap.h"
+#include "mmpriv.h"

 int mm_verbose = 1;
 int mm_dbg_flag = 0;
@@ -86,6 +86,8 @@ double cputime()

 	return kernelModeTime + userModeTime;
 }
+
+long peakrss(void) { return 0; }
 #else
 #include <sys/resource.h>
 #include <sys/time.h>
@@ -96,6 +98,18 @@ double cputime(void)
 	getrusage(RUSAGE_SELF, &r);
 	return r.ru_utime.tv_sec + r.ru_stime.tv_sec + 1e-6 * (r.ru_utime.tv_usec + r.ru_stime.tv_usec);
 }
+
+long peakrss(void)
+{
+	struct rusage r;
+	getrusage(RUSAGE_SELF, &r);
+#ifdef __linux__
+	return r.ru_maxrss * 1024;
+#else
+	return r.ru_maxrss;
+#endif
+}
+
 #endif /* WIN32 || _WIN32 */

 double realtime(void)
@@ -1,28 +1,179 @@
-The [K8 Javascript shell][k8] is needed to run Javascripts in this directory.
-Precompiled k8 binaries for Mac and Linux can be found at the [K8 release
-page][k8bin].
+## <a name="started"></a>Getting Started

-* [paf2aln.js](paf2aln.js): convert PAF to [MAF][maf] or BLAST-like output for
-  eyeballing. PAF has to be generated with minimap2 option `-S`, which writes
-  the aligned sequences to the `cs` tag. An example:
-  ```sh
-  ../minimap2 -S ../test/MT-*.fa | k8 paf2aln.js /dev/stdin
-  ```
+```sh
+# install minimap2
+git clone https://github.com/lh3/minimap2
+cd minimap2 && make
+# install the k8 javascript shell
+curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf -
+cp k8-0.2.4/k8-`uname -s` k8              # or copy it to a directory on your $PATH
+# export PATH="$PATH:`pwd`:`pwd`/misc"    # run this if k8, minimap2 or paftools.js not on your $PATH
+minimap2 --cs test/MT-human.fa test/MT-orang.fa | paftools.js view -     # view alignment
+minimap2 -c test/MT-human.fa test/MT-orang.fa | paftools.js stat -       # basic alignment statistics
+minimap2 -c --cs test/MT-human.fa test/MT-orang.fa \
+  | sort -k6,6 -k8,8n | paftools.js call -L15000 -     # calling variants from asm-to-ref alignment
+minimap2 -c test/MT-human.fa test/MT-orang.fa \
+  | paftools.js liftover -l10000 - <(echo -e "MT_orang\t2000\t5000")     # liftOver
+# no test data for the following examples
+paftools.js junceval -e anno.gtf splice.sam > out.txt  # compare splice junctions to annotations
+paftools.js splice2bed anno.gtf > anno.bed             # convert GTF/GFF3 to BED12
+```

-* [mapstat.js](mapstat.js): output basic statistics such as the number of
-  non-redundant mapped bases, number of split and secondary alignments and
-  number of long gaps. This scripts seamlessly works with both SAM and PAF.
+## Table of Contents

-* [sim-pbsim.js](sim-pbsim.js): convert reads simulated with [PBSIM][pbsim] to
-  FASTA and encode the true mapping positions to read names in a format like
-  `S1_33!chr1!225258409!225267761!-`.
+- [Getting Started](#started)
+- [Introduction](#intro)
+- [Evaluation](#eval)
+  - [Evaluating mapping accuracy with simulated reads](#mapeval)
+  - [Evaluating read overlap sensitivity](#oveval)
+- [Calling Variants from Assemblies](#asmvar)

-* [sim-eval.js](sim-eval.js): evaluate mapping accuracy for FASTA generated
-  with [sim-pbsim.js](sim-pbsim.js) or [sim-mason2.js](sim-mason2.js).
+## <a name="intro"></a>Introduction

-* [sam2paf.js](sam2paf.js): convert SAM to PAF.
+paftools.js is a script that processes alignments in the [PAF format][paf],
+such as converting between formats, evaluating mapping accuracy, lifting over
+BED files based on alignment, and calling variants from assembly-to-assembly
+alignment. This script *requires* the [k8 Javascript shell][k8] to run. On
+Linux or Mac, you can download the precompiled k8 binary with:

+```sh
+curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf -
+cp k8-0.2.4/k8-`uname -s` $HOME/bin/k8  # assuming $HOME/bin in your $PATH
+```
+
+It is highly recommended to copy the executable `k8` to a directory on your
+`$PATH` such as `/usr/bin/env` can find it. Like python scripts, once you
+install `k8`, you can launch paftools.js in one of the two ways:
+
+```sh
+path/to/paftools.js             # only if k8 is on your $PATH
+k8 path/to/paftools.js
+```
+
+In a nutshell, paftools.js has the following commands:
+
+```
+Usage: paftools.js <command> [arguments]
+Commands:
+  view       convert PAF to BLAST-like (for eyeballing) or MAF
+  splice2bed convert spliced alignment in PAF/SAM to BED12
+  sam2paf    convert SAM to PAF
+  delta2paf  convert MUMmer's delta to PAF
+  gff2bed    convert GTF/GFF3 to BED12
+
+  stat       collect basic mapping information in PAF/SAM
+  liftover   simplistic liftOver
+  call       call variants from asm-to-ref alignment with the cs tag
+  bedcov     compute the number of bases covered
+
+  mapeval    evaluate mapping accuracy using mason2/PBSIM-simulated FASTQ
+  mason2fq   convert mason2-simulated SAM to FASTQ
+  pbsim2fq   convert PBSIM-simulated MAF to FASTQ
+  junceval   evaluate splice junction consistency with known annotations
+  ov-eval    evaluate read overlap sensitivity using read-to-ref mapping
+```
+
+paftools.js seamlessly reads both plain text files and gzip'd text files.
+
+## <a name="eval"></a>Evaluation
+
+### <a name="mapeval"></a>Evaluating mapping accuracy with simulated reads
+
+The **pbsim2fq** command of paftools.js converts the MAF output of [pbsim][pbsim]
+to FASTQ and encodes the true mapping position in the read name in a format like
+`S1_33!chr1!225258409!225267761!-`. Similarly, the **mason2fq** command
+converts [mason2][mason2] simulated SAM to FASTQ.
+
+Command **mapeval** evaluates mapped SAM/PAF. Here is example output:
+
+```
+Q       60      32478   0       0.000000000     32478
+Q       22      16      1       0.000030775     32494
+Q       21      43      1       0.000061468     32537
+Q       19      73      1       0.000091996     32610
+Q       14      66      1       0.000122414     32676
+Q       10      27      3       0.000214048     32703
+Q       8       14      1       0.000244521     32717
+Q       7       13      2       0.000305530     32730
+Q       6       46      1       0.000335611     32776
+Q       3       10      1       0.000366010     32786
+Q       2       20      2       0.000426751     32806
+Q       1       248     94      0.003267381     33054
+Q       0       31      17      0.003778147     33085
+U       3
+```
+
+where each Q-line gives the quality threshold, the number of reads mapped with
+mapping quality equal to or greater than the threshold, number of wrong
+mappings, accumulative mapping error rate and the accumulative number of
+mapped reads. The U-line, if present, gives the number of unmapped reads if
+they are present in the SAM file.
+
+Suppose the reported mapping coordinate overlap with the true coordinate like
+the following:
+
+```
+truth:   --------------------
+mapper:           ----------------------
+         |<- l1 ->|<-- o -->|<-- l2 -->|
+```
+
+Let `r=o/(l1+o+l2)`. The reported mapping is considered correct if `r>0.1` by
+default.
+
+### <a name="oveval"></a>Evaluating read overlap sensitivity
+
+Command **ov-eval** takes *sorted* read-to-reference alignment and read
+overlaps in PAF as input, and evaluates the sensitivity. For example:
+
+```sh
+minimap2 -cx map-pb ref.fa reads.fq.gz | sort -k6,6 -k8,8n > reads-to-ref.paf
+minimap2 -x ava-pb reads.fq.gz reads.fq.gz > ovlp.paf
+k8 ov-eval.js reads-to-ref.paf ovlp.paf
+```
+
+## <a name="asmvar"></a>Calling Variants from Haploid Assemblies
+
+The **call** command of paftools.js calls variants from coordinate-sorted
+assembly-to-reference alignment. It calls variants from the [cs tag][cs] and
+identifies confident/callable regions as those covered by exactly one contig.
+Here are example command lines:
+
+```sh
+minimap2 -cx asm5 -t8 --cs ref.fa asm.fa > asm.paf  # keeping this file is recommended; --cs required!
+sort -k6,6 -k8,8n asm.paf > asm.srt.paf             # sort by reference start coordinate
+k8 paftools.js call asm.srt.paf > asm.var.txt
+```
+
+Here is sample output:
+
+```
+V   chr1    2276040 2276041 1   60  c   g   LJII01000171.1  1217409 1217410 +
+V   chr1    2280409 2280410 1   60  a   g   LJII01000171.1  1221778 1221779 +
+V   chr1    2280504 2280505 1   60  a   g   LJII01000171.1  1221873 1221874 +
+R   chr1    2325140 2436340
+V   chr1    2325287 2325287 1   60  -   ct  LJII01000171.1  1272894 1272896 +
+V   chr1    2325642 2325644 1   60  tt  -   LJII01000171.1  1273251 1273251 +
+V   chr1    2326051 2326052 1   60  c   t   LJII01000171.1  1273658 1273659 +
+V   chr1    2326287 2326288 1   60  c   t   LJII01000171.1  1273894 1273895 +
+```
+
+where a line starting with `R` gives regions covered by one query contig, and a
+V-line encodes a variant in the following format: chr, start, end, query depth,
+mapping quality, REF allele, ALT allele, query name, query start, end and the
+query orientation. Generally, you should only look at variants where column 5
+is one.
+
+By default, when calling variants, "paftools.js call" ignores alignments 50kb
+or shorter; when deriving callable regions, it ignores alignments 10kb or
+shorter.  It uses two thresholds to avoid edge effects. These defaults are
+designed for long-read assemblies. For short reads, both should be reduced.
+
+
+
+[paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
+[cs]: https://github.com/lh3/minimap2#cs
 [k8]: https://github.com/attractivechaos/k8
-[k8bin]: https://github.com/attractivechaos/k8/releases
 [maf]: https://genome.ucsc.edu/FAQ/FAQformat#format5
 [pbsim]: https://github.com/pfaucon/PBSIM-PacBio-Simulator
+[mason2]: https://github.com/seqan/seqan/tree/master/apps/mason2
@@ -1,258 +0,0 @@
-/*******************************
- * 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);
@@ -1,150 +0,0 @@
-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();
@@ -1,267 +0,0 @@
-/*******************************
- * 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.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 *
- *****************/
-
-var c, l_fuzzy = 0, print_ovlp = false, print_err_only = false, first_only = false;
-while ((c = getopt(arguments, "l:ep")) != null) {
-	if (c == 'l') l_fuzzy = parseInt(getopt.arg);
-	else if (c == 'e') print_err_only = print_ovlp = true;
-	else if (c == 'p') print_ovlp = true;
-}
-
-if (arguments.length - getopt.ind < 2) {
-	print("Usage: k8 intron-eval.js [options] <gene.gtf> <aln.sam>");
-	exit(1);
-}
-
-var file, buf = new Bytes();
-
-var tr = {};
-file = new File(arguments[getopt.ind]);
-while (file.readline(buf) >= 0) {
-	var m, t = buf.toString().split("\t");
-	if (t[0].charAt(0) == '#') continue;
-	if (t[2] != 'exon') continue;
-	var st = parseInt(t[3]) - 1;
-	var en = parseInt(t[4]);
-	if ((m = /transcript_id "(\S+)"/.exec(t[8])) == null) continue;
-	var tid = m[1];
-	if (tr[tid] == null) tr[tid] = [t[0], t[6], 0, 0, []];
-	tr[tid][4].push([st, en]);
-}
-file.close();
-
-var anno = {};
-for (var tid in tr) {
-	var t = tr[tid];
-	Interval.sort(t[4]);
-	t[2] = t[4][0][0];
-	t[3] = t[4][t[4].length - 1][1];
-	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])
-			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]]);
-	}
-}
-tr = null;
-
-for (var chr in anno) {
-	var e = anno[chr];
-	if (e.length == 0) continue;
-	Interval.sort(e);
-	var k = 0;
-	for (var i = 1; i < e.length; ++i) // dedup
-		if (e[i][0] != e[k][0] || e[i][1] != e[k][1])
-			e[++k] = e[i].slice(0);
-	e.length = k + 1;
-	Interval.index_end(e);
-}
-
-var n_pri = 0, n_unmapped = 0, n_mapped = 0;
-var n_sgl = 0, n_splice = 0, n_splice_hit = 0, n_splice_novel = 0;
-
-file = new File(arguments[getopt.ind+1]);
-var last_qname = null;
-var re_cigar = /(\d+)([MIDNSHX=])/g;
-while (file.readline(buf) >= 0) {
-	var m, t = buf.toString().split("\t");
-
-	if (t[0].charAt(0) == '@') continue;
-	var flag = parseInt(t[1]);
-	if (flag&0x100) continue;
-	if (first_only && last_qname == t[0]) continue;
-	if (t[2] == '*') {
-		++n_unmapped;
-		continue;
-	} else {
-		++n_pri;
-		if (last_qname != t[0]) ++n_mapped;
-	}
-
-	var pos = parseInt(t[3]) - 1, intron = [];
-	while ((m = re_cigar.exec(t[5])) != null) {
-		var len = parseInt(m[1]), op = m[2];
-		if (op == 'N') {
-			intron.push([pos, pos + len]);
-			pos += len;
-		} else if (op == 'M' || op == 'X' || op == '=' || op == 'D') pos += len;
-	}
-	if (intron.length == 0) {
-		++n_sgl;
-		continue;
-	}
-	n_splice += intron.length;
-
-	var chr = anno[t[2]];
-	if (chr != null) {
-		for (var i = 0; i < intron.length; ++i) {
-			var o = Interval.find_ovlp(chr, intron[i][0], intron[i][1]);
-			if (o.length > 0) {
-				var hit = false;
-				for (var j = 0; j < o.length; ++j) {
-					var st_diff = intron[i][0] - o[j][0];
-					var en_diff = intron[i][1] - o[j][1];
-					if (st_diff < 0) st_diff = -st_diff;
-					if (en_diff < 0) en_diff = -en_diff;
-					if (st_diff <= l_fuzzy && en_diff <= l_fuzzy)
-						++n_splice_hit, hit = true;
-					if (hit) break;
-				}
-				if (print_ovlp) {
-					var type = hit? 'C' : 'P';
-					if (hit && print_err_only) continue;
-					var x = '[';
-					for (var j = 0; j < o.length; ++j) {
-						if (j) x += ', ';
-						x += '(' + o[j][0] + "," + o[j][1] + ')';
-					}
-					x += ']';
-					print(type, t[0], i+1, t[2], intron[i][0], intron[i][1], x);
-				}
-			} else {
-				++n_splice_novel;
-				if (print_ovlp)
-					print('N', t[0], i+1, t[2], intron[i][0], intron[i][1]);
-			}
-		}
-	} else {
-		n_splice_novel += intron.length;
-	}
-	last_qname = t[0];
-}
-file.close();
-
-buf.destroy();
-
-if (!print_ovlp) {
-	print("# unmapped reads: " + n_unmapped);
-	print("# mapped reads: " + n_mapped);
-	print("# primary alignments: " + n_pri);
-	print("# singletons: " + n_sgl);
-	print("# predicted introns: " + n_splice);
-	print("# non-overlapping introns: " + n_splice_novel);
-	print("# correct introns: " + n_splice_hit + " (" + (n_splice_hit / n_splice * 100).toFixed(2) + "%)");
-}
@@ -1,183 +0,0 @@
-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, gap_out_len = null;
-while ((c = getopt(arguments, "l:")) != null)
-	if (c == 'l') gap_out_len = parseInt(getopt.arg);
-
-if (getopt.ind == arguments.length) {
-	print("Usage: k8 mapstat.js [-l gapOutLen] <in.sam>|<in.paf>");
-	exit(1);
-}
-
-var buf = new Bytes();
-var file = new File(arguments[getopt.ind]);
-var re = /(\d+)([MIDSHNX=])/g;
-
-var lineno = 0, n_pri = 0, n_2nd = 0, n_seq = 0, n_cigar_64k = 0, l_tot = 0, l_cov = 0;
-var n_gap = [[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]];
-
-function cov_len(regs)
-{
-	regs.sort(function(a,b) {return a[0]-b[0]});
-	var st = regs[0][0], en = regs[0][1], l = 0;
-	for (var i = 1; i < regs.length; ++i) {
-		if (regs[i][0] < en)
-			en = en > regs[i][1]? en : regs[i][1];
-		else l += en - st, st = regs[i][0], en = regs[i][1];
-	}
-	l += en - st;
-	return l;
-}
-
-var last = null, last_qlen = null, regs = [];
-while (file.readline(buf) >= 0) {
-	var line = buf.toString();
-	++lineno;
-	if (line.charAt(0) != '@') {
-		var t = line.split("\t", 12);
-		var m, rs, cigar = null, is_pri = false, is_sam = false, is_rev = false, tname = null;
-		var atlen = null, aqlen, qs, qe, mapq, ori_qlen;
-		if (t[4] == '+' || t[4] == '-') { // PAF
-			if (!/\ts2:i:\d+/.test(line)) {
-				++n_2nd;
-				continue;
-			}
-			if ((m = /\tcg:Z:(\S+)/.exec(line)) != null)
-				cigar = m[1];
-			if (cigar == null) {
-				warn("WARNING: no CIGAR at line " + lineno);
-				continue;
-			}
-			tname = t[5];
-			qs = parseInt(t[2]), qe = parseInt(t[3]);
-			aqlen = qe - qs;
-			is_rev = t[4] == '+'? false : true;
-			rs = parseInt(t[7]);
-			atlen = parseInt(t[8]) - rs;
-			mapq = parseInt(t[11]);
-			ori_qlen = parseInt(t[1]);
-		} else { // SAM
-			var flag = parseInt(t[1]);
-			if ((flag & 4) || t[2] == '*' || t[5] == '*') continue;
-			if (flag & 0x100) {
-				++n_2nd;
-				continue;
-			}
-			cigar = t[5];
-			tname = t[2];
-			rs = parseInt(t[3]) - 1;
-			mapq = parseInt(t[4]);
-			aqlen = t[9].length;
-			is_sam = true;
-			is_rev = !!(flag&0x10);
-		}
-		++n_pri;
-		if (last != t[0]) {
-			if (last != null) {
-				l_tot += last_qlen;
-				l_cov += cov_len(regs);
-			}
-			regs = [];
-			++n_seq, last = t[0];
-		}
-		var M = 0, tl = 0, ql = 0, clip = [0, 0], n_cigar = 0, sclip = 0;
-		while ((m = re.exec(cigar)) != null) {
-			var l = parseInt(m[1]);
-			++n_cigar;
-			if (m[2] == 'M' || m[2] == '=' || m[2] == 'X') {
-				tl += l, ql += l, M += l;
-			} else if (m[2] == 'I' || m[2] == 'D') {
-				var type;
-				if (l < 50) type = 0;
-				else if (l < 100) type = 1;
-				else if (l < 300) type = 2;
-				else if (l < 400) type = 3;
-				else if (l < 1000) type = 4;
-				else type = 5;
-				if (m[2] == 'I') ql += l, ++n_gap[0][type];
-				else tl += l, ++n_gap[1][type];
-				if (gap_out_len != null && l >= gap_out_len)
-					print(t[0], ql, is_rev? '-' : '+', tname, rs + tl, m[2], l);
-			} else if (m[2] == 'N') {
-				tl += l;
-			} else if (m[2] == 'S') {
-				clip[M == 0? 0 : 1] = l, sclip += l;
-			} else if (m[2] == 'H') {
-				clip[M == 0? 0 : 1] = l;
-			}
-		}
-		if (n_cigar > 65535) ++n_cigar_64k;
-		if (ql + sclip != aqlen)
-			warn("WARNING: aligned query length is inconsistent with CIGAR at line " + lineno + " (" + (ql+sclip) + " != " + aqlen + ")");
-		if (atlen != null && atlen != tl)
-			warn("WARNING: aligned reference length is inconsistent with CIGAR at line " + lineno);
-		if (is_sam) {
-			qs = clip[is_rev? 1 : 0], qe = qs + ql;
-			ori_qlen = clip[0] + ql + clip[1];
-		}
-		regs.push([qs, qe]);
-		last_qlen = ori_qlen;
-	}
-}
-l_tot += last_qlen;
-l_cov += cov_len(regs);
-
-file.close();
-buf.destroy();
-
-if (gap_out_len == null) {
-	print("Number of mapped sequences: " + n_seq);
-	print("Number of primary alignments: " + n_pri);
-	print("Number of secondary alignments: " + n_2nd);
-	print("Number of primary alignments with >65535 CIGAR operations: " + n_cigar_64k);
-	print("Number of bases in mapped sequences: " + l_tot);
-	print("Number of mapped bases: " + l_cov);
-	print("Number of insertions in [0,50): " + n_gap[0][0]);
-	print("Number of insertions in [50,100): " + n_gap[0][1]);
-	print("Number of insertions in [100,300): " + n_gap[0][2]);
-	print("Number of insertions in [300,400): " + n_gap[0][3]);
-	print("Number of insertions in [400,1000): " + n_gap[0][4]);
-	print("Number of insertions in [1000,inf): " + n_gap[0][5]);
-	print("Number of deletions in [0,50): " + n_gap[1][0]);
-	print("Number of deletions in [50,100): " + n_gap[1][1]);
-	print("Number of deletions in [100,300): " + n_gap[1][2]);
-	print("Number of deletions in [300,400): " + n_gap[1][3]);
-	print("Number of deletions in [400,1000): " + n_gap[1][4]);
-	print("Number of deletions in [1000,inf): " + n_gap[1][5]);
-}
@@ -1,105 +0,0 @@
-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");
@@ -1,196 +0,0 @@
-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, line_len = 80, fmt = "aln";
-while ((c = getopt(arguments, "f:l:")) != null) {
-	if (c == 'f') {
-		fmt = getopt.arg;
-		if (fmt != "aln" && fmt != "lastz-cigar" && fmt != "maf")
-			throw Error("format must be one of aln, lastz-cigar and maf");
-	} else if (c == 'l') line_len = parseInt(getopt.arg);
-}
-if (line_len == 0) line_len = 0x7fffffff;
-
-if (getopt.ind == arguments.length) {
-	print("Usage: k8 paf2aln.js [options] <in.paf>");
-	print("Options:");
-	print("  -f STR    output format: aln (BLAST-like), maf or lastz-cigar [aln]");
-	print("  -l INT    line length in BLAST-like output [80]");
-	exit(1);
-}
-
-function padding_str(x, len, right)
-{
-	var s = x.toString();
-	if (s.length < len) {
-		if (right) s += Array(len - s.length + 1).join(" ");
-		else s = Array(len - s.length + 1).join(" ") + s;
-	}
-	return s;
-}
-
-function update_aln(s_ref, s_qry, s_mid, type, seq, slen)
-{
-	var l = type == '*'? 1 : seq.length;
-	if (type == '=' || type == ':') {
-		s_ref.set(seq);
-		s_qry.set(seq);
-		s_mid.set(Array(l+1).join("|"));
-		slen[0] += l, slen[1] += l;
-	} else if (type == '*') {
-		s_ref.set(seq.charAt(0));
-		s_qry.set(seq.charAt(1));
-		s_mid.set(' ');
-		slen[0] += 1, slen[1] += 1;
-	} else if (type == '+') {
-		s_ref.set(Array(l+1).join("-"));
-		s_qry.set(seq);
-		s_mid.set(Array(l+1).join(" "));
-		slen[1] += l;
-	} else if (type == '-') {
-		s_ref.set(seq);
-		s_qry.set(Array(l+1).join("-"));
-		s_mid.set(Array(l+1).join(" "));
-		slen[0] += l;
-	}
-}
-
-function print_aln(rs, qs, strand, slen, elen, s_ref, s_qry, s_mid)
-{
-	print(["Ref+:", padding_str(rs + slen[0] + 1, 10, false), s_ref.toString(), padding_str(rs + elen[0], 10, true)].join(" "));
-	print("                 " + s_mid.toString());
-	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(" "));
-}
-
-var s_ref = new Bytes(), s_qry = new Bytes(), s_mid = new Bytes(); // these are used to show padded alignment
-var re_cs = /([:=\-\+\*])(\d+|[A-Za-z]+)/g;
-var re_cg = /(\d+)([MIDNSH])/g;
-
-var buf = new Bytes();
-var file = arguments[getopt.ind] == "-"? new File() : new File(arguments[getopt.ind]);
-var lineno = 0;
-if (fmt == "maf") print("##maf version=1\n");
-while (file.readline(buf) >= 0) {
-	var m, line = buf.toString();
-	var t = line.split("\t", 12);
-	++lineno;
-	s_ref.length = s_qry.length = s_mid.length = 0;
-	var slen = [0, 0], elen = [0, 0];
-	if (fmt == "lastz-cigar") { // LASTZ-cigar output
-		var cg = (m = /\tcg:Z:(\S+)/.exec(line)) != null? m[1] : null;
-		if (cg == null) {
-			warn("WARNING: converting to LASTZ-cigar format requires the 'cg' tag, which is absent on line " + lineno);
-			continue;
-		}
-		var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? m[1] : 0;
-		var out = ['cigar:', t[0], t[2], t[3], t[4], t[5], t[7], t[8], '+', score];
-		while ((m = re_cg.exec(cg)) != null)
-			out.push(m[2], m[1]);
-		print(out.join(" "));
-	} else if (fmt == "maf") { // MAF output
-		var cs = (m = /\tcs:Z:(\S+)/.exec(line)) != null? m[1] : null;
-		if (cs == null) {
-			warn("WARNING: converting to MAF requires the 'cs' tag, which is absent on line " + lineno);
-			continue;
-		}
-		while ((m = re_cs.exec(cs)) != null) {
-			if (m[1] == ':')
-				throw Error("converting to MAF only works with 'minimap2 --cs=long'");
-			update_aln(s_ref, s_qry, s_mid, m[1], m[2], elen);
-		}
-		var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : 0;
-		var len = t[0].length > t[5].length? t[0].length : t[5].length;
-		print("a " + score);
-		print(["s", padding_str(t[5], len, true), padding_str(t[7], 10, false), padding_str(parseInt(t[8]) - parseInt(t[7]), 10, false),
-				"+", padding_str(t[6], 10, false), s_ref.toString()].join(" "));
-		var qs, qe, ql = parseInt(t[1]);
-		if (t[4] == '+') {
-			qs = parseInt(t[2]);
-			qe = parseInt(t[3]);
-		} else {
-			qs = ql - parseInt(t[3]);
-			qe = ql - parseInt(t[2]);
-		}
-		print(["s", padding_str(t[0], len, true), padding_str(qs, 10, false), padding_str(qe - qs, 10, false),
-				t[4], padding_str(ql, 10, false), s_qry.toString()].join(" "));
-		print("");
-	} else { // BLAST-like output
-		var cs = (m = /\tcs:Z:(\S+)/.exec(line)) != null? m[1] : null;
-		if (cs == null) {
-			warn("WARNING: converting to BLAST-like alignment requires the 'cs' tag, which is absent on line " + lineno);
-			continue;
-		}
-		line = line.replace(/\tc[sg]:Z:\S+/g, ""); // get rid of cs or cg tags
-		print('>' + line);
-		var rs = parseInt(t[7]), qs = t[4] == '+'? parseInt(t[2]) : parseInt(t[3]);
-		var n_blocks = 0;
-		while ((m = re_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;
-				if (s_ref.length + rest >= line_len) {
-					l_proc = line_len - s_ref.length;
-					update_aln(s_ref, s_qry, s_mid, m[1], m[1] == '*'? m[2] : m[2].substr(start, l_proc), elen);
-					if (n_blocks > 0) print("");
-					print_aln(rs, qs, t[4], slen, elen, s_ref, s_qry, s_mid);
-					++n_blocks;
-					s_ref.length = s_qry.length = s_mid.length = 0;
-					slen[0] = elen[0], slen[1] = elen[1];
-				} else {
-					l_proc = rest;
-					update_aln(s_ref, s_qry, s_mid, m[1], m[1] == '*'? m[2] : m[2].substr(start, l_proc), elen);
-				}
-				rest -= l_proc, start += l_proc;
-			}
-		}
-		if (s_ref.length > 0) {
-			if (n_blocks > 0) print("");
-			print_aln(rs, qs, t[4], slen, elen, s_ref, s_qry, s_mid);
-			++n_blocks;
-		}
-		print("//");
-	}
-}
-file.close();
-buf.destroy();
-
-s_ref.destroy(); s_qry.destroy(); s_mid.destroy();
@@ -1,188 +0,0 @@
-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();
@@ -1,114 +0,0 @@
-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, pri_only = false;
-while ((c = getopt(arguments, "p")) != null)
-	if (c == 'p') pri_only = true;
-
-var file = arguments.length == getopt.ind? new File() : new File(arguments[getopt.ind]);
-var buf = new Bytes();
-var re = /(\d+)([MIDSHNX=])/g;
-
-var len = {}, lineno = 0;
-while (file.readline(buf) >= 0) {
-	var m, n_cigar = 0, line = buf.toString();
-	++lineno;
-	if (line.charAt(0) == '@') {
-		if (/^@SQ/.test(line)) {
-			var name = (m = /\tSN:(\S+)/.exec(line)) != null? m[1] : null;
-			var l = (m = /\tLN:(\d+)/.exec(line)) != null? parseInt(m[1]) : null;
-			if (name != null && l != null) len[name] = l;
-		}
-		continue;
-	}
-	var t = line.split("\t");
-	var flag = parseInt(t[1]);
-	if (t[9] != '*' && t[10] != '*' && t[9].length != t[10].length) throw Error("ERROR at line " + lineno + ": inconsistent SEQ and QUAL lengths - " + t[9].length + " != " + t[10].length);
-	if (t[2] == '*' || (flag&4)) continue;
-	if (pri_only && (flag&0x100)) continue;
-	var tlen = len[t[2]];
-	if (tlen == null) throw Error("ERROR at line " + lineno + ": can't find the length of contig " + t[2]);
-	var nn = (m = /\tnn:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : 0;
-	var NM = (m = /\tNM:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : null;
-	var have_NM = NM == null? false : true;
-	NM += nn;
-	var clip = [0, 0], I = [0, 0], D = [0, 0], M = 0, N = 0, ql = 0, tl = 0, mm = 0, ext_cigar = false;
-	while ((m = re.exec(t[5])) != null) {
-		var l = parseInt(m[1]);
-		if (m[2] == 'M') M += l, ql += l, tl += l, ext_cigar = false;
-		else if (m[2] == 'I') ++I[0], I[1] += l, ql += l;
-		else if (m[2] == 'D') ++D[0], D[1] += l, tl += l;
-		else if (m[2] == 'N') N += l, tl += l;
-		else if (m[2] == 'S') clip[M == 0? 0 : 1] = l, ql += l;
-		else if (m[2] == 'H') clip[M == 0? 0 : 1] = l;
-		else if (m[2] == '=') M += l, ql += l, tl += l, ext_cigar = true;
-		else if (m[2] == 'X') M += l, ql += l, tl += l, mm += l, ext_cigar = true;
-		++n_cigar;
-	}
-	if (n_cigar > 65535)
-		warn("WARNING at line " + lineno + ": " + n_cigar + " CIGAR operations");
-	if (tl + parseInt(t[3]) - 1 > tlen) {
-		warn("WARNING at line " + lineno + ": alignment end position larger than ref length; skipped");
-		continue;
-	}
-	if (t[9] != '*' && t[9].length != ql) {
-		warn("WARNING at line " + lineno + ": SEQ length inconsistent with CIGAR (" + t[9].length + " != " + ql + "); skipped");
-		continue;
-	}
-	if (!have_NM || ext_cigar) NM = I[1] + D[1] + mm;
-	if (NM < I[1] + D[1] + mm) {
-		warn("WARNING at line " + lineno + ": NM is less than the total number of gaps (" + NM + " < " + (I[1]+D[1]+mm) + ")");
-		NM = I[1] + D[1] + mm;
-	}
-	var extra = ["mm:i:"+(NM-I[1]-D[1]), "io:i:"+I[0], "in:i:"+I[1], "do:i:"+D[0], "dn:i:"+D[1]];
-	var match = M - (NM - I[1] - D[1]);
-	var blen = M + I[1] + D[1];
-	var qlen = M + I[1] + clip[0] + clip[1];
-	var qs, qe;
-	if (flag&16) qs = clip[1], qe = qlen - clip[0];
-	else qs = clip[0], qe = qlen - clip[1];
-	var ts = parseInt(t[3]) - 1, te = ts + M + D[1] + N;
-	var qname = t[0];
-	if ((flag&1) && (flag&0x40)) qname += '/1';
-	if ((flag&1) && (flag&0x80)) qname += '/2';
-	var a = [qname, qlen, qs, qe, flag&16? '-' : '+', t[2], tlen, ts, te, match, blen, t[4]];
-	print(a.join("\t"), extra.join("\t"));
-}
-
-buf.destroy();
-file.close();
@@ -1,193 +0,0 @@
-var getopt = function(args, ostr) {
-	var oli; // option letter list index
-	if (typeof(getopt.place) == 'undefined')
-		getopt.ind = 0, getopt.arg = null, getopt.place = -1;
-	if (getopt.place == -1) { // update scanning pointer
-		if (getopt.ind >= args.length || args[getopt.ind].charAt(getopt.place = 0) != '-') {
-			getopt.place = -1;
-			return null;
-		}
-		if (getopt.place + 1 < args[getopt.ind].length && args[getopt.ind].charAt(++getopt.place) == '-') { // found "--"
-			++getopt.ind;
-			getopt.place = -1;
-			return null;
-		}
-	}
-	var optopt = args[getopt.ind].charAt(getopt.place++); // character checked for validity
-	if (optopt == ':' || (oli = ostr.indexOf(optopt)) < 0) {
-		if (optopt == '-') return null; //  if the user didn't specify '-' as an option, assume it means null.
-		if (getopt.place < 0) ++getopt.ind;
-		return '?';
-	}
-	if (oli+1 >= ostr.length || ostr.charAt(++oli) != ':') { // don't need argument
-		getopt.arg = null;
-		if (getopt.place < 0 || getopt.place >= args[getopt.ind].length) ++getopt.ind, getopt.place = -1;
-	} else { // need an argument
-		if (getopt.place >= 0 && getopt.place < args[getopt.ind].length)
-			getopt.arg = args[getopt.ind].substr(getopt.place);
-		else if (args.length <= ++getopt.ind) { // no arg
-			getopt.place = -1;
-			if (ostr.length > 0 && ostr.charAt(0) == ':') return ':';
-			return '?';
-		} else getopt.arg = args[getopt.ind]; // white space
-		getopt.place = -1;
-		++getopt.ind;
-	}
-	return optopt;
-}
-
-var c, max_mapq = 60, mode = 0, err_out_q = 256, print_err = false, ovlp_ratio = 0.1, cap_short_mapq = false;
-while ((c = getopt(arguments, "Q:r:m:c")) != null) {
-	if (c == 'Q') err_out_q = parseInt(getopt.arg), print_err = true;
-	else if (c == 'r') ovlp_ratio = parseFloat(getopt.arg);
-	else if (c == 'm') mode = parseInt(getopt.arg);
-	else if (c == 'c') cap_short_mapq = true;
-}
-
-var file = arguments.length == getopt.ind? new File() : new File(arguments[getopt.ind]);
-var buf = new Bytes();
-
-var tot = [], err = [];
-for (var q = 0; q <= max_mapq; ++q)
-	tot[q] = err[q] = 0;
-
-function is_correct(s, b)
-{
-	if (s[0] != b[0] || s[3] != b[3]) return false;
-	var o, l;
-	if (s[1] < b[1]) {
-		if (s[2] <= b[1]) return false;
-		o = (s[2] < b[2]? s[2] : b[2]) - b[1];
-		l = (s[2] > b[2]? s[2] : b[2]) - s[1];
-	} else {
-		if (b[2] <= s[1]) return false;
-		o = (s[2] < b[2]? s[2] : b[2]) - s[1];
-		l = (s[2] > b[2]? s[2] : b[2]) - b[1];
-	}
-	return o/l > ovlp_ratio? true : false;
-}
-
-function count_err(qname, a, tot, err, mode)
-{
-	if (a.length == 0) return;
-
-	var m, s;
-	if ((m = /^(\S+)!(\S+)!(\d+)!(\d+)!([\+\-])$/.exec(qname)) != null) { // pbsim single-end reads
-		s = [m[1], m[2], parseInt(m[3]), parseInt(m[4]), m[5]];
-	} else if ((m = /^(\S+)!(\S+)!(\d+)_(\d+)!(\d+)_(\d+)!([\+\-])([\+\-])\/([12])$/.exec(qname)) != null) { // mason2 paired-end reads
-		if (m[9] == '1') {
-			s = [m[1], m[2], parseInt(m[3]), parseInt(m[5]), m[7]];
-		} else {
-			s = [m[1], m[2], parseInt(m[4]), parseInt(m[6]), m[8]];
-		}
-	} else throw Error("Failed to parse simulated read names '" + qname + "'");
-	s.shift(); // skip the orginal read name
-
-	if (mode == 0 || mode == 1) { // longest only or first only
-		var max_i = 0;
-		if (mode == 0) { // longest only
-			var max = 0;
-			for (var i = 0; i < a.length; ++i)
-				if (a[i][5] > max)
-					max = a[i][5], max_i = i;
-		}
-		var mapq = a[max_i][4];
-		++tot[mapq];
-		if (!is_correct(s, a[max_i])) {
-			if (mapq >= err_out_q)
-				print('E', qname, a[max_i].join("\t"));
-			++err[mapq];
-		}
-	} else if (mode == 2) { // all primary mode
-		var max_err_mapq = -1, max_mapq = 0, max_err_i = -1;
-		if (cap_short_mapq) {
-			var max = 0, max_q = 0;
-			for (var i = 0; i < a.length; ++i)
-				if (a[i][5] > max)
-					max = a[i][5], max_q = a[i][4];
-			for (var i = 0; i < a.length; ++i)
-				a[i][4] = max_q < a[i][4]? max_q : a[i][4];
-		}
-		for (var i = 0; i < a.length; ++i) {
-			max_mapq = max_mapq > a[i][4]? max_mapq : a[i][4];
-			if (!is_correct(s, a[i]))
-				if (a[i][4] > max_err_mapq)
-					max_err_mapq = a[i][4], max_err_i = i;
-		}
-		if (max_err_mapq >= 0) {
-			++tot[max_err_mapq], ++err[max_err_mapq];
-			if (max_err_mapq >= err_out_q)
-				print('E', qname, a[max_err_i].join("\t"));
-		} else ++tot[max_mapq];
-	}
-}
-
-var lineno = 0, last = null, a = [], n_unmapped = null;
-var re_cigar = /(\d+)([MIDSHN])/g;
-while (file.readline(buf) >= 0) {
-	var m, line = buf.toString();
-	++lineno;
-	if (line[0] != '@') {
-		var t = line.split("\t");
-		if (t[4] == '+' || t[4] == '-') { // PAF
-			if (last != t[0]) {
-				if (last != null) count_err(last, a, tot, err, mode);
-				a = [], last = t[0];
-			}
-			if (/\ts1:i:\d+/.test(line) && !/\ts2:i:\d+/.test(line)) // secondary alignment in minimap2 PAF
-				continue;
-			var mapq = parseInt(t[11]);
-			if (mapq > max_mapq) mapq = max_mapq;
-			a.push([t[5], parseInt(t[7]), parseInt(t[8]), t[4], mapq, parseInt(t[9])]);
-		} else { // SAM
-			var flag = parseInt(t[1]);
-			var read_no = flag>>6&0x3;
-			var qname = t[0];
-			if (!/\/[12]$/.test(qname))
-				qname = read_no == 1 || read_no == 2? t[0] + '/' + read_no : t[0];
-			if (last != qname) {
-				if (last != null) count_err(last, a, tot, err, mode);
-				a = [], last = qname;
-			}
-			if (flag&0x100) continue; // secondary alignment
-			if ((flag&0x4) || t[2] == '*') { // unmapped
-				if (n_unmapped == null) n_unmapped = 0;
-				++n_unmapped;
-				continue;
-			}
-			var mapq = parseInt(t[4]);
-			if (mapq > max_mapq) mapq = max_mapq;
-			var pos = parseInt(t[3]) - 1, pos_end = pos;
-			var n_gap = 0, mlen = 0;
-			while ((m = re_cigar.exec(t[5])) != null) {
-				var len = parseInt(m[1]);
-				if (m[2] == 'M') pos_end += len, mlen += len;
-				else if (m[2] == 'I') n_gap += len;
-				else if (m[2] == 'D') n_gap += len, pos_end += len;
-			}
-			var score = pos_end - pos;
-			if ((m = /\tNM:i:(\d+)/.exec(line)) != null) {
-				var NM = parseInt(m[1]);
-				if (NM >= n_gap) score = mlen - (NM - n_gap);
-			}
-			a.push([t[2], pos, pos_end, (flag&16)? '-' : '+', mapq, score]);
-		}
-	}
-}
-if (last != null) count_err(last, a, tot, err, mode);
-
-buf.destroy();
-file.close();
-
-var sum_tot = 0, sum_err = 0, q_out = -1, sum_tot2 = 0, sum_err2 = 0;
-for (var q = max_mapq; q >= 0; --q) {
-	if (tot[q] == 0) continue;
-	if (q_out < 0 || err[q] > 0) {
-		if (q_out >= 0) print('Q', q_out, sum_tot, sum_err, (sum_err2/sum_tot2).toFixed(9), sum_tot2);
-		sum_tot = sum_err = 0, q_out = q;
-	}
-	sum_tot += tot[q], sum_err += err[q];
-	sum_tot2 += tot[q], sum_err2 += err[q];
-}
-print('Q', q_out, sum_tot, sum_err, (sum_err2/sum_tot2).toFixed(9), sum_tot2);
-if (n_unmapped != null) print('U', n_unmapped);
@@ -1,105 +0,0 @@
-Bytes.prototype.reverse = function()
-{
-	for (var i = 0; i < this.length>>1; ++i) {
-		var tmp = this[i];
-		this[i] = this[this.length - i - 1];
-		this[this.length - i - 1] = tmp;
-	}
-}
-
-// reverse complement a DNA string
-Bytes.prototype.revcomp = function()
-{
-	if (Bytes.rctab == null) {
-		var s1 = 'WSATUGCYRKMBDHVNwsatugcyrkmbdhvn';
-		var s2 = 'WSTAACGRYMKVHDBNwstaacgrymkvhdbn';
-		Bytes.rctab = [];
-		for (var i = 0; i < 256; ++i) Bytes.rctab[i] = 0;
-		for (var i = 0; i < s1.length; ++i)
-			Bytes.rctab[s1.charCodeAt(i)] = s2.charCodeAt(i);
-	}
-	for (var i = 0; i < this.length>>1; ++i) {
-		var tmp = this[this.length - i - 1];
-		this[this.length - i - 1] = Bytes.rctab[this[i]];
-		this[i] = Bytes.rctab[tmp];
-	}
-	if (this.length&1)
-		this[this.length>>1] = Bytes.rctab[this[this.length>>1]];
-}
-
-if (arguments.length == 0) {
-	print("Usage: k8 sim-mason2.js <mason.sam>");
-	exit(1);
-}
-
-function print_se(a)
-{
-	print('@' + a.slice(0, 5).join("!") + " " + a[8]);
-	print(a[5]);
-	print("+");
-	print(a[6]);
-}
-
-var buf = new Bytes(), buf2 = new Bytes();
-var file = new File(arguments[0]);
-var re = /(\d+)([MIDSHN])/g;
-var last = null;
-while (file.readline(buf) >= 0) {
-	var t = buf.toString().split("\t");
-	if (t[0].charAt(0) == '@') continue;
-	var m, l_ref = 0;
-	while ((m = re.exec(t[5])) != null)
-		if (m[2] == 'D' || m[2] == 'M' || m[2] == 'N')
-			l_ref += parseInt(m[1]);
-	var flag = parseInt(t[1]);
-	var rev = !!(flag&16);
-	var seq, qual;
-	if (rev) {
-		buf2.length = 0;
-		buf2.set(t[9], 0);
-		buf2.revcomp();
-		seq = buf2.toString();
-		buf2.set(t[10], 0);
-		buf2.reverse();
-		qual = buf2.toString();
-	} else seq = t[9], qual = t[10];
-	var qname = t[0];
-	qname = qname.replace(/^simulated./, "");
-	var chr = t[2];
-	var pos = parseInt(t[3]) - 1;
-	var strand = (flag&16)? '-' : '+';
-	var read_no = flag&0xc0;
-	if (read_no == 0x40) read_no = 1;
-	else if (read_no == 0x80) read_no = 2;
-	else read_no = 0;
-	var err = 0, snp = 0, indel = 0;
-	for (var i = 11; i < t.length; ++i) {
-		if ((m = /^XE:i:(\d+)/.exec(t[i])) != null) err = m[1];
-		else if ((m = /^XS:i:(\d+)/.exec(t[i])) != null) snp = m[1];
-		else if ((m = /^XI:i:(\d+)/.exec(t[i])) != null) indel = m[1];
-	}
-	var comment = [err, snp, indel].join(":");
-	if (last == null) {
-		last = [qname, chr, pos, pos + l_ref, strand, seq, qual, read_no, comment];
-	} else if (last[0] != qname) {
-		print_se(last);
-		last = [qname, chr, pos, pos + l_ref, strand, seq, qual, read_no, comment];
-	} else {
-		if (read_no == 2) { // last[] is the first read
-			if (last[7] != 1) throw Error("ERROR: can't find read1");
-			var name = [qname, chr, last[2] + "_" + pos, last[3] + "_" + (pos + l_ref), last[4] + strand].join("!");
-			print('@' + name + '/1' + ' ' + last[8]); print(last[5]); print("+"); print(last[6]);
-			print('@' + name + '/2' + ' ' + comment); print(seq); print("+"); print(qual);
-		} else {
-			if (last[7] != 2) throw Error("ERROR: can't find read2");
-			var name = [qname, chr, pos + "_" + last[2], (pos + l_ref) + "_" + last[3], strand + last[4]].join("!");
-			print('@' + name + '/1' + ' ' + comment); print(seq); print("+"); print(qual);
-			print('@' + name + '/2' + ' ' + last[8]); print(last[5]); print("+"); print(last[6]);
-		}
-		last = null;
-	}
-}
-if (last != null) print_se(last);
-file.close();
-buf.destroy();
-buf2.destroy();
@@ -1,81 +0,0 @@
-Bytes.prototype.reverse = function()
-{
-	for (var i = 0; i < this.length>>1; ++i) {
-		var tmp = this[i];
-		this[i] = this[this.length - i - 1];
-		this[this.length - i - 1] = tmp;
-	}
-}
-
-// reverse complement a DNA string
-Bytes.prototype.revcomp = function()
-{
-	if (Bytes.rctab == null) {
-		var s1 = 'WSATUGCYRKMBDHVNwsatugcyrkmbdhvn';
-		var s2 = 'WSTAACGRYMKVHDBNwstaacgrymkvhdbn';
-		Bytes.rctab = [];
-		for (var i = 0; i < 256; ++i) Bytes.rctab[i] = 0;
-		for (var i = 0; i < s1.length; ++i)
-			Bytes.rctab[s1.charCodeAt(i)] = s2.charCodeAt(i);
-	}
-	for (var i = 0; i < this.length>>1; ++i) {
-		var tmp = this[this.length - i - 1];
-		this[this.length - i - 1] = Bytes.rctab[this[i]];
-		this[i] = Bytes.rctab[tmp];
-	}
-	if (this.length&1)
-		this[this.length>>1] = Bytes.rctab[this[this.length>>1]];
-}
-
-if (arguments.length < 2) {
-	print("Usage: k8 sim-pbsim.js <ref.fa.fai> <pbsim1.maf> [[pbsim2.maf] ...]");
-	exit(1);
-}
-
-var file, buf = new Bytes(), buf2 = new Bytes();
-file = new File(arguments[0]);
-var chr_list = [];
-while (file.readline(buf) >= 0) {
-	var t = buf.toString().split(/\s+/);
-	chr_list.push(t[0]);
-}
-file.close();
-
-for (var k = 1; k < arguments.length; ++k) {
-	var fn = arguments[k];
-	file = new File(fn);
-	var state = 0, reg;
-	while (file.readline(buf) >= 0) {
-		var line = buf.toString();
-		if (state == 0 && line.charAt(0) == 'a') {
-			state = 1;
-		} else if (state == 1 && line.charAt(0) == 's') {
-			var t = line.split(/\s+/);
-			var st = parseInt(t[2]);
-			reg = [st, st + parseInt(t[3])];
-			state = 2;
-		} else if (state == 2 && line.charAt(0) == 's') {
-			var m, t = line.split(/\s+/);
-			if ((m = /S(\d+)_\d+/.exec(t[1])) == null) throw Error("Failed to parse the read name");
-			var chr_id = parseInt(m[1]) - 1;
-			if (chr_id >= chr_list.length) throw Error("Index outside the chr list");
-			var name = [t[1], chr_list[chr_id], reg[0], reg[1], t[4]].join("!");
-			var seq = t[6].replace(/\-/g, "");
-			if (seq.length != parseInt(t[5])) throw Error("Inconsistent read length");
-			if (seq.indexOf("NN") < 0) {
-				if (t[4] == '-') {
-					buf2.set(seq, 0);
-					buf2.length = seq.length;
-					buf2.revcomp();
-					seq = buf2.toString();
-				}
-				print(">" + name);
-				print(seq);
-			}
-			state = 0;
-		}
-	}
-	file.close();
-}
-buf.destroy();
-buf2.destroy();
@@ -1,148 +0,0 @@
-var getopt = function(args, ostr) {
-	var oli; // option letter list index
-	if (typeof(getopt.place) == 'undefined')
-		getopt.ind = 0, getopt.arg = null, getopt.place = -1;
-	if (getopt.place == -1) { // update scanning pointer
-		if (getopt.ind >= args.length || args[getopt.ind].charAt(getopt.place = 0) != '-') {
-			getopt.place = -1;
-			return null;
-		}
-		if (getopt.place + 1 < args[getopt.ind].length && args[getopt.ind].charAt(++getopt.place) == '-') { // found "--"
-			++getopt.ind;
-			getopt.place = -1;
-			return null;
-		}
-	}
-	var optopt = args[getopt.ind].charAt(getopt.place++); // character checked for validity
-	if (optopt == ':' || (oli = ostr.indexOf(optopt)) < 0) {
-		if (optopt == '-') return null; //  if the user didn't specify '-' as an option, assume it means null.
-		if (getopt.place < 0) ++getopt.ind;
-		return '?';
-	}
-	if (oli+1 >= ostr.length || ostr.charAt(++oli) != ':') { // don't need argument
-		getopt.arg = null;
-		if (getopt.place < 0 || getopt.place >= args[getopt.ind].length) ++getopt.ind, getopt.place = -1;
-	} else { // need an argument
-		if (getopt.place >= 0 && getopt.place < args[getopt.ind].length)
-			getopt.arg = args[getopt.ind].substr(getopt.place);
-		else if (args.length <= ++getopt.ind) { // no arg
-			getopt.place = -1;
-			if (ostr.length > 0 && ostr.charAt(0) == ':') return ':';
-			return '?';
-		} else getopt.arg = args[getopt.ind]; // white space
-		getopt.place = -1;
-		++getopt.ind;
-	}
-	return optopt;
-}
-
-var colors = ["0,128,255", "255,0,0", "0,192,0"];
-
-function print_lines(a, fmt) {
-	if (a.length == 0) return;
-	if (fmt == "bed") {
-		var n_pri = 0;
-		for (var i = 0; i < a.length; ++i)
-			if (a[i][8] == 0) ++n_pri;
-		if (n_pri > 1) {
-			for (var i = 0; i < a.length; ++i)
-				if (a[i][8] == 0) a[i][8] = 1;
-		} else if (n_pri == 0) {
-			warn("Warning: " + a[0][3] + " doesn't have a primary alignment");
-		}
-		for (var i = 0; i < a.length; ++i) {
-			a[i][8] = colors[a[i][8]];
-			print(a[i].join("\t"));
-		}
-	}
-	a.length = 0;
-}
-
-function main(args) {
-	var re = /(\d+)([MIDNSH])/g;
-	var c, fmt = "bed", fn_name_conv = null;
-	while ((c = getopt(args, "f:n:")) != null) {
-		if (c == 'f') fmt = getopt.arg;
-		else if (c == 'n') fn_name_conv = getopt.arg;
-	}
-	if (getopt.ind == args.length) {
-		warn("Usage: k8 splice2bed.js <in.paf>");
-		exit(1);
-	}
-
-	var conv = null;
-	if (fn_name_conv != null) {
-		conv = new Map();
-		var file = new File(fn_name_conv);
-		var buf = new Bytes();
-		while (file.readline(buf) >= 0) {
-			var t = buf.toString().split("\t");
-			conv.put(t[0], t[1]);
-		}
-		buf.destroy();
-		file.close();
-	}
-
-	var file = new File(args[getopt.ind]);
-	var buf = new Bytes();
-	var a = [];
-	while (file.readline(buf) >= 0) {
-		var line = buf.toString();
-		if (line.charAt(0) == '@') continue; // skip SAM header lines
-		var t = line.split("\t");
-		var is_pri = false, cigar = null, a1;
-		var qname = conv != null? conv.get(t[0]) : null;
-		if (qname != null) t[0] = qname;
-		if (t.length >= 10 && t[4] != '+' && t[4] != '-' && /^\d+/.test(t[1])) { // SAM
-			var flag = parseInt(t[1]);
-			if (flag&1) t[0] += '/' + (flag>>6&3);
-		}
-		if (a.length && a[0][3] != t[0]) {
-			print_lines(a, fmt);
-			a = [];
-		}
-		if (t.length >= 12 && (t[4] == '+' || t[4] == '-')) { // PAF
-			for (var i = 12; i < t.length; ++i) {
-				if (t[i].substr(0, 5) == 'cg:Z:') {
-					cigar = t[i].substr(5);
-				} else if (t[i].substr(0, 5) == 's2:i:') {
-					is_pri = true;
-				}
-			}
-			a1 = [t[5], t[7], t[8], t[0], Math.floor(t[9]/t[10]*1000), t[4]];
-		} else if (t.length >= 10) { // SAM
-			var flag = parseInt(t[1]);
-			if ((flag&4) || a[2] == '*') continue;
-			cigar = t[5];
-			is_pri = (flag&0x100)? false : true;
-			a1 = [t[2], parseInt(t[3])-1, null, t[0], 1000, (flag&16)? '-' : '+'];
-		} else {
-			throw Error("unrecognized input format");
-		}
-		if (cigar == null) throw Error("missing CIGAR");
-		var m, x0 = 0, x = 0, bs = [], bl = [];
-		while ((m = re.exec(cigar)) != null) {
-			if (m[2] == 'M' || m[2] == 'D') {
-				x += parseInt(m[1]);
-			} else if (m[2] == 'N') {
-				bs.push(x0);
-				bl.push(x - x0);
-				x += parseInt(m[1]);
-				x0 = x;
-			}
-		}
-		bs.push(x0);
-		bl.push(x - x0);
-		// write the BED12 line
-		if (a1[2] == null) a1[2] = a1[1] + x;
-		a1.push(a1[1], a1[2]); // thick start/end is the same as start/end
-		a1.push(is_pri? 0 : 2, bs.length, bl.join(",")+",", bs.join(",")+",");
-		a.push(a1);
-	}
-	print_lines(a, fmt);
-	buf.destroy();
-	file.close();
-	if (conv != null) conv.destroy();
-}
-
-main(arguments);
@@ -16,6 +16,7 @@
 #define MM_SEED_LONG_JOIN  (1ULL<<40)
 #define MM_SEED_IGNORE     (1ULL<<41)
 #define MM_SEED_TANDEM     (1ULL<<42)
+#define MM_SEED_SELF       (1ULL<<43)

 #define MM_SEED_SEG_SHIFT  48
 #define MM_SEED_SEG_MASK   (0xffULL<<(MM_SEED_SEG_SHIFT))
@@ -47,6 +48,7 @@ typedef struct {

 double cputime(void);
 double realtime(void);
+long peakrss(void);

 void radix_sort_128x(mm128_t *beg, mm128_t *end);
 void radix_sort_64(uint64_t *beg, uint64_t *end);
@@ -61,7 +63,6 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se

 void mm_idxopt_init(mm_idxopt_t *opt);
 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, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
@@ -74,10 +75,10 @@ int mm_set_sam_pri(int n, mm_reg1_t *r);
 void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff);
 void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_, mm_reg1_t *r);
 void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int max_gap_ref, int min_diff, int best_n, int n_segs, const int *qlens, int *n_, mm_reg1_t *r);
-void mm_filter_regs(void *km, const mm_mapopt_t *opt, int *n_regs, mm_reg1_t *regs);
+void mm_filter_regs(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs);
 void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs, mm128_t *a);
 void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r);
-void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr);
+void mm_set_mapq(void *km, 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);

@@ -0,0 +1,155 @@
+#include <stdio.h>
+#include "mmpriv.h"
+
+void mm_idxopt_init(mm_idxopt_t *opt)
+{
+	memset(opt, 0, sizeof(mm_idxopt_t));
+	opt->k = 15, opt->w = 10, opt->flag = 0;
+	opt->bucket_bits = 14;
+	opt->mini_batch_size = 50000000;
+	opt->batch_size = 4000000000ULL;
+}
+
+void mm_mapopt_init(mm_mapopt_t *opt)
+{
+	memset(opt, 0, sizeof(mm_mapopt_t));
+	opt->seed = 11;
+	opt->mid_occ_frac = 2e-4f;
+	opt->sdust_thres = 0; // no SDUST masking
+
+	opt->min_cnt = 3;
+	opt->min_chain_score = 40;
+	opt->bw = 500;
+	opt->max_gap = 5000;
+	opt->max_gap_ref = -1;
+	opt->max_chain_skip = 25;
+
+	opt->mask_level = 0.5f;
+	opt->pri_ratio = 0.8f;
+	opt->best_n = 5;
+
+	opt->max_join_long = 20000;
+	opt->max_join_short = 2000;
+	opt->min_join_flank_sc = 1000;
+
+	opt->a = 2, opt->b = 4, opt->q = 4, opt->e = 2, opt->q2 = 24, opt->e2 = 1;
+	opt->zdrop = 400, opt->zdrop_inv = 200;
+	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->max_clip_ratio = 1.0f;
+	opt->mini_batch_size = 500000000;
+
+	opt->pe_ori = 0; // FF
+	opt->pe_bonus = 33;
+}
+
+void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)
+{
+	if ((opt->flag & MM_F_SPLICE_FOR) && (opt->flag & MM_F_SPLICE_REV))
+		opt->flag |= MM_F_SPLICE;
+	if (opt->mid_occ <= 0)
+		opt->mid_occ = mm_idx_cal_max_occ(mi, opt->mid_occ_frac);
+	if (mm_verbose >= 3)
+		fprintf(stderr, "[M::%s::%.3f*%.2f] mid_occ = %d\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), opt->mid_occ);
+}
+
+void mm_mapopt_max_intron_len(mm_mapopt_t *opt, int max_intron_len)
+{
+	if ((opt->flag & MM_F_SPLICE) && max_intron_len > 0)
+		opt->max_gap_ref = opt->bw = max_intron_len;
+}
+
+int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
+{
+	if (preset == 0) {
+		mm_idxopt_init(io);
+		mm_mapopt_init(mo);
+	} else if (strcmp(preset, "ava-ont") == 0) {
+		io->flag = 0, io->k = 15, io->w = 5;
+		mo->flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN;
+		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
+	} else if (strcmp(preset, "ava-pb") == 0) {
+		io->flag |= MM_I_HPC, io->k = 19, io->w = 5;
+		mo->flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN;
+		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
+	} else if (strcmp(preset, "map10k") == 0 || strcmp(preset, "map-pb") == 0) {
+		io->flag |= MM_I_HPC, io->k = 19;
+	} else if (strcmp(preset, "map-ont") == 0) {
+		io->flag = 0, io->k = 15;
+	} else if (strcmp(preset, "asm5") == 0) {
+		io->flag = 0, io->k = 19, io->w = 19;
+		mo->a = 1, mo->b = 19, mo->q = 39, mo->q2 = 81, mo->e = 3, mo->e2 = 1, mo->zdrop = mo->zdrop_inv = 200;
+		mo->min_dp_max = 200;
+		mo->best_n = 50;
+	} else if (strcmp(preset, "asm10") == 0) {
+		io->flag = 0, io->k = 19, io->w = 19;
+		mo->a = 1, mo->b = 9, mo->q = 16, mo->q2 = 41, mo->e = 2, mo->e2 = 1, mo->zdrop = mo->zdrop_inv = 200;
+		mo->min_dp_max = 200;
+		mo->best_n = 50;
+	} else if (strcmp(preset, "short") == 0 || strcmp(preset, "sr") == 0) {
+		io->flag = 0, io->k = 21, io->w = 11;
+		mo->flag |= MM_F_SR | MM_F_FRAG_MODE | MM_F_NO_PRINT_2ND | MM_F_2_IO_THREADS | MM_F_HEAP_SORT;
+		mo->pe_ori = 0<<1|1; // FR
+		mo->a = 2, mo->b = 8, mo->q = 12, mo->e = 2, mo->q2 = 24, mo->e2 = 1;
+		mo->zdrop = mo->zdrop_inv = 100;
+		mo->end_bonus = 10;
+		mo->max_frag_len = 800;
+		mo->max_gap = 100;
+		mo->bw = 100;
+		mo->pri_ratio = 0.5f;
+		mo->min_cnt = 2;
+		mo->min_chain_score = 25;
+		mo->min_dp_max = 40;
+		mo->best_n = 20;
+		mo->mid_occ = 1000;
+		mo->max_occ = 5000;
+		mo->mini_batch_size = 50000000;
+	} else if (strcmp(preset, "splice") == 0 || strcmp(preset, "cdna") == 0) {
+		io->flag = 0, io->k = 15, io->w = 5;
+		mo->flag |= MM_F_SPLICE | MM_F_SPLICE_FOR | MM_F_SPLICE_REV | MM_F_SPLICE_FLANK;
+		mo->max_gap = 2000, mo->max_gap_ref = mo->bw = 200000;
+		mo->a = 1, mo->b = 2, mo->q = 2, mo->e = 1, mo->q2 = 32, mo->e2 = 0;
+		mo->noncan = 9;
+		mo->zdrop = 200, mo->zdrop_inv = 100; // because mo->a is halved
+	} else return -1;
+	return 0;
+}
+
+int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
+{
+	if (mo->best_n < 0) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m -N must be no less than 0\033[0m\n");
+		return -4;
+	}
+	if (mo->best_n == 0 && mm_verbose >= 2)
+		fprintf(stderr, "[WARNING]\033[1;31m '-N 0' reduces mapping accuracy. Please use '--secondary=no' instead.\033[0m\n");
+	if (mo->pri_ratio < 0.0f || mo->pri_ratio > 1.0f) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m -p must be within 0 and 1 (including 0 and 1)\033[0m\n");
+		return -4;
+	}
+	if ((mo->flag & MM_F_FOR_ONLY) && (mo->flag & MM_F_REV_ONLY)) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m --for-only and --rev-only can't be applied at the same time\033[0m\n");
+		return -3;
+	}
+	if ((mo->q != mo->q2 || mo->e != mo->e2) && !(mo->e > mo->e2 && mo->q + mo->e < mo->q2 + mo->e2)) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m dual gap penalties violating E1>E2 and O1+E1<O2+E2\033[0m\n");
+		return -2;
+	}
+	if ((mo->q + mo->e) + (mo->q2 + mo->e2) > 127) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m scoring system violating ({-O}+{-E})+({-O2}+{-E2}) <= 127\033[0m\n");
+		return -1;
+	}
+	if (mo->zdrop < mo->zdrop_inv) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m Z-drop should not be less than inversion-Z-drop\033[0m\n");
+		return -5;
+	}
+	return 0;
+}
@@ -34,6 +34,8 @@ The following Python script demonstrates the key functionality of mappy:
 	import mappy as mp
 	a = mp.Aligner("test/MT-human.fa")  # load or build index
 	if not a: raise Exception("ERROR: failed to load/build index")
+	s = a.seq("MT_human", 100, 200)     # retrieve a subsequence from the index
+	print(mp.revcomp(s))                # reverse complement
 	for name, seq, qual in mp.fastx_read("test/MT-orang.fa"): # read a fasta/q sequence
 		for hit in a.map(seq): # traverse alignments
 			print("{}\t{}\t{}\t{}".format(hit.ctg, hit.r_st, hit.r_en, hit.cigar_str))
@@ -81,10 +83,21 @@ This constructor accepts the following arguments:

 .. code:: python

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

 This method aligns :code:`seq` against the index. It is a generator, *yielding*
-a series of :code:`mappy.Alignment` objects.
+a series of :code:`mappy.Alignment` objects. If :code:`seq2` is present, mappy
+performs paired-end alignment, assuming the two ends are in the FR orientation.
+Alignments of the two ends can be distinguished by the :code:`read_num` field
+(see Class mappy.Alignment below).
+
+.. code:: python
+
+	mappy.Aligner.seq(name, start=0, end=0x7fffffff)
+
+This method retrieves a (sub)sequence from the index and returns it as a Python
+string. :code:`None` is returned if :code:`name` is not present in the index or
+the start/end coordinates are invalid.

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

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

 * **cigar**: CIGAR returned as an array of shape :code:`(n_cigar,2)`. The two
@@ -133,13 +149,22 @@ the following format:
 It is effectively the PAF format without the QueryName and QueryLength columns
 (the first two columns in PAF).

-Function mappy.fastx_read
-~~~~~~~~~~~~~~~~~~~~~~~~~
+Miscellaneous Functions
+~~~~~~~~~~~~~~~~~~~~~~~

 .. code:: python

-	mappy.fastx_read(fn)
+	mappy.fastx_read(fn, read_comment=False)

 This generator function opens a FASTA/FASTQ file and *yields* a
 :code:`(name,seq,qual)` tuple for each sequence entry. The input file may be
-optionally gzip'd.
+optionally gzip'd. If :code:`read_comment` is True, this generator yields
+a :code:`(name,seq,qual,comment)` tuple instead.
+
+.. code:: python
+
+	mappy.revcomp(seq)
+
+Return the reverse complement of DNA string :code:`seq`. This function
+recognizes IUB code and preserves the letter cases. Uracil :code:`U` is
+complemented to :code:`A`.
@@ -15,6 +15,7 @@ typedef struct {
 	int32_t blen, mlen, NM, ctg_len;
 	uint8_t mapq, is_primary;
 	int8_t strand, trans_strand;
+	int32_t seg_id;
 	int32_t n_cigar32;
 	uint32_t *cigar32;
 } mm_hitpy_t;
@@ -32,6 +33,7 @@ static inline void mm_reg2hitpy(const mm_idx_t *mi, mm_reg1_t *r, mm_hitpy_t *h)
 	h->NM = r->blen - r->mlen + r->p->n_ambi;
 	h->trans_strand = r->p->trans_strand == 1? 1 : r->p->trans_strand == 2? -1 : 0;
 	h->is_primary = (r->id == r->parent);
+	h->seg_id = r->seg_id;
 	h->n_cigar32 = r->p->n_cigar;
 	h->cigar32 = r->p->cigar;
 }
@@ -68,4 +70,64 @@ static inline void mm_reset_timer(void)
 	mm_realtime0 = realtime();
 }

+extern unsigned char seq_comp_table[256];
+static inline mm_reg1_t *mm_map_aux(const mm_idx_t *mi, const char *seq1, const char *seq2, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt)
+{
+	if (seq2 == 0) {
+		return mm_map(mi, strlen(seq1), seq1, n_regs, b, opt, NULL);
+	} else {
+		int _n_regs[2];
+		mm_reg1_t *regs[2];
+		char *seq[2];
+		int i, len[2];
+		len[0] = strlen(seq1);
+		len[1] = strlen(seq2);
+		seq[0] = (char*)seq1;
+		seq[1] = strdup(seq2);
+		for (i = 0; i < len[1]>>1; ++i) {
+			int t = seq[1][len[1] - i - 1];
+			seq[1][len[1] - i - 1] = seq_comp_table[(uint8_t)seq[1][i]];
+			seq[1][i] = seq_comp_table[t];
+		}
+		if (len[1]&1) seq[1][len[1]>>1] = seq_comp_table[(uint8_t)seq[1][len[1]>>1]];
+		mm_map_frag(mi, 2, len, (const char**)seq, _n_regs, regs, b, opt, NULL);
+		for (i = 0; i < _n_regs[1]; ++i)
+			regs[1][i].rev = !regs[1][i].rev;
+		*n_regs = _n_regs[0] + _n_regs[1];
+		regs[0] = (mm_reg1_t*)realloc(regs[0], sizeof(mm_reg1_t) * (*n_regs));
+		memcpy(&regs[0][_n_regs[0]], regs[1], _n_regs[1] * sizeof(mm_reg1_t));
+		free(regs[1]);
+		return regs[0];
+	}
+}
+
+static inline char *mappy_revcomp(int len, const uint8_t *seq)
+{
+	int i;
+	char *rev;
+	rev = (char*)malloc(len + 1);
+	for (i = 0; i < len; ++i)
+		rev[len - i - 1] = seq_comp_table[seq[i]];
+	rev[len] = 0;
+	return rev;
+}
+
+static char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int en, int *len)
+{
+	int i, rid;
+	char *s;
+	*len = 0;
+	rid = mm_idx_name2id(mi, name);
+	if (rid < 0) return 0;
+	if (st >= mi->seq[i].len || st >= en) return 0;
+	if (en < 0 || en > mi->seq[i].len)
+		en = mi->seq[i].len;
+	s = (char*)malloc(en - st + 1);
+	*len = mm_idx_getseq(mi, rid, st, en, s);
+	for (i = 0; i < *len; ++i)
+		s[i] = "ACGTN"[(uint8_t)s[i]];
+	s[*len] = 0;
+	return s;
+}
+
 #endif
@@ -26,11 +26,12 @@ cdef extern from "minimap.h":
 		int min_join_flank_sc
 		int a, b, q, e, q2, e2
 		int noncan
-		int zdrop
+		int zdrop, zdrop_inv
 		int end_bonus
 		int min_dp_max
 		int min_ksw_len
 		int anchor_ext_len, anchor_ext_shift
+		float max_clip_ratio
 		int pe_ori, pe_bonus
 		float mid_occ_frac
 		int32_t mid_occ
@@ -57,7 +58,7 @@ cdef extern from "minimap.h":
 		mm_idx_seq_t *seq
 		uint32_t *S
 		mm_idx_bucket_t *B
-		void *km
+		void *km, *h

 	ctypedef struct mm_idx_reader_t:
 		pass
@@ -68,6 +69,8 @@ cdef extern from "minimap.h":
 	void mm_idx_destroy(mm_idx_t *mi)
 	void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)

+	int mm_idx_index_name(mm_idx_t *mi)
+
 	#
 	# Mapping (key struct defined in cmappy.h below)
 	#
@@ -79,7 +82,6 @@ cdef extern from "minimap.h":

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

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

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

 	ctypedef struct kstring_t:
 		unsigned l, m
@@ -114,5 +119,6 @@ cdef extern from "cmappy.h":
 	void mm_fastx_close(kseq_t *ks)
 	int kseq_read(kseq_t *seq)

+	char *mappy_revcomp(int l, const uint8_t *seq)
 	int mm_verbose_level(int v)
 	void mm_reset_timer()
@@ -1,6 +1,7 @@
 from libc.stdint cimport uint8_t, int8_t
 from libc.stdlib cimport free
 cimport cmappy
+import sys

 cmappy.mm_reset_timer()

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

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

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

+	@property
+	def read_num(self): return self._seg_id + 1
+
 	@property
 	def cigar_str(self):
 		return "".join(map(lambda x: str(x[0]) + 'MIDNSH'[x[1]], self._cigar))
@@ -106,7 +112,7 @@ cdef class Aligner:
 		if min_chain_score is not None: self.map_opt.min_chain_score = min_chain_score
 		if min_dp_score is not None: self.map_opt.min_dp_max = min_dp_score
 		if bw is not None: self.map_opt.bw = bw
-		if best_n is not None: self.best_n = best_n
+		if best_n is not None: self.map_opt.best_n = best_n

 		cdef cmappy.mm_idx_reader_t *r;
 		if fn_idx_out is None:
@@ -117,6 +123,7 @@ cdef class Aligner:
 			self._idx = cmappy.mm_idx_reader_read(r, n_threads) # NB: ONLY read the first part
 			cmappy.mm_idx_reader_close(r)
 			cmappy.mm_mapopt_update(&self.map_opt, self._idx)
+			cmappy.mm_idx_index_name(self._idx)

 	def __dealloc__(self):
 		if self._idx is not NULL:
@@ -125,7 +132,7 @@ cdef class Aligner:
 	def __bool__(self):
 		return (self._idx != NULL)

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

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

-def fastx_read(fn):
+	def seq(self, str name, int start=0, int end=0x7fffffff):
+		cdef int l
+		cdef char *s = cmappy.mappy_fetch_seq(self._idx, name.encode(), start, end, &l)
+		if l == 0: return None
+		r = s[:l] if isinstance(s, str) else s[:l].decode()
+		free(s)
+		return r
+
+	@property
+	def k(self): return self._idx.k
+
+	@property
+	def w(self): return self._idx.w
+
+	@property
+	def n_seq(self): return self._idx.n_seq
+
+def fastx_read(fn, read_comment=False):
 	cdef cmappy.kseq_t *ks
 	ks = cmappy.mm_fastx_open(str.encode(fn))
 	if ks is NULL: return None
@@ -155,9 +185,22 @@ def fastx_read(fn):
 		else: qual = None
 		name = ks.name.s if isinstance(ks.name.s, str) else ks.name.s.decode()
 		seq = ks.seq.s if isinstance(ks.seq.s, str) else ks.seq.s.decode()
-		yield name, seq, qual
+		if read_comment:
+			if ks.comment.l > 0: comment = ks.comment.s if isinstance(ks.comment.s, str) else ks.comment.s.decode()
+			else: comment = None
+			yield name, seq, qual, comment
+		else:
+			yield name, seq, qual
 	cmappy.mm_fastx_close(ks)

+def revcomp(seq):
+	l = len(seq)
+	bseq = seq if isinstance(seq, bytes) else seq.encode()
+	cdef char *s = cmappy.mappy_revcomp(l, bseq)
+	r = s[:l] if isinstance(s, str) else s[:l].decode()
+	free(s)
+	return r
+
 def verbose(v=None):
 	if v is None: v = -1
 	return cmappy.mm_verbose_level(v)
@@ -23,7 +23,7 @@ def readme():

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

 unsigned char seq_nt4_table[256] = {
 	0, 1, 2, 3,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
@@ -0,0 +1,127 @@
+>ref
+TGCGGAGGCTGAAGCAACTCCATCTTGGAAGCTAATCTACCATGTTGGCTTCTGATTAAC
+ATCAGTTCTGGGAAGGCTTGTAAGATTTCCTGTTTGTCTATTATTTCCTAGGTAAGAGCA
+GATACTTACTGTAAATCCTGCCCCTAGATTAAACAACCTTGGTGTTATCGTACTTCCATT
+GTCCTATACATCCCTTCGGAATCCCCCTTTCCCTATGGTCCTCAAGCCCTTGGTCTGGGG
+AGTAACAGCATAGGGATCAACCATCTCGTCTTGCCACTGCCCGAAATACAGACATGGCTT
+CTGTTCCTAAGTCCCTATTCAACTTTTCTTTCTAAGAAACTGGATTTGTCAGCCTCTTTC
+TTCACCTCTCAGCTTCCTTGGACTTTGGGGGTAGGTTTGCGTAGACATGCTCACCACAGA
+CACAATATCAGCTTCATTCTACAGATGAGGAAGGCAAGCCTTGGGGAGCTTAACCAACTT
+GTCGAGACTCATGTATATACCAACACTGAAAAGCAGATATTCCAGACTCCCAGTCATGCC
+ACAGGCACACCCCTCAGTGAGAGGTGGGGTTTGTAGTTGAGGCTATTTCCTGCCCAGGGA
+GCAGGGAGGCACTCTAGCTTCCCTGAGCTAACGTGGTTCTGCTTGTGTCTGACTTCCAGG
+TCTCTGCCCTTTCCAAGCTCACTAGGATGGGCTTCGGGTGTGTCAAATGCCTCAGACAGT
+ACAGATCCACACAGAATGGGCATATGCAACCAATCAGTGTCATAAAAAAGAAGGAAATGA
+CTCGGGCCCCCTGTGTGTTCAACATGTCGAAGGTATCTGTGCAGCAGAAGAAAGAGGGGC
+AAAAGCCCCCAGTGCCACAGGCCAGAGGCAGCAGCTTGGGCCCATGTGGGAGGGTTTGCT
+TTCCCCTGCCAAAGTGATGGGCTGCTGCAGCCTGGGGCTTGTGGGAATCCTTCCTGGGCC
+TGTGTGGGAAGTGTAGGCAGGGAGAGTGCTGCTTTCCCAAGCTCATCCCAGCTACAGCTA
+CCTTTGTGCTCTGGGATTCAGGACCCCCGAGGGGGCTGGCAGGAGAGTCTCTGTTCTCGG
+ATGGGTTGTCACCAGGGCATACATGGGAAGTGGGCTCTCTGGAGTCACCCTCCAGGGGAC
+AATGCCAATTCCAGACACATTTACTGGAACCCCTACACTGATGACCTTTTGTTGAGGGTT
+GAATTATGTCCCCAAAAAAGATACATTGAAGTCCAAACCTCTGGTGTCTATAAATGTGAT
+TTTATTTGAAAATGAGGTTTCTATGGACTAAATTGTGTCCCTCCCAAATTCATATTTTGA
+AGCCCTAGCCCCCAGTGTGACTATACCTAGAGACAGAGATCTTTAGGAGGTAATTAAGGT
+TCAATGAGGTCAGGTGGGTGGGGCCCTAAACCAACAGGAAGGACTGTGGCCTTACTAGAA
+AAGGAAGAAAAAGCATTTCCTCTCTTCTAGTATAAAAGGACACAGAAAGAAGGCAGATAT
+CTACAAGCCACGAAGAGAGACGTCACTGAGAACTGAATTTGTGTACATTGATCTGGAACT
+TCCAGCCTCCAGAACTTGAGAAATACATTTCTGTTGTTTATTTTTTTTTCATGTAATCAA
+TTCATTTATCATATATTTATTGAGTGCCTACTATGTGCCAGAGGATACAGCAGTAACAAA
+ACTAGGCAAAAATTGTGCCTAAAAGAGGGAAGATGACTTTTCTTAAAGTGTGGAATAAAG
+AAAAGTAAGATAGCGGATAGAAGCTTGAAGTGAAAGCAGGTTCACAGGAAGTTTCTTTGG
+TCATTTGTTTTGTTTTTAAATAGTGGAAAGATGTATATGTTTATGGAGAAAGATTGCCTT
+GAAGATGCAAGAGGAAGAGATGATCAAAATTCAAGAAGAAGCAGAAAGTGATAGAATAAA
+GAGCACAAGTGGAGAATTAGTGTTAATGAAAAGAAGGATGCTTCCTTTGATATGAAGTGA
+AGGAAGAGAGAATGAGTAAAGACCAAGACTTGAAGTCCCTAGTTTAATAGAGGGAGATTT
+CTTCTTTTGATAGCAACAATGGTATTCTGAATTATTTGAAGACATGTCATATTTCTCTTG
+TGCCATTTTCCTCCCAGTTTAAACATTCTCATAACCTCTATTCCTCACATGATGTTTTTC
+CAGGTCCTTTATTCTTTGGCACTCTCTTCTCTGGACACATTGTATTCTGTCATTGGTCCT
+AAAATTTAGATACCCACAATTGAACATACTCCTCTAGATATGGTCTAGCTAATGCAAAAG
+AACTGCTGCCTTCCAACTTGTTCAGACATCATATGTTTGTTGTCAAACGCTAAGTTGAGT
+TGTTATCTTTTAAGTTTTGTTTTTGTTTTTTTTTTTTTTTTTTAATTCCAAGAGGTGCCC
+ACGTTGGCTAAGTACCAAACAGGGTACTAGGGAATTTTACTTCTGAGTTAAATGCCATTC
+TAGTTGTTTTTTCTTCATCTCCAGTAAGGTTATCTTTATTCACCAGTTGTTACAATAGCT
+GTGGGTCTTGCTTCTCACAGTTTTATGCTGTCTGTGCTATTTTCTCTACTGATCATCACC
+ACAATCATTATTGCTTATCATAATTGTTATCTTTATTTTCTCCTTTAATCAAGAATCAGT
+CTTCCTTTATCTCATTATTCTCTTTTGCAGGCTTCAGGATAATTATGGTTGGAGTGCACT
+GGGGGAACCAGTGCAGCTAAGCTCTGACATCTTTGCATCCCTTTTCCATCTGCTGTTTTG
+GCACTCTGGTAGAATAGATAACCTAAAAACGACTTTAAAACATCTAGAAATTTTGGATAA
+AATATAACAAACATCCCTTTAAATGCACAACTGATCTTCCATGGAAGTCACAGAAATATA
+TAACGCCAAAAAGAAGGGAAGCTGAAACCCAGGGCTGTAAACATGAACATCATCTTCTCT
+CCCTTTTTCTTGTGACTTATCTTGTTTTTCTCAGCTTTGGTGCTACCAAGGCTTGACTTT
+AATAGGCATTTCCAATCAATGAGAGAATTTCTTTTGCTTTCATCAACAATTCAGTTATTG
+ATGTTAACATATATATCATTTGAGTACTTTTCTTTTTTTTATTATTATTATACTTTAAGT
+TTTAGGGTCCATGTGCACAATGTGCAGGTTAGTTACGTATGTATACATGTGCCATGCTGG
+TGTGCTGCACCCATTAACTCATCATTTAGCATTAGGTATATCTCCTAATGCTATCCCTTC
+CCCCTCTCCCCACCCCACAACAGTCCCCAGAGTGTTCCCCTTCCTGTGTCCATGTGTTCT
+CATTGTTCAATCCCCATCTATGAGTGAGAACATGCGGTGTTTGGTTTTTTGTCCTTGCAA
+TAGTTTACTGAGAATGATGATTTCTAATTTCATCCATGTCCCTAAAGAGCTTCTGCACAG
+CAAAAGAAACTACCATCAGAGTGAACAGGCAACCTACAAAATGGGAGAAAATTTTCACAA
+CCTGCTCATCTGACAAAGGGCTAATATCCAGAATCTACAATGAACTCAAACAAATTTACA
+AGAAAAAAACAAACAACCCCATCAAAAAGTGGGCAAAGGATATGAACAGACACTTCTCAA
+AAGAAGACATTTATGCAGCCAAAAGACACATGAAAAAATGCTCATCATCACTGGCCATCA
+GAGAAATGCAAACCAAAACCACAATGAGATACCATCTCACACCAGTTAAAATGGCAATCA
+TTAAAAAGTCAGGAAACAACAGGTGCTGGAGAGGATGTGGAGAAACAGGAACACTTTTAC
+ACTGTTGGTGGGACTGTAAACTAGTTCAACCATTGTGGAAGTCAGTGTGCTGATTCCTCA
+GGGATCTAGAACTAGAAATACCATTTGACCCAGCCATCCCATTACTGGGTATATACCCAA
+AGGACTATAAATCATGCTGCTATAAAGACACATGCACACGTATGTTTATTGCGGCACTAT
+TCACAATAGCAAAGACTTGGAACCAACCCAAATGTCCAACAATGATAGACTGGATTAAGA
+AAATGTGGCACATATACACCACGGAATACTGTGCAGCCATAAAAAATGATGAGTTCATGT
+CCTTTGTAGGGACACGGATGAAATTGGAAATCATTTCTGTTGTTTAAACCACGAAGTCTA
+TGGTATCTGGTTATGACAACCTGAGAATACTAACTCAAGGGTCTTTCGCAGATGTCATTA
+AGTTGTTAAAGTGAGGTCATTATGGTGGGTCCTAATCCAAGAGAAGAGATGCATGGACAG
+ACGTGCACAACGGGAGGACCAAGCCAAGACACACAGGGAGAATGGCCATGGGAAGATGGA
+GGCAGAGATCAAAGTGAGGCACCCACAAGCCAAGAAATGGCAGGAGCTACCAGCAGCTGG
+AAGATGCAGAGAAGCATTCCTTCTTAGAGGTTTCAGAGAGAGTATGGTGCTACTGACACC
+TTGATTTTGAACTTCTAGTCTCCAGAACTATGAGAGAATAAATTTCTGTTGGTTAAGCCA
+TCGAGTTTGTGTAAGTTTGTTATAAGAGCCCTAGGAAATAAACATATCCATTTATTCAGG
+AAAGCCTGCTAGAGTGCAAATATTTGGAAAAGATACTACTATGCAAATGTTTGAAAAAGA
+TATTGCTCTTGATTCTGCCTTATGGGTTTTTCATTTCTGTAAGCTATTCTCAAAGTTTTG
+TTCTTGGACTACTATTGGTAATTAAGACTGCAACATGTTTGGCAACATCAGTTGAGAACT
+GTTGCTCTGGGAACGTTTTCGGCAAGCCTCAGCCCTTCTTTTCCCTTGGCTTGCATTGAG
+GAGTTAGGTGATACTCTGCTGCTCAGGCCCAGCACCTTTATGGACCGTATTCCCCTGGTG
+GAATGACCATCTCTGCTTGCTCTGATTGGCTGTTGGGGTTTTCTAGCATGCCCTATTTAA
+TATGTATGATTTATCTCTTACTTCAGTTGGAAGGTACAGTTGCTCTGTAGTTGGCATGCA
+GTCATGGTGACTATGAAAATATAAAATAATGTTTTGGTTTACAGACACTTAGAAATAAGT
+TGTGTCTCAAAATTGGGTGACTATTCTAGTTATCTGCTACTCAATATCCTTGTGCGAGCC
+CTCTTTACCCAGAATCAAACTAAACCATGAGGGGCACTATAGAATGTCACCCCTGGGTCC
+AGGATACTATGGGGACTCAGAAGCCAAGCTCCCACTGGGGGATCTAGGGCATGCCCCCAA
+GGTAAGATTCCCACCTCTTTGTTCAGCAGGAAGCACCCATCACACAAGGAGGTAGGAATA
+AACAAGCATTCGTCAAGAACAAAAGATACAGATGTTCTGCTGGAGCTTGGATACATAGCA
+TAAGAGGGAACAGTTCTCACAGGTAAGAGTAAGTTTTCCTCTGGTGGTGACAGTGGGACC
+TGTGGGGGAGAGAATTGGGAGTACTGACAGGAAGGCAGAGTGGCTGTCCAAATGAACGGA
+TTGTTTGCACATGGCCTTTAGGGCACGTTGTGTTAGCCTTCCATTGCTGCTTATATTAGT
+CTGTTTTCACACTGCCCATAAATGCATACCTGAGACTGGATAATTTATAAAGAAAAAGAG
+CCTTAATGTACTCATAGTTGCATGTGGCTGGGGAGGCCTCACAATCATGGCAGAAGGTGA
+AAGGCACATCTTACATGGAAGCAGACAAGAGAGAATTGAGGACCAAGTGAAAGGGGTTTC
+CCCTTATAAAACCATCAGATCACATGAGACTTTTTCACCACCATGAGAACAGTAAGGGGA
+AAACTATGCTCATGATTCAATTGTCTCCCACTGGATTCCTCCCACAACACATAGGAATTA
+TGGGAGCTAAAATTCAAGATGAGATTTGGGTGAGGACACAGCCAAACCCTATCACTGCTG
+TAATCAATTCCCACCAACTTAGTGGCTCGAAACATCACAGATTTATGATCTTATGACGGT
+GGAGGTCCCCAAATGGATCTTCTAGGTCTAGAATCAAGGTATCAGCAGACCACTTCTTTT
+GGAGGCTCTGGTGGAGAAACCATTTCCTCGCCTTTTCCAGCTTCTAGAGGCTGCCCTTCT
+CATTCCTTGGTTCACGGCCACACTCATTTCCATCTCTGCTTCCACTGTGACAACTTCTCT
+GCCTCAGACCCTCCTGCTTTGCCTTTGTAAGGACCCTTGTGATGAGATCAGGCCCATCCA
+GGATTATCCCTCATCTCAAGACCTTTACCTTAATCACATTTGCAAGGTCTCTTCCACTGT
+GTCAGGTAACATTTTCACAGGTTCCAGGGATTAGGGTGTGGACATCTTGGGGAGCTGGAG
+GATATTATTTCATCTACCACACACATCTCTACCTTGTACAGGCAAGCACTTGCAAAGTGC
+AATGTGATCCTCTGGAGCCACTGTCCTCCCAGAGCTTATATATACTCTGAAAGTCAACTC
+TCAGACCACAGCCTCCTGTCCATGCACCACTCTCATCAACACCCCCACCCGAAACACTTT
+CACTCCACCCTCTTTGTCCCCTAACTCATGGAGAAGAAAATCTAATTAGTAGGAGTGGAA
+TTTGGCTTTCATCTTTACCAGTACTAGAAATATGGTGTGTGTCTTTTTGTAAAAATTCTC
+TCAACTAAATTGTTTTTATTAATTTCTGCAAAATGTGAACATCAACTCCCTTCATGTGAA
+TGTCAATAAGATTAAATGAGCTGTCTCAGCTCCTAGCCTGTGCAAGCTAACAGCTCAGGA
+GATGTTTATTTCTTTCCCTCTTCTTTCCTTAATGAAGCCCTCTCCTTTGACATCTTCAAT
+TCTGGAGCGCTTCTTTTCTGAGGCCTTGGCTCCCCCACATTGCCCACCCTTTTCCTGCTC
+GTCCACATTTCTGGCTTCTATTCTCTTGTCTTTACCATCTCCCTGAACAATGTTATCCGT
+TCCAATGACTTCAACAGTCTCTCCGCTTACATATGATGCCTCTCAAACTCTGATCTCCAA
+CTCTTCCAAAGAGCTCTGGACCTTTGTTCCAATTACCTGAAAAACATCTTCTTGGATGTC
+CCATTAGCACTGTTAAATCAAACAAGAATTTCCCTCCCTCCTGCCTTGCTGTAGTTCCCC
+TAGGGATTCGGTTGTGTGGGAAGATGTGTGGAGAGCTCTTAGTTGACTCCCTTCTCTGCA
+GTTCTACCTCTCTAGAGACTTGGAGGACCCACTGTTTCCGCCTCGCTTTTTCAGGCCTAG
+AGATTGCTCGCTCCTGGGCTGGCTGCTTCATAATTCCTTATTAGTAGTTTCCCAAGCTTA
+CATATCTGTAAATATTTACTTTAGTTAAATTCTCCCCAATTTCCACAATATGTTGGCTGC
+ACATGCTTTCTACTAGGAGTCACACAACTATGATAAGAACCAAGAAATATTAGTAAACGT
+TTTTTACCATTATTGGCCTATACCCTGGAATAGCCAACAATAACCTAGAACCTATGCAAC
+AAGAATATCCAACAAGAACCTAGAGACCTGTCAGTCTATAGGTGGGAACTACAGGATGAG
+A
@@ -313,3 +313,27 @@
 	Title = {Assembling large genomes with single-molecule sequencing and locality-sensitive hashing},
 	Volume = {33},
 	Year = {2015}}
+
+@article{Gurevich:2013aa,
+	Author = {Gurevich, Alexey and others},
+	Journal = {Bioinformatics},
+	Pages = {1072-5},
+	Title = {{QUAST}: quality assessment tool for genome assemblies},
+	Volume = {29},
+	Year = {2013}}
+
+@article{Li:2010fk,
+	Author = {Li, Heng and Durbin, Richard},
+	Journal = {Bioinformatics},
+	Pages = {589-95},
+	Title = {Fast and accurate long-read alignment with {Burrows-Wheeler} transform},
+	Volume = {26},
+	Year = {2010}}
+
+@article{Marcais:2018aa,
+	Author = {Mar{\c c}ais, Guillaume and others},
+	Journal = {PLoS Comput Biol},
+	Pages = {e1005944},
+	Title = {{MUMmer4}: A fast and versatile genome alignment system},
+	Volume = {14},
+	Year = {2018}}
@@ -1,6 +1,6 @@
 \documentclass{bioinfo}
-\copyrightyear{2017}
-\pubyear{2017}
+\copyrightyear{2018}
+\pubyear{2018}

 \usepackage{graphicx}
 \usepackage{hyperref}
@@ -40,9 +40,10 @@ full-length noisy Direct RNA or cDNA reads, and assembly contigs or closely
 related full chromosomes of hundreds of megabases in length. Minimap2 does
 split-read alignment, employs concave gap cost for long insertions and
 deletions (INDELs) and introduces new heuristics to reduce spurious alignments.
-It is 3--4 times faster than mainstream short-read mappers at comparable
-accuracy and $\ge$30 times faster at higher accuracy for both genomic and mRNA
-reads, surpassing most aligners specialized in one type of alignment.
+It is 3--4 times as fast as mainstream short-read mappers at comparable
+accuracy, and is $\ge$30 times faster than long-read genomic or cDNA
+mappers at higher accuracy, surpassing most aligners specialized in one type of
+alignment.

 \section{Availability and implementation:}
 \href{https://github.com/lh3/minimap2}{https://github.com/lh3/minimap2}
@@ -88,11 +89,11 @@ the versatility of minimap2.
 Minimap2 follows a typical seed-chain-align procedure as is used by most
 full-genome aligners. It collects minimizers~\citep{Roberts:2004fv} of the
 reference sequences and indexes them in a hash table. Then for each query
-sequence, minimap2 takes query minimizers as \emph{seeds}, finds matches to the
-reference, and identifies sets of colinear seeds, which are called
+sequence, minimap2 takes query minimizers as \emph{seeds}, finds exact matches
+(i.e. \emph{anchors}) to the reference, and identifies sets of colinear anchors as
 \emph{chains}. If base-level alignment is requested, minimap2 applies dynamic
-programming (DP) to extend from the ends of chains and to close unseeded
-regions between adjacent seeds in chains.
+programming (DP) to extend from the ends of chains and to close
+regions between adjacent anchors in chains.

 Minimap2 uses indexing and seeding algorithms similar to
 minimap~\citep{Li:2016aa}, and furthers the predecessor with more accurate
@@ -117,12 +118,12 @@ distance between two anchors is too large); otherwise
 \begin{equation}\label{eq:chain-gap}
 \beta(j,i)=\gamma_c\big((y_i-y_j)-(x_i-x_j)\big)
 \end{equation}
-In implementation, a gap of length $l$ costs
+In implementation, a gap of length $l\not=0$ costs
 \[
 \gamma_c(l)=0.01\cdot \bar{w}\cdot|l|+0.5\log_2|l|
 \]
-where $\bar{w}$ is the average seed length. For $m$ anchors, directly computing all $f(\cdot)$ with
-Eq.~(\ref{eq:chain}) takes $O(m^2)$ time. Although theoretically faster
+where $\bar{w}$ is the average seed length. For $N$ anchors, directly computing all $f(\cdot)$ with
+Eq.~(\ref{eq:chain}) takes $O(N^2)$ time. Although theoretically faster
 chaining algorithms exist~\citep{Abouelhoda:2005aa}, they
 are inapplicable to generic gap cost, complex to implement and usually
 associated with a large constant. We introduced a simple heuristic to
@@ -132,7 +133,7 @@ We note that if anchor $i$ is chained to $j$, chaining $i$ to a predecessor
 of $j$ is likely to yield a lower score. When evaluating Eq.~(\ref{eq:chain}),
 we start from anchor $i-1$ and stop the process if we cannot find a better
 score after up to $h$ iterations. This approach reduces the average time to
-$O(h\cdot m)$. In practice, we can almost always find the optimal chain with
+$O(hN)$. In practice, we can almost always find the optimal chain with
 $h=50$; even if the heuristic fails, the optimal chain is often close.

 \subsubsection{Backtracking}
@@ -146,9 +147,11 @@ in more than one chains.
 \subsubsection{Identifying primary chains}\label{sec:primary}
 In the absence of copy number changes, each query segment should not be mapped
 to two places in the reference. However, chains found at the previous step may
-have significant or complete overlaps due to repeats in the reference.
+have significant or complete overlaps due to repeats in the reference~\citep{Li:2010fk}.
 Minimap2 used the following procedure to identify \emph{primary chains} that do
-not greatly overlap on the query. Let $Q$ be an empty set initially. For each
+not greatly overlap on the query.
+
+Let $Q$ be an empty set initially. For each
 chain from the best to the worst according to their chaining scores: if on the
 query, the chain overlaps with a chain in $Q$ by 50\% or higher percentage of
 the shorter chain, mark the chain as secondary to the chain in $Q$; otherwise,
@@ -156,6 +159,16 @@ 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.

+For each primary chain, minimap2 estimates its mapping quality with an
+empirical formula:
+\[
+{\rm mapQ}=40\cdot (1-f_2/f_1)\cdot\min\{1,m/10\}\cdot\log f_1
+\]
+where $m$ is the number of anchors on the primary chain, $f_1$ is the chaining
+score, and $f_2\le f_1$ is the score of the best chain that is secondary to the
+primary chain. Intuitively, a chain is assigned to a higher mapping quality if
+it is long and its best secondary chain is weak.
+
 \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$
@@ -186,7 +199,7 @@ $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)
+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.

@@ -199,9 +212,9 @@ 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
+alignment. With HPC 19-mers per 5bp window, 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.
+performance. HPC-based indexing reduces the sensitivity for current ONT reads, though.

 \subsection{Aligning genomic DNA}\label{sec:genomic}

@@ -246,8 +259,8 @@ In case of 2-piece gap cost, define
 \[
 \left\{\begin{array}{ll}
 u_{ij}\triangleq H_{ij}-H_{i-1,j} & v_{ij}\triangleq H_{ij}-H_{i,j-1} \\
-x_{ij}\triangleq E_{i+1,j}-H_{ij} & \tilde{x}_{ij}\triangleq \tilde{E}_{i+1,j}-\tilde{H}_{ij} \\
-y_{ij}\triangleq F_{i,j+1}-H_{ij} & \tilde{y}_{ij}\triangleq \tilde{F}_{i,j+1}-\tilde{H}_{ij}
+x_{ij}\triangleq E_{i+1,j}-H_{ij} & \tilde{x}_{ij}\triangleq \tilde{E}_{i+1,j}-H_{ij} \\
+y_{ij}\triangleq F_{i,j+1}-H_{ij} & \tilde{y}_{ij}\triangleq \tilde{F}_{i,j+1}-H_{ij}
 \end{array}\right.
 \]
 We can transform Eq.~(\ref{eq:ae86}) to
@@ -311,7 +324,7 @@ the same anti-diagonal in one inner loop. It also simplifies banded alignment,
 which would be difficult with striped vectorization~\citep{Farrar:2007hs}.

 On the condition that $q+e<\tilde{q}+\tilde{e}$ and $e>\tilde{e}$, the initial
-values in the diagonal-antidiagonal formuation is
+values in the diagonal-antidiagonal formuation are
 \[
 \left\{\begin{array}{l}
 x_{r-1,-1}=y_{r-1,r}=-q-e\\
@@ -330,6 +343,13 @@ r\cdot(e-\tilde{e})-(\tilde{q}-q)-\tilde{e} & (r=\lceil\frac{\tilde{q}-q}{e-\til
 \]
 These can be derived from the initial values for Eq.~(\ref{eq:ae86}).

+When performing global alignment, we do not need to compute $H_{rt}$ in each cell.
+We use 16-way vectorization throughout the alignment process. When extending
+alignments from ends of chains, we need to find the cell $(r,t)$ where $H_{rt}$
+reaches the maximum. We resort to 4-way vectorization to compute
+$H_{rt}=H_{r-1,t}+u_{rt}$. Because this computation is simple,
+Eq.~(\ref{eq:suzuki}) is still the dominant performance bottleneck.
+
 In practice, our 16-way vectorized implementation of global alignment is three
 times as fast as Parasail's 4-way vectorization~\citep{Daily:2016aa}.  Without
 banding, our implementation is slower than Edlib~\citep{Sosic:2017aa}, but with
@@ -359,6 +379,16 @@ alignment between the two subsequences involved in the global alignment, but
 this time with the one subsequence reverse complemented. This additional
 alignment step may identify short inversions that are missed during chaining.

+\subsubsection{Filtering out misplaced anchors}
+Due to sequencing errors and local homology, some anchors in a chain may be
+wrong. If we blindly align regions between two misplaced anchors, we will
+produce a suboptimal alignment. To reduce this artifact, we filter out
+anchors that lead to a $>$10bp insertion and a $>$10bp deletion at the same
+time, and filter out terminal anchors that lead to a long gap towards the ends
+of a chain. These heuristics greatly alleviate the issues with misplaced
+anchors, but they are unable to fix all such errors. Local misalignment is a
+limitation of minimap2 which we hope to address in future.
+
 \subsection{Aligning spliced sequences}

 The algorithm described above can be adapted to spliced alignment. In this
@@ -397,7 +427,7 @@ p/2 & \mbox{if $T[i+1,i+3]$ is ${\tt GTC}$ or ${\tt GTT}$} \\
 p & \mbox{otherwise}
 \end{array}\right.\]
 where $T[i,j]$ extracts a substring of $T$ between $i$ and $j$ inclusively.
-$d(i)$ penalizes non-canonical donor sites with $p$ and less frequent Eukayotic
+$d(i)$ penalizes non-canonical donor sites with $p$ and less frequent Eukaryotic
 splicing signal ${\tt GT[C/T]}$ with $p/2$~\citep{Irimia:2008aa}. Similarly,
 \[a(i)=\left\{\begin{array}{ll}
 0 & \mbox{if $T[i-2,i]$ is ${\tt CAG}$ or ${\tt TAG}$} \\
@@ -424,13 +454,13 @@ alignment.

 \subsection{Aligning short paired-end reads}

-During chainging, minimap2 takes a pair of reads as one fragment with a gap of
+During chaining, minimap2 takes a pair of reads as one fragment with a gap of
 unknown length in the middle. It applies a normal gap cost between seeds on the
 same read but is a more permissive gap cost between seeds on different reads.
-More precisely, the gap cost during chaining is:
+More precisely, the gap cost during chaining is ($l\not=0$):
 \[
 \gamma_c(l)=\left\{\begin{array}{ll}
-0.01\cdot\bar{w}\cdot l+0.5\log_2 l & \mbox{if two seeds on the same read} \\
+0.01\cdot\bar{w}\cdot |l|+0.5\log_2 |l| & \mbox{if two seeds on the same read} \\
 \min\{0.01\cdot\bar{w}\cdot|l|,\log_2|l|\} & \mbox{otherwise}
 \end{array}\right.
 \]
@@ -443,6 +473,17 @@ consistent paired-end alignments.

 \section{Results}

+Minimap2 is implemented in the C programming language and comes with APIs in
+both C and Python. It is distributed under the MIT license, free to both
+commercial and academic uses. Minimap2 uses the same base algorithm for all
+applications, but it has to apply different sets of parameters depending on
+input data types. Similar to BWA-MEM, minimap2 introduces `presets' that
+modify multiple parameters with a simple invokation. Detailed settings
+and command-line options can be found in the minimap2 manpage. In addition to
+the applications evaluated in the following sections, minimap2 also retains
+minimap's functionality to find overlaps between long reads and to search
+against large multi-species databases such as \emph{nt} from NCBI.
+
 \subsection{Aligning long genomic reads}\label{sec:long-genomic}

 \begin{figure}[!tb]
@@ -450,19 +491,22 @@ consistent paired-end alignments.
 \includegraphics[width=.5\textwidth]{roc-color.pdf}
 \caption{Evaluation on aligning simulated reads. Simulated reads were mapped 
 to the primary assembly of human genome GRCh38. A read is considered correctly
-mapped if the true position overlaps with the best mapping position by 10\% of
-the read length. Read alignments are sorted by mapping quality in the
-descending order. For each mapping quality threshold, the fraction of
-alignments with mapping quality above the threshold and their error rate are
+mapped if its longest alignment overlaps with the true interval, and the
+overlap length is $\ge$10\% of the true interval length.  Read alignments are
+sorted by mapping quality in the descending order. For each mapping quality
+threshold, the fraction of alignments (out of the number of input reads) with
+mapping quality above the threshold and their error rate are
 plotted along the curve. (a) long-read alignment evaluation. 33,088 $\ge$1000bp
 reads were simulated using pbsim~\citep{Ono:2013aa} with error profile sampled
 from file `m131017\_060208\_42213\_*.1.*' downloaded at
 \href{http://bit.ly/chm1p5c3}{http://bit.ly/chm1p5c3}. The N50 read length is
 11,628. Aligners were run under the default setting for SMRT reads.
-(b) short-read alignment evaluation. 10 million pairs of 150bp reads were
-simulated using mason2~\citep{Holtgrewe:2010aa} with option
-`\mbox{--illumina-prob-mismatch-scale 2.5}'. Short-read aligners were run under the
-default setting except for changing the maximum fragment length to
+Kart outputted all alignments at mapping quality 60, so is not shown in the
+figure. It mapped nearly all reads with 4.1\% of alignments being wrong, less
+accurate than others.  (b) short-read alignment evaluation. 10 million pairs of
+150bp reads were simulated using mason2~\citep{Holtgrewe:2010aa} with option
+`\mbox{--illumina-prob-mismatch-scale 2.5}'. Short-read aligners were run under
+the default setting except for changing the maximum fragment length to
 800bp.}\label{fig:eval}
 \end{figure}

@@ -471,7 +515,7 @@ BLASR~(v1.MC.rc64; \citealp{Chaisson:2012aa}),
 BWA-MEM~(v0.7.15; \citealp{Li:2013aa}),
 GraphMap~(v0.5.2; \citealp{Sovic:2016aa}),
 Kart~(v2.2.5; \citealp{Lin:2017aa}),
-minialign~(v0.5.3; \citealp{Suzuki:2016}) and
+minialign~(v0.5.3; \href{https://github.com/ocxtal/minialign}{https://github.com/ocxtal/minialign}) and
 NGMLR~(v0.2.5; \citealp{Sedlazeck169557}). We excluded rHAT~\citep{Liu:2016ab}
 and LAMSA~\citep{Liu:2017aa} because they either
 crashed or produced malformatted output. In this evaluation, minimap2 has
@@ -480,11 +524,11 @@ higher mapping accuracy (Fig.~\ref{fig:eval}a). Minimap2 and
 NGMLR provide better mapping quality estimate: they rarely give repetitive hits
 high mapping quality.  Apparently, other aligners may
 occasionally miss close suboptimal hits and be overconfident in wrong mappings.
-On run time, minialign is slightly faster than minimap2 and Kart. They are over
-30 times faster than the rest.  Minimap2 consumed 6.1GB memory at the peak,
-more than BWA-MEM but less than others.
+On run time, minimap2 took 200 CPU seconds, comparable to minialign and Kart, and is over
+30 times faster than the rest.  Minimap2 consumed 6.8GB memory at the peak,
+more than BWA-MEM (5.4GB), similar to NGMLR and less than others.

-On real human SMRT reads, the relative performance and sensitivity of
+On real human SMRT reads, the relative performance and fraction of mapped reads reported by
 these aligners are broadly similar to the metrics on simulated data. We are
 unable to provide a good estimate of mapping error rate due to the lack of the
 truth.  On ONT $\sim$100kb human reads~\citep{Jain128835}, BWA-MEM failed.
@@ -526,7 +570,7 @@ Peak RAM (GByte)          & 8.9      & 14.5     & 3.2      & 29.2\vspace{1em}\\
 \% approx. introns        & 91.8\%   & 96.9\%   & 92.5\%   & 82.4\% \\
 \botrule
 \end{tabular}
-}{Mouse reads (AC:SRR5286960) were mapped to the primary assembly of mouse
+}{Mouse reads (AC:SRR5286960; R9.4 chemistry) were mapped to the primary assembly of mouse
 genome GRCm38 with the following tools and command options: minimap2 (`-ax
 splice'); GMAP (`-n 0 --min-intronlength 30 --cross-species'); SpAln (`-Q7 -LS
 -S3'); STARlong (according to
@@ -592,7 +636,7 @@ simulated data set than Bowtie2 and SNAP but less accurate than BWA-MEM
 (Fig.~\ref{fig:eval}b). Closer investigation reveals that BWA-MEM achieves
 a higher accuracy partly because it tries to locally align a read in a small
 region close to its mate. If we disable this feature, BWA-MEM becomes slightly
-less accurate than minimap2. We might consider to implement a similar heuristic
+less accurate than minimap2. We might implement a similar heuristic
 in minimap2 in future.

 To evaluate the accuracy of minimap2 on real data, we aligned human reads
@@ -608,14 +652,24 @@ million bases (FPPM; 3.0 vs 3.9), lower 2--50bp INDEL FNR (7.3\% vs 7.5\%) and
 similar INDEL FPPM (both 1.0). Minimap2 is broadly similar to BWA-MEM in the
 context of small variant calling.

-\subsection{Other applications}
+\subsection{Aligning long-read assemblies}

-Minimap2 retains minimap's functionality to find overlaps between long reads
-and to search against large multi-species databases such as \emph{nt} from
-NCBI. Minimap2 can also align similar genomes or different assemblies of the
-same species. It took 7 wall-clock minutes over 8 CPU cores to align a human
-SMRT assembly (AC:GCA\_001297185.1) to GRCh38, over 20 times faster
-MUMmer4~\citep{Kurtz:2004zr}.
+Minimap2 can align a SMRT assembly (AC:GCA\_001297185.1) against GRCh38 in 7
+minutes using 8 CPU cores, over 20 times faster than nucmer from
+MUMmer4~\citep{Marcais:2018aa}. With the paftools.js script from the minimap2
+package, we called 2.67 million single-base substitutions out of 2.78Gbp
+genomic regions. The transition-to-transversion ratio (ts/tv) is 2.01. In
+comparison, using MUMmer4's dnadiff pipeline, we called 2.86 million
+substitutions in 2.83Gbp at ts/tv=1.87. Given that ts/tv averaged across the
+human genome is about 2 but ts/tv averaged over random errors is 0.5, the
+minimap2 callset arguably has higher precision at lower sensitivity.
+
+The sample being assembled is a female. Minimap2 still called 201 substitutions
+on the Y chromosome. These substitutions all come from one contig aligned at
+96.8\% sequence identity. The contig could be a segmental duplication
+absent from GRCh38. In constrast, dnadiff called 9070 substitutions on the Y
+chromosome across 73 SMRT contigs. This again implies our minimap2-based
+pipeline has higher precision.

 \section{Discussions}

@@ -635,7 +689,7 @@ chaining alone is more accurate than all the other long-read mappers in
 Fig.~\ref{fig:eval}a (data not shown). This accuracy helps to reduce downstream
 base-level alignment of candidate chains, which is still times slower than
 chaining even with the Suzuki-Kasahara improvement. In addition, taking a
-general form, minimap2 chaining can be adapted to non-typical data types such
+general form, minimap2 chaining can be adapted to non-typical data types such as
 spliced reads and multiple reads per fragment. This gives us the opportunity to
 extend the same base algorithm to a variety of use cases.

@@ -647,8 +701,9 @@ 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
 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.
+alleviates the concern with the seeding uniqueness. At the same time, at low
+sequence identity, it is rare to see long seeds anyway. Hash table is the ideal
+data structure for mapping long noisy sequences.

 \section*{Acknowledgements}
 We owe a debt of gratitude to H. Suzuki and M. Kasahara for releasing their
Author	SHA1	Message	Date
Heng Li	eeb314edd6	Release minimap2-2.9 (r720)	2018-02-24 09:31:09 -05:00
Heng Li	83c57a9d98	r719: fixed bad memory access	2018-02-23 17:27:41 -05:00
Heng Li	24a4808826	r718: retrieve sequence from the index	2018-02-23 10:18:26 -05:00
Heng Li	29ed675ee5	bugfix: in-place revcomp() not working	2018-02-20 11:05:54 -05:00
Heng Li	7dc7097208	added reverse complement	2018-02-20 09:41:25 -05:00
Heng Li	f434653432	added peakrss(); not used for now	2018-02-17 20:40:31 -05:00
Heng Li	090361c25b	fixed a typo in mappy (#117 )	2018-02-16 21:10:57 -05:00
Heng Li	e1f18690f6	paftools-r713: fixed bug for cov1 regions	2018-02-16 12:20:46 -05:00
Heng Li	54a42aafe5	keep documents in sync	2018-02-15 17:21:09 -05:00
Heng Li	8fc5f8dc90	r711: assign proper mapq to primary inversions	2018-02-15 14:34:59 -05:00
Heng Li	a0d62519c1	r710: fixed incorrect inversion coordinate (#112 )	2018-02-15 14:23:42 -05:00
Heng Li	b71c01b316	added test data for inversions	2018-02-15 11:04:27 -05:00
Heng Li	1372977a37	r708: implemented double Z-drop thresholds (#112 ) When aligning long reads, we would prefer to align through low-quality regions. This requires a large Z-drop threshold. However, to find small inversions, we need to use a small Z-drop. This commit address this conflict with two Z-drop thresholds. When Z-drop exceeds the smaller threshold, we perform a local alignment to check if there is a potential inversion. If there is one, we break the alignment; otherwise we break the alignment only if Z-drop excess the larger threshold. This commit also fixes a bug that reported wrong coordinates when the inversion is on the forward strand (#112).	2018-02-15 10:50:49 -05:00
Heng Li	c0e0d5d84b	r707: bugfix for inversions on rev strand (#112 )	2018-02-14 14:09:03 -05:00
Heng Li	b328795051	r706: don't segfault upon wrong FASTA/Q (#111 ) The lack of robustness cost me several hours to identify.	2018-02-13 10:00:22 -05:00
Heng Li	3b17d62ccd	minor doc improvements	2018-02-12 22:27:02 -05:00
Heng Li	874b8c4795	rework misc/README; progressing but unfinished	2018-02-12 22:06:05 -05:00
Heng Li	7ef5490884	r703: added --max-clip-ratio still testing the option	2018-02-12 13:29:18 -05:00
Heng Li	f66de7df59	updated paftools revision number to r702	2018-02-12 11:39:04 -05:00
Heng Li	fbbd4e0968	changed output messages for clarity	2018-02-12 11:32:32 -05:00
Heng Li	8c89ba005e	merged cnt-feat into paftools	2018-02-12 11:29:42 -05:00
Heng Li	50775a1e6f	added liftOver tested on a few toy examples	2018-02-12 11:12:23 -05:00
Heng Li	87cf650168	r697: fixed wrong #mapped reads in junceval	2018-02-09 18:45:57 -05:00
Heng Li	6e65c5e631	Revision 1	2018-02-09 17:39:02 -05:00
Heng Li	a58b05a61b	extend the section on genome alignment	2018-02-09 13:59:00 -05:00
Heng Li	42dab6319b	don't call SNPs involving 'n'	2018-02-09 13:27:17 -05:00
Heng Li	8428809369	convert MUMmer's delta to PAF	2018-02-09 12:14:28 -05:00
Heng Li	642af5591e	merged intron-eval into paftools	2018-02-09 11:07:14 -05:00
Heng Li	2025a6279a	merged gff2bed and mapstat to paftools	2018-02-09 10:47:12 -05:00
Heng Li	3465b04724	merged a few converters to paftools	2018-02-09 10:23:30 -05:00
Heng Li	66c5e71fa8	merging scripts into one long script	2018-02-09 10:00:35 -05:00
Heng Li	01560f1db0	Revision v1, to be submitted	2018-02-07 11:07:08 -05:00
Heng Li	39535565ee	first round of revision	2018-02-06 16:19:22 -05:00
Heng Li	a8d476c6ad	r686: end seed trimming don't go over long join	2018-02-06 11:31:32 -05:00
Heng Li	29b4a1786c	r685: tune end seed filter again	2018-02-05 11:48:22 -05:00
Heng Li	dbf284b2d9	r684: separate end score from min_chain_score	2018-02-05 11:40:38 -05:00
Heng Li	3df5015668	General doc improvements	2018-02-02 21:45:46 -05:00
Heng Li	756379bf83	Documented paf2diff.js (more to come later)	2018-02-02 14:43:58 -05:00
Heng Li	41fd8a966a	Documented ov-eval.js	2018-02-02 14:18:18 -05:00
Heng Li	86e4933b1a	Added TOC	2018-02-02 14:12:30 -05:00
Heng Li	a633a744b6	more doc in misc	2018-02-02 14:07:17 -05:00
Heng Li	ddc31f57ba	Documented some k8 scripts; more coming	2018-02-02 13:57:08 -05:00
Heng Li	35d3e064bf	r677: reduce the change of missing hits that are close to end of alignments. It is still possible to create examples that fail the heuristic.	2018-02-02 10:35:33 -05:00
Heng Li	997ab9bb2e	Fixed the description of --dual	2018-02-01 15:11:51 -05:00
Heng Li	53ce317e59	Release minimap2-2.8 (r672)	2018-02-01 12:50:20 -05:00
Heng Li	da6947cfa3	r671: cleanup command line options	2018-01-31 13:59:52 -05:00
Heng Li	46d6349af4	r670: added PE support to mappy and minor code cleanup	2018-01-31 11:33:08 -05:00
Heng Li	12a5a5fa3c	r669: improved self chain extension (#10 ) This has not fully resolved #10, only alleviated the issue.	2018-01-30 20:05:02 -05:00
Heng Li	ad18fa490d	fixed typos	2018-01-30 10:11:50 -05:00
Heng Li	43bfa6199d	r667: warn if one query file has fewer records #92	2018-01-28 17:36:21 -05:00
Heng Li	72b9b0e3b6	r666: report if >=3 query files in SR mode #92	2018-01-28 17:15:57 -05:00
Heng Li	6205fa6f21	document --heap-sort	2018-01-26 15:15:40 -05:00
Heng Li	d676a5314b	r664: use --heat-sort for sr by default	2018-01-26 12:25:42 -05:00
Heng Li	dfc78b39d3	refactor the old sorting	2018-01-26 09:37:48 -05:00
Heng Li	7b57c9a619	heap sort working on MT	2018-01-26 09:21:45 -05:00
Heng Li	123bc1d91d	put option operations in another file	2018-01-26 08:38:37 -05:00
Heng Li	dd18307e66	code backup	2018-01-25 21:52:49 -05:00
Heng Li	543fa12e68	r659: for C++ compatibility	2018-01-19 10:40:18 -05:00
Heng Li	af1a871270	r658: gives a warning if -N0 is used	2018-01-19 08:33:20 -05:00
Heng Li	2b71181a37	r657: check -p (#96 ) Well, in principle, every option should be checked. Will do when someone raise issues...	2018-01-19 01:03:38 -05:00
Heng Li	0454e6be91	explain -M in the manpage	2018-01-18 11:47:11 -05:00
Heng Li	33f8157961	r655: options to map to one strand of the ref #91	2018-01-16 10:34:30 -05:00