Starting the work on next package update.
[samtools.git] / 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