Release minimap2-2.23 (r1111)

See #816 and #823

Makefile.simde

@@ -1,7 +1,7 @@
 CFLAGS=		-g -Wall -O2 -Wc++-compat #-Wextra
 CPPFLAGS=	-DHAVE_KALLOC -DUSE_SIMDE -DSIMDE_ENABLE_NATIVE_ALIASES
 INCLUDES=	-Ilib/simde
-OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o chain.o align.o hit.o map.o format.o pe.o esterr.o splitidx.o \
+OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o lchain.o align.o hit.o map.o format.o pe.o seed.o esterr.o splitidx.o \
 			ksw2_extz2_simde.o ksw2_extd2_simde.o ksw2_exts2_simde.o ksw2_ll_simde.o
 PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite

NEWS.md

@@ -1,3 +1,46 @@
+Release 2.23-r1111 (18 November 2021)
+-------------------------------------
+
+Notable changes:
+
+ * Bugfix: fixed missing alignments around long inversions (#806 and #816).
+   This bug affected v2.19 through v2.22.
+
+ * Improvement: avoid extremely long mapping time for pathologic reads with
+   highly repeated k-mers not in the reference (#771). Use --q-occ-frac=0
+   to disable the new heuristic.
+
+ * Change: use --cap-kalloc=1g by default.
+
+(2.23: 18 November 2021, r1111)
+
+
+
+Release 2.22-r1101 (7 August 2021)
+----------------------------------
+
+When choosing the best alignment, this release uses logarithm gap penalty and
+query-specific mismatch penalty. It improves the sensitivity to long INDELs in
+repetitive regions.
+
+Other notable changes:
+
+ * Bugfix: fixed an indirect memory leak that may waste a large amount of
+   memory given highly repetitive reference such as a 16S RNA database (#749).
+   All versions of minimap2 have this issue.
+
+ * New feature: added --cap-kalloc to reduce the peak memory. This option is
+   not enabled by default but may become the default in future releases.
+
+Known issue:
+
+ * Minimap2 may take a long time to map a read (#771). So far it is not clear
+   if this happens to v2.18 and earlier versions.
+
+(2.22: 7 August 2021, r1101)
+
+
+
 Release 2.21-r1071 (6 July 2021)
 --------------------------------
 
@@ -5,7 +48,7 @@ This release fixed a regression in short-read mapping introduced in v2.19
 (#776). It also fixed invalid comparisons of uninitialized variables, though
 these are harmless (#752). Long-read alignment should be identical to v2.20.
 
-(2.21: 6 July 2021)
+(2.21: 6 July 2021, r1071)
 
 
 

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.21/minimap2-2.21_x64-linux.tar.bz2 | tar -jxvf -
+curl -L https://github.com/lh3/minimap2/releases/download/v2.23/minimap2-2.23_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.21_x64-linux/minimap2
+./minimap2-2.23_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

@@ -237,10 +237,11 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 	r->p = p;
 }
 
-static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e, int is_eqx)
+static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e, int is_eqx, int log_gap)
 {
 	uint32_t k, l;
-	int32_t s = 0, max = 0, qshift, tshift, toff = 0, qoff = 0;
+	int32_t qshift, tshift, toff = 0, qoff = 0;
+	double s = 0.0, max = 0.0;
 	mm_extra_t *p = r->p;
 	if (p == 0) return;
 	mm_fix_cigar(r, qseq, tseq, &qshift, &tshift);
@@ -265,7 +266,8 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
 			for (l = 0; l < len; ++l)
 				if (qseq[qoff + l] > 3) ++n_ambi;
 			r->blen += len - n_ambi, p->n_ambi += n_ambi;
-			s -= q + e * len;
+			if (log_gap) s -= q + (double)e * mg_log2(1.0 + len);
+			else s -= q + e;
 			if (s < 0) s = 0;
 			qoff += len;
 		} else if (op == MM_CIGAR_DEL) {
@@ -273,14 +275,15 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
 			for (l = 0; l < len; ++l)
 				if (tseq[toff + l] > 3) ++n_ambi;
 			r->blen += len - n_ambi, p->n_ambi += n_ambi;
-			s -= q + e * len;
+			if (log_gap) s -= q + (double)e * mg_log2(1.0 + len);
+			else s -= q + e;
 			if (s < 0) s = 0;
 			toff += len;
 		} else if (op == MM_CIGAR_N_SKIP) {
 			toff += len;
 		}
 	}
-	p->dp_max = max;
+	p->dp_max = (int32_t)(max + .499);
 	assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
 	if (is_eqx) mm_update_cigar_eqx(r, qseq, tseq); // NB: it has to be called here as changes to qseq and tseq are not returned
 }
@@ -533,8 +536,13 @@ static int mm_seed_ext_score(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 	re = re + ext_len < (int32_t)mi->seq[rid].len? re + ext_len : mi->seq[rid].len;
 	qe = qe + ext_len < qlen? qe + ext_len : qlen;
 	tseq = (uint8_t*)kmalloc(km, re - rs);
-	mm_idx_getseq(mi, rid, rs, re, tseq);
+	if (opt->flag & MM_F_QSTRAND) {
-	qseq = qseq0[a->x>>63] + qs;
+		qseq = qseq0[0] + qs;
+		mm_idx_getseq2(mi, a->x>>63, rid, rs, re, tseq);
+	} else {
+		qseq = qseq0[a->x>>63] + qs;
+		mm_idx_getseq(mi, rid, rs, re, tseq);
+	}
 	qp = ksw_ll_qinit(km, 2, qe - qs, qseq, 5, mat);
 	score = ksw_ll_i16(qp, re - rs, tseq, opt->q, opt->e, &q_off, &t_off);
 	kfree(km, tseq);
@@ -690,8 +698,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	junc = (uint8_t*)kmalloc(km, re0 - rs0);
 
 	if (qs > 0 && rs > 0) { // left extension; probably the condition can be changed to "qs > qs0 && rs > rs0"
-		qseq = &qseq0[rev][qs0];
+		if (opt->flag & MM_F_QSTRAND) {
-		mm_idx_getseq(mi, rid, rs0, rs, tseq);
+			qseq = &qseq0[0][qs0];
+			mm_idx_getseq2(mi, rev, rid, rs0, rs, tseq);
+		} else {
+			qseq = &qseq0[rev][qs0];
+			mm_idx_getseq(mi, rid, rs0, rs, tseq);
+		}
 		mm_idx_bed_junc(mi, rid, rs0, rs, junc);
 		mm_seq_rev(qs - qs0, qseq);
 		mm_seq_rev(rs - rs0, tseq);
@@ -720,8 +733,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			if (a[as1+i].y & MM_SEED_LONG_JOIN)
 				bw1 = qe - qs > re - rs? qe - qs : re - rs;
 			// perform alignment
-			qseq = &qseq0[rev][qs];
+			if (opt->flag & MM_F_QSTRAND) {
-			mm_idx_getseq(mi, rid, rs, re, tseq);
+				qseq = &qseq0[0][qs];
+				mm_idx_getseq2(mi, rev, rid, rs, re, tseq);
+			} else {
+				qseq = &qseq0[rev][qs];
+				mm_idx_getseq(mi, rid, rs, re, tseq);
+			}
 			mm_idx_bed_junc(mi, rid, rs, re, junc);
 			if (is_sr) { // perform ungapped alignment
 				assert(qe - qs == re - rs);
@@ -757,7 +775,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				re1 = rs + (ez->max_t + 1);
 				qe1 = qs + (ez->max_q + 1);
 				if (cnt1 - (j + 1) >= opt->min_cnt) {
-					mm_split_reg(r, r2, as1 + j + 1 - r->as, qlen, a);
+					mm_split_reg(r, r2, as1 + j + 1 - r->as, qlen, a, !!(opt->flag&MM_F_QSTRAND));
 					if (zdrop_code == 2) r2->split_inv = 1;
 				}
 				break;
@@ -767,8 +785,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	}
 
 	if (!dropped && qe < qe0 && re < re0) { // right extension
-		qseq = &qseq0[rev][qe];
+		if (opt->flag & MM_F_QSTRAND) {
-		mm_idx_getseq(mi, rid, re, re0, tseq);
+			qseq = &qseq0[0][qe];
+			mm_idx_getseq2(mi, rev, rid, re, re0, tseq);
+		} else {
+			qseq = &qseq0[rev][qe];
+			mm_idx_getseq(mi, rid, re, re0, tseq);
+		}
 		mm_idx_bed_junc(mi, rid, re, re0, junc);
 		mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, junc, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
 		if (ez->n_cigar > 0) {
@@ -781,13 +804,19 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	assert(qe1 <= qlen);
 
 	r->rs = rs1, r->re = re1;
-	if (rev) r->qs = qlen - qe1, r->qe = qlen - qs1;
+	if (!rev || (opt->flag & MM_F_QSTRAND)) r->qs = qs1, r->qe = qe1;
-	else r->qs = qs1, r->qe = qe1;
+	else r->qs = qlen - qe1, r->qe = qlen - qs1;
 
 	assert(re1 - rs1 <= re0 - rs0);
 	if (r->p) {
-		mm_idx_getseq(mi, rid, rs1, re1, tseq);
+		if (opt->flag & MM_F_QSTRAND) {
-		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e, opt->flag & MM_F_EQX);
+			mm_idx_getseq2(mi, r->rev, rid, rs1, re1, tseq);
+			qseq = &qseq0[0][qs1];
+		} else {
+			mm_idx_getseq(mi, rid, rs1, re1, tseq);
+			qseq = &qseq0[r->rev][qs1];
+		}
+		mm_update_extra(r, qseq, tseq, mat, opt->q, opt->e, opt->flag & MM_F_EQX, !(opt->flag & MM_F_SR));
 		if (rev && r->p->trans_strand)
 			r->p->trans_strand ^= 3; // flip to the read strand
 	}
@@ -797,7 +826,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 }
 
 static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez)
-{
+{ // NB: this doesn't work with the qstrand mode
 	int tl, ql, score, ret = 0, q_off, t_off;
 	uint8_t *tseq, *qseq;
 	int8_t mat[25];
@@ -846,7 +875,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	}
 	r_inv->rs = r1->re + t_off;
 	r_inv->re = r_inv->rs + ez->max_t + 1;
-	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, opt->flag & MM_F_EQX);
+	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, opt->flag & MM_F_EQX, !(opt->flag & MM_F_SR));
 	ret = 1;
 end_align1_inv:
 	kfree(km, tseq);
@@ -863,6 +892,71 @@ static inline mm_reg1_t *mm_insert_reg(const mm_reg1_t *r, int i, int *n_regs, m
 	return regs;
 }
 
+static inline void mm_count_gaps(const mm_reg1_t *r, int32_t *n_gap_, int32_t *n_gapo_)
+{
+	uint32_t i;
+	int32_t n_gapo = 0, n_gap = 0;
+	*n_gap_ = *n_gapo_ = -1;
+	if (r->p == 0) return;
+	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)
+			++n_gapo, n_gap += len;
+	}
+	*n_gap_ = n_gap, *n_gapo_ = n_gapo;
+}
+
+double mm_event_identity(const mm_reg1_t *r)
+{
+	int32_t n_gap, n_gapo;
+	if (r->p == 0) return -1.0f;
+	mm_count_gaps(r, &n_gap, &n_gapo);
+	return (double)r->mlen / (r->blen + r->p->n_ambi - n_gap + n_gapo);
+}
+
+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;
+	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_mis = r->blen + r->p->n_ambi - r->mlen - n_gap;
+	return (int32_t)(match_sc * (r->mlen - b2 * n_mis - gap_cost) + .499);
+}
+
+void mm_update_dp_max(int qlen, int n_regs, mm_reg1_t *regs, float frac, int a, int b)
+{
+	int32_t max = -1, max2 = -1, i, max_i = -1;
+	double div, b2;
+	if (n_regs < 2) return;
+	for (i = 0; i < n_regs; ++i) {
+		mm_reg1_t *r = &regs[i];
+		if (r->p == 0) continue;
+		if (r->p->dp_max > max) max2 = max, max = r->p->dp_max, max_i = i;
+		else if (r->p->dp_max > max2) max2 = r->p->dp_max;
+	}
+	if (max_i < 0 || max < 0 || max2 < 0) return;
+	if (regs[max_i].qe - regs[max_i].qs < (double)qlen * frac) return;
+	if (max2 < (double)max * frac) return;
+	div = 1. - mm_event_identity(&regs[max_i]);
+	if (div < 0.02) div = 0.02;
+	b2 = 0.5 / div; // max value: 25
+	if (b2 * a < b) b2 = (double)a / b;
+	for (i = 0; i < n_regs; ++i) {
+		mm_reg1_t *r = &regs[i];
+		if (r->p == 0) continue;
+		r->p->dp_max = mm_recal_max_dp(r, b2, a);
+		if (r->p->dp_max < 0) r->p->dp_max = 0;
+	}
+}
+
 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)
 {
 	extern unsigned char seq_nt4_table[256];
@@ -906,7 +1000,7 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 				regs[i].p->trans_strand = opt->flag&MM_F_SPLICE_FOR? 1 : 2;
 		}
 		if (r2.cnt > 0) regs = mm_insert_reg(&r2, i, &n_regs, regs);
-		if (i > 0 && regs[i].split_inv) {
+		if (i > 0 && regs[i].split_inv && !(opt->flag & MM_F_NO_INV)) {
 			if (mm_align1_inv(km, opt, mi, qlen, qseq0, &regs[i-1], &regs[i], &r2, &ez)) {
 				regs = mm_insert_reg(&r2, i, &n_regs, regs);
 				++i; // skip the inserted INV alignment
@@ -917,6 +1011,10 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 	kfree(km, qseq0[0]);
 	kfree(km, ez.cigar);
 	mm_filter_regs(opt, qlen, n_regs_, regs);
+	if (!(opt->flag&MM_F_SR) && !opt->split_prefix && qlen >= opt->rank_min_len) {
+		mm_update_dp_max(qlen, *n_regs_, regs, opt->rank_frac, opt->a, opt->b);
+		mm_filter_regs(opt, qlen, n_regs_, regs);
+	}
 	mm_hit_sort(km, n_regs_, regs, opt->alt_drop);
 	return regs;
 }

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.21/minimap2-2.21_x64-linux.tar.bz2 | tar jxf -
+curl -L https://github.com/lh3/minimap2/releases/download/v2.23/minimap2-2.23_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.21_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
+cp minimap2-2.23_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 -
@@ -80,12 +80,12 @@ where a `U`-line gives the number of unmapped reads (for SAM input only); a
 5. Accumulative number of mappings
 
 For `paftools.js mapeval` to work, you need to encode the true read positions
-in read names in the right format. For [PBSIM][pbsim] and [mason2][mason2], we
+in read names in the right format. For [pbsim2][pbsim] and [mason2][mason2], we
 provide scripts to generate the right format. Simulated reads in this cookbook
 were created with the following command lines:
 ```sh
-# in PBSIM source code directory:
+# in the pbsim2 source code directory:
-src/pbsim ../ecoli_ref.fa --depth 1 --sample-fastq sample/sample.fastq
+src/pbsim --depth 1 --length-min 5000 --length-mean 20000 --accuracy-mean 0.95 --hmm_model data/R94.model ../ecoli_ref.fa
 paftools.js pbsim2fq ../ecoli_ref.fa.fai sd_0001.maf > ../ecoli_pbsim.fa
 
 # mason2 simulation
@@ -237,7 +237,7 @@ with `-x ava-pb` (99% vs 93% with `-x ava-ont`).
 
 
 
-[pbsim]: https://github.com/pfaucon/PBSIM-PacBio-Simulator
+[pbsim]: https://github.com/yukiteruono/pbsim2
 [mason2]: https://github.com/seqan/seqan/tree/master/apps/mason2
 [paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
 [v2.10]: https://github.com/lh3/minimap2/releases/tag/v2.10

format.c

@@ -217,7 +217,7 @@ static void write_MD_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_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)
+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, int is_qstrand)
 {
 	extern unsigned char seq_nt4_table[256];
 	int i;
@@ -227,14 +227,20 @@ static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_
 	qseq = (uint8_t*)kmalloc(km, r->qe - r->qs);
 	tseq = (uint8_t*)kmalloc(km, r->re - r->rs);
 	tmp = (char*)kmalloc(km, r->re - r->rs > r->qe - r->qs? r->re - r->rs + 1 : r->qe - r->qs + 1);
-	mm_idx_getseq(mi, r->rid, r->rs, r->re, tseq);
+	if (is_qstrand) {
-	if (!r->rev) {
+		mm_idx_getseq2(mi, r->rev, r->rid, r->rs, r->re, tseq);
 		for (i = r->qs; i < r->qe; ++i)
 			qseq[i - r->qs] = seq_nt4_table[(uint8_t)t->seq[i]];
 	} else {
-		for (i = r->qs; i < r->qe; ++i) {
+		mm_idx_getseq(mi, r->rid, r->rs, r->re, tseq);
-			uint8_t c = seq_nt4_table[(uint8_t)t->seq[i]];
+		if (!r->rev) {
-			qseq[r->qe - i - 1] = c >= 4? 4 : 3 - c;
+			for (i = r->qs; i < r->qe; ++i)
+				qseq[i - r->qs] = seq_nt4_table[(uint8_t)t->seq[i]];
+		} else {
+			for (i = r->qs; i < r->qe; ++i) {
+				uint8_t c = seq_nt4_table[(uint8_t)t->seq[i]];
+				qseq[r->qe - i - 1] = c >= 4? 4 : 3 - c;
+			}
 		}
 	}
 	if (is_MD) write_MD_core(s, tseq, qseq, r, tmp, write_tag);
@@ -242,14 +248,14 @@ static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_
 	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 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)
 {
 	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);
+	write_cs_or_MD(km, &str, mi, &t, r, no_iden, is_MD, 0, is_qstrand);
 	*max_len = str.m;
 	*buf = str.s;
 	return str.l;
@@ -257,24 +263,12 @@ int mm_gen_cs_or_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, cons
 
 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);
+	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 0, 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);
+	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 1, 0, 0);
-}
-
-double mm_event_identity(const mm_reg1_t *r)
-{
-	int32_t i, n_gapo = 0, n_gap = 0;
-	if (r->p == 0) return -1.0f;
-	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)
-			++n_gapo, n_gap += len;
-	}
-	return (double)r->mlen / (r->blen + r->p->n_ambi - n_gap + n_gapo);
 }
 
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
@@ -305,7 +299,7 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 	if (r->split) mm_sprintf_lite(s, "\tzd:i:%d", r->split);
 }
 
-void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag, int rep_len)
+void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int64_t opt_flag, int rep_len)
 {
 	s->l = 0;
 	if (r == 0) {
@@ -316,7 +310,11 @@ void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const
 	mm_sprintf_lite(s, "%s\t%d\t%d\t%d\t%c\t", t->name, t->l_seq, r->qs, r->qe, "+-"[r->rev]);
 	if (mi->seq[r->rid].name) mm_sprintf_lite(s, "%s", mi->seq[r->rid].name);
 	else mm_sprintf_lite(s, "%d", r->rid);
-	mm_sprintf_lite(s, "\t%d\t%d\t%d", mi->seq[r->rid].len, r->rs, r->re);
+	mm_sprintf_lite(s, "\t%d", mi->seq[r->rid].len);
+	if ((opt_flag & MM_F_QSTRAND) && r->rev)
+		mm_sprintf_lite(s, "\t%d\t%d", mi->seq[r->rid].len - r->re, mi->seq[r->rid].len - r->rs);
+	else
+		mm_sprintf_lite(s, "\t%d\t%d", r->rs, r->re);
 	mm_sprintf_lite(s, "\t%d\t%d", r->mlen, r->blen);
 	mm_sprintf_lite(s, "\t%d", r->mapq);
 	write_tags(s, r);
@@ -328,12 +326,12 @@ 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_MD)))
-		write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1);
+		write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1, !!(opt_flag&MM_F_QSTRAND));
 	if ((opt_flag & MM_F_COPY_COMMENT) && t->comment)
 		mm_sprintf_lite(s, "\t%s", t->comment);
 }
 
-void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag)
+void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int64_t opt_flag)
 {
 	mm_write_paf3(s, mi, t, r, km, opt_flag, -1);
 }
@@ -362,7 +360,7 @@ static inline const mm_reg1_t *get_sam_pri(int n_regs, const mm_reg1_t *regs)
 	return NULL;
 }
 
-static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, const mm_reg1_t *r, int opt_flag)
+static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, const mm_reg1_t *r, int64_t opt_flag)
 {
 	if (r->p == 0) {
 		mm_sprintf_lite(s, "*");
@@ -388,7 +386,7 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
 	}
 }
 
-void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag, int rep_len)
+void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int64_t opt_flag, int rep_len)
 {
 	const int max_bam_cigar_op = 65535;
 	int flag, n_regs = n_regss[seg_idx], cigar_in_tag = 0;
@@ -535,7 +533,7 @@ void mm_write_sam3(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, 1);
+			write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1, 0);
 		if (cigar_in_tag)
 			write_sam_cigar(s, flag, 1, t->l_seq, r, opt_flag);
 	}
@@ -547,7 +545,7 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 	s->s[s->l] = 0; // we always have room for an extra byte (see str_enlarge)
 }
 
-void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag)
+void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int64_t opt_flag)
 {
 	mm_write_sam3(s, mi, t, seg_idx, reg_idx, n_seg, n_regss, regss, km, opt_flag, -1);
 }

hit.c

@@ -20,14 +20,14 @@ static inline void mm_cal_fuzzy_len(mm_reg1_t *r, const mm128_t *a)
 	}
 }
 
-static inline void mm_reg_set_coor(mm_reg1_t *r, int32_t qlen, const mm128_t *a)
+static inline void mm_reg_set_coor(mm_reg1_t *r, int32_t qlen, const mm128_t *a, int is_qstrand)
 { // NB: r->as and r->cnt MUST BE set correctly for this function to work
 	int32_t k = r->as, q_span = (int32_t)(a[k].y>>32&0xff);
 	r->rev = a[k].x>>63;
 	r->rid = a[k].x<<1>>33;
 	r->rs = (int32_t)a[k].x + 1 > q_span? (int32_t)a[k].x + 1 - q_span : 0; // NB: target span may be shorter, so this test is necessary
 	r->re = (int32_t)a[k + r->cnt - 1].x + 1;
-	if (!r->rev) {
+	if (!r->rev || is_qstrand) {
 		r->qs = (int32_t)a[k].y + 1 - q_span;
 		r->qe = (int32_t)a[k + r->cnt - 1].y + 1;
 	} else {
@@ -49,7 +49,7 @@ static inline uint64_t hash64(uint64_t key)
 	return key;
 }
 
-mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a) // convert chains to hits
+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) // convert chains to hits
 {
 	mm128_t *z, tmp;
 	mm_reg1_t *r;
@@ -81,7 +81,7 @@ mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u,
 		ri->cnt = (int32_t)z[i].y;
 		ri->as = z[i].y >> 32;
 		ri->div = -1.0f;
-		mm_reg_set_coor(ri, qlen, a);
+		mm_reg_set_coor(ri, qlen, a, is_qstrand);
 	}
 	kfree(km, z);
 	return r;
@@ -103,7 +103,7 @@ static inline int mm_alt_score(int score, float alt_diff_frac)
 	return score > 0? score : 1;
 }
 
-void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a)
+void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a, int is_qstrand)
 {
 	if (n <= 0 || n >= r->cnt) return;
 	*r2 = *r;
@@ -115,10 +115,10 @@ void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a)
 	r2->score = (int32_t)(r->score * ((float)r2->cnt / r->cnt) + .499);
 	r2->as = r->as + n;
 	if (r->parent == r->id) r2->parent = MM_PARENT_TMP_PRI;
-	mm_reg_set_coor(r2, qlen, a);
+	mm_reg_set_coor(r2, qlen, a, is_qstrand);
 	r->cnt -= r2->cnt;
 	r->score -= r2->score;
-	mm_reg_set_coor(r, qlen, a);
+	mm_reg_set_coor(r, qlen, a, is_qstrand);
 	r->split |= 1, r2->split |= 2;
 }
 
@@ -252,7 +252,7 @@ void mm_sync_regs(void *km, int n_regs, mm_reg1_t *regs) // keep mm_reg1_t::{id,
 	mm_set_sam_pri(n_regs, regs);
 }
 
-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(void *km, float pri_ratio, int min_diff, int best_n, int check_strand, int min_strand_sc, int *n_, mm_reg1_t *r)
 {
 	if (pri_ratio > 0.0f && *n_ > 0) {
 		int i, k, n = *n_, n_2nd = 0;
@@ -264,6 +264,9 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *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);
+			} 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);
 		}
 		if (k != n) mm_sync_regs(km, k, r); // removing hits requires sync()
@@ -271,6 +274,19 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_,
 	}
 }
 
+int mm_filter_strand_retained(int n_regs, mm_reg1_t *r)
+{
+	int i, k;
+	for (i = k = 0; i < n_regs; ++i) {
+		int p = r[i].parent;
+		if (!r[i].strand_retained || r[i].div < r[p].div * 5.0f) {
+			if (k < i) r[k++] = r[i];
+			else ++k;
+		}
+	}
+	return k;
+}
+
 void mm_filter_regs(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs)
 { // NB: after this call, mm_reg1_t::parent can be -1 if its parent filtered out
 	int i, k;
@@ -358,7 +374,7 @@ mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int
 		}
 	}
 	for (s = 0; s < n_segs; ++s) {
-		regs[s] = mm_gen_regs(km, hash, qlens[s], seg[s].n_u, seg[s].u, seg[s].a);
+		regs[s] = mm_gen_regs(km, hash, qlens[s], seg[s].n_u, seg[s].u, seg[s].a, 0);
 		n_regs[s] = seg[s].n_u;
 		for (i = 0; i < n_regs[s]; ++i) {
 			regs[s][i].seg_split = 1;

index.c

@@ -161,6 +161,28 @@ int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, ui
 	return en - st;
 }
 
+int mm_idx_getseq_rev(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq)
+{
+	uint64_t i, st1, en1;
+	const mm_idx_seq_t *s;
+	if (rid >= mi->n_seq || st >= mi->seq[rid].len) return -1;
+	s = &mi->seq[rid];
+	if (en > s->len) en = s->len;
+	st1 = s->offset + (s->len - en);
+	en1 = s->offset + (s->len - st);
+	for (i = st1; i < en1; ++i) {
+		uint8_t c = mm_seq4_get(mi->S, i);
+		seq[en1 - i - 1] = c < 4? 3 - c : c;
+	}
+	return en - st;
+}
+
+int mm_idx_getseq2(const mm_idx_t *mi, int is_rev, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq)
+{
+	if (is_rev) return mm_idx_getseq_rev(mi, rid, st, en, seq);
+	else return mm_idx_getseq(mi, rid, st, en, seq);
+}
+
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f)
 {
 	int i;

lchain.c

@@ -6,16 +6,6 @@
 #include "kalloc.h"
 #include "krmq.h"
 
-static inline float mg_log2(float x) // NB: this doesn't work when x<2
-{
-	union { float f; uint32_t i; } z = { x };
-	float log_2 = ((z.i >> 23) & 255) - 128;
-	z.i &= ~(255 << 23);
-	z.i += 127 << 23;
-	log_2 += (-0.34484843f * z.f + 2.02466578f) * z.f - 0.67487759f;
-	return log_2;
-}
-
 uint64_t *mg_chain_backtrack(void *km, int64_t n, const int32_t *f, const int64_t *p, int32_t *v, int32_t *t, int32_t min_cnt, int32_t min_sc, int32_t *n_u_, int32_t *n_v_)
 {
 	mm128_t *z;

main.c

@@ -7,7 +7,7 @@
 #include "mmpriv.h"
 #include "ketopt.h"
 
-#define MM_VERSION "2.21-r1071"
+#define MM_VERSION "2.23-r1111"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -72,6 +72,10 @@ static ko_longopt_t long_options[] = {
 	{ "alt-drop",       ko_required_argument, 345 },
 	{ "mask-len",       ko_required_argument, 346 },
 	{ "rmq",            ko_optional_argument, 347 },
+	{ "qstrand",        ko_no_argument,       348 },
+	{ "cap-kalloc",     ko_required_argument, 349 },
+	{ "q-occ-frac",     ko_required_argument, 350 },
+	{ "chain-skip-scale",ko_required_argument,351 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
 	{ "version",        ko_no_argument,       'V' },
@@ -100,7 +104,7 @@ static inline int64_t mm_parse_num(const char *str)
 	return mm_parse_num2(str, 0);
 }
 
-static inline void yes_or_no(mm_mapopt_t *opt, int flag, int long_idx, const char *arg, int yes_to_set)
+static inline void yes_or_no(mm_mapopt_t *opt, int64_t flag, int long_idx, const char *arg, int yes_to_set)
 {
 	if (yes_to_set) {
 		if (strcmp(arg, "yes") == 0 || strcmp(arg, "y") == 0) opt->flag |= flag;
@@ -222,9 +226,13 @@ int main(int argc, char *argv[])
 		else if (c == 341) opt.junc_bonus = atoi(o.arg); // --junc-bonus
 		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
 		else if (c == 344) alt_list = o.arg; // --alt
 		else if (c == 345) opt.alt_drop = atof(o.arg); // --alt-drop
 		else if (c == 346) opt.mask_len = mm_parse_num(o.arg); // --mask-len
+		else if (c == 348) opt.flag |= MM_F_QSTRAND | MM_F_NO_INV; // --qstrand
+		else if (c == 349) opt.cap_kalloc = mm_parse_num(o.arg); // --cap-kalloc
+		else if (c == 350) opt.q_occ_frac = atof(o.arg); // --q-occ-frac
 		else if (c == 330) {
 			fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n");
 		} else if (c == 314) { // --frag
@@ -326,7 +334,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -N INT       retain at most INT secondary alignments [%d]\n", opt.best_n);
 		fprintf(fp_help, "  Alignment:\n");
 		fprintf(fp_help, "    -A INT       matching score [%d]\n", opt.a);
-		fprintf(fp_help, "    -B INT       mismatch penalty [%d]\n", opt.b);
+		fprintf(fp_help, "    -B INT       mismatch penalty (larger value for lower divergence) [%d]\n", opt.b);
 		fprintf(fp_help, "    -O INT[,INT] gap open penalty [%d,%d]\n", opt.q, opt.q2);
 		fprintf(fp_help, "    -E INT[,INT] gap extension penalty; a k-long gap costs min{O1+k*E1,O2+k*E2} [%d,%d]\n", opt.e, opt.e2);
 		fprintf(fp_help, "    -z INT[,INT] Z-drop score and inversion Z-drop score [%d,%d]\n", opt.zdrop, opt.zdrop_inv);
@@ -406,7 +414,10 @@ int main(int argc, char *argv[])
 		if (mm_verbose >= 3) mm_idx_stat(mi);
 		if (junc_bed) mm_idx_bed_read(mi, junc_bed, 1);
 		if (alt_list) mm_idx_alt_read(mi, alt_list);
-		if (argc - (o.ind + 1) == 0) continue; // no query files
+		if (argc - (o.ind + 1) == 0) {
+			mm_idx_destroy(mi);
+			continue; // no query files
+		}
 		ret = 0;
 		if (!(opt.flag & MM_F_FRAG_MODE)) {
 			for (i = o.ind + 1; i < argc; ++i) {

map.c

@@ -190,9 +190,13 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ,
 			if ((r[k]&1) == (q->q_pos&1)) { // forward strand
 				p->x = (r[k]&0xffffffff00000000ULL) | rpos;
 				p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
-			} else { // reverse strand
+			} else if (!(opt->flag & MM_F_QSTRAND)) { // reverse strand and not in the query-strand mode
 				p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | rpos;
 				p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
+			} else { // reverse strand; query-strand
+				int32_t len = mi->seq[r[k]>>32].len;
+				p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | (len - (rpos + 1 - q->q_span) - 1); // coordinate only accurate for non-HPC seeds
+				p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
 			}
 			p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
 			if (q->is_tandem) p->y |= MM_SEED_TANDEM;
@@ -208,7 +212,7 @@ static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_i
 {
 	if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
 		mm_set_parent(km, opt->mask_level, opt->mask_len, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
-		if (n_segs <= 1) mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
+		if (n_segs <= 1) mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, 1, opt->max_gap * 0.8, 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);
 	}
 }
@@ -219,7 +223,7 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k
 	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, opt->mask_len, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
-		mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
+		mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, 0, opt->max_gap * 0.8, n_regs, regs);
 		mm_set_sam_pri(*n_regs, regs);
 	}
 	return regs;
@@ -236,6 +240,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	mm128_v mv = {0,0,0};
 	mm_reg1_t *regs0;
 	km_stat_t kmst;
+	float chn_pen_gap, chn_pen_skip;
 
 	for (i = 0, qlen_sum = 0; i < n_segs; ++i)
 		qlen_sum += qlens[i], n_regs[i] = 0, regs[i] = 0;
@@ -248,6 +253,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	hash  = __ac_Wang_hash(hash);
 
 	collect_minimizers(b->km, opt, mi, n_segs, qlens, seqs, &mv);
+	if (opt->q_occ_frac > 0.0f) mm_seed_mz_flt(b->km, &mv, opt->mid_occ, opt->q_occ_frac);
 	if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 	else a = collect_seed_hits(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 
@@ -269,12 +275,14 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		if (max_chain_gap_ref < opt->max_gap) max_chain_gap_ref = opt->max_gap;
 	} else max_chain_gap_ref = opt->max_gap;
 
+	chn_pen_gap  = opt->chain_gap_scale * 0.01 * mi->k;
+	chn_pen_skip = opt->chain_skip_scale * 0.01 * mi->k;
 	if (opt->flag & MM_F_RMQ) {
 		a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
-						  opt->chain_gap_scale * 0.01 * mi->k, 0.0f, n_a, a, &n_regs0, &u, b->km);
+						  chn_pen_gap, chn_pen_skip, n_a, a, &n_regs0, &u, b->km);
 	} else {
 		a = mg_lchain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score,
-						 opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
+						 chn_pen_gap, chn_pen_skip, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 	}
 
 	if (opt->bw_long > opt->bw && (opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN)) == 0 && n_segs == 1 && n_regs0 > 1) { // re-chain/long-join for long sequences
@@ -285,7 +293,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 			kfree(b->km, u);
 			radix_sort_128x(a, a + n_a);
 			a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw_long, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
-							  opt->chain_gap_scale * 0.01 * mi->k, 0.0f, n_a, a, &n_regs0, &u, b->km);
+							  chn_pen_gap, chn_pen_skip, n_a, a, &n_regs0, &u, b->km);
 		}
 	} else if (opt->max_occ > opt->mid_occ && rep_len > 0 && !(opt->flag & MM_F_RMQ)) { // re-chain, mostly for short reads
 		int rechain = 0;
@@ -308,13 +316,13 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 			if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 			else a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 			a = mg_lchain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score,
-							 opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
+							 chn_pen_gap, chn_pen_skip, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 		}
 	}
 	b->frag_gap = max_chain_gap_ref;
 	b->rep_len = rep_len;
 
-	regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a);
+	regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a, !!(opt->flag&MM_F_QSTRAND));
 	if (mi->n_alt) {
 		mm_mark_alt(mi, n_regs0, regs0);
 		mm_hit_sort(b->km, &n_regs0, regs0, opt->alt_drop); // this step can be merged into mm_gen_regs(); will do if this shows up in profile
@@ -327,10 +335,14 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 						i == regs0[j].as? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));
 
 	chain_post(opt, max_chain_gap_ref, mi, b->km, qlen_sum, n_segs, qlens, &n_regs0, regs0, a);
-	if (!is_sr) mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos);
+	if (!is_sr && !(opt->flag&MM_F_QSTRAND)) {
+		mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos);
+		n_regs0 = mm_filter_strand_retained(n_regs0, regs0);
+	}
 
 	if (n_segs == 1) { // uni-segment
 		regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], &n_regs0, regs0, a);
+		regs0 = (mm_reg1_t*)realloc(regs0, sizeof(*regs0) * n_regs0);
 		mm_set_mapq(b->km, n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr);
 		n_regs[0] = n_regs0, regs[0] = regs0;
 	} else { // multi-segment
@@ -357,7 +369,9 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
 			fprintf(stderr, "QM\t%s\t%d\tcap=%ld,nCore=%ld,largest=%ld\n", qname, qlen_sum, kmst.capacity, kmst.n_cores, kmst.largest);
 		assert(kmst.n_blocks == kmst.n_cores); // otherwise, there is a memory leak
-		if (kmst.largest > 1U<<28) {
+		if (kmst.largest > 1U<<28 || (opt->cap_kalloc > 0 && kmst.capacity > opt->cap_kalloc)) {
+			if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
+				fprintf(stderr, "[W::%s] reset thread-local memory after read %s\n", __func__, qname);
 			km_destroy(b->km);
 			b->km = km_init();
 		}
@@ -402,10 +416,13 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
     step_t *s = (step_t*)_data;
 	int qlens[MM_MAX_SEG], j, off = s->seg_off[i], pe_ori = s->p->opt->pe_ori;
 	const char *qseqs[MM_MAX_SEG];
+	double t = 0.0;
 	mm_tbuf_t *b = s->buf[tid];
 	assert(s->n_seg[i] <= MM_MAX_SEG);
-	if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
+	if (mm_dbg_flag & MM_DBG_PRINT_QNAME) {
 		fprintf(stderr, "QR\t%s\t%d\t%d\n", s->seq[off].name, tid, s->seq[off].l_seq);
+		t = realtime();
+	}
 	for (j = 0; j < s->n_seg[i]; ++j) {
 		if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1))))
 			mm_revcomp_bseq(&s->seq[off + j]);
@@ -437,6 +454,8 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
 				r->rev = !r->rev;
 			}
 		}
+	if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
+		fprintf(stderr, "QT\t%s\t%d\t%.6f\n", s->seq[off].name, tid, realtime() - t);
 }
 
 static void merge_hits(step_t *s)
@@ -481,10 +500,18 @@ static void merge_hits(step_t *s)
 					}
 				}
 			}
+			if (!(opt->flag&MM_F_SR) && s->seq[k].l_seq >= opt->rank_min_len)
+				mm_update_dp_max(s->seq[k].l_seq, s->n_reg[k], s->reg[k], opt->rank_frac, opt->a, opt->b);
+			for (j = 0; j < s->n_reg[k]; ++j) {
+				mm_reg1_t *r = &s->reg[k][j];
+				if (r->p) r->p->dp_max2 = 0; // reset ->dp_max2 as mm_set_parent() doesn't clear it; necessary with mm_update_dp_max()
+				r->subsc = 0; // this may not be necessary
+				r->n_sub = 0; // n_sub will be an underestimate as we don't see all the chains now, but it can't be accurate anyway
+			}
 			mm_hit_sort(km, &s->n_reg[k], s->reg[k], opt->alt_drop);
 			mm_set_parent(km, opt->mask_level, opt->mask_len, s->n_reg[k], s->reg[k], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
 			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_select_sub(km, opt->pri_ratio, s->p->mi->k*2, opt->best_n, 0, opt->max_gap * 0.8, &s->n_reg[k], s->reg[k]);
 				mm_set_sam_pri(s->n_reg[k], s->reg[k]);
 			}
 			mm_set_mapq(km, s->n_reg[k], s->reg[k], opt->min_chain_score, opt->a, rep_len, !!(opt->flag & MM_F_SR));

minimap.h

@@ -36,7 +36,9 @@
 #define MM_F_NO_END_FLT    0x10000000
 #define MM_F_HARD_MLEVEL   0x20000000
 #define MM_F_SAM_HIT_ONLY  0x40000000
-#define MM_F_RMQ           0x80000000LL
+#define MM_F_RMQ           (0x80000000LL)
+#define MM_F_QSTRAND       (0x100000000LL)
+#define MM_F_NO_INV        (0x200000000LL)
 
 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2
@@ -106,7 +108,7 @@ typedef struct {
 	int32_t mlen, blen;     // seeded exact match length; seeded alignment block length
 	int32_t n_sub;          // number of suboptimal mappings
 	int32_t score0;         // initial chaining score (before chain merging/spliting)
-	uint32_t mapq:8, split:2, rev:1, inv:1, sam_pri:1, proper_frag:1, pe_thru:1, seg_split:1, seg_id:8, split_inv:1, is_alt:1, dummy:6;
+	uint32_t mapq:8, split:2, rev:1, inv:1, sam_pri:1, proper_frag:1, pe_thru:1, seg_split:1, seg_id:8, split_inv:1, is_alt:1, strand_retained:1, dummy:5;
 	uint32_t hash;
 	float div;
 	mm_extra_t *p;
@@ -133,6 +135,7 @@ typedef struct {
 	int min_cnt;         // min number of minimizers on each chain
 	int min_chain_score; // min chaining score
 	float chain_gap_scale;
+	float chain_skip_scale;
 	int rmq_size_cap, rmq_inner_dist;
 	int rmq_rescue_size;
 	float rmq_rescue_ratio;
@@ -155,14 +158,19 @@ typedef struct {
 	int anchor_ext_len, anchor_ext_shift;
 	float max_clip_ratio; // drop an alignment if BOTH ends are clipped above this ratio
 
+	int rank_min_len;
+	float rank_frac;
+
 	int pe_ori, pe_bonus;
 
 	float mid_occ_frac;  // only used by mm_mapopt_update(); see below
+	float q_occ_frac;
 	int32_t min_mid_occ, max_mid_occ;
 	int32_t mid_occ;     // ignore seeds with occurrences above this threshold
 	int32_t max_occ, max_max_occ, occ_dist;
 	int64_t mini_batch_size; // size of a batch of query bases to process in parallel
 	int64_t max_sw_mat;
+	int64_t cap_kalloc;
 
 	const char *split_prefix;
 } mm_mapopt_t;

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "6 July 2021" "minimap2-2.21 (r1071)" "Bioinformatics tools"
+.TH minimap2 1 "18 November 2021" "minimap2-2.23 (r1111)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -151,10 +151,16 @@ Lower and upper bounds of k-mer occurrences [10,1000000]. The final k-mer occurr
 .BR -f }}.
 This option prevents excessively small or large
 .B -f
-estimated from the input reference. It deprecates
+estimated from the input reference. Available since r1034 and deprecating
 .B --min-occ-floor
 in earlier versions of minimap2.
 .TP
+.BI --q-occ-frac \ FLOAT
+Discard a query minimizer if its occurrence is higher than
+.I FLOAT
+fraction of query minimizers and than the reference occurrence threshold
+[0.01]. Set 0 to disable. Available since r1105.
+.TP
 .BI -e \ INT
 Sample a high-frequency minimizer every
 .I INT
@@ -423,6 +429,11 @@ alignment.
 Skip alignment if the DP matrix size is above
 .IR NUM .
 Set 0 to disable [100m].
+.TP
+.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].
 .SS Input/output options
 .TP 10
 .B -a
@@ -573,7 +584,7 @@ Up to 20% sequence divergence.
 .B splice
 Long-read spliced alignment
 .RB ( -k15
-.B -w5 --splice -g2k -G200k -A1 -B2 -O2,32 -E1,0 -C9 -z200 -ub --junc-bonus=9 --cap-sw-mem=0
+.B -w5 --splice -g2k -G200k -A1 -B2 -O2,32 -E1,0 -b0 -C9 -z200 -ub --junc-bonus=9 --cap-sw-mem=0
 .BR --splice-flank=yes ).
 In the splice mode, 1) long deletions are taken as introns and represented as
 the
@@ -592,7 +603,7 @@ Long-read splice alignment for PacBio CCS reads
 .B sr
 Short single-end reads without splicing
 .RB ( -k21
-.B -w11 --sr --frag=yes -A2 -B8 -O12,32 -E2,1 -r100 -p.5 -N20 -f1000,5000 -n2 -m20
+.B -w11 --sr --frag=yes -A2 -B8 -O12,32 -E2,1 -b0 -r100 -p.5 -N20 -f1000,5000 -n2 -m20
 .B -s40 -g100 -2K50m --heap-sort=yes
 .BR --secondary=no ).
 .TP

misc/mmphase.js

@@ -0,0 +1,335 @@
+#!/usr/bin/env k8
+
+var getopt = function(args, ostr) {
+	var oli; // option letter list index
+	if (typeof(getopt.place) == 'undefined')
+		getopt.ind = 0, getopt.arg = null, getopt.place = -1;
+	if (getopt.place == -1) { // update scanning pointer
+		if (getopt.ind >= args.length || args[getopt.ind].charAt(getopt.place = 0) != '-') {
+			getopt.place = -1;
+			return null;
+		}
+		if (getopt.place + 1 < args[getopt.ind].length && args[getopt.ind].charAt(++getopt.place) == '-') { // found "--"
+			++getopt.ind;
+			getopt.place = -1;
+			return null;
+		}
+	}
+	var optopt = args[getopt.ind].charAt(getopt.place++); // character checked for validity
+	if (optopt == ':' || (oli = ostr.indexOf(optopt)) < 0) {
+		if (optopt == '-') return null; //  if the user didn't specify '-' as an option, assume it means null.
+		if (getopt.place < 0) ++getopt.ind;
+		return '?';
+	}
+	if (oli+1 >= ostr.length || ostr.charAt(++oli) != ':') { // don't need argument
+		getopt.arg = null;
+		if (getopt.place < 0 || getopt.place >= args[getopt.ind].length) ++getopt.ind, getopt.place = -1;
+	} else { // need an argument
+		if (getopt.place >= 0 && getopt.place < args[getopt.ind].length)
+			getopt.arg = args[getopt.ind].substr(getopt.place);
+		else if (args.length <= ++getopt.ind) { // no arg
+			getopt.place = -1;
+			if (ostr.length > 0 && ostr.charAt(0) == ':') return ':';
+			return '?';
+		} else getopt.arg = args[getopt.ind]; // white space
+		getopt.place = -1;
+		++getopt.ind;
+	}
+	return optopt;
+}
+
+function read_fastx(file, buf)
+{
+	if (file.readline(buf) < 0) return null;
+	var m, line = buf.toString();
+	if ((m = /^([>@])(\S+)/.exec(line)) == null)
+		throw Error("wrong fastx format");
+	var is_fq = (m[1] == '@');
+	var name = m[2];
+	if (file.readline(buf) < 0)
+		throw Error("missing sequence line");
+	var seq = buf.toString();
+	if (is_fq) { // skip quality
+		file.readline(buf);
+		file.readline(buf);
+	}
+	return [name, seq];
+}
+
+function filter_paf(a, opt)
+{
+	if (a.length == 0) return;
+	var k = 0;
+	for (var i = 0; i < a.length; ++i) {
+		var ai = a[i];
+		if (ai[10] < opt.min_blen) continue;
+		if (ai[9] < ai[10] * opt.min_iden) continue;
+		var clip = [0, 0];
+		if (ai[4] == '+') {
+			clip[0] = ai[2] < ai[7]? ai[2] : ai[7];
+			clip[1] = ai[1] - ai[3] < ai[6] - ai[8]? ai[1] - ai[3] : ai[6] - ai[8];
+		} else {
+			clip[0] = ai[2] < ai[6] - ai[8]? ai[2] : ai[6] - ai[8];
+			clip[1] = ai[1] - ai[3] < ai[7]? ai[1] - ai[3] : ai[7];
+		}
+		if (clip[0] > opt.max_clip_len || clip[1] > opt.max_clip_len) continue;
+		a[k++] = ai;
+	}
+	a.length = k;
+}
+
+function parse_events(t, ev, id, buf)
+{
+	var re = /(:(\d+))|(([\+\-\*])([a-z]+))/g;
+	var m, cs = null;
+	for (var j = 12; j < t.length; ++j) {
+		if ((m = /^cs:Z:(\S+)/.exec(t[j])) != null) {
+			cs = m[1].toLowerCase();
+			break;
+		}
+	}
+	if (cs == null) {
+		warn("Warning: no cs tag for read '" + t[0] + "'");
+		return;
+	}
+	var st = t[2], en = t[3];
+	var x = st;
+	while ((m = re.exec(cs)) != null) {
+		var l;
+		if (m[2] != null) { // an identitcal match ":\d+"
+			l = parseInt(m[2]);
+			// [start, end, type, index, changed_base]
+			ev.push([x, x + l, 0, id]);
+		} else {
+			if (m[4] == '*') {
+				l = 1;
+				ev.push([x, x + 1, 1, id, m[5][0]]);
+			} else if (m[4] == '+') {
+				l = m[5].length;
+				ev.push([x, x + l, 2, id]);
+			} else if (m[4] == '-') {
+				l = 0;
+				ev.push([x, x, -1, id, m[5]]);
+			}
+		}
+		x += l;
+	}
+	if (x != en)
+		throw Error("inconsistent cs for read '" + t[0] + "'");
+}
+
+function find_het_sub(ev, a, opt)
+{
+	var n = a.length, last0_i = -1, h = [], d = [];
+	for (var i = 0; i < n; ++i) h[i] = [], d[i] = [];
+	for (var i = 0; i < ev.length; ++i) {
+		if (ev[i][2] == 0) {
+			if (last0_i < 0 || ev[i][0] != ev[last0_i][0]) last0_i = i;
+			else if (ev[i][1] > ev[last0_i][1])
+				last0_i = i;
+		} else if (ev[i][2] == 1 && last0_i >= 0 && ev[i][0] < ev[last0_i][1]) {
+			if (ev[last0_i][1] - ev[last0_i][0] >= opt.min_mlen) {
+				if (opt.dbg_ev) print("EV", ev[last0_i].join("\t"), "|", ev[i].join("\t"));
+				var e0 = ev[last0_i], hl = h[e0[3]];
+				if (hl.length == 0 || hl[hl.length-1][0] != e0[0])
+					hl.push([e0[0], e0[1]]);
+				d[ev[i][3]].push([ev[i][0], e0[1] - e0[0]]);
+			}
+		}
+	}
+	var b = [];
+	for (var i = 0; i < n; ++i) {
+		var sh = 0, dh = 0;
+		for (var j = 0; j < h[i].length; ++j)
+			sh += h[i][j][1] - h[i][j][0];
+		for (var j = 0; j < d[i].length; ++j)
+			dh += d[i][j][1];
+		// [start, end, index, #consistent, lenConsistent, #conflictive, lenConflictive, identity, mlen]
+		b[i] = [a[i][2], a[i][3], i, h[i].length, sh, d[i].length, dh, a[i][9] / a[i][10], a[i][9]];
+	}
+	return b;
+}
+
+function flt_utg_for_ec(b, opt)
+{
+	var k = 0;
+	for (var i = 0; i < b.length; ++i) {
+		var bi = b[i];
+		if (bi[4] == 0 && bi[6] == 0) b[k++] = bi; // entirely ambiguous
+		else if (bi[6] < (bi[4] + bi[6]) * opt.max_ratio0) b[k++] = bi;
+	}
+	b.length = k;
+	if (b.length == 0) return;
+	// find the longest contiguous segment
+	b.sort(function(x,y) { return x[0]-y[0] });
+	var st = b[0][0], en = b[0][1], max_st = 0, max_en = 0, max_max_en = en;
+	for (var i = 1; i < b.length; ++i) {
+		if (b[i][0] > en) {
+			if (en - st > max_en - max_st)
+				max_st = st, max_en = en;
+			st = b[i][0], en = b[i][1];
+		} else {
+			en = en > b[i][1]? en : b[i][1];
+		}
+		max_max_en = max_max_en > b[i][1]? max_max_en : b[i][1];
+	}
+	if (en - st > max_en - max_st)
+		max_st = st, max_en = en;
+	if (max_max_en != en || st != b[0][0]) {
+		var k = 0;
+		for (var i = 0; i < b.length; ++i)
+			if (b[i][0] < max_en && b[i][1] > max_st)
+				b[k++] = b[i];
+		b.length = k;
+	}
+}
+
+function flt_utg_for_bin(b, opt) // filter out alignments clearly on the wrong phase
+{
+	var k = 0;
+	for (var i = 0; i < b.length; ++i) {
+		var bi = b[i];
+		if (bi[4] + bi[6] == 0 || bi[4] >= (bi[4] + bi[6]) * opt.max_ratio0) b[k++] = bi;
+	}
+	b.length = k;
+}
+
+function ec_core(b, n_a, ev, buf, ecb) // error correction
+{
+	var intv = [];
+	for (var i = 0; i < n_a; ++i)
+		intv[i] = null;
+	intv[b[0][2]] = [b[0][0], b[0][1]];
+	var en = b[0][1];
+	for (var i = 1; i < b.length; ++i) {
+		if (b[i][1] <= en) continue;
+		intv[b[i][2]] = [en, b[i][1]];
+		en = b[i][1];
+	}
+	var k = 0;
+	ecb.capacity = buf.capacity;
+	ecb.length = 0;
+	for (var i = 0; i < ev.length; ++i) {
+		var e = ev[i], I = intv[e[3]];
+		if (I == null) continue;
+		if (e[0] >= I[0] && e[0] < I[1]) { // this is to reduce duplicated events around junctions
+			//print("X", e.join("\t"));
+			if (e[2] == 0) {
+				ecb.length += e[1] - e[0];
+				for (var j = e[0]; j < e[1]; ++j)
+					ecb[k++] = buf[j];
+			} else if (e[2] == 1) {
+				++ecb.length;
+				ecb[k++] = e[4].charCodeAt(0);
+			} else if (e[2] < 0) {
+				ecb.length += e[4].length;
+				for (var j = 0; j < e[4].length; ++j)
+					ecb[k++] = e[4].charCodeAt(j);
+			} // else, skip e[2] == 2
+		}
+	}
+	if (ecb.length != k) throw Error("BUG!");
+}
+
+function process_paf(a, opt, fp_seq, buf, ecb)
+{
+	if (a.length == 0) return;
+	var len = a[0][1], name = a[0][0], seq = null;
+	if (len < opt.min_rlen) return;
+	if (fp_seq) {
+		var ret;
+		while ((ret = read_fastx(fp_seq, buf)) != null)
+			if (ret[0] == a[0][0])
+				break;
+		if (ret == null)
+			throw Error("failed to find sequence for read '" + a[0][0] + "'");
+		name = ret[0], seq = ret[1];
+		if (seq.length != len)
+			throw Error("inconsistent length for read '" + name + "'");
+	}
+	filter_paf(a, opt);
+	if (a.length == 0) return;
+	var ev = [];
+	for (var i = 0; i < a.length; ++i)
+		parse_events(a[i], ev, i, buf);
+	ev.sort(function(x,y) { return x[0]!=y[0]? x[0]-y[0] : x[2]-y[2] });
+	if (seq == null) print("SQ", name, a[0][1], a.length);
+	var b = find_het_sub(ev, a, opt);
+	if (opt.ec) flt_utg_for_ec(b, opt);
+	else flt_utg_for_bin(b, opt);
+	if (seq == null) {
+		for (var i = 0; i < b.length; ++i) {
+			var m, ai = a[b[i][2]], score = 0;
+			for (var j = 10; j < ai.length; ++j)
+				if ((m = /^AS:i:(\d+)/.exec(ai[j])) != null)
+					score = m[1];
+			print("TS", b[i][2], b[i][0], b[i][1], ai.slice(5, 9).join("\t"), b[i].slice(3, 7).join("\t"), score);
+		}
+		print("//");
+	} else { // error correction
+		if (b.length == 0) return;
+		buf.set(seq, 0);
+		ec_core(b, a.length, ev, buf, ecb);
+		print(">" + name);
+		print(ecb);
+	}
+}
+
+function main(args)
+{
+	var c, opt = { min_rlen:5000, min_blen:5000, min_iden:0.8, min_mlen:5, max_clip_len:500, max_ratio0:0.25, dbg_ev:false };
+	while ((c = getopt(args, "l:b:d:m:c:r:E")) != null) {
+		if (c == 'l') opt.min_rlen = parseInt(getopt.arg);
+		else if (c == 'b') opt.min_blen = parseInt(getopt.arg);
+		else if (c == 'd') opt.min_iden = parseFloat(getopt.arg);
+		else if (c == 'm') opt.min_slen = parseInt(getopt.arg);
+		else if (c == 'c') opt.max_clip_len = parseInt(getopt.arg);
+		else if (c == 'r') opt.max_ratio0 = parseFloat(getopt.arg);
+		else if (c == 'E') opt.dbg_ev = true;
+	}
+	if (args.length - getopt.ind < 1) {
+		print("Usage: mmphase.js [options] <map-with-cs.paf> [reads.fa]");
+		print("Options:");
+		print("  -l INT    min read length [" + opt.min_rlen + "]");
+		print("  -b INT    min alignment length [" + opt.min_blen + "]");
+		print("  -d FLOAT  min identity [" + opt.min_iden + "]");
+		print("  -s INT    min match length [" + opt.min_mlen + "]");
+		print("  -c INT    max clip length [" + opt.max_clip_len + "]");
+		print("  -r FLOAT  initial ratio for haplotype filtering [" + opt.max_ratio0 + "]");
+		return 0;
+	}
+
+	opt.ec = args.length - getopt.ind < 2? false : true;
+	if (!opt.ec) {
+		print("CC");
+		print("CC", "SQ  qName  qLen    nHits");
+		print("CC", "TS  index  qStart  qEnd  tName  tLen  tStart  tEnd  nConsistent  lCons  nConflictive  lConf  score");
+		print("CC");
+	}
+
+	var buf = new Bytes(), ecb = new Bytes();
+	var fp_paf = new File(args[getopt.ind]);
+	var fp_seq = args.length - getopt.ind >= 2? new File(args[getopt.ind+1]) : null;
+	var a = [];
+	while (fp_paf.readline(buf) >= 0) {
+		var t = buf.toString().split("\t");
+		if (a.length > 0 && a[0][0] != t[0]) {
+			process_paf(a, opt, fp_seq, buf, ecb);
+			a.length = 0;
+		}
+		for (var i = 1; i <= 3; ++i) t[i] = parseInt(t[i]);
+		if (t[1] < opt.min_rlen) continue;
+		for (var i = 6; i <= 10; ++i) t[i] = parseInt(t[i]);
+		if (t[10] < opt.min_blen) continue;
+		a.push(t);
+	}
+	if (a.length >= 0)
+		process_paf(a, opt, fp_seq, buf, ecb);
+	if (fp_seq) fp_seq.close();
+	fp_paf.close();
+	ecb.destroy();
+	buf.destroy();
+}
+
+var ret = main(arguments)
+exit(ret)

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = '2.21-r1071';
+var paftools_version = '2.23-r1111';
 
 /*****************************
  ***** Library functions *****
@@ -977,7 +977,7 @@ function paf_stat(args)
 	var re = /(\d+)([MIDSHNX=])/g;
 
 	var lineno = 0, n_pri = 0, n_2nd = 0, n_seq = 0, n_cigar_64k = 0, l_tot = 0, l_cov = 0;
-	var n_gap = [[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]];
+	var n_gap = [[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]], n_sub = 0;
 
 	function cov_len(regs)
 	{
@@ -999,7 +999,7 @@ function paf_stat(args)
 		if (line.charAt(0) != '@') {
 			var t = line.split("\t", 12);
 			var m, rs, cigar = null, is_pri = false, is_sam = false, is_rev = false, tname = null;
-			var atlen = null, aqlen, qs, qe, mapq, ori_qlen, NM = null;
+			var atlen = null, aqlen, qs, qe, mapq, ori_qlen, NM = null, nn = 0;
 			if (t.length < 2) continue;
 			if (t[4] == '+' || t[4] == '-' || t[4] == '*') { // PAF
 				if (t[4] == '*') continue; // unmapped
@@ -1009,6 +1009,8 @@ function paf_stat(args)
 				}
 				if ((m = /\tNM:i:(\d+)/.exec(line)) != null)
 					NM = parseInt(m[1]);
+				if ((m = /\tnn:i:(\d+)/.exec(line)) != null)
+					nn = parseInt(m[1]);
 				if ((m = /\tcg:Z:(\S+)/.exec(line)) != null)
 					cigar = m[1];
 				if (cigar == null) {
@@ -1032,6 +1034,8 @@ function paf_stat(args)
 				}
 				if ((m = /\tNM:i:(\d+)/.exec(line)) != null)
 					NM = parseInt(m[1]);
+				if ((m = /\tnn:i:(\d+)/.exec(line)) != null)
+					nn = parseInt(m[1]);
 				cigar = t[5];
 				tname = t[2];
 				rs = parseInt(t[3]) - 1;
@@ -1078,6 +1082,12 @@ function paf_stat(args)
 					clip[M == 0? 0 : 1] = l;
 				}
 			}
+			if (NM != null) {
+				var tmp = NM - n_gap_all - nn;
+				if (tmp < 0 && nn == 0) warn("WARNING: NM is smaller than the number of gaps at line " + lineno + ": NM=" + NM + ", nn=" + nn + ", G=" + n_gap_all);
+				if (tmp < 0) tmp = 0;
+				n_sub += tmp;
+			}
 			if (n_cigar > 65535) ++n_cigar_64k;
 			if (ql + sclip != aqlen)
 				warn("WARNING: aligned query length is inconsistent with CIGAR at line " + lineno + " (" + (ql+sclip) + " != " + aqlen + ")");
@@ -1112,6 +1122,7 @@ function paf_stat(args)
 		print("Number of primary alignments with >65535 CIGAR operations: " + n_cigar_64k);
 		print("Number of bases in mapped sequences: " + l_tot);
 		print("Number of mapped bases: " + l_cov);
+		print("Number of substitutions: " + n_sub);
 		print("Number of insertions in [0,50): " + n_gap[0][0]);
 		print("Number of insertions in [50,100): " + n_gap[0][1]);
 		print("Number of insertions in [100,300): " + n_gap[0][2]);
@@ -1475,8 +1486,14 @@ function paf_view(args)
 				warn("WARNING: converting to BLAST-like alignment requires the 'cs' tag, which is absent on line " + lineno);
 				continue;
 			}
+			var n_mm = 0, n_oi = 0, n_od = 0, n_ei = 0, n_ed = 0;
+			while ((m = re_cs.exec(cs)) != null) {
+				if (m[1] == '*') ++n_mm;
+				else if (m[1] == '+') ++n_oi, n_ei += m[2].length;
+				else if (m[1] == '-') ++n_od, n_ed += m[2].length;
+			}
 			line = line.replace(/\tc[sg]:Z:\S+/g, ""); // get rid of cs or cg tags
-			print('>' + line);
+			print('>' + line + "\tmm:i:"+n_mm + "\toi:i:"+n_oi + "\tei:i:"+n_ei + "\tod:i:"+n_od + "\ted:i:"+n_ed);
 			var rs = parseInt(t[7]), qs = t[4] == '+'? parseInt(t[2]) : parseInt(t[3]);
 			var n_blocks = 0;
 			while ((m = re_cs.exec(cs)) != null) {
@@ -1515,12 +1532,13 @@ function paf_view(args)
 
 function paf_gff2bed(args)
 {
-	var c, fn_ucsc_fai = null, is_short = false, keep_gff = false, print_junc = false;
+	var c, fn_ucsc_fai = null, is_short = false, keep_gff = false, print_junc = false, output_gene = false;
-	while ((c = getopt(args, "u:sgj")) != null) {
+	while ((c = getopt(args, "u:sgjG")) != null) {
 		if (c == 'u') fn_ucsc_fai = getopt.arg;
 		else if (c == 's') is_short = true;
 		else if (c == 'g') keep_gff = true;
 		else if (c == 'j') print_junc = true;
+		else if (c == 'G') output_gene = true;
 	}
 
 	if (getopt.ind == args.length) {
@@ -1588,8 +1606,10 @@ function paf_gff2bed(args)
 		print(a[0][0], st, en, name, 1000, a[0][3], cds_st, cds_en, color, a.length, sizes.join(",") + ",", starts.join(",") + ",");
 	}
 
-	var re_gtf = /\b(transcript_id|transcript_type|transcript_biotype|gene_name|gene_id|gbkey|transcript_name) "([^"]+)";/g;
+	var re_gtf  = /\b(transcript_id|transcript_type|transcript_biotype|gene_name|gene_id|gbkey|transcript_name) "([^"]+)";/g;
 	var re_gff3 = /\b(transcript_id|transcript_type|transcript_biotype|gene_name|gene_id|gbkey|transcript_name)=([^;]+)/g;
+	var re_gtf_gene  = /\b(gene_id|gene_type|gene_name) "([^;]+)";/g;
+	var re_gff3_gene = /\b(gene_id|gene_type|source_gene|gene_biotype|gene_name)=([^;]+);/g;
 	var buf = new Bytes();
 	var file = args[getopt.ind] == '-'? new File() : new File(args[getopt.ind]);
 
@@ -1603,6 +1623,26 @@ function paf_gff2bed(args)
 			continue;
 		}
 		if (t[0].charAt(0) == '#') continue;
+		if (output_gene) {
+			var id = null, src = null, biotype = null, type = "", name = "N/A";
+			if (t[2] != "gene") continue;
+			while ((m = re_gtf_gene.exec(t[8])) != null) {
+				if (m[1] == "gene_id") id = m[2];
+				else if (m[1] == "gene_type") type = m[2];
+				else if (m[1] == "gene_name") name = m[2];
+			}
+			while ((m = re_gff3_gene.exec(t[8])) != null) {
+				if (m[1] == "gene_id") id = m[2];
+				else if (m[1] == "source_gene") src = m[2];
+				else if (m[1] == "gene_type") type = m[2];
+				else if (m[1] == "gene_biotype") biotype = m[2];
+				else if (m[1] == "gene_name") name = m[2];
+			}
+			if (src != null) id = src;
+			if (type == "" && biotype != null) type = biotype;
+			print(t[0], parseInt(t[3]) - 1, t[4], [id, type, name].join("|"), 1000, t[6]);
+			continue;
+		}
 		if (t[2] != "CDS" && t[2] != "exon") continue;
 		t[3] = parseInt(t[3]) - 1;
 		t[4] = parseInt(t[4]);
@@ -2930,6 +2970,120 @@ function paf_vcfsel(args)
 	buf.destroy();
 }
 
+function paf_pafcmp(args)
+{
+	var c, opt = { min_len:5000, min_mapq:10, min_ovlp:0.5 };
+	while ((c = getopt(args, "q:")) != null) {
+		if (c == 'q') opt.min_mapq = parseInt(getopt.arg);
+	}
+
+	var buf = new Bytes();
+	if (args.length - getopt.ind < 2) {
+		print("Usage: paftools.js pafcmp [options] <base.paf> <test.paf>");
+		print("Options:");
+		print("  -q INT    min mapping quality [" + opt.min_mapq + "]");
+		return 1;
+	}
+
+	var eval = { n_base:0, n_test:0, n_out_high:0, n_out_low:0, n_hit:0, n_wrong:0, n_miss:0 };
+
+	function process_base(base, a) {
+		if (a.length != 1) return;
+		for (var i = 1; i < 4; ++i)
+			a[0][i] = parseInt(a[0][i]);
+		for (var i = 6; i < 12; ++i)
+			a[0][i] = parseInt(a[0][i]);
+		if (a[0][1] < opt.min_len) return;
+		if (a[0][11] >= opt.min_mapq) ++eval.n_base;
+		base[a[0][0]] = [a[0][5], a[0][7], a[0][8], a[0][11], 0, 0];
+	}
+
+	var file = new File(args[getopt.ind]);
+	warn("Reading " + args[getopt.ind] + "...");
+	var a = [], base = {};
+	while (file.readline(buf) >= 0) {
+		var line = buf.toString();
+		var t = line.split("\t");
+		if (/\ttp:A:S/.test(line)) continue;
+		if (a.length > 0 && a[0][0] != t[0]) {
+			process_base(base, a);
+			a = [];
+		}
+		a.push(t);
+	}
+	process_base(base, a);
+	file.close();
+
+	function process_test(base, a) {
+		for (var i = 1; i < 4; ++i)
+			a[0][i] = parseInt(a[0][i]);
+		for (var i = 6; i < 12; ++i)
+			a[0][i] = parseInt(a[0][i]);
+		if (a[0][1] < opt.min_len) return;
+		if (a[0][11] >= opt.min_mapq) ++eval.n_test;
+		var c = [a[0][5], a[0][7], a[0][8], a[0][11]];
+		if (base[a[0][0]] == null) {
+			if (c[3] >= opt.min_mapq) ++opt.n_out_high;
+			else ++opt.n_out_low;
+		} else {
+			var b = base[a[0][0]];
+			var inter = 0, union = (b[2] - b[1]) + (c[2] - c[1]);
+			if (b[0] == c[0]) { // same chr
+				if (b[1] < c[1]) {
+					if (b[2] > c[1])
+						inter = b[2] - c[1], union = c[2] - b[1];
+				} else { // c[1] < b[1]
+					if (c[2] > b[1])
+						inter = c[2] - b[1], union = b[2] - c[1];
+				}
+			}
+			if (inter >= union * opt.min_ovlp) {
+				if (b[3] >= opt.min_mapq) ++eval.n_hit;
+				++b[4];
+			} else {
+				if (b[3] >= opt.min_mapq) {
+					print("W", a[0][0], b.slice(0, 4).join("\t"), c.join("\t"));
+					++eval.n_wrong;
+				}
+				++b[5];
+			}
+		}
+	}
+
+	file = new File(args[getopt.ind+1]);
+	warn("Reading " + args[getopt.ind+1] + "...");
+	a = [];
+	while (file.readline(buf) >= 0) {
+		var line = buf.toString();
+		var t = line.split("\t");
+		if (/\ttp:A:S/.test(line)) continue;
+		if (a.length > 0 && a[0][0] != t[0]) {
+			process_test(base, a);
+			a = [];
+		}
+		a.push(t);
+	}
+	process_test(base, a);
+	file.close();
+
+	for (var r in base) {
+		var b = base[r];
+		if (b[3] >= opt.min_mapq && b[4] == 0 && b[5] == 0) {
+			++eval.n_miss;
+			print("M", r, b.slice(0, 4).join("\t"));
+		}
+	}
+
+	print("X", eval.n_base     + " base alignments with mapQ>=" + opt.min_mapq);
+//	print("X", eval.n_test     + " test alignments with mapQ>=" + opt.min_mapq);
+	print("X", eval.n_hit      + " base alignments correctly mapped by test");
+	print("X", eval.n_wrong    + " wrong test alignment");
+	print("X", eval.n_miss     + " base alignments missing");
+	print("X", eval.n_out_high + " additional test alignments with mapQ>=" + opt.min_mapq);
+
+	buf.destroy();
+}
+
 /*************************
  ***** main function *****
  *************************/
@@ -2957,6 +3111,7 @@ function main(args)
 		print("  version    print paftools.js version");
 		print("");
 		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("  pbsim2fq   convert PBSIM-simulated MAF to FASTQ");
 		print("  junceval   evaluate splice junction consistency with known annotations");
@@ -2978,6 +3133,7 @@ function main(args)
 	else if (cmd == 'vcfpair') paf_vcfpair(args);
 	else if (cmd == 'call') paf_call(args);
 	else if (cmd == 'mapeval') paf_mapeval(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 == 'pbsim2fq') paf_pbsim2fq(args);

mmpriv.h

@@ -61,18 +61,22 @@ uint32_t ks_ksmall_uint32_t(size_t n, uint32_t arr[], size_t kk);
 void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, int is_hpc, mm128_v *p);
 
 mm_seed_t *mm_collect_matches(void *km, int *_n_m, int qlen, int max_occ, int max_max_occ, int dist, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos);
+void mm_seed_mz_flt(void *km, mm128_v *mv, int32_t q_occ_max, float q_occ_frac);
 
+double mm_event_identity(const mm_reg1_t *r);
 int mm_write_sam_hdr(const mm_idx_t *mi, const char *rg, const char *ver, int argc, char *argv[]);
-void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag);
+void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int64_t opt_flag);
-void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag, int rep_len);
+void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int64_t opt_flag, int rep_len);
 void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int n_regs, const mm_reg1_t *regs);
-void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regs, const mm_reg1_t *const* regs, void *km, int opt_flag);
+void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regs, const mm_reg1_t *const* regs, void *km, int64_t opt_flag);
-void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag, int rep_len);
+void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int64_t opt_flag, int rep_len);
 
 void mm_idxopt_init(mm_idxopt_t *opt);
 const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f);
+int mm_idx_getseq2(const mm_idx_t *mi, int is_rev, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
 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);
@@ -81,18 +85,19 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
 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,
 					   int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
 
-mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a);
 void mm_mark_alt(const mm_idx_t *mi, int n, mm_reg1_t *r);
-void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a);
+void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a, int is_qstrand);
 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 mask_len, int n, mm_reg1_t *r, int sub_diff, int hard_mask_level, float alt_diff_frac);
-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(void *km, float pri_ratio, int min_diff, int best_n, int check_strand, int min_strand_sc, 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);
+int mm_filter_strand_retained(int n_regs, mm_reg1_t *r);
 void mm_filter_regs(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs);
 void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r, float alt_diff_frac);
 void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr);
+void mm_update_dp_max(int qlen, int n_regs, mm_reg1_t *regs, float frac, int a, int b);
 
 void mm_est_err(const mm_idx_t *mi, int qlen, int n_regs, mm_reg1_t *regs, const mm128_t *a, int32_t n, const uint64_t *mini_pos);
 
@@ -109,6 +114,16 @@ void mm_err_puts(const char *str);
 void mm_err_fwrite(const void *p, size_t size, size_t nitems, FILE *fp);
 void mm_err_fread(void *p, size_t size, size_t nitems, FILE *fp);
 
+static inline float mg_log2(float x) // NB: this doesn't work when x<2
+{
+	union { float f; uint32_t i; } z = { x };
+	float log_2 = ((z.i >> 23) & 255) - 128;
+	z.i &= ~(255 << 23);
+	z.i += 127 << 23;
+	log_2 += (-0.34484843f * z.f + 2.02466578f) * z.f - 0.67487759f;
+	return log_2;
+}
+
 #ifdef __cplusplus
 }
 #endif

options.c

@@ -19,6 +19,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->min_mid_occ = 10;
 	opt->max_mid_occ = 1000000;
 	opt->sdust_thres = 0; // no SDUST masking
+	opt->q_occ_frac = 0.01f;
 
 	opt->min_cnt = 3;
 	opt->min_chain_score = 40;
@@ -32,6 +33,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->rmq_rescue_size = 1000;
 	opt->rmq_rescue_ratio = 0.1f;
 	opt->chain_gap_scale = 0.8f;
+	opt->chain_skip_scale = 0.0f;
 	opt->max_max_occ = 4095;
 	opt->occ_dist = 500;
 
@@ -52,6 +54,10 @@ 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->rank_min_len = 500;
+	opt->rank_frac = 0.9f;
 
 	opt->pe_ori = 0; // FF
 	opt->pe_bonus = 33;
@@ -218,5 +224,10 @@ int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
 			fprintf(stderr, "[ERROR]\033[1;31m -X/-P and --secondary=no can't be applied at the same time\033[0m\n");
 		return -5;
 	}
+	if ((mo->flag & MM_F_QSTRAND) && ((mo->flag & (MM_F_OUT_SAM|MM_F_SPLICE|MM_F_FRAG_MODE)) || (io->flag & MM_I_HPC))) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m --qstrand doesn't work with -a, -H, --frag or --splice\033[0m\n");
+		return -5;
+	}
 	return 0;
 }

python/cmappy.pxd

@@ -23,6 +23,7 @@ cdef extern from "minimap.h":
 		int min_cnt
 		int min_chain_score
 		float chain_gap_scale
+		float chain_skip_scale
 		int rmq_size_cap, rmq_inner_dist
 		int rmq_rescue_size
 		float rmq_rescue_ratio
@@ -45,14 +46,19 @@ cdef extern from "minimap.h":
 		int anchor_ext_len, anchor_ext_shift
 		float max_clip_ratio
 
+		int rank_min_len
+		float rank_frac
+
 		int pe_ori, pe_bonus
 
 		float mid_occ_frac
+		float q_occ_frac
 		int32_t min_mid_occ
 		int32_t mid_occ
 		int32_t max_occ
 		int64_t mini_batch_size
 		int64_t max_sw_mat
+		int64_t cap_kalloc
 
 		const char *split_prefix
 

python/mappy.pyx

@@ -3,7 +3,7 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
-__version__ = '2.21'
+__version__ = '2.23'
 
 cmappy.mm_reset_timer()
 

seed.c

@@ -2,6 +2,31 @@
 #include "kalloc.h"
 #include "ksort.h"
 
+void mm_seed_mz_flt(void *km, mm128_v *mv, int32_t q_occ_max, float q_occ_frac)
+{
+	mm128_t *a;
+	size_t i, j, st;
+	if (mv->n <= q_occ_max || q_occ_frac <= 0.0f || q_occ_max <= 0) return;
+	KMALLOC(km, a, mv->n);
+	for (i = 0; i < mv->n; ++i)
+		a[i].x = mv->a[i].x, a[i].y = i;
+	radix_sort_128x(a, a + mv->n);
+	for (st = 0, i = 1; i <= mv->n; ++i) {
+		if (i == mv->n || a[i].x != a[st].x) {
+			int32_t cnt = i - st;
+			if (cnt > q_occ_max && cnt > mv->n * q_occ_frac)
+				for (j = st; j < i; ++j)
+					mv->a[a[j].y].x = 0;
+			st = i;
+		}
+	}
+	kfree(km, a);
+	for (i = j = 0; i < mv->n; ++i)
+		if (mv->a[i].x != 0)
+			mv->a[j++] = mv->a[i];
+	mv->n = j;
+}
+
 mm_seed_t *mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv, int32_t *n_m_)
 {
 	mm_seed_t *m;

setup.py

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

tex/minimap2.bib

@@ -338,3 +338,123 @@
 	Title = {Introducing difference recurrence relations for faster semi-global alignment of long sequences},
 	Volume = {19},
 	Year = {2018}}
+
+@article{Li:2018ab,
+	Author = {Li, Heng},
+	Journal = {Bioinformatics},
+	Pages = {3094-3100},
+	Title = {Minimap2: pairwise alignment for nucleotide sequences},
+	Volume = {34},
+	Year = {2018}}
+
+@article{Jain:2020aa,
+	Author = {Jain, Chirag and others},
+	Journal = {Bioinformatics},
+	Pages = {i111-i118},
+	Title = {Weighted minimizer sampling improves long read mapping},
+	Volume = {36},
+	Year = {2020}}
+
+@article{Miga:2020aa,
+	Author = {Miga, Karen H and others},
+	Journal = {Nature},
+	Pages = {79-84},
+	Title = {Telomere-to-telomere assembly of a complete human {X} chromosome},
+	Volume = {585},
+	Year = {2020}}
+
+@article {Jain2020.11.01.363887,
+	author = {Jain, Chirag and others},
+	title = {A long read mapping method for highly repetitive reference sequences},
+	elocation-id = {2020.11.01.363887},
+	year = {2020},
+	doi = {10.1101/2020.11.01.363887},
+	publisher = {Cold Spring Harbor Laboratory},
+	URL = {https://www.biorxiv.org/content/early/2020/11/02/2020.11.01.363887},
+	eprint = {https://www.biorxiv.org/content/early/2020/11/02/2020.11.01.363887.full.pdf},
+	journal = {bioRxiv}
+}
+
+@article{Li:2020aa,
+	Author = {Li, Heng and others},
+	Journal = {Genome Biol},
+	Pages = {265},
+	Title = {The design and construction of reference pangenome graphs with minigraph},
+	Volume = {21},
+	Year = {2020}}
+
+@article{Ren:2021aa,
+	Author = {Ren, Jingwen and Chaisson, Mark J P},
+	Journal = {PLoS Comput Biol},
+	Pages = {e1009078},
+	Title = {lra: A long read aligner for sequences and contigs},
+	Volume = {17},
+	Year = {2021}}
+
+@inproceedings{DBLP:conf/wabi/AbouelhodaO03,
+	Author = {Mohamed Ibrahim Abouelhoda and Enno Ohlebusch},
+	Booktitle = {Algorithms in Bioinformatics, Third International Workshop, {WABI} 2003, Budapest, Hungary, September 15-20, 2003, Proceedings},
+	Crossref = {DBLP:conf/wabi/2003},
+	Pages = {1--16},
+	Title = {A Local Chaining Algorithm and Its Applications in Comparative Genomics},
+	Year = {2003}}
+
+@article{Ono:2021aa,
+	Author = {Ono, Yukiteru and others},
+	Journal = {Bioinformatics},
+	Pages = {589-595},
+	Title = {{PBSIM2}: a simulator for long-read sequencers with a novel generative model of quality scores},
+	Volume = {37},
+	Year = {2021}}
+
+@article{Sedlazeck:2018ab,
+	Author = {Sedlazeck, Fritz J and others},
+	Journal = {Nat Methods},
+	Pages = {461-468},
+	Title = {Accurate detection of complex structural variations using single-molecule sequencing},
+	Volume = {15},
+	Year = {2018}}
+
+@article{Jeffares:2017aa,
+	Author = {Jeffares, Daniel C and others},
+	Journal = {Nat Commun},
+	Pages = {14061},
+	Title = {Transient structural variations have strong effects on quantitative traits and reproductive isolation in fission yeast},
+	Volume = {8},
+	Year = {2017}}
+
+@article{Zook:2020aa,
+	Author = {Zook, Justin M and others},
+	Journal = {Nat Biotechnol},
+	Pages = {1347-1355},
+	Title = {A robust benchmark for detection of germline large deletions and insertions},
+	Volume = {38},
+	Year = {2020}}
+
+@article{Harpak:2017aa,
+	Author = {Harpak, Arbel and others},
+	Journal = {Proc Natl Acad Sci U S A},
+	Pages = {12779-12784},
+	Title = {Frequent nonallelic gene conversion on the human lineage and its effect on the divergence of gene duplicates},
+	Volume = {114},
+	Year = {2017}}
+
+@article{Li:2018aa,
+	Author = {Li, Heng and others},
+	Journal = {Nat Methods},
+	Month = {Aug},
+	Number = {8},
+	Pages = {595-597},
+	Title = {A synthetic-diploid benchmark for accurate variant-calling evaluation},
+	Volume = {15},
+	Year = {2018}}
+
+@article{Gu:1995wt,
+	author = {Gu, X and Li, W H},
+	journal = {J Mol Evol},
+	month = {Apr},
+	number = {4},
+	pages = {464-73},
+	title = {The size distribution of insertions and deletions in human and rodent pseudogenes suggests the logarithmic gap penalty for sequence alignment},
+	volume = {40},
+	year = {1995}}

tex/mm2-update.tex

@@ -0,0 +1,240 @@
+\documentclass{bioinfo}
+\copyrightyear{2021}
+\pubyear{2021}
+
+\usepackage{graphicx}
+\usepackage{hyperref}
+\usepackage{url}
+\usepackage{amsmath}
+\usepackage[ruled,vlined]{algorithm2e}
+\newcommand\mycommfont[1]{\footnotesize\rmfamily{\it #1}}
+\SetCommentSty{mycommfont}
+\SetKwComment{Comment}{$\triangleright$\ }{}
+
+\usepackage{natbib}
+\bibliographystyle{apalike}
+
+\DeclareMathOperator*{\argmax}{argmax}
+
+\begin{document}
+\firstpage{1}
+
+\title[Improvements to minimap2]{New strategies to improve minimap2 alignment accuracy}
+\author[Li]{Heng Li$^{1,2}$}
+\address{$^1$Dana-Farber Cancer Institute, 450 Brookline Ave, Boston, MA 02215, USA,
+$^2$Harvard Medical School, 10 Shattuck St, Boston, MA 02215, USA}
+
+\maketitle
+
+\begin{abstract}
+
+\section{Summary:} We present several recent improvements to minimap2, a
+versatile pairwise aligner for nucleotide sequences. Now minimap2 v2.22 can
+more accurately map long reads to highly repetitive regions and align through
+insertions or deletions up to 100kb by default, addressing major weakness in
+minimap2 v2.18 or earlier.
+
+\section{Availability and implementation:}
+\href{https://github.com/lh3/minimap2}{https://github.com/lh3/minimap2}
+
+\section{Contact:} hli@ds.dfci.harvard.edu
+\end{abstract}
+
+\section{Introduction}
+Minimap2~\citep{Li:2018ab} is widely used for maping long sequence
+reads and assembly contigs. \citet{Jain:2020aa} found minimap2 v2.18 or earlier occasionally
+misaligned reads from highly repetitive regions as minimap2 ignored seeds of
+high occurrence. They also noticed minimap2 may misplace reads with structural
+variations (SVs) in such regions~\citep{Jain2020.11.01.363887}. These
+misalignments have become a pressing issue in the advent of
+temolere-to-telomore human assembly~\citep{Miga:2020aa}. Meanwhile, old minimap2
+was unable to efficiently align long insertions/deletions (INDELs) and often
+breaks an alignment around variable-number tandem repeats (VNTRs). This has
+inspired new chaining algorithms~\citep{Li:2020aa,Ren:2021aa} which are not
+integrated into minimap2. Here we will describe recent efforts implemented
+in v2.19 through v2.22 to improve mapping results.
+
+\begin{methods}
+\section{Methods}
+
+\subsection{Rescuing high-occurrence $k$-mers}\label{sec:high-occ}
+Minimap2 keeps all $k$-mer minimizers~\citep{Roberts:2004fv} during indexing. Its original
+implementation only selected low-occurrence minimizers during mapping. The
+cutoff is a few hundred for mapping long reads against a human genome. If a
+read habors only a few or even no low-occurrence minimizers, it will fail
+chaining due to insufficient anchors.
+
+To resolve this issue, we implemented a new heuristic to add additional
+minimizers. Suppose we are looking at two adjacent low-occurence $k$-mers
+located at position $x_1$ and $x_2$, respectively. If $|x_1-x_2|\ge L$,
+minimap2 v2.22 additionally selects $\lfloor|x_1-x_2|/L\rfloor$ minimizers
+of the lowest occurrence among minimizers between $x_1$ and $x_2$. Here
+parameter $L$ controls the frequency of sampling. It defaults to 500.
+This strategy adds necessary anchors at the cost of increasing total alignment
+time by a few percent on real data.
+
+\subsection{Aligning through longer INDELs}
+The original minimap2 may fail to align long INDELs due to its chaining
+heuristics. Briefly, minimap2 applies dynamic programming (DP) to chain
+minimizer anchors. This is a quadratic algorithm, slow for chaining
+contigs. For acceptable performance, the original minimap2 uses a 500bp band by
+default, which means a gap longer than 500bp will stop chaining.
+To align through longer gaps, older minimap2 implemented a long-join heurstic as follows.
+If there is an INDEL longer than 500bp and the two chains around the INDEL
+have no overlaps on either the query or the reference sequence, minimap2 may
+join the two short chains later.
+This heuristic may fail around VNTRs because short chains
+often have overlaps in VNTRs. More subtly, minimap2 may escape the inner DP
+loop early, again for performance, if the chaining result is not improved for
+50 iterations. When there is a copy number change in a long segmental
+duplication, the early escape may break around the event even if users
+specify a large band.
+
+In minigraph~\citep{Li:2020aa}, we developed a new chaining algorithm that
+finds up to 1kb INDELs with DP-based chaining and goes through longer INDELs with a
+subquadratic algorithm~\citep{DBLP:conf/wabi/AbouelhodaO03}. We ported the same
+algorithm to minimap2 for contig mapping. For long-read mapping, the minigraph
+algorithm is slower. Minimap2 v2.22 still uses the DP-based algorithm to
+find short chains and then invokes the minigraph algorithm to rechain anchors in
+these short chains. The rechaining step achieves the same goal as long-join
+but is more reliable because it can resolve overlaps between short chains. The old
+long-join heuristic has since been removed.
+
+\subsection{Properly mapping long reads with SVs}
+The original minimap2 ranks an alignment by its Smith-Waterman score and
+outputs the best scoring alignment. However, when there are SVs on the read,
+the best scoring alignment is sometimes not the correct alignment.
+\citet{Jain2020.11.01.363887} resolved this dilemma by altering the mapping
+algorithm.
+
+In our view, this problem is rooted in inapropriate scoring: affine-gap penalty
+over-penalizes a long INDEL that was often evolutionarily created in one event.
+We should not penalize a SV by a function linear in the SV length. Minimap2 v2.22 instead rescores
+an alignment with the following scoring function. Suppose an alignment consists
+of $M$ matching bases, $N$ substitutions and $G$ gap opens, we empirically
+score the alignment with
+$$
+S=M-\frac{N+G}{2d}-\sum_{i=1}^G\log_2(1+g_i)
+$$
+where $g_i\ge1$ is the length of the $i$-th gap and
+$$
+d=\max\left\{\frac{N+G}{M+N+G},0.02\right\}
+$$
+It approximates per-base sequence divergence except with the smallest value set
+to 2\%. As an analogy to affine-gap scoring, the matching score in our scheme
+is 1, the mismatch and gap open penalties are both $1/2d$ and the gap extension
+penalty is a logarithm function of the gap length~\citep{Gu:1995wt}. Our scoring gives a long SV
+a much milder penalty. In terms of time complexity, scoring an alignment is
+linear in the length of the alignment. The time spent on rescoring is negligible in
+practice.
+
+%If we assume sequences evolve under a duplication-mutation model, we may have a
+%better way to choose the best alignment. If a long read can be mapped to $n$
+%loci, we can take the read as the template and build a
+%pseudo-multi-sequence-alignment (pMSA) of $n+1$ sequences. In this pMSA, we say
+%a site on the read is informative if the $n$ reference subsequences differ at
+%the position.
+
+\end{methods}
+
+\section{Results}
+
+\begin{table}
+\processtable{Evaluation of minimap2 v2.22}
+{\footnotesize\label{tab:1}\begin{tabular}{p{4.2cm}rrrr}
+\toprule
+$[$Benchmark$]$ Metric & v2.22 & v2.18 & Winno & lra \\
+\midrule
+$[$sim-map$]$ \% mapped reads at Q10      & 97.9       & 97.6      & {\bf 99.0}& 97.3 \\
+$[$sim-map$]$ err. rate at Q10 (phredQ)   & {\bf 52}   & {\bf 52}  & 38        & 24 \\
+$[$winno-cmp$]$ rate of diff. (phredQ)    & {\bf 41}   & 37        & truth     & 18 \\
+$[$winno-cmp$]$ CPU time (hour)           & {\bf 5.0}  & 5.3       & 71.8      & 13.1 \\
+$[$winno-cmp$]$ peak RAM (Gb)             & 17.1       & 14.4      & {\bf 9.6} & 12.4 \\
+$[$sim-sv$]$  \% false negative rate      & {\bf 0.5}  & 2.0       & {\bf 0.5} & 1.4  \\
+$[$sim-sv$]$  \% false discovery rate     & {\bf 0.0}  & 0.1       & {\bf 0.0} & 0.1  \\
+$[$real-sv-1k$]$ \% false negative rate   & {\bf 7.3}  & 20.0      & 13.0      & N/A \\
+$[$real-sv-1k$]$ \% false discovery rate  & 2.7        & {\bf 2.4} & 2.7       & N/A \\
+\botrule
+\end{tabular}}
+{In $[$sim-map$]$, 152,713 reads were simulated from the CHM13 telomere-to-telomere assembly v1.1
+(AC: GCA\_009914755.3) with pbsim2~\citep{Ono:2021aa}: ``pbsim2 -{}-hmm\_model R94.model -{}-length-min
+5000 -{}-length-mean 20000 -{}-accuracy-mean 0.95''. Alignments of mapping quality
+10 or higher were evaluated by ``paftools.js mapeval''. The mapping error rate
+is measured in the phred scale: if the error rate is $e$, $-10\log_{10}e$ is
+reported in the table. In $[$winno-cmp$]$, 1.39 million CHM13 HiFi reads from
+SRR11292121 were mapped against the same CHM13 assembly. 99.3\% of them were mapped by Winnowmap2
+at mapping quality 10 or higher and were taken as ground truth to evaluate
+minimap2 and lra with ``paftools.js pafcmp''. $[$sim-sv$]$ simulated 1,000
+50bp to 1000bp INDELs from chr8 in CHM13 using SURVIVOR~\citep{Jeffares:2017aa} and simulated Nanopore
+reads at 30-fold coverage with the same pbsim2 command line. SVs were called with
+``sniffles -q 10''~\citep{Sedlazeck:2018ab} and compared to the simulated truth with ``SURVIVOR eval
+call.vcf truth.bed 50''. In $[$real-sv-1k$]$, small and long variants were
+called by dipcall-0.3~\citep{Li:2018aa} for HG002 assemblies (AC: GCA\_018852605.1 and
+GCA\_018852615.1) and compared to the GIAB truth~\citep{Zook:2020aa} using ``truvari -r 2000 -s
+1000 -S 400 -{}-multimatch -{}-passonly'' which sets the minimum INDEL size to 1kb in evaluation. }
+\end{table}
+
+We evaluated minimap2 v2.22 along with v2.18, Winnowmap2 v2.03 and lra v1.3.2
+(Table~\ref{tab:1}), using the default setting of each mapper according to the input data types.
+Both versions of minimap2 achieved high mapping accuracy on
+simulated Nanopore reads (sim-map). Winnowmap2 aligned more reads at mapping
+quality 10 or higher (mapQ10). However, it may occasionally assign a high mapping
+quality to a read with multiple identical best alignments. This reduced its
+mapping accuracy.
+
+In lack of groud truth for real data, we took Winnowmap2 mapping as ground
+truth to evaluate other mappers (winno-cmp in Table~\ref{tab:1}). Out of 1,378,092 reads with mapQ10
+alignments by Winnowmap2, minimap2 v2.22 could map all of them. 118 reads, less
+than 0.01\% of all reads, were mapped differently by v2.22. 51 of them have
+multiple identical best alignments. We believe these are more likely to be
+Winnowmap2 errors. Most of the remaining 67 (=118-51) reads have multiple
+highly similar but not identical alignments.
+Minimap2 v2.18 is less consistent with 275 differences including 30 unmapped
+reads mappable by both Winnowmap2 and v2.22.
+
+For the minimizer rescuing parameter $L$ in Section~\ref{sec:high-occ},
+we set its default to 500 such that v2.22 has comparable performance to v2.18 given simulated PacBio and Nanopore human reads.
+To see the effect of this parameter on real data, we tried several different $L$ values.
+v2.22 gave 99 mapping differences at $L=200$,
+118 at $L=500$ (default), 167 at $L=750$ and 224 differences at $L=1000$ in comparison to Winnowmap2.
+$L=200$ is 28\% slower than the default while $L=1000$ is 9\% faster.
+Changing the default minimizer window size (option ``-w'')
+and the initial minimizer occurrence cutoff (option ``-f'')
+also affects performance and accuracy to a similar magnitude.
+
+The two benchmarks above only evaluate read mappings when there are no variations between the reads and the reference.
+To measure the mapping accuracy in the presence of SVs (sim-sv), we reproduced
+the results by~\citep{Jain2020.11.01.363887}. Minimap2 v2.22 is as good as
+Winnowmap2 now. Note that we were setting the Sniffles mapping quality
+threshold to 10 in consistent with the benchmarks above. If we used the
+default threshold 20, v2.22 would miss additional five SVs (accounting for
+0.5\% of simulated SVs). For four out of these five missing SVs, minimap2 v2.22
+mapped more variant reads than Winnowmap2. Sniffles did not call these SVs
+because minimap2 tended to give them conservative mapping quality. It is worth
+noting that the simulation here only considers a simple scenario in evolution.
+Non-allelic gene conversions, which happen often in segmental
+duplications~\citep{Harpak:2017aa}, would obscure the optimal mapping
+strategies. How much such simple SV simulation informs real-world SV calling
+remains a question.
+
+To see if minimap2 v2.22 could improve long INDEL alignment, we ran dipcall on
+contig-to-reference alignments and focused on INDELs longer than 1kb
+(real-sv-1k). v2.22 is more sensitive at comparable specificity, confirming its
+advantage in more contiguous alignment. We could not get dipcall to work well with lra,
+so did not report the numbers.
+
+Minimap2 spends most computing time on base alignment. As recent improvements
+in v2.22 incur little additional computing and do not change the base alignment
+algorithm, the new version has similar performance to older versions. It is
+consistently faster than Winnowmap2 by several times. Sometimes simple
+heuristics can be as effective as more sophisticated yet slower solutions.
+
+\section*{Acknowledgements}
+We thank Arang Rhie and Chirag Jain for providing motivating examples for which
+older minimap2 underperforms.
+
+\paragraph{Funding\textcolon} This work is funded by NHGRI grant R01HG010040.
+
+\bibliography{minimap2}
+
+\end{document}