diff --git a/CommandLines.cpp b/CommandLines.cpp index 90f605a..2d27297 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -24,7 +24,8 @@ static ko_longopt_t long_options[] = { { "purge-cov", ko_required_argument, 309 }, { "pri-range", ko_required_argument, 310 }, { "high-het", ko_no_argument, 311 }, - { "lowQ", ko_required_argument, 312 }, + { "lowQ", ko_required_argument, 312 }, + { "min-hist-cnt", ko_required_argument, 313 }, { 0, 0, 0 } }; @@ -107,8 +108,7 @@ void init_opt(hifiasm_opt_t* asm_opt) 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; + asm_opt->min_hist_kmer_cnt = 5; asm_opt->load_index_from_disk = 1; asm_opt->write_index_to_disk = 1; asm_opt->number_of_round = 3; @@ -449,6 +449,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) } else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET; else if (c == 312) asm_opt->bed_inconsist_rate = atoi(opt.arg); + else if (c == 313) asm_opt->min_hist_kmer_cnt = atoi(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 b130949..86b351b 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -41,8 +41,7 @@ typedef struct { 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; + int min_hist_kmer_cnt; int load_index_from_disk; int write_index_to_disk; int number_of_round; diff --git a/hifiasm.1 b/hifiasm.1 index e30c661..2093c6f 100644 --- a/hifiasm.1 +++ b/hifiasm.1 @@ -111,6 +111,11 @@ assembly. .BI -r \ INT Rounds of haplotype-aware error corrections [3]. This option affects all outputs of hifiasm. +.TP +.BI --min-hist-cnt \ INT +When analyzing the k-mer spectrum, ignore counts below +.IR INT . + .SS Assembly options .TP diff --git a/hist.cpp b/hist.cpp index 9a1e31e..8e34806 100644 --- a/hist.cpp +++ b/hist.cpp @@ -1,4 +1,5 @@ #include +#include #include "htab.h" static void ha_hist_line(int c, int x, int exceed, int64_t cnt) @@ -11,16 +12,19 @@ static void ha_hist_line(int c, int x, int exceed, int64_t cnt) fprintf(stderr, " %lld\n", (long long)cnt); } -int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het) +int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het) { const int hist_max = 100; int i, start, low_i, max_i, max2_i, max3_i; int64_t max, max2, max3, min; - // find the low point from the left + // determine the start point + assert(n_cnt > start_cnt); *peak_het = -1; start = cnt[1] > 0? 1 : 2; - low_i = start; + + // find the low point from the left + low_i = start > start_cnt? start : start_cnt; for (i = low_i + 1; i < n_cnt; ++i) if (cnt[i] > cnt[i-1]) break; low_i = i - 1; diff --git a/htab.cpp b/htab.cpp index dcbf9c9..7082582 100644 --- a/htab.cpp +++ b/htab.cpp @@ -888,7 +888,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i if(!(ex_flag & HAF_SKIP_READ)) { ha_ct_hist(h, cnt, asm_opt->thread_num); - peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het); + peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, cnt, &peak_het); if (hom_cov) *hom_cov = peak_hom; if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); ///in default, asm_opt->high_factor = 5.0 @@ -925,7 +925,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f yak_realtime(), yak_cpu_usage(), (long)ct->tot); ha_ct_hist(ct, cnt, asm_opt->thread_num); 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); + peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, 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); @@ -1184,4 +1184,4 @@ int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_op free(gfa_name); return 1; -} \ No newline at end of file +} diff --git a/htab.h b/htab.h index c6feba6..bb53bb6 100644 --- a/htab.h +++ b/htab.h @@ -61,7 +61,7 @@ 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); void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct); -int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het); +int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het); static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function {