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 banded glocal alignment *
 376  *****************************************/
 377
 378 static float g_qual2prob[256];
 379
 380 #define set_u(u, b, i, k) { int x=(i)-(b); x=x>0?x:0; (u)=((k)-x+1)*3; }
 381
 382 ka_probpar_t ka_probpar_def = { 0.001, 0.1, 10 };
 383
 384 /*
 385   The topology of the profile HMM:
 386
 387            /\             /\        /\             /\
 388            I[1]           I[k-1]    I[k]           I[L]
 389             ^   \      \    ^    \   ^   \      \   ^
 390             |    \      \   |     \  |    \      \  |
 391     M[0]   M[1] -> ... -> M[k-1] -> M[k] -> ... -> M[L]   M[L+1]
 392                 \      \/        \/      \/      /
 393                  \     /\        /\      /\     /
 394                        -> D[k-1] -> D[k] ->
 395                           \/        \/
 396
 397    M[0] points to every {M,I}[k] and every {M,I}[k] points M[L+1].
 398  */
 399 int ka_prob_glocal(const uint8_t *_ref, int l_ref, const uint8_t *_query, int l_query, const uint8_t *iqual,
 400                                    const ka_probpar_t *c, int *state, uint8_t *q)
 401 {
 402         double **f, **b, *s, m[9], sI, sM, bI, bM, pb;
 403         float *qual, *_qual;
 404         const uint8_t *ref, *query;
 405         int bw, bw2, i, k, is_diff = 0;
 406
 407         /*** initialization ***/
 408         ref = _ref - 1; query = _query - 1; // change to 1-based coordinate
 409         bw = l_ref > l_query? l_ref : l_query;
 410         if (bw > c->bw) bw = c->bw;
 411         if (bw < abs(l_ref - l_query)) bw = abs(l_ref - l_query);
 412         bw2 = bw * 2 + 1;
 413         // allocate the forward and backward matrices f[][] and b[][] and the scaling array s[]
 414         f = calloc(l_query+1, sizeof(void*));
 415         b = calloc(l_query+1, sizeof(void*));
 416         for (i = 0; i <= l_query; ++i) {
 417                 f[i] = calloc(bw2 * 3 + 6, sizeof(double)); // FIXME: this is over-allocated for very short seqs
 418                 b[i] = calloc(bw2 * 3 + 6, sizeof(double));
 419         }
 420         s = calloc(l_query+2, sizeof(double)); // s[] is the scaling factor to avoid underflow
 421         // initialize qual
 422         _qual = calloc(l_query, sizeof(float));
 423         if (g_qual2prob[0] == 0)
 424                 for (i = 0; i < 256; ++i)
 425                         g_qual2prob[i] = pow(10, -i/10.);
 426         for (i = 0; i < l_query; ++i) _qual[i] = g_qual2prob[iqual? iqual[i] : 30];
 427         qual = _qual - 1;
 428         // initialize transition probability
 429         sM = sI = 1. / (2 * l_query + 2); // the value here seems not to affect results; FIXME: need proof
 430         m[0*3+0] = (1 - c->d - c->d) * (1 - sM); m[0*3+1] = m[0*3+2] = c->d * (1 - sM);
 431         m[1*3+0] = (1 - c->e) * (1 - sI); m[1*3+1] = c->e * (1 - sI); m[1*3+2] = 0.;
 432         m[2*3+0] = 1 - c->e; m[2*3+1] = 0.; m[2*3+2] = c->e;
 433         bM = (1 - c->d) / l_query; bI = c->d / l_query; // (bM+bI)*l_query==1
 434         /*** forward ***/
 435         // f[0]
 436         set_u(k, bw, 0, 0);
 437         f[0][k] = s[0] = 1.;
 438         { // f[1]
 439                 double *fi = f[1], sum;
 440                 int beg = 1, end = l_ref < bw + 1? l_ref : bw + 1, _beg, _end;
 441                 for (k = beg, sum = 0.; k <= end; ++k) {
 442                         int u;
 443                         double e = (ref[k] > 3 || query[1] > 3)? 1. : ref[k] == query[1]? 1. - qual[1] : qual[1] / 3.;
 444                         set_u(u, bw, 1, k);
 445                         fi[u+0] = e * bM; fi[u+1] = .25 * bI;
 446                         sum += fi[u] + fi[u+1];
 447                 }
 448                 // rescale
 449                 s[1] = sum;
 450                 set_u(_beg, bw, 1, beg); set_u(_end, bw, 1, end); _end += 2;
 451                 for (k = _beg; k <= _end; ++k) fi[k] /= sum;
 452         }
 453         // f[2..l_query]
 454         for (i = 2; i <= l_query; ++i) {
 455                 double *fi = f[i], *fi1 = f[i-1], sum, qli = qual[i];
 456                 int beg = 1, end = l_ref, x, _beg, _end;
 457                 uint8_t qyi = query[i];
 458                 x = i - bw; beg = beg > x? beg : x; // band start
 459                 x = i + bw; end = end < x? end : x; // band end
 460                 for (k = beg, sum = 0.; k <= end; ++k) {
 461                         int u, v11, v01, v10;
 462                         double e;
 463                         e = (ref[k] > 3 || qyi > 3)? 1. : ref[k] == qyi? 1. - qli : qli / 3.;
 464                         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);
 465                         fi[u+0] = e * (m[0] * fi1[v11+0] + m[3] * fi1[v11+1] + m[6] * fi1[v11+2]);
 466                         fi[u+1] = .25 * (m[1] * fi1[v10+0] + m[4] * fi1[v10+1]);
 467                         fi[u+2] = m[2] * fi[v01+0] + m[8] * fi[v01+2];
 468                         sum += fi[u] + fi[u+1] + fi[u+2];
 469 //                      fprintf(stderr, "F (%d,%d;%d): %lg,%lg,%lg\n", i, k, u, fi[u], fi[u+1], fi[u+2]); // DEBUG
 470                 }
 471                 // rescale
 472                 s[i] = sum;
 473                 set_u(_beg, bw, i, beg); set_u(_end, bw, i, end); _end += 2;
 474                 for (k = _beg, sum = 1./sum; k <= _end; ++k) fi[k] *= sum;
 475         }
 476         { // f[l_query+1]
 477                 double sum;
 478                 for (k = 1, sum = 0.; k <= l_ref; ++k) {
 479                         int u;
 480                         set_u(u, bw, l_query, k);
 481                         if (u < 3 || u >= bw2*3+3) continue;
 482                     sum += f[l_query][u+0] * sM + f[l_query][u+1] * sI;
 483                 }
 484                 s[l_query+1] = sum; // the last scaling factor
 485         }
 486         /*** backward ***/
 487         // 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)
 488         for (k = 1; k <= l_ref; ++k) {
 489                 int u;
 490                 double *bi = b[l_query];
 491                 set_u(u, bw, l_query, k);
 492                 if (u < 3 || u >= bw2*3+3) continue;
 493                 bi[u+0] = sM / s[l_query] / s[l_query+1]; bi[u+1] = sI / s[l_query] / s[l_query+1];
 494         }
 495         // b[l_query-1..1]
 496         for (i = l_query - 1; i >= 1; --i) {
 497                 int beg = 1, end = l_ref, x, _beg, _end;
 498                 double *bi = b[i], *bi1 = b[i+1], y = (i > 1), qli1 = qual[i+1];
 499                 uint8_t qyi1 = query[i+1];
 500                 x = i - bw; beg = beg > x? beg : x;
 501                 x = i + bw; end = end < x? end : x;
 502                 for (k = end; k >= beg; --k) {
 503                         int u, v11, v01, v10;
 504                         double e;
 505                         e = k >= l_ref? 0 : (ref[k+1] > 3 || qyi1 > 3)? 1. : ref[k+1] == qyi1? 1. - qli1 : qli1 / 3.;
 506                         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);
 507                         bi[u+0] = e * m[0] * bi1[v11+0] + .25 * m[1] * bi1[v10+1] + m[2] * bi[v01+2];
 508                         bi[u+1] = e * m[3] * bi1[v11+0] + .25 * m[4] * bi1[v10+1];
 509                         bi[u+2] = (e * m[6] * bi1[v11+0] + m[8] * bi[v01+2]) * y;
 510 //                      fprintf(stderr, "B (%d,%d;%d): %lg,%lg,%lg\n", i, k, u, bi[u], bi[u+1], bi[u+2]); // DEBUG
 511                 }
 512                 // rescale
 513                 set_u(_beg, bw, i, beg); set_u(_end, bw, i, end); _end += 2;
 514                 for (k = _beg, y = 1./s[i]; k <= _end; ++k) bi[k] *= y;
 515         }
 516         { // b[0]
 517                 int beg = 1, end = l_ref < bw + 1? l_ref : bw + 1;
 518                 double sum = 0.;
 519                 for (k = end; k >= beg; --k) {
 520                         int u;
 521                         double e = (ref[k] > 3 || query[1] > 3)? 1. : ref[k] == query[1]? 1. - qual[1] : qual[1] / 3.;
 522                         set_u(u, bw, 1, k);
 523                         if (u < 3 || u >= bw2*3+3) continue;
 524                     sum += e * b[1][u+0] * bM + .25 * b[1][u+1] * bI;
 525                 }
 526                 set_u(k, bw, 0, 0);
 527                 pb = b[0][k] = sum / s[0]; // if everything works as is expected, pb == 1.0
 528         }
 529         is_diff = fabs(pb - 1.) > 1e-7? 1 : 0;
 530         /*** MAP ***/
 531         for (i = 1; i <= l_query; ++i) {
 532                 double sum = 0., *fi = f[i], *bi = b[i], max = 0.;
 533                 int beg = 1, end = l_ref, x, max_k = -1;
 534                 x = i - bw; beg = beg > x? beg : x;
 535                 x = i + bw; end = end < x? end : x;
 536                 for (k = beg; k <= end; ++k) {
 537                         int u;
 538                         double z;
 539                         set_u(u, bw, i, k);
 540                         z = fi[u+0] * bi[u+0]; if (z > max) max = z, max_k = (k-1)<<2 | 0; sum += z;
 541                         z = fi[u+1] * bi[u+1]; if (z > max) max = z, max_k = (k-1)<<2 | 1; sum += z;
 542                 }
 543                 max /= sum; sum *= s[i]; // if everything works as is expected, sum == 1.0
 544                 if (state) state[i-1] = max_k;
 545                 if (q) k = -4.343 * log(1. - max), q[i-1] = k > 100? 99 : k;
 546 #ifdef _MAIN
 547                 fprintf(stderr, "(%.10lg,%.10lg) (%d,%d:%d)~%lg\n", pb, sum, i-1, max_k>>2, max_k&3, max); // DEBUG
 548 #endif
 549         }
 550         /*** free ***/
 551         for (i = 0; i <= l_query; ++i) {
 552                 free(f[i]); free(b[i]);
 553         }
 554         free(f); free(b); free(s); free(_qual);
 555         return 0;
 556 }
 557
 558 #ifdef _MAIN
 559 int main()
 560 {
 561         int l_ref = 5, l_query = 4;
 562         uint8_t *ref = (uint8_t*)"\0\1\3\3\1";
 563         uint8_t *query = (uint8_t*)"\0\3\3\1";
 564 //      uint8_t *query = (uint8_t*)"\1\3\3\1";
 565         static uint8_t qual[4] = {20, 20, 20, 20};
 566         ka_prob_glocal(ref, l_ref, query, l_query, qual, &ka_probpar_def, 0, 0);
 567         return 0;
 568 }
 569 #endif