Release minimap2-2.1.1 (r341)

old code does not check AVX2

.gitignore

@@ -1,3 +1,5 @@
+.cproject
+.project
 .*.swp
 *.a
 *.o

Makefile

@@ -1,15 +1,16 @@
 CC=			gcc
 CFLAGS=		-g -Wall -O2 -Wc++-compat
 CPPFLAGS=	-DHAVE_KALLOC
-INCLUDES=	-I.
-OBJS=		kthread.o kalloc.o ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_exts2_sse.o ksw2_ll_sse.o \
-			misc.o bseq.o sketch.o sdust.o index.o chain.o align.o hit.o map.o format.o
+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 ksw2_ll_sse.o
 PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite
 LIBS=		-lm -lz -lpthread
 
 ifeq ($(sse2only),)
-	CFLAGS+=-msse4
+	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
 
 .SUFFIXES:.c .o
@@ -21,8 +22,8 @@ all:$(PROG)
 
 extra:all $(PROG_EXTRA)
 
-minimap2:main.o libminimap2.a
-		$(CC) $(CFLAGS) $< -o $@ -L. -lminimap2 $(LIBS)
+minimap2:main.o getopt.o libminimap2.a
+		$(CC) $(CFLAGS) main.o getopt.o -o $@ -L. -lminimap2 $(LIBS)
 
 minimap2-lite:example.o libminimap2.a
 		$(CC) $(CFLAGS) $< -o $@ -L. -lminimap2 $(LIBS)
@@ -30,8 +31,29 @@ minimap2-lite:example.o libminimap2.a
 libminimap2.a:$(OBJS)
 		$(AR) -csru $@ $(OBJS)
 
-sdust:sdust.c kalloc.o kalloc.h kdq.h kvec.h kseq.h sdust.h
-		$(CC) -D_SDUST_MAIN $(CFLAGS) $< kalloc.o -o $@ -lz
+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
+
+ksw2_extz2_sse41.o:ksw2_extz2_sse.c ksw2.h kalloc.h
+		$(CC) -c -msse4 $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
+
+ksw2_extz2_sse2.o:ksw2_extz2_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@
+
+ksw2_extd2_sse41.o:ksw2_extd2_sse.c ksw2.h kalloc.h
+		$(CC) -c -msse4 $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
+
+ksw2_extd2_sse2.o:ksw2_extd2_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@
+
+ksw2_exts2_sse41.o:ksw2_exts2_sse.c ksw2.h kalloc.h
+		$(CC) -c -msse4 $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
+
+ksw2_exts2_sse2.o:ksw2_exts2_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@
+
+ksw2_dispatch.o:ksw2_dispatch.c ksw2.h
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
 
 clean:
 		rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM session*
@@ -46,6 +68,7 @@ bseq.o: bseq.h kseq.h
 chain.o: minimap.h mmpriv.h bseq.h kalloc.h
 example.o: minimap.h kseq.h
 format.o: kalloc.h mmpriv.h minimap.h bseq.h
+getopt.o: getopt.h
 hit.o: mmpriv.h minimap.h bseq.h kalloc.h
 index.o: kthread.h bseq.h minimap.h mmpriv.h kvec.h kalloc.h khash.h
 kalloc.o: kalloc.h
@@ -53,7 +76,7 @@ ksw2_extd2_sse.o: ksw2.h kalloc.h
 ksw2_exts2_sse.o: ksw2.h kalloc.h
 ksw2_extz2_sse.o: ksw2.h kalloc.h
 ksw2_ll_sse.o: ksw2.h kalloc.h
-main.o: bseq.h minimap.h mmpriv.h
+main.o: bseq.h minimap.h mmpriv.h getopt.h
 map.o: kthread.h kvec.h kalloc.h sdust.h mmpriv.h minimap.h bseq.h
 misc.o: minimap.h ksort.h
 sdust.o: kalloc.h kdq.h kvec.h sdust.h

NEWS.md

@@ -1,3 +1,26 @@
+Release 2.1.1-r341 (6 September 2017)
+-------------------------------------
+
+This is a maintenance release that is expected to output identical alignment to
+v2.1. Detailed changes include:
+
+ * Support CPU dispatch. By default, minimap2 is compiled with both SSE2 and
+   SSE4 based implementation of alignment and automatically chooses the right
+   one at runtime. This avoids unexpected errors on older CPUs (#21).
+
+ * Improved Windows support as is requested by Oxford Nanopore (#19). Minimap2
+   now avoids variable-length stacked arrays, eliminates alloca(), ships with
+   getopt_long() and provides timing functions implemented with Windows APIs.
+
+ * Fixed a potential segmentation fault when specifying -k/-w/-H with
+   multi-part index (#23).
+
+ * Fixed two memory leaks in example.c
+
+(2.1.1: 6 September 2017, r341)
+
+
+
 Release 2.1-r311 (25 August 2017)
 ---------------------------------
 

README.md

@@ -10,6 +10,8 @@ cd minimap2 && make
 ./minimap2 -ax map10k MT-human.mmi test/MT-orang.fa > test.sam
 # long-read overlap (no test data)
 ./minimap2 -x ava-pb your-reads.fa your-reads.fa > overlaps.paf
+# spliced alignment (no test data)
+./minimap2 -ax splice ref.fa rna-seq-reads.fa > spliced.sam
 # man page
 man ./minimap2.1
 ```

align.c

@@ -138,8 +138,10 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
 		int i;
 		fprintf(stderr, "===> q=(%d,%d), e=(%d,%d), bw=%d, flag=%d, zdrop=%d <===\n", opt->q, opt->q2, opt->e, opt->e2, w, flag, opt->zdrop);
-		for (i = 0; i < tlen; ++i) fputc("ACGTN"[tseq[i]], stderr); fputc('\n', stderr);
-		for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr); fputc('\n', stderr);
+		for (i = 0; i < tlen; ++i) fputc("ACGTN"[tseq[i]], stderr);
+		fputc('\n', stderr);
+		for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr);
+		fputc('\n', stderr);
 	}
 	if (opt->flag & MM_F_SPLICE)
 		ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, opt->zdrop, flag, ez);

example.c

@@ -35,13 +35,13 @@ int main(int argc, char *argv[])
 	opt.flag |= MM_F_CIGAR; // perform alignment
 	mm_tbuf_t *tbuf = mm_tbuf_init(); // thread buffer; for multi-threading, allocate one tbuf for each thread
 	while (kseq_read(ks) >= 0) { // each kseq_read() call reads one query sequence
-		const mm_reg1_t *reg;
+		mm_reg1_t *reg;
 		int j, i, n_reg;
 		// get all hits for the query
 		reg = mm_map(mi, ks->seq.l, ks->seq.s, &n_reg, tbuf, &opt, 0);
 		// traverse hits and print them out
 		for (j = 0; j < n_reg; ++j) {
-			const mm_reg1_t *r = &reg[j];
+			mm_reg1_t *r = &reg[j];
 			assert(r->p); // with MM_F_CIGAR, this should not be NULL
 			printf("%s\t%d\t%d\t%d\t%c\t", ks->name.s, ks->seq.l, r->qs, r->qe, "+-"[r->rev]);
 			printf("%s\t%d\t%d\t%d\t%d\t%d\t%d\tcg:Z:", mi->seq[r->rid].name, mi->seq[r->rid].len, r->rs, r->re,
@@ -49,7 +49,9 @@ int main(int argc, char *argv[])
 			for (i = 0; i < r->p->n_cigar; ++i) // IMPORTANT: this gives the CIGAR in the aligned regions. NO soft/hard clippings!
 				printf("%d%c", r->p->cigar[i]>>4, "MIDSHN"[r->p->cigar[i]&0xf]);
 			putchar('\n');
+			free(r->p);
 		}
+		free(reg);
 	}
 	mm_tbuf_destroy(tbuf);
 

getopt.c

@@ -0,0 +1,216 @@
+#include <stddef.h>
+#include <stdio.h>
+#include <string.h>
+#include "getopt.h"
+
+char *optarg;
+int optind=1, opterr=1, optopt, __optpos, optreset=0;
+
+#define optpos __optpos
+
+static void __getopt_msg(const char *a, const char *b, const char *c, size_t l)
+{
+	FILE *f = stderr;
+#if !defined(WIN32) && !defined(_WIN32)
+	flockfile(f);
+#endif
+	fputs(a, f);
+	fwrite(b, strlen(b), 1, f);
+	fwrite(c, 1, l, f);
+	fputc('\n', f);
+#if !defined(WIN32) && !defined(_WIN32)
+	funlockfile(f);
+#endif
+}
+
+int getopt(int argc, char * const argv[], const char *optstring)
+{
+	int i, c, d;
+	int k, l;
+	char *optchar;
+
+	if (!optind || optreset) {
+		optreset = 0;
+		__optpos = 0;
+		optind = 1;
+	}
+
+	if (optind >= argc || !argv[optind])
+		return -1;
+
+	if (argv[optind][0] != '-') {
+		if (optstring[0] == '-') {
+			optarg = argv[optind++];
+			return 1;
+		}
+		return -1;
+	}
+
+	if (!argv[optind][1])
+		return -1;
+
+	if (argv[optind][1] == '-' && !argv[optind][2])
+		return optind++, -1;
+
+	if (!optpos) optpos++;
+	c = argv[optind][optpos], k = 1;
+	optchar = argv[optind]+optpos;
+	optopt = c;
+	optpos += k;
+
+	if (!argv[optind][optpos]) {
+		optind++;
+		optpos = 0;
+	}
+
+	if (optstring[0] == '-' || optstring[0] == '+')
+		optstring++;
+
+	i = 0;
+	d = 0;
+	do {
+		d = optstring[i], l = 1;
+		if (l>0) i+=l; else i++;
+	} while (l && d != c);
+
+	if (d != c) {
+		if (optstring[0] != ':' && opterr)
+			__getopt_msg(argv[0], ": unrecognized option: ", optchar, k);
+		return '?';
+	}
+	if (optstring[i] == ':') {
+		if (optstring[i+1] == ':') optarg = 0;
+		else if (optind >= argc) {
+			if (optstring[0] == ':') return ':';
+			if (opterr) __getopt_msg(argv[0],
+				": option requires an argument: ",
+				optchar, k);
+			return '?';
+		}
+		if (optstring[i+1] != ':' || optpos) {
+			optarg = argv[optind++] + optpos;
+			optpos = 0;
+		}
+	}
+	return c;
+}
+
+static void permute(char *const *argv, int dest, int src)
+{
+	char **av = (char **)argv;
+	char *tmp = av[src];
+	int i;
+	for (i=src; i>dest; i--)
+		av[i] = av[i-1];
+	av[dest] = tmp;
+}
+
+static int __getopt_long_core(int argc, char *const *argv, const char *optstring, const struct option *longopts, int *idx, int longonly)
+{
+	optarg = 0;
+	if (longopts && argv[optind][0] == '-' &&
+		((longonly && argv[optind][1] && argv[optind][1] != '-') ||
+		 (argv[optind][1] == '-' && argv[optind][2])))
+	{
+		int colon = optstring[optstring[0]=='+'||optstring[0]=='-']==':';
+		int i, cnt, match;
+		char *opt;
+		for (cnt=i=0; longopts[i].name; i++) {
+			const char *name = longopts[i].name;
+			opt = argv[optind]+1;
+			if (*opt == '-') opt++;
+			for (; *name && *name == *opt; name++, opt++);
+			if (*opt && *opt != '=') continue;
+			match = i;
+			if (!*name) {
+				cnt = 1;
+				break;
+			}
+			cnt++;
+		}
+		if (cnt==1) {
+			i = match;
+			optind++;
+			optopt = longopts[i].val;
+			if (*opt == '=') {
+				if (!longopts[i].has_arg) {
+					if (colon || !opterr)
+						return '?';
+					__getopt_msg(argv[0],
+						": option does not take an argument: ",
+						longopts[i].name,
+						strlen(longopts[i].name));
+					return '?';
+				}
+				optarg = opt+1;
+			} else if (longopts[i].has_arg == required_argument) {
+				if (!(optarg = argv[optind])) {
+					if (colon) return ':';
+					if (!opterr) return '?';
+					__getopt_msg(argv[0],
+						": option requires an argument: ",
+						longopts[i].name,
+						strlen(longopts[i].name));
+					return '?';
+				}
+				optind++;
+			}
+			if (idx) *idx = i;
+			if (longopts[i].flag) {
+				*longopts[i].flag = longopts[i].val;
+				return 0;
+			}
+			return longopts[i].val;
+		}
+		if (argv[optind][1] == '-') {
+			if (!colon && opterr)
+				__getopt_msg(argv[0], cnt ?
+					": option is ambiguous: " :
+					": unrecognized option: ",
+					argv[optind]+2,
+					strlen(argv[optind]+2));
+			optind++;
+			return '?';
+		}
+	}
+	return getopt(argc, argv, optstring);
+}
+
+static int __getopt_long(int argc, char *const *argv, const char *optstring, const struct option *longopts, int *idx, int longonly)
+{
+	int ret, skipped, resumed;
+	if (!optind || optreset) {
+		optreset = 0;
+		__optpos = 0;
+		optind = 1;
+	}
+	if (optind >= argc || !argv[optind]) return -1;
+	skipped = optind;
+	if (optstring[0] != '+' && optstring[0] != '-') {
+		int i;
+		for (i=optind; ; i++) {
+			if (i >= argc || !argv[i]) return -1;
+			if (argv[i][0] == '-' && argv[i][1]) break;
+		}
+		optind = i;
+	}
+	resumed = optind;
+	ret = __getopt_long_core(argc, argv, optstring, longopts, idx, longonly);
+	if (resumed > skipped) {
+		int i, cnt = optind-resumed;
+		for (i=0; i<cnt; i++)
+			permute(argv, skipped, optind-1);
+		optind = skipped + cnt;
+	}
+	return ret;
+}
+
+int getopt_long(int argc, char *const *argv, const char *optstring, const struct option *longopts, int *idx)
+{
+	return __getopt_long(argc, argv, optstring, longopts, idx, 0);
+}
+
+int getopt_long_only(int argc, char *const *argv, const char *optstring, const struct option *longopts, int *idx)
+{
+	return __getopt_long(argc, argv, optstring, longopts, idx, 1);
+}

getopt.h

@@ -0,0 +1,53 @@
+/*
+  Copyright 2005-2014 Rich Felker, et al.
+  
+  Permission is hereby granted, free of charge, to any person obtaining
+  a copy of this software and associated documentation files (the
+  "Software"), to deal in the Software without restriction, including
+  without limitation the rights to use, copy, modify, merge, publish,
+  distribute, sublicense, and/or sell copies of the Software, and to
+  permit persons to whom the Software is furnished to do so, subject to
+  the following conditions:
+  
+  The above copyright notice and this permission notice shall be
+  included in all copies or substantial portions of the Software.
+  
+  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+*/
+
+#ifndef _GETOPT_H
+#define _GETOPT_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int getopt(int, char * const [], const char *);
+extern char *optarg;
+extern int optind, opterr, optopt, optreset;
+
+struct option {
+	const char *name;
+	int has_arg;
+	int *flag;
+	int val;
+};
+
+int getopt_long(int, char *const *, const char *, const struct option *, int *);
+int getopt_long_only(int, char *const *, const char *, const struct option *, int *);
+
+#define no_argument        0
+#define required_argument  1
+#define optional_argument  2
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

index.c

@@ -1,6 +1,10 @@
 #include <stdlib.h>
 #include <assert.h>
+#if defined(WIN32) || defined(_WIN32)
+#include <io.h> // for open(2)
+#else
 #include <unistd.h>
+#endif
 #include <fcntl.h>
 #include <stdio.h>
 #include "kthread.h"

kalloc.c

@@ -46,7 +46,7 @@ static size_t *morecore(kmem_t *km, size_t nu)
 	up = (size_t*)malloc(rnu * sizeof(size_t));
 	if (!up) { /* fail to allocate memory */
 		km_stat(km);
-		fprintf(stderr, "[morecore] %lu bytes requested but not available.\n", rnu * sizeof(size_t));
+		fprintf(stderr, "[morecore] %lu bytes requested but not available.\n", (unsigned long)rnu * sizeof(size_t));
 		exit(1);
 	}
 	/* put the pointer in km->list_head */
@@ -210,5 +210,5 @@ void km_stat(const void *_km)
 	--n_blocks;
 	frag = 1.0/1024.0 * n_units * sizeof(size_t) / n_blocks;
 	fprintf(stderr, "[kr_stat] tot=%lu, free=%lu, n_block=%u, max_block=%lu, frag_len=%.3fK\n",
-			km->total_allocated, n_units * sizeof(size_t), n_blocks, max_block * sizeof(size_t), frag);
+			(unsigned long)km->total_allocated, (unsigned long)n_units * sizeof(size_t), n_blocks, (unsigned long)max_block * sizeof(size_t), frag);
 }

kdq.h

@@ -3,11 +3,12 @@
 
 #include <stdlib.h>
 #include <string.h>
+#include <stdint.h>
 #include "kalloc.h"
 
 #define __KDQ_TYPE(type) \
 	typedef struct { \
-		size_t front:58, bits:6, count, mask; \
+		uint64_t front:58, bits:6, count, mask; \
 		type *a; \
 		void *km; \
 	} kdq_##type##_t;

ksort.h

@@ -30,6 +30,7 @@
 
 #include <stdlib.h>
 #include <string.h>
+#include <assert.h>
 
 typedef struct {
 	void *left, *right;
@@ -78,6 +79,7 @@ typedef const char *ksstr_t;
 #define KSORT_INIT_STR KSORT_INIT(str, ksstr_t, ks_lt_str)
 
 #define RS_MIN_SIZE 64
+#define RS_MAX_BITS 8
 
 #define KRADIX_SORT_INIT(name, rstype_t, rskey, sizeof_key) \
 	typedef struct { \
@@ -98,7 +100,8 @@ typedef const char *ksstr_t;
 	{ \
 		rstype_t *i; \
 		int size = 1<<n_bits, m = size - 1; \
-		rsbucket_##name##_t *k, b[size], *be = b + size; \
+		rsbucket_##name##_t *k, b[1<<RS_MAX_BITS], *be = b + size; \
+		assert(n_bits <= RS_MAX_BITS); \
 		for (k = b; k != be; ++k) k->b = k->e = beg; \
 		for (i = beg; i != end; ++i) ++b[rskey(*i)>>s&m].e; \
 		for (k = b + 1; k != be; ++k) \
@@ -127,7 +130,7 @@ typedef const char *ksstr_t;
 	void radix_sort_##name(rstype_t *beg, rstype_t *end) \
 	{ \
 		if (end - beg <= RS_MIN_SIZE) rs_insertsort_##name(beg, end); \
-		else rs_sort_##name(beg, end, 8, sizeof_key * 8 - 8); \
+		else rs_sort_##name(beg, end, RS_MAX_BITS, (sizeof_key - 1) * RS_MAX_BITS); \
 	}
 
 #endif

ksw2.h

@@ -169,5 +169,4 @@ static inline int ksw_apply_zdrop(ksw_extz_t *ez, int is_rot, int32_t H, int a,
 	}
 	return 0;
 }
-
 #endif

ksw2_dispatch.c

@@ -0,0 +1,97 @@
+#ifdef KSW_CPU_DISPATCH
+#include <stdlib.h>
+#include "ksw2.h"
+
+#define SIMD_SSE     0x1
+#define SIMD_SSE2    0x2
+#define SIMD_SSE3    0x4
+#define SIMD_SSSE3   0x8
+#define SIMD_SSE4_1  0x10
+#define SIMD_SSE4_2  0x20
+#define SIMD_AVX     0x40
+#define SIMD_AVX2    0x80
+#define SIMD_AVX512F 0x100
+
+#ifndef _MSC_VER
+// adapted from https://github.com/01org/linux-sgx/blob/master/common/inc/internal/linux/cpuid_gnu.h
+void __cpuidex(int cpuid[4], int func_id, int subfunc_id)
+{
+#if defined(__x86_64__)
+	asm volatile ("cpuid"
+			: "=a" (cpuid[0]), "=b" (cpuid[1]), "=c" (cpuid[2]), "=d" (cpuid[3])
+			: "0" (func_id), "2" (subfunc_id));
+#else // on 32bit, ebx can NOT be used as PIC code
+	asm volatile ("xchgl %%ebx, %1; cpuid; xchgl %%ebx, %1"
+			: "=a" (cpuid[0]), "=r" (cpuid[1]), "=c" (cpuid[2]), "=d" (cpuid[3])
+			: "0" (func_id), "2" (subfunc_id));
+#endif
+}
+#endif
+
+int x86_simd(void)
+{
+	int flag = 0, cpuid[4], max_id;
+	__cpuidex(cpuid, 0, 0);
+	max_id = cpuid[0];
+	if (max_id == 0) return 0;
+	__cpuidex(cpuid, 1, 0);
+	if (cpuid[3]>>25&1) flag |= SIMD_SSE;
+	if (cpuid[3]>>26&1) flag |= SIMD_SSE2;
+	if (cpuid[2]>>0 &1) flag |= SIMD_SSE3;
+	if (cpuid[2]>>9 &1) flag |= SIMD_SSSE3;
+	if (cpuid[2]>>19&1) flag |= SIMD_SSE4_1;
+	if (cpuid[2]>>20&1) flag |= SIMD_SSE4_2;
+	if (cpuid[2]>>28&1) flag |= SIMD_AVX;
+	if (max_id >= 7) {
+		__cpuidex(cpuid, 7, 0);
+		if (cpuid[1]>>5 &1) flag |= SIMD_AVX2;
+		if (cpuid[1]>>16&1) flag |= SIMD_AVX512F;
+	}
+	return flag;
+}
+
+void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez)
+{
+	extern void ksw_extz2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez);
+	extern void ksw_extz2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez);
+	unsigned simd;
+	simd = x86_simd();
+	if (simd & SIMD_SSE4_1)
+		ksw_extz2_sse41(km, qlen, query, tlen, target, m, mat, q, e, w, zdrop, flag, ez);
+	else if (simd & SIMD_SSE2)
+		ksw_extz2_sse2(km, qlen, query, tlen, target, m, mat, q, e, w, zdrop, flag, ez);
+	else abort();
+}
+
+void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int flag, ksw_extz_t *ez)
+{
+	extern void ksw_extd2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int flag, ksw_extz_t *ez);
+	extern void ksw_extd2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int flag, ksw_extz_t *ez);
+	unsigned simd;
+	simd = x86_simd();
+	if (simd & SIMD_SSE4_1)
+		ksw_extd2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, flag, ez);
+	else if (simd & SIMD_SSE2)
+		ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, flag, ez);
+	else abort();
+}
+
+void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez)
+{
+	extern void ksw_exts2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
+	extern void ksw_exts2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
+	unsigned simd;
+	simd = x86_simd();
+	if (simd & SIMD_SSE4_1)
+		ksw_exts2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
+	else if (simd & SIMD_SSE2)
+		ksw_exts2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
+	else abort();
+}
+#endif

ksw2_extd2_sse.c

@@ -6,12 +6,26 @@
 #ifdef __SSE2__
 #include <emmintrin.h>
 
+#ifdef KSW_SSE2_ONLY
+#undef __SSE4_1__
+#endif
+
 #ifdef __SSE4_1__
 #include <smmintrin.h>
 #endif
 
+#ifdef KSW_CPU_DISPATCH
+#ifdef __SSE4_1__
+void ksw_extd2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int flag, ksw_extz_t *ez)
+#else
+void ksw_extd2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int flag, ksw_extz_t *ez)
+#endif
+#else
 void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
 				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int flag, ksw_extz_t *ez)
+#endif // ~KSW_CPU_DISPATCH
 {
 #define __dp_code_block1 \
 	z = _mm_load_si128(&s[t]); \

ksw2_exts2_sse.c

@@ -6,12 +6,26 @@
 #ifdef __SSE2__
 #include <emmintrin.h>
 
+#ifdef KSW_SSE2_ONLY
+#undef __SSE4_1__
+#endif
+
 #ifdef __SSE4_1__
 #include <smmintrin.h>
 #endif
 
+#ifdef KSW_CPU_DISPATCH
+#ifdef __SSE4_1__
+void ksw_exts2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez)
+#else
+void ksw_exts2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez)
+#endif
+#else
 void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
 				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez)
+#endif // ~KSW_CPU_DISPATCH
 {
 #define __dp_code_block1 \
 	z = _mm_load_si128(&s[t]); \

ksw2_extz2_sse.c

@@ -5,11 +5,23 @@
 #ifdef __SSE2__
 #include <emmintrin.h>
 
+#ifdef KSW_SSE2_ONLY
+#undef __SSE4_1__
+#endif
+
 #ifdef __SSE4_1__
 #include <smmintrin.h>
 #endif
 
+#ifdef KSW_CPU_DISPATCH
+#ifdef __SSE4_1__
+void ksw_extz2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez)
+#else
+void ksw_extz2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez)
+#endif
+#else
 void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez)
+#endif // ~KSW_CPU_DISPATCH
 {
 #define __dp_code_block1 \
 	z = _mm_add_epi8(_mm_load_si128(&s[t]), qe2_); \

kthread.c

@@ -1,6 +1,11 @@
 #include <pthread.h>
 #include <stdlib.h>
 #include <limits.h>
+#include <stdint.h>
+
+#if (defined(WIN32) || defined(_WIN32)) && defined(_MSC_VER)
+#define __sync_fetch_and_add(ptr, addend)     _InterlockedExchangeAdd((void*)ptr, addend)
+#endif
 
 /************
  * kt_for() *
@@ -52,12 +57,13 @@ void kt_for(int n_threads, void (*func)(void*,long,int), void *data, long n)
 		kt_for_t t;
 		pthread_t *tid;
 		t.func = func, t.data = data, t.n_threads = n_threads, t.n = n;
-		t.w = (ktf_worker_t*)alloca(n_threads * sizeof(ktf_worker_t));
-		tid = (pthread_t*)alloca(n_threads * sizeof(pthread_t));
+		t.w = (ktf_worker_t*)calloc(n_threads, sizeof(ktf_worker_t));
+		tid = (pthread_t*)calloc(n_threads, sizeof(pthread_t));
 		for (i = 0; i < n_threads; ++i)
 			t.w[i].t = &t, t.w[i].i = i;
 		for (i = 0; i < n_threads; ++i) pthread_create(&tid[i], 0, ktf_worker, &t.w[i]);
 		for (i = 0; i < n_threads; ++i) pthread_join(tid[i], 0);
+		free(tid); free(t.w);
 	} else {
 		long j;
 		for (j = 0; j < n; ++j) func(data, j, 0);
@@ -135,16 +141,17 @@ void kt_pipeline(int n_threads, void *(*func)(void*, int, void*), void *shared_d
 	pthread_mutex_init(&aux.mutex, 0);
 	pthread_cond_init(&aux.cv, 0);
 
-	aux.workers = (ktp_worker_t*)alloca(n_threads * sizeof(ktp_worker_t));
+	aux.workers = (ktp_worker_t*)calloc(n_threads, sizeof(ktp_worker_t));
 	for (i = 0; i < n_threads; ++i) {
 		ktp_worker_t *w = &aux.workers[i];
 		w->step = 0; w->pl = &aux; w->data = 0;
 		w->index = aux.index++;
 	}
 
-	tid = (pthread_t*)alloca(n_threads * sizeof(pthread_t));
+	tid = (pthread_t*)calloc(n_threads, sizeof(pthread_t));
 	for (i = 0; i < n_threads; ++i) pthread_create(&tid[i], 0, ktp_worker, &aux.workers[i]);
 	for (i = 0; i < n_threads; ++i) pthread_join(tid[i], 0);
+	free(tid); free(aux.workers);
 
 	pthread_mutex_destroy(&aux.mutex);
 	pthread_cond_destroy(&aux.cv);

main.c

@@ -1,24 +1,26 @@
-#include <getopt.h>
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
-#include <sys/resource.h>
-#include <sys/time.h>
 #include "bseq.h"
 #include "minimap.h"
 #include "mmpriv.h"
+#include "getopt.h"
 
-#define MM_VERSION "2.1-r311"
+#define MM_VERSION "2.1.1-r341"
 
+#ifdef __linux__
+#include <sys/resource.h>
+#include <sys/time.h>
 void liftrlimit()
 {
-#ifdef __linux__
 	struct rlimit r;
 	getrlimit(RLIMIT_AS, &r);
 	r.rlim_cur = r.rlim_max;
 	setrlimit(RLIMIT_AS, &r);
-#endif
 }
+#else
+void liftrlimit() {}
+#endif
 
 static struct option long_options[] = {
 	{ "bucket-bits",    required_argument, 0, 0 },
@@ -240,6 +242,7 @@ int main(int argc, char *argv[])
 		mm_idx_t *mi;
 		if (fpr) {
 			mi = mm_idx_load(fpr);
+			if (mi == 0) break;
 			if (idx_par_set && mm_verbose >= 2 && (mi->k != k || mi->w != w || mi->is_hpc != is_hpc))
 				fprintf(stderr, "[WARNING] \033[1;31mIndexing parameters on the command line (-k/-w/-H) overridden by parameters in the prebuilt index.\033[0m\n");
 		} else {

map.c

@@ -382,7 +382,11 @@ int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int
 	pipeline_t pl;
 	memset(&pl, 0, sizeof(pipeline_t));
 	pl.fp = mm_bseq_open(fn);
-	if (pl.fp == 0) return -1;
+	if (pl.fp == 0) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "ERROR: failed to open file '%s'\n", fn);
+		return -1;
+	}
 	pl.opt = opt, pl.mi = idx;
 	pl.n_threads = n_threads, pl.mini_batch_size = mini_batch_size;
 	if ((opt->flag & MM_F_OUT_SAM) && !(opt->flag & MM_F_NO_SAM_SQ))

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "25 August 2017" "minimap2-2.1-r311" "Bioinformatics tools"
+.TH minimap2 1 "6 September 2017" "minimap2-2.1.1-r341" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences

misc.c

@@ -1,17 +1,102 @@
-#include <sys/resource.h>
-#include <sys/time.h>
 #include "minimap.h"
 
 int mm_verbose = 3;
 int mm_dbg_flag = 0;
 double mm_realtime0;
 
+#if defined(WIN32) || defined(_WIN32)
+#include <windows.h>
+
+struct timezone
+{
+  __int32  tz_minuteswest; /* minutes W of Greenwich */
+  int  tz_dsttime;     /* type of dst correction */
+};
+
+/*
+ * gettimeofday.c
+ *    Win32 gettimeofday() replacement
+ *    taken from PostgreSQL, according to
+ *    https://stackoverflow.com/questions/1676036/what-should-i-use-to-replace-gettimeofday-on-windows
+ *
+ * src/port/gettimeofday.c
+ *
+ * Copyright (c) 2003 SRA, Inc.
+ * Copyright (c) 2003 SKC, Inc.
+ *
+ * Permission to use, copy, modify, and distribute this software and
+ * its documentation for any purpose, without fee, and without a
+ * written agreement is hereby granted, provided that the above
+ * copyright notice and this paragraph and the following two
+ * paragraphs appear in all copies.
+ *
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE TO ANY PARTY FOR DIRECT,
+ * INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING
+ * LOST PROFITS, ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS
+ * DOCUMENTATION, EVEN IF THE UNIVERSITY OF CALIFORNIA HAS BEEN ADVISED
+ * OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ * THE AUTHOR SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ * A PARTICULAR PURPOSE.  THE SOFTWARE PROVIDED HEREUNDER IS ON AN "AS
+ * IS" BASIS, AND THE AUTHOR HAS NO OBLIGATIONS TO PROVIDE MAINTENANCE,
+ * SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
+ */
+
+/* FILETIME of Jan 1 1970 00:00:00. */
+static const unsigned __int64 epoch = ((unsigned __int64) 116444736000000000ULL);
+
+/*
+ * timezone information is stored outside the kernel so tzp isn't used anymore.
+ *
+ * Note: this function is not for Win32 high precision timing purpose. See
+ * elapsed_time().
+ */
+int gettimeofday(struct timeval * tp, struct timezone *tzp)
+{
+    FILETIME    file_time;
+    SYSTEMTIME  system_time;
+    ULARGE_INTEGER ularge;
+
+    GetSystemTime(&system_time);
+    SystemTimeToFileTime(&system_time, &file_time);
+    ularge.LowPart = file_time.dwLowDateTime;
+    ularge.HighPart = file_time.dwHighDateTime;
+
+    tp->tv_sec = (long) ((ularge.QuadPart - epoch) / 10000000L);
+    tp->tv_usec = (long) (system_time.wMilliseconds * 1000);
+
+    return 0;
+}
+
+// taken from https://stackoverflow.com/questions/5272470/c-get-cpu-usage-on-linux-and-windows
+double cputime()
+{
+	HANDLE hProcess = GetCurrentProcess();
+	FILETIME ftCreation, ftExit, ftKernel, ftUser;
+	SYSTEMTIME stKernel;
+	SYSTEMTIME stUser;
+
+	GetProcessTimes(hProcess, &ftCreation, &ftExit, &ftKernel, &ftUser);
+	FileTimeToSystemTime(&ftKernel, &stKernel);
+	FileTimeToSystemTime(&ftUser, &stUser);
+
+	double kernelModeTime = ((stKernel.wHour * 60.) + stKernel.wMinute * 60.) + stKernel.wSecond * 1. + stKernel.wMilliseconds / 1000.;
+	double userModeTime = ((stUser.wHour * 60.) + stUser.wMinute * 60.) + stUser.wSecond * 1. + stUser.wMilliseconds / 1000.;
+
+	return kernelModeTime + userModeTime;
+}
+#else
+#include <sys/resource.h>
+#include <sys/time.h>
+
 double cputime()
 {
 	struct rusage r;
 	getrusage(RUSAGE_SELF, &r);
 	return r.ru_utime.tv_sec + r.ru_stime.tv_sec + 1e-6 * (r.ru_utime.tv_usec + r.ru_stime.tv_usec);
 }
+#endif /* WIN32 || _WIN32 */
 
 double realtime()
 {

sdust.c

@@ -176,7 +176,7 @@ uint64_t *sdust(void *km, const uint8_t *seq, int l_seq, int T, int W, int *n)
 #ifdef _SDUST_MAIN
 #include <zlib.h>
 #include <stdio.h>
-#include <unistd.h>
+#include "getopt.h"
 #include "kseq.h"
 KSEQ_INIT(gzFile, gzread)
 

sketch.c

@@ -77,11 +77,10 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i
 {
 	uint64_t shift1 = 2 * (k - 1), mask = (1ULL<<2*k) - 1, kmer[2] = {0,0};
 	int i, j, l, buf_pos, min_pos, kmer_span = 0;
-	mm128_t *buf, min = { UINT64_MAX, UINT64_MAX };
+	mm128_t buf[256], min = { UINT64_MAX, UINT64_MAX };
 	tiny_queue_t tq;
 
-	assert(len > 0 && w > 0 && k > 0 && k <= 28); // 56 bits for k-mer; could use long k-mers, but 28 enough in practice
-	buf = (mm128_t*)alloca(w * 16);
+	assert(len > 0 && (w > 0 && w < 256) && (k > 0 && k <= 28)); // 56 bits for k-mer; could use long k-mers, but 28 enough in practice
 	memset(buf, 0xff, w * 16);
 	memset(&tq, 0, sizeof(tiny_queue_t));
 	kv_resize(mm128_t, km, *p, p->n + len/w);

tex/minimap2.tex

@@ -177,8 +177,8 @@ based on Eq.~(\ref{eq:ae86}) can achieve 16-way parallelization for short
 sequences, but only 4-way parallelization when the peak alignment score reaches
 32767. Long sequence alignment may exceed this threshold. Inspired by
 \citet{Wu:1996aa} and the following work, \citet{Suzuki:2016} proposed a
-difference-based formulation that lifted this limitation. In case of 2-piece
-gap cost, define
+difference-based formulation that lifted this limitation.
+In case of 2-piece gap cost, define
 \[
 \left\{\begin{array}{ll}
 u_{ij}\triangleq H_{ij}-H_{i-1,j} & v_{ij}\triangleq H_{ij}-H_{i,j-1} \\
@@ -199,9 +199,10 @@ y_{ij}&=&\max\{0,y_{i,j-1}+u_{i,j-1}-z_{ij}+q\}-q-e\\
 \tilde{y}_{ij}&=&\max\{0,\tilde{y}_{i,j-1}+u_{i,j-1}-z_{ij}+\tilde{q}\}-\tilde{q}-\tilde{e}
 \end{array}\right.
 \end{equation}
-where $z_{ij}$ is a temporary variable that does not need to be stored.  An
-important property of Eq.~(\ref{eq:suzuki}) is that all values are bounded. To
-see that,
+where $z_{ij}$ is a temporary variable that does not need to be stored.
+
+An important property of Eq.~(\ref{eq:suzuki}) is that all values are bounded
+by scoring parameters. To see that,
 \[
 x_{ij}=E_{i+1,j}-H_{ij}=\max\{-q,E_{ij}-H_{ij}\}-e
 \]
@@ -245,8 +246,8 @@ each other. This allows us to fully vectorize the computation of all cells on
 the same anti-diagonal in one inner loop. It also simplifies banded alignment,
 which would be difficult with striped vectorization~\citep{Farrar:2007hs}.
 
-On the condition that $q+e<\tilde{q}+\tilde{e}$ and $e>\tilde{e}$, the boundary
-condition of the equation above is
+On the condition that $q+e<\tilde{q}+\tilde{e}$ and $e>\tilde{e}$, the initial
+values in the diagonal-antidiagonal formuation is
 \[
 \left\{\begin{array}{l}
 x_{r-1,-1}=y_{r-1,r}=-q-e\\
@@ -263,7 +264,7 @@ r\cdot(e-\tilde{e})-(\tilde{q}-q)-\tilde{e} & (r=\lceil\frac{\tilde{q}-q}{e-\til
 -\tilde{e} & (r>\lceil\frac{\tilde{q}-q}{e-\tilde{e}}-1\rceil)
 \end{array}\right.
 \]
-These can be derived from the initial conditions of Eq.~(\ref{eq:ae86}).
+These can be derived from the initial values for Eq.~(\ref{eq:ae86}).
 
 In practice, our 16-way vectorized implementation of global alignment is three
 times as fast as Parasail's 4-way vectorization~\citep{Daily:2016aa}.  Without
@@ -285,9 +286,9 @@ $j'<j$, such that
 S(i',j')-S(i,j)>Z+e\cdot|(i-i')-(j-j')|
 \]
 where $e$ is the gap extension cost and $Z$ is an arbitrary threshold.
-This strategy is similar to X-drop employed in BLAST~\citep{Altschul:1997vn}.
-However, unlike X-drop, it would not break the alignment in the presence of a
-single long gap.
+This strategy is first used in BWA-MEM. It is similar to X-drop employed in
+BLAST~\citep{Altschul:1997vn}, but unlike X-drop, it would not break the
+alignment in the presence of a single long gap.
 
 When minimap2 breaks a global alignment between two anchors, it performs local
 alignment between the two subsequences involved in the global alignment, but
@@ -326,7 +327,7 @@ F_{i,j+1}= \max\{H_{ij}-q,F_{ij}\}-e\\
 \end{array}\right.
 \end{equation}
 Let $T$ be the reference sequence. $d(i)$ is the cost of a non-canonical donor
-site, which takes 0 if $T[i+1,i+2]={\tt GT}$, or a postive number $p$
+site, which takes 0 if $T[i+1,i+2]={\tt GT}$, or a positive number $p$
 otherwise. Similarly, $a(i)$ is the cost of a non-canonical acceptor site, which
 takes 0 if $T[i-1,i]={\tt AG}$, or $p$ otherwise. Eq.~(\ref{eq:splice}) is
 almost equivalent to the equation used by EXALIN~\citep{Zhang:2006aa} except
@@ -361,18 +362,18 @@ alignment.
 \centering
 \includegraphics[width=.5\textwidth]{roc-color.pdf}
 \caption{Evaluation on simulated SMRT reads aligned against human genome
-GRCh38. (a) ROC-like curve. Alignments are sorted by mapping quality in the
-descending order. For each mapping quality threshold, the fraction of
-alignments with mapping quality above the threshold and their error rate
-are plotted. (b) Accumulative mapping error rate as a function of mapping
-quality. 33,088 $\ge$1000bp reads were simulated using pbsim~\citep{Ono:2013aa}
+GRCh38.  33,088 $\ge$1000bp reads were simulated using pbsim~\citep{Ono:2013aa}
 with error profile sampled from file `m131017\_060208\_42213\_*.1.*' downloaded
 at \href{http://bit.ly/chm1p5c3}{http://bit.ly/chm1p5c3}. The N50 read length
 is 11,628. A read is considered correctly mapped if the true position overlaps
 with the best mapping position by 10\% of the read length. All aligners were
-run under the default setting for SMRT reads. Kart outputted all alignments at
-mapping quality 60, so is not shown in the figure. It mapped nearly all reads
-with 4.1\% of alignments being wrong, less accurate than others.}\label{fig:eval}
+run under the default setting for SMRT reads.  (a) ROC-like curve. Alignments
+are sorted by mapping quality in the descending order. For each mapping quality
+threshold, the fraction of alignments with mapping quality above the threshold
+and their error rate are plotted.  Kart outputted all alignments at mapping
+quality 60, so is not shown in the figure. It mapped nearly all reads with
+4.1\% of alignments being wrong, less accurate than others.  (b) Accumulative
+mapping error rate as a function of mapping quality.}\label{fig:eval}
 \end{figure}
 
 As a sanity check, we evaluated minimap2 on simulated human reads along with
@@ -383,7 +384,7 @@ Kart~(v2.2.5; \citealp{Lin:2017aa}),
 minialign~(v0.5.3; \citealp{Suzuki:2016}) and
 NGMLR~(v0.2.5; \citealp{Sedlazeck169557}). We excluded rHAT~\citep{Liu:2016ab}
 and LAMSA~\citep{Liu:2017aa} because they either
-crashed or produced malformatted output. In this evaluation, Minimap2 has
+crashed or produced malformatted output. In this evaluation, minimap2 has
 higher power to distinguish unique and repetitive hits, and achieves overall
 higher mapping accuracy (Fig.~\ref{fig:eval}a). It is still the most accurate
 even if we skip DP-based alignment (data not shown), confirming chaining alone
@@ -500,8 +501,8 @@ necessary to justify the use of minimap2 for such applications.
 We owe a debt of gratitude to Hajime Suzuki for releasing his 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 fix various
-issues.
+grateful to early minimap2 testers who have greatly helped to suggest features
+and to fix various issues.
 
 \bibliography{minimap2}