From aa164949f27438c26f260a402d5a729b8e64cc9c Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Mon, 31 May 2021 04:19:05 -0400 Subject: [PATCH] backup phasing --- Assembly.cpp | 2 +- Overlaps.cpp | 31 ++++++++++----------- Overlaps.h | 2 ++ hic.cpp | 14 ++++------ rcut.cpp | 77 ++++++++++++++++++++++++++++++++-------------------- 5 files changed, 71 insertions(+), 55 deletions(-) diff --git a/Assembly.cpp b/Assembly.cpp index c508c4f..358bbe2 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1660,7 +1660,7 @@ void ug_idx_build(ma_ug_t *ug, int hap_n) int ha_assemble(void) { // debug_mc_g_t(MC_NAME); - debug_mc_gg_t(MC_NAME, 0, 0); + // debug_mc_gg_t(MC_NAME, 0, 0); extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o); int r, hom_cov = -1, ovlp_loaded = 0; if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) { diff --git a/Overlaps.cpp b/Overlaps.cpp index 9972050..64b7988 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -12837,7 +12837,7 @@ void set_r_het_flag(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, ma_hit_t_all // dip_thre_max *= 0.75; dip_thre_max = (double)(dip_thre_max) - (((double)(dip_thre_max)*0.5)/asm_opt.polyploidy); - fprintf(stderr, "dip_thre_max: %lu\n", dip_thre_max); + // fprintf(stderr, "dip_thre_max: %lu\n", dip_thre_max); for (m = 0; m < ug->g->n_seq; m++) { @@ -13069,7 +13069,7 @@ void write_trans_chain(trans_chain* t_ch, const char *fn) } -trans_chain* load_hc_hits(const char *fn) +trans_chain* load_hc_trans(const char *fn) { uint64_t flag = 0; char *buf = (char*)calloc(strlen(fn) + 25, 1); @@ -13148,7 +13148,7 @@ long long gap_fuzz, bub_label_t* b_mask_t) hap_cov_t *cov = NULL; trans_chain* t_ch = NULL; - if((asm_opt.flag & HA_F_VERBOSE_GFA)) t_ch = load_hc_hits(output_file_name); + if((asm_opt.flag & HA_F_VERBOSE_GFA)) t_ch = load_hc_trans(output_file_name); if(!t_ch) { @@ -13308,7 +13308,7 @@ long long gap_fuzz, bub_label_t* b_mask_t) fclose(output_file); free(gfa_name); - if((asm_opt.flag & HA_F_VERBOSE_GFA)) t_ch = load_hc_hits(output_file_name); + if((asm_opt.flag & HA_F_VERBOSE_GFA)) t_ch = load_hc_trans(output_file_name); if(!t_ch) { t_ch = get_hic_polyploid_trans_chain(ug, sg, coverage_cut, sources, reverse_sources, max_hang, min_ovlp, ruIndex, b_mask_t); @@ -13317,7 +13317,7 @@ long long gap_fuzz, bub_label_t* b_mask_t) if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_trans_chain(t_ch, output_file_name); } - hic_analysis(ug, sg, t_ch, &opt, 1); + hic_analysis(ug, sg, t_ch, &opt, 0); destory_trans_chain(&t_ch); ma_ug_destroy(ug); asg_cleanup(sg); @@ -14101,7 +14101,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) { - fprintf(stderr, "******0******\n"); kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); ma_ug_t *ug = NULL; @@ -14112,13 +14111,13 @@ bub_label_t* b_mask_t) asg_t *copy_sg = copy_read_graph(sg); ma_ug_t *copy_ug = copy_untig_graph(ug); /*******************************for debug************************************/ - // 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); - adjust_utg_advance(copy_sg, copy_ug, reverse_sources, ruIndex, b_mask_t); - get_utg_ovlp(©_ug, copy_sg, sources, reverse_sources, coverage_cut, - ruIndex, max_hang, min_ovlp, &new_rtg_edges, b_mask_t, NULL); - exit(1); + 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); + // adjust_utg_advance(copy_sg, copy_ug, reverse_sources, ruIndex, b_mask_t); + // get_utg_ovlp(©_ug, copy_sg, sources, reverse_sources, coverage_cut, + // ruIndex, max_hang, min_ovlp, &new_rtg_edges, b_mask_t, NULL); + // exit(1); /*******************************for debug************************************/ print_utg(copy_ug, copy_sg, coverage_cut, output_file_name, sources, ruIndex, max_hang, min_ovlp, &new_rtg_edges); @@ -30816,10 +30815,10 @@ ma_sub_t **coverage_cut_ptr, int debug_g) else if(ha_opt_hic(&asm_opt)) { if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION; - // output_hic_graph(sg, coverage_cut, o_file, sources, reverse_sources, (asm_opt.max_short_tip*2), - // 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, gap_fuzz, &b_mask_t); - output_hic_graph_polyploid(sg, coverage_cut, o_file, sources, reverse_sources, (asm_opt.max_short_tip*2), + output_hic_graph(sg, coverage_cut, o_file, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, gap_fuzz, &b_mask_t); + // output_hic_graph_polyploid(sg, coverage_cut, o_file, sources, reverse_sources, (asm_opt.max_short_tip*2), + // 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, gap_fuzz, &b_mask_t); } else if((asm_opt.flag & HA_F_PARTITION) && (asm_opt.purge_level_primary > 0)) { diff --git a/Overlaps.h b/Overlaps.h index 2151679..bba4914 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -992,6 +992,8 @@ uint32_t set_utg_offset(uint32_t *a, uint32_t a_n, ma_ug_t *ug, asg_t *read_sg, uint32_t only_len); uint64_t get_utg_cov(ma_ug_t *ug, uint32_t uID, asg_t* read_g, const ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag); +trans_chain* load_hc_trans(const char *fn); +char *get_outfile_name(char* output_file_name); #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/hic.cpp b/hic.cpp index d42c54f..4284a23 100644 --- a/hic.cpp +++ b/hic.cpp @@ -785,9 +785,7 @@ ha_ug_index* build_unitig_index(ma_ug_t *ug, int k, uint64_t up_occ, uint64_t lo ha_ug_index* idx = NULL; CALLOC(idx, 1); pldat_t pl; pl.h = idx; pl.is_cnt = 1; double index_time = yak_realtime(), beg_time; - fprintf(stderr, "sa-0\n"); init_ha_ug_index_opt(idx, ug, k, &pl, up_occ, low_occ, thread_num); - fprintf(stderr, "sa-1\n"); beg_time = yak_realtime(); pl.is_cnt = 1; @@ -15733,7 +15731,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o skip_flipping: verbose_het_stat(&bub); - print_kv_weight(&k_trans); + // print_kv_weight(&k_trans); // horder_t *ho = init_horder_t(&sl.hits, idx->uID_bits, idx->pos_mode, idx->read_g, idx->ug, &bub, &(idx->t_ch->k_trans), opt, 3); @@ -15865,11 +15863,11 @@ int hic_short_align_poy(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, resolve_tangles_hic(idx, &bub, &sl.hits, &k_trans); measure_distance(idx, idx->ug, &sl.hits, &link, &bub, &(idx->t_ch->k_trans)); - // if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_psg_t(&s, asm_opt.output_file_name)) - // { - // bub.round_id = bub.n_round; - // goto skip_flipping; - // } + if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_psg_t(&s, asm_opt.output_file_name)) + { + bub.round_id = bub.n_round; + goto skip_flipping; + } s = init_psg_t(11, idx->ug, idx->read_g, opt); for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++) { diff --git a/rcut.cpp b/rcut.cpp index 9e57e99..dbe80a5 100644 --- a/rcut.cpp +++ b/rcut.cpp @@ -611,13 +611,13 @@ mb_nodes_t *update_mb_nodes_t(kv_u_trans_t *ref, mc_match_t* ma, uint32_t occ) return x; } -mb_g_t *init_mb_g_t(mc_g_t *mg, kv_u_trans_t *ref, uint32_t is_sys) +mb_g_t *init_mb_g_t(mc_match_t* e, kv_u_trans_t *ref, uint32_t is_sys) { mc_edge_t *o = NULL; uint32_t i, k, m, n, a_n[2], *a[2], qn, tn, qb, tb; mb_edge_t *t = NULL; mb_g_t *p = NULL; CALLOC(p, 1); - p->u = update_mb_nodes_t(ref, mg->e, mg->ug->g->n_seq); + p->u = update_mb_nodes_t(ref, e, e->n_seq); p->e = NULL; CALLOC(p->e, 1); p->e->n_seq = p->u->u.n; kv_init(p->e->ma); kv_init(p->e->idx); @@ -629,8 +629,8 @@ mb_g_t *init_mb_g_t(mc_g_t *mg, kv_u_trans_t *ref, uint32_t is_sys) for (k = 0; k < a_n[0]; k++) { qn = a[0][k]; - o = pt_a(*(mg->e), qn); - n = pt_n(*(mg->e), qn); + o = pt_a(*e, qn); + n = pt_n(*e, qn); for (m = 0; m < n; m++) { tn = ma_y(o[m]); @@ -650,8 +650,8 @@ mb_g_t *init_mb_g_t(mc_g_t *mg, kv_u_trans_t *ref, uint32_t is_sys) for (k = 0; k < a_n[1]; k++) { qn = a[1][k]; - o = pt_a(*(mg->e), qn); - n = pt_n(*(mg->e), qn); + o = pt_a(*e, qn); + n = pt_n(*e, qn); for (m = 0; m < n; m++) { tn = ma_y(o[m]); @@ -2152,10 +2152,7 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3 } } - fprintf(stderr, "\nBeg-[M::%s::score->%f]\n", __func__, mc_score(mg->e, b)); - print_mc_node(mg->e, b, 3838); - print_mc_node(mg->e, b, 36880); - + // print_mc_node(mg->e, b, 36880); // print_sc(opt, mg, b, sc_opt, n_iter); // mc_reset_z_debug(mg->e, b); // print_sc(opt, mg->e, b, sc_opt, n_iter); @@ -2189,10 +2186,6 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3 } sc_opt = sc; } - - fprintf(stderr, "\n"); - print_mc_node(mg->e, b, 3838); - print_mc_node(mg->e, b, 36880); } for (j = 0; j < b->cc_size; ++j) @@ -2200,8 +2193,7 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3 b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]]; b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]]; } - - fprintf(stderr, "End-[M::%s::score->%f]\n", __func__, mc_score(mg->e, b)); + return n_iter; } @@ -2594,7 +2586,7 @@ void mb_solve_core(mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t is_sys if(!ref) return; double index_time = yak_realtime(); uint32_t st, i; - mb_g_t *mbg = init_mb_g_t(mg, ref, is_sys); + mb_g_t *mbg = init_mb_g_t(mg->e, ref, is_sys); mb_svaux_t *bb; /**************************init**************************/ if(VERBOSE_CUT) @@ -2859,7 +2851,7 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u mc_g_t *mg = init_mc_g_t(ug, read_g, s, renew_s); update_mc_edges(mg, ovlp, ta, t_ch, f_rate, is_sys); - fprintf(stderr, "[M::%s:: # edges: %u]\n", __func__, (uint32_t)mg->e->ma.n); + // fprintf(stderr, "[M::%s:: # edges: %u]\n", __func__, (uint32_t)mg->e->ma.n); mb_solve_core(&opt, mg, ref, is_sys); ///debug_mc_g_t(mg); @@ -2873,7 +2865,7 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u if(ovlp) clean_ovlp_by_mc(mg, ovlp); - print_hap_s(s, ug->u.n); + // print_hap_s(s, ug->u.n); destory_mc_g_t(&mg); } @@ -3237,19 +3229,19 @@ t_w_t mcgg_init_spin(const mc_match_t *ma, mcgg_svaux_t *b) uint32_t n1 = ma_x(*e), n2 = ma_y(*e); if (b->s[n1].s == 0 && b->s[n2].s == 0) { b->s[n1].s = kr_drand_node(n1, b, &(b->s[n1].hc)); - debug_hapM(&(b->s[n1]), "s0"); + // debug_hapM(&(b->s[n1]), "s0"); b->s[n2].s = kr_drand_node_ref(n2, b, &(b->s[n2].hc), b->s[n1].s, e->w>0?1:0); - debug_hapM(&(b->s[n2]), "s1"); + // debug_hapM(&(b->s[n2]), "s1"); } else if(b->s[n1].s == 0) { b->s[n1].s = kr_drand_node_ref(n1, b, &(b->s[n1].hc), b->s[n2].s, e->w>0?1:0); - debug_hapM(&(b->s[n1]), "s2"); + // debug_hapM(&(b->s[n1]), "s2"); } else if(b->s[n2].s == 0) { b->s[n2].s = kr_drand_node_ref(n2, b, &(b->s[n2].hc), b->s[n1].s, e->w>0?1:0); - debug_hapM(&(b->s[n2]), "s3"); + // debug_hapM(&(b->s[n2]), "s3"); } } @@ -3354,19 +3346,40 @@ void inline back_status(mcgg_svaux_t *b, uint32_t id, uint32_t to_opt) } } - +/** static inline mcg_node_t kr_drand_node_new(uint64_t id, mcgg_svaux_t *b) { + if(id == 10070) fprintf(stderr, "id-%lu, h[0]-%u, h[1]-%u\n", id, b->s[id].h[0], b->s[id].h[1]); uint32_t hn = b->s[id].h[0]-1, is_old = 1, k; if((b->s[id].h[1] != (uint16_t)-1) && (kr_drand_r(&b->x) > b->s[id].hw[0])) { hn = b->s[id].h[1]-1; is_old = 0; } + if(id == 10070) fprintf(stderr, "id-%lu, hn-%u, b->m[hn].n-%u, is_old-%u\n", id, hn, b->m[hn].n, is_old); b->x = kr_splitmix64(b->x); k = b->x%(b->m[hn].n-is_old); + if(id == 10070) fprintf(stderr, "id-%lu, hn-%u, k-%u\n", id, hn, k); if(b->m[hn].a[k] == b->s[id].s) k = b->m[hn].n-1; return b->m[hn].a[k]; } +**/ + +static inline mcg_node_t kr_drand_node_new(uint64_t id, mcgg_svaux_t *b) +{ + uint32_t hn = b->s[id].h[0]-1, k; + if((b->s[id].h[1] != (uint16_t)-1) && (kr_drand_r(&b->x) > b->s[id].hw[0])) + { + hn = b->s[id].h[1]-1; + } + b->x = kr_splitmix64(b->x); + k = b->x%(b->m[hn].n); + if(b->m[hn].a[k] == b->s[id].s && b->m[hn].n > 1) + { + k = b->x%(b->m[hn].n-1); + if(b->m[hn].a[k] == b->s[id].s) k = b->m[hn].n-1; + } + return b->m[hn].a[k]; +} static void mcgg_perturb(const mc_opt_t *opt, const mc_match_t *ma, mcgg_svaux_t *b) { @@ -3448,8 +3461,8 @@ uint32_t mcgg_solve_cc(const mc_opt_t *opt, const mc_gg_t *mg, mcgg_svaux_t *b, for (j = 0; j < b->cc_size; ++j) back_status(b, b->cc_node[j], 0); } fprintf(stderr, "\nBeg-[M::%s::score->%f]\n", __func__, mcgg_score(mg->e, b)); - print_mcgg_node(mg->e, b, 3838); - print_mcgg_node(mg->e, b, 36880); + // print_mcgg_node(mg->e, b, 3838); + // print_mcgg_node(mg->e, b, 36880); // tt0 = tt1 = 0; for (k = 0; k < (uint32_t)opt->n_perturb; ++k) { @@ -3458,7 +3471,6 @@ uint32_t mcgg_solve_cc(const mc_opt_t *opt, const mc_gg_t *mg, mcgg_svaux_t *b, else mcgg_perturb_node(opt, mg->e, b, 3); // tt0 += yak_realtime()-t0; - // fprintf(stderr, "\n++(%u) sc_pre: %f\n", k, mcgg_score(mg->e, b)); // t1 = yak_realtime(); sc = mcgg_optimize_local(opt, mg->e, b, &n_iter); // tt1 += yak_realtime()-t1; @@ -3479,9 +3491,9 @@ uint32_t mcgg_solve_cc(const mc_opt_t *opt, const mc_gg_t *mg, mcgg_svaux_t *b, for (j = 0; j < b->cc_size; ++j) back_status(b, b->cc_node[j], 1); sc_opt = sc; } - fprintf(stderr, "\n"); - print_mcgg_node(mg->e, b, 3838); - print_mcgg_node(mg->e, b, 36880); + // fprintf(stderr, "\n"); + // print_mcgg_node(mg->e, b, 3838); + // print_mcgg_node(mg->e, b, 36880); // if((k&31)==0)fprintf(stderr, "+++(%u) sc: %f, sc_opt: %f, tt0: %.3f, tt1: %.3f\n", k, sc, sc_opt, tt0, tt1); } @@ -3680,6 +3692,9 @@ void debug_mc_gg_t(const char* fn, uint32_t update_ta, uint32_t convert_mc_g_t) mc_opt_init(&opt, asm_opt.n_perturb, asm_opt.f_perturb, asm_opt.seed); mc_g_t *mg = NULL; mg = to_mc_g_t(ta, s, un); + char *o_file = get_outfile_name(asm_opt.output_file_name); + trans_chain* t_ch = load_hc_trans(o_file); + mb_solve_core(&opt, mg, &(t_ch->k_trans), 1); mc_solve_core(&opt, mg, NULL); } else @@ -3699,6 +3714,8 @@ void clean_solve_general_ovlp(kv_u_trans_t *ta, uint32_t un, kv_gg_status *s) void mc_solve_general(kv_u_trans_t *ta, uint32_t un, kv_gg_status *s, uint16_t hapN, uint16_t update_ta, uint16_t write_dump) { + print_hap_p(s); + exit(1); mc_opt_t opt; mc_opt_init(&opt, asm_opt.n_perturb, asm_opt.f_perturb, asm_opt.seed); mc_gg_t *mg = init_mc_gg_t(un, s, hapN);