From: Heng Li <lh3@live.co.uk>
Date: Wed, 27 Oct 2010 13:58:56 +0000 (+0000)
Subject: improve pileup a little bit
X-Git-Url: https://git.donarmstrong.com/?p=samtools.git;a=commitdiff_plain;h=6ee6fd65332260af63fbfbf0a625ea3847371eda

improve pileup a little bit
---

diff --git a/bam_pileup.c b/bam_pileup.c
index 62d5fa1..a49f795 100644
--- 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 cf98f69..f0e63ce 100644
--- 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 c6d0570..1aff220 100644
--- 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	*