Release minimap2-2.5 (r572)

Fix typo in strand property

NEWS.md

@@ -1,3 +1,64 @@
+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 +87,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,4 +1,4 @@
-[![Release](https://img.shields.io/badge/Release-v2.3-blue.svg?style=flat)](https://github.com/lh3/minimap2/releases)
+[![Release](https://img.shields.io/badge/Release-v2.5-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)
 [![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)
@@ -9,15 +9,20 @@
 ```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 +39,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 +71,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.5/minimap2-2.5_x64-linux.tar.bz2 \
   | tar -jxvf -
-./minimap2-2.2_x64-linux/minimap2
+./minimap2-2.5_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 +137,30 @@ 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.
+
+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 +200,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 +212,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`.

align.c

@@ -110,8 +110,11 @@ 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);
 	}
@@ -130,27 +133,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 +156,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;
@@ -381,6 +377,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

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];

format.c

@@ -203,7 +203,6 @@ 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);
 	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)
@@ -258,7 +257,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 +266,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 +335,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 +344,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 +381,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 +429,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,7 +76,7 @@ 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;
@@ -264,45 +264,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 +403,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

@@ -267,7 +267,6 @@ static void *worker_pipeline(void *shared, int step, void *in)
 				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
 					seq->name = (char*)kmalloc(p->mi->km, strlen(s->seq[i].name) + 1);
 					strcpy(seq->name, s->seq[i].name);
 				} else seq->name = 0;

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.5-r572"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -34,13 +34,15 @@ 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 },
 	{ "help",           no_argument,       0, 'h' },
 	{ "max-intron-len", required_argument, 0, 'G' },
 	{ "version",        no_argument,       0, 'V' },
@@ -65,7 +67,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;
@@ -100,7 +102,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 +109,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 +118,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 +134,10 @@ 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 == 14) { // --frag
 			if (optarg == 0 || strcmp(optarg, "yes") == 0 || strcmp(optarg, "y") == 0)
 				opt.flag |= MM_F_FRAG_MODE;
@@ -154,6 +157,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)
@@ -219,7 +226,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 +233,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);

map.c

@@ -105,10 +105,10 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		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;
+		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;
@@ -270,7 +270,7 @@ 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 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;
@@ -298,7 +298,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) {
@@ -345,7 +345,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 
 	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,15 +354,12 @@ 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, a);
 	kfree(b->km, u);
@@ -414,7 +411,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;

minimap.h

@@ -22,6 +22,9 @@
 #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_IDX_MAGIC   "MMI\2"
 
@@ -57,14 +60,12 @@ 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 cnt:28, rev:1, seg_split:1, sam_pri:1, proper_frag: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
@@ -72,7 +73,7 @@ typedef struct {
 	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;
+	uint32_t pe_thru:1, score0:31;
 	uint32_t hash;
 	mm_extra_t *p;
 } mm_reg1_t;
@@ -99,7 +100,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;

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "22 October 2017" "minimap2-2.2-dirty (r531)" "Bioinformatics tools"
+.TH minimap2 1 "11 November 2017" "minimap2-2.5 (r572)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -219,6 +219,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 +244,26 @@ 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.
 .SS Input/output options
 .TP 10
 .B -a
@@ -261,7 +281,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 +302,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 +394,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 '

mmpriv.h

@@ -75,10 +75,9 @@ 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);
 
 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

@@ -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

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):

setup.py

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

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;
 			}

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