Mercurial > repos > youngkim > ezbamqc

comparison ezBAMQC/src/htslib/cram/cram_decode.c @ 0:dfa3745e5fd8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author youngkim
date Thu, 24 Mar 2016 17:12:52 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:dfa3745e5fd8
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 /*
32 * - In-memory decoding of CRAM data structures.
33 * - Iterator for reading CRAM record by record.
34 */
35
36 #ifdef HAVE_CONFIG_H
37 #include "io_lib_config.h"
38 #endif
39
40 #include <stdio.h>
41 #include <errno.h>
42 #include <assert.h>
43 #include <stdlib.h>
44 #include <string.h>
45 #include <zlib.h>
46 #include <sys/types.h>
47 #include <sys/stat.h>
48 #include <math.h>
49 #include <ctype.h>
50
51 #include "cram/cram.h"
52 #include "cram/os.h"
53 #include "cram/md5.h"
54
55 //Whether CIGAR has just M or uses = and X to indicate match and mismatch
56 //#define USE_X
57
58 /* ----------------------------------------------------------------------
59 * CRAM compression headers
60 */
61
62 /*
63 * Decodes the Tag Dictionary record in the preservation map
64 * Updates the cram compression header.
65 *
66 * Returns number of bytes decoded on success
67 * -1 on failure
68 */
69 int cram_decode_TD(char *cp, cram_block_compression_hdr *h) {
70 char *op = cp;
71 unsigned char *dat;
72 cram_block *b;
73 int32_t blk_size;
74 int nTL, i, sz;
75
76 if (!(b = cram_new_block(0, 0)))
77 return -1;
78 h->TD_blk = b;
79
80 /* Decode */
81 cp += itf8_get(cp, &blk_size);
82 if (!blk_size) {
83 h->nTL = 0;
84 h->TL = NULL;
85 cram_free_block(b);
86 return cp - op;
87 }
88
89 BLOCK_APPEND(b, cp, blk_size);
90 cp += blk_size;
91 sz = cp - op;
92
93 // Force nul termination if missing
94 if (BLOCK_DATA(b)[BLOCK_SIZE(b)-1])
95 BLOCK_APPEND_CHAR(b, '\0');
96
97 /* Set up TL lookup table */
98 dat = BLOCK_DATA(b);
99
100 // Count
101 for (nTL = i = 0; i < BLOCK_SIZE(b); i++) {
102 nTL++;
103 while (dat[i])
104 i++;
105 }
106
107 // Copy
108 h->nTL = nTL;
109 if (!(h->TL = calloc(h->nTL, sizeof(unsigned char *))))
110 return -1;
111 for (nTL = i = 0; i < BLOCK_SIZE(b); i++) {
112 h->TL[nTL++] = &dat[i];
113 while (dat[i])
114 i++;
115 }
116
117 return sz;
118 }
119
120 /*
121 * Decodes a CRAM block compression header.
122 * Returns header ptr on success
123 * NULL on failure
124 */
125 cram_block_compression_hdr *cram_decode_compression_header(cram_fd *fd,
126 cram_block *b) {
127 char *cp, *cp_copy;
128 cram_block_compression_hdr *hdr = calloc(1, sizeof(*hdr));
129 int i;
130 int32_t map_size, map_count;
131
132 if (!hdr)
133 return NULL;
134
135 if (b->method != RAW) {
136 if (cram_uncompress_block(b)) {
137 free(hdr);
138 return NULL;
139 }
140 }
141
142 cp = (char *)b->data;
143
144 if (CRAM_MAJOR_VERS(fd->version) == 1) {
145 cp += itf8_get(cp, &hdr->ref_seq_id);
146 cp += itf8_get(cp, &hdr->ref_seq_start);
147 cp += itf8_get(cp, &hdr->ref_seq_span);
148 cp += itf8_get(cp, &hdr->num_records);
149 cp += itf8_get(cp, &hdr->num_landmarks);
150 if (!(hdr->landmark = malloc(hdr->num_landmarks * sizeof(int32_t)))) {
151 free(hdr);
152 return NULL;
153 }
154 for (i = 0; i < hdr->num_landmarks; i++) {
155 cp += itf8_get(cp, &hdr->landmark[i]);
156 }
157 }
158
159 hdr->preservation_map = kh_init(map);
160
161 memset(hdr->rec_encoding_map, 0,
162 CRAM_MAP_HASH * sizeof(hdr->rec_encoding_map[0]));
163 memset(hdr->tag_encoding_map, 0,
164 CRAM_MAP_HASH * sizeof(hdr->tag_encoding_map[0]));
165
166 if (!hdr->preservation_map) {
167 cram_free_compression_header(hdr);
168 return NULL;
169 }
170
171 /* Initialise defaults for preservation map */
172 hdr->mapped_qs_included = 0;
173 hdr->unmapped_qs_included = 0;
174 hdr->unmapped_placed = 0;
175 hdr->qs_included = 0;
176 hdr->read_names_included = 0;
177 hdr->AP_delta = 1;
178 memcpy(hdr->substitution_matrix, "CGTNAGTNACTNACGNACGT", 20);
179
180 /* Preservation map */
181 cp += itf8_get(cp, &map_size); cp_copy = cp;
182 cp += itf8_get(cp, &map_count);
183 for (i = 0; i < map_count; i++) {
184 pmap_t hd;
185 khint_t k;
186 int r;
187
188 cp += 2;
189 switch(CRAM_KEY(cp[-2],cp[-1])) {
190 case CRAM_KEY('M','I'):
191 hd.i = *cp++;
192 k = kh_put(map, hdr->preservation_map, "MI", &r);
193 if (-1 == r) {
194 cram_free_compression_header(hdr);
195 return NULL;
196 }
197
198 kh_val(hdr->preservation_map, k) = hd;
199 hdr->mapped_qs_included = hd.i;
200 break;
201
202 case CRAM_KEY('U','I'):
203 hd.i = *cp++;
204 k = kh_put(map, hdr->preservation_map, "UI", &r);
205 if (-1 == r) {
206 cram_free_compression_header(hdr);
207 return NULL;
208 }
209
210 kh_val(hdr->preservation_map, k) = hd;
211 hdr->unmapped_qs_included = hd.i;
212 break;
213
214 case CRAM_KEY('P','I'):
215 hd.i = *cp++;
216 k = kh_put(map, hdr->preservation_map, "PI", &r);
217 if (-1 == r) {
218 cram_free_compression_header(hdr);
219 return NULL;
220 }
221
222 kh_val(hdr->preservation_map, k) = hd;
223 hdr->unmapped_placed = hd.i;
224 break;
225
226 case CRAM_KEY('R','N'):
227 hd.i = *cp++;
228 k = kh_put(map, hdr->preservation_map, "RN", &r);
229 if (-1 == r) {
230 cram_free_compression_header(hdr);
231 return NULL;
232 }
233
234 kh_val(hdr->preservation_map, k) = hd;
235 hdr->read_names_included = hd.i;
236 break;
237
238 case CRAM_KEY('A','P'):
239 hd.i = *cp++;
240 k = kh_put(map, hdr->preservation_map, "AP", &r);
241 if (-1 == r) {
242 cram_free_compression_header(hdr);
243 return NULL;
244 }
245
246 kh_val(hdr->preservation_map, k) = hd;
247 hdr->AP_delta = hd.i;
248 break;
249
250 case CRAM_KEY('R','R'):
251 hd.i = *cp++;
252 k = kh_put(map, hdr->preservation_map, "RR", &r);
253 if (-1 == r) {
254 cram_free_compression_header(hdr);
255 return NULL;
256 }
257
258 kh_val(hdr->preservation_map, k) = hd;
259 fd->no_ref = !hd.i;
260 break;
261
262 case CRAM_KEY('S','M'):
263 hdr->substitution_matrix[0][(cp[0]>>6)&3] = 'C';
264 hdr->substitution_matrix[0][(cp[0]>>4)&3] = 'G';
265 hdr->substitution_matrix[0][(cp[0]>>2)&3] = 'T';
266 hdr->substitution_matrix[0][(cp[0]>>0)&3] = 'N';
267
268 hdr->substitution_matrix[1][(cp[1]>>6)&3] = 'A';
269 hdr->substitution_matrix[1][(cp[1]>>4)&3] = 'G';
270 hdr->substitution_matrix[1][(cp[1]>>2)&3] = 'T';
271 hdr->substitution_matrix[1][(cp[1]>>0)&3] = 'N';
272
273 hdr->substitution_matrix[2][(cp[2]>>6)&3] = 'A';
274 hdr->substitution_matrix[2][(cp[2]>>4)&3] = 'C';
275 hdr->substitution_matrix[2][(cp[2]>>2)&3] = 'T';
276 hdr->substitution_matrix[2][(cp[2]>>0)&3] = 'N';
277
278 hdr->substitution_matrix[3][(cp[3]>>6)&3] = 'A';
279 hdr->substitution_matrix[3][(cp[3]>>4)&3] = 'C';
280 hdr->substitution_matrix[3][(cp[3]>>2)&3] = 'G';
281 hdr->substitution_matrix[3][(cp[3]>>0)&3] = 'N';
282
283 hdr->substitution_matrix[4][(cp[4]>>6)&3] = 'A';
284 hdr->substitution_matrix[4][(cp[4]>>4)&3] = 'C';
285 hdr->substitution_matrix[4][(cp[4]>>2)&3] = 'G';
286 hdr->substitution_matrix[4][(cp[4]>>0)&3] = 'T';
287
288 hd.p = cp;
289 cp += 5;
290
291 k = kh_put(map, hdr->preservation_map, "SM", &r);
292 if (-1 == r) {
293 cram_free_compression_header(hdr);
294 return NULL;
295 }
296 kh_val(hdr->preservation_map, k) = hd;
297 break;
298
299 case CRAM_KEY('T','D'): {
300 int sz = cram_decode_TD(cp, hdr); // tag dictionary
301 if (sz < 0) {
302 cram_free_compression_header(hdr);
303 return NULL;
304 }
305
306 hd.p = cp;
307 cp += sz;
308
309 k = kh_put(map, hdr->preservation_map, "TD", &r);
310 if (-1 == r) {
311 cram_free_compression_header(hdr);
312 return NULL;
313 }
314 kh_val(hdr->preservation_map, k) = hd;
315 break;
316 }
317
318 default:
319 fprintf(stderr, "Unrecognised preservation map key %c%c\n",
320 cp[-2], cp[-1]);
321 // guess byte;
322 cp++;
323 break;
324 }
325 }
326 if (cp - cp_copy != map_size) {
327 cram_free_compression_header(hdr);
328 return NULL;
329 }
330
331 /* Record encoding map */
332 cp += itf8_get(cp, &map_size); cp_copy = cp;
333 cp += itf8_get(cp, &map_count);
334 for (i = 0; i < map_count; i++) {
335 char *key = cp;
336 int32_t encoding;
337 int32_t size;
338 cram_map *m = malloc(sizeof(*m)); // FIXME: use pooled_alloc
339
340 if (!m) {
341 cram_free_compression_header(hdr);
342 return NULL;
343 }
344
345 cp += 2;
346 cp += itf8_get(cp, &encoding);
347 cp += itf8_get(cp, &size);
348
349 // Fill out cram_map purely for cram_dump to dump out.
350 m->key = (key[0]<<8)|key[1];
351 m->encoding = encoding;
352 m->size = size;
353 m->offset = cp - (char *)b->data;
354 m->codec = NULL;
355
356 if (m->encoding == E_NULL)
357 continue;
358
359 //printf("%s codes for %.2s\n", cram_encoding2str(encoding), key);
360
361 /*
362 * For CRAM1.0 CF and BF are Byte and not Int.
363 * Practically speaking it makes no difference unless we have a
364 * 1.0 format file that stores these in EXTERNAL as only then
365 * does Byte vs Int matter.
366 *
367 * Neither this C code nor Java reference implementations did this,
368 * so we gloss over it and treat them as int.
369 */
370
371 if (key[0] == 'B' && key[1] == 'F') {
372 if (!(hdr->codecs[DS_BF] = cram_decoder_init(encoding, cp, size,
373 E_INT,
374 fd->version))) {
375 cram_free_compression_header(hdr);
376 return NULL;
377 }
378 } else if (key[0] == 'C' && key[1] == 'F') {
379 if (!(hdr->codecs[DS_CF] = cram_decoder_init(encoding, cp, size,
380 E_INT,
381 fd->version))) {
382 cram_free_compression_header(hdr);
383 return NULL;
384 }
385 } else if (key[0] == 'R' && key[1] == 'I') {
386 if (!(hdr->codecs[DS_RI] = cram_decoder_init(encoding, cp, size,
387 E_INT,
388 fd->version))) {
389 cram_free_compression_header(hdr);
390 return NULL;
391 }
392 } else if (key[0] == 'R' && key[1] == 'L') {
393 if (!(hdr->codecs[DS_RL] = cram_decoder_init(encoding, cp, size,
394 E_INT,
395 fd->version))) {
396 cram_free_compression_header(hdr);
397 return NULL;
398 }
399 } else if (key[0] == 'A' && key[1] == 'P') {
400 if (!(hdr->codecs[DS_AP] = cram_decoder_init(encoding, cp, size,
401 E_INT,
402 fd->version))) {
403 cram_free_compression_header(hdr);
404 return NULL;
405 }
406 } else if (key[0] == 'R' && key[1] == 'G') {
407 if (!(hdr->codecs[DS_RG] = cram_decoder_init(encoding, cp, size,
408 E_INT,
409 fd->version))) {
410 cram_free_compression_header(hdr);
411 return NULL;
412 }
413 } else if (key[0] == 'M' && key[1] == 'F') {
414 if (!(hdr->codecs[DS_MF] = cram_decoder_init(encoding, cp, size,
415 E_INT,
416 fd->version))) {
417 cram_free_compression_header(hdr);
418 return NULL;
419 }
420 } else if (key[0] == 'N' && key[1] == 'S') {
421 if (!(hdr->codecs[DS_NS] = cram_decoder_init(encoding, cp, size,
422 E_INT,
423 fd->version))) {
424 cram_free_compression_header(hdr);
425 return NULL;
426 }
427 } else if (key[0] == 'N' && key[1] == 'P') {
428 if (!(hdr->codecs[DS_NP] = cram_decoder_init(encoding, cp, size,
429 E_INT,
430 fd->version))) {
431 cram_free_compression_header(hdr);
432 return NULL;
433 }
434 } else if (key[0] == 'T' && key[1] == 'S') {
435 if (!(hdr->codecs[DS_TS] = cram_decoder_init(encoding, cp, size,
436 E_INT,
437 fd->version))) {
438 cram_free_compression_header(hdr);
439 return NULL;
440 }
441 } else if (key[0] == 'N' && key[1] == 'F') {
442 if (!(hdr->codecs[DS_NF] = cram_decoder_init(encoding, cp, size,
443 E_INT,
444 fd->version))) {
445 cram_free_compression_header(hdr);
446 return NULL;
447 }
448 } else if (key[0] == 'T' && key[1] == 'C') {
449 if (!(hdr->codecs[DS_TC] = cram_decoder_init(encoding, cp, size,
450 E_BYTE,
451 fd->version))) {
452 cram_free_compression_header(hdr);
453 return NULL;
454 }
455 } else if (key[0] == 'T' && key[1] == 'N') {
456 if (!(hdr->codecs[DS_TN] = cram_decoder_init(encoding, cp, size,
457 E_INT,
458 fd->version))) {
459 cram_free_compression_header(hdr);
460 return NULL;
461 }
462 } else if (key[0] == 'F' && key[1] == 'N') {
463 if (!(hdr->codecs[DS_FN] = cram_decoder_init(encoding, cp, size,
464 E_INT,
465 fd->version))) {
466 cram_free_compression_header(hdr);
467 return NULL;
468 }
469 } else if (key[0] == 'F' && key[1] == 'C') {
470 if (!(hdr->codecs[DS_FC] = cram_decoder_init(encoding, cp, size,
471 E_BYTE,
472 fd->version))) {
473 cram_free_compression_header(hdr);
474 return NULL;
475 }
476 } else if (key[0] == 'F' && key[1] == 'P') {
477 if (!(hdr->codecs[DS_FP] = cram_decoder_init(encoding, cp, size,
478 E_INT,
479 fd->version))) {
480 cram_free_compression_header(hdr);
481 return NULL;
482 }
483 } else if (key[0] == 'B' && key[1] == 'S') {
484 if (!(hdr->codecs[DS_BS] = cram_decoder_init(encoding, cp, size,
485 E_BYTE,
486 fd->version))) {
487 cram_free_compression_header(hdr);
488 return NULL;
489 }
490 } else if (key[0] == 'I' && key[1] == 'N') {
491 if (!(hdr->codecs[DS_IN] = cram_decoder_init(encoding, cp, size,
492 E_BYTE_ARRAY,
493 fd->version))) {
494 cram_free_compression_header(hdr);
495 return NULL;
496 }
497 } else if (key[0] == 'S' && key[1] == 'C') {
498 if (!(hdr->codecs[DS_SC] = cram_decoder_init(encoding, cp, size,
499 E_BYTE_ARRAY,
500 fd->version))) {
501 cram_free_compression_header(hdr);
502 return NULL;
503 }
504 } else if (key[0] == 'D' && key[1] == 'L') {
505 if (!(hdr->codecs[DS_DL] = cram_decoder_init(encoding, cp, size,
506 E_INT,
507 fd->version))) {
508 cram_free_compression_header(hdr);
509 return NULL;
510 }
511 } else if (key[0] == 'B' && key[1] == 'A') {
512 if (!(hdr->codecs[DS_BA] = cram_decoder_init(encoding, cp, size,
513 E_BYTE,
514 fd->version))) {
515 cram_free_compression_header(hdr);
516 return NULL;
517 }
518 } else if (key[0] == 'B' && key[1] == 'B') {
519 if (!(hdr->codecs[DS_BB] = cram_decoder_init(encoding, cp, size,
520 E_BYTE_ARRAY,
521 fd->version))) {
522 cram_free_compression_header(hdr);
523 return NULL;
524 }
525 } else if (key[0] == 'R' && key[1] == 'S') {
526 if (!(hdr->codecs[DS_RS] = cram_decoder_init(encoding, cp, size,
527 E_INT,
528 fd->version))) {
529 cram_free_compression_header(hdr);
530 return NULL;
531 }
532 } else if (key[0] == 'P' && key[1] == 'D') {
533 if (!(hdr->codecs[DS_PD] = cram_decoder_init(encoding, cp, size,
534 E_INT,
535 fd->version))) {
536 cram_free_compression_header(hdr);
537 return NULL;
538 }
539 } else if (key[0] == 'H' && key[1] == 'C') {
540 if (!(hdr->codecs[DS_HC] = cram_decoder_init(encoding, cp, size,
541 E_INT,
542 fd->version))) {
543 cram_free_compression_header(hdr);
544 return NULL;
545 }
546 } else if (key[0] == 'M' && key[1] == 'Q') {
547 if (!(hdr->codecs[DS_MQ] = cram_decoder_init(encoding, cp, size,
548 E_INT,
549 fd->version))) {
550 cram_free_compression_header(hdr);
551 return NULL;
552 }
553 } else if (key[0] == 'R' && key[1] == 'N') {
554 if (!(hdr->codecs[DS_RN] = cram_decoder_init(encoding, cp, size,
555 E_BYTE_ARRAY_BLOCK,
556 fd->version))) {
557 cram_free_compression_header(hdr);
558 return NULL;
559 }
560 } else if (key[0] == 'Q' && key[1] == 'S') {
561 if (!(hdr->codecs[DS_QS] = cram_decoder_init(encoding, cp, size,
562 E_BYTE,
563 fd->version))) {
564 cram_free_compression_header(hdr);
565 return NULL;
566 }
567 } else if (key[0] == 'Q' && key[1] == 'Q') {
568 if (!(hdr->codecs[DS_QQ] = cram_decoder_init(encoding, cp, size,
569 E_BYTE_ARRAY,
570 fd->version))) {
571 cram_free_compression_header(hdr);
572 return NULL;
573 }
574 } else if (key[0] == 'T' && key[1] == 'L') {
575 if (!(hdr->codecs[DS_TL] = cram_decoder_init(encoding, cp, size,
576 E_INT,
577 fd->version))) {
578 cram_free_compression_header(hdr);
579 return NULL;
580 }
581 } else if (key[0] == 'T' && key[1] == 'M') {
582 } else if (key[0] == 'T' && key[1] == 'V') {
583 } else
584 fprintf(stderr, "Unrecognised key: %.2s\n", key);
585
586 cp += size;
587
588 m->next = hdr->rec_encoding_map[CRAM_MAP(key[0], key[1])];
589 hdr->rec_encoding_map[CRAM_MAP(key[0], key[1])] = m;
590 }
591 if (cp - cp_copy != map_size) {
592 cram_free_compression_header(hdr);
593 return NULL;
594 }
595
596 /* Tag encoding map */
597 cp += itf8_get(cp, &map_size); cp_copy = cp;
598 cp += itf8_get(cp, &map_count);
599 for (i = 0; i < map_count; i++) {
600 int32_t encoding;
601 int32_t size;
602 cram_map *m = malloc(sizeof(*m)); // FIXME: use pooled_alloc
603 char *key = cp+1;
604
605 if (!m) {
606 cram_free_compression_header(hdr);
607 return NULL;
608 }
609
610 m->key = (key[0]<<16)|(key[1]<<8)|key[2];
611
612 cp += 4; // Strictly ITF8, but this suffices
613 cp += itf8_get(cp, &encoding);
614 cp += itf8_get(cp, &size);
615
616 m->encoding = encoding;
617 m->size = size;
618 m->offset = cp - (char *)b->data;
619 if (!(m->codec = cram_decoder_init(encoding, cp, size,
620 E_BYTE_ARRAY_BLOCK, fd->version))) {
621 cram_free_compression_header(hdr);
622 free(m);
623 return NULL;
624 }
625
626 cp += size;
627
628 m->next = hdr->tag_encoding_map[CRAM_MAP(key[0],key[1])];
629 hdr->tag_encoding_map[CRAM_MAP(key[0],key[1])] = m;
630 }
631 if (cp - cp_copy != map_size) {
632 cram_free_compression_header(hdr);
633 return NULL;
634 }
635
636 return hdr;
637 }
638
639 /*
640 * Note we also need to scan through the record encoding map to
641 * see which data series share the same block, either external or
642 * CORE. For example if we need the BF data series but MQ and CF
643 * are also encoded in the same block then we need to add those in
644 * as a dependency in order to correctly decode BF.
645 *
646 * Returns 0 on success
647 * -1 on failure
648 */
649 int cram_dependent_data_series(cram_fd *fd,
650 cram_block_compression_hdr *hdr,
651 cram_slice *s) {
652 int *block_used;
653 int core_used = 0;
654 int i;
655 static int i_to_id[] = {
656 DS_BF, DS_AP, DS_FP, DS_RL, DS_DL, DS_NF, DS_BA, DS_QS,
657 DS_FC, DS_FN, DS_BS, DS_IN, DS_RG, DS_MQ, DS_TL, DS_RN,
658 DS_NS, DS_NP, DS_TS, DS_MF, DS_CF, DS_RI, DS_RS, DS_PD,
659 DS_HC, DS_SC, DS_BB, DS_QQ,
660 };
661 uint32_t orig_ds;
662
663 /*
664 * Set the data_series bit field based on fd->required_fields
665 * contents.
666 */
667 if (fd->required_fields && fd->required_fields != INT_MAX) {
668 hdr->data_series = 0;
669
670 if (fd->required_fields & SAM_QNAME)
671 hdr->data_series |= CRAM_RN;
672
673 if (fd->required_fields & SAM_FLAG)
674 hdr->data_series |= CRAM_BF;
675
676 if (fd->required_fields & SAM_RNAME)
677 hdr->data_series |= CRAM_RI | CRAM_BF;
678
679 if (fd->required_fields & SAM_POS)
680 hdr->data_series |= CRAM_AP | CRAM_BF;
681
682 if (fd->required_fields & SAM_MAPQ)
683 hdr->data_series |= CRAM_MQ;
684
685 if (fd->required_fields & SAM_CIGAR)
686 hdr->data_series |= CRAM_CIGAR;
687
688 if (fd->required_fields & SAM_RNEXT)
689 hdr->data_series |= CRAM_CF | CRAM_NF | CRAM_RI | CRAM_NS |CRAM_BF;
690
691 if (fd->required_fields & SAM_PNEXT)
692 hdr->data_series |= CRAM_CF | CRAM_NF | CRAM_AP | CRAM_NP | CRAM_BF;
693
694 if (fd->required_fields & SAM_TLEN)
695 hdr->data_series |= CRAM_CF | CRAM_NF | CRAM_AP | CRAM_TS |
696 CRAM_BF | CRAM_MF | CRAM_RI | CRAM_CIGAR;
697
698 if (fd->required_fields & SAM_SEQ)
699 hdr->data_series |= CRAM_SEQ;
700
701 if (!(fd->required_fields & SAM_AUX))
702 // No easy way to get MD/NM without other tags at present
703 fd->decode_md = 0;
704
705 if (fd->required_fields & SAM_QUAL)
706 hdr->data_series |= CRAM_SEQ;
707
708 if (fd->required_fields & SAM_AUX)
709 hdr->data_series |= CRAM_RG | CRAM_TL | CRAM_aux;
710
711 if (fd->required_fields & SAM_RGAUX)
712 hdr->data_series |= CRAM_RG | CRAM_BF;
713
714 // Always uncompress CORE block
715 if (cram_uncompress_block(s->block[0]))
716 return -1;
717 } else {
718 hdr->data_series = CRAM_ALL;
719
720 for (i = 0; i < s->hdr->num_blocks; i++) {
721 if (cram_uncompress_block(s->block[i]))
722 return -1;
723 }
724
725 return 0;
726 }
727
728 block_used = calloc(s->hdr->num_blocks+1, sizeof(int));
729 if (!block_used)
730 return -1;
731
732 do {
733 /*
734 * Also set data_series based on code prerequisites. Eg if we need
735 * CRAM_QS then we also need to know CRAM_RL so we know how long it
736 * is, or if we need FC/FP then we also need FN (number of features).
737 *
738 * It's not reciprocal though. We may be needing to decode FN
739 * but have no need to decode FC, FP and cigar ops.
740 */
741 if (hdr->data_series & CRAM_RS) hdr->data_series |= CRAM_FC|CRAM_FP;
742 if (hdr->data_series & CRAM_PD) hdr->data_series |= CRAM_FC|CRAM_FP;
743 if (hdr->data_series & CRAM_HC) hdr->data_series |= CRAM_FC|CRAM_FP;
744 if (hdr->data_series & CRAM_QS) hdr->data_series |= CRAM_FC|CRAM_FP;
745 if (hdr->data_series & CRAM_IN) hdr->data_series |= CRAM_FC|CRAM_FP;
746 if (hdr->data_series & CRAM_SC) hdr->data_series |= CRAM_FC|CRAM_FP;
747 if (hdr->data_series & CRAM_BS) hdr->data_series |= CRAM_FC|CRAM_FP;
748 if (hdr->data_series & CRAM_DL) hdr->data_series |= CRAM_FC|CRAM_FP;
749 if (hdr->data_series & CRAM_BA) hdr->data_series |= CRAM_FC|CRAM_FP;
750 if (hdr->data_series & CRAM_BB) hdr->data_series |= CRAM_FC|CRAM_FP;
751 if (hdr->data_series & CRAM_QQ) hdr->data_series |= CRAM_FC|CRAM_FP;
752
753 // cram_decode_seq() needs seq[] array
754 if (hdr->data_series & (CRAM_SEQ|CRAM_CIGAR)) hdr->data_series |= CRAM_RL;
755
756 if (hdr->data_series & CRAM_FP) hdr->data_series |= CRAM_FC;
757 if (hdr->data_series & CRAM_FC) hdr->data_series |= CRAM_FN;
758 if (hdr->data_series & CRAM_aux) hdr->data_series |= CRAM_TL;
759 if (hdr->data_series & CRAM_MF) hdr->data_series |= CRAM_CF;
760 if (hdr->data_series & CRAM_MQ) hdr->data_series |= CRAM_BF;
761 if (hdr->data_series & CRAM_BS) hdr->data_series |= CRAM_RI;
762 if (hdr->data_series & (CRAM_MF |CRAM_NS |CRAM_NP |CRAM_TS |CRAM_NF))
763 hdr->data_series |= CRAM_CF;
764 if (!hdr->read_names_included && hdr->data_series & CRAM_RN)
765 hdr->data_series |= CRAM_CF | CRAM_NF;
766 if (hdr->data_series & (CRAM_BA | CRAM_QS | CRAM_BB | CRAM_QQ))
767 hdr->data_series |= CRAM_BF | CRAM_CF | CRAM_RL;
768
769 orig_ds = hdr->data_series;
770
771 // Find which blocks are in use.
772 for (i = 0; i < sizeof(i_to_id)/sizeof(*i_to_id); i++) {
773 int bnum1, bnum2, j;
774 cram_codec *c = hdr->codecs[i_to_id[i]];
775
776 if (!(hdr->data_series & (1<<i)))
777 continue;
778
779 if (!c)
780 continue;
781
782 bnum1 = cram_codec_to_id(c, &bnum2);
783
784 for (;;) {
785 switch (bnum1) {
786 case -2:
787 break;
788
789 case -1:
790 core_used = 1;
791 break;
792
793 default:
794 for (j = 0; j < s->hdr->num_blocks; j++) {
795 if (s->block[j]->content_type == EXTERNAL &&
796 s->block[j]->content_id == bnum1) {
797 block_used[j] = 1;
798 if (cram_uncompress_block(s->block[j])) {
799 free(block_used);
800 return -1;
801 }
802 }
803 }
804 break;
805 }
806
807 if (bnum2 == -2 || bnum1 == bnum2)
808 break;
809
810 bnum1 = bnum2; // 2nd pass
811 }
812 }
813
814 // Tags too
815 if ((fd->required_fields & SAM_AUX) ||
816 (hdr->data_series & CRAM_aux)) {
817 for (i = 0; i < CRAM_MAP_HASH; i++) {
818 int bnum1, bnum2, j;
819 cram_map *m = hdr->tag_encoding_map[i];
820
821 while (m) {
822 cram_codec *c = m->codec;
823 if (!c)
824 continue;
825
826 bnum1 = cram_codec_to_id(c, &bnum2);
827
828 for (;;) {
829 switch (bnum1) {
830 case -2:
831 break;
832
833 case -1:
834 core_used = 1;
835 break;
836
837 default:
838 for (j = 0; j < s->hdr->num_blocks; j++) {
839 if (s->block[j]->content_type == EXTERNAL &&
840 s->block[j]->content_id == bnum1) {
841 block_used[j] = 1;
842 if (cram_uncompress_block(s->block[j])) {
843 free(block_used);
844 return -1;
845 }
846 }
847 }
848 break;
849 }
850
851 if (bnum2 == -2 || bnum1 == bnum2)
852 break;
853
854 bnum1 = bnum2; // 2nd pass
855 }
856
857 m = m->next;
858 }
859 }
860 }
861
862 // We now know which blocks are in used, so repeat and find
863 // which other data series need to be added.
864 for (i = 0; i < sizeof(i_to_id)/sizeof(*i_to_id); i++) {
865 int bnum1, bnum2, j;
866 cram_codec *c = hdr->codecs[i_to_id[i]];
867
868 if (!c)
869 continue;
870
871 bnum1 = cram_codec_to_id(c, &bnum2);
872
873 for (;;) {
874 switch (bnum1) {
875 case -2:
876 break;
877
878 case -1:
879 if (core_used) {
880 //printf(" + data series %08x:\n", 1<<i);
881 hdr->data_series |= 1<<i;
882 }
883 break;
884
885 default:
886 for (j = 0; j < s->hdr->num_blocks; j++) {
887 if (s->block[j]->content_type == EXTERNAL &&
888 s->block[j]->content_id == bnum1) {
889 if (block_used[j]) {
890 //printf(" + data series %08x:\n", 1<<i);
891 hdr->data_series |= 1<<i;
892 }
893 }
894 }
895 break;
896 }
897
898 if (bnum2 == -2 || bnum1 == bnum2)
899 break;
900
901 bnum1 = bnum2; // 2nd pass
902 }
903 }
904
905 // Tags too
906 for (i = 0; i < CRAM_MAP_HASH; i++) {
907 int bnum1, bnum2, j;
908 cram_map *m = hdr->tag_encoding_map[i];
909
910 while (m) {
911 cram_codec *c = m->codec;
912 if (!c)
913 continue;
914
915 bnum1 = cram_codec_to_id(c, &bnum2);
916
917 for (;;) {
918 switch (bnum1) {
919 case -2:
920 break;
921
922 case -1:
923 //printf(" + data series %08x:\n", CRAM_aux);
924 hdr->data_series |= CRAM_aux;
925 break;
926
927 default:
928 for (j = 0; j < s->hdr->num_blocks; j++) {
929 if (s->block[j]->content_type &&
930 s->block[j]->content_id == bnum1) {
931 if (block_used[j]) {
932 //printf(" + data series %08x:\n",
933 // CRAM_aux);
934 hdr->data_series |= CRAM_aux;
935 }
936 }
937 }
938 break;
939 }
940
941 if (bnum2 == -2 || bnum1 == bnum2)
942 break;
943
944 bnum1 = bnum2; // 2nd pass
945 }
946
947 m = m->next;
948 }
949 }
950 } while (orig_ds != hdr->data_series);
951
952 free(block_used);
953 return 0;
954 }
955
956 /* ----------------------------------------------------------------------
957 * CRAM slices
958 */
959
960 /*
961 * Decodes a CRAM (un)mapped slice header block.
962 * Returns slice header ptr on success
963 * NULL on failure
964 */
965 cram_block_slice_hdr *cram_decode_slice_header(cram_fd *fd, cram_block *b) {
966 cram_block_slice_hdr *hdr;
967 char *cp = (char *)b->data;
968 int i;
969
970 if (b->content_type != MAPPED_SLICE &&
971 b->content_type != UNMAPPED_SLICE)
972 return NULL;
973
974 if (!(hdr = calloc(1, sizeof(*hdr))))
975 return NULL;
976
977 hdr->content_type = b->content_type;
978
979 if (b->content_type == MAPPED_SLICE) {
980 cp += itf8_get(cp, &hdr->ref_seq_id);
981 cp += itf8_get(cp, &hdr->ref_seq_start);
982 cp += itf8_get(cp, &hdr->ref_seq_span);
983 }
984 cp += itf8_get(cp, &hdr->num_records);
985 hdr->record_counter = 0;
986 if (CRAM_MAJOR_VERS(fd->version) == 2) {
987 int32_t i32;
988 cp += itf8_get(cp, &i32);
989 hdr->record_counter = i32;
990 } else if (CRAM_MAJOR_VERS(fd->version) >= 3) {
991 cp += ltf8_get(cp, &hdr->record_counter);
992 }
993
994 cp += itf8_get(cp, &hdr->num_blocks);
995
996 cp += itf8_get(cp, &hdr->num_content_ids);
997 hdr->block_content_ids = malloc(hdr->num_content_ids * sizeof(int32_t));
998 if (!hdr->block_content_ids) {
999 free(hdr);
1000 return NULL;
1001 }
1002
1003 for (i = 0; i < hdr->num_content_ids; i++) {
1004 cp += itf8_get(cp, &hdr->block_content_ids[i]);
1005 }
1006
1007 if (b->content_type == MAPPED_SLICE) {
1008 cp += itf8_get(cp, &hdr->ref_base_id);
1009 }
1010
1011 if (CRAM_MAJOR_VERS(fd->version) != 1) {
1012 memcpy(hdr->md5, cp, 16);
1013 } else {
1014 memset(hdr->md5, 0, 16);
1015 }
1016
1017 return hdr;
1018 }
1019
1020
1021 #if 0
1022 /* Returns the number of bits set in val; it the highest bit used */
1023 static int nbits(int v) {
1024 static const int MultiplyDeBruijnBitPosition[32] = {
1025 1, 10, 2, 11, 14, 22, 3, 30, 12, 15, 17, 19, 23, 26, 4, 31,
1026 9, 13, 21, 29, 16, 18, 25, 8, 20, 28, 24, 7, 27, 6, 5, 32
1027 };
1028
1029 v |= v >> 1; // first up to set all bits 1 after the first 1 */
1030 v |= v >> 2;
1031 v |= v >> 4;
1032 v |= v >> 8;
1033 v |= v >> 16;
1034
1035 // DeBruijn magic to find top bit
1036 return MultiplyDeBruijnBitPosition[(uint32_t)(v * 0x07C4ACDDU) >> 27];
1037 }
1038 #endif
1039
1040 #if 0
1041 static int sort_freqs(const void *vp1, const void *vp2) {
1042 const int i1 = *(const int *)vp1;
1043 const int i2 = *(const int *)vp2;
1044 return i1-i2;
1045 }
1046 #endif
1047
1048 /* ----------------------------------------------------------------------
1049 * Primary CRAM sequence decoder
1050 */
1051
1052 /*
1053 * Internal part of cram_decode_slice().
1054 * Generates the sequence, quality and cigar components.
1055 */
1056 static int cram_decode_seq(cram_fd *fd, cram_container *c, cram_slice *s,
1057 cram_block *blk, cram_record *cr, SAM_hdr *bfd,
1058 int cf, char *seq, char *qual) {
1059 int prev_pos = 0, f, r = 0, out_sz = 1;
1060 int seq_pos = 1;
1061 int cig_len = 0, ref_pos = cr->apos;
1062 int32_t fn, i32;
1063 enum cigar_op cig_op = BAM_CMATCH;
1064 uint32_t *cigar = s->cigar;
1065 uint32_t ncigar = s->ncigar;
1066 uint32_t cigar_alloc = s->cigar_alloc;
1067 uint32_t nm = 0, md_dist = 0;
1068 int orig_aux = 0;
1069 int decode_md = fd->decode_md && s->ref;
1070 uint32_t ds = c->comp_hdr->data_series;
1071
1072 if ((ds & CRAM_QS) && !(cf & CRAM_FLAG_PRESERVE_QUAL_SCORES)) {
1073 memset(qual, 30, cr->len);
1074 }
1075
1076 if (decode_md) {
1077 orig_aux = BLOCK_SIZE(s->aux_blk);
1078 BLOCK_APPEND(s->aux_blk, "MDZ", 3);
1079 }
1080
1081 if (ds & CRAM_FN) {
1082 if (!c->comp_hdr->codecs[DS_FN]) return -1;
1083 r |= c->comp_hdr->codecs[DS_FN]->decode(s,c->comp_hdr->codecs[DS_FN],
1084 blk, (char *)&fn, &out_sz);
1085 } else {
1086 fn = 0;
1087 }
1088
1089 ref_pos--; // count from 0
1090 cr->cigar = ncigar;
1091
1092 if (!(ds & (CRAM_FC | CRAM_FP)))
1093 goto skip_cigar;
1094
1095 for (f = 0; f < fn; f++) {
1096 int32_t pos = 0;
1097 char op;
1098
1099 if (ncigar+2 >= cigar_alloc) {
1100 cigar_alloc = cigar_alloc ? cigar_alloc*2 : 1024;
1101 s->cigar = cigar;
1102 if (!(cigar = realloc(cigar, cigar_alloc * sizeof(*cigar))))
1103 return -1;
1104 }
1105
1106 if (ds & CRAM_FC) {
1107 if (!c->comp_hdr->codecs[DS_FC]) return -1;
1108 r |= c->comp_hdr->codecs[DS_FC]->decode(s,
1109 c->comp_hdr->codecs[DS_FC],
1110 blk,
1111 &op, &out_sz);
1112 }
1113
1114 if (!(ds & CRAM_FP))
1115 continue;
1116
1117 if (!c->comp_hdr->codecs[DS_FP]) return -1;
1118 r |= c->comp_hdr->codecs[DS_FP]->decode(s,
1119 c->comp_hdr->codecs[DS_FP],
1120 blk,
1121 (char *)&pos, &out_sz);
1122 pos += prev_pos;
1123
1124 if (pos > seq_pos) {
1125 if (pos > cr->len+1)
1126 return -1;
1127
1128 if (s->ref && cr->ref_id >= 0) {
1129 if (ref_pos + pos - seq_pos > bfd->ref[cr->ref_id].len) {
1130 static int whinged = 0;
1131 if (!whinged)
1132 fprintf(stderr, "Ref pos outside of ref "
1133 "sequence boundary\n");
1134 whinged = 1;
1135 } else {
1136 memcpy(&seq[seq_pos-1], &s->ref[ref_pos - s->ref_start +1],
1137 pos - seq_pos);
1138 }
1139 }
1140 #ifdef USE_X
1141 if (cig_len && cig_op != BAM_CBASE_MATCH) {
1142 cigar[ncigar++] = (cig_len<<4) + cig_op;
1143 cig_len = 0;
1144 }
1145 cig_op = BAM_CBASE_MATCH;
1146 #else
1147 if (cig_len && cig_op != BAM_CMATCH) {
1148 cigar[ncigar++] = (cig_len<<4) + cig_op;
1149 cig_len = 0;
1150 }
1151 cig_op = BAM_CMATCH;
1152 #endif
1153 cig_len += pos - seq_pos;
1154 ref_pos += pos - seq_pos;
1155 md_dist += pos - seq_pos;
1156 seq_pos = pos;
1157 }
1158
1159 prev_pos = pos;
1160
1161 if (!(ds & CRAM_FC))
1162 goto skip_cigar;
1163
1164 if (!(ds & CRAM_FC))
1165 continue;
1166
1167 switch(op) {
1168 case 'S': { // soft clip: IN
1169 int32_t out_sz2 = 1;
1170
1171 if (cig_len) {
1172 cigar[ncigar++] = (cig_len<<4) + cig_op;
1173 cig_len = 0;
1174 }
1175 if (ds & CRAM_IN) {
1176 switch (CRAM_MAJOR_VERS(fd->version)) {
1177 case 1:
1178 r |= c->comp_hdr->codecs[DS_IN]
1179 ? c->comp_hdr->codecs[DS_IN]
1180 ->decode(s, c->comp_hdr->codecs[DS_IN],
1181 blk, &seq[pos-1], &out_sz2)
1182 : (seq[pos-1] = 'N', out_sz2 = 1, 0);
1183 break;
1184
1185 case 2:
1186 default:
1187 r |= c->comp_hdr->codecs[DS_SC]
1188 ? c->comp_hdr->codecs[DS_SC]
1189 ->decode(s, c->comp_hdr->codecs[DS_SC],
1190 blk, &seq[pos-1], &out_sz2)
1191 : (seq[pos-1] = 'N', out_sz2 = 1, 0);
1192 break;
1193
1194 // default:
1195 // r |= c->comp_hdr->codecs[DS_BB]
1196 // ? c->comp_hdr->codecs[DS_BB]
1197 // ->decode(s, c->comp_hdr->codecs[DS_BB],
1198 // blk, &seq[pos-1], &out_sz2)
1199 // : (seq[pos-1] = 'N', out_sz2 = 1, 0);
1200 }
1201 cigar[ncigar++] = (out_sz2<<4) + BAM_CSOFT_CLIP;
1202 cig_op = BAM_CSOFT_CLIP;
1203 seq_pos += out_sz2;
1204 }
1205 break;
1206 }
1207
1208 case 'X': { // Substitution; BS
1209 unsigned char base;
1210 #ifdef USE_X
1211 if (cig_len && cig_op != BAM_CBASE_MISMATCH) {
1212 cigar[ncigar++] = (cig_len<<4) + cig_op;
1213 cig_len = 0;
1214 }
1215 if (ds & CRAM_BS) {
1216 if (!c->comp_hdr->codecs[DS_BS]) return -1;
1217 r |= c->comp_hdr->codecs[DS_BS]
1218 ->decode(s, c->comp_hdr->codecs[DS_BS], blk,
1219 (char *)&base, &out_sz);
1220 seq[pos-1] = 'N'; // FIXME look up BS=base value
1221 }
1222 cig_op = BAM_CBASE_MISMATCH;
1223 #else
1224 int ref_base;
1225 if (cig_len && cig_op != BAM_CMATCH) {
1226 cigar[ncigar++] = (cig_len<<4) + cig_op;
1227 cig_len = 0;
1228 }
1229 if (ds & CRAM_BS) {
1230 if (!c->comp_hdr->codecs[DS_BS]) return -1;
1231 r |= c->comp_hdr->codecs[DS_BS]
1232 ->decode(s, c->comp_hdr->codecs[DS_BS], blk,
1233 (char *)&base, &out_sz);
1234 if (ref_pos >= bfd->ref[cr->ref_id].len || !s->ref) {
1235 seq[pos-1] = 'N';
1236 } else {
1237 ref_base = fd->L1[(uc)s->ref[ref_pos - s->ref_start +1]];
1238 seq[pos-1] = c->comp_hdr->
1239 substitution_matrix[ref_base][base];
1240 if (decode_md) {
1241 BLOCK_APPEND_UINT(s->aux_blk, md_dist);
1242 BLOCK_APPEND_CHAR(s->aux_blk,
1243 s->ref[ref_pos-s->ref_start +1]);
1244 md_dist = 0;
1245 }
1246 }
1247 }
1248 cig_op = BAM_CMATCH;
1249 #endif
1250 nm++;
1251 cig_len++;
1252 seq_pos++;
1253 ref_pos++;
1254 break;
1255 }
1256
1257 case 'D': { // Deletion; DL
1258 if (cig_len && cig_op != BAM_CDEL) {
1259 cigar[ncigar++] = (cig_len<<4) + cig_op;
1260 cig_len = 0;
1261 }
1262 if (ds & CRAM_DL) {
1263 if (!c->comp_hdr->codecs[DS_DL]) return -1;
1264 r |= c->comp_hdr->codecs[DS_DL]
1265 ->decode(s, c->comp_hdr->codecs[DS_DL], blk,
1266 (char *)&i32, &out_sz);
1267 if (decode_md) {
1268 BLOCK_APPEND_UINT(s->aux_blk, md_dist);
1269 BLOCK_APPEND_CHAR(s->aux_blk, '^');
1270 BLOCK_APPEND(s->aux_blk,
1271 &s->ref[ref_pos - s->ref_start +1],
1272 i32);
1273 md_dist = 0;
1274 }
1275 cig_op = BAM_CDEL;
1276 cig_len += i32;
1277 ref_pos += i32;
1278 nm += i32;
1279 //printf(" %d: DL = %d (ret %d)\n", f, i32, r);
1280 }
1281 break;
1282 }
1283
1284 case 'I': { // Insertion (several bases); IN
1285 int32_t out_sz2 = 1;
1286
1287 if (cig_len && cig_op != BAM_CINS) {
1288 cigar[ncigar++] = (cig_len<<4) + cig_op;
1289 cig_len = 0;
1290 }
1291
1292 if (ds & CRAM_IN) {
1293 if (!c->comp_hdr->codecs[DS_IN]) return -1;
1294 r |= c->comp_hdr->codecs[DS_IN]
1295 ->decode(s, c->comp_hdr->codecs[DS_IN], blk,
1296 &seq[pos-1], &out_sz2);
1297 cig_op = BAM_CINS;
1298 cig_len += out_sz2;
1299 seq_pos += out_sz2;
1300 nm += out_sz2;
1301 //printf(" %d: IN(I) = %.*s (ret %d, out_sz %d)\n", f, out_sz2, dat, r, out_sz2);
1302 }
1303 break;
1304 }
1305
1306 case 'i': { // Insertion (single base); BA
1307 if (cig_len && cig_op != BAM_CINS) {
1308 cigar[ncigar++] = (cig_len<<4) + cig_op;
1309 cig_len = 0;
1310 }
1311 if (ds & CRAM_BA) {
1312 if (!c->comp_hdr->codecs[DS_BA]) return -1;
1313 r |= c->comp_hdr->codecs[DS_BA]
1314 ->decode(s, c->comp_hdr->codecs[DS_BA], blk,
1315 (char *)&seq[pos-1], &out_sz);
1316 }
1317 cig_op = BAM_CINS;
1318 cig_len++;
1319 seq_pos++;
1320 nm++;
1321 break;
1322 }
1323
1324 case 'b': { // Several bases
1325 int32_t len = 1;
1326
1327 if (cig_len && cig_op != BAM_CMATCH) {
1328 cigar[ncigar++] = (cig_len<<4) + cig_op;
1329 cig_len = 0;
1330 }
1331
1332 if (ds & CRAM_BB) {
1333 if (!c->comp_hdr->codecs[DS_BB]) return -1;
1334 r |= c->comp_hdr->codecs[DS_BB]
1335 ->decode(s, c->comp_hdr->codecs[DS_BB], blk,
1336 (char *)&seq[pos-1], &len);
1337 }
1338
1339 cig_op = BAM_CMATCH;
1340
1341 cig_len+=len;
1342 seq_pos+=len;
1343 ref_pos+=len;
1344 //prev_pos+=len;
1345 break;
1346 }
1347
1348 case 'q': { // Several quality values
1349 int32_t len = 1;
1350
1351 if (cig_len && cig_op != BAM_CMATCH) {
1352 cigar[ncigar++] = (cig_len<<4) + cig_op;
1353 cig_len = 0;
1354 }
1355
1356 if (ds & CRAM_QQ) {
1357 if (!c->comp_hdr->codecs[DS_QQ]) return -1;
1358 r |= c->comp_hdr->codecs[DS_QQ]
1359 ->decode(s, c->comp_hdr->codecs[DS_QQ], blk,
1360 (char *)&qual[pos-1], &len);
1361 }
1362
1363 cig_op = BAM_CMATCH;
1364
1365 cig_len+=len;
1366 seq_pos+=len;
1367 ref_pos+=len;
1368 //prev_pos+=len;
1369 break;
1370 }
1371
1372 case 'B': { // Read base; BA, QS
1373 #ifdef USE_X
1374 if (cig_len && cig_op != BAM_CBASE_MISMATCH) {
1375 cigar[ncigar++] = (cig_len<<4) + cig_op;
1376 cig_len = 0;
1377 }
1378 #else
1379 if (cig_len && cig_op != BAM_CMATCH) {
1380 cigar[ncigar++] = (cig_len<<4) + cig_op;
1381 cig_len = 0;
1382 }
1383 #endif
1384 if (ds & CRAM_BA) {
1385 if (!c->comp_hdr->codecs[DS_BA]) return -1;
1386 r |= c->comp_hdr->codecs[DS_BA]
1387 ->decode(s, c->comp_hdr->codecs[DS_BA], blk,
1388 (char *)&seq[pos-1], &out_sz);
1389 }
1390 if (ds & CRAM_QS) {
1391 if (!c->comp_hdr->codecs[DS_QS]) return -1;
1392 r |= c->comp_hdr->codecs[DS_QS]
1393 ->decode(s, c->comp_hdr->codecs[DS_QS], blk,
1394 (char *)&qual[pos-1], &out_sz);
1395 }
1396 #ifdef USE_X
1397 cig_op = BAM_CBASE_MISMATCH;
1398 #else
1399 cig_op = BAM_CMATCH;
1400 #endif
1401 cig_len++;
1402 seq_pos++;
1403 ref_pos++;
1404 //printf(" %d: BA/QS(B) = %c/%d (ret %d)\n", f, i32, qc, r);
1405 break;
1406 }
1407
1408 case 'Q': { // Quality score; QS
1409 if (ds & CRAM_QS) {
1410 if (!c->comp_hdr->codecs[DS_QS]) return -1;
1411 r |= c->comp_hdr->codecs[DS_QS]
1412 ->decode(s, c->comp_hdr->codecs[DS_QS], blk,
1413 (char *)&qual[pos-1], &out_sz);
1414 //printf(" %d: QS = %d (ret %d)\n", f, qc, r);
1415 }
1416 break;
1417 }
1418
1419 case 'H': { // hard clip; HC
1420 if (cig_len && cig_op != BAM_CHARD_CLIP) {
1421 cigar[ncigar++] = (cig_len<<4) + cig_op;
1422 cig_len = 0;
1423 }
1424 if (ds & CRAM_HC) {
1425 if (!c->comp_hdr->codecs[DS_HC]) return -1;
1426 r |= c->comp_hdr->codecs[DS_HC]
1427 ->decode(s, c->comp_hdr->codecs[DS_HC], blk,
1428 (char *)&i32, &out_sz);
1429 cig_op = BAM_CHARD_CLIP;
1430 cig_len += i32;
1431 nm += i32;
1432 }
1433 break;
1434 }
1435
1436 case 'P': { // padding; PD
1437 if (cig_len && cig_op != BAM_CPAD) {
1438 cigar[ncigar++] = (cig_len<<4) + cig_op;
1439 cig_len = 0;
1440 }
1441 if (ds & CRAM_PD) {
1442 if (!c->comp_hdr->codecs[DS_PD]) return -1;
1443 r |= c->comp_hdr->codecs[DS_PD]
1444 ->decode(s, c->comp_hdr->codecs[DS_PD], blk,
1445 (char *)&i32, &out_sz);
1446 cig_op = BAM_CPAD;
1447 cig_len += i32;
1448 nm += i32;
1449 }
1450 break;
1451 }
1452
1453 case 'N': { // Ref skip; RS
1454 if (cig_len && cig_op != BAM_CREF_SKIP) {
1455 cigar[ncigar++] = (cig_len<<4) + cig_op;
1456 cig_len = 0;
1457 }
1458 if (ds & CRAM_RS) {
1459 if (!c->comp_hdr->codecs[DS_RS]) return -1;
1460 r |= c->comp_hdr->codecs[DS_RS]
1461 ->decode(s, c->comp_hdr->codecs[DS_RS], blk,
1462 (char *)&i32, &out_sz);
1463 cig_op = BAM_CREF_SKIP;
1464 cig_len += i32;
1465 ref_pos += i32;
1466 nm += i32;
1467 }
1468 break;
1469 }
1470
1471 default:
1472 abort();
1473 }
1474 }
1475
1476 if (!(ds & CRAM_FC))
1477 goto skip_cigar;
1478
1479 /* An implement match op for any unaccounted for bases */
1480 if ((ds & CRAM_FN) && cr->len >= seq_pos) {
1481 if (s->ref) {
1482 if (ref_pos + cr->len - seq_pos + 1 > bfd->ref[cr->ref_id].len) {
1483 static int whinged = 0;
1484 if (!whinged)
1485 fprintf(stderr, "Ref pos outside of ref sequence boundary\n");
1486 whinged = 1;
1487 } else {
1488 memcpy(&seq[seq_pos-1], &s->ref[ref_pos - s->ref_start +1],
1489 cr->len - seq_pos + 1);
1490 ref_pos += cr->len - seq_pos + 1;
1491 md_dist += cr->len - seq_pos + 1;
1492 }
1493 }
1494
1495 if (ncigar+1 >= cigar_alloc) {
1496 cigar_alloc = cigar_alloc ? cigar_alloc*2 : 1024;
1497 s->cigar = cigar;
1498 if (!(cigar = realloc(cigar, cigar_alloc * sizeof(*cigar))))
1499 return -1;
1500 }
1501 #ifdef USE_X
1502 if (cig_len && cig_op != BAM_CBASE_MATCH) {
1503 cigar[ncigar++] = (cig_len<<4) + cig_op;
1504 cig_len = 0;
1505 }
1506 cig_op = BAM_CBASE_MATCH;
1507 #else
1508 if (cig_len && cig_op != BAM_CMATCH) {
1509 cigar[ncigar++] = (cig_len<<4) + cig_op;
1510 cig_len = 0;
1511 }
1512 cig_op = BAM_CMATCH;
1513 #endif
1514 cig_len += cr->len - seq_pos+1;
1515 }
1516
1517 skip_cigar:
1518
1519 if ((ds & CRAM_FN) && decode_md) {
1520 BLOCK_APPEND_UINT(s->aux_blk, md_dist);
1521 }
1522
1523 if (cig_len) {
1524 if (ncigar >= cigar_alloc) {
1525 cigar_alloc = cigar_alloc ? cigar_alloc*2 : 1024;
1526 s->cigar = cigar;
1527 if (!(cigar = realloc(cigar, cigar_alloc * sizeof(*cigar))))
1528 return -1;
1529 }
1530
1531 cigar[ncigar++] = (cig_len<<4) + cig_op;
1532 }
1533
1534 cr->ncigar = ncigar - cr->cigar;
1535 cr->aend = ref_pos;
1536
1537 //printf("2: %.*s %d .. %d\n", cr->name_len, DSTRING_STR(name_ds) + cr->name, cr->apos, ref_pos);
1538
1539 if (ds & CRAM_MQ) {
1540 if (!c->comp_hdr->codecs[DS_MQ]) return -1;
1541 r |= c->comp_hdr->codecs[DS_MQ]
1542 ->decode(s, c->comp_hdr->codecs[DS_MQ], blk,
1543 (char *)&cr->mqual, &out_sz);
1544 } else {
1545 cr->mqual = 40;
1546 }
1547
1548 if ((ds & CRAM_QS) && (cf & CRAM_FLAG_PRESERVE_QUAL_SCORES)) {
1549 int32_t out_sz2 = cr->len;
1550
1551 if (ds & CRAM_QS) {
1552 if (!c->comp_hdr->codecs[DS_QS]) return -1;
1553 r |= c->comp_hdr->codecs[DS_QS]
1554 ->decode(s, c->comp_hdr->codecs[DS_QS], blk,
1555 qual, &out_sz2);
1556 }
1557 }
1558
1559 s->cigar = cigar;
1560 s->cigar_alloc = cigar_alloc;
1561 s->ncigar = ncigar;
1562
1563 if (decode_md) {
1564 char buf[7];
1565 BLOCK_APPEND_CHAR(s->aux_blk, '\0'); // null terminate MD:Z:
1566 cr->aux_size += BLOCK_SIZE(s->aux_blk) - orig_aux;
1567 buf[0] = 'N'; buf[1] = 'M'; buf[2] = 'I';
1568 buf[3] = (nm>> 0) & 0xff;
1569 buf[4] = (nm>> 8) & 0xff;
1570 buf[5] = (nm>>16) & 0xff;
1571 buf[6] = (nm>>24) & 0xff;
1572 BLOCK_APPEND(s->aux_blk, buf, 7);
1573 cr->aux_size += 7;
1574 }
1575
1576 return r;
1577 }
1578
1579 /*
1580 * Quick and simple hash lookup for cram_map arrays
1581 */
1582 static cram_map *map_find(cram_map **map, unsigned char *key, int id) {
1583 cram_map *m;
1584
1585 m = map[CRAM_MAP(key[0],key[1])];
1586 while (m && m->key != id)
1587 m= m->next;
1588
1589 return m;
1590 }
1591
1592 //#define map_find(M,K,I) M[CRAM_MAP(K[0],K[1])];while (m && m->key != I);m= m->next
1593
1594
1595 static int cram_decode_aux_1_0(cram_container *c, cram_slice *s,
1596 cram_block *blk, cram_record *cr) {
1597 int i, r = 0, out_sz = 1;
1598 unsigned char ntags;
1599
1600 if (!c->comp_hdr->codecs[DS_TC]) return -1;
1601 r |= c->comp_hdr->codecs[DS_TC]->decode(s, c->comp_hdr->codecs[DS_TC], blk,
1602 (char *)&ntags, &out_sz);
1603 cr->ntags = ntags;
1604
1605 //printf("TC=%d\n", cr->ntags);
1606 cr->aux_size = 0;
1607 cr->aux = BLOCK_SIZE(s->aux_blk);
1608
1609 for (i = 0; i < cr->ntags; i++) {
1610 int32_t id, out_sz = 1;
1611 unsigned char tag_data[3];
1612 cram_map *m;
1613
1614 //printf("Tag %d/%d\n", i+1, cr->ntags);
1615 if (!c->comp_hdr->codecs[DS_TN]) return -1;
1616 r |= c->comp_hdr->codecs[DS_TN]->decode(s, c->comp_hdr->codecs[DS_TN],
1617 blk, (char *)&id, &out_sz);
1618 if (out_sz == 3) {
1619 tag_data[0] = ((char *)&id)[0];
1620 tag_data[1] = ((char *)&id)[1];
1621 tag_data[2] = ((char *)&id)[2];
1622 } else {
1623 tag_data[0] = (id>>16) & 0xff;
1624 tag_data[1] = (id>>8) & 0xff;
1625 tag_data[2] = id & 0xff;
1626 }
1627
1628 m = map_find(c->comp_hdr->tag_encoding_map, tag_data, id);
1629 if (!m)
1630 return -1;
1631 BLOCK_APPEND(s->aux_blk, (char *)tag_data, 3);
1632
1633 if (!m->codec) return -1;
1634 r |= m->codec->decode(s, m->codec, blk, (char *)s->aux_blk, &out_sz);
1635
1636 cr->aux_size += out_sz + 3;
1637 }
1638
1639 return r;
1640 }
1641
1642 static int cram_decode_aux(cram_container *c, cram_slice *s,
1643 cram_block *blk, cram_record *cr) {
1644 int i, r = 0, out_sz = 1;
1645 int32_t TL = 0;
1646 unsigned char *TN;
1647 uint32_t ds = c->comp_hdr->data_series;
1648
1649 if (!(ds & (CRAM_TL|CRAM_aux))) {
1650 cr->aux = 0;
1651 cr->aux_size = 0;
1652 return 0;
1653 }
1654
1655 if (!c->comp_hdr->codecs[DS_TL]) return -1;
1656 r |= c->comp_hdr->codecs[DS_TL]->decode(s, c->comp_hdr->codecs[DS_TL], blk,
1657 (char *)&TL, &out_sz);
1658 if (r || TL < 0 || TL >= c->comp_hdr->nTL)
1659 return -1;
1660
1661 TN = c->comp_hdr->TL[TL];
1662 cr->ntags = strlen((char *)TN)/3; // optimise to remove strlen
1663
1664 //printf("TC=%d\n", cr->ntags);
1665 cr->aux_size = 0;
1666 cr->aux = BLOCK_SIZE(s->aux_blk);
1667
1668 if (!(ds & CRAM_aux))
1669 return 0;
1670
1671 for (i = 0; i < cr->ntags; i++) {
1672 int32_t id, out_sz = 1;
1673 unsigned char tag_data[3];
1674 cram_map *m;
1675
1676 //printf("Tag %d/%d\n", i+1, cr->ntags);
1677 tag_data[0] = *TN++;
1678 tag_data[1] = *TN++;
1679 tag_data[2] = *TN++;
1680 id = (tag_data[0]<<16) | (tag_data[1]<<8) | tag_data[2];
1681
1682 m = map_find(c->comp_hdr->tag_encoding_map, tag_data, id);
1683 if (!m)
1684 return -1;
1685 BLOCK_APPEND(s->aux_blk, (char *)tag_data, 3);
1686
1687 if (!m->codec) return -1;
1688 r |= m->codec->decode(s, m->codec, blk, (char *)s->aux_blk, &out_sz);
1689 cr->aux_size += out_sz + 3;
1690 }
1691
1692 return r;
1693 }
1694
1695 /* Resolve mate pair cross-references between recs within this slice */
1696 static void cram_decode_slice_xref(cram_slice *s, int required_fields) {
1697 int rec;
1698
1699 if (!(required_fields & (SAM_RNEXT | SAM_PNEXT | SAM_TLEN))) {
1700 for (rec = 0; rec < s->hdr->num_records; rec++) {
1701 cram_record *cr = &s->crecs[rec];
1702
1703 cr->tlen = 0;
1704 cr->mate_pos = 0;
1705 cr->mate_ref_id = -1;
1706 }
1707
1708 return;
1709 }
1710
1711 for (rec = 0; rec < s->hdr->num_records; rec++) {
1712 cram_record *cr = &s->crecs[rec];
1713
1714 if (cr->mate_line >= 0) {
1715 if (cr->mate_line < s->hdr->num_records) {
1716 /*
1717 * On the first read, loop through computing lengths.
1718 * It's not perfect as we have one slice per reference so we
1719 * cannot detect when TLEN should be zero due to seqs that
1720 * map to multiple references.
1721 *
1722 * We also cannot set tlen correct when it spans a slice for
1723 * other reasons. This may make tlen too small. Should we
1724 * fix this by forcing TLEN to be stored verbatim in such cases?
1725 *
1726 * Or do we just admit defeat and output 0 for tlen? It's the
1727 * safe option...
1728 */
1729 if (cr->tlen == INT_MIN) {
1730 int id1 = rec, id2 = rec;
1731 int aleft = cr->apos, aright = cr->aend;
1732 int tlen;
1733 int ref = cr->ref_id;
1734
1735 // number of segments starting at the same point.
1736 int left_cnt = 0;
1737
1738 do {
1739 if (aleft > s->crecs[id2].apos)
1740 aleft = s->crecs[id2].apos, left_cnt = 1;
1741 else if (aleft == s->crecs[id2].apos)
1742 left_cnt++;
1743 if (aright < s->crecs[id2].aend)
1744 aright = s->crecs[id2].aend;
1745 if (s->crecs[id2].mate_line == -1) {
1746 s->crecs[id2].mate_line = rec;
1747 break;
1748 }
1749 assert(s->crecs[id2].mate_line > id2);
1750 id2 = s->crecs[id2].mate_line;
1751
1752 if (s->crecs[id2].ref_id != ref)
1753 ref = -1;
1754 } while (id2 != id1);
1755
1756 if (ref != -1) {
1757 tlen = aright - aleft + 1;
1758 id1 = id2 = rec;
1759
1760 /*
1761 * When we have two seqs with identical start and
1762 * end coordinates, set +/- tlen based on 1st/last
1763 * bit flags instead, as a tie breaker.
1764 */
1765 if (s->crecs[id2].apos == aleft) {
1766 if (left_cnt == 1 ||
1767 (s->crecs[id2].flags & BAM_FREAD1))
1768 s->crecs[id2].tlen = tlen;
1769 else
1770 s->crecs[id2].tlen = -tlen;
1771 } else {
1772 s->crecs[id2].tlen = -tlen;
1773 }
1774
1775 id2 = s->crecs[id2].mate_line;
1776 while (id2 != id1) {
1777 if (s->crecs[id2].apos == aleft) {
1778 if (left_cnt == 1 ||
1779 (s->crecs[id2].flags & BAM_FREAD1))
1780 s->crecs[id2].tlen = tlen;
1781 else
1782 s->crecs[id2].tlen = -tlen;
1783 } else {
1784 s->crecs[id2].tlen = -tlen;
1785 }
1786 id2 = s->crecs[id2].mate_line;
1787 }
1788 } else {
1789 id1 = id2 = rec;
1790
1791 s->crecs[id2].tlen = 0;
1792 id2 = s->crecs[id2].mate_line;
1793 while (id2 != id1) {
1794 s->crecs[id2].tlen = 0;
1795 id2 = s->crecs[id2].mate_line;
1796 }
1797 }
1798 }
1799
1800 cr->mate_pos = s->crecs[cr->mate_line].apos;
1801 cr->mate_ref_id = s->crecs[cr->mate_line].ref_id;
1802
1803 // paired
1804 cr->flags |= BAM_FPAIRED;
1805
1806 // set mate unmapped if needed
1807 if (s->crecs[cr->mate_line].flags & BAM_FUNMAP) {
1808 cr->flags |= BAM_FMUNMAP;
1809 cr->tlen = 0;
1810 }
1811 if (cr->flags & BAM_FUNMAP) {
1812 cr->tlen = 0;
1813 }
1814
1815 // set mate reversed if needed
1816 if (s->crecs[cr->mate_line].flags & BAM_FREVERSE)
1817 cr->flags |= BAM_FMREVERSE;
1818 } else {
1819 fprintf(stderr, "Mate line out of bounds: %d vs [0, %d]\n",
1820 cr->mate_line, s->hdr->num_records-1);
1821 }
1822
1823 /* FIXME: construct read names here too if needed */
1824 } else {
1825 if (cr->mate_flags & CRAM_M_REVERSE) {
1826 cr->flags |= BAM_FPAIRED | BAM_FMREVERSE;
1827 }
1828 if (cr->mate_flags & CRAM_M_UNMAP) {
1829 cr->flags |= BAM_FMUNMAP;
1830 //cr->mate_ref_id = -1;
1831 }
1832 if (!(cr->flags & BAM_FPAIRED))
1833 cr->mate_ref_id = -1;
1834 }
1835
1836 if (cr->tlen == INT_MIN)
1837 cr->tlen = 0; // Just incase
1838 }
1839 }
1840
1841 static char *md5_print(unsigned char *md5, char *out) {
1842 int i;
1843 for (i = 0; i < 16; i++) {
1844 out[i*2+0] = "0123456789abcdef"[md5[i]>>4];
1845 out[i*2+1] = "0123456789abcdef"[md5[i]&15];
1846 }
1847 out[32] = 0;
1848
1849 return out;
1850 }
1851
1852 /*
1853 * Decode an entire slice from container blocks. Fills out s->crecs[] array.
1854 * Returns 0 on success
1855 * -1 on failure
1856 */
1857 int cram_decode_slice(cram_fd *fd, cram_container *c, cram_slice *s,
1858 SAM_hdr *bfd) {
1859 cram_block *blk = s->block[0];
1860 int32_t bf, ref_id;
1861 unsigned char cf;
1862 int out_sz, r = 0;
1863 int rec;
1864 char *seq = NULL, *qual = NULL;
1865 int unknown_rg = -1;
1866 int embed_ref;
1867 char **refs = NULL;
1868 uint32_t ds;
1869
1870 if (cram_dependent_data_series(fd, c->comp_hdr, s) != 0)
1871 return -1;
1872
1873 ds = c->comp_hdr->data_series;
1874
1875 blk->bit = 7; // MSB first
1876
1877 /* Look for unknown RG, added as last by Java CRAM? */
1878 if (bfd->nrg > 0 &&
1879 !strcmp(bfd->rg[bfd->nrg-1].name, "UNKNOWN"))
1880 unknown_rg = bfd->nrg-1;
1881
1882 if (blk->content_type != CORE)
1883 return -1;
1884
1885 if (s->crecs)
1886 free(s->crecs);
1887 if (!(s->crecs = malloc(s->hdr->num_records * sizeof(*s->crecs))))
1888 return -1;
1889
1890 ref_id = s->hdr->ref_seq_id;
1891 embed_ref = s->hdr->ref_base_id >= 0 ? 1 : 0;
1892
1893 if (ref_id >= 0) {
1894 if (embed_ref) {
1895 cram_block *b;
1896 if (s->hdr->ref_base_id < 0) {
1897 fprintf(stderr, "No reference specified and "
1898 "no embedded reference is available.\n");
1899 return -1;
1900 }
1901 if (!s->block_by_id ||
1902 !(b = s->block_by_id[s->hdr->ref_base_id]))
1903 return -1;
1904 cram_uncompress_block(b);
1905 s->ref = (char *)BLOCK_DATA(b);
1906 s->ref_start = s->hdr->ref_seq_start;
1907 s->ref_end = s->hdr->ref_seq_start + s->hdr->ref_seq_span-1;
1908 } else if (!fd->no_ref) {
1909 //// Avoid Java cramtools bug by loading entire reference seq
1910 //s->ref = cram_get_ref(fd, s->hdr->ref_seq_id, 1, 0);
1911 //s->ref_start = 1;
1912
1913 if (fd->required_fields & SAM_SEQ)
1914 s->ref =
1915 cram_get_ref(fd, s->hdr->ref_seq_id,
1916 s->hdr->ref_seq_start,
1917 s->hdr->ref_seq_start + s->hdr->ref_seq_span -1);
1918 s->ref_start = s->hdr->ref_seq_start;
1919 s->ref_end = s->hdr->ref_seq_start + s->hdr->ref_seq_span-1;
1920
1921 /* Sanity check */
1922 if (s->ref_start < 0) {
1923 fprintf(stderr, "Slice starts before base 1.\n");
1924 s->ref_start = 0;
1925 }
1926 pthread_mutex_lock(&fd->ref_lock);
1927 pthread_mutex_lock(&fd->refs->lock);
1928 if ((fd->required_fields & SAM_SEQ) &&
1929 s->ref_end > fd->refs->ref_id[ref_id]->length) {
1930 fprintf(stderr, "Slice ends beyond reference end.\n");
1931 s->ref_end = fd->refs->ref_id[ref_id]->length;
1932 }
1933 pthread_mutex_unlock(&fd->refs->lock);
1934 pthread_mutex_unlock(&fd->ref_lock);
1935 }
1936 }
1937
1938 if ((fd->required_fields & SAM_SEQ) &&
1939 s->ref == NULL && s->hdr->ref_seq_id >= 0 && !fd->no_ref) {
1940 fprintf(stderr, "Unable to fetch reference #%d %d..%d\n",
1941 s->hdr->ref_seq_id, s->hdr->ref_seq_start,
1942 s->hdr->ref_seq_start + s->hdr->ref_seq_span-1);
1943 return -1;
1944 }
1945
1946 if (CRAM_MAJOR_VERS(fd->version) != 1
1947 && (fd->required_fields & SAM_SEQ)
1948 && s->hdr->ref_seq_id >= 0
1949 && !fd->ignore_md5
1950 && memcmp(s->hdr->md5, "\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0", 16)) {
1951 MD5_CTX md5;
1952 unsigned char digest[16];
1953
1954 if (s->ref && s->hdr->ref_seq_id >= 0) {
1955 int start, len;
1956
1957 if (s->hdr->ref_seq_start >= s->ref_start) {
1958 start = s->hdr->ref_seq_start - s->ref_start;
1959 } else {
1960 fprintf(stderr, "Slice starts before base 1.\n");
1961 start = 0;
1962 }
1963
1964 if (s->hdr->ref_seq_span <= s->ref_end - s->ref_start + 1) {
1965 len = s->hdr->ref_seq_span;
1966 } else {
1967 fprintf(stderr, "Slice ends beyond reference end.\n");
1968 len = s->ref_end - s->ref_start + 1;
1969 }
1970
1971 MD5_Init(&md5);
1972 if (start + len > s->ref_end - s->ref_start + 1)
1973 len = s->ref_end - s->ref_start + 1 - start;
1974 if (len >= 0)
1975 MD5_Update(&md5, s->ref + start, len);
1976 MD5_Final(digest, &md5);
1977 } else if (!s->ref && s->hdr->ref_base_id >= 0) {
1978 cram_block *b;
1979 if (s->block_by_id && (b = s->block_by_id[s->hdr->ref_base_id])) {
1980 MD5_Init(&md5);
1981 MD5_Update(&md5, b->data, b->uncomp_size);
1982 MD5_Final(digest, &md5);
1983 }
1984 }
1985
1986 if ((!s->ref && s->hdr->ref_base_id < 0)
1987 || memcmp(digest, s->hdr->md5, 16) != 0) {
1988 char M[33];
1989 fprintf(stderr, "ERROR: md5sum reference mismatch for ref "
1990 "%d pos %d..%d\n", ref_id, s->ref_start, s->ref_end);
1991 fprintf(stderr, "CRAM: %s\n", md5_print(s->hdr->md5, M));
1992 fprintf(stderr, "Ref : %s\n", md5_print(digest, M));
1993 return -1;
1994 }
1995 }
1996
1997 if (ref_id == -2) {
1998 pthread_mutex_lock(&fd->ref_lock);
1999 pthread_mutex_lock(&fd->refs->lock);
2000 refs = calloc(fd->refs->nref, sizeof(char *));
2001 pthread_mutex_unlock(&fd->refs->lock);
2002 pthread_mutex_unlock(&fd->ref_lock);
2003 if (!refs)
2004 return -1;
2005 }
2006
2007 for (rec = 0; rec < s->hdr->num_records; rec++) {
2008 cram_record *cr = &s->crecs[rec];
2009
2010 //fprintf(stderr, "Decode seq %d, %d/%d\n", rec, blk->byte, blk->bit);
2011
2012 cr->s = s;
2013
2014 out_sz = 1; /* decode 1 item */
2015 if (ds & CRAM_BF) {
2016 if (!c->comp_hdr->codecs[DS_BF]) return -1;
2017 r |= c->comp_hdr->codecs[DS_BF]
2018 ->decode(s, c->comp_hdr->codecs[DS_BF], blk,
2019 (char *)&bf, &out_sz);
2020 if (bf < 0 ||
2021 bf >= sizeof(fd->bam_flag_swap)/sizeof(*fd->bam_flag_swap))
2022 return -1;
2023 bf = fd->bam_flag_swap[bf];
2024 cr->flags = bf;
2025 } else {
2026 cr->flags = bf = 0x4; // unmapped
2027 }
2028
2029 if (ds & CRAM_CF) {
2030 if (CRAM_MAJOR_VERS(fd->version) == 1) {
2031 /* CF is byte in 1.0, int32 in 2.0 */
2032 if (!c->comp_hdr->codecs[DS_CF]) return -1;
2033 r |= c->comp_hdr->codecs[DS_CF]
2034 ->decode(s, c->comp_hdr->codecs[DS_CF], blk,
2035 (char *)&cf, &out_sz);
2036 cr->cram_flags = cf;
2037 } else {
2038 if (!c->comp_hdr->codecs[DS_CF]) return -1;
2039 r |= c->comp_hdr->codecs[DS_CF]
2040 ->decode(s, c->comp_hdr->codecs[DS_CF], blk,
2041 (char *)&cr->cram_flags,
2042 &out_sz);
2043 cf = cr->cram_flags;
2044 }
2045 }
2046
2047 if (CRAM_MAJOR_VERS(fd->version) != 1 && ref_id == -2) {
2048 if (ds & CRAM_RI) {
2049 if (!c->comp_hdr->codecs[DS_RI]) return -1;
2050 r |= c->comp_hdr->codecs[DS_RI]
2051 ->decode(s, c->comp_hdr->codecs[DS_RI], blk,
2052 (char *)&cr->ref_id, &out_sz);
2053 if ((fd->required_fields & (SAM_SEQ|SAM_TLEN))
2054 && cr->ref_id >= 0) {
2055 if (!fd->no_ref) {
2056 if (!refs[cr->ref_id])
2057 refs[cr->ref_id] = cram_get_ref(fd, cr->ref_id,
2058 1, 0);
2059 s->ref = refs[cr->ref_id];
2060 }
2061 s->ref_start = 1;
2062 pthread_mutex_lock(&fd->ref_lock);
2063 pthread_mutex_lock(&fd->refs->lock);
2064 s->ref_end = fd->refs->ref_id[cr->ref_id]->length;
2065 pthread_mutex_unlock(&fd->refs->lock);
2066 pthread_mutex_unlock(&fd->ref_lock);
2067 }
2068 } else {
2069 cr->ref_id = 0;
2070 }
2071 } else {
2072 cr->ref_id = ref_id; // Forced constant in CRAM 1.0
2073 }
2074
2075
2076 if (ds & CRAM_RL) {
2077 if (!c->comp_hdr->codecs[DS_RL]) return -1;
2078 r |= c->comp_hdr->codecs[DS_RL]
2079 ->decode(s, c->comp_hdr->codecs[DS_RL], blk,
2080 (char *)&cr->len, &out_sz);
2081 }
2082
2083 if (ds & CRAM_AP) {
2084 if (!c->comp_hdr->codecs[DS_AP]) return -1;
2085 r |= c->comp_hdr->codecs[DS_AP]
2086 ->decode(s, c->comp_hdr->codecs[DS_AP], blk,
2087 (char *)&cr->apos, &out_sz);
2088 if (c->comp_hdr->AP_delta)
2089 cr->apos += s->last_apos;
2090 s->last_apos= cr->apos;
2091 } else {
2092 cr->apos = c->ref_seq_start;
2093 }
2094
2095 if (ds & CRAM_RG) {
2096 if (!c->comp_hdr->codecs[DS_RG]) return -1;
2097 r |= c->comp_hdr->codecs[DS_RG]
2098 ->decode(s, c->comp_hdr->codecs[DS_RG], blk,
2099 (char *)&cr->rg, &out_sz);
2100 if (cr->rg == unknown_rg)
2101 cr->rg = -1;
2102 } else {
2103 cr->rg = -1;
2104 }
2105
2106 cr->name_len = 0;
2107
2108 if (c->comp_hdr->read_names_included) {
2109 int32_t out_sz2 = 1;
2110
2111 // Read directly into name cram_block
2112 cr->name = BLOCK_SIZE(s->name_blk);
2113 if (ds & CRAM_RN) {
2114 if (!c->comp_hdr->codecs[DS_RN]) return -1;
2115 r |= c->comp_hdr->codecs[DS_RN]
2116 ->decode(s, c->comp_hdr->codecs[DS_RN], blk,
2117 (char *)s->name_blk, &out_sz2);
2118 cr->name_len = out_sz2;
2119 }
2120 }
2121
2122 cr->mate_pos = 0;
2123 cr->mate_line = -1;
2124 cr->mate_ref_id = -1;
2125 if ((ds & CRAM_CF) && (cf & CRAM_FLAG_DETACHED)) {
2126 if (ds & CRAM_MF) {
2127 if (CRAM_MAJOR_VERS(fd->version) == 1) {
2128 /* MF is byte in 1.0, int32 in 2.0 */
2129 unsigned char mf;
2130 if (!c->comp_hdr->codecs[DS_MF]) return -1;
2131 r |= c->comp_hdr->codecs[DS_MF]
2132 ->decode(s, c->comp_hdr->codecs[DS_MF],
2133 blk, (char *)&mf, &out_sz);
2134 cr->mate_flags = mf;
2135 } else {
2136 if (!c->comp_hdr->codecs[DS_MF]) return -1;
2137 r |= c->comp_hdr->codecs[DS_MF]
2138 ->decode(s, c->comp_hdr->codecs[DS_MF],
2139 blk,
2140 (char *)&cr->mate_flags,
2141 &out_sz);
2142 }
2143 } else {
2144 cr->mate_flags = 0;
2145 }
2146
2147 if (!c->comp_hdr->read_names_included) {
2148 int32_t out_sz2 = 1;
2149
2150 // Read directly into name cram_block
2151 cr->name = BLOCK_SIZE(s->name_blk);
2152 if (ds & CRAM_RN) {
2153 if (!c->comp_hdr->codecs[DS_RN]) return -1;
2154 r |= c->comp_hdr->codecs[DS_RN]
2155 ->decode(s, c->comp_hdr->codecs[DS_RN],
2156 blk, (char *)s->name_blk,
2157 &out_sz2);
2158 cr->name_len = out_sz2;
2159 }
2160 }
2161
2162 if (ds & CRAM_NS) {
2163 if (!c->comp_hdr->codecs[DS_NS]) return -1;
2164 r |= c->comp_hdr->codecs[DS_NS]
2165 ->decode(s, c->comp_hdr->codecs[DS_NS], blk,
2166 (char *)&cr->mate_ref_id, &out_sz);
2167 }
2168
2169 // Skip as mate_ref of "*" is legit. It doesn't mean unmapped, just unknown.
2170 // if (cr->mate_ref_id == -1 && cr->flags & 0x01) {
2171 // /* Paired, but unmapped */
2172 // cr->flags |= BAM_FMUNMAP;
2173 // }
2174
2175 if (ds & CRAM_NP) {
2176 if (!c->comp_hdr->codecs[DS_NP]) return -1;
2177 r |= c->comp_hdr->codecs[DS_NP]
2178 ->decode(s, c->comp_hdr->codecs[DS_NP], blk,
2179 (char *)&cr->mate_pos, &out_sz);
2180 }
2181
2182 if (ds & CRAM_TS) {
2183 if (!c->comp_hdr->codecs[DS_TS]) return -1;
2184 r |= c->comp_hdr->codecs[DS_TS]
2185 ->decode(s, c->comp_hdr->codecs[DS_TS], blk,
2186 (char *)&cr->tlen, &out_sz);
2187 } else {
2188 cr->tlen = INT_MIN;
2189 }
2190 } else if ((ds & CRAM_CF) && (cf & CRAM_FLAG_MATE_DOWNSTREAM)) {
2191 if (ds & CRAM_NF) {
2192 if (!c->comp_hdr->codecs[DS_NF]) return -1;
2193 r |= c->comp_hdr->codecs[DS_NF]
2194 ->decode(s, c->comp_hdr->codecs[DS_NF], blk,
2195 (char *)&cr->mate_line, &out_sz);
2196 cr->mate_line += rec + 1;
2197
2198 //cr->name_len = sprintf(name, "%d", name_id++);
2199 //cr->name = DSTRING_LEN(name_ds);
2200 //dstring_nappend(name_ds, name, cr->name_len);
2201
2202 cr->mate_ref_id = -1;
2203 cr->tlen = INT_MIN;
2204 cr->mate_pos = 0;
2205 } else {
2206 cr->mate_flags = 0;
2207 cr->tlen = INT_MIN;
2208 }
2209 } else {
2210 cr->mate_flags = 0;
2211 cr->tlen = INT_MIN;
2212 }
2213 /*
2214 else if (!name[0]) {
2215 //name[0] = '?'; name[1] = 0;
2216 //cr->name_len = 1;
2217 //cr->name= DSTRING_LEN(s->name_ds);
2218 //dstring_nappend(s->name_ds, "?", 1);
2219
2220 cr->mate_ref_id = -1;
2221 cr->tlen = 0;
2222 cr->mate_pos = 0;
2223 }
2224 */
2225
2226 /* Auxiliary tags */
2227 if (CRAM_MAJOR_VERS(fd->version) == 1)
2228 r |= cram_decode_aux_1_0(c, s, blk, cr);
2229 else
2230 r |= cram_decode_aux(c, s, blk, cr);
2231
2232 /* Fake up dynamic string growth and appending */
2233 if (ds & CRAM_RL) {
2234 cr->seq = BLOCK_SIZE(s->seqs_blk);
2235 BLOCK_GROW(s->seqs_blk, cr->len);
2236 seq = (char *)BLOCK_END(s->seqs_blk);
2237 BLOCK_SIZE(s->seqs_blk) += cr->len;
2238
2239 if (!seq)
2240 return -1;
2241
2242 cr->qual = BLOCK_SIZE(s->qual_blk);
2243 BLOCK_GROW(s->qual_blk, cr->len);
2244 qual = (char *)BLOCK_END(s->qual_blk);
2245 BLOCK_SIZE(s->qual_blk) += cr->len;
2246
2247 if (!s->ref)
2248 memset(seq, '=', cr->len);
2249 }
2250
2251 if (!(bf & BAM_FUNMAP)) {
2252 /* Decode sequence and generate CIGAR */
2253 if (ds & (CRAM_SEQ | CRAM_MQ)) {
2254 r |= cram_decode_seq(fd, c, s, blk, cr, bfd, cf, seq, qual);
2255 } else {
2256 cr->cigar = 0;
2257 cr->ncigar = 0;
2258 cr->aend = cr->apos;
2259 cr->mqual = 0;
2260 }
2261 } else {
2262 int out_sz2 = cr->len;
2263
2264 //puts("Unmapped");
2265 cr->cigar = 0;
2266 cr->ncigar = 0;
2267 cr->aend = cr->apos;
2268 cr->mqual = 0;
2269
2270 if (ds & CRAM_BA) {
2271 if (!c->comp_hdr->codecs[DS_BA]) return -1;
2272 r |= c->comp_hdr->codecs[DS_BA]
2273 ->decode(s, c->comp_hdr->codecs[DS_BA], blk,
2274 (char *)seq, &out_sz2);
2275 }
2276
2277 if ((ds & CRAM_CF) && (cf & CRAM_FLAG_PRESERVE_QUAL_SCORES)) {
2278 out_sz2 = cr->len;
2279 if (ds & CRAM_QS) {
2280 if (!c->comp_hdr->codecs[DS_QS]) return -1;
2281 r |= c->comp_hdr->codecs[DS_QS]
2282 ->decode(s, c->comp_hdr->codecs[DS_QS],
2283 blk, qual, &out_sz2);
2284 }
2285 } else {
2286 if (ds & CRAM_RL)
2287 memset(qual, 30, cr->len);
2288 }
2289 }
2290 }
2291
2292 pthread_mutex_lock(&fd->ref_lock);
2293 if (refs) {
2294 int i;
2295 for (i = 0; i < fd->refs->nref; i++) {
2296 if (refs[i])
2297 cram_ref_decr(fd->refs, i);
2298 }
2299 free(refs);
2300 } else if (ref_id >= 0 && s->ref != fd->ref_free) {
2301 cram_ref_decr(fd->refs, ref_id);
2302 }
2303 pthread_mutex_unlock(&fd->ref_lock);
2304
2305 /* Resolve mate pair cross-references between recs within this slice */
2306 cram_decode_slice_xref(s, fd->required_fields);
2307
2308 return r;
2309 }
2310
2311 typedef struct {
2312 cram_fd *fd;
2313 cram_container *c;
2314 cram_slice *s;
2315 SAM_hdr *h;
2316 int exit_code;
2317 } cram_decode_job;
2318
2319 void *cram_decode_slice_thread(void *arg) {
2320 cram_decode_job *j = (cram_decode_job *)arg;
2321
2322 j->exit_code = cram_decode_slice(j->fd, j->c, j->s, j->h);
2323
2324 return j;
2325 }
2326
2327 /*
2328 * Spawn a multi-threaded version of cram_decode_slice().
2329 */
2330 int cram_decode_slice_mt(cram_fd *fd, cram_container *c, cram_slice *s,
2331 SAM_hdr *bfd) {
2332 cram_decode_job *j;
2333 int nonblock;
2334
2335 if (!fd->pool)
2336 return cram_decode_slice(fd, c, s, bfd);
2337
2338 if (!(j = malloc(sizeof(*j))))
2339 return -1;
2340
2341 j->fd = fd;
2342 j->c = c;
2343 j->s = s;
2344 j->h = bfd;
2345
2346 nonblock = t_pool_results_queue_sz(fd->rqueue) ? 1 : 0;
2347
2348 if (-1 == t_pool_dispatch2(fd->pool, fd->rqueue, cram_decode_slice_thread,
2349 j, nonblock)) {
2350 /* Would block */
2351 fd->job_pending = j;
2352 } else {
2353 fd->job_pending = NULL;
2354 }
2355
2356 // flush too
2357 return 0;
2358 }
2359
2360
2361 /* ----------------------------------------------------------------------
2362 * CRAM sequence iterators.
2363 */
2364
2365 /*
2366 * Converts a cram in-memory record into a bam in-memory record. We
2367 * pass a pointer to a bam_seq_t pointer along with the a pointer to
2368 * the allocated size. These can initially be pointers to NULL and zero.
2369 *
2370 * This function will reallocate the bam buffer as required and update
2371 * (*bam)->alloc accordingly, allowing it to be used within a loop
2372 * efficiently without needing to allocate new bam objects over and
2373 * over again.
2374 *
2375 * Returns the used size of the bam record on success
2376 * -1 on failure.
2377 */
2378 static int cram_to_bam(SAM_hdr *bfd, cram_fd *fd, cram_slice *s,
2379 cram_record *cr, int rec, bam_seq_t **bam) {
2380 int bam_idx, rg_len;
2381 char name_a[1024], *name;
2382 int name_len;
2383 char *aux, *aux_orig;
2384 char *seq, *qual;
2385
2386 /* Assign names if not explicitly set */
2387 if (fd->required_fields & SAM_QNAME) {
2388 if (cr->name_len) {
2389 name = (char *)BLOCK_DATA(s->name_blk) + cr->name;
2390 name_len = cr->name_len;
2391 } else {
2392 name = name_a;
2393 name_len = strlen(fd->prefix);
2394 memcpy(name, fd->prefix, name_len);
2395 name += name_len;
2396 *name++ = ':';
2397 if (cr->mate_line >= 0 && cr->mate_line < rec)
2398 name = (char *)append_uint64((unsigned char *)name,
2399 s->hdr->record_counter +
2400 cr->mate_line + 1);
2401 else
2402 name = (char *)append_uint64((unsigned char *)name,
2403 s->hdr->record_counter +
2404 rec + 1);
2405 name_len = name - name_a;
2406 name = name_a;
2407 }
2408 } else {
2409 name = "?";
2410 name_len = 1;
2411 }
2412
2413 /* Generate BAM record */
2414 if (cr->rg < -1 || cr->rg >= bfd->nrg)
2415 return -1;
2416 rg_len = (cr->rg != -1) ? bfd->rg[cr->rg].name_len + 4 : 0;
2417
2418 if (fd->required_fields & (SAM_SEQ | SAM_QUAL)) {
2419 if (!BLOCK_DATA(s->seqs_blk))
2420 return -1;
2421 seq = (char *)BLOCK_DATA(s->seqs_blk) + cr->seq;
2422 } else {
2423 seq = "*";
2424 cr->len = 1;
2425 }
2426
2427
2428 if (fd->required_fields & SAM_QUAL) {
2429 if (!BLOCK_DATA(s->qual_blk))
2430 return -1;
2431 qual = (char *)BLOCK_DATA(s->qual_blk) + cr->qual;
2432 } else {
2433 qual = NULL;
2434 }
2435
2436 bam_idx = bam_construct_seq(bam, cr->aux_size + rg_len,
2437 name, name_len,
2438 cr->flags,
2439 cr->ref_id,
2440 cr->apos,
2441 cr->aend,
2442 cr->mqual,
2443 cr->ncigar, &s->cigar[cr->cigar],
2444 cr->mate_ref_id,
2445 cr->mate_pos,
2446 cr->tlen,
2447 cr->len,
2448 seq,
2449 qual);
2450 if (bam_idx == -1)
2451 return -1;
2452
2453 aux = aux_orig = (char *)bam_aux(*bam);
2454
2455 /* Auxiliary strings */
2456 if (cr->aux_size != 0) {
2457 memcpy(aux, BLOCK_DATA(s->aux_blk) + cr->aux, cr->aux_size);
2458 aux += cr->aux_size;
2459 }
2460
2461 /* RG:Z: */
2462 if (cr->rg != -1) {
2463 int len = bfd->rg[cr->rg].name_len;
2464 *aux++ = 'R'; *aux++ = 'G'; *aux++ = 'Z';
2465 memcpy(aux, bfd->rg[cr->rg].name, len);
2466 aux += len;
2467 *aux++ = 0;
2468 }
2469
2470 return bam_idx + (aux - aux_orig);
2471 }
2472
2473 /*
2474 * Here be dragons! The multi-threading code in this is crufty beyond belief.
2475 */
2476 static cram_slice *cram_next_slice(cram_fd *fd, cram_container **cp) {
2477 cram_container *c;
2478 cram_slice *s = NULL;
2479
2480 if (!(c = fd->ctr)) {
2481 // Load first container.
2482 do {
2483 if (!(c = fd->ctr = cram_read_container(fd)))
2484 return NULL;
2485 } while (c->length == 0);
2486
2487 /*
2488 * The first container may be a result of a sub-range query.
2489 * In which case it may still not be the optimal starting point
2490 * due to skipped containers/slices in the index.
2491 */
2492 if (fd->range.refid != -2) {
2493 while (c->ref_seq_id != -2 &&
2494 (c->ref_seq_id < fd->range.refid ||
2495 c->ref_seq_start + c->ref_seq_span-1 < fd->range.start)) {
2496 if (0 != cram_seek(fd, c->length, SEEK_CUR))
2497 return NULL;
2498 cram_free_container(fd->ctr);
2499 do {
2500 if (!(c = fd->ctr = cram_read_container(fd)))
2501 return NULL;
2502 } while (c->length == 0);
2503 }
2504
2505 if (c->ref_seq_id != -2 && c->ref_seq_id != fd->range.refid)
2506 return NULL;
2507 }
2508
2509 if (!(c->comp_hdr_block = cram_read_block(fd)))
2510 return NULL;
2511 if (c->comp_hdr_block->content_type != COMPRESSION_HEADER)
2512 return NULL;
2513
2514 c->comp_hdr = cram_decode_compression_header(fd, c->comp_hdr_block);
2515 if (!c->comp_hdr)
2516 return NULL;
2517 if (!c->comp_hdr->AP_delta) {
2518 pthread_mutex_lock(&fd->ref_lock);
2519 fd->unsorted = 1;
2520 pthread_mutex_unlock(&fd->ref_lock);
2521 }
2522 }
2523
2524 if ((s = c->slice)) {
2525 c->slice = NULL;
2526 cram_free_slice(s);
2527 s = NULL;
2528 }
2529
2530 if (c->curr_slice == c->max_slice) {
2531 cram_free_container(c);
2532 c = NULL;
2533 }
2534
2535 /* Sorry this is so contorted! */
2536 for (;;) {
2537 if (fd->job_pending) {
2538 cram_decode_job *j = (cram_decode_job *)fd->job_pending;
2539 c = j->c;
2540 s = j->s;
2541 free(fd->job_pending);
2542 fd->job_pending = NULL;
2543 } else if (!fd->ooc) {
2544 empty_container:
2545 if (!c || c->curr_slice == c->max_slice) {
2546 // new container
2547 do {
2548 if (!(c = fd->ctr = cram_read_container(fd))) {
2549 if (fd->pool) {
2550 fd->ooc = 1;
2551 break;
2552 }
2553
2554 return NULL;
2555 }
2556 } while (c->length == 0);
2557 if (fd->ooc)
2558 break;
2559
2560 /* Skip containers not yet spanning our range */
2561 if (fd->range.refid != -2 && c->ref_seq_id != -2) {
2562 fd->required_fields |= SAM_POS;
2563
2564 if (c->ref_seq_id != fd->range.refid) {
2565 cram_free_container(c);
2566 fd->ctr = NULL;
2567 fd->ooc = 1;
2568 fd->eof = 1;
2569 break;
2570 }
2571
2572 if (c->ref_seq_start > fd->range.end) {
2573 cram_free_container(c);
2574 fd->ctr = NULL;
2575 fd->ooc = 1;
2576 fd->eof = 1;
2577 break;
2578 }
2579
2580 if (c->ref_seq_start + c->ref_seq_span-1 <
2581 fd->range.start) {
2582 c->curr_rec = c->max_rec;
2583 c->curr_slice = c->max_slice;
2584 cram_seek(fd, c->length, SEEK_CUR);
2585 cram_free_container(c);
2586 c = NULL;
2587 continue;
2588 }
2589 }
2590
2591 if (!(c->comp_hdr_block = cram_read_block(fd)))
2592 return NULL;
2593 if (c->comp_hdr_block->content_type != COMPRESSION_HEADER)
2594 return NULL;
2595
2596 c->comp_hdr =
2597 cram_decode_compression_header(fd, c->comp_hdr_block);
2598 if (!c->comp_hdr)
2599 return NULL;
2600
2601 if (!c->comp_hdr->AP_delta) {
2602 pthread_mutex_lock(&fd->ref_lock);
2603 fd->unsorted = 1;
2604 pthread_mutex_unlock(&fd->ref_lock);
2605 }
2606 }
2607
2608 if (c->num_records == 0) {
2609 cram_free_container(c); c = NULL;
2610 goto empty_container;
2611 }
2612
2613
2614 if (!(s = c->slice = cram_read_slice(fd)))
2615 return NULL;
2616 c->curr_slice++;
2617 c->curr_rec = 0;
2618 c->max_rec = s->hdr->num_records;
2619
2620 s->last_apos = s->hdr->ref_seq_start;
2621
2622 /* Skip slices not yet spanning our range */
2623 if (fd->range.refid != -2 && s->hdr->ref_seq_id != -2) {
2624 if (s->hdr->ref_seq_id != fd->range.refid) {
2625 fd->eof = 1;
2626 cram_free_slice(s);
2627 c->slice = NULL;
2628 return NULL;
2629 }
2630
2631 if (s->hdr->ref_seq_start > fd->range.end) {
2632 fd->eof = 1;
2633 cram_free_slice(s);
2634 c->slice = NULL;
2635 return NULL;
2636 }
2637
2638 if (s->hdr->ref_seq_start + s->hdr->ref_seq_span-1 <
2639 fd->range.start) {
2640 cram_free_slice(s);
2641 c->slice = NULL;
2642 cram_free_container(c);
2643 c = NULL;
2644 continue;
2645 }
2646 }
2647 }
2648
2649 /* Test decoding of 1st seq */
2650 if (!c || !s)
2651 break;
2652
2653 if (cram_decode_slice_mt(fd, c, s, fd->header) != 0) {
2654 // if (cram_decode_slice(fd, c, s, fd->header) != 0) {
2655 fprintf(stderr, "Failure to decode slice\n");
2656 cram_free_slice(s);
2657 c->slice = NULL;
2658 return NULL;
2659 }
2660
2661 if (!fd->pool || fd->job_pending)
2662 break;
2663
2664 // Push it a bit far, to qsize in queue rather than pending arrival,
2665 // as cram tends to be a bit bursty in decode timings.
2666 if (t_pool_results_queue_len(fd->rqueue) > fd->pool->qsize)
2667 break;
2668 }
2669
2670 if (fd->pool) {
2671 t_pool_result *res;
2672 cram_decode_job *j;
2673
2674 // fprintf(stderr, "Thread pool len = %d, %d\n",
2675 // t_pool_results_queue_len(fd->rqueue),
2676 // t_pool_results_queue_sz(fd->rqueue));
2677
2678 if (fd->ooc && t_pool_results_queue_empty(fd->rqueue))
2679 return NULL;
2680
2681 res = t_pool_next_result_wait(fd->rqueue);
2682
2683 if (!res || !res->data) {
2684 fprintf(stderr, "t_pool_next_result failure\n");
2685 return NULL;
2686 }
2687
2688 j = (cram_decode_job *)res->data;
2689 c = j->c;
2690 s = j->s;
2691
2692 fd->ctr = c;
2693
2694 t_pool_delete_result(res, 1);
2695 }
2696
2697 *cp = c;
2698 return s;
2699 }
2700
2701 /*
2702 * Read the next cram record and return it.
2703 * Note that to decode cram_record the caller will need to look up some data
2704 * in the current slice, pointed to by fd->ctr->slice. This is valid until
2705 * the next call to cram_get_seq (which may invalidate it).
2706 *
2707 * Returns record pointer on success (do not free)
2708 * NULL on failure
2709 */
2710 cram_record *cram_get_seq(cram_fd *fd) {
2711 cram_container *c;
2712 cram_slice *s;
2713
2714 for (;;) {
2715 c = fd->ctr;
2716 if (c && c->slice && c->curr_rec < c->max_rec) {
2717 s = c->slice;
2718 } else {
2719 if (!(s = cram_next_slice(fd, &c)))
2720 return NULL;
2721 }
2722
2723 if (fd->range.refid != -2) {
2724 if (s->crecs[c->curr_rec].ref_id < fd->range.refid) {
2725 c->curr_rec++;
2726 continue;
2727 }
2728
2729 if (s->crecs[c->curr_rec].ref_id != fd->range.refid) {
2730 fd->eof = 1;
2731 cram_free_slice(s);
2732 c->slice = NULL;
2733 return NULL;
2734 }
2735
2736 if (s->crecs[c->curr_rec].apos > fd->range.end) {
2737 fd->eof = 1;
2738 cram_free_slice(s);
2739 c->slice = NULL;
2740 return NULL;
2741 }
2742
2743 if (s->crecs[c->curr_rec].aend < fd->range.start) {
2744 c->curr_rec++;
2745 continue;
2746 }
2747 }
2748
2749 break;
2750 }
2751
2752 fd->ctr = c;
2753 c->slice = s;
2754 return &s->crecs[c->curr_rec++];
2755 }
2756
2757 /*
2758 * Read the next cram record and convert it to a bam_seq_t struct.
2759 *
2760 * Returns 0 on success
2761 * -1 on EOF or failure (check fd->err)
2762 */
2763 int cram_get_bam_seq(cram_fd *fd, bam_seq_t **bam) {
2764 cram_record *cr;
2765 cram_container *c;
2766 cram_slice *s;
2767
2768 if (!(cr = cram_get_seq(fd)))
2769 return -1;
2770
2771 c = fd->ctr;
2772 s = c->slice;
2773
2774 return cram_to_bam(fd->header, fd, s, cr, c->curr_rec-1, bam);
2775 }