|
0
|
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 }
|
