diff --git a/gfa_ut.cpp b/gfa_ut.cpp index a909f67..c3595ef 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -5467,7 +5467,8 @@ void gen_integer_normalize(ul_resolve_t *uidx) void clip_integer_chimeric(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul2ul_idx_t *ul2, integer_t *buf, int64_t min_dp) { - uint64_t k, is_srt = 0, is_del = 0; assert(o->id == qid); + uint64_t k, is_srt = 0, is_del = 0, rid, zs, ze; assert(o->id == qid); + ul_str_t *str = &(uidx->pstr.str.a[qid]); int64_t z, z_n; for (k = 0; k < o->cn; k++) { if(o->a[k].is_del) break; if(k > 0 && o->a[k].hid < o->a[k-1].hid) break; @@ -5534,6 +5535,68 @@ void clip_integer_chimeric(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul if(is_left && is_right && is_middle) is_del = 1; } + if(!is_del) { + for (k = 0; k < str->cn; k++) { + rid = (((uint32_t)str->a[k])>>1); + if(!IF_HOM(rid, *(uidx->bub))) break; + } + + if(k < str->cn) {///at least a hom node covered + for (k = 0, buf->u.n = 0; k < o->cn; k++) { + z = o->a[k].qs_k; z_n = o->a[k].qe_k; + for (; z < z_n; z++) { + rid = (((uint32_t)str->a[z])>>1); + if(!IF_HOM(rid, *(uidx->bub))) break; + } + if(z >= z_n) continue; + zs = z; + + for (z = z_n - 1; z >= (int64_t)zs; z--) { + rid = (((uint32_t)str->a[z])>>1); + if(!IF_HOM(rid, *(uidx->bub))) break; + } + ze = z + 1; + assert(zs < ze); + + kv_push(uint64_t, buf->u, (zs<<1)); + kv_push(uint64_t, buf->u, (ze<<1)|1); + + } + radix_sort_srt64(buf->u.a, buf->u.a + buf->u.n); + + b_n = buf->u.n; + for (k = 0, dp = 0, start = 0; k < b_n; ++k) { + old_dp = dp; + ///if a[j] is qe + if (buf->u.a[k]&1) --dp; + else ++dp; + + if (old_dp < min_dp && dp >= min_dp) {///old_dp < dp, b.a[j] is qs + start = buf->u.a[k]>>1; + } else if (old_dp >= min_dp && dp < min_dp) {///old_dp > min_dp, b.a[j] is qe + end = buf->u.a[k]>>1; + kv_push(uint64_t, buf->u, ((start<<32)|(end))); + } + } + + is_del = 0; + if(buf->u.n == b_n) { + is_del = 1; + } else { + uint32_t is_left = 0, is_right = 0, is_middle = 0; + for (k = b_n; k < buf->u.n; ++k) { + start = buf->u.a[k]>>32; end = (uint32_t)buf->u.a[k]; + if(start == 0) is_left = 1; + else is_middle = 1; + if(end == uidx->idx->a[qid].rlen) is_right = 1; + else is_middle = 1; + } + if(is_left == 0 && is_right == 0) is_del = 1; + if(is_left && is_right && is_middle) is_del = 1; + } + } + } + if(is_del) { o->is_del = 1; for (k = 0; k < o->cn; k++) o->a[k].is_del = 1; @@ -6396,8 +6459,8 @@ inline void get_iug_u_raw_occ(ul_resolve_t *uidx, uint32_t id, uint32_t *ul_occ, void ma_integer_ug_print0(const ma_ug_t *ug, ul_resolve_t *uidx, int print_seq, const char* prefix, FILE *fp) { - uint32_t i, j, l, x; ma_utg_t *p, *s; - char name[32]; + uint32_t i, j, l, x; ma_utg_t *p, *s; ul2ul_idx_t *idx = &(uidx->uovl); + char name[32]; uinfo_srt_warp_t *seq; for (i = 0; i < ug->u.n; ++i) { // the Segment lines in GFA p = &ug->u.a[i]; if(p->m == 0) continue; @@ -6423,6 +6486,14 @@ void ma_integer_ug_print0(const ma_ug_t *ug, ul_resolve_t *uidx, int print_seq, } l += (uint32_t)p->a[j]; } + + seq = &(idx->cc.iug_a[i]); + for (j = 0; j < seq->n; j++) { + x = seq->a[j].v>>1; s = &(uidx->init_ug->u.a[x]); + fprintf(fp, "U\t%s\t%c\tutg%.6d%c\t%d\t%d\tHG:A:*\n", + name, "+-"[seq->a[j].v&1], x + 1, "lc"[s->circ], seq->a[j].s, seq->a[j].e); + } + } if(ug->g) @@ -7859,6 +7930,7 @@ static void worker_update_ul_arc_supports(void *data, long i, int tid) // callba (*x) |= (w_v<<32); } + uint32_t ulg_arc_cut_supports(ul_resolve_t *uidx, ma_ug_t *ug, int32_t max_ext, uint32_t max_ext_hifi, float len_rat, uint32_t is_trio, uint32_t topo_level, uint32_t skip_hom, uint32_t *max_drop_len, asg64_v *in, asg64_v *ib) { @@ -7983,7 +8055,7 @@ uint64_t ulg_bub_pop_cut_aux(ul_resolve_t *uidx, ma_ug_t *ug, uint32_t v0, buf_t n_ext += ul_occ; kv_push(uint64_t, *b, x->b.a[i]); } - if(n_ext < max_ext) { + if(n_ext <= max_ext) { if(get_remove_hifi_occ(uidx, max_ext_hifi, b->a + bn, b->n - bn, ub, NULL)) r = 1; } b->n = bn; @@ -8223,6 +8295,229 @@ uint64_t ulg_pop_bubble(ul_resolve_t *uidx, ma_ug_t *ug, uint64_t* i_max_dist, u return n_pop; } +/** +void infer_reliable_regions(ul_resolve_t *uidx) +{ + ul2ul_idx_t *idx = &(uidx->uovl); ma_ug_t *iug = idx->i_ug; + uint64_t k, z, v; ma_utg_t *iu; uinfo_srt_warp_t *seq; + for (k = 0; k < iug->u.n; k++) { + seq = &(idx->cc.iug_a[k]); + for (z = 0; z < seq->n; z++) { + + } + + } + +} + +void fill_u2g(ul_resolve_t *uidx) +{ + renew_ul2_utg(uidx); + infer_reliable_regions(uidx); +} +**/ + +uint64_t get_ug_integer_seq_occ(ul_resolve_t *uidx, uint64_t *u_a, uint64_t u_n, asg64_v *b) +{ + uint32_t bn = b->n, k, occ; ma_ug_t *raw = uidx->l1_ug; + gen_ug_integer_seq_on_fly(uidx, u_a, u_n, b); + for (k = bn, occ = 0; k < b->n; k++) occ += raw->u.a[b->a[k]>>1].n; + b->n = bn; + return occ; +} + +uint32_t ul_occ_check(ul_resolve_t *uidx, ma_ug_t *ug, uint32_t qocc_ul, uint32_t qocc_hifi, uint32_t tv, float occ_rate, asg64_v *b, asg64_v *ub) +{ + uint32_t bn = b->n, z, tocc_ul, tocc_hifi, ul; + get_ul_path_info(uidx, ug, tv, NULL, NULL, NULL, NULL, NULL, b); + for (z = bn, tocc_ul = 0; z < b->n; z++) { + get_iug_u_raw_occ(uidx, b->a[z]>>1, &ul, NULL); tocc_ul += ul; + } + tocc_hifi = get_ug_integer_seq_occ(uidx, b->a + bn, b->n - bn, ub); + b->n = bn; + if((qocc_ul <= (tocc_ul*occ_rate)) && (qocc_hifi <= (tocc_hifi*occ_rate))) return 1; + return 0; +} + +uint32_t ul_homo_path_check(ul_resolve_t *uidx, ma_ug_t *ug, uint32_t v, uint32_t w, uint32_t raw_ug_occ, float match_rate, asg64_v *b, asg64_v *rb) +{ + bubble_type *bub = uidx->bub; + uint64_t k, l, z, bn = b->n, vn = 0, wn = 0, *va, *wa, rbn = rb->n, *rva, *rwa, rvn, rwn, rid, x; + get_ul_path_info(uidx, ug, v, NULL, NULL, NULL, NULL, NULL, b); vn = b->n - bn; + get_ul_path_info(uidx, ug, w, NULL, NULL, NULL, NULL, NULL, b); wn = b->n - bn - vn; + va = b->a + bn; wa = b->a + bn + vn; + + gen_ug_integer_seq_on_fly(uidx, va, vn, rb); rvn = rb->n - rbn; + gen_ug_integer_seq_on_fly(uidx, wa, wn, rb); rwn = rb->n - rbn - rvn; + rva = rb->a + rbn; rwa = rb->a + rbn + rvn; + if(rvn > raw_ug_occ) rvn = raw_ug_occ; + if(rwn > raw_ug_occ) rwn = raw_ug_occ; + + b->n = bn; + for (k = 0; k < rvn; k++) { + rid = rva[k]>>1; + if(IF_BUB(rid, *bub)) { + x = bub->index[rid]; x |= ((uint64_t)(0x80000000)); x <<= 32; + } else { + x = rid; x <<= 32; + } + kv_push(uint64_t, *b, x); + } + + for (k = 0; k < rwn; k++) { + rid = rwa[k]>>1; + if(IF_BUB(rid, *bub)) { + x = bub->index[rid]; x |= ((uint64_t)(0x80000000)); x <<= 32; + } else { + x = rid; x <<= 32; + } + x |= 1; kv_push(uint64_t, *b, x); + } + radix_sort_srt64(b->a + bn, b->a + b->n); + + uint64_t o[2]; + for (l = bn, k = bn + 1, o[0] = o[1] = 0; k <= b->n; k++) { + if (k == b->n || (b->a[k]>>32) != (b->a[l]>>32)) { + if((k - l > 1) && (((uint32_t)b->a[l]) != ((uint32_t)b->a[k-1]))) { + for (z = l; z < k; z++) o[(uint32_t)b->a[z]]++; + } + l = k; + } + } + b->n = bn; rb->n = rbn; + + if(o[0] >= (rvn*match_rate)) return 1; + if(o[1] >= (rwn*match_rate)) return 1; + + return 0; +} + +///small_occ_rate = 0.15; len_rat = 1.5 +uint32_t ulg_arc_cut_z(ul_resolve_t *uidx, ma_ug_t *ug, uint32_t max_ext, uint32_t max_ext_hifi, +float len_rat, float small_occ_rate, uint32_t raw_ug_occ, float raw_match_rate, uint32_t is_trio, +uint32_t skip_hom, uint32_t *max_drop_len, asg64_v *in, asg64_v *ib) +{ + asg64_v tx = {0,0,0}, tb = {0,0,0}, *b = NULL, *ub = NULL; asg_t *g = ug->g; + uint32_t v, w, wt, z, i, k, kv, nv, kw, kwt, nw, cnt = 0, n_vtx = g->n_seq<<1, ul, vp[2], wp[2]; + asg_arc_t *av, *aw, *ve, *we, *wte; uint64_t w_q, w_t, pb, raw_ul, raw_hifi; + + b = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); + for (v = b->n = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + if (nv < 2) continue; + for (i = kv = 0; i < nv && kv <= 2; ++i) { + if(av[i].del) continue; kv++; + } + if(kv != 2) continue; + for (i = 0; i < nv; ++i) { + if(av[i].del) continue; + kw = get_arcs(ug->g, av[i].v^1, NULL, 0); + if(kw == 2) { + kv_push(uint64_t, *b, ((uint64_t)(av-g->arc+i))); + } else if(kw == 1) { + ub->n = 0; + if(get_ul_path_info(uidx, ug, av[i].v, NULL, NULL, NULL, NULL, NULL, ub)==TWO_INPUT) { + for (z = raw_ul = 0; z < ub->n; z++) { + get_iug_u_raw_occ(uidx, ub->a[z]>>1, &ul, NULL); raw_ul += ul; + } + kv_push(uint64_t, *b, ((raw_ul<<32)|((uint64_t)(av-g->arc+i)))); + } + } + } + } + + radix_sort_srt64(b->a, b->a + b->n); + for (k = 0; k < b->n; k++) { + if(g->arc[(uint32_t)b->a[k]].del) continue; + v = g->arc[(uint32_t)b->a[k]].ul>>32; w = g->arc[(uint32_t)b->a[k]].v^1; + if(g->seq[v>>1].del || g->seq[w>>1].del) continue; + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + nw = asg_arc_n(g, w); aw = asg_arc_a(g, w); + if(nv <= 1 && nw <= 1) continue; + + vp[0] = v^1; vp[1] = (uint32_t)-1; + ve = &(g->arc[(uint32_t)b->a[k]]); + for (i = kv = 0; i < nv; ++i) { + if(av[i].del) continue; + if(av[i].v != (w^1)) vp[1] = av[i].v; + kv++; + } + + wp[0] = w^1; wp[1] = (uint32_t)-1; + for (i = kw = 0, we = NULL; i < nw; ++i) { + if (aw[i].del) continue; + if (aw[i].v == (v^1)) we = &(aw[i]); + else wp[1] = aw[i].v; + kw++; + } + if(kv <= 1 && kw <= 1) continue; + if(kv != 2 || kw > 2) continue; + raw_ul = 0; wt = w; wte = we; kwt = kw; pb = b->n; ub->n = 0; + if(kw == 1) { + if(get_ul_path_info(uidx, ug, w^1, &wt, NULL, NULL, NULL, NULL, b)==TWO_INPUT) { + for (z = pb, raw_ul = 0; z < b->n; z++) { + get_iug_u_raw_occ(uidx, b->a[z]>>1, &ul, NULL); raw_ul += ul; + } + ul = wt; get_arcs(ug->g, wt, &wt, 1); wt = ug->g->arc[wt].v^1; + wp[0] = wt^1; wp[1] = (uint32_t)-1; + nw = asg_arc_n(g, wt); aw = asg_arc_a(g, wt); + for (i = kwt = 0; i < nw; ++i) { + if (aw[i].del) continue; + if (aw[i].v == (ul^1)) wte = &(aw[i]); + else wp[1] = aw[i].v; + kwt++; + } + nw = asg_arc_n(g, w); aw = asg_arc_a(g, w); + } else { + b->n = pb; + continue; + } + } + + assert(kv == 2 && kwt == 2); + if(raw_ul <= max_ext && get_remove_hifi_occ(uidx, max_ext_hifi, b->a + pb, b->n - pb, ub, NULL)) { + raw_hifi = get_ug_integer_seq_occ(uidx, b->a + pb, b->n - pb, ub); + + b->n = pb; + if(raw_ul > 0) { + if(!ul_occ_check(uidx, ug, raw_ul, raw_hifi, vp[0], small_occ_rate, b, ub)) continue; + if(!ul_occ_check(uidx, ug, raw_ul, raw_hifi, vp[1], small_occ_rate, b, ub)) continue; + if(!ul_occ_check(uidx, ug, raw_ul, raw_hifi, wp[0], small_occ_rate, b, ub)) continue; + if(!ul_occ_check(uidx, ug, raw_ul, raw_hifi, wp[1], small_occ_rate, b, ub)) continue; + } + + + pb = b->n; ub->n = 0; + get_ul_arc_supports(uidx, ve, b, ub, skip_hom, &w_q, &w_t); + b->n = pb; ub->n = 0; + if((w_q == (uint64_t)-1) || (w_q > w_t*len_rat)) continue; + + pb = b->n; ub->n = 0; + get_ul_arc_supports(uidx, wte, b, ub, skip_hom, &w_q, &w_t); + b->n = pb; ub->n = 0; + if((w_q == (uint64_t)-1) || (w_q > w_t*len_rat)) continue; + + if(!ul_homo_path_check(uidx, ug, vp[0], wp[1], raw_ug_occ, raw_match_rate, b, ub)) continue; + if(!ul_homo_path_check(uidx, ug, wp[0], vp[1], raw_ug_occ, raw_match_rate, b, ub)) continue; + + if(kw == 1) { + get_ul_path_info(uidx, ug, w^1, &wt, NULL, NULL, NULL, NULL, b); + for (z = pb; z < b->n; z++) ulg_seq_del(ug, (b->a[z]>>1)); + } + ve->del = we->del = 1; cnt++; + } + + b->n = pb; + } + + + + if(!in) free(tx.a); if(!ib) free(tb.a); + if (cnt > 0) asg_cleanup(g); + return cnt; +} + void u2g_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt) { ul2ul_idx_t *idx = &(uidx->uovl); asg64_v bu = {0,0,0}, uu = {0,0,0}; @@ -8263,6 +8558,10 @@ void u2g_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt) ulg_pop_bubble(uidx, iug, NULL, ((int64_t)0x7fffffff), ulopt->max_tip_hifi, 1, &bu, &uu); + while(ulg_arc_cut_z(uidx, iug, ((int64_t)0x7fffffff), ulopt->max_tip_hifi, 1.5, 0.15, 100, 0.8, ulopt->is_trio, 1, NULL, &bu, &uu)); + + // fill_u2g(uidx); + // while (cnt) { // for (i = cnt = 0, mm_tip = ulopt->max_tip, drop = ulopt->min_ovlp_drop_ratio; i < ulopt->clean_round; i++, drop += step) { // if(drop > ulopt->max_ovlp_drop_ratio) drop = ulopt->max_ovlp_drop_ratio; @@ -8479,6 +8778,10 @@ void renew_u2g_cov(ul_resolve_t *uidx) { ul2ul_idx_t *idx = &(uidx->uovl); uinfo_srt_warp_t *x; ma_ug_t *i_ug = idx->i_ug, *raw = uidx->l1_ug; uint64_t k, z, iug_occ, m, l, *a, a_n; + free(idx->cc.uc); free(idx->cc.hc); free(idx->cc.raw_uc); + free(idx->cc.iug_a); free(idx->cc.iug_idx); free(idx->cc.iug_b); + memset(&(idx->cc), 0, sizeof(idx->cc)); + MALLOC(idx->cc.uc, i_ug->u.n); MALLOC(idx->cc.hc, i_ug->u.n); MALLOC(idx->cc.raw_uc, i_ug->u.n); @@ -8555,6 +8858,7 @@ void renew_u2g_bg(ul_resolve_t *uidx) ma_ug_t *iug = idx->i_ug, *raw = uidx->l1_ug; bubble_type *bub = uidx->bub; uint64_t k, z, l, *raw_a, raw_n, raw_id, iug_id, iug_off, v, w, nv, nw, n_vtx; uinfo_srt_warp_t *x; asg64_v buf = {0,0,0}; int64_t s, s_n; asg_arc_t *p, *av, *aw; + asg_destroy(bg->bg); free(bg->w_n); free(bg->a_n); memset(bg, 0, sizeof((*bg))); bg->bg = asg_init(); bg->bg->n_seq = 0; bg->bg->m_seq = raw->g->n_seq; MALLOC(bg->bg->seq, bg->bg->m_seq); @@ -8618,19 +8922,11 @@ void renew_u2g_bg(ul_resolve_t *uidx) av = asg_arc_a(bg->bg, v); nv = asg_arc_n(bg->bg, v); for (k = 0; k < nv && k < 2; k++) { v_occ[k] = av[k].ol; v_occ[k] <<= 32; v_occ[k] += av[k].v; - // if((v>>1) == 172 && (w>>1) == 168) { - // fprintf(stderr, "+[M::%s::] k::%lu, v>>1::%lu, av[k].v>>1:%u, av[k].ol::%u\n", - // __func__, k, v>>1, av[k].v>>1, av[k].ol); - // } } aw = asg_arc_a(bg->bg, (w^1)); nw = asg_arc_n(bg->bg, (w^1)); for (k = 0; k < nw && k < 2; k++) { w_occ[k] = aw[k].ol; w_occ[k] <<= 32; w_occ[k] += aw[k].v; - // if((v>>1) == 172 && (w>>1) == 168) { - // fprintf(stderr, "+[M::%s::] k::%lu, w>>1::%lu, aw[k].v>>1:%u, aw[k].ol::%u\n", - // __func__, k, w>>1, aw[k].v>>1, aw[k].ol); - // } } if((((uint32_t)v_occ[0]) == w) && (((uint32_t)w_occ[0]) == (v^1))) { @@ -8669,10 +8965,6 @@ void renew_u2g_bg(ul_resolve_t *uidx) assert(l != bg_unavailable); if(l == bg_wrong) bg->w_n[k]++; if(l == bg_ambiguous) bg->a_n[k]++; - // if(k == 79 || k == 80) { - // fprintf(stderr, "+[M::%s::] z::%lu, nv:%lu, pv>>1::%lu, cv>>1::%lu, l::%lu\n", - // __func__, z, nv, v_occ[(nv-2)&1]>>1, v_occ[(nv-1)&1]>>1, l); - // } } // fprintf(stderr, "-[M::%s::] k::%lu, x->n::%u, w_n[k]::%u, a_n[k]::%u\n", // __func__, k, (uint32_t)x->n, bg->w_n[k], bg->a_n[k]); @@ -8715,6 +9007,15 @@ void update_ul_tra_idx_t(ul_resolve_t *uidx) kt_for(uidx->str_b.n_thread, worker_update_ul_tra_idx, uidx, n_vtx);///all ul + ug } **/ +void renew_ul2_utg(ul_resolve_t *uidx) +{ + ul2ul_idx_t *z = &(uidx->uovl); + renew_utg(&(z->i_ug), z->i_g, NULL); + free(z->i_ug->g->seq_vis); CALLOC(z->i_ug->g->seq_vis, z->i_ug->g->n_seq*2); + renew_u2g_cov(uidx); + renew_u2g_bg(uidx); +} + ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt) { @@ -8758,7 +9059,8 @@ ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt) // output_integer_graph(uidx, z->i_ug, asm_opt.output_file_name); u2g_clean(uidx, ulopt); - // renew_utg(&(z->i_ug), z->i_g, NULL); + // renew_ul2_utg(uidx); + output_integer_graph(uidx, z->i_ug, asm_opt.output_file_name); return z; }