Release minimap2-2.13 (r850)

r849: option -P doesn't work
I don't know why I haven't found it at the beginning.
2026-09-24 16:48:12 +08:00 · 2018-10-11 13:18:31 -04:00 · 2018-10-04 16:11:59 -04:00 · 2018-09-27 15:13:44 -04:00 · 2018-09-27 14:46:26 -04:00 · 2018-09-23 20:27:49 -04:00
39 changed files with 2238 additions and 810 deletions
@@ -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
@@ -1,21 +1,25 @@
-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

-ifeq ($(arm_neon),)
-ifeq ($(sse2only),)
+ifeq ($(arm_neon),) # if arm_neon is not defined
+ifeq ($(sse2only),) # if sse2only is not defined
 	OBJS+=ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o
-else
+else                # if sse2only is defined
 	OBJS+=ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_exts2_sse.o
 endif
-else
-    OBJS+=ksw2_extz2_neon.o ksw2_extd2_neon.o ksw2_exts2_neon.o
-    CFLAGS+=-D_FILE_OFFSET_BITS=64 -mfpu=neon -fsigned-char
-    INCLUDES+=-I sse2neon
+else				# if arm_neon is defined
+	OBJS+=ksw2_extz2_neon.o ksw2_extd2_neon.o ksw2_exts2_neon.o
+    INCLUDES+=-Isse2neon
+ifeq ($(aarch64),)	#if aarch64 is not defined
+	CFLAGS+=-D_FILE_OFFSET_BITS=64 -mfpu=neon -fsigned-char
+else				#if aarch64 is defined
+	CFLAGS+=-D_FILE_OFFSET_BITS=64 -fsigned-char
+endif
 endif

 .PHONY:all extra clean depend
@@ -28,8 +32,8 @@ all:$(PROG)

 extra:all $(PROG_EXTRA)

-minimap2:main.o getopt.o libminimap2.a
-		$(CC) $(CFLAGS) main.o getopt.o -o $@ -L. -lminimap2 $(LIBS)
+minimap2:main.o libminimap2.a
+		$(CC) $(CFLAGS) main.o -o $@ -L. -lminimap2 $(LIBS)

 minimap2-lite:example.o libminimap2.a
 		$(CC) $(CFLAGS) $< -o $@ -L. -lminimap2 $(LIBS)
@@ -37,31 +41,36 @@ minimap2-lite:example.o libminimap2.a
 libminimap2.a:$(OBJS)
 		$(AR) -csru $@ $(OBJS)

-sdust:sdust.c getopt.o kalloc.o kalloc.h kdq.h kvec.h kseq.h sdust.h
-		$(CC) -D_SDUST_MAIN $(CFLAGS) $< getopt.o kalloc.o -o $@ -lz
+sdust:sdust.c kalloc.o kalloc.h kdq.h kvec.h kseq.h ketopt.h sdust.h
+		$(CC) -D_SDUST_MAIN $(CFLAGS) $< kalloc.o -o $@ -lz

 # SSE-specific targets on x86/x86_64

+ifeq ($(arm_neon),)   # if arm_neon is defined, compile this target with the default setting (i.e. no -msse2)
+ksw2_ll_sse.o:ksw2_ll_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) -msse2 $(CPPFLAGS) $(INCLUDES) $< -o $@
+endif
+
 ksw2_extz2_sse41.o:ksw2_extz2_sse.c ksw2.h kalloc.h
-		$(CC) -c -msse4 $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@

 ksw2_extz2_sse2.o:ksw2_extz2_sse.c ksw2.h kalloc.h
-		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse2 -mno-sse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@

 ksw2_extd2_sse41.o:ksw2_extd2_sse.c ksw2.h kalloc.h
-		$(CC) -c -msse4 $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@

 ksw2_extd2_sse2.o:ksw2_extd2_sse.c ksw2.h kalloc.h
-		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse2 -mno-sse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@

 ksw2_exts2_sse41.o:ksw2_exts2_sse.c ksw2.h kalloc.h
-		$(CC) -c -msse4 $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@

 ksw2_exts2_sse2.o:ksw2_exts2_sse.c ksw2.h kalloc.h
-		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse2 -mno-sse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH -DKSW_SSE2_ONLY $(INCLUDES) $< -o $@

 ksw2_dispatch.o:ksw2_dispatch.c ksw2.h
-		$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
+		$(CC) -c $(CFLAGS) -msse4.1 $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@

 # NEON-specific targets on ARM

@@ -90,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
@@ -99,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
@@ -1,3 +1,152 @@
+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)
+---------------------------------
+
+Changes to minimap2:
+
+ * Optionally output the MD tag for compatibility with existing tools (#63,
+   #118 and #137).
+
+ * Use SSE compiler flags more precisely to prevent compiling errors on certain
+   machines (#127).
+
+ * Added option --min-occ-floor to set a minimum occurrence threshold. Presets
+   intended for assembly-to-reference alignment set this option to 100. This
+   option alleviates issues with regions having high copy numbers (#107).
+
+ * Exit with non-zero code on file writing errors (e.g. disk full; #103 and
+   #132).
+
+ * Added option -y to copy FASTA/FASTQ comments in query sequences to the
+   output (#136).
+
+ * Added the asm20 preset for alignments between genomes at 5-10% sequence
+   divergence.
+
+ * Changed the band-width in the ava-ont preset from 500 to 2000. Oxford
+   Nanopore reads may contain long deletion sequencing errors that break
+   chaining.
+
+Changes to mappy, the Python binding:
+
+ * Fixed a typo in Align.seq() (#126).
+
+Changes to paftools.js, the companion script:
+
+ * Command sam2paf now converts the MD tag to cs.
+
+ * Support VCF output for assembly-to-reference variant calling (#109).
+
+This version should produce identical alignment for read overlapping, RNA-seq
+read mapping, and genomic read mapping. We have also added a cook book to show
+the variety uses of minimap2 on real datasets. Please see cookbook.md in the
+minimap2 source code directory.
+
+(2.10: 27 March 2017, r761)
+
+
+
 Release 2.9-r720 (23 February 2018)
 -----------------------------------

@@ -14,9 +14,11 @@ cd minimap2 && make
 # use presets (no test data)
 ./minimap2 -ax map-pb ref.fa pacbio.fq.gz > aln.sam       # PacBio genomic reads
 ./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam         # Oxford Nanopore genomic reads
+./minimap2 -ax asm20 ref.fa pacbio-ccs.fq.gz > aln.sam    # PacBio CCS genomic reads
 ./minimap2 -ax sr ref.fa read1.fa read2.fa > aln.sam      # short genomic paired-end reads
-./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads
-./minimap2 -ax splice -k14 -uf ref.fa reads.fa > aln.sam  # Nanopore Direct RNA-seq
+./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads (strand unknown)
+./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam  # noisy Nanopore Direct RNA-seq
+./minimap2 -ax splice -uf -C5 ref.fa query.fa > aln.sam   # Final PacBio Iso-seq or traditional cDNA
 ./minimap2 -cx asm5 asm1.fa asm2.fa > aln.paf             # intra-species asm-to-asm alignment
 ./minimap2 -x ava-pb reads.fa reads.fa > overlaps.paf     # PacBio read overlap
 ./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf    # Nanopore read overlap
@@ -61,16 +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.9/minimap2-2.9_x64-linux.tar.bz2 \
-  | tar -jxvf -
-./minimap2-2.9_x64-linux/minimap2
+curl -L https://github.com/lh3/minimap2/releases/download/v2.13/minimap2-2.13_x64-linux.tar.bz2 | tar -jxvf -
+./minimap2-2.13_x64-linux/minimap2
 ```
 If you want to compile from the source, you need to have a C compiler, GNU make
 and zlib development files installed. Then type `make` in the source code
@@ -78,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

@@ -137,7 +139,7 @@ Nanopore reads.
 #### <a name="map-long-splice"></a>Map long mRNA/cDNA reads

 ```sh
-minimap2 -ax splice -uf ref.fa iso-seq.fq > aln.sam          # PacBio Iso-seq/traditional cDNA
+minimap2 -ax splice -uf -C5 ref.fa iso-seq.fq > aln.sam      # PacBio Iso-seq/traditional cDNA
 minimap2 -ax splice ref.fa nanopore-cdna.fa > aln.sam        # Nanopore 2D cDNA-seq
 minimap2 -ax splice -uf -k14 ref.fa direct-rna.fq > aln.sam  # Nanopore Direct RNA-seq
 minimap2 -ax splice --splice-flank=no SIRV.fa SIRV-seq.fa    # mapping against SIRV control
@@ -229,7 +231,7 @@ sorting) still work with such BAM records; tools that read CIGAR will
 effectively ignore these records. It has been decided that future tools will
 will seamlessly recognize long-cigar records generated by option `-L`.

-**TD;DR**: if you work with ultra-long reads and use tools that only process
+**TL;DR**: if you work with ultra-long reads and use tools that only process
 BAM files, please add option `-L`.

 #### <a name="cs"></a>The cs optional tag
@@ -254,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

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

@@ -346,6 +351,10 @@ mappy` or [from BioConda][mappyconda] via `conda install -c bioconda mappy`.
  possible to add non-SIMD support, but it would make minimap2 slower by
  several times.

+* Minimap2 does not work with a single query or database sequence ~2
+  billion bases or longer (2,147,483,647 to be exact). The total length of all
+  sequences can well exceed this threshold.
+


 [paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
@@ -362,3 +371,5 @@ mappy` or [from BioConda][mappyconda] via `conda install -c bioconda mappy`.
 [issue]: https://github.com/lh3/minimap2/issues
 [k8]: https://github.com/attractivechaos/k8
 [manpage]: https://lh3.github.io/minimap2/minimap2.html
+[manpage-cs]: https://lh3.github.io/minimap2/minimap2.html#10
+[doi]: https://doi.org/10.1093/bioinformatics/bty191
@@ -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
@@ -147,8 +149,8 @@ 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)
 {
-	uint32_t k, l, toff = 0, qoff = 0;
-	int32_t s = 0, max = 0, qshift, tshift;
+	uint32_t k, l;
+	int32_t s = 0, max = 0, qshift, tshift, toff = 0, qoff = 0;
 	mm_extra_t *p = r->p;
 	if (p == 0) return;
 	mm_fix_cigar(r, qseq, tseq, &qshift, &tshift);
@@ -197,12 +199,12 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) //
 	mm_extra_t *p;
 	if (n_cigar == 0) return;
 	if (r->p == 0) {
-		uint32_t capacity = n_cigar + sizeof(mm_extra_t);
+		uint32_t capacity = n_cigar + sizeof(mm_extra_t)/4;
 		kroundup32(capacity);
 		r->p = (mm_extra_t*)calloc(capacity, 4);
 		r->p->capacity = capacity;
-	} else if (r->p->n_cigar + n_cigar + sizeof(mm_extra_t) > r->p->capacity) {
-		r->p->capacity = r->p->n_cigar + n_cigar + sizeof(mm_extra_t);
+	} else if (r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4 > r->p->capacity) {
+		r->p->capacity = r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4;
 		kroundup32(r->p->capacity);
 		r->p = (mm_extra_t*)realloc(r->p, r->p->capacity * 4);
 	}
@@ -217,6 +219,77 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) //
 	}
 }

+static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq) // written by @armintoepfer
+{
+	uint32_t n_EQX = 0;
+	uint32_t k, l, m, cap, toff = 0, qoff = 0, n_M = 0;
+	mm_extra_t *p;
+	if (r->p == 0) return;
+	for (k = 0; k < r->p->n_cigar; ++k) {
+		uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4;
+		if (op == 0) {
+			while (len > 0) {
+				for (l = 0; l < len && qseq[qoff + l] == tseq[toff + l]; ++l) {} // run of "="; TODO: N<=>N is converted to "="
+				if (l > 0) { ++n_EQX; len -= l; toff += l; qoff += l; }
+
+				for (l = 0; l < len && qseq[qoff + l] != tseq[toff + l]; ++l) {} // run of "X"
+				if (l > 0) { ++n_EQX; len -= l; toff += l; qoff += l; }
+			}
+			++n_M;
+		} else if (op == 1) { // insertion
+			qoff += len;
+		} else if (op == 2) { // deletion
+			toff += len;
+		} else if (op == 3) { // intron
+			toff += len;
+		}
+	}
+	// update in-place if we can
+	if (n_EQX == n_M) {
+		for (k = 0; k < r->p->n_cigar; ++k) {
+			uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4;
+			if (op == 0) r->p->cigar[k] = len << 4 | 7;
+		}
+		return;
+	}
+	// allocate new storage
+	cap = r->p->n_cigar + (n_EQX - n_M) + sizeof(mm_extra_t);
+	kroundup32(cap);
+	p = (mm_extra_t*)calloc(cap, 4);
+	memcpy(p, r->p, sizeof(mm_extra_t));
+	p->capacity = cap;
+	// update cigar while copying
+	toff = qoff = m = 0;
+	for (k = 0; k < r->p->n_cigar; ++k) {
+		uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4;
+		if (op == 0) { // match/mismatch
+			while (len > 0) {
+				// match
+				for (l = 0; l < len && qseq[qoff + l] == tseq[toff + l]; ++l) {}
+				if (l > 0) p->cigar[m++] = l << 4 | 7;
+				len -= l;
+				toff += l, qoff += l;
+				// mismatch
+				for (l = 0; l < len && qseq[qoff + l] != tseq[toff + l]; ++l) {}
+				if (l > 0) p->cigar[m++] = l << 4 | 8;
+				len -= l;
+				toff += l, qoff += l;
+			}
+			continue;
+		} else if (op == 1) { // insertion
+			qoff += len;
+		} else if (op == 2) { // deletion
+			toff += len;
+		} else if (op == 3) { // intron
+			toff += len;
+		}
+		p->cigar[m++] = r->p->cigar[k];
+	}
+	p->n_cigar = m;
+	free(r->p);
+	r->p = p;
+}
+
 static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
 {
 	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
@@ -268,20 +341,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 +376,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 +384,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 +397,42 @@ static void mm_filter_bad_seeds(void *km, int as1, int cnt1, mm128_t *a, int min
 	kfree(km, K);
 }

+static void mm_filter_bad_seeds_alt(void *km, int as1, int cnt1, mm128_t *a, int min_gap, int max_ext)
+{
+	int n, k, *K;
+	K = collect_long_gaps(km, as1, cnt1, a, min_gap, &n);
+	if (K == 0) return;
+	for (k = 0; k < n;) {
+		int i = K[k], l;
+		int gap1 = ((int32_t)a[as1 + i].y - (int32_t)a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - (int32_t)a[as1 + i - 1].x);
+		int re1 = (int32_t)a[as1 + i].x;
+		int qe1 = (int32_t)a[as1 + i].y;
+		gap1 = gap1 > 0? gap1 : -gap1;
+		for (l = k + 1; l < n; ++l) {
+			int j = K[l], gap2, q_span_pre, rs2, qs2, m;
+			if ((int32_t)a[as1 + j].y - qe1 > max_ext || (int32_t)a[as1 + j].x - re1 > max_ext) break;
+			gap2 = ((int32_t)a[as1 + j].y - (int32_t)a[as1 + j - 1].y) - (int32_t)(a[as1 + j].x - a[as1 + j - 1].x);
+			q_span_pre = a[as1 + j - 1].y >> 32 & 0xff;
+			rs2 = (int32_t)a[as1 + j - 1].x + q_span_pre;
+			qs2 = (int32_t)a[as1 + j - 1].x + q_span_pre;
+			m = rs2 - re1 < qs2 - qe1? rs2 - re1 : qs2 - qe1;
+			gap2 = gap2 > 0? gap2 : -gap2;
+			if (m > gap1 + gap2) break;
+			re1 = (int32_t)a[as1 + j].x;
+			qe1 = (int32_t)a[as1 + j].y;
+			gap1 = gap2;
+		}
+		if (l > k + 1) {
+			int j, end = K[l - 1];
+			for (j = K[k]; j < end; ++j)
+				a[as1 + j].y |= MM_SEED_IGNORE;
+			a[as1 + end].y |= MM_SEED_LONG_JOIN;
+		}
+		k = l;
+	}
+	kfree(km, K);
+}
+
 static void mm_fix_bad_ends(const mm_reg1_t *r, const mm128_t *a, int bw, int min_match, int32_t *as, int32_t *cnt)
 {
 	int32_t i, l, m;
@@ -350,7 +469,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 +507,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 +553,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) {
-			mm_fix_bad_ends_splice(km, opt, mi, r, mat, qlen, qseq0, a, &as1, &cnt1);
-		} else {
-			mm_fix_bad_ends(r, a, opt->bw, opt->min_chain_score * 2, &as1, &cnt1);
-		}
+		if (!(opt->flag & MM_F_NO_END_FLT)) {
+			if (is_splice)
+				mm_fix_bad_ends_splice(km, opt, mi, r, mat, qlen, qseq0, a, &as1, &cnt1);
+			else
+				mm_fix_bad_ends(r, a, opt->bw, opt->min_chain_score * 2, &as1, &cnt1);
+		} else as1 = r->as, cnt1 = r->cnt;
 		mm_filter_bad_seeds(km, as1, cnt1, a, 10, 40, opt->max_gap>>1, 10);
+		mm_filter_bad_seeds_alt(km, as1, cnt1, a, 30, opt->max_gap>>1);
 		mm_adjust_minier(mi, qseq0, &a[as1], &rs, &qs);
 		mm_adjust_minier(mi, qseq0, &a[as1 + cnt1 - 1], &re, &qe);
 	}
@@ -473,7 +594,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);
@@ -517,13 +638,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;
@@ -628,6 +749,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	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);
+		if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r, &qseq0[r->rev][qs1], tseq);
 		if (rev && r->p->trans_strand)
 			r->p->trans_strand ^= 3; // flip to the read strand
 	}
@@ -652,7 +774,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];
@@ -686,6 +808,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);
+	if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r_inv, &qseq[q_off], &tseq[t_off]);
 	ret = 1;
 end_align1_inv:
 	kfree(km, tseq);
@@ -755,7 +878,7 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 	*n_regs_ = n_regs;
 	kfree(km, qseq0[0]);
 	kfree(km, ez.cigar);
-	mm_filter_regs(km, opt, qlen, n_regs_, regs);
-	mm_hit_sort_by_dp(km, n_regs_, regs);
+	mm_filter_regs(opt, qlen, n_regs_, regs);
+	mm_hit_sort(km, n_regs_, regs);
 	return regs;
 }
@@ -62,19 +62,20 @@ static inline char *kstrdup(const kstring_t *s)
 	return t;
 }

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

-mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int chunk_size, int with_qual, int frag_mode, int *n_)
+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;
 	kvec_t(mm_bseq1_t) a = {0,0,0};
@@ -91,12 +92,12 @@ mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
 		assert(ks->seq.l <= INT32_MAX);
 		if (a.m == 0) kv_resize(mm_bseq1_t, 0, a, 256);
 		kv_pushp(mm_bseq1_t, 0, a, &s);
-		kseq2bseq(ks, s, with_qual);
+		kseq2bseq(ks, s, with_qual, with_comment);
 		size += s->l_seq;
 		if (size >= chunk_size) {
 			if (frag_mode && a.a[a.n-1].l_seq < CHECK_PAIR_THRES) {
 				while (kseq_read(ks) >= 0) {
-					kseq2bseq(ks, &fp->s, with_qual);
+					kseq2bseq(ks, &fp->s, with_qual, with_comment);
 					if (mm_qname_same(fp->s.name, a.a[a.n-1].name)) {
 						kv_push(mm_bseq1_t, 0, a, fp->s);
 						memset(&fp->s, 0, sizeof(mm_bseq1_t));
@@ -110,12 +111,17 @@ mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
 	return a.a;
 }

+mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int chunk_size, int with_qual, int frag_mode, int *n_)
+{
+	return mm_bseq_read3(fp, chunk_size, with_qual, 0, frag_mode, n_);
+}
+
 mm_bseq1_t *mm_bseq_read(mm_bseq_file_t *fp, int chunk_size, int with_qual, int *n_)
 {
 	return mm_bseq_read2(fp, chunk_size, with_qual, 0, n_);
 }

-mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int *n_)
+mm_bseq1_t *mm_bseq_read_frag2(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int with_comment, int *n_)
 {
 	int i;
 	int64_t size = 0;
@@ -136,7 +142,7 @@ mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int
 		for (i = 0; i < n_fp; ++i) {
 			mm_bseq1_t *s;
 			kv_pushp(mm_bseq1_t, 0, a, &s);
-			kseq2bseq(fp[i]->ks, s, with_qual);
+			kseq2bseq(fp[i]->ks, s, with_qual, with_comment);
 			size += s->l_seq;
 		}
 		if (size >= chunk_size) break;
@@ -145,6 +151,11 @@ mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int
 	return a.a;
 }

+mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int *n_)
+{
+	return mm_bseq_read_frag2(n_fp, fp, chunk_size, with_qual, 0, n_);
+}
+
 int mm_bseq_eof(mm_bseq_file_t *fp)
 {
 	return (ks_eof(fp->ks->f) && fp->s.seq == 0);
@@ -13,13 +13,15 @@ typedef struct mm_bseq_file_s mm_bseq_file_t;

 typedef struct {
 	int l_seq, rid;
-	char *name, *seq, *qual;
+	char *name, *seq, *qual, *comment;
 } mm_bseq1_t;

 mm_bseq_file_t *mm_bseq_open(const char *fn);
 void mm_bseq_close(mm_bseq_file_t *fp);
+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_);
 mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int chunk_size, int with_qual, int frag_mode, int *n_);
 mm_bseq1_t *mm_bseq_read(mm_bseq_file_t *fp, int chunk_size, int with_qual, int *n_);
+mm_bseq1_t *mm_bseq_read_frag2(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int with_comment, int *n_);
 mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int *n_);
 int mm_bseq_eof(mm_bseq_file_t *fp);

@@ -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;
@@ -0,0 +1,243 @@
+## Table of Contents
+
+- [Introduction & Installation](#intro)
+- [Mapping Genomic Reads](#map-reads)
+  * [Mapping long reads](#map-pb)
+  * [Mapping Illumina paired-end reads](#map-sr)
+  * [Evaluating mapping accuracy with simulated reads (for developers)](#mapeval)
+- [Mapping Long RNA-seq Reads](#map-rna)
+  * [Mapping Nanopore 2D cDNA reads](#map-ont-cdna-2d)
+  * [Mapping Nanopore direct-RNA reads](#map-direct-rna)
+  * [Mapping PacBio Iso-seq reads](#map-iso-seq)
+- [Full-Genome Alignment](#genome-aln)
+  * [Intra-species assembly alignment](#asm-to-ref)
+  * [Cross-species full-genome alignment](#x-species)
+  * [Eyeballing alignment](#view-aln)
+  * [Calling variants from assembly-to-reference alignment](#asm-var)
+  * [Constructing self-homology map](#hom-map)
+  * [Lift Over (for developers)](#liftover)
+- [Read Overlap](#read-overlap)
+  * [Long-read overlap](#long-read-overlap)
+  * [Evaluating overlap sensitivity (for developers)](#ov-eval)
+
+## <a name="intro"></a>Introduction & Installation
+
+This cookbook walks you through a variety of applications of minimap2 and its
+companion script `paftools.js`. All data here are freely available from the
+minimap2 release page at version tag [v2.10][v2.10]. Some examples only work
+with v2.10 or later.
+
+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.13/minimap2-2.13_x64-linux.tar.bz2 | tar jxf -
+cp minimap2-2.13_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
+export PATH="$PATH:"`pwd`                               # put the current directory on PATH
+# download example datasets
+curl -L https://github.com/lh3/minimap2/releases/download/v2.10/cookbook-data.tgz | tar zxf -
+```
+
+## <a name="map-reads"></a>Mapping Genomic Reads
+
+### <a name="map-pb"></a>Mapping long reads
+```sh
+minimap2 -ax map-pb -t4 ecoli_ref.fa ecoli_p6_25x_canu.fa > mapped.sam
+```
+Alternatively, you can create a minimap2 index first and then map:
+```sh
+minimap2 -x map-pb -d ecoli-pb.mmi ecoli_ref.fa                      # create an index
+minimap2 -ax map-pb ecoli-pb.mmi ecoli_p6_25x_canu.fa > mapped.sam
+```
+This will save you a couple of minutes when you map against the human genome.
+**HOWEVER**, key algorithm parameters such as the k-mer length and window
+size can't be changed after indexing. Minimap2 will give you a warning if
+parameters used in a pre-built index doesn't match parameters on the command
+line. **Please always make sure you are using an intended pre-built index.**
+
+### <a name="map-sr"></a>Mapping Illumina paired-end reads:
+```sh
+minimap2 -ax sr -t4 ecoli_ref.fa ecoli_mason_1.fq ecoli_mason_2.fq > mapped-sr.sam
+```
+
+### <a name="mapeval"></a>Evaluating mapping accuracy with simulated reads (for developers)
+```sh
+minimap2 -ax sr ecoli_ref.fa ecoli_mason_1.fq ecoli_mason_2.fq | paftools.js mapeval -
+```
+The output is:
+```
+Q       60      19712   0       0.000000000     19712
+Q       0       282     219     0.010953286     19994
+U       6
+```
+where a `U`-line gives the number of unmapped reads (for SAM input only); a
+`Q`-line gives:
+
+1. Mapping quality (mapQ) threshold
+2. Number of mapped reads between this threshold and the previous mapQ threshold.
+3. Number of wrong mappings in the same mapQ interval
+4. Accumulative mapping error rate
+5. Accumulative number of mappings
+
+For `paftools.js mapeval` to work, you need to encode the true read positions
+in read names in the right format. For [PBSIM][pbsim] and [mason2][mason2], we
+provide scripts to generate the right format. Simulated reads in this cookbook
+were created with the following command lines:
+```sh
+# in PBSIM source code directory:
+src/pbsim ../ecoli_ref.fa --depth 1 --sample-fastq sample/sample.fastq
+paftools.js pbsim2fq ../ecoli_ref.fa.fai sd_0001.maf > ../ecoli_pbsim.fa
+
+# mason2 simulation
+mason_simulator --illumina-prob-mismatch-scale 2.5 -ir ecoli_ref.fa -n 10000 -o tmp-l.fq -or tmp-r.fq -oa tmp.sam
+paftools.js mason2fq tmp.sam | seqtk seq -1 > ecoli_mason_1.fq
+paftools.js mason2fq tmp.sam | seqtk seq -2 > ecoli_mason_2.fq
+```
+
+
+
+## <a name="map-rna"></a>Mapping Long RNA-seq Reads
+
+### <a name="map-ont-cdna-2d"></a>Mapping Nanopore 2D cDNA reads
+```sh
+minimap2 -ax splice SIRV_E2.fa SIRV_ont-cdna.fa > aln.sam
+```
+You can compare the alignment to the true annotations with:
+```sh
+paftools.js junceval SIRV_E2C.gtf aln.sam
+```
+It gives the percentage of introns found in the annotation. For SIRV data, it
+is possible to achieve higher junction accuracy with
+```sh
+minimap2 -ax splice --splice-flank=no SIRV_E2.fa SIRV_ont-cdna.fa | paftools.js junceval SIRV_E2C.gtf
+```
+This is because minimap2 models one additional evolutionarily conserved base
+around a canonical junction, but SIRV doesn't honor this signal. Option
+`--splice-flank=no` asks minimap2 no to model this additional base.
+
+In the output a tag `ts:A:+` indicates that the read strand is the same as the
+transcript strand; `ts:A:-` indicates the read strand is opposite to the
+transcript strand. This tag is inferred from the GT-AG signal and is thus only
+available to spliced reads.
+
+### <a name="map-direct-rna"></a>Mapping Nanopore direct-RNA reads
+```sh
+minimap2 -ax splice -k14 -uf SIRV_E2.fa SIRV_ont-drna.fa > aln.sam
+```
+Direct-RNA reads are noisier, so we use a shorter k-mer for improved
+sensitivity. Here, option `-uf` forces minimap2 to map reads to the forward
+transcript strand only because direct-RNA reads are stranded. Again, applying
+`--splice-flank=no` helps junction accuracy for SIRV data.
+
+### <a name="map-iso-seq"></a>Mapping PacBio Iso-seq reads
+```sh
+minimap2 -ax splice -uf -C5 SIRV_E2.fa SIRV_iso-seq.fq > aln.sam
+```
+Option `-C5` reduces the penalty on non-canonical splicing sites. It helps
+to align such sites correctly for data with low error rate such as Iso-seq
+reads and traditional cDNAs. On this example, minimap2 makes one junction
+error. Applying `--splice-flank=no` fixes this alignment error.
+
+Note that the command line above is optimized for the final Iso-seq reads.
+PacBio's Iso-seq pipeline produces intermediate sequences at varying quality.
+For example, some intermediate reads are not stranded. For these reads, option
+`-uf` will lead to more errors. Please revise the minimap2 command line
+accordingly.
+
+
+
+## <a name="genome-aln"></a>Full-Genome Alignment
+
+### <a name="asm-to-ref"></a>Intra-species assembly alignment
+```sh
+# option "--cs" is recommended as paftools.js may need it
+minimap2 -cx asm5 --cs ecoli_ref.fa ecoli_canu.fa > ecoli_canu.paf
+```
+Here `ecoli_canu.fa` is the Canu assembly of `ecoli_p6_25x_canu.fa`. This
+command line outputs alignments in the [PAF format][paf]. Use `-a` instead of
+`-c` to get output in the SAM format.
+
+### <a name="x-species"></a>Cross-species full-genome alignment
+```sh
+minimap2 -cx asm20 --cs ecoli_ref.fa ecoli_O104:H4.fa > ecoli_O104:H4.paf
+sort -k6,6 -k8,8n ecoli_O104:H4.paf | paftools.js call -f ecoli_ref.fa -L10000 -l1000 - > out.vcf
+```
+Minimap2 has three presets for full-genome alignment: "asm5" for sequence
+divergence below 1%, "asm10" for divergence around a couple of percent and
+"asm20" for divergence not more than 10%. In theory, with the right setting,
+minimap2 should work for sequence pairs with sequence divergence up to ~15%,
+but this has not been carefully evaluated.
+
+### <a name="view-aln"></a>Eyeballing alignment
+```sh
+# option "--cs" required; minimap2-r741 or higher required for the "asm20" preset
+minimap2 -cx asm20 --cs ecoli_ref.fa ecoli_O104:H4.fa | paftools.js view - | less -S
+```
+This prints the alignment in a BLAST-like format.
+
+### <a name="asm-var"></a>Calling variants from assembly-to-reference alignment
+```sh
+# don't forget the "--cs" option; otherwise it doesn't work
+minimap2 -cx asm5 --cs ecoli_ref.fa ecoli_canu.fa \
+  | sort -k6,6 -k8,8n \
+  | paftools.js call -f ecoli_ref.fa - > out.vcf
+```
+Without option `-f`, `paftools.js call` outputs in a custom format. In this
+format, lines starting with `R` give the regions covered by one contig only.
+This information is not available in the VCF output.
+
+### <a name="hom-map"></a>Constructing self-homology map
+```sh
+minimap2 -DP -k19 -w19 -m200 ecoli_ref.fa ecoli_ref.fa > out.paf
+```
+Option `-D` asks minimap2 to ignore anchors from perfect self match and `-P`
+outputs all chains. For large nomes, we don't recommend to perform base-level
+alignment (with `-c`, `-a` or `--cs`) when `-P` is applied. This is because
+base-alignment is slow and occasionally gives wrong alignments close to the
+diagonal of a dotter plot. For E. coli, though, base-alignment is still fast.
+
+### <a name="liftover"></a>Lift over (for developers)
+```sh
+minimap2 -cx asm5 --cs ecoli_ref.fa ecoli_canu.fa > ecoli_canu.paf
+echo -e 'tig00000001\t200000\t300000' | paftools.js liftover ecoli_canu.paf -
+```
+This lifts over a region on query sequences to one or multiple regions on
+reference sequences. Note that this paftools.js command may not be efficient
+enough to lift millions of regions.
+
+
+
+## <a name="read-overlap"></a>Read Overlap
+
+### <a name="long-read-overlap"></a>Long read overlap
+```sh
+# For pacbio reads:
+minimap2 -x ava-pb ecoli_p6_25x_canu.fa ecoli_p6_25x_canu.fa > overlap.paf
+# For Nanopore reads (ava-ont also works with PacBio but not as good):
+minimap2 -x ava-ont -r 10000 ecoli_p6_25x_canu.fa ecoli_p6_25x_canu.fa > overlap.paf
+# If you have miniasm installed:
+miniasm -f ecoli_p6_25x_canu.fa overlap.paf > asm.gfa
+```
+Here we explicitly applied `-r 10000`. We are considering to set this as the
+default for the `ava-ont` mode as this seems to improve the contiguity for
+nanopore read assembly (Loman, personal communication).
+
+*Minimap2 doesn't work well with short-read overlap.*
+
+### <a name="ov-eval"></a>Evaluating overlap sensitivity (for developers)
+
+```sh
+# read to reference mapping
+minimap2 -cx map-pb ecoli_ref.fa ecoli_p6_25x_canu.fa > to-ref.paf
+# evaluate overlap sensitivity
+sort -k6,6 -k8,8n to-ref.paf | paftools.js ov-eval - overlap.paf
+```
+You can see that for PacBio reads, minimap2 achieves higher overlap sensitivity
+with `-x ava-pb` (99% vs 93% with `-x ava-ont`).
+
+
+
+[pbsim]: https://github.com/pfaucon/PBSIM-PacBio-Simulator
+[mason2]: https://github.com/seqan/seqan/tree/master/apps/mason2
+[paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
+[v2.10]: https://github.com/lh3/minimap2/releases/tag/v2.10
@@ -118,7 +118,7 @@ void mm_write_sam_hdr(const mm_idx_t *idx, const char *rg, const char *ver, int
 	if (idx) {
 		uint32_t i;
 		for (i = 0; i < idx->n_seq; ++i)
-			printf("@SQ\tSN:%s\tLN:%d\n", idx->seq[i].name, idx->seq[i].len);
+			mm_sprintf_lite(&str, "@SQ\tSN:%s\tLN:%d\n", idx->seq[i].name, idx->seq[i].len);
 	}
 	if (rg) sam_write_rg_line(&str, rg);
 	mm_sprintf_lite(&str, "@PG\tID:minimap2\tPN:minimap2");
@@ -129,36 +129,18 @@ void mm_write_sam_hdr(const mm_idx_t *idx, const char *rg, const char *ver, int
 		for (i = 1; i < argc; ++i)
 			mm_sprintf_lite(&str, " %s", argv[i]);
 	}
-	mm_sprintf_lite(&str, "\n");
-	fputs(str.s, stdout);
+	mm_err_puts(str.s);
 	free(str.s);
 }

-static void write_cs(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, 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)
 {
-	extern unsigned char seq_nt4_table[256];
 	int i, q_off, t_off;
-	uint8_t *qseq, *tseq;
-	char *tmp;
-	if (r->p == 0) return;
-	mm_sprintf_lite(s, "\tcs:Z:");
-	qseq = (uint8_t*)kmalloc(km, r->qe - r->qs);
-	tseq = (uint8_t*)kmalloc(km, r->re - r->rs);
-	tmp = (char*)kmalloc(km, r->re - r->rs > r->qe - r->qs? r->re - r->rs + 1 : r->qe - r->qs + 1);
-	mm_idx_getseq(mi, r->rid, r->rs, r->re, tseq);
-	if (!r->rev) {
-		for (i = r->qs; i < r->qe; ++i)
-			qseq[i - r->qs] = seq_nt4_table[(uint8_t)t->seq[i]];
-	} else {
-		for (i = r->qs; i < r->qe; ++i) {
-			uint8_t c = seq_nt4_table[(uint8_t)t->seq[i]];
-			qseq[r->qe - i - 1] = c >= 4? 4 : 3 - c;
-		}
-	}
-	for (i = q_off = t_off = 0; i < r->p->n_cigar; ++i) {
+	if (write_tag) mm_sprintf_lite(s, "\tcs:Z:");
+	for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
 		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
 		assert(op >= 0 && op <= 3);
-		if (op == 0) {
+		if (op == 0) { // match
 			int l_tmp = 0;
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
@@ -179,17 +161,17 @@ static void write_cs(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_
 				} else mm_sprintf_lite(s, ":%d", l_tmp);
 			}
 			q_off += len, t_off += len;
-		} else if (op == 1) {
+		} else if (op == 1) { // insertion to ref
 			for (j = 0, tmp[len] = 0; j < len; ++j)
 				tmp[j] = "acgtn"[qseq[q_off + j]];
 			mm_sprintf_lite(s, "+%s", tmp);
 			q_off += len;
-		} else if (op == 2) {
+		} else if (op == 2) { // deletion from ref
 			for (j = 0, tmp[len] = 0; j < len; ++j)
 				tmp[j] = "acgtn"[tseq[t_off + j]];
 			mm_sprintf_lite(s, "-%s", tmp);
 			t_off += len;
-		} else {
+		} else { // intron
 			assert(len >= 2);
 			mm_sprintf_lite(s, "~%c%c%d%c%c", "acgtn"[tseq[t_off]], "acgtn"[tseq[t_off+1]],
 				len, "acgtn"[tseq[t_off+len-2]], "acgtn"[tseq[t_off+len-1]]);
@@ -197,9 +179,87 @@ static void write_cs(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_
 		}
 	}
 	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, int write_tag)
+{
+	int i, q_off, t_off, l_MD = 0;
+	if (write_tag) mm_sprintf_lite(s, "\tMD:Z:");
+	for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
+		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
+		assert(op >= 0 && op <= 3);
+		if (op == 0) { // 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]]);
+					l_MD = 0;
+				} else ++l_MD;
+			}
+			q_off += len, t_off += len;
+		} else if (op == 1) { // insertion to ref
+			q_off += len;
+		} else if (op == 2) { // deletion from ref
+			for (j = 0, tmp[len] = 0; j < len; ++j)
+				tmp[j] = "ACGTN"[tseq[t_off + j]];
+			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, int write_tag)
+{
+	extern unsigned char seq_nt4_table[256];
+	int i;
+	uint8_t *qseq, *tseq;
+	char *tmp;
+	if (r->p == 0) return;
+	qseq = (uint8_t*)kmalloc(km, r->qe - r->qs);
+	tseq = (uint8_t*)kmalloc(km, r->re - r->rs);
+	tmp = (char*)kmalloc(km, r->re - r->rs > r->qe - r->qs? r->re - r->rs + 1 : r->qe - r->qs + 1);
+	mm_idx_getseq(mi, r->rid, r->rs, r->re, tseq);
+	if (!r->rev) {
+		for (i = r->qs; i < r->qe; ++i)
+			qseq[i - r->qs] = seq_nt4_table[(uint8_t)t->seq[i]];
+	} else {
+		for (i = r->qs; i < r->qe; ++i) {
+			uint8_t c = seq_nt4_table[(uint8_t)t->seq[i]];
+			qseq[r->qe - i - 1] = c >= 4? 4 : 3 - c;
+		}
+	}
+	if (is_MD) write_MD_core(s, tseq, qseq, r, tmp, write_tag);
+	else write_cs_core(s, tseq, qseq, r, tmp, no_iden, write_tag);
 	kfree(km, qseq); kfree(km, tseq); kfree(km, tmp);
 }

+int mm_gen_cs_or_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int is_MD, int no_iden)
+{
+	mm_bseq1_t t;
+	kstring_t str;
+	str.s = *buf, str.l = 0, str.m = *max_len;
+	t.l_seq = strlen(seq);
+	t.seq = (char*)seq;
+	write_cs_or_MD(km, &str, mi, &t, r, no_iden, is_MD, 0);
+	*max_len = str.m;
+	*buf = str.s;
+	return str.l;
+}
+
+int mm_gen_cs(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int no_iden)
+{
+	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 0, no_iden);
+}
+
+int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq)
+{
+	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 1, 0);
+}
+
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 {
 	int type;
@@ -224,6 +284,10 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag)
 {
 	s->l = 0;
+	if (r == 0) {
+		mm_sprintf_lite(s, "%s\t%d", t->name, t->l_seq);
+		return;
+	}
 	mm_sprintf_lite(s, "%s\t%d\t%d\t%d\t%c\t", t->name, t->l_seq, r->qs, r->qe, "+-"[r->rev]);
 	if (mi->seq[r->rid].name) mm_sprintf_lite(s, "%s", mi->seq[r->rid].name);
 	else mm_sprintf_lite(s, "%d", r->rid);
@@ -235,10 +299,12 @@ 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))
-		write_cs(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG));
+	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, 1);
+	if ((opt_flag & MM_F_COPY_COMMENT) && t->comment)
+		mm_sprintf_lite(s, "\t%s", t->comment);
 }

 static void sam_write_sq(kstring_t *s, char *seq, int l, int rev, int comp)
@@ -284,7 +350,7 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
 			int clip_char = (sam_flag&0x800) && !(opt_flag&MM_F_SOFTCLIP)? 'H' : 'S';
 			if (clip_len[0]) mm_sprintf_lite(s, "%d%c", clip_len[0], clip_char);
 			for (k = 0; k < r->p->n_cigar; ++k)
-				mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDN"[r->p->cigar[k]&0xf]);
+				mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDNSHP=XB"[r->p->cigar[k]&0xf]);
 			if (clip_len[1]) mm_sprintf_lite(s, "%d%c", clip_len[1], clip_char);
 		}
 	}
@@ -353,9 +419,11 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 				cigar_in_tag = 1;
 		}
 		if (cigar_in_tag) {
-			if (flag & 0x100) mm_sprintf_lite(s, "0S"); // secondary alignment
-			else if (flag & 0x800) mm_sprintf_lite(s, "%dS", r->re - r->rs); // supplementary alignment
-			else mm_sprintf_lite(s, "%dS", t->l_seq);
+			int slen;
+			if ((flag & 0x900) == 0 || (opt_flag & MM_F_SOFTCLIP)) slen = t->l_seq;
+			else if (flag & 0x100) slen = 0;
+			else slen = r->qe - r->qs;
+			mm_sprintf_lite(s, "%dS%dN", slen, r->re - r->rs);
 		} else write_sam_cigar(s, flag, 0, t->l_seq, r, opt_flag);
 	}

@@ -434,12 +502,15 @@ 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))
-			write_cs(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG));
+		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, 1);
 		if (cigar_in_tag)
 			write_sam_cigar(s, flag, 1, t->l_seq, r, opt_flag);
 	}

+	if ((opt_flag & MM_F_COPY_COMMENT) && t->comment)
+		mm_sprintf_lite(s, "\t%s", t->comment);
+
 	s->s[s->l] = 0; // we always have room for an extra byte (see str_enlarge)
 }

@@ -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);
-}
@@ -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
@@ -106,7 +106,7 @@ void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a)
 	r->split |= 1, r2->split |= 2;
 }

-void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff) // and compute mm_reg1_t::subsc
+void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff, int hard_mask_level) // and compute mm_reg1_t::subsc
 {
 	int i, j, k, *w;
 	uint64_t *cov;
@@ -118,6 +118,7 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff
 	for (i = 1, k = 1; i < n; ++i) {
 		mm_reg1_t *ri = &r[i];
 		int si = ri->qs, ei = ri->qe, n_cov = 0, uncov_len = 0;
+		if (hard_mask_level) goto skip_uncov;
 		for (j = 0; j < k; ++j) { // traverse existing primary hits to find overlapping hits
 			mm_reg1_t *rp = &r[w[j]];
 			int sj = rp->qs, ej = rp->qe;
@@ -132,18 +133,19 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff
 			int j, x = si;
 			radix_sort_64(cov, cov + n_cov);
 			for (j = 0; j < n_cov; ++j) {
-				if (cov[j]>>32 > x) uncov_len += (cov[j]>>32) - x;
+				if ((int)(cov[j]>>32) > x) uncov_len += (cov[j]>>32) - x;
 				x = (int32_t)cov[j] > x? (int32_t)cov[j] : x;
 			}
 			if (ei > x) uncov_len += ei - x;
 		}
+skip_uncov:
 		for (j = 0; j < k; ++j) { // traverse existing primary hits again
 			mm_reg1_t *rp = &r[w[j]];
 			int sj = rp->qs, ej = rp->qe, min, max, ol;
 			if (ej <= si || sj >= ei) continue; // no overlap
 			min = ej - sj < ei - si? ej - sj : ei - si;
 			max = ej - sj > ei - si? ej - sj : ei - si;
-			ol = si < sj? (ei < sj? 0 : ei < ej? ei - sj : ej - sj) : (ej < si? 0 : ej < ei? ej - si : ei - si); // overlap length
+			ol = si < sj? (ei < sj? 0 : ei < ej? ei - sj : ej - sj) : (ej < si? 0 : ej < ei? ej - si : ei - si); // overlap length; TODO: this can be simplified
 			if ((float)ol / min - (float)uncov_len / max > mask_level) {
 				int cnt_sub = 0;
 				ri->parent = rp->parent;
@@ -164,9 +166,9 @@ set_parent_test:
 	kfree(km, w);
 }

-void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r)
+void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r)
 {
-	int32_t i, n_aux, n = *n_regs;
+	int32_t i, n_aux, n = *n_regs, has_cigar = 0, no_cigar = 0;
 	mm128_t *aux;
 	mm_reg1_t *t;

@@ -175,14 +177,20 @@ void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r)
 	t = (mm_reg1_t*)kmalloc(km, n * sizeof(mm_reg1_t));
 	for (i = n_aux = 0; i < n; ++i) {
 		if (r[i].inv || r[i].cnt > 0) { // squeeze out elements with cnt==0 (soft deleted)
-			assert(r[i].p);
-			aux[n_aux].x = (uint64_t)r[i].p->dp_max << 32 | r[i].hash;
+			if (r[i].p) {
+				aux[n_aux].x = (uint64_t)r[i].p->dp_max << 32 | r[i].hash;
+				has_cigar = 1;
+			} else {
+				aux[n_aux].x = (uint64_t)r[i].score << 32 | r[i].hash;
+				no_cigar = 1;
+			}
 			aux[n_aux++].y = i;
 		} else if (r[i].p) {
 			free(r[i].p);
 			r[i].p = 0;
 		}
 	}
+	assert(has_cigar + no_cigar == 1);
 	radix_sort_128x(aux, aux + n_aux);
 	for (i = n_aux - 1; i >= 0; --i)
 		t[n_aux - 1 - i] = r[aux[i].y];
@@ -246,7 +254,7 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_,
 	}
 }

-void mm_filter_regs(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs)
+void mm_filter_regs(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs)
 { // NB: after this call, mm_reg1_t::parent can be -1 if its parent filtered out
 	int i, k;
 	for (i = k = 0; i < *n_regs; ++i) {
@@ -304,7 +312,7 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
 	for (i = n_aux - 1; i >= 1; --i) {
 		mm_reg1_t *r0 = &regs[(int32_t)aux[i-1]], *r1 = &regs[(int32_t)aux[i]];
 		mm128_t *a0e, *a1s;
-		int max_gap, min_gap, sc_thres;
+		int max_gap, min_gap, sc_thres, min_flank_len;

 		// test
 		if (r0->as + r0->cnt != r1->as) continue; // not adjacent in a[]
@@ -313,13 +321,14 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
 		a1s = &a[r1->as];
 		if (a1s->x <= a0e->x || (int32_t)a1s->y <= (int32_t)a0e->y) continue; // keep colinearity
 		max_gap = min_gap = (int32_t)a1s->y - (int32_t)a0e->y;
-		max_gap = max_gap > a1s->x - a0e->x? max_gap : a1s->x - a0e->x;
-		min_gap = min_gap < a1s->x - a0e->x? min_gap : a1s->x - a0e->x;
+		max_gap = a0e->x + max_gap > a1s->x? max_gap : a1s->x - a0e->x;
+		min_gap = a0e->x + min_gap < a1s->x? min_gap : a1s->x - a0e->x;
 		if (max_gap > opt->max_join_long || min_gap > opt->max_join_short) continue;
 		sc_thres = (int)((float)opt->min_join_flank_sc / opt->max_join_long * max_gap + .499);
 		if (r0->score < sc_thres || r1->score < sc_thres) continue; // require good flanking chains
-		if (r0->re - r0->rs < max_gap>>1 || r0->qe - r0->qs < max_gap>>1) continue; // require enough flanking length
-		if (r1->re - r1->rs < max_gap>>1 || r1->qe - r1->qs < max_gap>>1) continue;
+		min_flank_len = (int)(max_gap * opt->min_join_flank_ratio);
+		if (r0->re - r0->rs < min_flank_len || r0->qe - r0->qs < min_flank_len) continue; // require enough flanking length
+		if (r1->re - r1->rs < min_flank_len || r1->qe - r1->qs < min_flank_len) continue;

 		// all conditions satisfied; join
 		a[r1->as].y |= MM_SEED_LONG_JOIN;
@@ -339,7 +348,7 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r
 					r->parent = regs[r->parent].parent;
 			}
 		}
-		mm_filter_regs(km, opt, qlen, n_regs_, regs);
+		mm_filter_regs(opt, qlen, n_regs_, regs);
 		mm_sync_regs(km, *n_regs_, regs);
 	}
 }
@@ -411,23 +420,23 @@ void mm_seg_free(void *km, int n_segs, mm_seg_t *segs)
 static void mm_set_inv_mapq(void *km, int n_regs, mm_reg1_t *regs)
 {
 	int i, n_aux;
-	uint64_t *aux;
+	mm128_t *aux;
 	if (n_regs < 3) return;
 	for (i = 0; i < n_regs; ++i)
 		if (regs[i].inv) break;
 	if (i == n_regs) return; // no inversion hits

-	aux = (uint64_t*)kmalloc(km, n_regs * 8);
+	aux = (mm128_t*)kmalloc(km, n_regs * 16);
 	for (i = n_aux = 0; i < n_regs; ++i)
 		if (regs[i].parent == i || regs[i].parent < 0)
-			aux[n_aux++] = (uint64_t)regs[i].as << 32 | i;
-	radix_sort_64(aux, aux + n_aux);
+			aux[n_aux].y = i, aux[n_aux++].x = (uint64_t)regs[i].rid << 32 | regs[i].rs;
+	radix_sort_128x(aux, aux + n_aux);

 	for (i = 1; i < n_aux - 1; ++i) {
-		mm_reg1_t *inv = &regs[(int32_t)aux[i]];
+		mm_reg1_t *inv = &regs[aux[i].y];
 		if (inv->inv) {
-			mm_reg1_t *l = &regs[(int32_t)aux[i-1]];
-			mm_reg1_t *r = &regs[(int32_t)aux[i+1]];
+			mm_reg1_t *l = &regs[aux[i-1].y];
+			mm_reg1_t *r = &regs[aux[i+1].y];
 			inv->mapq = l->mapq < r->mapq? l->mapq : r->mapq;
 		}
 	}
@@ -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) {
-		free(mi->B[i].p);
-		free(mi->B[i].a.a);
-		kh_destroy(idx, (idxhash_t*)mi->B[i].h);
+	if (mi->B) {
+		for (i = 0; i < 1U<<mi->b; ++i) {
+			free(mi->B[i].p);
+			free(mi->B[i].a.a);
+			kh_destroy(idx, (idxhash_t*)mi->B[i].h);
+		}
 	}
 	if (!mi->km) {
 		for (i = 0; i < mi->n_seq; ++i)
@@ -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];
-	int i;
+	uint32_t x[5], i;

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

-#define MM_VERSION "2.9-r720"
+#define MM_VERSION "2.13-r850"

 #ifdef __linux__
 #include <sys/resource.h>
@@ -22,51 +22,60 @@ void liftrlimit()
 void liftrlimit() {}
 #endif

-static struct option long_options[] = {
-	{ "bucket-bits",    required_argument, 0, 0 },
-	{ "mb-size",        required_argument, 0, 'K' },
-	{ "seed",           required_argument, 0, 0 },
-	{ "no-kalloc",      no_argument,       0, 0 },
-	{ "print-qname",    no_argument,       0, 0 },
-	{ "no-self",        no_argument,       0, 'D' },
-	{ "print-seeds",    no_argument,       0, 0 },
-	{ "max-chain-skip", required_argument, 0, 0 },
-	{ "min-dp-len",     required_argument, 0, 0 },
-	{ "print-aln-seq",  no_argument,       0, 0 },
-	{ "splice",         no_argument,       0, 0 },
-	{ "cost-non-gt-ag", required_argument, 0, 'C' },
-	{ "no-long-join",   no_argument,       0, 0 },
-	{ "sr",             no_argument,       0, 0 },
-	{ "frag",           required_argument, 0, 0 },
-	{ "secondary",      required_argument, 0, 0 },
-	{ "cs",             optional_argument, 0, 0 },
-	{ "end-bonus",      required_argument, 0, 0 },
-	{ "no-pairing",     no_argument,       0, 0 },
-	{ "splice-flank",   required_argument, 0, 0 },
-	{ "idx-no-seq",     no_argument,       0, 0 },
-	{ "end-seed-pen",   required_argument, 0, 0 },   // 21
-	{ "for-only",       no_argument,       0, 0 },   // 22
-	{ "rev-only",       no_argument,       0, 0 },   // 23
-	{ "heap-sort",      required_argument, 0, 0 },   // 24
-	{ "all-chain",      no_argument,       0, 'P' },
-	{ "dual",           required_argument, 0, 0 },   // 26
-	{ "max-clip-ratio", required_argument, 0, 0 },   // 27
-	{ "help",           no_argument,       0, 'h' },
-	{ "max-intron-len", required_argument, 0, 'G' },
-	{ "version",        no_argument,       0, 'V' },
-	{ "min-count",      required_argument, 0, 'n' },
-	{ "min-chain-score",required_argument, 0, 'm' },
-	{ "mask-level",     required_argument, 0, 'M' },
-	{ "min-dp-score",   required_argument, 0, 's' },
-	{ "sam",            no_argument,       0, 'a' },
-	{ 0, 0, 0, 0}
+static ko_longopt_t long_options[] = {
+	{ "bucket-bits",    ko_required_argument, 300 },
+	{ "mb-size",        ko_required_argument, 'K' },
+	{ "seed",           ko_required_argument, 302 },
+	{ "no-kalloc",      ko_no_argument,       303 },
+	{ "print-qname",    ko_no_argument,       304 },
+	{ "no-self",        ko_no_argument,       'D' },
+	{ "print-seeds",    ko_no_argument,       306 },
+	{ "max-chain-skip", ko_required_argument, 307 },
+	{ "min-dp-len",     ko_required_argument, 308 },
+	{ "print-aln-seq",  ko_no_argument,       309 },
+	{ "splice",         ko_no_argument,       310 },
+	{ "cost-non-gt-ag", ko_required_argument, 'C' },
+	{ "no-long-join",   ko_no_argument,       312 },
+	{ "sr",             ko_no_argument,       313 },
+	{ "frag",           ko_required_argument, 314 },
+	{ "secondary",      ko_required_argument, 315 },
+	{ "cs",             ko_optional_argument, 316 },
+	{ "end-bonus",      ko_required_argument, 317 },
+	{ "no-pairing",     ko_no_argument,       318 },
+	{ "splice-flank",   ko_required_argument, 319 },
+	{ "idx-no-seq",     ko_no_argument,       320 },
+	{ "end-seed-pen",   ko_required_argument, 321 },
+	{ "for-only",       ko_no_argument,       322 },
+	{ "rev-only",       ko_no_argument,       323 },
+	{ "heap-sort",      ko_required_argument, 324 },
+	{ "all-chain",      ko_no_argument,       'P' },
+	{ "dual",           ko_required_argument, 326 },
+	{ "max-clip-ratio", ko_required_argument, 327 },
+	{ "min-occ-floor",  ko_required_argument, 328 },
+	{ "MD",             ko_no_argument,       329 },
+	{ "lj-min-ratio",   ko_required_argument, 330 },
+	{ "score-N",        ko_required_argument, 331 },
+	{ "eqx",            ko_no_argument,       332 },
+	{ "paf-no-hit",     ko_no_argument,       333 },
+	{ "split-prefix",   ko_required_argument, 334 },
+	{ "no-end-flt",     ko_no_argument,       335 },
+	{ "hard-mask-level",ko_no_argument,       336 },
+	{ "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;
@@ -88,10 +97,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:";
+	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;
@@ -102,30 +112,36 @@ int main(int argc, char *argv[])
 	mm_realtime0 = realtime();
 	mm_set_opt(0, &ipt, &opt);

-	while ((c = getopt_long(argc, argv, opt_str, long_options, &long_idx)) >= 0) // apply option -x/preset first
+	while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) { // test command line options and apply option -x/preset first
 		if (c == 'x') {
-			if (mm_set_opt(optarg, &ipt, &opt) < 0) {
-				fprintf(stderr, "[ERROR] unknown preset '%s'\n", optarg);
+			if (mm_set_opt(o.arg, &ipt, &opt) < 0) {
+				fprintf(stderr, "[ERROR] unknown preset '%s'\n", o.arg);
 				return 1;
 			}
-			break;
+		} else if (c == ':') {
+			fprintf(stderr, "[ERROR] missing option argument\n");
+			return 1;
+		} else if (c == '?') {
+			fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[o.i]);
+			return 1;
 		}
-	optreset = 1;
+	}
+	o = KETOPT_INIT;

-	while ((c = getopt_long(argc, argv, opt_str, long_options, &long_idx)) >= 0) {
-		if (c == 'w') ipt.w = atoi(optarg);
-		else if (c == 'k') ipt.k = atoi(optarg);
+	while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) {
+		if (c == 'w') ipt.w = atoi(o.arg);
+		else if (c == 'k') ipt.k = atoi(o.arg);
 		else if (c == 'H') ipt.flag |= MM_I_HPC;
-		else if (c == 'd') fnw = optarg; // the above are indexing related options, except -I
-		else if (c == 'r') opt.bw = (int)mm_parse_num(optarg);
-		else if (c == 't') n_threads = atoi(optarg);
-		else if (c == 'v') mm_verbose = atoi(optarg);
-		else if (c == 'g') opt.max_gap = (int)mm_parse_num(optarg);
-		else if (c == 'G') mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(optarg));
-		else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(optarg);
-		else if (c == 'N') opt.best_n = atoi(optarg);
-		else if (c == 'p') opt.pri_ratio = atof(optarg);
-		else if (c == 'M') opt.mask_level = atof(optarg);
+		else if (c == 'd') fnw = o.arg; // the above are indexing related options, except -I
+		else if (c == 'r') opt.bw = (int)mm_parse_num(o.arg);
+		else if (c == 't') n_threads = atoi(o.arg);
+		else if (c == 'v') mm_verbose = atoi(o.arg);
+		else if (c == 'g') opt.max_gap = (int)mm_parse_num(o.arg);
+		else if (c == 'G') mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(o.arg));
+		else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(o.arg);
+		else if (c == 'N') old_best_n = opt.best_n, opt.best_n = atoi(o.arg);
+		else if (c == 'p') opt.pri_ratio = atof(o.arg);
+		else if (c == 'M') opt.mask_level = atof(o.arg);
 		else if (c == 'c') opt.flag |= MM_F_OUT_CG | MM_F_CIGAR;
 		else if (c == 'D') opt.flag |= MM_F_NO_DIAG;
 		else if (c == 'P') opt.flag |= MM_F_ALL_CHAINS;
@@ -134,57 +150,67 @@ int main(int argc, char *argv[])
 		else if (c == 'Q') opt.flag |= MM_F_NO_QUAL;
 		else if (c == 'Y') opt.flag |= MM_F_SOFTCLIP;
 		else if (c == 'L') opt.flag |= MM_F_LONG_CIGAR;
-		else if (c == 'T') opt.sdust_thres = atoi(optarg);
-		else if (c == 'n') opt.min_cnt = atoi(optarg);
-		else if (c == 'm') opt.min_chain_score = atoi(optarg);
-		else if (c == 'A') opt.a = atoi(optarg);
-		else if (c == 'B') opt.b = atoi(optarg);
-		else if (c == 's') opt.min_dp_max = atoi(optarg);
-		else if (c == 'C') opt.noncan = atoi(optarg);
-		else if (c == 'I') ipt.batch_size = mm_parse_num(optarg);
-		else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(optarg);
-		else if (c == 'R') rg = optarg;
+		else if (c == 'y') opt.flag |= MM_F_COPY_COMMENT;
+		else if (c == 'T') opt.sdust_thres = atoi(o.arg);
+		else if (c == 'n') opt.min_cnt = atoi(o.arg);
+		else if (c == 'm') opt.min_chain_score = atoi(o.arg);
+		else if (c == 'A') opt.a = atoi(o.arg);
+		else if (c == 'B') opt.b = atoi(o.arg);
+		else if (c == 's') opt.min_dp_max = atoi(o.arg);
+		else if (c == 'C') opt.noncan = atoi(o.arg);
+		else if (c == 'I') ipt.batch_size = mm_parse_num(o.arg);
+		else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(o.arg);
+		else if (c == 'R') rg = o.arg;
 		else if (c == 'h') fp_help = stdout;
 		else if (c == '2') opt.flag |= MM_F_2_IO_THREADS;
-		else if (c == 0 && long_idx == 0) ipt.bucket_bits = atoi(optarg); // --bucket-bits
-		else if (c == 0 && long_idx == 2) opt.seed = atoi(optarg); // --seed
-		else if (c == 0 && long_idx == 3) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
-		else if (c == 0 && long_idx == 4) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
-		else if (c == 0 && long_idx == 6) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
-		else if (c == 0 && long_idx == 7) opt.max_chain_skip = atoi(optarg); // --max-chain-skip
-		else if (c == 0 && long_idx == 8) opt.min_ksw_len = atoi(optarg); // --min-dp-len
-		else if (c == 0 && long_idx == 9) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
-		else if (c == 0 && long_idx ==10) opt.flag |= MM_F_SPLICE; // --splice
-		else if (c == 0 && long_idx ==12) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
-		else if (c == 0 && long_idx ==13) opt.flag |= MM_F_SR; // --sr
-		else if (c == 0 && long_idx ==17) opt.end_bonus = atoi(optarg); // --end-bonus
-		else if (c == 0 && long_idx ==18) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
-		else if (c == 0 && long_idx ==20) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
-		else if (c == 0 && long_idx ==21) opt.anchor_ext_shift = atoi(optarg); // --end-seed-pen
-		else if (c == 0 && long_idx ==22) opt.flag |= MM_F_FOR_ONLY; // --for-only
-		else if (c == 0 && long_idx ==23) opt.flag |= MM_F_REV_ONLY; // --rev-only
-		else if (c == 0 && long_idx ==27) opt.max_clip_ratio = atof(optarg); // --max-clip-ratio
-		else if (c == 0 && long_idx == 14) { // --frag
-			yes_or_no(&opt, MM_F_FRAG_MODE, long_idx, optarg, 1);
-		} else if (c == 0 && long_idx == 15) { // --secondary
-			yes_or_no(&opt, MM_F_NO_PRINT_2ND, long_idx, optarg, 0);
-		} else if (c == 0 && long_idx == 16) { // --cs
+		else if (c == 300) ipt.bucket_bits = atoi(o.arg); // --bucket-bits
+		else if (c == 302) opt.seed = atoi(o.arg); // --seed
+		else if (c == 303) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
+		else if (c == 304) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
+		else if (c == 306) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
+		else if (c == 307) opt.max_chain_skip = atoi(o.arg); // --max-chain-skip
+		else if (c == 308) opt.min_ksw_len = atoi(o.arg); // --min-dp-len
+		else if (c == 309) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
+		else if (c == 310) opt.flag |= MM_F_SPLICE; // --splice
+		else if (c == 312) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
+		else if (c == 313) opt.flag |= MM_F_SR; // --sr
+		else if (c == 317) opt.end_bonus = atoi(o.arg); // --end-bonus
+		else if (c == 318) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
+		else if (c == 320) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
+		else if (c == 321) opt.anchor_ext_shift = atoi(o.arg); // --end-seed-pen
+		else if (c == 322) opt.flag |= MM_F_FOR_ONLY; // --for-only
+		else if (c == 323) opt.flag |= MM_F_REV_ONLY; // --rev-only
+		else if (c == 327) opt.max_clip_ratio = atof(o.arg); // --max-clip-ratio
+		else if (c == 328) opt.min_mid_occ = atoi(o.arg); // --min-occ-floor
+		else if (c == 329) opt.flag |= MM_F_OUT_MD; // --MD
+		else if (c == 330) opt.min_join_flank_ratio = atof(o.arg); // --lj-min-ratio
+		else if (c == 331) opt.sc_ambi = atoi(o.arg); // --score-N
+		else if (c == 332) opt.flag |= MM_F_EQX; // --eqx
+		else if (c == 333) opt.flag |= MM_F_PAF_NO_HIT; // --paf-no-hit
+		else if (c == 334) opt.split_prefix = o.arg; // --split-prefix
+		else if (c == 335) opt.flag |= MM_F_NO_END_FLT; // --no-end-flt
+		else if (c == 336) opt.flag |= MM_F_HARD_MLEVEL; // --hard-mask-level
+		else if (c == 314) { // --frag
+			yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1);
+		} else if (c == 315) { // --secondary
+			yes_or_no(&opt, MM_F_NO_PRINT_2ND, o.longidx, o.arg, 0);
+		} else if (c == 316) { // --cs
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR;
-			if (optarg == 0 || strcmp(optarg, "short") == 0) {
+			if (o.arg == 0 || strcmp(o.arg, "short") == 0) {
 				opt.flag &= ~MM_F_OUT_CS_LONG;
-			} else if (strcmp(optarg, "long") == 0) {
+			} else if (strcmp(o.arg, "long") == 0) {
 				opt.flag |= MM_F_OUT_CS_LONG;
-			} else if (strcmp(optarg, "none") == 0) {
+			} else if (strcmp(o.arg, "none") == 0) {
 				opt.flag &= ~MM_F_OUT_CS;
 			} else if (mm_verbose >= 2) {
 				fprintf(stderr, "[WARNING]\033[1;31m --cs only takes 'short' or 'long'. Invalid values are assumed to be 'short'.\033[0m\n");
 			}
-		} else if (c == 0 && long_idx == 19) { // --splice-flank
-			yes_or_no(&opt, MM_F_SPLICE_FLANK, long_idx, optarg, 1);
-		} else if (c == 0 && long_idx == 24) { // --heap-sort
-			yes_or_no(&opt, MM_F_HEAP_SORT, long_idx, optarg, 1);
-		} else if (c == 0 && long_idx == 26) { // --dual
-			yes_or_no(&opt, MM_F_NO_DUAL, long_idx, optarg, 0);
+		} else if (c == 319) { // --splice-flank
+			yes_or_no(&opt, MM_F_SPLICE_FLANK, o.longidx, o.arg, 1);
+		} else if (c == 324) { // --heap-sort
+			yes_or_no(&opt, MM_F_HEAP_SORT, o.longidx, o.arg, 1);
+		} else if (c == 326) { // --dual
+			yes_or_no(&opt, MM_F_NO_DUAL, o.longidx, o.arg, 0);
 		} else if (c == 'S') {
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR | MM_F_OUT_CS_LONG;
 			if (mm_verbose >= 2)
@@ -195,27 +221,27 @@ int main(int argc, char *argv[])
 		} else if (c == 'f') {
 			double x;
 			char *p;
-			x = strtod(optarg, &p);
+			x = strtod(o.arg, &p);
 			if (x < 1.0) opt.mid_occ_frac = x, opt.mid_occ = 0;
 			else opt.mid_occ = (int)(x + .499);
 			if (*p == ',') opt.max_occ = (int)(strtod(p+1, &p) + .499);
 		} else if (c == 'u') {
-			if (*optarg == 'b') opt.flag |= MM_F_SPLICE_FOR|MM_F_SPLICE_REV; // both strands
-			else if (*optarg == 'f') opt.flag |= MM_F_SPLICE_FOR, opt.flag &= ~MM_F_SPLICE_REV; // match GT-AG
-			else if (*optarg == 'r') opt.flag |= MM_F_SPLICE_REV, opt.flag &= ~MM_F_SPLICE_FOR; // match CT-AC (reverse complement of GT-AG)
-			else if (*optarg == 'n') opt.flag &= ~(MM_F_SPLICE_FOR|MM_F_SPLICE_REV); // don't try to match the GT-AG signal
+			if (*o.arg == 'b') opt.flag |= MM_F_SPLICE_FOR|MM_F_SPLICE_REV; // both strands
+			else if (*o.arg == 'f') opt.flag |= MM_F_SPLICE_FOR, opt.flag &= ~MM_F_SPLICE_REV; // match GT-AG
+			else if (*o.arg == 'r') opt.flag |= MM_F_SPLICE_REV, opt.flag &= ~MM_F_SPLICE_FOR; // match CT-AC (reverse complement of GT-AG)
+			else if (*o.arg == 'n') opt.flag &= ~(MM_F_SPLICE_FOR|MM_F_SPLICE_REV); // don't try to match the GT-AG signal
 			else {
 				fprintf(stderr, "[ERROR]\033[1;31m unrecognized cDNA direction\033[0m\n");
 				return 1;
 			}
 		} else if (c == 'z') {
-			opt.zdrop = opt.zdrop_inv = strtol(optarg, &s, 10);
+			opt.zdrop = opt.zdrop_inv = strtol(o.arg, &s, 10);
 			if (*s == ',') opt.zdrop_inv = strtol(s + 1, &s, 10);
 		} else if (c == 'O') {
-			opt.q = opt.q2 = strtol(optarg, &s, 10);
+			opt.q = opt.q2 = strtol(o.arg, &s, 10);
 			if (*s == ',') opt.q2 = strtol(s + 1, &s, 10);
 		} else if (c == 'E') {
-			opt.e = opt.e2 = strtol(optarg, &s, 10);
+			opt.e = opt.e2 = strtol(o.arg, &s, 10);
 			if (*s == ',') opt.e2 = strtol(s + 1, &s, 10);
 		}
 	}
@@ -227,12 +253,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");
@@ -264,35 +294,34 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -R STR       SAM read group line in a format like '@RG\\tID:foo\\tSM:bar' []\n");
 		fprintf(fp_help, "    -c           output CIGAR in PAF\n");
 		fprintf(fp_help, "    --cs[=STR]   output the cs tag; STR is 'short' (if absent) or 'long' [none]\n");
+		fprintf(fp_help, "    --MD         output the MD tag\n");
+		fprintf(fp_help, "    --eqx        write =/X CIGAR operators\n");
 		fprintf(fp_help, "    -Y           use soft clipping for supplementary alignments\n");
 		fprintf(fp_help, "    -t INT       number of threads [%d]\n", n_threads);
 		fprintf(fp_help, "    -K NUM       minibatch size for mapping [500M]\n");
 //		fprintf(fp_help, "    -v INT       verbose level [%d]\n", mm_verbose);
 		fprintf(fp_help, "    --version    show version number\n");
 		fprintf(fp_help, "  Preset:\n");
-		fprintf(fp_help, "    -x STR       preset (always applied before other options) []\n");
-		fprintf(fp_help, "                 map-pb: -Hk19 (PacBio vs reference mapping)\n");
-		fprintf(fp_help, "                 map-ont: -k15 (Oxford Nanopore vs reference mapping)\n");
-		fprintf(fp_help, "                 asm5: -k19 -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200 (asm to ref mapping; break at 5%% div.)\n");
-		fprintf(fp_help, "                 asm10: -k19 -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200 (asm to ref mapping; break at 10%% div.)\n");
-		fprintf(fp_help, "                 ava-pb: -Hk19 -Xw5 -m100 -g10000 --max-chain-skip 25 (PacBio read overlap)\n");
-		fprintf(fp_help, "                 ava-ont: -k15 -Xw5 -m100 -g10000 --max-chain-skip 25 (ONT read overlap)\n");
-		fprintf(fp_help, "                 splice: long-read spliced alignment (see minimap2.1 for details)\n");
-		fprintf(fp_help, "                 sr: short single-end reads without splicing (see minimap2.1 for details)\n");
-		fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of command-line options.\n");
+		fprintf(fp_help, "    -x STR       preset (always applied before other options; see minimap2.1 for details) []\n");
+		fprintf(fp_help, "                 - map-pb/map-ont: PacBio/Nanopore vs reference mapping\n");
+		fprintf(fp_help, "                 - ava-pb/ava-ont: PacBio/Nanopore read overlap\n");
+		fprintf(fp_help, "                 - asm5/asm10/asm20: asm-to-ref mapping, for ~0.1/1/5%% sequence divergence\n");
+		fprintf(fp_help, "                 - splice: long-read spliced alignment\n");
+		fprintf(fp_help, "                 - sr: genomic short-read mapping\n");
+		fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of these and other advanced command-line options.\n");
 		return fp_help == stdout? 0 : 1;
 	}

-	if ((opt.flag & MM_F_SR) && argc - optind > 3) {
+	if ((opt.flag & MM_F_SR) && argc - o.ind > 3) {
 		fprintf(stderr, "[ERROR] incorrect input: in the sr mode, please specify no more than two query files.\n");
 		return 1;
 	}
-	idx_rdr = mm_idx_reader_open(argv[optind], &ipt, fnw);
+	idx_rdr = mm_idx_reader_open(argv[o.ind], &ipt, fnw);
 	if (idx_rdr == 0) {
-		fprintf(stderr, "[ERROR] failed to open file '%s'\n", argv[optind]);
+		fprintf(stderr, "[ERROR] failed to open file '%s'\n", argv[o.ind]);
 		return 1;
 	}
-	if (!idx_rdr->is_idx && fnw == 0 && argc - optind < 2) {
+	if (!idx_rdr->is_idx && fnw == 0 && argc - o.ind < 2) {
 		fprintf(stderr, "[ERROR] missing input: please specify a query file to map or option -d to keep the index\n");
 		mm_idx_reader_close(idx_rdr);
 		return 1;
@@ -311,29 +340,40 @@ int main(int argc, char *argv[])
 				mm_write_sam_hdr(mi, rg, MM_VERSION, argc, argv);
 			} else {
 				mm_write_sam_hdr(0, rg, MM_VERSION, argc, argv);
-				if (mm_verbose >= 2)
-					fprintf(stderr, "[WARNING]\033[1;31m For a multi-part index, no @SQ lines will be outputted.\033[0m\n");
+				if (opt.split_prefix == 0 && mm_verbose >= 2)
+					fprintf(stderr, "[WARNING]\033[1;31m For a multi-part index, no @SQ lines will be outputted. Please use --split-prefix.\033[0m\n");
 			}
 		}
 		if (mm_verbose >= 3)
 			fprintf(stderr, "[M::%s::%.3f*%.2f] loaded/built the index for %d target sequence(s)\n",
 					__func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), mi->n_seq);
-		if (argc != optind + 1) mm_mapopt_update(&opt, mi);
+		if (argc != o.ind + 1) mm_mapopt_update(&opt, mi);
 		if (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);

-	fprintf(stderr, "[M::%s] Version: %s\n", __func__, MM_VERSION);
-	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());
+	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);
+	}
+
+	if (mm_verbose >= 3) {
+		fprintf(stderr, "[M::%s] Version: %s\n", __func__, MM_VERSION);
+		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; Peak RSS: %.3f GB\n", __func__, realtime() - mm_realtime0, cputime(), peakrss() / 1024.0 / 1024.0 / 1024.0);
+	}
 	return 0;
 }
@@ -11,6 +11,7 @@

 struct mm_tbuf_s {
 	void *km;
+	int rep_len, frag_gap;
 };

 mm_tbuf_t *mm_tbuf_init(void)
@@ -28,6 +29,11 @@ void mm_tbuf_destroy(mm_tbuf_t *b)
 	free(b);
 }

+void *mm_tbuf_get_km(mm_tbuf_t *b)
+{
+	return b->km;
+}
+
 static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *seq, int sdust_thres)
 {
 	int n_dreg, j, k, u = 0;
@@ -39,12 +45,12 @@ static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *se
 	for (j = k = 0; j < n; ++j) { // squeeze out minimizers that significantly overlap with LCRs
 		int32_t qpos = (uint32_t)a[j].y>>1, span = a[j].x&0xff;
 		int32_t s = qpos - (span - 1), e = s + span;
-		while (u < n_dreg && (uint32_t)dreg[u] <= s) ++u;
-		if (u < n_dreg && dreg[u]>>32 < e) {
+		while (u < n_dreg && (int32_t)dreg[u] <= s) ++u;
+		if (u < n_dreg && (int32_t)(dreg[u]>>32) < e) {
 			int v, l = 0;
-			for (v = u; v < n_dreg && dreg[v]>>32 < e; ++v) { // iterate over LCRs overlapping this minimizer
-				int ss = s > dreg[v]>>32? s : dreg[v]>>32;
-				int ee = e < (uint32_t)dreg[v]? e : (uint32_t)dreg[v];
+			for (v = u; v < n_dreg && (int32_t)(dreg[v]>>32) < e; ++v) { // iterate over LCRs overlapping this minimizer
+				int ss = s > (int32_t)(dreg[v]>>32)? s : dreg[v]>>32;
+				int ee = e < (int32_t)dreg[v]? e : (uint32_t)dreg[v];
 				l += ee - ss;
 			}
 			if (l <= span>>1) a[k++] = a[j]; // keep the minimizer if less than half of it falls in masked region
@@ -56,9 +62,10 @@ static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *se

 static void collect_minimizers(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, mm128_v *mv)
 {
-	int i, j, n, sum = 0;
+	int i, n, sum = 0;
 	mv->n = 0;
 	for (i = n = 0; i < n_segs; ++i) {
+		size_t j;
 		mm_sketch(km, seqs[i], qlens[i], mi->w, mi->k, i, mi->flag&MM_I_HPC, mv);
 		for (j = n; j < mv->n; ++j)
 			mv->a[j].y += sum << 1;
@@ -81,12 +88,13 @@ typedef struct {

 static mm_match_t *collect_matches(void *km, int *_n_m, int max_occ, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos)
 {
-	int i, rep_st = 0, rep_en = 0, n_m;
+	int rep_st = 0, rep_en = 0, n_m;
+	size_t i;
 	mm_match_t *m;
 	*n_mini_pos = 0;
 	*mini_pos = (uint64_t*)kmalloc(km, mv->n * sizeof(uint64_t));
 	m = (mm_match_t*)kmalloc(km, mv->n * sizeof(mm_match_t));
-	for (i = n_m = 0, *rep_len = 0, *n_a = 0; i < mv->n; ++i) {
+	for (i = 0, n_m = 0, *rep_len = 0, *n_a = 0; i < mv->n; ++i) {
 		const uint64_t *cr;
 		mm128_t *p = &mv->a[i];
 		uint32_t q_pos = (uint32_t)p->y, q_span = p->x & 0xff;
@@ -120,7 +128,7 @@ static inline int skip_seed(int flag, uint64_t r, const mm_match_t *q, const cha
 		const mm_idx_seq_t *s = &mi->seq[r>>32];
 		int cmp;
 		cmp = strcmp(qname, s->name);
-		if ((flag&MM_F_NO_DIAG) && cmp == 0 && s->len == qlen) {
+		if ((flag&MM_F_NO_DIAG) && cmp == 0 && (int)s->len == qlen) {
 			if ((uint32_t)r>>1 == (q->q_pos>>1)) return 1; // avoid the diagnonal anchors
 			if ((r&1) == (q->q_pos&1)) *is_self = 1; // this flag is used to avoid spurious extension on self chain
 		}
@@ -163,19 +171,20 @@ static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max
 		mm128_t *p;
 		uint64_t r = heap->x;
 		int32_t is_self, rpos = (uint32_t)r >> 1;
-		if (skip_seed(opt->flag, r, q, qname, qlen, mi, &is_self)) continue;
-		if ((r&1) == (q->q_pos&1)) { // forward strand
-			p = &a[n_for++];
-			p->x = (r&0xffffffff00000000ULL) | rpos;
-			p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
-		} else { // reverse strand
-			p = &a[(*n_a) - (++n_rev)];
-			p->x = 1ULL<<63 | (r&0xffffffff00000000ULL) | rpos;
-			p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
+		if (!skip_seed(opt->flag, r, q, qname, qlen, mi, &is_self)) {
+			if ((r&1) == (q->q_pos&1)) { // forward strand
+				p = &a[n_for++];
+				p->x = (r&0xffffffff00000000ULL) | rpos;
+				p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
+			} else { // reverse strand
+				p = &a[(*n_a) - (++n_rev)];
+				p->x = 1ULL<<63 | (r&0xffffffff00000000ULL) | rpos;
+				p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
+			}
+			p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
+			if (q->is_tandem) p->y |= MM_SEED_TANDEM;
+			if (is_self) p->y |= MM_SEED_SELF;
 		}
-		p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
-		if (q->is_tandem) p->y |= MM_SEED_TANDEM;
-		if (is_self) p->y |= MM_SEED_SELF;
 		// update the heap
 		if ((uint32_t)heap->y < q->n - 1) {
 			++heap[0].y;
@@ -205,7 +214,7 @@ static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max
 static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, const mm128_v *mv, int qlen, int64_t *n_a, int *rep_len,
 								  int *n_mini_pos, uint64_t **mini_pos)
 {
-	int i, k, n_m;
+	int i, n_m;
 	mm_match_t *m;
 	mm128_t *a;
 	m = collect_matches(km, &n_m, max_occ, mi, mv, n_a, rep_len, n_mini_pos, mini_pos);
@@ -213,6 +222,7 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ,
 	for (i = 0, *n_a = 0; i < n_m; ++i) {
 		mm_match_t *q = &m[i];
 		const uint64_t *r = q->cr;
+		uint32_t k;
 		for (k = 0; k < q->n; ++k) {
 			int32_t is_self, rpos = (uint32_t)r[k] >> 1;
 			mm128_t *p;
@@ -238,7 +248,7 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ,
 static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_idx_t *mi, void *km, int qlen, int n_segs, const int *qlens, int *n_regs, mm_reg1_t *regs, mm128_t *a)
 {
 	if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
-		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b);
+		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL);
 		if (n_segs <= 1) mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 		else mm_select_sub_multi(km, opt->pri_ratio, 0.2f, 0.7f, max_chain_gap_ref, mi->k*2, opt->best_n, n_segs, qlens, n_regs, regs);
 		if (!(opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN))) // long join not working well without primary chains
@@ -251,7 +261,7 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k
 	if (!(opt->flag & MM_F_CIGAR)) return regs;
 	regs = mm_align_skeleton(km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs()
 	if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
-		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b);
+		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL);
 		mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 		mm_set_sam_pri(*n_regs, regs);
 	}
@@ -308,10 +318,10 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		if (n_regs0 > 0) { // test if the best chain has all the segments
 			int n_chained_segs = 1, max = 0, max_i = -1, max_off = -1, off = 0;
 			for (i = 0; i < n_regs0; ++i) { // find the best chain
-				if (max < u[i]>>32) max = u[i]>>32, max_i = i, max_off = off;
+				if (max < (int)(u[i]>>32)) max = u[i]>>32, max_i = i, max_off = off;
 				off += (uint32_t)u[i];
 			}
-			for (i = 1; i < (uint32_t)u[max_i]; ++i) // count the number of segments in the best chain
+			for (i = 1; i < (int32_t)u[max_i]; ++i) // count the number of segments in the best chain
 				if ((a[max_off+i].y&MM_SEED_SEG_MASK) != (a[max_off+i-1].y&MM_SEED_SEG_MASK))
 					++n_chained_segs;
 			if (n_chained_segs < n_segs)
@@ -326,6 +336,8 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 			a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 		}
 	}
+	b->frag_gap = max_chain_gap_ref;
+	b->rep_len = rep_len;

 	regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a);

@@ -347,7 +359,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		seg = mm_seg_gen(b->km, hash, n_segs, qlens, n_regs0, regs0, n_regs, regs, a); // split fragment chain to separate segment chains
 		free(regs0);
 		for (i = 0; i < n_segs; ++i) {
-			mm_set_parent(b->km, opt->mask_level, n_regs[i], regs[i], opt->a * 2 + opt->b); // update mm_reg1_t::parent
+			mm_set_parent(b->km, opt->mask_level, n_regs[i], regs[i], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL); // update mm_reg1_t::parent
 			regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a);
 			mm_set_mapq(b->km, n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr);
 		}
@@ -390,13 +402,17 @@ typedef struct {
 	mm_bseq_file_t **fp;
 	const mm_idx_t *mi;
 	kstring_t str;
+
+	int n_parts;
+	uint32_t *rid_shift;
+	FILE *fp_split, **fp_parts;
 } pipeline_t;

 typedef struct {
 	const pipeline_t *p;
    int n_seq, n_frag;
 	mm_bseq1_t *seq;
-	int *n_reg, *seg_off, *n_seg;
+	int *n_reg, *seg_off, *n_seg, *rep_len, *frag_gap;
 	mm_reg1_t **reg;
 	mm_tbuf_t **buf;
 } step_t;
@@ -417,10 +433,17 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
 		qseqs[j] = s->seq[off + j].seq;
 	}
 	if (s->p->opt->flag & MM_F_INDEPEND_SEG) {
-		for (j = 0; j < s->n_seg[i]; ++j)
+		for (j = 0; j < s->n_seg[i]; ++j) {
 			mm_map_frag(s->p->mi, 1, &qlens[j], &qseqs[j], &s->n_reg[off+j], &s->reg[off+j], b, s->p->opt, s->seq[off+j].name);
+			s->rep_len[off + j] = b->rep_len;
+			s->frag_gap[off + j] = b->frag_gap;
+		}
 	} else {
 		mm_map_frag(s->p->mi, s->n_seg[i], qlens, qseqs, &s->n_reg[off], &s->reg[off], b, s->p->opt, s->seq[off].name);
+		for (j = 0; j < s->n_seg[i]; ++j) {
+			s->rep_len[off + j] = b->rep_len;
+			s->frag_gap[off + j] = b->frag_gap;
+		}
 	}
 	for (j = 0; j < s->n_seg[i]; ++j) // flip the query strand and coordinate to the original read strand
 		if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1)))) {
@@ -436,17 +459,75 @@ 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;
    pipeline_t *p = (pipeline_t*)shared;
    if (step == 0) { // step 0: read sequences
 		int with_qual = (!!(p->opt->flag & MM_F_OUT_SAM) && !(p->opt->flag & MM_F_NO_QUAL));
+		int with_comment = !!(p->opt->flag & MM_F_COPY_COMMENT);
 		int frag_mode = (p->n_fp > 1 || !!(p->opt->flag & MM_F_FRAG_MODE));
        step_t *s;
        s = (step_t*)calloc(1, sizeof(step_t));
-		if (p->n_fp > 1) s->seq = mm_bseq_read_frag(p->n_fp, p->fp, p->mini_batch_size, with_qual, &s->n_seq);
-		else s->seq = mm_bseq_read2(p->fp[0], p->mini_batch_size, with_qual, frag_mode, &s->n_seq);
+		if (p->n_fp > 1) s->seq = mm_bseq_read_frag2(p->n_fp, p->fp, p->mini_batch_size, with_qual, with_comment, &s->n_seq);
+		else s->seq = mm_bseq_read3(p->fp[0], p->mini_batch_size, with_qual, with_comment, frag_mode, &s->n_seq);
 		if (s->seq) {
 			s->p = p;
 			for (i = 0; i < s->n_seq; ++i)
@@ -454,9 +535,11 @@ static void *worker_pipeline(void *shared, int step, void *in)
 			s->buf = (mm_tbuf_t**)calloc(p->n_threads, sizeof(mm_tbuf_t*));
 			for (i = 0; i < p->n_threads; ++i)
 				s->buf[i] = mm_tbuf_init();
-			s->n_reg = (int*)calloc(3 * s->n_seq, sizeof(int));
-			s->seg_off = s->n_reg + s->n_seq; // seg_off and n_seg are allocated together with n_reg
+			s->n_reg = (int*)calloc(5 * s->n_seq, sizeof(int));
+			s->seg_off = s->n_reg + s->n_seq; // seg_off, n_seg, rep_len and frag_gap are allocated together with n_reg
 			s->n_seg = s->seg_off + s->n_seq;
+			s->rep_len = s->n_seg + s->n_seq;
+			s->frag_gap = s->rep_len + s->n_seq;
 			s->reg = (mm_reg1_t**)calloc(s->n_seq, sizeof(mm_reg1_t*));
 			for (i = 1, j = 0; i <= s->n_seq; ++i)
 				if (i == s->n_seq || !frag_mode || !mm_qname_same(s->seq[i-1].name, s->seq[i].name)) {
@@ -467,7 +550,8 @@ static void *worker_pipeline(void *shared, int step, void *in)
 			return s;
 		} else free(s);
    } else if (step == 1) { // step 1: map
-		kt_for(p->n_threads, worker_for, in, ((step_t*)in)->n_frag);
+		if (p->n_parts > 0) merge_hits((step_t*)in);
+		else kt_for(p->n_threads, worker_for, in, ((step_t*)in)->n_frag);
 		return in;
    } else if (step == 2) { // step 2: output
 		void *km = 0;
@@ -480,20 +564,36 @@ static void *worker_pipeline(void *shared, int step, void *in)
 			int seg_st = s->seg_off[k], seg_en = s->seg_off[k] + s->n_seg[k];
 			for (i = seg_st; i < seg_en; ++i) {
 				mm_bseq1_t *t = &s->seq[i];
-				for (j = 0; j < s->n_reg[i]; ++j) {
-					mm_reg1_t *r = &s->reg[i][j];
-					assert(!r->sam_pri || r->id == r->parent);
-					if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent)
-						continue;
+				if (p->opt->split_prefix && p->n_parts == 0) { // then write to temporary files
+					mm_err_fwrite(&s->n_reg[i],    sizeof(int), 1, p->fp_split);
+					mm_err_fwrite(&s->rep_len[i],  sizeof(int), 1, p->fp_split);
+					mm_err_fwrite(&s->frag_gap[i], sizeof(int), 1, p->fp_split);
+					for (j = 0; j < s->n_reg[i]; ++j) {
+						mm_reg1_t *r = &s->reg[i][j];
+						mm_err_fwrite(r, sizeof(mm_reg1_t), 1, p->fp_split);
+						if (p->opt->flag & MM_F_CIGAR) {
+							mm_err_fwrite(&r->p->capacity, 4, 1, p->fp_split);
+							mm_err_fwrite(r->p, r->p->capacity, 4, p->fp_split);
+						}
+					}
+				} else if (s->n_reg[i] > 0) { // the query has at least one hit
+					for (j = 0; j < s->n_reg[i]; ++j) {
+						mm_reg1_t *r = &s->reg[i][j];
+						assert(!r->sam_pri || r->id == r->parent);
+						if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent)
+							continue;
+						if (p->opt->flag & MM_F_OUT_SAM)
+							mm_write_sam2(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
+						else
+							mm_write_paf(&p->str, mi, t, r, km, p->opt->flag);
+						mm_err_puts(p->str.s);
+					}
+				} else if (p->opt->flag & (MM_F_OUT_SAM|MM_F_PAF_NO_HIT)) { // output an empty hit, if requested
 					if (p->opt->flag & MM_F_OUT_SAM)
-						mm_write_sam2(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
+						mm_write_sam2(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
 					else
-						mm_write_paf(&p->str, mi, t, r, km, p->opt->flag);
-					puts(p->str.s);
-				}
-				if (s->n_reg[i] == 0 && (p->opt->flag & MM_F_OUT_SAM)) {
-					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);
-					puts(p->str.s);
+						mm_write_paf(&p->str, mi, t, 0, 0, p->opt->flag);
+					mm_err_puts(p->str.s);
 				}
 			}
 			for (i = seg_st; i < seg_en; ++i) {
@@ -501,9 +601,10 @@ static void *worker_pipeline(void *shared, int step, void *in)
 				free(s->reg[i]);
 				free(s->seq[i].seq); free(s->seq[i].name);
 				if (s->seq[i].qual) free(s->seq[i].qual);
+				if (s->seq[i].comment) free(s->seq[i].comment);
 			}
 		}
-		free(s->reg); free(s->n_reg); free(s->seq); // seg_off and n_seg were allocated with reg; no memory leak here
+		free(s->reg); free(s->n_reg); free(s->seq); // seg_off, n_seg, rep_len and frag_gap were allocated with reg; no memory leak here
 		km_destroy(km);
 		if (mm_verbose >= 3)
 			fprintf(stderr, "[M::%s::%.3f*%.2f] mapped %d sequences\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), s->n_seq);
@@ -512,32 +613,44 @@ static void *worker_pipeline(void *shared, int step, void *in)
    return 0;
 }

+static mm_bseq_file_t **open_bseqs(int n, const char **fn)
+{
+	mm_bseq_file_t **fp;
+	int i, j;
+	fp = (mm_bseq_file_t**)calloc(n, sizeof(mm_bseq_file_t*));
+	for (i = 0; i < n; ++i) {
+		if ((fp[i] = mm_bseq_open(fn[i])) == 0) {
+			if (mm_verbose >= 1)
+				fprintf(stderr, "ERROR: failed to open file '%s'\n", fn[i]);
+			for (j = 0; j < i; ++j)
+				mm_bseq_close(fp[j]);
+			free(fp);
+			return 0;
+		}
+	}
+	return fp;
+}
+
 int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_mapopt_t *opt, int n_threads)
 {
-	int i, j, pl_threads;
+	int i, pl_threads;
 	pipeline_t pl;
 	if (n_segs < 1) return -1;
 	memset(&pl, 0, sizeof(pipeline_t));
 	pl.n_fp = n_segs;
-	pl.fp = (mm_bseq_file_t**)calloc(n_segs, sizeof(mm_bseq_file_t*));
-	for (i = 0; i < n_segs; ++i) {
-		pl.fp[i] = mm_bseq_open(fn[i]);
-		if (pl.fp[i] == 0) {
-			if (mm_verbose >= 1)
-				fprintf(stderr, "ERROR: failed to open file '%s'\n", fn[i]);
-			for (j = 0; j < i; ++j)
-				mm_bseq_close(pl.fp[j]);
-			free(pl.fp);
-			return -1;
-		}
-	}
+	pl.fp = open_bseqs(pl.n_fp, fn);
+	if (pl.fp == 0) return -1;
 	pl.opt = opt, pl.mi = idx;
 	pl.n_threads = n_threads > 1? n_threads : 1;
 	pl.mini_batch_size = opt->mini_batch_size;
+	if (opt->split_prefix)
+		pl.fp_split = mm_split_init(opt->split_prefix, idx);
 	pl_threads = n_threads == 1? 1 : (opt->flag&MM_F_2_IO_THREADS)? 3 : 2;
 	kt_pipeline(pl_threads, worker_pipeline, &pl, 3);
+
 	free(pl.str.s);
-	for (i = 0; i < n_segs; ++i)
+	if (pl.fp_split) fclose(pl.fp_split);
+	for (i = 0; i < pl.n_fp; ++i)
 		mm_bseq_close(pl.fp[i]);
 	free(pl.fp);
 	return 0;
@@ -547,3 +660,48 @@ int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int
 {
 	return mm_map_file_frag(idx, 1, &fn, opt, n_threads);
 }
+
+int mm_split_merge(int n_segs, const char **fn, const mm_mapopt_t *opt, int n_split_idx)
+{
+	int i;
+	pipeline_t pl;
+	mm_idx_t *mi;
+	if (n_segs < 1 || n_split_idx < 1) return -1;
+	memset(&pl, 0, sizeof(pipeline_t));
+	pl.n_fp = n_segs;
+	pl.fp = open_bseqs(pl.n_fp, fn);
+	if (pl.fp == 0) return -1;
+	pl.opt = opt;
+	pl.mini_batch_size = opt->mini_batch_size;
+
+	pl.n_parts = n_split_idx;
+	pl.fp_parts  = CALLOC(FILE*, pl.n_parts);
+	pl.rid_shift = CALLOC(uint32_t, pl.n_parts);
+	pl.mi = mi = mm_split_merge_prep(opt->split_prefix, n_split_idx, pl.fp_parts, pl.rid_shift);
+	if (pl.mi == 0) {
+		free(pl.fp_parts);
+		free(pl.rid_shift);
+		return -1;
+	}
+	for (i = n_split_idx - 1; i > 0; --i)
+		pl.rid_shift[i] = pl.rid_shift[i - 1];
+	for (pl.rid_shift[0] = 0, i = 1; i < n_split_idx; ++i)
+		pl.rid_shift[i] += pl.rid_shift[i - 1];
+	if (opt->flag & MM_F_OUT_SAM)
+		for (i = 0; i < (int32_t)pl.mi->n_seq; ++i)
+			printf("@SQ\tSN:%s\tLN:%d\n", pl.mi->seq[i].name, pl.mi->seq[i].len);
+
+	kt_pipeline(2, worker_pipeline, &pl, 3);
+
+	free(pl.str.s);
+	mm_idx_destroy(mi);
+	free(pl.rid_shift);
+	for (i = 0; i < n_split_idx; ++i)
+		fclose(pl.fp_parts[i]);
+	free(pl.fp_parts);
+	for (i = 0; i < pl.n_fp; ++i)
+		mm_bseq_close(pl.fp[i]);
+	free(pl.fp);
+	mm_split_rm_tmp(opt->split_prefix, n_split_idx);
+	return 0;
+}
@@ -29,6 +29,12 @@
 #define MM_F_REV_ONLY      0x200000
 #define MM_F_HEAP_SORT     0x400000
 #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
@@ -56,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)
@@ -113,8 +120,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;
@@ -126,9 +135,12 @@ typedef struct {
 	int pe_ori, pe_bonus;

 	float mid_occ_frac;  // only used by mm_mapopt_update(); see below
+	int32_t min_mid_occ;
 	int32_t mid_occ;     // ignore seeds with occurrences above this threshold
 	int32_t max_occ;
 	int mini_batch_size; // size of a batch of query bases to process in parallel
+
+	const char *split_prefix;
 } mm_mapopt_t;

 // index reader
@@ -213,6 +225,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
 *
@@ -261,6 +303,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
 *
@@ -297,6 +341,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);
@@ -1,4 +1,4 @@
-.TH minimap2 1 "24 February 2018" "minimap2-2.9 (r720)" "Bioinformatics tools"
+.TH minimap2 1 "11 October 2018" "minimap2-2.13 (r850)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -123,10 +123,27 @@ provided as the target sequences, options
 will be effectively overridden by the options stored in the index file.
 .SS Mapping options
 .TP 10
-.BI -f \ FLOAT
-Ignore top
+.BI -f \ FLOAT | INT1 [, INT2 ]
+If fraction, ignore top
 .I FLOAT
-fraction of most frequent minimizers [0.0002]
+fraction of most frequent minimizers [0.0002]. If integer,
+ignore minimizers occuring more than
+.I INT1
+times.
+.I INT2
+is only effective in the
+.B --sr
+or
+.B -xsr
+mode, which sets the threshold for a second round of seeding.
+.TP
+.BI --min-occ-floor \ INT
+Force minimap2 to always use k-mers occurring
+.I INT
+times or less [0]. In effect, the max occurrence threshold is set to
+the
+.RI max{ INT ,
+.BR -f }.
 .TP
 .BI -g \ INT
 Stop chain enlongation if there are no minimizers within
@@ -176,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 ,
@@ -209,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
@@ -222,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
@@ -231,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
@@ -305,6 +334,9 @@ no attempt to match GT-AG [n]
 .BI --end-bonus \ INT
 Score bonus when alignment extends to the end of the query sequence [0].
 .TP
+.BI --score-N \ INT
+Score of a mismatch involving ambiguous bases [1].
+.TP
 .BR --splice-flank = yes | no
 Assume the next base to a
 .B GT
@@ -333,6 +365,10 @@ the length of the terminal gap in the chain. This option is only effective
 with
 .BR --splice .
 It helps to avoid tiny terminal exons. [6]
+.TP
+.B --no-end-flt
+Don't filter seeds towards the ends of chains before performing base-level
+alignment.
 .SS Input/output options
 .TP 10
 .B -a
@@ -353,6 +389,9 @@ SAM read group line in a format like
 .B @RG\\\\tID:foo\\\\tSM:bar
 [].
 .TP
+.B -y
+Copy input FASTA/Q comments to output.
+.TP
 .B -c
 Generate CIGAR. In PAF, the CIGAR is written to the `cg' custom tag.
 .TP
@@ -371,6 +410,12 @@ is given,
 .I short
 is assumed. [none]
 .TP
+.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
@@ -433,18 +478,25 @@ is determined by the sequencing error mode.
 .B asm5
 Long assembly to reference mapping
 .RB ( -k19
-.B -w19 -A1 -B19 -O39,81 -E3,1 -s200
-.BR -z200 ).
+.B -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200
+.BR --min-occ-floor=100 ).
 Typically, the alignment will not extend to regions with 5% or higher sequence
 divergence. Only use this preset if the average divergence is far below 5%.
 .TP
 .B asm10
 Long assembly to reference mapping
 .RB ( -k19
-.B -w19 -A1 -B9 -O16,41 -E2,1 -s200
-.BR -z200 ).
+.B -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200
+.BR --min-occ-floor=100 ).
 Up to 10% sequence divergence.
 .TP
+.B asm20
+Long assembly to reference mapping
+.RB ( -k19
+.B -w10 -A1 -B6 -O6,26 -E2,1 -s200 -z200
+.BR --min-occ-floor=100 ).
+Up to 20% sequence divergence.
+.TP
 .B ava-pb
 PacBio all-vs-all overlap mapping
 .RB ( -Hk19
@@ -454,7 +506,7 @@ PacBio all-vs-all overlap mapping
 .B ava-ont
 Oxford Nanopore all-vs-all overlap mapping
 .RB ( -k15
-.B -Xw5 -m100 -g10000 --max-chain-skip
+.B -Xw5 -m100 -g10000 -r2000 --max-chain-skip
 .BR 25 ).
 Similarly, the major difference from
 .B ava-pb
@@ -535,12 +587,14 @@ cm	i	Number of minimizers on the chain
 s1	i	Chaining score
 s2	i	Chaining score of the best secondary chain
 NM	i	Total number of mismatches and gaps in the alignment
+MD	Z	To generate the ref sequence in the alignment
 AS	i	DP alignment score
 ms	i	DP score of the max scoring segment in the alignment
 nn	i	Number of ambiguous bases in the alignment
 ts	A	Transcript strand (splice mode only)
 cg	Z	CIGAR string (only in PAF)
 cs	Z	Difference string
+dv	f	Approximate per-base sequence divergence
 .TE

 .PP
@@ -1,3 +1,4 @@
+#include <stdlib.h>
 #include "mmpriv.h"

 int mm_verbose = 1;
@@ -120,6 +121,36 @@ double realtime(void)
 	return tp.tv_sec + tp.tv_usec * 1e-6;
 }

+void mm_err_puts(const char *str)
+{
+	int ret;
+	ret = puts(str);
+	if (ret == EOF) {
+		fprintf(stderr, "[ERROR] failed to write the results\n");
+		exit(EXIT_FAILURE);
+	}
+}
+
+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)
@@ -1,6 +1,6 @@
 #!/usr/bin/env k8

-var paftools_version = 'r713';
+var paftools_version = '2.13-r850';

 /*****************************
 ***** Library functions *****
@@ -131,6 +131,40 @@ Interval.find_ovlp = function(a, st, en)
 * Reverse and reverse complement *
 **********************************/

+function fasta_read(fn)
+{
+	var h = {}, gt = '>'.charCodeAt(0);
+	var file = fn == '-'? new File() : new File(fn);
+	var buf = new Bytes(), seq = null, name = null, seqlen = [];
+	while (file.readline(buf) >= 0) {
+		if (buf[0] == gt) {
+			if (seq != null && name != null) {
+				seqlen.push([name, seq.length]);
+				h[name] = seq;
+				name = seq = null;
+			}
+			var m, line = buf.toString();
+			if ((m = /^>(\S+)/.exec(line)) != null) {
+				name = m[1];
+				seq = new Bytes();
+			}
+		} else seq.set(buf);
+	}
+	if (seq != null && name != null) {
+		seqlen.push([name, seq.length]);
+		h[name] = seq;
+	}
+	buf.destroy();
+	file.close();
+	return [h, seqlen];
+}
+
+function fasta_free(fa)
+{
+	for (var name in fa)
+		fa[name].destroy();
+}
+
 Bytes.prototype.reverse = function()
 {
 	for (var i = 0; i < this.length>>1; ++i) {
@@ -305,14 +339,19 @@ function paf_liftover(args)
 // variant calling
 function paf_call(args)
 {
-	var re_cs = /([:=*+-])(\d+|[A-Za-z]+)/g;
-	var c, min_cov_len = 10000, min_var_len = 50000, gap_thres = 50, min_mapq = 5;
-	while ((c = getopt(args, "l:L:g:q:B:")) != null) {
+	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, gap_thres_long = 1000, min_mapq = 5;
+	var fa_tmp = null, fa, fa_lens, is_vcf = false, sample_name = "sample";
+	while ((c = getopt(args, "l:L:g:q:B:f:s:")) != null) {
 		if (c == 'l') min_cov_len = parseInt(getopt.arg);
 		else if (c == 'L') min_var_len = parseInt(getopt.arg);
 		else if (c == 'g') gap_thres = parseInt(getopt.arg);
+		else if (c == 'G') gap_thres_long = parseInt(getopt.arg);
 		else if (c == 'q') min_mapq = parseInt(getopt.arg);
+		else if (c == 'f') fa_tmp = fasta_read(getopt.arg, fa_lens);
+		else if (c == 's') sample_name = getopt.arg;
 	}
+	if (fa_tmp != null) fa = fa_tmp[0], fa_lens = fa_tmp[1], is_vcf = true;

 	if (args.length == getopt.ind) {
 		print("Usage: sort -k6,6 -k8,8n <with-cs.paf> | paftools.js call [options] -");
@@ -321,12 +360,35 @@ function paf_call(args)
 		print("  -L INT    min alignment length to call variants ["+min_var_len+"]");
 		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)
+	{
+		var v = null;
+		if (o[3] != 1) return; // coverage is one; skip
+		if (o[5] == '-' && o[6] == '-') return;
+		if (o[5] != '-' && o[6] != '-') { // snp
+			v = [o[0], o[1] + 1, '.', o[5].toUpperCase(), o[6].toUpperCase()];
+		} else if (o[1] > 0) { // shouldn't happen in theory
+			if (fa[o[0]] == null) throw Error('sequence "' + o[0] + '" is absent from the reference FASTA');
+			if (o[1] >= fa[o[0]].length) throw Error('position ' + o[1] + ' exceeds the length of sequence "' + o[0] + '"');
+			var ref = String.fromCharCode(fa[o[0]][o[1]-1]).toUpperCase();
+			if (o[5] == '-') // insertion
+				v = [o[0], o[1], '.', ref, ref + o[6].toUpperCase()];
+			else // deletion
+				v = [o[0], o[1], '.', ref + o[5].toUpperCase(), ref];
+		}
+		v.push(o[4], '.', 'QNAME=' + o[7] + ';QSTART=' + (o[8]+1) + ';QSTRAND=' + (rev? '-' : '+'), 'GT', '1/1');
+		if (v == null) throw Error("unexpected variant: [" + o.join(",") + "]");
+		print(v.join("\t"));
+	}

 	function count_var(o)
 	{
@@ -337,55 +399,77 @@ 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];
+			var s = (o[5] + o[6]).toLowerCase();
 			if (s == 'ag' || s == 'ga' || s == 'ct' || s == 'tc')
 				++n_sub[1];
 			else ++n_sub[2];
 		}
 	}

+	if (is_vcf) {
+		print('##fileformat=VCFv4.1');
+		for (var i = 0; i < fa_lens.length; ++i)
+			print('##contig=<ID=' + fa_lens[i][0] + ',length=' + fa_lens[i][1] + '>');
+		print('##INFO=<ID=QNAME,Number=1,Type=String,Description="Query name">');
+		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_name);
+	}
+
 	var a = [], out = [];
 	var c1_ctg = null, c1_start = 0, c1_end = 0, c1_counted = false, c1_len = 0;
 	while (file.readline(buf) >= 0) {
 		var line = buf.toString();
-		if (!/\ts2:i:/.test(line)) continue; // skip secondary alignments
 		var m, t = line.split("\t", 12);
 		for (var i = 6; i <= 11; ++i)
 			t[i] = parseInt(t[i]);
 		if (t[10] < min_cov_len || t[11] < min_mapq) continue;
-		print(t[0], t[7], t[8], c1_start, c1_end);
+		//print(t[0], t[7], t[8], c1_start, c1_end);
 		for (var i = 1; i <= 3; ++i)
 			t[i] = parseInt(t[i]);
 		var ctg = t[5], x = t[7], end = t[8];
 		var query = t[0], rev = (t[4] == '-'), y = rev? t[3] : t[2];
+		// collect tags
+		var cs = null, tp = null, have_s1 = false, have_s2 = false;
+		while ((m = re_tag.exec(line)) != null) {
+			if (m[1] == 'cs:Z') cs = m[2];
+			else if (m[1] == 'tp:A') tp = m[2];
+			else if (m[1] == 's1:i') have_s1 = true;
+			else if (m[1] == 's2:i') have_s2 = true;
+		}
+		if (have_s1 && !have_s2) continue;
+		if (tp != null && (tp == 'S' || tp == 'i')) continue;
 		// compute regions covered by 1 contig
 		if (ctg != c1_ctg || x >= c1_end) {
 			if (c1_counted && c1_end > c1_start) {
 				c1_len += c1_end - c1_start;
-				print('R', c1_ctg, c1_start, c1_end);
+				if (!is_vcf) print('R', c1_ctg, c1_start, c1_end);
 			}
 			c1_ctg = ctg, c1_start = x, c1_end = end;
 			c1_counted = (t[10] >= min_var_len);
 		} else if (end > c1_end) { // overlap
 			if (c1_counted && x > c1_start) {
 				c1_len += x - c1_start;
-				print('R', c1_ctg, c1_start, x);
+				if (!is_vcf) print('R', c1_ctg, c1_start, x);
 			}
 			c1_start = c1_end, c1_end = end;
 			c1_counted = (t[10] >= min_var_len);
 		} else if (end > c1_start) { // contained
 			if (c1_counted && x > c1_start) {
 				c1_len += x - c1_start;
-				print('R', c1_ctg, c1_start, x);
+				if (!is_vcf) print('R', c1_ctg, c1_start, x);
 			}
 			c1_start = end;
 		} // else, the alignment precedes the cov1 region; do nothing
@@ -393,7 +477,8 @@ function paf_call(args)
 		while (out.length) {
 			if (out[0][0] != ctg || out[0][2] <= x) {
 				count_var(out[0]);
-				print('V', out[0].join("\t"));
+				if (is_vcf) print_vcf(out[0], fa);
+				else print('V', out[0].join("\t"));
 				out.shift();
 			} else break;
 		}
@@ -404,20 +489,19 @@ 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
 		if (t[10] >= min_var_len) {
-			if ((m = /\tcs:Z:(\S+)/.exec(line)) == null) continue; // no cs tag
-			var cs = m[1];
+			if (cs == null) continue; // no cs tag
 			var blen = 0, n_diff = 0;
 			tot_len += t[10];
 			while ((m = re_cs.exec(cs)) != null) {
 				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]);
@@ -450,11 +534,12 @@ function paf_call(args)
 	}
 	if (c1_counted && c1_end > c1_start) {
 		c1_len += c1_end - c1_start;
-		print('R', c1_ctg, c1_start, c1_end);
+		if (!is_vcf) print('R', c1_ctg, c1_start, c1_end);
 	}
 	while (out.length) {
 		count_var(out[0]);
-		print('V', out[0].join("\t"));
+		if (is_vcf) print_vcf(out[0], fa);
+		else print('V', out[0].join("\t"));
 		out.shift();
 	}

@@ -467,11 +552,194 @@ function paf_call(args)
 	warn(n_ins[1] + " 2bp insertions");
 	warn(n_del[2] + " [3,"+gap_thres+") deletions");
 	warn(n_ins[2] + " [3,"+gap_thres+") insertions");
-	warn(n_del[3] + " >="+gap_thres+" deletions");
-	warn(n_ins[3] + " >="+gap_thres+" insertions");
+	warn(n_del[3] + " ["+gap_thres+","+gap_thres_long+") deletions");
+	warn(n_ins[3] + " ["+gap_thres+","+gap_thres_long+") insertions");
+	warn(n_del[4] + " >=" + gap_thres_long + " deletions");
+	warn(n_ins[4] + " >=" + gap_thres_long + " insertions");

 	buf.destroy();
 	file.close();
+	if (fa != null) fasta_free(fa);
+}
+
+function paf_asmstat(args)
+{
+	var c, min_seg_len = 10000, max_diff = 0.01;
+	while ((c = getopt(args, "l:d:")) != null) {
+		if (c == 'l') min_seg_len = parseInt(getopt.arg);
+		else if (c == 'd') max_diff = parseFloat(getopt.arg);
+	}
+	if (getopt.ind == args.length) {
+		print("Usage: paftools.js asmstat [options] <ref.fa.fai> <asm1.paf> [...]");
+		print("Options:");
+		print("  -l INT     min alignment block length [" + min_seg_len + "]");
+		print("  -d FLOAT   max gap-compressed sequence divergence [" + max_diff + "]");
+		exit(1);
+	}
+
+	var file, buf = new Bytes();
+
+	var ref_len = 0;
+	file = new File(args[getopt.ind]);
+	while (file.readline(buf) >= 0) {
+		var t = buf.toString().split("\t");
+		ref_len += parseInt(t[1]);
+	}
+	file.close();
+
+	function process_query(qblocks, qblock_len, bp) {
+		qblocks.sort(function(a,b) { return a[0]-b[0]; });
+		var last_k = null, last_blen = null, st = -1, en = -1, qcov = 0;
+		for (var k = 0; k < qblocks.length; ++k) {
+			var blen = qblocks[k][1] - qblocks[k][0];
+			if (k > 0 && qblocks[k][0] < qblocks[k-1][1]) {
+				if (qblocks[k][1] < qblocks[k-1][1]) continue;
+				blen = qblocks[k][1] - qblocks[k-1][1];
+			}
+			qblock_len.push(blen);
+			if (qblocks[k][0] > en) {
+				qcov += en - st;
+				st = qblocks[k][0];
+				en = qblocks[k][1];
+			} else en = en > qblocks[k][1]? en : qblocks[k][1];
+			if (last_k != null) {
+				var gap = 1000000000;
+				if (qblocks[k][2] == qblocks[last_k][2] && qblocks[k][3] == qblocks[last_k][3]) { // same chr and strand
+					var g1 = qblocks[k][0] - qblocks[last_k][1];
+					var g2 = qblocks[k][2] == '+'? qblocks[k][4] - qblocks[last_k][5] : qblocks[last_k][4] - qblocks[k][5];
+					gap = g1 > g2? g1 - g2 : g2 - g1;
+				}
+				var min = blen < last_blen? blen : last_blen;
+				var flank = k == 0? min : blen;
+				bp.push([flank, gap]);
+			}
+			last_k = k, last_blen = blen;
+		}
+		qcov += en - st;
+		return qcov;
+	}
+
+	function N50(lens, tot, quantile) {
+		lens.sort(function(a,b) { return b - a; });
+		if (tot == null) {
+			tot = 0;
+			for (var k = 0; k < lens.length; ++k)
+				tot += lens[k];
+		}
+		var sum = 0;
+		for (var k = 0; k < lens.length; ++k) {
+			if (sum <= quantile * tot && sum + lens[k] > quantile * tot)
+				return lens[k];
+			sum += lens[k];
+		}
+	}
+
+	function count_bp(bp, min_blen, min_gap) {
+		var n_bp = 0;
+		for (var k = 0; k < bp.length; ++k)
+			if (bp[k][0] >= min_blen && bp[k][1] >= min_gap)
+				++n_bp;
+		return n_bp;
+	}
+
+	function compute_diff(cigar, NM) {
+		var m, re = /(\d+)([MID])/g;
+		var n_M = 0, n_gapo = 0, n_gaps = 0;
+		while ((m = re.exec(cigar)) != null) {
+			var len = parseInt(m[1]);
+			if (m[2] == 'M') n_M += len;
+			else ++n_gapo, n_gaps += len;
+		}
+		if (NM < n_gaps) throw Error('NM is smaller the number of gaps');
+		return (NM - n_gaps + n_gapo) / (n_M + n_gapo);
+	}
+
+	var labels = ['Length', 'NG50', 'Coverage', 'Qcov', 'NGA50', '#breaks', 'bp(' + min_seg_len + ',0)', 'bp(' + min_seg_len + ',10k)'];
+	var rst = [];
+	for (var i = 0; i < labels.length; ++i)
+		rst[i] = [];
+
+	var n_asm = args.length - (getopt.ind + 1);
+	var header = ["Metric"];
+	for (var i = 0; i < n_asm; ++i) {
+		var n_breaks = 0, qcov = 0;
+		var fn = args[getopt.ind + 1 + i];
+		header.push(fn.replace(/.paf(.gz)?$/, ""));
+		var ref_blocks = [], qblock_len = [], qblocks = [], bp = [];
+		var query = {};
+		var last_qname = null;
+		file = new File(fn);
+		while (file.readline(buf) >= 0) {
+			var m, line = buf.toString();
+			var t = line.split("\t");
+			t[1] = parseInt(t[1]);
+			if (t.length >= 2) query[t[0]] = t[1];
+			if (t.length < 9) continue;
+			if (!/\ttp:A:[PI]/.test(line)) continue;
+			if ((m = /\tcg:Z:(\S+)/.exec(line)) == null) continue;
+			var cigar = m[1];
+			if ((m = /\tNM:i:(\d+)/.exec(line)) == null) continue;
+			var NM = parseInt(m[1]);
+			var diff = compute_diff(cigar, NM);
+			t[2] = parseInt(t[2]);
+			t[3] = parseInt(t[3]);
+			t[7] = parseInt(t[7]);
+			t[8] = parseInt(t[8]);
+			if (t[0] == last_qname) ++n_breaks;
+			if (diff > max_diff) continue;
+			if (t[3] - t[2] < min_seg_len) continue;
+			if (t[0] != last_qname) {
+				if (last_qname != null)
+					qcov += process_query(qblocks, qblock_len, bp);
+				qblocks = [];
+				last_qname = t[0];
+			}
+			ref_blocks.push([t[5], t[7], t[8]]);
+			qblocks.push([t[2], t[3], t[4], t[5], t[7], t[8]]);
+		}
+		if (last_qname != null)
+			qcov += process_query(qblocks, qblock_len, bp);
+		file.close();
+
+		// compute NG50
+		var asm_len = 0, asm_lens = []
+		for (var ctg in query) {
+			asm_len += query[ctg];
+			asm_lens.push(query[ctg]);
+		}
+		rst[0][i] = asm_len;
+		rst[1][i] = N50(asm_lens, ref_len, 0.5);
+
+		// compute coverage
+		var l_cov = 0;
+		ref_blocks.sort(function(a, b) { return a[0] > b[0]? 1 : a[0] < b[0]? -1 : a[1] - b[1]; });
+		var last_ref = null, st = -1, en = -1;
+		for (var j = 0; j < ref_blocks.length; ++j) {
+			if (ref_blocks[j][0] != last_ref || ref_blocks[j][1] > en) {
+				l_cov += en - st;
+				last_ref = ref_blocks[j][0];
+				st = ref_blocks[j][1];
+				en = ref_blocks[j][2];
+			} else en = en > ref_blocks[j][2]? en : ref_blocks[j][2];
+		}
+		l_cov += en - st;
+		rst[2][i] = (100.0 * (l_cov / ref_len)).toFixed(2) + '%';
+		rst[3][i] = (100.0 * (qcov  / asm_len)).toFixed(2) + '%';
+
+		// compute NGA50
+		rst[4][i] = N50(qblock_len, ref_len, 0.5);
+
+		// compute break points
+		rst[5][i] = n_breaks;
+		rst[6][i] = count_bp(bp, 500, 0);
+		rst[7][i] = count_bp(bp, 500, 10000);
+	}
+
+	print(header.join("\t"));
+	for (var i = 0; i < labels.length; ++i)
+		print(labels[i], rst[i].join("\t"));
+
+	buf.destroy();
 }

 function paf_stat(args)
@@ -595,8 +863,10 @@ function paf_stat(args)
 			last_qlen = ori_qlen;
 		}
 	}
-	l_tot += last_qlen;
-	l_cov += cov_len(regs);
+	if (regs.length) {
+		l_tot += last_qlen;
+		l_cov += cov_len(regs);
+	}

 	file.close();
 	buf.destroy();
@@ -912,14 +1182,15 @@ function paf_view(args)

 function paf_gff2bed(args)
 {
-	var c, fn_ucsc_fai = null, is_short = false;
-	while ((c = getopt(args, "u:s")) != null) {
+	var c, fn_ucsc_fai = null, is_short = false, keep_gff = false;
+	while ((c = getopt(args, "u:sg")) != null) {
 		if (c == 'u') fn_ucsc_fai = getopt.arg;
 		else if (c == 's') is_short = true;
+		else if (c == 'g') keep_gff = true;
 	}

 	if (getopt.ind == args.length) {
-		print("Usage: paftools.js gff2bed [-u ucsc-genome.fa.fai] <in.gff>");
+		print("Usage: paftools.js gff2bed [-g] [-u ucsc-genome.fa.fai] <in.gff>");
 		exit(1);
 	}

@@ -980,6 +1251,12 @@ function paf_gff2bed(args)
 	var exons = [], cds_st = 1<<30, cds_en = 0, last_id = null;
 	while (file.readline(buf) >= 0) {
 		var t = buf.toString().split("\t");
+		if (keep_gff) {
+			if (t[0].charAt(0) != '#' && ens2ucsc[t[0]] != null)
+				t[0] = ens2ucsc[t[0]];
+			print(t.join("\t"));
+			continue;
+		}
 		if (t[0].charAt(0) == '#') continue;
 		if (t[2] != "CDS" && t[2] != "exon") continue;
 		t[3] = parseInt(t[3]) - 1;
@@ -1028,15 +1305,19 @@ function paf_gff2bed(args)

 function paf_sam2paf(args)
 {
-	var c, pri_only = false;
+	var c, pri_only = false, use_eq = false;
 	while ((c = getopt(args, "p")) != null)
 		if (c == 'p') pri_only = true;
+	if (args.length == getopt.ind) {
+		print("Usage: paftools.js sam2paf [-p] <in.sam>");
+		exit(1);
+	}

-	var file = args.length == getopt.ind || args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]);
+	var file = args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]);
 	var buf = new Bytes();
-	var re = /(\d+)([MIDSHNX=])/g;
+	var re = /(\d+)([MIDSHNX=])/g, re_MD = /(\d+)|(\^[A-Za-z]+)|([A-Za-z])/g, re_tag = /\t(\S\S:[AZif]):(\S+)/g;

-	var len = {}, lineno = 0;
+	var ctg_len = {}, lineno = 0;
 	while (file.readline(buf) >= 0) {
 		var m, n_cigar = 0, line = buf.toString();
 		++lineno;
@@ -1044,37 +1325,52 @@ function paf_sam2paf(args)
 			if (/^@SQ/.test(line)) {
 				var name = (m = /\tSN:(\S+)/.exec(line)) != null? m[1] : null;
 				var l = (m = /\tLN:(\d+)/.exec(line)) != null? parseInt(m[1]) : null;
-				if (name != null && l != null) len[name] = l;
+				if (name != null && l != null) ctg_len[name] = l;
 			}
 			continue;
 		}
-		var t = line.split("\t");
+		var t = line.split("\t", 11);
 		var flag = parseInt(t[1]);
-		if (t[9] != '*' && t[10] != '*' && t[9].length != t[10].length) throw Error("ERROR at line " + lineno + ": inconsistent SEQ and QUAL lengths - " + t[9].length + " != " + t[10].length);
-		if (t[2] == '*' || (flag&4)) continue;
+		if (t[9] != '*' && t[10] != '*' && t[9].length != t[10].length)
+			throw Error("at line " + lineno + ": inconsistent SEQ and QUAL lengths - " + t[9].length + " != " + t[10].length);
+		if (t[2] == '*' || (flag&4) || t[5] == '*') continue;
 		if (pri_only && (flag&0x100)) continue;
-		var tlen = len[t[2]];
-		if (tlen == null) throw Error("ERROR at line " + lineno + ": can't find the length of contig " + t[2]);
-		var nn = (m = /\tnn:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : 0;
-		var NM = (m = /\tNM:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : null;
-		var have_NM = NM == null? false : true;
-		NM += nn;
-		var clip = [0, 0], I = [0, 0], D = [0, 0], M = 0, N = 0, ql = 0, tl = 0, mm = 0, ext_cigar = false;
-		while ((m = re.exec(t[5])) != null) {
-			var l = parseInt(m[1]);
-			if (m[2] == 'M') M += l, ql += l, tl += l, ext_cigar = false;
-			else if (m[2] == 'I') ++I[0], I[1] += l, ql += l;
-			else if (m[2] == 'D') ++D[0], D[1] += l, tl += l;
-			else if (m[2] == 'N') N += l, tl += l;
-			else if (m[2] == 'S') clip[M == 0? 0 : 1] = l, ql += l;
-			else if (m[2] == 'H') clip[M == 0? 0 : 1] = l;
-			else if (m[2] == '=') M += l, ql += l, tl += l, ext_cigar = true;
-			else if (m[2] == 'X') M += l, ql += l, tl += l, mm += l, ext_cigar = true;
-			++n_cigar;
+		var tlen = ctg_len[t[2]];
+		if (tlen == null) throw Error("at line " + lineno + ": can't find the length of contig " + t[2]);
+		// find tags
+		var nn = 0, NM = null, MD = null, md_list = [];
+		while ((m = re_tag.exec(line)) != null) {
+			if (m[1] == "NM:i") NM = parseInt(m[2]);
+			else if (m[1] == "nn:i") nn = parseInt(m[2]);
+			else if (m[1] == "MD:Z") MD = m[2];
 		}
+		if (t[9] == '*') MD = null;
+		// infer various lengths from CIGAR
+		var clip = [0, 0], soft_clip = 0, I = [0, 0], D = [0, 0], M = 0, N = 0, mm = 0, have_M = false, have_ext = false, cigar = [];
+		while ((m = re.exec(t[5])) != null) {
+			var l = parseInt(m[1]), op = m[2];
+			if (op == 'M') M += l, have_M = true;
+			else if (op == 'I') ++I[0], I[1] += l;
+			else if (op == 'D') ++D[0], D[1] += l;
+			else if (op == 'N') N += l;
+			else if (op == 'S') clip[n_cigar == 0? 0 : 1] = l, soft_clip += l;
+			else if (op == 'H') clip[n_cigar == 0? 0 : 1] = l;
+			else if (op == '=') M += l, have_ext = true, op = 'M';
+			else if (op == 'X') M += l, mm += l, have_ext = true, op = 'M';
+			++n_cigar;
+			if (MD != null && op != 'H') {
+				if (cigar.length > 0 && cigar[cigar.length-1][1] == op)
+					cigar[cigar.length-1][0] += l;
+				else cigar.push([l, op]);
+			}
+		}
+		var ql = M + I[1] + soft_clip;
+		var tl = M + D[1] + N;
+		var ts = parseInt(t[3]) - 1, te = ts + tl;
+		// checking coordinate and length consistencies
 		if (n_cigar > 65535)
 			warn("WARNING at line " + lineno + ": " + n_cigar + " CIGAR operations");
-		if (tl + parseInt(t[3]) - 1 > tlen) {
+		if (te > tlen) {
 			warn("WARNING at line " + lineno + ": alignment end position larger than ref length; skipped");
 			continue;
 		}
@@ -1082,24 +1378,78 @@ function paf_sam2paf(args)
 			warn("WARNING at line " + lineno + ": SEQ length inconsistent with CIGAR (" + t[9].length + " != " + ql + "); skipped");
 			continue;
 		}
-		if (!have_NM || ext_cigar) NM = I[1] + D[1] + mm;
-		if (NM < I[1] + D[1] + mm) {
-			warn("WARNING at line " + lineno + ": NM is less than the total number of gaps (" + NM + " < " + (I[1]+D[1]+mm) + ")");
-			NM = I[1] + D[1] + mm;
+		// parse MD
+		var cs = [];
+		if (MD != null) {
+			var k = 0, cx = 0, cy = 0, mx = 0, my = 0;
+			while ((m = re_MD.exec(MD)) != null) {
+				if (m[2] != null) { // deletion from the reference
+					var len = m[2].length - 1;
+					cs.push('-', m[2].substr(1));
+					mx += len, cx += len, ++k;
+				} else { // copy or mismatch
+					var ml = m[1] != null? parseInt(m[1]) : 1;
+					while (k < cigar.length && cigar[k][1] != 'D') {
+						var cl = cigar[k][0], op = cigar[k][1];
+						if (op == 'M') {
+							if (my + ml < cy + cl) {
+								if (ml > 0) {
+									if (m[3] != null) cs.push('*', m[3], t[9][my]);
+									else cs.push(':', ml);
+								}
+								mx += ml, my += ml, ml = 0;
+								break;
+							} else {
+								var dl = cy + cl - my;
+								cs.push(':', dl);
+								cx += cl, cy += cl, ++k;
+								mx += dl, my += dl, ml -= dl;
+							}
+						} else if (op == 'I') {
+							cs.push('+', t[9].substr(cy, cl));
+							cy += cl, my += cl, ++k;
+						} else if (op == 'S') {
+							cy += cl, my += cl, ++k;
+						} else throw Error("at line " + lineno + ": inconsistent MD tag");
+					}
+					if (ml != 0) throw Error("at line " + lineno + ": inconsistent MD tag");
+				}
+			}
+			if (cx != mx || cy != my) throw Error("at line " + lineno + ": inconsistent MD tag");
 		}
-		var extra = ["mm:i:"+(NM-I[1]-D[1]), "io:i:"+I[0], "in:i:"+I[1], "do:i:"+D[0], "dn:i:"+D[1]];
-		var match = M - (NM - I[1] - D[1]);
+		// compute matching length, block length and calibrate NM
+		if (have_ext && !have_M) { // extended CIGAR
+			if (NM != null && NM != I[1] + D[1] + mm)
+				warn("WARNING at line " + lineno + ": NM is different from sum of gaps and mismatches");
+			NM = I[1] + D[1] + mm;
+		} else if (NM != null) { // standard CIGAR; NM present
+			if (NM < I[1] + D[1]) {
+				warn("WARNING at line " + lineno + ": NM is less than the total number of gaps (" + NM + " < " + (I[1]+D[1]) + ")");
+				NM = I[1] + D[1];
+			}
+			mm = NM - (I[1] + D[1]);
+		} else { // no way to compute mm
+			warn("WARNING at line " + lineno + ": unable to find the number of mismatches; assuming zero");
+			mm = 0;
+		}
+		var mlen = M - mm;
 		var blen = M + I[1] + D[1];
+		// find query name, start and end
 		var qlen = M + I[1] + clip[0] + clip[1];
-		var qs, qe;
-		if (flag&16) qs = clip[1], qe = qlen - clip[0];
-		else qs = clip[0], qe = qlen - clip[1];
-		var ts = parseInt(t[3]) - 1, te = ts + M + D[1] + N;
-		var qname = t[0];
+		var qname = t[0], qs, qe;
 		if ((flag&1) && (flag&0x40)) qname += '/1';
 		if ((flag&1) && (flag&0x80)) qname += '/2';
-		var a = [qname, qlen, qs, qe, flag&16? '-' : '+', t[2], tlen, ts, te, match, blen, t[4]];
-		print(a.join("\t"), extra.join("\t"));
+		if (flag&16) qs = clip[1], qe = qlen - clip[0];
+		else qs = clip[0], qe = qlen - clip[1];
+		// optional tags
+		var type = flag&0x100? 'S' : 'P';
+		var tags = ["tp:A:" + type];
+		if (NM != null) tags.push("mm:i:"+mm);
+		tags.push("gn:i:"+(I[1]+D[1]), "go:i:"+(I[0]+D[0]), "cg:Z:" + t[5].replace(/\d+[SH]/g, ''));
+		if (cs.length > 0) tags.push("cs:Z:" + cs.join(""));
+		// print out
+		var a = [qname, qlen, qs, qe, flag&16? '-' : '+', t[2], tlen, ts, te, mlen, blen, t[4]];
+		print(a.join("\t"), tags.join("\t"));
 	}

 	buf.destroy();
@@ -1597,26 +1947,28 @@ function paf_pbsim2fq(args)

 function paf_junceval(args)
 {
-	var c, l_fuzzy = 0, print_ovlp = false, print_err_only = false, first_only = false;
-	while ((c = getopt(args, "l:ep")) != null) {
+	var c, l_fuzzy = 0, print_ovlp = false, print_err_only = false, first_only = false, chr_only = false;
+	while ((c = getopt(args, "l:epc")) != null) {
 		if (c == 'l') l_fuzzy = parseInt(getopt.arg);
 		else if (c == 'e') print_err_only = print_ovlp = true;
 		else if (c == 'p') print_ovlp = true;
+		else if (c == 'c') chr_only = true;
 	}

-	if (args.length - getopt.ind < 2) {
+	if (args.length - getopt.ind < 1) {
 		print("Usage: paftools.js junceval [options] <gene.gtf> <aln.sam>");
 		print("Options:");
 		print("  -l INT    tolerance of junction positions (0 for exact) [0]");
 		print("  -p        print overlapping introns");
 		print("  -e        print erroreous overlapping introns");
+		print("  -c        only consider alignments to /^(chr)?([0-9]+|X|Y)$/");
 		exit(1);
 	}

 	var file, buf = new Bytes();

 	var tr = {};
-	file = new File(args[getopt.ind]);
+	file = args[getopt.ind] == '-'? new File() : new File(args[getopt.ind]);
 	while (file.readline(buf) >= 0) {
 		var m, t = buf.toString().split("\t");
 		if (t[0].charAt(0) == '#') continue;
@@ -1661,13 +2013,14 @@ function paf_junceval(args)
 	var n_pri = 0, n_unmapped = 0, n_mapped = 0;
 	var n_sgl = 0, n_splice = 0, n_splice_hit = 0, n_splice_novel = 0;

-	file = new File(args[getopt.ind+1]);
+	file = getopt.ind+1 >= args.length || args[getopt.ind+1] == '-'? new File() : new File(args[getopt.ind+1]);
 	var last_qname = null;
 	var re_cigar = /(\d+)([MIDNSHX=])/g;
 	while (file.readline(buf) >= 0) {
 		var m, t = buf.toString().split("\t");

 		if (t[0].charAt(0) == '@') continue;
+		if (chr_only && !/^(chr)?([0-9]+|X|Y)$/.test(t[2])) continue;
 		var flag = parseInt(t[1]);
 		if (flag&0x100) continue;
 		if (first_only && last_qname == t[0]) continue;
@@ -1835,6 +2188,7 @@ function main(args)
 		print("  gff2bed    convert GTF/GFF3 to BED12");
 		print("");
 		print("  stat       collect basic mapping information in PAF/SAM");
+		print("  asmstat    collect basic assembly information");
 		print("  liftover   simplistic liftOver");
 		print("  call       call variants from asm-to-ref alignment with the cs tag");
 		print("  bedcov     compute the number of bases covered");
@@ -1855,6 +2209,7 @@ function main(args)
 	else if (cmd == 'splice2bed') paf_splice2bed(args);
 	else if (cmd == 'gff2bed') paf_gff2bed(args);
 	else if (cmd == 'stat') paf_stat(args);
+	else if (cmd == 'asmstat') paf_asmstat(args);
 	else if (cmd == 'liftover' || cmd == 'liftOver') paf_liftover(args);
 	else if (cmd == 'call') paf_call(args);
 	else if (cmd == 'mapeval') paf_mapeval(args);
@@ -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,6 +89,15 @@ 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
 }
 #endif
@@ -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,10 +49,12 @@ 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);
+	if (opt->mid_occ < opt->min_mid_occ)
+		opt->mid_occ = opt->min_mid_occ;
 	if (mm_verbose >= 3)
 		fprintf(stderr, "[M::%s::%.3f*%.2f] mid_occ = %d\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), opt->mid_occ);
 }
@@ -70,6 +74,7 @@ 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;
@@ -81,11 +86,19 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 	} else if (strcmp(preset, "asm5") == 0) {
 		io->flag = 0, io->k = 19, io->w = 19;
 		mo->a = 1, mo->b = 19, mo->q = 39, mo->q2 = 81, mo->e = 3, mo->e2 = 1, mo->zdrop = mo->zdrop_inv = 200;
+		mo->min_mid_occ = 100;
 		mo->min_dp_max = 200;
 		mo->best_n = 50;
 	} else if (strcmp(preset, "asm10") == 0) {
 		io->flag = 0, io->k = 19, io->w = 19;
 		mo->a = 1, mo->b = 9, mo->q = 16, mo->q2 = 41, mo->e = 2, mo->e2 = 1, mo->zdrop = mo->zdrop_inv = 200;
+		mo->min_mid_occ = 100;
+		mo->min_dp_max = 200;
+		mo->best_n = 50;
+	} else if (strcmp(preset, "asm20") == 0) {
+		io->flag = 0, io->k = 19, io->w = 10;
+		mo->a = 1, mo->b = 4, mo->q = 6, mo->q2 = 26, mo->e = 2, mo->e2 = 1, mo->zdrop = mo->zdrop_inv = 200;
+		mo->min_mid_occ = 100;
 		mo->min_dp_max = 200;
 		mo->best_n = 50;
 	} else if (strcmp(preset, "short") == 0 || strcmp(preset, "sr") == 0) {
@@ -119,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");
@@ -151,5 +174,10 @@ int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
 			fprintf(stderr, "[ERROR]\033[1;31m Z-drop should not be less than inversion-Z-drop\033[0m\n");
 		return -5;
 	}
+	if ((mo->flag & MM_F_NO_PRINT_2ND) && (mo->flag & MM_F_ALL_CHAINS)) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m -X/-P and --secondary=no can't be applied at the same time\033[0m\n");
+		return -5;
+	}
 	return 0;
 }
@@ -105,7 +105,7 @@ void mm_pair(void *km, int max_gap_ref, int pe_bonus, int sub_diff, int match_sc
 	max = -1;
 	max_idx[0] = max_idx[1] = -1;
 	last[0] = last[1] = -1;
-	kv_resize(uint64_t, km, sc, n);
+	kv_resize(uint64_t, km, sc, (size_t)n);
 	for (i = 0; i < n; ++i) {
 		if (a[i].key & 1) { // reverse first read or forward second read
 			mm_reg1_t *q, *r;
@@ -151,8 +151,8 @@ void mm_pair(void *km, int max_gap_ref, int pe_bonus, int sub_diff, int match_sc
 			}
 		}
 		mapq_pe = r[0]->mapq > r[1]->mapq? r[0]->mapq : r[1]->mapq;
-		for (i = 0; i < sc.n; ++i)
-			if ((sc.a[i]>>32) + sub_diff >= max>>32)
+		for (i = 0; i < (int)sc.n; ++i)
+			if ((sc.a[i]>>32) + sub_diff >= (uint64_t)max>>32)
 				++n_sub;
 		if (sc.n > 1) {
 			int mapq_pe_alt;
@@ -164,7 +164,7 @@ void mm_pair(void *km, int max_gap_ref, int pe_bonus, int sub_diff, int match_sc
 		if (sc.n == 1) {
 			if (r[0]->mapq < 2) r[0]->mapq = 2;
 			if (r[1]->mapq < 2) r[1]->mapq = 2;
-		} else if (max>>32 > sc.a[sc.n - 2]>>32) {
+		} else if ((uint64_t)max>>32 > sc.a[sc.n - 2]>>32) {
 			if (r[0]->mapq < 1) r[0]->mapq = 1;
 			if (r[1]->mapq < 1) r[1]->mapq = 1;
 		}
@@ -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:

@@ -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,15 +126,27 @@ 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[i].len || st >= en) return 0;
-	if (en < 0 || en > mi->seq[i].len)
-		en = mi->seq[i].len;
+	if ((uint32_t)st >= mi->seq[rid].len || st >= en) return 0;
+	if (en < 0 || (uint32_t)en > mi->seq[rid].len)
+		en = mi->seq[rid].len;
 	s = (char*)malloc(en - st + 1);
-	*len = mm_idx_getseq(mi, rid, st, en, s);
+	*len = mm_idx_getseq(mi, rid, st, en, (uint8_t*)s);
 	for (i = 0; i < *len; ++i)
 		s[i] = "ACGTN"[(uint8_t)s[i]];
 	s[*len] = 0;
 	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
@@ -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
@@ -34,9 +36,11 @@ cdef extern from "minimap.h":
 		float max_clip_ratio
 		int pe_ori, pe_bonus
 		float mid_occ_frac
+		int32_t min_mid_occ
 		int32_t mid_occ
 		int32_t max_occ
 		int mini_batch_size
+		const char *split_prefix

 	int mm_set_opt(char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 	int mm_verbose
@@ -58,7 +62,8 @@ cdef extern from "minimap.h":
 		mm_idx_seq_t *seq
 		uint32_t *S
 		mm_idx_bucket_t *B
-		void *km, *h
+		void *km
+		void *h

 	ctypedef struct mm_idx_reader_t:
 		pass
@@ -82,6 +87,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)
@@ -102,6 +110,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
@@ -3,6 +3,8 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys

+__version__ = '2.13'
+
 cmappy.mm_reset_timer()

 cdef class Alignment:
@@ -12,9 +14,9 @@ cdef class Alignment:
 	cdef int8_t _strand, _trans_strand
 	cdef uint8_t _mapq, _is_primary
 	cdef int _seg_id
-	cdef _ctg, _cigar # these are python objects
+	cdef _ctg, _cigar, _cs, _MD # these are python objects

-	def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand, seg_id):
+	def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand, seg_id, cs_str, MD_str):
 		self._ctg = ctg if isinstance(ctg, str) else ctg.decode()
 		self._ctg_len, self._r_st, self._r_en = cl, cs, ce
 		self._strand, self._q_st, self._q_en = strand, qs, qe
@@ -24,6 +26,8 @@ cdef class Alignment:
 		self._is_primary = is_primary
 		self._trans_strand = trans_strand
 		self._seg_id = seg_id
+		self._cs = cs_str
+		self._MD = MD_str

 	@property
 	def ctg(self): return self._ctg
@@ -70,6 +74,12 @@ cdef class Alignment:
 	@property
 	def read_num(self): return self._seg_id + 1

+	@property
+	def cs(self): return self._cs
+
+	@property
+	def MD(self): return self._MD
+
 	@property
 	def cigar_str(self):
 		return "".join(map(lambda x: str(x[0]) + 'MIDNSH'[x[1]], self._cigar))
@@ -83,8 +93,10 @@ cdef class Alignment:
 		if self._trans_strand > 0: ts = 'ts:A:+'
 		elif self._trans_strand < 0: ts = 'ts:A:-'
 		else: ts = 'ts:A:.'
-		return "\t".join([str(self._q_st), str(self._q_en), strand, self._ctg, str(self._ctg_len), str(self._r_st), str(self._r_en),
-				str(self._mlen), str(self._blen), str(self._mapq), tp, ts, "cg:Z:" + self.cigar_str])
+		a = [str(self._q_st), str(self._q_en), strand, self._ctg, str(self._ctg_len), str(self._r_st), str(self._r_en),
+			str(self._mlen), str(self._blen), str(self._mapq), tp, ts, "cg:Z:" + self.cigar_str]
+		if self._cs != "": a.append("cs:Z:" + self._cs)
+		return "\t".join(a)

 cdef class ThreadBuffer:
 	cdef cmappy.mm_tbuf_t *_b
@@ -100,7 +112,7 @@ cdef class Aligner:
 	cdef cmappy.mm_idxopt_t idx_opt
 	cdef cmappy.mm_mapopt_t map_opt

-	def __cinit__(self, fn_idx_in, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, best_n=None, n_threads=3, fn_idx_out=None):
+	def __cinit__(self, fn_idx_in=None, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, best_n=None, n_threads=3, fn_idx_out=None, max_frag_len=None, extra_flags=None, seq=None, scoring=None):
 		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)
-		if r is not NULL:
-			self._idx = cmappy.mm_idx_reader_read(r, n_threads) # NB: ONLY read the first part
-			cmappy.mm_idx_reader_close(r)
+			self._idx = cmappy.mappy_idx_seq(self.idx_opt.w, self.idx_opt.k, self.idx_opt.flag&1, self.idx_opt.bucket_bits, str.encode(seq), len(seq))
 			cmappy.mm_mapopt_update(&self.map_opt, self._idx)
-			cmappy.mm_idx_index_name(self._idx)
+			self.map_opt.mid_occ = 1000 # don't filter high-occ seeds

 	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 = []
-			for k in range(h.n_cigar32):
+			cigar, _cs, _MD = [], '', ''
+			for k in range(h.n_cigar32): # convert the 32-bit CIGAR encoding to Python array
 				c = h.cigar32[k]
 				cigar.append([c>>4, c&0xf])
-			yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.mlen, h.blen, h.NM, h.trans_strand, h.seg_id)
+			if cs or MD: # generate the cs and/or the MD tag, if requested
+				if cs:
+					l_cs_str = cmappy.mm_gen_cs(km, &cs_str, &m_cs_str, self._idx, &regs[i], _seq, 1)
+					_cs = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
+				if MD:
+					l_cs_str = cmappy.mm_gen_MD(km, &cs_str, &m_cs_str, self._idx, &regs[i], _seq)
+					_MD = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
+			yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.mlen, h.blen, h.NM, h.trans_strand, h.seg_id, _cs, _MD)
 			cmappy.mm_free_reg1(&regs[i])
 		free(regs)
+		free(cs_str)

 	def seq(self, str name, int start=0, int end=0x7fffffff):
 		cdef int l
@@ -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__":
@@ -70,10 +70,10 @@ void sdust_buf_destroy(sdust_buf_t *buf)
 static inline void shift_window(int t, kdq_t(int) *w, int T, int W, int *L, int *rw, int *rv, int *cw, int *cv)
 {
 	int s;
-	if (kdq_size(w) >= W - SD_WLEN + 1) { // TODO: is this right for SD_WLEN!=3?
+	if ((int)kdq_size(w) >= W - SD_WLEN + 1) { // TODO: is this right for SD_WLEN!=3?
 		s = *kdq_shift(int, w);
 		*rw -= --cw[s];
-		if (*L > kdq_size(w))
+		if (*L > (int)kdq_size(w))
 			--*L, *rv -= --cv[s];
 	}
 	kdq_push(int, w, t);
@@ -114,7 +114,7 @@ static void find_perfect(void *km, perf_intv_v *P, const kdq_t(int) *w, int T, i
 		r += c[t]++;
 		new_r = r, new_l = kdq_size(w) - i - 1;
 		if (new_r * 10 > T * new_l) {
-			for (j = 0; j < P->n && P->a[j].start >= i + start; ++j) { // find insertion position
+			for (j = 0; j < (int)P->n && P->a[j].start >= i + start; ++j) { // find insertion position
 				perf_intv_t *p = &P->a[j];
 				if (max_r == 0 || p->r * max_l > max_r * p->l)
 					max_r = p->r, max_l = p->l;
@@ -177,7 +177,7 @@ uint64_t *sdust(void *km, const uint8_t *seq, int l_seq, int T, int W, int *n)
 #ifdef _SDUST_MAIN
 #include <zlib.h>
 #include <stdio.h>
-#include "getopt.h"
+#include "ketopt.h"
 #include "kseq.h"
 KSEQ_INIT(gzFile, gzread)

@@ -186,16 +186,17 @@ int main(int argc, char *argv[])
 	gzFile fp;
 	kseq_t *ks;
 	int W = 64, T = 20, c;
+	ketopt_t o = KETOPT_INIT;

-	while ((c = getopt(argc, argv, "w:t:")) >= 0) {
-		if (c == 'w') W = atoi(optarg);
-		else if (c == 't') T = atoi(optarg);
+	while ((c = ketopt(&o, argc, argv, 1, "w:t:", 0)) >= 0) {
+		if (c == 'w') W = atoi(o.arg);
+		else if (c == 't') T = atoi(o.arg);
 	}
-	if (optind == argc) {
+	if (o.ind == argc) {
 		fprintf(stderr, "Usage: sdust [-w %d] [-t %d] <in.fa>\n", W, T);
 		return 1;
 	}
-	fp = strcmp(argv[optind], "-")? gzopen(argv[optind], "r") : gzdopen(fileno(stdin), "r");
+	fp = strcmp(argv[o.ind], "-")? gzopen(argv[o.ind], "r") : gzdopen(fileno(stdin), "r");
 	ks = kseq_init(fp);
 	while (kseq_read(ks) >= 0) {
 		uint64_t *r;
@@ -14,16 +14,26 @@ else: # with Cython
 	module_src = 'python/mappy.pyx'
 	cmdclass['build_ext'] = build_ext

-import sys
+import sys, platform
+
 sys.path.append('python')

+extra_compile_args = ['-DHAVE_KALLOC']
+include_dirs = ["."]
+
+if platform.machine() in ["aarch64", "arm64"]:
+	include_dirs.append("sse2neon/")
+	extra_compile_args.extend(['-ftree-vectorize', '-DKSW_SSE2_ONLY', '-D__SSE2__'])
+else:
+	extra_compile_args.append('-msse4.1') # WARNING: ancient x86_64 CPUs don't have SSE4
+
 def readme():
-    with open('python/README.rst') as f:
-        return f.read()
+	with open('python/README.rst') as f:
+		return f.read()

 setup(
 	name = 'mappy',
-	version = '2.9',
+	version = '2.13',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),
@@ -35,12 +45,12 @@ setup(
    ext_modules = [Extension('mappy',
 		sources = [module_src, 'align.c', 'bseq.c', 'chain.c', 'format.c', 'hit.c', 'index.c', 'pe.c', 'options.c',
 				   'ksw2_extd2_sse.c', 'ksw2_exts2_sse.c', 'ksw2_extz2_sse.c', 'ksw2_ll_sse.c',
-				   'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c', 'esterr.c'],
+				   'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c', 'esterr.c', 'splitidx.c'],
 		depends = ['minimap.h', 'bseq.h', 'kalloc.h', 'kdq.h', 'khash.h', 'kseq.h', 'ksort.h',
 				   'ksw2.h', 'kthread.h', 'kvec.h', 'mmpriv.h', 'sdust.h',
 				   'python/cmappy.h', 'python/cmappy.pxd'],
-		extra_compile_args = ['-DHAVE_KALLOC', '-msse4'], # WARNING: ancient x86_64 CPUs don't have SSE4
-		include_dirs = ['.'],
+		extra_compile_args = extra_compile_args,
+		include_dirs = include_dirs,
 		libraries = ['z', 'm', 'pthread'])],
 	classifiers = [
 		'Development Status :: 5 - Production/Stable',
@@ -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);
+}
@@ -1,4 +1,4 @@
->MT_orang
+>MT_orang co:Z:comment
 GTTTATGTAGCTTATTCTATCCAAAGCAATGCACTGAAAATGTCTCGACGGGCCCACACG
 CCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTGAGGTTACACAT
 GCAAGCATCCCCGCCCCAGTGAGTCGCCCTCCAAGTCACTCTGACTAAGAGGAGCAAGCA
@@ -61,13 +61,6 @@
 	Volume = {32},
 	Year = {2016}}

-@misc{Suzuki:2016,
-	title = {Fast and accurate alignment tool for PacBio and Nanopore long reads},
-	author = {Hajime Suzuki},
-	journal = {Unpublished},
-	howpublished = {\href{https://github.com/ocxtal/minialign}{https://github.com/ocxtal/minialign}},
-	year = {2016}}
-
@misc{Ruan:2016,
 	title = {Ultra-fast de novo assembler using long noisy reads},
 	author = {Jue Ruan},
@@ -172,14 +165,6 @@
 	Volume = {29},
 	Year = {2011}}

-@article {Suzuki130633,
-	author = {Suzuki, Hajime and Kasahara, Masahiro},
-	title = {Acceleration Of Nucleotide Semi-Global Alignment With Adaptive Banded Dynamic Programming},
-	year = {2017},
-	note = {doi:10.1101/130633},
-	publisher = {Cold Spring Harbor Labs Journals},
-	journal = {bioRxiv}}
-
@article{Gotoh:1982aa,
 	Author = {Gotoh, O},
 	Journal = {J Mol Biol},
@@ -337,3 +322,19 @@
 	Title = {{MUMmer4}: A fast and versatile genome alignment system},
 	Volume = {14},
 	Year = {2018}}
+
+@article{Li:2009ys,
+	Author = {Li, Heng and others},
+	Journal = {Bioinformatics},
+	Pages = {2078-9},
+	Title = {The {Sequence Alignment/Map format and SAMtools}},
+	Volume = {25},
+	Year = {2009}}
+
+@article{Suzuki:2018aa,
+	Author = {Suzuki, Hajime and Kasahara, Masahiro},
+	Journal = {BMC Bioinformatics},
+	Pages = {45},
+	Title = {Introducing difference recurrence relations for faster semi-global alignment of long sequences},
+	Volume = {19},
+	Year = {2018}}
@@ -19,7 +19,7 @@
 \begin{document}
 \firstpage{1}

-\title[Aligning nucleotide sequences with minimap2]{Minimap2: versatile pairwise alignment for nucleotide sequences}
+\title[Aligning nucleotide sequences with minimap2]{Minimap2: pairwise alignment for nucleotide sequences}
 \author[Li]{Heng Li}
 \address{Broad Institute, 415 Main Street, Cambridge, MA 02142, USA}

@@ -64,7 +64,7 @@ the thought that 10kb long sequences should be easier to map than 100bp reads
 because we can more effectively skip repetitive regions, which are often the
 bottleneck of short-read alignment. We confirmed our speculation by achieving
 approximate mapping 50 times faster than BWA-MEM~\citep{Li:2016aa}.
-\citet{Suzuki130633} extended our work with a fast and novel algorithm on
+\citet{Suzuki:2018aa} extended our work with a fast and novel algorithm on
 generating base-level alignment, which in turn inspired us to develop minimap2
 with added functionality.

@@ -88,7 +88,9 @@ the versatility of minimap2.

 Minimap2 follows a typical seed-chain-align procedure as is used by most
 full-genome aligners. It collects minimizers~\citep{Roberts:2004fv} of the
-reference sequences and indexes them in a hash table. Then for each query
+reference sequences and indexes them in a hash table, with the key being the
+hash of a minimizer and the value being a list of locations of the minimizer
+copies. Then for each query
 sequence, minimap2 takes query minimizers as \emph{seeds}, finds exact matches
 (i.e. \emph{anchors}) to the reference, and identifies sets of colinear anchors as
 \emph{chains}. If base-level alignment is requested, minimap2 applies dynamic
@@ -118,9 +120,12 @@ distance between two anchors is too large); otherwise
 \begin{equation}\label{eq:chain-gap}
 \beta(j,i)=\gamma_c\big((y_i-y_j)-(x_i-x_j)\big)
 \end{equation}
-In implementation, a gap of length $l\not=0$ costs
+In implementation, a gap of length $l$ costs
 \[
-\gamma_c(l)=0.01\cdot \bar{w}\cdot|l|+0.5\log_2|l|
+\gamma_c(l)=\left\{\begin{array}{ll}
+0.01\cdot \bar{w}\cdot|l|+0.5\log_2|l| & (l\not=0) \\
+0 & (l=0)
+\end{array}\right.
 \]
 where $\bar{w}$ is the average seed length. For $N$ anchors, directly computing all $f(\cdot)$ with
 Eq.~(\ref{eq:chain}) takes $O(N^2)$ time. Although theoretically faster
@@ -164,7 +169,7 @@ empirical formula:
 \[
 {\rm mapQ}=40\cdot (1-f_2/f_1)\cdot\min\{1,m/10\}\cdot\log f_1
 \]
-where $m$ is the number of anchors on the primary chain, $f_1$ is the chaining
+where $\log$ denotes natural logarithm, $m$ is the number of anchors on the primary chain, $f_1$ is the chaining
 score, and $f_2\le f_1$ is the score of the best chain that is secondary to the
 primary chain. Intuitively, a chain is assigned to a higher mapping quality if
 it is long and its best secondary chain is weak.
@@ -253,7 +258,7 @@ performance of minimap2. Traditional SSE implementations~\citep{Farrar:2007hs}
 based on Eq.~(\ref{eq:ae86}) can achieve 16-way parallelization for short
 sequences, but only 4-way parallelization when the peak alignment score reaches
 32767. Long sequence alignment may exceed this threshold. Inspired by
-\citet{Wu:1996aa} and the following work, \citet{Suzuki130633} proposed a
+\citet{Wu:1996aa} and the following work, \citet{Suzuki:2018aa} proposed a
 difference-based formulation that lifted this limitation.
 In case of 2-piece gap cost, define
 \[
@@ -320,7 +325,7 @@ y_{rt}&=&\max\{0,y_{r-1,t}+u_{r-1,t}-z_{rt}+q\}-q-e\\
 \end{equation*}
 In this formulation, cells with the same diagonal index $r$ are independent of
 each other. This allows us to fully vectorize the computation of all cells on
-the same anti-diagonal in one inner loop. It also simplifies banded alignment,
+the same anti-diagonal in one inner loop. It also simplifies banded alignment (500bp band width by default),
 which would be difficult with striped vectorization~\citep{Farrar:2007hs}.

 On the condition that $q+e<\tilde{q}+\tilde{e}$ and $e>\tilde{e}$, the initial
@@ -355,7 +360,7 @@ times as fast as Parasail's 4-way vectorization~\citep{Daily:2016aa}.  Without
 banding, our implementation is slower than Edlib~\citep{Sosic:2017aa}, but with
 a 1000bp band, it is considerably faster. When performing global alignment
 between anchors, we expect the alignment to stay close to the diagonal of the
-DP matrix. Banding is applicable most of time.
+DP matrix. Banding is applicable most of the time.

 \subsubsection{The Z-drop heuristic}

@@ -478,7 +483,7 @@ both C and Python. It is distributed under the MIT license, free to both
 commercial and academic uses. Minimap2 uses the same base algorithm for all
 applications, but it has to apply different sets of parameters depending on
 input data types. Similar to BWA-MEM, minimap2 introduces `presets' that
-modify multiple parameters with a simple invokation. Detailed settings
+modify multiple parameters with a simple invocation. Detailed settings
 and command-line options can be found in the minimap2 manpage. In addition to
 the applications evaluated in the following sections, minimap2 also retains
 minimap's functionality to find overlaps between long reads and to search
@@ -570,7 +575,7 @@ Peak RAM (GByte)          & 8.9      & 14.5     & 3.2      & 29.2\vspace{1em}\\
 \% approx. introns        & 91.8\%   & 96.9\%   & 92.5\%   & 82.4\% \\
 \botrule
 \end{tabular}
-}{Mouse reads (AC:SRR5286960; R9.4 chemistry) were mapped to the primary assembly of mouse
+}{Mouse cDNA reads (AC:SRR5286960; R9.4 chemistry) were mapped to the primary assembly of mouse
 genome GRCm38 with the following tools and command options: minimap2 (`-ax
 splice'); GMAP (`-n 0 --min-intronlength 30 --cross-species'); SpAln (`-Q7 -LS
 -S3'); STARlong (according to
@@ -579,7 +584,7 @@ compared to the EnsEMBL gene annotation, release 89. A predicted intron
 is \emph{novel} if it has no overlaps with any annotated introns. An intron
 is \emph{exact} if it is identical to an annotated intron. An intron is
 \emph{approximate} if both its 5'- and 3'-end are within 10bp around the ends
-of an annotated intron.}
+of an annotated intron. Chimeric alignments are defined in the SAM spec~\citep{Li:2009ys}.}
 \end{table}

 We next aligned real mouse reads~\citep{Byrne:2017aa} with GMAP~(v2017-06-20;
@@ -647,9 +652,9 @@ across the whole genome and have been \emph{de novo} assembled with SMRT reads
 to high quality. This allowed us to construct an independent truth variant
 dataset~\citep{Li223297} for
 ERR1341796. In this evaluation, minimap2 has higher SNP false negative rate
-(FNR; 2.5\% of minimap2 vs 2.2\% of BWA-MEM), but fewer false positive SNPs per
-million bases (FPPM; 3.0 vs 3.9), lower 2--50bp INDEL FNR (7.3\% vs 7.5\%) and
-similar INDEL FPPM (both 1.0). Minimap2 is broadly similar to BWA-MEM in the
+(FNR; 2.6\% of minimap2 vs 2.3\% of BWA-MEM), but fewer false positive SNPs per
+million bases (FPPM; 7.0 vs 8.8), similar INDEL FNR (11.2\% vs 11.3\%) and
+similar INDEL FPPM (6.4 vs 6.5). Minimap2 is broadly comparable to BWA-MEM in the
 context of small variant calling.

 \subsection{Aligning long-read assemblies}
@@ -687,7 +692,7 @@ involving $>$100kb introns, which was impractically slow ten years ago.  The
 minimap2 chaining algorithm is fast and highly accurate by itself.  In fact,
 chaining alone is more accurate than all the other long-read mappers in
 Fig.~\ref{fig:eval}a (data not shown). This accuracy helps to reduce downstream
-base-level alignment of candidate chains, which is still times slower than
+base-level alignment of candidate chains, which is still several times slower than
 chaining even with the Suzuki-Kasahara improvement. In addition, taking a
 general form, minimap2 chaining can be adapted to non-typical data types such as
 spliced reads and multiple reads per fragment. This gives us the opportunity to
@@ -712,6 +717,8 @@ Schatz, P. Rescheneder and F. Sedlazeck for pointing out the limitation of
 BWA-MEM. We are also grateful to minimap2 users who have greatly helped to
 suggest features and to fix various issues.

+\paragraph{Funding\textcolon} NHGRI 1R01HG010040-01
+
 \bibliography{minimap2}

 \end{document}
Author	SHA1	Message	Date
Heng Li	5d5d392c02	Release minimap2-2.13 (r850)	2018-10-11 13:18:31 -04:00
Heng Li	170863e553	r849: option -P doesn't work I don't know why I haven't found it at the beginning.	2018-10-04 16:11:59 -04:00
Heng Li	97f97306a4	r847: guard against -N0	2018-09-27 15:13:44 -04:00
Heng Li	1077b7ddc8	r846: added --hard-mask-level for #244	2018-09-27 14:46:26 -04:00
Heng Li	c57b59f02f	r845: log peak memory	2018-09-23 20:27:49 -04:00
Torsten Seemann	34be359e25	Add -s sample option for VCF output	2018-09-19 18:24:47 -04:00
Heng Li	8b12da8b0f	removed a useless dependency (not in repo)	2018-09-17 13:59:25 -04:00
Heng Li	c63a33904f	r836: fixed an integer overflow Forgot this one.	2018-09-14 23:29:31 -04:00
Heng Li	70b0fede64	r835: improved help message. Resolved #232	2018-09-14 22:29:25 -04:00
Heng Li	0b681e51e7	Merge remote-tracking branch 'remotes/origin/master'	2018-09-14 22:22:20 -04:00
Heng Li	7d80d6de4a	r832: fixed outdated -L. Resolved #231 and #233	2018-09-14 22:21:33 -04:00
Chris Rands	791e89ce0f	PEP8 and other minor styles changes	2018-09-12 10:58:56 -04:00
Heng Li	98c48a1c45	Fixed wrong version links in the cookbook	2018-09-03 09:07:34 -04:00
Heng Li	63d397120a	removed getopt.c from MANIFEST	2018-09-01 21:30:46 -04:00
Heng Li	7998fe9906	r829: replaced musl's getopt with ketopt	2018-09-01 21:18:02 -04:00
Heng Li	3a119d606f	r828: --MD to support spliced alignment	2018-08-22 10:47:45 -04:00
Heng Li	a5eafb75f9	Release minimap2-2.12 (r827)	2018-08-06 12:44:39 -04:00
Heng Li	9a567e4b37	allow mappy to change scoring	2018-08-06 09:52:13 -04:00
Heng Li	8e606bcc06	added a `-C5` example to Getting Started	2018-08-06 09:08:33 -04:00
Heng Li	a1b7219b5d	explain added parameters	2018-08-05 21:32:05 -04:00
Heng Li	b0f39a1a61	r823: mappy to index a single sequence	2018-08-05 20:57:05 -04:00
Heng Li	5ab6538757	r822: added option --no-end-flt	2018-08-05 19:42:12 -04:00
Heng Li	b32296e18f	r821: fixed memory when -y is used	2018-07-31 15:14:37 -04:00
Heng Li	99ecdf7b5d	mappy to support arm64 (#203 )	2018-07-24 23:53:02 -04:00
Heng Li	ff9917a1c4	r819: mappy to support cs/MD	2018-07-24 23:29:55 -04:00
Mark Bicknell	8c064a5f29	Fixed mm_idx_is_idx to return the correct result on Windows.	2018-07-17 09:16:33 -04:00
Heng Li	0e137670fc	Merge branch 'split'	2018-07-15 22:24:15 -04:00
Heng Li	395c8d678a	r815: fixed a memory leak	2018-07-15 22:11:32 -04:00
Heng Li	830da7fa27	r814: resumed versioning	2018-07-15 11:48:14 -04:00
Heng Li	a655cbef86	print SAM header; remove tmp files	2018-07-15 11:03:18 -04:00
Heng Li	4b707aac92	working with toy examples	2018-07-15 10:55:00 -04:00
Heng Li	951c0d1d35	apparently mm_append_cigar() wastes some memory	2018-07-14 23:47:44 -04:00
Heng Li	3545e35a42	pairing in the split-idx mode	2018-07-14 23:43:34 -04:00
Heng Li	f3417da838	bugfix: unmapped records are duplicated in output	2018-07-14 22:54:05 -04:00
Heng Li	e5277dbf5c	code backup	2018-07-14 22:52:36 -04:00
Heng Li	1a55227d5a	write hits to tmp files (unfinished)	2018-07-14 12:15:10 -04:00
Heng Li	5cfa621b2d	use unmapped records	2018-07-07 12:49:15 -05:00
Heng Li	a609a07f8c	optionally output unmapped query in PAF	2018-07-07 10:26:08 -05:00
Heng Li	bcf92b3c46	compute query coverage	2018-07-06 21:52:01 -05:00
Heng Li	097378ab90	reworked break point counting	2018-07-06 09:46:26 -04:00
Heng Li	10bbbe28c5	added asmstat	2018-07-05 13:22:32 -04:00
Hyeshik Chang	c92a6866f3	Release the GIL to allow native Python threading.	2018-07-05 08:00:27 -04:00
Maël Kerbiriou	6908dc59a5	--splice implied when searching for splicing sites on single strand	2018-07-05 07:41:26 -04:00
Heng Li	50dae10421	paftools call to show statistics on longer indels	2018-07-03 14:28:09 -04:00
Heng Li	0517972d02	Release minimap2-2.11 (r797)	2018-06-21 00:04:08 -04:00
Heng Li	d46e68e6ad	r796: don't use ssize_t	2018-06-20 12:45:27 -04:00
Heng Li	2584a4149a	r295: use -r2000 for ava-ont, NOT for ava-pb	2018-06-20 12:24:43 -04:00
Heng Li	66674afd09	r794: fixed a bug in seed filtering	2018-06-20 10:26:29 -04:00
Heng Li	e9ca0c9dab	added __version__; resolved #165 Not sure if this is the right way. Apparently working.	2018-06-19 15:40:26 -04:00
Heng Li	7e6e8ca73f	r792: fixed -Wextra warnings and resolved #184	2018-06-19 15:26:58 -04:00
Ilya Kolpakov	408e098859	fix deserialization of zero-length reference names	2018-06-19 14:46:05 -04:00
Ilya Kolpakov	57f37551f8	expose mm_idx_is_idx, mm_idx_load and mm_idx_dump	2018-06-19 14:46:05 -04:00
Ilya Kolpakov	4c66b689c3	fix serialization of empty names in mm_idx_dump	2018-06-19 14:46:05 -04:00
mvdbeek	1bde2cf076	Allow setting max_frag_len on a per alignment level	2018-06-19 14:39:30 -04:00
mvdbeek	31fc0f218a	Allow setting max_frag_len parameter in Aligner class	2018-06-19 14:39:30 -04:00
Aaron Wenger	3d3bcc29a8	Fix CIGAR reallocation with --eqx Fix the logic that calculates the number of CIGAR entries when match "M" entries are expanded into "=" and "X". The number of entries depends not on the number of mismatches but rather on the number of transitions between "=" to "X".	2018-06-19 14:37:41 -04:00
Hasindu Gamaarachchi	99dcd75f64	added support for 64 bit ARM architectures	2018-06-19 14:28:45 -04:00
Heng Li	154d2caf5b	r784: support the =/X CIGAR operators (#156 )	2018-05-30 16:11:22 -04:00
Heng Li	a3afeec0b2	r783: reverted to r781 (#155 )	2018-05-30 15:25:34 -04:00
Heng Li	3573784b4d	r782: no mask a chain having long ref ovlp (#155 )	2018-05-30 13:53:45 -04:00
Heng Li	872f300955	r781: fixed the buggy heapmerge (resolves #166 )	2018-05-30 11:55:14 -04:00
Heng Li	d7b61a039e	updated citation	2018-05-30 11:11:55 -04:00
Heng Li	248158a3e1	Resolved #168 : citing the manpage for cs	2018-05-30 11:03:16 -04:00
Heng Li	9f4309c376	r777: avoid skipping too many seeds	2018-05-11 10:25:18 -04:00
Heng Li	463f9309f9	Merge branch 'hot-fix' into fix-long-gap	2018-05-11 10:13:19 -04:00
Heng Li	abe989e355	the previous fix on int overflow is incomplete	2018-05-11 10:12:57 -04:00
Heng Li	881b4ca3a2	r774: Merge branch 'hot-fix' into fix-long-gap	2018-05-11 10:02:17 -04:00
Heng Li	10c6dd2551	r773: fixed an integer overflow	2018-05-11 10:01:23 -04:00
Heng Li	7ec6721c44	r772: option -Y not working	2018-05-11 10:00:11 -04:00
Heng Li	e61812ee55	reduced gap len to trigger bad seed filtering	2018-05-01 16:17:21 -04:00
Heng Li	734ac379bb	r770: matching N bases not working properly (#155 )	2018-04-30 19:55:23 -04:00
Heng Li	759f8e4ac9	r769: filter out seeds breaking long gaps	2018-04-24 15:37:37 -04:00
Heng Li	aef7b0744c	r768: shortened preset; added dv tag (#25 ) Also added asm20 to command line help (#151)	2018-04-24 12:48:54 -04:00
Heng Li	39f836eac8	r767: don't crash when there is no "cg" (#153 )	2018-04-24 12:32:37 -04:00
Heng Li	cbeb86dad6	Merge remote-tracking branch 'origin/master'	2018-04-10 09:12:26 -04:00
Heng Li	372c90ceb5	r764: fixed incorrect inversion mapq (#148 )	2018-04-10 09:11:49 -04:00
apregier	2e2e69107c	Small bugfix for paftools.js The bug resulted in the wrong coverage being printed in some cases.	2018-04-04 16:29:54 -04:00
Heng Li	ee4cd089f7	r763: fine control long join flank len (#128 )	2018-03-29 14:16:58 -04:00
Heng Li	2d7ec75d50	Release minimap2-2.10 (r761)	2018-03-27 11:45:44 -04:00
Heng Li	7938ed4893	fixed two mappy warnings	2018-03-26 15:12:54 -04:00
Heng Li	4740423afa	Merge remote-tracking branch 'origin/master'	2018-03-26 14:37:17 -04:00
Heng Li	1776311a9b	updated cookbook	2018-03-26 14:37:02 -04:00
Heng Li	c1a3e05cb0	Merge pull request #138 from zingdle/patch-1 Fix typo in README.md	2018-03-26 07:50:39 -04:00
zingdle	ecb6703d4a	Update README.md Fix typo. TD;DR -> TL;DR	2018-03-26 14:54:52 +08:00
Heng Li	0bf97a367b	r755: junceval to only evaluate chromosomes	2018-03-24 23:10:06 -04:00
Heng Li	1c504d72e4	r754: option to only change name	2018-03-23 12:27:11 -04:00
Heng Li	5ef9580b17	r753: change bandwidth in ava-ont to 2000bp	2018-03-23 10:15:23 -04:00
Heng Li	08bd2123b6	r752: option to copy comments to output (#136 )	2018-03-23 10:04:33 -04:00
Heng Li	8766d286df	r751: optionally output MD (#118 )	2018-03-22 14:15:33 -04:00
Heng Li	623b5d9d48	r750: check puts() return (#132 & #103 )	2018-03-22 11:31:58 -04:00
Heng Li	18659118cd	r749: don't print version etc at low verbose	2018-03-22 11:10:55 -04:00
Heng Li	d1050f4eaf	r748: optionally to use system getopt() (#134 )	2018-03-19 11:18:26 -04:00
Heng Li	b81d45510e	Revision v2	2018-03-16 11:18:06 -04:00
Heng Li	d135feb1a5	response to reviewers' comments, round 2	2018-03-15 21:59:57 -04:00
Heng Li	242ff4e91d	r745: junceval - recognize "-" as file name	2018-03-15 12:51:46 -04:00
Heng Li	7a0c1316ce	added more examples to cookbook	2018-03-12 22:24:22 -04:00
Heng Li	77ebd479f4	r743: added VCF output to "call"	2018-03-12 21:51:31 -04:00
Heng Li	e3f226a9d9	working on section "Full-Genome Alignment" found an apparent bug in paftools.js call. To debug...	2018-03-12 16:12:18 -04:00
Heng Li	bdc615c1d4	r741: added --min-occ-floor to improve #107	2018-03-12 14:32:27 -04:00
Heng Li	ad1beaf255	backup; not finished	2018-03-12 13:20:25 -04:00
Heng Li	de0480ac5b	finished section "Read Overlap"	2018-03-12 13:05:51 -04:00
Heng Li	f2866533a8	finished "mapping genomic reads"	2018-03-12 12:50:21 -04:00
Heng Li	0173850ef0	don't use bullets - they are hard to read	2018-03-12 12:42:13 -04:00
Heng Li	acea3594fb	updated cookbook	2018-03-12 12:39:57 -04:00
Heng Li	f78a247749	Started to work on a cookbook	2018-03-12 12:00:44 -04:00
Heng Li	ccaf12e1a2	r734: removed debugging print in call()	2018-03-09 23:46:21 -05:00
Heng Li	96b132c97d	fixed a bug/typo in Aligner.seq() (#126 )	2018-03-08 19:12:12 -05:00
Heng Li	70428ca3a8	speedup tag parsing for sam2paf	2018-03-02 10:17:55 -05:00
Heng Li	9aea79d621	faster tag parsing	2018-03-02 10:09:14 -05:00
Heng Li	1770988627	make calling work with sam2paf output	2018-03-02 10:03:06 -05:00
Heng Li	2bfdad34bb	minor cleanup to last commit	2018-03-02 01:07:37 -05:00
Heng Li	953766cedd	generate cs; not carefully tested	2018-03-02 00:14:55 -05:00
Heng Li	dc61301d9f	sam2paf cleanup	2018-03-01 21:58:12 -05:00
Heng Li	19e05a099d	for clarity	2018-03-01 18:25:29 -05:00
Heng Li	0238caa8b1	clarify minimap2 not working for >2Gbp seq #129	2018-03-01 18:23:40 -05:00
Heng Li	a22ebb9836	use SSE compiler flags more precisely (#127 )	2018-02-26 09:51:01 -05:00