Mercurial > repos > youngkim > ezbamqc
view ezBAMQC/src/htslib/sam.c @ 8:82bb8c455761
Uploaded
| author | cshl-bsr |
|---|---|
| date | Wed, 30 Mar 2016 12:03:10 -0400 |
| parents | dfa3745e5fd8 |
| children |
line wrap: on
line source
/* sam.c -- SAM and BAM file I/O and manipulation. Copyright (C) 2008-2010, 2012-2014 Genome Research Ltd. Copyright (C) 2010, 2012, 2013 Broad Institute. Author: Heng Li <lh3@sanger.ac.uk> Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #include <stdio.h> #include <stdlib.h> #include <string.h> #include <errno.h> #include <ctype.h> #include <zlib.h> #include "htslib/sam.h" #include "htslib/bgzf.h" #include "cram/cram.h" #include "htslib/hfile.h" #include "htslib/khash.h" KHASH_DECLARE(s2i, kh_cstr_t, int64_t) typedef khash_t(s2i) sdict_t; /********************** *** BAM header I/O *** **********************/ bam_hdr_t *bam_hdr_init() { return (bam_hdr_t*)calloc(1, sizeof(bam_hdr_t)); } void bam_hdr_destroy(bam_hdr_t *h) { int32_t i; if (h == NULL) return; if (h->target_name) { for (i = 0; i < h->n_targets; ++i) free(h->target_name[i]); free(h->target_name); free(h->target_len); } free(h->text); free(h->cigar_tab); if (h->sdict) kh_destroy(s2i, (sdict_t*)h->sdict); free(h); } bam_hdr_t *bam_hdr_dup(const bam_hdr_t *h0) { if (h0 == NULL) return NULL; bam_hdr_t *h; if ((h = bam_hdr_init()) == NULL) return NULL; // copy the simple data h->n_targets = h0->n_targets; h->ignore_sam_err = h0->ignore_sam_err; h->l_text = h0->l_text; // Then the pointery stuff h->cigar_tab = NULL; h->sdict = NULL; h->text = (char*)calloc(h->l_text + 1, 1); memcpy(h->text, h0->text, h->l_text); h->target_len = (uint32_t*)calloc(h->n_targets, sizeof(uint32_t)); h->target_name = (char**)calloc(h->n_targets, sizeof(char*)); int i; for (i = 0; i < h->n_targets; ++i) { h->target_len[i] = h0->target_len[i]; h->target_name[i] = strdup(h0->target_name[i]); } return h; } static bam_hdr_t *hdr_from_dict(sdict_t *d) { bam_hdr_t *h; khint_t k; h = bam_hdr_init(); h->sdict = d; h->n_targets = kh_size(d); h->target_len = (uint32_t*)malloc(sizeof(uint32_t) * h->n_targets); h->target_name = (char**)malloc(sizeof(char*) * h->n_targets); for (k = kh_begin(d); k != kh_end(d); ++k) { if (!kh_exist(d, k)) continue; h->target_name[kh_val(d, k)>>32] = (char*)kh_key(d, k); h->target_len[kh_val(d, k)>>32] = kh_val(d, k)<<32>>32; kh_val(d, k) >>= 32; } return h; } bam_hdr_t *bam_hdr_read(BGZF *fp) { bam_hdr_t *h; char buf[4]; int magic_len, has_EOF; int32_t i = 1, name_len; // check EOF has_EOF = bgzf_check_EOF(fp); if (has_EOF < 0) { perror("[W::sam_hdr_read] bgzf_check_EOF"); } else if (has_EOF == 0 && hts_verbose >= 2) fprintf(stderr, "[W::%s] EOF marker is absent. The input is probably truncated.\n", __func__); // read "BAM1" magic_len = bgzf_read(fp, buf, 4); if (magic_len != 4 || strncmp(buf, "BAM\1", 4)) { if (hts_verbose >= 1) fprintf(stderr, "[E::%s] invalid BAM binary header\n", __func__); return 0; } h = bam_hdr_init(); // read plain text and the number of reference sequences bgzf_read(fp, &h->l_text, 4); if (fp->is_be) ed_swap_4p(&h->l_text); h->text = (char*)malloc(h->l_text + 1); h->text[h->l_text] = 0; // make sure it is NULL terminated bgzf_read(fp, h->text, h->l_text); bgzf_read(fp, &h->n_targets, 4); if (fp->is_be) ed_swap_4p(&h->n_targets); // read reference sequence names and lengths h->target_name = (char**)calloc(h->n_targets, sizeof(char*)); h->target_len = (uint32_t*)calloc(h->n_targets, sizeof(uint32_t)); for (i = 0; i != h->n_targets; ++i) { bgzf_read(fp, &name_len, 4); if (fp->is_be) ed_swap_4p(&name_len); h->target_name[i] = (char*)calloc(name_len, 1); bgzf_read(fp, h->target_name[i], name_len); bgzf_read(fp, &h->target_len[i], 4); if (fp->is_be) ed_swap_4p(&h->target_len[i]); } return h; } int bam_hdr_write(BGZF *fp, const bam_hdr_t *h) { char buf[4]; int32_t i, name_len, x; // write "BAM1" strncpy(buf, "BAM\1", 4); bgzf_write(fp, buf, 4); // write plain text and the number of reference sequences if (fp->is_be) { x = ed_swap_4(h->l_text); bgzf_write(fp, &x, 4); if (h->l_text) bgzf_write(fp, h->text, h->l_text); x = ed_swap_4(h->n_targets); bgzf_write(fp, &x, 4); } else { bgzf_write(fp, &h->l_text, 4); if (h->l_text) bgzf_write(fp, h->text, h->l_text); bgzf_write(fp, &h->n_targets, 4); } // write sequence names and lengths for (i = 0; i != h->n_targets; ++i) { char *p = h->target_name[i]; name_len = strlen(p) + 1; if (fp->is_be) { x = ed_swap_4(name_len); bgzf_write(fp, &x, 4); } else bgzf_write(fp, &name_len, 4); bgzf_write(fp, p, name_len); if (fp->is_be) { x = ed_swap_4(h->target_len[i]); bgzf_write(fp, &x, 4); } else bgzf_write(fp, &h->target_len[i], 4); } bgzf_flush(fp); return 0; } int bam_name2id(bam_hdr_t *h, const char *ref) { sdict_t *d = (sdict_t*)h->sdict; khint_t k; if (h->sdict == 0) { int i, absent; d = kh_init(s2i); for (i = 0; i < h->n_targets; ++i) { k = kh_put(s2i, d, h->target_name[i], &absent); kh_val(d, k) = i; } h->sdict = d; } k = kh_get(s2i, d, ref); return k == kh_end(d)? -1 : kh_val(d, k); } /************************* *** BAM alignment I/O *** *************************/ bam1_t *bam_init1() { return (bam1_t*)calloc(1, sizeof(bam1_t)); } void bam_destroy1(bam1_t *b) { if (b == 0) return; free(b->data); free(b); } bam1_t *bam_copy1(bam1_t *bdst, const bam1_t *bsrc) { uint8_t *data = bdst->data; int m_data = bdst->m_data; // backup data and m_data if (m_data < bsrc->l_data) { // double the capacity m_data = bsrc->l_data; kroundup32(m_data); data = (uint8_t*)realloc(data, m_data); } memcpy(data, bsrc->data, bsrc->l_data); // copy var-len data *bdst = *bsrc; // copy the rest // restore the backup bdst->m_data = m_data; bdst->data = data; return bdst; } bam1_t *bam_dup1(const bam1_t *bsrc) { if (bsrc == NULL) return NULL; bam1_t *bdst = bam_init1(); if (bdst == NULL) return NULL; return bam_copy1(bdst, bsrc); } int bam_cigar2qlen(int n_cigar, const uint32_t *cigar) { int k, l; for (k = l = 0; k < n_cigar; ++k) if (bam_cigar_type(bam_cigar_op(cigar[k]))&1) l += bam_cigar_oplen(cigar[k]); return l; } int bam_cigar2rlen(int n_cigar, const uint32_t *cigar) { int k, l; for (k = l = 0; k < n_cigar; ++k) if (bam_cigar_type(bam_cigar_op(cigar[k]))&2) l += bam_cigar_oplen(cigar[k]); return l; } int32_t bam_endpos(const bam1_t *b) { if (!(b->core.flag & BAM_FUNMAP) && b->core.n_cigar > 0) return b->core.pos + bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b)); else return b->core.pos + 1; } static inline int aux_type2size(uint8_t type) { switch (type) { case 'A': case 'c': case 'C': return 1; case 's': case 'S': return 2; case 'i': case 'I': case 'f': return 4; case 'd': return 8; case 'Z': case 'H': case 'B': return type; default: return 0; } } static void swap_data(const bam1_core_t *c, int l_data, uint8_t *data, int is_host) { uint8_t *s; uint32_t *cigar = (uint32_t*)(data + c->l_qname); uint32_t i, n; s = data + c->n_cigar*4 + c->l_qname + c->l_qseq + (c->l_qseq + 1)/2; for (i = 0; i < c->n_cigar; ++i) ed_swap_4p(&cigar[i]); while (s < data + l_data) { int size; s += 2; // skip key size = aux_type2size(*s); ++s; // skip type switch (size) { case 1: ++s; break; case 2: ed_swap_2p(s); s += 2; break; case 4: ed_swap_4p(s); s += 4; break; case 8: ed_swap_8p(s); s += 8; break; case 'Z': case 'H': while (*s) ++s; ++s; break; case 'B': size = aux_type2size(*s); ++s; if (is_host) memcpy(&n, s, 4), ed_swap_4p(s); else ed_swap_4p(s), memcpy(&n, s, 4); s += 4; switch (size) { case 1: s += n; break; case 2: for (i = 0; i < n; ++i, s += 2) ed_swap_2p(s); break; case 4: for (i = 0; i < n; ++i, s += 4) ed_swap_4p(s); break; case 8: for (i = 0; i < n; ++i, s += 8) ed_swap_8p(s); break; } break; } } } int bam_read1(BGZF *fp, bam1_t *b) { bam1_core_t *c = &b->core; int32_t block_len, ret, i; uint32_t x[8]; if ((ret = bgzf_read(fp, &block_len, 4)) != 4) { if (ret == 0) return -1; // normal end-of-file else return -2; // truncated } if (bgzf_read(fp, x, 32) != 32) return -3; if (fp->is_be) { ed_swap_4p(&block_len); for (i = 0; i < 8; ++i) ed_swap_4p(x + i); } c->tid = x[0]; c->pos = x[1]; c->bin = x[2]>>16; c->qual = x[2]>>8&0xff; c->l_qname = x[2]&0xff; c->flag = x[3]>>16; c->n_cigar = x[3]&0xffff; c->l_qseq = x[4]; c->mtid = x[5]; c->mpos = x[6]; c->isize = x[7]; b->l_data = block_len - 32; if (b->l_data < 0 || c->l_qseq < 0) return -4; if ((char *)bam_get_aux(b) - (char *)b->data > b->l_data) return -4; if (b->m_data < b->l_data) { b->m_data = b->l_data; kroundup32(b->m_data); b->data = (uint8_t*)realloc(b->data, b->m_data); if (!b->data) return -4; } if (bgzf_read(fp, b->data, b->l_data) != b->l_data) return -4; //b->l_aux = b->l_data - c->n_cigar * 4 - c->l_qname - c->l_qseq - (c->l_qseq+1)/2; if (fp->is_be) swap_data(c, b->l_data, b->data, 0); return 4 + block_len; } int bam_write1(BGZF *fp, const bam1_t *b) { const bam1_core_t *c = &b->core; uint32_t x[8], block_len = b->l_data + 32, y; int i, ok; x[0] = c->tid; x[1] = c->pos; x[2] = (uint32_t)c->bin<<16 | c->qual<<8 | c->l_qname; x[3] = (uint32_t)c->flag<<16 | c->n_cigar; x[4] = c->l_qseq; x[5] = c->mtid; x[6] = c->mpos; x[7] = c->isize; ok = (bgzf_flush_try(fp, 4 + block_len) >= 0); if (fp->is_be) { for (i = 0; i < 8; ++i) ed_swap_4p(x + i); y = block_len; if (ok) ok = (bgzf_write(fp, ed_swap_4p(&y), 4) >= 0); swap_data(c, b->l_data, b->data, 1); } else { if (ok) ok = (bgzf_write(fp, &block_len, 4) >= 0); } if (ok) ok = (bgzf_write(fp, x, 32) >= 0); if (ok) ok = (bgzf_write(fp, b->data, b->l_data) >= 0); if (fp->is_be) swap_data(c, b->l_data, b->data, 0); return ok? 4 + block_len : -1; } /******************** *** BAM indexing *** ********************/ static hts_idx_t *bam_index(BGZF *fp, int min_shift) { int n_lvls, i, fmt; bam1_t *b; hts_idx_t *idx; bam_hdr_t *h; h = bam_hdr_read(fp); if (min_shift > 0) { int64_t max_len = 0, s; for (i = 0; i < h->n_targets; ++i) if (max_len < h->target_len[i]) max_len = h->target_len[i]; max_len += 256; for (n_lvls = 0, s = 1<<min_shift; max_len > s; ++n_lvls, s <<= 3); fmt = HTS_FMT_CSI; } else min_shift = 14, n_lvls = 5, fmt = HTS_FMT_BAI; idx = hts_idx_init(h->n_targets, fmt, bgzf_tell(fp), min_shift, n_lvls); bam_hdr_destroy(h); b = bam_init1(); while (bam_read1(fp, b) >= 0) { int l, ret; l = bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b)); if (l == 0) l = 1; // no zero-length records ret = hts_idx_push(idx, b->core.tid, b->core.pos, b->core.pos + l, bgzf_tell(fp), !(b->core.flag&BAM_FUNMAP)); if (ret < 0) { // unsorted bam_destroy1(b); hts_idx_destroy(idx); return NULL; } } hts_idx_finish(idx, bgzf_tell(fp)); bam_destroy1(b); return idx; } int bam_index_build(const char *fn, int min_shift) { hts_idx_t *idx; htsFile *fp; int ret = 0; if ((fp = hts_open(fn, "r")) == 0) return -1; switch (fp->format.format) { case cram: ret = cram_index_build(fp->fp.cram, fn); break; case bam: idx = bam_index(fp->fp.bgzf, min_shift); if (idx) { hts_idx_save(idx, fn, (min_shift > 0)? HTS_FMT_CSI : HTS_FMT_BAI); hts_idx_destroy(idx); } else ret = -1; break; default: ret = -1; break; } hts_close(fp); return ret; } static int bam_readrec(BGZF *fp, void *ignored, void *bv, int *tid, int *beg, int *end) { bam1_t *b = bv; int ret; if ((ret = bam_read1(fp, b)) >= 0) { *tid = b->core.tid; *beg = b->core.pos; *end = b->core.pos + (b->core.n_cigar? bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b)) : 1); } return ret; } // This is used only with read_rest=1 iterators, so need not set tid/beg/end. static int cram_readrec(BGZF *ignored, void *fpv, void *bv, int *tid, int *beg, int *end) { htsFile *fp = fpv; bam1_t *b = bv; return cram_get_bam_seq(fp->fp.cram, &b); } // This is used only with read_rest=1 iterators, so need not set tid/beg/end. static int sam_bam_cram_readrec(BGZF *bgzfp, void *fpv, void *bv, int *tid, int *beg, int *end) { htsFile *fp = fpv; bam1_t *b = bv; switch (fp->format.format) { case bam: return bam_read1(bgzfp, b); case cram: return cram_get_bam_seq(fp->fp.cram, &b); default: // TODO Need headers available to implement this for SAM files fprintf(stderr, "[sam_bam_cram_readrec] Not implemented for SAM files -- Exiting\n"); abort(); } } // The CRAM implementation stores the loaded index within the cram_fd rather // than separately as is done elsewhere in htslib. So if p is a pointer to // an hts_idx_t with p->fmt == HTS_FMT_CRAI, then it actually points to an // hts_cram_idx_t and should be cast accordingly. typedef struct hts_cram_idx_t { int fmt; cram_fd *cram; } hts_cram_idx_t; hts_idx_t *sam_index_load(samFile *fp, const char *fn) { switch (fp->format.format) { case bam: return bam_index_load(fn); case cram: { if (cram_index_load(fp->fp.cram, fn) < 0) return NULL; // Cons up a fake "index" just pointing at the associated cram_fd: hts_cram_idx_t *idx = malloc(sizeof (hts_cram_idx_t)); if (idx == NULL) return NULL; idx->fmt = HTS_FMT_CRAI; idx->cram = fp->fp.cram; return (hts_idx_t *) idx; } default: return NULL; // TODO Would use tbx_index_load if it returned hts_idx_t } } static hts_itr_t *cram_itr_query(const hts_idx_t *idx, int tid, int beg, int end, hts_readrec_func *readrec) { const hts_cram_idx_t *cidx = (const hts_cram_idx_t *) idx; hts_itr_t *iter = (hts_itr_t *) calloc(1, sizeof(hts_itr_t)); if (iter == NULL) return NULL; // Cons up a dummy iterator for which hts_itr_next() will simply invoke // the readrec function: iter->read_rest = 1; iter->off = NULL; iter->bins.a = NULL; iter->readrec = readrec; if (tid >= 0) { cram_range r = { tid, beg+1, end }; if (cram_set_option(cidx->cram, CRAM_OPT_RANGE, &r) != 0) { free(iter); return NULL; } iter->curr_off = 0; // The following fields are not required by hts_itr_next(), but are // filled in in case user code wants to look at them. iter->tid = tid; iter->beg = beg; iter->end = end; } else switch (tid) { case HTS_IDX_REST: iter->curr_off = 0; break; case HTS_IDX_NONE: iter->curr_off = 0; iter->finished = 1; break; default: fprintf(stderr, "[cram_itr_query] tid=%d not implemented for CRAM files -- Exiting\n", tid); abort(); break; } return iter; } hts_itr_t *sam_itr_queryi(const hts_idx_t *idx, int tid, int beg, int end) { const hts_cram_idx_t *cidx = (const hts_cram_idx_t *) idx; if (idx == NULL) return hts_itr_query(NULL, tid, beg, end, sam_bam_cram_readrec); else if (cidx->fmt == HTS_FMT_CRAI) return cram_itr_query(idx, tid, beg, end, cram_readrec); else return hts_itr_query(idx, tid, beg, end, bam_readrec); } static int cram_name2id(void *fdv, const char *ref) { cram_fd *fd = (cram_fd *) fdv; return sam_hdr_name2ref(fd->header, ref); } hts_itr_t *sam_itr_querys(const hts_idx_t *idx, bam_hdr_t *hdr, const char *region) { const hts_cram_idx_t *cidx = (const hts_cram_idx_t *) idx; if (cidx->fmt == HTS_FMT_CRAI) return hts_itr_querys(idx, region, cram_name2id, cidx->cram, cram_itr_query, cram_readrec); else return hts_itr_querys(idx, region, (hts_name2id_f)(bam_name2id), hdr, hts_itr_query, bam_readrec); } /********************** *** SAM header I/O *** **********************/ #include "htslib/kseq.h" #include "htslib/kstring.h" bam_hdr_t *sam_hdr_parse(int l_text, const char *text) { const char *q, *r, *p; khash_t(s2i) *d; d = kh_init(s2i); for (p = text; *p; ++p) { if (strncmp(p, "@SQ", 3) == 0) { char *sn = 0; int ln = -1; for (q = p + 4;; ++q) { if (strncmp(q, "SN:", 3) == 0) { q += 3; for (r = q; *r != '\t' && *r != '\n'; ++r); sn = (char*)calloc(r - q + 1, 1); strncpy(sn, q, r - q); q = r; } else if (strncmp(q, "LN:", 3) == 0) ln = strtol(q + 3, (char**)&q, 10); while (*q != '\t' && *q != '\n') ++q; if (*q == '\n') break; } p = q; if (sn && ln >= 0) { khint_t k; int absent; k = kh_put(s2i, d, sn, &absent); if (!absent) { if (hts_verbose >= 2) fprintf(stderr, "[W::%s] duplicated sequence '%s'\n", __func__, sn); free(sn); } else kh_val(d, k) = (int64_t)(kh_size(d) - 1)<<32 | ln; } } while (*p != '\n') ++p; } return hdr_from_dict(d); } bam_hdr_t *sam_hdr_read(htsFile *fp) { switch (fp->format.format) { case bam: return bam_hdr_read(fp->fp.bgzf); case cram: return cram_header_to_bam(fp->fp.cram->header); case sam: { kstring_t str; bam_hdr_t *h; int has_SQ = 0; str.l = str.m = 0; str.s = 0; while (hts_getline(fp, KS_SEP_LINE, &fp->line) >= 0) { if (fp->line.s[0] != '@') break; if (fp->line.l > 3 && strncmp(fp->line.s,"@SQ",3) == 0) has_SQ = 1; kputsn(fp->line.s, fp->line.l, &str); kputc('\n', &str); } if (! has_SQ && fp->fn_aux) { char line[2048]; FILE *f = fopen(fp->fn_aux, "r"); if (f == NULL) return NULL; while (fgets(line, sizeof line, f)) { const char *name = strtok(line, "\t"); const char *length = strtok(NULL, "\t"); ksprintf(&str, "@SQ\tSN:%s\tLN:%s\n", name, length); } fclose(f); } if (str.l == 0) kputsn("", 0, &str); h = sam_hdr_parse(str.l, str.s); h->l_text = str.l; h->text = str.s; return h; } default: abort(); } } int sam_hdr_write(htsFile *fp, const bam_hdr_t *h) { switch (fp->format.format) { case binary_format: fp->format.category = sequence_data; fp->format.format = bam; /* fall-through */ case bam: bam_hdr_write(fp->fp.bgzf, h); break; case cram: { cram_fd *fd = fp->fp.cram; if (cram_set_header(fd, bam_header_to_cram((bam_hdr_t *)h)) < 0) return -1; if (fp->fn_aux) cram_load_reference(fd, fp->fn_aux); if (cram_write_SAM_hdr(fd, fd->header) < 0) return -1; } break; case text_format: fp->format.category = sequence_data; fp->format.format = sam; /* fall-through */ case sam: { char *p; hputs(h->text, fp->fp.hfile); p = strstr(h->text, "@SQ\t"); // FIXME: we need a loop to make sure "@SQ\t" does not match something unwanted!!! if (p == 0) { int i; for (i = 0; i < h->n_targets; ++i) { fp->line.l = 0; kputsn("@SQ\tSN:", 7, &fp->line); kputs(h->target_name[i], &fp->line); kputsn("\tLN:", 4, &fp->line); kputw(h->target_len[i], &fp->line); kputc('\n', &fp->line); if ( hwrite(fp->fp.hfile, fp->line.s, fp->line.l) != fp->line.l ) return -1; } } if ( hflush(fp->fp.hfile) != 0 ) return -1; } break; default: abort(); } return 0; } /********************** *** SAM record I/O *** **********************/ int sam_parse1(kstring_t *s, bam_hdr_t *h, bam1_t *b) { #define _read_token(_p) (_p); for (; *(_p) && *(_p) != '\t'; ++(_p)); if (*(_p) != '\t') goto err_ret; *(_p)++ = 0 #define _read_token_aux(_p) (_p); for (; *(_p) && *(_p) != '\t'; ++(_p)); *(_p)++ = 0 // this is different in that it does not test *(_p)=='\t' #define _get_mem(type_t, _x, _s, _l) ks_resize((_s), (_s)->l + (_l)); *(_x) = (type_t*)((_s)->s + (_s)->l); (_s)->l += (_l) #define _parse_err(cond, msg) do { if ((cond) && hts_verbose >= 1) { fprintf(stderr, "[E::%s] " msg "\n", __func__); goto err_ret; } } while (0) #define _parse_warn(cond, msg) if ((cond) && hts_verbose >= 2) fprintf(stderr, "[W::%s] " msg "\n", __func__) uint8_t *t; char *p = s->s, *q; int i; kstring_t str; bam1_core_t *c = &b->core; str.l = b->l_data = 0; str.s = (char*)b->data; str.m = b->m_data; memset(c, 0, 32); if (h->cigar_tab == 0) { h->cigar_tab = (int8_t*) malloc(128); for (i = 0; i < 128; ++i) h->cigar_tab[i] = -1; for (i = 0; BAM_CIGAR_STR[i]; ++i) h->cigar_tab[(int)BAM_CIGAR_STR[i]] = i; } // qname q = _read_token(p); kputsn_(q, p - q, &str); c->l_qname = p - q; // flag c->flag = strtol(p, &p, 0); if (*p++ != '\t') goto err_ret; // malformated flag // chr q = _read_token(p); if (strcmp(q, "*")) { _parse_err(h->n_targets == 0, "missing SAM header"); c->tid = bam_name2id(h, q); _parse_warn(c->tid < 0, "urecognized reference name; treated as unmapped"); } else c->tid = -1; // pos c->pos = strtol(p, &p, 10) - 1; if (*p++ != '\t') goto err_ret; if (c->pos < 0 && c->tid >= 0) { _parse_warn(1, "mapped query cannot have zero coordinate; treated as unmapped"); c->tid = -1; } if (c->tid < 0) c->flag |= BAM_FUNMAP; // mapq c->qual = strtol(p, &p, 10); if (*p++ != '\t') goto err_ret; // cigar if (*p != '*') { uint32_t *cigar; size_t n_cigar = 0; for (q = p; *p && *p != '\t'; ++p) if (!isdigit(*p)) ++n_cigar; if (*p++ != '\t') goto err_ret; _parse_err(n_cigar >= 65536, "too many CIGAR operations"); c->n_cigar = n_cigar; _get_mem(uint32_t, &cigar, &str, c->n_cigar<<2); for (i = 0; i < c->n_cigar; ++i, ++q) { int op; cigar[i] = strtol(q, &q, 10)<<BAM_CIGAR_SHIFT; op = (uint8_t)*q >= 128? -1 : h->cigar_tab[(int)*q]; _parse_err(op < 0, "unrecognized CIGAR operator"); cigar[i] |= op; } i = bam_cigar2rlen(c->n_cigar, cigar); } else { _parse_warn(!(c->flag&BAM_FUNMAP), "mapped query must have a CIGAR; treated as unmapped"); c->flag |= BAM_FUNMAP; q = _read_token(p); i = 1; } c->bin = hts_reg2bin(c->pos, c->pos + i, 14, 5); // mate chr q = _read_token(p); if (strcmp(q, "=") == 0) c->mtid = c->tid; else if (strcmp(q, "*") == 0) c->mtid = -1; else c->mtid = bam_name2id(h, q); // mpos c->mpos = strtol(p, &p, 10) - 1; if (*p++ != '\t') goto err_ret; if (c->mpos < 0 && c->mtid >= 0) { _parse_warn(1, "mapped mate cannot have zero coordinate; treated as unmapped"); c->mtid = -1; } // tlen c->isize = strtol(p, &p, 10); if (*p++ != '\t') goto err_ret; // seq q = _read_token(p); if (strcmp(q, "*")) { c->l_qseq = p - q - 1; i = bam_cigar2qlen(c->n_cigar, (uint32_t*)(str.s + c->l_qname)); _parse_err(c->n_cigar && i != c->l_qseq, "CIGAR and query sequence are of different length"); i = (c->l_qseq + 1) >> 1; _get_mem(uint8_t, &t, &str, i); memset(t, 0, i); for (i = 0; i < c->l_qseq; ++i) t[i>>1] |= seq_nt16_table[(int)q[i]] << ((~i&1)<<2); } else c->l_qseq = 0; // qual q = _read_token_aux(p); _get_mem(uint8_t, &t, &str, c->l_qseq); if (strcmp(q, "*")) { _parse_err(p - q - 1 != c->l_qseq, "SEQ and QUAL are of different length"); for (i = 0; i < c->l_qseq; ++i) t[i] = q[i] - 33; } else memset(t, 0xff, c->l_qseq); // aux // Note that (like the bam1_core_t fields) this aux data in b->data is // stored in host endianness; so there is no byte swapping needed here. while (p < s->s + s->l) { uint8_t type; q = _read_token_aux(p); // FIXME: can be accelerated for long 'B' arrays _parse_err(p - q - 1 < 6, "incomplete aux field"); kputsn_(q, 2, &str); q += 3; type = *q++; ++q; // q points to value if (type == 'A' || type == 'a' || type == 'c' || type == 'C') { kputc_('A', &str); kputc_(*q, &str); } else if (type == 'i' || type == 'I') { if (*q == '-') { long x = strtol(q, &q, 10); if (x >= INT8_MIN) { kputc_('c', &str); kputc_(x, &str); } else if (x >= INT16_MIN) { int16_t y = x; kputc_('s', &str); kputsn_((char*)&y, 2, &str); } else { int32_t y = x; kputc_('i', &str); kputsn_(&y, 4, &str); } } else { unsigned long x = strtoul(q, &q, 10); if (x <= UINT8_MAX) { kputc_('C', &str); kputc_(x, &str); } else if (x <= UINT16_MAX) { uint16_t y = x; kputc_('S', &str); kputsn_(&y, 2, &str); } else { uint32_t y = x; kputc_('I', &str); kputsn_(&y, 4, &str); } } } else if (type == 'f') { float x; x = strtod(q, &q); kputc_('f', &str); kputsn_(&x, 4, &str); } else if (type == 'd') { double x; x = strtod(q, &q); kputc_('d', &str); kputsn_(&x, 8, &str); } else if (type == 'Z' || type == 'H') { kputc_(type, &str);kputsn_(q, p - q, &str); // note that this include the trailing NULL } else if (type == 'B') { int32_t n; char *r; _parse_err(p - q - 1 < 3, "incomplete B-typed aux field"); type = *q++; // q points to the first ',' following the typing byte for (r = q, n = 0; *r; ++r) if (*r == ',') ++n; kputc_('B', &str); kputc_(type, &str); kputsn_(&n, 4, &str); // FIXME: to evaluate which is faster: a) aligned array and then memmove(); b) unaligned array; c) kputsn_() if (type == 'c') while (q + 1 < p) { int8_t x = strtol(q + 1, &q, 0); kputc_(x, &str); } else if (type == 'C') while (q + 1 < p) { uint8_t x = strtoul(q + 1, &q, 0); kputc_(x, &str); } else if (type == 's') while (q + 1 < p) { int16_t x = strtol(q + 1, &q, 0); kputsn_(&x, 2, &str); } else if (type == 'S') while (q + 1 < p) { uint16_t x = strtoul(q + 1, &q, 0); kputsn_(&x, 2, &str); } else if (type == 'i') while (q + 1 < p) { int32_t x = strtol(q + 1, &q, 0); kputsn_(&x, 4, &str); } else if (type == 'I') while (q + 1 < p) { uint32_t x = strtoul(q + 1, &q, 0); kputsn_(&x, 4, &str); } else if (type == 'f') while (q + 1 < p) { float x = strtod(q + 1, &q); kputsn_(&x, 4, &str); } else _parse_err(1, "unrecognized type"); } else _parse_err(1, "unrecognized type"); } b->data = (uint8_t*)str.s; b->l_data = str.l; b->m_data = str.m; return 0; #undef _parse_warn #undef _parse_err #undef _get_mem #undef _read_token_aux #undef _read_token err_ret: b->data = (uint8_t*)str.s; b->l_data = str.l; b->m_data = str.m; return -2; } int sam_read1(htsFile *fp, bam_hdr_t *h, bam1_t *b) { switch (fp->format.format) { case bam: { int r = bam_read1(fp->fp.bgzf, b); if (r >= 0) { if (b->core.tid >= h->n_targets || b->core.tid < -1 || b->core.mtid >= h->n_targets || b->core.mtid < -1) return -3; } return r; } case cram: return cram_get_bam_seq(fp->fp.cram, &b); case sam: { int ret; err_recover: if (fp->line.l == 0) { ret = hts_getline(fp, KS_SEP_LINE, &fp->line); if (ret < 0) return -1; } ret = sam_parse1(&fp->line, h, b); fp->line.l = 0; if (ret < 0) { if (hts_verbose >= 1) fprintf(stderr, "[W::%s] parse error at line %lld\n", __func__, (long long)fp->lineno); if (h->ignore_sam_err) goto err_recover; } return ret; } default: abort(); } } int sam_format1(const bam_hdr_t *h, const bam1_t *b, kstring_t *str) { int i; uint8_t *s; const bam1_core_t *c = &b->core; str->l = 0; kputsn(bam_get_qname(b), c->l_qname-1, str); kputc('\t', str); // query name kputw(c->flag, str); kputc('\t', str); // flag if (c->tid >= 0) { // chr kputs(h->target_name[c->tid] , str); kputc('\t', str); } else kputsn("*\t", 2, str); kputw(c->pos + 1, str); kputc('\t', str); // pos kputw(c->qual, str); kputc('\t', str); // qual if (c->n_cigar) { // cigar uint32_t *cigar = bam_get_cigar(b); for (i = 0; i < c->n_cigar; ++i) { kputw(bam_cigar_oplen(cigar[i]), str); kputc(bam_cigar_opchr(cigar[i]), str); } } else kputc('*', str); kputc('\t', str); if (c->mtid < 0) kputsn("*\t", 2, str); // mate chr else if (c->mtid == c->tid) kputsn("=\t", 2, str); else { kputs(h->target_name[c->mtid], str); kputc('\t', str); } kputw(c->mpos + 1, str); kputc('\t', str); // mate pos kputw(c->isize, str); kputc('\t', str); // template len if (c->l_qseq) { // seq and qual uint8_t *s = bam_get_seq(b); for (i = 0; i < c->l_qseq; ++i) kputc("=ACMGRSVTWYHKDBN"[bam_seqi(s, i)], str); kputc('\t', str); s = bam_get_qual(b); if (s[0] == 0xff) kputc('*', str); else for (i = 0; i < c->l_qseq; ++i) kputc(s[i] + 33, str); } else kputsn("*\t*", 3, str); s = bam_get_aux(b); // aux while (s+4 <= b->data + b->l_data) { uint8_t type, key[2]; key[0] = s[0]; key[1] = s[1]; s += 2; type = *s++; kputc('\t', str); kputsn((char*)key, 2, str); kputc(':', str); if (type == 'A') { kputsn("A:", 2, str); kputc(*s, str); ++s; } else if (type == 'C') { kputsn("i:", 2, str); kputw(*s, str); ++s; } else if (type == 'c') { kputsn("i:", 2, str); kputw(*(int8_t*)s, str); ++s; } else if (type == 'S') { if (s+2 <= b->data + b->l_data) { kputsn("i:", 2, str); kputw(*(uint16_t*)s, str); s += 2; } else return -1; } else if (type == 's') { if (s+2 <= b->data + b->l_data) { kputsn("i:", 2, str); kputw(*(int16_t*)s, str); s += 2; } else return -1; } else if (type == 'I') { if (s+4 <= b->data + b->l_data) { kputsn("i:", 2, str); kputuw(*(uint32_t*)s, str); s += 4; } else return -1; } else if (type == 'i') { if (s+4 <= b->data + b->l_data) { kputsn("i:", 2, str); kputw(*(int32_t*)s, str); s += 4; } else return -1; } else if (type == 'f') { if (s+4 <= b->data + b->l_data) { ksprintf(str, "f:%g", *(float*)s); s += 4; } else return -1; } else if (type == 'd') { if (s+8 <= b->data + b->l_data) { ksprintf(str, "d:%g", *(double*)s); s += 8; } else return -1; } else if (type == 'Z' || type == 'H') { kputc(type, str); kputc(':', str); while (s < b->data + b->l_data && *s) kputc(*s++, str); if (s >= b->data + b->l_data) return -1; ++s; } else if (type == 'B') { uint8_t sub_type = *(s++); int32_t n; memcpy(&n, s, 4); s += 4; // no point to the start of the array if (s + n >= b->data + b->l_data) return -1; kputsn("B:", 2, str); kputc(sub_type, str); // write the typing for (i = 0; i < n; ++i) { // FIXME: for better performance, put the loop after "if" kputc(',', str); if ('c' == sub_type) { kputw(*(int8_t*)s, str); ++s; } else if ('C' == sub_type) { kputw(*(uint8_t*)s, str); ++s; } else if ('s' == sub_type) { kputw(*(int16_t*)s, str); s += 2; } else if ('S' == sub_type) { kputw(*(uint16_t*)s, str); s += 2; } else if ('i' == sub_type) { kputw(*(int32_t*)s, str); s += 4; } else if ('I' == sub_type) { kputuw(*(uint32_t*)s, str); s += 4; } else if ('f' == sub_type) { ksprintf(str, "%g", *(float*)s); s += 4; } } } } return str->l; } int sam_write1(htsFile *fp, const bam_hdr_t *h, const bam1_t *b) { switch (fp->format.format) { case binary_format: fp->format.category = sequence_data; fp->format.format = bam; /* fall-through */ case bam: return bam_write1(fp->fp.bgzf, b); case cram: return cram_put_bam_seq(fp->fp.cram, (bam1_t *)b); case text_format: fp->format.category = sequence_data; fp->format.format = sam; /* fall-through */ case sam: if (sam_format1(h, b, &fp->line) < 0) return -1; kputc('\n', &fp->line); if ( hwrite(fp->fp.hfile, fp->line.s, fp->line.l) != fp->line.l ) return -1; return fp->line.l; default: abort(); } } /************************ *** Auxiliary fields *** ************************/ void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data) { int ori_len = b->l_data; b->l_data += 3 + len; if (b->m_data < b->l_data) { b->m_data = b->l_data; kroundup32(b->m_data); b->data = (uint8_t*)realloc(b->data, b->m_data); } b->data[ori_len] = tag[0]; b->data[ori_len + 1] = tag[1]; b->data[ori_len + 2] = type; memcpy(b->data + ori_len + 3, data, len); } static inline uint8_t *skip_aux(uint8_t *s) { int size = aux_type2size(*s); ++s; // skip type uint32_t n; switch (size) { case 'Z': case 'H': while (*s) ++s; return s + 1; case 'B': size = aux_type2size(*s); ++s; memcpy(&n, s, 4); s += 4; return s + size * n; case 0: abort(); break; default: return s + size; } } uint8_t *bam_aux_get(const bam1_t *b, const char tag[2]) { uint8_t *s; int y = tag[0]<<8 | tag[1]; s = bam_get_aux(b); while (s < b->data + b->l_data) { int x = (int)s[0]<<8 | s[1]; s += 2; if (x == y) return s; s = skip_aux(s); } return 0; } // s MUST BE returned by bam_aux_get() int bam_aux_del(bam1_t *b, uint8_t *s) { uint8_t *p, *aux; int l_aux = bam_get_l_aux(b); aux = bam_get_aux(b); p = s - 2; s = skip_aux(s); memmove(p, s, l_aux - (s - aux)); b->l_data -= s - p; return 0; } int32_t bam_aux2i(const uint8_t *s) { int type; type = *s++; if (type == 'c') return (int32_t)*(int8_t*)s; else if (type == 'C') return (int32_t)*(uint8_t*)s; else if (type == 's') return (int32_t)*(int16_t*)s; else if (type == 'S') return (int32_t)*(uint16_t*)s; else if (type == 'i' || type == 'I') return *(int32_t*)s; else return 0; } double bam_aux2f(const uint8_t *s) { int type; type = *s++; if (type == 'd') return *(double*)s; else if (type == 'f') return *(float*)s; else return 0.0; } char bam_aux2A(const uint8_t *s) { int type; type = *s++; if (type == 'A') return *(char*)s; else return 0; } char *bam_aux2Z(const uint8_t *s) { int type; type = *s++; if (type == 'Z' || type == 'H') return (char*)s; else return 0; } int sam_open_mode(char *mode, const char *fn, const char *format) { // TODO Parse "bam5" etc for compression level if (format == NULL) { // Try to pick a format based on the filename extension const char *ext = fn? strrchr(fn, '.') : NULL; if (ext == NULL || strchr(ext, '/')) return -1; return sam_open_mode(mode, fn, ext+1); } else if (strcmp(format, "bam") == 0) strcpy(mode, "b"); else if (strcmp(format, "cram") == 0) strcpy(mode, "c"); else if (strcmp(format, "sam") == 0) strcpy(mode, ""); else return -1; return 0; } #define STRNCMP(a,b,n) (strncasecmp((a),(b),(n)) || strlen(a)!=(n)) int bam_str2flag(const char *str) { char *end, *beg = (char*) str; long int flag = strtol(str, &end, 0); if ( end!=str ) return flag; // the conversion was successful flag = 0; while ( *str ) { end = beg; while ( *end && *end!=',' ) end++; if ( !STRNCMP("PAIRED",beg,end-beg) ) flag |= BAM_FPAIRED; else if ( !STRNCMP("PROPER_PAIR",beg,end-beg) ) flag |= BAM_FPROPER_PAIR; else if ( !STRNCMP("UNMAP",beg,end-beg) ) flag |= BAM_FUNMAP; else if ( !STRNCMP("MUNMAP",beg,end-beg) ) flag |= BAM_FMUNMAP; else if ( !STRNCMP("REVERSE",beg,end-beg) ) flag |= BAM_FREVERSE; else if ( !STRNCMP("MREVERSE",beg,end-beg) ) flag |= BAM_FMREVERSE; else if ( !STRNCMP("READ1",beg,end-beg) ) flag |= BAM_FREAD1; else if ( !STRNCMP("READ2",beg,end-beg) ) flag |= BAM_FREAD2; else if ( !STRNCMP("SECONDARY",beg,end-beg) ) flag |= BAM_FSECONDARY; else if ( !STRNCMP("QCFAIL",beg,end-beg) ) flag |= BAM_FQCFAIL; else if ( !STRNCMP("DUP",beg,end-beg) ) flag |= BAM_FDUP; else if ( !STRNCMP("SUPPLEMENTARY",beg,end-beg) ) flag |= BAM_FSUPPLEMENTARY; else return -1; if ( !*end ) break; beg = end + 1; } return flag; } char *bam_flag2str(int flag) { kstring_t str = {0,0,0}; if ( flag&BAM_FPAIRED ) ksprintf(&str,"%s%s", str.l?",":"","PAIRED"); if ( flag&BAM_FPROPER_PAIR ) ksprintf(&str,"%s%s", str.l?",":"","PROPER_PAIR"); if ( flag&BAM_FUNMAP ) ksprintf(&str,"%s%s", str.l?",":"","UNMAP"); if ( flag&BAM_FMUNMAP ) ksprintf(&str,"%s%s", str.l?",":"","MUNMAP"); if ( flag&BAM_FREVERSE ) ksprintf(&str,"%s%s", str.l?",":"","REVERSE"); if ( flag&BAM_FMREVERSE ) ksprintf(&str,"%s%s", str.l?",":"","MREVERSE"); if ( flag&BAM_FREAD1 ) ksprintf(&str,"%s%s", str.l?",":"","READ1"); if ( flag&BAM_FREAD2 ) ksprintf(&str,"%s%s", str.l?",":"","READ2"); if ( flag&BAM_FSECONDARY ) ksprintf(&str,"%s%s", str.l?",":"","SECONDARY"); if ( flag&BAM_FQCFAIL ) ksprintf(&str,"%s%s", str.l?",":"","QCFAIL"); if ( flag&BAM_FDUP ) ksprintf(&str,"%s%s", str.l?",":"","DUP"); if ( flag&BAM_FSUPPLEMENTARY ) ksprintf(&str,"%s%s", str.l?",":"","SUPPLEMENTARY"); if ( str.l == 0 ) kputsn("", 0, &str); return str.s; } /************************** *** Pileup and Mpileup *** **************************/ #if !defined(BAM_NO_PILEUP) #include <assert.h> /******************* *** Memory pool *** *******************/ typedef struct { int k, x, y, end; } cstate_t; static cstate_t g_cstate_null = { -1, 0, 0, 0 }; typedef struct __linkbuf_t { bam1_t b; int32_t beg, end; cstate_t s; struct __linkbuf_t *next; } lbnode_t; typedef struct { int cnt, n, max; lbnode_t **buf; } mempool_t; static mempool_t *mp_init(void) { mempool_t *mp; mp = (mempool_t*)calloc(1, sizeof(mempool_t)); return mp; } static void mp_destroy(mempool_t *mp) { int k; for (k = 0; k < mp->n; ++k) { free(mp->buf[k]->b.data); free(mp->buf[k]); } free(mp->buf); free(mp); } static inline lbnode_t *mp_alloc(mempool_t *mp) { ++mp->cnt; if (mp->n == 0) return (lbnode_t*)calloc(1, sizeof(lbnode_t)); else return mp->buf[--mp->n]; } static inline void mp_free(mempool_t *mp, lbnode_t *p) { --mp->cnt; p->next = 0; // clear lbnode_t::next here if (mp->n == mp->max) { mp->max = mp->max? mp->max<<1 : 256; mp->buf = (lbnode_t**)realloc(mp->buf, sizeof(lbnode_t*) * mp->max); } mp->buf[mp->n++] = p; } /********************** *** CIGAR resolver *** **********************/ /* s->k: the index of the CIGAR operator that has just been processed. s->x: the reference coordinate of the start of s->k s->y: the query coordiante of the start of s->k */ static inline int resolve_cigar2(bam_pileup1_t *p, int32_t pos, cstate_t *s) { #define _cop(c) ((c)&BAM_CIGAR_MASK) #define _cln(c) ((c)>>BAM_CIGAR_SHIFT) bam1_t *b = p->b; bam1_core_t *c = &b->core; uint32_t *cigar = bam_get_cigar(b); int k; // determine the current CIGAR operation // 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); if (s->k == -1) { // never processed if (c->n_cigar == 1) { // just one operation, save a loop 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; } else { // find the first match or deletion for (k = 0, s->x = c->pos, s->y = 0; k < c->n_cigar; ++k) { int op = _cop(cigar[k]); int l = _cln(cigar[k]); if (op == BAM_CMATCH || op == BAM_CDEL || op == BAM_CEQUAL || op == BAM_CDIFF) break; else if (op == BAM_CREF_SKIP) s->x += l; else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) s->y += l; } assert(k < c->n_cigar); s->k = k; } } else { // the read has been processed before int op, l = _cln(cigar[s->k]); if (pos - s->x >= l) { // jump to the next operation assert(s->k < c->n_cigar); // otherwise a bug: this function should not be called in this case op = _cop(cigar[s->k+1]); if (op == BAM_CMATCH || op == BAM_CDEL || op == BAM_CREF_SKIP || op == BAM_CEQUAL || op == BAM_CDIFF) { // jump to the next without a loop if (_cop(cigar[s->k]) == BAM_CMATCH|| _cop(cigar[s->k]) == BAM_CEQUAL || _cop(cigar[s->k]) == BAM_CDIFF) s->y += l; s->x += l; ++s->k; } else { // find the next M/D/N/=/X if (_cop(cigar[s->k]) == BAM_CMATCH|| _cop(cigar[s->k]) == BAM_CEQUAL || _cop(cigar[s->k]) == BAM_CDIFF) s->y += l; s->x += l; for (k = s->k + 1; k < c->n_cigar; ++k) { op = _cop(cigar[k]), l = _cln(cigar[k]); if (op == BAM_CMATCH || op == BAM_CDEL || op == BAM_CREF_SKIP || op == BAM_CEQUAL || op == BAM_CDIFF) break; else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) s->y += l; } s->k = k; } assert(s->k < c->n_cigar); // otherwise a bug } // else, do nothing } { // collect pileup information int op, l; op = _cop(cigar[s->k]); l = _cln(cigar[s->k]); p->is_del = p->indel = p->is_refskip = 0; if (s->x + l - 1 == pos && s->k + 1 < c->n_cigar) { // peek the next operation int op2 = _cop(cigar[s->k+1]); int l2 = _cln(cigar[s->k+1]); if (op2 == BAM_CDEL) p->indel = -(int)l2; else if (op2 == BAM_CINS) p->indel = l2; else if (op2 == BAM_CPAD && s->k + 2 < c->n_cigar) { // no working for adjacent padding int l3 = 0; for (k = s->k + 2; k < c->n_cigar; ++k) { op2 = _cop(cigar[k]); l2 = _cln(cigar[k]); if (op2 == BAM_CINS) l3 += l2; else if (op2 == BAM_CDEL || op2 == BAM_CMATCH || op2 == BAM_CREF_SKIP || op2 == BAM_CEQUAL || op2 == BAM_CDIFF) break; } if (l3 > 0) p->indel = l3; } } if (op == BAM_CMATCH || op == BAM_CEQUAL || op == BAM_CDIFF) { p->qpos = s->y + (pos - s->x); } else if (op == BAM_CDEL || op == BAM_CREF_SKIP) { p->is_del = 1; p->qpos = s->y; // FIXME: distinguish D and N!!!!! p->is_refskip = (op == BAM_CREF_SKIP); } // cannot be other operations; otherwise a bug p->is_head = (pos == c->pos); p->is_tail = (pos == s->end); } return 1; } /*********************** *** Pileup iterator *** ***********************/ // Dictionary of overlapping reads KHASH_MAP_INIT_STR(olap_hash, lbnode_t *) typedef khash_t(olap_hash) olap_hash_t; struct __bam_plp_t { mempool_t *mp; lbnode_t *head, *tail, *dummy; int32_t tid, pos, max_tid, max_pos; int is_eof, max_plp, error, maxcnt; uint64_t id; bam_pileup1_t *plp; // for the "auto" interface only bam1_t *b; bam_plp_auto_f func; void *data; olap_hash_t *overlaps; }; bam_plp_t bam_plp_init(bam_plp_auto_f func, void *data) { bam_plp_t iter; iter = (bam_plp_t)calloc(1, sizeof(struct __bam_plp_t)); iter->mp = mp_init(); iter->head = iter->tail = mp_alloc(iter->mp); iter->dummy = mp_alloc(iter->mp); iter->max_tid = iter->max_pos = -1; iter->maxcnt = 8000; if (func) { iter->func = func; iter->data = data; iter->b = bam_init1(); } return iter; } void bam_plp_init_overlaps(bam_plp_t iter) { iter->overlaps = kh_init(olap_hash); // hash for tweaking quality of bases in overlapping reads } void bam_plp_destroy(bam_plp_t iter) { if ( iter->overlaps ) kh_destroy(olap_hash, iter->overlaps); mp_free(iter->mp, iter->dummy); mp_free(iter->mp, iter->head); if (iter->mp->cnt != 0) fprintf(stderr, "[bam_plp_destroy] memory leak: %d. Continue anyway.\n", iter->mp->cnt); mp_destroy(iter->mp); if (iter->b) bam_destroy1(iter->b); free(iter->plp); free(iter); } //--------------------------------- //--- Tweak overlapping reads //--------------------------------- /** * cigar_iref2iseq_set() - find the first CMATCH setting the ref and the read index * cigar_iref2iseq_next() - get the next CMATCH base * @cigar: pointer to current cigar block (rw) * @cigar_max: pointer just beyond the last cigar block * @icig: position within the current cigar block (rw) * @iseq: position in the sequence (rw) * @iref: position with respect to the beginning of the read (iref_pos - b->core.pos) (rw) * * Returns BAM_CMATCH or -1 when there is no more cigar to process or the requested position is not covered. */ static inline int cigar_iref2iseq_set(uint32_t **cigar, uint32_t *cigar_max, int *icig, int *iseq, int *iref) { int pos = *iref; if ( pos < 0 ) return -1; *icig = 0; *iseq = 0; *iref = 0; while ( *cigar<cigar_max ) { int cig = (**cigar) & BAM_CIGAR_MASK; int ncig = (**cigar) >> BAM_CIGAR_SHIFT; if ( cig==BAM_CSOFT_CLIP ) { (*cigar)++; *iseq += ncig; *icig = 0; continue; } if ( cig==BAM_CHARD_CLIP || cig==BAM_CPAD ) { (*cigar)++; *icig = 0; continue; } if ( cig==BAM_CMATCH || cig==BAM_CEQUAL || cig==BAM_CDIFF ) { pos -= ncig; if ( pos < 0 ) { *icig = ncig + pos; *iseq += *icig; *iref += *icig; return BAM_CMATCH; } (*cigar)++; *iseq += ncig; *icig = 0; *iref += ncig; continue; } if ( cig==BAM_CINS ) { (*cigar)++; *iseq += ncig; *icig = 0; continue; } if ( cig==BAM_CDEL || cig==BAM_CREF_SKIP ) { pos -= ncig; if ( pos<0 ) pos = 0; (*cigar)++; *icig = 0; *iref += ncig; continue; } fprintf(stderr,"todo: cigar %d\n", cig); assert(0); } *iseq = -1; return -1; } static inline int cigar_iref2iseq_next(uint32_t **cigar, uint32_t *cigar_max, int *icig, int *iseq, int *iref) { while ( *cigar < cigar_max ) { int cig = (**cigar) & BAM_CIGAR_MASK; int ncig = (**cigar) >> BAM_CIGAR_SHIFT; if ( cig==BAM_CMATCH || cig==BAM_CEQUAL || cig==BAM_CDIFF ) { if ( *icig >= ncig - 1 ) { *icig = 0; (*cigar)++; continue; } (*iseq)++; (*icig)++; (*iref)++; return BAM_CMATCH; } if ( cig==BAM_CDEL || cig==BAM_CREF_SKIP ) { (*cigar)++; (*iref) += ncig; *icig = 0; continue; } if ( cig==BAM_CINS ) { (*cigar)++; *iseq += ncig; *icig = 0; continue; } if ( cig==BAM_CSOFT_CLIP ) { (*cigar)++; *iseq += ncig; *icig = 0; continue; } if ( cig==BAM_CHARD_CLIP || cig==BAM_CPAD ) { (*cigar)++; *icig = 0; continue; } fprintf(stderr,"todo: cigar %d\n", cig); assert(0); } *iseq = -1; *iref = -1; return -1; } static void tweak_overlap_quality(bam1_t *a, bam1_t *b) { uint32_t *a_cigar = bam_get_cigar(a), *a_cigar_max = a_cigar + a->core.n_cigar; uint32_t *b_cigar = bam_get_cigar(b), *b_cigar_max = b_cigar + b->core.n_cigar; int a_icig = 0, a_iseq = 0; int b_icig = 0, b_iseq = 0; uint8_t *a_qual = bam_get_qual(a), *b_qual = bam_get_qual(b); uint8_t *a_seq = bam_get_seq(a), *b_seq = bam_get_seq(b); int iref = b->core.pos; int a_iref = iref - a->core.pos; int b_iref = iref - b->core.pos; int a_ret = cigar_iref2iseq_set(&a_cigar, a_cigar_max, &a_icig, &a_iseq, &a_iref); if ( a_ret<0 ) return; // no overlap int b_ret = cigar_iref2iseq_set(&b_cigar, b_cigar_max, &b_icig, &b_iseq, &b_iref); if ( b_ret<0 ) return; // no overlap #if DBG 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, 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))); #endif while ( 1 ) { // Increment reference position while ( a_iref>=0 && a_iref < iref - a->core.pos ) a_ret = cigar_iref2iseq_next(&a_cigar, a_cigar_max, &a_icig, &a_iseq, &a_iref); if ( a_ret<0 ) break; // done if ( iref < a_iref + a->core.pos ) iref = a_iref + a->core.pos; while ( b_iref>=0 && b_iref < iref - b->core.pos ) b_ret = cigar_iref2iseq_next(&b_cigar, b_cigar_max, &b_icig, &b_iseq, &b_iref); if ( b_ret<0 ) break; // done if ( iref < b_iref + b->core.pos ) iref = b_iref + b->core.pos; iref++; if ( a_iref+a->core.pos != b_iref+b->core.pos ) continue; // only CMATCH positions, don't know what to do with indels if ( bam_seqi(a_seq,a_iseq) == bam_seqi(b_seq,b_iseq) ) { #if DBG fprintf(stderr,"%c",seq_nt16_str[bam_seqi(a_seq,a_iseq)]); #endif // we are very confident about this base int qual = a_qual[a_iseq] + b_qual[b_iseq]; a_qual[a_iseq] = qual>200 ? 200 : qual; b_qual[b_iseq] = 0; } else { if ( a_qual[a_iseq] >= b_qual[b_iseq] ) { #if DBG fprintf(stderr,"[%c/%c]",seq_nt16_str[bam_seqi(a_seq,a_iseq)],tolower(seq_nt16_str[bam_seqi(b_seq,b_iseq)])); #endif a_qual[a_iseq] = 0.8 * a_qual[a_iseq]; // not so confident about a_qual anymore given the mismatch b_qual[b_iseq] = 0; } else { #if DBG fprintf(stderr,"[%c/%c]",tolower(seq_nt16_str[bam_seqi(a_seq,a_iseq)]),seq_nt16_str[bam_seqi(b_seq,b_iseq)]); #endif b_qual[b_iseq] = 0.8 * b_qual[b_iseq]; a_qual[a_iseq] = 0; } } } #if DBG fprintf(stderr,"\n"); #endif } // Fix overlapping reads. Simple soft-clipping did not give good results. // Lowering qualities of unwanted bases is more selective and works better. // static void overlap_push(bam_plp_t iter, lbnode_t *node) { if ( !iter->overlaps ) return; // mapped mates and paired reads only if ( node->b.core.flag&BAM_FMUNMAP || !(node->b.core.flag&BAM_FPROPER_PAIR) ) return; // no overlap possible, unless some wild cigar if ( abs(node->b.core.isize) >= 2*node->b.core.l_qseq ) return; khiter_t kitr = kh_get(olap_hash, iter->overlaps, bam_get_qname(&node->b)); if ( kitr==kh_end(iter->overlaps) ) { int ret; kitr = kh_put(olap_hash, iter->overlaps, bam_get_qname(&node->b), &ret); kh_value(iter->overlaps, kitr) = node; } else { lbnode_t *a = kh_value(iter->overlaps, kitr); tweak_overlap_quality(&a->b, &node->b); kh_del(olap_hash, iter->overlaps, kitr); assert(a->end-1 == a->s.end); a->end = a->b.core.pos + bam_cigar2rlen(a->b.core.n_cigar, bam_get_cigar(&a->b)); a->s.end = a->end - 1; } } static void overlap_remove(bam_plp_t iter, const bam1_t *b) { if ( !iter->overlaps ) return; khiter_t kitr; if ( b ) { kitr = kh_get(olap_hash, iter->overlaps, bam_get_qname(b)); if ( kitr!=kh_end(iter->overlaps) ) kh_del(olap_hash, iter->overlaps, kitr); } else { // remove all for (kitr = kh_begin(iter->overlaps); kitr<kh_end(iter->overlaps); kitr++) if ( kh_exist(iter->overlaps, kitr) ) kh_del(olap_hash, iter->overlaps, kitr); } } // Prepares next pileup position in bam records collected by bam_plp_auto -> user func -> bam_plp_push. Returns // pointer to the piled records if next position is ready or NULL if there is not enough records in the // buffer yet (the current position is still the maximum position across all buffered reads). const bam_pileup1_t *bam_plp_next(bam_plp_t iter, int *_tid, int *_pos, int *_n_plp) { if (iter->error) { *_n_plp = -1; return 0; } *_n_plp = 0; if (iter->is_eof && iter->head->next == 0) return 0; while (iter->is_eof || iter->max_tid > iter->tid || (iter->max_tid == iter->tid && iter->max_pos > iter->pos)) { int n_plp = 0; lbnode_t *p, *q; // write iter->plp at iter->pos iter->dummy->next = iter->head; for (p = iter->head, q = iter->dummy; p->next; q = p, p = p->next) { if (p->b.core.tid < iter->tid || (p->b.core.tid == iter->tid && p->end <= iter->pos)) { // then remove overlap_remove(iter, &p->b); q->next = p->next; mp_free(iter->mp, p); p = q; } else if (p->b.core.tid == iter->tid && p->beg <= iter->pos) { // here: p->end > pos; then add to pileup if (n_plp == iter->max_plp) { // then double the capacity iter->max_plp = iter->max_plp? iter->max_plp<<1 : 256; iter->plp = (bam_pileup1_t*)realloc(iter->plp, sizeof(bam_pileup1_t) * iter->max_plp); } iter->plp[n_plp].b = &p->b; if (resolve_cigar2(iter->plp + n_plp, iter->pos, &p->s)) ++n_plp; // actually always true... } } iter->head = iter->dummy->next; // dummy->next may be changed *_n_plp = n_plp; *_tid = iter->tid; *_pos = iter->pos; // update iter->tid and iter->pos if (iter->head->next) { if (iter->tid > iter->head->b.core.tid) { fprintf(stderr, "[%s] unsorted input. Pileup aborts.\n", __func__); iter->error = 1; *_n_plp = -1; return 0; } } if (iter->tid < iter->head->b.core.tid) { // come to a new reference sequence iter->tid = iter->head->b.core.tid; iter->pos = iter->head->beg; // jump to the next reference } else if (iter->pos < iter->head->beg) { // here: tid == head->b.core.tid iter->pos = iter->head->beg; // jump to the next position } else ++iter->pos; // scan contiguously // return if (n_plp) return iter->plp; if (iter->is_eof && iter->head->next == 0) break; } return 0; } int bam_plp_push(bam_plp_t iter, const bam1_t *b) { if (iter->error) return -1; if (b) { if (b->core.tid < 0) { overlap_remove(iter, b); return 0; } // Skip only unmapped reads here, any additional filtering must be done in iter->func if (b->core.flag & BAM_FUNMAP) { overlap_remove(iter, b); return 0; } if (iter->tid == b->core.tid && iter->pos == b->core.pos && iter->mp->cnt > iter->maxcnt) { overlap_remove(iter, b); return 0; } bam_copy1(&iter->tail->b, b); overlap_push(iter, iter->tail); #ifndef BAM_NO_ID iter->tail->b.id = iter->id++; #endif iter->tail->beg = b->core.pos; iter->tail->end = b->core.pos + bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b)); iter->tail->s = g_cstate_null; iter->tail->s.end = iter->tail->end - 1; // initialize cstate_t if (b->core.tid < iter->max_tid) { fprintf(stderr, "[bam_pileup_core] the input is not sorted (chromosomes out of order)\n"); iter->error = 1; return -1; } if ((b->core.tid == iter->max_tid) && (iter->tail->beg < iter->max_pos)) { fprintf(stderr, "[bam_pileup_core] the input is not sorted (reads out of order)\n"); iter->error = 1; return -1; } iter->max_tid = b->core.tid; iter->max_pos = iter->tail->beg; if (iter->tail->end > iter->pos || iter->tail->b.core.tid > iter->tid) { iter->tail->next = mp_alloc(iter->mp); iter->tail = iter->tail->next; } } else iter->is_eof = 1; return 0; } const bam_pileup1_t *bam_plp_auto(bam_plp_t iter, int *_tid, int *_pos, int *_n_plp) { const bam_pileup1_t *plp; if (iter->func == 0 || iter->error) { *_n_plp = -1; return 0; } if ((plp = bam_plp_next(iter, _tid, _pos, _n_plp)) != 0) return plp; else { // no pileup line can be obtained; read alignments *_n_plp = 0; if (iter->is_eof) return 0; int ret; while ( (ret=iter->func(iter->data, iter->b)) >= 0) { if (bam_plp_push(iter, iter->b) < 0) { *_n_plp = -1; return 0; } if ((plp = bam_plp_next(iter, _tid, _pos, _n_plp)) != 0) return plp; // otherwise no pileup line can be returned; read the next alignment. } if ( ret < -1 ) { iter->error = ret; *_n_plp = -1; return 0; } bam_plp_push(iter, 0); if ((plp = bam_plp_next(iter, _tid, _pos, _n_plp)) != 0) return plp; return 0; } } void bam_plp_reset(bam_plp_t iter) { lbnode_t *p, *q; iter->max_tid = iter->max_pos = -1; iter->tid = iter->pos = 0; iter->is_eof = 0; for (p = iter->head; p->next;) { overlap_remove(iter, NULL); q = p->next; mp_free(iter->mp, p); p = q; } iter->head = iter->tail; } void bam_plp_set_maxcnt(bam_plp_t iter, int maxcnt) { iter->maxcnt = maxcnt; } /************************ *** Mpileup iterator *** ************************/ struct __bam_mplp_t { int n; uint64_t min, *pos; bam_plp_t *iter; int *n_plp; const bam_pileup1_t **plp; }; bam_mplp_t bam_mplp_init(int n, bam_plp_auto_f func, void **data) { int i; bam_mplp_t iter; iter = (bam_mplp_t)calloc(1, sizeof(struct __bam_mplp_t)); iter->pos = (uint64_t*)calloc(n, sizeof(uint64_t)); iter->n_plp = (int*)calloc(n, sizeof(int)); iter->plp = (const bam_pileup1_t**)calloc(n, sizeof(bam_pileup1_t*)); iter->iter = (bam_plp_t*)calloc(n, sizeof(bam_plp_t)); iter->n = n; iter->min = (uint64_t)-1; for (i = 0; i < n; ++i) { iter->iter[i] = bam_plp_init(func, data[i]); iter->pos[i] = iter->min; } return iter; } void bam_mplp_init_overlaps(bam_mplp_t iter) { int i; for (i = 0; i < iter->n; ++i) bam_plp_init_overlaps(iter->iter[i]); } void bam_mplp_set_maxcnt(bam_mplp_t iter, int maxcnt) { int i; for (i = 0; i < iter->n; ++i) iter->iter[i]->maxcnt = maxcnt; } void bam_mplp_destroy(bam_mplp_t iter) { int i; for (i = 0; i < iter->n; ++i) bam_plp_destroy(iter->iter[i]); free(iter->iter); free(iter->pos); free(iter->n_plp); free(iter->plp); free(iter); } int bam_mplp_auto(bam_mplp_t iter, int *_tid, int *_pos, int *n_plp, const bam_pileup1_t **plp) { int i, ret = 0; uint64_t new_min = (uint64_t)-1; for (i = 0; i < iter->n; ++i) { if (iter->pos[i] == iter->min) { int tid, pos; iter->plp[i] = bam_plp_auto(iter->iter[i], &tid, &pos, &iter->n_plp[i]); if ( iter->iter[i]->error ) return -1; iter->pos[i] = iter->plp[i] ? (uint64_t)tid<<32 | pos : 0; } if (iter->plp[i] && iter->pos[i] < new_min) new_min = iter->pos[i]; } iter->min = new_min; if (new_min == (uint64_t)-1) return 0; *_tid = new_min>>32; *_pos = (uint32_t)new_min; for (i = 0; i < iter->n; ++i) { if (iter->pos[i] == iter->min) { // FIXME: valgrind reports "uninitialised value(s) at this line" n_plp[i] = iter->n_plp[i], plp[i] = iter->plp[i]; ++ret; } else n_plp[i] = 0, plp[i] = 0; } return ret; } #endif // ~!defined(BAM_NO_PILEUP)