From 253dbf083d8a84721ae6ac63af786c8dcf566ecb Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 27 Oct 2010 04:14:43 +0000 Subject: [PATCH] * samtools-0.1.8-20 (r778) * speed up pileup, although I do not know how much is the improvement --- bam.h | 2 +- bam_pileup.c | 139 +++++++++++++++++++++++++++++---------------------- bam_plcmd.c | 41 +++++++++++---- bam_tview.c | 2 +- bamtk.c | 2 +- 5 files changed, 114 insertions(+), 72 deletions(-) diff --git a/bam.h b/bam.h index a8769e2..5943123 100644 --- a/bam.h +++ b/bam.h @@ -495,7 +495,7 @@ extern "C" { bam1_t *b; int32_t qpos; int indel, level; - uint32_t is_del:1, is_head:1, is_tail:1; + uint32_t is_del:1, is_head:1, is_tail:1, is_refskip:1; } bam_pileup1_t; typedef int (*bam_plp_auto_f)(void *data, bam1_t *b); diff --git a/bam_pileup.c b/bam_pileup.c index 266b966..62d5fa1 100644 --- a/bam_pileup.c +++ b/bam_pileup.c @@ -4,9 +4,16 @@ #include #include "sam.h" +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; uint32_t beg, end; + cstate_t s; struct __linkbuf_t *next; } lbnode_t; @@ -53,72 +60,83 @@ static inline void mp_free(mempool_t *mp, lbnode_t *p) /* --- BEGIN: Auxiliary functions */ -static inline int resolve_cigar(bam_pileup1_t *p, uint32_t pos) +/* 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, uint32_t pos, cstate_t *s) { - unsigned k; +#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 x = c->pos, y = 0; - int ret = 1, is_restart = 1, first_op, last_op; - - if (c->flag&BAM_FUNMAP) return 0; // unmapped read - assert(x <= pos); // otherwise a bug - p->qpos = -1; p->indel = 0; p->is_del = p->is_head = p->is_tail = 0; - first_op = bam1_cigar(b)[0] & BAM_CIGAR_MASK; - last_op = bam1_cigar(b)[c->n_cigar-1] & BAM_CIGAR_MASK; - for (k = 0; k < c->n_cigar; ++k) { - int op = bam1_cigar(b)[k] & BAM_CIGAR_MASK; // operation - int l = bam1_cigar(b)[k] >> BAM_CIGAR_SHIFT; // length - if (op == BAM_CMATCH) { // NOTE: this assumes the first and the last operation MUST BE a match or a clip - if (x + l > pos) { // overlap with pos - p->indel = p->is_del = 0; - p->qpos = y + (pos - x); - if (x == pos && is_restart && first_op != BAM_CDEL) p->is_head = 1; - if (x + l - 1 == pos) { // come to the end of a match - int has_next_match = 0; - unsigned i; - for (i = k + 1; i < c->n_cigar; ++i) { - uint32_t cigar = bam1_cigar(b)[i]; - int opi = cigar&BAM_CIGAR_MASK; - if (opi == BAM_CMATCH) { - has_next_match = 1; - break; - } else if (opi == BAM_CSOFT_CLIP || opi == BAM_CREF_SKIP || opi == BAM_CHARD_CLIP) break; - } - if (!has_next_match && last_op != BAM_CDEL) p->is_tail = 1; - if (k < c->n_cigar - 1 && has_next_match) { // there are additional operation(s) - uint32_t cigar = bam1_cigar(b)[k+1]; // next CIGAR - int op_next = cigar&BAM_CIGAR_MASK; // next CIGAR operation - if (op_next == BAM_CDEL) p->indel = -(int32_t)(cigar>>BAM_CIGAR_SHIFT); // del - else if (op_next == BAM_CINS) p->indel = cigar>>BAM_CIGAR_SHIFT; // ins - else if (op_next == BAM_CPAD && k + 2 < c->n_cigar) { // no working for adjacent padding - cigar = bam1_cigar(b)[k+2]; op_next = cigar&BAM_CIGAR_MASK; - if (op_next == BAM_CDEL) p->indel = -(int32_t)(cigar>>BAM_CIGAR_SHIFT); // del - else if (op_next == BAM_CINS) p->indel = cigar>>BAM_CIGAR_SHIFT; // ins - } - } + uint32_t *cigar = bam1_cigar(b); + int k, is_head = 0; + // determine the current CIGAR operation +// fprintf(stderr, "%s\tpos=%d\tend=%d\t(%d,%d,%d)\n", bam1_qname(b), pos, s->end, s->k, s->x, s->y); + if (s->k == -1) { // never processed + is_head = 1; + if (c->n_cigar == 1) { // just one operation, save a loop + if (_cop(cigar[0]) == BAM_CMATCH) s->k = 0, s->x = c->pos, s->y = 0; + } else { // find the first match + 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) break; + else if (op == BAM_CDEL || 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) { // jump to the next without a loop + if (_cop(cigar[s->k]) == BAM_CMATCH) s->y += l; + s->x += l; + ++s->k; + } else { // find the next M/D/N + if (_cop(cigar[s->k]) == BAM_CMATCH) 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) break; + else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) s->y += l; } + s->k = k; } - x += l; y += l; - } else if (op == BAM_CDEL) { // then set ->is_del - if (k == 0 && x == pos) p->is_head = 1; - else if (x + l - 1 == pos && k == c->n_cigar - 1) p->is_tail = 1; - if (x + l > pos) { - p->indel = 0; p->is_del = 1; - p->qpos = y; + 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 (op == BAM_CMATCH) { + p->qpos = s->y + (pos - s->x); + 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 + op2 = _cop(cigar[s->k+2]); l2 = _cln(cigar[s->k+2]); + if (op2 == BAM_CDEL) p->indel = -(int)l2; + else if (op2 == BAM_CINS) p->indel = l2; + } } - x += l; - } else if (op == BAM_CREF_SKIP) x += l; - else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) y += l; - if (is_restart) is_restart ^= (op == BAM_CMATCH); - else is_restart ^= (op == BAM_CREF_SKIP || op == BAM_CSOFT_CLIP || op == BAM_CHARD_CLIP); - if (x > pos) { - if (op == BAM_CREF_SKIP) ret = 0; // then do not put it into pileup at all - break; - } + } 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); } - assert(x > pos); // otherwise a bug - return ret; + return 1; } /* --- END: Auxiliary functions */ @@ -187,7 +205,7 @@ const bam_pileup1_t *bam_plp_next(bam_plp_t iter, int *_tid, int *_pos, int *_n_ 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 + 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 @@ -221,6 +239,7 @@ int bam_plp_push(bam_plp_t iter, const bam1_t *b) 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)); + 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; diff --git a/bam_plcmd.c b/bam_plcmd.c index 4bdc351..cf98f69 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -91,6 +91,21 @@ static khash_t(64) *load_pos(const char *fn, bam_header_t *h) return hash; } +static inline int printw(int c, FILE *fp) +{ + char buf[16]; + int l, x; + if (c == 0) return fputc('0', fp); + for (l = 0, x = c < 0? -c : c; x > 0; x /= 10) buf[l++] = x%10 + '0'; + if (c < 0) buf[l++] = '-'; + buf[l] = 0; + for (x = 0; x < l/2; ++x) { + int y = buf[x]; buf[x] = buf[l-1-x]; buf[l-1-x] = y; + } + fputs(buf, fp); + return 0; +} + // an analogy to pileup_func() below static int glt3_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pu, void *data) { @@ -162,7 +177,10 @@ static int glt3_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pu, static void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, const char *ref) { - if (p->is_head) printf("^%c", p->b->core.qual > 93? 126 : p->b->core.qual + 33); + if (p->is_head) { + putchar('^'); + putchar(p->b->core.qual > 93? 126 : p->b->core.qual + 33); + } if (!p->is_del) { int j, rb, c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos)]; rb = (ref && pos < ref_len)? ref[pos] : 'N'; @@ -170,19 +188,19 @@ static void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, const char else c = bam1_strand(p->b)? tolower(c) : toupper(c); putchar(c); if (p->indel > 0) { - printf("+%d", p->indel); + putchar('+'); printw(p->indel, stdout); for (j = 1; j <= p->indel; ++j) { c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos + j)]; putchar(bam1_strand(p->b)? tolower(c) : toupper(c)); } } else if (p->indel < 0) { - printf("%d", p->indel); + printw(p->indel, stdout); for (j = 1; j <= -p->indel; ++j) { c = (ref && (int)pos+j < ref_len)? ref[pos+j] : 'N'; putchar(bam1_strand(p->b)? tolower(c) : toupper(c)); } } - } else putchar('*'); + } else putchar(p->is_refskip? (bam1_strand(p->b)? '<' : '>') : '*'); if (p->is_tail) putchar('$'); } @@ -256,12 +274,17 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p } } // print the first 3 columns - printf("%s\t%d\t%c\t", d->h->target_name[tid], pos + 1, rb); + fputs(d->h->target_name[tid], stdout); putchar('\t'); + printw(pos+1, stdout); putchar('\t'); putchar(rb); putchar('\t'); // print consensus information if required - if (d->format & BAM_PLF_CNS) - printf("%c\t%d\t%d\t%d\t", bam_nt16_rev_table[cns>>28], cns>>8&0xff, cns&0xff, cns>>16&0xff); + if (d->format & BAM_PLF_CNS) { + putchar(bam_nt16_rev_table[cns>>28]); putchar('\t'); + printw(cns>>8&0xff, stdout); putchar('\t'); + printw(cns&0xff, stdout); putchar('\t'); + printw(cns>>16&0xff, stdout); putchar('\t'); + } // print pileup sequences - printf("%d\t", n); + printw(n, stdout); putchar('\t'); rms_aux = 0; // we need to recalculate rms_mapq when -c is not flagged on the command line for (i = 0; i < n; ++i) { const bam_pileup1_t *p = pu + i; @@ -310,7 +333,7 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p for (i = 0; i < n; ++i) { int x = pu[i].qpos; int l = pu[i].b->core.l_qseq; - printf("%d,", x < l/2? x+1 : -((l-1)-x+1)); + printw(x < l/2? x+1 : -((l-1)-x+1), stdout); putchar(','); } } putchar('\n'); diff --git a/bam_tview.c b/bam_tview.c index 24cd12a..e48afa7 100644 --- a/bam_tview.c +++ b/bam_tview.c @@ -109,7 +109,7 @@ int tv_pl_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void if (tv->is_dot && toupper(c) == toupper(rb)) c = bam1_strand(p->b)? ',' : '.'; } } - } else c = '*'; + } else c = p->is_refskip? (bam1_strand(p->b)? '<' : '>') : '*'; } else { // padding if (j > p->indel) c = '*'; else { // insertion diff --git a/bamtk.c b/bamtk.c index 4de99d7..fb0e388 100644 --- a/bamtk.c +++ b/bamtk.c @@ -9,7 +9,7 @@ #endif #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.1.8-19 (r777)" +#define PACKAGE_VERSION "0.1.8-20 (r778)" #endif int bam_taf2baf(int argc, char *argv[]); -- 2.39.2