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,39 @@
+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)
 -----------------------------------
 

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.14/minimap2-2.14_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.14_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

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) {
@@ -612,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;
 				}
 			}
@@ -625,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;
@@ -664,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);
@@ -751,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
 	}
@@ -810,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,

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.14/minimap2-2.14_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.14_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

@@ -261,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;
@@ -273,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[]) \

main.c

@@ -6,7 +6,7 @@
 #include "mmpriv.h"
 #include "ketopt.h"
 
-#define MM_VERSION "2.14-r883"
+#define MM_VERSION "2.15-r905"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -61,6 +61,7 @@ static ko_longopt_t long_options[] = {
 	{ "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' },
@@ -192,6 +193,7 @@ int main(int argc, char *argv[])
 		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;

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "5 November 2018" "minimap2-2.14 (r883)" "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
@@ -458,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
@@ -604,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/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = '2.14-r883';
+var paftools_version = '2.15-r905';
 
 /*****************************
  ***** Library functions *****
@@ -564,16 +564,18 @@ function paf_call(args)
 
 function paf_asmstat(args)
 {
-	var c, min_seg_len = 10000, max_diff = 0.01, bp_flank_len = 0, bp_gap_len = 0;
-	while ((c = getopt(args, "l:d:b:g:")) != 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);
@@ -657,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] = [];
@@ -677,6 +679,7 @@ function paf_asmstat(args)
 			var m, line = buf.toString();
 			var t = line.split("\t");
 			t[1] = parseInt(t[1]);
+			if (t[1] < min_query_len) continue;
 			if (t.length >= 2) {
 				query[t[0]] = t[1];
 				if (qinfo[t[0]] == null) qinfo[t[0]] = {};
@@ -717,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;
@@ -732,18 +736,42 @@ 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
+		// nb-plot; NOT USED
+		/*
 		var qa = [];
 		for (var qn in qinfo)
 			qa.push([qinfo[qn].len, qinfo[qn].bp]);
@@ -760,6 +788,7 @@ function paf_asmstat(args)
 				if (next_quantile >= 1.0) break;
 			}
 		}
+		*/
 	}
 	buf.destroy();
 
@@ -772,11 +801,12 @@ function paf_asmstat(args)
 
 function paf_asmgene(args)
 {
-	var c, opt = { min_cov:0.99, min_iden:0.99 }, print_err = false;
-	while ((c = getopt(args, "i:c:e")) != null)
+	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) {
@@ -784,6 +814,7 @@ function paf_asmgene(args)
 		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);
 	}
@@ -829,7 +860,7 @@ function paf_asmgene(args)
 			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[t[0]][i - getopt.ind] = process_query(opt, a);
+				gene[a[0][0]][i - getopt.ind] = process_query(opt, a);
 				a = [];
 			}
 			a.push([t[0], ql, qs, qe, mlen, blen]);
@@ -866,6 +897,7 @@ function paf_asmgene(args)
 	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]++;
@@ -1168,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 *
  **********************/
@@ -2367,6 +2498,7 @@ function main(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);

python/mappy.pyx

@@ -3,7 +3,7 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
-__version__ = '2.14'
+__version__ = '2.15'
 
 cmappy.mm_reset_timer()
 

setup.py

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