diff --git a/cli.c b/cli.c index 83aa984..c4806ca 100644 --- a/cli.c +++ b/cli.c @@ -1,7 +1,10 @@ #include #include #include +#include #include "ksw2.h" +#include "kseq.h" +KSEQ_INIT(gzFile, gzread) unsigned char seq_nt4_table[256] = { 0, 1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, @@ -36,17 +39,35 @@ static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b) mat[(m - 1) * m + j] = 0; } +static void global_aln(const char *algo, void *km, char *qseq, char *tseq, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez) +{ + int i, qlen, tlen; + qlen = strlen(qseq); + tlen = strlen(tseq); + for (i = 0; i < qlen; ++i) + qseq[i] = seq_nt4_table[(uint8_t)qseq[i]]; + for (i = 0; i < tlen; ++i) + tseq[i] = seq_nt4_table[(uint8_t)tseq[i]]; + if (w < 0) w = qlen > tlen? qlen : tlen; + if (strcmp(algo, "gg") == 0) ez->score = ksw_gg(0, qlen, (uint8_t*)qseq, tlen, (uint8_t*)tseq, 5, mat, q, e, w, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + else if (strcmp(algo, "gg2") == 0) ez->score = ksw_gg2(0, qlen, (uint8_t*)qseq, tlen, (uint8_t*)tseq, 5, mat, q, e, w, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + else if (strcmp(algo, "gg2_sse") == 0) ez->score = ksw_gg2_sse(0, qlen, (uint8_t*)qseq, tlen, (uint8_t*)tseq, 5, mat, q, e, w, &ez->m_cigar, &ez->n_cigar, &ez->cigar); + else if (strcmp(algo, "extz2_sse") == 0) ksw_extz2_sse(0, qlen, (uint8_t*)qseq, tlen, (uint8_t*)tseq, 5, mat, q, e, w, zdrop, flag, ez); + else abort(); +} + int main(int argc, char *argv[]) { int8_t a = 1, b = 1, q = 1, e = 1; - int c, i, qlen, tlen, n_cigar, score; - char *qseq, *tseq, *algo = "gg2"; - uint32_t *cigar; + int c, i, pair = 1, w = -1; + char *algo = "gg2"; int8_t mat[25]; ksw_extz_t ez; + gzFile fp[2]; - while ((c = getopt(argc, argv, "t:")) >= 0) { + while ((c = getopt(argc, argv, "t:w:")) >= 0) { if (c == 't') algo = optarg; + else if (c == 'w') w = atoi(optarg); } if (argc - optind < 2) { fprintf(stderr, "Usage: ksw2-global \n"); @@ -54,25 +75,31 @@ int main(int argc, char *argv[]) } memset(&ez, 0, sizeof(ksw_extz_t)); ksw_gen_simple_mat(5, mat, a, -b); - tseq = argv[optind], qseq = argv[optind+1]; - tlen = strlen(tseq); - qlen = strlen(qseq); - for (i = 0; i < qlen; ++i) - qseq[i] = seq_nt4_table[(uint8_t)qseq[i]]; - for (i = 0; i < tlen; ++i) - tseq[i] = seq_nt4_table[(uint8_t)tseq[i]]; - if (strcmp(algo, "gg") == 0) score = ksw_gg(0, qlen, (uint8_t*)qseq, tlen, (uint8_t*)tseq, 5, mat, q, e, qlen > tlen? qlen : tlen, &n_cigar, &cigar); - else if (strcmp(algo, "gg2") == 0) score = ksw_gg2(0, qlen, (uint8_t*)qseq, tlen, (uint8_t*)tseq, 5, mat, q, e, qlen > tlen? qlen : tlen, &n_cigar, &cigar); - else if (strcmp(algo, "gg2_sse") == 0) score = ksw_gg2_sse(0, qlen, (uint8_t*)qseq, tlen, (uint8_t*)tseq, 5, mat, q, e, qlen > tlen? qlen : tlen, &n_cigar, &cigar); - else if (strcmp(algo, "extz_sse") == 0) { - ksw_extz_sse(0, qlen, (uint8_t*)qseq, tlen, (uint8_t*)tseq, 5, mat, q, e, qlen > tlen? qlen : tlen, 100, 0, &ez); - n_cigar = ez.n_cigar, cigar = ez.cigar, score = ez.score; - printf("max: %d; (%d,%d)\n", ez.max, ez.max_t, ez.max_q); - printf("mqe: %d; %d\n", ez.mqe, ez.mqe_t); - } else abort(); - printf("%d\t", score); - for (i = 0; i < n_cigar; ++i) - printf("%d%c", cigar[i]>>4, "MID"[cigar[i]&0xf]); - printf("\n"); + fp[0] = gzopen(argv[optind], "r"); + fp[1] = gzopen(argv[optind+1], "r"); + + if (fp[0] == 0 && fp[1] == 0) { + global_aln(algo, 0, argv[optind+1], argv[optind], 5, mat, q, e, w, 100, 0, &ez); + printf("%d\t", ez.score); + for (i = 0; i < ez.n_cigar; ++i) + printf("%d%c", ez.cigar[i]>>4, "MID"[ez.cigar[i]&0xf]); + printf("\n"); + } else if (fp[0] && fp[1]) { + kseq_t *ks[2]; + ks[0] = kseq_init(fp[0]); + ks[1] = kseq_init(fp[1]); + if (pair) { + while (kseq_read(ks[0]) > 0) { + if (kseq_read(ks[1]) <= 0) break; + global_aln(algo, 0, ks[0]->seq.s, ks[1]->seq.s, 5, mat, q, e, w, 100, 0, &ez); + printf("%s\t%s\t%d\t", ks[0]->name.s, ks[1]->name.s, ez.score); + for (i = 0; i < ez.n_cigar; ++i) + printf("%d%c", ez.cigar[i]>>4, "MID"[ez.cigar[i]&0xf]); + printf("\n"); + } + } + kseq_destroy(ks[0]); + kseq_destroy(ks[1]); + } return 0; } diff --git a/kseq.h b/kseq.h new file mode 100644 index 0000000..2f94a64 --- /dev/null +++ b/kseq.h @@ -0,0 +1,242 @@ +/* The MIT License + + Copyright (c) 2008, 2009, 2011 Attractive Chaos + + Permission is hereby granted, free of charge, to any person obtaining + a copy of this software and associated documentation files (the + "Software"), to deal in the Software without restriction, including + without limitation the rights to use, copy, modify, merge, publish, + distribute, sublicense, and/or sell copies of the Software, and to + permit persons to whom the Software is furnished to do so, subject to + the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS + BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN + ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + SOFTWARE. +*/ + +/* Last Modified: 05MAR2012 */ + +#ifndef AC_KSEQ_H +#define AC_KSEQ_H + +#include +#include +#include + +#define KS_SEP_SPACE 0 // isspace(): \t, \n, \v, \f, \r +#define KS_SEP_TAB 1 // isspace() && !' ' +#define KS_SEP_LINE 2 // line separator: "\n" (Unix) or "\r\n" (Windows) +#define KS_SEP_MAX 2 + +#define __KS_TYPE(type_t) \ + typedef struct __kstream_t { \ + unsigned char *buf; \ + int begin, end, is_eof; \ + type_t f; \ + } kstream_t; + +#define ks_err(ks) ((ks)->end == -1) +#define ks_eof(ks) ((ks)->is_eof && (ks)->begin >= (ks)->end) +#define ks_rewind(ks) ((ks)->is_eof = (ks)->begin = (ks)->end = 0) + +#define __KS_BASIC(type_t, __bufsize) \ + static inline kstream_t *ks_init(type_t f) \ + { \ + kstream_t *ks = (kstream_t*)calloc(1, sizeof(kstream_t)); \ + ks->f = f; \ + ks->buf = (unsigned char*)malloc(__bufsize); \ + return ks; \ + } \ + static inline void ks_destroy(kstream_t *ks) \ + { \ + if (ks) { \ + free(ks->buf); \ + free(ks); \ + } \ + } + +#define __KS_GETC(__read, __bufsize) \ + static inline int ks_getc(kstream_t *ks) \ + { \ + if (ks_err(ks)) return -3; \ + if (ks->is_eof && ks->begin >= ks->end) return -1; \ + if (ks->begin >= ks->end) { \ + ks->begin = 0; \ + ks->end = __read(ks->f, ks->buf, __bufsize); \ + if (ks->end == 0) { ks->is_eof = 1; return -1;} \ + if (ks->end == -1) { ks->is_eof = 1; return -3;}\ + } \ + return (int)ks->buf[ks->begin++]; \ + } + +#ifndef KSTRING_T +#define KSTRING_T kstring_t +typedef struct __kstring_t { + size_t l, m; + char *s; +} kstring_t; +#endif + +#ifndef kroundup32 +#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) +#endif + +#define __KS_GETUNTIL(__read, __bufsize) \ + static int ks_getuntil2(kstream_t *ks, int delimiter, kstring_t *str, int *dret, int append) \ + { \ + int gotany = 0; \ + if (dret) *dret = 0; \ + str->l = append? str->l : 0; \ + for (;;) { \ + int i; \ + if (ks_err(ks)) return -3; \ + if (ks->begin >= ks->end) { \ + if (!ks->is_eof) { \ + ks->begin = 0; \ + ks->end = __read(ks->f, ks->buf, __bufsize); \ + if (ks->end == 0) { ks->is_eof = 1; break; } \ + if (ks->end == -1) { ks->is_eof = 1; return -3; } \ + } else break; \ + } \ + if (delimiter == KS_SEP_LINE) { \ + for (i = ks->begin; i < ks->end; ++i) \ + if (ks->buf[i] == '\n') break; \ + } else if (delimiter > KS_SEP_MAX) { \ + for (i = ks->begin; i < ks->end; ++i) \ + if (ks->buf[i] == delimiter) break; \ + } else if (delimiter == KS_SEP_SPACE) { \ + for (i = ks->begin; i < ks->end; ++i) \ + if (isspace(ks->buf[i])) break; \ + } else if (delimiter == KS_SEP_TAB) { \ + for (i = ks->begin; i < ks->end; ++i) \ + if (isspace(ks->buf[i]) && ks->buf[i] != ' ') break; \ + } else i = 0; /* never come to here! */ \ + if (str->m - str->l < (size_t)(i - ks->begin + 1)) { \ + str->m = str->l + (i - ks->begin) + 1; \ + kroundup32(str->m); \ + str->s = (char*)realloc(str->s, str->m); \ + } \ + gotany = 1; \ + memcpy(str->s + str->l, ks->buf + ks->begin, i - ks->begin); \ + str->l = str->l + (i - ks->begin); \ + ks->begin = i + 1; \ + if (i < ks->end) { \ + if (dret) *dret = ks->buf[i]; \ + break; \ + } \ + } \ + if (!gotany && ks_eof(ks)) return -1; \ + if (str->s == 0) { \ + str->m = 1; \ + str->s = (char*)calloc(1, 1); \ + } else if (delimiter == KS_SEP_LINE && str->l > 1 && str->s[str->l-1] == '\r') --str->l; \ + str->s[str->l] = '\0'; \ + return str->l; \ + } \ + static inline int ks_getuntil(kstream_t *ks, int delimiter, kstring_t *str, int *dret) \ + { return ks_getuntil2(ks, delimiter, str, dret, 0); } + +#define KSTREAM_INIT(type_t, __read, __bufsize) \ + __KS_TYPE(type_t) \ + __KS_BASIC(type_t, __bufsize) \ + __KS_GETC(__read, __bufsize) \ + __KS_GETUNTIL(__read, __bufsize) + +#define kseq_rewind(ks) ((ks)->last_char = (ks)->f->is_eof = (ks)->f->begin = (ks)->f->end = 0) + +#define __KSEQ_BASIC(SCOPE, type_t) \ + SCOPE kseq_t *kseq_init(type_t fd) \ + { \ + kseq_t *s = (kseq_t*)calloc(1, sizeof(kseq_t)); \ + s->f = ks_init(fd); \ + return s; \ + } \ + SCOPE void kseq_destroy(kseq_t *ks) \ + { \ + if (!ks) return; \ + free(ks->name.s); free(ks->comment.s); free(ks->seq.s); free(ks->qual.s); \ + ks_destroy(ks->f); \ + free(ks); \ + } + +/* Return value: + >=0 length of the sequence (normal) + -1 end-of-file + -2 truncated quality string + -3 error reading stream + */ +#define __KSEQ_READ(SCOPE) \ + SCOPE int kseq_read(kseq_t *seq) \ + { \ + int c,r; \ + kstream_t *ks = seq->f; \ + if (seq->last_char == 0) { /* then jump to the next header line */ \ + while ((c = ks_getc(ks)) >= 0 && c != '>' && c != '@'); \ + if (c < 0) return c; /* end of file or error*/ \ + seq->last_char = c; \ + } /* else: the first header char has been read in the previous call */ \ + seq->comment.l = seq->seq.l = seq->qual.l = 0; /* reset all members */ \ + if ((r=ks_getuntil(ks, 0, &seq->name, &c)) < 0) return r; /* normal exit: EOF or error */ \ + if (c != '\n') ks_getuntil(ks, KS_SEP_LINE, &seq->comment, 0); /* read FASTA/Q comment */ \ + if (seq->seq.s == 0) { /* we can do this in the loop below, but that is slower */ \ + seq->seq.m = 256; \ + seq->seq.s = (char*)malloc(seq->seq.m); \ + } \ + while ((c = ks_getc(ks)) >= 0 && c != '>' && c != '+' && c != '@') { \ + if (c == '\n') continue; /* skip empty lines */ \ + seq->seq.s[seq->seq.l++] = c; /* this is safe: we always have enough space for 1 char */ \ + ks_getuntil2(ks, KS_SEP_LINE, &seq->seq, 0, 1); /* read the rest of the line */ \ + } \ + if (c == '>' || c == '@') seq->last_char = c; /* the first header char has been read */ \ + if (seq->seq.l + 1 >= seq->seq.m) { /* seq->seq.s[seq->seq.l] below may be out of boundary */ \ + seq->seq.m = seq->seq.l + 2; \ + kroundup32(seq->seq.m); /* rounded to the next closest 2^k */ \ + seq->seq.s = (char*)realloc(seq->seq.s, seq->seq.m); \ + } \ + seq->seq.s[seq->seq.l] = 0; /* null terminated string */ \ + if (c != '+') return seq->seq.l; /* FASTA */ \ + if (seq->qual.m < seq->seq.m) { /* allocate memory for qual in case insufficient */ \ + seq->qual.m = seq->seq.m; \ + seq->qual.s = (char*)realloc(seq->qual.s, seq->qual.m); \ + } \ + while ((c = ks_getc(ks)) >= 0 && c != '\n'); /* skip the rest of '+' line */ \ + if (c == -1) return -2; /* error: no quality string */ \ + while ((c = ks_getuntil2(ks, KS_SEP_LINE, &seq->qual, 0, 1) >= 0 && seq->qual.l < seq->seq.l)); \ + if (c == -3) return -3; /* stream error */ \ + seq->last_char = 0; /* we have not come to the next header line */ \ + if (seq->seq.l != seq->qual.l) return -2; /* error: qual string is of a different length */ \ + return seq->seq.l; \ + } + +#define __KSEQ_TYPE(type_t) \ + typedef struct { \ + kstring_t name, comment, seq, qual; \ + int last_char; \ + kstream_t *f; \ + } kseq_t; + +#define KSEQ_INIT2(SCOPE, type_t, __read) \ + KSTREAM_INIT(type_t, __read, 16384) \ + __KSEQ_TYPE(type_t) \ + __KSEQ_BASIC(SCOPE, type_t) \ + __KSEQ_READ(SCOPE) + +#define KSEQ_INIT(type_t, __read) KSEQ_INIT2(static, type_t, __read) + +#define KSEQ_DECLARE(type_t) \ + __KS_TYPE(type_t) \ + __KSEQ_TYPE(type_t) \ + extern kseq_t *kseq_init(type_t fd); \ + void kseq_destroy(kseq_t *ks); \ + int kseq_read(kseq_t *seq); + +#endif diff --git a/ksw2.h b/ksw2.h index cdcef2b..b3cc8d1 100644 --- a/ksw2.h +++ b/ksw2.h @@ -38,20 +38,21 @@ extern "C" { * @param flag flag (see KSW_EZ_* macros) * @param ez (out) scores and cigar */ -void ksw_extz_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez); +void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez); /** * Global alignment * * (first 10 parameters identical to ksw_extz_sse()) + * @param m_cigar (modified) max CIGAR length; feed 0 if cigar==0 * @param n_cigar (out) number of CIGAR elements * @param cigar (out) BAM-encoded CIGAR; caller need to deallocate with kfree(km, ) * * @return score of the alignment */ -int ksw_gg(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *n_cigar_, uint32_t **cigar_); -int ksw_gg2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *n_cigar_, uint32_t **cigar_); -int ksw_gg2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *n_cigar_, uint32_t **cigar_); +int ksw_gg(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_); +int ksw_gg2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_); +int ksw_gg2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_); #ifdef __cplusplus } diff --git a/ksw2_ext.c b/ksw2_ext.c index 4c1ea15..1cb6424 100644 --- a/ksw2_ext.c +++ b/ksw2_ext.c @@ -20,9 +20,9 @@ int ksw_ext(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *t } // fill the first row - eh[0].h = h0, eh[0].e = h0 - gapoe - gapo > 0? h0 - gapoe - gapo : 0; + eh[0].h = h0, eh[0].e = h0 - gapoe - gapoe > 0? h0 - gapoe - gapoe : 0; for (j = 1; j <= qlen && j <= w; ++j) { - eh[j].h = h0 - (gapo + gape * j), eh[j].e = h0 - (gapoe + gapo + gape * j); + eh[j].h = h0 - (gapoe + gape * j), eh[j].e = h0 - (gapoe + gapoe + gape * j); if (eh[j].e < 0) eh[j].e = 0; if (eh[j].h < 0) { eh[j].h = 0; @@ -49,9 +49,9 @@ int ksw_ext(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *t if (en > max_j0 + w + 1) en = max_j0 + w + 1; // compute the first column if (st == 0) { - h1 = h0 - (gapo + gape * i); + h1 = h0 - (gapoe + gape * i); if (h1 < 0) h1 = 0; - f = h0 - (gapoe + gapo + gape * i); + f = h0 - (gapoe + gapoe + gape * i); if (f < 0) f = 0; } else h1 = 0; for (j = st; j < en; ++j) { diff --git a/ksw2_extz2_sse.c b/ksw2_extz2_sse.c index b94f75f..f2f7534 100644 --- a/ksw2_extz2_sse.c +++ b/ksw2_extz2_sse.c @@ -8,7 +8,7 @@ #include #endif -void ksw_extz_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez) +void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int flag, ksw_extz_t *ez) { #define __dp_code_block1 \ z = _mm_add_epi8(_mm_loadu_si128((__m128i*)&s[t]), qe2_); \ diff --git a/ksw2_gg.c b/ksw2_gg.c index 095b3aa..998e0db 100644 --- a/ksw2_gg.c +++ b/ksw2_gg.c @@ -3,7 +3,7 @@ typedef struct { int32_t h, e; } eh_t; -int ksw_gg(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *n_cigar_, uint32_t **cigar_) +int ksw_gg(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_) { eh_t *eh; int8_t *qp; // query profile @@ -100,8 +100,8 @@ int ksw_gg(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *ta // backtrack score = eh[qlen].h; - if (n_cigar_ && cigar_) { - int n_cigar = 0, m_cigar = 0, which = 0; + if (m_cigar_ && n_cigar_ && cigar_) { + int n_cigar = 0, m_cigar = *m_cigar_, which = 0; uint32_t *cigar = 0, tmp; i = tlen - 1, k = last_en - 1; // (i,k) points to the last cell; FIXME: with a moving band, we need to take care of last deletion/insertion!!! while (i >= 0 && k >= 0) { @@ -117,7 +117,7 @@ int ksw_gg(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *ta if (k >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, k + 1); // first insertion for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR tmp = cigar[i], cigar[i] = cigar[n_cigar-1-i], cigar[n_cigar-1-i] = tmp; - *n_cigar_ = n_cigar, *cigar_ = cigar; + *m_cigar_ = m_cigar, *n_cigar_ = n_cigar, *cigar_ = cigar; } kfree(km, qp); kfree(km, eh); diff --git a/ksw2_gg2.c b/ksw2_gg2.c index 294b342..e15bad4 100644 --- a/ksw2_gg2.c +++ b/ksw2_gg2.c @@ -1,7 +1,7 @@ #include // for debugging only #include "ksw2.h" -int ksw_gg2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int *n_cigar_, uint32_t **cigar_) +int ksw_gg2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_) { int qe = q + e, qe2 = qe + qe, r, t, n_col, *off, score; int8_t *u, *v, *x, *y, *s; @@ -78,7 +78,7 @@ int ksw_gg2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *t } kfree(km, u); kfree(km, v); kfree(km, x); kfree(km, y); kfree(km, s); kfree(km, qr); { // backtrack - int n_cigar = 0, m_cigar = 0, which = 0, i, j, k, l; + int n_cigar = 0, m_cigar = *m_cigar_, which = 0, i, j, k, l; uint32_t *cigar = 0, tmp; i = tlen - 1, j = qlen - 1; while (i >= 0 && j >= 0) { @@ -95,7 +95,7 @@ int ksw_gg2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *t if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, j + 1); // first insertion for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR tmp = cigar[i], cigar[i] = cigar[n_cigar-1-i], cigar[n_cigar-1-i] = tmp; - *n_cigar_ = n_cigar, *cigar_ = cigar; + *m_cigar_ = m_cigar, *n_cigar_ = n_cigar, *cigar_ = cigar; // compute score for (k = 0, score = 0, i = j = 0; k < n_cigar; ++k) { diff --git a/ksw2_gg2_sse.c b/ksw2_gg2_sse.c index 0d1401d..60ec947 100644 --- a/ksw2_gg2_sse.c +++ b/ksw2_gg2_sse.c @@ -8,7 +8,7 @@ #include #endif -int ksw_gg2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int *n_cigar_, uint32_t **cigar_) +int ksw_gg2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int *m_cigar_, int *n_cigar_, uint32_t **cigar_) { int r, t, n_col, *off, score, tlen16; int8_t *u, *v, *x, *y, *s; @@ -111,7 +111,7 @@ int ksw_gg2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_ } kfree(km, mem); kfree(km, qr); { // backtrack - int n_cigar = 0, m_cigar = 0, which = 0, i, j, k, l; + int n_cigar = 0, m_cigar = *m_cigar_, which = 0, i, j, k, l; uint32_t *cigar = 0, tmp; i = tlen - 1, j = qlen - 1; while (i >= 0 && j >= 0) { @@ -128,7 +128,7 @@ int ksw_gg2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_ if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, j + 1); // first insertion for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR tmp = cigar[i], cigar[i] = cigar[n_cigar-1-i], cigar[n_cigar-1-i] = tmp; - *n_cigar_ = n_cigar, *cigar_ = cigar; + *m_cigar_ = m_cigar, *n_cigar_ = n_cigar, *cigar_ = cigar; // compute score for (k = 0, score = 0, i = j = 0; k < n_cigar; ++k) {