diff --git a/CommandLines.cpp b/CommandLines.cpp index 70594b8..f39bc98 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -10,11 +10,13 @@ hifiasm_opt_t asm_opt; static ko_longopt_t long_options[] = { - { "version", ko_no_argument, 300 }, - { "dbg-gfa", ko_no_argument, 301 }, - { "write-paf", ko_no_argument, 302 }, - { "write-ec", ko_no_argument, 303 }, - { "skip-triobin", ko_no_argument, 304 }, + { "version", ko_no_argument, 300 }, + { "dbg-gfa", ko_no_argument, 301 }, + { "write-paf", ko_no_argument, 302 }, + { "write-ec", ko_no_argument, 303 }, + { "skip-triobin", ko_no_argument, 304 }, + { "max-od-ec", ko_no_argument, 305 }, + { "max-od-final", ko_no_argument, 306 }, { 0, 0, 0 } }; @@ -83,7 +85,9 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->k_mer_length = 51; asm_opt->mz_win = 51; asm_opt->bf_shift = 37; - asm_opt->high_factor = 5.0f; + asm_opt->high_factor = 5.0; + asm_opt->max_ov_diff_ec = 0.04; + asm_opt->max_ov_diff_final = 0.03; asm_opt->hom_cov = 20; asm_opt->max_n_chain = 100; asm_opt->k_mer_min_freq = 3; @@ -264,6 +268,15 @@ int check_option(hifiasm_opt_t* asm_opt) return 0; } + if (asm_opt->max_ov_diff_ec < asm_opt->max_ov_diff_final) { + fprintf(stderr, "[ERROR] max_ov_diff_ec shouldn't be smaller than max_ov_diff_final\n"); + return 0; + } + + if (asm_opt->max_ov_diff_ec < HA_MIN_OV_DIFF) { + fprintf(stderr, "[ERROR] max_ov_diff_ec shouldn't be smaller than %g\n", HA_MIN_OV_DIFF); + return 0; + } if(asm_opt->max_short_tip < 0) { @@ -385,6 +398,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 302) asm_opt->flag |= HA_F_WRITE_PAF; else if (c == 303) asm_opt->flag |= HA_F_WRITE_EC; else if (c == 304) asm_opt->flag |= HA_F_SKIP_TRIOBIN; + else if (c == 305) asm_opt->max_ov_diff_ec = atof(opt.arg); + else if (c == 306) asm_opt->max_ov_diff_final = atof(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); diff --git a/CommandLines.h b/CommandLines.h index 98cce08..1147c42 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,18 +3,20 @@ #include -#define HA_VERSION "0.5-dirty-r248" +#define HA_VERSION "0.5-dirty-r249" #define VERBOSE 0 -#define HA_F_NO_HPC 0x1 -#define HA_F_NO_KMER_FLT 0x2 -#define HA_F_VERBOSE_GFA 0x4 -#define HA_F_WRITE_EC 0x8 -#define HA_F_WRITE_PAF 0x10 -#define HA_F_SKIP_TRIOBIN 0x20 +#define HA_F_NO_HPC 0x1 +#define HA_F_NO_KMER_FLT 0x2 +#define HA_F_VERBOSE_GFA 0x4 +#define HA_F_WRITE_EC 0x8 +#define HA_F_WRITE_PAF 0x10 +#define HA_F_SKIP_TRIOBIN 0x20 #define HA_F_PURGE_CONTAIN 0x40 -#define HA_F_PURGE_JOIN 0x80 +#define HA_F_PURGE_JOIN 0x80 + +#define HA_MIN_OV_DIFF 0.02 // min sequence divergence in an overlap typedef struct { int flag; @@ -28,7 +30,9 @@ typedef struct { int k_mer_length; int mz_win; int bf_shift; - float high_factor; // coverage cutoff set to high_factor*hom_cov + double high_factor; // coverage cutoff set to high_factor*hom_cov + double max_ov_diff_ec; + double max_ov_diff_final; int hom_cov; int max_n_chain; // fall-back max number of chains to consider int k_mer_min_freq; diff --git a/Correct.cpp b/Correct.cpp index 1733e94..8fd285a 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -260,7 +260,7 @@ All_reads* R_INF) x_len = x_end - x_start + 1; - threshold = x_len * THRESHOLD_RATE; + threshold = x_len * asm_opt.max_ov_diff_ec; /****************************may have bugs********************************/ threshold = Adjust_Threshold(threshold, x_len); /****************************may have bugs********************************/ @@ -483,7 +483,7 @@ char* r_string) ///overlap length between [window_start, window_end] x_len = x_end - x_start + 1; - threshold = x_len * THRESHOLD_RATE; + threshold = x_len * asm_opt.max_ov_diff_ec; /****************************may have bugs********************************/ threshold = Adjust_Threshold(threshold, x_len); /****************************may have bugs********************************/ @@ -548,8 +548,8 @@ inline double trim_error_rate(overlap_region_alloc* overlap_list, long long ID) } else { - ///tError += (Adjust_Threshold(subWinLen*THRESHOLD_RATE, subWinLen) * 2); - tError += (Adjust_Threshold(subWinLen*THRESHOLD_RATE, subWinLen) * 3); + ///tError += (Adjust_Threshold(subWinLen*asm_opt.max_ov_diff_ec, subWinLen) * 2); + tError += Adjust_Threshold(subWinLen * asm_opt.max_ov_diff_ec, subWinLen) * 3; } } } @@ -566,8 +566,8 @@ inline double trim_error_rate(overlap_region_alloc* overlap_list, long long ID) } else { - ///tError += (Adjust_Threshold(subWinLen*THRESHOLD_RATE, subWinLen) * 2); - tError += (Adjust_Threshold(subWinLen*THRESHOLD_RATE, subWinLen) * 3); + ///tError += (Adjust_Threshold(subWinLen*asm_opt.max_ov_diff_ec, subWinLen) * 2); + tError += Adjust_Threshold(subWinLen * asm_opt.max_ov_diff_ec, subWinLen) * 3; } } } @@ -2118,7 +2118,7 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN x_end = overlap_list->list[j].w_list[i].x_end; x_len = x_end - x_start + 1; /****************************may have bugs********************************/ - ///threshold = x_len * THRESHOLD_RATE; + ///threshold = x_len * asm_opt.max_ov_diff_ec; threshold = overlap_list->list[j].w_list[i].error_threshold; /****************************may have bugs********************************/ /****************************may have bugs********************************/ @@ -2415,7 +2415,7 @@ UC_Read* g_read) // } - threshold = xLen * THRESHOLD_RATE; + threshold = xLen * asm_opt.max_ov_diff_ec; threshold = Adjust_Threshold(threshold, xLen); threshold = double_error_threshold(threshold, xLen); @@ -2931,7 +2931,7 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea x_end = overlap_list->list[j].w_list[i].x_end; x_len = x_end - x_start + 1; /****************************may have bugs********************************/ - ///threshold = x_len * THRESHOLD_RATE; + ///threshold = x_len * asm_opt.max_ov_diff_ec; threshold = overlap_list->list[j].w_list[i].error_threshold; /****************************may have bugs********************************/ /****************************may have bugs********************************/ @@ -2998,19 +2998,16 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea } } - ///error_rate = trim_error_rate(overlap_list, j); error_rate = non_trim_error_rate(overlap_list, j, R_INF, dumy, g_read); - ///if(error_rate <= 0.015) - ///if(error_rate <= 0.03) - if(error_rate <= FINAL_OVERLAP_ERROR_RATE) + if (error_rate <= asm_opt.max_ov_diff_final) { overlap_list->mapped_overlaps_length += overlap_length; overlap_list->list[j].is_match = 1; calculate_boundary_cigars(&(overlap_list->list[j]), R_INF, dumy, g_read); } - else if(error_rate <= 0.045) + else if (error_rate <= asm_opt.max_ov_diff_final * 1.5) { overlap_list->list[j].is_match = 3; } @@ -4552,7 +4549,7 @@ Cigar_record* current_cigar, long long uncorrected_window_start, Round2_alignmen x_start = corrected_window_start; x_end = corrected_window_end; x_len = x_end - x_start + 1; - threshold = x_len * THRESHOLD_RATE; + threshold = x_len * asm_opt.max_ov_diff_ec; /****************************may have bugs********************************/ threshold = Adjust_Threshold(threshold, x_len); /****************************may have bugs********************************/ @@ -4741,7 +4738,7 @@ void generate_consensus(overlap_region_alloc* overlap_list, All_reads* R_INF, Window_Pool w_inf; - init_Window_Pool(&w_inf, g_read->length, WINDOW, TAIL_LENGTH); + init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec)); int flag = 0; ///for last window @@ -4885,7 +4882,7 @@ All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy, Graph* g, int* abnormal) (*abnormal) = 0; Window_Pool w_inf; - init_Window_Pool(&w_inf, g_read->length, WINDOW, TAIL_LENGTH); + init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec)); int flag = 0; long long realLen = 0, tmpLen = 0; @@ -6594,7 +6591,7 @@ Correct_dumy* dumy) { long long window_start, window_end; Window_Pool w_inf; - init_Window_Pool(&w_inf, read_length, WINDOW, TAIL_LENGTH); + init_Window_Pool(&w_inf, read_length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec)); int flag = 0; long long realLen = 0, realLen_100 = 0; int to_recover = 0; @@ -7046,7 +7043,7 @@ void partition_overlaps(overlap_region_alloc* overlap_list, All_reads* R_INF, long long num_availiable_win = 0; Window_Pool w_inf; - init_Window_Pool(&w_inf, g_read->length, WINDOW, TAIL_LENGTH); + init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec)); int flag = 0; while(get_Window(&w_inf, &window_start, &window_end) && flag != -2) @@ -7135,7 +7132,7 @@ void partition_overlaps_advance(overlap_region_alloc* overlap_list, All_reads* R long long num_availiable_win = 0; Window_Pool w_inf; - init_Window_Pool(&w_inf, g_read->length, WINDOW, TAIL_LENGTH); + init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec)); int flag = 0; while(get_Window(&w_inf, &window_start, &window_end) && flag != -2) @@ -7219,7 +7216,7 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, Window_Pool w_inf; - init_Window_Pool(&w_inf, g_read->length, WINDOW, TAIL_LENGTH); + init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec)); int flag = 0; diff --git a/Correct.h b/Correct.h index 01f3c1c..944583a 100644 --- a/Correct.h +++ b/Correct.h @@ -16,6 +16,7 @@ #define INSERTION 2 #define DELETION 3 +#define WINDOW_MAX_SIZE (WINDOW + (int)(1.0 / HA_MIN_OV_DIFF) + 3) // TODO: why 1/max_ov_diff? ///#define FLAG_THRE 0 diff --git a/Hash_Table.h b/Hash_Table.h index e9bebe2..f7c3a3f 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -6,22 +6,13 @@ #define MAX_SUFFIX_BITS 64 #define MODE_VALUE 101 -///#define WINDOW 350 -///#define THRESHOLD 14 - #define WINDOW 375 -//#define WINDOW_BOUNDARY 150 #define WINDOW_BOUNDARY 375 ///for one side, the first or last WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY bases should not be corrected #define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25 #define THRESHOLD 15 -#define THRESHOLD_RATE 0.04 -#define TAIL_LENGTH int(1/THRESHOLD_RATE) -///#define OVERLAP_THRESHOLD 0.9 #define OVERLAP_THRESHOLD_FILTER 0.9 -#define WINDOW_MAX_SIZE WINDOW + TAIL_LENGTH + 3 #define THRESHOLD_MAX_SIZE 31 -#define FINAL_OVERLAP_ERROR_RATE 0.03 #define GROUP_SIZE 4 ///the max cigar likes 10M10D10M10D10M diff --git a/Overlaps.cpp b/Overlaps.cpp index 9361f46..5efb639 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -1725,10 +1725,9 @@ ma_sub_t* coverage_cut, float shift_rate) collect_contain(&(paf[i]), NULL, rLen, &max_left, &max_right, 0.1); ///collect_contain(&(paf[i]), &(rev_paf[i]), rLen, &max_left, &max_right, 0.1); - ////shift_rate should be (FINAL_OVERLAP_ERROR_RATE*2) + ////shift_rate should be (asm_opt.max_ov_diff_final*2) ///this read is a normal read - if(max_left.e > max_right.s && - (max_left.e - max_right.s >= rLen * shift_rate)) + if (max_left.e > max_right.s && (max_left.e - max_right.s >= rLen * shift_rate)) { continue; } @@ -26062,7 +26061,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) ///ma_hit_sub is just use to init coverage_cut, ///it seems we do not need ma_hit_cut & ma_hit_flt ma_hit_sub(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut); - detect_chimeric_reads(sources, n_read, readLen, coverage_cut, FINAL_OVERLAP_ERROR_RATE*2); + detect_chimeric_reads(sources, n_read, readLen, coverage_cut, asm_opt.max_ov_diff_final * 2.0); ma_hit_cut(sources, n_read, readLen, mini_overlap_length, &coverage_cut); ///print_binned_reads(sources, n_read, coverage_cut); ma_hit_flt(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);