Merge branch 'master' into hifiasm-v0.14

merge README

Assembly.cpp

@@ -12,7 +12,7 @@
 #include "kthread.h"
 
 void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, 
-int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar);
+int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct);
 void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list);
 
 All_reads R_INF;
@@ -299,30 +299,6 @@ void push_overlaps(ma_hit_t_alloc* paf, overlap_region_alloc* overlap_list, int
     }
 }
 
-int if_exact_match(char* x, long long xLen, char* y, long long yLen, long long xBeg, long long xEnd, long long yBeg, long long yEnd)
-{
-    long long overlapLen = xEnd - xBeg + 1;
-
-    if(yEnd - yBeg + 1 == overlapLen)
-    {
-        long long i;
-
-        for (i = 0; i < overlapLen; i++)
-        {
-            if(x[xBeg + i] != y[yBeg + i])
-            {
-                break;
-            }
-        }
-
-        if(i == overlapLen)
-        {
-            return 1;
-        }
-    }
-
-    return 0;
-}
 
 long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, overlap_region_alloc* overlap_list, int flag)
 {
@@ -561,7 +537,7 @@ static void worker_ovec(void *data, long i, int tid)
 	int fully_cov, abnormal;
 
     ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 
-    0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region));
+    0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL);
 
 	clear_Cigar_record(&b->cigar1);
 	clear_Round2_alignment(&b->round2);
@@ -606,7 +582,6 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
     uint64_t k, queryNameLen;
     for (k = 0; k < R_INF_FLAG.query_num; k++)
     {
-        
         queryNameLen = strlen(R_INF_FLAG.read_name[k]);
         if (queryNameLen != Get_NAME_LENGTH((R_INF),i)) continue;
         if (memcmp(R_INF_FLAG.read_name[k], Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
@@ -617,10 +592,10 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
 
     if(k < R_INF_FLAG.query_num)
     {
-        int fully_cov, abnormal;
+        int fully_cov, abnormal, q_idx = k;
 
         ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 
-        0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region));
+        0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), &(R_INF_FLAG.candidate_count[q_idx]));
 
         clear_Cigar_record(&b->cigar1);
         clear_Round2_alignment(&b->round2);
@@ -649,6 +624,18 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
         fprintf(R_INF_FLAG.fp, "\n>%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i));
         fprintf(R_INF_FLAG.fp, "%d-th round, len: %lu, hom_cov: %d, max_n_chain: %d\n", 
             asm_opt.number_of_round, Get_READ_LENGTH(R_INF, i), asm_opt.hom_cov, asm_opt.max_n_chain);
+
+        fprintf(R_INF_FLAG.fp, "***************************k-mer counts (%d)***************************\n", (int)(R_INF_FLAG.candidate_count[q_idx].a.n));
+        
+        sort_kvec_t_u64_warp(&(R_INF_FLAG.candidate_count[q_idx]), 0);
+        for (k = 0; k < R_INF_FLAG.candidate_count[q_idx].a.n; k++) 
+        {
+            fprintf(R_INF_FLAG.fp, "[%lu] Count(%u): %lu, filtered: %lu\n", k, 
+            (uint32_t)R_INF_FLAG.candidate_count[q_idx].a.a[k], R_INF_FLAG.candidate_count[q_idx].a.a[k]>>33, 
+            (R_INF_FLAG.candidate_count[q_idx].a.a[k]>>32)&(uint64_t)1);
+        }
+
+        
         fprintf(R_INF_FLAG.fp, "***************************forward ovlp***************************\n");
         for (k = 0; k < b->olist.length; k++) 
         {
@@ -692,6 +679,7 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
         // }
         
         fprintf(R_INF_FLAG.fp, "R_INF.trio_flag[%ld]: %u\n", i, R_INF.trio_flag[i]);
+    
 
 	    pthread_mutex_unlock(&R_INF_FLAG.OutputMutex);
     }
@@ -860,14 +848,15 @@ void ha_overlap_and_correct(int round)
 		b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1));
     if(ha_idx) hom_cov = asm_opt.hom_cov;
 	if(ha_idx == NULL) ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
-	if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
+	///debug_adapter(&asm_opt, &R_INF);
+    if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
 		ha_opt_update_cov(&asm_opt, hom_cov);
 	if (asm_opt.required_read_name)
 		kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads);
 	else
 		kt_for(asm_opt.thread_num, worker_ovec, b, R_INF.total_reads);
 
-    if (r_out) write_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
+    if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
 	ha_pt_destroy(ha_idx);
 	ha_idx = NULL;
 
@@ -1298,7 +1287,7 @@ static void worker_ov_final(void *data, long i, int tid)
 
 	//get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
     ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 0.001, 
-    asm_opt.max_n_chain, 0, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region));
+    asm_opt.max_n_chain, 0, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL);
 
 	overlap_region_sort_y_id(b->olist.list, b->olist.length);
 	ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length);
@@ -1364,7 +1353,7 @@ static void worker_ov_final_high_het(void *data, long i, int tid)
     ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
 
     ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, HIGH_HET_ERROR_RATE, 
-    asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region));
+    asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL);
 
     overlap_region_sort_y_id(b->olist.list, b->olist.length);
     ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length);
@@ -1612,6 +1601,8 @@ 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);
@@ -1631,7 +1622,7 @@ int ha_assemble(void)
 	}
 	if (!ovlp_loaded) {
         ha_flt_tab = ha_idx = NULL;
-        if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
+        if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name), load_ct_index(&ha_ct_table, asm_opt.output_file_name);
 
 		// construct hash table for high occurrence k-mers
 		if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL) 

CommandLines.cpp

@@ -24,6 +24,12 @@ 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 },
 	{ 0, 0, 0 }
 };
 
@@ -59,6 +65,11 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "    -x FLOAT    max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate);
     fprintf(stderr, "    -y FLOAT    min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate);
     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");
 
@@ -80,6 +91,9 @@ 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");
@@ -94,8 +108,12 @@ 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;
@@ -104,8 +122,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;
@@ -126,18 +143,40 @@ 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->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]);
 }
 
 void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round)
@@ -177,6 +216,16 @@ 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)
@@ -319,12 +368,57 @@ int check_option(hifiasm_opt_t* asm_opt)
         return 0;
     }
 
+    if(asm_opt->bed_inconsist_rate < 0 || asm_opt->bed_inconsist_rate > 100)
+    {
+        fprintf(stderr, "[ERROR] inconsistency rate should be [0, 100] (--lowQ)\n");
+        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;
     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);
@@ -374,6 +468,61 @@ 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)
 {
@@ -438,6 +587,12 @@ 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);
@@ -456,15 +611,16 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 		}
     }
 
-
+    
     if (argc == opt.ind)
     {
         Print_H(asm_opt);
         return 0;
     }
-    ///fprintf(stderr, "max_ov_diff_ec: %f, max_ov_diff_final: %f\n", asm_opt->max_ov_diff_ec, asm_opt->max_ov_diff_final);
 
     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.12-r304"
+#define HA_VERSION "0.14-r309"
 
 #define VERBOSE 0
 
@@ -21,6 +21,11 @@
 
 #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;
@@ -30,9 +35,13 @@ 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
@@ -40,9 +49,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 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;
@@ -59,14 +68,18 @@ 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;
@@ -92,4 +105,9 @@ 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 
@@ -1175,6 +1175,20 @@ long long get_affine_gap_score(overlap_region* ovc, UC_Read* g_read, UC_Read* ov
 uint8_t* y_num, uint64_t EstimateXOlen, uint64_t EstimateYOlen);
 int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u8_warp* k_flag, float hp_rate, int rlen, FILE* fp);
 
+inline int if_exact_match(char* x, long long xLen, char* y, long long yLen, long long xBeg, long long xEnd, long long yBeg, long long yEnd)
+{
+    long long overlapLen = xEnd - xBeg + 1;
+
+    if(yEnd - yBeg + 1 == overlapLen)
+    {
+        if(memcmp(x+xBeg, y+yBeg, overlapLen)==0) return 1;
+    }
+
+    return 0;
+}
+
+
+
 #define FORWARD_KSW 0
 #define BACKWARD_KSW 1
 #define MATCH_SCORE_KSW 2

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
+			htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm
 
@@ -71,3 +71,4 @@ 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,6 +21,7 @@
 #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
@@ -106,6 +107,11 @@ 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
@@ -155,6 +161,7 @@ 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 {
@@ -394,7 +401,8 @@ 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;
@@ -442,11 +450,8 @@ 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, 
@@ -455,18 +460,6 @@ 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)))
 
@@ -481,14 +474,10 @@ 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);
+void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig, uint8_t* flag);
 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);
 
@@ -743,13 +732,12 @@ 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;
@@ -770,7 +758,6 @@ 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****************************************/
@@ -781,12 +768,11 @@ 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);
+				set_R_to_U(ruIndex, qn, (b_max.b_0->b.a[b_max.untigI]>>1), 1, &(read_sg->seq[qn].c));
 			}
 		}
 		/*****************************label all unitigs****************************************/
 
-
 		///each unitig
 		for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++)
 		{
@@ -820,7 +806,6 @@ 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++)
 		{
@@ -833,7 +818,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 			}
 		}
 		/*****************************label all unitigs****************************************/
-		
 	}
 	else
 	{
@@ -841,7 +825,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);
+			set_R_to_U(ruIndex, qn, 1, 1, &(read_sg->seq[qn].c));
 		}
 		/*****************************label all reads****************************************/
 
@@ -896,8 +880,6 @@ 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)
@@ -1047,8 +1029,6 @@ 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);
@@ -1061,9 +1041,77 @@ 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);
+uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes);
 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
 

Process_Read.cpp

@@ -338,93 +338,6 @@ void init_UC_Read(UC_Read* r)
 	}
 }
 
-void recover_UC_Read_sub_region_begin_end(char* r, long long start_pos, long long length, uint8_t strand,
-										  All_reads* R_INF, long long ID, int extra_begin, int extra_end)
-{
-	long long readLen = Get_READ_LENGTH((*R_INF), ID);
-	uint8_t* src = Get_READ((*R_INF), ID);
-
-	long long i;
-	long long copyLen;
-	long long end_pos = start_pos + length - 1;
-
-	if (strand == 0)
-	{
-		i = start_pos;
-		copyLen = 0;
-
-		long long initLen = start_pos % 4;
-
-		if (initLen != 0)
-		{
-			memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen);
-			copyLen = copyLen + 4 - initLen;
-			i = i + copyLen;
-		}
-		while (copyLen < length)
-		{
-			memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4);
-			copyLen = copyLen + 4;
-			i = i + 4;
-		}
-
-		if (R_INF->N_site[ID])
-		{
-			for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
-			{
-				if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)R_INF->N_site[ID][i] <= end_pos)
-				{
-					r[R_INF->N_site[ID][i] - start_pos] = 'N';
-				}
-				else if((long long)R_INF->N_site[ID][i] > end_pos)
-				{
-					break;
-				}
-			}
-		}
-	}
-	else
-	{
-		start_pos = readLen - start_pos - 1;
-		end_pos = readLen - end_pos - 1;
-
-		///start_pos > end_pos
-		i = start_pos;
-		copyLen = 0;
-		long long initLen = (start_pos + 1) % 4;
-
-		if (initLen != 0)
-		{
-			memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen);
-			copyLen = copyLen + initLen;
-			i = i - initLen;
-		}
-
-		while (copyLen < length)
-		{
-			memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4);
-			copyLen = copyLen + 4;
-			i = i - 4;
-		}
-
-		if (R_INF->N_site[ID])
-		{
-			long long offset = readLen - start_pos - 1;
-
-			for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
-			{
-				if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos)
-				{
-					r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N';
-				}
-				else if((long long)R_INF->N_site[ID][i] > start_pos)
-				{
-					break;
-				}
-			}
-		}
-	}
-}
 
 void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID)
 {
@@ -515,6 +428,104 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length,
 }
 
 
+void recover_UC_sub_Read(UC_Read* i_r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID)
+{
+	i_r->length = length;i_r->RID = ID;
+	if (i_r->length + 8 > i_r->size)
+	{
+		i_r->size = i_r->length + 4;
+		i_r->seq = (char*)realloc(i_r->seq,sizeof(char)*(i_r->size));
+	}
+	char* r = i_r->seq;
+	long long readLen = Get_READ_LENGTH((*R_INF), ID);
+	uint8_t* src = Get_READ((*R_INF), ID);
+
+	long long i;
+	long long copyLen;
+	long long end_pos = start_pos + length - 1;
+
+	if (strand == 0)
+	{
+		i = start_pos;
+		copyLen = 0;
+
+		long long initLen = start_pos % 4;
+
+		if (initLen != 0)
+		{
+			memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen);
+			copyLen = copyLen + 4 - initLen;
+			i = i + copyLen;
+		}
+
+		while (copyLen < length)
+		{
+			memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4);
+			copyLen = copyLen + 4;
+			i = i + 4;
+		}
+
+		if (R_INF->N_site[ID])
+		{
+			for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
+			{
+				if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)R_INF->N_site[ID][i] <= end_pos)
+				{
+					r[R_INF->N_site[ID][i] - start_pos] = 'N';
+				}
+				else if((long long)R_INF->N_site[ID][i] > end_pos)
+				{
+					break;
+				}
+			}
+		}
+	}
+	else
+	{
+		start_pos = readLen - start_pos - 1;
+		end_pos = readLen - end_pos - 1;
+
+		///start_pos > end_pos
+		i = start_pos;
+		copyLen = 0;
+		long long initLen = (start_pos + 1) % 4;
+
+		if (initLen != 0)
+		{
+			memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen);
+			copyLen = copyLen + initLen;
+			i = i - initLen;
+		}
+
+		while (copyLen < length)
+		{
+			memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4);
+			copyLen = copyLen + 4;
+			i = i - 4;
+		}
+
+		if (R_INF->N_site[ID])
+		{
+			long long offset = readLen - start_pos - 1;
+
+			for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
+			{
+				if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos)
+				{
+					r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N';
+				}
+				else if((long long)R_INF->N_site[ID][i] > start_pos)
+				{
+					break;
+				}
+			}
+		}
+	}
+}
+
+
+
+
 void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID)
 {
 	r->length = Get_READ_LENGTH((*R_INF), ID);
@@ -690,6 +701,7 @@ void init_Debug_reads(Debug_reads* x, const char* file)
 		x->query_num++;
 	}
 	x->read_name = (char**)malloc(sizeof(char*)*x->query_num);
+	x->candidate_count = (kvec_t_u64_warp*)malloc(sizeof(kvec_t_u64_warp)*x->query_num);
 	fseek(x->fp, 0, SEEK_SET);
 
 	i = 0;
@@ -699,6 +711,7 @@ void init_Debug_reads(Debug_reads* x, const char* file)
 		x->read_name[i] = (char*)malloc(sizeof(char)*(nameLen+1));
 		memcpy(x->read_name[i], Name_Buffer, sizeof(char)*nameLen);
 		x->read_name[i][nameLen] = '\0';
+		kv_init(x->candidate_count[i].a);
 		i++;
 	}
 	
@@ -716,6 +729,7 @@ void destory_Debug_reads(Debug_reads* x)
 	for (i = 0; i < x->query_num; i++)
 	{
 		free(x->read_name[i]);
+		kv_destroy(x->candidate_count[i].a);
 	}
 	
 	free(x->read_name);

Process_Read.h

@@ -148,6 +148,7 @@ typedef struct
 {
     char** read_name;
     uint64_t query_num;
+    kvec_t_u64_warp* candidate_count;
     FILE* fp;
     pthread_mutex_t OutputMutex;
 } Debug_reads;
@@ -168,4 +169,5 @@ void destory_All_reads(All_reads* r);
 int destory_read_bin(All_reads* r);
 void init_Debug_reads(Debug_reads* x, const char* file);
 void destory_Debug_reads(Debug_reads* x);
+void recover_UC_sub_Read(UC_Read* i_r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID);
 #endif

Purge_Dups.cpp

@@ -7,6 +7,7 @@
 #include "Correct.h"
 #include "kthread.h"
 #include "kdq.h"
+#include "hic.h"
 
 KDQ_INIT(uint64_t)
 
@@ -3291,7 +3292,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);
+			n_pop += asg_bub_pop1_primary_trio(purge_g, NULL, v, max_dist, &b, (uint32_t)-1, DROP, 1, NULL, NULL);
 	}
 	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);
@@ -3925,10 +3926,57 @@ 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)
+kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit,
+hc_links* link)
 {
     uint32_t v, n_vtx = purge_g->n_seq * 2, beg, end;
     long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen;
@@ -3948,6 +3996,8 @@ kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edg
         {
             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);
     }
@@ -4168,11 +4218,10 @@ 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)
+uint32_t just_contain, uint32_t just_coverage, hc_links* link)
 {
     asg_t *purge_g = NULL;
     purge_g = asg_init();
@@ -4229,7 +4278,7 @@ uint32_t just_contain, uint32_t just_coverage)
         for (i = 0, offset = 0; i < reads->n; i++)
         {
             rId = reads->a[i]>>33;
-            set_R_to_U(ruIndex, rId, uId, 1);
+            set_R_to_U(ruIndex, rId, uId, 1, &(read_g->seq[rId].c));
 
             position_index[rId] = offset;
             position_index[rId] = position_index[rId] << 32;
@@ -4301,6 +4350,7 @@ uint32_t just_contain, uint32_t just_coverage)
                 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)
@@ -4309,6 +4359,7 @@ uint32_t just_contain, uint32_t just_coverage)
                 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]));
         }
@@ -4360,7 +4411,7 @@ uint32_t just_contain, uint32_t just_coverage)
 
         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);
+        max_hang, min_ovlp, edge, hap_buf.buf[0].visit, link);
     }
 
     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);
+uint32_t just_contain, uint32_t just_coverage, hc_links* link);
 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);

README.md

@@ -8,7 +8,7 @@ cd hifiasm && make
 # Run on test data (use -f0 for small datasets)
 wget https://github.com/chhylp123/hifiasm/releases/download/v0.7/chr11-2M.fa.gz
 ./hifiasm -o test -t4 -f0 chr11-2M.fa.gz 2> test.log
-awk '/^S/{print ">"$1;print $2}' test.p_ctg.gfa > test.p_ctg.fa  # get primary contigs in FASTA
+awk '/^S/{print ">"$2;print $3}' test.p_ctg.gfa > test.p_ctg.fa  # get primary contigs in FASTA
 
 # Assemble inbred/homozygous genomes (-l0 disables duplication purging)
 hifiasm -o CHM13.asm -t32 -l0 CHM13-HiFi.fa.gz 2> CHM13.asm.log
@@ -30,6 +30,29 @@ produce primary/alternate assemblies of quality competitive with the best
 assemblers. It also introduces a new graph binning algorithm and achieves
 the best haplotype-resolved assembly given trio data.
 
+## Why Hifiasm?
+
+* Hifiasm delivers high-quality assemblies. It tends to generate longer contigs
+  and resolve more segmental duplications than other assemblers.
+
+* Given sequence reads from the parents, hifiasm can produce overall the best
+  haplotype-resolved assembly so far. It is the assembler of choice by the
+  [Human Pangenome Project][hpp] for the first batch of samples.
+
+* Hifiasm can purge duplications between haplotigs without relying on
+  third-party tools such as purge\_dups. Hifiasm does not need polishing tools
+  like pilon or racon, either. This simplifies the assembly pipeline and saves
+  running time.
+
+* Hifiasm is fast. It can assemble a human genome in half a day and assemble a
+  ~30Gb redwood genome in three days. No genome is too large for hifiasm.
+
+* Hifiasm is trivial to install and easy to use. It does not required python,
+  R or C++11 compilers and can be compiled into a single executable. The
+  default setting works well with a variety of genomes.
+
+[hpp]: https://humanpangenome.org
+
 ## Usage
 
 A typical hifiasm command line looks like:
@@ -144,12 +167,12 @@ redwood genome in a few days on a single machine. For trio binning assembly:
 [NA12891-data]: https://www.ebi.ac.uk/ena/data/view/ERR194160
 [NA12892-data]: https://www.ebi.ac.uk/ena/data/view/ERR194161
 
-Except NA12878, the assemblies above were produced by hifiasm v0.7 and can be
+Except NA12878, the assemblies above were produced by hifiasm v0.12 and can be
 downloaded at
 ```txt
-ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm/submission/v0.7/
+ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm/submission/hifiasm-0.12/
 ```
-NA12878 was assembled with a more recent version of hifiasm and is available at
+NA12878 was assembled with an older version of hifiasm and is available at
 ```txt
 ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm/NA12878-r253/
 ```
@@ -169,4 +192,4 @@ page](https://github.com/chhylp123/hifiasm/issues).
 
 ## Limitations
 
-1. Purging haplotig duplications may introduce misassemblies.
+1. Purging haplotig duplications may introduce misassemblies.

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];
+	uint32_t t[3], f_tmp = 0;
 	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,33 +92,36 @@ static yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ..
 		fclose(fp);
 		return 0;
 	}
-	fread(t, 4, 3, fp);
+	f_tmp += 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;
-		fread(t, 4, 2, fp);
+		f_tmp += fread(t, 4, 2, fp);
+		///t[0] = kh_capacity(h), t[1] = kh_size(h);
 		if (ch0 == 0) yak_ht_resize(h, t[0]);
 		for (j = 0; j < (int)t[1]; ++j) {
 			uint64_t key;
-			fread(&key, 8, 1, fp);
+			f_tmp += 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);

anchor.cpp

@@ -58,7 +58,7 @@ int ha_ov_type(const overlap_region *r, uint32_t len)
 }
 
 void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag,
-kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar)
+kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct)
 {
 	uint32_t i, rlen;
 	uint64_t k, l;
@@ -73,7 +73,7 @@ kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f
 	rlen = Get_READ_LENGTH(R_INF, rid); // read length
 
 	// get the list of anchors
-	ha_sketch_query(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag);
+	ha_sketch_query(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag, dbg_ct);
 	// minimizer of queried read
 	if (ab->mz.m > ab->old_mz_m) {
 		ab->old_mz_m = ab->mz.m;
@@ -205,14 +205,14 @@ void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf)
 
 
 void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, 
-int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar)
+int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct)
 {
 	extern void *ha_flt_tab;
 	extern ha_pt_t *ha_idx;
 	extern void *ha_flt_tab_hp;
 	extern ha_pt_t *ha_idx_hp;
 
-	ha_get_new_candidates(ab, rid, ucr, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar);
+	ha_get_new_candidates(ab, rid, ucr, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar, dbg_ct);
 
 	if(ha_idx_hp)
 	{
@@ -242,7 +242,7 @@ int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp*
 		overlap_list->length = k;
 
 
-		ha_get_new_candidates(ab, rid, ucr, overlap_list_hp, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab_hp, ha_idx_hp, f_cigar);	
+		ha_get_new_candidates(ab, rid, ucr, overlap_list_hp, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab_hp, ha_idx_hp, f_cigar, dbg_ct);	
 		
 		if(overlap_list->length + overlap_list_hp->length > overlap_list->size)
 		{

hic.h

@@ -0,0 +1,72 @@
+#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,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
@@ -198,7 +203,14 @@ Min and max coverage cutoff of primary contigs.
 Keep contigs with coverage in this range at p_ctg.gfa. 
 Inferred automatically in default.
 If INT2 is not specified, it is set to infinity. 
-Set -1 to disable 
+Set -1 to disable.
+
+.TP
+.BI --lowQ \ INT
+Output contig regions with >=INT% inconsistency to the bed file 
+with suffix
+.B lowQ.bed
+[70]. Set 0 to disable. 
 
 .SS Trio-partition options
 

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

@@ -37,6 +37,7 @@ void *ha_flt_tab;
 ha_pt_t *ha_idx;
 void *ha_flt_tab_hp;
 ha_pt_t *ha_idx_hp;
+void *ha_ct_table;
 
 /***************************
  * Yak specific parameters *
@@ -48,6 +49,7 @@ typedef struct {
 	int32_t pre;
 	int32_t n_thread;
 	int64_t chunk_size;
+	int adaLen;
 } yak_copt_t;
 
 void yak_copt_init(yak_copt_t *o)
@@ -132,8 +134,8 @@ typedef struct {
 	ha_ct1_t *h;
 } ha_ct_t;
 
-///for 0-th counting, k = 51, pre = 12, n_hash = 4, n_shift = 0
-///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 37
+///for 0-th counting, k = 51, pre = 12, n_hash = 4, n_shift = 37
+///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0
 static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift)
 {
 	ha_ct_t *h;
@@ -146,8 +148,7 @@ static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift)
 	///it seems there is a large hash table h, consisting 4096 small hash tables
 	for (i = 0; i < 1<<h->pre; ++i)
 		h->h[i].h = yak_ct_init();
-	///for 0-th counting, don't enter here
-	///seems used for minimzer
+	///for 0-th counting, enter here; used for bloom filter
 	if (n_hash > 0 && n_shift > h->pre) {
 		h->n_hash = n_hash, h->n_shift = n_shift;
 		for (i = 0; i < 1<<h->pre; ++i)
@@ -598,7 +599,9 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 			}
 		} else {
 			while ((ret = kseq_read(p->ks)) >= 0) {
-				int l = p->ks->seq.l;
+				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);
@@ -612,9 +615,9 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 						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]] >= 4)
+							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, l, &p->rs_out->N_site[p->n_seq], 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);
 					}
 				}
@@ -625,7 +628,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 					REALLOC(s->seq, s->m_seq);
 				}
 				MALLOC(s->seq[s->n_seq], l);
-				memcpy(s->seq[s->n_seq], p->ks->seq.s, 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;
@@ -720,6 +723,48 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 	return 0;
 }
 
+void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs)
+{
+	int ret;
+    uint32_t i, m, pass, unpass;
+    gzFile fp = 0;
+    kseq_t *ks = NULL;
+	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++;
+		}
+        kseq_destroy(ks);
+        gzclose(fp);
+		ks = NULL;
+		fp = 0;
+    }
+
+	destory_UC_Read(&ucr);
+
+	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)
 {
 	///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN|HAF_CREATE_NEW
@@ -749,10 +794,10 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 	} else if (c0) {
 		pl.ct = c0, pl.create_new = !!(flag&HAF_CREATE_NEW);
 		assert(c0->k == opt->k && c0->pre == opt->pre);
-	} else {///for 0-th counting and 1-th counting, go into here
+	} else {///for ft-th counting and 1-th counting, go into here
 		pl.create_new = 1; // alware create new elements if the count table is empty
-		///for 0-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0
-		///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 37
+		///for 0-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 37
+		///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0
 		///building a large hash table consisting of 4096 small hash tables
 		pl.ct = ha_ct_init(opt->k, opt->pre, opt->bf_n_hash, opt->bf_shift);
 	}
@@ -785,13 +830,13 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
 	opt.k = asm_opt->k_mer_length;
 	///always 0
 	opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC);
-	///for 0-th counting, shoud be 1
-	///for 1-th counting, shoud be 51
+	///for ft-counting, shoud be 1
 	opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win;
-	///for 0-th counting, shoud be 0
-	///for 1-th counting, shoud be 37
+	///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);
@@ -839,6 +884,32 @@ void ha_ft_destroy(void *h)
 	if (h) yak_ft_destroy((yak_ft_t*)h);
 }
 
+
+void debug_ct_index(void* q_ct_idx, void* r_ct_idx)
+{
+	ha_ct_t* ct_idx = (ha_ct_t*)q_ct_idx;
+	yak_ct_t *g = NULL;
+	uint64_t i;
+	khint_t k;
+	for (i = 0; (int)i < 1<<ct_idx->pre; i++)
+	{
+		g = ct_idx->h[i].h;
+		for (k = 0; k < kh_end(g); ++k) 
+		{
+			if (kh_exist(g, k)) 
+			{
+				int c = kh_key(g, k) & YAK_MAX_COUNT;
+				uint64_t hash = ((kh_key(g, k) >> ct_idx->pre)<<ct_idx->pre) | i;
+				int q = query_ct_index(r_ct_idx, hash);
+				if(q!=c)
+				{
+					fprintf(stderr, "ERROR:c: %d, q: %d\n", c, q);
+				}
+			}
+		}
+	}
+}
+
 /*************************
  * High-level interfaces *
  *************************/
@@ -851,10 +922,19 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i
 	if(is_hp_mode) ex_flag = HAF_RS_READ|HAF_SKIP_READ;
 	ha_ct_t *h;
 	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs);
+	if((asm_opt->flag & HA_F_VERBOSE_GFA))
+	{
+		write_ct_index((void*)h, asm_opt->output_file_name);
+		// load_ct_index(&ha_ct_table, asm_opt->output_file_name);
+		// debug_ct_index((void*)h, ha_ct_table);
+		// debug_ct_index(ha_ct_table, (void*)h);
+		// ha_ct_destroy((ha_ct_t *)ha_ct_table);
+	}
+	
 	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
@@ -891,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, 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);
@@ -917,9 +997,86 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 	return pt;
 }
 
+int query_ct_index(void* ct_idx, uint64_t hash)
+{
+	ha_ct1_t *g = &(((ha_ct_t*)ct_idx)->h[hash & ((1ULL<<((ha_ct_t*)ct_idx)->pre) - 1)]);
+	khint_t k;
+	k = yak_ct_get(g->h, hash);
+	if (k == kh_end(g->h)) return 0;
+	return kh_key(g->h, k)&YAK_MAX_COUNT;
+}
 
 
-int write_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
+int write_ct_index(void *i_ct_idx, char* file_name)
+{
+	char* gfa_name = (char*)malloc(strlen(file_name)+25);
+    sprintf(gfa_name, "%s.ct_flt", file_name);
+    FILE* fp = fopen(gfa_name, "w");
+	if (!fp) {
+		free(gfa_name);
+        return 0;
+    }
+	ha_ct_t* ct_idx = (ha_ct_t*)i_ct_idx;
+	int i;
+	ha_ct1_t *g;
+	fwrite(&ct_idx->k, sizeof(ct_idx->k), 1, fp);
+	fwrite(&ct_idx->pre, sizeof(ct_idx->pre), 1, fp);
+	fwrite(&ct_idx->n_hash, sizeof(ct_idx->n_hash), 1, fp);
+	fwrite(&ct_idx->n_shift, sizeof(ct_idx->n_shift), 1, fp);
+	fwrite(&ct_idx->tot, sizeof(ct_idx->tot), 1, fp);
+	for (i = 0; i < 1<<ct_idx->pre; i++)
+	{
+		g = &(ct_idx->h[i]);
+		yak_ct_save(g->h, fp);
+	}
+
+
+	fprintf(stderr, "[M::%s] Index has been written.\n", __func__);
+    free(gfa_name);
+	fclose(fp);
+	return 1;
+}
+
+int load_ct_index(void **i_ct_idx, char* file_name)
+{
+	char* gfa_name = (char*)malloc(strlen(file_name)+25);
+    sprintf(gfa_name, "%s.ct_flt", file_name);
+    FILE* fp = fopen(gfa_name, "r");
+	if (!fp) {
+		free(gfa_name);
+        return 0;
+    }
+	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);
+	CALLOC(h->h, 1<<h->pre);
+
+
+	index_time = yak_realtime();
+	for (i = 0; i < 1<<h->pre; ++i) 
+	{
+		g = &(h->h[i]);
+		yak_ct_load(&(g->h), fp);
+	}
+
+	(*ct_idx) = h;
+	fprintf(stderr, "[M::%s::%.3f] ==> Loaded count table\n", __func__, yak_realtime() - index_time);
+	fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__);
+	free(gfa_name);
+	return 1;
+}
+
+int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
 {
     char* gfa_name = (char*)malloc(strlen(file_name)+25);
     sprintf(gfa_name, "%s.pt_flt", file_name);
@@ -970,7 +1127,7 @@ int write_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt
 	return 1;
 }
 
-int load_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
+int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name)
 {
 	char* gfa_name = (char*)malloc(strlen(file_name)+25);
     sprintf(gfa_name, "%s.pt_flt", file_name);
@@ -1074,4 +1231,4 @@ int load_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t
 
 	free(gfa_name);
 	return 1;
-}
+}

htab.h

@@ -29,6 +29,8 @@ extern void *ha_flt_tab;
 extern ha_pt_t *ha_idx;
 extern void *ha_flt_tab_hp;
 extern ha_pt_t *ha_idx_hp;
+extern void *ha_ct_table;
+
 
 void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode);
 int ha_ft_isflt(const void *hh, uint64_t y);
@@ -38,8 +40,11 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 void ha_pt_destroy(ha_pt_t *h);
 const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);
 
-int write_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
-int load_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
+int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
+int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
+int write_ct_index(void *ct_idx, char* file_name);
+int load_ct_index(void **ct_idx, char* file_name);
+int query_ct_index(void* ct_idx, uint64_t hash);
 
 ha_abuf_t *ha_abuf_init(void);
 void ha_abuf_destroy(ha_abuf_t *ab);
@@ -55,8 +60,9 @@ double yak_cpu_usage(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);
-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);
-int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het);
+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);
+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,13 +147,14 @@ 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; \
-		fread(&n_buckets, sizeof(n_buckets), 1, fp); \
-		fread(&(*h)->bits, sizeof((*h)->bits), 1, fp); \
-		fread(&(*h)->count, sizeof((*h)->count), 1, fp); \
+		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); \
 		(*h)->used = (khint32_t*)kmalloc(__kh_fsize(n_buckets) * sizeof(khint32_t)); \
 		(*h)->keys = (khkey_t*)kmalloc(n_buckets * sizeof(khkey_t)); \
-		fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
-		fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); \
+		flag += fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
+		flag += fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); \
 		return 1; \
 	} \
 

sketch.cpp

@@ -139,12 +139,15 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 
 
 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_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct)
 {	///in default, w = 51, k = 51, is_hpc = 1
 	/**
 	 uint64_t x;
 	 uint64_t rid:28, pos:27, rev:1, span:8;
 	 **/
+	extern void *ha_ct_table;
+	if(dbg_ct != NULL) dbg_ct->a.n = 0;
+
 	static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 };
 	uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0}, filtered;
 	int i, j, l, buf_pos, min_pos, kmer_span = 0;
@@ -180,10 +183,12 @@ kvec_t_u8_warp* k_flag)
 				}
 				tq_push(&tq, skip_len);
 				kmer_span += skip_len;
+				///how many bases that are covered by this HPC k-mer
+				///kmer_span includes at most k HPC elements
 				if (tq.count > k) kmer_span -= tq_shift(&tq);
 			} else kmer_span = l + 1 < k? l + 1 : k; 
 			///kmer_span should be used for HPC k-mer
-			///so for non-HPC k-mer, kmer_span should be k in any case?
+			///non-HPC k-mer, kmer_span should be k
 			///kmer_span is used to calculate anchor pos on reverse complementary strand
 
 			if(k_flag != NULL) k_flag->a.a[i] = 1;///lable all useful base, which are not ignored by HPC
@@ -201,6 +206,7 @@ kvec_t_u8_warp* k_flag)
 
 				filtered = 0;
 				if(hf != 0) filtered = ha_ft_isflt(hf, y);
+				if(dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i));
 				///if (hf == 0 || ha_ft_isflt(hf, y) == 0)
 				if(filtered == 0)
 					info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;