mappy to support BED

index.c

@@ -611,7 +611,7 @@ 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_read_bed_list(const mm_idx_t *mi, const char *fn, uint32_t *n_)
+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;
@@ -656,20 +656,20 @@ mm_idx_bed_t *mm_idx_read_bed_list(const mm_idx_t *mi, const char *fn, uint32_t
 	return r;
 }
 
-int mm_idx_attach_bed(mm_idx_t *mi, uint32_t n, mm_idx_bed_t *r) // TODO: check errors
+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_read_bed(mm_idx_t *mi, const char *fn)
+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_read_bed_list(mi, fn, &n);
-	return mm_idx_attach_bed(mi, n, r);
+	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)

main.c

@@ -344,7 +344,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_read_bed(mi, fn_bed);
+		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

@@ -368,13 +368,15 @@ 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);
 
-int mm_idx_read_bed(mm_idx_t *mi, const char *fn);
-int mm_idx_bed_query(const mm_idx_t *mi, uint64_t x);
-
 #ifdef __cplusplus
 }
 #endif

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

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