|
0
|
1 /* tbx.c -- tabix API functions.
|
|
|
2
|
|
|
3 Copyright (C) 2009, 2010, 2012-2014 Genome Research Ltd.
|
|
|
4 Copyright (C) 2010-2012 Broad Institute.
|
|
|
5
|
|
|
6 Author: Heng Li <lh3@sanger.ac.uk>
|
|
|
7
|
|
|
8 Permission is hereby granted, free of charge, to any person obtaining a copy
|
|
|
9 of this software and associated documentation files (the "Software"), to deal
|
|
|
10 in the Software without restriction, including without limitation the rights
|
|
|
11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|
|
12 copies of the Software, and to permit persons to whom the Software is
|
|
|
13 furnished to do so, subject to the following conditions:
|
|
|
14
|
|
|
15 The above copyright notice and this permission notice shall be included in
|
|
|
16 all copies or substantial portions of the Software.
|
|
|
17
|
|
|
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
|
19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
|
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
|
|
|
21 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
|
22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|
|
23 FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
|
|
|
24 DEALINGS IN THE SOFTWARE. */
|
|
|
25
|
|
|
26 #include <stdlib.h>
|
|
|
27 #include <string.h>
|
|
|
28 #include <ctype.h>
|
|
|
29 #include <stdio.h>
|
|
|
30 #include <assert.h>
|
|
|
31 #include "htslib/tbx.h"
|
|
|
32 #include "htslib/bgzf.h"
|
|
|
33
|
|
|
34 #include "htslib/khash.h"
|
|
|
35 KHASH_DECLARE(s2i, kh_cstr_t, int64_t)
|
|
|
36
|
|
|
37 tbx_conf_t tbx_conf_gff = { 0, 1, 4, 5, '#', 0 };
|
|
|
38 tbx_conf_t tbx_conf_bed = { TBX_UCSC, 1, 2, 3, '#', 0 };
|
|
|
39 tbx_conf_t tbx_conf_psltbl = { TBX_UCSC, 15, 17, 18, '#', 0 };
|
|
|
40 tbx_conf_t tbx_conf_sam = { TBX_SAM, 3, 4, 0, '@', 0 };
|
|
|
41 tbx_conf_t tbx_conf_vcf = { TBX_VCF, 1, 2, 0, '#', 0 };
|
|
|
42
|
|
|
43 typedef struct {
|
|
|
44 int64_t beg, end;
|
|
|
45 char *ss, *se;
|
|
|
46 int tid;
|
|
|
47 } tbx_intv_t;
|
|
|
48
|
|
|
49 static inline int get_tid(tbx_t *tbx, const char *ss, int is_add)
|
|
|
50 {
|
|
|
51 khint_t k;
|
|
|
52 khash_t(s2i) *d;
|
|
|
53 if (tbx->dict == 0) tbx->dict = kh_init(s2i);
|
|
|
54 d = (khash_t(s2i)*)tbx->dict;
|
|
|
55 if (is_add) {
|
|
|
56 int absent;
|
|
|
57 k = kh_put(s2i, d, ss, &absent);
|
|
|
58 if (absent) {
|
|
|
59 kh_key(d, k) = strdup(ss);
|
|
|
60 kh_val(d, k) = kh_size(d) - 1;
|
|
|
61 }
|
|
|
62 } else k = kh_get(s2i, d, ss);
|
|
|
63 return k == kh_end(d)? -1 : kh_val(d, k);
|
|
|
64 }
|
|
|
65
|
|
|
66 int tbx_name2id(tbx_t *tbx, const char *ss)
|
|
|
67 {
|
|
|
68 return get_tid(tbx, ss, 0);
|
|
|
69 }
|
|
|
70
|
|
|
71 int tbx_parse1(const tbx_conf_t *conf, int len, char *line, tbx_intv_t *intv)
|
|
|
72 {
|
|
|
73 int i, b = 0, id = 1, ncols = 0;
|
|
|
74 char *s;
|
|
|
75 intv->ss = intv->se = 0; intv->beg = intv->end = -1;
|
|
|
76 for (i = 0; i <= len; ++i) {
|
|
|
77 if (line[i] == '\t' || line[i] == 0) {
|
|
|
78 ++ncols;
|
|
|
79 if (id == conf->sc) {
|
|
|
80 intv->ss = line + b; intv->se = line + i;
|
|
|
81 } else if (id == conf->bc) {
|
|
|
82 // here ->beg is 0-based.
|
|
|
83 intv->beg = intv->end = strtol(line + b, &s, 0);
|
|
|
84 if ( s==line+b ) return -1; // expected int
|
|
|
85 if (!(conf->preset&TBX_UCSC)) --intv->beg;
|
|
|
86 else ++intv->end;
|
|
|
87 if (intv->beg < 0) intv->beg = 0;
|
|
|
88 if (intv->end < 1) intv->end = 1;
|
|
|
89 } else {
|
|
|
90 if ((conf->preset&0xffff) == TBX_GENERIC) {
|
|
|
91 if (id == conf->ec)
|
|
|
92 {
|
|
|
93 intv->end = strtol(line + b, &s, 0);
|
|
|
94 if ( s==line+b ) return -1; // expected int
|
|
|
95 }
|
|
|
96 } else if ((conf->preset&0xffff) == TBX_SAM) {
|
|
|
97 if (id == 6) { // CIGAR
|
|
|
98 int l = 0, op;
|
|
|
99 char *t;
|
|
|
100 for (s = line + b; s < line + i;) {
|
|
|
101 long x = strtol(s, &t, 10);
|
|
|
102 op = toupper(*t);
|
|
|
103 if (op == 'M' || op == 'D' || op == 'N') l += x;
|
|
|
104 s = t + 1;
|
|
|
105 }
|
|
|
106 if (l == 0) l = 1;
|
|
|
107 intv->end = intv->beg + l;
|
|
|
108 }
|
|
|
109 } else if ((conf->preset&0xffff) == TBX_VCF) {
|
|
|
110 if (id == 4) {
|
|
|
111 if (b < i) intv->end = intv->beg + (i - b);
|
|
|
112 } else if (id == 8) { // look for "END="
|
|
|
113 int c = line[i];
|
|
|
114 line[i] = 0;
|
|
|
115 s = strstr(line + b, "END=");
|
|
|
116 if (s == line + b) s += 4;
|
|
|
117 else if (s) {
|
|
|
118 s = strstr(line + b, ";END=");
|
|
|
119 if (s) s += 5;
|
|
|
120 }
|
|
|
121 if (s) intv->end = strtol(s, &s, 0);
|
|
|
122 line[i] = c;
|
|
|
123 }
|
|
|
124 }
|
|
|
125 }
|
|
|
126 b = i + 1;
|
|
|
127 ++id;
|
|
|
128 }
|
|
|
129 }
|
|
|
130 if (intv->ss == 0 || intv->se == 0 || intv->beg < 0 || intv->end < 0) return -1;
|
|
|
131 return 0;
|
|
|
132 }
|
|
|
133
|
|
|
134 static inline int get_intv(tbx_t *tbx, kstring_t *str, tbx_intv_t *intv, int is_add)
|
|
|
135 {
|
|
|
136 if (tbx_parse1(&tbx->conf, str->l, str->s, intv) == 0) {
|
|
|
137 int c = *intv->se;
|
|
|
138 *intv->se = '\0'; intv->tid = get_tid(tbx, intv->ss, is_add); *intv->se = c;
|
|
|
139 return (intv->tid >= 0 && intv->beg >= 0 && intv->end >= 0)? 0 : -1;
|
|
|
140 } else {
|
|
|
141 char *type = NULL;
|
|
|
142 switch (tbx->conf.preset&0xffff)
|
|
|
143 {
|
|
|
144 case TBX_SAM: type = "TBX_SAM"; break;
|
|
|
145 case TBX_VCF: type = "TBX_VCF"; break;
|
|
|
146 case TBX_UCSC: type = "TBX_UCSC"; break;
|
|
|
147 default: type = "TBX_GENERIC"; break;
|
|
|
148 }
|
|
|
149 fprintf(stderr, "[E::%s] failed to parse %s, was wrong -p [type] used?\nThe offending line was: \"%s\"\n", __func__, type, str->s);
|
|
|
150 return -1;
|
|
|
151 }
|
|
|
152 }
|
|
|
153
|
|
|
154 int tbx_readrec(BGZF *fp, void *tbxv, void *sv, int *tid, int *beg, int *end)
|
|
|
155 {
|
|
|
156 tbx_t *tbx = (tbx_t *) tbxv;
|
|
|
157 kstring_t *s = (kstring_t *) sv;
|
|
|
158 int ret;
|
|
|
159 if ((ret = bgzf_getline(fp, '\n', s)) >= 0) {
|
|
|
160 tbx_intv_t intv;
|
|
|
161 get_intv(tbx, s, &intv, 0);
|
|
|
162 *tid = intv.tid; *beg = intv.beg; *end = intv.end;
|
|
|
163 }
|
|
|
164 return ret;
|
|
|
165 }
|
|
|
166
|
|
|
167 void tbx_set_meta(tbx_t *tbx)
|
|
|
168 {
|
|
|
169 int i, l = 0, l_nm;
|
|
|
170 uint32_t x[7];
|
|
|
171 char **name;
|
|
|
172 uint8_t *meta;
|
|
|
173 khint_t k;
|
|
|
174 khash_t(s2i) *d = (khash_t(s2i)*)tbx->dict;
|
|
|
175
|
|
|
176 memcpy(x, &tbx->conf, 24);
|
|
|
177 name = (char**)malloc(sizeof(char*) * kh_size(d));
|
|
|
178 for (k = kh_begin(d), l = 0; k != kh_end(d); ++k) {
|
|
|
179 if (!kh_exist(d, k)) continue;
|
|
|
180 name[kh_val(d, k)] = (char*)kh_key(d, k);
|
|
|
181 l += strlen(kh_key(d, k)) + 1; // +1 to include '\0'
|
|
|
182 }
|
|
|
183 l_nm = x[6] = l;
|
|
|
184 meta = (uint8_t*)malloc(l_nm + 28);
|
|
|
185 if (ed_is_big())
|
|
|
186 for (i = 0; i < 7; ++i)
|
|
|
187 x[i] = ed_swap_4(x[i]);
|
|
|
188 memcpy(meta, x, 28);
|
|
|
189 for (l = 28, i = 0; i < (int)kh_size(d); ++i) {
|
|
|
190 int x = strlen(name[i]) + 1;
|
|
|
191 memcpy(meta + l, name[i], x);
|
|
|
192 l += x;
|
|
|
193 }
|
|
|
194 free(name);
|
|
|
195 hts_idx_set_meta(tbx->idx, l, meta, 0);
|
|
|
196 }
|
|
|
197
|
|
|
198 tbx_t *tbx_index(BGZF *fp, int min_shift, const tbx_conf_t *conf)
|
|
|
199 {
|
|
|
200 tbx_t *tbx;
|
|
|
201 kstring_t str;
|
|
|
202 int ret, first = 0, n_lvls, fmt;
|
|
|
203 int64_t lineno = 0;
|
|
|
204 uint64_t last_off = 0;
|
|
|
205 tbx_intv_t intv;
|
|
|
206
|
|
|
207 str.s = 0; str.l = str.m = 0;
|
|
|
208 tbx = (tbx_t*)calloc(1, sizeof(tbx_t));
|
|
|
209 tbx->conf = *conf;
|
|
|
210 if (min_shift > 0) n_lvls = (TBX_MAX_SHIFT - min_shift + 2) / 3, fmt = HTS_FMT_CSI;
|
|
|
211 else min_shift = 14, n_lvls = 5, fmt = HTS_FMT_TBI;
|
|
|
212 while ((ret = bgzf_getline(fp, '\n', &str)) >= 0) {
|
|
|
213 ++lineno;
|
|
|
214 if (lineno <= tbx->conf.line_skip || str.s[0] == tbx->conf.meta_char) {
|
|
|
215 last_off = bgzf_tell(fp);
|
|
|
216 continue;
|
|
|
217 }
|
|
|
218 if (first == 0) {
|
|
|
219 tbx->idx = hts_idx_init(0, fmt, last_off, min_shift, n_lvls);
|
|
|
220 first = 1;
|
|
|
221 }
|
|
|
222 get_intv(tbx, &str, &intv, 1);
|
|
|
223 ret = hts_idx_push(tbx->idx, intv.tid, intv.beg, intv.end, bgzf_tell(fp), 1);
|
|
|
224 if (ret < 0)
|
|
|
225 {
|
|
|
226 free(str.s);
|
|
|
227 tbx_destroy(tbx);
|
|
|
228 return NULL;
|
|
|
229 }
|
|
|
230 }
|
|
|
231 if ( !tbx->idx ) tbx->idx = hts_idx_init(0, fmt, last_off, min_shift, n_lvls); // empty file
|
|
|
232 if ( !tbx->dict ) tbx->dict = kh_init(s2i);
|
|
|
233 hts_idx_finish(tbx->idx, bgzf_tell(fp));
|
|
|
234 tbx_set_meta(tbx);
|
|
|
235 free(str.s);
|
|
|
236 return tbx;
|
|
|
237 }
|
|
|
238
|
|
|
239 void tbx_destroy(tbx_t *tbx)
|
|
|
240 {
|
|
|
241 khash_t(s2i) *d = (khash_t(s2i)*)tbx->dict;
|
|
|
242 if (d != NULL)
|
|
|
243 {
|
|
|
244 khint_t k;
|
|
|
245 for (k = kh_begin(d); k != kh_end(d); ++k)
|
|
|
246 if (kh_exist(d, k)) free((char*)kh_key(d, k));
|
|
|
247 }
|
|
|
248 hts_idx_destroy(tbx->idx);
|
|
|
249 kh_destroy(s2i, d);
|
|
|
250 free(tbx);
|
|
|
251 }
|
|
|
252
|
|
|
253 int tbx_index_build(const char *fn, int min_shift, const tbx_conf_t *conf)
|
|
|
254 {
|
|
|
255 tbx_t *tbx;
|
|
|
256 BGZF *fp;
|
|
|
257 if ( bgzf_is_bgzf(fn)!=1 ) { fprintf(stderr,"Not a BGZF file: %s\n", fn); return -1; }
|
|
|
258 if ((fp = bgzf_open(fn, "r")) == 0) return -1;
|
|
|
259 if ( !fp->is_compressed ) { bgzf_close(fp); return -1; }
|
|
|
260 tbx = tbx_index(fp, min_shift, conf);
|
|
|
261 bgzf_close(fp);
|
|
|
262 if ( !tbx ) return -1;
|
|
|
263 hts_idx_save(tbx->idx, fn, min_shift > 0? HTS_FMT_CSI : HTS_FMT_TBI);
|
|
|
264 tbx_destroy(tbx);
|
|
|
265 return 0;
|
|
|
266 }
|
|
|
267
|
|
|
268 tbx_t *tbx_index_load(const char *fn)
|
|
|
269 {
|
|
|
270 tbx_t *tbx;
|
|
|
271 uint8_t *meta;
|
|
|
272 char *nm, *p;
|
|
|
273 uint32_t x[7];
|
|
|
274 int l_meta, l_nm;
|
|
|
275 tbx = (tbx_t*)calloc(1, sizeof(tbx_t));
|
|
|
276 tbx->idx = hts_idx_load(fn, HTS_FMT_TBI);
|
|
|
277 if ( !tbx->idx )
|
|
|
278 {
|
|
|
279 free(tbx);
|
|
|
280 return NULL;
|
|
|
281 }
|
|
|
282 meta = hts_idx_get_meta(tbx->idx, &l_meta);
|
|
|
283 if ( !meta )
|
|
|
284 {
|
|
|
285 free(tbx);
|
|
|
286 return NULL;
|
|
|
287 }
|
|
|
288 memcpy(x, meta, 28);
|
|
|
289 memcpy(&tbx->conf, x, 24);
|
|
|
290 p = nm = (char*)meta + 28;
|
|
|
291 l_nm = x[6];
|
|
|
292 for (; p - nm < l_nm; p += strlen(p) + 1) get_tid(tbx, p, 1);
|
|
|
293 return tbx;
|
|
|
294 }
|
|
|
295
|
|
|
296 const char **tbx_seqnames(tbx_t *tbx, int *n)
|
|
|
297 {
|
|
|
298 khash_t(s2i) *d = (khash_t(s2i)*)tbx->dict;
|
|
|
299 if (d == NULL)
|
|
|
300 {
|
|
|
301 *n = 0;
|
|
|
302 return NULL;
|
|
|
303 }
|
|
|
304 int tid, m = kh_size(d);
|
|
|
305 const char **names = (const char**) calloc(m,sizeof(const char*));
|
|
|
306 khint_t k;
|
|
|
307 for (k=kh_begin(d); k<kh_end(d); k++)
|
|
|
308 {
|
|
|
309 if ( !kh_exist(d,k) ) continue;
|
|
|
310 tid = kh_val(d,k);
|
|
|
311 assert( tid<m );
|
|
|
312 names[tid] = kh_key(d,k);
|
|
|
313 }
|
|
|
314 // sanity check: there should be no gaps
|
|
|
315 for (tid=0; tid<m; tid++)
|
|
|
316 assert(names[tid]);
|
|
|
317 *n = m;
|
|
|
318 return names;
|
|
|
319 }
|
|
|
320
|
