Update README.md

.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

.travis.yml

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

Makefile

@@ -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,13 +45,14 @@ PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite
 LIBS=		-lm -lz -lpthread
 
-ifneq ($(aarch64),)
-	arm_neon=1
+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
@@ -30,12 +68,12 @@ endif
 
 ifneq ($(asan),)
 	CFLAGS+=-fsanitize=address
-	LIBS+=-fsanitize=address -ldl
+	LIBS+=-fsanitize=address
 endif
 
 ifneq ($(tsan),)
 	CFLAGS+=-fsanitize=thread
-	LIBS+=-fsanitize=thread -ldl
+	LIBS+=-fsanitize=thread
 endif
 
 .PHONY:all extra clean depend
@@ -60,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)

NEWS.md

@@ -1,149 +1,3 @@
-Release 2.28-r1209 (27 March 2024)
-----------------------------------
-
-Notable changes to minimap2:
-
- * Bugfix: `--MD` was not working properly due to the addition of `--ds` in the
-   last release (#1181 and #1182).
-
- * New feature: added an experimental preset `lq:hqae` for aligning accurate
-   long reads back to their assembly. It has been observed that `map-hifi` and
-   `lr:hq` may produce many wrong alignments around centromeres when accurate
-   long reads (PacBio HiFi or Nanopore duplex/Q20+) are mapped to a diploid
-   assembly constructed from them. This new preset produces much more accurate
-   alignment. It is still experimental and may be subjective to changes in
-   future.
-
- * Change: reduced the default `--cap-kalloc` to 500m to lower the peak
-   memory consumption (#855).
-
-Notable changes to mappy:
-
- * Bugfix: mappy option struct was out of sync with minimap2 (#1177).
-
-Minimap2 should output identical alignments to v2.27.
-
-(2.28: 27 March 2024, r1209)
-
-
-
-Release 2.27-r1193 (12 March 2024)
-----------------------------------
-
-Notable changes to minimap2:
-
- * New feature: added the `lr:hq` preset for accurate long reads at ~1% error
-   rate. This was suggested by Oxford Nanopore developers (#1127). It is not
-   clear if this preset also works well for PacBio HiFi reads.
-
- * New feature: added the `map-iclr` preset for Illumina Complete Long Reads
-   (#1069), provided by Illumina developers.
-
- * New feature: added option `-b` to specify mismatch penalty for base
-   transitions (i.e. A-to-G or C-to-T changes).
-
- * New feature: added option `--ds` to generate a new `ds:Z` tag that
-   indicates uncertainty in INDEL positions. It is an extension to `cs`. The
-   `mgutils-es6.js` script in minigraph parses `ds`.
-
- * Bugfix: avoided a NULL pointer dereference (#1154). This would not have an
-   effect on most systems but would still be good to fix.
-
- * Bugfix: reverted the value of `ms:i` to pre-2.22 versions (#1146). This was
-   an oversight. See fcd4df2 for details.
-
-Notable changes to paftools.js and mappy:
-
- * New feature: expose `bw_long` to mappy's Aligner class (#1124).
-
- * Bugfix: fixed several compatibility issues with k8 v1.0 (#1161 and #1166).
-   Subcommands "call", "pbsim2fq" and "mason2fq" were not working with v1.0.
-
-Minimap2 should output identical alignments to v2.26, except the ms tag.
-
-(2.27: 12 March 2024, r1193)
-
-
-
-Release 2.26-r1175 (29 April 2023)
-----------------------------------
-
-Fixed the broken Python package. This is the only change.
-
-(2.26: 25 April 2023, r1173)
-
-
-
-Release 2.25-r1173 (25 April 2023)
-----------------------------------
-
-Notable changes:
-
- * Improvement: use the miniprot splice model for RNA-seq alignment by default.
-   This model considers non-GT-AG splice sites and leads to slightly higher
-   (<0.1%) accuracy and sensitivity on real human data.
-
- * Change: increased the default `-I` to `8G` such that minimap2 would create a
-   uni-part index for a pair of mammalian genomes. This change may increase the
-   memory for all-vs-all read overlap alignment given large datasets.
-
- * New feature: output the sequences in secondary alignments with option
-   `--secondary-seq` (#687).
-
- * Bugfix: --rmq was not parsed correctly (#1010)
-
- * Bugfix: possibly incorrect coordinate when applying end bonus to the target
-   sequence (#1025). This is a ksw2 bug. It does not affect minimap2 as
-   minimap2 is not using the affected feature.
-
- * Improvement: incorporated several changes for better compatibility with
-   Windows (#1051) and for minimap2 integration at Oxford Nanopore Technologies
-   (#1048 and #1033).
-
- * Improvement: output the HD-line in SAM output (#1019).
-
- * Improvement: check minimap2 index file in mappy to prevent segmentation
-   fault for certain indices (#1008).
-
-For genomic sequences, minimap2 should give identical output to v2.24.
-Long-read RNA-seq alignment may occasionally differ from previous versions.
-
-(2.25: 25 April 2023, r1173)
-
-
-
-Release 2.24-r1122 (26 December 2021)
--------------------------------------
-
-This release improves alignment around long poorly aligned regions. Older
-minimap2 may chain through such regions in rare cases which may result in
-missing alignments later. The issue has become worse since the the change of
-the chaining algorithm in v2.19. v2.23 implements an incomplete remedy. This
-release provides a better solution with a X-drop-like heuristic and by enabling
-two-bandwidth chaining in the assembly mode.
-
-(2.24: 26 December 2021, r1122)
-
-
-
-Release 2.23-r1111 (18 November 2021)
--------------------------------------
-
-Notable changes:
-
- * Bugfix: fixed missing alignments around long inversions (#806 and #816).
-   This bug affected v2.19 through v2.22.
-
- * Improvement: avoid extremely long mapping time for pathologic reads with
-   highly repeated k-mers not in the reference (#771). Use --q-occ-frac=0
-   to disable the new heuristic.
-
- * Change: use --cap-kalloc=1g by default.
-
-(2.23: 18 November 2021, r1111)
-
-
-
 Release 2.22-r1101 (7 August 2021)
 ----------------------------------
 

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

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

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.28/minimap2-2.28_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.28_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 -

format.c

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

hit.c

@@ -252,7 +252,7 @@ void mm_sync_regs(void *km, int n_regs, mm_reg1_t *regs) // keep mm_reg1_t::{id,
 	mm_set_sam_pri(n_regs, regs);
 }
 
-void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int check_strand, int min_strand_sc, int *n_, mm_reg1_t *r)
+void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_, mm_reg1_t *r)
 {
 	if (pri_ratio > 0.0f && *n_ > 0) {
 		int i, k, n = *n_, n_2nd = 0;
@@ -264,9 +264,6 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int chec
 				if (!(r[i].qs == r[p].qs && r[i].qe == r[p].qe && r[i].rid == r[p].rid && r[i].rs == r[p].rs && r[i].re == r[p].re)) // not identical hits
 					r[k++] = r[i], ++n_2nd;
 				else if (r[i].p) free(r[i].p);
-			} else if (check_strand && n_2nd < best_n && r[i].score > min_strand_sc && r[i].rev != r[p].rev) {
-				r[i].strand_retained = 1;
-				r[k++] = r[i], ++n_2nd;
 			} else if (r[i].p) free(r[i].p);
 		}
 		if (k != n) mm_sync_regs(km, k, r); // removing hits requires sync()
@@ -274,19 +271,6 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int chec
 	}
 }
 
-int mm_filter_strand_retained(int n_regs, mm_reg1_t *r)
-{
-	int i, k;
-	for (i = k = 0; i < n_regs; ++i) {
-		int p = r[i].parent;
-		if (!r[i].strand_retained || r[i].div < r[p].div * 5.0f || r[i].div < 0.01f) {
-			if (k < i) r[k++] = r[i];
-			else ++k;
-		}
-	}
-	return k;
-}
-
 void mm_filter_regs(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs)
 { // NB: after this call, mm_reg1_t::parent can be -1 if its parent filtered out
 	int i, k;

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

kalloc.c

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

kalloc.h

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

ksw2.h

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

ksw2_extd2_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);

ksw2_extd2_sse.c

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

ksw2_exts2_sse.c

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

ksw2_extz2_sse.c

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

lchain.c

@@ -5,26 +5,17 @@
 #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 int64_t mg_chain_bk_end(int32_t max_drop, const mm128_t *z, const int32_t *f, const int64_t *p, int32_t *t, int64_t k)
-{
-	int64_t i = z[k].y, end_i = -1, max_i = i;
-	int32_t max_s = 0;
-	if (i < 0 || t[i] != 0) return i;
-	do {
-		int32_t s;
-		t[i] = 2;
-		end_i = i = p[i];
-		s = i < 0? z[k].x : (int32_t)z[k].x - f[i];
-		if (s > max_s) max_s = s, max_i = i;
-		else if (max_s - s > max_drop) break;
-	} while (i >= 0 && t[i] == 0);
-	for (i = z[k].y; i >= 0 && i != end_i; i = p[i]) // reset modified t[]
-		t[i] = 0;
-	return max_i;
-}
+#ifdef MANUAL_PROFILING
+extern uint64_t dp_time, rmq_time, rmq_t1, rmq_t2, rmq_t3, rmq_t4;
+#endif
 
-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 max_drop, int32_t *n_u_, int32_t *n_v_)
+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;
 	uint64_t *u;
@@ -35,39 +26,33 @@ uint64_t *mg_chain_backtrack(void *km, int64_t n, const int32_t *f, const int64_
 	for (i = 0, n_z = 0; i < n; ++i) // precompute n_z
 		if (f[i] >= min_sc) ++n_z;
 	if (n_z == 0) return 0;
-	z = Kmalloc(km, mm128_t, n_z);
+	KMALLOC(km, z, n_z);
 	for (i = 0, k = 0; i < n; ++i) // populate z[]
 		if (f[i] >= min_sc) z[k].x = f[i], z[k++].y = i;
 	radix_sort_128x(z, z + n_z);
 
 	memset(t, 0, n * 4);
 	for (k = n_z - 1, n_v = n_u = 0; k >= 0; --k) { // precompute n_u
-		if (t[z[k].y] == 0) {
-			int64_t n_v0 = n_v, end_i;
-			int32_t sc;
-			end_i = mg_chain_bk_end(max_drop, z, f, p, t, k);
-			for (i = z[k].y; i != end_i; i = p[i])
-				++n_v, t[i] = 1;
-			sc = i < 0? z[k].x : (int32_t)z[k].x - f[i];
-			if (sc >= min_sc && n_v > n_v0 && n_v - n_v0 >= min_cnt)
-				++n_u;
-			else n_v = n_v0;
-		}
+		int64_t n_v0 = n_v;
+		int32_t sc;
+		for (i = z[k].y; i >= 0 && t[i] == 0; i = p[i])
+			++n_v, t[i] = 1;
+		sc = i < 0? z[k].x : (int32_t)z[k].x - f[i];
+		if (sc >= min_sc && n_v > n_v0 && n_v - n_v0 >= min_cnt)
+			++n_u;
+		else n_v = n_v0;
 	}
-	u = Kmalloc(km, uint64_t, n_u);
+	KMALLOC(km, u, n_u);
 	memset(t, 0, n * 4);
 	for (k = n_z - 1, n_v = n_u = 0; k >= 0; --k) { // populate u[]
-		if (t[z[k].y] == 0) {
-			int64_t n_v0 = n_v, end_i;
-			int32_t sc;
-			end_i = mg_chain_bk_end(max_drop, z, f, p, t, k);
-			for (i = z[k].y; i != end_i; i = p[i])
-				v[n_v++] = i, t[i] = 1;
-			sc = i < 0? z[k].x : (int32_t)z[k].x - f[i];
-			if (sc >= min_sc && n_v > n_v0 && n_v - n_v0 >= min_cnt)
-				u[n_u++] = (uint64_t)sc << 32 | (n_v - n_v0);
-			else n_v = n_v0;
-		}
+		int64_t n_v0 = n_v;
+		int32_t sc;
+		for (i = z[k].y; i >= 0 && t[i] == 0; i = p[i])
+			v[n_v++] = i, t[i] = 1;
+		sc = i < 0? z[k].x : (int32_t)z[k].x - f[i];
+		if (sc >= min_sc && n_v > n_v0 && n_v - n_v0 >= min_cnt)
+			u[n_u++] = (uint64_t)sc << 32 | (n_v - n_v0);
+		else n_v = n_v0;
 	}
 	kfree(km, z);
 	assert(n_v < INT32_MAX);
@@ -82,7 +67,7 @@ static mm128_t *compact_a(void *km, int32_t n_u, uint64_t *u, int32_t n_v, int32
 	int64_t i, j, k;
 
 	// write the result to b[]
-	b = Kmalloc(km, mm128_t, n_v);
+	KMALLOC(km, b, n_v);
 	for (i = 0, k = 0; i < n_u; ++i) {
 		int32_t k0 = k, ni = (int32_t)u[i];
 		for (j = 0; j < ni; ++j)
@@ -91,13 +76,13 @@ static mm128_t *compact_a(void *km, int32_t n_u, uint64_t *u, int32_t n_v, int32
 	kfree(km, v);
 
 	// sort u[] and a[] by the target position, such that adjacent chains may be joined
-	w = Kmalloc(km, mm128_t, n_u);
+	KMALLOC(km, w, n_u);
 	for (i = k = 0; i < n_u; ++i) {
 		w[i].x = b[k].x, w[i].y = (uint64_t)k<<32|i;
 		k += (int32_t)u[i];
 	}
 	radix_sort_128x(w, w + n_u);
-	u2 = Kmalloc(km, uint64_t, n_u);
+	KMALLOC(km, u2, n_u);
 	for (i = k = 0; i < n_u; ++i) {
 		int32_t j = (int32_t)w[i].y, n = (int32_t)u[j];
 		u2[i] = u[j];
@@ -112,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;
@@ -134,11 +169,18 @@ static inline int32_t comput_sc(const mm128_t *ai, const mm128_t *aj, int32_t ma
 			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;
 }
 
 /* Input:
- *   a[].x: rev<<63 | tid<<32 | tpos
+ *   a[].x: tid<<33 | rev<<32 | tpos
  *   a[].y: flags<<40 | q_span<<32 | q_pos
  * Output:
  *   n_u: #chains
@@ -148,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, max_drop = bw;
-	int64_t *p, i, j, max_ii, st = 0;
-	uint64_t *u;
+	///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);
@@ -159,11 +209,49 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
 	}
 	if (max_dist_x < bw) max_dist_x = bw;
 	if (max_dist_y < bw && !is_cdna) max_dist_y = bw;
-	if (is_cdna) max_drop = INT32_MAX;
-	p = Kmalloc(km, int64_t, n);
-	f = Kmalloc(km, int32_t, n);
-	v = Kmalloc(km, int32_t, n);
-	t = Kcalloc(km, int32_t, n);
+	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) {
@@ -171,9 +259,11 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
 		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);
+			++n_iter;
 			if (sc == INT32_MIN) continue;
 			sc += f[j];
 			if (sc > max_f) {
@@ -186,33 +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;
-		//fprintf(stderr, "X1\t%ld\t%ld:%d\t%ld\t%ld:%d\t%ld\t%ld\n", (long)i, (long)(a[i].x>>32), (int32_t)a[i].x, (long)max_j, max_j<0?-1L:(long)(a[max_j].x>>32), max_j<0?-1:(int32_t)a[max_j].x, (long)max_f, (long)v[i]);
 	}
 
-	u = mg_chain_backtrack(km, n, f, p, v, t, min_cnt, min_sc, max_drop, &n_u, &n_v);
+}
+//#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);
 }
 
@@ -250,8 +369,13 @@ 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)
 {
-	int32_t *f,*t, *v, n_u, n_v, mmax_f = 0, max_rmq_size = 0, max_drop = bw;
-	int64_t *p, i, i0, st = 0, st_inner = 0;
+#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;
 	lc_elem_t *root = 0, *root_inner = 0;
 	void *mem_mp = 0;
@@ -263,12 +387,11 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 		return 0;
 	}
 	if (max_dist < bw) max_dist = bw;
-	if (max_dist_inner < 0) max_dist_inner = 0;
-	if (max_dist_inner > max_dist) max_dist_inner = max_dist;
-	p = Kmalloc(km, int64_t, n);
-	f = Kmalloc(km, int32_t, n);
-	t = Kcalloc(km, int32_t, n);
-	v = Kmalloc(km, int32_t, n);
+	if (max_dist_inner <= 0 || max_dist_inner >= max_dist) max_dist_inner = 0;
+	KMALLOC(km, p, n);
+	KMALLOC(km, f, n);
+	KCALLOC(km, t, n);
+	KMALLOC(km, v, n);
 	mem_mp = km_init2(km, 0x10000);
 	mp = kmp_init_rmq(mem_mp);
 
@@ -278,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) {
@@ -292,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) {
@@ -301,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;
@@ -311,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;
@@ -326,11 +467,15 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 				krmq_interval(lc_elem, root_inner, &s, &lo, &hi);
 				if (lo) {
 					const lc_elem_t *q;
-					int32_t width;
+					int32_t width, n_rmq_iter = 0;
 					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;
 						sc = f[j] + comput_sc_simple(&a[i], &a[j], chn_pen_gap, chn_pen_skip, 0, &width);
 						if (width <= bw) {
@@ -344,10 +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;
@@ -357,12 +507,15 @@ mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_ski
 	}
 	km_destroy(mem_mp);
 
-	u = mg_chain_backtrack(km, n, f, p, v, t, min_cnt, min_sc, max_drop, &n_u, &n_v);
+	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);
 	if (n_u == 0) {
 		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,6 +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"
+
+//#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,14 +143,6 @@ static ko_longopt_t long_options[] = {
 	{ "rmq",            ko_optional_argument, 347 },
 	{ "qstrand",        ko_no_argument,       348 },
 	{ "cap-kalloc",     ko_required_argument, 349 },
-	{ "q-occ-frac",     ko_required_argument, 350 },
-	{ "chain-skip-scale",ko_required_argument,351 },
-	{ "print-chains",   ko_no_argument,       352 },
-	{ "no-hash-name",   ko_no_argument,       353 },
-	{ "secondary-seq",  ko_no_argument,       354 },
-	{ "ds",             ko_no_argument,       355 },
-	{ "rmq-inner",      ko_required_argument, 356 },
-	{ "dbg-seed-occ",   ko_no_argument,       501 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
 	{ "version",        ko_no_argument,       'V' },
@@ -123,7 +186,13 @@ static inline void yes_or_no(mm_mapopt_t *opt, int64_t flag, int long_idx, const
 
 int main(int argc, char *argv[])
 {
-	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:b:O:E:m:N:Qu:R:hF:LC:yYPo:e:U:J:";
+	// 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;
 	mm_idxopt_t ipt;
@@ -137,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;
@@ -181,7 +252,6 @@ int main(int argc, char *argv[])
 		else if (c == 'm') opt.min_chain_score = atoi(o.arg);
 		else if (c == 'A') opt.a = atoi(o.arg);
 		else if (c == 'B') opt.b = atoi(o.arg);
-		else if (c == 'b') opt.transition = atoi(o.arg);
 		else if (c == 's') opt.min_dp_max = atoi(o.arg);
 		else if (c == 'C') opt.noncan = atoi(o.arg);
 		else if (c == 'I') ipt.batch_size = mm_parse_num(o.arg);
@@ -190,12 +260,7 @@ int main(int argc, char *argv[])
 		else if (c == 'R') rg = o.arg;
 		else if (c == 'h') fp_help = stdout;
 		else if (c == '2') opt.flag |= MM_F_2_IO_THREADS;
-		else if (c == 'J') {
-			int t;
-			t = atoi(o.arg);
-			if (t == 0) opt.flag |= MM_F_SPLICE_OLD;
-			else if (t == 1) opt.flag &= ~MM_F_SPLICE_OLD;
-		} else if (c == 'o') {
+		else if (c == 'o') {
 			if (strcmp(o.arg, "-") != 0) {
 				if (freopen(o.arg, "wb", stdout) == NULL) {
 					fprintf(stderr, "[ERROR]\033[1;31m failed to write the output to file '%s'\033[0m: %s\n", o.arg, strerror(errno));
@@ -236,19 +301,11 @@ int main(int argc, char *argv[])
 		else if (c == 341) opt.junc_bonus = atoi(o.arg); // --junc-bonus
 		else if (c == 342) opt.flag |= MM_F_SAM_HIT_ONLY; // --sam-hit-only
 		else if (c == 343) opt.chain_gap_scale = atof(o.arg); // --chain-gap-scale
-		else if (c == 351) opt.chain_skip_scale = atof(o.arg); // --chain-skip-scale
 		else if (c == 344) alt_list = o.arg; // --alt
 		else if (c == 345) opt.alt_drop = atof(o.arg); // --alt-drop
 		else if (c == 346) opt.mask_len = mm_parse_num(o.arg); // --mask-len
 		else if (c == 348) opt.flag |= MM_F_QSTRAND | MM_F_NO_INV; // --qstrand
 		else if (c == 349) opt.cap_kalloc = mm_parse_num(o.arg); // --cap-kalloc
-		else if (c == 350) opt.q_occ_frac = atof(o.arg); // --q-occ-frac
-		else if (c == 352) mm_dbg_flag |= MM_DBG_PRINT_CHAIN; // --print-chains
-		else if (c == 353) opt.flag |= MM_F_NO_HASH_NAME; // --no-hash-name
-		else if (c == 354) opt.flag |= MM_F_SECONDARY_SEQ; // --secondary-seq
-		else if (c == 355) opt.flag |= MM_F_OUT_DS; // --ds
-		else if (c == 356) opt.rmq_inner_dist = mm_parse_num(o.arg); // --rmq-inner
-		else if (c == 501) mm_dbg_flag |= MM_DBG_SEED_FREQ; // --dbg-seed-occ
 		else if (c == 330) {
 			fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n");
 		} else if (c == 314) { // --frag
@@ -273,8 +330,7 @@ int main(int argc, char *argv[])
 		} else if (c == 326) { // --dual
 			yes_or_no(&opt, MM_F_NO_DUAL, o.longidx, o.arg, 0);
 		} else if (c == 347) { // --rmq
-			if (o.arg) yes_or_no(&opt, MM_F_RMQ, o.longidx, o.arg, 1);
-			else opt.flag |= MM_F_RMQ;
+			yes_or_no(&opt, MM_F_RMQ, o.longidx, o.arg, 1);
 		} else if (c == 'S') {
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR | MM_F_OUT_CS_LONG;
 			if (mm_verbose >= 2)
@@ -335,7 +391,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -H           use homopolymer-compressed k-mer (preferrable for PacBio)\n");
 		fprintf(fp_help, "    -k INT       k-mer size (no larger than 28) [%d]\n", ipt.k);
 		fprintf(fp_help, "    -w INT       minimizer window size [%d]\n", ipt.w);
-		fprintf(fp_help, "    -I NUM       split index for every ~NUM input bases [8G]\n");
+		fprintf(fp_help, "    -I NUM       split index for every ~NUM input bases [4G]\n");
 		fprintf(fp_help, "    -d FILE      dump index to FILE []\n");
 		fprintf(fp_help, "  Mapping:\n");
 		fprintf(fp_help, "    -f FLOAT     filter out top FLOAT fraction of repetitive minimizers [%g]\n", opt.mid_occ_frac);
@@ -357,7 +413,6 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -z INT[,INT] Z-drop score and inversion Z-drop score [%d,%d]\n", opt.zdrop, opt.zdrop_inv);
 		fprintf(fp_help, "    -s INT       minimal peak DP alignment score [%d]\n", opt.min_dp_max);
 		fprintf(fp_help, "    -u CHAR      how to find GT-AG. f:transcript strand, b:both strands, n:don't match GT-AG [n]\n");
-		fprintf(fp_help, "    -J INT       splice mode. 0: original minimap2 model; 1: miniprot model [1]\n");
 		fprintf(fp_help, "  Input/Output:\n");
 		fprintf(fp_help, "    -a           output in the SAM format (PAF by default)\n");
 		fprintf(fp_help, "    -o FILE      output alignments to FILE [stdout]\n");
@@ -365,7 +420,6 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -R STR       SAM read group line in a format like '@RG\\tID:foo\\tSM:bar' []\n");
 		fprintf(fp_help, "    -c           output CIGAR in PAF\n");
 		fprintf(fp_help, "    --cs[=STR]   output the cs tag; STR is 'short' (if absent) or 'long' [none]\n");
-		fprintf(fp_help, "    --ds         output the ds tag, which is an extension to cs\n");
 		fprintf(fp_help, "    --MD         output the MD tag\n");
 		fprintf(fp_help, "    --eqx        write =/X CIGAR operators\n");
 		fprintf(fp_help, "    -Y           use soft clipping for supplementary alignments\n");
@@ -375,12 +429,12 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    --version    show version number\n");
 		fprintf(fp_help, "  Preset:\n");
 		fprintf(fp_help, "    -x STR       preset (always applied before other options; see minimap2.1 for details) []\n");
-		fprintf(fp_help, "                 - lr:hq - accurate long reads (error rate <1%%) against a reference genome\n");
-		fprintf(fp_help, "                 - splice/splice:hq - spliced alignment for long reads/accurate long reads\n");
+		fprintf(fp_help, "                 - map-pb/map-ont - PacBio CLR/Nanopore vs reference mapping\n");
+		fprintf(fp_help, "                 - map-hifi - PacBio HiFi reads vs reference mapping\n");
+		fprintf(fp_help, "                 - ava-pb/ava-ont - PacBio/Nanopore read overlap\n");
 		fprintf(fp_help, "                 - asm5/asm10/asm20 - asm-to-ref mapping, for ~0.1/1/5%% sequence divergence\n");
-		fprintf(fp_help, "                 - sr - short reads against a reference\n");
-		fprintf(fp_help, "                 - map-pb/map-hifi/map-ont/map-iclr - CLR/HiFi/Nanopore/ICLR vs reference mapping\n");
-		fprintf(fp_help, "                 - ava-pb/ava-ont - PacBio CLR/Nanopore read overlap\n");
+		fprintf(fp_help, "                 - splice/splice:hq - long-read/Pacbio-CCS spliced alignment\n");
+		fprintf(fp_help, "                 - sr - genomic short-read mapping\n");
 		fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of these and other advanced command-line options.\n");
 		return fp_help == stdout? 0 : 1;
 	}
@@ -389,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));
@@ -431,6 +486,9 @@ 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) {
@@ -438,6 +496,16 @@ int main(int argc, char *argv[])
 			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);
@@ -446,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);
 
@@ -470,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,17 @@
 #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;
+	int rep_len, frag_gap;
+};
 
 mm_tbuf_t *mm_tbuf_init(void)
 {
@@ -168,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;
@@ -200,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;
 }
 
@@ -207,7 +224,7 @@ static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_i
 {
 	if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
 		mm_set_parent(km, opt->mask_level, opt->mask_len, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
-		if (n_segs <= 1) mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, 1, opt->max_gap * 0.8, n_regs, regs);
+		if (n_segs <= 1) mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 		else mm_select_sub_multi(km, opt->pri_ratio, 0.2f, 0.7f, max_chain_gap_ref, mi->k*2, opt->best_n, n_segs, qlens, n_regs, regs);
 	}
 }
@@ -218,7 +235,7 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k
 	regs = mm_align_skeleton(km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs()
 	if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
 		mm_set_parent(km, opt->mask_level, opt->mask_len, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
-		mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, 0, opt->max_gap * 0.8, n_regs, regs);
+		mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 		mm_set_sam_pri(*n_regs, regs);
 	}
 	return regs;
@@ -235,7 +252,6 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	mm128_v mv = {0,0,0};
 	mm_reg1_t *regs0;
 	km_stat_t kmst;
-	float chn_pen_gap, chn_pen_skip;
 
 	for (i = 0, qlen_sum = 0; i < n_segs; ++i)
 		qlen_sum += qlens[i], n_regs[i] = 0, regs[i] = 0;
@@ -243,12 +259,11 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	if (qlen_sum == 0 || n_segs <= 0 || n_segs > MM_MAX_SEG) return;
 	if (opt->max_qlen > 0 && qlen_sum > opt->max_qlen) return;
 
-	hash  = qname && !(opt->flag & MM_F_NO_HASH_NAME)? __ac_X31_hash_string(qname) : 0;
+	hash  = qname? __ac_X31_hash_string(qname) : 0;
 	hash ^= __ac_Wang_hash(qlen_sum) + __ac_Wang_hash(opt->seed);
 	hash  = __ac_Wang_hash(hash);
 
 	collect_minimizers(b->km, opt, mi, n_segs, qlens, seqs, &mv);
-	if (opt->q_occ_frac > 0.0f) mm_seed_mz_flt(b->km, &mv, opt->mid_occ, opt->q_occ_frac);
 	if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 	else a = collect_seed_hits(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 
@@ -270,25 +285,39 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		if (max_chain_gap_ref < opt->max_gap) max_chain_gap_ref = opt->max_gap;
 	} else max_chain_gap_ref = opt->max_gap;
 
-	chn_pen_gap  = opt->chain_gap_scale * 0.01 * mi->k;
-	chn_pen_skip = opt->chain_skip_scale * 0.01 * mi->k;
 	if (opt->flag & MM_F_RMQ) {
 		a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
-						  chn_pen_gap, chn_pen_skip, n_a, a, &n_regs0, &u, b->km);
+						  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,
-						 chn_pen_gap, chn_pen_skip, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
+						 opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 	}
 
 	if (opt->bw_long > opt->bw && (opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN)) == 0 && n_segs == 1 && n_regs0 > 1) { // re-chain/long-join for long sequences
 		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,
-							  chn_pen_gap, chn_pen_skip, n_a, a, &n_regs0, &u, b->km);
+							  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;
@@ -311,7 +340,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 			if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 			else a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 			a = mg_lchain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score,
-							 chn_pen_gap, chn_pen_skip, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
+							 opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 		}
 	}
 	b->frag_gap = max_chain_gap_ref;
@@ -323,17 +352,15 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		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
 	}
 
-	if (mm_dbg_flag & (MM_DBG_PRINT_SEED|MM_DBG_PRINT_CHAIN))
+	if (mm_dbg_flag & MM_DBG_PRINT_SEED)
 		for (j = 0; j < n_regs0; ++j)
 			for (i = regs0[j].as; i < regs0[j].as + regs0[j].cnt; ++i)
 				fprintf(stderr, "CN\t%d\t%s\t%d\t%c\t%d\t%d\t%d\n", j, mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff),
 						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 && !(opt->flag&MM_F_QSTRAND)) {
+	if (!is_sr && !(opt->flag&MM_F_QSTRAND))
 		mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos);
-		n_regs0 = mm_filter_strand_retained(n_regs0, regs0);
-	}
 
 	if (n_segs == 1) { // uni-segment
 		regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], &n_regs0, regs0, a);
@@ -506,7 +533,7 @@ static void merge_hits(step_t *s)
 			mm_hit_sort(km, &s->n_reg[k], s->reg[k], opt->alt_drop);
 			mm_set_parent(km, opt->mask_level, opt->mask_len, s->n_reg[k], s->reg[k], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
 			if (!(opt->flag & MM_F_ALL_CHAINS)) {
-				mm_select_sub(km, opt->pri_ratio, s->p->mi->k*2, opt->best_n, 0, opt->max_gap * 0.8, &s->n_reg[k], s->reg[k]);
+				mm_select_sub(km, opt->pri_ratio, s->p->mi->k*2, opt->best_n, &s->n_reg[k], s->reg[k]);
 				mm_set_sam_pri(s->n_reg[k], s->reg[k]);
 			}
 			mm_set_mapq(km, s->n_reg[k], s->reg[k], opt->min_chain_score, opt->a, rep_len, !!(opt->flag & MM_F_SR));
@@ -564,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
@@ -598,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

@@ -5,46 +5,40 @@
 #include <stdio.h>
 #include <sys/types.h>
 
-#define MM_VERSION "2.28-r1209"
-
-#define MM_F_NO_DIAG       (0x001LL) // no exact diagonal hit
-#define MM_F_NO_DUAL       (0x002LL) // skip pairs where query name is lexicographically larger than target name
-#define MM_F_CIGAR         (0x004LL)
-#define MM_F_OUT_SAM       (0x008LL)
-#define MM_F_NO_QUAL       (0x010LL)
-#define MM_F_OUT_CG        (0x020LL)
-#define MM_F_OUT_CS        (0x040LL)
-#define MM_F_SPLICE        (0x080LL) // splice mode
-#define MM_F_SPLICE_FOR    (0x100LL) // match GT-AG
-#define MM_F_SPLICE_REV    (0x200LL) // match CT-AC, the reverse complement of GT-AG
-#define MM_F_NO_LJOIN      (0x400LL)
-#define MM_F_OUT_CS_LONG   (0x800LL)
-#define MM_F_SR            (0x1000LL)
-#define MM_F_FRAG_MODE     (0x2000LL)
-#define MM_F_NO_PRINT_2ND  (0x4000LL)
-#define MM_F_2_IO_THREADS  (0x8000LL)
-#define MM_F_LONG_CIGAR    (0x10000LL)
-#define MM_F_INDEPEND_SEG  (0x20000LL)
-#define MM_F_SPLICE_FLANK  (0x40000LL)
-#define MM_F_SOFTCLIP      (0x80000LL)
-#define MM_F_FOR_ONLY      (0x100000LL)
-#define MM_F_REV_ONLY      (0x200000LL)
-#define MM_F_HEAP_SORT     (0x400000LL)
-#define MM_F_ALL_CHAINS    (0x800000LL)
-#define MM_F_OUT_MD        (0x1000000LL)
-#define MM_F_COPY_COMMENT  (0x2000000LL)
-#define MM_F_EQX           (0x4000000LL) // use =/X instead of M
-#define MM_F_PAF_NO_HIT    (0x8000000LL) // output unmapped reads to PAF
-#define MM_F_NO_END_FLT    (0x10000000LL)
-#define MM_F_HARD_MLEVEL   (0x20000000LL)
-#define MM_F_SAM_HIT_ONLY  (0x40000000LL)
+#define MM_F_NO_DIAG       0x001 // no exact diagonal hit
+#define MM_F_NO_DUAL       0x002 // skip pairs where query name is lexicographically larger than target name
+#define MM_F_CIGAR         0x004
+#define MM_F_OUT_SAM       0x008
+#define MM_F_NO_QUAL       0x010
+#define MM_F_OUT_CG        0x020
+#define MM_F_OUT_CS        0x040
+#define MM_F_SPLICE        0x080 // splice mode
+#define MM_F_SPLICE_FOR    0x100 // match GT-AG
+#define MM_F_SPLICE_REV    0x200 // match CT-AC, the reverse complement of GT-AG
+#define MM_F_NO_LJOIN      0x400
+#define MM_F_OUT_CS_LONG   0x800
+#define MM_F_SR            0x1000
+#define MM_F_FRAG_MODE     0x2000
+#define MM_F_NO_PRINT_2ND  0x4000
+#define MM_F_2_IO_THREADS  0x8000
+#define MM_F_LONG_CIGAR    0x10000
+#define MM_F_INDEPEND_SEG  0x20000
+#define MM_F_SPLICE_FLANK  0x40000
+#define MM_F_SOFTCLIP      0x80000
+#define MM_F_FOR_ONLY      0x100000
+#define MM_F_REV_ONLY      0x200000
+#define MM_F_HEAP_SORT     0x400000
+#define MM_F_ALL_CHAINS    0x800000
+#define MM_F_OUT_MD        0x1000000
+#define MM_F_COPY_COMMENT  0x2000000
+#define MM_F_EQX           0x4000000 // use =/X instead of M
+#define MM_F_PAF_NO_HIT    0x8000000 // output unmapped reads to PAF
+#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_QSTRAND       (0x100000000LL)
 #define MM_F_NO_INV        (0x200000000LL)
-#define MM_F_NO_HASH_NAME  (0x400000000LL)
-#define MM_F_SPLICE_OLD    (0x800000000LL)
-#define MM_F_SECONDARY_SEQ (0x1000000000LL)	//output SEQ field for seqondary alignments using hard clipping
-#define MM_F_OUT_DS        (0x2000000000LL)
 
 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2
@@ -98,7 +92,6 @@ typedef struct {
 typedef struct {
 	uint32_t capacity;                  // the capacity of cigar[]
 	int32_t dp_score, dp_max, dp_max2;  // DP score; score of the max-scoring segment; score of the best alternate mappings
-	int32_t dp_max0;                    // DP score before mm_update_dp_max() adjustment
 	uint32_t n_ambi:30, trans_strand:2; // number of ambiguous bases; transcript strand: 0 for unknown, 1 for +, 2 for -
 	uint32_t n_cigar;                   // number of cigar operations in cigar[]
 	uint32_t cigar[];
@@ -115,7 +108,7 @@ typedef struct {
 	int32_t mlen, blen;     // seeded exact match length; seeded alignment block length
 	int32_t n_sub;          // number of suboptimal mappings
 	int32_t score0;         // initial chaining score (before chain merging/spliting)
-	uint32_t mapq:8, split:2, rev:1, inv:1, sam_pri:1, proper_frag:1, pe_thru:1, seg_split:1, seg_id:8, split_inv:1, is_alt:1, strand_retained:1, dummy:5;
+	uint32_t mapq:8, split:2, rev:1, inv:1, sam_pri:1, proper_frag:1, pe_thru:1, seg_split:1, seg_id:8, split_inv:1, is_alt:1, dummy:6;
 	uint32_t hash;
 	float div;
 	mm_extra_t *p;
@@ -142,7 +135,6 @@ typedef struct {
 	int min_cnt;         // min number of minimizers on each chain
 	int min_chain_score; // min chaining score
 	float chain_gap_scale;
-	float chain_skip_scale;
 	int rmq_size_cap, rmq_inner_dist;
 	int rmq_rescue_size;
 	float rmq_rescue_ratio;
@@ -155,7 +147,6 @@ typedef struct {
 	float alt_drop;
 
 	int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties
-	int transition; // transition mismatch score (A:G, C:T)
 	int sc_ambi; // score when one or both bases are "N"
 	int noncan;      // cost of non-canonical splicing sites
 	int junc_bonus;
@@ -172,7 +163,6 @@ typedef struct {
 	int pe_ori, pe_bonus;
 
 	float mid_occ_frac;  // only used by mm_mapopt_update(); see below
-	float q_occ_frac;
 	int32_t min_mid_occ, max_mid_occ;
 	int32_t mid_occ;     // ignore seeds with occurrences above this threshold
 	int32_t max_occ, max_max_occ, occ_dist;
@@ -196,11 +186,6 @@ typedef struct {
 } mm_idx_reader_t;
 
 // memory buffer for thread-local storage during mapping
-struct mm_tbuf_s {
-	void *km;
-	int rep_len, frag_gap;
-};
-
 typedef struct mm_tbuf_s mm_tbuf_t;
 
 // global variables
@@ -300,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
  *
@@ -328,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 "12 March 2024" "minimap2-2.28 (r1209)" "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
@@ -77,21 +77,8 @@ SAM format.
 Minimizer k-mer length [15]
 .TP
 .BI -w \ INT
-Minimizer window size [10]. A minimizer is the smallest k-mer
+Minimizer window size [2/3 of k-mer length]. A minimizer is the smallest k-mer
 in a window of w consecutive k-mers.
-.TP
-.BI -j \ INT
-Syncmer submer size [10]. Option
-.B -j
-and
-.B -w
-will override each: if
-.B -w
-is applied after
-.BR -j ,
-.B -j
-will have no effect, and vice versa.
-
 .TP
 .B -H
 Use homopolymer-compressed (HPC) minimizers. An HPC sequence is constructed by
@@ -101,17 +88,16 @@ on the HPC sequence.
 .BI -I \ NUM
 Load at most
 .I NUM
-target bases into RAM for indexing [8G]. If there are more than
+target bases into RAM for indexing [4G]. If there are more than
 .I NUM
 bases in
 .IR target.fa ,
 minimap2 needs to read
 .I query.fa
-multiple times to map it against each batch of target sequences. This would create a multi-part index.
+multiple times to map it against each batch of target sequences.
 .I NUM
 may be ending with k/K/m/M/g/G. NB: mapping quality is incorrect given a
-multi-part index. See also option
-.BR --split-prefix .
+multi-part index.
 .TP
 .B --idx-no-seq
 Don't store target sequences in the index. It saves disk space and memory but
@@ -165,16 +151,10 @@ Lower and upper bounds of k-mer occurrences [10,1000000]. The final k-mer occurr
 .BR -f }}.
 This option prevents excessively small or large
 .B -f
-estimated from the input reference. Available since r1034 and deprecating
+estimated from the input reference. It deprecates
 .B --min-occ-floor
 in earlier versions of minimap2.
 .TP
-.BI --q-occ-frac \ FLOAT
-Discard a query minimizer if its occurrence is higher than
-.I FLOAT
-fraction of query minimizers and than the reference occurrence threshold
-[0.01]. Set 0 to disable. Available since r1105.
-.TP
 .BI -e \ INT
 Sample a high-frequency minimizer every
 .I INT
@@ -268,11 +248,6 @@ or more of the shorter chain [0.5]
 Use the minigraph chaining algorithm [no]. The minigraph algorithm is better
 for aligning contigs through long INDELs.
 .TP
-.BI --rmq-inner \ NUM
-Apply full dynamic programming for anchors within distance
-.I NUM
-[1000].
-.TP
 .B --hard-mask-level
 Honor option
 .B -M
@@ -337,9 +312,6 @@ faster for short reads, but slower for long reads. [no]
 .B --no-pairing
 Treat two reads in a pair as independent reads. The mate related fields in SAM
 are still properly populated.
-.TP
-.B --no-hash-name
-Produce the same alignment for identical sequences regardless of their sequence names.
 .SS Alignment options
 .TP 10
 .BI -A \ INT
@@ -348,10 +320,6 @@ Matching score [2]
 .BI -B \ INT
 Mismatching penalty [4]
 .TP
-.BI -b \ INT
-Mismatching penalty for transitions [same as
-.BR -B ].
-.TP
 .BI -O \ INT1[,INT2]
 Gap open penalty [4,24]. If
 .I INT2
@@ -365,19 +333,10 @@ costs
 .RI min{ O1 + k * E1 , O2 + k * E2 }.
 In the splice mode, the second gap penalties are not used.
 .TP
-.BI -J \ INT
-Splice model [1]. 0 for the original minimap2 splice model that always penalizes non-GT-AG splicing;
-1 for the miniprot model that considers non-GT-AG. Option
-.B -C
-has no effect with the default
-.BR -J1 .
-.BR -J0 .
-.TP
 .BI -C \ INT
 Cost for a non-canonical GT-AG splicing (effective with
-.B --splice
-.BR -J0 )
-[0].
+.BR --splice )
+[0]
 .TP
 .BI -z \ INT1[,INT2]
 Truncate an alignment if the running alignment score drops too quickly along
@@ -468,7 +427,7 @@ Set 0 to disable [100m].
 .BI --cap-kalloc \ NUM
 Free thread-local kalloc memory reservoir if after the alignment the size of the reservoir above
 .IR NUM .
-Set 0 to disable [500m].
+Set 0 to disable [0].
 .SS Input/output options
 .TP 10
 .B -a
@@ -524,9 +483,6 @@ Output =/X CIGAR operators for sequence match/mismatch.
 .B -Y
 In SAM output, use soft clipping for supplementary alignments.
 .TP
-.B --secondary-seq
-In SAM output, show query sequences for secondary alignments.
-.TP
 .BI --seed \ INT
 Integer seed for randomizing equally best hits. Minimap2 hashes
 .I INT
@@ -587,70 +543,42 @@ are:
 Align noisy long reads of ~10% error rate to a reference genome. This is the
 default mode.
 .TP
-.B lr:hq
-Align accurate long reads (error rate <1%) to a reference genome
-.RB ( -k19
-.B -w19 -U50,500
-.BR -g10k ).
-This was recommended by ONT developers for recent Nanopore reads
-produced with chemistry v14 that can reach ~99% in accuracy.
-It was shown to work better for accurate Nanopore reads
-than
-.BR map-hifi .
-.TP
 .B map-hifi
 Align PacBio high-fidelity (HiFi) reads to a reference genome
-.RB ( -xlr:hq
-.B -A1 -B4 -O6,26 -E2,1
+.RB ( -k19
+.B -w19 -U50,500 -g10k -A1 -B4 -O6,26 -E2,1
 .BR -s200 ).
-It differs from
-.B lr:hq
-only in scoring. It has not been tested whether
-.B lr:hq
-would work better for PacBio HiFi reads.
 .TP
 .B map-pb
 Align older PacBio continuous long (CLR) reads to a reference genome
 .RB ( -Hk19 ).
-Note that this data type is effectively deprecated by HiFi.
-Unless you work on very old data, you probably want to use
-.B map-hifi
-or
-.BR lr:hq .
-.TP
-.B map-iclr
-Align Illumina Complete Long Reads (ICLR) to a reference genome
-.RB ( -k19
-.B -B6 -b4
-.BR -O10,50 ).
-This was recommended by Illumina developers.
 .TP
 .B asm5
 Long assembly to reference mapping
 .RB ( -k19
-.B -w19 -U50,500 --rmq -r1k,100k -g10k -A1 -B19 -O39,81 -E3,1 -s200 -z200
+.B -w19 -U50,500 --rmq -r100k -g10k -A1 -B19 -O39,81 -E3,1 -s200 -z200
 .BR -N50 ).
 Typically, the alignment will not extend to regions with 5% or higher sequence
-divergence. Use this preset if the average divergence is not much higher than 0.1%.
+divergence. Only use this preset if the average divergence is far below 5%.
 .TP
 .B asm10
 Long assembly to reference mapping
 .RB ( -k19
-.B -w19 -U50,500 --rmq -r1k,100k -g10k -A1 -B9 -O16,41 -E2,1 -s200 -z200
+.B -w19 -U50,500 --rmq -r100k -g10k -A1 -B9 -O16,41 -E2,1 -s200 -z200
 .BR -N50 ).
-Use this if the average divergence is around 1%.
+Up to 10% sequence divergence.
 .TP
 .B asm20
 Long assembly to reference mapping
 .RB ( -k19
-.B -w10 -U50,500 --rmq -r1k,100k -g10k -A1 -B4 -O6,26 -E2,1 -s200 -z200
+.B -w10 -U50,500 --rmq -r100k -g10k -A1 -B4 -O6,26 -E2,1 -s200 -z200
 .BR -N50 ).
-Use this if the average divergence is around several percent.
+Up to 20% sequence divergence.
 .TP
 .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
@@ -661,15 +589,15 @@ costs are different during chaining; 4) the computation of the
 tag ignores introns to demote hits to pseudogenes.
 .TP
 .B splice:hq
-Spliced alignment for accurate long RNA-seq reads such as PacBio iso-seq
+Long-read splice alignment for PacBio CCS reads
 .RB ( -xsplice
 .B -C5 -O6,24
 .BR -B4 ).
 .TP
 .B sr
-Short-read alignment without splicing
+Short single-end reads without splicing
 .RB ( -k21
-.B -w11 --sr --frag=yes -A2 -B8 -O12,32 -E2,1 -b0 -r100 -p.5 -N20 -f1000,5000 -n2 -m25
+.B -w11 --sr --frag=yes -A2 -B8 -O12,32 -E2,1 -b0 -r100 -p.5 -N20 -f1000,5000 -n2 -m20
 .B -s40 -g100 -2K50m --heap-sort=yes
 .BR --secondary=no ).
 .TP

misc/paftools.js

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

mmpriv.h

@@ -13,8 +13,6 @@
 #define MM_DBG_PRINT_QNAME   0x2
 #define MM_DBG_PRINT_SEED    0x4
 #define MM_DBG_PRINT_ALN_SEQ 0x8
-#define MM_DBG_PRINT_CHAIN   0x10
-#define MM_DBG_SEED_FREQ     0x20
 
 #define MM_SEED_LONG_JOIN  (1ULL<<40)
 #define MM_SEED_IGNORE     (1ULL<<41)
@@ -63,7 +61,6 @@ uint32_t ks_ksmall_uint32_t(size_t n, uint32_t arr[], size_t kk);
 void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, int is_hpc, mm128_v *p);
 
 mm_seed_t *mm_collect_matches(void *km, int *_n_m, int qlen, int max_occ, int max_max_occ, int dist, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos);
-void mm_seed_mz_flt(void *km, mm128_v *mv, int32_t q_occ_max, float q_occ_frac);
 
 double mm_event_identity(const mm_reg1_t *r);
 int mm_write_sam_hdr(const mm_idx_t *mi, const char *rg, const char *ver, int argc, char *argv[]);
@@ -80,6 +77,8 @@ int mm_idx_getseq2(const mm_idx_t *mi, int is_rev, uint32_t rid, uint32_t st, ui
 mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
 mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a, int is_qstrand);
 
+mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float gap_scale,
+					 int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
 mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip,
 					  int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
 mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_skip, int cap_rmq_size, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip,
@@ -91,9 +90,8 @@ void mm_sync_regs(void *km, int n_regs, mm_reg1_t *regs);
 int mm_squeeze_a(void *km, int n_regs, mm_reg1_t *regs, mm128_t *a);
 int mm_set_sam_pri(int n, mm_reg1_t *r);
 void mm_set_parent(void *km, float mask_level, int mask_len, int n, mm_reg1_t *r, int sub_diff, int hard_mask_level, float alt_diff_frac);
-void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int check_strand, int min_strand_sc, int *n_, mm_reg1_t *r);
+void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_, mm_reg1_t *r);
 void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int max_gap_ref, int min_diff, int best_n, int n_segs, const int *qlens, int *n_, mm_reg1_t *r);
-int mm_filter_strand_retained(int n_regs, mm_reg1_t *r);
 void mm_filter_regs(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs);
 void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r, float alt_diff_frac);
 void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr);

options.c

@@ -1,14 +1,14 @@
 #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));
 	opt->k = 15, opt->w = 10, opt->flag = 0;
 	opt->bucket_bits = 14;
 	opt->mini_batch_size = 50000000;
-	opt->batch_size = 8000000000ULL;
+	opt->batch_size = 4000000000ULL;
 }
 
 void mm_mapopt_init(mm_mapopt_t *opt)
@@ -19,7 +19,6 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->min_mid_occ = 10;
 	opt->max_mid_occ = 1000000;
 	opt->sdust_thres = 0; // no SDUST masking
-	opt->q_occ_frac = 0.01f;
 
 	opt->min_cnt = 3;
 	opt->min_chain_score = 40;
@@ -33,7 +32,6 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->rmq_rescue_size = 1000;
 	opt->rmq_rescue_ratio = 0.1f;
 	opt->chain_gap_scale = 0.8f;
-	opt->chain_skip_scale = 0.0f;
 	opt->max_max_occ = 4095;
 	opt->occ_dist = 500;
 
@@ -45,7 +43,6 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->alt_drop = 0.15f;
 
 	opt->a = 2, opt->b = 4, opt->q = 4, opt->e = 2, opt->q2 = 24, opt->e2 = 1;
-	opt->transition = 0;
 	opt->sc_ambi = 1;
 	opt->zdrop = 400, opt->zdrop_inv = 200;
 	opt->end_bonus = -1;
@@ -55,7 +52,6 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->max_clip_ratio = 1.0f;
 	opt->mini_batch_size = 500000000;
 	opt->max_sw_mat = 100000000;
-	opt->cap_kalloc = 500000000;
 
 	opt->rank_min_len = 500;
 	opt->rank_frac = 0.9f;
@@ -75,7 +71,6 @@ void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)
 		if (opt->max_mid_occ > opt->min_mid_occ && opt->mid_occ > opt->max_mid_occ)
 			opt->mid_occ = opt->max_mid_occ;
 	}
-	if (opt->bw_long < opt->bw) opt->bw_long = opt->bw;
 	if (mm_verbose >= 3)
 		fprintf(stderr, "[M::%s::%.3f*%.2f] mid_occ = %d\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), opt->mid_occ);
 }
@@ -91,7 +86,7 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 	if (preset == 0) {
 		mm_idxopt_init(io);
 		mm_mapopt_init(mo);
-	} else if (strcmp(preset, "lr") == 0 || strcmp(preset, "map-ont") == 0) { // this is the same as the default
+	} else if (strcmp(preset, "map-ont") == 0) { // this is the same as the default
 	} else if (strcmp(preset, "ava-ont") == 0) {
 		io->flag = 0, io->k = 15, io->w = 5;
 		mo->flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN;
@@ -99,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;
@@ -106,33 +104,19 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_chain_skip = 25;
 		mo->bw_long = mo->bw;
 		mo->occ_dist = 0;
-	} else if (strcmp(preset, "lr:hq") == 0 || strcmp(preset, "map-hifi") == 0 || strcmp(preset, "map-ccs") == 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;
+		mo->occ_dist = 500;
 		mo->min_mid_occ = 50, mo->max_mid_occ = 500;
-		if (strcmp(preset, "map-hifi") == 0 || strcmp(preset, "map-ccs") == 0) {
-			mo->a = 1, mo->b = 4, mo->q = 6, mo->q2 = 26, mo->e = 2, mo->e2 = 1;
-			mo->min_dp_max = 200;
-		}
-	} else if (strcmp(preset, "lr:hqae") == 0) { // high-quality assembly evaluation
-		io->flag = 0, io->k = 25, io->w = 51;
-		mo->flag |= MM_F_RMQ;
-		mo->min_mid_occ = 50, mo->max_mid_occ = 500;
-		mo->rmq_inner_dist = 5000;
-		mo->occ_dist = 200;
-		mo->best_n = 100;
-		mo->chain_gap_scale = 5.0f;
-	} else if (strcmp(preset, "map-iclr-prerender") == 0) {
-		io->flag = 0, io->k = 15;
-		mo->b = 6, mo->transition = 1;
-		mo->q = 10, mo->q2 = 50;
-	} else if (strcmp(preset, "map-iclr") == 0) {
-		io->flag = 0, io->k = 19;
-		mo->b = 6, mo->transition = 4;
-		mo->q = 10, mo->q2 = 50;
+		mo->min_dp_max = 200;
 	} else if (strncmp(preset, "asm", 3) == 0) {
 		io->flag = 0, io->k = 19, io->w = 19;
-		mo->bw = 1000, mo->bw_long = 100000;
+		mo->bw = mo->bw_long = 100000;
 		mo->max_gap = 10000;
 		mo->flag |= MM_F_RMQ;
 		mo->min_mid_occ = 50, mo->max_mid_occ = 500;
@@ -174,7 +158,7 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		mo->junc_bonus = 9;
 		mo->zdrop = 200, mo->zdrop_inv = 100; // because mo->a is halved
 		if (strcmp(preset, "splice:hq") == 0)
-			mo->noncan = 5, mo->b = 4, mo->q = 6, mo->q2 = 24;
+			mo->junc_bonus = 5, mo->b = 4, mo->q = 6, mo->q2 = 24;
 	} else return -1;
 	return 0;
 }

pyproject.toml

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

python/README.rst

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

python/cmappy.pxd

@@ -23,7 +23,6 @@ cdef extern from "minimap.h":
 		int min_cnt
 		int min_chain_score
 		float chain_gap_scale
-		float chain_skip_scale
 		int rmq_size_cap, rmq_inner_dist
 		int rmq_rescue_size
 		float rmq_rescue_ratio
@@ -36,7 +35,6 @@ cdef extern from "minimap.h":
 		float alt_drop
 
 		int a, b, q, e, q2, e2
-		int transition
 		int sc_ambi
 		int noncan
 		int junc_bonus
@@ -53,7 +51,6 @@ cdef extern from "minimap.h":
 		int pe_ori, pe_bonus
 
 		float mid_occ_frac
-		float q_occ_frac
 		int32_t min_mid_occ
 		int32_t mid_occ
 		int32_t max_occ

python/mappy.pyx

@@ -3,7 +3,7 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
-__version__ = '2.28'
+__version__ = '2.22'
 
 cmappy.mm_reset_timer()
 
@@ -96,7 +96,6 @@ cdef class Alignment:
 		a = [str(self._q_st), str(self._q_en), strand, self._ctg, str(self._ctg_len), str(self._r_st), str(self._r_en),
 			str(self._mlen), str(self._blen), str(self._mapq), tp, ts, "cg:Z:" + self.cigar_str]
 		if self._cs != "": a.append("cs:Z:" + self._cs)
-		if self._MD != "": a.append("MD:Z:" + self._MD)
 		return "\t".join(a)
 
 cdef class ThreadBuffer:
@@ -113,7 +112,7 @@ cdef class Aligner:
 	cdef cmappy.mm_idxopt_t idx_opt
 	cdef cmappy.mm_mapopt_t map_opt
 
-	def __cinit__(self, fn_idx_in=None, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, bw_long=None, best_n=None, n_threads=3, fn_idx_out=None, max_frag_len=None, extra_flags=None, seq=None, scoring=None):
+	def __cinit__(self, fn_idx_in=None, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, best_n=None, n_threads=3, fn_idx_out=None, max_frag_len=None, extra_flags=None, seq=None, scoring=None):
 		self._idx = NULL
 		cmappy.mm_set_opt(NULL, &self.idx_opt, &self.map_opt) # set the default options
 		if preset is not None:
@@ -126,7 +125,6 @@ cdef class Aligner:
 		if min_chain_score is not None: self.map_opt.min_chain_score = min_chain_score
 		if min_dp_score is not None: self.map_opt.min_dp_max = min_dp_score
 		if bw is not None: self.map_opt.bw = bw
-		if bw_long is not None: self.map_opt.bw_long = bw_long
 		if best_n is not None: self.map_opt.best_n = best_n
 		if max_frag_len is not None: self.map_opt.max_frag_len = max_frag_len
 		if extra_flags is not None: self.map_opt.flag |= extra_flags
@@ -174,7 +172,6 @@ cdef class Aligner:
 		cdef cmappy.mm_mapopt_t map_opt
 
 		if self._idx == NULL: return
-		if ((self.map_opt.flag & 4) and (self._idx.flag & 2)): return
 		map_opt = self.map_opt
 		if max_frag_len is not None: map_opt.max_frag_len = max_frag_len
 		if extra_flags is not None: map_opt.flag |= extra_flags
@@ -220,7 +217,6 @@ cdef class Aligner:
 		cdef int l
 		cdef char *s
 		if self._idx == NULL: return
-		if ((self.map_opt.flag & 4) and (self._idx.flag & 2)): return
 		s = cmappy.mappy_fetch_seq(self._idx, name.encode(), start, end, &l)
 		if l == 0: return None
 		r = s[:l] if isinstance(s, str) else s[:l].decode()

python/minimap2.py

@@ -5,7 +5,7 @@ import getopt
 import mappy as mp
 
 def main(argv):
-	opts, args = getopt.getopt(argv[1:], "x:n:m:k:w:r:cM")
+	opts, args = getopt.getopt(argv[1:], "x:n:m:k:w:r:c")
 	if len(args) < 2:
 		print("Usage: minimap2.py [options] <ref.fa>|<ref.mmi> <query.fq>")
 		print("Options:")
@@ -16,11 +16,10 @@ def main(argv):
 		print("  -w INT      minimizer window length")
 		print("  -r INT      band width")
 		print("  -c          output the cs tag")
-		print("  -M          output the MD tag")
 		sys.exit(1)
 
 	preset = min_cnt = min_sc = k = w = bw = None
-	out_cs = out_MD = False
+	out_cs = False
 	for opt, arg in opts:
 		if opt == '-x': preset = arg
 		elif opt == '-n': min_cnt = int(arg)
@@ -29,12 +28,11 @@ def main(argv):
 		elif opt == '-k': k = int(arg)
 		elif opt == '-w': w = int(arg)
 		elif opt == '-c': out_cs = True
-		elif opt == '-M': out_MD = True
 
 	a = mp.Aligner(args[0], preset=preset, min_cnt=min_cnt, min_chain_score=min_sc, k=k, w=w, bw=bw)
 	if not a: raise Exception("ERROR: failed to load/build index file '{}'".format(args[0]))
 	for name, seq, qual in mp.fastx_read(args[1]): # read one sequence
-		for h in a.map(seq, cs=out_cs, MD=out_MD): # traverse hits
+		for h in a.map(seq, cs=out_cs): # traverse hits
 			print('{}\t{}\t{}'.format(name, len(seq), h))
 
 if __name__ == "__main__":

seed.c

@@ -1,34 +1,40 @@
 #include "mmpriv.h"
 #include "kalloc.h"
 #include "ksort.h"
+#include <stdlib.h>
+#include<algorithm>
+#include <x86intrin.h>
 
-void mm_seed_mz_flt(void *km, mm128_v *mv, int32_t q_occ_max, float q_occ_frac)
-{
-	mm128_t *a;
-	size_t i, j, st;
-	if (mv->n <= q_occ_max || q_occ_frac <= 0.0f || q_occ_max <= 0) return;
-	a = Kmalloc(km, mm128_t, mv->n);
-	for (i = 0; i < mv->n; ++i)
-		a[i].x = mv->a[i].x, a[i].y = i;
-	radix_sort_128x(a, a + mv->n);
-	for (st = 0, i = 1; i <= mv->n; ++i) {
-		if (i == mv->n || a[i].x != a[st].x) {
-			int32_t cnt = i - st;
-			if (cnt > q_occ_max && cnt > mv->n * q_occ_frac)
-				for (j = st; j < i; ++j)
-					mv->a[a[j].y].x = 0;
-			st = i;
-		}
-	}
-	kfree(km, a);
-	for (i = j = 0; i < mv->n; ++i)
-		if (mv->a[i].x != 0)
-			mv->a[j++] = mv->a[i];
-	mv->n = j;
-}
+#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;
@@ -39,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;
@@ -47,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;
 }
 
@@ -112,8 +132,7 @@ mm_seed_t *mm_collect_matches(void *km, int *_n_m, int qlen, int max_occ, int ma
 	}
 	for (i = 0, n_m = 0, *rep_len = 0, *n_a = 0; i < n_m0; ++i) {
 		mm_seed_t *q = &m[i];
-		if (mm_dbg_flag & MM_DBG_SEED_FREQ)
-			fprintf(stderr, "SF\t%d\t%d\t%d\n", q->q_pos>>1, q->n, q->flt);
+		//fprintf(stderr, "X\t%d\t%d\t%d\n", q->q_pos>>1, q->n, q->flt);
 		if (q->flt) {
 			int en = (q->q_pos >> 1) + 1, st = en - q->q_span;
 			if (st > rep_en) {

setup.py

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

tex/minimap2.bib

@@ -370,6 +370,7 @@
 	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}
@@ -448,13 +449,3 @@
 	Title = {A synthetic-diploid benchmark for accurate variant-calling evaluation},
 	Volume = {15},
 	Year = {2018}}
-
-@article{Gu:1995wt,
-	author = {Gu, X and Li, W H},
-	journal = {J Mol Evol},
-	month = {Apr},
-	number = {4},
-	pages = {464-73},
-	title = {The size distribution of insertions and deletions in human and rodent pseudogenes suggests the logarithmic gap penalty for sequence alignment},
-	volume = {40},
-	year = {1995}}

tex/mm2-update.tex

@@ -57,8 +57,8 @@ in v2.19 through v2.22 to improve mapping results.
 \begin{methods}
 \section{Methods}
 
-\subsection{Rescuing high-occurrence $k$-mers}\label{sec:high-occ}
-Minimap2 keeps all $k$-mer minimizers~\citep{Roberts:2004fv} during indexing. Its original
+\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
@@ -66,24 +66,23 @@ chaining due to insufficient anchors.
 
 To resolve this issue, we implemented a new heuristic to add additional
 minimizers. Suppose we are looking at two adjacent low-occurence $k$-mers
-located at position $x_1$ and $x_2$, respectively. If $|x_1-x_2|\ge L$,
-minimap2 v2.22 additionally selects $\lfloor|x_1-x_2|/L\rfloor$ minimizers
-of the lowest occurrence among minimizers between $x_1$ and $x_2$. Here
-parameter $L$ controls the frequency of sampling. It defaults to 500.
+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, slow for chaining
+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, which means a gap longer than 500bp will stop chaining.
-To align through longer gaps, older minimap2 implemented a long-join heurstic as follows.
-If there is an INDEL longer than 500bp and the two chains around the INDEL
+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.
-This heuristic may fail around VNTRs because short chains
+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
@@ -91,13 +90,13 @@ duplication, the early escape may break around the event even if users
 specify a large band.
 
 In minigraph~\citep{Li:2020aa}, we developed a new chaining algorithm that
-finds up to 1kb INDELs with DP-based chaining and goes through longer INDELs with a
+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 still uses the DP-based algorithm to
+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 because it can resolve overlaps between short chains. The old
+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}
@@ -107,25 +106,25 @@ the best scoring alignment is sometimes not the correct alignment.
 \citet{Jain2020.11.01.363887} resolved this dilemma by altering the mapping
 algorithm.
 
-In our view, this problem is rooted in inapropriate scoring: affine-gap penalty
+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 by a function linear in the SV length. Minimap2 v2.22 instead rescores
+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
 $$
-S=M-\frac{N+G}{2d}-\sum_{i=1}^G\log_2(1+g_i)
+M-\frac{N+G}{2d}-\sum_{i=1}^G\log_2(1+g_i)
 $$
 where $g_i\ge1$ is the length of the $i$-th gap and
 $$
 d=\max\left\{\frac{N+G}{M+N+G},0.02\right\}
 $$
-It approximates per-base sequence divergence except with the smallest value set
+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~\citep{Gu:1995wt}. Our scoring gives a long SV
+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. The time spent on rescoring is negligible in
+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
@@ -145,15 +144,13 @@ practice.
 \toprule
 $[$Benchmark$]$ Metric & v2.22 & v2.18 & Winno & lra \\
 \midrule
-$[$sim-map$]$ \% mapped reads at Q10      & 97.9       & 97.6      & {\bf 99.0}& 97.3 \\
-$[$sim-map$]$ err. rate at Q10 (phredQ)   & {\bf 52}   & {\bf 52}  & 38        & 24 \\
-$[$winno-cmp$]$ rate of diff. (phredQ)    & {\bf 41}   & 37        & truth     & 18 \\
-$[$winno-cmp$]$ CPU time (hour)           & {\bf 5.0}  & 5.3       & 71.8      & 13.1 \\
-$[$winno-cmp$]$ peak RAM (Gb)             & 17.1       & 14.4      & {\bf 9.6} & 12.4 \\
-$[$sim-sv$]$  \% false negative rate      & {\bf 0.5}  & 2.0       & {\bf 0.5} & 1.4  \\
-$[$sim-sv$]$  \% false discovery rate     & {\bf 0.0}  & 0.1       & {\bf 0.0} & 0.1  \\
-$[$real-sv-1k$]$ \% false negative rate   & {\bf 7.3}  & 20.0      & 13.0      & N/A \\
-$[$real-sv-1k$]$ \% false discovery rate  & 2.7        & {\bf 2.4} & 2.7       & N/A \\
+$[$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
@@ -162,11 +159,11 @@ $[$real-sv-1k$]$ \% false discovery rate  & 2.7        & {\bf 2.4} & 2.7       &
 10 or higher were evaluated by ``paftools.js mapeval''. The mapping error rate
 is measured in the phred scale: if the error rate is $e$, $-10\log_{10}e$ is
 reported in the table. In $[$winno-cmp$]$, 1.39 million CHM13 HiFi reads from
-SRR11292121 were mapped against the same CHM13 assembly. 99.3\% of them were mapped by Winnowmap2
+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-fold coverage with the same pbsim2 command line. SVs were called with
+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
@@ -175,44 +172,31 @@ GCA\_018852615.1) and compared to the GIAB truth~\citep{Zook:2020aa} using ``tru
 \end{table}
 
 We evaluated minimap2 v2.22 along with v2.18, Winnowmap2 v2.03 and lra v1.3.2
-(Table~\ref{tab:1}), using the default setting of each mapper according to the input data types.
-Both versions of minimap2 achieved high mapping accuracy on
+(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, we took Winnowmap2 mapping as ground
-truth to evaluate other mappers (winno-cmp in Table~\ref{tab:1}). Out of 1,378,092 reads with mapQ10
+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.
-Minimap2 v2.18 is less consistent with 275 differences including 30 unmapped
-reads mappable by both Winnowmap2 and v2.22.
+highly similar but not identical alignments. We are not sure what are real
+mapping errors.
 
-For the minimizer rescuing parameter $L$ in Section~\ref{sec:high-occ},
-we set its default to 500 such that v2.22 has comparable performance to v2.18 given simulated PacBio and Nanopore human reads.
-To see the effect of this parameter on real data, we tried several different $L$ values.
-v2.22 gave 99 mapping differences at $L=200$,
-118 at $L=500$ (default), 167 at $L=750$ and 224 differences at $L=1000$ in comparison to Winnowmap2.
-$L=200$ is 28\% slower than the default while $L=1000$ is 9\% faster.
-Changing the default minimizer window size (option ``-w'')
-and the initial minimizer occurrence cutoff (option ``-f'')
-also affects performance and accuracy to a similar magnitude.
-
-The two benchmarks above only evaluate read mappings when there are no variations between the reads and the reference.
+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 five SVs (accounting for
-0.5\% of simulated SVs). For four out of these five missing SVs, minimap2 v2.22
-mapped more variant reads than Winnowmap2. Sniffles did not call these SVs
-because minimap2 tended to give them conservative mapping quality. It is worth
-noting that the simulation here only considers a simple scenario in evolution.
-Non-allelic gene conversions, which happen often in segmental
+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.
@@ -220,17 +204,18 @@ remains a question.
 To see if minimap2 v2.22 could improve long INDEL alignment, we ran dipcall on
 contig-to-reference alignments and focused on INDELs longer than 1kb
 (real-sv-1k). v2.22 is more sensitive at comparable specificity, confirming its
-advantage in more contiguous alignment. We could not get dipcall to work well with lra,
-so did not report the numbers.
+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 versions. It is
+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 for which
+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.