kaln.c

   1 /* The MIT License
   2
   3    Copyright (c) 2003-2006, 2008, 2009, by Heng Li <lh3lh3@gmail.com>
   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 "kaln.h"
  32
  33 #define FROM_M 0
  34 #define FROM_I 1
  35 #define FROM_D 2
  36
  37 typedef struct {
  38         int i, j;
  39         unsigned char ctype;
  40 } path_t;
  41
  42 int aln_sm_blosum62[] = {
  43 /*       A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  *  X */
  44          4,-1,-2,-2, 0,-1,-1, 0,-2,-1,-1,-1,-1,-2,-1, 1, 0,-3,-2, 0,-4, 0,
  45         -1, 5, 0,-2,-3, 1, 0,-2, 0,-3,-2, 2,-1,-3,-2,-1,-1,-3,-2,-3,-4,-1,
  46         -2, 0, 6, 1,-3, 0, 0, 0, 1,-3,-3, 0,-2,-3,-2, 1, 0,-4,-2,-3,-4,-1,
  47         -2,-2, 1, 6,-3, 0, 2,-1,-1,-3,-4,-1,-3,-3,-1, 0,-1,-4,-3,-3,-4,-1,
  48          0,-3,-3,-3, 9,-3,-4,-3,-3,-1,-1,-3,-1,-2,-3,-1,-1,-2,-2,-1,-4,-2,
  49         -1, 1, 0, 0,-3, 5, 2,-2, 0,-3,-2, 1, 0,-3,-1, 0,-1,-2,-1,-2,-4,-1,
  50         -1, 0, 0, 2,-4, 2, 5,-2, 0,-3,-3, 1,-2,-3,-1, 0,-1,-3,-2,-2,-4,-1,
  51          0,-2, 0,-1,-3,-2,-2, 6,-2,-4,-4,-2,-3,-3,-2, 0,-2,-2,-3,-3,-4,-1,
  52         -2, 0, 1,-1,-3, 0, 0,-2, 8,-3,-3,-1,-2,-1,-2,-1,-2,-2, 2,-3,-4,-1,
  53         -1,-3,-3,-3,-1,-3,-3,-4,-3, 4, 2,-3, 1, 0,-3,-2,-1,-3,-1, 3,-4,-1,
  54         -1,-2,-3,-4,-1,-2,-3,-4,-3, 2, 4,-2, 2, 0,-3,-2,-1,-2,-1, 1,-4,-1,
  55         -1, 2, 0,-1,-3, 1, 1,-2,-1,-3,-2, 5,-1,-3,-1, 0,-1,-3,-2,-2,-4,-1,
  56         -1,-1,-2,-3,-1, 0,-2,-3,-2, 1, 2,-1, 5, 0,-2,-1,-1,-1,-1, 1,-4,-1,
  57         -2,-3,-3,-3,-2,-3,-3,-3,-1, 0, 0,-3, 0, 6,-4,-2,-2, 1, 3,-1,-4,-1,
  58         -1,-2,-2,-1,-3,-1,-1,-2,-2,-3,-3,-1,-2,-4, 7,-1,-1,-4,-3,-2,-4,-2,
  59          1,-1, 1, 0,-1, 0, 0, 0,-1,-2,-2, 0,-1,-2,-1, 4, 1,-3,-2,-2,-4, 0,
  60          0,-1, 0,-1,-1,-1,-1,-2,-2,-1,-1,-1,-1,-2,-1, 1, 5,-2,-2, 0,-4, 0,
  61         -3,-3,-4,-4,-2,-2,-3,-2,-2,-3,-2,-3,-1, 1,-4,-3,-2,11, 2,-3,-4,-2,
  62         -2,-2,-2,-3,-2,-1,-2,-3, 2,-1,-1,-2,-1, 3,-3,-2,-2, 2, 7,-1,-4,-1,
  63          0,-3,-3,-3,-1,-2,-2,-3,-3, 3, 1,-2, 1,-1,-2,-2, 0,-3,-1, 4,-4,-1,
  64         -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4, 1,-4,
  65          0,-1,-1,-1,-2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-2, 0, 0,-2,-1,-1,-4,-1
  66 };
  67
  68 int aln_sm_blast[] = {
  69         1, -3, -3, -3, -2,
  70         -3, 1, -3, -3, -2,
  71         -3, -3, 1, -3, -2,
  72         -3, -3, -3, 1, -2,
  73         -2, -2, -2, -2, -2
  74 };
  75
  76 ka_param_t ka_param_blast = {  5,  2,   5, 2, aln_sm_blast, 5, 50 };
  77 ka_param_t ka_param_aa2aa = { 10,  2,  10, 2, aln_sm_blosum62, 22, 50 };
  78
  79 static uint32_t *ka_path2cigar32(const path_t *path, int path_len, int *n_cigar)
  80 {
  81         int i, n;
  82         uint32_t *cigar;
  83         unsigned char last_type;
  84
  85         if (path_len == 0 || path == 0) {
  86                 *n_cigar = 0;
  87                 return 0;
  88         }
  89
  90         last_type = path->ctype;
  91         for (i = n = 1; i < path_len; ++i) {
  92                 if (last_type != path[i].ctype) ++n;
  93                 last_type = path[i].ctype;
  94         }
  95         *n_cigar = n;
  96         cigar = (uint32_t*)calloc(*n_cigar, 4);
  97
  98         cigar[0] = 1u << 4 | path[path_len-1].ctype;
  99         last_type = path[path_len-1].ctype;
 100         for (i = path_len - 2, n = 0; i >= 0; --i) {
 101                 if (path[i].ctype == last_type) cigar[n] += 1u << 4;
 102                 else {
 103                         cigar[++n] = 1u << 4 | path[i].ctype;
 104                         last_type = path[i].ctype;
 105                 }
 106         }
 107
 108         return cigar;
 109 }
 110
 111 /***************************/
 112 /* START OF common_align.c */
 113 /***************************/
 114
 115 #define SET_INF(s) (s).M = (s).I = (s).D = MINOR_INF;
 116
 117 #define set_M(MM, cur, p, sc)                                                   \
 118 {                                                                                                               \
 119         if ((p)->M >= (p)->I) {                                                         \
 120                 if ((p)->M >= (p)->D) {                                                 \
 121                         (MM) = (p)->M + (sc); (cur)->Mt = FROM_M;       \
 122                 } else {                                                                                \
 123                         (MM) = (p)->D + (sc); (cur)->Mt = FROM_D;       \
 124                 }                                                                                               \
 125         } else {                                                                                        \
 126                 if ((p)->I > (p)->D) {                                                  \
 127                         (MM) = (p)->I + (sc); (cur)->Mt = FROM_I;       \
 128                 } else {                                                                                \
 129                         (MM) = (p)->D + (sc); (cur)->Mt = FROM_D;       \
 130                 }                                                                                               \
 131         }                                                                                                       \
 132 }
 133 #define set_I(II, cur, p)                                                               \
 134 {                                                                                                               \
 135         if ((p)->M - gap_open > (p)->I) {                                       \
 136                 (cur)->It = FROM_M;                                                             \
 137                 (II) = (p)->M - gap_open - gap_ext;                             \
 138         } else {                                                                                        \
 139                 (cur)->It = FROM_I;                                                             \
 140                 (II) = (p)->I - gap_ext;                                                \
 141         }                                                                                                       \
 142 }
 143 #define set_end_I(II, cur, p)                                                   \
 144 {                                                                                                               \
 145         if (gap_end_ext >= 0) {                                                         \
 146                 if ((p)->M - gap_end_open > (p)->I) {                   \
 147                         (cur)->It = FROM_M;                                                     \
 148                         (II) = (p)->M - gap_end_open - gap_end_ext;     \
 149                 } else {                                                                                \
 150                         (cur)->It = FROM_I;                                                     \
 151                         (II) = (p)->I - gap_end_ext;                            \
 152                 }                                                                                               \
 153         } else set_I(II, cur, p);                                                       \
 154 }
 155 #define set_D(DD, cur, p)                                                               \
 156 {                                                                                                               \
 157         if ((p)->M - gap_open > (p)->D) {                                       \
 158                 (cur)->Dt = FROM_M;                                                             \
 159                 (DD) = (p)->M - gap_open - gap_ext;                             \
 160         } else {                                                                                        \
 161                 (cur)->Dt = FROM_D;                                                             \
 162                 (DD) = (p)->D - gap_ext;                                                \
 163         }                                                                                                       \
 164 }
 165 #define set_end_D(DD, cur, p)                                                   \
 166 {                                                                                                               \
 167         if (gap_end_ext >= 0) {                                                         \
 168                 if ((p)->M - gap_end_open > (p)->D) {                   \
 169                         (cur)->Dt = FROM_M;                                                     \
 170                         (DD) = (p)->M - gap_end_open - gap_end_ext;     \
 171                 } else {                                                                                \
 172                         (cur)->Dt = FROM_D;                                                     \
 173                         (DD) = (p)->D - gap_end_ext;                            \
 174                 }                                                                                               \
 175         } else set_D(DD, cur, p);                                                       \
 176 }
 177
 178 typedef struct {
 179         uint8_t Mt:3, It:2, Dt:2;
 180 } dpcell_t;
 181
 182 typedef struct {
 183         int M, I, D;
 184 } dpscore_t;
 185
 186 /***************************
 187  * banded global alignment *
 188  ***************************/
 189 uint32_t *ka_global_core(uint8_t *seq1, int len1, uint8_t *seq2, int len2, const ka_param_t *ap, int *_score, int *n_cigar)
 190 {
 191         int i, j;
 192         dpcell_t **dpcell, *q;
 193         dpscore_t *curr, *last, *s;
 194         int b1, b2, tmp_end;
 195         int *mat, end, max = 0;
 196         uint8_t type, ctype;
 197         uint32_t *cigar = 0;
 198
 199         int gap_open, gap_ext, gap_end_open, gap_end_ext, b;
 200         int *score_matrix, N_MATRIX_ROW;
 201
 202         /* initialize some align-related parameters. just for compatibility */
 203         gap_open = ap->gap_open;
 204         gap_ext = ap->gap_ext;
 205         gap_end_open = ap->gap_end_open;
 206         gap_end_ext = ap->gap_end_ext;
 207         b = ap->band_width;
 208         score_matrix = ap->matrix;
 209         N_MATRIX_ROW = ap->row;
 210
 211         *n_cigar = 0;
 212         if (len1 == 0 || len2 == 0) return 0;
 213
 214         /* calculate b1 and b2 */
 215         if (len1 > len2) {
 216                 b1 = len1 - len2 + b;
 217                 b2 = b;
 218         } else {
 219                 b1 = b;
 220                 b2 = len2 - len1 + b;
 221         }
 222         if (b1 > len1) b1 = len1;
 223         if (b2 > len2) b2 = len2;
 224         --seq1; --seq2;
 225
 226         /* allocate memory */
 227         end = (b1 + b2 <= len1)? (b1 + b2 + 1) : (len1 + 1);
 228         dpcell = (dpcell_t**)malloc(sizeof(dpcell_t*) * (len2 + 1));
 229         for (j = 0; j <= len2; ++j)
 230                 dpcell[j] = (dpcell_t*)malloc(sizeof(dpcell_t) * end);
 231         for (j = b2 + 1; j <= len2; ++j)
 232                 dpcell[j] -= j - b2;
 233         curr = (dpscore_t*)malloc(sizeof(dpscore_t) * (len1 + 1));
 234         last = (dpscore_t*)malloc(sizeof(dpscore_t) * (len1 + 1));
 235
 236         /* set first row */
 237         SET_INF(*curr); curr->M = 0;
 238         for (i = 1, s = curr + 1; i < b1; ++i, ++s) {
 239                 SET_INF(*s);
 240                 set_end_D(s->D, dpcell[0] + i, s - 1);
 241         }
 242         s = curr; curr = last; last = s;
 243
 244         /* core dynamic programming, part 1 */
 245         tmp_end = (b2 < len2)? b2 : len2 - 1;
 246         for (j = 1; j <= tmp_end; ++j) {
 247                 q = dpcell[j]; s = curr; SET_INF(*s);
 248                 set_end_I(s->I, q, last);
 249                 end = (j + b1 <= len1 + 1)? (j + b1 - 1) : len1;
 250                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 251                 ++s; ++q;
 252                 for (i = 1; i != end; ++i, ++s, ++q) {
 253                         set_M(s->M, q, last + i - 1, mat[seq1[i]]); /* this will change s->M ! */
 254                         set_I(s->I, q, last + i);
 255                         set_D(s->D, q, s - 1);
 256                 }
 257                 set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 258                 set_D(s->D, q, s - 1);
 259                 if (j + b1 - 1 > len1) { /* bug fixed, 040227 */
 260                         set_end_I(s->I, q, last + i);
 261                 } else s->I = MINOR_INF;
 262                 s = curr; curr = last; last = s;
 263         }
 264         /* last row for part 1, use set_end_D() instead of set_D() */
 265         if (j == len2 && b2 != len2 - 1) {
 266                 q = dpcell[j]; s = curr; SET_INF(*s);
 267                 set_end_I(s->I, q, last);
 268                 end = (j + b1 <= len1 + 1)? (j + b1 - 1) : len1;
 269                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 270                 ++s; ++q;
 271                 for (i = 1; i != end; ++i, ++s, ++q) {
 272                         set_M(s->M, q, last + i - 1, mat[seq1[i]]); /* this will change s->M ! */
 273                         set_I(s->I, q, last + i);
 274                         set_end_D(s->D, q, s - 1);
 275                 }
 276                 set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 277                 set_end_D(s->D, q, s - 1);
 278                 if (j + b1 - 1 > len1) { /* bug fixed, 040227 */
 279                         set_end_I(s->I, q, last + i);
 280                 } else s->I = MINOR_INF;
 281                 s = curr; curr = last; last = s;
 282                 ++j;
 283         }
 284
 285         /* core dynamic programming, part 2 */
 286         for (; j <= len2 - b2 + 1; ++j) {
 287                 SET_INF(curr[j - b2]);
 288                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 289                 end = j + b1 - 1;
 290                 for (i = j - b2 + 1, q = dpcell[j] + i, s = curr + i; i != end; ++i, ++s, ++q) {
 291                         set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 292                         set_I(s->I, q, last + i);
 293                         set_D(s->D, q, s - 1);
 294                 }
 295                 set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 296                 set_D(s->D, q, s - 1);
 297                 s->I = MINOR_INF;
 298                 s = curr; curr = last; last = s;
 299         }
 300
 301         /* core dynamic programming, part 3 */
 302         for (; j < len2; ++j) {
 303                 SET_INF(curr[j - b2]);
 304                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 305                 for (i = j - b2 + 1, q = dpcell[j] + i, s = curr + i; i < len1; ++i, ++s, ++q) {
 306                         set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 307                         set_I(s->I, q, last + i);
 308                         set_D(s->D, q, s - 1);
 309                 }
 310                 set_M(s->M, q, last + len1 - 1, mat[seq1[i]]);
 311                 set_end_I(s->I, q, last + i);
 312                 set_D(s->D, q, s - 1);
 313                 s = curr; curr = last; last = s;
 314         }
 315         /* last row */
 316         if (j == len2) {
 317                 SET_INF(curr[j - b2]);
 318                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 319                 for (i = j - b2 + 1, q = dpcell[j] + i, s = curr + i; i < len1; ++i, ++s, ++q) {
 320                         set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 321                         set_I(s->I, q, last + i);
 322                         set_end_D(s->D, q, s - 1);
 323                 }
 324                 set_M(s->M, q, last + len1 - 1, mat[seq1[i]]);
 325                 set_end_I(s->I, q, last + i);
 326                 set_end_D(s->D, q, s - 1);
 327                 s = curr; curr = last; last = s;
 328         }
 329
 330         *_score = last[len1].M;
 331         if (n_cigar) { /* backtrace */
 332                 path_t *p, *path = (path_t*)malloc(sizeof(path_t) * (len1 + len2 + 2));
 333                 i = len1; j = len2;
 334                 q = dpcell[j] + i;
 335                 s = last + len1;
 336                 max = s->M; type = q->Mt; ctype = FROM_M;
 337                 if (s->I > max) { max = s->I; type = q->It; ctype = FROM_I; }
 338                 if (s->D > max) { max = s->D; type = q->Dt; ctype = FROM_D; }
 339
 340                 p = path;
 341                 p->ctype = ctype; p->i = i; p->j = j; /* bug fixed 040408 */
 342                 ++p;
 343                 do {
 344                         switch (ctype) {
 345                         case FROM_M: --i; --j; break;
 346                         case FROM_I: --j; break;
 347                         case FROM_D: --i; break;
 348                         }
 349                         q = dpcell[j] + i;
 350                         ctype = type;
 351                         switch (type) {
 352                         case FROM_M: type = q->Mt; break;
 353                         case FROM_I: type = q->It; break;
 354                         case FROM_D: type = q->Dt; break;
 355                         }
 356                         p->ctype = ctype; p->i = i; p->j = j;
 357                         ++p;
 358                 } while (i || j);
 359                 cigar = ka_path2cigar32(path, p - path - 1, n_cigar);
 360                 free(path);
 361         }
 362
 363         /* free memory */
 364         for (j = b2 + 1; j <= len2; ++j)
 365                 dpcell[j] += j - b2;
 366         for (j = 0; j <= len2; ++j)
 367                 free(dpcell[j]);
 368         free(dpcell);
 369         free(curr); free(last);
 370
 371         return cigar;
 372 }
 373
 374 /***************************
 375  * Probabilistic extension *
 376  ***************************/
 377
 378 #define set_u(u, b, i, k) { int x=(i)-(b); x=x>0?x:0; (u)=((k)-x+1)*3; }
 379
 380 ka_probpar_t ka_probpar_def = { 0.0001, 0.1, 10 };
 381
 382 /*
 383   The profile HMM is:
 384
 385     /\      /\             /\        /\             /\
 386     I[0]    I[1]           I[k-1]    I[k]           I[L]
 387      ^   \   ^   \      \    ^    \   ^   \      \   ^
 388      |    \  |    \      \   |     \  |    \      \  |
 389     M[0] -> M[1] -> ... -> M[k-1] -> M[k] -> ... -> M[L]   M[L+1]
 390          \       \/     \/        \/      \/      /
 391           \      /\     /\        /\      /\     /
 392             D[1] ->     -> D[k-1] -> D[k] ->
 393             \/             \/        \/
 394
 395   Every {M,I}[k], k=0..L connects M[L+1] at the same probability, while
 396   no D[k], k=1..L-1 connects M[L+1]. This means an alignment can end up
 397   with M or I but not D.
 398
 399   Deletions are dumb states which do not emit residues. Frankly, I am
 400   not sure if they are handled properly in the following
 401   implementation. This is a potential concern to be resolved in future.
 402  */
 403 int ka_prob_extend(uint8_t *_ref, int l_ref, uint8_t *_query, int l_query, float *_qual,
 404                                    const ka_probpar_t *c, int *state, uint8_t *q)
 405 {
 406         double **f, **b, m[9], f_sink, sI, sM, pf, pb, pD;
 407         float *qual;
 408         uint8_t *ref, *query;
 409         int bw, bw2, i, k, is_diff = 0;
 410         /*** initialization ***/
 411         ref = _ref - 1;
 412         query = _query - 1;
 413         qual = _qual - 1;
 414         bw = c->bw;
 415         bw2 = c->bw * 2 + 1;
 416         f = calloc(l_query+1, sizeof(void*));
 417         b = calloc(l_query+1, sizeof(void*));
 418         for (i = 0; i <= l_query; ++i) {
 419                 f[i] = calloc(bw2 * 3 + 6, sizeof(double));
 420                 b[i] = calloc(bw2 * 3 + 6, sizeof(double));
 421         }
 422         // initialize m
 423         sM = sI = 1. / (2 * l_query + 2);
 424         m[0*3+0] = (1 - c->d - c->d) * (1 - sM); m[0*3+1] = m[0*3+2] = c->d * (1 - sM);
 425         m[1*3+0] = (1 - c->e) * (1 - sI); m[1*3+1] = c->e * (1 - sI); m[1*3+2] = 0.;
 426         m[2*3+0] = 1 - c->e; m[2*3+1] = 0.; m[2*3+2] = c->e;
 427         /*** forward ***/
 428         // f[0]
 429         set_u(k, bw, 0, 0);
 430         f[0][k] = 1.;
 431         {
 432                 int beg = 1, end = l_ref, x;
 433                 x = i - bw; beg = beg > x? beg : x;
 434                 x = i + bw; end = end < x? end : x;
 435                 for (k = beg; k <= end; ++k) {
 436                         int u, v01;
 437                         set_u(u, bw, 0, k); set_u(v01, bw, 0, k-1);
 438                         f[0][u+2] = m[2] * f[0][v01+0] + m[8] * f[0][v01+2];
 439                 }
 440         }
 441         // f[1..l_query]; core loop
 442         for (i = 1; i <= l_query; ++i) {
 443                 double *fi = f[i], *fi1 = f[i-1];
 444                 int beg = 0, end = l_ref, x;
 445                 x = i - bw; beg = beg > x? beg : x;
 446                 x = i + bw; end = end < x? end : x;
 447                 for (k = beg; k <= end; ++k) {
 448                         int u, v11, v01, v10;
 449                         double e;
 450                         e = (ref[k] > 3 || query[i] > 3)? 1. : ref[k] == query[i]? 1. - qual[i] : qual[i] / 3.;
 451                         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);
 452                         fi[u+0] = e * (m[0] * fi1[v11+0] + m[3] * fi1[v11+1] + m[6] * fi1[v11+2]);
 453                         fi[u+1] = .25 * (m[1] * fi1[v10+0] + m[4] * fi1[v10+1]);
 454                         fi[u+2] = m[2] * fi[v01+0] + m[8] * fi[v01+2];
 455 //                      fprintf(stderr, "F (%d,%d;%d): %lg,%lg,%lg\n", i, k, u, fi[u], fi[u+1], fi[u+2]); // DEBUG
 456                 }
 457         }
 458         // sink
 459         for (k = 0, f_sink = 0.; k <= l_ref; ++k) {
 460                 int u;
 461                 set_u(u, bw, l_query, k);
 462                 if (u < 3 || u >= bw2*3+3) continue;
 463                 f_sink += f[l_query][u+0] * sM + f[l_query][u+1] * sI;
 464         }
 465         pf = f_sink;
 466         /*** backward ***/
 467         // sink
 468         for (k = 0; k <= l_ref; ++k) {
 469                 int u;
 470                 double *bi = b[l_query];
 471                 set_u(u, bw, l_query, k);
 472                 if (u < 3 || u >= bw2*3+3) continue;
 473                 bi[u+0] = sM; bi[u+1] = sI;
 474         }
 475         // b[l_query-1..1]; core loop
 476         for (i = l_query - 1; i >= 0; --i) {
 477                 int beg = 0, end = l_ref, x;
 478                 double *bi = b[i], *bi1 = b[i+1];
 479                 x = i - bw; beg = beg > x? beg : x;
 480                 x = i + bw; end = end < x? end : x;
 481                 for (k = end; k >= beg; --k) {
 482                         int u, v11, v01, v10;
 483                         double e;
 484                         e = k >= l_ref? 0 : (ref[k+1] > 3 || query[i+1] > 3)? 1. : ref[k+1] == query[i+1]? 1. - qual[i+1] : qual[i+1] / 3.;
 485                         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);
 486                         bi[u+0] = e * m[0] * bi1[v11+0] + .25 * m[1] * bi1[v10+1] + m[2] * bi[v01+2];
 487                         bi[u+1] = e * m[3] * bi1[v11+0] + .25 * m[4] * bi1[v10+1];
 488                         bi[u+2] = e * m[6] * bi1[v11+0] + m[8] * bi[v01+2];
 489 //                      fprintf(stderr, "B (%d,%d;%d): %lg,%lg,%lg\n", i, k, u, bi[u], bi[u+1], bi[u+2]); // DEBUG
 490                 }
 491         }
 492         set_u(k, bw, 0, 0);
 493         pb = b[0][k];
 494         pD = (pf + pb) / 2.;
 495         is_diff = fabs(pb/pf - 1.) > 1e-7? 1 : 0;
 496         /*** MAP ***/
 497         for (i = 1; i <= l_query; ++i) {
 498                 double sum = 0., *fi = f[i], *bi = b[i], max = 0.;
 499                 int beg = 0, end = l_ref, x, max_k = -1;
 500                 x = i - bw; beg = beg > x? beg : x;
 501                 x = i + bw; end = end < x? end : x;
 502                 for (k = beg; k <= end; ++k) {
 503                         int u;
 504                         double z;
 505                         set_u(u, bw, i, k);
 506                         z = fi[u+0] * bi[u+0]; if (z > max) max = z, max_k = k<<2 | 0; sum += z;
 507                         z = fi[u+1] * bi[u+1]; if (z > max) max = z, max_k = k<<2 | 1; sum += z;
 508                 }
 509                 max /= sum; sum /= pD; // if everything works as is expected, sum == 1.0
 510                 if (state) state[i-1] = max_k;
 511                 if (q) q[i-1] = -4.343 * log(1. - max);
 512 #ifdef _MAIN
 513                 fprintf(stderr, "[%d],%d,%lg,%d:%d,%lg\n", is_diff, i, sum, max_k>>2, max_k&3, max); // DEBUG
 514 #endif
 515         }
 516         /*** free ***/
 517         for (i = 0; i <= l_query; ++i) {
 518                 free(f[i]); free(b[i]);
 519         }
 520         free(f); free(b);
 521         return 0;
 522 }
 523
 524 #ifdef _MAIN
 525 int main()
 526 {
 527         int l_ref = 5, l_query = 4;
 528         uint8_t *ref = (uint8_t*)"\0\1\3\3\1";
 529 //      uint8_t *query = (uint8_t*)"\0\3\3\1";
 530         uint8_t *query = (uint8_t*)"\1\3\3\1"; // FIXME: the output is not so right given this input!!!
 531         static float qual[4] = {.01, .01, .01, .01};
 532         ka_prob_extend(ref, l_ref, query, l_query, qual, &ka_probpar_def, 0, 0);
 533         return 0;
 534 }
 535 #endif