From b7f4d8a0f411b1b0e13fbed6fdf884ab9807af2d Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sun, 2 May 2021 18:55:37 -0400 Subject: [PATCH] removed the old minimap2 chaining --- Makefile | 4 +- chain.c | 164 ------------------------------------------------------- map.c | 2 - 3 files changed, 1 insertion(+), 169 deletions(-) delete mode 100644 chain.c diff --git a/Makefile b/Makefile index d540940..4118616 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ CFLAGS= -g -Wall -O2 -Wc++-compat #-Wextra CPPFLAGS= -DHAVE_KALLOC INCLUDES= OBJS= kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o \ - chain.o lchain.o align.o hit.o seed.o map.o format.o pe.o esterr.o splitidx.o \ + lchain.o align.o hit.o seed.o map.o format.o pe.o esterr.o splitidx.o \ ksw2_ll_sse.o PROG= minimap2 PROG_EXTRA= sdust minimap2-lite @@ -107,7 +107,6 @@ depend: align.o: minimap.h mmpriv.h bseq.h kseq.h ksw2.h kalloc.h bseq.o: bseq.h kvec.h kalloc.h kseq.h -chain.o: minimap.h mmpriv.h bseq.h kseq.h kalloc.h esterr.o: mmpriv.h minimap.h bseq.h kseq.h example.o: minimap.h kseq.h format.o: kalloc.h mmpriv.h minimap.h bseq.h kseq.h @@ -129,6 +128,5 @@ options.o: mmpriv.h minimap.h bseq.h kseq.h pe.o: mmpriv.h minimap.h bseq.h kseq.h kvec.h kalloc.h ksort.h sdust.o: kalloc.h kdq.h kvec.h sdust.h seed.o: mmpriv.h minimap.h bseq.h kseq.h kalloc.h ksort.h -self-chain.o: minimap.h kseq.h sketch.o: kvec.h kalloc.h mmpriv.h minimap.h bseq.h kseq.h splitidx.o: mmpriv.h minimap.h bseq.h kseq.h diff --git a/chain.c b/chain.c deleted file mode 100644 index a2f7ac5..0000000 --- a/chain.c +++ /dev/null @@ -1,164 +0,0 @@ -#include -#include -#include -#include "minimap.h" -#include "mmpriv.h" -#include "kalloc.h" - -static const char LogTable256[256] = { -#define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n - -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, - LT(4), LT(5), LT(5), LT(6), LT(6), LT(6), LT(6), - LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7) -}; - -static inline int ilog2_32(uint32_t v) -{ - uint32_t t, tt; - if ((tt = v>>16)) return (t = tt>>8) ? 24 + LogTable256[t] : 16 + LogTable256[tt]; - return (t = v>>8) ? 8 + LogTable256[t] : LogTable256[v]; -} - -mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float gap_scale, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km) -{ // TODO: make sure this works when n has more than 32 bits - int32_t k, *f, *p, *t, *v, n_u, n_v; - int64_t i, j, st = 0; - uint64_t *u, *u2, sum_qspan = 0; - float avg_qspan; - mm128_t *b, *w; - - if (_u) *_u = 0, *n_u_ = 0; - if (n == 0 || a == 0) { - kfree(km, a); - return 0; - } - f = (int32_t*)kmalloc(km, n * 4); - p = (int32_t*)kmalloc(km, n * 4); - t = (int32_t*)kmalloc(km, n * 4); - v = (int32_t*)kmalloc(km, n * 4); - memset(t, 0, n * 4); - - for (i = 0; i < n; ++i) sum_qspan += a[i].y>>32&0xff; - avg_qspan = (float)sum_qspan / n; - - // fill the score and backtrack arrays - for (i = 0; i < n; ++i) { - uint64_t ri = a[i].x; - int64_t max_j = -1; - int32_t qi = (int32_t)a[i].y, q_span = a[i].y>>32&0xff; // NB: only 8 bits of span is used!!! - int32_t max_f = q_span, n_skip = 0, min_d; - int32_t sidi = (a[i].y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT; - while (st < i && ri > a[st].x + max_dist_x) ++st; - if (i - st > max_iter) st = i - max_iter; - for (j = i - 1; j >= st; --j) { - int64_t dr = ri - a[j].x; - int32_t dq = qi - (int32_t)a[j].y, dd, sc, log_dd, gap_cost; - int32_t sidj = (a[j].y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT; - if ((sidi == sidj && dr == 0) || dq <= 0) continue; // don't skip if an anchor is used by multiple segments; see below - if ((sidi == sidj && dq > max_dist_y) || dq > max_dist_x) continue; - dd = dr > dq? dr - dq : dq - dr; - if (sidi == sidj && dd > bw) continue; - if (n_segs > 1 && !is_cdna && sidi == sidj && dr > max_dist_y) continue; - min_d = dq < dr? dq : dr; - sc = min_d > q_span? q_span : dq < dr? dq : dr; - log_dd = dd? ilog2_32(dd) : 0; - gap_cost = 0; - if (is_cdna || sidi != sidj) { - int c_log, c_lin; - c_lin = (int)(dd * .01 * avg_qspan); - c_log = log_dd; - if (sidi != sidj && dr == 0) ++sc; // possibly due to overlapping paired ends; give a minor bonus - else if (dr > dq || sidi != sidj) gap_cost = c_lin < c_log? c_lin : c_log; - else gap_cost = c_lin + (c_log>>1); - } else gap_cost = (int)(dd * .01 * avg_qspan) + (log_dd>>1); - sc -= (int)((double)gap_cost * gap_scale + .499); - sc += f[j]; - if (sc > max_f) { - max_f = sc, max_j = j; - if (n_skip > 0) --n_skip; - } else if (t[j] == i) { - if (++n_skip > max_skip) - break; - } - if (p[j] >= 0) t[p[j]] = i; - } - f[i] = max_f, p[i] = max_j; - v[i] = max_j >= 0 && v[max_j] > max_f? v[max_j] : max_f; // v[] keeps the peak score up to i; f[] is the score ending at i, not always the peak - } - - // find the ending positions of chains - memset(t, 0, n * 4); - for (i = 0; i < n; ++i) - if (p[i] >= 0) t[p[i]] = 1; - for (i = n_u = 0; i < n; ++i) - if (t[i] == 0 && v[i] >= min_sc) - ++n_u; - if (n_u == 0) { - kfree(km, a); kfree(km, f); kfree(km, p); kfree(km, t); kfree(km, v); - return 0; - } - u = (uint64_t*)kmalloc(km, n_u * 8); - for (i = n_u = 0; i < n; ++i) { - if (t[i] == 0 && v[i] >= min_sc) { - j = i; - while (j >= 0 && f[j] < v[j]) j = p[j]; // find the peak that maximizes f[] - if (j < 0) j = i; // TODO: this should really be assert(j>=0) - u[n_u++] = (uint64_t)f[j] << 32 | j; - } - } - radix_sort_64(u, u + n_u); - for (i = 0; i < n_u>>1; ++i) { // reverse, s.t. the highest scoring chain is the first - uint64_t t = u[i]; - u[i] = u[n_u - i - 1], u[n_u - i - 1] = t; - } - - // backtrack - memset(t, 0, n * 4); - for (i = n_v = k = 0; i < n_u; ++i) { // starting from the highest score - int32_t n_v0 = n_v, k0 = k; - j = (int32_t)u[i]; - do { - v[n_v++] = j; - t[j] = 1; - j = p[j]; - } while (j >= 0 && t[j] == 0); - if (j < 0) { - if (n_v - n_v0 >= min_cnt) u[k++] = u[i]>>32<<32 | (n_v - n_v0); - } else if ((int32_t)(u[i]>>32) - f[j] >= min_sc) { - if (n_v - n_v0 >= min_cnt) u[k++] = ((u[i]>>32) - f[j]) << 32 | (n_v - n_v0); - } - if (k0 == k) n_v = n_v0; // no new chain added, reset - } - *n_u_ = n_u = k, *_u = u; // NB: note that u[] may not be sorted by score here - - // free temporary arrays - kfree(km, f); kfree(km, p); kfree(km, t); - - // write the result to b[] - b = (mm128_t*)kmalloc(km, n_v * sizeof(mm128_t)); - for (i = 0, k = 0; i < n_u; ++i) { - int32_t k0 = k, ni = (int32_t)u[i]; - for (j = 0; j < ni; ++j) - b[k] = a[v[k0 + (ni - j - 1)]], ++k; - } - kfree(km, v); - - // sort u[] and a[] by a[].x, such that adjacent chains may be joined (required by mm_join_long) - w = (mm128_t*)kmalloc(km, n_u * sizeof(mm128_t)); - for (i = k = 0; i < n_u; ++i) { - w[i].x = b[k].x, w[i].y = (uint64_t)k<<32|i; - k += (int32_t)u[i]; - } - radix_sort_128x(w, w + n_u); - u2 = (uint64_t*)kmalloc(km, n_u * 8); - for (i = k = 0; i < n_u; ++i) { - int32_t j = (int32_t)w[i].y, n = (int32_t)u[j]; - u2[i] = u[j]; - memcpy(&a[k], &b[w[i].y>>32], n * sizeof(mm128_t)); - k += n; - } - if (n_u) memcpy(u, u2, n_u * 8); - if (k) memcpy(b, a, k * sizeof(mm128_t)); // write _a_ to _b_ and deallocate _a_ because _a_ is oversized, sometimes a lot - kfree(km, a); kfree(km, w); kfree(km, u2); - return b; -} diff --git a/map.c b/map.c index cf730d9..48d1f09 100644 --- a/map.c +++ b/map.c @@ -271,7 +271,6 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** if (max_chain_gap_ref < opt->max_gap) max_chain_gap_ref = opt->max_gap; } else max_chain_gap_ref = opt->max_gap; -// a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score, opt->chain_gap_scale, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km); a = mg_lchain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score, opt->chain_gap_scale * 0.2f, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km); @@ -295,7 +294,6 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** kfree(b->km, mini_pos); if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos); else a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos); -// a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score, opt->chain_gap_scale, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km); a = mg_lchain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score, opt->chain_gap_scale * 0.2f, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km); }