6 #include "bam_endian.h"
12 #define TAD_MIN_CHUNK_GAP 32768
13 // 1<<14 is the size of minimum bin.
14 #define TAD_LIDX_SHIFT 14
20 #define pair64_lt(a,b) ((a).u < (b).u)
21 KSORT_INIT(off, pair64_t, pair64_lt)
33 KHASH_MAP_INIT_INT(i, ti_binlist_t)
34 KHASH_MAP_INIT_STR(s, int)
45 int tid, beg, end, bin;
48 ti_conf_t ti_conf_gff = { 0, 1, 4, 5, '#', 0 };
49 ti_conf_t ti_conf_bed = { TI_FLAG_UCSC, 1, 2, 3, '#', 0 };
50 ti_conf_t ti_conf_psltbl = { TI_FLAG_UCSC, 15, 17, 18, '#', 0 };
51 ti_conf_t ti_conf_sam = { TI_PRESET_SAM, 3, 4, 0, '@', 0 };
52 ti_conf_t ti_conf_vcf = { TI_PRESET_VCF, 1, 2, 0, '#', 0 };
59 int ti_readline(BGZF *fp, kstring_t *str)
63 while ((c = bgzf_getc(fp)) >= 0 && c != '\n') {
65 if (c != '\r') kputc(c, str);
67 if (c < 0 && l == 0) return -1; // end of file
72 /* Below is a faster implementation largely equivalent to the one
73 * commented out above. */
74 int ti_readline(BGZF *fp, kstring_t *str)
77 unsigned char *buf = (unsigned char*)fp->uncompressed_block;
80 if (fp->block_offset >= fp->block_length) {
81 if (bgzf_read_block(fp) != 0) { state = -2; break; }
82 if (fp->block_length == 0) { state = -1; break; }
84 for (l = fp->block_offset; l < fp->block_length && buf[l] != '\n'; ++l);
85 if (l < fp->block_length) state = 1;
86 l -= fp->block_offset;
87 if (str->l + l + 1 >= str->m) {
88 str->m = str->l + l + 2;
90 str->s = (char*)realloc(str->s, str->m);
92 memcpy(str->s + str->l, buf + fp->block_offset, l);
94 fp->block_offset += l + 1;
95 if (fp->block_offset >= fp->block_length) {
97 fp->block_address = knet_tell(fp->x.fpr);
99 fp->block_address = ftello(fp->file);
101 fp->block_offset = 0;
102 fp->block_length = 0;
104 } while (state == 0);
105 if (str->l == 0 && state < 0) return state;
110 /*************************************
111 * get the interval from a data line *
112 *************************************/
114 static inline int ti_reg2bin(uint32_t beg, uint32_t end)
117 if (beg>>14 == end>>14) return 4681 + (beg>>14);
118 if (beg>>17 == end>>17) return 585 + (beg>>17);
119 if (beg>>20 == end>>20) return 73 + (beg>>20);
120 if (beg>>23 == end>>23) return 9 + (beg>>23);
121 if (beg>>26 == end>>26) return 1 + (beg>>26);
125 static int get_tid(ti_index_t *idx, const char *ss)
129 k = kh_get(s, idx->tname, ss);
130 if (k == kh_end(idx->tname)) { // a new target sequence
132 // update idx->n, ->max, ->index and ->index2
133 if (idx->n == idx->max) {
134 idx->max = idx->max? idx->max<<1 : 8;
135 idx->index = realloc(idx->index, idx->max * sizeof(void*));
136 idx->index2 = realloc(idx->index2, idx->max * sizeof(ti_lidx_t));
138 memset(&idx->index2[idx->n], 0, sizeof(ti_lidx_t));
139 idx->index[idx->n++] = kh_init(i);
141 tid = size = kh_size(idx->tname);
142 k = kh_put(s, idx->tname, strdup(ss), &ret);
143 kh_value(idx->tname, k) = size;
144 assert(idx->n == kh_size(idx->tname));
145 } else tid = kh_value(idx->tname, k);
149 static int get_intv(ti_index_t *idx, kstring_t *str, ti_intv_t *intv)
151 int i, b = 0, id = 1;
153 intv->tid = intv->beg = intv->end = intv->bin = -1;
154 for (i = 0; i <= str->l; ++i) {
155 if (str->s[i] == '\t' || str->s[i] == 0) {
156 if (id == idx->conf.sc) {
158 intv->tid = get_tid(idx, str->s + b);
159 if (i != str->l) str->s[i] = '\t';
160 } else if (id == idx->conf.bc) {
161 // here ->beg is 0-based.
162 intv->beg = intv->end = strtol(str->s + b, &s, 0);
163 if (!(idx->conf.preset&TI_FLAG_UCSC)) --intv->beg;
165 if ((idx->conf.preset&0xffff) == TI_PRESET_GENERIC) {
166 if (id == idx->conf.ec) intv->end = strtol(str->s + b, &s, 0);
167 } else if ((idx->conf.preset&0xffff) == TI_PRESET_SAM) {
168 if (id == 6) { // CIGAR
171 for (s = str->s + b; s < str->s + i;) {
172 long x = strtol(s, &t, 10);
174 if (op == 'M' || op == 'D' || op == 'N') l += x;
177 intv->end = intv->beg + l;
179 } else if ((idx->conf.preset&0xffff) == TI_PRESET_VCF) {
180 // FIXME: the following is NOT tested and is likely to be buggy
181 if (id == 5) { // ALT
184 for (s = str->s + b; s < str->s + i;) {
186 long x = strtol(s + 1, &t, 10);
187 if (x > max) max = x;
191 intv->end = intv->beg + max;
199 if (intv->tid < 0 || intv->beg < 0 || intv->end < 0) return -1;
200 intv->bin = ti_reg2bin(intv->beg, intv->end);
208 // requirement: len <= LEN_MASK
209 static inline void insert_offset(khash_t(i) *h, int bin, uint64_t beg, uint64_t end)
214 k = kh_put(i, h, bin, &ret);
216 if (ret) { // not present
218 l->list = (pair64_t*)calloc(l->m, 16);
222 l->list = (pair64_t*)realloc(l->list, l->m * 16);
224 l->list[l->n].u = beg; l->list[l->n++].v = end;
227 static inline void insert_offset2(ti_lidx_t *index2, int _beg, int _end, uint64_t offset)
230 beg = _beg >> TAD_LIDX_SHIFT;
231 end = (_end - 1) >> TAD_LIDX_SHIFT;
232 if (index2->m < end + 1) {
233 int old_m = index2->m;
235 kroundup32(index2->m);
236 index2->offset = (uint64_t*)realloc(index2->offset, index2->m * 8);
237 memset(index2->offset + old_m, 0, 8 * (index2->m - old_m));
239 for (i = beg + 1; i <= end; ++i)
240 if (index2->offset[i] == 0) index2->offset[i] = offset;
244 static void merge_chunks(ti_index_t *idx)
249 for (i = 0; i < idx->n; ++i) {
250 index = idx->index[i];
251 for (k = kh_begin(index); k != kh_end(index); ++k) {
253 if (!kh_exist(index, k)) continue;
254 p = &kh_value(index, k);
256 for (l = 1; l < p->n; ++l) {
257 if (p->list[m].v>>16 == p->list[l].u>>16) p->list[m].v = p->list[l].v;
258 else p->list[++m] = p->list[l];
265 ti_index_t *ti_index_core(BGZF *fp, const ti_conf_t *conf)
269 uint32_t last_bin, save_bin;
270 int32_t last_coor, last_tid, save_tid;
271 uint64_t save_off, last_off, lineno = 0;
274 str = calloc(1, sizeof(kstring_t));
276 idx = (ti_index_t*)calloc(1, sizeof(ti_index_t));
278 idx->n = idx->max = 0;
279 idx->tname = kh_init(s);
283 save_bin = save_tid = last_tid = last_bin = 0xffffffffu;
284 save_off = last_off = bgzf_tell(fp); last_coor = 0xffffffffu;
285 while ((ret = ti_readline(fp, str)) >= 0) {
288 if (lineno <= idx->conf.line_skip || str->s[0] == idx->conf.meta_char) {
289 last_off = bgzf_tell(fp);
292 get_intv(idx, str, &intv);
293 if (last_tid != intv.tid) { // change of chromosomes
295 last_bin = 0xffffffffu;
296 } else if (last_coor > intv.beg) {
297 fprintf(stderr, "[ti_index_core] the file is not sorted at line %llu\n", (unsigned long long)lineno);
300 if (intv.bin < 4681) insert_offset2(&idx->index2[intv.tid], intv.beg, intv.end, last_off);
301 if (intv.bin != last_bin) { // then possibly write the binning index
302 if (save_bin != 0xffffffffu) // save_bin==0xffffffffu only happens to the first record
303 insert_offset(idx->index[save_tid], save_bin, save_off, last_off);
305 save_bin = last_bin = intv.bin;
307 if (save_tid < 0) break;
309 if (bgzf_tell(fp) <= last_off) {
310 fprintf(stderr, "[ti_index_core] bug in BGZF: %llx < %llx\n",
311 (unsigned long long)bgzf_tell(fp), (unsigned long long)last_off);
314 last_off = bgzf_tell(fp);
315 last_coor = intv.beg;
317 if (save_tid >= 0) insert_offset(idx->index[save_tid], save_bin, save_off, bgzf_tell(fp));
320 free(str->s); free(str);
324 void ti_index_destroy(ti_index_t *idx)
328 if (idx == 0) return;
329 // destroy the name hash table
330 for (k = kh_begin(idx->tname); k != kh_end(idx->tname); ++k) {
331 if (kh_exist(idx->tname, k))
332 free((char*)kh_key(idx->tname, k));
334 kh_destroy(s, idx->tname);
335 // destroy the binning index
336 for (i = 0; i < idx->n; ++i) {
337 khash_t(i) *index = idx->index[i];
338 ti_lidx_t *index2 = idx->index2 + i;
339 for (k = kh_begin(index); k != kh_end(index); ++k) {
340 if (kh_exist(index, k))
341 free(kh_value(index, k).list);
343 kh_destroy(i, index);
344 free(index2->offset);
347 // destroy the linear index
356 void ti_index_save(const ti_index_t *idx, BGZF *fp)
358 int32_t i, size, ti_is_be;
360 ti_is_be = bam_is_big_endian();
361 bgzf_write(fp, "TBI\1", 4);
364 bgzf_write(fp, bam_swap_endian_4p(&x), 4);
365 } else bgzf_write(fp, &idx->n, 4);
366 assert(sizeof(ti_conf_t) == 24);
367 if (ti_is_be) { // write ti_conf_t;
369 memcpy(x, &idx->conf, 24);
370 for (i = 0; i < 6; ++i) bgzf_write(fp, bam_swap_endian_4p(&x[i]), 4);
371 } else bgzf_write(fp, &idx->conf, sizeof(ti_conf_t));
372 { // write target names
375 name = calloc(kh_size(idx->tname), sizeof(void*));
376 for (k = kh_begin(idx->tname); k != kh_end(idx->tname); ++k)
377 if (kh_exist(idx->tname, k))
378 name[kh_value(idx->tname, k)] = (char*)kh_key(idx->tname, k);
379 for (i = 0; i < kh_size(idx->tname); ++i)
380 l += strlen(name[i]) + 1;
381 if (ti_is_be) bgzf_write(fp, bam_swap_endian_4p(&l), 4);
382 else bgzf_write(fp, &l, 4);
383 for (i = 0; i < kh_size(idx->tname); ++i)
384 bgzf_write(fp, name[i], strlen(name[i]) + 1);
387 for (i = 0; i < idx->n; ++i) {
388 khash_t(i) *index = idx->index[i];
389 ti_lidx_t *index2 = idx->index2 + i;
390 // write binning index
391 size = kh_size(index);
392 if (ti_is_be) { // big endian
394 bgzf_write(fp, bam_swap_endian_4p(&x), 4);
395 } else bgzf_write(fp, &size, 4);
396 for (k = kh_begin(index); k != kh_end(index); ++k) {
397 if (kh_exist(index, k)) {
398 ti_binlist_t *p = &kh_value(index, k);
399 if (ti_is_be) { // big endian
401 x = kh_key(index, k); bgzf_write(fp, bam_swap_endian_4p(&x), 4);
402 x = p->n; bgzf_write(fp, bam_swap_endian_4p(&x), 4);
403 for (x = 0; (int)x < p->n; ++x) {
404 bam_swap_endian_8p(&p->list[x].u);
405 bam_swap_endian_8p(&p->list[x].v);
407 bgzf_write(fp, p->list, 16 * p->n);
408 for (x = 0; (int)x < p->n; ++x) {
409 bam_swap_endian_8p(&p->list[x].u);
410 bam_swap_endian_8p(&p->list[x].v);
413 bgzf_write(fp, &kh_key(index, k), 4);
414 bgzf_write(fp, &p->n, 4);
415 bgzf_write(fp, p->list, 16 * p->n);
419 // write linear index (index2)
422 bgzf_write(fp, bam_swap_endian_4p(&x), 4);
423 } else bgzf_write(fp, &index2->n, 4);
424 if (ti_is_be) { // big endian
426 for (x = 0; (int)x < index2->n; ++x)
427 bam_swap_endian_8p(&index2->offset[x]);
428 bgzf_write(fp, index2->offset, 8 * index2->n);
429 for (x = 0; (int)x < index2->n; ++x)
430 bam_swap_endian_8p(&index2->offset[x]);
431 } else bgzf_write(fp, index2->offset, 8 * index2->n);
435 static ti_index_t *ti_index_load_core(BGZF *fp)
440 ti_is_be = bam_is_big_endian();
442 fprintf(stderr, "[ti_index_load_core] fail to load index.\n");
445 bgzf_read(fp, magic, 4);
446 if (strncmp(magic, "TBI\1", 4)) {
447 fprintf(stderr, "[ti_index_load] wrong magic number.\n");
450 idx = (ti_index_t*)calloc(1, sizeof(ti_index_t));
451 bgzf_read(fp, &idx->n, 4);
452 if (ti_is_be) bam_swap_endian_4p(&idx->n);
453 idx->tname = kh_init(s);
454 idx->index = (khash_t(i)**)calloc(idx->n, sizeof(void*));
455 idx->index2 = (ti_lidx_t*)calloc(idx->n, sizeof(ti_lidx_t));
457 bgzf_read(fp, &idx->conf, sizeof(ti_conf_t));
459 bam_swap_endian_4p(&idx->conf.preset);
460 bam_swap_endian_4p(&idx->conf.sc);
461 bam_swap_endian_4p(&idx->conf.bc);
462 bam_swap_endian_4p(&idx->conf.ec);
463 bam_swap_endian_4p(&idx->conf.meta_char);
464 bam_swap_endian_4p(&idx->conf.line_skip);
466 { // read target names
471 bgzf_read(fp, &l, 4);
472 if (ti_is_be) bam_swap_endian_4p(&l);
474 bgzf_read(fp, buf, l);
475 str = calloc(1, sizeof(kstring_t));
476 for (i = j = 0; i < l; ++i) {
478 khint_t k = kh_put(s, idx->tname, strdup(str->s), &ret);
479 kh_value(idx->tname, k) = j++;
481 } else kputc(buf[i], str);
483 free(str->s); free(str); free(buf);
485 for (i = 0; i < idx->n; ++i) {
487 ti_lidx_t *index2 = idx->index2 + i;
492 index = idx->index[i] = kh_init(i);
493 // load binning index
494 bgzf_read(fp, &size, 4);
495 if (ti_is_be) bam_swap_endian_4p(&size);
496 for (j = 0; j < (int)size; ++j) {
497 bgzf_read(fp, &key, 4);
498 if (ti_is_be) bam_swap_endian_4p(&key);
499 k = kh_put(i, index, key, &ret);
500 p = &kh_value(index, k);
501 bgzf_read(fp, &p->n, 4);
502 if (ti_is_be) bam_swap_endian_4p(&p->n);
504 p->list = (pair64_t*)malloc(p->m * 16);
505 bgzf_read(fp, p->list, 16 * p->n);
508 for (x = 0; x < p->n; ++x) {
509 bam_swap_endian_8p(&p->list[x].u);
510 bam_swap_endian_8p(&p->list[x].v);
515 bgzf_read(fp, &index2->n, 4);
516 if (ti_is_be) bam_swap_endian_4p(&index2->n);
517 index2->m = index2->n;
518 index2->offset = (uint64_t*)calloc(index2->m, 8);
519 bgzf_read(fp, index2->offset, index2->n * 8);
521 for (j = 0; j < index2->n; ++j) bam_swap_endian_8p(&index2->offset[j]);
526 ti_index_t *ti_index_load_local(const char *_fn)
531 if (strstr(_fn, "ftp://") == _fn || strstr(_fn, "http://") == _fn) {
534 for (p = _fn + l - 1; p >= _fn; --p)
535 if (*p == '/') break;
537 } else fn = strdup(_fn);
538 fnidx = (char*)calloc(strlen(fn) + 5, 1);
539 strcpy(fnidx, fn); strcat(fnidx, ".tbi");
540 fp = bgzf_open(fnidx, "r");
541 free(fnidx); free(fn);
543 ti_index_t *idx = ti_index_load_core(fp);
550 static void download_from_remote(const char *url)
552 const int buf_size = 1 * 1024 * 1024;
558 if (strstr(url, "ftp://") != url && strstr(url, "http://") != url) return;
560 for (fn = (char*)url + l - 1; fn >= url; --fn)
561 if (*fn == '/') break;
562 ++fn; // fn now points to the file name
563 fp_remote = knet_open(url, "r");
564 if (fp_remote == 0) {
565 fprintf(stderr, "[download_from_remote] fail to open remote file.\n");
568 if ((fp = fopen(fn, "w")) == 0) {
569 fprintf(stderr, "[download_from_remote] fail to create file in the working directory.\n");
570 knet_close(fp_remote);
573 buf = (uint8_t*)calloc(buf_size, 1);
574 while ((l = knet_read(fp_remote, buf, buf_size)) != 0)
575 fwrite(buf, 1, l, fp);
578 knet_close(fp_remote);
581 static void download_from_remote(const char *url)
587 ti_index_t *ti_index_load(const char *fn)
590 idx = ti_index_load_local(fn);
591 if (idx == 0 && (strstr(fn, "ftp://") == fn || strstr(fn, "http://") == fn)) {
592 char *fnidx = calloc(strlen(fn) + 5, 1);
593 strcat(strcpy(fnidx, fn), ".tbi");
594 fprintf(stderr, "[ti_index_load] attempting to download the remote index file.\n");
595 download_from_remote(fnidx);
596 idx = ti_index_load_local(fn);
598 if (idx == 0) fprintf(stderr, "[ti_index_load] fail to load BAM index.\n");
602 int ti_index_build2(const char *fn, const ti_conf_t *conf, const char *_fnidx)
607 if ((fp = bgzf_open(fn, "r")) == 0) {
608 fprintf(stderr, "[ti_index_build2] fail to open the BAM file.\n");
611 idx = ti_index_core(fp, conf);
614 fnidx = (char*)calloc(strlen(fn) + 5, 1);
615 strcpy(fnidx, fn); strcat(fnidx, ".tbi");
616 } else fnidx = strdup(_fnidx);
617 fpidx = bgzf_open(fnidx, "w");
619 fprintf(stderr, "[ti_index_build2] fail to create the index file.\n");
623 ti_index_save(idx, fpidx);
624 ti_index_destroy(idx);
630 int ti_index_build(const char *fn, const ti_conf_t *conf)
632 return ti_index_build2(fn, conf, 0);
635 /********************************************
636 * parse a region in the format chr:beg-end *
637 ********************************************/
639 int ti_parse_region(ti_index_t *idx, const char *str, int *tid, int *begin, int *end)
644 khash_t(s) *h = idx->tname;
646 p = s = (char*)malloc(l+1);
647 /* squeeze out "," */
648 for (i = k = 0; i != l; ++i)
649 if (str[i] != ',' && !isspace(str[i])) s[k++] = str[i];
651 for (i = 0; i != k; ++i) if (s[i] == ':') break;
653 iter = kh_get(s, h, s); /* get the tid */
654 if (iter == kh_end(h)) { // name not found
658 *tid = kh_value(h, iter);
659 if (i == k) { /* dump the whole sequence */
660 *begin = 0; *end = 1<<29; free(s);
663 for (p = s + i + 1; i != k; ++i) if (s[i] == '-') break;
669 if (*begin > 0) --*begin;
671 if (*begin > *end) return -1;
675 /*******************************
676 * retrieve a specified region *
677 *******************************/
679 #define MAX_BIN 37450 // =(8^6-1)/7+1
681 static inline int reg2bins(uint32_t beg, uint32_t end, uint16_t list[MAX_BIN])
686 for (k = 1 + (beg>>26); k <= 1 + (end>>26); ++k) list[i++] = k;
687 for (k = 9 + (beg>>23); k <= 9 + (end>>23); ++k) list[i++] = k;
688 for (k = 73 + (beg>>20); k <= 73 + (end>>20); ++k) list[i++] = k;
689 for (k = 585 + (beg>>17); k <= 585 + (end>>17); ++k) list[i++] = k;
690 for (k = 4681 + (beg>>14); k <= 4681 + (end>>14); ++k) list[i++] = k;
694 static inline int is_overlap(int32_t beg, int32_t end, int rbeg, int rend)
696 return (rend > beg && rbeg < end);
699 // ti_fetch helper function retrieves
700 pair64_t *get_chunk_coordinates(const ti_index_t *idx, int tid, int beg, int end, int* cnt_off)
703 int i, n_bins, n_off;
709 bins = (uint16_t*)calloc(MAX_BIN, 2);
710 n_bins = reg2bins(beg, end, bins);
711 index = idx->index[tid];
712 min_off = (beg>>TAD_LIDX_SHIFT >= idx->index2[tid].n)? 0 : idx->index2[tid].offset[beg>>TAD_LIDX_SHIFT];
713 for (i = n_off = 0; i < n_bins; ++i) {
714 if ((k = kh_get(i, index, bins[i])) != kh_end(index))
715 n_off += kh_value(index, k).n;
718 free(bins); return 0;
720 off = (pair64_t*)calloc(n_off, 16);
721 for (i = n_off = 0; i < n_bins; ++i) {
722 if ((k = kh_get(i, index, bins[i])) != kh_end(index)) {
724 ti_binlist_t *p = &kh_value(index, k);
725 for (j = 0; j < p->n; ++j)
726 if (p->list[j].v > min_off) off[n_off++] = p->list[j];
732 ks_introsort(off, n_off, off);
733 // resolve completely contained adjacent blocks
734 for (i = 1, l = 0; i < n_off; ++i)
735 if (off[l].v < off[i].v)
738 // resolve overlaps between adjacent blocks; this may happen due to the merge in indexing
739 for (i = 1; i < n_off; ++i)
740 if (off[i-1].v >= off[i].u) off[i-1].v = off[i].u;
741 { // merge adjacent blocks
742 for (i = 1, l = 0; i < n_off; ++i) {
743 if (off[l].v>>16 == off[i].u>>16) off[l].v = off[i].v;
744 else off[++l] = off[i];
753 int ti_fetch(BGZF *fp, const ti_index_t *idx, int tid, int beg, int end, void *data, ti_fetch_f func)
756 pair64_t *off = get_chunk_coordinates(idx, tid, beg, end, &n_off);
757 if (off == 0) return 0;
759 // retrive alignments
762 kstring_t *str = calloc(1, sizeof(kstring_t));
763 n_seeks = 0; i = -1; curr_off = 0;
765 if (curr_off == 0 || curr_off >= off[i].v) { // then jump to the next chunk
766 if (i == n_off - 1) break; // no more chunks
767 if (i >= 0) assert(curr_off == off[i].v); // otherwise bug
768 if (i < 0 || off[i].v != off[i+1].u) { // not adjacent chunks; then seek
769 bgzf_seek(fp, off[i+1].u, SEEK_SET);
770 curr_off = bgzf_tell(fp);
775 if ((ret = ti_readline(fp, str)) >= 0) {
777 curr_off = bgzf_tell(fp);
778 if (str->s[0] == idx->conf.meta_char) continue;
779 get_intv((ti_index_t*)idx, str, &intv);
780 if (intv.tid != tid || intv.beg >= end) break; // no need to proceed
781 else if (is_overlap(beg, end, intv.beg, intv.end)) func(str->l, str->s, data);
782 } else break; // end of file
784 // fprintf(stderr, "[ti_fetch] # seek calls: %d\n", n_seeks);
785 free(str->s); free(str);