bedidx.c

   1 #include <stdlib.h>
   2 #include <stdint.h>
   3 #include <string.h>
   4 #include <stdio.h>
   5 #include <zlib.h>
   6
   7 #include "kseq.h"
   8 KSTREAM_INIT(gzFile, gzread, 8192)
   9
  10 typedef struct {
  11         int n, m;
  12         uint64_t *a;
  13         int *idx;
  14 } bed_reglist_t;
  15
  16 #include "khash.h"
  17 KHASH_MAP_INIT_STR(reg, bed_reglist_t)
  18
  19 #define LIDX_SHIFT 13
  20
  21 typedef kh_reg_t reghash_t;
  22
  23 int *bed_index_core(int n, uint64_t *a, int *n_idx)
  24 {
  25         int i, j, m, *idx;
  26         m = *n_idx = 0; idx = 0;
  27         for (i = 0; i < n; ++i) {
  28                 int beg, end;
  29                 beg = a[i]>>32 >> LIDX_SHIFT; end = ((uint32_t)a[i]) >> LIDX_SHIFT;
  30                 if (m < end + 1) {
  31                         int oldm = m;
  32                         m = end + 1;
  33                         kroundup32(m);
  34                         idx = realloc(idx, m * sizeof(int));
  35                         for (j = oldm; j < m; ++j) idx[j] = -1;
  36                 }
  37                 if (beg == end) {
  38                         if (idx[beg] < 0) idx[beg] = i;
  39                 } else {
  40                         for (j = beg; j <= end; ++j)
  41                                 if (idx[j] < 0) idx[j] = i;
  42                 }
  43                 *n_idx = end + 1;
  44         }
  45         return idx;
  46 }
  47
  48 void bed_index(void *_h)
  49 {
  50         reghash_t *h = (reghash_t*)_h;
  51         khint_t k;
  52         for (k = 0; k < kh_end(h); ++k) {
  53                 if (kh_exist(h, k)) {
  54                         bed_reglist_t *p = &kh_val(h, k);
  55                         if (p->idx) free(p->idx);
  56                         p->idx = bed_index_core(p->n, p->a, &p->m);
  57                 }
  58         }
  59 }
  60
  61 int bed_overlap_core(const bed_reglist_t *p, int beg, int end)
  62 {
  63         int i, min_off;
  64         if (p->n == 0) return 0;
  65         min_off = (beg>>LIDX_SHIFT >= p->n)? p->idx[p->n-1] : p->idx[beg>>LIDX_SHIFT];
  66         if (min_off < 0) { // TODO: this block can be improved, but speed should not matter too much here
  67                 int n = beg>>LIDX_SHIFT;
  68                 if (n > p->n) n = p->n;
  69                 for (i = n - 1; i >= 0; --i)
  70                         if (p->idx[i] >= 0) break;
  71                 min_off = i >= 0? p->idx[i] : 0;
  72         }
  73         for (i = min_off; i < p->n; ++i) {
  74                 if ((int)(p->a[i]>>32) >= end) break; // out of range; no need to proceed
  75                 if ((int32_t)p->a[i] > beg && (int32_t)(p->a[i]>>32) < end)
  76                         return 1; // find the overlap; return
  77         }
  78         return 0;
  79 }
  80
  81 int bed_overlap(const void *_h, const char *chr, int beg, int end)
  82 {
  83         const reghash_t *h = (const reghash_t*)_h;
  84         khint_t k;
  85         if (!h) return 0;
  86         k = kh_get(reg, h, chr);
  87         if (k == kh_end(h)) return 0;
  88         return bed_overlap_core(&kh_val(h, k), beg, end);
  89 }
  90
  91 void *bed_read(const char *fn)
  92 {
  93         reghash_t *h = kh_init(reg);
  94         gzFile fp;
  95         kstream_t *ks;
  96         int dret;
  97         kstring_t *str;
  98         // read the list
  99         fp = strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
 100         if (fp == 0) return 0;
 101         str = calloc(1, sizeof(kstring_t));
 102         ks = ks_init(fp);
 103         while (ks_getuntil(ks, 0, str, &dret) >= 0) { // read the chr name
 104                 int beg = -1, end = -1;
 105                 bed_reglist_t *p;
 106                 khint_t k = kh_get(reg, h, str->s);
 107                 if (k == kh_end(h)) { // absent from the hash table
 108                         int ret;
 109                         char *s = strdup(str->s);
 110                         k = kh_put(reg, h, s, &ret);
 111                         memset(&kh_val(h, k), 0, sizeof(bed_reglist_t));
 112                 }
 113                 p = &kh_val(h, k);
 114                 if (dret != '\n') { // if the lines has other characters
 115                         if (ks_getuntil(ks, 0, str, &dret) > 0 && isdigit(str->s[0])) {
 116                                 beg = atoi(str->s); // begin
 117                                 if (dret != '\n') {
 118                                         if (ks_getuntil(ks, 0, str, &dret) > 0 && isdigit(str->s[0]))
 119                                                 end = atoi(str->s); // end
 120                                 }
 121                         }
 122                 }
 123                 if (dret != '\n') while ((dret = ks_getc(ks)) > 0 && dret != '\n'); // skip the rest of the line
 124                 if (end < 0 && beg > 0) end = beg, beg = beg - 1; // if there is only one column
 125                 if (beg >= 0 && end > beg) {
 126                         if (p->n == p->m) {
 127                                 p->m = p->m? p->m<<1 : 4;
 128                                 p->a = realloc(p->a, p->m * 8);
 129                         }
 130                         p->a[p->n++] = (uint64_t)beg<<32 | end;
 131                 }
 132         }
 133         ks_destroy(ks);
 134         gzclose(fp);
 135         free(str->s); free(str);
 136         bed_index(h);
 137         return h;
 138 }
 139
 140 void bed_destroy(void *_h)
 141 {
 142         reghash_t *h = (reghash_t*)_h;
 143         khint_t k;
 144         for (k = 0; k < kh_end(h); ++k) {
 145                 if (kh_exist(h, k)) {
 146                         free(kh_val(h, k).a);
 147                         free(kh_val(h, k).idx);
 148                         free((char*)kh_key(h, k));
 149                 }
 150         }
 151         kh_destroy(reg, h);
 152 }