if (index2->n < end + 1) index2->n = end + 1;
}
-static void fill_missing(bcf_idx_t *idx)
-{
- int i, j;
- for (i = 0; i < idx->n; ++i) {
- bcf_lidx_t *idx2 = &idx->index2[i];
- for (j = 1; j < idx2->n; ++j)
- if (idx2->offset[j] == 0)
- idx2->offset[j] = idx2->offset[j-1];
- }
-}
-
-bcf_idx_t *bcf_idx_core(bcf_t *bp)
+bcf_idx_t *bcf_idx_core(bcf_t *bp, bcf_hdr_t *h)
{
bcf_idx_t *idx;
int32_t last_coor, last_tid;
b = calloc(1, sizeof(bcf1_t));
str = calloc(1, sizeof(kstring_t));
idx = (bcf_idx_t*)calloc(1, sizeof(bcf_idx_t));
- idx->n = bp->h.n_ref;
- idx->index2 = calloc(bp->h.n_ref, sizeof(bcf_lidx_t));
+ idx->n = h->n_ref;
+ idx->index2 = calloc(h->n_ref, sizeof(bcf_lidx_t));
last_tid = 0xffffffffu;
last_off = bgzf_tell(fp); last_coor = 0xffffffffu;
- while ((ret = bcf_read(bp, b)) > 0) {
+ while ((ret = bcf_read(bp, h, b)) > 0) {
int end, tmp;
if (last_tid != b->tid) { // change of chromosomes
last_tid = b->tid;
last_off = bgzf_tell(fp);
last_coor = b->pos;
}
- fill_missing(idx);
free(str->s); free(str); bcf_destroy(b);
return idx;
}
BGZF *fpidx;
bcf_t *bp;
bcf_idx_t *idx;
+ bcf_hdr_t *h;
if ((bp = bcf_open(fn, "r")) == 0) {
fprintf(stderr, "[bcf_idx_build2] fail to open the BAM file.\n");
return -1;
}
- idx = bcf_idx_core(bp);
+ h = bcf_hdr_read(bp);
+ idx = bcf_idx_core(bp, h);
bcf_close(bp);
if (_fnidx == 0) {
fnidx = (char*)calloc(strlen(fn) + 5, 1);
* retrieve a specified region *
*******************************/
-uint64_t bcf_idx_query(const bcf_idx_t *idx, int tid, int beg, int end)
+uint64_t bcf_idx_query(const bcf_idx_t *idx, int tid, int beg)
{
- uint64_t min_off;
+ uint64_t min_off, *offset;
+ int i;
if (beg < 0) beg = 0;
- if (end < beg) return 0xffffffffffffffffull;
- min_off = (beg>>TAD_LIDX_SHIFT >= idx->index2[tid].n)? idx->index2[tid].offset[idx->index2[tid].n-1]
- : idx->index2[tid].offset[beg>>TAD_LIDX_SHIFT];
- assert(min_off);
+ offset = idx->index2[tid].offset;
+ for (i = beg>>TAD_LIDX_SHIFT; i < idx->index2[tid].n && offset[i] == 0; ++i);
+ min_off = (i == idx->index2[tid].n)? offset[idx->index2[tid].n-1] : offset[i];
return min_off;
}