update k-mer peaks

Assembly.cpp

@@ -632,7 +632,7 @@ void Output_corrected_reads()
 
 void ha_overlap_and_correct(int round)
 {
-	int i, hom_cov;
+	int i, hom_cov, het_cov;
 	ha_ovec_buf_t **b;
 	ha_ecsave_buf_t *e;
 
@@ -640,7 +640,7 @@ void ha_overlap_and_correct(int round)
 	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));
-	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF, &hom_cov); // build the index
+	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF, &hom_cov, &het_cov); // build the index
 	if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
 		ha_opt_update_cov(&asm_opt, hom_cov);
 	if (asm_opt.required_read_name)
@@ -1141,20 +1141,55 @@ long long readNum, long long rescue_threshold, float cluster_threshold)
     fprintf(stderr, "[M::%s] took %0.2fs, revise edges #: %lld\n\n", __func__, Get_T()-startTime, revises);
 }
 
+void hap_recalculate_peaks(char* output_file_name)
+{
+    destory_read_bin(&R_INF);
+    destory_ma_hit_t_alloc(R_INF.paf);
+    destory_ma_hit_t_alloc(R_INF.reverse_paf);
+
+    char* gfa_name = (char*)malloc(strlen(output_file_name)+25);
+	sprintf(gfa_name, "%s.ec", output_file_name);
+
+    int hom_cov, het_cov;
+    // construct hash table for high occurrence k-mers
+    if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
+        ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov);
+        ha_opt_update_cov(&asm_opt, hom_cov);
+    }
+    free(R_INF.read_length);
+    free(R_INF.name_index);
+
+    load_All_reads(&R_INF, gfa_name);
+
+    ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index
+    asm_opt.hom_cov = hom_cov;
+    asm_opt.het_cov = het_cov;
+    ha_pt_destroy(ha_idx);
+    ha_idx = 0;
+
+    destory_read_bin(&R_INF);
+    free(gfa_name);
+
+    load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name); 
+    fprintf(stderr, "M::%s has done.\n", __func__);   
+}
+
 void ha_overlap_final(void)
 {
-	int i, hom_cov;
+	int i, hom_cov, het_cov;
 	ha_ovec_buf_t **b;
 	CALLOC(b, asm_opt.thread_num);
 	for (i = 0; i < asm_opt.thread_num; ++i)
 		b[i] = ha_ovec_init(1, 1);
-	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov); // build the index
+	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index
 	kt_for(asm_opt.thread_num, worker_ov_final, b, R_INF.total_reads);
 	ha_pt_destroy(ha_idx);
 	ha_idx = 0;
 	for (i = 0; i < asm_opt.thread_num; ++i)
 		ha_ovec_destroy(b[i]);
 	free(b);
+    asm_opt.hom_cov = hom_cov;
+    asm_opt.het_cov = het_cov;
 }
 
 int ha_assemble(void)
@@ -1172,6 +1207,7 @@ int ha_assemble(void)
         ///if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt);
 		if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();
 		if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
+        if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2;
 	}
 	if (!ovlp_loaded) {
 		// construct hash table for high occurrence k-mers
@@ -1201,6 +1237,7 @@ int ha_assemble(void)
 		if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
 		ha_triobin(&asm_opt);
 	}
+    if(ovlp_loaded == 2) ovlp_loaded = 0;
 
     build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf, 
         R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round, 

CommandLines.cpp

@@ -92,6 +92,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
 	asm_opt->max_ov_diff_ec = 0.04;
 	asm_opt->max_ov_diff_final = 0.03;
 	asm_opt->hom_cov = 20;
+    asm_opt->het_cov = -1024;
 	asm_opt->max_n_chain = 100;
     asm_opt->k_mer_min_freq = 3;
     asm_opt->k_mer_max_freq = 66;

CommandLines.h

@@ -38,6 +38,7 @@ typedef struct {
 	double max_ov_diff_ec;
 	double max_ov_diff_final;
 	int hom_cov;
+    int het_cov;
 	int max_n_chain; // fall-back max number of chains to consider
     int k_mer_min_freq;
     int k_mer_max_freq;

Overlaps.cpp

@@ -13489,6 +13489,7 @@ void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sou
 	write_ma_hit_ts(reverse_sources, RNF->total_reads, gfa_name);
 
 	free(gfa_name);
+    fprintf(stderr, "bin files have been written.\n");
 }
 
 int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, char* output_file_name)
@@ -26342,6 +26343,7 @@ long long bubble_dist, int read_graph, int write)
     asg_t *sg = NULL;
     ma_sub_t* coverage_cut = NULL;
 
+    
     // debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs", sources, reverse_sources, -1, "clean");
     
     ///actually min_thres = asm_opt.max_short_tip + 1 there are asm_opt.max_short_tip reads

Process_Read.cpp

@@ -98,6 +98,9 @@ void write_All_reads(All_reads* r, char* read_file_name)
 	fwrite(r->name, sizeof(char), r->total_name_length, fp);
 	fwrite(r->name_index, sizeof(uint64_t), r->name_index_size, fp);
 	fwrite(r->trio_flag, sizeof(uint8_t), r->total_reads, fp);
+	fwrite(&(asm_opt.hom_cov), sizeof(asm_opt.hom_cov), 1, fp);
+	fwrite(&(asm_opt.het_cov), sizeof(asm_opt.het_cov), 1, fp);
+
     free(index_name);    
 	fflush(fp);
     fclose(fp);
@@ -173,6 +176,8 @@ int load_All_reads(All_reads* r, char* read_file_name)
 	/****************************may have bugs********************************/
 	r->trio_flag = (uint8_t*)malloc(sizeof(uint8_t)*r->total_reads);
 	f_flag += fread(r->trio_flag, sizeof(uint8_t), r->total_reads, fp);
+	f_flag += fread(&(asm_opt.hom_cov), sizeof(asm_opt.hom_cov), 1, fp);
+    f_flag += fread(&(asm_opt.het_cov), sizeof(asm_opt.het_cov), 1, fp);
 	/****************************may have bugs********************************/
 
 	r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
@@ -195,6 +200,35 @@ int load_All_reads(All_reads* r, char* read_file_name)
 	return 1;
 }
 
+
+int destory_read_bin(All_reads* r)
+{
+
+	uint64_t i = 0;
+	for (i = 0; i < r->total_reads; i++)
+	{
+		if (r->N_site[i]) free(r->N_site[i]);
+		if (r->read_sperate[i]) free(r->read_sperate[i]);
+		if (r->cigars[i].record) free(r->cigars[i].record);
+		if (r->cigars[i].lost_base) free(r->cigars[i].lost_base);
+		if (r->second_round_cigar[i].record) free(r->second_round_cigar[i].record);
+		if (r->second_round_cigar[i].lost_base) free(r->second_round_cigar[i].lost_base);
+	}
+
+	free(r->N_site);
+	free(r->read_length);
+	free(r->read_size);
+	free(r->read_sperate);
+	free(r->name);
+	free(r->name_index);
+	free(r->trio_flag);
+	free(r->cigars);
+	free(r->second_round_cigar);
+	return 1;
+}
+
+
+
 void ha_insert_read_len(All_reads *r, int read_len, int name_len)
 {
 	r->total_reads++;

Process_Read.h

@@ -153,5 +153,6 @@ void reverse_complement(char* pattern, uint64_t length);
 void write_All_reads(All_reads* r, char* read_file_name);
 int load_All_reads(All_reads* r, char* read_file_name);
 void destory_All_reads(All_reads* r);
+int destory_read_bin(All_reads* r);
 
 #endif

Purge_Dups.cpp

@@ -173,7 +173,7 @@ void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo
         }
     }
 
-    print_peak(cov_buf, cov_buf_length, max_i);
+    ///print_peak(cov_buf, cov_buf_length, max_i);
 
     // look for smaller peak on the low end
 	max2 = -1; max2_i = -1;
@@ -190,7 +190,7 @@ void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo
 		}
 	}
 
-    fprintf(stderr, "***max2: %lld, max2_i: %lld\n", max2, max2_i);
+    ///fprintf(stderr, "***max2: %lld, max2_i: %lld\n", max2, max2_i);
 
     if (max2_i != -1 && max2_i > err_i && max2_i < max_i) 
     {
@@ -219,7 +219,7 @@ void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo
 		}
 	}
 
-    fprintf(stderr, "***max3: %lld, max3_i: %lld\n", max3, max3_i);
+    ///fprintf(stderr, "***max3: %lld, max3_i: %lld\n", max3, max3_i);
     
     //if found a peak
 	if (max3 != -1 && max3_i > max_i) 
@@ -245,7 +245,7 @@ void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo
             (*het_peak) = max_i;
         }
 
-        fprintf(stderr, "topo_peak: %lld, topo_peak_i: %lld\n", topo_peak, topo_peak_i);
+        ///fprintf(stderr, "topo_peak: %lld, topo_peak_i: %lld\n", topo_peak, topo_peak_i);
     }
     
 
@@ -262,10 +262,14 @@ void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo
         (*hom_peak) = max_i;
     }
 
-    fprintf(stderr, "max: %lld, max_i: %lld\n", max, max_i);
-    fprintf(stderr, "max2: %lld, max2_i: %lld\n", max2, max2_i);
-    fprintf(stderr, "max3: %lld, max3_i: %lld\n", max3, max3_i);
-    fprintf(stderr, "(*het_peak): %lld, (*hom_peak): %lld\n", (*het_peak), (*hom_peak));
+    // fprintf(stderr, "max: %lld, max_i: %lld\n", max, max_i);
+    // fprintf(stderr, "max2: %lld, max2_i: %lld\n", max2, max2_i);
+    // fprintf(stderr, "max3: %lld, max3_i: %lld\n", max3, max3_i);
+    fprintf(stderr, "[M::%s] Heterozygous k-mer peak: %d\n", __func__, asm_opt.het_cov);
+    fprintf(stderr, "[M::%s] Homozygous k-mer peak: %d\n", __func__, asm_opt.hom_cov);
+    fprintf(stderr, "[M::%s] Heterozygous coverage peak: %lld\n", __func__, (*het_peak));
+    fprintf(stderr, "[M::%s] Homozygous coverage peak: %lld\n", __func__, (*hom_peak));
+    fprintf(stderr, "[M::%s] Alter coverage peak: %lld\n", __func__, topo_peak_i);
 }
 
 
@@ -274,9 +278,12 @@ void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo
 long long get_alter_peak(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, uint64_t* position_index, 
 ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
 {
+    #define ALTER_COV_THRES 0.9
+    #define REAL_ALTER_THRES 0.1
+
     ma_utg_t* u = NULL;
     asg_t* nsg = ug->g;
-    uint64_t v, j, k, qn, n_vtx = nsg->n_seq;
+    uint64_t v, j, k, qn, n_vtx = nsg->n_seq, primary_bases = 0, alter_bases = 0;
     uint32_t tn, is_Unitig;
     long long* cov_buf = NULL;
     ma_hit_t *h;
@@ -295,6 +302,8 @@ ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
         {
             qn = u->a[k]>>33;
             position_index[qn] = 0;
+            R_bases = coverage_cut[qn].e - coverage_cut[qn].s;
+            primary_bases += R_bases;
         }
     }
 
@@ -305,6 +314,7 @@ ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
 
         C_bases = C_bases_primary = C_bases_alter = 0;
         R_bases = coverage_cut[qn].e - coverage_cut[qn].s;
+        alter_bases += R_bases;
         for (j = 0; j < (uint64_t)(sources[qn].length); j++)
         {
             h = &(sources[qn].buffer[j]);
@@ -335,16 +345,7 @@ ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
         // }
 
         C_bases = C_bases_primary + C_bases_alter;
-        if(C_bases_alter < C_bases * 0.8) continue;
-
-        ///fprintf(stderr, "******************qn: %lu\n",qn);
-        /**
-        if(qn == 1893151 || qn == 1929038)
-        {
-            fprintf(stderr, "******************qn: %lu\n",
-            qn);
-        }
-        **/
+        if(C_bases_alter < C_bases * ALTER_COV_THRES) continue;
 
         C_bases = C_bases/R_bases;
         if(C_bases < 0 || C_bases >= cov_buf_length) continue;
@@ -361,9 +362,9 @@ ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
         }
     }
 
-    fprintf(stderr, "alter max_i: %lld, max: %lld\n", max_i, max);
-
-    if(max_i < 5) max_i = max = -1;
+    ///fprintf(stderr, "alter max_i: %lld, max: %lld\n", max_i, max);
+    ///if(max_i < 5) max_i = max = -1;
+    if(alter_bases < primary_bases * REAL_ALTER_THRES) max_i = max = -1;
 
     free(cov_buf);
     memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);

htab.cpp

@@ -796,7 +796,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov)
 	return (void*)flt_tab;
 }
 
-ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov)
+ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov, int *het_cov)
 {
 	int64_t cnt[YAK_N_COUNTS], tot_cnt;
 	int peak_hom, peak_het, i, extra_flag1, extra_flag2;
@@ -818,6 +818,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 	fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]);
 	peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
 	if (hom_cov) *hom_cov = peak_hom;
+	if (het_cov) *het_cov = peak_het;
 	if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
 	if (flt_tab == 0) {
 		int cutoff = (int)(peak_hom * asm_opt->high_factor);

htab.h

@@ -30,7 +30,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov);
 int ha_ft_isflt(const void *hh, uint64_t y);
 void ha_ft_destroy(void *h);
 
-ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov);
+ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov, int *het_cov);
 void ha_pt_destroy(ha_pt_t *h);
 const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);