From a39f01f4d82f8dcb9c4b25d7b9f8e070e19c28d1 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Mon, 24 May 2021 19:43:50 -0400 Subject: [PATCH] haplotype popping --- Makefile | 4 +- Overlaps.cpp | 658 +++++----------------- Overlaps.h | 44 +- Purge_Dups.cpp | 30 +- Purge_Dups.h | 2 +- hic.cpp | 64 ++- hic.h | 5 +- tovlp.cpp | 1441 ++++++++++++++++++++++++++++++++++++++++++++++++ tovlp.h | 15 + 9 files changed, 1696 insertions(+), 567 deletions(-) create mode 100644 tovlp.cpp create mode 100644 tovlp.h diff --git a/Makefile b/Makefile index 4e9eb3b..030ab95 100644 --- a/Makefile +++ b/Makefile @@ -6,7 +6,8 @@ CPPFLAGS= INCLUDES= OBJS= CommandLines.o Process_Read.o Assembly.o Hash_Table.o \ POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \ - htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o rcut.o horder.o + htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o rcut.o horder.o \ + tovlp.o EXE= hifiasm LIBS= -lz -lpthread -lm @@ -74,3 +75,4 @@ sys.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h hic.o: hic.h rcut.o: rcut.h horder.o: horder.h +tovlp.o: tovlp.h diff --git a/Overlaps.cpp b/Overlaps.cpp index b9f61fe..42a8125 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -11,6 +11,7 @@ #include "Purge_Dups.h" #include "hic.h" #include "kthread.h" +#include "tovlp.h" uint32_t debug_purge_dup = 0; @@ -66,6 +67,13 @@ long long min_thres; uint32_t print_untig_by_read(ma_ug_t *g, const char* name, uint32_t in, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, const char* info); +int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov, utg_trans_t *o, uint32_t is_update_chain); +void get_utg_ovlp(ma_ug_t **ug, asg_t* read_g, float drop_rate, +ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, +long long tipsLen, float tip_drop_ratio, long long stops_threshold, +R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, +kvec_asg_arc_t_warp* new_rtg_edges, hap_cov_t **i_cov, bub_label_t* b_mask_t, +uint32_t collect_p_trans, uint32_t collect_p_trans_f); void init_bub_label_t(bub_label_t* x, uint32_t n_thres, uint32_t n_reads) { @@ -11798,7 +11806,7 @@ inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t igno return len; } -uint32_t set_utg_offset(uint32_t *a, uint32_t a_n, ma_ug_t *ug, asg_t *read_sg, hap_cov_t *cov, uint32_t is_clear, +uint32_t set_utg_offset(uint32_t *a, uint32_t a_n, ma_ug_t *ug, asg_t *read_sg, uint64_t* pos_idx, uint32_t is_clear, uint32_t only_len) { uint32_t ori, uid, v, nv, l, k; @@ -11840,13 +11848,13 @@ uint32_t only_len) if(is_clear == 1) { - cov->pos_idx[v>>1] = (uint64_t)-1; + pos_idx[v>>1] = (uint64_t)-1; } else { - cov->pos_idx[v>>1] = len; - cov->pos_idx[v>>1] <<= 32; - cov->pos_idx[v>>1] |= (uint64_t)v; + pos_idx[v>>1] = len; + pos_idx[v>>1] <<= 32; + pos_idx[v>>1] |= (uint64_t)v; } } } @@ -11879,11 +11887,6 @@ KRADIX_SORT_INIT(origin_trans_sort, asg_arc_t_offset, origin_trans_key, member_s KRADIX_SORT_INIT(origin_trans_el_sort, asg_arc_t_offset, origin_trans_el_key, 1) -inline uint32_t get_offset_adjust(uint32_t offset, uint32_t offsetLen, uint32_t targetLen) -{ - return ((double)(offset)/(double)(offsetLen))*targetLen; -} - void refine_u_trans_t(u_trans_hit_t *q, kv_ca_buf_t* cb) { ///already know [qScur, qEcur), [qSpre, qEpre) @@ -12320,10 +12323,10 @@ ma_ug_t *ug, uint32_t flag, double score, const char* cmd) ca_buf_t *tx = NULL; if(i_pri_len) pri_len = (*i_pri_len); - else pri_len = set_utg_offset(pri_a, pri_n, ug, read_sg, cov, 0, 1); + else pri_len = set_utg_offset(pri_a, pri_n, ug, read_sg, cov->pos_idx, 0, 1); if(i_aux_len) aux_len = (*i_aux_len); - else aux_len = set_utg_offset(aux_a, aux_n, ug, read_sg, cov, 0, 1); + else aux_len = set_utg_offset(aux_a, aux_n, ug, read_sg, cov->pos_idx, 0, 1); tt = NULL; t_ch->c_buf.n = 0; t_ch->k_t_b.n = 0; if(tailIndex->a.n > 0) tt = &(u_buffer->a.a[tailIndex->a.a[0]]); @@ -12475,8 +12478,8 @@ ma_ug_t *ug, asg_t *read_sg, hap_cov_t *cov) 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); + len_aux = set_utg_offset(aux->b.a, aux->b.n, ug, read_sg, cov->pos_idx, 0, 0); + chain_trans_ovlp(cov, NULL, ug, read_sg, pri, len_aux, &thre_pri); if(thre_pri > 0) { /*******************************for debug************************************/ @@ -12540,124 +12543,9 @@ ma_ug_t *ug, asg_t *read_sg, hap_cov_t *cov) if(occ >= thre_pri) break; } } - set_utg_offset(aux->b.a, aux->b.n, ug, read_sg, cov, 1, 0); + set_utg_offset(aux->b.a, aux->b.n, ug, read_sg, cov->pos_idx, 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) @@ -14026,7 +13914,10 @@ bub_label_t* b_mask_t) hap_cov_t *cov = NULL; asg_t *copy_sg = copy_read_graph(sg); ma_ug_t *copy_ug = copy_untig_graph(ug); - adjust_utg_by_primary(©_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, + // adjust_utg_by_primary(©_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, + // tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, + // max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 1); + get_utg_ovlp(©_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 1); print_utg(copy_ug, copy_sg, coverage_cut, output_file_name, sources, ruIndex, max_hang, @@ -14724,350 +14615,9 @@ void label_r_set(buf_t* b, R_to_U* ruIndex, ma_ug_t *ug, uint32_t flag) } } -/** -void get_trans_n(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint32_t uid, -uint32_t matchLen, uint32_t matchOcc, uint32_t *max_count, uint32_t *min_count) -{ - uint32_t i, j, qn, tn, is_Unitig, uId; - (*max_count) = (*min_count) = 0; - ma_utg_t *u = NULL; - u = &(ug->u.a[uid>>1]); - for (i = 0; i < u->n; i++) - { - qn = u->a[i]>>33; - if(reverse_sources[qn].length > 0) (*min_count)++; - for (j = 0; j < (long long)reverse_sources[qn].length; j++) - { - tn = Get_tn(reverse_sources[qn].buffer[j]); - if(rg->seq[tn].del == 1) - { - get_R_to_U(ruIndex, tn, &tn, &is_Unitig); - if(tn == (uint32_t)-1 || is_Unitig == 1 || rg->seq[tn].del == 1) continue; - } - - get_R_to_U(ruIndex, tn, &uId, &is_Unitig); - if(uId!=(uint32_t)-1 && is_Unitig == 1) - { - (*max_count)++; - break; - } - } - } -} - -void dfs_path(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint32_t v0, uint32_t fv, -kv_tip_t *tb, asg_t *g, uint32_t max_dist, uint32_t tigLen, uint32_t tigOcc) -{ - asg_arc_t *av = NULL; - tip_t *src = NULL, *t = NULL; - uint32_t i, v, nv, w, l, tot, ma, to_replace; - long long cur_w, max_w; - tb->st.n = tb->r.n = 0; - tb->b[v0].d = tb->b[v0].ma = tb->b[v0].tot = 0; tb->b[v0].p = (uint32_t)-1; tb->b[v0].in = 0; - kv_push(uint32_t, tb->st, v0); - kv_push(uint32_t, tb->r, v0); - - while (tb->st.n > 0) - { - tb->st.n--; - v = tb->st.a[tb->st.n]; - src = &(tb->b[v]); - - av = asg_arc_a(g, v); - nv = asg_arc_n(g, v); - for (i = 0; i < nv; ++i) - { - if (av[i].del) continue; - if (av[i].v == fv) continue; - if (av[i].v == v0) goto dfs_end; - w = av[i].v; - l = (uint32_t)av[i].ul; - t = &(tb->b[w]); - if(src->d + l > max_dist) continue; - get_trans_n(ug, rg, reverse_sources, ruIndex, w, tigLen, tigOcc, &ma, &tot); - if(t->vis == 0) - { - kv_push(uint32_t, tb->r, w); - t->p = v, t->vis = 1, t->d = src->d + l; - t->ma = ma + src->ma; - t->tot = tot + src->tot; - } - else - { - to_replace = 0; - cur_w = (long long)((src->ma + ma)*2) - (long long)(src->tot + tot); - max_w = (long long)(t->ma*2) - (long long)(t->tot); - if(cur_w > max_w) to_replace = 1; - else if(cur_w == max_w && (src->d + l) > t->d) to_replace = 1; - - if(to_replace) - { - t->p = v; - t->d = src->d + l; - t->ma = ma + src->ma; - t->tot = tot + src->tot; - } - } - } - } - - dfs_end: -} -**/ - -/** -int asg_arc_complex_tip(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, -R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov) -{ - double startTime = Get_T(); - uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, in, flag, operation; - uint32_t return_flag, convex_i, k; - long long ll, tmp, max_stop_nodeLen, max_stop_baseLen; - kv_tip_t tb; kv_init(tb); - - buf_t b, b_0, b_1; - memset(&b, 0, sizeof(buf_t)); - memset(&b_0, 0, sizeof(buf_t)); - memset(&b_1, 0, sizeof(buf_t)); - for (v = 0; v < n_vtx; ++v) - { - uint32_t i; - if(g->seq[v>>1].del || g->seq[v>>1].c == ALTER_LABLE) continue; - ///tip - if (get_real_length(g, v, NULL) != 0) continue; - if(get_real_length(g, v^1, NULL) != 1) continue; - - b.b.n = 0; - return_flag = get_unitig(g, ug, v^1, &convex, &ll, &tmp, &max_stop_nodeLen, - &max_stop_baseLen, 1, &b); - - if(return_flag != MUL_INPUT) continue; - in = convex^1; - get_real_length(g, convex, &convex); - convex = convex^1; - uint32_t n_convex = asg_arc_n(g, convex), convexLen = ll; - asg_arc_t *a_convex = asg_arc_a(g, convex); - for (i = 0; i < n_convex; i++) - { - if(a_convex[i].del) continue; - if(a_convex[i].v == in) break; - } - convex_i = i; - - label_r_set(&b, ruIndex, ug, 1); - tb.n = 0; - - - - - - - label_r_set(&b, ruIndex, ug, (uint32_t)-1); - - for (i = 0; i < n_convex; i++) - { - if(a_convex[i].del) continue; - if(i == convex_i) continue; - - return_flag = get_unitig(g, ug, a_convex[i].v, &convex, &ll, &tmp, &max_stop_nodeLen, - &max_stop_baseLen, stops_threshold, NULL); - - if(convexLen < ll*drop_ratio && max_stop_nodeLen >= ll*MAX_STOP_RATE) - { - n_reduced++; - operation = TRIM; - flag = check_different_haps_base(g, ug, read_sg, a_convex[convex_i].v, a_convex[i].v, - reverse_sources, &b_0, &b_1, ruIndex, min_edge_length, stops_threshold); - // #define UNAVAILABLE (uint32_t)-1 - // #define PLOID 0 - // #define NON_PLOID 1 - if(flag == NON_PLOID) operation = CUT; - - for (k = 0; k < b.b.n; k++) - { - g->seq[b.b.a[k]>>1].c = ALTER_LABLE; - } - - for (k = 0; k < b.b.n; k++) - { - asg_seq_drop(g, b.b.a[k]>>1); - } - - if(operation == CUT) - { - for (k = 0; k < b.b.n; k++) - { - g->seq[b.b.a[k]>>1].c = CUT; - } - } - - ///note: we need to remove b_0, insetad of b_1 here - if(cov && operation != CUT) - { - collect_trans_cov(__func__, &b_1, a_convex[i].ol, &b_0, a_convex[convex_i].ol, ug, read_sg, cov); - } - - - break; - } - } - } - - - if(VERBOSE >= 1) - { - fprintf(stderr, "[M::%s] removed %d long tips\n", - __func__, n_reduced); - fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime); - } - - asg_cleanup(g); - asg_symm(g); - free(b.b.a); - free(b_0.b.a); - free(b_1.b.a); - kv_destory(tb); - - return n_reduced; -} -**/ -uint32_t dfs_set(asg_t *g, uint32_t v0, uint32_t fbv, kvec_t_u32_warp *stack, kvec_t_u32_warp *tmp, uint8_t *vis, uint8_t flag) -{ - uint32_t cur, ncur, i, occ = 0;; - asg_arc_t *acur = NULL; - stack->a.n = 0; - kv_push(uint32_t, stack->a, v0); - // vis[v0] = 0; - - while (stack->a.n > 0) - { - stack->a.n--; - cur = stack->a.a[stack->a.n]; - if(vis[cur] && (!(vis[cur]&flag))) - { - vis[cur] |= flag; - kv_push(uint32_t, tmp->a, cur); - occ++; - } - if(vis[cur]) continue; - else kv_push(uint32_t, tmp->a, cur); - - vis[cur] |= flag; - ncur = asg_arc_n(g, cur); - acur = asg_arc_a(g, cur); - for (i = 0; i < ncur; i++) - { - if(acur[i].del) continue; - if((acur[i].v>>1) == fbv) continue; - if(vis[acur[i].v] && (!(vis[acur[i].v]&flag))) - { - vis[acur[i].v] |= flag; - kv_push(uint32_t, tmp->a, acur[i].v); - occ++; - continue; - } - - kv_push(uint32_t, stack->a, acur[i].v); - } - } - - return occ; -} - - -int asg_arc_decompress(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, -R_to_U* ruIndex, hap_cov_t *cov) -{ - double startTime = Get_T(); - asg_arc_t *as = NULL; - uint32_t v, s, sv, nc, ns, n_vtx = g->n_seq * 2, n_reduced = 0, convex; - uint32_t return_flag, k, i, p_n, a_n, *a_a = NULL;; - uint8_t fp = 1, fa = 2, found; - long long ll, tmp, max_stop_nodeLen, max_stop_baseLen; - kvec_t_u32_warp tt; kv_init(tt.a); - kvec_t_u32_warp stack; kv_init(stack.a); - uint8_t *vis = NULL; CALLOC(vis, n_vtx); - buf_t b; - memset(&b, 0, sizeof(buf_t)); - pdq pq_p, pq_a; - init_pdq(&pq_p, g->n_seq<<1); - init_pdq(&pq_a, g->n_seq<<1); - - for (v = 0; v < n_vtx; v++) - { - if(g->seq[v>>1].del || g->seq[v>>1].c == ALTER_LABLE) continue; - if(asg_arc_n(g, v) == 0 || get_real_length(g, v, NULL) != 1) continue; - get_real_length(g, v, &s); - if(get_real_length(g, s^1, NULL) < 2) continue; - return_flag = get_unitig(g, ug, v^1, &convex, &ll, &tmp, &max_stop_nodeLen,&max_stop_baseLen, 1, NULL); - if(return_flag == LOOP) continue; - get_real_length(g, convex, &convex); - if(get_real_length(g, convex^1, NULL) < 2) continue; - tt.a.n = 0; s^=1; sv = v^1; - dfs_set(g, sv, s>>1, &stack, &tt, vis, fp); - p_n = tt.a.n; - as = asg_arc_a(g, s); ns = asg_arc_n(g, s); found = 0; - - for (i = 0; i < ns; i++) - { - if(as[i].del || as[i].v == sv) continue; - nc = dfs_set(g, as[i].v, s>>1, &stack, &tt, vis, fa); - // if((s>>1)==9882) fprintf(stderr, "s-%u, as[i].v-%u, sv-%u, nc-%u\n", s, as[i].v, sv, nc); - if(nc && check_trans_relation_by_path(sv, as[i].v, &pq_p, &pq_a, g, - vis, fp+fa, nc, NULL, 0.45)) - { - found = 1; - } - a_a = tt.a.a + p_n; a_n = tt.a.n - p_n; - for (k = 0; k < a_n; k++) - { - if(vis[a_a[k]]&fa) vis[a_a[k]] -= fa; - } - tt.a.n = p_n; - if(found) break; - } - if(found) - { - b.b.n = 0; - get_unitig(g, ug, sv, &convex, &ll, &tmp, &max_stop_nodeLen,&max_stop_baseLen, 1, &b); - for (k = 0; k < b.b.n; k++) - { - g->seq[b.b.a[k]>>1].c = ALTER_LABLE; - } - - for (k = 0; k < b.b.n; k++) - { - asg_seq_drop(g, b.b.a[k]>>1); - } - - // fprintf(stderr, "++++++++utg%.6ul\n", (sv>>1)+1); - } - a_a = tt.a.a; a_n = p_n; - for (k = 0; k < a_n; k++) vis[a_a[k]] = 0; - // fprintf(stderr, "----------utg%.6ul (len: %u)\n", (sv>>1)+1, g->seq[sv>>1].len); - } - - if(VERBOSE >= 1) - { - fprintf(stderr, "[M::%s] removed %d long tips\n", - __func__, n_reduced); - fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime); - } - - asg_cleanup(g); - asg_symm(g); - free(b.b.a); - kv_destroy(tt.a); - kv_destroy(stack.a); - free(vis); - destory_pdq(&pq_p); - destory_pdq(&pq_a); - - return n_reduced; -} int asg_arc_cut_trio_long_tip_primary(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, -R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov) +R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov, utg_trans_t *o) { double startTime = Get_T(); ///the reason is that each read has two direction (query->target, target->query) @@ -15160,6 +14710,11 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov) { collect_trans_cov(__func__, &b_0, av[v_maxLen_i].ol, &b_1, av[i].ol, ug, read_sg, cov); } + + if(o && operation != CUT) + { + collect_trans_ovlp(__func__, &b_0, av[v_maxLen_i].ol, &b_1, av[i].ol, ug, o); + } } } @@ -15183,7 +14738,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov) } int asg_arc_cut_trio_long_tip_primary_complex(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, -R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_threshold, hap_cov_t *cov) +R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_threshold, hap_cov_t *cov, utg_trans_t *o) { double startTime = Get_T(); uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, in, flag, operation; @@ -15263,6 +14818,11 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_thre { collect_trans_cov(__func__, &b_1, a_convex[i].ol, &b_0, a_convex[convex_i].ol, ug, read_sg, cov); } + + if(o && operation != CUT) + { + collect_trans_ovlp(__func__, &b_1, a_convex[i].ol, &b_0, a_convex[convex_i].ol, ug, o); + } break; @@ -15288,7 +14848,8 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_thre } void renew_longest_tip_by_drop(asg_t *g, ma_ug_t *ug, asg_arc_t *av, uint32_t nv, -long long* base_maxLen, long long* base_maxLen_i, uint32_t stops_threshold, buf_t* b, uint32_t trio_flag) +long long* base_maxLen, long long* base_maxLen_i, uint32_t stops_threshold, buf_t* b, +uint32_t trio_flag) { if(trio_flag != FATHER && trio_flag != MOTHER) return; Trio_counter max, cur; @@ -15363,7 +14924,7 @@ long long* base_maxLen, long long* base_maxLen_i, uint32_t stops_threshold, buf_ int asg_arc_cut_trio_long_equal_tips_assembly(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_t trio_flag, -hap_cov_t *cov) +hap_cov_t *cov, utg_trans_t *o) { double startTime = Get_T(); uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, flag, is_hap, n_tips, return_flag, k; @@ -15451,6 +15012,11 @@ hap_cov_t *cov) { collect_trans_cov(__func__, &b_0, av[base_maxLen_i].ol, &b_1, av[i].ol, ug, read_sg, cov); } + + if(o) + { + collect_trans_ovlp(__func__, &b_0, av[base_maxLen_i].ol, &b_1, av[i].ol, ug, o); + } is_hap++; @@ -15682,7 +15248,8 @@ R_to_U* ruIndex, uint32_t positive_flag, float drop_rate) } int asg_arc_cut_trio_long_equal_tips_assembly_complex(asg_t *g, ma_ug_t *ug, asg_t *read_sg, -ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_t stops_threshold, hap_cov_t *cov) +ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_t stops_threshold, +hap_cov_t *cov, utg_trans_t *o) { double startTime = Get_T(); uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, in, flag; @@ -15764,6 +15331,11 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_ { collect_trans_cov(__func__, &b_1, a_convex[i].ol, &b_0, a_convex[convex_i].ol, ug, read_sg, cov); } + + if(o) + { + collect_trans_ovlp(__func__, &b_1, a_convex[i].ol, &b_0, a_convex[convex_i].ol, ug, o); + } ///lable the primary one @@ -15800,7 +15372,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_ int detect_chimeric_by_topo(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, long long miniedgeLen, uint32_t stops_threshold, float drop_rate, -R_to_U* ruIndex) +R_to_U* ruIndex, utg_trans_t *o) { double startTime = Get_T(); uint32_t i, k, v_i, v_beg, v_end, selfLen, w1, w2, wv, nw, n_vtx = g->n_seq * 2, n_reduced = 0, convex, convex_T, read_num; @@ -15947,6 +15519,7 @@ R_to_U* ruIndex) for (k = 0; k < b_0.b.n; k++) { asg_seq_drop(g, b_0.b.a[k]>>1); + if(o) asg_seq_del(o->cug->g, b_0.b.a[k]>>1); } if(read_num <= CHIMERIC_TRIM_THRES) @@ -15959,6 +15532,7 @@ R_to_U* ruIndex) } + if(o) asg_cleanup(o->cug->g); asg_cleanup(g); free(b_0.b.a); free(b_1.b.a); @@ -16165,7 +15739,7 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov) redo: ///print_untig((ug), 61955, "i-0:", 0); - asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov, 1); + asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov, NULL, 1); magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate); /**********debug**********/ if(just_bubble_pop == 0) @@ -16179,16 +15753,16 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov) { pre_cons = get_graph_statistic(g); ///need consider tangles - asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov, 1); + asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov, NULL, 1); /**********debug**********/ if(just_bubble_pop == 0) { ///need consider tangles - asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov); - asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, cov); - asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov); - asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov); - detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex); + asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, NULL); + asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, cov, NULL); + asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL); + asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL); + detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL); ///need consider tangles ///note we need both the read graph and the untig graph } @@ -16245,7 +15819,7 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov) redo: ///print_graph_statistic(g, "beg"); ///print_debug_gfa(read_g, ug, coverage_cut, "debug_trans_ovlp_hg002", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); - asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, 1); + asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, NULL, 1); if(just_bubble_pop == 0) { cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 2); @@ -16256,20 +15830,16 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov) while(pre_cons != cur_cons) { pre_cons = get_graph_statistic(g); - asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, 1); + asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, NULL, 1); if(just_bubble_pop == 0) { ///need consider tangles - asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov); - asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, cov); - asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov); - asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov); - detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex); - - if(asm_opt.polyploidy > 2) - { - asg_arc_decompress(g, ug, read_g, reverse_sources, ruIndex, cov); - } + asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, NULL); + asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, cov, NULL); + asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL); + asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL); + detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL); + if(round != T_ROUND) { unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, @@ -16285,7 +15855,7 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov) } resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, (uint32_t)-1, drop_ratio); drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex); - 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); ///print_graph_statistic(g, "end"); if(round > 0) @@ -16300,45 +15870,60 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov) } -void topo_ovlp_collect(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* sources, + +utg_trans_t *topo_ovlp_collect(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, long long tipsLen, float tip_drop_ratio, -long long stops_threshold, R_to_U* ruIndex, float chimeric_rate, float drop_ratio, hap_cov_t *cov) +long long stops_threshold, R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, +int min_ovlp, hap_cov_t *cov) { - // kv_u_trans_t *k_trans + utg_trans_t *o = init_utg_trans_t(ug, reverse_sources, coverage_cut, ruIndex, read_g, max_hang, min_ovlp); #define T_ROUND 2 asg_t *g = ug->g; int round = T_ROUND; + // print_debug_gfa(read_g, ug, coverage_cut, "debug_init", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); + redo: - asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, 1); + asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, o, 1); cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 2); long long pre_cons = get_graph_statistic(g); long long cur_cons = 0; while(pre_cons != cur_cons) - { - pre_cons = get_graph_statistic(g); - asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, 1); - - ///need consider tangles - asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov); - asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, cov); - asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov); - asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov); - detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex); - - if(asm_opt.polyploidy > 2) + { + while(pre_cons != cur_cons) { - asg_arc_decompress(g, ug, read_g, reverse_sources, ruIndex, cov); + while(pre_cons != cur_cons) + { + pre_cons = get_graph_statistic(g); + asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, o, 1); + + ///need consider tangles + asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, o); + asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, cov, o); + asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, o); + asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, o); + detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, o); + + cur_cons = get_graph_statistic(g); + } + + // if(asm_opt.polyploidy > 2) asg_arc_decompress(g, ug, read_g, reverse_sources, ruIndex, o); + if(asm_opt.polyploidy > 2) + { + asg_arc_decompress_mul(g, ug, read_g, (uint32_t)-1, DROP, reverse_sources, ruIndex, o); + } + + cur_cons = get_graph_statistic(g); } if(round != T_ROUND) { unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, reverse_sources, 0, 1); } - cur_cons = get_graph_statistic(g); - } + } + cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 2); @@ -16356,7 +15941,8 @@ long long stops_threshold, R_to_U* ruIndex, float chimeric_rate, float drop_rati } round--; goto redo; - } + } + return o; } void set_drop_trio_flag(ma_ug_t *ug) @@ -17565,8 +17151,8 @@ void chain_origin_trans_uid_s_bubble(buf_t *pri, buf_t* aux, uint32_t beg, uint3 uint64_t pri_len, aux_len; priBeg = priEnd = auxBeg = auxEnd = (uint32_t)-1; - pri_len = set_utg_offset(pri->b.a, pri->b.n, ug, cov->read_g, cov, 0, 1); - aux_len = set_utg_offset(aux->b.a, aux->b.n, ug, cov->read_g, cov, 0, 1); + pri_len = set_utg_offset(pri->b.a, pri->b.n, ug, cov->read_g, cov->pos_idx, 0, 1); + aux_len = set_utg_offset(aux->b.a, aux->b.n, ug, cov->read_g, cov->pos_idx, 0, 1); pri_v = pri->b.a[0]; aux_v = aux->b.a[0]; av = asg_arc_a(ug->g, beg); @@ -18356,7 +17942,7 @@ int max_hang, int min_ovlp, long long gap_fuzz) for (i = n_arc = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs if (!av[i].del) ++n_arc; if (n_arc < 2) continue; - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, max_dist, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0)) + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, max_dist, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) { //beg is v, end is b.S.a[0] //note b.b include end, does not include beg @@ -18528,7 +18114,7 @@ uint32_t positive_flag, uint32_t negative_flag, uint32_t found) buf_t b_new; memset(&b_new, 0, sizeof(buf_t)); b_new.a = (binfo_t*)calloc(g->n_seq * 2, sizeof(binfo_t)); - uint32_t n_pop = asg_bub_pop1_primary_trio(g, utg, v0, max_dist, &b_new, positive_flag, negative_flag, 0, NULL, NULL, NULL, 0, 0); + uint32_t n_pop = asg_bub_pop1_primary_trio(g, utg, v0, max_dist, &b_new, positive_flag, negative_flag, 0, NULL, NULL, NULL, 0, 0, NULL); if(n_pop != found) fprintf(stderr, "ERROR\n"); free(b_new.a); free(b_new.S.a); free(b_new.T.a); free(b_new.b.a); free(b_new.e.a); @@ -18536,7 +18122,7 @@ uint32_t positive_flag, uint32_t negative_flag, uint32_t found) uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, uint64_t max_dist, buf_t *b, uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes, -hap_cov_t *cov, uint32_t is_update_chain, uint32_t keep_d) +hap_cov_t *cov, uint32_t is_update_chain, uint32_t keep_d, utg_trans_t *o) { uint32_t i, n_pending = 0, is_first = 1, cur_m, cur_c, cur_np, cur_nc, to_replace, n_tips, tip_end; uint64_t n_pop = 0; @@ -18733,6 +18319,7 @@ hap_cov_t *cov, uint32_t is_update_chain, uint32_t keep_d) ///if(keep_d != 0) debug_asg_bub_pop1_primary_trio(g, utg, v0, max_dist, b, positive_flag, negative_flag, 1); /****************************may have bugs********************************/ if(cov && utg) asg_bub_backtrack_primary_cov(utg, v0, b, cov, is_update_chain); + if(o && utg) asg_bub_collect_ovlp(utg, v0, b, o); if(is_pop) asg_bub_backtrack_primary(g, v0, b); if(path_base_len || path_nodes) asg_bub_backtrack_primary_length(g, utg, v0, b, path_base_len, path_nodes); @@ -19198,7 +18785,7 @@ uint64_t get_s_bub_pop_max_dist_advance(asg_t *g, buf_s_t *b) // pop bubbles -int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov, uint32_t is_update_chain) +int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov, utg_trans_t *o, uint32_t is_update_chain) { asg_t *g = ug->g; uint32_t v, n_vtx = g->n_seq * 2, n_arc, nv, i; @@ -19226,7 +18813,7 @@ int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t posi for (i = n_arc = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs if (!av[i].del) ++n_arc; if (n_arc < 2) continue; - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, max_dist, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0)) + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, max_dist, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) { //beg is v, end is b.S.a[0] //note b.b include end, does not include beg @@ -19250,7 +18837,7 @@ int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t posi for (i = n_arc = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs if (!av[i].del) ++n_arc; if (n_arc > 1) - n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, cov, is_update_chain, 0); + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, cov, is_update_chain, 0, o); } if(VERBOSE >= 1) @@ -22341,13 +21928,13 @@ uint32_t positive_flag, uint32_t negative_flag) v = beg; if((!g->seq[v>>1].del)&&(g->seq[v>>1].c!=ALTER_LABLE)&&get_real_length(g, v, NULL)>=2) { - n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, NULL, 0, 1); + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, NULL, 0, 1, NULL); } v = end^1; if((!g->seq[v>>1].del)&&(g->seq[v>>1].c!=ALTER_LABLE)&&get_real_length(g, v, NULL)>=2) { - n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, NULL, 0, 1); + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, NULL, 0, 1, NULL); } @@ -22362,7 +21949,7 @@ uint32_t positive_flag, uint32_t negative_flag) { v = v|k; if(get_real_length(g, v, NULL)<=1) continue; - n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, NULL, 0, 1); + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, NULL, 0, 1, NULL); } } @@ -24858,6 +24445,7 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov kvec_asg_arc_t_warp* new_rtg_edges, hap_cov_t **i_cov, bub_label_t* b_mask_t, uint32_t collect_p_trans, uint32_t collect_p_trans_f) { + fprintf(stderr, "******1******\n"); asg_t* nsg = (*ug)->g; uint32_t v, n_vtx = nsg->n_seq, k, rId, just_contain; ma_utg_t* u = NULL; @@ -24876,7 +24464,7 @@ uint32_t collect_p_trans, uint32_t collect_p_trans_f) } topo_ovlp_collect(*ug, read_g, sources, reverse_sources, coverage_cut, tipsLen, tip_drop_ratio, - stops_threshold, ruIndex, chimeric_rate, drop_ratio, cov); + stops_threshold, ruIndex, chimeric_rate, drop_ratio, max_hang, min_ovlp, cov); delete_useless_nodes(ug); renew_utg(ug, read_g, new_rtg_edges); if(i_cov && collect_p_trans == 0) goto skip_purge; @@ -24985,7 +24573,7 @@ R_to_U* ruIndex, int max_hang, int min_ovlp) if(bubble_dist > 0) { - asg_pop_bubble_primary_trio(ug, &bubble_dist, (uint32_t)-1, DROP, NULL, 0); + asg_pop_bubble_primary_trio(ug, &bubble_dist, (uint32_t)-1, DROP, NULL, NULL, 0); delete_useless_nodes(&ug); renew_utg(&ug, sg, &new_rtg_edges); } @@ -30561,7 +30149,7 @@ void flat_bubbles(asg_t *sg, uint8_t* r_het) if(bs_flag[v] == 1) continue; if(bs_flag[v] == 0) bs_flag[v] = 1; - if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0)) + if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) { //note b.b include end, does not include beg for (i = path = 0; i < b.b.n; i++) @@ -30650,7 +30238,7 @@ void flat_bubbles(asg_t *sg, uint8_t* r_het) if(is_het_b > path && is_het_s > path && (is_het_b+is_het_s)>(path<<2)) { - asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0); + asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0, NULL); n_pop++; } } @@ -30741,7 +30329,7 @@ void flat_bubbles_advance(asg_t *sg, ma_hit_t_alloc* sources, R_to_U* ruIndex, u if(bs_flag[v] == 1) continue; if(bs_flag[v] == 0) bs_flag[v] = 1; - if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0)) + if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) { bs_flag[v] = 2; bs_flag[b.S.a[0]^1] = 2; //note b.b include end, does not include beg @@ -30794,7 +30382,7 @@ void flat_bubbles_advance(asg_t *sg, ma_hit_t_alloc* sources, R_to_U* ruIndex, u if((C_bases/R_bases) <= het_thres) { // fprintf(stderr, "s-utg%.6ul\te-utg%.6ul\tC_bases:%lu\tR_bases:%lu\n", (v>>1)+1, (b.S.a[0]>>1)+1, C_bases, R_bases); - asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0); + asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0, NULL); n_pop++; } diff --git a/Overlaps.h b/Overlaps.h index 8a661cd..1a46c5b 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -120,6 +120,7 @@ typedef struct { uint32_t occ; double nw; uint8_t f:6, rev:1, del:1; + uint8_t qo:4, to:4; } u_trans_t; typedef struct { @@ -401,7 +402,6 @@ typedef struct { kvec_t(uint32_t) e; // visited edges/arcs } buf_t; - typedef struct { kvec_t(uint64_t) Nodes; kvec_t(uint64_t) Edges; @@ -858,13 +858,44 @@ typedef struct{ hc_edge *a; }hc_edge_warp; +typedef struct { + uint32_t qs, qe, qn, qus, que; + uint32_t ts, te, tn, tus, tue; +} utg_thit_t; + +typedef struct { + size_t n, m; + utg_thit_t* a; +} kv_utg_thit_t_t; + +typedef struct { + ma_hit_t_alloc* reverse_sources; + ma_sub_t *coverage_cut; + R_to_U* ruIndex; + asg_t *read_g; + kvec_asg_arc_t_offset u_buffer; + kvec_t_i32_warp tailIndex; + kvec_t_i32_warp prevIndex; + kv_utg_thit_t_t k_t_b; + kv_ca_buf_t c_buf; + kv_u_trans_t k_trans; + uint64_t *pos_idx, rn; + kvec_t(uint32_t) topo_res; + ma_ug_t *cug; + int max_hang; + int min_ovlp; + + ma_utg_v u; + kv_u_trans_t t; + buf_t b0, b1; +} utg_trans_t; + void init_hc_links(hc_links* link, uint64_t ug_num, trans_chain* t_ch); void destory_hc_links(hc_links* link); uint64_t get_bub_pop_max_dist(asg_t *g, buf_t *b); uint64_t get_bub_pop_max_dist_advance(asg_t *g, buf_t *b); -int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov, uint32_t is_update_chain); uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, uint64_t max_dist, buf_t *b, uint32_t positive_flag, -uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes, hap_cov_t *cov, uint32_t is_update_chain, uint32_t keep_d); +uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes, hap_cov_t *cov, uint32_t is_update_chain, uint32_t keep_d, utg_trans_t *o); void adjust_utg_by_primary(ma_ug_t **ug, asg_t* read_g, float drop_rate, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, @@ -951,6 +982,13 @@ typedef struct {///[cBeg, cEnd) void reset_u_trans_hit_idx(u_trans_hit_idx *t, uint32_t* i_x_a, uint32_t i_x_n, ma_ug_t *i_ug, asg_t *i_read_sg, trans_chain* i_t_ch, uint32_t i_cBeg, uint32_t i_cEnd); uint32_t get_u_trans_hit(u_trans_hit_idx *t, u_trans_hit_t *hit); +inline uint32_t get_offset_adjust(uint32_t offset, uint32_t offsetLen, uint32_t targetLen) +{ + return ((double)(offset)/(double)(offsetLen))*targetLen; +} + +uint32_t set_utg_offset(uint32_t *a, uint32_t a_n, ma_ug_t *ug, asg_t *read_sg, uint64_t* pos_idx, uint32_t is_clear, +uint32_t only_len); #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index a96e16c..ba6d81d 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -9,6 +9,7 @@ #include "kdq.h" #include "hic.h" #include "rcut.h" +#include "tovlp.h" KDQ_INIT(uint64_t) KSORT_INIT_GENERIC(uint64_t) @@ -1689,15 +1690,16 @@ uint32_t v_in_pos, uint32_t w_in_pos, uint32_t xUnitigLen, uint32_t yUnitigLen, return tmp; } -void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd) +void chain_trans_ovlp(hap_cov_t *cover, utg_trans_t *o, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd) { - ma_hit_t_alloc* reverse_sources = cov->reverse_sources; - ma_sub_t *coverage_cut = cov->coverage_cut; - int max_hang = cov->max_hang; - int min_ovlp = cov->min_ovlp; - kvec_asg_arc_t_offset* u_buffer = &(cov->u_buffer); - kvec_t_i32_warp* tailIndex = &(cov->tailIndex); - kvec_t_i32_warp* prevIndex = &(cov->prevIndex); + ma_hit_t_alloc* reverse_sources = (o? o->reverse_sources:cover->reverse_sources); + ma_sub_t *coverage_cut = (o? o->coverage_cut:cover->coverage_cut); + int max_hang = (o? o->max_hang:cover->max_hang); + int min_ovlp = (o? o->min_ovlp:cover->min_ovlp); + kvec_asg_arc_t_offset* u_buffer = (o? &(o->u_buffer):&(cover->u_buffer)); + kvec_t_i32_warp* tailIndex = (o? &(o->tailIndex):&(cover->tailIndex)); + kvec_t_i32_warp* prevIndex = (o? &(o->prevIndex):&(cover->prevIndex)); + uint64_t *pos_idx = (o? o->pos_idx:cover->pos_idx); ma_hit_t_alloc *xR = NULL; ma_hit_t *h = NULL; ma_sub_t *sq = NULL, *st = NULL; @@ -1804,11 +1806,11 @@ void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads rId = t.v>>1; if(read_sg->seq[rId].del == 1) continue; - if(cov->pos_idx[rId] == (uint64_t)-1) continue; - w = (uint32_t)(cov->pos_idx[rId]); + if(pos_idx[rId] == (uint64_t)-1) continue; + w = (uint32_t)(pos_idx[rId]); if(rId != (w>>1)) continue; - tmp = get_xy_pos_by_pos(read_sg, &t, v, w, len, cov->pos_idx[w>>1]>>32, + tmp = get_xy_pos_by_pos(read_sg, &t, v, w, len, pos_idx[w>>1]>>32, (uint32_t)-1, targetBaseLen, &(t.el)); if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue; if(t.el) continue; ///must @@ -1905,11 +1907,11 @@ void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads rId = t.v>>1; if(read_sg->seq[rId].del == 1) continue; - if(cov->pos_idx[rId] == (uint64_t)-1) continue; - w = (uint32_t)(cov->pos_idx[rId]); + if(pos_idx[rId] == (uint64_t)-1) continue; + w = (uint32_t)(pos_idx[rId]); if(rId != (w>>1)) continue; - tmp = get_xy_pos_by_pos(read_sg, &t, v, w, len, cov->pos_idx[w>>1]>>32, + tmp = get_xy_pos_by_pos(read_sg, &t, v, w, len, pos_idx[w>>1]>>32, (uint32_t)-1, targetBaseLen, &(t.el)); if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue; if(t.el) continue; ///must diff --git a/Purge_Dups.h b/Purge_Dups.h index 5c4f273..55107c4 100644 --- a/Purge_Dups.h +++ b/Purge_Dups.h @@ -75,7 +75,7 @@ void enable_debug_mode(uint32_t mode); hap_cov_t* init_hap_cov_t(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* sources, R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, int max_hang, int min_ovlp, uint32_t is_collect_trans); void destory_hap_cov_t(hap_cov_t **x); -void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd); +void chain_trans_ovlp(hap_cov_t *cov, utg_trans_t *o, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd); int get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn); void set_reverse_hap_overlap(hap_overlaps* dest, hap_overlaps* source, uint32_t* types); void print_hap_paf(ma_ug_t *ug, hap_overlaps* ovlp); diff --git a/hic.cpp b/hic.cpp index 02821f7..88d2f14 100644 --- a/hic.cpp +++ b/hic.cpp @@ -2317,7 +2317,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t if(ug->g->seq[v>>1].del) continue; if(asg_arc_n(ug->g, v) < 2) continue; if((bub->index[v]&(uint32_t)3) != 0) continue; - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0)) + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) { //beg is v, end is b.S.a[0] //note b.b include end, does not include beg @@ -2338,7 +2338,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t for (v = 0; v < n_vtx; ++v) { if((bub->index[v]&(uint32_t)3) !=2) continue; - if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL, 0, 0)) + if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL, 0, 0, NULL)) { //note b.b include end, does not include beg i = b.b.n + 1; @@ -2370,7 +2370,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t if((bub->num.a[k]>>31) == 0) bub->s_bub++; v = (bub->num.a[k]<<1)>>1; bub->num.a[k] = bub->list.n; - if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL, 0, 0)) + if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL, 0, 0, NULL)) { kv_push(uint64_t, bub->pathLen, pathLen); //beg is v, end is b.S.a[0] @@ -2948,7 +2948,7 @@ uint32_t get_specific_shortest_path(pdq_spec *p) asg_arc_t *av = NULL; if(p->v == (uint64_t)-1) { - reset_pdq(p->pq); + // reset_pdq(p->pq); p->pq->dis.a[p->src] = 0; if(p->pre) p->pre[p->src] = (uint32_t)-1; push_pdq(p->pq, p->src, 0); @@ -2984,17 +2984,40 @@ uint32_t get_specific_shortest_path(pdq_spec *p) return 0; } - -uint32_t check_trans_relation_by_path(uint32_t v, uint32_t w, pdq* pqv, pdq* pqw, -asg_t *sg, uint8_t *dest, uint8_t df, uint32_t df_occ, uint32_t* pre, double rate) +void get_utg_path(uint64_t s, uint64_t e, uint32_t *path, buf_t *res) +{ + uint64_t p = path[e], i; + res->b.n = 0; + kv_push(uint32_t, res->b, e); + while (p != s) + { + kv_push(uint32_t, res->b, p); + p = path[p]; + } + kv_push(uint32_t, res->b, s); + for (i = 0; i < (res->b.n>>1); i++) + { + p = res->b.a[i]; + res->b.a[i] = res->b.a[res->b.n-i-1]; + res->b.a[res->b.n-i-1] = p; + } + res->b.n--; +} +uint32_t check_trans_relation_by_path(uint32_t v, uint32_t w, pdq* pqv, uint32_t* path_v, buf_t *resv, +pdq* pqw, uint32_t* path_w, buf_t *resw, asg_t *sg, uint8_t *dest, uint8_t df, uint32_t df_occ, +double rate, long long *dis) { uint64_t r1, r2, d1, d2; pdq_spec p1, p2, *p = NULL; p1.v = (uint64_t)-1; p1.src = v; p1.pq = pqv; p1.sg = sg; p1.df_occ = df_occ; - p1.occ = df_occ; p1.flag = df; p1.dest = dest; p1.pre = pre; + p1.occ = df_occ; p1.flag = df; p1.dest = dest; p1.pre = path_v; if(resv) resv->b.n = 0; + reset_pdq(p1.pq); p2.v = (uint64_t)-1; p2.src = w; p2.pq = pqw; p2.sg = sg; p2.df_occ = df_occ; - p2.occ = df_occ; p2.flag = df; p2.dest = dest; p2.pre = pre; + p2.occ = df_occ; p2.flag = df; p2.dest = dest; p2.pre = path_w; if(resw) resw->b.n = 0; + reset_pdq(p2.pq); + + if(dis) (*dis) = 0; while (1) { @@ -3007,11 +3030,22 @@ asg_t *sg, uint8_t *dest, uint8_t df, uint32_t df_occ, uint32_t* pre, double rat { if(d2 <= d1*(1+rate) && d2 >= d1*(1-rate)) { + if(resv) get_utg_path(p1.src, p->v, p1.pre, resv); + if(resw) get_utg_path(p2.src, p->v, p2.pre, resw); + if(dis) + { + (*dis) = d1 + d2; + (*dis) -= (sg->seq[p1.src>>1].len + sg->seq[p2.src>>1].len); + } + + // fprintf(stderr, "p-utg%.6ul(d1-%lu), a-utg%.6ul(d2-%lu), conver-utg%.6lul\n", + // (v>>1) + 1, d1, (w>>1) + 1, d2, (p->v>>1) + 1); + + return 1; } } } - p = &p2; r2 = get_specific_shortest_path(p); if(r2 && p1.pq->vis.a[p->v] && p2.pq->vis.a[p->v]) @@ -3021,11 +3055,19 @@ asg_t *sg, uint8_t *dest, uint8_t df, uint32_t df_occ, uint32_t* pre, double rat { if(d2 <= d1*(1+rate) && d2 >= d1*(1-rate)) { + if(resv) get_utg_path(p1.src, p->v, p1.pre, resv); + if(resw) get_utg_path(p2.src, p->v, p2.pre, resw); + if(dis) + { + (*dis) = d1 + d2; + (*dis) -= (sg->seq[p1.src>>1].len + sg->seq[p2.src>>1].len); + } + // fprintf(stderr, "p-utg%.6ul(d1-%lu), a-utg%.6ul(d2-%lu), conver-utg%.6lul\n", + // (v>>1) + 1, d1, (w>>1) + 1, d2, (p->v>>1) + 1); return 1; } } } - if(!r1 && !r2) return 0; } diff --git a/hic.h b/hic.h index 31f3d3d..f78b590 100644 --- a/hic.h +++ b/hic.h @@ -101,7 +101,8 @@ void des_ug_idx(); uint64_t count_unique_k_mers(char *r, uint64_t len, uint64_t query, uint64_t target, uint64_t *all, uint64_t *found); void init_pdq(pdq* q, uint64_t utg_num); void destory_pdq(pdq* q); -uint32_t check_trans_relation_by_path(uint32_t v, uint32_t w, pdq* pqv, pdq* pqw, -asg_t *sg, uint8_t *dest, uint8_t df, uint32_t df_occ, uint32_t* pre, double rate); +uint32_t check_trans_relation_by_path(uint32_t v, uint32_t w, pdq* pqv, uint32_t* path_v, buf_t *resv, +pdq* pqw, uint32_t* path_w, buf_t *resw, asg_t *sg, uint8_t *dest, uint8_t df, uint32_t df_occ, double rate, +long long *dis); #endif diff --git a/tovlp.cpp b/tovlp.cpp new file mode 100644 index 0000000..2ba1e69 --- /dev/null +++ b/tovlp.cpp @@ -0,0 +1,1441 @@ +#define __STDC_LIMIT_MACROS +#include "float.h" +#include "horder.h" +#include +#include "hic.h" +#include "htab.h" +#include "assert.h" +#include "Overlaps.h" +#include "Hash_Table.h" +#include "Correct.h" +#include "Purge_Dups.h" +#include "rcut.h" +#include "khashl.h" +#include "kthread.h" +#include "ksort.h" +#include "kseq.h" // FASTA/Q parser +#include "kdq.h" +#include "tovlp.h" +#include "hic.h" +KSEQ_INIT(gzFile, gzread) +KDQ_INIT(uint64_t) + + +typedef struct { + kvec_t_u32_warp tt; + kvec_t_u32_warp stack; + uint8_t *vis; + pdq pq_p; + pdq pq_a; + buf_t b; + uint32_t min_v, offset1, offset2; + long long min_d; + int is_b; +} clean_t; + +typedef struct { + clean_t *a; + size_t n, m; + uint64_t max_dist; + uint8_t *bs_flag, fp, fa; + asg_t *g; + ma_ug_t *ug; +} clean_mul_t; + +clean_mul_t *init_clean_mul_t(asg_t *g, ma_ug_t *ug, uint64_t n_threads, uint8_t fp, uint8_t fa) +{ + uint32_t i; + clean_t *kt = NULL; + clean_mul_t *p = NULL; CALLOC(p, 1); + p->g = g; p->ug = ug; + CALLOC(p->bs_flag, p->g->n_seq<<1); + p->n = p->m = n_threads; + CALLOC(p->a, p->n); + for (i = 0, kt = NULL; i < p->n; i++) + { + kt = &(p->a[i]); + kv_init(kt->tt.a); kv_init(kt->stack.a); + CALLOC(kt->vis, p->g->n_seq<<1); + init_pdq(&kt->pq_p, p->g->n_seq<<1); + init_pdq(&kt->pq_a, p->g->n_seq<<1); + memset(&kt->b, 0, sizeof(kt->b)); + kt->b.a = (binfo_t*)calloc(p->g->n_seq<<1, sizeof(binfo_t)); + } + p->max_dist = get_bub_pop_max_dist_advance(g, &kt->b); + p->fp = fp; p->fa = fa; + return p; +} + +void destroy_clean_mul_t(clean_mul_t **p) +{ + uint32_t i; + clean_t *kt = NULL; + free((*p)->bs_flag); + for (i = 0; i < (*p)->n; i++) + { + kt = &((*p)->a[i]); + kv_destroy(kt->tt.a); + kv_destroy(kt->stack.a); + free(kt->vis); + destory_pdq(&kt->pq_p); + destory_pdq(&kt->pq_a); + free(kt->b.a); free(kt->b.S.a); free(kt->b.T.a); free(kt->b.b.a); free(kt->b.e.a); + } + free((*p)->a); + free((*p)); +} + +utg_trans_t *init_utg_trans_t(ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, R_to_U* ruIndex, asg_t *read_g, int max_hang, int min_ovlp) +{ + utg_trans_t *o = NULL; + CALLOC(o, 1); + o->reverse_sources = reverse_sources; + o->coverage_cut = coverage_cut; + o->ruIndex = ruIndex; + o->read_g = read_g; + o->rn = read_g->n_seq; + memset(&(o->b0), 0, sizeof(o->b0)); + memset(&(o->b1), 0, sizeof(o->b1)); + MALLOC(o->pos_idx, o->rn); memset(o->pos_idx, -1, o->rn*sizeof(uint64_t)); + o->cug = copy_untig_graph(ug); + o->max_hang = max_hang; + o->min_ovlp = min_ovlp; + return o; +} + +void destroy_utg_trans_t(utg_trans_t **o) +{ + ; +} + +uint64_t get_utg_chain_length(ma_ug_t *ug, uint32_t *a, uint32_t a_n) +{ + uint64_t i, k, l, qn, v, w, nv, offset; + asg_arc_t *av = NULL; + for (i = 0, offset = 0; i < a_n; i++) + { + qn = a[i]>>1; + l = ug->g->seq[qn].len; + if(i + 1 < a_n) + { + v = a[i]; w = a[i+1]; + av = asg_arc_a(ug->g, v); + nv = asg_arc_n(ug->g, v); + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + if(av[k].v == w) + { + l = asg_arc_len(av[k]); + break; + } + } + if(k >= nv) fprintf(stderr, "ERROR-mc\n"); + } + offset += l; + } + return offset; +} + +void refine_utg_thit_t(utg_thit_t *q, kv_ca_buf_t* cb) +{ + ///already know [qScur, qEcur), [qSpre, qEpre) + uint32_t s, e, i, si, ei; + s = q->qs; e = q->qe;///[s, e) + for (i = 0, si = ei = cb->n; i < cb->n; i++) + { + if(cb->a[i].c_x_p > s && si == cb->n) + { + si = i; + } + + if(cb->a[i].c_x_p > (e-1) && ei == cb->n) + { + ei = i; + } + + if(si != cb->n && ei != cb->n) break; + } + + if(si == 0 || ei == 0) fprintf(stderr, "ERROR-si-ei-0\n"); + if(si >= cb->n || ei >= cb->n) fprintf(stderr, "ERROR-si-ei-1\n"); + si--; ei--; + + if(s < cb->a[si].c_x_p || ((si + 1) < cb->n && s >= cb->a[si + 1].c_x_p)) + { + fprintf(stderr, "ERROR3\n"); + } + if(e < cb->a[ei].c_x_p || ((ei + 1) < cb->n && e > cb->a[ei + 1].c_x_p)) + { + fprintf(stderr, "ERROR4\n"); + } + ///si and ei must be less than (cb->n-1) + q->ts = cb->a[si].c_y_p + + get_offset_adjust(s-cb->a[si].c_x_p, cb->a[si+1].c_x_p-cb->a[si].c_x_p, cb->a[si+1].c_y_p-cb->a[si].c_y_p); + + q->te = cb->a[ei].c_y_p + + get_offset_adjust(e-cb->a[ei].c_x_p, cb->a[ei+1].c_x_p-cb->a[ei].c_x_p, cb->a[ei+1].c_y_p-cb->a[ei].c_y_p); + + ///might be equal + if(q->ts > q->te) fprintf(stderr, "ERROR5\n"); + // if(q->tScur >= q->tEcur) + // { + // fprintf(stderr, "\n###q->tScur: %u, s: %u, si: %u, q->tEcur: %u, e: %u, ei: %u\n", + // q->tScur, s, si, q->tEcur, e, ei); + + // fprintf(stderr, "###cb->a[si].c_x_p: %u, cb->a[si].c_y_p: %u, cb->a[si+1].c_x_p: %u, cb->a[si+1].c_y_p: %u\n", + // cb->a[si].c_x_p, cb->a[si].c_y_p, cb->a[si+1].c_x_p, cb->a[si+1].c_y_p); + + // fprintf(stderr, "###cb->a[ei].c_x_p: %u, cb->a[ei].c_y_p: %u, cb->a[ei+1].c_x_p: %u, cb->a[ei+1].c_y_p: %u\n", + // cb->a[ei].c_x_p, cb->a[ei].c_y_p, cb->a[ei+1].c_x_p, cb->a[ei+1].c_y_p); + + // fprintf(stderr, "ERROR5\n"); + // } +} + +///[ts, te) +void extract_sub_overlaps_utg_thit(uint32_t i_ts, uint32_t i_te, uint32_t i_tus, uint32_t i_tue, +uint32_t tn, kv_utg_thit_t_t* ktb, uint32_t bn) +{ + uint32_t i, ovlp, found, beg, end, offS, offE; + utg_thit_t *q = NULL, x; + for (i = found = 0; i < bn; i++) + { + q = &(ktb->a[i]);///for q, already know [qs, qe), [qus, que), [ts, te) + + ovlp = ((MIN(i_te, q->te) > MAX(i_ts, q->ts))? + MIN(i_te, q->te) - MAX(i_ts, q->ts):0); + if(found == 1 && ovlp == 0) break; + if(ovlp > 0) found = 1; + if(ovlp == 0) continue; + + + beg = MAX(i_ts, q->ts); end = MIN(i_te, q->te); + offS = beg - q->ts; offE = q->te - end; + x.ts = q->ts + offS; + x.te = q->te - offE; + //x.qs = q->qs + offS; + x.qs = q->qs + get_offset_adjust(offS, q->te-q->ts, q->qe-q->qs); + ///x.qe = q->qe - offE; + x.qe = q->qe - get_offset_adjust(offE, q->te-q->ts, q->qe-q->qs); + + x.qn = q->qn; + offS = beg - q->ts; offE = q->te - end; + if((x.qn&1) == 0) + { + // x.qus = q->qus + offS; + x.qus = q->qus + get_offset_adjust(offS, q->te-q->ts, q->que-q->qus); + // x.que = q->que - offE; + x.que = q->que - get_offset_adjust(offE, q->te-q->ts, q->que-q->qus); + } + else + { + // x.qus = q->qus + offE; + x.qus = q->qus + get_offset_adjust(offE, q->te-q->ts, q->que-q->qus); + // x.que = q->que - offS; + x.que = q->que - get_offset_adjust(offS, q->te-q->ts, q->que-q->qus); + } + + x.tn = tn; + offS = beg - i_ts; offE = i_te - end; + if((x.tn&1) == 0) + { + // x.tus = i_tus + offS; + x.tus = i_tus + get_offset_adjust(offS, i_te-i_ts, i_tue-i_tus); + // x.tue = i_tue - offE; + x.tue = i_tue - get_offset_adjust(offE, i_te-i_ts, i_tue-i_tus); + } + else + { + // x.tus = i_tus + offE; + x.tus = i_tus + get_offset_adjust(offE, i_te-i_ts, i_tue-i_tus); + // x.tue = i_tue - offS; + x.tue = i_tue - get_offset_adjust(offS, i_te-i_ts, i_tue-i_tus); + } + + kv_push(utg_thit_t, *ktb, x); + + // if(x.tus >= x.tue || x.qus >= x.que) + // { + // fprintf(stderr, "\n*********x.qn: %u, x.tn: %u\n", x.qn, x.tn); + // fprintf(stderr, "x.qus: %u, x.que: %u, x.tus: %u, x.tue: %u\n", + // x.qus, x.que, x.tus, x.tue); + // fprintf(stderr, "q->qs: %u, q->qe: %u, q->qus: %u, q->que: %u\n", + // q->qs, q->qe, q->qus, q->que); + // fprintf(stderr, "q->ts: %u, q->te: %u, q->tus: %u, q->tue: %u\n", + // q->ts, q->te, q->tus, q->tue); + // fprintf(stderr, "i_ts: %u, i_te: %u, i_tus: %u, i_tue: %u\n", + // i_ts, i_te, i_tus, i_tue); + // } + } +} +typedef struct {///[cBeg, cEnd) + uint32_t ui, len, cBeg, cEnd; + uint32_t *a, an; + ma_ug_t *ug; + utg_trans_t *o; +} utg_trans_hit_idx; + +uint32_t get_utg_trans_hit(utg_trans_hit_idx *t, utg_thit_t *hit) +{ + uint32_t r_beg, r_end, ovlp; + uint32_t k, l, qn, v, w, nv; + uint32_t *a = t->a, a_n = t->an; + asg_arc_t *av = NULL; + hit->qs = hit->qe = hit->qn = hit->qus = hit->que = (uint32_t)-1; + hit->ts = hit->te = hit->tn = hit->tus = hit->tue = (uint32_t)-1; + + while (t->ui < a_n) + { + qn = a[t->ui]; + l = t->ug->g->seq[qn>>1].len; + if(t->ui + 1 < a_n) + { + v = a[t->ui]; w = a[t->ui+1]; + av = asg_arc_a(t->ug->g, v); + nv = asg_arc_n(t->ug->g, v); + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + if(av[k].v == w) + { + l = asg_arc_len(av[k]); + break; + } + } + if(k >= nv) fprintf(stderr, "ERROR-mc-1\n"); + } + + r_beg = t->len; r_end = t->len + t->ug->g->seq[qn>>1].len; + t->len += l; t->ui++; + + ovlp = ((MIN(t->cEnd, r_end) > MAX(t->cBeg, r_beg))? (MIN(t->cEnd, r_end) - MAX(t->cBeg, r_beg)) : 0); + if(ovlp == 0 && r_beg >= t->cEnd) return 0; + if(ovlp == 0) continue; + + hit->qn = qn; hit->qs = r_beg; hit->qe = r_end; hit->qus = 0; hit->que = t->ug->g->seq[qn>>1].len; + if(hit->qs < t->cBeg) + { + hit->qus += (t->cBeg - hit->qs); + hit->qs = t->cBeg; + } + + if(hit->qe > t->cEnd) + { + hit->que -= (hit->qe - t->cEnd); + hit->qe = t->cEnd; + } + + return 1; + } + return 0; +} + +void reset_utg_trans_hit_idx(utg_trans_hit_idx *t, uint32_t* i_x_a, uint32_t i_x_n, ma_ug_t *i_ug, +utg_trans_t *i_o, uint32_t i_cBeg, uint32_t i_cEnd) +{ + t->a = i_x_a; + t->an = i_x_n; + t->ug = i_ug; + t->o = i_o; + t->cBeg = i_cBeg; + t->cEnd = i_cEnd; + t->ui = t->len = 0; +} + +void chain_trans_d(utg_trans_t *o, +uint32_t *pri_a, uint32_t pri_n, uint32_t pri_beg, uint64_t *i_pri_len, +uint32_t *aux_a, uint32_t aux_n, uint32_t aux_beg, uint64_t *i_aux_len, +ma_ug_t *ug, uint32_t flag, double score, const char* cmd) +{ + uint32_t i, len, bn; + uint64_t pri_len, aux_len; + kvec_asg_arc_t_offset* u_buffer = &(o->u_buffer); + kvec_t_i32_warp* tailIndex = &(o->tailIndex); + asg_arc_t_offset *tt = NULL; + ca_buf_t *tx = NULL; + + if(i_pri_len) pri_len = (*i_pri_len); + else pri_len = get_utg_chain_length(ug, pri_a, pri_n); + + if(i_aux_len) aux_len = (*i_aux_len); + else aux_len = get_utg_chain_length(ug, aux_a, aux_n); + + tt = NULL; o->c_buf.n = 0; o->k_t_b.n = 0; + if(tailIndex->a.n > 0) tt = &(u_buffer->a.a[tailIndex->a.a[0]]); + if(!tt || (pri_beg < (tt->Off>>32) && aux_beg < ((uint32_t)tt->Off))) + { + kv_pushp(ca_buf_t, o->c_buf, &tx); + tx->c_x_p = pri_beg; + tx->c_y_p = aux_beg; + + for (i = 0; i < tailIndex->a.n; i++) + { + tt = &(u_buffer->a.a[tailIndex->a.a[i]]); + kv_pushp(ca_buf_t, o->c_buf, &tx); + + tx->c_x_p = tt->Off>>32; + tx->c_y_p = (uint32_t)tt->Off; + } + } + else if(tailIndex->a.n == 1)//1 ele in chain + { + kv_pushp(ca_buf_t, o->c_buf, &tx); + tx->c_x_p = pri_beg; tx->c_y_p = aux_beg; + } + else if(tailIndex->a.n > 0) + { + uint32_t cx, cy, ax, ay, found = 0; + for (i = 0; i < tailIndex->a.n; i++) + { + cx = cy = ax = ay = (uint32_t)-1; + + cx = u_buffer->a.a[tailIndex->a.a[i]].Off>>32; + cy = (uint32_t)u_buffer->a.a[tailIndex->a.a[i]].Off; + if((i + 1) < tailIndex->a.n) + { + ax = u_buffer->a.a[tailIndex->a.a[i+1]].Off>>32; + ay = (uint32_t)u_buffer->a.a[tailIndex->a.a[i+1]].Off; + } + + kv_pushp(ca_buf_t, o->c_buf, &tx); + tx->c_x_p = cx; tx->c_y_p = cy; + + if(found) continue; + if(pri_beg > cx && pri_beg < ax && aux_beg > cy && aux_beg < ay) + { + kv_pushp(ca_buf_t, o->c_buf, &tx); + tx->c_x_p = pri_beg; tx->c_y_p = aux_beg; + found = 1; + } + } + } + + + // fprintf(stderr, "\ncmd-%s\n", cmd); + // fprintf(stderr, "pri_beg=%u, pri_len=%lu\n", pri_beg, pri_len); + // fprintf(stderr, "aux_beg=%u, aux_len=%lu\n", aux_beg, aux_len); + + // print_buf_t(ug, pri, "pri"); + // print_buf_t(ug, aux, "aux"); + + + tx = &(o->c_buf.a[o->c_buf.n-1]); + len = MIN(pri_len - tx->c_x_p, aux_len - tx->c_y_p); + if(len > 0)///insert boundary + { + kv_pushp(ca_buf_t, o->c_buf, &tx); + tx->c_x_p = o->c_buf.a[o->c_buf.n-2].c_x_p + len; + tx->c_y_p = o->c_buf.a[o->c_buf.n-2].c_y_p + len; + } + + tx = &(o->c_buf.a[0]);///insert boundary + if(tx->c_x_p != 0 && tx->c_y_p != 0)///already at boundary + { + len = MIN(tx->c_x_p, tx->c_y_p);///offset + kv_pushp(ca_buf_t, o->c_buf, &tx); + for (i = 0; (i + 1)< o->c_buf.n; i++) + { + o->c_buf.a[o->c_buf.n - i - 1] = o->c_buf.a[o->c_buf.n - i - 2]; + } + o->c_buf.a[0].c_x_p = o->c_buf.a[1].c_x_p - len; + o->c_buf.a[0].c_y_p = o->c_buf.a[1].c_y_p - len; + } + + ///chain is [s, e) + if(o->c_buf.a[0].c_x_p != 0 && o->c_buf.a[0].c_y_p != 0) fprintf(stderr, "ERROR1\n"); + if(o->c_buf.a[o->c_buf.n-1].c_x_p!= pri_len && + o->c_buf.a[o->c_buf.n-1].c_y_p!= aux_len) + { + fprintf(stderr, "ERROR2\n"); + } + + utg_trans_hit_idx iter; + utg_thit_t hit, *kh = NULL; + ////////prx + reset_utg_trans_hit_idx(&iter, pri_a, pri_n, ug, o, o->c_buf.a[0].c_x_p, o->c_buf.a[o->c_buf.n-1].c_x_p); + while(get_utg_trans_hit(&iter, &hit))//get [qs, qe), [qus, que) + { + refine_utg_thit_t(&hit, &(o->c_buf)); ///get [ts, te) + kv_push(utg_thit_t, o->k_t_b, hit); + } + bn = o->k_t_b.n; + + ////////aux + reset_utg_trans_hit_idx(&iter, aux_a, aux_n, ug, o, o->c_buf.a[0].c_y_p, o->c_buf.a[o->c_buf.n-1].c_y_p); + while(get_utg_trans_hit(&iter, &hit)) + { + extract_sub_overlaps_utg_thit(hit.qs, hit.qe, hit.qus, hit.que, hit.qn, &(o->k_t_b), bn); + } + + if(o->k_t_b.n - bn == 0) fprintf(stderr, "ERROR7\n"); + + // fprintf(stderr, "\n******o->k_t_b.n: %u, bn: %u\n", (uint32_t)o->k_t_b.n, bn); + // for (i = 0; i < pri_n; i++) + // { + // fprintf(stderr, "p-utg%.6ul\n", (pri_a[i]>>1) + 1); + // } + // for (i = 0; i < aux_n; i++) + // { + // fprintf(stderr, "a-utg%.6ul\n", (aux_a[i]>>1) + 1); + // } + + u_trans_t *kt = NULL; + double x_score, y_score; + for (i = bn; i < o->k_t_b.n; i++) + { + kh = &(o->k_t_b.a[i]); + if(kh->que <= kh->qus) continue; + if(kh->tue <= kh->tus) continue; + kv_pushp(u_trans_t, o->k_trans, &kt); + kt->f = flag; kt->rev = ((kh->qn ^ kh->tn) & 1); kt->del = 0; + kt->qn = kh->qn>>1; kt->qs = kh->qus; kt->qe = kh->que; kt->qo = kh->qn&1; + kt->tn = kh->tn>>1; kt->ts = kh->tus; kt->te = kh->tue; kt->to = kh->tn&1; + if(score < 0) + { + kt->nw = (MIN((kt->qe - kt->qs), (kt->te - kt->ts)))*CHAIN_MATCH; + } + else + { + x_score = ((double)(kt->qe-kt->qs)/(double)(o->c_buf.a[o->c_buf.n-1].c_x_p-o->c_buf.a[0].c_x_p))*score; + y_score = ((double)(kt->te-kt->ts)/(double)(o->c_buf.a[o->c_buf.n-1].c_y_p-o->c_buf.a[0].c_y_p))*score; + kt->nw = MIN(x_score, y_score); + } + } +} + + +void chain_bubble(buf_t *pri, uint64_t pri_len, buf_t* aux, uint64_t aux_len, +uint32_t beg, uint32_t sink, ma_ug_t *ug, utg_trans_t *o) +{ + if(pri->b.n == 0 || aux->b.n == 0) return; + asg_arc_t *av = NULL; + uint32_t pri_v, aux_v, nv, i, priBeg, priEnd, auxBeg, auxEnd; + + priBeg = priEnd = auxBeg = auxEnd = (uint32_t)-1; + pri_v = pri->b.a[0]; aux_v = aux->b.a[0]; + av = asg_arc_a(ug->g, beg); + nv = asg_arc_n(ug->g, beg); + for (i = 0; i < nv; ++i) + { + if(av[i].del) continue; + if(av[i].v == pri_v) priBeg = av[i].ol; + if(av[i].v == aux_v) auxBeg = av[i].ol; + } + + pri_v = pri->b.a[pri->b.n-1]^1; aux_v = aux->b.a[aux->b.n-1]^1; + av = asg_arc_a(ug->g, sink); + nv = asg_arc_n(ug->g, sink); + for (i = 0; i < nv; ++i) + { + if(av[i].del) continue; + if(av[i].v == pri_v) priEnd = ((pri_len > av[i].ol)? (pri_len - av[i].ol - 1) : 0); + if(av[i].v == aux_v) auxEnd = ((aux_len > av[i].ol)? (aux_len - av[i].ol - 1) : 0); + } + + if(priBeg == (uint32_t)-1 || priEnd == (uint32_t)-1 || auxBeg == (uint32_t)-1 || auxEnd == (uint32_t)-1) + { + fprintf(stderr, "ERROR-s_bubble\n"); + } + + o->u_buffer.a.n = o->tailIndex.a.n = 0; + + kv_resize(asg_arc_t_offset, o->u_buffer.a, 1); + o->u_buffer.a.n = 1; + o->u_buffer.a.a[0].Off = priEnd; + o->u_buffer.a.a[0].Off <<= 32; + o->u_buffer.a.a[0].Off |= auxEnd; + + kv_resize(int32_t, o->tailIndex.a, 1); + o->tailIndex.a.n = 1; + o->tailIndex.a.a[0] = 0; + + chain_trans_d(o, pri->b.a, pri->b.n, priBeg, &pri_len, aux->b.a, aux->b.n, auxBeg, &aux_len, ug, RC_0, -1024, __func__); +} + +void chain_c_bubble(uint32_t query, buf_t *target, buf_t *idx, ma_ug_t *ug, utg_trans_t *o) +{ + if(target->b.n == 0) return; + uint32_t qs, qe, ts, te, i, v, ovlp; + qs = idx->a[query].d; qe = qs + ug->g->seq[query>>1].len; + uint64_t qlen = ug->g->seq[query>>1].len, tlen; + o->u_buffer.a.n = o->tailIndex.a.n = 0; + + for (i = 0; i < target->b.n; ++i) + { + v = target->b.a[i]; + if(v < query) continue; //avoid dup + ts = idx->a[v].d; te = ts + ug->g->seq[v>>1].len; tlen = ug->g->seq[v>>1].len; + ovlp = ((MIN(qe, te) > MAX(qs, ts))? (MIN(qe, te) - MAX(qs, ts)) : 0); + if(ovlp == 0) continue; + chain_trans_d(o, &query, 1, MAX(qs, ts) - qs, &qlen, + &v, 1, MAX(qs, ts) - ts, &tlen, ug, RC_0, -1024, __func__); + } + // fprintf(stderr, "\nocc: %u\n", (uint32_t)(cov->t_ch->k_trans.n - i_n)); + // for (i = i_n; i < cov->t_ch->k_trans.n; i++) + // { + // fprintf(stderr, "s-utg%.6ul\t%u\t%u\td-utg%.6ul\t%u\t%u\trev(%u)\n", + // cov->t_ch->k_trans.a[i].qn+1, cov->t_ch->k_trans.a[i].qs, cov->t_ch->k_trans.a[i].qe, + // cov->t_ch->k_trans.a[i].tn+1, cov->t_ch->k_trans.a[i].ts, cov->t_ch->k_trans.a[i].te, + // cov->t_ch->k_trans.a[i].rev); + // } +} + +void tpSort(asg_t *g, utg_trans_t *o, uint32_t beg, uint32_t sink) +{ + buf_t *b = &(o->b0); + uint64_t *visited = o->pos_idx; + uint32_t v = beg, nv, kv, i; + asg_arc_t *av = NULL; + + b->b.n = 0; o->topo_res.n = 0; + kv_push(uint32_t, b->b, v); + while (b->b.n > 0) + { + ///b->b.n--; + v = b->b.a[b->b.n-1]; + if(visited[v>>1] == (uint64_t)-1) + { + visited[v>>1] = 0; + } + + nv = asg_arc_n(g, v); + av = asg_arc_a(g, v); + for (i = kv = 0; i < nv; i++) + { + if(av[i].del) continue; + if((av[i].v>>1) == (beg>>1) || (av[i].v>>1) == (sink>>1)) continue; + if(visited[av[i].v>>1] != (uint64_t)-1) continue; + kv_push(uint32_t, b->b, av[i].v); + kv++; + } + + if(kv != 0) continue; + b->b.n--; + if(visited[v>>1] != 1) + { + kv_push(uint32_t, o->topo_res, v); + visited[v>>1] = 1; + } + } + for (i = 0; i < o->topo_res.n; ++i) + { + visited[o->topo_res.a[i]>>1] = (uint64_t)-1; + } + + o->topo_res.n--;//remove beg + for (i = 0; i < (o->topo_res.n>>1); ++i) + { + v = o->topo_res.a[i]; + o->topo_res.a[i] = o->topo_res.a[o->topo_res.n - i - 1]; + o->topo_res.a[o->topo_res.n - i - 1] = v; + } + + /*******************************for debug************************************/ + // for (i = 0; i < cov->n; ++i) + // { + // if(visited[i] != (uint64_t)-1) fprintf(stderr, "ERROR-2\n"); + // } + // debug_topo_sorting(g, cov, beg, sink); + /*******************************for debug************************************/ +} + +void dfs_trans_chain(asg_t *g, utg_trans_t *o, uint32_t v, uint32_t beg, uint32_t sink) +{ + buf_t *b = &(o->b0); + b->b.n = 0; + if(v == beg || v == sink) return; + uint64_t *flag = o->pos_idx; + asg_arc_t *acur = NULL; + uint32_t cur, ncur, i; + v = v << 1; + kv_push(uint32_t, b->b, v); + while (b->b.n > 0) + { + b->b.n--; + cur = b->b.a[b->b.n]; + if(flag[cur>>1] == 0 && (cur>>1) != (v>>1)) continue; + flag[cur>>1] = 0; + + ncur = asg_arc_n(g, cur); + acur = asg_arc_a(g, cur); + for (i = 0; i < ncur; i++) + { + if(acur[i].del) continue; + if((acur[i].v>>1) == beg || (acur[i].v>>1) == sink) continue; + if(flag[acur[i].v>>1] == 0) continue; + kv_push(uint32_t, b->b, acur[i].v); + } + } + + v = v + 1; + kv_push(uint32_t, b->b, v); + while (b->b.n > 0) + { + b->b.n--; + cur = b->b.a[b->b.n]; + if(flag[cur>>1] == 0 && (cur>>1) != (v>>1)) continue; + flag[cur>>1] = 0; + + ncur = asg_arc_n(g, cur); + acur = asg_arc_a(g, cur); + for (i = 0; i < ncur; i++) + { + if(acur[i].del) continue; + if((acur[i].v>>1) == beg || (acur[i].v>>1) == sink) continue; + if(flag[acur[i].v>>1] == 0) continue; + kv_push(uint32_t, b->b, acur[i].v); + } + } + + b->b.n = 0; + for (i = 0; i < o->topo_res.n; ++i) //has been sorted + { + if(flag[o->topo_res.a[i]>>1] == 0) + { + flag[o->topo_res.a[i]>>1] = (uint64_t)-1; + } + else + { + kv_push(uint32_t, b->b, o->topo_res.a[i]); + } + } + + + /*******************************for debug************************************/ + // for (i = 0; i < cov->n; ++i) + // { + // if(flag[i] != (uint64_t)-1) fprintf(stderr, "ERROR-0\n"); + // } + /*******************************for debug************************************/ +} + +void asg_bub_collect_ovlp(ma_ug_t *ug, uint32_t v0, buf_t *b, utg_trans_t *o) +{ + uint32_t i, uId; + ///b->S.a[0] is the sink of this bubble + if(get_real_length(ug->g, v0, NULL) == 2 && get_real_length(ug->g, b->S.a[0]^1, NULL) == 2) + { + long long tmp, max_stop_nodeLen, max_stop_baseLen, bch_occ[2], len[2]; + uint32_t bch[2], convex[2]; + get_real_length(ug->g, v0, bch); + + ///in rare cases, one side of a bubble might be empty + if((bch[0]>>1)!=(b->S.a[0]>>1) && (bch[1]>>1)!=(b->S.a[0]>>1)) + { + get_unitig(ug->g, NULL, bch[0], &convex[0], &bch_occ[0], &tmp, + &max_stop_nodeLen, &max_stop_baseLen, 1, NULL); + get_unitig(ug->g, NULL, bch[1], &convex[1], &bch_occ[1], &tmp, + &max_stop_nodeLen, &max_stop_baseLen, 1, NULL); + if(((bch_occ[0] + bch_occ[1] + 1) == (uint32_t)b->b.n) && + get_real_length(ug->g, convex[0], NULL) == 1 && get_real_length(ug->g, convex[1], NULL) == 1) + { + get_real_length(ug->g, convex[0], &convex[0]); + get_real_length(ug->g, convex[1], &convex[1]); + if(convex[0] == b->S.a[0] && convex[1] == b->S.a[0]) + { + o->b0.b.n = 0; + get_unitig(ug->g, NULL, bch[0], &convex[0], &bch_occ[0], &len[0], &max_stop_nodeLen, &max_stop_baseLen, 1, &(o->b0)); + o->b1.b.n = 0; + get_unitig(ug->g, NULL, bch[1], &convex[1], &bch_occ[1], &len[1], &max_stop_nodeLen, &max_stop_baseLen, 1, &(o->b1)); + + chain_bubble(&(o->b0), len[0], &(o->b1), len[1], v0, b->S.a[0]^1, ug, o); + return; + } + } + } + } + + tpSort(ug->g, o, v0, b->S.a[0]); + ///if(o->topo_res.n != b->b.n - 1) fprintf(stderr, "ERROR-4\n"); + if(o->topo_res.n == 0) return; + for (i = 0; i < o->topo_res.n; ++i) + { + uId = o->topo_res.a[i]>>1; + if(uId == (b->S.a[0]>>1)) continue; + + dfs_trans_chain(ug->g, o, uId, v0>>1, b->S.a[0]>>1); + + if(o->b0.b.n == 0) continue; + chain_c_bubble(o->topo_res.a[i], &(o->b0), b, ug, o); + } +} + + +void collect_trans_ovlp(const char* cmd, buf_t* pri, uint64_t pri_offset, buf_t* aux, uint64_t aux_offset, +ma_ug_t *ug, utg_trans_t *o) +{ + uint32_t thre_pri; + uint64_t len_aux; + if(pri->b.n == 0 || aux->b.n == 0) return; + + len_aux = set_utg_offset(aux->b.a, aux->b.n, ug, o->read_g, o->pos_idx, 0, 0); + chain_trans_ovlp(NULL, o, ug, o->read_g, pri, len_aux, &thre_pri); + + + if(thre_pri > 0) + { + chain_trans_d(o, pri->b.a, pri->b.n, pri_offset, NULL, aux->b.a, aux->b.n, aux_offset, &len_aux, + ug, RC_1, -1024, __func__); + // 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"); + } + set_utg_offset(aux->b.a, aux->b.n, ug, o->read_g, o->pos_idx, 1, 0); +} + + +uint32_t dfs_set(asg_t *g, uint32_t v0, uint32_t fbv, kvec_t_u32_warp *stack, kvec_t_u32_warp *tmp, uint8_t *vis, uint8_t flag) +{ + uint32_t cur, ncur, i, occ = 0;; + asg_arc_t *acur = NULL; + stack->a.n = 0; + kv_push(uint32_t, stack->a, v0); + // vis[v0] = 0; + + while (stack->a.n > 0) + { + stack->a.n--; + cur = stack->a.a[stack->a.n]; + if(vis[cur] && (!(vis[cur]&flag))) + { + vis[cur] |= flag; + kv_push(uint32_t, tmp->a, cur); + occ++; + } + if(vis[cur]) continue; + else kv_push(uint32_t, tmp->a, cur); + + vis[cur] |= flag; + ncur = asg_arc_n(g, cur); + acur = asg_arc_a(g, cur); + for (i = 0; i < ncur; i++) + { + if(acur[i].del) continue; + if((acur[i].v>>1) == fbv) continue; + if(vis[acur[i].v] && (!(vis[acur[i].v]&flag))) + { + vis[acur[i].v] |= flag; + kv_push(uint32_t, tmp->a, acur[i].v); + occ++; + continue; + } + + kv_push(uint32_t, stack->a, acur[i].v); + } + } + + return occ; +} + + +uint32_t dfs_reach(asg_t *g, uint32_t src, uint32_t dest, kvec_t_u32_warp *stack, kvec_t_u32_warp *tmp, uint8_t *vis) +{ + uint32_t cur, ncur, i, occ = 0;; + asg_arc_t *acur = NULL; + stack->a.n = 0; tmp->a.n = 0; + kv_push(uint32_t, stack->a, src); + + while (stack->a.n > 0) + { + stack->a.n--; + cur = stack->a.a[stack->a.n]; + if(vis[cur]) continue; + else kv_push(uint32_t, tmp->a, cur); + vis[cur] = 1; + if(cur == dest) + { + occ = 1; + break; + } + + + ncur = asg_arc_n(g, cur); + acur = asg_arc_a(g, cur); + for (i = 0; i < ncur; i++) + { + if(acur[i].del) continue; + if(vis[acur[i].v]) continue; + kv_push(uint32_t, stack->a, acur[i].v); + if(acur[i].v == dest) + { + occ = 1; + break; + } + } + if(occ) break; + } + + for (i = 0; i < tmp->a.n; i++) + { + vis[tmp->a.a[i]] = 0; + } + + stack->a.n = 0; tmp->a.n = 0; + return occ; +} + +int is_hap_ovlp(ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, + uint32_t *a, uint32_t a_n, uint32_t *b, uint32_t b_n) +{ + uint32_t i, k, j, qn, tn, is_Unitig, uId, found = 0; + ma_utg_t *u = NULL; + for (i = 0; i < b_n; i++) + { + u = &(ug->u.a[b[i]]); + for (k = 0; k < u->n; k++) + { + qn = (u->a[k]>>33); + set_R_to_U(ruIndex, qn, b[i], 1, NULL); + } + } + + + for (i = 0; i < a_n; i++) + { + u = &(ug->u.a[a[i]]); + for (k = 0; k < u->n; k++) + { + qn = (u->a[k]>>33); + for (j = 0; j < reverse_sources[qn].length; j++) + { + tn = Get_tn(reverse_sources[qn].buffer[j]); + if(read_sg->seq[tn].del == 1) + { + get_R_to_U(ruIndex, tn, &tn, &is_Unitig); + if(tn == (uint32_t)-1 || is_Unitig == 1 || read_sg->seq[tn].del == 1) continue; + } + + get_R_to_U(ruIndex, tn, &uId, &is_Unitig); + if(uId!=(uint32_t)-1 && is_Unitig == 1) + { + found = 1; + break; + } + } + if(found) break; + } + if(found) break; + } + + + + + for (i = 0; i < b_n; i++) + { + u = &(ug->u.a[b[i]]); + for (k = 0; k < u->n; k++) + { + qn = (u->a[k]>>33); + ruIndex->index[qn] = (uint32_t)-1; + } + } + + return found; +} + +uint32_t get_dir_v(uint32_t x, buf_t *b) +{ + uint32_t i; + for (i = 0; i < b->b.n; i++) + { + if((b->b.a[i]>>1) == x) return b->b.a[i]; + } + + return (uint32_t)-1; +} + +int is_pop_unitig(uint32_t v, asg_t *g, ma_ug_t *ug, kvec_t_u32_warp *tt, kvec_t_u32_warp *stack, uint8_t *vis, +pdq *pq_p, pdq *pq_a, uint8_t fp, uint8_t fa, long long *d) +{ + asg_arc_t *as = NULL; + uint32_t s, sv, return_flag, convex, ns, i, k, nc, p_n, found, *a_a = NULL, a_n; + long long ll, tmp, max_stop_nodeLen, max_stop_baseLen; + if(g->seq[v>>1].del || g->seq[v>>1].c == ALTER_LABLE) return 0; + if(asg_arc_n(g, v) == 0 || get_real_length(g, v, NULL) != 1) return 0; + get_real_length(g, v, &s); + if(get_real_length(g, s^1, NULL) < 2) return 0; + return_flag = get_unitig(g, ug, v^1, &convex, &ll, &tmp, &max_stop_nodeLen,&max_stop_baseLen, 1, NULL); + if(return_flag == LOOP) return 0; + + as = asg_arc_a(g, convex); ns = asg_arc_n(g, convex); + for (i = nc = 0; i < ns; i++) + { + if(as[i].del) continue; + if(get_real_length(g, as[i].v^1, NULL) < 2) break; + nc++; + } + if(nc == 0 || i < ns) return 0; + + tt->a.n = 0; s^=1; sv = v^1; + dfs_set(g, sv, s>>1, stack, tt, vis, fp); + p_n = tt->a.n; + as = asg_arc_a(g, s); ns = asg_arc_n(g, s); found = 0; + for (i = 0; i < ns; i++) + { + if(as[i].del || as[i].v == sv) continue; + nc = dfs_set(g, as[i].v, s>>1, stack, tt, vis, fa); + if(nc && check_trans_relation_by_path(sv, as[i].v, pq_p, NULL, NULL, pq_a, NULL, NULL, g, + vis, fp+fa, nc, 0.45, d)) + { + found = 1; + } + a_a = tt->a.a + p_n; a_n = tt->a.n - p_n; + for (k = 0; k < a_n; k++) + { + if(vis[a_a[k]]&fa) vis[a_a[k]] -= fa; + } + tt->a.n = p_n; + if(found) break; + } + a_a = tt->a.a; a_n = p_n; + for (k = 0; k < a_n; k++) vis[a_a[k]] = 0; + return found; +} + +static void unitig_iden_worker(void *_data, long eid, int tid) +{ + long long d; + clean_mul_t *buf = (clean_mul_t*)_data; + clean_t *b = &(buf->a[tid]); + if(is_pop_unitig(eid, buf->g, buf->ug, &(b->tt), &(b->stack), b->vis, &(b->pq_p), &(b->pq_a), buf->fp, buf->fa, &d)) + { + if(b->min_v == (uint32_t)-1 || (b->min_d > d) || + (b->min_d == d && buf->g->seq[b->min_v>>1].len > buf->g->seq[eid>>1].len) || + (b->min_d == d && buf->g->seq[b->min_v>>1].len == buf->g->seq[eid>>1].len && b->min_v > eid)) + { + b->min_v = eid; b->min_d = d; + } + } +} + +int is_pop_bub(uint32_t v, asg_t *g, ma_ug_t *ug, buf_t *b, uint64_t max_dist, uint8_t *bs_flag) +{ + bs_flag[v] = 4; + asg_arc_t *av = NULL; + uint32_t nv, i, n_arc; + int occ = 0; + nv = asg_arc_n(g, v); + av = asg_arc_a(g, v); + ///some node could be deleted + if (nv < 2 || g->seq[v>>1].del || g->seq[v>>1].c == ALTER_LABLE) return occ; + ///some edges could be deleted + for (i = n_arc = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs + if (!av[i].del) ++n_arc; + if (n_arc < 2) return occ; + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, max_dist, b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) + { + //beg is v, end is b.S.a[0] + //note b.b include end, does not include beg + bs_flag[v] = 0; + occ = 1; + } + return occ; +} + +static void bub_iden_worker(void *_data, long eid, int tid) +{ + clean_mul_t *buf = (clean_mul_t*)_data; + clean_t *b = &(buf->a[tid]); + if(is_pop_bub(eid, buf->g, buf->ug, &(b->b), buf->max_dist, buf->bs_flag)) + { + b->is_b = 1; + } +} + +int asg_pop_bubble_primary_trio(asg_t *g, ma_ug_t *ug, uint8_t *bs_flag, buf_t *b, +uint64_t max_dist, uint32_t positive_flag, uint32_t negative_flag, utg_trans_t *o) +{ + uint32_t v, n_vtx = g->n_seq * 2, n_arc, nv, i; + uint64_t n_pop = 0; + asg_arc_t *av = NULL; + + if(max_dist > 0) + { + for (v = 0; v < n_vtx; ++v) + { + if(bs_flag[v] != 0) continue; + nv = asg_arc_n(g, v); + av = asg_arc_a(g, v); + ///some node could be deleted + if (nv < 2 || g->seq[v>>1].del || g->seq[v>>1].c == ALTER_LABLE) continue; + ///some edges could be deleted + for (i = n_arc = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs + if (!av[i].del) ++n_arc; + if (n_arc < 2) continue; + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, max_dist, b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) + { + //beg is v, end is b.S.a[0] + //note b.b include end, does not include beg + for (i = 0; i < b->b.n; i++) + { + if(b->b.a[i]==v || b->b.a[i]==b->S.a[0]) continue; + bs_flag[b->b.a[i]] = bs_flag[b->b.a[i]^1] = 1; + } + bs_flag[v] = 2; bs_flag[b->S.a[0]^1] = 3; + } + } + + //traverse all node with two directions + for (v = 0; v < n_vtx; ++v) { + if(bs_flag[v] !=2) continue; + nv = asg_arc_n(g, v); + av = asg_arc_a(g, v); + ///some node could be deleted + if (nv < 2 || g->seq[v>>1].del || g->seq[v>>1].c == ALTER_LABLE) continue; + ///some edges could be deleted + for (i = n_arc = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs + if (!av[i].del) ++n_arc; + if (n_arc > 1) + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, b, positive_flag, negative_flag, 1, NULL, NULL, NULL, 0, 0, o); + } + + if(VERBOSE >= 1) + { + fprintf(stderr, "[M::%s] popped %lu bubbles\n", __func__, (unsigned long)n_pop); + } + } + + if (n_pop) asg_cleanup(g); + return n_pop; +} + +int get_min_dec(clean_mul_t *cl, uint32_t positive_flag, uint32_t negative_flag, utg_trans_t *o, uint32_t *v) +{ + uint32_t i, n_pop; + clean_t *p = NULL; + (*v) = (uint32_t)-1; + for (i = 0; i < cl->n; i++) cl->a[i].is_b = 0, cl->a[i].min_v = (uint32_t)-1; + kt_for(cl->n, bub_iden_worker, cl, cl->g->n_seq<<1); + for (i = 0; i < cl->n; i++) + { + if(cl->a[i].is_b) ///pop bubble + { + n_pop = asg_pop_bubble_primary_trio(cl->g, cl->ug, cl->bs_flag, &(cl->a[i].b), + cl->max_dist, positive_flag, negative_flag, o); + if(n_pop ==0) fprintf(stderr, "ERROR-n_pop\n"); + break; + } + } + kt_for(cl->n, unitig_iden_worker, cl, cl->g->n_seq<<1); + for (i = 0; i < cl->n; i++) + { + if(cl->a[i].min_v == (uint32_t)-1) continue; + // if(b->min_v == (uint32_t)-1 || (b->min_d > d) || (b->min_d == d && b->min_v < eid)) + if(!p || p->min_d > cl->a[i].min_d || + (p->min_d == cl->a[i].min_d && cl->g->seq[p->min_v>>1].len > cl->g->seq[cl->a[i].min_v>>1].len) || + (p->min_d == cl->a[i].min_d && cl->g->seq[p->min_v>>1].len == cl->g->seq[cl->a[i].min_v>>1].len && p->min_v > cl->a[i].min_v)) + { + p = &(cl->a[i]); + } + } + if(p) (*v) = p->min_v; + return p?1:0; +} + +int asg_arc_decompress_mul(asg_t *g, ma_ug_t *ug, asg_t *read_sg, uint32_t positive_flag, uint32_t negative_flag, +ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, utg_trans_t *o) +{ + double startTime = Get_T(); + asg_arc_t *as = NULL, *pm = NULL; + u_trans_t *kt = NULL; + uint32_t v, s, sv, nc, ns, n_vtx = g->n_seq * 2, n_reduced = 0, convex, pi, qn, tn; + uint32_t return_flag, k, m, i, p_n, a_n, *a_a = NULL; + uint8_t fp = 1, fa = 2, found; + long long ll, tmp, max_stop_nodeLen, max_stop_baseLen; + kvec_t_u32_warp tt; kv_init(tt.a); + kvec_t_u32_warp stack; kv_init(stack.a); + uint8_t *vis = NULL; CALLOC(vis, n_vtx); + buf_t b; + memset(&b, 0, sizeof(buf_t)); + pdq pq_p, pq_a; + init_pdq(&pq_p, g->n_seq<<1); + init_pdq(&pq_a, g->n_seq<<1); + uint32_t *path_p = NULL, *path_q = NULL; + CALLOC(path_p, g->n_seq<<1); + CALLOC(path_q, g->n_seq<<1); + clean_mul_t *cl = init_clean_mul_t(g, ug, asm_opt.thread_num, fp, fa); + + while (get_min_dec(cl, positive_flag, negative_flag, o, &v)) + { + if(g->seq[v>>1].del || g->seq[v>>1].c == ALTER_LABLE) continue; + if(asg_arc_n(g, v) == 0 || get_real_length(g, v, NULL) != 1) continue; + get_real_length(g, v, &s); + if(get_real_length(g, s^1, NULL) < 2) continue; + return_flag = get_unitig(g, ug, v^1, &convex, &ll, &tmp, &max_stop_nodeLen,&max_stop_baseLen, 1, NULL); + if(return_flag == LOOP) continue; + + as = asg_arc_a(g, convex); ns = asg_arc_n(g, convex); + for (i = nc = 0; i < ns; i++) + { + if(as[i].del) continue; + if(get_real_length(g, as[i].v^1, NULL) < 2) break; + nc++; + } + if(nc == 0 || i < ns) continue; + + tt.a.n = 0; s^=1; sv = v^1; + dfs_set(g, sv, s>>1, &stack, &tt, vis, fp); + p_n = tt.a.n; + as = asg_arc_a(g, s); ns = asg_arc_n(g, s); found = 0; + for (i = 0, pm = NULL; i < ns; i++) + { + if(as[i].del || as[i].v != sv) continue; + pm = &(as[i]); + break; + } + + for (i = 0; i < ns; i++) + { + if(as[i].del || as[i].v == sv) continue; + nc = dfs_set(g, as[i].v, s>>1, &stack, &tt, vis, fa); + if(nc && check_trans_relation_by_path(sv, as[i].v, &pq_p, path_p, &(o->b0), &pq_a, path_q, &(o->b1), g, + vis, fp+fa, nc, 0.45, NULL)) + { + found = 1; + } + a_a = tt.a.a + p_n; a_n = tt.a.n - p_n; + for (k = 0; k < a_n; k++) + { + if(vis[a_a[k]]&fa) vis[a_a[k]] -= fa; + } + tt.a.n = p_n; + if(found) break; + } + a_a = tt.a.a; a_n = p_n; + for (k = 0; k < a_n; k++) vis[a_a[k]] = 0; + + if(found == 0) fprintf(stderr, "ERROR-found\n"); + if(found) + { + b.b.n = 0; + get_unitig(g, ug, sv, &convex, &ll, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 1, &b); + for (k = 0; k < b.b.n; k++) + { + g->seq[b.b.a[k]>>1].c = ALTER_LABLE; + } + + for (k = 0; k < b.b.n; k++) + { + asg_seq_drop(g, b.b.a[k]>>1); + } + + // uint32_t ui; + // fprintf(stderr, "\n"); + // for (ui = 0; ui < o->b0.b.n; ui++) + // { + // fprintf(stderr, "p-utg%.6ul\n", (o->b0.b.a[ui]>>1)+1); + // } + // for (ui = 0; ui < o->b1.b.n; ui++) + // { + // fprintf(stderr, "a-utg%.6ul\n", (o->b1.b.a[ui]>>1)+1); + // } + + + if(o->b0.b.n < b.b.n) fprintf(stderr, "ERROR-ll\n"); + o->b0.b.n = b.b.n; pi = o->k_trans.n; + collect_trans_ovlp(__func__, &(o->b0), pm->ol, &(o->b1), as[i].ol, ug, o); + + for (k = m = pi; k < o->k_trans.n; k++) + { + kt = &(o->k_trans.a[k]); + if(is_hap_ovlp(ug, read_sg, reverse_sources, ruIndex, &(kt->qn), 1, &(kt->tn), 1) || + is_hap_ovlp(ug, read_sg, reverse_sources, ruIndex, &(kt->tn), 1, &(kt->qn), 1)) + { + qn = get_dir_v(kt->qn, &(o->b0)); + tn = get_dir_v(kt->tn, &(o->b1)); + + if(dfs_reach(o->cug->g, qn, tn, &stack, &tt, vis) || dfs_reach(o->cug->g, tn, qn, &stack, &tt, vis)) + { + // fprintf(stderr, "<<<<<qn + 1, kt->qs, kt->qe, kt->tn + 1, kt->ts, kt->te); + continue; + } + + // fprintf(stderr, "q-utg%.6ul (qs: %u, qe: %u), t-utg%.6ul (ts: %u, te: %u)\n", + // kt->qn + 1, kt->qs, kt->qe, kt->tn + 1, kt->ts, kt->te); + o->k_trans.a[m] = *kt; + m++; + } + // else + // { + // fprintf(stderr, "******delete: q-utg%.6ul (qs: %u, qe: %u), t-utg%.6ul (ts: %u, te: %u)\n", + // kt->qn + 1, kt->qs, kt->qe, kt->tn + 1, kt->ts, kt->te); + // } + } + o->k_trans.n = m; + } + } + + if(VERBOSE >= 1) + { + fprintf(stderr, "[M::%s] removed %d long tips\n", + __func__, n_reduced); + } + fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime); + asg_cleanup(g); + asg_symm(g); + free(b.b.a); + kv_destroy(tt.a); + kv_destroy(stack.a); + free(vis); + destory_pdq(&pq_p); + destory_pdq(&pq_a); + free(path_p); + free(path_q); + destroy_clean_mul_t(&cl); + return n_reduced; +} + + +int asg_arc_decompress(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, +R_to_U* ruIndex, utg_trans_t *o) +{ + double startTime = Get_T(); + asg_arc_t *as = NULL, *pm = NULL; + u_trans_t *kt = NULL; + uint32_t v, s, sv, nc, ns, n_vtx = g->n_seq * 2, n_reduced = 0, convex, pi, qn, tn; + uint32_t return_flag, k, m, i, p_n, a_n, *a_a = NULL; + uint8_t fp = 1, fa = 2, found; + long long ll, tmp, max_stop_nodeLen, max_stop_baseLen; + kvec_t_u32_warp tt; kv_init(tt.a); + kvec_t_u32_warp stack; kv_init(stack.a); + uint8_t *vis = NULL; CALLOC(vis, n_vtx); + buf_t b; + memset(&b, 0, sizeof(buf_t)); + pdq pq_p, pq_a; + init_pdq(&pq_p, g->n_seq<<1); + init_pdq(&pq_a, g->n_seq<<1); + uint32_t *path_p = NULL, *path_q = NULL; + CALLOC(path_p, g->n_seq<<1); + CALLOC(path_q, g->n_seq<<1); + + for (v = 0; v < n_vtx; v++) + { + if(g->seq[v>>1].del || g->seq[v>>1].c == ALTER_LABLE) continue; + if(asg_arc_n(g, v) == 0 || get_real_length(g, v, NULL) != 1) continue; + get_real_length(g, v, &s); + if(get_real_length(g, s^1, NULL) < 2) continue; + return_flag = get_unitig(g, ug, v^1, &convex, &ll, &tmp, &max_stop_nodeLen,&max_stop_baseLen, 1, NULL); + if(return_flag == LOOP) continue; + + as = asg_arc_a(g, convex); ns = asg_arc_n(g, convex); + for (i = nc = 0; i < ns; i++) + { + if(as[i].del) continue; + if(get_real_length(g, as[i].v^1, NULL) < 2) break; + nc++; + } + if(nc == 0 || i < ns) continue; + + tt.a.n = 0; s^=1; sv = v^1; + dfs_set(g, sv, s>>1, &stack, &tt, vis, fp); + p_n = tt.a.n; + as = asg_arc_a(g, s); ns = asg_arc_n(g, s); found = 0; + for (i = 0, pm = NULL; i < ns; i++) + { + if(as[i].del || as[i].v != sv) continue; + pm = &(as[i]); + break; + } + + for (i = 0; i < ns; i++) + { + if(as[i].del || as[i].v == sv) continue; + nc = dfs_set(g, as[i].v, s>>1, &stack, &tt, vis, fa); + if(nc && check_trans_relation_by_path(sv, as[i].v, &pq_p, path_p, &(o->b0), &pq_a, path_q, &(o->b1), g, + vis, fp+fa, nc, 0.45, NULL)) + { + found = 1; + } + a_a = tt.a.a + p_n; a_n = tt.a.n - p_n; + for (k = 0; k < a_n; k++) + { + if(vis[a_a[k]]&fa) vis[a_a[k]] -= fa; + } + tt.a.n = p_n; + if(found) break; + } + a_a = tt.a.a; a_n = p_n; + for (k = 0; k < a_n; k++) vis[a_a[k]] = 0; + + if(found) + { + b.b.n = 0; + get_unitig(g, ug, sv, &convex, &ll, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 1, &b); + for (k = 0; k < b.b.n; k++) + { + g->seq[b.b.a[k]>>1].c = ALTER_LABLE; + } + + for (k = 0; k < b.b.n; k++) + { + asg_seq_drop(g, b.b.a[k]>>1); + } + + // uint32_t ui; + // fprintf(stderr, "\n"); + // for (ui = 0; ui < o->b0.b.n; ui++) + // { + // fprintf(stderr, "p-utg%.6ul\n", (o->b0.b.a[ui]>>1)+1); + // } + // for (ui = 0; ui < o->b1.b.n; ui++) + // { + // fprintf(stderr, "a-utg%.6ul\n", (o->b1.b.a[ui]>>1)+1); + // } + + + if(o->b0.b.n < b.b.n) fprintf(stderr, "ERROR-ll\n"); + o->b0.b.n = b.b.n; pi = o->k_trans.n; + collect_trans_ovlp(__func__, &(o->b0), pm->ol, &(o->b1), as[i].ol, ug, o); + + for (k = m = pi; k < o->k_trans.n; k++) + { + kt = &(o->k_trans.a[k]); + if(is_hap_ovlp(ug, read_sg, reverse_sources, ruIndex, &(kt->qn), 1, &(kt->tn), 1) || + is_hap_ovlp(ug, read_sg, reverse_sources, ruIndex, &(kt->tn), 1, &(kt->qn), 1)) + { + qn = get_dir_v(kt->qn, &(o->b0)); + tn = get_dir_v(kt->tn, &(o->b1)); + + if(dfs_reach(o->cug->g, qn, tn, &stack, &tt, vis) || dfs_reach(o->cug->g, tn, qn, &stack, &tt, vis)) + { + // fprintf(stderr, "<<<<<qn + 1, kt->qs, kt->qe, kt->tn + 1, kt->ts, kt->te); + continue; + } + + // fprintf(stderr, "q-utg%.6ul (qs: %u, qe: %u), t-utg%.6ul (ts: %u, te: %u)\n", + // kt->qn + 1, kt->qs, kt->qe, kt->tn + 1, kt->ts, kt->te); + o->k_trans.a[m] = *kt; + m++; + } + // else + // { + // fprintf(stderr, "******delete: q-utg%.6ul (qs: %u, qe: %u), t-utg%.6ul (ts: %u, te: %u)\n", + // kt->qn + 1, kt->qs, kt->qe, kt->tn + 1, kt->ts, kt->te); + // } + } + o->k_trans.n = m; + } + } + + if(VERBOSE >= 1) + { + fprintf(stderr, "[M::%s] removed %d long tips\n", + __func__, n_reduced); + fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime); + } + + asg_cleanup(g); + asg_symm(g); + free(b.b.a); + kv_destroy(tt.a); + kv_destroy(stack.a); + free(vis); + destory_pdq(&pq_p); + destory_pdq(&pq_a); + free(path_p); + free(path_q); + return n_reduced; +} \ No newline at end of file diff --git a/tovlp.h b/tovlp.h new file mode 100644 index 0000000..4b9107b --- /dev/null +++ b/tovlp.h @@ -0,0 +1,15 @@ +#ifndef __TOVLP__ +#define __TOVLP__ +#include +#include "Overlaps.h" + +utg_trans_t *init_utg_trans_t(ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, R_to_U* ruIndex, asg_t *read_g, int max_hang, int min_ovlp); +void destroy_utg_trans_t(utg_trans_t **o); +void asg_bub_collect_ovlp(ma_ug_t *ug, uint32_t v0, buf_t *b, utg_trans_t *o); +void collect_trans_ovlp(const char* cmd, buf_t* pri, uint64_t pri_offset, buf_t* aux, uint64_t aux_offset, +ma_ug_t *ug, utg_trans_t *o); +int asg_arc_decompress(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, +R_to_U* ruIndex, utg_trans_t *o); +int asg_arc_decompress_mul(asg_t *g, ma_ug_t *ug, asg_t *read_sg, uint32_t positive_flag, uint32_t negative_flag, +ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, utg_trans_t *o); +#endif