diff --git a/CommandLines.h b/CommandLines.h index 84de370..1e4034e 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.19.0-r555" +#define HA_VERSION "0.19.0-r558" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index 92cf3bf..976e1a8 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -13528,7 +13528,7 @@ void hic_clean_adv(asg_t *sg, ug_opt_t *uopt) { uint32_t i, k, m, z, v, w, mk; ma_utg_t *u = NULL; uint32_t *ba, bn, n_vtx, beg, end, n0, n1; ma_utg_t *mz = NULL; ma_ug_t *ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); n_vtx = ug->g->n_seq<<1; double bub_rate = 0.1; - uint8_t *bf = NULL; bubble_type *bub = gen_bubble_chain(sg, ug, uopt, &bf); + uint8_t *bf = NULL; bubble_type *bub = gen_bubble_chain(sg, ug, uopt, &bf, ((asm_opt.polyploidy>2)?1:0)); uint64_t tLen, vocc, socc, pocc; buf_t b; memset(&b, 0, sizeof(buf_t)); CALLOC(b.a, n_vtx); REALLOC(bf, n_vtx); memset(bf, 0, sizeof((*bf))*n_vtx); kvec_t(uint64_t) buf; kv_init(buf); n0 = n1 = 0; @@ -15605,7 +15605,7 @@ void clean_u_trans_t_idx(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g) void gen_bp_phasing(ug_opt_t *opt, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *sg) { uint8_t *bf = NULL; - bubble_type *bub = gen_bubble_chain(sg, ug, opt, &bf); free(bf); + bubble_type *bub = gen_bubble_chain(sg, ug, opt, &bf, 0); free(bf); filter_u_trans(ta, asm_opt.is_bub_trans, asm_opt.is_topo_trans, asm_opt.is_read_trans, asm_opt.is_base_trans); // dbg_prt_utg_trans(ta, ug, "pre"); if(asm_opt.is_base_trans) { @@ -33734,7 +33734,7 @@ bub_label_t *b_mask_t, uint32_t is_trio, int32_t ul_aln_round, char *o_file, con int32_t k, strl = strlen(bin_file)+1, kt, cl, sl; char *id = NULL; renew_g(src, rev_src, n_read, readLen, cov, ruIndex, sg, mini_overlap_length, max_hang_length, uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, asm_opt.max_short_tip, b_mask_t, - is_trio, o_file, bin_file, (ul_aln_round<=1)?1:0, 1, /**((is_trio)?(0):(1))**/1); + is_trio, o_file, bin_file, (ul_aln_round<=1)?1:0, 1, /**((is_trio)?(0):(1))**/((asm_opt.polyploidy<=2)?1:0)); gen_ug_opt_t(uopt, *src, *rev_src, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, *readLen, *cov, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, b_mask_t); ug_ext_gfa(uopt, *sg, ug_ext_len); diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 06c884c..e18a5de 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -2995,7 +2995,7 @@ void ug_ext_gfa(ug_opt_t *uopt, asg_t *sg, uint32_t max_len) } -bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het) +bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het, uint8_t avoid_het) { kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); @@ -3013,7 +3013,14 @@ bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t ** // fprintf(stderr, "1[M::%s]\n", __func__); CALLOC(bub, 1); (*ir_het) = cov->t_ch->ir_het; cov->t_ch->ir_het = NULL; cov->is_r_het = NULL; - identify_bubbles(ug, bub, (*ir_het), NULL); + if(!avoid_het) { + identify_bubbles(ug, bub, (*ir_het), NULL); + } else { + uint8_t *pr_het; CALLOC(pr_het, sg->n_seq); + identify_bubbles_recal_poy(sg, ug, bub, pr_het, uopt->sources, uopt->ruIndex, NULL); + free(pr_het); + } + // fprintf(stderr, "2[M::%s]\n", __func__); kv_destroy(new_rtg_edges.a); destory_hap_cov_t(&cov); if(asm_opt.purge_level_primary == 0) {///all nodes are het @@ -13621,7 +13628,12 @@ void renew_usg_t_bub(ul_resolve_t *uidx, usg_t *ng, uint32_t *id_map, uint8_t *f // fprintf(stderr, "[M::%s] homozygous read coverage threshold: %d\n", __func__, asm_opt.hom_global_coverage_set? // asm_opt.hom_global_coverage:(int)(((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE))); // identify_bubbles(ug, uidx->bub, uidx->r_het, NULL); - identify_bubbles_recal(uidx->sg, ug, uidx->bub, uidx->r_het, uidx->uopt->sources, uidx->uopt->ruIndex, NULL); + if(asm_opt.polyploidy <= 2) { + identify_bubbles_recal(uidx->sg, ug, uidx->bub, uidx->r_het, uidx->uopt->sources, uidx->uopt->ruIndex, NULL); + } else { + identify_bubbles_recal_poy(uidx->sg, ug, uidx->bub, uidx->r_het, uidx->uopt->sources, uidx->uopt->ruIndex, NULL); + } + // fprintf(stderr, "0[M::%s::] f[51]::%u\n", __func__, ff[51]); for (i = 0; i < ng->n; i++) { if(id_map[i] != (uint32_t)-1) { @@ -17234,18 +17246,18 @@ ul_renew_t *ropt, const char *bin_file, uint64_t free_uld, uint64_t is_bridg, ui // fprintf(stderr, "2[M::%s]\n", __func__); // exit(1); - // char* gfa_name = NULL; MALLOC(gfa_name, strlen(o_file)+strlen(bin_file)+50); - // sprintf(gfa_name, "%s.%s", o_file, bin_file); - // print_debug_gfa(sg, init_ug, uopt->coverage_cut, gfa_name, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0); - // print_debug_gfa(sg, init_ug, uopt->coverage_cut, gfa_name, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1); - // free(gfa_name); + char* gfa_name = NULL; MALLOC(gfa_name, strlen(o_file)+strlen(bin_file)+50); + sprintf(gfa_name, "%s.%s", o_file, bin_file); + print_debug_gfa(sg, init_ug, uopt->coverage_cut, gfa_name, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0); + print_debug_gfa(sg, init_ug, uopt->coverage_cut, gfa_name, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1); + free(gfa_name); filter_sg_by_ug(sg, init_ug, uopt); // fprintf(stderr, "-0-[M::%s]\tUL_INF.a[25].rlen::%u\n", __func__, UL_INF.a[25].rlen); // print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-0"); - bub = gen_bubble_chain(sg, init_ug, uopt, &r_het); + bub = gen_bubble_chain(sg, init_ug, uopt, &r_het, ((asm_opt.polyploidy>2)?1:0)); // fprintf(stderr, "4[M::%s]\n", __func__); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-1"); ul_resolve_t *uidx = init_ul_resolve_t(sg, init_ug, bub, &UL_INF, uopt, r_het); diff --git a/gfa_ut.h b/gfa_ut.h index 6b47493..32c5675 100644 --- a/gfa_ut.h +++ b/gfa_ut.h @@ -35,7 +35,7 @@ void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd); asg_t *gen_ng(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **cov, R_to_U *ruI, uint64_t scaffold_len); void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover); // void print_raw_u2rgfa_seq(all_ul_t *aln, R_to_U* rI, uint32_t is_detail); -bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het); +bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het, uint8_t avoid_het); void filter_sg_by_ug(asg_t *rg, ma_ug_t *ug, ug_opt_t *uopt); void ug_ext_gfa(ug_opt_t *uopt, asg_t *sg, uint32_t max_len); void update_sg_uo(asg_t *g, ma_hit_t_alloc *src); diff --git a/hic.cpp b/hic.cpp index d60847c..59ec0c1 100644 --- a/hic.cpp +++ b/hic.cpp @@ -2520,7 +2520,8 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t // fprintf(stderr, "-bub->index[18759]: %u, bub->num.n: %u\n", (uint32_t)bub->index[18759], bub->num.n); } -uint32_t get_unitig_het_fly(ma_ug_t* ug, uint32_t uid, asg_t* sg, int64_t het_cov_thres, +uint32_t get_unitig_het_fly(ma_ug_t* ug, uint32_t uid, asg_t* sg, /**int64_t het_cov_thres,**/ +int64_t m_het_cov_thres, int64_t m_hom_cov_thres, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag, uint32_t m_het_occ, uint32_t m_het_label, uint32_t p_het_label, uint32_t n_het_label) { @@ -2596,8 +2597,10 @@ uint32_t p_het_label, uint32_t n_het_label) // fprintf(stderr, "[M::%s::uid->%u] u->n::%u, C_bases::%ld, R_bases::%ld, het_cov_thres::%ld, m_het_occ::%u\n", // __func__, uid, (uint32_t)u->n, C_bases, R_bases, het_cov_thres, m_het_occ); // } - if((cov <= (het_cov_thres*1.333333)) && (u->n >= m_het_occ)) return m_het_label; ///must het - if((cov >= (het_cov_thres*1.6))) return n_het_label; ///hom + // if((cov <= (het_cov_thres*1.333333)) && (u->n >= m_het_occ)) return m_het_label; ///must het + // if((cov >= (het_cov_thres*1.6))) return n_het_label; ///hom + if((cov <= m_het_cov_thres) && (u->n >= m_het_occ)) return m_het_label; ///must het + if(cov >= m_hom_cov_thres) return n_het_label; ///hom return p_het_label; ///potential het } @@ -2608,9 +2611,20 @@ kv_u_trans_t *ref) if (!ug->g->is_symm) asg_symm(ug->g); uint32_t v, n_vtx = ug->g->n_seq * 2, i, k, mode = (((uint32_t)-1)<<2); uint32_t beg, sink, n, *a, n_occ; - uint64_t pathLen; + uint64_t pathLen, hom_cov, het_cov, m_het_cov, m_hom_cov; bub->ug = ug; bub->b_bub = bub->b_end_bub = bub->tangle_bub = bub->cross_bub = bub->mess_bub = 0; + + if(asm_opt.hom_global_coverage_set) { + hom_cov = asm_opt.hom_global_coverage; + } else { + hom_cov = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE); + } + het_cov = hom_cov/asm_opt.polyploidy; + m_het_cov = hom_cov - het_cov + (het_cov*0.333333); + m_hom_cov = hom_cov - het_cov + (het_cov*0.6); + // fprintf(stderr, "hom_cov::%lu, het_cov::%lu, m_het_cov::%lu, m_hom_cov::%lu\n", hom_cov, het_cov, m_het_cov, m_hom_cov); + if(bub->round_id == 0) { buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); @@ -2697,20 +2711,11 @@ kv_u_trans_t *ref) kv_push(uint32_t, bub->num, bub->list.n); free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); bub->f_bub = bub->num.n - 1; ///bub->s_bub = bub->num.n - 1; - - uint64_t dip_thre_max; memset(r_het_flag, 0, sizeof((*r_het_flag))*sg->n_seq); - if(asm_opt.hom_global_coverage_set) { - dip_thre_max = asm_opt.hom_global_coverage; - } else { - dip_thre_max = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE); - } - // dip_thre_max = (double)(dip_thre_max) - (((double)(dip_thre_max)*0.5)/asm_opt.polyploidy); - dip_thre_max = (double)(dip_thre_max) - ((double)(dip_thre_max)/asm_opt.polyploidy); for (i = 0; i < ug->g->n_seq; i++) { - bub->index[i] = get_unitig_het_fly(ug, i, sg, dip_thre_max, sources, ruIndex, - r_het_flag, 20, M_het(*bub), P_het(*bub), (uint32_t)-1); + bub->index[i] = get_unitig_het_fly(ug, i, sg, m_het_cov, m_hom_cov, + sources, ruIndex, r_het_flag, 20, M_het(*bub), P_het(*bub), (uint32_t)-1); } for (i = 0; i < bub->f_bub; i++) @@ -2782,6 +2787,301 @@ kv_u_trans_t *ref) // fprintf(stderr, "-bub->index[18759]: %u, bub->num.n: %u\n", (uint32_t)bub->index[18759], bub->num.n); } +void reset_inner_bub_het_poy(asg_t* sg, ma_ug_t* ug, bubble_type* bub, uint32_t bid, uint64_t tLen, buf_t *b, +uint64_t m_het_cov, uint64_t m_hom_cov, uint8_t *r_het_flag, ma_hit_t_alloc* sources, R_to_U* ruIndex) +{ + uint32_t beg, sink, *ba, bn, m, v, z, socc; + get_bubbles(bub, bid, &beg, &sink, &ba, &bn, NULL); + for (m = 0; m < bn; m++) { + v = ba[m]; + if(IF_HOM((v>>1), *bub)) continue; + if(ug->g->seq[v>>1].del) continue; + if(bub->index[v>>1] == (uint32_t)-1) continue; + + if(get_unitig_het_fly(ug, v>>1, sg, m_het_cov, m_hom_cov, sources, ruIndex, r_het_flag, 20, 1, 0, (uint32_t)-1) == 1) { + continue; + } + + if(asg_arc_n(ug->g, v) < 2) continue; + if(get_real_length(ug->g, v, NULL) < 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, NULL)) { + //beg is v, end is b.S.a[0] + //note b.b include end, does not include beg + for (z = socc = 0; z < b->b.n; z++) { + if(b->b.a[z]==v || b->b.a[z]==b->S.a[0]) continue; + socc += ug->u.a[b->b.a[z]>>1].n; + } + if((socc > 10) && (socc > (ug->u.a[v>>1].n*3))) { + bub->index[v>>1] = (uint32_t)-1; + } + } + + + v ^= 1; + if(bub->index[v>>1] == (uint32_t)-1) continue; + if(asg_arc_n(ug->g, v) < 2) continue; + if(get_real_length(ug->g, v, NULL) < 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, NULL)) { + //beg is v, end is b.S.a[0] + //note b.b include end, does not include beg + for (z = socc = 0; z < b->b.n; z++) { + if(b->b.a[z]==v || b->b.a[z]==b->S.a[0]) continue; + socc += ug->u.a[b->b.a[z]>>1].n; + } + if((socc > 10) && (socc > (ug->u.a[v>>1].n*3))) { + bub->index[v>>1] = (uint32_t)-1; + } + } + } +} + +uint32_t deter_unitig_het_fly(ma_ug_t* ug, uint32_t uid, asg_t* sg, int64_t het_cov_thres, +ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag) +{ + ma_utg_t *u = &(ug->u.a[uid]); + uint32_t k, i, j, rId, nv, tn, is_Unitig; + asg_arc_t *av = NULL; ma_hit_t *h; + int64_t R_bases = 0, C_bases = 0, cov; + + ///set + u = &(ug->u.a[uid]); + for (k = 0; k < u->n; k++) { + rId = u->a[k]>>33; + r_flag[rId] = 1; + } + for (i = 0; i < 2; i++) { + nv = asg_arc_n(ug->g, (uid<<1)+i); + av = asg_arc_a(ug->g, (uid<<1)+i); + for (j = 0; j < nv; j++) { + u = &(ug->u.a[av[j].v>>1]); + for (k = 0; k < u->n; k++) { + rId = u->a[k]>>33; + r_flag[rId] = 2; + } + } + } + + u = &(ug->u.a[uid]); + for (k = 0; k < u->n; k++) { + if(u->a[k] == (uint64_t)-1) continue; + rId = u->a[k]>>33; + R_bases += sg->seq[rId].len; + for (j = 0; j < (uint64_t)(sources[rId].length); j++) { + h = &(sources[rId].buffer[j]); + if(h->el != 1) continue; + tn = Get_tn((*h)); + if(sg->seq[tn].del == 1) { + ///get the id of read that contains it + get_R_to_U(ruIndex, tn, &tn, &is_Unitig); + if(tn == (uint32_t)-1 || is_Unitig == 1 || sg->seq[tn].del == 1) continue; + } + if(sg->seq[tn].del == 1) continue; + if(r_flag[tn] == 0) continue; + if(r_flag[tn] == 1) { + C_bases += (Get_qe((*h)) - Get_qs((*h))); + } + if(r_flag[tn] == 2) { + C_bases += ((Get_qe((*h)) - Get_qs((*h)))/2); + } + } + } + + ///reset + u = &(ug->u.a[uid]); + for (k = 0; k < u->n; k++) { + rId = u->a[k]>>33; + r_flag[rId] = 0; + } + for (i = 0; i < 2; i++) { + nv = asg_arc_n(ug->g, (uid<<1)+i); + av = asg_arc_a(ug->g, (uid<<1)+i); + for (j = 0; j < nv; j++) { + u = &(ug->u.a[av[j].v>>1]); + for (k = 0; k < u->n; k++) { + rId = u->a[k]>>33; + r_flag[rId] = 0; + } + } + } + + u = &(ug->u.a[uid]); cov = 0; + if(R_bases > 0) cov = C_bases/R_bases; + if(cov <= het_cov_thres) return 1; + return 0; +} + +void identify_bubbles_recal_poy(asg_t* sg, ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, ma_hit_t_alloc* sources, R_to_U* ruIndex, +kv_u_trans_t *ref) +{ + asg_cleanup(ug->g); + if (!ug->g->is_symm) asg_symm(ug->g); + uint32_t v, n_vtx = ug->g->n_seq * 2, i, k, mode = (((uint32_t)-1)<<2); + uint32_t beg, sink, n, *a, n_occ; + uint64_t pathLen, hom_cov, het_cov, m_het_cov, m_hom_cov; + bub->ug = ug; + bub->b_bub = bub->b_end_bub = bub->tangle_bub = bub->cross_bub = bub->mess_bub = 0; + if(asm_opt.hom_global_coverage_set) { + hom_cov = asm_opt.hom_global_coverage; + } else { + hom_cov = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE); + } + het_cov = hom_cov/asm_opt.polyploidy; + m_het_cov = het_cov + (het_cov*0.333333);///hom_cov - het_cov + (het_cov*0.333333); + m_hom_cov = het_cov + (het_cov*0.6);///hom_cov - het_cov + (het_cov*0.6); + // fprintf(stderr, "hom_cov::%lu, het_cov::%lu, m_het_cov::%lu, m_hom_cov::%lu\n", hom_cov, het_cov, m_het_cov, m_hom_cov); + + buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); + uint64_t tLen = get_bub_pop_max_dist_advance(ug->g, &b); + kv_init(bub->list); kv_init(bub->num); kv_init(bub->pathLen); + kv_init(bub->b_s_idx); kv_malloc(bub->b_s_idx, ug->g->n_seq); + bub->b_ug = NULL; kv_init(bub->chain_weight); + bub->b_s_idx.n = ug->g->n_seq; + memset(bub->b_s_idx.a, -1, bub->b_s_idx.n * sizeof(uint64_t)); + CALLOC(bub->index, n_vtx); + for (i = 0; i < ug->g->n_seq; i++) ug->g->seq[i].c = 0; + 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((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, 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; + bub->index[b.b.a[i]] &= mode; bub->index[b.b.a[i]] += 1; + bub->index[b.b.a[i]^1] &= mode; bub->index[b.b.a[i]^1] += 1; + } + bub->index[v] &= mode; bub->index[v] += 2; + bub->index[b.S.a[0]^1] &= mode; bub->index[b.S.a[0]^1] += 3; + } + } + + kvec_t_u32_warp stack, result; + kv_init(stack.a); kv_init(result.a); + 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, NULL)) + { + //note b.b include end, does not include beg + i = b.b.n + 1; + if(b.b.n == 2 || b.b.n == 3 || b.b.n == 5) + { + for (i = 0; i < b.b.n; i++) + { + if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue; + dfs_bubble(ug->g, &stack, &result, b.b.a[i]>>1, v>>1, b.S.a[0]>>1); + if((result.a.n + 3) != b.b.n && (result.a.n + 2) != b.b.n) break; + } + } + + if(i == b.b.n) + { + kv_push(uint32_t, bub->num, v); + } + else + { + kv_push(uint32_t, bub->num, v + (1<<31)); + } + } + } + kv_destroy(stack.a); kv_destroy(result.a); + radix_sort_u32(bub->num.a, bub->num.a + bub->num.n); + bub->s_bub = 0; + for (k = 0; k < bub->num.n; k++) + { + 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, NULL)) + { + kv_push(uint64_t, bub->pathLen, pathLen); + //beg is v, end is b.S.a[0] + kv_push(uint32_t, bub->list, v); + kv_push(uint32_t, bub->list, b.S.a[0]^1); + + //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; + kv_push(uint32_t, bub->list, b.b.a[i]); + } + } + } + kv_push(uint32_t, bub->num, bub->list.n); + // free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); + bub->f_bub = bub->num.n - 1; ///bub->s_bub = bub->num.n - 1; + + + memset(r_het_flag, 0, sizeof((*r_het_flag))*sg->n_seq); + for (i = 0; i < ug->g->n_seq; i++) { + bub->index[i] = get_unitig_het_fly(ug, i, sg, m_het_cov, m_hom_cov, + sources, ruIndex, r_het_flag, 20, M_het(*bub), P_het(*bub), (uint32_t)-1); + } + + for (i = 0; i < bub->f_bub; i++) { + get_bubbles(bub, i, &beg, &sink, &a, &n, &pathLen); + for (v = n_occ = 0; v < n; v++) { + bub->index[(a[v]>>1)] = i; + n_occ += ug->u.a[a[v]>>1].n; + } + + if(n_occ > 3) { + if(bub->index[(beg>>1)] != M_het(*bub)) bub->index[(beg>>1)] = (uint32_t)-1; + if(bub->index[(sink>>1)] != M_het(*bub)) bub->index[(sink>>1)] = (uint32_t)-1; + } + + + v = beg>>1; + if(bub->b_s_idx.a[v] == (uint64_t)-1) { + bub->b_s_idx.a[v] <<= 32; + bub->b_s_idx.a[v] |= i; + } else if((bub->b_s_idx.a[v] & 0xffffffff00000000) == 0xffffffff00000000) { + bub->b_s_idx.a[v] <<= 32; + bub->b_s_idx.a[v] |= i; + } + + + v = sink>>1; + if(bub->b_s_idx.a[v] == (uint64_t)-1) { + bub->b_s_idx.a[v] <<= 32; + bub->b_s_idx.a[v] |= i; + } else if((bub->b_s_idx.a[v] & 0xffffffff00000000) == 0xffffffff00000000) { + bub->b_s_idx.a[v] <<= 32; + bub->b_s_idx.a[v] |= i; + } + } + for (i = 0; i < ug->g->n_seq; i++) { + if(bub->index[i] == M_het(*bub)) bub->index[i] = P_het(*bub); + } + + bub->b_g = NULL; + bub->b_ug = NULL; + build_bub_graph(ug, bub); + + + ///make het nodes to be hom + ma_utg_t *u = NULL; + for (i = 0; i < bub->b_ug->u.n; i++) { + u = &(bub->b_ug->u.a[i]); + if(u->n == 0) continue; + for (k = 0; k < u->n; k++) { + reset_inner_bub_het_poy(sg, ug, bub, u->a[k]>>33, tLen, &b, m_het_cov, m_hom_cov, r_het_flag, sources, ruIndex); + } + } + free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); + + for (i = 0; i < ug->g->n_seq; i++) { + if(IF_HOM(i, *bub)) continue; + if(deter_unitig_het_fly(ug, i, sg, het_cov*1.15, sources, ruIndex, r_het_flag)) continue; + bub->index[i] = (uint32_t)-1; + } + // fprintf(stderr, "-bub->index[18759]: %u, bub->num.n: %u\n", (uint32_t)bub->index[18759], bub->num.n); +} + void print_bubbles(ma_ug_t* ug, bubble_type* bub, kvec_pe_hit* hits, hc_links* link, ha_ug_index* idx) { uint64_t tLen, t_utg, i, k; diff --git a/hic.h b/hic.h index 1c3226f..321c044 100644 --- a/hic.h +++ b/hic.h @@ -95,6 +95,8 @@ void destory_bubbles(bubble_type* bub); void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_trans_t *ref); void identify_bubbles_recal(asg_t* sg, ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, ma_hit_t_alloc* sources, R_to_U* ruIndex, kv_u_trans_t *ref); +void identify_bubbles_recal_poy(asg_t* sg, ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, ma_hit_t_alloc* sources, R_to_U* ruIndex, +kv_u_trans_t *ref); void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub); uint32_t connect_bub_occ(bubble_type* bub, uint32_t root_id, uint32_t check_het); void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, uint32_t check_het); diff --git a/inter.cpp b/inter.cpp index 4b54085..eb49c7d 100644 --- a/inter.cpp +++ b/inter.cpp @@ -19599,7 +19599,7 @@ int32_t mini_cut, int32_t chain_cut, int32_t keep_unsymm_arc, ma_ug_t *ug, asg_t if(bub) { opt->bub = bub; } else { - opt->bub = gen_bubble_chain(sg, ug, uopt, &bf); free(bf); + opt->bub = gen_bubble_chain(sg, ug, uopt, &bf, 0); free(bf); } }