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