Release minimap2-2.15 (r905)

This was caused by mm_fix_cigar() which may change query/target offset in very
rare cases. Generating EQX has to beware of this change.

Resolves #266

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
+else				# if arm_neon is defined
-    OBJS+=ksw2_extz2_neon.o ksw2_extd2_neon.o ksw2_exts2_neon.o
+	OBJS+=ksw2_extz2_neon.o ksw2_extd2_neon.o ksw2_exts2_neon.o
-    CFLAGS+=-D_FILE_OFFSET_BITS=64 -mfpu=neon -fsigned-char
     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
+minimap2:main.o libminimap2.a
-		$(CC) $(CFLAGS) main.o getopt.o -o $@ -L. -lminimap2 $(LIBS)
+		$(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
+sdust:sdust.c kalloc.o kalloc.h kdq.h kvec.h kseq.h ketopt.h sdust.h
-		$(CC) -D_SDUST_MAIN $(CFLAGS) $< getopt.o kalloc.o -o $@ -lz
+		$(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,167 @@
+Release 2.15-r905 (10 January 2019)
+-----------------------------------
+
+Changes to minimap2:
+
+ * Fixed a rare segmentation fault when option -H is in use (#307). This may
+   happen when there are very long homopolymers towards the 5'-end of a read.
+
+ * Fixed wrong CIGARs when option --eqx is used (#266).
+
+ * Fixed a typo in the base encoding table (#264). This should have no
+   practical effect.
+
+ * Fixed a typo in the example code (#265).
+
+ * Improved the C++ compatibility by removing "register" (#261). However,
+   minimap2 still can't be compiled in the pedantic C++ mode (#306).
+
+ * Output a new "de" tag for gap-compressed sequence divergence.
+
+Changes to paftools.js:
+
+ * Added "asmgene" to evaluate the completeness of an assembly by measuring the
+   uniquely mapped single-copy genes. This command learns the idea of BUSCO.
+
+ * Added "vcfpair" to call a phased VCF from phased whole-genome assemblies. An
+   earlier version of this script is used to produce the ground truth for the
+   syndip benchmark [PMID:30013044].
+
+This release produces identical alignment coordinates and CIGARs in comparison
+to v2.14. Users are advised to upgrade due to the several bug fixes.
+
+(2.15: 10 Janurary 2019, r905)
+
+
+
+Release 2.14-r883 (5 November 2018)
+-----------------------------------
+
+Notable changes:
+
+ * Fixed two minor bugs caused by typos (#254 and #266).
+
+ * Fixed a bug that made minimap2 abort when --eqx was used together with --MD
+   or --cs (#257).
+
+ * Added --cap-sw-mem to cap the size of DP matrices (#259). Base alignment may
+   take a lot of memory in the splicing mode. This may lead to issues when we
+   run minimap2 on a cluster with a hard memory limit. The new option avoids
+   unlimited memory usage at the cost of missing a few long introns.
+
+ * Conforming to C99 and C11 when possible (#261).
+
+ * Warn about malformatted FASTA or FASTQ (#252 and #255).
+
+This release occasionally produces base alignments different from v2.13. The
+overall alignment accuracy remain similar.
+
+(2.14: 5 November 2018, r883)
+
+
+
+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 ref.fa rna-reads.fa > aln.sam       # spliced long reads (strand unknown)
-./minimap2 -ax splice -k14 -uf ref.fa reads.fa > aln.sam  # Nanopore Direct RNA-seq
+./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 -
+curl -L https://github.com/lh3/minimap2/releases/download/v2.15/minimap2-2.15_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.10_x64-linux/minimap2
+./minimap2-2.15_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
 
@@ -349,6 +355,8 @@ mappy` or [from BioConda][mappyconda] via `conda install -c bioconda mappy`.
   billion bases or longer (2,147,483,647 to be exact). The total length of all
   sequences can well exceed this threshold.
 
+* Minimap2 often misses small exons.
+
 
 
 [paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
@@ -365,3 +373,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)
@@ -90,7 +91,8 @@ static int mm_test_zdrop(void *km, const mm_mapopt_t *opt, const uint8_t *qseq,
 static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int *qshift, int *tshift)
 {
 	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;
 	for (k = 0; k < p->n_cigar; ++k) { // indel left alignment
@@ -145,10 +147,81 @@ 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_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq) // written by @armintoepfer
 {
-	uint32_t k, l, toff = 0, qoff = 0;
+	uint32_t n_EQX = 0;
-	int32_t s = 0, max = 0, qshift, tshift;
+	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_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 is_eqx)
+{
+	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);
@@ -190,6 +263,7 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
 	}
 	p->dp_max = max;
 	assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
+	if (is_eqx) mm_update_cigar_eqx(r, qseq, tseq); // NB: it has to be called here as changes to qseq and tseq are not returned
 }
 
 static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) // TODO: this calls the libc realloc()
@@ -197,12 +271,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) {
+	} 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);
+		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);
 	}
@@ -227,7 +301,10 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 		for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr);
 		fputc('\n', stderr);
 	}
-	if (opt->flag & MM_F_SPLICE)
+	if (opt->max_sw_mat > 0 && (int64_t)tlen * qlen > opt->max_sw_mat) {
+		ksw_reset_extz(ez);
+		ez->zdropped = 1;
+	} else if (opt->flag & MM_F_SPLICE)
 		ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, flag, ez);
 	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);
@@ -268,20 +345,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 +380,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 +388,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 +401,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].y + 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 +473,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 +511,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);
@@ -434,22 +557,24 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 
 	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) {
+		if (!(opt->flag & MM_F_NO_END_FLT)) {
-			mm_fix_bad_ends_splice(km, opt, mi, r, mat, qlen, qseq0, a, &as1, &cnt1);
+			if (is_splice)
-		} else {
+				mm_fix_bad_ends_splice(km, opt, mi, r, mat, qlen, qseq0, a, &as1, &cnt1);
-			mm_fix_bad_ends(r, a, opt->bw, opt->min_chain_score * 2, &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);
 	}
@@ -473,7 +598,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 		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;
+		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);
@@ -488,6 +613,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				if (++l > opt->min_cnt) {
 					l = rs0 - x > qs0 - y? rs0 - x : qs0 - y;
 					rs1 = rs0 - l, qs1 = qs0 - l;
+					if (rs1 < 0) rs1 = 0; // not strictly necessary; better have this guard for explicit
 					break;
 				}
 			}
@@ -501,6 +627,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			l = l < rs? l : rs;
 			rs1 = rs1 > rs - l? rs1 : rs - l;
 			rs0 = rs0 < rs1? rs0 : rs1;
+			rs0 = rs0 < rs? rs0 : rs;
 		} else rs0 = rs, qs0 = qs;
 		// compute re0 and qe0
 		re0 = (int32_t)a[r->as + r->cnt - 1].x + 1;
@@ -517,13 +644,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;
@@ -540,7 +667,7 @@ 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);
 
-	if (qs > 0 && rs > 0) { // left extension
+	if (qs > 0 && rs > 0) { // left extension; probably the condition can be changed to "qs > qs0 && rs > rs0"
 		qseq = &qseq0[rev][qs0];
 		mm_idx_getseq(mi, rid, rs0, rs, tseq);
 		mm_seq_rev(qs - qs0, qseq);
@@ -627,7 +754,7 @@ 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);
+		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e, opt->flag & MM_F_EQX);
 		if (rev && r->p->trans_strand)
 			r->p->trans_strand ^= 3; // flip to the read strand
 	}
@@ -652,7 +779,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];
@@ -685,7 +812,7 @@ 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, opt->flag & MM_F_EQX);
 	ret = 1;
 end_align1_inv:
 	kfree(km, tseq);
@@ -755,7 +882,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_filter_regs(opt, qlen, n_regs_, regs);
-	mm_hit_sort_by_dp(km, n_regs_, regs);
+	mm_hit_sort(km, n_regs_, regs);
 	return regs;
 }

bseq.c

@@ -15,7 +15,7 @@ unsigned char seq_comp_table[256] = {
 	 48,  49,  50,	51,	 52,  53,  54,	55,	 56,  57,  58,	59,	 60,  61,  62,	63,
 	 64, 'T', 'V', 'G', 'H', 'E', 'F', 'C', 'D', 'I', 'J', 'M', 'L', 'K', 'N', 'O',
 	'P', 'Q', 'Y', 'S', 'A', 'A', 'B', 'W', 'X', 'R', 'Z',	91,	 92,  93,  94,	95,
-	 64, 't', 'v', 'g', 'h', 'e', 'f', 'c', 'd', 'i', 'j', 'm', 'l', 'k', 'n', 'o',
+	 96, 't', 'v', 'g', 'h', 'e', 'f', 'c', 'd', 'i', 'j', 'm', 'l', 'k', 'n', 'o',
 	'p', 'q', 'y', 's', 'a', 'a', 'b', 'w', 'x', 'r', 'z', 123, 124, 125, 126, 127,
 	128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
 	144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159,
@@ -39,7 +39,7 @@ mm_bseq_file_t *mm_bseq_open(const char *fn)
 {
 	mm_bseq_file_t *fp;
 	gzFile f;
-	f = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
+	f = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(0, "r");
 	if (f == 0) return 0;
 	fp = (mm_bseq_file_t*)calloc(1, sizeof(mm_bseq_file_t));
 	fp->fp = f;
@@ -65,9 +65,11 @@ static inline char *kstrdup(const kstring_t *s)
 static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual, int with_comment)
 {
 	int i;
+	if (ks->name.l == 0)
+		fprintf(stderr, "[WARNING]\033[1;31m empty sequence name in the input.\033[0m\n");
 	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;
@@ -78,6 +80,7 @@ static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual, int with_
 mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int with_comment, int frag_mode, int *n_)
 {
 	int64_t size = 0;
+	int ret;
 	kvec_t(mm_bseq1_t) a = {0,0,0};
 	kseq_t *ks = fp->ks;
 	*n_ = 0;
@@ -87,7 +90,7 @@ mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
 		size = fp->s.l_seq;
 		memset(&fp->s, 0, sizeof(mm_bseq1_t));
 	}
-	while (kseq_read(ks) >= 0) {
+	while ((ret = kseq_read(ks)) >= 0) {
 		mm_bseq1_t *s;
 		assert(ks->seq.l <= INT32_MAX);
 		if (a.m == 0) kv_resize(mm_bseq1_t, 0, a, 256);
@@ -107,6 +110,8 @@ mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
 			break;
 		}
 	}
+	if (ret < -1)
+		fprintf(stderr, "[WARNING]\033[1;31m wrong FASTA/FASTQ record. Continue anyway.\033[0m\n");
 	*n_ = a.n;
 	return a.a;
 }

chain.c

@@ -14,7 +14,7 @@ static const char LogTable256[256] = {
 
 static inline int ilog2_32(uint32_t v)
 {
-	register uint32_t t, tt;
+	uint32_t t, tt;
 	if ((tt = v>>16)) return (t = tt>>8) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
 	return (t = v>>8) ? 8 + LogTable256[t] : LogTable256[v];
 }
@@ -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 -
+curl -L https://github.com/lh3/minimap2/releases/download/v2.15/minimap2-2.15_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.10_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
+cp minimap2-2.15_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 -

example.c

@@ -45,7 +45,7 @@ int main(int argc, char *argv[])
 				printf("%s\t%d\t%d\t%d\t%c\t", ks->name.s, ks->seq.l, r->qs, r->qe, "+-"[r->rev]);
 				printf("%s\t%d\t%d\t%d\t%d\t%d\t%d\tcg:Z:", mi->seq[r->rid].name, mi->seq[r->rid].len, r->rs, r->re, r->mlen, r->blen, r->mapq);
 				for (i = 0; i < r->p->n_cigar; ++i) // IMPORTANT: this gives the CIGAR in the aligned regions. NO soft/hard clippings!
-					printf("%d%c", r->p->cigar[i]>>4, "MIDSHN"[r->p->cigar[i]&0xf]);
+					printf("%d%c", r->p->cigar[i]>>4, "MIDNSH"[r->p->cigar[i]&0xf]);
 				putchar('\n');
 				free(r->p);
 			}

format.c

@@ -92,7 +92,8 @@ static void sam_write_rg_line(kstring_t *str, const char *s)
 		if (mm_verbose >= 1) fprintf(stderr, "[ERROR] the read group line contained literal <tab> characters -- replace with escaped tabs: \\t\n");
 		goto err_set_rg;
 	}
-	rg_line = strdup(s);
+	rg_line = (char*)malloc(strlen(s) + 1);
+	strcpy(rg_line, s);
 	mm_escape(rg_line);
 	if ((p = strstr(rg_line, "\tID:")) == 0) {
 		if (mm_verbose >= 1) fprintf(stderr, "[ERROR] no ID within the read group line\n");
@@ -133,14 +134,14 @@ 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:");
+	if (write_tag) mm_sprintf_lite(s, "\tcs:Z:");
-	for (i = q_off = t_off = 0; i < r->p->n_cigar; ++i) {
+	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);
+		assert((op >= 0 && op <= 3) || op == 7 || op == 8);
-		if (op == 0) { // match
+		if (op == 0 || op == 7 || op == 8) { // match
 			int l_tmp = 0;
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
@@ -181,14 +182,14 @@ 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:");
+	if (write_tag) mm_sprintf_lite(s, "\tMD:Z:");
-	for (i = q_off = t_off = 0; i < r->p->n_cigar; ++i) {
+	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) || op == 7 || op == 8);
-		if (op == 0) { // match
+		if (op == 0 || op == 7 || op == 8) { // match
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
 					mm_sprintf_lite(s, "%d%c", l_MD, "ACGTN"[tseq[t_off + j]]);
@@ -204,13 +205,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 +233,46 @@ 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);
+	if (is_MD) write_MD_core(s, tseq, qseq, r, tmp, write_tag);
-	else write_cs_core(s, tseq, qseq, r, tmp, no_iden);
+	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);
+}
+
+double mm_event_identity(const mm_reg1_t *r)
+{
+	int32_t i, n_gapo = 0, n_gap = 0;
+	if (r->p == 0) return -1.0f;
+	for (i = 0; i < r->p->n_cigar; ++i) {
+		int32_t op = r->p->cigar[i] & 0xf, len = r->p->cigar[i] >> 4;
+		if (op == 1 || op == 2)
+			++n_gapo, n_gap += len;
+	}
+	return (double)r->mlen / (r->blen - r->p->n_ambi - n_gap + n_gapo);
+}
+
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 {
 	int type;
@@ -247,10 +285,17 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 	}
 	mm_sprintf_lite(s, "\ttp:A:%c\tcm:i:%d\ts1:i:%d", type, r->cnt, r->score);
 	if (r->parent == r->id) mm_sprintf_lite(s, "\ts2:i:%d", r->subsc);
-	if (r->div >= 0.0f && r->div <= 1.0f) {
+	if (r->p) {
-		char buf[8];
+		char buf[16];
+		double div;
+		div = 1.0 - mm_event_identity(r);
+		if (div == 0.0) buf[0] = '0', buf[1] = 0;
+		else snprintf(buf, 16, "%.4f", 1.0 - mm_event_identity(r));
+		mm_sprintf_lite(s, "\tde:f:%s", buf);
+	} else if (r->div >= 0.0f && r->div <= 1.0f) {
+		char buf[16];
 		if (r->div == 0.0f) buf[0] = '0', buf[1] = 0;
-		else sprintf(buf, "%.4f", r->div);
+		else snprintf(buf, 16, "%.4f", r->div);
 		mm_sprintf_lite(s, "\tdv:f:%s", buf);
 	}
 	if (r->split) mm_sprintf_lite(s, "\tzd:i:%d", r->split);
@@ -259,6 +304,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 +319,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 +370,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 +439,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
+			int slen;
-			else if (flag & 0x800) mm_sprintf_lite(s, "%dS", r->re - r->rs); // supplementary alignment
+			if ((flag & 0x900) == 0 || (opt_flag & MM_F_SOFTCLIP)) slen = t->l_seq;
-			else mm_sprintf_lite(s, "%dS", 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 +523,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);
+			if (r[i].p) {
-			aux[n_aux].x = (uint64_t)r[i].p->dp_max << 32 | r[i].hash;
+				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;
+		max_gap = a0e->x + max_gap > a1s->x? max_gap : a1s->x - a0e->x;
-		min_gap = min_gap < a1s->x - a0e->x? min_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
+		min_flank_len = (int)(max_gap * opt->min_join_flank_ratio);
-		if (r1->re - r1->rs < max_gap>>1 || r1->qe - r1->qs < max_gap>>1) continue;
+		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;
+			aux[n_aux].y = i, aux[n_aux++].x = (uint64_t)regs[i].rid << 32 | regs[i].rs;
-	radix_sort_64(aux, aux + n_aux);
+	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 *l = &regs[aux[i-1].y];
-			mm_reg1_t *r = &regs[(int32_t)aux[i+1]];
+			mm_reg1_t *r = &regs[aux[i+1].y];
 			inv->mapq = l->mapq < r->mapq? l->mapq : r->mapq;
 		}
 	}

index.c

@@ -45,13 +45,15 @@ 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) {
+	if (mi->B) {
-		free(mi->B[i].p);
+		for (i = 0; i < 1U<<mi->b; ++i) {
-		free(mi->B[i].a.a);
+			free(mi->B[i].p);
-		kh_destroy(idx, (idxhash_t*)mi->B[i].h);
+			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)
@@ -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;
@@ -172,7 +175,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 +204,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 +220,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 +340,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 +372,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 +385,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 +423,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];
+	uint32_t x[5], i;
-	int 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;
+		if (mi->seq[i].name) {
-		l = strlen(mi->seq[i].name);
+			uint8_t l = strlen(mi->seq[i].name);
-		fwrite(&l, 1, 1, fp);
+			fwrite(&l, 1, 1, fp);
-		fwrite(mi->seq[i].name, 1, l, 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 +463,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 +478,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);
+		if (l) {
-		fread(s->name, 1, l, fp);
+			s->name = (char*)kmalloc(mi->km, l + 1);
-		s->name[l] = 0;
+			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 +518,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 +576,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;
 }

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((char*)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

ksort.h

@@ -37,7 +37,7 @@ typedef struct {
 	int depth;
 } ks_isort_stack_t;
 
-#define KSORT_SWAP(type_t, a, b) { register type_t t=(a); (a)=(b); (b)=t; }
+#define KSORT_SWAP(type_t, a, b) { type_t t=(a); (a)=(b); (b)=t; }
 
 #define KSORT_INIT(name, type_t, __sort_lt) \
 	void ks_heapdown_##name(size_t i, size_t n, type_t l[]) \

ksw2.h

@@ -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

@@ -17,18 +17,20 @@
 void __cpuidex(int cpuid[4], int func_id, int subfunc_id)
 {
 #if defined(__x86_64__)
-	asm volatile ("cpuid"
+	__asm__ volatile ("cpuid"
 			: "=a" (cpuid[0]), "=b" (cpuid[1]), "=c" (cpuid[2]), "=d" (cpuid[3])
 			: "0" (func_id), "2" (subfunc_id));
 #else // on 32bit, ebx can NOT be used as PIC code
-	asm volatile ("xchgl %%ebx, %1; cpuid; xchgl %%ebx, %1"
+	__asm__ volatile ("xchgl %%ebx, %1; cpuid; xchgl %%ebx, %1"
 			: "=a" (cpuid[0]), "=r" (cpuid[1]), "=c" (cpuid[2]), "=d" (cpuid[3])
 			: "0" (func_id), "2" (subfunc_id));
 #endif
 }
 #endif
 
-int x86_simd(void)
+static int ksw_simd = -1;
+
+static int x86_simd(void)
 {
 	int flag = 0, cpuid[4], max_id;
 	__cpuidex(cpuid, 0, 0);
@@ -54,11 +56,10 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 {
 	extern void ksw_extz2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
 	extern void ksw_extz2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
-	unsigned simd;
+	if (ksw_simd < 0) ksw_simd = x86_simd();
-	simd = x86_simd();
+	if (ksw_simd & SIMD_SSE4_1)
-	if (simd & SIMD_SSE4_1)
 		ksw_extz2_sse41(km, qlen, query, tlen, target, m, mat, q, e, w, zdrop, end_bonus, flag, ez);
-	else if (simd & SIMD_SSE2)
+	else if (ksw_simd & SIMD_SSE2)
 		ksw_extz2_sse2(km, qlen, query, tlen, target, m, mat, q, e, w, zdrop, end_bonus, flag, ez);
 	else abort();
 }
@@ -70,11 +71,10 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				   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);
 	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);
-	unsigned simd;
+	if (ksw_simd < 0) ksw_simd = x86_simd();
-	simd = x86_simd();
+	if (ksw_simd & SIMD_SSE4_1)
-	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);
-	else if (simd & SIMD_SSE2)
+	else if (ksw_simd & SIMD_SSE2)
 		ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, ez);
 	else abort();
 }
@@ -86,11 +86,10 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
 	extern void ksw_exts2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
 				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
-	unsigned simd;
+	if (ksw_simd < 0) ksw_simd = x86_simd();
-	simd = x86_simd();
+	if (ksw_simd & SIMD_SSE4_1)
-	if (simd & SIMD_SSE4_1)
 		ksw_exts2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
-	else if (simd & SIMD_SSE2)
+	else if (ksw_simd & SIMD_SSE2)
 		ksw_exts2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
 	else abort();
 }

ksw2_extd2_sse.c

@@ -76,7 +76,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	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;
@@ -111,7 +111,7 @@ void ksw_extd2_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;
@@ -218,7 +218,7 @@ void ksw_extd2_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, a2, b2, xt1, x2t1, vt1, ut, tmp;
@@ -265,7 +265,7 @@ void ksw_extd2_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, a2, b2, xt1, x2t1, vt1, ut, tmp;
@@ -382,7 +382,7 @@ void ksw_extd2_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) {

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 "ketopt.h"
-#include <getopt.h>
-#else
-#include "getopt.h"
-#endif
 
-#define MM_VERSION "2.10-r761"
+#define MM_VERSION "2.15-r905"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -26,53 +22,62 @@ void liftrlimit()
 void liftrlimit() {}
 #endif
 
-static struct option long_options[] = {
+static ko_longopt_t long_options[] = {
-	{ "bucket-bits",    required_argument, 0, 0 },
+	{ "bucket-bits",    ko_required_argument, 300 },
-	{ "mb-size",        required_argument, 0, 'K' },
+	{ "mb-size",        ko_required_argument, 'K' },
-	{ "seed",           required_argument, 0, 0 },
+	{ "seed",           ko_required_argument, 302 },
-	{ "no-kalloc",      no_argument,       0, 0 },
+	{ "no-kalloc",      ko_no_argument,       303 },
-	{ "print-qname",    no_argument,       0, 0 },
+	{ "print-qname",    ko_no_argument,       304 },
-	{ "no-self",        no_argument,       0, 'D' },
+	{ "no-self",        ko_no_argument,       'D' },
-	{ "print-seeds",    no_argument,       0, 0 },
+	{ "print-seeds",    ko_no_argument,       306 },
-	{ "max-chain-skip", required_argument, 0, 0 },
+	{ "max-chain-skip", ko_required_argument, 307 },
-	{ "min-dp-len",     required_argument, 0, 0 },
+	{ "min-dp-len",     ko_required_argument, 308 },
-	{ "print-aln-seq",  no_argument,       0, 0 },
+	{ "print-aln-seq",  ko_no_argument,       309 },
-	{ "splice",         no_argument,       0, 0 },
+	{ "splice",         ko_no_argument,       310 },
-	{ "cost-non-gt-ag", required_argument, 0, 'C' },
+	{ "cost-non-gt-ag", ko_required_argument, 'C' },
-	{ "no-long-join",   no_argument,       0, 0 },
+	{ "no-long-join",   ko_no_argument,       312 },
-	{ "sr",             no_argument,       0, 0 },
+	{ "sr",             ko_no_argument,       313 },
-	{ "frag",           required_argument, 0, 0 },
+	{ "frag",           ko_required_argument, 314 },
-	{ "secondary",      required_argument, 0, 0 },
+	{ "secondary",      ko_required_argument, 315 },
-	{ "cs",             optional_argument, 0, 0 },
+	{ "cs",             ko_optional_argument, 316 },
-	{ "end-bonus",      required_argument, 0, 0 },
+	{ "end-bonus",      ko_required_argument, 317 },
-	{ "no-pairing",     no_argument,       0, 0 },
+	{ "no-pairing",     ko_no_argument,       318 },
-	{ "splice-flank",   required_argument, 0, 0 },
+	{ "splice-flank",   ko_required_argument, 319 },
-	{ "idx-no-seq",     no_argument,       0, 0 },
+	{ "idx-no-seq",     ko_no_argument,       320 },
-	{ "end-seed-pen",   required_argument, 0, 0 },   // 21
+	{ "end-seed-pen",   ko_required_argument, 321 },
-	{ "for-only",       no_argument,       0, 0 },   // 22
+	{ "for-only",       ko_no_argument,       322 },
-	{ "rev-only",       no_argument,       0, 0 },   // 23
+	{ "rev-only",       ko_no_argument,       323 },
-	{ "heap-sort",      required_argument, 0, 0 },   // 24
+	{ "heap-sort",      ko_required_argument, 324 },
-	{ "all-chain",      no_argument,       0, 'P' },
+	{ "all-chain",      ko_no_argument,       'P' },
-	{ "dual",           required_argument, 0, 0 },   // 26
+	{ "dual",           ko_required_argument, 326 },
-	{ "max-clip-ratio", required_argument, 0, 0 },   // 27
+	{ "max-clip-ratio", ko_required_argument, 327 },
-	{ "min-occ-floor",  required_argument, 0, 0 },   // 28
+	{ "min-occ-floor",  ko_required_argument, 328 },
-	{ "MD",             no_argument,       0, 0 },   // 29
+	{ "MD",             ko_no_argument,       329 },
-	{ "help",           no_argument,       0, 'h' },
+	{ "lj-min-ratio",   ko_required_argument, 330 },
-	{ "max-intron-len", required_argument, 0, 'G' },
+	{ "score-N",        ko_required_argument, 331 },
-	{ "version",        no_argument,       0, 'V' },
+	{ "eqx",            ko_no_argument,       332 },
-	{ "min-count",      required_argument, 0, 'n' },
+	{ "paf-no-hit",     ko_no_argument,       333 },
-	{ "min-chain-score",required_argument, 0, 'm' },
+	{ "split-prefix",   ko_required_argument, 334 },
-	{ "mask-level",     required_argument, 0, 'M' },
+	{ "no-end-flt",     ko_no_argument,       335 },
-	{ "min-dp-score",   required_argument, 0, 's' },
+	{ "hard-mask-level",ko_no_argument,       336 },
-	{ "sam",            no_argument,       0, 'a' },
+	{ "cap-sw-mem",     ko_required_argument, 337 },
-	{ 0, 0, 0, 0}
+	{ "max-qlen",       ko_required_argument, 338 },
+	{ "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,10 +99,11 @@ 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;
+	int i, c, n_threads = 3, n_parts, old_best_n = -1;
 	char *fnw = 0, *rg = 0, *s;
 	FILE *fp_help = stderr;
 	mm_idx_reader_t *idx_rdr;
@@ -108,30 +114,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) {
+			if (mm_set_opt(o.arg, &ipt, &opt) < 0) {
-				fprintf(stderr, "[ERROR] unknown preset '%s'\n", optarg);
+				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 - 1]);
+			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) {
+	while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) {
-		if (c == 'w') ipt.w = atoi(optarg);
+		if (c == 'w') ipt.w = atoi(o.arg);
-		else if (c == 'k') ipt.k = atoi(optarg);
+		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 == 'd') fnw = o.arg; // the above are indexing related options, except -I
-		else if (c == 'r') opt.bw = (int)mm_parse_num(optarg);
+		else if (c == 'r') opt.bw = (int)mm_parse_num(o.arg);
-		else if (c == 't') n_threads = atoi(optarg);
+		else if (c == 't') n_threads = atoi(o.arg);
-		else if (c == 'v') mm_verbose = atoi(optarg);
+		else if (c == 'v') mm_verbose = atoi(o.arg);
-		else if (c == 'g') opt.max_gap = (int)mm_parse_num(optarg);
+		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(optarg));
+		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(optarg);
+		else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(o.arg);
-		else if (c == 'N') opt.best_n = atoi(optarg);
+		else if (c == 'N') old_best_n = opt.best_n, opt.best_n = atoi(o.arg);
-		else if (c == 'p') opt.pri_ratio = atof(optarg);
+		else if (c == 'p') opt.pri_ratio = atof(o.arg);
-		else if (c == 'M') opt.mask_level = atof(optarg);
+		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 +153,68 @@ 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 == 'T') opt.sdust_thres = atoi(o.arg);
-		else if (c == 'n') opt.min_cnt = atoi(optarg);
+		else if (c == 'n') opt.min_cnt = atoi(o.arg);
-		else if (c == 'm') opt.min_chain_score = atoi(optarg);
+		else if (c == 'm') opt.min_chain_score = atoi(o.arg);
-		else if (c == 'A') opt.a = atoi(optarg);
+		else if (c == 'A') opt.a = atoi(o.arg);
-		else if (c == 'B') opt.b = atoi(optarg);
+		else if (c == 'B') opt.b = atoi(o.arg);
-		else if (c == 's') opt.min_dp_max = atoi(optarg);
+		else if (c == 's') opt.min_dp_max = atoi(o.arg);
-		else if (c == 'C') opt.noncan = atoi(optarg);
+		else if (c == 'C') opt.noncan = atoi(o.arg);
-		else if (c == 'I') ipt.batch_size = mm_parse_num(optarg);
+		else if (c == 'I') ipt.batch_size = mm_parse_num(o.arg);
-		else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(optarg);
+		else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(o.arg);
-		else if (c == 'R') rg = optarg;
+		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 == 300) ipt.bucket_bits = atoi(o.arg); // --bucket-bits
-		else if (c == 0 && long_idx == 2) opt.seed = atoi(optarg); // --seed
+		else if (c == 302) opt.seed = atoi(o.arg); // --seed
-		else if (c == 0 && long_idx == 3) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
+		else if (c == 303) 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 == 304) 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 == 306) 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 == 307) opt.max_chain_skip = atoi(o.arg); // --max-chain-skip
-		else if (c == 0 && long_idx == 8) opt.min_ksw_len = atoi(optarg); // --min-dp-len
+		else if (c == 308) opt.min_ksw_len = atoi(o.arg); // --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 == 309) 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 == 310) 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 == 312) 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 == 313) opt.flag |= MM_F_SR; // --sr
-		else if (c == 0 && long_idx ==17) opt.end_bonus = atoi(optarg); // --end-bonus
+		else if (c == 317) opt.end_bonus = atoi(o.arg); // --end-bonus
-		else if (c == 0 && long_idx ==18) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
+		else if (c == 318) 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 == 320) 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 == 321) opt.anchor_ext_shift = atoi(o.arg); // --end-seed-pen
-		else if (c == 0 && long_idx ==22) opt.flag |= MM_F_FOR_ONLY; // --for-only
+		else if (c == 322) 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 == 323) 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 == 327) opt.max_clip_ratio = atof(o.arg); // --max-clip-ratio
-		else if (c == 0 && long_idx ==28) opt.min_mid_occ = atoi(optarg); // --min-occ-floor
+		else if (c == 328) opt.min_mid_occ = atoi(o.arg); // --min-occ-floor
-		else if (c == 0 && long_idx ==29) opt.flag |= MM_F_OUT_MD; // --MD
+		else if (c == 329) opt.flag |= MM_F_OUT_MD; // --MD
-		else if (c == 0 && long_idx == 14) { // --frag
+		else if (c == 330) opt.min_join_flank_ratio = atof(o.arg); // --lj-min-ratio
-			yes_or_no(&opt, MM_F_FRAG_MODE, long_idx, optarg, 1);
+		else if (c == 331) opt.sc_ambi = atoi(o.arg); // --score-N
-		} else if (c == 0 && long_idx == 15) { // --secondary
+		else if (c == 332) opt.flag |= MM_F_EQX; // --eqx
-			yes_or_no(&opt, MM_F_NO_PRINT_2ND, long_idx, optarg, 0);
+		else if (c == 333) opt.flag |= MM_F_PAF_NO_HIT; // --paf-no-hit
-		} else if (c == 0 && long_idx == 16) { // --cs
+		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) opt.max_sw_mat = mm_parse_num(o.arg); // --cap-sw-mat
+		else if (c == 338) opt.max_qlen = mm_parse_num(o.arg); // --max-qlen
+		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
+		} else if (c == 319) { // --splice-flank
-			yes_or_no(&opt, MM_F_SPLICE_FLANK, long_idx, optarg, 1);
+			yes_or_no(&opt, MM_F_SPLICE_FLANK, o.longidx, o.arg, 1);
-		} else if (c == 0 && long_idx == 24) { // --heap-sort
+		} else if (c == 324) { // --heap-sort
-			yes_or_no(&opt, MM_F_HEAP_SORT, long_idx, optarg, 1);
+			yes_or_no(&opt, MM_F_HEAP_SORT, o.longidx, o.arg, 1);
-		} else if (c == 0 && long_idx == 26) { // --dual
+		} else if (c == 326) { // --dual
-			yes_or_no(&opt, MM_F_NO_DUAL, long_idx, optarg, 0);
+			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 +225,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
+			if (*o.arg == '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 (*o.arg == '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 (*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 (*optarg == 'n') opt.flag &= ~(MM_F_SPLICE_FOR|MM_F_SPLICE_REV); // don't try to match the GT-AG signal
+			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 +257,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 +299,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, "    -x STR       preset (always applied before other options; see minimap2.1 for details) []\n");
-		fprintf(fp_help, "                 map-pb: -Hk19 (PacBio vs reference mapping)\n");
+		fprintf(fp_help, "                 - map-pb/map-ont: PacBio/Nanopore vs reference mapping\n");
-		fprintf(fp_help, "                 map-ont: -k15 (Oxford Nanopore vs reference mapping)\n");
+		fprintf(fp_help, "                 - ava-pb/ava-ont: PacBio/Nanopore read overlap\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, "                 - asm5/asm10/asm20: asm-to-ref mapping, for ~0.1/1/5%% sequence divergence\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, "                 - splice: long-read spliced alignment\n");
-		fprintf(fp_help, "                 ava-pb: -Hk19 -Xw5 -m100 -g10000 --max-chain-skip 25 (PacBio read overlap)\n");
+		fprintf(fp_help, "                 - sr: genomic short-read mapping\n");
-		fprintf(fp_help, "                 ava-ont: -k15 -Xw5 -m100 -g10000 -r2000 --max-chain-skip 25 (ONT read overlap)\n");
+		fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of these and other advanced command-line options.\n");
-		fprintf(fp_help, "                 splice: long-read spliced alignment (see minimap2.1 for details)\n");
-		fprintf(fp_help, "                 sr: short single-end reads without splicing (see minimap2.1 for details)\n");
-		fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of command-line options.\n");
 		return fp_help == stdout? 0 : 1;
 	}
 
-	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 +344,29 @@ 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)
+				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.\033[0m\n");
+					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 (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 +377,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;
+		while (u < n_dreg && (int32_t)dreg[u] <= s) ++u;
-		if (u < n_dreg && dreg[u]>>32 < e) {
+		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
+			for (v = u; v < n_dreg && (int32_t)(dreg[v]>>32) < e; ++v) { // iterate over LCRs overlapping this minimizer
-				int ss = s > dreg[v]>>32? s : dreg[v]>>32;
+				int ss = s > (int32_t)(dreg[v]>>32)? s : dreg[v]>>32;
-				int ee = e < (uint32_t)dreg[v]? e : (uint32_t)dreg[v];
+				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 (!skip_seed(opt->flag, r, q, qname, qlen, mi, &is_self)) {
-		if ((r&1) == (q->q_pos&1)) { // forward strand
+			if ((r&1) == (q->q_pos&1)) { // forward strand
-			p = &a[n_for++];
+				p = &a[n_for++];
-			p->x = (r&0xffffffff00000000ULL) | rpos;
+				p->x = (r&0xffffffff00000000ULL) | rpos;
-			p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
+				p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
-		} else { // reverse strand
+			} else { // reverse strand
-			p = &a[(*n_a) - (++n_rev)];
+				p = &a[(*n_a) - (++n_rev)];
-			p->x = 1ULL<<63 | (r&0xffffffff00000000ULL) | rpos;
+				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->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);
 	}
@@ -274,6 +284,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		qlen_sum += qlens[i], n_regs[i] = 0, regs[i] = 0;
 
 	if (qlen_sum == 0 || n_segs <= 0 || n_segs > MM_MAX_SEG) return;
+	if (opt->max_qlen > 0 && qlen_sum > opt->max_qlen) return;
 
 	hash  = qname? __ac_X31_hash_string(qname) : 0;
 	hash ^= __ac_Wang_hash(qlen_sum) + __ac_Wang_hash(opt->seed);
@@ -308,10 +319,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 +337,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 +360,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 +403,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 +434,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 +460,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 +536,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->n_reg = (int*)calloc(5 * 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->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 +551,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 +565,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) {
+				if (p->opt->split_prefix && p->n_parts == 0) { // then write to temporary files
-					mm_reg1_t *r = &s->reg[i][j];
+					mm_err_fwrite(&s->n_reg[i],    sizeof(int), 1, p->fp_split);
-					assert(!r->sam_pri || r->id == r->parent);
+					mm_err_fwrite(&s->rep_len[i],  sizeof(int), 1, p->fp_split);
-					if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent)
+					mm_err_fwrite(&s->frag_gap[i], sizeof(int), 1, p->fp_split);
-						continue;
+					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_write_paf(&p->str, mi, t, 0, 0, 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_err_puts(p->str.s);
 				}
 			}
@@ -502,9 +602,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 +614,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*));
+	pl.fp = open_bseqs(pl.n_fp, fn);
-	for (i = 0; i < n_segs; ++i) {
+	if (pl.fp == 0) return -1;
-		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.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 +661,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
@@ -58,6 +62,7 @@ typedef struct {
 typedef struct {
 	int32_t b, w, k, flag;
 	uint32_t n_seq;            // number of reference sequences
+	int32_t index;
 	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)
@@ -102,6 +107,8 @@ typedef struct {
 	int sdust_thres; // score threshold for SDUST; 0 to disable
 	int flag;        // see MM_F_* macros
 
+	int max_qlen;    // max query length
+
 	int bw;          // bandwidth
 	int max_gap, max_gap_ref; // break a chain if there are no minimizers in a max_gap window
 	int max_frag_len;
@@ -115,8 +122,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 +141,9 @@ 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
+	int64_t max_sw_mat;
+
+	const char *split_prefix;
 } mm_mapopt_t;
 
 // index reader
@@ -216,6 +228,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 +306,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,6 +344,22 @@ 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);

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "27 March 2018" "minimap2-2.10 (r761)" "Bioinformatics tools"
+.TH minimap2 1 "10 January 2019" "minimap2-2.15 (r905)" "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
@@ -262,6 +274,10 @@ Only map to the reverse complement strand of the reference sequences.
 .BR --heap-sort = no | yes
 If yes, sort anchors with heap merge, instead of radix sort. Heap merge is
 faster for short reads, but slower for long reads. [no]
+.TP
+.B --no-pairing
+Treat two reads in a pair as independent reads. The mate related fields in SAM
+are still properly populated.
 .SS Alignment options
 .TP 10
 .BI -A \ INT
@@ -322,6 +338,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 +369,15 @@ 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.
+.TP
+.BI --cap-sw-mem \ NUM
+Skip alignment if the DP matrix size is above
+.IR NUM .
+Set 0 to disable [0].
 .SS Input/output options
 .TP 10
 .B -a
@@ -394,6 +422,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
@@ -427,6 +458,13 @@ memory.
 .BR --secondary = yes | no
 Whether to output secondary alignments [yes]
 .TP
+.BI --max-qlen \ NUM
+Filter out query sequences longer than
+.IR NUM .
+.TP
+.B --paf-no-hit
+In PAF, output query name and length for an unmapped sequence.
+.TP
 .B --version
 Print version number to stdout
 .SS Preset options
@@ -471,7 +509,7 @@ Up to 10% sequence divergence.
 .B asm20
 Long assembly to reference mapping
 .RB ( -k19
-.B -w10 -A1 -B6 -O6,26 -E2,1 -s200 -z200
+.B -w10 -A1 -B4 -O6,26 -E2,1 -s200 -z200
 .BR --min-occ-floor=100 ).
 Up to 20% sequence divergence.
 .TP
@@ -572,6 +610,8 @@ 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
+de	f	Gap-compressed per-base sequence divergence
 .TE
 
 .PP

misc.c

@@ -116,8 +116,7 @@ long peakrss(void)
 double realtime(void)
 {
 	struct timeval tp;
-	struct timezone tzp;
+	gettimeofday(&tp, NULL);
-	gettimeofday(&tp, &tzp);
 	return tp.tv_sec + tp.tv_usec * 1e-6;
 }
 
@@ -131,6 +130,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.15-r905';
 
 /*****************************
  ***** 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 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;
+	var fa_tmp = null, fa, fa_lens, is_vcf = false, sample_name = "sample";
-	while ((c = getopt(args, "l:L:g:q:B:f:")) != null) {
+	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,382 @@ 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_del[3] + " ["+gap_thres+","+gap_thres_long+") deletions");
-	warn(n_ins[3] + " >="+gap_thres+" insertions");
+	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_query_len = 0, min_seg_len = 10000, max_diff = 0.01, bp_flank_len = 0, bp_gap_len = 0;
+	while ((c = getopt(args, "l:d:b:g:q:")) != null) {
+		if (c == 'l') min_seg_len = parseInt(getopt.arg);
+		else if (c == 'd') max_diff = parseFloat(getopt.arg);
+		else if (c == 'b') bp_flank_len = parseInt(getopt.arg);
+		else if (c == 'g') bp_gap_len = parseInt(getopt.arg);
+		else if (c == 'q') min_query_len = parseInt(getopt.arg);
+	}
+	if (getopt.ind == args.length) {
+		print("Usage: paftools.js asmstat [options] <ref.fa.fai> <asm1.paf> [...]");
+		print("Options:");
+		print("  -q INT     ignore query shorter than INT [0]");
+		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, qi) {
+		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]);
+				qi.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', 'l_cov', 'Rcov', 'Rdup', 'Qcov', 'NG75', 'NG50', '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];
+		var label = fn.replace(/.paf(.gz)?$/, "");
+		header.push(label);
+		var ref_blocks = [], qblock_len = [], qblocks = [], bp = [];
+		var query = {}, qinfo = {};
+		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[1] < min_query_len) continue;
+			if (t.length >= 2) {
+				query[t[0]] = t[1];
+				if (qinfo[t[0]] == null) qinfo[t[0]] = {};
+				qinfo[t[0]].len = t[1];
+				qinfo[t[0]].bp = [];
+			}
+			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, qinfo[last_qname]);
+				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, qinfo[last_qname]);
+		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[5][i] = N50(asm_lens, ref_len, 0.75);
+		rst[6][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[1][i] = l_cov;
+		rst[2][i] = (100.0 * (l_cov / ref_len)).toFixed(2) + '%';
+		rst[4][i] = (100.0 * (qcov  / asm_len)).toFixed(2) + '%';
+
+		// compute cov1 and cov2+ lengths; see paf_call() for details
+		var c1_ctg = null, c1_start = 0, c1_end = 0, c1_len = 0;
+		for (var j = 0; j < ref_blocks.length; ++j) {
+			if (ref_blocks[j][0] != c1_ctg || ref_blocks[j][1] >= c1_end) {
+				if (c1_end > c1_start)
+					c1_len += c1_end - c1_start;
+				c1_ctg = ref_blocks[j][0], c1_start = ref_blocks[j][1], c1_end = ref_blocks[j][2];
+			} else if (ref_blocks[j][2] > c1_end) { // overlap
+				if (ref_blocks[j][1] > c1_start)
+					c1_len += ref_blocks[j][1] - c1_start;
+				c1_start = c1_end, c1_end = ref_blocks[j][2];
+			} else if (ref_blocks[j][2] > c1_start) { // contained
+				if (ref_blocks[j][1] > c1_start)
+					c1_len += ref_blocks[j][1] - c1_start;
+				c1_start = ref_blocks[j][2];
+			}
+			//print(ref_blocks[j][0], ref_blocks[j][1], ref_blocks[j][2], c1_start, c1_end, c1_len);
+		}
+		if (c1_end > c1_start)
+			c1_len += c1_end - c1_start;
+		rst[3][i] = (100 * (l_cov - c1_len) / l_cov).toFixed(2) + '%';
+
+		// compute NGA50
+		rst[7][i] = N50(qblock_len, ref_len, 0.5);
+
+		// compute break points
+		rst[8][i] = n_breaks;
+		rst[9][i] = count_bp(bp, 500, 0);
+		rst[10][i] = count_bp(bp, 500, 10000);
+
+		// nb-plot; NOT USED
+		/*
+		var qa = [];
+		for (var qn in qinfo)
+			qa.push([qinfo[qn].len, qinfo[qn].bp]);
+		qa = qa.sort(function(a, b) { return b[0] - a[0] });
+		var sum = 0, n_bp = 0, next_quantile = 0.1;
+		for (var j = 0; j < qa.length; ++j) {
+			sum += qa[j][0];
+			for (var k = 0; k < qa[j][1].length; ++k)
+				if (qa[j][1][k][0] >= bp_flank_len && qa[j][1][k][1] >= bp_gap_len)
+					++n_bp;
+			if (sum >= ref_len * next_quantile) {
+				print(label, Math.floor(next_quantile * 100 + .5), qa[j][0], (sum / n_bp).toFixed(0), n_bp);
+				next_quantile += 0.1;
+				if (next_quantile >= 1.0) break;
+			}
+		}
+		*/
+	}
+	buf.destroy();
+
+	if (bp_flank_len <= 0) {
+		print(header.join("\t"));
+		for (var i = 0; i < labels.length; ++i)
+			print(labels[i], rst[i].join("\t"));
+	}
+}
+
+function paf_asmgene(args)
+{
+	var c, opt = { min_cov:0.99, min_iden:0.99 }, print_err = false, auto_only = false;
+	while ((c = getopt(args, "i:c:ea")) != null)
+		if (c == 'i') opt.min_iden = parseFloat(getopt.arg);
+		else if (c == 'c') opt.min_cov = parseFloat(getopt.arg);
+		else if (c == 'e') print_err = true;
+		else if (c == 'a') auto_only = true;
+
+	var n_fn = args.length - getopt.ind;
+	if (n_fn < 2) {
+		print("Usage: paftools.js asmgene [options] <ref-splice.paf> <asm-splice.paf> [...]");
+		print("Options:");
+		print("  -i FLOAT     min identity [" + opt.min_iden + "]");
+		print("  -c FLOAT     min coverage [" + opt.min_cov + "]");
+		print("  -a           only evaluate genes mapped to the autosomes");
+		print("  -e           print fragmented/missing genes");
+		exit(1);
+	}
+
+	function process_query(opt, a) {
+		var b = [], cnt = [0, 0, 0];
+		for (var j = 0; j < a.length; ++j) {
+			if (a[j][4] < a[j][5] * opt.min_iden)
+				continue;
+			b.push(a[j].slice(0));
+		}
+		if (b.length == 0) return cnt;
+		// count full
+		var n_full = 0;
+		for (var j = 0; j < b.length; ++j)
+			if (b[j][3] - b[j][2] >= b[j][1] * opt.min_cov)
+				++n_full;
+		cnt[0] = n_full;
+		// compute coverage
+		b = b.sort(function(x, y) { return x[2] - y[2] });
+		var l_cov = 0, st = b[0][2], en = b[0][3];
+		for (var j = 1; j < b.length; ++j) {
+			if (b[j][2] <= en)
+				en = b[j][3] > en? b[j][3] : en;
+			else l_cov += en - st;
+		}
+		l_cov += en - st;
+		cnt[1] = l_cov / b[0][1];
+		cnt[2] = b.length;
+		return cnt;
+	}
+
+	var buf = new Bytes();
+	var gene = {}, header = [], refpos = {};
+	for (var i = getopt.ind; i < args.length; ++i) {
+		var fn = args[i];
+		var label = fn.replace(/.paf(.gz)?$/, "");
+		header.push(label);
+		var file = new File(fn), a = [];
+		while (file.readline(buf) >= 0) {
+			var t = buf.toString().split("\t");
+			var ql = parseInt(t[1]), qs = parseInt(t[2]), qe = parseInt(t[3]), mlen = parseInt(t[9]), blen = parseInt(t[10]), mapq = parseInt(t[11]);
+			if (i == getopt.ind) refpos[t[0]] = [t[0], t[1], t[5], t[7], t[8]];
+			if (gene[t[0]] == null) gene[t[0]] = [];
+			if (a.length && t[0] != a[0][0]) {
+				gene[a[0][0]][i - getopt.ind] = process_query(opt, a);
+				a = [];
+			}
+			a.push([t[0], ql, qs, qe, mlen, blen]);
+		}
+		if (a.length)
+			gene[t[0]][i - getopt.ind] = process_query(opt, a);
+		file.close();
+	}
+
+	// select the longest genes (not optimal, but should be good enough)
+	var gene_list = [], gene_nr = {};
+	for (var g in refpos)
+		gene_list.push(refpos[g]);
+	gene_list = gene_list.sort(function(a, b) { return a[2] < b[2]? -1 : a[2] > b[2]? 1 : a[3] - b[3] });
+	var last = 0;
+	for (var j = 1; j < gene_list.length; ++j) {
+		if (gene_list[j][2] != gene_list[last][2] || gene_list[j][3] >= gene_list[last][4]) {
+			gene_nr[gene_list[last][0]] = 1;
+			last = j;
+		} else if (gene_list[j][1] > gene_list[last][1]) {
+			last = j;
+		}
+	}
+	gene_nr[gene_list[last][0]] = 1;
+
+	// count and print
+	var col1 = ["full_sgl", "full_dup", "frag", "part50+", "part10+", "part10-"];
+	var rst = [];
+	for (var k = 0; k < col1.length; ++k) {
+		rst[k] = [];
+		for (var i = 0; i < n_fn; ++i)
+			rst[k][i] = 0;
+	}
+	for (var g in gene) {
+		if (gene[g][0] == null || gene[g][0][0] != 1) continue;
+		if (gene_nr[g] == null) continue;
+		if (auto_only && /^(chr)?[XY]$/.test(refpos[g][2])) continue;
+		for (var i = 0; i < n_fn; ++i) {
+			if (gene[g][i] == null) {
+				rst[4][i]++;
+				if (print_err) print('M', header[i], refpos[g].join("\t"));
+			} else if (gene[g][i][0] == 1) rst[0][i]++;
+			else if (gene[g][i][0] > 1) {
+				rst[1][i]++;
+				if (print_err) print('D', header[i], refpos[g].join("\t"));
+			} else if (gene[g][i][1] >= opt.min_cov) {
+				rst[2][i]++;
+				if (print_err) print('F', header[i], refpos[g].join("\t"));
+			} else if (gene[g][i][1] >= 0.5) {
+				rst[3][i]++;
+				if (print_err) print('5', header[i], refpos[g].join("\t"));
+			} else if (gene[g][i][1] >= 0.1) {
+				rst[4][i]++;
+				if (print_err) print('1', header[i], refpos[g].join("\t"));
+			} else {
+				rst[5][i]++;
+				if (print_err) print('0', header[i], refpos[g].join("\t")); // TODO: reduce code duplicates...
+			}
+		}
+	}
+	print('H', 'Metric', header.join("\t"));
+	for (var k = 0; k < rst.length; ++k) {
+		print('X', col1[k], rst[k].join("\t"));
+	}
+	buf.destroy();
+}
+
 function paf_stat(args)
 {
 	var c, gap_out_len = null;
@@ -676,8 +1049,10 @@ function paf_stat(args)
 			last_qlen = ori_qlen;
 		}
 	}
-	l_tot += last_qlen;
+	if (regs.length) {
-	l_cov += cov_len(regs);
+		l_tot += last_qlen;
+		l_cov += cov_len(regs);
+	}
 
 	file.close();
 	buf.destroy();
@@ -825,6 +1200,105 @@ function paf_bedcov(args)
 	warn("# target bases overlapping regions: " + hit_len + ' (' + (100.0 * hit_len / tot_len).toFixed(2) + '%)');
 }
 
+function paf_vcfpair(args)
+{
+	var c, is_male = false, sample = 'syndip', hgver = null;
+	var PAR = { '37':[[0, 2699520], [154931043, 155260560]] };
+	while ((c = getopt(args, "ms:g:")) != null) {
+		if (c == 'm') is_male = true;
+		else if (c == 's') sample = getopt.arg;
+		else if (c == 'g') hgver = getopt.arg;
+	}
+	if (is_male && (hgver == null || PAR[hgver] == null))
+		throw("for a male, -g must be specified to properly handle PARs on chrX");
+
+	if (getopt.ind == args.length) {
+		print("Usage: paftools.js vcfpair [options] <in.pair.vcf>");
+		print("Options:");
+		print("  -m       the sample is male");
+		print("  -g STR   human genome version '37' []");
+		print("  -s STR   sample name [" + sample + "]");
+		exit(1);
+	}
+
+	var re_ctg = is_male? /^(chr)?([0-9]+|X|Y)$/ : /^(chr)?([0-9]+|X)$/;
+	var label = ['1', '2'];
+	var buf = new Bytes();
+	var file = args[getopt.ind] == '-'? new File() : new File(args[getopt.ind]);
+	while (file.readline(buf) >= 0) {
+		var m, line = buf.toString();
+		if (line.charAt(0) == '#') {
+			if (/^##(source|reference)=/.test(line)) continue;
+			if ((m = /^##contig=.*ID=([^\s,]+)/.exec(line)) != null) {
+				if (!re_ctg.test(m[1])) continue;
+			} else if (/^#CHROM/.test(line)) {
+				var t = line.split("\t");
+				--t.length;
+				t[t.length-1] = sample;
+				line = t.join("\t");
+				print('##FILTER=<ID=HET1,Description="Heterozygous in the first haplotype">');
+				print('##FILTER=<ID=HET2,Description="Heterozygous in the second haplotype">');
+				print('##FILTER=<ID=GAP1,Description="Uncalled in the first haplotype">');
+				print('##FILTER=<ID=GAP2,Description="Uncalled in the second haplotype">');
+			}
+			print(line);
+			continue;
+		}
+		var t = line.split("\t");
+		if (!re_ctg.test(t[0])) continue;
+		var GT = null, AD = null, FILTER = [], HT = [null, null];
+		for (var i = 0; i < 2; ++i) {
+			if ((m = /^(\.|[0-9]+)\/(\.|[0-9]+):(\S+)/.exec(t[9+i])) == null) {
+				warn(line);
+				throw Error("malformatted VCF");
+			}
+			var s = m[3].split(",");
+			if (AD == null) {
+				AD = [];
+				for (var j = 0; j < s.length; ++j)
+					AD[j] = 0;
+			}
+			for (var j = 0; j < s.length; ++j)
+				AD[j] += parseInt(s[j]);
+			if (m[1] == '.') {
+				FILTER.push('GAP' + label[i]);
+				HT[i] = '.';
+			} else if (m[1] != m[2]) {
+				FILTER.push('HET' + label[i]);
+				HT[i] = '.';
+			} else HT[i] = m[1];
+		}
+		--t.length;
+		// test if this is in a haploid region
+		var hap = 0, st = parseInt(t[1]), en = st + t[3].length;
+		if (is_male) {
+			if (/^(chr)?X/.test(t[0])) {
+				if (hgver != null && PAR[hgver] != null) {
+					var r = PAR[hgver], in_par = false;
+					for (var i = 0; i < r.length; ++i)
+						if (r[i][0] <= st && en <= r[i][1])
+							in_par = true;
+					hap = in_par? 0 : 2;
+				}
+			} else if (/^(chr)?Y/.test(t[0])) {
+				hap = 1;
+			}
+		}
+		// special treatment for haploid regions
+		if (hap > 0 && FILTER.length == 1) {
+			if ((hap == 2 && FILTER[0] == "GAP1") || (hap == 1 && FILTER[0] == "GAP2"))
+				FILTER.length = 0;
+		}
+		// update VCF
+		t[5] = 30; // fake QUAL
+		t[6] = FILTER.length? FILTER.join(";") : ".";
+		t[9] = HT.join("|") + ":" + AD.join(",");
+		print(t.join("\t"));
+	}
+	file.close();
+	buf.destroy();
+}
+
 /**********************
  * Conversion related *
  **********************/
@@ -1999,6 +2473,8 @@ 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("  asmgene    evaluate gene completeness (EXPERIMENTAL)");
 		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,7 +2495,10 @@ 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 == 'asmgene') paf_asmgene(args);
 	else if (cmd == 'liftover' || cmd == 'liftOver') paf_liftover(args);
+	else if (cmd == 'vcfpair') paf_vcfpair(args);
 	else if (cmd == 'call') paf_call(args);
 	else if (cmd == 'mapeval') paf_mapeval(args);
 	else if (cmd == 'bedcov') paf_bedcov(args);

mmpriv.h

@@ -28,6 +28,9 @@
 #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)))
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -72,12 +75,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 +89,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");
@@ -147,6 +159,11 @@ int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
 			fprintf(stderr, "[ERROR]\033[1;31m --for-only and --rev-only can't be applied at the same time\033[0m\n");
 		return -3;
 	}
+	if (mo->e <= 0 || mo->q <= 0) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m -O and -E must be positive\033[0m\n");
+		return -1;
+	}
 	if ((mo->q != mo->q2 || mo->e != mo->e2) && !(mo->e > mo->e2 && mo->q + mo->e < mo->q2 + mo->e2)) {
 		if (mm_verbose >= 1)
 			fprintf(stderr, "[ERROR]\033[1;31m dual gap penalties violating E1>E2 and O1+E1<O2+E2\033[0m\n");
@@ -162,5 +179,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

@@ -54,7 +54,7 @@ void mm_set_pe_thru(const int *qlens, int *n_regs, mm_reg1_t **regs)
 	if (n_pri[0] == 1 && n_pri[1] == 1) {
 		mm_reg1_t *p = &regs[0][pri[0]];
 		mm_reg1_t *q = &regs[1][pri[1]];
-		if (p->rid == q->rid && p->rev == q->rev && abs(p->rs - q->rs) < 3 && abs(p->re - p->re) < 3
+		if (p->rid == q->rid && p->rev == q->rev && abs(p->rs - q->rs) < 3 && abs(p->re - q->re) < 3
 			&& ((p->qs == 0 && qlens[1] - q->qe == 0) || (q->qs == 0 && qlens[0] - p->qe == 0)))
 		{
 			p->pe_thru = q->pe_thru = 1;
@@ -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)
+		for (i = 0; i < (int)sc.n; ++i)
-			if ((sc.a[i]>>32) + sub_diff >= max>>32)
+			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

@@ -73,13 +73,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 +101,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 +126,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 ((uint32_t)st >= mi->seq[rid].len || st >= en) return 0;
-	if (en < 0 || en > mi->seq[rid].len)
+	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 +137,16 @@ 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;
+}
+
 #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,8 @@ cdef extern from "minimap.h":
 		int32_t mid_occ
 		int32_t max_occ
 		int mini_batch_size
+		int64_t max_sw_mat
+		const char *split_prefix
 
 	int mm_set_opt(char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 	int mm_verbose
@@ -84,6 +88,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 +111,7 @@ 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 kstring_t:
 		unsigned l, m

python/mappy.pyx

@@ -3,6 +3,8 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
+__version__ = '2.15'
+
 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),
+		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])
+			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):
 		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,33 @@ 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)
+		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)
+			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))
-		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)
+			self.map_opt.mid_occ = 1000 # don't filter high-occ seeds
 
 	def __dealloc__(self):
 		if self._idx is not NULL:
@@ -132,32 +160,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 = []
+			cigar, _cs, _MD = [], '', ''
-			for k in range(h.n_cigar32):
+			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) {
+	while ((c = ketopt(&o, argc, argv, 1, "w:t:", 0)) >= 0) {
-		if (c == 'w') W = atoi(optarg);
+		if (c == 'w') W = atoi(o.arg);
-		else if (c == 't') T = atoi(optarg);
+		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:
+	with open('python/README.rst') as f:
-        return f.read()
+		return f.read()
 
 setup(
 	name = 'mappy',
-	version = '2.10',
+	version = '2.15',
 	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
+		extra_compile_args = extra_compile_args,
-		include_dirs = ['.'],
+		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);
+}