From f32bfc904a4034656aae44bce94ef7f7001f9166 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Fri, 10 Apr 2020 00:12:19 -0400 Subject: [PATCH] r215: prepare for list-based binning --- Assembly.cpp | 4 ++-- CommandLines.cpp | 21 +++++++++++---------- CommandLines.h | 11 ++++++++--- Overlaps.cpp | 16 ++++++++-------- Trio.cpp | 29 ++++++++++++++++++++--------- htab.h | 2 +- 6 files changed, 50 insertions(+), 33 deletions(-) diff --git a/Assembly.cpp b/Assembly.cpp index 577978b..716f60a 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1156,7 +1156,7 @@ int ha_assemble(void) if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) { ovlp_loaded = 1; fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage()); - if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) trio_partition(); + 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(); } @@ -1187,7 +1187,7 @@ int ha_assemble(void) ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads); ha_ft_destroy(ha_flt_tab); if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF(); - trio_partition(); + ha_triobin(&asm_opt); } 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, diff --git a/CommandLines.cpp b/CommandLines.cpp index fd4f849..4502785 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -61,14 +61,13 @@ 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->coverage = -1; asm_opt->num_reads = 0; asm_opt->read_file_names = NULL; asm_opt->output_file_name = (char*)(DEFAULT_OUTPUT); asm_opt->required_read_name = NULL; - asm_opt->pat_index = NULL; - asm_opt->mat_index = NULL; asm_opt->thread_num = 1; asm_opt->k_mer_length = 51; asm_opt->mz_win = 51; @@ -121,7 +120,7 @@ void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov) fprintf(stderr, "[M::%s] updated max_n_chain to %d\n", __func__, opt->max_n_chain); } -int check_file(char* name, const char* opt) +static int check_file(char* name, const char* opt) { if(!name) { @@ -256,8 +255,10 @@ int check_option(hifiasm_opt_t* asm_opt) } - if(asm_opt->pat_index != NULL && check_file(asm_opt->pat_index, "P") == 0) return 0; - if(asm_opt->mat_index != NULL && check_file(asm_opt->mat_index, "M") == 0) return 0; + if(asm_opt->fn_bin_yak[0] != NULL && check_file(asm_opt->fn_bin_yak[0], "YAK1") == 0) return 0; + if(asm_opt->fn_bin_yak[1] != NULL && check_file(asm_opt->fn_bin_yak[1], "YAK2") == 0) return 0; + if(asm_opt->fn_bin_list[0] != NULL && check_file(asm_opt->fn_bin_list[0], "LIST1") == 0) return 0; + if(asm_opt->fn_bin_list[1] != NULL && check_file(asm_opt->fn_bin_list[1], "LIST2") == 0) return 0; // fprintf(stderr, "input file num: %d\n", asm_opt->num_reads); // fprintf(stderr, "output file: %s\n", asm_opt->output_file_name); @@ -272,8 +273,6 @@ int check_option(hifiasm_opt_t* asm_opt) // fprintf(stderr, "size of popped large bubbles: %lld\n", asm_opt->large_pop_bubble_size); // fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip); // fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip); - // fprintf(stderr, "pat_index: %s\n", asm_opt->pat_index); - // fprintf(stderr, "mat_index: %s\n", asm_opt->mat_index); // fprintf(stderr, "min_cnt: %d\n", asm_opt->min_cnt); // fprintf(stderr, "mid_cnt: %d\n", asm_opt->mid_cnt); @@ -313,7 +312,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) int c; - while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:w:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:1:2:", long_options)) >= 0) { + while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:w:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:1:2:3:4:", long_options)) >= 0) { if (c == 'h') { Print_H(asm_opt); @@ -339,8 +338,10 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 'b') asm_opt->required_read_name = opt.arg; else if (c == 'c') asm_opt->min_cnt = atoi(opt.arg); else if (c == 'd') asm_opt->mid_cnt = atoi(opt.arg); - else if (c == '1' || c == 'P') asm_opt->pat_index = opt.arg; // -P/-M reserved for backward compatibility - else if (c == '2' || c == 'M') asm_opt->mat_index = opt.arg; + else if (c == '1' || c == 'P') asm_opt->fn_bin_yak[0] = opt.arg; // -P/-M reserved for backward compatibility + else if (c == '2' || c == 'M') asm_opt->fn_bin_yak[1] = opt.arg; + else if (c == '3') asm_opt->fn_bin_list[0] = opt.arg; + else if (c == '4') asm_opt->fn_bin_list[1] = opt.arg; else if (c == 'x') asm_opt->max_drop_rate = atof(opt.arg); else if (c == 'y') asm_opt->min_drop_rate = atof(opt.arg); else if (c == 'p') asm_opt->small_pop_bubble_size = atoll(opt.arg); diff --git a/CommandLines.h b/CommandLines.h index 7e251df..06f7202 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.3.0-dirty-r213" +#define HA_VERSION "0.3.0-dirty-r214" #define VERBOSE 0 @@ -20,8 +20,8 @@ typedef struct { char** read_file_names; char* output_file_name; char* required_read_name; - char* pat_index; - char* mat_index; + char *fn_bin_yak[2]; + char *fn_bin_list[2]; int thread_num; int k_mer_length; int mz_win; @@ -66,4 +66,9 @@ void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov); int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt); double Get_T(void); +static inline int ha_opt_triobin(const hifiasm_opt_t *opt) +{ + return ((opt->fn_bin_yak[0] && opt->fn_bin_yak[1]) || (opt->fn_bin_list[0] && opt->fn_bin_list[1])); +} + #endif diff --git a/Overlaps.cpp b/Overlaps.cpp index a715cf6..a8dca53 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -22786,7 +22786,7 @@ void pre_clean(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *sg, long ///remove isoloated single read tri_flag += asg_arc_del_single_node_directly(sg, asm_opt.max_short_tip, sources); - if (!((asm_opt.pat_index && asm_opt.mat_index))) + if (!ha_opt_triobin(&asm_opt)) { tri_flag += asg_arc_del_triangular_advance(sg, bubble_dist); ///remove the cross at the bubble carefully, just remove inexact cross @@ -24105,7 +24105,7 @@ kvec_t_u32_warp* new_rtg_nodes) int tri_flag = 0; ///remove very simple circle tri_flag += asg_arc_del_simple_circle_untig(sources, coverage_cut, r_g, 100, 0); - if (!((asm_opt.pat_index && asm_opt.mat_index))) + if (!ha_opt_triobin(&asm_opt)) { ///remove isoloated single read tri_flag += asg_arc_del_triangular_advance(r_g, bubble_dist); @@ -24376,7 +24376,7 @@ kvec_asg_arc_t_warp* new_rtg_edges) int tri_flag = 0; ///remove very simple circle tri_flag += asg_arc_del_simple_circle_untig(sources, coverage_cut, r_g, 100, 0); - if (!((asm_opt.pat_index && asm_opt.mat_index))) + if (!ha_opt_triobin(&asm_opt)) { ///remove isoloated single read tri_flag += asg_arc_del_triangular_advance(r_g, bubble_dist); @@ -24754,7 +24754,7 @@ uint32_t is_bubble_check, uint32_t is_primary_check) int tri_flag = 0; ///remove very simple circle tri_flag += asg_arc_del_simple_circle_untig(sources, coverage_cut, r_g, 100, 0); - if (!((asm_opt.pat_index && asm_opt.mat_index))) + if (!ha_opt_triobin(&asm_opt)) { ///remove isoloated single read tri_flag += asg_arc_del_triangular_advance(r_g, bubble_dist); @@ -26312,7 +26312,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) // -1, "clean"); - if ((asm_opt.pat_index && asm_opt.mat_index)) + if (ha_opt_triobin(&asm_opt)) { drop_edges_by_trio(sources, n_read); } @@ -26401,7 +26401,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) ///asg_arc_identify_simple_bubbles_multi(sg, 1); asg_arc_identify_simple_bubbles_multi(sg, 0); ///asg_arc_del_short_diploid_unclean_exact(sg, drop_ratio, sources); - if((asm_opt.pat_index && asm_opt.mat_index)) + if (ha_opt_triobin(&asm_opt)) { asg_arc_del_short_diploid_by_exact_trio(sg, asm_opt.max_short_tip, sources); } @@ -26414,7 +26414,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) asg_arc_identify_simple_bubbles_multi(sg, 1); - if((asm_opt.pat_index && asm_opt.mat_index)) + if (ha_opt_triobin(&asm_opt)) { asg_arc_del_short_diploid_by_length_trio(sg, drop_ratio, asm_opt.max_short_tip, reverse_sources, asm_opt.max_short_tip, 1, 1, 0, 0, ruIndex); @@ -26493,7 +26493,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) // rescue_no_coverage_aggressive(sg, sources, reverse_sources, &coverage_cut, ruIndex, max_hang_length, // mini_overlap_length, bubble_dist, 10); - if ((asm_opt.pat_index && asm_opt.mat_index)) + if (ha_opt_triobin(&asm_opt)) { char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); sprintf(buf, "%s.dip", output_file_name); diff --git a/Trio.cpp b/Trio.cpp index adef200..aeffeef 100644 --- a/Trio.cpp +++ b/Trio.cpp @@ -24,6 +24,9 @@ #define yak_ch_hash(a) ((a)>>YAK_COUNTER_BITS) KHASHL_SET_INIT(static klib_unused, yak_ht_t, yak_ht, uint64_t, yak_ch_hash, yak_ch_eq) +typedef const char *ha_cstr_t; +KHASHL_MAP_INIT(static klib_unused, cstr_ht_t, cstr_ht, ha_cstr_t, int64_t, kh_hash_str, kh_eq_str) + typedef struct { struct yak_ht_t *h; } yak_ch1_t; @@ -248,23 +251,18 @@ static void tb_worker(void *_data, long k, int tid) if(type == 'm') aux->seq->trio_flag[k] = MOTHER; } -void trio_partition(void) +static void ha_triobin_yak(const hifiasm_opt_t *opt) { - if (asm_opt.pat_index == NULL || asm_opt.mat_index == NULL) { - memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t)); - return; - } - yak_ch_t *ch; int i /**, min_cnt = 2, mid_cnt = 5**/; tb_shared_t aux; memset(&aux, 0, sizeof(tb_shared_t)); - aux.n_threads = asm_opt.thread_num, aux.print_diff = 0; + aux.n_threads = opt->thread_num, aux.print_diff = 0; aux.ratio_thres = 0.33; aux.seq = &R_INF; - ch = yak_ch_restore_core(0, asm_opt.pat_index, YAK_LOAD_TRIOBIN1, asm_opt.min_cnt, asm_opt.mid_cnt); - ch = yak_ch_restore_core(ch, asm_opt.mat_index, YAK_LOAD_TRIOBIN2, asm_opt.min_cnt, asm_opt.mid_cnt); + ch = yak_ch_restore_core(0, opt->fn_bin_yak[0], YAK_LOAD_TRIOBIN1, opt->min_cnt, opt->mid_cnt); + ch = yak_ch_restore_core(ch, opt->fn_bin_yak[1], YAK_LOAD_TRIOBIN2, opt->min_cnt, opt->mid_cnt); aux.k = ch->k; aux.ch = ch; @@ -285,3 +283,16 @@ void trio_partition(void) fprintf(stderr, "[M::%s::%.3f*%.2f] ==> partitioned reads using yak dumps\n", __func__, yak_realtime(), yak_cpu_usage()); } + +static void ha_triobin_list(const hifiasm_opt_t *opt) +{ +} + +void ha_triobin(const hifiasm_opt_t *opt) +{ + memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads * sizeof(uint8_t)); + if (opt->fn_bin_list[0] && opt->fn_bin_list[1]) + ha_triobin_list(opt); + else if (opt->fn_bin_yak[0] && opt->fn_bin_yak[1]) + ha_triobin_yak(opt); +} diff --git a/htab.h b/htab.h index 0355745..e0d4d0b 100644 --- a/htab.h +++ b/htab.h @@ -45,7 +45,7 @@ long yak_peakrss(void); double yak_peakrss_in_gb(void); double yak_cpu_usage(void); -void trio_partition(void); +void ha_triobin(const hifiasm_opt_t *opt); void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf); int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het);