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 = { 1e-4, 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    On input, _ref is the reference sequence and _query is the query
 400    sequence. Both are sequences of 0/1/2/3/4 where 4 stands for an
 401    ambiguous residue. iqual is the base quality. c sets the gap open
 402    probability, gap extension probability and band width.
 403
 404    On output, state and q are arrays of length l_query. The higher 30
 405    bits give the reference position the query base is matched to and the
 406    lower two bits can be 0 (an alignment match) or 1 (an
 407    insertion). q[i] gives the phred scaled posterior probability of
 408    state[i] being wrong.
 409  */
 410 int ka_prob_glocal(const uint8_t *_ref, int l_ref, const uint8_t *_query, int l_query, const uint8_t *iqual,
 411                                    const ka_probpar_t *c, int *state, uint8_t *q)
 412 {
 413         double **f, **b, *s, m[9], sI, sM, bI, bM, pb;
 414         float *qual, *_qual;
 415         const uint8_t *ref, *query;
 416         int bw, bw2, i, k, is_diff = 0;
 417
 418         /*** initialization ***/
 419         ref = _ref - 1; query = _query - 1; // change to 1-based coordinate
 420         bw = l_ref > l_query? l_ref : l_query;
 421         if (bw > c->bw) bw = c->bw;
 422         if (bw < abs(l_ref - l_query)) bw = abs(l_ref - l_query);
 423         bw2 = bw * 2 + 1;
 424         // allocate the forward and backward matrices f[][] and b[][] and the scaling array s[]
 425         f = calloc(l_query+1, sizeof(void*));
 426         b = calloc(l_query+1, sizeof(void*));
 427         for (i = 0; i <= l_query; ++i) {
 428                 f[i] = calloc(bw2 * 3 + 6, sizeof(double)); // FIXME: this is over-allocated for very short seqs
 429                 b[i] = calloc(bw2 * 3 + 6, sizeof(double));
 430         }
 431         s = calloc(l_query+2, sizeof(double)); // s[] is the scaling factor to avoid underflow
 432         // initialize qual
 433         _qual = calloc(l_query, sizeof(float));
 434         if (g_qual2prob[0] == 0)
 435                 for (i = 0; i < 256; ++i)
 436                         g_qual2prob[i] = pow(10, -i/10.);
 437         for (i = 0; i < l_query; ++i) _qual[i] = g_qual2prob[iqual? iqual[i] : 30];
 438         qual = _qual - 1;
 439         // initialize transition probability
 440         sM = sI = 1. / (2 * l_query + 2); // the value here seems not to affect results; FIXME: need proof
 441         m[0*3+0] = (1 - c->d - c->d) * (1 - sM); m[0*3+1] = m[0*3+2] = c->d * (1 - sM);
 442         m[1*3+0] = (1 - c->e) * (1 - sI); m[1*3+1] = c->e * (1 - sI); m[1*3+2] = 0.;
 443         m[2*3+0] = 1 - c->e; m[2*3+1] = 0.; m[2*3+2] = c->e;
 444         bM = (1 - c->d) / l_query; bI = c->d / l_query; // (bM+bI)*l_query==1
 445         /*** forward ***/
 446         // f[0]
 447         set_u(k, bw, 0, 0);
 448         f[0][k] = s[0] = 1.;
 449         { // f[1]
 450                 double *fi = f[1], sum;
 451                 int beg = 1, end = l_ref < bw + 1? l_ref : bw + 1, _beg, _end;
 452                 for (k = beg, sum = 0.; k <= end; ++k) {
 453                         int u;
 454                         double e = (ref[k] > 3 || query[1] > 3)? 1. : ref[k] == query[1]? 1. - qual[1] : qual[1] / 3.;
 455                         set_u(u, bw, 1, k);
 456                         fi[u+0] = e * bM; fi[u+1] = .25 * bI;
 457                         sum += fi[u] + fi[u+1];
 458                 }
 459                 // rescale
 460                 s[1] = sum;
 461                 set_u(_beg, bw, 1, beg); set_u(_end, bw, 1, end); _end += 2;
 462                 for (k = _beg; k <= _end; ++k) fi[k] /= sum;
 463         }
 464         // f[2..l_query]
 465         for (i = 2; i <= l_query; ++i) {
 466                 double *fi = f[i], *fi1 = f[i-1], sum, qli = qual[i];
 467                 int beg = 1, end = l_ref, x, _beg, _end;
 468                 uint8_t qyi = query[i];
 469                 x = i - bw; beg = beg > x? beg : x; // band start
 470                 x = i + bw; end = end < x? end : x; // band end
 471                 for (k = beg, sum = 0.; k <= end; ++k) {
 472                         int u, v11, v01, v10;
 473                         double e;
 474                         e = (ref[k] > 3 || qyi > 3)? 1. : ref[k] == qyi? 1. - qli : qli / 3.;
 475                         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);
 476                         fi[u+0] = e * (m[0] * fi1[v11+0] + m[3] * fi1[v11+1] + m[6] * fi1[v11+2]);
 477                         fi[u+1] = .25 * (m[1] * fi1[v10+0] + m[4] * fi1[v10+1]);
 478                         fi[u+2] = m[2] * fi[v01+0] + m[8] * fi[v01+2];
 479                         sum += fi[u] + fi[u+1] + fi[u+2];
 480 //                      fprintf(stderr, "F (%d,%d;%d): %lg,%lg,%lg\n", i, k, u, fi[u], fi[u+1], fi[u+2]); // DEBUG
 481                 }
 482                 // rescale
 483                 s[i] = sum;
 484                 set_u(_beg, bw, i, beg); set_u(_end, bw, i, end); _end += 2;
 485                 for (k = _beg, sum = 1./sum; k <= _end; ++k) fi[k] *= sum;
 486         }
 487         { // f[l_query+1]
 488                 double sum;
 489                 for (k = 1, sum = 0.; k <= l_ref; ++k) {
 490                         int u;
 491                         set_u(u, bw, l_query, k);
 492                         if (u < 3 || u >= bw2*3+3) continue;
 493                     sum += f[l_query][u+0] * sM + f[l_query][u+1] * sI;
 494                 }
 495                 s[l_query+1] = sum; // the last scaling factor
 496         }
 497         /*** backward ***/
 498         // 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)
 499         for (k = 1; k <= l_ref; ++k) {
 500                 int u;
 501                 double *bi = b[l_query];
 502                 set_u(u, bw, l_query, k);
 503                 if (u < 3 || u >= bw2*3+3) continue;
 504                 bi[u+0] = sM / s[l_query] / s[l_query+1]; bi[u+1] = sI / s[l_query] / s[l_query+1];
 505         }
 506         // b[l_query-1..1]
 507         for (i = l_query - 1; i >= 1; --i) {
 508                 int beg = 1, end = l_ref, x, _beg, _end;
 509                 double *bi = b[i], *bi1 = b[i+1], y = (i > 1), qli1 = qual[i+1];
 510                 uint8_t qyi1 = query[i+1];
 511                 x = i - bw; beg = beg > x? beg : x;
 512                 x = i + bw; end = end < x? end : x;
 513                 for (k = end; k >= beg; --k) {
 514                         int u, v11, v01, v10;
 515                         double e;
 516                         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);
 517                         e = (k >= l_ref? 0 : (ref[k+1] > 3 || qyi1 > 3)? 1. : ref[k+1] == qyi1? 1. - qli1 : qli1 / 3.) * bi1[v11];
 518                         bi[u+0] = e * m[0] + .25 * m[1] * bi1[v10+1] + m[2] * bi[v01+2]; // bi1[v11] has been foled into e.
 519                         bi[u+1] = e * m[3] + .25 * m[4] * bi1[v10+1];
 520                         bi[u+2] = (e * m[6] + m[8] * bi[v01+2]) * y;
 521 //                      fprintf(stderr, "B (%d,%d;%d): %lg,%lg,%lg\n", i, k, u, bi[u], bi[u+1], bi[u+2]); // DEBUG
 522                 }
 523                 // rescale
 524                 set_u(_beg, bw, i, beg); set_u(_end, bw, i, end); _end += 2;
 525                 for (k = _beg, y = 1./s[i]; k <= _end; ++k) bi[k] *= y;
 526         }
 527         { // b[0]
 528                 int beg = 1, end = l_ref < bw + 1? l_ref : bw + 1;
 529                 double sum = 0.;
 530                 for (k = end; k >= beg; --k) {
 531                         int u;
 532                         double e = (ref[k] > 3 || query[1] > 3)? 1. : ref[k] == query[1]? 1. - qual[1] : qual[1] / 3.;
 533                         set_u(u, bw, 1, k);
 534                         if (u < 3 || u >= bw2*3+3) continue;
 535                     sum += e * b[1][u+0] * bM + .25 * b[1][u+1] * bI;
 536                 }
 537                 set_u(k, bw, 0, 0);
 538                 pb = b[0][k] = sum / s[0]; // if everything works as is expected, pb == 1.0
 539         }
 540         is_diff = fabs(pb - 1.) > 1e-7? 1 : 0;
 541         /*** MAP ***/
 542         for (i = 1; i <= l_query; ++i) {
 543                 double sum = 0., *fi = f[i], *bi = b[i], max = 0.;
 544                 int beg = 1, end = l_ref, x, max_k = -1;
 545                 x = i - bw; beg = beg > x? beg : x;
 546                 x = i + bw; end = end < x? end : x;
 547                 for (k = beg; k <= end; ++k) {
 548                         int u;
 549                         double z;
 550                         set_u(u, bw, i, k);
 551                         z = fi[u+0] * bi[u+0]; if (z > max) max = z, max_k = (k-1)<<2 | 0; sum += z;
 552                         z = fi[u+1] * bi[u+1]; if (z > max) max = z, max_k = (k-1)<<2 | 1; sum += z;
 553                 }
 554                 max /= sum; sum *= s[i]; // if everything works as is expected, sum == 1.0
 555                 if (state) state[i-1] = max_k;
 556                 if (q) k = (int)(-4.343 * log(1. - max) + .499), q[i-1] = k > 100? 99 : k;
 557 #ifdef _MAIN
 558                 fprintf(stderr, "(%.10lg,%.10lg) (%d,%d:%d)~%lg\n", pb, sum, i-1, max_k>>2, max_k&3, max); // DEBUG
 559 #endif
 560         }
 561         /*** free ***/
 562         for (i = 0; i <= l_query; ++i) {
 563                 free(f[i]); free(b[i]);
 564         }
 565         free(f); free(b); free(s); free(_qual);
 566         return 0;
 567 }
 568
 569 #ifdef _MAIN
 570 int main()
 571 {
 572         int l_ref = 5, l_query = 4;
 573         uint8_t *ref = (uint8_t*)"\0\1\3\3\1";
 574         uint8_t *query = (uint8_t*)"\0\3\3\1";
 575 //      uint8_t *query = (uint8_t*)"\1\3\3\1";
 576         static uint8_t qual[4] = {20, 20, 20, 20};
 577         ka_prob_glocal(ref, l_ref, query, l_query, qual, &ka_probpar_def, 0, 0);
 578         return 0;
 579 }
 580 #endif