add "n-weight"

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);
 }

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;

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", 
+    // 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);
+    // 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);
 

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

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);
+	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));
+		fprintf(stderr, "*********before-[M::%s::mc_score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
-	/*******************************for debug************************************/
+		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)