Line | Count | Source |
1 | | /* tbx.c -- tabix API functions. |
2 | | |
3 | | Copyright (C) 2009, 2010, 2012-2015, 2017-2020, 2022-2023, 2025-2026 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 | | #define HTS_BUILDING_LIBRARY // Enables HTSLIB_EXPORT, see htslib/hts_defs.h |
27 | | #include <config.h> |
28 | | |
29 | | #include <stdlib.h> |
30 | | #include <string.h> |
31 | | #include <stdio.h> |
32 | | #include <assert.h> |
33 | | #include <errno.h> |
34 | | #include "htslib/tbx.h" |
35 | | #include "htslib/bgzf.h" |
36 | | #include "htslib/hts_alloc.h" |
37 | | #include "htslib/hts_endian.h" |
38 | | #include "hts_internal.h" |
39 | | #include "bgzf_internal.h" |
40 | | |
41 | | #include "htslib/khash.h" |
42 | | KHASH_DECLARE(s2i, kh_cstr_t, int64_t) |
43 | | |
44 | | HTSLIB_EXPORT |
45 | | const tbx_conf_t tbx_conf_gff = { 0, 1, 4, 5, '#', 0 }; |
46 | | |
47 | | HTSLIB_EXPORT |
48 | | const tbx_conf_t tbx_conf_bed = { TBX_UCSC, 1, 2, 3, '#', 0 }; |
49 | | |
50 | | HTSLIB_EXPORT |
51 | | const tbx_conf_t tbx_conf_psltbl = { TBX_UCSC, 15, 17, 18, '#', 0 }; |
52 | | |
53 | | HTSLIB_EXPORT |
54 | | const tbx_conf_t tbx_conf_sam = { TBX_SAM, 3, 4, 0, '@', 0 }; |
55 | | |
56 | | HTSLIB_EXPORT |
57 | | const tbx_conf_t tbx_conf_vcf = { TBX_VCF, 1, 2, 0, '#', 0 }; |
58 | | const tbx_conf_t tbx_conf_gaf = { TBX_GAF, 1, 6, 0, '#', 0 }; |
59 | | |
60 | | typedef struct { |
61 | | int64_t beg, end; |
62 | | char *ss, *se; |
63 | | int tid; |
64 | | } tbx_intv_t; |
65 | | |
66 | | static inline int get_tid(tbx_t *tbx, const char *ss, int is_add) |
67 | 0 | { |
68 | 0 | khint_t k; |
69 | 0 | khash_t(s2i) *d; |
70 | 0 | if ((tbx->conf.preset&0xffff) == TBX_GAF) return(0); |
71 | 0 | if (tbx->dict == 0) tbx->dict = kh_init(s2i); |
72 | 0 | if (!tbx->dict) return -1; // Out of memory |
73 | 0 | d = (khash_t(s2i)*)tbx->dict; |
74 | 0 | if (is_add) { |
75 | 0 | int absent; |
76 | 0 | k = kh_put(s2i, d, ss, &absent); |
77 | 0 | if (absent < 0) { |
78 | 0 | return -1; // Out of memory |
79 | 0 | } else if (absent) { |
80 | 0 | char *ss_dup = strdup(ss); |
81 | 0 | if (ss_dup) { |
82 | 0 | kh_key(d, k) = ss_dup; |
83 | 0 | kh_val(d, k) = kh_size(d) - 1; |
84 | 0 | } else { |
85 | 0 | kh_del(s2i, d, k); |
86 | 0 | return -1; // Out of memory |
87 | 0 | } |
88 | 0 | } |
89 | 0 | } else k = kh_get(s2i, d, ss); |
90 | 0 | return k == kh_end(d)? -1 : kh_val(d, k); |
91 | 0 | } |
92 | | |
93 | | int tbx_name2id(tbx_t *tbx, const char *ss) |
94 | 0 | { |
95 | 0 | return get_tid(tbx, ss, 0); |
96 | 0 | } |
97 | | |
98 | | int tbx_parse1(const tbx_conf_t *conf, size_t len, char *line, tbx_intv_t *intv) |
99 | 0 | { |
100 | 0 | size_t i, b = 0; |
101 | 0 | int id = 1, getlen = 0, alcnt = 0, use_svlen = 0, lenpos = -1; |
102 | 0 | char *s, *t; |
103 | 0 | uint8_t svlenals[8192]; |
104 | 0 | int64_t reflen = 0, svlen = 0, fmtlen = 0, tmp = 0; |
105 | |
|
106 | 0 | intv->ss = intv->se = 0; intv->beg = intv->end = -1; |
107 | 0 | for (i = 0; i <= len; ++i) { |
108 | 0 | if (line[i] == '\t' || line[i] == 0) { |
109 | 0 | if (id == conf->sc) { |
110 | 0 | intv->ss = line + b; intv->se = line + i; |
111 | 0 | } else if (id == conf->bc) { |
112 | | // here ->beg is 0-based. |
113 | 0 | if ((conf->preset&0xffff) == TBX_GAF){ |
114 | | // if gaf find the smallest and largest node id |
115 | 0 | char *t; |
116 | 0 | int64_t nodeid = -1; |
117 | 0 | for (s = line + b + 1; s < line + i;) { |
118 | 0 | nodeid = strtoll(s, &t, 0); |
119 | 0 | if(intv->beg == -1){ |
120 | 0 | intv->beg = intv->end = nodeid; |
121 | 0 | } else { |
122 | 0 | if(nodeid < intv->beg){ |
123 | 0 | intv->beg = nodeid; |
124 | 0 | } |
125 | |
|
126 | 0 | if(nodeid > intv->end){ |
127 | 0 | intv->end = nodeid; |
128 | 0 | } |
129 | 0 | } |
130 | 0 | s = t + 1; |
131 | 0 | } |
132 | 0 | } else { |
133 | 0 | intv->beg = strtoll(line + b, &s, 0); |
134 | |
|
135 | 0 | if (conf->bc <= conf->ec) // don't overwrite an already set end point |
136 | 0 | intv->end = intv->beg; |
137 | |
|
138 | 0 | if ( s==line+b ) return -1; // expected int |
139 | | |
140 | 0 | if (!(conf->preset&TBX_UCSC)) |
141 | 0 | --intv->beg; |
142 | 0 | else if (conf->bc <= conf->ec) |
143 | 0 | ++intv->end; |
144 | |
|
145 | 0 | if (intv->beg < 0) { |
146 | 0 | hts_log_warning("Coordinate <= 0 detected. " |
147 | 0 | "Did you forget to use the -0 option?"); |
148 | 0 | intv->beg = 0; |
149 | 0 | } |
150 | 0 | if (intv->end < 1) intv->end = 1; |
151 | 0 | } |
152 | 0 | } else { |
153 | 0 | if ((conf->preset&0xffff) == TBX_GENERIC) { |
154 | 0 | if (id == conf->ec) |
155 | 0 | { |
156 | 0 | intv->end = strtoll(line + b, &s, 0); |
157 | 0 | if ( s==line+b ) return -1; // expected int |
158 | 0 | } |
159 | 0 | } else if ((conf->preset&0xffff) == TBX_SAM) { |
160 | 0 | if (id == 6) { // CIGAR |
161 | 0 | int l = 0; |
162 | 0 | char *t; |
163 | 0 | for (s = line + b; s < line + i;) { |
164 | 0 | long x = strtol(s, &t, 10); |
165 | 0 | char op = toupper_c(*t); |
166 | 0 | if (op == 'M' || op == 'D' || op == 'N') l += x; |
167 | 0 | s = t + 1; |
168 | 0 | } |
169 | 0 | if (l == 0) l = 1; |
170 | 0 | intv->end = intv->beg + l; |
171 | 0 | } |
172 | 0 | } else if ((conf->preset&0xffff) == TBX_VCF) { |
173 | 0 | if (id == 4) { //ref allele |
174 | 0 | if (b < i) intv->end = intv->beg + (i - b); |
175 | 0 | ++alcnt; |
176 | 0 | reflen = i - b; |
177 | 0 | } if (id == 5) { //alt allele |
178 | 0 | int lastbyte = 0, c = line[i]; |
179 | 0 | svlenals[lastbyte] = 0; |
180 | 0 | line[i] = 0; |
181 | 0 | s = line + b; |
182 | 0 | do { |
183 | 0 | t = strchr(s, ','); |
184 | 0 | if (alcnt >> 3 != lastbyte) { //initialize insals |
185 | 0 | lastbyte = alcnt >> 3; |
186 | 0 | svlenals[lastbyte] = 0; |
187 | 0 | } |
188 | 0 | ++alcnt; |
189 | 0 | if (t) { |
190 | 0 | *t = 0; |
191 | 0 | } |
192 | 0 | if (svlen_on_ref_for_vcf_alt(s, -1)) { |
193 | | // Need to check SVLEN for this ALT |
194 | 0 | svlenals[lastbyte] |= 1 << ((alcnt - 1) & 7); |
195 | 0 | use_svlen = 1; |
196 | 0 | } else if (!strcmp("<*>", s) || |
197 | 0 | !strcmp("<NON_REF>", s)) { //note gvcf |
198 | 0 | getlen = 1; |
199 | 0 | } |
200 | 0 | if (t) { |
201 | 0 | *t = ','; |
202 | 0 | s = t + 1; |
203 | 0 | } |
204 | 0 | } while (t && alcnt < 65536); //max allcnt is 65535 |
205 | 0 | line[i] = c; |
206 | 0 | } else if (id == 8) { //INFO, look for "END=" / "SVLEN" |
207 | 0 | int c = line[i], d = 1; |
208 | 0 | line[i] = 0; |
209 | 0 | s = strstr(line + b, "END="); |
210 | 0 | if (s == line + b) s += 4; |
211 | 0 | else if (s) { |
212 | 0 | s = strstr(line + b, ";END="); |
213 | 0 | if (s) s += 5; |
214 | 0 | } |
215 | 0 | if (s && *s != '.') { |
216 | 0 | long long end = strtoll(s, &s, 0); |
217 | 0 | if (end <= intv->beg) { |
218 | 0 | static int reported = 0; |
219 | 0 | if (!reported) { |
220 | 0 | int l = intv->ss ? (int) (intv->se - intv->ss) : 0; |
221 | 0 | hts_log_warning("VCF INFO/END=%lld is smaller than POS at %.*s:%"PRIhts_pos"\n" |
222 | 0 | "This tag will be ignored. " |
223 | 0 | "Note: only one invalid END tag will be reported.", |
224 | 0 | end, l >= 0 ? l : 0, |
225 | 0 | intv->ss ? intv->ss : "", |
226 | 0 | intv->beg); |
227 | 0 | reported = 1; |
228 | 0 | } |
229 | 0 | } else { |
230 | 0 | intv->end = end; |
231 | 0 | } |
232 | 0 | } |
233 | 0 | s = strstr(line + b, "SVLEN="); |
234 | 0 | if (s == line + b) s += 6; //at start of info |
235 | 0 | else if (s) { //not at the start |
236 | 0 | s = strstr(line + b, ";SVLEN="); |
237 | 0 | if (s) s += 7; |
238 | 0 | } |
239 | 0 | while (s && d < alcnt) { |
240 | 0 | t = strchr(s, ','); |
241 | 0 | if ((use_svlen) && (svlenals[d >> 3] & (1 << (d & 7)))) { |
242 | | // <DEL> symbolic allele |
243 | 0 | tmp = atoll(s); |
244 | 0 | tmp = tmp < 0 ? llabs(tmp) : tmp; |
245 | 0 | } else { |
246 | 0 | tmp = 1; |
247 | 0 | } |
248 | 0 | svlen = svlen < tmp ? tmp : svlen; |
249 | 0 | s = t ? t + 1 : NULL; |
250 | 0 | ++d; |
251 | 0 | } |
252 | 0 | line[i] = c; |
253 | 0 | } else if (getlen && id == 9 ) { //FORMAT |
254 | 0 | int c = line[i], pos = -1; |
255 | 0 | line[i] = 0; |
256 | 0 | s = line + b; |
257 | 0 | while (s) { |
258 | 0 | ++pos; |
259 | 0 | if (!(t = strchr(s, ':'))) { //no further fields |
260 | 0 | if (!strcmp(s, "LEN")) { |
261 | 0 | lenpos = pos; |
262 | 0 | } |
263 | 0 | break; //not present at all! |
264 | 0 | } else { |
265 | 0 | *t = '\0'; |
266 | 0 | if (!strcmp(s, "LEN")) { |
267 | 0 | lenpos = pos; |
268 | 0 | *t = ':'; |
269 | 0 | break; |
270 | 0 | } |
271 | 0 | *t = ':'; |
272 | 0 | s = t + 1; //check next one |
273 | 0 | } |
274 | 0 | } |
275 | 0 | line[i] = c; |
276 | 0 | if (lenpos == -1) { //not present |
277 | 0 | break; |
278 | 0 | } |
279 | 0 | } else if (id > 9 && getlen && lenpos != -1) { |
280 | | //get LEN from sample |
281 | 0 | int c = line[i], d = 0; |
282 | 0 | line[i] = 0; tmp = 0; |
283 | 0 | s = line + b; |
284 | 0 | for (d = 0; d <= lenpos; ++d) { |
285 | 0 | if (d == lenpos) { |
286 | 0 | tmp = atoll(s); |
287 | 0 | break; |
288 | 0 | } |
289 | 0 | if ((t = strchr(s, ':'))) { |
290 | 0 | s = t + 1; |
291 | 0 | } else { |
292 | 0 | break; //not in sync with fmt def! |
293 | 0 | } |
294 | 0 | } |
295 | 0 | fmtlen = fmtlen < tmp ? tmp : fmtlen; |
296 | 0 | line[i] = c; |
297 | 0 | } |
298 | 0 | } |
299 | 0 | } |
300 | 0 | b = i + 1; //beginning if current field |
301 | 0 | ++id; |
302 | 0 | } |
303 | 0 | } |
304 | 0 | if ((conf->preset&0xffff) == TBX_VCF) { |
305 | 0 | tmp = reflen < svlen ? |
306 | 0 | svlen < fmtlen ? fmtlen : svlen : |
307 | 0 | reflen < fmtlen ? fmtlen : reflen ; |
308 | 0 | tmp += intv->beg; |
309 | 0 | intv->end = intv->end < tmp ? tmp : intv->end; |
310 | | |
311 | | //NOTE: 'end' calculation be in sync with end/rlen in vcf.c:get_rlen |
312 | 0 | } |
313 | 0 | if (intv->ss == 0 || intv->se == 0 || intv->beg < 0 || intv->end < 0) return -1; |
314 | 0 | return 0; |
315 | 0 | } |
316 | | |
317 | | static inline int get_intv(tbx_t *tbx, kstring_t *str, tbx_intv_t *intv, int is_add) |
318 | 0 | { |
319 | 0 | if (tbx_parse1(&tbx->conf, str->l, str->s, intv) == 0) { |
320 | 0 | int c = *intv->se; |
321 | 0 | *intv->se = '\0'; |
322 | 0 | if ((tbx->conf.preset&0xffff) == TBX_GAF){ |
323 | 0 | intv->tid = 0; |
324 | 0 | } else { |
325 | 0 | intv->tid = get_tid(tbx, intv->ss, is_add); |
326 | 0 | } |
327 | 0 | *intv->se = c; |
328 | 0 | if (intv->tid < 0) return -2; // get_tid out of memory |
329 | 0 | return (intv->beg >= 0 && intv->end >= 0)? 0 : -1; |
330 | 0 | } else { |
331 | 0 | char *type = NULL; |
332 | 0 | switch (tbx->conf.preset&0xffff) |
333 | 0 | { |
334 | 0 | case TBX_SAM: type = "TBX_SAM"; break; |
335 | 0 | case TBX_VCF: type = "TBX_VCF"; break; |
336 | 0 | case TBX_GAF: type = "TBX_GAF"; break; |
337 | 0 | case TBX_UCSC: type = "TBX_UCSC"; break; |
338 | 0 | default: type = "TBX_GENERIC"; break; |
339 | 0 | } |
340 | 0 | if (hts_is_utf16_text(str)) |
341 | 0 | hts_log_error("Failed to parse %s: offending line appears to be encoded as UTF-16", type); |
342 | 0 | else |
343 | 0 | hts_log_error("Failed to parse %s: was wrong -p [type] used?\nThe offending line was: \"%s\"", |
344 | 0 | type, str->s); |
345 | 0 | return -1; |
346 | 0 | } |
347 | 0 | } |
348 | | |
349 | | /* |
350 | | * Called by tabix iterator to read the next record. |
351 | | * Returns >= 0 on success |
352 | | * -1 on EOF |
353 | | * <= -2 on error |
354 | | */ |
355 | | int tbx_readrec(BGZF *fp, void *tbxv, void *sv, int *tid, hts_pos_t *beg, hts_pos_t *end) |
356 | 0 | { |
357 | 0 | tbx_t *tbx = (tbx_t *) tbxv; |
358 | 0 | kstring_t *s = (kstring_t *) sv; |
359 | 0 | int ret; |
360 | | |
361 | | // Get a line until either EOF or a non-meta character |
362 | 0 | do { |
363 | 0 | ret = bgzf_getline(fp, '\n', s); |
364 | 0 | } while (ret >= 0 && s->l && *s->s == tbx->conf.meta_char); |
365 | | |
366 | | // Parse line |
367 | 0 | if (ret >= 0) { |
368 | 0 | tbx_intv_t intv; |
369 | 0 | if (get_intv(tbx, s, &intv, 0) < 0) |
370 | 0 | return -2; |
371 | 0 | *tid = intv.tid; *beg = intv.beg; *end = intv.end; |
372 | 0 | } |
373 | | |
374 | 0 | return ret; |
375 | 0 | } |
376 | | |
377 | | /* |
378 | | Wrapper to get the tbx_t struct to tbx_readrec() when using the |
379 | | multi-region iterator interface. This is required to deal with |
380 | | differences between the single- and multi-region iterator interfaces. |
381 | | |
382 | | In particular, the multi-region one lacks a way to directly pass the tbx_t |
383 | | structure to the tbx_readrec() function. By using the structure below, |
384 | | tbx_itr_next1() can parcel a tbx_t pointer up along with one to the output |
385 | | buffer, then tbx_multi_readrec() can unwrap them to pass on to tbx_readrec(). |
386 | | */ |
387 | | |
388 | | typedef struct tbx_wrapper { |
389 | | void *result_ptr; |
390 | | tbx_t *tbx; |
391 | | } tbx_wrapper; |
392 | | |
393 | | int tbx_itr_next1(htsFile *htsfp, tbx_t *tbx, hts_itr_t *iter, void *r) |
394 | 0 | { |
395 | 0 | if (!htsfp->is_bgzf) { |
396 | 0 | hts_log_error("Only bgzf compressed files can be used with iterators"); |
397 | 0 | errno = EINVAL; |
398 | 0 | return -2; |
399 | 0 | } |
400 | | |
401 | 0 | if (iter->multi) { |
402 | 0 | tbx_wrapper tmp = { r, tbx }; |
403 | 0 | return hts_itr_multi_next(htsfp, iter, &tmp); |
404 | 0 | } else { |
405 | 0 | return hts_itr_next(htsfp->fp.bgzf, iter, r, tbx); |
406 | 0 | } |
407 | 0 | } |
408 | | |
409 | | static int tbx_multi_readrec(BGZF *fp, void *fpv, void *r, |
410 | | int *tid, hts_pos_t *beg, hts_pos_t *end) |
411 | 0 | { |
412 | 0 | tbx_wrapper *tmp = (tbx_wrapper *) r; |
413 | 0 | return tbx_readrec(fp, tmp->tbx, tmp->result_ptr, tid, beg, end); |
414 | 0 | } |
415 | | |
416 | | static int tbx_set_meta(tbx_t *tbx) |
417 | 0 | { |
418 | 0 | int i, l = 0, l_nm; |
419 | 0 | uint32_t x[7]; |
420 | 0 | char **name; |
421 | 0 | uint8_t *meta; |
422 | 0 | khint_t k; |
423 | 0 | khash_t(s2i) *d = (khash_t(s2i)*)tbx->dict; |
424 | |
|
425 | 0 | memcpy(x, &tbx->conf, 24); |
426 | 0 | name = hts_malloc_p(sizeof(char*), kh_size(d)); |
427 | 0 | if (!name) return -1; |
428 | 0 | for (k = kh_begin(d), l = 0; k != kh_end(d); ++k) { |
429 | 0 | if (!kh_exist(d, k)) continue; |
430 | 0 | name[kh_val(d, k)] = (char*)kh_key(d, k); |
431 | 0 | l += strlen(kh_key(d, k)) + 1; // +1 to include '\0' |
432 | 0 | } |
433 | 0 | l_nm = x[6] = l; |
434 | 0 | meta = hts_malloc_ps(sizeof(*meta), l_nm, 28); |
435 | 0 | if (!meta) { free(name); return -1; } |
436 | 0 | if (ed_is_big()) |
437 | 0 | for (i = 0; i < 7; ++i) |
438 | 0 | x[i] = ed_swap_4(x[i]); |
439 | 0 | memcpy(meta, x, 28); |
440 | 0 | for (l = 28, i = 0; i < (int)kh_size(d); ++i) { |
441 | 0 | int x = strlen(name[i]) + 1; |
442 | 0 | memcpy(meta + l, name[i], x); |
443 | 0 | l += x; |
444 | 0 | } |
445 | 0 | free(name); |
446 | 0 | hts_idx_set_meta(tbx->idx, l, meta, 0); |
447 | 0 | return 0; |
448 | 0 | } |
449 | | |
450 | | // Minimal effort parser to extract reference length out of VCF header line |
451 | | // This is used only used to adjust the number of levels if necessary, |
452 | | // so not a major problem if it doesn't always work. |
453 | | static void adjust_max_ref_len_vcf(const char *str, int64_t *max_ref_len) |
454 | 0 | { |
455 | 0 | const char *ptr; |
456 | 0 | int64_t len; |
457 | 0 | if (strncmp(str, "##contig", 8) != 0) return; |
458 | 0 | ptr = strstr(str + 8, "length"); |
459 | 0 | if (!ptr) return; |
460 | 0 | for (ptr += 6; *ptr == ' ' || *ptr == '='; ptr++) {} |
461 | 0 | len = strtoll(ptr, NULL, 10); |
462 | 0 | if (*max_ref_len < len) *max_ref_len = len; |
463 | 0 | } |
464 | | |
465 | | // Same for sam files |
466 | | static void adjust_max_ref_len_sam(const char *str, int64_t *max_ref_len) |
467 | 0 | { |
468 | 0 | const char *ptr; |
469 | 0 | int64_t len; |
470 | 0 | if (strncmp(str, "@SQ", 3) != 0) return; |
471 | 0 | ptr = strstr(str + 3, "\tLN:"); |
472 | 0 | if (!ptr) return; |
473 | 0 | ptr += 4; |
474 | 0 | len = strtoll(ptr, NULL, 10); |
475 | 0 | if (*max_ref_len < len) *max_ref_len = len; |
476 | 0 | } |
477 | | |
478 | | tbx_t *tbx_index(BGZF *fp, int min_shift, const tbx_conf_t *conf) |
479 | 0 | { |
480 | 0 | tbx_t *tbx; |
481 | 0 | kstring_t str; |
482 | 0 | int ret, first = 0, n_lvls, fmt; |
483 | 0 | int64_t lineno = 0; |
484 | 0 | uint64_t last_off = 0; |
485 | 0 | tbx_intv_t intv; |
486 | 0 | int64_t max_ref_len = 0; |
487 | |
|
488 | 0 | str.s = 0; str.l = str.m = 0; |
489 | 0 | tbx = (tbx_t*)calloc(1, sizeof(tbx_t)); |
490 | 0 | if (!tbx) return NULL; |
491 | 0 | tbx->conf = *conf; |
492 | 0 | if (min_shift > 0) n_lvls = (TBX_MAX_SHIFT - min_shift + 2) / 3, fmt = HTS_FMT_CSI; |
493 | 0 | else min_shift = 14, n_lvls = 5, fmt = HTS_FMT_TBI; |
494 | 0 | while ((ret = bgzf_getline(fp, '\n', &str)) >= 0) { |
495 | 0 | ++lineno; |
496 | 0 | if (str.s[0] == tbx->conf.meta_char && fmt == HTS_FMT_CSI) { |
497 | 0 | switch (tbx->conf.preset) { |
498 | 0 | case TBX_SAM: |
499 | 0 | adjust_max_ref_len_sam(str.s, &max_ref_len); break; |
500 | 0 | case TBX_VCF: |
501 | 0 | adjust_max_ref_len_vcf(str.s, &max_ref_len); break; |
502 | 0 | default: |
503 | 0 | break; |
504 | 0 | } |
505 | 0 | } |
506 | 0 | if (lineno <= tbx->conf.line_skip || str.s[0] == tbx->conf.meta_char) { |
507 | 0 | last_off = bgzf_tell(fp); |
508 | 0 | continue; |
509 | 0 | } |
510 | 0 | if (first == 0) { |
511 | 0 | if (fmt == HTS_FMT_CSI) { |
512 | 0 | if (max_ref_len) { |
513 | 0 | hts_adjust_csi_settings(max_ref_len, &min_shift, &n_lvls); |
514 | 0 | } else { |
515 | | // This will give a maximum reference length of at |
516 | | // least 100Gbases for min_shift >= 10, and the |
517 | | // maximum possible for min_shift < 10. |
518 | 0 | const int max_n_lvls = 9; // To prevent bin number overflow |
519 | 0 | n_lvls = (min_shift < 10 |
520 | 0 | ? max_n_lvls |
521 | 0 | : (min_shift < 25 |
522 | 0 | ? max_n_lvls - (min_shift - 10) / 3 |
523 | 0 | : 4)); |
524 | 0 | } |
525 | 0 | } |
526 | 0 | tbx->idx = hts_idx_init(0, fmt, last_off, min_shift, n_lvls); |
527 | 0 | if (!tbx->idx) goto fail; |
528 | 0 | first = 1; |
529 | 0 | } |
530 | 0 | ret = get_intv(tbx, &str, &intv, 1); |
531 | 0 | if (ret < 0) goto fail; // Out of memory or unparsable lines |
532 | 0 | if (hts_idx_push(tbx->idx, intv.tid, intv.beg, intv.end, |
533 | 0 | bgzf_tell(fp), 1) < 0) { |
534 | 0 | goto fail; |
535 | 0 | } |
536 | 0 | } |
537 | 0 | if (ret < -1) goto fail; |
538 | 0 | if ( !tbx->idx ) tbx->idx = hts_idx_init(0, fmt, last_off, min_shift, n_lvls); // empty file |
539 | 0 | if (!tbx->idx) goto fail; |
540 | 0 | if ( !tbx->dict ) tbx->dict = kh_init(s2i); |
541 | 0 | if (!tbx->dict) goto fail; |
542 | 0 | if (hts_idx_finish(tbx->idx, bgzf_tell(fp)) != 0) goto fail; |
543 | 0 | if (tbx_set_meta(tbx) != 0) goto fail; |
544 | 0 | free(str.s); |
545 | 0 | return tbx; |
546 | | |
547 | 0 | fail: |
548 | 0 | free(str.s); |
549 | 0 | tbx_destroy(tbx); |
550 | 0 | return NULL; |
551 | 0 | } |
552 | | |
553 | | void tbx_destroy(tbx_t *tbx) |
554 | 0 | { |
555 | 0 | khash_t(s2i) *d = (khash_t(s2i)*)tbx->dict; |
556 | 0 | if (d != NULL) |
557 | 0 | { |
558 | 0 | khint_t k; |
559 | 0 | for (k = kh_begin(d); k != kh_end(d); ++k) |
560 | 0 | if (kh_exist(d, k)) free((char*)kh_key(d, k)); |
561 | 0 | } |
562 | 0 | hts_idx_destroy(tbx->idx); |
563 | 0 | kh_destroy(s2i, d); |
564 | 0 | free(tbx); |
565 | 0 | } |
566 | | |
567 | | int tbx_index_build3(const char *fn, const char *fnidx, int min_shift, int n_threads, const tbx_conf_t *conf) |
568 | 0 | { |
569 | 0 | tbx_t *tbx; |
570 | 0 | BGZF *fp; |
571 | 0 | int ret; |
572 | 0 | if ((fp = bgzf_open(fn, "r")) == 0) return -1; |
573 | 0 | if ( n_threads ) bgzf_mt(fp, n_threads, 256); |
574 | 0 | if ( bgzf_compression(fp) != bgzf ) { bgzf_close(fp); return -2; } |
575 | 0 | tbx = tbx_index(fp, min_shift, conf); |
576 | 0 | bgzf_close(fp); |
577 | 0 | if ( !tbx ) return -1; |
578 | 0 | ret = hts_idx_save_as(tbx->idx, fn, fnidx, min_shift > 0? HTS_FMT_CSI : HTS_FMT_TBI); |
579 | 0 | tbx_destroy(tbx); |
580 | 0 | return ret; |
581 | 0 | } |
582 | | |
583 | | int tbx_index_build2(const char *fn, const char *fnidx, int min_shift, const tbx_conf_t *conf) |
584 | 0 | { |
585 | 0 | return tbx_index_build3(fn, fnidx, min_shift, 0, conf); |
586 | 0 | } |
587 | | |
588 | | int tbx_index_build(const char *fn, int min_shift, const tbx_conf_t *conf) |
589 | 0 | { |
590 | 0 | return tbx_index_build3(fn, NULL, min_shift, 0, conf); |
591 | 0 | } |
592 | | |
593 | | static tbx_t *index_load(const char *fn, const char *fnidx, int flags) |
594 | 3.96k | { |
595 | 3.96k | tbx_t *tbx; |
596 | 3.96k | uint8_t *meta; |
597 | 3.96k | char *nm, *p; |
598 | 3.96k | uint32_t l_meta, l_nm; |
599 | 3.96k | tbx = (tbx_t*)calloc(1, sizeof(tbx_t)); |
600 | 3.96k | if (!tbx) |
601 | 0 | return NULL; |
602 | 3.96k | tbx->idx = hts_idx_load3(fn, fnidx, HTS_FMT_TBI, flags); |
603 | 3.96k | if ( !tbx->idx ) |
604 | 3.96k | { |
605 | 3.96k | free(tbx); |
606 | 3.96k | return NULL; |
607 | 3.96k | } |
608 | 0 | meta = hts_idx_get_meta(tbx->idx, &l_meta); |
609 | 0 | if ( !meta || l_meta < 28) goto invalid; |
610 | | |
611 | 0 | tbx->conf.preset = le_to_i32(&meta[0]); |
612 | 0 | tbx->conf.sc = le_to_i32(&meta[4]); |
613 | 0 | tbx->conf.bc = le_to_i32(&meta[8]); |
614 | 0 | tbx->conf.ec = le_to_i32(&meta[12]); |
615 | 0 | tbx->conf.meta_char = le_to_i32(&meta[16]); |
616 | 0 | tbx->conf.line_skip = le_to_i32(&meta[20]); |
617 | 0 | l_nm = le_to_u32(&meta[24]); |
618 | 0 | if (l_nm > l_meta - 28) goto invalid; |
619 | | |
620 | 0 | p = nm = (char*)meta + 28; |
621 | | // This assumes meta is NUL-terminated, so we can merrily strlen away. |
622 | | // hts_idx_load_local() assures this for us by adding a NUL on the end |
623 | | // of whatever it reads. |
624 | 0 | for (; p - nm < l_nm; p += strlen(p) + 1) { |
625 | 0 | if (get_tid(tbx, p, 1) < 0) { |
626 | 0 | hts_log_error("%s", strerror(errno)); |
627 | 0 | goto fail; |
628 | 0 | } |
629 | 0 | } |
630 | 0 | return tbx; |
631 | | |
632 | 0 | invalid: |
633 | 0 | hts_log_error("Invalid index header for %s", fnidx ? fnidx : fn); |
634 | |
|
635 | 0 | fail: |
636 | 0 | tbx_destroy(tbx); |
637 | 0 | return NULL; |
638 | 0 | } |
639 | | |
640 | | tbx_t *tbx_index_load3(const char *fn, const char *fnidx, int flags) |
641 | 3.96k | { |
642 | 3.96k | return index_load(fn, fnidx, flags); |
643 | 3.96k | } |
644 | | |
645 | | tbx_t *tbx_index_load2(const char *fn, const char *fnidx) |
646 | 0 | { |
647 | 0 | return index_load(fn, fnidx, 1); |
648 | 0 | } |
649 | | |
650 | | tbx_t *tbx_index_load(const char *fn) |
651 | 0 | { |
652 | 0 | return index_load(fn, NULL, 1); |
653 | 0 | } |
654 | | |
655 | | const char **tbx_seqnames(tbx_t *tbx, int *n) |
656 | 0 | { |
657 | 0 | khash_t(s2i) *d = (khash_t(s2i)*)tbx->dict; |
658 | 0 | if (d == NULL) |
659 | 0 | { |
660 | 0 | *n = 0; |
661 | 0 | return calloc(1, sizeof(char *)); |
662 | 0 | } |
663 | 0 | int tid, m = kh_size(d); |
664 | 0 | const char **names = (const char**) calloc(m,sizeof(const char*)); |
665 | 0 | khint_t k; |
666 | 0 | if (!names) { |
667 | 0 | *n = 0; |
668 | 0 | return NULL; |
669 | 0 | } |
670 | 0 | for (k=kh_begin(d); k<kh_end(d); k++) |
671 | 0 | { |
672 | 0 | if ( !kh_exist(d,k) ) continue; |
673 | 0 | tid = kh_val(d,k); |
674 | 0 | assert( tid<m ); |
675 | 0 | names[tid] = kh_key(d,k); |
676 | 0 | } |
677 | | // sanity check: there should be no gaps |
678 | 0 | for (tid=0; tid<m; tid++) |
679 | 0 | assert(names[tid]); |
680 | 0 | *n = m; |
681 | 0 | return names; |
682 | 0 | } |
683 | | |
684 | | // Wrap around tbx_name2id() to get the right signature for hts_name2id_f |
685 | | static int tbx_name2id_wrapper(void *vhdr, const char *ref) |
686 | 0 | { |
687 | 0 | return tbx_name2id((tbx_t *) vhdr, ref); |
688 | 0 | } |
689 | | |
690 | | hts_itr_t *tbx_itr_querys1(tbx_t *tbx, const char *region) |
691 | 0 | { |
692 | 0 | return hts_itr_querys(tbx->idx, region, tbx_name2id_wrapper, tbx, |
693 | 0 | hts_itr_query, tbx_readrec); |
694 | 0 | } |
695 | | |
696 | | hts_itr_t *tbx_itr_regarray(tbx_t *tbx, char **regarray, unsigned int regcount) |
697 | 0 | { |
698 | 0 | hts_itr_t *itr = NULL; |
699 | 0 | hts_reglist_t *r_list = NULL; |
700 | 0 | int r_count = 0; |
701 | |
|
702 | 0 | r_list = hts_reglist_create(regarray, regcount, &r_count, tbx, |
703 | 0 | tbx_name2id_wrapper); |
704 | 0 | if (!r_list) |
705 | 0 | return NULL; |
706 | | |
707 | 0 | itr = hts_itr_regions(tbx->idx, r_list, r_count, tbx_name2id_wrapper, tbx, |
708 | 0 | hts_itr_multi_bam, tbx_multi_readrec, |
709 | 0 | bgzf_pseek, bgzf_ptell); |
710 | 0 | if (!itr) |
711 | 0 | hts_reglist_free(r_list, r_count); |
712 | |
|
713 | 0 | return itr; |
714 | |
|
715 | 0 | } |
716 | | |
717 | | hts_itr_t *tbx_itr_regions(const tbx_t *tbx, hts_reglist_t *reglist, |
718 | | unsigned int regcount) |
719 | 0 | { |
720 | 0 | return hts_itr_regions(tbx->idx, reglist, regcount, tbx_name2id_wrapper, |
721 | 0 | (void *) tbx, hts_itr_multi_bam, tbx_multi_readrec, |
722 | 0 | bgzf_pseek, bgzf_ptell); |
723 | 0 | } |