diff --git a/format.c b/format.c index b292857..c81aac1 100644 --- a/format.c +++ b/format.c @@ -125,7 +125,7 @@ void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const bseq1_t *t, const mm_r } else mm_sprintf_lite(s, "*"); mm_sprintf_lite(s, "\t*\t0\t0\t"); if ((flag & 0x900) == 0) sam_write_sq(s, t->seq, t->l_seq, r->rev, r->rev); - else if (flag & 0x100) mm_sprintf_lite(s, "\t*"); + else if (flag & 0x100) mm_sprintf_lite(s, "*"); else sam_write_sq(s, t->seq + r->qs, r->qe - r->qs, r->rev, r->rev); mm_sprintf_lite(s, "\t*"); // quality write_tags(s, r); diff --git a/hit.c b/hit.c index cfae504..4c84191 100644 --- a/hit.c +++ b/hit.c @@ -3,22 +3,37 @@ #include "mmpriv.h" #include "kalloc.h" -mm_reg1_t *mm_gen_regs(int qlen, int n_u, uint64_t *u, mm128_t *a) // convert chains to hits +mm_reg1_t *mm_gen_regs(void *km, int qlen, int n_u, uint64_t *u, mm128_t *a) // convert chains to hits { + mm128_t *z, tmp; mm_reg1_t *r; int i, k; - r = (mm_reg1_t*)calloc(n_u, sizeof(mm_reg1_t)); + + if (n_u == 0) return 0; + + // sort by score + z = (mm128_t*)kmalloc(km, n_u * 16); for (i = k = 0; i < n_u; ++i) { + z[i].x = u[i] >> 32; + z[i].y = (uint64_t)k << 32 | (int32_t)u[i]; + k += (int32_t)u[i]; + } + radix_sort_128x(z, z + n_u); + for (i = 0; i < n_u>>1; ++i) // reverse, s.t. larger score first + tmp = z[i], z[i] = z[n_u-1-i], z[n_u-1-i] = tmp; + + // populate r[] + r = (mm_reg1_t*)calloc(n_u, sizeof(mm_reg1_t)); + for (i = 0; i < n_u; ++i) { mm_reg1_t *ri = &r[i]; ri->id = i; ri->parent = MM_PARENT_UNSET; - ri->subsc = 0; - ri->score = u[i]>>32; - ri->cnt = (int32_t)u[i]; - ri->as = k; + ri->score = z[i].x; + ri->cnt = (int32_t)z[i].y; + ri->as = z[i].y >> 32; mm_reg_set_coor(ri, qlen, a); - k += ri->cnt; } + kfree(km, z); return r; } diff --git a/main.c b/main.c index 5f95cc9..c0d03a5 100644 --- a/main.c +++ b/main.c @@ -10,7 +10,7 @@ #include "minimap.h" #include "mmpriv.h" -#define MM_VERSION "2.0-r117-pre" +#define MM_VERSION "2.0-r118-pre" void liftrlimit() { diff --git a/map.c b/map.c index df98af4..f6c0ef9 100644 --- a/map.c +++ b/map.c @@ -235,7 +235,7 @@ mm_reg1_t *mm_map_frag(const mm_mapopt_t *opt, const mm_idx_t *mi, mm_tbuf_t *b, kfree(b->km, m); n_u = mm_chain_dp(opt->max_gap, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, n_a, a, &u, b->km); - regs = mm_gen_regs(qlen, n_u, u, a); + regs = mm_gen_regs(b->km, qlen, n_u, u, a); *n_regs = n_u; if (!(opt->flag & MM_F_AVA)) { // don't choose primary mapping(s) for read overlap mm_select_sub(b->km, opt->mask_level, opt->pri_ratio, n_regs, regs); diff --git a/mmpriv.h b/mmpriv.h index ee19d0f..65f82a1 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -36,7 +36,7 @@ void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const bseq1_t *t, const mm_r int mm_chain_dp(int max_dist, int bw, int max_skip, int min_cnt, int min_sc, int64_t n, mm128_t *a, uint64_t **_u, void *km); mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a); -mm_reg1_t *mm_gen_regs(int qlen, int n_u, uint64_t *u, mm128_t *a); +mm_reg1_t *mm_gen_regs(void *km, int qlen, int n_u, uint64_t *u, mm128_t *a); void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a); void mm_sync_regs(void *km, int n_regs, mm_reg1_t *regs); void mm_select_sub(void *km, float mask_level, float pri_ratio, int *n_, mm_reg1_t *r);