diff --git a/hit.c b/hit.c index 25cf538..b6585d1 100644 --- a/hit.c +++ b/hit.c @@ -79,6 +79,35 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r) // and compu kfree(km, w); } +void mm_update_parent(void *km, float mask_level, int n, mm_reg1_t *r) // due to changes to r.{qs,qe} after DP extension +{ + int i, j, k, *w, n_pri = 0; + if (n <= 0) return; + for (i = 0; i < n; ++i) + if (r[i].id == r[i].parent) ++n_pri; + if (n_pri <= 1) return; + w = (int*)kmalloc(km, n_pri * sizeof(int)); + for (i = j = 0; i < n; ++i) // find the first primary + if (r[i].id == r[i].parent) break; + for (w[0] = i, i = i + 1, k = 1; i < n; ++i) { + int si = r[i].qs, ei = r[i].qe; + if (r[i].id != r[i].parent) continue; // only check primary + for (j = 0; j < k; ++j) { + int sj = r[w[j]].qs, ej = r[w[j]].qe; + int min = ej - sj < ei - si? ej - sj : ei - si; + int ol = si < sj? (ei < sj? 0 : ei < ej? ei - sj : ej - sj) : (ej < si? 0 : ej < ei? ej - si : ei - si); + if (ol > mask_level * min) { + r[i].parent = r[w[j]].parent; + if (r[w[j]].subsc < r[i].score) + r[w[j]].subsc = r[i].score; + break; + } + } + if (j == k) w[k++] = i; + } + kfree(km, w); +} + void mm_sync_regs(void *km, int n_regs, mm_reg1_t *regs) // keep mm_reg1_t::{id,parent} in sync; also reset id { int *tmp, i, max_id = -1, n_tmp, n_pri; @@ -109,7 +138,6 @@ void mm_select_sub(void *km, float mask_level, float pri_ratio, int *n_, mm_reg1 { if (pri_ratio > 0.0f && *n_ > 0) { int i, k, n = *n_; - mm_set_parent(km, mask_level, n, r); for (i = k = 0; i < n; ++i) if (r[i].parent == i || r[i].score >= r[r[i].parent].score * pri_ratio) r[k++] = r[i]; diff --git a/main.c b/main.c index d24c952..dfe4ddc 100644 --- a/main.c +++ b/main.c @@ -10,7 +10,7 @@ #include "minimap.h" #include "mmpriv.h" -#define MM_VERSION "2.0-r121-pre" +#define MM_VERSION "2.0-r122-pre" void liftrlimit() { diff --git a/map.c b/map.c index 1ec5fec..1f4f9c9 100644 --- a/map.c +++ b/map.c @@ -238,11 +238,17 @@ mm_reg1_t *mm_map_frag(const mm_mapopt_t *opt, const mm_idx_t *mi, mm_tbuf_t *b, 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_set_parent(b->km, opt->mask_level, *n_regs, regs); mm_select_sub(b->km, opt->mask_level, opt->pri_ratio, n_regs, regs); mm_join_long(b->km, opt, qlen, *n_regs, regs, a); // TODO: this can be applied to all-vs-all in principle } - if (opt->flag & MM_F_CIGAR) + if (opt->flag & MM_F_CIGAR) { regs = mm_align_skeleton(b->km, opt, mi, qlen, seq, n_regs, regs, a); + if (!(opt->flag & MM_F_AVA)) { + mm_update_parent(b->km, opt->mask_level, *n_regs, regs); + mm_select_sub(b->km, opt->mask_level, opt->pri_ratio, n_regs, regs); + } + } mm_set_mapq(*n_regs, regs); // free diff --git a/mmpriv.h b/mmpriv.h index 65f82a1..bcbb3c7 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -39,6 +39,8 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m 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_set_parent(void *km, float mask_level, int n, mm_reg1_t *r); +void mm_update_parent(void *km, float mask_level, int n, mm_reg1_t *r); void mm_select_sub(void *km, float mask_level, float pri_ratio, int *n_, mm_reg1_t *r); void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int n_regs, mm_reg1_t *regs, mm128_t *a); void mm_set_mapq(int n_regs, mm_reg1_t *regs);