diff --git a/CommandLines.cpp b/CommandLines.cpp index 5aa19e4..71de953 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -37,6 +37,7 @@ static ko_longopt_t long_options[] = { { "n-perturb", ko_required_argument, 323 }, { "f-perturb", ko_required_argument, 324 }, { "n-hap", ko_required_argument, 325 }, + { "n-weight", ko_required_argument, 326 }, { 0, 0, 0 } }; @@ -113,6 +114,10 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " --h1 FILEs file names of Hi-C R1 [r1_1.fq,r1_2.fq,...]\n"); fprintf(stderr, " --h2 FILEs file names of Hi-C R2 [r2_1.fq,r2_2.fq,...]\n"); fprintf(stderr, " --seed INT RNG seed [%lu]\n", asm_opt->seed); + + + fprintf(stderr, " --n-weight INT\n"); + fprintf(stderr, " rounds of reweighting Hi-C links [%d]\n", asm_opt->n_weight); fprintf(stderr, " --n-perturb INT\n"); fprintf(stderr, " rounds of perturbation [%d]\n", asm_opt->n_perturb); fprintf(stderr, " --f-perturb FLOAT\n"); @@ -188,8 +193,9 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->polyploidy = 2; asm_opt->trio_flag_occ_thres = 60; asm_opt->seed = 11; - asm_opt->n_perturb = 50000; + asm_opt->n_perturb = 10000; asm_opt->f_perturb = 0.1; + asm_opt->n_weight = 3; } void destory_enzyme(enzyme* f) @@ -651,6 +657,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 323) asm_opt->n_perturb = atoi(opt.arg); else if (c == 324) asm_opt->f_perturb = atof(opt.arg); else if (c == 325) asm_opt->polyploidy = atoi(opt.arg); + else if (c == 326) asm_opt->n_weight = atoi(opt.arg); else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg); @@ -682,6 +689,5 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) get_queries(argc, argv, &opt, asm_opt); - return check_option(asm_opt); } diff --git a/CommandLines.h b/CommandLines.h index aac7329..f848f13 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -101,6 +101,7 @@ typedef struct { uint64_t seed; int32_t n_perturb; double f_perturb; + int32_t n_weight; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/hic.cpp b/hic.cpp index f674fbd..aefa275 100644 --- a/hic.cpp +++ b/hic.cpp @@ -2284,10 +2284,6 @@ uint32_t m_het_label, uint32_t p_het_label, uint32_t n_het_label) hom_occ++; } } - // if(uid == 18759) - // { - // fprintf(stderr, "het_occ: %u, hom_occ: %u\n", het_occ, hom_occ); - // } if((het_occ+hom_occ) == 0) return n_het_label; ///hom if((het_occ > ((het_occ+hom_occ)*0.8)) && ((het_occ+hom_occ) > m_het_occ)) return m_het_label; ///must het @@ -9271,7 +9267,7 @@ H_partition* hap, int8_t *s, trans_idx* dis) dis->med = ((buf.n&1)?buf.a[buf.n>>1]:((buf.a[buf.n>>1]+buf.a[(buf.n>>1)-1])/2)); } - fprintf(stderr, "dis->med: %lu\n", dis->med); + // fprintf(stderr, "dis->med: %lu\n", dis->med); kv_destroy(buf); return 1; @@ -13017,8 +13013,8 @@ void update_trans_g(ha_ug_index* idx, kv_u_trans_t *ta, bubble_type* bub) // append_boundary_chain_hic(idx->ug, ta, bub); - fprintf(stderr, "s_bub: %lu, f_bub: %lu, b_bub: %lu, b_end_bub: %lu, tangle_bub: %lu, cross_bub: %lu, mess_bub: %lu\n", - bub->s_bub, bub->f_bub, bub->b_bub, bub->b_end_bub, bub->tangle_bub, bub->cross_bub, bub->mess_bub); + // fprintf(stderr, "s_bub: %lu, f_bub: %lu, b_bub: %lu, b_end_bub: %lu, tangle_bub: %lu, cross_bub: %lu, mess_bub: %lu\n", + // bub->s_bub, bub->f_bub, bub->b_bub, bub->b_end_bub, bub->tangle_bub, bub->cross_bub, bub->mess_bub); ///reorder_bubbles(bub, ta, idx->ug->g->n_seq); // fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time); @@ -14206,6 +14202,18 @@ void destory_ps_t(ps_t **s) free((*s)); } +void verbose_het_stat(bubble_type *bub) +{ + uint64_t i, hetBase = 0, homBase = 0; + for (i = 0; i < bub->ug->g->n_seq; i++) + { + if(IF_HOM(i, *bub)) homBase += bub->ug->g->seq[i].len; + else hetBase += bub->ug->g->seq[i].len; + } + + fprintf(stderr, "[M::stat] # heterozygous bases: %lu; # homozygous bases: %lu\n", hetBase, homBase); +} + int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) { double index_time = yak_realtime(); @@ -14241,7 +14249,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) mb_nodes_t u; kv_init(u.bid); kv_init(u.idx); kv_init(u.u); memset(&bub, 0, sizeof(bubble_type)); - bub.round_id = 0; bub.n_round = 2; + bub.round_id = 0; bub.n_round = asm_opt.n_weight; for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++) { identify_bubbles(idx->ug, &bub, idx->t_ch->is_r_het, &(idx->t_ch->k_trans)); @@ -14267,6 +14275,8 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) **/ } + verbose_het_stat(&bub); + ///print_hc_links(&link, 0, &hap); // print_kv_u_trans(&k_trans, &link, s->s); diff --git a/hifiasm.1 b/hifiasm.1 index 75914e0..a1afd23 100644 --- a/hifiasm.1 +++ b/hifiasm.1 @@ -332,14 +332,19 @@ File names of input Hi-C R1 [r1_1.fq,r1_2.fq,...]. .BI --h2 \ FILEs File names of input Hi-C R2 [r2_1.fq,r2_2.fq,...]. +.TP +.BI --n-weight \ INT +Rounds of reweighting Hi-C links [3]. Increasing this may improves +phasing results but takes longer time. + .TP .BI --n-perturb \ INT -Rounds of perturbation [50000]. Increasing this improves +Rounds of perturbation [10000]. Increasing this may improves phasing results but takes longer time. .TP .BI --f-perturb \ FLOAT -Fraction to flip for perturbation [0.1]. Increasing this improves +Fraction to flip for perturbation [0.1]. Increasing this may improves phasing results but takes longer time. .TP diff --git a/rcut.cpp b/rcut.cpp index 66657a0..46ced5e 100644 --- a/rcut.cpp +++ b/rcut.cpp @@ -9,6 +9,8 @@ #include "kthread.h" #include "hic.h" +#define VERBOSE_CUT 0 + #define mc_edge_key(e) ((e).x) KRADIX_SORT_INIT(mce, mc_edge_t, mc_edge_key, member_size(mc_edge_t, x)) #define mb_edge_key(e) ((e).x) @@ -2565,7 +2567,11 @@ void mb_solve_core(mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t is_sys mb_g_t *mbg = init_mb_g_t(mg, ref, is_sys); mb_svaux_t *bb; /**************************init**************************/ - fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); + if(VERBOSE_CUT) + { + fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); + } + mc_svaux_t *b; mc_g_cc(mg->e); b = mc_svaux_init(mg, opt->seed); @@ -2576,11 +2582,17 @@ void mb_solve_core(mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t is_sys /**************************init**************************/ mb_g_cc(mbg); bb = mb_svaux_init(mbg, opt->seed); - /*******************************for debug************************************/ - // print_mb_g_blcok(mbg); - fprintf(stderr, "*********before-[M::%s::mc_score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); - fprintf(stderr, "*********before-[M::%s::mb_score->%f] ==> Partition\n", __func__, mb_score_all_advance(mg->e, mbg)); - /*******************************for debug************************************/ + + if(VERBOSE_CUT) + { + fprintf(stderr, "*********before-[M::%s::mc_score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); + fprintf(stderr, "*********before-[M::%s::mb_score->%f] ==> Partition\n", __func__, mb_score_all_advance(mg->e, mbg)); + /*******************************for debug************************************/ + // print_mb_g_blcok(mbg); + /*******************************for debug************************************/ + } + + opt->n_perturb = opt->n_b_perturb; for (st = 0, i = 1; i <= mbg->e->n_seq; ++i) { if (i == mbg->e->n_seq || mbg->e->cc[st]>>32 != mbg->e->cc[i]>>32) { @@ -2589,13 +2601,15 @@ void mb_solve_core(mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t is_sys } } opt->n_perturb = opt->n_s_perturb - opt->n_b_perturb; - /*******************************for debug************************************/ - // debug_mb_solve_core(mbg); - fprintf(stderr, "*********after-[M::%s::mc_score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); - fprintf(stderr, "*********after-[M::%s::mb_score->%f] ==> Partition\n", __func__, mb_score_all_advance(mg->e, mbg)); - /*******************************for debug************************************/ + mc_set_by_mbg(mg, mbg); - fprintf(stderr, "##############end-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); + if(VERBOSE_CUT) + { + /*******************************for debug************************************/ + // debug_mb_solve_core(mbg); + /*******************************for debug************************************/ + fprintf(stderr, "##############end-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); + } destory_mb_g_t(&mbg); mb_svaux_destroy(bb); fprintf(stderr, "[M::%s::%.3f] ==> Partition\n", __func__, yak_realtime()-index_time); @@ -2614,14 +2628,22 @@ void mc_solve_core(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub) if(bp) mc_init_spin_all(opt, mg, NULL, b); if(bp) mc_solve_bp(bp); /*******************************for debug************************************/ - fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); + if(VERBOSE_CUT) + { + fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); + } + for (st = 0, i = 1; i <= mg->e->n_seq; ++i) { if (i == mg->e->n_seq || mg->e->cc[st]>>32 != mg->e->cc[i]>>32) { mc_solve_cc(opt, mg, b, st, i - st); st = i; } } - fprintf(stderr, "##############end-[---M::%s::score->%f] ==> Partition\n", __func__, mc_score_all(mg->e, b)); + + if(VERBOSE_CUT) + { + fprintf(stderr, "##############end-[---M::%s::score->%f] ==> Partition\n", __func__, mc_score_all(mg->e, b)); + } if(bp) mc_solve_bp(bp); ///mc_write_info(g, b); mc_svaux_destroy(b); @@ -2796,7 +2818,7 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u mb_solve_core(&opt, mg, ref, is_sys); ///debug_mc_g_t(mg); - if(renew_s == 0) write_mc_g_t(&opt, mg, MC_NAME); + // if(renew_s == 0) write_mc_g_t(&opt, mg, MC_NAME); mc_solve_core(&opt, mg, bub); if((asm_opt.flag & HA_F_PARTITION) && t_ch)