index.c

   1 #include <ctype.h>
   2 #include <assert.h>
   3 #include "khash.h"
   4 #include "ksort.h"
   5 #include "kstring.h"
   6 #include "bam_endian.h"
   7 #ifdef _USE_KNETFILE
   8 #include "knetfile.h"
   9 #endif
  10 #include "tabix.h"
  11
  12 #define TAD_MIN_CHUNK_GAP 32768
  13 // 1<<14 is the size of minimum bin.
  14 #define TAD_LIDX_SHIFT    14
  15
  16 typedef struct {
  17         uint64_t u, v;
  18 } pair64_t;
  19
  20 #define pair64_lt(a,b) ((a).u < (b).u)
  21 KSORT_INIT(off, pair64_t, pair64_lt)
  22
  23 typedef struct {
  24         uint32_t m, n;
  25         pair64_t *list;
  26 } ti_binlist_t;
  27
  28 typedef struct {
  29         int32_t n, m;
  30         uint64_t *offset;
  31 } ti_lidx_t;
  32
  33 KHASH_MAP_INIT_INT(i, ti_binlist_t)
  34 KHASH_MAP_INIT_STR(s, int)
  35
  36 struct __ti_index_t {
  37         ti_conf_t conf;
  38         int32_t n, max;
  39         khash_t(s) *tname;
  40         khash_t(i) **index;
  41         ti_lidx_t *index2;
  42 };
  43
  44 typedef struct {
  45         int tid, beg, end, bin;
  46 } ti_intv_t;
  47
  48 ti_conf_t ti_conf_gff = { 0, 1, 4, 5, '#', 0 };
  49 ti_conf_t ti_conf_bed = { TI_FLAG_UCSC, 1, 2, 3, '#', 0 };
  50 ti_conf_t ti_conf_psltbl = { TI_FLAG_UCSC, 15, 17, 18, '#', 0 };
  51 ti_conf_t ti_conf_sam = { TI_PRESET_SAM, 3, 4, 0, '@', 0 };
  52 ti_conf_t ti_conf_vcf = { TI_PRESET_VCF, 1, 2, 0, '#', 0 };
  53
  54 /***************
  55  * read a line *
  56  ***************/
  57
  58 /*
  59 int ti_readline(BGZF *fp, kstring_t *str)
  60 {
  61         int c, l = 0;
  62         str->l = 0;
  63         while ((c = bgzf_getc(fp)) >= 0 && c != '\n') {
  64                 ++l;
  65                 if (c != '\r') kputc(c, str);
  66         }
  67         if (c < 0 && l == 0) return -1; // end of file
  68         return str->l;
  69 }
  70 */
  71
  72 /* Below is a faster implementation largely equivalent to the one
  73  * commented out above. */
  74 int ti_readline(BGZF *fp, kstring_t *str)
  75 {
  76         int l, state = 0;
  77         unsigned char *buf = (unsigned char*)fp->uncompressed_block;
  78         str->l = 0;
  79         do {
  80                 if (fp->block_offset >= fp->block_length) {
  81                         if (bgzf_read_block(fp) != 0) { state = -2; break; }
  82                         if (fp->block_length == 0) { state = -1; break; }
  83                 }
  84                 for (l = fp->block_offset; l < fp->block_length && buf[l] != '\n'; ++l);
  85                 if (l < fp->block_length) state = 1;
  86                 l -= fp->block_offset;
  87                 if (str->l + l + 1 >= str->m) {
  88                         str->m = str->l + l + 2;
  89                         kroundup32(str->m);
  90                         str->s = (char*)realloc(str->s, str->m);
  91                 }
  92                 memcpy(str->s + str->l, buf + fp->block_offset, l);
  93                 str->l += l;
  94                 fp->block_offset += l + 1;
  95                 if (fp->block_offset >= fp->block_length) {
  96 #ifdef _USE_KNETFILE
  97                         fp->block_address = knet_tell(fp->x.fpr);
  98 #else
  99                         fp->block_address = ftello(fp->file);
 100 #endif
 101                         fp->block_offset = 0;
 102                         fp->block_length = 0;
 103                 }
 104         } while (state == 0);
 105         if (str->l == 0 && state < 0) return state;
 106         str->s[str->l] = 0;
 107         return str->l;
 108 }
 109
 110 /*************************************
 111  * get the interval from a data line *
 112  *************************************/
 113
 114 static inline int ti_reg2bin(uint32_t beg, uint32_t end)
 115 {
 116         --end;
 117         if (beg>>14 == end>>14) return 4681 + (beg>>14);
 118         if (beg>>17 == end>>17) return  585 + (beg>>17);
 119         if (beg>>20 == end>>20) return   73 + (beg>>20);
 120         if (beg>>23 == end>>23) return    9 + (beg>>23);
 121         if (beg>>26 == end>>26) return    1 + (beg>>26);
 122         return 0;
 123 }
 124
 125 static int get_tid(ti_index_t *idx, const char *ss)
 126 {
 127         khint_t k;
 128         int tid;
 129         k = kh_get(s, idx->tname, ss);
 130         if (k == kh_end(idx->tname)) { // a new target sequence
 131                 int ret, size;
 132                 // update idx->n, ->max, ->index and ->index2
 133                 if (idx->n == idx->max) {
 134                         idx->max = idx->max? idx->max<<1 : 8;
 135                         idx->index = realloc(idx->index, idx->max * sizeof(void*));
 136                         idx->index2 = realloc(idx->index2, idx->max * sizeof(ti_lidx_t));
 137                 }
 138                 memset(&idx->index2[idx->n], 0, sizeof(ti_lidx_t));
 139                 idx->index[idx->n++] = kh_init(i);
 140                 // update ->tname
 141                 tid = size = kh_size(idx->tname);
 142                 k = kh_put(s, idx->tname, strdup(ss), &ret);
 143                 kh_value(idx->tname, k) = size;
 144                 assert(idx->n == kh_size(idx->tname));
 145         } else tid = kh_value(idx->tname, k);
 146         return tid;
 147 }
 148
 149 static int get_intv(ti_index_t *idx, kstring_t *str, ti_intv_t *intv)
 150 {
 151         int i, b = 0, id = 1;
 152         char *s;
 153         intv->tid = intv->beg = intv->end = intv->bin = -1;
 154         for (i = 0; i <= str->l; ++i) {
 155                 if (str->s[i] == '\t' || str->s[i] == 0) {
 156                         if (id == idx->conf.sc) {
 157                                 str->s[i] = 0;
 158                                 intv->tid = get_tid(idx, str->s + b);
 159                                 if (i != str->l) str->s[i] = '\t';
 160                         } else if (id == idx->conf.bc) {
 161                                 // here ->beg is 0-based.
 162                                 intv->beg = intv->end = strtol(str->s + b, &s, 0);
 163                                 if (!(idx->conf.preset&TI_FLAG_UCSC)) --intv->beg;
 164                         } else {
 165                                 if ((idx->conf.preset&0xffff) == TI_PRESET_GENERIC) {
 166                                         if (id == idx->conf.ec) intv->end = strtol(str->s + b, &s, 0);
 167                                 } else if ((idx->conf.preset&0xffff) == TI_PRESET_SAM) {
 168                                         if (id == 6) { // CIGAR
 169                                                 int l = 0, op;
 170                                                 char *t;
 171                                                 for (s = str->s + b; s < str->s + i;) {
 172                                                         long x = strtol(s, &t, 10);
 173                                                         op = toupper(*t);
 174                                                         if (op == 'M' || op == 'D' || op == 'N') l += x;
 175                                                         s = t + 1;
 176                                                 }
 177                                                 intv->end = intv->beg + l;
 178                                         }
 179                                 } else if ((idx->conf.preset&0xffff) == TI_PRESET_VCF) {
 180                                         // FIXME: the following is NOT tested and is likely to be buggy
 181                                         if (id == 5) { // ALT
 182                                                 char *t;
 183                                                 int max = 1;
 184                                                 for (s = str->s + b; s < str->s + i;) {
 185                                                         if (s[i] == 'D') {
 186                                                                 long x = strtol(s + 1, &t, 10);
 187                                                                 if (x > max) max = x;
 188                                                                 s = t + 1;
 189                                                         } else ++s;
 190                                                 }
 191                                                 intv->end = intv->beg + max;
 192                                         }
 193                                 }
 194                         }
 195                         b = i + 1;
 196                         ++id;
 197                 }
 198         }
 199         if (intv->tid < 0 || intv->beg < 0 || intv->end < 0) return -1;
 200         intv->bin = ti_reg2bin(intv->beg, intv->end);
 201         return 0;
 202 }
 203
 204 /************
 205  * indexing *
 206  ************/
 207
 208 // requirement: len <= LEN_MASK
 209 static inline void insert_offset(khash_t(i) *h, int bin, uint64_t beg, uint64_t end)
 210 {
 211         khint_t k;
 212         ti_binlist_t *l;
 213         int ret;
 214         k = kh_put(i, h, bin, &ret);
 215         l = &kh_value(h, k);
 216         if (ret) { // not present
 217                 l->m = 1; l->n = 0;
 218                 l->list = (pair64_t*)calloc(l->m, 16);
 219         }
 220         if (l->n == l->m) {
 221                 l->m <<= 1;
 222                 l->list = (pair64_t*)realloc(l->list, l->m * 16);
 223         }
 224         l->list[l->n].u = beg; l->list[l->n++].v = end;
 225 }
 226
 227 static inline void insert_offset2(ti_lidx_t *index2, int _beg, int _end, uint64_t offset)
 228 {
 229         int i, beg, end;
 230         beg = _beg >> TAD_LIDX_SHIFT;
 231         end = (_end - 1) >> TAD_LIDX_SHIFT;
 232         if (index2->m < end + 1) {
 233                 int old_m = index2->m;
 234                 index2->m = end + 1;
 235                 kroundup32(index2->m);
 236                 index2->offset = (uint64_t*)realloc(index2->offset, index2->m * 8);
 237                 memset(index2->offset + old_m, 0, 8 * (index2->m - old_m));
 238         }
 239         for (i = beg + 1; i <= end; ++i)
 240                 if (index2->offset[i] == 0) index2->offset[i] = offset;
 241         index2->n = end + 1;
 242 }
 243
 244 static void merge_chunks(ti_index_t *idx)
 245 {
 246         khash_t(i) *index;
 247         int i, l, m;
 248         khint_t k;
 249         for (i = 0; i < idx->n; ++i) {
 250                 index = idx->index[i];
 251                 for (k = kh_begin(index); k != kh_end(index); ++k) {
 252                         ti_binlist_t *p;
 253                         if (!kh_exist(index, k)) continue;
 254                         p = &kh_value(index, k);
 255                         m = 0;
 256                         for (l = 1; l < p->n; ++l) {
 257                                 if (p->list[m].v>>16 == p->list[l].u>>16) p->list[m].v = p->list[l].v;
 258                                 else p->list[++m] = p->list[l];
 259                         } // ~for(l)
 260                         p->n = m + 1;
 261                 } // ~for(k)
 262         } // ~for(i)
 263 }
 264
 265 ti_index_t *ti_index_core(BGZF *fp, const ti_conf_t *conf)
 266 {
 267         int ret;
 268         ti_index_t *idx;
 269         uint32_t last_bin, save_bin;
 270         int32_t last_coor, last_tid, save_tid;
 271         uint64_t save_off, last_off, lineno = 0;
 272         kstring_t *str;
 273
 274         str = calloc(1, sizeof(kstring_t));
 275
 276         idx = (ti_index_t*)calloc(1, sizeof(ti_index_t));
 277         idx->conf = *conf;
 278         idx->n = idx->max = 0;
 279         idx->tname = kh_init(s);
 280         idx->index = 0;
 281         idx->index2 = 0;
 282
 283         save_bin = save_tid = last_tid = last_bin = 0xffffffffu;
 284         save_off = last_off = bgzf_tell(fp); last_coor = 0xffffffffu;
 285         while ((ret = ti_readline(fp, str)) >= 0) {
 286                 ti_intv_t intv;
 287                 ++lineno;
 288                 if (lineno <= idx->conf.line_skip || str->s[0] == idx->conf.meta_char) {
 289                         last_off = bgzf_tell(fp);
 290                         continue;
 291                 }
 292                 get_intv(idx, str, &intv);
 293                 if (last_tid != intv.tid) { // change of chromosomes
 294                         last_tid = intv.tid;
 295                         last_bin = 0xffffffffu;
 296                 } else if (last_coor > intv.beg) {
 297                         fprintf(stderr, "[ti_index_core] the file is not sorted at line %llu\n", (unsigned long long)lineno);
 298                         exit(1);
 299                 }
 300                 if (intv.bin < 4681) insert_offset2(&idx->index2[intv.tid], intv.beg, intv.end, last_off);
 301                 if (intv.bin != last_bin) { // then possibly write the binning index
 302                         if (save_bin != 0xffffffffu) // save_bin==0xffffffffu only happens to the first record
 303                                 insert_offset(idx->index[save_tid], save_bin, save_off, last_off);
 304                         save_off = last_off;
 305                         save_bin = last_bin = intv.bin;
 306                         save_tid = intv.tid;
 307                         if (save_tid < 0) break;
 308                 }
 309                 if (bgzf_tell(fp) <= last_off) {
 310                         fprintf(stderr, "[ti_index_core] bug in BGZF: %llx < %llx\n",
 311                                         (unsigned long long)bgzf_tell(fp), (unsigned long long)last_off);
 312                         exit(1);
 313                 }
 314                 last_off = bgzf_tell(fp);
 315                 last_coor = intv.beg;
 316         }
 317         if (save_tid >= 0) insert_offset(idx->index[save_tid], save_bin, save_off, bgzf_tell(fp));
 318         merge_chunks(idx);
 319
 320         free(str->s); free(str);
 321         return idx;
 322 }
 323
 324 void ti_index_destroy(ti_index_t *idx)
 325 {
 326         khint_t k;
 327         int i;
 328         if (idx == 0) return;
 329         // destroy the name hash table
 330         for (k = kh_begin(idx->tname); k != kh_end(idx->tname); ++k) {
 331                 if (kh_exist(idx->tname, k))
 332                         free((char*)kh_key(idx->tname, k));
 333         }
 334         kh_destroy(s, idx->tname);
 335         // destroy the binning index
 336         for (i = 0; i < idx->n; ++i) {
 337                 khash_t(i) *index = idx->index[i];
 338                 ti_lidx_t *index2 = idx->index2 + i;
 339                 for (k = kh_begin(index); k != kh_end(index); ++k) {
 340                         if (kh_exist(index, k))
 341                                 free(kh_value(index, k).list);
 342                 }
 343                 kh_destroy(i, index);
 344                 free(index2->offset);
 345         }
 346         free(idx->index);
 347         // destroy the linear index
 348         free(idx->index2);
 349         free(idx);
 350 }
 351
 352 /******************
 353  * index file I/O *
 354  ******************/
 355
 356 void ti_index_save(const ti_index_t *idx, BGZF *fp)
 357 {
 358         int32_t i, size, ti_is_be;
 359         khint_t k;
 360         ti_is_be = bam_is_big_endian();
 361         bgzf_write(fp, "TBI\1", 4);
 362         if (ti_is_be) {
 363                 uint32_t x = idx->n;
 364                 bgzf_write(fp, bam_swap_endian_4p(&x), 4);
 365         } else bgzf_write(fp, &idx->n, 4);
 366         assert(sizeof(ti_conf_t) == 24);
 367         if (ti_is_be) { // write ti_conf_t;
 368                 uint32_t x[6];
 369                 memcpy(x, &idx->conf, 24);
 370                 for (i = 0; i < 6; ++i) bgzf_write(fp, bam_swap_endian_4p(&x[i]), 4);
 371         } else bgzf_write(fp, &idx->conf, sizeof(ti_conf_t));
 372         { // write target names
 373                 char **name;
 374                 int32_t l = 0;
 375                 name = calloc(kh_size(idx->tname), sizeof(void*));
 376                 for (k = kh_begin(idx->tname); k != kh_end(idx->tname); ++k)
 377                         if (kh_exist(idx->tname, k))
 378                                 name[kh_value(idx->tname, k)] = (char*)kh_key(idx->tname, k);
 379                 for (i = 0; i < kh_size(idx->tname); ++i)
 380                         l += strlen(name[i]) + 1;
 381                 if (ti_is_be) bgzf_write(fp, bam_swap_endian_4p(&l), 4);
 382                 else bgzf_write(fp, &l, 4);
 383                 for (i = 0; i < kh_size(idx->tname); ++i)
 384                         bgzf_write(fp, name[i], strlen(name[i]) + 1);
 385                 free(name);
 386         }
 387         for (i = 0; i < idx->n; ++i) {
 388                 khash_t(i) *index = idx->index[i];
 389                 ti_lidx_t *index2 = idx->index2 + i;
 390                 // write binning index
 391                 size = kh_size(index);
 392                 if (ti_is_be) { // big endian
 393                         uint32_t x = size;
 394                         bgzf_write(fp, bam_swap_endian_4p(&x), 4);
 395                 } else bgzf_write(fp, &size, 4);
 396                 for (k = kh_begin(index); k != kh_end(index); ++k) {
 397                         if (kh_exist(index, k)) {
 398                                 ti_binlist_t *p = &kh_value(index, k);
 399                                 if (ti_is_be) { // big endian
 400                                         uint32_t x;
 401                                         x = kh_key(index, k); bgzf_write(fp, bam_swap_endian_4p(&x), 4);
 402                                         x = p->n; bgzf_write(fp, bam_swap_endian_4p(&x), 4);
 403                                         for (x = 0; (int)x < p->n; ++x) {
 404                                                 bam_swap_endian_8p(&p->list[x].u);
 405                                                 bam_swap_endian_8p(&p->list[x].v);
 406                                         }
 407                                         bgzf_write(fp, p->list, 16 * p->n);
 408                                         for (x = 0; (int)x < p->n; ++x) {
 409                                                 bam_swap_endian_8p(&p->list[x].u);
 410                                                 bam_swap_endian_8p(&p->list[x].v);
 411                                         }
 412                                 } else {
 413                                         bgzf_write(fp, &kh_key(index, k), 4);
 414                                         bgzf_write(fp, &p->n, 4);
 415                                         bgzf_write(fp, p->list, 16 * p->n);
 416                                 }
 417                         }
 418                 }
 419                 // write linear index (index2)
 420                 if (ti_is_be) {
 421                         int x = index2->n;
 422                         bgzf_write(fp, bam_swap_endian_4p(&x), 4);
 423                 } else bgzf_write(fp, &index2->n, 4);
 424                 if (ti_is_be) { // big endian
 425                         int x;
 426                         for (x = 0; (int)x < index2->n; ++x)
 427                                 bam_swap_endian_8p(&index2->offset[x]);
 428                         bgzf_write(fp, index2->offset, 8 * index2->n);
 429                         for (x = 0; (int)x < index2->n; ++x)
 430                                 bam_swap_endian_8p(&index2->offset[x]);
 431                 } else bgzf_write(fp, index2->offset, 8 * index2->n);
 432         }
 433 }
 434
 435 static ti_index_t *ti_index_load_core(BGZF *fp)
 436 {
 437         int i, ti_is_be;
 438         char magic[4];
 439         ti_index_t *idx;
 440         ti_is_be = bam_is_big_endian();
 441         if (fp == 0) {
 442                 fprintf(stderr, "[ti_index_load_core] fail to load index.\n");
 443                 return 0;
 444         }
 445         bgzf_read(fp, magic, 4);
 446         if (strncmp(magic, "TBI\1", 4)) {
 447                 fprintf(stderr, "[ti_index_load] wrong magic number.\n");
 448                 return 0;
 449         }
 450         idx = (ti_index_t*)calloc(1, sizeof(ti_index_t));
 451         bgzf_read(fp, &idx->n, 4);
 452         if (ti_is_be) bam_swap_endian_4p(&idx->n);
 453         idx->tname = kh_init(s);
 454         idx->index = (khash_t(i)**)calloc(idx->n, sizeof(void*));
 455         idx->index2 = (ti_lidx_t*)calloc(idx->n, sizeof(ti_lidx_t));
 456         // read idx->conf
 457         bgzf_read(fp, &idx->conf, sizeof(ti_conf_t));
 458         if (ti_is_be) {
 459                 bam_swap_endian_4p(&idx->conf.preset);
 460                 bam_swap_endian_4p(&idx->conf.sc);
 461                 bam_swap_endian_4p(&idx->conf.bc);
 462                 bam_swap_endian_4p(&idx->conf.ec);
 463                 bam_swap_endian_4p(&idx->conf.meta_char);
 464                 bam_swap_endian_4p(&idx->conf.line_skip);
 465         }
 466         { // read target names
 467                 int j, ret;
 468                 kstring_t *str;
 469                 int32_t l;
 470                 uint8_t *buf;
 471                 bgzf_read(fp, &l, 4);
 472                 if (ti_is_be) bam_swap_endian_4p(&l);
 473                 buf = calloc(l, 1);
 474                 bgzf_read(fp, buf, l);
 475                 str = calloc(1, sizeof(kstring_t));
 476                 for (i = j = 0; i < l; ++i) {
 477                         if (buf[i] == 0) {
 478                                 khint_t k = kh_put(s, idx->tname, strdup(str->s), &ret);
 479                                 kh_value(idx->tname, k) = j++;
 480                                 str->l = 0;
 481                         } else kputc(buf[i], str);
 482                 }
 483                 free(str->s); free(str); free(buf);
 484         }
 485         for (i = 0; i < idx->n; ++i) {
 486                 khash_t(i) *index;
 487                 ti_lidx_t *index2 = idx->index2 + i;
 488                 uint32_t key, size;
 489                 khint_t k;
 490                 int j, ret;
 491                 ti_binlist_t *p;
 492                 index = idx->index[i] = kh_init(i);
 493                 // load binning index
 494                 bgzf_read(fp, &size, 4);
 495                 if (ti_is_be) bam_swap_endian_4p(&size);
 496                 for (j = 0; j < (int)size; ++j) {
 497                         bgzf_read(fp, &key, 4);
 498                         if (ti_is_be) bam_swap_endian_4p(&key);
 499                         k = kh_put(i, index, key, &ret);
 500                         p = &kh_value(index, k);
 501                         bgzf_read(fp, &p->n, 4);
 502                         if (ti_is_be) bam_swap_endian_4p(&p->n);
 503                         p->m = p->n;
 504                         p->list = (pair64_t*)malloc(p->m * 16);
 505                         bgzf_read(fp, p->list, 16 * p->n);
 506                         if (ti_is_be) {
 507                                 int x;
 508                                 for (x = 0; x < p->n; ++x) {
 509                                         bam_swap_endian_8p(&p->list[x].u);
 510                                         bam_swap_endian_8p(&p->list[x].v);
 511                                 }
 512                         }
 513                 }
 514                 // load linear index
 515                 bgzf_read(fp, &index2->n, 4);
 516                 if (ti_is_be) bam_swap_endian_4p(&index2->n);
 517                 index2->m = index2->n;
 518                 index2->offset = (uint64_t*)calloc(index2->m, 8);
 519                 bgzf_read(fp, index2->offset, index2->n * 8);
 520                 if (ti_is_be)
 521                         for (j = 0; j < index2->n; ++j) bam_swap_endian_8p(&index2->offset[j]);
 522         }
 523         return idx;
 524 }
 525
 526 ti_index_t *ti_index_load_local(const char *_fn)
 527 {
 528         BGZF *fp;
 529         char *fnidx, *fn;
 530
 531         if (strstr(_fn, "ftp://") == _fn || strstr(_fn, "http://") == _fn) {
 532                 const char *p;
 533                 int l = strlen(_fn);
 534                 for (p = _fn + l - 1; p >= _fn; --p)
 535                         if (*p == '/') break;
 536                 fn = strdup(p + 1);
 537         } else fn = strdup(_fn);
 538         fnidx = (char*)calloc(strlen(fn) + 5, 1);
 539         strcpy(fnidx, fn); strcat(fnidx, ".tbi");
 540         fp = bgzf_open(fnidx, "r");
 541         free(fnidx); free(fn);
 542         if (fp) {
 543                 ti_index_t *idx = ti_index_load_core(fp);
 544                 bgzf_close(fp);
 545                 return idx;
 546         } else return 0;
 547 }
 548
 549 #ifdef _USE_KNETFILE
 550 static void download_from_remote(const char *url)
 551 {
 552         const int buf_size = 1 * 1024 * 1024;
 553         char *fn;
 554         FILE *fp;
 555         uint8_t *buf;
 556         knetFile *fp_remote;
 557         int l;
 558         if (strstr(url, "ftp://") != url && strstr(url, "http://") != url) return;
 559         l = strlen(url);
 560         for (fn = (char*)url + l - 1; fn >= url; --fn)
 561                 if (*fn == '/') break;
 562         ++fn; // fn now points to the file name
 563         fp_remote = knet_open(url, "r");
 564         if (fp_remote == 0) {
 565                 fprintf(stderr, "[download_from_remote] fail to open remote file.\n");
 566                 return;
 567         }
 568         if ((fp = fopen(fn, "w")) == 0) {
 569                 fprintf(stderr, "[download_from_remote] fail to create file in the working directory.\n");
 570                 knet_close(fp_remote);
 571                 return;
 572         }
 573         buf = (uint8_t*)calloc(buf_size, 1);
 574         while ((l = knet_read(fp_remote, buf, buf_size)) != 0)
 575                 fwrite(buf, 1, l, fp);
 576         free(buf);
 577         fclose(fp);
 578         knet_close(fp_remote);
 579 }
 580 #else
 581 static void download_from_remote(const char *url)
 582 {
 583         return;
 584 }
 585 #endif
 586
 587 ti_index_t *ti_index_load(const char *fn)
 588 {
 589         ti_index_t *idx;
 590         idx = ti_index_load_local(fn);
 591         if (idx == 0 && (strstr(fn, "ftp://") == fn || strstr(fn, "http://") == fn)) {
 592                 char *fnidx = calloc(strlen(fn) + 5, 1);
 593                 strcat(strcpy(fnidx, fn), ".tbi");
 594                 fprintf(stderr, "[ti_index_load] attempting to download the remote index file.\n");
 595                 download_from_remote(fnidx);
 596                 idx = ti_index_load_local(fn);
 597         }
 598         if (idx == 0) fprintf(stderr, "[ti_index_load] fail to load BAM index.\n");
 599         return idx;
 600 }
 601
 602 int ti_index_build2(const char *fn, const ti_conf_t *conf, const char *_fnidx)
 603 {
 604         char *fnidx;
 605         BGZF *fp, *fpidx;
 606         ti_index_t *idx;
 607         if ((fp = bgzf_open(fn, "r")) == 0) {
 608                 fprintf(stderr, "[ti_index_build2] fail to open the BAM file.\n");
 609                 return -1;
 610         }
 611         idx = ti_index_core(fp, conf);
 612         bgzf_close(fp);
 613         if (_fnidx == 0) {
 614                 fnidx = (char*)calloc(strlen(fn) + 5, 1);
 615                 strcpy(fnidx, fn); strcat(fnidx, ".tbi");
 616         } else fnidx = strdup(_fnidx);
 617         fpidx = bgzf_open(fnidx, "w");
 618         if (fpidx == 0) {
 619                 fprintf(stderr, "[ti_index_build2] fail to create the index file.\n");
 620                 free(fnidx);
 621                 return -1;
 622         }
 623         ti_index_save(idx, fpidx);
 624         ti_index_destroy(idx);
 625         bgzf_close(fpidx);
 626         free(fnidx);
 627         return 0;
 628 }
 629
 630 int ti_index_build(const char *fn, const ti_conf_t *conf)
 631 {
 632         return ti_index_build2(fn, conf, 0);
 633 }
 634
 635 /********************************************
 636  * parse a region in the format chr:beg-end *
 637  ********************************************/
 638
 639 int ti_parse_region(ti_index_t *idx, const char *str, int *tid, int *begin, int *end)
 640 {
 641         char *s, *p;
 642         int i, l, k;
 643         khiter_t iter;
 644         khash_t(s) *h = idx->tname;
 645         l = strlen(str);
 646         p = s = (char*)malloc(l+1);
 647         /* squeeze out "," */
 648         for (i = k = 0; i != l; ++i)
 649                 if (str[i] != ',' && !isspace(str[i])) s[k++] = str[i];
 650         s[k] = 0;
 651         for (i = 0; i != k; ++i) if (s[i] == ':') break;
 652         s[i] = 0;
 653         iter = kh_get(s, h, s); /* get the tid */
 654         if (iter == kh_end(h)) { // name not found
 655                 *tid = -1; free(s);
 656                 return -1;
 657         }
 658         *tid = kh_value(h, iter);
 659         if (i == k) { /* dump the whole sequence */
 660                 *begin = 0; *end = 1<<29; free(s);
 661                 return 0;
 662         }
 663         for (p = s + i + 1; i != k; ++i) if (s[i] == '-') break;
 664         *begin = atoi(p);
 665         if (i < k) {
 666                 p = s + i + 1;
 667                 *end = atoi(p);
 668         } else *end = 1<<29;
 669         if (*begin > 0) --*begin;
 670         free(s);
 671         if (*begin > *end) return -1;
 672         return 0;
 673 }
 674
 675 /*******************************
 676  * retrieve a specified region *
 677  *******************************/
 678
 679 #define MAX_BIN 37450 // =(8^6-1)/7+1
 680
 681 static inline int reg2bins(uint32_t beg, uint32_t end, uint16_t list[MAX_BIN])
 682 {
 683         int i = 0, k;
 684         --end;
 685         list[i++] = 0;
 686         for (k =    1 + (beg>>26); k <=    1 + (end>>26); ++k) list[i++] = k;
 687         for (k =    9 + (beg>>23); k <=    9 + (end>>23); ++k) list[i++] = k;
 688         for (k =   73 + (beg>>20); k <=   73 + (end>>20); ++k) list[i++] = k;
 689         for (k =  585 + (beg>>17); k <=  585 + (end>>17); ++k) list[i++] = k;
 690         for (k = 4681 + (beg>>14); k <= 4681 + (end>>14); ++k) list[i++] = k;
 691         return i;
 692 }
 693
 694 static inline int is_overlap(int32_t beg, int32_t end, int rbeg, int rend)
 695 {
 696         return (rend > beg && rbeg < end);
 697 }
 698
 699 // ti_fetch helper function retrieves
 700 pair64_t *get_chunk_coordinates(const ti_index_t *idx, int tid, int beg, int end, int* cnt_off)
 701 {
 702         uint16_t *bins;
 703         int i, n_bins, n_off;
 704         pair64_t *off;
 705         khint_t k;
 706         khash_t(i) *index;
 707         uint64_t min_off;
 708
 709         bins = (uint16_t*)calloc(MAX_BIN, 2);
 710         n_bins = reg2bins(beg, end, bins);
 711         index = idx->index[tid];
 712         min_off = (beg>>TAD_LIDX_SHIFT >= idx->index2[tid].n)? 0 : idx->index2[tid].offset[beg>>TAD_LIDX_SHIFT];
 713         for (i = n_off = 0; i < n_bins; ++i) {
 714                 if ((k = kh_get(i, index, bins[i])) != kh_end(index))
 715                         n_off += kh_value(index, k).n;
 716         }
 717         if (n_off == 0) {
 718                 free(bins); return 0;
 719         }
 720         off = (pair64_t*)calloc(n_off, 16);
 721         for (i = n_off = 0; i < n_bins; ++i) {
 722                 if ((k = kh_get(i, index, bins[i])) != kh_end(index)) {
 723                         int j;
 724                         ti_binlist_t *p = &kh_value(index, k);
 725                         for (j = 0; j < p->n; ++j)
 726                                 if (p->list[j].v > min_off) off[n_off++] = p->list[j];
 727                 }
 728         }
 729         free(bins);
 730         {
 731                 int l;
 732                 ks_introsort(off, n_off, off);
 733                 // resolve completely contained adjacent blocks
 734                 for (i = 1, l = 0; i < n_off; ++i)
 735                         if (off[l].v < off[i].v)
 736                                 off[++l] = off[i];
 737                 n_off = l + 1;
 738                 // resolve overlaps between adjacent blocks; this may happen due to the merge in indexing
 739                 for (i = 1; i < n_off; ++i)
 740                         if (off[i-1].v >= off[i].u) off[i-1].v = off[i].u;
 741                 { // merge adjacent blocks
 742                         for (i = 1, l = 0; i < n_off; ++i) {
 743                                 if (off[l].v>>16 == off[i].u>>16) off[l].v = off[i].v;
 744                                 else off[++l] = off[i];
 745                         }
 746                         n_off = l + 1;
 747                 }
 748         }
 749         *cnt_off = n_off;
 750         return off;
 751 }
 752
 753 int ti_fetch(BGZF *fp, const ti_index_t *idx, int tid, int beg, int end, void *data, ti_fetch_f func)
 754 {
 755         int n_off;
 756         pair64_t *off = get_chunk_coordinates(idx, tid, beg, end, &n_off);
 757         if (off == 0) return 0;
 758         {
 759                 // retrive alignments
 760                 uint64_t curr_off;
 761                 int i, ret, n_seeks;
 762                 kstring_t *str = calloc(1, sizeof(kstring_t));
 763                 n_seeks = 0; i = -1; curr_off = 0;
 764                 for (;;) {
 765                         if (curr_off == 0 || curr_off >= off[i].v) { // then jump to the next chunk
 766                                 if (i == n_off - 1) break; // no more chunks
 767                                 if (i >= 0) assert(curr_off == off[i].v); // otherwise bug
 768                                 if (i < 0 || off[i].v != off[i+1].u) { // not adjacent chunks; then seek
 769                                         bgzf_seek(fp, off[i+1].u, SEEK_SET);
 770                                         curr_off = bgzf_tell(fp);
 771                                         ++n_seeks;
 772                                 }
 773                                 ++i;
 774                         }
 775                         if ((ret = ti_readline(fp, str)) >= 0) {
 776                                 ti_intv_t intv;
 777                                 curr_off = bgzf_tell(fp);
 778                                 if (str->s[0] == idx->conf.meta_char) continue;
 779                                 get_intv((ti_index_t*)idx, str, &intv);
 780                                 if (intv.tid != tid || intv.beg >= end) break; // no need to proceed
 781                                 else if (is_overlap(beg, end, intv.beg, intv.end)) func(str->l, str->s, data);
 782                         } else break; // end of file
 783                 }
 784 //              fprintf(stderr, "[ti_fetch] # seek calls: %d\n", n_seeks);
 785                 free(str->s); free(str);
 786         }
 787         free(off);
 788         return 0;
 789 }