Release minimap2-2.17 (r941)

for high-quality CCS/mRNA splice alignment

NEWS.md

@@ -1,3 +1,80 @@
+Release 2.17-r941 (4 May 2019)
+------------------------------
+
+Changes since the last release:
+
+ * Fixed flawed CIGARs like `5I6D7I` (#392).
+
+ * Bugfix: TLEN should be 0 when either end is unmapped (#373 and #365).
+
+ * Bugfix: mappy is unable to write index (#372).
+
+ * Added option `--junc-bed` to load known gene annotations in the BED12
+   format. Minimap2 prefers annotated junctions over novel junctions (#197 and
+   #348). GTF can be converted to BED12 with `paftools.js gff2bed`.
+
+ * Added option `--sam-hit-only` to suppress unmapped hits in SAM (#377).
+
+ * Added preset `splice:hq` for high-quality CCS or mRNA sequences. It applies
+   better scoring and improves the sensitivity to small exons. This preset may
+   introduce false small introns, but the overall accuracy should be higher.
+
+This version produces nearly identical alignments to v2.16, except for CIGARs
+affected by the bug mentioned above.
+
+(2.17: 5 May 2019, r941)
+
+
+
+Release 2.16-r922 (28 February 2019)
+------------------------------------
+
+This release is 50% faster for mapping ultra-long nanopore reads at comparable
+accuracy. For short-read mapping, long-read overlapping and ordinary long-read
+mapping, the performance and accuracy remain similar. This speedup is achieved
+with a new heuristic to limit the number of chaining iterations (#324). Users
+can disable the heuristic by increasing a new option `--max-chain-iter` to a
+huge number.
+
+Other changes to minimap2:
+
+ * Implemented option `--paf-no-hit` to output unmapped query sequences in PAF.
+   The strand and reference name columns are both `*` at an unmapped line. The
+   hidden option is available in earlier minimap2 but had a different 2-column
+   output format instead of PAF.
+
+ * Fixed a bug that leads to wrongly calculated `de` tags when ambiguous bases
+   are involved (#309). This bug only affects v2.15.
+
+ * Fixed a bug when parsing command-line option `--splice` (#344). This bug was
+   introduced in v2.13.
+
+ * Fixed two division-by-zero cases (#326). They don't affect final alignments
+   because the results of the divisions are not used in both case.
+
+ * Added an option `-o` to output alignments to a specified file. It is still
+   recommended to use UNIX pipes for on-the-fly conversion or compression.
+
+ * Output a new `rl` tag to give the length of query regions harboring
+   repetitive seeds.
+
+Changes to paftool.js:
+
+ * Added a new option to convert the MD tag to the long form of the cs tag.
+
+Changes to mappy:
+
+ * Added the `mappy.Aligner.seq_names` method to return sequence names (#312).
+
+For NA12878 ultra-long reads, this release changes the alignments of <0.1% of
+reads in comparison to v2.15. All these reads have highly fragmented alignments
+and are likely to be problematic anyway. For shorter or well aligned reads,
+this release should produce mostly identical alignments to v2.15.
+
+(2.16: 28 February 2019, r922)
+
+
+
 Release 2.15-r905 (10 January 2019)
 -----------------------------------
 

README.md

@@ -9,8 +9,8 @@ cd minimap2 && make
 # 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
+./minimap2 -x map-ont -d MT-human-ont.mmi test/MT-human.fa
+./minimap2 -a MT-human-ont.mmi test/MT-orang.fa > test.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
@@ -18,7 +18,7 @@ cd minimap2 && make
 ./minimap2 -ax sr ref.fa read1.fa read2.fa > aln.sam      # short genomic paired-end reads
 ./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads (strand unknown)
 ./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam  # noisy Nanopore Direct RNA-seq
-./minimap2 -ax splice -uf -C5 ref.fa query.fa > aln.sam   # Final PacBio Iso-seq or traditional cDNA
+./minimap2 -ax splice:hq -uf ref.fa query.fa > aln.sam    # Final PacBio Iso-seq or traditional cDNA
 ./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
@@ -71,8 +71,8 @@ Detailed evaluations are available from the [minimap2 paper][doi] or the
 Minimap2 is optimized for x86-64 CPUs. You can acquire precompiled binaries from
 the [release page][release] with:
 ```sh
-curl -L https://github.com/lh3/minimap2/releases/download/v2.15/minimap2-2.15_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.15_x64-linux/minimap2
+curl -L https://github.com/lh3/minimap2/releases/download/v2.17/minimap2-2.17_x64-linux.tar.bz2 | tar -jxvf -
+./minimap2-2.17_x64-linux/minimap2
 ```
 If you want to compile from the source, you need to have a C compiler, GNU make
 and zlib development files installed. Then type `make` in the source code
@@ -139,7 +139,7 @@ Nanopore reads.
 #### <a name="map-long-splice"></a>Map long mRNA/cDNA reads
 
 ```sh
-minimap2 -ax splice -uf -C5 ref.fa iso-seq.fq > aln.sam      # PacBio Iso-seq/traditional cDNA
+minimap2 -ax splice:hq -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
@@ -324,7 +324,7 @@ There is not a specific mailing list for the time being.
 If you use minimap2 in your work, please cite:
 
 > Li, H. (2018). Minimap2: pairwise alignment for nucleotide sequences.
-> Bioinformatics. [doi:10.1093/bioinformatics/bty191][doi]
+> *Bioinformatics*, **34**:3094-3100. [doi:10.1093/bioinformatics/bty191][doi]
 
 ## <a name="dguide"></a>Developers' Guide
 

align.c

@@ -123,6 +123,25 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 		}
 	}
 	assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
+	for (k = 0; k < p->n_cigar - 2; ++k) { // fix CIGAR like 5I6D7I
+		if ((p->cigar[k]&0xf) > 0 && (p->cigar[k]&0xf) + (p->cigar[k+1]&0xf) == 3) {
+			uint32_t l, s[3] = {0,0,0};
+			for (l = k; l < p->n_cigar; ++l) { // count number of adjacent I and D
+				uint32_t op = p->cigar[l]&0xf;
+				if (op == 1 || op == 2 || p->cigar[l]>>4 == 0)
+					s[op] += p->cigar[l] >> 4;
+				else break;
+			}
+			if (s[1] > 0 && s[2] > 0 && l - k > 2) { // turn to a single I and a single D
+				p->cigar[k]   = s[1]<<4|1;
+				p->cigar[k+1] = s[2]<<4|2;
+				for (k += 2; k < l; ++k)
+					p->cigar[k] &= 0xf;
+				to_shrink = 1;
+			}
+			k = l;
+		}
+	}
 	if (to_shrink) { // squeeze out zero-length operations
 		int32_t l = 0;
 		for (k = 0; k < p->n_cigar; ++k) // squeeze out zero-length operations
@@ -291,7 +310,7 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) //
 	}
 }
 
-static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
+static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const uint8_t *junc, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
 {
 	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
 		int i;
@@ -305,7 +324,7 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 		ksw_reset_extz(ez);
 		ez->zdropped = 1;
 	} else if (opt->flag & MM_F_SPLICE)
-		ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, flag, ez);
+		ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, opt->junc_bonus, flag, junc, ez);
 	else if (opt->q == opt->q2 && opt->e == opt->e2)
 		ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, flag, ez);
 	else
@@ -547,7 +566,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 {
 	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;
-	uint8_t *tseq, *qseq;
+	uint8_t *tseq, *qseq, *junc;
 	int32_t i, l, bw, dropped = 0, extra_flag = 0, rs0, re0, qs0, qe0;
 	int32_t rs, re, qs, qe;
 	int32_t rs1, qs1, re1, qe1;
@@ -666,13 +685,16 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 
 	assert(re0 > rs0);
 	tseq = (uint8_t*)kmalloc(km, re0 - rs0);
+	junc = (uint8_t*)kmalloc(km, re0 - rs0);
 
 	if (qs > 0 && rs > 0) { // left extension; probably the condition can be changed to "qs > qs0 && rs > rs0"
 		qseq = &qseq0[rev][qs0];
 		mm_idx_getseq(mi, rid, rs0, rs, tseq);
+		mm_idx_bed_junc(mi, rid, rs0, rs, junc);
 		mm_seq_rev(qs - qs0, qseq);
 		mm_seq_rev(rs - rs0, tseq);
-		mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, mat, bw, opt->end_bonus, r->split_inv? opt->zdrop_inv : opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez);
+		mm_seq_rev(rs - rs0, junc);
+		mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, junc, mat, bw, opt->end_bonus, r->split_inv? opt->zdrop_inv : opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez);
 		if (ez->n_cigar > 0) {
 			mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			r->p->dp_score += ez->max;
@@ -698,6 +720,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			// perform alignment
 			qseq = &qseq0[rev][qs];
 			mm_idx_getseq(mi, rid, rs, re, tseq);
+			mm_idx_bed_junc(mi, rid, rs, re, junc);
 			if (is_sr) { // perform ungapped alignment
 				assert(qe - qs == re - rs);
 				ksw_reset_extz(ez);
@@ -707,11 +730,11 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				}
 				ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, 0, qe - qs);
 			} else { // perform normal gapped alignment
-				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
+				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
 			}
 			// test Z-drop and inversion Z-drop
 			if ((zdrop_code = mm_test_zdrop(km, opt, qseq, tseq, ez->n_cigar, ez->cigar, mat)) != 0)
-				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, zdrop_code == 2? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate
+				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, zdrop_code == 2? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate
 			// update CIGAR
 			if (ez->n_cigar > 0)
 				mm_append_cigar(r, ez->n_cigar, ez->cigar);
@@ -737,7 +760,8 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	if (!dropped && qe < qe0 && re < re0) { // right extension
 		qseq = &qseq0[rev][qe];
 		mm_idx_getseq(mi, rid, re, re0, tseq);
-		mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
+		mm_idx_bed_junc(mi, rid, re, re0, junc);
+		mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, junc, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
 		if (ez->n_cigar > 0) {
 			mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			r->p->dp_score += ez->max;
@@ -760,6 +784,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	}
 
 	kfree(km, tseq);
+	kfree(km, junc);
 }
 
 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)
@@ -793,7 +818,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	mm_seq_rev(tl, tseq);
 	if (score < opt->min_dp_max) goto end_align1_inv;
 	q_off = ql - (q_off + 1), t_off = tl - (t_off + 1);
-	mm_align_pair(km, opt, ql - q_off, qseq + q_off, tl - t_off, tseq + t_off, mat, (int)(opt->bw * 1.5), -1, opt->zdrop, KSW_EZ_EXTZ_ONLY, ez);
+	mm_align_pair(km, opt, ql - q_off, qseq + q_off, tl - t_off, tseq + t_off, 0, mat, (int)(opt->bw * 1.5), -1, opt->zdrop, KSW_EZ_EXTZ_ONLY, ez);
 	if (ez->n_cigar == 0) goto end_align1_inv; // should never be here
 	mm_append_cigar(r_inv, ez->n_cigar, ez->cigar);
 	r_inv->p->dp_score = ez->max;

chain.c

@@ -19,7 +19,7 @@ static inline int ilog2_32(uint32_t v)
 	return (t = v>>8) ? 8 + LogTable256[t] : LogTable256[v];
 }
 
-mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km)
+mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km)
 { // TODO: make sure this works when n has more than 32 bits
 	int32_t k, *f, *p, *t, *v, n_u, n_v;
 	int64_t i, j, st = 0;
@@ -28,6 +28,10 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 	mm128_t *b, *w;
 
 	if (_u) *_u = 0, *n_u_ = 0;
+	if (n == 0 || a == 0) {
+		kfree(km, a);
+		return 0;
+	}
 	f = (int32_t*)kmalloc(km, n * 4);
 	p = (int32_t*)kmalloc(km, n * 4);
 	t = (int32_t*)kmalloc(km, n * 4);
@@ -45,6 +49,7 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 		int32_t max_f = q_span, n_skip = 0, min_d;
 		int32_t sidi = (a[i].y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
 		while (st < i && ri > a[st].x + max_dist_x) ++st;
+		if (i - st > max_iter) st = i - max_iter;
 		for (j = i - 1; j >= st; --j) {
 			int64_t dr = ri - a[j].x;
 			int32_t dq = qi - (int32_t)a[j].y, dd, sc, log_dd;

cookbook.md

@@ -31,8 +31,8 @@ To acquire the data used in this cookbook and to install minimap2 and paftools,
 please follow the command lines below:
 ```sh
 # install minimap2 executables
-curl -L https://github.com/lh3/minimap2/releases/download/v2.15/minimap2-2.15_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.15_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
+curl -L https://github.com/lh3/minimap2/releases/download/v2.17/minimap2-2.17_x64-linux.tar.bz2 | tar jxf -
+cp minimap2-2.17_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
 export PATH="$PATH:"`pwd`                               # put the current directory on PATH
 # download example datasets
 curl -L https://github.com/lh3/minimap2/releases/download/v2.10/cookbook-data.tgz | tar zxf -

format.c

@@ -270,7 +270,7 @@ double mm_event_identity(const mm_reg1_t *r)
 		if (op == 1 || op == 2)
 			++n_gapo, n_gap += len;
 	}
-	return (double)r->mlen / (r->blen - r->p->n_ambi - n_gap + n_gapo);
+	return (double)r->mlen / (r->blen - n_gap + n_gapo);
 }
 
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
@@ -301,11 +301,12 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 	if (r->split) mm_sprintf_lite(s, "\tzd:i:%d", r->split);
 }
 
-void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag)
+void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag, int rep_len)
 {
 	s->l = 0;
 	if (r == 0) {
-		mm_sprintf_lite(s, "%s\t%d", t->name, t->l_seq);
+		mm_sprintf_lite(s, "%s\t%d\t0\t0\t*\t*\t0\t0\t0\t0\t0\t0", t->name, t->l_seq);
+		if (rep_len >= 0) mm_sprintf_lite(s, "\trl:i:%d", rep_len);
 		return;
 	}
 	mm_sprintf_lite(s, "%s\t%d\t%d\t%d\t%c\t", t->name, t->l_seq, r->qs, r->qe, "+-"[r->rev]);
@@ -315,6 +316,7 @@ void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m
 	mm_sprintf_lite(s, "\t%d\t%d", r->mlen, r->blen);
 	mm_sprintf_lite(s, "\t%d", r->mapq);
 	write_tags(s, r);
+	if (rep_len >= 0) mm_sprintf_lite(s, "\trl:i:%d", rep_len);
 	if (r->p && (opt_flag & MM_F_OUT_CG)) {
 		uint32_t k;
 		mm_sprintf_lite(s, "\tcg:Z:");
@@ -327,6 +329,11 @@ void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m
 		mm_sprintf_lite(s, "\t%s", t->comment);
 }
 
+void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag)
+{
+	mm_write_paf3(s, mi, t, r, km, opt_flag, -1);
+}
+
 static void sam_write_sq(kstring_t *s, char *seq, int l, int rev, int comp)
 {
 	extern unsigned char seq_comp_table[256];
@@ -368,6 +375,7 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
 			if (clip_len[1]) mm_sprintf_lite(s, ",%u", clip_len[1]<<4|clip_char);
 		} else {
 			int clip_char = (sam_flag&0x800) && !(opt_flag&MM_F_SOFTCLIP)? 'H' : 'S';
+			assert(clip_len[0] < qlen && clip_len[1] < qlen);
 			if (clip_len[0]) mm_sprintf_lite(s, "%d%c", clip_len[0], clip_char);
 			for (k = 0; k < r->p->n_cigar; ++k)
 				mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDNSHP=XB"[r->p->cigar[k]&0xf]);
@@ -376,7 +384,7 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
 	}
 }
 
-void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag)
+void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag, int rep_len)
 {
 	const int max_bam_cigar_op = 65535;
 	int flag, n_regs = n_regss[seg_idx], cigar_in_tag = 0;
@@ -452,17 +460,17 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 		int tlen = 0;
 		if (this_rid >= 0 && r_next) {
 			if (this_rid == r_next->rid) {
-				int this_pos5 = r && r->rev? r->re - 1 : this_pos;
-				int next_pos5 = r_next->rev? r_next->re - 1 : r_next->rs;
-				tlen = next_pos5 - this_pos5;
+				if (r) {
+					int this_pos5 = r->rev? r->re - 1 : this_pos;
+					int next_pos5 = r_next->rev? r_next->re - 1 : r_next->rs;
+					tlen = next_pos5 - this_pos5;
+				}
 				mm_sprintf_lite(s, "\t=\t");
 			} else mm_sprintf_lite(s, "\t%s\t", mi->seq[r_next->rid].name);
 			mm_sprintf_lite(s, "%d\t", r_next->rs + 1);
 		} else if (r_next) { // && this_rid < 0
 			mm_sprintf_lite(s, "\t%s\t%d\t", mi->seq[r_next->rid].name, r_next->rs + 1);
 		} else if (this_rid >= 0) { // && r_next == NULL
-			int this_pos5 = this_rev? r->re - 1 : this_pos; // this_rev is only true when r != NULL
-			tlen = this_pos - this_pos5; // next_pos5 will be this_pos
 			mm_sprintf_lite(s, "\t=\t%d\t", this_pos + 1); // next segment will take r's coordinate
 		} else mm_sprintf_lite(s, "\t*\t0\t"); // neither has coordinates
 		if (tlen > 0) ++tlen;
@@ -527,6 +535,7 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 		if (cigar_in_tag)
 			write_sam_cigar(s, flag, 1, t->l_seq, r, opt_flag);
 	}
+	if (rep_len >= 0) mm_sprintf_lite(s, "\trl:i:%d", rep_len);
 
 	if ((opt_flag & MM_F_COPY_COMMENT) && t->comment)
 		mm_sprintf_lite(s, "\t%s", t->comment);
@@ -534,6 +543,11 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 	s->s[s->l] = 0; // we always have room for an extra byte (see str_enlarge)
 }
 
+void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag)
+{
+	mm_write_sam3(s, mi, t, seg_idx, reg_idx, n_seg, n_regss, regss, km, opt_flag, -1);
+}
+
 void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int n_regs, const mm_reg1_t *regs)
 {
 	int i;

hit.c

@@ -449,6 +449,7 @@ void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int ma
 	int64_t sum_sc = 0;
 	float uniq_ratio;
 	int i;
+	if (n_regs == 0) return;
 	for (i = 0; i < n_regs; ++i)
 		if (regs[i].parent == regs[i].id)
 			sum_sc += regs[i].score;

index.c

@@ -31,6 +31,16 @@ typedef struct mm_idx_bucket_s {
 	void *h;     // hash table indexing _p_ and minimizers appearing once
 } mm_idx_bucket_t;
 
+typedef struct {
+	int32_t st, en, max; // max is not used for now
+	int32_t score:30, strand:2;
+} mm_idx_intv1_t;
+
+typedef struct mm_idx_intv_s {
+	int32_t n, m;
+	mm_idx_intv1_t *a;
+} mm_idx_intv_t;
+
 mm_idx_t *mm_idx_init(int w, int k, int b, int flag)
 {
 	mm_idx_t *mi;
@@ -55,6 +65,11 @@ void mm_idx_destroy(mm_idx_t *mi)
 			kh_destroy(idx, (idxhash_t*)mi->B[i].h);
 		}
 	}
+	if (mi->I) {
+		for (i = 0; i < mi->n_seq; ++i)
+			free(mi->I[i].a);
+		free(mi->I);
+	}
 	if (!mi->km) {
 		for (i = 0; i < mi->n_seq; ++i)
 			free(mi->seq[i].name);
@@ -585,3 +600,127 @@ int mm_idx_reader_eof(const mm_idx_reader_t *r) // TODO: in extremely rare cases
 {
 	return r->is_idx? (feof(r->fp.idx) || ftell(r->fp.idx) == r->idx_size) : mm_bseq_eof(r->fp.seq);
 }
+
+#include <ctype.h>
+#include <zlib.h>
+#include "ksort.h"
+#include "kseq.h"
+KSTREAM_DECLARE(gzFile, gzread)
+
+#define sort_key_bed(a) ((a).st)
+KRADIX_SORT_INIT(bed, mm_idx_intv1_t, sort_key_bed, 4)
+
+mm_idx_intv_t *mm_idx_read_bed(const mm_idx_t *mi, const char *fn, int read_junc)
+{
+	gzFile fp;
+	kstream_t *ks;
+	kstring_t str = {0,0,0};
+	mm_idx_intv_t *I;
+
+	fp = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
+	if (fp == 0) return 0;
+	I = (mm_idx_intv_t*)calloc(mi->n_seq, sizeof(*I));
+	ks = ks_init(fp);
+	while (ks_getuntil(ks, KS_SEP_LINE, &str, 0) >= 0) {
+		mm_idx_intv_t *r;
+		mm_idx_intv1_t t = {-1,-1,-1,-1,0};
+		char *p, *q, *bl, *bs;
+		int32_t i, id = -1, n_blk = 0;
+		for (p = q = str.s, i = 0;; ++p) {
+			if (*p == 0 || isspace(*p)) {
+				int32_t c = *p;
+				*p = 0;
+				if (i == 0) { // chr
+					id = mm_idx_name2id(mi, q);
+					if (id < 0) break; // unknown name; TODO: throw a warning
+				} else if (i == 1) { // start
+					t.st = atol(q); // TODO: watch out integer overflow!
+					if (t.st < 0) break;
+				} else if (i == 2) { // end
+					t.en = atol(q);
+					if (t.en < 0) break;
+				} else if (i == 4) { // BED score
+					t.score = atol(q);
+				} else if (i == 5) { // strand
+					t.strand = *q == '+'? 1 : *q == '-'? -1 : 0;
+				} else if (i == 9) {
+					if (!isdigit(*q)) break;
+					n_blk = atol(q);
+				} else if (i == 10) {
+					bl = q;
+				} else if (i == 11) {
+					bs = q;
+					break;
+				}
+				if (c == 0) break;
+				++i, q = p + 1;
+			}
+		}
+		if (id < 0 || t.st < 0 || t.st >= t.en) continue;
+		r = &I[id];
+		if (i >= 11 && read_junc) { // BED12
+			int32_t st, sz, en;
+			st = strtol(bs, &bs, 10); ++bs;
+			sz = strtol(bl, &bl, 10); ++bl;
+			en = t.st + st + sz;
+			for (i = 1; i < n_blk; ++i) {
+				mm_idx_intv1_t s = t;
+				if (r->n == r->m) {
+					r->m = r->m? r->m + (r->m>>1) : 16;
+					r->a = (mm_idx_intv1_t*)realloc(r->a, sizeof(*r->a) * r->m);
+				}
+				st = strtol(bs, &bs, 10); ++bs;
+				sz = strtol(bl, &bl, 10); ++bl;
+				s.st = en, s.en = t.st + st;
+				en = t.st + st + sz;
+				if (s.en > s.st) r->a[r->n++] = s;
+			}
+		} else {
+			if (r->n == r->m) {
+				r->m = r->m? r->m + (r->m>>1) : 16;
+				r->a = (mm_idx_intv1_t*)realloc(r->a, sizeof(*r->a) * r->m);
+			}
+			r->a[r->n++] = t;
+		}
+	}
+	free(str.s);
+	ks_destroy(ks);
+	gzclose(fp);
+	return I;
+}
+
+int mm_idx_bed_read(mm_idx_t *mi, const char *fn, int read_junc)
+{
+	int32_t i;
+	if (mi->h == 0) mm_idx_index_name(mi);
+	mi->I = mm_idx_read_bed(mi, fn, read_junc);
+	if (mi->I == 0) return -1;
+	for (i = 0; i < mi->n_seq; ++i) // TODO: eliminate redundant intervals
+		radix_sort_bed(mi->I[i].a, mi->I[i].a + mi->I[i].n);
+	return 0;
+}
+
+int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uint8_t *s)
+{
+	int32_t i, left, right;
+	mm_idx_intv_t *r;
+	memset(s, 0, en - st);
+	if (mi->I == 0 || ctg < 0 || ctg >= mi->n_seq) return -1;
+	r = &mi->I[ctg];
+	left = 0, right = r->n;
+	while (right > left) {
+		int32_t mid = left + ((right - left) >> 1);
+		if (r->a[mid].st >= st) right = mid;
+		else left = mid + 1;
+	}
+	for (i = left; i < r->n; ++i) {
+		if (st <= r->a[i].st && en >= r->a[i].en && r->a[i].strand != 0) {
+			if (r->a[i].strand > 0) {
+				s[r->a[i].st - st] |= 1, s[r->a[i].en - 1 - st] |= 2;
+			} else {
+				s[r->a[i].st - st] |= 8, s[r->a[i].en - 1 - st] |= 4;
+			}
+		}
+	}
+	return left;
+}

ketopt.h

@@ -73,13 +73,17 @@ static int ketopt(ketopt_t *s, int argc, char *argv[], int permute, const char *
 		}
 		s->opt = 0, opt = '?', s->pos = -1;
 		if (longopts) { /* parse long options */
-			int k, n_matches = 0;
-			const ko_longopt_t *o = 0;
+			int k, n_exact = 0, n_partial = 0;
+			const ko_longopt_t *o = 0, *o_exact = 0, *o_partial = 0;
 			for (j = 2; argv[s->i][j] != '\0' && argv[s->i][j] != '='; ++j) {} /* find the end of the option name */
 			for (k = 0; longopts[k].name != 0; ++k)
-				if (strncmp(&argv[s->i][2], longopts[k].name, j - 2) == 0)
-					++n_matches, o = &longopts[k];
-			if (n_matches == 1) {
+				if (strncmp(&argv[s->i][2], longopts[k].name, j - 2) == 0) {
+					if (longopts[k].name[j - 2] == 0) ++n_exact, o_exact = &longopts[k];
+					else ++n_partial, o_partial = &longopts[k];
+				}
+			if (n_exact > 1 || (n_exact == 0 && n_partial > 1)) return '?';
+			o = n_exact == 1? o_exact : n_partial == 1? o_partial : 0;
+			if (o) {
 				s->opt = opt = o->val, s->longidx = o - longopts;
 				if (argv[s->i][j] == '=') s->arg = &argv[s->i][j + 1];
 				if (o->has_arg == 1 && argv[s->i][j] == '\0') {

kseq.h

@@ -37,6 +37,14 @@
 #define KS_SEP_LINE  2 // line separator: "\n" (Unix) or "\r\n" (Windows)
 #define KS_SEP_MAX   2
 
+#ifndef klib_unused
+#if (defined __clang__ && __clang_major__ >= 3) || (defined __GNUC__ && __GNUC__ >= 3)
+#define klib_unused __attribute__ ((__unused__))
+#else
+#define klib_unused
+#endif
+#endif /* klib_unused */
+
 #define __KS_TYPE(type_t) \
 	typedef struct __kstream_t { \
 		int begin, end; \
@@ -64,7 +72,7 @@
 	}
 
 #define __KS_INLINED(__read) \
-	static inline int ks_getc(kstream_t *ks) \
+	static inline klib_unused int ks_getc(kstream_t *ks) \
 	{ \
 		if (ks->is_eof && ks->begin >= ks->end) return -1; \
 		if (ks->begin >= ks->end) { \

ksw2.h

@@ -61,7 +61,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
 
 void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
+				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez);
 
 void ksw_extf2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t mch, int8_t mis, int8_t e, int w, int xdrop, ksw_extz_t *ez);
 

ksw2_dispatch.c

@@ -80,17 +80,17 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 }
 
 void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez)
 {
 	extern void ksw_exts2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez);
 	extern void ksw_exts2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez);
 	if (ksw_simd < 0) ksw_simd = x86_simd();
 	if (ksw_simd & SIMD_SSE4_1)
-		ksw_exts2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
+		ksw_exts2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, junc_bonus, flag, junc, ez);
 	else if (ksw_simd & SIMD_SSE2)
-		ksw_exts2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
+		ksw_exts2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, junc_bonus, flag, junc, ez);
 	else abort();
 }
 #endif

ksw2_exts2_sse.c

@@ -17,14 +17,14 @@
 #ifdef KSW_CPU_DISPATCH
 #ifdef __SSE4_1__
 void ksw_exts2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez)
 #else
 void ksw_exts2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez)
 #endif
 #else
 void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez)
 #endif // ~KSW_CPU_DISPATCH
 {
 #define __dp_code_block1 \
@@ -113,20 +113,53 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	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 - 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 can_type = 0;
-			if ((flag & KSW_EZ_SPLICE_FOR) && target[t-1] == 0 && target[t] == 2) can_type = 1; // ...yAG
-			if ((flag & KSW_EZ_SPLICE_REV) && target[t-1] == 0 && target[t] == 1) can_type = 1; // ...yAC
-			if (can_type && (target[t-2] == 1 || target[t-2] == 3)) can_type = 2;
-			if (can_type) ((int8_t*)acceptor)[t] = can_type == 2? 0 : semi_cost;
+		if (!(flag & KSW_EZ_REV_CIGAR)) {
+			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;
+			}
+			if (junc)
+				for (t = 0; t < tlen - 1; ++t)
+					if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t+1]&1)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t+1]&8)))
+						((int8_t*)donor)[t] += junc_bonus;
+			for (t = 2; t < tlen; ++t) {
+				int can_type = 0;
+				if ((flag & KSW_EZ_SPLICE_FOR) && target[t-1] == 0 && target[t] == 2) can_type = 1; // ...yAG
+				if ((flag & KSW_EZ_SPLICE_REV) && target[t-1] == 0 && target[t] == 1) can_type = 1; // ...yAC
+				if (can_type && (target[t-2] == 1 || target[t-2] == 3)) can_type = 2;
+				if (can_type) ((int8_t*)acceptor)[t] = can_type == 2? 0 : semi_cost;
+			}
+			if (junc)
+				for (t = 0; t < tlen; ++t)
+					if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t]&2)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t]&4)))
+						((int8_t*)acceptor)[t] += junc_bonus;
+		} else {
+			for (t = 0; t < tlen - 4; ++t) {
+				int can_type = 0; // type of canonical site: 0=none, 1=GT/AG only, 2=GTr/yAG
+				if ((flag & KSW_EZ_SPLICE_FOR) && target[t+1] == 2 && target[t+2] == 0) can_type = 1; // GAy...
+				if ((flag & KSW_EZ_SPLICE_REV) && target[t+1] == 1 && target[t+2] == 0) can_type = 1; // CAy...
+				if (can_type && (target[t+3] == 1 || target[t+3] == 3)) can_type = 2;
+				if (can_type) ((int8_t*)donor)[t] = can_type == 2? 0 : semi_cost;
+			}
+			if (junc)
+				for (t = 0; t < tlen - 1; ++t)
+					if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t+1]&2)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t+1]&4)))
+						((int8_t*)donor)[t] += junc_bonus;
+			for (t = 2; t < tlen; ++t) {
+				int can_type = 0;
+				if ((flag & KSW_EZ_SPLICE_FOR) && target[t-1] == 3 && target[t] == 2) can_type = 1; // ...rTG
+				if ((flag & KSW_EZ_SPLICE_REV) && target[t-1] == 3 && target[t] == 1) can_type = 1; // ...rTC
+				if (can_type && (target[t-2] == 0 || target[t-2] == 2)) can_type = 2;
+				if (can_type) ((int8_t*)acceptor)[t] = can_type == 2? 0 : semi_cost;
+			}
+			if (junc)
+				for (t = 0; t < tlen; ++t)
+					if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t]&1)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t]&8)))
+						((int8_t*)acceptor)[t] += junc_bonus;
 		}
 	}
 

main.c

@@ -6,7 +6,7 @@
 #include "mmpriv.h"
 #include "ketopt.h"
 
-#define MM_VERSION "2.15-r905"
+#define MM_VERSION "2.17-r941"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -62,6 +62,10 @@ static ko_longopt_t long_options[] = {
 	{ "hard-mask-level",ko_no_argument,       336 },
 	{ "cap-sw-mem",     ko_required_argument, 337 },
 	{ "max-qlen",       ko_required_argument, 338 },
+	{ "max-chain-iter", ko_required_argument, 339 },
+	{ "junc-bed",       ko_required_argument, 340 },
+	{ "junc-bonus",     ko_required_argument, 341 },
+	{ "sam-hit-only",   ko_no_argument,       342 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
 	{ "version",        ko_no_argument,       'V' },
@@ -99,12 +103,12 @@ static inline void yes_or_no(mm_mapopt_t *opt, int flag, int long_idx, const cha
 
 int main(int argc, char *argv[])
 {
-	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:yYP";
+	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:yYPo:";
 	ketopt_t o = KETOPT_INIT;
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
 	int i, c, n_threads = 3, n_parts, old_best_n = -1;
-	char *fnw = 0, *rg = 0, *s;
+	char *fnw = 0, *rg = 0, *junc_bed = 0, *s;
 	FILE *fp_help = stderr;
 	mm_idx_reader_t *idx_rdr;
 	mm_idx_t *mi;
@@ -165,12 +169,21 @@ int main(int argc, char *argv[])
 		else if (c == 'R') rg = o.arg;
 		else if (c == 'h') fp_help = stdout;
 		else if (c == '2') opt.flag |= MM_F_2_IO_THREADS;
+		else if (c == 'o') {
+			if (strcmp(o.arg, "-") != 0) {
+				if (freopen(o.arg, "wb", stdout) == NULL) {
+					fprintf(stderr, "[ERROR]\033[1;31m failed to write the output to file '%s'\033[0m\n", o.arg);
+					exit(1);
+				}
+			}
+		}
 		else if (c == 300) ipt.bucket_bits = atoi(o.arg); // --bucket-bits
 		else if (c == 302) opt.seed = atoi(o.arg); // --seed
 		else if (c == 303) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
 		else if (c == 304) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
 		else if (c == 306) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
 		else if (c == 307) opt.max_chain_skip = atoi(o.arg); // --max-chain-skip
+		else if (c == 339) opt.max_chain_iter = atoi(o.arg); // --max-chain-iter
 		else if (c == 308) opt.min_ksw_len = atoi(o.arg); // --min-dp-len
 		else if (c == 309) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
 		else if (c == 310) opt.flag |= MM_F_SPLICE; // --splice
@@ -194,6 +207,8 @@ int main(int argc, char *argv[])
 		else if (c == 336) opt.flag |= MM_F_HARD_MLEVEL; // --hard-mask-level
 		else if (c == 337) opt.max_sw_mat = mm_parse_num(o.arg); // --cap-sw-mat
 		else if (c == 338) opt.max_qlen = mm_parse_num(o.arg); // --max-qlen
+		else if (c == 340) junc_bed = o.arg; // --junc-bed
+		else if (c == 342) opt.flag |= MM_F_SAM_HIT_ONLY; // --sam-hit-only
 		else if (c == 314) { // --frag
 			yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1);
 		} else if (c == 315) { // --secondary
@@ -268,7 +283,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "  Indexing:\n");
 		fprintf(fp_help, "    -H           use homopolymer-compressed k-mer (preferrable for PacBio)\n");
 		fprintf(fp_help, "    -k INT       k-mer size (no larger than 28) [%d]\n", ipt.k);
-		fprintf(fp_help, "    -w INT       minizer window size [%d]\n", ipt.w);
+		fprintf(fp_help, "    -w INT       minimizer window size [%d]\n", ipt.w);
 		fprintf(fp_help, "    -I NUM       split index for every ~NUM input bases [4G]\n");
 		fprintf(fp_help, "    -d FILE      dump index to FILE []\n");
 		fprintf(fp_help, "  Mapping:\n");
@@ -293,7 +308,7 @@ int main(int argc, char *argv[])
 		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, "  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");
+		fprintf(fp_help, "    -o FILE      output alignments to FILE [stdout]\n");
 		fprintf(fp_help, "    -L           write CIGAR with >65535 ops at the CG tag\n");
 		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");
@@ -353,6 +368,7 @@ int main(int argc, char *argv[])
 					__func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), mi->n_seq);
 		if (argc != o.ind + 1) mm_mapopt_update(&opt, mi);
 		if (mm_verbose >= 3) mm_idx_stat(mi);
+		if (junc_bed) mm_idx_bed_read(mi, junc_bed, 1);
 		if (!(opt.flag & MM_F_FRAG_MODE)) {
 			for (i = o.ind + 1; i < argc; ++i)
 				mm_map_file(mi, argv[i], &opt, n_threads);

map.c

@@ -312,7 +312,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		if (max_chain_gap_ref < opt->max_gap) max_chain_gap_ref = opt->max_gap;
 	} else max_chain_gap_ref = opt->max_gap;
 
-	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);
+	a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 
 	if (opt->max_occ > opt->mid_occ && rep_len > 0) {
 		int rechain = 0;
@@ -334,7 +334,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 			kfree(b->km, mini_pos);
 			if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 			else a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
-			a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
+			a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 		}
 	}
 	b->frag_gap = max_chain_gap_ref;
@@ -584,16 +584,16 @@ static void *worker_pipeline(void *shared, int step, void *in)
 						if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent)
 							continue;
 						if (p->opt->flag & MM_F_OUT_SAM)
-							mm_write_sam2(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
+							mm_write_sam3(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag, s->rep_len[i]);
 						else
-							mm_write_paf(&p->str, mi, t, r, km, p->opt->flag);
+							mm_write_paf3(&p->str, mi, t, r, km, p->opt->flag, s->rep_len[i]);
 						mm_err_puts(p->str.s);
 					}
-				} else if (p->opt->flag & (MM_F_OUT_SAM|MM_F_PAF_NO_HIT)) { // output an empty hit, if requested
+				} else if ((p->opt->flag & MM_F_PAF_NO_HIT) || ((p->opt->flag & MM_F_OUT_SAM) && !(p->opt->flag & MM_F_SAM_HIT_ONLY))) { // output an empty hit, if requested
 					if (p->opt->flag & MM_F_OUT_SAM)
-						mm_write_sam2(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
+						mm_write_sam3(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag, s->rep_len[i]);
 					else
-						mm_write_paf(&p->str, mi, t, 0, 0, p->opt->flag);
+						mm_write_paf3(&p->str, mi, t, 0, 0, p->opt->flag, s->rep_len[i]);
 					mm_err_puts(p->str.s);
 				}
 			}

minimap.h

@@ -35,6 +35,7 @@
 #define MM_F_PAF_NO_HIT    0x8000000 // output unmapped reads to PAF
 #define MM_F_NO_END_FLT    0x10000000
 #define MM_F_HARD_MLEVEL   0x20000000
+#define MM_F_SAM_HIT_ONLY  0x40000000
 
 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2
@@ -66,6 +67,7 @@ typedef struct {
 	mm_idx_seq_t *seq;         // sequence name, length and offset
 	uint32_t *S;               // 4-bit packed sequence
 	struct mm_idx_bucket_s *B; // index (hidden)
+	struct mm_idx_intv_s *I;   // intervals (hidden)
 	void *km, *h;
 } mm_idx_t;
 
@@ -103,16 +105,16 @@ typedef struct {
 } mm_idxopt_t;
 
 typedef struct {
+	int64_t flag;    // see MM_F_* macros
 	int seed;
 	int sdust_thres; // score threshold for SDUST; 0 to disable
-	int flag;        // see MM_F_* macros
 
 	int max_qlen;    // max query length
 
 	int bw;          // bandwidth
 	int max_gap, max_gap_ref; // break a chain if there are no minimizers in a max_gap window
 	int max_frag_len;
-	int max_chain_skip;
+	int max_chain_skip, max_chain_iter;
 	int min_cnt;         // min number of minimizers on each chain
 	int min_chain_score; // min chaining score
 
@@ -127,6 +129,7 @@ typedef struct {
 	int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties
 	int sc_ambi; // score when one or both bases are "N"
 	int noncan;      // cost of non-canonical splicing sites
+	int junc_bonus;
 	int zdrop, zdrop_inv;   // break alignment if alignment score drops too fast along the diagonal
 	int end_bonus;
 	int min_dp_max;  // drop an alignment if the score of the max scoring segment is below this threshold
@@ -365,6 +368,9 @@ int mm_idx_index_name(mm_idx_t *mi);
 int mm_idx_name2id(const mm_idx_t *mi, const char *name);
 int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
 
+int mm_idx_bed_read(mm_idx_t *mi, const char *fn, int read_junc);
+int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uint8_t *s);
+
 // 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 flag, int n_threads);

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "10 January 2019" "minimap2-2.15 (r905)" "Bioinformatics tools"
+.TH minimap2 1 "4 May 2019" "minimap2-2.17 (r941)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -232,13 +232,19 @@ Honor option
 and disable a heurstic to save unmapped subsequences.
 .TP
 .BI --max-chain-skip \ INT
-A heuristics that stops chaining early [50]. Minimap2 uses dynamic programming
+A heuristics that stops chaining early [25]. Minimap2 uses dynamic programming
 for chaining. The time complexity is quadratic in the number of seeds. This
 option makes minimap2 exits the inner loop if it repeatedly sees seeds already
 on chains. Set
 .I INT
 to a large number to switch off this heurstics.
 .TP
+.BI --max-chain-iter \ INT
+Check up to
+.I INT
+partial chains during chaining [5000]. This is a heuristic to avoid quadratic
+time complexity in the worst case.
+.TP
 .B --no-long-join
 Disable the long gap patching heuristic. When this option is applied, the
 maximum alignment gap is mostly controlled by
@@ -358,6 +364,17 @@ on SIRV data, please add
 .B --splice-flank=no
 to the command line.
 .TP
+.BR --junc-bed \ FILE
+Gene annotations in the BED12 format (aka 12-column BED), or intron positions
+in 5-column BED. With this option, minimap2 prefers splicing in annotations.
+BED12 file can be converted from GTF/GFF3 with `paftools.js gff2bed anno.gtf'
+[].
+.TP
+.BR --junc-bonus \ INT
+Score bonus for a splice donor or acceptor found in annotation (effective with
+.BR --junc-bed )
+[0].
+.TP
 .BI --end-seed-pen \ INT
 Drop a terminal anchor if
 .IR s <log( g )+ INT ,
@@ -384,6 +401,11 @@ Set 0 to disable [0].
 Generate CIGAR and output alignments in the SAM format. Minimap2 outputs in PAF
 by default.
 .TP
+.BI -o \ FILE
+Output alignments to
+.I FILE
+[stdout].
+.TP
 .B -Q
 Ignore base quality in the input file.
 .TP
@@ -463,7 +485,12 @@ Filter out query sequences longer than
 .IR NUM .
 .TP
 .B --paf-no-hit
-In PAF, output query name and length for an unmapped sequence.
+In PAF, output unmapped queries; the strand and the reference name fields are
+set to `*'. Warning: some paftools.js commands may not work with such output
+for the moment.
+.TP
+.B --sam-hit-only
+In SAM, don't output unmapped reads.
 .TP
 .B --version
 Print version number to stdout
@@ -494,7 +521,7 @@ is determined by the sequencing error mode.
 .B asm5
 Long assembly to reference mapping
 .RB ( -k19
-.B -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200
+.B -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200 -N50
 .BR --min-occ-floor=100 ).
 Typically, the alignment will not extend to regions with 5% or higher sequence
 divergence. Only use this preset if the average divergence is far below 5%.
@@ -502,14 +529,14 @@ divergence. Only use this preset if the average divergence is far below 5%.
 .B asm10
 Long assembly to reference mapping
 .RB ( -k19
-.B -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200
+.B -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200 -N50
 .BR --min-occ-floor=100 ).
 Up to 10% sequence divergence.
 .TP
 .B asm20
 Long assembly to reference mapping
 .RB ( -k19
-.B -w10 -A1 -B4 -O6,26 -E2,1 -s200 -z200
+.B -w10 -A1 -B4 -O6,26 -E2,1 -s200 -z200 -N50
 .BR --min-occ-floor=100 ).
 Up to 20% sequence divergence.
 .TP
@@ -531,7 +558,7 @@ 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 -C9 -z200 -ub
+.B -w5 --splice -g2000 -G200k -A1 -B2 -O2,32 -E1,0 -C9 -z200 -ub --junc-bonus=9
 .BR --splice-flank=yes ).
 In the splice mode, 1) long deletions are taken as introns and represented as
 the
@@ -541,6 +568,12 @@ costs are different during chaining; 4) the computation of the
 .RB ` ms '
 tag ignores introns to demote hits to pseudogenes.
 .TP
+.B splice:hq
+Long-read splice alignment for PacBio CCS reads
+.RB ( -xsplice
+.B -C5 -O6,24
+.BR -B4 ).
+.TP
 .B sr
 Short single-end reads without splicing
 .RB ( -k21
@@ -612,6 +645,7 @@ cg	Z	CIGAR string (only in PAF)
 cs	Z	Difference string
 dv	f	Approximate per-base sequence divergence
 de	f	Gap-compressed per-base sequence divergence
+rl	i	Length of query regions harboring repetitive seeds
 .TE
 
 .PP

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = '2.15-r905';
+var paftools_version = '2.17-r941';
 
 /*****************************
  ***** Library functions *****
@@ -433,6 +433,7 @@ function paf_call(args)
 	while (file.readline(buf) >= 0) {
 		var line = buf.toString();
 		var m, t = line.split("\t", 12);
+		if (t.length < 12 || t[5] == '*') continue; // unmapped
 		for (var i = 6; i <= 11; ++i)
 			t[i] = parseInt(t[i]);
 		if (t[10] < min_cov_len || t[11] < min_mapq) continue;
@@ -680,13 +681,12 @@ function paf_asmstat(args)
 			var t = line.split("\t");
 			t[1] = parseInt(t[1]);
 			if (t[1] < min_query_len) continue;
-			if (t.length >= 2) {
-				query[t[0]] = t[1];
-				if (qinfo[t[0]] == null) qinfo[t[0]] = {};
-				qinfo[t[0]].len = t[1];
-				qinfo[t[0]].bp = [];
-			}
-			if (t.length < 9) continue;
+			if (t.length < 2) continue;
+			query[t[0]] = t[1];
+			if (qinfo[t[0]] == null) qinfo[t[0]] = {};
+			qinfo[t[0]].len = t[1];
+			qinfo[t[0]].bp = [];
+			if (t.length < 9 || t[5] == "*") continue;
 			if (!/\ttp:A:[PI]/.test(line)) continue;
 			if ((m = /\tcg:Z:(\S+)/.exec(line)) == null) continue;
 			var cigar = m[1];
@@ -967,7 +967,9 @@ function paf_stat(args)
 			var t = line.split("\t", 12);
 			var m, rs, cigar = null, is_pri = false, is_sam = false, is_rev = false, tname = null;
 			var atlen = null, aqlen, qs, qe, mapq, ori_qlen;
-			if (t[4] == '+' || t[4] == '-') { // PAF
+			if (t.length < 2) continue;
+			if (t[4] == '+' || t[4] == '-' || t[4] == '*') { // PAF
+				if (t[4] == '*') continue; // unmapped
 				if (!/\ts2:i:\d+/.test(line)) {
 					++n_2nd;
 					continue;
@@ -1467,15 +1469,21 @@ function paf_view(args)
 
 function paf_gff2bed(args)
 {
-	var c, fn_ucsc_fai = null, is_short = false, keep_gff = false;
-	while ((c = getopt(args, "u:sg")) != null) {
+	var c, fn_ucsc_fai = null, is_short = false, keep_gff = false, print_junc = false;
+	while ((c = getopt(args, "u:sgj")) != null) {
 		if (c == 'u') fn_ucsc_fai = getopt.arg;
 		else if (c == 's') is_short = true;
 		else if (c == 'g') keep_gff = true;
+		else if (c == 'j') print_junc = true;
 	}
 
 	if (getopt.ind == args.length) {
-		print("Usage: paftools.js gff2bed [-g] [-u ucsc-genome.fa.fai] <in.gff>");
+		print("Usage: paftools.js gff2bed [options] <in.gff>");
+		print("Options:");
+		print("  -j       Output junction BED");
+		print("  -s       Print names in the short form");
+		print("  -u FILE  hg38.fa.fai for chr name conversion");
+		print("  -g       Output GFF (used with -u)");
 		exit(1);
 	}
 
@@ -1507,11 +1515,16 @@ function paf_gff2bed(args)
 		'misc_RNA':'0,192,0'
 	};
 
-	function print_bed12(exons, cds_st, cds_en, is_short)
+	function print_bed12(exons, cds_st, cds_en, is_short, print_junc)
 	{
 		if (exons.length == 0) return;
 		var name = is_short? exons[0][7] + "|" + exons[0][5] : exons[0].slice(4, 7).join("|");
 		var a = exons.sort(function(a,b) {return a[1]-b[1]});
+		if (print_junc) {
+		for (var i = 1; i < a.length; ++i)
+			print(a[i][0], a[i-1][2], a[i][1], name, 1000, a[i][3]);
+			return;
+		}
 		var sizes = [], starts = [], st, en;
 		st = a[0][1];
 		en = a[a.length - 1][2];
@@ -1564,7 +1577,7 @@ function paf_gff2bed(args)
 		if (type == "" && biotype != "") type = biotype;
 		if (id == null) throw Error("No transcript_id");
 		if (id != last_id) {
-			print_bed12(exons, cds_st, cds_en, is_short);
+			print_bed12(exons, cds_st, cds_en, is_short, print_junc);
 			exons = [], cds_st = 1<<30, cds_en = 0;
 			last_id = id;
 		}
@@ -1582,7 +1595,7 @@ function paf_gff2bed(args)
 		}
 	}
 	if (last_id != null)
-		print_bed12(exons, cds_st, cds_en, is_short);
+		print_bed12(exons, cds_st, cds_en, is_short, print_junc);
 
 	file.close();
 	buf.destroy();
@@ -1590,11 +1603,16 @@ function paf_gff2bed(args)
 
 function paf_sam2paf(args)
 {
-	var c, pri_only = false, use_eq = false;
-	while ((c = getopt(args, "p")) != null)
+	var c, pri_only = false, long_cs = false;
+	while ((c = getopt(args, "pL")) != null) {
 		if (c == 'p') pri_only = true;
+		else if (c == 'L') long_cs = true;
+	}
 	if (args.length == getopt.ind) {
-		print("Usage: paftools.js sam2paf [-p] <in.sam>");
+		print("Usage: paftools.js sam2paf [options] <in.sam>");
+		print("Options:");
+		print("  -p      convert primary or supplementary alignments only");
+		print("  -L      output the cs tag in the long form");
 		exit(1);
 	}
 
@@ -1623,13 +1641,14 @@ function paf_sam2paf(args)
 		var tlen = ctg_len[t[2]];
 		if (tlen == null) throw Error("at line " + lineno + ": can't find the length of contig " + t[2]);
 		// find tags
-		var nn = 0, NM = null, MD = null, md_list = [];
+		var nn = 0, NM = null, MD = null, cs_str = null, md_list = [];
 		while ((m = re_tag.exec(line)) != null) {
 			if (m[1] == "NM:i") NM = parseInt(m[2]);
 			else if (m[1] == "nn:i") nn = parseInt(m[2]);
 			else if (m[1] == "MD:Z") MD = m[2];
+			else if (m[1] == "cs:Z") cs_str = m[2];
 		}
-		if (t[9] == '*') MD = null;
+		if (t[9] == '*') MD = cs_str = null;
 		// infer various lengths from CIGAR
 		var clip = [0, 0], soft_clip = 0, I = [0, 0], D = [0, 0], M = 0, N = 0, mm = 0, have_M = false, have_ext = false, cigar = [];
 		while ((m = re.exec(t[5])) != null) {
@@ -1665,8 +1684,8 @@ function paf_sam2paf(args)
 		}
 		// parse MD
 		var cs = [];
-		if (MD != null) {
-			var k = 0, cx = 0, cy = 0, mx = 0, my = 0;
+		if (MD != null && cs_str == null && t[9] != "*") {
+			var k = 0, cx = 0, cy = 0, mx = 0, my = 0; // cx: cigar ref position; cy: cigar query; mx: MD ref; my: MD query
 			while ((m = re_MD.exec(MD)) != null) {
 				if (m[2] != null) { // deletion from the reference
 					var len = m[2].length - 1;
@@ -1680,13 +1699,15 @@ function paf_sam2paf(args)
 							if (my + ml < cy + cl) {
 								if (ml > 0) {
 									if (m[3] != null) cs.push('*', m[3], t[9][my]);
+									else if (long_cs) cs.push('=', t[9].substr(my, ml));
 									else cs.push(':', ml);
 								}
 								mx += ml, my += ml, ml = 0;
 								break;
 							} else {
 								var dl = cy + cl - my;
-								cs.push(':', dl);
+								if (long_cs) cs.push('=', t[9].substr(my, dl));
+								else cs.push(':', dl);
 								cx += cl, cy += cl, ++k;
 								mx += dl, my += dl, ml -= dl;
 							}
@@ -1731,7 +1752,8 @@ function paf_sam2paf(args)
 		var tags = ["tp:A:" + type];
 		if (NM != null) tags.push("mm:i:"+mm);
 		tags.push("gn:i:"+(I[1]+D[1]), "go:i:"+(I[0]+D[0]), "cg:Z:" + t[5].replace(/\d+[SH]/g, ''));
-		if (cs.length > 0) tags.push("cs:Z:" + cs.join(""));
+		if (cs_str != null) tags.push("cs:Z:" + cs_str);
+		else if (cs.length > 0) tags.push("cs:Z:" + cs.join(""));
 		// print out
 		var a = [qname, qlen, qs, qe, flag&16? '-' : '+', t[2], tlen, ts, te, mlen, blen, t[4]];
 		print(a.join("\t"), tags.join("\t"));

mmpriv.h

@@ -61,13 +61,15 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i
 
 void mm_write_sam_hdr(const mm_idx_t *mi, const char *rg, const char *ver, int argc, char *argv[]);
 void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag);
+void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag, int rep_len);
 void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int n_regs, const mm_reg1_t *regs);
 void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regs, const mm_reg1_t *const* regs, void *km, int opt_flag);
+void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag, int rep_len);
 
 void mm_idxopt_init(mm_idxopt_t *opt);
 const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f);
-mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
+mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
 mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
 
 mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a);

options.c

@@ -23,6 +23,7 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->max_gap = 5000;
 	opt->max_gap_ref = -1;
 	opt->max_chain_skip = 25;
+	opt->max_chain_iter = 5000;
 
 	opt->mask_level = 0.5f;
 	opt->pri_ratio = 0.8f;
@@ -119,13 +120,16 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		mo->mid_occ = 1000;
 		mo->max_occ = 5000;
 		mo->mini_batch_size = 50000000;
-	} else if (strcmp(preset, "splice") == 0 || strcmp(preset, "cdna") == 0) {
+	} else if (strncmp(preset, "splice", 6) == 0 || strcmp(preset, "cdna") == 0) {
 		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 = 9;
+		mo->junc_bonus = 9;
 		mo->zdrop = 200, mo->zdrop_inv = 100; // because mo->a is halved
+		if (strcmp(preset, "splice:hq") == 0)
+			mo->junc_bonus = 5, mo->b = 4, mo->q = 6, mo->q2 = 24;
 	} else return -1;
 	return 0;
 }

python/README.rst

@@ -114,6 +114,12 @@ This method retrieves a (sub)sequence from the index and returns it as a Python
 string. :code:`None` is returned if :code:`name` is not present in the index or
 the start/end coordinates are invalid.
 
+.. code:: python
+
+	mappy.Aligner.seq_names
+
+This property gives the array of sequence names in the index.
+
 Class mappy.Alignment
 ~~~~~~~~~~~~~~~~~~~~~
 

python/cmappy.pxd

@@ -10,13 +10,14 @@ cdef extern from "minimap.h":
 		uint64_t batch_size
 
 	ctypedef struct mm_mapopt_t:
+		int64_t flag
 		int seed
 		int sdust_thres
-		int flag
+		int max_qlen
 		int bw
 		int max_gap, max_gap_ref
 		int max_frag_len
-		int max_chain_skip
+		int max_chain_skip, max_chain_iter
 		int min_cnt
 		int min_chain_score
 		float mask_level
@@ -24,10 +25,11 @@ cdef extern from "minimap.h":
 		int best_n
 		int max_join_long, max_join_short
 		int min_join_flank_sc
-		float min_join_flank_ratio;
+		float min_join_flank_ratio
 		int a, b, q, e, q2, e2
 		int sc_ambi
 		int noncan
+		int junc_bonus
 		int zdrop, zdrop_inv
 		int end_bonus
 		int min_dp_max

python/mappy.pyx

@@ -3,7 +3,7 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
-__version__ = '2.15'
+__version__ = '2.17'
 
 cmappy.mm_reset_timer()
 
@@ -142,7 +142,7 @@ cdef class Aligner:
 			if fn_idx_out is None:
 				r = cmappy.mm_idx_reader_open(str.encode(fn_idx_in), &self.idx_opt, NULL)
 			else:
-				r = cmappy.mm_idx_reader_open(str.encode(fn_idx_in), &self.idx_opt, fn_idx_out)
+				r = cmappy.mm_idx_reader_open(str.encode(fn_idx_in), &self.idx_opt, str.encode(fn_idx_out))
 			if r is not NULL:
 				self._idx = cmappy.mm_idx_reader_read(r, n_threads) # NB: ONLY read the first part
 				cmappy.mm_idx_reader_close(r)
@@ -221,6 +221,16 @@ cdef class Aligner:
 	@property
 	def n_seq(self): return self._idx.n_seq
 
+	@property
+	def seq_names(self):
+		cdef char *p
+		sn = []
+		for i in range(self._idx.n_seq):
+			p = self._idx.seq[i].name
+			s = p if isinstance(p, str) else p.decode()
+			sn.append(s)
+		return sn
+
 def fastx_read(fn, read_comment=False):
 	cdef cmappy.kseq_t *ks
 	ks = cmappy.mm_fastx_open(str.encode(fn))

setup.py

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

splitidx.c

@@ -11,8 +11,11 @@ FILE *mm_split_init(const char *prefix, const mm_idx_t *mi)
 	uint32_t i, k = mi->k;
 	fn = (char*)calloc(strlen(prefix) + 10, 1);
 	sprintf(fn, "%s.%.4d.tmp", prefix, mi->index);
-	fp = fopen(fn, "wb");
-	assert(fp);
+	if ((fp = fopen(fn, "wb")) == NULL) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m failed to write to temporary file '%s'\033[0m\n", fn);
+		exit(1);
+	}
 	mm_err_fwrite(&k, 4, 1, fp);
 	mm_err_fwrite(&mi->n_seq, 4, 1, fp);
 	for (i = 0; i < mi->n_seq; ++i) {