Release minimap2-2.2 (r409)

can't get cython working...

.gitignore

@@ -1,4 +1,8 @@
+.cproject
+.project
 .*.swp
 *.a
 *.o
 *.dSYM
+minimap2
+mappy.c

.travis.yml

@@ -1,5 +1,24 @@
-language: c
+matrix:
-compiler:
+  include:
-  - gcc
+    - language: c
-  - clang
+      compiler: gcc
-script: make
+      script: make
+    - language: c
+      compiler: clang
+      script: make
+    - language: python
+      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
+      script: python setup.py build_ext
+    - language: python
+      python: "3.6"
+      before_install: pip install cython
+      script: python setup.py build_ext

MANIFEST.in

@@ -0,0 +1,11 @@
+include *.h
+include Makefile
+include ksw2_dispatch.c
+include getopt.c
+include main.c
+include README.md
+include python/mappy.c
+include python/cmappy.h
+include python/cmappy.pxd
+include python/mappy.pyx
+include python/README.rst

Makefile

@@ -1,15 +1,15 @@
-CC=			gcc
 CFLAGS=		-g -Wall -O2 -Wc++-compat
 CPPFLAGS=	-DHAVE_KALLOC
-INCLUDES=	-I.
+INCLUDES=
-OBJS=		kthread.o kalloc.o ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_exts2_sse.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 ksw2_ll_sse.o
-			misc.o bseq.o sketch.o sdust.o index.o chain.o align.o hit.o map.o format.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 +21,8 @@ all:$(PROG)
 
 extra:all $(PROG_EXTRA)
 
-minimap2:main.o libminimap2.a
+minimap2:main.o getopt.o libminimap2.a
-		$(CC) $(CFLAGS) $< -o $@ -L. -lminimap2 $(LIBS)
+		$(CC) $(CFLAGS) main.o getopt.o -o $@ -L. -lminimap2 $(LIBS)
 
 minimap2-lite:example.o libminimap2.a
 		$(CC) $(CFLAGS) $< -o $@ -L. -lminimap2 $(LIBS)
@@ -30,11 +30,32 @@ 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
+sdust:sdust.c getopt.o kalloc.o kalloc.h kdq.h kvec.h kseq.h sdust.h
-		$(CC) -D_SDUST_MAIN $(CFLAGS) $< kalloc.o -o $@ -lz
+		$(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*
+		rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist mappy.so mappy.c python/mappy.c mappy.egg*
 
 depend:
 		(LC_ALL=C; export LC_ALL; makedepend -Y -- $(CFLAGS) $(CPPFLAGS) -- *.c)
@@ -46,6 +67,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 +75,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,55 @@
+Release 2.2-r409 (17 September 2017)
+------------------------------------
+
+This is a feature release. It improves single-end short-read alignment and
+comes with Python bindings. Detailed changes include:
+
+ * Added the **sr** preset for single-end short-read alignment. In this mode,
+   minimap2 runs faster than BWA-MEM, but is slightly less accurate on
+   simulated data sets. Paired-end alignment is not supported as of now.
+
+ * Improved mapping quality estimate with more accurate identification of
+   repetitive hits. This mainly helps short-read alignment.
+
+ * Implemented **mappy**, a Python binding for minimap2, which is available
+   from PyPI and can be installed with `pip install --user mappy`. Python users
+   can perform read alignment without the minimap2 executable.
+
+ * Restructured the indexing APIs and documented key minimap2 APIs in the
+   header file minimap.h. Updated example.c with the new APIs. Old APIs still
+   work but may become deprecated in future.
+
+This release may output alignments different from the previous version, though
+the overall alignment statistics, such as the number of aligned bases and long
+gaps, remain close.
+
+(2.2: 17 September 2017, r409)
+
+
+
+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

@@ -1,3 +1,8 @@
+[![Release](https://img.shields.io/badge/Release-v2.1.1-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)
+[![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)
 ## Getting Started
 ```sh
@@ -10,6 +15,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,14 @@ 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 < tlen; ++i) fputc("ACGTN"[tseq[i]], stderr);
-		for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr); fputc('\n', stderr);
+		fputc('\n', stderr);
+		for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr);
+		fputc('\n', stderr);
+	}
+	if (opt->flag & MM_F_APPROX_EXT) {
+		flag |= KSW_EZ_APPROX_MAX;
+		if (flag & KSW_EZ_EXTZ_ONLY) flag |= KSW_EZ_APPROX_DROP;
 	}
 	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);

bseq.c

@@ -5,7 +5,7 @@
 #include <assert.h>
 #include "bseq.h"
 #include "kseq.h"
-KSEQ_INIT(gzFile, gzread)
+KSEQ_INIT2(, gzFile, gzread)
 
 struct mm_bseq_file_s {
 	gzFile fp;

chain.c

@@ -45,20 +45,21 @@ int mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cn
 		while (st < i && ri - a[st].x > max_dist_x) ++st;
 		for (j = i - 1; j >= st; --j) {
 			int64_t dr = ri - a[j].x;
-			int32_t dq = qi - (int32_t)a[j].y, dd, sc;
+			int32_t dq = qi - (int32_t)a[j].y, dd, sc, log_dd;
 			if (dr == 0 || dq <= 0 || dq > max_dist_y) continue;
 			dd = dr > dq? dr - dq : dq - dr;
 			if (dd > bw) continue;
 			max_f_past = max_f_past > f[j]? max_f_past : f[j];
 			min_d = dq < dr? dq : dr;
 			sc = min_d > q_span? q_span : dq < dr? dq : dr;
+			log_dd = dd? ilog2_32(dd) : 0;
 			if (is_cdna) {
 				int c_log, c_lin;
 				c_lin = (int)(dd * .01 * avg_qspan);
-				c_log = ilog2_32(dd);
+				c_log = log_dd;
 				if (dr > dq) sc -= c_lin < c_log? c_lin : c_log;
 				else sc -= c_lin + (c_log>>1);
-			} else sc -= (int)(dd * .01 * avg_qspan) + (ilog2_32(dd)>>1);
+			} else sc -= (int)(dd * .01 * avg_qspan) + (log_dd>>1);
 			sc += f[j];
 			if (sc > max_f) {
 				max_f = sc, max_j = j;

example.c

@@ -11,7 +11,13 @@ KSEQ_INIT(gzFile, gzread)
 
 int main(int argc, char *argv[])
 {
+	mm_idxopt_t iopt;
+	mm_mapopt_t mopt;
+	int n_threads = 3;
+
 	mm_verbose = 2; // disable message output to stderr
+	mm_set_opt(0, &iopt, &mopt);
+	mopt.flag |= MM_F_CIGAR; // perform alignment
 
 	if (argc < 3) {
 		fprintf(stderr, "Usage: minimap2-lite <target.fa> <query.fa>\n");
@@ -23,39 +29,34 @@ int main(int argc, char *argv[])
 	assert(f);
 	kseq_t *ks = kseq_init(f);
 
-	// create index for target; we are creating one index for all target sequence
+	// open index reader
-	int n_threads = 4, w = 10, k = 15, is_hpc = 0;
+	mm_idx_reader_t *r = mm_idx_reader_open(argv[1], &iopt, 0);
-	mm_idx_t *mi = mm_idx_build(argv[1], w, k, is_hpc, n_threads);
+	mm_idx_t *mi;
-	assert(mi);
+	while ((mi = mm_idx_reader_read(r, n_threads)) != 0) { // traverse each part of the index
-
+		mm_mapopt_update(&mopt, mi); // this sets the maximum minimizer occurrence; TODO: set a better default in mm_mapopt_init()!
-	// mapping
+		mm_tbuf_t *tbuf = mm_tbuf_init(); // thread buffer; for multi-threading, allocate one tbuf for each thread
-	mm_mapopt_t opt;
+		while (kseq_read(ks) >= 0) { // each kseq_read() call reads one query sequence
-	mm_mapopt_init(&opt); // initialize mapping parameters
+			mm_reg1_t *reg;
-	mm_mapopt_update(&opt, mi); // this sets the maximum minimizer occurrence; TODO: set a better default in mm_mapopt_init()!
+			int j, i, n_reg;
-	opt.flag |= MM_F_CIGAR; // perform alignment
+			reg = mm_map(mi, ks->seq.l, ks->seq.s, &n_reg, tbuf, &mopt, 0); // get all hits for the query
-	mm_tbuf_t *tbuf = mm_tbuf_init(); // thread buffer; for multi-threading, allocate one tbuf for each thread
+			for (j = 0; j < n_reg; ++j) { // traverse hits and print them out
-	while (kseq_read(ks) >= 0) { // each kseq_read() call reads one query sequence
+				mm_reg1_t *r = &reg[j];
-		const mm_reg1_t *reg;
+				assert(r->p); // with MM_F_CIGAR, this should not be NULL
-		int j, i, n_reg;
+				printf("%s\t%d\t%d\t%d\t%c\t", ks->name.s, ks->seq.l, r->qs, r->qe, "+-"[r->rev]);
-		// get all hits for the query
+				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,
-		reg = mm_map(mi, ks->seq.l, ks->seq.s, &n_reg, tbuf, &opt, 0);
+						r->p->blen - r->p->n_ambi - r->p->n_diff, r->p->blen, r->mapq);
-		// traverse hits and print them out
+				for (i = 0; i < r->p->n_cigar; ++i) // IMPORTANT: this gives the CIGAR in the aligned regions. NO soft/hard clippings!
-		for (j = 0; j < n_reg; ++j) {
+					printf("%d%c", r->p->cigar[i]>>4, "MIDSHN"[r->p->cigar[i]&0xf]);
-			const mm_reg1_t *r = &reg[j];
+				putchar('\n');
-			assert(r->p); // with MM_F_CIGAR, this should not be NULL
+				free(r->p);
-			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,
+			free(reg);
-				   r->p->blen - r->p->n_ambi - r->p->n_diff, r->p->blen, r->mapq);
-			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');
 		}
+		mm_tbuf_destroy(tbuf);
+		mm_idx_destroy(mi);
 	}
-	mm_tbuf_destroy(tbuf);
+	mm_idx_reader_close(r); // close the index reader
-
+	kseq_destroy(ks); // close the query file
-	// deallocate index and close the query file
-	mm_idx_destroy(mi);
-	kseq_destroy(ks);
 	gzclose(f);
 	return 0;
 }

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 = -1;
+		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

hit.c

@@ -87,7 +87,7 @@ void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a)
 	r->split |= 1, r2->split |= 2;
 }
 
-void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r) // and compute mm_reg1_t::subsc
+void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff) // and compute mm_reg1_t::subsc
 {
 	int i, j, k, *w;
 	if (n <= 0) return;
@@ -103,14 +103,19 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r) // and compu
 			int min = ej - sj < ei - si? ej - sj : ei - si;
 			int ol = si < sj? (ei < sj? 0 : ei < ej? ei - sj : ej - sj) : (ej < si? 0 : ej < ei? ej - si : ei - si);
 			if (ol > mask_level * min) {
+				int cnt_sub = 0;
 				ri->parent = rp->parent;
 				rp->subsc = rp->subsc > ri->score? rp->subsc : ri->score;
-				if (rp->p && ri->p)
+				if (ri->cnt >= rp->cnt) cnt_sub = 1;
+				if (rp->p && ri->p) {
 					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;
+				}
+				if (cnt_sub) ++rp->n_sub;
 				break;
 			}
 		}
-		if (j == k) w[k++] = i, ri->parent = i;
+		if (j == k) w[k++] = i, ri->parent = i, ri->n_sub = 0;
 	}
 	kfree(km, w);
 }
@@ -288,9 +293,9 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
 	}
 }
 
-void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc)
+void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len)
 {
-	static const float q_coef = 30.0f;
+	static const float q_coef = 40.0f;
 	int i;
 	for (i = 0; i < n_regs; ++i) {
 		mm_reg1_t *r = &regs[i];
@@ -298,12 +303,21 @@ void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc)
 			r->mapq = 0;
 		} else if (r->parent == r->id) {
 			int mapq, subsc;
-			float pen_cm = r->cnt >= 10? 1.0f : 0.1f * r->cnt;
+			float pen_s1 = r->score > 100? 1.0f : 0.01f * r->score;
+			float pen_cm = r->cnt > 10? 1.0f : 0.1f * r->cnt;
+			if (r->score <= 100 && rep_len > 0) {
+				pen_s1 = 0.01f * (r->score - rep_len);
+				pen_s1 = pen_s1 > 0.1f? pen_s1 : 0.1f;
+			}
+			pen_cm = pen_s1 < pen_cm? pen_s1 : pen_cm;
 			subsc = r->subsc > min_chain_sc? r->subsc : min_chain_sc;
 			if (r->p && r->p->dp_max2 > 0 && r->p->dp_max > 0) {
 				float identity = (float)(r->p->blen - r->p->n_diff - r->p->n_ambi) / (r->p->blen - r->p->n_ambi);
-				mapq = (int)(identity * pen_cm * q_coef * (1. - (float)r->p->dp_max2 * subsc / r->p->dp_max / r->score) * logf(r->score));
+				int mapq_alt = (int)(6.02f * identity * identity * (r->p->dp_max - r->p->dp_max2) / match_sc + .499f); // BWA-MEM like mapQ, mostly for short reads
+				mapq = (int)(identity * pen_cm * q_coef * (1. - (float)r->p->dp_max2 * subsc / r->p->dp_max / r->score) * logf(r->score)); // more for long reads
+				mapq = mapq < mapq_alt? mapq : mapq_alt; // in case the long-read heuristic fails
 			} else mapq = (int)(pen_cm * q_coef * (1. - (float)subsc / r->score) * logf(r->score));
+			mapq -= (int)(4.343f * logf(r->n_sub + 1) + .499f);
 			mapq = mapq > 0? mapq : 0;
 			r->mapq = mapq < 60? mapq : 60;
 		} else r->mapq = 0;

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"
@@ -17,6 +21,22 @@ typedef khash_t(idx) idxhash_t;
 
 #define kroundup64(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, (x)|=(x)>>32, ++(x))
 
+typedef struct mm_idx_bucket_s {
+	mm128_v a;   // (minimizer, position) array
+	int32_t n;   // size of the _p_ array
+	uint64_t *p; // position array for minimizers appearing >1 times
+	void *h;     // hash table indexing _p_ and minimizers appearing once
+} mm_idx_bucket_t;
+
+void mm_idxopt_init(mm_idxopt_t *opt)
+{
+	memset(opt, 0, sizeof(mm_idxopt_t));
+	opt->k = 15, opt->w = 10, opt->is_hpc = 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 *mi;
@@ -100,13 +120,13 @@ int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, ui
 	return en - st;
 }
 
-uint32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f)
+int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f)
 {
 	int i;
 	size_t n = 0;
 	uint32_t thres;
 	khint_t *a, k;
-	if (f <= 0.) return UINT32_MAX;
+	if (f <= 0.) return INT32_MAX;
 	for (i = 0; i < 1<<mi->b; ++i)
 		if (mi->B[i].h) n += kh_size((idxhash_t*)mi->B[i].h);
 	a = (uint32_t*)malloc(n * 4);
@@ -310,7 +330,7 @@ mm_idx_t *mm_idx_build(const char *fn, int w, int k, int is_hpc, int n_threads)
 	mm_idx_t *mi;
 	fp = mm_bseq_open(fn);
 	if (fp == 0) return 0;
-	mi = mm_idx_gen(fp, w, k, MM_IDX_DEF_B, is_hpc, 1<<18, n_threads, UINT64_MAX, 1);
+	mi = mm_idx_gen(fp, w, k, 14, is_hpc, 1<<18, n_threads, UINT64_MAX, 1);
 	mm_bseq_close(fp);
 	return mi;
 }
@@ -425,3 +445,43 @@ int mm_idx_is_idx(const char *fn)
 	close(fd);
 	return is_idx;
 }
+
+mm_idx_reader_t *mm_idx_reader_open(const char *fn, const mm_idxopt_t *opt, const char *fn_out)
+{
+	int is_idx;
+	mm_idx_reader_t *r;
+	is_idx = mm_idx_is_idx(fn);
+	if (is_idx < 0) return 0; // failed to open the index
+	r = (mm_idx_reader_t*)calloc(1, sizeof(mm_idx_reader_t));
+	r->is_idx = is_idx;
+	if (opt) r->opt = *opt;
+	else mm_idxopt_init(&r->opt);
+	if (r->is_idx) r->fp.idx = fopen(fn, "rb");
+	else r->fp.seq = mm_bseq_open(fn);
+	if (fn_out) r->fp_out = fopen(fn_out, "wb");
+	return r;
+}
+
+void mm_idx_reader_close(mm_idx_reader_t *r)
+{
+	if (r->is_idx) fclose(r->fp.idx);
+	else mm_bseq_close(r->fp.seq);
+	if (r->fp_out) fclose(r->fp_out);
+	free(r);
+}
+
+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))
+			fprintf(stderr, "[WARNING] Indexing parameters (-k, -w or -H) overridden by parameters used in the prebuilt index.\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);
+	if (mi) {
+		if (r->fp_out) mm_idx_dump(r->fp_out, mi);
+		++r->n_parts;
+	}
+	return mi;
+}

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));
+		t.w = (ktf_worker_t*)calloc(n_threads, sizeof(ktf_worker_t));
-		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)
 			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,39 +1,42 @@
-#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.2-r409"
 
+#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 },
 	{ "mb-size",        required_argument, 0, 'K' },
-	{ "int-rname",      no_argument,       0, 0 },
+	{ "int-rname",      no_argument,       0, 0 }, // obsolete; kept as a placeholder
 	{ "no-kalloc",      no_argument,       0, 0 },
 	{ "print-qname",    no_argument,       0, 0 },
 	{ "no-self",        no_argument,       0, 0 },
-	{ "print-seed",     no_argument,       0, 0 },
+	{ "print-seeds",    no_argument,       0, 0 },
 	{ "max-chain-skip", required_argument, 0, 0 },
 	{ "min-dp-len",     required_argument, 0, 0 },
 	{ "print-aln-seq",  no_argument,       0, 0 },
 	{ "splice",         no_argument,       0, 0 },
 	{ "cost-non-gt-ag", required_argument, 0, 0 },
 	{ "no-sam-sq",      no_argument,       0, 0 },
+	{ "approx-ext",     no_argument,       0, 0 },
 	{ "help",           no_argument,       0, 'h' },
 	{ "max-intron-len", required_argument, 0, 'G' },
 	{ "version",        no_argument,       0, 'V' },
@@ -59,24 +62,24 @@ static inline int64_t mm_parse_num(const char *str)
 int main(int argc, char *argv[])
 {
 	mm_mapopt_t opt;
-	int i, c, k = 15, w = -1, bucket_bits = MM_IDX_DEF_B, n_threads = 3, keep_name = 1, is_idx, is_hpc = 0, long_idx, idx_par_set = 0, max_intron_len = 0, n_idx_part = 0;
+	mm_idxopt_t ipt;
-	int minibatch_size = 200000000;
+	int i, c, n_threads = 3, long_idx, max_intron_len = 0;
-	uint64_t batch_size = 4000000000ULL;
-	mm_bseq_file_t *fp = 0;
 	char *fnw = 0, *rg = 0, *s;
-	FILE *fpr = 0, *fpw = 0, *fp_help = stderr;
+	FILE *fp_help = stderr;
+	mm_idx_reader_t *idx_rdr;
+	mm_idx_t *mi;
 
+	mm_verbose = 3;
 	liftrlimit();
 	mm_realtime0 = realtime();
-	mm_mapopt_init(&opt);
+	mm_set_opt(0, &ipt, &opt);
 
 	while ((c = getopt_long(argc, argv, "aSw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:h", long_options, &long_idx)) >= 0) {
-		if (c == 'w') w = atoi(optarg), idx_par_set = 1;
+		if (c == 'w') ipt.w = atoi(optarg);
-		else if (c == 'k') k = atoi(optarg), idx_par_set = 1;
+		else if (c == 'k') ipt.k = atoi(optarg);
-		else if (c == 'H') is_hpc = 1, idx_par_set = 1;
+		else if (c == 'H') ipt.is_hpc = 1;
 		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 == 'f') opt.mid_occ_frac = atof(optarg);
 		else if (c == 't') n_threads = atoi(optarg);
 		else if (c == 'v') mm_verbose = atoi(optarg);
 		else if (c == 'g') opt.max_gap = (int)mm_parse_num(optarg);
@@ -96,12 +99,11 @@ int main(int argc, char *argv[])
 		else if (c == 'B') opt.b = atoi(optarg);
 		else if (c == 'z') opt.zdrop = atoi(optarg);
 		else if (c == 's') opt.min_dp_max = atoi(optarg);
-		else if (c == 'I') batch_size = mm_parse_num(optarg);
+		else if (c == 'I') ipt.batch_size = mm_parse_num(optarg);
-		else if (c == 'K') minibatch_size = (int)mm_parse_num(optarg);
+		else if (c == 'K') ipt.mini_batch_size = (int)mm_parse_num(optarg);
 		else if (c == 'R') rg = optarg;
 		else if (c == 'h') fp_help = stdout;
-		else if (c == 0 && long_idx == 0) bucket_bits = atoi(optarg); // --bucket-bits
+		else if (c == 0 && long_idx == 0) ipt.bucket_bits = atoi(optarg); // --bucket-bits
-		else if (c == 0 && long_idx == 2) keep_name = 0; // --int-rname
 		else if (c == 0 && long_idx == 3) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
 		else if (c == 0 && long_idx == 4) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
 		else if (c == 0 && long_idx == 5) opt.flag |= MM_F_NO_SELF; // --no-self
@@ -112,9 +114,15 @@ int main(int argc, char *argv[])
 		else if (c == 0 && long_idx ==10) opt.flag |= MM_F_SPLICE; // --splice
 		else if (c == 0 && long_idx ==11) opt.noncan = atoi(optarg); // --cost-non-gt-ag
 		else if (c == 0 && long_idx ==12) opt.flag |= MM_F_NO_SAM_SQ; // --no-sam-sq
+		else if (c == 0 && long_idx ==13) opt.flag |= MM_F_APPROX_EXT; // --approx-ext
 		else if (c == 'V') {
 			puts(MM_VERSION);
 			return 0;
+		} else if (c == 'f') {
+			double x;
+			x = atof(optarg);
+			if (x < 1.0) opt.mid_occ_frac = x, opt.mid_occ = 0;
+			else opt.mid_occ = (int)(x + .499);
 		} else if (c == 'u') {
 			if (*optarg == 'b') opt.flag |= MM_F_SPLICE_FOR|MM_F_SPLICE_REV;
 			else if (*optarg == 'B') opt.flag |= MM_F_SPLICE_BOTH;
@@ -132,42 +140,12 @@ int main(int argc, char *argv[])
 			opt.e = opt.e2 = strtol(optarg, &s, 10);
 			if (*s == ',') opt.e2 = strtol(s + 1, &s, 10);
 		} else if (c == 'x') {
-			if (strcmp(optarg, "ava-ont") == 0) {
+			if (mm_set_opt(optarg, &ipt, &opt) < 0) {
-				opt.flag |= MM_F_AVA | MM_F_NO_SELF;
-				opt.min_chain_score = 100, opt.pri_ratio = 0.0f, opt.max_gap = 10000, opt.max_chain_skip = 25;
-				minibatch_size = 500000000;
-				k = 15, w = 5;
-			} else if (strcmp(optarg, "ava-pb") == 0) {
-				opt.flag |= MM_F_AVA | MM_F_NO_SELF;
-				opt.min_chain_score = 100, opt.pri_ratio = 0.0f, opt.max_gap = 10000, opt.max_chain_skip = 25;
-				minibatch_size = 500000000;
-				is_hpc = 1, k = 19, w = 5;
-			} else if (strcmp(optarg, "map10k") == 0 || strcmp(optarg, "map-pb") == 0) {
-				is_hpc = 1, k = 19;
-			} else if (strcmp(optarg, "map-ont") == 0) {
-				is_hpc = 0, k = 15;
-			} else if (strcmp(optarg, "asm5") == 0) {
-				k = 19, w = 19;
-				opt.a = 1, opt.b = 19, opt.q = 39, opt.q2 = 81, opt.e = 3, opt.e2 = 1, opt.zdrop = 200;
-				opt.min_dp_max = 200;
-			} else if (strcmp(optarg, "asm10") == 0) {
-				k = 19, w = 19;
-				opt.a = 1, opt.b = 9, opt.q = 16, opt.q2 = 41, opt.e = 2, opt.e2 = 1, opt.zdrop = 200;
-				opt.min_dp_max = 200;
-			} else if (strcmp(optarg, "splice") == 0 || strcmp(optarg, "cdna") == 0) {
-				k = 15, w = 5;
-				opt.flag |= MM_F_SPLICE | MM_F_SPLICE_FOR | MM_F_SPLICE_REV;
-				opt.max_gap = 2000, opt.max_gap_ref = opt.bw = 200000;
-				opt.a = 1, opt.b = 2, opt.q = 2, opt.e = 1, opt.q2 = 32, opt.e2 = 0;
-				opt.noncan = 5;
-				opt.zdrop = 200;
-			} else {
 				fprintf(stderr, "[E::%s] unknown preset '%s'\n", __func__, optarg);
 				return 1;
 			}
 		}
 	}
-	if (w < 0) w = (int)(.6666667 * k + .499);
 	if ((opt.flag & MM_F_SPLICE) && max_intron_len > 0)
 		opt.max_gap_ref = opt.bw = max_intron_len;
 
@@ -176,8 +154,8 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "Options:\n");
 		fprintf(fp_help, "  Indexing:\n");
 		fprintf(fp_help, "    -H           use homopolymer-compressed k-mer\n");
-		fprintf(fp_help, "    -k INT       k-mer size (no larger than 28) [%d]\n", k);
+		fprintf(fp_help, "    -k INT       k-mer size (no larger than 28) [%d]\n", ipt.k);
-		fprintf(fp_help, "    -w INT       minizer window size [{-k}*2/3]\n");
+		fprintf(fp_help, "    -w INT       minizer window size [%d]\n", ipt.w);
 		fprintf(fp_help, "    -I NUM       split index for every ~NUM input bases [4G]\n");
 		fprintf(fp_help, "    -d FILE      dump index to FILE []\n");
 		fprintf(fp_help, "  Mapping:\n");
@@ -206,7 +184,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -c           output CIGAR in PAF\n");
 		fprintf(fp_help, "    -S           output the cs tag in PAF (cs encodes both query and ref sequences)\n");
 		fprintf(fp_help, "    -t INT       number of threads [%d]\n", n_threads);
-		fprintf(fp_help, "    -K NUM       minibatch size [200M]\n");
+		fprintf(fp_help, "    -K NUM       minibatch size for mapping [200M]\n");
 //		fprintf(fp_help, "    -v INT       verbose level [%d]\n", mm_verbose);
 		fprintf(fp_help, "    --version    show version number\n");
 		fprintf(fp_help, "  Preset:\n");
@@ -218,54 +196,35 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "                 ava-pb: -Hk19 -w5 -Xp0 -m100 -g10000 -K500m --max-chain-skip 25 (PacBio read overlap)\n");
 		fprintf(fp_help, "                 ava-ont: -k15 -w5 -Xp0 -m100 -g10000 -K500m --max-chain-skip 25 (ONT read overlap)\n");
 		fprintf(fp_help, "                 splice: long-read spliced alignment (see minimap2.1 for details)\n");
+		fprintf(fp_help, "                 sr: short single-end reads without splicing (see minimap2.1 for details)\n");
 		fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of command-line options.\n");
 		return fp_help == stdout? 0 : 1;
 	}
 
-	is_idx = mm_idx_is_idx(argv[optind]);
+	idx_rdr = mm_idx_reader_open(argv[optind], &ipt, fnw);
-	if (is_idx < 0) {
+	if (idx_rdr == 0) {
 		fprintf(stderr, "[ERROR] failed to open file '%s'\n", argv[optind]);
 		return 1;
 	}
-	if (!is_idx && fnw == 0 && argc - optind < 2) {
+	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");
 		return 1;
 	}
-	if (is_idx) fpr = fopen(argv[optind], "rb");
-	else fp = mm_bseq_open(argv[optind]);
-	if (fnw) fpw = fopen(fnw, "wb");
 	if (opt.flag & MM_F_OUT_SAM)
 		mm_write_sam_hdr_no_SQ(rg, MM_VERSION, argc, argv);
-	for (;;) {
+	while ((mi = mm_idx_reader_read(idx_rdr, n_threads)) != 0) {
-		mm_idx_t *mi;
+		if (mm_verbose >= 2 && idx_rdr->n_parts > 1 && (opt.flag&MM_F_OUT_SAM) && !(opt.flag&MM_F_NO_SAM_SQ))
-		if (fpr) {
-			mi = mm_idx_load(fpr);
-			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 {
-			mi = mm_idx_gen(fp, w, k, bucket_bits, is_hpc, minibatch_size, n_threads, batch_size, keep_name);
-		}
-		if (mi == 0) break;
-		++n_idx_part;
-		if (mm_verbose >= 2 && n_idx_part > 1 && (opt.flag&MM_F_OUT_SAM) && !(opt.flag&MM_F_NO_SAM_SQ))
 			fprintf(stderr, "[WARNING] \033[1;31mSAM output is malformated due to internal @SQ lines. Please add option --no-sam-sq or filter afterwards.\033[0m\n");
 		if (mm_verbose >= 3)
 			fprintf(stderr, "[M::%s::%.3f*%.2f] loaded/built the index for %d target sequence(s)\n",
 					__func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), mi->n_seq);
-		if (fpw) {
-			mm_idx_dump(fpw, mi);
-			if (mm_verbose >= 3)
-				fprintf(stderr, "[M::%s::%.3f*%.2f] dumpped the (partial) index to disk\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0));
-		}
 		if (argc != optind + 1) mm_mapopt_update(&opt, mi);
 		if (mm_verbose >= 3) mm_idx_stat(mi);
 		for (i = optind + 1; i < argc; ++i)
-			mm_map_file(mi, argv[i], &opt, n_threads, minibatch_size);
+			mm_map_file(mi, argv[i], &opt, n_threads);
 		mm_idx_destroy(mi);
 	}
-	if (fpw) fclose(fpw);
+	mm_idx_reader_close(idx_rdr);
-	if (fpr) fclose(fpr);
-	if (fp)  mm_bseq_close(fp);
 
 	fprintf(stderr, "[M::%s] Version: %s\n", __func__, MM_VERSION);
 	fprintf(stderr, "[M::%s] CMD:", __func__);

map.c

@@ -11,7 +11,6 @@
 void mm_mapopt_init(mm_mapopt_t *opt)
 {
 	memset(opt, 0, sizeof(mm_mapopt_t));
-	opt->max_occ_frac = 1e-5f;
 	opt->mid_occ_frac = 2e-4f;
 	opt->sdust_thres = 0;
 
@@ -34,21 +33,72 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->zdrop = 400;
 	opt->min_dp_max = opt->min_chain_score * opt->a;
 	opt->min_ksw_len = 200;
+	opt->mini_batch_size = 200000000;
 }
 
 void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)
 {
 	if (opt->flag & MM_F_SPLICE_BOTH)
 		opt->flag &= ~(MM_F_SPLICE_FOR|MM_F_SPLICE_REV);
-	opt->max_occ = mm_idx_cal_max_occ(mi, opt->max_occ_frac);
+	if (opt->mid_occ <= 0)
-	opt->mid_occ = mm_idx_cal_max_occ(mi, opt->mid_occ_frac);
+		opt->mid_occ = mm_idx_cal_max_occ(mi, opt->mid_occ_frac);
 	if (mm_verbose >= 3)
-		fprintf(stderr, "[M::%s::%.3f*%.2f] mid_occ = %d; max_occ = %d\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0),
+		fprintf(stderr, "[M::%s::%.3f*%.2f] mid_occ = %d\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), opt->mid_occ);
-				opt->mid_occ, opt->max_occ);
+}
+
+int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
+{
+	if (preset == 0) {
+		mm_idxopt_init(io);
+		mm_mapopt_init(mo);
+	} else if (strcmp(preset, "ava-ont") == 0) {
+		io->is_hpc = 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;
+		mo->mini_batch_size = 500000000;
+	} else if (strcmp(preset, "ava-pb") == 0) {
+		io->is_hpc = 1, 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;
+		mo->mini_batch_size = 500000000;
+	} else if (strcmp(preset, "map10k") == 0 || strcmp(preset, "map-pb") == 0) {
+		io->is_hpc = 1, io->k = 19;
+	} else if (strcmp(preset, "map-ont") == 0) {
+		io->is_hpc = 0, io->k = 15;
+	} else if (strcmp(preset, "asm5") == 0) {
+		io->is_hpc = 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;
+	} else if (strcmp(preset, "asm10") == 0) {
+		io->is_hpc = 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;
+	} else if (strcmp(preset, "short") == 0 || strcmp(preset, "sr") == 0) {
+		io->is_hpc = 0, io->k = 21, io->w = 11;
+		mo->flag |= MM_F_APPROX_EXT;
+		mo->a = 2, mo->b = 8, mo->q = 12, mo->e = 2, mo->q2 = 32, mo->e2 = 1;
+		mo->max_gap = 100;
+		mo->pri_ratio = 0.5f;
+		mo->min_cnt = 2;
+		mo->min_chain_score = 20;
+		mo->min_dp_max = 40;
+		mo->best_n = 20;
+		mo->bw = 50;
+		mo->mid_occ = 1000;
+		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;
+		mo->flag |= MM_F_SPLICE | MM_F_SPLICE_FOR | MM_F_SPLICE_REV;
+		mo->max_gap = 2000, mo->max_gap_ref = mo->bw = 200000;
+		mo->a = 1, mo->b = 2, mo->q = 2, mo->e = 1, mo->q2 = 32, mo->e2 = 0;
+		mo->noncan = 5;
+		mo->zdrop = 200;
+	} else return -1;
+	return 0;
 }
 
 typedef struct {
-	uint32_t n:31, is_alloc:1;
+	uint32_t n;
 	uint32_t qpos;
 	union {
 		const uint64_t *cr;
@@ -102,149 +152,85 @@ static void mm_dust_minier(mm128_v *mini, int l_seq, const char *seq, int sdust_
 	}
 	mini->n = k;
 }
-#if 0
-int mm_pair_thin_core(mm_tbuf_t *b, uint64_t x, int radius, int rel, int st0, int n, const uint64_t *z, uint64_v *a)
-{
-	int i, st = st0, en = n, mid = en - 1;
-	while (st < en) {
-		uint64_t y;
-		mid = st + ((en - st) >> 1);
-		y = z[mid];
-		if (y < x && (x - y)>>1 > radius) st = mid + 1;
-		else if (y >= x && (y - x)>>1 > radius) en = mid;
-		else break;
-	}
-	if (st < en) {
-		for (en = mid + 1; en < n; ++en)
-			if (z[en] > x && (z[en] - x)>>1 > radius)
-				break;
-		for (st = mid - 1; st >= st0; --st)
-			if (z[st] < x && (x - z[st])>>1 > radius)
-				break;
-		++st;
-		for (i = st; i < en; ++i) {
-			uint64_t y = z[i];
-			if (((x ^ y) & 1) == rel) {
-//				printf("* %d,%d\n", (uint32_t)x>>1, (uint32_t)y>>1);
-				kv_push(uint64_t, b->km, *a, y);
-			}
-		}
-		return en;
-	} else return st < n && z[st] < x? st + 1 : en;
-}
 
-void mm_pair_thin(mm_tbuf_t *b, int radius, mm_match_t *m1, mm_match_t *m2)
+mm_reg1_t *mm_map(const mm_idx_t *mi, int qlen, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
 {
-	mm_match_t *m[2];
+	int i, n, j, n_u, max_gap_ref, rep_st = 0, rep_en = 0, rep_len = 0;
-	const uint64_t *z[2];
-	uint64_v a[2];
-	int i, n[2], k[2], u = 0, rel = (m1->qpos ^ m2->qpos) & 1;
-
-	m[0] = m1, m[1] = m2;
-	for (i = 0; i < 2; ++i) {
-		n[i] = m[i]->n;
-		z[i] = m[i]->x.cr;
-		k[i] = 0;
-		kv_init(a[i]);
-		kv_resize(uint64_t, b->km, a[i], 256);
-	}
-	while (k[0] < n[0] && k[1] < n[1]) {
-		//printf("%d; %d,%d\n", u, k[0], k[1]);
-		int v = u^1, dist = (int)(m[v]->qpos>>1) - (int)(m[u]->qpos>>1);
-		uint64_t x = z[u][k[u]];
-		int uori = (x ^ m[u]->qpos) & 1, last;
-		int64_t tpos = x>>1 & 0x7fffffff;
-		tpos = uori == 0? tpos + dist : tpos - dist;
-		if (tpos < 0) tpos = 0;
-		x = x>>32<<32 | tpos<<1 | (x&1);
-		last = a[v].n;
-		k[v] = mm_pair_thin_core(b, x, radius, rel, k[v], n[v], z[v], &a[v]);
-		if (a[v].n > last) kv_push(uint64_t, b->km, a[u], z[u][k[u]]);
-		++k[u];
-		u ^= 1;
-	}
-	for (i = 0; i < 2; ++i)
-		m[i]->n = a[i].n, m[i]->x.r = a[i].a, m[i]->is_alloc = 1;
-//	printf("%d,%d; %d,%d\n", m[0]->qpos>>1, m[1]->qpos>>1, m[0]->n, m[1]->n);
-}
-#endif
-mm_reg1_t *mm_map_frag(const mm_mapopt_t *opt, const mm_idx_t *mi, mm_tbuf_t *b, uint32_t m_st, uint32_t m_en, const char *qname, int qlen, const char *seq, int *n_regs)
-{
-	int i, n = m_en - m_st, j, n_u, max_gap_ref;
 	int64_t n_a;
 	uint64_t *u;
 	mm_match_t *m;
 	mm128_t *a;
 	mm_reg1_t *regs;
 
+	// collect minimizers
+	b->mini.n = 0;
+	mm_sketch(b->km, seq, qlen, mi->w, mi->k, 0, mi->is_hpc, &b->mini);
+	n = b->mini.n;
+
+	if (opt->sdust_thres > 0)
+		mm_dust_minier(&b->mini, qlen, seq, opt->sdust_thres, b->sdb);
+
 	// convert to local representation
 	m = (mm_match_t*)kmalloc(b->km, n * sizeof(mm_match_t));
 	for (i = 0; i < n; ++i) {
 		int t;
-		mm128_t *p = &b->mini.a[i + m_st];
+		mm128_t *p = &b->mini.a[i];
-		m[i].is_alloc = 0;
 		m[i].qpos = (uint32_t)p->y;
 		m[i].x.cr = mm_idx_get(mi, p->x>>8, &t);
 		m[i].n = t;
 	}
-#if 0
+
-	int last = -1, last2 = -1;
-	// pair k-mer thinning
-	for (i = 0; i < n; ++i) {
-		if (m[i].n >= opt->mid_occ && m[i].n < opt->max_occ) {
-			if (last2 < 0) last2 = i;
-			if (last < 0 || m[last].n < m[i].n) last = i;
-			if (last >= 0 && (m[last].qpos>>1) + (m[last].span>>1) <= m[i].qpos>>1) {
-				mm_pair_thin(b, opt->bw, &m[last], &m[i]);
-				last2 = last = -1;
-			} else if (last2 >= 0 && (m[last2].qpos>>1) + (m[last2].span>>1) <= m[i].qpos>>1) {
-				mm_pair_thin(b, opt->bw, &m[last2], &m[i]);
-				last2 = last = -1;
-			}
-		}
-	}
-#endif
 	// fill the _a_ array
 	for (i = 0, n_a = 0; i < n; ++i) // find the length of a[]
 		if (m[i].n < opt->mid_occ) n_a += m[i].n;
 	a = (mm128_t*)kmalloc(b->km, n_a * sizeof(mm128_t));
 	for (i = j = 0; i < n; ++i) {
-		mm128_t *p = &b->mini.a[i + m_st];
+		mm128_t *p = &b->mini.a[i];
 		mm_match_t *q = &m[i];
 		const uint64_t *r = q->x.cr;
 		int k, q_span = p->x & 0xff, is_tandem = 0;
-		if (q->n >= opt->mid_occ) continue;
+		if (q->n >= opt->mid_occ) {
-		if (i > 0 && p->x>>8 == b->mini.a[m_st + i - 1].x>>8) is_tandem = 1;
+			int en = (q->qpos>>1) + 1, st = en - q_span;
-		if (i < n - 1 && p->x>>8 == b->mini.a[m_st + i + 1].x>>8) is_tandem = 1;
+			if (st > rep_en) {
+				rep_len += rep_en - rep_st;
+				rep_st = st, rep_en = en;
+			} else rep_en = en;
+			continue;
+		}
+		if (i > 0 && p->x>>8 == b->mini.a[i - 1].x>>8) is_tandem = 1;
+		if (i < n - 1 && p->x>>8 == b->mini.a[i + 1].x>>8) is_tandem = 1;
 		for (k = 0; k < q->n; ++k) {
-			const char *tname = mi->seq[r[k]>>32].name;
 			int32_t rpos = (uint32_t)r[k] >> 1;
 			mm128_t *p;
-			if (qname && (opt->flag&MM_F_NO_SELF) && strcmp(qname, tname) == 0 && rpos == (q->qpos>>1)) // avoid the diagonal
+			if (qname && (opt->flag&(MM_F_NO_SELF|MM_F_AVA))) {
-				continue;
+				const char *tname = mi->seq[r[k]>>32].name;
-			if (qname && (opt->flag&MM_F_AVA) && strcmp(qname, tname) > 0) // all-vs-all mode: map once
+				if ((opt->flag&MM_F_NO_SELF) && strcmp(qname, tname) == 0 && rpos == (q->qpos>>1)) // avoid the diagonal
-				continue;
+					continue;
+				if ((opt->flag&MM_F_AVA) && strcmp(qname, tname) > 0) // all-vs-all mode: map once
+					continue;
+			}
 			p = &a[j++];
 			if ((r[k]&1) == (q->qpos&1)) { // forward strand
-				p->x = (r[k]&0xffffffff00000000ULL) | (uint32_t)r[k]>>1;
+				p->x = (r[k]&0xffffffff00000000ULL) | rpos;
 				p->y = (uint64_t)q_span << 32 | q->qpos >> 1;
 			} else { // reverse strand
-				p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | (uint32_t)r[k]>>1;
+				p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | rpos;
 				p->y = (uint64_t)q_span << 32 | (qlen - ((q->qpos>>1) + 1 - q_span) - 1);
 			}
 			if (is_tandem) p->y |= MM_SEED_TANDEM;
 		}
 	}
+	rep_len += rep_en - rep_st;
 	n_a = j;
 	radix_sort_128x(a, a + n_a);
-	for (i = 0; i < n; ++i)
-		if (m[i].is_alloc) kfree(b->km, m[i].x.r);
 	kfree(b->km, m);
 
-	if (mm_dbg_flag & MM_DBG_PRINT_SEED)
+	if (mm_dbg_flag & MM_DBG_PRINT_SEED) {
+		fprintf(stderr, "RS\t%d\n", rep_len);
 		for (i = 0; i < n_a; ++i)
 			fprintf(stderr, "SD\t%s\t%d\t%c\t%d\t%d\t%d\n", mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff),
 					i == 0? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));
+	}
 
 	max_gap_ref = opt->max_gap_ref >= 0? opt->max_gap_ref : opt->max_gap;
 	n_u = mm_chain_dp(max_gap_ref, opt->max_gap, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, !!(opt->flag&MM_F_SPLICE), n_a, a, &u, b->km);
@@ -258,7 +244,7 @@ mm_reg1_t *mm_map_frag(const mm_mapopt_t *opt, const mm_idx_t *mi, mm_tbuf_t *b,
 						i == regs[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));
 
 	if (!(opt->flag & MM_F_AVA)) { // don't choose primary mapping(s) for read overlap
-		mm_set_parent(b->km, opt->mask_level, *n_regs, regs);
+		mm_set_parent(b->km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b);
 		mm_select_sub(b->km, opt->mask_level, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 		if (!(opt->flag & MM_F_SPLICE))
 			mm_join_long(b->km, opt, qlen, n_regs, regs, a); // TODO: this can be applied to all-vs-all in principle
@@ -266,12 +252,12 @@ mm_reg1_t *mm_map_frag(const mm_mapopt_t *opt, const mm_idx_t *mi, mm_tbuf_t *b,
 	if (opt->flag & MM_F_CIGAR) {
 		regs = mm_align_skeleton(b->km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs()
 		if (!(opt->flag & MM_F_AVA)) {
-			mm_set_parent(b->km, opt->mask_level, *n_regs, regs);
+			mm_set_parent(b->km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b);
 			mm_select_sub(b->km, opt->mask_level, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 			mm_set_sam_pri(*n_regs, regs);
 		}
 	}
-	mm_set_mapq(*n_regs, regs, opt->min_chain_score);
+	mm_set_mapq(*n_regs, regs, opt->min_chain_score, opt->a, rep_len);
 
 	// free
 	kfree(b->km, a);
@@ -279,17 +265,6 @@ mm_reg1_t *mm_map_frag(const mm_mapopt_t *opt, const mm_idx_t *mi, mm_tbuf_t *b,
 	return regs;
 }
 
-mm_reg1_t *mm_map(const mm_idx_t *mi, int l_seq, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
-{
-	mm_reg1_t *regs;
-	b->mini.n = 0;
-	mm_sketch(b->km, seq, l_seq, mi->w, mi->k, 0, mi->is_hpc, &b->mini);
-	if (opt->sdust_thres > 0)
-		mm_dust_minier(&b->mini, l_seq, seq, opt->sdust_thres, b->sdb);
-	regs = mm_map_frag(opt, mi, b, 0, b->mini.n, qname, l_seq, seq, n_regs);
-	return regs;
-}
-
 /**************************
  * Multi-threaded mapping *
  **************************/
@@ -377,14 +352,18 @@ static void *worker_pipeline(void *shared, int step, void *in)
     return 0;
 }
 
-int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int n_threads, int mini_batch_size)
+int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int n_threads)
 {
 	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;
+	pl.n_threads = n_threads, pl.mini_batch_size = opt->mini_batch_size;
 	if ((opt->flag & MM_F_OUT_SAM) && !(opt->flag & MM_F_NO_SAM_SQ))
 		mm_write_sam_SQ(idx);
 	kt_pipeline(n_threads == 1? 1 : 2, worker_pipeline, &pl, 3);

minimap.h

@@ -5,8 +5,6 @@
 #include <stdio.h>
 #include <sys/types.h>
 
-#define MM_IDX_DEF_B     14
-
 #define MM_F_NO_SELF     0x001
 #define MM_F_AVA         0x002
 #define MM_F_CIGAR       0x004
@@ -19,6 +17,7 @@
 #define MM_F_SPLICE_REV  0x200
 #define MM_F_SPLICE_BOTH 0x400
 #define MM_F_NO_SAM_SQ   0x800
+#define MM_F_APPROX_EXT  0x1000
 
 #define MM_IDX_MAGIC   "MMI\2"
 
@@ -26,21 +25,11 @@
 extern "C" {
 #endif
 
-typedef struct {
+// emulate 128-bit integers and arrays
- 	uint64_t x, y;
+typedef struct { uint64_t x, y; } mm128_t;
-} mm128_t;
-
 typedef struct { size_t n, m; mm128_t *a; } mm128_v;
-typedef struct { size_t n, m; uint64_t *a; } uint64_v;
-typedef struct { size_t n, m; uint32_t *a; } uint32_v;
-
-typedef struct {
-	mm128_v a;   // (minimizer, position) array
-	int32_t n;   // size of the _p_ array
-	uint64_t *p; // position array for minimizers appearing >1 times
-	void *h;     // hash table indexing _p_ and minimizers appearing once
-} mm_idx_bucket_t;
 
+// minimap2 index
 typedef struct {
 	char *name;      // name of the db sequence
 	uint64_t offset; // offset in mm_idx_t::S
@@ -49,103 +38,216 @@ typedef struct {
 
 typedef struct {
 	int32_t b, w, k, is_hpc;
-	uint32_t n_seq;     // number of reference sequences
+	uint32_t n_seq;            // number of reference sequences
-	mm_idx_seq_t *seq;  // sequence name, length and offset
+	mm_idx_seq_t *seq;         // sequence name, length and offset
-	uint32_t *S;        // 4-bit packed sequence
+	uint32_t *S;               // 4-bit packed sequence
-	mm_idx_bucket_t *B; // index
+	struct mm_idx_bucket_s *B; // index (hidden)
 	void *km;
 } mm_idx_t;
 
+// minimap2 alignment
 typedef struct {
-	uint32_t capacity;
+	uint32_t capacity;                  // the capacity of cigar[]
-	int32_t dp_score, dp_max, dp_max2;
+	int32_t dp_score, dp_max, dp_max2;  // DP score; score of the max-scoring segment; score of the best alternate mappings
-	uint32_t blen;
+	uint32_t blen;                      // block length
-	uint32_t n_diff;
+	uint32_t n_diff;                    // number of differences, including ambiguous bases
-	uint32_t n_ambi:30, trans_strand:2;
+	uint32_t n_ambi:30, trans_strand:2; // number of ambiguous bases; transcript strand: 0 for unknown, 1 for +, 2 for -
-	uint32_t n_cigar;
+	uint32_t n_cigar;                   // number of cigar operations in cigar[]
 	uint32_t cigar[];
 } mm_extra_t;
 
 typedef struct {
-	int32_t id;
+	int32_t id;                     // ID for internal uses (see also parent below)
-	uint32_t cnt:31, rev:1;
+	uint32_t cnt:31, rev:1;         // number of minimizers; if on the reverse strand
-	uint32_t rid:31, inv:1;
+	uint32_t rid:31, inv:1;         // reference index; if this is an alignment from inversion rescue
-	int32_t score;
+	int32_t score;                  // DP alignment score
-	int32_t qs, qe, rs, re;
+	int32_t qs, qe, rs, re;         // query start and end; reference start and end
-	int32_t parent, subsc;
+	int32_t parent, subsc;          // parent==id if primary; best alternate mapping score
-	int32_t as;
+	int32_t as;                     // offset in the a[] array (for internal uses only)
-	int32_t fuzzy_mlen, fuzzy_blen;
+	int32_t fuzzy_mlen, fuzzy_blen; // seeded exact match length; seeded alignment block length (approximate)
-	uint32_t mapq:8, split:2, sam_pri:1, n_sub:21; // TODO: n_sub is not used for now
+	uint32_t mapq:8, split:2, sam_pri:1, n_sub:21; // mapQ; split pattern; if SAM primary; number of suboptimal mappings
 	mm_extra_t *p;
 } mm_reg1_t;
 
+// indexing and mapping options
 typedef struct {
-	float max_occ_frac;
+	short k, w, is_hpc, bucket_bits;
-	float mid_occ_frac;
+	int mini_batch_size;
-	int sdust_thres;  // score threshold for SDUST; 0 to disable
+	uint64_t batch_size;
-	int flag;    // see MM_F_* macros
+} mm_idxopt_t;
 
-	int bw;  // bandwidth
+typedef struct {
+	int sdust_thres; // score threshold for SDUST; 0 to disable
+	int flag;        // see MM_F_* macros
+
+	int bw;          // bandwidth
 	int max_gap, max_gap_ref; // break a chain if there are no minimizers in a max_gap window
 	int max_chain_skip;
-	int min_cnt;
+	int min_cnt;         // min number of minimizers on each chain
-	int min_chain_score;
+	int min_chain_score; // min chaining score
 
 	float mask_level;
 	float pri_ratio;
-	int best_n;
+	int best_n;      // top best_n chains are subjected to DP alignment
 
 	int max_join_long, max_join_short;
 	int min_join_flank_sc;
 
 	int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties
-	int noncan;
+	int noncan;      // cost of non-canonical splicing sites
-	int zdrop;
+	int zdrop;       // break alignment if alignment score drops too fast along the diagonal
-	int min_dp_max;
+	int min_dp_max;  // drop an alignment if the score of the max scoring segment is below this threshold
 	int min_ksw_len;
 
-	int max_occ;
+	float mid_occ_frac;  // only used by mm_mapopt_update(); see below
-	int mid_occ;
+	int32_t mid_occ;     // ignore seeds with occurrences above this threshold
+	int mini_batch_size; // size of a batch of query bases to process in parallel
 } mm_mapopt_t;
 
-extern int mm_verbose, mm_dbg_flag;
+// index reader
-extern double mm_realtime0;
+typedef struct {
+	int is_idx, n_parts;
+	mm_idxopt_t opt;
+	FILE *fp_out;
+	union {
+		struct mm_bseq_file_s *seq;
+		FILE *idx;
+	} fp;
+} mm_idx_reader_t;
 
-struct mm_tbuf_s;
+// memory buffer for thread-local storage during mapping
 typedef struct mm_tbuf_s mm_tbuf_t;
 
-struct mm_bseq_file_s;
+// global variables
+extern int mm_verbose, mm_dbg_flag; // verbose level: 0 for no info, 1 for error, 2 for warning, 3 for message (default); debugging flag
+extern double mm_realtime0; // wall-clock timer
 
-#define mm_seq4_set(s, i, c) ((s)[(i)>>3] |= (uint32_t)(c) << (((i)&7)<<2))
+/**
-#define mm_seq4_get(s, i)    ((s)[(i)>>3] >> (((i)&7)<<2) & 0xf)
+ * Set default or preset parameters
+ *
+ * @param preset     NULL to set all parameters as default; otherwise apply preset to affected parameters
+ * @param io         pointer to indexing parameters
+ * @param mo         pointer to mapping parameters
+ *
+ * @return 0 if success; -1 if _present_ unknown
+ */
+int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo);
 
-// compute minimizers
+/**
-void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, int is_hpc, mm128_v *p);
+ * Update mm_mapopt_t::mid_occ via mm_mapopt_t::mid_occ_frac
-
+ *
-// minimizer indexing
+ * If mm_mapopt_t::mid_occ is 0, this function sets it to a number such that no
-mm_idx_t *mm_idx_init(int w, int k, int b, int is_hpc);
+ * more than mm_mapopt_t::mid_occ_frac of minimizers in the index have a higher
-void mm_idx_destroy(mm_idx_t *mi);
+ * occurrence.
-mm_idx_t *mm_idx_gen(struct mm_bseq_file_s *fp, int w, int k, int b, int is_hpc, int mini_batch_size, int n_threads, uint64_t batch_size, int keep_name);
+ *
-uint32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f);
+ * @param opt        mapping parameters
-void mm_idx_stat(const mm_idx_t *idx);
+ * @param mi         minimap2 index
-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);
-
-mm_idx_t *mm_idx_build(const char *fn, int w, int k, int is_hpc, int n_threads);
-int mm_idx_is_idx(const char *fn);
-
-// minimizer index I/O
-void mm_idx_dump(FILE *fp, const mm_idx_t *mi);
-mm_idx_t *mm_idx_load(FILE *fp);
-
-// mapping
-void mm_mapopt_init(mm_mapopt_t *opt);
 void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi);
+
+/**
+ * Initialize an index reader
+ *
+ * @param fn         index or fasta/fastq file name (this function tests the file type)
+ * @param opt        indexing parameters
+ * @param fn_out     if not NULL, write built index to this file
+ *
+ * @return an index reader on success; NULL if fail to open _fn_
+ */
+mm_idx_reader_t *mm_idx_reader_open(const char *fn, const mm_idxopt_t *opt, const char *fn_out);
+
+/**
+ * Read/build an index
+ *
+ * If the input file is an index file, this function reads one part of the
+ * index and returns. If the input file is a sequence file (fasta or fastq),
+ * this function constructs the index for about mm_idxopt_t::batch_size bases.
+ * Importantly, for a huge collection of sequences, this function may only
+ * return an index for part of sequences. It needs to be repeatedly called
+ * to traverse the entire index/sequence file.
+ *
+ * @param r          index reader
+ * @param n_threads  number of threads for constructing index
+ *
+ * @return an index on success; NULL if reaching the end of the input file
+ */
+mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads);
+
+/**
+ * Destroy/deallocate an index reader
+ *
+ * @param r          index reader
+ */
+void mm_idx_reader_close(mm_idx_reader_t *r);
+
+/**
+ * Print index statistics to stderr
+ *
+ * @param mi         minimap2 index
+ */
+void mm_idx_stat(const mm_idx_t *idx);
+
+/**
+ * Destroy/deallocate an index
+ *
+ * @param r          minimap2 index
+ */
+void mm_idx_destroy(mm_idx_t *mi);
+
+/**
+ * Initialize a thread-local buffer for mapping
+ *
+ * Each mapping thread requires a buffer specific to the thread (see mm_map()
+ * below). The primary purpose of this buffer is to reduce frequent heap
+ * allocations across threads. A buffer shall not be used by two or more
+ * threads.
+ *
+ * @return pointer to a thread-local buffer
+ */
 mm_tbuf_t *mm_tbuf_init(void);
+
+/**
+ * Destroy/deallocate a thread-local buffer for mapping
+ *
+ * @param b          the buffer
+ */
 void mm_tbuf_destroy(mm_tbuf_t *b);
+
+/**
+ * Align a query sequence against an index
+ *
+ * This function possibly finds multiple alignments of the query sequence.
+ * The returned array and the mm_reg1_t::p field of each element are allocated
+ * with malloc().
+ *
+ * @param mi         minimap2 index
+ * @param l_seq      length of the query sequence
+ * @param seq        the query sequence
+ * @param n_regs     number of hits (out)
+ * @param b          thread-local buffer; two mm_map() calls shall not use one buffer at the same time!
+ * @param opt        mapping parameters
+ * @param name       query name, used for all-vs-all overlapping and debugging
+ *
+ * @return an array of hits which need to be deallocated with free() together
+ *         with mm_reg1_t::p of each element. The size is written to _n_regs_.
+ */
 mm_reg1_t *mm_map(const mm_idx_t *mi, int l_seq, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *name);
 
-int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int n_threads, int tbatch_size);
+/**
+ * Align a fasta/fastq file and print alignments to stdout
+ *
+ * @param idx        minimap2 index
+ * @param fn         fasta/fastq file name
+ * @param opt        mapping parameters
+ * @param n_threads  number of threads
+ *
+ * @return 0 on success; -1 if _fn_ can't be read
+ */
+int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int n_threads);
+
+// 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);
 
 #ifdef __cplusplus
 }

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "25 August 2017" "minimap2-2.1-r311" "Bioinformatics tools"
+.TH minimap2 1 "17 September 2017" "minimap2-2.2 (r409)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -333,6 +333,12 @@ CIGAR operator; 2) long insertions are disabled; 3) deletion and insertion gap
 costs are different during chaining; 4) the computation of the
 .RB ` ms '
 tag ignores introns to demote hits to pseudogenes.
+.TP
+.B sr
+Short single-end reads without splicing
+.RB ( -k21
+.B -w11 -A2 -B8 -O12,32 -E2,1 -r50 -p.5 -N20 -f1000 -n2 -m20 -s40 -g100 -K50m
+.BR --approx-ext ).
 .RE
 .SS Miscellaneous options
 .TP 10
@@ -392,6 +398,7 @@ NM	i	Total number of mismatches and gaps in the alignment
 AS	i	DP alignment score
 ms	i	DP score of the max scoring segment in the alignment
 nn	i	Number of ambiguous bases in the alignment
+ts	A	Transcript strand (splice mode only)
 cg	Z	CIGAR string (only in PAF)
 .TE
 

misc.c

@@ -1,19 +1,104 @@
-#include <sys/resource.h>
-#include <sys/time.h>
 #include "minimap.h"
 
-int mm_verbose = 3;
+int mm_verbose = 1;
 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(void)
 {
 	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()
+double realtime(void)
 {
 	struct timeval tp;
 	struct timezone tzp;

misc/sim-eval.js

@@ -181,11 +181,11 @@ var sum_tot = 0, sum_err = 0, q_out = -1, sum_tot2 = 0, sum_err2 = 0;
 for (var q = max_mapq; q >= 0; --q) {
 	if (tot[q] == 0) continue;
 	if (q_out < 0 || err[q] > 0) {
-		if (q_out >= 0) print('Q', q_out, sum_tot, sum_err, (sum_err2/sum_tot2).toFixed(9));
+		if (q_out >= 0) print('Q', q_out, sum_tot, sum_err, (sum_err2/sum_tot2).toFixed(9), sum_tot2);
 		sum_tot = sum_err = 0, q_out = q;
 	}
 	sum_tot += tot[q], sum_err += err[q];
 	sum_tot2 += tot[q], sum_err2 += err[q];
 }
-print('Q', q_out, sum_tot, sum_err, (sum_err2/sum_tot2).toFixed(9));
+print('Q', q_out, sum_tot, sum_err, (sum_err2/sum_tot2).toFixed(9), sum_tot2);
 if (n_unmapped != null) print('U', n_unmapped);

mmpriv.h

@@ -21,6 +21,9 @@
 #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
 #endif
 
+#define mm_seq4_set(s, i, c) ((s)[(i)>>3] |= (uint32_t)(c) << (((i)&7)<<2))
+#define mm_seq4_get(s, i)    ((s)[(i)>>3] >> (((i)&7)<<2) & 0xf)
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -40,10 +43,17 @@ void radix_sort_128x(mm128_t *beg, mm128_t *end);
 void radix_sort_64(uint64_t *beg, uint64_t *end);
 uint32_t ks_ksmall_uint32_t(size_t n, uint32_t arr[], size_t kk);
 
+void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, int is_hpc, mm128_v *p);
+
 void mm_write_sam_SQ(const mm_idx_t *idx);
 void mm_write_sam_hdr_no_SQ(const char *rg, const char *ver, int argc, char *argv[]);
 void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag);
 void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int n_regs, const mm_reg1_t *regs);
+
+void mm_idxopt_init(mm_idxopt_t *opt);
+const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);
+int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
+int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f);
 int 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, int64_t n, mm128_t *a, uint64_t **_u, void *km);
 mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
 
@@ -51,12 +61,12 @@ mm_reg1_t *mm_gen_regs(void *km, 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);
 void mm_sync_regs(void *km, int n_regs, mm_reg1_t *regs);
 int mm_set_sam_pri(int n, mm_reg1_t *r);
-void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r);
+void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff);
 void mm_select_sub(void *km, float mask_level, float pri_ratio, int min_diff, int best_n, int *n_, mm_reg1_t *r);
 void mm_filter_regs(void *km, const mm_mapopt_t *opt, int *n_regs, mm_reg1_t *regs);
 void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs, mm128_t *a);
 void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r);
-void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc);
+void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len);
 
 #ifdef __cplusplus
 }

python/README.rst

@@ -0,0 +1,141 @@
+==============================
+Mappy: Minimap2 Python Binding
+==============================
+
+Mappy provides a convenient interface to `minimap2
+<https://github.com/lh3/minimap2>`_, a fast and accurate C program to align
+genomic and transcribe nucleotide sequences.
+
+Installation
+------------
+
+Mappy depends on `zlib <http://zlib.net>`_. It can be installed with `pip
+<https://en.wikipedia.org/wiki/Pip_(package_manager)>`_:
+
+.. code:: shell
+
+	pip install --user mappy
+
+or from the minimap2 github repo (`Cython <http://cython.org>`_ required):
+
+.. code:: shell
+
+	git clone https://github.com/lh3/minimap2
+	cd minimap2
+	python setup.py install
+
+Usage
+-----
+
+The following Python script demonstrates the key functionality of mappy:
+
+.. code:: python
+
+	import mappy as mp
+	a = mp.Aligner("test/MT-human.fa")  # load or build index
+	if not a: raise Exception("ERROR: failed to load/build index")
+	for name, seq, qual in mp.fastx_read("test/MT-orang.fa"): # read a fasta/q sequence
+		for hit in a.map(seq): # traverse alignments
+			print("{}\t{}\t{}\t{}".format(hit.ctg, hit.r_st, hit.r_en, hit.cigar_str))
+
+APIs
+----
+
+Mappy implements two classes and one global function.
+
+Class mappy.Aligner
+~~~~~~~~~~~~~~~~~~~
+
+.. code:: python
+
+	mappy.Aligner(fn_idx_in, preset=None, ...)
+
+This constructor accepts the following arguments:
+
+* **fn_idx_in**: index or sequence file name. Minimap2 automatically tests the
+  file type. If a sequence file is provided, minimap2 builds an index. The
+  sequence file can be optionally gzip'd.
+
+* **preset**: minimap2 preset. Currently, minimap2 supports the following
+  presets: **sr** for single-end short reads; **map-pb** for PacBio
+  read-to-reference mapping; **map-ont** for Oxford Nanopore read mapping;
+  **splice** for long-read spliced alignment; **asm5** for assembly-to-assembly
+  alignment; **asm10** for full genome alignment of closely related species. Note
+  that the Python module does not support all-vs-all read overlapping.
+
+* **k**: k-mer length, no larger than 28
+
+* **w**: minimizer window size, no larger than 255
+
+* **min_cnt**: mininum number of minimizers on a chain
+
+* **min_chain_score**: minimum chaing score
+
+* **bw**: chaining and alignment band width
+
+* **best_n**: max number of alignments to return
+
+* **n_threads**: number of indexing threads; 3 by default
+
+* **fn_idx_out**: name of file to which the index is written
+
+.. code:: python
+
+	mappy.Aligner.map(seq)
+
+This method aligns :code:`seq` against the index. It is a generator, *yielding*
+a series of :code:`mappy.Alignment` objects.
+
+Class mappy.Alignment
+~~~~~~~~~~~~~~~~~~~~~
+
+This class describes an alignment. An object of this class has the following
+properties:
+
+* **ctg**: name of the reference sequence the query is mapped to
+
+* **ctg_len**: total length of the reference sequence
+
+* **r_st** and **r_en**: start and end positions on the reference
+
+* **q_st** and **q_en**: start and end positions on the query
+
+* **strand**: +1 if on the forward strand; -1 if on the reverse strand
+
+* **mapq**: mapping quality
+
+* **NM**: number of mismatches and gaps in the alignment
+
+* **blen**: length of the alignment, including both alignment matches and gaps
+
+* **trans_strand**: transcript strand. +1 if on the forward strand; -1 if on the
+  reverse strand; 0 if unknown
+
+* **is_primary**: if the alignment is primary (typically the best and the first
+  to generate)
+
+* **cigar_str**: CIGAR string
+
+* **cigar**: CIGAR returned as an array of shape :code:`(n_cigar,2)`. The two
+  numbers give the length and the operator of each CIGAR operation.
+
+An :code:`Alignment` object can be converted to a string with :code:`str()` in
+the following format:
+
+::
+
+	q_st  q_en  strand  ctg  ctg_len  r_st  r_en  blen-NM  blen  mapq  cg:Z:cigar_str
+
+It is effectively the PAF format without the QueryName and QueryLength columns
+(the first two columns in PAF).
+
+Function mappy.fastx_read
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. code:: python
+
+	mappy.fastx_read(fn)
+
+This generator function opens a FASTA/FASTQ file and *yields* a
+:code:`(name,seq,qual)` tuple for each sequence entry. The input file may be
+optionally gzip'd.

python/cmappy.h

@@ -0,0 +1,70 @@
+#ifndef CMAPPY_H
+#define CMAPPY_H
+
+#include <stdlib.h>
+#include <string.h>
+#include <zlib.h>
+#include "minimap.h"
+#include "kseq.h"
+KSEQ_DECLARE(gzFile)
+
+typedef struct {
+	const char *ctg;
+	int32_t ctg_start, ctg_end;
+	int32_t qry_start, qry_end;
+	int32_t blen, NM, ctg_len;
+	uint8_t mapq, is_primary;
+	int8_t strand, trans_strand;
+	int32_t n_cigar32;
+	uint32_t *cigar32;
+} mm_hitpy_t;
+
+static inline void mm_reg2hitpy(const mm_idx_t *mi, mm_reg1_t *r, mm_hitpy_t *h)
+{
+	h->ctg = mi->seq[r->rid].name;
+	h->ctg_len = mi->seq[r->rid].len;
+	h->ctg_start = r->rs, h->ctg_end = r->re;
+	h->qry_start = r->qs, h->qry_end = r->qe;
+	h->strand = r->rev? -1 : 1;
+	h->mapq = r->mapq;
+	h->blen = r->p->blen;
+	h->NM = r->p->n_diff;
+	h->trans_strand = r->p->trans_strand == 1? 1 : r->p->trans_strand == 2? -1 : 0;
+	h->is_primary = (r->id == r->parent);
+	h->n_cigar32 = r->p->n_cigar;
+	h->cigar32 = r->p->cigar;
+}
+
+static inline void mm_free_reg1(mm_reg1_t *r)
+{
+	free(r->p);
+}
+
+static inline kseq_t *mm_fastx_open(const char *fn)
+{
+	gzFile fp;
+	fp = fn && strcmp(fn, "-") != 0? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
+	return kseq_init(fp);
+}
+
+static inline void mm_fastx_close(kseq_t *ks)
+{
+	gzFile fp;
+	fp = ks->f->f;
+	kseq_destroy(ks);
+	gzclose(fp);
+}
+
+static inline int mm_verbose_level(int v)
+{
+	if (v >= 0) mm_verbose = v;
+	return mm_verbose;
+}
+
+static inline void mm_reset_timer(void)
+{
+	extern double realtime(void);
+	mm_realtime0 = realtime();
+}
+
+#endif

python/cmappy.pxd

@@ -0,0 +1,112 @@
+from libc.stdint cimport int8_t, uint8_t, int32_t, int64_t, uint32_t, uint64_t
+
+cdef extern from "minimap.h":
+	#
+	# Options
+	#
+	ctypedef struct mm_idxopt_t:
+		short k, w, is_hpc, bucket_bits
+		int mini_batch_size
+		uint64_t batch_size
+
+	ctypedef struct mm_mapopt_t:
+		int sdust_thres
+		int flag
+		int bw
+		int max_gap, max_gap_ref
+		int max_chain_skip
+		int min_cnt
+		int min_chain_score
+		float mask_level
+		float pri_ratio
+		int best_n
+		int max_join_long, max_join_short
+		int min_join_flank_sc
+		int a, b, q, e, q2, e2
+		int noncan
+		int zdrop
+		int min_dp_max
+		int min_ksw_len
+		float mid_occ_frac
+		int32_t mid_occ
+		int mini_batch_size
+
+	int mm_set_opt(char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
+	int mm_verbose
+
+	#
+	# Indexing
+	#
+	ctypedef struct mm_idx_seq_t:
+		char *name
+		uint64_t offset
+		uint32_t len
+
+	ctypedef struct mm_idx_bucket_t:
+		pass
+
+	ctypedef struct mm_idx_t:
+		int32_t b, w, k, is_hpc
+		uint32_t n_seq
+		mm_idx_seq_t *seq
+		uint32_t *S
+		mm_idx_bucket_t *B
+		void *km
+
+	ctypedef struct mm_idx_reader_t:
+		pass
+
+	mm_idx_reader_t *mm_idx_reader_open(const char *fn, const mm_idxopt_t *opt, const char *fn_out)
+	mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads)
+	void mm_idx_reader_close(mm_idx_reader_t *r)
+	void mm_idx_destroy(mm_idx_t *mi)
+	void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)
+
+	#
+	# Mapping (key struct defined in cmappy.h below)
+	#
+	ctypedef struct mm_reg1_t:
+		pass
+
+	ctypedef struct mm_tbuf_t:
+		pass
+
+	mm_tbuf_t *mm_tbuf_init()
+	void mm_tbuf_destroy(mm_tbuf_t *b)
+	mm_reg1_t *mm_map(const mm_idx_t *mi, int l_seq, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *name)
+
+#
+# Helper header (because it is hard to expose mm_reg1_t with Cython)
+#
+cdef extern from "cmappy.h":
+	ctypedef struct mm_hitpy_t:
+		const char *ctg
+		int32_t ctg_start, ctg_end
+		int32_t qry_start, qry_end
+		int32_t blen, NM, ctg_len
+		uint8_t mapq, is_primary
+		int8_t strand, trans_strand
+		int32_t n_cigar32
+		uint32_t *cigar32
+
+	void mm_reg2hitpy(const mm_idx_t *mi, mm_reg1_t *r, mm_hitpy_t *h)
+	void mm_free_reg1(mm_reg1_t *r)
+
+	ctypedef struct kstring_t:
+		unsigned l, m
+		char *s
+
+	ctypedef struct kstream_t:
+		pass
+
+	ctypedef struct kseq_t:
+		kstring_t name, comment, seq, qual
+		int last_char
+		kstream_t *f
+
+	kseq_t *mm_fastx_open(const char *fn)
+	void mm_fastx_close(kseq_t *ks)
+	int kseq_read(kseq_t *seq)
+
+	int mm_verbose_level(int v)
+	void mm_reset_timer()

python/mappy.pyx

@@ -0,0 +1,154 @@
+from libc.stdint cimport uint8_t, int8_t
+from libc.stdlib cimport free
+cimport cmappy
+
+cmappy.mm_reset_timer()
+
+cdef class Alignment:
+	cdef int _ctg_len, _r_st, _r_en
+	cdef int _q_st, _q_en
+	cdef int _NM, _blen
+	cdef int8_t _strand, _trans_strand
+	cdef uint8_t _mapq, _is_primary
+	cdef _ctg, _cigar # these are python objects
+
+	def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, blen, NM, trans_strand):
+		self._ctg, self._ctg_len, self._r_st, self._r_en = str(ctg), cl, cs, ce
+		self._strand, self._q_st, self._q_en = strand, qs, qe
+		self._NM, self._blen = NM, blen
+		self._mapq = mapq
+		self._cigar = cigar
+		self._is_primary = is_primary
+		self._trans_strand = trans_strand
+
+	@property
+	def ctg(self): return self._ctg
+
+	@property
+	def ctg_len(self): return self._ctg_len
+
+	@property
+	def r_st(self): return self._r_st
+
+	@property
+	def r_en(self): return self._r_en
+
+	@property
+	def strand(self): return self.strand
+
+	@property
+	def trans_strand(self): return self._trans_strand
+
+	@property
+	def NM(self): return self._NM
+
+	@property
+	def is_primary(self): return (self._is_primary != 0)
+
+	@property
+	def q_st(self): return self._q_st
+
+	@property
+	def q_en(self): return self._q_en
+
+	@property
+	def mapq(self): return self._mapq
+
+	@property
+	def cigar(self): return self._cigar
+
+	@property
+	def cigar_str(self):
+		return "".join(map(lambda x: str(x[0]) + 'MIDNSH'[x[1]], self._cigar))
+
+	def __str__(self):
+		if self._strand > 0: strand = '+'
+		elif self._strand < 0: strand = '-'
+		else: strand = '?'
+		if self._is_primary != 0: tp = 'tp:A:P'
+		else: tp = 'tp:A:S'
+		if self._trans_strand > 0: ts = 'ts:A:+'
+		elif self._trans_strand < 0: ts = 'ts:A:-'
+		else: ts = 'ts:A:.'
+		return "\t".join([str(self._q_st), str(self._q_en), strand, self._ctg, str(self._ctg_len), str(self._r_st), str(self._r_en),
+				str(self._blen - self._NM), str(self._blen), str(self._mapq), tp, ts, "cg:Z:" + self.cigar_str])
+
+cdef class ThreadBuffer:
+	cdef cmappy.mm_tbuf_t *_b
+
+	def __cinit__(self):
+		self._b = cmappy.mm_tbuf_init()
+
+	def __dealloc__(self):
+		cmappy.mm_tbuf_destroy(self._b)
+
+cdef class Aligner:
+	cdef cmappy.mm_idx_t *_idx
+	cdef cmappy.mm_idxopt_t idx_opt
+	cdef cmappy.mm_mapopt_t map_opt
+
+	def __cinit__(self, fn_idx_in, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, best_n=None, n_threads=3, fn_idx_out=None):
+		cmappy.mm_set_opt(NULL, &self.idx_opt, &self.map_opt) # set the default options
+		if preset is not None:
+			cmappy.mm_set_opt(str.encode(preset), &self.idx_opt, &self.map_opt) # apply preset
+		self.map_opt.flag |= 4 # always perform alignment
+		self.idx_opt.batch_size = 0x7fffffffffffffffL # always build a uni-part index
+		if k is not None: self.idx_opt.k = k
+		if w is not None: self.idx_opt.w = w
+		if min_cnt is not None: self.map_opt.min_cnt = min_cnt
+		if min_chain_score is not None: self.map_opt.min_chain_score = min_chain_score
+		if min_dp_score is not None: self.map_opt.min_dp_max = min_dp_score
+		if bw is not None: self.map_opt.bw = bw
+		if best_n is not None: self.best_n = best_n
+
+		cdef cmappy.mm_idx_reader_t *r;
+		if fn_idx_out is None:
+			r = cmappy.mm_idx_reader_open(str.encode(fn_idx_in), &self.idx_opt, NULL)
+		else:
+			r = cmappy.mm_idx_reader_open(str.encode(fn_idx_in), &self.idx_opt, fn_idx_out)
+		if r is not NULL:
+			self._idx = cmappy.mm_idx_reader_read(r, n_threads) # NB: ONLY read the first part
+			cmappy.mm_idx_reader_close(r)
+			cmappy.mm_mapopt_update(&self.map_opt, self._idx)
+
+	def __dealloc__(self):
+		if self._idx is not NULL:
+			cmappy.mm_idx_destroy(self._idx)
+
+	def __bool__(self):
+		return (self._idx != NULL)
+
+	def map(self, seq, buf=None):
+		cdef cmappy.mm_reg1_t *regs
+		cdef cmappy.mm_hitpy_t h
+		cdef ThreadBuffer b
+		cdef int n_regs
+
+		if self._idx is NULL: return None
+		if buf is None: b = ThreadBuffer()
+		else: b = buf
+		regs = cmappy.mm_map(self._idx, len(seq), str.encode(seq), &n_regs, b._b, &self.map_opt, NULL)
+
+		for i in range(n_regs):
+			cmappy.mm_reg2hitpy(self._idx, &regs[i], &h)
+			cigar = []
+			for k in range(h.n_cigar32):
+				c = h.cigar32[k]
+				cigar.append([c>>4, c&0xf])
+			yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.blen, h.NM, h.trans_strand)
+			cmappy.mm_free_reg1(&regs[i])
+		free(regs)
+
+def fastx_read(fn):
+	cdef cmappy.kseq_t *ks
+	ks = cmappy.mm_fastx_open(str.encode(fn))
+	if ks is NULL: return None
+	while cmappy.kseq_read(ks) >= 0:
+		if ks.qual.l > 0: qual = str(ks.qual.s)
+		else: qual = None
+		yield str(ks.name.s), str(ks.seq.s), qual
+	cmappy.mm_fastx_close(ks)
+
+def verbose(v=None):
+	if v is None: v = -1
+	return cmappy.mm_verbose_level(v)

python/minimap2.py

@@ -0,0 +1,35 @@
+#!/usr/bin/env python
+
+import sys, getopt
+import mappy as mp
+
+def main(argv):
+	opts, args = getopt.getopt(argv[1:], "x:n:m:k:w:r:")
+	if len(args) < 2:
+		print("Usage: minimap2.py [options] <ref.fa>|<ref.mmi> <query.fq>")
+		print("Options:")
+		print("  -x STR      preset: sr, map-pb, map-ont, asm5, asm10 or splice")
+		print("  -n INT      mininum number of minimizers")
+		print("  -m INT      mininum chaining score")
+		print("  -k INT      k-mer length")
+		print("  -w INT      minimizer window length")
+		print("  -r INT      band width")
+		sys.exit(1)
+
+	preset, min_cnt, min_sc, k, w, bw = None, None, None, None, None, None
+	for opt, arg in opts:
+		if opt == '-x': preset = arg
+		elif opt == '-n': min_cnt = int(arg)
+		elif opt == '-m': min_chain_score = int(arg)
+		elif opt == '-r': bw = int(arg)
+		elif opt == '-k': k = int(arg)
+		elif opt == '-w': w = int(arg)
+
+	a = mp.Aligner(args[0], preset=preset, min_cnt=min_cnt, min_chain_score=min_sc, k=k, w=w, bw=bw)
+	if not a: raise Exception("ERROR: failed to load/build index file '{}'".format(args[0]))
+	for name, seq, qual in mp.fastx_read(args[1]): # read one sequence
+		for h in a.map(seq): # traverse hits
+			print('{}\t{}\t{}'.format(name, len(seq), h))
+
+if __name__ == "__main__":
+	main(sys.argv)

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)
 

setup.py

@@ -0,0 +1,55 @@
+try:
+	from setuptools import setup, Extension
+except ImportError:
+	from distutils.core import setup
+	from distutils.extension import Extension
+
+cmdclass = {}
+
+try:
+	from Cython.Build import build_ext
+except ImportError: # without Cython
+	module_src = 'python/mappy.c'
+else: # with Cython
+	module_src = 'python/mappy.pyx'
+	cmdclass['build_ext'] = build_ext
+
+import sys
+sys.path.append('python')
+
+def readme():
+    with open('python/README.rst') as f:
+        return f.read()
+
+setup(
+	name = 'mappy',
+	version = '2.2',
+	url = 'https://github.com/lh3/minimap2',
+	description = 'Minimap2 python binding',
+	long_description = readme(),
+	author = 'Heng Li',
+	author_email = 'lh3@me.com',
+	license = 'MIT',
+	keywords = 'sequence-alignment',
+	scripts = ['python/minimap2.py'],
+    ext_modules = [Extension('mappy',
+		sources = [module_src, 'align.c', 'bseq.c', 'chain.c', 'format.c', 'hit.c', 'index.c',
+				   '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'],
+		depends = ['minimap.h', 'bseq.h', 'kalloc.h', 'kdq.h', 'khash.h', 'kseq.h', 'ksort.h',
+				   'ksw2.h', 'kthread.h', 'kvec.h', 'mmpriv.h', 'sdust.h',
+				   'python/cmappy.h', 'python/cmappy.pxd'],
+		extra_compile_args = ['-msse4'], # WARNING: ancient x86_64 CPUs don't have SSE4
+		include_dirs = ['.'],
+		libraries = ['z', 'm', 'pthread'])],
+	classifiers = [
+		'Development Status :: 4 - Beta',
+		'License :: OSI Approved :: MIT License',
+		'Operating System :: POSIX',
+		'Programming Language :: C',
+		'Programming Language :: Cython',
+		'Programming Language :: Python :: 2.7',
+		'Programming Language :: Python :: 3',
+		'Intended Audience :: Science/Research',
+		'Topic :: Scientific/Engineering :: Bio-Informatics'],
+	cmdclass = cmdclass)

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
+	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
-	buf = (mm128_t*)alloca(w * 16);
 	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
+difference-based formulation that lifted this limitation.
-gap cost, define
+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
+where $z_{ij}$ is a temporary variable that does not need to be stored.
-important property of Eq.~(\ref{eq:suzuki}) is that all values are bounded. To
+
-see that,
+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
+On the condition that $q+e<\tilde{q}+\tilde{e}$ and $e>\tilde{e}$, the initial
-condition of the equation above is
+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}.
+This strategy is first used in BWA-MEM. It is similar to X-drop employed in
-However, unlike X-drop, it would not break the alignment in the presence of a
+BLAST~\citep{Altschul:1997vn}, but unlike X-drop, it would not break the
-single long gap.
+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
+GRCh38.  33,088 $\ge$1000bp reads were simulated using pbsim~\citep{Ono:2013aa}
-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}
 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
+run under the default setting for SMRT reads.  (a) ROC-like curve. Alignments
-mapping quality 60, so is not shown in the figure. It mapped nearly all reads
+are sorted by mapping quality in the descending order. For each mapping quality
-with 4.1\% of alignments being wrong, less accurate than others.}\label{fig:eval}
+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
+grateful to early minimap2 testers who have greatly helped to suggest features
-issues.
+and to fix various issues.
 
 \bibliography{minimap2}