Release minimap2-2.15 (r905)

This was caused by mm_fix_cigar() which may change query/target offset in very
rare cases. Generating EQX has to beware of this change.

Resolves #266

NEWS.md

@@ -1,3 +1,65 @@
+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

@@ -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.15/minimap2-2.15_x64-linux.tar.bz2 | tar -jxvf -
+./minimap2-2.15_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
@@ -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;
@@ -591,11 +595,9 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 		qs0 = 0, qe0 = qlen;
 		l = qs;
 		l += l * opt->a + opt->end_bonus > opt->q? (l * opt->a + opt->end_bonus - opt->q) / opt->e : 0;
-		l = l < opt->bw? l : opt->bw;
 		rs0 = rs - l > 0? rs - l : 0;
 		l = qlen - qe;
 		l += l * opt->a + opt->end_bonus > opt->q? (l * opt->a + opt->end_bonus - opt->q) / opt->e : 0;
-		l = l < opt->bw? l : opt->bw;
 		re0 = re + l < (int32_t)mi->seq[rid].len? re + l : mi->seq[rid].len;
 	} else {
 		// compute rs0 and qs0
@@ -611,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;
 				}
 			}
@@ -624,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;
@@ -663,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);
@@ -750,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
 	}
@@ -809,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,7 +14,7 @@ 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];
 }

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.15/minimap2-2.15_x64-linux.tar.bz2 | tar jxf -
+cp minimap2-2.15_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 - r->p->n_ambi - n_gap + n_gapo);
+}
+
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 {
 	int type;
@@ -272,10 +285,17 @@ 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);

ketopt.h

@@ -92,7 +92,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-r852-dirty"
+#define MM_VERSION "2.15-r905"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -60,6 +60,8 @@ 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 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
 	{ "version",        ko_no_argument,       'V' },
@@ -122,7 +124,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;
 		}
 	}
@@ -190,6 +192,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

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

minimap.h

@@ -107,6 +107,8 @@ 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;
@@ -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 "10 January 2019" "minimap2-2.15 (r905)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -274,6 +274,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,6 +373,11 @@ 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
@@ -449,6 +458,13 @@ 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 query name and length for an unmapped sequence.
+.TP
 .B --version
 Print version number to stdout
 .SS Preset options
@@ -493,7 +509,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 +611,7 @@ 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
 .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.15-r905';
 
 /*****************************
  ***** Library functions *****
@@ -564,14 +564,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 +591,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 +616,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 +659,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,16 +669,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[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 (!/\ttp:A:[PI]/.test(line)) continue;
 			if ((m = /\tcg:Z:(\S+)/.exec(line)) == null) continue;
@@ -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();
 }
 
@@ -1014,6 +1200,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 *
  **********************/
@@ -2189,6 +2474,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 +2496,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);

options.c

@@ -159,6 +159,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/cmappy.pxd

@@ -40,6 +40,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.15'
 
 cmappy.mm_reset_timer()
 

setup.py

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