3 Copyright (c) 2003-2006, 2008-2010, by Heng Li <lh3lh3@live.co.uk>
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:
13 The above copyright notice and this permission notice shall be
14 included in all copies or substantial portions of the Software.
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
33 /*****************************************
34 * Probabilistic banded glocal alignment *
35 *****************************************/
38 #define EM .33333333333
40 static float g_qual2prob[256];
42 #define set_u(u, b, i, k) { int x=(i)-(b); x=x>0?x:0; (u)=((k)-x+1)*3; }
44 kpa_par_t kpa_par_def = { 0.001, 0.1, 10 };
47 The topology of the profile HMM:
53 M[0] M[1] -> ... -> M[k-1] -> M[k] -> ... -> M[L] M[L+1]
58 M[0] points to every {M,I}[k] and every {M,I}[k] points M[L+1].
60 On input, _ref is the reference sequence and _query is the query
61 sequence. Both are sequences of 0/1/2/3/4 where 4 stands for an
62 ambiguous residue. iqual is the base quality. c sets the gap open
63 probability, gap extension probability and band width.
65 On output, state and q are arrays of length l_query. The higher 30
66 bits give the reference position the query base is matched to and the
67 lower two bits can be 0 (an alignment match) or 1 (an
68 insertion). q[i] gives the phred scaled posterior probability of
71 int kpa_glocal(const uint8_t *_ref, int l_ref, const uint8_t *_query, int l_query, const uint8_t *iqual,
72 const kpa_par_t *c, int *state, uint8_t *q)
74 double **f, **b, *s, m[9], sI, sM, bI, bM, pb;
76 const uint8_t *ref, *query;
77 int bw, bw2, i, k, is_diff = 0;
79 /*** initialization ***/
80 ref = _ref - 1; query = _query - 1; // change to 1-based coordinate
81 bw = l_ref > l_query? l_ref : l_query;
82 if (bw > c->bw) bw = c->bw;
83 if (bw < abs(l_ref - l_query)) bw = abs(l_ref - l_query);
85 // allocate the forward and backward matrices f[][] and b[][] and the scaling array s[]
86 f = calloc(l_query+1, sizeof(void*));
87 b = calloc(l_query+1, sizeof(void*));
88 for (i = 0; i <= l_query; ++i) {
89 f[i] = calloc(bw2 * 3 + 6, sizeof(double)); // FIXME: this is over-allocated for very short seqs
90 b[i] = calloc(bw2 * 3 + 6, sizeof(double));
92 s = calloc(l_query+2, sizeof(double)); // s[] is the scaling factor to avoid underflow
94 _qual = calloc(l_query, sizeof(float));
95 if (g_qual2prob[0] == 0)
96 for (i = 0; i < 256; ++i)
97 g_qual2prob[i] = pow(10, -i/10.);
98 for (i = 0; i < l_query; ++i) _qual[i] = g_qual2prob[iqual? iqual[i] : 30];
100 // initialize transition probability
101 sM = sI = 1. / (2 * l_query + 2); // the value here seems not to affect results; FIXME: need proof
102 m[0*3+0] = (1 - c->d - c->d) * (1 - sM); m[0*3+1] = m[0*3+2] = c->d * (1 - sM);
103 m[1*3+0] = (1 - c->e) * (1 - sI); m[1*3+1] = c->e * (1 - sI); m[1*3+2] = 0.;
104 m[2*3+0] = 1 - c->e; m[2*3+1] = 0.; m[2*3+2] = c->e;
105 bM = (1 - c->d) / l_query; bI = c->d / l_query; // (bM+bI)*l_query==1
111 double *fi = f[1], sum;
112 int beg = 1, end = l_ref < bw + 1? l_ref : bw + 1, _beg, _end;
113 for (k = beg, sum = 0.; k <= end; ++k) {
115 double e = (ref[k] > 3 || query[1] > 3)? 1. : ref[k] == query[1]? 1. - qual[1] : qual[1] * EM;
117 fi[u+0] = e * bM; fi[u+1] = EI * bI;
118 sum += fi[u] + fi[u+1];
122 set_u(_beg, bw, 1, beg); set_u(_end, bw, 1, end); _end += 2;
123 for (k = _beg; k <= _end; ++k) fi[k] /= sum;
126 for (i = 2; i <= l_query; ++i) {
127 double *fi = f[i], *fi1 = f[i-1], sum, qli = qual[i];
128 int beg = 1, end = l_ref, x, _beg, _end;
129 uint8_t qyi = query[i];
130 x = i - bw; beg = beg > x? beg : x; // band start
131 x = i + bw; end = end < x? end : x; // band end
132 for (k = beg, sum = 0.; k <= end; ++k) {
133 int u, v11, v01, v10;
135 e = (ref[k] > 3 || qyi > 3)? 1. : ref[k] == qyi? 1. - qli : qli * EM;
136 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);
137 fi[u+0] = e * (m[0] * fi1[v11+0] + m[3] * fi1[v11+1] + m[6] * fi1[v11+2]);
138 fi[u+1] = EI * (m[1] * fi1[v10+0] + m[4] * fi1[v10+1]);
139 fi[u+2] = m[2] * fi[v01+0] + m[8] * fi[v01+2];
140 sum += fi[u] + fi[u+1] + fi[u+2];
141 // fprintf(stderr, "F (%d,%d;%d): %lg,%lg,%lg\n", i, k, u, fi[u], fi[u+1], fi[u+2]); // DEBUG
145 set_u(_beg, bw, i, beg); set_u(_end, bw, i, end); _end += 2;
146 for (k = _beg, sum = 1./sum; k <= _end; ++k) fi[k] *= sum;
150 for (k = 1, sum = 0.; k <= l_ref; ++k) {
152 set_u(u, bw, l_query, k);
153 if (u < 3 || u >= bw2*3+3) continue;
154 sum += f[l_query][u+0] * sM + f[l_query][u+1] * sI;
156 s[l_query+1] = sum; // the last scaling factor
159 // 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)
160 for (k = 1; k <= l_ref; ++k) {
162 double *bi = b[l_query];
163 set_u(u, bw, l_query, k);
164 if (u < 3 || u >= bw2*3+3) continue;
165 bi[u+0] = sM / s[l_query] / s[l_query+1]; bi[u+1] = sI / s[l_query] / s[l_query+1];
168 for (i = l_query - 1; i >= 1; --i) {
169 int beg = 1, end = l_ref, x, _beg, _end;
170 double *bi = b[i], *bi1 = b[i+1], y = (i > 1), qli1 = qual[i+1];
171 uint8_t qyi1 = query[i+1];
172 x = i - bw; beg = beg > x? beg : x;
173 x = i + bw; end = end < x? end : x;
174 for (k = end; k >= beg; --k) {
175 int u, v11, v01, v10;
177 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);
178 e = (k >= l_ref? 0 : (ref[k+1] > 3 || qyi1 > 3)? 1. : ref[k+1] == qyi1? 1. - qli1 : qli1 * EM) * bi1[v11];
179 bi[u+0] = e * m[0] + EI * m[1] * bi1[v10+1] + m[2] * bi[v01+2]; // bi1[v11] has been foled into e.
180 bi[u+1] = e * m[3] + EI * m[4] * bi1[v10+1];
181 bi[u+2] = (e * m[6] + m[8] * bi[v01+2]) * y;
182 // fprintf(stderr, "B (%d,%d;%d): %lg,%lg,%lg\n", i, k, u, bi[u], bi[u+1], bi[u+2]); // DEBUG
185 set_u(_beg, bw, i, beg); set_u(_end, bw, i, end); _end += 2;
186 for (k = _beg, y = 1./s[i]; k <= _end; ++k) bi[k] *= y;
189 int beg = 1, end = l_ref < bw + 1? l_ref : bw + 1;
191 for (k = end; k >= beg; --k) {
193 double e = (ref[k] > 3 || query[1] > 3)? 1. : ref[k] == query[1]? 1. - qual[1] : qual[1] * EM;
195 if (u < 3 || u >= bw2*3+3) continue;
196 sum += e * b[1][u+0] * bM + EI * b[1][u+1] * bI;
199 pb = b[0][k] = sum / s[0]; // if everything works as is expected, pb == 1.0
201 is_diff = fabs(pb - 1.) > 1e-7? 1 : 0;
203 for (i = 1; i <= l_query; ++i) {
204 double sum = 0., *fi = f[i], *bi = b[i], max = 0.;
205 int beg = 1, end = l_ref, x, max_k = -1;
206 x = i - bw; beg = beg > x? beg : x;
207 x = i + bw; end = end < x? end : x;
208 for (k = beg; k <= end; ++k) {
212 z = fi[u+0] * bi[u+0]; if (z > max) max = z, max_k = (k-1)<<2 | 0; sum += z;
213 z = fi[u+1] * bi[u+1]; if (z > max) max = z, max_k = (k-1)<<2 | 1; sum += z;
215 max /= sum; sum *= s[i]; // if everything works as is expected, sum == 1.0
216 if (state) state[i-1] = max_k;
217 if (q) k = (int)(-4.343 * log(1. - max) + .499), q[i-1] = k > 100? 99 : k;
219 fprintf(stderr, "(%.10lg,%.10lg) (%d,%d:%c,%c:%d) %lg\n", pb, sum, i-1, max_k>>2,
220 "ACGT"[query[i]], "ACGT"[ref[(max_k>>2)+1]], max_k&3, max); // DEBUG
224 /* // compute the posterior of a gap, but I do not know how to use it...
227 a = calloc(3 * (l_ref + 2), sizeof(double));
228 for (i = 1; i < l_query; ++i) {
229 double sum = 0., *fi = f[i], *bi1 = b[i+1], qli1 = qual[i+1];
230 int beg = 1, end = l_ref, x;
231 uint8_t qyi1 = query[i+1];
232 x = i - bw; beg = beg > x? beg : x;
233 x = i + bw; end = end < x? end : x;
234 for (k = beg; k <= end; ++k) {
235 double *ak = a + 3 * k;
236 int u, v11, v01, v10;
238 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);
239 e = k >= l_ref? 0 : (ref[k+1] > 3 || qyi1 > 3)? 1. : ref[k+1] == qyi1? 1. - qli1 : qli1 * EM;
240 ak[0] += fi[u] * bi1[v11] * m[0] * e;
241 ak[1] += fi[u] * bi1[v10+1] * m[1] * EI;
242 ak[2] += fi[u] * bi1[v01+2] * m[2];
245 for (k = 1; k < l_ref; ++k) {
246 double sum = 0., *ak = a + 3 * k;
247 sum += 1. / (ak[0] + ak[1] + ak[2]);
248 ak[0] *= sum; ak[1] *= sum; ak[2] *= sum;
249 fprintf(stderr, "%d: %lf, %lf, %lf\n", k, ak[0], ak[1], ak[2]);
255 for (i = 0; i <= l_query; ++i) {
256 free(f[i]); free(b[i]);
258 free(f); free(b); free(s); free(_qual);
264 int main(int argc, char *argv[])
266 uint8_t conv[256], *iqual, *ref, *query;
267 int c, l_ref, l_query, i, q = 30, b = 10;
268 while ((c = getopt(argc, argv, "b:q:")) >= 0) {
270 case 'b': b = atoi(optarg); break;
271 case 'q': q = atoi(optarg); break;
274 if (optind + 2 > argc) {
275 fprintf(stderr, "Usage: %s [-q %d] [-b %d] <ref> <query>\n", argv[0], q, b); // example: acttc attc
278 memset(conv, 4, 256);
279 conv['a'] = conv['A'] = 0; conv['c'] = conv['C'] = 1;
280 conv['g'] = conv['G'] = 2; conv['t'] = conv['T'] = 3;
281 ref = (uint8_t*)argv[optind]; query = (uint8_t*)argv[optind+1];
282 l_ref = strlen((char*)ref); l_query = strlen((char*)query);
283 for (i = 0; i < l_ref; ++i) ref[i] = conv[ref[i]];
284 for (i = 0; i < l_query; ++i) query[i] = conv[query[i]];
285 iqual = malloc(l_query);
286 memset(iqual, q, l_query);
288 kpa_glocal(ref, l_ref, query, l_query, iqual, &kpa_par_def, 0, 0);