debugging code

We may use a large --end-bonus to mimic end-to-end alignment. In the short-read
mode, the candidate alignment region may be out of the band, which leads to
truncated alignment.

LICENSE.txt

@@ -1,6 +1,7 @@
 The MIT License
 
-Copyright (c) 2017 Broad Institute, Inc.
+Copyright (c) 2018-     Dana-Farber Cancer Institute
+              2017-2018 Broad Institute, Inc.
 
 Permission is hereby granted, free of charge, to any person obtaining
 a copy of this software and associated documentation files (the

MANIFEST.in

@@ -1,9 +1,9 @@
 include *.h
 include Makefile
 include ksw2_dispatch.c
-include getopt.c
 include main.c
 include README.md
+include sse2neon/emmintrin.h
 include python/mappy.c
 include python/cmappy.h
 include python/cmappy.pxd

Makefile

@@ -1,7 +1,7 @@
-CFLAGS=		-g -Wall -O2 -Wc++-compat
+CFLAGS=		-g -Wall -O2 -Wc++-compat #-Wextra
 CPPFLAGS=	-DHAVE_KALLOC
 INCLUDES=
-OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o chain.o align.o hit.o map.o format.o pe.o esterr.o ksw2_ll_sse.o
+OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o chain.o align.o hit.o map.o format.o pe.o esterr.o splitidx.o ksw2_ll_sse.o
 PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite
 LIBS=		-lm -lz -lpthread
@@ -12,10 +12,14 @@ ifeq ($(sse2only),) # if sse2only is not defined
 else                # if sse2only is defined
 	OBJS+=ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_exts2_sse.o
 endif
-else                # if arm_neon is defined
-    OBJS+=ksw2_extz2_neon.o ksw2_extd2_neon.o ksw2_exts2_neon.o
-    CFLAGS+=-D_FILE_OFFSET_BITS=64 -mfpu=neon -fsigned-char
+else				# if arm_neon is defined
+	OBJS+=ksw2_extz2_neon.o ksw2_extd2_neon.o ksw2_exts2_neon.o
     INCLUDES+=-Isse2neon
+ifeq ($(aarch64),)	#if aarch64 is not defined
+	CFLAGS+=-D_FILE_OFFSET_BITS=64 -mfpu=neon -fsigned-char
+else				#if aarch64 is defined
+	CFLAGS+=-D_FILE_OFFSET_BITS=64 -fsigned-char
+endif
 endif
 
 .PHONY:all extra clean depend
@@ -28,8 +32,8 @@ all:$(PROG)
 
 extra:all $(PROG_EXTRA)
 
-minimap2:main.o getopt.o libminimap2.a
-		$(CC) $(CFLAGS) main.o getopt.o -o $@ -L. -lminimap2 $(LIBS)
+minimap2:main.o libminimap2.a
+		$(CC) $(CFLAGS) main.o -o $@ -L. -lminimap2 $(LIBS)
 
 minimap2-lite:example.o libminimap2.a
 		$(CC) $(CFLAGS) $< -o $@ -L. -lminimap2 $(LIBS)
@@ -37,8 +41,8 @@ minimap2-lite:example.o libminimap2.a
 libminimap2.a:$(OBJS)
 		$(AR) -csru $@ $(OBJS)
 
-sdust:sdust.c getopt.o kalloc.o kalloc.h kdq.h kvec.h kseq.h sdust.h
-		$(CC) -D_SDUST_MAIN $(CFLAGS) $< getopt.o kalloc.o -o $@ -lz
+sdust:sdust.c kalloc.o kalloc.h kdq.h kvec.h kseq.h ketopt.h sdust.h
+		$(CC) -D_SDUST_MAIN $(CFLAGS) $< kalloc.o -o $@ -lz
 
 # SSE-specific targets on x86/x86_64
 
@@ -95,7 +99,6 @@ chain.o: minimap.h mmpriv.h bseq.h kalloc.h
 esterr.o: mmpriv.h minimap.h bseq.h
 example.o: minimap.h kseq.h
 format.o: kalloc.h mmpriv.h minimap.h bseq.h
-getopt.o: getopt.h
 hit.o: mmpriv.h minimap.h bseq.h kalloc.h khash.h
 index.o: kthread.h bseq.h minimap.h mmpriv.h kvec.h kalloc.h khash.h
 kalloc.o: kalloc.h
@@ -104,11 +107,12 @@ ksw2_exts2_sse.o: ksw2.h kalloc.h
 ksw2_extz2_sse.o: ksw2.h kalloc.h
 ksw2_ll_sse.o: ksw2.h kalloc.h
 kthread.o: kthread.h
-main.o: bseq.h minimap.h mmpriv.h getopt.h
+main.o: bseq.h minimap.h mmpriv.h ketopt.h
 map.o: kthread.h kvec.h kalloc.h sdust.h mmpriv.h minimap.h bseq.h khash.h
 map.o: ksort.h
 misc.o: mmpriv.h minimap.h bseq.h ksort.h
 options.o: mmpriv.h minimap.h bseq.h
 pe.o: mmpriv.h minimap.h bseq.h kvec.h kalloc.h ksort.h
-sdust.o: kalloc.h kdq.h kvec.h sdust.h
+sdust.o: kalloc.h kdq.h kvec.h ketopt.h sdust.h
 sketch.o: kvec.h kalloc.h mmpriv.h minimap.h bseq.h
+splitidx.o: mmpriv.h minimap.h bseq.h

NEWS.md

@@ -1,3 +1,105 @@
+Release 2.13-r850 (11 October 2018)
+-----------------------------------
+
+Changes to minimap2:
+
+ * Fixed wrongly formatted SAM when -L is in use (#231 and #233).
+
+ * Fixed an integer overflow in rare cases.
+
+ * Added --hard-mask-level to fine control split alignments (#244).
+
+ * Made --MD work with spliced alignment (#139).
+
+ * Replaced musl's getopt with ketopt for portability.
+
+ * Log peak memory usage on exit.
+
+This release should produce alignments identical to v2.12 and v2.11.
+
+(2.13: 11 October 2018, r850)
+
+
+
+Release 2.12-r827 (6 August 2018)
+---------------------------------
+
+Changes to minimap2:
+
+ * Added option --split-prefix to write proper alignments (correct mapping
+   quality and clustered query sequences) given a multi-part index (#141 and
+   #189; mostly by @hasindu2008).
+
+ * Fixed a memory leak when option -y is in use.
+
+Changes to mappy:
+
+ * Support the MD/cs tag (#183 and #203).
+
+ * Allow mappy to index a single sequence, to add extra flags and to change the
+   scoring system.
+
+Minimap2 should produce alignments identical to v2.11.
+
+(2.12: 6 August 2018, r827)
+
+
+
+Release 2.11-r797 (20 June 2018)
+--------------------------------
+
+Changes to minimap2:
+
+ * Improved alignment accuracy in low-complexity regions for SV calling. Thank
+   @armintoepfer for multiple offline examples.
+
+ * Added option --eqx to encode sequence match/mismatch with the =/X CIGAR
+   operators (#156, #157 and #175).
+
+ * When compiled with VC++, minimap2 generated wrong alignments due to a
+   comparison between a signed integer and an unsigned integer (#184). Also
+   fixed warnings reported by "clang -Wextra".
+
+ * Fixed incorrect anchor filtering due to a missing 64- to 32-bit cast.
+
+ * Fixed incorrect mapping quality for inversions (#148).
+
+ * Fixed incorrect alignment involving ambiguous bases (#155).
+
+ * Fixed incorrect presets: option `-r 2000` is intended to be used with
+   ava-ont, not ava-pb. The bug was introduced in 2.10.
+
+ * Fixed a bug when --for-only/--rev-only is used together with --sr or
+   --heap-sort=yes (#166).
+
+ * Fixed option -Y that was not working in the previous releases.
+
+ * Added option --lj-min-ratio to fine control the alignment of long gaps
+   found by the "long-join" heuristic (#128).
+
+ * Exposed `mm_idx_is_idx`, `mm_idx_load` and `mm_idx_dump` C APIs (#177).
+   Also fixed a bug when indexing without reference names (this feature is not
+   exposed to the command line).
+
+Changes to mappy:
+
+ * Added `__version__` (#165).
+
+ * Exposed the maximum fragment length parameter to mappy (#174).
+
+Changes to paftools:
+
+ * Don't crash when there is no "cg" tag (#153).
+
+ * Fixed wrong coverage report by "paftools.js call" (#145).
+
+This version may produce slightly different base-level alignment. The overall
+alignment statistics should remain similar.
+
+(2.11: 20 June 2018, r797)
+
+
+
 Release 2.10-r761 (27 March 2018)
 ---------------------------------
 

README.md

@@ -14,9 +14,11 @@ cd minimap2 && make
 # use presets (no test data)
 ./minimap2 -ax map-pb ref.fa pacbio.fq.gz > aln.sam       # PacBio genomic reads
 ./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam         # Oxford Nanopore genomic reads
+./minimap2 -ax asm20 ref.fa pacbio-ccs.fq.gz > aln.sam    # PacBio CCS genomic reads
 ./minimap2 -ax sr ref.fa read1.fa read2.fa > aln.sam      # short genomic paired-end reads
-./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads
-./minimap2 -ax splice -k14 -uf ref.fa reads.fa > aln.sam  # Nanopore Direct RNA-seq
+./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads (strand unknown)
+./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam  # noisy Nanopore Direct RNA-seq
+./minimap2 -ax splice -uf -C5 ref.fa query.fa > aln.sam   # Final PacBio Iso-seq or traditional cDNA
 ./minimap2 -cx asm5 asm1.fa asm2.fa > aln.paf             # intra-species asm-to-asm alignment
 ./minimap2 -x ava-pb reads.fa reads.fa > overlaps.paf     # PacBio read overlap
 ./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf    # Nanopore read overlap
@@ -61,15 +63,16 @@ mainstream long-read mappers such as BLASR, BWA-MEM, NGMLR and GMAP. It is more
 accurate on simulated long reads and produces biologically meaningful alignment
 ready for downstream analyses. For >100bp Illumina short reads, minimap2 is
 three times as fast as BWA-MEM and Bowtie2, and as accurate on simulated data.
-Detailed evaluations are available from the [minimap2 preprint][preprint].
+Detailed evaluations are available from the [minimap2 paper][doi] or the
+[preprint][preprint].
 
 ### <a name="install"></a>Installation
 
 Minimap2 is optimized for x86-64 CPUs. You can acquire precompiled binaries from
 the [release page][release] with:
 ```sh
-curl -L https://github.com/lh3/minimap2/releases/download/v2.10/minimap2-2.10_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.10_x64-linux/minimap2
+curl -L https://github.com/lh3/minimap2/releases/download/v2.13/minimap2-2.13_x64-linux.tar.bz2 | tar -jxvf -
+./minimap2-2.13_x64-linux/minimap2
 ```
 If you want to compile from the source, you need to have a C compiler, GNU make
 and zlib development files installed. Then type `make` in the source code
@@ -77,7 +80,7 @@ directory to compile. If you see compilation errors, try `make sse2only=1`
 to disable SSE4 code, which will make minimap2 slightly slower.
 
 Minimap2 also works with ARM CPUs supporting the NEON instruction sets. To
-compile, use `make arm_neon=1`.
+compile for 32 bit ARM architectures (such as ARMv7), use `make arm_neon=1`. To compile for for 64 bit ARM architectures (such as ARMv8), use `make arm_neon=1 aarch64=1`.
 
 ### <a name="general"></a>General usage
 
@@ -253,7 +256,9 @@ similar to the `MD` SAM tag but is standalone and easier to parse.
 If `--cs=long` is used, the `cs` string also contains identical sequences in
 the alignment. The above example will become
 `=CGATCG-ata=AATAGAGTAG+gtc=GAAT*at=GCA`. The long form of `cs` encodes both
-reference and query sequences in one string.
+reference and query sequences in one string. The `cs` tag also encodes intron
+positions and splicing signals (see the [minimap2 manpage][manpage-cs] for
+details).
 
 #### <a name="paftools"></a>Working with the PAF format
 
@@ -316,9 +321,10 @@ There is not a specific mailing list for the time being.
 
 ### <a name="cite"></a>Citing minimap2
 
-If you use minimap2 in your work, please consider to cite:
+If you use minimap2 in your work, please cite:
 
-> Li, H. (2017). Minimap2: fast pairwise alignment for long nucleotide sequences. [arXiv:1708.01492][preprint]
+> Li, H. (2018). Minimap2: pairwise alignment for nucleotide sequences.
+> Bioinformatics. [doi:10.1093/bioinformatics/bty191][doi]
 
 ## <a name="dguide"></a>Developers' Guide
 
@@ -365,3 +371,5 @@ mappy` or [from BioConda][mappyconda] via `conda install -c bioconda mappy`.
 [issue]: https://github.com/lh3/minimap2/issues
 [k8]: https://github.com/attractivechaos/k8
 [manpage]: https://lh3.github.io/minimap2/minimap2.html
+[manpage-cs]: https://lh3.github.io/minimap2/minimap2.html#10
+[doi]: https://doi.org/10.1093/bioinformatics/bty191

align.c

@@ -6,18 +6,19 @@
 #include "mmpriv.h"
 #include "ksw2.h"
 
-static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b)
+static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t sc_ambi)
 {
 	int i, j;
 	a = a < 0? -a : a;
 	b = b > 0? -b : b;
+	sc_ambi = sc_ambi > 0? -sc_ambi : sc_ambi;
 	for (i = 0; i < m - 1; ++i) {
 		for (j = 0; j < m - 1; ++j)
 			mat[i * m + j] = i == j? a : b;
-		mat[i * m + m - 1] = 0;
+		mat[i * m + m - 1] = sc_ambi;
 	}
 	for (j = 0; j < m; ++j)
-		mat[(m - 1) * m + j] = 0;
+		mat[(m - 1) * m + j] = sc_ambi;
 }
 
 static inline void mm_seq_rev(uint32_t len, uint8_t *seq)
@@ -87,12 +88,14 @@ static int mm_test_zdrop(void *km, const mm_mapopt_t *opt, const uint8_t *qseq,
 	return max_zdrop > opt->zdrop? 1 : 0;
 }
 
-static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int *qshift, int *tshift)
+static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int *qshift, int *tshift, int left_aln)
 {
 	mm_extra_t *p = r->p;
-	int32_t k, toff = 0, qoff = 0, to_shrink = 0;
+	int32_t toff = 0, qoff = 0, to_shrink = 0;
+	uint32_t k;
 	*qshift = *tshift = 0;
 	if (p->n_cigar <= 1) return;
+	if (!left_aln) goto end_left_aln;
 	for (k = 0; k < p->n_cigar; ++k) { // indel left alignment
 		uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4;
 		if (len == 0) to_shrink = 1;
@@ -121,6 +124,7 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 		}
 	}
 	assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
+end_left_aln:
 	if (to_shrink) { // squeeze out zero-length operations
 		int32_t l = 0;
 		for (k = 0; k < p->n_cigar; ++k) // squeeze out zero-length operations
@@ -145,13 +149,13 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 	}
 }
 
-static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e)
+static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e, int left_aln)
 {
-	uint32_t k, l, toff = 0, qoff = 0;
-	int32_t s = 0, max = 0, qshift, tshift;
+	uint32_t k, l;
+	int32_t s = 0, max = 0, qshift, tshift, toff = 0, qoff = 0;
 	mm_extra_t *p = r->p;
 	if (p == 0) return;
-	mm_fix_cigar(r, qseq, tseq, &qshift, &tshift);
+	mm_fix_cigar(r, qseq, tseq, &qshift, &tshift, left_aln);
 	qseq += qshift, tseq += tshift; // qseq and tseq may be shifted due to the removal of leading I/D
 	r->blen = r->mlen = 0;
 	for (k = 0; k < p->n_cigar; ++k) {
@@ -197,12 +201,12 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) //
 	mm_extra_t *p;
 	if (n_cigar == 0) return;
 	if (r->p == 0) {
-		uint32_t capacity = n_cigar + sizeof(mm_extra_t);
+		uint32_t capacity = n_cigar + sizeof(mm_extra_t)/4;
 		kroundup32(capacity);
 		r->p = (mm_extra_t*)calloc(capacity, 4);
 		r->p->capacity = capacity;
-	} else if (r->p->n_cigar + n_cigar + sizeof(mm_extra_t) > r->p->capacity) {
-		r->p->capacity = r->p->n_cigar + n_cigar + sizeof(mm_extra_t);
+	} else if (r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4 > r->p->capacity) {
+		r->p->capacity = r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4;
 		kroundup32(r->p->capacity);
 		r->p = (mm_extra_t*)realloc(r->p, r->p->capacity * 4);
 	}
@@ -217,7 +221,79 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) //
 	}
 }
 
-static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
+static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq) // written by @armintoepfer
+{
+	uint32_t n_EQX = 0;
+	uint32_t k, l, m, cap, toff = 0, qoff = 0, n_M = 0;
+	mm_extra_t *p;
+	if (r->p == 0) return;
+	for (k = 0; k < r->p->n_cigar; ++k) {
+		uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4;
+		if (op == 0) {
+			while (len > 0) {
+				for (l = 0; l < len && qseq[qoff + l] == tseq[toff + l]; ++l) {} // run of "="; TODO: N<=>N is converted to "="
+				if (l > 0) { ++n_EQX; len -= l; toff += l; qoff += l; }
+
+				for (l = 0; l < len && qseq[qoff + l] != tseq[toff + l]; ++l) {} // run of "X"
+				if (l > 0) { ++n_EQX; len -= l; toff += l; qoff += l; }
+			}
+			++n_M;
+		} else if (op == 1) { // insertion
+			qoff += len;
+		} else if (op == 2) { // deletion
+			toff += len;
+		} else if (op == 3) { // intron
+			toff += len;
+		}
+	}
+	// update in-place if we can
+	if (n_EQX == n_M) {
+		for (k = 0; k < r->p->n_cigar; ++k) {
+			uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4;
+			if (op == 0) r->p->cigar[k] = len << 4 | 7;
+		}
+		return;
+	}
+	// allocate new storage
+	cap = r->p->n_cigar + (n_EQX - n_M) + sizeof(mm_extra_t);
+	kroundup32(cap);
+	p = (mm_extra_t*)calloc(cap, 4);
+	memcpy(p, r->p, sizeof(mm_extra_t));
+	p->capacity = cap;
+	// update cigar while copying
+	toff = qoff = m = 0;
+	for (k = 0; k < r->p->n_cigar; ++k) {
+		uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4;
+		if (op == 0) { // match/mismatch
+			while (len > 0) {
+				// match
+				for (l = 0; l < len && qseq[qoff + l] == tseq[toff + l]; ++l) {}
+				if (l > 0) p->cigar[m++] = l << 4 | 7;
+				len -= l;
+				toff += l, qoff += l;
+				// mismatch
+				for (l = 0; l < len && qseq[qoff + l] != tseq[toff + l]; ++l) {}
+				if (l > 0) p->cigar[m++] = l << 4 | 8;
+				len -= l;
+				toff += l, qoff += l;
+			}
+			continue;
+		} else if (op == 1) { // insertion
+			qoff += len;
+		} else if (op == 2) { // deletion
+			toff += len;
+		} else if (op == 3) { // intron
+			toff += len;
+		}
+		p->cigar[m++] = r->p->cigar[k];
+	}
+	p->n_cigar = m;
+	free(r->p);
+	r->p = p;
+}
+
+static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const int8_t *ob, const int8_t *mat,
+						  int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
 {
 	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
 		int i;
@@ -232,7 +308,7 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 	else if (opt->q == opt->q2 && opt->e == opt->e2)
 		ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, flag, ez);
 	else
-		ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez);
+		ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ob, ez);
 	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
 		int i;
 		fprintf(stderr, "score=%d, cigar=", ez->score);
@@ -268,20 +344,30 @@ static inline void mm_adjust_minier(const mm_idx_t *mi, uint8_t *const qseq0[2],
 	}
 }
 
-static void mm_filter_bad_seeds(void *km, int as1, int cnt1, mm128_t *a, int min_gap, int diff_thres, int max_ext_len, int max_ext_cnt)
+static int *collect_long_gaps(void *km, int as1, int cnt1, mm128_t *a, int min_gap, int *n_)
 {
-	int max_st, max_en, n, i, k, max, *K;
+	int i, n, *K;
+	*n_ = 0;
 	for (i = 1, n = 0; i < cnt1; ++i) { // count the number of gaps longer than min_gap
 		int gap = ((int32_t)a[as1 + i].y - a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - a[as1 + i - 1].x);
 		if (gap < -min_gap || gap > min_gap) ++n;
 	}
-	if (n <= 1) return;
+	if (n <= 1) return 0;
 	K = (int*)kmalloc(km, n * sizeof(int));
 	for (i = 1, n = 0; i < cnt1; ++i) { // store the positions of long gaps
 		int gap = ((int32_t)a[as1 + i].y - a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - a[as1 + i - 1].x);
 		if (gap < -min_gap || gap > min_gap)
 			K[n++] = i;
 	}
+	*n_ = n;
+	return K;
+}
+
+static void mm_filter_bad_seeds(void *km, int as1, int cnt1, mm128_t *a, int min_gap, int diff_thres, int max_ext_len, int max_ext_cnt)
+{
+	int max_st, max_en, n, i, k, max, *K;
+	K = collect_long_gaps(km, as1, cnt1, a, min_gap, &n);
+	if (K == 0) return;
 	max = 0, max_st = max_en = -1;
 	for (k = 0;; ++k) { // traverse long gaps
 		int gap, l, n_ins = 0, n_del = 0, qs, rs, max_diff = 0, max_diff_l = -1;
@@ -293,7 +379,7 @@ static void mm_filter_bad_seeds(void *km, int as1, int cnt1, mm128_t *a, int min
 			if (k == n) break;
 		}
 		i = K[k];
-		gap = ((int32_t)a[as1 + i].y - a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - a[as1 + i - 1].x);
+		gap = ((int32_t)a[as1 + i].y - (int32_t)a[as1 + i - 1].y) - (int32_t)(a[as1 + i].x - a[as1 + i - 1].x);
 		if (gap > 0) n_ins += gap;
 		else n_del += -gap;
 		qs = (int32_t)a[as1 + i - 1].y;
@@ -301,7 +387,7 @@ static void mm_filter_bad_seeds(void *km, int as1, int cnt1, mm128_t *a, int min
 		for (l = k + 1; l < n && l <= k + max_ext_cnt; ++l) {
 			int j = K[l], diff;
 			if ((int32_t)a[as1 + j].y - qs > max_ext_len || (int32_t)a[as1 + j].x - rs > max_ext_len) break;
-			gap = ((int32_t)a[as1 + j].y - (int32_t)a[as1 + j - 1].y) - (a[as1 + j].x - a[as1 + j - 1].x);
+			gap = ((int32_t)a[as1 + j].y - (int32_t)a[as1 + j - 1].y) - (int32_t)(a[as1 + j].x - a[as1 + j - 1].x);
 			if (gap > 0) n_ins += gap;
 			else n_del += -gap;
 			diff = n_ins + n_del - abs(n_ins - n_del);
@@ -314,6 +400,42 @@ static void mm_filter_bad_seeds(void *km, int as1, int cnt1, mm128_t *a, int min
 	kfree(km, K);
 }
 
+static void mm_filter_bad_seeds_alt(void *km, int as1, int cnt1, mm128_t *a, int min_gap, int max_ext)
+{
+	int n, k, *K;
+	K = collect_long_gaps(km, as1, cnt1, a, min_gap, &n);
+	if (K == 0) return;
+	for (k = 0; k < n;) {
+		int i = K[k], l;
+		int gap1 = ((int32_t)a[as1 + i].y - (int32_t)a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - (int32_t)a[as1 + i - 1].x);
+		int re1 = (int32_t)a[as1 + i].x;
+		int qe1 = (int32_t)a[as1 + i].y;
+		gap1 = gap1 > 0? gap1 : -gap1;
+		for (l = k + 1; l < n; ++l) {
+			int j = K[l], gap2, q_span_pre, rs2, qs2, m;
+			if ((int32_t)a[as1 + j].y - qe1 > max_ext || (int32_t)a[as1 + j].x - re1 > max_ext) break;
+			gap2 = ((int32_t)a[as1 + j].y - (int32_t)a[as1 + j - 1].y) - (int32_t)(a[as1 + j].x - a[as1 + j - 1].x);
+			q_span_pre = a[as1 + j - 1].y >> 32 & 0xff;
+			rs2 = (int32_t)a[as1 + j - 1].x + q_span_pre;
+			qs2 = (int32_t)a[as1 + j - 1].x + q_span_pre;
+			m = rs2 - re1 < qs2 - qe1? rs2 - re1 : qs2 - qe1;
+			gap2 = gap2 > 0? gap2 : -gap2;
+			if (m > gap1 + gap2) break;
+			re1 = (int32_t)a[as1 + j].x;
+			qe1 = (int32_t)a[as1 + j].y;
+			gap1 = gap2;
+		}
+		if (l > k + 1) {
+			int j, end = K[l - 1];
+			for (j = K[k]; j < end; ++j)
+				a[as1 + j].y |= MM_SEED_IGNORE;
+			a[as1 + end].y |= MM_SEED_LONG_JOIN;
+		}
+		k = l;
+	}
+	kfree(km, K);
+}
+
 static void mm_fix_bad_ends(const mm_reg1_t *r, const mm128_t *a, int bw, int min_match, int32_t *as, int32_t *cnt)
 {
 	int32_t i, l, m;
@@ -350,7 +472,7 @@ static void mm_fix_bad_ends(const mm_reg1_t *r, const mm128_t *a, int bw, int mi
 	}
 }
 
-static void mm_max_stretch(const mm_mapopt_t *opt, const mm_reg1_t *r, const mm128_t *a, int32_t *as, int32_t *cnt)
+static void mm_max_stretch(const mm_reg1_t *r, const mm128_t *a, int32_t *as, int32_t *cnt)
 {
 	int32_t i, score, max_score, len, max_i, max_len;
 
@@ -388,7 +510,7 @@ static int mm_seed_ext_score(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 	qe = (uint32_t)a->y + 1, qs = qe - q_span;
 	rs = rs - ext_len > 0? rs - ext_len : 0;
 	qs = qs - ext_len > 0? qs - ext_len : 0;
-	re = re + ext_len < mi->seq[rid].len? re + ext_len : mi->seq[rid].len;
+	re = re + ext_len < (int32_t)mi->seq[rid].len? re + ext_len : mi->seq[rid].len;
 	qe = qe + ext_len < qlen? qe + ext_len : qlen;
 	tseq = (uint8_t*)kmalloc(km, re - rs);
 	mm_idx_getseq(mi, rid, rs, re, tseq);
@@ -428,28 +550,30 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	int32_t i, l, bw, dropped = 0, extra_flag = 0, rs0, re0, qs0, qe0;
 	int32_t rs, re, qs, qe;
 	int32_t rs1, qs1, re1, qe1;
-	int8_t mat[25];
+	int8_t mat[25], *ob;
 
 	if (is_sr) assert(!(mi->flag & MM_I_HPC)); // HPC won't work with SR because with HPC we can't easily tell if there is a gap
 
 	r2->cnt = 0;
 	if (r->cnt == 0) return;
-	ksw_gen_simple_mat(5, mat, opt->a, opt->b);
+	ksw_gen_simple_mat(5, mat, opt->a, opt->b, opt->sc_ambi);
 	bw = (int)(opt->bw * 1.5 + 1.);
 
 	if (is_sr && !(mi->flag & MM_I_HPC)) {
-		mm_max_stretch(opt, r, a, &as1, &cnt1);
+		mm_max_stretch(r, a, &as1, &cnt1);
 		rs = (int32_t)a[as1].x + 1 - (int32_t)(a[as1].y>>32&0xff);
 		qs = (int32_t)a[as1].y + 1 - (int32_t)(a[as1].y>>32&0xff);
 		re = (int32_t)a[as1+cnt1-1].x + 1;
 		qe = (int32_t)a[as1+cnt1-1].y + 1;
 	} else {
-		if (is_splice) {
-			mm_fix_bad_ends_splice(km, opt, mi, r, mat, qlen, qseq0, a, &as1, &cnt1);
-		} else {
-			mm_fix_bad_ends(r, a, opt->bw, opt->min_chain_score * 2, &as1, &cnt1);
-		}
+		if (!(opt->flag & MM_F_NO_END_FLT)) {
+			if (is_splice)
+				mm_fix_bad_ends_splice(km, opt, mi, r, mat, qlen, qseq0, a, &as1, &cnt1);
+			else
+				mm_fix_bad_ends(r, a, opt->bw, opt->min_chain_score * 2, &as1, &cnt1);
+		} else as1 = r->as, cnt1 = r->cnt;
 		mm_filter_bad_seeds(km, as1, cnt1, a, 10, 40, opt->max_gap>>1, 10);
+		mm_filter_bad_seeds_alt(km, as1, cnt1, a, 30, opt->max_gap>>1);
 		mm_adjust_minier(mi, qseq0, &a[as1], &rs, &qs);
 		mm_adjust_minier(mi, qseq0, &a[as1 + cnt1 - 1], &re, &qe);
 	}
@@ -470,10 +594,12 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 		qs0 = 0, qe0 = qlen;
 		l = qs;
 		l += l * opt->a + opt->end_bonus > opt->q? (l * opt->a + opt->end_bonus - opt->q) / opt->e : 0;
+		l = l < opt->bw? l : opt->bw;
 		rs0 = rs - l > 0? rs - l : 0;
 		l = qlen - qe;
 		l += l * opt->a + opt->end_bonus > opt->q? (l * opt->a + opt->end_bonus - opt->q) / opt->e : 0;
-		re0 = re + l < mi->seq[rid].len? re + l : mi->seq[rid].len;
+		l = l < opt->bw? l : opt->bw;
+		re0 = re + l < (int32_t)mi->seq[rid].len? re + l : mi->seq[rid].len;
 	} else {
 		// compute rs0 and qs0
 		rs0 = (int32_t)a[r->as].x + 1 - (int32_t)(a[r->as].y>>32&0xff);
@@ -517,13 +643,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				}
 			}
 		}
-		if (qe < qlen && re < mi->seq[rid].len) {
+		if (qe < qlen && re < (int32_t)mi->seq[rid].len) {
 			l = qlen - qe < opt->max_gap? qlen - qe : opt->max_gap;
 			qe1 = qe1 < qe + l? qe1 : qe + l;
 			qe0 = qe0 > qe1? qe0 : qe1; // at least include qe0
 			l += l * opt->a > opt->q? (l * opt->a - opt->q) / opt->e : 0;
 			l = l < opt->max_gap? l : opt->max_gap;
-			l = l < mi->seq[rid].len - re? l : mi->seq[rid].len - re;
+			l = l < (int32_t)mi->seq[rid].len - re? l : mi->seq[rid].len - re;
 			re1 = re1 < re + l? re1 : re + l;
 			re0 = re0 > re1? re0 : re1;
 		} else re0 = re, qe0 = qe;
@@ -539,13 +665,14 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 
 	assert(re0 > rs0);
 	tseq = (uint8_t*)kmalloc(km, re0 - rs0);
+	ob = (int8_t*)kmalloc(km, re0 - rs0);
 
 	if (qs > 0 && rs > 0) { // left extension
 		qseq = &qseq0[rev][qs0];
-		mm_idx_getseq(mi, rid, rs0, rs, tseq);
+		mm_idx_getseq2(mi, rid, rs0, rs, tseq, ob);
 		mm_seq_rev(qs - qs0, qseq);
 		mm_seq_rev(rs - rs0, tseq);
-		mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, mat, bw, opt->end_bonus, r->split_inv? opt->zdrop_inv : opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez);
+		mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, ob, mat, bw, opt->end_bonus, r->split_inv? opt->zdrop_inv : opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez);
 		if (ez->n_cigar > 0) {
 			mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			r->p->dp_score += ez->max;
@@ -570,7 +697,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				bw1 = qe - qs > re - rs? qe - qs : re - rs;
 			// perform alignment
 			qseq = &qseq0[rev][qs];
-			mm_idx_getseq(mi, rid, rs, re, tseq);
+			mm_idx_getseq2(mi, rid, rs, re, tseq, ob);
 			if (is_sr) { // perform ungapped alignment
 				assert(qe - qs == re - rs);
 				ksw_reset_extz(ez);
@@ -580,11 +707,11 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				}
 				ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, 0, qe - qs);
 			} else { // perform normal gapped alignment
-				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
+				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, ob, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
 			}
 			// test Z-drop and inversion Z-drop
 			if ((zdrop_code = mm_test_zdrop(km, opt, qseq, tseq, ez->n_cigar, ez->cigar, mat)) != 0)
-				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, zdrop_code == 2? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate
+				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, ob, mat, bw1, -1, zdrop_code == 2? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate
 			// update CIGAR
 			if (ez->n_cigar > 0)
 				mm_append_cigar(r, ez->n_cigar, ez->cigar);
@@ -609,8 +736,8 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 
 	if (!dropped && qe < qe0 && re < re0) { // right extension
 		qseq = &qseq0[rev][qe];
-		mm_idx_getseq(mi, rid, re, re0, tseq);
-		mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
+		mm_idx_getseq2(mi, rid, re, re0, tseq, ob);
+		mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, ob, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
 		if (ez->n_cigar > 0) {
 			mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			r->p->dp_score += ez->max;
@@ -626,13 +753,16 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 
 	assert(re1 - rs1 <= re0 - rs0);
 	if (r->p) {
-		mm_idx_getseq(mi, rid, rs1, re1, tseq);
-		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e);
+		int left_aln = (opt->flag & MM_F_SPLICE) || mi->n_R == 0? 1 : 0;
+		mm_idx_getseq2(mi, rid, rs1, re1, tseq, ob);
+		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e, left_aln);
+		if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r, &qseq0[r->rev][qs1], tseq);
 		if (rev && r->p->trans_strand)
 			r->p->trans_strand ^= 3; // flip to the read strand
 	}
 
 	kfree(km, tseq);
+	kfree(km, ob);
 }
 
 static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez)
@@ -652,7 +782,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	if (ql < opt->min_chain_score || ql > opt->max_gap) return 0;
 	if (tl < opt->min_chain_score || tl > opt->max_gap) return 0;
 
-	ksw_gen_simple_mat(5, mat, opt->a, opt->b);
+	ksw_gen_simple_mat(5, mat, opt->a, opt->b, opt->sc_ambi);
 	tseq = (uint8_t*)kmalloc(km, tl);
 	mm_idx_getseq(mi, r1->rid, r1->re, r2->rs, tseq);
 	qseq = r1->rev? &qseq0[0][r2->qe] : &qseq0[1][qlen - r2->qs];
@@ -666,7 +796,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	mm_seq_rev(tl, tseq);
 	if (score < opt->min_dp_max) goto end_align1_inv;
 	q_off = ql - (q_off + 1), t_off = tl - (t_off + 1);
-	mm_align_pair(km, opt, ql - q_off, qseq + q_off, tl - t_off, tseq + t_off, mat, (int)(opt->bw * 1.5), -1, opt->zdrop, KSW_EZ_EXTZ_ONLY, ez);
+	mm_align_pair(km, opt, ql - q_off, qseq + q_off, tl - t_off, tseq + t_off, 0, mat, (int)(opt->bw * 1.5), -1, opt->zdrop, KSW_EZ_EXTZ_ONLY, ez);
 	if (ez->n_cigar == 0) goto end_align1_inv; // should never be here
 	mm_append_cigar(r_inv, ez->n_cigar, ez->cigar);
 	r_inv->p->dp_score = ez->max;
@@ -685,7 +815,8 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	}
 	r_inv->rs = r1->re + t_off;
 	r_inv->re = r_inv->rs + ez->max_t + 1;
-	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e);
+	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, 0);
+	if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r_inv, &qseq[q_off], &tseq[t_off]);
 	ret = 1;
 end_align1_inv:
 	kfree(km, tseq);
@@ -755,7 +886,7 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 	*n_regs_ = n_regs;
 	kfree(km, qseq0[0]);
 	kfree(km, ez.cigar);
-	mm_filter_regs(km, opt, qlen, n_regs_, regs);
-	mm_hit_sort_by_dp(km, n_regs_, regs);
+	mm_filter_regs(opt, qlen, n_regs_, regs);
+	mm_hit_sort(km, n_regs_, regs);
 	return regs;
 }

bseq.c

@@ -67,7 +67,7 @@ static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual, int with_
 	int i;
 	s->name = kstrdup(&ks->name);
 	s->seq = kstrdup(&ks->seq);
-	for (i = 0; i < ks->seq.l; ++i) // convert U to T
+	for (i = 0; i < (int)ks->seq.l; ++i) // convert U to T
 		if (s->seq[i] == 'u' || s->seq[i] == 'U')
 			--s->seq[i];
 	s->qual = with_qual && ks->qual.l? kstrdup(&ks->qual) : 0;

chain.c

@@ -44,7 +44,7 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 		int32_t qi = (int32_t)a[i].y, q_span = a[i].y>>32&0xff; // NB: only 8 bits of span is used!!!
 		int32_t max_f = q_span, n_skip = 0, min_d;
 		int32_t sidi = (a[i].y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
-		while (st < i && ri - a[st].x > max_dist_x) ++st;
+		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, log_dd;

cookbook.md

@@ -31,8 +31,8 @@ To acquire the data used in this cookbook and to install minimap2 and paftools,
 please follow the command lines below:
 ```sh
 # install minimap2 executables
-curl -L https://github.com/lh3/minimap2/releases/download/v2.10/minimap2-2.10_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.10_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
+curl -L https://github.com/lh3/minimap2/releases/download/v2.13/minimap2-2.13_x64-linux.tar.bz2 | tar jxf -
+cp minimap2-2.13_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
 export PATH="$PATH:"`pwd`                               # put the current directory on PATH
 # download example datasets
 curl -L https://github.com/lh3/minimap2/releases/download/v2.10/cookbook-data.tgz | tar zxf -

format.c

@@ -133,11 +133,11 @@ void mm_write_sam_hdr(const mm_idx_t *idx, const char *rg, const char *ver, int
 	free(str.s);
 }
 
-static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp, int no_iden)
+static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp, int no_iden, int write_tag)
 {
 	int i, q_off, t_off;
-	mm_sprintf_lite(s, "\tcs:Z:");
-	for (i = q_off = t_off = 0; i < r->p->n_cigar; ++i) {
+	if (write_tag) mm_sprintf_lite(s, "\tcs:Z:");
+	for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
 		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
 		assert(op >= 0 && op <= 3);
 		if (op == 0) { // match
@@ -181,13 +181,13 @@ static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 	assert(t_off == r->re - r->rs && q_off == r->qe - r->qs);
 }
 
-static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp)
+static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp, int write_tag)
 {
 	int i, q_off, t_off, l_MD = 0;
-	mm_sprintf_lite(s, "\tMD:Z:");
-	for (i = q_off = t_off = 0; i < r->p->n_cigar; ++i) {
+	if (write_tag) mm_sprintf_lite(s, "\tMD:Z:");
+	for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
 		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
-		assert(op >= 0 && op <= 2); // introns (aka reference skips) are not supported
+		assert(op >= 0 && op <= 3);
 		if (op == 0) { // match
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
@@ -204,13 +204,15 @@ static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 			mm_sprintf_lite(s, "%d^%s", l_MD, tmp);
 			l_MD = 0;
 			t_off += len;
+		} else if (op == 3) { // reference skip
+			t_off += len;
 		}
 	}
 	if (l_MD > 0) mm_sprintf_lite(s, "%d", l_MD);
 	assert(t_off == r->re - r->rs && q_off == r->qe - r->qs);
 }
 
-static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int no_iden, int is_MD)
+static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int no_iden, int is_MD, int write_tag)
 {
 	extern unsigned char seq_nt4_table[256];
 	int i;
@@ -230,11 +232,34 @@ static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_
 			qseq[r->qe - i - 1] = c >= 4? 4 : 3 - c;
 		}
 	}
-	if (is_MD) write_MD_core(s, tseq, qseq, r, tmp);
-	else write_cs_core(s, tseq, qseq, r, tmp, no_iden);
+	if (is_MD) write_MD_core(s, tseq, qseq, r, tmp, write_tag);
+	else write_cs_core(s, tseq, qseq, r, tmp, no_iden, write_tag);
 	kfree(km, qseq); kfree(km, tseq); kfree(km, tmp);
 }
 
+int mm_gen_cs_or_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int is_MD, int no_iden)
+{
+	mm_bseq1_t t;
+	kstring_t str;
+	str.s = *buf, str.l = 0, str.m = *max_len;
+	t.l_seq = strlen(seq);
+	t.seq = (char*)seq;
+	write_cs_or_MD(km, &str, mi, &t, r, no_iden, is_MD, 0);
+	*max_len = str.m;
+	*buf = str.s;
+	return str.l;
+}
+
+int mm_gen_cs(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int no_iden)
+{
+	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 0, no_iden);
+}
+
+int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq)
+{
+	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 1, 0);
+}
+
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 {
 	int type;
@@ -259,6 +284,10 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 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)
 {
 	s->l = 0;
+	if (r == 0) {
+		mm_sprintf_lite(s, "%s\t%d", t->name, t->l_seq);
+		return;
+	}
 	mm_sprintf_lite(s, "%s\t%d\t%d\t%d\t%c\t", t->name, t->l_seq, r->qs, r->qe, "+-"[r->rev]);
 	if (mi->seq[r->rid].name) mm_sprintf_lite(s, "%s", mi->seq[r->rid].name);
 	else mm_sprintf_lite(s, "%d", r->rid);
@@ -270,10 +299,10 @@ void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m
 		uint32_t k;
 		mm_sprintf_lite(s, "\tcg:Z:");
 		for (k = 0; k < r->p->n_cigar; ++k)
-			mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDN"[r->p->cigar[k]&0xf]);
+			mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDNSHP=XB"[r->p->cigar[k]&0xf]);
 	}
 	if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_MD)))
-		write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD);
+		write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1);
 	if ((opt_flag & MM_F_COPY_COMMENT) && t->comment)
 		mm_sprintf_lite(s, "\t%s", t->comment);
 }
@@ -321,7 +350,7 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
 			int clip_char = (sam_flag&0x800) && !(opt_flag&MM_F_SOFTCLIP)? 'H' : 'S';
 			if (clip_len[0]) mm_sprintf_lite(s, "%d%c", clip_len[0], clip_char);
 			for (k = 0; k < r->p->n_cigar; ++k)
-				mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDN"[r->p->cigar[k]&0xf]);
+				mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDNSHP=XB"[r->p->cigar[k]&0xf]);
 			if (clip_len[1]) mm_sprintf_lite(s, "%d%c", clip_len[1], clip_char);
 		}
 	}
@@ -390,9 +419,11 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 				cigar_in_tag = 1;
 		}
 		if (cigar_in_tag) {
-			if (flag & 0x100) mm_sprintf_lite(s, "0S"); // secondary alignment
-			else if (flag & 0x800) mm_sprintf_lite(s, "%dS", r->re - r->rs); // supplementary alignment
-			else mm_sprintf_lite(s, "%dS", t->l_seq);
+			int slen;
+			if ((flag & 0x900) == 0 || (opt_flag & MM_F_SOFTCLIP)) slen = t->l_seq;
+			else if (flag & 0x100) slen = 0;
+			else slen = r->qe - r->qs;
+			mm_sprintf_lite(s, "%dS%dN", slen, r->re - r->rs);
 		} else write_sam_cigar(s, flag, 0, t->l_seq, r, opt_flag);
 	}
 
@@ -472,7 +503,7 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 			}
 		}
 		if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_MD)))
-			write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD);
+			write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1);
 		if (cigar_in_tag)
 			write_sam_cigar(s, flag, 1, t->l_seq, r, opt_flag);
 	}

getopt.c

@@ -1,216 +0,0 @@
-#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

@@ -1,53 +0,0 @@
-/*
-  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

@@ -106,7 +106,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, int sub_diff) // and compute mm_reg1_t::subsc
+void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff, int hard_mask_level) // and compute mm_reg1_t::subsc
 {
 	int i, j, k, *w;
 	uint64_t *cov;
@@ -118,6 +118,7 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff
 	for (i = 1, k = 1; i < n; ++i) {
 		mm_reg1_t *ri = &r[i];
 		int si = ri->qs, ei = ri->qe, n_cov = 0, uncov_len = 0;
+		if (hard_mask_level) goto skip_uncov;
 		for (j = 0; j < k; ++j) { // traverse existing primary hits to find overlapping hits
 			mm_reg1_t *rp = &r[w[j]];
 			int sj = rp->qs, ej = rp->qe;
@@ -132,18 +133,19 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff
 			int j, x = si;
 			radix_sort_64(cov, cov + n_cov);
 			for (j = 0; j < n_cov; ++j) {
-				if (cov[j]>>32 > x) uncov_len += (cov[j]>>32) - x;
+				if ((int)(cov[j]>>32) > x) uncov_len += (cov[j]>>32) - x;
 				x = (int32_t)cov[j] > x? (int32_t)cov[j] : x;
 			}
 			if (ei > x) uncov_len += ei - x;
 		}
+skip_uncov:
 		for (j = 0; j < k; ++j) { // traverse existing primary hits again
 			mm_reg1_t *rp = &r[w[j]];
 			int sj = rp->qs, ej = rp->qe, min, max, ol;
 			if (ej <= si || sj >= ei) continue; // no overlap
 			min = ej - sj < ei - si? ej - sj : ei - si;
 			max = ej - sj > ei - si? ej - sj : ei - si;
-			ol = si < sj? (ei < sj? 0 : ei < ej? ei - sj : ej - sj) : (ej < si? 0 : ej < ei? ej - si : ei - si); // overlap length
+			ol = si < sj? (ei < sj? 0 : ei < ej? ei - sj : ej - sj) : (ej < si? 0 : ej < ei? ej - si : ei - si); // overlap length; TODO: this can be simplified
 			if ((float)ol / min - (float)uncov_len / max > mask_level) {
 				int cnt_sub = 0;
 				ri->parent = rp->parent;
@@ -164,9 +166,9 @@ set_parent_test:
 	kfree(km, w);
 }
 
-void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r)
+void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r)
 {
-	int32_t i, n_aux, n = *n_regs;
+	int32_t i, n_aux, n = *n_regs, has_cigar = 0, no_cigar = 0;
 	mm128_t *aux;
 	mm_reg1_t *t;
 
@@ -175,14 +177,20 @@ void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r)
 	t = (mm_reg1_t*)kmalloc(km, n * sizeof(mm_reg1_t));
 	for (i = n_aux = 0; i < n; ++i) {
 		if (r[i].inv || r[i].cnt > 0) { // squeeze out elements with cnt==0 (soft deleted)
-			assert(r[i].p);
-			aux[n_aux].x = (uint64_t)r[i].p->dp_max << 32 | r[i].hash;
+			if (r[i].p) {
+				aux[n_aux].x = (uint64_t)r[i].p->dp_max << 32 | r[i].hash;
+				has_cigar = 1;
+			} else {
+				aux[n_aux].x = (uint64_t)r[i].score << 32 | r[i].hash;
+				no_cigar = 1;
+			}
 			aux[n_aux++].y = i;
 		} else if (r[i].p) {
 			free(r[i].p);
 			r[i].p = 0;
 		}
 	}
+	assert(has_cigar + no_cigar == 1);
 	radix_sort_128x(aux, aux + n_aux);
 	for (i = n_aux - 1; i >= 0; --i)
 		t[n_aux - 1 - i] = r[aux[i].y];
@@ -246,7 +254,7 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_,
 	}
 }
 
-void mm_filter_regs(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs)
+void mm_filter_regs(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs)
 { // NB: after this call, mm_reg1_t::parent can be -1 if its parent filtered out
 	int i, k;
 	for (i = k = 0; i < *n_regs; ++i) {
@@ -304,7 +312,7 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
 	for (i = n_aux - 1; i >= 1; --i) {
 		mm_reg1_t *r0 = &regs[(int32_t)aux[i-1]], *r1 = &regs[(int32_t)aux[i]];
 		mm128_t *a0e, *a1s;
-		int max_gap, min_gap, sc_thres;
+		int max_gap, min_gap, sc_thres, min_flank_len;
 
 		// test
 		if (r0->as + r0->cnt != r1->as) continue; // not adjacent in a[]
@@ -313,13 +321,14 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
 		a1s = &a[r1->as];
 		if (a1s->x <= a0e->x || (int32_t)a1s->y <= (int32_t)a0e->y) continue; // keep colinearity
 		max_gap = min_gap = (int32_t)a1s->y - (int32_t)a0e->y;
-		max_gap = max_gap > a1s->x - a0e->x? max_gap : a1s->x - a0e->x;
-		min_gap = min_gap < a1s->x - a0e->x? min_gap : a1s->x - a0e->x;
+		max_gap = a0e->x + max_gap > a1s->x? max_gap : a1s->x - a0e->x;
+		min_gap = a0e->x + min_gap < a1s->x? min_gap : a1s->x - a0e->x;
 		if (max_gap > opt->max_join_long || min_gap > opt->max_join_short) continue;
 		sc_thres = (int)((float)opt->min_join_flank_sc / opt->max_join_long * max_gap + .499);
 		if (r0->score < sc_thres || r1->score < sc_thres) continue; // require good flanking chains
-		if (r0->re - r0->rs < max_gap>>1 || r0->qe - r0->qs < max_gap>>1) continue; // require enough flanking length
-		if (r1->re - r1->rs < max_gap>>1 || r1->qe - r1->qs < max_gap>>1) continue;
+		min_flank_len = (int)(max_gap * opt->min_join_flank_ratio);
+		if (r0->re - r0->rs < min_flank_len || r0->qe - r0->qs < min_flank_len) continue; // require enough flanking length
+		if (r1->re - r1->rs < min_flank_len || r1->qe - r1->qs < min_flank_len) continue;
 
 		// all conditions satisfied; join
 		a[r1->as].y |= MM_SEED_LONG_JOIN;
@@ -339,7 +348,7 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
 					r->parent = regs[r->parent].parent;
 			}
 		}
-		mm_filter_regs(km, opt, qlen, n_regs_, regs);
+		mm_filter_regs(opt, qlen, n_regs_, regs);
 		mm_sync_regs(km, *n_regs_, regs);
 	}
 }
@@ -411,23 +420,23 @@ void mm_seg_free(void *km, int n_segs, mm_seg_t *segs)
 static void mm_set_inv_mapq(void *km, int n_regs, mm_reg1_t *regs)
 {
 	int i, n_aux;
-	uint64_t *aux;
+	mm128_t *aux;
 	if (n_regs < 3) return;
 	for (i = 0; i < n_regs; ++i)
 		if (regs[i].inv) break;
 	if (i == n_regs) return; // no inversion hits
 
-	aux = (uint64_t*)kmalloc(km, n_regs * 8);
+	aux = (mm128_t*)kmalloc(km, n_regs * 16);
 	for (i = n_aux = 0; i < n_regs; ++i)
 		if (regs[i].parent == i || regs[i].parent < 0)
-			aux[n_aux++] = (uint64_t)regs[i].as << 32 | i;
-	radix_sort_64(aux, aux + n_aux);
+			aux[n_aux].y = i, aux[n_aux++].x = (uint64_t)regs[i].rid << 32 | regs[i].rs;
+	radix_sort_128x(aux, aux + n_aux);
 
 	for (i = 1; i < n_aux - 1; ++i) {
-		mm_reg1_t *inv = &regs[(int32_t)aux[i]];
+		mm_reg1_t *inv = &regs[aux[i].y];
 		if (inv->inv) {
-			mm_reg1_t *l = &regs[(int32_t)aux[i-1]];
-			mm_reg1_t *r = &regs[(int32_t)aux[i+1]];
+			mm_reg1_t *l = &regs[aux[i-1].y];
+			mm_reg1_t *r = &regs[aux[i+1].y];
 			inv->mapq = l->mapq < r->mapq? l->mapq : r->mapq;
 		}
 	}

index.c

@@ -45,20 +45,22 @@ mm_idx_t *mm_idx_init(int w, int k, int b, int flag)
 
 void mm_idx_destroy(mm_idx_t *mi)
 {
-	int i;
+	uint32_t i;
 	if (mi == 0) return;
 	if (mi->h) kh_destroy(str, (khash_t(str)*)mi->h);
-	for (i = 0; i < 1<<mi->b; ++i) {
-		free(mi->B[i].p);
-		free(mi->B[i].a.a);
-		kh_destroy(idx, (idxhash_t*)mi->B[i].h);
+	if (mi->B) {
+		for (i = 0; i < 1U<<mi->b; ++i) {
+			free(mi->B[i].p);
+			free(mi->B[i].a.a);
+			kh_destroy(idx, (idxhash_t*)mi->B[i].h);
+		}
 	}
 	if (!mi->km) {
 		for (i = 0; i < mi->n_seq; ++i)
 			free(mi->seq[i].name);
 		free(mi->seq);
 	} else km_destroy(mi->km);
-	free(mi->B); free(mi->S); free(mi);
+	free(mi->R); free(mi->B); free(mi->S); free(mi);
 }
 
 const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n)
@@ -82,14 +84,15 @@ const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n)
 
 void mm_idx_stat(const mm_idx_t *mi)
 {
-	int i, n = 0, n1 = 0;
+	int n = 0, n1 = 0;
+	uint32_t i;
 	uint64_t sum = 0, len = 0;
 	fprintf(stderr, "[M::%s] kmer size: %d; skip: %d; is_hpc: %d; #seq: %d\n", __func__, mi->k, mi->w, mi->flag&MM_I_HPC, mi->n_seq);
 	for (i = 0; i < mi->n_seq; ++i)
 		len += mi->seq[i].len;
-	for (i = 0; i < 1<<mi->b; ++i)
+	for (i = 0; i < 1U<<mi->b; ++i)
 		if (mi->B[i].h) n += kh_size((idxhash_t*)mi->B[i].h);
-	for (i = 0; i < 1<<mi->b; ++i) {
+	for (i = 0; i < 1U<<mi->b; ++i) {
 		idxhash_t *h = (idxhash_t*)mi->B[i].h;
 		khint_t k;
 		if (h == 0) continue;
@@ -131,7 +134,7 @@ int mm_idx_name2id(const mm_idx_t *mi, const char *name)
 	return k == kh_end(h)? -1 : kh_val(h, k);
 }
 
-int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq)
+int mm_idx_getseq2(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq, int8_t *b)
 {
 	uint64_t i, st1, en1;
 	if (rid >= mi->n_seq || st >= mi->seq[rid].len) return -1;
@@ -140,9 +143,26 @@ int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, ui
 	en1 = mi->seq[rid].offset + en;
 	for (i = st1; i < en1; ++i)
 		seq[i - st1] = mm_seq4_get(mi->S, i);
+	if (b) memset(b, 0, en - st);
+	if (b && mi->R) {
+		uint32_t i, z;
+		memset(b, 0, en - st);
+		z = mm_idx_bed_query(mi, (uint64_t)rid << 32 | st);
+		for (i = z < 0? 0 : z; i < mi->n_R; ++i) {
+			uint32_t j, rr, rs, re;
+			rr = mi->R[i].x >> 32, rs = (uint32_t)mi->R[i].x, re = mi->R[i].end;
+			if (rr > rid || rs >= en) break;
+			if (rr < rid) continue;
+			re = re < en? re : en;
+			for (j = st > rs? st : rs; j < re; ++j)
+				b[j - st] = mi->R[i].score;
+		}
+	}
 	return en - st;
 }
 
+int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq) { return mm_idx_getseq2(mi, rid, st, en, seq, 0); }
+
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f)
 {
 	int i;
@@ -172,7 +192,8 @@ int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f)
 
 static void worker_post(void *g, long i, int tid)
 {
-	int j, start_a, start_p, n, n_keys;
+	int n, n_keys;
+	size_t j, start_a, start_p;
 	idxhash_t *h;
 	mm_idx_t *mi = (mm_idx_t*)g;
 	mm_idx_bucket_t *b = &mi->B[i];
@@ -200,7 +221,7 @@ static void worker_post(void *g, long i, int tid)
 			int absent;
 			mm128_t *p = &b->a.a[j-1];
 			itr = kh_put(idx, h, p->x>>8>>mi->b<<1, &absent);
-			assert(absent && j - start_a == n);
+			assert(absent && j == start_a + n);
 			if (n == 1) {
 				kh_key(h, itr) |= 1;
 				kh_val(h, itr) = p->y;
@@ -216,7 +237,7 @@ static void worker_post(void *g, long i, int tid)
 		} else ++n;
 	}
 	b->h = h;
-	assert(b->n == start_p);
+	assert(b->n == (int32_t)start_p);
 
 	// deallocate and clear b->a
 	kfree(0, b->a.a);
@@ -336,7 +357,7 @@ mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int flag, int mini
 	pipeline_t pl;
 	if (fp == 0 || mm_bseq_eof(fp)) return 0;
 	memset(&pl, 0, sizeof(pipeline_t));
-	pl.mini_batch_size = mini_batch_size < batch_size? mini_batch_size : batch_size;
+	pl.mini_batch_size = (uint64_t)mini_batch_size < batch_size? mini_batch_size : batch_size;
 	pl.batch_size = batch_size;
 	pl.fp = fp;
 	pl.mi = mm_idx_init(w, k, b, flag);
@@ -368,7 +389,9 @@ mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const cha
 	uint64_t sum_len = 0;
 	mm128_v a = {0,0,0};
 	mm_idx_t *mi;
+	khash_t(str) *h;
 	int i, flag = 0;
+
 	if (n <= 0) return 0;
 	for (i = 0; i < n; ++i) // get the total length
 		sum_len += strlen(seq[i]);
@@ -379,13 +402,17 @@ mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const cha
 	mi->n_seq = n;
 	mi->seq = (mm_idx_seq_t*)kcalloc(mi->km, n, sizeof(mm_idx_seq_t)); // ->seq is allocated from km
 	mi->S = (uint32_t*)calloc((sum_len + 7) / 8, 4);
+	mi->h = h = kh_init(str);
 	for (i = 0, sum_len = 0; i < n; ++i) {
 		const char *s = seq[i];
 		mm_idx_seq_t *p = &mi->seq[i];
 		uint32_t j;
 		if (name && name[i]) {
+			int absent;
 			p->name = (char*)kmalloc(mi->km, strlen(name[i]) + 1);
 			strcpy(p->name, name[i]);
+			kh_put(str, h, p->name, &absent);
+			assert(absent);
 		}
 		p->offset = sum_len;
 		p->len = strlen(s);
@@ -413,17 +440,20 @@ mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const cha
 void mm_idx_dump(FILE *fp, const mm_idx_t *mi)
 {
 	uint64_t sum_len = 0;
-	uint32_t x[5];
-	int i;
+	uint32_t x[5], i;
 
 	x[0] = mi->w, x[1] = mi->k, x[2] = mi->b, x[3] = mi->n_seq, x[4] = mi->flag;
 	fwrite(MM_IDX_MAGIC, 1, 4, fp);
 	fwrite(x, 4, 5, fp);
 	for (i = 0; i < mi->n_seq; ++i) {
-		uint8_t l;
-		l = strlen(mi->seq[i].name);
-		fwrite(&l, 1, 1, fp);
-		fwrite(mi->seq[i].name, 1, l, fp);
+		if (mi->seq[i].name) {
+			uint8_t l = strlen(mi->seq[i].name);
+			fwrite(&l, 1, 1, fp);
+			fwrite(mi->seq[i].name, 1, l, fp);
+		} else {
+			uint8_t l = 0;
+			fwrite(&l, 1, 1, fp);
+		}
 		fwrite(&mi->seq[i].len, 4, 1, fp);
 		sum_len += mi->seq[i].len;
 	}
@@ -450,9 +480,8 @@ void mm_idx_dump(FILE *fp, const mm_idx_t *mi)
 
 mm_idx_t *mm_idx_load(FILE *fp)
 {
-	int i;
 	char magic[4];
-	uint32_t x[5];
+	uint32_t x[5], i;
 	uint64_t sum_len = 0;
 	mm_idx_t *mi;
 
@@ -466,9 +495,11 @@ mm_idx_t *mm_idx_load(FILE *fp)
 		uint8_t l;
 		mm_idx_seq_t *s = &mi->seq[i];
 		fread(&l, 1, 1, fp);
-		s->name = (char*)kmalloc(mi->km, l + 1);
-		fread(s->name, 1, l, fp);
-		s->name[l] = 0;
+		if (l) {
+			s->name = (char*)kmalloc(mi->km, l + 1);
+			fread(s->name, 1, l, fp);
+			s->name[l] = 0;
+		}
 		fread(&s->len, 4, 1, fp);
 		s->offset = sum_len;
 		sum_len += s->len;
@@ -504,14 +535,19 @@ mm_idx_t *mm_idx_load(FILE *fp)
 int64_t mm_idx_is_idx(const char *fn)
 {
 	int fd, is_idx = 0;
-	off_t ret, off_end;
+	int64_t ret, off_end;
 	char magic[4];
 
 	if (strcmp(fn, "-") == 0) return 0; // read from pipe; not an index
 	fd = open(fn, O_RDONLY);
 	if (fd < 0) return -1; // error
+#ifdef WIN32
+	if ((off_end = _lseeki64(fd, 0, SEEK_END)) >= 4) {
+		_lseeki64(fd, 0, SEEK_SET);
+#else
 	if ((off_end = lseek(fd, 0, SEEK_END)) >= 4) {
 		lseek(fd, 0, SEEK_SET);
+#endif // WIN32
 		ret = read(fd, magic, 4);
 		if (ret == 4 && strncmp(magic, MM_IDX_MAGIC, 4) == 0)
 			is_idx = 1;
@@ -557,7 +593,7 @@ mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads)
 		mi = mm_idx_gen(r->fp.seq, r->opt.w, r->opt.k, r->opt.bucket_bits, r->opt.flag, r->opt.mini_batch_size, n_threads, r->opt.batch_size);
 	if (mi) {
 		if (r->fp_out) mm_idx_dump(r->fp_out, mi);
-		++r->n_parts;
+		mi->index = r->n_parts++;
 	}
 	return mi;
 }
@@ -566,3 +602,89 @@ int mm_idx_reader_eof(const mm_idx_reader_t *r) // TODO: in extremely rare cases
 {
 	return r->is_idx? (feof(r->fp.idx) || ftell(r->fp.idx) == r->idx_size) : mm_bseq_eof(r->fp.seq);
 }
+
+#include <ctype.h>
+#include <zlib.h>
+#include "ksort.h"
+#include "kseq.h"
+KSTREAM_DECLARE(gzFile, gzread)
+
+#define sort_key_bed(a) ((a).x)
+KRADIX_SORT_INIT(bed, mm_idx_bed_t, sort_key_bed, 8)
+
+mm_idx_bed_t *mm_idx_bed_read_list(const mm_idx_t *mi, const char *fn, uint32_t *n_)
+{
+	gzFile fp;
+	kstream_t *ks;
+	kstring_t str = {0,0,0};
+	uint32_t n = 0, m = 0;
+	mm_idx_bed_t *r = 0;
+
+	fp = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
+	if (fp == 0) return 0;
+	ks = ks_init(fp);
+	while (ks_getuntil(ks, KS_SEP_LINE, &str, 0) >= 0) {
+		mm_idx_bed_t t;
+		char *p, *q;
+		int i, id = -1, st = -1, en = -1, sc = -1;
+		for (p = q = str.s, i = 0;; ++p) {
+			if (*p == 0 || isspace(*p)) {
+				int32_t c = *p;
+				*p = 0;
+				if (i == 0) { // chr
+					id = mm_idx_name2id(mi, q);
+					if (id < 0) break; // unknown name; TODO: throw a warning
+				} else if (i == 1) { // start
+					st = atoi(q);
+					if (st < 0) break;
+				} else if (i == 2) { // end
+					en = atoi(q);
+					if (en < 0) break;
+				} else if (i == 3) { // name; do nothing
+				} else if (i == 4) { // BED score
+					sc = atoi(q);
+					assert(sc >= 0 && sc <= 127);
+				} else break;
+				if (c == 0) break;
+				++i, q = p + 1;
+			}
+		}
+		if (en < 0) en = st + 1;
+		if (st < 0 || st >= en) continue;
+		if (m == n) EXPAND(r, m);
+		t.x = (uint64_t)id << 32 | st, t.end = en, t.score = sc >= 0? sc : 0, t.idx = -1;
+		r[n++] = t;
+	}
+	ks_destroy(ks);
+	gzclose(fp);
+	*n_ = n;
+	return r;
+}
+
+int mm_idx_bed_attach(mm_idx_t *mi, uint32_t n, mm_idx_bed_t *r) // TODO: check errors
+{
+	radix_sort_bed(r, r + n);
+	mi->R = r, mi->n_R = n;
+	return 0;
+}
+
+int mm_idx_bed_read(mm_idx_t *mi, const char *fn)
+{
+	mm_idx_bed_t *r;
+	uint32_t n;
+	if (mi->h == 0) mm_idx_index_name(mi);
+	r = mm_idx_bed_read_list(mi, fn, &n);
+	return mm_idx_bed_attach(mi, n, r);
+}
+
+int mm_idx_bed_query(const mm_idx_t *mi, uint64_t x)
+{
+	int32_t left = -1, right = mi->n_R;
+	while (right - left > 1) {
+		int32_t mid = left + ((right - left) >> 1);
+		if (mi->R[mid].x > x) right = mid;
+		else if (mi->R[mid].x < x) left = mid;
+		else return mid;
+	}
+	return left;
+}

ketopt.h

@@ -0,0 +1,116 @@
+#ifndef KETOPT_H
+#define KETOPT_H
+
+#include <string.h> /* for strchr() and strncmp() */
+
+#define ko_no_argument       0
+#define ko_required_argument 1
+#define ko_optional_argument 2
+
+typedef struct {
+	int ind;   /* equivalent to optind */
+	int opt;   /* equivalent to optopt */
+	char *arg; /* equivalent to optarg */
+	int longidx; /* index of a long option; or -1 if short */
+	/* private variables not intended for external uses */
+	int i, pos, n_args;
+} ketopt_t;
+
+typedef struct {
+	char *name;
+	int has_arg;
+	int val;
+} ko_longopt_t;
+
+static ketopt_t KETOPT_INIT = { 1, 0, 0, -1, 1, 0, 0 };
+
+static void ketopt_permute(char *argv[], int j, int n) /* move argv[j] over n elements to the left */
+{
+	int k;
+	char *p = argv[j];
+	for (k = 0; k < n; ++k)
+		argv[j - k] = argv[j - k - 1];
+	argv[j - k] = p;
+}
+
+/**
+ * Parse command-line options and arguments
+ *
+ * This fuction has a similar interface to GNU's getopt_long(). Each call
+ * parses one option and returns the option name.  s->arg points to the option
+ * argument if present. The function returns -1 when all command-line arguments
+ * are parsed. In this case, s->ind is the index of the first non-option
+ * argument.
+ *
+ * @param s         status; shall be initialized to KETOPT_INIT on the first call
+ * @param argc      length of argv[]
+ * @param argv      list of command-line arguments; argv[0] is ignored
+ * @param permute   non-zero to move options ahead of non-option arguments
+ * @param ostr      option string
+ * @param longopts  long options
+ *
+ * @return ASCII for a short option; ko_longopt_t::val for a long option; -1 if
+ *         argv[] is fully processed; '?' for an unknown option or an ambiguous
+ *         long option; ':' if an option argument is missing
+ */
+static int ketopt(ketopt_t *s, int argc, char *argv[], int permute, const char *ostr, const ko_longopt_t *longopts)
+{
+	int opt = -1, i0, j;
+	if (permute) {
+		while (s->i < argc && (argv[s->i][0] != '-' || argv[s->i][1] == '\0'))
+			++s->i, ++s->n_args;
+	}
+	s->arg = 0, s->longidx = -1, i0 = s->i;
+	if (s->i >= argc || argv[s->i][0] != '-' || argv[s->i][1] == '\0') {
+		s->ind = s->i - s->n_args;
+		return -1;
+	}
+	if (argv[s->i][0] == '-' && argv[s->i][1] == '-') { /* "--" or a long option */
+		if (argv[s->i][2] == '\0') { /* a bare "--" */
+			ketopt_permute(argv, s->i, s->n_args);
+			++s->i, s->ind = s->i - s->n_args;
+			return -1;
+		}
+		s->opt = 0, opt = '?', s->pos = -1;
+		if (longopts) { /* parse long options */
+			int k, n_matches = 0;
+			const ko_longopt_t *o = 0;
+			for (j = 2; argv[s->i][j] != '\0' && argv[s->i][j] != '='; ++j) {} /* find the end of the option name */
+			for (k = 0; longopts[k].name != 0; ++k)
+				if (strncmp(&argv[s->i][2], longopts[k].name, j - 2) == 0)
+					++n_matches, o = &longopts[k];
+			if (n_matches == 1) {
+				s->opt = opt = o->val, s->longidx = o - longopts;
+				if (argv[s->i][j] == '=') s->arg = &argv[s->i][j + 1];
+				if (o->has_arg == 1 && argv[s->i][j] == '\0') {
+					if (s->i < argc - 1) s->arg = argv[++s->i];
+					else opt = ':'; /* missing option argument */
+				}
+			}
+		}
+	} else { /* a short option */
+		char *p;
+		if (s->pos == 0) s->pos = 1;
+		opt = s->opt = argv[s->i][s->pos++];
+		p = strchr(ostr, opt);
+		if (p == 0) {
+			opt = '?'; /* unknown option */
+		} else if (p[1] == ':') {
+			if (argv[s->i][s->pos] == 0) {
+				if (s->i < argc - 1) s->arg = argv[++s->i];
+				else opt = ':'; /* missing option argument */
+			} else s->arg = &argv[s->i][s->pos];
+			s->pos = -1;
+		}
+	}
+	if (s->pos < 0 || argv[s->i][s->pos] == 0) {
+		++s->i, s->pos = 0;
+		if (s->n_args > 0) /* permute */
+			for (j = i0; j < s->i; ++j)
+				ketopt_permute(argv, j, s->n_args);
+	}
+	s->ind = s->i - s->n_args;
+	return opt;
+}
+
+#endif

kseq.h

@@ -37,6 +37,14 @@
 #define KS_SEP_LINE  2 // line separator: "\n" (Unix) or "\r\n" (Windows)
 #define KS_SEP_MAX   2
 
+#ifndef klib_unused
+#if (defined __clang__ && __clang_major__ >= 3) || (defined __GNUC__ && __GNUC__ >= 3)
+#define klib_unused __attribute__ ((__unused__))
+#else
+#define klib_unused
+#endif
+#endif /* klib_unused */
+
 #define __KS_TYPE(type_t) \
 	typedef struct __kstream_t { \
 		int begin, end; \
@@ -64,7 +72,7 @@
 	}
 
 #define __KS_INLINED(__read) \
-	static inline int ks_getc(kstream_t *ks) \
+	static inline klib_unused int ks_getc(kstream_t *ks) \
 	{ \
 		if (ks->is_eof && ks->begin >= ks->end) return -1; \
 		if (ks->begin >= ks->end) { \

ksw2.h

@@ -58,7 +58,7 @@ void ksw_extd(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t
 			  int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int flag, ksw_extz_t *ez);
 
 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 gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
+				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez);
 
 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 gapo, int8_t gape, int8_t gapo2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
@@ -127,11 +127,11 @@ static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int min_intro
 			r = i + j;
 			if (i < off[r]) force_state = 2;
 			if (off_end && i > off_end[r]) force_state = 1;
-			tmp = force_state < 0? p[r * n_col + i - off[r]] : 0;
+			tmp = force_state < 0? p[(size_t)r * n_col + i - off[r]] : 0;
 		} else {
 			if (j < off[i]) force_state = 2;
 			if (off_end && j > off_end[i]) force_state = 1;
-			tmp = force_state < 0? p[i * n_col + j - off[i]] : 0;
+			tmp = force_state < 0? p[(size_t)i * n_col + j - off[i]] : 0;
 		}
 		if (state == 0) state = tmp & 7; // if requesting the H state, find state one maximizes it.
 		else if (!(tmp >> (state + 2) & 1)) state = 0; // if requesting other states, _state_ stays the same if it is a continuation; otherwise, set to H

ksw2_dispatch.c

@@ -64,18 +64,18 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 }
 
 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 end_bonus, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, 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 end_bonus, int flag, ksw_extz_t *ez);
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, 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 end_bonus, int flag, ksw_extz_t *ez);
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, 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, end_bonus, flag, ez);
+		ksw_extd2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, qd, ez);
 	else if (simd & SIMD_SSE2)
-		ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, ez);
+		ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, qd, ez);
 	else abort();
 }
 

ksw2_extd2_sse.c

@@ -17,17 +17,18 @@
 #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 end_bonus, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, 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 end_bonus, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, 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 end_bonus, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez)
 #endif // ~KSW_CPU_DISPATCH
 {
 #define __dp_code_block1 \
+	dt = _mm_load_si128(&dv[t]); \
 	z = _mm_load_si128(&s[t]); \
 	xt1 = _mm_load_si128(&x[t]);                     /* xt1 <- x[r-1][t..t+15] */ \
 	tmp = _mm_srli_si128(xt1, 15);                   /* tmp <- x[r-1][t+15] */ \
@@ -50,10 +51,10 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 #define __dp_code_block2 \
 	_mm_store_si128(&u[t], _mm_sub_epi8(z, vt1));    /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \
 	_mm_store_si128(&v[t], _mm_sub_epi8(z, ut));     /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \
-	tmp = _mm_sub_epi8(z, q_); \
+	tmp = _mm_sub_epi8(z, _mm_add_epi8(dt, q_)); \
 	a = _mm_sub_epi8(a, tmp); \
 	b = _mm_sub_epi8(b, tmp); \
-	tmp = _mm_sub_epi8(z, q2_); \
+	tmp = _mm_sub_epi8(z, _mm_add_epi8(dt, q2_)); \
 	a2= _mm_sub_epi8(a2, tmp); \
 	b2= _mm_sub_epi8(b2, tmp);
 
@@ -61,8 +62,8 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX);
 	int32_t *H = 0, H0 = 0, last_H0_t = 0;
 	uint8_t *qr, *sf, *mem, *mem2 = 0;
-	__m128i q_, q2_, qe_, qe2_, zero_, sc_mch_, sc_mis_, m1_, sc_N_;
-	__m128i *u, *v, *x, *y, *x2, *y2, *s, *p = 0;
+	__m128i q_, q2_, qe_, qe2_, e_, e2_, zero_, sc_mch_, sc_mis_, m1_, sc_N_;
+	__m128i *u, *v, *x, *y, *x2, *y2, *s, *p = 0, *dv;
 
 	ksw_reset_extz(ez);
 	if (m <= 1 || qlen <= 0 || tlen <= 0) return;
@@ -72,11 +73,13 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	zero_   = _mm_set1_epi8(0);
 	q_      = _mm_set1_epi8(q);
 	q2_     = _mm_set1_epi8(q2);
+	e_      = _mm_set1_epi8(e);
+	e2_     = _mm_set1_epi8(e2);
 	qe_     = _mm_set1_epi8(q + e);
 	qe2_    = _mm_set1_epi8(q2 + e2);
 	sc_mch_ = _mm_set1_epi8(mat[0]);
 	sc_mis_ = _mm_set1_epi8(mat[1]);
-	sc_N_   = _mm_set1_epi8(-e2);
+	sc_N_   = mat[m*m-1] == 0? _mm_set1_epi8(-e2) : _mm_set1_epi8(mat[m*m-1]);
 	m1_     = _mm_set1_epi8(m - 1); // wildcard
 
 	if (w < 0) w = tlen > qlen? tlen : qlen;
@@ -96,22 +99,29 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		++long_thres;
 	long_diff = long_thres * (e - e2) - (q2 - q) - e2;
 
-	mem = (uint8_t*)kcalloc(km, tlen_ * 8 + qlen_ + 1, 16);
+	mem = (uint8_t*)kcalloc(km, tlen_ * 9 + qlen_ + 1, 16);
 	u = (__m128i*)(((size_t)mem + 15) >> 4 << 4); // 16-byte aligned
-	v = u + tlen_, x = v + tlen_, y = x + tlen_, x2 = y + tlen_, y2 = x2 + tlen_;
-	s = y2 + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16;
+	v = u + tlen_, x = v + tlen_, y = x + tlen_, x2 = y + tlen_, y2 = x2 + tlen_, dv = y2 + tlen_;
+	s = dv + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16;
 	memset(u,  -q  - e,  tlen_ * 16);
 	memset(v,  -q  - e,  tlen_ * 16);
 	memset(x,  -q  - e,  tlen_ * 16);
 	memset(y,  -q  - e,  tlen_ * 16);
 	memset(x2, -q2 - e2, tlen_ * 16);
 	memset(y2, -q2 - e2, tlen_ * 16);
+	if (qd) {
+		int8_t *tmp = (int8_t*)dv;
+		for (t = 0; t < tlen; ++t) tmp[t] = -qd[t];
+		fprintf(stderr, "%d\t%d\t%x\n", tlen, qlen, flag&KSW_EZ_RIGHT);
+		for (t = 0; t < tlen; ++t) fputc("ACGTN"[target[t]], stderr); fputc('\n', stderr);
+		for (t = 0; t < tlen; ++t) fputc('0' + qd[t], stderr); fputc('\n', stderr);
+	}
 	if (!approx_max) {
 		H = (int32_t*)kmalloc(km, tlen_ * 16 * 4);
 		for (t = 0; t < tlen_ * 16; ++t) H[t] = KSW_NEG_INF;
 	}
 	if (with_cigar) {
-		mem2 = (uint8_t*)kmalloc(km, ((qlen + tlen - 1) * n_col_ + 1) * 16);
+		mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 16);
 		p = (__m128i*)(((size_t)mem2 + 15) >> 4 << 4);
 		off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2);
 		off_end = off + qlen + tlen - 1;
@@ -182,7 +192,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		assert(en_ - st_ + 1 <= n_col_);
 		if (!with_cigar) { // score only
 			for (t = st_; t <= en_; ++t) {
-				__m128i z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp;
+				__m128i z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp, dt;
 				__dp_code_block1;
 #ifdef __SSE4_1__
 				z = _mm_max_epi8(z, a);
@@ -191,10 +201,10 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				z = _mm_max_epi8(z, b2);
 				z = _mm_min_epi8(z, sc_mch_);
 				__dp_code_block2; // save u[] and v[]; update a, b, a2 and b2
-				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_max_epi8(a,  zero_), qe_));
-				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_max_epi8(b,  zero_), qe_));
-				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_max_epi8(a2, zero_), qe2_));
-				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_max_epi8(b2, zero_), qe2_));
+				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_max_epi8(a,  zero_), _mm_add_epi8(dt, qe_)));
+				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_max_epi8(b,  zero_), _mm_add_epi8(dt, qe_)));
+				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_max_epi8(a2, zero_), _mm_add_epi8(dt, qe2_)));
+				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_max_epi8(b2, zero_), _mm_add_epi8(dt, qe2_)));
 #else
 				tmp = _mm_cmpgt_epi8(a,  z);
 				z = _mm_or_si128(_mm_andnot_si128(tmp, z), _mm_and_si128(tmp, a));
@@ -208,20 +218,20 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				z = _mm_or_si128(_mm_and_si128(tmp, sc_mch_), _mm_andnot_si128(tmp, z));
 				__dp_code_block2;
 				tmp = _mm_cmpgt_epi8(a, zero_);
-				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_and_si128(tmp, a),  qe_));
+				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_and_si128(tmp, a),  _mm_add_epi8(dt, qe_)));
 				tmp = _mm_cmpgt_epi8(b, zero_);
-				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_and_si128(tmp, b),  qe_));
+				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_and_si128(tmp, b),  _mm_add_epi8(dt, qe_)));
 				tmp = _mm_cmpgt_epi8(a2, zero_);
-				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), qe2_));
+				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), _mm_add_epi8(dt, qe2_)));
 				tmp = _mm_cmpgt_epi8(b2, zero_);
-				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), qe2_));
+				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), _mm_add_epi8(dt, qe2_)));
 #endif
 			}
 		} else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment
-			__m128i *pr = p + r * n_col_ - st_;
+			__m128i *pr = p + (size_t)r * n_col_ - st_;
 			off[r] = st, off_end[r] = en;
 			for (t = st_; t <= en_; ++t) {
-				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp;
+				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp, dt;
 				__dp_code_block1;
 #ifdef __SSE4_1__
 				d = _mm_and_si128(_mm_cmpgt_epi8(a, z), _mm_set1_epi8(1));       // d = a  > z? 1 : 0
@@ -251,24 +261,24 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 #endif
 				__dp_code_block2;
 				tmp = _mm_cmpgt_epi8(a, zero_);
-				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_and_si128(tmp, a),  qe_));
+				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_and_si128(tmp, a),  _mm_add_epi8(dt, qe_)));
 				d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x08))); // d = a > 0? 1<<3 : 0
 				tmp = _mm_cmpgt_epi8(b, zero_);
-				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_and_si128(tmp, b),  qe_));
+				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_and_si128(tmp, b),  _mm_add_epi8(dt, qe_)));
 				d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x10))); // d = b > 0? 1<<4 : 0
 				tmp = _mm_cmpgt_epi8(a2, zero_);
-				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), qe2_));
+				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), _mm_add_epi8(dt, qe2_)));
 				d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x20))); // d = a > 0? 1<<5 : 0
 				tmp = _mm_cmpgt_epi8(b2, zero_);
-				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), qe2_));
+				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), _mm_add_epi8(dt, qe2_)));
 				d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x40))); // d = b > 0? 1<<6 : 0
 				_mm_store_si128(&pr[t], d);
 			}
 		} else { // gap right-alignment
-			__m128i *pr = p + r * n_col_ - st_;
+			__m128i *pr = p + (size_t)r * n_col_ - st_;
 			off[r] = st, off_end[r] = en;
 			for (t = st_; t <= en_; ++t) {
-				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp;
+				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp, dt;
 				__dp_code_block1;
 #ifdef __SSE4_1__
 				d = _mm_andnot_si128(_mm_cmpgt_epi8(z, a), _mm_set1_epi8(1));    // d = z > a?  0 : 1
@@ -298,16 +308,16 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 #endif
 				__dp_code_block2;
 				tmp = _mm_cmpgt_epi8(zero_, a);
-				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_andnot_si128(tmp, a),  qe_));
+				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_andnot_si128(tmp, a),  _mm_add_epi8(dt, qe_)));
 				d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x08))); // d = a > 0? 1<<3 : 0
 				tmp = _mm_cmpgt_epi8(zero_, b);
-				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_andnot_si128(tmp, b),  qe_));
+				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_andnot_si128(tmp, b),  _mm_add_epi8(dt, qe_)));
 				d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x10))); // d = b > 0? 1<<4 : 0
 				tmp = _mm_cmpgt_epi8(zero_, a2);
-				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, a2), qe2_));
+				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, a2), _mm_add_epi8(dt, qe2_)));
 				d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x20))); // d = a > 0? 1<<5 : 0
 				tmp = _mm_cmpgt_epi8(zero_, b2);
-				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, b2), qe2_));
+				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, b2), _mm_add_epi8(dt, qe2_)));
 				d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x40))); // d = b > 0? 1<<6 : 0
 				_mm_store_si128(&pr[t], d);
 			}
@@ -376,13 +386,14 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		last_st = st, last_en = en;
 		//for (t = st0; t <= en0; ++t) printf("(%d,%d)\t(%d,%d,%d,%d)\t%d\n", r, t, ((int8_t*)u)[t], ((int8_t*)v)[t], ((int8_t*)x)[t], ((int8_t*)y)[t], H[t]); // for debugging
 	}
+	fprintf(stderr, "score: %d\n", ez->score);
 	kfree(km, mem);
 	if (!approx_max) kfree(km, H);
 	if (with_cigar) { // backtrack
 		int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR);
 		if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY)) {
 			ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
-		} else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > ez->max) {
+		} else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > (int)ez->max) {
 			ez->reach_end = 1;
 			ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, ez->mqe_t, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
 		} else if (ez->max_t >= 0 && ez->max_q >= 0) {

ksw2_exts2_sse.c

@@ -71,7 +71,7 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	qe_     = _mm_set1_epi8(q + e);
 	sc_mch_ = _mm_set1_epi8(mat[0]);
 	sc_mis_ = _mm_set1_epi8(mat[1]);
-	sc_N_   = _mm_set1_epi8(-e);
+	sc_N_   = mat[m*m-1] == 0? _mm_set1_epi8(-e) : _mm_set1_epi8(mat[m*m-1]);
 	m1_     = _mm_set1_epi8(m - 1); // wildcard
 
 	tlen_ = (tlen + 15) / 16;
@@ -100,7 +100,7 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		for (t = 0; t < tlen_ * 16; ++t) H[t] = KSW_NEG_INF;
 	}
 	if (with_cigar) {
-		mem2 = (uint8_t*)kmalloc(km, ((qlen + tlen - 1) * n_col_ + 1) * 16);
+		mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 16);
 		p = (__m128i*)(((size_t)mem2 + 15) >> 4 << 4);
 		off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2);
 		off_end = off + qlen + tlen - 1;

ksw2_extz2_sse.c

@@ -65,7 +65,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	flag16_ = _mm_set1_epi8(0x10);
 	sc_mch_ = _mm_set1_epi8(mat[0]);
 	sc_mis_ = _mm_set1_epi8(mat[1]);
-	sc_N_   = _mm_set1_epi8(-e);
+	sc_N_   = mat[m*m-1] == 0? _mm_set1_epi8(-e) : _mm_set1_epi8(mat[m*m-1]);
 	m1_     = _mm_set1_epi8(m - 1); // wildcard
 	max_sc_ = _mm_set1_epi8(mat[0] + (q + e) * 2);
 
@@ -89,7 +89,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		for (t = 0; t < tlen_ * 16; ++t) H[t] = KSW_NEG_INF;
 	}
 	if (with_cigar) {
-		mem2 = (uint8_t*)kmalloc(km, ((qlen + tlen - 1) * n_col_ + 1) * 16);
+		mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 16);
 		p = (__m128i*)(((size_t)mem2 + 15) >> 4 << 4);
 		off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2);
 		off_end = off + qlen + tlen - 1;
@@ -169,7 +169,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 #endif
 			}
 		} else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment
-			__m128i *pr = p + r * n_col_ - st_;
+			__m128i *pr = p + (size_t)r * n_col_ - st_;
 			off[r] = st, off_end[r] = en;
 			for (t = st_; t <= en_; ++t) {
 				__m128i d, z, a, b, xt1, vt1, ut, tmp;
@@ -195,7 +195,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				_mm_store_si128(&pr[t], d);
 			}
 		} else { // gap right-alignment
-			__m128i *pr = p + r * n_col_ - st_;
+			__m128i *pr = p + (size_t)r * n_col_ - st_;
 			off[r] = st, off_end[r] = en;
 			for (t = st_; t <= en_; ++t) {
 				__m128i d, z, a, b, xt1, vt1, ut, tmp;
@@ -293,7 +293,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		int rev_cigar = !!(flag & KSW_EZ_REV_CIGAR);
 		if (!ez->zdropped && !(flag&KSW_EZ_EXTZ_ONLY)) {
 			ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, tlen-1, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
-		} else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > ez->max) {
+		} else if (!ez->zdropped && (flag&KSW_EZ_EXTZ_ONLY) && ez->mqe + end_bonus > (int)ez->max) {
 			ez->reach_end = 1;
 			ksw_backtrack(km, 1, rev_cigar, 0, (uint8_t*)p, off, off_end, n_col_*16, ez->mqe_t, qlen-1, &ez->m_cigar, &ez->n_cigar, &ez->cigar);
 		} else if (ez->max_t >= 0 && ez->max_q >= 0) {

main.c

@@ -4,13 +4,9 @@
 #include "bseq.h"
 #include "minimap.h"
 #include "mmpriv.h"
-#ifdef HAVE_GETOPT
-#include <getopt.h>
-#else
-#include "getopt.h"
-#endif
+#include "ketopt.h"
 
-#define MM_VERSION "2.10-r761"
+#define MM_VERSION "2.13-r852-dirty"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -26,53 +22,61 @@ void liftrlimit()
 void liftrlimit() {}
 #endif
 
-static struct option long_options[] = {
-	{ "bucket-bits",    required_argument, 0, 0 },
-	{ "mb-size",        required_argument, 0, 'K' },
-	{ "seed",           required_argument, 0, 0 },
-	{ "no-kalloc",      no_argument,       0, 0 },
-	{ "print-qname",    no_argument,       0, 0 },
-	{ "no-self",        no_argument,       0, 'D' },
-	{ "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, 'C' },
-	{ "no-long-join",   no_argument,       0, 0 },
-	{ "sr",             no_argument,       0, 0 },
-	{ "frag",           required_argument, 0, 0 },
-	{ "secondary",      required_argument, 0, 0 },
-	{ "cs",             optional_argument, 0, 0 },
-	{ "end-bonus",      required_argument, 0, 0 },
-	{ "no-pairing",     no_argument,       0, 0 },
-	{ "splice-flank",   required_argument, 0, 0 },
-	{ "idx-no-seq",     no_argument,       0, 0 },
-	{ "end-seed-pen",   required_argument, 0, 0 },   // 21
-	{ "for-only",       no_argument,       0, 0 },   // 22
-	{ "rev-only",       no_argument,       0, 0 },   // 23
-	{ "heap-sort",      required_argument, 0, 0 },   // 24
-	{ "all-chain",      no_argument,       0, 'P' },
-	{ "dual",           required_argument, 0, 0 },   // 26
-	{ "max-clip-ratio", required_argument, 0, 0 },   // 27
-	{ "min-occ-floor",  required_argument, 0, 0 },   // 28
-	{ "MD",             no_argument,       0, 0 },   // 29
-	{ "help",           no_argument,       0, 'h' },
-	{ "max-intron-len", required_argument, 0, 'G' },
-	{ "version",        no_argument,       0, 'V' },
-	{ "min-count",      required_argument, 0, 'n' },
-	{ "min-chain-score",required_argument, 0, 'm' },
-	{ "mask-level",     required_argument, 0, 'M' },
-	{ "min-dp-score",   required_argument, 0, 's' },
-	{ "sam",            no_argument,       0, 'a' },
-	{ 0, 0, 0, 0}
+static ko_longopt_t long_options[] = {
+	{ "bucket-bits",    ko_required_argument, 300 },
+	{ "mb-size",        ko_required_argument, 'K' },
+	{ "seed",           ko_required_argument, 302 },
+	{ "no-kalloc",      ko_no_argument,       303 },
+	{ "print-qname",    ko_no_argument,       304 },
+	{ "no-self",        ko_no_argument,       'D' },
+	{ "print-seeds",    ko_no_argument,       306 },
+	{ "max-chain-skip", ko_required_argument, 307 },
+	{ "min-dp-len",     ko_required_argument, 308 },
+	{ "print-aln-seq",  ko_no_argument,       309 },
+	{ "splice",         ko_no_argument,       310 },
+	{ "cost-non-gt-ag", ko_required_argument, 'C' },
+	{ "no-long-join",   ko_no_argument,       312 },
+	{ "sr",             ko_no_argument,       313 },
+	{ "frag",           ko_required_argument, 314 },
+	{ "secondary",      ko_required_argument, 315 },
+	{ "cs",             ko_optional_argument, 316 },
+	{ "end-bonus",      ko_required_argument, 317 },
+	{ "no-pairing",     ko_no_argument,       318 },
+	{ "splice-flank",   ko_required_argument, 319 },
+	{ "idx-no-seq",     ko_no_argument,       320 },
+	{ "end-seed-pen",   ko_required_argument, 321 },
+	{ "for-only",       ko_no_argument,       322 },
+	{ "rev-only",       ko_no_argument,       323 },
+	{ "heap-sort",      ko_required_argument, 324 },
+	{ "all-chain",      ko_no_argument,       'P' },
+	{ "dual",           ko_required_argument, 326 },
+	{ "max-clip-ratio", ko_required_argument, 327 },
+	{ "min-occ-floor",  ko_required_argument, 328 },
+	{ "MD",             ko_no_argument,       329 },
+	{ "lj-min-ratio",   ko_required_argument, 330 },
+	{ "score-N",        ko_required_argument, 331 },
+	{ "eqx",            ko_no_argument,       332 },
+	{ "paf-no-hit",     ko_no_argument,       333 },
+	{ "split-prefix",   ko_required_argument, 334 },
+	{ "no-end-flt",     ko_no_argument,       335 },
+	{ "hard-mask-level",ko_no_argument,       336 },
+	{ "bed",            ko_required_argument, 337 },
+	{ "help",           ko_no_argument,       'h' },
+	{ "max-intron-len", ko_required_argument, 'G' },
+	{ "version",        ko_no_argument,       'V' },
+	{ "min-count",      ko_required_argument, 'n' },
+	{ "min-chain-score",ko_required_argument, 'm' },
+	{ "mask-level",     ko_required_argument, 'M' },
+	{ "min-dp-score",   ko_required_argument, 's' },
+	{ "sam",            ko_no_argument,       'a' },
+	{ 0, 0, 0 }
 };
 
 static inline int64_t mm_parse_num(const char *str)
 {
 	double x;
 	char *p;
-	x = strtod(optarg, &p);
+	x = strtod(str, &p);
 	if (*p == 'G' || *p == 'g') x *= 1e9;
 	else if (*p == 'M' || *p == 'm') x *= 1e6;
 	else if (*p == 'K' || *p == 'k') x *= 1e3;
@@ -94,11 +98,12 @@ static inline void yes_or_no(mm_mapopt_t *opt, int flag, int long_idx, const cha
 
 int main(int argc, char *argv[])
 {
-	const char *opt_str = "2aSDw: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:hF:LC:y";
+	const char *opt_str = "2aSDw: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:hF:LC:yYP";
+	ketopt_t o = KETOPT_INIT;
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
-	int i, c, n_threads = 3, long_idx;
-	char *fnw = 0, *rg = 0, *s;
+	int i, c, n_threads = 3, n_parts, old_best_n = -1;
+	char *fnw = 0, *fn_bed = 0, *rg = 0, *s;
 	FILE *fp_help = stderr;
 	mm_idx_reader_t *idx_rdr;
 	mm_idx_t *mi;
@@ -108,30 +113,36 @@ int main(int argc, char *argv[])
 	mm_realtime0 = realtime();
 	mm_set_opt(0, &ipt, &opt);
 
-	while ((c = getopt_long(argc, argv, opt_str, long_options, &long_idx)) >= 0) // apply option -x/preset first
+	while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) { // test command line options and apply option -x/preset first
 		if (c == 'x') {
-			if (mm_set_opt(optarg, &ipt, &opt) < 0) {
-				fprintf(stderr, "[ERROR] unknown preset '%s'\n", optarg);
+			if (mm_set_opt(o.arg, &ipt, &opt) < 0) {
+				fprintf(stderr, "[ERROR] unknown preset '%s'\n", o.arg);
 				return 1;
 			}
-			break;
+		} else if (c == ':') {
+			fprintf(stderr, "[ERROR] missing option argument\n");
+			return 1;
+		} else if (c == '?') {
+			fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[o.i]);
+			return 1;
 		}
-	optind = 0; // for musl getopt, optind=0 has the same effect as optreset=1; older libc doesn't have optreset
+	}
+	o = KETOPT_INIT;
 
-	while ((c = getopt_long(argc, argv, opt_str, long_options, &long_idx)) >= 0) {
-		if (c == 'w') ipt.w = atoi(optarg);
-		else if (c == 'k') ipt.k = atoi(optarg);
+	while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) {
+		if (c == 'w') ipt.w = atoi(o.arg);
+		else if (c == 'k') ipt.k = atoi(o.arg);
 		else if (c == 'H') ipt.flag |= MM_I_HPC;
-		else if (c == 'd') fnw = optarg; // the above are indexing related options, except -I
-		else if (c == 'r') opt.bw = (int)mm_parse_num(optarg);
-		else if (c == 't') n_threads = atoi(optarg);
-		else if (c == 'v') mm_verbose = atoi(optarg);
-		else if (c == 'g') opt.max_gap = (int)mm_parse_num(optarg);
-		else if (c == 'G') mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(optarg));
-		else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(optarg);
-		else if (c == 'N') opt.best_n = atoi(optarg);
-		else if (c == 'p') opt.pri_ratio = atof(optarg);
-		else if (c == 'M') opt.mask_level = atof(optarg);
+		else if (c == 'd') fnw = o.arg; // the above are indexing related options, except -I
+		else if (c == 'r') opt.bw = (int)mm_parse_num(o.arg);
+		else if (c == 't') n_threads = atoi(o.arg);
+		else if (c == 'v') mm_verbose = atoi(o.arg);
+		else if (c == 'g') opt.max_gap = (int)mm_parse_num(o.arg);
+		else if (c == 'G') mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(o.arg));
+		else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(o.arg);
+		else if (c == 'N') old_best_n = opt.best_n, opt.best_n = atoi(o.arg);
+		else if (c == 'p') opt.pri_ratio = atof(o.arg);
+		else if (c == 'M') opt.mask_level = atof(o.arg);
 		else if (c == 'c') opt.flag |= MM_F_OUT_CG | MM_F_CIGAR;
 		else if (c == 'D') opt.flag |= MM_F_NO_DIAG;
 		else if (c == 'P') opt.flag |= MM_F_ALL_CHAINS;
@@ -141,59 +152,67 @@ int main(int argc, char *argv[])
 		else if (c == 'Y') opt.flag |= MM_F_SOFTCLIP;
 		else if (c == 'L') opt.flag |= MM_F_LONG_CIGAR;
 		else if (c == 'y') opt.flag |= MM_F_COPY_COMMENT;
-		else if (c == 'T') opt.sdust_thres = atoi(optarg);
-		else if (c == 'n') opt.min_cnt = atoi(optarg);
-		else if (c == 'm') opt.min_chain_score = atoi(optarg);
-		else if (c == 'A') opt.a = atoi(optarg);
-		else if (c == 'B') opt.b = atoi(optarg);
-		else if (c == 's') opt.min_dp_max = atoi(optarg);
-		else if (c == 'C') opt.noncan = atoi(optarg);
-		else if (c == 'I') ipt.batch_size = mm_parse_num(optarg);
-		else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(optarg);
-		else if (c == 'R') rg = optarg;
+		else if (c == 'T') opt.sdust_thres = atoi(o.arg);
+		else if (c == 'n') opt.min_cnt = atoi(o.arg);
+		else if (c == 'm') opt.min_chain_score = atoi(o.arg);
+		else if (c == 'A') opt.a = atoi(o.arg);
+		else if (c == 'B') opt.b = atoi(o.arg);
+		else if (c == 's') opt.min_dp_max = atoi(o.arg);
+		else if (c == 'C') opt.noncan = atoi(o.arg);
+		else if (c == 'I') ipt.batch_size = mm_parse_num(o.arg);
+		else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(o.arg);
+		else if (c == 'R') rg = o.arg;
 		else if (c == 'h') fp_help = stdout;
 		else if (c == '2') opt.flag |= MM_F_2_IO_THREADS;
-		else if (c == 0 && long_idx == 0) ipt.bucket_bits = atoi(optarg); // --bucket-bits
-		else if (c == 0 && long_idx == 2) opt.seed = atoi(optarg); // --seed
-		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 == 6) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
-		else if (c == 0 && long_idx == 7) opt.max_chain_skip = atoi(optarg); // --max-chain-skip
-		else if (c == 0 && long_idx == 8) opt.min_ksw_len = atoi(optarg); // --min-dp-len
-		else if (c == 0 && long_idx == 9) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
-		else if (c == 0 && long_idx ==10) opt.flag |= MM_F_SPLICE; // --splice
-		else if (c == 0 && long_idx ==12) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
-		else if (c == 0 && long_idx ==13) opt.flag |= MM_F_SR; // --sr
-		else if (c == 0 && long_idx ==17) opt.end_bonus = atoi(optarg); // --end-bonus
-		else if (c == 0 && long_idx ==18) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
-		else if (c == 0 && long_idx ==20) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
-		else if (c == 0 && long_idx ==21) opt.anchor_ext_shift = atoi(optarg); // --end-seed-pen
-		else if (c == 0 && long_idx ==22) opt.flag |= MM_F_FOR_ONLY; // --for-only
-		else if (c == 0 && long_idx ==23) opt.flag |= MM_F_REV_ONLY; // --rev-only
-		else if (c == 0 && long_idx ==27) opt.max_clip_ratio = atof(optarg); // --max-clip-ratio
-		else if (c == 0 && long_idx ==28) opt.min_mid_occ = atoi(optarg); // --min-occ-floor
-		else if (c == 0 && long_idx ==29) opt.flag |= MM_F_OUT_MD; // --MD
-		else if (c == 0 && long_idx == 14) { // --frag
-			yes_or_no(&opt, MM_F_FRAG_MODE, long_idx, optarg, 1);
-		} else if (c == 0 && long_idx == 15) { // --secondary
-			yes_or_no(&opt, MM_F_NO_PRINT_2ND, long_idx, optarg, 0);
-		} else if (c == 0 && long_idx == 16) { // --cs
+		else if (c == 300) ipt.bucket_bits = atoi(o.arg); // --bucket-bits
+		else if (c == 302) opt.seed = atoi(o.arg); // --seed
+		else if (c == 303) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
+		else if (c == 304) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
+		else if (c == 306) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
+		else if (c == 307) opt.max_chain_skip = atoi(o.arg); // --max-chain-skip
+		else if (c == 308) opt.min_ksw_len = atoi(o.arg); // --min-dp-len
+		else if (c == 309) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
+		else if (c == 310) opt.flag |= MM_F_SPLICE; // --splice
+		else if (c == 312) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
+		else if (c == 313) opt.flag |= MM_F_SR; // --sr
+		else if (c == 317) opt.end_bonus = atoi(o.arg); // --end-bonus
+		else if (c == 318) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
+		else if (c == 320) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
+		else if (c == 321) opt.anchor_ext_shift = atoi(o.arg); // --end-seed-pen
+		else if (c == 322) opt.flag |= MM_F_FOR_ONLY; // --for-only
+		else if (c == 323) opt.flag |= MM_F_REV_ONLY; // --rev-only
+		else if (c == 327) opt.max_clip_ratio = atof(o.arg); // --max-clip-ratio
+		else if (c == 328) opt.min_mid_occ = atoi(o.arg); // --min-occ-floor
+		else if (c == 329) opt.flag |= MM_F_OUT_MD; // --MD
+		else if (c == 330) opt.min_join_flank_ratio = atof(o.arg); // --lj-min-ratio
+		else if (c == 331) opt.sc_ambi = atoi(o.arg); // --score-N
+		else if (c == 332) opt.flag |= MM_F_EQX; // --eqx
+		else if (c == 333) opt.flag |= MM_F_PAF_NO_HIT; // --paf-no-hit
+		else if (c == 334) opt.split_prefix = o.arg; // --split-prefix
+		else if (c == 335) opt.flag |= MM_F_NO_END_FLT; // --no-end-flt
+		else if (c == 336) opt.flag |= MM_F_HARD_MLEVEL; // --hard-mask-level
+		else if (c == 337) fn_bed = o.arg; // --bed-prefer
+		else if (c == 314) { // --frag
+			yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1);
+		} else if (c == 315) { // --secondary
+			yes_or_no(&opt, MM_F_NO_PRINT_2ND, o.longidx, o.arg, 0);
+		} else if (c == 316) { // --cs
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR;
-			if (optarg == 0 || strcmp(optarg, "short") == 0) {
+			if (o.arg == 0 || strcmp(o.arg, "short") == 0) {
 				opt.flag &= ~MM_F_OUT_CS_LONG;
-			} else if (strcmp(optarg, "long") == 0) {
+			} else if (strcmp(o.arg, "long") == 0) {
 				opt.flag |= MM_F_OUT_CS_LONG;
-			} else if (strcmp(optarg, "none") == 0) {
+			} else if (strcmp(o.arg, "none") == 0) {
 				opt.flag &= ~MM_F_OUT_CS;
 			} else if (mm_verbose >= 2) {
 				fprintf(stderr, "[WARNING]\033[1;31m --cs only takes 'short' or 'long'. Invalid values are assumed to be 'short'.\033[0m\n");
 			}
-		} else if (c == 0 && long_idx == 19) { // --splice-flank
-			yes_or_no(&opt, MM_F_SPLICE_FLANK, long_idx, optarg, 1);
-		} else if (c == 0 && long_idx == 24) { // --heap-sort
-			yes_or_no(&opt, MM_F_HEAP_SORT, long_idx, optarg, 1);
-		} else if (c == 0 && long_idx == 26) { // --dual
-			yes_or_no(&opt, MM_F_NO_DUAL, long_idx, optarg, 0);
+		} else if (c == 319) { // --splice-flank
+			yes_or_no(&opt, MM_F_SPLICE_FLANK, o.longidx, o.arg, 1);
+		} else if (c == 324) { // --heap-sort
+			yes_or_no(&opt, MM_F_HEAP_SORT, o.longidx, o.arg, 1);
+		} else if (c == 326) { // --dual
+			yes_or_no(&opt, MM_F_NO_DUAL, o.longidx, o.arg, 0);
 		} else if (c == 'S') {
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR | MM_F_OUT_CS_LONG;
 			if (mm_verbose >= 2)
@@ -204,27 +223,27 @@ int main(int argc, char *argv[])
 		} else if (c == 'f') {
 			double x;
 			char *p;
-			x = strtod(optarg, &p);
+			x = strtod(o.arg, &p);
 			if (x < 1.0) opt.mid_occ_frac = x, opt.mid_occ = 0;
 			else opt.mid_occ = (int)(x + .499);
 			if (*p == ',') opt.max_occ = (int)(strtod(p+1, &p) + .499);
 		} else if (c == 'u') {
-			if (*optarg == 'b') opt.flag |= MM_F_SPLICE_FOR|MM_F_SPLICE_REV; // both strands
-			else if (*optarg == 'f') opt.flag |= MM_F_SPLICE_FOR, opt.flag &= ~MM_F_SPLICE_REV; // match GT-AG
-			else if (*optarg == 'r') opt.flag |= MM_F_SPLICE_REV, opt.flag &= ~MM_F_SPLICE_FOR; // match CT-AC (reverse complement of GT-AG)
-			else if (*optarg == 'n') opt.flag &= ~(MM_F_SPLICE_FOR|MM_F_SPLICE_REV); // don't try to match the GT-AG signal
+			if (*o.arg == 'b') opt.flag |= MM_F_SPLICE_FOR|MM_F_SPLICE_REV; // both strands
+			else if (*o.arg == 'f') opt.flag |= MM_F_SPLICE_FOR, opt.flag &= ~MM_F_SPLICE_REV; // match GT-AG
+			else if (*o.arg == 'r') opt.flag |= MM_F_SPLICE_REV, opt.flag &= ~MM_F_SPLICE_FOR; // match CT-AC (reverse complement of GT-AG)
+			else if (*o.arg == 'n') opt.flag &= ~(MM_F_SPLICE_FOR|MM_F_SPLICE_REV); // don't try to match the GT-AG signal
 			else {
 				fprintf(stderr, "[ERROR]\033[1;31m unrecognized cDNA direction\033[0m\n");
 				return 1;
 			}
 		} else if (c == 'z') {
-			opt.zdrop = opt.zdrop_inv = strtol(optarg, &s, 10);
+			opt.zdrop = opt.zdrop_inv = strtol(o.arg, &s, 10);
 			if (*s == ',') opt.zdrop_inv = strtol(s + 1, &s, 10);
 		} else if (c == 'O') {
-			opt.q = opt.q2 = strtol(optarg, &s, 10);
+			opt.q = opt.q2 = strtol(o.arg, &s, 10);
 			if (*s == ',') opt.q2 = strtol(s + 1, &s, 10);
 		} else if (c == 'E') {
-			opt.e = opt.e2 = strtol(optarg, &s, 10);
+			opt.e = opt.e2 = strtol(o.arg, &s, 10);
 			if (*s == ',') opt.e2 = strtol(s + 1, &s, 10);
 		}
 	}
@@ -236,12 +255,16 @@ int main(int argc, char *argv[])
 		ipt.flag |= MM_I_NO_SEQ;
 	if (mm_check_opt(&ipt, &opt) < 0)
 		return 1;
+	if (opt.best_n == 0) {
+		fprintf(stderr, "[WARNING]\033[1;31m changed '-N 0' to '-N %d --secondary=no'.\033[0m\n", old_best_n);
+		opt.best_n = old_best_n, opt.flag |= MM_F_NO_PRINT_2ND;
+	}
 
-	if (argc == optind || fp_help == stdout) {
+	if (argc == o.ind || fp_help == stdout) {
 		fprintf(fp_help, "Usage: minimap2 [options] <target.fa>|<target.idx> [query.fa] [...]\n");
 		fprintf(fp_help, "Options:\n");
 		fprintf(fp_help, "  Indexing:\n");
-		fprintf(fp_help, "    -H           use homopolymer-compressed k-mer\n");
+		fprintf(fp_help, "    -H           use homopolymer-compressed k-mer (preferrable for PacBio)\n");
 		fprintf(fp_help, "    -k INT       k-mer size (no larger than 28) [%d]\n", ipt.k);
 		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");
@@ -274,35 +297,33 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -c           output CIGAR in PAF\n");
 		fprintf(fp_help, "    --cs[=STR]   output the cs tag; STR is 'short' (if absent) or 'long' [none]\n");
 		fprintf(fp_help, "    --MD         output the MD tag\n");
+		fprintf(fp_help, "    --eqx        write =/X CIGAR operators\n");
 		fprintf(fp_help, "    -Y           use soft clipping for supplementary alignments\n");
 		fprintf(fp_help, "    -t INT       number of threads [%d]\n", n_threads);
 		fprintf(fp_help, "    -K NUM       minibatch size for mapping [500M]\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");
-		fprintf(fp_help, "    -x STR       preset (always applied before other options) []\n");
-		fprintf(fp_help, "                 map-pb: -Hk19 (PacBio vs reference mapping)\n");
-		fprintf(fp_help, "                 map-ont: -k15 (Oxford Nanopore vs reference mapping)\n");
-		fprintf(fp_help, "                 asm5: -k19 -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200 (asm to ref mapping; break at 5%% div.)\n");
-		fprintf(fp_help, "                 asm10: -k19 -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200 (asm to ref mapping; break at 10%% div.)\n");
-		fprintf(fp_help, "                 ava-pb: -Hk19 -Xw5 -m100 -g10000 --max-chain-skip 25 (PacBio read overlap)\n");
-		fprintf(fp_help, "                 ava-ont: -k15 -Xw5 -m100 -g10000 -r2000 --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");
+		fprintf(fp_help, "    -x STR       preset (always applied before other options; see minimap2.1 for details) []\n");
+		fprintf(fp_help, "                 - map-pb/map-ont: PacBio/Nanopore vs reference mapping\n");
+		fprintf(fp_help, "                 - ava-pb/ava-ont: PacBio/Nanopore read overlap\n");
+		fprintf(fp_help, "                 - asm5/asm10/asm20: asm-to-ref mapping, for ~0.1/1/5%% sequence divergence\n");
+		fprintf(fp_help, "                 - splice: long-read spliced alignment\n");
+		fprintf(fp_help, "                 - sr: genomic short-read mapping\n");
+		fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of these and other advanced command-line options.\n");
 		return fp_help == stdout? 0 : 1;
 	}
 
-	if ((opt.flag & MM_F_SR) && argc - optind > 3) {
+	if ((opt.flag & MM_F_SR) && argc - o.ind > 3) {
 		fprintf(stderr, "[ERROR] incorrect input: in the sr mode, please specify no more than two query files.\n");
 		return 1;
 	}
-	idx_rdr = mm_idx_reader_open(argv[optind], &ipt, fnw);
+	idx_rdr = mm_idx_reader_open(argv[o.ind], &ipt, fnw);
 	if (idx_rdr == 0) {
-		fprintf(stderr, "[ERROR] failed to open file '%s'\n", argv[optind]);
+		fprintf(stderr, "[ERROR] failed to open file '%s'\n", argv[o.ind]);
 		return 1;
 	}
-	if (!idx_rdr->is_idx && fnw == 0 && argc - optind < 2) {
+	if (!idx_rdr->is_idx && fnw == 0 && argc - o.ind < 2) {
 		fprintf(stderr, "[ERROR] missing input: please specify a query file to map or option -d to keep the index\n");
 		mm_idx_reader_close(idx_rdr);
 		return 1;
@@ -321,25 +342,30 @@ int main(int argc, char *argv[])
 				mm_write_sam_hdr(mi, rg, MM_VERSION, argc, argv);
 			} else {
 				mm_write_sam_hdr(0, rg, MM_VERSION, argc, argv);
-				if (mm_verbose >= 2)
-					fprintf(stderr, "[WARNING]\033[1;31m For a multi-part index, no @SQ lines will be outputted.\033[0m\n");
+				if (opt.split_prefix == 0 && mm_verbose >= 2)
+					fprintf(stderr, "[WARNING]\033[1;31m For a multi-part index, no @SQ lines will be outputted. Please use --split-prefix.\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 (argc != optind + 1) mm_mapopt_update(&opt, mi);
+		if (argc != o.ind + 1) mm_mapopt_update(&opt, mi);
+		if (fn_bed) mm_idx_bed_read(mi, fn_bed);
 		if (mm_verbose >= 3) mm_idx_stat(mi);
 		if (!(opt.flag & MM_F_FRAG_MODE)) {
-			for (i = optind + 1; i < argc; ++i)
+			for (i = o.ind + 1; i < argc; ++i)
 				mm_map_file(mi, argv[i], &opt, n_threads);
 		} else {
-			mm_map_file_frag(mi, argc - (optind + 1), (const char**)&argv[optind + 1], &opt, n_threads);
+			mm_map_file_frag(mi, argc - (o.ind + 1), (const char**)&argv[o.ind + 1], &opt, n_threads);
 		}
 		mm_idx_destroy(mi);
 	}
+	n_parts = idx_rdr->n_parts;
 	mm_idx_reader_close(idx_rdr);
 
+	if (opt.split_prefix)
+		mm_split_merge(argc - (o.ind + 1), (const char**)&argv[o.ind + 1], &opt, n_parts);
+
 	if (fflush(stdout) == EOF) {
 		fprintf(stderr, "[ERROR] failed to write the results\n");
 		exit(EXIT_FAILURE);
@@ -350,7 +376,7 @@ int main(int argc, char *argv[])
 		fprintf(stderr, "[M::%s] CMD:", __func__);
 		for (i = 0; i < argc; ++i)
 			fprintf(stderr, " %s", argv[i]);
-		fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec\n", __func__, realtime() - mm_realtime0, cputime());
+		fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, realtime() - mm_realtime0, cputime(), peakrss() / 1024.0 / 1024.0 / 1024.0);
 	}
 	return 0;
 }

map.c

@@ -11,6 +11,7 @@
 
 struct mm_tbuf_s {
 	void *km;
+	int rep_len, frag_gap;
 };
 
 mm_tbuf_t *mm_tbuf_init(void)
@@ -28,6 +29,11 @@ void mm_tbuf_destroy(mm_tbuf_t *b)
 	free(b);
 }
 
+void *mm_tbuf_get_km(mm_tbuf_t *b)
+{
+	return b->km;
+}
+
 static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *seq, int sdust_thres)
 {
 	int n_dreg, j, k, u = 0;
@@ -39,12 +45,12 @@ static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *se
 	for (j = k = 0; j < n; ++j) { // squeeze out minimizers that significantly overlap with LCRs
 		int32_t qpos = (uint32_t)a[j].y>>1, span = a[j].x&0xff;
 		int32_t s = qpos - (span - 1), e = s + span;
-		while (u < n_dreg && (uint32_t)dreg[u] <= s) ++u;
-		if (u < n_dreg && dreg[u]>>32 < e) {
+		while (u < n_dreg && (int32_t)dreg[u] <= s) ++u;
+		if (u < n_dreg && (int32_t)(dreg[u]>>32) < e) {
 			int v, l = 0;
-			for (v = u; v < n_dreg && dreg[v]>>32 < e; ++v) { // iterate over LCRs overlapping this minimizer
-				int ss = s > dreg[v]>>32? s : dreg[v]>>32;
-				int ee = e < (uint32_t)dreg[v]? e : (uint32_t)dreg[v];
+			for (v = u; v < n_dreg && (int32_t)(dreg[v]>>32) < e; ++v) { // iterate over LCRs overlapping this minimizer
+				int ss = s > (int32_t)(dreg[v]>>32)? s : dreg[v]>>32;
+				int ee = e < (int32_t)dreg[v]? e : (uint32_t)dreg[v];
 				l += ee - ss;
 			}
 			if (l <= span>>1) a[k++] = a[j]; // keep the minimizer if less than half of it falls in masked region
@@ -56,9 +62,10 @@ static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *se
 
 static void collect_minimizers(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, mm128_v *mv)
 {
-	int i, j, n, sum = 0;
+	int i, n, sum = 0;
 	mv->n = 0;
 	for (i = n = 0; i < n_segs; ++i) {
+		size_t j;
 		mm_sketch(km, seqs[i], qlens[i], mi->w, mi->k, i, mi->flag&MM_I_HPC, mv);
 		for (j = n; j < mv->n; ++j)
 			mv->a[j].y += sum << 1;
@@ -81,12 +88,13 @@ typedef struct {
 
 static mm_match_t *collect_matches(void *km, int *_n_m, int max_occ, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos)
 {
-	int i, rep_st = 0, rep_en = 0, n_m;
+	int rep_st = 0, rep_en = 0, n_m;
+	size_t i;
 	mm_match_t *m;
 	*n_mini_pos = 0;
 	*mini_pos = (uint64_t*)kmalloc(km, mv->n * sizeof(uint64_t));
 	m = (mm_match_t*)kmalloc(km, mv->n * sizeof(mm_match_t));
-	for (i = n_m = 0, *rep_len = 0, *n_a = 0; i < mv->n; ++i) {
+	for (i = 0, n_m = 0, *rep_len = 0, *n_a = 0; i < mv->n; ++i) {
 		const uint64_t *cr;
 		mm128_t *p = &mv->a[i];
 		uint32_t q_pos = (uint32_t)p->y, q_span = p->x & 0xff;
@@ -120,7 +128,7 @@ static inline int skip_seed(int flag, uint64_t r, const mm_match_t *q, const cha
 		const mm_idx_seq_t *s = &mi->seq[r>>32];
 		int cmp;
 		cmp = strcmp(qname, s->name);
-		if ((flag&MM_F_NO_DIAG) && cmp == 0 && s->len == qlen) {
+		if ((flag&MM_F_NO_DIAG) && cmp == 0 && (int)s->len == qlen) {
 			if ((uint32_t)r>>1 == (q->q_pos>>1)) return 1; // avoid the diagnonal anchors
 			if ((r&1) == (q->q_pos&1)) *is_self = 1; // this flag is used to avoid spurious extension on self chain
 		}
@@ -163,19 +171,20 @@ static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max
 		mm128_t *p;
 		uint64_t r = heap->x;
 		int32_t is_self, rpos = (uint32_t)r >> 1;
-		if (skip_seed(opt->flag, r, q, qname, qlen, mi, &is_self)) continue;
-		if ((r&1) == (q->q_pos&1)) { // forward strand
-			p = &a[n_for++];
-			p->x = (r&0xffffffff00000000ULL) | rpos;
-			p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
-		} else { // reverse strand
-			p = &a[(*n_a) - (++n_rev)];
-			p->x = 1ULL<<63 | (r&0xffffffff00000000ULL) | rpos;
-			p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
+		if (!skip_seed(opt->flag, r, q, qname, qlen, mi, &is_self)) {
+			if ((r&1) == (q->q_pos&1)) { // forward strand
+				p = &a[n_for++];
+				p->x = (r&0xffffffff00000000ULL) | rpos;
+				p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
+			} else { // reverse strand
+				p = &a[(*n_a) - (++n_rev)];
+				p->x = 1ULL<<63 | (r&0xffffffff00000000ULL) | rpos;
+				p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
+			}
+			p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
+			if (q->is_tandem) p->y |= MM_SEED_TANDEM;
+			if (is_self) p->y |= MM_SEED_SELF;
 		}
-		p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
-		if (q->is_tandem) p->y |= MM_SEED_TANDEM;
-		if (is_self) p->y |= MM_SEED_SELF;
 		// update the heap
 		if ((uint32_t)heap->y < q->n - 1) {
 			++heap[0].y;
@@ -205,7 +214,7 @@ static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max
 static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, const mm128_v *mv, int qlen, int64_t *n_a, int *rep_len,
 								  int *n_mini_pos, uint64_t **mini_pos)
 {
-	int i, k, n_m;
+	int i, n_m;
 	mm_match_t *m;
 	mm128_t *a;
 	m = collect_matches(km, &n_m, max_occ, mi, mv, n_a, rep_len, n_mini_pos, mini_pos);
@@ -213,6 +222,7 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ,
 	for (i = 0, *n_a = 0; i < n_m; ++i) {
 		mm_match_t *q = &m[i];
 		const uint64_t *r = q->cr;
+		uint32_t k;
 		for (k = 0; k < q->n; ++k) {
 			int32_t is_self, rpos = (uint32_t)r[k] >> 1;
 			mm128_t *p;
@@ -238,7 +248,7 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ,
 static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_idx_t *mi, void *km, int qlen, int n_segs, const int *qlens, int *n_regs, mm_reg1_t *regs, mm128_t *a)
 {
 	if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
-		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b);
+		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL);
 		if (n_segs <= 1) mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 		else mm_select_sub_multi(km, opt->pri_ratio, 0.2f, 0.7f, max_chain_gap_ref, mi->k*2, opt->best_n, n_segs, qlens, n_regs, regs);
 		if (!(opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN))) // long join not working well without primary chains
@@ -251,7 +261,7 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k
 	if (!(opt->flag & MM_F_CIGAR)) return regs;
 	regs = mm_align_skeleton(km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs()
 	if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
-		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b);
+		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL);
 		mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 		mm_set_sam_pri(*n_regs, regs);
 	}
@@ -308,10 +318,10 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		if (n_regs0 > 0) { // test if the best chain has all the segments
 			int n_chained_segs = 1, max = 0, max_i = -1, max_off = -1, off = 0;
 			for (i = 0; i < n_regs0; ++i) { // find the best chain
-				if (max < u[i]>>32) max = u[i]>>32, max_i = i, max_off = off;
+				if (max < (int)(u[i]>>32)) max = u[i]>>32, max_i = i, max_off = off;
 				off += (uint32_t)u[i];
 			}
-			for (i = 1; i < (uint32_t)u[max_i]; ++i) // count the number of segments in the best chain
+			for (i = 1; i < (int32_t)u[max_i]; ++i) // count the number of segments in the best chain
 				if ((a[max_off+i].y&MM_SEED_SEG_MASK) != (a[max_off+i-1].y&MM_SEED_SEG_MASK))
 					++n_chained_segs;
 			if (n_chained_segs < n_segs)
@@ -326,6 +336,8 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 			a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 		}
 	}
+	b->frag_gap = max_chain_gap_ref;
+	b->rep_len = rep_len;
 
 	regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a);
 
@@ -347,7 +359,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		seg = mm_seg_gen(b->km, hash, n_segs, qlens, n_regs0, regs0, n_regs, regs, a); // split fragment chain to separate segment chains
 		free(regs0);
 		for (i = 0; i < n_segs; ++i) {
-			mm_set_parent(b->km, opt->mask_level, n_regs[i], regs[i], opt->a * 2 + opt->b); // update mm_reg1_t::parent
+			mm_set_parent(b->km, opt->mask_level, n_regs[i], regs[i], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL); // update mm_reg1_t::parent
 			regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a);
 			mm_set_mapq(b->km, n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr);
 		}
@@ -390,13 +402,17 @@ typedef struct {
 	mm_bseq_file_t **fp;
 	const mm_idx_t *mi;
 	kstring_t str;
+
+	int n_parts;
+	uint32_t *rid_shift;
+	FILE *fp_split, **fp_parts;
 } pipeline_t;
 
 typedef struct {
 	const pipeline_t *p;
     int n_seq, n_frag;
 	mm_bseq1_t *seq;
-	int *n_reg, *seg_off, *n_seg;
+	int *n_reg, *seg_off, *n_seg, *rep_len, *frag_gap;
 	mm_reg1_t **reg;
 	mm_tbuf_t **buf;
 } step_t;
@@ -417,10 +433,17 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
 		qseqs[j] = s->seq[off + j].seq;
 	}
 	if (s->p->opt->flag & MM_F_INDEPEND_SEG) {
-		for (j = 0; j < s->n_seg[i]; ++j)
+		for (j = 0; j < s->n_seg[i]; ++j) {
 			mm_map_frag(s->p->mi, 1, &qlens[j], &qseqs[j], &s->n_reg[off+j], &s->reg[off+j], b, s->p->opt, s->seq[off+j].name);
+			s->rep_len[off + j] = b->rep_len;
+			s->frag_gap[off + j] = b->frag_gap;
+		}
 	} else {
 		mm_map_frag(s->p->mi, s->n_seg[i], qlens, qseqs, &s->n_reg[off], &s->reg[off], b, s->p->opt, s->seq[off].name);
+		for (j = 0; j < s->n_seg[i]; ++j) {
+			s->rep_len[off + j] = b->rep_len;
+			s->frag_gap[off + j] = b->frag_gap;
+		}
 	}
 	for (j = 0; j < s->n_seg[i]; ++j) // flip the query strand and coordinate to the original read strand
 		if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1)))) {
@@ -436,6 +459,63 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
 		}
 }
 
+static void merge_hits(step_t *s)
+{
+	int f, i, k0, k, max_seg = 0, *n_reg_part, *rep_len_part, *frag_gap_part, *qlens;
+	void *km;
+	FILE **fp = s->p->fp_parts;
+	const mm_mapopt_t *opt = s->p->opt;
+
+	km = km_init();
+	for (f = 0; f < s->n_frag; ++f)
+		max_seg = max_seg > s->n_seg[f]? max_seg : s->n_seg[f];
+	qlens = CALLOC(int, max_seg + s->p->n_parts * 3);
+	n_reg_part = qlens + max_seg;
+	rep_len_part = n_reg_part + s->p->n_parts;
+	frag_gap_part = rep_len_part + s->p->n_parts;
+	for (f = 0, k = k0 = 0; f < s->n_frag; ++f) {
+		k0 = k;
+		for (i = 0; i < s->n_seg[f]; ++i, ++k) {
+			int j, l, t, rep_len = 0;
+			qlens[i] = s->seq[k].l_seq;
+			for (j = 0, s->n_reg[k] = 0; j < s->p->n_parts; ++j) {
+				mm_err_fread(&n_reg_part[j],    sizeof(int), 1, fp[j]);
+				mm_err_fread(&rep_len_part[j],  sizeof(int), 1, fp[j]);
+				mm_err_fread(&frag_gap_part[j], sizeof(int), 1, fp[j]);
+				s->n_reg[k] += n_reg_part[j];
+				if (rep_len < rep_len_part[j])
+					rep_len = rep_len_part[j];
+			}
+			s->reg[k] = CALLOC(mm_reg1_t, s->n_reg[k]);
+			for (j = 0, l = 0; j < s->p->n_parts; ++j) {
+				for (t = 0; t < n_reg_part[j]; ++t, ++l) {
+					mm_reg1_t *r = &s->reg[k][l];
+					uint32_t capacity;
+					mm_err_fread(r, sizeof(mm_reg1_t), 1, fp[j]);
+					r->rid += s->p->rid_shift[j];
+					if (opt->flag & MM_F_CIGAR) {
+						mm_err_fread(&capacity, 4, 1, fp[j]);
+						r->p = (mm_extra_t*)calloc(capacity, 4);
+						r->p->capacity = capacity;
+						mm_err_fread(r->p, r->p->capacity, 4, fp[j]);
+					}
+				}
+			}
+			mm_hit_sort(km, &s->n_reg[k], s->reg[k]);
+			mm_set_parent(km, opt->mask_level, s->n_reg[k], s->reg[k], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL);
+			if (!(opt->flag & MM_F_ALL_CHAINS)) {
+				mm_select_sub(km, opt->pri_ratio, s->p->mi->k*2, opt->best_n, &s->n_reg[k], s->reg[k]);
+				mm_set_sam_pri(s->n_reg[k], s->reg[k]);
+			}
+			mm_set_mapq(km, s->n_reg[k], s->reg[k], opt->min_chain_score, opt->a, rep_len, !!(opt->flag & MM_F_SR));
+		}
+		if (s->n_seg[f] == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR))
+			mm_pair(km, frag_gap_part[0], opt->pe_bonus, opt->a * 2 + opt->b, opt->a, qlens, &s->n_reg[k0], &s->reg[k0]);
+	}
+	free(qlens);
+	km_destroy(km);
+}
+
 static void *worker_pipeline(void *shared, int step, void *in)
 {
 	int i, j, k;
@@ -455,9 +535,11 @@ static void *worker_pipeline(void *shared, int step, void *in)
 			s->buf = (mm_tbuf_t**)calloc(p->n_threads, sizeof(mm_tbuf_t*));
 			for (i = 0; i < p->n_threads; ++i)
 				s->buf[i] = mm_tbuf_init();
-			s->n_reg = (int*)calloc(3 * s->n_seq, sizeof(int));
-			s->seg_off = s->n_reg + s->n_seq; // seg_off and n_seg are allocated together with n_reg
+			s->n_reg = (int*)calloc(5 * s->n_seq, sizeof(int));
+			s->seg_off = s->n_reg + s->n_seq; // seg_off, n_seg, rep_len and frag_gap are allocated together with n_reg
 			s->n_seg = s->seg_off + s->n_seq;
+			s->rep_len = s->n_seg + s->n_seq;
+			s->frag_gap = s->rep_len + s->n_seq;
 			s->reg = (mm_reg1_t**)calloc(s->n_seq, sizeof(mm_reg1_t*));
 			for (i = 1, j = 0; i <= s->n_seq; ++i)
 				if (i == s->n_seq || !frag_mode || !mm_qname_same(s->seq[i-1].name, s->seq[i].name)) {
@@ -468,7 +550,8 @@ static void *worker_pipeline(void *shared, int step, void *in)
 			return s;
 		} else free(s);
     } else if (step == 1) { // step 1: map
-		kt_for(p->n_threads, worker_for, in, ((step_t*)in)->n_frag);
+		if (p->n_parts > 0) merge_hits((step_t*)in);
+		else kt_for(p->n_threads, worker_for, in, ((step_t*)in)->n_frag);
 		return in;
     } else if (step == 2) { // step 2: output
 		void *km = 0;
@@ -481,19 +564,35 @@ static void *worker_pipeline(void *shared, int step, void *in)
 			int seg_st = s->seg_off[k], seg_en = s->seg_off[k] + s->n_seg[k];
 			for (i = seg_st; i < seg_en; ++i) {
 				mm_bseq1_t *t = &s->seq[i];
-				for (j = 0; j < s->n_reg[i]; ++j) {
-					mm_reg1_t *r = &s->reg[i][j];
-					assert(!r->sam_pri || r->id == r->parent);
-					if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent)
-						continue;
+				if (p->opt->split_prefix && p->n_parts == 0) { // then write to temporary files
+					mm_err_fwrite(&s->n_reg[i],    sizeof(int), 1, p->fp_split);
+					mm_err_fwrite(&s->rep_len[i],  sizeof(int), 1, p->fp_split);
+					mm_err_fwrite(&s->frag_gap[i], sizeof(int), 1, p->fp_split);
+					for (j = 0; j < s->n_reg[i]; ++j) {
+						mm_reg1_t *r = &s->reg[i][j];
+						mm_err_fwrite(r, sizeof(mm_reg1_t), 1, p->fp_split);
+						if (p->opt->flag & MM_F_CIGAR) {
+							mm_err_fwrite(&r->p->capacity, 4, 1, p->fp_split);
+							mm_err_fwrite(r->p, r->p->capacity, 4, p->fp_split);
+						}
+					}
+				} else if (s->n_reg[i] > 0) { // the query has at least one hit
+					for (j = 0; j < s->n_reg[i]; ++j) {
+						mm_reg1_t *r = &s->reg[i][j];
+						assert(!r->sam_pri || r->id == r->parent);
+						if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent)
+							continue;
+						if (p->opt->flag & MM_F_OUT_SAM)
+							mm_write_sam2(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
+						else
+							mm_write_paf(&p->str, mi, t, r, km, p->opt->flag);
+						mm_err_puts(p->str.s);
+					}
+				} else if (p->opt->flag & (MM_F_OUT_SAM|MM_F_PAF_NO_HIT)) { // output an empty hit, if requested
 					if (p->opt->flag & MM_F_OUT_SAM)
-						mm_write_sam2(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
+						mm_write_sam2(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
 					else
-						mm_write_paf(&p->str, mi, t, r, km, p->opt->flag);
-					mm_err_puts(p->str.s);
-				}
-				if (s->n_reg[i] == 0 && (p->opt->flag & MM_F_OUT_SAM)) { // write an unmapped record
-					mm_write_sam2(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
+						mm_write_paf(&p->str, mi, t, 0, 0, p->opt->flag);
 					mm_err_puts(p->str.s);
 				}
 			}
@@ -502,9 +601,10 @@ static void *worker_pipeline(void *shared, int step, void *in)
 				free(s->reg[i]);
 				free(s->seq[i].seq); free(s->seq[i].name);
 				if (s->seq[i].qual) free(s->seq[i].qual);
+				if (s->seq[i].comment) free(s->seq[i].comment);
 			}
 		}
-		free(s->reg); free(s->n_reg); free(s->seq); // seg_off and n_seg were allocated with reg; no memory leak here
+		free(s->reg); free(s->n_reg); free(s->seq); // seg_off, n_seg, rep_len and frag_gap were allocated with reg; no memory leak here
 		km_destroy(km);
 		if (mm_verbose >= 3)
 			fprintf(stderr, "[M::%s::%.3f*%.2f] mapped %d sequences\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), s->n_seq);
@@ -513,32 +613,44 @@ static void *worker_pipeline(void *shared, int step, void *in)
     return 0;
 }
 
+static mm_bseq_file_t **open_bseqs(int n, const char **fn)
+{
+	mm_bseq_file_t **fp;
+	int i, j;
+	fp = (mm_bseq_file_t**)calloc(n, sizeof(mm_bseq_file_t*));
+	for (i = 0; i < n; ++i) {
+		if ((fp[i] = mm_bseq_open(fn[i])) == 0) {
+			if (mm_verbose >= 1)
+				fprintf(stderr, "ERROR: failed to open file '%s'\n", fn[i]);
+			for (j = 0; j < i; ++j)
+				mm_bseq_close(fp[j]);
+			free(fp);
+			return 0;
+		}
+	}
+	return fp;
+}
+
 int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_mapopt_t *opt, int n_threads)
 {
-	int i, j, pl_threads;
+	int i, pl_threads;
 	pipeline_t pl;
 	if (n_segs < 1) return -1;
 	memset(&pl, 0, sizeof(pipeline_t));
 	pl.n_fp = n_segs;
-	pl.fp = (mm_bseq_file_t**)calloc(n_segs, sizeof(mm_bseq_file_t*));
-	for (i = 0; i < n_segs; ++i) {
-		pl.fp[i] = mm_bseq_open(fn[i]);
-		if (pl.fp[i] == 0) {
-			if (mm_verbose >= 1)
-				fprintf(stderr, "ERROR: failed to open file '%s'\n", fn[i]);
-			for (j = 0; j < i; ++j)
-				mm_bseq_close(pl.fp[j]);
-			free(pl.fp);
-			return -1;
-		}
-	}
+	pl.fp = open_bseqs(pl.n_fp, fn);
+	if (pl.fp == 0) return -1;
 	pl.opt = opt, pl.mi = idx;
 	pl.n_threads = n_threads > 1? n_threads : 1;
 	pl.mini_batch_size = opt->mini_batch_size;
+	if (opt->split_prefix)
+		pl.fp_split = mm_split_init(opt->split_prefix, idx);
 	pl_threads = n_threads == 1? 1 : (opt->flag&MM_F_2_IO_THREADS)? 3 : 2;
 	kt_pipeline(pl_threads, worker_pipeline, &pl, 3);
+
 	free(pl.str.s);
-	for (i = 0; i < n_segs; ++i)
+	if (pl.fp_split) fclose(pl.fp_split);
+	for (i = 0; i < pl.n_fp; ++i)
 		mm_bseq_close(pl.fp[i]);
 	free(pl.fp);
 	return 0;
@@ -548,3 +660,48 @@ int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int
 {
 	return mm_map_file_frag(idx, 1, &fn, opt, n_threads);
 }
+
+int mm_split_merge(int n_segs, const char **fn, const mm_mapopt_t *opt, int n_split_idx)
+{
+	int i;
+	pipeline_t pl;
+	mm_idx_t *mi;
+	if (n_segs < 1 || n_split_idx < 1) return -1;
+	memset(&pl, 0, sizeof(pipeline_t));
+	pl.n_fp = n_segs;
+	pl.fp = open_bseqs(pl.n_fp, fn);
+	if (pl.fp == 0) return -1;
+	pl.opt = opt;
+	pl.mini_batch_size = opt->mini_batch_size;
+
+	pl.n_parts = n_split_idx;
+	pl.fp_parts  = CALLOC(FILE*, pl.n_parts);
+	pl.rid_shift = CALLOC(uint32_t, pl.n_parts);
+	pl.mi = mi = mm_split_merge_prep(opt->split_prefix, n_split_idx, pl.fp_parts, pl.rid_shift);
+	if (pl.mi == 0) {
+		free(pl.fp_parts);
+		free(pl.rid_shift);
+		return -1;
+	}
+	for (i = n_split_idx - 1; i > 0; --i)
+		pl.rid_shift[i] = pl.rid_shift[i - 1];
+	for (pl.rid_shift[0] = 0, i = 1; i < n_split_idx; ++i)
+		pl.rid_shift[i] += pl.rid_shift[i - 1];
+	if (opt->flag & MM_F_OUT_SAM)
+		for (i = 0; i < (int32_t)pl.mi->n_seq; ++i)
+			printf("@SQ\tSN:%s\tLN:%d\n", pl.mi->seq[i].name, pl.mi->seq[i].len);
+
+	kt_pipeline(2, worker_pipeline, &pl, 3);
+
+	free(pl.str.s);
+	mm_idx_destroy(mi);
+	free(pl.rid_shift);
+	for (i = 0; i < n_split_idx; ++i)
+		fclose(pl.fp_parts[i]);
+	free(pl.fp_parts);
+	for (i = 0; i < pl.n_fp; ++i)
+		mm_bseq_close(pl.fp[i]);
+	free(pl.fp);
+	mm_split_rm_tmp(opt->split_prefix, n_split_idx);
+	return 0;
+}

minimap.h

@@ -31,6 +31,10 @@
 #define MM_F_ALL_CHAINS    0x800000
 #define MM_F_OUT_MD        0x1000000
 #define MM_F_COPY_COMMENT  0x2000000
+#define MM_F_EQX           0x4000000 // use =/X instead of M
+#define MM_F_PAF_NO_HIT    0x8000000 // output unmapped reads to PAF
+#define MM_F_NO_END_FLT    0x10000000
+#define MM_F_HARD_MLEVEL   0x20000000
 
 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2
@@ -55,12 +59,21 @@ typedef struct {
 	uint32_t len;    // length
 } mm_idx_seq_t;
 
+typedef struct {
+	uint64_t x;
+	int32_t end, idx;
+	int32_t score; // NB: wasting 4 bytes due to memory alignment
+} mm_idx_bed_t;
+
 typedef struct {
 	int32_t b, w, k, flag;
 	uint32_t n_seq;            // number of reference sequences
+	int32_t index;
+	uint32_t n_R;
 	mm_idx_seq_t *seq;         // sequence name, length and offset
 	uint32_t *S;               // 4-bit packed sequence
 	struct mm_idx_bucket_s *B; // index (hidden)
+	mm_idx_bed_t *R;
 	void *km, *h;
 } mm_idx_t;
 
@@ -115,8 +128,10 @@ typedef struct {
 
 	int max_join_long, max_join_short;
 	int min_join_flank_sc;
+	float min_join_flank_ratio;
 
 	int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties
+	int sc_ambi; // score when one or both bases are "N"
 	int noncan;      // cost of non-canonical splicing sites
 	int zdrop, zdrop_inv;   // break alignment if alignment score drops too fast along the diagonal
 	int end_bonus;
@@ -132,6 +147,8 @@ typedef struct {
 	int32_t mid_occ;     // ignore seeds with occurrences above this threshold
 	int32_t max_occ;
 	int mini_batch_size; // size of a batch of query bases to process in parallel
+
+	const char *split_prefix;
 } mm_mapopt_t;
 
 // index reader
@@ -216,6 +233,36 @@ void mm_idx_reader_close(mm_idx_reader_t *r);
 
 int mm_idx_reader_eof(const mm_idx_reader_t *r);
 
+/**
+ * Check whether the file contains a minimap2 index
+ *
+ * @param fn         file name
+ *
+ * @return the file size if fn is an index file; 0 if fn is not.
+ */
+int64_t mm_idx_is_idx(const char *fn);
+
+/**
+ * Load a part of an index
+ *
+ * Given a uni-part index, this function loads the entire index into memory.
+ * Given a multi-part index, it loads one part only and places the file pointer
+ * at the end of that part.
+ *
+ * @param fp         pointer to FILE object
+ *
+ * @return minimap2 index read from fp
+ */
+mm_idx_t *mm_idx_load(FILE *fp);
+
+/**
+ * Append an index (or one part of a full index) to file
+ *
+ * @param fp         pointer to FILE object
+ * @param mi         minimap2 index
+ */
+void mm_idx_dump(FILE *fp, const mm_idx_t *mi);
+
 /**
  * Create an index from strings in memory
  *
@@ -264,6 +311,8 @@ mm_tbuf_t *mm_tbuf_init(void);
  */
 void mm_tbuf_destroy(mm_tbuf_t *b);
 
+void *mm_tbuf_get_km(mm_tbuf_t *b);
+
 /**
  * Align a query sequence against an index
  *
@@ -300,11 +349,32 @@ int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int
 
 int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_mapopt_t *opt, int n_threads);
 
+/**
+ * Generate the cs tag (new in 2.12)
+ *
+ * @param km         memory blocks; set to NULL if unsure
+ * @param buf        buffer to write the cs/MD tag; typicall NULL on the first call
+ * @param max_len    max length of the buffer; typically set to 0 on the first call
+ * @param mi         index
+ * @param r          alignment
+ * @param seq        query sequence
+ * @param no_iden    true to use : instead of =
+ *
+ * @return the length of cs
+ */
+int mm_gen_cs(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int no_iden);
+int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq);
+
 // query sequence name and sequence in the minimap2 index
 int mm_idx_index_name(mm_idx_t *mi);
 int mm_idx_name2id(const mm_idx_t *mi, const char *name);
 int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
 
+// BED operations
+int mm_idx_bed_read(mm_idx_t *mi, const char *fn);
+int mm_idx_bed_attach(mm_idx_t *mi, uint32_t n, mm_idx_bed_t *r);
+int mm_idx_bed_query(const mm_idx_t *mi, uint64_t x);
+
 // 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 flag, int n_threads);

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "27 March 2018" "minimap2-2.10 (r761)" "Bioinformatics tools"
+.TH minimap2 1 "11 October 2018" "minimap2-2.13 (r850)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -193,7 +193,7 @@ Primarily used for all-vs-all read overlapping.
 .BI -p \ FLOAT
 Minimal secondary-to-primary score ratio to output secondary mappings [0.8].
 Between two chains overlaping over half of the shorter chain (controlled by
-.BR --mask-level ),
+.BR -M ),
 the chain with a lower score is secondary to the chain with a higher score.
 If the ratio of the scores is below
 .IR FLOAT ,
@@ -226,6 +226,11 @@ Mark as secondary a chain that overlaps with a better chain by
 .I FLOAT
 or more of the shorter chain [0.5]
 .TP
+.B --hard-mask-level
+Honor option
+.B -M
+and disable a heurstic to save unmapped subsequences.
+.TP
 .BI --max-chain-skip \ INT
 A heuristics that stops chaining early [50]. Minimap2 uses dynamic programming
 for chaining. The time complexity is quadratic in the number of seeds. This
@@ -239,6 +244,10 @@ Disable the long gap patching heuristic. When this option is applied, the
 maximum alignment gap is mostly controlled by
 .BR -r .
 .TP
+.BI --lj-min-ratio \ FLOAT
+Fraction of query sequence length required to bridge a long gap [0.5]. A
+smaller value helps to recover longer gaps, at the cost of more false gaps.
+.TP
 .B --splice
 Enable the splice alignment mode.
 .TP
@@ -248,6 +257,9 @@ applies a second round of chaining with a higher minimizer occurrence threshold
 if no good chain is found. In addition, minimap2 attempts to patch gaps between
 seeds with ungapped alignment.
 .TP
+.BI --split-prefix \ STR
+Prefix to create temporary files. Typically used for a multi-part index.
+.TP
 .BR --frag = no | yes
 Whether to enable the fragment mode [no]
 .TP
@@ -322,6 +334,9 @@ no attempt to match GT-AG [n]
 .BI --end-bonus \ INT
 Score bonus when alignment extends to the end of the query sequence [0].
 .TP
+.BI --score-N \ INT
+Score of a mismatch involving ambiguous bases [1].
+.TP
 .BR --splice-flank = yes | no
 Assume the next base to a
 .B GT
@@ -350,6 +365,10 @@ the length of the terminal gap in the chain. This option is only effective
 with
 .BR --splice .
 It helps to avoid tiny terminal exons. [6]
+.TP
+.B --no-end-flt
+Don't filter seeds towards the ends of chains before performing base-level
+alignment.
 .SS Input/output options
 .TP 10
 .B -a
@@ -394,6 +413,9 @@ is assumed. [none]
 .B --MD
 Output the MD tag (see the SAM spec).
 .TP
+.B --eqx
+Output =/X CIGAR operators for sequence match/mismatch.
+.TP
 .B -Y
 In SAM output, use soft clipping for supplementary alignments.
 .TP
@@ -572,6 +594,7 @@ nn	i	Number of ambiguous bases in the alignment
 ts	A	Transcript strand (splice mode only)
 cg	Z	CIGAR string (only in PAF)
 cs	Z	Difference string
+dv	f	Approximate per-base sequence divergence
 .TE
 
 .PP

misc.c

@@ -131,6 +131,26 @@ void mm_err_puts(const char *str)
 	}
 }
 
+void mm_err_fwrite(const void *p, size_t size, size_t nitems, FILE *fp)
+{
+	int ret;
+	ret = fwrite(p, size, nitems, fp);
+	if (ret == EOF) {
+		fprintf(stderr, "[ERROR] failed to write data\n");
+		exit(EXIT_FAILURE);
+	}
+}
+
+void mm_err_fread(void *p, size_t size, size_t nitems, FILE *fp)
+{
+	int ret;
+	ret = fread(p, size, nitems, fp);
+	if (ret == EOF) {
+		fprintf(stderr, "[ERROR] failed to read data\n");
+		exit(EXIT_FAILURE);
+	}
+}
+
 #include "ksort.h"
 
 #define sort_key_128x(a) ((a).x)

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = 'r755';
+var paftools_version = '2.13-r850';
 
 /*****************************
  ***** Library functions *****
@@ -340,14 +340,16 @@ function paf_liftover(args)
 function paf_call(args)
 {
 	var re_cs = /([:=*+-])(\d+|[A-Za-z]+)/g, re_tag = /\t(\S\S:[AZif]):(\S+)/g;
-	var c, min_cov_len = 10000, min_var_len = 50000, gap_thres = 50, min_mapq = 5;
-	var fa_tmp = null, fa, fa_lens, is_vcf = false;
-	while ((c = getopt(args, "l:L:g:q:B:f:")) != null) {
+	var c, min_cov_len = 10000, min_var_len = 50000, gap_thres = 50, gap_thres_long = 1000, min_mapq = 5;
+	var fa_tmp = null, fa, fa_lens, is_vcf = false, sample_name = "sample";
+	while ((c = getopt(args, "l:L:g:q:B:f:s:")) != null) {
 		if (c == 'l') min_cov_len = parseInt(getopt.arg);
 		else if (c == 'L') min_var_len = parseInt(getopt.arg);
 		else if (c == 'g') gap_thres = parseInt(getopt.arg);
+		else if (c == 'G') gap_thres_long = parseInt(getopt.arg);
 		else if (c == 'q') min_mapq = parseInt(getopt.arg);
 		else if (c == 'f') fa_tmp = fasta_read(getopt.arg, fa_lens);
+		else if (c == 's') sample_name = getopt.arg;
 	}
 	if (fa_tmp != null) fa = fa_tmp[0], fa_lens = fa_tmp[1], is_vcf = true;
 
@@ -359,12 +361,13 @@ function paf_call(args)
 		print("  -q INT    min mapping quality ["+min_mapq+"]");
 		print("  -g INT    short/long gap threshold (for statistics only) ["+gap_thres+"]");
 		print("  -f FILE   reference sequences (enabling VCF output) [null]");
+		print("  -s NAME   sample name in VCF header ["+sample_name+"]");
 		exit(1);
 	}
 
 	var file = args[getopt.ind] == '-'? new File() : new File(args[getopt.ind]);
 	var buf = new Bytes();
-	var tot_len = 0, n_sub = [0, 0, 0], n_ins = [0, 0, 0, 0], n_del = [0, 0, 0, 0];
+	var tot_len = 0, n_sub = [0, 0, 0], n_ins = [0, 0, 0, 0, 0], n_del = [0, 0, 0, 0, 0];
 
 	function print_vcf(o, fa)
 	{
@@ -396,13 +399,15 @@ function paf_call(args)
 			if (l == 1) ++n_ins[0];
 			else if (l == 2) ++n_ins[1];
 			else if (l < gap_thres) ++n_ins[2];
-			else ++n_ins[3];
+			else if (l < gap_thres_long) ++n_ins[3];
+			else ++n_ins[4];
 		} else if (o[6] == '-') { // deletion
 			var l = o[5].length;
 			if (l == 1) ++n_del[0];
 			else if (l == 2) ++n_del[1];
 			else if (l < gap_thres) ++n_del[2];
-			else ++n_del[3];
+			else if (l < gap_thres_long) ++n_del[3];
+			else ++n_del[4];
 		} else {
 			++n_sub[0];
 			var s = (o[5] + o[6]).toLowerCase();
@@ -420,7 +425,7 @@ function paf_call(args)
 		print('##INFO=<ID=QSTART,Number=1,Type=Integer,Description="Query start">');
 		print('##INFO=<ID=QSTRAND,Number=1,Type=String,Description="Query strand">');
 		print('##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">');
-		print('#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	sample');
+		print('#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	'+sample_name);
 	}
 
 	var a = [], out = [];
@@ -484,7 +489,7 @@ function paf_call(args)
 		// drop alignments that don't overlap with the current one
 		var k = 0;
 		for (var i = 0; i < a.length; ++i)
-			if (a[0][0] == ctg && a[0][2] > x)
+			if (a[i][0] == ctg && a[i][2] > x)
 				a[k++] = a[i];
 		a.length = k;
 		// core loop
@@ -496,7 +501,7 @@ function paf_call(args)
 				var cov = 1;
 				if (m[1] == '*' || m[1] == '+' || m[1] == '-')
 					for (var i = 0; i < a.length; ++i)
-						if (a[0][2] > x) ++cov;
+						if (a[i][2] > x) ++cov;
 				var qs, qe;
 				if (m[1] == '=' || m[1] == ':') {
 					var l = m[1] == '='? m[2].length : parseInt(m[2]);
@@ -547,14 +552,196 @@ function paf_call(args)
 	warn(n_ins[1] + " 2bp insertions");
 	warn(n_del[2] + " [3,"+gap_thres+") deletions");
 	warn(n_ins[2] + " [3,"+gap_thres+") insertions");
-	warn(n_del[3] + " >="+gap_thres+" deletions");
-	warn(n_ins[3] + " >="+gap_thres+" insertions");
+	warn(n_del[3] + " ["+gap_thres+","+gap_thres_long+") deletions");
+	warn(n_ins[3] + " ["+gap_thres+","+gap_thres_long+") insertions");
+	warn(n_del[4] + " >=" + gap_thres_long + " deletions");
+	warn(n_ins[4] + " >=" + gap_thres_long + " insertions");
 
 	buf.destroy();
 	file.close();
 	if (fa != null) fasta_free(fa);
 }
 
+function paf_asmstat(args)
+{
+	var c, min_seg_len = 10000, max_diff = 0.01;
+	while ((c = getopt(args, "l:d:")) != null) {
+		if (c == 'l') min_seg_len = parseInt(getopt.arg);
+		else if (c == 'd') max_diff = parseFloat(getopt.arg);
+	}
+	if (getopt.ind == args.length) {
+		print("Usage: paftools.js asmstat [options] <ref.fa.fai> <asm1.paf> [...]");
+		print("Options:");
+		print("  -l INT     min alignment block length [" + min_seg_len + "]");
+		print("  -d FLOAT   max gap-compressed sequence divergence [" + max_diff + "]");
+		exit(1);
+	}
+
+	var file, buf = new Bytes();
+
+	var ref_len = 0;
+	file = new File(args[getopt.ind]);
+	while (file.readline(buf) >= 0) {
+		var t = buf.toString().split("\t");
+		ref_len += parseInt(t[1]);
+	}
+	file.close();
+
+	function process_query(qblocks, qblock_len, bp) {
+		qblocks.sort(function(a,b) { return a[0]-b[0]; });
+		var last_k = null, last_blen = null, st = -1, en = -1, qcov = 0;
+		for (var k = 0; k < qblocks.length; ++k) {
+			var blen = qblocks[k][1] - qblocks[k][0];
+			if (k > 0 && qblocks[k][0] < qblocks[k-1][1]) {
+				if (qblocks[k][1] < qblocks[k-1][1]) continue;
+				blen = qblocks[k][1] - qblocks[k-1][1];
+			}
+			qblock_len.push(blen);
+			if (qblocks[k][0] > en) {
+				qcov += en - st;
+				st = qblocks[k][0];
+				en = qblocks[k][1];
+			} else en = en > qblocks[k][1]? en : qblocks[k][1];
+			if (last_k != null) {
+				var gap = 1000000000;
+				if (qblocks[k][2] == qblocks[last_k][2] && qblocks[k][3] == qblocks[last_k][3]) { // same chr and strand
+					var g1 = qblocks[k][0] - qblocks[last_k][1];
+					var g2 = qblocks[k][2] == '+'? qblocks[k][4] - qblocks[last_k][5] : qblocks[last_k][4] - qblocks[k][5];
+					gap = g1 > g2? g1 - g2 : g2 - g1;
+				}
+				var min = blen < last_blen? blen : last_blen;
+				var flank = k == 0? min : blen;
+				bp.push([flank, gap]);
+			}
+			last_k = k, last_blen = blen;
+		}
+		qcov += en - st;
+		return qcov;
+	}
+
+	function N50(lens, tot, quantile) {
+		lens.sort(function(a,b) { return b - a; });
+		if (tot == null) {
+			tot = 0;
+			for (var k = 0; k < lens.length; ++k)
+				tot += lens[k];
+		}
+		var sum = 0;
+		for (var k = 0; k < lens.length; ++k) {
+			if (sum <= quantile * tot && sum + lens[k] > quantile * tot)
+				return lens[k];
+			sum += lens[k];
+		}
+	}
+
+	function count_bp(bp, min_blen, min_gap) {
+		var n_bp = 0;
+		for (var k = 0; k < bp.length; ++k)
+			if (bp[k][0] >= min_blen && bp[k][1] >= min_gap)
+				++n_bp;
+		return n_bp;
+	}
+
+	function compute_diff(cigar, NM) {
+		var m, re = /(\d+)([MID])/g;
+		var n_M = 0, n_gapo = 0, n_gaps = 0;
+		while ((m = re.exec(cigar)) != null) {
+			var len = parseInt(m[1]);
+			if (m[2] == 'M') n_M += len;
+			else ++n_gapo, n_gaps += len;
+		}
+		if (NM < n_gaps) throw Error('NM is smaller the number of gaps');
+		return (NM - n_gaps + n_gapo) / (n_M + n_gapo);
+	}
+
+	var labels = ['Length', 'NG50', 'Coverage', 'Qcov', 'NGA50', '#breaks', 'bp(' + min_seg_len + ',0)', 'bp(' + min_seg_len + ',10k)'];
+	var rst = [];
+	for (var i = 0; i < labels.length; ++i)
+		rst[i] = [];
+
+	var n_asm = args.length - (getopt.ind + 1);
+	var header = ["Metric"];
+	for (var i = 0; i < n_asm; ++i) {
+		var n_breaks = 0, qcov = 0;
+		var fn = args[getopt.ind + 1 + i];
+		header.push(fn.replace(/.paf(.gz)?$/, ""));
+		var ref_blocks = [], qblock_len = [], qblocks = [], bp = [];
+		var query = {};
+		var last_qname = null;
+		file = new File(fn);
+		while (file.readline(buf) >= 0) {
+			var m, line = buf.toString();
+			var t = line.split("\t");
+			t[1] = parseInt(t[1]);
+			if (t.length >= 2) query[t[0]] = t[1];
+			if (t.length < 9) continue;
+			if (!/\ttp:A:[PI]/.test(line)) continue;
+			if ((m = /\tcg:Z:(\S+)/.exec(line)) == null) continue;
+			var cigar = m[1];
+			if ((m = /\tNM:i:(\d+)/.exec(line)) == null) continue;
+			var NM = parseInt(m[1]);
+			var diff = compute_diff(cigar, NM);
+			t[2] = parseInt(t[2]);
+			t[3] = parseInt(t[3]);
+			t[7] = parseInt(t[7]);
+			t[8] = parseInt(t[8]);
+			if (t[0] == last_qname) ++n_breaks;
+			if (diff > max_diff) continue;
+			if (t[3] - t[2] < min_seg_len) continue;
+			if (t[0] != last_qname) {
+				if (last_qname != null)
+					qcov += process_query(qblocks, qblock_len, bp);
+				qblocks = [];
+				last_qname = t[0];
+			}
+			ref_blocks.push([t[5], t[7], t[8]]);
+			qblocks.push([t[2], t[3], t[4], t[5], t[7], t[8]]);
+		}
+		if (last_qname != null)
+			qcov += process_query(qblocks, qblock_len, bp);
+		file.close();
+
+		// compute NG50
+		var asm_len = 0, asm_lens = []
+		for (var ctg in query) {
+			asm_len += query[ctg];
+			asm_lens.push(query[ctg]);
+		}
+		rst[0][i] = asm_len;
+		rst[1][i] = N50(asm_lens, ref_len, 0.5);
+
+		// compute coverage
+		var l_cov = 0;
+		ref_blocks.sort(function(a, b) { return a[0] > b[0]? 1 : a[0] < b[0]? -1 : a[1] - b[1]; });
+		var last_ref = null, st = -1, en = -1;
+		for (var j = 0; j < ref_blocks.length; ++j) {
+			if (ref_blocks[j][0] != last_ref || ref_blocks[j][1] > en) {
+				l_cov += en - st;
+				last_ref = ref_blocks[j][0];
+				st = ref_blocks[j][1];
+				en = ref_blocks[j][2];
+			} else en = en > ref_blocks[j][2]? en : ref_blocks[j][2];
+		}
+		l_cov += en - st;
+		rst[2][i] = (100.0 * (l_cov / ref_len)).toFixed(2) + '%';
+		rst[3][i] = (100.0 * (qcov  / asm_len)).toFixed(2) + '%';
+
+		// compute NGA50
+		rst[4][i] = N50(qblock_len, ref_len, 0.5);
+
+		// compute break points
+		rst[5][i] = n_breaks;
+		rst[6][i] = count_bp(bp, 500, 0);
+		rst[7][i] = count_bp(bp, 500, 10000);
+	}
+
+	print(header.join("\t"));
+	for (var i = 0; i < labels.length; ++i)
+		print(labels[i], rst[i].join("\t"));
+
+	buf.destroy();
+}
+
 function paf_stat(args)
 {
 	var c, gap_out_len = null;
@@ -676,8 +863,10 @@ function paf_stat(args)
 			last_qlen = ori_qlen;
 		}
 	}
-	l_tot += last_qlen;
-	l_cov += cov_len(regs);
+	if (regs.length) {
+		l_tot += last_qlen;
+		l_cov += cov_len(regs);
+	}
 
 	file.close();
 	buf.destroy();
@@ -1999,6 +2188,7 @@ function main(args)
 		print("  gff2bed    convert GTF/GFF3 to BED12");
 		print("");
 		print("  stat       collect basic mapping information in PAF/SAM");
+		print("  asmstat    collect basic assembly information");
 		print("  liftover   simplistic liftOver");
 		print("  call       call variants from asm-to-ref alignment with the cs tag");
 		print("  bedcov     compute the number of bases covered");
@@ -2019,6 +2209,7 @@ function main(args)
 	else if (cmd == 'splice2bed') paf_splice2bed(args);
 	else if (cmd == 'gff2bed') paf_gff2bed(args);
 	else if (cmd == 'stat') paf_stat(args);
+	else if (cmd == 'asmstat') paf_asmstat(args);
 	else if (cmd == 'liftover' || cmd == 'liftOver') paf_liftover(args);
 	else if (cmd == 'call') paf_call(args);
 	else if (cmd == 'mapeval') paf_mapeval(args);

mmpriv.h

@@ -28,6 +28,16 @@
 #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)
 
+#define MALLOC(type, len) ((type*)malloc((len) * sizeof(type)))
+#define CALLOC(type, len) ((type*)calloc((len), sizeof(type)))
+
+#define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr))))
+
+#define EXPAND(a, m) do { \
+		(m) = (m)? (m) + ((m)>>1) : 16; \
+		REALLOC((a), (m)); \
+	} while (0)
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -64,6 +74,7 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 void mm_idxopt_init(mm_idxopt_t *opt);
 const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f);
+int mm_idx_getseq2(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq, int8_t *b);
 mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
 mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
 
@@ -72,12 +83,12 @@ 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_squeeze_a(void *km, int n_regs, mm_reg1_t *regs, mm128_t *a);
 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, int sub_diff);
+void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff, int hard_mask_level);
 void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_, mm_reg1_t *r);
 void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int max_gap_ref, int min_diff, int best_n, int n_segs, const int *qlens, int *n_, mm_reg1_t *r);
-void mm_filter_regs(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs);
+void mm_filter_regs(const mm_mapopt_t *opt, int qlen, 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_hit_sort(void *km, int *n_regs, mm_reg1_t *r);
 void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr);
 
 void mm_est_err(const mm_idx_t *mi, int qlen, int n_regs, mm_reg1_t *regs, const mm128_t *a, int32_t n, const uint64_t *mini_pos);
@@ -86,7 +97,14 @@ mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int
 void mm_seg_free(void *km, int n_segs, mm_seg_t *segs);
 void mm_pair(void *km, int max_gap_ref, int dp_bonus, int sub_diff, int match_sc, const int *qlens, int *n_regs, mm_reg1_t **regs);
 
+FILE *mm_split_init(const char *prefix, const mm_idx_t *mi);
+mm_idx_t *mm_split_merge_prep(const char *prefix, int n_splits, FILE **fp, uint32_t *n_seq_part);
+int mm_split_merge(int n_segs, const char **fn, const mm_mapopt_t *opt, int n_split_idx);
+void mm_split_rm_tmp(const char *prefix, int n_splits);
+
 void mm_err_puts(const char *str);
+void mm_err_fwrite(const void *p, size_t size, size_t nitems, FILE *fp);
+void mm_err_fread(void *p, size_t size, size_t nitems, FILE *fp);
 
 #ifdef __cplusplus
 }

options.c

@@ -31,8 +31,10 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->max_join_long = 20000;
 	opt->max_join_short = 2000;
 	opt->min_join_flank_sc = 1000;
+	opt->min_join_flank_ratio = 0.5f;
 
 	opt->a = 2, opt->b = 4, opt->q = 4, opt->e = 2, opt->q2 = 24, opt->e2 = 1;
+	opt->sc_ambi = 1;
 	opt->zdrop = 400, opt->zdrop_inv = 200;
 	opt->end_bonus = -1;
 	opt->min_dp_max = opt->min_chain_score * opt->a;
@@ -47,7 +49,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 
 void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)
 {
-	if ((opt->flag & MM_F_SPLICE_FOR) && (opt->flag & MM_F_SPLICE_REV))
+	if ((opt->flag & MM_F_SPLICE_FOR) || (opt->flag & MM_F_SPLICE_REV))
 		opt->flag |= MM_F_SPLICE;
 	if (opt->mid_occ <= 0)
 		opt->mid_occ = mm_idx_cal_max_occ(mi, opt->mid_occ_frac);
@@ -72,11 +74,11 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		io->flag = 0, io->k = 15, io->w = 5;
 		mo->flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN;
 		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
+		mo->bw = 2000;
 	} else if (strcmp(preset, "ava-pb") == 0) {
 		io->flag |= MM_I_HPC, io->k = 19, io->w = 5;
 		mo->flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN;
 		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
-		mo->bw = 2000;
 	} else if (strcmp(preset, "map10k") == 0 || strcmp(preset, "map-pb") == 0) {
 		io->flag |= MM_I_HPC, io->k = 19;
 	} else if (strcmp(preset, "map-ont") == 0) {
@@ -130,6 +132,16 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 
 int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
 {
+	if (mo->split_prefix && (mo->flag & (MM_F_OUT_CS|MM_F_OUT_MD))) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m --cs or --MD doesn't work with --split-prefix\033[0m\n");
+		return -6;
+	}
+	if (io->k <= 0 || io->w <= 0) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m -k and -w must be positive\033[0m\n");
+		return -5;
+	}
 	if (mo->best_n < 0) {
 		if (mm_verbose >= 1)
 			fprintf(stderr, "[ERROR]\033[1;31m -N must be no less than 0\033[0m\n");
@@ -162,5 +174,10 @@ int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
 			fprintf(stderr, "[ERROR]\033[1;31m Z-drop should not be less than inversion-Z-drop\033[0m\n");
 		return -5;
 	}
+	if ((mo->flag & MM_F_NO_PRINT_2ND) && (mo->flag & MM_F_ALL_CHAINS)) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m -X/-P and --secondary=no can't be applied at the same time\033[0m\n");
+		return -5;
+	}
 	return 0;
 }

pe.c

@@ -105,7 +105,7 @@ void mm_pair(void *km, int max_gap_ref, int pe_bonus, int sub_diff, int match_sc
 	max = -1;
 	max_idx[0] = max_idx[1] = -1;
 	last[0] = last[1] = -1;
-	kv_resize(uint64_t, km, sc, n);
+	kv_resize(uint64_t, km, sc, (size_t)n);
 	for (i = 0; i < n; ++i) {
 		if (a[i].key & 1) { // reverse first read or forward second read
 			mm_reg1_t *q, *r;
@@ -151,8 +151,8 @@ void mm_pair(void *km, int max_gap_ref, int pe_bonus, int sub_diff, int match_sc
 			}
 		}
 		mapq_pe = r[0]->mapq > r[1]->mapq? r[0]->mapq : r[1]->mapq;
-		for (i = 0; i < sc.n; ++i)
-			if ((sc.a[i]>>32) + sub_diff >= max>>32)
+		for (i = 0; i < (int)sc.n; ++i)
+			if ((sc.a[i]>>32) + sub_diff >= (uint64_t)max>>32)
 				++n_sub;
 		if (sc.n > 1) {
 			int mapq_pe_alt;
@@ -164,7 +164,7 @@ void mm_pair(void *km, int max_gap_ref, int pe_bonus, int sub_diff, int match_sc
 		if (sc.n == 1) {
 			if (r[0]->mapq < 2) r[0]->mapq = 2;
 			if (r[1]->mapq < 2) r[1]->mapq = 2;
-		} else if (max>>32 > sc.a[sc.n - 2]>>32) {
+		} else if ((uint64_t)max>>32 > sc.a[sc.n - 2]>>32) {
 			if (r[0]->mapq < 1) r[0]->mapq = 1;
 			if (r[1]->mapq < 1) r[1]->mapq = 1;
 		}

python/README.rst

@@ -43,20 +43,24 @@ The following Python script demonstrates the key functionality of mappy:
 APIs
 ----
 
-Mappy implements two classes and one global function.
+Mappy implements two classes and two global function.
 
 Class mappy.Aligner
 ~~~~~~~~~~~~~~~~~~~
 
 .. code:: python
 
-	mappy.Aligner(fn_idx_in, preset=None, ...)
+	mappy.Aligner(fn_idx_in=None, 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.
+  sequence file can be optionally gzip'd. This option has no effect if **seq**
+  is set.
+
+* **seq**: a single sequence to index. The sequence name will be set to
+  :code:`N/A`.
 
 * **preset**: minimap2 preset. Currently, minimap2 supports the following
   presets: **sr** for single-end short reads; **map-pb** for PacBio
@@ -79,17 +83,28 @@ This constructor accepts the following arguments:
 
 * **n_threads**: number of indexing threads; 3 by default
 
-* **fn_idx_out**: name of file to which the index is written
+* **extra_flags**: additional flags defined in minimap.h
+
+* **fn_idx_out**: name of file to which the index is written. This parameter
+  has no effect if **seq** is set.
+
+* **scoring**: scoring system. It is a tuple/list consisting of 4, 6 or 7
+  positive integers. The first 4 elements specify match scoring, mismatch
+  penalty, gap open and gap extension penalty. The 5th and 6th elements, if
+  present, set long-gap open and long-gap extension penalty. The 7th sets a
+  mismatch penalty involving ambiguous bases.
 
 .. code:: python
 
-	mappy.Aligner.map(seq, seq2=None)
+	mappy.Aligner.map(seq, seq2=None, cs=False, MD=False)
 
 This method aligns :code:`seq` against the index. It is a generator, *yielding*
 a series of :code:`mappy.Alignment` objects. If :code:`seq2` is present, mappy
 performs paired-end alignment, assuming the two ends are in the FR orientation.
 Alignments of the two ends can be distinguished by the :code:`read_num` field
-(see Class mappy.Alignment below).
+(see Class mappy.Alignment below). Argument :code:`cs` asks mappy to generate
+the :code:`cs` tag; :code:`MD` is similar. These two arguments might slightly
+degrade performance and are not enabled by default.
 
 .. code:: python
 
@@ -139,6 +154,11 @@ properties:
 * **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.
 
+* **MD**: the :code:`MD` tag as in the SAM format. It is an empty string unless
+  the :code:`MD` argument is applied when calling :code:`mappy.Aligner.map()`.
+
+* **cs**: the :code:`cs` tag.
+
 An :code:`Alignment` object can be converted to a string with :code:`str()` in
 the following format:
 

python/cmappy.h

@@ -20,6 +20,11 @@ typedef struct {
 	uint32_t *cigar32;
 } mm_hitpy_t;
 
+typedef struct {
+	int32_t n, m;
+	mm_idx_bed_t *r;
+} mm_bedpy_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;
@@ -73,13 +78,17 @@ static inline void mm_reset_timer(void)
 extern unsigned char seq_comp_table[256];
 static inline mm_reg1_t *mm_map_aux(const mm_idx_t *mi, const char *seq1, const char *seq2, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt)
 {
+	mm_reg1_t *r;
+
+	Py_BEGIN_ALLOW_THREADS
 	if (seq2 == 0) {
-		return mm_map(mi, strlen(seq1), seq1, n_regs, b, opt, NULL);
+		r = mm_map(mi, strlen(seq1), seq1, n_regs, b, opt, NULL);
 	} else {
 		int _n_regs[2];
 		mm_reg1_t *regs[2];
 		char *seq[2];
 		int i, len[2];
+
 		len[0] = strlen(seq1);
 		len[1] = strlen(seq2);
 		seq[0] = (char*)seq1;
@@ -97,8 +106,11 @@ static inline mm_reg1_t *mm_map_aux(const mm_idx_t *mi, const char *seq1, const
 		regs[0] = (mm_reg1_t*)realloc(regs[0], sizeof(mm_reg1_t) * (*n_regs));
 		memcpy(&regs[0][_n_regs[0]], regs[1], _n_regs[1] * sizeof(mm_reg1_t));
 		free(regs[1]);
-		return regs[0];
+		r = regs[0];
 	}
+	Py_END_ALLOW_THREADS
+
+	return r;
 }
 
 static inline char *mappy_revcomp(int len, const uint8_t *seq)
@@ -119,8 +131,8 @@ static char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int e
 	*len = 0;
 	rid = mm_idx_name2id(mi, name);
 	if (rid < 0) return 0;
-	if (st >= mi->seq[rid].len || st >= en) return 0;
-	if (en < 0 || en > mi->seq[rid].len)
+	if ((uint32_t)st >= mi->seq[rid].len || st >= en) return 0;
+	if (en < 0 || (uint32_t)en > mi->seq[rid].len)
 		en = mi->seq[rid].len;
 	s = (char*)malloc(en - st + 1);
 	*len = mm_idx_getseq(mi, rid, st, en, (uint8_t*)s);
@@ -130,4 +142,45 @@ static char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int e
 	return s;
 }
 
+static mm_idx_t *mappy_idx_seq(int w, int k, int is_hpc, int bucket_bits, const char *seq, int len)
+{
+	const char *fake_name = "N/A";
+	char *s;
+	mm_idx_t *mi;
+	s = (char*)calloc(len + 1, 1);
+	memcpy(s, seq, len);
+	mi = mm_idx_str(w, k, is_hpc, bucket_bits, 1, (const char**)&s, (const char**)&fake_name);
+	free(s);
+	return mi;
+}
+
+static mm_bedpy_t *mappy_bed_new(void)
+{
+	return (mm_bedpy_t*)calloc(1, sizeof(mm_bedpy_t));
+}
+
+static int mappy_bed_add(mm_bedpy_t *bed, mm_idx_t *mi, const char *name, uint32_t st, uint32_t en)
+{
+	mm_idx_bed_t *b;
+	int id;
+	if (mi->h == 0) mm_idx_index_name(mi);
+	if (bed->n == bed->m) {
+		bed->m = bed->m? bed->m + (bed->m>>1) : 16;
+		bed->r = (mm_idx_bed_t*)realloc(bed->r, sizeof(mm_idx_bed_t) * bed->m);
+	}
+	id = mm_idx_name2id(mi, name);
+	if (id < 0 || st >= en) return -1;
+	if (en > mi->seq[id].len) en = mi->seq[id].len;
+	b = &bed->r[bed->n++];
+	b->x = (uint64_t)id << 32 | st;
+	b->end = en, b->idx = -1;
+	return 0;
+}
+
+static void mappy_bed_finalize(mm_bedpy_t *bed, mm_idx_t *mi)
+{
+	mm_idx_bed_attach(mi, bed->n, bed->r); // bed->r is now owned by mi and will be deallocated with it
+	free(bed);
+}
+
 #endif

python/cmappy.pxd

@@ -24,7 +24,9 @@ cdef extern from "minimap.h":
 		int best_n
 		int max_join_long, max_join_short
 		int min_join_flank_sc
+		float min_join_flank_ratio;
 		int a, b, q, e, q2, e2
+		int sc_ambi
 		int noncan
 		int zdrop, zdrop_inv
 		int end_bonus
@@ -38,6 +40,7 @@ cdef extern from "minimap.h":
 		int32_t mid_occ
 		int32_t max_occ
 		int mini_batch_size
+		const char *split_prefix
 
 	int mm_set_opt(char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 	int mm_verbose
@@ -50,6 +53,10 @@ cdef extern from "minimap.h":
 		uint64_t offset
 		uint32_t len
 
+	ctypedef struct mm_idx_bed_t:
+		uint64_t x
+		int32_t end, idx
+
 	ctypedef struct mm_idx_bucket_t:
 		pass
 
@@ -59,6 +66,7 @@ cdef extern from "minimap.h":
 		mm_idx_seq_t *seq
 		uint32_t *S
 		mm_idx_bucket_t *B
+		mm_idx_bed_t *R
 		void *km
 		void *h
 
@@ -84,6 +92,9 @@ cdef extern from "minimap.h":
 
 	mm_tbuf_t *mm_tbuf_init()
 	void mm_tbuf_destroy(mm_tbuf_t *b)
+	void *mm_tbuf_get_km(mm_tbuf_t *b)
+	int mm_gen_cs(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int no_iden)
+	int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq)
 
 #
 # Helper header (because it is hard to expose mm_reg1_t with Cython)
@@ -104,6 +115,15 @@ cdef extern from "cmappy.h":
 	void mm_free_reg1(mm_reg1_t *r)
 	mm_reg1_t *mm_map_aux(const mm_idx_t *mi, const char *seq1, const char *seq2, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt)
 	char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int en, int *l)
+	mm_idx_t *mappy_idx_seq(int w, int k, int is_hpc, int bucket_bits, const char *seq, int l)
+
+	ctypedef struct mm_bedpy_t:
+		int32_t n, m
+		mm_idx_bed_t *r
+
+	mm_bedpy_t *mappy_bed_new()
+	int mappy_bed_add(mm_bedpy_t *bed, mm_idx_t *mi, const char *name, uint32_t st, uint32_t en)
+	void mappy_bed_finalize(mm_bedpy_t *bed, mm_idx_t *mi)
 
 	ctypedef struct kstring_t:
 		unsigned l, m

python/mappy.pyx

@@ -3,6 +3,8 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
+__version__ = '2.13'
+
 cmappy.mm_reset_timer()
 
 cdef class Alignment:
@@ -12,9 +14,9 @@ cdef class Alignment:
 	cdef int8_t _strand, _trans_strand
 	cdef uint8_t _mapq, _is_primary
 	cdef int _seg_id
-	cdef _ctg, _cigar # these are python objects
+	cdef _ctg, _cigar, _cs, _MD # these are python objects
 
-	def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand, seg_id):
+	def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand, seg_id, cs_str, MD_str):
 		self._ctg = ctg if isinstance(ctg, str) else ctg.decode()
 		self._ctg_len, self._r_st, self._r_en = cl, cs, ce
 		self._strand, self._q_st, self._q_en = strand, qs, qe
@@ -24,6 +26,8 @@ cdef class Alignment:
 		self._is_primary = is_primary
 		self._trans_strand = trans_strand
 		self._seg_id = seg_id
+		self._cs = cs_str
+		self._MD = MD_str
 
 	@property
 	def ctg(self): return self._ctg
@@ -70,6 +74,12 @@ cdef class Alignment:
 	@property
 	def read_num(self): return self._seg_id + 1
 
+	@property
+	def cs(self): return self._cs
+
+	@property
+	def MD(self): return self._MD
+
 	@property
 	def cigar_str(self):
 		return "".join(map(lambda x: str(x[0]) + 'MIDNSH'[x[1]], self._cigar))
@@ -83,8 +93,10 @@ cdef class Alignment:
 		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._mlen), str(self._blen), str(self._mapq), tp, ts, "cg:Z:" + self.cigar_str])
+		a = [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._mlen), str(self._blen), str(self._mapq), tp, ts, "cg:Z:" + self.cigar_str]
+		if self._cs != "": a.append("cs:Z:" + self._cs)
+		return "\t".join(a)
 
 cdef class ThreadBuffer:
 	cdef cmappy.mm_tbuf_t *_b
@@ -100,7 +112,7 @@ cdef class Aligner:
 	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):
+	def __cinit__(self, fn_idx_in=None, 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, max_frag_len=None, extra_flags=None, seq=None, scoring=None, bed=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
@@ -113,17 +125,42 @@ cdef class Aligner:
 		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.map_opt.best_n = best_n
+		if max_frag_len is not None: self.map_opt.max_frag_len = max_frag_len
+		if extra_flags is not None: self.map_opt.flag |= extra_flags
+		if scoring is not None and len(scoring) >= 4:
+			self.map_opt.a, self.map_opt.b = scoring[0], scoring[1]
+			self.map_opt.q, self.map_opt.e = scoring[2], scoring[3]
+			self.map_opt.q2, self.map_opt.e2 = self.map_opt.q, self.map_opt.e
+			if len(scoring) >= 6:
+				self.map_opt.q2, self.map_opt.e2 = scoring[4], scoring[5]
+				if len(scoring) >= 7:
+					self.map_opt.sc_ambi = scoring[6]
 
 		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)
+		cdef cmappy.mm_bedpy_t *bed_agg;
+
+		if seq is None:
+			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)
+				cmappy.mm_idx_index_name(self._idx)
 		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)
+			self._idx = cmappy.mappy_idx_seq(self.idx_opt.w, self.idx_opt.k, self.idx_opt.flag&1, self.idx_opt.bucket_bits, str.encode(seq), len(seq))
 			cmappy.mm_mapopt_update(&self.map_opt, self._idx)
-			cmappy.mm_idx_index_name(self._idx)
+			self.map_opt.mid_occ = 1000 # don't filter high-occ seeds
+
+		if bed is not None:
+			bed_agg = cmappy.mappy_bed_new()
+			for b in bed:
+				if len(b) < 3: en = int(b[1]) + 1
+				else: en = int(b[2])
+				cmappy.mappy_bed_add(bed_agg, self._idx, str.encode(b[0]), int(b[1]), en)
+			cmappy.mappy_bed_finalize(bed_agg, self._idx)
 
 	def __dealloc__(self):
 		if self._idx is not NULL:
@@ -132,32 +169,49 @@ cdef class Aligner:
 	def __bool__(self):
 		return (self._idx != NULL)
 
-	def map(self, seq, seq2=None, buf=None):
+	def map(self, seq, seq2=None, buf=None, cs=False, MD=False, max_frag_len=None, extra_flags=None):
 		cdef cmappy.mm_reg1_t *regs
 		cdef cmappy.mm_hitpy_t h
 		cdef ThreadBuffer b
 		cdef int n_regs
+		cdef char *cs_str = NULL
+		cdef int l_cs_str, m_cs_str = 0
+		cdef void *km
+		cdef cmappy.mm_mapopt_t map_opt
+
+		map_opt = self.map_opt
+		if max_frag_len is not None: map_opt.max_frag_len = max_frag_len
+		if extra_flags is not None: map_opt.flag |= extra_flags
 
 		if self._idx is NULL: return None
 		if buf is None: b = ThreadBuffer()
 		else: b = buf
+		km = cmappy.mm_tbuf_get_km(b._b)
 
 		_seq = seq if isinstance(seq, bytes) else seq.encode()
 		if seq2 is None:
-			regs = cmappy.mm_map_aux(self._idx, _seq, NULL,  &n_regs, b._b, &self.map_opt)
+			regs = cmappy.mm_map_aux(self._idx, _seq, NULL,  &n_regs, b._b, &map_opt)
 		else:
 			_seq2 = seq2 if isinstance(seq2, bytes) else seq2.encode()
-			regs = cmappy.mm_map_aux(self._idx, _seq, _seq2, &n_regs, b._b, &self.map_opt)
+			regs = cmappy.mm_map_aux(self._idx, _seq, _seq2, &n_regs, b._b, &map_opt)
 
 		for i in range(n_regs):
 			cmappy.mm_reg2hitpy(self._idx, &regs[i], &h)
-			cigar = []
-			for k in range(h.n_cigar32):
+			cigar, _cs, _MD = [], '', ''
+			for k in range(h.n_cigar32): # convert the 32-bit CIGAR encoding to Python array
 				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.mlen, h.blen, h.NM, h.trans_strand, h.seg_id)
+			if cs or MD: # generate the cs and/or the MD tag, if requested
+				if cs:
+					l_cs_str = cmappy.mm_gen_cs(km, &cs_str, &m_cs_str, self._idx, &regs[i], _seq, 1)
+					_cs = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
+				if MD:
+					l_cs_str = cmappy.mm_gen_MD(km, &cs_str, &m_cs_str, self._idx, &regs[i], _seq)
+					_MD = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
+			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.mlen, h.blen, h.NM, h.trans_strand, h.seg_id, _cs, _MD)
 			cmappy.mm_free_reg1(&regs[i])
 		free(regs)
+		free(cs_str)
 
 	def seq(self, str name, int start=0, int end=0x7fffffff):
 		cdef int l

python/minimap2.py

@@ -1,10 +1,11 @@
 #!/usr/bin/env python
 
-import sys, getopt
+import sys
+import getopt
 import mappy as mp
 
 def main(argv):
-	opts, args = getopt.getopt(argv[1:], "x:n:m:k:w:r:")
+	opts, args = getopt.getopt(argv[1:], "x:n:m:k:w:r:c")
 	if len(args) < 2:
 		print("Usage: minimap2.py [options] <ref.fa>|<ref.mmi> <query.fq>")
 		print("Options:")
@@ -14,9 +15,11 @@ def main(argv):
 		print("  -k INT      k-mer length")
 		print("  -w INT      minimizer window length")
 		print("  -r INT      band width")
+		print("  -c          output the cs tag")
 		sys.exit(1)
 
-	preset, min_cnt, min_sc, k, w, bw = None, None, None, None, None, None
+	preset = min_cnt = min_sc = k = w = bw = None
+	out_cs = False
 	for opt, arg in opts:
 		if opt == '-x': preset = arg
 		elif opt == '-n': min_cnt = int(arg)
@@ -24,11 +27,12 @@ def main(argv):
 		elif opt == '-r': bw = int(arg)
 		elif opt == '-k': k = int(arg)
 		elif opt == '-w': w = int(arg)
+		elif opt == '-c': out_cs = True
 
 	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
+		for h in a.map(seq, cs=out_cs): # traverse hits
 			print('{}\t{}\t{}'.format(name, len(seq), h))
 
 if __name__ == "__main__":

sdust.c

@@ -70,10 +70,10 @@ void sdust_buf_destroy(sdust_buf_t *buf)
 static inline void shift_window(int t, kdq_t(int) *w, int T, int W, int *L, int *rw, int *rv, int *cw, int *cv)
 {
 	int s;
-	if (kdq_size(w) >= W - SD_WLEN + 1) { // TODO: is this right for SD_WLEN!=3?
+	if ((int)kdq_size(w) >= W - SD_WLEN + 1) { // TODO: is this right for SD_WLEN!=3?
 		s = *kdq_shift(int, w);
 		*rw -= --cw[s];
-		if (*L > kdq_size(w))
+		if (*L > (int)kdq_size(w))
 			--*L, *rv -= --cv[s];
 	}
 	kdq_push(int, w, t);
@@ -114,7 +114,7 @@ static void find_perfect(void *km, perf_intv_v *P, const kdq_t(int) *w, int T, i
 		r += c[t]++;
 		new_r = r, new_l = kdq_size(w) - i - 1;
 		if (new_r * 10 > T * new_l) {
-			for (j = 0; j < P->n && P->a[j].start >= i + start; ++j) { // find insertion position
+			for (j = 0; j < (int)P->n && P->a[j].start >= i + start; ++j) { // find insertion position
 				perf_intv_t *p = &P->a[j];
 				if (max_r == 0 || p->r * max_l > max_r * p->l)
 					max_r = p->r, max_l = p->l;
@@ -177,7 +177,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 "getopt.h"
+#include "ketopt.h"
 #include "kseq.h"
 KSEQ_INIT(gzFile, gzread)
 
@@ -186,16 +186,17 @@ int main(int argc, char *argv[])
 	gzFile fp;
 	kseq_t *ks;
 	int W = 64, T = 20, c;
+	ketopt_t o = KETOPT_INIT;
 
-	while ((c = getopt(argc, argv, "w:t:")) >= 0) {
-		if (c == 'w') W = atoi(optarg);
-		else if (c == 't') T = atoi(optarg);
+	while ((c = ketopt(&o, argc, argv, 1, "w:t:", 0)) >= 0) {
+		if (c == 'w') W = atoi(o.arg);
+		else if (c == 't') T = atoi(o.arg);
 	}
-	if (optind == argc) {
+	if (o.ind == argc) {
 		fprintf(stderr, "Usage: sdust [-w %d] [-t %d] <in.fa>\n", W, T);
 		return 1;
 	}
-	fp = strcmp(argv[optind], "-")? gzopen(argv[optind], "r") : gzdopen(fileno(stdin), "r");
+	fp = strcmp(argv[o.ind], "-")? gzopen(argv[o.ind], "r") : gzdopen(fileno(stdin), "r");
 	ks = kseq_init(fp);
 	while (kseq_read(ks) >= 0) {
 		uint64_t *r;

setup.py

@@ -14,16 +14,26 @@ else: # with Cython
 	module_src = 'python/mappy.pyx'
 	cmdclass['build_ext'] = build_ext
 
-import sys
+import sys, platform
+
 sys.path.append('python')
 
+extra_compile_args = ['-DHAVE_KALLOC']
+include_dirs = ["."]
+
+if platform.machine() in ["aarch64", "arm64"]:
+	include_dirs.append("sse2neon/")
+	extra_compile_args.extend(['-ftree-vectorize', '-DKSW_SSE2_ONLY', '-D__SSE2__'])
+else:
+	extra_compile_args.append('-msse4.1') # WARNING: ancient x86_64 CPUs don't have SSE4
+
 def readme():
-    with open('python/README.rst') as f:
-        return f.read()
+	with open('python/README.rst') as f:
+		return f.read()
 
 setup(
 	name = 'mappy',
-	version = '2.10',
+	version = '2.13',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),
@@ -35,12 +45,12 @@ setup(
     ext_modules = [Extension('mappy',
 		sources = [module_src, 'align.c', 'bseq.c', 'chain.c', 'format.c', 'hit.c', 'index.c', 'pe.c', 'options.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', 'esterr.c'],
+				   'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c', 'esterr.c', 'splitidx.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 = ['-DHAVE_KALLOC', '-msse4.1'], # WARNING: ancient x86_64 CPUs don't have SSE4
-		include_dirs = ['.'],
+		extra_compile_args = extra_compile_args,
+		include_dirs = include_dirs,
 		libraries = ['z', 'm', 'pthread'])],
 	classifiers = [
 		'Development Status :: 5 - Production/Stable',

splitidx.c

@@ -0,0 +1,80 @@
+#include <string.h>
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "mmpriv.h"
+
+FILE *mm_split_init(const char *prefix, const mm_idx_t *mi)
+{
+	char *fn;
+	FILE *fp;
+	uint32_t i, k = mi->k;
+	fn = (char*)calloc(strlen(prefix) + 10, 1);
+	sprintf(fn, "%s.%.4d.tmp", prefix, mi->index);
+	fp = fopen(fn, "wb");
+	assert(fp);
+	mm_err_fwrite(&k, 4, 1, fp);
+	mm_err_fwrite(&mi->n_seq, 4, 1, fp);
+	for (i = 0; i < mi->n_seq; ++i) {
+		uint8_t l;
+		l = strlen(mi->seq[i].name);
+		mm_err_fwrite(&l, 1, 1, fp);
+		mm_err_fwrite(mi->seq[i].name, 1, l, fp);
+		mm_err_fwrite(&mi->seq[i].len, 4, 1, fp);
+	}
+	free(fn);
+	return fp;
+}
+
+mm_idx_t *mm_split_merge_prep(const char *prefix, int n_splits, FILE **fp, uint32_t *n_seq_part)
+{
+	mm_idx_t *mi = 0;
+	char *fn;
+	int i, j;
+
+	if (n_splits < 1) return 0;
+	fn = CALLOC(char, strlen(prefix) + 10);
+	for (i = 0; i < n_splits; ++i) {
+		sprintf(fn, "%s.%.4d.tmp", prefix, i);
+		if ((fp[i] = fopen(fn, "rb")) == 0) {
+			if (mm_verbose >= 1)
+				fprintf(stderr, "ERROR: failed to open temporary file '%s'\n", fn);
+			for (j = 0; j < i; ++j)
+				fclose(fp[j]);
+			free(fn);
+			return 0;
+		}
+	}
+	free(fn);
+
+	mi = CALLOC(mm_idx_t, 1);
+	for (i = 0; i < n_splits; ++i) {
+		mm_err_fread(&mi->k, 4, 1, fp[i]); // TODO: check if k is all the same
+		mm_err_fread(&n_seq_part[i], 4, 1, fp[i]);
+		mi->n_seq += n_seq_part[i];
+	}
+	mi->seq = CALLOC(mm_idx_seq_t, mi->n_seq);
+	for (i = j = 0; i < n_splits; ++i) {
+		uint32_t k;
+		for (k = 0; k < n_seq_part[i]; ++k, ++j) {
+			uint8_t l;
+			mm_err_fread(&l, 1, 1, fp[i]);
+			mi->seq[j].name = (char*)calloc(l + 1, 1);
+			mm_err_fread(mi->seq[j].name, 1, l, fp[i]);
+			mm_err_fread(&mi->seq[j].len, 4, 1, fp[i]);
+		}
+	}
+	return mi;
+}
+
+void mm_split_rm_tmp(const char *prefix, int n_splits)
+{
+	int i;
+	char *fn;
+	fn = CALLOC(char, strlen(prefix) + 10);
+	for (i = 0; i < n_splits; ++i) {
+		sprintf(fn, "%s.%.4d.tmp", prefix, i);
+		remove(fn);
+	}
+	free(fn);
+}