improve pileup a little bit

diff --git a/bam_pileup.c b/bam_pileup.c
index 62d5fa1d2b3d9b2171b98eccb9f583755e379822..a49f795d73d4f0ed1ad3a9592108f3be2a917579 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);

diff --git a/bam_plcmd.c b/bam_plcmd.c
index cf98f69565ee6d52176df5ca5c9d14a44d2ef3e6..f0e63ce62cd91cc6cada955eaa19c9b2493e1872 100644 (file)
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -175,32 +175,33 @@ static int glt3_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pu,
        return 0;
 }
 
-static void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, const char *ref)
+static inline void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, const char *ref)
 {
+       int j;
        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)];
+               int rb, c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos)];
                rb = (ref && pos < ref_len)? ref[pos] : 'N';
                if (c == '=' || toupper(c) == toupper(rb)) c = bam1_strand(p->b)? ',' : '.';
                else c = bam1_strand(p->b)? tolower(c) : toupper(c);
                putchar(c);
-               if (p->indel > 0) {
-                       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) {
-                       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(p->is_refskip? (bam1_strand(p->b)? '<' : '>') : '*');
+       if (p->indel > 0) {
+               putchar('+'); printw(p->indel, stdout);
+               for (j = 1; j <= p->indel; ++j) {
+                       int 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) {
+               printw(p->indel, stdout);
+               for (j = 1; j <= -p->indel; ++j) {
+                       int c = (ref && (int)pos+j < ref_len)? ref[pos+j] : 'N';
+                       putchar(bam1_strand(p->b)? tolower(c) : toupper(c));
+               }
+       }
        if (p->is_tail) putchar('$');
 }
 

diff --git a/examples/toy.sam b/examples/toy.sam
index c6d0570277cdaf3877c0388a37f88c58baa1659b..1aff2204e4078b9adc9943a9533a2626cca887c8 100644 (file)
--- a/examples/toy.sam
+++ b/examples/toy.sam
@@ -1,7 +1,7 @@
 @SQ    SN:ref  LN:45
 @SQ    SN:ref2 LN:40
-r001   163     ref     7       30      8M2I4M1D3M      =       37      39      TTAGATAAAGGATACTG       *
-r002   0       ref     9       30      1S2I6M1P1I4M2I  *       0       0       AAAAGATAAGGATAAA        *
+r001   163     ref     7       30      8M4I4M1D3M      =       37      39      TTAGATAAAGAGGATACTG     *
+r002   0       ref     9       30      1S2I6M1P1I1P1I4M2I      *       0       0       AAAAGATAAGGGATAAA       *
 r003   0       ref     9       30      5H6M    *       0       0       AGCTAA  *
 r004   0       ref     16      30      6M14N1I5M       *       0       0       ATAGCTCTCAGC    *
 r003   16      ref     29      30      6H5M    *       0       0       TAGGC   *