Release minimap2-2.28 (r1209)

Such that minimap2 can give the same alignment in other modes

NEWS.md

@@ -1,3 +1,32 @@
+Release 2.28-r1209 (27 March 2024)
+----------------------------------
+
+Notable changes to minimap2:
+
+ * Bugfix: `--MD` was not working properly due to the addition of `--ds` in the
+   last release (#1181 and #1182).
+
+ * New feature: added an experimental preset `lq:hqae` for aligning accurate
+   long reads back to their assembly. It has been observed that `map-hifi` and
+   `lr:hq` may produce many wrong alignments around centromeres when accurate
+   long reads (PacBio HiFi or Nanopore duplex/Q20+) are mapped to a diploid
+   assembly constructed from them. This new preset produces much more accurate
+   alignment. It is still experimental and may be subjective to changes in
+   future.
+
+ * Change: reduced the default `--cap-kalloc` to 500m to lower the peak
+   memory consumption (#855).
+
+Notable changes to mappy:
+
+ * Bugfix: mappy option struct was out of sync with minimap2 (#1177).
+
+Minimap2 should output identical alignments to v2.27.
+
+(2.28: 27 March 2024, r1209)
+
+
+
 Release 2.27-r1193 (12 March 2024)
 ----------------------------------
 

README.md

@@ -74,8 +74,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.27/minimap2-2.27_x64-linux.tar.bz2 | tar -jxvf -
+curl -L https://github.com/lh3/minimap2/releases/download/v2.28/minimap2-2.28_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.27_x64-linux/minimap2
+./minimap2-2.28_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

@@ -933,14 +933,14 @@ double mm_event_identity(const mm_reg1_t *r)
 static int32_t mm_recal_max_dp(const mm_reg1_t *r, double b2, int32_t match_sc)
 {
 	uint32_t i;
-	int32_t n_gap = 0, n_gapo = 0, n_mis;
+	int32_t n_gap = 0, n_mis;
 	double gap_cost = 0.0;
 	if (r->p == 0) return -1;
 	for (i = 0; i < r->p->n_cigar; ++i) {
 		int32_t op = r->p->cigar[i] & 0xf, len = r->p->cigar[i] >> 4;
 		if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL) {
 			gap_cost += b2 + (double)mg_log2(1.0 + len);
-			++n_gapo, n_gap += len;
+			n_gap += len;
 		}
 	}
 	n_mis = r->blen + r->p->n_ambi - r->mlen - n_gap;

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.27/minimap2-2.27_x64-linux.tar.bz2 | tar jxf -
+curl -L https://github.com/lh3/minimap2/releases/download/v2.28/minimap2-2.28_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.27_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
+cp minimap2-2.28_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

@@ -310,7 +310,7 @@ static void write_cs_ds_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const
 			}
 		}
 	}
-	if (is_MD == 1) write_MD_core(s, tseq, qseq, r, tmp, write_tag);
+	if (is_MD) write_MD_core(s, tseq, qseq, r, tmp, write_tag);
 	else write_cs_ds_core(s, tseq, qseq, r, tmp, no_iden, is_ds, write_tag);
 	kfree(km, qseq); kfree(km, tseq); kfree(km, tmp);
 }
@@ -393,7 +393,7 @@ void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const
 			mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, MM_CIGAR_STR[r->p->cigar[k]&0xf]);
 	}
 	if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_DS|MM_F_OUT_MD)))
-		write_cs_ds_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, !!(opt_flag&MM_F_OUT_DS), 1, !!(opt_flag&MM_F_QSTRAND));
+		write_cs_ds_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), !!(opt_flag&MM_F_OUT_MD), !!(opt_flag&MM_F_OUT_DS), 1, !!(opt_flag&MM_F_QSTRAND));
 	if ((opt_flag & MM_F_COPY_COMMENT) && t->comment)
 		mm_sprintf_lite(s, "\t%s", t->comment);
 }

lchain.c

@@ -149,7 +149,7 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
 					  int is_cdna, int n_seg, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km)
 { // TODO: make sure this works when n has more than 32 bits
 	int32_t *f, *t, *v, n_u, n_v, mmax_f = 0, max_drop = bw;
-	int64_t *p, i, j, max_ii, st = 0, n_iter = 0;
+	int64_t *p, i, j, max_ii, st = 0;
 	uint64_t *u;
 
 	if (_u) *_u = 0, *n_u_ = 0;
@@ -174,7 +174,6 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
 		for (j = i - 1; j >= st; --j) {
 			int32_t sc;
 			sc = comput_sc(&a[i], &a[j], max_dist_x, max_dist_y, bw, chn_pen_gap, chn_pen_skip, is_cdna, n_seg);
-			++n_iter;
 			if (sc == INT32_MIN) continue;
 			sc += f[j];
 			if (sc > max_f) {
@@ -204,6 +203,7 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
 		if (max_ii < 0 || (a[i].x - a[max_ii].x <= (int64_t)max_dist_x && f[max_ii] < f[i]))
 			max_ii = i;
 		if (mmax_f < max_f) mmax_f = max_f;
+		//fprintf(stderr, "X1\t%ld\t%ld:%d\t%ld\t%ld:%d\t%ld\t%ld\n", (long)i, (long)(a[i].x>>32), (int32_t)a[i].x, (long)max_j, max_j<0?-1L:(long)(a[max_j].x>>32), max_j<0?-1:(int32_t)a[max_j].x, (long)max_f, (long)v[i]);
 	}
 
 	u = mg_chain_backtrack(km, n, f, p, v, t, min_cnt, min_sc, max_drop, &n_u, &n_v);
@@ -263,7 +263,8 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 		return 0;
 	}
 	if (max_dist < bw) max_dist = bw;
-	if (max_dist_inner <= 0 || max_dist_inner >= max_dist) max_dist_inner = 0;
+	if (max_dist_inner < 0) max_dist_inner = 0;
+	if (max_dist_inner > max_dist) max_dist_inner = max_dist;
 	p = Kmalloc(km, int64_t, n);
 	f = Kmalloc(km, int32_t, n);
 	t = Kcalloc(km, int32_t, n);
@@ -325,12 +326,11 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 				krmq_interval(lc_elem, root_inner, &s, &lo, &hi);
 				if (lo) {
 					const lc_elem_t *q;
-					int32_t width, n_rmq_iter = 0;
+					int32_t width;
 					krmq_itr_t(lc_elem) itr;
 					krmq_itr_find(lc_elem, root_inner, lo, &itr);
 					while ((q = krmq_at(&itr)) != 0) {
 						if (q->y < (int32_t)a[i].y - max_dist_inner) break;
-						++n_rmq_iter;
 						j = q->i;
 						sc = f[j] + comput_sc_simple(&a[i], &a[j], chn_pen_gap, chn_pen_skip, 0, &width);
 						if (width <= bw) {

main.c

@@ -78,6 +78,8 @@ static ko_longopt_t long_options[] = {
 	{ "no-hash-name",   ko_no_argument,       353 },
 	{ "secondary-seq",  ko_no_argument,       354 },
 	{ "ds",             ko_no_argument,       355 },
+	{ "rmq-inner",      ko_required_argument, 356 },
+	{ "dbg-seed-occ",   ko_no_argument,       501 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
 	{ "version",        ko_no_argument,       'V' },
@@ -245,6 +247,8 @@ int main(int argc, char *argv[])
 		else if (c == 353) opt.flag |= MM_F_NO_HASH_NAME; // --no-hash-name
 		else if (c == 354) opt.flag |= MM_F_SECONDARY_SEQ; // --secondary-seq
 		else if (c == 355) opt.flag |= MM_F_OUT_DS; // --ds
+		else if (c == 356) opt.rmq_inner_dist = mm_parse_num(o.arg); // --rmq-inner
+		else if (c == 501) mm_dbg_flag |= MM_DBG_SEED_FREQ; // --dbg-seed-occ
 		else if (c == 330) {
 			fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n");
 		} else if (c == 314) { // --frag

minimap.h

@@ -5,7 +5,7 @@
 #include <stdio.h>
 #include <sys/types.h>
 
-#define MM_VERSION "2.27-r1193"
+#define MM_VERSION "2.28-r1209"
 
 #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

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "12 March 2024" "minimap2-2.27 (r1193)" "Bioinformatics tools"
+.TH minimap2 1 "12 March 2024" "minimap2-2.28 (r1209)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -268,6 +268,11 @@ or more of the shorter chain [0.5]
 Use the minigraph chaining algorithm [no]. The minigraph algorithm is better
 for aligning contigs through long INDELs.
 .TP
+.BI --rmq-inner \ NUM
+Apply full dynamic programming for anchors within distance
+.I NUM
+[1000].
+.TP
 .B --hard-mask-level
 Honor option
 .B -M
@@ -463,7 +468,7 @@ Set 0 to disable [100m].
 .BI --cap-kalloc \ NUM
 Free thread-local kalloc memory reservoir if after the alignment the size of the reservoir above
 .IR NUM .
-Set 0 to disable [0].
+Set 0 to disable [500m].
 .SS Input/output options
 .TP 10
 .B -a

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = '2.27-r1193';
+var paftools_version = '2.28-r1209';
 
 /*****************************
  ***** Library functions *****
@@ -1740,15 +1740,17 @@ function paf_gff2bed(args)
 
 function paf_sam2paf(args)
 {
-	var c, pri_only = false, long_cs = false;
+	var c, pri_only = false, long_cs = false, pri_pri_only = false;
-	while ((c = getopt(args, "pL")) != null) {
+	while ((c = getopt(args, "pPL")) != null) {
 		if (c == 'p') pri_only = true;
+		else if (c == 'P') pri_pri_only = pri_only = true;
 		else if (c == 'L') long_cs = true;
 	}
 	if (args.length == getopt.ind) {
 		print("Usage: paftools.js sam2paf [options] <in.sam>");
 		print("Options:");
 		print("  -p      convert primary or supplementary alignments only");
+		print("  -P      convert primary alignments only");
 		print("  -L      output the cs tag in the long form");
 		exit(1);
 	}
@@ -1775,6 +1777,7 @@ function paf_sam2paf(args)
 			throw Error("at line " + lineno + ": inconsistent SEQ and QUAL lengths - " + t[9].length + " != " + t[10].length);
 		if (t[2] == '*' || (flag&4) || t[5] == '*') continue;
 		if (pri_only && (flag&0x100)) continue;
+		if (pri_pri_only && (flag&0x900)) continue;
 		var tlen = ctg_len[t[2]];
 		if (tlen == null) throw Error("at line " + lineno + ": can't find the length of contig " + t[2]);
 		// find tags
@@ -1887,7 +1890,10 @@ function paf_sam2paf(args)
 		// optional tags
 		var type = flag&0x100? 'S' : 'P';
 		var tags = ["tp:A:" + type];
-		if (NM != null) tags.push("mm:i:"+mm);
+		if (NM != null) {
+			tags.push("NM:i:"+NM);
+			tags.push("mm:i:"+mm);
+		}
 		tags.push("gn:i:"+(I[1]+D[1]), "go:i:"+(I[0]+D[0]), "cg:Z:" + t[5].replace(/\d+[SH]/g, ''));
 		if (cs_str != null) tags.push("cs:Z:" + cs_str);
 		else if (cs.length > 0) tags.push("cs:Z:" + cs.join(""));

mmpriv.h

@@ -14,6 +14,7 @@
 #define MM_DBG_PRINT_SEED    0x4
 #define MM_DBG_PRINT_ALN_SEQ 0x8
 #define MM_DBG_PRINT_CHAIN   0x10
+#define MM_DBG_SEED_FREQ     0x20
 
 #define MM_SEED_LONG_JOIN  (1ULL<<40)
 #define MM_SEED_IGNORE     (1ULL<<41)
@@ -79,8 +80,6 @@ int mm_idx_getseq2(const mm_idx_t *mi, int is_rev, uint32_t rid, uint32_t st, ui
 mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
 mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a, int is_qstrand);
 
-mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float gap_scale,
-					 int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
 mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip,
 					  int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
 mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_skip, int cap_rmq_size, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip,

options.c

@@ -55,7 +55,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->max_clip_ratio = 1.0f;
 	opt->mini_batch_size = 500000000;
 	opt->max_sw_mat = 100000000;
-	opt->cap_kalloc = 1000000000;
+	opt->cap_kalloc = 500000000;
 
 	opt->rank_min_len = 500;
 	opt->rank_frac = 0.9f;
@@ -114,6 +114,14 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 			mo->a = 1, mo->b = 4, mo->q = 6, mo->q2 = 26, mo->e = 2, mo->e2 = 1;
 			mo->min_dp_max = 200;
 		}
+	} else if (strcmp(preset, "lr:hqae") == 0) { // high-quality assembly evaluation
+		io->flag = 0, io->k = 25, io->w = 51;
+		mo->flag |= MM_F_RMQ;
+		mo->min_mid_occ = 50, mo->max_mid_occ = 500;
+		mo->rmq_inner_dist = 5000;
+		mo->occ_dist = 200;
+		mo->best_n = 100;
+		mo->chain_gap_scale = 5.0f;
 	} else if (strcmp(preset, "map-iclr-prerender") == 0) {
 		io->flag = 0, io->k = 15;
 		mo->b = 6, mo->transition = 1;

python/cmappy.pxd

@@ -36,6 +36,7 @@ cdef extern from "minimap.h":
 		float alt_drop
 
 		int a, b, q, e, q2, e2
+		int transition
 		int sc_ambi
 		int noncan
 		int junc_bonus

python/mappy.pyx

@@ -3,7 +3,7 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
-__version__ = '2.27'
+__version__ = '2.28'
 
 cmappy.mm_reset_timer()
 
@@ -96,6 +96,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._MD != "": a.append("MD:Z:" + self._MD)
 		return "\t".join(a)
 
 cdef class ThreadBuffer:

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

seed.c

@@ -112,7 +112,8 @@ mm_seed_t *mm_collect_matches(void *km, int *_n_m, int qlen, int max_occ, int ma
 	}
 	for (i = 0, n_m = 0, *rep_len = 0, *n_a = 0; i < n_m0; ++i) {
 		mm_seed_t *q = &m[i];
-		//fprintf(stderr, "X\t%d\t%d\t%d\n", q->q_pos>>1, q->n, q->flt);
+		if (mm_dbg_flag & MM_DBG_SEED_FREQ)
+			fprintf(stderr, "SF\t%d\t%d\t%d\n", q->q_pos>>1, q->n, q->flt);
 		if (q->flt) {
 			int en = (q->q_pos >> 1) + 1, st = en - q->q_span;
 			if (st > rep_en) {

setup.py

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