Release minimap2-2.2 (r409)

can't get cython working...

.gitignore

@@ -4,3 +4,5 @@
 *.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,4 +1,3 @@
-CC=			gcc
 CFLAGS=		-g -Wall -O2 -Wc++-compat
 CPPFLAGS=	-DHAVE_KALLOC
 INCLUDES=
@@ -56,7 +55,7 @@ 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)

NEWS.md

@@ -1,3 +1,32 @@
+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)
 -------------------------------------
 

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

align.c

@@ -143,6 +143,10 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 		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);
 	else if (opt->q == opt->q2 && opt->e == opt->e2)

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,53 +11,52 @@ 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");
 		return 1;
 	}
-	
+
 	// open query file for reading; you may use your favorite FASTA/Q parser
 	gzFile f = gzopen(argv[2], "r");
 	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];
-		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]);
-			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');
-			free(r->p);
 		}
-		free(reg);
+		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

@@ -113,7 +113,7 @@ static int __getopt_long_core(int argc, char *const *argv, const char *optstring
 		 (argv[optind][1] == '-' && argv[optind][2])))
 	{
 		int colon = optstring[optstring[0]=='+'||optstring[0]=='-']==':';
-		int i, cnt, match;
+		int i, cnt, match = -1;
 		char *opt;
 		for (cnt=i=0; longopts[i].name; i++) {
 			const char *name = longopts[i].name;

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

@@ -21,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;
@@ -104,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);
@@ -314,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;
 }
@@ -429,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;
+}

main.c

@@ -6,7 +6,7 @@
 #include "mmpriv.h"
 #include "getopt.h"
 
-#define MM_VERSION "2.1.1-r341"
+#define MM_VERSION "2.2-r409"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -25,17 +25,18 @@ void liftrlimit() {}
 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' },
@@ -61,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);
@@ -98,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
@@ -114,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;
@@ -134,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;
 
@@ -178,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");
@@ -208,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");
@@ -220,55 +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 (mi == 0) break;
-			if (idx_par_set && mm_verbose >= 2 && (mi->k != k || mi->w != w || mi->is_hpc != is_hpc))
-				fprintf(stderr, "[WARNING] \033[1;31mIndexing parameters on the command line (-k/-w/-H) overridden by parameters in the prebuilt index.\033[0m\n");
-		} else {
-			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,7 +352,7 @@ 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));
@@ -388,7 +363,7 @@ int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int
 		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 "6 September 2017" "minimap2-2.1.1-r341" "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,6 +1,6 @@
 #include "minimap.h"
 
-int mm_verbose = 3;
+int mm_verbose = 1;
 int mm_dbg_flag = 0;
 double mm_realtime0;
 
@@ -90,7 +90,7 @@ double cputime()
 #include <sys/resource.h>
 #include <sys/time.h>
 
-double cputime()
+double cputime(void)
 {
 	struct rusage r;
 	getrusage(RUSAGE_SELF, &r);
@@ -98,7 +98,7 @@ double cputime()
 }
 #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)

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)