Release minimap2-2.28 (r1209)

Such that minimap2 can give the same alignment in other modes

.travis.yml

@@ -1,24 +0,0 @@
-matrix:
-  include:
-    - language: c
-      compiler: gcc
-      script: make
-    - language: c
-      compiler: clang
-      script: make
-    - arch: arm64
-      language: c
-      compiler: gcc
-      script: make arm_neon=1 aarch64=1
-    - language: python
-      python: "2.7"
-      before_install: pip install cython
-      script: python setup.py build_ext
-    - language: python
-      python: "3.5"
-      before_install: pip install cython
-      script: python setup.py build_ext
-    - language: python
-      python: "3.9"
-      before_install: pip install cython
-      script: python setup.py build_ext

Makefile

@@ -8,6 +8,10 @@ PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite
 LIBS=		-lm -lz -lpthread
 
+ifneq ($(aarch64),)
+	arm_neon=1
+endif
+
 ifeq ($(arm_neon),) # if arm_neon is not defined
 ifeq ($(sse2only),) # if sse2only is not defined
 	OBJS+=ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o
@@ -26,12 +30,12 @@ endif
 
 ifneq ($(asan),)
 	CFLAGS+=-fsanitize=address
-	LIBS+=-fsanitize=address
+	LIBS+=-fsanitize=address -ldl
 endif
 
 ifneq ($(tsan),)
 	CFLAGS+=-fsanitize=thread
-	LIBS+=-fsanitize=thread
+	LIBS+=-fsanitize=thread -ldl
 endif
 
 .PHONY:all extra clean depend

NEWS.md

@@ -1,3 +1,117 @@
+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)
+----------------------------------
+
+Notable changes to minimap2:
+
+ * New feature: added the `lr:hq` preset for accurate long reads at ~1% error
+   rate. This was suggested by Oxford Nanopore developers (#1127). It is not
+   clear if this preset also works well for PacBio HiFi reads.
+
+ * New feature: added the `map-iclr` preset for Illumina Complete Long Reads
+   (#1069), provided by Illumina developers.
+
+ * New feature: added option `-b` to specify mismatch penalty for base
+   transitions (i.e. A-to-G or C-to-T changes).
+
+ * New feature: added option `--ds` to generate a new `ds:Z` tag that
+   indicates uncertainty in INDEL positions. It is an extension to `cs`. The
+   `mgutils-es6.js` script in minigraph parses `ds`.
+
+ * Bugfix: avoided a NULL pointer dereference (#1154). This would not have an
+   effect on most systems but would still be good to fix.
+
+ * Bugfix: reverted the value of `ms:i` to pre-2.22 versions (#1146). This was
+   an oversight. See fcd4df2 for details.
+
+Notable changes to paftools.js and mappy:
+
+ * New feature: expose `bw_long` to mappy's Aligner class (#1124).
+
+ * Bugfix: fixed several compatibility issues with k8 v1.0 (#1161 and #1166).
+   Subcommands "call", "pbsim2fq" and "mason2fq" were not working with v1.0.
+
+Minimap2 should output identical alignments to v2.26, except the ms tag.
+
+(2.27: 12 March 2024, r1193)
+
+
+
+Release 2.26-r1175 (29 April 2023)
+----------------------------------
+
+Fixed the broken Python package. This is the only change.
+
+(2.26: 25 April 2023, r1173)
+
+
+
+Release 2.25-r1173 (25 April 2023)
+----------------------------------
+
+Notable changes:
+
+ * Improvement: use the miniprot splice model for RNA-seq alignment by default.
+   This model considers non-GT-AG splice sites and leads to slightly higher
+   (<0.1%) accuracy and sensitivity on real human data.
+
+ * Change: increased the default `-I` to `8G` such that minimap2 would create a
+   uni-part index for a pair of mammalian genomes. This change may increase the
+   memory for all-vs-all read overlap alignment given large datasets.
+
+ * New feature: output the sequences in secondary alignments with option
+   `--secondary-seq` (#687).
+
+ * Bugfix: --rmq was not parsed correctly (#1010)
+
+ * Bugfix: possibly incorrect coordinate when applying end bonus to the target
+   sequence (#1025). This is a ksw2 bug. It does not affect minimap2 as
+   minimap2 is not using the affected feature.
+
+ * Improvement: incorporated several changes for better compatibility with
+   Windows (#1051) and for minimap2 integration at Oxford Nanopore Technologies
+   (#1048 and #1033).
+
+ * Improvement: output the HD-line in SAM output (#1019).
+
+ * Improvement: check minimap2 index file in mappy to prevent segmentation
+   fault for certain indices (#1008).
+
+For genomic sequences, minimap2 should give identical output to v2.24.
+Long-read RNA-seq alignment may occasionally differ from previous versions.
+
+(2.25: 25 April 2023, r1173)
+
+
+
 Release 2.24-r1122 (26 December 2021)
 -------------------------------------
 

README.md

@@ -15,7 +15,7 @@ cd minimap2 && make
 ./minimap2 -ax map-pb ref.fa pacbio.fq.gz > aln.sam       # PacBio CLR genomic reads
 ./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam         # Oxford Nanopore genomic reads
 ./minimap2 -ax map-hifi ref.fa pacbio-ccs.fq.gz > aln.sam # PacBio HiFi/CCS genomic reads (v2.19 or later)
-./minimap2 -ax asm20 ref.fa pacbio-ccs.fq.gz > aln.sam    # PacBio HiFi/CCS genomic reads (v2.18 or earlier)
+./minimap2 -ax lr:hq ref.fa ont-Q20.fq.gz > aln.sam       # Nanopore Q20 genomic reads (v2.27 or later)
 ./minimap2 -ax sr ref.fa read1.fa read2.fa > aln.sam      # short genomic paired-end reads
 ./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads (strand unknown)
 ./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam  # noisy Nanopore Direct RNA-seq
@@ -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.24/minimap2-2.24_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.24_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
@@ -139,12 +139,15 @@ parameters at the same time. The default setting is the same as `map-ont`.
 ```sh
 minimap2 -ax map-pb  ref.fa pacbio-reads.fq > aln.sam   # for PacBio CLR reads
 minimap2 -ax map-ont ref.fa ont-reads.fq > aln.sam      # for Oxford Nanopore reads
+minimap2 -ax map-iclr ref.fa iclr-reads.fq > aln.sam    # for Illumina Complete Long Reads
 ```
 The difference between `map-pb` and `map-ont` is that `map-pb` uses
 homopolymer-compressed (HPC) minimizers as seeds, while `map-ont` uses ordinary
-minimizers as seeds. Emperical evaluation suggests HPC minimizers improve
+minimizers as seeds. Empirical evaluation suggests HPC minimizers improve
 performance and sensitivity when aligning PacBio CLR reads, but hurt when aligning
-Nanopore reads.
+Nanopore reads. `map-iclr` uses an adjusted alignment scoring matrix that
+accounts for the low overall error rate in the reads, with transversion errors
+being less frequent than transitions.
 
 #### <a name="map-long-splice"></a>Map long mRNA/cDNA reads
 
@@ -350,6 +353,11 @@ If you use minimap2 in your work, please cite:
 > Li, H. (2018). Minimap2: pairwise alignment for nucleotide sequences.
 > *Bioinformatics*, **34**:3094-3100. [doi:10.1093/bioinformatics/bty191][doi]
 
+and/or:
+
+> Li, H. (2021). New strategies to improve minimap2 alignment accuracy.
+> *Bioinformatics*, **37**:4572-4574. [doi:10.1093/bioinformatics/btab705][doi2]
+
 ## <a name="dguide"></a>Developers' Guide
 
 Minimap2 is not only a command line tool, but also a programming library.
@@ -399,5 +407,6 @@ mappy` or [from BioConda][mappyconda] via `conda install -c bioconda mappy`.
 [manpage]: https://lh3.github.io/minimap2/minimap2.html
 [manpage-cs]: https://lh3.github.io/minimap2/minimap2.html#10
 [doi]: https://doi.org/10.1093/bioinformatics/bty191
-[smide]: https://github.com/nemequ/simde
+[doi2]: https://doi.org/10.1093/bioinformatics/btab705
+[simde]: https://github.com/nemequ/simde
 [unimap]: https://github.com/lh3/unimap

align.c

@@ -21,6 +21,18 @@ static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t sc
 		mat[(m - 1) * m + j] = sc_ambi;
 }
 
+static void ksw_gen_ts_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t transition, int8_t sc_ambi)
+{
+	assert(m == 5);
+	ksw_gen_simple_mat(m, mat, a, b, sc_ambi);
+	if (transition == 0 || transition == b) return;
+	transition = transition > 0? -transition : transition;
+	mat[0 * m + 2] = transition;  // A->G
+	mat[1 * m + 3] = transition;  // C->T
+	mat[2 * m + 0] = transition;  // G->A
+	mat[3 * m + 1] = transition;  // T->C
+}
+
 static inline void mm_seq_rev(uint32_t len, uint8_t *seq)
 {
 	uint32_t i;
@@ -283,7 +295,7 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
 			toff += len;
 		}
 	}
-	p->dp_max = (int32_t)(max + .499);
+	p->dp_max = p->dp_max0 = (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
 }
@@ -323,12 +335,16 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 		for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr);
 		fputc('\n', stderr);
 	}
+	if (opt->transition != 0 && opt->b != opt->transition)
+		flag |= KSW_EZ_GENERIC_SC;
 	if (opt->max_sw_mat > 0 && (int64_t)tlen * qlen > opt->max_sw_mat) {
 		ksw_reset_extz(ez);
 		ez->zdropped = 1;
-	} else if (opt->flag & MM_F_SPLICE)
+	} else if (opt->flag & MM_F_SPLICE) {
-		ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, opt->junc_bonus, flag, junc, ez);
+		int flag_tmp = flag;
-	else if (opt->q == opt->q2 && opt->e == opt->e2)
+		if (!(opt->flag & MM_F_SPLICE_OLD)) flag_tmp |= KSW_EZ_SPLICE_CMPLX;
+		ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, opt->junc_bonus, flag_tmp, junc, ez);
+	} else if (opt->q == opt->q2 && opt->e == opt->e2)
 		ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, flag, ez);
 	else
 		ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez);
@@ -584,7 +600,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 
 	r2->cnt = 0;
 	if (r->cnt == 0) return;
-	ksw_gen_simple_mat(5, mat, opt->a, opt->b, opt->sc_ambi);
+	ksw_gen_ts_mat(5, mat, opt->a, opt->b, opt->transition, opt->sc_ambi);
 	bw = (int)(opt->bw * 1.5 + 1.);
 	bw_long = (int)(opt->bw_long * 1.5 + 1.);
 	if (bw_long < bw) bw_long = bw;
@@ -842,7 +858,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	if (ql < opt->min_chain_score || ql > opt->max_gap) return 0;
 	if (tl < opt->min_chain_score || tl > opt->max_gap) return 0;
 
-	ksw_gen_simple_mat(5, mat, opt->a, opt->b, opt->sc_ambi);
+	ksw_gen_ts_mat(5, mat, opt->a, opt->b, opt->transition, opt->sc_ambi);
 	tseq = (uint8_t*)kmalloc(km, tl);
 	mm_idx_getseq(mi, r1->rid, r1->re, r2->rs, tseq);
 	qseq = r1->rev? &qseq0[0][r2->qe] : &qseq0[1][qlen - r2->qs];
@@ -917,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.24/minimap2-2.24_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.24_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

@@ -119,6 +119,7 @@ int mm_write_sam_hdr(const mm_idx_t *idx, const char *rg, const char *ver, int a
 {
 	kstring_t str = {0,0,0};
 	int ret = 0;
+	mm_sprintf_lite(&str, "@HD\tVN:1.6\tSO:unsorted\tGO:query\n");
 	if (idx) {
 		uint32_t i;
 		for (i = 0; i < idx->n_seq; ++i)
@@ -138,10 +139,48 @@ int mm_write_sam_hdr(const mm_idx_t *idx, const char *rg, const char *ver, int a
 	return ret;
 }
 
-static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp, int no_iden, int write_tag)
+static void write_indel_ds(kstring_t *str, int64_t len, const uint8_t *seq, int64_t ll, int64_t lr) // write an indel to ds; adapted from minigraph
 {
-	int i, q_off, t_off;
+	int64_t i;
-	if (write_tag) mm_sprintf_lite(s, "\tcs:Z:");
+	if (ll + lr >= len) {
+		mm_sprintf_lite(str, "[");
+		for (i = 0; i < len; ++i)
+			mm_sprintf_lite(str, "%c", "acgtn"[seq[i]]);
+		mm_sprintf_lite(str, "]");
+	} else {
+		int64_t k = 0;
+		if (ll > 0) {
+			mm_sprintf_lite(str, "[");
+			for (i = 0; i < ll; ++i)
+				mm_sprintf_lite(str, "%c", "acgtn"[seq[k+i]]);
+			mm_sprintf_lite(str, "]");
+			k += ll;
+		}
+		for (i = 0; i < len - lr - ll; ++i)
+			mm_sprintf_lite(str, "%c", "acgtn"[seq[k+i]]);
+		k += len - lr - ll;
+		if (lr > 0) {
+			mm_sprintf_lite(str, "[");
+			for (i = 0; i < lr; ++i)
+				mm_sprintf_lite(str, "%c", "acgtn"[seq[k+i]]);
+			mm_sprintf_lite(str, "]");
+		}
+	}
+}
+
+static void write_cs_ds_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp, int no_iden, int is_ds, int write_tag)
+{
+	int i, q_off, t_off, q_len = 0, t_len = 0;
+	if (write_tag) mm_sprintf_lite(s, "\t%cs:Z:", is_ds? 'd' : 'c');
+	for (i = 0; i < (int)r->p->n_cigar; ++i) {
+		int op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
+		if (op == MM_CIGAR_MATCH || op == MM_CIGAR_EQ_MATCH || op == MM_CIGAR_X_MISMATCH)
+			q_len += len, t_len += len;
+		else if (op == MM_CIGAR_INS)
+			q_len += len;
+		else if (op == MM_CIGAR_DEL || op == MM_CIGAR_N_SKIP)
+			t_len += len;
+	}
 	for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
 		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
 		assert((op >= MM_CIGAR_MATCH && op <= MM_CIGAR_N_SKIP) || op == MM_CIGAR_EQ_MATCH || op == MM_CIGAR_X_MISMATCH);
@@ -167,14 +206,42 @@ static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 			}
 			q_off += len, t_off += len;
 		} else if (op == MM_CIGAR_INS) {
-			for (j = 0, tmp[len] = 0; j < len; ++j)
+			if (is_ds) {
-				tmp[j] = "acgtn"[qseq[q_off + j]];
+				int z, ll, lr, y = q_off;
-			mm_sprintf_lite(s, "+%s", tmp);
+				for (z = 1; z <= len; ++z)
+					if (y - z < 0 || qseq[y + len - z] != qseq[y - z])
+						break;
+				lr = z - 1;
+				for (z = 0; z < len; ++z)
+					if (y + len + z >= q_len || qseq[y + len + z] != qseq[y + z])
+						break;
+				ll = z;
+				mm_sprintf_lite(s, "+");
+				write_indel_ds(s, len, &qseq[y], ll, lr);
+			} else {
+				for (j = 0, tmp[len] = 0; j < len; ++j)
+					tmp[j] = "acgtn"[qseq[q_off + j]];
+				mm_sprintf_lite(s, "+%s", tmp);
+			}
 			q_off += len;
 		} else if (op == MM_CIGAR_DEL) {
-			for (j = 0, tmp[len] = 0; j < len; ++j)
+			if (is_ds) {
-				tmp[j] = "acgtn"[tseq[t_off + j]];
+				int z, ll, lr, x = t_off;
-			mm_sprintf_lite(s, "-%s", tmp);
+				for (z = 1; z <= len; ++z)
+					if (x - z < 0 || tseq[x + len - z] != tseq[x - z])
+						break;
+				lr = z - 1;
+				for (z = 0; z < len; ++z)
+					if (x + len + z >= t_len || tseq[x + z] != tseq[x + len + z])
+						break;
+				ll = z;
+				mm_sprintf_lite(s, "-");
+				write_indel_ds(s, len, &tseq[x], ll, lr);
+			} else {
+				for (j = 0, tmp[len] = 0; j < len; ++j)
+					tmp[j] = "acgtn"[tseq[t_off + j]];
+				mm_sprintf_lite(s, "-%s", tmp);
+			}
 			t_off += len;
 		} else { // intron
 			assert(len >= 2);
@@ -217,7 +284,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, int is_qstrand)
+static void write_cs_ds_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 is_ds, int write_tag, int is_qstrand)
 {
 	extern unsigned char seq_nt4_table[256];
 	int i;
@@ -244,7 +311,7 @@ static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_
 		}
 	}
 	if (is_MD) write_MD_core(s, tseq, qseq, r, tmp, write_tag);
-	else write_cs_core(s, tseq, qseq, r, tmp, no_iden, write_tag);
+	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);
 }
 
@@ -255,7 +322,7 @@ int mm_gen_cs_or_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, cons
 	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, is_qstrand);
+	write_cs_ds_or_MD(km, &str, mi, &t, r, no_iden, is_MD, 0, 0, is_qstrand);
 	*max_len = str.m;
 	*buf = str.s;
 	return str.l;
@@ -277,7 +344,7 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 	if (r->id == r->parent) type = r->inv? 'I' : 'P';
 	else type = r->inv? 'i' : 'S';
 	if (r->p) {
-		mm_sprintf_lite(s, "\tNM:i:%d\tms:i:%d\tAS:i:%d\tnn:i:%d", r->blen - r->mlen + r->p->n_ambi, r->p->dp_max, r->p->dp_score, r->p->n_ambi);
+		mm_sprintf_lite(s, "\tNM:i:%d\tms:i:%d\tAS:i:%d\tnn:i:%d", r->blen - r->mlen + r->p->n_ambi, r->p->dp_max0, r->p->dp_score, r->p->n_ambi);
 		if (r->p->trans_strand == 1 || r->p->trans_strand == 2)
 			mm_sprintf_lite(s, "\tts:A:%c", "?+-?"[r->p->trans_strand]);
 	}
@@ -325,8 +392,8 @@ void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const
 		for (k = 0; k < r->p->n_cigar; ++k)
 			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)))
+	if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_DS|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, !!(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);
 }
@@ -369,14 +436,16 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
 		clip_len[0] = r->rev? qlen - r->qe : r->qs;
 		clip_len[1] = r->rev? r->qs : qlen - r->qe;
 		if (in_tag) {
-			int clip_char = (sam_flag&0x800) && !(opt_flag&MM_F_SOFTCLIP)? 5 : 4;
+			int clip_char = (((sam_flag&0x800) || ((sam_flag&0x100) && (opt_flag&MM_F_SECONDARY_SEQ))) &&
+							 !(opt_flag&MM_F_SOFTCLIP)) ? 5 : 4;
 			mm_sprintf_lite(s, "\tCG:B:I");
 			if (clip_len[0]) mm_sprintf_lite(s, ",%u", clip_len[0]<<4|clip_char);
 			for (k = 0; k < r->p->n_cigar; ++k)
 				mm_sprintf_lite(s, ",%u", r->p->cigar[k]);
 			if (clip_len[1]) mm_sprintf_lite(s, ",%u", clip_len[1]<<4|clip_char);
 		} else {
-			int clip_char = (sam_flag&0x800) && !(opt_flag&MM_F_SOFTCLIP)? 'H' : 'S';
+			int clip_char = (((sam_flag&0x800) || ((sam_flag&0x100) && (opt_flag&MM_F_SECONDARY_SEQ))) &&
+							 !(opt_flag&MM_F_SOFTCLIP)) ? 'H' : 'S';
 			assert(clip_len[0] < qlen && clip_len[1] < qlen);
 			if (clip_len[0]) mm_sprintf_lite(s, "%d%c", clip_len[0], clip_char);
 			for (k = 0; k < r->p->n_cigar; ++k)
@@ -451,7 +520,7 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 		if (cigar_in_tag) {
 			int slen;
 			if ((flag & 0x900) == 0 || (opt_flag & MM_F_SOFTCLIP)) slen = t->l_seq;
-			else if (flag & 0x100) slen = 0;
+			else if ((flag & 0x100) && !(opt_flag & MM_F_SECONDARY_SEQ)) slen = 0;
 			else slen = r->qe - r->qs;
 			mm_sprintf_lite(s, "%dS%dN", slen, r->re - r->rs);
 		} else write_sam_cigar(s, flag, 0, t->l_seq, r, opt_flag);
@@ -492,7 +561,7 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 			mm_sprintf_lite(s, "\t");
 			if (t->qual) sam_write_sq(s, t->qual, t->l_seq, r->rev, 0);
 			else mm_sprintf_lite(s, "*");
-		} else if (flag & 0x100) {
+		} else if ((flag & 0x100) && !(opt_flag & MM_F_SECONDARY_SEQ)){
 			mm_sprintf_lite(s, "*\t*");
 		} else {
 			sam_write_sq(s, t->seq + r->qs, r->qe - r->qs, r->rev, r->rev);
@@ -532,8 +601,8 @@ 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)))
+		if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_DS|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, 0);
+			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, 0);
 		if (cigar_in_tag)
 			write_sam_cigar(s, flag, 1, t->l_seq, r, opt_flag);
 	}

index.c

@@ -192,6 +192,7 @@ int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f)
 	if (f <= 0.) return INT32_MAX;
 	for (i = 0; i < 1<<mi->b; ++i)
 		if (mi->B[i].h) n += kh_size((idxhash_t*)mi->B[i].h);
+	if (n == 0) return INT32_MAX;
 	a = (uint32_t*)malloc(n * 4);
 	for (i = n = 0; i < 1<<mi->b; ++i) {
 		idxhash_t *h = (idxhash_t*)mi->B[i].h;

kalloc.c

@@ -40,7 +40,8 @@ void *km_init2(void *km_par, size_t min_core_size)
 	kmem_t *km;
 	km = (kmem_t*)kcalloc(km_par, 1, sizeof(kmem_t));
 	km->par = km_par;
-	km->min_core_size = min_core_size > 0? min_core_size : 0x80000;
+	if (km_par) km->min_core_size = min_core_size > 0? min_core_size : ((kmem_t*)km_par)->min_core_size - 2;
+	else km->min_core_size = min_core_size > 0? min_core_size : 0x80000;
 	return (void*)km;
 }
 
@@ -183,6 +184,16 @@ void *krealloc(void *_km, void *ap, size_t n_bytes) // TODO: this can be made mo
 	return q;
 }
 
+void *krelocate(void *km, void *ap, size_t n_bytes)
+{
+	void *p;
+	if (km == 0 || ap == 0) return ap;
+	p = kmalloc(km, n_bytes);
+	memcpy(p, ap, n_bytes);
+	kfree(km, ap);
+	return p;
+}
+
 void km_stat(const void *_km, km_stat_t *s)
 {
 	kmem_t *km = (kmem_t*)_km;
@@ -203,3 +214,11 @@ void km_stat(const void *_km, km_stat_t *s)
 		s->largest = s->largest > size? s->largest : size;
 	}
 }
+
+void km_stat_print(const void *km)
+{
+	km_stat_t st;
+	km_stat(km, &st);
+	fprintf(stderr, "[km_stat] cap=%ld, avail=%ld, largest=%ld, n_core=%ld, n_block=%ld\n",
+			st.capacity, st.available, st.largest, st.n_blocks, st.n_cores);
+}

kalloc.h

@@ -13,6 +13,7 @@ typedef struct {
 
 void *kmalloc(void *km, size_t size);
 void *krealloc(void *km, void *ptr, size_t size);
+void *krelocate(void *km, void *ap, size_t n_bytes);
 void *kcalloc(void *km, size_t count, size_t size);
 void kfree(void *km, void *ptr);
 
@@ -20,11 +21,21 @@ void *km_init(void);
 void *km_init2(void *km_par, size_t min_core_size);
 void km_destroy(void *km);
 void km_stat(const void *_km, km_stat_t *s);
+void km_stat_print(const void *km);
 
 #ifdef __cplusplus
 }
 #endif
 
+#define Kmalloc(km, type, cnt)       ((type*)kmalloc((km), (cnt) * sizeof(type)))
+#define Kcalloc(km, type, cnt)       ((type*)kcalloc((km), (cnt), sizeof(type)))
+#define Krealloc(km, type, ptr, cnt) ((type*)krealloc((km), (ptr), (cnt) * sizeof(type)))
+
+#define Kexpand(km, type, a, m) do { \
+		(m) = (m) >= 4? (m) + ((m)>>1) : 16; \
+		(a) = Krealloc(km, type, (a), (m)); \
+	} while (0)
+
 #define KMALLOC(km, ptr, len) ((ptr) = (__typeof__(ptr))kmalloc((km), (len) * sizeof(*(ptr))))
 #define KCALLOC(km, ptr, len) ((ptr) = (__typeof__(ptr))kcalloc((km), (len), sizeof(*(ptr))))
 #define KREALLOC(km, ptr, len) ((ptr) = (__typeof__(ptr))krealloc((km), (ptr), (len) * sizeof(*(ptr))))
@@ -50,7 +61,7 @@ void km_stat(const void *_km, km_stat_t *s);
 	} kmp_##name##_t; \
 	SCOPE kmp_##name##_t *kmp_init_##name(void *km) { \
 		kmp_##name##_t *mp; \
-		KCALLOC(km, mp, 1); \
+		mp = Kcalloc(km, kmp_##name##_t, 1); \
 		mp->km = km; \
 		return mp; \
 	} \
@@ -66,7 +77,7 @@ void km_stat(const void *_km, km_stat_t *s);
 	} \
 	SCOPE void kmp_free_##name(kmp_##name##_t *mp, kmptype_t *p) { \
 		--mp->cnt; \
-		if (mp->n == mp->max) KEXPAND(mp->km, mp->buf, mp->max); \
+		if (mp->n == mp->max) Kexpand(mp->km, kmptype_t*, mp->buf, mp->max); \
 		mp->buf[mp->n++] = p; \
 	}
 

ksw2.h

@@ -15,6 +15,7 @@
 #define KSW_EZ_SPLICE_FOR  0x100
 #define KSW_EZ_SPLICE_REV  0x200
 #define KSW_EZ_SPLICE_FLANK 0x400
+#define KSW_EZ_SPLICE_CMPLX 0x800
 
 // The subset of CIGAR operators used by ksw code.
 // Use MM_CIGAR_* from minimap.h if you need the full list.

ksw2_extd2_sse.c

@@ -358,7 +358,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 			} else H[0] = v8[0] - qe, max_H = H[0], max_t = 0; // special casing r==0
 			// update ez
 			if (en0 == tlen - 1 && H[en0] > ez->mte)
-				ez->mte = H[en0], ez->mte_q = r - en;
+				ez->mte = H[en0], ez->mte_q = r - en0;
 			if (r - st0 == qlen - 1 && H[st0] > ez->mqe)
 				ez->mqe = H[st0], ez->mqe_t = st0;
 			if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, e2)) break;

ksw2_exts2_sse.c

@@ -71,6 +71,7 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 
 	ksw_reset_extz(ez);
 	if (m <= 1 || qlen <= 0 || tlen <= 0 || q2 <= q + e) return;
+	assert((flag & KSW_EZ_SPLICE_FOR) == 0 || (flag & KSW_EZ_SPLICE_REV) == 0); // can't be both set
 
 	zero_   = _mm_set1_epi8(0);
 	q_      = _mm_set1_epi8(q);
@@ -118,55 +119,93 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 
 	// set the donor and acceptor arrays. TODO: this assumes 0/1/2/3 encoding!
 	if (flag & (KSW_EZ_SPLICE_FOR|KSW_EZ_SPLICE_REV)) {
-		int semi_cost = flag&KSW_EZ_SPLICE_FLANK? -noncan/2 : 0; // GTr or yAG is worth 0.5 bit; see PMID:18688272
+		const int sp0[4] = { 8, 15, 21, 30 };
-		memset(donor, -noncan, tlen_ * 16);
+		int sp[4];
-		memset(acceptor, -noncan, tlen_ * 16);
+		if (flag & KSW_EZ_SPLICE_CMPLX) {
+			for (t = 0; t < 4; ++t)
+				sp[t] = (int)((double)sp0[t] / 3. + .499);
+		} else {
+			sp[0] = flag&KSW_EZ_SPLICE_FLANK? noncan / 2 : 0;
+			sp[1] = sp[2] = sp[3] = noncan;
+		}
+		memset(donor,    -sp[3], tlen_ * 16);
+		memset(acceptor, -sp[3], tlen_ * 16);
 		if (!(flag & KSW_EZ_REV_CIGAR)) {
 			for (t = 0; t < tlen - 4; ++t) {
-				int can_type = 0; // type of canonical site: 0=none, 1=GT/AG only, 2=GTr/yAG
+				int z = 3;
-				if ((flag & KSW_EZ_SPLICE_FOR) && target[t+1] == 2 && target[t+2] == 3) can_type = 1; // GTr...
+				if (flag & KSW_EZ_SPLICE_FOR) {
-				if ((flag & KSW_EZ_SPLICE_REV) && target[t+1] == 1 && target[t+2] == 3) can_type = 1; // CTr...
+					if (target[t+1] == 2 && target[t+2] == 3)             // |GT.
-				if (can_type && (target[t+3] == 0 || target[t+3] == 2)) can_type = 2;
+						z = target[t+3] == 0 || target[t+3] == 2? -1 : 0; // |GTr or not
-				if (can_type) ((int8_t*)donor)[t] = can_type == 2? 0 : semi_cost;
+					else if (target[t+1] == 2 && target[t+2] == 1) z = 1; // |GC.
+					else if (target[t+1] == 0 && target[t+2] == 3) z = 2; // |AT.
+				} else if (flag & KSW_EZ_SPLICE_REV) {
+					if (target[t+1] == 1 && target[t+2] == 3)             // |CT. (revcomp of .AG|)
+						z = target[t+3] == 0 || target[t+3] == 2? -1 : 0;
+					else if (target[t+1] == 2 && target[t+2] == 3) z = 2; // |GT. (revcomp of .AC|)
+				}
+				((int8_t*)donor)[t] = z < 0? 0 : -sp[z];
 			}
-			if (junc)
-				for (t = 0; t < tlen - 1; ++t)
-					if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t+1]&1)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t+1]&8)))
-						((int8_t*)donor)[t] += junc_bonus;
 			for (t = 2; t < tlen; ++t) {
-				int can_type = 0;
+				int z = 3;
-				if ((flag & KSW_EZ_SPLICE_FOR) && target[t-1] == 0 && target[t] == 2) can_type = 1; // ...yAG
+				if (flag & KSW_EZ_SPLICE_FOR) {
-				if ((flag & KSW_EZ_SPLICE_REV) && target[t-1] == 0 && target[t] == 1) can_type = 1; // ...yAC
+					if (target[t-1] == 0 && target[t] == 2)               // .AG|
-				if (can_type && (target[t-2] == 1 || target[t-2] == 3)) can_type = 2;
+						z = target[t-2] == 1 || target[t-2] == 3? -1 : 0; // yAG| or not
-				if (can_type) ((int8_t*)acceptor)[t] = can_type == 2? 0 : semi_cost;
+					else if (target[t-1] == 0 && target[t] == 1) z = 2;   // .AC|
+				} else if (flag & KSW_EZ_SPLICE_REV) {
+					if (target[t-1] == 0 && target[t] == 1)               // .AC| (revcomp of |GT.)
+						z = target[t-2] == 1 || target[t-2] == 3? -1 : 0; // yAC| or not
+					else if (target[t-1] == 2 && target[t] == 1) z = 1;   // .GC| (revcomp of |GC.)
+					else if (target[t-1] == 0 && target[t] == 3) z = 2;   // .AT| (revcomp of |AT.)
+				}
+				((int8_t*)acceptor)[t] = z < 0? 0 : -sp[z];
 			}
-			if (junc)
-				for (t = 0; t < tlen; ++t)
-					if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t]&2)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t]&4)))
-						((int8_t*)acceptor)[t] += junc_bonus;
 		} else {
 			for (t = 0; t < tlen - 4; ++t) {
-				int can_type = 0; // type of canonical site: 0=none, 1=GT/AG only, 2=GTr/yAG
+				int z = 3;
-				if ((flag & KSW_EZ_SPLICE_FOR) && target[t+1] == 2 && target[t+2] == 0) can_type = 1; // GAy...
+				if (flag & KSW_EZ_SPLICE_FOR) {
-				if ((flag & KSW_EZ_SPLICE_REV) && target[t+1] == 1 && target[t+2] == 0) can_type = 1; // CAy...
+					if (target[t+1] == 2 && target[t+2] == 0)             // |GA. (rev of .AG|)
-				if (can_type && (target[t+3] == 1 || target[t+3] == 3)) can_type = 2;
+						z = target[t+3] == 1 || target[t+3] == 3? -1 : 0;
-				if (can_type) ((int8_t*)donor)[t] = can_type == 2? 0 : semi_cost;
+					else if (target[t+1] == 1 && target[t+2] == 0) z = 2; // |CA. (rev of .AC|)
+				} else if (flag & KSW_EZ_SPLICE_REV) {
+					if (target[t+1] == 1 && target[t+2] == 0)             // |CA. (comp of |GT.)
+						z = target[t+3] == 1 || target[t+3] == 3? -1 : 0;
+					else if (target[t+1] == 1 && target[t+2] == 2) z = 1; // |CG. (comp of |GC.)
+					else if (target[t+1] == 3 && target[t+2] == 0) z = 2; // |TA. (comp of |AT.)
+				}
+				((int8_t*)donor)[t] = z < 0? 0 : -sp[z];
 			}
-			if (junc)
-				for (t = 0; t < tlen - 1; ++t)
-					if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t+1]&2)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t+1]&4)))
-						((int8_t*)donor)[t] += junc_bonus;
 			for (t = 2; t < tlen; ++t) {
-				int can_type = 0;
+				int z = 3;
-				if ((flag & KSW_EZ_SPLICE_FOR) && target[t-1] == 3 && target[t] == 2) can_type = 1; // ...rTG
+				if (flag & KSW_EZ_SPLICE_FOR) {
-				if ((flag & KSW_EZ_SPLICE_REV) && target[t-1] == 3 && target[t] == 1) can_type = 1; // ...rTC
+					if (target[t-1] == 3 && target[t] == 2)               // .TG| (rev of |GT.)
-				if (can_type && (target[t-2] == 0 || target[t-2] == 2)) can_type = 2;
+						z = target[t-2] == 0 || target[t-2] == 2? -1 : 0;
-				if (can_type) ((int8_t*)acceptor)[t] = can_type == 2? 0 : semi_cost;
+					else if (target[t-1] == 1 && target[t] == 2) z = 1;   // .CG| (rev of |GC.)
+					else if (target[t-1] == 3 && target[t] == 0) z = 2;   // .TA| (rev of |AT.)
+				} else if (flag & KSW_EZ_SPLICE_REV) {
+					if (target[t-1] == 3 && target[t] == 1)               // .TC| (comp of .AG|)
+						z = target[t-2] == 0 || target[t-2] == 2? -1 : 0;
+					else if (target[t-1] == 3 && target[t] == 2) z = 2;   // .TG| (comp of .AC|)
+				}
+				((int8_t*)acceptor)[t] = z < 0? 0 : -sp[z];
 			}
-			if (junc)
+		}
-				for (t = 0; t < tlen; ++t)
+	}
-					if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t]&1)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t]&8)))
+
-						((int8_t*)acceptor)[t] += junc_bonus;
+	if (junc) {
+		if (!(flag & KSW_EZ_REV_CIGAR)) {
+			for (t = 0; t < tlen - 1; ++t)
+				if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t+1]&1)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t+1]&8)))
+					((int8_t*)donor)[t] += junc_bonus;
+			for (t = 0; t < tlen; ++t)
+				if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t]&2)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t]&4)))
+					((int8_t*)acceptor)[t] += junc_bonus;
+		} else {
+			for (t = 0; t < tlen - 1; ++t)
+				if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t+1]&2)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t+1]&4)))
+					((int8_t*)donor)[t] += junc_bonus;
+			for (t = 0; t < tlen; ++t)
+				if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t]&1)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t]&8)))
+					((int8_t*)acceptor)[t] += junc_bonus;
 		}
 	}
 
@@ -376,7 +415,7 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 			} else H[0] = v8[0] - qe, max_H = H[0], max_t = 0; // special casing r==0
 			// update ez
 			if (en0 == tlen - 1 && H[en0] > ez->mte)
-				ez->mte = H[en0], ez->mte_q = r - en;
+				ez->mte = H[en0], ez->mte_q = r - en0;
 			if (r - st0 == qlen - 1 && H[st0] > ez->mqe)
 				ez->mqe = H[st0], ez->mqe_t = st0;
 			if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, 0)) break;

ksw2_extz2_sse.c

@@ -269,7 +269,7 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 			} else H[0] = v8[0] - qe - qe, max_H = H[0], max_t = 0; // special casing r==0
 			// update ez
 			if (en0 == tlen - 1 && H[en0] > ez->mte)
-				ez->mte = H[en0], ez->mte_q = r - en;
+				ez->mte = H[en0], ez->mte_q = r - en0;
 			if (r - st0 == qlen - 1 && H[st0] > ez->mqe)
 				ez->mqe = H[st0], ez->mqe_t = st0;
 			if (ksw_apply_zdrop(ez, 1, max_H, r, max_t, zdrop, e)) break;

lchain.c

@@ -35,7 +35,7 @@ uint64_t *mg_chain_backtrack(void *km, int64_t n, const int32_t *f, const int64_
 	for (i = 0, n_z = 0; i < n; ++i) // precompute n_z
 		if (f[i] >= min_sc) ++n_z;
 	if (n_z == 0) return 0;
-	KMALLOC(km, z, n_z);
+	z = Kmalloc(km, mm128_t, n_z);
 	for (i = 0, k = 0; i < n; ++i) // populate z[]
 		if (f[i] >= min_sc) z[k].x = f[i], z[k++].y = i;
 	radix_sort_128x(z, z + n_z);
@@ -54,7 +54,7 @@ uint64_t *mg_chain_backtrack(void *km, int64_t n, const int32_t *f, const int64_
 			else n_v = n_v0;
 		}
 	}
-	KMALLOC(km, u, n_u);
+	u = Kmalloc(km, uint64_t, n_u);
 	memset(t, 0, n * 4);
 	for (k = n_z - 1, n_v = n_u = 0; k >= 0; --k) { // populate u[]
 		if (t[z[k].y] == 0) {
@@ -82,7 +82,7 @@ static mm128_t *compact_a(void *km, int32_t n_u, uint64_t *u, int32_t n_v, int32
 	int64_t i, j, k;
 
 	// write the result to b[]
-	KMALLOC(km, b, n_v);
+	b = Kmalloc(km, mm128_t, n_v);
 	for (i = 0, k = 0; i < n_u; ++i) {
 		int32_t k0 = k, ni = (int32_t)u[i];
 		for (j = 0; j < ni; ++j)
@@ -91,13 +91,13 @@ static mm128_t *compact_a(void *km, int32_t n_u, uint64_t *u, int32_t n_v, int32
 	kfree(km, v);
 
 	// sort u[] and a[] by the target position, such that adjacent chains may be joined
-	KMALLOC(km, w, n_u);
+	w = Kmalloc(km, mm128_t, n_u);
 	for (i = k = 0; i < n_u; ++i) {
 		w[i].x = b[k].x, w[i].y = (uint64_t)k<<32|i;
 		k += (int32_t)u[i];
 	}
 	radix_sort_128x(w, w + n_u);
-	KMALLOC(km, u2, n_u);
+	u2 = Kmalloc(km, uint64_t, n_u);
 	for (i = k = 0; i < n_u; ++i) {
 		int32_t j = (int32_t)w[i].y, n = (int32_t)u[j];
 		u2[i] = u[j];
@@ -138,7 +138,7 @@ static inline int32_t comput_sc(const mm128_t *ai, const mm128_t *aj, int32_t ma
 }
 
 /* Input:
- *   a[].x: tid<<33 | rev<<32 | tpos
+ *   a[].x: rev<<63 | tid<<32 | tpos
  *   a[].y: flags<<40 | q_span<<32 | q_pos
  * Output:
  *   n_u: #chains
@@ -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;
@@ -160,10 +160,10 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
 	if (max_dist_x < bw) max_dist_x = bw;
 	if (max_dist_y < bw && !is_cdna) max_dist_y = bw;
 	if (is_cdna) max_drop = INT32_MAX;
-	KMALLOC(km, p, n);
+	p = Kmalloc(km, int64_t, n);
-	KMALLOC(km, f, n);
+	f = Kmalloc(km, int32_t, n);
-	KMALLOC(km, v, n);
+	v = Kmalloc(km, int32_t, n);
-	KCALLOC(km, t, n);
+	t = Kcalloc(km, int32_t, n);
 
 	// fill the score and backtrack arrays
 	for (i = 0, max_ii = -1; i < n; ++i) {
@@ -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);
@@ -251,7 +251,7 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 					   int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km)
 {
 	int32_t *f,*t, *v, n_u, n_v, mmax_f = 0, max_rmq_size = 0, max_drop = bw;
-	int64_t *p, i, i0, st = 0, st_inner = 0, n_iter = 0;
+	int64_t *p, i, i0, st = 0, st_inner = 0;
 	uint64_t *u;
 	lc_elem_t *root = 0, *root_inner = 0;
 	void *mem_mp = 0;
@@ -263,11 +263,12 @@ 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;
-	KMALLOC(km, p, n);
+	if (max_dist_inner > max_dist) max_dist_inner = max_dist;
-	KMALLOC(km, f, n);
+	p = Kmalloc(km, int64_t, n);
-	KCALLOC(km, t, n);
+	f = Kmalloc(km, int32_t, n);
-	KMALLOC(km, v, n);
+	t = Kcalloc(km, int32_t, n);
+	v = Kmalloc(km, int32_t, n);
 	mem_mp = km_init2(km, 0x10000);
 	mp = kmp_init_rmq(mem_mp);
 
@@ -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) {
@@ -345,7 +345,6 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 						}
 						if (!krmq_itr_prev(lc_elem, &itr)) break;
 					}
-					n_iter += n_rmq_iter;
 				}
 			}
 		}

main.c

@@ -7,8 +7,6 @@
 #include "mmpriv.h"
 #include "ketopt.h"
 
-#define MM_VERSION "2.24-r1122"
-
 #ifdef __linux__
 #include <sys/resource.h>
 #include <sys/time.h>
@@ -78,6 +76,10 @@ static ko_longopt_t long_options[] = {
 	{ "chain-skip-scale",ko_required_argument,351 },
 	{ "print-chains",   ko_no_argument,       352 },
 	{ "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' },
@@ -121,7 +123,7 @@ static inline void yes_or_no(mm_mapopt_t *opt, int64_t flag, int long_idx, const
 
 int main(int argc, char *argv[])
 {
-	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:yYPo:e:U:";
+	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:b:O:E:m:N:Qu:R:hF:LC:yYPo:e:U:J:";
 	ketopt_t o = KETOPT_INIT;
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
@@ -179,6 +181,7 @@ int main(int argc, char *argv[])
 		else if (c == 'm') opt.min_chain_score = atoi(o.arg);
 		else if (c == 'A') opt.a = atoi(o.arg);
 		else if (c == 'B') opt.b = atoi(o.arg);
+		else if (c == 'b') opt.transition = atoi(o.arg);
 		else if (c == 's') opt.min_dp_max = atoi(o.arg);
 		else if (c == 'C') opt.noncan = atoi(o.arg);
 		else if (c == 'I') ipt.batch_size = mm_parse_num(o.arg);
@@ -187,7 +190,12 @@ int main(int argc, char *argv[])
 		else if (c == 'R') rg = o.arg;
 		else if (c == 'h') fp_help = stdout;
 		else if (c == '2') opt.flag |= MM_F_2_IO_THREADS;
-		else if (c == 'o') {
+		else if (c == 'J') {
+			int t;
+			t = atoi(o.arg);
+			if (t == 0) opt.flag |= MM_F_SPLICE_OLD;
+			else if (t == 1) opt.flag &= ~MM_F_SPLICE_OLD;
+		} else if (c == 'o') {
 			if (strcmp(o.arg, "-") != 0) {
 				if (freopen(o.arg, "wb", stdout) == NULL) {
 					fprintf(stderr, "[ERROR]\033[1;31m failed to write the output to file '%s'\033[0m: %s\n", o.arg, strerror(errno));
@@ -237,6 +245,10 @@ int main(int argc, char *argv[])
 		else if (c == 350) opt.q_occ_frac = atof(o.arg); // --q-occ-frac
 		else if (c == 352) mm_dbg_flag |= MM_DBG_PRINT_CHAIN; // --print-chains
 		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
@@ -261,7 +273,8 @@ int main(int argc, char *argv[])
 		} else if (c == 326) { // --dual
 			yes_or_no(&opt, MM_F_NO_DUAL, o.longidx, o.arg, 0);
 		} else if (c == 347) { // --rmq
-			yes_or_no(&opt, MM_F_RMQ, o.longidx, o.arg, 1);
+			if (o.arg) yes_or_no(&opt, MM_F_RMQ, o.longidx, o.arg, 1);
+			else opt.flag |= MM_F_RMQ;
 		} else if (c == 'S') {
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR | MM_F_OUT_CS_LONG;
 			if (mm_verbose >= 2)
@@ -322,7 +335,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -H           use homopolymer-compressed k-mer (preferrable for PacBio)\n");
 		fprintf(fp_help, "    -k INT       k-mer size (no larger than 28) [%d]\n", ipt.k);
 		fprintf(fp_help, "    -w INT       minimizer window size [%d]\n", ipt.w);
-		fprintf(fp_help, "    -I NUM       split index for every ~NUM input bases [4G]\n");
+		fprintf(fp_help, "    -I NUM       split index for every ~NUM input bases [8G]\n");
 		fprintf(fp_help, "    -d FILE      dump index to FILE []\n");
 		fprintf(fp_help, "  Mapping:\n");
 		fprintf(fp_help, "    -f FLOAT     filter out top FLOAT fraction of repetitive minimizers [%g]\n", opt.mid_occ_frac);
@@ -344,6 +357,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -z INT[,INT] Z-drop score and inversion Z-drop score [%d,%d]\n", opt.zdrop, opt.zdrop_inv);
 		fprintf(fp_help, "    -s INT       minimal peak DP alignment score [%d]\n", opt.min_dp_max);
 		fprintf(fp_help, "    -u CHAR      how to find GT-AG. f:transcript strand, b:both strands, n:don't match GT-AG [n]\n");
+		fprintf(fp_help, "    -J INT       splice mode. 0: original minimap2 model; 1: miniprot model [1]\n");
 		fprintf(fp_help, "  Input/Output:\n");
 		fprintf(fp_help, "    -a           output in the SAM format (PAF by default)\n");
 		fprintf(fp_help, "    -o FILE      output alignments to FILE [stdout]\n");
@@ -351,6 +365,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -R STR       SAM read group line in a format like '@RG\\tID:foo\\tSM:bar' []\n");
 		fprintf(fp_help, "    -c           output CIGAR in PAF\n");
 		fprintf(fp_help, "    --cs[=STR]   output the cs tag; STR is 'short' (if absent) or 'long' [none]\n");
+		fprintf(fp_help, "    --ds         output the ds tag, which is an extension to cs\n");
 		fprintf(fp_help, "    --MD         output the MD tag\n");
 		fprintf(fp_help, "    --eqx        write =/X CIGAR operators\n");
 		fprintf(fp_help, "    -Y           use soft clipping for supplementary alignments\n");
@@ -360,12 +375,12 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    --version    show version number\n");
 		fprintf(fp_help, "  Preset:\n");
 		fprintf(fp_help, "    -x STR       preset (always applied before other options; see minimap2.1 for details) []\n");
-		fprintf(fp_help, "                 - map-pb/map-ont - PacBio CLR/Nanopore vs reference mapping\n");
+		fprintf(fp_help, "                 - lr:hq - accurate long reads (error rate <1%%) against a reference genome\n");
-		fprintf(fp_help, "                 - map-hifi - PacBio HiFi reads vs reference mapping\n");
+		fprintf(fp_help, "                 - splice/splice:hq - spliced alignment for long reads/accurate long reads\n");
-		fprintf(fp_help, "                 - ava-pb/ava-ont - PacBio/Nanopore read overlap\n");
 		fprintf(fp_help, "                 - asm5/asm10/asm20 - asm-to-ref mapping, for ~0.1/1/5%% sequence divergence\n");
-		fprintf(fp_help, "                 - splice/splice:hq - long-read/Pacbio-CCS spliced alignment\n");
+		fprintf(fp_help, "                 - sr - short reads against a reference\n");
-		fprintf(fp_help, "                 - sr - genomic short-read mapping\n");
+		fprintf(fp_help, "                 - map-pb/map-hifi/map-ont/map-iclr - CLR/HiFi/Nanopore/ICLR vs reference mapping\n");
+		fprintf(fp_help, "                 - ava-pb/ava-ont - PacBio CLR/Nanopore read overlap\n");
 		fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of these and other advanced command-line options.\n");
 		return fp_help == stdout? 0 : 1;
 	}

map.c

@@ -10,11 +10,6 @@
 #include "bseq.h"
 #include "khash.h"
 
-struct mm_tbuf_s {
-	void *km;
-	int rep_len, frag_gap;
-};
-
 mm_tbuf_t *mm_tbuf_init(void)
 {
 	mm_tbuf_t *b;

minimap.h

@@ -5,41 +5,46 @@
 #include <stdio.h>
 #include <sys/types.h>
 
-#define MM_F_NO_DIAG       0x001 // no exact diagonal hit
+#define MM_VERSION "2.28-r1209"
-#define MM_F_NO_DUAL       0x002 // skip pairs where query name is lexicographically larger than target name
+
-#define MM_F_CIGAR         0x004
+#define MM_F_NO_DIAG       (0x001LL) // no exact diagonal hit
-#define MM_F_OUT_SAM       0x008
+#define MM_F_NO_DUAL       (0x002LL) // skip pairs where query name is lexicographically larger than target name
-#define MM_F_NO_QUAL       0x010
+#define MM_F_CIGAR         (0x004LL)
-#define MM_F_OUT_CG        0x020
+#define MM_F_OUT_SAM       (0x008LL)
-#define MM_F_OUT_CS        0x040
+#define MM_F_NO_QUAL       (0x010LL)
-#define MM_F_SPLICE        0x080 // splice mode
+#define MM_F_OUT_CG        (0x020LL)
-#define MM_F_SPLICE_FOR    0x100 // match GT-AG
+#define MM_F_OUT_CS        (0x040LL)
-#define MM_F_SPLICE_REV    0x200 // match CT-AC, the reverse complement of GT-AG
+#define MM_F_SPLICE        (0x080LL) // splice mode
-#define MM_F_NO_LJOIN      0x400
+#define MM_F_SPLICE_FOR    (0x100LL) // match GT-AG
-#define MM_F_OUT_CS_LONG   0x800
+#define MM_F_SPLICE_REV    (0x200LL) // match CT-AC, the reverse complement of GT-AG
-#define MM_F_SR            0x1000
+#define MM_F_NO_LJOIN      (0x400LL)
-#define MM_F_FRAG_MODE     0x2000
+#define MM_F_OUT_CS_LONG   (0x800LL)
-#define MM_F_NO_PRINT_2ND  0x4000
+#define MM_F_SR            (0x1000LL)
-#define MM_F_2_IO_THREADS  0x8000
+#define MM_F_FRAG_MODE     (0x2000LL)
-#define MM_F_LONG_CIGAR    0x10000
+#define MM_F_NO_PRINT_2ND  (0x4000LL)
-#define MM_F_INDEPEND_SEG  0x20000
+#define MM_F_2_IO_THREADS  (0x8000LL)
-#define MM_F_SPLICE_FLANK  0x40000
+#define MM_F_LONG_CIGAR    (0x10000LL)
-#define MM_F_SOFTCLIP      0x80000
+#define MM_F_INDEPEND_SEG  (0x20000LL)
-#define MM_F_FOR_ONLY      0x100000
+#define MM_F_SPLICE_FLANK  (0x40000LL)
-#define MM_F_REV_ONLY      0x200000
+#define MM_F_SOFTCLIP      (0x80000LL)
-#define MM_F_HEAP_SORT     0x400000
+#define MM_F_FOR_ONLY      (0x100000LL)
-#define MM_F_ALL_CHAINS    0x800000
+#define MM_F_REV_ONLY      (0x200000LL)
-#define MM_F_OUT_MD        0x1000000
+#define MM_F_HEAP_SORT     (0x400000LL)
-#define MM_F_COPY_COMMENT  0x2000000
+#define MM_F_ALL_CHAINS    (0x800000LL)
-#define MM_F_EQX           0x4000000 // use =/X instead of M
+#define MM_F_OUT_MD        (0x1000000LL)
-#define MM_F_PAF_NO_HIT    0x8000000 // output unmapped reads to PAF
+#define MM_F_COPY_COMMENT  (0x2000000LL)
-#define MM_F_NO_END_FLT    0x10000000
+#define MM_F_EQX           (0x4000000LL) // use =/X instead of M
-#define MM_F_HARD_MLEVEL   0x20000000
+#define MM_F_PAF_NO_HIT    (0x8000000LL) // output unmapped reads to PAF
-#define MM_F_SAM_HIT_ONLY  0x40000000
+#define MM_F_NO_END_FLT    (0x10000000LL)
+#define MM_F_HARD_MLEVEL   (0x20000000LL)
+#define MM_F_SAM_HIT_ONLY  (0x40000000LL)
 #define MM_F_RMQ           (0x80000000LL)
 #define MM_F_QSTRAND       (0x100000000LL)
 #define MM_F_NO_INV        (0x200000000LL)
 #define MM_F_NO_HASH_NAME  (0x400000000LL)
+#define MM_F_SPLICE_OLD    (0x800000000LL)
+#define MM_F_SECONDARY_SEQ (0x1000000000LL)	//output SEQ field for seqondary alignments using hard clipping
+#define MM_F_OUT_DS        (0x2000000000LL)
 
 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2
@@ -93,6 +98,7 @@ typedef struct {
 typedef struct {
 	uint32_t capacity;                  // the capacity of cigar[]
 	int32_t dp_score, dp_max, dp_max2;  // DP score; score of the max-scoring segment; score of the best alternate mappings
+	int32_t dp_max0;                    // DP score before mm_update_dp_max() adjustment
 	uint32_t n_ambi:30, trans_strand:2; // number of ambiguous bases; transcript strand: 0 for unknown, 1 for +, 2 for -
 	uint32_t n_cigar;                   // number of cigar operations in cigar[]
 	uint32_t cigar[];
@@ -149,6 +155,7 @@ typedef struct {
 	float alt_drop;
 
 	int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties
+	int transition; // transition mismatch score (A:G, C:T)
 	int sc_ambi; // score when one or both bases are "N"
 	int noncan;      // cost of non-canonical splicing sites
 	int junc_bonus;
@@ -189,6 +196,11 @@ typedef struct {
 } mm_idx_reader_t;
 
 // memory buffer for thread-local storage during mapping
+struct mm_tbuf_s {
+	void *km;
+	int rep_len, frag_gap;
+};
+
 typedef struct mm_tbuf_s mm_tbuf_t;
 
 // global variables

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "18 December 2021" "minimap2-2.24 (r1122)" "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
@@ -79,6 +79,19 @@ Minimizer k-mer length [15]
 .BI -w \ INT
 Minimizer window size [10]. A minimizer is the smallest k-mer
 in a window of w consecutive k-mers.
+.TP
+.BI -j \ INT
+Syncmer submer size [10]. Option
+.B -j
+and
+.B -w
+will override each: if
+.B -w
+is applied after
+.BR -j ,
+.B -j
+will have no effect, and vice versa.
+
 .TP
 .B -H
 Use homopolymer-compressed (HPC) minimizers. An HPC sequence is constructed by
@@ -88,16 +101,17 @@ on the HPC sequence.
 .BI -I \ NUM
 Load at most
 .I NUM
-target bases into RAM for indexing [4G]. If there are more than
+target bases into RAM for indexing [8G]. If there are more than
 .I NUM
 bases in
 .IR target.fa ,
 minimap2 needs to read
 .I query.fa
-multiple times to map it against each batch of target sequences.
+multiple times to map it against each batch of target sequences. This would create a multi-part index.
 .I NUM
 may be ending with k/K/m/M/g/G. NB: mapping quality is incorrect given a
-multi-part index.
+multi-part index. See also option
+.BR --split-prefix .
 .TP
 .B --idx-no-seq
 Don't store target sequences in the index. It saves disk space and memory but
@@ -254,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
@@ -329,6 +348,10 @@ Matching score [2]
 .BI -B \ INT
 Mismatching penalty [4]
 .TP
+.BI -b \ INT
+Mismatching penalty for transitions [same as
+.BR -B ].
+.TP
 .BI -O \ INT1[,INT2]
 Gap open penalty [4,24]. If
 .I INT2
@@ -342,10 +365,19 @@ costs
 .RI min{ O1 + k * E1 , O2 + k * E2 }.
 In the splice mode, the second gap penalties are not used.
 .TP
+.BI -J \ INT
+Splice model [1]. 0 for the original minimap2 splice model that always penalizes non-GT-AG splicing;
+1 for the miniprot model that considers non-GT-AG. Option
+.B -C
+has no effect with the default
+.BR -J1 .
+.BR -J0 .
+.TP
 .BI -C \ INT
 Cost for a non-canonical GT-AG splicing (effective with
-.BR --splice )
+.B --splice
-[0]
+.BR -J0 )
+[0].
 .TP
 .BI -z \ INT1[,INT2]
 Truncate an alignment if the running alignment score drops too quickly along
@@ -436,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
@@ -492,6 +524,9 @@ Output =/X CIGAR operators for sequence match/mismatch.
 .B -Y
 In SAM output, use soft clipping for supplementary alignments.
 .TP
+.B --secondary-seq
+In SAM output, show query sequences for secondary alignments.
+.TP
 .BI --seed \ INT
 Integer seed for randomizing equally best hits. Minimap2 hashes
 .I INT
@@ -552,15 +587,43 @@ are:
 Align noisy long reads of ~10% error rate to a reference genome. This is the
 default mode.
 .TP
+.B lr:hq
+Align accurate long reads (error rate <1%) to a reference genome
+.RB ( -k19
+.B -w19 -U50,500
+.BR -g10k ).
+This was recommended by ONT developers for recent Nanopore reads
+produced with chemistry v14 that can reach ~99% in accuracy.
+It was shown to work better for accurate Nanopore reads
+than
+.BR map-hifi .
+.TP
 .B map-hifi
 Align PacBio high-fidelity (HiFi) reads to a reference genome
-.RB ( -k19
+.RB ( -xlr:hq
-.B -w19 -U50,500 -g10k -A1 -B4 -O6,26 -E2,1
+.B -A1 -B4 -O6,26 -E2,1
 .BR -s200 ).
+It differs from
+.B lr:hq
+only in scoring. It has not been tested whether
+.B lr:hq
+would work better for PacBio HiFi reads.
 .TP
 .B map-pb
 Align older PacBio continuous long (CLR) reads to a reference genome
 .RB ( -Hk19 ).
+Note that this data type is effectively deprecated by HiFi.
+Unless you work on very old data, you probably want to use
+.B map-hifi
+or
+.BR lr:hq .
+.TP
+.B map-iclr
+Align Illumina Complete Long Reads (ICLR) to a reference genome
+.RB ( -k19
+.B -B6 -b4
+.BR -O10,50 ).
+This was recommended by Illumina developers.
 .TP
 .B asm5
 Long assembly to reference mapping
@@ -568,26 +631,26 @@ Long assembly to reference mapping
 .B -w19 -U50,500 --rmq -r1k,100k -g10k -A1 -B19 -O39,81 -E3,1 -s200 -z200
 .BR -N50 ).
 Typically, the alignment will not extend to regions with 5% or higher sequence
-divergence. Only use this preset if the average divergence is far below 5%.
+divergence. Use this preset if the average divergence is not much higher than 0.1%.
 .TP
 .B asm10
 Long assembly to reference mapping
 .RB ( -k19
 .B -w19 -U50,500 --rmq -r1k,100k -g10k -A1 -B9 -O16,41 -E2,1 -s200 -z200
 .BR -N50 ).
-Up to 10% sequence divergence.
+Use this if the average divergence is around 1%.
 .TP
 .B asm20
 Long assembly to reference mapping
 .RB ( -k19
 .B -w10 -U50,500 --rmq -r1k,100k -g10k -A1 -B4 -O6,26 -E2,1 -s200 -z200
 .BR -N50 ).
-Up to 20% sequence divergence.
+Use this if the average divergence is around several percent.
 .TP
 .B splice
 Long-read spliced alignment
 .RB ( -k15
-.B -w5 --splice -g2k -G200k -A1 -B2 -O2,32 -E1,0 -b0 -C9 -z200 -ub --junc-bonus=9 --cap-sw-mem=0
+.B -w5 --splice -g2k -G200k -A1 -B2 -O2,32 -E1,0 -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
@@ -598,13 +661,13 @@ costs are different during chaining; 4) the computation of the
 tag ignores introns to demote hits to pseudogenes.
 .TP
 .B splice:hq
-Long-read splice alignment for PacBio CCS reads
+Spliced alignment for accurate long RNA-seq reads such as PacBio iso-seq
 .RB ( -xsplice
 .B -C5 -O6,24
 .BR -B4 ).
 .TP
 .B sr
-Short single-end reads without splicing
+Short-read alignment without splicing
 .RB ( -k21
 .B -w11 --sr --frag=yes -A2 -B8 -O12,32 -E2,1 -b0 -r100 -p.5 -N20 -f1000,5000 -n2 -m25
 .B -s40 -g100 -2K50m --heap-sort=yes

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = '2.24-r1122';
+var paftools_version = '2.28-r1209';
 
 /*****************************
  ***** Library functions *****
@@ -133,26 +133,50 @@ Interval.find_ovlp = function(a, st, en)
 
 function fasta_read(fn)
 {
-	var h = {}, gt = '>'.charCodeAt(0);
+	var h = {}, seqlen = [];
+	var buf = new Bytes();
 	var file = fn == '-'? new File() : new File(fn);
-	var buf = new Bytes(), seq = null, name = null, seqlen = [];
+	if (typeof k8_version == "undefined") { // for k8-0.x
-	while (file.readline(buf) >= 0) {
+		var seq = null, name = null, gt = '>'.charCodeAt(0);
-		if (buf[0] == gt) {
+		while (file.readline(buf) >= 0) {
-			if (seq != null && name != null) {
+			if (buf[0] == gt) {
-				seqlen.push([name, seq.length]);
+				if (seq != null && name != null) {
-				h[name] = seq;
+					seqlen.push([name, seq.length]);
-				name = seq = null;
+					h[name] = seq;
-			}
+					name = seq = null;
-			var m, line = buf.toString();
+				}
-			if ((m = /^>(\S+)/.exec(line)) != null) {
+				var m, line = buf.toString();
-				name = m[1];
+				if ((m = /^>(\S+)/.exec(line)) != null) {
-				seq = new Bytes();
+					name = m[1];
-			}
+					seq = new Bytes();
-		} else seq.set(buf);
+				}
-	}
+			} else seq.set(buf);
-	if (seq != null && name != null) {
+		}
-		seqlen.push([name, seq.length]);
+		if (seq != null && name != null) {
-		h[name] = seq;
+			seqlen.push([name, seq.length]);
+			h[name] = seq;
+		}
+	} else { // for k8-1.x
+		var seq = null, name = null;
+		while (file.readline(buf) >= 0) {
+			var line = buf.toString();
+			if (line[0] == ">") {
+				if (seq != null && name != null) {
+					seqlen.push([name, seq.length]);
+					h[name] = new Uint8Array(seq.buffer);
+					name = seq = null;
+				}
+				var m;
+				if ((m = /^>(\S+)/.exec(line)) != null) {
+					name = m[1];
+					seq = new Bytes();
+				}
+			} else seq.set(line);
+		}
+		if (seq != null && name != null) {
+			seqlen.push([name, seq.length]);
+			h[name] = new Uint8Array(seq.buffer);
+		}
 	}
 	buf.destroy();
 	file.close();
@@ -161,16 +185,27 @@ function fasta_read(fn)
 
 function fasta_free(fa)
 {
-	for (var name in fa)
+	if (typeof k8_version == "undefined")
-		fa[name].destroy();
+		for (var name in fa)
+			fa[name].destroy();
+	// FIXME: for k8-1.0, sequences are not freed. This is ok for now but not general.
 }
 
 Bytes.prototype.reverse = function()
 {
-	for (var i = 0; i < this.length>>1; ++i) {
+	if (typeof k8_version === "undefined") { // k8-0.x
-		var tmp = this[i];
+		for (var i = 0; i < this.length>>1; ++i) {
-		this[i] = this[this.length - i - 1];
+			var tmp = this[i];
-		this[this.length - i - 1] = tmp;
+			this[i] = this[this.length - i - 1];
+			this[this.length - i - 1] = tmp;
+		}
+	} else { // k8-1.x
+		var buf = new Uint8Array(this.buffer);
+		for (var i = 0; i < buf.length>>1; ++i) {
+			var tmp = buf[i];
+			buf[i] = buf[buf.length - i - 1];
+			buf[buf.length - i - 1] = tmp;
+		}
 	}
 }
 
@@ -185,13 +220,24 @@ Bytes.prototype.revcomp = function()
 		for (var i = 0; i < s1.length; ++i)
 			Bytes.rctab[s1.charCodeAt(i)] = s2.charCodeAt(i);
 	}
-	for (var i = 0; i < this.length>>1; ++i) {
+	if (typeof k8_version === "undefined") { // k8-0.x
-		var tmp = this[this.length - i - 1];
+		for (var i = 0; i < this.length>>1; ++i) {
-		this[this.length - i - 1] = Bytes.rctab[this[i]];
+			var tmp = this[this.length - i - 1];
-		this[i] = Bytes.rctab[tmp];
+			this[this.length - i - 1] = Bytes.rctab[this[i]];
+			this[i] = Bytes.rctab[tmp];
+		}
+		if (this.length&1)
+			this[this.length>>1] = Bytes.rctab[this[this.length>>1]];
+	} else { // k8-1.x
+		var buf = new Uint8Array(this.buffer);
+		for (var i = 0; i < buf.length>>1; ++i) {
+			var tmp = buf[buf.length - i - 1];
+			buf[buf.length - i - 1] = Bytes.rctab[buf[i]];
+			buf[i] = Bytes.rctab[tmp];
+		}
+		if (buf.length&1)
+			buf[buf.length>>1] = Bytes.rctab[buf[buf.length>>1]];
 	}
-	if (this.length&1)
-		this[this.length>>1] = Bytes.rctab[this[this.length>>1]];
 }
 
 /********************
@@ -1532,22 +1578,24 @@ function paf_view(args)
 
 function paf_gff2bed(args)
 {
-	var c, fn_ucsc_fai = null, is_short = false, keep_gff = false, print_junc = false, output_gene = false;
+	var c, fn_ucsc_fai = null, is_short = false, keep_gff = false, print_junc = false, output_gene = false, ens_canon_only = false;
-	while ((c = getopt(args, "u:sgjG")) != null) {
+	while ((c = getopt(args, "u:sgjGe")) != 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;
+		else if (c == 'e') ens_canon_only = true;
 	}
 
 	if (getopt.ind == args.length) {
 		print("Usage: paftools.js gff2bed [options] <in.gff>");
 		print("Options:");
-		print("  -j       Output junction BED");
+		print("  -j       output junction BED");
-		print("  -s       Print names in the short form");
+		print("  -s       print names in the short form");
 		print("  -u FILE  hg38.fa.fai for chr name conversion");
-		print("  -g       Output GFF (used with -u)");
+		print("  -e       only show transcript tagged with 'Ensembl_canonical'");
+		print("  -g       output GFF (used with -u)");
 		exit(1);
 	}
 
@@ -1606,7 +1654,7 @@ 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|tag) "([^"]+)";/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;
@@ -1646,13 +1694,14 @@ function paf_gff2bed(args)
 		if (t[2] != "CDS" && t[2] != "exon") continue;
 		t[3] = parseInt(t[3]) - 1;
 		t[4] = parseInt(t[4]);
-		var id = null, type = "", name = "N/A", biotype = "", m, tname = "N/A";
+		var id = null, type = "", name = "N/A", biotype = "", m, tname = "N/A", ens_canonical = false;
 		while ((m = re_gtf.exec(t[8])) != null) {
 			if (m[1] == "transcript_id") id = m[2];
 			else if (m[1] == "transcript_type") type = m[2];
 			else if (m[1] == "transcript_biotype" || m[1] == "gbkey") biotype = m[2];
 			else if (m[1] == "gene_name" || m[1] == "gene_id") name = m[2];
 			else if (m[1] == "transcript_name") tname = m[2];
+			else if (m[1] == "tag" && m[2] == "Ensembl_canonical") ens_canonical = true;
 		}
 		while ((m = re_gff3.exec(t[8])) != null) {
 			if (m[1] == "transcript_id") id = m[2];
@@ -1661,6 +1710,7 @@ function paf_gff2bed(args)
 			else if (m[1] == "gene_name" || m[1] == "gene_id") name = m[2];
 			else if (m[1] == "transcript_name") tname = m[2];
 		}
+		if (ens_canon_only && !ens_canonical) continue;
 		if (type == "" && biotype != "") type = biotype;
 		if (id == null) throw Error("No transcript_id");
 		if (id != last_id) {
@@ -1690,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);
 	}
@@ -1725,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
@@ -1837,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(""));
@@ -2047,7 +2103,7 @@ function paf_mapeval(args)
 		warn("Usage: paftools.js mapeval [options] <in.paf>|<in.sam>");
 		warn("Options:");
 		warn("  -r FLOAT   mapping correct if overlap_length/union_length>FLOAT [" + ovlp_ratio + "]");
-		warn("  -Q INT     print wrong mappings with mapQ>INT [don't print]");
+		warn("  -Q INT     print wrong mappings with mapQ>=INT [don't print]");
 		warn("  -m INT     0: eval the longest aln only; 1: first aln only; 2: all primary aln [0]");
 		exit(1);
 	}
@@ -2341,12 +2397,15 @@ function paf_pbsim2fq(args)
 
 function paf_junceval(args)
 {
-	var c, l_fuzzy = 0, print_ovlp = false, print_err_only = false, first_only = false, chr_only = false;
+	var c, l_fuzzy = 0, print_ovlp = false, print_err_only = false, first_only = false, chr_only = false, aa = false, is_bed = false;
-	while ((c = getopt(args, "l:epc")) != null) {
+	while ((c = getopt(args, "l:epcab1")) != null) {
 		if (c == 'l') l_fuzzy = parseInt(getopt.arg);
 		else if (c == 'e') print_err_only = print_ovlp = true;
 		else if (c == 'p') print_ovlp = true;
 		else if (c == 'c') chr_only = true;
+		else if (c == 'a') aa = true;
+		else if (c == 'b') is_bed = true;
+		else if (c == '1') first_only = true;
 	}
 
 	if (args.length - getopt.ind < 1) {
@@ -2356,6 +2415,9 @@ function paf_junceval(args)
 		print("  -p        print overlapping introns");
 		print("  -e        print erroreous overlapping introns");
 		print("  -c        only consider alignments to /^(chr)?([0-9]+|X|Y)$/");
+		print("  -a        miniprot PAF as input");
+		print("  -b        BED as input");
+		print("  -1        only process the first alignment of each query");
 		exit(1);
 	}
 
@@ -2409,13 +2471,17 @@ function paf_junceval(args)
 
 	file = getopt.ind+1 >= args.length || args[getopt.ind+1] == '-'? new File() : new File(args[getopt.ind+1]);
 	var last_qname = null;
-	var re_cigar = /(\d+)([MIDNSHP=X])/g;
+	var re_cigar = /(\d+)([MIDNSHP=XFGUV])/g;
 	while (file.readline(buf) >= 0) {
 		var m, t = buf.toString().split("\t");
-		var ctg_name = null, cigar = null, pos = null, qname = t[0];
+		var ctg_name = null, cigar = null, pos = null, qname;
 
 		if (t[0].charAt(0) == '@') continue;
-		if (t[4] == '+' || t[4] == '-' || t[4] == '*') { // PAF
+		if (t[0] == "##PAF") t.shift();
+		qname = t[0];
+		if (is_bed) {
+			ctg_name = t[0], pos = parseInt(t[1]), cigar == null;
+		} else if (t[4] == '+' || t[4] == '-' || t[4] == '*') { // PAF
 			ctg_name = t[5], pos = parseInt(t[7]);
 			var type = 'P';
 			for (i = 12; i < t.length; ++i) {
@@ -2445,12 +2511,43 @@ function paf_junceval(args)
 		}
 
 		var intron = [];
-		while ((m = re_cigar.exec(cigar)) != null) {
+		if (is_bed) {
-			var len = parseInt(m[1]), op = m[2];
+			intron.push([pos, parseInt(t[2])]);
-			if (op == 'N') {
+		} else if (aa) {
-				intron.push([pos, pos + len]);
+			var tmp_junc = [], tmp = 0;
-				pos += len;
+			while ((m = re_cigar.exec(cigar)) != null) {
-			} else if (op == 'M' || op == 'X' || op == '=' || op == 'D') pos += len;
+				var len = parseInt(m[1]), op = m[2];
+				if (op == 'N') {
+					tmp_junc.push([tmp, tmp + len]);
+					tmp += len;
+				} else if (op == 'U') {
+					tmp_junc.push([tmp + 1, tmp + len - 2]);
+					tmp += len;
+				} else if (op == 'V') {
+					tmp_junc.push([tmp + 2, tmp + len - 1]);
+					tmp += len;
+				} else if (op == 'M' || op == 'X' || op == '=' || op == 'D') {
+					tmp += len * 3;
+				} else if (op == 'F' || op == 'G') {
+					tmp += len;
+				}
+			}
+			if (t[4] == '+') {
+				for (var i = 0; i < tmp_junc.length; ++i)
+					intron.push([pos + tmp_junc[i][0], pos + tmp_junc[i][1]]);
+			} else if (t[4] == '-') {
+				var glen = parseInt(t[8]) - parseInt(t[7]);
+				for (var i = tmp_junc.length - 1; i >= 0; --i)
+					intron.push([pos + (glen - tmp_junc[i][1]), pos + (glen - tmp_junc[i][0])]);
+			}
+		} else {
+			while ((m = re_cigar.exec(cigar)) != null) {
+				var len = parseInt(m[1]), op = m[2];
+				if (op == 'N') {
+					intron.push([pos, pos + len]);
+					pos += len;
+				} else if (op == 'M' || op == 'X' || op == '=' || op == 'D') pos += len;
+			}
 		}
 		if (intron.length == 0) {
 			++n_sgl;
@@ -2509,6 +2606,276 @@ function paf_junceval(args)
 	}
 }
 
+function paf_exoneval(args) // adapted from paf_junceval()
+{
+	var c, l_fuzzy = 0, print_ovlp = false, print_err_only = false, first_only = false, chr_only = false, aa = false, is_bed = false, use_cds = false, eval_base = false;
+	while ((c = getopt(args, "l:epcab1ds")) != null) {
+		if (c == 'l') l_fuzzy = parseInt(getopt.arg);
+		else if (c == 'e') print_err_only = print_ovlp = true;
+		else if (c == 'p') print_ovlp = true;
+		else if (c == 'c') chr_only = true;
+		else if (c == 'a') aa = true, use_cds = true;
+		else if (c == 'b') is_bed = true;
+		else if (c == '1') first_only = true;
+		else if (c == 'd') use_cds = true;
+		else if (c == 's') eval_base = true;
+	}
+
+	if (args.length - getopt.ind < 1) {
+		print("Usage: paftools.js exoneval [options] <gene.gtf> <aln.sam>");
+		print("Options:");
+		print("  -l INT    tolerance of junction positions (0 for exact) [0]");
+		print("  -d        evaluate coding regions only (exon regions by default)");
+		print("  -a        miniprot PAF as input (force -d)");
+		print("  -p        print overlapping exons");
+		print("  -e        print erroreous overlapping exons");
+		print("  -c        only consider alignments to /^(chr)?([0-9]+|X|Y)$/");
+		print("  -1        only process the first alignment of each query");
+		print("  -b        BED as input");
+		print("  -s        compute base Sn and Sp (more memory)");
+		exit(1);
+	}
+
+	var file, buf = new Bytes();
+
+	warn("Reading reference GTF...");
+	var tr = {};
+	file = args[getopt.ind] == '-'? new File() : new File(args[getopt.ind]);
+	while (file.readline(buf) >= 0) {
+		var m, t = buf.toString().split("\t");
+		if (t[0].charAt(0) == '#') continue;
+		if (use_cds) {
+			if (t[2] != "cds" && t[2] != "CDS") continue;
+		} else {
+			if (t[2] != 'exon') continue;
+		}
+		var st = parseInt(t[3]) - 1;
+		var en = parseInt(t[4]);
+		if ((m = /transcript_id "(\S+)"/.exec(t[8])) == null) continue;
+		var tid = m[1];
+		if (tr[tid] == null) tr[tid] = [t[0], t[6], 0, 0, []];
+		tr[tid][4].push([st, en]); // this keeps transcript
+	}
+	file.close();
+
+	var anno = {};
+	for (var tid in tr) { // traverse each transcript
+		var t = tr[tid];
+		Interval.sort(t[4]);
+		t[2] = t[4][0][0];
+		t[3] = t[4][t[4].length - 1][1];
+		if (anno[t[0]] == null) anno[t[0]] = [];
+		var s = t[4];
+		for (var i = 0; i < s.length; ++i) // traverse each exon
+			anno[t[0]].push([s[i][0], s[i][1]]);
+	}
+	tr = null;
+
+	for (var chr in anno) { // index exons
+		var e = anno[chr];
+		if (e.length == 0) continue;
+		Interval.sort(e);
+		var k = 0;
+		for (var i = 1; i < e.length; ++i) // dedup
+			if (e[i][0] != e[k][0] || e[i][1] != e[k][1])
+				e[++k] = e[i].slice(0);
+		e.length = k + 1;
+		Interval.index_end(e);
+	}
+
+	var n_pri = 0, n_unmapped = 0, n_mapped = 0;
+	var n_exon = 0, n_exon_hit = 0, n_exon_novel = 0;
+
+	file = getopt.ind+1 >= args.length || args[getopt.ind+1] == '-'? new File() : new File(args[getopt.ind+1]);
+	var last_qname = null, qexon = {};
+	var re_cigar = /(\d+)([MIDNSHP=XFGUV])/g;
+
+	warn("Evaluating alignments...");
+	while (file.readline(buf) >= 0) {
+		var m, t = buf.toString().split("\t");
+		var ctg_name = null, cigar = null, pos = null, qname;
+
+		if (t[0].charAt(0) == '@') continue;
+		if (t[0] == "##PAF") t.shift();
+		qname = t[0];
+		if (is_bed) {
+			ctg_name = t[0], pos = parseInt(t[1]), cigar == null;
+		} else if (t[4] == '+' || t[4] == '-' || t[4] == '*') { // PAF
+			ctg_name = t[5], pos = parseInt(t[7]);
+			var type = 'P';
+			for (i = 12; i < t.length; ++i) {
+				if ((m = /^(tp:A|cg:Z):(\S+)/.exec(t[i])) != null) {
+					if (m[1] == 'tp:A') type = m[2];
+					else cigar = m[2];
+				}
+			}
+			if (type == 'S') continue; // secondary
+		} else { // SAM
+			ctg_name = t[2], pos = parseInt(t[3]) - 1, cigar = t[5];
+			var flag = parseInt(t[1]);
+			if (flag&0x100) continue; // secondary
+		}
+
+		if (chr_only && !/^(chr)?([0-9]+|X|Y)$/.test(ctg_name)) continue;
+		if (first_only && last_qname == qname) continue;
+		if (ctg_name == '*') { // unmapped
+			++n_unmapped;
+			continue;
+		} else {
+			++n_pri;
+			if (last_qname != qname) {
+				++n_mapped;
+				last_qname = qname;
+			}
+		}
+
+		var exon = [];
+		if (is_bed) { // BED
+			exon.push([pos, parseInt(t[2])]);
+		} else if (aa) {
+			var tmp_exon = [], tmp = 0, tmp_st = 0;
+			while ((m = re_cigar.exec(cigar)) != null) {
+				var len = parseInt(m[1]), op = m[2];
+				if (op == 'N') {
+					tmp_exon.push([tmp_st, tmp]);
+					tmp_st = tmp + len, tmp += len;
+				} else if (op == 'U') {
+					tmp_exon.push([tmp_st, tmp + 1]);
+					tmp_st = tmp + len - 2, tmp += len;
+				} else if (op == 'V') {
+					tmp_exon.push([tmp_st, tmp + 2]);
+					tmp_st = tmp + len - 1, tmp += len;
+				} else if (op == 'M' || op == 'X' || op == '=' || op == 'D') {
+					tmp += len * 3;
+				} else if (op == 'F' || op == 'G') {
+					tmp += len;
+				}
+			}
+			tmp_exon.push([tmp_st, tmp]);
+			if (t[4] == '+') {
+				for (var i = 0; i < tmp_exon.length; ++i)
+					exon.push([pos + tmp_exon[i][0], pos + tmp_exon[i][1]]);
+			} else if (t[4] == '-') { // For protein-to-genome alignment, the coordinates are on the query strand. Need to flip them.
+				var glen = parseInt(t[8]) - parseInt(t[7]);
+				for (var i = tmp_exon.length - 1; i >= 0; --i)
+					exon.push([pos + (glen - tmp_exon[i][1]), pos + (glen - tmp_exon[i][0])]);
+			}
+		} else {
+			var tmp_st = pos;
+			while ((m = re_cigar.exec(cigar)) != null) {
+				var len = parseInt(m[1]), op = m[2];
+				if (op == 'N') {
+					exon.push([tmp_st, pos]);
+					tmp_st = pos + len, pos += len;
+				} else if (op == 'M' || op == 'X' || op == '=' || op == 'D') pos += len;
+			}
+			exon.push([tmp_st, pos]);
+		}
+		n_exon += exon.length;
+
+		var chr = anno[ctg_name];
+		if (chr != null) {
+			for (var i = 0; i < exon.length; ++i) {
+				if (eval_base) {
+					if (qexon[ctg_name] == null) qexon[ctg_name] = [];
+					qexon[ctg_name].push([exon[i][0], exon[i][1]]);
+				}
+				var o = Interval.find_ovlp(chr, exon[i][0], exon[i][1]);
+				if (o.length > 0) {
+					var hit = false;
+					for (var j = 0; j < o.length; ++j) {
+						var st_diff = exon[i][0] - o[j][0];
+						var en_diff = exon[i][1] - o[j][1];
+						if (st_diff < 0) st_diff = -st_diff;
+						if (en_diff < 0) en_diff = -en_diff;
+						if (st_diff <= l_fuzzy && en_diff <= l_fuzzy)
+							++n_exon_hit, hit = true;
+						if (hit) break;
+					}
+					if (print_ovlp) {
+						var type = hit? 'C' : 'P';
+						if (hit && print_err_only) continue;
+						var x = '[';
+						for (var j = 0; j < o.length; ++j) {
+							if (j) x += ', ';
+							x += '(' + o[j][0] + "," + o[j][1] + ')';
+						}
+						x += ']';
+						print(type, qname, i+1, ctg_name, exon[i][0], exon[i][1], x);
+					}
+				} else {
+					++n_exon_novel;
+					if (print_ovlp)
+						print('N',  qname, i+1, ctg_name, exon[i][0], exon[i][1]);
+				}
+			}
+		} else {
+			n_exon_novel += exon.length;
+		}
+	}
+	file.close();
+
+	buf.destroy();
+
+	if (!print_ovlp) {
+		print("# unmapped reads: " + n_unmapped);
+		print("# mapped reads: " + n_mapped);
+		print("# primary alignments: " + n_pri);
+		print("# predicted exons: " + n_exon);
+		print("# non-overlapping exons: " + n_exon_novel);
+		print("# correct exons: " + n_exon_hit + " (" + (n_exon_hit / n_exon * 100).toFixed(2) + "%)");
+	}
+
+	function merge_and_index(ex) {
+		for (var chr in ex) {
+			var a = [];
+			e = ex[chr];
+			Interval.sort(e);
+			var st = e[0][0], en = e[0][1];
+			for (var i = 1; i < e.length; ++i) { // merge
+				if (e[i][0] > en) {
+					a.push([st, en]);
+					st = e[i][0], en = e[i][1];
+				} else {
+					en = en > e[i][1]? en : e[i][1];
+				}
+			}
+			a.push([st, en]);
+			Interval.index_end(a);
+			ex[chr] = a;
+		}
+	}
+
+	function cal_sn(a0, a1) {
+		var tot = 0, cov = 0;
+		for (var chr in a1) {
+			var e0 = a0[chr], e1 = a1[chr];
+			for (var i = 0; i < e1.length; ++i)
+				tot += e1[i][1] - e1[i][0];
+			if (e0 == null) continue;
+			for (var i = 0; i < e1.length; ++i) {
+				var o = Interval.find_ovlp(e0, e1[i][0], e1[i][1]);
+				for (var j = 0; j < o.length; ++j) { // this only works when there are no overlaps between intervals
+					var st = e1[i][0] > o[j][0]? e1[i][0] : o[j][0];
+					var en = e1[i][1] < o[j][1]? e1[i][1] : o[j][1];
+					cov += en - st;
+				}
+			}
+		}
+		return [tot, cov];
+	}
+
+	if (eval_base) {
+		warn("Computing base Sn and Sp...");
+		merge_and_index(qexon);
+		merge_and_index(anno);
+		var sn = cal_sn(qexon, anno);
+		var sp = cal_sn(anno, qexon);
+		print("Base Sn: " + sn[1] + " / " + sn[0] + " = " + (sn[1] / sn[0] * 100).toFixed(2) + "%");
+		print("Base Sp: " + sp[1] + " / " + sp[0] + " = " + (sp[1] / sp[0] * 100).toFixed(2) + "%");
+	}
+}
+
 // evaluate overlap sensitivity
 function paf_ov_eval(args)
 {
@@ -2704,6 +3071,23 @@ function paf_misjoin(args)
 		return len < (en - st) * cen_ratio? false : true;
 	}
 
+	function test_cen_point(cen, chr, x) {
+		var b = cen[chr];
+		if (b == null) return false;
+		for (var j = 0; j < b.length; ++j)
+			if (x >= b[j][0] && x < b[j][1])
+				return true;
+		return false;
+	}
+
+	if (show_err || show_long) {
+		print("C\tJ  inter-chromosomal misjoin");
+		print("C\tj  inter-chromosomal misjoin with both breakpoints ending in centromeres");
+		print("C\tG  long gap on the reference genome");
+		print("C\tg  long gap on the reference genome with both breakpoints ending in centromeres");
+		print("C\tM  closed inversion");
+		print("C");
+	}
 	function process(a) {
 		var k = 0;
 		for (var i = 0; i < a.length; ++i) {
@@ -2716,14 +3100,17 @@ function paf_misjoin(args)
 		a = a.sort(function(x,y){return x[2]-y[2]});
 		if (show_long) for (var i = 0; i < a.length; ++i) print(a[i].join("\t"));
 		for (var i = 1; i < a.length; ++i) {
-			var ov = [false, false];
+			var ov = [false, false], end_cen = [false, false];
 			ov[0] = test_cen(cen, a[i-1][5], a[i-1][7], a[i-1][8]);
 			ov[1] = test_cen(cen, a[i][5], a[i][7], a[i][8]);
+			end_cen[0] = test_cen_point(cen, a[i-1][5], a[i-1][4] == '+'? a[i-1][8] : a[i-1][7]);
+			end_cen[1] = test_cen_point(cen, a[i][5],   a[i][4] == '+'?   a[i][7]   : a[i][8]);
 			if (a[i-1][5] != a[i][5]) { // different chr
 				if (ov[0] || ov[1]) ++n_diff[1];
 				else if (show_err) {
-					print("J", a[i-1].slice(0, 12).join("\t"));
+					var label = end_cen[0] && end_cen[1]? 'j' : 'J';
-					print("J", a[i].slice(0, 12).join("\t"));
+					print(label, a[i-1].slice(0, 12).join("\t"));
+					print(label, a[i].slice(0, 12).join("\t"));
 				}
 				++n_diff[0];
 			} else if (a[i-1][4] == a[i][4]) { // a gap
@@ -2733,8 +3120,9 @@ function paf_misjoin(args)
 				if (gap > max_gap) {
 					if (ov[0] || ov[1]) ++n_gap[1];
 					else if (show_err) {
-						print("G", a[i-1].slice(0, 12).join("\t"));
+						var label = end_cen[0] && end_cen[1]? 'g' : 'G';
-						print("G", a[i].slice(0, 12).join("\t"));
+						print(label, a[i-1].slice(0, 12).join("\t"));
+						print(label, a[i].slice(0, 12).join("\t"));
 					}
 					++n_gap[0];
 				}
@@ -3084,6 +3472,183 @@ function paf_pafcmp(args)
 	buf.destroy();
 }
 
+function paf_longcs2seq(args) {
+	var c, opt = { query:false };
+	while ((c = getopt(args, "q")) != null)
+		if (c == 'q') opt.query = true;
+	if (args.length == getopt.ind) {
+		print("Usage: paftools.js longcs2seq [-q] <long-cs.paf>");
+		return;
+	}
+	var re_cs = /([:=*+-])(\d+|[A-Za-z]+)/g
+	var buf = new Bytes();
+	var file = args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]);
+	while (file.readline(buf) >= 0) {
+		var m, cs = null, t = buf.toString().split("\t");
+		for (var i = 12; i < t.length; ++i)
+			if ((m = /^cs:Z:(\S+)/.exec(t[i])) != null) {
+				cs = m[1];
+				break;
+			}
+		if (cs == null) continue;
+		var ts = "", qs = "";
+		while ((m = re_cs.exec(cs)) != null) {
+			if (m[1] == "=") ts += m[2], qs += m[2];
+			else if (m[1] == "+") qs += m[2].toUpperCase();
+			else if (m[1] == "-") ts += m[2].toUpperCase();
+			else if (m[1] == "*") ts += m[2][0].toUpperCase(), qs += m[2][1].toUpperCase();
+			else if (m[1] == ":") throw Error("Long cs is required");
+		}
+		if (opt.query) {
+			print(">" + t[0] + "_" + t[2] + "_" + t[3]);
+			print(qs);
+		} else {
+			print(">" + t[5] + "_" + t[7] + "_" + t[8]);
+			print(ts);
+		}
+	}
+	file.close();
+	buf.destroy();
+}
+
+function paf_paf2gff(args) {
+	var c, opt = { aa:false };
+	var re_cigar = /(\d+)([A-Z=])/g;
+	while ((c = getopt(args, "a")) != null) {
+		if (c == 'a') opt.aa = true;
+	}
+	if (args.length == getopt.ind) {
+		print("Usage: paftools.js paf2gff [-a] <in.paf>");
+		return;
+	}
+	var buf = new Bytes();
+	var file = args[getopt.ind] == '-'? new File() : new File(args[getopt.ind]);
+	var hid = 1, last_name = null;
+	while (file.readline(buf) >= 0) {
+		var m, t = buf.toString().split("\t");
+		if (t[5] == '*') continue; // skip unmapped lines
+
+		if (t[0] != last_name) last_name = t[0], hid = 1;
+		else ++hid;
+		for (var i = 1; i <= 3; ++i) t[i] = parseInt(t[i]);
+		for (var i = 6; i <= 11; ++i) t[i] = parseInt(t[i]);
+		var cigar = null, score = null, np = null, dist_stop = null, dist_start = null;
+		for (var i = 12; i < t.length; ++i) {
+			if ((m = /^(cg:Z|AS:i|np:i|da:i|do:i):(\S+)/.exec(t[i])) != null) {
+				if (m[1] == 'cg:Z') cigar = m[2];
+				else if (m[1] == 'AS:i') score = parseInt(m[2]);
+				else if (m[1] == 'np:i') np = parseInt(m[2]);
+				else if (m[1] == 'do:i') dist_stop = parseInt(m[2]);
+				else if (m[1] == 'da:i') dist_start = parseInt(m[2]);
+			}
+		}
+		if (cigar == null) throw Error("failed to find the cg:Z tag");
+		if (score == null) throw Error("failed to find the AS:i tag");
+
+		var st = 0, en = 0, phase = 0, pseudo = false, fs = 0, a = [];
+		if (dist_start != null && dist_start == 0)
+			a.push([t[5], 'paf2gff', 'start_codon', 0, 3, 0, t[4], '.', 0]);
+		while ((m = re_cigar.exec(cigar)) != null) {
+			var len = parseInt(m[1]);
+			if (m[2] == 'M' || m[2] == 'D') {
+				en += opt.aa? len * 3 : len;
+			} else if (m[2] == 'F' || m[2] == 'G' || m[2] == 'R') {
+				en += len, pseudo = true, fs = 1;
+			} else if (m[2] == 'N') {
+				a.push([t[5], 'paf2gff', 'exon', st, en, 0, t[4], phase, fs]);
+				st = en + len, en += len, phase = 0, fs = 0;
+			} else if (m[2] == 'U') { // ...xGT...AGxx...
+				a.push([t[5], 'paf2gff', 'exon', st, en + 1, 0, t[4], phase, fs]);
+				st = en + len - 2, en += len, phase = 2, fs = 0;
+			} else if (m[2] == 'V') { // ...xxGT...AGx...
+				a.push([t[5], 'paf2gff', 'exon', st, en + 2, 0, t[4], phase, fs]);
+				st = en + len - 1, en += len, phase = 1, fs = 0;
+			}
+		}
+		a.push([t[5], 'paf2gff', 'exon', st, en, 0, t[4], phase, fs]);
+		if (en != t[8] - t[7]) throw Error("inconsistent cigar");
+		if (dist_stop != null && dist_stop == 0)
+			a.push([t[5], 'paf2gff', 'stop_codon', en, en + 3, 0, t[4], '.', 0]);
+		var type = pseudo? 'pseudogene' : 'protein_coding';
+		var attr = ['transcript_id=' + t[0] + '#' + hid, 'transcript_type=' + type].join(";");
+		var trans_attr = 'identity=' + (t[9] / t[10]).toFixed(4);
+		if (np != null) trans_attr += ';positive=' + (np * 3 / t[10]).toFixed(4);
+		trans_attr += ';aa_start=' + t[2];
+		trans_attr += ';aa_end=' + (t[1] - t[3]);
+		if (dist_start != null && dist_start >= 0) trans_attr += ';dist_start_codon=' + dist_start;
+		if (dist_stop != null && dist_stop >= 0) trans_attr += ';dist_stop_codon=' + dist_stop;
+		var trans_st = t[7], trans_en = t[8];
+		if (dist_stop != null && dist_stop == 0) {
+			if (t[4] == '-') trans_st -= 3;
+			else trans_en += 3;
+		}
+		print([t[5], 'paf2gff', 'transcript', trans_st + 1, trans_en, score, t[4], '.', attr + ';' + trans_attr].join("\t"));
+		if (opt.aa && t[4] == '-') {
+			var b = [], len = t[8] - t[7];
+			for (var i = a.length - 1; i >= 0; --i) {
+				var x = len - a[i][3];
+				a[i][3] = len - a[i][4];
+				a[i][4] = x;
+				//a[i][7] = a[i][7] == 0? 0 : 3 - a[i][7]; // not sure if this line is needed
+				b.push(a[i]);
+			}
+			a = b;
+		}
+		for (var i = 0; i < a.length; ++i) {
+			if (!pseudo && a[i][2] == "exon") a[i][2] = "CDS";
+			a[i][3] += t[7] + 1;
+			a[i][4] += t[7];
+			a[i][8] = attr + ";frameshift=" + a[i][8];
+			print(a[i].join("\t"));
+		}
+	}
+	file.close();
+	buf.destroy();
+}
+
+function paf_gff2junc(args) {
+	var c, feat = "CDS";
+	while ((c = getopt(args, "f:")) != null) {
+		if (c == 'f') feat = getopt.arg;
+	}
+	if (getopt.ind == args.length) {
+		print("Usage: paftools.js gff2junc [-f feature] <in.gff3>");
+		return;
+	}
+	var buf = new Bytes();
+	var file = args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]);
+
+	function process_a(a) {
+		if (a.length < 2) return;
+		a = a.sort(function(x, y) { return x[4] - y[4] });
+		for (var i = 1; i < a.length; ++i)
+			print([a[i][1], a[i-1][5], a[i][4], a[i][0], 0, a[i][7]].join("\t"));
+	}
+
+	var a = [];
+	while (file.readline(buf) >= 0) {
+		var m, t = buf.toString().split("\t");
+		if (t[0][0] == '#') continue;
+		if (t[2].toLowerCase() != feat.toLowerCase()) continue;
+		//print(t.join("\t"));
+		if ((m = /\bParent=([^;]+)/.exec(t[8])) == null) {
+			warn("Can't find Parent");
+			continue;
+		}
+		t[3] = parseInt(t[3]) - 1;
+		t[4] = parseInt(t[4]);
+		t.unshift(m[1]);
+		if (a.length > 0 && a[0][0] != m[1]) {
+			process_a(a);
+			a.length = 0;
+			a.push(t);
+		} else a.push(t);
+	}
+	process_a(a);
+	file.close();
+	buf.destroy();
+}
+
 /*************************
  ***** main function *****
  *************************/
@@ -3098,6 +3663,9 @@ function main(args)
 		print("  sam2paf    convert SAM to PAF");
 		print("  delta2paf  convert MUMmer's delta to PAF");
 		print("  gff2bed    convert GTF/GFF3 to BED12");
+		print("  gff2junc   convert GFF3 to junction BED");
+		print("  longcs2seq convert long-cs PAF to sequences");
+//		print("  paf2gff    convert PAF to GFF3 (tested for miniprot only)");
 		print("");
 		print("  stat       collect basic mapping information in PAF/SAM");
 		print("  asmstat    collect basic assembly information");
@@ -3115,6 +3683,7 @@ function main(args)
 		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");
+		print("  exoneval   evaluate exon-level consistency with known annotations");
 		print("  ov-eval    evaluate read overlap sensitivity using read-to-ref mapping");
 		exit(1);
 	}
@@ -3125,6 +3694,7 @@ function main(args)
 	else if (cmd == 'delta2paf') paf_delta2paf(args);
 	else if (cmd == 'splice2bed') paf_splice2bed(args);
 	else if (cmd == 'gff2bed') paf_gff2bed(args);
+	else if (cmd == 'gff2junc') paf_gff2junc(args);
 	else if (cmd == 'stat') paf_stat(args);
 	else if (cmd == 'asmstat') paf_asmstat(args);
 	else if (cmd == 'asmgene') paf_asmgene(args);
@@ -3138,10 +3708,13 @@ function main(args)
 	else if (cmd == 'mason2fq') paf_mason2fq(args);
 	else if (cmd == 'pbsim2fq') paf_pbsim2fq(args);
 	else if (cmd == 'junceval') paf_junceval(args);
+	else if (cmd == 'exoneval') paf_exoneval(args);
 	else if (cmd == 'ov-eval') paf_ov_eval(args);
 	else if (cmd == 'vcfstat') paf_vcfstat(args);
 	else if (cmd == 'sveval') paf_sveval(args);
 	else if (cmd == 'vcfsel') paf_vcfsel(args);
+	else if (cmd == 'longcs2seq') paf_longcs2seq(args);
+	else if (cmd == 'paf2gff') paf_paf2gff(args);
 	else if (cmd == 'version') print(paftools_version);
 	else throw Error("unrecognized command: " + cmd);
 }

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

@@ -8,7 +8,7 @@ void mm_idxopt_init(mm_idxopt_t *opt)
 	opt->k = 15, opt->w = 10, opt->flag = 0;
 	opt->bucket_bits = 14;
 	opt->mini_batch_size = 50000000;
-	opt->batch_size = 4000000000ULL;
+	opt->batch_size = 8000000000ULL;
 }
 
 void mm_mapopt_init(mm_mapopt_t *opt)
@@ -45,6 +45,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->alt_drop = 0.15f;
 
 	opt->a = 2, opt->b = 4, opt->q = 4, opt->e = 2, opt->q2 = 24, opt->e2 = 1;
+	opt->transition = 0;
 	opt->sc_ambi = 1;
 	opt->zdrop = 400, opt->zdrop_inv = 200;
 	opt->end_bonus = -1;
@@ -54,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;
@@ -90,7 +91,7 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 	if (preset == 0) {
 		mm_idxopt_init(io);
 		mm_mapopt_init(mo);
-	} else if (strcmp(preset, "map-ont") == 0) { // this is the same as the default
+	} else if (strcmp(preset, "lr") == 0 || strcmp(preset, "map-ont") == 0) { // this is the same as the default
 	} else if (strcmp(preset, "ava-ont") == 0) {
 		io->flag = 0, io->k = 15, io->w = 5;
 		mo->flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN;
@@ -105,13 +106,30 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_chain_skip = 25;
 		mo->bw_long = mo->bw;
 		mo->occ_dist = 0;
-	} else if (strcmp(preset, "map-hifi") == 0 || strcmp(preset, "map-ccs") == 0) {
+	} else if (strcmp(preset, "lr:hq") == 0 || strcmp(preset, "map-hifi") == 0 || strcmp(preset, "map-ccs") == 0) {
 		io->flag = 0, io->k = 19, io->w = 19;
 		mo->max_gap = 10000;
-		mo->a = 1, mo->b = 4, mo->q = 6, mo->q2 = 26, mo->e = 2, mo->e2 = 1;
-		mo->occ_dist = 500;
 		mo->min_mid_occ = 50, mo->max_mid_occ = 500;
-		mo->min_dp_max = 200;
+		if (strcmp(preset, "map-hifi") == 0 || strcmp(preset, "map-ccs") == 0) {
+			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;
+		mo->q = 10, mo->q2 = 50;
+	} else if (strcmp(preset, "map-iclr") == 0) {
+		io->flag = 0, io->k = 19;
+		mo->b = 6, mo->transition = 4;
+		mo->q = 10, mo->q2 = 50;
 	} else if (strncmp(preset, "asm", 3) == 0) {
 		io->flag = 0, io->k = 19, io->w = 19;
 		mo->bw = 1000, mo->bw_long = 100000;
@@ -156,7 +174,7 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		mo->junc_bonus = 9;
 		mo->zdrop = 200, mo->zdrop_inv = 100; // because mo->a is halved
 		if (strcmp(preset, "splice:hq") == 0)
-			mo->junc_bonus = 5, mo->b = 4, mo->q = 6, mo->q2 = 24;
+			mo->noncan = 5, mo->b = 4, mo->q = 6, mo->q2 = 24;
 	} else return -1;
 	return 0;
 }

pyproject.toml

@@ -0,0 +1,2 @@
+[build-system]
+requires = ["setuptools", "wheel", "Cython"]

python/README.rst

@@ -77,7 +77,9 @@ This constructor accepts the following arguments:
 
 * **min_chain_score**: minimum chaing score
 
-* **bw**: chaining and alignment band width
+* **bw**: chaining and alignment band width (initial chaining and extension)
+
+* **bw_long**: chaining and alignment band width (RMQ-based rechaining and closing gaps)
 
 * **best_n**: max number of alignments to return
 

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.24'
+__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:
@@ -112,7 +113,7 @@ cdef class Aligner:
 	cdef cmappy.mm_idxopt_t idx_opt
 	cdef cmappy.mm_mapopt_t map_opt
 
-	def __cinit__(self, fn_idx_in=None, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, best_n=None, n_threads=3, fn_idx_out=None, max_frag_len=None, extra_flags=None, seq=None, scoring=None):
+	def __cinit__(self, fn_idx_in=None, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, bw_long=None, best_n=None, n_threads=3, fn_idx_out=None, max_frag_len=None, extra_flags=None, seq=None, scoring=None):
 		self._idx = NULL
 		cmappy.mm_set_opt(NULL, &self.idx_opt, &self.map_opt) # set the default options
 		if preset is not None:
@@ -125,6 +126,7 @@ cdef class Aligner:
 		if min_chain_score is not None: self.map_opt.min_chain_score = min_chain_score
 		if min_dp_score is not None: self.map_opt.min_dp_max = min_dp_score
 		if bw is not None: self.map_opt.bw = bw
+		if bw_long is not None: self.map_opt.bw_long = bw_long
 		if best_n is not None: self.map_opt.best_n = best_n
 		if max_frag_len is not None: self.map_opt.max_frag_len = max_frag_len
 		if extra_flags is not None: self.map_opt.flag |= extra_flags
@@ -172,6 +174,7 @@ cdef class Aligner:
 		cdef cmappy.mm_mapopt_t map_opt
 
 		if self._idx == NULL: return
+		if ((self.map_opt.flag & 4) and (self._idx.flag & 2)): return
 		map_opt = self.map_opt
 		if max_frag_len is not None: map_opt.max_frag_len = max_frag_len
 		if extra_flags is not None: map_opt.flag |= extra_flags
@@ -217,6 +220,7 @@ cdef class Aligner:
 		cdef int l
 		cdef char *s
 		if self._idx == NULL: return
+		if ((self.map_opt.flag & 4) and (self._idx.flag & 2)): return
 		s = cmappy.mappy_fetch_seq(self._idx, name.encode(), start, end, &l)
 		if l == 0: return None
 		r = s[:l] if isinstance(s, str) else s[:l].decode()

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

@@ -7,7 +7,7 @@ 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);
+	a = Kmalloc(km, mm128_t, 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);
@@ -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.24',
+	version = '2.28',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),