diff --git a/index.c b/index.c index 05420d8..5be102c 100644 --- a/index.c +++ b/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) @@ -585,3 +585,85 @@ 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 +#include +#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_read_bed_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; + for (p = q = str.s, i = 0;; ++p) { + if (*p == 0 || isspace(*p)) { + int32_t c = *p; + *p = 0; + if (i == 0) { + id = mm_idx_name2id(mi, q); + if (id < 0) break; // unknown name; TODO: throw a warning + } else if (i == 1) { + st = atoi(q); + if (st < 0) break; + } else if (i == 2) { + en = atoi(q); + if (en < 0) break; + } else break; + if (c == 0) break; + ++i, q = p + 1; + } + } + if (i == 2) en = st + 1; + if (i < 2 || st >= en) continue; + if (m == n) EXPAND(r, m); + t.x = (uint64_t)id << 32 | st, t.end = en, t.idx = -1; + r[n++] = t; + } + ks_destroy(ks); + gzclose(fp); + *n_ = n; + return r; +} + +int mm_idx_attach_bed(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) +{ + 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); +} + +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; +} diff --git a/kseq.h b/kseq.h index d301ddc..8021e56 100644 --- a/kseq.h +++ b/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) { \ diff --git a/main.c b/main.c index 0e0812b..07795a9 100644 --- a/main.c +++ b/main.c @@ -59,6 +59,7 @@ static ko_longopt_t long_options[] = { { "paf-no-hit", ko_no_argument, 333 }, { "split-prefix", ko_required_argument, 334 }, { "no-end-flt", ko_no_argument, 335 }, + { "bed", ko_required_argument, 336 }, { "help", ko_no_argument, 'h' }, { "max-intron-len", ko_required_argument, 'G' }, { "version", ko_no_argument, 'V' }, @@ -101,7 +102,7 @@ int main(int argc, char *argv[]) mm_mapopt_t opt; mm_idxopt_t ipt; int i, c, n_threads = 3, n_parts; - 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; @@ -188,6 +189,7 @@ int main(int argc, char *argv[]) else if (c == 333) opt.flag |= MM_F_PAF_NO_HIT; // --paf-no-hit 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) 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 @@ -342,6 +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 (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 27ba911..6bae3fb 100644 --- a/minimap.h +++ b/minimap.h @@ -58,13 +58,20 @@ typedef struct { uint32_t len; // length } mm_idx_seq_t; +typedef struct { + uint64_t x; + int32_t end, idx; +} 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; @@ -365,6 +372,9 @@ int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, ui 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/mmpriv.h b/mmpriv.h index f847445..9e4a95b 100644 --- a/mmpriv.h +++ b/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