bamcheck: new stats, number of filtered vs raw sequences

[samtools.git] / bam_color.c

diff --git a/bam_color.c b/bam_color.c
index 75aedd603a1798e85a40c982dd9392ee2a3e5b3f..8b86e2fb6fb0a328b4a903e6dd2d463c3a0114c9 100644 (file)
--- a/bam_color.c
+++ b/bam_color.c
@@ -19,8 +19,14 @@ char bam_aux_getCSi(bam1_t *b, int i)
 
        cs = bam_aux2Z(c);
        // adjust for strandedness and leading adaptor
-       if(bam1_strand(b)) i = strlen(cs) - 1 - i;
-       else i++;
+       if(bam1_strand(b)) {
+           i = strlen(cs) - 1 - i;
+           // adjust for leading hard clip
+           uint32_t cigar = bam1_cigar(b)[0];
+           if((cigar & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) {
+               i -= cigar >> BAM_CIGAR_SHIFT;
+           }
+       } else { i++; }
        return cs[i];
 }
 
@@ -42,7 +48,14 @@ char bam_aux_getCQi(bam1_t *b, int i)
 
        cq = bam_aux2Z(c);
        // adjust for strandedness
-       if(bam1_strand(b)) i = strlen(cq) - 1 - i;
+       if(bam1_strand(b)) {
+           i = strlen(cq) - 1 - i;
+           // adjust for leading hard clip
+           uint32_t cigar = bam1_cigar(b)[0];
+           if((cigar & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) {
+               i -= (cigar >> BAM_CIGAR_SHIFT);
+           }
+       }
        return cq[i];
 }
 
@@ -98,10 +111,15 @@ char bam_aux_getCEi(bam1_t *b, int i)
        // adjust for strandedness and leading adaptor
        if(bam1_strand(b)) { //reverse strand
                cs_i = strlen(cs) - 1 - i;
+               // adjust for leading hard clip
+               uint32_t cigar = bam1_cigar(b)[0];
+               if((cigar & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) {
+                       cs_i -= cigar >> BAM_CIGAR_SHIFT;
+               }
                // get current color
                cur_color = cs[cs_i];
-               // get previous base
-               prev_b = (0 == cs_i) ? cs[0] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i+1)];
+               // get previous base.  Note: must rc adaptor
+               prev_b = (cs_i == 1) ? "TGCAN"[(int)bam_aux_nt2int(cs[0])] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i+1)];
                // get current base
                cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)]; 
        }