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 };
45 kpa_par_t kpa_par_alt = { 0.0001, 0.01, 10 };
48 The topology of the profile HMM:
54 M[0] M[1] -> ... -> M[k-1] -> M[k] -> ... -> M[L] M[L+1]
59 M[0] points to every {M,I}[k] and every {M,I}[k] points M[L+1].
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.
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
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)
75 double **f, **b = 0, *s, m[9], sI, sM, bI, bM, pb;
77 const uint8_t *ref, *query;
78 int bw, bw2, i, k, is_diff = 0, is_backward = 1, Pr;
80 if ( l_ref<=0 || l_query<=0 ) return 0; // FIXME: this may not be an ideal fix, just prevents sefgault
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);
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));
96 s = calloc(l_query+2, sizeof(double)); // s[] is the scaling factor to avoid underflow
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];
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
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) {
119 double e = (ref[k] > 3 || query[1] > 3)? 1. : ref[k] == query[1]? 1. - qual[1] : qual[1] * EM;
121 fi[u+0] = e * bM; fi[u+1] = EI * bI;
122 sum += fi[u] + fi[u+1];
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;
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;
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
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;
154 for (k = 1, sum = 0.; k <= l_ref; ++k) {
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;
160 s[l_query+1] = sum; // the last scaling factor
162 { // compute likelihood
163 double p = 1., Pr1 = 0.;
164 for (i = 0; i <= l_query + 1; ++i) {
166 if (p < 1e-100) Pr1 += -4.343 * log(p), p = 1.;
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);
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) {
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];
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;
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
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;
207 int beg = 1, end = l_ref < bw + 1? l_ref : bw + 1;
209 for (k = end; k >= beg; --k) {
211 double e = (ref[k] > 3 || query[1] > 3)? 1. : ref[k] == query[1]? 1. - qual[1] : qual[1] * EM;
213 if (u < 3 || u >= bw2*3+3) continue;
214 sum += e * b[1][u+0] * bM + EI * b[1][u+1] * bI;
217 pb = b[0][k] = sum / s[0]; // if everything works as is expected, pb == 1.0
219 is_diff = fabs(pb - 1.) > 1e-7? 1 : 0;
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) {
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;
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;
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
242 for (i = 0; i <= l_query; ++i) {
243 free(f[i]); free(b[i]);
245 free(f); free(b); free(s); free(_qual);
251 int main(int argc, char *argv[])
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) {
257 case 'b': b = atoi(optarg); break;
258 case 'q': q = atoi(optarg); break;
261 if (optind + 2 > argc) {
262 fprintf(stderr, "Usage: %s [-q %d] [-b %d] <ref> <query>\n", argv[0], q, b); // example: acttc attc
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);
275 P = kpa_glocal(ref, l_ref, query, l_query, iqual, &kpa_par_alt, 0, 0);
276 fprintf(stderr, "%d\n", P);