Update README.md

This is more apparent when there are many candidate chains. Although only a
small numbers of them are extended, they are still occupying memory. A
realloc() solves this problem. This is a long existiing issue.

.gitmodules

@@ -1,3 +1,6 @@
 [submodule "lib/simde"]
 	path = lib/simde
 	url = https://github.com/nemequ/simde.git
+[submodule "ext/TAL"]
+	path = ext/TAL
+	url = https://github.com/IntelLabs/Trans-Omics-Acceleration-Library.git

Makefile

@@ -1,6 +1,43 @@
 CFLAGS=		-g -Wall -O2 -Wc++-compat #-Wextra
-CPPFLAGS=	-DHAVE_KALLOC
-INCLUDES=
+CPPFLAGS=	-DHAVE_KALLOC #-march=native #-DALIGN_AVX -DPARALLEL_CHAINING #-DMANUAL_PROFILING
+COMP_FLAG = -march=native
+
+ifeq ($(avx2_compile), 1)
+	COMP_FLAG = -mavx2
+endif
+
+#CPPFLAGS=	-DHAVE_KALLOC -mavx2 -DALIGN_AVX -DAPPLY_AVX2 -DPARALLEL_CHAINING #-DLISA_HASH -DUINT64 -DVECTORIZE #-DMANUAL_PROFILING
+#CPPFLAGS=	-DHAVE_KALLOC -mavx2 -DPARALLEL_CHAINING  #-DMANUAL_PROFILING
+
+OPT_FLAGS= -DPARALLEL_CHAINING -DALIGN_AVX -DAPPLY_AVX2 
+OPT_FLAGS+=$(COMP_FLAG)
+ifeq ($(lhash_index), 1)	
+	CPPFLAGS+=	-DLISA_INDEX 
+endif
+ifeq ($(lhash), 1)	
+	OPT_FLAGS+=	-DLISA_HASH -DUINT64 -DVECTORIZE 
+endif
+ifeq ($(manual_profile), 1)
+	CPPFLAGS+= -DMANUAL_PROFILING 
+endif
+
+#ifeq ($(use_avx2), 1)
+#	OPT_FLAGS+= -DAPPLY_AVX2
+#endif
+
+ifeq ($(disable_output), 1)
+	CPPFLAGS+= -DDISABLE_OUTPUT
+endif
+
+ifeq ($(no_opt),)
+	CPPFLAGS+= $(OPT_FLAGS)
+endif
+
+
+
+#INCLUDES=
+#INCLUDES=	-I./ext/TAL_offline/src/LISA-hash #-I./ext/TAL/src/dynamic-programming 
+INCLUDES=	-I./ext/TAL/src/LISA-hash -I./ext/TAL/src/dynamic-programming 
 OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o \
 			lchain.o align.o hit.o seed.o map.o format.o pe.o esterr.o splitidx.o \
 			ksw2_ll_sse.o
@@ -8,9 +45,14 @@ PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite
 LIBS=		-lm -lz -lpthread
 
+CC=$(CXX)
+ifeq ($(CC), g++)
+	CC=g++ -std=c++11
+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
+	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 ksw2_extd2_avx.o
 else                # if sse2only is defined
 	OBJS+=ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_exts2_sse.o
 endif
@@ -56,6 +98,17 @@ libminimap2.a:$(OBJS)
 sdust:sdust.c kalloc.o kalloc.h kdq.h kvec.h kseq.h ketopt.h sdust.h
 		$(CC) -D_SDUST_MAIN $(CFLAGS) $< kalloc.o -o $@ -lz
 
+multi:
+		$(MAKE) clean
+		$(MAKE)
+		mv minimap2 mm2-fast
+		$(MAKE) clean
+		$(MAKE) lhash=1
+		mv minimap2 mm2-fast-lhash
+		$(MAKE) clean
+		$(MAKE) no_opt=1
+		mv minimap2 mm2-fast-no-opt
+
 # SSE-specific targets on x86/x86_64
 
 ifeq ($(arm_neon),)   # if arm_neon is defined, compile this target with the default setting (i.e. no -msse2)

Makefile.simde

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

NEWS.md

@@ -1,3 +1,39 @@
+Release 2.22-r1101 (7 August 2021)
+----------------------------------
+
+When choosing the best alignment, this release uses logarithm gap penalty and
+query-specific mismatch penalty. It improves the sensitivity to long INDELs in
+repetitive regions.
+
+Other notable changes:
+
+ * Bugfix: fixed an indirect memory leak that may waste a large amount of
+   memory given highly repetitive reference such as a 16S RNA database (#749).
+   All versions of minimap2 have this issue.
+
+ * New feature: added --cap-kalloc to reduce the peak memory. This option is
+   not enabled by default but may become the default in future releases.
+
+Known issue:
+
+ * Minimap2 may take a long time to map a read (#771). So far it is not clear
+   if this happens to v2.18 and earlier versions.
+
+(2.22: 7 August 2021, r1101)
+
+
+
+Release 2.21-r1071 (6 July 2021)
+--------------------------------
+
+This release fixed a regression in short-read mapping introduced in v2.19
+(#776). It also fixed invalid comparisons of uninitialized variables, though
+these are harmless (#752). Long-read alignment should be identical to v2.20.
+
+(2.21: 6 July 2021, r1071)
+
+
+
 Release 2.20-r1061 (27 May 2021)
 --------------------------------
 

README.md

@@ -1,3 +1,76 @@
+## mm2-fast
+### Introduction
+mm2-fast is an accelerated implementation of minimap2 on modern CPUs. mm2-fast accelerates all the three major modules of minimap2: (a) seeding, (b) chaining, and (c) pairwise alignment, achieving up to 1.8x speedup using AVX512 over minimap2. 
+mm2-fast is a drop-in replacement of minimap2, providing the same functionality with the exact same output.
+In the current version, all the modules are optimized using **AVX-512** and **AVX2** vectorization. Detailed benchmark results are available in our [preprint](https://doi.org/10.1101/2021.07.21.453294).
+
+### System requirement
+Operating System: Linux   
+mm2-fast was tested using g++ (GCC) 9.2.0 and icpc version 19.1.3.304   
+Architecture: x86\_64 CPUs with [AVX512, AVX2](https://en.wikipedia.org/wiki/Advanced_Vector_Extensions)   
+Memory requirement: ~30GB for human genome   
+
+### Installation
+Clone the *fast-contrib-v2.22* branch from minimap2 github page. The source code can be compiled by using *make* command. It only takes a few seconds.
+```
+git clone --recursive https://github.com/lh3/minimap2.git -b fast-contrib-v2.22 mm2-fast   
+cd mm2-fast
+make 
+```
+
+### Usage
+The usage of mm2-fast is same as minimap2. Here is an example of mapping ONT reads with test data.
+```sh
+./minimap2 -ax map-ont test/MT-human.fa test/MT-orang.fa > mm2-fast_output
+```
+
+### Accuracy evaluation
+As mm2-fast is an accelerated version of minimap2-v2.22, the output of mm2-fast can be verified against minimap2-v2.22. Note that the optimized chaining in mm2-fast is strictly required to be run with a chaining parameter *max-chain-skip=infinity*. Note that having parameter *max-chain-skip=infinity* leads to higher chaining precision. Therefore, for correctness verification, minimap2 should run with a larger value of *max-chain-skip* parameter. Follow the below steps to verify the accuracy of mm2-fast.
+```sh
+git clone --recursive https://github.com/lh3/minimap2.git -b fast-contrib-v2.22 mm2-fast   
+cd mm2-fast && make 
+./minimap2 -ax map-ont test/MT-human.fa test/MT-orang.fa --max-chain-skip=1000000 > mm2-fast_output
+```
+```sh
+git clone https://github.com/lh3/minimap2.git -b v2.22
+cd minimap2 && make
+./minimap2 -ax map-ont test/MT-human.fa test/MT-orang.fa --max-chain-skip=1000000 > minimap2_output
+```
+The output generated by minimap2 and mm2-fast should match. 
+```sh
+diff minimap2_output mm2-fast_output > diff_result
+```
+The file ```diff_result``` should be empty, meaning a difference of 0 lines.
+
+### Advanced options
+The default compilation using make applies two optimizations: vectorized chaining and sequence alignment. The learned-indexes based seeding is disabled by default as it requires availability of [Rust](https://en.wikipedia.org/wiki/Rust_(programming_language)). This is because the learned hash-table uses an external training library that runs on Rust. Rust is trivial to install, see https://rustup.rs/ and add its path to .bashrc file. Rust installation only takes a few seconds. Following are the steps to enable learned hash table optimization in mm2-fast:
+```sh
+# Start by building learned hash table index for optimized seeding module 
+./build_rmi.sh test/MT-human.fa map-ont               ##Takes two arguments: 1. path-to-reference-seq-file 2. preset. 
+						      ##For human genome, this step should take around 2-3 minutes to finish.    
+
+# Next, compile and run the mapping phase  
+make clean && make lhash=1
+./minimap2 -ax map-ont test/MT-human.fa test/MT-orang.fa > mm2-fast-lhash_output
+```
+To compile mm2-fast with all optimizations turned off and switch back to default minimap2, use the following command during compilation. This could be useful for debugging.
+```sh
+make clean && make no_opt=1
+```
+
+### Performance
+We have observed up to 1.8x speedup across datasets (please refer to the paper for more details). For example, for the randomly sampled 100K reads from ["HG002\_GM24385\_1\_2\_3\_Guppy\_3.6.0\_prom.fastq.gz"](https://precision.fda.gov/challenges/10/view), minimap2 takes 92 seconds, while mm2-fast takes 54 seconds to map against the human genome on a 28 cores Intel® Xeon® Platinum 8280 CPUs.  Our sampled datasets with 100K reads are available [here](https://drive.google.com/drive/folders/1131j7ejHdT7QZnjxLcTLi5qqwYcfFbuv).
+
+### Future Plans
+
+
+### Citations
+["Accelerating long-read analysis on modern CPUs"](https://doi.org/10.1101/2021.07.21.453294); Saurabh Kalikar, Chirag Jain, Vasimuddin Md, Sanchit Misra; BioRxiv 2021
+
+---
+The original README content of minimap2 follows.
+
+
 [![GitHub Downloads](https://img.shields.io/github/downloads/lh3/minimap2/total.svg?style=social&logo=github&label=Download)](https://github.com/lh3/minimap2/releases)
 [![BioConda Install](https://img.shields.io/conda/dn/bioconda/minimap2.svg?style=flag&label=BioConda%20install)](https://anaconda.org/bioconda/minimap2)
 [![PyPI](https://img.shields.io/pypi/v/mappy.svg?style=flat)](https://pypi.python.org/pypi/mappy)
@@ -14,7 +87,7 @@ cd minimap2 && make
 # use presets (no test data)
 ./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 lateer)
+./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 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)
@@ -74,8 +147,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.20/minimap2-2.20_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.20_x64-linux/minimap2
+curl -L https://github.com/lh3/minimap2/releases/download/v2.22/minimap2-2.22_x64-linux.tar.bz2 | tar -jxvf -
+./minimap2-2.22_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
@@ -96,7 +169,7 @@ with the ARM related command lines given above.
 
 Without any options, minimap2 takes a reference database and a query sequence
 file as input and produce approximate mapping, without base-level alignment
-(i.e. no CIGAR), in the [PAF format][paf]:
+(i.e. coordinates are only approximate and no CIGAR in output), in the [PAF format][paf]:
 ```sh
 minimap2 ref.fa query.fq > approx-mapping.paf
 ```

align.c

@@ -5,6 +5,13 @@
 #include "minimap.h"
 #include "mmpriv.h"
 #include "ksw2.h"
+#include "ksw2_extd2_avx.h"
+#include <x86intrin.h>
+extern uint64_t avg;
+extern uint64_t alignment_time;
+extern void *km1;
+extern uint64_t km_size;// = 500000000; // 500 MB
+extern int km_top;
 
 static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t sc_ambi)
 {
@@ -53,16 +60,16 @@ static int mm_test_zdrop(void *km, const mm_mapopt_t *opt, const uint8_t *qseq,
 	// find the score and the region where score drops most along diagonal
 	for (k = 0, score = 0; k < n_cigar; ++k) {
 		uint32_t l, op = cigar[k]&0xf, len = cigar[k]>>4;
-		if (op == 0) {
+		if (op == MM_CIGAR_MATCH) {
 			for (l = 0; l < len; ++l) {
 				score += mat[tseq[i + l] * 5 + qseq[j + l]];
 				update_max_zdrop(score, i+l, j+l, &max, &max_i, &max_j, opt->e, &max_zdrop, pos);
 			}
 			i += len, j += len;
-		} else if (op == 1 || op == 2 || op == 3) {
+		} else if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL || op == MM_CIGAR_N_SKIP) {
 			score -= opt->q + opt->e * len;
-			if (op == 1) j += len; // insertion
-			else i += len;         // deletion
+			if (op == MM_CIGAR_INS) j += len;
+			else i += len;
 			update_max_zdrop(score, i, j, &max, &max_i, &max_j, opt->e, &max_zdrop, pos);
 		}
 	}
@@ -98,12 +105,12 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 	for (k = 0; k < p->n_cigar; ++k) { // indel left alignment
 		uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4;
 		if (len == 0) to_shrink = 1;
-		if (op == 0) {
+		if (op == MM_CIGAR_MATCH) {
 			toff += len, qoff += len;
-		} else if (op == 1 || op == 2) { // insertion or deletion
+		} else if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL) {
 			if (k > 0 && k < p->n_cigar - 1 && (p->cigar[k-1]&0xf) == 0 && (p->cigar[k+1]&0xf) == 0) {
 				int l, prev_len = p->cigar[k-1] >> 4;
-				if (op == 1) {
+				if (op == MM_CIGAR_INS) {
 					for (l = 0; l < prev_len; ++l)
 						if (qseq[qoff - 1 - l] != qseq[qoff + len - 1 - l])
 							break;
@@ -116,9 +123,9 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 					p->cigar[k-1] -= l<<4, p->cigar[k+1] += l<<4, qoff -= l, toff -= l;
 				if (l == prev_len) to_shrink = 1;
 			}
-			if (op == 1) qoff += len;
+			if (op == MM_CIGAR_INS) qoff += len;
 			else toff += len;
-		} else if (op == 3) {
+		} else if (op == MM_CIGAR_N_SKIP) {
 			toff += len;
 		}
 	}
@@ -128,13 +135,13 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 			uint32_t l, s[3] = {0,0,0};
 			for (l = k; l < p->n_cigar; ++l) { // count number of adjacent I and D
 				uint32_t op = p->cigar[l]&0xf;
-				if (op == 1 || op == 2 || p->cigar[l]>>4 == 0)
+				if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL || p->cigar[l]>>4 == 0)
 					s[op] += p->cigar[l] >> 4;
 				else break;
 			}
 			if (s[1] > 0 && s[2] > 0 && l - k > 2) { // turn to a single I and a single D
-				p->cigar[k]   = s[1]<<4|1;
-				p->cigar[k+1] = s[2]<<4|2;
+				p->cigar[k]   = s[1]<<4|MM_CIGAR_INS;
+				p->cigar[k+1] = s[2]<<4|MM_CIGAR_DEL;
 				for (k += 2; k < l; ++k)
 					p->cigar[k] &= 0xf;
 				to_shrink = 1;
@@ -154,9 +161,9 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 			else p->cigar[k+1] += p->cigar[k]>>4<<4; // add length to the next CIGAR operator
 		p->n_cigar = l;
 	}
-	if ((p->cigar[0]&0xf) == 1 || (p->cigar[0]&0xf) == 2) { // get rid of leading I or D
+	if ((p->cigar[0]&0xf) == MM_CIGAR_INS || (p->cigar[0]&0xf) == MM_CIGAR_DEL) { // get rid of leading I or D
 		int32_t l = p->cigar[0] >> 4;
-		if ((p->cigar[0]&0xf) == 1) {
+		if ((p->cigar[0]&0xf) == MM_CIGAR_INS) {
 			if (r->rev) r->qe -= l;
 			else r->qs += l;
 			*qshift = l;
@@ -174,7 +181,7 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 	if (r->p == 0) return;
 	for (k = 0; k < r->p->n_cigar; ++k) {
 		uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4;
-		if (op == 0) {
+		if (op == MM_CIGAR_MATCH) {
 			while (len > 0) {
 				for (l = 0; l < len && qseq[qoff + l] == tseq[toff + l]; ++l) {} // run of "="; TODO: N<=>N is converted to "="
 				if (l > 0) { ++n_EQX; len -= l; toff += l; qoff += l; }
@@ -183,11 +190,11 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 				if (l > 0) { ++n_EQX; len -= l; toff += l; qoff += l; }
 			}
 			++n_M;
-		} else if (op == 1) { // insertion
+		} else if (op == MM_CIGAR_INS) {
 			qoff += len;
-		} else if (op == 2) { // deletion
+		} else if (op == MM_CIGAR_DEL) {
 			toff += len;
-		} else if (op == 3) { // intron
+		} else if (op == MM_CIGAR_N_SKIP) {
 			toff += len;
 		}
 	}
@@ -195,7 +202,7 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 	if (n_EQX == n_M) {
 		for (k = 0; k < r->p->n_cigar; ++k) {
 			uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4;
-			if (op == 0) r->p->cigar[k] = len << 4 | 7;
+			if (op == MM_CIGAR_MATCH) r->p->cigar[k] = len << 4 | MM_CIGAR_EQ_MATCH;
 		}
 		return;
 	}
@@ -209,25 +216,25 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 	toff = qoff = m = 0;
 	for (k = 0; k < r->p->n_cigar; ++k) {
 		uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4;
-		if (op == 0) { // match/mismatch
+		if (op == MM_CIGAR_MATCH) {
 			while (len > 0) {
 				// match
 				for (l = 0; l < len && qseq[qoff + l] == tseq[toff + l]; ++l) {}
-				if (l > 0) p->cigar[m++] = l << 4 | 7;
+				if (l > 0) p->cigar[m++] = l << 4 | MM_CIGAR_EQ_MATCH;
 				len -= l;
 				toff += l, qoff += l;
 				// mismatch
 				for (l = 0; l < len && qseq[qoff + l] != tseq[toff + l]; ++l) {}
-				if (l > 0) p->cigar[m++] = l << 4 | 8;
+				if (l > 0) p->cigar[m++] = l << 4 | MM_CIGAR_X_MISMATCH;
 				len -= l;
 				toff += l, qoff += l;
 			}
 			continue;
-		} else if (op == 1) { // insertion
+		} else if (op == MM_CIGAR_INS) {
 			qoff += len;
-		} else if (op == 2) { // deletion
+		} else if (op == MM_CIGAR_DEL) {
 			toff += len;
-		} else if (op == 3) { // intron
+		} else if (op == MM_CIGAR_N_SKIP) {
 			toff += len;
 		}
 		p->cigar[m++] = r->p->cigar[k];
@@ -237,10 +244,11 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 	r->p = p;
 }
 
-static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e, int is_eqx)
+static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e, int is_eqx, int log_gap)
 {
 	uint32_t k, l;
-	int32_t s = 0, max = 0, qshift, tshift, toff = 0, qoff = 0;
+	int32_t qshift, tshift, toff = 0, qoff = 0;
+	double s = 0.0, max = 0.0;
 	mm_extra_t *p = r->p;
 	if (p == 0) return;
 	mm_fix_cigar(r, qseq, tseq, &qshift, &tshift);
@@ -248,7 +256,7 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
 	r->blen = r->mlen = 0;
 	for (k = 0; k < p->n_cigar; ++k) {
 		uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4;
-		if (op == 0) { // match/mismatch
+		if (op == MM_CIGAR_MATCH) {
 			int n_ambi = 0, n_diff = 0;
 			for (l = 0; l < len; ++l) {
 				int cq = qseq[qoff + l], ct = tseq[toff + l];
@@ -260,27 +268,29 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
 			}
 			r->blen += len - n_ambi, r->mlen += len - (n_ambi + n_diff), p->n_ambi += n_ambi;
 			toff += len, qoff += len;
-		} else if (op == 1) { // insertion
+		} else if (op == MM_CIGAR_INS) {
 			int n_ambi = 0;
 			for (l = 0; l < len; ++l)
 				if (qseq[qoff + l] > 3) ++n_ambi;
 			r->blen += len - n_ambi, p->n_ambi += n_ambi;
-			s -= q + e * len;
+			if (log_gap) s -= q + (double)e * mg_log2(1.0 + len);
+			else s -= q + e;
 			if (s < 0) s = 0;
 			qoff += len;
-		} else if (op == 2) { // deletion
+		} else if (op == MM_CIGAR_DEL) {
 			int n_ambi = 0;
 			for (l = 0; l < len; ++l)
 				if (tseq[toff + l] > 3) ++n_ambi;
 			r->blen += len - n_ambi, p->n_ambi += n_ambi;
-			s -= q + e * len;
+			if (log_gap) s -= q + (double)e * mg_log2(1.0 + len);
+			else s -= q + e;
 			if (s < 0) s = 0;
 			toff += len;
-		} else if (op == 3) { // intron
+		} else if (op == MM_CIGAR_N_SKIP) {
 			toff += len;
 		}
 	}
-	p->dp_max = max;
+	p->dp_max = (int32_t)(max + .499);
 	assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
 	if (is_eqx) mm_update_cigar_eqx(r, qseq, tseq); // NB: it has to be called here as changes to qseq and tseq are not returned
 }
@@ -310,6 +320,7 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) //
 	}
 }
 
+#if 0
 static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const uint8_t *junc, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
 {
 	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
@@ -333,11 +344,69 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 		int i;
 		fprintf(stderr, "score=%d, cigar=", ez->score);
 		for (i = 0; i < ez->n_cigar; ++i)
-			fprintf(stderr, "%d%c", ez->cigar[i]>>4, "MIDN"[ez->cigar[i]&0xf]);
+			fprintf(stderr, "%d%c", ez->cigar[i]>>4, MM_CIGAR_STR[ez->cigar[i]&0xf]);
 		fprintf(stderr, "\n");
 	}
 }
+#endif
+#if 1
+static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const uint8_t *junc, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
+{
+#ifdef MANUAL_PROFILING
+	uint64_t align_start = __rdtsc();
+#endif
 
+	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
+		int i;
+		fprintf(stderr, "===> q=(%d,%d), e=(%d,%d), bw=%d, flag=%d, zdrop=%d <===\n", opt->q, opt->q2, opt->e, opt->e2, w, flag, opt->zdrop);
+		for (i = 0; i < tlen; ++i) fputc("ACGTN"[tseq[i]], stderr);
+		fputc('\n', stderr);
+		for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr);
+		fputc('\n', stderr);
+	}
+	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)
+		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);
+	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{
+#if defined (ALIGN_AVX) && (defined(__AVX512BW__) || (defined(__AVX2__) && defined(APPLY_AVX2)))
+#ifdef __AVX512BW__
+		
+            	ksw_extd2_avx512(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez);
+#elif __AVX2__
+	avg = 0;
+//	uint64_t *ptr_km = (uint64_t *) km1;
+//	for(uint64_t itr = 0; itr < km_size/512; itr++){
+//		avg+=ptr_km[itr];
+//	}
+
+//#ifdef MANUAL_PROFILING
+//	uint64_t align_start = __rdtsc();
+//#endif
+		ksw_extd2_avx2(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez);
+//#ifdef MANUAL_PROFILING
+//	alignment_time += (__rdtsc() - align_start);
+//#endif
+#endif
+#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);
+#endif
+	   }
+	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
+		int i;
+		fprintf(stderr, "score=%d, cigar=", ez->score);
+		for (i = 0; i < ez->n_cigar; ++i)
+			fprintf(stderr, "%d%c", ez->cigar[i]>>4, "MIDN"[ez->cigar[i]&0xf]);
+		fprintf(stderr, "\n");
+	}
+#ifdef MANUAL_PROFILING
+	alignment_time += (__rdtsc() - align_start);
+#endif
+}
+#endif
 static inline int mm_get_hplen_back(const mm_idx_t *mi, uint32_t rid, uint32_t x)
 {
 	int64_t i, off0 = mi->seq[rid].offset, off = off0 + x;
@@ -533,8 +602,13 @@ static int mm_seed_ext_score(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 	re = re + ext_len < (int32_t)mi->seq[rid].len? re + ext_len : mi->seq[rid].len;
 	qe = qe + ext_len < qlen? qe + ext_len : qlen;
 	tseq = (uint8_t*)kmalloc(km, re - rs);
-	mm_idx_getseq(mi, rid, rs, re, tseq);
-	qseq = qseq0[a->x>>63] + qs;
+	if (opt->flag & MM_F_QSTRAND) {
+		qseq = qseq0[0] + qs;
+		mm_idx_getseq2(mi, a->x>>63, rid, rs, re, tseq);
+	} else {
+		qseq = qseq0[a->x>>63] + qs;
+		mm_idx_getseq(mi, rid, rs, re, tseq);
+	}
 	qp = ksw_ll_qinit(km, 2, qe - qs, qseq, 5, mat);
 	score = ksw_ll_i16(qp, re - rs, tseq, opt->q, opt->e, &q_off, &t_off);
 	kfree(km, tseq);
@@ -690,8 +764,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	junc = (uint8_t*)kmalloc(km, re0 - rs0);
 
 	if (qs > 0 && rs > 0) { // left extension; probably the condition can be changed to "qs > qs0 && rs > rs0"
-		qseq = &qseq0[rev][qs0];
-		mm_idx_getseq(mi, rid, rs0, rs, tseq);
+		if (opt->flag & MM_F_QSTRAND) {
+			qseq = &qseq0[0][qs0];
+			mm_idx_getseq2(mi, rev, rid, rs0, rs, tseq);
+		} else {
+			qseq = &qseq0[rev][qs0];
+			mm_idx_getseq(mi, rid, rs0, rs, tseq);
+		}
 		mm_idx_bed_junc(mi, rid, rs0, rs, junc);
 		mm_seq_rev(qs - qs0, qseq);
 		mm_seq_rev(rs - rs0, tseq);
@@ -720,8 +799,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			if (a[as1+i].y & MM_SEED_LONG_JOIN)
 				bw1 = qe - qs > re - rs? qe - qs : re - rs;
 			// perform alignment
-			qseq = &qseq0[rev][qs];
-			mm_idx_getseq(mi, rid, rs, re, tseq);
+			if (opt->flag & MM_F_QSTRAND) {
+				qseq = &qseq0[0][qs];
+				mm_idx_getseq2(mi, rev, rid, rs, re, tseq);
+			} else {
+				qseq = &qseq0[rev][qs];
+				mm_idx_getseq(mi, rid, rs, re, tseq);
+			}
 			mm_idx_bed_junc(mi, rid, rs, re, junc);
 			if (is_sr) { // perform ungapped alignment
 				assert(qe - qs == re - rs);
@@ -730,7 +814,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 					if (qseq[j] >= 4 || tseq[j] >= 4) ez->score += opt->e2;
 					else ez->score += qseq[j] == tseq[j]? opt->a : -opt->b;
 				}
-				ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, 0, qe - qs);
+				ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, MM_CIGAR_MATCH, qe - qs);
 			} else { // perform normal gapped alignment
 				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
 			}
@@ -757,7 +841,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				re1 = rs + (ez->max_t + 1);
 				qe1 = qs + (ez->max_q + 1);
 				if (cnt1 - (j + 1) >= opt->min_cnt) {
-					mm_split_reg(r, r2, as1 + j + 1 - r->as, qlen, a);
+					mm_split_reg(r, r2, as1 + j + 1 - r->as, qlen, a, !!(opt->flag&MM_F_QSTRAND));
 					if (zdrop_code == 2) r2->split_inv = 1;
 				}
 				break;
@@ -767,8 +851,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	}
 
 	if (!dropped && qe < qe0 && re < re0) { // right extension
-		qseq = &qseq0[rev][qe];
-		mm_idx_getseq(mi, rid, re, re0, tseq);
+		if (opt->flag & MM_F_QSTRAND) {
+			qseq = &qseq0[0][qe];
+			mm_idx_getseq2(mi, rev, rid, re, re0, tseq);
+		} else {
+			qseq = &qseq0[rev][qe];
+			mm_idx_getseq(mi, rid, re, re0, tseq);
+		}
 		mm_idx_bed_junc(mi, rid, re, re0, junc);
 		mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, junc, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
 		if (ez->n_cigar > 0) {
@@ -781,13 +870,19 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	assert(qe1 <= qlen);
 
 	r->rs = rs1, r->re = re1;
-	if (rev) r->qs = qlen - qe1, r->qe = qlen - qs1;
-	else r->qs = qs1, r->qe = qe1;
+	if (!rev || (opt->flag & MM_F_QSTRAND)) r->qs = qs1, r->qe = qe1;
+	else r->qs = qlen - qe1, r->qe = qlen - qs1;
 
 	assert(re1 - rs1 <= re0 - rs0);
 	if (r->p) {
-		mm_idx_getseq(mi, rid, rs1, re1, tseq);
-		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e, opt->flag & MM_F_EQX);
+		if (opt->flag & MM_F_QSTRAND) {
+			mm_idx_getseq2(mi, r->rev, rid, rs1, re1, tseq);
+			qseq = &qseq0[0][qs1];
+		} else {
+			mm_idx_getseq(mi, rid, rs1, re1, tseq);
+			qseq = &qseq0[r->rev][qs1];
+		}
+		mm_update_extra(r, qseq, tseq, mat, opt->q, opt->e, opt->flag & MM_F_EQX, !(opt->flag & MM_F_SR));
 		if (rev && r->p->trans_strand)
 			r->p->trans_strand ^= 3; // flip to the read strand
 	}
@@ -797,7 +892,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 }
 
 static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez)
-{
+{ // NB: this doesn't work with the qstrand mode
 	int tl, ql, score, ret = 0, q_off, t_off;
 	uint8_t *tseq, *qseq;
 	int8_t mat[25];
@@ -846,7 +941,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	}
 	r_inv->rs = r1->re + t_off;
 	r_inv->re = r_inv->rs + ez->max_t + 1;
-	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, opt->flag & MM_F_EQX);
+	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, opt->flag & MM_F_EQX, !(opt->flag & MM_F_SR));
 	ret = 1;
 end_align1_inv:
 	kfree(km, tseq);
@@ -863,6 +958,71 @@ static inline mm_reg1_t *mm_insert_reg(const mm_reg1_t *r, int i, int *n_regs, m
 	return regs;
 }
 
+static inline void mm_count_gaps(const mm_reg1_t *r, int32_t *n_gap_, int32_t *n_gapo_)
+{
+	uint32_t i;
+	int32_t n_gapo = 0, n_gap = 0;
+	*n_gap_ = *n_gapo_ = -1;
+	if (r->p == 0) return;
+	for (i = 0; i < r->p->n_cigar; ++i) {
+		int32_t op = r->p->cigar[i] & 0xf, len = r->p->cigar[i] >> 4;
+		if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL)
+			++n_gapo, n_gap += len;
+	}
+	*n_gap_ = n_gap, *n_gapo_ = n_gapo;
+}
+
+double mm_event_identity(const mm_reg1_t *r)
+{
+	int32_t n_gap, n_gapo;
+	if (r->p == 0) return -1.0f;
+	mm_count_gaps(r, &n_gap, &n_gapo);
+	return (double)r->mlen / (r->blen + r->p->n_ambi - n_gap + n_gapo);
+}
+
+static int32_t mm_recal_max_dp(const mm_reg1_t *r, double b2, int32_t match_sc)
+{
+	uint32_t i;
+	int32_t n_gap = 0, n_gapo = 0, n_mis;
+	double gap_cost = 0.0;
+	if (r->p == 0) return -1;
+	for (i = 0; i < r->p->n_cigar; ++i) {
+		int32_t op = r->p->cigar[i] & 0xf, len = r->p->cigar[i] >> 4;
+		if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL) {
+			gap_cost += b2 + (double)mg_log2(1.0 + len);
+			++n_gapo, n_gap += len;
+		}
+	}
+	n_mis = r->blen + r->p->n_ambi - r->mlen - n_gap;
+	return (int32_t)(match_sc * (r->mlen - b2 * n_mis - gap_cost) + .499);
+}
+
+void mm_update_dp_max(int qlen, int n_regs, mm_reg1_t *regs, float frac, int a, int b)
+{
+	int32_t max = -1, max2 = -1, i, max_i = -1;
+	double div, b2;
+	if (n_regs < 2) return;
+	for (i = 0; i < n_regs; ++i) {
+		mm_reg1_t *r = &regs[i];
+		if (r->p == 0) continue;
+		if (r->p->dp_max > max) max2 = max, max = r->p->dp_max, max_i = i;
+		else if (r->p->dp_max > max2) max2 = r->p->dp_max;
+	}
+	if (max_i < 0 || max < 0 || max2 < 0) return;
+	if (regs[max_i].qe - regs[max_i].qs < (double)qlen * frac) return;
+	if (max2 < (double)max * frac) return;
+	div = 1. - mm_event_identity(&regs[max_i]);
+	if (div < 0.02) div = 0.02;
+	b2 = 0.5 / div; // max value: 25
+	if (b2 * a < b) b2 = (double)a / b;
+	for (i = 0; i < n_regs; ++i) {
+		mm_reg1_t *r = &regs[i];
+		if (r->p == 0) continue;
+		r->p->dp_max = mm_recal_max_dp(r, b2, a);
+		if (r->p->dp_max < 0) r->p->dp_max = 0;
+	}
+}
+
 mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a)
 {
 	extern unsigned char seq_nt4_table[256];
@@ -906,7 +1066,7 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 				regs[i].p->trans_strand = opt->flag&MM_F_SPLICE_FOR? 1 : 2;
 		}
 		if (r2.cnt > 0) regs = mm_insert_reg(&r2, i, &n_regs, regs);
-		if (i > 0 && regs[i].split_inv) {
+		if (i > 0 && regs[i].split_inv && !(opt->flag & MM_F_NO_INV)) {
 			if (mm_align1_inv(km, opt, mi, qlen, qseq0, &regs[i-1], &regs[i], &r2, &ez)) {
 				regs = mm_insert_reg(&r2, i, &n_regs, regs);
 				++i; // skip the inserted INV alignment
@@ -917,6 +1077,10 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 	kfree(km, qseq0[0]);
 	kfree(km, ez.cigar);
 	mm_filter_regs(opt, qlen, n_regs_, regs);
+	if (!(opt->flag&MM_F_SR) && !opt->split_prefix && qlen >= opt->rank_min_len) {
+		mm_update_dp_max(qlen, *n_regs_, regs, opt->rank_frac, opt->a, opt->b);
+		mm_filter_regs(opt, qlen, n_regs_, regs);
+	}
 	mm_hit_sort(km, n_regs_, regs, opt->alt_drop);
 	return regs;
 }

build_rmi.sh

@@ -0,0 +1,16 @@
+ref_data=$1
+preset=$2
+
+make clean && make lhash_index=1
+touch temp_read.fastq
+./minimap2  -ax $2 $1 temp_read.fastq >/dev/null      
+
+kv_file=$1"_"$2"_minimizers_key_value_sorted"  
+
+full_path=`readlink -f $kv_file`
+
+cd ./ext/TAL
+make lisa_hash
+./build-lisa-hash-index $full_path
+
+rm ../../temp_read.fastq

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

example.c

@@ -47,7 +47,7 @@ int main(int argc, char *argv[])
 				printf("%s\t%d\t%d\t%d\t%c\t", ks->name.s, ks->seq.l, r->qs, r->qe, "+-"[r->rev]);
 				printf("%s\t%d\t%d\t%d\t%d\t%d\t%d\tcg:Z:", mi->seq[r->rid].name, mi->seq[r->rid].len, r->rs, r->re, r->mlen, r->blen, r->mapq);
 				for (i = 0; i < r->p->n_cigar; ++i) // IMPORTANT: this gives the CIGAR in the aligned regions. NO soft/hard clippings!
-					printf("%d%c", r->p->cigar[i]>>4, "MIDNSH"[r->p->cigar[i]&0xf]);
+					printf("%d%c", r->p->cigar[i]>>4, MM_CIGAR_STR[r->p->cigar[i]&0xf]);
 				putchar('\n');
 				free(r->p);
 			}

format.c

@@ -144,8 +144,8 @@ static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 	if (write_tag) mm_sprintf_lite(s, "\tcs:Z:");
 	for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
 		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
-		assert((op >= 0 && op <= 3) || op == 7 || op == 8);
-		if (op == 0 || op == 7 || op == 8) { // match
+		assert((op >= MM_CIGAR_MATCH && op <= MM_CIGAR_N_SKIP) || op == MM_CIGAR_EQ_MATCH || op == MM_CIGAR_X_MISMATCH);
+		if (op == MM_CIGAR_MATCH || op == MM_CIGAR_EQ_MATCH || op == MM_CIGAR_X_MISMATCH) {
 			int l_tmp = 0;
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
@@ -166,12 +166,12 @@ static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 				} else mm_sprintf_lite(s, ":%d", l_tmp);
 			}
 			q_off += len, t_off += len;
-		} else if (op == 1) { // insertion to ref
+		} else if (op == MM_CIGAR_INS) {
 			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 == 2) { // deletion from ref
+		} else if (op == MM_CIGAR_DEL) {
 			for (j = 0, tmp[len] = 0; j < len; ++j)
 				tmp[j] = "acgtn"[tseq[t_off + j]];
 			mm_sprintf_lite(s, "-%s", tmp);
@@ -192,8 +192,8 @@ static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 	if (write_tag) mm_sprintf_lite(s, "\tMD:Z:");
 	for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
 		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
-		assert((op >= 0 && op <= 3) || op == 7 || op == 8);
-		if (op == 0 || op == 7 || op == 8) { // match
+		assert((op >= MM_CIGAR_MATCH && op <= MM_CIGAR_N_SKIP) || op == MM_CIGAR_EQ_MATCH || op == MM_CIGAR_X_MISMATCH);
+		if (op == MM_CIGAR_MATCH || op == MM_CIGAR_EQ_MATCH || op == MM_CIGAR_X_MISMATCH) {
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
 					mm_sprintf_lite(s, "%d%c", l_MD, "ACGTN"[tseq[t_off + j]]);
@@ -201,15 +201,15 @@ static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 				} else ++l_MD;
 			}
 			q_off += len, t_off += len;
-		} else if (op == 1) { // insertion to ref
+		} else if (op == MM_CIGAR_INS) {
 			q_off += len;
-		} else if (op == 2) { // deletion from ref
+		} else if (op == MM_CIGAR_DEL) {
 			for (j = 0, tmp[len] = 0; j < len; ++j)
 				tmp[j] = "ACGTN"[tseq[t_off + j]];
 			mm_sprintf_lite(s, "%d^%s", l_MD, tmp);
 			l_MD = 0;
 			t_off += len;
-		} else if (op == 3) { // reference skip
+		} else if (op == MM_CIGAR_N_SKIP) {
 			t_off += len;
 		}
 	}
@@ -217,7 +217,7 @@ static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 	assert(t_off == r->re - r->rs && q_off == r->qe - r->qs);
 }
 
-static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int no_iden, int is_MD, int write_tag)
+static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int no_iden, int is_MD, int write_tag, int is_qstrand)
 {
 	extern unsigned char seq_nt4_table[256];
 	int i;
@@ -227,14 +227,20 @@ static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_
 	qseq = (uint8_t*)kmalloc(km, r->qe - r->qs);
 	tseq = (uint8_t*)kmalloc(km, r->re - r->rs);
 	tmp = (char*)kmalloc(km, r->re - r->rs > r->qe - r->qs? r->re - r->rs + 1 : r->qe - r->qs + 1);
-	mm_idx_getseq(mi, r->rid, r->rs, r->re, tseq);
-	if (!r->rev) {
+	if (is_qstrand) {
+		mm_idx_getseq2(mi, r->rev, r->rid, r->rs, r->re, tseq);
 		for (i = r->qs; i < r->qe; ++i)
 			qseq[i - r->qs] = seq_nt4_table[(uint8_t)t->seq[i]];
 	} else {
-		for (i = r->qs; i < r->qe; ++i) {
-			uint8_t c = seq_nt4_table[(uint8_t)t->seq[i]];
-			qseq[r->qe - i - 1] = c >= 4? 4 : 3 - c;
+		mm_idx_getseq(mi, r->rid, r->rs, r->re, tseq);
+		if (!r->rev) {
+			for (i = r->qs; i < r->qe; ++i)
+				qseq[i - r->qs] = seq_nt4_table[(uint8_t)t->seq[i]];
+		} else {
+			for (i = r->qs; i < r->qe; ++i) {
+				uint8_t c = seq_nt4_table[(uint8_t)t->seq[i]];
+				qseq[r->qe - i - 1] = c >= 4? 4 : 3 - c;
+			}
 		}
 	}
 	if (is_MD) write_MD_core(s, tseq, qseq, r, tmp, write_tag);
@@ -242,14 +248,14 @@ static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_
 	kfree(km, qseq); kfree(km, tseq); kfree(km, tmp);
 }
 
-int mm_gen_cs_or_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int is_MD, int no_iden)
+int mm_gen_cs_or_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int is_MD, int no_iden, int is_qstrand)
 {
 	mm_bseq1_t t;
 	kstring_t str;
 	str.s = *buf, str.l = 0, str.m = *max_len;
 	t.l_seq = strlen(seq);
 	t.seq = (char*)seq;
-	write_cs_or_MD(km, &str, mi, &t, r, no_iden, is_MD, 0);
+	write_cs_or_MD(km, &str, mi, &t, r, no_iden, is_MD, 0, is_qstrand);
 	*max_len = str.m;
 	*buf = str.s;
 	return str.l;
@@ -257,24 +263,12 @@ int mm_gen_cs_or_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, cons
 
 int mm_gen_cs(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int no_iden)
 {
-	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 0, no_iden);
+	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 0, no_iden, 0);
 }
 
 int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq)
 {
-	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 1, 0);
-}
-
-double mm_event_identity(const mm_reg1_t *r)
-{
-	int32_t i, n_gapo = 0, n_gap = 0;
-	if (r->p == 0) return -1.0f;
-	for (i = 0; i < r->p->n_cigar; ++i) {
-		int32_t op = r->p->cigar[i] & 0xf, len = r->p->cigar[i] >> 4;
-		if (op == 1 || op == 2)
-			++n_gapo, n_gap += len;
-	}
-	return (double)r->mlen / (r->blen + r->p->n_ambi - n_gap + n_gapo);
+	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 1, 0, 0);
 }
 
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
@@ -305,7 +299,7 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 	if (r->split) mm_sprintf_lite(s, "\tzd:i:%d", r->split);
 }
 
-void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag, int rep_len)
+void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int64_t opt_flag, int rep_len)
 {
 	s->l = 0;
 	if (r == 0) {
@@ -316,7 +310,11 @@ void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const
 	mm_sprintf_lite(s, "%s\t%d\t%d\t%d\t%c\t", t->name, t->l_seq, r->qs, r->qe, "+-"[r->rev]);
 	if (mi->seq[r->rid].name) mm_sprintf_lite(s, "%s", mi->seq[r->rid].name);
 	else mm_sprintf_lite(s, "%d", r->rid);
-	mm_sprintf_lite(s, "\t%d\t%d\t%d", mi->seq[r->rid].len, r->rs, r->re);
+	mm_sprintf_lite(s, "\t%d", mi->seq[r->rid].len);
+	if ((opt_flag & MM_F_QSTRAND) && r->rev)
+		mm_sprintf_lite(s, "\t%d\t%d", mi->seq[r->rid].len - r->re, mi->seq[r->rid].len - r->rs);
+	else
+		mm_sprintf_lite(s, "\t%d\t%d", r->rs, r->re);
 	mm_sprintf_lite(s, "\t%d\t%d", r->mlen, r->blen);
 	mm_sprintf_lite(s, "\t%d", r->mapq);
 	write_tags(s, r);
@@ -325,15 +323,15 @@ void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const
 		uint32_t k;
 		mm_sprintf_lite(s, "\tcg:Z:");
 		for (k = 0; k < r->p->n_cigar; ++k)
-			mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDNSHP=XB"[r->p->cigar[k]&0xf]);
+			mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, MM_CIGAR_STR[r->p->cigar[k]&0xf]);
 	}
 	if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_MD)))
-		write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1);
+		write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1, !!(opt_flag&MM_F_QSTRAND));
 	if ((opt_flag & MM_F_COPY_COMMENT) && t->comment)
 		mm_sprintf_lite(s, "\t%s", t->comment);
 }
 
-void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag)
+void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int64_t opt_flag)
 {
 	mm_write_paf3(s, mi, t, r, km, opt_flag, -1);
 }
@@ -362,7 +360,7 @@ static inline const mm_reg1_t *get_sam_pri(int n_regs, const mm_reg1_t *regs)
 	return NULL;
 }
 
-static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, const mm_reg1_t *r, int opt_flag)
+static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, const mm_reg1_t *r, int64_t opt_flag)
 {
 	if (r->p == 0) {
 		mm_sprintf_lite(s, "*");
@@ -382,13 +380,13 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
 			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)
-				mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDNSHP=XB"[r->p->cigar[k]&0xf]);
+				mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, MM_CIGAR_STR[r->p->cigar[k]&0xf]);
 			if (clip_len[1]) mm_sprintf_lite(s, "%d%c", clip_len[1], clip_char);
 		}
 	}
 }
 
-void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag, int rep_len)
+void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int64_t opt_flag, int rep_len)
 {
 	const int max_bam_cigar_op = 65535;
 	int flag, n_regs = n_regss[seg_idx], cigar_in_tag = 0;
@@ -535,7 +533,7 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 			}
 		}
 		if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_MD)))
-			write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1);
+			write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1, 0);
 		if (cigar_in_tag)
 			write_sam_cigar(s, flag, 1, t->l_seq, r, opt_flag);
 	}
@@ -547,7 +545,7 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 	s->s[s->l] = 0; // we always have room for an extra byte (see str_enlarge)
 }
 
-void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag)
+void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int64_t opt_flag)
 {
 	mm_write_sam3(s, mi, t, seg_idx, reg_idx, n_seg, n_regss, regss, km, opt_flag, -1);
 }

hit.c

@@ -20,14 +20,14 @@ static inline void mm_cal_fuzzy_len(mm_reg1_t *r, const mm128_t *a)
 	}
 }
 
-static inline void mm_reg_set_coor(mm_reg1_t *r, int32_t qlen, const mm128_t *a)
+static inline void mm_reg_set_coor(mm_reg1_t *r, int32_t qlen, const mm128_t *a, int is_qstrand)
 { // NB: r->as and r->cnt MUST BE set correctly for this function to work
 	int32_t k = r->as, q_span = (int32_t)(a[k].y>>32&0xff);
 	r->rev = a[k].x>>63;
 	r->rid = a[k].x<<1>>33;
 	r->rs = (int32_t)a[k].x + 1 > q_span? (int32_t)a[k].x + 1 - q_span : 0; // NB: target span may be shorter, so this test is necessary
 	r->re = (int32_t)a[k + r->cnt - 1].x + 1;
-	if (!r->rev) {
+	if (!r->rev || is_qstrand) {
 		r->qs = (int32_t)a[k].y + 1 - q_span;
 		r->qe = (int32_t)a[k + r->cnt - 1].y + 1;
 	} else {
@@ -49,7 +49,7 @@ static inline uint64_t hash64(uint64_t key)
 	return key;
 }
 
-mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a) // convert chains to hits
+mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a, int is_qstrand) // convert chains to hits
 {
 	mm128_t *z, tmp;
 	mm_reg1_t *r;
@@ -81,7 +81,7 @@ mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u,
 		ri->cnt = (int32_t)z[i].y;
 		ri->as = z[i].y >> 32;
 		ri->div = -1.0f;
-		mm_reg_set_coor(ri, qlen, a);
+		mm_reg_set_coor(ri, qlen, a, is_qstrand);
 	}
 	kfree(km, z);
 	return r;
@@ -103,7 +103,7 @@ static inline int mm_alt_score(int score, float alt_diff_frac)
 	return score > 0? score : 1;
 }
 
-void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a)
+void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a, int is_qstrand)
 {
 	if (n <= 0 || n >= r->cnt) return;
 	*r2 = *r;
@@ -115,10 +115,10 @@ void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a)
 	r2->score = (int32_t)(r->score * ((float)r2->cnt / r->cnt) + .499);
 	r2->as = r->as + n;
 	if (r->parent == r->id) r2->parent = MM_PARENT_TMP_PRI;
-	mm_reg_set_coor(r2, qlen, a);
+	mm_reg_set_coor(r2, qlen, a, is_qstrand);
 	r->cnt -= r2->cnt;
 	r->score -= r2->score;
-	mm_reg_set_coor(r, qlen, a);
+	mm_reg_set_coor(r, qlen, a, is_qstrand);
 	r->split |= 1, r2->split |= 2;
 }
 
@@ -358,7 +358,7 @@ mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int
 		}
 	}
 	for (s = 0; s < n_segs; ++s) {
-		regs[s] = mm_gen_regs(km, hash, qlens[s], seg[s].n_u, seg[s].u, seg[s].a);
+		regs[s] = mm_gen_regs(km, hash, qlens[s], seg[s].n_u, seg[s].u, seg[s].a, 0);
 		n_regs[s] = seg[s].n_u;
 		for (i = 0; i < n_regs[s]; ++i) {
 			regs[s][i].seg_split = 1;

index.c

@@ -14,6 +14,21 @@
 #include "mmpriv.h"
 #include "kvec.h"
 #include "khash.h"
+#include <map>
+#include <fstream>
+#include <vector>
+#include <algorithm>
+#include <x86intrin.h>
+using namespace std;
+
+
+extern uint64_t minimizer_lookup_time, alignment_time, dp_time, rmq_time, rmq_t1, rmq_t2, rmq_t3, rmq_t4;
+#ifdef LISA_HASH
+#include "lisa_hash.h"
+extern lisa_hash<uint64_t, uint64_t> *lh;
+#endif
+
+
 
 #define idx_hash(a) ((a)>>1)
 #define idx_eq(a, b) ((a)>>1 == (b)>>1)
@@ -52,9 +67,43 @@ mm_idx_t *mm_idx_init(int w, int k, int b, int flag)
 	if (!(mm_dbg_flag & 1)) mi->km = km_init();
 	return mi;
 }
+void mm_idx_destroy_mm_hash(mm_idx_t *mi)
+{
+	//fprintf(stderr, "mm_destroy_hash\n");
+	uint32_t i;
+	if (mi == 0) return;
+	if (mi->h) kh_destroy(str, (khash_t(str)*)mi->h);
+	if (mi->B) {
+		for (i = 0; i < 1U<<mi->b; ++i) {
+			free(mi->B[i].p);
+			free(mi->B[i].a.a);
+			kh_destroy(idx, (idxhash_t*)mi->B[i].h);
+		}
+	}
+}
+void mm_idx_destroy_seq(mm_idx_t *mi)
+{
+	//fprintf(stderr, "mm_destroy_seq\n");
+
+	uint32_t i;
+	if (mi == 0) return;
+	if (mi->I) {
+		for (i = 0; i < mi->n_seq; ++i)
+			free(mi->I[i].a);
+		free(mi->I);
+	}
+	if (!mi->km) {
+		for (i = 0; i < mi->n_seq; ++i)
+			free(mi->seq[i].name);
+		free(mi->seq);
+	} else km_destroy(mi->km);
+	free(mi->B); free(mi->S); free(mi);
+}
+
 
 void mm_idx_destroy(mm_idx_t *mi)
 {
+
 	uint32_t i;
 	if (mi == 0) return;
 	if (mi->h) kh_destroy(str, (khash_t(str)*)mi->h);
@@ -96,6 +145,317 @@ const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n)
 		return &b->p[kh_val(h, k)>>32];
 	}
 }
+//Output minimap2's hash table entries
+class hash_entry {
+	public: 
+	uint64_t key;
+	uint64_t n;
+	uint64_t *p;
+	hash_entry(uint64_t k, uint64_t n_, uint64_t *p_){
+		key = k;
+		n = n_;
+		p = p_;
+	}
+
+};
+bool key_sort( hash_entry i1, hash_entry i2)
+{
+    return (i1.key < i2.key);
+} 
+
+#if 0
+void mm_idx_load_key_value_lisa(const char* f_name, const mm_idx_t *mi)  
+{
+	uint64_t tic = __rdtsc();
+	std::vector<hash_entry> v_hash;
+
+	//ofstream f(f_name);
+	fprintf(stderr, "Building sorted key-val map\n");
+
+	uint32_t i,j;
+	uint64_t num_values = 0;
+	for (i = 0; i < 1U<<mi->b; ++i) {
+		
+		
+		//fprintf(stderr, "BucketID %lu \n", i);
+		idxhash_t *h = (idxhash_t*)mi->B[i].h;
+		khint_t k;
+		if (h == 0) continue;
+		for (k = 0; k < kh_end(h); ++k){
+			if (kh_exist(h, k)) {
+				uint64_t key = kh_key(h, k), bucket_id = i;
+				key = key>>1;
+				
+				key = key<<mi->b | bucket_id;
+				
+				if(kh_key(h, k)&1)
+				{
+					//print key value
+					//fprintf(stderr, "%llu %llu %llu\n", key, kh_val(h, k), 0);
+					v_hash.push_back(hash_entry(key, kh_val(h, k), NULL));
+				}
+				else
+				{	// print key
+					uint32_t n = (uint32_t)kh_val(h, k);
+					//fprintf(stderr, "%llu %llu %llu ", key, kh_val(h, k), n);
+					// for 0 to lsb 32 val
+					//  print b->p[msb 32 of val]
+					
+					v_hash.push_back(hash_entry(key, n, &mi->B[i].p[(kh_val(h, k)>>32) + 0]));
+				}
+			}
+		}
+	}
+	sort(v_hash.begin(), v_hash.end(), key_sort);
+	fprintf(stderr, "Sorted map building time = %lld \n", __rdtsc() - tic);
+	fprintf(stderr, "Storing hash to %s \n", f_name);
+	tic = __rdtsc();
+
+	int64_t itr_p = 0;
+	for( int i = 0; i < v_hash.size(); i++){
+	
+		if(v_hash[i].p == NULL){	
+			//f<<v_hash[i].key << " "<<1<<"\n"<<v_hash[i].n<<" \n";
+			
+			lh->p[itr_p++] = v_hash[i].n;
+			continue;
+		}
+		
+		//f<<v_hash[i].key << " "<<v_hash[i].n<<endl;
+		
+		for(int j = 0; j < v_hash[i].n; j++){
+		//	f<<v_hash[i].p[j]<<" ";
+			lh->p[itr_p++] = v_hash[i].p[j];
+			num_values++;
+		}
+		//f<<endl;
+	
+	}
+
+	//f.close();	
+
+	string size_file_name = (string) f_name + "_size";
+	ofstream size_f(size_file_name);
+	size_f<<v_hash.size()<<" "<<num_values;
+	size_f.close();
+
+	string prefix = (string)f_name + "_keys";	
+	string keys_bin_file_name = prefix + ".uint64";		
+	ofstream wf(keys_bin_file_name, ios::out | ios::binary);
+	wf.write((char*)&key_list[0], (key_list.size())*sizeof(uint64_t));		
+	wf.close();	
+
+	key_list.clear();
+
+	m.clear();
+	v_hash.clear();
+
+	fprintf(stderr, "Index store File IO time %lld \n", __rdtsc() - tic);
+
+}
+#endif
+
+void mm_idx_dump_hash(const char* f_name, const mm_idx_t *mi)  
+{
+	uint64_t tic = __rdtsc();
+	//std::map<uint64_t, vector<uint64_t>> m;
+	std::vector<hash_entry> v_hash;
+
+	//ofstream f(f_name);
+	fprintf(stderr, "Building sorted key-val map\n");
+
+	uint32_t i,j;
+	uint64_t num_values = 0;
+	for (i = 0; i < 1U<<mi->b; ++i) {
+		
+		
+		//fprintf(stderr, "BucketID %lu \n", i);
+		idxhash_t *h = (idxhash_t*)mi->B[i].h;
+		khint_t k;
+		if (h == 0) continue;
+		for (k = 0; k < kh_end(h); ++k){
+			if (kh_exist(h, k)) {
+				uint64_t key = kh_key(h, k), bucket_id = i;
+				key = key>>1;
+				
+				key = key<<mi->b | bucket_id;
+				
+				if(kh_key(h, k)&1)
+				{
+					//print key value
+					//fprintf(stderr, "%llu %llu %llu\n", key, kh_val(h, k), 0);
+					//m[key].push_back(kh_val(h, k));
+					v_hash.push_back(hash_entry(key, kh_val(h, k), NULL));
+				}
+				else
+				{	// print key
+					uint32_t n = (uint32_t)kh_val(h, k);
+					//fprintf(stderr, "%llu %llu %llu ", key, kh_val(h, k), n);
+					// for 0 to lsb 32 val
+					//  print b->p[msb 32 of val]
+					
+					v_hash.push_back(hash_entry(key, n, &mi->B[i].p[(kh_val(h, k)>>32) + 0]));
+				}
+			}
+		}
+	}
+	sort(v_hash.begin(), v_hash.end(), key_sort);
+	fprintf(stderr, "Sorted map building time = %lld \n", __rdtsc() - tic);
+	fprintf(stderr, "Storing hash to %s \n", f_name);
+	tic = __rdtsc();
+
+	vector<uint64_t> key_list;
+	vector<uint64_t> val_list;
+	vector<uint64_t> p_list;
+/*
+	key_list.push_back(m.size());	
+	for(auto k : m){
+		key_list.push_back(k.first);
+		f<<k.first << " "<<k.second.size()<<endl;
+		for(int j = 0; j < k.second.size(); j++){
+			f<<k.second[j]<<" ";
+			num_values++;
+		}
+		f<<endl;
+	}
+*/
+
+
+	key_list.push_back(v_hash.size());
+	int64_t itr_p = 0;
+	uint64_t sum_pos = 0;
+	string f1_name = (string)f_name + "_pos_bin";
+	string f2_name = (string)f_name + "_val_bin";
+	ofstream f1(f1_name, ios::out | ios::binary);
+	ofstream f2(f2_name, ios::out | ios::binary);
+	for( int i = 0; i < v_hash.size(); i++){
+		key_list.push_back(v_hash[i].key);
+		if(v_hash[i].p == NULL){	
+			//f<<v_hash[i].key << " "<<1<<"\n"<<v_hash[i].n<<" \n";
+			val_list.push_back(sum_pos<<32|(uint64_t)1);
+			sum_pos+=1;
+			p_list.push_back(v_hash[i].n);
+			num_values++;
+			continue;
+		}
+		
+		//f<<v_hash[i].key << " "<<v_hash[i].n<<endl;
+		val_list.push_back(sum_pos<<32|(uint64_t)v_hash[i].n);
+		sum_pos+=v_hash[i].n;	
+			
+		num_values+=v_hash[i].n;
+
+
+		for(int j = 0; j < v_hash[i].n; j++){
+			//f<<v_hash[i].p[j]<<" ";
+			p_list.push_back(v_hash[i].p[j]);
+		}
+//		f<<endl;
+	
+	}
+	f1.write((char*)&val_list[0], (val_list.size())*sizeof(uint64_t));		
+	f2.write((char*)&p_list[0], (p_list.size())*sizeof(uint64_t));		
+	f1.close();
+	f2.close();
+	fprintf(stderr, "Index sorted SoA time %lld \n", __rdtsc() - tic);
+
+	//f.close();	
+
+	string size_file_name = (string) f_name + "_size";
+	ofstream size_f(size_file_name);
+	size_f<<v_hash.size()<<" "<<num_values;
+	size_f.close();
+
+	string prefix = (string)f_name + "_keys";	
+	string keys_bin_file_name = prefix + ".uint64";		
+	ofstream wf(keys_bin_file_name, ios::out | ios::binary);
+	wf.write((char*)&key_list[0], (key_list.size())*sizeof(uint64_t));		
+	wf.close();	
+
+	key_list.clear();
+
+	//m.clear();
+	v_hash.clear();
+
+	fprintf(stderr, "Index store File IO time %lld \n", __rdtsc() - tic);
+
+}
+void mm_idx_dump_hash_1(const char* f_name, const mm_idx_t *mi)  
+{
+	uint64_t tic = __rdtsc();
+	std::map<uint64_t, vector<uint64_t>> m;
+
+	ofstream f(f_name);
+	fprintf(stderr, "Building sorted key-val map\n");
+
+	uint32_t i,j;
+	uint64_t num_values = 0;
+	for (i = 0; i < 1U<<mi->b; ++i) {
+		
+		
+		//fprintf(stderr, "BucketID %lu \n", i);
+		idxhash_t *h = (idxhash_t*)mi->B[i].h;
+		khint_t k;
+		if (h == 0) continue;
+		for (k = 0; k < kh_end(h); ++k){
+			if (kh_exist(h, k)) {
+				uint64_t key = kh_key(h, k), bucket_id = i;
+				key = key>>1;
+				
+				key = key<<mi->b | bucket_id;
+				
+				if(kh_key(h, k)&1)
+				{
+					//print key value
+					//fprintf(stderr, "%llu %llu %llu\n", key, kh_val(h, k), 0);
+					m[key].push_back(kh_val(h, k));
+				}
+				else
+				{	// print key
+					uint32_t n = (uint32_t)kh_val(h, k);
+					//fprintf(stderr, "%llu %llu %llu ", key, kh_val(h, k), n);
+					// for 0 to lsb 32 val
+					//  print b->p[msb 32 of val]
+					for(j = 0; j < n; j++)
+					{
+						//fprintf(stderr, "%llu ", mi->B[i].p[(kh_val(h, k)>>32) + j]);
+						m[key].push_back(mi->B[i].p[(kh_val(h, k)>>32) + j]);
+					}
+				}
+			}
+		}
+	}
+	fprintf(stderr, "Sorted map building time = %lld \n", __rdtsc() - tic);
+	fprintf(stderr, "Storing hash to %s \n", f_name);
+	tic = __rdtsc();
+	vector<uint64_t> key_list;
+	key_list.push_back(m.size());	
+	for(auto k : m){
+		key_list.push_back(k.first);
+		f<<k.first << " "<<k.second.size()<<endl;
+		for(int j = 0; j < k.second.size(); j++){
+			f<<k.second[j]<<" ";
+			num_values++;
+		}
+		f<<endl;
+	}
+	f.close();	
+	string size_file_name = (string) f_name + "_size";
+	ofstream size_f(size_file_name);
+	size_f<<m.size()<<" "<<num_values;
+	size_f.close();
+
+	string prefix = (string)f_name + "_keys";	
+	string keys_bin_file_name = prefix + ".uint64";		
+	ofstream wf(keys_bin_file_name, ios::out | ios::binary);
+	wf.write((char*)&key_list[0], (key_list.size())*sizeof(uint64_t));		
+	wf.close();	
+
+	key_list.clear();
+	m.clear();
+	fprintf(stderr, "Index store File IO time %lld \n", __rdtsc() - tic);
+
+}
 
 void mm_idx_stat(const mm_idx_t *mi)
 {
@@ -119,6 +479,7 @@ void mm_idx_stat(const mm_idx_t *mi)
 	}
 	fprintf(stderr, "[M::%s::%.3f*%.2f] distinct minimizers: %d (%.2f%% are singletons); average occurrences: %.3lf; average spacing: %.3lf; total length: %ld\n",
 			__func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), n, 100.0*n1/n, (double)sum / n, (double)len / sum, (long)len);
+	fprintf(stderr, "minimizer-lookup: %lld dp: %lld rmq: %lld rmq_t1: %lld rmq_t2: %lld rmq_t3: %lld rmq_t4: %lld alignment: %lld \n", minimizer_lookup_time, dp_time, rmq_time, rmq_t1, rmq_t2, rmq_t3, rmq_t4, alignment_time);
 }
 
 int mm_idx_index_name(mm_idx_t *mi)
@@ -161,6 +522,28 @@ int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, ui
 	return en - st;
 }
 
+int mm_idx_getseq_rev(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq)
+{
+	uint64_t i, st1, en1;
+	const mm_idx_seq_t *s;
+	if (rid >= mi->n_seq || st >= mi->seq[rid].len) return -1;
+	s = &mi->seq[rid];
+	if (en > s->len) en = s->len;
+	st1 = s->offset + (s->len - en);
+	en1 = s->offset + (s->len - st);
+	for (i = st1; i < en1; ++i) {
+		uint8_t c = mm_seq4_get(mi->S, i);
+		seq[en1 - i - 1] = c < 4? 3 - c : c;
+	}
+	return en - st;
+}
+
+int mm_idx_getseq2(const mm_idx_t *mi, int is_rev, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq)
+{
+	if (is_rev) return mm_idx_getseq_rev(mi, rid, st, en, seq);
+	else return mm_idx_getseq(mi, rid, st, en, seq);
+}
+
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f)
 {
 	int i;

krmq.h

@@ -247,7 +247,7 @@ int main(void) {
 		unsigned char dir[KRMQ_MAX_DEPTH]; \
 		int i, d = 0, cmp; \
 		unsigned cnt = 0; \
-		fake.__head.p[0] = *root_, fake.__head.p[1] = 0; \
+		fake = **root_, fake.__head.p[0] = *root_, fake.__head.p[1] = 0; \
 		if (cnt_) *cnt_ = 0; \
 		if (x) { \
 			for (cmp = -1, p = &fake; cmp; cmp = __cmp(x, p)) { \

ksw2.h

@@ -16,6 +16,13 @@
 #define KSW_EZ_SPLICE_REV  0x200
 #define KSW_EZ_SPLICE_FLANK 0x400
 
+// The subset of CIGAR operators used by ksw code.
+// Use MM_CIGAR_* from minimap.h if you need the full list.
+#define KSW_CIGAR_MATCH  0
+#define KSW_CIGAR_INS    1
+#define KSW_CIGAR_DEL    2
+#define KSW_CIGAR_N_SKIP 3
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -137,13 +144,13 @@ static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int min_intro
 		else if (!(tmp >> (state + 2) & 1)) state = 0; // if requesting other states, _state_ stays the same if it is a continuation; otherwise, set to H
 		if (state == 0) state = tmp & 7; // TODO: probably this line can be merged into the "else if" line right above; not 100% sure
 		if (force_state >= 0) state = force_state;
-		if (state == 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 0, 1), --i, --j; // match
-		else if (state == 1 || (state == 3 && min_intron_len <= 0)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i; // deletion
-		else if (state == 3 && min_intron_len > 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 3, 1), --i; // intron
-		else cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, 1), --j; // insertion
+		if (state == 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_MATCH, 1), --i, --j;
+		else if (state == 1 || (state == 3 && min_intron_len <= 0)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_DEL, 1), --i;
+		else if (state == 3 && min_intron_len > 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_N_SKIP, 1), --i;
+		else cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_INS, 1), --j;
 	}
-	if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, min_intron_len > 0 && i >= min_intron_len? 3 : 2, i + 1); // first deletion
-	if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, j + 1); // first insertion
+	if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, min_intron_len > 0 && i >= min_intron_len? KSW_CIGAR_N_SKIP : KSW_CIGAR_DEL, i + 1); // first deletion
+	if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_INS, j + 1); // first insertion
 	if (!is_rev)
 		for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR
 			tmp = cigar[i], cigar[i] = cigar[n_cigar-1-i], cigar[n_cigar-1-i] = tmp;

ksw2_extd2_avx.h

@@ -0,0 +1,42 @@
+/* The MIT License
+
+Copyright (c) 2018-     Dana-Farber Cancer Institute
+              2017-2018 Broad Institute, Inc.
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
+BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
+ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+Modified Copyright (C) 2021 Intel Corporation
+   Contacts: Saurabh Kalikar <saurabh.kalikar@intel.com>; 
+	Vasimuddin Md <vasimuddin.md@intel.com>; Sanchit Misra <sanchit.misra@intel.com>; 
+	Chirag Jain <chirag@iisc.ac.in>; Heng Li <hli@jimmy.harvard.edu>
+*/
+#include <string.h>
+#include <stdio.h>
+#include <assert.h>
+#include "ksw2.h"
+#include <immintrin.h>
+#include <x86intrin.h>
+#include <smmintrin.h>
+#include <emmintrin.h>
+void ksw_extd2_avx512(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+                   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
+
+void ksw_extd2_avx2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
+                   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);

lchain.c

@@ -5,17 +5,16 @@
 #include "mmpriv.h"
 #include "kalloc.h"
 #include "krmq.h"
+#include <x86intrin.h>
+//#include "simd_chain.h"
+//#include "parallel_chaining_32_bit.h"
+#include "parallel_chaining_v2_22.h"
 
-static inline float mg_log2(float x) // NB: this doesn't work when x<2
-{
-	union { float f; uint32_t i; } z = { x };
-	float log_2 = ((z.i >> 23) & 255) - 128;
-	z.i &= ~(255 << 23);
-	z.i += 127 << 23;
-	log_2 += (-0.34484843f * z.f + 2.02466578f) * z.f - 0.67487759f;
-	return log_2;
-}
+#ifdef MANUAL_PROFILING
+extern uint64_t dp_time, rmq_time, rmq_t1, rmq_t2, rmq_t3, rmq_t4;
+#endif
 
+extern bool enable_vect_dp_chaining;
 uint64_t *mg_chain_backtrack(void *km, int64_t n, const int32_t *f, const int64_t *p, int32_t *v, int32_t *t, int32_t min_cnt, int32_t min_sc, int32_t *n_u_, int32_t *n_v_)
 {
 	mm128_t *z;
@@ -98,15 +97,65 @@ static mm128_t *compact_a(void *km, int32_t n_u, uint64_t *u, int32_t n_v, int32
 
 static inline int32_t comput_sc(const mm128_t *ai, const mm128_t *aj, int32_t max_dist_x, int32_t max_dist_y, int32_t bw, float chn_pen_gap, float chn_pen_skip, int is_cdna, int n_seg)
 {
-	int32_t dq = (int32_t)ai->y - (int32_t)aj->y, dr, dd, dg, q_span, sc;
-	int32_t sidi = (ai->y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
-	int32_t sidj = (aj->y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
-	if (dq <= 0 || dq > max_dist_x) return INT32_MIN;
-	dr = (int32_t)(ai->x - aj->x);
-	if (sidi == sidj && (dr == 0 || dq > max_dist_y)) return INT32_MIN;
+
+	uint64_t ai_x, ai_y, aj_x, aj_y;
+	ai_x = ai->x; ai_y = ai->y; aj_x = aj->x; aj_y = aj->y;
+
+#ifdef CHAIN_DEBUG
+	int32_t sc_vect = obj.comput_sc_vectorized_avx2_caller(ai_x, ai_y, aj_x, aj_y, aj->y>>32&0xff);
+#endif
+
+	//if (sc_vect == 0) return INT32_MIN;
+	//else 
+	//return sc_vect;
+
+	//fprintf(stderr, "%lld %lld %lld %lld \n", ai_x, ai_y, aj_x, aj_y); 
+	//fprintf(stderr, "%lld %lld %lld %f %f %d %d\n", max_dist_x, max_dist_y, bw, chn_pen_gap, chn_pen_skip, is_cdna, n_seg);
+	int32_t dq = (int32_t)ai_y - (int32_t)aj_y, dr, dd, dg, q_span, sc;
+	int32_t sidi = (ai_y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
+	int32_t sidj = (aj_y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
+	if (dq <= 0 || dq > max_dist_x) { 
+
+#ifdef CHAIN_DEBUG
+		if(INT32_MIN != sc_vect){
+		//fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect);
+			fprintf(stderr, "int-min exit: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect);
+		}
+#endif
+	  return INT32_MIN; 
+	}
+	dr = (int32_t)(ai_x - aj_x);
+	if (sidi == sidj && (dr == 0 || dq > max_dist_y)) { 
+
+#ifdef CHAIN_DEBUG
+		if(INT32_MIN != sc_vect){
+		//fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect);
+			fprintf(stderr, "int-min exit: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect);
+		}
+#endif
+		return INT32_MIN;
+	}
 	dd = dr > dq? dr - dq : dq - dr;
-	if (sidi == sidj && dd > bw) return INT32_MIN;
-	if (n_seg > 1 && !is_cdna && sidi == sidj && dr > max_dist_y) return INT32_MIN;
+	if (sidi == sidj && dd > bw) {
+
+#ifdef CHAIN_DEBUG
+		if(INT32_MIN != sc_vect){
+		//fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect);
+			fprintf(stderr, "int-min exit: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect);
+		}
+#endif
+		return INT32_MIN; 
+	}
+	if (n_seg > 1 && !is_cdna && sidi == sidj && dr > max_dist_y) { 
+
+#ifdef CHAIN_DEBUG
+		if(INT32_MIN != sc_vect){
+		//fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect);
+			fprintf(stderr, "int-min exit: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect);
+		}
+#endif
+		return INT32_MIN;
+	}
 	dg = dr < dq? dr : dq;
 	q_span = aj->y>>32&0xff;
 	sc = q_span < dg? q_span : dg;
@@ -114,11 +163,19 @@ static inline int32_t comput_sc(const mm128_t *ai, const mm128_t *aj, int32_t ma
 		float lin_pen, log_pen;
 		lin_pen = chn_pen_gap * (float)dd + chn_pen_skip * (float)dg;
 		log_pen = dd >= 1? mg_log2(dd + 1) : 0.0f; // mg_log2() only works for dd>=2
-		if (is_cdna) {
-			if (dr > dq) sc -= (int)(lin_pen < log_pen? lin_pen : log_pen); // deletion or jump between paired ends
+		if (is_cdna || sidi != sidj) {
+			if (sidi != sidj && dr == 0) ++sc; // possibly due to overlapping paired ends; give a minor bonus
+			else if (dr > dq || sidi != sidj) sc -= (int)(lin_pen < log_pen? lin_pen : log_pen); // deletion or jump between paired ends
 			else sc -= (int)(lin_pen + .5f * log_pen);
 		} else sc -= (int)(lin_pen + .5f * log_pen);
 	}
+#ifdef CHAIN_DEBUG
+
+	if(sc != sc_vect ){
+		//fprintf(stderr, "score mismatch %d -- %d", sc , sc_vect);
+		fprintf(stderr, "outer: %llu, %llu, %llu, %llu : %d -- %d\n", ai_x, ai_y, aj_x, aj_y, sc, sc_vect);
+	}
+#endif
 	return sc;
 }
 
@@ -133,10 +190,18 @@ static inline int32_t comput_sc(const mm128_t *ai, const mm128_t *aj, int32_t ma
 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_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;
+	///fprintf(stderr, "chaining called\n");
+
+
+
+#ifdef MANUAL_PROFILING
+	uint64_t align_start = __rdtsc();
+#endif
+
+	int32_t *f, *t, *v, *v_1, *p_1, n_u, n_v, mmax_f = 0;
 	int64_t *p, i, j, max_ii, st = 0, n_iter = 0;
 	uint64_t *u;
-
+	uint32_t* f_1;
 	if (_u) *_u = 0, *n_u_ = 0;
 	if (n == 0 || a == 0) {
 		kfree(km, a);
@@ -145,16 +210,56 @@ 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;
 	KMALLOC(km, p, n);
+	KMALLOC(km, p_1, n);
 	KMALLOC(km, f, n);
+	KMALLOC(km, f_1, n);
 	KMALLOC(km, v, n);
+	KMALLOC(km, v_1, n);
 	KCALLOC(km, t, n);
 
+//#ifdef PARALLEL_CHAINING
+if(enable_vect_dp_chaining){
+	// Parallel chaining data-structures
+	anchor_t* anchors = (anchor_t*)malloc(n* sizeof(anchor_t));
+	for (i = 0; i < n; ++i) {
+		uint64_t ri = a[i].x;
+		int32_t qi = (int32_t)a[i].y, q_span = a[i].y>>32&0xff; // NB: only 8 bits of span is used!!!
+		anchors[i].r = ri;
+		anchors[i].q = qi;
+		anchors[i].l = q_span;
+	}
+	num_bits_t *anchor_r, *anchor_q, *anchor_l;
+	create_SoA_Anchors_32_bit(anchors, n, anchor_r, anchor_q, anchor_l);
+	dp_chain obj(max_dist_x, max_dist_y, bw, max_skip, max_iter, min_cnt, min_sc, chn_pen_gap, chn_pen_skip, is_cdna, n_seg);
+
+#ifdef PARALLEL_CHAINING
+	obj.mm_dp_vectorized(n, &anchors[0], anchor_r, anchor_q, anchor_l, f_1, p_1, v_1, max_dist_x, max_dist_y, NULL, NULL);
+#endif
+	// -16 is due to extra padding at the start of arrays
+	anchor_r -= 16; anchor_q -= 16; anchor_l -= 16;
+	free(anchor_r); 
+	free(anchor_q); 
+	free(anchor_l);
+	free(anchors);
+	for(int i = 0; i < n; i++){
+#if 1
+			f[i] = f_1[i];
+			p[i] = p_1[i];
+			v[i] = v_1[i];
+#endif
+	}
+
+//
+} else {
+//#else
+
 	// fill the score and backtrack arrays
 	for (i = 0, max_ii = -1; i < n; ++i) {
 		int64_t max_j = -1, end_j;
 		int32_t max_f = a[i].y>>32&0xff, n_skip = 0;
 		while (st < i && (a[i].x>>32 != a[st].x>>32 || a[i].x > a[st].x + max_dist_x)) ++st;
 		if (i - st > max_iter) st = i - max_iter;
+		int my_cnt = 0;
 		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);
@@ -171,32 +276,62 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
 			if (p[j] >= 0) t[p[j]] = i;
 		}
 		end_j = j;
+		int debug_iter = 2057329;
+
 		if (max_ii < 0 || a[i].x - a[max_ii].x > (int64_t)max_dist_x) {
 			int32_t max = INT32_MIN;
 			max_ii = -1;
-			for (j = i - 1; j >= st; --j)
-				if (max < f[j]) max = f[j], max_ii = j;
+			for (j = i - 1; j >= st; --j) {
+				if (max < (int32_t)f[j]) max = f[j], max_ii = j;
+			}
 		}
 		if (max_ii >= 0 && max_ii < end_j) {
 			int32_t tmp;
 			tmp = comput_sc(&a[i], &a[max_ii], max_dist_x, max_dist_y, bw, chn_pen_gap, chn_pen_skip, is_cdna, n_seg);
-			if (tmp != INT32_MIN && max_f < tmp + f[max_ii])
+		//	if (i == debug_iter) fprintf(stderr, "mm2: endj: %d max_ii: %d max_f: %d tmp_score: %d \n", end_j, max_ii, max_f, tmp);	
+		
+
+			if (tmp != INT32_MIN && max_f < tmp + f[max_ii]){
 				max_f = tmp + f[max_ii], max_j = max_ii;
+			}
 		}
 		f[i] = max_f, p[i] = max_j;
 		v[i] = max_j >= 0 && v[max_j] > max_f? v[max_j] : max_f; // v[] keeps the peak score up to i; f[] is the score ending at i, not always the peak
+
 		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;
 	}
 
+}
+//#endif
+
+#ifdef CHAIN_DEBUG
+
+		for(int i = 0; i < n; i++){
+			if(f[i] != f_1[i] || p[i] != p_1[i] || v[i] !=v_1[i])
+			{
+				fprintf(stderr, "i:%d %d %d %d %d %d %d\n",i, f[i], f_1[i], p[i], p_1[i], v[i], v_1[i] );
+			}
+#if 0
+			f[i] = f_1[i];
+			p[i] = p_1[i];
+			v[i] = v_1[i];
+#endif
+		}
+#endif
 	u = mg_chain_backtrack(km, n, f, p, v, t, min_cnt, min_sc, &n_u, &n_v);
 	*n_u_ = n_u, *_u = u; // NB: note that u[] may not be sorted by score here
-	kfree(km, p); kfree(km, f); kfree(km, t);
+	kfree(km, p); kfree(km, p_1); kfree(km, f); kfree(km, f_1); kfree(km, t); kfree(km, v_1);
 	if (n_u == 0) {
 		kfree(km, a); kfree(km, v);
 		return 0;
 	}
+
+
+#ifdef MANUAL_PROFILING
+	dp_time += __rdtsc() - align_start;
+#endif
 	return compact_a(km, n_u, u, n_v, v, a);
 }
 
@@ -234,6 +369,11 @@ static inline int32_t comput_sc_simple(const mm128_t *ai, const mm128_t *aj, flo
 mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_skip, int cap_rmq_size, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip,
 					   int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km)
 {
+#ifdef MANUAL_PROFILING
+	uint64_t start = __rdtsc();
+#endif
+	uint64_t tim;
+	//fprintf(stderr, "rmq call \n");
 	int32_t *f,*t, *v, n_u, n_v, mmax_f = 0, max_rmq_size = 0;
 	int64_t *p, i, i0, st = 0, st_inner = 0, n_iter = 0;
 	uint64_t *u;
@@ -261,6 +401,9 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 		int32_t q_span = a[i].y>>32&0xff, max_f = q_span;
 		lc_elem_t s, *q, *r, lo, hi;
 		// add in-range anchors
+#ifdef MANUAL_PROFILING_RMQ
+	tim = __rdtsc();
+#endif
 		if (i0 < i && a[i0].x != a[i].x) {
 			int64_t j;
 			for (j = i0; j < i; ++j) {
@@ -275,7 +418,13 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 			}
 			i0 = i;
 		}
+#ifdef MANUAL_PROFILING_RMQ
+	rmq_t1 += __rdtsc() - tim;
+#endif
 		// get rid of active chains out of range
+#ifdef MANUAL_PROFILING_RMQ
+	tim = __rdtsc();
+#endif
 		while (st < i && (a[i].x>>32 != a[st].x>>32 || a[i].x > a[st].x + max_dist || krmq_size(head, root) > cap_rmq_size)) {
 			s.y = (int32_t)a[st].y, s.i = st;
 			if ((q = krmq_find(lc_elem, root, &s, 0)) != 0) {
@@ -284,6 +433,12 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 			}
 			++st;
 		}
+#ifdef MANUAL_PROFILING_RMQ
+	rmq_t2 += __rdtsc() - tim;
+#endif
+#ifdef MANUAL_PROFILING_RMQ
+	tim = __rdtsc();
+#endif
 		if (max_dist_inner > 0)  { // similar to the block above, but applied to the inner tree
 			while (st_inner < i && (a[i].x>>32 != a[st_inner].x>>32 || a[i].x > a[st_inner].x + max_dist_inner || krmq_size(head, root_inner) > cap_rmq_size)) {
 				s.y = (int32_t)a[st_inner].y, s.i = st_inner;
@@ -294,6 +449,9 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 				++st_inner;
 			}
 		}
+#ifdef MANUAL_PROFILING_RMQ
+	rmq_t3 += __rdtsc() - tim;
+#endif
 		// RMQ
 		lo.i = INT32_MAX, lo.y = (int32_t)a[i].y - max_dist;
 		hi.i = 0, hi.y = (int32_t)a[i].y;
@@ -313,6 +471,9 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 					krmq_itr_t(lc_elem) itr;
 					krmq_itr_find(lc_elem, root_inner, lo, &itr);
 					while ((q = krmq_at(&itr)) != 0) {
+#ifdef MANUAL_PROFILING_RMQ
+	tim = __rdtsc();
+#endif
 						if (q->y < (int32_t)a[i].y - max_dist_inner) break;
 						++n_rmq_iter;
 						j = q->i;
@@ -328,11 +489,15 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 							if (p[j] >= 0) t[p[j]] = i;
 						}
 						if (!krmq_itr_prev(lc_elem, &itr)) break;
+#ifdef MANUAL_PROFILING_RMQ
+	rmq_t4 += __rdtsc() - tim;
+#endif
 					}
 					n_iter += n_rmq_iter;
 				}
 			}
 		}
+
 		// set max
 		assert(max_j < 0 || (a[max_j].x < a[i].x && (int32_t)a[max_j].y < (int32_t)a[i].y));
 		f[i] = max_f, p[i] = max_j;
@@ -349,5 +514,8 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 		kfree(km, a); kfree(km, v);
 		return 0;
 	}
+#ifdef MANUAL_PROFILING
+	rmq_time += __rdtsc() -  start;
+#endif
 	return compact_a(km, n_u, u, n_v, v, a);
 }

main.c

@@ -6,8 +6,77 @@
 #include "minimap.h"
 #include "mmpriv.h"
 #include "ketopt.h"
+#include <x86intrin.h>
+#include <immintrin.h>
+#include <sys/time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <string>
+#include <map>
+#include <errno.h>
+#include "bseq.h"
+#include "minimap.h"
+#include "mmpriv.h"
+#include "ketopt.h"
 
-#define MM_VERSION "2.20-r1061"
+//#include "profile.h"
+#include <stdint.h>
+#include <unistd.h>
+#include <x86intrin.h>
+
+using namespace std;
+uint64_t avg;
+uint64_t minimizer_lookup_time, alignment_time, dp_time, rmq_time, rmq_t1, rmq_t2, rmq_t3, rmq_t4;
+
+bool enable_vect_dp_chaining = false;
+
+#ifdef LISA_HASH
+#include "lisa_hash.h"
+lisa_hash<uint64_t, uint64_t> *lh;
+#endif
+
+// New memory allocation approach for alignment optimizations
+//
+void *km1;
+uint64_t km_size = 500000000; // 500 MB
+int km_top;
+/*
+void *kcalloc_(void* km, int count, int size)
+{
+    assert(count*size < km_size);
+    km_top += count*size + 1024;
+    memset(km, 0, count * size);
+        
+    // printf("km_top: %d\n", km_top);
+    return km;
+}
+
+void *kmalloc_(void* km, int count) {
+    if(km_top + count >= km_size)
+        printf("count: %d\n", count);
+    assert(km_top + count < km_size);
+    void *mem = (void*) ((int8_t*) km + km_top);
+    km_top += count + 1024;
+    // printf("km_top: %d\n", km_top);
+    return mem;
+}
+
+void kfree_all() { km_top = 0;}
+  */  
+
+// Memory for alignment end
+
+
+#ifndef __rdtsc
+#ifdef _rdtsc
+#define __rdtsc _rdtsc
+#else
+#define __rdtsc __builtin_ia32_rdtsc
+#endif
+#endif
+
+#define MM_VERSION "2.22-r1101"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -72,6 +141,8 @@ static ko_longopt_t long_options[] = {
 	{ "alt-drop",       ko_required_argument, 345 },
 	{ "mask-len",       ko_required_argument, 346 },
 	{ "rmq",            ko_optional_argument, 347 },
+	{ "qstrand",        ko_no_argument,       348 },
+	{ "cap-kalloc",     ko_required_argument, 349 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
 	{ "version",        ko_no_argument,       'V' },
@@ -100,7 +171,7 @@ static inline int64_t mm_parse_num(const char *str)
 	return mm_parse_num2(str, 0);
 }
 
-static inline void yes_or_no(mm_mapopt_t *opt, int flag, int long_idx, const char *arg, int yes_to_set)
+static inline void yes_or_no(mm_mapopt_t *opt, int64_t flag, int long_idx, const char *arg, int yes_to_set)
 {
 	if (yes_to_set) {
 		if (strcmp(arg, "yes") == 0 || strcmp(arg, "y") == 0) opt->flag |= flag;
@@ -115,6 +186,12 @@ static inline void yes_or_no(mm_mapopt_t *opt, int flag, int long_idx, const cha
 
 int main(int argc, char *argv[])
 {
+	// Memory allocation for alignment optimizations
+    //km1 = calloc(km_size, 1); // 10 MB init contg. alloc
+#ifdef PARALLEL_CHAINING
+	enable_vect_dp_chaining = true;
+#endif
+
 	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:";
 	ketopt_t o = KETOPT_INIT;
 	mm_mapopt_t opt;
@@ -129,9 +206,11 @@ int main(int argc, char *argv[])
 	liftrlimit();
 	mm_realtime0 = realtime();
 	mm_set_opt(0, &ipt, &opt);
+	string preset_arg = "";
 
 	while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) { // test command line options and apply option -x/preset first
 		if (c == 'x') {
+			preset_arg +=  (string) o.arg;
 			if (mm_set_opt(o.arg, &ipt, &opt) < 0) {
 				fprintf(stderr, "[ERROR] unknown preset '%s'\n", o.arg);
 				return 1;
@@ -225,6 +304,8 @@ int main(int argc, char *argv[])
 		else if (c == 344) alt_list = o.arg; // --alt
 		else if (c == 345) opt.alt_drop = atof(o.arg); // --alt-drop
 		else if (c == 346) opt.mask_len = mm_parse_num(o.arg); // --mask-len
+		else if (c == 348) opt.flag |= MM_F_QSTRAND | MM_F_NO_INV; // --qstrand
+		else if (c == 349) opt.cap_kalloc = mm_parse_num(o.arg); // --cap-kalloc
 		else if (c == 330) {
 			fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n");
 		} else if (c == 314) { // --frag
@@ -326,7 +407,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -N INT       retain at most INT secondary alignments [%d]\n", opt.best_n);
 		fprintf(fp_help, "  Alignment:\n");
 		fprintf(fp_help, "    -A INT       matching score [%d]\n", opt.a);
-		fprintf(fp_help, "    -B INT       mismatch penalty [%d]\n", opt.b);
+		fprintf(fp_help, "    -B INT       mismatch penalty (larger value for lower divergence) [%d]\n", opt.b);
 		fprintf(fp_help, "    -O INT[,INT] gap open penalty [%d,%d]\n", opt.q, opt.q2);
 		fprintf(fp_help, "    -E INT[,INT] gap extension penalty; a k-long gap costs min{O1+k*E1,O2+k*E2} [%d,%d]\n", opt.e, opt.e2);
 		fprintf(fp_help, "    -z INT[,INT] Z-drop score and inversion Z-drop score [%d,%d]\n", opt.zdrop, opt.zdrop_inv);
@@ -362,6 +443,7 @@ int main(int argc, char *argv[])
 		fprintf(stderr, "[ERROR] incorrect input: in the sr mode, please specify no more than two query files.\n");
 		return 1;
 	}
+	preset_arg = (string)argv[o.ind] + "_" + preset_arg + "_minimizers_key_value_sorted";
 	idx_rdr = mm_idx_reader_open(argv[o.ind], &ipt, fnw);
 	if (idx_rdr == 0) {
 		fprintf(stderr, "[ERROR] failed to open file '%s': %s\n", argv[o.ind], strerror(errno));
@@ -404,10 +486,26 @@ int main(int argc, char *argv[])
 					__func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), mi->n_seq);
 		if (argc != o.ind + 1) mm_mapopt_update(&opt, mi);
 		if (mm_verbose >= 3) mm_idx_stat(mi);
+#ifdef LISA_INDEX
+		mm_idx_dump_hash(preset_arg.c_str(), mi); 
+#endif	
 		if (junc_bed) mm_idx_bed_read(mi, junc_bed, 1);
 		if (alt_list) mm_idx_alt_read(mi, alt_list);
-		if (argc - (o.ind + 1) == 0) continue; // no query files
+		if (argc - (o.ind + 1) == 0) {
+			mm_idx_destroy(mi);
+			continue; // no query files
+		}
 		ret = 0;
+#ifdef LISA_HASH
+	fprintf(stderr, "Using LISA_HASH..\n");
+	mm_idx_destroy_mm_hash(mi);
+	char* prefix;
+	lh = new lisa_hash<uint64_t, uint64_t>(preset_arg, prefix);
+	fprintf(stderr, "Loading done.\n");
+//	total_time =  __rdtsc();
+//	fprintf(stderr, "\nIndexing Real time: %.3f sec;\n", realtime() - mapping_time);
+#endif
+	mm_realtime0 = realtime();
 		if (!(opt.flag & MM_F_FRAG_MODE)) {
 			for (i = o.ind + 1; i < argc; ++i) {
 				ret = mm_map_file(mi, argv[i], &opt, n_threads);
@@ -416,12 +514,17 @@ int main(int argc, char *argv[])
 		} else {
 			ret = mm_map_file_frag(mi, argc - (o.ind + 1), (const char**)&argv[o.ind + 1], &opt, n_threads);
 		}
-		mm_idx_destroy(mi);
+		//mm_idx_destroy(mi);
 		if (ret < 0) {
 			fprintf(stderr, "ERROR: failed to map the query file\n");
 			exit(EXIT_FAILURE);
 		}
 	}
+#ifdef LISA_HASH
+		mm_idx_destroy_seq(mi);
+#else
+		mm_idx_destroy(mi);
+#endif
 	n_parts = idx_rdr->n_parts;
 	mm_idx_reader_close(idx_rdr);
 
@@ -440,5 +543,10 @@ int main(int argc, char *argv[])
 			fprintf(stderr, " %s", argv[i]);
 		fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, realtime() - mm_realtime0, cputime(), peakrss() / 1024.0 / 1024.0 / 1024.0);
 	}
+
+	fprintf(stderr, "minimizer-lookup: %lld dp: %lld rmq: %lld rmq_t1: %lld rmq_t2: %lld rmq_t3: %lld rmq_t4: %lld alignment: %lld %lld\n", minimizer_lookup_time, dp_time, rmq_time, rmq_t1, rmq_t2, rmq_t3, rmq_t4, alignment_time, avg);
+#ifdef LISA_HASH
+	delete lh;
+#endif	
 	return 0;
 }

map.c

@@ -9,6 +9,12 @@
 #include "mmpriv.h"
 #include "bseq.h"
 #include "khash.h"
+#include <x86intrin.h>
+
+#ifdef MANUAL_PROFILING
+extern uint64_t minimizer_lookup_time;
+extern uint64_t rmq_time;
+#endif
 
 struct mm_tbuf_s {
 	void *km;
@@ -173,6 +179,9 @@ static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max
 static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, const mm128_v *mv, int qlen, int64_t *n_a, int *rep_len,
 								  int *n_mini_pos, uint64_t **mini_pos)
 {
+#ifdef MANUAL_PROFILING
+	uint64_t lookup_start = __rdtsc();
+#endif
 	int i, n_m;
 	mm_seed_t *m;
 	mm128_t *a;
@@ -190,9 +199,13 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ,
 			if ((r[k]&1) == (q->q_pos&1)) { // forward strand
 				p->x = (r[k]&0xffffffff00000000ULL) | rpos;
 				p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
-			} else { // reverse strand
+			} else if (!(opt->flag & MM_F_QSTRAND)) { // reverse strand and not in the query-strand mode
 				p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | rpos;
 				p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
+			} else { // reverse strand; query-strand
+				int32_t len = mi->seq[r[k]>>32].len;
+				p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | (len - (rpos + 1 - q->q_span) - 1); // coordinate only accurate for non-HPC seeds
+				p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
 			}
 			p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
 			if (q->is_tandem) p->y |= MM_SEED_TANDEM;
@@ -201,6 +214,9 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ,
 	}
 	kfree(km, m);
 	radix_sort_128x(a, a + (*n_a));
+#ifdef MANUAL_PROFILING
+	minimizer_lookup_time += __rdtsc() - lookup_start;
+#endif
 	return a;
 }
 
@@ -272,7 +288,12 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	if (opt->flag & MM_F_RMQ) {
 		a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
 						  opt->chain_gap_scale * 0.01 * mi->k, 0.0f, n_a, a, &n_regs0, &u, b->km);
+//		a = mg_lchain_dp(opt->max_gap, opt->rmq_inner_dist, opt->bw, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
+//						 	opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
+
+
 	} else {
+		//fprintf(stderr, "dp call - n_a = %lld\n", n_a);
 		a = mg_lchain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score,
 						 opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 	}
@@ -280,12 +301,23 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	if (opt->bw_long > opt->bw && (opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN)) == 0 && n_segs == 1 && n_regs0 > 1) { // re-chain/long-join for long sequences
 		int32_t st = (int32_t)a[0].y, en = (int32_t)a[(int32_t)u[0] - 1].y;
 		if (qlen_sum - (en - st) > opt->rmq_rescue_size || en - st > qlen_sum * opt->rmq_rescue_ratio) {
+#ifdef MANUAL_PROFILING
+//			uint64_t tim = __rdtsc();
+#endif
+	//		fprintf(stderr, "pre: rmq rechain call - n_a = %lld n_regs = %lld\n",n_a, n_regs0);
 			int32_t i;
+			int64_t prev_n_a = n_a;
 			for (i = 0, n_a = 0; i < n_regs0; ++i) n_a += (int32_t)u[i];
 			kfree(b->km, u);
 			radix_sort_128x(a, a + n_a);
+//			fprintf(stderr, "post: rmq rechain call - prev_n_a = %lld n_a = %lld n_regs = %lld\n",prev_n_a, n_a, n_regs0);
+//			a = mg_lchain_dp(opt->max_gap, opt->rmq_inner_dist, opt->bw_long, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
+//						 	opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 			a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw_long, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
 							  opt->chain_gap_scale * 0.01 * mi->k, 0.0f, n_a, a, &n_regs0, &u, b->km);
+#ifdef MANUAL_PROFILING
+//			rmq_time += __rdtsc() - tim;
+#endif
 		}
 	} else if (opt->max_occ > opt->mid_occ && rep_len > 0 && !(opt->flag & MM_F_RMQ)) { // re-chain, mostly for short reads
 		int rechain = 0;
@@ -314,7 +346,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	b->frag_gap = max_chain_gap_ref;
 	b->rep_len = rep_len;
 
-	regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a);
+	regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a, !!(opt->flag&MM_F_QSTRAND));
 	if (mi->n_alt) {
 		mm_mark_alt(mi, n_regs0, regs0);
 		mm_hit_sort(b->km, &n_regs0, regs0, opt->alt_drop); // this step can be merged into mm_gen_regs(); will do if this shows up in profile
@@ -327,10 +359,12 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 						i == regs0[j].as? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));
 
 	chain_post(opt, max_chain_gap_ref, mi, b->km, qlen_sum, n_segs, qlens, &n_regs0, regs0, a);
-	if (!is_sr) mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos);
+	if (!is_sr && !(opt->flag&MM_F_QSTRAND))
+		mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos);
 
 	if (n_segs == 1) { // uni-segment
 		regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], &n_regs0, regs0, a);
+		regs0 = (mm_reg1_t*)realloc(regs0, sizeof(*regs0) * n_regs0);
 		mm_set_mapq(b->km, n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr);
 		n_regs[0] = n_regs0, regs[0] = regs0;
 	} else { // multi-segment
@@ -357,7 +391,9 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
 			fprintf(stderr, "QM\t%s\t%d\tcap=%ld,nCore=%ld,largest=%ld\n", qname, qlen_sum, kmst.capacity, kmst.n_cores, kmst.largest);
 		assert(kmst.n_blocks == kmst.n_cores); // otherwise, there is a memory leak
-		if (kmst.largest > 1U<<28) {
+		if (kmst.largest > 1U<<28 || (opt->cap_kalloc > 0 && kmst.capacity > opt->cap_kalloc)) {
+			if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
+				fprintf(stderr, "[W::%s] reset thread-local memory after read %s\n", __func__, qname);
 			km_destroy(b->km);
 			b->km = km_init();
 		}
@@ -402,10 +438,13 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
     step_t *s = (step_t*)_data;
 	int qlens[MM_MAX_SEG], j, off = s->seg_off[i], pe_ori = s->p->opt->pe_ori;
 	const char *qseqs[MM_MAX_SEG];
+	double t = 0.0;
 	mm_tbuf_t *b = s->buf[tid];
 	assert(s->n_seg[i] <= MM_MAX_SEG);
-	if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
+	if (mm_dbg_flag & MM_DBG_PRINT_QNAME) {
 		fprintf(stderr, "QR\t%s\t%d\t%d\n", s->seq[off].name, tid, s->seq[off].l_seq);
+		t = realtime();
+	}
 	for (j = 0; j < s->n_seg[i]; ++j) {
 		if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1))))
 			mm_revcomp_bseq(&s->seq[off + j]);
@@ -437,6 +476,8 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
 				r->rev = !r->rev;
 			}
 		}
+	if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
+		fprintf(stderr, "QT\t%s\t%d\t%.6f\n", s->seq[off].name, tid, realtime() - t);
 }
 
 static void merge_hits(step_t *s)
@@ -481,6 +522,14 @@ static void merge_hits(step_t *s)
 					}
 				}
 			}
+			if (!(opt->flag&MM_F_SR) && s->seq[k].l_seq >= opt->rank_min_len)
+				mm_update_dp_max(s->seq[k].l_seq, s->n_reg[k], s->reg[k], opt->rank_frac, opt->a, opt->b);
+			for (j = 0; j < s->n_reg[k]; ++j) {
+				mm_reg1_t *r = &s->reg[k][j];
+				if (r->p) r->p->dp_max2 = 0; // reset ->dp_max2 as mm_set_parent() doesn't clear it; necessary with mm_update_dp_max()
+				r->subsc = 0; // this may not be necessary
+				r->n_sub = 0; // n_sub will be an underestimate as we don't see all the chains now, but it can't be accurate anyway
+			}
 			mm_hit_sort(km, &s->n_reg[k], s->reg[k], opt->alt_drop);
 			mm_set_parent(km, opt->mask_level, opt->mask_len, s->n_reg[k], s->reg[k], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
 			if (!(opt->flag & MM_F_ALL_CHAINS)) {
@@ -542,6 +591,7 @@ static void *worker_pipeline(void *shared, int step, void *in)
 		if ((p->opt->flag & MM_F_OUT_CS) && !(mm_dbg_flag & MM_DBG_NO_KALLOC)) km = km_init();
 		for (k = 0; k < s->n_frag; ++k) {
 			int seg_st = s->seg_off[k], seg_en = s->seg_off[k] + s->n_seg[k];
+#ifndef DISABLE_OUTPUT
 			for (i = seg_st; i < seg_en; ++i) {
 				mm_bseq1_t *t = &s->seq[i];
 				if (p->opt->split_prefix && p->n_parts == 0) { // then write to temporary files
@@ -576,6 +626,7 @@ static void *worker_pipeline(void *shared, int step, void *in)
 					mm_err_puts(p->str.s);
 				}
 			}
+#endif
 			for (i = seg_st; i < seg_en; ++i) {
 				for (j = 0; j < s->n_reg[i]; ++j) free(s->reg[i][j].p);
 				free(s->reg[i]);

minimap.h

@@ -36,7 +36,9 @@
 #define MM_F_NO_END_FLT    0x10000000
 #define MM_F_HARD_MLEVEL   0x20000000
 #define MM_F_SAM_HIT_ONLY  0x40000000
-#define MM_F_RMQ           0x80000000LL
+#define MM_F_RMQ           (0x80000000LL)
+#define MM_F_QSTRAND       (0x100000000LL)
+#define MM_F_NO_INV        (0x200000000LL)
 
 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2
@@ -46,6 +48,18 @@
 
 #define MM_MAX_SEG       255
 
+#define MM_CIGAR_MATCH      0
+#define MM_CIGAR_INS        1
+#define MM_CIGAR_DEL        2
+#define MM_CIGAR_N_SKIP     3
+#define MM_CIGAR_SOFTCLIP   4
+#define MM_CIGAR_HARDCLIP   5
+#define MM_CIGAR_PADDING    6
+#define MM_CIGAR_EQ_MATCH   7
+#define MM_CIGAR_X_MISMATCH 8
+
+#define MM_CIGAR_STR  "MIDNSHP=XB"
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -143,6 +157,9 @@ typedef struct {
 	int anchor_ext_len, anchor_ext_shift;
 	float max_clip_ratio; // drop an alignment if BOTH ends are clipped above this ratio
 
+	int rank_min_len;
+	float rank_frac;
+
 	int pe_ori, pe_bonus;
 
 	float mid_occ_frac;  // only used by mm_mapopt_update(); see below
@@ -151,6 +168,7 @@ typedef struct {
 	int32_t max_occ, max_max_occ, occ_dist;
 	int64_t mini_batch_size; // size of a batch of query bases to process in parallel
 	int64_t max_sw_mat;
+	int64_t cap_kalloc;
 
 	const char *split_prefix;
 } mm_mapopt_t;
@@ -267,6 +285,13 @@ mm_idx_t *mm_idx_load(FILE *fp);
  */
 void mm_idx_dump(FILE *fp, const mm_idx_t *mi);
 
+/**
+ * Store hash table from minimap2 index into a file
+ * @param f_name     File name for output file
+ * @param mi 	     minimap2 index 	
+ */
+void mm_idx_dump_hash(const char* f_name, const mm_idx_t *mi);
+
 /**
  * Create an index from strings in memory
  *
@@ -295,6 +320,19 @@ void mm_idx_stat(const mm_idx_t *idx);
  * @param r          minimap2 index
  */
 void mm_idx_destroy(mm_idx_t *mi);
+/**
+ * Destroy/deallocate an hash table index
+ *
+ * @param r          minimap2 index
+ */
+void mm_idx_destroy_mm_hash(mm_idx_t *mi);
+
+/**
+ * Destroy/deallocate target sequences
+ *
+ * @param r          minimap2 index
+ */
+void mm_idx_destroy_seq(mm_idx_t *mi);
 
 /**
  * Initialize a thread-local buffer for mapping

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "27 May 2021" "minimap2-2.20 (r1061)" "Bioinformatics tools"
+.TH minimap2 1 "7 August 2021" "minimap2-2.22 (r1101)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -423,6 +423,11 @@ alignment.
 Skip alignment if the DP matrix size is above
 .IR NUM .
 Set 0 to disable [100m].
+.TP
+.BI --cap-kalloc \ NUM
+Free thread-local kalloc memory reservoir if after the alignment the size of the reservoir above
+.IR NUM .
+Set 0 to disable [0].
 .SS Input/output options
 .TP 10
 .B -a
@@ -573,7 +578,7 @@ Up to 20% sequence divergence.
 .B splice
 Long-read spliced alignment
 .RB ( -k15
-.B -w5 --splice -g2k -G200k -A1 -B2 -O2,32 -E1,0 -C9 -z200 -ub --junc-bonus=9 --cap-sw-mem=0
+.B -w5 --splice -g2k -G200k -A1 -B2 -O2,32 -E1,0 -b0 -C9 -z200 -ub --junc-bonus=9 --cap-sw-mem=0
 .BR --splice-flank=yes ).
 In the splice mode, 1) long deletions are taken as introns and represented as
 the
@@ -592,8 +597,8 @@ Long-read splice alignment for PacBio CCS reads
 .B sr
 Short single-end reads without splicing
 .RB ( -k21
-.B -w11 --sr --frag=yes -A2 -B8 -O12,32 -E2,1 -r50 -p.5 -N20 -f1000,5000 -n2 -m20
-.B -s40 -g200 -2K50m --heap-sort=yes
+.B -w11 --sr --frag=yes -A2 -B8 -O12,32 -E2,1 -b0 -r100 -p.5 -N20 -f1000,5000 -n2 -m20
+.B -s40 -g100 -2K50m --heap-sort=yes
 .BR --secondary=no ).
 .TP
 .B ava-pb

misc/mmphase.js

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

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = '2.20-r1061';
+var paftools_version = '2.22-r1101';
 
 /*****************************
  ***** Library functions *****
@@ -977,7 +977,7 @@ function paf_stat(args)
 	var re = /(\d+)([MIDSHNX=])/g;
 
 	var lineno = 0, n_pri = 0, n_2nd = 0, n_seq = 0, n_cigar_64k = 0, l_tot = 0, l_cov = 0;
-	var n_gap = [[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]];
+	var n_gap = [[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]], n_sub = 0;
 
 	function cov_len(regs)
 	{
@@ -999,7 +999,7 @@ function paf_stat(args)
 		if (line.charAt(0) != '@') {
 			var t = line.split("\t", 12);
 			var m, rs, cigar = null, is_pri = false, is_sam = false, is_rev = false, tname = null;
-			var atlen = null, aqlen, qs, qe, mapq, ori_qlen, NM = null;
+			var atlen = null, aqlen, qs, qe, mapq, ori_qlen, NM = null, nn = 0;
 			if (t.length < 2) continue;
 			if (t[4] == '+' || t[4] == '-' || t[4] == '*') { // PAF
 				if (t[4] == '*') continue; // unmapped
@@ -1009,6 +1009,8 @@ function paf_stat(args)
 				}
 				if ((m = /\tNM:i:(\d+)/.exec(line)) != null)
 					NM = parseInt(m[1]);
+				if ((m = /\tnn:i:(\d+)/.exec(line)) != null)
+					nn = parseInt(m[1]);
 				if ((m = /\tcg:Z:(\S+)/.exec(line)) != null)
 					cigar = m[1];
 				if (cigar == null) {
@@ -1032,6 +1034,8 @@ function paf_stat(args)
 				}
 				if ((m = /\tNM:i:(\d+)/.exec(line)) != null)
 					NM = parseInt(m[1]);
+				if ((m = /\tnn:i:(\d+)/.exec(line)) != null)
+					nn = parseInt(m[1]);
 				cigar = t[5];
 				tname = t[2];
 				rs = parseInt(t[3]) - 1;
@@ -1078,6 +1082,12 @@ function paf_stat(args)
 					clip[M == 0? 0 : 1] = l;
 				}
 			}
+			if (NM != null) {
+				var tmp = NM - n_gap_all - nn;
+				if (tmp < 0 && nn == 0) warn("WARNING: NM is smaller than the number of gaps at line " + lineno + ": NM=" + NM + ", nn=" + nn + ", G=" + n_gap_all);
+				if (tmp < 0) tmp = 0;
+				n_sub += tmp;
+			}
 			if (n_cigar > 65535) ++n_cigar_64k;
 			if (ql + sclip != aqlen)
 				warn("WARNING: aligned query length is inconsistent with CIGAR at line " + lineno + " (" + (ql+sclip) + " != " + aqlen + ")");
@@ -1112,6 +1122,7 @@ function paf_stat(args)
 		print("Number of primary alignments with >65535 CIGAR operations: " + n_cigar_64k);
 		print("Number of bases in mapped sequences: " + l_tot);
 		print("Number of mapped bases: " + l_cov);
+		print("Number of substitutions: " + n_sub);
 		print("Number of insertions in [0,50): " + n_gap[0][0]);
 		print("Number of insertions in [50,100): " + n_gap[0][1]);
 		print("Number of insertions in [100,300): " + n_gap[0][2]);
@@ -1419,7 +1430,7 @@ function paf_view(args)
 
 	var s_ref = new Bytes(), s_qry = new Bytes(), s_mid = new Bytes(); // these are used to show padded alignment
 	var re_cs = /([:=\-\+\*])(\d+|[A-Za-z]+)/g;
-	var re_cg = /(\d+)([MIDNSH])/g;
+	var re_cg = /(\d+)([MIDNSHP=X])/g;
 
 	var buf = new Bytes();
 	var file = args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]);
@@ -1475,8 +1486,14 @@ function paf_view(args)
 				warn("WARNING: converting to BLAST-like alignment requires the 'cs' tag, which is absent on line " + lineno);
 				continue;
 			}
+			var n_mm = 0, n_oi = 0, n_od = 0, n_ei = 0, n_ed = 0;
+			while ((m = re_cs.exec(cs)) != null) {
+				if (m[1] == '*') ++n_mm;
+				else if (m[1] == '+') ++n_oi, n_ei += m[2].length;
+				else if (m[1] == '-') ++n_od, n_ed += m[2].length;
+			}
 			line = line.replace(/\tc[sg]:Z:\S+/g, ""); // get rid of cs or cg tags
-			print('>' + line);
+			print('>' + line + "\tmm:i:"+n_mm + "\toi:i:"+n_oi + "\tei:i:"+n_ei + "\tod:i:"+n_od + "\ted:i:"+n_ed);
 			var rs = parseInt(t[7]), qs = t[4] == '+'? parseInt(t[2]) : parseInt(t[3]);
 			var n_blocks = 0;
 			while ((m = re_cs.exec(cs)) != null) {
@@ -1899,7 +1916,7 @@ function paf_splice2bed(args)
 		a.length = 0;
 	}
 
-	var re = /(\d+)([MIDNSH])/g;
+	var re = /(\d+)([MIDNSHP=X])/g;
 	var c, fmt = "bed", fn_name_conv = null, keep_multi = false;
 	while ((c = getopt(args, "f:n:m")) != null) {
 		if (c == 'f') fmt = getopt.arg;
@@ -2369,7 +2386,7 @@ 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+)([MIDNSHX=])/g;
+	var re_cigar = /(\d+)([MIDNSHP=X])/g;
 	while (file.readline(buf) >= 0) {
 		var m, t = buf.toString().split("\t");
 		var ctg_name = null, cigar = null, pos = null, qname = t[0];
@@ -2930,6 +2947,120 @@ function paf_vcfsel(args)
 	buf.destroy();
 }
 
+function paf_pafcmp(args)
+{
+	var c, opt = { min_len:5000, min_mapq:10, min_ovlp:0.5 };
+	while ((c = getopt(args, "q:")) != null) {
+		if (c == 'q') opt.min_mapq = parseInt(opt.arg);
+	}
+
+	var buf = new Bytes();
+	if (args.length - getopt.ind < 2) {
+		print("Usage: paftools.js pafcmp [options] <base.paf> <test.paf>");
+		print("Options:");
+		print("  -q INT    min mapping quality [" + opt.min_mapq + "]");
+		return 1;
+	}
+
+	var eval = { n_base:0, n_test:0, n_out_high:0, n_out_low:0, n_hit:0, n_wrong:0, n_miss:0 };
+
+	function process_base(base, a) {
+		if (a.length != 1) return;
+		for (var i = 1; i < 4; ++i)
+			a[0][i] = parseInt(a[0][i]);
+		for (var i = 6; i < 12; ++i)
+			a[0][i] = parseInt(a[0][i]);
+		if (a[0][1] < opt.min_len) return;
+		if (a[0][11] >= opt.min_mapq) ++eval.n_base;
+		base[a[0][0]] = [a[0][5], a[0][7], a[0][8], a[0][11], 0, 0];
+	}
+
+	var file = new File(args[getopt.ind]);
+	warn("Reading " + args[getopt.ind] + "...");
+	var a = [], base = {};
+	while (file.readline(buf) >= 0) {
+		var line = buf.toString();
+		var t = line.split("\t");
+		if (/\ttp:A:S/.test(line)) continue;
+		if (a.length > 0 && a[0][0] != t[0]) {
+			process_base(base, a);
+			a = [];
+		}
+		a.push(t);
+	}
+	process_base(base, a);
+	file.close();
+
+	function process_test(base, a) {
+		for (var i = 1; i < 4; ++i)
+			a[0][i] = parseInt(a[0][i]);
+		for (var i = 6; i < 12; ++i)
+			a[0][i] = parseInt(a[0][i]);
+		if (a[0][1] < opt.min_len) return;
+		if (a[0][11] >= opt.min_mapq) ++eval.n_test;
+		var c = [a[0][5], a[0][7], a[0][8], a[0][11]];
+		if (base[a[0][0]] == null) {
+			if (c[3] >= opt.min_mapq) ++opt.n_out_high;
+			else ++opt.n_out_low;
+		} else {
+			var b = base[a[0][0]];
+			var inter = 0, union = (b[2] - b[1]) + (c[2] - c[1]);
+			if (b[0] == c[0]) { // same chr
+				if (b[1] < c[1]) {
+					if (b[2] > c[1])
+						inter = b[2] - c[1], union = c[2] - b[1];
+				} else { // c[1] < b[1]
+					if (c[2] > b[1])
+						inter = c[2] - b[1], union = b[2] - c[1];
+				}
+			}
+			if (inter >= union * opt.min_ovlp) {
+				if (b[3] >= opt.min_mapq) ++eval.n_hit;
+				++b[4];
+			} else {
+				if (b[3] >= opt.min_mapq) {
+					print("W", a[0][0], b.slice(0, 4).join("\t"), c.join("\t"));
+					++eval.n_wrong;
+				}
+				++b[5];
+			}
+		}
+	}
+
+	file = new File(args[getopt.ind+1]);
+	warn("Reading " + args[getopt.ind+1] + "...");
+	a = [];
+	while (file.readline(buf) >= 0) {
+		var line = buf.toString();
+		var t = line.split("\t");
+		if (/\ttp:A:S/.test(line)) continue;
+		if (a.length > 0 && a[0][0] != t[0]) {
+			process_test(base, a);
+			a = [];
+		}
+		a.push(t);
+	}
+	process_test(base, a);
+	file.close();
+
+	for (var r in base) {
+		var b = base[r];
+		if (b[3] >= opt.min_mapq && b[4] == 0 && b[5] == 0) {
+			++eval.n_miss;
+			print("M", r, b.slice(0, 4).join("\t"));
+		}
+	}
+
+	print("X", eval.n_base     + " base alignments with mapQ>=" + opt.min_mapq);
+//	print("X", eval.n_test     + " test alignments with mapQ>=" + opt.min_mapq);
+	print("X", eval.n_hit      + " base alignments correctly mapped by test");
+	print("X", eval.n_wrong    + " wrong test alignment");
+	print("X", eval.n_miss     + " base alignments missing");
+	print("X", eval.n_out_high + " additional test alignments with mapQ>=" + opt.min_mapq);
+
+	buf.destroy();
+}
+
 /*************************
  ***** main function *****
  *************************/
@@ -2957,6 +3088,7 @@ function main(args)
 		print("  version    print paftools.js version");
 		print("");
 		print("  mapeval    evaluate mapping accuracy using mason2/PBSIM-simulated FASTQ");
+		print("  pafcmp     compare two PAF files");
 		print("  mason2fq   convert mason2-simulated SAM to FASTQ");
 		print("  pbsim2fq   convert PBSIM-simulated MAF to FASTQ");
 		print("  junceval   evaluate splice junction consistency with known annotations");
@@ -2978,6 +3110,7 @@ function main(args)
 	else if (cmd == 'vcfpair') paf_vcfpair(args);
 	else if (cmd == 'call') paf_call(args);
 	else if (cmd == 'mapeval') paf_mapeval(args);
+	else if (cmd == 'pafcmp') paf_pafcmp(args);
 	else if (cmd == 'bedcov') paf_bedcov(args);
 	else if (cmd == 'mason2fq') paf_mason2fq(args);
 	else if (cmd == 'pbsim2fq') paf_pbsim2fq(args);

mmpriv.h

@@ -62,17 +62,20 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i
 
 mm_seed_t *mm_collect_matches(void *km, int *_n_m, int qlen, int max_occ, int max_max_occ, int dist, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos);
 
+double mm_event_identity(const mm_reg1_t *r);
 int mm_write_sam_hdr(const mm_idx_t *mi, const char *rg, const char *ver, int argc, char *argv[]);
-void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag);
-void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag, int rep_len);
+void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int64_t opt_flag);
+void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int64_t opt_flag, int rep_len);
 void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int n_regs, const mm_reg1_t *regs);
-void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regs, const mm_reg1_t *const* regs, void *km, int opt_flag);
-void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag, int rep_len);
+void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regs, const mm_reg1_t *const* regs, void *km, int64_t opt_flag);
+void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int64_t opt_flag, int rep_len);
 
 void mm_idxopt_init(mm_idxopt_t *opt);
 const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f);
+int mm_idx_getseq2(const mm_idx_t *mi, int is_rev, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
 mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
+mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a, int is_qstrand);
 
 mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float gap_scale,
 					 int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
@@ -81,9 +84,8 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
 mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_skip, int cap_rmq_size, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip,
 					   int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
 
-mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a);
 void mm_mark_alt(const mm_idx_t *mi, int n, mm_reg1_t *r);
-void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a);
+void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a, int is_qstrand);
 void mm_sync_regs(void *km, int n_regs, mm_reg1_t *regs);
 int mm_squeeze_a(void *km, int n_regs, mm_reg1_t *regs, mm128_t *a);
 int mm_set_sam_pri(int n, mm_reg1_t *r);
@@ -93,6 +95,7 @@ void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int
 void mm_filter_regs(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs);
 void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r, float alt_diff_frac);
 void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr);
+void mm_update_dp_max(int qlen, int n_regs, mm_reg1_t *regs, float frac, int a, int b);
 
 void mm_est_err(const mm_idx_t *mi, int qlen, int n_regs, mm_reg1_t *regs, const mm128_t *a, int32_t n, const uint64_t *mini_pos);
 
@@ -109,6 +112,16 @@ void mm_err_puts(const char *str);
 void mm_err_fwrite(const void *p, size_t size, size_t nitems, FILE *fp);
 void mm_err_fread(void *p, size_t size, size_t nitems, FILE *fp);
 
+static inline float mg_log2(float x) // NB: this doesn't work when x<2
+{
+	union { float f; uint32_t i; } z = { x };
+	float log_2 = ((z.i >> 23) & 255) - 128;
+	z.i &= ~(255 << 23);
+	z.i += 127 << 23;
+	log_2 += (-0.34484843f * z.f + 2.02466578f) * z.f - 0.67487759f;
+	return log_2;
+}
+
 #ifdef __cplusplus
 }
 #endif

options.c

@@ -1,7 +1,7 @@
 #include <stdio.h>
 #include <limits.h>
 #include "mmpriv.h"
-
+extern bool enable_vect_dp_chaining;
 void mm_idxopt_init(mm_idxopt_t *opt)
 {
 	memset(opt, 0, sizeof(mm_idxopt_t));
@@ -53,6 +53,9 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->mini_batch_size = 500000000;
 	opt->max_sw_mat = 100000000;
 
+	opt->rank_min_len = 500;
+	opt->rank_frac = 0.9f;
+
 	opt->pe_ori = 0; // FF
 	opt->pe_bonus = 33;
 }
@@ -91,6 +94,9 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		mo->bw = mo->bw_long = 2000;
 		mo->occ_dist = 0;
 	} else if (strcmp(preset, "map10k") == 0 || strcmp(preset, "map-pb") == 0) {
+		#if defined (PARALLEL_CHAINING) && (defined(__AVX2__)) && (!defined(__AVX512BW__))
+		enable_vect_dp_chaining = false;
+		#endif
 		io->flag |= MM_I_HPC, io->k = 19;
 	} else if (strcmp(preset, "ava-pb") == 0) {
 		io->flag |= MM_I_HPC, io->k = 19, io->w = 5;
@@ -99,6 +105,9 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		mo->bw_long = mo->bw;
 		mo->occ_dist = 0;
 	} else if (strcmp(preset, "map-hifi") == 0 || strcmp(preset, "map-ccs") == 0) {
+		#if defined (PARALLEL_CHAINING) && (defined(__AVX2__)) && (!defined(__AVX512BW__))
+		enable_vect_dp_chaining = false;
+		#endif
 		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;
@@ -218,5 +227,10 @@ int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
 			fprintf(stderr, "[ERROR]\033[1;31m -X/-P and --secondary=no can't be applied at the same time\033[0m\n");
 		return -5;
 	}
+	if ((mo->flag & MM_F_QSTRAND) && ((mo->flag & (MM_F_OUT_SAM|MM_F_SPLICE|MM_F_FRAG_MODE)) || (io->flag & MM_I_HPC))) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m --qstrand doesn't work with -a, -H, --frag or --splice\033[0m\n");
+		return -5;
+	}
 	return 0;
 }

python/cmappy.pxd

@@ -45,6 +45,9 @@ cdef extern from "minimap.h":
 		int anchor_ext_len, anchor_ext_shift
 		float max_clip_ratio
 
+		int rank_min_len
+		float rank_frac
+
 		int pe_ori, pe_bonus
 
 		float mid_occ_frac
@@ -53,6 +56,7 @@ cdef extern from "minimap.h":
 		int32_t max_occ
 		int64_t mini_batch_size
 		int64_t max_sw_mat
+		int64_t cap_kalloc
 
 		const char *split_prefix
 

python/mappy.pyx

@@ -3,7 +3,7 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
-__version__ = '2.20'
+__version__ = '2.22'
 
 cmappy.mm_reset_timer()
 
@@ -82,7 +82,7 @@ cdef class Alignment:
 
 	@property
 	def cigar_str(self):
-		return "".join(map(lambda x: str(x[0]) + 'MIDNSH'[x[1]], self._cigar))
+		return "".join(map(lambda x: str(x[0]) + 'MIDNSHP=XB'[x[1]], self._cigar))
 
 	def __str__(self):
 		if self._strand > 0: strand = '+'

seed.c

@@ -1,9 +1,40 @@
 #include "mmpriv.h"
 #include "kalloc.h"
 #include "ksort.h"
+#include <stdlib.h>
+#include<algorithm>
+#include <x86intrin.h>
+
+#ifdef LISA_HASH
+#include "lisa_hash.h"
+extern lisa_hash<uint64_t, uint64_t> *lh;
+#endif
+extern uint64_t minimizer_lookup_time;
 
 mm_seed_t *mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv, int32_t *n_m_)
 {
+//#ifdef MANUAL_PROFILING
+//	uint64_t lookup_start = __rdtsc();
+//#endif
+
+#ifdef LISA_HASH
+//-----------------------------------
+	uint64_t** cr_batch = (uint64_t**) malloc((mv->n)*sizeof(uint64_t*));
+	int* t_batch = (int*)malloc((mv->n)*sizeof(int));
+	uint64_t* minimizers = (uint64_t*) malloc((mv->n)*sizeof(uint64_t));
+	int64_t* lisa_pos = (int64_t*) malloc((max(32, (int)mv->n))* sizeof(int64_t));
+
+	for (size_t i = 0; i < mv->n; i++) {
+		mm128_t *p = &mv->a[i];
+		minimizers[i] = p->x>>8;
+	}
+	
+	lh->mm_idx_get_batched(minimizers, mv->n, lisa_pos, cr_batch, t_batch); 
+//-----------------------------------
+
+#endif
+
+
 	mm_seed_t *m;
 	size_t i;
 	int32_t k;
@@ -14,7 +45,12 @@ mm_seed_t *mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv,
 		mm128_t *p = &mv->a[i];
 		uint32_t q_pos = (uint32_t)p->y, q_span = p->x & 0xff;
 		int t;
+#ifdef LISA_HASH
+		t = t_batch[i];
+		cr = cr_batch[i];
+#else		
 		cr = mm_idx_get(mi, p->x>>8, &t);
+#endif
 		if (t == 0) continue;
 		q = &m[k++];
 		q->q_pos = q_pos, q->q_span = q_span, q->cr = cr, q->n = t, q->seg_id = p->y >> 32;
@@ -22,7 +58,16 @@ mm_seed_t *mm_seed_collect_all(void *km, const mm_idx_t *mi, const mm128_v *mv,
 		if (i > 0 && p->x>>8 == mv->a[i - 1].x>>8) q->is_tandem = 1;
 		if (i < mv->n - 1 && p->x>>8 == mv->a[i + 1].x>>8) q->is_tandem = 1;
 	}
+#ifdef LISA_HASH
+	free(cr_batch);
+	free(t_batch);
+	free(minimizers);
+	free(lisa_pos);
+#endif
 	*n_m_ = k;
+//#ifdef MANUAL_PROFILING
+//	minimizer_lookup_time += __rdtsc() - lookup_start;
+//#endif
 	return m;
 }
 
@@ -46,19 +91,20 @@ void mm_seed_select(int32_t n, mm_seed_t *a, int len, int max_occ, int max_max_o
 				int32_t pe = i == n? len : (uint32_t)a[i].q_pos>>1;
 				int32_t j, k, st = last0 + 1, en = i;
 				int32_t max_high_occ = (int32_t)((double)(pe - ps) / dist + .499);
-				//fprintf(stderr, "Y\t%d\t%d\n", ps, pe);
-				if (max_high_occ > MAX_MAX_HIGH_OCC)
-					max_high_occ = MAX_MAX_HIGH_OCC;
-				for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)
-					b[k] = (uint64_t)a[j].n<<32 | j;
-				ks_heapmake_uint64_t(k, b); // initialize the binomial heap
-				for (; j < en; ++j) { // if there are more, choose top max_high_occ
-					if (a[j].n < (int32_t)(b[0]>>32)) { // then update the heap
-						b[0] = (uint64_t)a[j].n<<32 | j;
-						ks_heapdown_uint64_t(0, k, b);
+				if (max_high_occ > 0) {
+					if (max_high_occ > MAX_MAX_HIGH_OCC)
+						max_high_occ = MAX_MAX_HIGH_OCC;
+					for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)
+						b[k] = (uint64_t)a[j].n<<32 | j;
+					ks_heapmake_uint64_t(k, b); // initialize the binomial heap
+					for (; j < en; ++j) { // if there are more, choose top max_high_occ
+						if (a[j].n < (int32_t)(b[0]>>32)) { // then update the heap
+							b[0] = (uint64_t)a[j].n<<32 | j;
+							ks_heapdown_uint64_t(0, k, b);
+						}
 					}
+					for (j = 0; j < k; ++j) a[(uint32_t)b[j]].flt = 1;
 				}
-				for (j = 0; j < k; ++j) a[(uint32_t)b[j]].flt = 1;
 				for (j = st; j < en; ++j) a[j].flt ^= 1;
 				for (j = st; j < en; ++j)
 					if (a[j].n > max_max_occ)

setup.py

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

tex/minimap2.bib

@@ -338,3 +338,114 @@
 	Title = {Introducing difference recurrence relations for faster semi-global alignment of long sequences},
 	Volume = {19},
 	Year = {2018}}
+
+@article{Li:2018ab,
+	Author = {Li, Heng},
+	Journal = {Bioinformatics},
+	Pages = {3094-3100},
+	Title = {Minimap2: pairwise alignment for nucleotide sequences},
+	Volume = {34},
+	Year = {2018}}
+
+@article{Jain:2020aa,
+	Author = {Jain, Chirag and others},
+	Journal = {Bioinformatics},
+	Pages = {i111-i118},
+	Title = {Weighted minimizer sampling improves long read mapping},
+	Volume = {36},
+	Year = {2020}}
+
+@article{Miga:2020aa,
+	Author = {Miga, Karen H and others},
+	Journal = {Nature},
+	Pages = {79-84},
+	Title = {Telomere-to-telomere assembly of a complete human {X} chromosome},
+	Volume = {585},
+	Year = {2020}}
+
+@article {Jain2020.11.01.363887,
+	author = {Jain, Chirag and others},
+	title = {A long read mapping method for highly repetitive reference sequences},
+	elocation-id = {2020.11.01.363887},
+	year = {2020},
+	doi = {10.1101/2020.11.01.363887},
+	publisher = {Cold Spring Harbor Laboratory},
+	abstract = {About 5-10\% of the human genome remains inaccessible for functional analysis due to the presence of repetitive sequences such as segmental duplications and tandem repeat arrays. To enable high-quality resequencing of personal genomes, it is crucial to support end-to-end genome variant discovery using repeat-aware read mapping methods. In this study, we highlight the fact that existing long read mappers often yield incorrect alignments and variant calls within long, near-identical repeats, as they remain vulnerable to allelic bias. In the presence of a non-reference allele within a repeat, a read sampled from that region could be mapped to an incorrect repeat copy because the standard pairwise sequence alignment scoring system penalizes true variants.To address the above problem, we propose a novel, long read mapping method that addresses allelic bias by making use of minimal confidently alignable substrings (MCASs). MCASs are formulated as minimal length substrings of a read that have unique alignments to a reference locus with sufficient mapping confidence (i.e., a mapping quality score above a user-specified threshold). This approach treats each read mapping as a collection of confident sub-alignments, which is more tolerant of structural variation and more sensitive to paralog-specific variants (PSVs) within repeats. We mathematically define MCASs and discuss an exact algorithm as well as a practical heuristic to compute them. The proposed method, referred to as Winnowmap2, is evaluated using simulated as well as real long read benchmarks using the recently completed gapless assemblies of human chromosomes X and 8 as a reference. We show that Winnowmap2 successfully addresses the issue of allelic bias, enabling more accurate downstream variant calls in repetitive sequences. As an example, using simulated PacBio HiFi reads and structural variants in chromosome 8, Winnowmap2 alignments achieved the lowest false-negative and false-positive rates (1.89\%, 1.89\%) for calling structural variants within near-identical repeats compared to minimap2 (39.62\%, 5.88\%) and NGMLR (56.60\%, 36.11\%) respectively.Winnowmap2 code is accessible at https://github.com/marbl/WinnowmapCompeting Interest StatementThe authors have declared no competing interest.},
+	URL = {https://www.biorxiv.org/content/early/2020/11/02/2020.11.01.363887},
+	eprint = {https://www.biorxiv.org/content/early/2020/11/02/2020.11.01.363887.full.pdf},
+	journal = {bioRxiv}
+}
+
+@article{Li:2020aa,
+	Author = {Li, Heng and others},
+	Journal = {Genome Biol},
+	Pages = {265},
+	Title = {The design and construction of reference pangenome graphs with minigraph},
+	Volume = {21},
+	Year = {2020}}
+
+@article{Ren:2021aa,
+	Author = {Ren, Jingwen and Chaisson, Mark J P},
+	Journal = {PLoS Comput Biol},
+	Pages = {e1009078},
+	Title = {lra: A long read aligner for sequences and contigs},
+	Volume = {17},
+	Year = {2021}}
+
+@inproceedings{DBLP:conf/wabi/AbouelhodaO03,
+	Author = {Mohamed Ibrahim Abouelhoda and Enno Ohlebusch},
+	Booktitle = {Algorithms in Bioinformatics, Third International Workshop, {WABI} 2003, Budapest, Hungary, September 15-20, 2003, Proceedings},
+	Crossref = {DBLP:conf/wabi/2003},
+	Pages = {1--16},
+	Title = {A Local Chaining Algorithm and Its Applications in Comparative Genomics},
+	Year = {2003}}
+
+@article{Ono:2021aa,
+	Author = {Ono, Yukiteru and others},
+	Journal = {Bioinformatics},
+	Pages = {589-595},
+	Title = {{PBSIM2}: a simulator for long-read sequencers with a novel generative model of quality scores},
+	Volume = {37},
+	Year = {2021}}
+
+@article{Sedlazeck:2018ab,
+	Author = {Sedlazeck, Fritz J and others},
+	Journal = {Nat Methods},
+	Pages = {461-468},
+	Title = {Accurate detection of complex structural variations using single-molecule sequencing},
+	Volume = {15},
+	Year = {2018}}
+
+@article{Jeffares:2017aa,
+	Author = {Jeffares, Daniel C and others},
+	Journal = {Nat Commun},
+	Pages = {14061},
+	Title = {Transient structural variations have strong effects on quantitative traits and reproductive isolation in fission yeast},
+	Volume = {8},
+	Year = {2017}}
+
+@article{Zook:2020aa,
+	Author = {Zook, Justin M and others},
+	Journal = {Nat Biotechnol},
+	Pages = {1347-1355},
+	Title = {A robust benchmark for detection of germline large deletions and insertions},
+	Volume = {38},
+	Year = {2020}}
+
+@article{Harpak:2017aa,
+	Author = {Harpak, Arbel and others},
+	Journal = {Proc Natl Acad Sci U S A},
+	Pages = {12779-12784},
+	Title = {Frequent nonallelic gene conversion on the human lineage and its effect on the divergence of gene duplicates},
+	Volume = {114},
+	Year = {2017}}
+
+@article{Li:2018aa,
+	Author = {Li, Heng and others},
+	Journal = {Nat Methods},
+	Month = {Aug},
+	Number = {8},
+	Pages = {595-597},
+	Title = {A synthetic-diploid benchmark for accurate variant-calling evaluation},
+	Volume = {15},
+	Year = {2018}}

tex/mm2-update.tex

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