Mercurial > repos > lsong10 > psiclass

comparison PsiCLASS-1.0.2/samtools-0.1.19/bam_color.c @ 0:903fc43d6227 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author lsong10
date Fri, 26 Mar 2021 16:52:45 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:903fc43d6227
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)) {
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;
28 }
29 } else { i++; }
30 return cs[i];
31 }
32
33 /*!
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
37 @return color quality
38
39 @discussion Returns 0 no color information is found.
40 */
41 char bam_aux_getCQi(bam1_t *b, int i)
42 {
43 uint8_t *c = bam_aux_get(b, "CQ");
44 char *cq = NULL;
45
46 // return the base if the tag was not found
47 if(0 == c) return 0;
48
49 cq = bam_aux2Z(c);
50 // adjust for strandedness
51 if(bam1_strand(b)) {
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);
57 }
58 }
59 return cq[i];
60 }
61
62 char bam_aux_nt2int(char a)
63 {
64 switch(toupper(a)) {
65 case 'A':
66 return 0;
67 break;
68 case 'C':
69 return 1;
70 break;
71 case 'G':
72 return 2;
73 break;
74 case 'T':
75 return 3;
76 break;
77 default:
78 return 4;
79 break;
80 }
81 }
82
83 char bam_aux_ntnt2cs(char a, char b)
84 {
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)];
89 }
90
91 /*!
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
95
96 @discussion Returns 0 no color information is found.
97 */
98 char bam_aux_getCEi(bam1_t *b, int i)
99 {
100 int cs_i;
101 uint8_t *c = bam_aux_get(b, "CS");
102 char *cs = NULL;
103 char prev_b, cur_b;
104 char cur_color, cor_color;
105
106 // return the base if the tag was not found
107 if(0 == c) return 0;
108
109 cs = bam_aux2Z(c);
110
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;
118 }
119 // get current color
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)];
123 // get current base
124 cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)];
125 }
126 else {
127 cs_i=i+1;
128 // get current color
129 cur_color = cs[cs_i];
130 // get previous base
131 prev_b = (0 == i) ? cs[0] : bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i-1)];
132 // get current base
133 cur_b = bam_nt16_rev_table[bam1_seqi(bam1_seq(b), i)];
134 }
135
136 // corrected color
137 cor_color = bam_aux_ntnt2cs(prev_b, cur_b);
138
139 if(cur_color == cor_color) {
140 return '-';
141 }
142 else {
143 return cur_color;
144 }
145 }