Release minimap2-2.16 (r922)

Resolves #324

LICENSE.txt

@@ -1,6 +1,7 @@
 The MIT License
 
-Copyright (c) 2017 Broad Institute, Inc.
+Copyright (c) 2018-     Dana-Farber Cancer Institute
+              2017-2018 Broad Institute, Inc.
 
 Permission is hereby granted, free of charge, to any person obtaining
 a copy of this software and associated documentation files (the

NEWS.md

@@ -1,3 +1,114 @@
+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)
+-----------------------------------
+
+Changes to minimap2:
+
+ * Fixed a rare segmentation fault when option -H is in use (#307). This may
+   happen when there are very long homopolymers towards the 5'-end of a read.
+
+ * Fixed wrong CIGARs when option --eqx is used (#266).
+
+ * Fixed a typo in the base encoding table (#264). This should have no
+   practical effect.
+
+ * Fixed a typo in the example code (#265).
+
+ * Improved the C++ compatibility by removing "register" (#261). However,
+   minimap2 still can't be compiled in the pedantic C++ mode (#306).
+
+ * Output a new "de" tag for gap-compressed sequence divergence.
+
+Changes to paftools.js:
+
+ * Added "asmgene" to evaluate the completeness of an assembly by measuring the
+   uniquely mapped single-copy genes. This command learns the idea of BUSCO.
+
+ * Added "vcfpair" to call a phased VCF from phased whole-genome assemblies. An
+   earlier version of this script is used to produce the ground truth for the
+   syndip benchmark [PMID:30013044].
+
+This release produces identical alignment coordinates and CIGARs in comparison
+to v2.14. Users are advised to upgrade due to the several bug fixes.
+
+(2.15: 10 Janurary 2019, r905)
+
+
+
+Release 2.14-r883 (5 November 2018)
+-----------------------------------
+
+Notable changes:
+
+ * Fixed two minor bugs caused by typos (#254 and #266).
+
+ * Fixed a bug that made minimap2 abort when --eqx was used together with --MD
+   or --cs (#257).
+
+ * Added --cap-sw-mem to cap the size of DP matrices (#259). Base alignment may
+   take a lot of memory in the splicing mode. This may lead to issues when we
+   run minimap2 on a cluster with a hard memory limit. The new option avoids
+   unlimited memory usage at the cost of missing a few long introns.
+
+ * Conforming to C99 and C11 when possible (#261).
+
+ * Warn about malformatted FASTA or FASTQ (#252 and #255).
+
+This release occasionally produces base alignments different from v2.13. The
+overall alignment accuracy remain similar.
+
+(2.14: 5 November 2018, r883)
+
+
+
 Release 2.13-r850 (11 October 2018)
 -----------------------------------
 

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
@@ -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.13/minimap2-2.13_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.13_x64-linux/minimap2
+curl -L https://github.com/lh3/minimap2/releases/download/v2.16/minimap2-2.16_x64-linux.tar.bz2 | tar -jxvf -
+./minimap2-2.16_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
@@ -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
 
@@ -355,6 +355,8 @@ mappy` or [from BioConda][mappyconda] via `conda install -c bioconda mappy`.
   billion bases or longer (2,147,483,647 to be exact). The total length of all
   sequences can well exceed this threshold.
 
+* Minimap2 often misses small exons.
+
 
 
 [paf]: https://github.com/lh3/miniasm/blob/master/PAF.md

align.c

@@ -147,78 +147,6 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 	}
 }
 
-static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e)
-{
-	uint32_t k, l;
-	int32_t s = 0, max = 0, qshift, tshift, toff = 0, qoff = 0;
-	mm_extra_t *p = r->p;
-	if (p == 0) return;
-	mm_fix_cigar(r, qseq, tseq, &qshift, &tshift);
-	qseq += qshift, tseq += tshift; // qseq and tseq may be shifted due to the removal of leading I/D
-	r->blen = r->mlen = 0;
-	for (k = 0; k < p->n_cigar; ++k) {
-		uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4;
-		if (op == 0) { // match/mismatch
-			int n_ambi = 0, n_diff = 0;
-			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;
-				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;
-			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;
-			s -= q + e * len;
-			if (s < 0) s = 0;
-			qoff += len;
-		} else if (op == 2) { // deletion
-			int n_ambi = 0;
-			for (l = 0; l < len; ++l)
-				if (tseq[toff + l] > 3) ++n_ambi;
-			r->blen += len - n_ambi, p->n_ambi += n_ambi;
-			s -= q + e * len;
-			if (s < 0) s = 0;
-			toff += len;
-		} else if (op == 3) { // intron
-			toff += len;
-		}
-	}
-	p->dp_max = max;
-	assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
-}
-
-static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) // TODO: this calls the libc realloc()
-{
-	mm_extra_t *p;
-	if (n_cigar == 0) return;
-	if (r->p == 0) {
-		uint32_t capacity = n_cigar + sizeof(mm_extra_t)/4;
-		kroundup32(capacity);
-		r->p = (mm_extra_t*)calloc(capacity, 4);
-		r->p->capacity = capacity;
-	} else if (r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4 > r->p->capacity) {
-		r->p->capacity = r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4;
-		kroundup32(r->p->capacity);
-		r->p = (mm_extra_t*)realloc(r->p, r->p->capacity * 4);
-	}
-	p = r->p;
-	if (p->n_cigar > 0 && (p->cigar[p->n_cigar-1]&0xf) == (cigar[0]&0xf)) { // same CIGAR op at the boundary
-		p->cigar[p->n_cigar-1] += cigar[0]>>4<<4;
-		if (n_cigar > 1) memcpy(p->cigar + p->n_cigar, cigar + 1, (n_cigar - 1) * 4);
-		p->n_cigar += n_cigar - 1;
-	} else {
-		memcpy(p->cigar + p->n_cigar, cigar, n_cigar * 4);
-		p->n_cigar += n_cigar;
-	}
-}
-
 static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq) // written by @armintoepfer
 {
 	uint32_t n_EQX = 0;
@@ -290,6 +218,79 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 	r->p = p;
 }
 
+static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e, int is_eqx)
+{
+	uint32_t k, l;
+	int32_t s = 0, max = 0, qshift, tshift, toff = 0, qoff = 0;
+	mm_extra_t *p = r->p;
+	if (p == 0) return;
+	mm_fix_cigar(r, qseq, tseq, &qshift, &tshift);
+	qseq += qshift, tseq += tshift; // qseq and tseq may be shifted due to the removal of leading I/D
+	r->blen = r->mlen = 0;
+	for (k = 0; k < p->n_cigar; ++k) {
+		uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4;
+		if (op == 0) { // match/mismatch
+			int n_ambi = 0, n_diff = 0;
+			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;
+				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;
+			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;
+			s -= q + e * len;
+			if (s < 0) s = 0;
+			qoff += len;
+		} else if (op == 2) { // deletion
+			int n_ambi = 0;
+			for (l = 0; l < len; ++l)
+				if (tseq[toff + l] > 3) ++n_ambi;
+			r->blen += len - n_ambi, p->n_ambi += n_ambi;
+			s -= q + e * len;
+			if (s < 0) s = 0;
+			toff += len;
+		} else if (op == 3) { // intron
+			toff += len;
+		}
+	}
+	p->dp_max = max;
+	assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
+	if (is_eqx) mm_update_cigar_eqx(r, qseq, tseq); // NB: it has to be called here as changes to qseq and tseq are not returned
+}
+
+static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) // TODO: this calls the libc realloc()
+{
+	mm_extra_t *p;
+	if (n_cigar == 0) return;
+	if (r->p == 0) {
+		uint32_t capacity = n_cigar + sizeof(mm_extra_t)/4;
+		kroundup32(capacity);
+		r->p = (mm_extra_t*)calloc(capacity, 4);
+		r->p->capacity = capacity;
+	} else if (r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4 > r->p->capacity) {
+		r->p->capacity = r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4;
+		kroundup32(r->p->capacity);
+		r->p = (mm_extra_t*)realloc(r->p, r->p->capacity * 4);
+	}
+	p = r->p;
+	if (p->n_cigar > 0 && (p->cigar[p->n_cigar-1]&0xf) == (cigar[0]&0xf)) { // same CIGAR op at the boundary
+		p->cigar[p->n_cigar-1] += cigar[0]>>4<<4;
+		if (n_cigar > 1) memcpy(p->cigar + p->n_cigar, cigar + 1, (n_cigar - 1) * 4);
+		p->n_cigar += n_cigar - 1;
+	} else {
+		memcpy(p->cigar + p->n_cigar, cigar, n_cigar * 4);
+		p->n_cigar += n_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)
 {
 	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
@@ -300,7 +301,10 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 		for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr);
 		fputc('\n', stderr);
 	}
-	if (opt->flag & MM_F_SPLICE)
+	if (opt->max_sw_mat > 0 && (int64_t)tlen * qlen > opt->max_sw_mat) {
+		ksw_reset_extz(ez);
+		ez->zdropped = 1;
+	} else if (opt->flag & MM_F_SPLICE)
 		ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, flag, 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);
@@ -414,7 +418,7 @@ static void mm_filter_bad_seeds_alt(void *km, int as1, int cnt1, mm128_t *a, int
 			gap2 = ((int32_t)a[as1 + j].y - (int32_t)a[as1 + j - 1].y) - (int32_t)(a[as1 + j].x - a[as1 + j - 1].x);
 			q_span_pre = a[as1 + j - 1].y >> 32 & 0xff;
 			rs2 = (int32_t)a[as1 + j - 1].x + q_span_pre;
-			qs2 = (int32_t)a[as1 + j - 1].x + q_span_pre;
+			qs2 = (int32_t)a[as1 + j - 1].y + q_span_pre;
 			m = rs2 - re1 < qs2 - qe1? rs2 - re1 : qs2 - qe1;
 			gap2 = gap2 > 0? gap2 : -gap2;
 			if (m > gap1 + gap2) break;
@@ -609,6 +613,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				if (++l > opt->min_cnt) {
 					l = rs0 - x > qs0 - y? rs0 - x : qs0 - y;
 					rs1 = rs0 - l, qs1 = qs0 - l;
+					if (rs1 < 0) rs1 = 0; // not strictly necessary; better have this guard for explicit
 					break;
 				}
 			}
@@ -622,6 +627,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			l = l < rs? l : rs;
 			rs1 = rs1 > rs - l? rs1 : rs - l;
 			rs0 = rs0 < rs1? rs0 : rs1;
+			rs0 = rs0 < rs? rs0 : rs;
 		} else rs0 = rs, qs0 = qs;
 		// compute re0 and qe0
 		re0 = (int32_t)a[r->as + r->cnt - 1].x + 1;
@@ -661,7 +667,7 @@ 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);
 
-	if (qs > 0 && rs > 0) { // left extension
+	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_seq_rev(qs - qs0, qseq);
@@ -748,8 +754,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	assert(re1 - rs1 <= re0 - rs0);
 	if (r->p) {
 		mm_idx_getseq(mi, rid, rs1, re1, tseq);
-		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e);
-		if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r, &qseq0[r->rev][qs1], tseq);
+		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e, opt->flag & MM_F_EQX);
 		if (rev && r->p->trans_strand)
 			r->p->trans_strand ^= 3; // flip to the read strand
 	}
@@ -807,8 +812,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	}
 	r_inv->rs = r1->re + t_off;
 	r_inv->re = r_inv->rs + ez->max_t + 1;
-	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e);
-	if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r_inv, &qseq[q_off], &tseq[t_off]);
+	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, opt->flag & MM_F_EQX);
 	ret = 1;
 end_align1_inv:
 	kfree(km, tseq);

bseq.c

@@ -15,7 +15,7 @@ unsigned char seq_comp_table[256] = {
 	 48,  49,  50,	51,	 52,  53,  54,	55,	 56,  57,  58,	59,	 60,  61,  62,	63,
 	 64, 'T', 'V', 'G', 'H', 'E', 'F', 'C', 'D', 'I', 'J', 'M', 'L', 'K', 'N', 'O',
 	'P', 'Q', 'Y', 'S', 'A', 'A', 'B', 'W', 'X', 'R', 'Z',	91,	 92,  93,  94,	95,
-	 64, 't', 'v', 'g', 'h', 'e', 'f', 'c', 'd', 'i', 'j', 'm', 'l', 'k', 'n', 'o',
+	 96, 't', 'v', 'g', 'h', 'e', 'f', 'c', 'd', 'i', 'j', 'm', 'l', 'k', 'n', 'o',
 	'p', 'q', 'y', 's', 'a', 'a', 'b', 'w', 'x', 'r', 'z', 123, 124, 125, 126, 127,
 	128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
 	144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159,
@@ -39,7 +39,7 @@ mm_bseq_file_t *mm_bseq_open(const char *fn)
 {
 	mm_bseq_file_t *fp;
 	gzFile f;
-	f = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
+	f = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(0, "r");
 	if (f == 0) return 0;
 	fp = (mm_bseq_file_t*)calloc(1, sizeof(mm_bseq_file_t));
 	fp->fp = f;
@@ -65,6 +65,8 @@ static inline char *kstrdup(const kstring_t *s)
 static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual, int with_comment)
 {
 	int i;
+	if (ks->name.l == 0)
+		fprintf(stderr, "[WARNING]\033[1;31m empty sequence name in the input.\033[0m\n");
 	s->name = kstrdup(&ks->name);
 	s->seq = kstrdup(&ks->seq);
 	for (i = 0; i < (int)ks->seq.l; ++i) // convert U to T
@@ -78,6 +80,7 @@ static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual, int with_
 mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int with_comment, int frag_mode, int *n_)
 {
 	int64_t size = 0;
+	int ret;
 	kvec_t(mm_bseq1_t) a = {0,0,0};
 	kseq_t *ks = fp->ks;
 	*n_ = 0;
@@ -87,7 +90,7 @@ mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
 		size = fp->s.l_seq;
 		memset(&fp->s, 0, sizeof(mm_bseq1_t));
 	}
-	while (kseq_read(ks) >= 0) {
+	while ((ret = kseq_read(ks)) >= 0) {
 		mm_bseq1_t *s;
 		assert(ks->seq.l <= INT32_MAX);
 		if (a.m == 0) kv_resize(mm_bseq1_t, 0, a, 256);
@@ -107,6 +110,8 @@ mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
 			break;
 		}
 	}
+	if (ret < -1)
+		fprintf(stderr, "[WARNING]\033[1;31m wrong FASTA/FASTQ record. Continue anyway.\033[0m\n");
 	*n_ = a.n;
 	return a.a;
 }

chain.c

@@ -14,12 +14,12 @@ static const char LogTable256[256] = {
 
 static inline int ilog2_32(uint32_t v)
 {
-	register uint32_t t, tt;
+	uint32_t t, tt;
 	if ((tt = v>>16)) return (t = tt>>8) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
 	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.13/minimap2-2.13_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.13_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
+curl -L https://github.com/lh3/minimap2/releases/download/v2.16/minimap2-2.16_x64-linux.tar.bz2 | tar jxf -
+cp minimap2-2.16_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 -

example.c

@@ -45,7 +45,7 @@ int main(int argc, char *argv[])
 				printf("%s\t%d\t%d\t%d\t%c\t", ks->name.s, ks->seq.l, r->qs, r->qe, "+-"[r->rev]);
 				printf("%s\t%d\t%d\t%d\t%d\t%d\t%d\tcg:Z:", mi->seq[r->rid].name, mi->seq[r->rid].len, r->rs, r->re, r->mlen, r->blen, r->mapq);
 				for (i = 0; i < r->p->n_cigar; ++i) // IMPORTANT: this gives the CIGAR in the aligned regions. NO soft/hard clippings!
-					printf("%d%c", r->p->cigar[i]>>4, "MIDSHN"[r->p->cigar[i]&0xf]);
+					printf("%d%c", r->p->cigar[i]>>4, "MIDNSH"[r->p->cigar[i]&0xf]);
 				putchar('\n');
 				free(r->p);
 			}

format.c

@@ -92,7 +92,8 @@ static void sam_write_rg_line(kstring_t *str, const char *s)
 		if (mm_verbose >= 1) fprintf(stderr, "[ERROR] the read group line contained literal <tab> characters -- replace with escaped tabs: \\t\n");
 		goto err_set_rg;
 	}
-	rg_line = strdup(s);
+	rg_line = (char*)malloc(strlen(s) + 1);
+	strcpy(rg_line, s);
 	mm_escape(rg_line);
 	if ((p = strstr(rg_line, "\tID:")) == 0) {
 		if (mm_verbose >= 1) fprintf(stderr, "[ERROR] no ID within the read group line\n");
@@ -139,8 +140,8 @@ static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 	if (write_tag) mm_sprintf_lite(s, "\tcs:Z:");
 	for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
 		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
-		assert(op >= 0 && op <= 3);
-		if (op == 0) { // match
+		assert((op >= 0 && op <= 3) || op == 7 || op == 8);
+		if (op == 0 || op == 7 || op == 8) { // match
 			int l_tmp = 0;
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
@@ -187,8 +188,8 @@ static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 	if (write_tag) mm_sprintf_lite(s, "\tMD:Z:");
 	for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
 		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
-		assert(op >= 0 && op <= 3);
-		if (op == 0) { // match
+		assert((op >= 0 && op <= 3) || op == 7 || op == 8);
+		if (op == 0 || op == 7 || op == 8) { // match
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
 					mm_sprintf_lite(s, "%d%c", l_MD, "ACGTN"[tseq[t_off + j]]);
@@ -260,6 +261,18 @@ int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_r
 	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 1, 0);
 }
 
+double mm_event_identity(const mm_reg1_t *r)
+{
+	int32_t i, n_gapo = 0, n_gap = 0;
+	if (r->p == 0) return -1.0f;
+	for (i = 0; i < r->p->n_cigar; ++i) {
+		int32_t op = r->p->cigar[i] & 0xf, len = r->p->cigar[i] >> 4;
+		if (op == 1 || op == 2)
+			++n_gapo, n_gap += len;
+	}
+	return (double)r->mlen / (r->blen - n_gap + n_gapo);
+}
+
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 {
 	int type;
@@ -272,20 +285,28 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 	}
 	mm_sprintf_lite(s, "\ttp:A:%c\tcm:i:%d\ts1:i:%d", type, r->cnt, r->score);
 	if (r->parent == r->id) mm_sprintf_lite(s, "\ts2:i:%d", r->subsc);
-	if (r->div >= 0.0f && r->div <= 1.0f) {
-		char buf[8];
+	if (r->p) {
+		char buf[16];
+		double div;
+		div = 1.0 - mm_event_identity(r);
+		if (div == 0.0) buf[0] = '0', buf[1] = 0;
+		else snprintf(buf, 16, "%.4f", 1.0 - mm_event_identity(r));
+		mm_sprintf_lite(s, "\tde:f:%s", buf);
+	} else if (r->div >= 0.0f && r->div <= 1.0f) {
+		char buf[16];
 		if (r->div == 0.0f) buf[0] = '0', buf[1] = 0;
-		else sprintf(buf, "%.4f", r->div);
+		else snprintf(buf, 16, "%.4f", r->div);
 		mm_sprintf_lite(s, "\tdv:f:%s", buf);
 	}
 	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]);
@@ -295,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:");
@@ -307,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];
@@ -348,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]);
@@ -356,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;
@@ -507,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);
@@ -514,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;

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') {
@@ -92,7 +96,7 @@ static int ketopt(ketopt_t *s, int argc, char *argv[], int permute, const char *
 		char *p;
 		if (s->pos == 0) s->pos = 1;
 		opt = s->opt = argv[s->i][s->pos++];
-		p = strchr(ostr, opt);
+		p = strchr((char*)ostr, opt);
 		if (p == 0) {
 			opt = '?'; /* unknown option */
 		} else if (p[1] == ':') {

ksort.h

@@ -37,7 +37,7 @@ typedef struct {
 	int depth;
 } ks_isort_stack_t;
 
-#define KSORT_SWAP(type_t, a, b) { register type_t t=(a); (a)=(b); (b)=t; }
+#define KSORT_SWAP(type_t, a, b) { type_t t=(a); (a)=(b); (b)=t; }
 
 #define KSORT_INIT(name, type_t, __sort_lt) \
 	void ks_heapdown_##name(size_t i, size_t n, type_t l[]) \

ksw2_dispatch.c

@@ -17,18 +17,20 @@
 void __cpuidex(int cpuid[4], int func_id, int subfunc_id)
 {
 #if defined(__x86_64__)
-	asm volatile ("cpuid"
+	__asm__ volatile ("cpuid"
 			: "=a" (cpuid[0]), "=b" (cpuid[1]), "=c" (cpuid[2]), "=d" (cpuid[3])
 			: "0" (func_id), "2" (subfunc_id));
 #else // on 32bit, ebx can NOT be used as PIC code
-	asm volatile ("xchgl %%ebx, %1; cpuid; xchgl %%ebx, %1"
+	__asm__ volatile ("xchgl %%ebx, %1; cpuid; xchgl %%ebx, %1"
 			: "=a" (cpuid[0]), "=r" (cpuid[1]), "=c" (cpuid[2]), "=d" (cpuid[3])
 			: "0" (func_id), "2" (subfunc_id));
 #endif
 }
 #endif
 
-int x86_simd(void)
+static int ksw_simd = -1;
+
+static int x86_simd(void)
 {
 	int flag = 0, cpuid[4], max_id;
 	__cpuidex(cpuid, 0, 0);
@@ -54,11 +56,10 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 {
 	extern void ksw_extz2_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, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
 	extern void ksw_extz2_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, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
-	unsigned simd;
-	simd = x86_simd();
-	if (simd & SIMD_SSE4_1)
+	if (ksw_simd < 0) ksw_simd = x86_simd();
+	if (ksw_simd & SIMD_SSE4_1)
 		ksw_extz2_sse41(km, qlen, query, tlen, target, m, mat, q, e, w, zdrop, end_bonus, flag, ez);
-	else if (simd & SIMD_SSE2)
+	else if (ksw_simd & SIMD_SSE2)
 		ksw_extz2_sse2(km, qlen, query, tlen, target, m, mat, q, e, w, zdrop, end_bonus, flag, ez);
 	else abort();
 }
@@ -70,11 +71,10 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
 	extern void ksw_extd2_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 e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
-	unsigned simd;
-	simd = x86_simd();
-	if (simd & SIMD_SSE4_1)
+	if (ksw_simd < 0) ksw_simd = x86_simd();
+	if (ksw_simd & SIMD_SSE4_1)
 		ksw_extd2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, ez);
-	else if (simd & SIMD_SSE2)
+	else if (ksw_simd & SIMD_SSE2)
 		ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, ez);
 	else abort();
 }
@@ -86,11 +86,10 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, 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);
-	unsigned simd;
-	simd = x86_simd();
-	if (simd & SIMD_SSE4_1)
+	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);
-	else if (simd & SIMD_SSE2)
+	else if (ksw_simd & SIMD_SSE2)
 		ksw_exts2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
 	else abort();
 }

main.c

@@ -6,7 +6,7 @@
 #include "mmpriv.h"
 #include "ketopt.h"
 
-#define MM_VERSION "2.13-r850"
+#define MM_VERSION "2.16-r922"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -60,6 +60,9 @@ static ko_longopt_t long_options[] = {
 	{ "split-prefix",   ko_required_argument, 334 },
 	{ "no-end-flt",     ko_no_argument,       335 },
 	{ "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 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
 	{ "version",        ko_no_argument,       'V' },
@@ -97,7 +100,7 @@ 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;
@@ -122,7 +125,7 @@ int main(int argc, char *argv[])
 			fprintf(stderr, "[ERROR] missing option argument\n");
 			return 1;
 		} else if (c == '?') {
-			fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[o.i]);
+			fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[o.i - 1]);
 			return 1;
 		}
 	}
@@ -163,12 +166,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
@@ -190,6 +202,8 @@ int main(int argc, char *argv[])
 		else if (c == 334) opt.split_prefix = o.arg; // --split-prefix
 		else if (c == 335) opt.flag |= MM_F_NO_END_FLT; // --no-end-flt
 		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 == 314) { // --frag
 			yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1);
 		} else if (c == 315) { // --secondary
@@ -264,7 +278,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");
@@ -289,7 +303,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");

map.c

@@ -284,6 +284,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		qlen_sum += qlens[i], n_regs[i] = 0, regs[i] = 0;
 
 	if (qlen_sum == 0 || n_segs <= 0 || n_segs > MM_MAX_SEG) return;
+	if (opt->max_qlen > 0 && qlen_sum > opt->max_qlen) return;
 
 	hash  = qname? __ac_X31_hash_string(qname) : 0;
 	hash ^= __ac_Wang_hash(qlen_sum) + __ac_Wang_hash(opt->seed);
@@ -311,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;
@@ -333,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;
@@ -583,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
 					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

@@ -107,10 +107,12 @@ typedef struct {
 	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
 
@@ -139,6 +141,7 @@ typedef struct {
 	int32_t mid_occ;     // ignore seeds with occurrences above this threshold
 	int32_t max_occ;
 	int mini_batch_size; // size of a batch of query bases to process in parallel
+	int64_t max_sw_mat;
 
 	const char *split_prefix;
 } mm_mapopt_t;

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "11 October 2018" "minimap2-2.13 (r850)" "Bioinformatics tools"
+.TH minimap2 1 "28 Feburary 2019" "minimap2-2.16-dirty (r922)" "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
@@ -274,6 +280,10 @@ Only map to the reverse complement strand of the reference sequences.
 .BR --heap-sort = no | yes
 If yes, sort anchors with heap merge, instead of radix sort. Heap merge is
 faster for short reads, but slower for long reads. [no]
+.TP
+.B --no-pairing
+Treat two reads in a pair as independent reads. The mate related fields in SAM
+are still properly populated.
 .SS Alignment options
 .TP 10
 .BI -A \ INT
@@ -369,12 +379,22 @@ It helps to avoid tiny terminal exons. [6]
 .B --no-end-flt
 Don't filter seeds towards the ends of chains before performing base-level
 alignment.
+.TP
+.BI --cap-sw-mem \ NUM
+Skip alignment if the DP matrix size is above
+.IR NUM .
+Set 0 to disable [0].
 .SS Input/output options
 .TP 10
 .B -a
 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
@@ -449,6 +469,15 @@ memory.
 .BR --secondary = yes | no
 Whether to output secondary alignments [yes]
 .TP
+.BI --max-qlen \ NUM
+Filter out query sequences longer than
+.IR NUM .
+.TP
+.B --paf-no-hit
+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 --version
 Print version number to stdout
 .SS Preset options
@@ -493,7 +522,7 @@ Up to 10% sequence divergence.
 .B asm20
 Long assembly to reference mapping
 .RB ( -k19
-.B -w10 -A1 -B6 -O6,26 -E2,1 -s200 -z200
+.B -w10 -A1 -B4 -O6,26 -E2,1 -s200 -z200
 .BR --min-occ-floor=100 ).
 Up to 20% sequence divergence.
 .TP
@@ -595,6 +624,8 @@ ts	A	Transcript strand (splice mode only)
 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.c

@@ -116,8 +116,7 @@ long peakrss(void)
 double realtime(void)
 {
 	struct timeval tp;
-	struct timezone tzp;
-	gettimeofday(&tp, &tzp);
+	gettimeofday(&tp, NULL);
 	return tp.tv_sec + tp.tv_usec * 1e-6;
 }
 

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = '2.13-r850';
+var paftools_version = '2.16-r922';
 
 /*****************************
  ***** 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;
@@ -564,14 +565,18 @@ function paf_call(args)
 
 function paf_asmstat(args)
 {
-	var c, min_seg_len = 10000, max_diff = 0.01;
-	while ((c = getopt(args, "l:d:")) != null) {
+	var c, min_query_len = 0, min_seg_len = 10000, max_diff = 0.01, bp_flank_len = 0, bp_gap_len = 0;
+	while ((c = getopt(args, "l:d:b:g:q:")) != null) {
 		if (c == 'l') min_seg_len = parseInt(getopt.arg);
 		else if (c == 'd') max_diff = parseFloat(getopt.arg);
+		else if (c == 'b') bp_flank_len = parseInt(getopt.arg);
+		else if (c == 'g') bp_gap_len = parseInt(getopt.arg);
+		else if (c == 'q') min_query_len = parseInt(getopt.arg);
 	}
 	if (getopt.ind == args.length) {
 		print("Usage: paftools.js asmstat [options] <ref.fa.fai> <asm1.paf> [...]");
 		print("Options:");
+		print("  -q INT     ignore query shorter than INT [0]");
 		print("  -l INT     min alignment block length [" + min_seg_len + "]");
 		print("  -d FLOAT   max gap-compressed sequence divergence [" + max_diff + "]");
 		exit(1);
@@ -587,7 +592,7 @@ function paf_asmstat(args)
 	}
 	file.close();
 
-	function process_query(qblocks, qblock_len, bp) {
+	function process_query(qblocks, qblock_len, bp, qi) {
 		qblocks.sort(function(a,b) { return a[0]-b[0]; });
 		var last_k = null, last_blen = null, st = -1, en = -1, qcov = 0;
 		for (var k = 0; k < qblocks.length; ++k) {
@@ -612,6 +617,7 @@ function paf_asmstat(args)
 				var min = blen < last_blen? blen : last_blen;
 				var flank = k == 0? min : blen;
 				bp.push([flank, gap]);
+				qi.bp.push([flank, gap]);
 			}
 			last_k = k, last_blen = blen;
 		}
@@ -654,7 +660,7 @@ function paf_asmstat(args)
 		return (NM - n_gaps + n_gapo) / (n_M + n_gapo);
 	}
 
-	var labels = ['Length', 'NG50', 'Coverage', 'Qcov', 'NGA50', '#breaks', 'bp(' + min_seg_len + ',0)', 'bp(' + min_seg_len + ',10k)'];
+	var labels = ['Length', 'l_cov', 'Rcov', 'Rdup', 'Qcov', 'NG75', 'NG50', 'NGA50', '#breaks', 'bp(' + min_seg_len + ',0)', 'bp(' + min_seg_len + ',10k)'];
 	var rst = [];
 	for (var i = 0; i < labels.length; ++i)
 		rst[i] = [];
@@ -664,17 +670,23 @@ function paf_asmstat(args)
 	for (var i = 0; i < n_asm; ++i) {
 		var n_breaks = 0, qcov = 0;
 		var fn = args[getopt.ind + 1 + i];
-		header.push(fn.replace(/.paf(.gz)?$/, ""));
+		var label = fn.replace(/.paf(.gz)?$/, "");
+		header.push(label);
 		var ref_blocks = [], qblock_len = [], qblocks = [], bp = [];
-		var query = {};
+		var query = {}, qinfo = {};
 		var last_qname = null;
 		file = new File(fn);
 		while (file.readline(buf) >= 0) {
 			var m, line = buf.toString();
 			var t = line.split("\t");
 			t[1] = parseInt(t[1]);
-			if (t.length >= 2) query[t[0]] = t[1];
-			if (t.length < 9) continue;
+			if (t[1] < min_query_len) 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];
@@ -690,7 +702,7 @@ function paf_asmstat(args)
 			if (t[3] - t[2] < min_seg_len) continue;
 			if (t[0] != last_qname) {
 				if (last_qname != null)
-					qcov += process_query(qblocks, qblock_len, bp);
+					qcov += process_query(qblocks, qblock_len, bp, qinfo[last_qname]);
 				qblocks = [];
 				last_qname = t[0];
 			}
@@ -698,7 +710,7 @@ function paf_asmstat(args)
 			qblocks.push([t[2], t[3], t[4], t[5], t[7], t[8]]);
 		}
 		if (last_qname != null)
-			qcov += process_query(qblocks, qblock_len, bp);
+			qcov += process_query(qblocks, qblock_len, bp, qinfo[last_qname]);
 		file.close();
 
 		// compute NG50
@@ -708,7 +720,8 @@ function paf_asmstat(args)
 			asm_lens.push(query[ctg]);
 		}
 		rst[0][i] = asm_len;
-		rst[1][i] = N50(asm_lens, ref_len, 0.5);
+		rst[5][i] = N50(asm_lens, ref_len, 0.75);
+		rst[6][i] = N50(asm_lens, ref_len, 0.5);
 
 		// compute coverage
 		var l_cov = 0;
@@ -723,22 +736,195 @@ function paf_asmstat(args)
 			} else en = en > ref_blocks[j][2]? en : ref_blocks[j][2];
 		}
 		l_cov += en - st;
+		rst[1][i] = l_cov;
 		rst[2][i] = (100.0 * (l_cov / ref_len)).toFixed(2) + '%';
-		rst[3][i] = (100.0 * (qcov  / asm_len)).toFixed(2) + '%';
+		rst[4][i] = (100.0 * (qcov  / asm_len)).toFixed(2) + '%';
+
+		// compute cov1 and cov2+ lengths; see paf_call() for details
+		var c1_ctg = null, c1_start = 0, c1_end = 0, c1_len = 0;
+		for (var j = 0; j < ref_blocks.length; ++j) {
+			if (ref_blocks[j][0] != c1_ctg || ref_blocks[j][1] >= c1_end) {
+				if (c1_end > c1_start)
+					c1_len += c1_end - c1_start;
+				c1_ctg = ref_blocks[j][0], c1_start = ref_blocks[j][1], c1_end = ref_blocks[j][2];
+			} else if (ref_blocks[j][2] > c1_end) { // overlap
+				if (ref_blocks[j][1] > c1_start)
+					c1_len += ref_blocks[j][1] - c1_start;
+				c1_start = c1_end, c1_end = ref_blocks[j][2];
+			} else if (ref_blocks[j][2] > c1_start) { // contained
+				if (ref_blocks[j][1] > c1_start)
+					c1_len += ref_blocks[j][1] - c1_start;
+				c1_start = ref_blocks[j][2];
+			}
+			//print(ref_blocks[j][0], ref_blocks[j][1], ref_blocks[j][2], c1_start, c1_end, c1_len);
+		}
+		if (c1_end > c1_start)
+			c1_len += c1_end - c1_start;
+		rst[3][i] = (100 * (l_cov - c1_len) / l_cov).toFixed(2) + '%';
 
 		// compute NGA50
-		rst[4][i] = N50(qblock_len, ref_len, 0.5);
+		rst[7][i] = N50(qblock_len, ref_len, 0.5);
 
 		// compute break points
-		rst[5][i] = n_breaks;
-		rst[6][i] = count_bp(bp, 500, 0);
-		rst[7][i] = count_bp(bp, 500, 10000);
+		rst[8][i] = n_breaks;
+		rst[9][i] = count_bp(bp, 500, 0);
+		rst[10][i] = count_bp(bp, 500, 10000);
+
+		// nb-plot; NOT USED
+		/*
+		var qa = [];
+		for (var qn in qinfo)
+			qa.push([qinfo[qn].len, qinfo[qn].bp]);
+		qa = qa.sort(function(a, b) { return b[0] - a[0] });
+		var sum = 0, n_bp = 0, next_quantile = 0.1;
+		for (var j = 0; j < qa.length; ++j) {
+			sum += qa[j][0];
+			for (var k = 0; k < qa[j][1].length; ++k)
+				if (qa[j][1][k][0] >= bp_flank_len && qa[j][1][k][1] >= bp_gap_len)
+					++n_bp;
+			if (sum >= ref_len * next_quantile) {
+				print(label, Math.floor(next_quantile * 100 + .5), qa[j][0], (sum / n_bp).toFixed(0), n_bp);
+				next_quantile += 0.1;
+				if (next_quantile >= 1.0) break;
+			}
+		}
+		*/
+	}
+	buf.destroy();
+
+	if (bp_flank_len <= 0) {
+		print(header.join("\t"));
+		for (var i = 0; i < labels.length; ++i)
+			print(labels[i], rst[i].join("\t"));
+	}
+}
+
+function paf_asmgene(args)
+{
+	var c, opt = { min_cov:0.99, min_iden:0.99 }, print_err = false, auto_only = false;
+	while ((c = getopt(args, "i:c:ea")) != null)
+		if (c == 'i') opt.min_iden = parseFloat(getopt.arg);
+		else if (c == 'c') opt.min_cov = parseFloat(getopt.arg);
+		else if (c == 'e') print_err = true;
+		else if (c == 'a') auto_only = true;
+
+	var n_fn = args.length - getopt.ind;
+	if (n_fn < 2) {
+		print("Usage: paftools.js asmgene [options] <ref-splice.paf> <asm-splice.paf> [...]");
+		print("Options:");
+		print("  -i FLOAT     min identity [" + opt.min_iden + "]");
+		print("  -c FLOAT     min coverage [" + opt.min_cov + "]");
+		print("  -a           only evaluate genes mapped to the autosomes");
+		print("  -e           print fragmented/missing genes");
+		exit(1);
 	}
 
-	print(header.join("\t"));
-	for (var i = 0; i < labels.length; ++i)
-		print(labels[i], rst[i].join("\t"));
+	function process_query(opt, a) {
+		var b = [], cnt = [0, 0, 0];
+		for (var j = 0; j < a.length; ++j) {
+			if (a[j][4] < a[j][5] * opt.min_iden)
+				continue;
+			b.push(a[j].slice(0));
+		}
+		if (b.length == 0) return cnt;
+		// count full
+		var n_full = 0;
+		for (var j = 0; j < b.length; ++j)
+			if (b[j][3] - b[j][2] >= b[j][1] * opt.min_cov)
+				++n_full;
+		cnt[0] = n_full;
+		// compute coverage
+		b = b.sort(function(x, y) { return x[2] - y[2] });
+		var l_cov = 0, st = b[0][2], en = b[0][3];
+		for (var j = 1; j < b.length; ++j) {
+			if (b[j][2] <= en)
+				en = b[j][3] > en? b[j][3] : en;
+			else l_cov += en - st;
+		}
+		l_cov += en - st;
+		cnt[1] = l_cov / b[0][1];
+		cnt[2] = b.length;
+		return cnt;
+	}
 
+	var buf = new Bytes();
+	var gene = {}, header = [], refpos = {};
+	for (var i = getopt.ind; i < args.length; ++i) {
+		var fn = args[i];
+		var label = fn.replace(/.paf(.gz)?$/, "");
+		header.push(label);
+		var file = new File(fn), a = [];
+		while (file.readline(buf) >= 0) {
+			var t = buf.toString().split("\t");
+			var ql = parseInt(t[1]), qs = parseInt(t[2]), qe = parseInt(t[3]), mlen = parseInt(t[9]), blen = parseInt(t[10]), mapq = parseInt(t[11]);
+			if (i == getopt.ind) refpos[t[0]] = [t[0], t[1], t[5], t[7], t[8]];
+			if (gene[t[0]] == null) gene[t[0]] = [];
+			if (a.length && t[0] != a[0][0]) {
+				gene[a[0][0]][i - getopt.ind] = process_query(opt, a);
+				a = [];
+			}
+			a.push([t[0], ql, qs, qe, mlen, blen]);
+		}
+		if (a.length)
+			gene[t[0]][i - getopt.ind] = process_query(opt, a);
+		file.close();
+	}
+
+	// select the longest genes (not optimal, but should be good enough)
+	var gene_list = [], gene_nr = {};
+	for (var g in refpos)
+		gene_list.push(refpos[g]);
+	gene_list = gene_list.sort(function(a, b) { return a[2] < b[2]? -1 : a[2] > b[2]? 1 : a[3] - b[3] });
+	var last = 0;
+	for (var j = 1; j < gene_list.length; ++j) {
+		if (gene_list[j][2] != gene_list[last][2] || gene_list[j][3] >= gene_list[last][4]) {
+			gene_nr[gene_list[last][0]] = 1;
+			last = j;
+		} else if (gene_list[j][1] > gene_list[last][1]) {
+			last = j;
+		}
+	}
+	gene_nr[gene_list[last][0]] = 1;
+
+	// count and print
+	var col1 = ["full_sgl", "full_dup", "frag", "part50+", "part10+", "part10-"];
+	var rst = [];
+	for (var k = 0; k < col1.length; ++k) {
+		rst[k] = [];
+		for (var i = 0; i < n_fn; ++i)
+			rst[k][i] = 0;
+	}
+	for (var g in gene) {
+		if (gene[g][0] == null || gene[g][0][0] != 1) continue;
+		if (gene_nr[g] == null) continue;
+		if (auto_only && /^(chr)?[XY]$/.test(refpos[g][2])) continue;
+		for (var i = 0; i < n_fn; ++i) {
+			if (gene[g][i] == null) {
+				rst[4][i]++;
+				if (print_err) print('M', header[i], refpos[g].join("\t"));
+			} else if (gene[g][i][0] == 1) rst[0][i]++;
+			else if (gene[g][i][0] > 1) {
+				rst[1][i]++;
+				if (print_err) print('D', header[i], refpos[g].join("\t"));
+			} else if (gene[g][i][1] >= opt.min_cov) {
+				rst[2][i]++;
+				if (print_err) print('F', header[i], refpos[g].join("\t"));
+			} else if (gene[g][i][1] >= 0.5) {
+				rst[3][i]++;
+				if (print_err) print('5', header[i], refpos[g].join("\t"));
+			} else if (gene[g][i][1] >= 0.1) {
+				rst[4][i]++;
+				if (print_err) print('1', header[i], refpos[g].join("\t"));
+			} else {
+				rst[5][i]++;
+				if (print_err) print('0', header[i], refpos[g].join("\t")); // TODO: reduce code duplicates...
+			}
+		}
+	}
+	print('H', 'Metric', header.join("\t"));
+	for (var k = 0; k < rst.length; ++k) {
+		print('X', col1[k], rst[k].join("\t"));
+	}
 	buf.destroy();
 }
 
@@ -781,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;
@@ -1014,6 +1202,105 @@ function paf_bedcov(args)
 	warn("# target bases overlapping regions: " + hit_len + ' (' + (100.0 * hit_len / tot_len).toFixed(2) + '%)');
 }
 
+function paf_vcfpair(args)
+{
+	var c, is_male = false, sample = 'syndip', hgver = null;
+	var PAR = { '37':[[0, 2699520], [154931043, 155260560]] };
+	while ((c = getopt(args, "ms:g:")) != null) {
+		if (c == 'm') is_male = true;
+		else if (c == 's') sample = getopt.arg;
+		else if (c == 'g') hgver = getopt.arg;
+	}
+	if (is_male && (hgver == null || PAR[hgver] == null))
+		throw("for a male, -g must be specified to properly handle PARs on chrX");
+
+	if (getopt.ind == args.length) {
+		print("Usage: paftools.js vcfpair [options] <in.pair.vcf>");
+		print("Options:");
+		print("  -m       the sample is male");
+		print("  -g STR   human genome version '37' []");
+		print("  -s STR   sample name [" + sample + "]");
+		exit(1);
+	}
+
+	var re_ctg = is_male? /^(chr)?([0-9]+|X|Y)$/ : /^(chr)?([0-9]+|X)$/;
+	var label = ['1', '2'];
+	var buf = new Bytes();
+	var file = args[getopt.ind] == '-'? new File() : new File(args[getopt.ind]);
+	while (file.readline(buf) >= 0) {
+		var m, line = buf.toString();
+		if (line.charAt(0) == '#') {
+			if (/^##(source|reference)=/.test(line)) continue;
+			if ((m = /^##contig=.*ID=([^\s,]+)/.exec(line)) != null) {
+				if (!re_ctg.test(m[1])) continue;
+			} else if (/^#CHROM/.test(line)) {
+				var t = line.split("\t");
+				--t.length;
+				t[t.length-1] = sample;
+				line = t.join("\t");
+				print('##FILTER=<ID=HET1,Description="Heterozygous in the first haplotype">');
+				print('##FILTER=<ID=HET2,Description="Heterozygous in the second haplotype">');
+				print('##FILTER=<ID=GAP1,Description="Uncalled in the first haplotype">');
+				print('##FILTER=<ID=GAP2,Description="Uncalled in the second haplotype">');
+			}
+			print(line);
+			continue;
+		}
+		var t = line.split("\t");
+		if (!re_ctg.test(t[0])) continue;
+		var GT = null, AD = null, FILTER = [], HT = [null, null];
+		for (var i = 0; i < 2; ++i) {
+			if ((m = /^(\.|[0-9]+)\/(\.|[0-9]+):(\S+)/.exec(t[9+i])) == null) {
+				warn(line);
+				throw Error("malformatted VCF");
+			}
+			var s = m[3].split(",");
+			if (AD == null) {
+				AD = [];
+				for (var j = 0; j < s.length; ++j)
+					AD[j] = 0;
+			}
+			for (var j = 0; j < s.length; ++j)
+				AD[j] += parseInt(s[j]);
+			if (m[1] == '.') {
+				FILTER.push('GAP' + label[i]);
+				HT[i] = '.';
+			} else if (m[1] != m[2]) {
+				FILTER.push('HET' + label[i]);
+				HT[i] = '.';
+			} else HT[i] = m[1];
+		}
+		--t.length;
+		// test if this is in a haploid region
+		var hap = 0, st = parseInt(t[1]), en = st + t[3].length;
+		if (is_male) {
+			if (/^(chr)?X/.test(t[0])) {
+				if (hgver != null && PAR[hgver] != null) {
+					var r = PAR[hgver], in_par = false;
+					for (var i = 0; i < r.length; ++i)
+						if (r[i][0] <= st && en <= r[i][1])
+							in_par = true;
+					hap = in_par? 0 : 2;
+				}
+			} else if (/^(chr)?Y/.test(t[0])) {
+				hap = 1;
+			}
+		}
+		// special treatment for haploid regions
+		if (hap > 0 && FILTER.length == 1) {
+			if ((hap == 2 && FILTER[0] == "GAP1") || (hap == 1 && FILTER[0] == "GAP2"))
+				FILTER.length = 0;
+		}
+		// update VCF
+		t[5] = 30; // fake QUAL
+		t[6] = FILTER.length? FILTER.join(";") : ".";
+		t[9] = HT.join("|") + ":" + AD.join(",");
+		print(t.join("\t"));
+	}
+	file.close();
+	buf.destroy();
+}
+
 /**********************
  * Conversion related *
  **********************/
@@ -1305,11 +1592,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);
 	}
 
@@ -1338,13 +1630,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) {
@@ -1380,8 +1673,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;
@@ -1395,13 +1688,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;
 							}
@@ -1446,7 +1741,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"));
@@ -2189,6 +2485,7 @@ function main(args)
 		print("");
 		print("  stat       collect basic mapping information in PAF/SAM");
 		print("  asmstat    collect basic assembly information");
+		print("  asmgene    evaluate gene completeness (EXPERIMENTAL)");
 		print("  liftover   simplistic liftOver");
 		print("  call       call variants from asm-to-ref alignment with the cs tag");
 		print("  bedcov     compute the number of bases covered");
@@ -2210,7 +2507,9 @@ function main(args)
 	else if (cmd == 'gff2bed') paf_gff2bed(args);
 	else if (cmd == 'stat') paf_stat(args);
 	else if (cmd == 'asmstat') paf_asmstat(args);
+	else if (cmd == 'asmgene') paf_asmgene(args);
 	else if (cmd == 'liftover' || cmd == 'liftOver') paf_liftover(args);
+	else if (cmd == 'vcfpair') paf_vcfpair(args);
 	else if (cmd == 'call') paf_call(args);
 	else if (cmd == 'mapeval') paf_mapeval(args);
 	else if (cmd == 'bedcov') paf_bedcov(args);

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;
@@ -159,6 +160,11 @@ int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
 			fprintf(stderr, "[ERROR]\033[1;31m --for-only and --rev-only can't be applied at the same time\033[0m\n");
 		return -3;
 	}
+	if (mo->e <= 0 || mo->q <= 0) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m -O and -E must be positive\033[0m\n");
+		return -1;
+	}
 	if ((mo->q != mo->q2 || mo->e != mo->e2) && !(mo->e > mo->e2 && mo->q + mo->e < mo->q2 + mo->e2)) {
 		if (mm_verbose >= 1)
 			fprintf(stderr, "[ERROR]\033[1;31m dual gap penalties violating E1>E2 and O1+E1<O2+E2\033[0m\n");

pe.c

@@ -54,7 +54,7 @@ void mm_set_pe_thru(const int *qlens, int *n_regs, mm_reg1_t **regs)
 	if (n_pri[0] == 1 && n_pri[1] == 1) {
 		mm_reg1_t *p = &regs[0][pri[0]];
 		mm_reg1_t *q = &regs[1][pri[1]];
-		if (p->rid == q->rid && p->rev == q->rev && abs(p->rs - q->rs) < 3 && abs(p->re - p->re) < 3
+		if (p->rid == q->rid && p->rev == q->rev && abs(p->rs - q->rs) < 3 && abs(p->re - q->re) < 3
 			&& ((p->qs == 0 && qlens[1] - q->qe == 0) || (q->qs == 0 && qlens[0] - p->qe == 0)))
 		{
 			p->pe_thru = q->pe_thru = 1;

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

@@ -13,10 +13,11 @@ cdef extern from "minimap.h":
 		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,7 +25,7 @@ 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
@@ -40,6 +41,7 @@ cdef extern from "minimap.h":
 		int32_t mid_occ
 		int32_t max_occ
 		int mini_batch_size
+		int64_t max_sw_mat
 		const char *split_prefix
 
 	int mm_set_opt(char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)

python/mappy.pyx

@@ -3,7 +3,7 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
-__version__ = '2.13'
+__version__ = '2.16'
 
 cmappy.mm_reset_timer()
 
@@ -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.13',
+	version = '2.16',
 	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) {