|
0
|
1 /* The MIT License
|
|
|
2
|
|
|
3 Copyright (c) 2003-2006, 2008-2010, by Heng Li <lh3lh3@live.co.uk>
|
|
|
4
|
|
|
5 Permission is hereby granted, free of charge, to any person obtaining
|
|
|
6 a copy of this software and associated documentation files (the
|
|
|
7 "Software"), to deal in the Software without restriction, including
|
|
|
8 without limitation the rights to use, copy, modify, merge, publish,
|
|
|
9 distribute, sublicense, and/or sell copies of the Software, and to
|
|
|
10 permit persons to whom the Software is furnished to do so, subject to
|
|
|
11 the following conditions:
|
|
|
12
|
|
|
13 The above copyright notice and this permission notice shall be
|
|
|
14 included in all copies or substantial portions of the Software.
|
|
|
15
|
|
|
16 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|
|
17 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
|
|
|
18 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|
|
19 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
|
|
|
20 BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
|
|
|
21 ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
|
|
|
22 CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
|
|
23 SOFTWARE.
|
|
|
24 */
|
|
|
25
|
|
|
26 #include <stdlib.h>
|
|
|
27 #include <stdio.h>
|
|
|
28 #include <string.h>
|
|
|
29 #include <stdint.h>
|
|
|
30 #include <math.h>
|
|
|
31 #include "kprobaln.h"
|
|
|
32
|
|
|
33 /*****************************************
|
|
|
34 * Probabilistic banded glocal alignment *
|
|
|
35 *****************************************/
|
|
|
36
|
|
|
37 #define EI .25
|
|
|
38 #define EM .33333333333
|
|
|
39
|
|
|
40 static float g_qual2prob[256];
|
|
|
41
|
|
|
42 #define set_u(u, b, i, k) { int x=(i)-(b); x=x>0?x:0; (u)=((k)-x+1)*3; }
|
|
|
43
|
|
|
44 kpa_par_t kpa_par_def = { 0.001, 0.1, 10 };
|
|
|
45 kpa_par_t kpa_par_alt = { 0.0001, 0.01, 10 };
|
|
|
46
|
|
|
47 /*
|
|
|
48 The topology of the profile HMM:
|
|
|
49
|
|
|
50 /\ /\ /\ /\
|
|
|
51 I[1] I[k-1] I[k] I[L]
|
|
|
52 ^ \ \ ^ \ ^ \ \ ^
|
|
|
53 | \ \ | \ | \ \ |
|
|
|
54 M[0] M[1] -> ... -> M[k-1] -> M[k] -> ... -> M[L] M[L+1]
|
|
|
55 \ \/ \/ \/ /
|
|
|
56 \ /\ /\ /\ /
|
|
|
57 -> D[k-1] -> D[k] ->
|
|
|
58
|
|
|
59 M[0] points to every {M,I}[k] and every {M,I}[k] points M[L+1].
|
|
|
60
|
|
|
61 On input, _ref is the reference sequence and _query is the query
|
|
|
62 sequence. Both are sequences of 0/1/2/3/4 where 4 stands for an
|
|
|
63 ambiguous residue. iqual is the base quality. c sets the gap open
|
|
|
64 probability, gap extension probability and band width.
|
|
|
65
|
|
|
66 On output, state and q are arrays of length l_query. The higher 30
|
|
|
67 bits give the reference position the query base is matched to and the
|
|
|
68 lower two bits can be 0 (an alignment match) or 1 (an
|
|
|
69 insertion). q[i] gives the phred scaled posterior probability of
|
|
|
70 state[i] being wrong.
|
|
|
71 */
|
|
|
72 int kpa_glocal(const uint8_t *_ref, int l_ref, const uint8_t *_query, int l_query, const uint8_t *iqual,
|
|
|
73 const kpa_par_t *c, int *state, uint8_t *q)
|
|
|
74 {
|
|
|
75 double **f, **b = 0, *s, m[9], sI, sM, bI, bM, pb;
|
|
|
76 float *qual, *_qual;
|
|
|
77 const uint8_t *ref, *query;
|
|
|
78 int bw, bw2, i, k, is_diff = 0, is_backward = 1, Pr;
|
|
|
79
|
|
|
80 if ( l_ref<=0 || l_query<=0 ) return 0; // FIXME: this may not be an ideal fix, just prevents sefgault
|
|
|
81
|
|
|
82 /*** initialization ***/
|
|
|
83 is_backward = state && q? 1 : 0;
|
|
|
84 ref = _ref - 1; query = _query - 1; // change to 1-based coordinate
|
|
|
85 bw = l_ref > l_query? l_ref : l_query;
|
|
|
86 if (bw > c->bw) bw = c->bw;
|
|
|
87 if (bw < abs(l_ref - l_query)) bw = abs(l_ref - l_query);
|
|
|
88 bw2 = bw * 2 + 1;
|
|
|
89 // allocate the forward and backward matrices f[][] and b[][] and the scaling array s[]
|
|
|
90 f = calloc(l_query+1, sizeof(void*));
|
|
|
91 if (is_backward) b = calloc(l_query+1, sizeof(void*));
|
|
|
92 for (i = 0; i <= l_query; ++i) { // FIXME: this will lead in segfault for l_query==0
|
|
|
93 f[i] = calloc(bw2 * 3 + 6, sizeof(double)); // FIXME: this is over-allocated for very short seqs
|
|
|
94 if (is_backward) b[i] = calloc(bw2 * 3 + 6, sizeof(double));
|
|
|
95 }
|
|
|
96 s = calloc(l_query+2, sizeof(double)); // s[] is the scaling factor to avoid underflow
|
|
|
97 // initialize qual
|
|
|
98 _qual = calloc(l_query, sizeof(float));
|
|
|
99 if (g_qual2prob[0] == 0)
|
|
|
100 for (i = 0; i < 256; ++i)
|
|
|
101 g_qual2prob[i] = pow(10, -i/10.);
|
|
|
102 for (i = 0; i < l_query; ++i) _qual[i] = g_qual2prob[iqual? iqual[i] : 30];
|
|
|
103 qual = _qual - 1;
|
|
|
104 // initialize transition probability
|
|
|
105 sM = sI = 1. / (2 * l_query + 2); // the value here seems not to affect results; FIXME: need proof
|
|
|
106 m[0*3+0] = (1 - c->d - c->d) * (1 - sM); m[0*3+1] = m[0*3+2] = c->d * (1 - sM);
|
|
|
107 m[1*3+0] = (1 - c->e) * (1 - sI); m[1*3+1] = c->e * (1 - sI); m[1*3+2] = 0.;
|
|
|
108 m[2*3+0] = 1 - c->e; m[2*3+1] = 0.; m[2*3+2] = c->e;
|
|
|
109 bM = (1 - c->d) / l_ref; bI = c->d / l_ref; // (bM+bI)*l_ref==1
|
|
|
110 /*** forward ***/
|
|
|
111 // f[0]
|
|
|
112 set_u(k, bw, 0, 0);
|
|
|
113 f[0][k] = s[0] = 1.;
|
|
|
114 { // f[1]
|
|
|
115 double *fi = f[1], sum;
|
|
|
116 int beg = 1, end = l_ref < bw + 1? l_ref : bw + 1, _beg, _end;
|
|
|
117 for (k = beg, sum = 0.; k <= end; ++k) {
|
|
|
118 int u;
|
|
|
119 double e = (ref[k] > 3 || query[1] > 3)? 1. : ref[k] == query[1]? 1. - qual[1] : qual[1] * EM;
|
|
|
120 set_u(u, bw, 1, k);
|
|
|
121 fi[u+0] = e * bM; fi[u+1] = EI * bI;
|
|
|
122 sum += fi[u] + fi[u+1];
|
|
|
123 }
|
|
|
124 // rescale
|
|
|
125 s[1] = sum;
|
|
|
126 set_u(_beg, bw, 1, beg); set_u(_end, bw, 1, end); _end += 2;
|
|
|
127 for (k = _beg; k <= _end; ++k) fi[k] /= sum;
|
|
|
128 }
|
|
|
129 // f[2..l_query]
|
|
|
130 for (i = 2; i <= l_query; ++i) {
|
|
|
131 double *fi = f[i], *fi1 = f[i-1], sum, qli = qual[i];
|
|
|
132 int beg = 1, end = l_ref, x, _beg, _end;
|
|
|
133 uint8_t qyi = query[i];
|
|
|
134 x = i - bw; beg = beg > x? beg : x; // band start
|
|
|
135 x = i + bw; end = end < x? end : x; // band end
|
|
|
136 for (k = beg, sum = 0.; k <= end; ++k) {
|
|
|
137 int u, v11, v01, v10;
|
|
|
138 double e;
|
|
|
139 e = (ref[k] > 3 || qyi > 3)? 1. : ref[k] == qyi? 1. - qli : qli * EM;
|
|
|
140 set_u(u, bw, i, k); set_u(v11, bw, i-1, k-1); set_u(v10, bw, i-1, k); set_u(v01, bw, i, k-1);
|
|
|
141 fi[u+0] = e * (m[0] * fi1[v11+0] + m[3] * fi1[v11+1] + m[6] * fi1[v11+2]);
|
|
|
142 fi[u+1] = EI * (m[1] * fi1[v10+0] + m[4] * fi1[v10+1]);
|
|
|
143 fi[u+2] = m[2] * fi[v01+0] + m[8] * fi[v01+2];
|
|
|
144 sum += fi[u] + fi[u+1] + fi[u+2];
|
|
|
145 // fprintf(stderr, "F (%d,%d;%d): %lg,%lg,%lg\n", i, k, u, fi[u], fi[u+1], fi[u+2]); // DEBUG
|
|
|
146 }
|
|
|
147 // rescale
|
|
|
148 s[i] = sum;
|
|
|
149 set_u(_beg, bw, i, beg); set_u(_end, bw, i, end); _end += 2;
|
|
|
150 for (k = _beg, sum = 1./sum; k <= _end; ++k) fi[k] *= sum;
|
|
|
151 }
|
|
|
152 { // f[l_query+1]
|
|
|
153 double sum;
|
|
|
154 for (k = 1, sum = 0.; k <= l_ref; ++k) {
|
|
|
155 int u;
|
|
|
156 set_u(u, bw, l_query, k);
|
|
|
157 if (u < 3 || u >= bw2*3+3) continue;
|
|
|
158 sum += f[l_query][u+0] * sM + f[l_query][u+1] * sI;
|
|
|
159 }
|
|
|
160 s[l_query+1] = sum; // the last scaling factor
|
|
|
161 }
|
|
|
162 { // compute likelihood
|
|
|
163 double p = 1., Pr1 = 0.;
|
|
|
164 for (i = 0; i <= l_query + 1; ++i) {
|
|
|
165 p *= s[i];
|
|
|
166 if (p < 1e-100) Pr1 += -4.343 * log(p), p = 1.;
|
|
|
167 }
|
|
|
168 Pr1 += -4.343 * log(p * l_ref * l_query);
|
|
|
169 Pr = (int)(Pr1 + .499);
|
|
|
170 if (!is_backward) { // skip backward and MAP
|
|
|
171 for (i = 0; i <= l_query; ++i) free(f[i]);
|
|
|
172 free(f); free(s); free(_qual);
|
|
|
173 return Pr;
|
|
|
174 }
|
|
|
175 }
|
|
|
176 /*** backward ***/
|
|
|
177 // b[l_query] (b[l_query+1][0]=1 and thus \tilde{b}[][]=1/s[l_query+1]; this is where s[l_query+1] comes from)
|
|
|
178 for (k = 1; k <= l_ref; ++k) {
|
|
|
179 int u;
|
|
|
180 double *bi = b[l_query];
|
|
|
181 set_u(u, bw, l_query, k);
|
|
|
182 if (u < 3 || u >= bw2*3+3) continue;
|
|
|
183 bi[u+0] = sM / s[l_query] / s[l_query+1]; bi[u+1] = sI / s[l_query] / s[l_query+1];
|
|
|
184 }
|
|
|
185 // b[l_query-1..1]
|
|
|
186 for (i = l_query - 1; i >= 1; --i) {
|
|
|
187 int beg = 1, end = l_ref, x, _beg, _end;
|
|
|
188 double *bi = b[i], *bi1 = b[i+1], y = (i > 1), qli1 = qual[i+1];
|
|
|
189 uint8_t qyi1 = query[i+1];
|
|
|
190 x = i - bw; beg = beg > x? beg : x;
|
|
|
191 x = i + bw; end = end < x? end : x;
|
|
|
192 for (k = end; k >= beg; --k) {
|
|
|
193 int u, v11, v01, v10;
|
|
|
194 double e;
|
|
|
195 set_u(u, bw, i, k); set_u(v11, bw, i+1, k+1); set_u(v10, bw, i+1, k); set_u(v01, bw, i, k+1);
|
|
|
196 e = (k >= l_ref? 0 : (ref[k+1] > 3 || qyi1 > 3)? 1. : ref[k+1] == qyi1? 1. - qli1 : qli1 * EM) * bi1[v11];
|
|
|
197 bi[u+0] = e * m[0] + EI * m[1] * bi1[v10+1] + m[2] * bi[v01+2]; // bi1[v11] has been foled into e.
|
|
|
198 bi[u+1] = e * m[3] + EI * m[4] * bi1[v10+1];
|
|
|
199 bi[u+2] = (e * m[6] + m[8] * bi[v01+2]) * y;
|
|
|
200 // fprintf(stderr, "B (%d,%d;%d): %lg,%lg,%lg\n", i, k, u, bi[u], bi[u+1], bi[u+2]); // DEBUG
|
|
|
201 }
|
|
|
202 // rescale
|
|
|
203 set_u(_beg, bw, i, beg); set_u(_end, bw, i, end); _end += 2;
|
|
|
204 for (k = _beg, y = 1./s[i]; k <= _end; ++k) bi[k] *= y;
|
|
|
205 }
|
|
|
206 { // b[0]
|
|
|
207 int beg = 1, end = l_ref < bw + 1? l_ref : bw + 1;
|
|
|
208 double sum = 0.;
|
|
|
209 for (k = end; k >= beg; --k) {
|
|
|
210 int u;
|
|
|
211 double e = (ref[k] > 3 || query[1] > 3)? 1. : ref[k] == query[1]? 1. - qual[1] : qual[1] * EM;
|
|
|
212 set_u(u, bw, 1, k);
|
|
|
213 if (u < 3 || u >= bw2*3+3) continue;
|
|
|
214 sum += e * b[1][u+0] * bM + EI * b[1][u+1] * bI;
|
|
|
215 }
|
|
|
216 set_u(k, bw, 0, 0);
|
|
|
217 pb = b[0][k] = sum / s[0]; // if everything works as is expected, pb == 1.0
|
|
|
218 }
|
|
|
219 is_diff = fabs(pb - 1.) > 1e-7? 1 : 0;
|
|
|
220 /*** MAP ***/
|
|
|
221 for (i = 1; i <= l_query; ++i) {
|
|
|
222 double sum = 0., *fi = f[i], *bi = b[i], max = 0.;
|
|
|
223 int beg = 1, end = l_ref, x, max_k = -1;
|
|
|
224 x = i - bw; beg = beg > x? beg : x;
|
|
|
225 x = i + bw; end = end < x? end : x;
|
|
|
226 for (k = beg; k <= end; ++k) {
|
|
|
227 int u;
|
|
|
228 double z;
|
|
|
229 set_u(u, bw, i, k);
|
|
|
230 z = fi[u+0] * bi[u+0]; if (z > max) max = z, max_k = (k-1)<<2 | 0; sum += z;
|
|
|
231 z = fi[u+1] * bi[u+1]; if (z > max) max = z, max_k = (k-1)<<2 | 1; sum += z;
|
|
|
232 }
|
|
|
233 max /= sum; sum *= s[i]; // if everything works as is expected, sum == 1.0
|
|
|
234 if (state) state[i-1] = max_k;
|
|
|
235 if (q) k = (int)(-4.343 * log(1. - max) + .499), q[i-1] = k > 100? 99 : k;
|
|
|
236 #ifdef _MAIN
|
|
|
237 fprintf(stderr, "(%.10lg,%.10lg) (%d,%d:%c,%c:%d) %lg\n", pb, sum, i-1, max_k>>2,
|
|
|
238 "ACGT"[query[i]], "ACGT"[ref[(max_k>>2)+1]], max_k&3, max); // DEBUG
|
|
|
239 #endif
|
|
|
240 }
|
|
|
241 /*** free ***/
|
|
|
242 for (i = 0; i <= l_query; ++i) {
|
|
|
243 free(f[i]); free(b[i]);
|
|
|
244 }
|
|
|
245 free(f); free(b); free(s); free(_qual);
|
|
|
246 return Pr;
|
|
|
247 }
|
|
|
248
|
|
|
249 #ifdef _MAIN
|
|
|
250 #include <unistd.h>
|
|
|
251 int main(int argc, char *argv[])
|
|
|
252 {
|
|
|
253 uint8_t conv[256], *iqual, *ref, *query;
|
|
|
254 int c, l_ref, l_query, i, q = 30, b = 10, P;
|
|
|
255 while ((c = getopt(argc, argv, "b:q:")) >= 0) {
|
|
|
256 switch (c) {
|
|
|
257 case 'b': b = atoi(optarg); break;
|
|
|
258 case 'q': q = atoi(optarg); break;
|
|
|
259 }
|
|
|
260 }
|
|
|
261 if (optind + 2 > argc) {
|
|
|
262 fprintf(stderr, "Usage: %s [-q %d] [-b %d] <ref> <query>\n", argv[0], q, b); // example: acttc attc
|
|
|
263 return 1;
|
|
|
264 }
|
|
|
265 memset(conv, 4, 256);
|
|
|
266 conv['a'] = conv['A'] = 0; conv['c'] = conv['C'] = 1;
|
|
|
267 conv['g'] = conv['G'] = 2; conv['t'] = conv['T'] = 3;
|
|
|
268 ref = (uint8_t*)argv[optind]; query = (uint8_t*)argv[optind+1];
|
|
|
269 l_ref = strlen((char*)ref); l_query = strlen((char*)query);
|
|
|
270 for (i = 0; i < l_ref; ++i) ref[i] = conv[ref[i]];
|
|
|
271 for (i = 0; i < l_query; ++i) query[i] = conv[query[i]];
|
|
|
272 iqual = malloc(l_query);
|
|
|
273 memset(iqual, q, l_query);
|
|
|
274 kpa_par_def.bw = b;
|
|
|
275 P = kpa_glocal(ref, l_ref, query, l_query, iqual, &kpa_par_alt, 0, 0);
|
|
|
276 fprintf(stderr, "%d\n", P);
|
|
|
277 free(iqual);
|
|
|
278 return 0;
|
|
|
279 }
|
|
|
280 #endif
|
