debugging code

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;
@@ -305,7 +308,7 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 	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);
@@ -547,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
 
@@ -662,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;
@@ -693,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);
@@ -703,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);
@@ -732,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;
@@ -749,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)
@@ -790,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;
@@ -809,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:

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

@@ -64,18 +64,18 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 }
 
 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);
+				   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, ez);
+		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, ez);
+		ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, qd, 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

@@ -60,6 +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 },
+	{ "bed",            ko_required_argument, 337 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
 	{ "version",        ko_no_argument,       'V' },
@@ -102,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;
@@ -190,6 +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) 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
@@ -348,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;
 
@@ -362,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);

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

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

@@ -53,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
 
@@ -62,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
 
@@ -112,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

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