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