Merge pull request #420 from chhylp123/hifiasm_dev_debug

Hifiasm dev debug

Assembly.cpp

@@ -645,6 +645,44 @@ static void worker_ovec(void *data, long i, int tid)
 }
 
 
+static void worker_ovec_cal0(void *data, long i, int tid)
+{
+	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
+	int fully_cov, abnormal;
+    // if(i != 33) return;
+    // fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
+    // if (memcmp("m64012_190920_173625/88015004/ccs", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
+    //     fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i);
+    // } else {
+    //     return;
+    // }
+
+    ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 
+    0.02, asm_opt.max_n_chain, 1, NULL/**&(b->k_flag)**/, &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp));
+
+	clear_Cigar_record(&b->cigar1);
+	clear_Round2_alignment(&b->round2);
+
+	correct_overlap(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon,
+			&b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 0, &fully_cov, &abnormal);
+
+	b->num_read_base += b->self_read.length;
+	b->num_correct_base += b->correct.corrected_base;
+	b->num_recorrect_base += b->round2.dumy.corrected_base;
+
+	
+
+	R_INF.paf[i].is_fully_corrected = 0;
+	R_INF.paf[i].is_abnormal = abnormal;
+    R_INF.trio_flag[i] = AMBIGU;
+    
+    push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, &b->olist, 1);
+    push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
+
+    if(het_cnt) het_cnt[i] = get_het_cnt(&b->hap);
+}
+
+
 static void worker_ovec_related_reads(void *data, long i, int tid)
 {
 	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
@@ -982,6 +1020,45 @@ void ha_overlap_and_correct(int round)
     ///debug_print_pob_regions();
 }
 
+void ha_overlap_cal(int round)
+{
+	int i, hom_cov, het_cov;
+	ha_ovec_buf_t **b;
+    ha_flt_tab_hp = ha_idx_hp = NULL;
+
+	// overlap and correct reads
+	CALLOC(b, asm_opt.thread_num);
+	for (i = 0; i < asm_opt.thread_num; ++i)
+		b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1),0);
+    if(ha_idx) hom_cov = asm_opt.hom_cov;
+	if(ha_idx == NULL) ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
+	///debug_adapter(&asm_opt, &R_INF);
+    if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
+		ha_opt_update_cov(&asm_opt, hom_cov);
+    het_cnt = NULL;
+    if(round == asm_opt.number_of_round-1 && asm_opt.is_dbg_het_cnt) CALLOC(het_cnt, R_INF.total_reads);
+    // fprintf(stderr, "[M::%s-start]\n", __func__);
+	kt_for(asm_opt.thread_num, worker_ovec_cal0, b, R_INF.total_reads);///debug_for_fix
+    // fprintf(stderr, "[M::%s-end]\n", __func__);
+
+	ha_pt_destroy(ha_idx);
+	ha_idx = NULL;
+
+    if(het_cnt) {
+        print_het_cnt_log(het_cnt); free(het_cnt); het_cnt = NULL;
+    }
+
+	// collect statistics
+	for (i = 0; i < asm_opt.thread_num; ++i) {
+		asm_opt.num_bases += b[i]->num_read_base;
+		asm_opt.num_corrected_bases += b[i]->num_correct_base;
+		asm_opt.num_recorrected_bases += b[i]->num_recorrect_base;
+		asm_opt.mem_buf += ha_ovec_mem(b[i], NULL);
+		ha_ovec_destroy(b[i]);
+	}
+	free(b);
+}
+
 
 void update_overlaps(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf, 
 UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact)
@@ -1518,6 +1595,52 @@ void Output_PAF()
     fprintf(stderr, "PAF has been written.\n");
 }
 
+
+void Output_PAF0(ma_hit_t_alloc* sources, const char *prefix)
+{
+    fprintf(stderr, "Writing PAF to disk ...... \n");
+    char* paf_name = (char*)malloc(strlen(asm_opt.output_file_name)+strlen(prefix)+50);
+    sprintf(paf_name, "%s.%s.ovlp.paf", asm_opt.output_file_name, prefix);
+    FILE* output_file = fopen(paf_name, "w");
+    uint64_t i, j;
+
+    for (i = 0; i < R_INF.total_reads; i++)
+    {
+        for (j = 0; j < sources[i].length; j++)
+        {
+            fwrite(Get_NAME(R_INF, Get_qn(sources[i].buffer[j])), 1, 
+            Get_NAME_LENGTH(R_INF, Get_qn(sources[i].buffer[j])), output_file);
+            fwrite("\t", 1, 1, output_file);
+            fprintf(output_file, "%lu\t", (unsigned long)Get_READ_LENGTH(R_INF, Get_qn(sources[i].buffer[j])));
+            fprintf(output_file, "%d\t", Get_qs(sources[i].buffer[j]));
+            fprintf(output_file, "%d\t", Get_qe(sources[i].buffer[j]));
+            if(sources[i].buffer[j].rev)
+            {
+                fprintf(output_file, "-\t");
+            }
+            else
+            {
+                fprintf(output_file, "+\t");
+            }
+            fwrite(Get_NAME(R_INF, Get_tn(sources[i].buffer[j])), 1, 
+            Get_NAME_LENGTH(R_INF, Get_tn(sources[i].buffer[j])), output_file);
+            fwrite("\t", 1, 1, output_file);
+            fprintf(output_file, "%lu\t", (unsigned long)Get_READ_LENGTH(R_INF, Get_tn(sources[i].buffer[j])));
+            fprintf(output_file, "%d\t", Get_ts(sources[i].buffer[j]));
+            fprintf(output_file, "%d\t", Get_te(sources[i].buffer[j]));
+            fprintf(output_file, "%d\t", sources[i].buffer[j].ml);
+            fprintf(output_file, "%d\t", sources[i].buffer[j].bl);
+            fprintf(output_file, "255\n");
+
+        }
+    }
+
+    free(paf_name);
+    fclose(output_file);
+
+    fprintf(stderr, "PAF has been written.\n");
+}
+
 int check_cluster(uint64_t* list, long long listLen, ma_hit_t_alloc* paf, float threshold)
 {
     long long i, k;
@@ -1747,6 +1870,37 @@ void ug_idx_build(ma_ug_t *ug, int hap_n)
     exit(1);
 }
 
+int ha_assemble_ovec(void)
+{
+	extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
+	int r = 0, hom_cov = -1;
+	
+    ha_flt_tab = ha_idx = NULL;
+
+    // construct hash table for high occurrence k-mers
+    if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL) 
+    {
+        ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0);
+        ha_opt_update_cov(&asm_opt, hom_cov);
+    }
+    // error correction
+    assert(asm_opt.number_of_round > 0);
+    ha_opt_reset_to_round(&asm_opt, r); // this update asm_opt.roundID and a few other fields
+    ha_overlap_cal(r);
+    fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0);
+    
+    
+    ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads);
+    ha_ft_destroy(ha_flt_tab);
+
+    Output_PAF0(R_INF.paf, "0");
+    Output_PAF0(R_INF.reverse_paf, "1");
+    if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
+
+	destory_All_reads(&R_INF);
+	return 0;
+}
+
 int ha_assemble(void)
 {
     // debug_mc_g_t(MC_NAME);

Assembly.h

@@ -48,7 +48,7 @@ void ug_idx_build(ma_ug_t *ug, int hap_n);
 ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov, int is_ug);
 ha_ovec_buf_t *ha_ovec_buf_init(void *km, int is_final, int save_ov, int is_ug);
 void ha_ovec_destroy(ha_ovec_buf_t *b);
-int64_t ha_ovec_mem(const ha_ovec_buf_t *b, int64_t *mem_a)
+int64_t ha_ovec_mem(const ha_ovec_buf_t *b, int64_t *mem_a);
-;
+int ha_assemble_ovec(void);
 
 #endif

CommandLines.cpp

@@ -57,6 +57,7 @@ static ko_longopt_t long_options[] = {
     { "ul-round",     ko_required_argument, 342},
     { "prt-raw",     ko_no_argument, 343},
     { "integer-correct",     ko_required_argument, 344},
+    { "dbg-ovec",     ko_no_argument, 345},
 	{ 0, 0, 0 }
 };
 
@@ -269,6 +270,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->ul_clean_round = 1;
     asm_opt->prt_dbg_gfa = 0;
     asm_opt->integer_correct_round = 0;
+    asm_opt->dbg_ovec_cal = 0;
 }
 
 void destory_enzyme(enzyme* f)
@@ -804,6 +806,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 342) asm_opt->ul_clean_round = atol(opt.arg);
         else if (c == 343) asm_opt->prt_dbg_gfa = 1;
         else if (c == 344) asm_opt->integer_correct_round = atol(opt.arg);
+        else if (c == 345) asm_opt->dbg_ovec_cal = 1;
         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);
         }

CommandLines.h

@@ -5,7 +5,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.19.0-r554"
+#define HA_VERSION "0.19.0-r555"
 
 #define VERBOSE 0
 
@@ -141,6 +141,7 @@ typedef struct {
     int32_t ul_clean_round;
     int32_t prt_dbg_gfa;
     int32_t integer_correct_round;
+    uint8_t dbg_ovec_cal;
 } hifiasm_opt_t;
 
 extern hifiasm_opt_t asm_opt;

gfa_ut.cpp

@@ -16488,7 +16488,7 @@ ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt, uin
     u2g_clean(uidx, ulopt, keep_raw_utg, is_bridg);
     // renew_ul2_utg(uidx);
 
-    output_integer_graph(uidx, z->i_ug, asm_opt.output_file_name, 0);
+    // output_integer_graph(uidx, z->i_ug, asm_opt.output_file_name, 0);
     return z;
 }
 

main.cpp

@@ -61,8 +61,9 @@ int main(int argc, char *argv[])
 	// ed_band_cal_global((char*)"ACTTTTTT", 8, (char*)"AATTTT", 6, 3),
 	// ed_band_cal_global_128bit((char*)"ACTTTTTT", 8, (char*)"AATTTT", 6, 3));
 	// exit(1);
+	if(!(asm_opt.dbg_ovec_cal)) ret = ha_assemble();
+	else ret = ha_assemble_ovec();
 	
-	ret = ha_assemble();
     destory_opt(&asm_opt);
 	fprintf(stderr, "[M::%s] Version: %s\n", __func__, HA_VERSION);
 	fprintf(stderr, "[M::%s] CMD:", __func__);