From 236bc7b9ee173bf8d293101742f5493f107dc606 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Wed, 28 Sep 2022 11:15:29 -0400 Subject: [PATCH] backup for minimizer --- Correct.cpp | 299 +++++++++++++++++++++++++++++++++++++---- Levenshtein_distance.h | 8 +- inter.cpp | 37 ++--- 3 files changed, 289 insertions(+), 55 deletions(-) diff --git a/Correct.cpp b/Correct.cpp index d35e34a..0b73d00 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -12626,21 +12626,45 @@ int64_t qs, int64_t qe, int64_t thre, int64_t *ts, int64_t *te, int64_t *aux_beg return 1; } -int64_t cal_exz_infi(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, bit_extz_t *exz, char* qstr, UC_Read *tu, int64_t qs, int64_t qe, int64_t ts, int64_t te, int64_t thre, int64_t mode) + +void adjust_ext_offset(int64_t *qs, int64_t *qe, int64_t *ts, int64_t *te, int64_t ql, int64_t tl, int64_t thre, int64_t mode) +{ + int64_t qoff, toff; + if(mode == 1) {///forward extension + qoff = ql - (*qs); toff = tl - (*ts); + if(qoff <= toff) { + (*qe) = ql; (*te) = (*ts) + qoff + thre; + } else { + (*te) = tl; (*qe) = (*qs) + toff + thre; + } + } else if(mode == 2) {///backward extension + qoff = (*qe); toff = (*te); + if(qoff <= toff) { + (*qs) = 0; (*ts) = (*te) - qoff - thre; + } else { + (*ts) = 0; (*qs) = (*qe) - toff - thre; + } + } + if((*qs) < 0) (*qs) = 0; + if((*ts) < 0) (*ts) = 0; + if((*qe) > ql) (*qe) = ql; + if((*te) > tl) (*te) = tl; +} + +int64_t cal_exz_infi(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, bit_extz_t *exz, char* qstr, UC_Read *tu, int64_t qs, int64_t qe, int64_t ts, int64_t te, int64_t thre, int64_t q_tot_l, int64_t mode) { int64_t aux_beg = 0, bd = (((thre)<<1)+1), ql, tl, t_tot_l = -1; int32_t nword = ((bd>>bitw)+(!!(bd&bitz))); char *q_string, *t_string; int32_t rev = z->y_pos_strand, id = z->y_id; ql = qe - qs; + if(hpc_g) t_tot_l = hpc_len(*hpc_g, id); + else if(uref) t_tot_l = uref->ug->u.a[id].len; + else t_tot_l = Get_READ_LENGTH((*rref), id); + if(mode == 3) { update_semi_coord(uref, hpc_g, rref, z, qs, qe, thre, &ts, &te, &aux_beg); - } else if(mode == 2) { - ts = te - ql - thre; if(ts < 0) ts = 0; - } else if(mode == 1) { - te = ts + ql + thre; - if(hpc_g) t_tot_l = hpc_len(*hpc_g, id); - else if(uref) t_tot_l = uref->ug->u.a[id].len; - else t_tot_l = Get_READ_LENGTH((*rref), id); - if(te > t_tot_l) te = t_tot_l; + } else if(mode == 1 || mode == 2) { + adjust_ext_offset(&qs, &qe, &ts, &te, q_tot_l, t_tot_l, thre, mode); } + if((qe > qs) && (te > ts) && (ts != -1) && (te != -1)) { ql = qe - qs; q_string = qstr + qs; tl = te - ts; resize_UC_Read(tu, tl); @@ -12697,16 +12721,16 @@ int64_t cal_exz_infi(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_ void hc_aln_exz(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, int64_t qs, int64_t qe, int64_t ts, int64_t te, int64_t estimate_err, -int64_t mode, int64_t wl, bit_extz_t *exz, double e_rate) +int64_t mode, int64_t wl, bit_extz_t *exz, int64_t q_tot, double e_rate) { int64_t thre, ql = qe - qs, thre0; if(((ts == -1) && (te == -1))) mode = 3;///set to semi-global // fprintf(stderr, "[M::%s::ql::%ld] qs::%ld, qe::%ld, ts::%ld, te::%ld, mode::%ld, estimate_err::%ld, e_rate::%f", // __func__, ql, qs, qe, ts, te, mode, estimate_err, e_rate); - if(ql <= MAX_L) { + if(ql <= MAX_L && (estimate_err*1.2) <= MAX_E) { thre = scale_ed_thre(estimate_err); if(thre > ql) thre = ql; - if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, mode)) { + if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, q_tot, mode)) { // fprintf(stderr, ", err::%d, thre::%d, scale::%ld(+)\n", exz->err, exz->thre, thre); return; } @@ -12714,7 +12738,7 @@ int64_t mode, int64_t wl, bit_extz_t *exz, double e_rate) thre0 = thre; thre = ql*e_rate; thre = scale_ed_thre(thre); if(thre > ql) thre = ql; if(thre > thre0) { - if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, mode)) { + if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, q_tot, mode)) { // fprintf(stderr, ", err::%d, thre::%d, scale::%ld(-)\n", exz->err, exz->thre, thre); return; } @@ -12723,7 +12747,7 @@ int64_t mode, int64_t wl, bit_extz_t *exz, double e_rate) thre0 = thre; thre <<= 1; thre = scale_ed_thre(thre); if(thre > ql) thre = ql; if(thre > thre0) { - if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, mode)) { + if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, q_tot, mode)) { // fprintf(stderr, ", err::%d, thre::%d, scale::%ld(-)\n", exz->err, exz->thre, thre); return; } @@ -12732,19 +12756,19 @@ int64_t mode, int64_t wl, bit_extz_t *exz, double e_rate) thre0 = thre; thre = ql*0.51; thre = scale_ed_thre(thre); if(thre > ql) thre = ql; if(thre > thre0) { - if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, mode)) { - // fprintf(stderr, ", err::%d, thre::%d, scale::%ld(-)\n", exz->err, exz->thre, thre); + if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, q_tot, mode)) { + // fprintf(stderr, ", err::%d, thre::%d, scale::%ld(*)\n", exz->err, exz->thre, thre); return; } } - - // fprintf(stderr, ", err::%d, thre::%d\n", INT32_MAX, exz->thre); } + // fprintf(stderr, ", err::%d, thre::%d\n", INT32_MAX, exz->thre); } void sub_ciagar_gen(overlap_region *z, uint64_t s, uint64_t e, uint64_t wl, -const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate) +const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate, +int64_t ql, uint64_t rid) { uint64_t qs, qe, sid, eid, k, l, m, tot_e, c_e; int64_t q[2], t[2], mode; qs = (s/wl)*wl; if(qs < z->x_pos_s) qs = z->x_pos_s; if(qs > z->x_pos_e) return; @@ -12813,7 +12837,7 @@ const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bi // fprintf(stderr, "[M::%s::ql::%lu] qs::%lu, qe::%lu, ts::%lu, te::%lu, mode::%ld, tot_e::%lu\n", // __func__, q[1]-q[0], q[0], q[1], t[0], t[1], mode, tot_e); // } - hc_aln_exz(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], tot_e, mode, wl, exz, e_rate); + hc_aln_exz(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], tot_e, mode, wl, exz, ql, e_rate); } l = k; @@ -12822,27 +12846,248 @@ const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bi } void cigar_gen(overlap_region *z, ul_ov_t *ov, uint64_t on, uint64_t qn, uint64_t wl, -const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate) +const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate, int64_t ql, uint64_t rid) { - uint64_t i; + uint64_t i, qs = (uint64_t)-1, qe = (uint64_t)-1; for (i = 0; i < on && ov[i].qn == qn; i++) { - sub_ciagar_gen(z, ov[i].qs, ov[i].qe, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate); + assert((i<=0)||(ov[i].qs >= ov[i-1].qe)); + if(ov[i].qs <= qe && qe != (uint64_t)-1) { + qe = ov[i].qe; + } else { + if(qs != (uint64_t)-1) sub_ciagar_gen(z, qs, qe, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); + qs = ov[i].qs; qe = ov[i].qe; + } } - + if(qs != (uint64_t)-1) sub_ciagar_gen(z, qs, qe, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); } void ul_gap_filling(overlap_region_alloc* ol, kv_ul_ov_t *aln, uint64_t wl, -const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate) +const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, +double e_rate, int64_t ql, uint64_t rid) { int64_t i, on = ol->length; for (i = 0; i < on; i++) { if(ol->list[i].align_length == (uint32_t)-1) continue; cigar_gen(&(ol->list[i]), aln->a+ol->list[i].align_length, aln->n-ol->list[i].align_length, i, wl, - uref, hpc_g, rref, qstr, tu, exz, e_rate); + uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); } } +///[ys, ye) +uint64_t get_win_aln(overlap_region *z, uint64_t wid, int64_t *ys, int64_t *ye, int64_t *err) +{ + (*err) = -2; + if((wid > 0) && (z->w_list.a[wid].y_end != -1) && (z->w_list.a[wid-1].y_end != -1) && (z->w_list.a[wid].y_end > z->w_list.a[wid-1].y_end)) { + (*ys) = z->w_list.a[wid-1].y_end+1; + (*ye) = z->w_list.a[wid].y_end+1; + (*err) = z->w_list.a[wid].error; + return 1; + } + return 0; +} + +char* retrive_str_piece_exz(All_reads *rref, const ul_idx_t *uref, char *buf, int64_t s, int64_t l, int64_t rev, int64_t id) +{ + if(rref) recover_UC_Read_sub_region(buf, s, l, rev, rref, id); + else if(uref) retrieve_u_seq(NULL, buf, &(uref->ug->u.a[id]), rev, s, l, NULL); + else return NULL; + return buf; +} + +uint64_t dp_commen_sketch(kv_ul_ov_t *aln, overlap_region_alloc* ol, uint64_t *id_a, int64_t id_n, +uint64_t *win_a, int64_t win_n, uint64_t *dp, int64_t n_skip, int64_t wl) +{ + if(!win_n) return 0; + int64_t k, ws, ws0, i, m, wid, wid0, s, e, err, s0, e0, err0, sc, max_sc, p, long_sc, long_idx; + overlap_region *z; s = e = err = s0 = e0 = err0 = -1; + if(n_skip < win_n) { + long_sc = long_idx = -1; + for (k = 0; k < n_skip; k++) { + dp[k] = ((k>0)?(k-1):((uint32_t)-1)); dp[k] <<= 32; dp[k] |= k+1; max_sc = k+1; + if(long_sc < max_sc) { + long_sc = max_sc; long_idx = k; + } + } + + for (; k < win_n; k++) { + ws = win_a[k]; sc = 1; p = -1; max_sc = sc; + for (m = k-1; m >= 0; m--) { + ws0 = win_a[m]; + for (i = 0; i < id_n; i++) { + z = &(ol->list[aln->a[(uint32_t)id_a[i]].qn]); + wid = get_win_id_by_s(z, ws, wl, NULL); + wid0 = get_win_id_by_s(z, ws0, wl, NULL); + + get_win_aln(z, wid, &s, &e, &err); + get_win_aln(z, wid0, &s0, &e0, &err0); + if(e0 > s) break; + } + if((i >= id_n) && ((sc + ((uint32_t)dp[m])) > max_sc)) { + max_sc = (sc + ((uint32_t)dp[m])); p = m; + } + } + dp[k] = ((p>=0)?(p):((uint32_t)-1)); dp[k] <<= 32; dp[k] |= max_sc; + if(long_sc < max_sc) { + long_sc = max_sc; long_idx = k; + } + } + + k = long_idx; + while (k >= 0) { + win_a[k] |= ((uint64_t)0x8000000000000000); + k = ((dp[k]>>32)!=((uint32_t)-1))?(dp[k]>>32):(-1); + } + + for (k = 0, m = 0; k < win_n; k++) { + if(!(win_a[k]&((uint64_t)0x8000000000000000))) continue; + win_a[m++] = win_a[k]; + } + win_n = m; + } + + for (k = 0; k < win_n; k++) { + ws = win_a[k]; + for (i = 0; i < id_n; i++) { + z = &(ol->list[aln->a[(uint32_t)id_a[i]].qn]); + wid = get_win_id_by_s(z, ws, wl, NULL); + z->w_list.a[wid].extra_end = -1; + get_win_aln(z, wid, &s, &e, &err); + if(z->align_length < (uint32_t)e) z->align_length = e;///for conliner + } + } + return win_n; +} + +uint64_t gen_commen_sketch(All_reads *rref, const ul_idx_t *uref, overlap_region_alloc* ol, uint64_t *id_a, uint64_t id_n, uint64_t s, uint64_t e, uint64_t ql, uint64_t wl, +uint64_t *buf, uint64_t dp, char *str0, char *str1, kv_ul_ov_t *aln, asg64_v *trace)///[s, e) +{ + if(!id_n) return id_n; + uint64_t i, m, k, rm_n = 0, buf_n = 0, qs, qe, wid, co, occ; char *qstring, *tstring; + overlap_region *z; uint64_t ws, we; int64_t r_y[2], r_err, p_y[2], p_err; + ///shrink [qs, qe) + qs = (s/wl)*wl; if(qs < s) qs += wl; if(qs >= ql) return id_n; + qe = (e/wl)*wl; if(qe >= ql) qe = ql; + if(qs >= qe) return id_n; + //idx_a[] is sorted by aln[].qs + for (k = 0; k < id_n; k++) { + if(aln->a[id_a[k]].qs<=qs && aln->a[id_a[k]].qe>=qe) { + buf[buf_n++] = id_a[k]; + } + if(aln->a[id_a[k]].qe < e) rm_n++; + } + assert(buf_n == dp && buf_n > 1); + + if(buf_n > 0) { + ///fs = fe = (uint64_t)-1; + co = 1; trace->n = occ = 0; + for (k = qs; k < qe; k += wl) { + ws = k; we = ws + wl; if(we > qe) we = qe;//[ws, we) + ///first overlap + z = &(ol->list[aln->a[(uint32_t)buf[0]].qn]); + wid = get_win_id_by_s(z, ws, wl, NULL); + if(!get_win_aln(z, wid, &(r_y[0]), &(r_y[1]), &r_err)) continue; + if(r_y[0] < z->align_length) continue;///not co-linear + qstring = tstring = NULL; + for (i = 1; i < buf_n; i++) { + z = &(ol->list[aln->a[(uint32_t)buf[i]].qn]); + wid = get_win_id_by_s(z, ws, wl, NULL); + if(!get_win_aln(z, wid, &(p_y[0]), &(p_y[1]), &p_err)) break; + if(p_y[0] < z->align_length) continue;///not co-linear + + if(((r_y[1]-r_y[0]) != (p_y[1]-p_y[0])) || (r_err != p_err)) break; + if(r_err == 0) continue; + if(!qstring) { + qstring = retrive_str_piece_exz(rref, uref, str0, r_y[0], r_y[1]-r_y[0], + ol->list[aln->a[(uint32_t)buf[0]].qn].y_pos_strand, ol->list[aln->a[(uint32_t)buf[0]].qn].y_id); + } + tstring = retrive_str_piece_exz(rref, uref, str1, p_y[0], p_y[1]-p_y[0], z->y_pos_strand, z->y_id); + if(memcmp(str0, str1, (p_y[1]-p_y[0]))) { + // fprintf(stderr, "[M::%s::] qs::%ld, qe::%lu, ts::%ld, te::%lu, err::%ld, rts::%ld, rte::%lu, err::%ld\n", + // __func__, ws, we, p_y[0], p_y[1], r_err, r_y[0], r_y[1], p_err); + // fprintf(stderr, "str0::%.*s\n", (int32_t)(r_y[1]-r_y[0]), str0); + // fprintf(stderr, "str1::%.*s\n", (int32_t)(p_y[1]-p_y[0]), str1); + break; + } + } + if(i < buf_n) continue; + for (i = 0, p_y[0] = p_y[1] = -1; i < buf_n; i++) { + z = &(ol->list[aln->a[(uint32_t)buf[i]].qn]); + wid = get_win_id_by_s(z, ws, wl, NULL); + + if(co) { + get_win_aln(z, wid, &(p_y[0]), &(p_y[1]), &p_err); + if((buf[i]>>32) > (uint64_t)p_y[0]) co = 0; + m = p_y[1]; m <<= 32; m |= (uint32_t)buf[i]; buf[i] = m; + } + } + if(co) occ++; + kv_push(uint64_t, *trace, ws); + } + + if(trace->n) { + if(!co) kv_resize(uint64_t, *trace, trace->n<<1); + trace->n = dp_commen_sketch(aln, ol, buf, buf_n, trace->a, trace->n, trace->a + trace->n, occ, wl); + } + } + if(rm_n) { + for (i = m = 0; i < id_n; i++) { + if(aln->a[id_a[i]].qe < e) continue; + id_a[m++] = id_a[i]; + } + id_n = m; + } + return id_n; +} + +void update_sketch_trace(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t *uopt, +All_reads *rref, UC_Read* tu, asg64_v* idx, asg64_v *b0, asg64_v *b1, int64_t ql, int64_t wl, kv_ul_ov_t *aln, uint64_t rid) +{ + if(!aln->n) return; + uint64_t i, k, own, srt_n; int64_t dp, old_dp, beg, end; overlap_region *z; + for (i = 0; i < ol->length; i++) { + z = &(ol->list[i]); append_unmatched_wins(z, wl); + own = z->w_list.n; z->align_length = (uint32_t)-1; + for (k = 0; k < own; k++) { + if(z->w_list.a[k].extra_end < 0) z->w_list.a[k].extra_end = 0; + } + } + + kv_resize(uint64_t, *idx, (aln->n<<1)); kv_resize(uint64_t, *b0, aln->n); + for (i = srt_n = 0; i < aln->n; i++) { + // if(i == 0 || aln->a[i].qn != aln->a[i-1].qn) ol->list[aln->a[i].qn].align_length = i; + ol->list[aln->a[i].qn].align_length = 0;///for co-linear + idx->a[srt_n] = aln->a[i].qs<<1; idx->a[srt_n] <<= 32; idx->a[srt_n] += i; srt_n++; + idx->a[srt_n] = ((aln->a[i].qe-1)<<1)+1; idx->a[srt_n] <<= 32; idx->a[srt_n] += i; srt_n++; + aln->a[i].el = 1; + } + + radix_sort_bc64(idx->a, idx->a+srt_n); idx->n = srt_n; resize_UC_Read(tu, (wl<<1)); + for (i = 0, dp = 0, beg = 0, end = -1; i < srt_n; ++i) {///[beg, end] + old_dp = dp; + ///if idx->a.a[] is qe + if ((idx->a[i]>>32)&1) { + --dp; end = (idx->a[i]>>33)+1; + }else { + //meet a new overlap; the overlaps are pushed by the x_pos_s + ++dp; end = (idx->a[i]>>33); + kv_push(uint64_t, *idx, ((uint32_t)idx->a[i])); + } + + if((end > beg) && (end - beg > wl) && (old_dp >= 2) ) { + idx->n = srt_n + + gen_commen_sketch(rref, uref, ol, idx->a+srt_n, idx->n-srt_n, beg, end, ql, wl, b0->a, old_dp, tu->seq, tu->seq+wl, aln, b1); + } + beg = end; + } + + for (i = srt_n = 0; i < aln->n; i++) { + if(i == 0 || aln->a[i].qn != aln->a[i-1].qn) ol->list[aln->a[i].qn].align_length = i; + } + + return; +} + void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, char *qstr, uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, bit_extz_t *exz, @@ -12898,7 +13143,7 @@ void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur } else { // fprintf(stderr, "-[M::%s] on::%lu\n", __func__, ol->length); if(ol->length <= 1) return; - ul_gap_filling(ol, aln, wl, uref, NULL, NULL, qu->seq, tu, exz, err); + ul_gap_filling(ol, aln, wl, uref, NULL, NULL, qu->seq, tu, exz, err, ql, sid); // for (i = 0; (i < ol->length) && (ol->list[i].is_match == 1); i++); on = i; // if(on <= 1) return; // kv_resize(uint64_t, v_idx->a, (on<<1)); v_idx->a.n = on; diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index 6755b5b..7938777 100644 --- a/Levenshtein_distance.h +++ b/Levenshtein_distance.h @@ -538,8 +538,12 @@ inline uint32_t pop_trace(asg16_v *res, uint32_t i, uint16_t *c, uint32_t *len) } ///511 -> 16 64-bits -#define MAX_E 511 -#define MAX_L 2500 +// #define MAX_E 511 +// #define MAX_L 2500 + +///511 -> 32 64-bits +#define MAX_E 1023 +#define MAX_L 8192 typedef uint64_t w_sig; #define bitw (6) diff --git a/inter.cpp b/inter.cpp index dd7b6e8..6a3172c 100644 --- a/inter.cpp +++ b/inter.cpp @@ -5872,6 +5872,7 @@ uint64_t *buf, uint64_t dp, char *str0, char *str1, kv_ul_ov_t *aln)///[s, e) return id_n; } + void update_shared_intervals(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t *uopt, All_reads *rref, UC_Read* tu, kvec_t_u64_warp* idx, st_mt_t *sps, int64_t ql, int64_t wl, kv_ul_ov_t *aln, uint64_t rid) { @@ -5916,30 +5917,6 @@ All_reads *rref, UC_Read* tu, kvec_t_u64_warp* idx, st_mt_t *sps, int64_t ql, in beg = end; } - - - - // for (i = 1, j = 0; i <= aln->n; i++) { - // if (i == aln->n || aln->a[i].qn != aln->a[j].qn) { - // z = &(ol->list[aln->a[j].qn]); own = z->w_list.n; dp = old_dp = 0; - // for (k = j; k < i; k++) { - // old_dp += aln->a[k].qe-aln->a[k].qs; - // assert(k <= j || aln->a[k].qs >= aln->a[k-1].qe); - // } - // for (k = 0; k < own; k++) { - // if(z->w_list.a[k].extra_end < 0) dp += z->w_list.a[k].x_end+1-z->w_list.a[k].x_start; - // } - // fprintf(stderr, "[M::%s::utg%.6dl] pre_len::%ld, cur_len::%ld\n", __func__, (int32_t)z->y_id+1, old_dp, dp); - // j = i; - // } - // } - - // if(aln->n > aln_n) {///this function works - // idx->a.n = aln->n - aln_n; kv_resize(uint64_t, idx->a, idx->a.n); - // for (i = aln_n; i < aln->n; i++) { - // fprintf(stderr, "***[M::%s::shared] q[%u, %u)\n", __func__, aln->a[i].qs, aln->a[i].qe); - // } - // } return; } @@ -6054,6 +6031,7 @@ static void worker_for_ul_scall_alignment(void *data, long i, int tid) // callba } + static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // callback for kt_for() { utepdat_t *s = (utepdat_t*)data; @@ -6061,6 +6039,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call glchain_t *bl = &(s->ll[tid]); int64_t /**rid = s->id+i,**/ winLen = MIN((((double)THRESHOLD_MAX_SIZE)/s->opt->diff_ec_ul), WINDOW), ton = 0; uint32_t high_occ = 2, phase = 1; + asg64_v b0, b1, b2; // uint64_t align = 0; // if(UL_INF.a[s->id+i].rlen != s->len[i]) { @@ -6071,7 +6050,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // if(s->id+i!=41927 && s->id+i!=47072 && s->id+i!=67641 && s->id+i!=90305 && s->id+i!=698342 && s->id+i!=329421) { // return; // } - // if((s->id+i!=2936) /**&& (s->id+i!=44) && (s->id+i!=948)**/) return; + // if((s->id+i!=230) /**&& (s->id+i!=44) && (s->id+i!=948)**/) return; // fprintf(stderr, "\n[M::%s] rid::%ld, len::%lu, name::%.*s\n", __func__, s->id+i, s->len[i], // (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a); @@ -6105,7 +6084,13 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a); if(gen_shared_intervals(&b->olist, s->uu, s->uopt, winLen, &b->r_buf, &(bl->lo))) { filter_topN(&b->olist, &(bl->lo), s->len[i], winLen, UL_TOPN, bl); - update_shared_intervals(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b->r_buf, &(s->sps[tid]), s->len[i], winLen, &(bl->lo), s->id+i); + // update_shared_intervals(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b->r_buf, &(s->sps[tid]), s->len[i], winLen, &(bl->lo), s->id+i); + + copy_asg_arr(b0, b->r_buf.a); copy_asg_arr(b1, s->sps[tid]); copy_asg_arr(b2, b->hap.snp_srt); + // update_shared_intervals(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b0, &b1, &b2, s->len[i], winLen, &(bl->lo), s->id+i); + update_sketch_trace(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b0, &b1, &b2, s->len[i], winLen, &(bl->lo), s->id+i); + copy_asg_arr(b->r_buf.a, b0); copy_asg_arr(s->sps[tid], b1); copy_asg_arr(b->hap.snp_srt, b2); + ul_lalign(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, &b->correct, &b->exz, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, &(bl->lo), s->id+i, NULL); // ul_lalign_old_ed(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read,