Released minimap2-2.7 (r654)

This happened when the last operation during backtracking is an intron.

.travis.yml

@@ -10,10 +10,6 @@ matrix:
       python: "2.7"
       before_install: pip install cython
       script: python setup.py build_ext
-    - language: python
-      python: "3.3"
-      before_install: pip install cython
-      script: python setup.py build_ext
     - language: python
       python: "3.5"
       before_install: pip install cython

Makefile

@@ -1,17 +1,24 @@
 CFLAGS=		-g -Wall -O2 -Wc++-compat
 CPPFLAGS=	-DHAVE_KALLOC
 INCLUDES=
-OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o index.o chain.o align.o hit.o map.o format.o pe.o ksw2_ll_sse.o
+OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o index.o chain.o align.o hit.o map.o format.o pe.o esterr.o ksw2_ll_sse.o
 PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite
 LIBS=		-lm -lz -lpthread
 
+ifeq ($(arm_neon),)
 ifeq ($(sse2only),)
 	OBJS+=ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o
 else
 	OBJS+=ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_exts2_sse.o
 endif
+else
+    OBJS+=ksw2_extz2_neon.o ksw2_extd2_neon.o ksw2_exts2_neon.o
+    CFLAGS+=-D_FILE_OFFSET_BITS=64 -mfpu=neon -fsigned-char
+    INCLUDES+=-I sse2neon
+endif
 
+.PHONY:all extra clean depend
 .SUFFIXES:.c .o
 
 .c.o:
@@ -33,6 +40,8 @@ libminimap2.a:$(OBJS)
 sdust:sdust.c getopt.o kalloc.o kalloc.h kdq.h kvec.h kseq.h sdust.h
 		$(CC) -D_SDUST_MAIN $(CFLAGS) $< getopt.o kalloc.o -o $@ -lz
 
+# SSE-specific targets on x86/x86_64
+
 ksw2_extz2_sse41.o:ksw2_extz2_sse.c ksw2.h kalloc.h
 		$(CC) -c -msse4 $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
 
@@ -54,6 +63,19 @@ ksw2_exts2_sse2.o:ksw2_exts2_sse.c ksw2.h kalloc.h
 ksw2_dispatch.o:ksw2_dispatch.c ksw2.h
 		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
 
+# NEON-specific targets on ARM
+
+ksw2_extz2_neon.o:ksw2_extz2_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_SSE2_ONLY -D__SSE2__ $(INCLUDES) $< -o $@
+
+ksw2_extd2_neon.o:ksw2_extd2_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_SSE2_ONLY -D__SSE2__ $(INCLUDES) $< -o $@
+
+ksw2_exts2_neon.o:ksw2_exts2_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_SSE2_ONLY -D__SSE2__ $(INCLUDES) $< -o $@
+
+# other non-file targets
+
 clean:
 		rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist mappy.so mappy.c python/mappy.c mappy.egg*
 
@@ -65,6 +87,7 @@ depend:
 align.o: minimap.h mmpriv.h bseq.h ksw2.h kalloc.h
 bseq.o: bseq.h kvec.h kalloc.h kseq.h
 chain.o: minimap.h mmpriv.h bseq.h kalloc.h
+esterr.o: mmpriv.h minimap.h bseq.h
 example.o: minimap.h kseq.h
 format.o: kalloc.h mmpriv.h minimap.h bseq.h
 getopt.o: getopt.h

NEWS.md

@@ -1,3 +1,60 @@
+Release 2.7-r654 (9 January 2018)
+---------------------------------
+
+This release fixed a bug in the splice mode and added a few minor features:
+
+ * Fixed a bug that occasionally takes an intron as a long deletion in the
+   splice mode. This was caused by wrong backtracking at the last CIGAR
+   operator. The current fix eliminates the error, but it is not optimal in
+   that it often produces a wrong junction when the last operator is an intron.
+   A future version of minimap2 may improve upon this.
+
+ * Support high-end ARM CPUs that implement the NEON instruction set (#81).
+   This enables minimap2 to work on Raspberry Pi 3 and Odroid XU4.
+
+ * Added a C API to construct a minimizer index from a set of C strings (#80).
+
+ * Check scoring specified on the command line (#79). Due to the 8-bit limit,
+   excessively large score penalties fail minimap2.
+
+For genomic sequences, minimap2 should give identical alignments to v2.6.
+
+(2.7: 9 January 2018, r654)
+
+
+
+Release 2.6-r623 (12 December 2017)
+-----------------------------------
+
+This release adds several features and fixes two minor bugs:
+
+ * Optionally build an index without sequences. This helps to reduce the
+   peak memory for read overlapping and is automatically applied when
+   base-level alignment is not requested.
+
+ * Approximately estimate per-base sequence divergence (i.e. 1-identity)
+   without performing base-level alignment, using a MashMap-like method. The
+   estimate is written to a new dv:f tag.
+
+ * Reduced the number of tiny terminal exons in RNA-seq alignment. The current
+   setting is conservative. Increase --end-seed-pen to drop more such exons.
+
+ * Reduced the peak memory when aligning long query sequences.
+
+ * Fixed a bug that is caused by HPC minimizers longer than 256bp. This should
+   have no effect in practice, but it is recommended to rebuild HPC indices if
+   possible.
+
+ * Fixed a bug when identifying identical hits (#71). This should only affect
+   artifactual reference consisting of near identical sequences.
+
+For genomic sequences, minimap2 should give nearly identical alignments to
+v2.5, except the new dv:f tag.
+
+(2.6: 12 December 2017, r623)
+
+
+
 Release 2.5-r572 (11 November 2017)
 -----------------------------------
 

README.md

@@ -1,10 +1,7 @@
-[![Release](https://img.shields.io/badge/Release-v2.5-blue.svg?style=flat)](https://github.com/lh3/minimap2/releases)
-[![BioConda](https://img.shields.io/conda/vn/bioconda/minimap2.svg?style=flat)](https://anaconda.org/bioconda/minimap2)
+[![GitHub Downloads](https://img.shields.io/github/downloads/lh3/minimap2/total.svg?style=social&logo=github&label=Download)](https://github.com/lh3/minimap2/releases)
+[![BioConda Install](https://img.shields.io/conda/dn/bioconda/minimap2.svg?style=flag&label=BioConda%20install)](https://anaconda.org/bioconda/minimap2)
 [![PyPI](https://img.shields.io/pypi/v/mappy.svg?style=flat)](https://pypi.python.org/pypi/mappy)
-[![Python Version](https://img.shields.io/pypi/pyversions/mappy.svg?style=flat)](https://pypi.python.org/pypi/mappy)
-[![License](https://img.shields.io/badge/License-MIT-blue.svg?style=flat)](LICENSE.txt)
 [![Build Status](https://travis-ci.org/lh3/minimap2.svg?branch=master)](https://travis-ci.org/lh3/minimap2)
-[![Downloads](https://img.shields.io/github/downloads/lh3/minimap2/total.svg?style=flat)](https://github.com/lh3/minimap2/releases)
 ## <a name="started"></a>Getting Started
 ```sh
 git clone https://github.com/lh3/minimap2
@@ -68,18 +65,21 @@ Detailed evaluations are available from the [minimap2 preprint][preprint].
 
 ### <a name="install"></a>Installation
 
-Minimap2 only works on x86-64 CPUs. You can acquire precompiled binaries from
+Minimap2 is optimized for x86-64 CPUs. You can acquire precompiled binaries from
 the [release page][release] with:
 ```sh
-curl -L https://github.com/lh3/minimap2/releases/download/v2.5/minimap2-2.5_x64-linux.tar.bz2 \
+curl -L https://github.com/lh3/minimap2/releases/download/v2.7/minimap2-2.7_x64-linux.tar.bz2 \
   | tar -jxvf -
-./minimap2-2.5_x64-linux/minimap2
+./minimap2-2.7_x64-linux/minimap2
 ```
 If you want to compile from the source, you need to have a C compiler, GNU make
 and zlib development files installed. Then type `make` in the source code
 directory to compile. If you see compilation errors, try `make sse2only=1`
 to disable SSE4 code, which will make minimap2 slightly slower.
 
+Minimap2 also works with ARM CPUs supporting the NEON instruction sets. To
+compile, use `make arm_neon=1`.
+
 ### <a name="general"></a>General usage
 
 Without any options, minimap2 takes a reference database and a query sequence
@@ -154,6 +154,20 @@ slight improvement to accuracy. For noisy Nanopore Direct RNA-seq reads, it is
 recommended to use a smaller k-mer size for increased sensitivity to the first
 or the last exons.
 
+Minimap2 rates an alignment by the score of the max-scoring sub-segment,
+*excluding* introns, and marks the best alignment as primary in SAM. When a
+spliced gene also has unspliced pseudogenes, minimap2 does not intentionally
+prefer spliced alignment, though in practice it more often marks the spliced
+alignment as the primary. By default, minimap2 outputs up to five secondary
+alignments (i.e. likely pseudogenes in the context of RNA-seq mapping). This
+can be tuned with option **-N**.
+
+For long RNA-seq reads, minimap2 may produce chimeric alignments potentially
+caused by gene fusions/structural variations or by an intron longer than the
+max intron length **-G** (200k by default). For now, it is not recommended to
+apply an excessively large **-G** as this slows down minimap2 and sometimes
+leads to false alignments.
+
 It is worth noting that by default `-x splice` prefers GT[A/G]..[C/T]AG
 over GT[C/T]..[A/G]AG, and then over other splicing signals. Considering
 one additional base improves the junction accuracy for noisy reads, but
@@ -309,7 +323,7 @@ highlighted in bold. The description may help to tune minimap2 parameters.
 
 ### <a name="help"></a>Getting help
 
-Manpage [minimap2.1](minimap2.1) provides detailed description of minimap2
+Manpage [minimap2.1][manpage] provides detailed description of minimap2
 command line options and optional tags. If you encounter bugs or have further
 questions or requests, you can raise an issue at the [issue page][issue].
 There is not a specific mailing list for the time being.
@@ -341,8 +355,9 @@ mappy` or [from BioConda][mappyconda] via `conda install -c bioconda mappy`.
   regions where seed positions may be suboptimal. This should not be a big
   concern because even the optimal alignment may be wrong in such regions.
 
-* Minimap2 requires SSE2 instructions to compile. It is possible to add
-  non-SSE2 support, but it would make minimap2 slower by several times.
+* Minimap2 requires SSE2 instructions on x86 CPUs or NEON on ARM CPUs. It is
+  possible to add non-SIMD support, but it would make minimap2 slower by
+  several times.
 
 In general, minimap2 is a young project with most code written since June, 2017.
 It may have bugs and room for improvements. Bug reports and suggestions are
@@ -363,3 +378,4 @@ warmly welcomed.
 [mappyconda]: https://anaconda.org/bioconda/mappy
 [issue]: https://github.com/lh3/minimap2/issues
 [k8]: https://github.com/attractivechaos/k8
+[manpage]: https://lh3.github.io/minimap2/minimap2.html

align.c

@@ -1,6 +1,7 @@
 #include <assert.h>
 #include <string.h>
 #include <stdlib.h>
+#include <math.h>
 #include "minimap.h"
 #include "mmpriv.h"
 #include "ksw2.h"
@@ -120,7 +121,7 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 	}
 }
 
-static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *qual, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e)
+static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e)
 {
 	uint32_t k, l, toff = 0, qoff = 0;
 	int32_t s = 0, max = 0, qshift, tshift;
@@ -209,6 +210,13 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 		ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, flag, ez);
 	else
 		ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez);
+	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
+		int i;
+		fprintf(stderr, "score=%d, cigar=", ez->score);
+		for (i = 0; i < ez->n_cigar; ++i)
+			fprintf(stderr, "%d%c", ez->cigar[i]>>4, "MIDN"[ez->cigar[i]&0xf]);
+		fprintf(stderr, "\n");
+	}
 }
 
 static inline int mm_get_hplen_back(const mm_idx_t *mi, uint32_t rid, uint32_t x)
@@ -222,7 +230,7 @@ static inline int mm_get_hplen_back(const mm_idx_t *mi, uint32_t rid, uint32_t x
 
 static inline void mm_adjust_minier(const mm_idx_t *mi, uint8_t *const qseq0[2], mm128_t *a, int32_t *r, int32_t *q)
 {
-	if (mi->is_hpc) {
+	if (mi->flag & MM_I_HPC) {
 		const uint8_t *qseq = qseq0[a->x>>63];
 		int i, c;
 		*q = (int32_t)a->y;
@@ -341,7 +349,49 @@ static void mm_max_stretch(const mm_mapopt_t *opt, const mm_reg1_t *r, const mm1
 	*as = max_i, *cnt = max_len;
 }
 
-static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], uint8_t *qual0[2], mm_reg1_t *r, mm_reg1_t *r2, int n_a, mm128_t *a, ksw_extz_t *ez, int splice_flag)
+static int mm_seed_ext_score(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, const int8_t mat[25], int qlen, uint8_t *qseq0[2], const mm128_t *a)
+{
+	uint8_t *qseq, *tseq;
+	int q_span = a->y>>32&0xff, qs, qe, rs, re, rid, score, q_off, t_off, ext_len = opt->anchor_ext_len;
+	void *qp;
+	rid = a->x<<1>>33;
+	re = (uint32_t)a->x + 1, rs = re - q_span;
+	qe = (uint32_t)a->y + 1, qs = qe - q_span;
+	rs = rs - ext_len > 0? rs - ext_len : 0;
+	qs = qs - ext_len > 0? qs - ext_len : 0;
+	re = re + ext_len < mi->seq[rid].len? re + ext_len : mi->seq[rid].len;
+	qe = qe + ext_len < qlen? qe + ext_len : qlen;
+	tseq = (uint8_t*)kmalloc(km, re - rs);
+	mm_idx_getseq(mi, rid, rs, re, tseq);
+	qseq = qseq0[a->x>>63] + qs;
+	qp = ksw_ll_qinit(km, 2, qe - qs, qseq, 5, mat);
+	score = ksw_ll_i16(qp, re - rs, tseq, opt->q, opt->e, &q_off, &t_off);
+	kfree(km, tseq);
+	kfree(km, qp);
+	return score;
+}
+
+static void mm_fix_bad_ends_splice(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, const mm_reg1_t *r, const int8_t mat[25], int qlen, uint8_t *qseq0[2], const mm128_t *a, int *as1, int *cnt1)
+{ // this assumes a very crude k-mer based mode; it is not necessary to use a good model just for filtering bounary exons
+	int score;
+	double log_gap;
+	*as1 = r->as, *cnt1 = r->cnt;
+	if (r->cnt < 3) return;
+	log_gap = log((int32_t)a[r->as + 1].x - (int32_t)a[r->as].x);
+	if ((a[r->as].y>>32&0xff) < log_gap + opt->anchor_ext_shift) {
+		score = mm_seed_ext_score(km, opt, mi, mat, qlen, qseq0, &a[r->as]);
+		if ((double)score / mat[0] < log_gap + opt->anchor_ext_shift) // a more exact format is "score < log_4(gap) + shift"
+			++(*as1), --(*cnt1);
+	}
+	log_gap = log((int32_t)a[r->as + r->cnt - 1].x - (int32_t)a[r->as + r->cnt - 2].x);
+	if ((a[r->as + r->cnt - 1].y>>32&0xff) < log_gap + opt->anchor_ext_shift) {
+		score = mm_seed_ext_score(km, opt, mi, mat, qlen, qseq0, &a[r->as + r->cnt - 1]);
+		if ((double)score / mat[0] < log_gap + opt->anchor_ext_shift)
+			--(*cnt1);
+	}
+}
+
+static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], mm_reg1_t *r, mm_reg1_t *r2, int n_a, mm128_t *a, ksw_extz_t *ez, int splice_flag)
 {
 	int is_sr = !!(opt->flag & MM_F_SR), is_splice = !!(opt->flag & MM_F_SPLICE);
 	int32_t rid = a[r->as].x<<1>>33, rev = a[r->as].x>>63, as1, cnt1;
@@ -351,23 +401,25 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	int32_t rs1, qs1, re1, qe1;
 	int8_t mat[25];
 
-	if (is_sr) assert(!mi->is_hpc); // HPC won't work with SR because with HPC we can't easily tell if there is a gap
+	if (is_sr) assert(!(mi->flag & MM_I_HPC)); // HPC won't work with SR because with HPC we can't easily tell if there is a gap
 
 	r2->cnt = 0;
 	if (r->cnt == 0) return;
 	ksw_gen_simple_mat(5, mat, opt->a, opt->b);
 	bw = (int)(opt->bw * 1.5 + 1.);
 
-	if (is_sr && !mi->is_hpc) {
+	if (is_sr && !(mi->flag & MM_I_HPC)) {
 		mm_max_stretch(opt, r, a, &as1, &cnt1);
 		rs = (int32_t)a[as1].x + 1 - (int32_t)(a[as1].y>>32&0xff);
 		qs = (int32_t)a[as1].y + 1 - (int32_t)(a[as1].y>>32&0xff);
 		re = (int32_t)a[as1+cnt1-1].x + 1;
 		qe = (int32_t)a[as1+cnt1-1].y + 1;
 	} else {
-		if (!is_splice)
+		if (is_splice) {
+			mm_fix_bad_ends_splice(km, opt, mi, r, mat, qlen, qseq0, a, &as1, &cnt1);
+		} else {
 			mm_fix_bad_ends(r, a, opt->bw, &as1, &cnt1);
-		else as1 = r->as, cnt1 = r->cnt;
+		}
 		mm_filter_bad_seeds(km, as1, cnt1, a, 10, 40, opt->max_gap>>1, 10);
 		mm_adjust_minier(mi, qseq0, &a[as1], &rs, &qs);
 		mm_adjust_minier(mi, qseq0, &a[as1 + cnt1 - 1], &re, &qe);
@@ -470,7 +522,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 
 	for (i = is_sr? cnt1 - 1 : 1; i < cnt1; ++i) { // gap filling
 		if ((a[as1+i].y & (MM_SEED_IGNORE|MM_SEED_TANDEM)) && i != cnt1 - 1) continue;
-		if (is_sr && !mi->is_hpc) {
+		if (is_sr && !(mi->flag & MM_I_HPC)) {
 			re = (int32_t)a[as1 + i].x + 1;
 			qe = (int32_t)a[as1 + i].y + 1;
 		} else mm_adjust_minier(mi, qseq0, &a[as1 + i], &re, &qe);
@@ -533,7 +585,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	assert(re1 - rs1 <= re0 - rs0);
 	if (r->p) {
 		mm_idx_getseq(mi, rid, rs1, re1, tseq);
-		mm_update_extra(r, &qseq0[r->rev][qs1], qual0[r->rev]? &qual0[r->rev][qs1] : 0, tseq, mat, opt->q, opt->e);
+		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e);
 		if (rev && r->p->trans_strand)
 			r->p->trans_strand ^= 3; // flip to the read strand
 	}
@@ -541,10 +593,10 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	kfree(km, tseq);
 }
 
-static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], uint8_t *qual0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez)
+static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez)
 {
 	int tl, ql, score, ret = 0, q_off, t_off;
-	uint8_t *tseq, *qseq, *qual;
+	uint8_t *tseq, *qseq;
 	int8_t mat[25];
 	void *qp;
 
@@ -562,7 +614,6 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	tseq = (uint8_t*)kmalloc(km, tl);
 	mm_idx_getseq(mi, r1->rid, r1->re, r2->rs, tseq);
 	qseq = &qseq0[!r1->rev][qlen - r2->qs];
-	qual = qual0[!r1->rev]? &qseq0[!r1->rev][qlen - r2->qs] : 0;
 
 	mm_seq_rev(ql, qseq);
 	mm_seq_rev(tl, tseq);
@@ -582,9 +633,10 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	r_inv->inv = 1;
 	r_inv->rev = !r1->rev;
 	r_inv->rid = r1->rid;
+	r_inv->div = -1.0f;
 	r_inv->qs = r1->qe + q_off, r_inv->qe = r_inv->qs + ez->max_q + 1;
 	r_inv->rs = r1->re + t_off, r_inv->re = r_inv->rs + ez->max_t + 1;
-	mm_update_extra(r_inv, &qseq[q_off], qual? &qual[q_off] : 0, &tseq[t_off], mat, opt->q, opt->e);
+	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e);
 	ret = 1;
 end_align1_inv:
 	kfree(km, tseq);
@@ -601,11 +653,11 @@ static inline mm_reg1_t *mm_insert_reg(const mm_reg1_t *r, int i, int *n_regs, m
 	return regs;
 }
 
-mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, const char *qual, int *n_regs_, mm_reg1_t *regs, mm128_t *a)
+mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a)
 {
 	extern unsigned char seq_nt4_table[256];
 	int32_t i, n_regs = *n_regs_, n_a;
-	uint8_t *qseq0[2], *qual0[2];
+	uint8_t *qseq0[2];
 	ksw_extz_t ez;
 
 	// encode the query sequence
@@ -615,12 +667,6 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 		qseq0[0][i] = seq_nt4_table[(uint8_t)qstr[i]];
 		qseq0[1][qlen - 1 - i] = qseq0[0][i] < 4? 3 - qseq0[0][i] : 4;
 	}
-	if (qual) {
-		qual0[0] = (uint8_t*)kmalloc(km, qlen * 2);
-		qual0[1] = qual0[0] + qlen;
-		for (i = 0; i < qlen; ++i)
-			qual0[0][i] = qual0[1][qlen - 1 - i] = qual[i] - 33;
-	} else qual0[0] = qual0[1] = 0;
 
 	// align through seed hits
 	n_a = mm_squeeze_a(km, n_regs, regs, a);
@@ -631,8 +677,8 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 			mm_reg1_t s[2], s2[2];
 			int which, trans_strand;
 			s[0] = s[1] = regs[i];
-			mm_align1(km, opt, mi, qlen, qseq0, qual0, &s[0], &s2[0], n_a, a, &ez, MM_F_SPLICE_FOR);
-			mm_align1(km, opt, mi, qlen, qseq0, qual0, &s[1], &s2[1], n_a, a, &ez, MM_F_SPLICE_REV);
+			mm_align1(km, opt, mi, qlen, qseq0, &s[0], &s2[0], n_a, a, &ez, MM_F_SPLICE_FOR);
+			mm_align1(km, opt, mi, qlen, qseq0, &s[1], &s2[1], n_a, a, &ez, MM_F_SPLICE_REV);
 			if (s[0].p->dp_score > s[1].p->dp_score) which = 0, trans_strand = 1;
 			else if (s[0].p->dp_score < s[1].p->dp_score) which = 1, trans_strand = 2;
 			else trans_strand = 3, which = (qlen + s[0].p->dp_score) & 1; // randomly choose a strand, effectively
@@ -645,13 +691,13 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 			}
 			regs[i].p->trans_strand = trans_strand;
 		} else { // one round of alignment
-			mm_align1(km, opt, mi, qlen, qseq0, qual0, &regs[i], &r2, n_a, a, &ez, opt->flag);
+			mm_align1(km, opt, mi, qlen, qseq0, &regs[i], &r2, n_a, a, &ez, opt->flag);
 			if (opt->flag&MM_F_SPLICE)
 				regs[i].p->trans_strand = opt->flag&MM_F_SPLICE_FOR? 1 : 2;
 		}
 		if (r2.cnt > 0) regs = mm_insert_reg(&r2, i, &n_regs, regs);
 		if (!(opt->flag&MM_F_SPLICE) && !(opt->flag&MM_F_SR) && i > 0) { // don't try inversion alignment for -xsplice or -xsr
-			if (mm_align1_inv(km, opt, mi, qlen, qseq0, qual0, &regs[i-1], &regs[i], &r2, &ez)) {
+			if (mm_align1_inv(km, opt, mi, qlen, qseq0, &regs[i-1], &regs[i], &r2, &ez)) {
 				regs = mm_insert_reg(&r2, i, &n_regs, regs);
 				++i; // skip the inserted INV alignment
 			}
@@ -659,7 +705,6 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 	}
 	*n_regs_ = n_regs;
 	kfree(km, qseq0[0]);
-	if (qual0[0]) kfree(km, qual0[0]);
 	kfree(km, ez.cigar);
 	mm_filter_regs(km, opt, n_regs_, regs);
 	mm_hit_sort_by_dp(km, n_regs_, regs);

esterr.c

@@ -0,0 +1,64 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+#include "mmpriv.h"
+
+static inline int32_t get_for_qpos(int32_t qlen, const mm128_t *a)
+{
+	int32_t x = (int32_t)a->y;
+	int32_t q_span = a->y>>32 & 0xff;
+	if (a->x>>63)
+		x = qlen - 1 - (x + 1 - q_span); // revert the position to the forward strand of query
+	return x;
+}
+
+static int get_mini_idx(int qlen, const mm128_t *a, int32_t n, const uint64_t *mini_pos)
+{
+	int32_t x, L = 0, R = n - 1;
+	x = get_for_qpos(qlen, a);
+	while (L <= R) { // binary search
+		int32_t m = ((uint64_t)L + R) >> 1;
+		int32_t y = (int32_t)mini_pos[m];
+		if (y < x) L = m + 1;
+		else if (y > x) R = m - 1;
+		else return m;
+	}
+	return -1;
+}
+
+void mm_est_err(const mm_idx_t *mi, int qlen, int n_regs, mm_reg1_t *regs, const mm128_t *a, int32_t n, const uint64_t *mini_pos)
+{
+	int i;
+	uint64_t sum_k = 0;
+	float avg_k;
+
+	if (n == 0) return;
+	for (i = 0; i < n; ++i)
+		sum_k += mini_pos[i] >> 32 & 0xff;
+	avg_k = (float)sum_k / n;
+
+	for (i = 0; i < n_regs; ++i) {
+		mm_reg1_t *r = &regs[i];
+		int32_t st, en, j, k, n_match, n_tot, l_ref;
+		r->div = -1.0f;
+		if (r->cnt == 0) continue;
+		st = en = get_mini_idx(qlen, r->rev? &a[r->as + r->cnt - 1] : &a[r->as], n, mini_pos);
+		if (st < 0) {
+			if (mm_verbose >= 2)
+				fprintf(stderr, "[WARNING] logic inconsistency in mm_est_err(). Please contact the developer.\n");
+			continue;
+		}
+		l_ref = mi->seq[r->rid].len;
+		for (k = 1, j = st + 1, n_match = 1; j < n && k < r->cnt; ++j) {
+			int32_t x;
+			x = get_for_qpos(qlen, r->rev? &a[r->as + r->cnt - 1 - k] : &a[r->as + k]);
+			if (x == (int32_t)mini_pos[j])
+				++k, en = j, ++n_match;
+		}
+		n_tot = en - st + 1;
+		if (r->qs > avg_k && r->rs > avg_k) ++n_tot;
+		if (qlen - r->qs > avg_k && l_ref - r->re > avg_k) ++n_tot;
+		r->div = logf((float)n_tot / n_match) / avg_k;
+	}
+}

format.c

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

hit.c

@@ -80,6 +80,7 @@ mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u,
 		ri->hash = (uint32_t)z[i].x;
 		ri->cnt = (int32_t)z[i].y;
 		ri->as = z[i].y >> 32;
+		ri->div = -1.0f;
 		mm_reg_set_coor(ri, qlen, a);
 	}
 	kfree(km, z);
@@ -147,7 +148,7 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff
 				ri->parent = rp->parent;
 				rp->subsc = rp->subsc > ri->score? rp->subsc : ri->score;
 				if (ri->cnt >= rp->cnt) cnt_sub = 1;
-				if (rp->p && ri->p && (rp->rs != ri->rs || rp->re != ri->re || ol != min)) { // the last condition excludes identical hits after DP
+				if (rp->p && ri->p && (rp->rid != ri->rid || rp->rs != ri->rs || rp->re != ri->re || ol != min)) { // the last condition excludes identical hits after DP
 					rp->p->dp_max2 = rp->p->dp_max2 > ri->p->dp_max? rp->p->dp_max2 : ri->p->dp_max;
 					if (rp->p->dp_max - ri->p->dp_max <= sub_diff) cnt_sub = 1;
 				}
@@ -234,7 +235,7 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_,
 			if (p == i || r[i].inv) { // primary or inversion
 				r[k++] = r[i];
 			} else if ((r[i].score >= r[p].score * pri_ratio || r[i].score + min_diff >= r[p].score) && n_2nd < best_n) {
-				if (!(r[i].qs == r[p].qs && r[i].qe == r[p].qe && r[i].rs == r[p].rs && r[i].re == r[p].re)) // not identical hits
+				if (!(r[i].qs == r[p].qs && r[i].qe == r[p].qe && r[i].rid == r[p].rid && r[i].rs == r[p].rs && r[i].re == r[p].re)) // not identical hits
 					r[k++] = r[i], ++n_2nd;
 				else if (r[i].p) free(r[i].p);
 			} else if (r[i].p) free(r[i].p);

index.c

@@ -31,19 +31,19 @@ typedef struct mm_idx_bucket_s {
 void mm_idxopt_init(mm_idxopt_t *opt)
 {
 	memset(opt, 0, sizeof(mm_idxopt_t));
-	opt->k = 15, opt->w = 10, opt->is_hpc = 0;
+	opt->k = 15, opt->w = 10, opt->flag = 0;
 	opt->bucket_bits = 14;
 	opt->mini_batch_size = 50000000;
 	opt->batch_size = 4000000000ULL;
 }
 
-mm_idx_t *mm_idx_init(int w, int k, int b, int is_hpc)
+mm_idx_t *mm_idx_init(int w, int k, int b, int flag)
 {
 	mm_idx_t *mi;
 	if (k*2 < b) b = k * 2;
 	if (w < 1) w = 1;
 	mi = (mm_idx_t*)calloc(1, sizeof(mm_idx_t));
-	mi->w = w, mi->k = k, mi->b = b, mi->is_hpc = is_hpc;
+	mi->w = w, mi->k = k, mi->b = b, mi->flag = flag;
 	mi->B = (mm_idx_bucket_t*)calloc(1<<b, sizeof(mm_idx_bucket_t));
 	if (!(mm_dbg_flag & 1)) mi->km = km_init();
 	return mi;
@@ -89,7 +89,7 @@ void mm_idx_stat(const mm_idx_t *mi)
 {
 	int i, n = 0, n1 = 0;
 	uint64_t sum = 0, len = 0;
-	fprintf(stderr, "[M::%s] kmer size: %d; skip: %d; is_HPC: %d; #seq: %d\n", __func__, mi->k, mi->w, mi->is_hpc, mi->n_seq);
+	fprintf(stderr, "[M::%s] kmer size: %d; skip: %d; is_hpc: %d; #seq: %d\n", __func__, mi->k, mi->w, mi->flag&MM_I_HPC, mi->n_seq);
 	for (i = 0; i < mi->n_seq; ++i)
 		len += mi->seq[i].len;
 	for (i = 0; i < 1<<mi->b; ++i)
@@ -214,7 +214,7 @@ static void mm_idx_post(mm_idx_t *mi, int n_threads)
 #include "bseq.h"
 
 typedef struct {
-	int mini_batch_size, keep_name;
+	int mini_batch_size;
 	uint64_t batch_size, sum_len;
 	mm_bseq_file_t *fp;
 	mm_idx_t *mi;
@@ -246,7 +246,6 @@ static void *worker_pipeline(void *shared, int step, void *in)
 		s->seq = mm_bseq_read(p->fp, p->mini_batch_size, 0, &s->n_seq); // read a mini-batch
 		if (s->seq) {
 			uint32_t old_m, m;
-			uint64_t sum_len, old_max_len, max_len;
 			assert((uint64_t)p->mi->n_seq + s->n_seq <= UINT32_MAX); // to prevent integer overflow
 			// make room for p->mi->seq
 			old_m = p->mi->n_seq, m = p->mi->n_seq + s->n_seq;
@@ -254,29 +253,34 @@ static void *worker_pipeline(void *shared, int step, void *in)
 			if (old_m != m)
 				p->mi->seq = (mm_idx_seq_t*)krealloc(p->mi->km, p->mi->seq, m * sizeof(mm_idx_seq_t));
 			// make room for p->mi->S
-			for (i = 0, sum_len = 0; i < s->n_seq; ++i) sum_len += s->seq[i].l_seq;
-			old_max_len = (p->sum_len + 7) / 8;
-			max_len = (p->sum_len + sum_len + 7) / 8;
-			kroundup64(old_max_len); kroundup64(max_len);
-			if (old_max_len != max_len) {
-				p->mi->S = (uint32_t*)realloc(p->mi->S, max_len * 4);
-				memset(&p->mi->S[old_max_len], 0, 4 * (max_len - old_max_len));
+			if (!(p->mi->flag & MM_I_NO_SEQ)) {
+				uint64_t sum_len, old_max_len, max_len;
+				for (i = 0, sum_len = 0; i < s->n_seq; ++i) sum_len += s->seq[i].l_seq;
+				old_max_len = (p->sum_len + 7) / 8;
+				max_len = (p->sum_len + sum_len + 7) / 8;
+				kroundup64(old_max_len); kroundup64(max_len);
+				if (old_max_len != max_len) {
+					p->mi->S = (uint32_t*)realloc(p->mi->S, max_len * 4);
+					memset(&p->mi->S[old_max_len], 0, 4 * (max_len - old_max_len));
+				}
 			}
 			// populate p->mi->seq
 			for (i = 0; i < s->n_seq; ++i) {
 				mm_idx_seq_t *seq = &p->mi->seq[p->mi->n_seq];
 				uint32_t j;
-				if (p->keep_name) {
+				if (!(p->mi->flag & MM_I_NO_NAME)) {
 					seq->name = (char*)kmalloc(p->mi->km, strlen(s->seq[i].name) + 1);
 					strcpy(seq->name, s->seq[i].name);
 				} else seq->name = 0;
 				seq->len = s->seq[i].l_seq;
 				seq->offset = p->sum_len;
 				// copy the sequence
-				for (j = 0; j < seq->len; ++j) { // TODO: this is not the fastest way, but let's first see if speed matters here
-					uint64_t o = p->sum_len + j;
-					int c = seq_nt4_table[(uint8_t)s->seq[i].seq[j]];
-					mm_seq4_set(p->mi->S, o, c);
+				if (!(p->mi->flag & MM_I_NO_SEQ)) {
+					for (j = 0; j < seq->len; ++j) { // TODO: this is not the fastest way, but let's first see if speed matters here
+						uint64_t o = p->sum_len + j;
+						int c = seq_nt4_table[(uint8_t)s->seq[i].seq[j]];
+						mm_seq4_set(p->mi->S, o, c);
+					}
 				}
 				// update p->sum_len and p->mi->n_seq
 				p->sum_len += seq->len;
@@ -288,7 +292,10 @@ static void *worker_pipeline(void *shared, int step, void *in)
         step_t *s = (step_t*)in;
 		for (i = 0; i < s->n_seq; ++i) {
 			mm_bseq1_t *t = &s->seq[i];
-			mm_sketch(0, t->seq, t->l_seq, p->mi->w, p->mi->k, t->rid, p->mi->is_hpc, &s->a);
+			if (t->l_seq > 0)
+				mm_sketch(0, t->seq, t->l_seq, p->mi->w, p->mi->k, t->rid, p->mi->flag&MM_I_HPC, &s->a);
+			else if (mm_verbose >= 2)
+				fprintf(stderr, "[WARNING] the length database sequence '%s' is 0\n", t->name);
 			free(t->seq); free(t->name);
 		}
 		free(s->seq); s->seq = 0;
@@ -301,16 +308,15 @@ static void *worker_pipeline(void *shared, int step, void *in)
     return 0;
 }
 
-mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int is_hpc, int mini_batch_size, int n_threads, uint64_t batch_size, int keep_name)
+mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int flag, int mini_batch_size, int n_threads, uint64_t batch_size)
 {
 	pipeline_t pl;
 	if (fp == 0 || mm_bseq_eof(fp)) return 0;
 	memset(&pl, 0, sizeof(pipeline_t));
 	pl.mini_batch_size = mini_batch_size < batch_size? mini_batch_size : batch_size;
-	pl.keep_name = keep_name;
 	pl.batch_size = batch_size;
 	pl.fp = fp;
-	pl.mi = mm_idx_init(w, k, b, is_hpc);
+	pl.mi = mm_idx_init(w, k, b, flag);
 
 	kt_pipeline(n_threads < 3? n_threads : 3, worker_pipeline, &pl, 3);
 	if (mm_verbose >= 3)
@@ -323,17 +329,60 @@ mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int is_hpc, int mi
 	return pl.mi;
 }
 
-mm_idx_t *mm_idx_build(const char *fn, int w, int k, int is_hpc, int n_threads) // a simpler interface
+mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads) // a simpler interface; deprecated
 {
 	mm_bseq_file_t *fp;
 	mm_idx_t *mi;
 	fp = mm_bseq_open(fn);
 	if (fp == 0) return 0;
-	mi = mm_idx_gen(fp, w, k, 14, is_hpc, 1<<18, n_threads, UINT64_MAX, 1);
+	mi = mm_idx_gen(fp, w, k, 14, flag, 1<<18, n_threads, UINT64_MAX);
 	mm_bseq_close(fp);
 	return mi;
 }
 
+mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const char **seq, const char **name)
+{
+	uint64_t sum_len = 0;
+	mm128_v a = {0,0,0};
+	mm_idx_t *mi;
+	int i, flag = 0;
+	if (n <= 0) return 0;
+	for (i = 0; i < n; ++i) // get the total length
+		sum_len += strlen(seq[i]);
+	if (is_hpc) flag |= MM_I_HPC;
+	if (name == 0) flag |= MM_I_NO_NAME;
+	if (bucket_bits < 0) bucket_bits = 14;
+	mi = mm_idx_init(w, k, bucket_bits, flag);
+	mi->n_seq = n;
+	mi->seq = (mm_idx_seq_t*)kcalloc(mi->km, n, sizeof(mm_idx_seq_t)); // ->seq is allocated from km
+	mi->S = (uint32_t*)calloc((sum_len + 7) / 8, 4);
+	for (i = 0, sum_len = 0; i < n; ++i) {
+		const char *s = seq[i];
+		mm_idx_seq_t *p = &mi->seq[i];
+		uint32_t j;
+		if (name && name[i]) {
+			p->name = (char*)kmalloc(mi->km, strlen(name[i]) + 1);
+			strcpy(p->name, name[i]);
+		}
+		p->offset = sum_len;
+		p->len = strlen(s);
+		for (j = 0; j < p->len; ++j) {
+			int c = seq_nt4_table[(uint8_t)s[j]];
+			uint64_t o = sum_len + j;
+			mm_seq4_set(mi->S, o, c);
+		}
+		sum_len += p->len;
+		if (p->len > 0) {
+			a.n = 0;
+			mm_sketch(0, s, p->len, w, k, i, is_hpc, &a);
+			mm_idx_add(mi, a.n, a.a);
+		}
+	}
+	free(a.a);
+	mm_idx_post(mi, 1);
+	return mi;
+}
+
 /*************
  * index I/O *
  *************/
@@ -344,7 +393,7 @@ void mm_idx_dump(FILE *fp, const mm_idx_t *mi)
 	uint32_t x[5];
 	int i;
 
-	x[0] = mi->w, x[1] = mi->k, x[2] = mi->b, x[3] = mi->n_seq, x[4] = mi->is_hpc;
+	x[0] = mi->w, x[1] = mi->k, x[2] = mi->b, x[3] = mi->n_seq, x[4] = mi->flag;
 	fwrite(MM_IDX_MAGIC, 1, 4, fp);
 	fwrite(x, 4, 5, fp);
 	for (i = 0; i < mi->n_seq; ++i) {
@@ -371,7 +420,8 @@ void mm_idx_dump(FILE *fp, const mm_idx_t *mi)
 			fwrite(x, 8, 2, fp);
 		}
 	}
-	fwrite(mi->S, 4, (sum_len + 7) / 8, fp);
+	if (!(mi->flag & MM_I_NO_SEQ))
+		fwrite(mi->S, 4, (sum_len + 7) / 8, fp);
 	fflush(fp);
 }
 
@@ -421,8 +471,10 @@ mm_idx_t *mm_idx_load(FILE *fp)
 			kh_val(h, k) = x[1];
 		}
 	}
-	mi->S = (uint32_t*)malloc((sum_len + 7) / 8 * 4);
-	fread(mi->S, 4, (sum_len + 7) / 8, fp);
+	if (!(mi->flag & MM_I_NO_SEQ)) {
+		mi->S = (uint32_t*)malloc((sum_len + 7) / 8 * 4);
+		fread(mi->S, 4, (sum_len + 7) / 8, fp);
+	}
 	return mi;
 }
 
@@ -476,10 +528,10 @@ mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads)
 	mm_idx_t *mi;
 	if (r->is_idx) {
 		mi = mm_idx_load(r->fp.idx);
-		if (mi && mm_verbose >= 2 && (mi->k != r->opt.k || mi->w != r->opt.w || mi->is_hpc != r->opt.is_hpc))
+		if (mi && mm_verbose >= 2 && (mi->k != r->opt.k || mi->w != r->opt.w || (mi->flag&MM_I_HPC) != (r->opt.flag&MM_I_HPC)))
 			fprintf(stderr, "[WARNING]\033[1;31m Indexing parameters (-k, -w or -H) overridden by parameters used in the prebuilt index.\033[0m\n");
 	} else
-		mi = mm_idx_gen(r->fp.seq, r->opt.w, r->opt.k, r->opt.bucket_bits, r->opt.is_hpc, r->opt.mini_batch_size, n_threads, r->opt.batch_size, 1);
+		mi = mm_idx_gen(r->fp.seq, r->opt.w, r->opt.k, r->opt.bucket_bits, r->opt.flag, r->opt.mini_batch_size, n_threads, r->opt.batch_size);
 	if (mi) {
 		if (r->fp_out) mm_idx_dump(r->fp_out, mi);
 		++r->n_parts;

kalloc.c

@@ -189,6 +189,10 @@ void km_stat(const void *_km, km_stat_t *s)
 			panic("[km_stat] The end of a free block enters another free block.");
 		if (p->ptr == km->loop_head) break;
 	}
-	for (p = km->core_head; p != NULL; p = p->ptr)
-		++s->n_cores, s->capacity += p->size * sizeof(header_t);
+	for (p = km->core_head; p != NULL; p = p->ptr) {
+		size_t size = p->size * sizeof(header_t);
+		++s->n_cores;
+		s->capacity += size;
+		s->largest = s->largest > size? s->largest : size;
+	}
 }

kalloc.h

@@ -8,7 +8,7 @@ extern "C" {
 #endif
 
 typedef struct {
-	size_t capacity, available, n_blocks, n_cores;
+	size_t capacity, available, n_blocks, n_cores, largest;
 } km_stat_t;
 
 void *kmalloc(void *km, size_t size);

ksw2.h

@@ -116,7 +116,7 @@ static inline uint32_t *ksw_push_cigar(void *km, int *n_cigar, int *m_cigar, uin
 //   bit 0-2: which type gets the max - 0 for H, 1 for E, 2 for F, 3 for \tilde{E} and 4 for \tilde{F}
 //   bit 3/0x08: 1 if a continuation on the E state (bit 5/0x20 for a continuation on \tilde{E})
 //   bit 4/0x10: 1 if a continuation on the F state (bit 6/0x40 for a continuation on \tilde{F})
-static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int with_N, const uint8_t *p, const int *off, const int *off_end, int n_col, int i0, int j0,
+static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int min_intron_len, const uint8_t *p, const int *off, const int *off_end, int n_col, int i0, int j0,
 								 int *m_cigar_, int *n_cigar_, uint32_t **cigar_)
 { // p[] - lower 3 bits: which type gets the max; bit
 	int n_cigar = 0, m_cigar = *m_cigar_, i = i0, j = j0, r, state = 0;
@@ -138,11 +138,11 @@ static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int with_N, c
 		if (state == 0) state = tmp & 7; // TODO: probably this line can be merged into the "else if" line right above; not 100% sure
 		if (force_state >= 0) state = force_state;
 		if (state == 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 0, 1), --i, --j; // match
-		else if (state == 1 || (state == 3 && !with_N)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i; // deletion
-		else if (state == 3 && with_N) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 3, 1), --i; // intron
+		else if (state == 1 || (state == 3 && min_intron_len <= 0)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i; // deletion
+		else if (state == 3 && min_intron_len > 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 3, 1), --i; // intron
 		else cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, 1), --j; // insertion
 	}
-	if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, i + 1); // first deletion
+	if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, min_intron_len > 0 && i >= min_intron_len? 3 : 2, i + 1); // first deletion
 	if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, j + 1); // first insertion
 	if (!is_rev)
 		for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR

ksw2_exts2_sse.c

@@ -367,9 +367,9 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	if (with_cigar) { // backtrack
 		int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR);
 		if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY))
-			ksw_backtrack(km, 1, rev_cigar, 1, (uint8_t*)p, off, off_end, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+			ksw_backtrack(km, 1, rev_cigar, long_thres, (uint8_t*)p, off, off_end, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
 		else if (ez->max_t >= 0 && ez->max_q >= 0)
-			ksw_backtrack(km, 1, rev_cigar, 1, (uint8_t*)p, off, off_end, n_col_*16, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
+			ksw_backtrack(km, 1, rev_cigar, long_thres, (uint8_t*)p, off, off_end, n_col_*16, ez->max_t, ez->max_q, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
 		kfree(km, mem2); kfree(km, off);
 	}
 }

main.c

@@ -6,7 +6,7 @@
 #include "mmpriv.h"
 #include "getopt.h"
 
-#define MM_VERSION "2.5-r572"
+#define MM_VERSION "2.7-r654"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -43,6 +43,8 @@ static struct option long_options[] = {
 	{ "end-bonus",      required_argument, 0, 0 },
 	{ "no-pairing",     no_argument,       0, 0 },
 	{ "splice-flank",   optional_argument, 0, 0 },
+	{ "idx-no-seq",     no_argument,       0, 0 },
+	{ "end-seed-pen",   required_argument, 0, 0 },   // 21
 	{ "help",           no_argument,       0, 'h' },
 	{ "max-intron-len", required_argument, 0, 'G' },
 	{ "version",        no_argument,       0, 'V' },
@@ -94,7 +96,7 @@ int main(int argc, char *argv[])
 	while ((c = getopt_long(argc, argv, opt_str, long_options, &long_idx)) >= 0) {
 		if (c == 'w') ipt.w = atoi(optarg);
 		else if (c == 'k') ipt.k = atoi(optarg);
-		else if (c == 'H') ipt.is_hpc = 1;
+		else if (c == 'H') ipt.flag |= MM_I_HPC;
 		else if (c == 'd') fnw = optarg; // the above are indexing related options, except -I
 		else if (c == 'r') opt.bw = (int)mm_parse_num(optarg);
 		else if (c == 't') n_threads = atoi(optarg);
@@ -132,12 +134,14 @@ int main(int argc, char *argv[])
 		else if (c == 0 && long_idx == 6) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
 		else if (c == 0 && long_idx == 7) opt.max_chain_skip = atoi(optarg); // --max-chain-skip
 		else if (c == 0 && long_idx == 8) opt.min_ksw_len = atoi(optarg); // --min-dp-len
-		else if (c == 0 && long_idx == 9) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ; // --print-aln-seq
+		else if (c == 0 && long_idx == 9) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
 		else if (c == 0 && long_idx ==10) opt.flag |= MM_F_SPLICE; // --splice
 		else if (c == 0 && long_idx ==12) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
 		else if (c == 0 && long_idx ==13) opt.flag |= MM_F_SR; // --sr
 		else if (c == 0 && long_idx ==17) opt.end_bonus = atoi(optarg); // --end-bonus
 		else if (c == 0 && long_idx ==18) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
+		else if (c == 0 && long_idx ==20) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
+		else if (c == 0 && long_idx ==21) opt.anchor_ext_shift = atoi(optarg); // --end-seed-pen
 		else if (c == 0 && long_idx == 14) { // --frag
 			if (optarg == 0 || strcmp(optarg, "yes") == 0 || strcmp(optarg, "y") == 0)
 				opt.flag |= MM_F_FRAG_MODE;
@@ -196,6 +200,10 @@ int main(int argc, char *argv[])
 		fprintf(stderr, "[ERROR]\033[1;31m --splice and --frag should not be specified at the same time.\033[0m\n");
 		return 1;
 	}
+	if (!fnw && !(opt.flag&MM_F_CIGAR))
+		ipt.flag |= MM_I_NO_SEQ;
+	if (mm_check_opt(&ipt, &opt) < 0)
+		return 1;
 
 	if (argc == optind || fp_help == stdout) {
 		fprintf(fp_help, "Usage: minimap2 [options] <target.fa>|<target.idx> [query.fa] [...]\n");
@@ -259,11 +267,18 @@ int main(int argc, char *argv[])
 	}
 	if (!idx_rdr->is_idx && fnw == 0 && argc - optind < 2) {
 		fprintf(stderr, "[ERROR] missing input: please specify a query file to map or option -d to keep the index\n");
+		mm_idx_reader_close(idx_rdr);
 		return 1;
 	}
 	if (opt.best_n == 0 && (opt.flag&MM_F_CIGAR) && mm_verbose >= 2)
 		fprintf(stderr, "[WARNING]\033[1;31m `-N 0' reduces alignment accuracy. Please use --secondary=no to suppress secondary alignments.\033[0m\n");
 	while ((mi = mm_idx_reader_read(idx_rdr, n_threads)) != 0) {
+		if ((opt.flag & MM_F_CIGAR) && (mi->flag & MM_I_NO_SEQ)) {
+			fprintf(stderr, "[ERROR] the prebuilt index doesn't contain sequences.\n");
+			mm_idx_destroy(mi);
+			mm_idx_reader_close(idx_rdr);
+			return 1;
+		}
 		if ((opt.flag & MM_F_OUT_SAM) && idx_rdr->n_parts == 1) {
 			if (mm_idx_reader_eof(idx_rdr)) {
 				mm_write_sam_hdr(mi, rg, MM_VERSION, argc, argv);

map.c

@@ -36,6 +36,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->end_bonus = -1;
 	opt->min_dp_max = opt->min_chain_score * opt->a;
 	opt->min_ksw_len = 200;
+	opt->anchor_ext_len = 20, opt->anchor_ext_shift = 6;
 	opt->mini_batch_size = 500000000;
 
 	opt->pe_ori = 0; // FF
@@ -64,29 +65,29 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		mm_idxopt_init(io);
 		mm_mapopt_init(mo);
 	} else if (strcmp(preset, "ava-ont") == 0) {
-		io->is_hpc = 0, io->k = 15, io->w = 5;
+		io->flag = 0, io->k = 15, io->w = 5;
 		mo->flag |= MM_F_AVA | MM_F_NO_SELF;
 		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
 	} else if (strcmp(preset, "ava-pb") == 0) {
-		io->is_hpc = 1, io->k = 19, io->w = 5;
+		io->flag |= MM_I_HPC, io->k = 19, io->w = 5;
 		mo->flag |= MM_F_AVA | MM_F_NO_SELF;
 		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
 	} else if (strcmp(preset, "map10k") == 0 || strcmp(preset, "map-pb") == 0) {
-		io->is_hpc = 1, io->k = 19;
+		io->flag |= MM_I_HPC, io->k = 19;
 	} else if (strcmp(preset, "map-ont") == 0) {
-		io->is_hpc = 0, io->k = 15;
+		io->flag = 0, io->k = 15;
 	} else if (strcmp(preset, "asm5") == 0) {
-		io->is_hpc = 0, io->k = 19, io->w = 19;
+		io->flag = 0, io->k = 19, io->w = 19;
 		mo->a = 1, mo->b = 19, mo->q = 39, mo->q2 = 81, mo->e = 3, mo->e2 = 1, mo->zdrop = 200;
 		mo->min_dp_max = 200;
 		mo->best_n = 50;
 	} else if (strcmp(preset, "asm10") == 0) {
-		io->is_hpc = 0, io->k = 19, io->w = 19;
+		io->flag = 0, io->k = 19, io->w = 19;
 		mo->a = 1, mo->b = 9, mo->q = 16, mo->q2 = 41, mo->e = 2, mo->e2 = 1, mo->zdrop = 200;
 		mo->min_dp_max = 200;
 		mo->best_n = 50;
 	} else if (strcmp(preset, "short") == 0 || strcmp(preset, "sr") == 0) {
-		io->is_hpc = 0, io->k = 21, io->w = 11;
+		io->flag = 0, io->k = 21, io->w = 11;
 		mo->flag |= MM_F_SR | MM_F_FRAG_MODE | MM_F_NO_PRINT_2ND | MM_F_2_IO_THREADS;
 		mo->pe_ori = 0<<1|1; // FR
 		mo->a = 2, mo->b = 8, mo->q = 12, mo->e = 2, mo->q2 = 24, mo->e2 = 1;
@@ -104,7 +105,7 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		mo->max_occ = 5000;
 		mo->mini_batch_size = 50000000;
 	} else if (strcmp(preset, "splice") == 0 || strcmp(preset, "cdna") == 0) {
-		io->is_hpc = 0, io->k = 15, io->w = 5;
+		io->flag = 0, io->k = 15, io->w = 5;
 		mo->flag |= MM_F_SPLICE | MM_F_SPLICE_FOR | MM_F_SPLICE_REV | MM_F_SPLICE_FLANK;
 		mo->max_gap = 2000, mo->max_gap_ref = mo->bw = 200000;
 		mo->a = 1, mo->b = 2, mo->q = 2, mo->e = 1, mo->q2 = 32, mo->e2 = 0;
@@ -114,6 +115,21 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 	return 0;
 }
 
+int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
+{
+	if ((mo->q != mo->q2 || mo->e != mo->e2) && !(mo->e > mo->e2 && mo->q + mo->e < mo->q2 + mo->e2)) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m dual gap penalties violating E1>E2 and O1+E1<O2+E2\033[0m\n");
+		return -2;
+	}
+	if ((mo->q + mo->e) + (mo->q2 + mo->e2) > 127) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m scoring system violating ({-O}+{-E})+({-O2}+{-E2}) <= 127\033[0m\n");
+		return -1;
+	}
+	return 0;
+}
+
 typedef struct {
 	uint32_t n;
 	uint32_t qpos;
@@ -122,8 +138,6 @@ typedef struct {
 } mm_match_t;
 
 struct mm_tbuf_s {
-	sdust_buf_t *sdb;
-	mm128_v mini;
 	void *km;
 };
 
@@ -132,24 +146,23 @@ mm_tbuf_t *mm_tbuf_init(void)
 	mm_tbuf_t *b;
 	b = (mm_tbuf_t*)calloc(1, sizeof(mm_tbuf_t));
 	if (!(mm_dbg_flag & 1)) b->km = km_init();
-	b->sdb = sdust_buf_init(b->km);
 	return b;
 }
 
 void mm_tbuf_destroy(mm_tbuf_t *b)
 {
 	if (b == 0) return;
-	kfree(b->km, b->mini.a);
-	sdust_buf_destroy(b->sdb);
 	km_destroy(b->km);
 	free(b);
 }
 
-static int mm_dust_minier(int n, mm128_t *a, int l_seq, const char *seq, int sdust_thres, sdust_buf_t *sdb)
+static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *seq, int sdust_thres)
 {
 	int n_dreg, j, k, u = 0;
 	const uint64_t *dreg;
-	if (sdust_thres <= 0 || sdb == 0) return n;
+	sdust_buf_t *sdb;
+	if (sdust_thres <= 0) return n;
+	sdb = sdust_buf_init(km);
 	dreg = sdust_core((const uint8_t*)seq, l_seq, sdust_thres, 64, &n_dreg, sdb);
 	for (j = k = 0; j < n; ++j) { // squeeze out minimizers that significantly overlap with LCRs
 		int32_t qpos = (uint32_t)a[j].y>>1, span = a[j].x&0xff;
@@ -163,45 +176,49 @@ static int mm_dust_minier(int n, mm128_t *a, int l_seq, const char *seq, int sdu
 				l += ee - ss;
 			}
 			if (l <= span>>1) a[k++] = a[j]; // keep the minimizer if less than half of it falls in masked region
-		}
+		} else a[k++] = a[j];
 	}
+	sdust_buf_destroy(sdb);
 	return k; // the new size
 }
 
-static void collect_minimizers(const mm_mapopt_t *opt, const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, mm_tbuf_t *b)
+static void collect_minimizers(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, mm128_v *mv)
 {
 	int i, j, n, sum = 0;
-	b->mini.n = 0;
+	mv->n = 0;
 	for (i = n = 0; i < n_segs; ++i) {
-		mm_sketch(b->km, seqs[i], qlens[i], mi->w, mi->k, i, mi->is_hpc, &b->mini);
-		for (j = n; j < b->mini.n; ++j)
-			b->mini.a[j].y += sum << 1;
+		mm_sketch(km, seqs[i], qlens[i], mi->w, mi->k, i, mi->flag&MM_I_HPC, mv);
+		for (j = n; j < mv->n; ++j)
+			mv->a[j].y += sum << 1;
 		if (opt->sdust_thres > 0) // mask low-complexity minimizers
-			b->mini.n = n + mm_dust_minier(b->mini.n - n, b->mini.a + n, qlens[i], seqs[i], opt->sdust_thres, b->sdb);
-		sum += qlens[i], n = b->mini.n;
+			mv->n = n + mm_dust_minier(km, mv->n - n, mv->a + n, qlens[i], seqs[i], opt->sdust_thres);
+		sum += qlens[i], n = mv->n;
 	}
 }
 
-static mm128_t *collect_seed_hits(const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, int qlen, int64_t *n_a, int *rep_len, mm_tbuf_t *b)
+static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, const mm128_v *mv, int qlen, int64_t *n_a, int *rep_len,
+								  int *n_mini_pos, uint64_t **mini_pos)
 {
 	int rep_st = 0, rep_en = 0, i;
 	mm_match_t *m;
 	mm128_t *a;
 
-	m = (mm_match_t*)kmalloc(b->km, b->mini.n * sizeof(mm_match_t));
-	for (i = 0; i < b->mini.n; ++i) {
+	*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) {
 		int t;
-		mm128_t *p = &b->mini.a[i];
+		mm128_t *p = &mv->a[i];
 		m[i].qpos = (uint32_t)p->y;
 		m[i].cr = mm_idx_get(mi, p->x>>8, &t);
 		m[i].n = t;
 		m[i].seg_id = p->y >> 32;
 	}
-	for (i = 0, *n_a = 0; i < b->mini.n; ++i) // find the length of a[]
+	for (i = 0, *n_a = 0; i < mv->n; ++i) // find the length of a[]
 		if (m[i].n < max_occ) *n_a += m[i].n;
-	a = (mm128_t*)kmalloc(b->km, *n_a * sizeof(mm128_t));
-	for (i = *rep_len = 0, *n_a = 0; i < b->mini.n; ++i) {
-		mm128_t *p = &b->mini.a[i];
+	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;
@@ -213,8 +230,9 @@ static mm128_t *collect_seed_hits(const mm_mapopt_t *opt, int max_occ, const mm_
 			} else rep_en = en;
 			continue;
 		}
-		if (i > 0 && p->x>>8 == b->mini.a[i - 1].x>>8) is_tandem = 1;
-		if (i < b->mini.n - 1 && p->x>>8 == b->mini.a[i + 1].x>>8) is_tandem = 1;
+		(*mini_pos)[(*n_mini_pos)++] = (uint64_t)q_span<<32 | q->qpos>>1;
+		if (i > 0 && p->x>>8 == mv->a[i - 1].x>>8) is_tandem = 1;
+		if (i < mv->n - 1 && p->x>>8 == mv->a[i + 1].x>>8) is_tandem = 1;
 		for (k = 0; k < q->n; ++k) {
 			int32_t rpos = (uint32_t)r[k] >> 1;
 			mm128_t *p;
@@ -240,7 +258,7 @@ static mm128_t *collect_seed_hits(const mm_mapopt_t *opt, int max_occ, const mm_
 		}
 	}
 	*rep_len += rep_en - rep_st;
-	kfree(b->km, m);
+	kfree(km, m);
 	return a;
 }
 
@@ -258,7 +276,7 @@ static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_i
 static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *km, int qlen, const char *seq, const char *qual, int *n_regs, mm_reg1_t *regs, mm128_t *a)
 {
 	if (!(opt->flag & MM_F_CIGAR)) return regs;
-	regs = mm_align_skeleton(km, opt, mi, qlen, seq, qual, n_regs, regs, a); // this calls mm_filter_regs()
+	regs = mm_align_skeleton(km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs()
 	if (!(opt->flag & MM_F_AVA)) {
 		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b);
 		mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
@@ -269,13 +287,15 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k
 
 void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, const char **quals, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
 {
-	int i, j, rep_len, qlen_sum, n_regs0;
+	int i, j, rep_len, qlen_sum, n_regs0, n_mini_pos;
 	int max_chain_gap_qry, max_chain_gap_ref, is_splice = !!(opt->flag & MM_F_SPLICE), is_sr = !!(opt->flag & MM_F_SR);
 	uint32_t hash;
 	int64_t n_a;
-	uint64_t *u;
+	uint64_t *u, *mini_pos;
 	mm128_t *a;
+	mm128_v mv = {0,0,0};
 	mm_reg1_t *regs0;
+	km_stat_t kmst;
 
 	for (i = 0, qlen_sum = 0; i < n_segs; ++i)
 		qlen_sum += qlens[i], n_regs[i] = 0, regs[i] = 0;
@@ -286,8 +306,8 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	hash ^= __ac_Wang_hash(qlen_sum) + __ac_Wang_hash(opt->seed);
 	hash  = __ac_Wang_hash(hash);
 
-	collect_minimizers(opt, mi, n_segs, qlens, seqs, b);
-	a = collect_seed_hits(opt, opt->mid_occ, mi, qname, qlen_sum, &n_a, &rep_len, b);
+	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 (mm_dbg_flag & MM_DBG_PRINT_SEED) {
@@ -327,7 +347,8 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		if (rechain) { // redo chaining with a higher max_occ threshold
 			kfree(b->km, a);
 			kfree(b->km, u);
-			a = collect_seed_hits(opt, opt->max_occ, mi, qname, qlen_sum, &n_a, &rep_len, b);
+			kfree(b->km, mini_pos);
+			a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 			radix_sort_128x(a, a + n_a);
 			a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 		}
@@ -342,6 +363,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 						i == regs0[j].as? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));
 
 	chain_post(opt, max_chain_gap_ref, mi, b->km, qlen_sum, n_segs, qlens, &n_regs0, regs0, a);
+	if (!is_sr) mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos);
 
 	if (n_segs == 1) { // uni-segment
 		regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], quals? quals[0] : 0, &n_regs0, regs0, a);
@@ -361,8 +383,21 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 			mm_pair(b->km, max_chain_gap_ref, opt->pe_bonus, opt->a * 2 + opt->b, opt->a, qlens, n_regs, regs); // pairing
 	}
 
+	kfree(b->km, mv.a);
 	kfree(b->km, a);
 	kfree(b->km, u);
+	kfree(b->km, mini_pos);
+
+	if (b->km) {
+		km_stat(b->km, &kmst);
+		if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
+			fprintf(stderr, "QM\t%s\t%d\tcap=%ld,nCore=%ld,largest=%ld\n", qname, qlen_sum, kmst.capacity, kmst.n_cores, kmst.largest);
+		assert(kmst.n_blocks == kmst.n_cores); // otherwise, there is a memory leak
+		if (kmst.largest > 1U<<28) {
+			km_destroy(b->km);
+			b->km = km_init();
+		}
+	}
 }
 
 mm_reg1_t *mm_map(const mm_idx_t *mi, int qlen, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
@@ -396,14 +431,13 @@ typedef struct {
 static void worker_for(void *_data, long i, int tid) // kt_for() callback
 {
     step_t *s = (step_t*)_data;
-	int *qlens, j, off = s->seg_off[i], pe_ori = s->p->opt->pe_ori, is_sr = !!(s->p->opt->flag & MM_F_SR);
-	const char **qseqs, **quals = 0;
+	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];
 	mm_tbuf_t *b = s->buf[tid];
+	assert(s->n_seg[i] <= MM_MAX_SEG);
+	memset(quals, 0, sizeof(char*) * MM_MAX_SEG);
 	if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
-		fprintf(stderr, "QR\t%s\t%d\n", s->seq[off].name, tid);
-	qlens = (int*)kmalloc(b->km, s->n_seg[i] * sizeof(int));
-	qseqs = (const char**)kmalloc(b->km, s->n_seg[i] * sizeof(const char**));
-	quals = (const char**)kmalloc(b->km, s->n_seg[i] * sizeof(const char**));
+		fprintf(stderr, "QR\t%s\t%d\t%d\n", s->seq[off].name, tid, s->seq[off].l_seq);
 	for (j = 0; j < s->n_seg[i]; ++j) {
 		if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1))))
 			mm_revcomp_bseq(&s->seq[off + j]);
@@ -429,9 +463,6 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
 				r->rev = !r->rev;
 			}
 		}
-	kfree(b->km, qlens);
-	kfree(b->km, qseqs);
-	kfree(b->km, quals);
 }
 
 static void *worker_pipeline(void *shared, int step, void *in)

minimap.h

@@ -26,6 +26,10 @@
 #define MM_F_SPLICE_FLANK  0x40000
 #define MM_F_SOFTCLIP      0x80000
 
+#define MM_I_HPC          0x1
+#define MM_I_NO_SEQ       0x2
+#define MM_I_NO_NAME      0x4
+
 #define MM_IDX_MAGIC   "MMI\2"
 
 #define MM_MAX_SEG       255
@@ -46,7 +50,7 @@ typedef struct {
 } mm_idx_seq_t;
 
 typedef struct {
-	int32_t b, w, k, is_hpc;
+	int32_t b, w, k, flag;
 	uint32_t n_seq;            // number of reference sequences
 	mm_idx_seq_t *seq;         // sequence name, length and offset
 	uint32_t *S;               // 4-bit packed sequence
@@ -64,23 +68,25 @@ typedef struct {
 } mm_extra_t;
 
 typedef struct {
-	int32_t id;                     // ID for internal uses (see also parent below)
-	uint32_t cnt:28, rev:1, seg_split:1, sam_pri:1, proper_frag:1; // number of minimizers; if on the reverse strand
-	uint32_t rid:31, inv:1;         // reference index; if this is an alignment from inversion rescue
-	int32_t score;                  // DP alignment score
-	int32_t qs, qe, rs, re;         // query start and end; reference start and end
-	int32_t parent, subsc;          // parent==id if primary; best alternate mapping score
-	int32_t as;                     // offset in the a[] array (for internal uses only)
-	int32_t mlen, blen;             // seeded exact match length; seeded alignment block length
-	uint32_t mapq:8, split:2, n_sub:22; // mapQ; split pattern; number of suboptimal mappings
-	uint32_t pe_thru:1, score0:31;
+	int32_t id;             // ID for internal uses (see also parent below)
+	int32_t cnt;            // number of minimizers; if on the reverse strand
+	int32_t rid;            // reference index; if this is an alignment from inversion rescue
+	int32_t score;          // DP alignment score
+	int32_t qs, qe, rs, re; // query start and end; reference start and end
+	int32_t parent, subsc;  // parent==id if primary; best alternate mapping score
+	int32_t as;             // offset in the a[] array (for internal uses only)
+	int32_t mlen, blen;     // seeded exact match length; seeded alignment block length
+	int32_t n_sub;          // number of suboptimal mappings
+	int32_t score0;         // initial chaining score (before chain merging/spliting)
+	uint32_t mapq:8, split:2, rev:1, inv:1, sam_pri:1, proper_frag:1, pe_thru:1, seg_split:1, dummy:16;
 	uint32_t hash;
+	float div;
 	mm_extra_t *p;
 } mm_reg1_t;
 
 // indexing and mapping options
 typedef struct {
-	short k, w, is_hpc, bucket_bits;
+	short k, w, flag, bucket_bits;
 	int mini_batch_size;
 	uint64_t batch_size;
 } mm_idxopt_t;
@@ -110,6 +116,7 @@ typedef struct {
 	int end_bonus;
 	int min_dp_max;  // drop an alignment if the score of the max scoring segment is below this threshold
 	int min_ksw_len;
+	int anchor_ext_len, anchor_ext_shift;
 
 	int pe_ori, pe_bonus;
 
@@ -148,6 +155,7 @@ extern double mm_realtime0; // wall-clock timer
  * @return 0 if success; -1 if _present_ unknown
  */
 int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo);
+int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo);
 
 /**
  * Update mm_mapopt_t::mid_occ via mm_mapopt_t::mid_occ_frac
@@ -200,6 +208,21 @@ void mm_idx_reader_close(mm_idx_reader_t *r);
 
 int mm_idx_reader_eof(const mm_idx_reader_t *r);
 
+/**
+ * Create an index from strings in memory
+ *
+ * @param w            minimizer window size
+ * @param k            minimizer k-mer size
+ * @param is_hpc       use HPC k-mer if true
+ * @param bucket_bits  number of bits for the first level of the hash table
+ * @param n            number of sequences
+ * @param seq          sequences in A/C/G/T
+ * @param name         sequence names; could be NULL
+ *
+ * @return minimap2 index
+ */
+mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const char **seq, const char **name);
+
 /**
  * Print index statistics to stderr
  *
@@ -269,7 +292,7 @@ int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_
 
 // deprecated APIs for backward compatibility
 void mm_mapopt_init(mm_mapopt_t *opt);
-mm_idx_t *mm_idx_build(const char *fn, int w, int k, int is_hpc, int n_threads);
+mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads);
 
 #ifdef __cplusplus
 }

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "11 November 2017" "minimap2-2.5 (r572)" "Bioinformatics tools"
+.TH minimap2 1 "9 January 2018" "minimap2-2.7 (r654)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -99,6 +99,14 @@ multiple times to map it against each batch of target sequences.
 may be ending with k/K/m/M/g/G. NB: mapping quality is incorrect given a
 multi-part index.
 .TP
+.B --idx-no-seq
+Don't store target sequences in the index. It saves disk space and memory but
+the index generated with this option will not work with
+.B -a
+or
+.BR -c .
+When base-level alignment is not requested, this option is automatically applied.
+.TP
 .BI -d \ FILE
 Save the minimizer index of
 .I target.fa
@@ -264,6 +272,18 @@ However, the SIRV control does not honor this trend
 on SIRV data, please add
 .B --splice-flank=no
 to the command line.
+.TP
+.BI --end-seed-pen \ INT
+Drop a terminal anchor if
+.IR s <log( g )+ INT ,
+where
+.I s
+is the local alignment score around the anchor and
+.I g
+the length of the terminal gap in the chain. This option is only effective
+with
+.BR --splice .
+It helps to avoid tiny terminal exons. [6]
 .SS Input/output options
 .TP 10
 .B -a
@@ -461,7 +481,7 @@ cb | cb | cb
 r | c | l .
 Tag	Type	Description
 _
-tp	A	Type of aln: P/primary, S/secondary and I/inversion
+tp	A	Type of aln: P/primary, S/secondary and I,i/inversion
 cm	i	Number of minimizers on the chain
 s1	i	Chaining score
 s2	i	Chaining score of the best secondary chain

misc/cnt-feat.js

@@ -0,0 +1,258 @@
+/*******************************
+ * Command line option parsing *
+ *******************************/
+
+var getopt = function(args, ostr) {
+	var oli; // option letter list index
+	if (typeof(getopt.place) == 'undefined')
+		getopt.ind = 0, getopt.arg = null, getopt.place = -1;
+	if (getopt.place == -1) { // update scanning pointer
+		if (getopt.ind >= args.length || args[getopt.ind].charAt(getopt.place = 0) != '-') {
+			getopt.place = -1;
+			return null;
+		}
+		if (getopt.place + 1 < args[getopt.ind].length && args[getopt.ind].charAt(++getopt.place) == '-') { // found "--"
+			++getopt.ind;
+			getopt.place = -1;
+			return null;
+		}
+	}
+	var optopt = args[getopt.ind].charAt(getopt.place++); // character checked for validity
+	if (optopt == ':' || (oli = ostr.indexOf(optopt)) < 0) {
+		if (optopt == '-') return null; //  if the user didn't specify '-' as an option, assume it means null.
+		if (getopt.place < 0) ++getopt.ind;
+		return '?';
+	}
+	if (oli+1 >= ostr.length || ostr.charAt(++oli) != ':') { // don't need argument
+		getopt.arg = null;
+		if (getopt.place < 0 || getopt.place >= args[getopt.ind].length) ++getopt.ind, getopt.place = -1;
+	} else { // need an argument
+		if (getopt.place >= 0 && getopt.place < args[getopt.ind].length)
+			getopt.arg = args[getopt.ind].substr(getopt.place);
+		else if (args.length <= ++getopt.ind) { // no arg
+			getopt.place = -1;
+			if (ostr.length > 0 && ostr.charAt(0) == ':') return ':';
+			return '?';
+		} else getopt.arg = args[getopt.ind]; // white space
+		getopt.place = -1;
+		++getopt.ind;
+	}
+	return optopt;
+}
+
+/***********************
+ * Interval operations *
+ ***********************/
+
+Interval = {};
+
+Interval.sort = function(a)
+{
+	if (typeof a[0] == 'number')
+		a.sort(function(x, y) { return x - y });
+	else a.sort(function(x, y) { return x[0] != y[0]? x[0] - y[0] : x[1] - y[1] });
+}
+
+Interval.merge = function(a, sorted)
+{
+	if (typeof sorted == 'undefined') sorted = true;
+	if (!sorted) Interval.sort(a);
+	var k = 0;
+	for (var i = 1; i < a.length; ++i) {
+		if (a[k][1] >= a[i][0])
+			a[k][1] = a[k][1] > a[i][1]? a[k][1] : a[i][1];
+		else a[++k] = a[i].slice(0);
+	}
+	a.length = k + 1;
+}
+
+Interval.dedup = function(a, sorted)
+{
+	if (typeof sorted == 'undefined') sorted = true;
+	if (!sorted) Interval.sort(a);
+	var k = 0;
+	for (var i = 1; i < a.length; ++i)
+		if (a[k][0] != a[i][0] || a[k][1] != a[i][1])
+			a[++k] = a[i].slice(0);
+	a.length = k + 1;
+}
+
+Interval.index_end = function(a, sorted)
+{
+	if (a.length == 0) return;
+	if (typeof sorted == 'undefined') sorted = true;
+	if (!sorted) Interval.sort(a);
+	a[0].push(0);
+	var k = 0, k_en = a[0][1];
+	for (var i = 1; i < a.length; ++i) {
+		if (k_en <= a[i][0]) {
+			for (++k; k < i; ++k)
+				if (a[k][1] > a[i][0])
+					break;
+			k_en = a[k][1];
+		}
+		a[i].push(k);
+	}
+}
+
+Interval.find_intv = function(a, x)
+{
+	var left = -1, right = a.length;
+	if (typeof a[0] == 'number') {
+		while (right - left > 1) {
+			var mid = left + ((right - left) >> 1);
+			if (a[mid] > x) right = mid;
+			else if (a[mid] < x) left = mid;
+			else return mid;
+		}
+	} else {
+		while (right - left > 1) {
+			var mid = left + ((right - left) >> 1);
+			if (a[mid][0] > x) right = mid;
+			else if (a[mid][0] < x) left = mid;
+			else return mid;
+		}
+	}
+	return left;
+}
+
+Interval.find_ovlp = function(a, st, en)
+{
+	if (a.length == 0 || st >= en) return [];
+	var l = Interval.find_intv(a, st);
+	var k = l < 0? 0 : a[l][a[l].length - 1];
+	var b = [];
+	for (var i = k; i < a.length; ++i) {
+		if (a[i][0] >= en) break;
+		else if (st < a[i][1])
+			b.push(a[i]);
+	}
+	return b;
+}
+
+/*****************
+ * Main function *
+ *****************/
+
+function read_bed(fn, to_merge, to_dedup)
+{
+	var file = new File(fn);
+	var buf = new Bytes();
+	var h = {};
+	while (file.readline(buf) >= 0) {
+		var t = buf.toString().split("\t");
+		if (h[t[0]] == null)
+			h[t[0]] = [];
+		var bst = parseInt(t[1]);
+		var ben = parseInt(t[2]);
+		if (t.length >= 12 && /^\d+$/.test(t[9])) {
+			t[9] = parseInt(t[9]);
+			var sz = t[10].split(",");
+			var st = t[11].split(",");
+			for (var i = 0; i < t[9]; ++i) {
+				st[i] = parseInt(st[i]);
+				sz[i] = parseInt(sz[i]);
+				h[t[0]].push([bst + st[i], bst + st[i] + sz[i], 0, 0, 0]);
+			}
+		} else {
+			h[t[0]].push([bst, ben, 0, 0, 0]);
+		}
+	}
+	buf.destroy();
+	file.close();
+	for (var chr in h) {
+		if (to_merge) Interval.merge(h[chr], false);
+		else if (to_dedup) Interval.dedup(h[chr], false);
+		else Interval.sort(h[chr]);
+		Interval.index_end(h[chr]);
+	}
+	return h;
+}
+
+function main(args)
+{
+	var c, print_len = false, to_merge = true, to_dedup = false, fn_excl = null;
+	while ((c = getopt(args, "pde:")) != null) {
+		if (c == 'p') print_len = true;
+		else if (c == 'd') to_dedup = true, to_merge = false;
+		else if (c == 'e') fn_excl = getopt.arg;
+	}
+
+	if (args.length - getopt.ind < 2) {
+		print("Usage: k8 cnt-feat.js [options] <target.bed> <feature.bed>");
+		print("Options:");
+		print("  -e FILE    exclude features overlapping regions in BED FILE []");
+		print("  -p         print number of covered bases for each feature");
+		exit(1);
+	}
+
+	var excl = fn_excl != null? read_bed(fn_excl, true, false) : null;
+	var target = read_bed(args[getopt.ind], to_merge, to_dedup);
+
+	var file, buf = new Bytes();
+	var tot_len = 0, hit_len = 0;
+	file = args[getopt.ind+1] != '-'? new File(args[getopt.ind+1]) : new File();
+	while (file.readline(buf) >= 0) {
+		var t = buf.toString().split("\t");
+		var a = [];
+		var bst = parseInt(t[1]);
+		var ben = parseInt(t[2]);
+		if (t.length >= 12 && /^\d+$/.test(t[9])) { // BED12
+			t[9] = parseInt(t[9]);
+			var sz = t[10].split(",");
+			var st = t[11].split(",");
+			for (var i = 0; i < t[9]; ++i) {
+				st[i] = parseInt(st[i]);
+				sz[i] = parseInt(sz[i]);
+				a.push([bst + st[i], bst + st[i] + sz[i], false]);
+			}
+		} else a.push([bst, ben, false]); // 3-column BED
+		var feat_len = 0;
+		for (var i = 0; i < a.length; ++i) {
+			if (excl != null && excl[t[0]] != null) {
+				var oe = Interval.find_ovlp(excl[t[0]], a[i][0], a[i][1]);
+				if (oe.length > 0)
+					continue;
+			}
+			a[i][2] = true;
+			feat_len += a[i][1] - a[i][0];
+		}
+		tot_len += feat_len;
+		if (target[t[0]] == null) continue;
+		var b = [];
+		for (var i = 0; i < a.length; ++i) {
+			if (!a[i][2]) continue;
+			var o = Interval.find_ovlp(target[t[0]], a[i][0], a[i][1]);
+			for (var j = 0; j < o.length; ++j) {
+				var max_st = o[j][0] > a[i][0]? o[j][0] : a[i][0];
+				var min_en = o[j][1] < a[i][1]? o[j][1] : a[i][1];
+				b.push([max_st, min_en]);
+				o[j][2] += min_en - max_st;
+				++o[j][3];
+				if (max_st == o[j][0] && min_en == o[j][1])
+					++o[j][4];
+			}
+		}
+		// find the length covered
+		var feat_hit_len = 0;
+		if (b.length > 0) {
+			b.sort(function(a,b) {return a[0]-b[0]});
+			var st = b[0][0], en = b[0][1];
+			for (var i = 1; i < b.length; ++i) {
+				if (b[i][0] <= en) en = en > b[i][1]? en : b[i][1];
+				else feat_hit_len += en - st, st = b[i][0], en = b[i][1];
+			}
+			feat_hit_len += en - st;
+		}
+		hit_len += feat_hit_len;
+		if (print_len) print('F', t.slice(0, 4).join("\t"), feat_len, feat_hit_len);
+	}
+	file.close();
+
+	buf.destroy();
+
+	warn("# feature bases: " + tot_len);
+	warn("# feature bases overlapping targets: " + hit_len + ' (' + (100.0 * hit_len / tot_len).toFixed(2) + '%)');
+}
+
+main(arguments);

misc/gff2bed.js

@@ -0,0 +1,150 @@
+var getopt = function(args, ostr) {
+	var oli; // option letter list index
+	if (typeof(getopt.place) == 'undefined')
+		getopt.ind = 0, getopt.arg = null, getopt.place = -1;
+	if (getopt.place == -1) { // update scanning pointer
+		if (getopt.ind >= args.length || args[getopt.ind].charAt(getopt.place = 0) != '-') {
+			getopt.place = -1;
+			return null;
+		}
+		if (getopt.place + 1 < args[getopt.ind].length && args[getopt.ind].charAt(++getopt.place) == '-') { // found "--"
+			++getopt.ind;
+			getopt.place = -1;
+			return null;
+		}
+	}
+	var optopt = args[getopt.ind].charAt(getopt.place++); // character checked for validity
+	if (optopt == ':' || (oli = ostr.indexOf(optopt)) < 0) {
+		if (optopt == '-') return null; //  if the user didn't specify '-' as an option, assume it means null.
+		if (getopt.place < 0) ++getopt.ind;
+		return '?';
+	}
+	if (oli+1 >= ostr.length || ostr.charAt(++oli) != ':') { // don't need argument
+		getopt.arg = null;
+		if (getopt.place < 0 || getopt.place >= args[getopt.ind].length) ++getopt.ind, getopt.place = -1;
+	} else { // need an argument
+		if (getopt.place >= 0 && getopt.place < args[getopt.ind].length)
+			getopt.arg = args[getopt.ind].substr(getopt.place);
+		else if (args.length <= ++getopt.ind) { // no arg
+			getopt.place = -1;
+			if (ostr.length > 0 && ostr.charAt(0) == ':') return ':';
+			return '?';
+		} else getopt.arg = args[getopt.ind]; // white space
+		getopt.place = -1;
+		++getopt.ind;
+	}
+	return optopt;
+}
+
+var c, fn_ucsc_fai = null, is_short = false;
+while ((c = getopt(arguments, "u:s")) != null) {
+	if (c == 'u') fn_ucsc_fai = getopt.arg;
+	else if (c == 's') is_short = true;
+}
+
+if (getopt.ind == arguments.length) {
+	print("Usage: k8 gff2bed.js [-u ucsc-genome.fa.fai] <in.gff>");
+	exit(1);
+}
+
+var ens2ucsc = {};
+if (fn_ucsc_fai != null) {
+	var buf = new Bytes();
+	var file = new File(fn_ucsc_fai);
+	while (file.readline(buf) >= 0) {
+		var t = buf.toString().split("\t");
+		var s = t[0];
+		if (/_(random|alt|decoy)$/.test(s)) {
+			s = s.replace(/_(random|alt|decoy)$/, '');
+			s = s.replace(/^chr\S+_/, '');
+		} else {
+			s = s.replace(/^chrUn_/, '');
+		}
+		s = s.replace(/v(\d+)/, ".$1");
+		if (s != t[0]) ens2ucsc[s] = t[0];
+	}
+	file.close();
+	buf.destroy();
+}
+
+var colors = {
+		'protein_coding':'0,128,255',
+		'lincRNA':'0,192,0',
+		'snRNA':'0,192,0',
+		'miRNA':'0,192,0',
+		'misc_RNA':'0,192,0'
+	};
+
+function print_bed12(exons, cds_st, cds_en, is_short)
+{
+	if (exons.length == 0) return;
+	var name = is_short? exons[0][7] + "|" + exons[0][5] : exons[0].slice(4, 7).join("|");
+	var a = exons.sort(function(a,b) {return a[1]-b[1]});
+	var sizes = [], starts = [], st, en;
+	st = a[0][1];
+	en = a[a.length - 1][2];
+	if (cds_st == 1<<30) cds_st = st;
+	if (cds_en == 0) cds_en = en;
+	if (cds_st < st || cds_en > en)
+		throw Error("inconsistent thick start or end for transcript " + a[0][4]);
+	for (var i = 0; i < a.length; ++i) {
+		sizes.push(a[i][2] - a[i][1]);
+		starts.push(a[i][1] - st);
+	}
+	var color = colors[a[0][5]];
+	if (color == null) color = '196,196,196';
+	print(a[0][0], st, en, name, 1000, a[0][3], cds_st, cds_en, color, a.length, sizes.join(",") + ",", starts.join(",") + ",");
+}
+
+var re_gtf = /(transcript_id|transcript_type|transcript_biotype|gene_name|transcript_name) "([^"]+)";/g;
+var re_gff3 = /(transcript_id|transcript_type|transcript_biotype|gene_name|transcript_name)=([^;]+)/g;
+var buf = new Bytes();
+var file = new File(arguments[getopt.ind]);
+
+var exons = [], cds_st = 1<<30, cds_en = 0, last_id = null;
+while (file.readline(buf) >= 0) {
+	var t = buf.toString().split("\t");
+	if (t[0].charAt(0) == '#') continue;
+	if (t[2] != "CDS" && t[2] != "exon") continue;
+	t[3] = parseInt(t[3]) - 1;
+	t[4] = parseInt(t[4]);
+	var id = null, type = "", gname = "N/A", biotype = "", m, tname = "N/A";
+	while ((m = re_gtf.exec(t[8])) != null) {
+		if (m[1] == "transcript_id") id = m[2];
+		else if (m[1] == "transcript_type") type = m[2];
+		else if (m[1] == "transcript_biotype") biotype = m[2];
+		else if (m[1] == "gene_name") name = m[2];
+		else if (m[1] == "transcript_name") tname = m[2];
+	}
+	while ((m = re_gff3.exec(t[8])) != null) {
+		if (m[1] == "transcript_id") id = m[2];
+		else if (m[1] == "transcript_type") type = m[2];
+		else if (m[1] == "transcript_biotype") biotype = m[2];
+		else if (m[1] == "gene_name") name = m[2];
+		else if (m[1] == "transcript_name") tname = m[2];
+	}
+	if (type == "" && biotype != "") type = biotype;
+	if (id == null) throw Error("No transcript_id");
+	if (id != last_id) {
+		print_bed12(exons, cds_st, cds_en, is_short);
+		exons = [], cds_st = 1<<30, cds_en = 0;
+		last_id = id;
+	}
+	if (t[2] == "CDS") {
+		cds_st = cds_st < t[3]? cds_st : t[3];
+		cds_en = cds_en > t[4]? cds_en : t[4];
+	} else if (t[2] == "exon") {
+		if (fn_ucsc_fai != null) {
+			if (ens2ucsc[t[0]] != null)
+				t[0] = ens2ucsc[t[0]];
+			else if (/^[A-Z]+\d+\.\d+$/.test(t[0]))
+				t[0] = t[0].replace(/([A-Z]+\d+)\.(\d+)/, "chrUn_$1v$2");
+		}
+		exons.push([t[0], t[3], t[4], t[6], id, type, name, tname]);
+	}
+}
+if (last_id != null)
+	print_bed12(exons, cds_st, cds_en, is_short);
+
+file.close();
+buf.destroy();

misc/intron-eval.js

@@ -161,7 +161,8 @@ for (var tid in tr) {
 	if (anno[t[0]] == null) anno[t[0]] = [];
 	var s = t[4];
 	for (var i = 0; i < s.length - 1; ++i) {
-		if (s[i][1] >= s[i+1][0]) throw Error("ERROR: wrong annotation!");
+		if (s[i][1] >= s[i+1][0])
+			warn("WARNING: incorrect annotation for transcript "+tid+" ("+s[i][1]+" >= "+s[i+1][0]+")")
 		anno[t[0]].push([s[i][1], s[i+1][0]]);
 	}
 }

misc/ov-eval.js

@@ -0,0 +1,105 @@
+var getopt = function(args, ostr) {
+	var oli; // option letter list index
+	if (typeof(getopt.place) == 'undefined')
+		getopt.ind = 0, getopt.arg = null, getopt.place = -1;
+	if (getopt.place == -1) { // update scanning pointer
+		if (getopt.ind >= args.length || args[getopt.ind].charAt(getopt.place = 0) != '-') {
+			getopt.place = -1;
+			return null;
+		}
+		if (getopt.place + 1 < args[getopt.ind].length && args[getopt.ind].charAt(++getopt.place) == '-') { // found "--"
+			++getopt.ind;
+			getopt.place = -1;
+			return null;
+		}
+	}
+	var optopt = args[getopt.ind].charAt(getopt.place++); // character checked for validity
+	if (optopt == ':' || (oli = ostr.indexOf(optopt)) < 0) {
+		if (optopt == '-') return null; //  if the user didn't specify '-' as an option, assume it means null.
+		if (getopt.place < 0) ++getopt.ind;
+		return '?';
+	}
+	if (oli+1 >= ostr.length || ostr.charAt(++oli) != ':') { // don't need argument
+		getopt.arg = null;
+		if (getopt.place < 0 || getopt.place >= args[getopt.ind].length) ++getopt.ind, getopt.place = -1;
+	} else { // need an argument
+		if (getopt.place >= 0 && getopt.place < args[getopt.ind].length)
+			getopt.arg = args[getopt.ind].substr(getopt.place);
+		else if (args.length <= ++getopt.ind) { // no arg
+			getopt.place = -1;
+			if (ostr.length > 0 && ostr.charAt(0) == ':') return ':';
+			return '?';
+		} else getopt.arg = args[getopt.ind]; // white space
+		getopt.place = -1;
+		++getopt.ind;
+	}
+	return optopt;
+}
+
+var c, min_ovlp = 2000, min_frac = 0.95, min_mapq = 10;
+while ((c = getopt(arguments, "q:l:f:")) != null) {
+	if (c == 'q') min_mapq = parseInt(getopt.arg);
+	else if (c == 'l') min_ovlp = parseInt(getopt.arg);
+	else if (c == 'f') min_frac = parseFloat(getopt.arg);
+}
+if (arguments.length - getopt.ind < 2) {
+	print("Usage: sort -k6,6 -k8,8n to-ref.paf | k8 ov-eval.js [options] - <ovlp.paf>");
+	print("Options:");
+	print("  -l INT     min overlap length [2000]");
+	print("  -q INT     min mapping quality [10]");
+	print("  -f FLOAT   min fraction of mapped length [0.95]");
+	exit(1);
+}
+
+var buf = new Bytes();
+var file = arguments[getopt.ind] == '-'? new File() : new File(arguments[getopt.ind]);
+var a = [], h = {};
+while (file.readline(buf) >= 0) {
+	var t = buf.toString().split("\t");
+	var is_pri = false;
+	if (parseInt(t[11]) < min_mapq) continue;
+	for (var i = 12; i < t.length; ++i)
+		if (t[i] == 'tp:A:P')
+			is_pri = true;
+	if (!is_pri) continue;
+	for (var i = 1; i <= 3; ++i)
+		t[i] = parseInt(t[i]);
+	for (var i = 6; i <= 8; ++i)
+		t[i] = parseInt(t[i]);
+	if (t[3] - t[2] < min_ovlp || t[8] - t[7] < min_ovlp || (t[3] - t[2]) / t[1] < min_frac)
+		continue;
+	var ctg = t[5], st = t[7], en = t[8];
+	while (a.length > 0) {
+		if (a[0][0] == ctg && a[0][2] > st)
+			break;
+		else a.shift();
+	}
+	for (var j = 0; j < a.length; ++j) {
+		if (a[j][3] == t[0]) continue;
+		var len = (en > a[j][2]? a[j][2] : en) - st;
+		if (len >= min_ovlp) {
+			var key = a[j][3] < t[0]? a[j][3] + "\t" + t[0] : t[0] + "\t" + a[j][3];
+			h[key] = len;
+		}
+	}
+	a.push([ctg, st, en, t[0]]);
+}
+file.close();
+
+file = new File(arguments[getopt.ind + 1]);
+while (file.readline(buf) >= 0) {
+	var t = buf.toString().split("\t");
+	var key = t[0] < t[5]? t[0] + "\t" + t[5] : t[5] + "\t" + t[0];
+	if (h[key] > 0) h[key] = -h[key];
+}
+file.close();
+buf.destroy();
+
+var n_ovlp = 0, n_missing = 0;
+for (var key in h) {
+	++n_ovlp;
+	if (h[key] > 0) ++n_missing;
+}
+print(n_ovlp + " overlaps inferred from the reference mapping");
+print(n_missing + " missed by the read overlapper");
+print((100 * (1 - n_missing / n_ovlp)).toFixed(2) + "% sensitivity");

misc/paf2aln.js

@@ -36,20 +36,21 @@ var getopt = function(args, ostr) {
 	return optopt;
 }
 
-var c, maf_out = false, line_len = 80;
-while ((c = getopt(arguments, "ml:")) != null) {
-	if (c == 'm') maf_out = true;
-	else if (c == 'l') line_len = parseInt(getopt.arg); // TODO: not implemented yet
+var c, line_len = 80, fmt = "aln";
+while ((c = getopt(arguments, "f:l:")) != null) {
+	if (c == 'f') {
+		fmt = getopt.arg;
+		if (fmt != "aln" && fmt != "lastz-cigar" && fmt != "maf")
+			throw Error("format must be one of aln, lastz-cigar and maf");
+	} else if (c == 'l') line_len = parseInt(getopt.arg);
 }
 if (line_len == 0) line_len = 0x7fffffff;
 
 if (getopt.ind == arguments.length) {
-	print("Usage: k8 paf2aln.js [options] <with-cs.paf>");
+	print("Usage: k8 paf2aln.js [options] <in.paf>");
 	print("Options:");
-	print("  -m        MAF output (BLAST-like output by default)");
+	print("  -f STR    output format: aln (BLAST-like), maf or lastz-cigar [aln]");
 	print("  -l INT    line length in BLAST-like output [80]");
-	print("");
-	print("Note: this script only works when minimap2 is run with option '-S'");
 	exit(1);
 }
 
@@ -66,7 +67,7 @@ function padding_str(x, len, right)
 function update_aln(s_ref, s_qry, s_mid, type, seq, slen)
 {
 	var l = type == '*'? 1 : seq.length;
-	if (type == '=') {
+	if (type == '=' || type == ':') {
 		s_ref.set(seq);
 		s_qry.set(seq);
 		s_mid.set(Array(l+1).join("|"));
@@ -96,49 +97,73 @@ function print_aln(rs, qs, strand, slen, elen, s_ref, s_qry, s_mid)
 	var st, en;
 	if (strand == '+') st = qs + slen[1] + 1, en = qs + elen[1];
 	else st = qs - slen[1], en = qs - elen[1] + 1;
-	print(["Qry" + strand + ":", padding_str(st,               10, false), s_qry.toString(), padding_str(en          , 10, true)].join(" "));
+	print(["Qry" + strand + ":", padding_str(st, 10, false), s_qry.toString(), padding_str(en, 10, true)].join(" "));
 }
 
-var s_ref = new Bytes(), s_qry = new Bytes(), s_mid = new Bytes();
-var re = /([=\-\+\*])([A-Za-z]+)/g;
+var s_ref = new Bytes(), s_qry = new Bytes(), s_mid = new Bytes(); // these are used to show padded alignment
+var re_cs = /([:=\-\+\*])(\d+|[A-Za-z]+)/g;
+var re_cg = /(\d+)([MIDNSH])/g;
 
 var buf = new Bytes();
-var file = new File(arguments[getopt.ind]);
-if (maf_out) print("##maf version=1\n");
+var file = arguments[getopt.ind] == "-"? new File() : new File(arguments[getopt.ind]);
+var lineno = 0;
+if (fmt == "maf") print("##maf version=1\n");
 while (file.readline(buf) >= 0) {
 	var m, line = buf.toString();
 	var t = line.split("\t", 12);
-	if ((m = /\tcs:Z:(\S+)/.exec(line)) == null) continue;
-	var cs = m[1];
+	++lineno;
 	s_ref.length = s_qry.length = s_mid.length = 0;
 	var slen = [0, 0], elen = [0, 0];
-	if (maf_out) {
-		while ((m = re.exec(cs)) != null)
-			update_aln(s_ref, s_qry, s_mid, m[1], m[2], elen);
-		if (maf_out) {
-			var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : 0;
-			var len = t[0].length > t[5].length? t[0].length : t[5].length;
-			print("a " + score);
-			print(["s", padding_str(t[5], len, true), padding_str(t[7], 10, false), padding_str(parseInt(t[8]) - parseInt(t[7]), 10, false),
-				   "+", padding_str(t[6], 10, false), s_ref.toString()].join(" "));
-			var qs, qe, ql = parseInt(t[1]);
-			if (t[4] == '+') {
-				qs = parseInt(t[2]);
-				qe = parseInt(t[3]);
-			} else {
-				qs = ql - parseInt(t[3]);
-				qe = ql - parseInt(t[2]);
-			}
-			print(["s", padding_str(t[0], len, true), padding_str(qs, 10, false), padding_str(qe - qs, 10, false),
-				   t[4], padding_str(ql, 10, false), s_qry.toString()].join(" "));
-			print("");
+	if (fmt == "lastz-cigar") { // LASTZ-cigar output
+		var cg = (m = /\tcg:Z:(\S+)/.exec(line)) != null? m[1] : null;
+		if (cg == null) {
+			warn("WARNING: converting to LASTZ-cigar format requires the 'cg' tag, which is absent on line " + lineno);
+			continue;
 		}
-	} else {
-		line = line.replace(/\tc[sg]:Z:\S+/g, "");
+		var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? m[1] : 0;
+		var out = ['cigar:', t[0], t[2], t[3], t[4], t[5], t[7], t[8], '+', score];
+		while ((m = re_cg.exec(cg)) != null)
+			out.push(m[2], m[1]);
+		print(out.join(" "));
+	} else if (fmt == "maf") { // MAF output
+		var cs = (m = /\tcs:Z:(\S+)/.exec(line)) != null? m[1] : null;
+		if (cs == null) {
+			warn("WARNING: converting to MAF requires the 'cs' tag, which is absent on line " + lineno);
+			continue;
+		}
+		while ((m = re_cs.exec(cs)) != null) {
+			if (m[1] == ':')
+				throw Error("converting to MAF only works with 'minimap2 --cs=long'");
+			update_aln(s_ref, s_qry, s_mid, m[1], m[2], elen);
+		}
+		var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : 0;
+		var len = t[0].length > t[5].length? t[0].length : t[5].length;
+		print("a " + score);
+		print(["s", padding_str(t[5], len, true), padding_str(t[7], 10, false), padding_str(parseInt(t[8]) - parseInt(t[7]), 10, false),
+				"+", padding_str(t[6], 10, false), s_ref.toString()].join(" "));
+		var qs, qe, ql = parseInt(t[1]);
+		if (t[4] == '+') {
+			qs = parseInt(t[2]);
+			qe = parseInt(t[3]);
+		} else {
+			qs = ql - parseInt(t[3]);
+			qe = ql - parseInt(t[2]);
+		}
+		print(["s", padding_str(t[0], len, true), padding_str(qs, 10, false), padding_str(qe - qs, 10, false),
+				t[4], padding_str(ql, 10, false), s_qry.toString()].join(" "));
+		print("");
+	} else { // BLAST-like output
+		var cs = (m = /\tcs:Z:(\S+)/.exec(line)) != null? m[1] : null;
+		if (cs == null) {
+			warn("WARNING: converting to BLAST-like alignment requires the 'cs' tag, which is absent on line " + lineno);
+			continue;
+		}
+		line = line.replace(/\tc[sg]:Z:\S+/g, ""); // get rid of cs or cg tags
 		print('>' + line);
 		var rs = parseInt(t[7]), qs = t[4] == '+'? parseInt(t[2]) : parseInt(t[3]);
 		var n_blocks = 0;
-		while ((m = re.exec(cs)) != null) {
+		while ((m = re_cs.exec(cs)) != null) {
+			if (m[1] == ':') m[2] = Array(parseInt(m[2]) + 1).join("=");
 			var start = 0, rest = m[1] == '*'? 1 : m[2].length;
 			while (rest > 0) {
 				var l_proc;

misc/paf2diff.js

@@ -0,0 +1,188 @@
+var getopt = function(args, ostr) {
+	var oli; // option letter list index
+	if (typeof(getopt.place) == 'undefined')
+		getopt.ind = 0, getopt.arg = null, getopt.place = -1;
+	if (getopt.place == -1) { // update scanning pointer
+		if (getopt.ind >= args.length || args[getopt.ind].charAt(getopt.place = 0) != '-') {
+			getopt.place = -1;
+			return null;
+		}
+		if (getopt.place + 1 < args[getopt.ind].length && args[getopt.ind].charAt(++getopt.place) == '-') { // found "--"
+			++getopt.ind;
+			getopt.place = -1;
+			return null;
+		}
+	}
+	var optopt = args[getopt.ind].charAt(getopt.place++); // character checked for validity
+	if (optopt == ':' || (oli = ostr.indexOf(optopt)) < 0) {
+		if (optopt == '-') return null; //  if the user didn't specify '-' as an option, assume it means null.
+		if (getopt.place < 0) ++getopt.ind;
+		return '?';
+	}
+	if (oli+1 >= ostr.length || ostr.charAt(++oli) != ':') { // don't need argument
+		getopt.arg = null;
+		if (getopt.place < 0 || getopt.place >= args[getopt.ind].length) ++getopt.ind, getopt.place = -1;
+	} else { // need an argument
+		if (getopt.place >= 0 && getopt.place < args[getopt.ind].length)
+			getopt.arg = args[getopt.ind].substr(getopt.place);
+		else if (args.length <= ++getopt.ind) { // no arg
+			getopt.place = -1;
+			if (ostr.length > 0 && ostr.charAt(0) == ':') return ':';
+			return '?';
+		} else getopt.arg = args[getopt.ind]; // white space
+		getopt.place = -1;
+		++getopt.ind;
+	}
+	return optopt;
+}
+
+var re_cs = /([:=*+-])(\d+|[A-Za-z]+)/g;
+var c, min_cov_len = 10000, min_var_len = 50000, gap_thres = 50, min_mapq = 5;
+while ((c = getopt(arguments, "l:L:g:q:")) != null) {
+	if (c == 'l') min_cov_len = parseInt(getopt.arg);
+	else if (c == 'L') min_var_len = parseInt(optarg.arg);
+	else if (c == 'g') gap_thres = parseInt(optarg.arg);
+	else if (c == 'q') min_mapq = parseInt(optarg.arg);
+}
+
+if (arguments.length == getopt.ind) {
+	print("Usage: k8 paf2diff.js [options] <with-cs.paf>");
+	print("Options:");
+	print("  -l INT    min alignment length to compute coverage ["+min_cov_len+"]");
+	print("  -L INT    min alignment length to call variants ["+min_var_len+"]");
+	print("  -q INT    min mapping quality ["+min_mapq+"]");
+	print("  -g INT    short/long gap threshold (for statistics only) ["+gap_thres+"]");
+	exit(1);
+}
+
+var file = new File(arguments[getopt.ind]);
+var buf = new Bytes();
+var tot_len = 0, n_sub = [0, 0, 0], n_ins = [0, 0, 0, 0], n_del = [0, 0, 0, 0];
+
+function count_var(o)
+{
+	if (o[3] > 1) return;
+	if (o[5] == '-' && o[6] == '-') return;
+	if (o[5] == '-') { // insertion
+		var l = o[6].length;
+		if (l == 1) ++n_ins[0];
+		else if (l == 2) ++n_ins[1];
+		else if (l < gap_thres) ++n_ins[2];
+		else ++n_ins[3];
+	} else if (o[6] == '-') { // deletion
+		var l = o[5].length;
+		if (l == 1) ++n_del[0];
+		else if (l == 2) ++n_del[1];
+		else if (l < gap_thres) ++n_del[2];
+		else ++n_del[3];
+	} else {
+		++n_sub[0];
+		var s = o[5] + o[6];
+		if (s == 'ag' || s == 'ga' || s == 'ct' || s == 'tc')
+			++n_sub[1];
+		else ++n_sub[2];
+	}
+}
+
+var a = [], out = [];
+var c1_ctg = null, c1_start = 0, c1_end = 0, c1_counted = false, c1_len = 0;
+while (file.readline(buf) >= 0) {
+	var line = buf.toString();
+	if (!/\ts2:i:/.test(line)) continue; // skip secondary alignments
+	var m, t = line.split("\t", 12);
+	for (var i = 6; i <= 11; ++i)
+		t[i] = parseInt(t[i]);
+	if (t[10] < min_cov_len || t[11] < min_mapq) continue;
+	var ctg = t[5], x = t[7], end = t[8];
+	// compute regions covered by 1 contig
+	if (ctg != c1_ctg || x >= c1_end) {
+		if (c1_counted && c1_end > c1_start) {
+			c1_len += c1_end - c1_start;
+			print('R', c1_ctg, c1_start, c1_end);
+		}
+		c1_ctg = ctg, c1_start = x, c1_end = end;
+		c1_counted = (t[10] >= min_var_len);
+	} else if (end > c1_end) { // overlap
+		if (c1_counted && x > c1_start) {
+			c1_len += x - c1_start;
+			print('R', c1_ctg, c1_start, x);
+		}
+		c1_start = c1_end, c1_end = end;
+		c1_counted = (t[10] >= min_var_len);
+	} else { // contained
+		if (c1_counted && x > c1_start) {
+			c1_len += x - c1_start;
+			print('R', c1_ctg, c1_start, x);
+		}
+		c1_start = end;
+	}
+	// output variants ahead of this alignment
+	while (out.length) {
+		if (out[0][0] != ctg || out[0][2] <= x) {
+			count_var(out[0]);
+			print('V', out[0].join("\t"));
+			out.shift();
+		} else break;
+	}
+	// update coverage
+	for (var i = 0; i < out.length; ++i)
+		if (out[i][1] >= x && out[i][2] <= end)
+			++out[i][3];
+	// drop alignments that don't overlap with the current one
+	var k = 0;
+	for (var i = 0; i < a.length; ++i)
+		if (a[0][0] == ctg && a[0][2] > x)
+			a[k++] = a[i];
+	a.length = k;
+	// core loop
+	if (t[10] >= min_var_len) {
+		if ((m = /\tcs:Z:(\S+)/.exec(line)) == null) continue; // no cs tag
+		var cs = m[1];
+		var blen = 0, n_diff = 0;
+		tot_len += t[10];
+		while ((m = re_cs.exec(cs)) != null) {
+			var cov = 1;
+			if (m[1] == '*' || m[1] == '+' || m[1] == '-')
+				for (var i = 0; i < a.length; ++i)
+					if (a[0][2] > x) ++cov;
+			if (m[1] == '=' || m[1] == ':') {
+				var l = m[1] == '='? m[2].length : parseInt(m[2]);
+				x += l, blen += l;
+			} else if (m[1] == '*') {
+				out.push([t[5], x, x+1, cov, t[11], m[2].charAt(0), m[2].charAt(1)]);
+				++x, ++blen, ++n_diff;
+			} else if (m[1] == '+') {
+				out.push([t[5], x, x, cov, t[11], '-', m[2]]);
+				++blen, ++n_diff;
+			} else if (m[1] == '-') {
+				out.push([t[5], x, x + m[2].length, cov, t[11], m[2], '-']);
+				x += m[2].length, ++blen, ++n_diff;
+			}
+		}
+	}
+	a.push([t[5], t[7], t[8]]);
+}
+if (c1_counted && c1_end > c1_start) {
+	c1_len += c1_end - c1_start;
+	print('R', c1_ctg, c1_start, c1_end);
+}
+while (out.length) {
+	count_var(out[0]);
+	print('V', out[0].join("\t"));
+	out.shift();
+}
+
+//warn(tot_len + " alignment columns considered in calling");
+warn(c1_len + " reference bases covered by exactly one contig");
+warn(n_sub[0] + " substitutions; ts/tv = " + (n_sub[1]/n_sub[2]).toFixed(3));
+warn(n_del[0] + " 1bp deletions");
+warn(n_ins[0] + " 1bp insertions");
+warn(n_del[1] + " 2bp deletions");
+warn(n_ins[1] + " 2bp insertions");
+warn(n_del[2] + " [3,"+gap_thres+") deletions");
+warn(n_ins[2] + " [3,"+gap_thres+") insertions");
+warn(n_del[3] + " >="+gap_thres+" deletions");
+warn(n_ins[3] + " >="+gap_thres+" insertions");
+
+buf.destroy();
+file.close();

misc/splice2bed.js

@@ -0,0 +1,148 @@
+var getopt = function(args, ostr) {
+	var oli; // option letter list index
+	if (typeof(getopt.place) == 'undefined')
+		getopt.ind = 0, getopt.arg = null, getopt.place = -1;
+	if (getopt.place == -1) { // update scanning pointer
+		if (getopt.ind >= args.length || args[getopt.ind].charAt(getopt.place = 0) != '-') {
+			getopt.place = -1;
+			return null;
+		}
+		if (getopt.place + 1 < args[getopt.ind].length && args[getopt.ind].charAt(++getopt.place) == '-') { // found "--"
+			++getopt.ind;
+			getopt.place = -1;
+			return null;
+		}
+	}
+	var optopt = args[getopt.ind].charAt(getopt.place++); // character checked for validity
+	if (optopt == ':' || (oli = ostr.indexOf(optopt)) < 0) {
+		if (optopt == '-') return null; //  if the user didn't specify '-' as an option, assume it means null.
+		if (getopt.place < 0) ++getopt.ind;
+		return '?';
+	}
+	if (oli+1 >= ostr.length || ostr.charAt(++oli) != ':') { // don't need argument
+		getopt.arg = null;
+		if (getopt.place < 0 || getopt.place >= args[getopt.ind].length) ++getopt.ind, getopt.place = -1;
+	} else { // need an argument
+		if (getopt.place >= 0 && getopt.place < args[getopt.ind].length)
+			getopt.arg = args[getopt.ind].substr(getopt.place);
+		else if (args.length <= ++getopt.ind) { // no arg
+			getopt.place = -1;
+			if (ostr.length > 0 && ostr.charAt(0) == ':') return ':';
+			return '?';
+		} else getopt.arg = args[getopt.ind]; // white space
+		getopt.place = -1;
+		++getopt.ind;
+	}
+	return optopt;
+}
+
+var colors = ["0,128,255", "255,0,0", "0,192,0"];
+
+function print_lines(a, fmt) {
+	if (a.length == 0) return;
+	if (fmt == "bed") {
+		var n_pri = 0;
+		for (var i = 0; i < a.length; ++i)
+			if (a[i][8] == 0) ++n_pri;
+		if (n_pri > 1) {
+			for (var i = 0; i < a.length; ++i)
+				if (a[i][8] == 0) a[i][8] = 1;
+		} else if (n_pri == 0) {
+			warn("Warning: " + a[0][3] + " doesn't have a primary alignment");
+		}
+		for (var i = 0; i < a.length; ++i) {
+			a[i][8] = colors[a[i][8]];
+			print(a[i].join("\t"));
+		}
+	}
+	a.length = 0;
+}
+
+function main(args) {
+	var re = /(\d+)([MIDNSH])/g;
+	var c, fmt = "bed", fn_name_conv = null;
+	while ((c = getopt(args, "f:n:")) != null) {
+		if (c == 'f') fmt = getopt.arg;
+		else if (c == 'n') fn_name_conv = getopt.arg;
+	}
+	if (getopt.ind == args.length) {
+		warn("Usage: k8 splice2bed.js <in.paf>");
+		exit(1);
+	}
+
+	var conv = null;
+	if (fn_name_conv != null) {
+		conv = new Map();
+		var file = new File(fn_name_conv);
+		var buf = new Bytes();
+		while (file.readline(buf) >= 0) {
+			var t = buf.toString().split("\t");
+			conv.put(t[0], t[1]);
+		}
+		buf.destroy();
+		file.close();
+	}
+
+	var file = new File(args[getopt.ind]);
+	var buf = new Bytes();
+	var a = [];
+	while (file.readline(buf) >= 0) {
+		var line = buf.toString();
+		if (line.charAt(0) == '@') continue; // skip SAM header lines
+		var t = line.split("\t");
+		var is_pri = false, cigar = null, a1;
+		var qname = conv != null? conv.get(t[0]) : null;
+		if (qname != null) t[0] = qname;
+		if (t.length >= 10 && t[4] != '+' && t[4] != '-' && /^\d+/.test(t[1])) { // SAM
+			var flag = parseInt(t[1]);
+			if (flag&1) t[0] += '/' + (flag>>6&3);
+		}
+		if (a.length && a[0][3] != t[0]) {
+			print_lines(a, fmt);
+			a = [];
+		}
+		if (t.length >= 12 && (t[4] == '+' || t[4] == '-')) { // PAF
+			for (var i = 12; i < t.length; ++i) {
+				if (t[i].substr(0, 5) == 'cg:Z:') {
+					cigar = t[i].substr(5);
+				} else if (t[i].substr(0, 5) == 's2:i:') {
+					is_pri = true;
+				}
+			}
+			a1 = [t[5], t[7], t[8], t[0], Math.floor(t[9]/t[10]*1000), t[4]];
+		} else if (t.length >= 10) { // SAM
+			var flag = parseInt(t[1]);
+			if ((flag&4) || a[2] == '*') continue;
+			cigar = t[5];
+			is_pri = (flag&0x100)? false : true;
+			a1 = [t[2], parseInt(t[3])-1, null, t[0], 1000, (flag&16)? '-' : '+'];
+		} else {
+			throw Error("unrecognized input format");
+		}
+		if (cigar == null) throw Error("missing CIGAR");
+		var m, x0 = 0, x = 0, bs = [], bl = [];
+		while ((m = re.exec(cigar)) != null) {
+			if (m[2] == 'M' || m[2] == 'D') {
+				x += parseInt(m[1]);
+			} else if (m[2] == 'N') {
+				bs.push(x0);
+				bl.push(x - x0);
+				x += parseInt(m[1]);
+				x0 = x;
+			}
+		}
+		bs.push(x0);
+		bl.push(x - x0);
+		// write the BED12 line
+		if (a1[2] == null) a1[2] = a1[1] + x;
+		a1.push(a1[1], a1[2]); // thick start/end is the same as start/end
+		a1.push(is_pri? 0 : 2, bs.length, bl.join(",")+",", bs.join(",")+",");
+		a.push(a1);
+	}
+	print_lines(a, fmt);
+	buf.destroy();
+	file.close();
+	if (conv != null) conv.destroy();
+}
+
+main(arguments);

mmpriv.h

@@ -64,7 +64,7 @@ const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);
 int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f);
 mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
-mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, const char *qual, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
+mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
 
 mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a);
 void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a);
@@ -79,6 +79,8 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_re
 void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r);
 void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr);
 
+void mm_est_err(const mm_idx_t *mi, int qlen, int n_regs, mm_reg1_t *regs, const mm128_t *a, int32_t n, const uint64_t *mini_pos);
+
 mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int n_regs0, const mm_reg1_t *regs0, int *n_regs, mm_reg1_t **regs, const mm128_t *a);
 void mm_seg_free(void *km, int n_segs, mm_seg_t *segs);
 void mm_pair(void *km, int max_gap_ref, int dp_bonus, int sub_diff, int match_sc, const int *qlens, int *n_regs, mm_reg1_t **regs);

python/cmappy.pxd

@@ -5,7 +5,7 @@ cdef extern from "minimap.h":
 	# Options
 	#
 	ctypedef struct mm_idxopt_t:
-		short k, w, is_hpc, bucket_bits
+		short k, w, flag, bucket_bits
 		int mini_batch_size
 		uint64_t batch_size
 
@@ -30,6 +30,7 @@ cdef extern from "minimap.h":
 		int end_bonus
 		int min_dp_max
 		int min_ksw_len
+		int anchor_ext_len, anchor_ext_shift
 		int pe_ori, pe_bonus
 		float mid_occ_frac
 		int32_t mid_occ
@@ -51,7 +52,7 @@ cdef extern from "minimap.h":
 		pass
 
 	ctypedef struct mm_idx_t:
-		int32_t b, w, k, is_hpc
+		int32_t b, w, k, flag
 		uint32_t n_seq
 		mm_idx_seq_t *seq
 		uint32_t *S

sdust.c

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

setup.py

@@ -23,7 +23,7 @@ def readme():
 
 setup(
 	name = 'mappy',
-	version = '2.5',
+	version = '2.7',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),
@@ -35,7 +35,7 @@ setup(
     ext_modules = [Extension('mappy',
 		sources = [module_src, 'align.c', 'bseq.c', 'chain.c', 'format.c', 'hit.c', 'index.c', 'pe.c',
 				   'ksw2_extd2_sse.c', 'ksw2_exts2_sse.c', 'ksw2_extz2_sse.c', 'ksw2_ll_sse.c',
-				   'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c'],
+				   'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c', 'esterr.c'],
 		depends = ['minimap.h', 'bseq.h', 'kalloc.h', 'kdq.h', 'khash.h', 'kseq.h', 'ksort.h',
 				   'ksw2.h', 'kthread.h', 'kvec.h', 'mmpriv.h', 'sdust.h',
 				   'python/cmappy.h', 'python/cmappy.pxd'],

sketch.c

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

tex/minimap2.bib

@@ -268,18 +268,10 @@
 	Title = {Mason -- a read simulator for second generation sequencing data},
 	Year = {2010}}
 
-@article{Langmead:2012fk,
-	Author = {Langmead, Ben and Salzberg, Steven L},
-	Journal = {Nat Methods},
-	Pages = {357-9},
-	Title = {Fast gapped-read alignment with Bowtie 2},
-	Volume = {9},
-	Year = {2012}}
-
 @article{Zaharia:2011aa,
 	Author = {Zaharia, Matei and others},
 	Journal = {arXiv:1111:5572},
-	Title = {Faster and More Accurate Sequence Alignment with SNAP},
+	Title = {Faster and More Accurate Sequence Alignment with {SNAP}},
 	Year = {2011}}
 
 @article{Irimia:2008aa,
@@ -305,3 +297,19 @@
 	Title = {Versatile and open software for comparing large genomes},
 	Volume = {5},
 	Year = {2004}}
+
+@article {Li223297,
+	author = {Li, Heng and others},
+	title = {New synthetic-diploid benchmark for accurate variant calling evaluation},
+	year = {2017},
+	note = {doi:10.1101/223297},
+	journal = {bioRxiv}
+}
+
+@article{Berlin:2015xy,
+	Author = {Berlin, Konstantin and others},
+	Journal = {Nat Biotechnol},
+	Pages = {623-30},
+	Title = {Assembling large genomes with single-molecule sequencing and locality-sensitive hashing},
+	Volume = {33},
+	Year = {2015}}

tex/minimap2.tex

@@ -13,7 +13,6 @@
 
 \usepackage{natbib}
 \bibliographystyle{apalike}
-\usepackage{hyperref}
 
 \DeclareMathOperator*{\argmax}{argmax}
 
@@ -138,7 +137,7 @@ $h=50$; even if the heuristic fails, the optimal chain is often close.
 
 \subsubsection{Backtracking}
 Let $P(i)$ be the index of the best predecessor of anchor $i$. It equals 0 if
-$f(i)=w_i$ or $\argmax_j\{f(j)+\eta(j,i)-\gamma(j,i)\}$ otherwise. For each
+$f(i)=w_i$ or $\argmax_j\{f(j)+\alpha(j,i)-\beta(j,i)\}$ otherwise. For each
 anchor $i$ in the descending order of $f(i)$, we apply $P(\cdot)$ repeatedly to
 find its predecessor and mark each visited $i$ as `used', until $P(i)=0$ or we
 reach an already `used' $i$. This way we find all chains with no anchors used
@@ -157,6 +156,53 @@ add the chain to $Q$. In the end, $Q$ contains all the primary chains. We did
 not choose a more sophisticated data structure (e.g. range tree or k-d tree)
 because this step is not the performance bottleneck.
 
+\subsubsection{Estimating per-base sequence divergence}
+Suppose a query sequence harbors $n$ seeds of length $k$, $m$ of which are
+present in a chain. We want to estimate the sequence divergence $\epsilon$
+between the query and the reference sequences in the chain. This is useful
+when base-level alignment is too expensive to perform.
+
+If we model substitutions with a homogeneous Poisson process along the query
+sequence, the probablity of seeing $k$ consecutive bases without substitutions
+is $e^{-k\epsilon}$. On the assumption that all $k$-mers are independent of
+each other, the likelihood function of $\epsilon$ is
+\[
+\mathcal{L}(\epsilon|n,m,k)=e^{-m\cdot k\epsilon}(1-e^{-k\epsilon})^{n-m}
+\]
+The maximum likelihood estimate of $\epsilon$ is
+\[
+\hat{\epsilon}=\frac{1}{k}\log\frac{n}{m}
+\]
+In reality, sequencing errors are sometimes clustered and $k$-mers are not
+independent of each other, especially when we take minimizers as seeds. These
+violate the assumptions in the derivation above. As a result, $\hat{\epsilon}$
+is only approximate and can be biased. It also ignores long deletions from the
+reference sequence. In practice, fortunately, $\hat{\epsilon}$ is often close
+to and strongly correlated with the sequence divergence estimated from
+base-level alignments. On the several datasets used in
+Section~\ref{sec:long-genomic}, the Spearman correlation coefficient is around
+$0.9$.
+
+\subsubsection{Indexing with homopolymer compressed $k$-mers}
+SmartDenovo
+(\href{https://github.com/ruanjue/smartdenovo}{https://github.com/ruanjue/smartdenovo};
+J Ruan, personal communication) indexes reads with homopolymer-compressed (HPC)
+$k$-mers and finds the strategy improves overlap sensitivity for SMRT reads.
+Minimap2 adopts the same heuristic.
+
+The HPC string of a string $s$, denoted by ${\rm HPC}(s)$, is constructed by
+contracting homopolymers in $s$ to a single base.  An HPC $k$-mer of $s$ is a
+$k$-long substring of ${\rm HPC}(s)$. For example, suppose $s={\tt GGATTTTCCA}$,
+${\rm HPC}(s)={\tt GATCA}$ and the first HPC 4-mer is ${\tt GATC}$.
+
+To demonstrate the effectiveness of HPC $k$-mers, we performed read overlapping
+for the example {\it E. coli} SMRT reads from PBcR~\citep{Berlin:2015xy}, using
+different types of $k$-mers. With normal 15bp minimizers per 5bp window,
+minimap2 finds 90.9\% of $\ge$2kb overlaps inferred from the read-to-reference
+alignment. With HPC 19-mers, minimap2 finds 97.4\% of overlaps. It achieves this
+higher sensitivity by indexing 1/3 fewer minimizers, which further helps
+performance. HPC-based indexing reduces the sensitivity for ONT reads, though.
+
 \subsection{Aligning genomic DNA}\label{sec:genomic}
 
 \subsubsection{Alignment with 2-piece affine gap cost}
@@ -397,7 +443,7 @@ consistent paired-end alignments.
 
 \section{Results}
 
-\subsection{Aligning long genomic reads}
+\subsection{Aligning long genomic reads}\label{sec:long-genomic}
 
 \begin{figure}[!tb]
 \centering
@@ -500,14 +546,18 @@ more junctions with a higher percentage being exactly or approximately correct.
 Minimap2 is over 40 times faster than GMAP and SpAln. While STAR is close to
 minimap2 in speed, it does not work well with noisy reads.
 
-We have also evaluated spliced aligners on public Iso-Seq data (human Alzheimer
-brain from \href{http://bit.ly/isoseqpub}{http://bit.ly/isoseqpub}). The
-observation is similar: minimap2 is faster at higher junction accuracy.
-On a private Nanopore Direct RNA data set with $\sim$17\% sequencing error rate
-(N. Loman, personal communication), minimap2 aligned 96\,467 introns
-from 37\,068 mapped reads with 95.4\% of them consistent with human gene
-annotations. In comparison, only 74.8\% of GMAP introns found in known gene
-annotations.
+We have also evaluated spliced aligners on a human Nanopore Direct RNA-seq
+dataset (\href{http://bit.ly/na12878ont}{http://bit.ly/na12878ont}). Minimap2
+aligned 10 million reads in $<$1 wall-clock hour using 16 CPU cores. 94.2\% of
+aligned splice junctions consistent with gene annotations. In comparison,
+GMAP under option `-k 14 -n 0 --min-intronlength 30 --cross-species' is 160
+times slower; 68.7\% of GMAP junctions are found in known gene annotations. The
+percentage increases to 84.1\% if an aligned junction within 10bp from an
+annotated junction is considered to be correct. On a public Iso-Seq dataset
+(human Alzheimer brain from
+\href{http://bit.ly/isoseqpub}{http://bit.ly/isoseqpub}), minimap2 is also
+faster at higher junction accuracy in comparison to other aligners in
+Table~\ref{tab:intron}.
 
 We noted that GMAP and SpAln have not been optimized for noisy reads. We are
 showing the best setting we have experimented, but their developers should be
@@ -551,8 +601,7 @@ with GATK HaplotypeCaller v3.5~\citep{Depristo:2011vn}. This run was sequenced
 from experimentally mixed CHM1 and CHM13 cell lines. Both of them are homozygous
 across the whole genome and have been \emph{de novo} assembled with SMRT reads
 to high quality. This allowed us to construct an independent truth variant
-data set
-(\href{https://github.com/lh3/CHM-eval}{https://github.com/lh3/CHM-eval}) for
+dataset~\citep{Li223297} for
 ERR1341796. In this evaluation, minimap2 has higher SNP false negative rate
 (FNR; 2.5\% of minimap2 vs 2.2\% of BWA-MEM), but fewer false positive SNPs per
 million bases (FPPM; 3.0 vs 3.9), lower 2--50bp INDEL FNR (7.3\% vs 7.5\%) and
@@ -597,7 +646,7 @@ uniqueness and reduce unsuccessful extensions. Minimap2 indexes reference
 k-mers with a hash table instead. Such fixed-length seeds are inferior to
 variable-length seeds in theory, but can be computed much more efficiently in
 practice. When a query sequence has multiple seed hits, we can afford to skip
-some highly repetitive seeds without affecting the final accuracy. This further
+highly repetitive seeds without affecting the final accuracy. This further
 alleviates the concern with the uniqueness of seeds. Hash table is the ideal
 data structure for mapping long query sequences.
 
@@ -605,8 +654,8 @@ data structure for mapping long query sequences.
 We owe a debt of gratitude to H. Suzuki and M. Kasahara for releasing their
 masterpiece and insightful notes before formal publication. We thank M.
 Schatz, P. Rescheneder and F. Sedlazeck for pointing out the limitation of
-BWA-MEM. We are also grateful to early minimap2 testers who have greatly helped
-to suggest features and to fix various issues.
+BWA-MEM. We are also grateful to minimap2 users who have greatly helped to
+suggest features and to fix various issues.
 
 \bibliography{minimap2}