|
0
|
1 /* sam.c -- SAM and BAM file I/O and manipulation.
|
|
|
2
|
|
|
3 Copyright (C) 2008-2010, 2012-2014 Genome Research Ltd.
|
|
|
4 Copyright (C) 2010, 2012, 2013 Broad Institute.
|
|
|
5
|
|
|
6 Author: Heng Li <lh3@sanger.ac.uk>
|
|
|
7
|
|
|
8 Permission is hereby granted, free of charge, to any person obtaining a copy
|
|
|
9 of this software and associated documentation files (the "Software"), to deal
|
|
|
10 in the Software without restriction, including without limitation the rights
|
|
|
11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|
|
12 copies of the Software, and to permit persons to whom the Software is
|
|
|
13 furnished to do so, subject to the following conditions:
|
|
|
14
|
|
|
15 The above copyright notice and this permission notice shall be included in
|
|
|
16 all copies or substantial portions of the Software.
|
|
|
17
|
|
|
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
|
19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
|
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
|
|
|
21 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
|
22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|
|
23 FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
|
|
|
24 DEALINGS IN THE SOFTWARE. */
|
|
|
25
|
|
|
26 #include <stdio.h>
|
|
|
27 #include <stdlib.h>
|
|
|
28 #include <string.h>
|
|
|
29 #include <errno.h>
|
|
|
30 #include <ctype.h>
|
|
|
31 #include <zlib.h>
|
|
|
32 #include "htslib/sam.h"
|
|
|
33 #include "htslib/bgzf.h"
|
|
|
34 #include "cram/cram.h"
|
|
|
35 #include "htslib/hfile.h"
|
|
|
36
|
|
|
37 #include "htslib/khash.h"
|
|
|
38 KHASH_DECLARE(s2i, kh_cstr_t, int64_t)
|
|
|
39
|
|
|
40 typedef khash_t(s2i) sdict_t;
|
|
|
41
|
|
|
42 /**********************
|
|
|
43 *** BAM header I/O ***
|
|
|
44 **********************/
|
|
|
45
|
|
|
46 bam_hdr_t *bam_hdr_init()
|
|
|
47 {
|
|
|
48 return (bam_hdr_t*)calloc(1, sizeof(bam_hdr_t));
|
|
|
49 }
|
|
|
50
|
|
|
51 void bam_hdr_destroy(bam_hdr_t *h)
|
|
|
52 {
|
|
|
53 int32_t i;
|
|
|
54 if (h == NULL) return;
|
|
|
55 if (h->target_name) {
|
|
|
56 for (i = 0; i < h->n_targets; ++i)
|
|
|
57 free(h->target_name[i]);
|
|
|
58 free(h->target_name);
|
|
|
59 free(h->target_len);
|
|
|
60 }
|
|
|
61 free(h->text); free(h->cigar_tab);
|
|
|
62 if (h->sdict) kh_destroy(s2i, (sdict_t*)h->sdict);
|
|
|
63 free(h);
|
|
|
64 }
|
|
|
65
|
|
|
66 bam_hdr_t *bam_hdr_dup(const bam_hdr_t *h0)
|
|
|
67 {
|
|
|
68 if (h0 == NULL) return NULL;
|
|
|
69 bam_hdr_t *h;
|
|
|
70 if ((h = bam_hdr_init()) == NULL) return NULL;
|
|
|
71 // copy the simple data
|
|
|
72 h->n_targets = h0->n_targets;
|
|
|
73 h->ignore_sam_err = h0->ignore_sam_err;
|
|
|
74 h->l_text = h0->l_text;
|
|
|
75 // Then the pointery stuff
|
|
|
76 h->cigar_tab = NULL;
|
|
|
77 h->sdict = NULL;
|
|
|
78 h->text = (char*)calloc(h->l_text + 1, 1);
|
|
|
79 memcpy(h->text, h0->text, h->l_text);
|
|
|
80 h->target_len = (uint32_t*)calloc(h->n_targets, sizeof(uint32_t));
|
|
|
81 h->target_name = (char**)calloc(h->n_targets, sizeof(char*));
|
|
|
82 int i;
|
|
|
83 for (i = 0; i < h->n_targets; ++i) {
|
|
|
84 h->target_len[i] = h0->target_len[i];
|
|
|
85 h->target_name[i] = strdup(h0->target_name[i]);
|
|
|
86 }
|
|
|
87 return h;
|
|
|
88 }
|
|
|
89
|
|
|
90
|
|
|
91 static bam_hdr_t *hdr_from_dict(sdict_t *d)
|
|
|
92 {
|
|
|
93 bam_hdr_t *h;
|
|
|
94 khint_t k;
|
|
|
95 h = bam_hdr_init();
|
|
|
96 h->sdict = d;
|
|
|
97 h->n_targets = kh_size(d);
|
|
|
98 h->target_len = (uint32_t*)malloc(sizeof(uint32_t) * h->n_targets);
|
|
|
99 h->target_name = (char**)malloc(sizeof(char*) * h->n_targets);
|
|
|
100 for (k = kh_begin(d); k != kh_end(d); ++k) {
|
|
|
101 if (!kh_exist(d, k)) continue;
|
|
|
102 h->target_name[kh_val(d, k)>>32] = (char*)kh_key(d, k);
|
|
|
103 h->target_len[kh_val(d, k)>>32] = kh_val(d, k)<<32>>32;
|
|
|
104 kh_val(d, k) >>= 32;
|
|
|
105 }
|
|
|
106 return h;
|
|
|
107 }
|
|
|
108
|
|
|
109 bam_hdr_t *bam_hdr_read(BGZF *fp)
|
|
|
110 {
|
|
|
111 bam_hdr_t *h;
|
|
|
112 char buf[4];
|
|
|
113 int magic_len, has_EOF;
|
|
|
114 int32_t i = 1, name_len;
|
|
|
115 // check EOF
|
|
|
116 has_EOF = bgzf_check_EOF(fp);
|
|
|
117 if (has_EOF < 0) {
|
|
|
118 perror("[W::sam_hdr_read] bgzf_check_EOF");
|
|
|
119 } else if (has_EOF == 0 && hts_verbose >= 2)
|
|
|
120 fprintf(stderr, "[W::%s] EOF marker is absent. The input is probably truncated.\n", __func__);
|
|
|
121 // read "BAM1"
|
|
|
122 magic_len = bgzf_read(fp, buf, 4);
|
|
|
123 if (magic_len != 4 || strncmp(buf, "BAM\1", 4)) {
|
|
|
124 if (hts_verbose >= 1) fprintf(stderr, "[E::%s] invalid BAM binary header\n", __func__);
|
|
|
125 return 0;
|
|
|
126 }
|
|
|
127 h = bam_hdr_init();
|
|
|
128 // read plain text and the number of reference sequences
|
|
|
129 bgzf_read(fp, &h->l_text, 4);
|
|
|
130 if (fp->is_be) ed_swap_4p(&h->l_text);
|
|
|
131 h->text = (char*)malloc(h->l_text + 1);
|
|
|
132 h->text[h->l_text] = 0; // make sure it is NULL terminated
|
|
|
133 bgzf_read(fp, h->text, h->l_text);
|
|
|
134 bgzf_read(fp, &h->n_targets, 4);
|
|
|
135 if (fp->is_be) ed_swap_4p(&h->n_targets);
|
|
|
136 // read reference sequence names and lengths
|
|
|
137 h->target_name = (char**)calloc(h->n_targets, sizeof(char*));
|
|
|
138 h->target_len = (uint32_t*)calloc(h->n_targets, sizeof(uint32_t));
|
|
|
139 for (i = 0; i != h->n_targets; ++i) {
|
|
|
140 bgzf_read(fp, &name_len, 4);
|
|
|
141 if (fp->is_be) ed_swap_4p(&name_len);
|
|
|
142 h->target_name[i] = (char*)calloc(name_len, 1);
|
|
|
143 bgzf_read(fp, h->target_name[i], name_len);
|
|
|
144 bgzf_read(fp, &h->target_len[i], 4);
|
|
|
145 if (fp->is_be) ed_swap_4p(&h->target_len[i]);
|
|
|
146 }
|
|
|
147 return h;
|
|
|
148 }
|
|
|
149
|
|
|
150 int bam_hdr_write(BGZF *fp, const bam_hdr_t *h)
|
|
|
151 {
|
|
|
152 char buf[4];
|
|
|
153 int32_t i, name_len, x;
|
|
|
154 // write "BAM1"
|
|
|
155 strncpy(buf, "BAM\1", 4);
|
|
|
156 bgzf_write(fp, buf, 4);
|
|
|
157 // write plain text and the number of reference sequences
|
|
|
158 if (fp->is_be) {
|
|
|
159 x = ed_swap_4(h->l_text);
|
|
|
160 bgzf_write(fp, &x, 4);
|
|
|
161 if (h->l_text) bgzf_write(fp, h->text, h->l_text);
|
|
|
162 x = ed_swap_4(h->n_targets);
|
|
|
163 bgzf_write(fp, &x, 4);
|
|
|
164 } else {
|
|
|
165 bgzf_write(fp, &h->l_text, 4);
|
|
|
166 if (h->l_text) bgzf_write(fp, h->text, h->l_text);
|
|
|
167 bgzf_write(fp, &h->n_targets, 4);
|
|
|
168 }
|
|
|
169 // write sequence names and lengths
|
|
|
170 for (i = 0; i != h->n_targets; ++i) {
|
|
|
171 char *p = h->target_name[i];
|
|
|
172 name_len = strlen(p) + 1;
|
|
|
173 if (fp->is_be) {
|
|
|
174 x = ed_swap_4(name_len);
|
|
|
175 bgzf_write(fp, &x, 4);
|
|
|
176 } else bgzf_write(fp, &name_len, 4);
|
|
|
177 bgzf_write(fp, p, name_len);
|
|
|
178 if (fp->is_be) {
|
|
|
179 x = ed_swap_4(h->target_len[i]);
|
|
|
180 bgzf_write(fp, &x, 4);
|
|
|
181 } else bgzf_write(fp, &h->target_len[i], 4);
|
|
|
182 }
|
|
|
183 bgzf_flush(fp);
|
|
|
184 return 0;
|
|
|
185 }
|
|
|
186
|
|
|
187 int bam_name2id(bam_hdr_t *h, const char *ref)
|
|
|
188 {
|
|
|
189 sdict_t *d = (sdict_t*)h->sdict;
|
|
|
190 khint_t k;
|
|
|
191 if (h->sdict == 0) {
|
|
|
192 int i, absent;
|
|
|
193 d = kh_init(s2i);
|
|
|
194 for (i = 0; i < h->n_targets; ++i) {
|
|
|
195 k = kh_put(s2i, d, h->target_name[i], &absent);
|
|
|
196 kh_val(d, k) = i;
|
|
|
197 }
|
|
|
198 h->sdict = d;
|
|
|
199 }
|
|
|
200 k = kh_get(s2i, d, ref);
|
|
|
201 return k == kh_end(d)? -1 : kh_val(d, k);
|
|
|
202 }
|
|
|
203
|
|
|
204 /*************************
|
|
|
205 *** BAM alignment I/O ***
|
|
|
206 *************************/
|
|
|
207
|
|
|
208 bam1_t *bam_init1()
|
|
|
209 {
|
|
|
210 return (bam1_t*)calloc(1, sizeof(bam1_t));
|
|
|
211 }
|
|
|
212
|
|
|
213 void bam_destroy1(bam1_t *b)
|
|
|
214 {
|
|
|
215 if (b == 0) return;
|
|
|
216 free(b->data); free(b);
|
|
|
217 }
|
|
|
218
|
|
|
219 bam1_t *bam_copy1(bam1_t *bdst, const bam1_t *bsrc)
|
|
|
220 {
|
|
|
221 uint8_t *data = bdst->data;
|
|
|
222 int m_data = bdst->m_data; // backup data and m_data
|
|
|
223 if (m_data < bsrc->l_data) { // double the capacity
|
|
|
224 m_data = bsrc->l_data; kroundup32(m_data);
|
|
|
225 data = (uint8_t*)realloc(data, m_data);
|
|
|
226 }
|
|
|
227 memcpy(data, bsrc->data, bsrc->l_data); // copy var-len data
|
|
|
228 *bdst = *bsrc; // copy the rest
|
|
|
229 // restore the backup
|
|
|
230 bdst->m_data = m_data;
|
|
|
231 bdst->data = data;
|
|
|
232 return bdst;
|
|
|
233 }
|
|
|
234
|
|
|
235 bam1_t *bam_dup1(const bam1_t *bsrc)
|
|
|
236 {
|
|
|
237 if (bsrc == NULL) return NULL;
|
|
|
238 bam1_t *bdst = bam_init1();
|
|
|
239 if (bdst == NULL) return NULL;
|
|
|
240 return bam_copy1(bdst, bsrc);
|
|
|
241 }
|
|
|
242
|
|
|
243 int bam_cigar2qlen(int n_cigar, const uint32_t *cigar)
|
|
|
244 {
|
|
|
245 int k, l;
|
|
|
246 for (k = l = 0; k < n_cigar; ++k)
|
|
|
247 if (bam_cigar_type(bam_cigar_op(cigar[k]))&1)
|
|
|
248 l += bam_cigar_oplen(cigar[k]);
|
|
|
249 return l;
|
|
|
250 }
|
|
|
251
|
|
|
252 int bam_cigar2rlen(int n_cigar, const uint32_t *cigar)
|
|
|
253 {
|
|
|
254 int k, l;
|
|
|
255 for (k = l = 0; k < n_cigar; ++k)
|
|
|
256 if (bam_cigar_type(bam_cigar_op(cigar[k]))&2)
|
|
|
257 l += bam_cigar_oplen(cigar[k]);
|
|
|
258 return l;
|
|
|
259 }
|
|
|
260
|
|
|
261 int32_t bam_endpos(const bam1_t *b)
|
|
|
262 {
|
|
|
263 if (!(b->core.flag & BAM_FUNMAP) && b->core.n_cigar > 0)
|
|
|
264 return b->core.pos + bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b));
|
|
|
265 else
|
|
|
266 return b->core.pos + 1;
|
|
|
267 }
|
|
|
268
|
|
|
269 static inline int aux_type2size(uint8_t type)
|
|
|
270 {
|
|
|
271 switch (type) {
|
|
|
272 case 'A': case 'c': case 'C':
|
|
|
273 return 1;
|
|
|
274 case 's': case 'S':
|
|
|
275 return 2;
|
|
|
276 case 'i': case 'I': case 'f':
|
|
|
277 return 4;
|
|
|
278 case 'd':
|
|
|
279 return 8;
|
|
|
280 case 'Z': case 'H': case 'B':
|
|
|
281 return type;
|
|
|
282 default:
|
|
|
283 return 0;
|
|
|
284 }
|
|
|
285 }
|
|
|
286
|
|
|
287 static void swap_data(const bam1_core_t *c, int l_data, uint8_t *data, int is_host)
|
|
|
288 {
|
|
|
289 uint8_t *s;
|
|
|
290 uint32_t *cigar = (uint32_t*)(data + c->l_qname);
|
|
|
291 uint32_t i, n;
|
|
|
292 s = data + c->n_cigar*4 + c->l_qname + c->l_qseq + (c->l_qseq + 1)/2;
|
|
|
293 for (i = 0; i < c->n_cigar; ++i) ed_swap_4p(&cigar[i]);
|
|
|
294 while (s < data + l_data) {
|
|
|
295 int size;
|
|
|
296 s += 2; // skip key
|
|
|
297 size = aux_type2size(*s); ++s; // skip type
|
|
|
298 switch (size) {
|
|
|
299 case 1: ++s; break;
|
|
|
300 case 2: ed_swap_2p(s); s += 2; break;
|
|
|
301 case 4: ed_swap_4p(s); s += 4; break;
|
|
|
302 case 8: ed_swap_8p(s); s += 8; break;
|
|
|
303 case 'Z':
|
|
|
304 case 'H':
|
|
|
305 while (*s) ++s;
|
|
|
306 ++s;
|
|
|
307 break;
|
|
|
308 case 'B':
|
|
|
309 size = aux_type2size(*s); ++s;
|
|
|
310 if (is_host) memcpy(&n, s, 4), ed_swap_4p(s);
|
|
|
311 else ed_swap_4p(s), memcpy(&n, s, 4);
|
|
|
312 s += 4;
|
|
|
313 switch (size) {
|
|
|
314 case 1: s += n; break;
|
|
|
315 case 2: for (i = 0; i < n; ++i, s += 2) ed_swap_2p(s); break;
|
|
|
316 case 4: for (i = 0; i < n; ++i, s += 4) ed_swap_4p(s); break;
|
|
|
317 case 8: for (i = 0; i < n; ++i, s += 8) ed_swap_8p(s); break;
|
|
|
318 }
|
|
|
319 break;
|
|
|
320 }
|
|
|
321 }
|
|
|
322 }
|
|
|
323
|
|
|
324 int bam_read1(BGZF *fp, bam1_t *b)
|
|
|
325 {
|
|
|
326 bam1_core_t *c = &b->core;
|
|
|
327 int32_t block_len, ret, i;
|
|
|
328 uint32_t x[8];
|
|
|
329 if ((ret = bgzf_read(fp, &block_len, 4)) != 4) {
|
|
|
330 if (ret == 0) return -1; // normal end-of-file
|
|
|
331 else return -2; // truncated
|
|
|
332 }
|
|
|
333 if (bgzf_read(fp, x, 32) != 32) return -3;
|
|
|
334 if (fp->is_be) {
|
|
|
335 ed_swap_4p(&block_len);
|
|
|
336 for (i = 0; i < 8; ++i) ed_swap_4p(x + i);
|
|
|
337 }
|
|
|
338 c->tid = x[0]; c->pos = x[1];
|
|
|
339 c->bin = x[2]>>16; c->qual = x[2]>>8&0xff; c->l_qname = x[2]&0xff;
|
|
|
340 c->flag = x[3]>>16; c->n_cigar = x[3]&0xffff;
|
|
|
341 c->l_qseq = x[4];
|
|
|
342 c->mtid = x[5]; c->mpos = x[6]; c->isize = x[7];
|
|
|
343 b->l_data = block_len - 32;
|
|
|
344 if (b->l_data < 0 || c->l_qseq < 0) return -4;
|
|
|
345 if ((char *)bam_get_aux(b) - (char *)b->data > b->l_data)
|
|
|
346 return -4;
|
|
|
347 if (b->m_data < b->l_data) {
|
|
|
348 b->m_data = b->l_data;
|
|
|
349 kroundup32(b->m_data);
|
|
|
350 b->data = (uint8_t*)realloc(b->data, b->m_data);
|
|
|
351 if (!b->data)
|
|
|
352 return -4;
|
|
|
353 }
|
|
|
354 if (bgzf_read(fp, b->data, b->l_data) != b->l_data) return -4;
|
|
|
355 //b->l_aux = b->l_data - c->n_cigar * 4 - c->l_qname - c->l_qseq - (c->l_qseq+1)/2;
|
|
|
356 if (fp->is_be) swap_data(c, b->l_data, b->data, 0);
|
|
|
357 return 4 + block_len;
|
|
|
358 }
|
|
|
359
|
|
|
360 int bam_write1(BGZF *fp, const bam1_t *b)
|
|
|
361 {
|
|
|
362 const bam1_core_t *c = &b->core;
|
|
|
363 uint32_t x[8], block_len = b->l_data + 32, y;
|
|
|
364 int i, ok;
|
|
|
365 x[0] = c->tid;
|
|
|
366 x[1] = c->pos;
|
|
|
367 x[2] = (uint32_t)c->bin<<16 | c->qual<<8 | c->l_qname;
|
|
|
368 x[3] = (uint32_t)c->flag<<16 | c->n_cigar;
|
|
|
369 x[4] = c->l_qseq;
|
|
|
370 x[5] = c->mtid;
|
|
|
371 x[6] = c->mpos;
|
|
|
372 x[7] = c->isize;
|
|
|
373 ok = (bgzf_flush_try(fp, 4 + block_len) >= 0);
|
|
|
374 if (fp->is_be) {
|
|
|
375 for (i = 0; i < 8; ++i) ed_swap_4p(x + i);
|
|
|
376 y = block_len;
|
|
|
377 if (ok) ok = (bgzf_write(fp, ed_swap_4p(&y), 4) >= 0);
|
|
|
378 swap_data(c, b->l_data, b->data, 1);
|
|
|
379 } else {
|
|
|
380 if (ok) ok = (bgzf_write(fp, &block_len, 4) >= 0);
|
|
|
381 }
|
|
|
382 if (ok) ok = (bgzf_write(fp, x, 32) >= 0);
|
|
|
383 if (ok) ok = (bgzf_write(fp, b->data, b->l_data) >= 0);
|
|
|
384 if (fp->is_be) swap_data(c, b->l_data, b->data, 0);
|
|
|
385 return ok? 4 + block_len : -1;
|
|
|
386 }
|
|
|
387
|
|
|
388 /********************
|
|
|
389 *** BAM indexing ***
|
|
|
390 ********************/
|
|
|
391
|
|
|
392 static hts_idx_t *bam_index(BGZF *fp, int min_shift)
|
|
|
393 {
|
|
|
394 int n_lvls, i, fmt;
|
|
|
395 bam1_t *b;
|
|
|
396 hts_idx_t *idx;
|
|
|
397 bam_hdr_t *h;
|
|
|
398 h = bam_hdr_read(fp);
|
|
|
399 if (min_shift > 0) {
|
|
|
400 int64_t max_len = 0, s;
|
|
|
401 for (i = 0; i < h->n_targets; ++i)
|
|
|
402 if (max_len < h->target_len[i]) max_len = h->target_len[i];
|
|
|
403 max_len += 256;
|
|
|
404 for (n_lvls = 0, s = 1<<min_shift; max_len > s; ++n_lvls, s <<= 3);
|
|
|
405 fmt = HTS_FMT_CSI;
|
|
|
406 } else min_shift = 14, n_lvls = 5, fmt = HTS_FMT_BAI;
|
|
|
407 idx = hts_idx_init(h->n_targets, fmt, bgzf_tell(fp), min_shift, n_lvls);
|
|
|
408 bam_hdr_destroy(h);
|
|
|
409 b = bam_init1();
|
|
|
410 while (bam_read1(fp, b) >= 0) {
|
|
|
411 int l, ret;
|
|
|
412 l = bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b));
|
|
|
413 if (l == 0) l = 1; // no zero-length records
|
|
|
414 ret = hts_idx_push(idx, b->core.tid, b->core.pos, b->core.pos + l, bgzf_tell(fp), !(b->core.flag&BAM_FUNMAP));
|
|
|
415 if (ret < 0)
|
|
|
416 {
|
|
|
417 // unsorted
|
|
|
418 bam_destroy1(b);
|
|
|
419 hts_idx_destroy(idx);
|
|
|
420 return NULL;
|
|
|
421 }
|
|
|
422 }
|
|
|
423 hts_idx_finish(idx, bgzf_tell(fp));
|
|
|
424 bam_destroy1(b);
|
|
|
425 return idx;
|
|
|
426 }
|
|
|
427
|
|
|
428 int bam_index_build(const char *fn, int min_shift)
|
|
|
429 {
|
|
|
430 hts_idx_t *idx;
|
|
|
431 htsFile *fp;
|
|
|
432 int ret = 0;
|
|
|
433
|
|
|
434 if ((fp = hts_open(fn, "r")) == 0) return -1;
|
|
|
435 switch (fp->format.format) {
|
|
|
436 case cram:
|
|
|
437 ret = cram_index_build(fp->fp.cram, fn);
|
|
|
438 break;
|
|
|
439
|
|
|
440 case bam:
|
|
|
441 idx = bam_index(fp->fp.bgzf, min_shift);
|
|
|
442 if (idx) {
|
|
|
443 hts_idx_save(idx, fn, (min_shift > 0)? HTS_FMT_CSI : HTS_FMT_BAI);
|
|
|
444 hts_idx_destroy(idx);
|
|
|
445 }
|
|
|
446 else ret = -1;
|
|
|
447 break;
|
|
|
448
|
|
|
449 default:
|
|
|
450 ret = -1;
|
|
|
451 break;
|
|
|
452 }
|
|
|
453 hts_close(fp);
|
|
|
454
|
|
|
455 return ret;
|
|
|
456 }
|
|
|
457
|
|
|
458 static int bam_readrec(BGZF *fp, void *ignored, void *bv, int *tid, int *beg, int *end)
|
|
|
459 {
|
|
|
460 bam1_t *b = bv;
|
|
|
461 int ret;
|
|
|
462 if ((ret = bam_read1(fp, b)) >= 0) {
|
|
|
463 *tid = b->core.tid; *beg = b->core.pos;
|
|
|
464 *end = b->core.pos + (b->core.n_cigar? bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b)) : 1);
|
|
|
465 }
|
|
|
466 return ret;
|
|
|
467 }
|
|
|
468
|
|
|
469 // This is used only with read_rest=1 iterators, so need not set tid/beg/end.
|
|
|
470 static int cram_readrec(BGZF *ignored, void *fpv, void *bv, int *tid, int *beg, int *end)
|
|
|
471 {
|
|
|
472 htsFile *fp = fpv;
|
|
|
473 bam1_t *b = bv;
|
|
|
474 return cram_get_bam_seq(fp->fp.cram, &b);
|
|
|
475 }
|
|
|
476
|
|
|
477 // This is used only with read_rest=1 iterators, so need not set tid/beg/end.
|
|
|
478 static int sam_bam_cram_readrec(BGZF *bgzfp, void *fpv, void *bv, int *tid, int *beg, int *end)
|
|
|
479 {
|
|
|
480 htsFile *fp = fpv;
|
|
|
481 bam1_t *b = bv;
|
|
|
482 switch (fp->format.format) {
|
|
|
483 case bam: return bam_read1(bgzfp, b);
|
|
|
484 case cram: return cram_get_bam_seq(fp->fp.cram, &b);
|
|
|
485 default:
|
|
|
486 // TODO Need headers available to implement this for SAM files
|
|
|
487 fprintf(stderr, "[sam_bam_cram_readrec] Not implemented for SAM files -- Exiting\n");
|
|
|
488 abort();
|
|
|
489 }
|
|
|
490 }
|
|
|
491
|
|
|
492 // The CRAM implementation stores the loaded index within the cram_fd rather
|
|
|
493 // than separately as is done elsewhere in htslib. So if p is a pointer to
|
|
|
494 // an hts_idx_t with p->fmt == HTS_FMT_CRAI, then it actually points to an
|
|
|
495 // hts_cram_idx_t and should be cast accordingly.
|
|
|
496 typedef struct hts_cram_idx_t {
|
|
|
497 int fmt;
|
|
|
498 cram_fd *cram;
|
|
|
499 } hts_cram_idx_t;
|
|
|
500
|
|
|
501 hts_idx_t *sam_index_load(samFile *fp, const char *fn)
|
|
|
502 {
|
|
|
503 switch (fp->format.format) {
|
|
|
504 case bam:
|
|
|
505 return bam_index_load(fn);
|
|
|
506
|
|
|
507 case cram: {
|
|
|
508 if (cram_index_load(fp->fp.cram, fn) < 0) return NULL;
|
|
|
509 // Cons up a fake "index" just pointing at the associated cram_fd:
|
|
|
510 hts_cram_idx_t *idx = malloc(sizeof (hts_cram_idx_t));
|
|
|
511 if (idx == NULL) return NULL;
|
|
|
512 idx->fmt = HTS_FMT_CRAI;
|
|
|
513 idx->cram = fp->fp.cram;
|
|
|
514 return (hts_idx_t *) idx;
|
|
|
515 }
|
|
|
516
|
|
|
517 default:
|
|
|
518 return NULL; // TODO Would use tbx_index_load if it returned hts_idx_t
|
|
|
519 }
|
|
|
520 }
|
|
|
521
|
|
|
522 static hts_itr_t *cram_itr_query(const hts_idx_t *idx, int tid, int beg, int end, hts_readrec_func *readrec)
|
|
|
523 {
|
|
|
524 const hts_cram_idx_t *cidx = (const hts_cram_idx_t *) idx;
|
|
|
525 hts_itr_t *iter = (hts_itr_t *) calloc(1, sizeof(hts_itr_t));
|
|
|
526 if (iter == NULL) return NULL;
|
|
|
527
|
|
|
528 // Cons up a dummy iterator for which hts_itr_next() will simply invoke
|
|
|
529 // the readrec function:
|
|
|
530 iter->read_rest = 1;
|
|
|
531 iter->off = NULL;
|
|
|
532 iter->bins.a = NULL;
|
|
|
533 iter->readrec = readrec;
|
|
|
534
|
|
|
535 if (tid >= 0) {
|
|
|
536 cram_range r = { tid, beg+1, end };
|
|
|
537 if (cram_set_option(cidx->cram, CRAM_OPT_RANGE, &r) != 0) { free(iter); return NULL; }
|
|
|
538 iter->curr_off = 0;
|
|
|
539 // The following fields are not required by hts_itr_next(), but are
|
|
|
540 // filled in in case user code wants to look at them.
|
|
|
541 iter->tid = tid;
|
|
|
542 iter->beg = beg;
|
|
|
543 iter->end = end;
|
|
|
544 }
|
|
|
545 else switch (tid) {
|
|
|
546 case HTS_IDX_REST:
|
|
|
547 iter->curr_off = 0;
|
|
|
548 break;
|
|
|
549 case HTS_IDX_NONE:
|
|
|
550 iter->curr_off = 0;
|
|
|
551 iter->finished = 1;
|
|
|
552 break;
|
|
|
553 default:
|
|
|
554 fprintf(stderr, "[cram_itr_query] tid=%d not implemented for CRAM files -- Exiting\n", tid);
|
|
|
555 abort();
|
|
|
556 break;
|
|
|
557 }
|
|
|
558
|
|
|
559 return iter;
|
|
|
560 }
|
|
|
561
|
|
|
562 hts_itr_t *sam_itr_queryi(const hts_idx_t *idx, int tid, int beg, int end)
|
|
|
563 {
|
|
|
564 const hts_cram_idx_t *cidx = (const hts_cram_idx_t *) idx;
|
|
|
565 if (idx == NULL)
|
|
|
566 return hts_itr_query(NULL, tid, beg, end, sam_bam_cram_readrec);
|
|
|
567 else if (cidx->fmt == HTS_FMT_CRAI)
|
|
|
568 return cram_itr_query(idx, tid, beg, end, cram_readrec);
|
|
|
569 else
|
|
|
570 return hts_itr_query(idx, tid, beg, end, bam_readrec);
|
|
|
571 }
|
|
|
572
|
|
|
573 static int cram_name2id(void *fdv, const char *ref)
|
|
|
574 {
|
|
|
575 cram_fd *fd = (cram_fd *) fdv;
|
|
|
576 return sam_hdr_name2ref(fd->header, ref);
|
|
|
577 }
|
|
|
578
|
|
|
579 hts_itr_t *sam_itr_querys(const hts_idx_t *idx, bam_hdr_t *hdr, const char *region)
|
|
|
580 {
|
|
|
581 const hts_cram_idx_t *cidx = (const hts_cram_idx_t *) idx;
|
|
|
582 if (cidx->fmt == HTS_FMT_CRAI)
|
|
|
583 return hts_itr_querys(idx, region, cram_name2id, cidx->cram, cram_itr_query, cram_readrec);
|
|
|
584 else
|
|
|
585 return hts_itr_querys(idx, region, (hts_name2id_f)(bam_name2id), hdr, hts_itr_query, bam_readrec);
|
|
|
586 }
|
|
|
587
|
|
|
588 /**********************
|
|
|
589 *** SAM header I/O ***
|
|
|
590 **********************/
|
|
|
591
|
|
|
592 #include "htslib/kseq.h"
|
|
|
593 #include "htslib/kstring.h"
|
|
|
594
|
|
|
595 bam_hdr_t *sam_hdr_parse(int l_text, const char *text)
|
|
|
596 {
|
|
|
597 const char *q, *r, *p;
|
|
|
598 khash_t(s2i) *d;
|
|
|
599 d = kh_init(s2i);
|
|
|
600 for (p = text; *p; ++p) {
|
|
|
601 if (strncmp(p, "@SQ", 3) == 0) {
|
|
|
602 char *sn = 0;
|
|
|
603 int ln = -1;
|
|
|
604 for (q = p + 4;; ++q) {
|
|
|
605 if (strncmp(q, "SN:", 3) == 0) {
|
|
|
606 q += 3;
|
|
|
607 for (r = q; *r != '\t' && *r != '\n'; ++r);
|
|
|
608 sn = (char*)calloc(r - q + 1, 1);
|
|
|
609 strncpy(sn, q, r - q);
|
|
|
610 q = r;
|
|
|
611 } else if (strncmp(q, "LN:", 3) == 0)
|
|
|
612 ln = strtol(q + 3, (char**)&q, 10);
|
|
|
613 while (*q != '\t' && *q != '\n') ++q;
|
|
|
614 if (*q == '\n') break;
|
|
|
615 }
|
|
|
616 p = q;
|
|
|
617 if (sn && ln >= 0) {
|
|
|
618 khint_t k;
|
|
|
619 int absent;
|
|
|
620 k = kh_put(s2i, d, sn, &absent);
|
|
|
621 if (!absent) {
|
|
|
622 if (hts_verbose >= 2)
|
|
|
623 fprintf(stderr, "[W::%s] duplicated sequence '%s'\n", __func__, sn);
|
|
|
624 free(sn);
|
|
|
625 } else kh_val(d, k) = (int64_t)(kh_size(d) - 1)<<32 | ln;
|
|
|
626 }
|
|
|
627 }
|
|
|
628 while (*p != '\n') ++p;
|
|
|
629 }
|
|
|
630 return hdr_from_dict(d);
|
|
|
631 }
|
|
|
632
|
|
|
633 bam_hdr_t *sam_hdr_read(htsFile *fp)
|
|
|
634 {
|
|
|
635 switch (fp->format.format) {
|
|
|
636 case bam:
|
|
|
637 return bam_hdr_read(fp->fp.bgzf);
|
|
|
638
|
|
|
639 case cram:
|
|
|
640 return cram_header_to_bam(fp->fp.cram->header);
|
|
|
641
|
|
|
642 case sam: {
|
|
|
643 kstring_t str;
|
|
|
644 bam_hdr_t *h;
|
|
|
645 int has_SQ = 0;
|
|
|
646 str.l = str.m = 0; str.s = 0;
|
|
|
647 while (hts_getline(fp, KS_SEP_LINE, &fp->line) >= 0) {
|
|
|
648 if (fp->line.s[0] != '@') break;
|
|
|
649 if (fp->line.l > 3 && strncmp(fp->line.s,"@SQ",3) == 0) has_SQ = 1;
|
|
|
650 kputsn(fp->line.s, fp->line.l, &str);
|
|
|
651 kputc('\n', &str);
|
|
|
652 }
|
|
|
653 if (! has_SQ && fp->fn_aux) {
|
|
|
654 char line[2048];
|
|
|
655 FILE *f = fopen(fp->fn_aux, "r");
|
|
|
656 if (f == NULL) return NULL;
|
|
|
657 while (fgets(line, sizeof line, f)) {
|
|
|
658 const char *name = strtok(line, "\t");
|
|
|
659 const char *length = strtok(NULL, "\t");
|
|
|
660 ksprintf(&str, "@SQ\tSN:%s\tLN:%s\n", name, length);
|
|
|
661 }
|
|
|
662 fclose(f);
|
|
|
663 }
|
|
|
664 if (str.l == 0) kputsn("", 0, &str);
|
|
|
665 h = sam_hdr_parse(str.l, str.s);
|
|
|
666 h->l_text = str.l; h->text = str.s;
|
|
|
667 return h;
|
|
|
668 }
|
|
|
669
|
|
|
670 default:
|
|
|
671 abort();
|
|
|
672 }
|
|
|
673 }
|
|
|
674
|
|
|
675 int sam_hdr_write(htsFile *fp, const bam_hdr_t *h)
|
|
|
676 {
|
|
|
677 switch (fp->format.format) {
|
|
|
678 case binary_format:
|
|
|
679 fp->format.category = sequence_data;
|
|
|
680 fp->format.format = bam;
|
|
|
681 /* fall-through */
|
|
|
682 case bam:
|
|
|
683 bam_hdr_write(fp->fp.bgzf, h);
|
|
|
684 break;
|
|
|
685
|
|
|
686 case cram: {
|
|
|
687 cram_fd *fd = fp->fp.cram;
|
|
|
688 if (cram_set_header(fd, bam_header_to_cram((bam_hdr_t *)h)) < 0) return -1;
|
|
|
689 if (fp->fn_aux)
|
|
|
690 cram_load_reference(fd, fp->fn_aux);
|
|
|
691 if (cram_write_SAM_hdr(fd, fd->header) < 0) return -1;
|
|
|
692 }
|
|
|
693 break;
|
|
|
694
|
|
|
695 case text_format:
|
|
|
696 fp->format.category = sequence_data;
|
|
|
697 fp->format.format = sam;
|
|
|
698 /* fall-through */
|
|
|
699 case sam: {
|
|
|
700 char *p;
|
|
|
701 hputs(h->text, fp->fp.hfile);
|
|
|
702 p = strstr(h->text, "@SQ\t"); // FIXME: we need a loop to make sure "@SQ\t" does not match something unwanted!!!
|
|
|
703 if (p == 0) {
|
|
|
704 int i;
|
|
|
705 for (i = 0; i < h->n_targets; ++i) {
|
|
|
706 fp->line.l = 0;
|
|
|
707 kputsn("@SQ\tSN:", 7, &fp->line); kputs(h->target_name[i], &fp->line);
|
|
|
708 kputsn("\tLN:", 4, &fp->line); kputw(h->target_len[i], &fp->line); kputc('\n', &fp->line);
|
|
|
709 if ( hwrite(fp->fp.hfile, fp->line.s, fp->line.l) != fp->line.l ) return -1;
|
|
|
710 }
|
|
|
711 }
|
|
|
712 if ( hflush(fp->fp.hfile) != 0 ) return -1;
|
|
|
713 }
|
|
|
714 break;
|
|
|
715
|
|
|
716 default:
|
|
|
717 abort();
|
|
|
718 }
|
|
|
719 return 0;
|
|
|
720 }
|
|
|
721
|
|
|
722 /**********************
|
|
|
723 *** SAM record I/O ***
|
|
|
724 **********************/
|
|
|
725
|
|
|
726 int sam_parse1(kstring_t *s, bam_hdr_t *h, bam1_t *b)
|
|
|
727 {
|
|
|
728 #define _read_token(_p) (_p); for (; *(_p) && *(_p) != '\t'; ++(_p)); if (*(_p) != '\t') goto err_ret; *(_p)++ = 0
|
|
|
729 #define _read_token_aux(_p) (_p); for (; *(_p) && *(_p) != '\t'; ++(_p)); *(_p)++ = 0 // this is different in that it does not test *(_p)=='\t'
|
|
|
730 #define _get_mem(type_t, _x, _s, _l) ks_resize((_s), (_s)->l + (_l)); *(_x) = (type_t*)((_s)->s + (_s)->l); (_s)->l += (_l)
|
|
|
731 #define _parse_err(cond, msg) do { if ((cond) && hts_verbose >= 1) { fprintf(stderr, "[E::%s] " msg "\n", __func__); goto err_ret; } } while (0)
|
|
|
732 #define _parse_warn(cond, msg) if ((cond) && hts_verbose >= 2) fprintf(stderr, "[W::%s] " msg "\n", __func__)
|
|
|
733
|
|
|
734 uint8_t *t;
|
|
|
735 char *p = s->s, *q;
|
|
|
736 int i;
|
|
|
737 kstring_t str;
|
|
|
738 bam1_core_t *c = &b->core;
|
|
|
739
|
|
|
740 str.l = b->l_data = 0;
|
|
|
741 str.s = (char*)b->data; str.m = b->m_data;
|
|
|
742 memset(c, 0, 32);
|
|
|
743 if (h->cigar_tab == 0) {
|
|
|
744 h->cigar_tab = (int8_t*) malloc(128);
|
|
|
745 for (i = 0; i < 128; ++i)
|
|
|
746 h->cigar_tab[i] = -1;
|
|
|
747 for (i = 0; BAM_CIGAR_STR[i]; ++i)
|
|
|
748 h->cigar_tab[(int)BAM_CIGAR_STR[i]] = i;
|
|
|
749 }
|
|
|
750 // qname
|
|
|
751 q = _read_token(p);
|
|
|
752 kputsn_(q, p - q, &str);
|
|
|
753 c->l_qname = p - q;
|
|
|
754 // flag
|
|
|
755 c->flag = strtol(p, &p, 0);
|
|
|
756 if (*p++ != '\t') goto err_ret; // malformated flag
|
|
|
757 // chr
|
|
|
758 q = _read_token(p);
|
|
|
759 if (strcmp(q, "*")) {
|
|
|
760 _parse_err(h->n_targets == 0, "missing SAM header");
|
|
|
761 c->tid = bam_name2id(h, q);
|
|
|
762 _parse_warn(c->tid < 0, "urecognized reference name; treated as unmapped");
|
|
|
763 } else c->tid = -1;
|
|
|
764 // pos
|
|
|
765 c->pos = strtol(p, &p, 10) - 1;
|
|
|
766 if (*p++ != '\t') goto err_ret;
|
|
|
767 if (c->pos < 0 && c->tid >= 0) {
|
|
|
768 _parse_warn(1, "mapped query cannot have zero coordinate; treated as unmapped");
|
|
|
769 c->tid = -1;
|
|
|
770 }
|
|
|
771 if (c->tid < 0) c->flag |= BAM_FUNMAP;
|
|
|
772 // mapq
|
|
|
773 c->qual = strtol(p, &p, 10);
|
|
|
774 if (*p++ != '\t') goto err_ret;
|
|
|
775 // cigar
|
|
|
776 if (*p != '*') {
|
|
|
777 uint32_t *cigar;
|
|
|
778 size_t n_cigar = 0;
|
|
|
779 for (q = p; *p && *p != '\t'; ++p)
|
|
|
780 if (!isdigit(*p)) ++n_cigar;
|
|
|
781 if (*p++ != '\t') goto err_ret;
|
|
|
782 _parse_err(n_cigar >= 65536, "too many CIGAR operations");
|
|
|
783 c->n_cigar = n_cigar;
|
|
|
784 _get_mem(uint32_t, &cigar, &str, c->n_cigar<<2);
|
|
|
785 for (i = 0; i < c->n_cigar; ++i, ++q) {
|
|
|
786 int op;
|
|
|
787 cigar[i] = strtol(q, &q, 10)<<BAM_CIGAR_SHIFT;
|
|
|
788 op = (uint8_t)*q >= 128? -1 : h->cigar_tab[(int)*q];
|
|
|
789 _parse_err(op < 0, "unrecognized CIGAR operator");
|
|
|
790 cigar[i] |= op;
|
|
|
791 }
|
|
|
792 i = bam_cigar2rlen(c->n_cigar, cigar);
|
|
|
793 } else {
|
|
|
794 _parse_warn(!(c->flag&BAM_FUNMAP), "mapped query must have a CIGAR; treated as unmapped");
|
|
|
795 c->flag |= BAM_FUNMAP;
|
|
|
796 q = _read_token(p);
|
|
|
797 i = 1;
|
|
|
798 }
|
|
|
799 c->bin = hts_reg2bin(c->pos, c->pos + i, 14, 5);
|
|
|
800 // mate chr
|
|
|
801 q = _read_token(p);
|
|
|
802 if (strcmp(q, "=") == 0) c->mtid = c->tid;
|
|
|
803 else if (strcmp(q, "*") == 0) c->mtid = -1;
|
|
|
804 else c->mtid = bam_name2id(h, q);
|
|
|
805 // mpos
|
|
|
806 c->mpos = strtol(p, &p, 10) - 1;
|
|
|
807 if (*p++ != '\t') goto err_ret;
|
|
|
808 if (c->mpos < 0 && c->mtid >= 0) {
|
|
|
809 _parse_warn(1, "mapped mate cannot have zero coordinate; treated as unmapped");
|
|
|
810 c->mtid = -1;
|
|
|
811 }
|
|
|
812 // tlen
|
|
|
813 c->isize = strtol(p, &p, 10);
|
|
|
814 if (*p++ != '\t') goto err_ret;
|
|
|
815 // seq
|
|
|
816 q = _read_token(p);
|
|
|
817 if (strcmp(q, "*")) {
|
|
|
818 c->l_qseq = p - q - 1;
|
|
|
819 i = bam_cigar2qlen(c->n_cigar, (uint32_t*)(str.s + c->l_qname));
|
|
|
820 _parse_err(c->n_cigar && i != c->l_qseq, "CIGAR and query sequence are of different length");
|
|
|
821 i = (c->l_qseq + 1) >> 1;
|
|
|
822 _get_mem(uint8_t, &t, &str, i);
|
|
|
823 memset(t, 0, i);
|
|
|
824 for (i = 0; i < c->l_qseq; ++i)
|
|
|
825 t[i>>1] |= seq_nt16_table[(int)q[i]] << ((~i&1)<<2);
|
|
|
826 } else c->l_qseq = 0;
|
|
|
827 // qual
|
|
|
828 q = _read_token_aux(p);
|
|
|
829 _get_mem(uint8_t, &t, &str, c->l_qseq);
|
|
|
830 if (strcmp(q, "*")) {
|
|
|
831 _parse_err(p - q - 1 != c->l_qseq, "SEQ and QUAL are of different length");
|
|
|
832 for (i = 0; i < c->l_qseq; ++i) t[i] = q[i] - 33;
|
|
|
833 } else memset(t, 0xff, c->l_qseq);
|
|
|
834 // aux
|
|
|
835 // Note that (like the bam1_core_t fields) this aux data in b->data is
|
|
|
836 // stored in host endianness; so there is no byte swapping needed here.
|
|
|
837 while (p < s->s + s->l) {
|
|
|
838 uint8_t type;
|
|
|
839 q = _read_token_aux(p); // FIXME: can be accelerated for long 'B' arrays
|
|
|
840 _parse_err(p - q - 1 < 6, "incomplete aux field");
|
|
|
841 kputsn_(q, 2, &str);
|
|
|
842 q += 3; type = *q++; ++q; // q points to value
|
|
|
843 if (type == 'A' || type == 'a' || type == 'c' || type == 'C') {
|
|
|
844 kputc_('A', &str);
|
|
|
845 kputc_(*q, &str);
|
|
|
846 } else if (type == 'i' || type == 'I') {
|
|
|
847 if (*q == '-') {
|
|
|
848 long x = strtol(q, &q, 10);
|
|
|
849 if (x >= INT8_MIN) {
|
|
|
850 kputc_('c', &str); kputc_(x, &str);
|
|
|
851 } else if (x >= INT16_MIN) {
|
|
|
852 int16_t y = x;
|
|
|
853 kputc_('s', &str); kputsn_((char*)&y, 2, &str);
|
|
|
854 } else {
|
|
|
855 int32_t y = x;
|
|
|
856 kputc_('i', &str); kputsn_(&y, 4, &str);
|
|
|
857 }
|
|
|
858 } else {
|
|
|
859 unsigned long x = strtoul(q, &q, 10);
|
|
|
860 if (x <= UINT8_MAX) {
|
|
|
861 kputc_('C', &str); kputc_(x, &str);
|
|
|
862 } else if (x <= UINT16_MAX) {
|
|
|
863 uint16_t y = x;
|
|
|
864 kputc_('S', &str); kputsn_(&y, 2, &str);
|
|
|
865 } else {
|
|
|
866 uint32_t y = x;
|
|
|
867 kputc_('I', &str); kputsn_(&y, 4, &str);
|
|
|
868 }
|
|
|
869 }
|
|
|
870 } else if (type == 'f') {
|
|
|
871 float x;
|
|
|
872 x = strtod(q, &q);
|
|
|
873 kputc_('f', &str); kputsn_(&x, 4, &str);
|
|
|
874 } else if (type == 'd') {
|
|
|
875 double x;
|
|
|
876 x = strtod(q, &q);
|
|
|
877 kputc_('d', &str); kputsn_(&x, 8, &str);
|
|
|
878 } else if (type == 'Z' || type == 'H') {
|
|
|
879 kputc_(type, &str);kputsn_(q, p - q, &str); // note that this include the trailing NULL
|
|
|
880 } else if (type == 'B') {
|
|
|
881 int32_t n;
|
|
|
882 char *r;
|
|
|
883 _parse_err(p - q - 1 < 3, "incomplete B-typed aux field");
|
|
|
884 type = *q++; // q points to the first ',' following the typing byte
|
|
|
885 for (r = q, n = 0; *r; ++r)
|
|
|
886 if (*r == ',') ++n;
|
|
|
887 kputc_('B', &str); kputc_(type, &str); kputsn_(&n, 4, &str);
|
|
|
888 // FIXME: to evaluate which is faster: a) aligned array and then memmove(); b) unaligned array; c) kputsn_()
|
|
|
889 if (type == 'c') while (q + 1 < p) { int8_t x = strtol(q + 1, &q, 0); kputc_(x, &str); }
|
|
|
890 else if (type == 'C') while (q + 1 < p) { uint8_t x = strtoul(q + 1, &q, 0); kputc_(x, &str); }
|
|
|
891 else if (type == 's') while (q + 1 < p) { int16_t x = strtol(q + 1, &q, 0); kputsn_(&x, 2, &str); }
|
|
|
892 else if (type == 'S') while (q + 1 < p) { uint16_t x = strtoul(q + 1, &q, 0); kputsn_(&x, 2, &str); }
|
|
|
893 else if (type == 'i') while (q + 1 < p) { int32_t x = strtol(q + 1, &q, 0); kputsn_(&x, 4, &str); }
|
|
|
894 else if (type == 'I') while (q + 1 < p) { uint32_t x = strtoul(q + 1, &q, 0); kputsn_(&x, 4, &str); }
|
|
|
895 else if (type == 'f') while (q + 1 < p) { float x = strtod(q + 1, &q); kputsn_(&x, 4, &str); }
|
|
|
896 else _parse_err(1, "unrecognized type");
|
|
|
897 } else _parse_err(1, "unrecognized type");
|
|
|
898 }
|
|
|
899 b->data = (uint8_t*)str.s; b->l_data = str.l; b->m_data = str.m;
|
|
|
900 return 0;
|
|
|
901
|
|
|
902 #undef _parse_warn
|
|
|
903 #undef _parse_err
|
|
|
904 #undef _get_mem
|
|
|
905 #undef _read_token_aux
|
|
|
906 #undef _read_token
|
|
|
907 err_ret:
|
|
|
908 b->data = (uint8_t*)str.s; b->l_data = str.l; b->m_data = str.m;
|
|
|
909 return -2;
|
|
|
910 }
|
|
|
911
|
|
|
912 int sam_read1(htsFile *fp, bam_hdr_t *h, bam1_t *b)
|
|
|
913 {
|
|
|
914 switch (fp->format.format) {
|
|
|
915 case bam: {
|
|
|
916 int r = bam_read1(fp->fp.bgzf, b);
|
|
|
917 if (r >= 0) {
|
|
|
918 if (b->core.tid >= h->n_targets || b->core.tid < -1 ||
|
|
|
919 b->core.mtid >= h->n_targets || b->core.mtid < -1)
|
|
|
920 return -3;
|
|
|
921 }
|
|
|
922 return r;
|
|
|
923 }
|
|
|
924
|
|
|
925 case cram:
|
|
|
926 return cram_get_bam_seq(fp->fp.cram, &b);
|
|
|
927
|
|
|
928 case sam: {
|
|
|
929 int ret;
|
|
|
930 err_recover:
|
|
|
931 if (fp->line.l == 0) {
|
|
|
932 ret = hts_getline(fp, KS_SEP_LINE, &fp->line);
|
|
|
933 if (ret < 0) return -1;
|
|
|
934 }
|
|
|
935 ret = sam_parse1(&fp->line, h, b);
|
|
|
936 fp->line.l = 0;
|
|
|
937 if (ret < 0) {
|
|
|
938 if (hts_verbose >= 1)
|
|
|
939 fprintf(stderr, "[W::%s] parse error at line %lld\n", __func__, (long long)fp->lineno);
|
|
|
940 if (h->ignore_sam_err) goto err_recover;
|
|
|
941 }
|
|
|
942 return ret;
|
|
|
943 }
|
|
|
944
|
|
|
945 default:
|
|
|
946 abort();
|
|
|
947 }
|
|
|
948 }
|
|
|
949
|
|
|
950 int sam_format1(const bam_hdr_t *h, const bam1_t *b, kstring_t *str)
|
|
|
951 {
|
|
|
952 int i;
|
|
|
953 uint8_t *s;
|
|
|
954 const bam1_core_t *c = &b->core;
|
|
|
955
|
|
|
956 str->l = 0;
|
|
|
957 kputsn(bam_get_qname(b), c->l_qname-1, str); kputc('\t', str); // query name
|
|
|
958 kputw(c->flag, str); kputc('\t', str); // flag
|
|
|
959 if (c->tid >= 0) { // chr
|
|
|
960 kputs(h->target_name[c->tid] , str);
|
|
|
961 kputc('\t', str);
|
|
|
962 } else kputsn("*\t", 2, str);
|
|
|
963 kputw(c->pos + 1, str); kputc('\t', str); // pos
|
|
|
964 kputw(c->qual, str); kputc('\t', str); // qual
|
|
|
965 if (c->n_cigar) { // cigar
|
|
|
966 uint32_t *cigar = bam_get_cigar(b);
|
|
|
967 for (i = 0; i < c->n_cigar; ++i) {
|
|
|
968 kputw(bam_cigar_oplen(cigar[i]), str);
|
|
|
969 kputc(bam_cigar_opchr(cigar[i]), str);
|
|
|
970 }
|
|
|
971 } else kputc('*', str);
|
|
|
972 kputc('\t', str);
|
|
|
973 if (c->mtid < 0) kputsn("*\t", 2, str); // mate chr
|
|
|
974 else if (c->mtid == c->tid) kputsn("=\t", 2, str);
|
|
|
975 else {
|
|
|
976 kputs(h->target_name[c->mtid], str);
|
|
|
977 kputc('\t', str);
|
|
|
978 }
|
|
|
979 kputw(c->mpos + 1, str); kputc('\t', str); // mate pos
|
|
|
980 kputw(c->isize, str); kputc('\t', str); // template len
|
|
|
981 if (c->l_qseq) { // seq and qual
|
|
|
982 uint8_t *s = bam_get_seq(b);
|
|
|
983 for (i = 0; i < c->l_qseq; ++i) kputc("=ACMGRSVTWYHKDBN"[bam_seqi(s, i)], str);
|
|
|
984 kputc('\t', str);
|
|
|
985 s = bam_get_qual(b);
|
|
|
986 if (s[0] == 0xff) kputc('*', str);
|
|
|
987 else for (i = 0; i < c->l_qseq; ++i) kputc(s[i] + 33, str);
|
|
|
988 } else kputsn("*\t*", 3, str);
|
|
|
989 s = bam_get_aux(b); // aux
|
|
|
990 while (s+4 <= b->data + b->l_data) {
|
|
|
991 uint8_t type, key[2];
|
|
|
992 key[0] = s[0]; key[1] = s[1];
|
|
|
993 s += 2; type = *s++;
|
|
|
994 kputc('\t', str); kputsn((char*)key, 2, str); kputc(':', str);
|
|
|
995 if (type == 'A') {
|
|
|
996 kputsn("A:", 2, str);
|
|
|
997 kputc(*s, str);
|
|
|
998 ++s;
|
|
|
999 } else if (type == 'C') {
|
|
|
1000 kputsn("i:", 2, str);
|
|
|
1001 kputw(*s, str);
|
|
|
1002 ++s;
|
|
|
1003 } else if (type == 'c') {
|
|
|
1004 kputsn("i:", 2, str);
|
|
|
1005 kputw(*(int8_t*)s, str);
|
|
|
1006 ++s;
|
|
|
1007 } else if (type == 'S') {
|
|
|
1008 if (s+2 <= b->data + b->l_data) {
|
|
|
1009 kputsn("i:", 2, str);
|
|
|
1010 kputw(*(uint16_t*)s, str);
|
|
|
1011 s += 2;
|
|
|
1012 } else return -1;
|
|
|
1013 } else if (type == 's') {
|
|
|
1014 if (s+2 <= b->data + b->l_data) {
|
|
|
1015 kputsn("i:", 2, str);
|
|
|
1016 kputw(*(int16_t*)s, str);
|
|
|
1017 s += 2;
|
|
|
1018 } else return -1;
|
|
|
1019 } else if (type == 'I') {
|
|
|
1020 if (s+4 <= b->data + b->l_data) {
|
|
|
1021 kputsn("i:", 2, str);
|
|
|
1022 kputuw(*(uint32_t*)s, str);
|
|
|
1023 s += 4;
|
|
|
1024 } else return -1;
|
|
|
1025 } else if (type == 'i') {
|
|
|
1026 if (s+4 <= b->data + b->l_data) {
|
|
|
1027 kputsn("i:", 2, str);
|
|
|
1028 kputw(*(int32_t*)s, str);
|
|
|
1029 s += 4;
|
|
|
1030 } else return -1;
|
|
|
1031 } else if (type == 'f') {
|
|
|
1032 if (s+4 <= b->data + b->l_data) {
|
|
|
1033 ksprintf(str, "f:%g", *(float*)s);
|
|
|
1034 s += 4;
|
|
|
1035 } else return -1;
|
|
|
1036
|
|
|
1037 } else if (type == 'd') {
|
|
|
1038 if (s+8 <= b->data + b->l_data) {
|
|
|
1039 ksprintf(str, "d:%g", *(double*)s);
|
|
|
1040 s += 8;
|
|
|
1041 } else return -1;
|
|
|
1042 } else if (type == 'Z' || type == 'H') {
|
|
|
1043 kputc(type, str); kputc(':', str);
|
|
|
1044 while (s < b->data + b->l_data && *s) kputc(*s++, str);
|
|
|
1045 if (s >= b->data + b->l_data)
|
|
|
1046 return -1;
|
|
|
1047 ++s;
|
|
|
1048 } else if (type == 'B') {
|
|
|
1049 uint8_t sub_type = *(s++);
|
|
|
1050 int32_t n;
|
|
|
1051 memcpy(&n, s, 4);
|
|
|
1052 s += 4; // no point to the start of the array
|
|
|
1053 if (s + n >= b->data + b->l_data)
|
|
|
1054 return -1;
|
|
|
1055 kputsn("B:", 2, str); kputc(sub_type, str); // write the typing
|
|
|
1056 for (i = 0; i < n; ++i) { // FIXME: for better performance, put the loop after "if"
|
|
|
1057 kputc(',', str);
|
|
|
1058 if ('c' == sub_type) { kputw(*(int8_t*)s, str); ++s; }
|
|
|
1059 else if ('C' == sub_type) { kputw(*(uint8_t*)s, str); ++s; }
|
|
|
1060 else if ('s' == sub_type) { kputw(*(int16_t*)s, str); s += 2; }
|
|
|
1061 else if ('S' == sub_type) { kputw(*(uint16_t*)s, str); s += 2; }
|
|
|
1062 else if ('i' == sub_type) { kputw(*(int32_t*)s, str); s += 4; }
|
|
|
1063 else if ('I' == sub_type) { kputuw(*(uint32_t*)s, str); s += 4; }
|
|
|
1064 else if ('f' == sub_type) { ksprintf(str, "%g", *(float*)s); s += 4; }
|
|
|
1065 }
|
|
|
1066 }
|
|
|
1067 }
|
|
|
1068 return str->l;
|
|
|
1069 }
|
|
|
1070
|
|
|
1071 int sam_write1(htsFile *fp, const bam_hdr_t *h, const bam1_t *b)
|
|
|
1072 {
|
|
|
1073 switch (fp->format.format) {
|
|
|
1074 case binary_format:
|
|
|
1075 fp->format.category = sequence_data;
|
|
|
1076 fp->format.format = bam;
|
|
|
1077 /* fall-through */
|
|
|
1078 case bam:
|
|
|
1079 return bam_write1(fp->fp.bgzf, b);
|
|
|
1080
|
|
|
1081 case cram:
|
|
|
1082 return cram_put_bam_seq(fp->fp.cram, (bam1_t *)b);
|
|
|
1083
|
|
|
1084 case text_format:
|
|
|
1085 fp->format.category = sequence_data;
|
|
|
1086 fp->format.format = sam;
|
|
|
1087 /* fall-through */
|
|
|
1088 case sam:
|
|
|
1089 if (sam_format1(h, b, &fp->line) < 0) return -1;
|
|
|
1090 kputc('\n', &fp->line);
|
|
|
1091 if ( hwrite(fp->fp.hfile, fp->line.s, fp->line.l) != fp->line.l ) return -1;
|
|
|
1092 return fp->line.l;
|
|
|
1093
|
|
|
1094 default:
|
|
|
1095 abort();
|
|
|
1096 }
|
|
|
1097 }
|
|
|
1098
|
|
|
1099 /************************
|
|
|
1100 *** Auxiliary fields ***
|
|
|
1101 ************************/
|
|
|
1102
|
|
|
1103 void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data)
|
|
|
1104 {
|
|
|
1105 int ori_len = b->l_data;
|
|
|
1106 b->l_data += 3 + len;
|
|
|
1107 if (b->m_data < b->l_data) {
|
|
|
1108 b->m_data = b->l_data;
|
|
|
1109 kroundup32(b->m_data);
|
|
|
1110 b->data = (uint8_t*)realloc(b->data, b->m_data);
|
|
|
1111 }
|
|
|
1112 b->data[ori_len] = tag[0]; b->data[ori_len + 1] = tag[1];
|
|
|
1113 b->data[ori_len + 2] = type;
|
|
|
1114 memcpy(b->data + ori_len + 3, data, len);
|
|
|
1115 }
|
|
|
1116
|
|
|
1117 static inline uint8_t *skip_aux(uint8_t *s)
|
|
|
1118 {
|
|
|
1119 int size = aux_type2size(*s); ++s; // skip type
|
|
|
1120 uint32_t n;
|
|
|
1121 switch (size) {
|
|
|
1122 case 'Z':
|
|
|
1123 case 'H':
|
|
|
1124 while (*s) ++s;
|
|
|
1125 return s + 1;
|
|
|
1126 case 'B':
|
|
|
1127 size = aux_type2size(*s); ++s;
|
|
|
1128 memcpy(&n, s, 4); s += 4;
|
|
|
1129 return s + size * n;
|
|
|
1130 case 0:
|
|
|
1131 abort();
|
|
|
1132 break;
|
|
|
1133 default:
|
|
|
1134 return s + size;
|
|
|
1135 }
|
|
|
1136 }
|
|
|
1137
|
|
|
1138 uint8_t *bam_aux_get(const bam1_t *b, const char tag[2])
|
|
|
1139 {
|
|
|
1140 uint8_t *s;
|
|
|
1141 int y = tag[0]<<8 | tag[1];
|
|
|
1142 s = bam_get_aux(b);
|
|
|
1143 while (s < b->data + b->l_data) {
|
|
|
1144 int x = (int)s[0]<<8 | s[1];
|
|
|
1145 s += 2;
|
|
|
1146 if (x == y) return s;
|
|
|
1147 s = skip_aux(s);
|
|
|
1148 }
|
|
|
1149 return 0;
|
|
|
1150 }
|
|
|
1151 // s MUST BE returned by bam_aux_get()
|
|
|
1152 int bam_aux_del(bam1_t *b, uint8_t *s)
|
|
|
1153 {
|
|
|
1154 uint8_t *p, *aux;
|
|
|
1155 int l_aux = bam_get_l_aux(b);
|
|
|
1156 aux = bam_get_aux(b);
|
|
|
1157 p = s - 2;
|
|
|
1158 s = skip_aux(s);
|
|
|
1159 memmove(p, s, l_aux - (s - aux));
|
|
|
1160 b->l_data -= s - p;
|
|
|
1161 return 0;
|
|
|
1162 }
|
|
|
1163
|
|
|
1164 int32_t bam_aux2i(const uint8_t *s)
|
|
|
1165 {
|
|
|
1166 int type;
|
|
|
1167 type = *s++;
|
|
|
1168 if (type == 'c') return (int32_t)*(int8_t*)s;
|
|
|
1169 else if (type == 'C') return (int32_t)*(uint8_t*)s;
|
|
|
1170 else if (type == 's') return (int32_t)*(int16_t*)s;
|
|
|
1171 else if (type == 'S') return (int32_t)*(uint16_t*)s;
|
|
|
1172 else if (type == 'i' || type == 'I') return *(int32_t*)s;
|
|
|
1173 else return 0;
|
|
|
1174 }
|
|
|
1175
|
|
|
1176 double bam_aux2f(const uint8_t *s)
|
|
|
1177 {
|
|
|
1178 int type;
|
|
|
1179 type = *s++;
|
|
|
1180 if (type == 'd') return *(double*)s;
|
|
|
1181 else if (type == 'f') return *(float*)s;
|
|
|
1182 else return 0.0;
|
|
|
1183 }
|
|
|
1184
|
|
|
1185 char bam_aux2A(const uint8_t *s)
|
|
|
1186 {
|
|
|
1187 int type;
|
|
|
1188 type = *s++;
|
|
|
1189 if (type == 'A') return *(char*)s;
|
|
|
1190 else return 0;
|
|
|
1191 }
|
|
|
1192
|
|
|
1193 char *bam_aux2Z(const uint8_t *s)
|
|
|
1194 {
|
|
|
1195 int type;
|
|
|
1196 type = *s++;
|
|
|
1197 if (type == 'Z' || type == 'H') return (char*)s;
|
|
|
1198 else return 0;
|
|
|
1199 }
|
|
|
1200
|
|
|
1201 int sam_open_mode(char *mode, const char *fn, const char *format)
|
|
|
1202 {
|
|
|
1203 // TODO Parse "bam5" etc for compression level
|
|
|
1204 if (format == NULL) {
|
|
|
1205 // Try to pick a format based on the filename extension
|
|
|
1206 const char *ext = fn? strrchr(fn, '.') : NULL;
|
|
|
1207 if (ext == NULL || strchr(ext, '/')) return -1;
|
|
|
1208 return sam_open_mode(mode, fn, ext+1);
|
|
|
1209 }
|
|
|
1210 else if (strcmp(format, "bam") == 0) strcpy(mode, "b");
|
|
|
1211 else if (strcmp(format, "cram") == 0) strcpy(mode, "c");
|
|
|
1212 else if (strcmp(format, "sam") == 0) strcpy(mode, "");
|
|
|
1213 else return -1;
|
|
|
1214
|
|
|
1215 return 0;
|
|
|
1216 }
|
|
|
1217
|
|
|
1218 #define STRNCMP(a,b,n) (strncasecmp((a),(b),(n)) || strlen(a)!=(n))
|
|
|
1219 int bam_str2flag(const char *str)
|
|
|
1220 {
|
|
|
1221 char *end, *beg = (char*) str;
|
|
|
1222 long int flag = strtol(str, &end, 0);
|
|
|
1223 if ( end!=str ) return flag; // the conversion was successful
|
|
|
1224 flag = 0;
|
|
|
1225 while ( *str )
|
|
|
1226 {
|
|
|
1227 end = beg;
|
|
|
1228 while ( *end && *end!=',' ) end++;
|
|
|
1229 if ( !STRNCMP("PAIRED",beg,end-beg) ) flag |= BAM_FPAIRED;
|
|
|
1230 else if ( !STRNCMP("PROPER_PAIR",beg,end-beg) ) flag |= BAM_FPROPER_PAIR;
|
|
|
1231 else if ( !STRNCMP("UNMAP",beg,end-beg) ) flag |= BAM_FUNMAP;
|
|
|
1232 else if ( !STRNCMP("MUNMAP",beg,end-beg) ) flag |= BAM_FMUNMAP;
|
|
|
1233 else if ( !STRNCMP("REVERSE",beg,end-beg) ) flag |= BAM_FREVERSE;
|
|
|
1234 else if ( !STRNCMP("MREVERSE",beg,end-beg) ) flag |= BAM_FMREVERSE;
|
|
|
1235 else if ( !STRNCMP("READ1",beg,end-beg) ) flag |= BAM_FREAD1;
|
|
|
1236 else if ( !STRNCMP("READ2",beg,end-beg) ) flag |= BAM_FREAD2;
|
|
|
1237 else if ( !STRNCMP("SECONDARY",beg,end-beg) ) flag |= BAM_FSECONDARY;
|
|
|
1238 else if ( !STRNCMP("QCFAIL",beg,end-beg) ) flag |= BAM_FQCFAIL;
|
|
|
1239 else if ( !STRNCMP("DUP",beg,end-beg) ) flag |= BAM_FDUP;
|
|
|
1240 else if ( !STRNCMP("SUPPLEMENTARY",beg,end-beg) ) flag |= BAM_FSUPPLEMENTARY;
|
|
|
1241 else return -1;
|
|
|
1242 if ( !*end ) break;
|
|
|
1243 beg = end + 1;
|
|
|
1244 }
|
|
|
1245 return flag;
|
|
|
1246 }
|
|
|
1247
|
|
|
1248 char *bam_flag2str(int flag)
|
|
|
1249 {
|
|
|
1250 kstring_t str = {0,0,0};
|
|
|
1251 if ( flag&BAM_FPAIRED ) ksprintf(&str,"%s%s", str.l?",":"","PAIRED");
|
|
|
1252 if ( flag&BAM_FPROPER_PAIR ) ksprintf(&str,"%s%s", str.l?",":"","PROPER_PAIR");
|
|
|
1253 if ( flag&BAM_FUNMAP ) ksprintf(&str,"%s%s", str.l?",":"","UNMAP");
|
|
|
1254 if ( flag&BAM_FMUNMAP ) ksprintf(&str,"%s%s", str.l?",":"","MUNMAP");
|
|
|
1255 if ( flag&BAM_FREVERSE ) ksprintf(&str,"%s%s", str.l?",":"","REVERSE");
|
|
|
1256 if ( flag&BAM_FMREVERSE ) ksprintf(&str,"%s%s", str.l?",":"","MREVERSE");
|
|
|
1257 if ( flag&BAM_FREAD1 ) ksprintf(&str,"%s%s", str.l?",":"","READ1");
|
|
|
1258 if ( flag&BAM_FREAD2 ) ksprintf(&str,"%s%s", str.l?",":"","READ2");
|
|
|
1259 if ( flag&BAM_FSECONDARY ) ksprintf(&str,"%s%s", str.l?",":"","SECONDARY");
|
|
|
1260 if ( flag&BAM_FQCFAIL ) ksprintf(&str,"%s%s", str.l?",":"","QCFAIL");
|
|
|
1261 if ( flag&BAM_FDUP ) ksprintf(&str,"%s%s", str.l?",":"","DUP");
|
|
|
1262 if ( flag&BAM_FSUPPLEMENTARY ) ksprintf(&str,"%s%s", str.l?",":"","SUPPLEMENTARY");
|
|
|
1263 if ( str.l == 0 ) kputsn("", 0, &str);
|
|
|
1264 return str.s;
|
|
|
1265 }
|
|
|
1266
|
|
|
1267
|
|
|
1268 /**************************
|
|
|
1269 *** Pileup and Mpileup ***
|
|
|
1270 **************************/
|
|
|
1271
|
|
|
1272 #if !defined(BAM_NO_PILEUP)
|
|
|
1273
|
|
|
1274 #include <assert.h>
|
|
|
1275
|
|
|
1276 /*******************
|
|
|
1277 *** Memory pool ***
|
|
|
1278 *******************/
|
|
|
1279
|
|
|
1280 typedef struct {
|
|
|
1281 int k, x, y, end;
|
|
|
1282 } cstate_t;
|
|
|
1283
|
|
|
1284 static cstate_t g_cstate_null = { -1, 0, 0, 0 };
|
|
|
1285
|
|
|
1286 typedef struct __linkbuf_t {
|
|
|
1287 bam1_t b;
|
|
|
1288 int32_t beg, end;
|
|
|
1289 cstate_t s;
|
|
|
1290 struct __linkbuf_t *next;
|
|
|
1291 } lbnode_t;
|
|
|
1292
|
|
|
1293 typedef struct {
|
|
|
1294 int cnt, n, max;
|
|
|
1295 lbnode_t **buf;
|
|
|
1296 } mempool_t;
|
|
|
1297
|
|
|
1298 static mempool_t *mp_init(void)
|
|
|
1299 {
|
|
|
1300 mempool_t *mp;
|
|
|
1301 mp = (mempool_t*)calloc(1, sizeof(mempool_t));
|
|
|
1302 return mp;
|
|
|
1303 }
|
|
|
1304 static void mp_destroy(mempool_t *mp)
|
|
|
1305 {
|
|
|
1306 int k;
|
|
|
1307 for (k = 0; k < mp->n; ++k) {
|
|
|
1308 free(mp->buf[k]->b.data);
|
|
|
1309 free(mp->buf[k]);
|
|
|
1310 }
|
|
|
1311 free(mp->buf);
|
|
|
1312 free(mp);
|
|
|
1313 }
|
|
|
1314 static inline lbnode_t *mp_alloc(mempool_t *mp)
|
|
|
1315 {
|
|
|
1316 ++mp->cnt;
|
|
|
1317 if (mp->n == 0) return (lbnode_t*)calloc(1, sizeof(lbnode_t));
|
|
|
1318 else return mp->buf[--mp->n];
|
|
|
1319 }
|
|
|
1320 static inline void mp_free(mempool_t *mp, lbnode_t *p)
|
|
|
1321 {
|
|
|
1322 --mp->cnt; p->next = 0; // clear lbnode_t::next here
|
|
|
1323 if (mp->n == mp->max) {
|
|
|
1324 mp->max = mp->max? mp->max<<1 : 256;
|
|
|
1325 mp->buf = (lbnode_t**)realloc(mp->buf, sizeof(lbnode_t*) * mp->max);
|
|
|
1326 }
|
|
|
1327 mp->buf[mp->n++] = p;
|
|
|
1328 }
|
|
|
1329
|
|
|
1330 /**********************
|
|
|
1331 *** CIGAR resolver ***
|
|
|
1332 **********************/
|
|
|
1333
|
|
|
1334 /* s->k: the index of the CIGAR operator that has just been processed.
|
|
|
1335 s->x: the reference coordinate of the start of s->k
|
|
|
1336 s->y: the query coordiante of the start of s->k
|
|
|
1337 */
|
|
|
1338 static inline int resolve_cigar2(bam_pileup1_t *p, int32_t pos, cstate_t *s)
|
|
|
1339 {
|
|
|
1340 #define _cop(c) ((c)&BAM_CIGAR_MASK)
|
|
|
1341 #define _cln(c) ((c)>>BAM_CIGAR_SHIFT)
|
|
|
1342
|
|
|
1343 bam1_t *b = p->b;
|
|
|
1344 bam1_core_t *c = &b->core;
|
|
|
1345 uint32_t *cigar = bam_get_cigar(b);
|
|
|
1346 int k;
|
|
|
1347 // determine the current CIGAR operation
|
|
|
1348 // fprintf(stderr, "%s\tpos=%d\tend=%d\t(%d,%d,%d)\n", bam_get_qname(b), pos, s->end, s->k, s->x, s->y);
|
|
|
1349 if (s->k == -1) { // never processed
|
|
|
1350 if (c->n_cigar == 1) { // just one operation, save a loop
|
|
|
1351 if (_cop(cigar[0]) == BAM_CMATCH || _cop(cigar[0]) == BAM_CEQUAL || _cop(cigar[0]) == BAM_CDIFF) s->k = 0, s->x = c->pos, s->y = 0;
|
|
|
1352 } else { // find the first match or deletion
|
|
|
1353 for (k = 0, s->x = c->pos, s->y = 0; k < c->n_cigar; ++k) {
|
|
|
1354 int op = _cop(cigar[k]);
|
|
|
1355 int l = _cln(cigar[k]);
|
|
|
1356 if (op == BAM_CMATCH || op == BAM_CDEL || op == BAM_CEQUAL || op == BAM_CDIFF) break;
|
|
|
1357 else if (op == BAM_CREF_SKIP) s->x += l;
|
|
|
1358 else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) s->y += l;
|
|
|
1359 }
|
|
|
1360 assert(k < c->n_cigar);
|
|
|
1361 s->k = k;
|
|
|
1362 }
|
|
|
1363 } else { // the read has been processed before
|
|
|
1364 int op, l = _cln(cigar[s->k]);
|
|
|
1365 if (pos - s->x >= l) { // jump to the next operation
|
|
|
1366 assert(s->k < c->n_cigar); // otherwise a bug: this function should not be called in this case
|
|
|
1367 op = _cop(cigar[s->k+1]);
|
|
|
1368 if (op == BAM_CMATCH || op == BAM_CDEL || op == BAM_CREF_SKIP || op == BAM_CEQUAL || op == BAM_CDIFF) { // jump to the next without a loop
|
|
|
1369 if (_cop(cigar[s->k]) == BAM_CMATCH|| _cop(cigar[s->k]) == BAM_CEQUAL || _cop(cigar[s->k]) == BAM_CDIFF) s->y += l;
|
|
|
1370 s->x += l;
|
|
|
1371 ++s->k;
|
|
|
1372 } else { // find the next M/D/N/=/X
|
|
|
1373 if (_cop(cigar[s->k]) == BAM_CMATCH|| _cop(cigar[s->k]) == BAM_CEQUAL || _cop(cigar[s->k]) == BAM_CDIFF) s->y += l;
|
|
|
1374 s->x += l;
|
|
|
1375 for (k = s->k + 1; k < c->n_cigar; ++k) {
|
|
|
1376 op = _cop(cigar[k]), l = _cln(cigar[k]);
|
|
|
1377 if (op == BAM_CMATCH || op == BAM_CDEL || op == BAM_CREF_SKIP || op == BAM_CEQUAL || op == BAM_CDIFF) break;
|
|
|
1378 else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) s->y += l;
|
|
|
1379 }
|
|
|
1380 s->k = k;
|
|
|
1381 }
|
|
|
1382 assert(s->k < c->n_cigar); // otherwise a bug
|
|
|
1383 } // else, do nothing
|
|
|
1384 }
|
|
|
1385 { // collect pileup information
|
|
|
1386 int op, l;
|
|
|
1387 op = _cop(cigar[s->k]); l = _cln(cigar[s->k]);
|
|
|
1388 p->is_del = p->indel = p->is_refskip = 0;
|
|
|
1389 if (s->x + l - 1 == pos && s->k + 1 < c->n_cigar) { // peek the next operation
|
|
|
1390 int op2 = _cop(cigar[s->k+1]);
|
|
|
1391 int l2 = _cln(cigar[s->k+1]);
|
|
|
1392 if (op2 == BAM_CDEL) p->indel = -(int)l2;
|
|
|
1393 else if (op2 == BAM_CINS) p->indel = l2;
|
|
|
1394 else if (op2 == BAM_CPAD && s->k + 2 < c->n_cigar) { // no working for adjacent padding
|
|
|
1395 int l3 = 0;
|
|
|
1396 for (k = s->k + 2; k < c->n_cigar; ++k) {
|
|
|
1397 op2 = _cop(cigar[k]); l2 = _cln(cigar[k]);
|
|
|
1398 if (op2 == BAM_CINS) l3 += l2;
|
|
|
1399 else if (op2 == BAM_CDEL || op2 == BAM_CMATCH || op2 == BAM_CREF_SKIP || op2 == BAM_CEQUAL || op2 == BAM_CDIFF) break;
|
|
|
1400 }
|
|
|
1401 if (l3 > 0) p->indel = l3;
|
|
|
1402 }
|
|
|
1403 }
|
|
|
1404 if (op == BAM_CMATCH || op == BAM_CEQUAL || op == BAM_CDIFF) {
|
|
|
1405 p->qpos = s->y + (pos - s->x);
|
|
|
1406 } else if (op == BAM_CDEL || op == BAM_CREF_SKIP) {
|
|
|
1407 p->is_del = 1; p->qpos = s->y; // FIXME: distinguish D and N!!!!!
|
|
|
1408 p->is_refskip = (op == BAM_CREF_SKIP);
|
|
|
1409 } // cannot be other operations; otherwise a bug
|
|
|
1410 p->is_head = (pos == c->pos); p->is_tail = (pos == s->end);
|
|
|
1411 }
|
|
|
1412 return 1;
|
|
|
1413 }
|
|
|
1414
|
|
|
1415 /***********************
|
|
|
1416 *** Pileup iterator ***
|
|
|
1417 ***********************/
|
|
|
1418
|
|
|
1419 // Dictionary of overlapping reads
|
|
|
1420 KHASH_MAP_INIT_STR(olap_hash, lbnode_t *)
|
|
|
1421 typedef khash_t(olap_hash) olap_hash_t;
|
|
|
1422
|
|
|
1423 struct __bam_plp_t {
|
|
|
1424 mempool_t *mp;
|
|
|
1425 lbnode_t *head, *tail, *dummy;
|
|
|
1426 int32_t tid, pos, max_tid, max_pos;
|
|
|
1427 int is_eof, max_plp, error, maxcnt;
|
|
|
1428 uint64_t id;
|
|
|
1429 bam_pileup1_t *plp;
|
|
|
1430 // for the "auto" interface only
|
|
|
1431 bam1_t *b;
|
|
|
1432 bam_plp_auto_f func;
|
|
|
1433 void *data;
|
|
|
1434 olap_hash_t *overlaps;
|
|
|
1435 };
|
|
|
1436
|
|
|
1437 bam_plp_t bam_plp_init(bam_plp_auto_f func, void *data)
|
|
|
1438 {
|
|
|
1439 bam_plp_t iter;
|
|
|
1440 iter = (bam_plp_t)calloc(1, sizeof(struct __bam_plp_t));
|
|
|
1441 iter->mp = mp_init();
|
|
|
1442 iter->head = iter->tail = mp_alloc(iter->mp);
|
|
|
1443 iter->dummy = mp_alloc(iter->mp);
|
|
|
1444 iter->max_tid = iter->max_pos = -1;
|
|
|
1445 iter->maxcnt = 8000;
|
|
|
1446 if (func) {
|
|
|
1447 iter->func = func;
|
|
|
1448 iter->data = data;
|
|
|
1449 iter->b = bam_init1();
|
|
|
1450 }
|
|
|
1451 return iter;
|
|
|
1452 }
|
|
|
1453
|
|
|
1454 void bam_plp_init_overlaps(bam_plp_t iter)
|
|
|
1455 {
|
|
|
1456 iter->overlaps = kh_init(olap_hash); // hash for tweaking quality of bases in overlapping reads
|
|
|
1457 }
|
|
|
1458
|
|
|
1459 void bam_plp_destroy(bam_plp_t iter)
|
|
|
1460 {
|
|
|
1461 if ( iter->overlaps ) kh_destroy(olap_hash, iter->overlaps);
|
|
|
1462 mp_free(iter->mp, iter->dummy);
|
|
|
1463 mp_free(iter->mp, iter->head);
|
|
|
1464 if (iter->mp->cnt != 0)
|
|
|
1465 fprintf(stderr, "[bam_plp_destroy] memory leak: %d. Continue anyway.\n", iter->mp->cnt);
|
|
|
1466 mp_destroy(iter->mp);
|
|
|
1467 if (iter->b) bam_destroy1(iter->b);
|
|
|
1468 free(iter->plp);
|
|
|
1469 free(iter);
|
|
|
1470 }
|
|
|
1471
|
|
|
1472
|
|
|
1473 //---------------------------------
|
|
|
1474 //--- Tweak overlapping reads
|
|
|
1475 //---------------------------------
|
|
|
1476
|
|
|
1477 /**
|
|
|
1478 * cigar_iref2iseq_set() - find the first CMATCH setting the ref and the read index
|
|
|
1479 * cigar_iref2iseq_next() - get the next CMATCH base
|
|
|
1480 * @cigar: pointer to current cigar block (rw)
|
|
|
1481 * @cigar_max: pointer just beyond the last cigar block
|
|
|
1482 * @icig: position within the current cigar block (rw)
|
|
|
1483 * @iseq: position in the sequence (rw)
|
|
|
1484 * @iref: position with respect to the beginning of the read (iref_pos - b->core.pos) (rw)
|
|
|
1485 *
|
|
|
1486 * Returns BAM_CMATCH or -1 when there is no more cigar to process or the requested position is not covered.
|
|
|
1487 */
|
|
|
1488 static inline int cigar_iref2iseq_set(uint32_t **cigar, uint32_t *cigar_max, int *icig, int *iseq, int *iref)
|
|
|
1489 {
|
|
|
1490 int pos = *iref;
|
|
|
1491 if ( pos < 0 ) return -1;
|
|
|
1492 *icig = 0;
|
|
|
1493 *iseq = 0;
|
|
|
1494 *iref = 0;
|
|
|
1495 while ( *cigar<cigar_max )
|
|
|
1496 {
|
|
|
1497 int cig = (**cigar) & BAM_CIGAR_MASK;
|
|
|
1498 int ncig = (**cigar) >> BAM_CIGAR_SHIFT;
|
|
|
1499
|
|
|
1500 if ( cig==BAM_CSOFT_CLIP ) { (*cigar)++; *iseq += ncig; *icig = 0; continue; }
|
|
|
1501 if ( cig==BAM_CHARD_CLIP || cig==BAM_CPAD ) { (*cigar)++; *icig = 0; continue; }
|
|
|
1502 if ( cig==BAM_CMATCH || cig==BAM_CEQUAL || cig==BAM_CDIFF )
|
|
|
1503 {
|
|
|
1504 pos -= ncig;
|
|
|
1505 if ( pos < 0 ) { *icig = ncig + pos; *iseq += *icig; *iref += *icig; return BAM_CMATCH; }
|
|
|
1506 (*cigar)++; *iseq += ncig; *icig = 0; *iref += ncig;
|
|
|
1507 continue;
|
|
|
1508 }
|
|
|
1509 if ( cig==BAM_CINS ) { (*cigar)++; *iseq += ncig; *icig = 0; continue; }
|
|
|
1510 if ( cig==BAM_CDEL || cig==BAM_CREF_SKIP )
|
|
|
1511 {
|
|
|
1512 pos -= ncig;
|
|
|
1513 if ( pos<0 ) pos = 0;
|
|
|
1514 (*cigar)++; *icig = 0; *iref += ncig;
|
|
|
1515 continue;
|
|
|
1516 }
|
|
|
1517 fprintf(stderr,"todo: cigar %d\n", cig);
|
|
|
1518 assert(0);
|
|
|
1519 }
|
|
|
1520 *iseq = -1;
|
|
|
1521 return -1;
|
|
|
1522 }
|
|
|
1523 static inline int cigar_iref2iseq_next(uint32_t **cigar, uint32_t *cigar_max, int *icig, int *iseq, int *iref)
|
|
|
1524 {
|
|
|
1525 while ( *cigar < cigar_max )
|
|
|
1526 {
|
|
|
1527 int cig = (**cigar) & BAM_CIGAR_MASK;
|
|
|
1528 int ncig = (**cigar) >> BAM_CIGAR_SHIFT;
|
|
|
1529
|
|
|
1530 if ( cig==BAM_CMATCH || cig==BAM_CEQUAL || cig==BAM_CDIFF )
|
|
|
1531 {
|
|
|
1532 if ( *icig >= ncig - 1 ) { *icig = 0; (*cigar)++; continue; }
|
|
|
1533 (*iseq)++; (*icig)++; (*iref)++;
|
|
|
1534 return BAM_CMATCH;
|
|
|
1535 }
|
|
|
1536 if ( cig==BAM_CDEL || cig==BAM_CREF_SKIP ) { (*cigar)++; (*iref) += ncig; *icig = 0; continue; }
|
|
|
1537 if ( cig==BAM_CINS ) { (*cigar)++; *iseq += ncig; *icig = 0; continue; }
|
|
|
1538 if ( cig==BAM_CSOFT_CLIP ) { (*cigar)++; *iseq += ncig; *icig = 0; continue; }
|
|
|
1539 if ( cig==BAM_CHARD_CLIP || cig==BAM_CPAD ) { (*cigar)++; *icig = 0; continue; }
|
|
|
1540 fprintf(stderr,"todo: cigar %d\n", cig);
|
|
|
1541 assert(0);
|
|
|
1542 }
|
|
|
1543 *iseq = -1;
|
|
|
1544 *iref = -1;
|
|
|
1545 return -1;
|
|
|
1546 }
|
|
|
1547
|
|
|
1548 static void tweak_overlap_quality(bam1_t *a, bam1_t *b)
|
|
|
1549 {
|
|
|
1550 uint32_t *a_cigar = bam_get_cigar(a), *a_cigar_max = a_cigar + a->core.n_cigar;
|
|
|
1551 uint32_t *b_cigar = bam_get_cigar(b), *b_cigar_max = b_cigar + b->core.n_cigar;
|
|
|
1552 int a_icig = 0, a_iseq = 0;
|
|
|
1553 int b_icig = 0, b_iseq = 0;
|
|
|
1554 uint8_t *a_qual = bam_get_qual(a), *b_qual = bam_get_qual(b);
|
|
|
1555 uint8_t *a_seq = bam_get_seq(a), *b_seq = bam_get_seq(b);
|
|
|
1556
|
|
|
1557 int iref = b->core.pos;
|
|
|
1558 int a_iref = iref - a->core.pos;
|
|
|
1559 int b_iref = iref - b->core.pos;
|
|
|
1560 int a_ret = cigar_iref2iseq_set(&a_cigar, a_cigar_max, &a_icig, &a_iseq, &a_iref);
|
|
|
1561 if ( a_ret<0 ) return; // no overlap
|
|
|
1562 int b_ret = cigar_iref2iseq_set(&b_cigar, b_cigar_max, &b_icig, &b_iseq, &b_iref);
|
|
|
1563 if ( b_ret<0 ) return; // no overlap
|
|
|
1564
|
|
|
1565 #if DBG
|
|
|
1566 fprintf(stderr,"tweak %s n_cigar=%d %d .. %d-%d vs %d-%d\n", bam_get_qname(a), a->core.n_cigar, b->core.n_cigar,
|
|
|
1567 a->core.pos+1,a->core.pos+bam_cigar2rlen(a->core.n_cigar,bam_get_cigar(a)), b->core.pos+1, b->core.pos+bam_cigar2rlen(b->core.n_cigar,bam_get_cigar(b)));
|
|
|
1568 #endif
|
|
|
1569
|
|
|
1570 while ( 1 )
|
|
|
1571 {
|
|
|
1572 // Increment reference position
|
|
|
1573 while ( a_iref>=0 && a_iref < iref - a->core.pos )
|
|
|
1574 a_ret = cigar_iref2iseq_next(&a_cigar, a_cigar_max, &a_icig, &a_iseq, &a_iref);
|
|
|
1575 if ( a_ret<0 ) break; // done
|
|
|
1576 if ( iref < a_iref + a->core.pos ) iref = a_iref + a->core.pos;
|
|
|
1577
|
|
|
1578 while ( b_iref>=0 && b_iref < iref - b->core.pos )
|
|
|
1579 b_ret = cigar_iref2iseq_next(&b_cigar, b_cigar_max, &b_icig, &b_iseq, &b_iref);
|
|
|
1580 if ( b_ret<0 ) break; // done
|
|
|
1581 if ( iref < b_iref + b->core.pos ) iref = b_iref + b->core.pos;
|
|
|
1582
|
|
|
1583 iref++;
|
|
|
1584 if ( a_iref+a->core.pos != b_iref+b->core.pos ) continue; // only CMATCH positions, don't know what to do with indels
|
|
|
1585
|
|
|
1586 if ( bam_seqi(a_seq,a_iseq) == bam_seqi(b_seq,b_iseq) )
|
|
|
1587 {
|
|
|
1588 #if DBG
|
|
|
1589 fprintf(stderr,"%c",seq_nt16_str[bam_seqi(a_seq,a_iseq)]);
|
|
|
1590 #endif
|
|
|
1591 // we are very confident about this base
|
|
|
1592 int qual = a_qual[a_iseq] + b_qual[b_iseq];
|
|
|
1593 a_qual[a_iseq] = qual>200 ? 200 : qual;
|
|
|
1594 b_qual[b_iseq] = 0;
|
|
|
1595 }
|
|
|
1596 else
|
|
|
1597 {
|
|
|
1598 if ( a_qual[a_iseq] >= b_qual[b_iseq] )
|
|
|
1599 {
|
|
|
1600 #if DBG
|
|
|
1601 fprintf(stderr,"[%c/%c]",seq_nt16_str[bam_seqi(a_seq,a_iseq)],tolower(seq_nt16_str[bam_seqi(b_seq,b_iseq)]));
|
|
|
1602 #endif
|
|
|
1603 a_qual[a_iseq] = 0.8 * a_qual[a_iseq]; // not so confident about a_qual anymore given the mismatch
|
|
|
1604 b_qual[b_iseq] = 0;
|
|
|
1605 }
|
|
|
1606 else
|
|
|
1607 {
|
|
|
1608 #if DBG
|
|
|
1609 fprintf(stderr,"[%c/%c]",tolower(seq_nt16_str[bam_seqi(a_seq,a_iseq)]),seq_nt16_str[bam_seqi(b_seq,b_iseq)]);
|
|
|
1610 #endif
|
|
|
1611 b_qual[b_iseq] = 0.8 * b_qual[b_iseq];
|
|
|
1612 a_qual[a_iseq] = 0;
|
|
|
1613 }
|
|
|
1614 }
|
|
|
1615 }
|
|
|
1616 #if DBG
|
|
|
1617 fprintf(stderr,"\n");
|
|
|
1618 #endif
|
|
|
1619 }
|
|
|
1620
|
|
|
1621 // Fix overlapping reads. Simple soft-clipping did not give good results.
|
|
|
1622 // Lowering qualities of unwanted bases is more selective and works better.
|
|
|
1623 //
|
|
|
1624 static void overlap_push(bam_plp_t iter, lbnode_t *node)
|
|
|
1625 {
|
|
|
1626 if ( !iter->overlaps ) return;
|
|
|
1627
|
|
|
1628 // mapped mates and paired reads only
|
|
|
1629 if ( node->b.core.flag&BAM_FMUNMAP || !(node->b.core.flag&BAM_FPROPER_PAIR) ) return;
|
|
|
1630
|
|
|
1631 // no overlap possible, unless some wild cigar
|
|
|
1632 if ( abs(node->b.core.isize) >= 2*node->b.core.l_qseq ) return;
|
|
|
1633
|
|
|
1634 khiter_t kitr = kh_get(olap_hash, iter->overlaps, bam_get_qname(&node->b));
|
|
|
1635 if ( kitr==kh_end(iter->overlaps) )
|
|
|
1636 {
|
|
|
1637 int ret;
|
|
|
1638 kitr = kh_put(olap_hash, iter->overlaps, bam_get_qname(&node->b), &ret);
|
|
|
1639 kh_value(iter->overlaps, kitr) = node;
|
|
|
1640 }
|
|
|
1641 else
|
|
|
1642 {
|
|
|
1643 lbnode_t *a = kh_value(iter->overlaps, kitr);
|
|
|
1644 tweak_overlap_quality(&a->b, &node->b);
|
|
|
1645 kh_del(olap_hash, iter->overlaps, kitr);
|
|
|
1646 assert(a->end-1 == a->s.end);
|
|
|
1647 a->end = a->b.core.pos + bam_cigar2rlen(a->b.core.n_cigar, bam_get_cigar(&a->b));
|
|
|
1648 a->s.end = a->end - 1;
|
|
|
1649 }
|
|
|
1650 }
|
|
|
1651
|
|
|
1652 static void overlap_remove(bam_plp_t iter, const bam1_t *b)
|
|
|
1653 {
|
|
|
1654 if ( !iter->overlaps ) return;
|
|
|
1655
|
|
|
1656 khiter_t kitr;
|
|
|
1657 if ( b )
|
|
|
1658 {
|
|
|
1659 kitr = kh_get(olap_hash, iter->overlaps, bam_get_qname(b));
|
|
|
1660 if ( kitr!=kh_end(iter->overlaps) )
|
|
|
1661 kh_del(olap_hash, iter->overlaps, kitr);
|
|
|
1662 }
|
|
|
1663 else
|
|
|
1664 {
|
|
|
1665 // remove all
|
|
|
1666 for (kitr = kh_begin(iter->overlaps); kitr<kh_end(iter->overlaps); kitr++)
|
|
|
1667 if ( kh_exist(iter->overlaps, kitr) ) kh_del(olap_hash, iter->overlaps, kitr);
|
|
|
1668 }
|
|
|
1669 }
|
|
|
1670
|
|
|
1671
|
|
|
1672
|
|
|
1673 // Prepares next pileup position in bam records collected by bam_plp_auto -> user func -> bam_plp_push. Returns
|
|
|
1674 // pointer to the piled records if next position is ready or NULL if there is not enough records in the
|
|
|
1675 // buffer yet (the current position is still the maximum position across all buffered reads).
|
|
|
1676 const bam_pileup1_t *bam_plp_next(bam_plp_t iter, int *_tid, int *_pos, int *_n_plp)
|
|
|
1677 {
|
|
|
1678 if (iter->error) { *_n_plp = -1; return 0; }
|
|
|
1679 *_n_plp = 0;
|
|
|
1680 if (iter->is_eof && iter->head->next == 0) return 0;
|
|
|
1681 while (iter->is_eof || iter->max_tid > iter->tid || (iter->max_tid == iter->tid && iter->max_pos > iter->pos)) {
|
|
|
1682 int n_plp = 0;
|
|
|
1683 lbnode_t *p, *q;
|
|
|
1684 // write iter->plp at iter->pos
|
|
|
1685 iter->dummy->next = iter->head;
|
|
|
1686 for (p = iter->head, q = iter->dummy; p->next; q = p, p = p->next) {
|
|
|
1687 if (p->b.core.tid < iter->tid || (p->b.core.tid == iter->tid && p->end <= iter->pos)) { // then remove
|
|
|
1688 overlap_remove(iter, &p->b);
|
|
|
1689 q->next = p->next; mp_free(iter->mp, p); p = q;
|
|
|
1690 } else if (p->b.core.tid == iter->tid && p->beg <= iter->pos) { // here: p->end > pos; then add to pileup
|
|
|
1691 if (n_plp == iter->max_plp) { // then double the capacity
|
|
|
1692 iter->max_plp = iter->max_plp? iter->max_plp<<1 : 256;
|
|
|
1693 iter->plp = (bam_pileup1_t*)realloc(iter->plp, sizeof(bam_pileup1_t) * iter->max_plp);
|
|
|
1694 }
|
|
|
1695 iter->plp[n_plp].b = &p->b;
|
|
|
1696 if (resolve_cigar2(iter->plp + n_plp, iter->pos, &p->s)) ++n_plp; // actually always true...
|
|
|
1697 }
|
|
|
1698 }
|
|
|
1699 iter->head = iter->dummy->next; // dummy->next may be changed
|
|
|
1700 *_n_plp = n_plp; *_tid = iter->tid; *_pos = iter->pos;
|
|
|
1701 // update iter->tid and iter->pos
|
|
|
1702 if (iter->head->next) {
|
|
|
1703 if (iter->tid > iter->head->b.core.tid) {
|
|
|
1704 fprintf(stderr, "[%s] unsorted input. Pileup aborts.\n", __func__);
|
|
|
1705 iter->error = 1;
|
|
|
1706 *_n_plp = -1;
|
|
|
1707 return 0;
|
|
|
1708 }
|
|
|
1709 }
|
|
|
1710 if (iter->tid < iter->head->b.core.tid) { // come to a new reference sequence
|
|
|
1711 iter->tid = iter->head->b.core.tid; iter->pos = iter->head->beg; // jump to the next reference
|
|
|
1712 } else if (iter->pos < iter->head->beg) { // here: tid == head->b.core.tid
|
|
|
1713 iter->pos = iter->head->beg; // jump to the next position
|
|
|
1714 } else ++iter->pos; // scan contiguously
|
|
|
1715 // return
|
|
|
1716 if (n_plp) return iter->plp;
|
|
|
1717 if (iter->is_eof && iter->head->next == 0) break;
|
|
|
1718 }
|
|
|
1719 return 0;
|
|
|
1720 }
|
|
|
1721
|
|
|
1722 int bam_plp_push(bam_plp_t iter, const bam1_t *b)
|
|
|
1723 {
|
|
|
1724 if (iter->error) return -1;
|
|
|
1725 if (b) {
|
|
|
1726 if (b->core.tid < 0) { overlap_remove(iter, b); return 0; }
|
|
|
1727 // Skip only unmapped reads here, any additional filtering must be done in iter->func
|
|
|
1728 if (b->core.flag & BAM_FUNMAP) { overlap_remove(iter, b); return 0; }
|
|
|
1729 if (iter->tid == b->core.tid && iter->pos == b->core.pos && iter->mp->cnt > iter->maxcnt)
|
|
|
1730 {
|
|
|
1731 overlap_remove(iter, b);
|
|
|
1732 return 0;
|
|
|
1733 }
|
|
|
1734 bam_copy1(&iter->tail->b, b);
|
|
|
1735 overlap_push(iter, iter->tail);
|
|
|
1736 #ifndef BAM_NO_ID
|
|
|
1737 iter->tail->b.id = iter->id++;
|
|
|
1738 #endif
|
|
|
1739 iter->tail->beg = b->core.pos;
|
|
|
1740 iter->tail->end = b->core.pos + bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b));
|
|
|
1741 iter->tail->s = g_cstate_null; iter->tail->s.end = iter->tail->end - 1; // initialize cstate_t
|
|
|
1742 if (b->core.tid < iter->max_tid) {
|
|
|
1743 fprintf(stderr, "[bam_pileup_core] the input is not sorted (chromosomes out of order)\n");
|
|
|
1744 iter->error = 1;
|
|
|
1745 return -1;
|
|
|
1746 }
|
|
|
1747 if ((b->core.tid == iter->max_tid) && (iter->tail->beg < iter->max_pos)) {
|
|
|
1748 fprintf(stderr, "[bam_pileup_core] the input is not sorted (reads out of order)\n");
|
|
|
1749 iter->error = 1;
|
|
|
1750 return -1;
|
|
|
1751 }
|
|
|
1752 iter->max_tid = b->core.tid; iter->max_pos = iter->tail->beg;
|
|
|
1753 if (iter->tail->end > iter->pos || iter->tail->b.core.tid > iter->tid) {
|
|
|
1754 iter->tail->next = mp_alloc(iter->mp);
|
|
|
1755 iter->tail = iter->tail->next;
|
|
|
1756 }
|
|
|
1757 } else iter->is_eof = 1;
|
|
|
1758 return 0;
|
|
|
1759 }
|
|
|
1760
|
|
|
1761 const bam_pileup1_t *bam_plp_auto(bam_plp_t iter, int *_tid, int *_pos, int *_n_plp)
|
|
|
1762 {
|
|
|
1763 const bam_pileup1_t *plp;
|
|
|
1764 if (iter->func == 0 || iter->error) { *_n_plp = -1; return 0; }
|
|
|
1765 if ((plp = bam_plp_next(iter, _tid, _pos, _n_plp)) != 0) return plp;
|
|
|
1766 else { // no pileup line can be obtained; read alignments
|
|
|
1767 *_n_plp = 0;
|
|
|
1768 if (iter->is_eof) return 0;
|
|
|
1769 int ret;
|
|
|
1770 while ( (ret=iter->func(iter->data, iter->b)) >= 0) {
|
|
|
1771 if (bam_plp_push(iter, iter->b) < 0) {
|
|
|
1772 *_n_plp = -1;
|
|
|
1773 return 0;
|
|
|
1774 }
|
|
|
1775 if ((plp = bam_plp_next(iter, _tid, _pos, _n_plp)) != 0) return plp;
|
|
|
1776 // otherwise no pileup line can be returned; read the next alignment.
|
|
|
1777 }
|
|
|
1778 if ( ret < -1 ) { iter->error = ret; *_n_plp = -1; return 0; }
|
|
|
1779 bam_plp_push(iter, 0);
|
|
|
1780 if ((plp = bam_plp_next(iter, _tid, _pos, _n_plp)) != 0) return plp;
|
|
|
1781 return 0;
|
|
|
1782 }
|
|
|
1783 }
|
|
|
1784
|
|
|
1785 void bam_plp_reset(bam_plp_t iter)
|
|
|
1786 {
|
|
|
1787 lbnode_t *p, *q;
|
|
|
1788 iter->max_tid = iter->max_pos = -1;
|
|
|
1789 iter->tid = iter->pos = 0;
|
|
|
1790 iter->is_eof = 0;
|
|
|
1791 for (p = iter->head; p->next;) {
|
|
|
1792 overlap_remove(iter, NULL);
|
|
|
1793 q = p->next;
|
|
|
1794 mp_free(iter->mp, p);
|
|
|
1795 p = q;
|
|
|
1796 }
|
|
|
1797 iter->head = iter->tail;
|
|
|
1798 }
|
|
|
1799
|
|
|
1800 void bam_plp_set_maxcnt(bam_plp_t iter, int maxcnt)
|
|
|
1801 {
|
|
|
1802 iter->maxcnt = maxcnt;
|
|
|
1803 }
|
|
|
1804
|
|
|
1805 /************************
|
|
|
1806 *** Mpileup iterator ***
|
|
|
1807 ************************/
|
|
|
1808
|
|
|
1809 struct __bam_mplp_t {
|
|
|
1810 int n;
|
|
|
1811 uint64_t min, *pos;
|
|
|
1812 bam_plp_t *iter;
|
|
|
1813 int *n_plp;
|
|
|
1814 const bam_pileup1_t **plp;
|
|
|
1815 };
|
|
|
1816
|
|
|
1817 bam_mplp_t bam_mplp_init(int n, bam_plp_auto_f func, void **data)
|
|
|
1818 {
|
|
|
1819 int i;
|
|
|
1820 bam_mplp_t iter;
|
|
|
1821 iter = (bam_mplp_t)calloc(1, sizeof(struct __bam_mplp_t));
|
|
|
1822 iter->pos = (uint64_t*)calloc(n, sizeof(uint64_t));
|
|
|
1823 iter->n_plp = (int*)calloc(n, sizeof(int));
|
|
|
1824 iter->plp = (const bam_pileup1_t**)calloc(n, sizeof(bam_pileup1_t*));
|
|
|
1825 iter->iter = (bam_plp_t*)calloc(n, sizeof(bam_plp_t));
|
|
|
1826 iter->n = n;
|
|
|
1827 iter->min = (uint64_t)-1;
|
|
|
1828 for (i = 0; i < n; ++i) {
|
|
|
1829 iter->iter[i] = bam_plp_init(func, data[i]);
|
|
|
1830 iter->pos[i] = iter->min;
|
|
|
1831 }
|
|
|
1832 return iter;
|
|
|
1833 }
|
|
|
1834
|
|
|
1835 void bam_mplp_init_overlaps(bam_mplp_t iter)
|
|
|
1836 {
|
|
|
1837 int i;
|
|
|
1838 for (i = 0; i < iter->n; ++i)
|
|
|
1839 bam_plp_init_overlaps(iter->iter[i]);
|
|
|
1840 }
|
|
|
1841
|
|
|
1842 void bam_mplp_set_maxcnt(bam_mplp_t iter, int maxcnt)
|
|
|
1843 {
|
|
|
1844 int i;
|
|
|
1845 for (i = 0; i < iter->n; ++i)
|
|
|
1846 iter->iter[i]->maxcnt = maxcnt;
|
|
|
1847 }
|
|
|
1848
|
|
|
1849 void bam_mplp_destroy(bam_mplp_t iter)
|
|
|
1850 {
|
|
|
1851 int i;
|
|
|
1852 for (i = 0; i < iter->n; ++i) bam_plp_destroy(iter->iter[i]);
|
|
|
1853 free(iter->iter); free(iter->pos); free(iter->n_plp); free(iter->plp);
|
|
|
1854 free(iter);
|
|
|
1855 }
|
|
|
1856
|
|
|
1857 int bam_mplp_auto(bam_mplp_t iter, int *_tid, int *_pos, int *n_plp, const bam_pileup1_t **plp)
|
|
|
1858 {
|
|
|
1859 int i, ret = 0;
|
|
|
1860 uint64_t new_min = (uint64_t)-1;
|
|
|
1861 for (i = 0; i < iter->n; ++i) {
|
|
|
1862 if (iter->pos[i] == iter->min) {
|
|
|
1863 int tid, pos;
|
|
|
1864 iter->plp[i] = bam_plp_auto(iter->iter[i], &tid, &pos, &iter->n_plp[i]);
|
|
|
1865 if ( iter->iter[i]->error ) return -1;
|
|
|
1866 iter->pos[i] = iter->plp[i] ? (uint64_t)tid<<32 | pos : 0;
|
|
|
1867 }
|
|
|
1868 if (iter->plp[i] && iter->pos[i] < new_min) new_min = iter->pos[i];
|
|
|
1869 }
|
|
|
1870 iter->min = new_min;
|
|
|
1871 if (new_min == (uint64_t)-1) return 0;
|
|
|
1872 *_tid = new_min>>32; *_pos = (uint32_t)new_min;
|
|
|
1873 for (i = 0; i < iter->n; ++i) {
|
|
|
1874 if (iter->pos[i] == iter->min) { // FIXME: valgrind reports "uninitialised value(s) at this line"
|
|
|
1875 n_plp[i] = iter->n_plp[i], plp[i] = iter->plp[i];
|
|
|
1876 ++ret;
|
|
|
1877 } else n_plp[i] = 0, plp[i] = 0;
|
|
|
1878 }
|
|
|
1879 return ret;
|
|
|
1880 }
|
|
|
1881
|
|
|
1882 #endif // ~!defined(BAM_NO_PILEUP)
|
