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 int aln_sm_qual[] = {
  77           0, -23, -23, -23, 0,
  78         -23,   0, -23, -23, 0,
  79         -23, -23,   0, -23, 0,
  80         -23, -23, -23,   0, 0,
  81           0,   0,   0,   0, 0
  82 };
  83
  84 ka_param_t ka_param_blast = {  5,  2,   5, 2, aln_sm_blast, 5, 50 };
  85 ka_param_t ka_param_aa2aa = { 10,  2,  10, 2, aln_sm_blosum62, 22, 50 };
  86
  87 ka_param2_t ka_param2_qual  = { 37, 11, 37, 11, 37, 11, 0, 0, aln_sm_qual, 5, 50 };
  88
  89 static uint32_t *ka_path2cigar32(const path_t *path, int path_len, int *n_cigar)
  90 {
  91         int i, n;
  92         uint32_t *cigar;
  93         unsigned char last_type;
  94
  95         if (path_len == 0 || path == 0) {
  96                 *n_cigar = 0;
  97                 return 0;
  98         }
  99
 100         last_type = path->ctype;
 101         for (i = n = 1; i < path_len; ++i) {
 102                 if (last_type != path[i].ctype) ++n;
 103                 last_type = path[i].ctype;
 104         }
 105         *n_cigar = n;
 106         cigar = (uint32_t*)calloc(*n_cigar, 4);
 107
 108         cigar[0] = 1u << 4 | path[path_len-1].ctype;
 109         last_type = path[path_len-1].ctype;
 110         for (i = path_len - 2, n = 0; i >= 0; --i) {
 111                 if (path[i].ctype == last_type) cigar[n] += 1u << 4;
 112                 else {
 113                         cigar[++n] = 1u << 4 | path[i].ctype;
 114                         last_type = path[i].ctype;
 115                 }
 116         }
 117
 118         return cigar;
 119 }
 120
 121 /***************************/
 122 /* START OF common_align.c */
 123 /***************************/
 124
 125 #define SET_INF(s) (s).M = (s).I = (s).D = MINOR_INF;
 126
 127 #define set_M(MM, cur, p, sc)                                                   \
 128 {                                                                                                               \
 129         if ((p)->M >= (p)->I) {                                                         \
 130                 if ((p)->M >= (p)->D) {                                                 \
 131                         (MM) = (p)->M + (sc); (cur)->Mt = FROM_M;       \
 132                 } else {                                                                                \
 133                         (MM) = (p)->D + (sc); (cur)->Mt = FROM_D;       \
 134                 }                                                                                               \
 135         } else {                                                                                        \
 136                 if ((p)->I > (p)->D) {                                                  \
 137                         (MM) = (p)->I + (sc); (cur)->Mt = FROM_I;       \
 138                 } else {                                                                                \
 139                         (MM) = (p)->D + (sc); (cur)->Mt = FROM_D;       \
 140                 }                                                                                               \
 141         }                                                                                                       \
 142 }
 143 #define set_I(II, cur, p)                                                               \
 144 {                                                                                                               \
 145         if ((p)->M - gap_open > (p)->I) {                                       \
 146                 (cur)->It = FROM_M;                                                             \
 147                 (II) = (p)->M - gap_open - gap_ext;                             \
 148         } else {                                                                                        \
 149                 (cur)->It = FROM_I;                                                             \
 150                 (II) = (p)->I - gap_ext;                                                \
 151         }                                                                                                       \
 152 }
 153 #define set_end_I(II, cur, p)                                                   \
 154 {                                                                                                               \
 155         if (gap_end_ext >= 0) {                                                         \
 156                 if ((p)->M - gap_end_open > (p)->I) {                   \
 157                         (cur)->It = FROM_M;                                                     \
 158                         (II) = (p)->M - gap_end_open - gap_end_ext;     \
 159                 } else {                                                                                \
 160                         (cur)->It = FROM_I;                                                     \
 161                         (II) = (p)->I - gap_end_ext;                            \
 162                 }                                                                                               \
 163         } else set_I(II, cur, p);                                                       \
 164 }
 165 #define set_D(DD, cur, p)                                                               \
 166 {                                                                                                               \
 167         if ((p)->M - gap_open > (p)->D) {                                       \
 168                 (cur)->Dt = FROM_M;                                                             \
 169                 (DD) = (p)->M - gap_open - gap_ext;                             \
 170         } else {                                                                                        \
 171                 (cur)->Dt = FROM_D;                                                             \
 172                 (DD) = (p)->D - gap_ext;                                                \
 173         }                                                                                                       \
 174 }
 175 #define set_end_D(DD, cur, p)                                                   \
 176 {                                                                                                               \
 177         if (gap_end_ext >= 0) {                                                         \
 178                 if ((p)->M - gap_end_open > (p)->D) {                   \
 179                         (cur)->Dt = FROM_M;                                                     \
 180                         (DD) = (p)->M - gap_end_open - gap_end_ext;     \
 181                 } else {                                                                                \
 182                         (cur)->Dt = FROM_D;                                                     \
 183                         (DD) = (p)->D - gap_end_ext;                            \
 184                 }                                                                                               \
 185         } else set_D(DD, cur, p);                                                       \
 186 }
 187
 188 typedef struct {
 189         uint8_t Mt:3, It:2, Dt:3;
 190 } dpcell_t;
 191
 192 typedef struct {
 193         int M, I, D;
 194 } dpscore_t;
 195
 196 /***************************
 197  * banded global alignment *
 198  ***************************/
 199 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)
 200 {
 201         int i, j;
 202         dpcell_t **dpcell, *q;
 203         dpscore_t *curr, *last, *s;
 204         int b1, b2, tmp_end;
 205         int *mat, end, max = 0;
 206         uint8_t type, ctype;
 207         uint32_t *cigar = 0;
 208
 209         int gap_open, gap_ext, gap_end_open, gap_end_ext, b;
 210         int *score_matrix, N_MATRIX_ROW;
 211
 212         /* initialize some align-related parameters. just for compatibility */
 213         gap_open = ap->gap_open;
 214         gap_ext = ap->gap_ext;
 215         gap_end_open = ap->gap_end_open;
 216         gap_end_ext = ap->gap_end_ext;
 217         b = ap->band_width;
 218         score_matrix = ap->matrix;
 219         N_MATRIX_ROW = ap->row;
 220
 221         if (n_cigar) *n_cigar = 0;
 222         if (len1 == 0 || len2 == 0) return 0;
 223
 224         /* calculate b1 and b2 */
 225         if (len1 > len2) {
 226                 b1 = len1 - len2 + b;
 227                 b2 = b;
 228         } else {
 229                 b1 = b;
 230                 b2 = len2 - len1 + b;
 231         }
 232         if (b1 > len1) b1 = len1;
 233         if (b2 > len2) b2 = len2;
 234         --seq1; --seq2;
 235
 236         /* allocate memory */
 237         end = (b1 + b2 <= len1)? (b1 + b2 + 1) : (len1 + 1);
 238         dpcell = (dpcell_t**)malloc(sizeof(dpcell_t*) * (len2 + 1));
 239         for (j = 0; j <= len2; ++j)
 240                 dpcell[j] = (dpcell_t*)malloc(sizeof(dpcell_t) * end);
 241         for (j = b2 + 1; j <= len2; ++j)
 242                 dpcell[j] -= j - b2;
 243         curr = (dpscore_t*)malloc(sizeof(dpscore_t) * (len1 + 1));
 244         last = (dpscore_t*)malloc(sizeof(dpscore_t) * (len1 + 1));
 245
 246         /* set first row */
 247         SET_INF(*curr); curr->M = 0;
 248         for (i = 1, s = curr + 1; i < b1; ++i, ++s) {
 249                 SET_INF(*s);
 250                 set_end_D(s->D, dpcell[0] + i, s - 1);
 251         }
 252         s = curr; curr = last; last = s;
 253
 254         /* core dynamic programming, part 1 */
 255         tmp_end = (b2 < len2)? b2 : len2 - 1;
 256         for (j = 1; j <= tmp_end; ++j) {
 257                 q = dpcell[j]; s = curr; SET_INF(*s);
 258                 set_end_I(s->I, q, last);
 259                 end = (j + b1 <= len1 + 1)? (j + b1 - 1) : len1;
 260                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 261                 ++s; ++q;
 262                 for (i = 1; i != end; ++i, ++s, ++q) {
 263                         set_M(s->M, q, last + i - 1, mat[seq1[i]]); /* this will change s->M ! */
 264                         set_I(s->I, q, last + i);
 265                         set_D(s->D, q, s - 1);
 266                 }
 267                 set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 268                 set_D(s->D, q, s - 1);
 269                 if (j + b1 - 1 > len1) { /* bug fixed, 040227 */
 270                         set_end_I(s->I, q, last + i);
 271                 } else s->I = MINOR_INF;
 272                 s = curr; curr = last; last = s;
 273         }
 274         /* last row for part 1, use set_end_D() instead of set_D() */
 275         if (j == len2 && b2 != len2 - 1) {
 276                 q = dpcell[j]; s = curr; SET_INF(*s);
 277                 set_end_I(s->I, q, last);
 278                 end = (j + b1 <= len1 + 1)? (j + b1 - 1) : len1;
 279                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 280                 ++s; ++q;
 281                 for (i = 1; i != end; ++i, ++s, ++q) {
 282                         set_M(s->M, q, last + i - 1, mat[seq1[i]]); /* this will change s->M ! */
 283                         set_I(s->I, q, last + i);
 284                         set_end_D(s->D, q, s - 1);
 285                 }
 286                 set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 287                 set_end_D(s->D, q, s - 1);
 288                 if (j + b1 - 1 > len1) { /* bug fixed, 040227 */
 289                         set_end_I(s->I, q, last + i);
 290                 } else s->I = MINOR_INF;
 291                 s = curr; curr = last; last = s;
 292                 ++j;
 293         }
 294
 295         /* core dynamic programming, part 2 */
 296         for (; j <= len2 - b2 + 1; ++j) {
 297                 SET_INF(curr[j - b2]);
 298                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 299                 end = j + b1 - 1;
 300                 for (i = j - b2 + 1, q = dpcell[j] + i, s = curr + i; i != end; ++i, ++s, ++q) {
 301                         set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 302                         set_I(s->I, q, last + i);
 303                         set_D(s->D, q, s - 1);
 304                 }
 305                 set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 306                 set_D(s->D, q, s - 1);
 307                 s->I = MINOR_INF;
 308                 s = curr; curr = last; last = s;
 309         }
 310
 311         /* core dynamic programming, part 3 */
 312         for (; j < len2; ++j) {
 313                 SET_INF(curr[j - b2]);
 314                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 315                 for (i = j - b2 + 1, q = dpcell[j] + i, s = curr + i; i < len1; ++i, ++s, ++q) {
 316                         set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 317                         set_I(s->I, q, last + i);
 318                         set_D(s->D, q, s - 1);
 319                 }
 320                 set_M(s->M, q, last + len1 - 1, mat[seq1[i]]);
 321                 set_end_I(s->I, q, last + i);
 322                 set_D(s->D, q, s - 1);
 323                 s = curr; curr = last; last = s;
 324         }
 325         /* last row */
 326         if (j == len2) {
 327                 SET_INF(curr[j - b2]);
 328                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 329                 for (i = j - b2 + 1, q = dpcell[j] + i, s = curr + i; i < len1; ++i, ++s, ++q) {
 330                         set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 331                         set_I(s->I, q, last + i);
 332                         set_end_D(s->D, q, s - 1);
 333                 }
 334                 set_M(s->M, q, last + len1 - 1, mat[seq1[i]]);
 335                 set_end_I(s->I, q, last + i);
 336                 set_end_D(s->D, q, s - 1);
 337                 s = curr; curr = last; last = s;
 338         }
 339
 340         *_score = last[len1].M;
 341         if (n_cigar) { /* backtrace */
 342                 path_t *p, *path = (path_t*)malloc(sizeof(path_t) * (len1 + len2 + 2));
 343                 i = len1; j = len2;
 344                 q = dpcell[j] + i;
 345                 s = last + len1;
 346                 max = s->M; type = q->Mt; ctype = FROM_M;
 347                 if (s->I > max) { max = s->I; type = q->It; ctype = FROM_I; }
 348                 if (s->D > max) { max = s->D; type = q->Dt; ctype = FROM_D; }
 349
 350                 p = path;
 351                 p->ctype = ctype; p->i = i; p->j = j; /* bug fixed 040408 */
 352                 ++p;
 353                 do {
 354                         switch (ctype) {
 355                         case FROM_M: --i; --j; break;
 356                         case FROM_I: --j; break;
 357                         case FROM_D: --i; break;
 358                         }
 359                         q = dpcell[j] + i;
 360                         ctype = type;
 361                         switch (type) {
 362                         case FROM_M: type = q->Mt; break;
 363                         case FROM_I: type = q->It; break;
 364                         case FROM_D: type = q->Dt; break;
 365                         }
 366                         p->ctype = ctype; p->i = i; p->j = j;
 367                         ++p;
 368                 } while (i || j);
 369                 cigar = ka_path2cigar32(path, p - path - 1, n_cigar);
 370                 free(path);
 371         }
 372
 373         /* free memory */
 374         for (j = b2 + 1; j <= len2; ++j)
 375                 dpcell[j] += j - b2;
 376         for (j = 0; j <= len2; ++j)
 377                 free(dpcell[j]);
 378         free(dpcell);
 379         free(curr); free(last);
 380
 381         return cigar;
 382 }
 383
 384 typedef struct {
 385         int M, I, D;
 386 } score_aux_t;
 387
 388 #define MINUS_INF -0x40000000
 389
 390 // matrix: len2 rows and len1 columns
 391 int ka_global_score(const uint8_t *_seq1, int len1, const uint8_t *_seq2, int len2, const ka_param2_t *ap)
 392 {
 393
 394 #define __score_aux(_p, _q0, _sc, _io, _ie, _do, _de) {                                 \
 395                 int t1, t2;                                                                                                             \
 396                 score_aux_t *_q;                                                                                                \
 397                 _q = _q0;                                                                                                               \
 398                 _p->M = _q->M >= _q->I? _q->M : _q->I;                                                  \
 399                 _p->M = _p->M >= _q->D? _p->M : _q->D;                                                  \
 400                 _p->M += (_sc);                                                                                                 \
 401                 ++_q;      t1 = _q->M - _io - _ie; t2 = _q->I - _ie; _p->I = t1 >= t2? t1 : t2; \
 402                 _q = _p-1; t1 = _q->M - _do - _de; t2 = _q->D - _de; _p->D = t1 >= t2? t1 : t2; \
 403         }
 404
 405         int i, j, bw, scmat_size = ap->row, *scmat = ap->matrix, ret;
 406         const uint8_t *seq1, *seq2;
 407         score_aux_t *curr, *last, *swap;
 408         bw = abs(len1 - len2) + ap->band_width;
 409         i = len1 > len2? len1 : len2;
 410         if (bw > i + 1) bw = i + 1;
 411         seq1 = _seq1 - 1; seq2 = _seq2 - 1;
 412         curr = calloc(len1 + 2, sizeof(score_aux_t));
 413         last = calloc(len1 + 2, sizeof(score_aux_t));
 414         { // the zero-th row
 415                 int x, end = len1;
 416                 score_aux_t *p;
 417                 j = 0;
 418                 x = j + bw; end = len1 < x? len1 : x; // band end
 419                 p = curr;
 420                 p->M = 0; p->I = p->D = MINUS_INF;
 421                 for (i = 1, p = &curr[1]; i <= end; ++i, ++p)
 422                         p->M = p->I = MINUS_INF, p->D = -(ap->edo + ap->ede * i);
 423                 p->M = p->I = p->D = MINUS_INF;
 424                 swap = curr; curr = last; last = swap;
 425         }
 426         for (j = 1; j < len2; ++j) {
 427                 int x, beg = 0, end = len1, *scrow, col_end;
 428                 score_aux_t *p;
 429                 x = j - bw; beg =    0 > x?    0 : x; // band start
 430                 x = j + bw; end = len1 < x? len1 : x; // band end
 431                 if (beg == 0) { // from zero-th column
 432                         p = curr;
 433                         p->M = p->D = MINUS_INF; p->I = -(ap->eio + ap->eie * j);
 434                         ++beg; // then beg = 1
 435                 }
 436                 scrow = scmat + seq2[j] * scmat_size;
 437                 if (end == len1) col_end = 1, --end;
 438                 else col_end = 0;
 439                 for (i = beg, p = &curr[beg]; i <= end; ++i, ++p)
 440                         __score_aux(p, &last[i-1], scrow[(int)seq1[i]], ap->iio, ap->iie, ap->ido, ap->ide);
 441                 if (col_end) {
 442                         __score_aux(p, &last[i-1], scrow[(int)seq1[i]], ap->eio, ap->eie, ap->ido, ap->ide);
 443                         ++p;
 444                 }
 445                 p->M = p->I = p->D = MINUS_INF;
 446 //              for (i = 0; i <= len1; ++i) printf("(%d,%d,%d) ", curr[i].M, curr[i].I, curr[i].D); putchar('\n');
 447                 swap = curr; curr = last; last = swap;
 448         }
 449         { // the last row
 450                 int x, beg = 0, *scrow;
 451                 score_aux_t *p;
 452                 j = len2;
 453                 x = j - bw; beg = 0 > x?    0 : x; // band start
 454                 if (beg == 0) { // from zero-th column
 455                         p = curr;
 456                         p->M = p->D = MINUS_INF; p->I = -(ap->eio + ap->eie * j);
 457                         ++beg; // then beg = 1
 458                 }
 459                 scrow = scmat + seq2[j] * scmat_size;
 460                 for (i = beg, p = &curr[beg]; i < len1; ++i, ++p)
 461                         __score_aux(p, &last[i-1], scrow[(int)seq1[i]], ap->iio, ap->iie, ap->edo, ap->ede);
 462                 __score_aux(p, &last[i-1], scrow[(int)seq1[i]], ap->eio, ap->eie, ap->edo, ap->ede);
 463 //              for (i = 0; i <= len1; ++i) printf("(%d,%d,%d) ", curr[i].M, curr[i].I, curr[i].D); putchar('\n');
 464         }
 465         ret = curr[len1].M >= curr[len1].I? curr[len1].M : curr[len1].I;
 466         ret = ret >= curr[len1].D? ret : curr[len1].D;
 467         free(curr); free(last);
 468         return ret;
 469 }
 470
 471 #ifdef _MAIN
 472 int main(int argc, char *argv[])
 473 {
 474 //      int len1 = 35, len2 = 35;
 475 //      uint8_t *seq1 = (uint8_t*)"\0\0\3\3\2\0\0\0\1\0\2\1\2\1\3\2\3\3\3\0\2\3\2\1\1\3\3\3\2\3\3\1\0\0\1";
 476 //      uint8_t *seq2 = (uint8_t*)"\0\0\3\3\2\0\0\0\1\0\2\1\2\1\3\2\3\3\3\0\2\3\2\1\1\3\3\3\2\3\3\1\0\1\0";
 477         int len1 = 4, len2 = 4;
 478         uint8_t *seq1 = (uint8_t*)"\1\0\0\1";
 479         uint8_t *seq2 = (uint8_t*)"\1\0\1\0";
 480         int sc;
 481 //      ka_global_core(seq1, 2, seq2, 1, &ka_param_qual, &sc, 0);
 482         sc = ka_global_score(seq1, len1, seq2, len2, &ka_param2_qual);
 483         printf("%d\n", sc);
 484         return 0;
 485 }
 486 #endif