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

@@ -32,8 +32,8 @@ all:$(PROG)
 
 extra:all $(PROG_EXTRA)
 
-minimap2:main.o getopt.o libminimap2.a
+minimap2:main.o libminimap2.a
-		$(CC) $(CFLAGS) main.o getopt.o -o $@ -L. -lminimap2 $(LIBS)
+		$(CC) $(CFLAGS) main.o -o $@ -L. -lminimap2 $(LIBS)
 
 minimap2-lite:example.o libminimap2.a
 		$(CC) $(CFLAGS) $< -o $@ -L. -lminimap2 $(LIBS)
@@ -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
+sdust:sdust.c kalloc.o kalloc.h kdq.h kvec.h kseq.h ketopt.h sdust.h
-		$(CC) -D_SDUST_MAIN $(CFLAGS) $< getopt.o kalloc.o -o $@ -lz
+		$(CC) -D_SDUST_MAIN $(CFLAGS) $< kalloc.o -o $@ -lz
 
 # SSE-specific targets on x86/x86_64
 
@@ -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 ref.fa rna-reads.fa > aln.sam       # spliced long reads (strand unknown)
-./minimap2 -ax splice -k14 -uf ref.fa reads.fa > aln.sam  # Nanopore Direct RNA-seq
+./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam  # noisy Nanopore Direct RNA-seq
+./minimap2 -ax splice -uf -C5 ref.fa query.fa > aln.sam   # Final PacBio Iso-seq or traditional cDNA
 ./minimap2 -cx asm5 asm1.fa asm2.fa > aln.paf             # intra-species asm-to-asm alignment
 ./minimap2 -x ava-pb reads.fa reads.fa > overlaps.paf     # PacBio read overlap
 ./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf    # Nanopore read overlap
@@ -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 -
+curl -L https://github.com/lh3/minimap2/releases/download/v2.13/minimap2-2.13_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.11_x64-linux/minimap2
+./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

@@ -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) {
+		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 {
+			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);

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
+minimap2 release page at version tag [v2.10][v2.10]. Some examples only work
-with v2.11 or later.
+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 -
+curl -L https://github.com/lh3/minimap2/releases/download/v2.13/minimap2-2.13_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.11_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
+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);
+	if (is_MD) write_MD_core(s, tseq, qseq, r, tmp, write_tag);
-	else write_cs_core(s, tseq, qseq, r, tmp, no_iden);
+	else write_cs_core(s, tseq, qseq, r, tmp, no_iden, write_tag);
 	kfree(km, qseq); kfree(km, tseq); kfree(km, tmp);
 }
 
+int mm_gen_cs_or_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int is_MD, int no_iden)
+{
+	mm_bseq1_t t;
+	kstring_t str;
+	str.s = *buf, str.l = 0, str.m = *max_len;
+	t.l_seq = strlen(seq);
+	t.seq = (char*)seq;
+	write_cs_or_MD(km, &str, mi, &t, r, no_iden, is_MD, 0);
+	*max_len = str.m;
+	*buf = str.s;
+	return str.l;
+}
+
+int mm_gen_cs(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int no_iden)
+{
+	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 0, no_iden);
+}
+
+int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq)
+{
+	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 1, 0);
+}
+
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 {
 	int type;
@@ -277,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);
 }
@@ -394,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
+			int slen;
-			else if (flag & 0x800) mm_sprintf_lite(s, "%dS", r->re - r->rs); // supplementary alignment
+			if ((flag & 0x900) == 0 || (opt_flag & MM_F_SOFTCLIP)) slen = t->l_seq;
-			else mm_sprintf_lite(s, "%dS", t->l_seq);
+			else if (flag & 0x100) slen = 0;
+			else slen = r->qe - r->qs;
+			mm_sprintf_lite(s, "%dS%dN", slen, r->re - r->rs);
 		} else write_sam_cigar(s, flag, 0, t->l_seq, r, opt_flag);
 	}
 
@@ -476,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;

index.c

@@ -372,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]);
@@ -383,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);
@@ -512,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;

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 "ketopt.h"
-#include <getopt.h>
-#else
-#include "getopt.h"
-#endif
 
-#define MM_VERSION "2.11-r815-dirty"
+#define MM_VERSION "2.13-r850"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -26,51 +22,53 @@ void liftrlimit()
 void liftrlimit() {}
 #endif
 
-static struct option long_options[] = {
+static ko_longopt_t long_options[] = {
-	{ "bucket-bits",    required_argument, 0, 0 },
+	{ "bucket-bits",    ko_required_argument, 300 },
-	{ "mb-size",        required_argument, 0, 'K' },
+	{ "mb-size",        ko_required_argument, 'K' },
-	{ "seed",           required_argument, 0, 0 },
+	{ "seed",           ko_required_argument, 302 },
-	{ "no-kalloc",      no_argument,       0, 0 },
+	{ "no-kalloc",      ko_no_argument,       303 },
-	{ "print-qname",    no_argument,       0, 0 },
+	{ "print-qname",    ko_no_argument,       304 },
-	{ "no-self",        no_argument,       0, 'D' },
+	{ "no-self",        ko_no_argument,       'D' },
-	{ "print-seeds",    no_argument,       0, 0 },
+	{ "print-seeds",    ko_no_argument,       306 },
-	{ "max-chain-skip", required_argument, 0, 0 },
+	{ "max-chain-skip", ko_required_argument, 307 },
-	{ "min-dp-len",     required_argument, 0, 0 },
+	{ "min-dp-len",     ko_required_argument, 308 },
-	{ "print-aln-seq",  no_argument,       0, 0 },
+	{ "print-aln-seq",  ko_no_argument,       309 },
-	{ "splice",         no_argument,       0, 0 },
+	{ "splice",         ko_no_argument,       310 },
-	{ "cost-non-gt-ag", required_argument, 0, 'C' },
+	{ "cost-non-gt-ag", ko_required_argument, 'C' },
-	{ "no-long-join",   no_argument,       0, 0 },
+	{ "no-long-join",   ko_no_argument,       312 },
-	{ "sr",             no_argument,       0, 0 },
+	{ "sr",             ko_no_argument,       313 },
-	{ "frag",           required_argument, 0, 0 },
+	{ "frag",           ko_required_argument, 314 },
-	{ "secondary",      required_argument, 0, 0 },
+	{ "secondary",      ko_required_argument, 315 },
-	{ "cs",             optional_argument, 0, 0 },
+	{ "cs",             ko_optional_argument, 316 },
-	{ "end-bonus",      required_argument, 0, 0 },
+	{ "end-bonus",      ko_required_argument, 317 },
-	{ "no-pairing",     no_argument,       0, 0 },
+	{ "no-pairing",     ko_no_argument,       318 },
-	{ "splice-flank",   required_argument, 0, 0 },
+	{ "splice-flank",   ko_required_argument, 319 },
-	{ "idx-no-seq",     no_argument,       0, 0 },
+	{ "idx-no-seq",     ko_no_argument,       320 },
-	{ "end-seed-pen",   required_argument, 0, 0 },   // 21
+	{ "end-seed-pen",   ko_required_argument, 321 },
-	{ "for-only",       no_argument,       0, 0 },   // 22
+	{ "for-only",       ko_no_argument,       322 },
-	{ "rev-only",       no_argument,       0, 0 },   // 23
+	{ "rev-only",       ko_no_argument,       323 },
-	{ "heap-sort",      required_argument, 0, 0 },   // 24
+	{ "heap-sort",      ko_required_argument, 324 },
-	{ "all-chain",      no_argument,       0, 'P' },
+	{ "all-chain",      ko_no_argument,       'P' },
-	{ "dual",           required_argument, 0, 0 },   // 26
+	{ "dual",           ko_required_argument, 326 },
-	{ "max-clip-ratio", required_argument, 0, 0 },   // 27
+	{ "max-clip-ratio", ko_required_argument, 327 },
-	{ "min-occ-floor",  required_argument, 0, 0 },   // 28
+	{ "min-occ-floor",  ko_required_argument, 328 },
-	{ "MD",             no_argument,       0, 0 },   // 29
+	{ "MD",             ko_no_argument,       329 },
-	{ "lj-min-ratio",   required_argument, 0, 0 },   // 30
+	{ "lj-min-ratio",   ko_required_argument, 330 },
-	{ "score-N",        required_argument, 0, 0 },   // 31
+	{ "score-N",        ko_required_argument, 331 },
-	{ "eqx",            no_argument,       0, 0 },   // 32
+	{ "eqx",            ko_no_argument,       332 },
-	{ "paf-no-hit",     no_argument,       0, 0 },   // 33
+	{ "paf-no-hit",     ko_no_argument,       333 },
-	{ "split-prefix",   required_argument, 0, 0 },   // 34
+	{ "split-prefix",   ko_required_argument, 334 },
-	{ "help",           no_argument,       0, 'h' },
+	{ "no-end-flt",     ko_no_argument,       335 },
-	{ "max-intron-len", required_argument, 0, 'G' },
+	{ "hard-mask-level",ko_no_argument,       336 },
-	{ "version",        no_argument,       0, 'V' },
+	{ "help",           ko_no_argument,       'h' },
-	{ "min-count",      required_argument, 0, 'n' },
+	{ "max-intron-len", ko_required_argument, 'G' },
-	{ "min-chain-score",required_argument, 0, 'm' },
+	{ "version",        ko_no_argument,       'V' },
-	{ "mask-level",     required_argument, 0, 'M' },
+	{ "min-count",      ko_required_argument, 'n' },
-	{ "min-dp-score",   required_argument, 0, 's' },
+	{ "min-chain-score",ko_required_argument, 'm' },
-	{ "sam",            no_argument,       0, 'a' },
+	{ "mask-level",     ko_required_argument, 'M' },
-	{ 0, 0, 0, 0}
+	{ "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)
@@ -99,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, n_parts, 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;
@@ -113,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) {
+			if (mm_set_opt(o.arg, &ipt, &opt) < 0) {
-				fprintf(stderr, "[ERROR] unknown preset '%s'\n", optarg);
+				fprintf(stderr, "[ERROR] unknown preset '%s'\n", o.arg);
 				return 1;
 			}
-			break;
+		} else if (c == ':') {
+			fprintf(stderr, "[ERROR] missing option argument\n");
+			return 1;
+		} else if (c == '?') {
+			fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[o.i]);
+			return 1;
 		}
-	optind = 0; // for musl getopt, optind=0 has the same effect as optreset=1; older libc doesn't have optreset
+	}
+	o = KETOPT_INIT;
 
-	while ((c = getopt_long(argc, argv, opt_str, long_options, &long_idx)) >= 0) {
+	while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) {
-		if (c == 'w') ipt.w = atoi(optarg);
+		if (c == 'w') ipt.w = atoi(o.arg);
-		else if (c == 'k') ipt.k = atoi(optarg);
+		else if (c == 'k') ipt.k = atoi(o.arg);
 		else if (c == 'H') ipt.flag |= MM_I_HPC;
-		else if (c == 'd') fnw = optarg; // the above are indexing related options, except -I
+		else if (c == 'd') fnw = o.arg; // the above are indexing related options, except -I
-		else if (c == 'r') opt.bw = (int)mm_parse_num(optarg);
+		else if (c == 'r') opt.bw = (int)mm_parse_num(o.arg);
-		else if (c == 't') n_threads = atoi(optarg);
+		else if (c == 't') n_threads = atoi(o.arg);
-		else if (c == 'v') mm_verbose = atoi(optarg);
+		else if (c == 'v') mm_verbose = atoi(o.arg);
-		else if (c == 'g') opt.max_gap = (int)mm_parse_num(optarg);
+		else if (c == 'g') opt.max_gap = (int)mm_parse_num(o.arg);
-		else if (c == 'G') mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(optarg));
+		else if (c == 'G') mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(o.arg));
-		else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(optarg);
+		else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(o.arg);
-		else if (c == 'N') opt.best_n = atoi(optarg);
+		else if (c == 'N') old_best_n = opt.best_n, opt.best_n = atoi(o.arg);
-		else if (c == 'p') opt.pri_ratio = atof(optarg);
+		else if (c == 'p') opt.pri_ratio = atof(o.arg);
-		else if (c == 'M') opt.mask_level = atof(optarg);
+		else if (c == 'M') opt.mask_level = atof(o.arg);
 		else if (c == 'c') opt.flag |= MM_F_OUT_CG | MM_F_CIGAR;
 		else if (c == 'D') opt.flag |= MM_F_NO_DIAG;
 		else if (c == 'P') opt.flag |= MM_F_ALL_CHAINS;
@@ -146,64 +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 == 'T') opt.sdust_thres = atoi(o.arg);
-		else if (c == 'n') opt.min_cnt = atoi(optarg);
+		else if (c == 'n') opt.min_cnt = atoi(o.arg);
-		else if (c == 'm') opt.min_chain_score = atoi(optarg);
+		else if (c == 'm') opt.min_chain_score = atoi(o.arg);
-		else if (c == 'A') opt.a = atoi(optarg);
+		else if (c == 'A') opt.a = atoi(o.arg);
-		else if (c == 'B') opt.b = atoi(optarg);
+		else if (c == 'B') opt.b = atoi(o.arg);
-		else if (c == 's') opt.min_dp_max = atoi(optarg);
+		else if (c == 's') opt.min_dp_max = atoi(o.arg);
-		else if (c == 'C') opt.noncan = atoi(optarg);
+		else if (c == 'C') opt.noncan = atoi(o.arg);
-		else if (c == 'I') ipt.batch_size = mm_parse_num(optarg);
+		else if (c == 'I') ipt.batch_size = mm_parse_num(o.arg);
-		else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(optarg);
+		else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(o.arg);
-		else if (c == 'R') rg = optarg;
+		else if (c == 'R') rg = o.arg;
 		else if (c == 'h') fp_help = stdout;
 		else if (c == '2') opt.flag |= MM_F_2_IO_THREADS;
-		else if (c == 0 && long_idx == 0) ipt.bucket_bits = atoi(optarg); // --bucket-bits
+		else if (c == 300) ipt.bucket_bits = atoi(o.arg); // --bucket-bits
-		else if (c == 0 && long_idx == 2) opt.seed = atoi(optarg); // --seed
+		else if (c == 302) opt.seed = atoi(o.arg); // --seed
-		else if (c == 0 && long_idx == 3) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
+		else if (c == 303) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
-		else if (c == 0 && long_idx == 4) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
+		else if (c == 304) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
-		else if (c == 0 && long_idx == 6) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
+		else if (c == 306) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
-		else if (c == 0 && long_idx == 7) opt.max_chain_skip = atoi(optarg); // --max-chain-skip
+		else if (c == 307) opt.max_chain_skip = atoi(o.arg); // --max-chain-skip
-		else if (c == 0 && long_idx == 8) opt.min_ksw_len = atoi(optarg); // --min-dp-len
+		else if (c == 308) opt.min_ksw_len = atoi(o.arg); // --min-dp-len
-		else if (c == 0 && long_idx == 9) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
+		else if (c == 309) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
-		else if (c == 0 && long_idx ==10) opt.flag |= MM_F_SPLICE; // --splice
+		else if (c == 310) opt.flag |= MM_F_SPLICE; // --splice
-		else if (c == 0 && long_idx ==12) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
+		else if (c == 312) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
-		else if (c == 0 && long_idx ==13) opt.flag |= MM_F_SR; // --sr
+		else if (c == 313) opt.flag |= MM_F_SR; // --sr
-		else if (c == 0 && long_idx ==17) opt.end_bonus = atoi(optarg); // --end-bonus
+		else if (c == 317) opt.end_bonus = atoi(o.arg); // --end-bonus
-		else if (c == 0 && long_idx ==18) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
+		else if (c == 318) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
-		else if (c == 0 && long_idx ==20) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
+		else if (c == 320) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
-		else if (c == 0 && long_idx ==21) opt.anchor_ext_shift = atoi(optarg); // --end-seed-pen
+		else if (c == 321) opt.anchor_ext_shift = atoi(o.arg); // --end-seed-pen
-		else if (c == 0 && long_idx ==22) opt.flag |= MM_F_FOR_ONLY; // --for-only
+		else if (c == 322) opt.flag |= MM_F_FOR_ONLY; // --for-only
-		else if (c == 0 && long_idx ==23) opt.flag |= MM_F_REV_ONLY; // --rev-only
+		else if (c == 323) opt.flag |= MM_F_REV_ONLY; // --rev-only
-		else if (c == 0 && long_idx ==27) opt.max_clip_ratio = atof(optarg); // --max-clip-ratio
+		else if (c == 327) opt.max_clip_ratio = atof(o.arg); // --max-clip-ratio
-		else if (c == 0 && long_idx ==28) opt.min_mid_occ = atoi(optarg); // --min-occ-floor
+		else if (c == 328) opt.min_mid_occ = atoi(o.arg); // --min-occ-floor
-		else if (c == 0 && long_idx ==29) opt.flag |= MM_F_OUT_MD; // --MD
+		else if (c == 329) opt.flag |= MM_F_OUT_MD; // --MD
-		else if (c == 0 && long_idx ==30) opt.min_join_flank_ratio = atof(optarg); // --lj-min-ratio
+		else if (c == 330) opt.min_join_flank_ratio = atof(o.arg); // --lj-min-ratio
-		else if (c == 0 && long_idx ==31) opt.sc_ambi = atoi(optarg); // --score-N
+		else if (c == 331) opt.sc_ambi = atoi(o.arg); // --score-N
-		else if (c == 0 && long_idx ==32) opt.flag |= MM_F_EQX; // --eqx
+		else if (c == 332) opt.flag |= MM_F_EQX; // --eqx
-		else if (c == 0 && long_idx ==33) opt.flag |= MM_F_PAF_NO_HIT; // --paf-no-hit
+		else if (c == 333) opt.flag |= MM_F_PAF_NO_HIT; // --paf-no-hit
-		else if (c == 0 && long_idx ==34) opt.split_prefix = optarg; // --split-prefix
+		else if (c == 334) opt.split_prefix = o.arg; // --split-prefix
-		else if (c == 0 && long_idx == 14) { // --frag
+		else if (c == 335) opt.flag |= MM_F_NO_END_FLT; // --no-end-flt
-			yes_or_no(&opt, MM_F_FRAG_MODE, long_idx, optarg, 1);
+		else if (c == 336) opt.flag |= MM_F_HARD_MLEVEL; // --hard-mask-level
-		} else if (c == 0 && long_idx == 15) { // --secondary
+		else if (c == 314) { // --frag
-			yes_or_no(&opt, MM_F_NO_PRINT_2ND, long_idx, optarg, 0);
+			yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1);
-		} else if (c == 0 && long_idx == 16) { // --cs
+		} else if (c == 315) { // --secondary
+			yes_or_no(&opt, MM_F_NO_PRINT_2ND, o.longidx, o.arg, 0);
+		} else if (c == 316) { // --cs
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR;
-			if (optarg == 0 || strcmp(optarg, "short") == 0) {
+			if (o.arg == 0 || strcmp(o.arg, "short") == 0) {
 				opt.flag &= ~MM_F_OUT_CS_LONG;
-			} else if (strcmp(optarg, "long") == 0) {
+			} else if (strcmp(o.arg, "long") == 0) {
 				opt.flag |= MM_F_OUT_CS_LONG;
-			} else if (strcmp(optarg, "none") == 0) {
+			} else if (strcmp(o.arg, "none") == 0) {
 				opt.flag &= ~MM_F_OUT_CS;
 			} else if (mm_verbose >= 2) {
 				fprintf(stderr, "[WARNING]\033[1;31m --cs only takes 'short' or 'long'. Invalid values are assumed to be 'short'.\033[0m\n");
 			}
-		} else if (c == 0 && long_idx == 19) { // --splice-flank
+		} else if (c == 319) { // --splice-flank
-			yes_or_no(&opt, MM_F_SPLICE_FLANK, long_idx, optarg, 1);
+			yes_or_no(&opt, MM_F_SPLICE_FLANK, o.longidx, o.arg, 1);
-		} else if (c == 0 && long_idx == 24) { // --heap-sort
+		} else if (c == 324) { // --heap-sort
-			yes_or_no(&opt, MM_F_HEAP_SORT, long_idx, optarg, 1);
+			yes_or_no(&opt, MM_F_HEAP_SORT, o.longidx, o.arg, 1);
-		} else if (c == 0 && long_idx == 26) { // --dual
+		} else if (c == 326) { // --dual
-			yes_or_no(&opt, MM_F_NO_DUAL, long_idx, optarg, 0);
+			yes_or_no(&opt, MM_F_NO_DUAL, o.longidx, o.arg, 0);
 		} else if (c == 'S') {
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR | MM_F_OUT_CS_LONG;
 			if (mm_verbose >= 2)
@@ -214,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
+			if (*o.arg == 'b') opt.flag |= MM_F_SPLICE_FOR|MM_F_SPLICE_REV; // both strands
-			else if (*optarg == 'f') opt.flag |= MM_F_SPLICE_FOR, opt.flag &= ~MM_F_SPLICE_REV; // match GT-AG
+			else if (*o.arg == 'f') opt.flag |= MM_F_SPLICE_FOR, opt.flag &= ~MM_F_SPLICE_REV; // match GT-AG
-			else if (*optarg == 'r') opt.flag |= MM_F_SPLICE_REV, opt.flag &= ~MM_F_SPLICE_FOR; // match CT-AC (reverse complement of GT-AG)
+			else if (*o.arg == 'r') opt.flag |= MM_F_SPLICE_REV, opt.flag &= ~MM_F_SPLICE_FOR; // match CT-AC (reverse complement of GT-AG)
-			else if (*optarg == 'n') opt.flag &= ~(MM_F_SPLICE_FOR|MM_F_SPLICE_REV); // don't try to match the GT-AG signal
+			else if (*o.arg == 'n') opt.flag &= ~(MM_F_SPLICE_FOR|MM_F_SPLICE_REV); // don't try to match the GT-AG signal
 			else {
 				fprintf(stderr, "[ERROR]\033[1;31m unrecognized cDNA direction\033[0m\n");
 				return 1;
 			}
 		} else if (c == 'z') {
-			opt.zdrop = opt.zdrop_inv = strtol(optarg, &s, 10);
+			opt.zdrop = opt.zdrop_inv = strtol(o.arg, &s, 10);
 			if (*s == ',') opt.zdrop_inv = strtol(s + 1, &s, 10);
 		} else if (c == 'O') {
-			opt.q = opt.q2 = strtol(optarg, &s, 10);
+			opt.q = opt.q2 = strtol(o.arg, &s, 10);
 			if (*s == ',') opt.q2 = strtol(s + 1, &s, 10);
 		} else if (c == 'E') {
-			opt.e = opt.e2 = strtol(optarg, &s, 10);
+			opt.e = opt.e2 = strtol(o.arg, &s, 10);
 			if (*s == ',') opt.e2 = strtol(s + 1, &s, 10);
 		}
 	}
@@ -246,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");
@@ -297,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;
@@ -336,13 +347,13 @@ int main(int argc, char *argv[])
 		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);
 	}
@@ -350,7 +361,7 @@ int main(int argc, char *argv[])
 	mm_idx_reader_close(idx_rdr);
 
 	if (opt.split_prefix)
-		mm_split_merge(argc - (optind + 1), (const char**)&argv[optind + 1], &opt, n_parts);
+		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");
@@ -362,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

@@ -29,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;
@@ -243,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
@@ -256,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);
 	}
@@ -354,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);
 		}
@@ -497,7 +502,7 @@ static void merge_hits(step_t *s)
 				}
 			}
 			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);
+			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]);
@@ -596,6 +601,7 @@ 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, n_seg, rep_len and frag_gap were allocated with reg; no memory leak here
@@ -682,7 +688,7 @@ int mm_split_merge(int n_segs, const char **fn, const mm_mapopt_t *opt, int n_sp
 	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 < pl.mi->n_seq; ++i)
+		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);

minimap.h

@@ -33,6 +33,8 @@
 #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
@@ -301,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
  *
@@ -337,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
@@ -360,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/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = 'r767';
+var paftools_version = '2.13-r850';
 
 /*****************************
  ***** Library functions *****
@@ -341,14 +341,15 @@ 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, gap_thres_long = 1000, min_mapq = 5;
-	var fa_tmp = null, fa, fa_lens, is_vcf = false;
+	var fa_tmp = null, fa, fa_lens, is_vcf = false, sample_name = "sample";
-	while ((c = getopt(args, "l:L:g:q:B:f:")) != null) {
+	while ((c = getopt(args, "l:L:g:q:B:f:s:")) != null) {
 		if (c == 'l') min_cov_len = parseInt(getopt.arg);
 		else if (c == 'L') min_var_len = parseInt(getopt.arg);
 		else if (c == 'g') gap_thres = parseInt(getopt.arg);
 		else if (c == 'G') gap_thres_long = parseInt(getopt.arg);
 		else if (c == 'q') min_mapq = parseInt(getopt.arg);
 		else if (c == 'f') fa_tmp = fasta_read(getopt.arg, fa_lens);
+		else if (c == 's') sample_name = getopt.arg;
 	}
 	if (fa_tmp != null) fa = fa_tmp[0], fa_lens = fa_tmp[1], is_vcf = true;
 
@@ -360,6 +361,7 @@ 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);
 	}
 
@@ -423,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 = [];

mmpriv.h

@@ -75,7 +75,7 @@ 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);

options.c

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

@@ -126,8 +126,8 @@ static char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int e
 	*len = 0;
 	rid = mm_idx_name2id(mi, name);
 	if (rid < 0) return 0;
-	if (st >= mi->seq[rid].len || st >= en) return 0;
+	if ((uint32_t)st >= mi->seq[rid].len || st >= en) return 0;
-	if (en < 0 || en > mi->seq[rid].len)
+	if (en < 0 || (uint32_t)en > mi->seq[rid].len)
 		en = mi->seq[rid].len;
 	s = (char*)malloc(en - st + 1);
 	*len = mm_idx_getseq(mi, rid, st, en, (uint8_t*)s);
@@ -137,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),
+		a = [str(self._q_st), str(self._q_en), strand, self._ctg, str(self._ctg_len), str(self._r_st), str(self._r_en),
-				str(self._mlen), str(self._blen), str(self._mapq), tp, ts, "cg:Z:" + self.cigar_str])
+			str(self._mlen), str(self._blen), str(self._mapq), tp, ts, "cg:Z:" + self.cigar_str]
+		if self._cs != "": a.append("cs:Z:" + self._cs)
+		return "\t".join(a)
 
 cdef class ThreadBuffer:
 	cdef cmappy.mm_tbuf_t *_b
@@ -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,8 +126,19 @@ 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 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:
@@ -127,6 +148,10 @@ cdef class Aligner:
 				cmappy.mm_idx_reader_close(r)
 				cmappy.mm_mapopt_update(&self.map_opt, self._idx)
 				cmappy.mm_idx_index_name(self._idx)
+		else:
+			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)
+			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 = []
+			cigar, _cs, _MD = [], '', ''
-			for k in range(h.n_cigar32):
+			for k in range(h.n_cigar32): # convert the 32-bit CIGAR encoding to Python array
 				c = h.cigar32[k]
 				cigar.append([c>>4, c&0xf])
-			yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.mlen, h.blen, h.NM, h.trans_strand, h.seg_id)
+			if cs or MD: # generate the cs and/or the MD tag, if requested
+				if cs:
+					l_cs_str = cmappy.mm_gen_cs(km, &cs_str, &m_cs_str, self._idx, &regs[i], _seq, 1)
+					_cs = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
+				if MD:
+					l_cs_str = cmappy.mm_gen_MD(km, &cs_str, &m_cs_str, self._idx, &regs[i], _seq)
+					_MD = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
+			yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.mlen, h.blen, h.NM, h.trans_strand, h.seg_id, _cs, _MD)
 			cmappy.mm_free_reg1(&regs[i])
 		free(regs)
+		free(cs_str)
 
 	def seq(self, str name, int start=0, int end=0x7fffffff):
 		cdef int l

python/minimap2.py

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

sdust.c

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

setup.py

@@ -14,16 +14,26 @@ else: # with Cython
 	module_src = 'python/mappy.pyx'
 	cmdclass['build_ext'] = build_ext
 
-import sys
+import sys, platform
+
 sys.path.append('python')
 
+extra_compile_args = ['-DHAVE_KALLOC']
+include_dirs = ["."]
+
+if platform.machine() in ["aarch64", "arm64"]:
+	include_dirs.append("sse2neon/")
+	extra_compile_args.extend(['-ftree-vectorize', '-DKSW_SSE2_ONLY', '-D__SSE2__'])
+else:
+	extra_compile_args.append('-msse4.1') # WARNING: ancient x86_64 CPUs don't have SSE4
+
 def readme():
 	with open('python/README.rst') as f:
 		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
+		extra_compile_args = extra_compile_args,
-		include_dirs = ['.'],
+		include_dirs = include_dirs,
 		libraries = ['z', 'm', 'pthread'])],
 	classifiers = [
 		'Development Status :: 5 - Production/Stable',