added --min-hist-cnt for #49

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);

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;

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

hist.cpp

@@ -1,4 +1,5 @@
 #include <stdio.h>
+#include <assert.h>
 #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;

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;
-}
+}

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
 {