-z version number

Assembly.cpp

@@ -1601,8 +1601,6 @@ void ha_overlap_final(void)
     asm_opt.het_cov = het_cov;
 }
 
-
-
 int ha_assemble(void)
 {
 	extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);

CommandLines.cpp

@@ -24,12 +24,7 @@ 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 },
-	{ "min-hist-cnt",  ko_required_argument, 313 },
-    { "h1",      ko_required_argument, 314 },
-    { "h2",      ko_required_argument, 315 },
-    { "enzyme",      ko_required_argument, 316 },
-    { "b-cov",      ko_required_argument, 317 },
+    { "lowQ",      ko_required_argument, 312 },
 	{ 0, 0, 0 }
 };
 
@@ -67,9 +62,6 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "    -u          disable post join contigs step which may improve N50\n");
     fprintf(stderr, "    --lowQ      INT\n");
     fprintf(stderr, "                output contig regions with >=INT%% inconsistency in BED format; 0 to disable [%d]\n", asm_opt->bed_inconsist_rate);
-    fprintf(stderr, "    --b-cov     INT\n");
-    fprintf(stderr, "                break contigs at breakpoints with coverage drop at <INT-fold coverage [%d]\n", asm_opt->break_cov);
-    
 //	fprintf(stderr, "    --pri-range INT1[,INT2]\n");
 //	fprintf(stderr, "                keep contigs with coverage in this range in p_ctg.gfa; -1 to disable [auto,inf]\n");
 
@@ -91,9 +83,6 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "                coverage upper bound of Purge-dups [auto]\n");
     fprintf(stderr, "    --high-het  enable this mode for high heterozygosity sample [experimental, not stable]\n");
     
-    fprintf(stderr, "  Hi-C-partition [experimental, not stable]:\n");
-    fprintf(stderr, "    --h1 FILEs   file names of Hi-C R1  [r1_1.fq,r1_2.fq,...]\n");
-    fprintf(stderr, "    --h2 FILEs   file names of Hi-C R2  [r2_1.fq,r2_2.fq,...]\n");
 
     fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
     fprintf(stderr, "See `man ./hifiasm.1' for detailed description of these command-line options.\n");
@@ -108,12 +97,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->read_file_names = NULL;
     asm_opt->output_file_name = (char*)(DEFAULT_OUTPUT);
     asm_opt->required_read_name = NULL;
-    asm_opt->hic_enzymes = NULL;
-    asm_opt->hic_reads[0] = NULL;
-    asm_opt->hic_reads[1] = NULL;
     asm_opt->thread_num = 1;
     asm_opt->k_mer_length = 51;
-    asm_opt->hic_mer_length = 31;
 	asm_opt->mz_win = 51;
 	asm_opt->bf_shift = 37;
 	asm_opt->high_factor = 5.0;
@@ -122,7 +107,8 @@ 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->min_hist_kmer_cnt = 5;
+    asm_opt->k_mer_min_freq = 3;
+    asm_opt->k_mer_max_freq = 66;
     asm_opt->load_index_from_disk = 1;
     asm_opt->write_index_to_disk = 1;
     asm_opt->number_of_round = 3;
@@ -143,40 +129,19 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->purge_level_primary = 2;
     asm_opt->purge_level_trio = 0;
     asm_opt->purge_simi_rate = 0.75;
-    asm_opt->purge_simi_rate_hic = 0.85;
     asm_opt->purge_overlap_len = 1;
-    asm_opt->purge_overlap_len_hic = 50;
     asm_opt->recover_atg_cov_min = -1024;
     asm_opt->recover_atg_cov_max = INT_MAX;
     asm_opt->hom_global_coverage = -1;
     asm_opt->bed_inconsist_rate = 70;
-    asm_opt->hic_inconsist_rate = 30;
-    ///asm_opt->bub_mer_length = 3;
-    asm_opt->bub_mer_length = 1000000;
-    asm_opt->break_cov = 0;
-}
-
-void destory_enzyme(enzyme* f)
-{
-    int i;
-    if(f != NULL)
-    {
-        for (i = 0; i < f->n; i++)
-        {
-            free(f->a[i]);
-        }
-        free(f->a);
-        free(f->l);
-        free(f);
-    }
 }
 
 void destory_opt(hifiasm_opt_t* asm_opt)
 {
-    if(asm_opt->read_file_names != NULL) free(asm_opt->read_file_names);
-    if(asm_opt->hic_enzymes != NULL) destory_enzyme(asm_opt->hic_enzymes);
-    if(asm_opt->hic_reads[0] != NULL) destory_enzyme(asm_opt->hic_reads[0]);
-    if(asm_opt->hic_reads[1] != NULL) destory_enzyme(asm_opt->hic_reads[1]);
+    if(asm_opt->read_file_names != NULL)
+    {
+        free(asm_opt->read_file_names);
+    }
 }
 
 void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round)
@@ -216,16 +181,6 @@ static int check_file(char* name, const char* opt)
     return 1;
 }
 
-static int check_hic_reads(enzyme* f, const char* opt)
-{
-    int i;
-    for (i = 0; i < f->n; i++)
-    {
-        if(check_file(f->a[i], opt) == 0) return 0;
-    }
-    return 1;
-}
-
 int check_option(hifiasm_opt_t* asm_opt)
 {
     if(asm_opt->read_file_names == NULL || asm_opt->num_reads == 0)
@@ -370,7 +325,7 @@ int check_option(hifiasm_opt_t* asm_opt)
 
     if(asm_opt->bed_inconsist_rate < 0 || asm_opt->bed_inconsist_rate > 100)
     {
-        fprintf(stderr, "[ERROR] inconsistency rate should be [0, 100] (--lowQ)\n");
+        fprintf(stderr, "[ERROR] inconsistency rate should be [0, 100] (--pb-range)\n");
         return 0;
     }
 
@@ -380,45 +335,6 @@ int check_option(hifiasm_opt_t* asm_opt)
     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;
     if(asm_opt->required_read_name != NULL && check_file(asm_opt->required_read_name, "b") == 0) return 0;
-    
-    if(asm_opt->hic_reads[0] != NULL && check_hic_reads(asm_opt->hic_reads[0], "HIC1") == 0) return 0;
-    if(asm_opt->hic_reads[1] != NULL && check_hic_reads(asm_opt->hic_reads[1], "HIC2") == 0) return 0;
-    if(asm_opt->hic_reads[0] != NULL && asm_opt->hic_reads[1] == NULL)
-    {
-        fprintf(stderr, "[ERROR] lack r2 of HiC reads (--h2)\n");
-        return 0;
-    }
-    if(asm_opt->hic_reads[1] != NULL && asm_opt->hic_reads[0] == NULL)
-    {
-        fprintf(stderr, "[ERROR] lack r1 of HiC reads (--h1)\n");
-        return 0;
-    }
-
-    if(asm_opt->hic_reads[0] != NULL && asm_opt->hic_reads[1] != NULL && 
-                            asm_opt->hic_reads[0]->n != asm_opt->hic_reads[1]->n)
-    {
-        fprintf(stderr, "[ERROR] wrong r1 and r2 of HiC reads (--h1 && --h2)\n");
-        return 0;
-    }
-
-    if(asm_opt->hic_enzymes != NULL && asm_opt->hic_enzymes->n == 0)
-    {
-        fprintf(stderr, "[ERROR] wrong HiC enzymes (--enzyme)\n");
-        return 0;
-    }
-
-    if(asm_opt->hic_reads[0] != NULL && asm_opt->hic_reads[0]->n == 0)
-    {
-        fprintf(stderr, "[ERROR] wrong r1 of HiC reads (--h1)\n");
-        return 0;
-    }
-
-    if(asm_opt->hic_reads[1] != NULL && asm_opt->hic_reads[1]->n == 0)
-    {
-        fprintf(stderr, "[ERROR] wrong r2 of HiC reads (--h2)\n");
-        return 0;
-    }
-
     // fprintf(stderr, "input file num: %d\n", asm_opt->num_reads);
     // fprintf(stderr, "output file: %s\n", asm_opt->output_file_name);
     // fprintf(stderr, "number of threads: %d\n", asm_opt->thread_num);
@@ -468,61 +384,6 @@ void get_queries(int argc, char *argv[], ketopt_t* opt, hifiasm_opt_t* asm_opt)
     }
 }
 
-void get_hic_enzymes(char *argv, enzyme** x, int check_name)
-{
-    int i, k, pre_i, len = strlen(argv);
-    (*x) = (enzyme*)calloc(1, sizeof(enzyme));
-    if(len == 0)
-    {
-        (*x)->n = 0; (*x)->l = NULL; (*x)->a = NULL;
-        return;
-    }
-
-
-    (*x)->n = 1;
-    for (i = pre_i = 0; i < len; i++)
-    {
-        if(argv[i] == ',')
-        {
-            (*x)->n++;
-            continue;
-        } 
-        
-        if(check_name)
-        {
-            if(argv[i] != 'A' && argv[i] != 'C' && argv[i] != 'G' && argv[i] != 'T' &&
-                argv[i] != 'a' && argv[i] != 'c' && argv[i] != 'g' && argv[i] != 't' && 
-                argv[i] != 'N' && argv[i] != 'n')
-            {
-                (*x)->n = 0;
-                (*x)->l = NULL;
-                (*x)->a = NULL;
-                return;
-            }
-        }
-        
-    }
-    (*x)->l = (int*)calloc((*x)->n, sizeof(int));
-    (*x)->a = (char**)calloc((*x)->n, sizeof(char*));
-
-    for (i = pre_i = k = 0; i < len; i++)
-    {
-        if(argv[i] == ',')
-        {
-            (*x)->l[k] = i - pre_i;
-            (*x)->a[k] = (char*)malloc(sizeof(char)*((*x)->l[k]+1));
-            memcpy((*x)->a[k], argv + pre_i, (*x)->l[k]);
-            (*x)->a[k][(*x)->l[k]] = '\0';
-            pre_i = i + 1;
-            k++;
-        }
-    }
-
-    (*x)->l[k] = i - pre_i;
-    (*x)->a[k] = (char*)malloc(sizeof(char)*((*x)->l[k]+1));
-    memcpy((*x)->a[k], argv + pre_i, (*x)->l[k]);
-    (*x)->a[k][(*x)->l[k]] = '\0';
-}
 
 int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 {
@@ -588,11 +449,6 @@ 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 == 314) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[0]), 0);
-        else if (c == 315) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[1]), 0);
-        else if (c == 316) get_hic_enzymes(opt.arg, &(asm_opt->hic_enzymes), 1);
-        else if (c == 317) asm_opt->break_cov = 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);
@@ -611,7 +467,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 		}
     }
 
-    
+
     if (argc == opt.ind)
     {
         Print_H(asm_opt);
@@ -620,7 +476,5 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 
     get_queries(argc, argv, &opt, asm_opt);
 
-
-
     return check_option(asm_opt);
 }

CommandLines.h

@@ -3,7 +3,7 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.14-r309"
+#define HA_VERSION "0.13-r308-z"
 
 #define VERBOSE 0
 
@@ -21,11 +21,6 @@
 
 #define HA_MIN_OV_DIFF       0.02 // min sequence divergence in an overlap
 
-typedef struct{
-    int *l, n; 
-    char **a;
-}enzyme;
-
 typedef struct {
 	int flag;
     int num_reads;
@@ -35,13 +30,9 @@ typedef struct {
 	char *fn_bin_yak[2];
 	char *fn_bin_list[2];
 	char *extract_list;
-    enzyme *hic_reads[2];
-    enzyme *hic_enzymes;
 	int extract_iter;
     int thread_num;
     int k_mer_length;
-    int hic_mer_length;
-    int bub_mer_length;
 	int mz_win;
 	int bf_shift;
 	double high_factor; // coverage cutoff set to high_factor*hom_cov
@@ -49,9 +40,9 @@ typedef struct {
 	double max_ov_diff_final;
 	int hom_cov;
     int het_cov;
-    int break_cov;
 	int max_n_chain; // fall-back max number of chains to consider
-	int min_hist_kmer_cnt;
+    int k_mer_min_freq;
+    int k_mer_max_freq;
     int load_index_from_disk;
     int write_index_to_disk;
     int number_of_round;
@@ -68,18 +59,15 @@ typedef struct {
     int purge_level_primary;
     int purge_level_trio;
     int purge_overlap_len;
-    int purge_overlap_len_hic;
     int recover_atg_cov_min;
     int recover_atg_cov_max;
     int hom_global_coverage;
     int bed_inconsist_rate;
-    int hic_inconsist_rate;
 
     float max_hang_rate;
     float min_drop_rate;
     float max_drop_rate;
     float purge_simi_rate;
-    float purge_simi_rate_hic;
 
     long long small_pop_bubble_size;
     long long large_pop_bubble_size;
@@ -105,9 +93,4 @@ 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]));
 }
 
-static inline int ha_opt_hic(const hifiasm_opt_t *opt)
-{
-    return ((opt->hic_reads[0] && opt->hic_reads[1]));
-}
-
 #endif

Correct.h

@@ -21,8 +21,8 @@
 
 ///#define FLAG_THRE 0
 
-#define MAX(x, y) (((x) >= (y))?(x):(y))
-#define MIN(x, y) (((x) <= (y))?(x):(y))
+#define MAX(x, y) ((x >= y)?(x):(y))
+#define MIN(x, y) ((x <= y)?(x):(y))
 #define DIFF(x, y) ((MAX((x), (y))) - (MIN((x), (y))))
 #define OVERLAP(x_start, x_end, y_start, y_end) (MIN(x_end, y_end) - MAX(x_start, y_start) + 1) 
 ///#define OVERLAP(x_start, x_end, y_start, y_end) MIN(x_end, y_end) - MAX(x_start, y_start) + 1 

Makefile

@@ -6,7 +6,7 @@ CPPFLAGS=
 INCLUDES=
 OBJS=		CommandLines.o Process_Read.o Assembly.o Hash_Table.o \
 			POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \
-			htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o
+			htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm
 
@@ -71,4 +71,3 @@ main.o: CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h Assembly.h
 main.o: Levenshtein_distance.h htab.h
 sketch.o: kvec.h htab.h Process_Read.h Overlaps.h kdq.h CommandLines.h
 sys.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
-hic.o: hic.h

Overlaps.h

@@ -21,7 +21,6 @@
 #define PRIMARY_LABLE 0
 #define ALTER_LABLE 1
 #define HAP_LABLE 2
-#define FAKE_LABLE 4
 #define TRIO_THRES 0.9
 #define DOUBLE_CHECK_THRES 0.1
 #define FINAL_DOUBLE_CHECK_THRES 0.2
@@ -107,11 +106,6 @@ typedef struct {
 	uint8_t no_l_indel;
 } asg_arc_t;
 
-typedef struct {
-	size_t n, m;
-	asg_arc_t* a;
-} kv_asg_arc_t;
-
 
 typedef struct {
 	uint32_t len:31, circ:1; // len: length of the unitig; circ: circular if non-zero
@@ -161,7 +155,6 @@ typedef struct { size_t n, m; ma_utg_t *a; } ma_utg_v;
 typedef struct {
 	ma_utg_v u;
 	asg_t *g;
-	kvec_t(uint64_t) occ;
 } ma_ug_t;
 
 typedef struct {
@@ -401,8 +394,7 @@ typedef struct {
 }kvec_asg_arc_t_warp;
 
 void sort_kvec_t_u64_warp(kvec_t_u64_warp* u_vecs, uint32_t is_descend);
-int asg_arc_del_multi(asg_t *g);
-int asg_arc_del_asymm(asg_t *g);
+
 
 typedef struct {
 	uint32_t q_pos;
@@ -450,8 +442,11 @@ long long max_hang_length, long long clean_round, long long gap_fuzz,
 float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, 
 long long bubble_dist, int read_graph, int write);
 
-void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int id, char* command);
+void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources, 
+ma_hit_t_alloc* reverse_sources, int id, char* command);
+
 void collect_abnormal_edges(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long readNum);
+
 void add_overlaps(ma_hit_t_alloc* source_paf, ma_hit_t_alloc* dest_paf, uint64_t* source_index, long long listLen);
 void remove_overlaps(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen);
 void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, 
@@ -460,6 +455,18 @@ uint64_t* source_index, long long listLen);
 #define EvaluateLen(U, id) ((U).a[(id)].start)
 #define IsMerge(U, id) ((U).a[(id)].end)
 #define kv_reuse(v, rn, rm, r) ((v).n = (rn), (v).m = (rm), (v).a = (r))
+#define long_tip(U, id, threshold) ((EvaluateLen((U), (id))>=(threshold))&&(!((U).a[(id)].circ)))
+///there are threee cases: 
+///1. if this untig is too long (>maxShortUntig), it must be not short untig/must be a long untig
+///2. if this untig is long (>minLongUntig && EvaluateLen(ug->u, av[i].v>>1) > (EvaluateLen(ug->u, v>>1)*l_untig_rate)), it might be a long tip
+#define check_long_tip(U, id, minLongUntig, maxShortUntig, ShortUntigRate, mainLen) \
+                        ((!((U).a[(id)].circ)) \
+                        && \
+                        ((EvaluateLen((U), (id)) > (maxShortUntig))\
+                        ||\
+                        ((long_tip((U), (id), (minLongUntig)))\
+                        &&\
+                        (EvaluateLen((U), (id)) > (ShortUntigRate)*(mainLen)))))
 #define Get_vis(visit, v, d) (((visit)[(v)>>1])&(((((v)<<(d))&1)+1)))
 #define Set_vis(visit, v, d) (((visit)[(v)>>1])|=(((((v)<<(d))&1)+1)))
 
@@ -474,10 +481,14 @@ typedef struct {
 
 void init_R_to_U(R_to_U* x, uint64_t len);
 void destory_R_to_U(R_to_U* x);
-void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig, uint8_t* flag);
+void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig);
 void get_R_to_U(R_to_U* x, uint32_t rID, uint32_t* uID, uint32_t* is_Unitig);
 void transfor_R_to_U(R_to_U* x);
 void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, int require_equal_nv, int test_tangle);
+void clean_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources,
+long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
+R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen, 
+uint32_t miniBiGraph, float chimeric_rate, int is_final_clean);
 int asg_pop_bubble_primary(asg_t *g, int max_dist);
 long long asg_arc_del_simple_circle_untig(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *g, long long circleLen, int is_drop);
 
@@ -732,12 +743,13 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 	stops_threshold, b_0) == LOOP)
 	{
 		return UNAVAILABLE;
-	}
+	} 
 	if(get_unitig(nsg, ug, v_1, &vEnd, &ELen_1, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 
 	stops_threshold, b_1) == LOOP)
 	{
 		return UNAVAILABLE;
 	}
+
 	if(ELen_0<=min_edge_length || ELen_1<=min_edge_length) return UNAVAILABLE;
 
 	rIdContig b_max, b_min;
@@ -758,6 +770,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 
 	uint32_t max_count = 0, min_count = 0;
 	ma_utg_t *node_min = NULL, *node_max = NULL;
+
 	if(ug != NULL)
 	{
 		/*****************************label all unitigs****************************************/
@@ -768,11 +781,12 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 			for (b_max.readI = 0; b_max.readI < node_max->n; b_max.readI++)
 			{
 				qn = (node_max->a[b_max.readI]>>33);
-				set_R_to_U(ruIndex, qn, (b_max.b_0->b.a[b_max.untigI]>>1), 1, &(read_sg->seq[qn].c));
+				set_R_to_U(ruIndex, qn, (b_max.b_0->b.a[b_max.untigI]>>1), 1);
 			}
 		}
 		/*****************************label all unitigs****************************************/
 
+
 		///each unitig
 		for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
 		{
@@ -806,6 +820,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 				}
 			}
 		}
+
 		/*****************************label all unitigs****************************************/
 		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
 		{
@@ -818,6 +833,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 			}
 		}
 		/*****************************label all unitigs****************************************/
+		
 	}
 	else
 	{
@@ -825,7 +841,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 		for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++)
 		{
 			qn = (b_max.b_0->b.a[b_max.untigI]>>1);
-			set_R_to_U(ruIndex, qn, 1, 1, &(read_sg->seq[qn].c));
+			set_R_to_U(ruIndex, qn, 1, 1);
 		}
 		/*****************************label all reads****************************************/
 
@@ -880,6 +896,8 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 	return NON_PLOID;
 }
 
+
+
 inline uint32_t check_different_haps_naive(asg_t *nsg, ma_ug_t *ug, asg_t *read_sg, 
 uint32_t v_0, uint32_t v_1, ma_hit_t_alloc* reverse_sources, buf_t* b_0, buf_t* b_1, 
 R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
@@ -1029,6 +1047,8 @@ uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp*
 void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
 R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t is_bubble_check, 
 uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges);
+void deduplicate(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig, 
+long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t resolve_tangle);
 void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, 
 ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate);
 void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
@@ -1041,77 +1061,9 @@ uint32_t is_bubble_check, uint32_t is_primary_check);
 uint32_t get_edge_from_source(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, 
 R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t query, uint32_t target, asg_arc_t* t);
 uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, int max_dist, buf_t *b, 
-uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes);
+uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop);
 int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio);
 
-
-typedef struct{
-    double weight;
-    uint32_t uID:31, del:1;
-    uint64_t dis;
-	uint64_t occ;
-	///uint64_t occ:63, scaff:1;
-    ///uint32_t enzyme;
-} hc_edge;
-
-typedef struct{
-    kvec_t(hc_edge) e;
-    kvec_t(hc_edge) f;//forbiden
-} hc_linkeage;
-
-typedef struct{
-	uint64_t beg, end;
-}bed_interval;
-
-typedef struct{
-	size_t n, m;
-	bed_interval* a;
-}bed_in;
-
-typedef struct{
-    kvec_t(hc_linkeage) a;
-    kvec_t(uint64_t) enzymes;
-	kvec_t(bed_in) bed;
-    uint32_t* u_idx;
-	uint64_t r_num;
-} hc_links;
-
-typedef struct{
-    ///kvec_t(hc_edge) a;
-    size_t n, m; 
-    hc_edge *a;
-}hc_edge_warp;
-
-void init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num);
-void destory_hc_links(hc_links* link);
-void clean_primary_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources,
-long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
-R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen, 
-uint32_t miniBiGraph, float chimeric_rate, int is_final_clean, int just_bubble_pop, 
-float drop_ratio, hc_links* link);
-void adjust_utg_by_primary(ma_ug_t **ug, asg_t* read_g, float drop_rate,
-ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, 
-long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
-R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp,
-kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link);
-void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug, asg_t *read_sg);
-ma_ug_t* copy_untig_graph(ma_ug_t *src);
-ma_ug_t* output_trio_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, 
-uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist, 
-long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, 
-float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench);
-asg_t* copy_read_graph(asg_t *src);
-ma_ug_t *ma_ug_gen(asg_t *g);
-void ma_ug_destroy(ma_ug_t *ug);
-
-inline int inter_interval(int a_s, int a_e, int b_s, int b_e, int* i_s, int* i_e)
-{
-    if(a_s > b_e || b_s > a_e) return 0;
-    if(i_s) (*i_s) = a_s >= b_s? a_s : b_s; ///MAX(a_s, b_s);
-    if(i_e) (*i_e) = a_e <= b_e? a_e : b_e; ///MIN(a_e, b_e);
-    return 1;
-}
-
 #define JUNK_COV 5
 #define DISCARD_RATE 0.8
 

Purge_Dups.cpp

@@ -7,7 +7,6 @@
 #include "Correct.h"
 #include "kthread.h"
 #include "kdq.h"
-#include "hic.h"
 
 KDQ_INIT(uint64_t)
 
@@ -3292,7 +3291,7 @@ int asg_pop_bubble_purge_graph(asg_t *purge_g, int max_dist)
 		for (i = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs
 			if (!av[i].del) ++n_arc;
 		if (n_arc > 1)
-			n_pop += asg_bub_pop1_primary_trio(purge_g, NULL, v, max_dist, &b, (uint32_t)-1, DROP, 1, NULL, NULL);
+			n_pop += asg_bub_pop1_primary_trio(purge_g, NULL, v, max_dist, &b, (uint32_t)-1, DROP, 1);
 	}
 	free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
 	if (n_pop) asg_cleanup(purge_g);
@@ -3926,57 +3925,10 @@ kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edg
     
 }
 
-
-
-void collect_reverse_unitig_pair(hc_links* link, ma_ug_t *ug, hap_overlaps* t)
-{
-    uint32_t i = 0, k = 0, rId_0, rId_1, pre_0, pre_1, b_0 = t->xUid, b_1 = t->yUid;
-    uint64_t d = RC_2;
-    ma_utg_t* u_b_0 = &(ug->u.a[b_0]);
-    ma_utg_t* u_b_1 = &(ug->u.a[b_1]);
-    if(u_b_0->n == 0) return;
-    if(u_b_1->n == 0) return;
-
-    for (i = t->x_beg_id, pre_0 = (uint32_t)-1; i < t->x_end_id; i++)
-    {
-        rId_0 = u_b_0->a[i]>>33;
-        if(link->u_idx[rId_0] == (uint32_t)-1) continue;
-        if(pre_0 == link->u_idx[rId_0]) continue;
-        pre_0 = link->u_idx[rId_0];
-        
-        for (k = t->y_beg_id, pre_1  = (uint32_t)-1; k < t->y_end_id; k++)
-        {
-            rId_1 = u_b_1->a[k]>>33;
-            if(link->u_idx[rId_1] == (uint32_t)-1) continue;
-            if(pre_1 == link->u_idx[rId_1]) continue;
-            pre_1 = link->u_idx[rId_1];
-            push_hc_edge(&(link->a.a[pre_0]), pre_1, 1, 1, &d);
-            push_hc_edge(&(link->a.a[pre_1]), pre_0, 1, 1, &d);
-        }
-    }
-
-}
-
-
-void collect_reverse_unitigs_purge(buf_t* b_0, hc_links* link, ma_ug_t *ug, hap_overlaps_list* all_ovlp)
-{
-    if(b_0->b.n <= 1) return;
-    uint32_t k;
-    int index = 0;
-    for (k = 0; k < b_0->b.n - 1; k++)
-    {
-        index = get_specific_hap_overlap(&(all_ovlp->x[b_0->b.a[k]>>1]), b_0->b.a[k]>>1, b_0->b.a[k+1]>>1);
-        if(index == -1) continue;
-        collect_reverse_unitig_pair(link, ug, &(all_ovlp->x[b_0->b.a[k]>>1].a.a[index]));
-    }
-}
-
-
 void link_unitigs(asg_t *purge_g, ma_ug_t *ug, hap_overlaps_list* all_ovlp,
 R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, 
 uint64_t* position_index, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, 
-kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit,
-hc_links* link)
+kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit)
 {
     uint32_t v, n_vtx = purge_g->n_seq * 2, beg, end;
     long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen;
@@ -3996,8 +3948,6 @@ hc_links* link)
         {
             continue;
         }
-
-        if(link) collect_reverse_unitigs_purge(&b_0, link, ug, all_ovlp);
         purge_merge(purge_g, ug, all_ovlp, &b_0, ruIndex, reverse_sources, coverage_cut, 
         read_g, position_index, u_buffer, tailIndex, prevIndex,max_hang, min_ovlp, edge, visit);
     }
@@ -4218,10 +4168,11 @@ uint32_t minLen, double purge_threshold)
     return 0;
 }
 
+
 void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, 
 ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, 
 uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, 
-uint32_t just_contain, uint32_t just_coverage, hc_links* link)
+uint32_t just_contain, uint32_t just_coverage)
 {
     asg_t *purge_g = NULL;
     purge_g = asg_init();
@@ -4278,7 +4229,7 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link)
         for (i = 0, offset = 0; i < reads->n; i++)
         {
             rId = reads->a[i]>>33;
-            set_R_to_U(ruIndex, rId, uId, 1, &(read_g->seq[rId].c));
+            set_R_to_U(ruIndex, rId, uId, 1);
 
             position_index[rId] = offset;
             position_index[rId] = position_index[rId] << 32;
@@ -4350,7 +4301,6 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link)
                 purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].c = ALTER_LABLE;
                 purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].del = 1;
                 all_ovlp.x[uId].a.a[i].status = DELETE;
-                if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i]));
             }
 
             if(all_ovlp.x[uId].a.a[i].type == XCY)
@@ -4359,7 +4309,6 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link)
                 purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].c = ALTER_LABLE;
                 purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].del = 1;
                 all_ovlp.x[uId].a.a[i].status = DELETE;
-                if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i]));
             }
             ///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
         }
@@ -4411,7 +4360,7 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link)
 
         link_unitigs(purge_g, ug, &all_ovlp, ruIndex, reverse_sources, coverage_cut, read_g, position_index, 
         &(hap_buf.buf[0].u_buffer), &(hap_buf.buf[0].u_buffer_tailIndex), &(hap_buf.buf[0].u_buffer_prevIndex),
-        max_hang, min_ovlp, edge, hap_buf.buf[0].visit, link);
+        max_hang, min_ovlp, edge, hap_buf.buf[0].visit);
     }
 
     for (v = 0; v < all_ovlp.num; v++)

Purge_Dups.h

@@ -15,7 +15,7 @@
 void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, 
 ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, 
 uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, 
-uint32_t just_contain, uint32_t just_coverage, hc_links* link);
+uint32_t just_contain, uint32_t just_coverage);
 void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge,
 uint32_t is_circle, uint64_t* rLen);
 void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen);

Trio.cpp

@@ -67,7 +67,7 @@ static yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ..
 {
 	va_list ap;
 	FILE *fp;
-	uint32_t t[3], f_tmp = 0;
+	uint32_t t[3];
 	char magic[4];
 	int i, j, absent, min_cnt = 0, mid_cnt = 0, mode_err = 0;
 	uint64_t mask = (1ULL<<YAK_COUNTER_BITS) - 1, n_ins = 0, n_new = 0;
@@ -92,36 +92,33 @@ static yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ..
 		fclose(fp);
 		return 0;
 	}
-	f_tmp += fread(t, 4, 3, fp);
+	fread(t, 4, 3, fp);
 	if (t[2] != YAK_COUNTER_BITS) {
 		fprintf(stderr, "ERROR: saved counter bits: %d; compile-time counter bits: %d\n", t[2], YAK_COUNTER_BITS);
 		fclose(fp);
 		return 0;
 	}
-	///t[0] = k; t[1] = pre, t[2] = YAK_COUNTER_BITS;
+
 	ch = ch0 == 0? yak_ch_init(t[0], t[1]) : ch0;
 	assert((int)t[0] == ch->k && (int)t[1] == ch->pre);
 	for (i = 0; i < 1<<ch->pre; ++i) {
 		yak_ht_t *h = ch->h[i].h;
-		f_tmp += fread(t, 4, 2, fp);
-		///t[0] = kh_capacity(h), t[1] = kh_size(h);
+		fread(t, 4, 2, fp);
 		if (ch0 == 0) yak_ht_resize(h, t[0]);
 		for (j = 0; j < (int)t[1]; ++j) {
 			uint64_t key;
-			f_tmp += fread(&key, 8, 1, fp);
+			fread(&key, 8, 1, fp);
 			if (mode == YAK_LOAD_ALL) {
 				++n_ins;
 				yak_ht_put(h, key, &absent);
 				if (absent) ++n_new;
 			} else if (mode == YAK_LOAD_TRIOBIN1 || mode == YAK_LOAD_TRIOBIN2) {
 				int cnt = key & mask, x, shift = mode == YAK_LOAD_TRIOBIN1? 0 : 2;
-				//1. filter singleton k-mer; 2. label non-repeat and repeat
 				if (cnt >= mid_cnt) x = 2<<shift;
 				else if (cnt >= min_cnt) x = 1<<shift;
 				else x = -1;
 				if (x >= 0) {
 					khint_t k;
-					///no need cnt at all
 					key = (key & ~mask) | x;
 					++n_ins;
 					k = yak_ht_put(h, key, &absent);

hic.h

@@ -1,72 +0,0 @@
-#ifndef __HIC__
-#define __HIC__
-#include <stdint.h>
-#include "Overlaps.h"
-
-#define kdq_clear(q) ((q)->count = (q)->front = 0)
-#define kv_malloc(v, s) ((v).n = 0, (v).m = (s), MALLOC((v).a, (s)))
-#define RC_0 0
-#define RC_1 1
-#define RC_2 2
-
-hc_edge* get_hc_edge(hc_links* link, uint64_t src, uint64_t dest, uint64_t dir);
-void push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t* d);
-void hic_analysis(ma_ug_t *ug, asg_t* read_g, hc_links* link);
-void hic_benchmark(ma_ug_t *ug, asg_t* read_g);
-
-typedef struct {
-    double w;
-    uint32_t id, occ;
-    ///uint32_t *bid, bid_n;
-    ma_utg_t *u;
-    uint64_t l_d, r_d;
-}chain_hic_w_type;
-
-typedef struct {
-    size_t n, m;
-    chain_hic_w_type* a;
-    uint32_t max_bub_id;
-    uint32_t *chain_idx, u_n;
-}chain_hic_warp;
-
-typedef struct {
-    long long g_occ, b_occ;
-    uint64_t id;
-    uint8_t del;
-}chain_w_type;
-
-typedef struct {
-    uint32_t* index, round_id, n_round;
-    ma_ug_t* ug;
-    kvec_t(uint32_t) list;
-    kvec_t(uint32_t) num;
-    kvec_t(uint64_t) pathLen;
-    kvec_t(uint64_t) b_s_idx;
-    uint64_t s_bub, f_bub, b_bub, b_end_bub, tangle_bub, cross_bub, mess_bub;
-    uint32_t check_het;
-    asg_t *b_g;
-    ma_ug_t* b_ug;
-    kvec_t(chain_w_type) chain_weight;
-    chain_hic_warp c_w;
-} bubble_type;
-#define P_het(B) ((B).num.n)
-#define M_het(B) ((B).num.n + 1)
-// #define IF_BUB(ID, B) ((B).index[(ID)] < (B).num.n)
-// #define IF_HET(ID, B) ((B).index[(ID)] == (B).num.n)
-// #define IF_HOM(ID, B) ((B).index[(ID)] > (B).num.n)
-#define IF_BUB(ID, B) ((B).index[(ID)] < (B).f_bub+1)
-#define IF_HET(ID, B) ((B).index[(ID)] == (B).f_bub+1)
-#define IF_HOM(ID, B) ((B).index[(ID)] > (B).f_bub+1)
-#define Get_bub_num(RECORD) ((RECORD).num.n-1)
-void get_bubbles(bubble_type* bub, uint64_t id, uint32_t* beg, uint32_t* sink, uint32_t** a, uint32_t* n, uint64_t* pathBase);
-int load_hc_links(hc_links* link, const char *fn);
-void write_hc_links(hc_links* link, const char *fn);
-void destory_bubbles(bubble_type* bub);
-void identify_bubbles(ma_ug_t* ug, bubble_type* bub, hc_links* link);
-void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub);
-uint32_t connect_bub_occ(bubble_type* bub, uint32_t root_id, uint32_t check_het);
-void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, uint32_t check_het);
-void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint32_t is_end);
-void set_b_utg_weight_flag(bubble_type* bub, buf_t* b, uint32_t v, uint8_t* vis_flag, uint32_t flag, uint32_t* occ);
-
-#endif

hifiasm.1

@@ -111,11 +111,6 @@ 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,5 +1,4 @@
 #include <stdio.h>
-#include <assert.h>
 #include "htab.h"
 
 static void ha_hist_line(int c, int x, int exceed, int64_t cnt)
@@ -12,19 +11,16 @@ 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, int start_cnt, const int64_t *cnt, int *peak_het)
+int ha_analyze_count(int n_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;
 
-	// determine the start point
-	assert(n_cnt > start_cnt);
+	// find the low point from the left
 	*peak_het = -1;
 	start = cnt[1] > 0? 1 : 2;
-
-	// find the low point from the left
-	low_i = start > start_cnt? start : start_cnt;
+	low_i = start;
 	for (i = low_i + 1; i < n_cnt; ++i)
 		if (cnt[i] > cnt[i-1]) break;
 	low_i = i - 1;

htab.cpp

@@ -566,203 +566,203 @@ static void worker_for_mz(void *data, long i, int tid)
 
 static void *worker_count(void *data, int step, void *in) // callback for kt_pipeline()
 {
-	pl_data_t *p = (pl_data_t*)data;
-	if (step == 0) { // step 1: read a block of sequences
-		int ret;
-		st_data_t *s;
-		CALLOC(s, 1);
-		s->p = p;
-		s->n_seq0 = p->n_seq;
-		if (p->rs_in && (p->flag & HAF_RS_READ)) {
-			while (p->n_seq < p->rs_in->total_reads) {
-				if((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU)
-				{
-					++p->n_seq;
-					continue;
-				}
-				int l;
-				recover_UC_Read(&p->ucr, p->rs_in, p->n_seq);
-				l = p->ucr.length;
-				if (s->n_seq == s->m_seq) {
-					s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
-					REALLOC(s->len, s->m_seq);
-					REALLOC(s->seq, s->m_seq);
-				}
-				MALLOC(s->seq[s->n_seq], l);
-				memcpy(s->seq[s->n_seq], p->ucr.seq, l);
-				s->len[s->n_seq++] = l;
-				++p->n_seq;
-				s->sum_len += l;
-				s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
-				if (s->sum_len >= p->opt->chunk_size)
-					break;
-			}
-		} else {
-			while ((ret = kseq_read(p->ks)) >= 0) {
-				int l = (int)(p->ks->seq.l) - (int)(p->opt->adaLen) - (int)(p->opt->adaLen);
-				if(l <= 0) continue;
+    pl_data_t *p = (pl_data_t*)data;
+    if (step == 0) { // step 1: read a block of sequences
+        int ret;
+        st_data_t *s;
+        CALLOC(s, 1);
+        s->p = p;
+        s->n_seq0 = p->n_seq;
+        if (p->rs_in && (p->flag & HAF_RS_READ)) {
+            while (p->n_seq < p->rs_in->total_reads) {
+                if((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU)
+                {
+                    ++p->n_seq;
+                    continue;
+                }
+                int l;
+                recover_UC_Read(&p->ucr, p->rs_in, p->n_seq);
+                l = p->ucr.length;
+                if (s->n_seq == s->m_seq) {
+                    s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
+                    REALLOC(s->len, s->m_seq);
+                    REALLOC(s->seq, s->m_seq);
+                }
+                MALLOC(s->seq[s->n_seq], l);
+                memcpy(s->seq[s->n_seq], p->ucr.seq, l);
+                s->len[s->n_seq++] = l;
+                ++p->n_seq;
+                s->sum_len += l;
+                s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
+                if (s->sum_len >= p->opt->chunk_size)
+                    break;
+            }
+        } else {
+            while ((ret = kseq_read(p->ks)) >= 0) {
+                int l = (int)(p->ks->seq.l) - (int)(p->opt->adaLen) - (int)(p->opt->adaLen);
+                if(l <= 0) continue;
 
-				if (p->n_seq >= 1<<28) {
-					fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);
-					exit(1);
-				}
-				if (p->rs_out) {
-					///for 0-th count, just insert read length to R_INF, instead of read
-					if (p->flag & HAF_RS_WRITE_LEN) {
-						assert(p->n_seq == p->rs_out->total_reads);
-						ha_insert_read_len(p->rs_out, l, p->ks->name.l);
-					} else if (p->flag & HAF_RS_WRITE_SEQ) {
-						int i, n_N;
-						assert(l == (int)p->rs_out->read_length[p->n_seq]);
-						for (i = n_N = 0; i < l; ++i) // count number of ambiguous bases
-							if (seq_nt4_table[(uint8_t)p->ks->seq.s[i+p->opt->adaLen]] >= 4)
-								++n_N;
-						ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s+p->opt->adaLen, l, &p->rs_out->N_site[p->n_seq], n_N);
-						memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);
-					}
-				}
-				///for 0-th count, insert both seq and length to local block
-				if (s->n_seq == s->m_seq) {
-					s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
-					REALLOC(s->len, s->m_seq);
-					REALLOC(s->seq, s->m_seq);
-				}
-				MALLOC(s->seq[s->n_seq], l);
-				memcpy(s->seq[s->n_seq], p->ks->seq.s+p->opt->adaLen, l);
-				s->len[s->n_seq++] = l;
-				++p->n_seq;
-				s->sum_len += l;
-				s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
-				///p->opt->chunk_size is the block max size
-				if (s->sum_len >= p->opt->chunk_size)
-					break;
-			}
-		}
-		if (s->sum_len == 0) free(s);
-		else return s;
-	} else if (step == 1) { // step 2: extract k-mers
-		///s is the block of reads
-		st_data_t *s = (st_data_t*)in;
-		///for 0-th counting, n_pre = 4096
-		int i, n_pre = 1<<p->opt->pre, m;
-		// allocate the k-mer buffer
-		CALLOC(s->buf, n_pre);
-		m = (int)(s->nk * 1.2 / n_pre) + 1;
-		//pre-allocate memory for each of 4096 buffer
-		for (i = 0; i < n_pre; ++i) {
-			s->buf[i].m = m;
-			///for 0-th counting, p->pt = NULL
-			if (p->pt) MALLOC(s->buf[i].b, m);
-			else MALLOC(s->buf[i].a, m);
-		}
-		// fill the buffer
-		///for 0-th counting, p->opt->w == 1
-		if (p->opt->w == 1) { // enumerate all k-mers
-			///scan all reads
-			for (i = 0; i < s->n_seq; ++i) {
-				if (p->opt->is_HPC)
-					count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
-				else
-					count_seq_buf(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
-				if (!p->is_store) free(s->seq[i]);
-			}
-		} else { // minimizers only
-			uint32_t j;
-			// compute minimizers
-			// s->n_seq is how many reads at this buffer
-			// s->mz && s->mz_buf are lists of minimzer vectors
-			CALLOC(s->mz, s->n_seq);
-			CALLOC(s->mz_buf, p->opt->n_thread);
-			///calculate minimzers for each read, each read corresponds to one thread
-			kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq);
-			for (i = 0; i < p->opt->n_thread; ++i)
-				free(s->mz_buf[i].a);
-			free(s->mz_buf);
-			// insert minimizers
-			if (p->pt) {///insert whole minimizer
-				for (i = 0; i < s->n_seq; ++i)
-					for (j = 0; j < s->mz[i].n; ++j)
-						pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]);
-			} else {///just insert the hash key of minimizer
-				for (i = 0; i < s->n_seq; ++i)
-					for (j = 0; j < s->mz[i].n; ++j)
-						ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);
-			}
-			for (i = 0; i < s->n_seq; ++i) {
-				free(s->mz[i].a);
-				if (!p->is_store) free(s->seq[i]);
-			}
-			free(s->mz);
-		}
-		///just clean seq
-		free(s->seq); free(s->len);
-		s->seq = 0, s->len = 0;
-		return s;
-	} else if (step == 2) { // step 3: insert k-mers to hash table
-		st_data_t *s = (st_data_t*)in;
-		int i, n = 1<<p->opt->pre;
-		uint64_t n_ins = 0;
-		///for 0-th counting, p->pt = NULL
-		kt_for(p->opt->n_thread, worker_for_insert, s, n);
-		///n_ins is number of distinct k-mers
-		for (i = 0; i < n; ++i) {
-			n_ins += s->buf[i].n_ins;
-			if (p->pt) free(s->buf[i].b);
-			else free(s->buf[i].a);
-		}
-		if (p->ct) p->ct->tot += n_ins;
-		if (p->pt) p->pt->tot_pos += n_ins;
-		free(s->buf);
-		#if 0
-		fprintf(stderr, "[M::%s::%.3f*%.2f] processed %ld sequences; %ld %s in the hash table\n", __func__,
-				yak_realtime(), yak_cpu_usage(), (long)s->n_seq0 + s->n_seq,
-				(long)(p->pt? p->pt->tot_pos : p->ct->tot), p->pt? "positions" : "distinct k-mers");
-		#endif
-		free(s);
-	}
-	return 0;
+                if (p->n_seq >= 1<<28) {
+                    fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);
+                    exit(1);
+                }
+                if (p->rs_out) {
+                    ///for 0-th count, just insert read length to R_INF, instead of read
+                    if (p->flag & HAF_RS_WRITE_LEN) {
+                        assert(p->n_seq == p->rs_out->total_reads);
+                        ha_insert_read_len(p->rs_out, l, p->ks->name.l);
+                    } else if (p->flag & HAF_RS_WRITE_SEQ) {
+                        int i, n_N;
+                        assert(l == (int)p->rs_out->read_length[p->n_seq]);
+                        for (i = n_N = 0; i < l; ++i) // count number of ambiguous bases
+                            if (seq_nt4_table[(uint8_t)p->ks->seq.s[i+p->opt->adaLen]] >= 4)
+                                ++n_N;
+                        ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s+p->opt->adaLen, l, &p->rs_out->N_site[p->n_seq], n_N);
+                        memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);
+                    }
+                }
+                ///for 0-th count, insert both seq and length to local block
+                if (s->n_seq == s->m_seq) {
+                    s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
+                    REALLOC(s->len, s->m_seq);
+                    REALLOC(s->seq, s->m_seq);
+                }
+                MALLOC(s->seq[s->n_seq], l);
+                memcpy(s->seq[s->n_seq], p->ks->seq.s+p->opt->adaLen, l);
+                s->len[s->n_seq++] = l;
+                ++p->n_seq;
+                s->sum_len += l;
+                s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
+                ///p->opt->chunk_size is the block max size
+                if (s->sum_len >= p->opt->chunk_size)
+                    break;
+            }
+        }
+        if (s->sum_len == 0) free(s);
+        else return s;
+    } else if (step == 1) { // step 2: extract k-mers
+        ///s is the block of reads
+        st_data_t *s = (st_data_t*)in;
+        ///for 0-th counting, n_pre = 4096
+        int i, n_pre = 1<<p->opt->pre, m;
+        // allocate the k-mer buffer
+        CALLOC(s->buf, n_pre);
+        m = (int)(s->nk * 1.2 / n_pre) + 1;
+        //pre-allocate memory for each of 4096 buffer
+        for (i = 0; i < n_pre; ++i) {
+            s->buf[i].m = m;
+            ///for 0-th counting, p->pt = NULL
+            if (p->pt) MALLOC(s->buf[i].b, m);
+            else MALLOC(s->buf[i].a, m);
+        }
+        // fill the buffer
+        ///for 0-th counting, p->opt->w == 1
+        if (p->opt->w == 1) { // enumerate all k-mers
+            ///scan all reads
+            for (i = 0; i < s->n_seq; ++i) {
+                if (p->opt->is_HPC)
+                    count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
+                else
+                    count_seq_buf(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
+                if (!p->is_store) free(s->seq[i]);
+            }
+        } else { // minimizers only
+            uint32_t j;
+            // compute minimizers
+            // s->n_seq is how many reads at this buffer
+            // s->mz && s->mz_buf are lists of minimzer vectors
+            CALLOC(s->mz, s->n_seq);
+            CALLOC(s->mz_buf, p->opt->n_thread);
+            ///calculate minimzers for each read, each read corresponds to one thread
+            kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq);
+            for (i = 0; i < p->opt->n_thread; ++i)
+                free(s->mz_buf[i].a);
+            free(s->mz_buf);
+            // insert minimizers
+            if (p->pt) {///insert whole minimizer
+                for (i = 0; i < s->n_seq; ++i)
+                    for (j = 0; j < s->mz[i].n; ++j)
+                        pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]);
+            } else {///just insert the hash key of minimizer
+                for (i = 0; i < s->n_seq; ++i)
+                    for (j = 0; j < s->mz[i].n; ++j)
+                        ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);
+            }
+            for (i = 0; i < s->n_seq; ++i) {
+                free(s->mz[i].a);
+                if (!p->is_store) free(s->seq[i]);
+            }
+            free(s->mz);
+        }
+        ///just clean seq
+        free(s->seq); free(s->len);
+        s->seq = 0, s->len = 0;
+        return s;
+    } else if (step == 2) { // step 3: insert k-mers to hash table
+        st_data_t *s = (st_data_t*)in;
+        int i, n = 1<<p->opt->pre;
+        uint64_t n_ins = 0;
+        ///for 0-th counting, p->pt = NULL
+        kt_for(p->opt->n_thread, worker_for_insert, s, n);
+        ///n_ins is number of distinct k-mers
+        for (i = 0; i < n; ++i) {
+            n_ins += s->buf[i].n_ins;
+            if (p->pt) free(s->buf[i].b);
+            else free(s->buf[i].a);
+        }
+        if (p->ct) p->ct->tot += n_ins;
+        if (p->pt) p->pt->tot_pos += n_ins;
+        free(s->buf);
+        #if 0
+        fprintf(stderr, "[M::%s::%.3f*%.2f] processed %ld sequences; %ld %s in the hash table\n", __func__,
+                yak_realtime(), yak_cpu_usage(), (long)s->n_seq0 + s->n_seq,
+                (long)(p->pt? p->pt->tot_pos : p->ct->tot), p->pt? "positions" : "distinct k-mers");
+        #endif
+        free(s);
+    }
+    return 0;
 }
 
 void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs)
 {
-	int ret;
+    int ret;
     uint32_t i, m, pass, unpass;
     gzFile fp = 0;
     kseq_t *ks = NULL;
-	UC_Read ucr;
-	init_UC_Read(&ucr);
+    UC_Read ucr;
+    init_UC_Read(&ucr);
 
     for (i = m = pass = unpass = 0; i < (uint32_t)asm_opt->num_reads; ++i)
     {
-		if ((fp = gzopen(asm_opt->read_file_names[i], "r")) == 0) continue;
-		ks = kseq_init(fp);
-		while ((ret = kseq_read(ks)) >= 0) 
-		{
-			int l = ks->seq.l;
-			if((l - asm_opt->adapterLen*2) <= 0) continue;
-			recover_UC_Read(&ucr, rs, m);
-			fprintf(stderr, "l: %d, ucr.length: %lld, asm_opt->adapterLen: %d\n", 
-			l, ucr.length, asm_opt->adapterLen);
-			if(memcmp(ucr.seq, ks->seq.s+asm_opt->adapterLen, ucr.length) == 0)
-			{
-				pass++;
-			}
-			else
-			{
-				unpass++;
-			}
-			m++;
-		}
+        if ((fp = gzopen(asm_opt->read_file_names[i], "r")) == 0) continue;
+        ks = kseq_init(fp);
+        while ((ret = kseq_read(ks)) >= 0) 
+        {
+            int l = ks->seq.l;
+            if((l - asm_opt->adapterLen*2) <= 0) continue;
+            recover_UC_Read(&ucr, rs, m);
+            fprintf(stderr, "l: %d, ucr.length: %lld, asm_opt->adapterLen: %d\n", 
+            l, ucr.length, asm_opt->adapterLen);
+            if(memcmp(ucr.seq, ks->seq.s+asm_opt->adapterLen, ucr.length) == 0)
+            {
+                pass++;
+            }
+            else
+            {
+                unpass++;
+            }
+            m++;
+        }
         kseq_destroy(ks);
         gzclose(fp);
-		ks = NULL;
-		fp = 0;
+        ks = NULL;
+        fp = 0;
     }
 
-	destory_UC_Read(&ucr);
+    destory_UC_Read(&ucr);
 
-	fprintf(stderr, "[M::%s::# reads: %u, # pass: %u, # unpass: %u\n]", __func__, m, pass, unpass);
-	exit(1);
+    fprintf(stderr, "[M::%s::# reads: %u, # pass: %u, # unpass: %u\n]", __func__, m, pass, unpass);
+    exit(1);
 }
 
 static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt_t *p0, ha_ct_t *c0, const void *flt_tab, All_reads *rs, int64_t *n_seq)
@@ -814,35 +814,35 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 
 ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const void *flt_tab, All_reads *rs)
 {
-	int i;
-	int64_t n_seq = 0;
-	yak_copt_t opt;
-	ha_ct_t *h = 0;
-	assert(!(flag & HAF_RS_WRITE_LEN) || !(flag & HAF_RS_WRITE_SEQ)); // not both
-	///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN
-	if (rs) {
-		if (flag & HAF_RS_WRITE_LEN)
-			init_All_reads(rs);
-		else if (flag & HAF_RS_WRITE_SEQ)
-			malloc_All_reads(rs);
-	}
-	yak_copt_init(&opt);
-	opt.k = asm_opt->k_mer_length;
-	///always 0
-	opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC);
-	///for ft-counting, shoud be 1
-	opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win;
-	///for ft-counting, shoud be 37
-	///for ha_pt_gen, shoud be 0
-	opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift;
-	opt.n_thread = asm_opt->thread_num;
-	opt.adaLen = asm_opt->adapterLen;
-	///asm_opt->num_reads is the number of fastq files
-	for (i = 0; i < asm_opt->num_reads; ++i)
-		h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, &n_seq);
-	if (h && opt.bf_shift > 0)
-		ha_ct_destroy_bf(h);
-	return h;
+    int i;
+    int64_t n_seq = 0;
+    yak_copt_t opt;
+    ha_ct_t *h = 0;
+    assert(!(flag & HAF_RS_WRITE_LEN) || !(flag & HAF_RS_WRITE_SEQ)); // not both
+    ///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN
+    if (rs) {
+        if (flag & HAF_RS_WRITE_LEN)
+            init_All_reads(rs);
+        else if (flag & HAF_RS_WRITE_SEQ)
+            malloc_All_reads(rs);
+    }
+    yak_copt_init(&opt);
+    opt.k = asm_opt->k_mer_length;
+    ///always 0
+    opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC);
+    ///for ft-counting, shoud be 1
+    opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win;
+    ///for ft-counting, shoud be 37
+    ///for ha_pt_gen, shoud be 0
+    opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift;
+    opt.n_thread = asm_opt->thread_num;
+    opt.adaLen = asm_opt->adapterLen;
+    ///asm_opt->num_reads is the number of fastq files
+    for (i = 0; i < asm_opt->num_reads; ++i)
+        h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, &n_seq);
+    if (h && opt.bf_shift > 0)
+        ha_ct_destroy_bf(h);
+    return h;
 }
 
 /***************************
@@ -934,7 +934,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, asm_opt->min_hist_kmer_cnt, cnt, &peak_het);
+		peak_hom = ha_analyze_count(YAK_N_COUNTS, 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
@@ -971,7 +971,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, asm_opt->min_hist_kmer_cnt, cnt, &peak_het);
+	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);
@@ -1048,17 +1048,16 @@ int load_ct_index(void **i_ct_idx, char* file_name)
     }
 	ha_ct_t** ct_idx = (ha_ct_t**)i_ct_idx;
 	double index_time = 0;
-	uint64_t flag = 0;
 	int i;
 	ha_ct_t *h = 0;
 	ha_ct1_t *g;
 	CALLOC(h, 1);
 
-	flag += fread(&h->k, sizeof(h->k), 1, fp);
-	flag += fread(&h->pre, sizeof(h->pre), 1, fp);
-	flag += fread(&h->n_hash, sizeof(h->n_hash), 1, fp);
-	flag += fread(&h->n_shift, sizeof(h->n_shift), 1, fp);
-	flag += fread(&h->tot, sizeof(h->tot), 1, fp);
+	fread(&h->k, sizeof(h->k), 1, fp);
+	fread(&h->pre, sizeof(h->pre), 1, fp);
+	fread(&h->n_hash, sizeof(h->n_hash), 1, fp);
+	fread(&h->n_shift, sizeof(h->n_shift), 1, fp);
+	fread(&h->tot, sizeof(h->tot), 1, fp);
 	CALLOC(h->h, 1<<h->pre);
 
 
@@ -1231,4 +1230,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, int start_cnt, const int64_t *cnt, int *peak_het);
+int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het);
 void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs);
 
 static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function

khashl.h

@@ -147,14 +147,13 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 	SCOPE khint_t prefix##_load(HType **h, FILE* fp) { \
 		(*h) = prefix##_init(); \
 		khint_t n_buckets; \
-		uint64_t flag = 0;\
-		flag += fread(&n_buckets, sizeof(n_buckets), 1, fp); \
-		flag += fread(&(*h)->bits, sizeof((*h)->bits), 1, fp); \
-		flag += fread(&(*h)->count, sizeof((*h)->count), 1, fp); \
+		fread(&n_buckets, sizeof(n_buckets), 1, fp); \
+		fread(&(*h)->bits, sizeof((*h)->bits), 1, fp); \
+		fread(&(*h)->count, sizeof((*h)->count), 1, fp); \
 		(*h)->used = (khint32_t*)kmalloc(__kh_fsize(n_buckets) * sizeof(khint32_t)); \
 		(*h)->keys = (khkey_t*)kmalloc(n_buckets * sizeof(khkey_t)); \
-		flag += fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
-		flag += fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); \
+		fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
+		fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); \
 		return 1; \
 	} \