Release minimap2-2.6 (r623)

* Don't preallocate sdust_buf or minizer list. kalloc should be fast enough -
  benchmarks needed to confirm.

* Fixed a memory leak caused by divergence estimate (post v2.5)

* Reset the kalloc buffer after mapping a long query. This reduces peak memory
  when large chunks of memory are allocated, at the cost of performance, though.

Makefile

@@ -1,7 +1,7 @@
 CFLAGS=		-g -Wall -O2 -Wc++-compat
 CPPFLAGS=	-DHAVE_KALLOC
 INCLUDES=
-OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o index.o chain.o align.o hit.o map.o format.o pe.o ksw2_ll_sse.o
+OBJS=		kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o index.o chain.o align.o hit.o map.o format.o pe.o esterr.o ksw2_ll_sse.o
 PROG=		minimap2
 PROG_EXTRA=	sdust minimap2-lite
 LIBS=		-lm -lz -lpthread

NEWS.md

@@ -1,3 +1,96 @@
+Release 2.6-r623 (12 December 2017)
+-----------------------------------
+
+This release adds several features and fixes two minor bugs:
+
+ * Optionally build an index without sequences. This helps to reduce the
+   peak memory for read overlapping and is automatically applied when
+   base-level alignment is not requested.
+
+ * Approximately estimate per-base sequence divergence (i.e. 1-identity)
+   without performing base-level alignment, using a MashMap-like method. The
+   estimate is written to a new dv:f tag.
+
+ * Reduced the number of tiny terminal exons in RNA-seq alignment. The current
+   setting is conservative. Increase --end-seed-pen to drop more such exons.
+
+ * Reduced the peak memory when aligning long query sequences.
+
+ * Fixed a bug that is caused by HPC minimizers longer than 256bp. This should
+   have no effect in practice, but it is recommended to rebuild HPC indices if
+   possible.
+
+ * Fixed a bug when identifying identical hits (#71). This should only affect
+   artifactual reference consisting of near identical sequences.
+
+For genomic sequences, minimap2 should give nearly identical alignments to
+v2.5, except the new dv:f tag.
+
+(2.6: 12 December 2017, r623)
+
+
+
+Release 2.5-r572 (11 November 2017)
+-----------------------------------
+
+This release fixes several bugs and brings a couple of minor improvements:
+
+ * Fixed a severe bug that leads to incorrect mapping coordinates in rare
+   corner cases.
+
+ * Fixed underestimated mapping quality for chimeric alignments when the whole
+   query sequence contain many repetitive minimizers, and for chimeric
+   alignments caused by Z-drop.
+
+ * Fixed two bugs in Python binding: incorrect strand field (#57) and incorrect
+   sequence names for Python3 (#55).
+
+ * Improved mapping accuracy for highly overlapping paired ends.
+
+ * Added option -Y to use soft clipping for supplementary alignments (#56).
+
+(2.5: 11 November 2017, r572)
+
+
+
+Release 2.4-r555 (6 November 2017)
+----------------------------------
+
+As is planned, this release focuses on fine tuning the base algorithm. Notable
+changes include
+
+ * Changed the mapping quality scale to match the scale of BWA-MEM. This makes
+   minimap2 and BWA-MEM achieve similar sensitivity-specificity balance on real
+   short-read data.
+
+ * Improved the accuracy of splice alignment by modeling one additional base
+   close to the GT-AG signal. This model is used by default with `-x splice`.
+   For SIRV control data, however, it is recommended to add `--splice-flank=no`
+   to disable this feature as the SIRV splice signals are slightly different.
+
+ * Tuned the parameters for Nanopore Direct RNA reads. The recommended command
+   line is `-axsplice -k14 -uf` (#46).
+
+ * Fixed a segmentation fault when aligning PacBio reads (#47 and #48). This
+   bug is very rare but it affects all versions of minimap2. It is also
+   recommended to re-index reference genomes created with `map-pb`. For human,
+   two minimizers in an old index are wrong.
+
+ * Changed option `-L` in sync with the final decision of hts-specs: a fake
+   CIGAR takes the form of `<readLen>S<refLen>N`. Note that `-L` only enables
+   future tools to recognize long CIGARs. It is not possible for older tools to
+   work with such alignments in BAM (#43 and #51).
+
+ * Fixed a tiny issue whereby minimap2 may waste 8 bytes per candidate
+   alignment.
+
+The minimap2 technical note hosted at arXiv has also been updated to reflect
+recent changes.
+
+(2.4: 6 November 2017, r555)
+
+
+
 Release 2.3-r531 (22 October 2017)
 ----------------------------------
 
@@ -26,7 +119,7 @@ This release come with many improvements and bug fixes:
 
 This release has implemented all the major features I planned five months ago,
 with the addition of spliced long-read alignment. The next couple of releases
-will focus on fine tuning of base algorithms.
+will focus on fine tuning of the base algorithms.
 
 (2.3: 22 October 2017, r531)
 

README.md

@@ -1,23 +1,25 @@
-[![Release](https://img.shields.io/badge/Release-v2.3-blue.svg?style=flat)](https://github.com/lh3/minimap2/releases)
-[![BioConda](https://img.shields.io/conda/vn/bioconda/minimap2.svg?style=flat)](https://anaconda.org/bioconda/minimap2)
+[![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)
-[![Python Version](https://img.shields.io/pypi/pyversions/mappy.svg?style=flat)](https://pypi.python.org/pypi/mappy)
-[![License](https://img.shields.io/badge/License-MIT-blue.svg?style=flat)](LICENSE.txt)
 [![Build Status](https://travis-ci.org/lh3/minimap2.svg?branch=master)](https://travis-ci.org/lh3/minimap2)
-[![Downloads](https://img.shields.io/github/downloads/lh3/minimap2/total.svg?style=flat)](https://github.com/lh3/minimap2/releases)
 ## <a name="started"></a>Getting Started
 ```sh
 git clone https://github.com/lh3/minimap2
 cd minimap2 && make
-# long reads against a reference genome
+# long sequences against a reference genome
 ./minimap2 -a test/MT-human.fa test/MT-orang.fa > test.sam
 # create an index first and then map
 ./minimap2 -d MT-human.mmi test/MT-human.fa
 ./minimap2 -a MT-human.mmi test/MT-orang.fa > test.sam
-# long-read overlap (no test data)
-./minimap2 -x ava-pb your-reads.fa your-reads.fa > overlaps.paf
-# spliced alignment (no test data)
-./minimap2 -ax splice ref.fa rna-seq-reads.fa > spliced.sam
+# use presets (no test data)
+./minimap2 -ax map-pb ref.fa pacbio.fq.gz > aln.sam       # PacBio genomic reads
+./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam         # Oxford Nanopore genomic reads
+./minimap2 -ax sr ref.fa read1.fa read2.fa > aln.sam      # short genomic paired-end reads
+./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads
+./minimap2 -ax splice -k14 -uf ref.fa reads.fa > aln.sam  # Nanopore Direct RNA-seq
+./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
 ```
@@ -34,7 +36,7 @@ man ./minimap2.1
     - [Map short accurate genomic reads](#short-genomic)
     - [Full genome/assembly alignment](#full-genome)
   - [Advanced features](#advanced)
-    - [Working CIGARs with >65535 operations](#long-cigar)
+    - [Working with >65535 CIGAR operations](#long-cigar)
     - [The cs optional tag](#cs)
     - [Evaluation scripts](#eval)
   - [Algorithm overview](#algo)
@@ -66,9 +68,9 @@ Detailed evaluations are available from the [minimap2 preprint][preprint].
 Minimap2 only works on x86-64 CPUs. You can acquire precompiled binaries from
 the [release page][release] with:
 ```sh
-wget --no-check-certificate -O- https://github.com/lh3/minimap2/releases/download/v2.2/minimap2-2.2_x64-linux.tar.bz2 \
+curl -L https://github.com/lh3/minimap2/releases/download/v2.6/minimap2-2.6_x64-linux.tar.bz2 \
   | tar -jxvf -
-./minimap2-2.2_x64-linux/minimap2
+./minimap2-2.6_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
@@ -132,14 +134,44 @@ Nanopore reads.
 #### <a name="map-long-splice"></a>Map long mRNA/cDNA reads
 
 ```sh
-minimap2 -ax splice ref.fa spliced.fq > aln.sam      # strand unknown
-minimap2 -ax splice -uf ref.fa spliced.fq > aln.sam  # assuming transcript strand
+minimap2 -ax splice -uf ref.fa iso-seq.fq > aln.sam          # PacBio Iso-seq/traditional cDNA
+minimap2 -ax splice ref.fa nanopore-cdna.fa > aln.sam        # Nanopore 2D cDNA-seq
+minimap2 -ax splice -uf -k14 ref.fa direct-rna.fq > aln.sam  # Nanopore Direct RNA-seq
+minimap2 -ax splice --splice-flank=no SIRV.fa SIRV-seq.fa    # mapping against SIRV control
 ```
-This command line has been tested on PacBio Iso-Seq reads and Nanopore 2D cDNA
-reads, and been shown to work with Nanopore 1D Direct RNA reads by others. Like
-typical RNA-seq mappers, minimap2 represents an intron with the `N` CIGAR
-operator. For spliced reads, minimap2 will try to infer the strand relative to
-transcript and may write the strand to the `ts` SAM/PAF tag.
+There are different long-read RNA-seq technologies, including tranditional
+full-length cDNA, EST, PacBio Iso-seq, Nanopore 2D cDNA-seq and Direct RNA-seq.
+They produce data of varying quality and properties. By default, `-x splice`
+assumes the read orientation relative to the transcript strand is unknown. It
+tries two rounds of alignment to infer the orientation and write the strand to
+the `ts` SAM/PAF tag if possible. For Iso-seq, Direct RNA-seq and tranditional
+full-length cDNAs, it would be desired to apply `-u f` to force minimap2 to
+consider the forward transcript strand only. This speeds up alignment with
+slight improvement to accuracy. For noisy Nanopore Direct RNA-seq reads, it is
+recommended to use a smaller k-mer size for increased sensitivity to the first
+or the last exons.
+
+Minimap2 rates an alignment by the score of the max-scoring sub-segment,
+*excluding* introns, and marks the best alignment as primary in SAM. When a
+spliced gene also has unspliced pseudogenes, minimap2 does not intentionally
+prefer spliced alignment, though in practice it more often marks the spliced
+alignment as the primary. By default, minimap2 outputs up to five secondary
+alignments (i.e. likely pseudogenes in the context of RNA-seq mapping). This
+can be tuned with option **-N**.
+
+For long RNA-seq reads, minimap2 may produce chimeric alignments potentially
+caused by gene fusions/structural variations or by an intron longer than the
+max intron length **-G** (200k by default). For now, it is not recommended to
+apply an excessively large **-G** as this slows down minimap2 and sometimes
+leads to false alignments.
+
+It is worth noting that by default `-x splice` prefers GT[A/G]..[C/T]AG
+over GT[C/T]..[A/G]AG, and then over other splicing signals. Considering
+one additional base improves the junction accuracy for noisy reads, but
+reduces the accuracy when aligning against the widely used SIRV control data.
+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.
 
 #### <a name="long-overlap"></a>Find overlaps between long reads
 
@@ -179,7 +211,7 @@ according to the sequence divergence.
 
 ### <a name="advanced"></a>Advanced features
 
-#### <a name="long-cigar"></a>Working CIGARs with >65535 operations
+#### <a name="long-cigar"></a>Working with >65535 CIGAR operations
 
 Due to a design flaw, BAM does not work with CIGAR strings with >65535
 operations (SAM and CRAM work). However, for ultra-long nanopore reads minimap2
@@ -191,10 +223,8 @@ To avoid this issue, you can add option `-L` at the minimap2 command line.
 This option moves a long CIGAR to the `CG` tag and leaves a fully clipped CIGAR
 at the SAM CIGAR column. Current tools that don't read CIGAR (e.g. merging and
 sorting) still work with such BAM records; tools that read CIGAR will
-effectively ignore these records. I have pull requests to the SAM spec, htslib,
-htsjdk, bedtools2, Rsamtools and igv.js. If they are accepted, future versions
-of these tools will seamlessly recognize long-cigar records generated by option
-`-L`.
+effectively ignore these records. It has been decided that future tools will
+will seamlessly recognize long-cigar records generated by option `-L`.
 
 **TD;DR**: if you work with ultra-long reads and use tools that only process
 BAM files, please add option `-L`.
@@ -290,7 +320,7 @@ highlighted in bold. The description may help to tune minimap2 parameters.
 
 ### <a name="help"></a>Getting help
 
-Manpage [minimap2.1](minimap2.1) provides detailed description of minimap2
+Manpage [minimap2.1][manpage] provides detailed description of minimap2
 command line options and optional tags. If you encounter bugs or have further
 questions or requests, you can raise an issue at the [issue page][issue].
 There is not a specific mailing list for the time being.
@@ -344,3 +374,4 @@ warmly welcomed.
 [mappyconda]: https://anaconda.org/bioconda/mappy
 [issue]: https://github.com/lh3/minimap2/issues
 [k8]: https://github.com/attractivechaos/k8
+[manpage]: https://lh3.github.io/minimap2/minimap2.html

align.c

@@ -1,6 +1,7 @@
 #include <assert.h>
 #include <string.h>
 #include <stdlib.h>
+#include <math.h>
 #include "minimap.h"
 #include "mmpriv.h"
 #include "ksw2.h"
@@ -110,14 +111,17 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 	}
 	if ((p->cigar[0]&0xf) == 1 || (p->cigar[0]&0xf) == 2) { // get rid of leading I or D
 		int32_t l = p->cigar[0] >> 4;
-		if ((p->cigar[0]&0xf) == 1) r->qs += l, *qshift = l;
-		else r->rs += l, *tshift = l;
+		if ((p->cigar[0]&0xf) == 1) {
+			if (r->rev) r->qe -= l;
+			else r->qs += l;
+			*qshift = l;
+		} else r->rs += l, *tshift = l;
 		--p->n_cigar;
 		memmove(p->cigar, p->cigar + 1, p->n_cigar * 4);
 	}
 }
 
-static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *qual, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e)
+static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e)
 {
 	uint32_t k, l, toff = 0, qoff = 0;
 	int32_t s = 0, max = 0, qshift, tshift;
@@ -130,27 +134,21 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *qu
 		uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4;
 		if (op == 0) { // match/mismatch
 			int n_ambi = 0, n_diff = 0;
-			float n_diff2 = 0.0f;
 			for (l = 0; l < len; ++l) {
 				int cq = qseq[qoff + l], ct = tseq[toff + l];
 				if (ct > 3 || cq > 3) ++n_ambi;
-				else if (ct != cq) {
-					++n_diff;
-					n_diff2 += qual == 0 || qual[qoff + l] >= 20? 1.0f : .05f * qual[qoff + l];
-				}
+				else if (ct != cq) ++n_diff;
 				s += mat[ct * 5 + cq];
 				if (s < 0) s = 0;
 				else max = max > s? max : s;
 			}
 			r->blen += len - n_ambi, r->mlen += len - (n_ambi + n_diff), p->n_ambi += n_ambi;
-			p->n_diff2 += n_diff2, p->blen2 += len - n_ambi;
 			toff += len, qoff += len;
 		} else if (op == 1) { // insertion
 			int n_ambi = 0;
 			for (l = 0; l < len; ++l)
 				if (qseq[qoff + l] > 3) ++n_ambi;
 			r->blen += len - n_ambi, p->n_ambi += n_ambi;
-			p->n_diff2 += 1.0f, ++p->blen2;
 			s -= q + e * len;
 			if (s < 0) s = 0;
 			qoff += len;
@@ -159,7 +157,6 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *qu
 			for (l = 0; l < len; ++l)
 				if (tseq[toff + l] > 3) ++n_ambi;
 			r->blen += len - n_ambi, p->n_ambi += n_ambi;
-			p->n_diff2 += 1.0f, ++p->blen2;
 			s -= q + e * len;
 			if (s < 0) s = 0;
 			toff += len;
@@ -226,7 +223,7 @@ static inline int mm_get_hplen_back(const mm_idx_t *mi, uint32_t rid, uint32_t x
 
 static inline void mm_adjust_minier(const mm_idx_t *mi, uint8_t *const qseq0[2], mm128_t *a, int32_t *r, int32_t *q)
 {
-	if (mi->is_hpc) {
+	if (mi->flag & MM_I_HPC) {
 		const uint8_t *qseq = qseq0[a->x>>63];
 		int i, c;
 		*q = (int32_t)a->y;
@@ -345,7 +342,49 @@ static void mm_max_stretch(const mm_mapopt_t *opt, const mm_reg1_t *r, const mm1
 	*as = max_i, *cnt = max_len;
 }
 
-static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], uint8_t *qual0[2], mm_reg1_t *r, mm_reg1_t *r2, int n_a, mm128_t *a, ksw_extz_t *ez, int splice_flag)
+static int mm_seed_ext_score(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, const int8_t mat[25], int qlen, uint8_t *qseq0[2], const mm128_t *a)
+{
+	uint8_t *qseq, *tseq;
+	int q_span = a->y>>32&0xff, qs, qe, rs, re, rid, score, q_off, t_off, ext_len = opt->anchor_ext_len;
+	void *qp;
+	rid = a->x<<1>>33;
+	re = (uint32_t)a->x + 1, rs = re - q_span;
+	qe = (uint32_t)a->y + 1, qs = qe - q_span;
+	rs = rs - ext_len > 0? rs - ext_len : 0;
+	qs = qs - ext_len > 0? qs - ext_len : 0;
+	re = re + ext_len < mi->seq[rid].len? re + ext_len : mi->seq[rid].len;
+	qe = qe + ext_len < qlen? qe + ext_len : qlen;
+	tseq = (uint8_t*)kmalloc(km, re - rs);
+	mm_idx_getseq(mi, rid, rs, re, tseq);
+	qseq = qseq0[a->x>>63] + qs;
+	qp = ksw_ll_qinit(km, 2, qe - qs, qseq, 5, mat);
+	score = ksw_ll_i16(qp, re - rs, tseq, opt->q, opt->e, &q_off, &t_off);
+	kfree(km, tseq);
+	kfree(km, qp);
+	return score;
+}
+
+static void mm_fix_bad_ends_splice(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, const mm_reg1_t *r, const int8_t mat[25], int qlen, uint8_t *qseq0[2], const mm128_t *a, int *as1, int *cnt1)
+{ // this assumes a very crude k-mer based mode; it is not necessary to use a good model just for filtering bounary exons
+	int score;
+	double log_gap;
+	*as1 = r->as, *cnt1 = r->cnt;
+	if (r->cnt < 3) return;
+	log_gap = log((int32_t)a[r->as + 1].x - (int32_t)a[r->as].x);
+	if ((a[r->as].y>>32&0xff) < log_gap + opt->anchor_ext_shift) {
+		score = mm_seed_ext_score(km, opt, mi, mat, qlen, qseq0, &a[r->as]);
+		if ((double)score / mat[0] < log_gap + opt->anchor_ext_shift) // a more exact format is "score < log_4(gap) + shift"
+			++(*as1), --(*cnt1);
+	}
+	log_gap = log((int32_t)a[r->as + r->cnt - 1].x - (int32_t)a[r->as + r->cnt - 2].x);
+	if ((a[r->as + r->cnt - 1].y>>32&0xff) < log_gap + opt->anchor_ext_shift) {
+		score = mm_seed_ext_score(km, opt, mi, mat, qlen, qseq0, &a[r->as + r->cnt - 1]);
+		if ((double)score / mat[0] < log_gap + opt->anchor_ext_shift)
+			--(*cnt1);
+	}
+}
+
+static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], mm_reg1_t *r, mm_reg1_t *r2, int n_a, mm128_t *a, ksw_extz_t *ez, int splice_flag)
 {
 	int is_sr = !!(opt->flag & MM_F_SR), is_splice = !!(opt->flag & MM_F_SPLICE);
 	int32_t rid = a[r->as].x<<1>>33, rev = a[r->as].x>>63, as1, cnt1;
@@ -355,23 +394,25 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	int32_t rs1, qs1, re1, qe1;
 	int8_t mat[25];
 
-	if (is_sr) assert(!mi->is_hpc); // HPC won't work with SR because with HPC we can't easily tell if there is a gap
+	if (is_sr) assert(!(mi->flag & MM_I_HPC)); // HPC won't work with SR because with HPC we can't easily tell if there is a gap
 
 	r2->cnt = 0;
 	if (r->cnt == 0) return;
 	ksw_gen_simple_mat(5, mat, opt->a, opt->b);
 	bw = (int)(opt->bw * 1.5 + 1.);
 
-	if (is_sr && !mi->is_hpc) {
+	if (is_sr && !(mi->flag & MM_I_HPC)) {
 		mm_max_stretch(opt, r, a, &as1, &cnt1);
 		rs = (int32_t)a[as1].x + 1 - (int32_t)(a[as1].y>>32&0xff);
 		qs = (int32_t)a[as1].y + 1 - (int32_t)(a[as1].y>>32&0xff);
 		re = (int32_t)a[as1+cnt1-1].x + 1;
 		qe = (int32_t)a[as1+cnt1-1].y + 1;
 	} else {
-		if (!is_splice)
+		if (is_splice) {
+			mm_fix_bad_ends_splice(km, opt, mi, r, mat, qlen, qseq0, a, &as1, &cnt1);
+		} else {
 			mm_fix_bad_ends(r, a, opt->bw, &as1, &cnt1);
-		else as1 = r->as, cnt1 = r->cnt;
+		}
 		mm_filter_bad_seeds(km, as1, cnt1, a, 10, 40, opt->max_gap>>1, 10);
 		mm_adjust_minier(mi, qseq0, &a[as1], &rs, &qs);
 		mm_adjust_minier(mi, qseq0, &a[as1 + cnt1 - 1], &re, &qe);
@@ -381,6 +422,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	if (is_splice) {
 		if (splice_flag & MM_F_SPLICE_FOR) extra_flag |= rev? KSW_EZ_SPLICE_REV : KSW_EZ_SPLICE_FOR;
 		if (splice_flag & MM_F_SPLICE_REV) extra_flag |= rev? KSW_EZ_SPLICE_FOR : KSW_EZ_SPLICE_REV;
+		if (opt->flag & MM_F_SPLICE_FLANK) extra_flag |= KSW_EZ_SPLICE_FLANK;
 	}
 
 	/* Look for the start and end of regions to perform DP. This sounds easy
@@ -473,7 +515,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 
 	for (i = is_sr? cnt1 - 1 : 1; i < cnt1; ++i) { // gap filling
 		if ((a[as1+i].y & (MM_SEED_IGNORE|MM_SEED_TANDEM)) && i != cnt1 - 1) continue;
-		if (is_sr && !mi->is_hpc) {
+		if (is_sr && !(mi->flag & MM_I_HPC)) {
 			re = (int32_t)a[as1 + i].x + 1;
 			qe = (int32_t)a[as1 + i].y + 1;
 		} else mm_adjust_minier(mi, qseq0, &a[as1 + i], &re, &qe);
@@ -536,7 +578,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	assert(re1 - rs1 <= re0 - rs0);
 	if (r->p) {
 		mm_idx_getseq(mi, rid, rs1, re1, tseq);
-		mm_update_extra(r, &qseq0[r->rev][qs1], qual0[r->rev]? &qual0[r->rev][qs1] : 0, tseq, mat, opt->q, opt->e);
+		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e);
 		if (rev && r->p->trans_strand)
 			r->p->trans_strand ^= 3; // flip to the read strand
 	}
@@ -544,10 +586,10 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	kfree(km, tseq);
 }
 
-static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], uint8_t *qual0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez)
+static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez)
 {
 	int tl, ql, score, ret = 0, q_off, t_off;
-	uint8_t *tseq, *qseq, *qual;
+	uint8_t *tseq, *qseq;
 	int8_t mat[25];
 	void *qp;
 
@@ -565,7 +607,6 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	tseq = (uint8_t*)kmalloc(km, tl);
 	mm_idx_getseq(mi, r1->rid, r1->re, r2->rs, tseq);
 	qseq = &qseq0[!r1->rev][qlen - r2->qs];
-	qual = qual0[!r1->rev]? &qseq0[!r1->rev][qlen - r2->qs] : 0;
 
 	mm_seq_rev(ql, qseq);
 	mm_seq_rev(tl, tseq);
@@ -585,9 +626,10 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	r_inv->inv = 1;
 	r_inv->rev = !r1->rev;
 	r_inv->rid = r1->rid;
+	r_inv->div = -1.0f;
 	r_inv->qs = r1->qe + q_off, r_inv->qe = r_inv->qs + ez->max_q + 1;
 	r_inv->rs = r1->re + t_off, r_inv->re = r_inv->rs + ez->max_t + 1;
-	mm_update_extra(r_inv, &qseq[q_off], qual? &qual[q_off] : 0, &tseq[t_off], mat, opt->q, opt->e);
+	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e);
 	ret = 1;
 end_align1_inv:
 	kfree(km, tseq);
@@ -604,11 +646,11 @@ static inline mm_reg1_t *mm_insert_reg(const mm_reg1_t *r, int i, int *n_regs, m
 	return regs;
 }
 
-mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, const char *qual, int *n_regs_, mm_reg1_t *regs, mm128_t *a)
+mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a)
 {
 	extern unsigned char seq_nt4_table[256];
 	int32_t i, n_regs = *n_regs_, n_a;
-	uint8_t *qseq0[2], *qual0[2];
+	uint8_t *qseq0[2];
 	ksw_extz_t ez;
 
 	// encode the query sequence
@@ -618,12 +660,6 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 		qseq0[0][i] = seq_nt4_table[(uint8_t)qstr[i]];
 		qseq0[1][qlen - 1 - i] = qseq0[0][i] < 4? 3 - qseq0[0][i] : 4;
 	}
-	if (qual) {
-		qual0[0] = (uint8_t*)kmalloc(km, qlen * 2);
-		qual0[1] = qual0[0] + qlen;
-		for (i = 0; i < qlen; ++i)
-			qual0[0][i] = qual0[1][qlen - 1 - i] = qual[i] - 33;
-	} else qual0[0] = qual0[1] = 0;
 
 	// align through seed hits
 	n_a = mm_squeeze_a(km, n_regs, regs, a);
@@ -634,8 +670,8 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 			mm_reg1_t s[2], s2[2];
 			int which, trans_strand;
 			s[0] = s[1] = regs[i];
-			mm_align1(km, opt, mi, qlen, qseq0, qual0, &s[0], &s2[0], n_a, a, &ez, MM_F_SPLICE_FOR);
-			mm_align1(km, opt, mi, qlen, qseq0, qual0, &s[1], &s2[1], n_a, a, &ez, MM_F_SPLICE_REV);
+			mm_align1(km, opt, mi, qlen, qseq0, &s[0], &s2[0], n_a, a, &ez, MM_F_SPLICE_FOR);
+			mm_align1(km, opt, mi, qlen, qseq0, &s[1], &s2[1], n_a, a, &ez, MM_F_SPLICE_REV);
 			if (s[0].p->dp_score > s[1].p->dp_score) which = 0, trans_strand = 1;
 			else if (s[0].p->dp_score < s[1].p->dp_score) which = 1, trans_strand = 2;
 			else trans_strand = 3, which = (qlen + s[0].p->dp_score) & 1; // randomly choose a strand, effectively
@@ -648,13 +684,13 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 			}
 			regs[i].p->trans_strand = trans_strand;
 		} else { // one round of alignment
-			mm_align1(km, opt, mi, qlen, qseq0, qual0, &regs[i], &r2, n_a, a, &ez, opt->flag);
+			mm_align1(km, opt, mi, qlen, qseq0, &regs[i], &r2, n_a, a, &ez, opt->flag);
 			if (opt->flag&MM_F_SPLICE)
 				regs[i].p->trans_strand = opt->flag&MM_F_SPLICE_FOR? 1 : 2;
 		}
 		if (r2.cnt > 0) regs = mm_insert_reg(&r2, i, &n_regs, regs);
 		if (!(opt->flag&MM_F_SPLICE) && !(opt->flag&MM_F_SR) && i > 0) { // don't try inversion alignment for -xsplice or -xsr
-			if (mm_align1_inv(km, opt, mi, qlen, qseq0, qual0, &regs[i-1], &regs[i], &r2, &ez)) {
+			if (mm_align1_inv(km, opt, mi, qlen, qseq0, &regs[i-1], &regs[i], &r2, &ez)) {
 				regs = mm_insert_reg(&r2, i, &n_regs, regs);
 				++i; // skip the inserted INV alignment
 			}
@@ -662,7 +698,6 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
 	}
 	*n_regs_ = n_regs;
 	kfree(km, qseq0[0]);
-	if (qual0[0]) kfree(km, qual0[0]);
 	kfree(km, ez.cigar);
 	mm_filter_regs(km, opt, n_regs_, regs);
 	mm_hit_sort_by_dp(km, n_regs_, regs);

chain.c

@@ -49,7 +49,7 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 			int64_t dr = ri - a[j].x;
 			int32_t dq = qi - (int32_t)a[j].y, dd, sc, log_dd;
 			int32_t sidj = (a[j].y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
-			if (dr == 0 || dq <= 0) continue;
+			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;
 			dd = dr > dq? dr - dq : dq - dr;
 			if (sidi == sidj && dd > bw) continue;
@@ -61,7 +61,8 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 				int c_log, c_lin;
 				c_lin = (int)(dd * .01 * avg_qspan);
 				c_log = log_dd;
-				if (dr > dq || sidi != sidj) sc -= c_lin < c_log? c_lin : c_log;
+				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 sc -= c_lin + (c_log>>1);
 			} else sc -= (int)(dd * .01 * avg_qspan) + (log_dd>>1);
 			sc += f[j];

esterr.c

@@ -0,0 +1,64 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+#include "mmpriv.h"
+
+static inline int32_t get_for_qpos(int32_t qlen, const mm128_t *a)
+{
+	int32_t x = (int32_t)a->y;
+	int32_t q_span = a->y>>32 & 0xff;
+	if (a->x>>63)
+		x = qlen - 1 - (x + 1 - q_span); // revert the position to the forward strand of query
+	return x;
+}
+
+static int get_mini_idx(int qlen, const mm128_t *a, int32_t n, const uint64_t *mini_pos)
+{
+	int32_t x, L = 0, R = n - 1;
+	x = get_for_qpos(qlen, a);
+	while (L <= R) { // binary search
+		int32_t m = ((uint64_t)L + R) >> 1;
+		int32_t y = (int32_t)mini_pos[m];
+		if (y < x) L = m + 1;
+		else if (y > x) R = m - 1;
+		else return m;
+	}
+	return -1;
+}
+
+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)
+{
+	int i;
+	uint64_t sum_k = 0;
+	float avg_k;
+
+	if (n == 0) return;
+	for (i = 0; i < n; ++i)
+		sum_k += mini_pos[i] >> 32 & 0xff;
+	avg_k = (float)sum_k / n;
+
+	for (i = 0; i < n_regs; ++i) {
+		mm_reg1_t *r = &regs[i];
+		int32_t st, en, j, k, n_match, n_tot, l_ref;
+		r->div = -1.0f;
+		if (r->cnt == 0) continue;
+		st = en = get_mini_idx(qlen, r->rev? &a[r->as + r->cnt - 1] : &a[r->as], n, mini_pos);
+		if (st < 0) {
+			if (mm_verbose >= 2)
+				fprintf(stderr, "[WARNING] logic inconsistency in mm_est_err(). Please contact the developer.\n");
+			continue;
+		}
+		l_ref = mi->seq[r->rid].len;
+		for (k = 1, j = st + 1, n_match = 1; j < n && k < r->cnt; ++j) {
+			int32_t x;
+			x = get_for_qpos(qlen, r->rev? &a[r->as + r->cnt - 1 - k] : &a[r->as + k]);
+			if (x == (int32_t)mini_pos[j])
+				++k, en = j, ++n_match;
+		}
+		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;
+	}
+}

format.c

@@ -202,8 +202,9 @@ static void write_cs(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_
 
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 {
-	int type = r->inv? 'I' : r->id == r->parent? 'P' : 'S';
-	if (r->iden_flt) mm_sprintf_lite(s, "\tom:i:%d", r->mapq);
+	int type;
+	if (r->id == r->parent) type = r->inv? 'I' : 'P';
+	else type = r->inv? 'i' : 'S';
 	if (r->p) {
 		mm_sprintf_lite(s, "\tNM:i:%d\tms:i:%d\tAS:i:%d\tnn:i:%d", r->blen - r->mlen + r->p->n_ambi, r->p->dp_max, r->p->dp_score, r->p->n_ambi);
 		if (r->p->trans_strand == 1 || r->p->trans_strand == 2)
@@ -211,6 +212,12 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 	}
 	mm_sprintf_lite(s, "\ttp:A:%c\tcm:i:%d\ts1:i:%d", type, r->cnt, r->score);
 	if (r->parent == r->id) mm_sprintf_lite(s, "\ts2:i:%d", r->subsc);
+	if (r->div >= 0.0f && r->div <= 1.0f) {
+		char buf[8];
+		if (r->div == 0.0f) buf[0] = '0', buf[1] = 0;
+		else sprintf(buf, "%.4f", r->div);
+		mm_sprintf_lite(s, "\tdv:f:%s", buf);
+	}
 	if (r->split) mm_sprintf_lite(s, "\tzd:i:%d", r->split);
 }
 
@@ -258,7 +265,7 @@ static inline const mm_reg1_t *get_sam_pri(int n_regs, const mm_reg1_t *regs)
 	return NULL;
 }
 
-static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, const mm_reg1_t *r)
+static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, const mm_reg1_t *r, int opt_flag)
 {
 	if (r->p == 0) {
 		mm_sprintf_lite(s, "*");
@@ -267,14 +274,14 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
 		clip_len[0] = r->rev? qlen - r->qe : r->qs;
 		clip_len[1] = r->rev? r->qs : qlen - r->qe;
 		if (in_tag) {
-			int clip_char = (sam_flag&0x800)? 5 : 4;
+			int clip_char = (sam_flag&0x800) && !(opt_flag&MM_F_SOFTCLIP)? 5 : 4;
 			mm_sprintf_lite(s, "\tCG:B:I");
 			if (clip_len[0]) mm_sprintf_lite(s, ",%u", clip_len[0]<<4|clip_char);
 			for (k = 0; k < r->p->n_cigar; ++k)
 				mm_sprintf_lite(s, ",%u", r->p->cigar[k]);
 			if (clip_len[1]) mm_sprintf_lite(s, ",%u", clip_len[1]<<4|clip_char);
 		} else {
-			int clip_char = (sam_flag&0x800)? 'H' : 'S';
+			int clip_char = (sam_flag&0x800) && !(opt_flag&MM_F_SOFTCLIP)? 'H' : 'S';
 			if (clip_len[0]) mm_sprintf_lite(s, "%d%c", clip_len[0], clip_char);
 			for (k = 0; k < r->p->n_cigar; ++k)
 				mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDN"[r->p->cigar[k]&0xf]);
@@ -336,9 +343,8 @@ void mm_write_sam2(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 {
-		int mapq = !r->iden_flt? r->mapq : r->mapq < 3? r->mapq : 3;
 		this_rid = r->rid, this_pos = r->rs, this_rev = r->rev;
-		mm_sprintf_lite(s, "\t%s\t%d\t%d\t", mi->seq[r->rid].name, r->rs+1, mapq);
+		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;
 			if (r->qs != 0) ++n_cigar;
@@ -346,8 +352,11 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 			if (n_cigar > max_bam_cigar_op)
 				cigar_in_tag = 1;
 		}
-		if (cigar_in_tag) mm_sprintf_lite(s, "%dS", t->l_seq);
-		else write_sam_cigar(s, flag, 0, t->l_seq, r);
+		if (cigar_in_tag) {
+			if (flag & 0x100) mm_sprintf_lite(s, "0S"); // secondary alignment
+			else if (flag & 0x800) mm_sprintf_lite(s, "%dS", r->re - r->rs); // supplementary alignment
+			else mm_sprintf_lite(s, "%dS", t->l_seq);
+		} else write_sam_cigar(s, flag, 0, t->l_seq, r, opt_flag);
 	}
 
 	// write mate positions
@@ -380,7 +389,7 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 		if (t->qual) sam_write_sq(s, t->qual, t->l_seq, 0, 0);
 		else mm_sprintf_lite(s, "*");
 	} else {
-		if ((flag & 0x900) == 0) {
+		if ((flag & 0x900) == 0 || (opt_flag & MM_F_SOFTCLIP)) {
 			sam_write_sq(s, t->seq, t->l_seq, r->rev, r->rev);
 			mm_sprintf_lite(s, "\t");
 			if (t->qual) sam_write_sq(s, t->qual, t->l_seq, r->rev, 0);
@@ -428,7 +437,7 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 		if (r->p && (opt_flag & MM_F_OUT_CS))
 			write_cs(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG));
 		if (cigar_in_tag)
-			write_sam_cigar(s, flag, 1, t->l_seq, r);
+			write_sam_cigar(s, flag, 1, t->l_seq, r, opt_flag);
 	}
 
 	s->s[s->l] = 0; // we always have room for an extra byte (see str_enlarge)

hit.c

@@ -76,10 +76,11 @@ mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u,
 		mm_reg1_t *ri = &r[i];
 		ri->id = i;
 		ri->parent = MM_PARENT_UNSET;
-		ri->score = z[i].x >> 32;
+		ri->score = ri->score0 = z[i].x >> 32;
 		ri->hash = (uint32_t)z[i].x;
 		ri->cnt = (int32_t)z[i].y;
 		ri->as = z[i].y >> 32;
+		ri->div = -1.0f;
 		mm_reg_set_coor(ri, qlen, a);
 	}
 	kfree(km, z);
@@ -147,7 +148,7 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff
 				ri->parent = rp->parent;
 				rp->subsc = rp->subsc > ri->score? rp->subsc : ri->score;
 				if (ri->cnt >= rp->cnt) cnt_sub = 1;
-				if (rp->p && ri->p && (rp->rs != ri->rs || rp->re != ri->re || ol != min)) { // the last condition excludes identical hits after DP
+				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;
 					if (rp->p->dp_max - ri->p->dp_max <= sub_diff) cnt_sub = 1;
 				}
@@ -234,7 +235,7 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_,
 			if (p == i || r[i].inv) { // primary or inversion
 				r[k++] = r[i];
 			} else if ((r[i].score >= r[p].score * pri_ratio || r[i].score + min_diff >= r[p].score) && n_2nd < best_n) {
-				if (!(r[i].qs == r[p].qs && r[i].qe == r[p].qe && r[i].rs == r[p].rs && r[i].re == r[p].re)) // not identical hits
+				if (!(r[i].qs == r[p].qs && r[i].qe == r[p].qe && r[i].rid == r[p].rid && r[i].rs == r[p].rs && r[i].re == r[p].re)) // not identical hits
 					r[k++] = r[i], ++n_2nd;
 				else if (r[i].p) free(r[i].p);
 			} else if (r[i].p) free(r[i].p);
@@ -264,45 +265,6 @@ void mm_filter_regs(void *km, const mm_mapopt_t *opt, int *n_regs, mm_reg1_t *re
 	*n_regs = k;
 }
 
-void mm_filter_by_identity(void *km, int n_regs, mm_reg1_t *regs, float min_iden, int qlen, const char *qual) // TODO: make sure it is not beyond the ends of contigs
-{
-	int i, j, n_aux = 0, en, blen = 0;
-	uint64_t *aux;
-	float n_diff = 0.0f;
-	if (n_regs <= 0) return;
-	for (i = 0; i < n_regs; ++i)
-		if (regs[i].id == regs[i].parent && regs[i].pe_thru) // sequenced through the fragment; don't filter
-			return;
-	for (i = 0; i < n_regs; ++i)
-		if (regs[i].id == regs[i].parent)
-			++n_aux;
-	assert(n_aux >= 1);
-	aux = (uint64_t*)kmalloc(km, n_aux * 8);
-	for (i = 0, n_aux = 0; i < n_regs; ++i)
-		if (regs[i].id == regs[i].parent)
-			aux[n_aux++] = (uint64_t)regs[i].qs<<32 | i;
-	radix_sort_64(aux, aux + n_aux);
-	for (i = 0, en = 0; i < n_aux; ++i) {
-		mm_reg1_t *r = &regs[(int32_t)aux[i]];
-		if (r->qs > en) {
-			for (j = en; j < r->qs; ++j)
-				n_diff += qual == 0 || qual[j] >= 53? .25f : .05f * .25f * (qual[j] - 33);
-			blen += r->qs - en;
-		}
-		assert(r->p);
-		blen += r->p->blen2;
-		n_diff += r->p->n_diff2;
-		en = en > r->qe? en : r->qe;
-	}
-	for (j = en; j < qlen; ++j)
-		n_diff += qual == 0 || qual[j] >= 53? .25f : .05f * .25f * (qual[j] - 33);
-	blen += qlen - en;
-	kfree(km, aux);
-	if (1.0f - n_diff / blen < min_iden)
-		for (i = 0; i < n_regs; ++i)
-			regs[i].iden_flt = 1;
-}
-
 int mm_squeeze_a(void *km, int n_regs, mm_reg1_t *regs, mm128_t *a)
 { // squeeze out regions in a[] that are not referenced by regs[]
 	int i, as = 0;
@@ -442,26 +404,43 @@ void mm_seg_free(void *km, int n_segs, mm_seg_t *segs)
 	kfree(km, segs);
 }
 
-void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len)
+void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len, int is_sr)
 {
 	static const float q_coef = 40.0f;
+	int64_t sum_sc = 0;
+	float uniq_ratio;
 	int i;
+	for (i = 0; i < n_regs; ++i)
+		if (regs[i].parent == regs[i].id)
+			sum_sc += regs[i].score;
+	uniq_ratio = (float)sum_sc / (sum_sc + rep_len);
 	for (i = 0; i < n_regs; ++i) {
 		mm_reg1_t *r = &regs[i];
 		if (r->inv) {
 			r->mapq = 0;
 		} else if (r->parent == r->id) {
 			int mapq, subsc;
-			float pen_s1 = (r->score > 100? 1.0f : 0.01f * r->score) * ((float)r->score / (r->score + rep_len));
+			float pen_s1 = (r->score > 100? 1.0f : 0.01f * r->score) * uniq_ratio;
 			float pen_cm = r->cnt > 10? 1.0f : 0.1f * r->cnt;
 			pen_cm = pen_s1 < pen_cm? pen_s1 : pen_cm;
 			subsc = r->subsc > min_chain_sc? r->subsc : min_chain_sc;
 			if (r->p && r->p->dp_max2 > 0 && r->p->dp_max > 0) {
 				float identity = (float)r->mlen / r->blen;
-				int mapq_alt = (int)(6.02f * identity * identity * (r->p->dp_max - r->p->dp_max2) / match_sc + .499f); // BWA-MEM like mapQ, mostly for short reads
-				mapq = (int)(identity * pen_cm * q_coef * (1. - (float)r->p->dp_max2 * subsc / r->p->dp_max / r->score) * logf(r->score)); // more for long reads
-				mapq = mapq < mapq_alt? mapq : mapq_alt; // in case the long-read heuristic fails
-			} else mapq = (int)(pen_cm * q_coef * (1. - (float)subsc / r->score) * logf(r->score));
+				float x = (float)r->p->dp_max2 * subsc / r->p->dp_max / r->score0;
+				mapq = (int)(identity * pen_cm * q_coef * (1.0f - x * x) * logf((float)r->p->dp_max / match_sc));
+				if (!is_sr) {
+					int mapq_alt = (int)(6.02f * identity * identity * (r->p->dp_max - r->p->dp_max2) / match_sc + .499f); // BWA-MEM like mapQ, mostly for short reads
+					mapq = mapq < mapq_alt? mapq : mapq_alt; // in case the long-read heuristic fails
+				}
+			} else {
+				float x = (float)subsc / r->score0;
+				if (r->p) {
+					float identity = (float)r->mlen / r->blen;
+					mapq = (int)(identity * pen_cm * q_coef * (1.0f - x) * logf((float)r->p->dp_max / match_sc));
+				} else {
+					mapq = (int)(pen_cm * q_coef * (1.0f - x) * logf(r->score));
+				}
+			}
 			mapq -= (int)(4.343f * logf(r->n_sub + 1) + .499f);
 			mapq = mapq > 0? mapq : 0;
 			r->mapq = mapq < 60? mapq : 60;

index.c

@@ -31,19 +31,19 @@ typedef struct mm_idx_bucket_s {
 void mm_idxopt_init(mm_idxopt_t *opt)
 {
 	memset(opt, 0, sizeof(mm_idxopt_t));
-	opt->k = 15, opt->w = 10, opt->is_hpc = 0;
+	opt->k = 15, opt->w = 10, opt->flag = 0;
 	opt->bucket_bits = 14;
 	opt->mini_batch_size = 50000000;
 	opt->batch_size = 4000000000ULL;
 }
 
-mm_idx_t *mm_idx_init(int w, int k, int b, int is_hpc)
+mm_idx_t *mm_idx_init(int w, int k, int b, int flag)
 {
 	mm_idx_t *mi;
 	if (k*2 < b) b = k * 2;
 	if (w < 1) w = 1;
 	mi = (mm_idx_t*)calloc(1, sizeof(mm_idx_t));
-	mi->w = w, mi->k = k, mi->b = b, mi->is_hpc = is_hpc;
+	mi->w = w, mi->k = k, mi->b = b, mi->flag = flag;
 	mi->B = (mm_idx_bucket_t*)calloc(1<<b, sizeof(mm_idx_bucket_t));
 	if (!(mm_dbg_flag & 1)) mi->km = km_init();
 	return mi;
@@ -89,7 +89,7 @@ void mm_idx_stat(const mm_idx_t *mi)
 {
 	int i, n = 0, n1 = 0;
 	uint64_t sum = 0, len = 0;
-	fprintf(stderr, "[M::%s] kmer size: %d; skip: %d; is_HPC: %d; #seq: %d\n", __func__, mi->k, mi->w, mi->is_hpc, mi->n_seq);
+	fprintf(stderr, "[M::%s] kmer size: %d; skip: %d; is_hpc: %d; #seq: %d\n", __func__, mi->k, mi->w, mi->flag&MM_I_HPC, mi->n_seq);
 	for (i = 0; i < mi->n_seq; ++i)
 		len += mi->seq[i].len;
 	for (i = 0; i < 1<<mi->b; ++i)
@@ -214,7 +214,7 @@ static void mm_idx_post(mm_idx_t *mi, int n_threads)
 #include "bseq.h"
 
 typedef struct {
-	int mini_batch_size, keep_name;
+	int mini_batch_size;
 	uint64_t batch_size, sum_len;
 	mm_bseq_file_t *fp;
 	mm_idx_t *mi;
@@ -246,7 +246,6 @@ static void *worker_pipeline(void *shared, int step, void *in)
 		s->seq = mm_bseq_read(p->fp, p->mini_batch_size, 0, &s->n_seq); // read a mini-batch
 		if (s->seq) {
 			uint32_t old_m, m;
-			uint64_t sum_len, old_max_len, max_len;
 			assert((uint64_t)p->mi->n_seq + s->n_seq <= UINT32_MAX); // to prevent integer overflow
 			// make room for p->mi->seq
 			old_m = p->mi->n_seq, m = p->mi->n_seq + s->n_seq;
@@ -254,30 +253,34 @@ static void *worker_pipeline(void *shared, int step, void *in)
 			if (old_m != m)
 				p->mi->seq = (mm_idx_seq_t*)krealloc(p->mi->km, p->mi->seq, m * sizeof(mm_idx_seq_t));
 			// make room for p->mi->S
-			for (i = 0, sum_len = 0; i < s->n_seq; ++i) sum_len += s->seq[i].l_seq;
-			old_max_len = (p->sum_len + 7) / 8;
-			max_len = (p->sum_len + sum_len + 7) / 8;
-			kroundup64(old_max_len); kroundup64(max_len);
-			if (old_max_len != max_len) {
-				p->mi->S = (uint32_t*)realloc(p->mi->S, max_len * 4);
-				memset(&p->mi->S[old_max_len], 0, 4 * (max_len - old_max_len));
+			if (!(p->mi->flag & MM_I_NO_SEQ)) {
+				uint64_t sum_len, old_max_len, max_len;
+				for (i = 0, sum_len = 0; i < s->n_seq; ++i) sum_len += s->seq[i].l_seq;
+				old_max_len = (p->sum_len + 7) / 8;
+				max_len = (p->sum_len + sum_len + 7) / 8;
+				kroundup64(old_max_len); kroundup64(max_len);
+				if (old_max_len != max_len) {
+					p->mi->S = (uint32_t*)realloc(p->mi->S, max_len * 4);
+					memset(&p->mi->S[old_max_len], 0, 4 * (max_len - old_max_len));
+				}
 			}
 			// populate p->mi->seq
 			for (i = 0; i < s->n_seq; ++i) {
 				mm_idx_seq_t *seq = &p->mi->seq[p->mi->n_seq];
 				uint32_t j;
-				if (p->keep_name) {
-					assert(strlen(s->seq[i].name) <= 254); // a long query name breaks BAM
+				if (!(p->mi->flag & MM_I_NO_NAME)) {
 					seq->name = (char*)kmalloc(p->mi->km, strlen(s->seq[i].name) + 1);
 					strcpy(seq->name, s->seq[i].name);
 				} else seq->name = 0;
 				seq->len = s->seq[i].l_seq;
 				seq->offset = p->sum_len;
 				// copy the sequence
-				for (j = 0; j < seq->len; ++j) { // TODO: this is not the fastest way, but let's first see if speed matters here
-					uint64_t o = p->sum_len + j;
-					int c = seq_nt4_table[(uint8_t)s->seq[i].seq[j]];
-					mm_seq4_set(p->mi->S, o, c);
+				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
+						uint64_t o = p->sum_len + j;
+						int c = seq_nt4_table[(uint8_t)s->seq[i].seq[j]];
+						mm_seq4_set(p->mi->S, o, c);
+					}
 				}
 				// update p->sum_len and p->mi->n_seq
 				p->sum_len += seq->len;
@@ -289,7 +292,10 @@ static void *worker_pipeline(void *shared, int step, void *in)
         step_t *s = (step_t*)in;
 		for (i = 0; i < s->n_seq; ++i) {
 			mm_bseq1_t *t = &s->seq[i];
-			mm_sketch(0, t->seq, t->l_seq, p->mi->w, p->mi->k, t->rid, p->mi->is_hpc, &s->a);
+			if (t->l_seq > 0)
+				mm_sketch(0, t->seq, t->l_seq, p->mi->w, p->mi->k, t->rid, p->mi->flag&MM_I_HPC, &s->a);
+			else if (mm_verbose >= 2)
+				fprintf(stderr, "[WARNING] the length database sequence '%s' is 0\n", t->name);
 			free(t->seq); free(t->name);
 		}
 		free(s->seq); s->seq = 0;
@@ -302,16 +308,15 @@ static void *worker_pipeline(void *shared, int step, void *in)
     return 0;
 }
 
-mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int is_hpc, int mini_batch_size, int n_threads, uint64_t batch_size, int keep_name)
+mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int flag, int mini_batch_size, int n_threads, uint64_t batch_size)
 {
 	pipeline_t pl;
 	if (fp == 0 || mm_bseq_eof(fp)) return 0;
 	memset(&pl, 0, sizeof(pipeline_t));
 	pl.mini_batch_size = mini_batch_size < batch_size? mini_batch_size : batch_size;
-	pl.keep_name = keep_name;
 	pl.batch_size = batch_size;
 	pl.fp = fp;
-	pl.mi = mm_idx_init(w, k, b, is_hpc);
+	pl.mi = mm_idx_init(w, k, b, flag);
 
 	kt_pipeline(n_threads < 3? n_threads : 3, worker_pipeline, &pl, 3);
 	if (mm_verbose >= 3)
@@ -324,13 +329,13 @@ mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int is_hpc, int mi
 	return pl.mi;
 }
 
-mm_idx_t *mm_idx_build(const char *fn, int w, int k, int is_hpc, int n_threads) // a simpler interface
+mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads) // a simpler interface
 {
 	mm_bseq_file_t *fp;
 	mm_idx_t *mi;
 	fp = mm_bseq_open(fn);
 	if (fp == 0) return 0;
-	mi = mm_idx_gen(fp, w, k, 14, is_hpc, 1<<18, n_threads, UINT64_MAX, 1);
+	mi = mm_idx_gen(fp, w, k, 14, flag, 1<<18, n_threads, UINT64_MAX);
 	mm_bseq_close(fp);
 	return mi;
 }
@@ -345,7 +350,7 @@ void mm_idx_dump(FILE *fp, const mm_idx_t *mi)
 	uint32_t x[5];
 	int i;
 
-	x[0] = mi->w, x[1] = mi->k, x[2] = mi->b, x[3] = mi->n_seq, x[4] = mi->is_hpc;
+	x[0] = mi->w, x[1] = mi->k, x[2] = mi->b, x[3] = mi->n_seq, x[4] = mi->flag;
 	fwrite(MM_IDX_MAGIC, 1, 4, fp);
 	fwrite(x, 4, 5, fp);
 	for (i = 0; i < mi->n_seq; ++i) {
@@ -372,7 +377,8 @@ void mm_idx_dump(FILE *fp, const mm_idx_t *mi)
 			fwrite(x, 8, 2, fp);
 		}
 	}
-	fwrite(mi->S, 4, (sum_len + 7) / 8, fp);
+	if (!(mi->flag & MM_I_NO_SEQ))
+		fwrite(mi->S, 4, (sum_len + 7) / 8, fp);
 	fflush(fp);
 }
 
@@ -422,8 +428,10 @@ mm_idx_t *mm_idx_load(FILE *fp)
 			kh_val(h, k) = x[1];
 		}
 	}
-	mi->S = (uint32_t*)malloc((sum_len + 7) / 8 * 4);
-	fread(mi->S, 4, (sum_len + 7) / 8, fp);
+	if (!(mi->flag & MM_I_NO_SEQ)) {
+		mi->S = (uint32_t*)malloc((sum_len + 7) / 8 * 4);
+		fread(mi->S, 4, (sum_len + 7) / 8, fp);
+	}
 	return mi;
 }
 
@@ -477,10 +485,10 @@ mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads)
 	mm_idx_t *mi;
 	if (r->is_idx) {
 		mi = mm_idx_load(r->fp.idx);
-		if (mi && mm_verbose >= 2 && (mi->k != r->opt.k || mi->w != r->opt.w || mi->is_hpc != r->opt.is_hpc))
+		if (mi && mm_verbose >= 2 && (mi->k != r->opt.k || mi->w != r->opt.w || (mi->flag&MM_I_HPC) != (r->opt.flag&MM_I_HPC)))
 			fprintf(stderr, "[WARNING]\033[1;31m Indexing parameters (-k, -w or -H) overridden by parameters used in the prebuilt index.\033[0m\n");
 	} else
-		mi = mm_idx_gen(r->fp.seq, r->opt.w, r->opt.k, r->opt.bucket_bits, r->opt.is_hpc, r->opt.mini_batch_size, n_threads, r->opt.batch_size, 1);
+		mi = mm_idx_gen(r->fp.seq, r->opt.w, r->opt.k, r->opt.bucket_bits, r->opt.flag, r->opt.mini_batch_size, n_threads, r->opt.batch_size);
 	if (mi) {
 		if (r->fp_out) mm_idx_dump(r->fp_out, mi);
 		++r->n_parts;

kalloc.c

@@ -189,6 +189,10 @@ void km_stat(const void *_km, km_stat_t *s)
 			panic("[km_stat] The end of a free block enters another free block.");
 		if (p->ptr == km->loop_head) break;
 	}
-	for (p = km->core_head; p != NULL; p = p->ptr)
-		++s->n_cores, s->capacity += p->size * sizeof(header_t);
+	for (p = km->core_head; p != NULL; p = p->ptr) {
+		size_t size = p->size * sizeof(header_t);
+		++s->n_cores;
+		s->capacity += size;
+		s->largest = s->largest > size? s->largest : size;
+	}
 }

kalloc.h

@@ -8,7 +8,7 @@ extern "C" {
 #endif
 
 typedef struct {
-	size_t capacity, available, n_blocks, n_cores;
+	size_t capacity, available, n_blocks, n_cores, largest;
 } km_stat_t;
 
 void *kmalloc(void *km, size_t size);

ksw2.h

@@ -14,6 +14,7 @@
 #define KSW_EZ_REV_CIGAR   0x80 // reverse CIGAR in the output
 #define KSW_EZ_SPLICE_FOR  0x100
 #define KSW_EZ_SPLICE_REV  0x200
+#define KSW_EZ_SPLICE_FLANK 0x400
 
 #ifdef __cplusplus
 extern "C" {

ksw2_exts2_sse.c

@@ -111,19 +111,22 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 
 	// set the donor and acceptor arrays. TODO: this assumes 0/1/2/3 encoding!
 	if (flag & (KSW_EZ_SPLICE_FOR|KSW_EZ_SPLICE_REV)) {
+		int semi_cost = flag&KSW_EZ_SPLICE_FLANK? -noncan/2 : 0; // GTr or yAG is worth 0.5 bit; see PMID:18688272
 		memset(donor, -noncan, tlen_ * 16);
-		for (t = 0; t < tlen - 2; ++t) {
-			int is_can = 0; // is a canonical site
-			if ((flag & KSW_EZ_SPLICE_FOR) && target[t+1] == 2 && target[t+2] == 3) is_can = 1;
-			if ((flag & KSW_EZ_SPLICE_REV) && target[t+1] == 1 && target[t+2] == 3) is_can = 1;
-			if (is_can) ((int8_t*)donor)[t] = 0;
+		for (t = 0; t < tlen - 4; ++t) {
+			int can_type = 0; // type of canonical site: 0=none, 1=GT/AG only, 2=GTr/yAG
+			if ((flag & KSW_EZ_SPLICE_FOR) && target[t+1] == 2 && target[t+2] == 3) can_type = 1; // GTr...
+			if ((flag & KSW_EZ_SPLICE_REV) && target[t+1] == 1 && target[t+2] == 3) can_type = 1; // CTr...
+			if (can_type && (target[t+3] == 0 || target[t+3] == 2)) can_type = 2;
+			if (can_type) ((int8_t*)donor)[t] = can_type == 2? 0 : semi_cost;
 		}
 		memset(acceptor, -noncan, tlen_ * 16);
 		for (t = 2; t < tlen; ++t) {
-			int is_can = 0;
-			if ((flag & KSW_EZ_SPLICE_FOR) && target[t-1] == 0 && target[t] == 2) is_can = 1;
-			if ((flag & KSW_EZ_SPLICE_REV) && target[t-1] == 0 && target[t] == 1) is_can = 1;
-			if (is_can) ((int8_t*)acceptor)[t] = 0;
+			int can_type = 0;
+			if ((flag & KSW_EZ_SPLICE_FOR) && target[t-1] == 0 && target[t] == 2) can_type = 1; // ...yAG
+			if ((flag & KSW_EZ_SPLICE_REV) && target[t-1] == 0 && target[t] == 1) can_type = 1; // ...yAC
+			if (can_type && (target[t-2] == 1 || target[t-2] == 3)) can_type = 2;
+			if (can_type) ((int8_t*)acceptor)[t] = can_type == 2? 0 : semi_cost;
 		}
 	}
 

main.c

@@ -6,7 +6,7 @@
 #include "mmpriv.h"
 #include "getopt.h"
 
-#define MM_VERSION "2.3-r531"
+#define MM_VERSION "2.6-r623"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -34,13 +34,17 @@ static struct option long_options[] = {
 	{ "min-dp-len",     required_argument, 0, 0 },
 	{ "print-aln-seq",  no_argument,       0, 0 },
 	{ "splice",         no_argument,       0, 0 },
-	{ "cost-non-gt-ag", required_argument, 0, 0 },
+	{ "cost-non-gt-ag", required_argument, 0, 'C' },
 	{ "no-long-join",   no_argument,       0, 0 },
 	{ "sr",             no_argument,       0, 0 },
 	{ "frag",           optional_argument, 0, 0 },
 	{ "secondary",      optional_argument, 0, 0 },
 	{ "cs",             optional_argument, 0, 0 },
 	{ "end-bonus",      required_argument, 0, 0 },
+	{ "no-pairing",     no_argument,       0, 0 },
+	{ "splice-flank",   optional_argument, 0, 0 },
+	{ "idx-no-seq",     no_argument,       0, 0 },
+	{ "end-seed-pen",   required_argument, 0, 0 },   // 21
 	{ "help",           no_argument,       0, 'h' },
 	{ "max-intron-len", required_argument, 0, 'G' },
 	{ "version",        no_argument,       0, 'V' },
@@ -65,7 +69,7 @@ static inline int64_t mm_parse_num(const char *str)
 
 int main(int argc, char *argv[])
 {
-	const char *opt_str = "2aSw: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:i:L";
+	const char *opt_str = "2aSw: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:";
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
 	int i, c, n_threads = 3, long_idx;
@@ -92,7 +96,7 @@ int main(int argc, char *argv[])
 	while ((c = getopt_long(argc, argv, opt_str, long_options, &long_idx)) >= 0) {
 		if (c == 'w') ipt.w = atoi(optarg);
 		else if (c == 'k') ipt.k = atoi(optarg);
-		else if (c == 'H') ipt.is_hpc = 1;
+		else if (c == 'H') ipt.flag |= MM_I_HPC;
 		else if (c == 'd') fnw = optarg; // the above are indexing related options, except -I
 		else if (c == 'r') opt.bw = (int)mm_parse_num(optarg);
 		else if (c == 't') n_threads = atoi(optarg);
@@ -100,7 +104,6 @@ int main(int argc, char *argv[])
 		else if (c == 'g') opt.max_gap = (int)mm_parse_num(optarg);
 		else if (c == 'G') mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(optarg));
 		else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(optarg);
-		else if (c == 'i') opt.min_iden = atof(optarg);
 		else if (c == 'N') opt.best_n = atoi(optarg);
 		else if (c == 'p') opt.pri_ratio = atof(optarg);
 		else if (c == 'M') opt.mask_level = atof(optarg);
@@ -108,6 +111,7 @@ int main(int argc, char *argv[])
 		else if (c == 'X') opt.flag |= MM_F_AVA | MM_F_NO_SELF;
 		else if (c == 'a') opt.flag |= MM_F_OUT_SAM | MM_F_CIGAR;
 		else if (c == 'Q') opt.flag |= MM_F_NO_QUAL;
+		else if (c == 'Y') opt.flag |= MM_F_SOFTCLIP;
 		else if (c == 'L') opt.flag |= MM_F_LONG_CIGAR;
 		else if (c == 'T') opt.sdust_thres = atoi(optarg);
 		else if (c == 'n') opt.min_cnt = atoi(optarg);
@@ -116,6 +120,7 @@ int main(int argc, char *argv[])
 		else if (c == 'B') opt.b = atoi(optarg);
 		else if (c == 'z') opt.zdrop = atoi(optarg);
 		else if (c == 's') opt.min_dp_max = atoi(optarg);
+		else if (c == 'C') opt.noncan = atoi(optarg);
 		else if (c == 'I') ipt.batch_size = mm_parse_num(optarg);
 		else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(optarg);
 		else if (c == 'R') rg = optarg;
@@ -131,10 +136,12 @@ int main(int argc, char *argv[])
 		else if (c == 0 && long_idx == 8) opt.min_ksw_len = atoi(optarg); // --min-dp-len
 		else if (c == 0 && long_idx == 9) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ; // --print-aln-seq
 		else if (c == 0 && long_idx ==10) opt.flag |= MM_F_SPLICE; // --splice
-		else if (c == 0 && long_idx ==11) opt.noncan = atoi(optarg); // --cost-non-gt-ag
 		else if (c == 0 && long_idx ==12) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
 		else if (c == 0 && long_idx ==13) opt.flag |= MM_F_SR; // --sr
 		else if (c == 0 && long_idx ==17) opt.end_bonus = atoi(optarg); // --end-bonus
+		else if (c == 0 && long_idx ==18) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
+		else if (c == 0 && long_idx ==20) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
+		else if (c == 0 && long_idx ==21) opt.anchor_ext_shift = atoi(optarg); // --end-seed-pen
 		else if (c == 0 && long_idx == 14) { // --frag
 			if (optarg == 0 || strcmp(optarg, "yes") == 0 || strcmp(optarg, "y") == 0)
 				opt.flag |= MM_F_FRAG_MODE;
@@ -154,6 +161,10 @@ int main(int argc, char *argv[])
 			} else if (mm_verbose >= 2) {
 				fprintf(stderr, "[WARNING]\033[1;31m --cs only takes 'short' or 'long'. Invalid values are assumed to be 'short'.\033[0m\n");
 			}
+		} else if (c == 0 && long_idx == 19) { // --splice-flank
+			if (optarg == 0 || strcmp(optarg, "yes") == 0 || strcmp(optarg, "y") == 0)
+				opt.flag |= MM_F_SPLICE_FLANK;
+			else opt.flag &= ~MM_F_SPLICE_FLANK;
 		} else if (c == 'S') {
 			opt.flag |= MM_F_OUT_CS | MM_F_CIGAR | MM_F_OUT_CS_LONG;
 			if (mm_verbose >= 2)
@@ -189,6 +200,8 @@ int main(int argc, char *argv[])
 		fprintf(stderr, "[ERROR]\033[1;31m --splice and --frag should not be specified at the same time.\033[0m\n");
 		return 1;
 	}
+	if (!fnw && !(opt.flag&MM_F_CIGAR))
+		ipt.flag |= MM_I_NO_SEQ;
 
 	if (argc == optind || fp_help == stdout) {
 		fprintf(fp_help, "Usage: minimap2 [options] <target.fa>|<target.idx> [query.fa] [...]\n");
@@ -219,7 +232,6 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -z INT       Z-drop score [%d]\n", opt.zdrop);
 		fprintf(fp_help, "    -s INT       minimal peak DP alignment score [%d]\n", opt.min_dp_max);
 		fprintf(fp_help, "    -u CHAR      how to find GT-AG. f:transcript strand, b:both strands, n:don't match GT-AG [n]\n");
-		fprintf(fp_help, "    -i FLOAT     min identity (mapQ reduced to 0 if below) [0]\n");
 		fprintf(fp_help, "  Input/Output:\n");
 		fprintf(fp_help, "    -a           output in the SAM format (PAF by default)\n");
 		fprintf(fp_help, "    -Q           don't output base quality in SAM\n");
@@ -227,6 +239,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -R STR       SAM read group line in a format like '@RG\\tID:foo\\tSM:bar' []\n");
 		fprintf(fp_help, "    -c           output CIGAR in PAF\n");
 		fprintf(fp_help, "    --cs[=STR]   output the cs tag; STR is 'short' (if absent) or 'long' [none]\n");
+		fprintf(fp_help, "    -Y           use soft clipping for supplementary alignments\n");
 		fprintf(fp_help, "    -t INT       number of threads [%d]\n", n_threads);
 		fprintf(fp_help, "    -K NUM       minibatch size for mapping [500M]\n");
 //		fprintf(fp_help, "    -v INT       verbose level [%d]\n", mm_verbose);
@@ -252,11 +265,18 @@ int main(int argc, char *argv[])
 	}
 	if (!idx_rdr->is_idx && fnw == 0 && argc - optind < 2) {
 		fprintf(stderr, "[ERROR] missing input: please specify a query file to map or option -d to keep the index\n");
+		mm_idx_reader_close(idx_rdr);
 		return 1;
 	}
 	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) {
+		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);
+			mm_idx_reader_close(idx_rdr);
+			return 1;
+		}
 		if ((opt.flag & MM_F_OUT_SAM) && idx_rdr->n_parts == 1) {
 			if (mm_idx_reader_eof(idx_rdr)) {
 				mm_write_sam_hdr(mi, rg, MM_VERSION, argc, argv);

map.c

@@ -36,6 +36,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->end_bonus = -1;
 	opt->min_dp_max = opt->min_chain_score * opt->a;
 	opt->min_ksw_len = 200;
+	opt->anchor_ext_len = 20, opt->anchor_ext_shift = 6;
 	opt->mini_batch_size = 500000000;
 
 	opt->pe_ori = 0; // FF
@@ -64,29 +65,29 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		mm_idxopt_init(io);
 		mm_mapopt_init(mo);
 	} else if (strcmp(preset, "ava-ont") == 0) {
-		io->is_hpc = 0, io->k = 15, io->w = 5;
+		io->flag = 0, io->k = 15, io->w = 5;
 		mo->flag |= MM_F_AVA | MM_F_NO_SELF;
 		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
 	} else if (strcmp(preset, "ava-pb") == 0) {
-		io->is_hpc = 1, io->k = 19, io->w = 5;
+		io->flag |= MM_I_HPC, io->k = 19, io->w = 5;
 		mo->flag |= MM_F_AVA | MM_F_NO_SELF;
 		mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_gap = 10000, mo->max_chain_skip = 25;
 	} else if (strcmp(preset, "map10k") == 0 || strcmp(preset, "map-pb") == 0) {
-		io->is_hpc = 1, io->k = 19;
+		io->flag |= MM_I_HPC, io->k = 19;
 	} else if (strcmp(preset, "map-ont") == 0) {
-		io->is_hpc = 0, io->k = 15;
+		io->flag = 0, io->k = 15;
 	} else if (strcmp(preset, "asm5") == 0) {
-		io->is_hpc = 0, io->k = 19, io->w = 19;
+		io->flag = 0, io->k = 19, io->w = 19;
 		mo->a = 1, mo->b = 19, mo->q = 39, mo->q2 = 81, mo->e = 3, mo->e2 = 1, mo->zdrop = 200;
 		mo->min_dp_max = 200;
 		mo->best_n = 50;
 	} else if (strcmp(preset, "asm10") == 0) {
-		io->is_hpc = 0, io->k = 19, io->w = 19;
+		io->flag = 0, io->k = 19, io->w = 19;
 		mo->a = 1, mo->b = 9, mo->q = 16, mo->q2 = 41, mo->e = 2, mo->e2 = 1, mo->zdrop = 200;
 		mo->min_dp_max = 200;
 		mo->best_n = 50;
 	} else if (strcmp(preset, "short") == 0 || strcmp(preset, "sr") == 0) {
-		io->is_hpc = 0, io->k = 21, io->w = 11;
+		io->flag = 0, io->k = 21, io->w = 11;
 		mo->flag |= MM_F_SR | MM_F_FRAG_MODE | MM_F_NO_PRINT_2ND | MM_F_2_IO_THREADS;
 		mo->pe_ori = 0<<1|1; // FR
 		mo->a = 2, mo->b = 8, mo->q = 12, mo->e = 2, mo->q2 = 24, mo->e2 = 1;
@@ -104,11 +105,11 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		mo->max_occ = 5000;
 		mo->mini_batch_size = 50000000;
 	} else if (strcmp(preset, "splice") == 0 || strcmp(preset, "cdna") == 0) {
-		io->is_hpc = 0, io->k = 15, io->w = 5;
-		mo->flag |= MM_F_SPLICE | MM_F_SPLICE_FOR | MM_F_SPLICE_REV;
+		io->flag = 0, io->k = 15, io->w = 5;
+		mo->flag |= MM_F_SPLICE | MM_F_SPLICE_FOR | MM_F_SPLICE_REV | MM_F_SPLICE_FLANK;
 		mo->max_gap = 2000, mo->max_gap_ref = mo->bw = 200000;
 		mo->a = 1, mo->b = 2, mo->q = 2, mo->e = 1, mo->q2 = 32, mo->e2 = 0;
-		mo->noncan = 5;
+		mo->noncan = 9;
 		mo->zdrop = 200;
 	} else return -1;
 	return 0;
@@ -122,8 +123,6 @@ typedef struct {
 } mm_match_t;
 
 struct mm_tbuf_s {
-	sdust_buf_t *sdb;
-	mm128_v mini;
 	void *km;
 };
 
@@ -132,24 +131,23 @@ mm_tbuf_t *mm_tbuf_init(void)
 	mm_tbuf_t *b;
 	b = (mm_tbuf_t*)calloc(1, sizeof(mm_tbuf_t));
 	if (!(mm_dbg_flag & 1)) b->km = km_init();
-	b->sdb = sdust_buf_init(b->km);
 	return b;
 }
 
 void mm_tbuf_destroy(mm_tbuf_t *b)
 {
 	if (b == 0) return;
-	kfree(b->km, b->mini.a);
-	sdust_buf_destroy(b->sdb);
 	km_destroy(b->km);
 	free(b);
 }
 
-static int mm_dust_minier(int n, mm128_t *a, int l_seq, const char *seq, int sdust_thres, sdust_buf_t *sdb)
+static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *seq, int sdust_thres)
 {
 	int n_dreg, j, k, u = 0;
 	const uint64_t *dreg;
-	if (sdust_thres <= 0 || sdb == 0) return n;
+	sdust_buf_t *sdb;
+	if (sdust_thres <= 0) return n;
+	sdb = sdust_buf_init(km);
 	dreg = sdust_core((const uint8_t*)seq, l_seq, sdust_thres, 64, &n_dreg, sdb);
 	for (j = k = 0; j < n; ++j) { // squeeze out minimizers that significantly overlap with LCRs
 		int32_t qpos = (uint32_t)a[j].y>>1, span = a[j].x&0xff;
@@ -163,45 +161,49 @@ static int mm_dust_minier(int n, mm128_t *a, int l_seq, const char *seq, int sdu
 				l += ee - ss;
 			}
 			if (l <= span>>1) a[k++] = a[j]; // keep the minimizer if less than half of it falls in masked region
-		}
+		} else a[k++] = a[j];
 	}
+	sdust_buf_destroy(sdb);
 	return k; // the new size
 }
 
-static void collect_minimizers(const mm_mapopt_t *opt, const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, mm_tbuf_t *b)
+static void collect_minimizers(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, mm128_v *mv)
 {
 	int i, j, n, sum = 0;
-	b->mini.n = 0;
+	mv->n = 0;
 	for (i = n = 0; i < n_segs; ++i) {
-		mm_sketch(b->km, seqs[i], qlens[i], mi->w, mi->k, i, mi->is_hpc, &b->mini);
-		for (j = n; j < b->mini.n; ++j)
-			b->mini.a[j].y += sum << 1;
+		mm_sketch(km, seqs[i], qlens[i], mi->w, mi->k, i, mi->flag&MM_I_HPC, mv);
+		for (j = n; j < mv->n; ++j)
+			mv->a[j].y += sum << 1;
 		if (opt->sdust_thres > 0) // mask low-complexity minimizers
-			b->mini.n = n + mm_dust_minier(b->mini.n - n, b->mini.a + n, qlens[i], seqs[i], opt->sdust_thres, b->sdb);
-		sum += qlens[i], n = b->mini.n;
+			mv->n = n + mm_dust_minier(km, mv->n - n, mv->a + n, qlens[i], seqs[i], opt->sdust_thres);
+		sum += qlens[i], n = mv->n;
 	}
 }
 
-static mm128_t *collect_seed_hits(const mm_mapopt_t *opt, int max_occ, const mm_idx_t *mi, const char *qname, int qlen, int64_t *n_a, int *rep_len, mm_tbuf_t *b)
+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 rep_st = 0, rep_en = 0, i;
 	mm_match_t *m;
 	mm128_t *a;
 
-	m = (mm_match_t*)kmalloc(b->km, b->mini.n * sizeof(mm_match_t));
-	for (i = 0; i < b->mini.n; ++i) {
+	*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) {
 		int t;
-		mm128_t *p = &b->mini.a[i];
+		mm128_t *p = &mv->a[i];
 		m[i].qpos = (uint32_t)p->y;
 		m[i].cr = mm_idx_get(mi, p->x>>8, &t);
 		m[i].n = t;
 		m[i].seg_id = p->y >> 32;
 	}
-	for (i = 0, *n_a = 0; i < b->mini.n; ++i) // find the length of a[]
+	for (i = 0, *n_a = 0; i < mv->n; ++i) // find the length of a[]
 		if (m[i].n < max_occ) *n_a += m[i].n;
-	a = (mm128_t*)kmalloc(b->km, *n_a * sizeof(mm128_t));
-	for (i = *rep_len = 0, *n_a = 0; i < b->mini.n; ++i) {
-		mm128_t *p = &b->mini.a[i];
+	a = (mm128_t*)kmalloc(km, *n_a * sizeof(mm128_t));
+	for (i = *rep_len = 0, *n_a = 0; i < mv->n; ++i) {
+		mm128_t *p = &mv->a[i];
 		mm_match_t *q = &m[i];
 		const uint64_t *r = q->cr;
 		int k, q_span = p->x & 0xff, is_tandem = 0;
@@ -213,8 +215,9 @@ static mm128_t *collect_seed_hits(const mm_mapopt_t *opt, int max_occ, const mm_
 			} else rep_en = en;
 			continue;
 		}
-		if (i > 0 && p->x>>8 == b->mini.a[i - 1].x>>8) is_tandem = 1;
-		if (i < b->mini.n - 1 && p->x>>8 == b->mini.a[i + 1].x>>8) is_tandem = 1;
+		(*mini_pos)[(*n_mini_pos)++] = (uint64_t)q_span<<32 | q->qpos>>1;
+		if (i > 0 && p->x>>8 == mv->a[i - 1].x>>8) is_tandem = 1;
+		if (i < mv->n - 1 && p->x>>8 == mv->a[i + 1].x>>8) is_tandem = 1;
 		for (k = 0; k < q->n; ++k) {
 			int32_t rpos = (uint32_t)r[k] >> 1;
 			mm128_t *p;
@@ -240,7 +243,7 @@ static mm128_t *collect_seed_hits(const mm_mapopt_t *opt, int max_occ, const mm_
 		}
 	}
 	*rep_len += rep_en - rep_st;
-	kfree(b->km, m);
+	kfree(km, m);
 	return a;
 }
 
@@ -258,7 +261,7 @@ static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_i
 static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *km, int qlen, const char *seq, const char *qual, int *n_regs, mm_reg1_t *regs, mm128_t *a)
 {
 	if (!(opt->flag & MM_F_CIGAR)) return regs;
-	regs = mm_align_skeleton(km, opt, mi, qlen, seq, qual, n_regs, regs, a); // this calls mm_filter_regs()
+	regs = mm_align_skeleton(km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs()
 	if (!(opt->flag & MM_F_AVA)) {
 		mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b);
 		mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs);
@@ -269,13 +272,15 @@ 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, const char **quals, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
 {
-	int i, j, rep_len, qlen_sum, n_regs0;
-	int max_chain_gap_qry, max_chain_gap_ref, is_splice = !!(opt->flag & MM_F_SPLICE);
+	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;
 	int64_t n_a;
-	uint64_t *u;
+	uint64_t *u, *mini_pos;
 	mm128_t *a;
+	mm128_v mv = {0,0,0};
 	mm_reg1_t *regs0;
+	km_stat_t kmst;
 
 	for (i = 0, qlen_sum = 0; i < n_segs; ++i)
 		qlen_sum += qlens[i], n_regs[i] = 0, regs[i] = 0;
@@ -286,8 +291,8 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	hash ^= __ac_Wang_hash(qlen_sum) + __ac_Wang_hash(opt->seed);
 	hash  = __ac_Wang_hash(hash);
 
-	collect_minimizers(opt, mi, n_segs, qlens, seqs, b);
-	a = collect_seed_hits(opt, opt->mid_occ, mi, qname, qlen_sum, &n_a, &rep_len, b);
+	collect_minimizers(b->km, opt, mi, n_segs, qlens, seqs, &mv);
+	a = collect_seed_hits(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 	radix_sort_128x(a, a + n_a);
 
 	if (mm_dbg_flag & MM_DBG_PRINT_SEED) {
@@ -298,7 +303,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 	}
 
 	// set max chaining gap on the query and the reference sequence
-	if (opt->flag & MM_F_SR)
+	if (is_sr)
 		max_chain_gap_qry = qlen_sum > opt->max_gap? qlen_sum : opt->max_gap;
 	else max_chain_gap_qry = opt->max_gap;
 	if (opt->max_gap_ref > 0) {
@@ -327,7 +332,8 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		if (rechain) { // redo chaining with a higher max_occ threshold
 			kfree(b->km, a);
 			kfree(b->km, u);
-			a = collect_seed_hits(opt, opt->max_occ, mi, qname, qlen_sum, &n_a, &rep_len, b);
+			kfree(b->km, mini_pos);
+			a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 			radix_sort_128x(a, a + n_a);
 			a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 		}
@@ -342,10 +348,11 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 						i == regs0[j].as? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x));
 
 	chain_post(opt, max_chain_gap_ref, mi, b->km, qlen_sum, n_segs, qlens, &n_regs0, regs0, a);
+	if (!is_sr) mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos);
 
 	if (n_segs == 1) { // uni-segment
 		regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], quals? quals[0] : 0, &n_regs0, regs0, a);
-		mm_set_mapq(n_regs0, regs0, opt->min_chain_score, opt->a, rep_len);
+		mm_set_mapq(n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr);
 		n_regs[0] = n_regs0, regs[0] = regs0;
 	} else { // multi-segment
 		mm_seg_t *seg;
@@ -354,18 +361,28 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		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); // update mm_reg1_t::parent
 			regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], quals? quals[i] : 0, &n_regs[i], regs[i], seg[i].a);
-			mm_set_mapq(n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len);
+			mm_set_mapq(n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr);
 		}
 		mm_seg_free(b->km, n_segs, seg);
 		if (n_segs == 2 && opt->pe_ori >= 0 && (opt->flag&MM_F_CIGAR))
 			mm_pair(b->km, max_chain_gap_ref, opt->pe_bonus, opt->a * 2 + opt->b, opt->a, qlens, n_regs, regs); // pairing
 	}
-	if (opt->min_iden > 0.0f)
-		for (i = 0; i < n_segs; ++i)
-			mm_filter_by_identity(b->km, n_regs[i], regs[i], opt->min_iden, qlens[i], quals[i]);
 
+	kfree(b->km, mv.a);
 	kfree(b->km, a);
 	kfree(b->km, u);
+	kfree(b->km, mini_pos);
+
+	if (b->km) {
+		km_stat(b->km, &kmst);
+		if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
+			fprintf(stderr, "QM\t%s\t%d\tcap=%ld,nCore=%ld,largest=%ld\n", qname, qlen_sum, kmst.capacity, kmst.n_cores, kmst.largest);
+		assert(kmst.n_blocks == kmst.n_cores); // otherwise, there is a memory leak
+		if (kmst.largest > 1U<<28) {
+			km_destroy(b->km);
+			b->km = km_init();
+		}
+	}
 }
 
 mm_reg1_t *mm_map(const mm_idx_t *mi, int qlen, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname)
@@ -399,14 +416,13 @@ typedef struct {
 static void worker_for(void *_data, long i, int tid) // kt_for() callback
 {
     step_t *s = (step_t*)_data;
-	int *qlens, j, off = s->seg_off[i], pe_ori = s->p->opt->pe_ori, is_sr = !!(s->p->opt->flag & MM_F_SR);
-	const char **qseqs, **quals = 0;
+	int qlens[MM_MAX_SEG], j, off = s->seg_off[i], pe_ori = s->p->opt->pe_ori, is_sr = !!(s->p->opt->flag & MM_F_SR);
+	const char *qseqs[MM_MAX_SEG], *quals[MM_MAX_SEG];
 	mm_tbuf_t *b = s->buf[tid];
+	assert(s->n_seg[i] <= MM_MAX_SEG);
+	memset(quals, 0, sizeof(char*) * MM_MAX_SEG);
 	if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
-		fprintf(stderr, "QR\t%s\t%d\n", s->seq[off].name, tid);
-	qlens = (int*)kmalloc(b->km, s->n_seg[i] * sizeof(int));
-	qseqs = (const char**)kmalloc(b->km, s->n_seg[i] * sizeof(const char**));
-	quals = (const char**)kmalloc(b->km, s->n_seg[i] * sizeof(const char**));
+		fprintf(stderr, "QR\t%s\t%d\t%d\n", s->seq[off].name, tid, s->seq[off].l_seq);
 	for (j = 0; j < s->n_seg[i]; ++j) {
 		if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1))))
 			mm_revcomp_bseq(&s->seq[off + j]);
@@ -414,7 +430,12 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
 		qseqs[j] = s->seq[off + j].seq;
 		quals[j] = is_sr? s->seq[off + j].qual : 0;
 	}
-	mm_map_frag(s->p->mi, s->n_seg[i], qlens, qseqs, quals, &s->n_reg[off], &s->reg[off], b, s->p->opt, s->seq[off].name);
+	if (s->p->opt->flag & MM_F_INDEPEND_SEG) {
+		for (j = 0; j < s->n_seg[i]; ++j)
+			mm_map_frag(s->p->mi, 1, &qlens[j], &qseqs[j], &quals[j], &s->n_reg[off+j], &s->reg[off+j], b, s->p->opt, s->seq[off+j].name);
+	} else {
+		mm_map_frag(s->p->mi, s->n_seg[i], qlens, qseqs, quals, &s->n_reg[off], &s->reg[off], b, s->p->opt, s->seq[off].name);
+	}
 	for (j = 0; j < s->n_seg[i]; ++j) // flip the query strand and coordinate to the original read strand
 		if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1)))) {
 			int k, t;
@@ -427,9 +448,6 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
 				r->rev = !r->rev;
 			}
 		}
-	kfree(b->km, qlens);
-	kfree(b->km, qseqs);
-	kfree(b->km, quals);
 }
 
 static void *worker_pipeline(void *shared, int step, void *in)

minimap.h

@@ -22,6 +22,13 @@
 #define MM_F_NO_PRINT_2ND  0x4000
 #define MM_F_2_IO_THREADS  0x8000
 #define MM_F_LONG_CIGAR    0x10000
+#define MM_F_INDEPEND_SEG  0x20000
+#define MM_F_SPLICE_FLANK  0x40000
+#define MM_F_SOFTCLIP      0x80000
+
+#define MM_I_HPC          0x1
+#define MM_I_NO_SEQ       0x2
+#define MM_I_NO_NAME      0x4
 
 #define MM_IDX_MAGIC   "MMI\2"
 
@@ -43,7 +50,7 @@ typedef struct {
 } mm_idx_seq_t;
 
 typedef struct {
-	int32_t b, w, k, is_hpc;
+	int32_t b, w, k, flag;
 	uint32_t n_seq;            // number of reference sequences
 	mm_idx_seq_t *seq;         // sequence name, length and offset
 	uint32_t *S;               // 4-bit packed sequence
@@ -57,29 +64,29 @@ typedef struct {
 	int32_t dp_score, dp_max, dp_max2;  // DP score; score of the max-scoring segment; score of the best alternate mappings
 	uint32_t n_ambi:30, trans_strand:2; // number of ambiguous bases; transcript strand: 0 for unknown, 1 for +, 2 for -
 	uint32_t n_cigar;                   // number of cigar operations in cigar[]
-	float n_diff2;
-	uint32_t blen2;
 	uint32_t cigar[];
 } mm_extra_t;
 
 typedef struct {
-	int32_t id;                     // ID for internal uses (see also parent below)
-	uint32_t cnt:30, rev:1, seg_split:1; // number of minimizers; if on the reverse strand
-	uint32_t rid:31, inv:1;         // reference index; if this is an alignment from inversion rescue
-	int32_t score;                  // DP alignment score
-	int32_t qs, qe, rs, re;         // query start and end; reference start and end
-	int32_t parent, subsc;          // parent==id if primary; best alternate mapping score
-	int32_t as;                     // offset in the a[] array (for internal uses only)
-	int32_t mlen, blen;             // seeded exact match length; seeded alignment block length
-	uint32_t mapq:8, split:2, n_sub:22; // mapQ; split pattern; number of suboptimal mappings
-	uint32_t sam_pri:1, proper_frag:1, iden_flt:1, pe_thru:1, dummy:29;
+	int32_t id;             // ID for internal uses (see also parent below)
+	int32_t cnt;            // number of minimizers; if on the reverse strand
+	int32_t rid;            // reference index; if this is an alignment from inversion rescue
+	int32_t score;          // DP alignment score
+	int32_t qs, qe, rs, re; // query start and end; reference start and end
+	int32_t parent, subsc;  // parent==id if primary; best alternate mapping score
+	int32_t as;             // offset in the a[] array (for internal uses only)
+	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, dummy:16;
 	uint32_t hash;
+	float div;
 	mm_extra_t *p;
 } mm_reg1_t;
 
 // indexing and mapping options
 typedef struct {
-	short k, w, is_hpc, bucket_bits;
+	short k, w, flag, bucket_bits;
 	int mini_batch_size;
 	uint64_t batch_size;
 } mm_idxopt_t;
@@ -99,7 +106,6 @@ typedef struct {
 	float mask_level;
 	float pri_ratio;
 	int best_n;      // top best_n chains are subjected to DP alignment
-	float min_iden;
 
 	int max_join_long, max_join_short;
 	int min_join_flank_sc;
@@ -110,6 +116,7 @@ typedef struct {
 	int end_bonus;
 	int min_dp_max;  // drop an alignment if the score of the max scoring segment is below this threshold
 	int min_ksw_len;
+	int anchor_ext_len, anchor_ext_shift;
 
 	int pe_ori, pe_bonus;
 
@@ -269,7 +276,7 @@ int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_
 
 // deprecated APIs for backward compatibility
 void mm_mapopt_init(mm_mapopt_t *opt);
-mm_idx_t *mm_idx_build(const char *fn, int w, int k, int is_hpc, int n_threads);
+mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads);
 
 #ifdef __cplusplus
 }

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "22 October 2017" "minimap2-2.2-dirty (r531)" "Bioinformatics tools"
+.TH minimap2 1 "12 December 2017" "minimap2-2.6 (r623)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -99,6 +99,14 @@ multiple times to map it against each batch of target sequences.
 may be ending with k/K/m/M/g/G. NB: mapping quality is incorrect given a
 multi-part index.
 .TP
+.B --idx-no-seq
+Don't store target sequences in the index. It saves disk space and memory but
+the index generated with this option will not work with
+.B -a
+or
+.BR -c .
+When base-level alignment is not requested, this option is automatically applied.
+.TP
 .BI -d \ FILE
 Save the minimizer index of
 .I target.fa
@@ -219,6 +227,11 @@ costs
 .RI min{ O1 + k * E1 , O2 + k * E2 }.
 In the splice mode, the second gap penalties are not used.
 .TP
+.BI -C \ INT
+Cost for a non-canonical GT-AG splicing (effective with
+.BR --splice )
+[0]
+.TP
 .BI -z \ INT
 Break an alignment if the running score drops too quickly along the diagonal of
 the DP matrix (diagonal X-drop, or Z-drop) [400]. Increasing the value improves
@@ -239,11 +252,38 @@ both strands;
 .BR n :
 no attempt to match GT-AG [n]
 .TP
-.BI --cost-non-gt-ag \ INT
-Cost of non-canonical splicing sites [0].
-.TP
 .BI --end-bonus \ INT
-Score bonus when alignment extends to the end of the query sequence [10].
+Score bonus when alignment extends to the end of the query sequence [0].
+.TP
+.BR --splice-flank [= yes | no ]
+Assume the next base to a
+.B GT
+donor site tends to be A/G (91% in human and 92% in mouse) and the preceding
+base to a
+.B AG
+acceptor tends to be C/T [yes with
+.BR --splice ].
+This trend is evolutionarily conservative, all the way to S. cerevisiae
+(PMID:18688272). Specifying this option generally leads to higher junction
+accuracy by several percents, so it is applied by default with
+.BR --splice .
+However, the SIRV control does not honor this trend
+(only ~60%). This option reduces accuracy. If you are benchmarking minimap2
+on SIRV data, please add
+.B --splice-flank=no
+to the command line.
+.TP
+.BI --end-seed-pen \ INT
+Drop a terminal anchor if
+.IR s <log( g )+ INT ,
+where
+.I s
+is the local alignment score around the anchor and
+.I g
+the length of the terminal gap in the chain. This option is only effective
+with
+.BR --splice .
+It helps to avoid tiny terminal exons. [6]
 .SS Input/output options
 .TP 10
 .B -a
@@ -261,7 +301,7 @@ the real CIGAR in memory.
 .TP
 .BI -R \ STR
 SAM read group line in a format like
-.RB @RG\\\\tID:foo\\\\tSM:bar
+.B @RG\\\\tID:foo\\\\tSM:bar
 [].
 .TP
 .B -c
@@ -282,6 +322,9 @@ is given,
 .I short
 is assumed. [none]
 .TP
+.B -Y
+In SAM output, use soft clipping for supplementary alignments.
+.TP
 .BI --seed \ INT
 Integer seed for randomizing equally best hits. Minimap2 hashes
 .I INT
@@ -371,8 +414,8 @@ is that this preset is not using HPC minimizers.
 .B splice
 Long-read spliced alignment
 .RB ( -k15
-.B -w5 --splice -g2000 -G200k -A1 -B2 -O2,32 -E1,0 -z200 -ub --cost-non-gt-ag
-.BR 5 ).
+.B -w5 --splice -g2000 -G200k -A1 -B2 -O2,32 -E1,0 -C9 -z200 -ub
+.BR --splice-flank=yes ).
 In the splice mode, 1) long deletions are taken as introns and represented as
 the
 .RB ` N '
@@ -438,7 +481,7 @@ cb | cb | cb
 r | c | l .
 Tag	Type	Description
 _
-tp	A	Type of aln: P/primary, S/secondary and I/inversion
+tp	A	Type of aln: P/primary, S/secondary and I,i/inversion
 cm	i	Number of minimizers on the chain
 s1	i	Chaining score
 s2	i	Chaining score of the best secondary chain

misc/intron-eval.js

@@ -161,7 +161,8 @@ for (var tid in tr) {
 	if (anno[t[0]] == null) anno[t[0]] = [];
 	var s = t[4];
 	for (var i = 0; i < s.length - 1; ++i) {
-		if (s[i][1] >= s[i+1][0]) throw Error("ERROR: wrong annotation!");
+		if (s[i][1] >= s[i+1][0])
+			warn("WARNING: incorrect annotation for transcript "+tid+" ("+s[i][1]+" >= "+s[i+1][0]+")")
 		anno[t[0]].push([s[i][1], s[i+1][0]]);
 	}
 }

misc/paf2aln.js

@@ -36,20 +36,21 @@ var getopt = function(args, ostr) {
 	return optopt;
 }
 
-var c, maf_out = false, line_len = 80;
-while ((c = getopt(arguments, "ml:")) != null) {
-	if (c == 'm') maf_out = true;
-	else if (c == 'l') line_len = parseInt(getopt.arg); // TODO: not implemented yet
+var c, line_len = 80, fmt = "aln";
+while ((c = getopt(arguments, "f:l:")) != null) {
+	if (c == 'f') {
+		fmt = getopt.arg;
+		if (fmt != "aln" && fmt != "lastz-cigar" && fmt != "maf")
+			throw Error("format must be one of aln, lastz-cigar and maf");
+	} else if (c == 'l') line_len = parseInt(getopt.arg);
 }
 if (line_len == 0) line_len = 0x7fffffff;
 
 if (getopt.ind == arguments.length) {
-	print("Usage: k8 paf2aln.js [options] <with-cs.paf>");
+	print("Usage: k8 paf2aln.js [options] <in.paf>");
 	print("Options:");
-	print("  -m        MAF output (BLAST-like output by default)");
+	print("  -f STR    output format: aln (BLAST-like), maf or lastz-cigar [aln]");
 	print("  -l INT    line length in BLAST-like output [80]");
-	print("");
-	print("Note: this script only works when minimap2 is run with option '-S'");
 	exit(1);
 }
 
@@ -66,7 +67,7 @@ function padding_str(x, len, right)
 function update_aln(s_ref, s_qry, s_mid, type, seq, slen)
 {
 	var l = type == '*'? 1 : seq.length;
-	if (type == '=') {
+	if (type == '=' || type == ':') {
 		s_ref.set(seq);
 		s_qry.set(seq);
 		s_mid.set(Array(l+1).join("|"));
@@ -96,49 +97,73 @@ function print_aln(rs, qs, strand, slen, elen, s_ref, s_qry, s_mid)
 	var st, en;
 	if (strand == '+') st = qs + slen[1] + 1, en = qs + elen[1];
 	else st = qs - slen[1], en = qs - elen[1] + 1;
-	print(["Qry" + strand + ":", padding_str(st,               10, false), s_qry.toString(), padding_str(en          , 10, true)].join(" "));
+	print(["Qry" + strand + ":", padding_str(st, 10, false), s_qry.toString(), padding_str(en, 10, true)].join(" "));
 }
 
-var s_ref = new Bytes(), s_qry = new Bytes(), s_mid = new Bytes();
-var re = /([=\-\+\*])([A-Za-z]+)/g;
+var s_ref = new Bytes(), s_qry = new Bytes(), s_mid = new Bytes(); // these are used to show padded alignment
+var re_cs = /([:=\-\+\*])(\d+|[A-Za-z]+)/g;
+var re_cg = /(\d+)([MIDNSH])/g;
 
 var buf = new Bytes();
-var file = new File(arguments[getopt.ind]);
-if (maf_out) print("##maf version=1\n");
+var file = arguments[getopt.ind] == "-"? new File() : new File(arguments[getopt.ind]);
+var lineno = 0;
+if (fmt == "maf") print("##maf version=1\n");
 while (file.readline(buf) >= 0) {
 	var m, line = buf.toString();
 	var t = line.split("\t", 12);
-	if ((m = /\tcs:Z:(\S+)/.exec(line)) == null) continue;
-	var cs = m[1];
+	++lineno;
 	s_ref.length = s_qry.length = s_mid.length = 0;
 	var slen = [0, 0], elen = [0, 0];
-	if (maf_out) {
-		while ((m = re.exec(cs)) != null)
-			update_aln(s_ref, s_qry, s_mid, m[1], m[2], elen);
-		if (maf_out) {
-			var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : 0;
-			var len = t[0].length > t[5].length? t[0].length : t[5].length;
-			print("a " + score);
-			print(["s", padding_str(t[5], len, true), padding_str(t[7], 10, false), padding_str(parseInt(t[8]) - parseInt(t[7]), 10, false),
-				   "+", padding_str(t[6], 10, false), s_ref.toString()].join(" "));
-			var qs, qe, ql = parseInt(t[1]);
-			if (t[4] == '+') {
-				qs = parseInt(t[2]);
-				qe = parseInt(t[3]);
-			} else {
-				qs = ql - parseInt(t[3]);
-				qe = ql - parseInt(t[2]);
-			}
-			print(["s", padding_str(t[0], len, true), padding_str(qs, 10, false), padding_str(qe - qs, 10, false),
-				   t[4], padding_str(ql, 10, false), s_qry.toString()].join(" "));
-			print("");
+	if (fmt == "lastz-cigar") { // LASTZ-cigar output
+		var cg = (m = /\tcg:Z:(\S+)/.exec(line)) != null? m[1] : null;
+		if (cg == null) {
+			warn("WARNING: converting to LASTZ-cigar format requires the 'cg' tag, which is absent on line " + lineno);
+			continue;
 		}
-	} else {
-		line = line.replace(/\tc[sg]:Z:\S+/g, "");
+		var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? m[1] : 0;
+		var out = ['cigar:', t[0], t[2], t[3], t[4], t[5], t[7], t[8], '+', score];
+		while ((m = re_cg.exec(cg)) != null)
+			out.push(m[2], m[1]);
+		print(out.join(" "));
+	} else if (fmt == "maf") { // MAF output
+		var cs = (m = /\tcs:Z:(\S+)/.exec(line)) != null? m[1] : null;
+		if (cs == null) {
+			warn("WARNING: converting to MAF requires the 'cs' tag, which is absent on line " + lineno);
+			continue;
+		}
+		while ((m = re_cs.exec(cs)) != null) {
+			if (m[1] == ':')
+				throw Error("converting to MAF only works with 'minimap2 --cs=long'");
+			update_aln(s_ref, s_qry, s_mid, m[1], m[2], elen);
+		}
+		var score = (m = /\tAS:i:(\d+)/.exec(line)) != null? parseInt(m[1]) : 0;
+		var len = t[0].length > t[5].length? t[0].length : t[5].length;
+		print("a " + score);
+		print(["s", padding_str(t[5], len, true), padding_str(t[7], 10, false), padding_str(parseInt(t[8]) - parseInt(t[7]), 10, false),
+				"+", padding_str(t[6], 10, false), s_ref.toString()].join(" "));
+		var qs, qe, ql = parseInt(t[1]);
+		if (t[4] == '+') {
+			qs = parseInt(t[2]);
+			qe = parseInt(t[3]);
+		} else {
+			qs = ql - parseInt(t[3]);
+			qe = ql - parseInt(t[2]);
+		}
+		print(["s", padding_str(t[0], len, true), padding_str(qs, 10, false), padding_str(qe - qs, 10, false),
+				t[4], padding_str(ql, 10, false), s_qry.toString()].join(" "));
+		print("");
+	} else { // BLAST-like output
+		var cs = (m = /\tcs:Z:(\S+)/.exec(line)) != null? m[1] : null;
+		if (cs == null) {
+			warn("WARNING: converting to BLAST-like alignment requires the 'cs' tag, which is absent on line " + lineno);
+			continue;
+		}
+		line = line.replace(/\tc[sg]:Z:\S+/g, ""); // get rid of cs or cg tags
 		print('>' + line);
 		var rs = parseInt(t[7]), qs = t[4] == '+'? parseInt(t[2]) : parseInt(t[3]);
 		var n_blocks = 0;
-		while ((m = re.exec(cs)) != null) {
+		while ((m = re_cs.exec(cs)) != null) {
+			if (m[1] == ':') m[2] = Array(parseInt(m[2]) + 1).join("=");
 			var start = 0, rest = m[1] == '*'? 1 : m[2].length;
 			while (rest > 0) {
 				var l_proc;

misc/paf2diff.js

@@ -0,0 +1,188 @@
+var getopt = function(args, ostr) {
+	var oli; // option letter list index
+	if (typeof(getopt.place) == 'undefined')
+		getopt.ind = 0, getopt.arg = null, getopt.place = -1;
+	if (getopt.place == -1) { // update scanning pointer
+		if (getopt.ind >= args.length || args[getopt.ind].charAt(getopt.place = 0) != '-') {
+			getopt.place = -1;
+			return null;
+		}
+		if (getopt.place + 1 < args[getopt.ind].length && args[getopt.ind].charAt(++getopt.place) == '-') { // found "--"
+			++getopt.ind;
+			getopt.place = -1;
+			return null;
+		}
+	}
+	var optopt = args[getopt.ind].charAt(getopt.place++); // character checked for validity
+	if (optopt == ':' || (oli = ostr.indexOf(optopt)) < 0) {
+		if (optopt == '-') return null; //  if the user didn't specify '-' as an option, assume it means null.
+		if (getopt.place < 0) ++getopt.ind;
+		return '?';
+	}
+	if (oli+1 >= ostr.length || ostr.charAt(++oli) != ':') { // don't need argument
+		getopt.arg = null;
+		if (getopt.place < 0 || getopt.place >= args[getopt.ind].length) ++getopt.ind, getopt.place = -1;
+	} else { // need an argument
+		if (getopt.place >= 0 && getopt.place < args[getopt.ind].length)
+			getopt.arg = args[getopt.ind].substr(getopt.place);
+		else if (args.length <= ++getopt.ind) { // no arg
+			getopt.place = -1;
+			if (ostr.length > 0 && ostr.charAt(0) == ':') return ':';
+			return '?';
+		} else getopt.arg = args[getopt.ind]; // white space
+		getopt.place = -1;
+		++getopt.ind;
+	}
+	return optopt;
+}
+
+var re_cs = /([:=*+-])(\d+|[A-Za-z]+)/g;
+var c, min_cov_len = 10000, min_var_len = 50000, gap_thres = 50, min_mapq = 5;
+while ((c = getopt(arguments, "l:L:g:q:")) != null) {
+	if (c == 'l') min_cov_len = parseInt(getopt.arg);
+	else if (c == 'L') min_var_len = parseInt(optarg.arg);
+	else if (c == 'g') gap_thres = parseInt(optarg.arg);
+	else if (c == 'q') min_mapq = parseInt(optarg.arg);
+}
+
+if (arguments.length == getopt.ind) {
+	print("Usage: k8 paf2diff.js [options] <with-cs.paf>");
+	print("Options:");
+	print("  -l INT    min alignment length to compute coverage ["+min_cov_len+"]");
+	print("  -L INT    min alignment length to call variants ["+min_var_len+"]");
+	print("  -q INT    min mapping quality ["+min_mapq+"]");
+	print("  -g INT    short/long gap threshold (for statistics only) ["+gap_thres+"]");
+	exit(1);
+}
+
+var file = new File(arguments[getopt.ind]);
+var buf = new Bytes();
+var tot_len = 0, n_sub = [0, 0, 0], n_ins = [0, 0, 0, 0], n_del = [0, 0, 0, 0];
+
+function count_var(o)
+{
+	if (o[3] > 1) return;
+	if (o[5] == '-' && o[6] == '-') return;
+	if (o[5] == '-') { // insertion
+		var l = o[6].length;
+		if (l == 1) ++n_ins[0];
+		else if (l == 2) ++n_ins[1];
+		else if (l < gap_thres) ++n_ins[2];
+		else ++n_ins[3];
+	} else if (o[6] == '-') { // deletion
+		var l = o[5].length;
+		if (l == 1) ++n_del[0];
+		else if (l == 2) ++n_del[1];
+		else if (l < gap_thres) ++n_del[2];
+		else ++n_del[3];
+	} else {
+		++n_sub[0];
+		var s = o[5] + o[6];
+		if (s == 'ag' || s == 'ga' || s == 'ct' || s == 'tc')
+			++n_sub[1];
+		else ++n_sub[2];
+	}
+}
+
+var a = [], out = [];
+var c1_ctg = null, c1_start = 0, c1_end = 0, c1_counted = false, c1_len = 0;
+while (file.readline(buf) >= 0) {
+	var line = buf.toString();
+	if (!/\ts2:i:/.test(line)) continue; // skip secondary alignments
+	var m, t = line.split("\t", 12);
+	for (var i = 6; i <= 11; ++i)
+		t[i] = parseInt(t[i]);
+	if (t[10] < min_cov_len || t[11] < min_mapq) continue;
+	var ctg = t[5], x = t[7], end = t[8];
+	// compute regions covered by 1 contig
+	if (ctg != c1_ctg || x >= c1_end) {
+		if (c1_counted && c1_end > c1_start) {
+			c1_len += c1_end - c1_start;
+			print('R', c1_ctg, c1_start, c1_end);
+		}
+		c1_ctg = ctg, c1_start = x, c1_end = end;
+		c1_counted = (t[10] >= min_var_len);
+	} else if (end > c1_end) { // overlap
+		if (c1_counted && x > c1_start) {
+			c1_len += x - c1_start;
+			print('R', c1_ctg, c1_start, x);
+		}
+		c1_start = c1_end, c1_end = end;
+		c1_counted = (t[10] >= min_var_len);
+	} else { // contained
+		if (c1_counted && x > c1_start) {
+			c1_len += x - c1_start;
+			print('R', c1_ctg, c1_start, x);
+		}
+		c1_start = end;
+	}
+	// output variants ahead of this alignment
+	while (out.length) {
+		if (out[0][0] != ctg || out[0][2] <= x) {
+			count_var(out[0]);
+			print('V', out[0].join("\t"));
+			out.shift();
+		} else break;
+	}
+	// update coverage
+	for (var i = 0; i < out.length; ++i)
+		if (out[i][1] >= x && out[i][2] <= end)
+			++out[i][3];
+	// drop alignments that don't overlap with the current one
+	var k = 0;
+	for (var i = 0; i < a.length; ++i)
+		if (a[0][0] == ctg && a[0][2] > x)
+			a[k++] = a[i];
+	a.length = k;
+	// core loop
+	if (t[10] >= min_var_len) {
+		if ((m = /\tcs:Z:(\S+)/.exec(line)) == null) continue; // no cs tag
+		var cs = m[1];
+		var blen = 0, n_diff = 0;
+		tot_len += t[10];
+		while ((m = re_cs.exec(cs)) != null) {
+			var cov = 1;
+			if (m[1] == '*' || m[1] == '+' || m[1] == '-')
+				for (var i = 0; i < a.length; ++i)
+					if (a[0][2] > x) ++cov;
+			if (m[1] == '=' || m[1] == ':') {
+				var l = m[1] == '='? m[2].length : parseInt(m[2]);
+				x += l, blen += l;
+			} else if (m[1] == '*') {
+				out.push([t[5], x, x+1, cov, t[11], m[2].charAt(0), m[2].charAt(1)]);
+				++x, ++blen, ++n_diff;
+			} else if (m[1] == '+') {
+				out.push([t[5], x, x, cov, t[11], '-', m[2]]);
+				++blen, ++n_diff;
+			} else if (m[1] == '-') {
+				out.push([t[5], x, x + m[2].length, cov, t[11], m[2], '-']);
+				x += m[2].length, ++blen, ++n_diff;
+			}
+		}
+	}
+	a.push([t[5], t[7], t[8]]);
+}
+if (c1_counted && c1_end > c1_start) {
+	c1_len += c1_end - c1_start;
+	print('R', c1_ctg, c1_start, c1_end);
+}
+while (out.length) {
+	count_var(out[0]);
+	print('V', out[0].join("\t"));
+	out.shift();
+}
+
+//warn(tot_len + " alignment columns considered in calling");
+warn(c1_len + " reference bases covered by exactly one contig");
+warn(n_sub[0] + " substitutions; ts/tv = " + (n_sub[1]/n_sub[2]).toFixed(3));
+warn(n_del[0] + " 1bp deletions");
+warn(n_ins[0] + " 1bp insertions");
+warn(n_del[1] + " 2bp deletions");
+warn(n_ins[1] + " 2bp insertions");
+warn(n_del[2] + " [3,"+gap_thres+") deletions");
+warn(n_ins[2] + " [3,"+gap_thres+") insertions");
+warn(n_del[3] + " >="+gap_thres+" deletions");
+warn(n_ins[3] + " >="+gap_thres+" insertions");
+
+buf.destroy();
+file.close();

misc/splice2bed.js

@@ -0,0 +1,127 @@
+var getopt = function(args, ostr) {
+	var oli; // option letter list index
+	if (typeof(getopt.place) == 'undefined')
+		getopt.ind = 0, getopt.arg = null, getopt.place = -1;
+	if (getopt.place == -1) { // update scanning pointer
+		if (getopt.ind >= args.length || args[getopt.ind].charAt(getopt.place = 0) != '-') {
+			getopt.place = -1;
+			return null;
+		}
+		if (getopt.place + 1 < args[getopt.ind].length && args[getopt.ind].charAt(++getopt.place) == '-') { // found "--"
+			++getopt.ind;
+			getopt.place = -1;
+			return null;
+		}
+	}
+	var optopt = args[getopt.ind].charAt(getopt.place++); // character checked for validity
+	if (optopt == ':' || (oli = ostr.indexOf(optopt)) < 0) {
+		if (optopt == '-') return null; //  if the user didn't specify '-' as an option, assume it means null.
+		if (getopt.place < 0) ++getopt.ind;
+		return '?';
+	}
+	if (oli+1 >= ostr.length || ostr.charAt(++oli) != ':') { // don't need argument
+		getopt.arg = null;
+		if (getopt.place < 0 || getopt.place >= args[getopt.ind].length) ++getopt.ind, getopt.place = -1;
+	} else { // need an argument
+		if (getopt.place >= 0 && getopt.place < args[getopt.ind].length)
+			getopt.arg = args[getopt.ind].substr(getopt.place);
+		else if (args.length <= ++getopt.ind) { // no arg
+			getopt.place = -1;
+			if (ostr.length > 0 && ostr.charAt(0) == ':') return ':';
+			return '?';
+		} else getopt.arg = args[getopt.ind]; // white space
+		getopt.place = -1;
+		++getopt.ind;
+	}
+	return optopt;
+}
+
+var colors = ["0,128,255", "255,0,0", "0,192,0"];
+
+function print_lines(a, fmt) {
+	if (a.length == 0) return;
+	if (fmt == "bed") {
+		var n_pri = 0;
+		for (var i = 0; i < a.length; ++i)
+			if (a[i][8] == 0) ++n_pri;
+		if (n_pri > 1) {
+			for (var i = 0; i < a.length; ++i)
+				if (a[i][8] == 0) a[i][8] = 1;
+		} else if (n_pri == 0) {
+			warn("Warning: " + a[0][3] + " doesn't have a primary alignment");
+		}
+		for (var i = 0; i < a.length; ++i) {
+			a[i][8] = colors[a[i][8]];
+			print(a[i].join("\t"));
+		}
+	}
+	a.length = 0;
+}
+
+function main(args) {
+	var re = /(\d+)([MIDNSH])/g;
+	var c, fmt = "bed";
+	while ((c = getopt(args, "f:")) != null) {
+		if (c == 'f') fmt = getopt.arg;
+	}
+	if (getopt.ind == args.length) {
+		warn("Usage: k8 splice2bed.js <in.paf>");
+		exit(1);
+	}
+
+	var file = new File(args[getopt.ind]);
+	var buf = new Bytes();
+	var a = [];
+	while (file.readline(buf) >= 0) {
+		var line = buf.toString();
+		if (line.charAt(0) == '@') continue; // skip SAM header lines
+		var t = line.split("\t");
+		var is_pri = false, cigar = null, a1;
+		if (a.length && a[0][3] != t[0]) {
+			print_lines(a, fmt);
+			a = [];
+		}
+		if (t.length >= 12 && (t[4] == '+' || t[4] == '-')) {
+			for (var i = 12; i < t.length; ++i) {
+				if (t[i].substr(0, 5) == 'cg:Z:') {
+					cigar = t[i].substr(5);
+				} else if (t[i].substr(0, 5) == 's2:i:') {
+					is_pri = true;
+				}
+			}
+			a1 = [t[5], t[7], t[8], t[0], Math.floor(t[9]/t[10]*1000), t[4]];
+		} else if (t.length >= 10) {
+			var flag = parseInt(t[1]);
+			if ((flag&4) || a[2] == '*') continue;
+			cigar = t[5];
+			is_pri = (flag&0x100)? false : true;
+			a1 = [t[2], parseInt(t[3])-1, null, t[0], 1000, (flag&16)? '-' : '+'];
+		} else {
+			throw Error("unrecognized input format");
+		}
+		if (cigar == null) throw Error("missing CIGAR");
+		var m, x0 = 0, x = 0, bs = [], bl = [];
+		while ((m = re.exec(cigar)) != null) {
+			if (m[2] == 'M' || m[2] == 'D') {
+				x += parseInt(m[1]);
+			} else if (m[2] == 'N') {
+				bs.push(x0);
+				bl.push(x - x0);
+				x += parseInt(m[1]);
+				x0 = x;
+			}
+		}
+		bs.push(x0);
+		bl.push(x - x0);
+		// write the BED12 line
+		if (a1[2] == null) a1[2] = a1[1] + x;
+		a1.push(a1[1], a1[2]); // thick start/end is the same as start/end
+		a1.push(is_pri? 0 : 2, bs.length, bl.join(",")+",", bs.join(",")+",");
+		a.push(a1);
+	}
+	print_lines(a, fmt);
+	buf.destroy();
+	file.close();
+}
+
+main(arguments);

mmpriv.h

@@ -64,7 +64,7 @@ const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);
 int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
 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 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);
-mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, const char *qual, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
+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_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a);
@@ -75,10 +75,11 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff
 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(void *km, const mm_mapopt_t *opt, int *n_regs, mm_reg1_t *regs);
-void mm_filter_by_identity(void *km, int n_regs, mm_reg1_t *regs, float min_iden, int qlen, const char *qual);
 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_by_dp(void *km, int *n_regs, mm_reg1_t *r);
-void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len);
+void mm_set_mapq(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);
 
 mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int n_regs0, const mm_reg1_t *regs0, int *n_regs, mm_reg1_t **regs, const mm128_t *a);
 void mm_seg_free(void *km, int n_segs, mm_seg_t *segs);

pe.c

@@ -159,8 +159,8 @@ void mm_pair(void *km, int max_gap_ref, int pe_bonus, int sub_diff, int match_sc
 			mapq_pe_alt = (int)(6.02f * ((max>>32) - (sc.a[sc.n - 2]>>32)) / match_sc - 4.343f * logf(n_sub)); // n_sub > 0 because it counts the optimal, too
 			mapq_pe = mapq_pe < mapq_pe_alt? mapq_pe : mapq_pe_alt;
 		}
-		if (r[0]->mapq < mapq_pe) r[0]->mapq = (r[0]->mapq + mapq_pe) / 2;
-		if (r[1]->mapq < mapq_pe) r[1]->mapq = (r[1]->mapq + mapq_pe) / 2;
+		if (r[0]->mapq < mapq_pe) r[0]->mapq = (int)(.2f * r[0]->mapq + .8f * mapq_pe + .499f);
+		if (r[1]->mapq < mapq_pe) r[1]->mapq = (int)(.2f * r[1]->mapq + .8f * mapq_pe + .499f);
 		if (sc.n == 1) {
 			if (r[0]->mapq < 2) r[0]->mapq = 2;
 			if (r[1]->mapq < 2) r[1]->mapq = 2;

python/cmappy.pxd

@@ -5,7 +5,7 @@ cdef extern from "minimap.h":
 	# Options
 	#
 	ctypedef struct mm_idxopt_t:
-		short k, w, is_hpc, bucket_bits
+		short k, w, flag, bucket_bits
 		int mini_batch_size
 		uint64_t batch_size
 
@@ -22,7 +22,6 @@ cdef extern from "minimap.h":
 		float mask_level
 		float pri_ratio
 		int best_n
-		float min_iden
 		int max_join_long, max_join_short
 		int min_join_flank_sc
 		int a, b, q, e, q2, e2
@@ -31,6 +30,7 @@ cdef extern from "minimap.h":
 		int end_bonus
 		int min_dp_max
 		int min_ksw_len
+		int anchor_ext_len, anchor_ext_shift
 		int pe_ori, pe_bonus
 		float mid_occ_frac
 		int32_t mid_occ
@@ -52,7 +52,7 @@ cdef extern from "minimap.h":
 		pass
 
 	ctypedef struct mm_idx_t:
-		int32_t b, w, k, is_hpc
+		int32_t b, w, k, flag
 		uint32_t n_seq
 		mm_idx_seq_t *seq
 		uint32_t *S

python/mappy.pyx

@@ -13,7 +13,8 @@ cdef class Alignment:
 	cdef _ctg, _cigar # these are python objects
 
 	def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand):
-		self._ctg, self._ctg_len, self._r_st, self._r_en = str(ctg), cl, cs, ce
+		self._ctg = ctg if isinstance(ctg, str) else ctg.decode()
+		self._ctg_len, self._r_st, self._r_en = cl, cs, ce
 		self._strand, self._q_st, self._q_en = strand, qs, qe
 		self._NM, self._mlen, self._blen = NM, mlen, blen
 		self._mapq = mapq
@@ -34,7 +35,7 @@ cdef class Alignment:
 	def r_en(self): return self._r_en
 
 	@property
-	def strand(self): return self.strand
+	def strand(self): return self._strand
 
 	@property
 	def trans_strand(self): return self._trans_strand
@@ -150,9 +151,11 @@ def fastx_read(fn):
 	ks = cmappy.mm_fastx_open(str.encode(fn))
 	if ks is NULL: return None
 	while cmappy.kseq_read(ks) >= 0:
-		if ks.qual.l > 0: qual = str(ks.qual.s)
+		if ks.qual.l > 0: qual = ks.qual.s if isinstance(ks.qual.s, str) else ks.qual.s.decode()
 		else: qual = None
-		yield str(ks.name.s), str(ks.seq.s), qual
+		name = ks.name.s if isinstance(ks.name.s, str) else ks.name.s.decode()
+		seq = ks.seq.s if isinstance(ks.seq.s, str) else ks.seq.s.decode()
+		yield name, seq, qual
 	cmappy.mm_fastx_close(ks)
 
 def verbose(v=None):

sdust.c

@@ -56,6 +56,7 @@ sdust_buf_t *sdust_buf_init(void *km)
 	buf = (sdust_buf_t*)kcalloc(km, 1, sizeof(sdust_buf_t));
 	buf->km = km;
 	buf->w = kdq_init(int, buf->km);
+	kdq_resize(int, buf->w, 8);
 	return buf;
 }
 

setup.py

@@ -23,7 +23,7 @@ def readme():
 
 setup(
 	name = 'mappy',
-	version = '2.3',
+	version = '2.6',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),
@@ -35,7 +35,7 @@ setup(
     ext_modules = [Extension('mappy',
 		sources = [module_src, 'align.c', 'bseq.c', 'chain.c', 'format.c', 'hit.c', 'index.c', 'pe.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'],
+				   'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c', 'esterr.c'],
 		depends = ['minimap.h', 'bseq.h', 'kalloc.h', 'kdq.h', 'khash.h', 'kseq.h', 'ksort.h',
 				   'ksw2.h', 'kthread.h', 'kvec.h', 'mmpriv.h', 'sdust.h',
 				   'python/cmappy.h', 'python/cmappy.pxd'],

sketch.c

@@ -101,13 +101,13 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i
 				tq_push(&tq, skip_len);
 				kmer_span += skip_len;
 				if (tq.count > k) kmer_span -= tq_shift(&tq);
-				if (kmer_span >= 256) continue; // make sure $kmer_span does not take more than 8 bits
 			} else kmer_span = l + 1 < k? l + 1 : k;
 			kmer[0] = (kmer[0] << 2 | c) & mask;           // forward k-mer
 			kmer[1] = (kmer[1] >> 2) | (3ULL^c) << shift1; // reverse k-mer
 			if (kmer[0] == kmer[1]) continue; // skip "symmetric k-mers" as we don't know it strand
 			z = kmer[0] < kmer[1]? 0 : 1; // strand
-			if (++l >= k) {
+			++l;
+			if (l >= k && kmer_span < 256) {
 				info.x = hash64(kmer[z], mask) << 8 | kmer_span;
 				info.y = (uint64_t)rid<<32 | (uint32_t)i<<1 | z;
 			}
@@ -115,12 +115,12 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i
 		buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
 		if (l == w + k - 1) { // special case for the first window - because identical k-mers are not stored yet
 			for (j = buf_pos + 1; j < w; ++j)
-				if (min.x == buf[j].x && buf[j].y != min.y) kv_push(mm128_t, km, *p, buf[j]);
+				if (min.x == buf[j].x && buf[j].y != min.y && buf[j].y != UINT64_MAX) kv_push(mm128_t, km, *p, buf[j]);
 			for (j = 0; j < buf_pos; ++j)
-				if (min.x == buf[j].x && buf[j].y != min.y) kv_push(mm128_t, km, *p, buf[j]);
+				if (min.x == buf[j].x && buf[j].y != min.y && buf[j].y != UINT64_MAX) kv_push(mm128_t, km, *p, buf[j]);
 		}
 		if (info.x <= min.x) { // a new minimum; then write the old min
-			if (l >= w + k) kv_push(mm128_t, km, *p, min);
+			if (l >= w + k && min.y != UINT64_MAX) kv_push(mm128_t, km, *p, min);
 			min = info, min_pos = buf_pos;
 		} else if (buf_pos == min_pos) { // old min has moved outside the window
 			if (l >= w + k - 1) kv_push(mm128_t, km, *p, min);
@@ -130,9 +130,9 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i
 				if (min.x >= buf[j].x) min = buf[j], min_pos = j;
 			if (l >= w + k - 1) { // write identical k-mers
 				for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
-					if (min.x == buf[j].x && min.y != buf[j].y) kv_push(mm128_t, km, *p, buf[j]);
+					if (min.x == buf[j].x && min.y != buf[j].y && buf[j].y != UINT64_MAX) kv_push(mm128_t, km, *p, buf[j]);
 				for (j = 0; j <= buf_pos; ++j)
-					if (min.x == buf[j].x && min.y != buf[j].y) kv_push(mm128_t, km, *p, buf[j]);
+					if (min.x == buf[j].x && min.y != buf[j].y && buf[j].y != UINT64_MAX) kv_push(mm128_t, km, *p, buf[j]);
 			}
 		}
 		if (++buf_pos == w) buf_pos = 0;

tex/minimap2.bib

@@ -281,3 +281,27 @@
 	Journal = {arXiv:1111:5572},
 	Title = {Faster and More Accurate Sequence Alignment with SNAP},
 	Year = {2011}}
+
+@article{Irimia:2008aa,
+	Author = {Irimia, Manuel and Roy, Scott William},
+	Journal = {PLoS Genet},
+	Pages = {e1000148},
+	Title = {Evolutionary convergence on highly-conserved 3' intron structures in intron-poor eukaryotes and insights into the ancestral eukaryotic genome},
+	Volume = {4},
+	Year = {2008}}
+
+@article{Depristo:2011vn,
+	Author = {Depristo, Mark A and others},
+	Journal = {Nat Genet},
+	Pages = {491-8},
+	Title = {A framework for variation discovery and genotyping using next-generation {DNA} sequencing data},
+	Volume = {43},
+	Year = {2011}}
+
+@article{Kurtz:2004zr,
+	Author = {Kurtz, Stefan and others},
+	Journal = {Genome Biol},
+	Pages = {R12},
+	Title = {Versatile and open software for comparing large genomes},
+	Volume = {5},
+	Year = {2004}}

tex/minimap2.tex

@@ -31,10 +31,10 @@
 \section{Motivation:} Recent advances in sequencing technologies promise
 ultra-long reads of $\sim$100 kilo bases (kb) in average, full-length mRNA or
 cDNA reads in high throughput and genomic contigs over 100 mega bases (Mb) in
-length. Existing alignment tools are unable or inefficient to process such data
+length. Existing alignment programs are unable or inefficient to process such data
 at scale, which presses for the development of new alignment algorithms.
 
-\section{Results:} Minimap2 is a general-purpose aligner to map DNA or long
+\section{Results:} Minimap2 is a general-purpose alignment program to map DNA or long
 mRNA sequences against a large reference database. It works with accurate short
 reads of $\ge$100bp in length, $\ge$1kb genomic reads at error rate $\sim$15\%,
 full-length noisy Direct RNA or cDNA reads, and assembly contigs or closely
@@ -64,17 +64,24 @@ the thought that 10kb long sequences should be easier to map than 100bp reads
 because we can more effectively skip repetitive regions, which are often the
 bottleneck of short-read alignment. We confirmed our speculation by achieving
 approximate mapping 50 times faster than BWA-MEM~\citep{Li:2016aa}.
-\citet{Suzuki:2016} extended our work with a fast and novel algorithm on
+\citet{Suzuki130633} extended our work with a fast and novel algorithm on
 generating base-level alignment, which in turn inspired us to develop minimap2
-towards higher accuracy and more practical functionality.
+with added functionality.
 
 Both SMRT and ONT have been applied to the sequencing of spliced mRNAs (RNA-seq). While
 traditional mRNA aligners work~\citep{Wu:2005vn,Iwata:2012aa}, they are not
 optimized for long noisy sequence reads and are tens of times slower than
 dedicated long-read aligners. When developing minimap2 initially for aligning
-genomic DNA only, we realized minor modifications could make it competitive for
-aligning mRNAs as well. Minimap2 is a first RNA-seq aligner specifically
-designed for long noisy reads.
+genomic DNA only, we realized minor modifications could enable the base
+algorithm to map mRNAs as well. Minimap2 becomes a first RNA-seq aligner
+specifically designed for long noisy reads. We have also extended the original
+algorithm to map short reads at a speed faster than several mainstream
+short-read mappers.
+
+In this article, we will describe the minimap2 algorithm and its applications
+to different types of input sequences. We will evaluate the performance and
+accuracy of minimap2 on several simulated and real data sets and demonstrate
+the versatility of minimap2.
 
 \begin{methods}
 \section{Methods}
@@ -179,7 +186,7 @@ where $s(i,j)$ is the score between the $i$-th reference base and $j$-th query
 base. Eq.~(\ref{eq:ae86}) is a natural extension to the equation under affine
 gap cost~\citep{Gotoh:1982aa,Altschul:1986aa}.
 
-\subsubsection{Suzuki's formulation}
+\subsubsection{The Suzuki-Kasahara formulation}
 
 When we allow gaps longer than several hundred base pairs, nucleotide-level
 alignment is much slower than chaining. SSE acceleration is critical to the
@@ -187,7 +194,7 @@ performance of minimap2. Traditional SSE implementations~\citep{Farrar:2007hs}
 based on Eq.~(\ref{eq:ae86}) can achieve 16-way parallelization for short
 sequences, but only 4-way parallelization when the peak alignment score reaches
 32767. Long sequence alignment may exceed this threshold. Inspired by
-\citet{Wu:1996aa} and the following work, \citet{Suzuki:2016} proposed a
+\citet{Wu:1996aa} and the following work, \citet{Suzuki130633} proposed a
 difference-based formulation that lifted this limitation.
 In case of 2-piece gap cost, define
 \[
@@ -337,18 +344,24 @@ F_{i,j+1}= \max\{H_{ij}-q,F_{ij}\}-e\\
 \tilde{E}_{i+1,j}= \max\{H_{ij}-d(i)-\tilde{q},\tilde{E}_{ij}\}\\
 \end{array}\right.
 \end{equation}
-Let $T$ be the reference sequence. $d(i)$ is the cost of a non-canonical donor
-site, which takes 0 if $T[i+1,i+2]={\tt GT}$, or a positive number $p$
-otherwise. Similarly, $a(i)$ is the cost of a non-canonical acceptor site, which
-takes 0 if $T[i-1,i]={\tt AG}$, or $p$ otherwise. Eq.~(\ref{eq:splice}) is
-almost equivalent to the equation used by EXALIN~\citep{Zhang:2006aa} except
-that we allow insertions immediately followed by deletions and vice versa; in
-addition, we use Suzuki's diagonal formulation in actual implementation.
-
-%Given that $d_i$ and $a_i$
-%are a function of the reference sequence, it is possible to incorporate
-%splicing signals with more sophisticated models, such as positional weight
-%matrices. We have not tried this approach.
+Let $T$ be the reference sequence. $d(i)$ is computed as
+\[d(i)=\left\{\begin{array}{ll}
+0 & \mbox{if $T[i+1,i+3]$ is ${\tt GTA}$ or ${\tt GTG}$} \\
+p/2 & \mbox{if $T[i+1,i+3]$ is ${\tt GTC}$ or ${\tt GTT}$} \\
+p & \mbox{otherwise}
+\end{array}\right.\]
+where $T[i,j]$ extracts a substring of $T$ between $i$ and $j$ inclusively.
+$d(i)$ penalizes non-canonical donor sites with $p$ and less frequent Eukayotic
+splicing signal ${\tt GT[C/T]}$ with $p/2$~\citep{Irimia:2008aa}. Similarly,
+\[a(i)=\left\{\begin{array}{ll}
+0 & \mbox{if $T[i-2,i]$ is ${\tt CAG}$ or ${\tt TAG}$} \\
+p/2 & \mbox{if $T[i-2,i]$ is ${\tt AAG}$ or ${\tt GAG}$} \\
+p & \mbox{otherwise}
+\end{array}\right.\]
+models the acceptor signal. Eq.~(\ref{eq:splice}) is close to an equation in
+\citet{Zhang:2006aa} except that we allow insertions immediately followed by
+deletions and vice versa; in addition, we use the Suzuki-Kasahara diagonal
+formulation in actual implementation.
 
 If RNA-seq reads are not sequenced from stranded libraries, the read strand
 relative to the underlying transcript is unknown. By default, minimap2 aligns
@@ -360,12 +373,12 @@ reads that span canonical splicing sites.
 
 In the spliced alignment mode, minimap2 further increases the density of
 minimizers and disables banded alignment. Together with the two-round DP-based
-alignment, spliced alignment is several times slower than DNA sequence
+alignment, spliced alignment is several times slower than genomic DNA
 alignment.
 
 \subsection{Aligning short paired-end reads}
 
-During chainging, minimap2 takes a pair of reads as one read with a gap of
+During chainging, minimap2 takes a pair of reads as one fragment with a gap of
 unknown length in the middle. It applies a normal gap cost between seeds on the
 same read but is a more permissive gap cost between seeds on different reads.
 More precisely, the gap cost during chaining is:
@@ -417,9 +430,7 @@ NGMLR~(v0.2.5; \citealp{Sedlazeck169557}). We excluded rHAT~\citep{Liu:2016ab}
 and LAMSA~\citep{Liu:2017aa} because they either
 crashed or produced malformatted output. In this evaluation, minimap2 has
 higher power to distinguish unique and repetitive hits, and achieves overall
-higher mapping accuracy (Fig.~\ref{fig:eval}a). It is still the most accurate
-even if we skip DP-based alignment (data not shown), confirming chaining alone
-is sufficient to achieve high accuracy for approximate mapping. Minimap2 and
+higher mapping accuracy (Fig.~\ref{fig:eval}a). Minimap2 and
 NGMLR provide better mapping quality estimate: they rarely give repetitive hits
 high mapping quality.  Apparently, other aligners may
 occasionally miss close suboptimal hits and be overconfident in wrong mappings.
@@ -440,16 +451,16 @@ to the 2-piece affine gap cost.
 \subsection{Aligning long spliced reads}
 
 We evaluated minimap2 on SIRV control data~(AC:SRR5286959;
-\citealp{Byrne:2017aa}) where the truth is known. Minimap2 predicted 59\,916
-introns from 11\,017 reads. 93.0\% of splice juctions are precise. We examined
+\citealp{Byrne:2017aa}) where the truth is known. Minimap2 predicted 59\,918
+introns from 11\,018 reads. 93.8\% of splice juctions are precise. We examined
 wrongly predicted junctions and found the majority were caused by clustered
 splicing signals (e.g. two adjacent ${\tt GT}$ sites). When INDEL sequencing
 errors are frequent, it is difficult to find precise splicing sites in this
 case. If we allow up to 10bp distance from true splicing sites, 98.4\% of
-aligned introns are approximately correct.  Given this observation, we might be
-able to improve boundary detection by initializing $d(\cdot)$ and $a(\cdot)$ in
-Eq.~(\ref{eq:splice}) with position-specific scoring matrices or more
-sophisticated models. We have not tried this approach.
+aligned introns are approximately correct. It is worth noting that for SIRV, we
+asked minimap2 to model the ${\tt GT..AG}$ splicing signal only without extra
+bases. This is because SIRV does not honor the evolutionarily prevalent signal
+${\tt GT[A/G]..[C/T]AG}$~\citep{Irimia:2008aa}.
 
 \begin{table}[!tb]
 \processtable{Evaluation of junction accuracy on 2D ONT reads}
@@ -460,13 +471,13 @@ sophisticated models. We have not tried this approach.
 \midrule
 Run time (CPU min)        & 631      & 15.9     & 2\,076   & 33.9 \\
 Peak RAM (GByte)          & 8.9      & 14.5     & 3.2      & 29.2\vspace{1em}\\
-\# aligned reads          & 103\,669 & 104\,200 & 103\,711 & 26\,479 \\
+\# aligned reads          & 103\,669 & 104\,199 & 103\,711 & 26\,479 \\
 \# chimeric alignments    & 1\,904   & 1\,488   & 0        & 0 \\
-\# non-spliced alignments & 15\,854  & 14\,639  & 17\,033  & 10\,545\vspace{1em}\\
-\# aligned introns        & 692\,275 & 694\,103 & 692\,945 & 78\,603 \\
-\# novel introns          & 11\,239  & 3\,207   & 8\,550   & 1\,214 \\
-\% exact introns          & 83.8\%   & 91.7\%   & 87.9\%   & 55.2\% \\
-\% approx. introns        & 91.8\%   & 96.5\%   & 92.5\%   & 82.4\% \\
+\# non-spliced alignments & 15\,854  & 14\,798  & 17\,033  & 10\,545\vspace{1em}\\
+\# aligned introns        & 692\,275 & 693\,553 & 692\,945 & 78\,603 \\
+\# novel introns          & 11\,239  & 3\,113   & 8\,550   & 1\,214 \\
+\% exact introns          & 83.8\%   & 94.0\%   & 87.9\%   & 55.2\% \\
+\% approx. introns        & 91.8\%   & 96.9\%   & 92.5\%   & 82.4\% \\
 \botrule
 \end{tabular}
 }{Mouse reads (AC:SRR5286960) were mapped to the primary assembly of mouse
@@ -487,10 +498,16 @@ STAR~(v2.5.3a; \citealp{Dobin:2013kx}). In general, minimap2 is more
 consistent with existing annotations (Table~\ref{tab:intron}): it finds
 more junctions with a higher percentage being exactly or approximately correct.
 Minimap2 is over 40 times faster than GMAP and SpAln. While STAR is close to
-minimap2 in speed, it does not work well with noisy reads.  We have also
-evaluated spliced aligners on public Iso-Seq data (human Alzheimer brain
-from \href{http://bit.ly/isoseqpub}{http://bit.ly/isoseqpub}). The observation
-is similar: minimap2 is faster at higher junction accuracy.
+minimap2 in speed, it does not work well with noisy reads.
+
+We have also evaluated spliced aligners on public Iso-Seq data (human Alzheimer
+brain from \href{http://bit.ly/isoseqpub}{http://bit.ly/isoseqpub}). The
+observation is similar: minimap2 is faster at higher junction accuracy.
+On a private Nanopore Direct RNA data set with $\sim$17\% sequencing error rate
+(N. Loman, personal communication), minimap2 aligned 96\,467 introns
+from 37\,068 mapped reads with 95.4\% of them consistent with human gene
+annotations. In comparison, only 74.8\% of GMAP introns found in known gene
+annotations.
 
 We noted that GMAP and SpAln have not been optimized for noisy reads. We are
 showing the best setting we have experimented, but their developers should be
@@ -518,8 +535,8 @@ able to improve their accuracy further.
 
 \subsection{Aligning short genomic reads}
 
-We evaluated minimap2 along with Bowtie2~\citep{Langmead:2012fk}, BWA-MEM and
-SNAP~\citep{Zaharia:2011aa}. Minimap2 is 3--4 times as fast as Bowtie2 and
+We evaluated minimap2 along with Bowtie2~(v2.3.3; \citealt{Langmead:2012fk}), BWA-MEM and
+SNAP (v1.0beta23; \citealt{Zaharia:2011aa}). Minimap2 is 3--4 times as fast as Bowtie2 and
 BWA-MEM, but is 1.3 times slower than SNAP. Minimap2 is more accurate on this
 simulated data set than Bowtie2 and SNAP but less accurate than BWA-MEM
 (Fig.~\ref{fig:eval}b). Closer investigation reveals that BWA-MEM achieves
@@ -528,23 +545,68 @@ region close to its mate. If we disable this feature, BWA-MEM becomes slightly
 less accurate than minimap2. We might consider to implement a similar heuristic
 in minimap2 in future.
 
-\section{Conclusion}
+To evaluate the accuracy of minimap2 on real data, we aligned human reads
+(AC:ERR1341796) with BWA-MEM and minimap2, and called SNPs and small INDELs
+with GATK HaplotypeCaller v3.5~\citep{Depristo:2011vn}. This run was sequenced
+from experimentally mixed CHM1 and CHM13 cell lines. Both of them are homozygous
+across the whole genome and have been \emph{de novo} assembled with SMRT reads
+to high quality. This allowed us to construct an independent truth variant
+data set
+(\href{https://github.com/lh3/CHM-eval}{https://github.com/lh3/CHM-eval}) for
+ERR1341796. In this evaluation, minimap2 has higher SNP false negative rate
+(FNR; 2.5\% of minimap2 vs 2.2\% of BWA-MEM), but fewer false positive SNPs per
+million bases (FPPM; 3.0 vs 3.9), lower 2--50bp INDEL FNR (7.3\% vs 7.5\%) and
+similar INDEL FPPM (both 1.0). Minimap2 is broadly similar to BWA-MEM in the
+context of small variant calling.
 
-Minimap2 is a fast, accurate and versatile aligner for long nucleotide
-sequences.  In addition to reference-based read mapping, minimap2 inherits
-minimap's functionality to search against huge multi-species databases and to
-find read overlaps. On a few test data sets, minimap2 appears to yield slightly
-better miniasm assembly~\citep{Li:2016aa}. Minimap2 can also align similar
-genomes or different assemblies of the same species. However, full-genome
-alignment is an intricate research topic. More thorough evaluations would be
-necessary to justify the use of minimap2 for such applications.
+\subsection{Other applications}
+
+Minimap2 retains minimap's functionality to find overlaps between long reads
+and to search against large multi-species databases such as \emph{nt} from
+NCBI. Minimap2 can also align similar genomes or different assemblies of the
+same species. It took 7 wall-clock minutes over 8 CPU cores to align a human
+SMRT assembly (AC:GCA\_001297185.1) to GRCh38, over 20 times faster
+MUMmer4~\citep{Kurtz:2004zr}.
+
+\section{Discussions}
+
+Minimap2 is a versatile mapper and pairwise aligner for nucleotide sequences.
+It works with short reads, assembly contigs and long noisy genomic and RNA-seq
+reads, and can be used as a read mapper, long-read overlapper or a full-genome
+aligner. Minimap2 is also accurate and efficient, often outperforming other
+domain-specific alignment tools in terms of both speed and accuracy.
+
+The capability of minimap2 comes from a fast base-level alignment algorithm and
+an accurate chaining algorithm. When aligning long query sequences, base-level
+alignment is often the performance bottleneck. The Suzuki-Kasahara algorithm
+greatly alleviates the bottleneck and enables DP-based splice alignment
+involving $>$100kb introns, which was impractically slow ten years ago.  The
+minimap2 chaining algorithm is fast and highly accurate by itself.  In fact,
+chaining alone is more accurate than all the other long-read mappers in
+Fig.~\ref{fig:eval}a (data not shown). This accuracy helps to reduce downstream
+base-level alignment of candidate chains, which is still times slower than
+chaining even with the Suzuki-Kasahara improvement. In addition, taking a
+general form, minimap2 chaining can be adapted to non-typical data types such
+spliced reads and multiple reads per fragment. This gives us the opportunity to
+extend the same base algorithm to a variety of use cases.
+
+Modern mainstream aligners often use a full-text index, such as suffix array or
+FM-index, to index reference sequences. An advantage of this approach is that
+we can use exact seeds of arbitrary lengths, which helps to increase seed
+uniqueness and reduce unsuccessful extensions. Minimap2 indexes reference
+k-mers with a hash table instead. Such fixed-length seeds are inferior to
+variable-length seeds in theory, but can be computed much more efficiently in
+practice. When a query sequence has multiple seed hits, we can afford to skip
+some highly repetitive seeds without affecting the final accuracy. This further
+alleviates the concern with the uniqueness of seeds. Hash table is the ideal
+data structure for mapping long query sequences.
 
 \section*{Acknowledgements}
-We owe a debt of gratitude to Hajime Suzuki for releasing his masterpiece and
-insightful notes before formal publication. We thank M. Schatz, P. Rescheneder
-and F.  Sedlazeck for pointing out the limitation of BWA-MEM. We are also
-grateful to early minimap2 testers who have greatly helped to suggest features
-and to fix various issues.
+We owe a debt of gratitude to H. Suzuki and M. Kasahara for releasing their
+masterpiece and insightful notes before formal publication. We thank M.
+Schatz, P. Rescheneder and F. Sedlazeck for pointing out the limitation of
+BWA-MEM. We are also grateful to early minimap2 testers who have greatly helped
+to suggest features and to fix various issues.
 
 \bibliography{minimap2}
 

tex/mm2-s3.sam.eval

@@ -1,60 +1,62 @@
-Q	60	18345673	8	0.000000436	18345673
-Q	59	33966	4	0.000000653	18379639
-Q	58	34178	1	0.000000706	18413817
-Q	56	49138	1	0.000000758	18462955
-Q	54	22442	4	0.000000974	18485397
-Q	53	19070	2	0.000001081	18504467
-Q	52	14169	3	0.000001242	18518636
-Q	51	13233	4	0.000001457	18531869
-Q	50	12133	2	0.000001564	18544002
-Q	49	11138	4	0.000001778	18555140
-Q	48	11174	8	0.000002208	18566314
-Q	47	17139	4	0.000002422	18583453
-Q	46	20428	10	0.000002956	18603881
-Q	45	16503	3	0.000003115	18620384
-Q	44	11933	6	0.000003435	18632317
-Q	43	25392	11	0.000004020	18657709
-Q	42	16734	9	0.000004498	18674443
-Q	41	13826	10	0.000005030	18688269
-Q	40	13023	10	0.000005561	18701292
-Q	39	12686	10	0.000006092	18713978
-Q	38	17275	4	0.000006300	18731253
-Q	37	17241	2	0.000006401	18748494
-Q	36	12458	12	0.000007036	18760952
-Q	35	11981	5	0.000007298	18772933
-Q	34	12004	11	0.000007879	18784937
-Q	33	12111	7	0.000008246	18797048
-Q	32	11782	9	0.000008719	18808830
-Q	31	11811	7	0.000009086	18820641
-Q	30	33507	32	0.000010767	18854148
-Q	29	11243	21	0.000011874	18865391
-Q	28	10779	17	0.000012767	18876170
-Q	27	15733	24	0.000014027	18891903
-Q	26	16762	40	0.000016130	18908665
-Q	25	13811	49	0.000018708	18922476
-Q	24	14141	46	0.000021123	18936617
-Q	23	13429	55	0.000024010	18950046
-Q	22	13116	26	0.000025365	18963162
-Q	21	13436	46	0.000027771	18976598
-Q	20	13441	55	0.000030648	18990039
-Q	19	12988	53	0.000033416	19003027
-Q	18	13353	51	0.000036074	19016380
-Q	17	13782	77	0.000040094	19030162
-Q	16	14065	94	0.000045001	19044227
-Q	15	14044	124	0.000051474	19058271
-Q	14	14714	140	0.000058774	19072985
-Q	13	17459	197	0.000069040	19090444
-Q	12	17339	259	0.000082532	19107783
-Q	11	17381	280	0.000097097	19125164
-Q	10	17732	295	0.000112418	19142896
-Q	9	17959	416	0.000134023	19160855
-Q	8	18234	530	0.000161530	19179089
-Q	7	19048	514	0.000188143	19198137
-Q	6	19722	656	0.000222085	19217859
-Q	5	19753	775	0.000262143	19237612
-Q	4	19818	1030	0.000315359	19257430
-Q	3	17088	1100	0.000372149	19274518
-Q	2	43045	6708	0.000718569	19317563
-Q	1	126377	25255	0.002012761	19443940
-Q	0	554357	372087	0.020562901	19998297
+Q	60	18579866	27	0.000001453	18579866
+Q	59	27087	4	0.000001666	18606953
+Q	58	21435	1	0.000001718	18628388
+Q	57	45663	3	0.000001874	18674051
+Q	56	36031	2	0.000001978	18710082
+Q	55	18499	2	0.000002082	18728581
+Q	54	14754	2	0.000002187	18743335
+Q	53	25541	2	0.000002291	18768876
+Q	52	26397	5	0.000002554	18795273
+Q	51	15090	3	0.000002711	18810363
+Q	50	13425	11	0.000003294	18823788
+Q	49	15175	2	0.000003397	18838963
+Q	48	19407	4	0.000003606	18858370
+Q	47	11538	16	0.000004452	18869908
+Q	46	12558	17	0.000005349	18882466
+Q	45	40362	28	0.000006817	18922828
+Q	44	10465	13	0.000007500	18933293
+Q	43	10098	20	0.000008552	18943391
+Q	42	10682	19	0.000009549	18954073
+Q	41	9823	11	0.000010125	18963896
+Q	40	9685	16	0.000010963	18973581
+Q	39	10273	18	0.000011905	18983854
+Q	38	9515	18	0.000012847	18993369
+Q	37	9474	27	0.000014261	19002843
+Q	36	10430	25	0.000015568	19013273
+Q	35	9241	34	0.000017348	19022514
+Q	34	9162	31	0.000018968	19031676
+Q	33	10164	49	0.000021532	19041840
+Q	32	9152	55	0.000024408	19050992
+Q	31	9252	35	0.000026233	19060244
+Q	30	9872	55	0.000029103	19070116
+Q	29	8938	65	0.000032496	19079054
+Q	28	8951	73	0.000036306	19088005
+Q	27	9949	95	0.000041261	19097954
+Q	26	9784	97	0.000046316	19107738
+Q	25	10126	97	0.000051366	19117864
+Q	24	11260	123	0.000057765	19129124
+Q	23	10047	114	0.000063691	19139171
+Q	22	9661	123	0.000070083	19148832
+Q	21	10339	168	0.000078813	19159171
+Q	20	17928	193	0.000088804	19177099
+Q	19	9842	193	0.000098817	19186941
+Q	18	14737	247	0.000111605	19201678
+Q	17	10218	238	0.000123934	19211896
+Q	16	10271	242	0.000136457	19222167
+Q	15	12241	333	0.000153683	19234408
+Q	14	9189	336	0.000171070	19243597
+Q	13	9493	515	0.000197734	19253090
+Q	12	11502	743	0.000236185	19264592
+Q	11	8211	507	0.000262390	19272803
+Q	10	9133	606	0.000293695	19281936
+Q	9	10014	931	0.000341801	19291950
+Q	8	8436	698	0.000377816	19300386
+Q	7	8443	705	0.000414163	19308829
+Q	6	10203	944	0.000462808	19319032
+Q	5	6936	756	0.000501760	19325968
+Q	4	6732	843	0.000545190	19332700
+Q	3	8215	1104	0.000602040	19340915
+Q	2	21201	5440	0.000882342	19362116
+Q	1	82328	22186	0.002019600	19444444
+Q	0	553853	371953	0.020562901	19998297
 U	1703

tex/mm2.approx.eval

@@ -1,30 +1,12 @@
-Q	60	32066	0	0.000000000
-Q	40	32	1	0.000031155
-Q	38	19	1	0.000062272
-Q	36	11	1	0.000093376
-Q	35	32	1	0.000124378
-Q	33	15	1	0.000155400
-Q	32	58	1	0.000186145
-Q	27	11	1	0.000217095
-Q	26	80	1	0.000247494
-Q	21	19	2	0.000309186
-Q	20	16	1	0.000339936
-Q	19	19	1	0.000370622
-Q	18	22	2	0.000432099
-Q	17	37	5	0.000585751
-Q	15	24	2	0.000646930
-Q	14	18	3	0.000738939
-Q	13	30	6	0.000922821
-Q	12	18	1	0.000953054
-Q	11	29	2	0.001013638
-Q	10	30	1	0.001043393
-Q	9	20	5	0.001196099
-Q	8	25	8	0.001440348
-Q	7	28	6	0.001622830
-Q	6	35	12	0.001988132
-Q	5	34	12	0.002352725
-Q	4	29	8	0.002594865
-Q	3	36	14	0.003018937
-Q	2	46	15	0.003471482
-Q	1	69	36	0.004558162
-Q	0	167	94	0.007377173
+Q	60	32084	0	0.000000000	32084
+Q	24	318	2	0.000061725	32402
+Q	11	98	2	0.000123077	32500
+Q	8	37	2	0.000184405	32537
+Q	7	37	3	0.000276294	32574
+Q	6	40	3	0.000367940	32614
+Q	5	34	2	0.000428816	32648
+Q	4	37	5	0.000581306	32685
+Q	3	28	6	0.000764222	32713
+Q	2	38	6	0.000946536	32751
+Q	1	50	21	0.001585318	32801
+Q	0	286	150	0.006105117	33087

tex/mm2.eval

@@ -1,9 +1,13 @@
-Q	60	32226	0	0.000000000	32226
-Q	20	267	1	0.000030776	32493
-Q	10	34	1	0.000061487	32527
-Q	9	118	1	0.000091898	32645
-Q	5	27	2	0.000153036	32672
-Q	4	68	2	0.000213806	32740
-Q	1	314	101	0.003267381	33054
+Q	60	32477	0	0.000000000	32477
+Q	22	16	1	0.000030776	32493
+Q	21	44	1	0.000061468	32537
+Q	19	73	1	0.000091996	32610
+Q	14	66	1	0.000122414	32676
+Q	10	26	3	0.000214054	32702
+Q	8	14	1	0.000244529	32716
+Q	7	13	2	0.000305539	32729
+Q	6	47	1	0.000335611	32776
+Q	3	10	1	0.000366010	32786
+Q	2	20	2	0.000426751	32806
+Q	1	248	94	0.003267381	33054
 Q	0	31	17	0.003778147	33085
-U	3