bug fixing

CommandLines.cpp

@@ -23,7 +23,6 @@ static ko_longopt_t long_options[] = {
 	{ "ex-iter",       ko_required_argument, 308 },
     { "purge-cov",     ko_required_argument, 309 },
     { "pri-range",     ko_required_argument, 310 },
-    ///{ "high-het",      ko_no_argument, 311 },
     { "lowQ",          ko_required_argument, 312 },
 	{ "min-hist-cnt",  ko_required_argument, 313 },
     { "h1",      ko_required_argument, 314 },
@@ -32,7 +31,7 @@ static ko_longopt_t long_options[] = {
     { "b-cov",      ko_required_argument, 317 },
     { "h-cov",      ko_required_argument, 318 },
     { "m-rate",      ko_required_argument, 319 },
-    { "b-partition",      ko_no_argument, 320 },
+    { "primary",      ko_no_argument, 320 },
     { "t-occ",      ko_required_argument, 321 },
     { "seed",      ko_required_argument, 322 },
     { "n-perturb",      ko_required_argument, 323 },
@@ -82,6 +81,7 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "    --m-rate     FLOAT\n");
     fprintf(stderr, "                 break contigs at positions with <=FLOAT*coverage exact overlaps;\n");
     fprintf(stderr, "                 only work with '--b-cov' or '--h-cov'[%.2f]\n", asm_opt->m_rate);
+    fprintf(stderr, "    --primary    output a primary assembly and an alternate assembly\n");
     
 //	fprintf(stderr, "    --pri-range INT1[,INT2]\n");
 //	fprintf(stderr, "                keep contigs with coverage in this range in p_ctg.gfa; -1 to disable [auto,inf]\n");
@@ -126,7 +126,8 @@ void Print_H(hifiasm_opt_t* asm_opt)
 void init_opt(hifiasm_opt_t* asm_opt)
 {
 	memset(asm_opt, 0, sizeof(hifiasm_opt_t));
-	asm_opt->flag = 0;
+	///asm_opt->flag = 0;
+    asm_opt->flag = HA_F_PARTITION;
     asm_opt->coverage = -1;
     asm_opt->num_reads = 0;
     asm_opt->read_file_names = NULL;
@@ -483,7 +484,7 @@ int check_option(hifiasm_opt_t* asm_opt)
         fprintf(stderr, "[ERROR] [-s] must >= 0\n");
         return 0;
     }
-    
+
     return 1;
 }
 
@@ -644,7 +645,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 317) asm_opt->b_low_cov = atoi(opt.arg);
         else if (c == 318) asm_opt->b_high_cov = atoi(opt.arg);
         else if (c == 319) asm_opt->m_rate = atof(opt.arg);
-        else if (c == 320) asm_opt->flag |= HA_F_PARTITION;
+        else if (c == 320) asm_opt->flag -= HA_F_PARTITION;
         else if (c == 321) asm_opt->trio_flag_occ_thres = atoi(opt.arg);
         else if (c == 322) asm_opt->seed = atol(opt.arg);
         else if (c == 323) asm_opt->n_perturb = atoi(opt.arg);

Overlaps.cpp

@@ -12603,7 +12603,9 @@ bub_label_t* b_mask_t)
     
 
     hap_cov_t *cov = NULL;
-    trans_chain* t_ch = load_hc_hits(output_file_name);
+    trans_chain* t_ch = NULL;
+    if((asm_opt.flag & HA_F_VERBOSE_GFA)) t_ch = load_hc_hits(output_file_name);
+    
     if(!t_ch)
     {
         new_rtg_edges.a.n = 0;
@@ -12627,7 +12629,7 @@ bub_label_t* b_mask_t)
         ma_ug_print_bed(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, 
                 max_hang, min_ovlp, asm_opt.hic_inconsist_rate, NULL, NULL, cov);
 
-        write_trans_chain(cov->t_ch, output_file_name);
+        if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_trans_chain(cov->t_ch, output_file_name);
     }
 
     hic_analysis(ug, sg, cov?cov->t_ch:t_ch);
@@ -30000,11 +30002,13 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
 
     if (ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt))
     {
+        if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION;
         benchmark_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, &b_mask_t);
     }
     else if (ha_opt_triobin(&asm_opt))
     {   
+        if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION;
         output_trio_unitig_graph(sg, coverage_cut, o_file, FATHER, sources, 
         reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 
         0.05, 0.9, max_hang_length, mini_overlap_length, 0, &b_mask_t);
@@ -30014,16 +30018,18 @@ 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, &b_mask_t);
     }
-    else if(asm_opt.flag & HA_F_PARTITION)
+    else if((asm_opt.flag & HA_F_PARTITION) && (asm_opt.purge_level_primary > 0))
     {
         output_bp_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, &b_mask_t);
     }
     else
     {
+        if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION;
         output_contig_graph_primary(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, &b_mask_t);
 

hifiasm.1

@@ -235,7 +235,13 @@ Only work with
 and
 .B --h-cov.
 
-
+.TP
+.BI --primary
+Output a primary assembly and an alternate assembly.
+Hifiasm outputs two balanced assemblies and a primary
+assembly in default. Enable this option or 
+.B -l0
+outputs a primary assembly and an alternate assembly.
 
 
 .SS Trio-partition options
@@ -288,7 +294,7 @@ times in the other sample.
 Level of purge-dup. 0 to disable purge-dup, 1 to only purge contained haplotigs, 
 2 to purge all types of haplotigs, 3 to purge all types of haplotigs in most aggressive way 
 for high heterozygosity sample. 
-In default, [2] for non-trio assembly, [0] for trio assembly.
+In default, [3] for non-trio assembly, [0] for trio assembly.
 For trio assembly, only level 0 and level 1 are allowed.
 
 .TP
@@ -304,6 +310,10 @@ Min number of overlapped reads for duplicate haplotigs that should be purged [1]
 Coverage upper bound of Purge-dups, which is inferred automatically in default.
 If the coverage of a contig is higher than this bound, don't apply Purge-dups.
 
+.TP
+.BI   --n-hap \ INT
+Assumption of haplotype number.
+
 
 .SS Debugging options
 
@@ -316,12 +326,25 @@ Write additional files to speed up the debugging of graph cleaning.
 
 .TP
 .BI --h1 \ FILEs 
-File names of input Hi-C R1 [r1_1.fq,r1_2.fq,...]
+File names of input Hi-C R1 [r1_1.fq,r1_2.fq,...].
 
 .TP
 .BI --h2 \ FILEs 
-File names of input Hi-C R2 [r2_1.fq,r2_2.fq,...]
+File names of input Hi-C R2 [r2_1.fq,r2_2.fq,...].
 
+.TP
+.BI --n-perturb \ INT 
+Rounds of perturbation [50000]. Increasing this improves
+phasing results but takes longer time.
+
+.TP
+.BI --f-perturb \ FLOAT 
+Fraction to flip for perturbation [0.1]. Increasing this improves
+phasing results but takes longer time.
+
+.TP
+.BI --seed \ INT 
+RNG seed [11].
 
 .SH OUTPUTS
 

rcut.cpp

@@ -1305,20 +1305,25 @@ void mc_init_spin_all(const mc_opt_t *opt, mc_g_t *mg, mc_svaux_t *b)
 
 void mc_solve_core(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub)
 {
+	fprintf(stderr, "#######0#######\n");
 	double index_time = yak_realtime();
 	uint32_t st, i;
 	mc_svaux_t *b;
 	mc_bp_t *bp = NULL;
+	fprintf(stderr, "#######1#######\n");
 	mc_g_cc(mg->e);
-
+	fprintf(stderr, "#######2#######\n");
 	b = mc_svaux_init(mg, opt->seed);
+	fprintf(stderr, "#######3#######\n");
 	if(bub) bp = mc_bp_t_init(mg->e, b, bub, asm_opt.thread_num);
+	fprintf(stderr, "#######4#######\n");
 	/*******************************for debug************************************/
 	if(bp)
 	{
 		mc_init_spin_all(opt, mg, b);
 		mc_solve_bp(bp);
 	} 
+	fprintf(stderr, "#######5#######\n");
 	/*******************************for debug************************************/
 	// fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all(mg->e, b));
 	for (st = 0, i = 1; i <= mg->e->n_seq; ++i) {
@@ -1328,11 +1333,14 @@ void mc_solve_core(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub)
 		}
 	}
 	// fprintf(stderr, "##############end-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all(mg->e, b));
-
+	fprintf(stderr, "#######6#######\n");
 	if(bp) mc_solve_bp(bp);	
+	fprintf(stderr, "#######7#######\n");
 	///mc_write_info(g, b);
-	if(bp) mc_svaux_destroy(b);
+	mc_svaux_destroy(b);
-	destroy_mc_bp_t(&bp);
+	fprintf(stderr, "#######8#######\n");
+	if(bp) destroy_mc_bp_t(&bp);
+	fprintf(stderr, "#######9#######\n");
 	fprintf(stderr, "[M::%s::%.3f] ==> Partition\n", __func__, yak_realtime()-index_time);
 }
 
@@ -1442,11 +1450,16 @@ void p_nodes(mc_g_t *mg, trans_chain* t_ch, uint8_t* trio_flag)
 void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub)
 {
 	mc_opt_t opt;
+	fprintf(stderr, "*****0******\n");
 	mc_opt_init(&opt, asm_opt.n_perturb, asm_opt.f_perturb, asm_opt.seed);
+	fprintf(stderr, "*****1******\n");
 	mc_g_t *mg = init_mc_g_t(ug, read_g, s, renew_s);
+	fprintf(stderr, "*****2******\n");
 	update_mc_edges(mg, ovlp, ta, t_ch, f_rate, is_sys);
+	fprintf(stderr, "*****3******\n");
 	///debug_mc_g_t(mg);
 	mc_solve_core(&opt, mg, bub);
+	fprintf(stderr, "*****4******\n");
 
 	if((asm_opt.flag & HA_F_PARTITION) && t_ch)
 	{