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