diff --git a/index.c b/index.c index 6a3160a..9fe8e3c 100644 --- a/index.c +++ b/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) diff --git a/main.c b/main.c index 07795a9..4eb24b4 100644 --- a/main.c +++ b/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) diff --git a/minimap.h b/minimap.h index 6bae3fb..50e084f 100644 --- a/minimap.h +++ b/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 diff --git a/python/cmappy.h b/python/cmappy.h index 6bc5635..a9215be 100644 --- a/python/cmappy.h +++ b/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 diff --git a/python/cmappy.pxd b/python/cmappy.pxd index 28ce3f7..620fc9c 100644 --- a/python/cmappy.pxd +++ b/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 diff --git a/python/mappy.pyx b/python/mappy.pyx index 3516acb..9e7f842 100644 --- a/python/mappy.pyx +++ b/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)