0
|
1 #include <stdlib.h>
|
|
2 #include <ctype.h>
|
|
3 #include <string.h>
|
|
4 #include <stdio.h>
|
|
5 #include "bamlite.h"
|
|
6
|
|
7 /*********************
|
|
8 * from bam_endian.c *
|
|
9 *********************/
|
|
10
|
|
11 static inline int bam_is_big_endian()
|
|
12 {
|
|
13 long one= 1;
|
|
14 return !(*((char *)(&one)));
|
|
15 }
|
|
16 static inline uint16_t bam_swap_endian_2(uint16_t v)
|
|
17 {
|
|
18 return (uint16_t)(((v & 0x00FF00FFU) << 8) | ((v & 0xFF00FF00U) >> 8));
|
|
19 }
|
|
20 static inline void *bam_swap_endian_2p(void *x)
|
|
21 {
|
|
22 *(uint16_t*)x = bam_swap_endian_2(*(uint16_t*)x);
|
|
23 return x;
|
|
24 }
|
|
25 static inline uint32_t bam_swap_endian_4(uint32_t v)
|
|
26 {
|
|
27 v = ((v & 0x0000FFFFU) << 16) | (v >> 16);
|
|
28 return ((v & 0x00FF00FFU) << 8) | ((v & 0xFF00FF00U) >> 8);
|
|
29 }
|
|
30 static inline void *bam_swap_endian_4p(void *x)
|
|
31 {
|
|
32 *(uint32_t*)x = bam_swap_endian_4(*(uint32_t*)x);
|
|
33 return x;
|
|
34 }
|
|
35 static inline uint64_t bam_swap_endian_8(uint64_t v)
|
|
36 {
|
|
37 v = ((v & 0x00000000FFFFFFFFLLU) << 32) | (v >> 32);
|
|
38 v = ((v & 0x0000FFFF0000FFFFLLU) << 16) | ((v & 0xFFFF0000FFFF0000LLU) >> 16);
|
|
39 return ((v & 0x00FF00FF00FF00FFLLU) << 8) | ((v & 0xFF00FF00FF00FF00LLU) >> 8);
|
|
40 }
|
|
41 static inline void *bam_swap_endian_8p(void *x)
|
|
42 {
|
|
43 *(uint64_t*)x = bam_swap_endian_8(*(uint64_t*)x);
|
|
44 return x;
|
|
45 }
|
|
46
|
|
47 /**************
|
|
48 * from bam.c *
|
|
49 **************/
|
|
50
|
|
51 int bam_is_be;
|
|
52
|
|
53 bam_header_t *bam_header_init()
|
|
54 {
|
|
55 bam_is_be = bam_is_big_endian();
|
|
56 return (bam_header_t*)calloc(1, sizeof(bam_header_t));
|
|
57 }
|
|
58
|
|
59 void bam_header_destroy(bam_header_t *header)
|
|
60 {
|
|
61 int32_t i;
|
|
62 if (header == 0) return;
|
|
63 if (header->target_name) {
|
|
64 for (i = 0; i < header->n_targets; ++i)
|
|
65 free(header->target_name[i]);
|
|
66 free(header->target_name);
|
|
67 free(header->target_len);
|
|
68 }
|
|
69 free(header->text);
|
|
70 free(header);
|
|
71 }
|
|
72
|
|
73 bam_header_t *bam_header_read(bamFile fp)
|
|
74 {
|
|
75 bam_header_t *header;
|
|
76 char buf[4];
|
|
77 int magic_len;
|
|
78 int32_t i = 1, name_len;
|
|
79 // read "BAM1"
|
|
80 magic_len = bam_read(fp, buf, 4);
|
|
81 if (magic_len != 4 || strncmp(buf, "BAM\001", 4) != 0) {
|
|
82 fprintf(stderr, "[bam_header_read] invalid BAM binary header (this is not a BAM file).\n");
|
|
83 return 0;
|
|
84 }
|
|
85 header = bam_header_init();
|
|
86 // read plain text and the number of reference sequences
|
|
87 bam_read(fp, &header->l_text, 4);
|
|
88 if (bam_is_be) bam_swap_endian_4p(&header->l_text);
|
|
89 header->text = (char*)calloc(header->l_text + 1, 1);
|
|
90 bam_read(fp, header->text, header->l_text);
|
|
91 bam_read(fp, &header->n_targets, 4);
|
|
92 if (bam_is_be) bam_swap_endian_4p(&header->n_targets);
|
|
93 // read reference sequence names and lengths
|
|
94 header->target_name = (char**)calloc(header->n_targets, sizeof(char*));
|
|
95 header->target_len = (uint32_t*)calloc(header->n_targets, 4);
|
|
96 for (i = 0; i != header->n_targets; ++i) {
|
|
97 bam_read(fp, &name_len, 4);
|
|
98 if (bam_is_be) bam_swap_endian_4p(&name_len);
|
|
99 header->target_name[i] = (char*)calloc(name_len, 1);
|
|
100 bam_read(fp, header->target_name[i], name_len);
|
|
101 bam_read(fp, &header->target_len[i], 4);
|
|
102 if (bam_is_be) bam_swap_endian_4p(&header->target_len[i]);
|
|
103 }
|
|
104 return header;
|
|
105 }
|
|
106
|
|
107 static void swap_endian_data(const bam1_core_t *c, int data_len, uint8_t *data)
|
|
108 {
|
|
109 uint8_t *s;
|
|
110 uint32_t i, *cigar = (uint32_t*)(data + c->l_qname);
|
|
111 s = data + c->n_cigar*4 + c->l_qname + c->l_qseq + (c->l_qseq + 1)/2;
|
|
112 for (i = 0; i < c->n_cigar; ++i) bam_swap_endian_4p(&cigar[i]);
|
|
113 while (s < data + data_len) {
|
|
114 uint8_t type;
|
|
115 s += 2; // skip key
|
|
116 type = toupper(*s); ++s; // skip type
|
|
117 if (type == 'C' || type == 'A') ++s;
|
|
118 else if (type == 'S') { bam_swap_endian_2p(s); s += 2; }
|
|
119 else if (type == 'I' || type == 'F') { bam_swap_endian_4p(s); s += 4; }
|
|
120 else if (type == 'D') { bam_swap_endian_8p(s); s += 8; }
|
|
121 else if (type == 'Z' || type == 'H') { while (*s) ++s; ++s; }
|
|
122 }
|
|
123 }
|
|
124
|
|
125 int bam_read1(bamFile fp, bam1_t *b)
|
|
126 {
|
|
127 bam1_core_t *c = &b->core;
|
|
128 int32_t block_len, ret, i;
|
|
129 uint32_t x[8];
|
|
130
|
|
131 if ((ret = bam_read(fp, &block_len, 4)) != 4) {
|
|
132 if (ret == 0) return -1; // normal end-of-file
|
|
133 else return -2; // truncated
|
|
134 }
|
|
135 if (bam_read(fp, x, sizeof(bam1_core_t)) != sizeof(bam1_core_t)) return -3;
|
|
136 if (bam_is_be) {
|
|
137 bam_swap_endian_4p(&block_len);
|
|
138 for (i = 0; i < 8; ++i) bam_swap_endian_4p(x + i);
|
|
139 }
|
|
140 c->tid = x[0]; c->pos = x[1];
|
|
141 c->bin = x[2]>>16; c->qual = x[2]>>8&0xff; c->l_qname = x[2]&0xff;
|
|
142 c->flag = x[3]>>16; c->n_cigar = x[3]&0xffff;
|
|
143 c->l_qseq = x[4];
|
|
144 c->mtid = x[5]; c->mpos = x[6]; c->isize = x[7];
|
|
145 b->data_len = block_len - sizeof(bam1_core_t);
|
|
146 if (b->m_data < b->data_len) {
|
|
147 b->m_data = b->data_len;
|
|
148 kroundup32(b->m_data);
|
|
149 b->data = (uint8_t*)realloc(b->data, b->m_data);
|
|
150 }
|
|
151 if (bam_read(fp, b->data, b->data_len) != b->data_len) return -4;
|
|
152 b->l_aux = b->data_len - c->n_cigar * 4 - c->l_qname - c->l_qseq - (c->l_qseq+1)/2;
|
|
153 if (bam_is_be) swap_endian_data(c, b->data_len, b->data);
|
|
154 return 4 + block_len;
|
|
155 }
|