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

@@ -6,12 +6,19 @@ 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,28 @@
+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)
 -----------------------------------
 

README.md

@@ -65,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.6/minimap2-2.6_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.6_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
@@ -352,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

align.c

@@ -210,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)

index.c

@@ -329,7 +329,7 @@ mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int flag, int mini
 	return pl.mi;
 }
 
-mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, 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;
@@ -340,6 +340,49 @@ mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads) //
 	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 *
  *************/

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.6-r623"
+#define MM_VERSION "2.7-r654"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -134,7 +134,7 @@ int main(int argc, char *argv[])
 		else if (c == 0 && long_idx == 6) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, 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
@@ -202,6 +202,8 @@ int main(int argc, char *argv[])
 	}
 	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");

map.c

@@ -115,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;

minimap.h

@@ -155,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
@@ -207,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
  *

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "12 December 2017" "minimap2-2.6 (r623)" "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

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/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/splice2bed.js

@@ -60,15 +60,29 @@ function print_lines(a, fmt) {
 
 function main(args) {
 	var re = /(\d+)([MIDNSH])/g;
-	var c, fmt = "bed";
-	while ((c = getopt(args, "f:")) != null) {
+	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 = [];
@@ -77,11 +91,17 @@ function main(args) {
 		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] == '-')) {
+		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);
@@ -90,7 +110,7 @@ function main(args) {
 				}
 			}
 			a1 = [t[5], t[7], t[8], t[0], Math.floor(t[9]/t[10]*1000), t[4]];
-		} else if (t.length >= 10) {
+		} else if (t.length >= 10) { // SAM
 			var flag = parseInt(t[1]);
 			if ((flag&4) || a[2] == '*') continue;
 			cigar = t[5];
@@ -122,6 +142,7 @@ function main(args) {
 	print_lines(a, fmt);
 	buf.destroy();
 	file.close();
+	if (conv != null) conv.destroy();
 }
 
 main(arguments);

setup.py

@@ -23,7 +23,7 @@ def readme():
 
 setup(
 	name = 'mappy',
-	version = '2.6',
+	version = '2.7',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),

sketch.c

@@ -113,26 +113,26 @@ 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 && buf[j].y != UINT64_MAX) kv_push(mm128_t, km, *p, buf[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 && buf[j].y != UINT64_MAX) kv_push(mm128_t, km, *p, buf[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 && min.y != UINT64_MAX) 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 && buf[j].y != UINT64_MAX) kv_push(mm128_t, km, *p, buf[j]);
+					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)
-					if (min.x == buf[j].x && min.y != buf[j].y && buf[j].y != UINT64_MAX) kv_push(mm128_t, km, *p, buf[j]);
+					if (min.x == buf[j].x && min.y != buf[j].y) kv_push(mm128_t, km, *p, buf[j]);
 			}
 		}
 		if (++buf_pos == w) buf_pos = 0;

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}