comparison pyPRADA_1.2/tools/bwa-0.5.7-mh/cs2nt.c @ 0:acc2ca1a3ba4

Uploaded
author siyuan
date Thu, 20 Feb 2014 00:44:58 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:acc2ca1a3ba4
1 #include <string.h>
2 #include <stdint.h>
3 #include <stdlib.h>
4 #include "bwtaln.h"
5 #include "stdaln.h"
6
7 /*
8 Here is a delicate example. ref_nt=ATTAAC(RBRBG), read_cs=RBBOG. If we
9 decode as ATTGAC(RBGOG), there are one color change and one nt change;
10 if we decode as ATTAAC(RBRBG), there are two color changes.
11
12 In DP, if color quality is smaller than COLOR_MM, we will use COLOR_MM
13 as the penalty; otherwise, we will use color quality as the
14 penalty. This means we always prefer two consistent color changes over
15 a nt change, but if a color has high quality, we may prefer one nt
16 change.
17
18 In the above example, the penalties of the two types of decoding are
19 q(B)+25 and q(B)+q(O), respectively. If q(O)>25, we prefer the first;
20 otherwise the second. Note that no matter what we choose, the fourth
21 base will get a low nt quality.
22 */
23
24 #define COLOR_MM 19
25 #define NUCL_MM 25
26
27 static const int nst_ntnt2cs_table[] = { 4, 0, 0, 1, 0, 2, 3, 4, 0, 3, 2, 4, 1, 4, 4, 4 };
28
29 /*
30 {A,C,G,T,N} -> {0,1,2,3,4}
31 nt_ref[0..size]: nucleotide reference: 0/1/2/3/4
32 cs_read[0..size-1]: color read+qual sequence: base<<6|qual; qual==63 for N
33 nt_read[0..size]: nucleotide read sequence: 0/1/2/3 (returned)
34 btarray[0..4*size]: backtrack array (working space)
35 */
36 void cs2nt_DP(int size, const uint8_t *nt_ref, const uint8_t *cs_read, uint8_t *nt_read, uint8_t *btarray)
37 {
38 int h[8], curr, last;
39 int x, y, xmin, hmin, k;
40
41 // h[0..3] and h[4..7] are the current and last best score array, depending on curr and last
42
43 // recursion: initial value
44 if (nt_ref[0] >= 4) memset(h, 0, sizeof(int) << 2);
45 else {
46 for (x = 0; x != 4; ++x) h[x] = NUCL_MM;
47 h[nt_ref[0]] = 0;
48 }
49 // recursion: main loop
50 curr = 1; last = 0;
51 for (k = 1; k <= size; ++k) {
52 for (x = 0; x != 4; ++x) {
53 int min = 0x7fffffff, ymin = 0;
54 for (y = 0; y != 4; ++y) {
55 int s = h[last<<2|y];
56 if ((cs_read[k-1]&0x3f) != 63 && cs_read[k-1]>>6 != nst_ntnt2cs_table[1<<x|1<<y])
57 s += ((cs_read[k-1]&0x3f) < COLOR_MM)? COLOR_MM : (cs_read[k-1]&0x3f); // color mismatch
58 if (nt_ref[k] < 4 && nt_ref[k] != x) s += NUCL_MM; // nt mismatch
59 if (s < min) {
60 min = s; ymin = y;
61 }
62 }
63 h[curr<<2|x] = min; btarray[k<<2|x] = ymin;
64 }
65 last = curr; curr = 1 - curr; // swap
66 }
67 // back trace
68 hmin = 0x7fffffff; xmin = 0;
69 for (x = 0; x != 4; ++x) {
70 if (h[last<<2|x] < hmin) {
71 hmin = h[last<<2|x]; xmin = x;
72 }
73 }
74 nt_read[size] = xmin;
75 for (k = size - 1; k >= 0; --k)
76 nt_read[k] = btarray[(k+1)<<2 | nt_read[k+1]];
77 }
78 /*
79 nt_read[0..size]: nucleotide read sequence: 0/1/2/3
80 cs_read[0..size-1]: color read+qual sequence: base<<6|qual; qual==63 for N
81 tarray[0..size*2-1]: temporary array
82 */
83 uint8_t *cs2nt_nt_qual(int size, const uint8_t *nt_read, const uint8_t *cs_read, uint8_t *tarray)
84 {
85 int k, c1, c2;
86 uint8_t *t2array = tarray + size;
87 // get the color sequence of nt_read
88 c1 = nt_read[0];
89 for (k = 1; k <= size; ++k) {
90 c2 = nt_read[k]; // in principle, there is no 'N' in nt_read[]; just in case
91 tarray[k-1] = (c1 >= 4 || c2 >= 4)? 4 : nst_ntnt2cs_table[1<<c1 | 1<<c2];
92 c1 = c2;
93 }
94 for (k = 1; k != size; ++k) {
95 int q = 0;
96 if (tarray[k-1] == cs_read[k-1]>>6 && tarray[k] == cs_read[k]>>6) {
97 q = (int)(cs_read[k-1]&0x3f) + (int)(cs_read[k]&0x3f) + 10;
98 } else if (tarray[k-1] == cs_read[k-1]>>6) {
99 q = (int)(cs_read[k-1]&0x3f) - (int)(cs_read[k]&0x3f);
100 } else if (tarray[k] == cs_read[k]>>6) {
101 q = (int)(cs_read[k]&0x3f) - (int)(cs_read[k-1]&0x3f);
102 } // else, q = 0
103 if (q < 0) q = 0;
104 if (q > 60) q = 60;
105 t2array[k] = nt_read[k]<<6 | q;
106 if ((cs_read[k-1]&0x3f) == 63 || (cs_read[k]&0x3f) == 63) t2array[k] = 0;
107 }
108 return t2array + 1; // of size-2
109 }
110
111 // this function will be called when p->seq has been reversed by refine_gapped()
112 void bwa_cs2nt_core(bwa_seq_t *p, bwtint_t l_pac, ubyte_t *pac)
113 {
114 uint8_t *ta, *nt_read, *btarray, *tarray, *nt_ref, *cs_read, *new_nt_read;
115 int i, len;
116 uint8_t *seq;
117
118 // set temporary arrays
119 if (p->type == BWA_TYPE_NO_MATCH) return;
120 len = p->len + p->n_gapo + p->n_gape + 100; // leave enough space
121 ta = (uint8_t*)malloc(len * 7);
122 nt_ref = ta;
123 cs_read = nt_ref + len;
124 nt_read = cs_read + len;
125 btarray = nt_read + len;
126 tarray = nt_read + len;
127
128 #define __gen_csbase(_cs, _i, _seq) do { \
129 int q = p->qual[p->strand? p->len - 1 - (_i) : (_i)] - 33; \
130 if (q > 60) q = 60; \
131 if (_seq[_i] > 3) q = 63; \
132 (_cs) = _seq[_i]<<6 | q; \
133 } while (0)
134
135 // generate len, nt_ref[] and cs_read
136 seq = p->strand? p->rseq : p->seq;
137 nt_ref[0] = p->pos? bns_pac(pac, p->pos-1) : 4;
138 if (p->cigar == 0) { // no gap or clipping
139 len = p->len;
140 for (i = 0; i < p->len; ++i) {
141 __gen_csbase(cs_read[i], i, seq);
142 nt_ref[i+1] = bns_pac(pac, p->pos + i);
143 }
144 } else {
145 int k, z;
146 bwtint_t x, y;
147 x = p->pos; y = 0;
148 for (k = z = 0; k < p->n_cigar; ++k) {
149 int l = __cigar_len(p->cigar[k]);
150 if (__cigar_op(p->cigar[k]) == FROM_M) {
151 for (i = 0; i < l; ++i, ++x, ++y) {
152 __gen_csbase(cs_read[z], y, seq);
153 nt_ref[z+1] = bns_pac(pac, x);
154 ++z;
155 }
156 } else if (__cigar_op(p->cigar[k]) == FROM_I) {
157 for (i = 0; i < l; ++i, ++y) {
158 __gen_csbase(cs_read[z], y, seq);
159 nt_ref[z+1] = 4;
160 ++z;
161 }
162 } else if (__cigar_op(p->cigar[k]) == FROM_S) y += l;
163 else x += l;
164 }
165 len = z;
166 }
167
168 cs2nt_DP(len, nt_ref, cs_read, nt_read, btarray);
169 new_nt_read = cs2nt_nt_qual(len, nt_read, cs_read, tarray);
170
171 // update p
172 p->len = p->full_len = len - 1;
173 for (i = 0; i < p->len; ++i) {
174 if ((new_nt_read[i]&0x3f) == 63) {
175 p->qual[i] = 33; seq[i] = 4;
176 } else {
177 p->qual[i] = (new_nt_read[i]&0x3f) + 33;
178 seq[i] = new_nt_read[i]>>6;
179 }
180 }
181 p->qual[p->len] = seq[p->len] = 0;
182 if (p->strand) {
183 memcpy(p->seq, seq, p->len);
184 seq_reverse(p->len, p->seq, 1);
185 seq_reverse(p->len, p->qual, 0);
186 } else {
187 memcpy(p->rseq, seq, p->len);
188 seq_reverse(p->len, p->rseq, 1);
189 }
190 free(ta);
191 }