index.c

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