diff --git a/bwamem.c b/bwamem.c index b44e54e..68078d0 100644 --- a/bwamem.c +++ b/bwamem.c @@ -345,14 +345,18 @@ int mem_sort_and_dedup(int n, mem_alnreg_t *a) return m; } -int mem_choose_alnreg_se(const mem_opt_t *opt, int n, mem_alnreg_t *a) // IMPORTANT: must run mem_sort_and_dedup() before calling this function +void mem_mark_primary_se(const mem_opt_t *opt, int n, mem_alnreg_t *a) // IMPORTANT: must run mem_sort_and_dedup() before calling this function { // similar to the loop in mem_chain_flt() - int i, j, m, tmp; - if (n <= 1) return n; - for (i = 0; i < n; ++i) a[i].sub = 0; + int i, k, tmp; + kvec_t(int) z; + if (n == 0) return; + kv_init(z); + for (i = 0; i < n; ++i) a[i].sub = a[i].is_primary = 0; tmp = opt->a + opt->b > opt->q + opt->r? opt->a + opt->b : opt->q + opt->r; - for (i = 1, m = 1; i < n; ++i) { - for (j = 0; j < m; ++j) { + kv_push(int, z, 0); + for (i = 1; i < n; ++i) { + for (k = 0; k < z.n; ++k) { + int j = z.a[k]; int b_max = a[j].qb > a[i].qb? a[j].qb : a[i].qb; int e_min = a[j].qe < a[i].qe? a[j].qe : a[i].qe; if (e_min > b_max) { // have overlap @@ -364,9 +368,10 @@ int mem_choose_alnreg_se(const mem_opt_t *opt, int n, mem_alnreg_t *a) // IMPORT } } } - if (j == m) a[m++] = a[i]; + if (k == z.n) kv_push(int, z, i); } - return m; + for (k = 0; k < z.n; ++k) a[z.a[k]].is_primary = 1; + free(z.a); } /************************ @@ -601,10 +606,11 @@ void mem_sam_se(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, b int k; kstring_t str; str.l = str.m = 0; str.s = 0; - a->n = mem_choose_alnreg_se(opt, a->n, a->a); // NOTE: mem_sort_and_dedup() called in worker1() if (a->n > 0) { + mem_mark_primary_se(opt, a->n, a->a); // NOTE: mem_sort_and_dedup() called in worker1() for (k = 0; k < a->n; ++k) { bwahit_t h; + if (!a->a[k].is_primary) continue; mem_alnreg2hit(&a->a[k], &h); h.qual = approx_mapq_se(opt, &a->a[k]); bwa_hit2sam(&str, opt->mat, opt->q, opt->r, opt->w, bns, pac, s, &h, opt->is_hard); diff --git a/bwamem.h b/bwamem.h index 4e2e5ce..d511254 100644 --- a/bwamem.h +++ b/bwamem.h @@ -34,7 +34,7 @@ typedef struct { typedef struct { int64_t rb, re; int score, qb, qe, seedcov, sub, csub; // sub: suboptimal score; csub: suboptimal inside the chain - int sub_n; + int sub_n, is_primary; } mem_alnreg_t; typedef struct {