debugging code

NEWS.md

@@ -1,29 +1,3 @@
-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.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.13/minimap2-2.13_x64-linux.tar.bz2 | tar -jxvf -
+./minimap2-2.13_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,8 +355,6 @@ 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

@@ -88,13 +88,14 @@ static int mm_test_zdrop(void *km, const mm_mapopt_t *opt, const uint8_t *qseq,
 	return max_zdrop > opt->zdrop? 1 : 0;
 }
 
-static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int *qshift, int *tshift)
+static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int *qshift, int *tshift, int left_aln)
 {
 	mm_extra_t *p = r->p;
 	int32_t toff = 0, qoff = 0, to_shrink = 0;
 	uint32_t k;
 	*qshift = *tshift = 0;
 	if (p->n_cigar <= 1) return;
+	if (!left_aln) goto end_left_aln;
 	for (k = 0; k < p->n_cigar; ++k) { // indel left alignment
 		uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4;
 		if (len == 0) to_shrink = 1;
@@ -123,6 +124,7 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 		}
 	}
 	assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
+end_left_aln:
 	if (to_shrink) { // squeeze out zero-length operations
 		int32_t l = 0;
 		for (k = 0; k < p->n_cigar; ++k) // squeeze out zero-length operations
@@ -147,13 +149,13 @@ 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)
+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 left_aln)
 {
 	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);
+	mm_fix_cigar(r, qseq, tseq, &qshift, &tshift, left_aln);
 	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) {
@@ -290,7 +292,8 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 	r->p = p;
 }
 
-static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
+static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const int8_t *ob, const int8_t *mat,
+						  int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
 {
 	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
 		int i;
@@ -300,15 +303,12 @@ 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->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)
+	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);
 	else
-		ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez);
+		ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ob, ez);
 	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
 		int i;
 		fprintf(stderr, "score=%d, cigar=", ez->score);
@@ -417,7 +417,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].y + q_span_pre;
+			qs2 = (int32_t)a[as1 + j - 1].x + q_span_pre;
 			m = rs2 - re1 < qs2 - qe1? rs2 - re1 : qs2 - qe1;
 			gap2 = gap2 > 0? gap2 : -gap2;
 			if (m > gap1 + gap2) break;
@@ -550,7 +550,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 	int32_t i, l, bw, dropped = 0, extra_flag = 0, rs0, re0, qs0, qe0;
 	int32_t rs, re, qs, qe;
 	int32_t rs1, qs1, re1, qe1;
-	int8_t mat[25];
+	int8_t mat[25], *ob;
 
 	if (is_sr) assert(!(mi->flag & MM_I_HPC)); // HPC won't work with SR because with HPC we can't easily tell if there is a gap
 
@@ -594,9 +594,11 @@ 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
@@ -663,13 +665,14 @@ 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);
+	ob = (int8_t*)kmalloc(km, re0 - rs0);
 
 	if (qs > 0 && rs > 0) { // left extension
 		qseq = &qseq0[rev][qs0];
-		mm_idx_getseq(mi, rid, rs0, rs, tseq);
+		mm_idx_getseq2(mi, rid, rs0, rs, tseq, ob);
 		mm_seq_rev(qs - qs0, qseq);
 		mm_seq_rev(rs - rs0, tseq);
-		mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, mat, bw, opt->end_bonus, r->split_inv? opt->zdrop_inv : opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez);
+		mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, ob, mat, bw, opt->end_bonus, r->split_inv? opt->zdrop_inv : opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez);
 		if (ez->n_cigar > 0) {
 			mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			r->p->dp_score += ez->max;
@@ -694,7 +697,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				bw1 = qe - qs > re - rs? qe - qs : re - rs;
 			// perform alignment
 			qseq = &qseq0[rev][qs];
-			mm_idx_getseq(mi, rid, rs, re, tseq);
+			mm_idx_getseq2(mi, rid, rs, re, tseq, ob);
 			if (is_sr) { // perform ungapped alignment
 				assert(qe - qs == re - rs);
 				ksw_reset_extz(ez);
@@ -704,11 +707,11 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				}
 				ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, 0, qe - qs);
 			} else { // perform normal gapped alignment
-				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
+				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, ob, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
 			}
 			// test Z-drop and inversion Z-drop
 			if ((zdrop_code = mm_test_zdrop(km, opt, qseq, tseq, ez->n_cigar, ez->cigar, mat)) != 0)
-				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, zdrop_code == 2? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate
+				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, ob, mat, bw1, -1, zdrop_code == 2? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate
 			// update CIGAR
 			if (ez->n_cigar > 0)
 				mm_append_cigar(r, ez->n_cigar, ez->cigar);
@@ -733,8 +736,8 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 
 	if (!dropped && qe < qe0 && re < re0) { // right extension
 		qseq = &qseq0[rev][qe];
-		mm_idx_getseq(mi, rid, re, re0, tseq);
-		mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
+		mm_idx_getseq2(mi, rid, re, re0, tseq, ob);
+		mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, ob, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
 		if (ez->n_cigar > 0) {
 			mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			r->p->dp_score += ez->max;
@@ -750,14 +753,16 @@ 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);
+		int left_aln = (opt->flag & MM_F_SPLICE) || mi->n_R == 0? 1 : 0;
+		mm_idx_getseq2(mi, rid, rs1, re1, tseq, ob);
+		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e, left_aln);
 		if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r, &qseq0[r->rev][qs1], tseq);
 		if (rev && r->p->trans_strand)
 			r->p->trans_strand ^= 3; // flip to the read strand
 	}
 
 	kfree(km, tseq);
+	kfree(km, ob);
 }
 
 static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez)
@@ -791,7 +796,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	mm_seq_rev(tl, tseq);
 	if (score < opt->min_dp_max) goto end_align1_inv;
 	q_off = ql - (q_off + 1), t_off = tl - (t_off + 1);
-	mm_align_pair(km, opt, ql - q_off, qseq + q_off, tl - t_off, tseq + t_off, mat, (int)(opt->bw * 1.5), -1, opt->zdrop, KSW_EZ_EXTZ_ONLY, ez);
+	mm_align_pair(km, opt, ql - q_off, qseq + q_off, tl - t_off, tseq + t_off, 0, mat, (int)(opt->bw * 1.5), -1, opt->zdrop, KSW_EZ_EXTZ_ONLY, ez);
 	if (ez->n_cigar == 0) goto end_align1_inv; // should never be here
 	mm_append_cigar(r_inv, ez->n_cigar, ez->cigar);
 	r_inv->p->dp_score = ez->max;
@@ -810,7 +815,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);
+	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, 0);
 	if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r_inv, &qseq[q_off], &tseq[t_off]);
 	ret = 1;
 end_align1_inv:

bseq.c

@@ -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(0, "r");
+	f = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
 	if (f == 0) return 0;
 	fp = (mm_bseq_file_t*)calloc(1, sizeof(mm_bseq_file_t));
 	fp->fp = f;
@@ -65,8 +65,6 @@ 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
@@ -80,7 +78,6 @@ 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;
@@ -90,7 +87,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 ((ret = kseq_read(ks)) >= 0) {
+	while (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);
@@ -110,8 +107,6 @@ 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;
 }

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.13/minimap2-2.13_x64-linux.tar.bz2 | tar jxf -
+cp minimap2-2.13_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
 export PATH="$PATH:"`pwd`                               # put the current directory on PATH
 # download example datasets
 curl -L https://github.com/lh3/minimap2/releases/download/v2.10/cookbook-data.tgz | tar zxf -

format.c

@@ -92,8 +92,7 @@ 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 = (char*)malloc(strlen(s) + 1);
-	strcpy(rg_line, s);
+	rg_line = strdup(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");
@@ -140,8 +139,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) || op == 7 || op == 8);
-		if (op == 0 || op == 7 || op == 8) { // match
+		assert(op >= 0 && op <= 3);
+		if (op == 0) { // match
 			int l_tmp = 0;
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
@@ -188,8 +187,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) || op == 7 || op == 8);
-		if (op == 0 || op == 7 || op == 8) { // match
+		assert(op >= 0 && op <= 3);
+		if (op == 0) { // 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]]);

index.c

@@ -60,7 +60,7 @@ void mm_idx_destroy(mm_idx_t *mi)
 			free(mi->seq[i].name);
 		free(mi->seq);
 	} else km_destroy(mi->km);
-	free(mi->B); free(mi->S); free(mi);
+	free(mi->R); free(mi->B); free(mi->S); free(mi);
 }
 
 const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n)
@@ -134,7 +134,7 @@ int mm_idx_name2id(const mm_idx_t *mi, const char *name)
 	return k == kh_end(h)? -1 : kh_val(h, k);
 }
 
-int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq)
+int mm_idx_getseq2(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq, int8_t *b)
 {
 	uint64_t i, st1, en1;
 	if (rid >= mi->n_seq || st >= mi->seq[rid].len) return -1;
@@ -143,9 +143,26 @@ int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, ui
 	en1 = mi->seq[rid].offset + en;
 	for (i = st1; i < en1; ++i)
 		seq[i - st1] = mm_seq4_get(mi->S, i);
+	if (b) memset(b, 0, en - st);
+	if (b && mi->R) {
+		uint32_t i, z;
+		memset(b, 0, en - st);
+		z = mm_idx_bed_query(mi, (uint64_t)rid << 32 | st);
+		for (i = z < 0? 0 : z; i < mi->n_R; ++i) {
+			uint32_t j, rr, rs, re;
+			rr = mi->R[i].x >> 32, rs = (uint32_t)mi->R[i].x, re = mi->R[i].end;
+			if (rr > rid || rs >= en) break;
+			if (rr < rid) continue;
+			re = re < en? re : en;
+			for (j = st > rs? st : rs; j < re; ++j)
+				b[j - st] = mi->R[i].score;
+		}
+	}
 	return en - st;
 }
 
+int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq) { return mm_idx_getseq2(mi, rid, st, en, seq, 0); }
+
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f)
 {
 	int i;
@@ -585,3 +602,89 @@ int mm_idx_reader_eof(const mm_idx_reader_t *r) // TODO: in extremely rare cases
 {
 	return r->is_idx? (feof(r->fp.idx) || ftell(r->fp.idx) == r->idx_size) : mm_bseq_eof(r->fp.seq);
 }
+
+#include <ctype.h>
+#include <zlib.h>
+#include "ksort.h"
+#include "kseq.h"
+KSTREAM_DECLARE(gzFile, gzread)
+
+#define sort_key_bed(a) ((a).x)
+KRADIX_SORT_INIT(bed, mm_idx_bed_t, sort_key_bed, 8)
+
+mm_idx_bed_t *mm_idx_bed_read_list(const mm_idx_t *mi, const char *fn, uint32_t *n_)
+{
+	gzFile fp;
+	kstream_t *ks;
+	kstring_t str = {0,0,0};
+	uint32_t n = 0, m = 0;
+	mm_idx_bed_t *r = 0;
+
+	fp = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
+	if (fp == 0) return 0;
+	ks = ks_init(fp);
+	while (ks_getuntil(ks, KS_SEP_LINE, &str, 0) >= 0) {
+		mm_idx_bed_t t;
+		char *p, *q;
+		int i, id = -1, st = -1, en = -1, sc = -1;
+		for (p = q = str.s, i = 0;; ++p) {
+			if (*p == 0 || isspace(*p)) {
+				int32_t c = *p;
+				*p = 0;
+				if (i == 0) { // chr
+					id = mm_idx_name2id(mi, q);
+					if (id < 0) break; // unknown name; TODO: throw a warning
+				} else if (i == 1) { // start
+					st = atoi(q);
+					if (st < 0) break;
+				} else if (i == 2) { // end
+					en = atoi(q);
+					if (en < 0) break;
+				} else if (i == 3) { // name; do nothing
+				} else if (i == 4) { // BED score
+					sc = atoi(q);
+					assert(sc >= 0 && sc <= 127);
+				} else break;
+				if (c == 0) break;
+				++i, q = p + 1;
+			}
+		}
+		if (en < 0) en = st + 1;
+		if (st < 0 || st >= en) continue;
+		if (m == n) EXPAND(r, m);
+		t.x = (uint64_t)id << 32 | st, t.end = en, t.score = sc >= 0? sc : 0, t.idx = -1;
+		r[n++] = t;
+	}
+	ks_destroy(ks);
+	gzclose(fp);
+	*n_ = n;
+	return r;
+}
+
+int mm_idx_bed_attach(mm_idx_t *mi, uint32_t n, mm_idx_bed_t *r) // TODO: check errors
+{
+	radix_sort_bed(r, r + n);
+	mi->R = r, mi->n_R = n;
+	return 0;
+}
+
+int mm_idx_bed_read(mm_idx_t *mi, const char *fn)
+{
+	mm_idx_bed_t *r;
+	uint32_t n;
+	if (mi->h == 0) mm_idx_index_name(mi);
+	r = mm_idx_bed_read_list(mi, fn, &n);
+	return mm_idx_bed_attach(mi, n, r);
+}
+
+int mm_idx_bed_query(const mm_idx_t *mi, uint64_t x)
+{
+	int32_t left = -1, right = mi->n_R;
+	while (right - left > 1) {
+		int32_t mid = left + ((right - left) >> 1);
+		if (mi->R[mid].x > x) right = mid;
+		else if (mi->R[mid].x < x) left = mid;
+		else return mid;
+	}
+	return left;
+}

kseq.h

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

ksw2.h

@@ -58,7 +58,7 @@ void ksw_extd(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t
 			  int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int flag, ksw_extz_t *ez);
 
 void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
+				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez);
 
 void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
 				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);

ksw2_dispatch.c

@@ -17,20 +17,18 @@
 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
 
-static int ksw_simd = -1;
-
-static int x86_simd(void)
+int x86_simd(void)
 {
 	int flag = 0, cpuid[4], max_id;
 	__cpuidex(cpuid, 0, 0);
@@ -56,26 +54,28 @@ 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);
-	if (ksw_simd < 0) ksw_simd = x86_simd();
-	if (ksw_simd & SIMD_SSE4_1)
+	unsigned simd;
+	simd = x86_simd();
+	if (simd & SIMD_SSE4_1)
 		ksw_extz2_sse41(km, qlen, query, tlen, target, m, mat, q, e, w, zdrop, end_bonus, flag, ez);
-	else if (ksw_simd & SIMD_SSE2)
+	else if (simd & SIMD_SSE2)
 		ksw_extz2_sse2(km, qlen, query, tlen, target, m, mat, q, e, w, zdrop, end_bonus, flag, ez);
 	else abort();
 }
 
 void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez)
 {
 	extern void ksw_extd2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, 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);
-	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 (ksw_simd & SIMD_SSE2)
-		ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, ez);
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez);
+	unsigned simd;
+	simd = x86_simd();
+	if (simd & SIMD_SSE4_1)
+		ksw_extd2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, qd, ez);
+	else if (simd & SIMD_SSE2)
+		ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, qd, ez);
 	else abort();
 }
 
@@ -86,10 +86,11 @@ 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);
-	if (ksw_simd < 0) ksw_simd = x86_simd();
-	if (ksw_simd & SIMD_SSE4_1)
+	unsigned simd;
+	simd = x86_simd();
+	if (simd & SIMD_SSE4_1)
 		ksw_exts2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
-	else if (ksw_simd & SIMD_SSE2)
+	else if (simd & SIMD_SSE2)
 		ksw_exts2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
 	else abort();
 }

ksw2_extd2_sse.c

@@ -17,17 +17,18 @@
 #ifdef KSW_CPU_DISPATCH
 #ifdef __SSE4_1__
 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)
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez)
 #else
 void ksw_extd2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez)
 #endif
 #else
 void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez)
 #endif // ~KSW_CPU_DISPATCH
 {
 #define __dp_code_block1 \
+	dt = _mm_load_si128(&dv[t]); \
 	z = _mm_load_si128(&s[t]); \
 	xt1 = _mm_load_si128(&x[t]);                     /* xt1 <- x[r-1][t..t+15] */ \
 	tmp = _mm_srli_si128(xt1, 15);                   /* tmp <- x[r-1][t+15] */ \
@@ -50,10 +51,10 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 #define __dp_code_block2 \
 	_mm_store_si128(&u[t], _mm_sub_epi8(z, vt1));    /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \
 	_mm_store_si128(&v[t], _mm_sub_epi8(z, ut));     /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \
-	tmp = _mm_sub_epi8(z, q_); \
+	tmp = _mm_sub_epi8(z, _mm_add_epi8(dt, q_)); \
 	a = _mm_sub_epi8(a, tmp); \
 	b = _mm_sub_epi8(b, tmp); \
-	tmp = _mm_sub_epi8(z, q2_); \
+	tmp = _mm_sub_epi8(z, _mm_add_epi8(dt, q2_)); \
 	a2= _mm_sub_epi8(a2, tmp); \
 	b2= _mm_sub_epi8(b2, tmp);
 
@@ -61,8 +62,8 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX);
 	int32_t *H = 0, H0 = 0, last_H0_t = 0;
 	uint8_t *qr, *sf, *mem, *mem2 = 0;
-	__m128i q_, q2_, qe_, qe2_, zero_, sc_mch_, sc_mis_, m1_, sc_N_;
-	__m128i *u, *v, *x, *y, *x2, *y2, *s, *p = 0;
+	__m128i q_, q2_, qe_, qe2_, e_, e2_, zero_, sc_mch_, sc_mis_, m1_, sc_N_;
+	__m128i *u, *v, *x, *y, *x2, *y2, *s, *p = 0, *dv;
 
 	ksw_reset_extz(ez);
 	if (m <= 1 || qlen <= 0 || tlen <= 0) return;
@@ -72,6 +73,8 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	zero_   = _mm_set1_epi8(0);
 	q_      = _mm_set1_epi8(q);
 	q2_     = _mm_set1_epi8(q2);
+	e_      = _mm_set1_epi8(e);
+	e2_     = _mm_set1_epi8(e2);
 	qe_     = _mm_set1_epi8(q + e);
 	qe2_    = _mm_set1_epi8(q2 + e2);
 	sc_mch_ = _mm_set1_epi8(mat[0]);
@@ -96,16 +99,23 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		++long_thres;
 	long_diff = long_thres * (e - e2) - (q2 - q) - e2;
 
-	mem = (uint8_t*)kcalloc(km, tlen_ * 8 + qlen_ + 1, 16);
+	mem = (uint8_t*)kcalloc(km, tlen_ * 9 + qlen_ + 1, 16);
 	u = (__m128i*)(((size_t)mem + 15) >> 4 << 4); // 16-byte aligned
-	v = u + tlen_, x = v + tlen_, y = x + tlen_, x2 = y + tlen_, y2 = x2 + tlen_;
-	s = y2 + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16;
+	v = u + tlen_, x = v + tlen_, y = x + tlen_, x2 = y + tlen_, y2 = x2 + tlen_, dv = y2 + tlen_;
+	s = dv + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16;
 	memset(u,  -q  - e,  tlen_ * 16);
 	memset(v,  -q  - e,  tlen_ * 16);
 	memset(x,  -q  - e,  tlen_ * 16);
 	memset(y,  -q  - e,  tlen_ * 16);
 	memset(x2, -q2 - e2, tlen_ * 16);
 	memset(y2, -q2 - e2, tlen_ * 16);
+	if (qd) {
+		int8_t *tmp = (int8_t*)dv;
+		for (t = 0; t < tlen; ++t) tmp[t] = -qd[t];
+		fprintf(stderr, "%d\t%d\t%x\n", tlen, qlen, flag&KSW_EZ_RIGHT);
+		for (t = 0; t < tlen; ++t) fputc("ACGTN"[target[t]], stderr); fputc('\n', stderr);
+		for (t = 0; t < tlen; ++t) fputc('0' + qd[t], stderr); fputc('\n', stderr);
+	}
 	if (!approx_max) {
 		H = (int32_t*)kmalloc(km, tlen_ * 16 * 4);
 		for (t = 0; t < tlen_ * 16; ++t) H[t] = KSW_NEG_INF;
@@ -182,7 +192,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		assert(en_ - st_ + 1 <= n_col_);
 		if (!with_cigar) { // score only
 			for (t = st_; t <= en_; ++t) {
-				__m128i z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp;
+				__m128i z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp, dt;
 				__dp_code_block1;
 #ifdef __SSE4_1__
 				z = _mm_max_epi8(z, a);
@@ -191,10 +201,10 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				z = _mm_max_epi8(z, b2);
 				z = _mm_min_epi8(z, sc_mch_);
 				__dp_code_block2; // save u[] and v[]; update a, b, a2 and b2
-				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_max_epi8(a,  zero_), qe_));
-				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_max_epi8(b,  zero_), qe_));
-				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_max_epi8(a2, zero_), qe2_));
-				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_max_epi8(b2, zero_), qe2_));
+				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_max_epi8(a,  zero_), _mm_add_epi8(dt, qe_)));
+				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_max_epi8(b,  zero_), _mm_add_epi8(dt, qe_)));
+				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_max_epi8(a2, zero_), _mm_add_epi8(dt, qe2_)));
+				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_max_epi8(b2, zero_), _mm_add_epi8(dt, qe2_)));
 #else
 				tmp = _mm_cmpgt_epi8(a,  z);
 				z = _mm_or_si128(_mm_andnot_si128(tmp, z), _mm_and_si128(tmp, a));
@@ -208,20 +218,20 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				z = _mm_or_si128(_mm_and_si128(tmp, sc_mch_), _mm_andnot_si128(tmp, z));
 				__dp_code_block2;
 				tmp = _mm_cmpgt_epi8(a, zero_);
-				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_and_si128(tmp, a),  qe_));
+				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_and_si128(tmp, a),  _mm_add_epi8(dt, qe_)));
 				tmp = _mm_cmpgt_epi8(b, zero_);
-				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_and_si128(tmp, b),  qe_));
+				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_and_si128(tmp, b),  _mm_add_epi8(dt, qe_)));
 				tmp = _mm_cmpgt_epi8(a2, zero_);
-				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), qe2_));
+				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), _mm_add_epi8(dt, qe2_)));
 				tmp = _mm_cmpgt_epi8(b2, zero_);
-				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), qe2_));
+				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), _mm_add_epi8(dt, qe2_)));
 #endif
 			}
 		} else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment
 			__m128i *pr = p + (size_t)r * n_col_ - st_;
 			off[r] = st, off_end[r] = en;
 			for (t = st_; t <= en_; ++t) {
-				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp;
+				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp, dt;
 				__dp_code_block1;
 #ifdef __SSE4_1__
 				d = _mm_and_si128(_mm_cmpgt_epi8(a, z), _mm_set1_epi8(1));       // d = a  > z? 1 : 0
@@ -251,16 +261,16 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 #endif
 				__dp_code_block2;
 				tmp = _mm_cmpgt_epi8(a, zero_);
-				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_and_si128(tmp, a),  qe_));
+				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_and_si128(tmp, a),  _mm_add_epi8(dt, qe_)));
 				d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x08))); // d = a > 0? 1<<3 : 0
 				tmp = _mm_cmpgt_epi8(b, zero_);
-				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_and_si128(tmp, b),  qe_));
+				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_and_si128(tmp, b),  _mm_add_epi8(dt, qe_)));
 				d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x10))); // d = b > 0? 1<<4 : 0
 				tmp = _mm_cmpgt_epi8(a2, zero_);
-				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), qe2_));
+				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), _mm_add_epi8(dt, qe2_)));
 				d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x20))); // d = a > 0? 1<<5 : 0
 				tmp = _mm_cmpgt_epi8(b2, zero_);
-				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), qe2_));
+				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), _mm_add_epi8(dt, qe2_)));
 				d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x40))); // d = b > 0? 1<<6 : 0
 				_mm_store_si128(&pr[t], d);
 			}
@@ -268,7 +278,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 			__m128i *pr = p + (size_t)r * n_col_ - st_;
 			off[r] = st, off_end[r] = en;
 			for (t = st_; t <= en_; ++t) {
-				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp;
+				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp, dt;
 				__dp_code_block1;
 #ifdef __SSE4_1__
 				d = _mm_andnot_si128(_mm_cmpgt_epi8(z, a), _mm_set1_epi8(1));    // d = z > a?  0 : 1
@@ -298,16 +308,16 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 #endif
 				__dp_code_block2;
 				tmp = _mm_cmpgt_epi8(zero_, a);
-				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_andnot_si128(tmp, a),  qe_));
+				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_andnot_si128(tmp, a),  _mm_add_epi8(dt, qe_)));
 				d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x08))); // d = a > 0? 1<<3 : 0
 				tmp = _mm_cmpgt_epi8(zero_, b);
-				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_andnot_si128(tmp, b),  qe_));
+				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_andnot_si128(tmp, b),  _mm_add_epi8(dt, qe_)));
 				d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x10))); // d = b > 0? 1<<4 : 0
 				tmp = _mm_cmpgt_epi8(zero_, a2);
-				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, a2), qe2_));
+				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, a2), _mm_add_epi8(dt, qe2_)));
 				d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x20))); // d = a > 0? 1<<5 : 0
 				tmp = _mm_cmpgt_epi8(zero_, b2);
-				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, b2), qe2_));
+				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, b2), _mm_add_epi8(dt, qe2_)));
 				d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x40))); // d = b > 0? 1<<6 : 0
 				_mm_store_si128(&pr[t], d);
 			}
@@ -376,6 +386,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		last_st = st, last_en = en;
 		//for (t = st0; t <= en0; ++t) printf("(%d,%d)\t(%d,%d,%d,%d)\t%d\n", r, t, ((int8_t*)u)[t], ((int8_t*)v)[t], ((int8_t*)x)[t], ((int8_t*)y)[t], H[t]); // for debugging
 	}
+	fprintf(stderr, "score: %d\n", ez->score);
 	kfree(km, mem);
 	if (!approx_max) kfree(km, H);
 	if (with_cigar) { // backtrack

main.c

@@ -6,7 +6,7 @@
 #include "mmpriv.h"
 #include "ketopt.h"
 
-#define MM_VERSION "2.14-r883"
+#define MM_VERSION "2.13-r852-dirty"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -60,7 +60,7 @@ 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 },
+	{ "bed",            ko_required_argument, 337 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
 	{ "version",        ko_no_argument,       'V' },
@@ -103,7 +103,7 @@ int main(int argc, char *argv[])
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
 	int i, c, n_threads = 3, n_parts, old_best_n = -1;
-	char *fnw = 0, *rg = 0, *s;
+	char *fnw = 0, *fn_bed = 0, *rg = 0, *s;
 	FILE *fp_help = stderr;
 	mm_idx_reader_t *idx_rdr;
 	mm_idx_t *mi;
@@ -123,7 +123,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 - 1]);
+			fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[o.i]);
 			return 1;
 		}
 	}
@@ -191,7 +191,7 @@ 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 == 337) fn_bed = o.arg; // --bed-prefer
 		else if (c == 314) { // --frag
 			yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1);
 		} else if (c == 315) { // --secondary
@@ -350,6 +350,7 @@ int main(int argc, char *argv[])
 			fprintf(stderr, "[M::%s::%.3f*%.2f] loaded/built the index for %d target sequence(s)\n",
 					__func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), mi->n_seq);
 		if (argc != o.ind + 1) mm_mapopt_update(&opt, mi);
+		if (fn_bed) mm_idx_bed_read(mi, fn_bed);
 		if (mm_verbose >= 3) mm_idx_stat(mi);
 		if (!(opt.flag & MM_F_FRAG_MODE)) {
 			for (i = o.ind + 1; i < argc; ++i)

minimap.h

@@ -59,13 +59,21 @@ typedef struct {
 	uint32_t len;    // length
 } mm_idx_seq_t;
 
+typedef struct {
+	uint64_t x;
+	int32_t end, idx;
+	int32_t score; // NB: wasting 4 bytes due to memory alignment
+} mm_idx_bed_t;
+
 typedef struct {
 	int32_t b, w, k, flag;
 	uint32_t n_seq;            // number of reference sequences
 	int32_t index;
+	uint32_t n_R;
 	mm_idx_seq_t *seq;         // sequence name, length and offset
 	uint32_t *S;               // 4-bit packed sequence
 	struct mm_idx_bucket_s *B; // index (hidden)
+	mm_idx_bed_t *R;
 	void *km, *h;
 } mm_idx_t;
 
@@ -139,7 +147,6 @@ 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;
@@ -363,6 +370,11 @@ int mm_idx_index_name(mm_idx_t *mi);
 int mm_idx_name2id(const mm_idx_t *mi, const char *name);
 int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
 
+// BED operations
+int mm_idx_bed_read(mm_idx_t *mi, const char *fn);
+int mm_idx_bed_attach(mm_idx_t *mi, uint32_t n, mm_idx_bed_t *r);
+int mm_idx_bed_query(const mm_idx_t *mi, uint64_t x);
+
 // deprecated APIs for backward compatibility
 void mm_mapopt_init(mm_mapopt_t *opt);
 mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads);

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "5 November 2018" "minimap2-2.14 (r883)" "Bioinformatics tools"
+.TH minimap2 1 "11 October 2018" "minimap2-2.13 (r850)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -274,10 +274,6 @@ 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
@@ -373,11 +369,6 @@ 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
@@ -502,7 +493,7 @@ Up to 10% sequence divergence.
 .B asm20
 Long assembly to reference mapping
 .RB ( -k19
-.B -w10 -A1 -B4 -O6,26 -E2,1 -s200 -z200
+.B -w10 -A1 -B6 -O6,26 -E2,1 -s200 -z200
 .BR --min-occ-floor=100 ).
 Up to 20% sequence divergence.
 .TP

misc.c

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

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = '2.14-r883';
+var paftools_version = '2.13-r850';
 
 /*****************************
  ***** Library functions *****
@@ -564,12 +564,10 @@ 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_seg_len = 10000, max_diff = 0.01;
+	while ((c = getopt(args, "l:d:")) != 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);
 	}
 	if (getopt.ind == args.length) {
 		print("Usage: paftools.js asmstat [options] <ref.fa.fai> <asm1.paf> [...]");
@@ -589,7 +587,7 @@ function paf_asmstat(args)
 	}
 	file.close();
 
-	function process_query(qblocks, qblock_len, bp, qi) {
+	function process_query(qblocks, qblock_len, bp) {
 		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) {
@@ -614,7 +612,6 @@ 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;
 		}
@@ -667,22 +664,16 @@ 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];
-		var label = fn.replace(/.paf(.gz)?$/, "");
-		header.push(label);
+		header.push(fn.replace(/.paf(.gz)?$/, ""));
 		var ref_blocks = [], qblock_len = [], qblocks = [], bp = [];
-		var query = {}, qinfo = {};
+		var query = {};
 		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 (qinfo[t[0]] == null) qinfo[t[0]] = {};
-				qinfo[t[0]].len = t[1];
-				qinfo[t[0]].bp = [];
-			}
+			if (t.length >= 2) query[t[0]] = t[1];
 			if (t.length < 9) continue;
 			if (!/\ttp:A:[PI]/.test(line)) continue;
 			if ((m = /\tcg:Z:(\S+)/.exec(line)) == null) continue;
@@ -699,7 +690,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, qinfo[last_qname]);
+					qcov += process_query(qblocks, qblock_len, bp);
 				qblocks = [];
 				last_qname = t[0];
 			}
@@ -707,7 +698,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, qinfo[last_qname]);
+			qcov += process_query(qblocks, qblock_len, bp);
 		file.close();
 
 		// compute NG50
@@ -742,157 +733,12 @@ function paf_asmstat(args)
 		rst[5][i] = n_breaks;
 		rst[6][i] = count_bp(bp, 500, 0);
 		rst[7][i] = count_bp(bp, 500, 10000);
-
-		// nb-plot
-		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;
-	while ((c = getopt(args, "i:c:e")) != 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;
-
-	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("  -e           print fragmented/missing genes");
-		exit(1);
 	}
 
-	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;
-	}
+	print(header.join("\t"));
+	for (var i = 0; i < labels.length; ++i)
+		print(labels[i], rst[i].join("\t"));
 
-	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[t[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;
-		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();
 }
 
@@ -2343,7 +2189,6 @@ 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");
@@ -2365,7 +2210,6 @@ 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 == 'call') paf_call(args);
 	else if (cmd == 'mapeval') paf_mapeval(args);

mmpriv.h

@@ -31,6 +31,13 @@
 #define MALLOC(type, len) ((type*)malloc((len) * sizeof(type)))
 #define CALLOC(type, len) ((type*)calloc((len), sizeof(type)))
 
+#define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr))))
+
+#define EXPAND(a, m) do { \
+		(m) = (m)? (m) + ((m)>>1) : 16; \
+		REALLOC((a), (m)); \
+	} while (0)
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -67,6 +74,7 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 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);
+int mm_idx_getseq2(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq, int8_t *b);
 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);
 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);
 

options.c

@@ -159,11 +159,6 @@ 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 - q->re) < 3
+		if (p->rid == q->rid && p->rev == q->rev && abs(p->rs - q->rs) < 3 && abs(p->re - p->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.h

@@ -20,6 +20,11 @@ typedef struct {
 	uint32_t *cigar32;
 } mm_hitpy_t;
 
+typedef struct {
+	int32_t n, m;
+	mm_idx_bed_t *r;
+} mm_bedpy_t;
+
 static inline void mm_reg2hitpy(const mm_idx_t *mi, mm_reg1_t *r, mm_hitpy_t *h)
 {
 	h->ctg = mi->seq[r->rid].name;
@@ -149,4 +154,33 @@ static mm_idx_t *mappy_idx_seq(int w, int k, int is_hpc, int bucket_bits, const
 	return mi;
 }
 
+static mm_bedpy_t *mappy_bed_new(void)
+{
+	return (mm_bedpy_t*)calloc(1, sizeof(mm_bedpy_t));
+}
+
+static int mappy_bed_add(mm_bedpy_t *bed, mm_idx_t *mi, const char *name, uint32_t st, uint32_t en)
+{
+	mm_idx_bed_t *b;
+	int id;
+	if (mi->h == 0) mm_idx_index_name(mi);
+	if (bed->n == bed->m) {
+		bed->m = bed->m? bed->m + (bed->m>>1) : 16;
+		bed->r = (mm_idx_bed_t*)realloc(bed->r, sizeof(mm_idx_bed_t) * bed->m);
+	}
+	id = mm_idx_name2id(mi, name);
+	if (id < 0 || st >= en) return -1;
+	if (en > mi->seq[id].len) en = mi->seq[id].len;
+	b = &bed->r[bed->n++];
+	b->x = (uint64_t)id << 32 | st;
+	b->end = en, b->idx = -1;
+	return 0;
+}
+
+static void mappy_bed_finalize(mm_bedpy_t *bed, mm_idx_t *mi)
+{
+	mm_idx_bed_attach(mi, bed->n, bed->r); // bed->r is now owned by mi and will be deallocated with it
+	free(bed);
+}
+
 #endif

python/cmappy.pxd

@@ -40,7 +40,6 @@ 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)
@@ -54,6 +53,10 @@ cdef extern from "minimap.h":
 		uint64_t offset
 		uint32_t len
 
+	ctypedef struct mm_idx_bed_t:
+		uint64_t x
+		int32_t end, idx
+
 	ctypedef struct mm_idx_bucket_t:
 		pass
 
@@ -63,6 +66,7 @@ cdef extern from "minimap.h":
 		mm_idx_seq_t *seq
 		uint32_t *S
 		mm_idx_bucket_t *B
+		mm_idx_bed_t *R
 		void *km
 		void *h
 
@@ -113,6 +117,14 @@ cdef extern from "cmappy.h":
 	char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int en, int *l)
 	mm_idx_t *mappy_idx_seq(int w, int k, int is_hpc, int bucket_bits, const char *seq, int l)
 
+	ctypedef struct mm_bedpy_t:
+		int32_t n, m
+		mm_idx_bed_t *r
+
+	mm_bedpy_t *mappy_bed_new()
+	int mappy_bed_add(mm_bedpy_t *bed, mm_idx_t *mi, const char *name, uint32_t st, uint32_t en)
+	void mappy_bed_finalize(mm_bedpy_t *bed, mm_idx_t *mi)
+
 	ctypedef struct kstring_t:
 		unsigned l, m
 		char *s

python/mappy.pyx

@@ -3,7 +3,7 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
-__version__ = '2.14'
+__version__ = '2.13'
 
 cmappy.mm_reset_timer()
 
@@ -112,7 +112,7 @@ cdef class Aligner:
 	cdef cmappy.mm_idxopt_t idx_opt
 	cdef cmappy.mm_mapopt_t map_opt
 
-	def __cinit__(self, fn_idx_in=None, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, best_n=None, n_threads=3, fn_idx_out=None, max_frag_len=None, extra_flags=None, seq=None, scoring=None):
+	def __cinit__(self, fn_idx_in=None, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, best_n=None, n_threads=3, fn_idx_out=None, max_frag_len=None, extra_flags=None, seq=None, scoring=None, bed=None):
 		cmappy.mm_set_opt(NULL, &self.idx_opt, &self.map_opt) # set the default options
 		if preset is not None:
 			cmappy.mm_set_opt(str.encode(preset), &self.idx_opt, &self.map_opt) # apply preset
@@ -137,6 +137,7 @@ cdef class Aligner:
 					self.map_opt.sc_ambi = scoring[6]
 
 		cdef cmappy.mm_idx_reader_t *r;
+		cdef cmappy.mm_bedpy_t *bed_agg;
 
 		if seq is None:
 			if fn_idx_out is None:
@@ -153,6 +154,14 @@ cdef class Aligner:
 			cmappy.mm_mapopt_update(&self.map_opt, self._idx)
 			self.map_opt.mid_occ = 1000 # don't filter high-occ seeds
 
+		if bed is not None:
+			bed_agg = cmappy.mappy_bed_new()
+			for b in bed:
+				if len(b) < 3: en = int(b[1]) + 1
+				else: en = int(b[2])
+				cmappy.mappy_bed_add(bed_agg, self._idx, str.encode(b[0]), int(b[1]), en)
+			cmappy.mappy_bed_finalize(bed_agg, self._idx)
+
 	def __dealloc__(self):
 		if self._idx is not NULL:
 			cmappy.mm_idx_destroy(self._idx)

setup.py

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