bam_index.c

   1 #include <ctype.h>
   2 #include <assert.h>
   3 #include "bam.h"
   4 #include "khash.h"
   5 #include "ksort.h"
   6 #include "bam_endian.h"
   7 #ifdef _USE_KNETFILE
   8 #include "knetfile.h"
   9 #endif
  10
  11 /*!
  12   @header
  13
  14   Alignment indexing. Before indexing, BAM must be sorted based on the
  15   leftmost coordinate of alignments. In indexing, BAM uses two indices:
  16   a UCSC binning index and a simple linear index. The binning index is
  17   efficient for alignments spanning long distance, while the auxiliary
  18   linear index helps to reduce unnecessary seek calls especially for
  19   short alignments.
  20
  21   The UCSC binning scheme was suggested by Richard Durbin and Lincoln
  22   Stein and is explained by Kent et al. (2002). In this scheme, each bin
  23   represents a contiguous genomic region which can be fully contained in
  24   another bin; each alignment is associated with a bin which represents
  25   the smallest region containing the entire alignment. The binning
  26   scheme is essentially another representation of R-tree. A distinct bin
  27   uniquely corresponds to a distinct internal node in a R-tree. Bin A is
  28   a child of Bin B if region A is contained in B.
  29
  30   In BAM, each bin may span 2^29, 2^26, 2^23, 2^20, 2^17 or 2^14 bp. Bin
  31   0 spans a 512Mbp region, bins 1-8 span 64Mbp, 9-72 8Mbp, 73-584 1Mbp,
  32   585-4680 128Kbp and bins 4681-37449 span 16Kbp regions. If we want to
  33   find the alignments overlapped with a region [rbeg,rend), we need to
  34   calculate the list of bins that may be overlapped the region and test
  35   the alignments in the bins to confirm the overlaps. If the specified
  36   region is short, typically only a few alignments in six bins need to
  37   be retrieved. The overlapping alignments can be quickly fetched.
  38
  39  */
  40
  41 #define BAM_MIN_CHUNK_GAP 32768
  42 // 1<<14 is the size of minimum bin.
  43 #define BAM_LIDX_SHIFT    14
  44
  45 #define BAM_MAX_BIN 37450 // =(8^6-1)/7+1
  46
  47 typedef struct {
  48         uint64_t u, v;
  49 } pair64_t;
  50
  51 #define pair64_lt(a,b) ((a).u < (b).u)
  52 KSORT_INIT(off, pair64_t, pair64_lt)
  53
  54 typedef struct {
  55         uint32_t m, n;
  56         pair64_t *list;
  57 } bam_binlist_t;
  58
  59 typedef struct {
  60         int32_t n, m;
  61         uint64_t *offset;
  62 } bam_lidx_t;
  63
  64 KHASH_MAP_INIT_INT(i, bam_binlist_t)
  65
  66 struct __bam_index_t {
  67         int32_t n;
  68         uint64_t n_no_coor; // unmapped reads without coordinate
  69         khash_t(i) **index;
  70         bam_lidx_t *index2;
  71 };
  72
  73 // requirement: len <= LEN_MASK
  74 static inline void insert_offset(khash_t(i) *h, int bin, uint64_t beg, uint64_t end)
  75 {
  76         khint_t k;
  77         bam_binlist_t *l;
  78         int ret;
  79         k = kh_put(i, h, bin, &ret);
  80         l = &kh_value(h, k);
  81         if (ret) { // not present
  82                 l->m = 1; l->n = 0;
  83                 l->list = (pair64_t*)calloc(l->m, 16);
  84         }
  85         if (l->n == l->m) {
  86                 l->m <<= 1;
  87                 l->list = (pair64_t*)realloc(l->list, l->m * 16);
  88         }
  89         l->list[l->n].u = beg; l->list[l->n++].v = end;
  90 }
  91
  92 static inline void insert_offset2(bam_lidx_t *index2, bam1_t *b, uint64_t offset)
  93 {
  94         int i, beg, end;
  95         beg = b->core.pos >> BAM_LIDX_SHIFT;
  96         end = (bam_calend(&b->core, bam1_cigar(b)) - 1) >> BAM_LIDX_SHIFT;
  97         if (index2->m < end + 1) {
  98                 int old_m = index2->m;
  99                 index2->m = end + 1;
 100                 kroundup32(index2->m);
 101                 index2->offset = (uint64_t*)realloc(index2->offset, index2->m * 8);
 102                 memset(index2->offset + old_m, 0, 8 * (index2->m - old_m));
 103         }
 104         if (beg == end) {
 105                 if (index2->offset[beg] == 0) index2->offset[beg] = offset;
 106         } else {
 107                 for (i = beg; i <= end; ++i)
 108                         if (index2->offset[i] == 0) index2->offset[i] = offset;
 109         }
 110         index2->n = end + 1;
 111 }
 112
 113 static void merge_chunks(bam_index_t *idx)
 114 {
 115 #if defined(BAM_TRUE_OFFSET) || defined(BAM_VIRTUAL_OFFSET16)
 116         khash_t(i) *index;
 117         int i, l, m;
 118         khint_t k;
 119         for (i = 0; i < idx->n; ++i) {
 120                 index = idx->index[i];
 121                 for (k = kh_begin(index); k != kh_end(index); ++k) {
 122                         bam_binlist_t *p;
 123                         if (!kh_exist(index, k) || kh_key(index, k) == BAM_MAX_BIN) continue;
 124                         p = &kh_value(index, k);
 125                         m = 0;
 126                         for (l = 1; l < p->n; ++l) {
 127 #ifdef BAM_TRUE_OFFSET
 128                                 if (p->list[m].v + BAM_MIN_CHUNK_GAP > p->list[l].u) p->list[m].v = p->list[l].v;
 129 #else
 130                                 if (p->list[m].v>>16 == p->list[l].u>>16) p->list[m].v = p->list[l].v;
 131 #endif
 132                                 else p->list[++m] = p->list[l];
 133                         } // ~for(l)
 134                         p->n = m + 1;
 135                 } // ~for(k)
 136         } // ~for(i)
 137 #endif // defined(BAM_TRUE_OFFSET) || defined(BAM_BGZF)
 138 }
 139
 140 static void fill_missing(bam_index_t *idx)
 141 {
 142         int i, j;
 143         for (i = 0; i < idx->n; ++i) {
 144                 bam_lidx_t *idx2 = &idx->index2[i];
 145                 for (j = 1; j < idx2->n; ++j)
 146                         if (idx2->offset[j] == 0)
 147                                 idx2->offset[j] = idx2->offset[j-1];
 148         }
 149 }
 150
 151 bam_index_t *bam_index_core(bamFile fp)
 152 {
 153         bam1_t *b;
 154         bam_header_t *h;
 155         int i, ret;
 156         bam_index_t *idx;
 157         uint32_t last_bin, save_bin;
 158         int32_t last_coor, last_tid, save_tid;
 159         bam1_core_t *c;
 160         uint64_t save_off, last_off, n_mapped, n_unmapped, off_beg, off_end, n_no_coor;
 161
 162         idx = (bam_index_t*)calloc(1, sizeof(bam_index_t));
 163         b = (bam1_t*)calloc(1, sizeof(bam1_t));
 164         h = bam_header_read(fp);
 165         c = &b->core;
 166
 167         idx->n = h->n_targets;
 168         bam_header_destroy(h);
 169         idx->index = (khash_t(i)**)calloc(idx->n, sizeof(void*));
 170         for (i = 0; i < idx->n; ++i) idx->index[i] = kh_init(i);
 171         idx->index2 = (bam_lidx_t*)calloc(idx->n, sizeof(bam_lidx_t));
 172
 173         save_bin = save_tid = last_tid = last_bin = 0xffffffffu;
 174         save_off = last_off = bam_tell(fp); last_coor = 0xffffffffu;
 175         n_mapped = n_unmapped = n_no_coor = off_end = 0;
 176         off_beg = off_end = bam_tell(fp);
 177         while ((ret = bam_read1(fp, b)) >= 0) {
 178                 if (c->tid < 0) ++n_no_coor;
 179                 if (last_tid < c->tid || (last_tid >= 0 && c->tid < 0)) { // change of chromosomes
 180                         last_tid = c->tid;
 181                         last_bin = 0xffffffffu;
 182                 } else if ((uint32_t)last_tid > (uint32_t)c->tid) {
 183                         fprintf(stderr, "[bam_index_core] the alignment is not sorted (%s): %d-th chr > %d-th chr\n",
 184                                         bam1_qname(b), last_tid+1, c->tid+1);
 185                         return NULL;
 186                 } else if ((int32_t)c->tid >= 0 && last_coor > c->pos) {
 187                         fprintf(stderr, "[bam_index_core] the alignment is not sorted (%s): %u > %u in %d-th chr\n",
 188                                         bam1_qname(b), last_coor, c->pos, c->tid+1);
 189                         return NULL;
 190                 }
 191                 if (c->tid >= 0 && !(c->flag & BAM_FUNMAP)) insert_offset2(&idx->index2[b->core.tid], b, last_off);
 192                 if (c->bin != last_bin) { // then possibly write the binning index
 193                         if (save_bin != 0xffffffffu) // save_bin==0xffffffffu only happens to the first record
 194                                 insert_offset(idx->index[save_tid], save_bin, save_off, last_off);
 195                         if (last_bin == 0xffffffffu && save_tid != 0xffffffffu) { // write the meta element
 196                                 off_end = last_off;
 197                                 insert_offset(idx->index[save_tid], BAM_MAX_BIN, off_beg, off_end);
 198                                 insert_offset(idx->index[save_tid], BAM_MAX_BIN, n_mapped, n_unmapped);
 199                                 n_mapped = n_unmapped = 0;
 200                                 off_beg = off_end;
 201                         }
 202                         save_off = last_off;
 203                         save_bin = last_bin = c->bin;
 204                         save_tid = c->tid;
 205                         if (save_tid < 0) break;
 206                 }
 207                 if (bam_tell(fp) <= last_off) {
 208                         fprintf(stderr, "[bam_index_core] bug in BGZF/RAZF: %llx < %llx\n",
 209                                         (unsigned long long)bam_tell(fp), (unsigned long long)last_off);
 210                         return NULL;
 211                 }
 212                 if (c->flag & BAM_FUNMAP) ++n_unmapped;
 213                 else ++n_mapped;
 214                 last_off = bam_tell(fp);
 215                 last_coor = b->core.pos;
 216         }
 217         if (save_tid >= 0) {
 218                 insert_offset(idx->index[save_tid], save_bin, save_off, bam_tell(fp));
 219                 insert_offset(idx->index[save_tid], BAM_MAX_BIN, off_beg, bam_tell(fp));
 220                 insert_offset(idx->index[save_tid], BAM_MAX_BIN, n_mapped, n_unmapped);
 221         }
 222         merge_chunks(idx);
 223         fill_missing(idx);
 224         if (ret >= 0) {
 225                 while ((ret = bam_read1(fp, b)) >= 0) {
 226                         ++n_no_coor;
 227                         if (c->tid >= 0 && n_no_coor) {
 228                                 fprintf(stderr, "[bam_index_core] the alignment is not sorted: reads without coordinates prior to reads with coordinates.\n");
 229                                 return NULL;
 230                         }
 231                 }
 232         }
 233         if (ret < -1) fprintf(stderr, "[bam_index_core] truncated file? Continue anyway. (%d)\n", ret);
 234         free(b->data); free(b);
 235         idx->n_no_coor = n_no_coor;
 236         return idx;
 237 }
 238
 239 void bam_index_destroy(bam_index_t *idx)
 240 {
 241         khint_t k;
 242         int i;
 243         if (idx == 0) return;
 244         for (i = 0; i < idx->n; ++i) {
 245                 khash_t(i) *index = idx->index[i];
 246                 bam_lidx_t *index2 = idx->index2 + i;
 247                 for (k = kh_begin(index); k != kh_end(index); ++k) {
 248                         if (kh_exist(index, k))
 249                                 free(kh_value(index, k).list);
 250                 }
 251                 kh_destroy(i, index);
 252                 free(index2->offset);
 253         }
 254         free(idx->index); free(idx->index2);
 255         free(idx);
 256 }
 257
 258 void bam_index_save(const bam_index_t *idx, FILE *fp)
 259 {
 260         int32_t i, size;
 261         khint_t k;
 262         if (fwrite("BAI\1", 1, 4, fp) < 4) {
 263                 fprintf(stderr, "[%s] failed to write magic number.\n", __func__);
 264                 return;
 265         }
 266         if (bam_is_be) {
 267                 uint32_t x = idx->n;
 268                 if (fwrite(bam_swap_endian_4p(&x), 4, 1, fp) < 1) {
 269                         fprintf(stderr, "[%s] failed to write n_ref.\n", __func__);
 270                         return;
 271                 }
 272         }
 273         else {
 274                 if (fwrite(&idx->n, 4, 1, fp) < 1) {
 275                         fprintf(stderr, "[%s] failed to write n_ref.\n", __func__);
 276                         return;
 277                 }
 278         }
 279         for (i = 0; i < idx->n; ++i) {
 280                 khash_t(i) *index = idx->index[i];
 281                 bam_lidx_t *index2 = idx->index2 + i;
 282                 // write binning index
 283                 size = kh_size(index);
 284                 if (bam_is_be) { // big endian
 285                         uint32_t x = size;
 286                         if (fwrite(bam_swap_endian_4p(&x), 4, 1, fp) < 1) {
 287                                 fprintf(stderr, "[%s] failed to write n_bin.\n", __func__);
 288                                 return;
 289                         }
 290                 }
 291                 else {
 292                         if (fwrite(&size, 4, 1, fp) < 1) {
 293                                 fprintf(stderr, "[%s] failed to write n_bin.\n", __func__);
 294                                 return;
 295                         }
 296                 }
 297                 for (k = kh_begin(index); k != kh_end(index); ++k) {
 298                         if (kh_exist(index, k)) {
 299                                 bam_binlist_t *p = &kh_value(index, k);
 300                                 if (bam_is_be) { // big endian
 301                                         uint32_t x;
 302                                         x = kh_key(index, k);
 303                                         if (fwrite(bam_swap_endian_4p(&x), 4, 1, fp) < 1) {
 304                                                 fprintf(stderr, "[%s] failed to write bin.\n", __func__);
 305                                                 return;
 306                                         }
 307                                         x = p->n;
 308                                         if (fwrite(bam_swap_endian_4p(&x), 4, 1, fp) < 1) {
 309                                                 fprintf(stderr, "[%s] failed to write n_chunk.\n", __func__);
 310                                                 return;
 311                                         }
 312                                         for (x = 0; (int)x < p->n; ++x) {
 313                                                 bam_swap_endian_8p(&p->list[x].u);
 314                                                 bam_swap_endian_8p(&p->list[x].v);
 315                                         }
 316                                         if (fwrite(p->list, 16, p->n, fp) < p->n) {
 317                                                 fprintf(stderr, "[%s] failed to write %d chunks.\n", __func__, p->n);
 318                                                 return;
 319                                         }
 320                                         /*
 321                                         for (x = 0; (int)x < p->n; ++x) {
 322                                                 bam_swap_endian_8p(&p->list[x].u);
 323                                                 bam_swap_endian_8p(&p->list[x].v);
 324                                         } */
 325                                 }
 326                                 else {
 327                                         if (fwrite(&kh_key(index, k), 4, 1, fp) < 1) {
 328                                                 fprintf(stderr, "[%s] failed to write bin.\n", __func__);
 329                                                 return;
 330                                         }
 331                                         if (fwrite(&p->n, 4, 1, fp) < 1) {
 332                                                 fprintf(stderr, "[%s] failed to write n_chunk.\n", __func__);
 333                                                 return;
 334                                         }
 335                                         if (fwrite(p->list, 16, p->n, fp) < p->n) {
 336                                                 fprintf(stderr, "[%s] failed to write %d chunks.\n", __func__, p->n);
 337                                                 return;
 338                                         }
 339                                 }
 340                         }
 341                 }
 342                 // write linear index (index2)
 343                 if (bam_is_be) {
 344                         int x = index2->n;
 345                         if (fwrite(bam_swap_endian_4p(&x), 4, 1, fp) < 1) {
 346                                 fprintf(stderr, "[%s] failed to write n_intv.\n", __func__);
 347                                 return;
 348                         }
 349                 }
 350                 else {
 351                         if (fwrite(&index2->n, 4, 1, fp) < 1) {
 352                                 fprintf(stderr, "[%s] failed to write n_intv.\n", __func__);
 353                                 return;
 354                         }
 355                 }
 356                 if (bam_is_be) { // big endian
 357                         int x;
 358                         for (x = 0; (int)x < index2->n; ++x)
 359                                 bam_swap_endian_8p(&index2->offset[x]);
 360                         if (fwrite(index2->offset, 8, index2->n, fp) < index2->n) {
 361                                 fprintf(stderr, "[%s] failed to write ioffset.\n", __func__);
 362                                 return;
 363                         }
 364                         for (x = 0; (int)x < index2->n; ++x)
 365                                 bam_swap_endian_8p(&index2->offset[x]);
 366                 }
 367                 else {
 368                         if (fwrite(index2->offset, 8, index2->n, fp) < index2->n) {
 369                                 fprintf(stderr, "[%s] failed to write ioffset.\n", __func__);
 370                                 return;
 371                         }
 372                 }
 373         }
 374         { // write the number of reads coor-less records.
 375                 uint64_t x = idx->n_no_coor;
 376                 if (bam_is_be) bam_swap_endian_8p(&x);
 377                 if (fwrite(&x, 8, 1, fp) < 1) {
 378                         fprintf(stderr, "[%s] failed to write n_no_coor.\n", __func__);
 379                 }
 380         }
 381         fflush(fp);
 382 }
 383
 384 static bam_index_t *bam_index_load_core(FILE *fp)
 385 {
 386         int i;
 387         char magic[4];
 388         bam_index_t *idx;
 389         if (fp == 0) {
 390                 fprintf(stderr, "[bam_index_load_core] fail to load index.\n");
 391                 return 0;
 392         }
 393         if (fread(magic, 1, 4, fp) < 4) {
 394                 fprintf(stderr, "[%s] failed to read magic number.\n", __func__);
 395                 fclose(fp);
 396                 return 0;
 397         }
 398         if (strncmp(magic, "BAI\1", 4)) {
 399                 fprintf(stderr, "[bam_index_load] wrong magic number.\n");
 400                 fclose(fp);
 401                 return 0;
 402         }
 403         idx = (bam_index_t*)calloc(1, sizeof(bam_index_t));
 404         if (fread(&idx->n, 4, 1, fp) < 1) {
 405                 fprintf(stderr, "[%s] failed to read n_ref.\n", __func__);
 406                 fclose(fp);
 407                 return 0;
 408         }
 409         if (bam_is_be) bam_swap_endian_4p(&idx->n);
 410         idx->index = (khash_t(i)**)calloc(idx->n, sizeof(void*));
 411         idx->index2 = (bam_lidx_t*)calloc(idx->n, sizeof(bam_lidx_t));
 412         for (i = 0; i < idx->n; ++i) {
 413                 khash_t(i) *index;
 414                 bam_lidx_t *index2 = idx->index2 + i;
 415                 uint32_t key, size;
 416                 khint_t k;
 417                 int j, ret;
 418                 bam_binlist_t *p;
 419                 index = idx->index[i] = kh_init(i);
 420                 // load binning index
 421                 if (fread(&size, 4, 1, fp) < 1) {
 422                         fprintf(stderr, "[%s] failed to read n_bin.\n", __func__);
 423                         fclose(fp);
 424                         return 0;
 425                 }
 426                 if (bam_is_be) bam_swap_endian_4p(&size);
 427                 for (j = 0; j < (int)size; ++j) {
 428                         if (fread(&key, 4, 1, fp) < 1) {
 429                                 fprintf(stderr, "[%s] failed to read bin.\n", __func__);
 430                                 fclose(fp);
 431                                 return 0;
 432                         }
 433                         if (bam_is_be) bam_swap_endian_4p(&key);
 434                         k = kh_put(i, index, key, &ret);
 435                         p = &kh_value(index, k);
 436                         if (fread(&p->n, 4, 1, fp) < 1) {
 437                                 fprintf(stderr, "[%s] failed to read n_chunk.\n", __func__);
 438                                 fclose(fp);
 439                                 return 0;
 440                         }
 441                         if (bam_is_be) bam_swap_endian_4p(&p->n);
 442                         p->m = p->n;
 443                         p->list = (pair64_t*)malloc(p->m * 16);
 444                         if (fread(p->list, 16, p->n, fp) < p->n) {
 445                                 fprintf(stderr, "[%s] failed to read %d chunks.\n", __func__, p->n);
 446                                 fclose(fp);
 447                                 return 0;
 448                         }
 449                         if (bam_is_be) {
 450                                 int x;
 451                                 for (x = 0; x < p->n; ++x) {
 452                                         bam_swap_endian_8p(&p->list[x].u);
 453                                         bam_swap_endian_8p(&p->list[x].v);
 454                                 }
 455                         }
 456                 }
 457                 // load linear index
 458                 if (fread(&index2->n, 4, 1, fp) < 1) {
 459                         fprintf(stderr, "[%s] failed to read n_intv", __func__);
 460                         fclose(fp);
 461                         return 0;
 462                 }
 463                 if (bam_is_be) bam_swap_endian_4p(&index2->n);
 464                 if (index2->n > 0) {
 465                         index2->m = index2->n;
 466                         index2->offset = (uint64_t*)calloc(index2->m, 8);
 467                         if (fread(index2->offset, index2->n, 8, fp) < index2->n) {
 468                                 fprintf(stderr, "[%s] failed to read %d intervals.", __func__, index2->n);
 469                                 fclose(fp);
 470                                 return 0;
 471                         }
 472                         if (bam_is_be) {
 473                                 for (j = 0; j < index2->n; ++j) {
 474                                         bam_swap_endian_8p(&index2->offset[j]);
 475                                 }
 476                         }
 477                 }
 478         }
 479         if (fread(&idx->n_no_coor, 8, 1, fp) == 0) idx->n_no_coor = 0;
 480         if (bam_is_be) bam_swap_endian_8p(&idx->n_no_coor);
 481         return idx;
 482 }
 483
 484 bam_index_t *bam_index_load_local(const char *_fn)
 485 {
 486         FILE *fp;
 487         char *fnidx, *fn;
 488
 489         if (strstr(_fn, "ftp://") == _fn || strstr(_fn, "http://") == _fn) {
 490                 const char *p;
 491                 int l = strlen(_fn);
 492                 for (p = _fn + l - 1; p >= _fn; --p)
 493                         if (*p == '/') break;
 494                 fn = strdup(p + 1);
 495         } else fn = strdup(_fn);
 496         fnidx = (char*)calloc(strlen(fn) + 5, 1);
 497         strcpy(fnidx, fn); strcat(fnidx, ".bai");
 498         fp = fopen(fnidx, "rb");
 499         if (fp == 0) { // try "{base}.bai"
 500                 char *s = strstr(fn, "bam");
 501                 if (s == fn + strlen(fn) - 3) {
 502                         strcpy(fnidx, fn);
 503                         fnidx[strlen(fn)-1] = 'i';
 504                         fp = fopen(fnidx, "rb");
 505                 }
 506         }
 507         free(fnidx); free(fn);
 508         if (fp) {
 509                 bam_index_t *idx = bam_index_load_core(fp);
 510                 fclose(fp);
 511                 return idx;
 512         } else return 0;
 513 }
 514
 515 #ifdef _USE_KNETFILE
 516 static void download_from_remote(const char *url)
 517 {
 518         const int buf_size = 1 * 1024 * 1024;
 519         char *fn;
 520         FILE *fp;
 521         uint8_t *buf;
 522         knetFile *fp_remote;
 523         int l;
 524         if (strstr(url, "ftp://") != url && strstr(url, "http://") != url) return;
 525         l = strlen(url);
 526         for (fn = (char*)url + l - 1; fn >= url; --fn)
 527                 if (*fn == '/') break;
 528         ++fn; // fn now points to the file name
 529         fp_remote = knet_open(url, "r");
 530         if (fp_remote == 0) {
 531                 fprintf(stderr, "[download_from_remote] fail to open remote file.\n");
 532                 return;
 533         }
 534         if ((fp = fopen(fn, "wb")) == 0) {
 535                 fprintf(stderr, "[download_from_remote] fail to create file in the working directory.\n");
 536                 knet_close(fp_remote);
 537                 return;
 538         }
 539         buf = (uint8_t*)calloc(buf_size, 1);
 540         while ((l = knet_read(fp_remote, buf, buf_size)) != 0) {
 541                 if (fwrite(buf, 1, l, fp) < l) {
 542                         fprintf(stderr, "[%s] fail to write to destination file.\n", __func__);
 543                         break;
 544                 }
 545         }
 546         free(buf);
 547         fclose(fp);
 548         knet_close(fp_remote);
 549 }
 550 #else
 551 static void download_from_remote(const char *url)
 552 {
 553         return;
 554 }
 555 #endif
 556
 557 bam_index_t *bam_index_load(const char *fn)
 558 {
 559         bam_index_t *idx;
 560         idx = bam_index_load_local(fn);
 561         if (idx == 0 && (strstr(fn, "ftp://") == fn || strstr(fn, "http://") == fn)) {
 562                 char *fnidx = calloc(strlen(fn) + 5, 1);
 563                 strcat(strcpy(fnidx, fn), ".bai");
 564                 fprintf(stderr, "[bam_index_load] attempting to download the remote index file.\n");
 565                 download_from_remote(fnidx);
 566                 idx = bam_index_load_local(fn);
 567         }
 568         if (idx == 0) fprintf(stderr, "[bam_index_load] fail to load BAM index.\n");
 569         return idx;
 570 }
 571
 572 int bam_index_build2(const char *fn, const char *_fnidx)
 573 {
 574         char *fnidx;
 575         FILE *fpidx;
 576         bamFile fp;
 577         bam_index_t *idx;
 578         if ((fp = bam_open(fn, "r")) == 0) {
 579                 fprintf(stderr, "[bam_index_build2] fail to open the BAM file.\n");
 580                 return -1;
 581         }
 582         idx = bam_index_core(fp);
 583         bam_close(fp);
 584         if(idx == 0) {
 585                 fprintf(stderr, "[bam_index_build2] fail to index the BAM file.\n");
 586                 return -1;
 587         }
 588         if (_fnidx == 0) {
 589                 fnidx = (char*)calloc(strlen(fn) + 5, 1);
 590                 strcpy(fnidx, fn); strcat(fnidx, ".bai");
 591         } else fnidx = strdup(_fnidx);
 592         fpidx = fopen(fnidx, "wb");
 593         if (fpidx == 0) {
 594                 fprintf(stderr, "[bam_index_build2] fail to create the index file.\n");
 595                 free(fnidx);
 596                 return -1;
 597         }
 598         bam_index_save(idx, fpidx);
 599         bam_index_destroy(idx);
 600         fclose(fpidx);
 601         free(fnidx);
 602         return 0;
 603 }
 604
 605 int bam_index_build(const char *fn)
 606 {
 607         return bam_index_build2(fn, 0);
 608 }
 609
 610 int bam_index(int argc, char *argv[])
 611 {
 612         if (argc < 2) {
 613                 fprintf(stderr, "Usage: samtools index <in.bam> [out.index]\n");
 614                 return 1;
 615         }
 616         if (argc >= 3) bam_index_build2(argv[1], argv[2]);
 617         else bam_index_build(argv[1]);
 618         return 0;
 619 }
 620
 621 int bam_idxstats(int argc, char *argv[])
 622 {
 623         bam_index_t *idx;
 624         bam_header_t *header;
 625         bamFile fp;
 626         int i;
 627         if (argc < 2) {
 628                 fprintf(stderr, "Usage: samtools idxstats <in.bam>\n");
 629                 return 1;
 630         }
 631         fp = bam_open(argv[1], "r");
 632         if (fp == 0) { fprintf(stderr, "[%s] fail to open BAM.\n", __func__); return 1; }
 633         header = bam_header_read(fp);
 634         bam_close(fp);
 635         idx = bam_index_load(argv[1]);
 636         if (idx == 0) { fprintf(stderr, "[%s] fail to load the index.\n", __func__); return 1; }
 637         for (i = 0; i < idx->n; ++i) {
 638                 khint_t k;
 639                 khash_t(i) *h = idx->index[i];
 640                 printf("%s\t%d", header->target_name[i], header->target_len[i]);
 641                 k = kh_get(i, h, BAM_MAX_BIN);
 642                 if (k != kh_end(h))
 643                         printf("\t%llu\t%llu", (long long)kh_val(h, k).list[1].u, (long long)kh_val(h, k).list[1].v);
 644                 else printf("\t0\t0");
 645                 putchar('\n');
 646         }
 647         printf("*\t0\t0\t%llu\n", (long long)idx->n_no_coor);
 648         bam_header_destroy(header);
 649         bam_index_destroy(idx);
 650         return 0;
 651 }
 652
 653 static inline int reg2bins(uint32_t beg, uint32_t end, uint16_t list[BAM_MAX_BIN])
 654 {
 655         int i = 0, k;
 656         if (beg >= end) return 0;
 657         if (end >= 1u<<29) end = 1u<<29;
 658         --end;
 659         list[i++] = 0;
 660         for (k =    1 + (beg>>26); k <=    1 + (end>>26); ++k) list[i++] = k;
 661         for (k =    9 + (beg>>23); k <=    9 + (end>>23); ++k) list[i++] = k;
 662         for (k =   73 + (beg>>20); k <=   73 + (end>>20); ++k) list[i++] = k;
 663         for (k =  585 + (beg>>17); k <=  585 + (end>>17); ++k) list[i++] = k;
 664         for (k = 4681 + (beg>>14); k <= 4681 + (end>>14); ++k) list[i++] = k;
 665         return i;
 666 }
 667
 668 static inline int is_overlap(uint32_t beg, uint32_t end, const bam1_t *b)
 669 {
 670         uint32_t rbeg = b->core.pos;
 671         uint32_t rend = b->core.n_cigar? bam_calend(&b->core, bam1_cigar(b)) : b->core.pos + 1;
 672         return (rend > beg && rbeg < end);
 673 }
 674
 675 struct __bam_iter_t {
 676         int from_first; // read from the first record; no random access
 677         int tid, beg, end, n_off, i, finished;
 678         uint64_t curr_off;
 679         pair64_t *off;
 680 };
 681
 682 // bam_fetch helper function retrieves
 683 bam_iter_t bam_iter_query(const bam_index_t *idx, int tid, int beg, int end)
 684 {
 685         uint16_t *bins;
 686         int i, n_bins, n_off;
 687         pair64_t *off;
 688         khint_t k;
 689         khash_t(i) *index;
 690         uint64_t min_off;
 691         bam_iter_t iter = 0;
 692
 693         if (beg < 0) beg = 0;
 694         if (end < beg) return 0;
 695         // initialize iter
 696         iter = calloc(1, sizeof(struct __bam_iter_t));
 697         iter->tid = tid, iter->beg = beg, iter->end = end; iter->i = -1;
 698         //
 699         bins = (uint16_t*)calloc(BAM_MAX_BIN, 2);
 700         n_bins = reg2bins(beg, end, bins);
 701         index = idx->index[tid];
 702         if (idx->index2[tid].n > 0) {
 703                 min_off = (beg>>BAM_LIDX_SHIFT >= idx->index2[tid].n)? idx->index2[tid].offset[idx->index2[tid].n-1]
 704                         : idx->index2[tid].offset[beg>>BAM_LIDX_SHIFT];
 705                 if (min_off == 0) { // improvement for index files built by tabix prior to 0.1.4
 706                         int n = beg>>BAM_LIDX_SHIFT;
 707                         if (n > idx->index2[tid].n) n = idx->index2[tid].n;
 708                         for (i = n - 1; i >= 0; --i)
 709                                 if (idx->index2[tid].offset[i] != 0) break;
 710                         if (i >= 0) min_off = idx->index2[tid].offset[i];
 711                 }
 712         } else min_off = 0; // tabix 0.1.2 may produce such index files
 713         for (i = n_off = 0; i < n_bins; ++i) {
 714                 if ((k = kh_get(i, index, bins[i])) != kh_end(index))
 715                         n_off += kh_value(index, k).n;
 716         }
 717         if (n_off == 0) {
 718                 free(bins); return iter;
 719         }
 720         off = (pair64_t*)calloc(n_off, 16);
 721         for (i = n_off = 0; i < n_bins; ++i) {
 722                 if ((k = kh_get(i, index, bins[i])) != kh_end(index)) {
 723                         int j;
 724                         bam_binlist_t *p = &kh_value(index, k);
 725                         for (j = 0; j < p->n; ++j)
 726                                 if (p->list[j].v > min_off) off[n_off++] = p->list[j];
 727                 }
 728         }
 729         free(bins);
 730         if (n_off == 0) {
 731                 free(off); return iter;
 732         }
 733         {
 734                 bam1_t *b = (bam1_t*)calloc(1, sizeof(bam1_t));
 735                 int l;
 736                 ks_introsort(off, n_off, off);
 737                 // resolve completely contained adjacent blocks
 738                 for (i = 1, l = 0; i < n_off; ++i)
 739                         if (off[l].v < off[i].v)
 740                                 off[++l] = off[i];
 741                 n_off = l + 1;
 742                 // resolve overlaps between adjacent blocks; this may happen due to the merge in indexing
 743                 for (i = 1; i < n_off; ++i)
 744                         if (off[i-1].v >= off[i].u) off[i-1].v = off[i].u;
 745                 { // merge adjacent blocks
 746 #if defined(BAM_TRUE_OFFSET) || defined(BAM_VIRTUAL_OFFSET16)
 747                         for (i = 1, l = 0; i < n_off; ++i) {
 748 #ifdef BAM_TRUE_OFFSET
 749                                 if (off[l].v + BAM_MIN_CHUNK_GAP > off[i].u) off[l].v = off[i].v;
 750 #else
 751                                 if (off[l].v>>16 == off[i].u>>16) off[l].v = off[i].v;
 752 #endif
 753                                 else off[++l] = off[i];
 754                         }
 755                         n_off = l + 1;
 756 #endif
 757                 }
 758                 bam_destroy1(b);
 759         }
 760         iter->n_off = n_off; iter->off = off;
 761         return iter;
 762 }
 763
 764 pair64_t *get_chunk_coordinates(const bam_index_t *idx, int tid, int beg, int end, int *cnt_off)
 765 { // for pysam compatibility
 766         bam_iter_t iter;
 767         pair64_t *off;
 768         iter = bam_iter_query(idx, tid, beg, end);
 769         off = iter->off; *cnt_off = iter->n_off;
 770         free(iter);
 771         return off;
 772 }
 773
 774 void bam_iter_destroy(bam_iter_t iter)
 775 {
 776         if (iter) { free(iter->off); free(iter); }
 777 }
 778
 779 int bam_iter_read(bamFile fp, bam_iter_t iter, bam1_t *b)
 780 {
 781         int ret;
 782         if (iter && iter->finished) return -1;
 783         if (iter == 0 || iter->from_first) {
 784                 ret = bam_read1(fp, b);
 785                 if (ret < 0 && iter) iter->finished = 1;
 786                 return ret;
 787         }
 788         if (iter->off == 0) return -1;
 789         for (;;) {
 790                 if (iter->curr_off == 0 || iter->curr_off >= iter->off[iter->i].v) { // then jump to the next chunk
 791                         if (iter->i == iter->n_off - 1) { ret = -1; break; } // no more chunks
 792                         if (iter->i >= 0) assert(iter->curr_off == iter->off[iter->i].v); // otherwise bug
 793                         if (iter->i < 0 || iter->off[iter->i].v != iter->off[iter->i+1].u) { // not adjacent chunks; then seek
 794                                 bam_seek(fp, iter->off[iter->i+1].u, SEEK_SET);
 795                                 iter->curr_off = bam_tell(fp);
 796                         }
 797                         ++iter->i;
 798                 }
 799                 if ((ret = bam_read1(fp, b)) >= 0) {
 800                         iter->curr_off = bam_tell(fp);
 801                         if (b->core.tid != iter->tid || b->core.pos >= iter->end) { // no need to proceed
 802                                 ret = bam_validate1(NULL, b)? -1 : -5; // determine whether end of region or error
 803                                 break;
 804                         }
 805                         else if (is_overlap(iter->beg, iter->end, b)) return ret;
 806                 } else break; // end of file or error
 807         }
 808         iter->finished = 1;
 809         return ret;
 810 }
 811
 812 int bam_fetch(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_fetch_f func)
 813 {
 814         int ret;
 815         bam_iter_t iter;
 816         bam1_t *b;
 817         b = bam_init1();
 818         iter = bam_iter_query(idx, tid, beg, end);
 819         while ((ret = bam_iter_read(fp, iter, b)) >= 0) func(b, data);
 820         bam_iter_destroy(iter);
 821         bam_destroy1(b);
 822         return (ret == -1)? 0 : ret;
 823 }