Release minimap2-2.31 (r1302)

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

.github/workflows/ci.yaml

@@ -35,7 +35,7 @@ jobs:
       uses: actions/checkout@v4
 
     - name: Compile with ${{ matrix.compiler }}
-      uses: uraimo/run-on-arch-action@v2
+      uses: uraimo/run-on-arch-action@v3
       with:
         arch: aarch64
         distro: ubuntu22.04
@@ -65,4 +65,4 @@ jobs:
       run: |
         make CC=${{ matrix.compiler }} arm_neon=1 aarch64=1 -j
         file minimap2 | grep arm64
-        
+        

Makefile

@@ -102,7 +102,7 @@ ksw2_exts2_neon.o:ksw2_exts2_sse.c ksw2.h kalloc.h
 # other non-file targets
 
 clean:
-		rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist mappy*.so mappy.c python/mappy.c mappy.egg*
+		rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist mappy*.so mappy.c python/mappy.c mappy.egg* .eggs
 
 depend:
 		(LC_ALL=C; export LC_ALL; makedepend -Y -- $(CFLAGS) $(CPPFLAGS) -- *.c)

NEWS.md

@@ -1,3 +1,51 @@
+Release 2.31-r1302 (19 May 2026)
+--------------------------------
+
+Notable changes to minimap2:
+
+ * Bugfix: supplementary and secondary alignments were occasionally flagged
+   incorrectly.
+
+ * Bugfix: Smith-Waterman alignment for inversion alignment led to an
+   out-of-bound access in rare cases.
+
+Changes to paftools.js:
+
+ * New feature: new `sim2bed` subcommand to get a BED file from simulated
+   reads.
+
+ * New feature: new `badread2fa` subcommand to format reads simulated by
+   the Badread simulator.
+
+Change to the python binding:
+
+ * New feature: mappy optionally writes the `ds` tag.
+
+ * Bugfix: a use-after-free error (#1345)
+
+The two bugs in minimap2 had existed for years. They were caught by Jeremy Wang
+at UNC when he ported minimap2 to Rust. Due to the two bug fixes, this version
+occasionally produces alignment different from the last version.
+
+(2.31: 19 May 2026, r1302)
+
+
+
+Release 2.30-r1287 (15 June 2025)
+---------------------------------
+
+Notable changes:
+
+ * Improvement: consolidated `--spsc`.
+
+ * Deprecation: subcommands `splice2bed`, `gff2bed`, `gff2junc`, `junceval` and
+   `exoneval` in `paftools.js` are deprecated by minigff. They will remain
+   indefinitely for backward compatibility.
+
+(2.30: 15 June 2025, r1287)
+
+
+
 Release 2.29-r1283 (18 April 2025)
 ----------------------------------
 

README.md

@@ -3,6 +3,7 @@
 [![PyPI](https://img.shields.io/pypi/v/mappy.svg?style=flat)](https://pypi.python.org/pypi/mappy)
 [![Build Status](https://github.com/lh3/minimap2/actions/workflows/ci.yaml/badge.svg)](https://github.com/lh3/minimap2/actions)
 ## <a name="started"></a>Getting Started
+**ALERT:** `minimap2.com` is a [phishing site](https://github.com/lh3/minimap2/issues/1316). Please don't use anything from that website.
 ```sh
 git clone https://github.com/lh3/minimap2
 cd minimap2 && make
@@ -77,8 +78,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.29/minimap2-2.29_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.29_x64-linux/minimap2
+curl -L https://github.com/lh3/minimap2/releases/download/v2.31/minimap2-2.31_x64-linux.tar.bz2 | tar -jxvf -
+./minimap2-2.31_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

cookbook.md

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

format.c

@@ -361,24 +361,34 @@ static void write_cs_ds_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const
 	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, int is_qstrand)
+int mm_gen_cs_ds_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 is_ds, int no_iden, int is_qstrand)
 {
 	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_ds_or_MD(km, &str, mi, &t, r, no_iden, is_MD, 0, 0, is_qstrand);
+	write_cs_ds_or_MD(km, &str, mi, &t, r, no_iden, is_MD, is_ds, 0, is_qstrand);
 	*max_len = str.m;
 	*buf = str.s;
 	return str.l;
 }
 
+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, int is_qstrand)
+{
+	return mm_gen_cs_ds_or_MD(km, buf, max_len, mi, r, seq, is_MD, 0, no_iden, is_qstrand);
+}
+
 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, 0);
 }
 
+int mm_gen_ds(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_ds_or_MD(km, buf, max_len, mi, r, seq, 0, 1, no_iden, 0);
+}
+
 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, 0);

hit.c

@@ -256,19 +256,25 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int chec
 {
 	if (pri_ratio > 0.0f && *n_ > 0) {
 		int i, k, n = *n_, n_2nd = 0;
-		for (i = k = 0; i < n; ++i) {
+		uint8_t *keep = (uint8_t*)kmalloc(km, n);
+		for (i = 0; i < n; ++i) {
 			int p = r[i].parent;
+			keep[i] = 0;
 			if (p == i || r[i].inv) { // primary or inversion
-				r[k++] = r[i];
+				keep[i] = 1;
 			} else if ((r[i].score >= r[p].score * pri_ratio || r[i].score + min_diff >= r[p].score) && n_2nd < best_n) {
 				if (!(r[i].qs == r[p].qs && r[i].qe == r[p].qe && r[i].rid == r[p].rid && r[i].rs == r[p].rs && r[i].re == r[p].re)) // not identical hits
-					r[k++] = r[i], ++n_2nd;
-				else if (r[i].p) free(r[i].p);
+					keep[i] = 1, ++n_2nd;
 			} else if (check_strand && n_2nd < best_n && r[i].score > min_strand_sc && r[i].rev != r[p].rev) {
 				r[i].strand_retained = 1;
-				r[k++] = r[i], ++n_2nd;
-			} else if (r[i].p) free(r[i].p);
+				keep[i] = 1, ++n_2nd;
+			}
 		}
+		for (i = k = 0; i < n; ++i) {
+			if (keep[i]) r[k++] = r[i];
+			else if (r[i].p) free(r[i].p);
+		}
+		kfree(km, keep);
 		if (k != n) mm_sync_regs(km, k, r); // removing hits requires sync()
 		*n_ = k;
 	}
@@ -277,13 +283,18 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int chec
 int mm_filter_strand_retained(int n_regs, mm_reg1_t *r)
 {
 	int i, k;
-	for (i = k = 0; i < n_regs; ++i) {
+	uint8_t *keep = (uint8_t*)malloc(n_regs);
+	for (i = 0; i < n_regs; ++i) {
 		int p = r[i].parent;
-		if (!r[i].strand_retained || r[i].div < r[p].div * 5.0f || r[i].div < 0.01f) {
+		keep[i] = (!r[i].strand_retained || r[i].div < r[p].div * 5.0f || r[i].div < 0.01f);
+	}
+	for (i = k = 0; i < n_regs; ++i) {
+		if (keep[i]) {
 			if (k < i) r[k++] = r[i];
 			else ++k;
 		}
 	}
+	free(keep);
 	return k;
 }
 

index.c

@@ -967,7 +967,7 @@ typedef struct mm_idx_spsc_s {
 	uint64_t *a; // pos<<56 | score<<1 | acceptor
 } mm_idx_spsc_t;
 
-int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc)
+int32_t mm_idx_spsc_read2(mm_idx_t *idx, const char *fn, int32_t max_sc, float scale)
 {
 	gzFile fp;
 	kstring_t str = {0,0,0};
@@ -1007,6 +1007,8 @@ int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc)
 			}
 		}
 		if (i < 4) continue; // not enough fields
+		if (scale > 0.0f && scale < 1.0f)
+			score = score > 0.0f? (int)(score * scale + .499) : (int)(score * scale - .499);
 		if (score > max_sc) score = max_sc;
 		if (score < -max_sc) score = -max_sc;
 		cid = mm_idx_name2id(idx, name);
@@ -1030,6 +1032,11 @@ int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc)
 	return 0;
 }
 
+int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc)
+{
+	return mm_idx_spsc_read2(idx, fn, max_sc, 1.0f);
+}
+
 static int32_t mm_idx_find_intv(int32_t n, const uint64_t *a, int64_t x)
 {
 	int32_t s = 0, e = n;

ksw2_ll_sse.c

@@ -67,7 +67,7 @@ void *ksw_ll_qinit(void *km, int size, int qlen, const uint8_t *query, int m, co
 			const int8_t *ma = mat + a * m;
 			for (i = 0; i < slen; ++i)
 				for (k = i; k < nlen; k += slen) // p iterations
-					*t++ = (k >= qlen? 0 : ma[query[k]]) + q->shift;
+					*t++ = (k >= qlen? -1 : ma[query[k]]) + q->shift;
 		}
 	} else {
 		int16_t *t = (int16_t*)q->qp;
@@ -76,7 +76,7 @@ void *ksw_ll_qinit(void *km, int size, int qlen, const uint8_t *query, int m, co
 			const int8_t *ma = mat + a * m;
 			for (i = 0; i < slen; ++i)
 				for (k = i; k < nlen; k += slen) // p iterations
-					*t++ = (k >= qlen? 0 : ma[query[k]]);
+					*t++ = (k >= qlen? -1 : ma[query[k]]);
 		}
 	}
 	return q;

main.c

@@ -85,6 +85,8 @@ static ko_longopt_t long_options[] = {
 	{ "jump-min-match", ko_required_argument, 360 },
 	{ "write-junc",     ko_no_argument,       361 },
 	{ "pass1",          ko_required_argument, 362 },
+	{ "spsc-scale",     ko_required_argument, 363 },
+	{ "spsc0",          ko_required_argument, 364 },
 	{ "dbg-seed-occ",   ko_no_argument,       501 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
@@ -134,6 +136,7 @@ int main(int argc, char *argv[])
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
 	int i, c, n_threads = 3, n_parts, old_best_n = -1;
+	float spsc_scale = 0.7f;
 	char *fnw = 0, *rg = 0, *fn_bed_junc = 0, *fn_bed_jump = 0, *fn_bed_pass1 = 0, *fn_spsc = 0, *s, *alt_list = 0;
 	FILE *fp_help = stderr;
 	mm_idx_reader_t *idx_rdr;
@@ -240,7 +243,6 @@ int main(int argc, char *argv[])
 		else if (c == 338) opt.max_qlen = mm_parse_num(o.arg); // --max-qlen
 		else if (c == 340) fn_bed_junc = o.arg; // --junc-bed
 		else if (c == 341) opt.junc_bonus = atoi(o.arg); // --junc-bonus
-		else if (c == 358) opt.junc_pen = atoi(o.arg); // --junc-pen
 		else if (c == 342) opt.flag |= MM_F_SAM_HIT_ONLY; // --sam-hit-only
 		else if (c == 343) opt.chain_gap_scale = atof(o.arg); // --chain-gap-scale
 		else if (c == 351) opt.chain_skip_scale = atof(o.arg); // --chain-skip-scale
@@ -260,6 +262,8 @@ int main(int argc, char *argv[])
 		else if (c == 361) opt.flag |= MM_F_OUT_JUNC | MM_F_CIGAR; // --write-junc
 		else if (c == 362) fn_bed_pass1 = o.arg; // --jump-pass1
 		else if (c == 501) mm_dbg_flag |= MM_DBG_SEED_FREQ; // --dbg-seed-occ
+		else if (c == 363) spsc_scale = atof(o.arg); // --spsc-scale
+		else if (c == 358 || c == 364) opt.junc_pen = atoi(o.arg); // --junc-pen or --spsc0
 		else if (c == 330) {
 			fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n");
 		} else if (c == 313) { // --sr
@@ -476,7 +480,7 @@ int main(int argc, char *argv[])
 				fprintf(stderr, "[WARNING] failed to load the pass-1 jump BED file\n");
 		}
 		if (fn_spsc) {
-			mm_idx_spsc_read(mi, fn_spsc, mm_max_spsc_bonus(&opt));
+			mm_idx_spsc_read2(mi, fn_spsc, mm_max_spsc_bonus(&opt), spsc_scale);
 			if (mi->spsc == 0 && mm_verbose >= 2)
 				fprintf(stderr, "[WARNING] failed to load the splice score file\n");
 		}

minimap.h

@@ -5,7 +5,7 @@
 #include <stdio.h>
 #include <sys/types.h>
 
-#define MM_VERSION "2.29-r1283"
+#define MM_VERSION "2.31-r1302"
 
 #define MM_F_NO_DIAG       (0x001LL) // no exact diagonal hit
 #define MM_F_NO_DUAL       (0x002LL) // skip pairs where query name is lexicographically larger than target name
@@ -163,7 +163,8 @@ typedef struct {
 	int transition; // transition mismatch score (A:G, C:T)
 	int sc_ambi; // score when one or both bases are "N"
 	int noncan;      // cost of non-canonical splicing sites
-	int junc_bonus, junc_pen;
+	int junc_bonus;  // bonus for a splice site in annotation
+	int junc_pen;    // penalty for GT- or -AG not scored in --spsc
 	int zdrop, zdrop_inv;   // break alignment if alignment score drops too fast along the diagonal
 	int end_bonus;
 	int min_dp_max;  // drop an alignment if the score of the max scoring segment is below this threshold
@@ -407,6 +408,7 @@ int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_
  * @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_ds(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
@@ -420,6 +422,7 @@ int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uin
 
 int mm_max_spsc_bonus(const mm_mapopt_t *mo);
 int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc);
+int32_t mm_idx_spsc_read2(mm_idx_t *idx, const char *fn, int32_t max_sc, float scale);
 int64_t mm_idx_spsc_get(const mm_idx_t *db, int32_t cid, int64_t st0, int64_t en0, int32_t rev, uint8_t *sc);
 
 // deprecated APIs for backward compatibility

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "18 April 2025" "minimap2-2.29 (r1283)" "Bioinformatics tools"
+.TH minimap2 1 "19 May 2026" "minimap2-2.31 (r1302)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -443,14 +443,23 @@ line corresponds to a donor site and `A' for an acceptor site.
 A positive score suggests the junction is preferred and a negative score
 suggests the junction is not preferred.
 .TP
-.BR --junc-pen \ INT
-Penalty for a position not in FILE specified by
+.BR --spsc0 \ INT
+Penalty for positions not in
+.I FILE
+specified by
 .B --spsc
 [5]. Effective with
 .B --spsc
 but not
 .BR --junc-bed .
 .TP
+.BR --spsc-scale \ FLOAT
+Scale splice scores in
+.B --spsc
+by
+.IR FLOAT
+rounded to the nearest integer [0.7].
+.TP
 .BR --junc-bed \ FILE
 Junctions to prefer during base alignment [].
 Same format as

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = '2.29-r1283';
+var paftools_version = '2.31-r1302';
 
 /*****************************
  ***** Library functions *****
@@ -1740,10 +1740,11 @@ function paf_gff2bed(args)
 
 function paf_sam2paf(args)
 {
-	var c, pri_only = false, long_cs = false, pri_pri_only = false;
-	while ((c = getopt(args, "pPL")) != null) {
+	var c, pri_only = false, long_cs = false, pri_pri_only = false, allow_unmapped = false;
+	while ((c = getopt(args, "pPUL")) != null) {
 		if (c == 'p') pri_only = true;
 		else if (c == 'P') pri_pri_only = pri_only = true;
+		else if (c == 'U') allow_unmapped = true;
 		else if (c == 'L') long_cs = true;
 	}
 	if (args.length == getopt.ind) {
@@ -1751,6 +1752,7 @@ function paf_sam2paf(args)
 		print("Options:");
 		print("  -p      convert primary or supplementary alignments only");
 		print("  -P      convert primary alignments only");
+		print("  -U      convert unmapped reads as well");
 		print("  -L      output the cs tag in the long form");
 		exit(1);
 	}
@@ -1775,7 +1777,15 @@ function paf_sam2paf(args)
 		var flag = parseInt(t[1]);
 		if (t[9] != '*' && t[10] != '*' && t[9].length != t[10].length)
 			throw Error("at line " + lineno + ": inconsistent SEQ and QUAL lengths - " + t[9].length + " != " + t[10].length);
-		if (t[2] == '*' || (flag&4) || t[5] == '*') continue;
+        if (t[2] == '*' || (flag&4) || t[5] == '*') {
+			if (allow_unmapped) {
+            // emit an unmapped PAF line instead of skipping
+            // fields: qname, qlen, qstart, qend, strand, tname, tlen, tstart, tend, n_match, aln_len, mapq
+            var qlen_val = (t[9] == '*' ? 0 : t[9].length);
+            print([t[0], qlen_val, 0, 0, '*', '*', 0, 0, 0, 0, 0, 0].join("\t"));
+			}
+            continue;
+        }
 		if (pri_only && (flag&0x100)) continue;
 		if (pri_pri_only && (flag&0x900)) continue;
 		var tlen = ctg_len[t[2]];
@@ -2338,6 +2348,41 @@ function paf_mason2fq(args)
 	buf2.destroy();
 }
 
+// convert Mason read names to BED
+function paf_sim2bed(args)
+{
+	if (args.length == 0) {
+		print("Usage: paftools.js sim2bed <sim.txt>");
+		exit(1);
+	}
+	var buf = new Bytes();
+	var file = new File(args[0]);
+	while (file.readline(buf) >= 0) {
+		var line = buf.toString();
+		var t = line.split("!");
+		if (t.length < 5) continue;
+		var chr = t[1], st, en, strand;
+		if (t[2].indexOf("_") >= 0) { // mason paired-end
+			var pos = t[2].split("_");
+			var end = t[3].split("_");
+			var m = /^(.)(.)\/([12])$/.exec(t[4]);
+			if (m == null) continue;
+			strand = m[3] == "1" ? m[1] : m[2];
+			var read_no = parseInt(m[3]) - 1;
+			st = parseInt(pos[read_no]);
+			en = parseInt(end[read_no]);
+		} else { // badread/pbsim long reads
+			st = parseInt(t[2]);
+			en = parseInt(t[3]);
+			strand = t[4];
+		}
+		if (st > en) { var tmp = st; st = en; en = tmp; }
+		print([chr, st, en, line, 0, strand].join("\t"));
+	}
+	file.close();
+	buf.destroy();
+}
+
 // convert pbsim MAF to FASTQ
 function paf_pbsim2fq(args)
 {
@@ -2395,6 +2440,53 @@ function paf_pbsim2fq(args)
 	buf2.destroy();
 }
 
+function paf_badread2fa(args)
+{
+	if (args.length < 2) {
+		print("Usage: paftools.js badread2fa <ref.fa.fai> <badread.fq>");
+		exit(1);
+	}
+
+	var len = {}, file, buf = new Bytes();
+	file = new File(args[0]);
+	while (file.readline(buf) >= 0) {
+		var t = buf.toString().split("\t");
+		len[t[0]] = parseInt(t[1]);
+	}
+	file.close();
+
+	var id = 0, n_discard = 0;
+	file = new File(args[1]);
+	while (file.readline(buf) >= 0) {
+		var line = buf.toString();
+		var m, tag = '', a = null, is_fq = line[0] == '@'? true : false;
+		if (!/\schimera\s/.test(line) && (m = /\s(\S+),([+-])strand,(\d+)-(\d+).*read_identity=([0-9\.]+)%/.exec(line)) != null) {
+			if (len[m[1]] == null) throw Error("failed to find the contig length of " + m[1]);
+			m[3] = parseInt(m[3]);
+			m[4] = parseInt(m[4]);
+			if (m[2] == '+')
+				a = [ "S" + (id+1), m[1], m[3], m[4], m[2] ];
+			else
+				a = [ "S" + (id+1), m[1], len[m[1]] - m[4], len[m[1]] - m[3], m[2] ];
+			tag = "ri:f:" + m[5];
+		}
+		file.readline(buf);
+		var seq = buf.toString();
+		if (is_fq) {
+			file.readline(buf);
+			file.readline(buf);
+		}
+		if (a != null) {
+			print(">" + a.join("!"), tag);
+			print(seq);
+		} else ++n_discard;
+		++id;
+	}
+	file.close();
+	buf.destroy();
+	warn("WARNING: discarded " + n_discard + " reads");
+}
+
 function paf_junceval(args)
 {
 	var c, l_fuzzy = 0, print_ovlp = false, print_err_only = false, first_only = false, chr_only = false, aa = false, is_bed = false;
@@ -2613,7 +2705,8 @@ function paf_junceval(args)
 function paf_exoneval(args) // adapted from paf_junceval()
 {
 	var c, l_fuzzy = 0, print_ovlp = false, print_err_only = false, first_only = false, chr_only = false, aa = false, is_bed = false, use_cds = false, eval_base = false;
-	while ((c = getopt(args, "l:epcab1ds")) != null) {
+	var skip_start = false, skip_last = false;
+	while ((c = getopt(args, "l:epcab1dsft")) != null) {
 		if (c == 'l') l_fuzzy = parseInt(getopt.arg);
 		else if (c == 'e') print_err_only = print_ovlp = true;
 		else if (c == 'p') print_ovlp = true;
@@ -2623,6 +2716,8 @@ function paf_exoneval(args) // adapted from paf_junceval()
 		else if (c == '1') first_only = true;
 		else if (c == 'd') use_cds = true;
 		else if (c == 's') eval_base = true;
+		else if (c == 'f') skip_start = true;
+		else if (c == 't') skip_last = skip_start = true;
 	}
 
 	if (args.length - getopt.ind < 1) {
@@ -2635,6 +2730,8 @@ function paf_exoneval(args) // adapted from paf_junceval()
 		print("  -e        print erroreous overlapping exons");
 		print("  -c        only consider alignments to /^(chr)?([0-9]+|X|Y)$/");
 		print("  -1        only process the first alignment of each query");
+		print("  -f        skip the first exon in the miniprot mode");
+		print("  -t        skip the first and the last exons");
 		print("  -b        BED as input");
 		print("  -s        compute base Sn and Sp (more memory)");
 		exit(1);
@@ -2764,6 +2861,8 @@ function paf_exoneval(args) // adapted from paf_junceval()
 				for (var i = tmp_exon.length - 1; i >= 0; --i)
 					exon.push([pos + (glen - tmp_exon[i][1]), pos + (glen - tmp_exon[i][0])]);
 			}
+			if (skip_start) exon.shift();
+			if (skip_last) exon.pop();
 		} else {
 			var tmp_st = pos;
 			while ((m = re_cigar.exec(cigar)) != null) {
@@ -3686,7 +3785,9 @@ function main(args)
 		print("  mapeval    evaluate mapping accuracy using mason2/PBSIM-simulated FASTQ");
 		print("  pafcmp     compare two PAF files");
 		print("  mason2fq   convert mason2-simulated SAM to FASTQ");
+	print("  sim2bed    convert mason2-simulated read names to BED");
 		print("  pbsim2fq   convert PBSIM-simulated MAF to FASTQ");
+		print("  badread2fa convert Baderead FASTQ to FASTA");
 		print("  junceval   evaluate splice junction consistency with known annotations");
 		print("  exoneval   evaluate exon-level consistency with known annotations");
 		print("  ov-eval    evaluate read overlap sensitivity using read-to-ref mapping");
@@ -3711,7 +3812,9 @@ function main(args)
 	else if (cmd == 'pafcmp') paf_pafcmp(args);
 	else if (cmd == 'bedcov') paf_bedcov(args);
 	else if (cmd == 'mason2fq') paf_mason2fq(args);
+	else if (cmd == 'sim2bed') paf_sim2bed(args);
 	else if (cmd == 'pbsim2fq') paf_pbsim2fq(args);
+	else if (cmd == 'badread2fa') paf_badread2fa(args);
 	else if (cmd == 'junceval') paf_junceval(args);
 	else if (cmd == 'exoneval') paf_exoneval(args);
 	else if (cmd == 'ov-eval') paf_ov_eval(args);

pe.c

@@ -8,7 +8,8 @@ void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int
 	if (pri_ratio > 0.0f && *n_ > 0) {
 		int i, k, n = *n_, n_2nd = 0;
 		int max_dist = n_segs == 2? qlens[0] + qlens[1] + max_gap_ref : 0;
-		for (i = k = 0; i < n; ++i) {
+		uint8_t *keep = (uint8_t*)kmalloc(km, n);
+		for (i = 0; i < n; ++i) {
 			int to_keep = 0;
 			if (r[i].parent == i) { // primary
 				to_keep = 1;
@@ -34,9 +35,13 @@ void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int
 			if (to_keep && r[i].parent != i) {
 				if (n_2nd++ >= best_n) to_keep = 0; // don't keep if there are too many secondary hits
 			}
-			if (to_keep) r[k++] = r[i];
+			keep[i] = to_keep;
+		}
+		for (i = k = 0; i < n; ++i) {
+			if (keep[i]) r[k++] = r[i];
 			else if (r[i].p) free(r[i].p);
 		}
+		kfree(km, keep);
 		if (k != n) mm_sync_regs(km, k, r); // removing hits requires sync()
 		*n_ = k;
 	}

python/cmappy.pxd

@@ -112,6 +112,7 @@ cdef extern from "minimap.h":
 	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_ds(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)
 
 #

python/mappy.pyx

@@ -3,7 +3,7 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
-__version__ = '2.29'
+__version__ = '2.31'
 
 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, _cs, _MD # these are python objects
+	cdef _ctg, _cigar, _cs, _ds, _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, cs_str, MD_str):
+	def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand, seg_id, cs_str, ds_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
@@ -27,6 +27,7 @@ cdef class Alignment:
 		self._trans_strand = trans_strand
 		self._seg_id = seg_id
 		self._cs = cs_str
+		self._ds = ds_str
 		self._MD = MD_str
 
 	@property
@@ -77,6 +78,9 @@ cdef class Alignment:
 	@property
 	def cs(self): return self._cs
 
+	@property
+	def ds(self): return self._ds
+
 	@property
 	def MD(self): return self._MD
 
@@ -96,6 +100,7 @@ cdef class Alignment:
 		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)
+		if self._ds != "": a.append("ds:Z:" + self._ds)
 		if self._MD != "": a.append("MD:Z:" + self._MD)
 		return "\t".join(a)
 
@@ -165,7 +170,7 @@ cdef class Aligner:
 	def __bool__(self):
 		return (self._idx != NULL)
 
-	def map(self, seq, seq2=None, name=None, buf=None, cs=False, MD=False, max_frag_len=None, extra_flags=None):
+	def map(self, seq, seq2=None, name=None, buf=None, cs=False, ds=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
@@ -184,7 +189,6 @@ cdef class Aligner:
 		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 name is not None:
@@ -206,19 +210,23 @@ cdef class Aligner:
 			i = 0
 			while i < n_regs:
 				cmappy.mm_reg2hitpy(self._idx, &regs[i], &h)
-				cigar, _cs, _MD = [], '', ''
+				cigar, _cs, _ds, _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])
-				if cs or MD: # generate the cs and/or the MD tag, if requested
+				if cs or ds or MD: # generate the cs/ds and/or the MD tag, if requested
+					km = cmappy.mm_tbuf_get_km(b._b)
 					_cur_seq = _seq2 if h.seg_id > 0 and seq2 is not None else _seq
 					if cs:
 						l_cs_str = cmappy.mm_gen_cs(km, &cs_str, &m_cs_str, self._idx, &regs[i], _cur_seq, 1)
 						_cs = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
+					if ds:
+						l_cs_str = cmappy.mm_gen_ds(km, &cs_str, &m_cs_str, self._idx, &regs[i], _cur_seq, 1)
+						_ds = 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], _cur_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)
+				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, _ds, _MD)
 				cmappy.mm_free_reg1(&regs[i])
 				i += 1
 		finally:

python/minimap2.py

@@ -5,7 +5,7 @@ import getopt
 import mappy as mp
 
 def main(argv):
-	opts, args = getopt.getopt(argv[1:], "x:n:m:k:w:r:cM")
+	opts, args = getopt.getopt(argv[1:], "x:n:m:k:w:r:cdM")
 	if len(args) < 2:
 		print("Usage: minimap2.py [options] <ref.fa>|<ref.mmi> <query.fq>")
 		print("Options:")
@@ -16,11 +16,12 @@ def main(argv):
 		print("  -w INT      minimizer window length")
 		print("  -r INT      band width")
 		print("  -c          output the cs tag")
+		print("  -d          output the ds tag")
 		print("  -M          output the MD tag")
 		sys.exit(1)
 
 	preset = min_cnt = min_sc = k = w = bw = None
-	out_cs = out_MD = False
+	out_cs = out_ds = out_MD = False
 	for opt, arg in opts:
 		if opt == '-x': preset = arg
 		elif opt == '-n': min_cnt = int(arg)
@@ -29,12 +30,13 @@ def main(argv):
 		elif opt == '-k': k = int(arg)
 		elif opt == '-w': w = int(arg)
 		elif opt == '-c': out_cs = True
+		elif opt == '-d': out_ds = True
 		elif opt == '-M': out_MD = 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, cs=out_cs, MD=out_MD): # traverse hits
+		for h in a.map(seq, cs=out_cs, ds=out_ds, MD=out_MD): # traverse hits
 			print('{}\t{}\t{}'.format(name, len(seq), h))
 
 if __name__ == "__main__":

setup.py

@@ -23,7 +23,7 @@ def readme():
 
 setup(
 	name = 'mappy',
-	version = '2.29',
+	version = '2.31',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),

test/x3s-aln.txt

@@ -0,0 +1,5 @@
+mm2: TGTTATCCCTAGGGTAACTTGTTCCGTTGGTCAAGTTATTGGATCAATTGAGTATAGTAGTGCACTCAC......................................................................................................................................CACTTGGAGCCATTCATACAGGTCCCTATTTAAGGAACAAGTGATTATGCTACCTTTGCACGGTT
+     |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||                                                                                                                                      |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+ref: TGTTATCCCTAGGGTAACTTGTTCCGTTGGTCAAGTTATTGGATCAATTGAGTATAGTAGTGCACTCACctGCTTCGCTTTGACTGGTGAAGTCTTAGCATGTACTGCTCGGAGGTTGGGTTCTGCTCCGAGGTCGCCCCAACCGAAATTTTTAATGCAGGTTTGGTAGTTTAGGACCTGTGGGTTTGTTAGGCTAACCTCacCACTTGGAGCCATTCATACAGGTCCCTATTTAAGGAACAAGTGATTATGCTACCTTTGCACGGTT
+     ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||                                                                                                                                      ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+sta: TGTTATCCCTAGGGTAACTTGTTCCGTTGGTCAAGTTATTGGATCAATTGAGTATAGTAGTGCA......................................................................................................................................CTCACCACTTGGAGCCATTCATACAGGTCCCTATTTAAGGAACAAGTGATTATGCTACCTTTGCACGGTT

test/x3s-qry.fa

@@ -0,0 +1,5 @@
+>query
+AACCGTGCAAAGGTAGCATAATCACTTGTTCCTTAAATAGGGACCTGTATGAATGGCTCC
+AAGTG
+GTGAGTGCA
+CTACTATACTCAATTGATCCAATAACTTGACCAACGGAACAAGTTACCCTAGGGATAACA

test/x3s-ref.fa

@@ -0,0 +1,10 @@
+>ref
+TGATCCAACATCGAGGTCGTAAACCCTATTGTTGATATGGACTCTAGAATAGGATTGCGC
+TGTTATCCCTAGGGTAACTTGTTCCGTTGGTCAAGTTATTGGATCAATTGAGTATAGTAG
+TGCACTCAC
+ctGCTTCGCTTTGACTGGTGAAGTCTTAGCATGTACTGCTCGGAGGTTGGGTTCTGCTCC
+GAGGTCGCCCCAACCGAAATTTTTAATGCAGGTTTGGTAGTTTAGGACCTGTGGGTTTGT
+TAGGCTAACCTCac
+CACTTGGAGCCATTCATACAGGTCCCTATTTAAGGAACAAGTGATTATGCTACCTTTGCA
+CGGTTAGGGTACCGCGGCCGTTAAACATGTGTCACTGGGCAGGCGGTGCCTCTAATACTG
+GTGAT