From 788e9d1e3dad7c5477d075371af81f45f1ff55b9 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Mon, 4 Feb 2013 15:40:26 -0500 Subject: [PATCH] fixed a couple of leaks; buggy atm --- bwamem.c | 21 +++++++++++++++++++-- fastmap.c | 10 ++++++++-- ksw.c | 4 ++-- ksw.h | 2 +- 4 files changed, 30 insertions(+), 7 deletions(-) diff --git a/bwamem.c b/bwamem.c index f373ef0..032a54e 100644 --- a/bwamem.c +++ b/bwamem.c @@ -10,7 +10,7 @@ void mem_fill_scmat(int a, int b, int8_t mat[25]) { int i, j, k; - for (i = k = 0; i < 5; ++i) { + for (i = k = 0; i < 4; ++i) { for (j = 0; j < 4; ++j) mat[k++] = i == j? a : -b; mat[k++] = 0; // ambiguous base @@ -233,11 +233,28 @@ mem_aln_t mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, s = &c->seeds[0]; if (s->qbeg + s->len != l_query) { // right extension of the first seed int qle, tle, qe, re; + int16_t *qw = 0; qe = s->qbeg + s->len; re = s->rbeg + s->len - rmax[0]; for (j = 0; j < l_query - qe; ++j) putchar("ACGTN"[(int)query[j+qe]]); putchar('\n'); for (j = 0; j < rmax[1] - rmax[0] - re; ++j) putchar("ACGTN"[(int)rseq[j+re]]); putchar('\n'); - a.score = ksw_extend(l_query - qe, query + qe, rmax[1] - rmax[0] - re, rseq + re, 5, opt->mat, opt->q, opt->r, opt->w, a.score, 0, &qle, &tle); + if (c->n > 1) { // generate $qw + int l = rmax[1] - (s->rbeg + s->len); + assert(l >= 0 && l < 1000); + qw = malloc(l * 2); + for (i = 0; i < l; ++i) qw[i] = -1; // no constraint by default + for (i = 1; i < c->n; ++i) { + const mem_seed_t *t = &c->seeds[i]; + for (j = 0; j < t->len; ++j) { + int x = t->rbeg + j - (s->rbeg + s->len), y = t->qbeg + j - (s->qbeg + s->len); + assert(x < l); + if (qw[x] == -1) qw[x] = x > y? x - y : y - x; + else if (qw[x] >= 0) qw[x] = -2; // in a seed overlap, do not set any constraint + } + } + } + a.score = ksw_extend(l_query - qe, query + qe, rmax[1] - rmax[0] - re, rseq + re, 5, opt->mat, opt->q, opt->r, opt->w, a.score, qw, &qle, &tle); a.qe = qe + qle; a.re = rmax[0] + re + tle; + free(qw); } else a.qe = l_query, a.re = s->rbeg + s->len; a.is_all = 1; diff --git a/fastmap.c b/fastmap.c index 4a677c3..f3100c7 100644 --- a/fastmap.c +++ b/fastmap.c @@ -17,7 +17,7 @@ int main_mem(int argc, char *argv[]) bwt_t *bwt; bntseq_t *bns; int i, j, c; - gzFile *fp; + gzFile fp; kseq_t *seq; uint8_t *pac = 0; @@ -66,9 +66,15 @@ int main_mem(int argc, char *argv[]) putchar('\n'); } puts("//"); + for (i = 0; i < chain.n; ++i) free(chain.chains[i].seeds); + free(chain.chains); } - free(opt); + free(pac); free(opt); + bns_destroy(bns); + bwt_destroy(bwt); + kseq_destroy(seq); + gzclose(fp); return 0; } diff --git a/ksw.c b/ksw.c index 763c774..05f597d 100644 --- a/ksw.c +++ b/ksw.c @@ -315,7 +315,7 @@ typedef struct { int32_t h, e; } eh_t; -int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int h0, const int *qw, int *_qle, int *_tle) +int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int h0, const int16_t *qw, int *_qle, int *_tle) { eh_t *eh; int8_t *qp; @@ -350,7 +350,7 @@ int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, h1 = h0 - (gapo + gape * (i + 1)); if (h1 < 0) h1 = 0; // apply the band and the constraint (if provided) - t = (qw && qw[i] < w)? qw[i] : w; // this is the band width at $i + t = (qw && qw[i] >= 0 && qw[i] < w)? qw[i] : w; // this is the band width at $i if (beg < i - t) beg = i - t; if (end > i + t + 1) end = i + t + 1; if (end > qlen) end = qlen; diff --git a/ksw.h b/ksw.h index 3c9b959..220a8d7 100644 --- a/ksw.h +++ b/ksw.h @@ -49,7 +49,7 @@ extern "C" { /** Unified interface for ksw_sse2_8() and ksw_sse2_16() */ int ksw_sse2(ksw_query_t *q, int tlen, const uint8_t *target, ksw_aux_t *a); - int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int h0, const int *qw, int *_qle, int *_tle); + int ksw_extend(int qlen, const uint8_t *query, int tlen, const uint8_t *target, int m, const int8_t *mat, int gapo, int gape, int w, int h0, const int16_t *qw, int *_qle, int *_tle); #ifdef __cplusplus }