Release minimap2-2.13 (r850)

I don't know why I haven't found it at the beginning.

MANIFEST.in

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

Makefile

@@ -1,7 +1,7 @@
 CFLAGS=		-g -Wall -O2 -Wc++-compat #-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
@@ -32,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)
@@ -41,8 +41,8 @@ minimap2-lite:example.o libminimap2.a
 libminimap2.a:$(OBJS)
 		$(AR) -csru $@ $(OBJS)
 
-sdust:sdust.c getopt.o kalloc.o kalloc.h kdq.h kvec.h kseq.h sdust.h
-		$(CC) -D_SDUST_MAIN $(CFLAGS) $< getopt.o kalloc.o -o $@ -lz
+sdust:sdust.c kalloc.o kalloc.h kdq.h kvec.h kseq.h ketopt.h sdust.h
+		$(CC) -D_SDUST_MAIN $(CFLAGS) $< kalloc.o -o $@ -lz
 
 # SSE-specific targets on x86/x86_64
 
@@ -99,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
@@ -108,11 +107,12 @@ ksw2_exts2_sse.o: ksw2.h kalloc.h
 ksw2_extz2_sse.o: ksw2.h kalloc.h
 ksw2_ll_sse.o: ksw2.h kalloc.h
 kthread.o: kthread.h
-main.o: bseq.h minimap.h mmpriv.h getopt.h
+main.o: bseq.h minimap.h mmpriv.h ketopt.h
 map.o: kthread.h kvec.h kalloc.h sdust.h mmpriv.h minimap.h bseq.h khash.h
 map.o: ksort.h
 misc.o: mmpriv.h minimap.h bseq.h ksort.h
 options.o: mmpriv.h minimap.h bseq.h
 pe.o: mmpriv.h minimap.h bseq.h kvec.h kalloc.h ksort.h
-sdust.o: kalloc.h kdq.h kvec.h sdust.h
+sdust.o: kalloc.h kdq.h kvec.h ketopt.h sdust.h
 sketch.o: kvec.h kalloc.h mmpriv.h minimap.h bseq.h
+splitidx.o: mmpriv.h minimap.h bseq.h

NEWS.md

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

README.md

@@ -14,9 +14,11 @@ cd minimap2 && make
 # use presets (no test data)
 ./minimap2 -ax map-pb ref.fa pacbio.fq.gz > aln.sam       # PacBio genomic reads
 ./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam         # Oxford Nanopore genomic reads
+./minimap2 -ax asm20 ref.fa pacbio-ccs.fq.gz > aln.sam    # PacBio CCS genomic reads
 ./minimap2 -ax sr ref.fa read1.fa read2.fa > aln.sam      # short genomic paired-end reads
-./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads
-./minimap2 -ax splice -k14 -uf ref.fa reads.fa > aln.sam  # Nanopore Direct RNA-seq
+./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads (strand unknown)
+./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam  # noisy Nanopore Direct RNA-seq
+./minimap2 -ax splice -uf -C5 ref.fa query.fa > aln.sam   # Final PacBio Iso-seq or traditional cDNA
 ./minimap2 -cx asm5 asm1.fa asm2.fa > aln.paf             # intra-species asm-to-asm alignment
 ./minimap2 -x ava-pb reads.fa reads.fa > overlaps.paf     # PacBio read overlap
 ./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf    # Nanopore read overlap
@@ -69,8 +71,8 @@ Detailed evaluations are available from the [minimap2 paper][doi] or the
 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.11/minimap2-2.11_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.11_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

align.c

@@ -199,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); // TODO: should this be "n_cigar + sizeof(mm_extra_t)/4" instead?
+		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);
 	}
@@ -563,11 +563,12 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 		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);
@@ -878,6 +879,6 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 	kfree(km, qseq0[0]);
 	kfree(km, ez.cigar);
 	mm_filter_regs(opt, qlen, n_regs_, regs);
-	mm_hit_sort_by_dp(km, n_regs_, regs);
+	mm_hit_sort(km, n_regs_, regs);
 	return regs;
 }

cookbook.md

@@ -24,18 +24,18 @@
 
 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.11][v2.11]. Some examples only work
-with v2.11 or later.
+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.11/minimap2-2.11_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.11_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
+curl -L https://github.com/lh3/minimap2/releases/download/v2.13/minimap2-2.13_x64-linux.tar.bz2 | tar jxf -
+cp minimap2-2.13_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
 export PATH="$PATH:"`pwd`                               # put the current directory on PATH
 # download example datasets
-curl -L https://github.com/lh3/minimap2/releases/download/v2.11/cookbook-data.tgz | tar zxf -
+curl -L https://github.com/lh3/minimap2/releases/download/v2.10/cookbook-data.tgz | tar zxf -
 ```
 
 ## <a name="map-reads"></a>Mapping Genomic Reads
@@ -240,4 +240,4 @@ 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.11]: https://github.com/lh3/minimap2/releases/tag/v2.11
+[v2.10]: https://github.com/lh3/minimap2/releases/tag/v2.10

format.c

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

getopt.c

@@ -1,216 +0,0 @@
-#include <stddef.h>
-#include <stdio.h>
-#include <string.h>
-#include "getopt.h"
-
-char *optarg;
-int optind=1, opterr=1, optopt, __optpos, optreset=0;
-
-#define optpos __optpos
-
-static void __getopt_msg(const char *a, const char *b, const char *c, size_t l)
-{
-	FILE *f = stderr;
-#if !defined(WIN32) && !defined(_WIN32)
-	flockfile(f);
-#endif
-	fputs(a, f);
-	fwrite(b, strlen(b), 1, f);
-	fwrite(c, 1, l, f);
-	fputc('\n', f);
-#if !defined(WIN32) && !defined(_WIN32)
-	funlockfile(f);
-#endif
-}
-
-int getopt(int argc, char * const argv[], const char *optstring)
-{
-	int i, c, d;
-	int k, l;
-	char *optchar;
-
-	if (!optind || optreset) {
-		optreset = 0;
-		__optpos = 0;
-		optind = 1;
-	}
-
-	if (optind >= argc || !argv[optind])
-		return -1;
-
-	if (argv[optind][0] != '-') {
-		if (optstring[0] == '-') {
-			optarg = argv[optind++];
-			return 1;
-		}
-		return -1;
-	}
-
-	if (!argv[optind][1])
-		return -1;
-
-	if (argv[optind][1] == '-' && !argv[optind][2])
-		return optind++, -1;
-
-	if (!optpos) optpos++;
-	c = argv[optind][optpos], k = 1;
-	optchar = argv[optind]+optpos;
-	optopt = c;
-	optpos += k;
-
-	if (!argv[optind][optpos]) {
-		optind++;
-		optpos = 0;
-	}
-
-	if (optstring[0] == '-' || optstring[0] == '+')
-		optstring++;
-
-	i = 0;
-	d = 0;
-	do {
-		d = optstring[i], l = 1;
-		if (l>0) i+=l; else i++;
-	} while (l && d != c);
-
-	if (d != c) {
-		if (optstring[0] != ':' && opterr)
-			__getopt_msg(argv[0], ": unrecognized option: ", optchar, k);
-		return '?';
-	}
-	if (optstring[i] == ':') {
-		if (optstring[i+1] == ':') optarg = 0;
-		else if (optind >= argc) {
-			if (optstring[0] == ':') return ':';
-			if (opterr) __getopt_msg(argv[0],
-				": option requires an argument: ",
-				optchar, k);
-			return '?';
-		}
-		if (optstring[i+1] != ':' || optpos) {
-			optarg = argv[optind++] + optpos;
-			optpos = 0;
-		}
-	}
-	return c;
-}
-
-static void permute(char *const *argv, int dest, int src)
-{
-	char **av = (char **)argv;
-	char *tmp = av[src];
-	int i;
-	for (i=src; i>dest; i--)
-		av[i] = av[i-1];
-	av[dest] = tmp;
-}
-
-static int __getopt_long_core(int argc, char *const *argv, const char *optstring, const struct option *longopts, int *idx, int longonly)
-{
-	optarg = 0;
-	if (longopts && argv[optind][0] == '-' &&
-		((longonly && argv[optind][1] && argv[optind][1] != '-') ||
-		 (argv[optind][1] == '-' && argv[optind][2])))
-	{
-		int colon = optstring[optstring[0]=='+'||optstring[0]=='-']==':';
-		int i, cnt, match = -1;
-		char *opt;
-		for (cnt=i=0; longopts[i].name; i++) {
-			const char *name = longopts[i].name;
-			opt = argv[optind]+1;
-			if (*opt == '-') opt++;
-			for (; *name && *name == *opt; name++, opt++);
-			if (*opt && *opt != '=') continue;
-			match = i;
-			if (!*name) {
-				cnt = 1;
-				break;
-			}
-			cnt++;
-		}
-		if (cnt==1) {
-			i = match;
-			optind++;
-			optopt = longopts[i].val;
-			if (*opt == '=') {
-				if (!longopts[i].has_arg) {
-					if (colon || !opterr)
-						return '?';
-					__getopt_msg(argv[0],
-						": option does not take an argument: ",
-						longopts[i].name,
-						strlen(longopts[i].name));
-					return '?';
-				}
-				optarg = opt+1;
-			} else if (longopts[i].has_arg == required_argument) {
-				if (!(optarg = argv[optind])) {
-					if (colon) return ':';
-					if (!opterr) return '?';
-					__getopt_msg(argv[0],
-						": option requires an argument: ",
-						longopts[i].name,
-						strlen(longopts[i].name));
-					return '?';
-				}
-				optind++;
-			}
-			if (idx) *idx = i;
-			if (longopts[i].flag) {
-				*longopts[i].flag = longopts[i].val;
-				return 0;
-			}
-			return longopts[i].val;
-		}
-		if (argv[optind][1] == '-') {
-			if (!colon && opterr)
-				__getopt_msg(argv[0], cnt ?
-					": option is ambiguous: " :
-					": unrecognized option: ",
-					argv[optind]+2,
-					strlen(argv[optind]+2));
-			optind++;
-			return '?';
-		}
-	}
-	return getopt(argc, argv, optstring);
-}
-
-static int __getopt_long(int argc, char *const *argv, const char *optstring, const struct option *longopts, int *idx, int longonly)
-{
-	int ret, skipped, resumed;
-	if (!optind || optreset) {
-		optreset = 0;
-		__optpos = 0;
-		optind = 1;
-	}
-	if (optind >= argc || !argv[optind]) return -1;
-	skipped = optind;
-	if (optstring[0] != '+' && optstring[0] != '-') {
-		int i;
-		for (i=optind; ; i++) {
-			if (i >= argc || !argv[i]) return -1;
-			if (argv[i][0] == '-' && argv[i][1]) break;
-		}
-		optind = i;
-	}
-	resumed = optind;
-	ret = __getopt_long_core(argc, argv, optstring, longopts, idx, longonly);
-	if (resumed > skipped) {
-		int i, cnt = optind-resumed;
-		for (i=0; i<cnt; i++)
-			permute(argv, skipped, optind-1);
-		optind = skipped + cnt;
-	}
-	return ret;
-}
-
-int getopt_long(int argc, char *const *argv, const char *optstring, const struct option *longopts, int *idx)
-{
-	return __getopt_long(argc, argv, optstring, longopts, idx, 0);
-}
-
-int getopt_long_only(int argc, char *const *argv, const char *optstring, const struct option *longopts, int *idx)
-{
-	return __getopt_long(argc, argv, optstring, longopts, idx, 1);
-}

getopt.h

@@ -1,53 +0,0 @@
-/*
-  Copyright 2005-2014 Rich Felker, et al.
-  
-  Permission is hereby granted, free of charge, to any person obtaining
-  a copy of this software and associated documentation files (the
-  "Software"), to deal in the Software without restriction, including
-  without limitation the rights to use, copy, modify, merge, publish,
-  distribute, sublicense, and/or sell copies of the Software, and to
-  permit persons to whom the Software is furnished to do so, subject to
-  the following conditions:
-  
-  The above copyright notice and this permission notice shall be
-  included in all copies or substantial portions of the Software.
-  
-  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
-  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
-  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
-  IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
-  CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
-  TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
-  SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-*/
-
-#ifndef _GETOPT_H
-#define _GETOPT_H
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-int getopt(int, char * const [], const char *);
-extern char *optarg;
-extern int optind, opterr, optopt, optreset;
-
-struct option {
-	const char *name;
-	int has_arg;
-	int *flag;
-	int val;
-};
-
-int getopt_long(int, char *const *, const char *, const struct option *, int *);
-int getopt_long_only(int, char *const *, const char *, const struct option *, int *);
-
-#define no_argument        0
-#define required_argument  1
-#define optional_argument  2
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif

hit.c

@@ -106,7 +106,7 @@ void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a)
 	r->split |= 1, r2->split |= 2;
 }
 
-void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff) // and compute mm_reg1_t::subsc
+void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff, int hard_mask_level) // and compute mm_reg1_t::subsc
 {
 	int i, j, k, *w;
 	uint64_t *cov;
@@ -118,6 +118,7 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff
 	for (i = 1, k = 1; i < n; ++i) {
 		mm_reg1_t *ri = &r[i];
 		int si = ri->qs, ei = ri->qe, n_cov = 0, uncov_len = 0;
+		if (hard_mask_level) goto skip_uncov;
 		for (j = 0; j < k; ++j) { // traverse existing primary hits to find overlapping hits
 			mm_reg1_t *rp = &r[w[j]];
 			int sj = rp->qs, ej = rp->qe;
@@ -137,6 +138,7 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff
 			}
 			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;
@@ -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];

index.c

@@ -48,10 +48,12 @@ void mm_idx_destroy(mm_idx_t *mi)
 	uint32_t i;
 	if (mi == 0) return;
 	if (mi->h) kh_destroy(str, (khash_t(str)*)mi->h);
-	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->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)
@@ -370,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]);
@@ -381,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);
@@ -510,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;
@@ -563,7 +576,7 @@ mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads)
 		mi = mm_idx_gen(r->fp.seq, r->opt.w, r->opt.k, r->opt.bucket_bits, r->opt.flag, r->opt.mini_batch_size, n_threads, r->opt.batch_size);
 	if (mi) {
 		if (r->fp_out) mm_idx_dump(r->fp_out, mi);
-		++r->n_parts;
+		mi->index = r->n_parts++;
 	}
 	return mi;
 }

ketopt.h

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

ksw2_exts2_sse.c

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

main.c

@@ -4,13 +4,9 @@
 #include "bseq.h"
 #include "minimap.h"
 #include "mmpriv.h"
-#ifdef HAVE_GETOPT
-#include <getopt.h>
-#else
-#include "getopt.h"
-#endif
+#include "ketopt.h"
 
-#define MM_VERSION "2.11-r797"
+#define MM_VERSION "2.13-r850"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -26,49 +22,53 @@ void liftrlimit()
 void liftrlimit() {}
 #endif
 
-static struct option long_options[] = {
-	{ "bucket-bits",    required_argument, 0, 0 },
-	{ "mb-size",        required_argument, 0, 'K' },
-	{ "seed",           required_argument, 0, 0 },
-	{ "no-kalloc",      no_argument,       0, 0 },
-	{ "print-qname",    no_argument,       0, 0 },
-	{ "no-self",        no_argument,       0, 'D' },
-	{ "print-seeds",    no_argument,       0, 0 },
-	{ "max-chain-skip", required_argument, 0, 0 },
-	{ "min-dp-len",     required_argument, 0, 0 },
-	{ "print-aln-seq",  no_argument,       0, 0 },
-	{ "splice",         no_argument,       0, 0 },
-	{ "cost-non-gt-ag", required_argument, 0, 'C' },
-	{ "no-long-join",   no_argument,       0, 0 },
-	{ "sr",             no_argument,       0, 0 },
-	{ "frag",           required_argument, 0, 0 },
-	{ "secondary",      required_argument, 0, 0 },
-	{ "cs",             optional_argument, 0, 0 },
-	{ "end-bonus",      required_argument, 0, 0 },
-	{ "no-pairing",     no_argument,       0, 0 },
-	{ "splice-flank",   required_argument, 0, 0 },
-	{ "idx-no-seq",     no_argument,       0, 0 },
-	{ "end-seed-pen",   required_argument, 0, 0 },   // 21
-	{ "for-only",       no_argument,       0, 0 },   // 22
-	{ "rev-only",       no_argument,       0, 0 },   // 23
-	{ "heap-sort",      required_argument, 0, 0 },   // 24
-	{ "all-chain",      no_argument,       0, 'P' },
-	{ "dual",           required_argument, 0, 0 },   // 26
-	{ "max-clip-ratio", required_argument, 0, 0 },   // 27
-	{ "min-occ-floor",  required_argument, 0, 0 },   // 28
-	{ "MD",             no_argument,       0, 0 },   // 29
-	{ "lj-min-ratio",   required_argument, 0, 0 },   // 30
-	{ "score-N",        required_argument, 0, 0 },   // 31
-	{ "eqx",            no_argument,       0, 0 },   // 32
-	{ "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)
@@ -97,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:yY";
+	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;
@@ -111,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;
 		}
-	optind = 0; // for musl getopt, optind=0 has the same effect as optreset=1; older libc doesn't have optreset
+	}
+	o = KETOPT_INIT;
 
-	while ((c = getopt_long(argc, argv, opt_str, long_options, &long_idx)) >= 0) {
-		if (c == 'w') ipt.w = atoi(optarg);
-		else if (c == 'k') ipt.k = atoi(optarg);
+	while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) {
+		if (c == 'w') ipt.w = atoi(o.arg);
+		else if (c == 'k') ipt.k = atoi(o.arg);
 		else if (c == 'H') ipt.flag |= MM_I_HPC;
-		else if (c == 'd') fnw = optarg; // the above are indexing related options, except -I
-		else if (c == 'r') opt.bw = (int)mm_parse_num(optarg);
-		else if (c == 't') n_threads = atoi(optarg);
-		else if (c == 'v') mm_verbose = atoi(optarg);
-		else if (c == 'g') opt.max_gap = (int)mm_parse_num(optarg);
-		else if (c == 'G') mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(optarg));
-		else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(optarg);
-		else if (c == 'N') opt.best_n = atoi(optarg);
-		else if (c == 'p') opt.pri_ratio = atof(optarg);
-		else if (c == 'M') opt.mask_level = atof(optarg);
+		else if (c == 'd') fnw = o.arg; // the above are indexing related options, except -I
+		else if (c == 'r') opt.bw = (int)mm_parse_num(o.arg);
+		else if (c == 't') n_threads = atoi(o.arg);
+		else if (c == 'v') mm_verbose = atoi(o.arg);
+		else if (c == 'g') opt.max_gap = (int)mm_parse_num(o.arg);
+		else if (c == 'G') mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(o.arg));
+		else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(o.arg);
+		else if (c == 'N') old_best_n = opt.best_n, opt.best_n = atoi(o.arg);
+		else if (c == 'p') opt.pri_ratio = atof(o.arg);
+		else if (c == 'M') opt.mask_level = atof(o.arg);
 		else if (c == 'c') opt.flag |= MM_F_OUT_CG | MM_F_CIGAR;
 		else if (c == 'D') opt.flag |= MM_F_NO_DIAG;
 		else if (c == 'P') opt.flag |= MM_F_ALL_CHAINS;
@@ -144,62 +151,66 @@ int main(int argc, char *argv[])
 		else if (c == 'Y') opt.flag |= MM_F_SOFTCLIP;
 		else if (c == 'L') opt.flag |= MM_F_LONG_CIGAR;
 		else if (c == 'y') opt.flag |= MM_F_COPY_COMMENT;
-		else if (c == 'T') opt.sdust_thres = atoi(optarg);
-		else if (c == 'n') opt.min_cnt = atoi(optarg);
-		else if (c == 'm') opt.min_chain_score = atoi(optarg);
-		else if (c == 'A') opt.a = atoi(optarg);
-		else if (c == 'B') opt.b = atoi(optarg);
-		else if (c == 's') opt.min_dp_max = atoi(optarg);
-		else if (c == 'C') opt.noncan = atoi(optarg);
-		else if (c == 'I') ipt.batch_size = mm_parse_num(optarg);
-		else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(optarg);
-		else if (c == 'R') rg = optarg;
+		else if (c == 'T') opt.sdust_thres = atoi(o.arg);
+		else if (c == 'n') opt.min_cnt = atoi(o.arg);
+		else if (c == 'm') opt.min_chain_score = atoi(o.arg);
+		else if (c == 'A') opt.a = atoi(o.arg);
+		else if (c == 'B') opt.b = atoi(o.arg);
+		else if (c == 's') opt.min_dp_max = atoi(o.arg);
+		else if (c == 'C') opt.noncan = atoi(o.arg);
+		else if (c == 'I') ipt.batch_size = mm_parse_num(o.arg);
+		else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(o.arg);
+		else if (c == 'R') rg = o.arg;
 		else if (c == 'h') fp_help = stdout;
 		else if (c == '2') opt.flag |= MM_F_2_IO_THREADS;
-		else if (c == 0 && long_idx == 0) ipt.bucket_bits = atoi(optarg); // --bucket-bits
-		else if (c == 0 && long_idx == 2) opt.seed = atoi(optarg); // --seed
-		else if (c == 0 && long_idx == 3) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
-		else if (c == 0 && long_idx == 4) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
-		else if (c == 0 && long_idx == 6) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
-		else if (c == 0 && long_idx == 7) opt.max_chain_skip = atoi(optarg); // --max-chain-skip
-		else if (c == 0 && long_idx == 8) opt.min_ksw_len = atoi(optarg); // --min-dp-len
-		else if (c == 0 && long_idx == 9) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
-		else if (c == 0 && long_idx ==10) opt.flag |= MM_F_SPLICE; // --splice
-		else if (c == 0 && long_idx ==12) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
-		else if (c == 0 && long_idx ==13) opt.flag |= MM_F_SR; // --sr
-		else if (c == 0 && long_idx ==17) opt.end_bonus = atoi(optarg); // --end-bonus
-		else if (c == 0 && long_idx ==18) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
-		else if (c == 0 && long_idx ==20) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
-		else if (c == 0 && long_idx ==21) opt.anchor_ext_shift = atoi(optarg); // --end-seed-pen
-		else if (c == 0 && long_idx ==22) opt.flag |= MM_F_FOR_ONLY; // --for-only
-		else if (c == 0 && long_idx ==23) opt.flag |= MM_F_REV_ONLY; // --rev-only
-		else if (c == 0 && long_idx ==27) opt.max_clip_ratio = atof(optarg); // --max-clip-ratio
-		else if (c == 0 && long_idx ==28) opt.min_mid_occ = atoi(optarg); // --min-occ-floor
-		else if (c == 0 && long_idx ==29) opt.flag |= MM_F_OUT_MD; // --MD
-		else if (c == 0 && long_idx ==30) opt.min_join_flank_ratio = atof(optarg); // --lj-min-ratio
-		else if (c == 0 && long_idx ==31) opt.sc_ambi = atoi(optarg); // --score-N
-		else if (c == 0 && long_idx ==32) opt.flag |= MM_F_EQX; // --eqx
-		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)
@@ -210,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);
 		}
 	}
@@ -242,8 +253,12 @@ 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");
@@ -293,20 +308,20 @@ int main(int argc, char *argv[])
 		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 command-line options.\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;
@@ -325,25 +340,29 @@ int main(int argc, char *argv[])
 				mm_write_sam_hdr(mi, rg, MM_VERSION, argc, argv);
 			} else {
 				mm_write_sam_hdr(0, rg, MM_VERSION, argc, argv);
-				if (mm_verbose >= 2)
-					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);
 
+	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);
@@ -354,7 +373,7 @@ int main(int argc, char *argv[])
 		fprintf(stderr, "[M::%s] CMD:", __func__);
 		for (i = 0; i < argc; ++i)
 			fprintf(stderr, " %s", argv[i]);
-		fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec\n", __func__, realtime() - mm_realtime0, cputime());
+		fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, realtime() - mm_realtime0, cputime(), peakrss() / 1024.0 / 1024.0 / 1024.0);
 	}
 	return 0;
 }

map.c

@@ -11,6 +11,7 @@
 
 struct mm_tbuf_s {
 	void *km;
+	int rep_len, frag_gap;
 };
 
 mm_tbuf_t *mm_tbuf_init(void)
@@ -28,6 +29,11 @@ void mm_tbuf_destroy(mm_tbuf_t *b)
 	free(b);
 }
 
+void *mm_tbuf_get_km(mm_tbuf_t *b)
+{
+	return b->km;
+}
+
 static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *seq, int sdust_thres)
 {
 	int n_dreg, j, k, u = 0;
@@ -242,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
@@ -255,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);
 	}
@@ -330,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);
 
@@ -351,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);
 		}
@@ -394,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;
@@ -421,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)))) {
@@ -440,6 +459,63 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
 		}
 }
 
+static void merge_hits(step_t *s)
+{
+	int f, i, k0, k, max_seg = 0, *n_reg_part, *rep_len_part, *frag_gap_part, *qlens;
+	void *km;
+	FILE **fp = s->p->fp_parts;
+	const mm_mapopt_t *opt = s->p->opt;
+
+	km = km_init();
+	for (f = 0; f < s->n_frag; ++f)
+		max_seg = max_seg > s->n_seg[f]? max_seg : s->n_seg[f];
+	qlens = CALLOC(int, max_seg + s->p->n_parts * 3);
+	n_reg_part = qlens + max_seg;
+	rep_len_part = n_reg_part + s->p->n_parts;
+	frag_gap_part = rep_len_part + s->p->n_parts;
+	for (f = 0, k = k0 = 0; f < s->n_frag; ++f) {
+		k0 = k;
+		for (i = 0; i < s->n_seg[f]; ++i, ++k) {
+			int j, l, t, rep_len = 0;
+			qlens[i] = s->seq[k].l_seq;
+			for (j = 0, s->n_reg[k] = 0; j < s->p->n_parts; ++j) {
+				mm_err_fread(&n_reg_part[j],    sizeof(int), 1, fp[j]);
+				mm_err_fread(&rep_len_part[j],  sizeof(int), 1, fp[j]);
+				mm_err_fread(&frag_gap_part[j], sizeof(int), 1, fp[j]);
+				s->n_reg[k] += n_reg_part[j];
+				if (rep_len < rep_len_part[j])
+					rep_len = rep_len_part[j];
+			}
+			s->reg[k] = CALLOC(mm_reg1_t, s->n_reg[k]);
+			for (j = 0, l = 0; j < s->p->n_parts; ++j) {
+				for (t = 0; t < n_reg_part[j]; ++t, ++l) {
+					mm_reg1_t *r = &s->reg[k][l];
+					uint32_t capacity;
+					mm_err_fread(r, sizeof(mm_reg1_t), 1, fp[j]);
+					r->rid += s->p->rid_shift[j];
+					if (opt->flag & MM_F_CIGAR) {
+						mm_err_fread(&capacity, 4, 1, fp[j]);
+						r->p = (mm_extra_t*)calloc(capacity, 4);
+						r->p->capacity = capacity;
+						mm_err_fread(r->p, r->p->capacity, 4, fp[j]);
+					}
+				}
+			}
+			mm_hit_sort(km, &s->n_reg[k], s->reg[k]);
+			mm_set_parent(km, opt->mask_level, s->n_reg[k], s->reg[k], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL);
+			if (!(opt->flag & MM_F_ALL_CHAINS)) {
+				mm_select_sub(km, opt->pri_ratio, s->p->mi->k*2, opt->best_n, &s->n_reg[k], s->reg[k]);
+				mm_set_sam_pri(s->n_reg[k], s->reg[k]);
+			}
+			mm_set_mapq(km, s->n_reg[k], s->reg[k], opt->min_chain_score, opt->a, rep_len, !!(opt->flag & MM_F_SR));
+		}
+		if (s->n_seg[f] == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR))
+			mm_pair(km, frag_gap_part[0], opt->pe_bonus, opt->a * 2 + opt->b, opt->a, qlens, &s->n_reg[k0], &s->reg[k0]);
+	}
+	free(qlens);
+	km_destroy(km);
+}
+
 static void *worker_pipeline(void *shared, int step, void *in)
 {
 	int i, j, k;
@@ -459,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)) {
@@ -472,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;
@@ -485,19 +564,35 @@ static void *worker_pipeline(void *shared, int step, void *in)
 			int seg_st = s->seg_off[k], seg_en = s->seg_off[k] + s->n_seg[k];
 			for (i = seg_st; i < seg_en; ++i) {
 				mm_bseq1_t *t = &s->seq[i];
-				for (j = 0; j < s->n_reg[i]; ++j) {
-					mm_reg1_t *r = &s->reg[i][j];
-					assert(!r->sam_pri || r->id == r->parent);
-					if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent)
-						continue;
+				if (p->opt->split_prefix && p->n_parts == 0) { // then write to temporary files
+					mm_err_fwrite(&s->n_reg[i],    sizeof(int), 1, p->fp_split);
+					mm_err_fwrite(&s->rep_len[i],  sizeof(int), 1, p->fp_split);
+					mm_err_fwrite(&s->frag_gap[i], sizeof(int), 1, p->fp_split);
+					for (j = 0; j < s->n_reg[i]; ++j) {
+						mm_reg1_t *r = &s->reg[i][j];
+						mm_err_fwrite(r, sizeof(mm_reg1_t), 1, p->fp_split);
+						if (p->opt->flag & MM_F_CIGAR) {
+							mm_err_fwrite(&r->p->capacity, 4, 1, p->fp_split);
+							mm_err_fwrite(r->p, r->p->capacity, 4, p->fp_split);
+						}
+					}
+				} else if (s->n_reg[i] > 0) { // the query has at least one hit
+					for (j = 0; j < s->n_reg[i]; ++j) {
+						mm_reg1_t *r = &s->reg[i][j];
+						assert(!r->sam_pri || r->id == r->parent);
+						if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent)
+							continue;
+						if (p->opt->flag & MM_F_OUT_SAM)
+							mm_write_sam2(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
+						else
+							mm_write_paf(&p->str, mi, t, r, km, p->opt->flag);
+						mm_err_puts(p->str.s);
+					}
+				} else if (p->opt->flag & (MM_F_OUT_SAM|MM_F_PAF_NO_HIT)) { // output an empty hit, if requested
 					if (p->opt->flag & MM_F_OUT_SAM)
-						mm_write_sam2(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
+						mm_write_sam2(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
 					else
-						mm_write_paf(&p->str, mi, t, r, km, p->opt->flag);
-					mm_err_puts(p->str.s);
-				}
-				if (s->n_reg[i] == 0 && (p->opt->flag & MM_F_OUT_SAM)) { // write an unmapped record
-					mm_write_sam2(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
+						mm_write_paf(&p->str, mi, t, 0, 0, p->opt->flag);
 					mm_err_puts(p->str.s);
 				}
 			}
@@ -506,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);
@@ -517,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;
@@ -552,3 +660,48 @@ int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int
 {
 	return mm_map_file_frag(idx, 1, &fn, opt, n_threads);
 }
+
+int mm_split_merge(int n_segs, const char **fn, const mm_mapopt_t *opt, int n_split_idx)
+{
+	int i;
+	pipeline_t pl;
+	mm_idx_t *mi;
+	if (n_segs < 1 || n_split_idx < 1) return -1;
+	memset(&pl, 0, sizeof(pipeline_t));
+	pl.n_fp = n_segs;
+	pl.fp = open_bseqs(pl.n_fp, fn);
+	if (pl.fp == 0) return -1;
+	pl.opt = opt;
+	pl.mini_batch_size = opt->mini_batch_size;
+
+	pl.n_parts = n_split_idx;
+	pl.fp_parts  = CALLOC(FILE*, pl.n_parts);
+	pl.rid_shift = CALLOC(uint32_t, pl.n_parts);
+	pl.mi = mi = mm_split_merge_prep(opt->split_prefix, n_split_idx, pl.fp_parts, pl.rid_shift);
+	if (pl.mi == 0) {
+		free(pl.fp_parts);
+		free(pl.rid_shift);
+		return -1;
+	}
+	for (i = n_split_idx - 1; i > 0; --i)
+		pl.rid_shift[i] = pl.rid_shift[i - 1];
+	for (pl.rid_shift[0] = 0, i = 1; i < n_split_idx; ++i)
+		pl.rid_shift[i] += pl.rid_shift[i - 1];
+	if (opt->flag & MM_F_OUT_SAM)
+		for (i = 0; i < (int32_t)pl.mi->n_seq; ++i)
+			printf("@SQ\tSN:%s\tLN:%d\n", pl.mi->seq[i].name, pl.mi->seq[i].len);
+
+	kt_pipeline(2, worker_pipeline, &pl, 3);
+
+	free(pl.str.s);
+	mm_idx_destroy(mi);
+	free(pl.rid_shift);
+	for (i = 0; i < n_split_idx; ++i)
+		fclose(pl.fp_parts[i]);
+	free(pl.fp_parts);
+	for (i = 0; i < pl.n_fp; ++i)
+		mm_bseq_close(pl.fp[i]);
+	free(pl.fp);
+	mm_split_rm_tmp(opt->split_prefix, n_split_idx);
+	return 0;
+}

minimap.h

@@ -32,6 +32,9 @@
 #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
@@ -59,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)
@@ -135,6 +139,8 @@ typedef struct {
 	int32_t mid_occ;     // ignore seeds with occurrences above this threshold
 	int32_t max_occ;
 	int mini_batch_size; // size of a batch of query bases to process in parallel
+
+	const char *split_prefix;
 } mm_mapopt_t;
 
 // index reader
@@ -297,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
  *
@@ -333,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);

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "20 June 2018" "minimap2-2.11 (r797)" "Bioinformatics tools"
+.TH minimap2 1 "11 October 2018" "minimap2-2.13 (r850)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -193,7 +193,7 @@ Primarily used for all-vs-all read overlapping.
 .BI -p \ FLOAT
 Minimal secondary-to-primary score ratio to output secondary mappings [0.8].
 Between two chains overlaping over half of the shorter chain (controlled by
-.BR --mask-level ),
+.BR -M ),
 the chain with a lower score is secondary to the chain with a higher score.
 If the ratio of the scores is below
 .IR FLOAT ,
@@ -226,6 +226,11 @@ Mark as secondary a chain that overlaps with a better chain by
 .I FLOAT
 or more of the shorter chain [0.5]
 .TP
+.B --hard-mask-level
+Honor option
+.B -M
+and disable a heurstic to save unmapped subsequences.
+.TP
 .BI --max-chain-skip \ INT
 A heuristics that stops chaining early [50]. Minimap2 uses dynamic programming
 for chaining. The time complexity is quadratic in the number of seeds. This
@@ -239,7 +244,7 @@ Disable the long gap patching heuristic. When this option is applied, the
 maximum alignment gap is mostly controlled by
 .BR -r .
 .TP
-.B --lj-min-ratio \ FLOAT
+.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
@@ -252,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
@@ -357,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

misc.c

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

misc/paftools.js

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

mmpriv.h

@@ -28,6 +28,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(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs);
 void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs, mm128_t *a);
-void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r);
+void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r);
 void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr);
 
 void mm_est_err(const mm_idx_t *mi, int qlen, int n_regs, mm_reg1_t *regs, const mm128_t *a, int32_t n, const uint64_t *mini_pos);
@@ -86,7 +89,14 @@ mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int
 void mm_seg_free(void *km, int n_segs, mm_seg_t *segs);
 void mm_pair(void *km, int max_gap_ref, int dp_bonus, int sub_diff, int match_sc, const int *qlens, int *n_regs, mm_reg1_t **regs);
 
+FILE *mm_split_init(const char *prefix, const mm_idx_t *mi);
+mm_idx_t *mm_split_merge_prep(const char *prefix, int n_splits, FILE **fp, uint32_t *n_seq_part);
+int mm_split_merge(int n_segs, const char **fn, const mm_mapopt_t *opt, int n_split_idx);
+void mm_split_rm_tmp(const char *prefix, int n_splits);
+
 void mm_err_puts(const char *str);
+void mm_err_fwrite(const void *p, size_t size, size_t nitems, FILE *fp);
+void mm_err_fread(void *p, size_t size, size_t nitems, FILE *fp);
 
 #ifdef __cplusplus
 }

options.c

@@ -49,7 +49,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 
 void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)
 {
-	if ((opt->flag & MM_F_SPLICE_FOR) && (opt->flag & MM_F_SPLICE_REV))
+	if ((opt->flag & MM_F_SPLICE_FOR) || (opt->flag & MM_F_SPLICE_REV))
 		opt->flag |= MM_F_SPLICE;
 	if (opt->mid_occ <= 0)
 		opt->mid_occ = mm_idx_cal_max_occ(mi, opt->mid_occ_frac);
@@ -132,6 +132,11 @@ 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");
@@ -169,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;
 }

python/README.rst

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

python/cmappy.h

@@ -73,13 +73,17 @@ static inline void mm_reset_timer(void)
 extern unsigned char seq_comp_table[256];
 static inline mm_reg1_t *mm_map_aux(const mm_idx_t *mi, const char *seq1, const char *seq2, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt)
 {
+	mm_reg1_t *r;
+
+	Py_BEGIN_ALLOW_THREADS
 	if (seq2 == 0) {
-		return mm_map(mi, strlen(seq1), seq1, n_regs, b, opt, NULL);
+		r = mm_map(mi, strlen(seq1), seq1, n_regs, b, opt, NULL);
 	} else {
 		int _n_regs[2];
 		mm_reg1_t *regs[2];
 		char *seq[2];
 		int i, len[2];
+
 		len[0] = strlen(seq1);
 		len[1] = strlen(seq2);
 		seq[0] = (char*)seq1;
@@ -97,8 +101,11 @@ static inline mm_reg1_t *mm_map_aux(const mm_idx_t *mi, const char *seq1, const
 		regs[0] = (mm_reg1_t*)realloc(regs[0], sizeof(mm_reg1_t) * (*n_regs));
 		memcpy(&regs[0][_n_regs[0]], regs[1], _n_regs[1] * sizeof(mm_reg1_t));
 		free(regs[1]);
-		return regs[0];
+		r = regs[0];
 	}
+	Py_END_ALLOW_THREADS
+
+	return r;
 }
 
 static inline char *mappy_revcomp(int len, const uint8_t *seq)
@@ -119,8 +126,8 @@ static char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int e
 	*len = 0;
 	rid = mm_idx_name2id(mi, name);
 	if (rid < 0) return 0;
-	if (st >= mi->seq[rid].len || st >= en) return 0;
-	if (en < 0 || en > mi->seq[rid].len)
+	if ((uint32_t)st >= mi->seq[rid].len || st >= en) return 0;
+	if (en < 0 || (uint32_t)en > mi->seq[rid].len)
 		en = mi->seq[rid].len;
 	s = (char*)malloc(en - st + 1);
 	*len = mm_idx_getseq(mi, rid, st, en, (uint8_t*)s);
@@ -130,4 +137,16 @@ static char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int e
 	return s;
 }
 
+static mm_idx_t *mappy_idx_seq(int w, int k, int is_hpc, int bucket_bits, const char *seq, int len)
+{
+	const char *fake_name = "N/A";
+	char *s;
+	mm_idx_t *mi;
+	s = (char*)calloc(len + 1, 1);
+	memcpy(s, seq, len);
+	mi = mm_idx_str(w, k, is_hpc, bucket_bits, 1, (const char**)&s, (const char**)&fake_name);
+	free(s);
+	return mi;
+}
+
 #endif

python/cmappy.pxd

@@ -40,6 +40,7 @@ cdef extern from "minimap.h":
 		int32_t mid_occ
 		int32_t max_occ
 		int mini_batch_size
+		const char *split_prefix
 
 	int mm_set_opt(char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 	int mm_verbose
@@ -86,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)
@@ -106,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

python/mappy.pyx

@@ -3,7 +3,7 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
-__version__ = '2.11'
+__version__ = '2.13'
 
 cmappy.mm_reset_timer()
 
@@ -14,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
@@ -26,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
@@ -72,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))
@@ -85,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
@@ -102,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, max_frag_len=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
@@ -116,17 +126,32 @@ cdef class Aligner:
 		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:
@@ -135,18 +160,24 @@ cdef class Aligner:
 	def __bool__(self):
 		return (self._idx != NULL)
 
-	def map(self, seq, seq2=None, buf=None, max_frag_len=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:
@@ -157,13 +188,21 @@ cdef class Aligner:
 
 		for i in range(n_regs):
 			cmappy.mm_reg2hitpy(self._idx, &regs[i], &h)
-			cigar = []
-			for k in range(h.n_cigar32):
+			cigar, _cs, _MD = [], '', ''
+			for k in range(h.n_cigar32): # convert the 32-bit CIGAR encoding to Python array
 				c = h.cigar32[k]
 				cigar.append([c>>4, c&0xf])
-			yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.mlen, h.blen, h.NM, h.trans_strand, h.seg_id)
+			if cs or MD: # generate the cs and/or the MD tag, if requested
+				if cs:
+					l_cs_str = cmappy.mm_gen_cs(km, &cs_str, &m_cs_str, self._idx, &regs[i], _seq, 1)
+					_cs = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
+				if MD:
+					l_cs_str = cmappy.mm_gen_MD(km, &cs_str, &m_cs_str, self._idx, &regs[i], _seq)
+					_MD = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
+			yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.mlen, h.blen, h.NM, h.trans_strand, h.seg_id, _cs, _MD)
 			cmappy.mm_free_reg1(&regs[i])
 		free(regs)
+		free(cs_str)
 
 	def seq(self, str name, int start=0, int end=0x7fffffff):
 		cdef int l

python/minimap2.py

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

sdust.c

@@ -177,7 +177,7 @@ uint64_t *sdust(void *km, const uint8_t *seq, int l_seq, int T, int W, int *n)
 #ifdef _SDUST_MAIN
 #include <zlib.h>
 #include <stdio.h>
-#include "getopt.h"
+#include "ketopt.h"
 #include "kseq.h"
 KSEQ_INIT(gzFile, gzread)
 
@@ -186,16 +186,17 @@ int main(int argc, char *argv[])
 	gzFile fp;
 	kseq_t *ks;
 	int W = 64, T = 20, c;
+	ketopt_t o = KETOPT_INIT;
 
-	while ((c = getopt(argc, argv, "w:t:")) >= 0) {
-		if (c == 'w') W = atoi(optarg);
-		else if (c == 't') T = atoi(optarg);
+	while ((c = ketopt(&o, argc, argv, 1, "w:t:", 0)) >= 0) {
+		if (c == 'w') W = atoi(o.arg);
+		else if (c == 't') T = atoi(o.arg);
 	}
-	if (optind == argc) {
+	if (o.ind == argc) {
 		fprintf(stderr, "Usage: sdust [-w %d] [-t %d] <in.fa>\n", W, T);
 		return 1;
 	}
-	fp = strcmp(argv[optind], "-")? gzopen(argv[optind], "r") : gzdopen(fileno(stdin), "r");
+	fp = strcmp(argv[o.ind], "-")? gzopen(argv[o.ind], "r") : gzdopen(fileno(stdin), "r");
 	ks = kseq_init(fp);
 	while (kseq_read(ks) >= 0) {
 		uint64_t *r;

setup.py

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

splitidx.c

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