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