From: Heng Li Date: Sat, 12 Jun 2010 02:37:09 +0000 (+0000) Subject: * samtools-0.1.7-12 (r589) X-Git-Url: https://git.donarmstrong.com/?p=samtools.git;a=commitdiff_plain;h=bc8511255e4a4ad7eb3233c4e9b65dfbe70a87ed * samtools-0.1.7-12 (r589) * added iterator-like APIs for pileup --- diff --git a/bam.h b/bam.h index 86939b2..21a7598 100644 --- a/bam.h +++ b/bam.h @@ -463,11 +463,15 @@ extern "C" { uint32_t is_del:1, is_head:1, is_tail:1; } bam_pileup1_t; - struct __bam_plbuf_t; - /*! @abstract pileup buffer */ - typedef struct __bam_plbuf_t bam_plbuf_t; + struct __bam_plp_t; + typedef struct __bam_plp_t *bam_plp_t; - void bam_plbuf_set_mask(bam_plbuf_t *buf, int mask); + bam_plp_t bam_plp_init(); + int bam_plp_push(bam_plp_t iter, const bam1_t *b); + const bam_pileup1_t *bam_plp_next(bam_plp_t iter, int *_n_plp, int *_tid, int *_pos); + void bam_plp_set_mask(bam_plp_t iter, int mask); + void bam_plp_reset(bam_plp_t iter); + void bam_plp_destroy(bam_plp_t iter); /*! @typedef @abstract Type of function to be called by bam_plbuf_push(). @@ -480,44 +484,16 @@ extern "C" { */ typedef int (*bam_pileup_f)(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data); - /*! - @abstract Reset a pileup buffer for another pileup process - @param buf the pileup buffer to be reset - */ - void bam_plbuf_reset(bam_plbuf_t *buf); + typedef struct { + bam_plp_t iter; + bam_pileup_f func; + void *data; + } bam_plbuf_t; - /*! - @abstract Initialize a buffer for pileup. - @param func fucntion to be called by bam_pileup_core() - @param data user provided data - @return pointer to the pileup buffer - */ + void bam_plbuf_set_mask(bam_plbuf_t *buf, int mask); + void bam_plbuf_reset(bam_plbuf_t *buf); bam_plbuf_t *bam_plbuf_init(bam_pileup_f func, void *data); - - /*! - @abstract Destroy a pileup buffer. - @param buf pointer to the pileup buffer - */ void bam_plbuf_destroy(bam_plbuf_t *buf); - - /*! - @abstract Push an alignment to the pileup buffer. - @param b alignment to be pushed - @param buf pileup buffer - @see bam_plbuf_init() - @return always 0 currently - - @discussion If all the alignments covering a particular site have - been collected, this function will call the user defined function - as is provided to bam_plbuf_init(). The coordinate of the site and - all the alignments will be transferred to the user defined - function as function parameters. - - When all the alignments are pushed to the buffer, this function - needs to be called with b equal to NULL. This will flush the - buffer. A pileup buffer can only be reused when bam_plbuf_reset() - is called. - */ int bam_plbuf_push(const bam1_t *b, bam_plbuf_t *buf); int bam_pileup_file(bamFile fp, int mask, bam_pileup_f func, void *func_data); diff --git a/bam_pileup.c b/bam_pileup.c index f68f400..5a87746 100644 --- a/bam_pileup.c +++ b/bam_pileup.c @@ -108,119 +108,131 @@ static inline int resolve_cigar(bam_pileup1_t *p, uint32_t pos) /* --- END: Auxiliary functions */ -struct __bam_plbuf_t { +/******************* + * pileup iterator * + *******************/ + +struct __bam_plp_t { mempool_t *mp; lbnode_t *head, *tail, *dummy; - bam_pileup_f func; - void *func_data; int32_t tid, pos, max_tid, max_pos; - int max_pu, is_eof; - bam_pileup1_t *pu; - int flag_mask; + int is_eof, flag_mask, max_plp; + bam_pileup1_t *plp; }; -void bam_plbuf_reset(bam_plbuf_t *buf) +bam_plp_t bam_plp_init(void) { - lbnode_t *p, *q; - buf->max_tid = buf->max_pos = -1; - buf->tid = buf->pos = 0; - buf->is_eof = 0; - for (p = buf->head; p->next;) { - q = p->next; - mp_free(buf->mp, p); - p = q; - } - buf->head = buf->tail; + bam_plp_t iter; + iter = 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->flag_mask = BAM_DEF_MASK; + return iter; } -void bam_plbuf_set_mask(bam_plbuf_t *buf, int mask) +const bam_pileup1_t *bam_plp_next(bam_plp_t iter, int *_n_plp, int *_tid, int *_pos) { - if (mask < 0) buf->flag_mask = BAM_DEF_MASK; - else buf->flag_mask = BAM_FUNMAP | mask; -} - -bam_plbuf_t *bam_plbuf_init(bam_pileup_f func, void *data) -{ - bam_plbuf_t *buf; - buf = (bam_plbuf_t*)calloc(1, sizeof(bam_plbuf_t)); - buf->func = func; buf->func_data = data; - buf->mp = mp_init(); - buf->head = buf->tail = mp_alloc(buf->mp); - buf->dummy = mp_alloc(buf->mp); - buf->max_tid = buf->max_pos = -1; - buf->flag_mask = BAM_DEF_MASK; - return buf; -} - -void bam_plbuf_destroy(bam_plbuf_t *buf) -{ - mp_free(buf->mp, buf->dummy); - mp_free(buf->mp, buf->head); - if (buf->mp->cnt != 0) - fprintf(stderr, "[bam_plbuf_destroy] memory leak: %d. Continue anyway.\n", buf->mp->cnt); - mp_destroy(buf->mp); - free(buf->pu); - free(buf); + *_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 + 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_cigar(iter->plp + n_plp, iter->pos)) ++n_plp; // skip the read if we are looking at ref-skip + } + } + 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__); + *_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_plbuf_push(const bam1_t *b, bam_plbuf_t *buf) +int bam_plp_push(bam_plp_t iter, const bam1_t *b) { - if (b) { // fill buffer + if (b) { if (b->core.tid < 0) return 0; - if (b->core.flag & buf->flag_mask) return 0; - bam_copy1(&buf->tail->b, b); - buf->tail->beg = b->core.pos; buf->tail->end = bam_calend(&b->core, bam1_cigar(b)); - if (b->core.tid < buf->max_tid) { + if (b->core.flag & iter->flag_mask) return 0; + bam_copy1(&iter->tail->b, b); + iter->tail->beg = b->core.pos; iter->tail->end = bam_calend(&b->core, bam1_cigar(b)); + if (b->core.tid < iter->max_tid) { fprintf(stderr, "[bam_pileup_core] the input is not sorted (chromosomes out of order)\n"); return -1; } - if ((b->core.tid == buf->max_tid) && (buf->tail->beg < buf->max_pos)) { + 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"); return -1; } - buf->max_tid = b->core.tid; buf->max_pos = buf->tail->beg; - if (buf->tail->end > buf->pos || buf->tail->b.core.tid > buf->tid) { - buf->tail->next = mp_alloc(buf->mp); - buf->tail = buf->tail->next; - } - } else buf->is_eof = 1; - while (buf->is_eof || buf->max_tid > buf->tid || (buf->max_tid == buf->tid && buf->max_pos > buf->pos)) { - int n_pu = 0; - lbnode_t *p, *q; - buf->dummy->next = buf->head; - for (p = buf->head, q = buf->dummy; p->next; q = p, p = p->next) { - if (p->b.core.tid < buf->tid || (p->b.core.tid == buf->tid && p->end <= buf->pos)) { // then remove from the list - q->next = p->next; mp_free(buf->mp, p); p = q; - } else if (p->b.core.tid == buf->tid && p->beg <= buf->pos) { // here: p->end > pos; then add to pileup - if (n_pu == buf->max_pu) { // then double the capacity - buf->max_pu = buf->max_pu? buf->max_pu<<1 : 256; - buf->pu = (bam_pileup1_t*)realloc(buf->pu, sizeof(bam_pileup1_t) * buf->max_pu); - } - buf->pu[n_pu].b = &p->b; - if (resolve_cigar(buf->pu + n_pu, buf->pos)) ++n_pu; // skip the read if we are looking at BAM_CREF_SKIP - } + 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; } - buf->head = buf->dummy->next; // dummy->next may be changed - if (n_pu) { // then call user defined function - buf->func(buf->tid, buf->pos, n_pu, buf->pu, buf->func_data); - } - // update tid and pos - if (buf->head->next) { - if (buf->tid > buf->head->b.core.tid) { - fprintf(stderr, "[bam_plbuf_push] unsorted input. Pileup aborts.\n"); - return 1; - } - } - if (buf->tid < buf->head->b.core.tid) { // come to a new reference sequence - buf->tid = buf->head->b.core.tid; buf->pos = buf->head->beg; // jump to the next reference - } else if (buf->pos < buf->head->beg) { // here: tid == head->b.core.tid - buf->pos = buf->head->beg; // jump to the next position - } else ++buf->pos; // scan contiguously - if (buf->is_eof && buf->head->next == 0) break; - } + } else iter->is_eof = 1; 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;) { + q = p->next; + mp_free(iter->mp, p); + p = q; + } + iter->head = iter->tail; +} + +void bam_plp_set_mask(bam_plp_t iter, int mask) +{ + iter->flag_mask = mask < 0? BAM_DEF_MASK : (BAM_FUNMAP | mask); +} + +void bam_plp_destroy(bam_plp_t iter) +{ + 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); + free(iter->plp); + free(iter); +} + +/***************** + * callback APIs * + *****************/ + int bam_pileup_file(bamFile fp, int mask, bam_pileup_f func, void *func_data) { bam_plbuf_t *buf; @@ -236,3 +248,40 @@ int bam_pileup_file(bamFile fp, int mask, bam_pileup_f func, void *func_data) bam_destroy1(b); return 0; } + +void bam_plbuf_set_mask(bam_plbuf_t *buf, int mask) +{ + bam_plp_set_mask(buf->iter, mask); +} + +void bam_plbuf_reset(bam_plbuf_t *buf) +{ + bam_plp_reset(buf->iter); +} + +bam_plbuf_t *bam_plbuf_init(bam_pileup_f func, void *data) +{ + bam_plbuf_t *buf; + buf = calloc(1, sizeof(bam_plbuf_t)); + buf->iter = bam_plp_init(); + buf->func = func; + buf->data = data; + return buf; +} + +void bam_plbuf_destroy(bam_plbuf_t *buf) +{ + bam_plp_destroy(buf->iter); + free(buf); +} + +int bam_plbuf_push(const bam1_t *b, bam_plbuf_t *buf) +{ + int ret, n_plp, tid, pos; + const bam_pileup1_t *plp; + ret = bam_plp_push(buf->iter, b); + if (ret < 0) return ret; + while ((plp = bam_plp_next(buf->iter, &n_plp, &tid, &pos)) != 0) + buf->func(tid, pos, n_plp, plp, buf->data); + return 0; +} diff --git a/bamtk.c b/bamtk.c index 32357b3..417b64c 100644 --- a/bamtk.c +++ b/bamtk.c @@ -9,7 +9,7 @@ #endif #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.1.7-11 (r588)" +#define PACKAGE_VERSION "0.1.7-12 (r589)" #endif int bam_taf2baf(int argc, char *argv[]);