diff --git a/Overlaps.cpp b/Overlaps.cpp index da3000b..7fbd9b9 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -12712,104 +12712,6 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp, kvec_asg_a /*******************************for debug************************************/ } -void hic_clean(asg_t* read_g) -{ - uint32_t n_vtx, v, u; - uint64_t i, k, k_i, tLen, v_occ, u_occ, utg_occ; - double bub_rate = 0.1; - ma_ug_t *ug = NULL; - ug = ma_ug_gen_primary(read_g, PRIMARY_LABLE); - n_vtx = ug->g->n_seq * 2; - buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); - ///for (i = 0, tLen = 1; i < ug->u.n; i++) tLen += ug->u.a[i].len; - tLen = get_bub_pop_max_dist_advance(ug->g, &b); - uint8_t* bs_flag = (uint8_t*)calloc(n_vtx, 1); - kvec_t(uint32_t) ax; - kv_init(ax); - for (v = 0; v < ug->g->n_seq; ++v) - { - if(ug->g->seq[v].del) continue; - ug->g->seq[v].c = PRIMARY_LABLE; - EvaluateLen(ug->u, v) = ug->u.a[v].n; - } - - - for (v = 0; v < n_vtx; ++v) - { - if(ug->g->seq[v>>1].del) continue; - if(asg_arc_n(ug->g, v) < 2) continue; - if(bs_flag[v] != 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)) - { - //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; - } - } - - - for (v = 0; v < n_vtx; ++v) - { - if(bs_flag[v] !=2) 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)) - { - //note b.b include end, does not include beg - for (i = v_occ = ax.n = 0; i < b.b.n; i++) - { - if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue; - v_occ += ug->u.a[b.b.a[i]>>1].n; - kv_push(uint32_t, ax, b.b.a[i]>>1); - } - - for (i = 0; i < ax.n; i++) - { - for (k = 0; k < 2; k++) - { - u = (ax.a[i]<<1) + k; - if(asg_arc_n(ug->g, u) < 2) continue; - if(asg_bub_pop1_primary_trio(ug->g, NULL, u, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0)) - { - for (k_i = u_occ = utg_occ = 0; k_i < b.b.n; k_i++) - { - if(b.b.a[k_i]==u || b.b.a[k_i]==b.S.a[0]) continue; - u_occ += ug->u.a[b.b.a[k_i]>>1].n; - utg_occ++; - } - - if(u_occ >= v_occ*bub_rate) continue; - if(u_occ > 3) continue; - if(utg_occ > 2) continue; - asg_bub_pop1_primary_trio(ug->g, NULL, u, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0); - } - } - } - } - } - - ma_utg_t* m = NULL; - for (v = 0; v < ug->g->n_seq; ++v) - { - if(ug->g->seq[v].del) continue; - if(ug->g->seq[v].c != ALTER_LABLE) continue; - m = &(ug->u.a[v]); - if(m->m == 0) continue; - for (k = 0; k < m->n; k++) - { - asg_seq_del(read_g, m->a[k]>>33); - } - } - - asg_cleanup(read_g); - free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); free(bs_flag); - ma_ug_destroy(ug); - kv_destroy(ax); -} - void write_trans_chain(trans_chain* t_ch, const char *fn) { char *buf = (char*)calloc(strlen(fn) + 25, 1); @@ -12899,7 +12801,6 @@ 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, long long gap_fuzz, bub_label_t* b_mask_t) { - hic_clean(sg); ug_opt_t opt; memset(&opt, 0, sizeof(opt)); opt.coverage_cut = coverage_cut; opt.sources = sources; @@ -13675,8 +13576,6 @@ 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, bub_label_t* b_mask_t) { - hic_clean(sg); - kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); ma_ug_t *ug = NULL; @@ -29782,7 +29681,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, NULL, 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)) { //note b.b include end, does not include beg for (i = path = 0; i < b.b.n; i++) @@ -29871,7 +29770,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, NULL, 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); n_pop++; } } @@ -29962,7 +29861,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, NULL, 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)) { bs_flag[v] = 2; bs_flag[b.S.a[0]^1] = 2; //note b.b include end, does not include beg @@ -30015,7 +29914,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, NULL, 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); n_pop++; } diff --git a/hic.cpp b/hic.cpp index d9763c1..104755b 100644 --- a/hic.cpp +++ b/hic.cpp @@ -2344,7 +2344,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, NULL, 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)) { //note b.b include end, does not include beg i = b.b.n + 1; @@ -2376,7 +2376,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, NULL, 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)) { kv_push(uint64_t, bub->pathLen, pathLen); //beg is v, end is b.S.a[0]