5 @abstract Get the color encoding the previous and current base
6 @param b pointer to an alignment
7 @param i The i-th position, 0-based
10 @discussion Returns 0 no color information is found.
12 char bam_aux_getCSi(bam1_t *b, int i)
14 uint8_t *c = bam_aux_get(b, "CS");
17 // return the base if the tag was not found
21 // adjust for strandedness and leading adaptor
23 i = strlen(cs) - 1 - i;
24 // adjust for leading hard clip
25 uint32_t cigar = bam1_cigar(b)[0];
26 if((cigar & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) {
27 i -= cigar >> BAM_CIGAR_SHIFT;
34 @abstract Get the color quality of the color encoding the previous and current base
35 @param b pointer to an alignment
36 @param i The i-th position, 0-based
39 @discussion Returns 0 no color information is found.
41 char bam_aux_getCQi(bam1_t *b, int i)
43 uint8_t *c = bam_aux_get(b, "CQ");
46 // return the base if the tag was not found
50 // adjust for strandedness
52 i = strlen(cq) - 1 - i;
53 // adjust for leading hard clip
54 uint32_t cigar = bam1_cigar(b)[0];
55 if((cigar & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) {
56 i -= (cigar >> BAM_CIGAR_SHIFT);
62 char bam_aux_nt2int(char a)
83 char bam_aux_ntnt2cs(char a, char b)
85 a = bam_aux_nt2int(a);
86 b = bam_aux_nt2int(b);
87 if(4 == a || 4 == b) return '4';
88 return "0123"[(int)(a ^ b)];
92 @abstract Get the color error profile at the give position
93 @param b pointer to an alignment
94 @return the original color if the color was an error, '-' (dash) otherwise
96 @discussion Returns 0 no color information is found.
98 char bam_aux_getCEi(bam1_t *b, int i)
101 uint8_t *c = bam_aux_get(b, "CS");
104 char cur_color, cor_color;
106 // return the base if the tag was not found
111 // adjust for strandedness and leading adaptor
112 if(bam1_strand(b)) { //reverse strand
113 cs_i = strlen(cs) - 1 - i;
114 // adjust for leading hard clip
115 uint32_t cigar = bam1_cigar(b)[0];
116 if((cigar & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) {
117 cs_i -= cigar >> BAM_CIGAR_SHIFT;
120 cur_color = cs[cs_i];
121 // get previous base. Note: must rc adaptor
122 prev_b = (cs_i == 1) ? "TGCAN"[(int)bam_aux_nt2int(cs[0])] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i+1)];
124 cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)];
129 cur_color = cs[cs_i];
131 prev_b = (0 == i) ? cs[0] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i-1)];
133 cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)];
137 cor_color = bam_aux_ntnt2cs(prev_b, cur_b);
139 if(cur_color == cor_color) {