4 #include "bam_endian.h"
11 #define TAD_LIDX_SHIFT 13
27 static inline void insert_offset2(bcf_lidx_t *index2, int _beg, int _end, uint64_t offset)
30 beg = _beg >> TAD_LIDX_SHIFT;
31 end = (_end - 1) >> TAD_LIDX_SHIFT;
32 if (index2->m < end + 1) {
33 int old_m = index2->m;
35 kroundup32(index2->m);
36 index2->offset = (uint64_t*)realloc(index2->offset, index2->m * 8);
37 memset(index2->offset + old_m, 0, 8 * (index2->m - old_m));
40 if (index2->offset[beg] == 0) index2->offset[beg] = offset;
42 for (i = beg; i <= end; ++i)
43 if (index2->offset[i] == 0) index2->offset[i] = offset;
45 if (index2->n < end + 1) index2->n = end + 1;
48 static void fill_missing(bcf_idx_t *idx)
51 for (i = 0; i < idx->n; ++i) {
52 bcf_lidx_t *idx2 = &idx->index2[i];
53 for (j = 1; j < idx2->n; ++j)
54 if (idx2->offset[j] == 0)
55 idx2->offset[j] = idx2->offset[j-1];
59 bcf_idx_t *bcf_idx_core(bcf_t *bp, bcf_hdr_t *h)
62 int32_t last_coor, last_tid;
69 b = calloc(1, sizeof(bcf1_t));
70 str = calloc(1, sizeof(kstring_t));
71 idx = (bcf_idx_t*)calloc(1, sizeof(bcf_idx_t));
73 idx->index2 = calloc(h->n_ref, sizeof(bcf_lidx_t));
75 last_tid = 0xffffffffu;
76 last_off = bgzf_tell(fp); last_coor = 0xffffffffu;
77 while ((ret = bcf_read(bp, h, b)) > 0) {
79 if (last_tid != b->tid) { // change of chromosomes
81 } else if (last_coor > b->pos) {
82 fprintf(stderr, "[bcf_idx_core] the input is out of order\n");
83 free(str->s); free(str); free(idx); bcf_destroy(b);
87 end = b->pos + (tmp > 0? tmp : 1);
88 insert_offset2(&idx->index2[b->tid], b->pos, end, last_off);
89 last_off = bgzf_tell(fp);
93 free(str->s); free(str); bcf_destroy(b);
97 void bcf_idx_destroy(bcf_idx_t *idx)
100 if (idx == 0) return;
101 for (i = 0; i < idx->n; ++i) free(idx->index2[i].offset);
110 void bcf_idx_save(const bcf_idx_t *idx, BGZF *fp)
113 ti_is_be = bam_is_big_endian();
114 bgzf_write(fp, "BCI\4", 4);
117 bgzf_write(fp, bam_swap_endian_4p(&x), 4);
118 } else bgzf_write(fp, &idx->n, 4);
119 for (i = 0; i < idx->n; ++i) {
120 bcf_lidx_t *index2 = idx->index2 + i;
121 // write linear index (index2)
124 bgzf_write(fp, bam_swap_endian_4p(&x), 4);
125 } else bgzf_write(fp, &index2->n, 4);
126 if (ti_is_be) { // big endian
128 for (x = 0; (int)x < index2->n; ++x)
129 bam_swap_endian_8p(&index2->offset[x]);
130 bgzf_write(fp, index2->offset, 8 * index2->n);
131 for (x = 0; (int)x < index2->n; ++x)
132 bam_swap_endian_8p(&index2->offset[x]);
133 } else bgzf_write(fp, index2->offset, 8 * index2->n);
137 static bcf_idx_t *bcf_idx_load_core(BGZF *fp)
142 ti_is_be = bam_is_big_endian();
144 fprintf(stderr, "[%s] fail to load index.\n", __func__);
147 bgzf_read(fp, magic, 4);
148 if (strncmp(magic, "BCI\4", 4)) {
149 fprintf(stderr, "[%s] wrong magic number.\n", __func__);
152 idx = (bcf_idx_t*)calloc(1, sizeof(bcf_idx_t));
153 bgzf_read(fp, &idx->n, 4);
154 if (ti_is_be) bam_swap_endian_4p(&idx->n);
155 idx->index2 = (bcf_lidx_t*)calloc(idx->n, sizeof(bcf_lidx_t));
156 for (i = 0; i < idx->n; ++i) {
157 bcf_lidx_t *index2 = idx->index2 + i;
159 bgzf_read(fp, &index2->n, 4);
160 if (ti_is_be) bam_swap_endian_4p(&index2->n);
161 index2->m = index2->n;
162 index2->offset = (uint64_t*)calloc(index2->m, 8);
163 bgzf_read(fp, index2->offset, index2->n * 8);
165 for (j = 0; j < index2->n; ++j) bam_swap_endian_8p(&index2->offset[j]);
170 bcf_idx_t *bcf_idx_load_local(const char *fnidx)
173 fp = bgzf_open(fnidx, "r");
175 bcf_idx_t *idx = bcf_idx_load_core(fp);
182 static void download_from_remote(const char *url)
184 const int buf_size = 1 * 1024 * 1024;
190 if (strstr(url, "ftp://") != url && strstr(url, "http://") != url) return;
192 for (fn = (char*)url + l - 1; fn >= url; --fn)
193 if (*fn == '/') break;
194 ++fn; // fn now points to the file name
195 fp_remote = knet_open(url, "r");
196 if (fp_remote == 0) {
197 fprintf(stderr, "[download_from_remote] fail to open remote file.\n");
200 if ((fp = fopen(fn, "w")) == 0) {
201 fprintf(stderr, "[download_from_remote] fail to create file in the working directory.\n");
202 knet_close(fp_remote);
205 buf = (uint8_t*)calloc(buf_size, 1);
206 while ((l = knet_read(fp_remote, buf, buf_size)) != 0)
207 fwrite(buf, 1, l, fp);
210 knet_close(fp_remote);
213 static void download_from_remote(const char *url)
219 static char *get_local_version(const char *fn)
222 char *fnidx = (char*)calloc(strlen(fn) + 5, 1);
223 strcat(strcpy(fnidx, fn), ".bci");
224 if ((strstr(fnidx, "ftp://") == fnidx || strstr(fnidx, "http://") == fnidx)) {
226 int l = strlen(fnidx);
227 for (p = fnidx + l - 1; p >= fnidx; --p)
228 if (*p == '/') break;
229 url = fnidx; fnidx = strdup(p + 1);
230 if (stat(fnidx, &sbuf) == 0) {
234 fprintf(stderr, "[%s] downloading the index file...\n", __func__);
235 download_from_remote(url);
238 if (stat(fnidx, &sbuf) == 0) return fnidx;
239 free(fnidx); return 0;
242 bcf_idx_t *bcf_idx_load(const char *fn)
245 char *fname = get_local_version(fn);
246 if (fname == 0) return 0;
247 idx = bcf_idx_load_local(fname);
252 int bcf_idx_build2(const char *fn, const char *_fnidx)
259 if ((bp = bcf_open(fn, "r")) == 0) {
260 fprintf(stderr, "[bcf_idx_build2] fail to open the BAM file.\n");
263 h = bcf_hdr_read(bp);
264 idx = bcf_idx_core(bp, h);
267 fnidx = (char*)calloc(strlen(fn) + 5, 1);
268 strcpy(fnidx, fn); strcat(fnidx, ".bci");
269 } else fnidx = strdup(_fnidx);
270 fpidx = bgzf_open(fnidx, "w");
272 fprintf(stderr, "[bcf_idx_build2] fail to create the index file.\n");
276 bcf_idx_save(idx, fpidx);
277 bcf_idx_destroy(idx);
283 int bcf_idx_build(const char *fn)
285 return bcf_idx_build2(fn, 0);
288 /********************************************
289 * parse a region in the format chr:beg-end *
290 ********************************************/
292 int bcf_parse_region(void *str2id, const char *str, int *tid, int *begin, int *end)
297 p = s = (char*)malloc(l+1);
298 /* squeeze out "," */
299 for (i = k = 0; i != l; ++i)
300 if (str[i] != ',' && !isspace(str[i])) s[k++] = str[i];
302 for (i = 0; i != k; ++i) if (s[i] == ':') break;
304 if ((*tid = bcf_str2id(str2id, s)) < 0) {
308 if (i == k) { /* dump the whole sequence */
309 *begin = 0; *end = 1<<29; free(s);
312 for (p = s + i + 1; i != k; ++i) if (s[i] == '-') break;
318 if (*begin > 0) --*begin;
320 if (*begin > *end) return -1;
324 /*******************************
325 * retrieve a specified region *
326 *******************************/
328 uint64_t bcf_idx_query(const bcf_idx_t *idx, int tid, int beg, int end)
331 if (beg < 0) beg = 0;
332 if (end < beg) return 0xffffffffffffffffull;
333 min_off = (beg>>TAD_LIDX_SHIFT >= idx->index2[tid].n)? idx->index2[tid].offset[idx->index2[tid].n-1]
334 : idx->index2[tid].offset[beg>>TAD_LIDX_SHIFT];
339 int bcf_main_index(int argc, char *argv[])
342 fprintf(stderr, "Usage: bcftools index <in.bcf>\n");
345 bcf_idx_build(argv[1]);