diff --git a/Correct.cpp b/Correct.cpp index 8c5f75a..9ce9c96 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -13,6 +13,7 @@ #include "kalloc.h" #include "htab.h" #include "Overlaps.h" +#include "inter.h" #define A_L 16 #define ext_w 6 @@ -14732,6 +14733,36 @@ void debug_overlap_region(overlap_region *au, char* qstr, UC_Read *tu, const ul_ } } +void update_overlap_region(overlap_region *des, overlap_region *src, int64_t xl, int64_t yl) +{ + kv_resize(uint16_t, des->w_list.c, src->w_list.c.n); + des->w_list.c.n = src->w_list.c.n; + memcpy(des->w_list.c.a, src->w_list.c.a, src->w_list.c.n*(sizeof((*(src->w_list.c.a))))); + + kv_resize(window_list, des->w_list, src->w_list.n); + des->w_list.n = src->w_list.n; + memcpy(des->w_list.a, src->w_list.a, src->w_list.n*(sizeof((*(src->w_list.a))))); + + if(src->w_list.n) { + des->x_pos_s = src->w_list.a[0].x_start; des->x_pos_e = src->w_list.a[src->w_list.n-1].x_end; + des->y_pos_s = src->w_list.a[0].y_start; des->y_pos_e = src->w_list.a[src->w_list.n-1].y_end; + } + + int64_t xr, yr; + if(des->x_pos_s <= des->y_pos_s) { + des->y_pos_s -= des->x_pos_s; des->x_pos_s = 0; + } else { + des->x_pos_s -= des->y_pos_s; des->y_pos_s = 0; + } + + xr = xl-des->x_pos_e-1; yr = yl-des->y_pos_e-1; + if(xr <= yr) { + des->x_pos_e = xl-1; des->y_pos_e += xr; + } else { + des->y_pos_e = yl-1; des->x_pos_e += yr; + } +} + 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, int64_t h_khit) @@ -14780,6 +14811,9 @@ UC_Read *tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, radix_sort_window_list_xs_srt(aux_o->w_list.a, aux_o->w_list.a+aux_o->w_list.n); } + ///update z by aux_o + update_overlap_region(z, aux_o, ql, tl); + // debug_overlap_region(aux_o, qstr, tu, uref, hpc_g, rref); @@ -14797,29 +14831,529 @@ UC_Read *tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, // } } +#define set_bit_extz_t(x, z, id) do {\ + (x).cigar.a = (z).w_list.c.a+(z).w_list.a[(id)].cidx;\ + (x).cigar.n = (x).cigar.m = (z).w_list.a[(id)].clen;\ + (x).ts = (z).w_list.a[(id)].x_start;\ + (x).te = (z).w_list.a[(id)].x_end;\ + (x).ps = (z).w_list.a[(id)].y_start;\ + (x).pe = (z).w_list.a[(id)].y_end;\ + (x).err = (z).w_list.a[(id)].error;\ + } while (0) + +#define gen_err_unaligned(xl, yl) (((xl)<=FORCE_SIN_L)?(MAX((xl), (yl))):MAX((MIN((xl), (yl))), ((xl*0.51)+1))) + +#define ovlp_id(x) ((x).tn) +#define ovlp_min_wid(x) ((x).ts) +#define ovlp_max_wid(x) ((x).te) +#define ovlp_cur_wid(x) ((x).qn) +#define ovlp_cur_xoff(x) ((x).qs) +#define ovlp_cur_coff(x) ((x).qe) + + +int64_t retrieve_cigar_err(bit_extz_t *ez, int64_t s, int64_t e, int64_t *xk, int64_t *ck) +{ + if(!ez->cigar.n) return 0; + int64_t cn = ez->cigar.n, op, err = 0; int64_t ws, we, os, oe, ovlp; + if(((*ck) < 0) || ((*ck) > cn)) {//(*ck) == cn is allowed + (*ck) = 0; (*xk) = ez->ts; + } + + while ((*ck) > 0 && (*xk) > s) { + --(*ck); + op = ez->cigar.a[(*ck)]>>14; + if(op!=2) (*xk) -= (ez->cigar.a[(*ck)]&(0x3fff)); + } + + //some cigar will span s or e + while ((*ck) < cn && (*xk) < e) {//[s, e) + ws = (*xk); + op = ez->cigar.a[(*ck)]>>14; + if(op!=2) (*xk) += (ez->cigar.a[(*ck)]&(0x3fff)); + we = (*xk); + os = MAX(s, ws); oe = MIN(e, we); + ovlp = ((oe>os)? (oe-os):0); + if((op==2) && (ws>=s) && (wscigar.a[(*ck)]&(0x3fff)); + } + (*ck)++; + if((!ovlp) || (!op)) continue; + err += ovlp; + } + return err; +} +///[s, e) +int64_t extract_sub_cigar_err(overlap_region *z, int64_t s, int64_t e, ul_ov_t *p) +{ + int64_t wk = ovlp_cur_wid(*p), xk = ovlp_cur_xoff(*p), ck = ovlp_cur_coff(*p); + int64_t min_w = ovlp_min_wid(*p), max_w = ovlp_max_wid(*p);//[min_w, max_w] + bit_extz_t ez; window_list *m; + int64_t ws, we, os, oe, ovlp, err = 0, xl, yl, werr, tot = e - s; + if(wk < min_w || wk > max_w) wk = min_w; + for (; wk >= min_w && z->w_list.a[wk].x_start > s; wk--); + if(wk < min_w || wk > max_w) return -1; + for (; wk <= max_w && z->w_list.a[wk].x_end < s; wk++); + if(wk < min_w || wk > max_w) return -1; + //s >= w_list.a[wk].x_start && s <= w_list.a[wk].x_end + if(wk != ovlp_cur_wid(*p)) {//xk is global, while ck is local + xk = z->w_list.a[wk].x_start; ck = 0; + } + // fprintf(stderr, "[M::%s] wk::%ld, ck::%ld, xk::%ld\n", __func__, wk, ck, xk); + // fprintf(stderr, "+[M::%s] wk::%ld, ck::%ld, xk::%ld, w::[%d, %d), bound::[%ld, %ld)\n", + // __func__, wk, ck, xk, z->w_list.a[wk].x_start, z->w_list.a[wk].x_end+1, s, e); + + while(wk <= max_w && z->w_list.a[wk].x_start < e) {///[s, e) + m = &(z->w_list.a[wk]); + ws = m->x_start; we = m->x_end+1; + os = MAX(s, ws); oe = MIN(e, we); + ovlp = ((oe>os)? (oe-os):0); + // fprintf(stderr, "-[M::%s] ovlp::%ld, wk::%ld, ck::%ld, xk::%ld, w::[%ld, %ld), bound::[%ld, %ld)\n", + // __func__, ovlp, wk, ck, xk, ws, we, s, e); + if(ovlp) { + xl = m->x_end+1-m->x_start; + yl = m->y_end+1-m->y_start; + if((is_ualn_win((*m))) || (is_est_aln((*m)))) { + if(is_ualn_win((*m))) { //unmapped + werr = gen_err_unaligned(xl, yl); + } else { + werr = m->error;//shared window + } + if(ovlp < xl) { + werr = (((double)ovlp)/((double)xl))*((double)werr); + } + //skip the whole window + err += werr; xk = m->x_end+1; ck = m->clen; + } else { + if(ovlp == xl) { + //skip the whole window + err += m->error; xk = m->x_end+1; ck = m->clen; + } else { + set_bit_extz_t(ez, (*z), wk); + err += retrieve_cigar_err(&ez, os, oe, &xk, &ck); + } + } + } + tot -= ovlp; + if(xk >= e) break;//[min_w, max_w] && [s, e) + wk++; if(wk > max_w) break; + xk = z->w_list.a[wk].x_start; ck = 0;//reset + } + assert(!tot); + ovlp_cur_wid(*p) = wk; ovlp_cur_xoff(*p) = xk; ovlp_cur_coff(*p) = ck; + return err; +} + +#define bst_ov(x) ((x).misBase) +#define ov_dif(x) ((x).cov) +#define ov_id(x) ((x).overlapID) +#define ov_xoff(x) ((x).site) +#define var_id(x) ((x).overlapSite) + +#define var_s(x) ((x).site) +#define var_l(x) ((x).overlap_num) +#define var_occ(x) ((x).occ_0) +#define var_min_dif(x) ((x).score) +#define var_min_ovid(x) ((x).id) +#define var_h_idx(x) ((x).occ_1) + + +uint64_t query_gen_gov_idx(asg64_v *ovidx, uint64_t v, uint64_t w) +{ + uint64_t m, s, e; + if(v > w) { + m = v; v = w; w = m; + } + s = ovidx->a[v]>>32; e = s + (uint32_t)ovidx->a[v]; + for (m = s; m < e; m++) { + if(ovidx->a[m] == w) return 1; + } + return 0; +} + +///[s, e) +uint64_t gen_region_phase(overlap_region* ol, uint64_t *id_a, uint64_t id_n, uint64_t s, uint64_t e, uint64_t dp, ul_ov_t *c_idx, uint64_t *buf, asg64_v *ovidx) +{ + if(!id_n) return id_n; + uint64_t k, m, mn, q[2], buf_n, rm_n, i; int64_t err, msc, msc_k, msc_n; + overlap_region *z; ul_ov_t *p; + for (k = buf_n = rm_n = 0; k < id_n; k++) { + p = &(c_idx[id_a[k]]); + q[0] = ol[ovlp_id(*p)].w_list.a[ovlp_min_wid(*p)].x_start; + q[1] = ol[ovlp_id(*p)].w_list.a[ovlp_max_wid(*p)].x_end+1; + if(q[0]<=s && q[1]>=e) { + buf[buf_n++] = id_a[k]; + } + if(q[1] < e) rm_n++; + } + assert(buf_n == dp);//not right + + if(buf_n > 0) { + for (k = 0, msc = INT32_MAX, msc_k = -1, msc_n = 0; k < buf_n; k++) { + p = &(c_idx[(uint32_t)buf[k]]); z = &(ol[ovlp_id(*p)]); + // fprintf(stderr, "+++[M::%s::utg%.6dl] wid::%u, xoff::%u, coff::%u\n", __func__, + // (int32_t)ol[ovlp_id(*p)].y_id+1, ovlp_cur_wid(*p), ovlp_cur_xoff(*p), ovlp_cur_coff(*p)); + err = extract_sub_cigar_err(z, s, e, p); + // fprintf(stderr, "---[M::%s::utg%.6dl] wid::%u, xoff::%u, coff::%u, err::%ld\n", __func__, + // (int32_t)ol[ovlp_id(*p)].y_id+1, ovlp_cur_wid(*p), ovlp_cur_xoff(*p), ovlp_cur_coff(*p), err); + assert(err >= 0); + if(err < msc) { + msc = err; msc_k = k; msc_n = 1; + } else if(err == msc) { + msc_n++; + } + buf[k] |= (((uint64_t)err)<<32); + } + + if(msc_n == 1) { + p = &(c_idx[(uint32_t)buf[msc_k]]); + z = &(ol[ovlp_id(*p)]); mn = 1; + if(msc_k != 0) { + m = buf[msc_k]; + buf[msc_k] = buf[0]; + buf[0] = m; + } + } else { + for (k = mn = 0; k < buf_n && (int64_t)mn < msc_n; k++) { + p = &(c_idx[(uint32_t)buf[k]]); + z = &(ol[ovlp_id(*p)]); + if((buf[k]>>32) == (uint64_t)msc) { + if(mn != k) { + m = buf[k]; + buf[k] = buf[mn]; + buf[mn] = m; + } + mn++; + } + } + } + // fprintf(stderr, "[M::%s] buf_n::%ld, msc_n::%ld, mn::%lu\n", __func__, buf_n, msc_n, mn); + for (k = 0; k < buf_n; k++) { + buf[k] = ovlp_id((c_idx[(uint32_t)buf[k]])); + // if(k < mn) fprintf(stderr, "d::bst::[M::%s::utg%.6dl]\n", __func__, (int32_t)ol[buf[k]].y_id+1); + } + for (k = mn; k < buf_n; k++) { + z = &(ol[buf[k]]); + z->align_length -= e - s; + for (i = 0; i < mn; i++) { + if(query_gen_gov_idx(ovidx, buf[k], buf[i])) break; + } + if(i < mn) { + z->align_length += e - s; + // fprintf(stderr, "i::bst::[M::%s::utg%.6dl]\n", __func__, (int32_t)ol[buf[k]].y_id+1); + } + } + } + + + if(rm_n) { + for (k = m = 0; k < id_n; k++) { + p = &(c_idx[id_a[k]]); + q[1] = ol[ovlp_id(*p)].w_list.a[ovlp_max_wid(*p)].x_end+1; + if(q[1] < e) continue; + id_a[m++] = id_a[k]; + } + id_n = m; + } + return id_n; +} + + +int64_t infer_rovlp(ul_ov_t *li, ul_ov_t *lj, uc_block_t *bi, uc_block_t *bj, All_reads *ridx, ma_ug_t *ug) +{ + int64_t in, is, ie, irev, iqs, iqe, jn, js, je, jrev, jqs, jqe, ir, jr, ts, te, max_s, min_e, s_shift, e_shift; + + if(li) { + in = ug?ug->u.a[li->tn].len:Get_READ_LENGTH(R_INF, li->tn); + is = li->ts; ie = li->te; irev = li->rev; iqs = li->qs; iqe = li->qe; + } else if(bi) { + in = ug?ug->u.a[bi->hid].len:Get_READ_LENGTH(R_INF, bi->hid); + is = bi->ts; ie = bi->te; irev = bi->rev; iqs = bi->qs; iqe = bi->qe; + } else { + return 0; + } + + if(lj) { + jn = ug?ug->u.a[lj->tn].len:Get_READ_LENGTH(R_INF, lj->tn); + js = lj->ts; je = lj->te; jrev = lj->rev; jqs = lj->qs; jqe = lj->qe; + } else if(bj) { + jn = ug?ug->u.a[bj->hid].len:Get_READ_LENGTH(R_INF, bj->hid); + js = bj->ts; je = bj->te; jrev = bj->rev; jqs = bj->qs; jqe = bj->qe; + } else { + return 0; + } + + max_s = MAX(iqs, jqs); min_e = MIN(iqe, jqe); + if(min_e <= max_s) return 0; + s_shift = get_offset_adjust(max_s - iqs, iqe-iqs, ie-is); + e_shift = get_offset_adjust(iqe - min_e, iqe-iqs, ie-is); + if(irev) { + ts = s_shift; s_shift = e_shift; e_shift = ts; + } + is += s_shift; ie-= e_shift; + + // if(li && lj && li->tn == 324 && lj->tn == 319 && li->qs == 63841) { + // fprintf(stderr, "+++in:%ld, is:%ld, ie:%ld, irev:%ld, jn:%ld, js:%ld, je:%ld, jrev:%ld\n", in, is, ie, irev, jn, js, je, jrev); + // } + + s_shift = get_offset_adjust(max_s - jqs, jqe-jqs, je-js); + e_shift = get_offset_adjust(jqe - min_e, jqe-jqs, je-js); + if(jrev) { + ts = s_shift; s_shift = e_shift; e_shift = ts; + } + js += s_shift; je-= e_shift; + + if(irev) { + ts = in - ie; te = in - is; + is = ts; ie = te; + } + + if(jrev) { + ts = jn - je; te = jn - js; + js = ts; je = te; + } + + // if(li && lj && li->tn == 324 && lj->tn == 319 && li->qs == 63841) { + // fprintf(stderr, "---in:%ld, is:%ld, ie:%ld, irev:%ld, jn:%ld, js:%ld, je:%ld, jrev:%ld\n", in, is, ie, irev, jn, js, je, jrev); + // } + + if(is <= js) { + js -= is; is = 0; + } else { + is -= js; js = 0; + } + + ir = in - ie; jr = jn - je; + + if(ir <= jr){ + ie = in; je += ir; + } + else { + je = jn; ie += jr; + } + + ir = ie - is; jr = je - js; + return MAX(ir, jr); +} + +void convert_ul_ov_t(ul_ov_t *des, overlap_region *src, const ul_idx_t *uref) +{ + des->qn = (uint32_t)-1; des->qs = src->x_pos_s; des->qe = src->x_pos_e+1; + des->tn = src->y_id; des->el = 1; des->rev = src->y_pos_strand; + des->sec = src->non_homopolymer_errors; + if(des->rev) { + des->ts = uref->ug->u.a[des->tn].len - (src->y_pos_e+1); + des->te = uref->ug->u.a[des->tn].len - src->y_pos_s; + } else { + des->ts = src->y_pos_s; + des->te = src->y_pos_e+1; + } +} + +uint64_t check_connect_ug(const ul_idx_t *uref, uint32_t v, uint32_t w, int64_t bw, double diff_ec_ul, int64_t dq) +{ + const asg_t *g = uref?uref->ug->g:NULL; int64_t dt = -1; + uint32_t nv = asg_arc_n(g, v), i; asg_arc_t *av = asg_arc_a(g, v); + for (i = 0; i < nv; i++) { + if(av[i].del || av[i].v != w) continue; + dt = av[i].ol; + break; + } + if(dt < 0) return 0; + int64_t diff = (dq>dt? dq-dt:dt-dq), mm = MAX(dq, dt); + mm *= diff_ec_ul; if(mm < bw) mm = bw; + if(diff <= mm) return 1; + return 0; +} + +uint64_t check_connect_rg(const ul_idx_t *uref, const ug_opt_t *uopt, uint32_t uv, uint32_t uw, int64_t bw, double diff_ec_ul, int64_t dq) +{ + int64_t dt = -1; + if(uref->ug->u.a[uv>>1].circ || uref->ug->u.a[uw>>1].circ) return 0; + uint32_t rv = (uref->ug->u.a[uv>>1].a[(uv&1)?(0):(uref->ug->u.a[uv>>1].n-1)]>>32)^(uv&1); + uint32_t rw = (uref->ug->u.a[uw>>1].a[(uw&1)?(uref->ug->u.a[uw>>1].n-1):(0)]>>32)^(uw&1); + ma_hit_t_alloc* src = uopt->sources; + int64_t min_ovlp = uopt->min_ovlp; + int64_t max_hang = uopt->max_hang; + uint64_t z, qn, tn, x = rv>>1; int32_t r = 1; asg_arc_t e; + for (z = 0; z < src[x].length; z++) { + qn = Get_qn(src[x].buffer[z]); tn = Get_tn(src[x].buffer[z]); + if(tn != (rw>>1)) continue; + r = ma_hit2arc(&(src[x].buffer[z]), Get_READ_LENGTH(R_INF, qn), Get_READ_LENGTH(R_INF, tn), max_hang, asm_opt.max_hang_rate, min_ovlp, &e); + if(r < 0) continue; + if((e.ul>>32) != rv || e.v != rw) continue; + dt = e.ol; + break; + } + if(dt < 0) return 0; + int64_t diff = (dq>dt? dq-dt:dt-dq), mm = MAX(dq, dt); + mm *= diff_ec_ul; if(mm < bw) mm = bw; + if(diff <= mm) return 1; + return 0; +} + +uint32_t govlp_check(const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, ul_ov_t *li, ul_ov_t *lj) +{ + int64_t qo = infer_rovlp(li, lj, NULL, NULL, /**ridx**/NULL, uref->ug); ///overlap length in query (UL read) + + // fprintf(stderr, "+++[M::%s::utg%.6dl->utg%.6dl] qo::%ld\n", __func__, (int32_t)li->tn+1, (int32_t)lj->tn+1, qo); + if(check_connect_ug(uref, ((li->tn<<1)|li->rev)^1, ((lj->tn<<1)|lj->rev)^1, bw, diff_ec_ul, qo)) return 1; + // fprintf(stderr, "[M::%s::] check_connect_ug fail\n", __func__); + if(check_connect_rg(uref, uopt, ((li->tn<<1)|li->rev)^1, ((lj->tn<<1)|lj->rev)^1, bw, diff_ec_ul, qo)) return 1; + // fprintf(stderr, "[M::%s::] check_connect_rg fail\n", __func__); + return 0; +} + +void gen_gov_idx(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, asg64_v* idx) +{ + int64_t on = ol->length, k, i; uint64_t os, oe, ovlp; ul_ov_t p, q, *li, *lj; + kv_resize(uint64_t, *idx, (uint64_t)on); memset(idx->a, 0, sizeof(*(idx->a))*on); + for (k = 0, idx->n = on; k < on; k++) { + convert_ul_ov_t(&p, &(ol->list[k]), uref); p.qn = k; + idx->a[k] = idx->n; idx->a[k] <<= 32; + for (i = on - 1; i >= 0 && i > k && ol->list[i].x_pos_e >= ol->list[k].x_pos_s; i--) { + // if(k >= i) continue; + convert_ul_ov_t(&q, &(ol->list[i]), uref); q.qn = i; + + if(p.qe > q.qe) li = &p, lj = &q; + else if(p.qe == q.qe && p.qs >= q.qs) li = &p, lj = &q; + else lj = &p, li = &q; + os = MAX(li->qs, lj->qs), oe = MIN(li->qe, lj->qe); + ovlp = ((oe > os)? (oe - os):0); + if(!ovlp) continue;//no overlap + + if(lj->qs <= li->qs+G_CHAIN_INDEL) { + if(govlp_check(uref, uopt, bw, diff_ec_ul, li, lj)) { + idx->a[k]++; kv_push(uint64_t, *idx, i); + } + } else if((lj->qe+G_CHAIN_INDEL>=li->qe) && (lj->qs+G_CHAIN_INDEL>=li->qs)) { + if(govlp_check(uref, uopt, bw, diff_ec_ul, lj, li)) { + idx->a[k]++; kv_push(uint64_t, *idx, i); + } + } + } + } + + + // for (k = 0; k < on; k++) { + // int64_t s, e; + // s = idx->a[k]>>32; e = s + (uint32_t)idx->a[k]; + // for (i = s; i < e; i++) { + // fprintf(stderr, "k::%ld[M::%s::utg%.6dl] utg%.6dl\n", k, __func__, + // (int32_t)ol->list[k].y_id+1, (int32_t)ol->list[idx->a[i]].y_id+1); + // } + // } +} + +void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t *uopt, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, asg64_v* buf1) +{ + int64_t on = ol->length, k, i, zwn, q[2], t[2], w[2]; + uint64_t m; overlap_region *z; ul_ov_t *cp; + kv_resize(uint64_t, *idx, (ol->length<<1)); + kv_resize(ul_ov_t, *c_idx, ol->length); + for (k = idx->n = c_idx->n = 0; k < on; k++) { + z = &(ol->list[k]); zwn = z->w_list.n; + z->align_length = z->overlapLen = z->x_pos_e+1-z->x_pos_s; + if(!zwn) continue; + q[0] = q[1] = t[0] = t[1] = w[0] = w[1] = INT32_MIN; + for (i = 0; i < zwn; i++) { + if((z->w_list.a[i].x_start==(q[1]+1)) && ((z->w_list.a[i].y_start==(t[1]+1)))) { + q[1] = z->w_list.a[i].x_end; + t[1] = z->w_list.a[i].y_end; + w[1] = i; + } else { + if(q[0] != INT32_MIN) { + m = ((uint64_t)q[0])<<1; m <<= 32; + m += c_idx->n; kv_push(uint64_t, *idx, m); + m = (((uint64_t)q[1])<<1)+1; m <<= 32; + m += c_idx->n; kv_push(uint64_t, *idx, m); + + kv_pushp(ul_ov_t, *c_idx, &cp); + ovlp_id(*cp) = k; ///ovlp id + ovlp_min_wid(*cp) = w[0]; ///beg id of windows + ovlp_max_wid(*cp) = w[1]; ///end id of windows + ovlp_cur_wid(*cp) = w[0]; ///cur id of windows + ovlp_cur_xoff(*cp) = z->w_list.a[w[0]].x_start; ///cur xpos + ovlp_cur_coff(*cp) = 0; ///cur cigar off in cur window + } + + q[0] = z->w_list.a[i].x_start; q[1] = z->w_list.a[i].x_end; + t[0] = z->w_list.a[i].y_start; t[1] = z->w_list.a[i].y_end; + w[0] = i; w[1] = i; + } + } + if(q[0] != INT32_MIN) { + m = ((uint64_t)q[0])<<1; m <<= 32; + m += c_idx->n; kv_push(uint64_t, *idx, m); + m = (((uint64_t)q[1])<<1)+1; m <<= 32; + m += c_idx->n; kv_push(uint64_t, *idx, m); + + kv_pushp(ul_ov_t, *c_idx, &cp); + ovlp_id(*cp) = k; ///ovlp id + ovlp_min_wid(*cp) = w[0]; ///beg id of windows + ovlp_max_wid(*cp) = w[1]; ///end id of windows + ovlp_cur_wid(*cp) = w[0]; ///cur id of windows + ovlp_cur_xoff(*cp) = z->w_list.a[w[0]].x_start; ///cur xpos + ovlp_cur_coff(*cp) = 0; ///cur cigar off in cur window + } + } + radix_sort_bc64(idx->a, idx->a+idx->n); + kv_resize(uint64_t, *buf, idx->n); + gen_gov_idx(ol, uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, buf1); + // for (m = 0; m < c_idx->n; m++) { + // fprintf(stderr, "+++[M::%s::utg%.6dl] q[%d, %d), t[%d, %d)\n", __func__, + // (int32_t)ol->list[ovlp_id(c_idx->a[m])].y_id+1, + // ol->list[ovlp_id(c_idx->a[m])].w_list.a[ovlp_min_wid(c_idx->a[m])].x_start, + // ol->list[ovlp_id(c_idx->a[m])].w_list.a[ovlp_max_wid(c_idx->a[m])].x_end+1, + // ol->list[ovlp_id(c_idx->a[m])].w_list.a[ovlp_min_wid(c_idx->a[m])].y_start, + // ol->list[ovlp_id(c_idx->a[m])].w_list.a[ovlp_max_wid(c_idx->a[m])].y_end+1); + // } + + + int64_t srt_n = idx->n, dp, old_dp, beg, end; + for (i = k = 0, dp = old_dp = 0, beg = 0, end = -1; i < srt_n; ++i) {///[beg, end) but coordinates in idx is [, ] + ///if idx->a.a[] is qe + old_dp = dp; + 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])); + } + // fprintf(stderr, "\n[M::%s::] beg::%ld, end::%ld, old_dp::%ld\n", __func__, beg, end, old_dp); + if((end > beg) && (old_dp >= 2)) { + idx->n = srt_n + gen_region_phase(ol->list, idx->a+srt_n, idx->n-srt_n, beg, end, old_dp, c_idx->a, buf->a, buf1); + } + beg = end; + } + ///hap->length + + +} void ul_gap_filling_adv(overlap_region_alloc* ol, Candidates_list *cl, 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, overlap_region *aux_o, asg64_v* buf, asg64_v* iidx, double e_rate, int64_t ql, uint64_t rid, int64_t khit, int64_t base_chekc_k_hit, int64_t max_lgap) { - int64_t k, l, ch_n, a_n = aln->n; overlap_region *z; //k_mer_hit *ch_a; + int64_t k, l, ch_n, a_n = aln->n; uint64_t pqn, pk; overlap_region *z; //k_mer_hit *ch_a; // count_k_hits(rref, uref, qstr, tu, ol, cl, buf, khit, base_chekc_k_hit); count_k_hits_adv(rref, uref, qstr, tu, ol, cl, buf, &(cl->chainDP), e_rate, khit, base_chekc_k_hit); - for (k = 1, l = 0; k <= a_n; k++) { + for (k = 1, l = 0, pqn = 0; k <= a_n; k++) { if(k == a_n || aln->a[l].qn != aln->a[k].qn) { - z = &(ol->list[aln->a[l].qn]); assert(z->align_length == l); + z = &(ol->list[aln->a[l].qn]); assert(z->align_length == l); + + for (pk = pqn; pk < aln->a[l].qn; pk++) ol->list[pk].w_list.n = 0; + pqn = aln->a[l].qn+1; + ch_n = gen_cns_chain(z, cl, iidx, max_lgap, e_rate, 0); 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, 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); l = k; } } + for (pk = pqn; pk < ol->length; pk++) ol->list[pk].w_list.n = 0; } inline uint32_t ovlp_win_check(overlap_region *z, uint32_t id0, uint32_t id1, int64_t max_lgap, double small_bw_rate, int64_t min_small_bw) @@ -15054,13 +15588,13 @@ kv_ul_ov_t *aln, uint64_t rid, int64_t max_lgap, double sgap_rate) -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, - haplotype_evdience_alloc* hap, kvec_t_u64_warp* v_idx, overlap_region *aux_o, - double e_rate, int64_t wl, kv_ul_ov_t *aln, int64_t sid, uint64_t khit, void *km) +void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, const ug_opt_t *uopt, + char *qstr, uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, bit_extz_t *exz, haplotype_evdience_alloc* hap, + kvec_t_u64_warp* v_idx, overlap_region *aux_o, double e_rate, int64_t wl, kv_ul_ov_t *aln, int64_t sid, uint64_t khit, + st_mt_t *stb, void *km) { uint64_t i, bs, k, ovl/**, on**/; Window_Pool w; double err; - /**int64_t sc;**/ overlap_region t; overlap_region *z; asg64_v iidx, buf; + /**int64_t sc;**/ overlap_region t; overlap_region *z; asg64_v iidx, buf, buf1; ol->mapped_overlaps_length = 0; if(ol->length <= 0) return; @@ -15110,10 +15644,14 @@ void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur // fprintf(stderr, "-[M::%s] on::%lu\n", __func__, ol->length); if(ol->length <= 1) return; ///coordinates for all intervals with cov > 1 - copy_asg_arr(iidx, hap->snp_srt); copy_asg_arr(buf, v_idx->a); + copy_asg_arr(iidx, hap->snp_srt); copy_asg_arr(buf, v_idx->a); copy_asg_arr(buf1, (*stb)); // fprintf(stderr, "\n[M::%s] iidx_n::%ld\n", __func__, (int64_t)iidx.n); ul_gap_filling_adv(ol, cl, aln, wl, uref, NULL, NULL, qu->seq, tu, exz, aux_o, &buf, &iidx, err, ql, sid, khit, 1, MAX_LGAP(ql)); - copy_asg_arr(hap->snp_srt, iidx); copy_asg_arr(v_idx->a, buf); + copy_asg_arr(hap->snp_srt, iidx); copy_asg_arr(v_idx->a, buf); copy_asg_arr((*stb), buf1); + + copy_asg_arr(iidx, hap->snp_srt); copy_asg_arr(buf, v_idx->a); copy_asg_arr(buf1, (*stb)); + region_phase(ol, uref, uopt, aln, &iidx, &buf, &buf1); + copy_asg_arr(hap->snp_srt, iidx); copy_asg_arr(v_idx->a, buf); copy_asg_arr((*stb), buf1); // 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; @@ -15132,19 +15670,19 @@ void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur // ul_phase(ol->list, on, v_idx->a.a, v_idx->a.a+on); - for (i = 0; i < ol->length; i++) { - z = &(ol->list[i]); ovl = z->x_pos_e+1-z->x_pos_s; z->is_match = 1; - for (k = 0; k < z->w_list.n; k++) { - if(z->w_list.a[k].clen) continue; - gen_backtrace_adv_exz(&(z->w_list.a[k]), z, NULL, NULL, uref, qu->seq, tu->seq, exz, z->y_pos_strand, z->y_id); - } - // verify_aln(sid, z, qu, tu, NULL, NULL, uref); + // for (i = 0; i < ol->length; i++) { + // z = &(ol->list[i]); ovl = z->x_pos_e+1-z->x_pos_s; z->is_match = 1; + // for (k = 0; k < z->w_list.n; k++) { + // if(z->w_list.a[k].clen) continue; + // gen_backtrace_adv_exz(&(z->w_list.a[k]), z, NULL, NULL, uref, qu->seq, tu->seq, exz, z->y_pos_strand, z->y_id); + // } + // // verify_aln(sid, z, qu, tu, NULL, NULL, uref); - ol->mapped_overlaps_length += ovl; - append_unmatched_wins(z, w.window_length); - calculate_ul_boundary_cigars(z, uref, dumy, qu, err, w.window_length); - } - partition_ul_overlaps_advance(ol, uref, qu, tu, dumy, hap, 1, err, w.window_length, km); + // ol->mapped_overlaps_length += ovl; + // append_unmatched_wins(z, w.window_length); + // calculate_ul_boundary_cigars(z, uref, dumy, qu, err, w.window_length); + // } + // partition_ul_overlaps_advance(ol, uref, qu, tu, dumy, hap, 1, err, w.window_length, km); } } \ No newline at end of file diff --git a/Correct.h b/Correct.h index cf076d8..37b9e4c 100644 --- a/Correct.h +++ b/Correct.h @@ -1134,10 +1134,10 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref kvec_t_u64_warp* v_idx, window_list_alloc* win_ciagr_buf, int force_repeat, int is_consensus, int* fully_cov, int* abnormal, double max_ov_diff_ec, long long winLen, void *km); -void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, char *qstr, +void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, const ug_opt_t *uopt, char *qstr, uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, bit_extz_t *exz, haplotype_evdience_alloc* hap, kvec_t_u64_warp* v_idx, overlap_region *aux_o, - double e_rate, int64_t wl, kv_ul_ov_t *aln, int64_t sid, uint64_t hpc_k, void *km); + double e_rate, int64_t wl, kv_ul_ov_t *aln, int64_t sid, uint64_t hpc_k, st_mt_t *stb, void *km); void ul_lalign_old_ed(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, @@ -1330,6 +1330,12 @@ void update_sketch_trace(overlap_region_alloc* ol, const ul_idx_t *uref, const u 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, int64_t max_lgap, double sgap_rate); +int64_t infer_rovlp(ul_ov_t *li, ul_ov_t *lj, uc_block_t *bi, uc_block_t *bj, All_reads *ridx, ma_ug_t *ug); +void convert_ul_ov_t(ul_ov_t *des, overlap_region *src, const ul_idx_t *uref); +uint64_t check_connect_ug(const ul_idx_t *uref, uint32_t v, uint32_t w, int64_t bw, double diff_ec_ul, int64_t dq); +uint64_t check_connect_rg(const ul_idx_t *uref, const ug_opt_t *uopt, uint32_t uv, uint32_t uw, int64_t bw, double diff_ec_ul, int64_t dq); +uint32_t govlp_check(const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, ul_ov_t *li, ul_ov_t *lj); + #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)) diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index 706d33d..ad370fc 100644 --- a/Levenshtein_distance.h +++ b/Levenshtein_distance.h @@ -537,6 +537,15 @@ inline uint32_t pop_trace(asg16_v *res, uint32_t i, uint16_t *c, uint32_t *len) return i; } +inline int32_t pop_trace_back(asg16_v *res, int32_t i, uint16_t *c, uint32_t *len) +{ + (*c) = (res->a[i]>>14); (*len) = (res->a[i]&(0x3fff)); + for (i--; (i >= 0) && ((*c) == (res->a[i]>>14)); i--) { + (*len) += (res->a[i]&(0x3fff)); + } + return i; +} + ///511 -> 16 64-bits // #define MAX_E 511 // #define MAX_L 2500 diff --git a/inter.cpp b/inter.cpp index d3dabfa..a00dc38 100644 --- a/inter.cpp +++ b/inter.cpp @@ -18,6 +18,9 @@ #include "Assembly.h" KSEQ_INIT(gzFile, gzread) +#define oreg_xe_lt(a, b) (((uint64_t)(a).x_pos_e<<32|(a).x_pos_s) < ((uint64_t)(b).x_pos_e<<32|(b).x_pos_s)) +KSORT_INIT(or_xe, overlap_region, oreg_xe_lt) + void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t high_occ, void *km); void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, @@ -46,6 +49,9 @@ void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t #define generic_key(x) (x) KRADIX_SORT_INIT(gfa64, uint64_t, generic_key, 8) +#define generic_key(x) (x) +KRADIX_SORT_INIT(gfa64i, int64_t, generic_key, 8) + #define ul_ov_srt_qe_key(p) ((p).qe) KRADIX_SORT_INIT(ul_ov_srt_qe, ul_ov_t, ul_ov_srt_qe_key, member_size(ul_ov_t, qe)) @@ -2964,83 +2970,6 @@ ma_hit_t* query_ovlp_src(const ug_opt_t *uopt, uint32_t v, uint32_t w, int64_t o return NULL; } -int64_t infer_rovlp(ul_ov_t *li, ul_ov_t *lj, uc_block_t *bi, uc_block_t *bj, All_reads *ridx, ma_ug_t *ug) -{ - int64_t in, is, ie, irev, iqs, iqe, jn, js, je, jrev, jqs, jqe, ir, jr, ts, te, max_s, min_e, s_shift, e_shift; - - if(li) { - in = ug?ug->u.a[li->tn].len:Get_READ_LENGTH(R_INF, li->tn); - is = li->ts; ie = li->te; irev = li->rev; iqs = li->qs; iqe = li->qe; - } else if(bi) { - in = ug?ug->u.a[bi->hid].len:Get_READ_LENGTH(R_INF, bi->hid); - is = bi->ts; ie = bi->te; irev = bi->rev; iqs = bi->qs; iqe = bi->qe; - } else { - return 0; - } - - if(lj) { - jn = ug?ug->u.a[lj->tn].len:Get_READ_LENGTH(R_INF, lj->tn); - js = lj->ts; je = lj->te; jrev = lj->rev; jqs = lj->qs; jqe = lj->qe; - } else if(bj) { - jn = ug?ug->u.a[bj->hid].len:Get_READ_LENGTH(R_INF, bj->hid); - js = bj->ts; je = bj->te; jrev = bj->rev; jqs = bj->qs; jqe = bj->qe; - } else { - return 0; - } - - max_s = MAX(iqs, jqs); min_e = MIN(iqe, jqe); - if(min_e <= max_s) return 0; - s_shift = get_offset_adjust(max_s - iqs, iqe-iqs, ie-is); - e_shift = get_offset_adjust(iqe - min_e, iqe-iqs, ie-is); - if(irev) { - ts = s_shift; s_shift = e_shift; e_shift = ts; - } - is += s_shift; ie-= e_shift; - - // if(li && lj && li->tn == 324 && lj->tn == 319 && li->qs == 63841) { - // fprintf(stderr, "+++in:%ld, is:%ld, ie:%ld, irev:%ld, jn:%ld, js:%ld, je:%ld, jrev:%ld\n", in, is, ie, irev, jn, js, je, jrev); - // } - - s_shift = get_offset_adjust(max_s - jqs, jqe-jqs, je-js); - e_shift = get_offset_adjust(jqe - min_e, jqe-jqs, je-js); - if(jrev) { - ts = s_shift; s_shift = e_shift; e_shift = ts; - } - js += s_shift; je-= e_shift; - - if(irev) { - ts = in - ie; te = in - is; - is = ts; ie = te; - } - - if(jrev) { - ts = jn - je; te = jn - js; - js = ts; je = te; - } - - // if(li && lj && li->tn == 324 && lj->tn == 319 && li->qs == 63841) { - // fprintf(stderr, "---in:%ld, is:%ld, ie:%ld, irev:%ld, jn:%ld, js:%ld, je:%ld, jrev:%ld\n", in, is, ie, irev, jn, js, je, jrev); - // } - - if(is <= js) { - js -= is; is = 0; - } else { - is -= js; js = 0; - } - - ir = in - ie; jr = jn - je; - - if(ir <= jr){ - ie = in; je += ir; - } - else { - je = jn; ie += jr; - } - - ir = ie - is; jr = je - js; - return MAX(ir, jr); -} - void debug_infer_read_ovlp(const ug_opt_t *uopt, double diff_ec_ul, ul_ov_t *li, ul_ov_t *lj, ma_utg_t *u, uint32_t i_idx, uint32_t j_idx, All_reads *ridx, ma_ug_t *ug) { @@ -5102,6 +5031,472 @@ int64_t debug_i, int64_t tid, void *km) return 1; } +#define aln_sc(a, w) (((int64_t)((a).sec))-((int64_t)(((a).qe-(a).qs-(a).sec)*(w)))) + +void gen_gl_aln(overlap_region_alloc* olist, const ul_idx_t *uref, kv_ul_ov_t *res) +{ + uint64_t k; ul_ov_t *p = NULL; + res->n = 0; kv_resize(ul_ov_t, *res, olist->length); + for (k = 0; k < olist->length; k++) { + // fprintf(stderr, "+++[M::%s::utg%.6dl] q[%d, %d), t[%d, %d), tot::%u, cis::%u\n", __func__, + // (int32_t)olist->list[k].y_id+1, olist->list[k].x_pos_s, olist->list[k].x_pos_e+1, + // olist->list[k].y_pos_s, olist->list[k].y_pos_e+1, + // olist->list[k].overlapLen, olist->list[k].align_length); + p = &(res->a[res->n++]); + p->qn = k; p->qs = olist->list[k].x_pos_s; p->qe = olist->list[k].x_pos_e+1; + p->tn = olist->list[k].y_id; p->sec = olist->list[k].align_length; + p->rev = olist->list[k].y_pos_strand; p->el = (olist->list[k].is_match==1?1:0); + if(p->rev) { + p->ts = uref->ug->u.a[p->tn].len - (olist->list[k].y_pos_e+1); + p->te = uref->ug->u.a[p->tn].len - olist->list[k].y_pos_s; + } else { + p->ts = olist->list[k].y_pos_s; + p->te = olist->list[k].y_pos_e+1; + } + } +} + +void gen_gg_aln(overlap_region_alloc* olist, const ul_idx_t *uref, int64_t trans_sc, vec_mg_lchain_t *res) +{ + uint64_t k; mg_lchain_t *p = NULL; + res->n = 0; kv_resize(mg_lchain_t, *res, olist->length); + for (k = 0; k < olist->length; k++) { + p = &(res->a[res->n++]); memset(p, 0, sizeof((*p))); + p->v = ((olist->list[k].y_id<<1)|(olist->list[k].y_pos_strand)); p->off = k; + p->score = (((int64_t)(olist->list[k].align_length)) + -((int64_t)((olist->list[k].overlapLen-olist->list[k].align_length)*(trans_sc)))); + p->qs = olist->list[k].x_pos_s; p->qe = olist->list[k].x_pos_e+1; + if((p->v&1)) { + p->rs = uref->ug->u.a[p->v>>1].len - (olist->list[k].y_pos_e+1); + p->re = uref->ug->u.a[p->v>>1].len - olist->list[k].y_pos_s; + } else { + p->rs = olist->list[k].y_pos_s; + p->re = olist->list[k].y_pos_e+1; + } + } +} + + + +int64_t gl_chain_lin(kv_ul_ov_t *res, overlap_region *ol, ul_ov_t *ex, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, +double diff_ec_ul, int64_t qlen, int64_t max_skip, int64_t max_iter, int64_t max_dis, Chain_Data* dp, int64_t trans_sc, +uint64_t mode, All_reads *ridx, ma_ug_t *ug, int64_t need_srt) +{ + if(res->n == 0) return 0; + uint32_t li_v, lj_v, rev_n; int32_t *f, *c_n, *c_sc; int64_t *p, *t, res_n = res->n, st, max_ii, max; + int64_t mm_ovlp, x, i, j, k, sc, csc, mm_sc, mm_idx, qo, share, n_skip, end_j, plus; ul_ov_t *li, *lj, rev_t; + resize_Chain_Data(dp, res_n, NULL); + t = dp->tmp; f = dp->score; p = dp->pre; c_n = dp->occ; c_sc = dp->self_length; + if(need_srt) { + radix_sort_ul_ov_srt_qe(res->a, res->a + res_n); + for (i = 1, j = 0; i <= res_n; i++) { + if (i == res_n || res->a[i].qe != res->a[j].qe) { + if(i - j > 1) { + radix_sort_ul_ov_srt_qs(res->a+j, res->a+i); + } + j = i; + } + } + } + + memset(t, 0, (res_n*sizeof((*t)))); + for (i = st = plus = 0, max_ii = -1; i < res_n; ++i) { + li = &(res->a[i]); li_v = (li->tn<<1)|li->rev; + mm_ovlp = mode?max_ovlp_src(uopt, li_v^1):max_ovlp(uref->ug->g, li_v^1); + x = (li->qs + mm_ovlp)*diff_ec_ul; + if(x < bw) x = bw; + x += li->qs + mm_ovlp; + if (x > qlen+1) x = qlen+1; + x = find_ul_ov_max(i, res->a, x+G_CHAIN_INDEL); + csc = aln_sc((*li), trans_sc); + mm_sc = csc; mm_idx = -1; + + n_skip = 0; end_j = -1; + if ((x-st) > max_iter) st = x-max_iter; + for (j = x; j >= st; --j) { // collect potential destination vertices + lj = &(res->a[j]); lj_v = (lj->tn<<1)|lj->rev; + if(lj->qe+G_CHAIN_INDEL <= li->qs) break;//even this pair has a overlap, its length will be very small; just ignore + if(lj->qs >= li->qs) continue; + qo = infer_rovlp(li, lj, NULL, NULL, ridx, ug); ///overlap length in query (UL read) + if(li_v != lj_v && get_ecov_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, mode, &share)) { + sc = csc + f[j]; + if(sc > mm_sc) { + mm_sc = sc, mm_idx = 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; + } + } + + end_j = j; + if (max_ii < 0 || (res->a[i].qe>(res->a[max_ii].qe+max_dis))) {//too long + max = INT32_MIN; max_ii = -1; + for (j = i - 1; (j >= st) && (res->a[i].qe<=(max_dis+res->a[j].qe)); --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] + lj = &(res->a[max_ii]); lj_v = (lj->tn<<1)|lj->rev; + if(lj->qe+G_CHAIN_INDEL <= li->qs) break;//even this pair has a overlap, its length will be very small; just ignore + if(lj->qs >= li->qs) continue; + qo = infer_rovlp(li, lj, NULL, NULL, ridx, ug); ///overlap length in query (UL read) + if(li_v != lj_v && get_ecov_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, mode, &share)) { + sc = csc + f[j]; + if(sc > mm_sc) { + mm_sc = sc; mm_idx = max_ii; + } + } + } + if(mm_sc < 0) { + mm_sc = csc; mm_idx = -1; + } + f[i] = mm_sc; p[i] = mm_idx; + + if ((max_ii < 0) || ((res->a[i].qe<=max_dis+res->a[max_ii].qe) && (f[max_ii]= 0; --k) { + n_v0 = n_v; + for (i = ((uint32_t)t[k])>>1; i >= 0 && (t[i]&1) == 0; ) { + ex[n_v++] = res->a[i]; t[i] |= 1; i = p[i]; + } + if(n_v0 == n_v) continue; + sc = (i<0?(t[k]>>32):((t[k]>>32)-f[i])); + c_n[n_u] = n_v-n_v0; c_sc[n_u] = sc; n_u++; + } + + for (k = 0, n_v = n_v0 = 0; k < n_u; k++) { + n_v0 = n_v; n_v += c_n[k]; + res->a[k].qn = c_sc[k];//score + res->a[k].ts = n_v0; res->a[k].te = n_v;///idx + + rev_n = c_n[k]>>1; + ///we need to consider contained reads; so determining qs is not such easy + res->a[k].qs = (uint32_t)-1; res->a[k].qe = ex[n_v0].qe; + for (i = 0; i < rev_n; i++) { + rev_t = ex[n_v0+i]; ex[n_v0+i] = ex[n_v-i-1]; ex[n_v-i-1] = rev_t; + + if(res->a[k].qs > ex[n_v0+i].qs) res->a[k].qs = ex[n_v0+i].qs; + if(res->a[k].qs > ex[n_v-i-1].qs) res->a[k].qs = ex[n_v-i-1].qs; + ex[n_v0+i].sec = ex[n_v-i-1].sec = SEC_MODE; + } + if(c_n[k]&1) { + if(res->a[k].qs > ex[n_v0+i].qs) res->a[k].qs = ex[n_v0+i].qs; + ex[n_v0+i].sec = SEC_MODE; + } + } + res->n = n_u; + radix_sort_ul_ov_srt_qn(res->a, res->a + res->n);//sort by score + // fprintf(stderr, "---[M::%s] n_u:%ld, n_v:%ld\n", __func__, n_u, n_v); + return n_v; +} + +void set_ul_ov_t_by_mg_lchain_t(ul_ov_t *u, mg_lchain_t *l) +{ + u->tn = l->v>>1; u->rev = (l->v&1); + u->ts = l->rs; u->te = l->re; + u->qs = l->qs; u->qe = l->qe; +} + +int64_t find_mg_lchain_max(int64_t n, const mg_lchain_t *a, int32_t x) +{ + int64_t s = 0, e = n; + if (n == 0) return -1; + if (a[n-1].qe < x) return n - 1; + if (a[0].qe >= x) return -1; + while (e > s) { // TODO: finish this block + int64_t m = s + (e - s) / 2; + if (a[m].qe >= x) e = m; + else s = m + 1; + } + assert(s == e); + return s; +} + + +int64_t hc_target_len(asg_t *g, mg_lchain_t *s, mg_lchain_t *e) +{ + // int64_t ql = s->qe - e->qe, tp, tm; + // if((s->v^1)&1) tp = g->seq[s->v>>1].len - s->re; + // else tp = s->rs; + + // if((e->v^1)&1) tm = g->seq[e->v>>1].len - e->re; + // else tm = e->rs; + int64_t ql = (int64_t)s->qs - (int64_t)e->qe, tp, tm; + int64_t sts, ete; + sts = (s->v&1)?g->seq[s->v>>1].len-s->re:s->rs; tp = g->seq[s->v>>1].len - sts; + ete = (e->v&1)?g->seq[e->v>>1].len-e->rs:e->re; tm = g->seq[e->v>>1].len - ete; + // fprintf(stderr, "[M::%s::] ql:%ld, tp:%ld, tm:%ld, sts:%ld, ete:%ld\n", __func__, ql, tp, tm, sts, ete); + return ql + tp - tm; +} + +inline int32_t cal_gchain_sc(const mg_path_dst_t *dj, const mg_lchain_t *li, const mg_lchain_t *lc, int64_t *f, int64_t b_w, float diff_thre, float chn_pen_gap) +{ + // const mg_lchain_t *lj; + int32_t gap, sc; + float lin_pen, log_pen; + if (dj->n_path == 0) return INT32_MIN; + gap = dj->dist - dj->target_dist; + // lj = &lc[dj->meta]; + if (gap < 0) gap = -gap; + if ((gap > ((dj->target_dist)*diff_thre)) && (gap > b_w)) return INT32_MIN; + // if (lj->qe <= li->qs) sc = li->score; + // else sc = (int32_t)((double)(li->qe - lj->qe) / (li->qe - li->qs) * li->score + .499); // dealing with overlap on query + sc = li->score; + //sc += dj->mlen; // TODO: is this line the right thing to do? + // if (dj->is_0) sc += ref_bonus; + lin_pen = chn_pen_gap * (float)gap; + log_pen = gap >= 2? mg_log2(gap) : 0.0f; + sc -= (int32_t)(lin_pen + log_pen); + sc += f[dj->meta]; + return sc; +} + +int64_t gl_chain_graph(void *km, const ul_idx_t *uref, const ma_ug_t *ug, vec_mg_lchain_t *lc, +vec_mg_lchain_t *sw, vec_mg_path_dst_t *dst, vec_sp_node_t *out, vec_mg_pathv_t *path, +int64_t qlen, const ug_opt_t *uopt, int64_t bw, double diff_thre, uint64_t *srt, st_mt_t *bf, +Chain_Data* dp, int64_t max_skip, int64_t need_srt) +{ + bf->n = 0; + if(lc->n == 0) return 0; + int64_t i, j, lc_n = lc->n, n_ext, mm_ovlp, target_dist, max_target_dist, x, m_idx, m_sc, qo, sc; + int64_t max_f, max_j = -1, max_d = -1, max_inner = 0, share; uint32_t max_hash = 0; int64_t k, k0, n_u, n_v, ni; + mg_lchain_t *r, *li, *lj; mg_path_dst_t *q; asg_t *g = ug->g; uint64_t isolated, *u; ul_ov_t ui, uj; + if(!need_srt) { + for (i = n_ext = 0; i < lc_n; i++) { + r = &lc->a[i]; r->dist_pre = -1; isolated = 0;///dist_pre -> parent in graph chain + if((r->re < g->seq[r->v>>1].len) && (r->rs > 0)) isolated = 1;///UL contained in one vertice + if (!isolated) { + srt[n_ext] = r->qe; srt[n_ext] <<= 32; + srt[n_ext] |= (uint64_t)i; srt[n_ext] |= (isolated<<63); + ++n_ext; + } + } + j = n_ext; + if(j < lc_n) { + for (i = 0; i < lc_n; i++) { + r = &lc->a[i]; r->dist_pre = -1; isolated = 0;///dist_pre -> parent in graph chain + if((r->re < g->seq[r->v>>1].len) && (r->rs > 0)) isolated = 1;///UL contained in one vertice + if (isolated) { + srt[j] = r->qe; srt[j] <<= 32; srt[j] |= (uint64_t)i; srt[j] |= (isolated<<63); ++j; + } + } + } + assert(j == lc_n); + } else { + for (i = n_ext = 0; i < lc_n; i++) { + r = &lc->a[i]; r->dist_pre = -1; isolated = 0;///dist_pre -> parent in graph chain + if((r->re < g->seq[r->v>>1].len) && (r->rs > 0)) isolated = 1;///UL contained in one vertice + if (!isolated) ++n_ext; + srt[i] = r->qe; srt[i] <<= 32; srt[i] |= (uint64_t)i; srt[i] |= (isolated<<63); + } + radix_sort_gfa64(srt, srt+lc_n); + for (i = 1, j = 0; i <= lc_n; i++) { + if (i == lc_n || (srt[i]>>32) != (srt[j]>>32)) { + if(i - j > 1) { + for (x = j; x < i; x++) { + srt[x] <<= 32; srt[x] >>= 32; srt[x] |= ((uint64_t)lc->a[(uint32_t)srt[x]].qs)<<32; + } + radix_sort_gfa64(srt+j, srt+i); + } + j = i; + } + } + } + if((n_ext != lc_n) || (need_srt)) { + kv_resize(mg_lchain_t, *sw, (uint64_t)lc_n); sw->n = lc_n; + for (i = 0; i < lc_n; i++) sw->a[i] = lc->a[(uint32_t)srt[i]]; + memcpy(lc->a, sw->a, lc_n *sizeof((*(lc->a)))); + } + + + resize_Chain_Data(dp, lc_n, NULL); + int32_t *p; int64_t *f, *t, n_skip, dst_n, is_f, plus, n_v0; mg_path_dst_t *dj; + t = dp->tmp; p = dp->score; f = dp->pre; + memset(t, 0, (n_ext*sizeof((*t)))); + + for (i = plus = 0; i < n_ext; ++i) { // core loop + li = &lc->a[i]; set_ul_ov_t_by_mg_lchain_t(&ui, li); + mm_ovlp = max_ovlp(g, li->v^1); + x = (li->qs + mm_ovlp)*diff_thre; if(x < bw) x = bw; + x += li->qs + mm_ovlp; + if (x > qlen+1) x = qlen+1; + x = find_mg_lchain_max(i, lc->a, x+G_CHAIN_INDEL); + + n_skip = 0; is_f = 0; + // collect potential destination vertices + for (dst->n = 0, max_target_dist = -1, j = x; j >= 0; --j) { + lj = &lc->a[j]; ///extend_end_coord(lj, qlen, g->seq[lj->v>>1].len, &jqs, &jqe, &jrs, &jre); + //lj contained in li; actually in circle, this might happen; need to deal with it later + if(lj->qs >= li->qs/**+G_CHAIN_INDEL**/) continue; + ///if there is a circle, the two linear chains might be at the same vertice + target_dist = hc_target_len(g, li, lj); + if(target_dist < 0) continue; + kv_pushp(mg_path_dst_t, *dst, &q); + memset(q, 0, sizeof(*q)); + q->inner = 0;//we set q->inner = 0 to allow circles + q->v = lj->v^1;///must be v^1 instead of v + q->meta = j; + ///lj->qs************lj->qe + /// li->qs************li->qe + q->qlen = li->qs - lj->qe;///might be negative; this is the region that need to be checked in base-level + q->target_dist = target_dist;///cannot understand the target_dist + q->target_hash = 0; + q->check_hash = 0; + if(max_target_dist < target_dist) max_target_dist = target_dist; + if (t[j] == i) { + if (++n_skip > max_skip) + break; + } + if (p[j] >= 0) t[p[j]] = i; + if((!is_f) && (lj->qe+G_CHAIN_INDEL > li->qs)) { + set_ul_ov_t_by_mg_lchain_t(&uj, lj); + qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug); + if(li->v!=lj->v && get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, N_GCHAIN_RATE, qo, 0, &share)) { + is_f = 1; if(n_skip > 0) n_skip--; + if(n_skip < (max_skip>>1)) n_skip= (max_skip>>1); + } + } + } + + // confirm reach-ability + max_f = li->score, max_j = -1, max_d = -1, max_inner = 0; max_hash = 0; + if(dst->n) { + max_target_dist *= (1+diff_thre); if(max_target_dist < bw) max_target_dist = bw; + hc_shortest_k(km, g, li->v^1, dst->n, dst->a, max_target_dist, MG_MAX_SHORT_K, bf, srt, out, NULL, 1, 0, 0); + // remove unreachable destinations + //TODO: check sequence identity + dst_n = dst->n; + for (j = 0; j < dst_n; ++j) { + dj = &dst->a[j]; + if (dj->n_path == 0) continue; // unreachable + sc = cal_gchain_sc(dj, li, lc->a, f, bw, diff_thre, W_CHN_PEN_GAP); + + if (sc == INT32_MIN) continue; // out of band + // if (sc < 0) continue;// negative score + if (sc > max_f) { + max_f = sc, max_j = dj->meta, max_d = dj->dist, max_hash = dj->hash, max_inner = dj->inner; + } + } + } + if(max_f < 0) { + max_f = li->score; max_j = -1; + } + + f[i] = max_f; p[i] = max_j; + ///same time for gchain + li->dist_pre = max_d; + li->hash_pre = max_hash; + li->inner_pre = max_inner; + if(max_f < plus) plus = max_f;//minmun negative + } + + for (; i < lc_n; i++) { + li = &lc->a[i]; + max_f = li->score, max_j = -1, max_d = -1, max_inner = 0; max_hash = 0; + f[i] = max_f; p[i] = max_j; + ///same time for gchain + li->dist_pre = max_d; + li->hash_pre = max_hash; + li->inner_pre = max_inner; + if(max_f < plus) plus = max_f;//minmun negative + } + + for (i = 0; i < lc_n; ++i) { + f[i]-=plus; t[i] = f[i]<<32; t[i] += (i<<1); + } + + sw->n = 0; kv_resize(mg_lchain_t, *sw, (uint64_t)lc_n); + kv_resize(uint64_t, *bf, (uint64_t)lc_n); u = bf->a; + n_u = n_v = 0; radix_sort_gfa64i(t, t + lc_n); + for (k = lc_n-1, n_v = n_u = 0; k >= 0; --k) { + n_v0 = n_v; + for (i = ((uint32_t)t[k])>>1; i >= 0 && (t[i]&1) == 0; ) { + sw->a[n_v++] = lc->a[i]; t[i] |= 1; i = p[i]; + } + if(n_v0 == n_v) continue; + sc = (i<0?(t[k]>>32):((t[k]>>32)-f[i])); + u[n_u++] = (((uint64_t)sc)<<32) | ((uint64_t)(n_v-n_v0)); + } + + m_idx = m_sc = -1; + for (i = 0, k = 0; i < n_u; ++i) { + k0 = k, ni = (int32_t)u[i]; + for (j = 0; j < ni; ++j) { + lc->a[k++] = sw->a[k0 + (ni - j - 1)]; + } + if(m_idx < 0 || m_sc < ((int64_t)(u[i]>>32))) { + m_idx = i; m_sc = ((int64_t)(u[i]>>32)); + } + } + assert(k == n_v); bf->n = n_u; + return m_idx; +} + +int64_t gl_chain(mg_tbuf_t *b, ul_vec_t *rch, overlap_region_alloc* olist, +Chain_Data* dp, haplotype_evdience_alloc *hap, st_mt_t *sps, glchain_t *ll, +gdpchain_t *gdp, const ul_idx_t *uref, double diff_ec_ul, int64_t winLen, +int64_t qlen, const ug_opt_t *uopt, int64_t debug_i, int64_t tid, void *km) +{ + ll->tk.n = ll->lo.n = 0; + kv_ul_ov_t *idx = &(ll->lo); + gen_gl_aln(olist, uref, idx); + // fprintf(stderr, "0-[M::%s] idx->n::%lu\n", __func__, (uint64_t)idx->n); + if(idx->n == 0) return 0; + // fprintf(stderr, "(beg0) [M::%s::tid:%ld] debug_i:%ld, qlen:%ld, # cis:%lu, # trans:%lu\n", __func__, tid, debug_i, qlen, (uint64_t)idx->n, o2); + int64_t max_idx, occ = 0, f = 0; + + kv_resize(ul_ov_t, ll->tk, idx->n); + occ = gl_chain_lin(idx, olist->list, ll->tk.a, uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, qlen, UG_SKIP_N, UG_ITER_N, UG_DIS_N, dp, UG_TRANS_W, 0, NULL, uref->ug, 0); + + if(occ) { + if(ff_chain(idx, qlen, 0.99/**P_CHAIN_COV**/, -1/**G_CHAIN_TRANS_RATE**/, ll->tk.a, NULL, NULL, NULL, diff_ec_ul, winLen, km)) { + f = l2g_res_chain(uref->ug, ll->tk.a+idx->a[idx->n-1].ts, idx->a[idx->n-1].te-idx->a[idx->n-1].ts, &(gdp->swap), -1/**N_GCHAIN_RATE**/); + } + } + // fprintf(stderr, "1-[M::%s] f::%ld\n", __func__, f); + // fprintf(stderr, "(beg1) [M::%s] debug_i:%ld, qlen:%ld\n", __func__, debug_i, qlen); + if(!f) { + // gl_chain_gen(olist, uref, idx, 0, hap, km);///no trans + // l2g_chain(uref, idx, &(gdp->l)); ll->tk.n = 0; + gen_gg_aln(olist, uref, UG_TRANS_W, &(gdp->l)); ll->tk.n = 0; + ///buffer + // kv_resize(uint64_t, ll->srt.a, gdp->l.n); kv_resize(uint64_t, hap->snp_srt, gdp->l.n); + // kv_resize(uint64_t, gdp->v, gdp->l.n); kv_resize(int64_t, gdp->f, gdp->l.n); + // max_idx = hc_gchain1_dp(b->km, uref, uref->ug, &(gdp->l), &(gdp->swap), &(gdp->dst), &(gdp->out), &(gdp->path), rch->rlen, + // uopt, G_CHAIN_BW, diff_ec_ul, -1, ll->srt.a.a, sps, gdp->f.a, hap->snp_srt.a, gdp->v.a); + kv_resize(uint64_t, ll->srt.a, gdp->l.n); + max_idx = gl_chain_graph(b->km, uref, uref->ug, &(gdp->l), &(gdp->swap), &(gdp->dst), &(gdp->out), + &(gdp->path), rch->rlen, uopt, G_CHAIN_BW, diff_ec_ul, ll->srt.a.a, sps, dp, UG_SKIP_GRAPH_N, 0); + if(max_idx >= 0 && gen_max_gchain_adv(b->km, uref, debug_i, sps, &(gdp->l), &(ll->tk), NULL, rch->rlen, P_CHAIN_COV, 0.3/**P_FRAGEMENT_PRIMARY_CHAIN_COV**/, + 0.1/**P_FRAGEMENT_PRIMARY_SECOND_COV**/, PRIMARY_UL_CHAIN_MIN, uref->ug->g, &(gdp->dst_done), &(gdp->out), &(gdp->path), ll->srt.a.a, &(gdp->swap))) { + // print_raw_chains(&(gdp->swap), debug_i); + // f = check_trans_rate_gap(&(gdp->swap), &(ll->tk), olist, hap, uref, diff_ec_ul, winLen, G_CHAIN_TRANS_RATE); + f = 1; + } + } + // if(debug_i == 1756) fprintf(stderr, "[M::%s] ulid:%ld, qlen:%ld, f:%ld\n", __func__, debug_i, qlen, f); + if(f) update_ul_vec_t_ug(uref, rch, &(gdp->swap), debug_i); + // debug_ul_vec_t_chain(km, uref->ug->g, rch, &(gdp->dst_done), &(gdp->out)); + // fprintf(stderr, "(beg3) [M::%s::tid:%ld] debug_i:%ld, qlen:%ld\n", __func__, tid, debug_i, qlen); + return 1; +} + + int64_t comput_err_partial_cigar(int64_t ol, overlap_region *z, int64_t *rk) { int64_t k = 0, err = 0, e = z->x_pos_s+ol, wn = z->w_list.n; (*rk) = -1; @@ -5214,19 +5609,21 @@ int64_t comput_sc_partial_cigar(int64_t sc, int64_t ol, double err_sc_r, overlap ///mode: 0->ug; 1->read int64_t ed_dp_c(overlap_region_alloc *o, kv_ul_ov_t *res, ul_ov_t *ex, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, int64_t qlen, int64_t max_skip, uint64_t *srt, uint64_t *idx, uint64_t *track, double err_sc, -uint64_t mode, All_reads *ridx, ma_ug_t *ug) +uint64_t mode, All_reads *ridx, ma_ug_t *ug, uint32_t need_srt) { if(res->n == 0) return 0; uint32_t li_v, lj_v, rev_n; int64_t mm_ovlp, x, i, j, k, sc, csc, mm_sc, mm_idx, qo, qovl, share, minus_sc, pj, n_skip, wi, werr; ul_ov_t *li = NULL, *lj = NULL, rev_t; - radix_sort_ul_ov_srt_qe(res->a, res->a + res->n); - for (i = 1, j = 0; i <= (int64_t)res->n; i++) { - if (i == (int64_t)res->n || res->a[i].qe != res->a[j].qe) { - if(i - j > 1) radix_sort_ul_ov_srt_qs(res->a+j, res->a+i); - j = i; - } - } + if(need_srt) { + radix_sort_ul_ov_srt_qe(res->a, res->a + res->n); + for (i = 1, j = 0; i <= (int64_t)res->n; i++) { + if (i == (int64_t)res->n || res->a[i].qe != res->a[j].qe) { + if(i - j > 1) radix_sort_ul_ov_srt_qs(res->a+j, res->a+i); + j = i; + } + } + } ///res->a[0].qe: min_qe; res->a[res->n-1].qs: max_qs if(res->a[0].qe == qlen && res->a[res->n-1].qs == 0) {///all alignments are contained for (i = 0; i < (int64_t)res->n; ++i) { @@ -5488,13 +5885,14 @@ const ul_idx_t *uref, double diff_ec_ul, int64_t wl, int64_t ql, const ug_opt_t uint64_t k, nw; ul_ov_t *m, *p; int64_t occ, i, ovlp, idx_n; ll->tk.n = ll->lo.n = 0; kv_ul_ov_t *idx = &(ll->lo); + ks_introsort_or_xe(olist->length, olist->list); regen_ul_ov_t_lst(uref, olist, idx); if(idx->n == 0) return 0; kv_resize(uint64_t, ll->srt.a, idx->n); kv_resize(uint64_t, *sps, idx->n); kv_resize(ul_ov_t, ll->tk, idx->n); - occ = ed_dp_c(olist, idx, ll->tk.a, uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, ql, 75, dumy->overlapID, ll->srt.a.a, sps->a, ERROR_RATE, 0, NULL, uref->ug); + occ = ed_dp_c(olist, idx, ll->tk.a, uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, ql, 75, dumy->overlapID, ll->srt.a.a, sps->a, ERROR_RATE, 0, NULL, uref->ug, 0); if((!occ) || (!idx->n)) return 0; idx_n = idx->n; p = &(idx->a[idx_n-1]); if(idx_n <= 1) {//one chain; nothing to do @@ -5529,73 +5927,6 @@ const ul_idx_t *uref, double diff_ec_ul, int64_t wl, int64_t ql, const ug_opt_t return 1; } -void convert_ul_ov_t(ul_ov_t *des, overlap_region *src, const ul_idx_t *uref) -{ - des->qn = (uint32_t)-1; des->qs = src->x_pos_s; des->qe = src->x_pos_e+1; - des->tn = src->y_id; des->el = 1; des->rev = src->y_pos_strand; - des->sec = src->non_homopolymer_errors; - if(des->rev) { - des->ts = uref->ug->u.a[des->tn].len - (src->y_pos_e+1); - des->te = uref->ug->u.a[des->tn].len - src->y_pos_s; - } else { - des->ts = src->y_pos_s; - des->te = src->y_pos_e+1; - } -} - -uint64_t check_connect_ug(const ul_idx_t *uref, uint32_t v, uint32_t w, int64_t bw, double diff_ec_ul, int64_t dq) -{ - const asg_t *g = uref?uref->ug->g:NULL; int64_t dt = -1; - uint32_t nv = asg_arc_n(g, v), i; asg_arc_t *av = asg_arc_a(g, v); - for (i = 0; i < nv; i++) { - if(av[i].del || av[i].v != w) continue; - dt = av[i].ol; - break; - } - if(dt < 0) return 0; - int64_t diff = (dq>dt? dq-dt:dt-dq), mm = MAX(dq, dt); - mm *= diff_ec_ul; if(mm < bw) mm = bw; - if(diff <= mm) return 1; - return 0; -} - -uint64_t check_connect_rg(const ul_idx_t *uref, const ug_opt_t *uopt, uint32_t uv, uint32_t uw, int64_t bw, double diff_ec_ul, int64_t dq) -{ - int64_t dt = -1; - if(uref->ug->u.a[uv>>1].circ || uref->ug->u.a[uw>>1].circ) return 0; - uint32_t rv = (uref->ug->u.a[uv>>1].a[(uv&1)?(0):(uref->ug->u.a[uv>>1].n-1)]>>32)^(uv&1); - uint32_t rw = (uref->ug->u.a[uw>>1].a[(uw&1)?(uref->ug->u.a[uw>>1].n-1):(0)]>>32)^(uw&1); - ma_hit_t_alloc* src = uopt->sources; - int64_t min_ovlp = uopt->min_ovlp; - int64_t max_hang = uopt->max_hang; - uint64_t z, qn, tn, x = rv>>1; int32_t r = 1; asg_arc_t e; - for (z = 0; z < src[x].length; z++) { - qn = Get_qn(src[x].buffer[z]); tn = Get_tn(src[x].buffer[z]); - if(tn != (rw>>1)) continue; - r = ma_hit2arc(&(src[x].buffer[z]), Get_READ_LENGTH(R_INF, qn), Get_READ_LENGTH(R_INF, tn), max_hang, asm_opt.max_hang_rate, min_ovlp, &e); - if(r < 0) continue; - if((e.ul>>32) != rv || e.v != rw) continue; - dt = e.ol; - break; - } - if(dt < 0) return 0; - int64_t diff = (dq>dt? dq-dt:dt-dq), mm = MAX(dq, dt); - mm *= diff_ec_ul; if(mm < bw) mm = bw; - if(diff <= mm) return 1; - return 0; -} - -uint32_t govlp_check(const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, ul_ov_t *li, ul_ov_t *lj) -{ - int64_t qo = infer_rovlp(li, lj, NULL, NULL, /**ridx**/NULL, uref->ug); ///overlap length in query (UL read) - - // fprintf(stderr, "+++[M::%s::utg%.6dl->utg%.6dl] qo::%ld\n", __func__, (int32_t)li->tn+1, (int32_t)lj->tn+1, qo); - if(check_connect_ug(uref, ((li->tn<<1)|li->rev)^1, ((lj->tn<<1)|lj->rev)^1, bw, diff_ec_ul, qo)) return 1; - // fprintf(stderr, "[M::%s::] check_connect_ug fail\n", __func__); - if(check_connect_rg(uref, uopt, ((li->tn<<1)|li->rev)^1, ((lj->tn<<1)|lj->rev)^1, bw, diff_ec_ul, qo)) return 1; - // fprintf(stderr, "[M::%s::] check_connect_rg fail\n", __func__); - return 0; -} uint64_t gen_shared_trace(overlap_region_alloc* ol, uint64_t *id_a, uint64_t id_n, uint64_t s, uint64_t e, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, int64_t *is_srt, kv_ul_ov_t *res)///[s, e] @@ -6049,7 +6380,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call ha_ovec_buf_t *b = s->hab[tid]; 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; + uint32_t high_occ = 2, phase = 1, k; asg64_v b0, b1, b2; overlap_region *aux_o = NULL; // uint64_t align = 0; @@ -6084,8 +6415,8 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // &b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul, winLen, NULL); // memset(&b->self_read, 0, sizeof(b->self_read)); - 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, aux_o, s->opt->diff_ec_ul, winLen, NULL, s->id+i, s->opt->k, NULL); + ul_lalign(&b->olist, &b->clist, s->uu, s->uopt, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, + &b->correct, &b->exz, &b->hap, &b->r_buf, aux_o, s->opt->diff_ec_ul, winLen, NULL, s->id+i, s->opt->k, &(s->sps[tid]), NULL); // ul_lalign_old_ed(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, // &b->correct, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, 1, NULL); ton = b->olist.length;//all alignments pass similary check @@ -6104,12 +6435,20 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call 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, MAX_LGAP(s->len[i]), s->opt->diff_ec_ul); copy_asg_arr(b->hap.snp_srt, b0); copy_asg_arr(s->sps[tid], b1); copy_asg_arr(b->r_buf.a, 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, aux_o, s->opt->diff_ec_ul, winLen, &(bl->lo), s->id+i, s->opt->k, NULL); + ul_lalign(&b->olist, &b->clist, s->uu, s->uopt, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, + &b->correct, &b->exz, &b->hap, &b->r_buf, aux_o, s->opt->diff_ec_ul, winLen, &(bl->lo), s->id+i, s->opt->k, &(s->sps[tid]), NULL); // ul_lalign_old_ed(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, // &b->correct, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, 0, NULL); } } + ///recover alignments + for (k = b->olist.length; k < ton; k++) { + b->olist.list[k].align_length = 0; b->olist.list[k].is_match = 2; + b->olist.list[k].overlapLen = b->olist.list[k].x_pos_e+1-b->olist.list[k].x_pos_s; + } + b->olist.length = ton; + + gl_chain(s->buf[tid], &(UL_INF.a[s->id+i]), &b->olist, &(b->clist.chainDP), &b->hap, &(s->sps[tid]), bl, &(s->gdp[tid]), s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, s->id+i, tid, NULL); // exit(1); // uint64_t k; @@ -6120,7 +6459,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // gl_chain_refine(&b->olist, &b->correct, &b->hap, bl, s->uu, s->opt->diff_ec_ul, winLen, s->len[i], km); - gl_chain_refine_advance_combine(s->buf[tid], &(UL_INF.a[s->id+i]), &b->olist, &b->correct, &b->hap, &(s->sps[tid]), bl, &(s->gdp[tid]), s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, s->id+i, tid, NULL); + // gl_chain_refine_advance_combine(s->buf[tid], &(UL_INF.a[s->id+i]), &b->olist, &b->correct, &b->hap, &(s->sps[tid]), bl, &(s->gdp[tid]), s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, s->id+i, tid, NULL); // return; // b->num_read_base += b->self_read.length; // b->num_correct_base += b->correct.corrected_base; @@ -7264,60 +7603,6 @@ void dump_linear_chain(ma_ug_t *ug, kv_ul_ov_t *lidx, kv_ul_ov_t *autom, vec_mg_ } } -int64_t find_mg_lchain_max(int64_t n, const mg_lchain_t *a, int32_t x) -{ - int64_t s = 0, e = n; - if (n == 0) return -1; - if (a[n-1].qe < x) return n - 1; - if (a[0].qe >= x) return -1; - while (e > s) { // TODO: finish this block - int64_t m = s + (e - s) / 2; - if (a[m].qe >= x) e = m; - else s = m + 1; - } - assert(s == e); - return s; -} - - -int64_t hc_target_len(asg_t *g, mg_lchain_t *s, mg_lchain_t *e) -{ - // int64_t ql = s->qe - e->qe, tp, tm; - // if((s->v^1)&1) tp = g->seq[s->v>>1].len - s->re; - // else tp = s->rs; - - // if((e->v^1)&1) tm = g->seq[e->v>>1].len - e->re; - // else tm = e->rs; - int64_t ql = (int64_t)s->qs - (int64_t)e->qe, tp, tm; - int64_t sts, ete; - sts = (s->v&1)?g->seq[s->v>>1].len-s->re:s->rs; tp = g->seq[s->v>>1].len - sts; - ete = (e->v&1)?g->seq[e->v>>1].len-e->rs:e->re; tm = g->seq[e->v>>1].len - ete; - // fprintf(stderr, "[M::%s::] ql:%ld, tp:%ld, tm:%ld, sts:%ld, ete:%ld\n", __func__, ql, tp, tm, sts, ete); - return ql + tp - tm; -} - -inline int32_t cal_gchain_sc(const mg_path_dst_t *dj, const mg_lchain_t *li, const mg_lchain_t *lc, int64_t *f, int64_t b_w, float diff_thre, float chn_pen_gap) -{ - // const mg_lchain_t *lj; - int32_t gap, sc; - float lin_pen, log_pen; - if (dj->n_path == 0) return INT32_MIN; - gap = dj->dist - dj->target_dist; - // lj = &lc[dj->meta]; - if (gap < 0) gap = -gap; - if ((gap > ((dj->target_dist)*diff_thre)) && (gap > b_w)) return INT32_MIN; - // if (lj->qe <= li->qs) sc = li->score; - // else sc = (int32_t)((double)(li->qe - lj->qe) / (li->qe - li->qs) * li->score + .499); // dealing with overlap on query - sc = li->score; - //sc += dj->mlen; // TODO: is this line the right thing to do? - // if (dj->is_0) sc += ref_bonus; - lin_pen = chn_pen_gap * (float)gap; - log_pen = gap >= 2? mg_log2(gap) : 0.0f; - sc -= (int32_t)(lin_pen + log_pen); - sc += f[dj->meta]; - return sc; -} - void set_trans_arr(uint64_t *trans, vec_sp_node_t *out, int64_t idx) { int64_t i; @@ -7813,12 +8098,7 @@ int64_t *n_u_, int64_t *n_v_) return n_u; } -void set_ul_ov_t_by_mg_lchain_t(ul_ov_t *u, mg_lchain_t *l) -{ - u->tn = l->v>>1; u->rev = (l->v&1); - u->ts = l->rs; u->te = l->re; - u->qs = l->qs; u->qe = l->qe; -} + uint64_t primary_chain_check(uint64_t *idx, int64_t idx_n, mg_lchain_t *a) { diff --git a/inter.h b/inter.h index 8adc7c2..09ee587 100644 --- a/inter.h +++ b/inter.h @@ -13,6 +13,11 @@ #define G_CHAIN_GAP 0.1 #define UG_SKIP 5 #define RG_SKIP 25 +#define UG_SKIP_GRAPH_N 50 +#define UG_SKIP_N 100 +#define UG_ITER_N 5000 +#define UG_DIS_N 50000 +#define UG_TRANS_W 2 #define G_CHAIN_TRANS_RATE 0.25 #define G_CHAIN_TRANS_WEIGHT -1 #define G_CHAIN_INDEL 128