1 /*
2 Copyright (c) 2012-2014 Genome Research Ltd.
3 Author: James Bonfield <jkb@sanger.ac.uk>
4
5 Redistribution and use in source and binary forms, with or without
6 modification, are permitted provided that the following conditions are met:
7
8 1. Redistributions of source code must retain the above copyright notice,
9 this list of conditions and the following disclaimer.
10
11 2. Redistributions in binary form must reproduce the above copyright notice,
12 this list of conditions and the following disclaimer in the documentation
13 and/or other materials provided with the distribution.
14
15 3. Neither the names Genome Research Ltd and Wellcome Trust Sanger
16 Institute nor the names of its contributors may be used to endorse or promote
17 products derived from this software without specific prior written permission.
18
19 THIS SOFTWARE IS PROVIDED BY GENOME RESEARCH LTD AND CONTRIBUTORS "AS IS" AND
20 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
21 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
22 DISCLAIMED. IN NO EVENT SHALL GENOME RESEARCH LTD OR CONTRIBUTORS BE LIABLE
23 FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
25 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
26 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
27 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
28 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 */
30
31 /*! \file
32 * Include cram.h instead.
33 *
34 * This is an internal part of the CRAM system and is automatically included
35 * when you #include cram.h.
36 *
37 * Implements the low level CRAM I/O primitives.
38 * This includes basic data types such as byte, int, ITF-8,
39 * maps, bitwise I/O, etc.
40 */
41
42 #ifndef _CRAM_IO_H_
43 #define _CRAM_IO_H_
44
45 #include <stdint.h>
46 #include <cram/misc.h>
47
48 #ifdef __cplusplus
49 extern "C" {
50 #endif
51
52 /**@{ ----------------------------------------------------------------------
53 * ITF8 encoding and decoding.
54 *
55 * Also see the itf8_get and itf8_put macros.
56 */
57
58 /*! INTERNAL: Converts two characters into an integer for use in switch{} */
59 #define CRAM_KEY(a,b) (((a)<<8)|((b)))
60
61 /*! Reads an integer in ITF-8 encoding from 'fd' and stores it in
62 * *val.
63 *
64 * @return
65 * Returns the number of bytes read on success;
66 * -1 on failure
67 */
68 int itf8_decode(cram_fd *fd, int32_t *val);
69
itf8_get(char * cp,int32_t * val_p)70 static inline int itf8_get(char *cp, int32_t *val_p) {
71 unsigned char *up = (unsigned char *)cp;
72
73 if (up[0] < 0x80) {
74 *val_p = up[0];
75 return 1;
76 } else if (up[0] < 0xc0) {
77 *val_p = ((up[0] <<8) | up[1]) & 0x3fff;
78 return 2;
79 } else if (up[0] < 0xe0) {
80 *val_p = ((up[0]<<16) | (up[1]<< 8) | up[2]) & 0x1fffff;
81 return 3;
82 } else if (up[0] < 0xf0) {
83 *val_p = ((up[0]<<24) | (up[1]<<16) | (up[2]<<8) | up[3]) & 0x0fffffff;
84 return 4;
85 } else {
86 *val_p = ((up[0] & 0x0f)<<28) | (up[1]<<20) | (up[2]<<12) | (up[3]<<4) | (up[4] & 0x0f);
87 return 5;
88 }
89 }
90
91 /*
92 * Stores a value to memory in ITF-8 format.
93 *
94 * Returns the number of bytes required to store the number.
95 * This is a maximum of 5 bytes.
96 */
itf8_put(char * cp,int32_t val)97 static inline int itf8_put(char *cp, int32_t val) {
98 unsigned char *up = (unsigned char *)cp;
99 if (!(val & ~0x00000007f)) { // 1 byte
100 *up = val;
101 return 1;
102 } else if (!(val & ~0x00003fff)) { // 2 byte
103 *up++ = (val >> 8 ) | 0x80;
104 *up = val & 0xff;
105 return 2;
106 } else if (!(val & ~0x01fffff)) { // 3 byte
107 *up++ = (val >> 16) | 0xc0;
108 *up++ = (val >> 8 ) & 0xff;
109 *up = val & 0xff;
110 return 3;
111 } else if (!(val & ~0x0fffffff)) { // 4 byte
112 *up++ = (val >> 24) | 0xe0;
113 *up++ = (val >> 16) & 0xff;
114 *up++ = (val >> 8 ) & 0xff;
115 *up = val & 0xff;
116 return 4;
117 } else { // 5 byte
118 *up++ = 0xf0 | ((val>>28) & 0xff);
119 *up++ = (val >> 20) & 0xff;
120 *up++ = (val >> 12) & 0xff;
121 *up++ = (val >> 4 ) & 0xff;
122 *up = val & 0x0f;
123 return 5;
124 }
125 }
126
127
128 /* 64-bit itf8 variant */
ltf8_put(char * cp,int64_t val)129 static inline int ltf8_put(char *cp, int64_t val) {
130 unsigned char *up = (unsigned char *)cp;
131 if (!(val & ~((1LL<<7)-1))) {
132 *up = val;
133 return 1;
134 } else if (!(val & ~((1LL<<(6+8))-1))) {
135 *up++ = (val >> 8 ) | 0x80;
136 *up = val & 0xff;
137 return 2;
138 } else if (!(val & ~((1LL<<(5+2*8))-1))) {
139 *up++ = (val >> 16) | 0xc0;
140 *up++ = (val >> 8 ) & 0xff;
141 *up = val & 0xff;
142 return 3;
143 } else if (!(val & ~((1LL<<(4+3*8))-1))) {
144 *up++ = (val >> 24) | 0xe0;
145 *up++ = (val >> 16) & 0xff;
146 *up++ = (val >> 8 ) & 0xff;
147 *up = val & 0xff;
148 return 4;
149 } else if (!(val & ~((1LL<<(3+4*8))-1))) {
150 *up++ = (val >> 32) | 0xf0;
151 *up++ = (val >> 24) & 0xff;
152 *up++ = (val >> 16) & 0xff;
153 *up++ = (val >> 8 ) & 0xff;
154 *up = val & 0xff;
155 return 5;
156 } else if (!(val & ~((1LL<<(2+5*8))-1))) {
157 *up++ = (val >> 40) | 0xf8;
158 *up++ = (val >> 32) & 0xff;
159 *up++ = (val >> 24) & 0xff;
160 *up++ = (val >> 16) & 0xff;
161 *up++ = (val >> 8 ) & 0xff;
162 *up = val & 0xff;
163 return 6;
164 } else if (!(val & ~((1LL<<(1+6*8))-1))) {
165 *up++ = (val >> 48) | 0xfc;
166 *up++ = (val >> 40) & 0xff;
167 *up++ = (val >> 32) & 0xff;
168 *up++ = (val >> 24) & 0xff;
169 *up++ = (val >> 16) & 0xff;
170 *up++ = (val >> 8 ) & 0xff;
171 *up = val & 0xff;
172 return 7;
173 } else if (!(val & ~((1LL<<(7*8))-1))) {
174 *up++ = (val >> 56) | 0xfe;
175 *up++ = (val >> 48) & 0xff;
176 *up++ = (val >> 40) & 0xff;
177 *up++ = (val >> 32) & 0xff;
178 *up++ = (val >> 24) & 0xff;
179 *up++ = (val >> 16) & 0xff;
180 *up++ = (val >> 8 ) & 0xff;
181 *up = val & 0xff;
182 return 8;
183 } else {
184 *up++ = 0xff;
185 *up++ = (val >> 56) & 0xff;
186 *up++ = (val >> 48) & 0xff;
187 *up++ = (val >> 40) & 0xff;
188 *up++ = (val >> 32) & 0xff;
189 *up++ = (val >> 24) & 0xff;
190 *up++ = (val >> 16) & 0xff;
191 *up++ = (val >> 8 ) & 0xff;
192 *up = val & 0xff;
193 return 9;
194 }
195 }
196
ltf8_get(char * cp,int64_t * val_p)197 static inline int ltf8_get(char *cp, int64_t *val_p) {
198 unsigned char *up = (unsigned char *)cp;
199
200 if (up[0] < 0x80) {
201 *val_p = up[0];
202 return 1;
203 } else if (up[0] < 0xc0) {
204 *val_p = (((uint64_t)up[0]<< 8) |
205 (uint64_t)up[1]) & (((1LL<<(6+8)))-1);
206 return 2;
207 } else if (up[0] < 0xe0) {
208 *val_p = (((uint64_t)up[0]<<16) |
209 ((uint64_t)up[1]<< 8) |
210 (uint64_t)up[2]) & ((1LL<<(5+2*8))-1);
211 return 3;
212 } else if (up[0] < 0xf0) {
213 *val_p = (((uint64_t)up[0]<<24) |
214 ((uint64_t)up[1]<<16) |
215 ((uint64_t)up[2]<< 8) |
216 (uint64_t)up[3]) & ((1LL<<(4+3*8))-1);
217 return 4;
218 } else if (up[0] < 0xf8) {
219 *val_p = (((uint64_t)up[0]<<32) |
220 ((uint64_t)up[1]<<24) |
221 ((uint64_t)up[2]<<16) |
222 ((uint64_t)up[3]<< 8) |
223 (uint64_t)up[4]) & ((1LL<<(3+4*8))-1);
224 return 5;
225 } else if (up[0] < 0xfc) {
226 *val_p = (((uint64_t)up[0]<<40) |
227 ((uint64_t)up[1]<<32) |
228 ((uint64_t)up[2]<<24) |
229 ((uint64_t)up[3]<<16) |
230 ((uint64_t)up[4]<< 8) |
231 (uint64_t)up[5]) & ((1LL<<(2+5*8))-1);
232 return 6;
233 } else if (up[0] < 0xfe) {
234 *val_p = (((uint64_t)up[0]<<48) |
235 ((uint64_t)up[1]<<40) |
236 ((uint64_t)up[2]<<32) |
237 ((uint64_t)up[3]<<24) |
238 ((uint64_t)up[4]<<16) |
239 ((uint64_t)up[5]<< 8) |
240 (uint64_t)up[6]) & ((1LL<<(1+6*8))-1);
241 return 7;
242 } else if (up[0] < 0xff) {
243 *val_p = (((uint64_t)up[1]<<48) |
244 ((uint64_t)up[2]<<40) |
245 ((uint64_t)up[3]<<32) |
246 ((uint64_t)up[4]<<24) |
247 ((uint64_t)up[5]<<16) |
248 ((uint64_t)up[6]<< 8) |
249 (uint64_t)up[7]) & ((1LL<<(7*8))-1);
250 return 8;
251 } else {
252 *val_p = (((uint64_t)up[1]<<56) |
253 ((uint64_t)up[2]<<48) |
254 ((uint64_t)up[3]<<40) |
255 ((uint64_t)up[4]<<32) |
256 ((uint64_t)up[5]<<24) |
257 ((uint64_t)up[6]<<16) |
258 ((uint64_t)up[7]<< 8) |
259 (uint64_t)up[8]);
260 return 9;
261 }
262 }
263
264 #define itf8_size(v) ((!((v)&~0x7f))?1:(!((v)&~0x3fff))?2:(!((v)&~0x1fffff))?3:(!((v)&~0xfffffff))?4:5)
265
266
267 /* Version of itf8_get that checks it hasn't run out of input */
268
269 extern const int itf8_bytes[16];
270 extern const int ltf8_bytes[256];
271
safe_itf8_get(const char * cp,const char * endp,int32_t * val_p)272 static inline int safe_itf8_get(const char *cp, const char *endp,
273 int32_t *val_p) {
274 const unsigned char *up = (unsigned char *)cp;
275
276 if (endp - cp < 5 &&
277 (cp >= endp || endp - cp < itf8_bytes[up[0]>>4])) {
278 *val_p = 0;
279 return 0;
280 }
281
282 if (up[0] < 0x80) {
283 *val_p = up[0];
284 return 1;
285 } else if (up[0] < 0xc0) {
286 *val_p = ((up[0] <<8) | up[1]) & 0x3fff;
287 return 2;
288 } else if (up[0] < 0xe0) {
289 *val_p = ((up[0]<<16) | (up[1]<< 8) | up[2]) & 0x1fffff;
290 return 3;
291 } else if (up[0] < 0xf0) {
292 *val_p = (((uint32_t)up[0]<<24) | (up[1]<<16) | (up[2]<<8) | up[3]) & 0x0fffffff;
293 return 4;
294 } else {
295 uint32_t uv = (((uint32_t)up[0] & 0x0f)<<28) | (up[1]<<20) | (up[2]<<12) | (up[3]<<4) | (up[4] & 0x0f);
296 *val_p = uv < 0x80000000UL ? uv : -((int32_t) (0xffffffffUL - uv)) - 1;
297 return 5;
298 }
299 }
300
safe_ltf8_get(const char * cp,const char * endp,int64_t * val_p)301 static inline int safe_ltf8_get(const char *cp, const char *endp,
302 int64_t *val_p) {
303 unsigned char *up = (unsigned char *)cp;
304
305 if (endp - cp < 9 &&
306 (cp >= endp || endp - cp < ltf8_bytes[up[0]])) return 0;
307
308 if (up[0] < 0x80) {
309 *val_p = up[0];
310 return 1;
311 } else if (up[0] < 0xc0) {
312 *val_p = (((uint64_t)up[0]<< 8) |
313 (uint64_t)up[1]) & (((1LL<<(6+8)))-1);
314 return 2;
315 } else if (up[0] < 0xe0) {
316 *val_p = (((uint64_t)up[0]<<16) |
317 ((uint64_t)up[1]<< 8) |
318 (uint64_t)up[2]) & ((1LL<<(5+2*8))-1);
319 return 3;
320 } else if (up[0] < 0xf0) {
321 *val_p = (((uint64_t)up[0]<<24) |
322 ((uint64_t)up[1]<<16) |
323 ((uint64_t)up[2]<< 8) |
324 (uint64_t)up[3]) & ((1LL<<(4+3*8))-1);
325 return 4;
326 } else if (up[0] < 0xf8) {
327 *val_p = (((uint64_t)up[0]<<32) |
328 ((uint64_t)up[1]<<24) |
329 ((uint64_t)up[2]<<16) |
330 ((uint64_t)up[3]<< 8) |
331 (uint64_t)up[4]) & ((1LL<<(3+4*8))-1);
332 return 5;
333 } else if (up[0] < 0xfc) {
334 *val_p = (((uint64_t)up[0]<<40) |
335 ((uint64_t)up[1]<<32) |
336 ((uint64_t)up[2]<<24) |
337 ((uint64_t)up[3]<<16) |
338 ((uint64_t)up[4]<< 8) |
339 (uint64_t)up[5]) & ((1LL<<(2+5*8))-1);
340 return 6;
341 } else if (up[0] < 0xfe) {
342 *val_p = (((uint64_t)up[0]<<48) |
343 ((uint64_t)up[1]<<40) |
344 ((uint64_t)up[2]<<32) |
345 ((uint64_t)up[3]<<24) |
346 ((uint64_t)up[4]<<16) |
347 ((uint64_t)up[5]<< 8) |
348 (uint64_t)up[6]) & ((1LL<<(1+6*8))-1);
349 return 7;
350 } else if (up[0] < 0xff) {
351 *val_p = (((uint64_t)up[1]<<48) |
352 ((uint64_t)up[2]<<40) |
353 ((uint64_t)up[3]<<32) |
354 ((uint64_t)up[4]<<24) |
355 ((uint64_t)up[5]<<16) |
356 ((uint64_t)up[6]<< 8) |
357 (uint64_t)up[7]) & ((1LL<<(7*8))-1);
358 return 8;
359 } else {
360 *val_p = (((uint64_t)up[1]<<56) |
361 ((uint64_t)up[2]<<48) |
362 ((uint64_t)up[3]<<40) |
363 ((uint64_t)up[4]<<32) |
364 ((uint64_t)up[5]<<24) |
365 ((uint64_t)up[6]<<16) |
366 ((uint64_t)up[7]<< 8) |
367 (uint64_t)up[8]);
368 return 9;
369 }
370 }
371
372 /*! Pushes a value in ITF8 format onto the end of a block.
373 *
374 * This shouldn't be used for high-volume data as it is not the fastest
375 * method.
376 *
377 * @return
378 * Returns the number of bytes written
379 */
380 int itf8_put_blk(cram_block *blk, int val);
381
382 /*! Pulls a literal 32-bit value from a block.
383 *
384 * @returns the number of bytes decoded;
385 * -1 on failure.
386 */
387 int int32_get_blk(cram_block *b, int32_t *val);
388
389 /*! Pushes a literal 32-bit value onto the end of a block.
390 *
391 * @return
392 * Returns 0 on success;
393 * -1 on failure.
394 */
395 int int32_put_blk(cram_block *blk, int32_t val);
396
397
398 /**@}*/
399 /**@{ ----------------------------------------------------------------------
400 * CRAM blocks - the dynamically growable data block. We have code to
401 * create, update, (un)compress and read/write.
402 *
403 * These are derived from the deflate_interlaced.c blocks, but with the
404 * CRAM extension of content types and IDs.
405 */
406
407 /*! Allocates a new cram_block structure with a specified content_type and
408 * id.
409 *
410 * @return
411 * Returns block pointer on success;
412 * NULL on failure
413 */
414 cram_block *cram_new_block(enum cram_content_type content_type,
415 int content_id);
416
417 /*! Reads a block from a cram file.
418 *
419 * @return
420 * Returns cram_block pointer on success;
421 * NULL on failure
422 */
423 cram_block *cram_read_block(cram_fd *fd);
424
425 /*! Writes a CRAM block.
426 *
427 * @return
428 * Returns 0 on success;
429 * -1 on failure
430 */
431 int cram_write_block(cram_fd *fd, cram_block *b);
432
433 /*! Frees a CRAM block, deallocating internal data too.
434 */
435 void cram_free_block(cram_block *b);
436
437 /*! Uncompress a memory block using Zlib.
438 *
439 * @return
440 * Returns 0 on success;
441 * -1 on failure
442 */
443 char *zlib_mem_inflate(char *cdata, size_t csize, size_t *size);
444
445 /*! Uncompresses a CRAM block, if compressed.
446 *
447 * @return
448 * Returns 0 on success;
449 * -1 on failure
450 */
451 int cram_uncompress_block(cram_block *b);
452
453 /*! Compresses a block.
454 *
455 * Compresses a block using one of two different zlib strategies. If we only
456 * want one choice set strat2 to be -1.
457 *
458 * The logic here is that sometimes Z_RLE does a better job than Z_FILTERED
459 * or Z_DEFAULT_STRATEGY on quality data. If so, we'd rather use it as it is
460 * significantly faster.
461 *
462 * @return
463 * Returns 0 on success;
464 * -1 on failure
465 */
466 int cram_compress_block(cram_fd *fd, cram_block *b, cram_metrics *metrics,
467 int method, int level);
468
469 cram_metrics *cram_new_metrics(void);
470 char *cram_block_method2str(enum cram_block_method m);
471 char *cram_content_type2str(enum cram_content_type t);
472
473 /*
474 * Find an external block by its content_id
475 */
476
cram_get_block_by_id(cram_slice * slice,int id)477 static inline cram_block *cram_get_block_by_id(cram_slice *slice, int id) {
478 if (slice->block_by_id && id >= 0 && id < 1024) {
479 return slice->block_by_id[id];
480 } else {
481 int i;
482 for (i = 0; i < slice->hdr->num_blocks; i++) {
483 cram_block *b = slice->block[i];
484 if (b && b->content_type == EXTERNAL && b->content_id == id)
485 return b;
486 }
487 }
488 return NULL;
489 }
490
491 /* --- Accessor macros for manipulating blocks on a byte by byte basis --- */
492
493 /* Block size and data pointer. */
494 #define BLOCK_SIZE(b) ((b)->byte)
495 #define BLOCK_DATA(b) ((b)->data)
496
497 /* Returns the address one past the end of the block */
498 #define BLOCK_END(b) (&(b)->data[(b)->byte])
499
500 /* Request block to be at least 'l' bytes long */
501 #define BLOCK_RESIZE(b,l) \
502 do { \
503 while((b)->alloc <= (l)) { \
504 (b)->alloc = (b)->alloc ? (b)->alloc*1.5 : 1024; \
505 (b)->data = realloc((b)->data, (b)->alloc); \
506 } \
507 } while(0)
508
509 /* Make block exactly 'l' bytes long */
510 #define BLOCK_RESIZE_EXACT(b,l) \
511 do { \
512 (b)->alloc = (l); \
513 (b)->data = realloc((b)->data, (b)->alloc); \
514 } while(0)
515
516 /* Ensure the block can hold at least another 'l' bytes */
517 #define BLOCK_GROW(b,l) BLOCK_RESIZE((b), BLOCK_SIZE((b)) + (l))
518
519 /* Append string 's' of length 'l' */
520 #define BLOCK_APPEND(b,s,l) \
521 do { \
522 BLOCK_GROW((b),(l)); \
523 memcpy(BLOCK_END((b)), (s), (l)); \
524 BLOCK_SIZE((b)) += (l); \
525 } while (0)
526
527 /* Append as single character 'c' */
528 #define BLOCK_APPEND_CHAR(b,c) \
529 do { \
530 BLOCK_GROW((b),1); \
531 (b)->data[(b)->byte++] = (c); \
532 } while (0)
533
534 /* Append a single unsigned integer */
535 #define BLOCK_APPEND_UINT(b,i) \
536 do { \
537 unsigned char *cp; \
538 BLOCK_GROW((b),11); \
539 cp = &(b)->data[(b)->byte]; \
540 (b)->byte += append_uint32(cp, (i)) - cp; \
541 } while (0)
542
append_uint32(unsigned char * cp,uint32_t i)543 static inline unsigned char *append_uint32(unsigned char *cp, uint32_t i) {
544 uint32_t j;
545
546 if (i == 0) {
547 *cp++ = '0';
548 return cp;
549 }
550
551 if (i < 100) goto b1;
552 if (i < 10000) goto b3;
553 if (i < 1000000) goto b5;
554 if (i < 100000000) goto b7;
555
556 if ((j = i / 1000000000)) {*cp++ = j + '0'; i -= j*1000000000; goto x8;}
557 if ((j = i / 100000000)) {*cp++ = j + '0'; i -= j*100000000; goto x7;}
558 b7:if ((j = i / 10000000)) {*cp++ = j + '0'; i -= j*10000000; goto x6;}
559 if ((j = i / 1000000)) {*cp++ = j + '0', i -= j*1000000; goto x5;}
560 b5:if ((j = i / 100000)) {*cp++ = j + '0', i -= j*100000; goto x4;}
561 if ((j = i / 10000)) {*cp++ = j + '0', i -= j*10000; goto x3;}
562 b3:if ((j = i / 1000)) {*cp++ = j + '0', i -= j*1000; goto x2;}
563 if ((j = i / 100)) {*cp++ = j + '0', i -= j*100; goto x1;}
564 b1:if ((j = i / 10)) {*cp++ = j + '0', i -= j*10; goto x0;}
565 if (i) *cp++ = i + '0';
566 return cp;
567
568 x8: *cp++ = i / 100000000 + '0', i %= 100000000;
569 x7: *cp++ = i / 10000000 + '0', i %= 10000000;
570 x6: *cp++ = i / 1000000 + '0', i %= 1000000;
571 x5: *cp++ = i / 100000 + '0', i %= 100000;
572 x4: *cp++ = i / 10000 + '0', i %= 10000;
573 x3: *cp++ = i / 1000 + '0', i %= 1000;
574 x2: *cp++ = i / 100 + '0', i %= 100;
575 x1: *cp++ = i / 10 + '0', i %= 10;
576 x0: *cp++ = i + '0';
577
578 return cp;
579 }
580
append_sub32(unsigned char * cp,uint32_t i)581 static inline unsigned char *append_sub32(unsigned char *cp, uint32_t i) {
582 *cp++ = i / 100000000 + '0', i %= 100000000;
583 *cp++ = i / 10000000 + '0', i %= 10000000;
584 *cp++ = i / 1000000 + '0', i %= 1000000;
585 *cp++ = i / 100000 + '0', i %= 100000;
586 *cp++ = i / 10000 + '0', i %= 10000;
587 *cp++ = i / 1000 + '0', i %= 1000;
588 *cp++ = i / 100 + '0', i %= 100;
589 *cp++ = i / 10 + '0', i %= 10;
590 *cp++ = i + '0';
591
592 return cp;
593 }
594
append_uint64(unsigned char * cp,uint64_t i)595 static inline unsigned char *append_uint64(unsigned char *cp, uint64_t i) {
596 uint64_t j;
597
598 if (i <= 0xffffffff)
599 return append_uint32(cp, i);
600
601 if ((j = i/1000000000) > 1000000000) {
602 cp = append_uint32(cp, j/1000000000);
603 j %= 1000000000;
604 cp = append_sub32(cp, j);
605 } else {
606 cp = append_uint32(cp, i / 1000000000);
607 }
608 cp = append_sub32(cp, i % 1000000000);
609
610 return cp;
611 }
612
613 #define BLOCK_UPLEN(b) \
614 (b)->comp_size = (b)->uncomp_size = BLOCK_SIZE((b))
615
616 /**@}*/
617 /**@{ ----------------------------------------------------------------------
618 * Reference sequence handling
619 */
620
621 /*! Loads a reference set from fn and stores in the cram_fd.
622 *
623 * @return
624 * Returns 0 on success;
625 * -1 on failure
626 */
627 int cram_load_reference(cram_fd *fd, char *fn);
628
629 /*! Generates a lookup table in refs based on the SQ headers in SAM_hdr.
630 *
631 * Indexes references by the order they appear in a BAM file. This may not
632 * necessarily be the same order they appear in the fasta reference file.
633 *
634 * @return
635 * Returns 0 on success;
636 * -1 on failure
637 */
638 int refs2id(refs_t *r, SAM_hdr *bfd);
639
640 void refs_free(refs_t *r);
641
642 /*! Returns a portion of a reference sequence from start to end inclusive.
643 *
644 * The returned pointer is owned by the cram_file fd and should not be freed
645 * by the caller. It is valid only until the next cram_get_ref is called
646 * with the same fd parameter (so is thread-safe if given multiple files).
647 *
648 * To return the entire reference sequence, specify start as 1 and end
649 * as 0.
650 *
651 * @return
652 * Returns reference on success;
653 * NULL on failure
654 */
655 char *cram_get_ref(cram_fd *fd, int id, int start, int end);
656 void cram_ref_incr(refs_t *r, int id);
657 void cram_ref_decr(refs_t *r, int id);
658 /**@}*/
659 /**@{ ----------------------------------------------------------------------
660 * Containers
661 */
662
663 /*! Creates a new container, specifying the maximum number of slices
664 * and records permitted.
665 *
666 * @return
667 * Returns cram_container ptr on success;
668 * NULL on failure
669 */
670 cram_container *cram_new_container(int nrec, int nslice);
671 void cram_free_container(cram_container *c);
672
673 /*! Reads a container header.
674 *
675 * @return
676 * Returns cram_container on success;
677 * NULL on failure or no container left (fd->err == 0).
678 */
679 cram_container *cram_read_container(cram_fd *fd);
680
681 /*! Writes a container structure.
682 *
683 * @return
684 * Returns 0 on success;
685 * -1 on failure
686 */
687 int cram_write_container(cram_fd *fd, cram_container *h);
688
689 /*! Flushes a container to disk.
690 *
691 * Flushes a completely or partially full container to disk, writing
692 * container structure, header and blocks. This also calls the encoder
693 * functions.
694 *
695 * @return
696 * Returns 0 on success;
697 * -1 on failure
698 */
699 int cram_flush_container(cram_fd *fd, cram_container *c);
700 int cram_flush_container_mt(cram_fd *fd, cram_container *c);
701
702
703 /**@}*/
704 /**@{ ----------------------------------------------------------------------
705 * Compression headers; the first part of the container
706 */
707
708 /*! Creates a new blank container compression header
709 *
710 * @return
711 * Returns header ptr on success;
712 * NULL on failure
713 */
714 cram_block_compression_hdr *cram_new_compression_header(void);
715
716 /*! Frees a cram_block_compression_hdr */
717 void cram_free_compression_header(cram_block_compression_hdr *hdr);
718
719
720 /**@}*/
721 /**@{ ----------------------------------------------------------------------
722 * Slices and slice headers
723 */
724
725 /*! Frees a slice header */
726 void cram_free_slice_header(cram_block_slice_hdr *hdr);
727
728 /*! Frees a slice */
729 void cram_free_slice(cram_slice *s);
730
731 /*! Creates a new empty slice in memory, for subsequent writing to
732 * disk.
733 *
734 * @return
735 * Returns cram_slice ptr on success;
736 * NULL on failure
737 */
738 cram_slice *cram_new_slice(enum cram_content_type type, int nrecs);
739
740 /*! Loads an entire slice.
741 *
742 * FIXME: In 1.0 the native unit of slices within CRAM is broken
743 * as slices contain references to objects in other slices.
744 * To work around this while keeping the slice oriented outer loop
745 * we read all slices and stitch them together into a fake large
746 * slice instead.
747 *
748 * @return
749 * Returns cram_slice ptr on success;
750 * NULL on failure
751 */
752 cram_slice *cram_read_slice(cram_fd *fd);
753
754
755
756 /**@}*/
757 /**@{ ----------------------------------------------------------------------
758 * CRAM file definition (header)
759 */
760
761 /*! Reads a CRAM file definition structure.
762 *
763 * @return
764 * Returns file_def ptr on success;
765 * NULL on failure
766 */
767 cram_file_def *cram_read_file_def(cram_fd *fd);
768
769 /*! Writes a cram_file_def structure to cram_fd.
770 *
771 * @return
772 * Returns 0 on success;
773 * -1 on failure
774 */
775 int cram_write_file_def(cram_fd *fd, cram_file_def *def);
776
777 /*! Frees a cram_file_def structure. */
778 void cram_free_file_def(cram_file_def *def);
779
780
781 /**@}*/
782 /**@{ ----------------------------------------------------------------------
783 * SAM header I/O
784 */
785
786 /*! Reads the SAM header from the first CRAM data block.
787 *
788 * Also performs minimal parsing to extract read-group
789 * and sample information.
790 *
791 * @return
792 * Returns SAM hdr ptr on success;
793 * NULL on failure
794 */
795 SAM_hdr *cram_read_SAM_hdr(cram_fd *fd);
796
797 /*! Writes a CRAM SAM header.
798 *
799 * @return
800 * Returns 0 on success;
801 * -1 on failure
802 */
803 int cram_write_SAM_hdr(cram_fd *fd, SAM_hdr *hdr);
804
805
806 /**@}*/
807 /**@{ ----------------------------------------------------------------------
808 * The top-level cram opening, closing and option handling
809 */
810
811 /*! Opens a CRAM file for read (mode "rb") or write ("wb").
812 *
813 * The filename may be "-" to indicate stdin or stdout.
814 *
815 * @return
816 * Returns file handle on success;
817 * NULL on failure.
818 */
819 cram_fd *cram_open(const char *filename, const char *mode);
820
821 /*! Opens an existing stream for reading or writing.
822 *
823 * @return
824 * Returns file handle on success;
825 * NULL on failure.
826 */
827 cram_fd *cram_dopen(struct hFILE *fp, const char *filename, const char *mode);
828
829 /*! Closes a CRAM file.
830 *
831 * @return
832 * Returns 0 on success;
833 * -1 on failure
834 */
835 int cram_close(cram_fd *fd);
836
837 /*
838 * Seek within a CRAM file.
839 *
840 * Returns 0 on success
841 * -1 on failure
842 */
843 int cram_seek(cram_fd *fd, off_t offset, int whence);
844
845 int64_t cram_tell(cram_fd *fd);
846
847 /*
848 * Flushes a CRAM file.
849 * Useful for when writing to stdout without wishing to close the stream.
850 *
851 * Returns 0 on success
852 * -1 on failure
853 */
854 int cram_flush(cram_fd *fd);
855
856 /*! Checks for end of file on a cram_fd stream.
857 *
858 * @return
859 * Returns 0 if not at end of file
860 * 1 if we hit an expected EOF (end of range or EOF block)
861 * 2 for other EOF (end of stream without EOF block)
862 */
863 int cram_eof(cram_fd *fd);
864
865 /*! Sets options on the cram_fd.
866 *
867 * See CRAM_OPT_* definitions in cram_structs.h.
868 * Use this immediately after opening.
869 *
870 * @return
871 * Returns 0 on success;
872 * -1 on failure
873 */
874 int cram_set_option(cram_fd *fd, enum hts_fmt_option opt, ...);
875
876 /*! Sets options on the cram_fd.
877 *
878 * See CRAM_OPT_* definitions in cram_structs.h.
879 * Use this immediately after opening.
880 *
881 * @return
882 * Returns 0 on success;
883 * -1 on failure
884 */
885 int cram_set_voption(cram_fd *fd, enum hts_fmt_option opt, va_list args);
886
887 /*!
888 * Attaches a header to a cram_fd.
889 *
890 * This should be used when creating a new cram_fd for writing where
891 * we have an SAM_hdr already constructed (eg from a file we've read
892 * in).
893 *
894 * @return
895 * Returns 0 on success;
896 * -1 on failure
897 */
898 int cram_set_header(cram_fd *fd, SAM_hdr *hdr);
899
900 /*!
901 * Returns the hFILE connected to a cram_fd.
902 */
cram_hfile(cram_fd * fd)903 static inline struct hFILE *cram_hfile(cram_fd *fd) {
904 return fd->fp;
905 }
906
907 #ifdef __cplusplus
908 }
909 #endif
910
911 #endif /* _CRAM_IO_H_ */
912