From 475ebb8075ef1b3f7c88eea053871e89a2b1f766 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Tue, 20 Apr 2021 01:29:15 -0400 Subject: [PATCH] keep long range information --- Overlaps.cpp | 495 +++++++++++++-------------------------------------- Overlaps.h | 10 +- hic.cpp | 381 ++++++++++++++++++++++++++++++++++----- hic.h | 5 +- 4 files changed, 465 insertions(+), 426 deletions(-) diff --git a/Overlaps.cpp b/Overlaps.cpp index 542a888..86220fb 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -11917,6 +11917,121 @@ ma_ug_t *ug, asg_t *read_sg, hap_cov_t *cov) set_utg_offset(aux->b.a, aux->b.n, ug, read_sg, cov, 1, 0); } +// void collect_trans_cov(const char* cmd, buf_t* pri, uint64_t pri_offset, buf_t* aux, uint64_t aux_offset, +// ma_ug_t *ug, asg_t *read_sg, hap_cov_t *cov) +// { +// uint32_t i, k, rid, occ, ori, thre_pri; +// uint32_t p_uId, c_uId, x_occ, y_occ; +// uint64_t len_aux, uLen, uCov; +// ma_utg_t* u = NULL; +// trans_chain* t_ch = cov->t_ch; +// if(pri->b.n == 0 || aux->b.n == 0) return; + +// len_aux = set_utg_offset(aux->b.a, aux->b.n, ug, read_sg, cov, 0, 0); +// chain_trans_ovlp(cov, ug, read_sg, pri, len_aux, &thre_pri); +// if(thre_pri > 0) +// { +// /*******************************for debug************************************/ +// // fprintf(stderr, "\n%s, thre_pri: %u, len_aux: %lu\n", cmd, thre_pri, len_aux); +// // print_buf_t(ug, pri, "pri"); +// // print_buf_t(ug, aux, "aux"); +// /*******************************for debug************************************/ + +// if(t_ch) +// { +// chain_origin_trans_uid_by_distance(cov, read_sg, pri->b.a, pri->b.n, pri_offset, NULL, +// aux->b.a, aux->b.n, aux_offset, &len_aux, ug, RC_1, -1024, cmd); +// } + + +// for (i = uCov = 0, p_uId = (uint32_t)-1; i < aux->b.n; i++) +// { +// u = &(ug->u.a[aux->b.a[i]>>1]); +// if(u->n == 0) continue; +// ori = aux->b.a[i] & 1; +// for (k = 0; k < u->n; k++) +// { +// rid = (ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33)); +// uCov += cov->cov[rid]; + +// // if(t_ch) t_ch->is_r_het[(ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33))] |= P_HET; +// if(t_ch) +// { +// t_ch->is_r_het[(ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33))] |= P_HET; +// c_uId = get_origin_uid((ori == 1?((u->a[u->n-k-1]^(uint64_t)(0x100000000))>>32):(u->a[k]>>32)), +// t_ch, NULL, NULL); +// if(c_uId == (uint32_t)-1 || p_uId == c_uId) continue; +// p_uId = c_uId; +// kv_push(uint32_t, t_ch->st.uIDs, c_uId); +// } +// } +// } + +// if(t_ch) kv_push(uint32_t, t_ch->st.iDXs, t_ch->st.uIDs.n);///dedup_push_trans_chain(t_ch); + +// for (i = uLen = occ = 0; i < pri->b.n; i++) +// { +// u = &(ug->u.a[pri->b.a[i]>>1]); +// if(u->n == 0) continue; +// ori = pri->b.a[i] & 1; +// for (k = 0; k < u->n; k++, occ++) +// { +// if(occ >= thre_pri) break; +// rid = (ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33)); +// uLen += read_sg->seq[rid].len; + +// //if(t_ch) t_ch->is_r_het[(ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33))] |= P_HET; +// if(t_ch) +// { +// t_ch->is_r_het[(ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33))] |= P_HET; +// c_uId = get_origin_uid((ori == 1?((u->a[u->n-k-1]^(uint64_t)(0x100000000))>>32):(u->a[k]>>32)), +// t_ch, NULL, NULL); +// if(c_uId == (uint32_t)-1 || p_uId == c_uId) continue; +// p_uId = c_uId; +// kv_push(uint32_t, t_ch->st.uIDs, c_uId); +// } +// } +// if(occ >= thre_pri) break; +// } + +// if(t_ch) kv_push(uint32_t, t_ch->st.iDXs, t_ch->st.uIDs.n);///dedup_push_trans_chain(t_ch); + +// if(t_ch) +// { +// x_occ = y_occ = 0; +// get_chain_trans(t_ch, t_ch->chain_num, NULL, &x_occ, NULL, &y_occ); +// if(x_occ == 0 || y_occ == 0) +// { +// t_ch->uIDs.n -= (x_occ + y_occ); +// t_ch->iDXs.n -= 2; +// } +// else +// { +// t_ch->chain_num++; +// t_ch->l0_chain++; +// } +// } + +// uCov = (uLen == 0? 0 : uCov / uLen); + +// for (i = occ = 0; i < pri->b.n; i++) +// { +// u = &(ug->u.a[pri->b.a[i]>>1]); +// if(u->n == 0) continue; +// ori = pri->b.a[i] & 1; +// for (k = 0; k < u->n; k++, occ++) +// { +// if(occ >= thre_pri) break; +// ///rid = u->a[k]>>33; +// rid = (ori == 1?(u->a[u->n-k-1]>>33):(u->a[k]>>33)); +// cov->cov[rid] += (uCov * read_sg->seq[rid].len); +// } +// if(occ >= thre_pri) break; +// } +// } +// set_utg_offset(aux->b.a, aux->b.n, ug, read_sg, cov, 1, 0); +// } + int asg_arc_cut_long_equal_tips_assembly_complex(asg_t *g, ma_hit_t_alloc* reverse_sources, long long miniedgeLen, uint32_t stops_threshold, R_to_U* ruIndex) @@ -12295,6 +12410,10 @@ trans_chain* init_trans_chain(ma_ug_t *ug, uint64_t r_num) kv_malloc(x->bed, x->u_num); x->bed.n = x->u_num; for (k = 0; k < x->bed.n; k++) kv_init(x->bed.a[k]); + x->st.chain_num = 0; + kv_init(x->st.uIDs); + kv_init(x->st.iDXs); + kv_push(uint32_t, x->st.iDXs, 0); return x; } @@ -12315,6 +12434,8 @@ void destory_trans_chain(trans_chain **x) kv_destroy((*x)->topo_buf); kv_destroy((*x)->topo_res); kv_destroy((*x)->c_buf); + kv_destroy((*x)->st.uIDs); + kv_destroy((*x)->st.iDXs); ///free((*x)->uLen); free((*x)); } @@ -26452,156 +26573,6 @@ bub_label_t* b_mask_t) set_rtg_flag_by_bubble(bub, ug, r_g, m_bub, expect_vis, 0); } } - /** - for (i = beg_idx; i < beg_idx + occ; i++) - { - get_bubbles(bub, i, &beg_utg, &sink_utg, &a, &n, NULL); - if(beg_utg == (uint32_t)-1 || sink_utg == (uint32_t)-1) continue; - - l_bub = r_bub = (uint64_t)-1; m_bub = i; - get_bub_id(bub, beg_utg>>1, &bub_0, &bub_1, 0); - if(bub_0 != i) l_bub = bub_0; - if(bub_1 != i) l_bub = bub_1; - - get_bub_id(bub, sink_utg>>1, &bub_0, &bub_1, 0); - if(bub_0 != i) r_bub = bub_0; - if(bub_1 != i) r_bub = bub_1; - - - set_rtg_flag_by_bubble(bub, ug, r_g, l_bub, expect_vis, 1); - set_rtg_flag_by_bubble(bub, ug, r_g, r_bub, expect_vis, 1); - set_rtg_flag_by_bubble(bub, ug, r_g, m_bub, expect_vis, 1); - for (k_i = 0; k_i < n; k_i++) - { - uId = a[k_i]>>1; - nsu = &(ug->u.a[uId]); - if(nsu->m == 0) continue; - for (k_v = 0; k_v < 2; k_v++) - { - v = (uId<<1) + k_v; - if(get_real_length(nsg, v, NULL) != 0) continue; - ///fprintf(stderr, "++++++tig-utg%.6ul\n", uId+1); - - if(v&1) - { - endRid = nsu->start^1; - } - else - { - endRid = nsu->end^1; - } - - //x is the end read of a tip - ///find all overlap of x - x = &(sources[(endRid>>1)]); - ava_ol_max = ava_max = 0; ava_min_chain = (uint32_t)-1; - h_max = NULL; - for (k = 0; k < x->length; k++) - { - ///fprintf(stderr, "k: %u\n", k); - ///h is the edge of endRid - h = &(x->buffer[k]); - ///means we found a contained read - if(get_contained_reads_chain_by_broken_bub(h, sources, coverage_cut, ruIndex, ug, r_g, - max_hang, min_ovlp, endRid, uId, &chain_edges, &ava_cur, &test_oLen, &ava_chainLen, - expect_vis, circle_vis, utg_vis, chainLenThres, 1)) - { - is_update = 0; - - if(ava_cur > ava_max) - { - is_update = 1; - } - else if(ava_cur == ava_max) - { - if(ava_chainLen < ava_min_chain) - { - is_update = 1; - } - else if(ava_chainLen == ava_min_chain && test_oLen > ava_ol_max) - { - is_update = 1; - } - } - - if(is_update) - { - ava_min_chain = ava_chainLen; - ava_max = ava_cur; - ava_ol_max = test_oLen; - h_max = h; - } - } - } - - if(ava_max > 0) - { - ///fprintf(stderr, "ava_max: %u\n", ava_max); - get_contained_reads_chain_by_broken_bub(h_max, sources, coverage_cut, ruIndex, ug, r_g, - max_hang, min_ovlp, endRid, uId, &chain_edges, &ava_cur, &test_oLen, &ava_chainLen, - expect_vis, circle_vis, NULL, chainLenThres, 1); - if(chain_edges.a.n < 1) continue; - ///the last cantained read - t_max = chain_edges.a.a[chain_edges.a.n-1]; - - k = 0; rbub_edges.n = 0; - ///edges from the last contained read to other unitigs - - while(get_edge2existing_node_advance_by_broken_bub(ug, r_g, sources, coverage_cut, - ruIndex, expect_vis, max_hang, min_ovlp, t_max.v, &k, &r_edge, 1)) - { - kv_push(asg_arc_t, rbub_edges, r_edge); - } - - ///need to do transitive reduction - ///note here is different to standard transitive reduction - minor_transitive_reduction_r_g(r_g, rbub_edges.a, rbub_edges.n); - - for (k = 0; k < chain_edges.a.n; k++) - { - t_max = chain_edges.a.a[k]; - ///save all infor for reverting - get_R_to_U(ruIndex, t_max.v>>1, &contain_rId, &is_Unitig); - a_nodes=contain_rId; - a_nodes=a_nodes<<32; - a_nodes=a_nodes|((uint64_t)(t_max.v>>1)); - kv_push(uint64_t, u_vecs.a, a_nodes); - - r_g->seq[t_max.v>>1].del = 0; - coverage_cut[t_max.v>>1].del = 0; - coverage_cut[t_max.v>>1].c = PRIMARY_LABLE; - - get_edge_from_source(sources, coverage_cut, ruIndex, max_hang, min_ovlp, - (t_max.ul>>32), t_max.v, &t); - kv_push(asg_arc_t, new_edges, t); - - get_edge_from_source(sources, coverage_cut, ruIndex, max_hang, min_ovlp, - (t_max.v^1), ((t_max.ul>>32)^1), &t); - kv_push(asg_arc_t, new_edges, t); - } - - - for (k = 0; k < rbub_edges.n; k++) - { - t = rbub_edges.a[k]; - if(t.del) continue; - - kv_push(asg_arc_t, new_edges, t); - get_edge_from_source(sources, coverage_cut, ruIndex, max_hang, min_ovlp, - (t.v^1), ((t.ul>>32)^1), &t); - kv_push(asg_arc_t, new_edges, t); - } - - } - } - } - - - set_rtg_flag_by_bubble(bub, ug, r_g, l_bub, expect_vis, 0); - set_rtg_flag_by_bubble(bub, ug, r_g, r_bub, expect_vis, 0); - set_rtg_flag_by_bubble(bub, ug, r_g, m_bub, expect_vis, 0); - } - **/ asg_arc_t* p = NULL; for (k = 0; k < new_edges.n; k++) @@ -26779,80 +26750,6 @@ bub_label_t* b_mask_t) } } - /** - for (i = beg_idx; i < beg_idx + occ; i++) - { - get_bubbles(bub, i, &beg_utg, &sink_utg, &a, &n, NULL); - ///fprintf(stderr, "\n%u-th: bub-utg%.6ul, bub-utg%.6ul\n", i - beg_idx, (beg_utg>>1)+1, (sink_utg>>1)+1); - if(beg_utg == (uint32_t)-1 || sink_utg == (uint32_t)-1) continue; - - l_bub = r_bub = (uint64_t)-1; m_bub = i; - get_bub_id(bub, beg_utg>>1, &bub_0, &bub_1, 0); - if(bub_0 != i) l_bub = bub_0; - if(bub_1 != i) l_bub = bub_1; - - get_bub_id(bub, sink_utg>>1, &bub_0, &bub_1, 0); - if(bub_0 != i) r_bub = bub_0; - if(bub_1 != i) r_bub = bub_1; - - - set_rtg_flag_by_bubble(bub, ug, r_g, l_bub, expect_vis, 1); - set_rtg_flag_by_bubble(bub, ug, r_g, r_bub, expect_vis, 1); - set_rtg_flag_by_bubble(bub, ug, r_g, m_bub, expect_vis, 1); - for (k_i = 0; k_i < n; k_i++) - { - uId = a[k_i]>>1; - nsu = &(ug->u.a[uId]); - if(nsu->m == 0) continue; - for (k_v = 0; k_v < 2; k_v++) - { - v = (uId<<1) + k_v; - if(get_real_length(nsg, v, NULL) != 0) continue; - ///fprintf(stderr, "++++++tig-utg%.6ul\n", uId+1); - - if(v&1) - { - endRid = nsu->start^1; - } - else - { - endRid = nsu->end^1; - } - - - k = 0; rbub_edges.n = 0; - while(get_edge2existing_node_advance_by_broken_bub(ug, r_g, sources, coverage_cut, - ruIndex, expect_vis, max_hang, min_ovlp, endRid, &k, &r_edge, 1)) - { - kv_push(asg_arc_t, rbub_edges, r_edge); - } - - if(rbub_edges.n > 0) - { - ///need to do transitive reduction - ///note here is different to standard transitive reduction - minor_transitive_reduction_r_g(r_g, rbub_edges.a, rbub_edges.n); - for (k = 0; k < rbub_edges.n; k++) - { - t = rbub_edges.a[k]; - if(t.del) continue; - - kv_push(asg_arc_t, new_edges, t); - get_edge_from_source(sources, coverage_cut, ruIndex, max_hang, min_ovlp, - (t.v^1), ((t.ul>>32)^1), &t); - kv_push(asg_arc_t, new_edges, t); - } - } - } - } - - - set_rtg_flag_by_bubble(bub, ug, r_g, l_bub, expect_vis, 0); - set_rtg_flag_by_bubble(bub, ug, r_g, r_bub, expect_vis, 0); - set_rtg_flag_by_bubble(bub, ug, r_g, m_bub, expect_vis, 0); - } - **/ - asg_arc_t* p = NULL; for (k = 0; k < new_edges.n; k++) { @@ -27090,154 +26987,6 @@ R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t backward_steps, uint32_t b } } - /** - for (i = beg_idx; i < beg_idx + occ; i++) - { - get_bubbles(bub, i, &beg_utg, &sink_utg, &a, &n, NULL); - if(beg_utg == (uint32_t)-1 || sink_utg == (uint32_t)-1) continue; - - l_bub = r_bub = (uint64_t)-1; m_bub = i; - get_bub_id(bub, beg_utg>>1, &bub_0, &bub_1, 0); - if(bub_0 != i) l_bub = bub_0; - if(bub_1 != i) l_bub = bub_1; - - get_bub_id(bub, sink_utg>>1, &bub_0, &bub_1, 0); - if(bub_0 != i) r_bub = bub_0; - if(bub_1 != i) r_bub = bub_1; - - - set_rtg_flag_by_bubble(bub, ug, r_g, l_bub, expect_vis, 1); - set_rtg_flag_by_bubble(bub, ug, r_g, r_bub, expect_vis, 1); - set_rtg_flag_by_bubble(bub, ug, r_g, m_bub, expect_vis, 1); - for (k_i = 0; k_i < n; k_i++) - { - uId = a[k_i]>>1; - nsu = &(ug->u.a[uId]); - if(nsu->m == 0) continue; - - rbub_edges.n = round = 0; - for (k_v = 0; k_v < 2; k_v++) - { - if(rbub_edges.n > 0) - { - cur_backward_steps = nsu->n - round - 1; - if(cur_backward_steps > backward_steps) - { - cur_backward_steps = backward_steps; - } - } - else - { - cur_backward_steps = backward_steps; - } - - - v = (uId<<1) + k_v; - if(get_real_length(nsg, v, NULL) != 0) continue; - ///fprintf(stderr, "++++++tig-utg%.6ul\n", uId+1); - ///that means this unitig has been changed - // if(nsu->start!=((uint64_t)(nsu->a[0])>>32)) continue; - // if((nsu->end^1)!=((uint64_t)(nsu->a[nsu->n-1])>>32)) continue; - - if(v&1) - { - init = 0; - step = 1; - mode = 1; - } - else - { - init = nsu->n - 1; - step = -1; - mode = 0; - } - - rbub_edges.n = 0; - for (round = 0; round < cur_backward_steps && init >= 0 && init < (long long)nsu->n; - init = init + step, round++) - { - endRid = ((uint64_t)(nsu->a[init]))>>32; - endRid = endRid^mode; - - k = 0; rbub_edges.n = 0; - while(get_edge2existing_node_advance_by_broken_bub(ug, r_g, sources, coverage_cut, - ruIndex, expect_vis, max_hang, min_ovlp, endRid, &k, &r_edge, 1)) - { - kv_push(asg_arc_t, rbub_edges, r_edge); - } - - if(rbub_edges.n > 0) break; - } - - if(rbub_edges.n > 0) - { - //save for revert - tmp = mode; tmp = tmp <<31; tmp = tmp | (uint64_t)(init); tmp = tmp << 32; tmp = tmp | uId; - kv_push(uint64_t, u_vecs.a, tmp); - ///need to do transitive reduction - ///note here is different to standard transitive reduction - minor_transitive_reduction_r_g(r_g, rbub_edges.a, rbub_edges.n); - - - ///modify read graph - for (init = init - step; init >= 0 && init < (long long)nsu->n; init = init - step) - { - w = ((uint64_t)(nsu->a[init]))>>32; - nv = asg_arc_n(r_g, w); - av = asg_arc_a(r_g, w); - for (k = 0; k < nv; k++) - { - if(av[k].del) continue; - kv_push(asg_arc_t, recover_edges.a, av[k]); - if(asg_get_arc(r_g, av[k].v^1, av[k].ul>>32^1, &t)==0) - { - fprintf(stderr, "error\n"); - } - kv_push(asg_arc_t, recover_edges.a, t); - } - - - nv = asg_arc_n(r_g, w^1); - av = asg_arc_a(r_g, w^1); - for (k = 0; k < nv; k++) - { - if(av[k].del) continue; - kv_push(asg_arc_t, recover_edges.a, av[k]); - if(asg_get_arc(r_g, av[k].v^1, av[k].ul>>32^1, &t)==0) - { - fprintf(stderr, "error\n"); - } - kv_push(asg_arc_t, recover_edges.a, t); - } - - ///w = ((uint64_t)(nsu->a[init]))>>32; - asg_seq_del(r_g, w>>1); - expect_vis[w>>1] = 0; - } - - - for (k = 0; k < rbub_edges.n; k++) - { - t = rbub_edges.a[k]; - if(t.del) continue; - - kv_push(asg_arc_t, new_edges, t); - get_edge_from_source(sources, coverage_cut, ruIndex, max_hang, min_ovlp, - (t.v^1), ((t.ul>>32)^1), &t); - kv_push(asg_arc_t, new_edges, t); - } - } - - } - } - - - set_rtg_flag_by_bubble(bub, ug, r_g, l_bub, expect_vis, 0); - set_rtg_flag_by_bubble(bub, ug, r_g, r_bub, expect_vis, 0); - set_rtg_flag_by_bubble(bub, ug, r_g, m_bub, expect_vis, 0); - } - **/ - asg_arc_t* p = NULL; for (k = 0; k < new_edges.n; k++) { @@ -27346,7 +27095,7 @@ void reset_bub(bubble_type* bub, ma_ug_t *ug, trans_chain* back_ug_chain, kvec_a new_rtg_edges->a.n = 0; ///classify_untigs(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, max_hang, min_ovlp); - identify_bubbles(ug, bub, back_ug_chain->is_r_het); + identify_bubbles(ug, bub, back_ug_chain->is_r_het, NULL); update_bubble_chain(ug, bub, 0, 1); resolve_bubble_chain_tangle(ug, bub); // fprintf(stderr, "bub.f_bub: %lu, bub.b_bub: %lu, bub.b_end_bub: %lu, bub.tangle_bub: %lu, bub.cross_bub: %lu\n", diff --git a/Overlaps.h b/Overlaps.h index 6caf5c2..33d20c9 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1060,8 +1060,9 @@ void asg_bub_backtrack_primary(asg_t *g, uint32_t v0, buf_t *b); typedef struct{ double weight; - uint32_t uID:31, del:1; + uint32_t uID; uint64_t dis; + uint8_t is_cc:7, del:1; uint64_t occ; ///uint64_t occ:63, scaff:1; ///uint32_t enzyme; @@ -1102,6 +1103,12 @@ typedef struct { ca_buf_t* a; } kv_ca_buf_t; +typedef struct { + kvec_t(uint32_t) uIDs; + kvec_t(uint32_t) iDXs; + uint32_t chain_num; +} sub_tran_t; + typedef struct{ uint32_t* rUidx; uint64_t* rUpos; @@ -1115,6 +1122,7 @@ typedef struct{ kv_u_trans_t k_trans; kv_u_trans_hit_t k_t_b; kv_ca_buf_t c_buf; + sub_tran_t st; }trans_chain; typedef struct { diff --git a/hic.cpp b/hic.cpp index bebdb69..5d2a1b0 100644 --- a/hic.cpp +++ b/hic.cpp @@ -291,7 +291,7 @@ typedef struct{ typedef struct{ trans_p_t* a; size_t n, m; - uint64_t max; + uint64_t max, med; } trans_idx; @@ -2264,8 +2264,11 @@ void dfs_bubble(asg_t *g, kvec_t_u32_warp* stack, kvec_t_u32_warp* result, uint3 } } -uint32_t get_unitig_het_arb(ma_utg_t* u, uint8_t *r_het_flag, uint32_t m_het_label, uint32_t p_het_label, uint32_t n_het_label) +uint32_t get_unitig_het_arb(ma_ug_t* ug, uint32_t uid, uint8_t *r_het_flag, kv_u_trans_t *ref, uint32_t m_het_occ, +uint32_t m_het_label, uint32_t p_het_label, uint32_t n_het_label) { + if(u_trans_n(*ref, uid) > 0) return m_het_label; + ma_utg_t *u = &(ug->u.a[uid]); uint32_t k, rId; uint32_t het_occ, hom_occ; for (k = 0, het_occ = hom_occ = 0; k < u->n; k++) @@ -2281,14 +2284,18 @@ uint32_t get_unitig_het_arb(ma_utg_t* u, uint8_t *r_het_flag, uint32_t m_het_lab hom_occ++; } } + // if(uid == 18759) + // { + // fprintf(stderr, "het_occ: %u, hom_occ: %u\n", het_occ, hom_occ); + // } if((het_occ+hom_occ) == 0) return n_het_label; ///hom - if(het_occ > ((het_occ+hom_occ)*0.8)) return m_het_label; ///must het + if((het_occ > ((het_occ+hom_occ)*0.8)) && ((het_occ+hom_occ) > m_het_occ)) return m_het_label; ///must het if(het_occ >= hom_occ) return p_het_label; ///potential het return n_het_label; ///hom } void update_bub_b_s_idx(bubble_type* bub); -void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag) +void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_trans_t *ref) { asg_cleanup(ug->g); if (!ug->g->is_symm) asg_symm(ug->g); @@ -2391,7 +2398,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag) for (i = 0; i < ug->g->n_seq; i++) { - bub->index[i] = get_unitig_het_arb(&(ug->u.a[i]), r_het_flag, M_het(*bub), P_het(*bub), (uint32_t)-1); + bub->index[i] = get_unitig_het_arb(ug, i, r_het_flag, ref, 20, M_het(*bub), P_het(*bub), (uint32_t)-1); } for (i = 0; i < bub->f_bub; i++) @@ -2460,8 +2467,8 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag) } bub->b_g = NULL; bub->b_ug = NULL; - build_bub_graph(ug, bub); + // fprintf(stderr, "-bub->index[18759]: %u, bub->num.n: %u\n", (uint32_t)bub->index[18759], bub->num.n); } void print_bubbles(ma_ug_t* ug, bubble_type* bub, kvec_pe_hit* hits, hc_links* link, ha_ug_index* idx) @@ -2635,7 +2642,7 @@ void print_bubbles(ma_ug_t* ug, bubble_type* bub, kvec_pe_hit* hits, hc_links* l // fprintf(stderr, "************het utgs************\n"); } -void push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t* d) +hc_edge* push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t* d) { uint64_t k, n; hc_edge* a = NULL; @@ -2658,7 +2665,7 @@ void push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t { a[k].weight += weight; if(d) a[k].dis = (*d); - return; + return &(a[k]); } } @@ -2676,6 +2683,7 @@ void push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t p->uID = uID; p->weight = weight; if(d) p->dis = (*d); + return p; } long long get_enzyme_occ_debug(char* t, long long tlen, char* p, long long plen) @@ -3124,6 +3132,7 @@ static void worker_for_dis(void *data, long i, int tid) for (k = 0; k < t->e.n; k++) { if(t->e.a[k].del) continue; + if(i == t->e.a[k].uID) continue; u = t->e.a[k].uID; for (v = ((uint64_t)(i)<<1); v < ((uint64_t)(i+1)<<1); v++)///two directions @@ -3794,6 +3803,7 @@ void update_containment_distance(asg_t *sg, kv_u_trans_t *ta, hc_links* link) kv_destroy(buf.a); } + void collect_hc_links(const ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, MT* M) { double index_time = yak_realtime(); @@ -3840,6 +3850,249 @@ void collect_hc_links(const ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, fprintf(stderr, "[M::%s::%.3f] ==> Enzymes have been counted\n", __func__, yak_realtime()-index_time); } + +void update_dis_connected_gfa(asg_t *sg, hc_links *link, MT *M) +{ + uint32_t i, m, v, nv, x, y; + asg_arc_t *av = NULL; + kvec_t(uint32_t) stack; kv_init(stack); + uint32_t *flag = NULL; MALLOC(flag, sg->n_seq); + uint64_t *group = NULL; MALLOC(group, sg->n_seq); + + for (i = 0; i < link->a.n; i++) + { + for (m = 0; m < link->a.a[i].e.n; m++) + { + if(link->a.a[i].e.a[m].dis == (uint64_t)-1) + { + if(link->a.a[i].e.a[m].is_cc != 0) + { + fprintf(stderr, "ERROR 1\n"); + } + } + + if(link->a.a[i].e.a[m].occ == (uint64_t)-1) + { + if(link->a.a[i].e.a[m].is_cc != 0 || link->a.a[i].e.a[m].dis == (uint64_t)-1) + { + fprintf(stderr, "ERROR 2\n"); + } + } + } + } + + + memset(flag, -1, sg->n_seq*sizeof(uint32_t)); + // connected componets + for (i = 0; i < sg->n_seq; ++i) { + if (flag[i] != (uint32_t)-1) continue; + stack.n = 0; + kv_push(uint32_t, stack, i); + while (stack.n > 0) { + stack.n--; + flag[stack.a[stack.n]] = i;///group id + + v = (stack.a[stack.n])<<1; + av = asg_arc_a(sg, v); + nv = asg_arc_n(sg, v); + for (m = 0; m < nv; ++m) { + if (flag[av[m].v>>1] != (uint32_t)-1) continue; + kv_push(uint32_t, stack, av[m].v>>1); + } + + v++; + av = asg_arc_a(sg, v); + nv = asg_arc_n(sg, v); + for (m = 0; m < nv; ++m) { + if (flag[av[m].v>>1] != (uint32_t)-1) continue; + kv_push(uint32_t, stack, av[m].v>>1); + } + } + } + kv_destroy(stack); + + + + // precalculate the size of each group + for (i = 0; i < sg->n_seq; ++i) + group[i] = (uint64_t)flag[i] << 32 | i; + radix_sort_hc64(group, group + sg->n_seq); + for (i = 1, x = y = 0; i <= sg->n_seq; ++i) { + if (i == sg->n_seq || (group[i]>>32) != (group[x]>>32)) { + uint32_t j; + for (j = x; j < i; ++j) + group[j] = (uint64_t)y << 32 | (uint32_t)group[j];///(group id)|first element in this group + ++y, x = i; + } + } + + memset(flag, 0, sg->n_seq*sizeof(uint32_t)); + for (i = 1, x = y = 0; i <= sg->n_seq; ++i) + { + if (i == sg->n_seq || (group[i]>>32) != (group[x]>>32)) + { + x = i; + } + } + + + free(flag); +} + +void append_tig_link(uint64_t *cc, uint32_t cc_off, uint32_t cc_size, hc_links *link, asg_t *sg, kvec_t_u32_warp *buf) +{ + uint32_t i, k, id, tig_occ, b_cc_occ, *tig, *b_cc, qn, tn; + uint64_t t_d = (uint64_t)-1; + hc_edge *p = NULL; + for (i = 0; i < cc_size; ++i) ///how many nodes + { + id = (uint32_t)cc[cc_off + i];///node id + for (k = 0; k < link->a.a[id].e.n; k++) + { + if(link->a.a[id].e.a[k].is_cc == 0) break; + } + if(k < link->a.a[id].e.n) break; + } + if(i >= cc_size) return; + + buf->a.n = 0; tig_occ = 0; + for (i = 0; i < cc_size; ++i) ///how many nodes + { + id = (uint32_t)cc[cc_off + i];///node id + + if(asg_arc_n(sg, (id<<1)) == 0 || asg_arc_n(sg, ((id<<1)+1)) == 0)///tig + { + kv_push(uint32_t, buf->a, id); + tig_occ++; + } + } + + if(tig_occ == 0) return; + + b_cc_occ = 0; + for (i = 0; i < cc_size; ++i) ///how many nodes + { + id = (uint32_t)cc[cc_off + i];///node id + if(link->a.a[id].e.n == 0) continue; + + for (k = 0; k < link->a.a[id].e.n; k++) + { + if(link->a.a[id].e.a[k].is_cc == 0) break; + } + if(k < link->a.a[id].e.n)///disconnected + { + kv_push(uint32_t, buf->a, id); + b_cc_occ++; + } + } + + if(b_cc_occ == 0) return; + + tig = buf->a.a; + b_cc = buf->a.a + tig_occ; + for (i = 0; i < tig_occ; i++) + { + qn = tig[i]; + for (k = 0; k < b_cc_occ; k++) + { + tn = b_cc[k]; + if(qn == tn) continue; + p = get_hc_edge(link, qn, tn, 0); + if(p) continue; + + t_d = (uint64_t)-1; + p = push_hc_edge(&(link->a.a[qn]), tn, 0, 0, &t_d); + p->is_cc = 0; p->occ = (uint64_t)-1; + p = push_hc_edge(&(link->a.a[tn]), qn, 0, 0, &t_d); + p->is_cc = 0; p->occ = (uint64_t)-1; + } + } +} + +void update_ug_by_tigs(asg_t *sg, hc_links *link) +{ + uint32_t i, m, v, nv, x, y, qn, tn; + asg_arc_t *av = NULL; + kvec_t(uint32_t) stack; kv_init(stack); + uint32_t *flag = NULL; MALLOC(flag, sg->n_seq); + uint64_t *group = NULL; MALLOC(group, sg->n_seq); + kvec_t_u32_warp buf; kv_init(buf.a); + + + memset(flag, -1, sg->n_seq*sizeof(uint32_t)); + // connected componets + for (i = 0; i < sg->n_seq; ++i) { + if (flag[i] != (uint32_t)-1) continue; + stack.n = 0; + kv_push(uint32_t, stack, i); + while (stack.n > 0) { + stack.n--; + flag[stack.a[stack.n]] = i;///group id + + v = (stack.a[stack.n])<<1; + av = asg_arc_a(sg, v); + nv = asg_arc_n(sg, v); + for (m = 0; m < nv; ++m) { + if (flag[av[m].v>>1] != (uint32_t)-1) continue; + kv_push(uint32_t, stack, av[m].v>>1); + } + + v++; + av = asg_arc_a(sg, v); + nv = asg_arc_n(sg, v); + for (m = 0; m < nv; ++m) { + if (flag[av[m].v>>1] != (uint32_t)-1) continue; + kv_push(uint32_t, stack, av[m].v>>1); + } + } + } + kv_destroy(stack); + + + + // precalculate the size of each group + for (i = 0; i < sg->n_seq; ++i) + group[i] = (uint64_t)flag[i] << 32 | i; + radix_sort_hc64(group, group + sg->n_seq); + for (i = 1, x = y = 0; i <= sg->n_seq; ++i) { + if (i == sg->n_seq || (group[i]>>32) != (group[x]>>32)) { + uint32_t j; + for (j = x; j < i; ++j) + group[j] = (uint64_t)y << 32 | (uint32_t)group[j];///(group id)|first element in this group + ++y, x = i; + } + } + + + for (i = 0; i < link->a.n; i++) + { + qn = i; + for (m = 0; m < link->a.a[i].e.n; m++) + { + link->a.a[i].e.a[m].occ = 0; + link->a.a[i].e.a[m].is_cc = 0; + tn = link->a.a[i].e.a[m].uID; + if((group[qn]>>32) == (group[tn]>>32)) + { + link->a.a[i].e.a[m].is_cc = 1; + } + } + } + + + for (i = 1, x = 0; i <= sg->n_seq; ++i) + { + if (i == sg->n_seq || (group[i]>>32) != (group[x]>>32)) + { + append_tig_link(group, x, i - x, link, sg, &buf); + x = i; + } + } + + free(flag); kv_destroy(buf.a); +} + + void measure_distance(const ma_ug_t* ug, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, kv_u_trans_t *ta) { // double index_time = yak_realtime(); @@ -3877,16 +4130,16 @@ void measure_distance(const ma_ug_t* ug, kvec_pe_hit* hits, hc_links* link, bubb push_hc_edge(&(link->a.a[k]), i, 0, 0, &t_d); } } - } - + asg_t *copy_sg = copy_read_graph(ug->g); - update_ug_by_trans(copy_sg, ta); + // update_ug_by_tigs(copy_sg, link); all_pair_shortest_path(copy_sg, link, &M); fill_utg_distance_multi(copy_sg, link, &M, bub); update_containment_distance(copy_sg, ta, link); + // update_dis_connected_gfa(copy_sg, link, &M); asg_destroy(copy_sg); @@ -5634,6 +5887,8 @@ inline double get_trans(const ha_ug_index* idx, uint64_t x) inline double get_trans_weight_advance(const ha_ug_index* idx, uint64_t x, trans_idx* dis) { long double rate = 0; + + if(x == (uint64_t)-1) x = dis->med; if(x < dis->max) { uint64_t i; @@ -6106,6 +6361,7 @@ uint64_t rLEN, uint64_t r_hetLen, int beg_idx, int end_idx, uint32_t check_het) if(beg_idx > end_idx) return; uint64_t r_homLen = rLEN - r_hetLen; + ///pLen: total length, rLEN: beg/sink length if(pLen > 0 && rLEN > 0 && r_hetLen > 0 && rLEN < pLen*0.5 && r_hetLen < rLEN * 0.2) ///set het to hom { set_bub_idx(bub, bu, untig_sg, beg_idx, end_idx, 1, bub->check_het); @@ -6140,7 +6396,7 @@ void detect_bub_graph(bubble_type* bub, asg_t *untig_sg) t = NULL; if(k+1 < u->n) t = &(arc_first(bg, u->a[k]>>32)); - pLen += bg->seq[rId].len; + pLen += bg->seq[rId].len;///path length in bubble if(end_idx < beg_idx) ///first bubble { pLen += untig_sg->seq[r_root>>1].len; @@ -8764,6 +9020,8 @@ void build_bub_graph(ma_ug_t* ug, bubble_type* bub) // detect_bub_graph(bub, ug->g, 1); update_bubble_chain(ug, bub, 1, 0); ///print_bubble_chain(bub, "second round"); + update_bubble_chain(ug, bub, 0, 1); + resolve_bubble_chain_tangle(ug, bub); } void get_forward_distance(uint32_t src, uint32_t dest, asg_t *sg, hc_links* link, MT* M) @@ -8858,6 +9116,7 @@ H_partition* hap, int8_t *s, trans_idx* dis) } buf.n = MIN(r_idx, f_idx); + trans_p_t* p = NULL; dis->n = 0; uint64_t bin_size = MIN(2250, buf.n>>8), m; @@ -8968,7 +9227,7 @@ H_partition* hap, int8_t *s, trans_idx* dis) dis->max = dis->a[dis->n-1].end; - kv_destroy(buf); + if(idx->a < 0) idx->a = 0; if(idx->a == 0) { @@ -8982,6 +9241,39 @@ H_partition* hap, int8_t *s, trans_idx* dis) // fprintf(stderr, "idx->a: %f, idx->b: %f, idx->frac: %f, med: %lu\n", // (double)idx->a, (double)idx->b, (double)idx->frac, med); + + + + + buf.n = 0; + for (k = 0; k < hits->a.n; ++k) + { + beg = ((hits->a.a[k].s<<1)>>(64 - idx->uID_bits)); + end = ((hits->a.a[k].e<<1)>>(64 - idx->uID_bits)); + + if(IF_HOM(beg, *bub)) continue; + if(IF_HOM(end, *bub)) continue; + if(beg == end) continue; + + t_d = get_hic_distance(&(hits->a.a[k]), link, idx); + if(t_d == (uint64_t)-1) continue; + kv_push(uint64_t, buf, t_d); + } + + ///might have bias, we may not use right linkage larger than trans rc linkage + radix_sort_hc64(buf.a, buf.a+buf.n); + if(buf.n == 0) + { + dis->med = 0; + } + else + { + dis->med = ((buf.n&1)?buf.a[buf.n>>1]:((buf.a[buf.n>>1]+buf.a[(buf.n>>1)-1])/2)); + } + + fprintf(stderr, "dis->med: %lu\n", dis->med); + + kv_destroy(buf); return 1; } @@ -12715,17 +13007,20 @@ void reorder_bubbles(bubble_type* bub, kv_u_trans_t *ta, uint32_t n_utg) void update_trans_g(ha_ug_index* idx, kv_u_trans_t *ta, bubble_type* bub) { // double index_time = yak_realtime(); - update_bubble_chain(idx->ug, bub, 0, 1); + // update_bubble_chain(idx->ug, bub, 0, 1); - resolve_bubble_chain_tangle(idx->ug, bub); + // resolve_bubble_chain_tangle(idx->ug, bub); clean_bubble_chain_by_hic(idx->ug, ta, bub); // print_debug_bubble_graph(bub, idx->ug, "bub-2"); - append_boundary_chain_hic(idx->ug, ta, bub); + // append_boundary_chain_hic(idx->ug, ta, bub); - reorder_bubbles(bub, ta, idx->ug->g->n_seq); + fprintf(stderr, "s_bub: %lu, f_bub: %lu, b_bub: %lu, b_end_bub: %lu, tangle_bub: %lu, cross_bub: %lu, mess_bub: %lu\n", + bub->s_bub, bub->f_bub, bub->b_bub, bub->b_end_bub, bub->tangle_bub, bub->cross_bub, bub->mess_bub); + + ///reorder_bubbles(bub, ta, idx->ug->g->n_seq); // fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time); } @@ -13305,8 +13600,8 @@ void print_debug_bubble_graph(bubble_type* bub, ma_ug_t* ug, const char *fn) void print_bubble_chain(bubble_type* bub) { - uint32_t m, i; - uint32_t beg, sink; + uint32_t m, i, k; + uint32_t beg, sink, *a = NULL, n; uint64_t bid; ma_utg_t *u = NULL; for (m = 0; m < bub->chain_weight.n; m++) @@ -13317,9 +13612,16 @@ void print_bubble_chain(bubble_type* bub) for (i = 0; i < u->n; i++) { bid = u->a[i]>>33; - get_bubbles(bub, bid, &beg, &sink, NULL, NULL, NULL); + get_bubbles(bub, bid, &beg, &sink, &a, &n, NULL); fprintf(stderr, "btg%.6lu%c, beg-utg%.6ul, sink-utg%.6ul\n", bid, "fb"[bidf_bub?0:1], (beg>>1)+1, (sink>>1)+1); + for (k = 0; k < n; k++) + { + if(k != 0 && (k%5)==0) fprintf(stderr, "\n"); + fprintf(stderr, "m-utg%.6ul\t", (a[k]>>1)+1); + } + fprintf(stderr, "\n"); + } } } @@ -13444,13 +13746,13 @@ int alignment_worker_pipeline(sldat_t* sl, const enzyme *fn1, const enzyme *fn2) return 1; } -void debug_gfa_space(ma_ug_t* ug, trans_chain* t_ch) +void debug_gfa_space(ma_ug_t* ug, trans_chain* t_ch, kv_u_trans_t *ref) { bubble_type bub; memset(&bub, 0, sizeof(bubble_type)); bub.round_id = 0; bub.n_round = 2; - identify_bubbles(ug, &bub, t_ch->is_r_het); + identify_bubbles(ug, &bub, t_ch->is_r_het, ref); hc_links link; init_hc_links(&link, ug->g->n_seq, t_ch); @@ -13583,17 +13885,14 @@ kv_u_trans_t *ta, trans_idx* dis) if(IF_HOM(end, *bub)) continue; t_d = get_hic_distance(&(hits->a.a[k]), link, idx); - if(t_d == (uint64_t)-1) continue; + // if(t_d == (uint64_t)-1) continue; get_u_trans_spec(ta, beg, end, &e1, NULL); get_u_trans_spec(ta, end, beg, &e2, NULL); if(e1 == NULL || e2 == NULL) continue; weight = 1; - if(dis) - { - weight = get_trans_weight_advance(idx, t_d, dis); - } + if(dis) weight = get_trans_weight_advance(idx, t_d, dis); e1->nw -= weight; e1->occ++; e2->nw -= weight; e2->occ++; @@ -13656,7 +13955,7 @@ pe_hit *hits, uint32_t occ, uint32_t qid, uint32_t qs, uint32_t qe, uint32_t tid if(!(ts <= e_beg && te >= e_end)) continue; t_d = get_hic_distance(&hits[k], link, idx); - if(t_d == (uint64_t)-1) continue; + // if(t_d == (uint64_t)-1) continue; weight = 1; if(dis) weight = get_trans_weight_advance(idx, t_d, dis); @@ -13676,7 +13975,7 @@ pe_hit *hits, uint32_t occ, uint32_t qid, uint32_t qs, uint32_t qe, uint32_t tid if(!(ts <= e_beg && te >= e_end)) continue; t_d = get_hic_distance(&hits[k], link, idx); - if(t_d == (uint64_t)-1) continue; + // if(t_d == (uint64_t)-1) continue; weight = 1; if(dis) weight = get_trans_weight_advance(idx, t_d, dis); @@ -13760,7 +14059,7 @@ kv_u_trans_t *ta, kv_u_trans_t *ref, trans_idx* dis) ta->a[m].nw /= (double)(MIN(hits->occ.a[ta->a[m].qn], hits->occ.a[ta->a[m].tn])); } /*******************************for debug************************************/ - trans_checking_pass(bub, ref, ta->a[i].qn, ta->a[i].tn); + // trans_checking_pass(bub, ref, ta->a[i].qn, ta->a[i].tn); // if(trans_checking_pass(bub, ref, ta->a[i].qn, ta->a[i].tn)) ta->a[m].nw = 0; /*******************************for debug************************************/ @@ -13856,21 +14155,6 @@ ha_ug_index* idx, bubble_type* bub, int8_t *s, uint32_t ignore_dis) is_comples_weight = get_trans_rate_function_advance(idx, hits, lk, bub, NULL, s, &dis); } - hc_edge *e = NULL; - for (i = 0; i < lk->a.n; i++) - { - for (k = 0; k < lk->a.a[i].e.n; k++) - { - if(lk->a.a[i].e.a[k].del) continue; - if(lk->a.a[i].e.a[k].dis == (uint64_t)-1) - { - e = get_hc_edge(lk, lk->a.a[i].e.a[k].uID, i, 0); - if(!e) fprintf(stderr, "ERROR\n"); - e->del = lk->a.a[i].e.a[k].del = 1; - } - } - } - for (i = 0; i < lk->a.n; i++) { for (k = m = 0; k < lk->a.a[i].e.n; k++) @@ -13958,11 +14242,10 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) bub.round_id = 0; bub.n_round = 2; for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++) { - identify_bubbles(idx->ug, &bub, idx->t_ch->is_r_het); - if(bub.round_id == 0) + identify_bubbles(idx->ug, &bub, idx->t_ch->is_r_het, &(idx->t_ch->k_trans)); + if(bub.round_id == 0) { measure_distance(idx->ug, &sl.hits, &link, &bub, &(idx->t_ch->k_trans)); - collect_hc_reverse_links(&link, idx->ug, &bub); } renew_kv_u_trans(&k_trans, &link, &sl.hits, &(idx->t_ch->k_trans), idx, &bub, s->s, 0); @@ -13991,7 +14274,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) ///print_debug_bubble_graph(&bub, idx->ug, asm_opt.output_file_name); - // print_bubble_chain(&bub); + print_bubble_chain(&bub); // destory_contig_partition(&hap); kv_destroy(back_hc_edge.a); kv_destroy(sl.hits.a); diff --git a/hic.h b/hic.h index 93e6a72..fe1d166 100644 --- a/hic.h +++ b/hic.h @@ -10,7 +10,7 @@ #define RC_2 2 hc_edge* get_hc_edge(hc_links* link, uint64_t src, uint64_t dest, uint64_t dir); -void push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t* d); +hc_edge* push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t* d); void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch); void hic_benchmark(ma_ug_t *ug, asg_t* read_g); @@ -67,13 +67,12 @@ void get_bubbles(bubble_type* bub, uint64_t id, uint32_t* beg, uint32_t* sink, u int load_hc_links(hc_links* link, const char *fn); void write_hc_links(hc_links* link, const char *fn); void destory_bubbles(bubble_type* bub); -void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag); +void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_trans_t *ref); void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub); uint32_t connect_bub_occ(bubble_type* bub, uint32_t root_id, uint32_t check_het); void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, uint32_t check_het); void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint32_t is_end); void set_b_utg_weight_flag(bubble_type* bub, buf_t* b, uint32_t v, uint8_t* vis_flag, uint32_t flag, uint32_t* occ); -uint32_t get_unitig_het_arb(ma_utg_t* u, uint8_t *r_het_flag, uint32_t m_het_label, uint32_t p_het_label, uint32_t n_het_label); void debug_gfa_space(ma_ug_t* ug, hap_cov_t *cov); #endif