* samtools-0.1.9-2 (r787)

[samtools.git] / bam_pileup.c

diff --git a/bam_pileup.c b/bam_pileup.c
index 62d5fa1d2b3d9b2171b98eccb9f583755e379822..55e51e27804cfc58fc5576b40ae7ba10416a3e8f 100644 (file)
--- a/bam_pileup.c
+++ b/bam_pileup.c
@@ -114,22 +114,25 @@ static inline int resolve_cigar2(bam_pileup1_t *p, uint32_t pos, cstate_t *s)
        }
        { // collect pileup information
                int op, l;
-               op = _cop(cigar[s->k]);
-               l = _cln(cigar[s->k]);
+               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;
+               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) break;
                                }
+                               if (l3 > 0) p->indel = l3;
                        }
+               }
+               if (op == BAM_CMATCH) {
+                       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);
@@ -149,7 +152,7 @@ struct __bam_plp_t {
        mempool_t *mp;
        lbnode_t *head, *tail, *dummy;
        int32_t tid, pos, max_tid, max_pos;
-       int is_eof, flag_mask, max_plp, error;
+       int is_eof, flag_mask, max_plp, error, maxcnt;
        bam_pileup1_t *plp;
        // for the "auto" interface only
        bam1_t *b;
@@ -166,6 +169,7 @@ bam_plp_t bam_plp_init(bam_plp_auto_f func, void *data)
        iter->dummy = mp_alloc(iter->mp);
        iter->max_tid = iter->max_pos = -1;
        iter->flag_mask = BAM_DEF_MASK;
+       iter->maxcnt = 8000;
        if (func) {
                iter->func = func;
                iter->data = data;
@@ -237,6 +241,7 @@ int bam_plp_push(bam_plp_t iter, const bam1_t *b)
        if (b) {
                if (b->core.tid < 0) return 0;
                if (b->core.flag & iter->flag_mask) return 0;
+               if (iter->tid == b->core.tid && iter->pos == b->core.pos && iter->mp->cnt > iter->maxcnt) 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
@@ -264,7 +269,7 @@ const bam_pileup1_t *bam_plp_auto(bam_plp_t iter, int *_tid, int *_pos, int *_n_
        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 {
+       else { // no pileup line can be obtained; read alignments
                *_n_plp = 0;
                if (iter->is_eof) return 0;
                while (iter->func(iter->data, iter->b) >= 0) {
@@ -273,6 +278,7 @@ const bam_pileup1_t *bam_plp_auto(bam_plp_t iter, int *_tid, int *_pos, int *_n_
                                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.
                }
                bam_plp_push(iter, 0);
                if ((plp = bam_plp_next(iter, _tid, _pos, _n_plp)) != 0) return plp;
@@ -299,6 +305,11 @@ 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_set_maxcnt(bam_plp_t iter, int maxcnt)
+{
+       iter->maxcnt = maxcnt;
+}
+
 /*****************
  * callback APIs *
  *****************/
@@ -386,6 +397,13 @@ bam_mplp_t bam_mplp_init(int n, bam_plp_auto_f func, void **data)
        return iter;
 }
 
+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;