From 33afa25398b2ec5eb8899e393504fc1c5ecaabf3 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Sun, 24 Jul 2022 18:13:16 -0400 Subject: [PATCH] UL graph done --- Overlaps.cpp | 18 +- Overlaps.h | 4 +- gfa_ut.cpp | 3582 +++++++++++++++++++++++++++++++++++++++++++++++++- gfa_ut.h | 3 +- inter.cpp | 15 +- 5 files changed, 3581 insertions(+), 41 deletions(-) diff --git a/Overlaps.cpp b/Overlaps.cpp index b9878c8..991dd96 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -11180,7 +11180,6 @@ void debug_contig_end(ma_ug_t *ug, asg_t* read_g, kvec_asg_arc_t_warp* edge) } } -void renew_utg(ma_ug_t **ug, asg_t* read_g, kvec_asg_arc_t_warp* edge); void push_sub_unitig(ma_ug_t *n_ug, ma_utg_t *src_u, asg_t *read_g, kvec_asg_arc_t_warp* edge, uint32_t beg_idx, uint32_t occ) { @@ -14819,14 +14818,14 @@ const char* command) uint32_t print_debug_gfa(asg_t *read_g, ma_ug_t *ug, ma_sub_t* coverage_cut, const char* output_file_name, -ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) +ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp, int is_polish, int is_update_ou, int is_check_alter_lable) { kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); if(ug == NULL) { ug = ma_ug_gen(read_g); - } else { + } else if(is_check_alter_lable) { uint32_t i; for (i = 0; i < ug->u.n; ++i) { @@ -14844,8 +14843,8 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) } } - update_ug_ou(ug, read_g); - ma_ug_seq(ug, read_g, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 0); + if(is_update_ou) update_ug_ou(ug, read_g); + if(is_polish) ma_ug_seq(ug, read_g, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 0); fprintf(stderr, "Writing raw unitig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+25); @@ -16627,7 +16626,7 @@ int min_ovlp, hap_cov_t *cov) resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, cov->is_r_het, (uint32_t)-1, drop_ratio); drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het); - print_debug_gfa(read_g, ug, coverage_cut, "debug_clean_end", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); + // print_debug_gfa(read_g, ug, coverage_cut, "debug_clean_end", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, reverse_sources, 0, 1); if(round > 0) @@ -31549,8 +31548,11 @@ ma_sub_t **coverage_cut_ptr, int debug_g) gen_ug_opt_t(&uopt, sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, &b_mask_t); } - if(asm_opt.ar) ul_realignment_gfa(&uopt, sg); - print_debug_gfa(sg, NULL, coverage_cut, "UL.debug", sources, ruIndex, max_hang_length, mini_overlap_length); + if(asm_opt.ar) { + ul_realignment_gfa(&uopt, sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, + asm_opt.max_short_tip, &b_mask_t, ha_opt_triobin(&asm_opt)); + } + print_debug_gfa(sg, NULL, coverage_cut, "UL.debug", sources, ruIndex, max_hang_length, mini_overlap_length, 0, 0, 0); /** asg_cut_tip(sg, asm_opt.max_short_tip); ///debug_info_of_specfic_node("m64043_200505_112554/8849050/ccs", sg, "inner_1"); diff --git a/Overlaps.h b/Overlaps.h index b750254..56f9d51 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -259,6 +259,7 @@ typedef struct { // kv_ul_ov_t *ov; } ul_idx_t; +#define MA_HT_DOUBLE (-1024) #define MA_HT_INT (-1) #define MA_HT_QCONT (-2) #define MA_HT_TCONT (-3) @@ -1079,7 +1080,7 @@ int asg_arc_del_triangular_directly(asg_t *g, long long min_edge_length, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex); int asg_arc_del_short_diploid_by_exact(asg_t *g, int max_ext, ma_hit_t_alloc* sources); uint32_t print_debug_gfa(asg_t *read_g, ma_ug_t *ug, ma_sub_t* coverage_cut, const char* output_file_name, -ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp); +ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp, int is_polish, int is_update_ou, int is_check_alter_lable); void debug_info_of_specfic_node(const char* name, asg_t *g, R_to_U* ruIndex, const char* command); ma_ug_t *gen_polished_ug(const ug_opt_t *uopt, asg_t *sg); void output_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, @@ -1087,6 +1088,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp); void flat_soma_v(asg_t *sg, ma_hit_t_alloc* sources, R_to_U* ruIndex); void hic_clean(asg_t* read_g); int64_t count_edges_v_w(asg_t *g, uint32_t v, uint32_t w); +void renew_utg(ma_ug_t **ug, asg_t* read_g, kvec_asg_arc_t_warp* edge); #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/gfa_ut.cpp b/gfa_ut.cpp index c805bcf..a909f67 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -22,6 +22,85 @@ KRADIX_SORT_INIT(srt64, uint64_t, generic_key, 8) #define UL_TRAV_HERATE 0.2 #define UL_TRAV_FT_RATE 0.8 + +typedef struct{ + int64_t tipsLen; + float tip_drop_ratio; + int64_t stops_threshold; + float chimeric_rate; + float drop_ratio; + + bub_label_t* b_mask_t; + int64_t clean_round; + double min_ovlp_drop_ratio; + double max_ovlp_drop_ratio; + double hom_check_drop_rate; + int64_t max_tip, max_tip_hifi; + uint32_t is_trio; +}ulg_opt_t; + +typedef struct { + uint32_t hid; + uint32_t qs, qe, ts, te; + uint32_t qs_k, qe_k, ts_k, te_k; + uint8_t is_rev:6, is_del:1, is_ct:1; +} ul2ul_t; + +typedef struct { + ul2ul_t *a; + size_t n, m; + uint32_t id:31, is_del:1; + uint32_t cn; +} ul2ul_item_t; + +typedef struct { + uint32_t v, s, e, n; +} uinfo_srt_t; + +typedef struct { + size_t n, m; + uinfo_srt_t *a; +} uinfo_srt_warp_t; + +typedef struct { + uint32_t *uc, *hc, *raw_uc; + uinfo_srt_warp_t *iug_a; + uint32_t *iug_idx; + uint64_t *iug_b; +} ul_cov_t; + + +typedef struct { + asg_t *bg; + uint32_t *w_n, *a_n; +} ul_bg_t; + +// typedef struct { +// uint32_t v, n, wn; +// } ul_tra_t; + +// typedef struct { +// kvec_t(ul_tra_t) arc; +// kvec_t(uint32_t) idx; +// } ul_tra_idx_t; + +// #define iug_tra_arc_n(z, v) ((z)->idx.a[(v)+1]-(z)->idx.a[(v)]) +// #define iug_tra_arc_a(z, v) ((z)->arc.a + (z)->idx.a[(v)]) + +typedef struct { + ul2ul_item_t *a; + size_t n, m; + uint64_t uln, gn, tot; + uint32_t *item_idx; + asg_t *i_g; ma_ug_t *i_ug; + ul_cov_t cc; + ul_bg_t bg; + asg64_v *iug_tra; +} ul2ul_idx_t; + +#define ul2ul_srt_key(p) ((p).hid) +KRADIX_SORT_INIT(ul2ul_srt, ul2ul_t, ul2ul_srt_key, member_size(ul2ul_t, hid)) + typedef struct { asg_t *g; ma_hit_t_alloc *src; @@ -49,13 +128,10 @@ typedef struct { //s: state, s=0, this edge has not been visited, otherwise, s=1 } uinfo_t; -typedef struct { - uint64_t c; // max count of positive reads - uint32_t i; -} uinfo_srt_t; -#define uinfo_srt_t_c_key(p) ((p).c) -KRADIX_SORT_INIT(uinfo_srt_t_c, uinfo_srt_t, uinfo_srt_t_c_key, member_size(uinfo_srt_t, c)) + +// #define uinfo_srt_t_c_key(p) ((p).se) +// KRADIX_SORT_INIT(uinfo_srt_t_c, uinfo_srt_t, uinfo_srt_t_c_key, member_size(uinfo_srt_t, se)) typedef struct { ///all information for each node @@ -118,7 +194,6 @@ typedef struct{ #define integer_aln_t_vqk_key(x) ((x).tn_rev_qk) KRADIX_SORT_INIT(integer_aln_t_srt, integer_aln_t, integer_aln_t_vqk_key, member_size(integer_aln_t, tn_rev_qk)) - typedef struct { size_t n, m; integer_aln_t *a; @@ -228,6 +303,7 @@ typedef struct { } cul_g_t; typedef struct{ + ug_opt_t *uopt; ma_ug_t *init_ug; ma_ug_t *l0_ug; ma_ug_t *l1_ug; @@ -240,6 +316,7 @@ typedef struct{ ul_str_idx_t pstr; integer_ml_t str_b; cul_g_t *cg; + ul2ul_idx_t uovl; // ul_path_srt_t psrt; }ul_resolve_t; @@ -2597,10 +2674,10 @@ void init_ul_str_idx_t(ul_resolve_t *p) } } -ul_resolve_t *init_ul_resolve_t(asg_t *sg, ma_ug_t *init_ug, bubble_type* bub, all_ul_t *idx, uint8_t *r_het) +ul_resolve_t *init_ul_resolve_t(asg_t *sg, ma_ug_t *init_ug, bubble_type* bub, all_ul_t *idx, ug_opt_t *uopt, uint8_t *r_het) { ul_resolve_t *p = NULL; CALLOC(p, 1); - p->sg = sg; p->init_ug = init_ug; p->bub = bub; p->idx = idx; p->r_het = r_het; + p->sg = sg; p->init_ug = init_ug; p->bub = bub; p->idx = idx; p->r_het = r_het; p->uopt = uopt; p->l1_ug = copy_untig_graph(p->init_ug); init_integer_ml_t(&p->str_b, p, asm_opt.thread_num); // init_ul_path_srt_t(p); @@ -3624,6 +3701,180 @@ uint64_t *cns, uint64_t cns_occ) return 1; } + +void uc_block_qse_cutoff(uint32_t c_ts, uint32_t c_te, uc_block_t *x, uint32_t *nqs, uint32_t *nqe) +{ + uint32_t l, r; + assert(c_ts >= x->ts && c_te <= x->te); + if(!(x->rev)) { + l = c_ts - x->ts; r = x->te - c_te; + } else { + r = c_ts - x->ts; l = x->te - c_te; + } + + (*nqs) = x->qs + get_offset_adjust(l, x->te - x->ts, x->qe - x->qs); + (*nqe) = x->qe - get_offset_adjust(r, x->te - x->ts, x->qe - x->qs); + assert((*nqs) >= x->qs && (*nqe) <= x->qe); +} + +int64_t update_exact_ul_ovlps(all_ul_t *ul_idx, ma_ug_t *ug, ul_str_t *str, integer_aln_t *aln, ul_chain_t *idx, int64_t qid, +integer_t *buf, ul2ul_t *res) +{ + if(idx->e<=idx->s) return 0; + int64_t qk, tk, is_rev, tid, q_end, t_end; + ul_str_t *q_str, *t_str; integer_aln_t *x; + tid = aln[idx->s].tn_rev_qk>>33; is_rev = ((aln[idx->s].tn_rev_qk>>32)&1); + q_str = &(str[qid]); t_str = &(str[tid]); + uint64_t os, oe; uint32_t v, w, q_ns, q_ne, t_ns, t_ne, acc, chained, qbs, qbe, tbs, tbe; + // if(qid == 95 || qid == 36) { + // fprintf(stderr,"[M::%s::qid->%ld::tid->%ld] qs_idx::%u, qe_idx::%u, ts_idx::%u, te_idx::%u\n", __func__, qid, tid, + // (uint32_t)(aln[idx->s].tn_rev_qk), (uint32_t)(aln[idx->e-1].tn_rev_qk), + // aln[idx->s].tk, aln[idx->e-1].tk); + // } + //beg + x = &(aln[idx->s]); + qk = (uint32_t)(x->tn_rev_qk); tk = ((is_rev == 0)? (x->tk):(t_str->cn - x->tk - 1)); + if((qk > 0) && (x->tk > 0)) return 0; + idx->q_sidx = (uint32_t)(x->tn_rev_qk); idx->t_sidx = x->tk; + + //end + x = &(aln[idx->e-1]); + qk = (uint32_t)(x->tn_rev_qk); tk = ((is_rev == 0)? (x->tk):(t_str->cn - x->tk - 1)); + if((((uint32_t)qk + 1) < q_str->cn) && ((x->tk + 1) < t_str->cn)) return 0; + idx->q_eidx = (uint32_t)(x->tn_rev_qk)+1; idx->t_eidx = x->tk+1; + + if((idx->q_eidx - idx->q_sidx) != (idx->t_eidx - idx->t_sidx)) return 0; + if(idx->q_eidx <= idx->q_sidx) return 0; + + qk = idx->q_sidx; q_end = idx->q_eidx; + tk = idx->t_sidx; t_end = idx->t_eidx; + for (; qk < q_end && tk < t_end; qk++, tk++){ + v = ((uint32_t)q_str->a[qk]); + if(!is_rev) w = ((uint32_t)t_str->a[tk]); + else w = ((uint32_t)t_str->a[t_str->cn-tk-1])^1; + if(v != w) break; + } + if(qk!=q_end || tk!=t_end) return 0; + + uc_block_t *qi, *ti; + qk = idx->q_sidx; q_end = idx->q_eidx; + tk = idx->t_sidx; t_end = idx->t_eidx; + + buf->p.n = buf->f.n = buf->o.n = buf->u.n = 0; + kv_resize(int64_t, buf->p, idx->q_eidx - idx->q_sidx); int64_t *p = buf->p.a; + kv_resize(int64_t, buf->f, idx->q_eidx - idx->q_sidx); int64_t *f = buf->f.a; + kv_resize(uint64_t, buf->o, idx->q_eidx - idx->q_sidx); uint64_t *qc = buf->o.a; + kv_resize(uint64_t, buf->u, idx->q_eidx - idx->q_sidx); uint64_t *tc = buf->u.a; + int64_t max_f, max_z, csc, ps, sc, tf, tz, z; + for (acc = 0, chained = 1, tf = tz = -1; qk < q_end && tk < t_end; qk++, tk++) { + qi = &(ul_idx->a[qid].bb.a[q_str->a[qk]>>32]); + if(!is_rev) ti = &(ul_idx->a[tid].bb.a[t_str->a[tk]>>32]); + else ti = &(ul_idx->a[tid].bb.a[t_str->a[t_str->cn-tk-1]>>32]); + assert(qi->hid == ti->hid); + os = MAX(qi->ts, ti->ts); oe = MIN(qi->te, ti->te); + if(oe <= os) continue; + uc_block_qse_cutoff(os, oe, qi, &q_ns, &q_ne); + uc_block_qse_cutoff(os, oe, ti, &t_ns, &t_ne); + qc[acc] = q_ns; qc[acc] <<= 32; qc[acc] |= q_ne; + tc[acc] = t_ns; tc[acc] <<= 32; tc[acc] |= t_ne; + + z = acc; z -= 1; csc = q_ne - q_ns; if(csc <= 0) csc = 1; + max_f = csc; max_z = -1; ps = 0; + if(chained && z >= 0) { + qbs = qc[z]>>32; qbe = (uint32_t)qc[z]; + tbs = tc[z]>>32; tbe = (uint32_t)tc[z]; + if(!is_rev) { + if(qbs <= q_ns && qbe <= q_ne && tbs <= t_ns && tbe <= t_ne) ps = 1; + } else { + if(qbs <= q_ns && qbe <= q_ne && tbs >= t_ns && tbe >= t_ne) ps = 1; + } + if(ps) { + sc = csc + f[z]; + if(sc > max_f) { + max_f = sc; max_z = z; + } + } else { + chained = 0; + } + } + + if((!ps) || (!chained)) { + for (; z >= 0; z--) { + qbs = qc[z]>>32; qbe = (uint32_t)qc[z]; + tbs = tc[z]>>32; tbe = (uint32_t)tc[z]; + ps = 0; + if(!is_rev) { + if(qbs <= q_ns && qbe <= q_ne && tbs <= t_ns && tbe <= t_ne) ps = 1; + } else { + if(qbs <= q_ns && qbe <= q_ne && tbs >= t_ns && tbe >= t_ne) ps = 1; + } + if(!ps) continue; + sc = csc + f[z]; + if(sc > max_f) { + max_f = sc; max_z = z; + } + } + } + + f[acc] = max_f; p[acc] = max_z; + if(tf < max_f) { + tf = max_f; tz = acc; + } + acc++; + } + if(acc <= 0) return 0; + assert(tz >= 0); + res->hid = res->qs = res->qe = res->ts = res->te = (uint32_t)-1; + res->hid = tid; res->is_rev = !!is_rev; res->is_del = 0; res->is_ct = 0; + res->qs_k = idx->q_sidx; res->qe_k = idx->q_eidx; + if(!is_rev) { + res->ts_k = idx->t_sidx; res->te_k = idx->t_eidx; + } else { + res->ts_k = t_str->cn - idx->t_eidx; res->te_k = t_str->cn - idx->t_sidx; + } + + + res->qe = (uint32_t)qc[tz]; + if(!is_rev) res->te = (uint32_t)tc[tz]; + else res->ts = tc[tz]>>32; + + for (z = tz; z >= 0; z = p[z]) { + res->qs = qc[z]>>32; + if(!is_rev) res->ts = tc[z]>>32; + else res->te = (uint32_t)tc[z]; + } + assert(res->qs <= res->qe && res->ts <= res->te); + + + int64_t qs = 0, qe = 0, rs = 0, re = 0, qtail = 0, rtail = 0; + qs = res->qs; qe = res->qe; rs = res->ts; re = res->te; + if(is_rev) { + rs = ul_idx->a[tid].rlen - res->te; re = ul_idx->a[tid].rlen - res->ts; + } + if(qs <= rs) { + rs -= qs; qs = 0; + } else { + qs -= rs; rs = 0; + } + + qtail = ul_idx->a[qid].rlen - qe; rtail = ul_idx->a[tid].rlen - re; + if(qtail <= rtail) { + qe = ul_idx->a[qid].rlen; re += qtail; + } + else { + re = ul_idx->a[tid].rlen; qe += rtail; + } + res->qs = qs; res->qe = qe; + res->ts = rs; res->te = re; + if(is_rev) { + res->ts = ul_idx->a[tid].rlen - re; + res->te = ul_idx->a[tid].rlen - rs; + } + return 1; +} + + + void reset_poa_g_t(poa_g_t *g) { g->seq.n = g->arc.n = g->idx.n = 0; g->update_arc = g->update_seq = g->e_idx.n = 0; @@ -4734,7 +4985,7 @@ void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_ // print_ul_alignment(ug, &UL_INF, 27512, "inner-1"); radix_sort_integer_aln_t_srt(buf->b.a, buf->b.a + buf->b.n); b_n = buf->b.n; buf->sc.n = 0; - for (k = 1, z = 0; k < b_n; k++) { + for (k = 1, z = 0; k <= b_n; k++) { if(k == b_n || (buf->b.a[z].tn_rev_qk>>32) != (buf->b.a[k].tn_rev_qk>>32)) { if(integer_chain(qid, buf->b.a + z, k - z, z, buf, ug, str_idx, uidx->idx, &sc) && sc.v != (uint32_t)-1) { if((buf->sc.n > 0) && ((buf->sc.a[buf->sc.n-1].v>>1) == (sc.v>>1))) { @@ -4752,6 +5003,7 @@ void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_ uint64_t *o, o_n, cns_het, cns_het_occ, ref_cns_occ, corrected = 0; o_n = integer_chain_dp(uidx->bub, buf, str_idx->str.a, buf->b.a, buf->sc.a, buf->sc.n, qid, is_hom, 2, &corrected); assert(o_n <= str->cn); + // if(qid == 2062 || qid == 2093) fprintf(stderr,"[M::%s::] o_n::%lu, str->cn::%u, corrected::%lu\n", __func__, o_n, str->cn, corrected); if(o_n <= 0) return; if(corrected) return; // print_ul_alignment(ug, &UL_INF, 27512, "inner-3"); @@ -4774,7 +5026,7 @@ void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_ // fprintf(stderr, "\n"); // print_integer_seq(ug, str_idx->str.a, qid, 1); // print_aln_seq(ug, uidx->idx, qid, 1); - // print_cns_seq(ug, str, o, o_n); + // if(qid == 2062 || qid == 2093) print_cns_seq(ug, str, o, o_n); // print_ul_alignment(ug, &UL_INF, 27512, "inner-5"); for (k = m = 0; k < buf->sc.n; k++) { @@ -4812,6 +5064,363 @@ void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_ update_raw_integer_seq(&(buf->pg), ug, buf->pg.srt_b.res.a, buf->pg.srt_b.res.n, uidx->idx, str_idx->str.a, qid, buf, buf->sc.a, buf->sc.n); } +ul2ul_item_t *get_ul_ovlp(ul2ul_idx_t *z, uint64_t id, uint64_t is_ul) +{ + uint64_t x = (is_ul?(id):(id+z->uln)); + if(z->item_idx[x] == (uint32_t)-1) return NULL; + return &(z->a[z->item_idx[x]]); +} + +ul2ul_t *get_ul_o(ul2ul_idx_t *idx, uint32_t qid, uint32_t is_q_ul, uint32_t tid, uint32_t is_t_ul) +{ + ul2ul_item_t *z = get_ul_ovlp(idx, qid, is_q_ul); uint64_t k; + if(!is_t_ul) tid += idx->uln; + for (k = 0; k < z->n; k++) { + if(z->a[k].hid == tid) return &(z->a[k]); + } + return NULL; +} + +void integer_gen_ovlp(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, ul2ul_item_t *o, uint32_t ug_offset) +{ + ul_str_idx_t *str_idx = &(uidx->pstr); ma_ug_t *ug = uidx->l1_ug; ul2ul_t res; + ul_str_t *str = &(str_idx->str.a[qid]); integer_aln_t *p; ul_chain_t sc; + uint64_t k, z, *hid_a, hid_n, sck, b_n, m; uint32_t vk, vz; uc_block_t *xi; + o->n = o->cn = 0; o->id = qid; o->is_del = 0; + // if(qid == 95 || qid == 36) print_integer_seq(ug, str_idx->str.a, qid, 1); + if(str->cn < 2) return; + for (k = 0, buf->b.n = 0; k < str->cn; k++) { + xi = &(uidx->idx->a[qid].bb.a[str->a[k]>>32]); + assert(((xi->hid<<1)+xi->rev)==((uint32_t)str->a[k])); + sck = ug_occ_w(xi->ts, xi->te, &(ug->u.a[xi->hid])); + + vk = (uint32_t)str->a[k]; + hid_a = str_idx->occ.a + str_idx->idx.a[vk>>1]; + hid_n = str_idx->idx.a[(vk>>1)+1] - str_idx->idx.a[vk>>1]; + for (z = 0; z < hid_n; z++) { + // if(qid == 95 || qid == 36) { + // fprintf(stderr,"[M::%s::qid->%u::tid->%lu] k->%lu, tid_occ->%u\n", + // __func__, qid, (hid_a[z]>>32), k, str_idx->str.a[hid_a[z]>>32].cn); + // } + if((hid_a[z]>>32) == qid) continue; + if(str_idx->str.a[hid_a[z]>>32].cn < 2) continue; + vz = (uint32_t)(str_idx->str.a[hid_a[z]>>32].a[(uint32_t)hid_a[z]]); + assert((vk>>1) == (vz>>1)); + kv_pushp(integer_aln_t, buf->b, &p); + p->vq = vk; p->tk = (uint32_t)hid_a[z]; + if((vk^vz)&1) p->tk = str_idx->str.a[hid_a[z]>>32].cn - p->tk - 1;///rev + p->tn_rev_qk = (hid_a[z]>>32); p->tn_rev_qk <<= 1; p->tn_rev_qk |= ((vk^vz)&1); + p->tn_rev_qk <<= 32; p->tn_rev_qk += k; + ///set score of this pair + p->sc = sck; + xi = &(uidx->idx->a[hid_a[z]>>32].bb.a[(str_idx->str.a[hid_a[z]>>32].a[(uint32_t)hid_a[z]])>>32]); + assert(((xi->hid<<1)+xi->rev)==vz); m = ug_occ_w(xi->ts, xi->te, &(ug->u.a[xi->hid])); + if(p->sc > m) p->sc = m; + } + } + + radix_sort_integer_aln_t_srt(buf->b.a, buf->b.a + buf->b.n); + b_n = buf->b.n; buf->sc.n = 0; + for (k = 1, z = 0; k <= b_n; k++) { + if(k == b_n || (buf->b.a[z].tn_rev_qk>>32) != (buf->b.a[k].tn_rev_qk>>32)) { + // if(qid == 95 || qid == 36) { + // fprintf(stderr,"[M::%s::qid->%u::tid->%lu]\n", __func__, qid, buf->b.a[z].tn_rev_qk>>33); + // } + if(integer_chain(qid, buf->b.a + z, k - z, z, buf, ug, str_idx, uidx->idx, &sc) && sc.v != (uint32_t)-1) { + if((buf->sc.n > 0) && ((buf->sc.a[buf->sc.n-1].v>>1) == (sc.v>>1))) { + if(buf->sc.a[buf->sc.n-1].sc < sc.sc) { + buf->sc.a[buf->sc.n-1] = sc; + } + } else { + kv_push(ul_chain_t, buf->sc, sc); + } + } + z = k; + } + } + + for (k = o->n = 0; k < buf->sc.n; k++) { + if(update_exact_ul_ovlps(uidx->idx, ug, str_idx->str.a, buf->b.a, &(buf->sc.a[k]), qid, buf, &res)) { + kv_push(ul2ul_t, *o, res); + } + } + // if(str->cn > 0) { + // for (k = 0, buf->b.n = 0; k < str->cn; k++) { + // xi = &(uidx->idx->a[qid].bb.a[str->a[k]>>32]); + + // // if(xi->ts) + // } + // // uint32_t v, nv; asg_arc_t *av; asg_t *g = ug->g; ul2ul_t *r; + // // v = ((uint32_t)str->a[0])^1; + // // nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + // // for (k = 0; k < nv; k++) { + // // if(av[k].del) continue; + // // kv_pushp(ul2ul_t, *o, &r); + // // r->hid = (av[k].v>>1) + ug_offset; + // // r->is_rev = (av[k].v^v)&1; + // // r->qs = 0; r->qe = av[k].ol; + + // // } + // // v = ((uint32_t)str->a[str->cn-1]); + // // nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + // } + o->cn = o->n; + if(o->n <= 0) return; +} + +#define ulg_id(ul2, i) (((i)>=(ul2).uln)?((i)-(ul2).uln):(i)) +#define ulg_type(ul2, i) (((i)>=(ul2).uln)?(0):(1)) +#define ulg_len(uidx, i) (ulg_type((uidx).uovl,(i))?((uidx).idx->a[ulg_id((uidx).uovl,(i))].rlen):((uidx).l1_ug->u.a[ulg_id((uidx).uovl,(i))].len)) +#define ulg_occ(uidx, i) (ulg_type((uidx).uovl,(i))?((uidx).pstr.str.a[ulg_id((uidx).uovl,(i))].cn):(1)) + +void integer_normalize_ovlp(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul2ul_idx_t *ul2) +{ + ul2ul_t *z; uint64_t k; ul2ul_item_t *t; + for (k = 0; k < o->n; k++) { + assert(ul2->item_idx[o->a[k].hid] != (uint32_t)-1); + z = get_ul_o(ul2, ulg_id(*ul2, o->a[k].hid), ulg_type(*ul2, o->a[k].hid), + ulg_id(*ul2, qid), ulg_type(*ul2, qid)); + // assert(z && z->hid == qid); + // if(o->a[k].is_rev == z->is_rev && (*z).qs == o->a[k].ts && (*z).qe == o->a[k].te && + // (*z).ts == o->a[k].qs && (*z).te == o->a[k].qe) { + // fprintf(stderr, "good::[M::%s] fn::%u(%c), fqs::%u, fqe::%u, fts::%u, fte::%u, rn::%u(%c), rqs::%u, rqe::%u, rts::%u, rte::%u\n", __func__, + // o->a[k].hid, "+-"[o->a[k].is_rev], o->a[k].qs, o->a[k].qe, o->a[k].ts, o->a[k].te, + // z->hid, "+-"[z->is_rev], z->qs, z->qe, z->ts, z->te); + // } else { + // fprintf(stderr, "bad::[M::%s] fn::%u(%c), fqs::%u, fqe::%u, fts::%u, fte::%u, rn::%u(%c), rqs::%u, rqe::%u, rts::%u, rte::%u\n", __func__, + // o->a[k].hid, "+-"[o->a[k].is_rev], o->a[k].qs, o->a[k].qe, o->a[k].ts, o->a[k].te, + // z->hid, "+-"[z->is_rev], z->qs, z->qe, z->ts, z->te); + // } + // continue; + if(z && o->a[k].hid < qid) continue; + if(!z) { + t = get_ul_ovlp(ul2, ulg_id(*ul2, o->a[k].hid), ulg_type(*ul2, o->a[k].hid)); + assert(t); + kv_pushp(ul2ul_t, *t, &z); + (*z).hid = qid; + (*z).qs = o->a[k].ts; (*z).qe = o->a[k].te; + (*z).ts = o->a[k].qs; (*z).te = o->a[k].qe; + (*z).qs_k = o->a[k].ts_k; (*z).qe_k = o->a[k].te_k; + (*z).ts_k = o->a[k].qs_k; (*z).te_k = o->a[k].qe_k; + (*z).is_rev = o->a[k].is_rev; (*z).is_del = o->a[k].is_del; + (*z).is_ct = o->a[k].is_ct; + if(!((*z).is_del)) { + t->cn++; + } + } else { + if((o->a[k].qe - o->a[k].qs) >= (z->te - z->ts)) { + (*z).qs = o->a[k].ts; (*z).qe = o->a[k].te; + (*z).ts = o->a[k].qs; (*z).te = o->a[k].qe; + (*z).qs_k = o->a[k].ts_k; (*z).qe_k = o->a[k].te_k; + (*z).ts_k = o->a[k].qs_k; (*z).te_k = o->a[k].qe_k; + (*z).is_rev = o->a[k].is_rev; (*z).is_del = o->a[k].is_del; + (*z).is_ct = o->a[k].is_ct; + } else { + o->a[k].qs = (*z).ts; o->a[k].qe = (*z).te; + o->a[k].ts = (*z).qs; o->a[k].te = (*z).qe; + o->a[k].qs_k = (*z).ts_k; o->a[k].qe_k = (*z).te_k; + o->a[k].ts_k = (*z).qs_k; o->a[k].te_k = (*z).qe_k; + o->a[k].is_rev = (*z).is_rev; o->a[k].is_del = (*z).is_del; + o->a[k].is_ct = (*z).is_ct; + } + } + } +} + + +void integer_normalize_ovlp_purge(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul2ul_idx_t *ul2) +{ + if(o->is_del) return; + ul2ul_t *z; uint64_t k, is_del; + for (k = is_del = 0; k < o->n; k++) { + if(o->a[k].is_del) continue; + z = get_ul_o(ul2, o->a[k].hid, 1, qid, 1); + if((!z) || (z->is_del)) { + o->a[k].is_del = 1; is_del++; + } + } + + if(is_del) { + for (k = o->cn = 0; k < o->n; k++) { + if(o->a[k].is_del) o->a[k].hid |= ((uint32_t)(0x80000000)); + else o->cn++; + } + + radix_sort_ul2ul_srt(o->a, o->a + o->n); + + for (k = 0; k < o->n; k++) { + if(o->a[k].hid&((uint32_t)(0x80000000))) { + o->a[k].hid -= ((uint32_t)(0x80000000)); + } + } + } + + if(o->cn == 0) o->is_del = 1; +} + +static inline int64_t integer_hit2arc_idx_contain(const ul2ul_t *z, int64_t qn, int64_t tn, uint64_t *dir) +{ + int64_t tn5, tn3, extn5, extn3, qsn = z->qs_k, qen = qn - ((int64_t)z->qe_k); + if (z->is_rev) tn5 = tn - z->te_k, tn3 = z->ts_k; + else tn5 = z->ts_k, tn3 = tn - z->te_k; + (*dir) = (uint64_t)-1; + + extn5 = ((qsn 0 || extn3 > 0) return MA_HT_INT;///overhang + if (qsn < tn5 && qen < tn3) { // query contained in target + return MA_HT_QCONT; + } else if (qsn > tn5 && qen > tn3) { // target contained in query + return MA_HT_TCONT; + } else if(qsn == tn5 && qen == tn3) { + (*dir) = (uint64_t)-1; + } else if (qsn > tn5) { ///query-to-target overlap + (*dir) = 0; + } else if(qsn < tn5) { ///target-to-query overlaps + (*dir) = 1; + } else if(qen > tn3) { ///target-to-query overlaps + (*dir) = 1; + } else if(qen < tn3) { ///query-to-target overlap + (*dir) = 0; + } + return MA_HT_DOUBLE; +} + +static inline int integer_hit2arc(const ul2ul_t *z, int64_t ql, int64_t tl, int64_t qocc, int64_t tocc, uint64_t qid, uint64_t tid, +int64_t min_ovlp, asg_arc_t *p) +{ + int64_t tl5, tl3, ext5, ext3, qs = z->qs, rf; + uint64_t u, v, l, rr; // u: query end; v: target end; l: length from u to v + + ///if query and target are in different strand + if (z->is_rev) tl5 = tl - z->te, tl3 = z->ts; // tl5: 5'-end overhang (on the query strand); tl3: similar + else tl5 = z->ts, tl3 = tl - z->te; + + ///ext5 and ext3 is the hang on left side and right side, respectively + ext5 = qs < tl5? qs : tl5; + ext3 = (((ql - ((int64_t)z->qe)) < tl3)? (ql - ((int64_t)z->qe)) : tl3); + + if(ext5 > 0 || ext3 > 0) return MA_HT_INT;///overhang + if (qs <= tl5 && (ql - (int64_t)z->qe) <= tl3) { // query contained in target + return MA_HT_QCONT; + } else if (qs >= tl5 && (ql - (int64_t)z->qe) >= tl3) { // target contained in query + return MA_HT_TCONT; + } else if (qs > tl5) { ///u = 0 means query-to-target overlap, l is the length of node in string graph (not the overlap length) + u = 0, v = !!(z->is_rev), l = qs - tl5; + } else { ///u = 1 means target-to-query overlaps, l is the length of node in string graph (not the overlap length) + u = 1, v = !(z->is_rev), l = (ql - z->qe) - tl3; + } + if ((int64_t)z->qe - qs + ext5 + ext3 < min_ovlp || (int64_t)z->te - (int64_t)z->ts + ext5 + ext3 < min_ovlp) { + return MA_HT_SHORT_OVLP; // short overlap + } + rf = integer_hit2arc_idx_contain(z, qocc, tocc, &rr); + if(rf != MA_HT_DOUBLE) return rf; + if(rr != (uint64_t)-1 && rr != u) { + // fprintf(stderr, "[M::%s::] z->is_rev::%u, z->qs::%u, z->qe::%u, z->ts::%u, z->te::%u, z->qs_k::%u, z->qe_k::%u, z->ts_k::%u, z->te_k::%u, ql::%ld, tl::%ld, qocc::%ld, tocc::%ld\n", __func__, + // z->is_rev, z->qs, z->qe, z->ts, z->te, z->qs_k, z->qe_k, z->ts_k, z->te_k, ql, tl, qocc, tocc); + return MA_HT_INT;///overhang + } + ///u = 0 / 1 means query-to-target / target-to-query overlaps, + ///l is the length of node in string graph (not the overlap length between two reads) + u |= qid<<1, v |= tid<<1; + /** + p->ul: |____________31__________|__________1___________|______________32_____________| + qn direction of overlap length of this node (not overlap length) + (in the view of query) + p->v : |___________31___________|__________1___________| + tn reverse direction of overlap + (in the view of target) + p->ol: overlap length + **/ + if(p) { + p->ul = (uint64_t)u<<32 | l, p->v = v, p->ol = ql - l, p->del = 0; + ///l is the length of node in string graph (not the overlap length) + + p->strong = 1; p->el = 1; p->no_l_indel = 1; + } + return l; +} + + +void integer_append_ug_ovlp(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul2ul_idx_t *ul2) +{ + // if(o->is_del) return; + uint64_t k; uc_block_t *xi; ul2ul_t *z; ul2ul_item_t *t; int32_t r; + ma_ug_t *ug = uidx->l1_ug; ul_str_t *str = &(uidx->pstr.str.a[qid]); + for (k = 0; k < str->cn; k++) { + xi = &(uidx->idx->a[qid].bb.a[str->a[k]>>32]); + assert(((xi->hid<<1)+xi->rev)==((uint32_t)str->a[k])); + + ///ul side + kv_pushp(ul2ul_t, *o, &z); + z->hid = xi->hid + ul2->uln; z->is_rev = xi->rev; z->is_del = o->is_del; z->is_ct = 0; + z->qs = xi->qs; z->qe = xi->qe; z->ts = xi->ts; z->te = xi->te; + z->qs_k = k; z->qe_k = k + 1; z->ts_k = 0; z->te_k = 1; + + r = integer_hit2arc(z, uidx->idx->a[qid].rlen, ug->u.a[xi->hid].len, uidx->pstr.str.a[qid].cn, + 1, qid, z->hid, 0, NULL); + if(r == MA_HT_INT) { + o->n--; continue; + } + ///ug side + t = get_ul_ovlp(ul2, xi->hid, 0); + assert(t); + kv_pushp(ul2ul_t, *t, &z); + z->hid = qid; z->is_rev = xi->rev; z->is_del = o->is_del; z->is_ct = 0; + z->qs = xi->ts; z->qe = xi->te; z->ts = xi->qs; z->te = xi->qe; + z->qs_k = 0; z->qe_k = 1; z->ts_k = k; z->te_k = k + 1; + } +} + +void integer_node_del(ul2ul_idx_t *ul2, uint64_t id, uint64_t is_ct) +{ + ul2ul_item_t *o = get_ul_ovlp(ul2, ulg_id(*ul2, id), ulg_type(*ul2, id)); + if(o) { + uint64_t k; ul2ul_t *z; + for (k = 0; k < o->cn; k++) { + if(is_ct) o->a[k].is_ct = 1; + else o->a[k].is_del = 1; + z = get_ul_o(ul2, ulg_id(*ul2, o->a[k].hid), ulg_type(*ul2, o->a[k].hid), ulg_id(*ul2, id), ulg_type(*ul2, id)); + if(is_ct) z->is_ct = 1; + else z->is_del = 1; + } + o->is_del = 1; + } +} + +void integer_containment_purge(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *q, ul2ul_idx_t *ul2) +{ + if(q->is_del) return; + uint64_t k; int32_t r; ul2ul_item_t *t; ul2ul_t *z; assert(qid == q->id); + for (k = 0; k < q->cn; k++) { + if(q->a[k].is_del || q->a[k].is_ct) continue; + t = get_ul_ovlp(ul2, ulg_id(*ul2, q->a[k].hid), ulg_type(*ul2, q->a[k].hid)); + assert(t); + if(t->is_del) continue; + r = integer_hit2arc(&(q->a[k]), ulg_len(*uidx, qid), ulg_len(*uidx, q->a[k].hid), + ulg_occ(*uidx, qid), ulg_occ(*uidx, q->a[k].hid), qid, q->a[k].hid, 0, NULL); + // assert(r != MA_HT_INT); + if (r == MA_HT_QCONT) { + q->a[k].is_ct = 1; + z = get_ul_o(ul2, ulg_id(*ul2, q->a[k].hid), ulg_type(*ul2, q->a[k].hid), ulg_id(*ul2, qid), ulg_type(*ul2, qid)); + assert(z && (!z->is_del) && (!z->is_ct)); z->is_ct = 1; + integer_node_del(ul2, qid, 1); + } else if (r == MA_HT_TCONT) { + q->a[k].is_ct = 1; + z = get_ul_o(ul2, ulg_id(*ul2, q->a[k].hid), ulg_type(*ul2, q->a[k].hid), ulg_id(*ul2, qid), ulg_type(*ul2, qid)); + assert(z && (!z->is_del) && (!z->is_ct)); z->is_ct = 1; + integer_node_del(ul2, q->a[k].hid, 1); + } + } + + // if(ulg_type(*ul2, k)) { + // for (k = 0; k < q->cn; k++) { + // if((!(q->a[k].is_del)) && (!(q->a[k].is_ct))) break; + // } + // if(k >= q->cn) q->is_del = 1; + // } +} + static void worker_integer_correction(void *data, long i, int tid) // callback for kt_for() { @@ -4831,12 +5440,223 @@ static void worker_integer_postprecess(void *data, long i, int tid) // callback ul_resolve_t *uidx = (ul_resolve_t *)data; integer_ml_t *sl = &(uidx->str_b); integer_t *buf = &(sl->buf[tid]); + ul2ul_item_t *it = get_ul_ovlp(&(uidx->uovl), i, 1); + if(!it) return; + assert(uidx->pstr.str.a[i].cn > 1); // uc_block_t *uls; uint64_t uls_n; uint32_t v; // uint64_t *srt_a, srt_n, is_circle = ((uidx->psrt.idx.a[i]&((uint64_t)(0x100000000)))?1:0); // srt_a = uidx->psrt.srt.a + (uint32_t)uidx->psrt.idx.a[i]; srt_n = uidx->psrt.idx.a[i]>>33; // if(srt_n == 0) return; // integer_candidate(uidx, srt_a, srt_n, is_circle, buf); - integer_candidate(uidx, buf, i, (asm_opt.purge_level_primary == 0?1:0)); + integer_gen_ovlp(uidx, buf, i, it, uidx->uovl.uln); +} + + +void gen_integer_normalize(ul_resolve_t *uidx) +{ + uint64_t k; ul2ul_idx_t *u2o = &(uidx->uovl); ul2ul_item_t *it; + for (k = 0; k < u2o->uln; k++) { + it = get_ul_ovlp(u2o, k, 1); + if(!it) continue; + assert(uidx->pstr.str.a[k].cn > 1); + if(it->is_del) continue; + integer_normalize_ovlp(uidx, k, it, u2o); + } +} + + +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); + 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; + } + if(k < o->cn) { + is_srt = 1; + } else { + for (; k < o->n; k++) { + if(!(o->a[k].is_del)) break; + if(k > o->cn && o->a[k].hid < o->a[k-1].hid) break; + } + if(k < o->n) is_srt = 1; + } + + if(is_srt) { + for (k = o->cn = 0; k < o->n; k++) { + if(o->a[k].is_del) o->a[k].hid |= ((uint32_t)(0x80000000)); + else o->cn++; + } + + radix_sort_ul2ul_srt(o->a, o->a + o->n); + + for (k = 0; k < o->n; k++) { + if(o->a[k].hid&((uint32_t)(0x80000000))) { + o->a[k].hid -= ((uint32_t)(0x80000000)); + } + } + } + + for (k = 0, buf->u.n = 0; k < o->cn; k++) { + kv_push(uint64_t, buf->u, (o->a[k].qs<<1)); + kv_push(uint64_t, buf->u, (o->a[k].qe<<1)|1); + } + radix_sort_srt64(buf->u.a, buf->u.a + buf->u.n); + + int64_t dp, old_dp; uint64_t start, end, 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; + if(o->cn && o->n > o->cn) radix_sort_ul2ul_srt(o->a, o->a + o->n); + o->cn = 0; + } +} + +void clean_srt_integer(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul2ul_idx_t *ul2) +{ + uint64_t k, l, i, m; ul2ul_t *p; + radix_sort_ul2ul_srt(o->a, o->a + o->n); + for(k = 1, l = m = o->cn = 0; k <= o->n; k++) { + if(k == o->n || o->a[k].hid != o->a[l].hid) { + for (i = l, p = &(o->a[l]); i < k; i++) { + if(((p->is_del) && (!o->a[i].is_del)) || + ((o->a[i].qe - o->a[i].qs + o->a[i].te - o->a[i].ts) > (p->qe - p->qs + p->te - p->ts))) { + p = &(o->a[i]); + } + } + if(!(p->is_del)) o->cn++; + o->a[m++] = *p; + l = k; + } + } + o->n = m; + if(o->n == o->cn) return; + + for (k = o->cn = 0; k < o->n; k++) { + if(o->a[k].is_del) o->a[k].hid |= ((uint32_t)(0x80000000)); + else o->cn++; + } + radix_sort_ul2ul_srt(o->a, o->a + o->n); + + for (k = 0; k < o->n; k++) { + if(o->a[k].hid&((uint32_t)(0x80000000))) { + o->a[k].hid -= ((uint32_t)(0x80000000)); + } + } +} + + +static void worker_detect_chimeric(void *data, long i, int tid) // callback for kt_for() +{ + ul_resolve_t *uidx = (ul_resolve_t *)data; + integer_ml_t *sl = &(uidx->str_b); + integer_t *buf = &(sl->buf[tid]); + ul2ul_item_t *it = get_ul_ovlp(&(uidx->uovl), i, 1); + if(!it) return; + assert(uidx->pstr.str.a[i].cn > 1); + if(it->is_del) return; + clip_integer_chimeric(uidx, i, it, &(uidx->uovl), buf, 1); +} + + +void chimeric_integer_deal(ul_resolve_t *uidx) +{ + uint64_t k; ul2ul_idx_t *u2o = &(uidx->uovl); ul2ul_item_t *it; + kt_for(uidx->str_b.n_thread, worker_detect_chimeric, uidx, uidx->idx->n); + for (k = 0; k < u2o->uln; k++) { + it = get_ul_ovlp(u2o, k, 1); + if(!it) continue; + assert(uidx->pstr.str.a[k].cn > 1); + integer_normalize_ovlp_purge(uidx, k, it, u2o); + } +} + + + + +static void worker_integert_clean(void *data, long i, int tid) // callback for kt_for() +{ + ul_resolve_t *uidx = (ul_resolve_t *)data; + ul2ul_item_t *it = get_ul_ovlp(&(uidx->uovl), ulg_id(uidx->uovl, (uint32_t)i), + ulg_type(uidx->uovl, (uint32_t)i)); + if(!it) return; + if((uint32_t)i >= uidx->uovl.uln) it->id = i; + clean_srt_integer(uidx, i, it, &(uidx->uovl)); +} + +static void worker_integert_debug_sym(void *data, long i, int tid) // callback for kt_for() +{ + ul_resolve_t *uidx = (ul_resolve_t *)data; + ul2ul_item_t *o = get_ul_ovlp(&(uidx->uovl), ulg_id(uidx->uovl, (uint32_t)i), + ulg_type(uidx->uovl, (uint32_t)i)); + if(!o) return; + ul2ul_idx_t *ul2 = &(uidx->uovl); + // if(o->is_del) { + // assert(o->cn == 0); + // } + + ul2ul_t *z; uint64_t k, qid = i, ct_n = 0, del_n = 0; assert(o->id == qid); + for (k = 0; k < o->n; k++) { + if(o->a[k].is_del) del_n++; + else if(o->a[k].is_ct) ct_n++; + assert(ul2->item_idx[o->a[k].hid] != (uint32_t)-1); + z = get_ul_o(ul2, ulg_id(*ul2, o->a[k].hid), ulg_type(*ul2, o->a[k].hid), + ulg_id(*ul2, qid), ulg_type(*ul2, qid)); + assert(z && z->hid == qid && z->is_del == o->a[k].is_del && z->is_ct == o->a[k].is_ct); + + // if(!(o->a[k].is_rev == z->is_rev && (*z).qs == o->a[k].ts && (*z).qe == o->a[k].te && + // (*z).ts == o->a[k].qs && (*z).te == o->a[k].qe)){ + // fprintf(stderr, "uln::%u[M::%s] fn::%u(%c), fqs::%u, fqe::%u, fts::%u, fte::%u, fdel::%u, rn::%u(%c), rqs::%u, rqe::%u, rts::%u, rte::%u, rdel::%u\n", ul2->uln, __func__, + // o->a[k].hid, "+-"[o->a[k].is_rev], o->a[k].qs, o->a[k].qe, o->a[k].ts, o->a[k].te, o->a[k].is_del, + // z->hid, "+-"[z->is_rev], z->qs, z->qe, z->ts, z->te, z->is_del); + // } + // if(o->a[k].is_rev == z->is_rev && (*z).qs == o->a[k].ts && (*z).qe == o->a[k].te && + // (*z).ts == o->a[k].qs && (*z).te == o->a[k].qe) { + // fprintf(stderr, "good::[M::%s] fn::%u(%c), fqs::%u, fqe::%u, fts::%u, fte::%u, rn::%u(%c), rqs::%u, rqe::%u, rts::%u, rte::%u\n", __func__, + // o->a[k].hid, "+-"[o->a[k].is_rev], o->a[k].qs, o->a[k].qe, o->a[k].ts, o->a[k].te, + // z->hid, "+-"[z->is_rev], z->qs, z->qe, z->ts, z->te); + // } else { + // fprintf(stderr, "bad::[M::%s] fn::%u(%c), fqs::%u, fqe::%u, fts::%u, fte::%u, rn::%u(%c), rqs::%u, rqe::%u, rts::%u, rte::%u\n", __func__, + // o->a[k].hid, "+-"[o->a[k].is_rev], o->a[k].qs, o->a[k].qe, o->a[k].ts, o->a[k].te, + // z->hid, "+-"[z->is_rev], z->qs, z->qe, z->ts, z->te); + // } + assert(o->a[k].is_rev == z->is_rev && (*z).qs == o->a[k].ts && (*z).qe == o->a[k].te && + (*z).ts == o->a[k].qs && (*z).te == o->a[k].qe); + if(kcn) assert(!o->a[k].is_del); + else assert(o->a[k].is_del); + } + + if(o->is_del) assert((ct_n + del_n) == o->n); } void print_primary_ul_chain(ul_str_t *p_str, ul_vec_t *ul, ma_ug_t *ug) @@ -5222,8 +6042,224 @@ void rebuid_idx(ul_resolve_t *uidx) init_ul_str_idx_t(uidx); } +void shrink_1b(ma_ug_t *ug, uc_block_t *z, uint32_t is_forward) +{ + if(z->ts != 0 || z->te != ug->g->seq[z->hid].len) return; + uc_block_t bc = *z; + if(z->ts + 1 < z->te && z->qs + 1 < z->qe) { + uint32_t off = get_offset_adjust(1, z->te-z->ts, z->qe-z->qs); + if(is_forward) { + if((!z->rev)) { + z->ts += 1; z->qs += off; + } else { + z->te -= 1; z->qs += off; + } + } else { + if((!z->rev)) { + z->te -= 1; z->qe -= off; + } else { + z->ts += 1; z->qe -= off; + } + } + if((ugl_cover_check(bc.ts, bc.te, &(ug->u.a[bc.hid]))) && + (!ugl_cover_check(z->ts, z->te, &(ug->u.a[z->hid])))) { + *z = bc; + } + } + +} + +void renew_ul_vec_t(ul_vec_t *x, ma_ug_t *ug) +{ + if(x->bb.n <= 0) return; + int64_t l, k, z, bn = x->bb.n, dd; int64_t qs = x->bb.a[0].qs, qe; + for (l = 0, k = 1, qs = 0; k <= bn; k++) { + if(k == bn || x->bb.a[k].pidx == (uint32_t)-1) { + if(k - l > 0) { + for (z = k - 1, dd = 0; z >= l; z--) { + if(z > l) { + assert(x->bb.a[z].pidx == z-1); + dd += x->bb.a[z].pdis; + } else { + assert(x->bb.a[z].pidx == (uint32_t)-1); + dd += ug->g->seq[x->bb.a[z].hid].len; + } + dd -= (int64_t)(ug->g->seq[x->bb.a[z].hid].len - (x->bb.a[z].te - x->bb.a[z].ts)); + } + if(dd < 0) dd = 0; + qe = qs + dd; + if(k < bn) { + dd = qe + (int64_t)(x->bb.a[k].qs) - (int64_t)(x->bb.a[k-1].qe); + if(dd < 0) dd = 0; + } + for (z = k - 1; z >= l; z--) { + qs = qe - (int64_t)(x->bb.a[z].te - x->bb.a[z].ts); + if(qs < 0) qs = 0; + x->bb.a[z].qs = qs; x->bb.a[z].qe = qe; + if(z > l) { + qe += (int64_t)(ug->g->seq[x->bb.a[z].hid].len - (x->bb.a[z].te - x->bb.a[z].ts)); + qe -= (int64_t)(x->bb.a[z].pdis); + if(qe < 0) qe = 0; + } + } + + if(k < bn) qs = dd; + } + l = k; + } + } + x->rlen = x->bb.a[bn-1].qe; +} + +void shrink_ul0(all_ul_t *uls, ul_str_t *str, uint64_t id, integer_t *buf, ma_ug_t *ug) +{ + uint32_t k, c_k, p_k, cv, pv, bl, i; uc_block_t *xi; buf->u.n = 0; nid_t *np = NULL; + asg_arc_t *av; uint32_t nv, s, e, m, d, mm, is_conn; uint64_t *z; ul_vec_t *x; + if(str->cn < 2) return; + for (k = 0, bl = 0, c_k = p_k = pv = (uint32_t)-1; k < str->cn; k++) { + xi = &(uls->a[id].bb.a[str->a[k]>>32]); is_conn = 0; + c_k = k; cv = ((uint32_t)str->a[k])^1; + assert(((xi->hid<<1)+xi->rev)==((uint32_t)str->a[k])); + if(p_k != (uint32_t)-1) { + if(xi->pidx == (str->a[p_k]>>32)) { + av = asg_arc_a(ug->g, cv); nv = asg_arc_n(ug->g, cv); + for (i = 0; i < nv; i++) { + if(av[i].del) continue; + if(av[i].v == pv) break; + } + if(i < nv) { + s = (uint32_t)av[i].ul; + if(s == xi->pdis) { + is_conn = 1; + } else { + d = (s>=xi->pdis?s-xi->pdis:xi->pdis-s); + mm = MAX(s, xi->pdis); + if((d <= (mm*0.08)) || (d <= 512)) is_conn = 1; + } + } + } else if(xi->pidx == (uint32_t)-1) { + is_conn = 1; + } + } + if(is_conn) { + bl++; + } else { + if(bl > 0) { + kv_pushp(uint64_t, buf->u, &z); + (*z) = k - 1 - bl; (*z) <<= 32; (*z) += bl + 1; + } + bl = 0; + } + p_k = c_k; pv = cv; + } + + if(bl > 0) { + kv_pushp(uint64_t, buf->u, &z); + (*z) = k - 1 - bl; (*z) <<= 32; (*z) += bl + 1; + } + // fprintf(stderr, "[M::%s::] buf->u.n::%u, str->cn::%u\n", __func__, buf->u.n, str->cn); + if(buf->u.n <= 0) { + str->cn = str->n = 0; + } else if(buf->u.n > 0) { + uint32_t pidx, aidx, pdis; + for (i = 0; i + 1 < buf->u.n; i++) { + s = (buf->u.a[i]>>32); e = s + ((uint32_t)buf->u.a[i]); + ///new id + kv_pushp(nid_t, uls->nid, &np); + np->n = uls->nid.a[id].n; MALLOC(np->a, np->n+1); + memcpy(np->a, uls->nid.a[id].a, np->n); np->a[np->n] = '\0'; + ///new ovlps + kv_pushp(ul_vec_t, *uls, &x); memset(x, 0, sizeof(*x)); + MALLOC(x->bb.a, e - s); x->bb.n = x->bb.m = e - s; + + for (k = s, m = 0; k < e; k++, m++) { + pidx = uls->a[id].bb.a[str->a[k]>>32].pidx; + aidx = uls->a[id].bb.a[str->a[k]>>32].aidx; + pdis = uls->a[id].bb.a[str->a[k]>>32].pdis; + x->bb.a[m] = uls->a[id].bb.a[str->a[k]>>32]; + if(k == s) { + x->bb.a[m].pidx = x->bb.a[m].pdis = (uint32_t)-1; + } else if(pidx != (uint32_t)-1) { + x->bb.a[m].pidx = m - 1; x->bb.a[m].pdis = pdis; + } + if(k + 1 == e) { + x->bb.a[m].aidx = (uint32_t)-1; + } else if(aidx != (uint32_t)-1) { + x->bb.a[m].aidx = m + 1; + } + } + x->bb.n = m; + assert(x->bb.n > 1); + + + + shrink_1b(ug, &(x->bb.a[0]), 1); + shrink_1b(ug, &(x->bb.a[x->bb.n-1]), 0); + d = x->bb.a[0].qs; + for (k = 0; k < x->bb.n; k++) { + x->bb.a[k].qs -= d; x->bb.a[k].qe -= d; + } + x->rlen = x->bb.a[x->bb.n-1].qe; + renew_ul_vec_t(x, ug); + } + + ///the last one; update in-place + s = (buf->u.a[i]>>32); e = s + ((uint32_t)buf->u.a[i]); x = &(uls->a[id]); + for (k = s, m = 0; k < e; k++, m++) { + pidx = x->bb.a[str->a[k]>>32].pidx; + aidx = x->bb.a[str->a[k]>>32].aidx; + pdis = x->bb.a[str->a[k]>>32].pdis; + x->bb.a[m] = x->bb.a[str->a[k]>>32]; + if(k == s) { + x->bb.a[m].pidx = x->bb.a[m].pdis = (uint32_t)-1; + } else if(pidx != (uint32_t)-1) { + x->bb.a[m].pidx = m - 1; x->bb.a[m].pdis = pdis; + } + if(k + 1 == e) { + x->bb.a[m].aidx = (uint32_t)-1; + } else if(aidx != (uint32_t)-1) { + x->bb.a[m].aidx = m + 1; + } + } + x->bb.n = m; + assert(x->bb.n > 1); + shrink_1b(ug, &(x->bb.a[0]), 1); + shrink_1b(ug, &(x->bb.a[x->bb.n-1]), 0); + d = x->bb.a[0].qs; + for (k = 0; k < x->bb.n; k++) { + x->bb.a[k].qs -= d; x->bb.a[k].qe -= d; + } + x->rlen = x->bb.a[x->bb.n-1].qe; + renew_ul_vec_t(x, ug); + } +} + +void shrink_uls(ul_resolve_t *uidx) +{ + all_ul_t *uls = uidx->idx; uint64_t k, ul_n = uls->n; + for (k = 0; k < ul_n; k++) { + // fprintf(stderr, "+[M::%s::k->%lu] ul_n::%lu, uls->n::%u\n", __func__, k, ul_n, uls->n); + shrink_ul0(uls, &(uidx->pstr.str.a[k]), k, &(uidx->str_b.buf[0]), uidx->l1_ug); + // fprintf(stderr, "-[M::%s::k->%lu] ul_n::%lu, uls->n::%u\n", __func__, k, ul_n, uls->n); + } + rebuid_idx(uidx); +} + +void print_ul_seq(ul_resolve_t *uidx, uint64_t id) +{ + uint64_t i; ul_str_t *str = &(uidx->pstr.str.a[id]); ma_ug_t *ug = uidx->init_ug; + fprintf(stderr, "%.*s\tid::%lu\t", (int32_t)uidx->idx->nid.a[id].n, + uidx->idx->nid.a[id].a, id); + for (i = 0; i < str->cn; i++) { + fprintf(stderr, "utg%.6d%c(%c)\t", (((uint32_t)str->a[i])>>1)+1, + "lc"[ug->u.a[(((uint32_t)str->a[i])>>1)].circ], "+-"[(((uint32_t)str->a[i])&1)]); + } + fprintf(stderr,"\n"); +} + void ul_re_correct(ul_resolve_t *uidx, uint64_t n_r) { + // print_ul_seq(uidx, 2062); print_ul_seq(uidx, 2093); uint64_t k, occ, n_circle; ///uidx->str_b.n_thread = 1; for (k = 0; k < n_r; k++) { kt_for(uidx->str_b.n_thread, worker_integer_correction, uidx, uidx->idx->n); @@ -5235,9 +6271,2496 @@ void ul_re_correct(ul_resolve_t *uidx, uint64_t n_r) fprintf(stderr, "-[M::%s::round->%lu] # corrected UL reads::%lu, # circle UL reads::%lu\n", __func__, k, occ, n_circle); rebuid_idx(uidx); + // print_ul_seq(uidx, 2062); print_ul_seq(uidx, 2093); + // exit(1); } - // kt_for(uidx->str_b.n_thread, worker_integer_postprecess, uidx, uidx->idx->n); + shrink_uls(uidx); +} + +void append_utg_es(ul_resolve_t *uidx) +{ + uint64_t k; ul2ul_idx_t *u2o = &(uidx->uovl); ul2ul_item_t *it; + for (k = 0; k < u2o->uln; k++) { + it = get_ul_ovlp(u2o, k, 1); + if(!it) continue; + assert(uidx->pstr.str.a[k].cn > 1); + integer_append_ug_ovlp(uidx, k, it, u2o); + } + + kt_for(uidx->str_b.n_thread, worker_integert_clean, uidx, u2o->tot);///all ul + ug + + for (k = 0; k < u2o->tot; k++) { + it = get_ul_ovlp(u2o, ulg_id(*u2o, k), ulg_type(*u2o, k)); + if(!it) continue; + integer_normalize_ovlp(uidx, k, it, u2o); + } + + kt_for(uidx->str_b.n_thread, worker_integert_clean, uidx, u2o->tot);///all ul + ug +} + + +void remove_integert_containment(ul_resolve_t *uidx) +{ + uint64_t k; ul2ul_idx_t *u2o = &(uidx->uovl); ul2ul_item_t *o; + for (k = 0; k < u2o->tot; k++) { + o = get_ul_ovlp(u2o, ulg_id(*u2o, k), ulg_type(*u2o, k)); + if(!o) continue; + integer_containment_purge(uidx, k, o, u2o); + } + + kt_for(uidx->str_b.n_thread, worker_integert_clean, uidx, u2o->tot);///all ul + ug +} + +void print_integert_ovlp_stat(ul2ul_idx_t *ul2) +{ + uint64_t i, k, occ_r = 0, occ_o = 0; ul2ul_item_t *o; + for (i = 0; i < ul2->uln; i++) { + o = get_ul_ovlp(ul2, ulg_id(*ul2, i), ulg_type(*ul2, i)); + if((!o) || (o->is_del)) continue; + occ_r++; + for (k = 0; k < o->n; k++) { + if(o->a[k].is_del || o->a[k].is_ct || (!ulg_type(*ul2, o->a[k].hid))) continue; + occ_o++; + } + } + fprintf(stderr, "[M::%s::] # UL reads::%lu, # UL ovlps::%lu\n", __func__, occ_r, occ_o); +} + +asg_t *integer_sg_gen(ul_resolve_t *uidx, uint64_t min_ovlp) +{ + ul2ul_idx_t *ul2 = &(uidx->uovl); ma_ug_t *ug = uidx->l1_ug; asg_t *raw_g = ug->g; + uint64_t i, k, is_del, v, w, nv, z; ul2ul_item_t *o, *ow; int32_t r; asg_arc_t t, *p; asg_arc_t *av; + asg_t *g = asg_init(); + for (i = 0; i < ul2->tot; i++) { + is_del = 0; + o = get_ul_ovlp(ul2, ulg_id(*ul2, i), ulg_type(*ul2, i)); + if((!o) || (o->is_del)) is_del = 1; + asg_seq_set(g, i, ulg_len(*uidx, i), is_del); + g->seq[i].c = 0; + } + // CALLOC(g->seq_vis, g->n_seq*2); + + for (i = 0; i < ul2->tot; ++i) { + o = get_ul_ovlp(ul2, ulg_id(*ul2, i), ulg_type(*ul2, i)); + if((!o) || (o->is_del)) continue; + for (k = 0; k < o->n; k++) { + if(o->a[k].is_del || o->a[k].is_ct) continue; + r = integer_hit2arc(&(o->a[k]), ulg_len(*uidx, i), ulg_len(*uidx, o->a[k].hid), + ulg_occ(*uidx, i), ulg_occ(*uidx, o->a[k].hid), i, o->a[k].hid, min_ovlp, &t); + if (r >= 0) { + p = asg_arc_pushp(g); + *p = t; + } + } + if(i >= ul2->uln) {///is a node of ug + v = (i-ul2->uln)<<1; + nv = asg_arc_n(raw_g, v); av = asg_arc_a(raw_g, v); + for (z = 0; z < nv; z++) { + if(av[z].del) continue; + w = av[z].v + (ul2->uln<<1); + ow = get_ul_ovlp(ul2, ulg_id(*ul2, (w>>1)), ulg_type(*ul2, (w>>1))); + if((!ow) || (ow->is_del)) continue; + p = asg_arc_pushp(g); + *p = av[z]; p->ul += (((uint64_t)ul2->uln)<<33); p->v += (ul2->uln<<1); + } + + v = ((i-ul2->uln)<<1)+1; + nv = asg_arc_n(raw_g, v); av = asg_arc_a(raw_g, v); + for (z = 0; z < nv; z++) { + if(av[z].del) continue; + w = av[z].v + (ul2->uln<<1); + ow = get_ul_ovlp(ul2, ulg_id(*ul2, (w>>1)), ulg_type(*ul2, (w>>1))); + if((!ow) || (ow->is_del)) continue; + p = asg_arc_pushp(g); + *p = av[z]; p->ul += (((uint64_t)ul2->uln)<<33); p->v += (ul2->uln<<1); + } + } + } + asg_cleanup(g); + g->r_seq = g->n_seq; + // fprintf(stderr, "[M::%s::] # ig nodes::%u, # ig archs::%u\n", __func__, g->n_seq, g->n_arc); + return g; +} + +inline uint64_t get_ul_occ(ul_resolve_t *uidx, uint64_t id) +{ + return uidx->uovl.cc.uc[id]; +} + +inline void get_iug_u_raw_occ(ul_resolve_t *uidx, uint32_t id, uint32_t *ul_occ, uint32_t *raw_ug_occ) +{ + if(ul_occ) *ul_occ = uidx->uovl.cc.raw_uc[id]; + if(raw_ug_occ) *raw_ug_occ = uidx->uovl.i_ug->u.a[id].n - uidx->uovl.cc.raw_uc[id]; +} + +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]; + for (i = 0; i < ug->u.n; ++i) { // the Segment lines in GFA + p = &ug->u.a[i]; + if(p->m == 0) continue; + sprintf(name, "%s%.6d%c", prefix, i + 1, "lc"[p->circ]); + fprintf(fp, "S\t%s\t*\tLN:i:%d\trd:i:%lu\n", name, p->len, get_ul_occ(uidx, i)); + + for (j = l = 0; j < p->n; j++) { + if(p->a[j] != (uint64_t)-1) { + x = p->a[j]>>33; + if(ulg_type(uidx->uovl, x)) {//read + x = ulg_id(uidx->uovl, x); + fprintf(fp, "A\t%s\t%d\t%c\t%.*s\t%d\t%d\tid:i:%d\tHG:A:*\n", name, l, "+-"[p->a[j]>>32&1], + (int32_t)uidx->idx->nid.a[x].n, uidx->idx->nid.a[x].a, 0, uidx->idx->a[x].rlen, x); + } else { ///node + x = ulg_id(uidx->uovl, x); s = &(uidx->init_ug->u.a[x]); + fprintf(fp, "A\t%s\t%d\t%c\tutg%.6d%c\t%d\t%d\tid:i:%d\tHG:A:*\n", name, l, "+-"[p->a[j]>>32&1], + x + 1, "lc"[s->circ], 0, s->len, x); + } + } + else + { + fprintf(fp, "A\t%s\t%d\t*\t*\t*\t*\tid:i:*\tHG:A:*\n", name, l); + } + l += (uint32_t)p->a[j]; + } + } + + if(ug->g) + { + asg_arc_t* au = NULL; + uint32_t nu, u, v; + for (i = 0; i < ug->u.n; ++i) { + if(ug->u.a[i].m == 0) continue; + if(ug->u.a[i].circ) + { + fprintf(fp, "L\t%s%.6dc\t+\t%s%.6dc\t+\t%dM\tL1:i:%d\n", + prefix, i+1, prefix, i+1, 0, ug->u.a[i].len); + fprintf(fp, "L\t%s%.6dc\t-\t%s%.6dc\t-\t%dM\tL1:i:%d\n", + prefix, i+1, prefix, i+1, 0, ug->u.a[i].len); + } + u = i<<1; + au = asg_arc_a(ug->g, u); + nu = asg_arc_n(ug->g, u); + for (j = 0; j < nu; j++) + { + if(au[j].del) continue; + v = au[j].v; + fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\tL2:i:%u\n", + prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], + prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), au[j].ou); + } + + + u = (i<<1) + 1; + au = asg_arc_a(ug->g, u); + nu = asg_arc_n(ug->g, u); + for (j = 0; j < nu; j++) + { + if(au[j].del) continue; + v = au[j].v; + fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\tL2:i:%u\n", + prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], + prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), au[j].ou); + } + } + } +} + + +void output_integer_graph(ul_resolve_t *uidx, ma_ug_t *iug, const char *nn) +{ + char* gfa_name = NULL; MALLOC(gfa_name, strlen(nn)+50); + sprintf(gfa_name, "%s.integer.noseq.gfa", nn); + FILE* fp = fopen(gfa_name, "w"); free(gfa_name); + if (!fp) return; + ma_integer_ug_print0(iug, uidx, 0, "itg", fp); + fclose(fp); +} + + +void print_raw_uls_seq(ul_resolve_t *uidx, const char *nn) +{ + char* gfa_name = NULL; MALLOC(gfa_name, strlen(nn)+70); + sprintf(gfa_name, "%s.raw.integer.seq.log", nn); + FILE* fp = fopen(gfa_name, "w"); free(gfa_name); + if (!fp) return; + ma_ug_t *ug = uidx->init_ug; all_ul_t *aln = uidx->idx; + uint64_t id; uc_block_t *a = NULL; int64_t k, a_n; + for (id = 0; id < aln->n; id++) { + a = aln->a[id].bb.a; a_n = aln->a[id].bb.n; + if(a_n == 0) continue; + fprintf(fp,"%.*s\tid::%lu\t", (int32_t)aln->nid.a[id].n, aln->nid.a[id].a, id); + for (k = 0; k < a_n && ug_occ_w(a[k].ts, a[k].te, &(ug->u.a[a[k].hid])) == 0; k++); + for (; k < a_n; k++) { + if(ug_occ_w(a[k].ts, a[k].te, &(ug->u.a[a[k].hid])) == 0) break; + fprintf(fp, "utg%.6d%c(%c)\t", a[k].hid + 1, "lc"[ug->u.a[a[k].hid].circ], "+-"[a[k].rev]); + } + fprintf(fp,"\n"); + } + fclose(fp); +} + + +void print_uls_seq(ul_resolve_t *uidx, const char *nn) +{ + char* gfa_name = NULL; MALLOC(gfa_name, strlen(nn)+50); + sprintf(gfa_name, "%s.integer.seq.log", nn); + FILE* fp = fopen(gfa_name, "w"); free(gfa_name); + if (!fp) return; + uint64_t k, i; all_ul_t *uls = uidx->idx; ul2ul_item_t *o; + ul_str_idx_t *str_idx = &(uidx->pstr); ma_ug_t *ug = uidx->init_ug; + for (k = 0; k < uls->n; k++) { + if(str_idx->str.a[k].cn < 2) continue; + o = get_ul_ovlp(&(uidx->uovl), ulg_id(uidx->uovl, k), ulg_type(uidx->uovl, k)); + if(!o) continue; + fprintf(fp,"%.*s\tid::%lu\tdel::%u\t", (int32_t)uidx->idx->nid.a[k].n, uidx->idx->nid.a[k].a, k, o->is_del); + for (i = 0; i < str_idx->str.a[k].cn; i++) { + fprintf(fp, "utg%.6d%c(%c)\t", (((uint32_t)str_idx->str.a[k].a[i])>>1)+1, + "lc"[ug->u.a[(((uint32_t)str_idx->str.a[k].a[i])>>1)].circ], "+-"[(((uint32_t)str_idx->str.a[k].a[i])&1)]); + } + fprintf(fp,"\n"); + } + fclose(fp); +} + + +void print_uls_ovs(ul_resolve_t *uidx, const char *nn) +{ + char* gfa_name = NULL; MALLOC(gfa_name, strlen(nn)+50); + sprintf(gfa_name, "%s.integer.ovlp.log", nn); + FILE* fp = fopen(gfa_name, "w"); free(gfa_name); + if (!fp) return; + uint64_t i, k, qv, tv; ma_ug_t *ug = uidx->l1_ug; ul2ul_idx_t *ul2 = &(uidx->uovl); ul2ul_item_t *o; + for (i = 0; i < ul2->uln; ++i) { + o = get_ul_ovlp(ul2, ulg_id(*ul2, i), ulg_type(*ul2, i)); + if(!o) continue; + for (k = 0; k < o->n; k++) { + qv = i; tv = o->a[k].hid; + if(ulg_type(uidx->uovl, qv)) {//ul read + qv = ulg_id(uidx->uovl, qv); + fprintf(fp, "%.*s(id::%lu)\t%u\t%u(i::%u)\t%u(i::%u)\t%c\t", (int32_t)uidx->idx->nid.a[qv].n, uidx->idx->nid.a[qv].a, + i, uidx->idx->a[qv].rlen, o->a[k].qs, o->a[k].qs_k, o->a[k].qe, o->a[k].qe_k, "+-"[o->a[k].is_rev]); + } else {///ug node + qv = ulg_id(uidx->uovl, qv); + fprintf(fp, "utg%.6d%c(id::%lu)\t%u\t%u(i::%u)\t%u(i::%u)\t%c\t", (int32_t)qv + 1, "lc"[ug->u.a[qv].circ], + i, ug->u.a[qv].len, o->a[k].qs, o->a[k].qs_k, o->a[k].qe, o->a[k].qe_k, "+-"[o->a[k].is_rev]); + } + + if(ulg_type(uidx->uovl, tv)) {//ul read + tv = ulg_id(uidx->uovl, tv); + fprintf(fp, "%.*s(id::%u)\t%u\t%u(i::%u)\t%u(i::%u)\t", (int32_t)uidx->idx->nid.a[tv].n, uidx->idx->nid.a[tv].a, + o->a[k].hid, uidx->idx->a[tv].rlen, o->a[k].ts, o->a[k].ts_k, o->a[k].te, o->a[k].te_k); + + } else { + tv = ulg_id(uidx->uovl, tv); + fprintf(fp, "utg%.6d%c(id::%u)\t%u\t%u(i::%u)\t%u(i::%u)\t\t", (int32_t)tv + 1, "lc"[ug->u.a[tv].circ], + o->a[k].hid, ug->u.a[tv].len, o->a[k].ts, o->a[k].ts_k, o->a[k].te, o->a[k].te_k); + } + + fprintf(fp, "del::%u\tct::%u\n", o->a[k].is_del, o->a[k].is_ct); + } + } + fclose(fp); +} + + + +uint64_t gen_srt_cov_interval(uint64_t *b, uint64_t b_n) +{ + uint64_t i, m, start, end; int64_t dp, old_dp; + radix_sort_srt64(b, b + b_n); + for (i = m = 0, dp = 0, start = 0; i < b_n; ++i) { + old_dp = dp; + if (b[i]&1) --dp; + else ++dp; + + if (old_dp < 1 && dp >= 1) {///old_dp < dp, b.a[j] is qs + start = b[i]>>1; + } else if (old_dp >= 1 && dp < 1){ + end = b[i]>>1; + b[m] = start; b[m] <<= 32; b[m] += end; m++; + } + } + + return m; +} + +inline uint64_t get_remove_hifi_occ_back(ul_resolve_t *uidx, uint64_t thres, uint64_t *v_a, uint64_t v_n, asg64_v *buf) +{ + ul2ul_idx_t *idx = &(uidx->uovl); + ma_ug_t *iug = idx->i_ug, *raw = uidx->l1_ug; + uint64_t k, l, z, m, p, *raw_a, raw_n, raw_id, iug_id, iug_off, del_n, keep_ns[2], dup_raw, n_mask, pn, fn, full_del, is, ie; + uinfo_srt_warp_t *iu; uint64_t *s_a, s_n, ms, me; uinfo_srt_t *ps; + for (k = del_n = fn = 0, buf->n = dup_raw = 0; k < v_n; k++) { + iu = &(idx->cc.iug_a[v_a[k]]); ///integer unitigs + for (z = 0; z < iu->n; z++) { + raw_id = (iu->a[z].v>>1); + raw_a = idx->cc.iug_b + idx->cc.iug_idx[raw_id]; + raw_n = idx->cc.iug_idx[raw_id+1] - idx->cc.iug_idx[raw_id]; + assert(raw_n > 0); + for (m = keep_ns[0] = keep_ns[1] = 0, pn = buf->n, n_mask = 0; m < raw_n; m++) { + iug_id = raw_a[m]>>32; iug_off = (uint32_t)raw_a[m]; + if(iug->g->seq[iug_id].del) continue; + if(raw_a[m]&((uint64_t)(0x8000000000000000))) { + n_mask++; + continue; + } + if(iug_id == v_a[k] && iug_off == z) {///query itself + raw_a[m] |= ((uint64_t)(0x8000000000000000)); + p = (uint32_t)-1; p <<= 32; p += idx->cc.iug_idx[raw_id] + m; + kv_push(uint64_t, *buf, p); + } else { + if(!(idx->cc.iug_a[iug_id].a[iug_off].v&1)) { + keep_ns[0] = MAX(idx->cc.iug_a[iug_id].a[iug_off].n, keep_ns[0]); + } else { + keep_ns[1] = MAX(idx->cc.iug_a[iug_id].a[iug_off].n, keep_ns[1]); + } + } + } + assert(buf->n == pn + 1); + if(keep_ns[0] + keep_ns[1] >= raw->u.a[raw_id].n) { + keep_ns[0] = keep_ns[1] = raw->u.a[raw_id].n; + } + + is = keep_ns[0]; ie = raw->u.a[raw_id].n - keep_ns[1];///[is, ie) -> uncovered coordinates + if(ie > is) { + if(!(iu->a[z].v&1)) { + ms = 0; me = iu->a[z].n; + } else { + ms = raw->u.a[raw_id].n - iu->a[z].n; me = raw->u.a[raw_id].n; + } + if(MAX(is, ms) < MIN(ie, me)) { + del_n += MIN(ie, me) - MAX(is, ms); + buf->a[pn] <<= 32; buf->a[pn] >>= 32; buf->a[pn] |= (raw_id<<32); fn++; + if(n_mask) dup_raw = 1; + } + } + } + } + + if(del_n > thres && dup_raw) { + radix_sort_srt64(buf->a, buf->a + buf->n); del_n = 0; + for (l = 0, k = 1; k <= fn; k++) { + if(k == fn || (buf->a[k]>>32) != (buf->a[l]>>32)) { + raw_id = (buf->a[l]>>32); pn = buf->n; + raw_a = idx->cc.iug_b + idx->cc.iug_idx[raw_id]; + raw_n = idx->cc.iug_idx[raw_id+1] - idx->cc.iug_idx[raw_id]; + assert(raw_n > 0); + for (m = full_del = 0; m < raw_n; m++) { + iug_id = raw_a[m]>>32; iug_off = (uint32_t)raw_a[m]; + if(iug->g->seq[iug_id].del) continue; + if(raw_a[m]&((uint64_t)(0x8000000000000000))) { + kv_push(uint64_t, *buf, (idx->cc.iug_a[iug_id].a[iug_off].s<<1)); + kv_push(uint64_t, *buf, (idx->cc.iug_a[iug_id].a[iug_off].e<<1)+1); + if(idx->cc.iug_a[iug_id].a[iug_off].s == 0 && + idx->cc.iug_a[iug_id].a[iug_off].e == raw->u.a[raw_id].len) { + full_del = 1; break; + } + } + } + + s_a = buf->a + buf->n; s_n = buf->n - pn; + assert(s_n > 0); + + if(s_n > 2 && full_del == 0) { + s_n = gen_srt_cov_interval(s_a, s_n); + } else if(full_del) {///an interval has already cover the whole unitig + s_a[0] = raw->u.a[raw_id].len; s_n = 1; + } else {///sn == 2; + s_a[0] >>= 1; s_a[0] <<= 32; s_a[0] += (s_a[1]>>1); s_n = 1; + } + + + for (m = full_del = 0; m < raw_n && full_del < s_n; m++) { + iug_id = raw_a[m]>>32; iug_off = (uint32_t)raw_a[m]; + if(iug->g->seq[iug_id].del) continue; + if(raw_a[m]&((uint64_t)(0x8000000000000000))) continue; + ps = &(idx->cc.iug_a[iug_id].a[iug_off]); + for (z = 0; z < s_n; z++) { + if(s_a[z] == ((uint64_t)-1)) continue; + is = s_a[z]>>32; ie = (uint32_t)s_a[z]; + ms = MAX(is, ps->s); me = MIN(ie, ps->e); + if(ms >= me) continue; + ///[ms, me) is unlikely to be contained in the [is, ie) + assert(ms == is || me == ie); + if(ms == is) is = me; + else if(me == ie) ie = ms; + if(is >= ie) { + s_a[z] = ((uint64_t)-1); full_del++; + } else { + s_a[z] = is; s_a[z] <<= 32; s_a[z] += ie; + } + } + } + + for (z = 0; z < s_n; z++) { + if(s_a[z] == ((uint64_t)-1)) continue; + is = s_a[z]>>32; ie = (uint32_t)s_a[z]; + del_n += ug_occ_w(is, ie, &(raw->u.a[raw_id])); + } + if(del_n > thres) break; + l = k; buf->n = pn; + } + } + } + + for (k = 0; k < buf->n; k++) { + assert(idx->cc.iug_b[(uint32_t)buf->a[k]]&((uint64_t)(0x8000000000000000))); + idx->cc.iug_b[(uint32_t)buf->a[k]] <<= 1; idx->cc.iug_b[(uint32_t)buf->a[k]] >>= 1; + } + + + return ((del_n > thres)?0:1); +} + + +inline uint64_t get_remove_hifi_occ(ul_resolve_t *uidx, uint64_t thres, uint64_t *v_a, uint64_t v_n, asg64_v *buf, uint32_t *del_occ) +{ + ul2ul_idx_t *idx = &(uidx->uovl); + ma_ug_t *iug = idx->i_ug, *raw = uidx->l1_ug; + uint64_t k, l, z, m, p, *raw_a, raw_n, raw_id, iug_id, iug_off, del_n, keep_ns[2], del_ns[2], dup_raw, n_mask, pn, fn, is, ie; + uinfo_srt_warp_t *iu; uint64_t ms, me; + for (k = del_n = fn = 0, buf->n = dup_raw = 0; k < v_n; k++) { + iu = &(idx->cc.iug_a[v_a[k]>>1]); ///integer unitigs + for (z = 0; z < iu->n; z++) { + raw_id = (iu->a[z].v>>1); + raw_a = idx->cc.iug_b + idx->cc.iug_idx[raw_id]; + raw_n = idx->cc.iug_idx[raw_id+1] - idx->cc.iug_idx[raw_id]; + assert(raw_n > 0); + for (m = keep_ns[0] = keep_ns[1] = 0, pn = buf->n, n_mask = 0; m < raw_n; m++) { + iug_id = (raw_a[m]<<1)>>33; iug_off = (uint32_t)raw_a[m]; + // if(iug_id >= iug->g->n_seq){ + // fprintf(stderr, ">>>[M::%s::] v_a[%lu]>>1::%lu, raw_id::%lu, raw_n::%lu, m::%lu, raw_a[m]::%lu, uidx->uovl.cc.iug_b[455]::%lu, iug_id::%lu, iug_off::%lu, iug->g->n_seq::%u\n", __func__, + // k, v_a[k]>>1, raw_id, raw_n, m, raw_a[m], uidx->uovl.cc.iug_b[455], iug_id, iug_off, (uint32_t)iug->g->n_seq); + // } + if(iug->g->seq[iug_id].del) continue; + if(raw_a[m]&((uint64_t)(0x8000000000000000))) { + n_mask++; + continue; + } + if(iug_id == (v_a[k]>>1) && iug_off == z) {///query itself + raw_a[m] |= ((uint64_t)(0x8000000000000000)); + p = (uint32_t)-1; p <<= 32; p += idx->cc.iug_idx[raw_id] + m; + kv_push(uint64_t, *buf, p); + } else { + if(!(idx->cc.iug_a[iug_id].a[iug_off].v&1)) { + keep_ns[0] = MAX(idx->cc.iug_a[iug_id].a[iug_off].n, keep_ns[0]); + } else { + keep_ns[1] = MAX(idx->cc.iug_a[iug_id].a[iug_off].n, keep_ns[1]); + } + } + } + assert(buf->n == pn + 1); + if(keep_ns[0] + keep_ns[1] >= raw->u.a[raw_id].n) { + keep_ns[0] = keep_ns[1] = raw->u.a[raw_id].n; + } + + is = keep_ns[0]; ie = raw->u.a[raw_id].n - keep_ns[1];///[is, ie) -> uncovered coordinates + if(ie > is) { + if(!(iu->a[z].v&1)) { + ms = 0; me = iu->a[z].n; + } else { + ms = raw->u.a[raw_id].n - iu->a[z].n; me = raw->u.a[raw_id].n; + } + if(MAX(is, ms) < MIN(ie, me)) { + del_n += MIN(ie, me) - MAX(is, ms); + buf->a[pn] <<= 32; buf->a[pn] >>= 32; buf->a[pn] |= (raw_id<<32); fn++; + if(n_mask) dup_raw = 1; + } + } + } + } + + if(del_n > thres && dup_raw) { + radix_sort_srt64(buf->a, buf->a + buf->n); del_n = 0; + for (l = 0, k = 1; k <= fn; k++) { + if(k == fn || (buf->a[k]>>32) != (buf->a[l]>>32)) { + raw_id = (buf->a[l]>>32); + raw_a = idx->cc.iug_b + idx->cc.iug_idx[raw_id]; + raw_n = idx->cc.iug_idx[raw_id+1] - idx->cc.iug_idx[raw_id]; + assert(raw_n > 0); + for (m = del_ns[0] = del_ns[1] = keep_ns[0] = keep_ns[1] = 0; m < raw_n; m++) { + iug_id = (raw_a[m]<<1)>>33; iug_off = (uint32_t)raw_a[m]; + if(iug->g->seq[iug_id].del) continue; + if(raw_a[m]&((uint64_t)(0x8000000000000000))) { + if(!(idx->cc.iug_a[iug_id].a[iug_off].v&1)) { + del_ns[0] = MAX(idx->cc.iug_a[iug_id].a[iug_off].n, del_ns[0]); + } else { + del_ns[1] = MAX(idx->cc.iug_a[iug_id].a[iug_off].n, del_ns[1]); + } + } else { + if(!(idx->cc.iug_a[iug_id].a[iug_off].v&1)) { + keep_ns[0] = MAX(idx->cc.iug_a[iug_id].a[iug_off].n, keep_ns[0]); + } else { + keep_ns[1] = MAX(idx->cc.iug_a[iug_id].a[iug_off].n, keep_ns[1]); + } + } + } + assert(del_ns[0] + del_ns[1] > 0); + if(keep_ns[0] + keep_ns[1] >= raw->u.a[raw_id].n) { + keep_ns[0] = keep_ns[1] = raw->u.a[raw_id].n; + } + + + is = keep_ns[0]; ie = raw->u.a[raw_id].n - keep_ns[1];///[is, ie) -> uncovered coordinates + if(ie > is) { + if(del_ns[0] + del_ns[1] >= raw->u.a[raw_id].n) { + del_n += ie - is; + } else { + ms = 0; me = del_ns[0]; + if(ms < me && MAX(is, ms) < MIN(ie, me)) { + del_n += MIN(ie, me) - MAX(is, ms); + } + + ms = raw->u.a[raw_id].n - del_ns[1]; me = raw->u.a[raw_id].n; + if(ms < me && MAX(is, ms) < MIN(ie, me)) { + del_n += MIN(ie, me) - MAX(is, ms); + } + } + } + if(del_n > thres) break; + l = k; + } + } + } + + for (k = 0; k < buf->n; k++) { + assert(idx->cc.iug_b[(uint32_t)buf->a[k]]&((uint64_t)(0x8000000000000000))); + idx->cc.iug_b[(uint32_t)buf->a[k]] <<= 1; idx->cc.iug_b[(uint32_t)buf->a[k]] >>= 1; + } + // fprintf(stderr, "*[M::%s::] del_n::%lu\n", __func__, del_n); + if(del_occ) *del_occ = del_n; + return ((del_n > thres)?0:1); +} +#define bg_correct (1) +#define bg_wrong (2) +#define bg_ambiguous (3) +#define bg_unavailable ((uint32_t)-1) +inline uint32_t get_bg_flag(ul_resolve_t *uidx, uint32_t v, uint32_t w) +{ + asg_t *g = uidx->uovl.bg.bg; asg_arc_t *av; uint32_t k, nv; + if(g->seq[v>>1].del || g->seq[w>>1].del) return bg_unavailable; + + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + for (k = 0; k < nv; k++) { + if(av[k].v == w) break; + } + if(k >= nv) return bg_unavailable; + if(av[k].ou == 3) return bg_correct; + if(av[k].ou == 2) return bg_wrong; + return bg_ambiguous; +} + +inline void get_bridges(ul_resolve_t *uidx, uint64_t *v_a, uint64_t v_n, uint32_t *w_occ, uint32_t *am_occ) +{ + ul_bg_t *bg = &(uidx->uovl.bg); uint64_t k, w, am, l; uint32_t uv, uw, bv, bw; + for (k = w = am = 0, uv = uw = (uint32_t)-1; k < v_n; k++) { + w += bg->w_n[v_a[k]>>1]; am += bg->a_n[v_a[k]>>1]; uw = v_a[k]; + if(uv != (uint32_t)-1) { + bv = uidx->uovl.cc.iug_a[uv>>1].a[((uv&1)?(0):(uidx->uovl.cc.iug_a[uv>>1].n-1))].v; if(uv&1) bv ^= 1; + bw = uidx->uovl.cc.iug_a[uw>>1].a[((uw&1)?(uidx->uovl.cc.iug_a[uw>>1].n-1):(0))].v; if(uw&1) bw ^= 1; + if((ulg_type(uidx->uovl, (bv>>1))) && (ulg_type(uidx->uovl, (bw>>1)))) { + l = get_bg_flag(uidx, bv, bw); + if(l == bg_wrong) w++; + if(l == bg_ambiguous) am++; + } + } + uv = uw; + } + if(w_occ) *w_occ = w; + if(am_occ) *am_occ = am; +} + +static inline void ulg_seq_del(ma_ug_t *ug, uint32_t s) +{ + uint32_t k; asg_t *g = ug->g; + g->seq[s].del = 1; + for (k = 0; k < 2; ++k) { + uint32_t i, v = s<<1 | k; + uint32_t nv = asg_arc_n(g, v); + asg_arc_t *av = asg_arc_a(g, v); + for (i = 0; i < nv; ++i) { + av[i].del = 1; + asg_arc_del(g, av[i].v^1, v^1, 1); + } + } + free(ug->u.a[s].a); memset(&(ug->u.a[s]), 0, sizeof(ug->u.a[s])); +} + + +uint32_t ulg_arc_cut_tips(ul_resolve_t *uidx, ma_ug_t *ug, uint32_t max_ext, uint32_t max_ext_hifi, 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 n_vtx = g->n_seq<<1, v, w, i, k, cnt = 0, nv, kv, pb, w_occ, a_occ, ul_occ, del_occ; + asg_arc_t *av = NULL; uint64_t lw; + b = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); + for (v = 0, b->n = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + + av = asg_arc_a(g, v^1); nv = asg_arc_n(g, v^1); + for (i = kv = 0; i < nv; i++) { + if (av[i].del) continue; + kv++; break; + } + + if(kv) continue; + // kv = ug->u.a[v>>1].n;/// get_ul_occ(uidx, v>>1); + get_iug_u_raw_occ(uidx, v>>1, &ul_occ, NULL); kv = ul_occ; + for (i = 0, w = v; i < max_ext; i++) { + if(asg_end(g, w^1, &lw, NULL)!=0) break; + w = (uint32_t)lw; + // kv += ug->u.a[w>>1].n; ///get_ul_occ(uidx, w>>1); + get_iug_u_raw_occ(uidx, w>>1, &ul_occ, NULL); kv += ul_occ; + } + + if(kv <= max_ext) kv_push(uint64_t, *b, (((uint64_t)kv)<<32)|v); + } + + radix_sort_srt64(b->a, b->a + b->n); + + for (k = 0; k < b->n; k++) { + v = (uint32_t)(b->a[k]); + if (g->seq[v>>1].del) continue; + + av = asg_arc_a(g, v^1); nv = asg_arc_n(g, v^1); + for (i = kv = 0; i < nv; i++) { + if (av[i].del) continue; + kv++; break; + } + + if(kv) continue; + + pb = b->n; kv_push(uint64_t, *b, v); + // kv = ug->u.a[v>>1].n;/// get_ul_occ(uidx, v>>1); + get_iug_u_raw_occ(uidx, v>>1, &ul_occ, NULL); kv = ul_occ; + for (i = 0, w = v; i < max_ext; i++) { + if(asg_end(g, w^1, &lw, NULL)!=0) break; + w = (uint32_t)lw; kv_push(uint64_t, *b, w); + // kv += ug->u.a[w>>1].n; ///get_ul_occ(uidx, w>>1); + get_iug_u_raw_occ(uidx, w>>1, &ul_occ, NULL); kv += ul_occ; + } + + + + if(kv <= max_ext && get_remove_hifi_occ(uidx, max_ext_hifi, b->a + pb, b->n - pb, ub, &del_occ)) { + // fprintf(stderr, "*[M::%s::] k::%u, v>>1::%u, v&1::%u, kv::%u, max_ext::%u, max_ext_hifi::%u\n\n", + // __func__, k, v>>1, v&1, kv, max_ext, max_ext_hifi); + get_bridges(uidx, b->a + pb, b->n - pb, &w_occ, &a_occ); + if(w_occ > 0 || a_occ > 0 || kv == 0 || del_occ == 0) { + for (i = pb; i < b->n; i++) ulg_seq_del(ug, (b->a[i]>>1)); + cnt += b->n - pb; + } + } + b->n = pb; + } + + // stats_sysm(g); + if(!in) free(tx.a); if(!ib) free(tb.a); + if (cnt > 0) asg_cleanup(g); + + return cnt; +} + +uint32_t is_het_ulg_edge(ul_resolve_t *uidx, uint32_t uv, uint32_t uw) +{ + ma_ug_t *iug = uidx->uovl.i_ug; bubble_type *bub = uidx->bub; uint32_t qid, tid; + qid = iug->u.a[uv>>1].a[((uv&1)?(0):(iug->u.a[uv>>1].n-1))]>>33; + tid = iug->u.a[uw>>1].a[((uw&1)?(iug->u.a[uw>>1].n-1):(0))]>>33; + + if(!ulg_type(uidx->uovl, qid)) return (!IF_HOM(ulg_id(uidx->uovl, qid), *bub)); + if(!ulg_type(uidx->uovl, tid)) return (!IF_HOM(ulg_id(uidx->uovl, tid), *bub)); + + if(uidx->uovl.item_idx[qid] == (uint32_t)-1) return (uint32_t)-1; + ul2ul_item_t *o = &(uidx->uovl.a[uidx->uovl.item_idx[qid]]); + ul2ul_t *z; uint64_t k, n_hom, n_het; ul_str_t *str; + for (k = 0, z = NULL; k < o->cn; k++) { + if(o->a[k].hid != tid) continue; + z = &(o->a[k]); + break; + } + if(!z) return (uint32_t)-1; + str = &(uidx->pstr.str.a[ulg_id(uidx->uovl, qid)]); + + for (k = z->qs_k, n_hom = n_het = 0; k < z->qe_k; k++) { + if(IF_HOM((((uint32_t)str->a[k])>>1), *bub)) { + n_hom++; + } else { + n_het++; break; + } + } + + if(n_het) return 1; + if(n_hom) return 0; + return (uint32_t)-1; +} + + +int32_t usg_topocut_aux(ul_resolve_t *uidx, ma_ug_t *ug, uint32_t v, int32_t max_ext, int32_t max_ext_hifi, asg64_v *b, asg64_v *ub) +{ + int32_t n_ext; asg_arc_t *av; uint32_t w = v, nv, i, kv, ul_occ, pn = b->n; + for (n_ext = 0; n_ext < max_ext; v = w) { + av = asg_arc_a(ug->g, v^1); nv = asg_arc_n(ug->g, v^1); + for (i = kv = 0; i < nv && kv <= 1; i++) { + if (av[i].del) continue; + kv++; + } + if(kv!=1) break; + get_iug_u_raw_occ(uidx, v>>1, &ul_occ, NULL); n_ext += ul_occ; kv_push(uint64_t, *b, v); + + av = asg_arc_a(ug->g, v); nv = asg_arc_n(ug->g, v); + for (i = kv = 0; i < nv && kv <= 1; i++) { + if (av[i].del) continue; + kv++; w = av[i].v; + } + if(kv!=1) break; + } + + if(n_ext < max_ext) { + if(get_remove_hifi_occ(uidx, max_ext_hifi, b->a + pn, b->n - pn, ub, NULL)) { + b->n = pn; + return 1; + } + } + + b->n = pn; + return 0; +} + +uint32_t ulg_arc_cut_length(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 hom_check, 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, i, k, kv, kw, nv, nw, cnt = 0, n_vtx = g->n_seq<<1, n_het, n_hom, ff, ol_max, mm_ol, to_del; + asg_arc_t *av, *aw, *ve, *we, *vl_max, *wl_max; + + b = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); + for (v = b->n = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + if(g->seq_vis[v] == 0) { + 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; + if(hom_check) { + for (i = n_het = n_hom = 0; i < nv; ++i) { + if(av[i].del) continue; + ff = is_het_ulg_edge(uidx, av[i].ul>>32, av[i].v); assert(ff != (uint32_t)-1); + if(ff) n_het++; + else n_hom++; + if(n_het > 0 && n_hom > 0) break; + } + if(n_het == 0 || n_hom == 0) continue; + } + + for (i = 0; i < nv; ++i) { + if(av[i].del) continue; + kv_push(uint64_t, *b, (((uint64_t)av[i].ol)<<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); nw = asg_arc_n(g, w); + av = asg_arc_a(g, v); aw = asg_arc_a(g, w); + if(nv<=1 && nw <= 1) continue; + if(hom_check) { + ff = is_het_ulg_edge(uidx, g->arc[(uint32_t)b->a[k]].ul>>32, g->arc[(uint32_t)b->a[k]].v); + assert(ff != (uint32_t)-1); + if(ff) continue; + } + + ve = &(g->arc[(uint32_t)b->a[k]]); + for (i = 0; i < nw; ++i) { + if (aw[i].v == (v^1)) { + we = &(aw[i]); + break; + } + } + mm_ol = MIN(ve->ol, we->ol); ///ve and we are hom edges + + for (i = kv = ol_max = 0, vl_max = NULL; i < nv; ++i) { + if(av[i].del) continue; + kv++; + if(hom_check) { + ff = is_het_ulg_edge(uidx, av[i].ul>>32, av[i].v); assert(ff != (uint32_t)-1); + if(!ff) continue;///vl_max must be a het edge + } + if(ol_max < av[i].ol) ol_max = av[i].ol, vl_max = &(av[i]); + } + if (kv < 1 || (!vl_max)) continue; + if (kv >= 2) { + if (mm_ol > ol_max*len_rat) continue; + } + + for (i = kw = ol_max = 0, wl_max = NULL; i < nw; ++i) { + if(aw[i].del) continue; + kw++; + if(hom_check) { + ff = is_het_ulg_edge(uidx, aw[i].ul>>32, aw[i].v); assert(ff != (uint32_t)-1); + if(!ff) continue;///wl_max must be a het edge + } + if(ol_max < aw[i].ol) ol_max = aw[i].ol, wl_max = &(aw[i]); + } + if (kw < 1 || (!wl_max)) continue; + if (kw >= 2) { + if (mm_ol > ol_max*len_rat) continue; + } + + if (kv <= 1 && kw <= 1) continue; + + to_del = 0; + if(topo_level == 0) { + to_del = 1; + } else if(topo_level == 2) { + if (kv > 1 && kw > 1) to_del = 1; + } else { + if (kv > 1 && kw > 1) { + to_del = 1; + } else if (kw == 1) { + if (usg_topocut_aux(uidx, ug, w^1, max_ext, max_ext_hifi, b, ub)) to_del = 1; + } else if (kv == 1) { + if (usg_topocut_aux(uidx, ug, v^1, max_ext, max_ext_hifi, b, ub)) to_del = 1; + } + } + + if (to_del) { + ve->del = we->del = 1, ++cnt; + } + } + + + if(!in) free(tx.a); if(!ib) free(tb.a); + if (cnt > 0) asg_cleanup(g); + + return cnt; +} + +uint32_t is_ul_edge(ul_resolve_t *uidx, uint32_t uv, uint32_t uw) +{ + ma_ug_t *iug = uidx->uovl.i_ug; uint32_t qid, tid; + qid = iug->u.a[uv>>1].a[((uv&1)?(0):(iug->u.a[uv>>1].n-1))]>>33; + tid = iug->u.a[uw>>1].a[((uw&1)?(iug->u.a[uw>>1].n-1):(0))]>>33; + if((qid != tid) && (!ulg_type(uidx->uovl, qid)) && (!ulg_type(uidx->uovl, tid))) return 0; + return 1; +} + +uint32_t ulg_arc_cut_occ(ul_resolve_t *uidx, ma_ug_t *ug, int32_t max_ext, uint32_t max_ext_hifi, +uint32_t is_trio, uint32_t topo_level, asg64_v *in, asg64_v *ib) +{ + asg_t *g = ug->g; asg64_v tx = {0,0,0}, tb = {0,0,0}, *b = NULL, *ub = NULL; + uint32_t v, w, i, z, kv, kw, nv, nw, cnt = 0, n_vtx = g->n_seq<<1, to_del, n_ul, n_ug; + asg_arc_t *av, *aw; 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 = n_ul = n_ug = kv = 0; i < nv; ++i) { + if(av[i].del) continue; + if(is_ul_edge(uidx, av[i].ul>>32, av[i].v)) n_ul++; + else n_ug++; + kv++; + } + if(n_ul == 0 || n_ug == 0 || kv < 2) continue; + + for (i = 0; i < nv; ++i) { + if(av[i].del) continue; + if(is_ul_edge(uidx, av[i].ul>>32, av[i].v)) continue; + w = av[i].v^1; if(g->seq[w>>1].del) continue; + kw = get_arcs(g, w, NULL, 0); + if (kv <= 1 && kw <= 1) continue; + + to_del = 0; + if(topo_level == 0) { + to_del = 1; + } else if(topo_level == 2) { + if (kv > 1 && kw > 1) to_del = 1; + } else { + if (kv > 1 && kw > 1) { + to_del = 1; + } else if (kw == 1) { + if (usg_topocut_aux(uidx, ug, w^1, max_ext, max_ext_hifi, b, ub)) to_del = 1; + } else if (kv == 1) { + if (usg_topocut_aux(uidx, ug, v^1, max_ext, max_ext_hifi, b, ub)) to_del = 1; + } + } + + if (to_del) { + ++cnt; av[i].del = 1; + aw = asg_arc_a(g, w); nw = asg_arc_n(g, w); + for (z = 0; z < nw; ++z) { + if (aw[z].v == (v^1)) { + aw[z].del = 1; + break; + } + } + assert(z < nw); + } + } + } + + if(!in) free(tx.a); if(!ib) free(tb.a); + if (cnt > 0) asg_cleanup(g); + + return cnt; +} + + +uint32_t get_ul_path_info(ul_resolve_t *uidx, ma_ug_t *ug, uint32_t s, uint32_t *e, uint32_t *occ, +uint32_t *ul_cnt, uint32_t *bridge_w, uint32_t *bridge_am, asg64_v *b) +{ + uint32_t v = s, w = (uint32_t)-1, kv, kw, uv = (uint32_t)-1, uw = (uint32_t)-1, bv, bw, l; + ul_bg_t *bg = &(uidx->uovl.bg); + if(occ) (*occ) = 0; if(ul_cnt) (*ul_cnt) = 0; if(bridge_w) (*bridge_w) = 0; if(bridge_am) (*bridge_am) = 0; + + while (1) { + if(occ) (*occ)++; + kv = get_arcs(ug->g, v, &w, 1); + if(e) (*e) = v; + + ///if(b) kv_push(uint32_t, b->b, v>>1); + if(b) kv_push(uint64_t, *b, v); + if(ul_cnt) (*ul_cnt) += get_ul_occ(uidx, v>>1); + if(bridge_w || bridge_am) { + uw = v; + if(bridge_w) (*bridge_w) += bg->w_n[v>>1]; + if(bridge_am) (*bridge_am) += bg->a_n[v>>1]; + if(uv != (uint32_t)-1) { + bv = uidx->uovl.cc.iug_a[uv>>1].a[((uv&1)?(0):(uidx->uovl.cc.iug_a[uv>>1].n-1))].v; if(uv&1) bv ^= 1; + bw = uidx->uovl.cc.iug_a[uw>>1].a[((uw&1)?(uidx->uovl.cc.iug_a[uw>>1].n-1):(0))].v; if(uw&1) bw ^= 1; + if((ulg_type(uidx->uovl, (bv>>1))) && (ulg_type(uidx->uovl, (bw>>1)))) { + l = get_bg_flag(uidx, bv, bw); + if(l == bg_wrong && bridge_w) (*bridge_w)++; + if(l == bg_ambiguous && bridge_am) (*bridge_am)++; + } + } + uv = uw; + } + + if(kv == 0) return END_TIPS; + if(kv == 2) return TWO_OUTPUT; + if(kv > 2) return MUL_OUTPUT; + w = ug->g->arc[w].v; + ///up to here, kv=1 + ///kw must >= 1 + kw = get_arcs(ug->g, w^1, NULL, 0); + v = w; + + if(kw == 2) return TWO_INPUT; + if(kw > 2) return MUL_INPUT; + if(v == s) return LOOP; + } + + return LONG_TIPS; +} + +#define ul_path_w(ul_cnt, bridge_w, bridge_am) (((bridge_w)+(bridge_w))==0?((uint32_t)-1):((ul_cnt)/((bridge_w)+(bridge_w)))) + + +uint32_t usg_bridge_topocut_aux(ul_resolve_t *uidx, ma_ug_t *ug, uint32_t v, uint32_t max_ext, uint32_t max_ext_hifi, uint32_t topo_level, uint32_t double_check, asg64_v *b, asg64_v *ub) +{ + uint32_t k, z, bn = b->n, w, bridge_w = 1, bridge_am = 1, raw_ul, ul, del_occ, is_del = 0, kk; + asg_arc_t *av; uint32_t nv; + get_ul_path_info(uidx, ug, v, &w, NULL, NULL, (double_check?(&bridge_w):(NULL)), + (double_check?(&bridge_am):(NULL)), b); + for (k = bn, raw_ul = 0; k < b->n; k++) { + get_iug_u_raw_occ(uidx, b->a[k]>>1, &ul, NULL); raw_ul += ul; + } + + if(raw_ul <= max_ext && get_remove_hifi_occ(uidx, max_ext_hifi, b->a + bn, b->n - bn, ub, &del_occ)) { + b->n = bn; + if(topo_level == 0) { + is_del = 1; + } else if(double_check == 0 || bridge_w > 0 || bridge_am > 0 || raw_ul == 0 || del_occ == 0) { + av = asg_arc_a(ug->g, v^1); nv = asg_arc_n(ug->g, v^1); + for (z = 0; z < nv; z++) { + if(av[z].del) continue; + kk = get_arcs(ug->g, av[z].v^1, NULL, 0); + if(topo_level == 2 && kk > 1) continue; + ///==1, might be a tip + if (topo_level == 1 && usg_topocut_aux(uidx, ug, av[z].v, max_ext, max_ext_hifi, b, ub)) continue; + break; + } + + if(z >= nv) { + av = asg_arc_a(ug->g, w); nv = asg_arc_n(ug->g, w); + for (z = 0; z < nv; z++) { + if(av[z].del) continue; + kk = get_arcs(ug->g, av[z].v^1, NULL, 0); + if(topo_level == 2 && kk > 1) continue; + ///==1, might be a tip + if (topo_level == 1 && usg_topocut_aux(uidx, ug, av[z].v, max_ext, max_ext_hifi, b, ub)) continue; + break; + } + + if(z >= nv) is_del = 1; + } + } + } + b->n = bn; + + return is_del; +} + +uint32_t ulg_arc_cut_bridge(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 *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, i, k, kv, nv, cnt = 0, n_vtx = g->n_seq<<1, ff, ol_max, ul_max, mm_ol; + asg_arc_t *av, *ve, *vl_max; uint32_t ul_cnt, bridge_w, bridge_am, pb; + // fprintf(stderr, "+++[M::%s::] idx->cc.iug_b[455]::%lu\n", __func__, uidx->uovl.cc.iug_b[455]); + b = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); + for (v = b->n = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + if(g->seq_vis[v] == 0) { + 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; + get_ul_path_info(uidx, ug, av[i].v, NULL, NULL, &ul_cnt, &bridge_w, &bridge_am, NULL); + ff = ul_path_w(ul_cnt, bridge_w, bridge_am); + // if((v>>1) == 5) { + // fprintf(stderr, "+[M::%s::] v>>1::%u, v&1::%u, w>>1::%u, w&1::%u, ul_cnt::%u, bridge_w::%u, bridge_am::%u, ff::%u\n", __func__, + // v>>1, v&1, av[i].v>>1, av[i].v&1, ul_cnt, bridge_w, bridge_am, ff); + // } + kv_push(uint64_t, *b, (((uint64_t)ff)<<32) | ((uint64_t)(av-g->arc+i))); + } + } + } + + // fprintf(stderr, "---[M::%s::] idx->cc.iug_b[455]::%lu\n", __func__, uidx->uovl.cc.iug_b[455]); + 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); + if(nv<=1 && asg_arc_n(g, w) <= 1) continue; + if(get_arcs(ug->g, w, NULL, 0) != 1) continue; + + ve = &(g->arc[(uint32_t)b->a[k]]); + get_ul_path_info(uidx, ug, ve->v, NULL, NULL, &ul_cnt, &bridge_w, &bridge_am, NULL); + ff = ul_path_w(ul_cnt, bridge_w, bridge_am); + if(ff == (uint32_t)-1) continue; + mm_ol = ff; + + for (i = kv = ol_max = ul_max = 0, vl_max = NULL; i < nv; ++i) { + if(av[i].del) continue; + kv++; + get_ul_path_info(uidx, ug, av[i].v, NULL, NULL, &ul_cnt, &bridge_w, &bridge_am, NULL); + ff = ul_path_w(ul_cnt, bridge_w, bridge_am); + if((ol_max < ff) || (ol_max == ff && ul_max < ul_cnt)) { + ol_max = ff; ul_max = ul_cnt; vl_max = &(av[i]); + } + } + + // if((v>>1) == 5) { + // fprintf(stderr, "***[M::%s::] v>>1::%u, v&1::%u, w>>1::%u, w&1::%u, ul_cnt::%u, bridge_w::%u, bridge_am::%u, ff::%u, mm_ol::%u, ol_max::%u, kv::%u\n", __func__, + // v>>1, v&1, w>>1, w&1, ul_cnt, bridge_w, bridge_am, ff, mm_ol, ol_max, kv); + // } + + if (kv <= 1 || (!vl_max)) continue; + if (kv >= 2) { + if (mm_ol > ol_max*len_rat) continue; + } + + // if((v>>1) == 5) { + // fprintf(stderr, "###[M::%s::] v>>1::%u, v&1::%u, w>>1::%u, w&1::%u, ul_cnt::%u, bridge_w::%u, bridge_am::%u, ff::%u, mm_ol::%u, ol_max::%u\n", __func__, + // v>>1, v&1, w>>1, w&1, ul_cnt, bridge_w, bridge_am, ff, mm_ol, ol_max); + // } + + + if(usg_bridge_topocut_aux(uidx, ug, ve->v, max_ext, max_ext_hifi, topo_level, 1, b, ub)) { + pb = b->n; + get_ul_path_info(uidx, ug, ve->v, NULL, NULL, NULL, NULL, NULL, b); + for (i = pb; i < b->n; i++) ulg_seq_del(ug, (b->a[i]>>1)); + cnt += b->n - pb; b->n = pb; + } + } + + + if(!in) free(tx.a); if(!ib) free(tb.a); + if (cnt > 0) asg_cleanup(g); + return cnt; +} + + +inline asg_arc_t* get_specfic_edge(asg_t *g, uint32_t v, uint32_t w) +{ + asg_arc_t *av = asg_arc_a(g, v); uint32_t nv = asg_arc_n(g, v), k; + for (k = 0; k < nv; k++) { + if(av[k].del) continue; + if(av[k].v == w) break; + } + + if(k < nv) return (&av[k]); + return NULL; +} + +ul2ul_t* get_ul_spec_ovlp(ul2ul_idx_t *z, uint64_t qid, uint64_t tid) +{ + // uint64_t x = (is_qul?(qid):(qid+z->uln)); + if(z->item_idx[qid] == (uint32_t)-1) return NULL; + ul2ul_item_t *o = &(z->a[z->item_idx[qid]]); uint64_t k; + for (k = 0; k < o->cn; k++) { + if(o->a[k].hid != tid) continue; + return &(o->a[k]); + } + return NULL; +} + +uint64_t gen_ug_integer_seq_on_fly(ul_resolve_t *uidx, uint64_t *u_a, uint64_t u_n, asg64_v *res) +{ + ul2ul_idx_t *idx = &(uidx->uovl); ma_ug_t *raw = uidx->l1_ug; + uint64_t k, rn = res->n, nn, u_rev, uv, uw, bv, bw, is_bv_ul, is_bw_ul, x; int64_t z, zn; + ma_ug_t *iug = idx->i_ug; uinfo_srt_warp_t *seq; ul2ul_t *m; + for (k = 0; k < u_n; k++) { + seq = &(idx->cc.iug_a[u_a[k]>>1]); u_rev = u_a[k]&1; + // if(u_n == 6 && u_a[0] == 37683 && u_a[5] == 45894) { + // fprintf(stderr, "\n[M::%s::] k::%lu, u_a[k]>>1::%lu, u_a[k]&1::%lu\n", __func__, k, u_a[k]>>1, u_a[k]&1); + // for (z = 0, zn = seq->n; z < zn; z++) { + // fprintf(stderr, "[M::%s::] v>>1::%u, v&1::%u\n", __func__, seq->a[z].v>>1, seq->a[z].v&1); + // } + // } + if(k > 0) { + uv = u_a[k-1]; uw = u_a[k]; + bv = iug->u.a[uv>>1].a[((uv&1)?(0):(iug->u.a[uv>>1].n-1))]>>32; if(uv&1) bv ^= 1; //pre + bw = iug->u.a[uw>>1].a[((uw&1)?(iug->u.a[uw>>1].n-1):(0))]>>32; if(uw&1) bw ^= 1; //current + is_bv_ul = ulg_type((*idx), (bv>>1)); is_bw_ul = ulg_type((*idx), (bw>>1)); + if(is_bw_ul) { + m = get_ul_spec_ovlp(idx, bv>>1, bw>>1); + assert(m && (!m->is_del)); + nn = res->n + seq->n - (m->te_k - m->ts_k); + kv_resize(uint64_t, *res, nn); + nn = m->te_k - m->ts_k; ///skipped length + zn = seq->n; + // fprintf(stderr, "[M::%s::] bv>>1::%lu, bv&1::%lu, bw>>1::%lu, bw&1::%lu, m->qs_k::%u, m->qe_k::%u, m->ts_k::%u, m->te_k::%u\n", __func__, + // bv>>1, bv&1, bw>>1, bw&1, m->qs_k, m->qe_k, m->ts_k, m->te_k); + if(!u_rev) { + ///debug + for (z = 0; z < ((int64_t)nn); z++) { + // if(res->a[res->n-nn+z] != seq->a[z].v) { + // fprintf(stderr, "****[M::%s::] k::%lu, z::%ld, nn::%lu, res->n::%u, res->v>>1::%lu, res->v&1::%lu, seq->a[z].v>>1::%u, seq->a[z].v&1::%u\n", + // __func__, k, z, nn, (uint32_t)res->n, res->a[res->n-nn+z]>>1, res->a[res->n-nn+z]&1, seq->a[z].v>>1, seq->a[z].v&1); + // } + assert(res->a[res->n-nn+z] == seq->a[z].v); + } + + for (z = nn; z < zn; z++) { + res->a[res->n++] = seq->a[z].v; + } + } else { + ///debug + for (z = zn - 1; z >= zn - ((int64_t)nn); z--) { + // fprintf(stderr, "[M::%s::z->%ld::idx->%ld] nn::%lu, zn::%ld, res->n::%ld, res>>1::%lu, res&1::%lu, seq>>1::%u, seq&1::%u\n", + // __func__, z, (int64_t)(res->n-nn+(zn-1-z)), nn, zn, (int64_t)(res->n), res->a[res->n-nn+(zn-1-z)]>>1, res->a[res->n-nn+(zn-1-z)]&1, seq->a[z].v>>1, seq->a[z].v&1); + assert(res->a[res->n-nn+(zn-1-z)] == (seq->a[z].v^1)); + } + + for (z = zn - 1 - ((int64_t)nn); z >= 0; z--) { + res->a[res->n++] = seq->a[z].v^1; + } + } + } else { + x = ulg_id(uidx->uovl, (bw>>1)); x <<= 1; x += (bw&1); + if(!is_bv_ul) {///if previous ul is also a raw utg node + assert(get_specfic_edge(raw->g, res->a[res->n-1], x)); + } else { + // assert(((res->a[res->n-1].v) == x)); + if(((res->a[res->n-1]) == x)) { + res->n--; + } else { + assert(get_specfic_edge(raw->g, res->a[res->n-1], x)); + } + } + nn = res->n + seq->n; zn = seq->n; + kv_resize(uint64_t, *res, nn); + if(!u_rev) { + for (z = 0; z < zn; z++) { + res->a[res->n++] = seq->a[z].v; + } + } else { + for (z = zn - 1; z >= 0; z--) { + res->a[res->n++] = seq->a[z].v^1; + } + } + } + } else { + nn = res->n + seq->n; zn = seq->n; + kv_resize(uint64_t, *res, nn); + if(!u_rev) { + for (z = 0; z < zn; z++) { + // fprintf(stderr, ">+<[M::%s::res->n->%ld] a>>1::%u, a&1::%u\n", + // __func__, (int64_t)res->n, seq->a[z].v>>1, seq->a[z].v&1); + res->a[res->n++] = seq->a[z].v; + } + } else { + for (z = zn - 1; z >= 0; z--) { + // fprintf(stderr, ">-<[M::%s::res->n->%ld] a>>1::%u, a&1::%u\n", + // __func__, (int64_t)res->n, seq->a[z].v>>1, (seq->a[z].v^1)&1); + res->a[res->n++] = seq->a[z].v^1; + } + } + } + } + + return res->n - rn; +} + +uint32_t get_integer_seq_ovlps(ul_resolve_t *uidx, uint64_t *p_a, int64_t p_n, int64_t match_bound, uint64_t skip_hom, asg64_v *res, uint64_t *r_w) +{ + ul_str_idx_t *str_idx = &(uidx->pstr); uint32_t v = p_a[match_bound]; uc_block_t *xi; + uint64_t *hid_a, hid_n, z, vz, ps, pe, ww[2], sw[2], occ; ul_str_t *str; int64_t s_n, s, p; + bubble_type *bub = uidx->bub; + + hid_a = str_idx->occ.a + str_idx->idx.a[v>>1]; + hid_n = str_idx->idx.a[(v>>1)+1] - str_idx->idx.a[v>>1]; + for (z = ww[0] = ww[1] = (*r_w) = occ = 0; z < hid_n; z++) { + str = &(str_idx->str.a[hid_a[z]>>32]); s_n = str->cn; + if(s_n < 2) continue; + vz = (uint32_t)(str->a[(uint32_t)hid_a[z]]); + assert((v>>1) == (vz>>1)); ps = pe = (uint64_t)-1; + sw[0] = sw[1] = 0; + if(v == vz) { + s = ((uint32_t)hid_a[z]) + 1; p = match_bound + 1; + for (; (s < s_n) && (p < p_n) && ((uint32_t)(str->a[s]) == p_a[p]); s++, p++) { + xi = &(uidx->idx->a[hid_a[z]>>32].bb.a[str->a[s]>>32]); + assert(((xi->hid<<1)+xi->rev)==((uint32_t)str->a[s])); + if(skip_hom && IF_HOM(xi->hid, *bub)) continue; + sw[1] += ug_occ_w(xi->ts, xi->te, &(uidx->l1_ug->u.a[xi->hid])); + } + if(s < s_n && p < p_n) continue; + if(p <= match_bound + 1) continue;///bridging the two sides of the breakpoint + if(sw[1] == 0) continue; + pe = p; + + s = (uint32_t)hid_a[z]; p = match_bound; + for (; (s >= 0) && (p >= 0) && ((uint32_t)(str->a[s]) == p_a[p]); s--, p--) { + xi = &(uidx->idx->a[hid_a[z]>>32].bb.a[str->a[s]>>32]); + assert(((xi->hid<<1)+xi->rev)==((uint32_t)str->a[s])); + if(skip_hom && IF_HOM(xi->hid, *bub)) continue; + sw[0] += ug_occ_w(xi->ts, xi->te, &(uidx->l1_ug->u.a[xi->hid])); + } + assert(p < match_bound); + if(s >= 0 && p >= 0) continue; + if(sw[0] == 0) continue; + ps = p + 1; + } else { + s = ((int32_t)((uint32_t)hid_a[z]))-1; p = match_bound + 1; + for (; (s >= 0) && (p < p_n) && ((uint32_t)(str->a[s]) == (p_a[p]^1)); s--, p++) { + xi = &(uidx->idx->a[hid_a[z]>>32].bb.a[str->a[s]>>32]); + assert(((xi->hid<<1)+xi->rev)==((uint32_t)str->a[s])); + if(skip_hom && IF_HOM(xi->hid, *bub)) continue; + sw[1] += ug_occ_w(xi->ts, xi->te, &(uidx->l1_ug->u.a[xi->hid])); + } + if(s >= 0 && p < p_n) continue; + if(p <= match_bound + 1) continue;///bridging the two sides of the breakpoint + if(sw[1] == 0) continue; + pe = p; + + s = (uint32_t)hid_a[z]; p = match_bound; + for (; (s < s_n) && (p >= 0) && ((uint32_t)(str->a[s]) == (p_a[p]^1)); s++, p--) { + xi = &(uidx->idx->a[hid_a[z]>>32].bb.a[str->a[s]>>32]); + assert(((xi->hid<<1)+xi->rev)==((uint32_t)str->a[s])); + if(skip_hom && IF_HOM(xi->hid, *bub)) continue; + sw[0] += ug_occ_w(xi->ts, xi->te, &(uidx->l1_ug->u.a[xi->hid])); + } + assert(p < match_bound); + if(s < s_n && p >= 0) continue; + if(sw[0] == 0) continue; + ps = p + 1; + } + if(res) kv_push(uint64_t, *res, ((ps<<32)|(pe))); + occ++; + ww[0] += sw[0]; ww[1] += sw[1]; + // if(match_bound == 15 && p_n == 36) { + // fprintf(stderr, "[M::%s::] z::%lu, ulid::%lu, ps::%lu, pe::%lu\n", __func__, z, hid_a[z]>>32, ps, pe); + // } + } + + (*r_w) = MIN(ww[0], ww[1]); + return occ; +} + +uint32_t usg_misjoin_topocut_aux(ul_resolve_t *uidx, ma_ug_t *ug, uint32_t v, uint32_t ref_v, uint32_t beg_v, +asg64_v *b, asg64_v *ub) +{ + #define cut_rate 0.49999 + uint32_t k, bn = b->n, is_del = 0, kk, ks, ke; + uint64_t n_ref, n_v, *a_ref, *a_v, n_min, n_ref_ov, n_v_ov, *a_ref_ov, *a_v_ov, l_occ, r_occ, ref_occ, v_occ, ww[2]; + ub->n = 0; kv_push(uint64_t, *ub, beg_v); + get_ul_path_info(uidx, ug, ref_v, NULL, NULL, NULL, NULL, NULL, ub); + // fprintf(stderr, "+++[M::%s::] start..., v>>1::%u, v&1::%u, ref_v>>1::%u, ref_v&1::%u, beg_v>>1::%u, beg_v&1::%u\n", __func__, v>>1, v&1, ref_v>>1, ref_v&1, beg_v>>1, beg_v&1); + n_ref = gen_ug_integer_seq_on_fly(uidx, ub->a, ub->n, b); + // fprintf(stderr, "+++[M::%s::] done...\n", __func__); + + ub->n = 0; kv_push(uint64_t, *ub, beg_v); + get_ul_path_info(uidx, ug, v, NULL, NULL, NULL, NULL, NULL, ub); + // fprintf(stderr, "---[M::%s::] start...\n", __func__); + n_v = gen_ug_integer_seq_on_fly(uidx, ub->a, ub->n, b); + // fprintf(stderr, "---[M::%s::] done...\n", __func__); + + a_ref = b->a + bn; a_v = b->a + bn + n_ref; n_min = MIN(n_ref, n_v); ub->n = 0; + for (k = 0; k < n_min && a_ref[k] == a_v[k]; k++); + + // if((beg_v>>1) == 5958) { + // fprintf(stderr, "[M::%s::] beg_v>>1::%u, beg_v&1::%u, v>>1::%u, v&1::%u, n_v::%lu, ref_v>>1::%u, ref_v&1::%u, n_ref::%lu, # prefix::%u\n", + // __func__, beg_v>>1, beg_v&1, v>>1, v&1, n_v, ref_v>>1, ref_v&1, n_ref, k); + // for (kk = 0; kk < n_v; kk++) { + // fprintf(stderr, "+[M::%s::] a_v[%u]>>1::%lu, a_v[%u]&1::%lu\n", + // __func__, kk, a_v[kk]>>1, kk, a_v[kk]&1); + // } + + // for (kk = 0; kk < n_ref; kk++) { + // fprintf(stderr, "+[M::%s::] a_ref[%u]>>1::%lu, a_ref[%u]&1::%lu\n", + // __func__, kk, a_ref[kk]>>1, kk, a_ref[kk]&1); + // } + // } + + + if(k > 0) { + n_ref_ov = get_integer_seq_ovlps(uidx, a_ref, n_ref, k - 1, 1, ub, &ww[0]); + n_v_ov = get_integer_seq_ovlps(uidx, a_v, n_v, k - 1, 1, ub, &ww[1]); + if(n_v_ov <= (n_ref_ov*cut_rate)) { + is_del = 1; + } else { + a_ref_ov = ub->a; a_v_ov = ub->a + n_ref_ov; kk = k; + + for (k = ref_occ = 0; k < n_ref_ov; k++) { + ks = a_ref_ov[k]>>32; ke = (uint32_t)a_ref_ov[k]; + l_occ = kk - ks; r_occ = ke - kk; + ref_occ += MIN(l_occ, r_occ); + } + + for (k = v_occ = 0; k < n_v_ov; k++) { + ks = a_v_ov[k]>>32; ke = (uint32_t)a_v_ov[k]; + l_occ = kk - ks; r_occ = ke - kk; + v_occ += MIN(l_occ, r_occ); + } + + if(v_occ <= (ref_occ*cut_rate)) is_del = 1; + } + } + b->n = bn; + return is_del; +} + + +uint32_t ulg_arc_cut_misjoin(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 *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, i, k, kv, nv, kw, nw, cnt = 0, n_vtx = g->n_seq<<1, mm_ul, ul_max, to_del; + asg_arc_t *av, *aw, *ve, *we; uint32_t ul_cnt, pb; + + b = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); + for (v = b->n = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + + if(g->seq_vis[v] == 0) { + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + if (nv < 2) continue; + + for (i = kv = 0, ul_max = -1; i < nv; ++i) { + if(av[i].del) continue; kv++; + get_ul_path_info(uidx, ug, av[i].v, NULL, NULL, &ul_cnt, NULL, NULL, NULL); + if(ul_cnt > ul_max) ul_max = ul_cnt; + } + // if((v>>1) == 5958) { + // fprintf(stderr, "[M::%s::] v>>1::%u, v&1::%u, kv::%u, ul_max::%u\n", + // __func__, v>>1, v&1, kv, ul_max); + // } + if(kv < 2 || ul_max == 0) continue; //note: ul_cnt/ul_max might be 0 + + for (i = 0; i < nv; ++i) { + if(av[i].del) continue; + get_ul_path_info(uidx, ug, av[i].v, NULL, NULL, &ul_cnt, NULL, NULL, NULL); + // if((v>>1) == 5958) { + // fprintf(stderr, "+++[M::%s::] v>>1::%u, v&1::%u, av[i].v>>1::%u, av[i].v&1::%u, ul_cnt::%u\n", + // __func__, v>>1, v&1, av[i].v>>1, av[i].v&1, ul_cnt); + // } + if(ul_cnt > 0) ul_cnt = ul_max/ul_cnt; + else ul_cnt = ul_max<<2; + ul_cnt = ((uint32_t)-1) - ul_cnt; + kv_push(uint64_t, *b, (((uint64_t)ul_cnt)<<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); + if(nv<=1 && asg_arc_n(g, w) <= 1) continue; + + ve = &(g->arc[(uint32_t)b->a[k]]); pb = b->n; ub->n = 0; + get_ul_path_info(uidx, ug, ve->v, NULL, NULL, &mm_ul, NULL, NULL, NULL); + // if((v>>1) == 5958 && (w>>1) == 11105) { + // fprintf(stderr, ">>>[M::%s::] v>>1::%u, v&1::%u, w>>1::%u, w&1::%u, mm_ul::%u\n", + // __func__, v>>1, v&1, w>>1, w&1, mm_ul); + // } + if(mm_ul == 0) continue;///no UL read support this path + + for (i = kv = 0; i < nv; ++i) { + if(av[i].del) continue; + kv++; + if(av[i].v == ve->v) continue; + get_ul_path_info(uidx, ug, av[i].v, NULL, NULL, &ul_cnt, NULL, NULL, NULL); + // if((v>>1) == 5958 && (w>>1) == 11105) { + // fprintf(stderr, "*[M::%s::i->%u] v>>1::%u, v&1::%u, w>>1::%u, w&1::%u, av[i].v>>1::%u, av[i].v&1::%u, ul_cnt::%u\n", + // __func__, i, v>>1, v&1, w>>1, w&1, av[i].v>>1, av[i].v&1, ul_cnt); + // } + if (mm_ul <= ul_cnt*len_rat) { + ul_cnt = ((uint32_t)-1) - ul_cnt; + kv_push(uint64_t, *b, (((uint64_t)ul_cnt)<<32)|((uint64_t)(i))); + } + } + if(b->n == pb) continue; + to_del = 0; assert(kv >= 2); + + aw = asg_arc_a(g, w); nw = asg_arc_n(g, w); + for (i = kw = 0, we = NULL; i < nw; ++i) { + if (aw[i].del) continue; + if (aw[i].v == (v^1)) we = &(aw[i]); + kw++; + } + + if((kv > 1 && kw > 1) || (usg_bridge_topocut_aux(uidx, ug, ve->v, max_ext, max_ext_hifi, topo_level, 0, b, ub))) { + radix_sort_srt64(b->a + pb, b->a + b->n); + for (i = pb; i < b->n; i++) { + // if((v>>1) == 5958 && (w>>1) == 11105) { + // fprintf(stderr, "#[M::%s::srt_i->%u] v>>1::%u, v&1::%u, w>>1::%u, w&1::%u, av[(uint32_t)b->a[i]].v>>1::%u, av[(uint32_t)b->a[i]].v&1::%u\n", + // __func__, (uint32_t)b->a[i], v>>1, v&1, w>>1, w&1, av[(uint32_t)b->a[i]].v>>1, av[(uint32_t)b->a[i]].v&1); + // } + if(usg_misjoin_topocut_aux(uidx, ug, ve->v, av[(uint32_t)b->a[i]].v, v, b, ub)) { + to_del = 1; + break; + } + } + } + + b->n = pb; + if(to_del) { + if(kv > 1 && kw > 1) { + ve->del = we->del = 1; ++cnt; + } else { + pb = b->n; + get_ul_path_info(uidx, ug, ve->v, NULL, NULL, NULL, NULL, NULL, b); + for (i = pb; i < b->n; i++) ulg_seq_del(ug, (b->a[i]>>1)); + cnt += b->n - pb; b->n = pb; + } + } + } + + + if(!in) free(tx.a); if(!ib) free(tb.a); + if (cnt > 0) asg_cleanup(g); + return cnt; +} + +void get_ul_arc_supports(ul_resolve_t *uidx, asg_arc_t *ve, asg64_v *b_int, asg64_v *b_raw, uint64_t skip_hom, uint64_t *retrun_w_v, uint64_t *retrun_w_r) +{ + ma_ug_t *iug = uidx->uovl.i_ug; asg_t *g = iug->g; uint32_t v, k, z, zn, nv, n_pre, l_v, l_r; + asg_arc_t *av; v = ve->ul>>32; uint32_t b_int_s = b_int->n, b_raw_s = b_raw->n; + (*retrun_w_v) = (*retrun_w_r) = (uint64_t)-1; uint64_t *raw_v, *raw_r, w_v, w_r, min_w_v, min_w_r; + + get_ul_path_info(uidx, iug, v^1, NULL, NULL, NULL, NULL, NULL, b_int); nv = b_int->n-b_int_s; + for (k = 0, n_pre = b_int->n; k < (nv>>1); k++) { + v = b_int->a[k+b_int_s]; + b_int->a[k+b_int_s] = b_int->a[b_int_s+nv-k-1]^1; + b_int->a[b_int_s+nv-k-1] = v^1; + } + if(nv&1) b_int->a[k+b_int_s] ^= 1; + v = ve->ul>>32; + + get_ul_path_info(uidx, iug, ve->v, NULL, NULL, NULL, NULL, NULL, b_int); + // if(((ve->ul>>32) == 24093) && (ve->v == 61472)) { + // for (z = 0; z < b_int->n-b_int_s; z++) { + // fprintf(stderr, "[M::%s::] integ_v[%u]>>1:%lu, integ_v[%u]&1:%lu\n", __func__, + // z, b_int->a[b_int_s+z]>>1, z, b_int->a[b_int_s+z]&1); + // } + // } + l_v = gen_ug_integer_seq_on_fly(uidx, b_int->a+b_int_s, b_int->n-b_int_s, b_raw); + + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + for (k = 0, min_w_v = min_w_r = (uint64_t)-1; k < nv; k++) { + if(av[k].del || av[k].v == ve->v) continue; + + b_int->n = n_pre; b_raw->n = b_raw_s + l_v; + get_ul_path_info(uidx, iug, av[k].v, NULL, NULL, NULL, NULL, NULL, b_int); + // if(((ve->ul>>32) == 24093) && (ve->v == 61472)) { + // for (z = 0; z < b_int->n-b_int_s; z++) { + // fprintf(stderr, "[M::%s::k->%u] integ_r[%u]>>1:%lu, integ_r[%u]&1:%lu\n", __func__, + // k, z, b_int->a[b_int_s+z]>>1, z, b_int->a[b_int_s+z]&1); + // } + // } + + l_r = gen_ug_integer_seq_on_fly(uidx, b_int->a+b_int_s, b_int->n-b_int_s, b_raw); + + raw_v = b_raw->a + b_raw_s; raw_r = b_raw->a + b_raw_s + l_v; zn = MIN(l_v, l_r); + // if((v>>1) == 77 && (ve->v>>1) == 78) { + // fprintf(stderr, "[M::%s::] l_v::%u\n", __func__, l_v); + // for (z = 0; z < l_v; z++) { + // fprintf(stderr, "raw_v[%u]>>1::utg%.6dl, raw_v[%u]&1::%lu\n", + // z, (int32_t)(raw_v[z]>>1)+1, z , raw_v[z]&1); + // } + // fprintf(stderr, "[M::%s::] l_r::%u\n", __func__, l_r); + // for (z = 0; z < l_r; z++) { + // fprintf(stderr, "raw_r[%u]>>1::utg%.6dl, raw_r[%u]&1::%lu\n", + // z, (int32_t)(raw_r[z]>>1)+1, z, raw_r[z]&1); + // } + // } + for (z = 0; z < zn && raw_v[z] == raw_r[z]; z++); + assert(z > 0); + get_integer_seq_ovlps(uidx, raw_v, l_v, z - 1, skip_hom, NULL, &w_v); + get_integer_seq_ovlps(uidx, raw_r, l_r, z - 1, skip_hom, NULL, &w_r); + // if((v>>1) == 77 && (ve->v>>1) == 78) { + // fprintf(stderr, "[M::%s::] l_v::%u, l_r::%u, z::%u, w_v::%lu, w_r::%lu\n", __func__, l_v, l_r, z, w_v, w_r); + // } + if((min_w_v == (uint64_t)-1) || (min_w_v > w_v) || (min_w_v == w_v && min_w_r < w_r)) { + min_w_v = w_v; min_w_r = w_r; + } + } + + (*retrun_w_v) = min_w_v; (*retrun_w_r) = min_w_r; +} + +static void worker_update_ul_arc_supports(void *data, long i, int tid) // callback for kt_for() +{ + ul_resolve_t *uidx = (ul_resolve_t *)data; integer_t *buf = &(uidx->str_b.buf[tid]); + uint64_t *x = &(uidx->uovl.iug_tra->a[i]); asg_arc_t *ve = &(uidx->uovl.i_ug->g->arc[*x]); + asg64_v b_v, b_r; uint64_t w_v, w_r; + b_v.a = buf->u.a; b_v.n = buf->u.n; b_v.m = buf->u.m; + b_r.a = buf->o.a; b_r.n = buf->o.n; b_r.m = buf->o.m; + + b_v.n = b_r.n = 0; + get_ul_arc_supports(uidx, ve, &b_v, &b_r, 1, &w_v, &w_r); + + buf->u.a = b_v.a; buf->u.n = b_v.n; buf->u.m = b_v.m; + buf->o.a = b_r.a; buf->o.n = b_r.n; buf->o.m = b_r.m; + + (*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) +{ + asg64_v tx = {0,0,0}, tb = {0,0,0}, *b = NULL, *ub = NULL; asg_t *g = ug->g; + uint32_t v, w, i, k, kv, nv, kw, nw, cnt = 0, n_vtx = g->n_seq<<1, to_del; + asg_arc_t *av, *aw, *ve, *we; uint64_t w_q, w_t, pb; + + b = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); + for (v = b->n = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + + if(g->seq_vis[v] == 0) { + 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; + kv_push(uint64_t, *b, ((uint64_t)(av-g->arc+i))); + } + } + } + + // fprintf(stderr, "\n#[M::%s::] Starting...\n", __func__); + uidx->uovl.iug_tra = b; + kt_for(uidx->str_b.n_thread, worker_update_ul_arc_supports, uidx, b->n);///all ul + ug + uidx->uovl.iug_tra = NULL; + // fprintf(stderr, "#[M::%s::] Done\n", __func__); + + 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; + ve = &(g->arc[(uint32_t)b->a[k]]); + for (i = kv = 0; i < nv; ++i) { + if(av[i].del) continue; + kv++; + } + for (i = kw = 0, we = NULL; i < nw; ++i) { + if (aw[i].del) continue; + if (aw[i].v == (v^1)) we = &(aw[i]); + kw++; + } + if(kv <= 1 && kw <= 1) continue; + + // if((v>>1) == 78 && (w>>1) == 77) { + // fprintf(stderr, "\n#[M::%s::] v>>1::%u, v&1::%u, kv::%u, w>>1::%u, w&1::%u, kw::%u\n", + // __func__, v>>1, v&1, kv, w>>1, w&1, kw); + // } + + if(kv > 1) { + 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((v>>1) == 78 && (w>>1) == 77) { + // fprintf(stderr, "+[M::%s::] v>>1::%u, v&1::%u, kv::%u, w>>1::%u, w&1::%u, kw::%u, w_q::%lu, w_t::%lu\n", + // __func__, v>>1, v&1, kv, w>>1, w&1, kw, w_q, w_t); + // } + if(w_q == (uint64_t)-1) continue; + if(w_q > w_t*len_rat) continue; + } + + if(kw > 1) { + pb = b->n; ub->n = 0; + get_ul_arc_supports(uidx, we, b, ub, skip_hom, &w_q, &w_t); + b->n = pb; ub->n = 0; + // if((v>>1) == 78 && (w>>1) == 77) { + // fprintf(stderr, "-[M::%s::] v>>1::%u, v&1::%u, kv::%u, w>>1::%u, w&1::%u, kw::%u, w_q::%lu, w_t::%lu\n", + // __func__, v>>1, v&1, kv, w>>1, w&1, kw, w_q, w_t); + // } + if(w_q == (uint64_t)-1) continue; + if(w_q > w_t*len_rat) continue; + } + + to_del = 0; + if(topo_level == 0) { + to_del = 1; + } else if(topo_level == 2) { + if (kv > 1 && kw > 1) to_del = 1; + } else { + if (kv > 1 && kw > 1) { + to_del = 1; + } else if (kw == 1) { + if (usg_topocut_aux(uidx, ug, w^1, max_ext, max_ext_hifi, b, ub)) to_del = 1; + } else if (kv == 1) { + if (usg_topocut_aux(uidx, ug, v^1, max_ext, max_ext_hifi, b, ub)) to_del = 1; + } + } + + if (to_del) { + ve->del = we->del = 1, ++cnt; + } + } + + if(!in) free(tx.a); if(!ib) free(tb.a); + if (cnt > 0) asg_cleanup(g); + return cnt; +} + +uint64_t ulg_bub_pop_cut_aux(ul_resolve_t *uidx, ma_ug_t *ug, uint32_t v0, buf_t *x, uint32_t max_ext, uint32_t max_ext_hifi, asg64_v *b, asg64_v *ub) +{ + uint32_t i, v, u, bn = b->n, n_ext, ul_occ, r = 0; binfo_t *t; + v = x->S.a[0]; + do { + u = x->a[v].p; // u->v + x->a[v].d = (uint32_t)-1; + v = u; + } while (v != v0); + + for (i = n_ext = 0; i < x->b.n; ++i) { // clear the states of visited vertices + t = &x->a[x->b.a[i]]; + if(t->d == (uint32_t)-1) continue; + get_iug_u_raw_occ(uidx, x->b.a[i]>>1, &ul_occ, NULL); + n_ext += ul_occ; kv_push(uint64_t, *b, x->b.a[i]); + } + + 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; + return r; +} + +uint32_t ulg_bub_pop_backtrack(ma_ug_t *ug, uint32_t v0, buf_t *b) +{ + uint32_t i, v, u, cnt = b->e.n; asg_t *g = ug->g; + + ///b->S.a[0] is the sink of this bubble + for (i = 0; i < b->b.n; ++i) g->seq[b->b.a[i]>>1].del = 1; + + ///remove all edges (self/reverse for each edge) in this bubble + for (i = 0; i < b->e.n; ++i) { + g->arc[b->e.a[i]].del = 1; + asg_arc_del(g, g->arc[b->e.a[i]].v^1, (g->arc[b->e.a[i]].ul>>32)^1, 1); + } + + ///v is the sink of this bubble + v = b->S.a[0]; + do { + u = b->a[v].p; // u->v + g->seq[v>>1].del = 0; + asg_arc_del(g, u, v, 0); + asg_arc_del(g, v^1, u^1, 0); + cnt--; + v = u; + } while (v != v0); + + for (i = 0; i < b->b.n; ++i) { + if(!g->seq[b->b.a[i]>>1].del) continue; + ulg_seq_del(ug, (b->b.a[i]>>1)); + } + + return cnt; +} + +uint64_t ulg_bub_pop1(ul_resolve_t *uidx, ma_ug_t *ug, uint32_t v0, uint64_t max_dist, buf_t *x, +uint32_t max_ext, uint32_t max_ext_hifi, uint32_t check_bubble_only, uint32_t is_pop, uint32_t skip_hom, +asg64_v *b, asg64_v *ub, uint32_t *r_w_c, uint32_t *r_w_m) +{ + asg_t *g = ug->g; uint32_t pb = b->n; uint64_t w_q, w_t, wc, wm, ww; (*r_w_c) = (*r_w_m) = (uint32_t)-1; + uint32_t v, w, i, kv, nv, kw, cnt = 0, fail_b = 0, n_tips = 0, tip_end = (uint32_t)-1; + uint32_t l, d, c, m, n_pending = 0, z, to_replace; asg_arc_t *av, *ve, *we; binfo_t *t; + if(g->seq[v0>>1].del || get_arcs(g, v0, NULL, 0) < 2) return 0; // already deleted + // fprintf(stderr, "sbsbsbsbsbsbsb[M::%s::(v0>>1)->%u::(v0&1)->%u] check_bubble_only::%u, is_pop::%u\n", + // __func__, v0>>1, v0&1, check_bubble_only, is_pop); + + x->S.n = x->T.n = x->b.n = x->e.n = 0; + x->a[v0].c = x->a[v0].d = x->a[v0].m = x->a[v0].nc = x->a[v0].np = 0; + kv_push(uint32_t, x->S, v0); + + do { + v = kv_pop(x->S); d = x->a[v].d; c = x->a[v].c; m = x->a[v].m; + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); kv = get_arcs(g, v, NULL, 0); + for (i = 0; i < nv; ++i) { + if (av[i].del) continue; + w = av[i].v; t = &(x->a[w]); l = ((v == v0)?(0):((uint32_t)av[i].ul)); + if ((w>>1) == (v0>>1)) { + fail_b = 1; + break; + } + kv_push(uint32_t, x->e, (g->idx[v]>>32) + i); ///for backtracking + if (d + l > max_dist) { + fail_b = 1; + break; + } + kw = get_arcs(g, w^1, NULL, 0); + wc = 0; wm = get_ul_occ(uidx, w>>1); + if(!check_bubble_only) { + if(kv > 1) { + pb = b->n; ub->n = 0; ve = &(av[i]); + // fprintf(stderr, "[M::%s::] ve->v>>1::%lu, ve->v&1::%lu, ve->w>>1::%u, ve->w&1::%u\n", + // __func__, ve->ul>>33, (ve->ul>>32)&1, ve->v>>1, ve->v&1); + get_ul_arc_supports(uidx, ve, b, ub, skip_hom, &w_q, &w_t); + b->n = pb; ub->n = 0; + if(w_q < w_t) { + if(w_q == 0) { + ww = 10; + } else { + if(w_t != 0) ww = ((w_t - w_q)*10)/w_t; + else ww = 0; + } + if(ww > wc) wc = ww; + } + } + + if(kw > 1) { + pb = b->n; ub->n = 0; we = get_specfic_edge(g, w^1, v^1); assert(we); + // fprintf(stderr, "[M::%s::] we->v>>1::%lu, we->v&1::%lu, we->w>>1::%u, we->w&1::%u\n", + // __func__, we->ul>>33, (we->ul>>32)&1, we->v>>1, we->v&1); + get_ul_arc_supports(uidx, we, b, ub, skip_hom, &w_q, &w_t); + b->n = pb; ub->n = 0; + if(w_q < w_t) { + if(w_q == 0) { + ww = 10; + } else { + if(w_t != 0) ww = ((w_t - w_q)*10)/w_t; + else ww = 0; + } + if(ww > wc) wc = ww; + } + } + } + + if (t->s == 0) { + kv_push(uint32_t, x->b, w); + t->p = v, t->s = 1, t->d = d + l; + t->c = c + wc; t->m = m + wm; + t->r = kw; + ++n_pending; + } else { + to_replace = 0; + if((c + wc) < t->c) { + to_replace = 1; + } else if(((c + wc) == t->c) && (m + wm > t->m)) { + to_replace = 1; + } else if(((c + wc) == t->c) && (m + wm == t->m) && (d + l > t->d)) { + to_replace = 1; + } + if(to_replace) { + t->p = v; t->c = c + wc; t->m = m + wm; + } + if (d + l < t->d) t->d = d + l; // update dist + } + + if (--(t->r) == 0) { + z = get_arcs(g, w, NULL, 0); + if(z > 0) { + kv_push(uint32_t, x->S, w); + } + else { + ///at most one tip + if(n_tips != 0) { + fail_b = 1; + break; + } + n_tips++; tip_end = w; + } + --n_pending; + } + } + if(fail_b) break; + if(n_tips == 1) { + if(tip_end != (uint32_t)-1 && n_pending == 0 && x->S.n == 0) { + kv_push(uint32_t, x->S, tip_end); + break; + } + fail_b = 1; + break; + } + + if (i < nv || x->S.n == 0) { + fail_b = 1; + break; + } + } while (x->S.n > 1 || n_pending); + + if(!fail_b) {//there is a bubble + cnt = 1; (*r_w_c) = x->a[x->S.a[0]].c; (*r_w_m) = x->a[x->S.a[0]].m; + if(!check_bubble_only) { + if(is_pop) { + cnt = ulg_bub_pop_backtrack(ug, v0, x); + } else { + cnt = ulg_bub_pop_cut_aux(uidx, ug, v0, x, max_ext, max_ext_hifi, b, ub); + } + } + } + for (i = 0; i < x->b.n; ++i) { // clear the states of visited vertices + // if(v0 == 119) { + // fprintf(stderr, "-[M::%s::(v0>>1)->%u::(v0&1)->%u] v[%u]>>1::%u, v[%u]&1::%u, cnt::%u, is_pop::%u\n", + // __func__, v0>>1, v0&1, i, x->b.a[i]>>1, i, x->b.a[i]&1, cnt, is_pop); + // } + t = &x->a[x->b.a[i]]; + t->s = t->c = t->d = t->m = t->nc = t->np = 0; + } + if(!cnt) (*r_w_c) = (*r_w_m) = (uint32_t)-1; + return cnt; +} + +uint64_t ulg_pop_bubble(ul_resolve_t *uidx, ma_ug_t *ug, uint64_t* i_max_dist, uint32_t max_ext, uint32_t max_ext_hifi, uint32_t skip_hom, asg64_v *in, asg64_v *ib) +{ + fprintf(stderr, "[M::%s::] Starting...\n", __func__); + asg_t *g = ug->g; + asg64_v tx = {0,0,0}, tb = {0,0,0}, *ob = NULL, *ub = NULL; + uint32_t v, w, n_vtx = g->n_seq<<1, n_arc, nv, i, wc[2], wm[2], mm_c, mm_m, mm_v; + uint64_t n_pop = 0, max_dist; + asg_arc_t *av = NULL; if (!g->is_symm) asg_symm(g); + buf_t b; memset(&b, 0, sizeof(buf_t)); + b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); + if(i_max_dist) max_dist = (*i_max_dist); + else max_dist = get_bub_pop_max_dist_advance(g, &b); + ob = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); ob->n = ub->n = 0; + + if(max_dist > 0) { + for (v = 0; v < n_vtx; ++v) { + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + if (nv < 2 || g->seq[v>>1].del) continue; + for (i = n_arc = 0; i < nv; ++i) { + if (!av[i].del) ++n_arc; + } + if (n_arc < 2) continue; + ///find a bubble + ob->n = ub->n = 0; + if(ulg_bub_pop1(uidx, ug, v, max_dist, &b, max_ext, max_ext_hifi, 1, 0, skip_hom, ob, ub, &(wc[0]), &(wm[0]))) { + w = b.S.a[0]^1; mm_c = mm_m = mm_v = (uint32_t)-1; + + ob->n = ub->n = 0; + if(ulg_bub_pop1(uidx, ug, v, max_dist, &b, max_ext, max_ext_hifi, 0, 0, skip_hom, + ob, ub, &(wc[0]), &(wm[0]))) { + mm_c = wc[0]; mm_m = wm[0]; mm_v = v; + } + + ob->n = ub->n = 0; + if(ulg_bub_pop1(uidx, ug, w, max_dist, &b, max_ext, max_ext_hifi, 0, 0, skip_hom, + ob, ub, &(wc[1]), &(wm[1]))) { + if((wc[1] < mm_c) && (wc[1] == mm_c && wm[1] > mm_m)) { + mm_c = wc[1]; mm_m = wm[1]; mm_v = w; + } + } + + if(mm_v != (uint32_t)-1) { + ob->n = ub->n = 0; + w = ulg_bub_pop1(uidx, ug, mm_v, max_dist, &b, max_ext, max_ext_hifi, 0, 1, skip_hom, ob, ub, &(wc[0]), &(wm[0])); + assert(w); + n_pop += w; + } + } + } + } + + free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); + if(n_pop) asg_cleanup(g); + if(!in) free(tx.a); if(!ib) free(tb.a); + fprintf(stderr, "[M::%s::] Done...\n", __func__); + return n_pop; +} + +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}; + ma_ug_t *iug = idx->i_ug; int64_t i, mm_tip = ulopt->max_tip; uint64_t cnt = 1, topo_level; + double step = (ulopt->clean_round==1?ulopt->max_ovlp_drop_ratio: + ((ulopt->max_ovlp_drop_ratio-ulopt->min_ovlp_drop_ratio)/(ulopt->clean_round-1))); + double drop = ulopt->min_ovlp_drop_ratio; CALLOC(iug->g->seq_vis, iug->g->n_seq*2); + + // ulg_arc_cut_tips(uidx, iug, ulopt->max_tip, ulopt->max_tip_hifi, &bu, &uu); + + for (i = 0, drop = ulopt->min_ovlp_drop_ratio; i < ulopt->clean_round; i++, drop += step) { + fprintf(stderr, "\n[M::%s::] Starting round-%ld, drop::%f\n", __func__, i, drop); + cnt = 1; topo_level = 2; mm_tip = ulopt->max_tip; + while (cnt) { + cnt = 0; + asg_arc_identify_simple_bubbles_multi(iug->g, ulopt->b_mask_t, 0); + cnt += ulg_arc_cut_supports(uidx, iug, mm_tip, ulopt->max_tip_hifi, drop, ulopt->is_trio, topo_level, 1, NULL, &bu, &uu); + cnt += ulg_arc_cut_tips(uidx, iug, mm_tip, ulopt->max_tip_hifi, &bu, &uu); + } + + cnt = 1; topo_level = 1; mm_tip = ulopt->max_tip; + while (cnt) { + cnt = 0; + asg_arc_identify_simple_bubbles_multi(iug->g, ulopt->b_mask_t, 0); + cnt += ulg_arc_cut_supports(uidx, iug, mm_tip, ulopt->max_tip_hifi, drop, ulopt->is_trio, topo_level, 1, NULL, &bu, &uu); + cnt += ulg_arc_cut_tips(uidx, iug, mm_tip, ulopt->max_tip_hifi, &bu, &uu); + } + + cnt = 1; topo_level = 1; mm_tip = ((int64_t)0x7fffffff); + while (cnt) { + cnt = 0; + asg_arc_identify_simple_bubbles_multi(iug->g, ulopt->b_mask_t, 0); + cnt += ulg_arc_cut_supports(uidx, iug, mm_tip, ulopt->max_tip_hifi, drop, ulopt->is_trio, topo_level, 1, NULL, &bu, &uu); + cnt += ulg_arc_cut_tips(uidx, iug, mm_tip, ulopt->max_tip_hifi, &bu, &uu); + } + fprintf(stderr, "[M::%s::] Done round-%ld, drop::%f\n", __func__, i, drop); + } + + ulg_pop_bubble(uidx, iug, NULL, ((int64_t)0x7fffffff), ulopt->max_tip_hifi, 1, &bu, &uu); + + // 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; + // cnt += ulg_arc_cut_occ(uidx, iug, mm_tip, ulopt->max_tip_hifi, ulopt->is_trio, 2, &bu, &uu); + // asg_arc_identify_simple_bubbles_multi(iug->g, ulopt->b_mask_t, 0); + // cnt += ulg_arc_cut_length(uidx, iug, mm_tip, ulopt->max_tip_hifi, drop/**MIN(drop, ulopt->hom_check_drop_rate)**/, ulopt->is_trio, 2, 1, NULL, &bu, &uu); + // cnt += ulg_arc_cut_tips(uidx, iug, mm_tip, ulopt->max_tip_hifi, &bu, &uu); + + // asg_arc_identify_simple_bubbles_multi(iug->g, ulopt->b_mask_t, 0); + // cnt += ulg_arc_cut_bridge(uidx, iug, mm_tip, ulopt->max_tip_hifi, 0.5, ulopt->is_trio, 2, NULL, &bu, &uu); + + // asg_arc_identify_simple_bubbles_multi(iug->g, ulopt->b_mask_t, 0); + // cnt += ulg_arc_cut_misjoin(uidx, iug, mm_tip, ulopt->max_tip_hifi, 0.5, ulopt->is_trio, 2, NULL, &bu, &uu); + // } + + // for (i = cnt = 0, mm_tip = ((int64_t)0x7fffffff), 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; + // cnt += ulg_arc_cut_occ(uidx, iug, mm_tip, ulopt->max_tip_hifi, ulopt->is_trio, 2, &bu, &uu); + // asg_arc_identify_simple_bubbles_multi(iug->g, ulopt->b_mask_t, 0); + // cnt += ulg_arc_cut_length(uidx, iug, mm_tip, ulopt->max_tip_hifi, drop/**MIN(drop, ulopt->hom_check_drop_rate)**/, ulopt->is_trio, 2, 1, NULL, &bu, &uu); + // cnt += ulg_arc_cut_tips(uidx, iug, mm_tip, ulopt->max_tip_hifi, &bu, &uu); + + // asg_arc_identify_simple_bubbles_multi(iug->g, ulopt->b_mask_t, 0); + // cnt += ulg_arc_cut_bridge(uidx, iug, mm_tip, ulopt->max_tip_hifi, 0.5, ulopt->is_trio, 2, NULL, &bu, &uu); + + // asg_arc_identify_simple_bubbles_multi(iug->g, ulopt->b_mask_t, 0); + // cnt += ulg_arc_cut_misjoin(uidx, iug, mm_tip, ulopt->max_tip_hifi, 0.5, ulopt->is_trio, 2, NULL, &bu, &uu); + // } + // } + + free(bu.a); free(uu.a); +} + +void clc_contain(ul_resolve_t *uidx, uint64_t id, uint64_t is_ul, integer_t *buf) +{ + ma_ug_t *raw = uidx->l1_ug; ma_utg_t *ru; ug_opt_t *uopt = uidx->uopt; + uint64_t k, m, qn, tn; ma_hit_t_alloc *x; asg_arc_t e; ul2ul_item_t *o; + int64_t min_ovlp = uopt->min_ovlp, max_hang = uopt->max_hang, r; + + if(is_ul) { + o = get_ul_ovlp(&(uidx->uovl), id, 1); + assert(o && (!o->is_del)); + for (k = 0; k < o->cn; k++) { + if((o->a[k].is_del) || (!ulg_type(uidx->uovl, o->a[k].hid))) continue; + r = integer_hit2arc(&(o->a[k]), ulg_len(*uidx, id), ulg_len(*uidx, o->a[k].hid), + ulg_occ(*uidx, id), ulg_occ(*uidx, o->a[k].hid), id, o->a[k].hid, 0, NULL); + if(r != MA_HT_TCONT) continue; + kv_push(uint64_t, buf->o, o->a[k].hid); + } + } else { + ru = &(raw->u.a[id]); + for (k = 0; k < ru->n; k++) { + x = &(uopt->sources[ru->a[k]>>33]); + for (m = 0; m < x->length; m++) { + qn = Get_qn(x->buffer[m]); tn = Get_tn(x->buffer[m]); + r = ma_hit2arc(&(x->buffer[m]), Get_READ_LENGTH(R_INF, qn), Get_READ_LENGTH(R_INF, tn), + max_hang, asm_opt.max_hang_rate, min_ovlp, &e); + if(r != MA_HT_TCONT) continue; + kv_push(uint64_t, buf->u, tn); + } + } + } +} + +static void worker_renew_u2g_cov(void *data, long i, int tid) // callback for kt_for() +{ + ul_resolve_t *uidx = (ul_resolve_t *)data; + integer_t *buf = &(uidx->str_b.buf[tid]); + ma_utg_t *iu = &(uidx->uovl.i_ug->u.a[i]); + uint64_t z, ri, k, ul_n, ug_n, *a, a_n; ma_ug_t *raw = uidx->l1_ug; + ul_cov_t *idx = &(uidx->uovl.cc); + + idx->uc[i] = idx->raw_uc[i] = idx->hc[i] = ul_n = ug_n = 0; buf->u.n = buf->o.n = 0; + for (z = 0; z < iu->n; z++) { + ri = iu->a[z]>>33; + if(ulg_type(uidx->uovl, ri)) {///ul + clc_contain(uidx, ulg_id(uidx->uovl, ri), 1, buf); ul_n++; + } else {///ug node + clc_contain(uidx, ulg_id(uidx->uovl, ri), 0, buf); + ug_n += raw->u.a[ulg_id(uidx->uovl, ri)].n; + } + } + idx->raw_uc[i] = ul_n; + + ///buf->o: UL contained reads + ///buf->u: HiFi reads + if(buf->o.n > 0) { + a = buf->o.a; a_n = buf->o.n; + radix_sort_srt64(a, a + a_n); + for (z = 0, k = 1; k <= a_n; k++) { + if(k == a_n || a[z] != a[k]) { + ul_n++; z = k; + } + } + } + if(buf->u.n > 0) { + a = buf->u.a; a_n = buf->u.n; + radix_sort_srt64(a, a + a_n); + for (z = 0, k = 1; k <= a_n; k++) { + if(k == a_n || a[z] != a[k]) { + ug_n++; z = k; + } + } + } + idx->uc[i] = ul_n; idx->hc[i] = ug_n; +} + +ul2ul_t* get_ulg_spec_ovlp(ul_resolve_t *uidx, uint64_t uv, uint64_t uw) +{ + ma_ug_t *iug = uidx->uovl.i_ug; uint32_t qid, tid; + qid = iug->u.a[uv>>1].a[((uv&1)?(0):(iug->u.a[uv>>1].n-1))]>>33; + tid = iug->u.a[uw>>1].a[((uw&1)?(iug->u.a[uw>>1].n-1):(0))]>>33; + + if(uidx->uovl.item_idx[qid] == (uint32_t)-1) return NULL; + ul2ul_item_t *o = &(uidx->uovl.a[uidx->uovl.item_idx[qid]]); uint64_t k; + for (k = 0; k < o->cn; k++) { + if(o->a[k].hid != tid) continue; + return &(o->a[k]); + } + return NULL; +} + + + + + +void gen_raw_ug_seq(ul_resolve_t *uidx, ul_str_t *str, ma_utg_t *u, ma_ug_t *raw, uinfo_srt_warp_t *res, uint64_t iug_id) +{ + uint64_t k, cd, nd, rev, x, os, oe; uinfo_srt_t *p; ma_utg_t *ru; + uc_block_t *xi; ul2ul_t *z = NULL; ul_str_t *c_str; int64_t m, cs, ce; + for (k = 0, res->n = 0; k < u->n; k++) { + cd = u->a[k]>>33; rev = ((u->a[k]>>32)&1); ru = NULL; c_str = NULL; + if(ulg_type(uidx->uovl, cd)) c_str = &(str[ulg_id(uidx->uovl, cd)]); ///ul + else ru = &(raw->u.a[ulg_id(uidx->uovl, cd)]); //ug + + if(c_str) { + if(k + 1 < u->n) { + nd = u->a[k+1]>>33; + z = get_ul_spec_ovlp(&(uidx->uovl), cd, nd); + assert(z && (!z->is_del)); + if(!rev) { + cs = 0; ce = z->qs_k + 1; + } else { + cs = z->qe_k - 1; ce = c_str->cn; + } + } else { + cs = 0; ce = c_str->cn; + } + + if(!rev) { + xi = &(uidx->idx->a[cd].bb.a[c_str->a[cs]>>32]); + x = (uint32_t)c_str->a[cs]; + } else { + xi = &(uidx->idx->a[cd].bb.a[c_str->a[ce-1]>>32]); + x = (uint32_t)c_str->a[ce-1]; x ^= 1; + } + + if(res->n) { + assert((res->a[res->n-1].v) == x); + os = MAX(xi->ts, res->a[res->n-1].s); oe = MIN(xi->te, res->a[res->n-1].e); + assert(oe > os); + os = MIN(xi->ts, res->a[res->n-1].s); oe = MAX(xi->te, res->a[res->n-1].e); + res->a[res->n-1].s = os; res->a[res->n-1].e = oe; + if(!rev) cs++; + else ce--; + } + + if(!rev) { + for (m = cs; m < ce; m++) { + xi = &(uidx->idx->a[cd].bb.a[c_str->a[m]>>32]); + x = (uint32_t)c_str->a[m]; + kv_pushp(uinfo_srt_t, *res, &p); + p->v = x; p->s = xi->ts; p->e = xi->te; + } + } else { + for (m = ce-1; m >= cs; m--) { + xi = &(uidx->idx->a[cd].bb.a[c_str->a[m]>>32]); + x = ((uint32_t)c_str->a[m])^1; + kv_pushp(uinfo_srt_t, *res, &p); + p->v = x; p->s = xi->ts; p->e = xi->te; + } + } + } else if(ru) { + x = ulg_id(uidx->uovl, cd); x <<= 1; if(rev) x^=1; + if(res->n) { + if((!ulg_type(uidx->uovl, (u->a[k-1]>>33)))) {///in previous ul is also a raw utg node + assert(get_specfic_edge(raw->g, res->a[res->n-1].v, x)); + kv_pushp(uinfo_srt_t, *res, &p); + p->v = x; p->s = 0; p->e = ru->len; + } else { + // assert(((res->a[res->n-1].v) == x)); + if(((res->a[res->n-1].v) == x)) { + res->a[res->n-1].s = 0; res->a[res->n-1].e = ru->len; + } else { + assert(get_specfic_edge(raw->g, res->a[res->n-1].v, x)); + kv_pushp(uinfo_srt_t, *res, &p); + p->v = x; p->s = 0; p->e = ru->len; + } + } + } else { + kv_pushp(uinfo_srt_t, *res, &p); + p->v = x; p->s = 0; p->e = ru->len; + } + } + } + + + for (k = 0; k < res->n; k++) { + res->a[k].n = ug_occ_w(res->a[k].s, res->a[k].e, &(raw->u.a[res->a[k].v>>1])); + } +} + +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; + 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); + + kt_for(uidx->str_b.n_thread, worker_renew_u2g_cov, uidx, i_ug->u.n); + + CALLOC(idx->cc.iug_a, i_ug->u.n); CALLOC(idx->cc.iug_idx, raw->u.n+1); + for (k = iug_occ = 0; k < i_ug->u.n; k++) { + // fprintf(stderr, "-[M::%s::] k::%lu, i_ug->u.n:%u\n", __func__, k, (uint32_t)i_ug->u.n); + gen_raw_ug_seq(uidx, uidx->pstr.str.a, &(i_ug->u.a[k]), raw, &(idx->cc.iug_a[k]), k); + iug_occ += idx->cc.iug_a[k].n; x = &(idx->cc.iug_a[k]); + for (z = 0; z < x->n; z++) idx->cc.iug_idx[x->a[z].v>>1]++; + } + + MALLOC(idx->cc.iug_b, iug_occ); + for (k = l = 0; k < raw->u.n+1; k++) { + m = idx->cc.iug_idx[k]; + idx->cc.iug_idx[k] = l; + l += m; + if(m > 0) idx->cc.iug_b[l-1] = 0; + } + assert(l == iug_occ); + + for (k = 0; k < i_ug->u.n; k++) { + x = &(idx->cc.iug_a[k]); + for (z = 0; z < x->n; z++) { + a = idx->cc.iug_b + idx->cc.iug_idx[x->a[z].v>>1]; + a_n = idx->cc.iug_idx[(x->a[z].v>>1)+1] - idx->cc.iug_idx[x->a[z].v>>1]; + if(a_n) { + if(a[a_n-1] == a_n-1) a[a_n-1] = (k<<32)|z; + else a[a[a_n-1]++] = (k<<32)|z; + } + } + } +} + +static void worker_renew_integer_bridge(void *data, long i, int tid) // callback for kt_for() +{ + ul_resolve_t *uidx = (ul_resolve_t *)data; bubble_type *bub = uidx->bub; + asg_t *bg = uidx->uovl.bg.bg; ul_str_idx_t *str_idx = &(uidx->pstr); + asg_arc_t *p = &(bg->arc[i]); uint64_t *hid_a, hid_n, z; ul_str_t *str; + uint32_t v = p->ul>>32, w = p->v, vz, wz; int64_t s, s_n, occ = 0; + + hid_a = str_idx->occ.a + str_idx->idx.a[v>>1]; + hid_n = str_idx->idx.a[(v>>1)+1] - str_idx->idx.a[v>>1]; + for (z = 0; z < hid_n; z++) { + str = &(str_idx->str.a[hid_a[z]>>32]); s_n = str->cn; + if(s_n < 2) continue; + vz = (uint32_t)(str->a[(uint32_t)hid_a[z]]); + assert((v>>1) == (vz>>1)); + s = (uint32_t)hid_a[z]; s -= 1; vz = ((uint32_t)(str->a[(uint32_t)hid_a[z]]))^1; + for (; s >= 0; s--) { + wz = (uint32_t)(str->a[s]); wz ^= 1; + if(IF_HOM((wz>>1), *bub)) continue; + if(vz == v && wz == w) occ++; + else if(vz == (v^1) && wz == (w^1)) occ++; + break; + } + + s = (uint32_t)hid_a[z]; s += 1; vz = (uint32_t)(str->a[(uint32_t)hid_a[z]]); + for (; s < s_n; s++) { + wz = (uint32_t)(str->a[s]); + if(IF_HOM((wz>>1), *bub)) continue; + if(vz == v && wz == w) occ++; + else if(vz == (v^1) && wz == (w^1)) occ++; + break; + } + } + p->ol = occ; p->ul >>= 32; p->ul <<= 32; p->ul += (((uint32_t)-1) - p->ol); +} + +void renew_u2g_bg(ul_resolve_t *uidx) +{ + ul2ul_idx_t *idx = &(uidx->uovl); ul_bg_t *bg = &(idx->bg); + 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; + + 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); + for (k = 0; k < raw->g->n_seq; k++) { + raw_id = k; if(IF_HOM(raw_id, *bub)) continue; + raw_a = idx->cc.iug_b + idx->cc.iug_idx[raw_id]; + raw_n = idx->cc.iug_idx[raw_id+1] - idx->cc.iug_idx[raw_id]; + for (z = 0, buf.n = 0; z < raw_n; z++) { + iug_id = raw_a[z]>>32; iug_off = (uint32_t)raw_a[z]; + x = &(idx->cc.iug_a[iug_id]); s_n = x->n; + assert(raw_id == (x->a[iug_off].v>>1)); + if(iug->g->seq[iug_id].del) continue; + for (s = ((int64_t)iug_off) - 1, v = x->a[iug_off].v^1; s >= 0; s--) { + if(IF_HOM((x->a[s].v>>1), *bub)) continue; + kv_push(uint64_t, buf, ((v<<32)|((uint64_t)(x->a[s].v^1)))); + break; + } + for (s = ((int64_t)iug_off) + 1, v = x->a[iug_off].v; s < s_n; s++) { + if(IF_HOM((x->a[s].v>>1), *bub)) continue; + kv_push(uint64_t, buf, ((v<<32)|((uint64_t)(x->a[s].v)))); + break; + } + } + + asg_seq_set(bg->bg, k, raw->g->seq[k].len, ((buf.n>0)?0:1)); + bg->bg->seq[k].c = 0; + + radix_sort_srt64(buf.a, buf.a + buf.n); + for (l = 0, z = 1; z <= buf.n; z++) { + if(z == buf.n || buf.a[l] != buf.a[z]) { + p = asg_arc_pushp(bg->bg); memset(p, 0, sizeof(*p)); + p->v = (uint32_t)buf.a[l]; p->ol = z - l; + p->ul = buf.a[l]>>32; p->ul <<= 32; + p->ul += (((uint32_t)-1) - p->ol); + l = z; + } + } + } + + kt_for(uidx->str_b.n_thread, worker_renew_integer_bridge, uidx, bg->bg->n_arc); + asg_cleanup(bg->bg); bg->bg->r_seq = bg->bg->n_seq; + + /***********debug***********/ + for (k = 0; k < bg->bg->n_arc; k++) { + p = &(bg->bg->arc[k]); v = p->v^1; w = (p->ul>>32)^1; + av = asg_arc_a(bg->bg, v); nv = asg_arc_n(bg->bg, v); + for (z = 0; z < nv; z++) { + if(av[z].v == w) break; + } + assert(z < nv && p->ol == av[z].ol); + } + /***********debug***********/ + + uint64_t v_occ[2], w_occ[2]; double ss = 0.500001; + for (v = 0, n_vtx = bg->bg->n_seq<<1; v < n_vtx; v++) { + if(bg->bg->seq[v>>1].del) continue; + av = asg_arc_a(bg->bg, v); nv = asg_arc_n(bg->bg, v); + if(!nv) continue; w = av[0].v; + v_occ[0] = v_occ[1] = w_occ[0] = w_occ[1] = (uint32_t)-1; + + 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))) { + if(((((v_occ[0]>>32)*ss) >= (v_occ[1]>>32)) && (((w_occ[0]>>32)*ss) >= (w_occ[1]>>32))) || + ((((v_occ[0]>>32)+(w_occ[0]>>32))*ss) >= ((v_occ[1]>>32)+(w_occ[1]>>32)))) { + for (k = 0; k < nv; k++) av[k].ou = 2;///wrong + for (k = 0; k < nw; k++) aw[k].ou = 2;///wrong + av[0].ou = aw[0].ou = 3;//correct + } else { + for (k = 0; k < nv; k++) av[k].ou = 1;///ambg + for (k = 0; k < nw; k++) aw[k].ou = 1;///ambg + } + } + } + + for (k = 0; k < bg->bg->n_arc; k++) { + p = &(bg->bg->arc[k]); v = p->v^1; w = (p->ul>>32)^1; + av = asg_arc_a(bg->bg, v); nv = asg_arc_n(bg->bg, v); + for (z = 0; z < nv; z++) { + if(av[z].v == w) break; + } + assert(z < nv && p->ol == av[z].ol); + l = MAX(p->ou, av[z].ou); if(l == 0) l = 1; + p->ou = av[z].ou = l; + } + + + MALLOC(bg->w_n, iug->u.n); MALLOC(bg->a_n, iug->u.n); + for (k = 0; k < iug->u.n; k++) { + x = &(idx->cc.iug_a[k]); bg->w_n[k] = bg->a_n[k] = 0; + for (z = nv = 0; z < x->n; z++) { + if(IF_HOM((x->a[z].v>>1), *bub)) continue; + v_occ[nv&1] = x->a[z].v; nv++; + if(nv < 2) continue; + l = get_bg_flag(uidx, v_occ[(nv-2)&1], v_occ[(nv-1)&1]); + 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]); + } + + kv_destroy(buf); +} + +/** +static void worker_update_ul_tra_idx(void *data, long i, int tid) // callback for kt_for() +{ + ul_resolve_t *uidx = (ul_resolve_t *)data; integer_t *buf = &(uidx->str_b.buf[tid]); + ma_ug_t *iug = uidx->uovl.i_ug; ul_tra_idx_t *iug_tra = &(uidx->uovl.iug_tra); + uint32_t v = i, n_tra = iug_tra_arc_n(iug_tra, v), nv, k; + if(n_tra == 0) return; + ul_tra_t *a_tra = iug_tra_arc_a(iug_tra, v); asg_arc_t *av; + + av = asg_arc_a(iug->g, v); nv = asg_arc_n(iug->g, v); kv_resize(uint64_t, buf->u, 2); + for (k = 0; k < nv; k++) { + buf->u.n = 0; buf->u.a[buf->u.n++] = v; buf->u.a[buf->u.n++] = av[k].v; + } +} + +void update_ul_tra_idx_t(ul_resolve_t *uidx) +{ + ma_ug_t *iug = uidx->uovl.i_ug; ul_tra_idx_t *iug_tra = &(uidx->uovl.iug_tra); + uint64_t v, n_vtx = iug->g->n_seq<<1, l, m; + iug_tra->arc.n = iug_tra->idx.n = 0; + kv_resize(uint32_t, iug_tra->idx, n_vtx + 1); iug_tra->idx.n = n_vtx + 1; + memset(iug_tra->idx.a, 0, iug_tra->idx.n *sizeof(*(iug_tra->idx.a))); + for (v = l = 0; v < n_vtx; v++) { + m = asg_arc_n(iug->g, v); + if(m < 2) m = 0; + iug_tra->idx.a[v] = l; + l += m; + } + iug_tra->idx.a[v] = l; + kv_resize(ul_tra_t, iug_tra->arc, l); iug_tra->arc.n = l; + + kt_for(uidx->str_b.n_thread, worker_update_ul_tra_idx, uidx, n_vtx);///all ul + ug +} +**/ + +ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt) +{ + uint64_t k, m; + ma_ug_t *ug = uidx->l1_ug; all_ul_t *uls = uidx->idx; + ul2ul_idx_t *z = &(uidx->uovl); ul_str_idx_t *str_idx = &(uidx->pstr); + z->uln = uls->n; z->gn = ug->g->n_seq; + z->tot = uls->n + ug->g->n_seq; MALLOC(z->item_idx, z->tot); + for (k = m = 0; k < z->tot; k++) { + z->item_idx[k] = (uint32_t)-1; + if(k < z->uln) {///is a ul read + if(str_idx->str.a[k].cn > 1) { + z->item_idx[k] = m; m++; + } + } else {//is a node in graph + z->item_idx[k] = m; m++; + } + } + CALLOC(z->a, m); z->n = z->m = m; + + kt_for(uidx->str_b.n_thread, worker_integer_postprecess, uidx, uls->n); + gen_integer_normalize(uidx); + chimeric_integer_deal(uidx); + append_utg_es(uidx); + kt_for(uidx->str_b.n_thread, worker_integert_debug_sym, uidx, z->tot);///all ul + ug + print_integert_ovlp_stat(z); + + remove_integert_containment(uidx); + kt_for(uidx->str_b.n_thread, worker_integert_debug_sym, uidx, z->tot);///all ul + ug + print_integert_ovlp_stat(z); + // print_uls_seq(uidx, asm_opt.output_file_name); + // print_uls_ovs(uidx, asm_opt.output_file_name); + + z->i_g = integer_sg_gen(uidx, uopt->min_ovlp); + asg_arc_del_trans(z->i_g, uopt->gap_fuzz); + + z->i_ug = ma_ug_gen(z->i_g); + CALLOC(z->i_ug->g->seq_vis, z->i_ug->g->n_seq*2); + renew_u2g_cov(uidx); + renew_u2g_bg(uidx); + + // 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); + output_integer_graph(uidx, z->i_ug, asm_opt.output_file_name); + return z; } uint64_t str_occ_w(ul_str_t *str, ul_vec_t *raw, ma_ug_t *ug) @@ -5268,9 +8791,30 @@ void gen_cul_g_t(ul_resolve_t *uidx) } } -void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg) +void init_ulg_opt_t(ulg_opt_t *z, ug_opt_t *uopt, int64_t clean_round, double min_ovlp_drop_ratio, +double max_ovlp_drop_ratio, double hom_check_drop_rate, int64_t max_tip, int64_t max_tip_hifi, bub_label_t *b_mask_t, uint32_t is_trio) { - uint64_t i; uint8_t *r_het = NULL; bubble_type *bub = NULL; + z->tipsLen = uopt->tipsLen; + z->tip_drop_ratio = uopt->tip_drop_ratio; + z->stops_threshold = uopt->stops_threshold; + z->chimeric_rate = uopt->chimeric_rate; + z->drop_ratio = uopt->drop_ratio; + + + z->b_mask_t = b_mask_t; + z->clean_round = clean_round; + z->min_ovlp_drop_ratio = min_ovlp_drop_ratio; + z->max_ovlp_drop_ratio = max_ovlp_drop_ratio; + z->hom_check_drop_rate = hom_check_drop_rate; + z->max_tip = max_tip; + z->max_tip_hifi = max_tip_hifi; + z->is_trio = is_trio; +} + +void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, +double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio) +{ + uint64_t i; uint8_t *r_het = NULL; bubble_type *bub = NULL; ulg_opt_t uu; for (i = 0; i < sg->n_seq; ++i) { if(sg->seq[i].del) continue; sg->seq[i].c = PRIMARY_LABLE; @@ -5279,13 +8823,17 @@ void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg) ma_ug_t *init_ug = ul_realignment(uopt, sg, 0); // exit(1); filter_sg_by_ug(sg, init_ug, uopt); + // print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-0"); bub = gen_bubble_chain(sg, init_ug, uopt, &r_het); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-1"); - ul_resolve_t *uidx = init_ul_resolve_t(sg, init_ug, bub, &UL_INF, r_het); + ul_resolve_t *uidx = init_ul_resolve_t(sg, init_ug, bub, &UL_INF, uopt, r_het); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-2"); // exit(1); - ul_re_correct(uidx, 3); + print_raw_uls_seq(uidx, asm_opt.output_file_name); + ul_re_correct(uidx, 3); + init_ulg_opt_t(&uu, uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_tip<<1, b_mask_t, is_trio); + /**ul2ul_idx_t *u2o = **/gen_ul2ul(uidx, uopt, &uu); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-3"); // print_debug_ul("UL.debug", init_ug, sg, uopt->coverage_cut, uopt->sources, uopt->ruIndex, bub, &UL_INF); diff --git a/gfa_ut.h b/gfa_ut.h index 8382b9c..5c3b6c6 100644 --- a/gfa_ut.h +++ b/gfa_ut.h @@ -14,5 +14,6 @@ void asg_arc_cut_bub_links(asg_t *g, asg64_v *in, float len_rat, float sec_len_r void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float ou_rat, uint32_t is_ou, bub_label_t *b_mask_t); uint32_t asg_cut_large_indel(asg_t *g, asg64_v *in, int32_t max_ext, float ou_rat, uint32_t is_ou); uint32_t asg_cut_semi_circ(asg_t *g, uint32_t lim_len, uint32_t is_clean); -void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg); +void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, +double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio); #endif diff --git a/inter.cpp b/inter.cpp index 4177c87..7c3242c 100644 --- a/inter.cpp +++ b/inter.cpp @@ -8913,18 +8913,9 @@ static void filter_short_ulalignments(void *data, long i, int tid) // callback f const ma_ug_t *ug = (ma_ug_t *)data; uc_block_t *a = NULL; uc_block_t *p; int64_t k, a_n; uint32_t z, fz, lz, l, bz; a = UL_INF.a[i].bb.a; a_n = UL_INF.a[i].bb.n; - // if(i == 4) { - // fprintf(stderr, "[M::%s::i->%ld] a_n->%ld\n", __func__, i, a_n); - // } + for (k = a_n - 1; k >= 0; k--) { p = &(a[k]); - // if(i == 27512) { - // fprintf(stderr, "+[M::%s::k->%ld::a_n->%ld] p->ts::%u, p->te::%u, p->pchain::%u, p->pidx::%u, p->aidx::%u, p->qs::%u, p->qe::%u\n", - // __func__, k, a_n, p->ts, p->te, p->pchain, p->pidx, p->aidx, p->qs, p->qe); - // } - // if(p->qs > p->qe) { - // fprintf(stderr, "[M::%s] id->%ld, qs->%u, qe->%u\n", __func__, i, p->qs, p->qe); - // } if(p->base || (!p->el) || (!p->pchain)) continue; if(p->pidx == (uint32_t)-1) { if(!ugl_cover_check(p->ts, p->te, &(ug->u.a[p->hid]))) { @@ -8956,10 +8947,6 @@ static void filter_short_ulalignments(void *data, long i, int tid) // callback f for (k = a_n - 1; k >= 0; k--) { p = &(a[k]); - // if(i == 27512) { - // fprintf(stderr, "-[M::%s::k->%ld::a_n->%ld] p->ts::%u, p->te::%u, p->pchain::%u, p->pidx::%u, p->aidx::%u, p->qs::%u, p->qe::%u\n", - // __func__, k, a_n, p->ts, p->te, p->pchain, p->pidx, p->aidx, p->qs, p->qe); - // } if(p->base || (!p->el) || (!p->pchain)) continue; if(p->pidx != (uint32_t)-1) { assert(a[p->pidx].aidx == (uint32_t)k);