Merge commit 'upstream/0.1.18'
[samtools.git] / bam_color.c
1 #include <ctype.h>
2 #include "bam.h"
3
4 /*!
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
8  @return       color
9
10  @discussion   Returns 0 no color information is found.
11  */
12 char bam_aux_getCSi(bam1_t *b, int i)
13 {
14         uint8_t *c = bam_aux_get(b, "CS");
15         char *cs = NULL;
16
17         // return the base if the tag was not found
18         if(0 == c) return 0;
19
20         cs = bam_aux2Z(c);
21         // adjust for strandedness and leading adaptor
22         if(bam1_strand(b)) i = strlen(cs) - 1 - i;
23         else i++;
24         return cs[i];
25 }
26
27 /*!
28  @abstract     Get the color quality of the color encoding the previous and current base
29  @param b      pointer to an alignment
30  @param i      The i-th position, 0-based
31  @return       color quality
32
33  @discussion   Returns 0 no color information is found.
34  */
35 char bam_aux_getCQi(bam1_t *b, int i)
36 {
37         uint8_t *c = bam_aux_get(b, "CQ");
38         char *cq = NULL;
39         
40         // return the base if the tag was not found
41         if(0 == c) return 0;
42
43         cq = bam_aux2Z(c);
44         // adjust for strandedness
45         if(bam1_strand(b)) i = strlen(cq) - 1 - i;
46         return cq[i];
47 }
48
49 char bam_aux_nt2int(char a)
50 {
51         switch(toupper(a)) {
52                 case 'A':
53                         return 0;
54                         break;
55                 case 'C':
56                         return 1;
57                         break;
58                 case 'G':
59                         return 2;
60                         break;
61                 case 'T':
62                         return 3;
63                         break;
64                 default:
65                         return 4;
66                         break;
67         }
68 }
69
70 char bam_aux_ntnt2cs(char a, char b)
71 {
72         a = bam_aux_nt2int(a);
73         b = bam_aux_nt2int(b);
74         if(4 == a || 4 == b) return '4';
75         return "0123"[(int)(a ^ b)];
76 }
77
78 /*!
79  @abstract     Get the color error profile at the give position    
80  @param b      pointer to an alignment
81  @return       the original color if the color was an error, '-' (dash) otherwise
82
83  @discussion   Returns 0 no color information is found.
84  */
85 char bam_aux_getCEi(bam1_t *b, int i)
86 {
87         int cs_i;
88         uint8_t *c = bam_aux_get(b, "CS");
89         char *cs = NULL;
90         char prev_b, cur_b;
91         char cur_color, cor_color;
92
93         // return the base if the tag was not found
94         if(0 == c) return 0;
95         
96         cs = bam_aux2Z(c);
97
98         // adjust for strandedness and leading adaptor
99         if(bam1_strand(b)) { //reverse strand
100                 cs_i = strlen(cs) - 1 - i;
101                 // get current color
102                 cur_color = cs[cs_i];
103                 // get previous base.  Note: must rc adaptor
104                 prev_b = (cs_i == 1) ? "TGCAN"[(int)bam_aux_nt2int(cs[0])] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i+1)];
105                 // get current base
106                 cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)]; 
107         }
108         else {
109                 cs_i=i+1;
110                 // get current color
111                 cur_color = cs[cs_i];
112                 // get previous base
113                 prev_b = (0 == i) ? cs[0] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i-1)];
114                 // get current base
115                 cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)];
116         }
117
118         // corrected color
119         cor_color = bam_aux_ntnt2cs(prev_b, cur_b);
120
121         if(cur_color == cor_color) { 
122                 return '-';
123         }
124         else {
125                 return cur_color;
126         }
127 }