Estimated installation time

Based on DOI:10.1101/2021.01.15.426881. One minimap2 reviewer suggested the
right formula to me but I thought the difference would be insignificant. I was
wrong.

.gitmodules

@@ -0,0 +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

@@ -6,6 +6,10 @@ matrix:
     - 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
@@ -15,6 +19,6 @@ matrix:
       before_install: pip install cython
       script: python setup.py build_ext
     - language: python
-      python: "3.6"
+      python: "3.9"
       before_install: pip install cython
       script: python setup.py build_ext

Makefile

@@ -1,14 +1,74 @@
-CFLAGS=		-g -Wall -O2 -Wc++-compat #-Wextra
+
-CPPFLAGS=	-DHAVE_KALLOC
+## /* The MIT License
-INCLUDES=
+## 
+## 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>
+## */
+## 
+CFLAGS=		-Wall -O2 -Wc++-compat #-Wextra
+CPPFLAGS=	-DHAVE_KALLOC -march=native
+
+OPT_FLAGS= -DVECTORIZED_CHAINING -DALIGN_AVX  
+ 
+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=	-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 chain.o align.o hit.o map.o format.o pe.o esterr.o splitidx.o ksw2_ll_sse.o
 PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite
-LIBS=		-lm -lz -lpthread
+LIBS=		-lm -lz -lpthread 
+
+CC=$(CXX)
+ifeq ($(CC), g++)
+	CC=g++ -std=c++11
+endif
 
 ifeq ($(arm_neon),) # if arm_neon is not defined
 ifeq ($(sse2only),) # if sse2only is not defined
-	OBJS+=ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o
+	OBJS+=ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o ksw2_extd2_avx.o
 else                # if sse2only is defined
 	OBJS+=ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_exts2_sse.o
 endif
@@ -54,6 +114,17 @@ libminimap2.a:$(OBJS)
 sdust:sdust.c kalloc.o kalloc.h kdq.h kvec.h kseq.h ketopt.h sdust.h
 		$(CC) -D_SDUST_MAIN $(CFLAGS) $< kalloc.o -o $@ -lz
 
+multi:
+		$(MAKE) clean
+		$(MAKE)
+		mv minimap2 mm2-fast
+		$(MAKE) clean
+		$(MAKE) lhash=1
+		mv minimap2 mm2-fast-lhash
+		$(MAKE) clean
+		$(MAKE) no_opt=1
+		mv minimap2 mm2-fast-no-opt
+
 # SSE-specific targets on x86/x86_64
 
 ifeq ($(arm_neon),)   # if arm_neon is defined, compile this target with the default setting (i.e. no -msse2)

Makefile.simde

@@ -0,0 +1,97 @@
+CFLAGS=		-g -Wall -O2 -Wc++-compat #-Wextra
+CPPFLAGS=	-DHAVE_KALLOC -DUSE_SIMDE -DSIMDE_ENABLE_NATIVE_ALIASES
+INCLUDES=	-Ilib/simde
+OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o chain.o align.o hit.o map.o format.o pe.o esterr.o splitidx.o \
+			ksw2_extz2_simde.o ksw2_extd2_simde.o ksw2_exts2_simde.o ksw2_ll_simde.o
+PROG=		minimap2
+PROG_EXTRA=	sdust minimap2-lite
+LIBS=		-lm -lz -lpthread
+
+
+ifneq ($(arm_neon),) # if arm_neon is defined
+ifeq ($(aarch64),)   #if aarch64 is not defined
+	CFLAGS+=-D_FILE_OFFSET_BITS=64 -mfpu=neon -fsigned-char
+else                 #if aarch64 is defined
+	CFLAGS+=-D_FILE_OFFSET_BITS=64 -fsigned-char
+endif
+endif
+
+ifneq ($(asan),)
+	CFLAGS+=-fsanitize=address
+	LIBS+=-fsanitize=address
+endif
+
+ifneq ($(tsan),)
+	CFLAGS+=-fsanitize=thread
+	LIBS+=-fsanitize=thread
+endif
+
+.PHONY:all extra clean depend
+.SUFFIXES:.c .o
+
+.c.o:
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@
+
+all:$(PROG)
+
+extra:all $(PROG_EXTRA)
+
+minimap2:main.o libminimap2.a
+		$(CC) $(CFLAGS) main.o -o $@ -L. -lminimap2 $(LIBS)
+
+minimap2-lite:example.o libminimap2.a
+		$(CC) $(CFLAGS) $< -o $@ -L. -lminimap2 $(LIBS)
+
+libminimap2.a:$(OBJS)
+		$(AR) -csru $@ $(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
+
+ksw2_ll_simde.o:ksw2_ll_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) -msse2 $(CPPFLAGS) $(INCLUDES) $< -o $@
+
+ksw2_extz2_simde.o:ksw2_extz2_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) -msse4.1 $(CPPFLAGS) $(INCLUDES) $< -o $@
+
+ksw2_extd2_simde.o:ksw2_extd2_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) -msse4.1 $(CPPFLAGS) $(INCLUDES) $< -o $@
+
+ksw2_exts2_simde.o:ksw2_exts2_sse.c ksw2.h kalloc.h
+		$(CC) -c $(CFLAGS) -msse4.1 $(CPPFLAGS) $(INCLUDES) $< -o $@
+
+# other non-file targets
+
+clean:
+		rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist mappy*.so mappy.c python/mappy.c mappy.egg*
+
+depend:
+		(LC_ALL=C; export LC_ALL; makedepend -Y -- $(CFLAGS) $(CPPFLAGS) -- *.c)
+
+# DO NOT DELETE
+
+align.o: minimap.h mmpriv.h bseq.h kseq.h ksw2.h kalloc.h
+bseq.o: bseq.h kvec.h kalloc.h kseq.h
+chain.o: minimap.h mmpriv.h bseq.h kseq.h kalloc.h
+esterr.o: mmpriv.h minimap.h bseq.h kseq.h
+example.o: minimap.h kseq.h
+format.o: kalloc.h mmpriv.h minimap.h bseq.h kseq.h
+hit.o: mmpriv.h minimap.h bseq.h kseq.h kalloc.h khash.h
+index.o: kthread.h bseq.h minimap.h mmpriv.h kseq.h kvec.h kalloc.h khash.h
+index.o: ksort.h
+kalloc.o: kalloc.h
+ksw2_extd2_sse.o: ksw2.h kalloc.h
+ksw2_exts2_sse.o: ksw2.h kalloc.h
+ksw2_extz2_sse.o: ksw2.h kalloc.h
+ksw2_ll_sse.o: ksw2.h kalloc.h
+kthread.o: kthread.h
+main.o: bseq.h minimap.h mmpriv.h kseq.h ketopt.h
+map.o: kthread.h kvec.h kalloc.h sdust.h mmpriv.h minimap.h bseq.h kseq.h
+map.o: khash.h ksort.h
+misc.o: mmpriv.h minimap.h bseq.h kseq.h ksort.h
+options.o: mmpriv.h minimap.h bseq.h kseq.h
+pe.o: mmpriv.h minimap.h bseq.h kseq.h kvec.h kalloc.h ksort.h
+sdust.o: kalloc.h kdq.h kvec.h sdust.h
+self-chain.o: minimap.h kseq.h
+sketch.o: kvec.h kalloc.h mmpriv.h minimap.h bseq.h kseq.h
+splitidx.o: mmpriv.h minimap.h bseq.h kseq.h

NEWS.md

@@ -1,3 +1,89 @@
+Release 2.18-r1015 (9 April 2021)
+---------------------------------
+
+This release fixes multiple rare bugs in minimap2 and adds additional
+functionality to paftools.js.
+
+Changes to minimap2:
+
+ * Bugfix: a rare segfault caused by an off-by-one error (#489)
+
+ * Bugfix: minimap2 segfaulted due to an uninitilized variable (#622 and #625).
+
+ * Bugfix: minimap2 parsed spaces as field separators in BED (#721). This led
+   to issues when the BED name column contains spaces.
+
+ * Bugfix: minimap2 `--split-prefix` did not work with long reference names
+   (#394).
+
+ * Bugfix: option `--junc-bonus` didn't work (#513)
+
+ * Bugfix: minimap2 didn't return 1 on I/O errors (#532)
+
+ * Bugfix: the `de:f` tag (sequence divergence) could be negative if there were
+   ambiguous bases
+
+ * Bugfix: fixed two undefined behaviors caused by calling memcpy() on
+   zero-length blocks (#443)
+
+ * Bugfix: there were duplicated SAM @SQ lines if option `--split-prefix` is in
+   use (#400 and #527)
+
+ * Bugfix: option -K had to be smaller than 2 billion (#491). This was caused
+   by a 32-bit integer overflow.
+
+ * Improvement: optionally compile against SIMDe (#597). Minimap2 should work
+   with IBM POWER CPUs, though this has not been tested. To compile with SIMDe,
+   please use `make -f Makefile.simde`.
+
+ * Improvement: more informative error message for I/O errors (#454) and for
+   FASTQ parsing errors (#510)
+
+ * Improvement: abort given malformatted RG line (#541)
+
+ * Improvement: better formula to estimate the `dv:f` tag (approximate sequence
+   divergence). See DOI:10.1101/2021.01.15.426881.
+
+ * New feature: added the `--mask-len` option to fine control the removal of
+   redundant hits (#659). The default behavior is unchanged.
+
+Changes to mappy:
+
+ * Bugfix: mappy caused segmentation fault if the reference index is not
+   present (#413).
+
+ * Bugfix: fixed a memory leak via 238b6bb3
+
+ * Change: always require Cython to compile the mappy module (#723). Older
+   mappy packages at PyPI bundled the C source code generated by Cython such
+   that end users did not need to install Cython to compile mappy. However, as
+   Python 3.9 is breaking backward compatibility, older mappy does not work
+   with Python 3.9 anymore. We have to add this Cython dependency as a
+   workaround.
+
+Changes to paftools.js:
+
+ * Bugfix: the "part10-" line from asmgene was wrong (#581)
+
+ * Improvement: compatibility with GTF files from GenBank (#422)
+
+ * New feature: asmgene also checks missing multi-copy genes
+
+ * New feature: added the misjoin command to evaluate large-scale misjoins and
+   megabase-long inversions.
+
+Although given the many bug fixes and minor improvements, the core algorithm
+stays the same. This version of minimap2 produces nearly identical alignments
+to v2.17 except very rare corner cases.
+
+Now unimap is recommended over minimap2 for aligning long contigs against a
+reference genome. It often takes less wall-clock time and is much more
+sensitive to long insertions and deletions.
+
+(2.18: 9 April 2021, r1015)
+
+
+
 Release 2.17-r941 (4 May 2019)
 ------------------------------
 

README.md

@@ -1,3 +1,73 @@
+## 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 3.5x speedup 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** 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](https://en.wikipedia.org/wiki/AVX-512)   
+Memory requirement: ~30GB for human genome   
+
+### Installation
+Clone the *fast-contrib* branch from minimap2 github page. The source code can be compiled by simple using *make* command. It only takes a few seconds.
+```
+git clone --recursive https://github.com/lh3/minimap2.git -b fast-contrib 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.18, the output of mm2-fast can be verified against minimap2-v2.18. Note that AVX512-based chaining in mm2-fast by default runs with a chaining parameter *max-skip=infinity* for higher chaining precision. Therefore, for correctness verification, minimap2 should run with a larger value of *max-skip* parameter. Follow the below steps to verify the accuracy of mm2-fast. 
+```sh
+git clone https://github.com/lh3/minimap2.git -b v2.18
+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 show a clean-diff with the difference of 2 lines, i.e., the lines containing the command-line parameters for minimap2 and mm2-fast.
+
+### Advanced options
+The default compilation using make applies two optimizations: AVX512 vectorized chaining and alignment, and 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 20-30 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
+```
+mm2-fast includes preliminary support for AVX2 architecture. Currently, chaining step is not optimized for AVX2 but the seeding and alignment steps are available. To try mm2-fast on AVX2 systems, use the following command to compile.
+```sh
+make clean && make lhash=1 use_avx2=1
+```
+### Performance
+We have observed up to 3.5x 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 80 seconds, while mm2-fast takes 38 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
+The current version of mm2-fast is based on minimap2-v2.18. We are planning to apply our optimizations to minimap2 master branch.
+### 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)
@@ -19,12 +89,17 @@ cd minimap2 && make
 ./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
 ./minimap2 -ax splice:hq -uf ref.fa query.fa > aln.sam    # Final PacBio Iso-seq or traditional cDNA
+./minimap2 -ax splice --junc-bed anno.bed12 ref.fa query.fa > aln.sam  # prioritize on annotated junctions
 ./minimap2 -cx asm5 asm1.fa asm2.fa > aln.paf             # intra-species asm-to-asm alignment
 ./minimap2 -x ava-pb reads.fa reads.fa > overlaps.paf     # PacBio read overlap
 ./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf    # Nanopore read overlap
 # man page for detailed command line options
 man ./minimap2.1
 ```
+[Unimap][unimap] is recommended for aligning long contigs against a reference
+genome. It often takes less wall-clock time and is much more sensitive to long
+insertions and deletions.
+
 ## Table of Contents
 
 - [Getting Started](#started)
@@ -71,8 +146,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.17/minimap2-2.17_x64-linux.tar.bz2 | tar -jxvf -
+curl -L https://github.com/lh3/minimap2/releases/download/v2.18/minimap2-2.18_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.17_x64-linux/minimap2
+./minimap2-2.18_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
@@ -80,7 +155,14 @@ directory to compile. If you see compilation errors, try `make sse2only=1`
 to disable SSE4 code, which will make minimap2 slightly slower.
 
 Minimap2 also works with ARM CPUs supporting the NEON instruction sets. To
-compile for 32 bit ARM architectures (such as ARMv7), use `make arm_neon=1`. To compile for for 64 bit ARM architectures (such as ARMv8), use `make arm_neon=1 aarch64=1`.
+compile for 32 bit ARM architectures (such as ARMv7), use `make arm_neon=1`. To
+compile for for 64 bit ARM architectures (such as ARMv8), use `make arm_neon=1
+aarch64=1`.
+
+Minimap2 can use [SIMD Everywhere (SIMDe)][simde] library for porting
+implementation to the different SIMD instruction sets. To compile using SIMDe,
+use `make -f Makefile.simde`. To compile for ARM CPUs, use `Makefile.simde`
+with the ARM related command lines given above.
 
 ### <a name="general"></a>General usage
 
@@ -178,6 +260,19 @@ This is because SIRV does not honor the evolutionarily conservative splicing
 signal. If you are studying SIRV, you may apply `--splice-flank=no` to let
 minimap2 only model GT..AG, ignoring the additional base.
 
+Since v2.17, minimap2 can optionally take annotated genes as input and
+prioritize on annotated splice junctions. To use this feature, you can 
+```sh
+paftools.js gff2bed anno.gff > anno.bed
+minimap2 -ax splice --junc-bed anno.bed ref.fa query.fa > aln.sam
+```
+Here, `anno.gff` is the gene annotation in the GTF or GFF3 format (`gff2bed`
+automatically tests the format). The output of `gff2bed` is in the 12-column
+BED format, or the BED12 format. With the `--junc-bed` option, minimap2 adds a
+bonus score (tuned by `--junc-bonus`) if an aligned junction matches a junction
+in the annotation. Option `--junc-bed` also takes 5-column BED, including the
+strand field. In this case, each line indicates an oriented junction.
+
 #### <a name="long-overlap"></a>Find overlaps between long reads
 
 ```sh
@@ -376,3 +471,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
+[smide]: https://github.com/nemequ/simde
+[unimap]: https://github.com/lh3/unimap

align.c

@@ -1,3 +1,33 @@
+/* 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 <assert.h>
 #include <string.h>
 #include <stdlib.h>
@@ -5,7 +35,9 @@
 #include "minimap.h"
 #include "mmpriv.h"
 #include "ksw2.h"
-
+#include "ksw2_extd2_avx.h"
+#include <x86intrin.h>
+extern uint64_t alignment_time;
 static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t sc_ambi)
 {
 	int i, j;
@@ -38,8 +70,8 @@ static inline void update_max_zdrop(int32_t score, int i, int j, int32_t *max, i
 		int z = *max - score - diff * e;
 		if (z > *max_zdrop) {
 			*max_zdrop = z;
-			pos[0][0] = *max_i, pos[0][1] = i + 1;
+			pos[0][0] = *max_i, pos[0][1] = i;
-			pos[1][0] = *max_j, pos[1][1] = j + 1;
+			pos[1][0] = *max_j, pos[1][1] = j;
 		}
 	} else *max = score, *max_i = i, *max_j = j;
 }
@@ -312,6 +344,10 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) //
 
 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);
@@ -327,8 +363,18 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 		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
+	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__
+		ksw_extd2_avx2(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez);
+#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);
@@ -336,6 +382,9 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 			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
 }
 
 static inline int mm_get_hplen_back(const mm_idx_t *mi, uint32_t rid, uint32_t x)
@@ -739,6 +788,13 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			if (ez->n_cigar > 0)
 				mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			if (ez->zdropped) { // truncated by Z-drop; TODO: sometimes Z-drop kicks in because the next seed placement is wrong. This can be fixed in principle.
+				if (!r->p) {
+					assert(ez->n_cigar == 0);
+					uint32_t capacity = sizeof(mm_extra_t)/4;
+					kroundup32(capacity);
+					r->p = (mm_extra_t*)calloc(capacity, 4);
+					r->p->capacity = capacity;
+				}
 				for (j = i - 1; j >= 0; --j)
 					if ((int32_t)a[as1 + j].x <= rs + ez->max_t)
 						break;
@@ -908,6 +964,6 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 	kfree(km, qseq0[0]);
 	kfree(km, ez.cigar);
 	mm_filter_regs(opt, qlen, n_regs_, regs);
-	mm_hit_sort(km, n_regs_, regs);
+	mm_hit_sort(km, n_regs_, regs, opt->alt_drop);
 	return regs;
 }

bseq.c

@@ -77,7 +77,7 @@ static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual, int with_
 	s->l_seq = ks->seq.l;
 }
 
-mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int with_comment, int frag_mode, int *n_)
+mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int64_t chunk_size, int with_qual, int with_comment, int frag_mode, int *n_)
 {
 	int64_t size = 0;
 	int ret;
@@ -99,7 +99,7 @@ mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
 		size += s->l_seq;
 		if (size >= chunk_size) {
 			if (frag_mode && a.a[a.n-1].l_seq < CHECK_PAIR_THRES) {
-				while (kseq_read(ks) >= 0) {
+				while ((ret = kseq_read(ks)) >= 0) {
 					kseq2bseq(ks, &fp->s, with_qual, with_comment);
 					if (mm_qname_same(fp->s.name, a.a[a.n-1].name)) {
 						kv_push(mm_bseq1_t, 0, a, fp->s);
@@ -110,23 +110,25 @@ mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
 			break;
 		}
 	}
-	if (ret < -1)
+	if (ret < -1) {
-		fprintf(stderr, "[WARNING]\033[1;31m wrong FASTA/FASTQ record. Continue anyway.\033[0m\n");
+		if (a.n) fprintf(stderr, "[WARNING]\033[1;31m failed to parse the FASTA/FASTQ record next to '%s'. Continue anyway.\033[0m\n", a.a[a.n-1].name);
+		else fprintf(stderr, "[WARNING]\033[1;31m failed to parse the first FASTA/FASTQ record. Continue anyway.\033[0m\n");
+	}
 	*n_ = a.n;
 	return a.a;
 }
 
-mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int chunk_size, int with_qual, int frag_mode, int *n_)
+mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int64_t chunk_size, int with_qual, int frag_mode, int *n_)
 {
 	return mm_bseq_read3(fp, chunk_size, with_qual, 0, frag_mode, n_);
 }
 
-mm_bseq1_t *mm_bseq_read(mm_bseq_file_t *fp, int chunk_size, int with_qual, int *n_)
+mm_bseq1_t *mm_bseq_read(mm_bseq_file_t *fp, int64_t chunk_size, int with_qual, int *n_)
 {
 	return mm_bseq_read2(fp, chunk_size, with_qual, 0, n_);
 }
 
-mm_bseq1_t *mm_bseq_read_frag2(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int with_comment, int *n_)
+mm_bseq1_t *mm_bseq_read_frag2(int n_fp, mm_bseq_file_t **fp, int64_t chunk_size, int with_qual, int with_comment, int *n_)
 {
 	int i;
 	int64_t size = 0;
@@ -156,7 +158,7 @@ mm_bseq1_t *mm_bseq_read_frag2(int n_fp, mm_bseq_file_t **fp, int chunk_size, in
 	return a.a;
 }
 
-mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int *n_)
+mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int64_t chunk_size, int with_qual, int *n_)
 {
 	return mm_bseq_read_frag2(n_fp, fp, chunk_size, with_qual, 0, n_);
 }

bseq.h

@@ -18,11 +18,11 @@ typedef struct {
 
 mm_bseq_file_t *mm_bseq_open(const char *fn);
 void mm_bseq_close(mm_bseq_file_t *fp);
-mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int with_comment, int frag_mode, int *n_);
+mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int64_t chunk_size, int with_qual, int with_comment, int frag_mode, int *n_);
-mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int chunk_size, int with_qual, int frag_mode, int *n_);
+mm_bseq1_t *mm_bseq_read2(mm_bseq_file_t *fp, int64_t chunk_size, int with_qual, int frag_mode, int *n_);
-mm_bseq1_t *mm_bseq_read(mm_bseq_file_t *fp, int chunk_size, int with_qual, int *n_);
+mm_bseq1_t *mm_bseq_read(mm_bseq_file_t *fp, int64_t chunk_size, int with_qual, int *n_);
-mm_bseq1_t *mm_bseq_read_frag2(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int with_comment, int *n_);
+mm_bseq1_t *mm_bseq_read_frag2(int n_fp, mm_bseq_file_t **fp, int64_t chunk_size, int with_qual, int with_comment, int *n_);
-mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int chunk_size, int with_qual, int *n_);
+mm_bseq1_t *mm_bseq_read_frag(int n_fp, mm_bseq_file_t **fp, int64_t chunk_size, int with_qual, int *n_);
 int mm_bseq_eof(mm_bseq_file_t *fp);
 
 extern unsigned char seq_nt4_table[256];

build_rmi.sh

@@ -0,0 +1,16 @@
+ref_data=$1
+preset=$2
+
+make clean && make no_opt=1
+touch temp_read.fastq
+./minimap2  -ax $2 $1 temp_read.fastq -Z 1 >/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

chain.c

@@ -1,3 +1,32 @@
+/* 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 <stdint.h>
 #include <string.h>
 #include <stdio.h>
@@ -5,6 +34,10 @@
 #include "mmpriv.h"
 #include "kalloc.h"
 
+#if defined(VECTORIZED_CHAINING) && defined(__AVX512BW__)
+#include "parallel_chaining_32_bit.h"
+#endif
+
 static const char LogTable256[256] = {
 #define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
 	-1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
@@ -19,12 +52,13 @@ static inline int ilog2_32(uint32_t v)
 	return (t = v>>8) ? 8 + LogTable256[t] : LogTable256[v];
 }
 
-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, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km)
+
+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)
 { // TODO: make sure this works when n has more than 32 bits
-	int32_t k, *f, *p, *t, *v, n_u, n_v;
+	int32_t k, *p, *t, *v, n_u, n_v;
-	int64_t i, j, st = 0;
+	uint32_t *f;
-	uint64_t *u, *u2, sum_qspan = 0;
+	int64_t i, j;
-	float avg_qspan;
+	uint64_t *u, *u2;
 	mm128_t *b, *w;
 
 	if (_u) *_u = 0, *n_u_ = 0;
@@ -32,12 +66,40 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 		kfree(km, a);
 		return 0;
 	}
-	f = (int32_t*)kmalloc(km, n * 4);
+	f = (uint32_t*)kmalloc(km, n * 4);
 	p = (int32_t*)kmalloc(km, n * 4);
 	t = (int32_t*)kmalloc(km, n * 4);
 	v = (int32_t*)kmalloc(km, n * 4);
 	memset(t, 0, n * 4);
+#if defined(VECTORIZED_CHAINING) && defined(__AVX512BW__)
+	/* Allocation for debugging 
+	f_avx = (uint32_t*)kmalloc(km, n * 4);
+	p_avx = (int32_t*)kmalloc(km, n * 4);
+	*/
+	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, gap_scale, is_cdna, n_segs);
+	obj.mm_dp_vectorized(n, &anchors[0], anchor_r, anchor_q, anchor_l, f, p, v, max_dist_x, max_dist_y, NULL, NULL);	
+
+	// -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);
+#else
+	int64_t st = 0;
+	uint64_t sum_qspan = 0; 
+	float avg_qspan;
 	for (i = 0; i < n; ++i) sum_qspan += a[i].y>>32&0xff;
 	avg_qspan = (float)sum_qspan / n;
 
@@ -52,7 +114,7 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 		if (i - st > max_iter) st = i - max_iter;
 		for (j = i - 1; j >= st; --j) {
 			int64_t dr = ri - a[j].x;
-			int32_t dq = qi - (int32_t)a[j].y, dd, sc, log_dd;
+			int32_t dq = qi - (int32_t)a[j].y, dd, sc, log_dd, gap_cost;
 			int32_t sidj = (a[j].y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
 			if ((sidi == sidj && dr == 0) || dq <= 0) continue; // don't skip if an anchor is used by multiple segments; see below
 			if ((sidi == sidj && dq > max_dist_y) || dq > max_dist_x) continue;
@@ -62,14 +124,16 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 			min_d = dq < dr? dq : dr;
 			sc = min_d > q_span? q_span : dq < dr? dq : dr;
 			log_dd = dd? ilog2_32(dd) : 0;
+			gap_cost = 0;
 			if (is_cdna || sidi != sidj) {
 				int c_log, c_lin;
 				c_lin = (int)(dd * .01 * avg_qspan);
 				c_log = log_dd;
 				if (sidi != sidj && dr == 0) ++sc; // possibly due to overlapping paired ends; give a minor bonus
-				else if (dr > dq || sidi != sidj) sc -= c_lin < c_log? c_lin : c_log;
+				else if (dr > dq || sidi != sidj) gap_cost = c_lin < c_log? c_lin : c_log;
-				else sc -= c_lin + (c_log>>1);
+				else gap_cost = c_lin + (c_log>>1);
-			} else sc -= (int)(dd * .01 * avg_qspan) + (log_dd>>1);
+			} else gap_cost = (int)(dd * .01 * avg_qspan) + (log_dd>>1);
+			sc -= (int)((double)gap_cost * gap_scale + .499);
 			sc += f[j];
 			if (sc > max_f) {
 				max_f = sc, max_j = j;
@@ -83,7 +147,46 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 		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 0
+	for (i = 0; i < n; ++i) {
+		assert(f[i] == f_avx[i] && p[i] == p_avx[i]);
+	
+		//if(! (f[i] == f_avx[i] && p[i] == p_avx[i]))
+		{
+			#if 0
+			fprintf(stderr, "mm2-score:\n");
+			for (int itt = 0; itt < n; ++itt) {
+				fprintf(stderr, "%ld %ld \n", f[itt], p[itt]);
+			}
+			fprintf(stderr, "mm2-simd-score:\n");
+			for (int itt = 0; itt < n; ++itt) {
+				fprintf(stderr, "%ld %ld \n", f_avx[itt], p_avx[itt]);
+			}
+			fprintf(stderr, "anchors:\n");
+			fprintf(stderr, "%lld\n", n);
+			for (int itt = 0; itt < n; ++itt) {
+				uint64_t ri = a[itt].x;
+				int32_t qi = (int32_t)a[itt].y, q_span = a[itt].y>>32&0xff; // NB: only 8 bits of span is used!!!
+				fprintf(stderr, "%llu %ld %ld\n", ri, qi, q_span);
+			}
+			//exit(0);
+			#endif
+		}
+	
+	}
+	#if 0
+	fprintf(stderr, "%llu\n", n);
+	for (int itt = 0; itt < n; ++itt) {
+		uint64_t ri = a[itt].x;
+		int32_t qi = (int32_t)a[itt].y, q_span = a[itt].y>>32&0xff; // NB: only 8 bits of span is used!!!
+		fprintf(stderr, "%llu %ld %ld\n", ri, qi, q_span);
+	}
 
+	#endif
+	kfree(km, f_avx); kfree(km, p_avx);
+	#endif
+#endif
 	// find the ending positions of chains
 	memset(t, 0, n * 4);
 	for (i = 0; i < n; ++i)

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

esterr.c

@@ -59,6 +59,6 @@ void mm_est_err(const mm_idx_t *mi, int qlen, int n_regs, mm_reg1_t *regs, const
 		n_tot = en - st + 1;
 		if (r->qs > avg_k && r->rs > avg_k) ++n_tot;
 		if (qlen - r->qs > avg_k && l_ref - r->re > avg_k) ++n_tot;
-		r->div = logf((float)n_tot / n_match) / avg_k;
+		r->div = n_match >= n_tot? 0.0f : (float)(1.0 - pow((double)n_match / n_tot, 1.0 / avg_k));
 	}
 }

format.c

@@ -274,7 +274,7 @@ double mm_event_identity(const mm_reg1_t *r)
 		if (op == 1 || op == 2)
 			++n_gapo, n_gap += len;
 	}
-	return (double)r->mlen / (r->blen - n_gap + n_gapo);
+	return (double)r->mlen / (r->blen + r->p->n_ambi - n_gap + n_gapo);
 }
 
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
@@ -392,7 +392,7 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 {
 	const int max_bam_cigar_op = 65535;
 	int flag, n_regs = n_regss[seg_idx], cigar_in_tag = 0;
-	int this_rid = -1, this_pos = -1, this_rev = 0;
+	int this_rid = -1, this_pos = -1;
 	const mm_reg1_t *regs = regss[seg_idx], *r_prev = NULL, *r_next;
 	const mm_reg1_t *r = n_regs > 0 && reg_idx < n_regs && reg_idx >= 0? &regs[reg_idx] : NULL;
 
@@ -441,7 +441,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%s\t%d\t0\t*", mi->seq[this_rid].name, this_pos+1);
 		} else mm_sprintf_lite(s, "\t*\t0\t0\t*");
 	} else {
-		this_rid = r->rid, this_pos = r->rs, this_rev = r->rev;
+		this_rid = r->rid, this_pos = r->rs;
 		mm_sprintf_lite(s, "\t%s\t%d\t%d\t", mi->seq[r->rid].name, r->rs+1, r->mapq);
 		if ((opt_flag & MM_F_LONG_CIGAR) && r->p && r->p->n_cigar > max_bam_cigar_op - 2) {
 			int n_cigar = r->p->n_cigar;

hit.c

@@ -87,6 +87,22 @@ mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u,
 	return r;
 }
 
+void mm_mark_alt(const mm_idx_t *mi, int n, mm_reg1_t *r)
+{
+	int i;
+	if (mi->n_alt == 0) return;
+	for (i = 0; i < n; ++i)
+		if (mi->seq[r[i].rid].is_alt)
+			r[i].is_alt = 1;
+}
+
+static inline int mm_alt_score(int score, float alt_diff_frac)
+{
+	if (score < 0) return score;
+	score = (int)(score * (1.0 - alt_diff_frac) + .499);
+	return score > 0? score : 1;
+}
+
 void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a)
 {
 	if (n <= 0 || n >= r->cnt) return;
@@ -106,7 +122,7 @@ void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a)
 	r->split |= 1, r2->split |= 2;
 }
 
-void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff, int hard_mask_level) // and compute mm_reg1_t::subsc
+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) // and compute mm_reg1_t::subsc
 {
 	int i, j, k, *w;
 	uint64_t *cov;
@@ -146,13 +162,16 @@ skip_uncov:
 			min = ej - sj < ei - si? ej - sj : ei - si;
 			max = ej - sj > ei - si? ej - sj : ei - si;
 			ol = si < sj? (ei < sj? 0 : ei < ej? ei - sj : ej - sj) : (ej < si? 0 : ej < ei? ej - si : ei - si); // overlap length; TODO: this can be simplified
-			if ((float)ol / min - (float)uncov_len / max > mask_level) {
+			if ((float)ol / min - (float)uncov_len / max > mask_level && uncov_len <= mask_len) { // then this is a secondary hit
-				int cnt_sub = 0;
+				int cnt_sub = 0, sci = ri->score;
 				ri->parent = rp->parent;
-				rp->subsc = rp->subsc > ri->score? rp->subsc : ri->score;
+				if (!rp->is_alt && ri->is_alt) sci = mm_alt_score(sci, alt_diff_frac);
+				rp->subsc = rp->subsc > sci? rp->subsc : sci;
 				if (ri->cnt >= rp->cnt) cnt_sub = 1;
 				if (rp->p && ri->p && (rp->rid != ri->rid || rp->rs != ri->rs || rp->re != ri->re || ol != min)) { // the last condition excludes identical hits after DP
-					rp->p->dp_max2 = rp->p->dp_max2 > ri->p->dp_max? rp->p->dp_max2 : ri->p->dp_max;
+					sci = ri->p->dp_max;
+					if (!rp->is_alt && ri->is_alt) sci = mm_alt_score(sci, alt_diff_frac);
+					rp->p->dp_max2 = rp->p->dp_max2 > sci? rp->p->dp_max2 : sci;
 					if (rp->p->dp_max - ri->p->dp_max <= sub_diff) cnt_sub = 1;
 				}
 				if (cnt_sub) ++rp->n_sub;
@@ -166,7 +185,7 @@ set_parent_test:
 	kfree(km, w);
 }
 
-void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r)
+void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r, float alt_diff_frac)
 {
 	int32_t i, n_aux, n = *n_regs, has_cigar = 0, no_cigar = 0;
 	mm128_t *aux;
@@ -177,13 +196,11 @@ void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r)
 	t = (mm_reg1_t*)kmalloc(km, n * sizeof(mm_reg1_t));
 	for (i = n_aux = 0; i < n; ++i) {
 		if (r[i].inv || r[i].cnt > 0) { // squeeze out elements with cnt==0 (soft deleted)
-			if (r[i].p) {
+			int score;
-				aux[n_aux].x = (uint64_t)r[i].p->dp_max << 32 | r[i].hash;
+			if (r[i].p) score = r[i].p->dp_max, has_cigar = 1;
-				has_cigar = 1;
+			else score = r[i].score, no_cigar = 1;
-			} else {
+			if (r[i].is_alt) score = mm_alt_score(score, alt_diff_frac);
-				aux[n_aux].x = (uint64_t)r[i].score << 32 | r[i].hash;
+			aux[n_aux].x = (uint64_t)score << 32 | r[i].hash;
-				no_cigar = 1;
-			}
 			aux[n_aux++].y = i;
 		} else if (r[i].p) {
 			free(r[i].p);

index.c

@@ -1,4 +1,37 @@
+/* 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 <stdlib.h>
+#include<map>
+#include <vector>
+#include <fstream>
+using namespace std;
 #include <assert.h>
 #if defined(WIN32) || defined(_WIN32)
 #include <io.h> // for open(2)
@@ -53,6 +86,37 @@ mm_idx_t *mm_idx_init(int w, int k, int b, int flag)
 	return mi;
 }
 
+
+void mm_idx_destroy_mm_hash(mm_idx_t *mi)
+{
+	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)
+{
+	uint32_t i;
+	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;
@@ -97,6 +161,81 @@ const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n)
 	}
 }
 
+//Output minimap2's hash table entries 
+void mm_idx_dump_hash(const char* f_name, const mm_idx_t *mi)  
+{
+	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, "Storing hash to %s \n", f_name);
+	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();
+
+}
+
+
 void mm_idx_stat(const mm_idx_t *mi)
 {
 	int n = 0, n1 = 0;
@@ -316,6 +455,7 @@ static void *worker_pipeline(void *shared, int step, void *in)
 				} else seq->name = 0;
 				seq->len = s->seq[i].l_seq;
 				seq->offset = p->sum_len;
+				seq->is_alt = 0;
 				// copy the sequence
 				if (!(p->mi->flag & MM_I_NO_SEQ)) {
 					for (j = 0; j < seq->len; ++j) { // TODO: this is not the fastest way, but let's first see if speed matters here
@@ -414,6 +554,7 @@ mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const cha
 		}
 		p->offset = sum_len;
 		p->len = strlen(s);
+		p->is_alt = 0;
 		for (j = 0; j < p->len; ++j) {
 			int c = seq_nt4_table[(uint8_t)s[j]];
 			uint64_t o = sum_len + j;
@@ -500,6 +641,7 @@ mm_idx_t *mm_idx_load(FILE *fp)
 		}
 		fread(&s->len, 4, 1, fp);
 		s->offset = sum_len;
+		s->is_alt = 0;
 		sum_len += s->len;
 	}
 	for (i = 0; i < 1<<mi->b; ++i) {
@@ -607,6 +749,30 @@ int mm_idx_reader_eof(const mm_idx_reader_t *r) // TODO: in extremely rare cases
 #include "kseq.h"
 KSTREAM_DECLARE(gzFile, gzread)
 
+int mm_idx_alt_read(mm_idx_t *mi, const char *fn)
+{
+	int n_alt = 0;
+	gzFile fp;
+	kstream_t *ks;
+	kstring_t str = {0,0,0};
+	fp = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
+	if (fp == 0) return -1;
+	ks = ks_init(fp);
+	if (mi->h == 0) mm_idx_index_name(mi);
+	while (ks_getuntil(ks, KS_SEP_LINE, &str, 0) >= 0) {
+		char *p;
+		int id;
+		for (p = str.s; *p && !isspace(*p); ++p) { }
+		*p = 0;
+		id = mm_idx_name2id(mi, str.s);
+		if (id >= 0) mi->seq[id].is_alt = 1, ++n_alt;
+	}
+	mi->n_alt = n_alt;
+	if (mm_verbose >= 3)
+		fprintf(stderr, "[M::%s] found %d ALT contigs\n", __func__, n_alt);
+	return n_alt;
+}
+
 #define sort_key_bed(a) ((a).st)
 KRADIX_SORT_INIT(bed, mm_idx_intv1_t, sort_key_bed, 4)
 
@@ -627,7 +793,7 @@ mm_idx_intv_t *mm_idx_read_bed(const mm_idx_t *mi, const char *fn, int read_junc
 		char *p, *q, *bl, *bs;
 		int32_t i, id = -1, n_blk = 0;
 		for (p = q = str.s, i = 0;; ++p) {
-			if (*p == 0 || isspace(*p)) {
+			if (*p == 0 || *p == '\t') {
 				int32_t c = *p;
 				*p = 0;
 				if (i == 0) { // chr

kseq.h

@@ -89,7 +89,7 @@
 #ifndef KSTRING_T
 #define KSTRING_T kstring_t
 typedef struct __kstring_t {
-	unsigned l, m;
+	size_t l, m;
 	char *s;
 } kstring_t;
 #endif

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

@@ -4,15 +4,23 @@
 #include "ksw2.h"
 
 #ifdef __SSE2__
+#ifdef USE_SIMDE
+#include <simde/x86/sse2.h>
+#else
 #include <emmintrin.h>
+#endif
 
 #ifdef KSW_SSE2_ONLY
 #undef __SSE4_1__
 #endif
 
 #ifdef __SSE4_1__
+#ifdef USE_SIMDE
+#include <simde/x86/sse4.1.h>
+#else
 #include <smmintrin.h>
 #endif
+#endif
 
 #ifdef KSW_CPU_DISPATCH
 #ifdef __SSE4_1__

ksw2_exts2_sse.c

@@ -4,15 +4,22 @@
 #include "ksw2.h"
 
 #ifdef __SSE2__
+#ifdef USE_SIMDE
+#include <simde/x86/sse2.h>
+#else
 #include <emmintrin.h>
-
+#endif
 #ifdef KSW_SSE2_ONLY
 #undef __SSE4_1__
 #endif
 
 #ifdef __SSE4_1__
+#ifdef USE_SIMDE
+#include <simde/x86/sse4.1.h>
+#else
 #include <smmintrin.h>
 #endif
+#endif
 
 #ifdef KSW_CPU_DISPATCH
 #ifdef __SSE4_1__

ksw2_extz2_sse.c

@@ -3,15 +3,23 @@
 #include "ksw2.h"
 
 #ifdef __SSE2__
+#ifdef USE_SIMDE
+#include <simde/x86/sse2.h>
+#else
 #include <emmintrin.h>
+#endif
 
 #ifdef KSW_SSE2_ONLY
 #undef __SSE4_1__
 #endif
 
 #ifdef __SSE4_1__
+#ifdef USE_SIMDE
+#include <simde/x86/sse4.1.h>
+#else
 #include <smmintrin.h>
 #endif
+#endif
 
 #ifdef KSW_CPU_DISPATCH
 #ifdef __SSE4_1__

ksw2_ll_sse.c

@@ -1,9 +1,14 @@
 #include <stdlib.h>
 #include <stdint.h>
 #include <string.h>
-#include <emmintrin.h>
 #include "ksw2.h"
 
+#ifdef USE_SIMDE
+#include <simde/x86/sse2.h>
+#else
+#include <emmintrin.h>
+#endif
+
 #ifdef __GNUC__
 #define LIKELY(x) __builtin_expect((x),1)
 #define UNLIKELY(x) __builtin_expect((x),0)

main.c

@@ -1,13 +1,54 @@
+/* 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 <stdlib.h>
 #include <stdio.h>
 #include <string.h>
+#include <string>
 #include <errno.h>
 #include "bseq.h"
 #include "minimap.h"
 #include "mmpriv.h"
 #include "ketopt.h"
+#include <x86intrin.h>
 
-#define MM_VERSION "2.17-r963-dirty"
+#define MM_VERSION "2.18-r1015"
+using namespace std;
+
+#ifdef MANUAL_PROFILING
+uint64_t num_reads = 0, minimizer_hit_time = 0, dp_chaining_time = 0, alignment_time = 0;
+#endif
+
+#ifdef LISA_HASH
+#include "lisa_hash.h"
+lisa_hash<uint64_t, uint64_t> *lh;
+#endif
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -67,6 +108,10 @@ static ko_longopt_t long_options[] = {
 	{ "junc-bed",       ko_required_argument, 340 },
 	{ "junc-bonus",     ko_required_argument, 341 },
 	{ "sam-hit-only",   ko_no_argument,       342 },
+	{ "chain-gap-scale",ko_required_argument, 343 },
+	{ "alt",            ko_required_argument, 344 },
+	{ "alt-drop",       ko_required_argument, 345 },
+	{ "mask-len",       ko_required_argument, 346 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
 	{ "version",        ko_no_argument,       'V' },
@@ -104,12 +149,40 @@ static inline void yes_or_no(mm_mapopt_t *opt, int flag, int long_idx, const cha
 
 int main(int argc, char *argv[])
 {
-	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:yYPo:";
+#ifdef LISA_HASH
+	#if VECTORIZE && __AVX512BW__
+		fprintf(stderr, "Using LISA hash with AVX512-vectorized last-mile search.\n");
+	#else
+		fprintf(stderr, "Using LISA hash with sequential last-mile search.\n");
+	#endif	
+#else
+	fprintf(stderr, "Using default hash lookup.\n");
+#endif	
+
+#if defined(VECTORIZED_CHAINING) && defined(__AVX512BW__)
+	fprintf(stderr, "Using AVX512-vectorized chaining.\n");
+#else
+	fprintf(stderr, "Using default chaining.\n");
+#endif
+
+
+#if defined (ALIGN_AVX) && (defined(__AVX512BW__) || (defined(__AVX2__) && defined(APPLY_AVX2)))
+#ifdef __AVX512BW__
+	fprintf(stderr, "Using AVX512-vectorized alignment.\n");
+#elif __AVX2__
+	fprintf(stderr, "Using AVX2-vectorized alignment.\n");
+#endif
+#else
+	fprintf(stderr, "Using default SSE-vectorized alignment.\n");
+#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:Z:";
 	ketopt_t o = KETOPT_INIT;
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
 	int i, c, n_threads = 3, n_parts, old_best_n = -1;
-	char *fnw = 0, *rg = 0, *junc_bed = 0, *s;
+	uint64_t total_time = 0;
+	char *fnw = 0, *rg = 0, *junc_bed = 0, *s, *alt_list = 0;
 	FILE *fp_help = stderr;
 	mm_idx_reader_t *idx_rdr;
 	mm_idx_t *mi;
@@ -117,10 +190,13 @@ int main(int argc, char *argv[])
 	mm_verbose = 3;
 	liftrlimit();
 	mm_realtime0 = realtime();
+	double mapping_time = 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;
@@ -137,6 +213,7 @@ int main(int argc, char *argv[])
 
 	while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) {
 		if (c == 'w') ipt.w = atoi(o.arg);
+        	else if (c == 'Z') opt.L_hash = atoi(o.arg);
 		else if (c == 'k') ipt.k = atoi(o.arg);
 		else if (c == 'H') ipt.flag |= MM_I_HPC;
 		else if (c == 'd') fnw = o.arg; // the above are indexing related options, except -I
@@ -166,7 +243,7 @@ int main(int argc, char *argv[])
 		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);
-		else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(o.arg);
+		else if (c == 'K') opt.mini_batch_size = mm_parse_num(o.arg);
 		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;
@@ -211,6 +288,10 @@ int main(int argc, char *argv[])
 		else if (c == 340) junc_bed = o.arg; // --junc-bed
 		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 == 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 == 314) { // --frag
 			yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1);
 		} else if (c == 315) { // --secondary
@@ -337,7 +418,12 @@ 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);
+	total_time =  __rdtsc();
 	if (idx_rdr == 0) {
 		fprintf(stderr, "[ERROR] failed to open file '%s': %s\n", argv[o.ind], strerror(errno));
 		return 1;
@@ -350,6 +436,7 @@ int main(int argc, char *argv[])
 	if (opt.best_n == 0 && (opt.flag&MM_F_CIGAR) && mm_verbose >= 2)
 		fprintf(stderr, "[WARNING]\033[1;31m `-N 0' reduces alignment accuracy. Please use --secondary=no to suppress secondary alignments.\033[0m\n");
 	while ((mi = mm_idx_reader_read(idx_rdr, n_threads)) != 0) {
+		int ret;
 		if ((opt.flag & MM_F_CIGAR) && (mi->flag & MM_I_NO_SEQ)) {
 			fprintf(stderr, "[ERROR] the prebuilt index doesn't contain sequences.\n");
 			mm_idx_destroy(mi);
@@ -357,9 +444,11 @@ int main(int argc, char *argv[])
 			return 1;
 		}
 		if ((opt.flag & MM_F_OUT_SAM) && idx_rdr->n_parts == 1) {
-			int ret;
 			if (mm_idx_reader_eof(idx_rdr)) {
-				ret = mm_write_sam_hdr(mi, rg, MM_VERSION, argc, argv);
+				if (opt.split_prefix == 0)
+					ret = mm_write_sam_hdr(mi, rg, MM_VERSION, argc, argv);
+				else
+					ret = mm_write_sam_hdr(0, rg, MM_VERSION, argc, argv);
 			} else {
 				ret = mm_write_sam_hdr(0, rg, MM_VERSION, argc, argv);
 				if (opt.split_prefix == 0 && mm_verbose >= 2)
@@ -376,14 +465,42 @@ 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);
-		if (junc_bed) mm_idx_bed_read(mi, junc_bed, 1);
+		if(opt.L_hash == 1) {
-		if (!(opt.flag & MM_F_FRAG_MODE)) {
+			fprintf(stderr, "Generating lisa-hash..\n"); 
-			for (i = o.ind + 1; i < argc; ++i)
+			mm_idx_dump_hash(preset_arg.c_str(), mi); 
-				mm_map_file(mi, argv[i], &opt, n_threads);
+			fprintf(stderr, "Lisa-hash saving done.. \n");
-		} else {
+			exit(0); 
-			mm_map_file_frag(mi, argc - (o.ind + 1), (const char**)&argv[o.ind + 1], &opt, n_threads);
 		}
+		if (junc_bed) mm_idx_bed_read(mi, junc_bed, 1);
+		if (alt_list) mm_idx_alt_read(mi, alt_list);
+		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
+	mapping_time = 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);
+				if (ret < 0) break;
+			}
+		} else {
+			ret = mm_map_file_frag(mi, argc - (o.ind + 1), (const char**)&argv[o.ind + 1], &opt, n_threads);
+		}
+		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);
@@ -403,5 +520,14 @@ 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);
 	}
-	return 0;
+
+#ifdef MANUAL_PROFILING
+		fprintf(stderr, "\n Number of reads = %lld Minimizer hit time = %lld dp_chaining time = %lld alignment time = %lld total time = %lld \n", num_reads, minimizer_hit_time, dp_chaining_time, alignment_time, __rdtsc() - total_time);
+#endif
+	fprintf(stderr, "Total ticks: %lld \n",__rdtsc() - total_time);
+	fprintf(stderr, "\nMapping Real time: %.3f sec;\n", realtime() - mapping_time);
+#ifdef LISA_HASH
+	delete lh;
+#endif	
+return 0;
 }

map.c

@@ -1,3 +1,32 @@
+/* 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 <stdlib.h>
 #include <string.h>
 #include <assert.h>
@@ -9,6 +38,17 @@
 #include "mmpriv.h"
 #include "bseq.h"
 #include "khash.h"
+#include <x86intrin.h>
+
+#ifdef LISA_HASH
+#include "lisa_hash.h"
+extern lisa_hash<uint64_t, uint64_t> *lh;
+#endif
+
+#ifdef MANUAL_PROFILING
+extern uint64_t num_reads, minimizer_hit_time, dp_chaining_time, alignment_time;
+#endif
+
 
 struct mm_tbuf_s {
 	void *km;
@@ -87,6 +127,77 @@ typedef struct {
 	const uint64_t *cr;
 } mm_match_t;
 
+
+#ifdef LISA_HASH
+static mm_match_t *collect_matches_lisa_hash(void *km, int *_n_m, int max_occ, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos)
+{
+	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));
+	
+	int rep_st = 0, rep_en = 0, n_m;
+	size_t i;
+	mm_match_t *m;
+	*n_mini_pos = 0;
+	*mini_pos = (uint64_t*)kmalloc(km, mv->n * sizeof(uint64_t));
+	m = (mm_match_t*)kmalloc(km, mv->n * sizeof(mm_match_t));
+	
+	for (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); 
+	for (i = 0, n_m = 0, *rep_len = 0, *n_a = 0; i < mv->n; ++i) {
+		const uint64_t *cr;
+		mm128_t *p = &mv->a[i];
+		uint32_t q_pos = (uint32_t)p->y, q_span = p->x & 0xff;
+		int t;
+		cr = cr_batch[i]; t = t_batch[i];
+/*Correctness check for lisa_hash*/
+#ifdef LISA_HASH_ASSERT	
+		int t_minimap2_original;
+		const uint64_t *cr_minimap2_hash = mm_idx_get(mi, p->x>>8, &t);
+		cr_minimap2_hash = mm_idx_get(mi, p->x>>8, &t_minimap2_original);
+		assert(t == t_minimap2_original);	
+#endif
+		if (t >= max_occ) {
+
+			int en = (q_pos >> 1) + 1, st = en - q_span;
+			if (st > rep_en) {
+				*rep_len += rep_en - rep_st;
+				rep_st = st, rep_en = en;
+			} else rep_en = en;
+		} else {
+#ifdef LISA_HASH_ASSERT
+			//Correctness assertion
+			for(int itr = 0; itr < t; itr++){
+				assert((cr[itr] == cr_minimap2_hash[itr]));
+			}
+#endif	
+			mm_match_t *q = &m[n_m++];
+			q->q_pos = q_pos, q->q_span = q_span, q->cr = cr, q->n = t, q->seg_id = p->y >> 32;
+			q->is_tandem = 0;
+			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;
+			*n_a += q->n;
+			(*mini_pos)[(*n_mini_pos)++] = (uint64_t)q_span<<32 | q_pos>>1;
+		}
+	}
+
+	free(cr_batch);
+	free(t_batch);
+	free(minimizers);
+	free(lisa_pos);
+
+	*rep_len += rep_en - rep_st;
+	*_n_m = n_m;
+	return m;
+}
+#endif
+
+
 static mm_match_t *collect_matches(void *km, int *_n_m, int max_occ, const mm_idx_t *mi, const mm128_v *mv, int64_t *n_a, int *rep_len, int *n_mini_pos, uint64_t **mini_pos)
 {
 	int rep_st = 0, rep_en = 0, n_m;
@@ -146,6 +257,49 @@ static inline int skip_seed(int flag, uint64_t r, const mm_match_t *q, const cha
 	return 0;
 }
 
+
+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)
+{
+	int i, n_m;
+	mm_match_t *m;
+	mm128_t *a;
+#ifndef LISA_HASH
+	m = collect_matches(km, &n_m, max_occ, mi, mv, n_a, rep_len, n_mini_pos, mini_pos);
+#else
+	m = collect_matches_lisa_hash(km, &n_m, max_occ, mi, mv, n_a, rep_len, n_mini_pos, mini_pos);
+#endif
+
+	a = (mm128_t*)kmalloc(km, *n_a * sizeof(mm128_t));
+	for (i = 0, *n_a = 0; i < n_m; ++i) {
+		
+		mm_match_t *q = &m[i];
+		const uint64_t *r = q->cr;
+		uint32_t k;
+		for (k = 0; k < q->n; ++k) {
+			uint64_t r_k = r[k];
+			int32_t is_self, rpos = (uint32_t)r_k >> 1;
+			mm128_t *p;
+			if (skip_seed(opt->flag, r_k, q, qname, qlen, mi, &is_self)) continue;
+			p = &a[(*n_a)++];
+			if ((r_k&1) == (q->q_pos&1)) { // forward strand
+				p->x = (r_k & 0xffffffff00000000ULL) | rpos;
+				p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
+			} else { // reverse strand
+				p->x = 1ULL<<63 | (r_k & 0xffffffff00000000ULL) | rpos;
+				p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
+			}
+			p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
+			if (q->is_tandem) p->y |= MM_SEED_TANDEM;
+			if (is_self) p->y |= MM_SEED_SELF;
+		}
+	}
+	
+	kfree(km, m);
+	radix_sort_128x(a, a + (*n_a));
+	return a;
+}
+
 static mm128_t *collect_seed_hits_heap(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)
 {
@@ -212,44 +366,10 @@ static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max
 	return a;
 }
 
-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)
-{
-	int i, n_m;
-	mm_match_t *m;
-	mm128_t *a;
-	m = collect_matches(km, &n_m, max_occ, mi, mv, n_a, rep_len, n_mini_pos, mini_pos);
-	a = (mm128_t*)kmalloc(km, *n_a * sizeof(mm128_t));
-	for (i = 0, *n_a = 0; i < n_m; ++i) {
-		mm_match_t *q = &m[i];
-		const uint64_t *r = q->cr;
-		uint32_t k;
-		for (k = 0; k < q->n; ++k) {
-			int32_t is_self, rpos = (uint32_t)r[k] >> 1;
-			mm128_t *p;
-			if (skip_seed(opt->flag, r[k], q, qname, qlen, mi, &is_self)) continue;
-			p = &a[(*n_a)++];
-			if ((r[k]&1) == (q->q_pos&1)) { // forward strand
-				p->x = (r[k]&0xffffffff00000000ULL) | rpos;
-				p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
-			} else { // reverse strand
-				p->x = 1ULL<<63 | (r[k]&0xffffffff00000000ULL) | rpos;
-				p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos>>1) + 1 - q->q_span) - 1);
-			}
-			p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
-			if (q->is_tandem) p->y |= MM_SEED_TANDEM;
-			if (is_self) p->y |= MM_SEED_SELF;
-		}
-	}
-	kfree(km, m);
-	radix_sort_128x(a, a + (*n_a));
-	return a;
-}
-
 static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_idx_t *mi, void *km, int qlen, int n_segs, const int *qlens, int *n_regs, mm_reg1_t *regs, mm128_t *a)
 {
 	if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
-		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL);
+		mm_set_parent(km, opt->mask_level, opt->mask_len, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
 		if (n_segs <= 1) mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 		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);
 		if (!(opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN))) // long join not working well without primary chains
@@ -262,7 +382,7 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k
 	if (!(opt->flag & MM_F_CIGAR)) return regs;
 	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, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL);
+		mm_set_parent(km, opt->mask_level, opt->mask_len, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
 		mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
 		mm_set_sam_pri(*n_regs, regs);
 	}
@@ -271,6 +391,11 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k
 
 void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
 {
+	
+#ifdef MANUAL_PROFILING
+	num_reads++;
+#endif
+
 	int i, j, rep_len, qlen_sum, n_regs0, n_mini_pos;
 	int max_chain_gap_qry, max_chain_gap_ref, is_splice = !!(opt->flag & MM_F_SPLICE), is_sr = !!(opt->flag & MM_F_SR);
 	uint32_t hash;
@@ -293,8 +418,18 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 
 	collect_minimizers(b->km, opt, mi, n_segs, qlens, seqs, &mv);
 	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);
+	else {
 
+#ifdef MANUAL_PROFILING
+		uint64_t mm_hit_start = __rdtsc();
+#endif
+
+		a = collect_seed_hits(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
+
+#ifdef MANUAL_PROFILING
+		minimizer_hit_time += (__rdtsc() - mm_hit_start);
+#endif
+	}
 	if (mm_dbg_flag & MM_DBG_PRINT_SEED) {
 		fprintf(stderr, "RS\t%d\n", rep_len);
 		for (i = 0; i < n_a; ++i)
@@ -312,8 +447,15 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		max_chain_gap_ref = opt->max_frag_len - qlen_sum;
 		if (max_chain_gap_ref < opt->max_gap) max_chain_gap_ref = opt->max_gap;
 	} else max_chain_gap_ref = opt->max_gap;
+#ifdef MANUAL_PROFILING
+	uint64_t dp_start = __rdtsc();
+#endif
+	a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score, opt->chain_gap_scale, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 
-	a = mm_chain_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, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
+
+#ifdef MANUAL_PROFILING
+	dp_chaining_time += (__rdtsc() - dp_start);
+#endif
 
 	if (opt->max_occ > opt->mid_occ && rep_len > 0) {
 		int rechain = 0;
@@ -335,13 +477,18 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 			kfree(b->km, mini_pos);
 			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 = mm_chain_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, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
+
+			a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score, opt->chain_gap_scale, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 		}
 	}
 	b->frag_gap = max_chain_gap_ref;
 	b->rep_len = rep_len;
 
 	regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a);
+	if (mi->n_alt) {
+		mm_mark_alt(mi, n_regs0, regs0);
+		mm_hit_sort(b->km, &n_regs0, regs0, opt->alt_drop); // this step can be merged into mm_gen_regs(); will do if this shows up in profile
+	}
 
 	if (mm_dbg_flag & MM_DBG_PRINT_SEED)
 		for (j = 0; j < n_regs0; ++j)
@@ -361,7 +508,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		seg = mm_seg_gen(b->km, hash, n_segs, qlens, n_regs0, regs0, n_regs, regs, a); // split fragment chain to separate segment chains
 		free(regs0);
 		for (i = 0; i < n_segs; ++i) {
-			mm_set_parent(b->km, opt->mask_level, n_regs[i], regs[i], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL); // update mm_reg1_t::parent
+			mm_set_parent(b->km, opt->mask_level, opt->mask_len, n_regs[i], regs[i], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop); // update mm_reg1_t::parent
 			regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a);
 			mm_set_mapq(b->km, n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr);
 		}
@@ -399,7 +546,8 @@ mm_reg1_t *mm_map(const mm_idx_t *mi, int qlen, const char *seq, int *n_regs, mm
  **************************/
 
 typedef struct {
-	int mini_batch_size, n_processed, n_threads, n_fp;
+	int n_processed, n_threads, n_fp;
+	int64_t mini_batch_size;
 	const mm_mapopt_t *opt;
 	mm_bseq_file_t **fp;
 	const mm_idx_t *mi;
@@ -425,6 +573,7 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
 	int qlens[MM_MAX_SEG], j, off = s->seg_off[i], pe_ori = s->p->opt->pe_ori;
 	const char *qseqs[MM_MAX_SEG];
 	mm_tbuf_t *b = s->buf[tid];
+
 	assert(s->n_seg[i] <= MM_MAX_SEG);
 	if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
 		fprintf(stderr, "QR\t%s\t%d\t%d\n", s->seq[off].name, tid, s->seq[off].l_seq);
@@ -503,8 +652,8 @@ static void merge_hits(step_t *s)
 					}
 				}
 			}
-			mm_hit_sort(km, &s->n_reg[k], s->reg[k]);
+			mm_hit_sort(km, &s->n_reg[k], s->reg[k], opt->alt_drop);
-			mm_set_parent(km, opt->mask_level, s->n_reg[k], s->reg[k], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL);
+			mm_set_parent(km, opt->mask_level, opt->mask_len, s->n_reg[k], s->reg[k], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL, opt->alt_drop);
 			if (!(opt->flag & MM_F_ALL_CHAINS)) {
 				mm_select_sub(km, opt->pri_ratio, s->p->mi->k*2, opt->best_n, &s->n_reg[k], s->reg[k]);
 				mm_set_sam_pri(s->n_reg[k], s->reg[k]);
@@ -556,6 +705,7 @@ static void *worker_pipeline(void *shared, int step, void *in)
 		else kt_for(p->n_threads, worker_for, in, ((step_t*)in)->n_frag);
 		return in;
     } else if (step == 2) { // step 2: output
+
 		void *km = 0;
         step_t *s = (step_t*)in;
 		const mm_idx_t *mi = p->mi;
@@ -564,6 +714,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 +749,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

@@ -1,3 +1,32 @@
+/* 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>
+*/
 #ifndef MINIMAP2_H
 #define MINIMAP2_H
 
@@ -58,12 +87,14 @@ typedef struct {
 	char *name;      // name of the db sequence
 	uint64_t offset; // offset in mm_idx_t::S
 	uint32_t len;    // length
+	uint32_t is_alt;
 } mm_idx_seq_t;
 
 typedef struct {
 	int32_t b, w, k, flag;
 	uint32_t n_seq;            // number of reference sequences
 	int32_t index;
+	int32_t n_alt;
 	mm_idx_seq_t *seq;         // sequence name, length and offset
 	uint32_t *S;               // 4-bit packed sequence
 	struct mm_idx_bucket_s *B; // index (hidden)
@@ -91,7 +122,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, dummy:7;
+	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;
@@ -100,7 +131,7 @@ typedef struct {
 // indexing and mapping options
 typedef struct {
 	short k, w, flag, bucket_bits;
-	int mini_batch_size;
+	int64_t mini_batch_size;
 	uint64_t batch_size;
 } mm_idxopt_t;
 
@@ -117,8 +148,10 @@ typedef struct {
 	int max_chain_skip, max_chain_iter;
 	int min_cnt;         // min number of minimizers on each chain
 	int min_chain_score; // min chaining score
+	float chain_gap_scale;
 
 	float mask_level;
+	int mask_len;
 	float pri_ratio;
 	int best_n;      // top best_n chains are subjected to DP alignment
 
@@ -126,6 +159,8 @@ typedef struct {
 	int min_join_flank_sc;
 	float min_join_flank_ratio;
 
+	float alt_drop;
+
 	int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties
 	int sc_ambi; // score when one or both bases are "N"
 	int noncan;      // cost of non-canonical splicing sites
@@ -143,10 +178,12 @@ typedef struct {
 	int32_t min_mid_occ;
 	int32_t mid_occ;     // ignore seeds with occurrences above this threshold
 	int32_t max_occ;
-	int mini_batch_size; // size of a batch of query bases to process in parallel
+	int64_t mini_batch_size; // size of a batch of query bases to process in parallel
 	int64_t max_sw_mat;
 
 	const char *split_prefix;
+    	// Store minimizer hash to a file as key and list of values
+    	int L_hash;
 } mm_mapopt_t;
 
 // index reader
@@ -261,6 +298,14 @@ 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
  *
@@ -290,6 +335,21 @@ void mm_idx_stat(const mm_idx_t *idx);
  */
 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
  *
@@ -368,6 +428,7 @@ int mm_idx_index_name(mm_idx_t *mi);
 int mm_idx_name2id(const mm_idx_t *mi, const char *name);
 int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
 
+int mm_idx_alt_read(mm_idx_t *mi, const char *fn);
 int mm_idx_bed_read(mm_idx_t *mi, const char *fn, int read_junc);
 int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uint8_t *s);
 

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "4 May 2019" "minimap2-2.17 (r941)" "Bioinformatics tools"
+.TH minimap2 1 "9 April 2021" "minimap2-2.18 (r1015)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -121,6 +121,14 @@ provided as the target sequences, options
 .BR -w ,
 .B -I
 will be effectively overridden by the options stored in the index file.
+.TP
+.BI --alt \ FILE
+List of ALT contigs [null]
+.TP
+.BI --alt-drop \ FLOAT
+Drop ALT hits by
+.I FLOAT
+fraction when ranking and computing mapping quality [0.15]
 .SS Mapping options
 .TP 10
 .BI -f \ FLOAT | INT1 [, INT2 ]
@@ -229,7 +237,14 @@ or more of the shorter chain [0.5]
 .B --hard-mask-level
 Honor option
 .B -M
-and disable a heurstic to save unmapped subsequences.
+and disable a heurstic to save unmapped subsequences and disables
+.BR --mask-len .
+.TP
+.BI --mask-len \ NUM
+Keep an alignment if dropping it leaves an unaligned region on query longer than
+.IR INT
+[inf]. Effective without
+.BR --hard-mask-level .
 .TP
 .BI --max-chain-skip \ INT
 A heuristics that stops chaining early [25]. Minimap2 uses dynamic programming
@@ -245,6 +260,9 @@ Check up to
 partial chains during chaining [5000]. This is a heuristic to avoid quadratic
 time complexity in the worst case.
 .TP
+.BI --chain-gap-scale \ FLOAT
+Scale of gap cost during chaining [1.0]
+.TP
 .B --no-long-join
 Disable the long gap patching heuristic. When this option is applied, the
 maximum alignment gap is mostly controlled by
@@ -373,7 +391,7 @@ BED12 file can be converted from GTF/GFF3 with `paftools.js gff2bed anno.gtf'
 .BR --junc-bonus \ INT
 Score bonus for a splice donor or acceptor found in annotation (effective with
 .BR --junc-bed )
-[0].
+[9].
 .TP
 .BI --end-seed-pen \ INT
 Drop a terminal anchor if

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = '2.17-r949-dirty';
+var paftools_version = '2.18-r1015';
 
 /*****************************
  ***** Library functions *****
@@ -640,6 +640,23 @@ function paf_asmstat(args)
 		}
 	}
 
+	function AUN(lens, tot) {
+		lens.sort(function(a,b) { return b - a; });
+		if (tot == null) {
+			tot = 0;
+			for (var k = 0; k < lens.length; ++k)
+				tot += lens[k];
+		}
+		var x = 0, y = 0;
+		for (var k = 0; k < lens.length; ++k) {
+			var l = x + lens[k] <= tot? lens[k] : tot - x;
+			x += lens[k];
+			y += l * (l / tot);
+			if (x >= tot) break;
+		}
+		return y.toFixed(0);
+	}
+
 	function count_bp(bp, min_blen, min_gap) {
 		var n_bp = 0;
 		for (var k = 0; k < bp.length; ++k)
@@ -660,7 +677,7 @@ function paf_asmstat(args)
 		return (NM - n_gaps + n_gapo) / (n_M + n_gapo);
 	}
 
-	var labels = ['Length', 'l_cov', 'Rcov', 'Rdup', 'Qcov', 'NG75', 'NG50', 'NGA50', '#breaks', 'bp(' + min_seg_len + ',0)', 'bp(' + min_seg_len + ',10k)'];
+	var labels = ['Length', 'l_cov', 'Rcov', 'Rdup', 'Qcov', 'NG75', 'NG50', 'NGA50', 'AUNGA', '#breaks', 'bp(' + min_seg_len + ',0)', 'bp(' + min_seg_len + ',10k)'];
 	var rst = [];
 	for (var i = 0; i < labels.length; ++i)
 		rst[i] = [];
@@ -688,11 +705,9 @@ function paf_asmstat(args)
 			qinfo[t[0]].bp = [];
 			if (t.length < 9 || t[5] == "*") continue;
 			if (!/\ttp:A:[PI]/.test(line)) continue;
-			if ((m = /\tcg:Z:(\S+)/.exec(line)) == null) continue;
+			var cigar = (m = /\tcg:Z:(\S+)/.exec(line)) != null? m[1] : null;
-			var cigar = m[1];
+			var NM = (m = /\tNM:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : null;
-			if ((m = /\tNM:i:(\d+)/.exec(line)) == null) continue;
+			var diff = cigar != null && NM != null? compute_diff(cigar, NM) : 0;
-			var NM = parseInt(m[1]);
-			var diff = compute_diff(cigar, NM);
 			t[2] = parseInt(t[2]);
 			t[3] = parseInt(t[3]);
 			t[7] = parseInt(t[7]);
@@ -765,10 +780,13 @@ function paf_asmstat(args)
 		// compute NGA50
 		rst[7][i] = N50(qblock_len, ref_len, 0.5);
 
+		// compute AUNGA
+		rst[8][i] = AUN(qblock_len, ref_len);
+
 		// compute break points
-		rst[8][i] = n_breaks;
+		rst[9][i] = n_breaks;
-		rst[9][i] = count_bp(bp, 500, 0);
+		rst[10][i] = count_bp(bp, 500, 0);
-		rst[10][i] = count_bp(bp, 500, 10000);
+		rst[11][i] = count_bp(bp, 500, 10000);
 
 		// nb-plot; NOT USED
 		/*
@@ -887,23 +905,24 @@ function paf_asmgene(args)
 	gene_nr[gene_list[last][0]] = 1;
 
 	// count and print
-	var col1 = ["full_sgl", "full_dup", "frag", "part50+", "part10+", "part10-"];
+	var col1 = ["full_sgl", "full_dup", "frag", "part50+", "part10+", "part10-", "dup_cnt", "dup_sum"];
 	var rst = [];
 	for (var k = 0; k < col1.length; ++k) {
 		rst[k] = [];
 		for (var i = 0; i < n_fn; ++i)
 			rst[k][i] = 0;
 	}
-	for (var g in gene) {
+	for (var g in gene) { // count single-copy genes
 		if (gene[g][0] == null || gene[g][0][0] != 1) continue;
 		if (gene_nr[g] == null) continue;
 		if (auto_only && /^(chr)?[XY]$/.test(refpos[g][2])) continue;
 		for (var i = 0; i < n_fn; ++i) {
 			if (gene[g][i] == null) {
-				rst[4][i]++;
+				rst[5][i]++;
 				if (print_err) print('M', header[i], refpos[g].join("\t"));
-			} else if (gene[g][i][0] == 1) rst[0][i]++;
+			} else if (gene[g][i][0] == 1) {
-			else if (gene[g][i][0] > 1) {
+				rst[0][i]++;
+			} else if (gene[g][i][0] > 1) {
 				rst[1][i]++;
 				if (print_err) print('D', header[i], refpos[g].join("\t"));
 			} else if (gene[g][i][1] >= opt.min_cov) {
@@ -921,6 +940,19 @@ function paf_asmgene(args)
 			}
 		}
 	}
+	for (var g in gene) { // count multi-copy genes
+		if (gene[g][0] == null || gene[g][0][0] <= 1) continue;
+		if (gene_nr[g] == null) continue;
+		if (auto_only && /^(chr)?[XY]$/.test(refpos[g][2])) continue;
+		for (var i = 0; i < n_fn; ++i) {
+			if (gene[g][i] != null) rst[7][i] += gene[g][i][0];
+			if (gene[g][i] != null && gene[g][i][0] > 1) {
+				rst[6][i]++;
+			} else if (print_err) {
+				print('d', header[i], gene[g][0][0], refpos[g].join("\t"));
+			}
+		}
+	}
 	print('H', 'Metric', header.join("\t"));
 	for (var k = 0; k < rst.length; ++k) {
 		print('X', col1[k], rst[k].join("\t"));
@@ -930,12 +962,13 @@ function paf_asmgene(args)
 
 function paf_stat(args)
 {
-	var c, gap_out_len = null;
+	var c, gap_out_len = null, count_err = false;
-	while ((c = getopt(args, "l:")) != null)
+	while ((c = getopt(args, "cl:")) != null)
 		if (c == 'l') gap_out_len = parseInt(getopt.arg);
+		else if (c == 'c') count_err = true;
 
 	if (getopt.ind == args.length) {
-		print("Usage: paftools.js stat [-l gapOutLen] <in.sam>|<in.paf>");
+		print("Usage: paftools.js stat [-c] [-l gapOutLen] <in.sam>|<in.paf>");
 		exit(1);
 	}
 
@@ -966,7 +999,7 @@ function paf_stat(args)
 		if (line.charAt(0) != '@') {
 			var t = line.split("\t", 12);
 			var m, rs, cigar = null, is_pri = false, is_sam = false, is_rev = false, tname = null;
-			var atlen = null, aqlen, qs, qe, mapq, ori_qlen;
+			var atlen = null, aqlen, qs, qe, mapq, ori_qlen, NM = null;
 			if (t.length < 2) continue;
 			if (t[4] == '+' || t[4] == '-' || t[4] == '*') { // PAF
 				if (t[4] == '*') continue; // unmapped
@@ -974,6 +1007,8 @@ function paf_stat(args)
 					++n_2nd;
 					continue;
 				}
+				if ((m = /\tNM:i:(\d+)/.exec(line)) != null)
+					NM = parseInt(m[1]);
 				if ((m = /\tcg:Z:(\S+)/.exec(line)) != null)
 					cigar = m[1];
 				if (cigar == null) {
@@ -995,6 +1030,8 @@ function paf_stat(args)
 					++n_2nd;
 					continue;
 				}
+				if ((m = /\tNM:i:(\d+)/.exec(line)) != null)
+					NM = parseInt(m[1]);
 				cigar = t[5];
 				tname = t[2];
 				rs = parseInt(t[3]) - 1;
@@ -1013,11 +1050,13 @@ function paf_stat(args)
 				++n_seq, last = t[0];
 			}
 			var M = 0, tl = 0, ql = 0, clip = [0, 0], n_cigar = 0, sclip = 0;
+			var n_gapo = 0, n_gap_all = 0, l_match = 0;
 			while ((m = re.exec(cigar)) != null) {
 				var l = parseInt(m[1]);
 				++n_cigar;
 				if (m[2] == 'M' || m[2] == '=' || m[2] == 'X') {
 					tl += l, ql += l, M += l;
+					l_match += l;
 				} else if (m[2] == 'I' || m[2] == 'D') {
 					var type;
 					if (l < 50) type = 0;
@@ -1030,6 +1069,7 @@ function paf_stat(args)
 					else tl += l, ++n_gap[1][type];
 					if (gap_out_len != null && l >= gap_out_len)
 						print(t[0], ql, is_rev? '-' : '+', tname, rs + tl, m[2], l);
+					++n_gapo, n_gap_all += l;
 				} else if (m[2] == 'N') {
 					tl += l;
 				} else if (m[2] == 'S') {
@@ -1047,6 +1087,12 @@ function paf_stat(args)
 				qs = clip[is_rev? 1 : 0], qe = qs + ql;
 				ori_qlen = clip[0] + ql + clip[1];
 			}
+			if (count_err && NM != null) {
+				var n_mm = NM - n_gap_all;
+				if (n_mm < 0) warn("WARNING: NM is smaller than the number of gaps at line " + lineno);
+				if (n_mm < 0) n_mm = 0;
+				print(t[0], ori_qlen, t[11], ori_qlen - (qe - qs), NM, l_match + n_gap_all, n_mm + n_gapo, l_match + n_gapo);
+			}
 			regs.push([qs, qe]);
 			last_qlen = ori_qlen;
 		}
@@ -1059,7 +1105,7 @@ function paf_stat(args)
 	file.close();
 	buf.destroy();
 
-	if (gap_out_len == null) {
+	if (gap_out_len == null && !count_err) {
 		print("Number of mapped sequences: " + n_seq);
 		print("Number of primary alignments: " + n_pri);
 		print("Number of secondary alignments: " + n_2nd);
@@ -2480,6 +2526,310 @@ function paf_ov_eval(args)
 	print((100 * (1 - n_missing / n_ovlp)).toFixed(2) + "% sensitivity");
 }
 
+function paf_vcfstat(args)
+{
+	var c, ts = { "AG":1, "GA":1, "CT":1, "TC":1 };
+	while ((c = getopt(args, "")) != null) {
+	}
+	var buf = new Bytes();
+	var file = args.length == getopt.ind? new File() : new File(args[getopt.ind]);
+	var x = { sub:0, ts:0, tv:0, ins:0, del:0, ins1:0, del1:0, ins2:0, del2:0, ins50:0, del50:0, ins1k:0, del1k:0, ins7k:0, del7k:0, insinf:0, delinf:0 };
+	while (file.readline(buf) >= 0) {
+		var t = buf.toString().split("\t");
+		if (t[0][0] == '#') continue;
+		var alt = t[4].split(",");
+		var ref = t[3];
+		for (var i = 0; i < alt.length; ++i) {
+			var a = alt[i];
+			if (a[0] == '<' || a[1] == '>') continue;
+			var l = ref.length < a.length? ref.length : a.length;
+			for (var j = 0; j < l; ++j) {
+				if (ref[j] != a[j]) {
+					++x.sub;
+					if (ts[ref[j] + a[j]]) ++x.ts;
+					else ++x.tv;
+				}
+			}
+			var d = a.length - ref.length;
+			if (d > 0) {
+				++x.ins;
+				if (d == 1) ++x.ins1;
+				else if (d == 2) ++x.ins2;
+				else if (d < 50) ++x.ins50;
+				else if (d < 1000) ++x.ins1k;
+				else if (d < 7000) ++x.ins7k;
+				else ++x.insinf;
+			} else if (d < 0) {
+				d = -d;
+				++x.del;
+				if (d == 1) ++x.del1;
+				else if (d == 2) ++x.del2;
+				else if (d < 50) ++x.del50;
+				else if (d < 1000) ++x.del1k;
+				else if (d < 7000) ++x.del7k;
+				else ++x.delinf;
+			}
+		}
+	}
+	file.close();
+	buf.destroy();
+	print("# substitutions: " + x.sub);
+	print("ts/tv: " + (x.ts / x.tv).toFixed(3));
+	print("# insertions: " + x.ins);
+	print("# 1bp insertions: " + x.ins1);
+	print("# 2bp insertions: " + x.ins2);
+	print("# [3,50) insertions: " + x.ins50);
+	print("# [50,1000) insertions: " + x.ins1k);
+	print("# [1000,7000) insertions: " + x.ins7k);
+	print("# >=7000 insertions: " + x.insinf);
+	print("# deletions: " + x.del);
+	print("# 1bp deletions: " + x.del1);
+	print("# 2bp deletions: " + x.del2);
+	print("# [3,50) deletions: " + x.del50);
+	print("# [50,1000) deletions: " + x.del1k);
+	print("# [1000,7000) deletions: " + x.del7k);
+	print("# >=7000 deletions: " + x.delinf);
+}
+
+function paf_parseNum(s) {
+	var m, x = null;
+	if ((m = /^(\d*\.?\d*)([mMgGkK]?)/.exec(s)) != null) {
+		x = parseFloat(m[1]);
+		if (m[2] == 'k' || m[2] == 'K') x *= 1000;
+		else if (m[2] == 'm' || m[2] == 'M') x *= 1000000;
+		else if (m[2] == 'g' || m[2] == 'G') x *= 1000000000;
+	}
+	return Math.floor(x + .499);
+}
+
+function paf_misjoin(args)
+{
+	var c, min_seg_len = 1000000, max_gap = 1000000, fn_cen = null, show_long = false, show_err = false, cen_ratio = 0.5;
+	var n_diff = [0, 0], n_gap = [0, 0], n_inv = [0, 0], n_inv_end = [0, 0];
+	while ((c = getopt(args, "l:g:c:per:")) != null) {
+		if (c == 'l') min_seg_len = paf_parseNum(getopt.arg);
+		else if (c == 'g') max_gap = paf_parseNum(getopt.arg);
+		else if (c == 'c') fn_cen = getopt.arg;
+		else if (c == 'r') cen_ratio = parseFloat(getopt.arg);
+		else if (c == 'p') show_long = true;
+		else if (c == 'e') show_err = true;
+	}
+	if (args.length == getopt.ind) {
+		print("Usage: paftools.js misjoin [options] <in.paf>");
+		print("Options:");
+		print("  -c FILE   BED for centromeres []");
+		print("  -r FLOAT  count a centromeric event if overlap ratio > FLOAT [" + cen_ratio + "]");
+		print("  -l NUM    min alignment block length [1m]");
+		print("  -g NUM    max gap size [1m]");
+		print("  -e        output misjoins not involving centromeres");
+		print("  -p        output long alignment blocks for debugging");
+		return;
+	}
+	var cen = {};
+	var file, buf = new Bytes();
+	if (fn_cen != null) {
+		file = new File(fn_cen);
+		while (file.readline(buf) >= 0) {
+			var t = buf.toString().split("\t");
+			if (cen[t[0]] == null) cen[t[0]] = [];
+			cen[t[0]].push([parseInt(t[1]), parseInt(t[2])]);
+		}
+		file.close();
+	}
+
+	function test_cen(cen, chr, st, en) {
+		var b = cen[chr], len = 0;
+		if (b == null) return false;
+		for (var j = 0; j < b.length; ++j)
+			if (b[j][0] < en && b[j][1] > st) {
+				var s = b[j][0] > st? b[j][0] : st;
+				var e = b[j][1] < en? b[j][1] : en;
+				len += e - s;
+			}
+		return len < (en - st) * cen_ratio? false : true;
+	}
+
+	function process(a) {
+		var k = 0;
+		for (var i = 0; i < a.length; ++i) {
+			for (var j = 1; j <= 3;  ++j) a[i][j] = parseInt(a[i][j]);
+			for (var j = 6; j <= 11; ++j) a[i][j] = parseInt(a[i][j]);
+			if (a[i][10] >= min_seg_len) a[k++] = a[i];
+		}
+		a.length = k;
+		if (a.length == 1) return;
+		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];
+			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]);
+			if (a[i-1][5] != a[i][5]) { // different chr
+				if (ov[0] || ov[1]) ++n_diff[1];
+				else if (show_err) {
+					print("J", a[i-1].slice(0, 12).join("\t"));
+					print("J", a[i].slice(0, 12).join("\t"));
+				}
+				++n_diff[0];
+			} else if (a[i-1][4] == a[i][4]) { // a gap
+				var dq = a[i][2] - a[i-1][3];
+				var dr = a[i][4] == '+'? a[i][7] - a[i-1][8] : a[i-1][7] - a[i][8];
+				var gap = dr > dq? dr - dq : dq - dr;
+				if (gap > max_gap) {
+					if (ov[0] || ov[1]) ++n_gap[1];
+					else if (show_err) {
+						print("G", a[i-1].slice(0, 12).join("\t"));
+						print("G", a[i].slice(0, 12).join("\t"));
+					}
+					++n_gap[0];
+				}
+			} else if (i + 1 < a.length && a[i+1][4] == a[i-1][4]) { // bracketed inversion
+				if (ov[0] || ov[1]) ++n_inv[1];
+				else if (show_err) {
+					print("M", a[i-1].slice(0, 12).join("\t"));
+					print("M", a[i].slice(0, 12).join("\t"));
+					print("M", a[i+1].slice(0, 12).join("\t"));
+				}
+				++n_inv[0];
+				++i;
+			} else { // hanging inversion
+				if (ov[0] || ov[1]) ++n_inv_end[1];
+				++n_inv_end[0];
+			}
+		}
+	}
+
+	file = args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]);
+	var a = [];
+	while (file.readline(buf) >= 0) {
+		var t = buf.toString().split("\t");
+		if (a.length > 0 && a[0][0] != t[0]) {
+			process(a);
+			a.length = 0;
+		}
+		a.push(t);
+	}
+	if (a.length > 0) process(a);
+	file.close();
+	buf.destroy();
+	print("# inter-chromosomal misjoins: " + n_diff.join(","));
+	print("# intra-chromosomal gaps: " + n_gap.join(","));
+	print("# candidate inversions in the middle: " + n_inv.join(","));
+	print("# candidate inversions at contig ends: " + n_inv_end.join(","));
+}
+
+function paf_sveval(args)
+{
+	var c, min_flt = 30, min_size = 50, max_size = 10000, win_size = 500, print_err = false, bed_fn = null;
+	while ((c = getopt(args, "f:i:x:w:er:")) != null) {
+		if (c == 'f') min_flt = paf_parseNum(getopt.arg);
+		else if (c == 'i') min_size = paf_parseNum(getopt.arg);
+		else if (c == 'x') max_size = paf_parseNum(getopt.arg);
+		else if (c == 'w') win_size = paf_parseNum(getopt.arg);
+		else if (c == 'r') bed_fn = getopt.arg;
+		else if (c == 'e') print_err = true;
+	}
+	if (args.length - getopt.ind < 2) {
+		print("Usage: paftools.js sveval [options] <base.vcf> <call.vcf>");
+		print("Options:");
+		print("  -r FILE    confident region in BED []");
+		print("  -f INT     min length to discard [" + min_flt + "]");
+		print("  -i INT     min SV length [" + min_size + "]");
+		print("  -x INT     max SV length [" + max_size + "]");
+		print("  -w INT     fuzzy windown size [" + win_size + "]");
+		print("  -e         print errors");
+		return;
+	}
+
+	function read_bed(fn) {
+		var buf = new Bytes();
+		var file = new File(fn);
+		var bed = {};
+		while (file.readline(buf) >= 0) {
+			var t = buf.toString().split("\t");
+			if (bed[t[0]] == null) bed[t[0]] = [];
+			bed[t[0]].push([parseInt(t[1]), parseInt(t[2])]);
+		}
+		file.close();
+		buf.destroy();
+		for (var x in bed) {
+			Interval.sort(bed[x]);
+			Interval.merge(bed[x]);
+			Interval.index_end(bed[x]);
+		}
+		return bed;
+	}
+
+	var bed = bed_fn != null? read_bed(bed_fn) : null;
+
+	function read_vcf(fn, bed) {
+		var buf = new Bytes();
+		var file = new File(fn);
+		var v = {};
+		while (file.readline(buf) >= 0) {
+			var m, t = buf.toString().split("\t");
+			if (t[0][0] == '#') continue;
+			if (bed != null && bed[t[0]] == null) continue;
+			if (t[4] == '<INV>' || t[4] == '<INVDUP>') continue; // no inversion
+			if (/[\[\]]/.test(t[4])) continue; // no break points
+			var st = parseInt(t[1]) - 1, en = st + t[3].length;
+			if ((m = /((;END)|(^END))=(\d+)/.exec(t[7])) != null)
+				en = parseInt(m[4]);
+			if (bed != null && Interval.find_ovlp(bed[t[0]], st, en).length == 0) continue;
+			// determine svlen
+			var s = t[4].split(","), max_del = 0, max_ins = 0;
+			for (var i = 0; i < s.length; ++i) {
+				var l = s[i].length - t[3].length;
+				if (l > 0)
+					max_ins = max_ins > l? max_ins : l;
+				else if (l < 0)
+					max_del = max_del > -l? max_del : -l;
+			}
+			if (max_ins < min_flt && max_del < min_flt) continue;
+			var svlen = max_ins > max_del? max_ins : -max_del;
+			if ((m = /((;SVLEN)|(^SVLEN))=(\d+)/.exec(t[7])) != null)
+				svlen = parseInt(m[4]);
+			var abslen = svlen > 0? svlen : -svlen;
+			if (abslen < min_flt || abslen > max_size) continue;
+			// insert
+			if (v[t[0]] == null) v[t[0]] = [];
+			v[t[0]].push([st, en, svlen, abslen]);
+		}
+		file.close();
+		buf.destroy();
+		for (var x in v) {
+			Interval.sort(v[x]);
+			Interval.index_end(v[x]);
+		}
+		return v;
+	}
+
+	function compare_vcf(v0, v1, label) {
+		var m = 0, n = 0;
+		for (var x in v1) {
+			var a1 = v1[x], a0 = v0[x];
+			for (var i = 0; i < a1.length; ++i) {
+				if (a1[i][3] < min_size) continue;
+				++n;
+				if (a0 == null) continue;
+				var st = a1[i][0] > win_size? a1[i][0] - win_size : 0;
+				b = Interval.find_ovlp(a0, st, a1[i][1] + win_size);
+				if (b.length > 0) ++m;
+				else if (print_err) print(label, x, a1[i].slice(0, 3).join("\t"));
+			}
+		}
+		return [n, m];
+	}
+
+	var v_base = read_vcf(args[getopt.ind+0], bed);
+	var v_call = read_vcf(args[getopt.ind+1], bed);
+	var fn = compare_vcf(v_call, v_base, 'FN');
+	var fp = compare_vcf(v_base, v_call, 'FP');
+	print('SN', fn[0], fn[1], (fn[1] / fn[0]).toFixed(6));
+	print('PC', fp[0], fp[1], (fp[1] / fp[0]).toFixed(6));
+	print('F1', ((fn[1] / fn[0] + fp[1] / fp[0]) / 2).toFixed(6));
+}
+
 /*************************
  ***** main function *****
  *************************/
@@ -2497,10 +2847,13 @@ function main(args)
 		print("");
 		print("  stat       collect basic mapping information in PAF/SAM");
 		print("  asmstat    collect basic assembly information");
-		print("  asmgene    evaluate gene completeness (EXPERIMENTAL)");
+		print("  asmgene    evaluate gene completeness");
+		print("  misjoin    evaluate large-scale misjoins");
 		print("  liftover   simplistic liftOver");
 		print("  call       call variants from asm-to-ref alignment with the cs tag");
 		print("  bedcov     compute the number of bases covered");
+		print("  vcfstat    VCF statistics");
+		print("  sveval     compare two SV callsets in VCF");
 		print("  version    print paftools.js version");
 		print("");
 		print("  mapeval    evaluate mapping accuracy using mason2/PBSIM-simulated FASTQ");
@@ -2520,6 +2873,7 @@ function main(args)
 	else if (cmd == 'stat') paf_stat(args);
 	else if (cmd == 'asmstat') paf_asmstat(args);
 	else if (cmd == 'asmgene') paf_asmgene(args);
+	else if (cmd == 'misjoin') paf_misjoin(args);
 	else if (cmd == 'liftover' || cmd == 'liftOver') paf_liftover(args);
 	else if (cmd == 'vcfpair') paf_vcfpair(args);
 	else if (cmd == 'call') paf_call(args);
@@ -2529,6 +2883,8 @@ function main(args)
 	else if (cmd == 'pbsim2fq') paf_pbsim2fq(args);
 	else if (cmd == 'junceval') paf_junceval(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 == 'version') print(paftools_version);
 	else throw Error("unrecognized command: " + cmd);
 }

mmpriv.h

@@ -4,6 +4,7 @@
 #include <assert.h>
 #include "minimap.h"
 #include "bseq.h"
+#include "kseq.h"
 
 #define MM_PARENT_UNSET   (-1)
 #define MM_PARENT_TMP_PRI (-2)
@@ -35,14 +36,6 @@
 extern "C" {
 #endif
 
-#ifndef KSTRING_T
-#define KSTRING_T kstring_t
-typedef struct __kstring_t {
-	unsigned l, m;
-	char *s;
-} kstring_t;
-#endif
-
 typedef struct {
 	int n_u, n_a;
 	uint64_t *u;
@@ -69,20 +62,21 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 void mm_idxopt_init(mm_idxopt_t *opt);
 const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f);
-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, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
+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);
 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);
+void mm_mark_alt(const mm_idx_t *mi, int n, mm_reg1_t *r);
 void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a);
 void mm_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 n, mm_reg1_t *r, int sub_diff, int hard_mask_level);
+void mm_set_parent(void *km, float mask_level, int mask_len, int n, mm_reg1_t *r, int sub_diff, int hard_mask_level, float alt_diff_frac);
 void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_, mm_reg1_t *r);
 void mm_select_sub_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);
 void mm_filter_regs(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs);
 void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs, mm128_t *a);
-void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r);
+void mm_hit_sort(void *km, int *n_regs, mm_reg1_t *r, float alt_diff_frac);
 void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr);
 
 void mm_est_err(const mm_idx_t *mi, int qlen, int n_regs, mm_reg1_t *regs, const mm128_t *a, int32_t n, const uint64_t *mini_pos);

options.c

@@ -1,4 +1,5 @@
 #include <stdio.h>
+#include <limits.h>
 #include "mmpriv.h"
 
 void mm_idxopt_init(mm_idxopt_t *opt)
@@ -24,8 +25,10 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->max_gap_ref = -1;
 	opt->max_chain_skip = 25;
 	opt->max_chain_iter = 5000;
+	opt->chain_gap_scale = 1.0f;
 
 	opt->mask_level = 0.5f;
+	opt->mask_len = INT_MAX;
 	opt->pri_ratio = 0.8f;
 	opt->best_n = 5;
 
@@ -34,6 +37,8 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->min_join_flank_sc = 1000;
 	opt->min_join_flank_ratio = 0.5f;
 
+	opt->alt_drop = 0.15f;
+
 	opt->a = 2, opt->b = 4, opt->q = 4, opt->e = 2, opt->q2 = 24, opt->e2 = 1;
 	opt->sc_ambi = 1;
 	opt->zdrop = 400, opt->zdrop_inv = 200;

python/cmappy.pxd

@@ -6,7 +6,7 @@ cdef extern from "minimap.h":
 	#
 	ctypedef struct mm_idxopt_t:
 		short k, w, flag, bucket_bits
-		int mini_batch_size
+		int64_t mini_batch_size
 		uint64_t batch_size
 
 	ctypedef struct mm_mapopt_t:
@@ -20,12 +20,15 @@ cdef extern from "minimap.h":
 		int max_chain_skip, max_chain_iter
 		int min_cnt
 		int min_chain_score
+		float chain_gap_scale
 		float mask_level
+		int mask_len
 		float pri_ratio
 		int best_n
 		int max_join_long, max_join_short
 		int min_join_flank_sc
 		float min_join_flank_ratio
+		float alt_drop
 		int a, b, q, e, q2, e2
 		int sc_ambi
 		int noncan
@@ -41,7 +44,7 @@ cdef extern from "minimap.h":
 		int32_t min_mid_occ
 		int32_t mid_occ
 		int32_t max_occ
-		int mini_batch_size
+		int64_t mini_batch_size
 		int64_t max_sw_mat
 		const char *split_prefix
 

python/mappy.pyx

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

setup.py

@@ -4,16 +4,6 @@ except ImportError:
 	from distutils.core import setup
 	from distutils.extension import Extension
 
-cmdclass = {}
-
-try:
-	from Cython.Build import build_ext
-except ImportError: # without Cython
-	module_src = 'python/mappy.c'
-else: # with Cython
-	module_src = 'python/mappy.pyx'
-	cmdclass['build_ext'] = build_ext
-
 import sys, platform
 
 sys.path.append('python')
@@ -33,7 +23,7 @@ def readme():
 
 setup(
 	name = 'mappy',
-	version = '2.17',
+	version = '2.18',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),
@@ -42,8 +32,8 @@ setup(
 	license = 'MIT',
 	keywords = 'sequence-alignment',
 	scripts = ['python/minimap2.py'],
-    ext_modules = [Extension('mappy',
+	ext_modules = [Extension('mappy',
-		sources = [module_src, 'align.c', 'bseq.c', 'chain.c', 'format.c', 'hit.c', 'index.c', 'pe.c', 'options.c',
+		sources = ['python/mappy.pyx', 'align.c', 'bseq.c', 'chain.c', 'format.c', 'hit.c', 'index.c', 'pe.c', 'options.c',
 				   'ksw2_extd2_sse.c', 'ksw2_exts2_sse.c', 'ksw2_extz2_sse.c', 'ksw2_ll_sse.c',
 				   'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c', 'esterr.c', 'splitidx.c'],
 		depends = ['minimap.h', 'bseq.h', 'kalloc.h', 'kdq.h', 'khash.h', 'kseq.h', 'ksort.h',
@@ -62,4 +52,4 @@ setup(
 		'Programming Language :: Python :: 3',
 		'Intended Audience :: Science/Research',
 		'Topic :: Scientific/Engineering :: Bio-Informatics'],
-	cmdclass = cmdclass)
+	setup_requires=["cython"])