From 9943e5fdd0a4d9b4b29d7aa4c3d887e930b72fa1 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 20 Sep 2017 14:35:46 -0400 Subject: [PATCH] backup --- hit.c | 64 ++++++++++++++++++++++++++++++++----- map.c | 96 +++++++++++++++++++++++++++++++++----------------------- mmpriv.h | 9 ++++-- 3 files changed, 118 insertions(+), 51 deletions(-) diff --git a/hit.c b/hit.c index f775239..1a7b99a 100644 --- a/hit.c +++ b/hit.c @@ -293,19 +293,67 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs_, mm_r mm_sync_regs(km, *n_regs_, regs); } } -/* -void mm_sep_seg(void *km, int n_segs, int qlen_sum, int n_regs0, mm_reg1_t *regs0, mm128_t *a, mm_segreg_t *seg) -{ // TODO: this implement is inefficient given a large n_segs - int s, i, j; - memset(seg, 0, n_segs * sizeof(mm_segreg_t)); + +mm_seg_t *mm_seg_gen(void *km, int n_segs, const int *qlens, int n_regs0, const mm_reg1_t *regs0, int *n_regs, mm_reg1_t **regs, const mm128_t *a) +{ + int s, i, j, acc_qlen[MM_MAX_SEG+1], qlen_sum = 0; + mm_seg_t *seg; + + assert(n_segs <= MM_MAX_SEG); + for (s = 1, acc_qlen[0] = 0; s < n_segs; ++s) + acc_qlen[s] = acc_qlen[s-1] + qlens[s-1]; + qlen_sum = acc_qlen[n_segs - 1] + qlens[n_segs - 1]; + + seg = (mm_seg_t*)kcalloc(km, n_segs, sizeof(mm_seg_t)); + for (s = 0; s < n_segs; ++s) { + seg[s].u = (uint64_t*)kmalloc(km, n_regs0 * 8); + for (i = 0; i < n_regs0; ++i) + seg[s].u[i] = (uint64_t)regs0[i].score << 32; + } for (i = 0; i < n_regs0; ++i) { - mm_reg1_t *r = ®s0[i]; + const mm_reg1_t *r = ®s0[i]; for (j = 0; j < r->cnt; ++j) { - ++seg[(a[r->as + j].y&MM_SEED_SEG_MASK)>>MM_SEED_SEG_SHIFT].n_a; + int sid = (a[r->as + j].y&MM_SEED_SEG_MASK)>>MM_SEED_SEG_SHIFT; + ++seg[sid].u[i]; + ++seg[sid].n_a; } } + for (s = 0; s < n_segs; ++s) { + mm_seg_t *sr = &seg[s]; + for (i = 0, sr->n_u = 0; i < n_regs0; ++i) // squeeze out zero-length per-segment chains + if ((int32_t)sr->u[i] != 0) + sr->u[sr->n_u++] = sr->u[i]; + sr->a = (mm128_t*)kmalloc(km, sr->n_a * sizeof(mm128_t)); + sr->n_a = 0; + } + + for (i = 0; i < n_regs0; ++i) { + const mm_reg1_t *r = ®s0[i]; + for (j = 0; j < r->cnt; ++j) { + int sid = (a[r->as + j].y&MM_SEED_SEG_MASK)>>MM_SEED_SEG_SHIFT; + mm128_t a1 = a[r->as + j]; + // on reverse strand, the segment position is: + // x_for_cat = qlen_sum - 1 - (int32_t)a1.y - 1 + q_span + // (int32_t)new_a1.y = qlens[sid] - (x_for_cat - acc_qlen[sid] + 1 - q_span) - 1 = (int32_t)a1.y - (qlen_sum - (qlens[sid] + acc_qlen[sid])) + a1.y -= a1.x>>63? qlen_sum - (qlens[sid] + acc_qlen[sid]) : acc_qlen[sid]; + seg[sid].a[seg[sid].n_a++] = a1; + } + } + for (s = 0; s < n_segs; ++s) { + regs[s] = mm_gen_regs(km, qlens[s], seg[s].n_u, seg[s].u, seg[s].a); + n_regs[s] = seg[s].n_u; + } + return seg; } -*/ + +void mm_seg_free(void *km, int n_segs, mm_seg_t *segs) +{ + int i; + for (i = 0; i < n_segs; ++i) kfree(km, segs[i].u); + for (i = 0; i < n_segs; ++i) kfree(km, segs[i].a); + kfree(km, segs); +} + void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len) { static const float q_coef = 40.0f; diff --git a/map.c b/map.c index af0f78a..4778b07 100644 --- a/map.c +++ b/map.c @@ -154,20 +154,18 @@ static int mm_dust_minier(int n, mm128_t *a, int l_seq, const char *seq, int sdu return k; // the new size } -static int collect_minimizers(const mm_mapopt_t *opt, const mm_idx_t *mi, int n_segs, int *qlens, const char *seqs, mm_tbuf_t *b) +static void collect_minimizers(const mm_mapopt_t *opt, const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, mm_tbuf_t *b) { - int i, j, n; - const char *seq; - b->mini.n = 0, seq = seqs; + int i, j, n, sum = 0; + b->mini.n = 0; for (i = n = 0; i < n_segs; ++i) { - mm_sketch(b->km, seq, qlens[i], mi->w, mi->k, i, mi->is_hpc, &b->mini); + mm_sketch(b->km, seqs[i], qlens[i], mi->w, mi->k, i, mi->is_hpc, &b->mini); for (j = n; j < b->mini.n; ++j) - b->mini.a[j].x += (seq - seqs) << 1; + b->mini.a[j].x += sum << 1; if (opt->sdust_thres > 0) // mask low-complexity minimizers - b->mini.n = n + mm_dust_minier(b->mini.n - n, b->mini.a + n, qlens[i], seq, opt->sdust_thres, b->sdb); - seq += qlens[i], n = b->mini.n; + b->mini.n = n + mm_dust_minier(b->mini.n - n, b->mini.a + n, qlens[i], seqs[i], opt->sdust_thres, b->sdb); + sum += qlens[i], n = b->mini.n; } - return seq - seqs; } static mm128_t *collect_seed_hits(const mm_mapopt_t *opt, const mm_idx_t *mi, const char *qname, int qlen, int64_t *n_a, int *rep_len, mm_tbuf_t *b) @@ -230,20 +228,42 @@ static mm128_t *collect_seed_hits(const mm_mapopt_t *opt, const mm_idx_t *mi, co return a; } -mm_reg1_t *mm_map_seg(const mm_idx_t *mi, int n_segs, int *qlens, const char *seqs, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname) +static void chain_post(const mm_mapopt_t *opt, const mm_idx_t *mi, void *km, int qlen, int *n_regs, mm_reg1_t *regs, mm128_t *a) { - int i, j, n_u, max_gap_ref, rep_len, qlen_sum; + if (!(opt->flag & MM_F_AVA)) { // don't choose primary mapping(s) for read overlap + mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b); + mm_select_sub(km, opt->mask_level, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs); + if (!(opt->flag & MM_F_SPLICE) && !(opt->flag & MM_F_SR)) + mm_join_long(km, opt, qlen, n_regs, regs, a); // TODO: this can be applied to all-vs-all in principle + } +} + +static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *km, int qlen, const char *seq, int *n_regs, mm_reg1_t *regs, mm128_t *a) +{ + if (!(opt->flag & MM_F_CIGAR)) return regs; + regs = mm_align_skeleton(km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs() + if (!(opt->flag & MM_F_AVA)) { + mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b); + mm_select_sub(km, opt->mask_level, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs); + mm_set_sam_pri(*n_regs, regs); + } + return regs; +} + +void mm_map_multi(const mm_idx_t *mi, int n_segs, const int *qlens, const char **seqs, int *n_regs, mm_reg1_t **regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname) +{ + int i, j, n_u, max_gap_ref, rep_len, qlen_sum, n_regs0; int64_t n_a; uint64_t *u; mm128_t *a; - mm_reg1_t *regs; + mm_reg1_t *regs0; - if (n_segs > MM_MAX_SEG || n_segs <= 0) return 0; + for (i = 0, qlen_sum = 0; i < n_segs; ++i) + qlen_sum += qlens[i], n_regs[i] = 0, regs[i] = 0; - qlen_sum = collect_minimizers(opt, mi, n_segs, qlens, seqs, b); + if (qlen_sum == 0 || n_segs <= 0 || n_segs > MM_MAX_SEG) return; - *n_regs = 0; - if (qlen_sum == 0) return 0; + collect_minimizers(opt, mi, n_segs, qlens, seqs, b); a = collect_seed_hits(opt, mi, qname, qlen_sum, &n_a, &rep_len, b); radix_sort_128x(a, a + n_a); @@ -267,45 +287,41 @@ mm_reg1_t *mm_map_seg(const mm_idx_t *mi, int n_segs, int *qlens, const char *se kfree(b->km, a0); } } - regs = mm_gen_regs(b->km, qlen_sum, n_u, u, a); - *n_regs = n_u; + regs0 = mm_gen_regs(b->km, qlen_sum, n_u, u, a); + n_regs0 = n_u; if (mm_dbg_flag & MM_DBG_PRINT_SEED) for (j = 0; j < n_u; ++j) - for (i = regs[j].as; i < regs[j].as + regs[j].cnt; ++i) + for (i = regs0[j].as; i < regs0[j].as + regs0[j].cnt; ++i) fprintf(stderr, "CN\t%d\t%s\t%d\t%c\t%d\t%d\t%d\n", j, mi->seq[a[i].x<<1>>33].name, (int32_t)a[i].x, "+-"[a[i].x>>63], (int32_t)a[i].y, (int32_t)(a[i].y>>32&0xff), - i == regs[j].as? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x)); + i == regs0[j].as? 0 : ((int32_t)a[i].y - (int32_t)a[i-1].y) - ((int32_t)a[i].x - (int32_t)a[i-1].x)); - 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, opt->a * 2 + opt->b); - mm_select_sub(b->km, opt->mask_level, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs); - if (!(opt->flag & MM_F_SPLICE) && !(opt->flag & MM_F_SR)) - mm_join_long(b->km, opt, qlen_sum, n_regs, regs, a); // TODO: this can be applied to all-vs-all in principle - } + chain_post(opt, mi, b->km, qlen_sum, &n_regs0, regs0, a); - if (opt->flag & MM_F_CIGAR) { - if (n_segs == 1) { - regs = mm_align_skeleton(b->km, opt, mi, qlen_sum, seqs, n_regs, regs, a); // this calls mm_filter_regs() - if (!(opt->flag & MM_F_AVA)) { - mm_set_parent(b->km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b); - mm_select_sub(b->km, opt->mask_level, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs); - mm_set_sam_pri(*n_regs, regs); - } - } else { - abort(); + if (n_segs == 1) { + regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], &n_regs0, regs0, a); + mm_set_mapq(n_regs0, regs0, opt->min_chain_score, opt->a, rep_len); + n_regs[0] = n_regs0, regs[0] = regs0; + } else { + mm_seg_t *seg; + seg = mm_seg_gen(b->km, n_segs, qlens, n_regs0, regs0, n_regs, regs, a); + free(regs0); + for (i = 0; i < n_segs; ++i) { + regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a); + mm_set_mapq(n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len); } + mm_seg_free(b->km, n_segs, seg); } - mm_set_mapq(*n_regs, regs, opt->min_chain_score, opt->a, rep_len); - // free kfree(b->km, a); kfree(b->km, u); - return regs; } mm_reg1_t *mm_map(const mm_idx_t *mi, int qlen, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *qname) { - return mm_map_seg(mi, 1, &qlen, seq, n_regs, b, opt, qname); + mm_reg1_t *regs; + mm_map_multi(mi, 1, &qlen, &seq, n_regs, ®s, b, opt, qname); + return regs; } /************************** diff --git a/mmpriv.h b/mmpriv.h index 19caa0a..0afb432 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -40,10 +40,10 @@ typedef struct __kstring_t { #endif typedef struct { - int n_regs, n_a; - mm_reg1_t *regs; + int n_u, n_a; + uint64_t *u; mm128_t *a; -} mm_segreg_t; +} mm_seg_t; double cputime(void); double realtime(void); @@ -77,6 +77,9 @@ void mm_join_long(void *km, const mm_mapopt_t *opt, int qlen, int *n_regs, mm_re void mm_hit_sort_by_dp(void *km, int *n_regs, mm_reg1_t *r); void mm_set_mapq(int n_regs, mm_reg1_t *regs, int min_chain_sc, int match_sc, int rep_len); +mm_seg_t *mm_seg_gen(void *km, int n_segs, const int *qlens, int n_regs0, const mm_reg1_t *regs0, int *n_regs, mm_reg1_t **regs, const mm128_t *a); +void mm_seg_free(void *km, int n_segs, mm_seg_t *segs); + #ifdef __cplusplus } #endif