From 8a4d6a5f24c4a854ceefdc8f68d00d55f3304559 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Tue, 4 Oct 2022 15:31:23 -0400 Subject: [PATCH] done alignemnt --- Correct.cpp | 439 +++++++++++++++++++++++++++++++++++++++++++++---- Correct.h | 3 + Hash_Table.cpp | 263 ++++++++++++++++++++++++++++- Hash_Table.h | 7 + 4 files changed, 683 insertions(+), 29 deletions(-) diff --git a/Correct.cpp b/Correct.cpp index da021d9..8c5f75a 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -35,8 +35,14 @@ KRADIX_SORT_INIT(window_list_xs_srt, window_list, window_list_xs_key, member_siz #define uov_qs_key(p) ((p).qs) KRADIX_SORT_INIT(uov_srt_qs, ul_ov_t, uov_qs_key, member_size(ul_ov_t, qs)) -int ha_ov_type(const overlap_region *r, uint32_t len); +#define k_mer_hit_self_key(p) ((p).self_offset) +KRADIX_SORT_INIT(k_mer_hit_self, k_mer_hit, k_mer_hit_self_key, member_size(k_mer_hit, self_offset)) +#define k_mer_hit_off_key(p) ((p).offset) +KRADIX_SORT_INIT(k_mer_hit_off, k_mer_hit, k_mer_hit_off_key, member_size(k_mer_hit, offset)) + +int ha_ov_type(const overlap_region *r, uint32_t len); +void set_lchain_dp_op(uint32_t is_accurate, uint32_t mz_k, int64_t *max_skip, int64_t *max_iter, int64_t *max_dis, double *chn_pen_gap, double *chn_pen_skip, int64_t *quick_check); void clear_Round2_alignment(Round2_alignment* h) { @@ -13797,17 +13803,18 @@ void append_wcigar(window_list *idx, window_list_alloc *res, bit_extz_t *exz) void push_alnw(overlap_region *aux_o, bit_extz_t *exz) { - window_list *p = NULL; + window_list *p = NULL; int64_t t; if(aux_o->w_list.n > 0) { p = &(aux_o->w_list.a[aux_o->w_list.n-1]); // fprintf(stderr, "+[M::%s::wn->%d] px::[%d, %d], py::[%d, %d], pe::%d, exz->t::[%u, %u], exz->p::[%u, %u], exz->e::%d, clen::%u\n", // __func__, (int32_t)(aux_o->w_list.n), p->x_start, p->x_end, p->y_start, p->y_end, p->error, // exz->ts, exz->te, exz->ps, exz->pe, exz->err, p->clen); - assert((p->x_endts)&&(p->y_endps)); + // assert((p->x_endts)&&(p->y_endps)); if(p->clen > 0) { - if(((p->x_end+1) == exz->ts) && ((p->y_end+1) == exz->ps) - && ((p->error+exz->err) <= INT16_MAX)) { + t = ((int64_t)p->error) + ((int64_t)exz->err); + ///note: t cannot be equal to INT16_MAX; otherwise it is unable to distiguish unaligned regions + if(((p->x_end+1) == exz->ts) && ((p->y_end+1) == exz->ps) && (t < INT16_MAX)) { p->x_end = exz->te; p->y_end = exz->pe; p->error += exz->err; append_wcigar(p, &(aux_o->w_list), exz); // fprintf(stderr, "-[M::%s::wn->%d] px::[%d, %d], py::[%d, %d], pe::%d, exz->t::[%u, %u], exz->p::[%u, %u], exz->e::%d, clen::%u\n", @@ -13822,7 +13829,7 @@ void push_alnw(overlap_region *aux_o, bit_extz_t *exz) kv_pushp(window_list, aux_o->w_list, &p); p->x_start = exz->ts; p->x_end = exz->te; p->y_start = exz->ps; p->y_end = exz->pe; - p->error_threshold = 0; p->error = exz->err; + p->error_threshold = 0; p->error = exz->err;///single round of alignment cannot have INT16_MAX errors push_wcigar(p, &(aux_o->w_list), exz); } @@ -13954,7 +13961,7 @@ int64_t estimate_err, overlap_region *aux_o) } if(ql <= force_l) { - thre = maxe; + thre = maxe; if(cal_exz_infi_adv(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, &pts, &pte, thre, &pthre, q_tot, mode)) { // ref_cigar_check(qstr, tu, uref, hpc_g, rref, z->y_id, z->y_pos_strand, exz); // fprintf(stderr, ", err::%d, thre::%d, scale::%ld(*)\n", exz->err, exz->thre, thre); @@ -14259,7 +14266,7 @@ int64_t fusion_chain_ovlp(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, ul_o // } -int64_t ovlp_base_aln_all(overlap_region *z, Chain_Data *dp, k_mer_hit *ch_a, int64_t ch_n, +int64_t ovlp_base_aln_all(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, int64_t soff, int64_t eoff, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, ul_ov_t *ov, int64_t ql, int64_t tl, int64_t wl, bit_extz_t *exz, overlap_region *aux_o, double e_rate) @@ -14306,7 +14313,7 @@ overlap_region *aux_o, double e_rate) return 0; } -void ovlp_base_aln(overlap_region *z, Chain_Data *dp, k_mer_hit *ch_a, int64_t ch_n, +void ovlp_base_aln(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, ul_ov_t *ov, int64_t wl, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, uint64_t rid) { @@ -14341,8 +14348,7 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u } else { mode = 3;///no primary hit within [ibeg, iend] } - // fprintf(stderr, "000000[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld, l::%ld, i::%ld, ch_n::%ld\n", - // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode, l, i, ch_n); + if((mode == 0) && is_alnw(ch_a[l]) && is_alnw(ch_a[i]) && (ch_a[l].strand == 0) && (ch_a[i].strand == 1)) { is_done = 1; @@ -14356,26 +14362,383 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); is_done = hc_aln_exz_adv(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, ql, e_rate, MAX_CNS_L, MAX_CNS_E, FORCE_CNS_L, -1, aux_o); } - + if(!is_done) {///postprocess // fprintf(stderr, ">[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n", // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); - is_done = ovlp_base_aln_all(z, dp, ch_a, ch_n, l, i, uref, hpc_g, rref, qstr, tu, ov, ql, tl, wl, exz, aux_o, e_rate); + is_done = ovlp_base_aln_all(z, ch_a, ch_n, l, i, uref, hpc_g, rref, qstr, tu, ov, ql, tl, wl, exz, aux_o, e_rate); } + // if(rid == (uint64_t)-1) { + // fprintf(stderr, "-is_done::%ld[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld, l::%ld, i::%ld, ch_n::%ld\n", + // is_done, __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode, l, i, ch_n); + // } + // if(rid == (uint64_t)-1) { + // fprintf(stderr, "+is_done::%ld[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld, l::%ld, i::%ld, ch_n::%ld\n", + // is_done, __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode, l, i, ch_n); + // } l = i; } } } +inline void push_khit(Candidates_list *res, int32_t xs, int32_t ys, uint32_t len, uint32_t h_khit, uint32_t *ic) +{ + uint32_t p, c; k_mer_hit *z; + c = ((len >= h_khit)?1:2); if(ic) c = *ic; c <<= 8; + if(len > 0) { + while (len >= (0xffu)) { + p = (c + (0xffu)); + kv_pushp_cl(k_mer_hit, (*res), &z); + memset(z, 0, sizeof((*z))); + z->self_offset = xs; z->offset = ys; z->cnt = p; + len -= (0xffu); + } + if(len) { + p = (c + len); + kv_pushp_cl(k_mer_hit, (*res), &z); + memset(z, 0, sizeof((*z))); + z->self_offset = xs; z->offset = ys; z->cnt = p; + } + } else { + p = (c + len); + kv_pushp_cl(k_mer_hit, (*res), &z); + memset(z, 0, sizeof((*z))); + z->self_offset = xs; z->offset = ys; z->cnt = p; + } +} + +uint32_t extract_exact_cigar(asg16_v *ez, int32_t ps, int32_t ts, int32_t pmin, int32_t pmax, +int32_t tmin, int32_t tmax, Candidates_list *res, int32_t minl, int64_t min_w_l, int64_t h_khit) +{ + uint32_t ci = 0, cl, occ = 0; uint16_t c; + int32_t pi = ps, ti = ts, p[2], t[2], poff, toff, maxl; + int32_t pos, poe, tos, toe, l; poff = toff = maxl = -1; + while (ci < ez->n && pi < pmax && ti < tmax) { + ci = pop_trace(ez, ci, &c, &cl); + if(c == 0) { + p[0] = pi; p[1] = pi + cl; + t[0] = ti; t[1] = ti + cl; + pos = MAX(p[0], pmin); poe = MIN(p[1], pmax); + tos = MAX(t[0], tmin); toe = MIN(t[1], tmax); + if((poe > pos) && (toe > tos)) { + l = poe - pos; + if(l == (toe - tos)) { + poe--; toe--; + if(l > maxl) { + poff = poe; toff = toe; maxl = l; + } + if(l >= minl) { + push_khit(res, toe, poe, l, h_khit, NULL); occ++; + } + } + } + pi+=cl; ti+=cl; + } else if(c == 1) { + pi+=cl; ti+=cl; + } else if(c == 2) {///more p + pi+=cl; + } else if(c == 3) { + ti+=cl; + } + } + ///(ts >= tmin) && (ti >= (min_w_l + ts)): here is a whole window + if(maxl > 0 && maxl < minl && (ts >= tmin) && (ti >= (min_w_l + ts))) { + uint32_t w = 3; + push_khit(res, toff, poff, maxl, h_khit, &w); occ++; + } + return occ; +} + +int64_t debug_k_mer_hit_retrive(k_mer_hit *z, hpc_t *hpc_g, All_reads *rref, const ul_idx_t *uref, +char* qstr, UC_Read *tu, int64_t id, int64_t rev) +{ + int64_t qs, qe, ts, te; char *q_string, *t_string; + qe = z->self_offset+1; qs = qe - (z->cnt&(0xffu)); + te = z->offset+1; ts = te - (z->cnt&(0xffu)); + if(qe == qs) return 1; + q_string = qstr + qs; + t_string = retrieve_str_seq_exz(tu, ts, te-ts, -1, -1, rev, uref, hpc_g, rref, id); + fprintf(stderr, "[M::%s::] q_string::%.*s\n", __func__, (int32_t)(qe-qs), q_string); + fprintf(stderr, "[M::%s::] t_string::%.*s\n", __func__, (int32_t)(te-ts), t_string); + if(memcmp(q_string, t_string, qe-qs)) { + fprintf(stderr, "[M::%s::] qsite::%u, tsite::%u\n", __func__, z->self_offset, z->offset); + return 0; + } + return 1; +} + +int64_t gen_single_khit(Candidates_list *cl, int64_t ch_n, int64_t h_khit, int64_t mode, int64_t qs, int64_t qe, int64_t ts, int64_t te, int64_t max_skip, int64_t max_iter) +{ + // if(ch_n != 3 || mode != 0 || qs != 171728 || qe != 172258) return 0; + k_mer_hit *ch_a = cl->list + cl->length; int64_t k, i, j, occ, m, ncn, prefix, suffix, srt = 1; + prefix = suffix = 0; + if(mode == 0 || mode == 2) suffix = 1; + if(mode == 0 || mode == 1) prefix = 1; + // fprintf(stderr, "\n[M::%s::mode->%ld] ch_n::%ld, q::[%ld, %ld), t::[%ld, %ld)\n", + // __func__, mode, ch_n, qs, qe, ts, te); + for (k = occ = 0; k < ch_n; k++) { + // fprintf(stderr, "+i::%ld[M::%s::] x::[%u, %u), y::[%u, %u)\n", k, __func__, + // ch_a[k].self_offset+1-(ch_a[k].cnt&((uint32_t)(0xffu))), ch_a[k].self_offset+1, + // ch_a[k].offset+1-(ch_a[k].cnt&((uint32_t)(0xffu))), ch_a[k].offset+1); + if(!(ch_a[k].cnt&(0xffu))) continue; + occ++; + if((ch_a[k].cnt&(0xffu)) > 1) occ++; + } + occ += prefix + suffix; + + ncn = occ + cl->length; + if(cl->size < ncn) { + cl->size = ncn; + cl->list = (k_mer_hit*)realloc(cl->list, (sizeof((*(cl->list)))*cl->length)); + } + ch_a = cl->list + cl->length; assert((cl->length+occ)<= cl->size); + + k_mer_hit cht; + ///global or backward + if(suffix) { + cht.self_offset = qe; + cht.offset = te; + cht.cnt = 1; cht.readID = 1;//make it as primary chain + cht.strand = 0; + ch_a[--occ] = cht; + // fprintf(stderr, "occ::%ld[M::%s::] x::%u, y::%u, cnt::%u, cov::%u\n", occ, __func__, + // cht.self_offset, cht.offset, cht.cnt, cht.readID); + } + for (k = ch_n-1; k >= 0; k--) { + if(!(ch_a[k].cnt&(0xffu))) continue; + ///end + cht.self_offset = ch_a[k].self_offset+1; + cht.offset = ch_a[k].offset+1; + cht.strand = 0; + cht.cnt = cht.readID = (ch_a[k].cnt&(0xffu)); + //make it as non-primary chain + if((ch_a[k].cnt&(0xffu)) < h_khit) cht.readID = cht.cnt + 1; + ch_a[--occ] = cht; + // fprintf(stderr, "occ::%ld[M::%s::] x::%u, y::%u, cnt::%u, cov::%u\n", occ, __func__, + // cht.self_offset, cht.offset, cht.cnt, cht.readID); + + if((ch_a[k].cnt&(0xffu)) > 1) { + ///start + cht.self_offset = ch_a[k].self_offset+1-(ch_a[k].cnt&(0xffu)); + cht.offset = ch_a[k].offset+1-(ch_a[k].cnt&(0xffu)); + cht.strand = 0; + cht.cnt = cht.readID = (ch_a[k].cnt&(0xffu)); + //make it as non-primary chain + if((ch_a[k].cnt&(0xffu)) < h_khit) cht.readID = cht.cnt + 1; + ch_a[--occ] = cht; + // fprintf(stderr, "occ::%ld[M::%s::] x::%u, y::%u, cnt::%u, cov::%u\n", occ, __func__, + // cht.self_offset, cht.offset, cht.cnt, cht.readID); + } + } + + if(prefix) { ///global or forward + cht.self_offset = qs; + cht.offset = ts; + cht.cnt = 1; cht.readID = 1;//make it as primary chain + cht.strand = 0; + ch_a[--occ] = cht; + // fprintf(stderr, "occ::%ld[M::%s::] x::%u, y::%u, cnt::%u, cov::%u\n", occ, __func__, + // cht.self_offset, cht.offset, cht.cnt, cht.readID); + } + assert(occ == 0); + ch_n = occ = ncn - cl->length; + uint64_t q[2], t[2]; + q[0] = q[1] = t[0] = t[1] = (uint64_t)-1; + if(prefix) { + q[0] = qs; t[0] = ts; + } + if(suffix) { + q[1] = qe; t[1] = te; + } + + // for (k = 0; k < ch_n; k++) { + // fprintf(stderr, "0) && (ch_a[k].self_offset==q[0]) && (ch_a[k].offset=t[0])) continue; + if(((k+1)0) { + if((ch_a[k].self_offset>ch_a[m-1].self_offset) && (ch_a[k].offset>ch_a[m-1].offset)) { + occ++; + } + if(ch_a[k].self_offset<=ch_a[m-1].self_offset) srt = 0; + } else { + occ++; + } + ch_a[m++] = ch_a[k]; + } + ch_n = m; + if(occ == ch_n) return ch_n;///already colinear + if(!srt) { + radix_sort_k_mer_hit_self(ch_a, ch_a + ch_n); + for (i = 1, j = 0; i <= ch_n; i++) { + if (i == ch_n || ch_a[i].self_offset != ch_a[j].self_offset) { + if(i - j > 1) radix_sort_k_mer_hit_off(ch_a+j, ch_a+i); + j = i; + } + } + } + + // for (k = 0; k < ch_n; k++) { + // fprintf(stderr, ">i::%ld[M::%s::] x::%u, y::%u, cnt::%u, cov::%u\n", k, __func__, + // ch_a[k].self_offset, ch_a[k].offset, ch_a[k].cnt, ch_a[k].readID); + // } + occ = ch_n; + ch_n = lchain_simple(ch_a+prefix, ch_n-prefix-suffix, ch_a+prefix, &(cl->chainDP), max_skip, max_iter); + ch_n += prefix + suffix; if(suffix) ch_a[ch_n-1] = ch_a[occ-1]; + // for (k = 0; k < ch_n; k++) { + // fprintf(stderr, "-i::%ld[M::%s::] x::%u, y::%u, cnt::%u, cov::%u\n", k, __func__, + // ch_a[k].self_offset, ch_a[k].offset, ch_a[k].cnt, ch_a[k].readID); + // } + + return ch_n; +} + +///[qs, qe) && [ts, te) +int64_t gen_win_chain(overlap_region *z, Candidates_list *cl, int64_t qs, int64_t qe, int64_t ts, int64_t te, +int64_t wl, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, +int64_t ql, int64_t tl, double e_rate, int64_t h_khit, int64_t mode) +{ + assert(mode < 3); + int64_t k, ws, wsk, we, wek, rcn = cl->length, ncn, occ = 0; asg16_v ez; uint32_t w = 1; + ws = qs; if(ws < z->x_pos_s) ws = z->x_pos_s; + we = qe-1; if(we > z->x_pos_e) we = z->x_pos_e; + wsk = get_win_id_by_s(z, ((ws/wl)*wl), wl, NULL); + assert((ws>=z->w_list.a[wsk].x_start) && (ws<=z->w_list.a[wsk].x_end)); + wek = get_win_id_by_e(z, ((we/wl)*wl), wl, NULL); + assert((we>=z->w_list.a[wek].x_start) && (we<=z->w_list.a[wek].x_end)); + ///global or forward + if(mode == 0 || mode == 1) push_khit(cl, qs, ts, 0, 0, &w); + //[ws, we] && [wsk, wek]; [qs, qe) && [ts, te) + for (k = wsk; k <= wek; k++) { + if(z->w_list.a[k].y_end == -1) continue; + if(!(z->w_list.a[k].clen)) { + gen_backtrace_adv_exz(&(z->w_list.a[k]), z, NULL, NULL, uref, qstr, tu->seq, exz, z->y_pos_strand, z->y_id); + } + ez.a = z->w_list.c.a + z->w_list.a[k].cidx; + ez.n = ez.m = z->w_list.a[k].clen; + occ += extract_exact_cigar(&ez, z->w_list.a[k].y_start, z->w_list.a[k].x_start, ts, te, qs, qe, cl, 10, wl, h_khit); + } + ///global or backward + if(mode == 0 || mode == 2) push_khit(cl, qe-1, te-1, 0, 0, &w); + // fprintf(stderr, "[M::%s::] rcn::%ld, cl->length::%lld\n", __func__, rcn, cl->length); + if(!occ) { + cl->length = rcn; return 0; + } + ncn = cl->length; cl->length = rcn; + k_mer_hit *ch_a = cl->list + rcn; int64_t ch_n0 = ncn - rcn, ch_n; + int64_t max_skip, max_iter, max_dis, quick_check; double chn_pen_gap, chn_pen_skip; + set_lchain_dp_op(0, h_khit, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check); + max_dis = MAX_SIN_L>>1; + // for (k = 0; k < ch_n0; k++) { + // assert(debug_k_mer_hit_retrive(&(ch_a[k]), hpc_g, rref, uref, qstr, tu, z->y_id, z->y_pos_strand)); + // } + ch_n = lchain_qdp_fix(ch_a, ch_n0, &(cl->chainDP), max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, + e_rate, ql, tl, 1, ((mode==0)||(mode==1))?1:0, ((mode==0)||(mode==2))?1:0); + // fprintf(stderr, "[M::%s::] ch_n0::%ld, ch_n::%ld, mode::%ld, ql::%ld, tl::%ld\n", + // __func__, ch_n0, ch_n, mode, qe-qs, te-ts); + for (k = occ = 0; k < ch_n; k++) { + ch_a[k] = ch_a[cl->chainDP.tmp[k]]; + if((ch_a[k].cnt&(0xffu))) occ++; + // assert(debug_k_mer_hit_retrive(&(ch_a[k]), hpc_g, rref, uref, qstr, tu, z->y_id, z->y_pos_strand)); + } + if(occ <= 0) return 0; + ch_n = gen_single_khit(cl, ch_n, h_khit, mode, qs, qe, ts, te, max_skip, max_iter); + return ch_n; +} + + +void rechain_aln(overlap_region *z, Candidates_list *cl, overlap_region *aux_o, int64_t aux_i, int64_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, +int64_t ql, int64_t tl, int64_t h_khit, int64_t rid) +{ + int64_t rcn = cl->length, ch_n, qs, qe, ts, te, mode, an0, an, todo; + k_mer_hit *ch_a; ul_ov_t idx; uint8_t q[2], t[2]; + ///[qs, qe) && [ts, te) + qs = aux_o->w_list.a[aux_i].x_start; qe = aux_o->w_list.a[aux_i].x_end+1; + ts = aux_o->w_list.a[aux_i].y_start; te = aux_o->w_list.a[aux_i].y_end+1; + mode = aux_o->w_list.a[aux_i].error_threshold; + ch_n = gen_win_chain(z, cl, qs, qe, ts, te, wl, uref, hpc_g, rref, qstr, tu, exz, ql, tl, e_rate, h_khit, mode); + ch_a = cl->list + rcn; + if(ch_n) { + idx.ts = idx.te = (uint32_t)-1; idx.qs = 0; idx.qe = ql; todo = 1; + if(mode == 0) {//global + idx.qn = 0; idx.tn = ch_n - 1; + idx.qs = ch_a[idx.qn].self_offset; + idx.ts = ch_a[idx.qn].offset; + idx.qe = ch_a[idx.tn].self_offset; + idx.te = ch_a[idx.tn].offset; + assert(ch_a[0].self_offset == qs && ch_a[0].offset == ts); + assert(ch_a[ch_n-1].self_offset == qe && ch_a[ch_n-1].offset == te); + if(ch_n <= 2) todo = 0; + } else if(mode == 1) {//forward ext + idx.qn = 0; idx.tn = ch_n; + idx.qs = ch_a[idx.qn].self_offset; + idx.ts = ch_a[idx.qn].offset; + idx.qe = ql; + assert(ch_a[0].self_offset == qs && ch_a[0].offset == ts); + if(ch_n <= 1) todo = 0; + } else if(mode == 2) {///backward ext + idx.qn = (uint32_t)-1; idx.tn = ch_n-1; + idx.qs = 0; + idx.qe = ch_a[idx.tn].self_offset; + idx.te = ch_a[idx.tn].offset; + assert(ch_a[ch_n-1].self_offset == qe && ch_a[ch_n-1].offset == te); + if(ch_n <= 1) todo = 0; + } + if(todo) { + an0 = aux_o->w_list.n; + ovlp_base_aln(z, ch_a, ch_n, &idx, wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, (uint64_t)-1); + an = aux_o->w_list.n; q[0] = q[1] = t[0] = t[1] = 0; todo = 0; + // fprintf(stderr, "[M::%s::] awn0::%ld, awn::%lu\n", __func__, an0, an); + ///old unaligned window could be replaced by the new aligned window + if((an == (an0 + 1)) && (!(is_ualn_win(aux_o->w_list.a[an-1])))) { + if(aux_o->w_list.a[aux_i].x_start == aux_o->w_list.a[an-1].x_start) q[0] = 1; + if(aux_o->w_list.a[aux_i].x_end == aux_o->w_list.a[an-1].x_end) q[1] = 1; + if(aux_o->w_list.a[aux_i].y_start == aux_o->w_list.a[an-1].y_start) t[0] = 1; + if(aux_o->w_list.a[aux_i].y_end == aux_o->w_list.a[an-1].y_end) t[1] = 1; + if((mode == 0) && q[0] && q[1] && t[0] && t[1]) todo = 1; + if((mode == 1) && q[0] && t[0]) todo = 1; + if((mode == 2) && q[1] && t[1]) todo = 1; + if(todo) { + aux_o->w_list.a[aux_i] = aux_o->w_list.a[an-1]; aux_o->w_list.n--; + } + } + // if(an > an0) {///should always > 0 as there are unmapped windows + // } + // aux_o->w_list.n = an0; + } + } + cl->length = rcn;///must reset!!!! +} + +void debug_overlap_region(overlap_region *au, char* qstr, UC_Read *tu, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref) +{ + int64_t wn = au->w_list.n, k; bit_extz_t ez; + for (k = 0; k < wn; k++) { + assert((k<=0)||((au->w_list.a[k].x_start>au->w_list.a[k-1].x_end) + &&(au->w_list.a[k].y_start>au->w_list.a[k-1].y_end))); + assert(au->w_list.a[k].x_end>au->w_list.a[k].x_start); + assert(au->w_list.a[k].y_end>au->w_list.a[k].y_start); + if(is_ualn_win(au->w_list.a[k]) || is_est_aln(au->w_list.a[k])) continue; + ez.cigar.a = au->w_list.c.a + au->w_list.a[k].cidx; + ez.cigar.n = ez.cigar.m = au->w_list.a[k].clen; + ez.ts = au->w_list.a[k].x_start; ez.te = au->w_list.a[k].x_end; + ez.ps = au->w_list.a[k].y_start; ez.pe = au->w_list.a[k].y_end; + ez.err = au->w_list.a[k].error; + ref_cigar_check(qstr, tu, uref, hpc_g, rref, au->y_id, au->y_pos_strand, &ez); + } +} + void cigar_gen_by_chain_adv(overlap_region *z, Candidates_list *cl, int64_t ch_idx, int64_t ch_n, -ul_ov_t *ov, int64_t on, uint64_t wl, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, -char* qstr, UC_Read *tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, uint64_t rid) +ul_ov_t *ov, int64_t on, uint64_t wl, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, +UC_Read *tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, uint64_t rid, int64_t h_khit) { if(on <= 0) return; - int64_t i, tl, id = z->y_id; + int64_t i, tl, id = z->y_id, m; k_mer_hit *ch_a = cl->list + ch_idx; - Chain_Data *dp = &(cl->chainDP); if(hpc_g) tl = hpc_len(*hpc_g, id); else if(uref) tl = uref->ug->u.a[id].len; else tl = Get_READ_LENGTH((*rref), id); @@ -14392,19 +14755,35 @@ char* qstr, UC_Read *tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, // ov->qs, ov->qe, ov->ts, ov->te, // (ov->qn!=((uint32_t)-1))?(int32_t)ov->qn:-1, (int32_t)ov->tn); assert((i<=0)||(ov[i].qs>ov[i-1].qe)); - ovlp_base_aln(z, dp, ch_a, ch_n, &(ov[i]), wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, rid); + ovlp_base_aln(z, ch_a, ch_n, &(ov[i]), wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, rid); } int64_t aux_n = aux_o->w_list.n; for (i = 0; i < aux_n; i++) { - // fprintf(stderr, "[aln::i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, clen::%u\n", i, - // aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start, - // aux_o->w_list.a[i].x_start, aux_o->w_list.a[i].x_end+1, - // aux_o->w_list.a[i].y_start, aux_o->w_list.a[i].y_end+1, - // aux_o->w_list.a[i].error, aux_o->w_list.a[i].clen); + if(!(is_ualn_win(aux_o->w_list.a[i]))) continue; + // if((aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start) <= FORCE_CNS_L) { + // fprintf(stderr, "[aln::i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, clen::%u, mode::%d\n", i, + // aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start, + // aux_o->w_list.a[i].x_start, aux_o->w_list.a[i].x_end+1, + // aux_o->w_list.a[i].y_start, aux_o->w_list.a[i].y_end+1, + // aux_o->w_list.a[i].error, aux_o->w_list.a[i].clen, aux_o->w_list.a[i].error_threshold); + // } + //will overwrite ch_a; does not matter + rechain_aln(z, cl, aux_o, i, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, tl, h_khit, rid); + } + if(((int64_t)aux_o->w_list.n) > aux_n) { + for (i = m = 0; i < ((int64_t)aux_o->w_list.n); i++) { + if((i < aux_n) && (is_ualn_win(aux_o->w_list.a[i]))) continue; + aux_o->w_list.a[m++] = aux_o->w_list.a[i]; + } + aux_o->w_list.n = m; + radix_sort_window_list_xs_srt(aux_o->w_list.a, aux_o->w_list.a+aux_o->w_list.n); } - ch_a = cl->list + ch_idx; //update + // debug_overlap_region(aux_o, qstr, tu, uref, hpc_g, rref); + + + // ch_a = cl->list + ch_idx; //update // for (i = 0; i < wn; i++) z->w_list.a[i].clen = 0;///clean cigar // if(on > 1) { // fprintf(stderr, "[M::%s::] rid::%lu, on::%ld\n", __func__, rid, on); @@ -14434,7 +14813,7 @@ int64_t max_lgap) if(ch_n) { ///ch_a = cl->list + cl->length; // cigar_gen_by_chain(z, &(cl->chainDP), ch_a, ch_n, aln->a+l, k-l, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); - cigar_gen_by_chain_adv(z, cl, cl->length, ch_n, aln->a+l, k-l, wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, rid); + cigar_gen_by_chain_adv(z, cl, cl->length, ch_n, aln->a+l, k-l, wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, rid, khit); // m = fusion_coordinates(z, ch_a, ch_n, aln->a+l, k-l); } // m += cigar_gen_cns(z, cl, aln->a+l, k-l, i, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid, aln->a+m); @@ -14541,7 +14920,7 @@ int64_t max_lgap, double sgap_rate, int64_t sgap)///[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; + overlap_region *z; uint64_t ws, we; int64_t r_y[2], r_err, p_y[2], p_err, rxl, ryl, pxl; ///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; @@ -14566,13 +14945,17 @@ int64_t max_lgap, double sgap_rate, int64_t sgap)///[s, e) 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(!ovlp_win_check(z, z->align_length, wid, max_lgap, sgap_rate, sgap)) continue;///not co-linear + rxl = z->w_list.a[wid].x_end+1-z->w_list.a[wid].x_start; + ryl = r_y[1]-r_y[0]; + if(rxl > 1 && ryl > 1) continue;///length of window should be longer than 1 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(!ovlp_win_check(z, z->align_length, wid, max_lgap, sgap_rate, sgap)) continue;///not co-linear - + if(!ovlp_win_check(z, z->align_length, wid, max_lgap, sgap_rate, sgap)) break;///not co-linear + pxl = z->w_list.a[wid].x_end+1-z->w_list.a[wid].x_start; + if(pxl != rxl) break; 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) { diff --git a/Correct.h b/Correct.h index af73081..cf076d8 100644 --- a/Correct.h +++ b/Correct.h @@ -1331,6 +1331,9 @@ All_reads *rref, UC_Read* tu, asg64_v* idx, asg64_v *b0, asg64_v *b1, int64_t ql kv_ul_ov_t *aln, uint64_t rid, int64_t max_lgap, double sgap_rate); #define copy_asg_arr(des, src) ((des).a = (src).a, (des).n = (src).n, (des).m = (src).m) +#define is_ualn_win(a) (((a).error==INT16_MAX)&&((a).clen==0)&&((a).extra_end<0)) +#define is_exact_aln(a) (((a).error0)) +#define is_est_aln(a) (((a).errorcnt&(0xffu); sc = q_span < dg? q_span : dg; sc = normal_w(sc, ((int32_t)(ai->cnt>>8))); - if (dd || dg > q_span) { + if (dd || (dg > q_span && dg > 0)) { + double lin_pen, a_pen; + lin_pen = (chn_pen_gap*(double)dd); + a_pen = ((double)(sc))*((((double)dd)/((double)dg))/bw_rate); + if(lin_pen > a_pen) lin_pen = a_pen; + lin_pen += (chn_pen_skip*(double)dg); + sc -= (int32_t)lin_pen; + } + return sc; +} + +inline int32_t comput_sc_ff(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, double chn_pen_gap, double chn_pen_skip, int64_t sl, int64_t ol) +{ + ///ai is the suffix of aj + int32_t dq, dr, dd, dg, q_span, sc; + dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset); + if(dq < 0) return INT32_MIN; + dr = (int64_t)(ai->offset) - (int64_t)(aj->offset); + if(dr < 0) return INT32_MIN; + dd = dr > dq? dr - dq : dq - dr;//gap + // if((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, sl, ol))) return INT32_MIN; + dg = dr < dq? dr : dq;//len + q_span = ai->cnt&(0xffu); + sc = q_span < dg? q_span : dg; + sc = normal_w(sc, ((int32_t)(ai->cnt>>8))); + if (dd || (dg > q_span && dg > 0)) { double lin_pen, a_pen; lin_pen = (chn_pen_gap*(double)dd); a_pen = ((double)(sc))*((((double)dd)/((double)dg))/bw_rate); @@ -1816,6 +1841,178 @@ uint64_t lchain_qdp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, o } + + +#define rev_khit(an, xl, yl) do { \ + (an).self_offset = (xl)-1-((an).self_offset+1-((an).cnt&((uint32_t)(0xffu)))); \ + (an).offset = (yl)-1-((an).offset+1-((an).cnt&((uint32_t)(0xffu))));\ + } while (0) + +///it is unpossibale that (left_fix == 0 && right_fix == 0) +uint64_t lchain_qdp_fix(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_skip, + int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip, + double bw_rate, int64_t xl, int64_t yl, int64_t quick_check, + int64_t left_fix, int64_t right_fix) +{ + int64_t *p, *t, max_f, n_skip, st, max_j, end_j, sc, msc, msc_i, bw, max_ii, ovl, movl; + int32_t *f, max, tmp; int64_t i, j, ret, cL = 0, must_p = 1, is_reorder = 0; k_mer_hit z; + resize_Chain_Data(dp, a_n, NULL); + t = dp->tmp; f = dp->score; p = dp->pre; + bw = ((xl < yl)?xl:yl); bw *= bw_rate; + msc = msc_i = -1; movl = INT32_MAX; + + if(quick_check) { + ret = lchain_qcheck(a, a_n, dp, bw_rate); + if (ret > 0) { + a_n = ret; msc_i = a_n-1; msc = f[msc_i]; + goto skip_ldp; + } + } + + memset(t, 0, (a_n*sizeof((*t)))); + if((right_fix) && (!left_fix)) { + // fprintf(stderr, "\n[M::%s::] a_n::%ld\n", __func__, a_n); + // for (i = 0; i < a_n; ++i) { + // fprintf(stderr, "+[M::%s::] x::[%u, %u), y::[%u, %u)\n", __func__, + // a[i].self_offset+1-(a[i].cnt&((uint32_t)(0xffu))), a[i].self_offset+1, + // a[i].offset+1-(a[i].cnt&((uint32_t)(0xffu))), a[i].offset+1); + // } + n_skip = a_n>>1; + for (i=0; i max_iter) st = i-max_iter; + // if(a[0].self_offset == 171728) { + // fprintf(stderr, "i::%ld[M::%s::] x::[%u, %u), y::[%u, %u)\n", i, __func__, + // a[i].self_offset+1-(a[i].cnt&((uint32_t)(0xffu))), a[i].self_offset+1, + // a[i].offset+1-(a[i].cnt&((uint32_t)(0xffu))), a[i].offset+1); + // } + for (j = i - 1; j >= 0; --j) { + sc = comput_sc_ff(&a[i], &a[j], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl); + // if(a[0].self_offset == 171728) { + // fprintf(stderr, "j::%ld[M::%s::sc->%ld] x::[%u, %u), y::[%u, %u)\n", j, __func__, sc, + // a[j].self_offset+1-(a[j].cnt&((uint32_t)(0xffu))), a[j].self_offset+1, + // a[j].offset+1-(a[j].cnt&((uint32_t)(0xffu))), a[j].offset+1); + // } + if (sc == INT32_MIN) continue; + sc += f[j]; + if (sc > max_f) { + max_f = sc, max_j = j; + if (n_skip > 0) --n_skip; + } else if (t[j] == (int32_t)i) { + if ((++n_skip) > max_skip) { + if((max_j != -1) || (must_p == 0)) break; + } + } + if (p[j] >= 0) t[p[j]] = i; + ///put it here will allow at least one prefix no matter max_dis + ///this is special for gap filling, not for chaining + if (a[i].self_offset > (max_dis + a[j].self_offset)) { + if((max_j != -1)) break; + } + if (j < st) { + if((max_j != -1) || (must_p == 0)) break; + } + } + end_j = j; + + if (max_ii < 0 || ((int64_t)a[i].self_offset) - ((int64_t)a[max_ii].self_offset) > max_dis) { + max = INT32_MIN; max_ii = -1; + for (j = i - 1; (j >= st) && ((((int64_t)a[i].self_offset)-((int64_t)a[j].self_offset))<=max_dis); --j) { + if (max < f[j]) { + max = f[j], max_ii = j; + } + } + } + + if (max_ii >= 0 && max_ii < end_j) {///just have a try with a[i]<->a[max_ii] + tmp = comput_sc_ff(&a[i], &a[max_ii], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl); + if (tmp != INT32_MIN && max_f < tmp + f[max_ii]) + max_f = tmp + f[max_ii], max_j = max_ii; + } + if(max_j == -1) { + f[i] = 0; p[i] = max_j; + } else { + f[i] = max_f; p[i] = max_j; + } + + if ((max_ii < 0) || (((((int64_t)a[i].self_offset)-((int64_t)a[max_ii].self_offset))<=max_dis) && (f[max_ii]= msc) { + ovl = get_chainLen(a[i].self_offset, a[i].self_offset, xl, a[i].offset, a[i].offset, yl); + if(f[i] > msc || ovl < movl) { + msc = f[i]; msc_i = i; movl = ovl; + } + } + } + + skip_ldp: + if(right_fix && left_fix) msc_i = a_n-1; + ///a[] has been sorted by self_offset + i = msc_i; cL = 0; + while (i >= 0) { + t[cL++] = i; msc_i = i; i = p[i]; + } + // if((right_fix) && (!left_fix)) { + // fprintf(stderr, "[M::%s::] cL::%ld\n", __func__, cL); + // } + + if(is_reorder) { + n_skip = a_n>>1; + for (i=0; i[M::%s::] t[%ld]::%ld\n", __func__, i, t[i]); + } + // for (i = 0; i < a_n; ++i) { + // fprintf(stderr, ">[M::%s::] x::[%u, %u), y::[%u, %u)\n", __func__, + // a[i].self_offset+1-(a[i].cnt&((uint32_t)(0xffu))), a[i].self_offset+1, + // a[i].offset+1-(a[i].cnt&((uint32_t)(0xffu))), a[i].offset+1); + // } + } else { + n_skip = cL>>1; + for (i = 0; i < n_skip; i++) { + msc_i = t[i]; t[i] = t[cL-i-1]; t[cL-i-1] = msc_i; + } + } + // if(cL != a_n) { + // fprintf(stderr, "\n[M::%s::] a_n::%ld, left_fix::%ld, right_fix::%ld\n", __func__, a_n, left_fix, right_fix); + // for (i = 0; i < a_n; ++i) { + // fprintf(stderr, "+[M::%s::] x::[%u, %u), y::[%u, %u)\n", __func__, + // a[i].self_offset+1-(a[i].cnt&((uint32_t)(0xffu))), a[i].self_offset+1, + // a[i].offset+1-(a[i].cnt&((uint32_t)(0xffu))), a[i].offset+1); + // } + // for (i = 0; i < cL; i++) { + // fprintf(stderr, ">[M::%s::] t[%ld]::%ld\n", __func__, i, t[i]); + // } + // } + return cL; +} + + uint64_t lchain_refine(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, int64_t max_skip, int64_t max_iter, int64_t max_dis, int64_t long_gap) { @@ -1905,6 +2102,70 @@ uint64_t lchain_refine(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp return cL; } + +uint64_t lchain_simple(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, + int64_t max_skip, int64_t max_iter) +{ + if(a_n <= 0) return 0; + int64_t *p, *t, max_f, n_skip, st, max_j, sc, msc, msc_i; + int32_t *f; int64_t i, j, cL = 0; + resize_Chain_Data(dp, a_n, NULL); + t = dp->tmp; f = dp->score; p = dp->pre; msc = msc_i = -1; + + for (i=1, f[0]=a[0].cnt, p[0]=-1, msc_i=a_n-1; i a[j].self_offset)&&(a[i].offset > a[j].offset)) { + p[i] = j; f[i] = f[j]+a[i].cnt; + } else { + break; + } + } + if(i >= a_n) goto ss_kip; + + memset(t, 0, (a_n*sizeof((*t)))); + f[0]=a[0].cnt; p[0]=-1; msc = f[0]; msc_i = 0; + + for (i = 1, st = 0; i < a_n; ++i) { + max_f = INT32_MIN; n_skip = 0; max_j = -1; + if ((i-st) > max_iter) st = i-max_iter; + ///[st, i-2] + for (j=i-1; j >= st; --j) { + if((a[i].self_offset > a[j].self_offset)&&(a[i].offset > a[j].offset)) { + sc = f[j]+a[i].cnt; + if (sc > max_f) { + max_f = sc, max_j = j; + if (n_skip > 0) --n_skip; + } else if (t[j] == (int32_t)i) { + if (++n_skip > max_skip) + break; + } + if (p[j] >= 0) t[p[j]] = i; + } + } + f[i] = max_f; p[i] = max_j; + if(f[i] > msc) { + msc = f[i]; msc_i = i; + } + } + + ss_kip: + ///a[] has been sorted by self_offset + i = msc_i; + cL = 0; + while (i >= 0) { + t[cL++] = i; i = p[i]; + } + + n_skip = cL>>1; + for (i = 0; i < n_skip; i++) { + msc_i = t[i]; t[i] = t[cL-i-1]; t[cL-i-1] = msc_i; + } + if(des) { + for (i = 0; i < cL; i++) des[i] = a[t[i]]; + } + return cL; +} + inline int64_t hit_long_gap(k_mer_hit *a, k_mer_hit *b, int64_t max_lgap, double small_bw_rate, int64_t min_small_bw) { int64_t dq, dr, dd, dm; diff --git a/Hash_Table.h b/Hash_Table.h index 1080eab..928f860 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -214,4 +214,11 @@ uint64_t lchain_qdp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, o int ovlp_chain_qgen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int64_t yl, int64_t apend_be, k_mer_hit* hit, int64_t n_hit); uint64_t lchain_refine(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, int64_t max_skip, int64_t max_iter, int64_t max_dis, int64_t long_gap); +uint64_t lchain_qdp_fix(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_skip, + int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip, + double bw_rate, int64_t xl, int64_t yl, int64_t quick_check, + int64_t left_fix, int64_t right_fix); +uint64_t lchain_simple(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, + int64_t max_skip, int64_t max_iter); + #endif