Mercurial > repos > youngkim > ezbamqc
view ezBAMQC/src/htslib/hts.c @ 18:494b5cd02238
bash script
author | youngkim |
---|---|
date | Wed, 30 Mar 2016 13:39:05 -0400 |
parents | dfa3745e5fd8 |
children |
line wrap: on
line source
/* hts.c -- format-neutral I/O, indexing, and iterator API functions. Copyright (C) 2008, 2009, 2012-2015 Genome Research Ltd. Copyright (C) 2012, 2013 Broad Institute. Author: Heng Li <lh3@sanger.ac.uk> Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #include <zlib.h> #include <ctype.h> #include <stdio.h> #include <string.h> #include <stdlib.h> #include <limits.h> #include <fcntl.h> #include <errno.h> #include <sys/stat.h> #include "htslib/bgzf.h" #include "htslib/hts.h" #include "cram/cram.h" #include "htslib/hfile.h" #include "version.h" #include "htslib/kseq.h" #define KS_BGZF 1 #if KS_BGZF // bgzf now supports gzip-compressed files, the gzFile branch can be removed KSTREAM_INIT2(, BGZF*, bgzf_read, 65536) #else KSTREAM_INIT2(, gzFile, gzread, 16384) #endif #include "htslib/khash.h" KHASH_INIT2(s2i,, kh_cstr_t, int64_t, 1, kh_str_hash_func, kh_str_hash_equal) int hts_verbose = 3; const char *hts_version() { return HTS_VERSION; } const unsigned char seq_nt16_table[256] = { 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 1, 2, 4, 8, 15,15,15,15, 15,15,15,15, 15, 0 /*=*/,15,15, 15, 1,14, 2, 13,15,15, 4, 11,15,15,12, 15, 3,15,15, 15,15, 5, 6, 8,15, 7, 9, 15,10,15,15, 15,15,15,15, 15, 1,14, 2, 13,15,15, 4, 11,15,15,12, 15, 3,15,15, 15,15, 5, 6, 8,15, 7, 9, 15,10,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15, 15,15,15,15 }; const char seq_nt16_str[] = "=ACMGRSVTWYHKDBN"; const int seq_nt16_int[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 }; /********************** *** Basic file I/O *** **********************/ static enum htsFormatCategory format_category(enum htsExactFormat fmt) { switch (fmt) { case bam: case sam: case cram: return sequence_data; case vcf: case bcf: return variant_data; case bai: case crai: case csi: case gzi: case tbi: return index_file; case bed: return region_list; case unknown_format: case binary_format: case text_format: case format_maximum: break; } return unknown_category; } // Decompress up to ten or so bytes by peeking at the file, which must be // positioned at the start of a GZIP block. static size_t decompress_peek(hFILE *fp, unsigned char *dest, size_t destsize) { // Typically at most a couple of hundred bytes of input are required // to get a few bytes of output from inflate(), so hopefully this buffer // size suffices in general. unsigned char buffer[512]; z_stream zs; ssize_t npeek = hpeek(fp, buffer, sizeof buffer); if (npeek < 0) return 0; zs.zalloc = NULL; zs.zfree = NULL; zs.next_in = buffer; zs.avail_in = npeek; zs.next_out = dest; zs.avail_out = destsize; if (inflateInit2(&zs, 31) != Z_OK) return 0; while (zs.total_out < destsize) if (inflate(&zs, Z_SYNC_FLUSH) != Z_OK) break; destsize = zs.total_out; inflateEnd(&zs); return destsize; } // Parse "x.y" text, taking care because the string is not NUL-terminated // and filling in major/minor only when the digits are followed by a delimiter, // so we don't misread "1.10" as "1.1" due to reaching the end of the buffer. static void parse_version(htsFormat *fmt, const unsigned char *u, const unsigned char *ulim) { const char *str = (const char *) u; const char *slim = (const char *) ulim; const char *s; fmt->version.major = fmt->version.minor = -1; for (s = str; s < slim; s++) if (!isdigit(*s)) break; if (s < slim) { fmt->version.major = atoi(str); if (*s == '.') { str = &s[1]; for (s = str; s < slim; s++) if (!isdigit(*s)) break; if (s < slim) fmt->version.minor = atoi(str); } else fmt->version.minor = 0; } } int hts_detect_format(hFILE *hfile, htsFormat *fmt) { unsigned char s[21]; ssize_t len = hpeek(hfile, s, 18); if (len < 0) return -1; if (len >= 2 && s[0] == 0x1f && s[1] == 0x8b) { // The stream is either gzip-compressed or BGZF-compressed. // Determine which, and decompress the first few bytes. fmt->compression = (len >= 18 && (s[3] & 4) && memcmp(&s[12], "BC\2\0", 4) == 0)? bgzf : gzip; len = decompress_peek(hfile, s, sizeof s); } else { fmt->compression = no_compression; len = hpeek(hfile, s, sizeof s); } if (len < 0) return -1; fmt->compression_level = -1; fmt->specific = NULL; if (len >= 6 && memcmp(s,"CRAM",4) == 0 && s[4]>=1 && s[4]<=3 && s[5]<=1) { fmt->category = sequence_data; fmt->format = cram; fmt->version.major = s[4], fmt->version.minor = s[5]; fmt->compression = custom; return 0; } else if (len >= 4 && s[3] <= '\4') { if (memcmp(s, "BAM\1", 4) == 0) { fmt->category = sequence_data; fmt->format = bam; // TODO Decompress enough to pick version from @HD-VN header fmt->version.major = 1, fmt->version.minor = -1; return 0; } else if (memcmp(s, "BAI\1", 4) == 0) { fmt->category = index_file; fmt->format = bai; fmt->version.major = -1, fmt->version.minor = -1; return 0; } else if (memcmp(s, "BCF\4", 4) == 0) { fmt->category = variant_data; fmt->format = bcf; fmt->version.major = 1, fmt->version.minor = -1; return 0; } else if (memcmp(s, "BCF\2", 4) == 0) { fmt->category = variant_data; fmt->format = bcf; fmt->version.major = s[3]; fmt->version.minor = (len >= 5 && s[4] <= 2)? s[4] : 0; return 0; } else if (memcmp(s, "CSI\1", 4) == 0) { fmt->category = index_file; fmt->format = csi; fmt->version.major = 1, fmt->version.minor = -1; return 0; } else if (memcmp(s, "TBI\1", 4) == 0) { fmt->category = index_file; fmt->format = tbi; fmt->version.major = -1, fmt->version.minor = -1; return 0; } } else if (len >= 16 && memcmp(s, "##fileformat=VCF", 16) == 0) { fmt->category = variant_data; fmt->format = vcf; if (len >= 21 && s[16] == 'v') parse_version(fmt, &s[17], &s[len]); else fmt->version.major = fmt->version.minor = -1; return 0; } else if (len >= 4 && s[0] == '@' && (memcmp(s, "@HD\t", 4) == 0 || memcmp(s, "@SQ\t", 4) == 0 || memcmp(s, "@RG\t", 4) == 0 || memcmp(s, "@PG\t", 4) == 0)) { fmt->category = sequence_data; fmt->format = sam; // @HD-VN is not guaranteed to be the first tag, but then @HD is // not guaranteed to be present at all... if (len >= 9 && memcmp(s, "@HD\tVN:", 7) == 0) parse_version(fmt, &s[7], &s[len]); else fmt->version.major = 1, fmt->version.minor = -1; return 0; } else { // Various possibilities for tab-delimited text: // .crai (gzipped tab-delimited six columns: seqid 5*number) // .bed ([3..12] tab-delimited columns) // .bedpe (>= 10 tab-delimited columns) // .sam (tab-delimited >= 11 columns: seqid number seqid...) // FIXME For now, assume it's SAM fmt->category = sequence_data; fmt->format = sam; fmt->version.major = 1, fmt->version.minor = -1; return 0; } fmt->category = unknown_category; fmt->format = unknown_format; fmt->version.major = fmt->version.minor = -1; fmt->compression = no_compression; return 0; } char *hts_format_description(const htsFormat *format) { kstring_t str = { 0, 0, NULL }; switch (format->format) { case sam: kputs("SAM", &str); break; case bam: kputs("BAM", &str); break; case cram: kputs("CRAM", &str); break; case vcf: kputs("VCF", &str); break; case bcf: if (format->version.major == 1) kputs("Legacy BCF", &str); else kputs("BCF", &str); break; case bai: kputs("BAI", &str); break; case crai: kputs("CRAI", &str); break; case csi: kputs("CSI", &str); break; case tbi: kputs("Tabix", &str); break; default: kputs("unknown", &str); break; } if (format->version.major >= 0) { kputs(" version ", &str); kputw(format->version.major, &str); if (format->version.minor >= 0) { kputc('.', &str); kputw(format->version.minor, &str); } } switch (format->compression) { case custom: kputs(" compressed", &str); break; case gzip: kputs(" gzip-compressed", &str); break; case bgzf: switch (format->format) { case bam: case bcf: case csi: case tbi: // These are by definition BGZF, so just use the generic term kputs(" compressed", &str); break; default: kputs(" BGZF-compressed", &str); break; } break; default: break; } switch (format->category) { case sequence_data: kputs(" sequence", &str); break; case variant_data: kputs(" variant calling", &str); break; case index_file: kputs(" index", &str); break; case region_list: kputs(" genomic region", &str); break; default: break; } if (format->compression == no_compression) switch (format->format) { case sam: case crai: case vcf: case bed: kputs(" text", &str); break; default: kputs(" data", &str); break; } else kputs(" data", &str); return ks_release(&str); } htsFile *hts_open(const char *fn, const char *mode) { htsFile *fp = NULL; hFILE *hfile = hopen(fn, mode); if (hfile == NULL) goto error; fp = hts_hopen(hfile, fn, mode); if (fp == NULL) goto error; return fp; error: if (hts_verbose >= 2) fprintf(stderr, "[E::%s] fail to open file '%s'\n", __func__, fn); if (hfile) hclose_abruptly(hfile); return NULL; } htsFile *hts_hopen(struct hFILE *hfile, const char *fn, const char *mode) { htsFile *fp = (htsFile*)calloc(1, sizeof(htsFile)); if (fp == NULL) goto error; fp->fn = strdup(fn); fp->is_be = ed_is_big(); if (strchr(mode, 'r')) { if (hts_detect_format(hfile, &fp->format) < 0) goto error; } else if (strchr(mode, 'w') || strchr(mode, 'a')) { htsFormat *fmt = &fp->format; fp->is_write = 1; if (strchr(mode, 'b')) fmt->format = binary_format; else if (strchr(mode, 'c')) fmt->format = cram; else fmt->format = text_format; if (strchr(mode, 'z')) fmt->compression = bgzf; else if (strchr(mode, 'g')) fmt->compression = gzip; else if (strchr(mode, 'u')) fmt->compression = no_compression; else { // No compression mode specified, set to the default for the format switch (fmt->format) { case binary_format: fmt->compression = bgzf; break; case cram: fmt->compression = custom; break; case text_format: fmt->compression = no_compression; break; default: abort(); } } // Fill in category (if determinable; e.g. 'b' could be BAM or BCF) fmt->category = format_category(fmt->format); fmt->version.major = fmt->version.minor = -1; fmt->compression_level = -1; fmt->specific = NULL; } else goto error; switch (fp->format.format) { case binary_format: case bam: case bcf: fp->fp.bgzf = bgzf_hopen(hfile, mode); if (fp->fp.bgzf == NULL) goto error; fp->is_bin = 1; break; case cram: fp->fp.cram = cram_dopen(hfile, fn, mode); if (fp->fp.cram == NULL) goto error; if (!fp->is_write) cram_set_option(fp->fp.cram, CRAM_OPT_DECODE_MD, 1); fp->is_cram = 1; break; case text_format: case sam: case vcf: if (!fp->is_write) { #if KS_BGZF BGZF *gzfp = bgzf_hopen(hfile, mode); #else // TODO Implement gzip hFILE adaptor hclose(hfile); // This won't work, especially for stdin gzFile gzfp = strcmp(fn, "-")? gzopen(fn, "rb") : gzdopen(fileno(stdin), "rb"); #endif if (gzfp) fp->fp.voidp = ks_init(gzfp); else goto error; } else if (fp->format.compression != no_compression) { fp->fp.bgzf = bgzf_hopen(hfile, mode); if (fp->fp.bgzf == NULL) goto error; } else fp->fp.hfile = hfile; break; default: goto error; } return fp; error: if (hts_verbose >= 2) fprintf(stderr, "[E::%s] fail to open file '%s'\n", __func__, fn); if (fp) { free(fp->fn); free(fp->fn_aux); free(fp); } return NULL; } int hts_close(htsFile *fp) { int ret, save; switch (fp->format.format) { case binary_format: case bam: case bcf: ret = bgzf_close(fp->fp.bgzf); break; case cram: if (!fp->is_write) { switch (cram_eof(fp->fp.cram)) { case 0: fprintf(stderr, "[E::%s] Failed to decode sequence.\n", __func__); return -1; case 2: fprintf(stderr, "[W::%s] EOF marker is absent. The input is probably truncated.\n", __func__); break; default: /* case 1, expected EOF */ break; } } ret = cram_close(fp->fp.cram); break; case text_format: case sam: case vcf: if (!fp->is_write) { #if KS_BGZF BGZF *gzfp = ((kstream_t*)fp->fp.voidp)->f; ret = bgzf_close(gzfp); #else gzFile gzfp = ((kstream_t*)fp->fp.voidp)->f; ret = gzclose(gzfp); #endif ks_destroy((kstream_t*)fp->fp.voidp); } else if (fp->format.compression != no_compression) ret = bgzf_close(fp->fp.bgzf); else ret = hclose(fp->fp.hfile); break; default: ret = -1; break; } save = errno; free(fp->fn); free(fp->fn_aux); free(fp->line.s); free(fp); errno = save; return ret; } const htsFormat *hts_get_format(htsFile *fp) { return fp? &fp->format : NULL; } int hts_set_opt(htsFile *fp, enum cram_option opt, ...) { int r; va_list args; if (fp->format.format != cram) return 0; va_start(args, opt); r = cram_set_voption(fp->fp.cram, opt, args); va_end(args); return r; } int hts_set_threads(htsFile *fp, int n) { if (fp->format.compression == bgzf) { return bgzf_mt(fp->fp.bgzf, n, 256); } else if (fp->format.format == cram) { return hts_set_opt(fp, CRAM_OPT_NTHREADS, n); } else return 0; } int hts_set_fai_filename(htsFile *fp, const char *fn_aux) { free(fp->fn_aux); if (fn_aux) { fp->fn_aux = strdup(fn_aux); if (fp->fn_aux == NULL) return -1; } else fp->fn_aux = NULL; if (fp->format.format == cram) cram_set_option(fp->fp.cram, CRAM_OPT_REFERENCE, fp->fn_aux); return 0; } // For VCF/BCF backward sweeper. Not exposing these functions because their // future is uncertain. Things will probably have to change with hFILE... BGZF *hts_get_bgzfp(htsFile *fp) { if ( fp->is_bin ) return fp->fp.bgzf; else return ((kstream_t*)fp->fp.voidp)->f; } int hts_useek(htsFile *fp, long uoffset, int where) { if ( fp->is_bin ) return bgzf_useek(fp->fp.bgzf, uoffset, where); else { ks_rewind((kstream_t*)fp->fp.voidp); ((kstream_t*)fp->fp.voidp)->seek_pos = uoffset; return bgzf_useek(((kstream_t*)fp->fp.voidp)->f, uoffset, where); } } long hts_utell(htsFile *fp) { if ( fp->is_bin ) return bgzf_utell(fp->fp.bgzf); else return ((kstream_t*)fp->fp.voidp)->seek_pos; } int hts_getline(htsFile *fp, int delimiter, kstring_t *str) { int ret, dret; ret = ks_getuntil((kstream_t*)fp->fp.voidp, delimiter, str, &dret); ++fp->lineno; return ret; } char **hts_readlist(const char *string, int is_file, int *_n) { int m = 0, n = 0, dret; char **s = 0; if ( is_file ) { #if KS_BGZF BGZF *fp = bgzf_open(string, "r"); #else gzFile fp = gzopen(string, "r"); #endif if ( !fp ) return NULL; kstream_t *ks; kstring_t str; str.s = 0; str.l = str.m = 0; ks = ks_init(fp); while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) { if (str.l == 0) continue; n++; hts_expand(char*,n,m,s); s[n-1] = strdup(str.s); } ks_destroy(ks); #if KS_BGZF bgzf_close(fp); #else gzclose(fp); #endif free(str.s); } else { const char *q = string, *p = string; while ( 1 ) { if (*p == ',' || *p == 0) { n++; hts_expand(char*,n,m,s); s[n-1] = (char*)calloc(p - q + 1, 1); strncpy(s[n-1], q, p - q); q = p + 1; } if ( !*p ) break; p++; } } s = (char**)realloc(s, n * sizeof(char*)); *_n = n; return s; } char **hts_readlines(const char *fn, int *_n) { int m = 0, n = 0, dret; char **s = 0; #if KS_BGZF BGZF *fp = bgzf_open(fn, "r"); #else gzFile fp = gzopen(fn, "r"); #endif if ( fp ) { // read from file kstream_t *ks; kstring_t str; str.s = 0; str.l = str.m = 0; ks = ks_init(fp); while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) { if (str.l == 0) continue; if (m == n) { m = m? m<<1 : 16; s = (char**)realloc(s, m * sizeof(char*)); } s[n++] = strdup(str.s); } ks_destroy(ks); #if KS_BGZF bgzf_close(fp); #else gzclose(fp); #endif s = (char**)realloc(s, n * sizeof(char*)); free(str.s); } else if (*fn == ':') { // read from string const char *q, *p; for (q = p = fn + 1;; ++p) if (*p == ',' || *p == 0) { if (m == n) { m = m? m<<1 : 16; s = (char**)realloc(s, m * sizeof(char*)); } s[n] = (char*)calloc(p - q + 1, 1); strncpy(s[n++], q, p - q); q = p + 1; if (*p == 0) break; } } else return 0; s = (char**)realloc(s, n * sizeof(char*)); *_n = n; return s; } // DEPRECATED: To be removed in a future HTSlib release int hts_file_type(const char *fname) { int len = strlen(fname); if ( !strcasecmp(".vcf.gz",fname+len-7) ) return FT_VCF_GZ; if ( !strcasecmp(".vcf",fname+len-4) ) return FT_VCF; if ( !strcasecmp(".bcf",fname+len-4) ) return FT_BCF_GZ; if ( !strcmp("-",fname) ) return FT_STDIN; hFILE *f = hopen(fname, "r"); if (f == NULL) return 0; htsFormat fmt; if (hts_detect_format(f, &fmt) < 0) { hclose_abruptly(f); return 0; } if (hclose(f) < 0) return 0; switch (fmt.format) { case vcf: return (fmt.compression == no_compression)? FT_VCF : FT_VCF_GZ; case bcf: return (fmt.compression == no_compression)? FT_BCF : FT_BCF_GZ; default: return 0; } } /**************** *** Indexing *** ****************/ #define HTS_MIN_MARKER_DIST 0x10000 // Finds the special meta bin // ((1<<(3 * n_lvls + 3)) - 1) / 7 + 1 #define META_BIN(idx) ((idx)->n_bins + 1) #define pair64_lt(a,b) ((a).u < (b).u) #include "htslib/ksort.h" KSORT_INIT(_off, hts_pair64_t, pair64_lt) typedef struct { int32_t m, n; uint64_t loff; hts_pair64_t *list; } bins_t; #include "htslib/khash.h" KHASH_MAP_INIT_INT(bin, bins_t) typedef khash_t(bin) bidx_t; typedef struct { int32_t n, m; uint64_t *offset; } lidx_t; struct __hts_idx_t { int fmt, min_shift, n_lvls, n_bins; uint32_t l_meta; int32_t n, m; uint64_t n_no_coor; bidx_t **bidx; lidx_t *lidx; uint8_t *meta; struct { uint32_t last_bin, save_bin; int last_coor, last_tid, save_tid, finished; uint64_t last_off, save_off; uint64_t off_beg, off_end; uint64_t n_mapped, n_unmapped; } z; // keep internal states }; static inline void insert_to_b(bidx_t *b, int bin, uint64_t beg, uint64_t end) { khint_t k; bins_t *l; int absent; k = kh_put(bin, b, bin, &absent); l = &kh_value(b, k); if (absent) { l->m = 1; l->n = 0; l->list = (hts_pair64_t*)calloc(l->m, sizeof(hts_pair64_t)); } if (l->n == l->m) { l->m <<= 1; l->list = (hts_pair64_t*)realloc(l->list, l->m * sizeof(hts_pair64_t)); } l->list[l->n].u = beg; l->list[l->n++].v = end; } static inline void insert_to_l(lidx_t *l, int64_t _beg, int64_t _end, uint64_t offset, int min_shift) { int i, beg, end; beg = _beg >> min_shift; end = (_end - 1) >> min_shift; if (l->m < end + 1) { int old_m = l->m; l->m = end + 1; kroundup32(l->m); l->offset = (uint64_t*)realloc(l->offset, l->m * sizeof(uint64_t)); memset(l->offset + old_m, 0xff, 8 * (l->m - old_m)); // fill l->offset with (uint64_t)-1 } if (beg == end) { // to save a loop in this case if (l->offset[beg] == (uint64_t)-1) l->offset[beg] = offset; } else { for (i = beg; i <= end; ++i) if (l->offset[i] == (uint64_t)-1) l->offset[i] = offset; } if (l->n < end + 1) l->n = end + 1; } hts_idx_t *hts_idx_init(int n, int fmt, uint64_t offset0, int min_shift, int n_lvls) { hts_idx_t *idx; idx = (hts_idx_t*)calloc(1, sizeof(hts_idx_t)); if (idx == NULL) return NULL; idx->fmt = fmt; idx->min_shift = min_shift; idx->n_lvls = n_lvls; idx->n_bins = ((1<<(3 * n_lvls + 3)) - 1) / 7; idx->z.save_bin = idx->z.save_tid = idx->z.last_tid = idx->z.last_bin = 0xffffffffu; idx->z.save_off = idx->z.last_off = idx->z.off_beg = idx->z.off_end = offset0; idx->z.last_coor = 0xffffffffu; if (n) { idx->n = idx->m = n; idx->bidx = (bidx_t**)calloc(n, sizeof(bidx_t*)); if (idx->bidx == NULL) { free(idx); return NULL; } idx->lidx = (lidx_t*) calloc(n, sizeof(lidx_t)); if (idx->lidx == NULL) { free(idx->bidx); free(idx); return NULL; } } return idx; } static void update_loff(hts_idx_t *idx, int i, int free_lidx) { bidx_t *bidx = idx->bidx[i]; lidx_t *lidx = &idx->lidx[i]; khint_t k; int l; uint64_t offset0 = 0; if (bidx) { k = kh_get(bin, bidx, META_BIN(idx)); if (k != kh_end(bidx)) offset0 = kh_val(bidx, k).list[0].u; for (l = 0; l < lidx->n && lidx->offset[l] == (uint64_t)-1; ++l) lidx->offset[l] = offset0; } else l = 1; for (; l < lidx->n; ++l) // fill missing values if (lidx->offset[l] == (uint64_t)-1) lidx->offset[l] = lidx->offset[l-1]; if (bidx == 0) return; for (k = kh_begin(bidx); k != kh_end(bidx); ++k) // set loff if (kh_exist(bidx, k)) { if ( kh_key(bidx, k) < idx->n_bins ) { int bot_bin = hts_bin_bot(kh_key(bidx, k), idx->n_lvls); // disable linear index if bot_bin out of bounds kh_val(bidx, k).loff = bot_bin < lidx->n ? lidx->offset[bot_bin] : 0; } else kh_val(bidx, k).loff = 0; } if (free_lidx) { free(lidx->offset); lidx->m = lidx->n = 0; lidx->offset = 0; } } static void compress_binning(hts_idx_t *idx, int i) { bidx_t *bidx = idx->bidx[i]; khint_t k; int l, m; if (bidx == 0) return; // merge a bin to its parent if the bin is too small for (l = idx->n_lvls; l > 0; --l) { unsigned start = hts_bin_first(l); for (k = kh_begin(bidx); k != kh_end(bidx); ++k) { bins_t *p, *q; if (!kh_exist(bidx, k) || kh_key(bidx, k) >= idx->n_bins || kh_key(bidx, k) < start) continue; p = &kh_value(bidx, k); if (l < idx->n_lvls && p->n > 1) ks_introsort(_off, p->n, p->list); if ((p->list[p->n - 1].v>>16) - (p->list[0].u>>16) < HTS_MIN_MARKER_DIST) { khint_t kp; kp = kh_get(bin, bidx, hts_bin_parent(kh_key(bidx, k))); if (kp == kh_end(bidx)) continue; q = &kh_val(bidx, kp); if (q->n + p->n > q->m) { q->m = q->n + p->n; kroundup32(q->m); q->list = (hts_pair64_t*)realloc(q->list, q->m * sizeof(hts_pair64_t)); } memcpy(q->list + q->n, p->list, p->n * sizeof(hts_pair64_t)); q->n += p->n; free(p->list); kh_del(bin, bidx, k); } } } k = kh_get(bin, bidx, 0); if (k != kh_end(bidx)) ks_introsort(_off, kh_val(bidx, k).n, kh_val(bidx, k).list); // merge adjacent chunks that start from the same BGZF block for (k = kh_begin(bidx); k != kh_end(bidx); ++k) { bins_t *p; if (!kh_exist(bidx, k) || kh_key(bidx, k) >= idx->n_bins) continue; p = &kh_value(bidx, k); for (l = 1, m = 0; l < p->n; ++l) { if (p->list[m].v>>16 >= p->list[l].u>>16) { if (p->list[m].v < p->list[l].v) p->list[m].v = p->list[l].v; } else p->list[++m] = p->list[l]; } p->n = m + 1; } } void hts_idx_finish(hts_idx_t *idx, uint64_t final_offset) { int i; if (idx == NULL || idx->z.finished) return; // do not run this function on an empty index or multiple times if (idx->z.save_tid >= 0) { insert_to_b(idx->bidx[idx->z.save_tid], idx->z.save_bin, idx->z.save_off, final_offset); insert_to_b(idx->bidx[idx->z.save_tid], META_BIN(idx), idx->z.off_beg, final_offset); insert_to_b(idx->bidx[idx->z.save_tid], META_BIN(idx), idx->z.n_mapped, idx->z.n_unmapped); } for (i = 0; i < idx->n; ++i) { update_loff(idx, i, (idx->fmt == HTS_FMT_CSI)); compress_binning(idx, i); } idx->z.finished = 1; } int hts_idx_push(hts_idx_t *idx, int tid, int beg, int end, uint64_t offset, int is_mapped) { int bin; if (tid<0) beg = -1, end = 0; if (tid >= idx->m) { // enlarge the index int32_t oldm = idx->m; idx->m = idx->m? idx->m<<1 : 2; idx->bidx = (bidx_t**)realloc(idx->bidx, idx->m * sizeof(bidx_t*)); idx->lidx = (lidx_t*) realloc(idx->lidx, idx->m * sizeof(lidx_t)); memset(&idx->bidx[oldm], 0, (idx->m - oldm) * sizeof(bidx_t*)); memset(&idx->lidx[oldm], 0, (idx->m - oldm) * sizeof(lidx_t)); } if (idx->n < tid + 1) idx->n = tid + 1; if (idx->z.finished) return 0; if (idx->z.last_tid != tid || (idx->z.last_tid >= 0 && tid < 0)) { // change of chromosome if ( tid>=0 && idx->n_no_coor ) { if (hts_verbose >= 1) fprintf(stderr,"[E::%s] NO_COOR reads not in a single block at the end %d %d\n", __func__, tid,idx->z.last_tid); return -1; } if (tid>=0 && idx->bidx[tid] != 0) { if (hts_verbose >= 1) fprintf(stderr, "[E::%s] chromosome blocks not continuous\n", __func__); return -1; } idx->z.last_tid = tid; idx->z.last_bin = 0xffffffffu; } else if (tid >= 0 && idx->z.last_coor > beg) { // test if positions are out of order if (hts_verbose >= 1) fprintf(stderr, "[E::%s] unsorted positions\n", __func__); return -1; } if ( tid>=0 ) { if (idx->bidx[tid] == 0) idx->bidx[tid] = kh_init(bin); if ( is_mapped) insert_to_l(&idx->lidx[tid], beg, end, idx->z.last_off, idx->min_shift); // last_off points to the start of the current record } else idx->n_no_coor++; bin = hts_reg2bin(beg, end, idx->min_shift, idx->n_lvls); if ((int)idx->z.last_bin != bin) { // then possibly write the binning index if (idx->z.save_bin != 0xffffffffu) // save_bin==0xffffffffu only happens to the first record insert_to_b(idx->bidx[idx->z.save_tid], idx->z.save_bin, idx->z.save_off, idx->z.last_off); if (idx->z.last_bin == 0xffffffffu && idx->z.save_bin != 0xffffffffu) { // change of chr; keep meta information idx->z.off_end = idx->z.last_off; insert_to_b(idx->bidx[idx->z.save_tid], META_BIN(idx), idx->z.off_beg, idx->z.off_end); insert_to_b(idx->bidx[idx->z.save_tid], META_BIN(idx), idx->z.n_mapped, idx->z.n_unmapped); idx->z.n_mapped = idx->z.n_unmapped = 0; idx->z.off_beg = idx->z.off_end; } idx->z.save_off = idx->z.last_off; idx->z.save_bin = idx->z.last_bin = bin; idx->z.save_tid = tid; } if (is_mapped) ++idx->z.n_mapped; else ++idx->z.n_unmapped; idx->z.last_off = offset; idx->z.last_coor = beg; return 0; } void hts_idx_destroy(hts_idx_t *idx) { khint_t k; int i; if (idx == 0) return; // For HTS_FMT_CRAI, idx actually points to a different type -- see sam.c if (idx->fmt == HTS_FMT_CRAI) { free(idx); return; } for (i = 0; i < idx->m; ++i) { bidx_t *bidx = idx->bidx[i]; free(idx->lidx[i].offset); if (bidx == 0) continue; for (k = kh_begin(bidx); k != kh_end(bidx); ++k) if (kh_exist(bidx, k)) free(kh_value(bidx, k).list); kh_destroy(bin, bidx); } free(idx->bidx); free(idx->lidx); free(idx->meta); free(idx); } static inline long idx_read(int is_bgzf, void *fp, void *buf, long l) { if (is_bgzf) return bgzf_read((BGZF*)fp, buf, l); else return (long)fread(buf, 1, l, (FILE*)fp); } static inline long idx_write(int is_bgzf, void *fp, const void *buf, long l) { if (is_bgzf) return bgzf_write((BGZF*)fp, buf, l); else return (long)fwrite(buf, 1, l, (FILE*)fp); } static inline void swap_bins(bins_t *p) { int i; for (i = 0; i < p->n; ++i) { ed_swap_8p(&p->list[i].u); ed_swap_8p(&p->list[i].v); } } static void hts_idx_save_core(const hts_idx_t *idx, void *fp, int fmt) { int32_t i, size, is_be; int is_bgzf = (fmt != HTS_FMT_BAI); is_be = ed_is_big(); if (is_be) { uint32_t x = idx->n; idx_write(is_bgzf, fp, ed_swap_4p(&x), 4); } else idx_write(is_bgzf, fp, &idx->n, 4); if (fmt == HTS_FMT_TBI && idx->l_meta) idx_write(is_bgzf, fp, idx->meta, idx->l_meta); for (i = 0; i < idx->n; ++i) { khint_t k; bidx_t *bidx = idx->bidx[i]; lidx_t *lidx = &idx->lidx[i]; // write binning index size = bidx? kh_size(bidx) : 0; if (is_be) { // big endian uint32_t x = size; idx_write(is_bgzf, fp, ed_swap_4p(&x), 4); } else idx_write(is_bgzf, fp, &size, 4); if (bidx == 0) goto write_lidx; for (k = kh_begin(bidx); k != kh_end(bidx); ++k) { bins_t *p; if (!kh_exist(bidx, k)) continue; p = &kh_value(bidx, k); if (is_be) { // big endian uint32_t x; x = kh_key(bidx, k); idx_write(is_bgzf, fp, ed_swap_4p(&x), 4); if (fmt == HTS_FMT_CSI) { uint64_t y = kh_val(bidx, k).loff; idx_write(is_bgzf, fp, ed_swap_4p(&y), 8); } x = p->n; idx_write(is_bgzf, fp, ed_swap_4p(&x), 4); swap_bins(p); idx_write(is_bgzf, fp, p->list, 16 * p->n); swap_bins(p); } else { idx_write(is_bgzf, fp, &kh_key(bidx, k), 4); if (fmt == HTS_FMT_CSI) idx_write(is_bgzf, fp, &kh_val(bidx, k).loff, 8); //int j;for(j=0;j<p->n;++j)fprintf(stderr,"%d,%llx,%d,%llx:%llx\n",kh_key(bidx,k),kh_val(bidx, k).loff,j,p->list[j].u,p->list[j].v); idx_write(is_bgzf, fp, &p->n, 4); idx_write(is_bgzf, fp, p->list, p->n << 4); } } write_lidx: if (fmt != HTS_FMT_CSI) { if (is_be) { int32_t x = lidx->n; idx_write(is_bgzf, fp, ed_swap_4p(&x), 4); for (x = 0; x < lidx->n; ++x) ed_swap_8p(&lidx->offset[x]); idx_write(is_bgzf, fp, lidx->offset, lidx->n << 3); for (x = 0; x < lidx->n; ++x) ed_swap_8p(&lidx->offset[x]); } else { idx_write(is_bgzf, fp, &lidx->n, 4); idx_write(is_bgzf, fp, lidx->offset, lidx->n << 3); } } } if (is_be) { // write the number of reads without coordinates uint64_t x = idx->n_no_coor; idx_write(is_bgzf, fp, &x, 8); } else idx_write(is_bgzf, fp, &idx->n_no_coor, 8); } void hts_idx_save(const hts_idx_t *idx, const char *fn, int fmt) { char *fnidx; fnidx = (char*)calloc(1, strlen(fn) + 5); strcpy(fnidx, fn); if (fmt == HTS_FMT_CSI) { BGZF *fp; uint32_t x[3]; int is_be, i; is_be = ed_is_big(); fp = bgzf_open(strcat(fnidx, ".csi"), "w"); bgzf_write(fp, "CSI\1", 4); x[0] = idx->min_shift; x[1] = idx->n_lvls; x[2] = idx->l_meta; if (is_be) { for (i = 0; i < 3; ++i) bgzf_write(fp, ed_swap_4p(&x[i]), 4); } else bgzf_write(fp, &x, 12); if (idx->l_meta) bgzf_write(fp, idx->meta, idx->l_meta); hts_idx_save_core(idx, fp, HTS_FMT_CSI); bgzf_close(fp); } else if (fmt == HTS_FMT_TBI) { BGZF *fp; fp = bgzf_open(strcat(fnidx, ".tbi"), "w"); bgzf_write(fp, "TBI\1", 4); hts_idx_save_core(idx, fp, HTS_FMT_TBI); bgzf_close(fp); } else if (fmt == HTS_FMT_BAI) { FILE *fp; fp = fopen(strcat(fnidx, ".bai"), "w"); fwrite("BAI\1", 1, 4, fp); hts_idx_save_core(idx, fp, HTS_FMT_BAI); fclose(fp); } else abort(); free(fnidx); } static int hts_idx_load_core(hts_idx_t *idx, void *fp, int fmt) { int32_t i, n, is_be; int is_bgzf = (fmt != HTS_FMT_BAI); is_be = ed_is_big(); if (idx == NULL) return -4; for (i = 0; i < idx->n; ++i) { bidx_t *h; lidx_t *l = &idx->lidx[i]; uint32_t key; int j, absent; bins_t *p; h = idx->bidx[i] = kh_init(bin); if (idx_read(is_bgzf, fp, &n, 4) != 4) return -1; if (is_be) ed_swap_4p(&n); for (j = 0; j < n; ++j) { khint_t k; if (idx_read(is_bgzf, fp, &key, 4) != 4) return -1; if (is_be) ed_swap_4p(&key); k = kh_put(bin, h, key, &absent); if (absent <= 0) return -3; // Duplicate bin number p = &kh_val(h, k); if (fmt == HTS_FMT_CSI) { if (idx_read(is_bgzf, fp, &p->loff, 8) != 8) return -1; if (is_be) ed_swap_8p(&p->loff); } else p->loff = 0; if (idx_read(is_bgzf, fp, &p->n, 4) != 4) return -1; if (is_be) ed_swap_4p(&p->n); p->m = p->n; p->list = (hts_pair64_t*)malloc(p->m * sizeof(hts_pair64_t)); if (p->list == NULL) return -2; if (idx_read(is_bgzf, fp, p->list, p->n<<4) != p->n<<4) return -1; if (is_be) swap_bins(p); } if (fmt != HTS_FMT_CSI) { // load linear index int j; if (idx_read(is_bgzf, fp, &l->n, 4) != 4) return -1; if (is_be) ed_swap_4p(&l->n); l->m = l->n; l->offset = (uint64_t*)malloc(l->n * sizeof(uint64_t)); if (l->offset == NULL) return -2; if (idx_read(is_bgzf, fp, l->offset, l->n << 3) != l->n << 3) return -1; if (is_be) for (j = 0; j < l->n; ++j) ed_swap_8p(&l->offset[j]); for (j = 1; j < l->n; ++j) // fill missing values; may happen given older samtools and tabix if (l->offset[j] == 0) l->offset[j] = l->offset[j-1]; update_loff(idx, i, 1); } } if (idx_read(is_bgzf, fp, &idx->n_no_coor, 8) != 8) idx->n_no_coor = 0; if (is_be) ed_swap_8p(&idx->n_no_coor); return 0; } hts_idx_t *hts_idx_load_local(const char *fn, int fmt) { uint8_t magic[4]; int i, is_be; hts_idx_t *idx = NULL; is_be = ed_is_big(); if (fmt == HTS_FMT_CSI) { BGZF *fp; uint32_t x[3], n; uint8_t *meta = 0; if ((fp = bgzf_open(fn, "r")) == 0) return NULL; if (bgzf_read(fp, magic, 4) != 4) goto csi_fail; if (memcmp(magic, "CSI\1", 4) != 0) goto csi_fail; if (bgzf_read(fp, x, 12) != 12) goto csi_fail; if (is_be) for (i = 0; i < 3; ++i) ed_swap_4p(&x[i]); if (x[2]) { if ((meta = (uint8_t*)malloc(x[2])) == NULL) goto csi_fail; if (bgzf_read(fp, meta, x[2]) != x[2]) goto csi_fail; } if (bgzf_read(fp, &n, 4) != 4) goto csi_fail; if (is_be) ed_swap_4p(&n); if ((idx = hts_idx_init(n, fmt, 0, x[0], x[1])) == NULL) goto csi_fail; idx->l_meta = x[2]; idx->meta = meta; meta = NULL; if (hts_idx_load_core(idx, fp, HTS_FMT_CSI) < 0) goto csi_fail; bgzf_close(fp); return idx; csi_fail: bgzf_close(fp); hts_idx_destroy(idx); free(meta); return NULL; } else if (fmt == HTS_FMT_TBI) { BGZF *fp; uint32_t x[8]; if ((fp = bgzf_open(fn, "r")) == 0) return NULL; if (bgzf_read(fp, magic, 4) != 4) goto tbi_fail; if (memcmp(magic, "TBI\1", 4) != 0) goto tbi_fail; if (bgzf_read(fp, x, 32) != 32) goto tbi_fail; if (is_be) for (i = 0; i < 8; ++i) ed_swap_4p(&x[i]); if ((idx = hts_idx_init(x[0], fmt, 0, 14, 5)) == NULL) goto tbi_fail; idx->l_meta = 28 + x[7]; if ((idx->meta = (uint8_t*)malloc(idx->l_meta)) == NULL) goto tbi_fail; memcpy(idx->meta, &x[1], 28); if (bgzf_read(fp, idx->meta + 28, x[7]) != x[7]) goto tbi_fail; if (hts_idx_load_core(idx, fp, HTS_FMT_TBI) < 0) goto tbi_fail; bgzf_close(fp); return idx; tbi_fail: bgzf_close(fp); hts_idx_destroy(idx); return NULL; } else if (fmt == HTS_FMT_BAI) { uint32_t n; FILE *fp; if ((fp = fopen(fn, "rb")) == 0) return NULL; if (fread(magic, 1, 4, fp) != 4) goto bai_fail; if (memcmp(magic, "BAI\1", 4) != 0) goto bai_fail; if (fread(&n, 4, 1, fp) != 1) goto bai_fail; if (is_be) ed_swap_4p(&n); idx = hts_idx_init(n, fmt, 0, 14, 5); if (hts_idx_load_core(idx, fp, HTS_FMT_BAI) < 0) goto bai_fail; fclose(fp); return idx; bai_fail: fclose(fp); hts_idx_destroy(idx); return NULL; } else abort(); } void hts_idx_set_meta(hts_idx_t *idx, int l_meta, uint8_t *meta, int is_copy) { if (idx->meta) free(idx->meta); idx->l_meta = l_meta; if (is_copy) { idx->meta = (uint8_t*)malloc(l_meta); memcpy(idx->meta, meta, l_meta); } else idx->meta = meta; } uint8_t *hts_idx_get_meta(hts_idx_t *idx, int *l_meta) { *l_meta = idx->l_meta; return idx->meta; } const char **hts_idx_seqnames(const hts_idx_t *idx, int *n, hts_id2name_f getid, void *hdr) { if ( !idx->n ) { *n = 0; return NULL; } int tid = 0, i; const char **names = (const char**) calloc(idx->n,sizeof(const char*)); for (i=0; i<idx->n; i++) { bidx_t *bidx = idx->bidx[i]; if ( !bidx ) continue; names[tid++] = getid(hdr,i); } *n = tid; return names; } int hts_idx_get_stat(const hts_idx_t* idx, int tid, uint64_t* mapped, uint64_t* unmapped) { if ( idx->fmt == HTS_FMT_CRAI ) { *mapped = 0; *unmapped = 0; return -1; } bidx_t *h = idx->bidx[tid]; khint_t k = kh_get(bin, h, META_BIN(idx)); if (k != kh_end(h)) { *mapped = kh_val(h, k).list[1].u; *unmapped = kh_val(h, k).list[1].v; return 0; } else { *mapped = 0; *unmapped = 0; return -1; } } uint64_t hts_idx_get_n_no_coor(const hts_idx_t* idx) { return idx->n_no_coor; } /**************** *** Iterator *** ****************/ static inline int reg2bins(int64_t beg, int64_t end, hts_itr_t *itr, int min_shift, int n_lvls) { int l, t, s = min_shift + (n_lvls<<1) + n_lvls; if (beg >= end) return 0; if (end >= 1LL<<s) end = 1LL<<s; for (--end, l = 0, t = 0; l <= n_lvls; s -= 3, t += 1<<((l<<1)+l), ++l) { int b, e, n, i; b = t + (beg>>s); e = t + (end>>s); n = e - b + 1; if (itr->bins.n + n > itr->bins.m) { itr->bins.m = itr->bins.n + n; kroundup32(itr->bins.m); itr->bins.a = (int*)realloc(itr->bins.a, sizeof(int) * itr->bins.m); } for (i = b; i <= e; ++i) itr->bins.a[itr->bins.n++] = i; } return itr->bins.n; } hts_itr_t *hts_itr_query(const hts_idx_t *idx, int tid, int beg, int end, hts_readrec_func *readrec) { int i, n_off, l, bin; hts_pair64_t *off; khint_t k; bidx_t *bidx; uint64_t min_off; hts_itr_t *iter = 0; if (tid < 0) { int finished0 = 0; uint64_t off0 = (uint64_t)-1; khint_t k; switch (tid) { case HTS_IDX_START: // Find the smallest offset, note that sequence ids may not be ordered sequentially for (i=0; i<idx->n; i++) { bidx = idx->bidx[i]; k = kh_get(bin, bidx, META_BIN(idx)); if (k == kh_end(bidx)) continue; if ( off0 > kh_val(bidx, k).list[0].u ) off0 = kh_val(bidx, k).list[0].u; } if ( off0==(uint64_t)-1 && idx->n_no_coor ) off0 = 0; // only no-coor reads in this bam break; case HTS_IDX_NOCOOR: if ( idx->n>0 ) { bidx = idx->bidx[idx->n - 1]; k = kh_get(bin, bidx, META_BIN(idx)); if (k != kh_end(bidx)) off0 = kh_val(bidx, k).list[0].v; } if ( off0==(uint64_t)-1 && idx->n_no_coor ) off0 = 0; // only no-coor reads in this bam break; case HTS_IDX_REST: off0 = 0; break; case HTS_IDX_NONE: finished0 = 1; off0 = 0; break; default: return 0; } if (off0 != (uint64_t)-1) { iter = (hts_itr_t*)calloc(1, sizeof(hts_itr_t)); iter->read_rest = 1; iter->finished = finished0; iter->curr_off = off0; iter->readrec = readrec; return iter; } else return 0; } if (beg < 0) beg = 0; if (end < beg) return 0; if (tid >= idx->n || (bidx = idx->bidx[tid]) == NULL) return 0; iter = (hts_itr_t*)calloc(1, sizeof(hts_itr_t)); iter->tid = tid, iter->beg = beg, iter->end = end; iter->i = -1; iter->readrec = readrec; // compute min_off bin = hts_bin_first(idx->n_lvls) + (beg>>idx->min_shift); do { int first; k = kh_get(bin, bidx, bin); if (k != kh_end(bidx)) break; first = (hts_bin_parent(bin)<<3) + 1; if (bin > first) --bin; else bin = hts_bin_parent(bin); } while (bin); if (bin == 0) k = kh_get(bin, bidx, bin); min_off = k != kh_end(bidx)? kh_val(bidx, k).loff : 0; // retrieve bins reg2bins(beg, end, iter, idx->min_shift, idx->n_lvls); for (i = n_off = 0; i < iter->bins.n; ++i) if ((k = kh_get(bin, bidx, iter->bins.a[i])) != kh_end(bidx)) n_off += kh_value(bidx, k).n; if (n_off == 0) return iter; off = (hts_pair64_t*)calloc(n_off, sizeof(hts_pair64_t)); for (i = n_off = 0; i < iter->bins.n; ++i) { if ((k = kh_get(bin, bidx, iter->bins.a[i])) != kh_end(bidx)) { int j; bins_t *p = &kh_value(bidx, k); for (j = 0; j < p->n; ++j) if (p->list[j].v > min_off) off[n_off++] = p->list[j]; } } if (n_off == 0) { free(off); return iter; } ks_introsort(_off, n_off, off); // resolve completely contained adjacent blocks for (i = 1, l = 0; i < n_off; ++i) if (off[l].v < off[i].v) off[++l] = off[i]; n_off = l + 1; // resolve overlaps between adjacent blocks; this may happen due to the merge in indexing for (i = 1; i < n_off; ++i) if (off[i-1].v >= off[i].u) off[i-1].v = off[i].u; // merge adjacent blocks for (i = 1, l = 0; i < n_off; ++i) { if (off[l].v>>16 == off[i].u>>16) off[l].v = off[i].v; else off[++l] = off[i]; } n_off = l + 1; iter->n_off = n_off; iter->off = off; return iter; } void hts_itr_destroy(hts_itr_t *iter) { if (iter) { free(iter->off); free(iter->bins.a); free(iter); } } const char *hts_parse_reg(const char *s, int *beg, int *end) { int i, k, l, name_end; *beg = *end = -1; name_end = l = strlen(s); // determine the sequence name for (i = l - 1; i >= 0; --i) if (s[i] == ':') break; // look for colon from the end if (i >= 0) name_end = i; if (name_end < l) { // check if this is really the end int n_hyphen = 0; for (i = name_end + 1; i < l; ++i) { if (s[i] == '-') ++n_hyphen; else if (!isdigit(s[i]) && s[i] != ',') break; } if (i < l || n_hyphen > 1) name_end = l; // malformated region string; then take str as the name } // parse the interval if (name_end < l) { char *tmp; tmp = (char*)alloca(l - name_end + 1); for (i = name_end + 1, k = 0; i < l; ++i) if (s[i] != ',') tmp[k++] = s[i]; tmp[k] = 0; if ((*beg = strtol(tmp, &tmp, 10) - 1) < 0) *beg = 0; *end = *tmp? strtol(tmp + 1, &tmp, 10) : INT_MAX; if (*beg > *end) name_end = l; } if (name_end == l) *beg = 0, *end = INT_MAX; return s + name_end; } hts_itr_t *hts_itr_querys(const hts_idx_t *idx, const char *reg, hts_name2id_f getid, void *hdr, hts_itr_query_func *itr_query, hts_readrec_func *readrec) { int tid, beg, end; char *q, *tmp; if (strcmp(reg, ".") == 0) return itr_query(idx, HTS_IDX_START, 0, 0, readrec); else if (strcmp(reg, "*") != 0) { q = (char*)hts_parse_reg(reg, &beg, &end); tmp = (char*)alloca(q - reg + 1); strncpy(tmp, reg, q - reg); tmp[q - reg] = 0; if ((tid = getid(hdr, tmp)) < 0) tid = getid(hdr, reg); if (tid < 0) return 0; return itr_query(idx, tid, beg, end, readrec); } else return itr_query(idx, HTS_IDX_NOCOOR, 0, 0, readrec); } int hts_itr_next(BGZF *fp, hts_itr_t *iter, void *r, void *data) { int ret, tid, beg, end; if (iter == NULL || iter->finished) return -1; if (iter->read_rest) { if (iter->curr_off) { // seek to the start bgzf_seek(fp, iter->curr_off, SEEK_SET); iter->curr_off = 0; // only seek once } ret = iter->readrec(fp, data, r, &tid, &beg, &end); if (ret < 0) iter->finished = 1; iter->curr_tid = tid; iter->curr_beg = beg; iter->curr_end = end; return ret; } if (iter->off == 0) return -1; for (;;) { if (iter->curr_off == 0 || iter->curr_off >= iter->off[iter->i].v) { // then jump to the next chunk if (iter->i == iter->n_off - 1) { ret = -1; break; } // no more chunks if (iter->i < 0 || iter->off[iter->i].v != iter->off[iter->i+1].u) { // not adjacent chunks; then seek bgzf_seek(fp, iter->off[iter->i+1].u, SEEK_SET); iter->curr_off = bgzf_tell(fp); } ++iter->i; } if ((ret = iter->readrec(fp, data, r, &tid, &beg, &end)) >= 0) { iter->curr_off = bgzf_tell(fp); if (tid != iter->tid || beg >= iter->end) { // no need to proceed ret = -1; break; } else if (end > iter->beg && iter->end > beg) { iter->curr_tid = tid; iter->curr_beg = beg; iter->curr_end = end; return ret; } } else break; // end of file or error } iter->finished = 1; return ret; } /********************** *** Retrieve index *** **********************/ static char *test_and_fetch(const char *fn) { FILE *fp; if (hisremote(fn)) { const int buf_size = 1 * 1024 * 1024; hFILE *fp_remote; uint8_t *buf; int l; const char *p; for (p = fn + strlen(fn) - 1; p >= fn; --p) if (*p == '/') break; ++p; // p now points to the local file name // Attempt to open local file first if ((fp = fopen((char*)p, "rb")) != 0) { fclose(fp); return (char*)p; } // Attempt to open remote file. Stay quiet on failure, it is OK to fail when trying first .csi then .tbi index. if ((fp_remote = hopen(fn, "r")) == 0) return 0; if ((fp = fopen(p, "w")) == 0) { if (hts_verbose >= 1) fprintf(stderr, "[E::%s] fail to create file '%s' in the working directory\n", __func__, p); hclose_abruptly(fp_remote); return 0; } if (hts_verbose >= 3) fprintf(stderr, "[M::%s] downloading file '%s' to local directory\n", __func__, fn); buf = (uint8_t*)calloc(buf_size, 1); while ((l = hread(fp_remote, buf, buf_size)) > 0) fwrite(buf, 1, l, fp); free(buf); fclose(fp); if (hclose(fp_remote) != 0) fprintf(stderr, "[E::%s] fail to close remote file '%s'\n", __func__, fn); return (char*)p; } else { if ((fp = fopen(fn, "rb")) == 0) return 0; fclose(fp); return (char*)fn; } } char *hts_idx_getfn(const char *fn, const char *ext) { int i, l_fn, l_ext; char *fnidx, *ret; l_fn = strlen(fn); l_ext = strlen(ext); fnidx = (char*)calloc(l_fn + l_ext + 1, 1); strcpy(fnidx, fn); strcpy(fnidx + l_fn, ext); if ((ret = test_and_fetch(fnidx)) == 0) { for (i = l_fn - 1; i > 0; --i) if (fnidx[i] == '.') break; strcpy(fnidx + i, ext); ret = test_and_fetch(fnidx); } if (ret == 0) { free(fnidx); return 0; } l_fn = strlen(ret); memmove(fnidx, ret, l_fn + 1); return fnidx; } hts_idx_t *hts_idx_load(const char *fn, int fmt) { char *fnidx; hts_idx_t *idx; fnidx = hts_idx_getfn(fn, ".csi"); if (fnidx) fmt = HTS_FMT_CSI; else fnidx = hts_idx_getfn(fn, fmt == HTS_FMT_BAI? ".bai" : ".tbi"); if (fnidx == 0) return 0; // Check that the index file is up to date, the main file might have changed struct stat stat_idx,stat_main; if ( !stat(fn, &stat_main) && !stat(fnidx, &stat_idx) ) { if ( stat_idx.st_mtime < stat_main.st_mtime ) fprintf(stderr, "Warning: The index file is older than the data file: %s\n", fnidx); } idx = hts_idx_load_local(fnidx, fmt); free(fnidx); return idx; }