Merge branch 'master' into hifiasm-v0.14

merge README

Assembly.cpp

@@ -11,10 +11,12 @@
 #include "htab.h"
 #include "kthread.h"
 
-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);
+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, kvec_t_u64_warp* dbg_ct);
 void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list);
 
 All_reads R_INF;
+Debug_reads R_INF_FLAG;
 
 void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length, char* new_read, int* new_length)
 {
@@ -297,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)
 {
@@ -442,6 +420,7 @@ typedef struct {
 	UC_Read self_read, ovlp_read;
 	Candidates_list clist;
 	overlap_region_alloc olist;
+    overlap_region_alloc olist_hp;
 	ha_abuf_t *ab;
 	// error correction related buffers
 	int64_t num_read_base, num_correct_base, num_recorrect_base;
@@ -451,6 +430,10 @@ typedef struct {
 	Correct_dumy correct;
 	haplotype_evdience_alloc hap;
 	Round2_alignment round2;
+    kvec_t_u32_warp b_buf;
+    kvec_t_u64_warp r_buf;
+    kvec_t_u8_warp k_flag;
+    overlap_region tmp_region;
 } ha_ovec_buf_t;
 
 ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
@@ -462,6 +445,11 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
 	init_UC_Read(&b->ovlp_read);
 	init_Candidates_list(&b->clist);
 	init_overlap_region_alloc(&b->olist);
+    init_overlap_region_alloc(&b->olist_hp);
+    init_fake_cigar(&(b->tmp_region.f_cigar));
+    kv_init(b->b_buf.a);
+    kv_init(b->r_buf.a);
+    kv_init(b->k_flag.a);
 	b->ab = ha_abuf_init();
 	if (!b->is_final) {
 		init_Cigar_record(&b->cigar1);
@@ -480,7 +468,12 @@ void ha_ovec_destroy(ha_ovec_buf_t *b)
 	destory_UC_Read(&b->ovlp_read);
 	destory_Candidates_list(&b->clist);
 	destory_overlap_region_alloc(&b->olist);
+    destory_overlap_region_alloc(&b->olist_hp);
 	ha_abuf_destroy(b->ab);
+    destory_fake_cigar(&(b->tmp_region.f_cigar));
+    kv_destroy(b->b_buf.a);
+    kv_destroy(b->r_buf.a);
+    kv_destroy(b->k_flag.a);
 	if (!b->is_final) {
 		destory_Cigar_record(&b->cigar1);
 		destory_Graph(&b->POA_Graph);
@@ -510,6 +503,7 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b)
 {
 	int64_t i, mem = 0, mem_clist, mem_olist;
 	mem_clist = b->clist.size * sizeof(k_mer_hit) + b->clist.chainDP.size * 7 * 4;
+    
 	mem_olist = b->olist.size * sizeof(overlap_region);
 	for (i = 0; i < (int64_t)b->olist.size; ++i) {
 		const overlap_region *r = &b->olist.list[i];
@@ -517,6 +511,14 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b)
 		mem_olist += r->f_cigar.size * 8;
 		mem_olist += r->boundary_cigars.size * sizeof(window_list);
 	}
+    mem_olist += b->olist_hp.size * sizeof(overlap_region);
+	for (i = 0; i < (int64_t)b->olist_hp.size; ++i) {
+		const overlap_region *r = &b->olist_hp.list[i];
+		mem_olist += r->w_list_size * sizeof(window_list);
+		mem_olist += r->f_cigar.size * 8;
+		mem_olist += r->boundary_cigars.size * sizeof(window_list);
+	}
+
 	mem = ha_abuf_mem(b->ab) + mem_clist + mem_olist;
 	if (!b->is_final) {
 		mem += sizeof(Cigar_record) + b->cigar1.lost_base_size + b->cigar1.size * 4;
@@ -534,8 +536,8 @@ static void worker_ovec(void *data, long i, int tid)
 	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
 	int fully_cov, abnormal;
 
-	ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1);
+    ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 
-    ///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 1);
+    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);
@@ -557,7 +559,15 @@ static void worker_ovec(void *data, long i, int tid)
 	}
 	R_INF.paf[i].is_abnormal = abnormal;
 
-	if (b->save_ov) {
+    R_INF.trio_flag[i] = AMBIGU;
+    
+    ///need to be fixed in r305
+    // if(ha_idx_hp == NULL)
+    // {
+    //     R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &(b->k_flag), RESEED_HP_RATE, Get_READ_LENGTH(R_INF, i), NULL);
+    // }
+
+    if (R_INF.trio_flag[i] != AMBIGU || b->save_ov) {
 		int is_rev = (asm_opt.number_of_round % 2 == 0);
 		push_overlaps(&(R_INF.paf[i]), &b->olist, 1, &R_INF, is_rev);
 		push_overlaps(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, is_rev);
@@ -568,25 +578,111 @@ static void worker_ovec(void *data, long i, int tid)
 static void worker_ovec_related_reads(void *data, long i, int tid)
 {
 	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
-	int required_read_name_length = strlen(asm_opt.required_read_name);
-	uint64_t k;
-	if (required_read_name_length == (int)Get_NAME_LENGTH((R_INF),i)
-			&&
-			memcmp(asm_opt.required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
-	{
-		ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1);
-        ///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 1);
 
-		fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i));
+    uint64_t k, queryNameLen;
-		recover_UC_Read(&b->self_read, &R_INF, i);
+    for (k = 0; k < R_INF_FLAG.query_num; k++)
-		fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq);
+    {
+        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)
+        {
+            break;
+        }
+    }
 
-		for (k = 0; k < b->olist.length; k++) {
+    if(k < R_INF_FLAG.query_num)
-			fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
+    {
-			recover_UC_Read(&b->self_read, &R_INF, b->olist.list[k].y_id);
+        int fully_cov, abnormal, q_idx = k;
-			fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq);
+
+        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), &(R_INF_FLAG.candidate_count[q_idx]));
+
+        clear_Cigar_record(&b->cigar1);
+        clear_Round2_alignment(&b->round2);
+
+        correct_overlap(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon,
+                &b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal);
+
+        b->num_read_base += b->self_read.length;
+        b->num_correct_base += b->correct.corrected_base;
+        b->num_recorrect_base += b->round2.dumy.corrected_base;
+
+        push_cigar(R_INF.cigars, i, &b->cigar1);
+        push_cigar(R_INF.second_round_cigar, i, &b->round2.cigar);
+
+        R_INF.paf[i].is_fully_corrected = 0;
+        if (fully_cov) {
+            if (get_cigar_errors(&b->cigar1) == 0 && get_cigar_errors(&b->round2.cigar) == 0)
+                R_INF.paf[i].is_fully_corrected = 1;
+        }
+        R_INF.paf[i].is_abnormal = abnormal;
+
+        
+
+        pthread_mutex_lock(&R_INF_FLAG.OutputMutex);
+
+        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++) 
+        {
+            if(b->olist.list[k].is_match != 1) continue;
+			fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
+            fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n",
+            b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e, 
+            b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel,
+            Get_READ_LENGTH(R_INF, b->olist.list[k].y_id));
 		}
-	}
+
+        fprintf(R_INF_FLAG.fp, "***************************reverse ovlp***************************\n");
+        for (k = 0; k < b->olist.length; k++) 
+        {
+            if(b->olist.list[k].is_match != 2) continue;
+			fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
+            fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n",
+            b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e, 
+            b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel,
+            Get_READ_LENGTH(R_INF, b->olist.list[k].y_id));
+		}
+
+        fprintf(R_INF_FLAG.fp, "***************************unmatched ovlp***************************\n");
+        for (k = 0; k < b->olist.length; k++) 
+        {
+            if(b->olist.list[k].is_match == 1) continue;
+            if(b->olist.list[k].is_match == 2) continue;
+			fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id));
+            fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n",
+            b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e, 
+            b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel,
+            Get_READ_LENGTH(R_INF, b->olist.list[k].y_id));
+		}
+
+        R_INF.trio_flag[i] = AMBIGU;
+
+        ///need to be fixed in r305
+        // if(ha_idx_hp == NULL)
+        // {
+        //     R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &(b->k_flag), RESEED_HP_RATE, Get_READ_LENGTH(R_INF, i), R_INF_FLAG.fp);
+        // }
+        
+        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);
+    }
 }
 
 
@@ -693,25 +789,76 @@ void Output_corrected_reads()
     fclose(output_file);
 }
 
+void debug_print_pob_regions()
+{
+    uint64_t i, total = 0;
+    for (i = 0; i < R_INF.total_reads; i++)
+    {
+        if(R_INF.trio_flag[i]!=AMBIGU)
+        {
+            total++;
+            fprintf(stderr, "(%lu) %.*s\n", i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
+        }
+    }
+    fprintf(stderr, "total hp reads: %lu, R_INF.total_reads: %lu\n", total, R_INF.total_reads);
+    exit(1);
+}
+
+void rescue_hp_reads(ha_ovec_buf_t **b)
+{
+    int hom_cov, het_cov;
+    ha_flt_tab_hp = ha_idx_hp = NULL;
+    if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
+       ha_flt_tab_hp = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 1);
+    }
+    ha_idx_hp = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 1, &R_INF, &hom_cov, &het_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);
+
+
+
+
+    ha_ft_destroy(ha_flt_tab_hp); ha_flt_tab_hp = NULL;
+    ha_pt_destroy(ha_idx_hp); ha_idx_hp = NULL;
+}
+
+
+
+
 void ha_overlap_and_correct(int round)
 {
-	int i, hom_cov, het_cov;
+	int i, hom_cov, het_cov, r_out = 0;
 	ha_ovec_buf_t **b;
 	ha_ecsave_buf_t *e;
+    ha_flt_tab_hp = ha_idx_hp = NULL;
 
+    if((ha_idx == NULL)&&(asm_opt.flag & HA_F_VERBOSE_GFA)&&(round == asm_opt.number_of_round - 1))
+    {
+        r_out = 1;
+    }
+
+    if(asm_opt.required_read_name) init_Debug_reads(&R_INF_FLAG, asm_opt.required_read_name); // for debugging only
 	// overlap and correct reads
 	CALLOC(b, asm_opt.thread_num);
 	for (i = 0; i < asm_opt.thread_num; ++i)
 		b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1));
-	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF, &hom_cov, &het_cov); // build the index
+    if(ha_idx) hom_cov = asm_opt.hom_cov;
-	if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
+	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
+	///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_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
 	ha_pt_destroy(ha_idx);
-	ha_idx = 0;
+	ha_idx = NULL;
 
 	// collect statistics
 	for (i = 0; i < asm_opt.thread_num; ++i) {
@@ -723,8 +870,8 @@ void ha_overlap_and_correct(int round)
 	}
 	free(b);
 
-	if (asm_opt.required_read_name) exit(0); // for debugging only
+	if (asm_opt.required_read_name) destory_Debug_reads(&R_INF_FLAG), exit(0); // for debugging only
-
+    
 	// save corrected reads to R_INF
 	CALLOC(e, asm_opt.thread_num);
 	for (i = 0; i < asm_opt.thread_num; ++i) {
@@ -740,6 +887,7 @@ void ha_overlap_and_correct(int round)
 		free(e[i].second_round_read);
 	}
 	free(e);
+    ///debug_print_pob_regions();
 }
 
 
@@ -1138,15 +1286,8 @@ static void worker_ov_final(void *data, long i, int tid)
 	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
 
 	//get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
-	ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.001, asm_opt.max_n_chain, 0);
+    ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 0.001, 
-    ///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 0);
+    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);
-
-	/**
-	  correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon,
-	  &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1,
-	  &current_cigar, &hap, &second_round, 0, 0);
-	  push_final_overlaps(&(R_INF.paf[i]), &overlap_list);
-	 **/
 
 	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);
@@ -1211,7 +1352,8 @@ 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_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, HIGH_HET_ERROR_RATE, asm_opt.max_n_chain, 1);
+    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), 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);
@@ -1410,7 +1552,7 @@ void hap_recalculate_peaks(char* output_file_name)
     int hom_cov, het_cov;
     // construct hash table for high occurrence k-mers
     if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
-        ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov);
+        ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0);
         ha_opt_update_cov(&asm_opt, hom_cov);
     }
     free(R_INF.read_length);
@@ -1418,7 +1560,7 @@ void hap_recalculate_peaks(char* output_file_name)
 
     load_All_reads(&R_INF, gfa_name);
 
-    ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index
+    ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
     asm_opt.hom_cov = hom_cov;
     asm_opt.het_cov = het_cov;
     ha_pt_destroy(ha_idx);
@@ -1435,10 +1577,12 @@ void ha_overlap_final(void)
 {
 	int i, hom_cov, het_cov;
 	ha_ovec_buf_t **b;
+    ha_flt_tab_hp = ha_idx_hp = NULL;
+
 	CALLOC(b, asm_opt.thread_num);
 	for (i = 0; i < asm_opt.thread_num; ++i)
 		b[i] = ha_ovec_init(asm_opt.flag & HA_F_HIGH_HET, 1);///b[i] = ha_ovec_init(1, 1);
-	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index
+	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index
     if(asm_opt.flag & HA_F_HIGH_HET)
     {
         kt_for(asm_opt.thread_num, worker_ov_final_high_het, b, R_INF.total_reads);
@@ -1457,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);
@@ -1475,14 +1621,18 @@ int ha_assemble(void)
         if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2;
 	}
 	if (!ovlp_loaded) {
+        ha_flt_tab = ha_idx = NULL;
+        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)) {
+		if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL) 
-			ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov);
+        {
+			ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0);
 			ha_opt_update_cov(&asm_opt, hom_cov);
 		}
 		// error correction
 		assert(asm_opt.number_of_round > 0);
-		for (r = 0; r < asm_opt.number_of_round; ++r) {
+		for (r = ha_idx?asm_opt.number_of_round-1:0; r < asm_opt.number_of_round; ++r) {
 			ha_opt_reset_to_round(&asm_opt, r); // this update asm_opt.roundID and a few other fields
 			ha_overlap_and_correct(r);
 			fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(),

Assembly.h

@@ -8,6 +8,11 @@
 #define Get_Cigar_Type(RECORD) (RECORD&3)
 #define Get_Cigar_Length(RECORD) (RECORD>>2)
 
+#define RESEED_DP 4
+#define RESEED_PEAK_RATE 0.15
+#define RESEED_LEN 2000
+#define RESEED_HP_RATE 0.9
+
 int ha_assemble(void);
 
 #endif

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->min_hist_kmer_cnt = 5;
-    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;
@@ -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)
+    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);
-        free(asm_opt->read_file_names);
+    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,11 +368,56 @@ 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);
@@ -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.10-r299"
+#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 min_hist_kmer_cnt;
-    int k_mer_max_freq;
     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.cpp

@@ -7,8 +7,12 @@
 #include "Assembly.h"
 #include "CommandLines.h"
 #include "ksw2.h"
+#include "ksort.h"
 
+#define generic_key(x) (x)
+KRADIX_SORT_INIT(b32, uint32_t, generic_key, 4)
 
+int ha_ov_type(const overlap_region *r, uint32_t len);
 
 
 void clear_Round2_alignment(Round2_alignment* h)
@@ -7202,7 +7206,175 @@ void partition_overlaps_advance(overlap_region_alloc* overlap_list, All_reads* R
 }
 
 
+void collect_no_cov_regions(overlap_region_alloc* overlap_list, All_reads* R_INF, 
+kvec_t_u32_warp* b, kvec_t_u64_warp* r, int min_dp, int min_len)
+{
+    b->a.n = r->a.n = 0;
+    ///if(overlap_list->length == 0) return;
+    long long i = 0, xLen = Get_READ_LENGTH((*R_INF), overlap_list->list[0].x_id);
+    uint32_t qs, qe;
+    uint64_t tmp;
+    int dp, old_dp, s_start = 0, s_end = 0;
+    ///at least 1
+    if(min_len < 1) min_len = 1;
 
+    
+    for (i = 0; i < (long long)overlap_list->length; i++)
+    {
+        if (overlap_list->list[i].is_match != 1 && overlap_list->list[i].is_match != 2) continue;
+
+        qs = overlap_list->list[i].x_pos_s;
+        qe = overlap_list->list[i].x_pos_e + 1;
+        kv_push(uint32_t, b->a, qs<<1);
+        kv_push(uint32_t, b->a, qe<<1|1);
+    }
+
+
+    ///we can identify the qs and qe by the 0-th bit
+    radix_sort_b32(b->a.a, b->a.a + b->a.n);
+
+    for (i = 0, dp = 0; i < (long long)b->a.n; ++i) 
+    {
+        old_dp = dp;
+        //if a[j] is qe
+        if (b->a.a[i]&1) --dp;
+        else ++dp;
+        /**
+        min_dp is the coverage drop threshold
+        there are two cases: 
+            1. old_dp = dp + 1 (b.a[j] is qe); 2. old_dp = dp - 1 (b.a[j] is qs);
+        **/
+        if (old_dp < min_dp && dp >= min_dp) ///old_dp < dp, b.a[j] is qs
+        { 
+            ///case 2, a[j] is qs
+            s_end = b->a.a[i]>>1;
+            ///at least 1
+            if(s_end-s_start >= min_len)
+            {
+                tmp = s_start; tmp = tmp << 32; tmp = tmp | (uint64_t)(s_end-1);
+                kv_push(uint64_t, r->a, tmp);
+            }
+        }
+        else if (old_dp >= min_dp && dp < min_dp) ///old_dp > min_dp, b.a[j] is qe
+        {
+            s_start = b->a.a[i]>>1;
+        }
+    }
+
+
+    if(s_start < xLen && xLen-s_start >= min_len)
+    {
+        s_end = xLen;
+        tmp = s_start; tmp = tmp << 32; tmp = tmp | (uint64_t)(s_end-1);
+        kv_push(uint64_t, r->a, tmp);
+    }
+}
+
+int collect_hp_regions_back(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u32_warp* b, kvec_t_u64_warp* r, kvec_t_u8_warp* k_flag, float hp_rate, FILE* fp)
+{
+    int i, k, qs, qe, ava_k_mer = 0, hp_k_mer = 0, min_dp;
+    // min_dp = RESEED_DP;
+    // if(asm_opt.hom_cov > 0) min_dp = asm_opt.hom_cov * RESEED_PEAK_RATE;
+    // if(min_dp > RESEED_DP) min_dp = RESEED_DP;
+    min_dp = RESEED_DP;
+    if(asm_opt.hom_cov > 0) min_dp = asm_opt.hom_cov * RESEED_PEAK_RATE;
+    if(asm_opt.het_cov > 0) min_dp = asm_opt.het_cov * RESEED_PEAK_RATE;
+    collect_no_cov_regions(olist, R_INF, b, r, min_dp, RESEED_LEN);
+
+    for (i = 0; i < (int)r->a.n; i++)
+    {
+        ///[qs, qe]
+        qs = r->a.a[i]>>32;
+        qe = (r->a.a[i]<<32)>>32;
+        
+        for (k = qs; k <= qe; k++)
+        {
+            if(k_flag->a.a[k] > 1) ava_k_mer++;
+            if(k_flag->a.a[k] > 2) hp_k_mer++;
+        }
+
+        if(fp) fprintf(fp, "qs: %d, qe: %d, ava_k_mer: %d, hp_k_mer: %d\n", qs, qe, ava_k_mer, hp_k_mer);        
+    }
+
+    if(fp) fprintf(fp, "ava_k_mer: %d, hp_k_mer: %d, hp_rate: %f, min_dp: %d, a.n: %d\n", ava_k_mer, hp_k_mer, hp_rate, min_dp, (int)r->a.n);
+
+    // if(fp)
+    // {
+    //     for (k = 0; k < (int)k_flag->a.n; k++)
+    //     {
+    //         if(k_flag->a.a[k] > 0) fprintf(fp, "(%d) %u\n", k, k_flag->a.a[k]);
+    //     }
+    // }
+    
+    if(hp_k_mer > ava_k_mer*hp_rate) return 1; ///must use '>' instead of '>='
+    r->a.n = 0;
+    return 0;
+}
+
+
+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)
+{
+    int i, ava_k_mer = 0, hp_k_mer = 0, vLen, min_dp;
+    int32_t w, n[4];
+    n[0] = n[1] = n[2] = n[3] = 0;
+    min_dp = RESEED_DP;
+    if(asm_opt.hom_cov > 0) min_dp = asm_opt.hom_cov * RESEED_PEAK_RATE;
+    if(min_dp > RESEED_DP) min_dp = RESEED_DP;
+    overlap_region* ov = NULL;
+
+    for (i = 0; i < (long long)olist->length; i++)
+    {
+        ov = &(olist->list[i]);
+        if (ov->is_match != 1 && ov->is_match != 2) continue;
+
+        w = ha_ov_type(ov, rlen);
+		++n[w];
+    }
+
+    if(fp) fprintf(fp, "n[0]: %d, n[1]: %d, n[2]: %d, n[3]: %d\n", n[0], n[1], n[2], n[3]);
+
+    // n[0] += n[2];
+    // n[1] += n[2];
+
+    if(n[0] < min_dp)
+    {
+        ava_k_mer = hp_k_mer = 0;
+        vLen = MIN(k_flag->a.n, RESEED_LEN);
+        for (i = 0; i < vLen; i++)
+        {
+            if(k_flag->a.a[i] > 1) ava_k_mer++;
+            if(k_flag->a.a[i] > 2) hp_k_mer++;
+        }
+        if(hp_k_mer > ava_k_mer*hp_rate) return 1;
+    }
+
+
+    if(n[1] < min_dp)
+    {
+        ava_k_mer = hp_k_mer = 0;
+        vLen = MIN(k_flag->a.n, RESEED_LEN);
+        for (i = k_flag->a.n - vLen; i < (int)k_flag->a.n; i++)
+        {
+            if(k_flag->a.a[i] > 1) ava_k_mer++;
+            if(k_flag->a.a[i] > 2) hp_k_mer++;
+        }
+        if(hp_k_mer > ava_k_mer*hp_rate) return 1;
+    }
+
+
+    if(fp) fprintf(fp, "ava_k_mer: %d, hp_k_mer: %d, hp_rate: %f, min_dp: %d\n", ava_k_mer, hp_k_mer, hp_rate, min_dp);
+
+    // if(fp)
+    // {
+    //     for (k = 0; k < (int)k_flag->a.n; k++)
+    //     {
+    //         if(k_flag->a.a[k] > 0) fprintf(fp, "(%d) %u\n", k, k_flag->a.a[k]);
+    //     }
+    // }
+    
+    return 0;
+}
 
 void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, 
                         UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, 

Correct.h

@@ -21,8 +21,8 @@
 
 ///#define FLAG_THRE 0
 
-#define MAX(x, y) ((x >= y)?(x):(y))
+#define MAX(x, y) (((x) >= (y))?(x):(y))
-#define MIN(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 
@@ -1173,6 +1173,21 @@ void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_I
                         UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read);
 long long get_affine_gap_score(overlap_region* ovc, UC_Read* g_read, UC_Read* overlap_read, uint8_t* x_num,
 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

Hash_Table.cpp

@@ -305,6 +305,9 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp)
     }
 }
 
+
+
+
 long long get_chainLen(long long x_beg, long long x_end, long long xLen, 
                        long long y_beg, long long y_end, long long yLen)
 {
@@ -337,6 +340,50 @@ long long get_chainLen(long long x_beg, long long x_end, long long xLen,
     return x_end - x_beg + 1;
 }
 
+
+void debug_chain_single_site(k_mer_hit* a, long long a_n, Chain_Data* dp, int x_readLen, int y_readLen, int s_index)
+{
+    long long j, current_j = s_index;
+    long long selfLen = 0, indels = 0;
+    long long distance_self_pos, distance_pos, distance_gap;
+
+    j = s_index;
+    while (j >= 0)
+    {
+        current_j = j;
+
+        j = dp->pre[j];
+
+        if(j != -1)
+        {
+            distance_self_pos = a[current_j].self_offset - a[j].self_offset;
+            distance_pos = a[current_j].offset - a[j].offset;
+            distance_gap = distance_pos > distance_self_pos? distance_pos - distance_self_pos : distance_self_pos - distance_pos;
+
+            indels += distance_gap; 
+            selfLen += distance_self_pos;
+        }
+        fprintf(stderr, "j: %lld, score: %lld, occ: %d, pre_j: %lld\n", 
+                    current_j, (long long)dp->score[current_j], dp->occ[current_j], j);
+    }
+
+    fprintf(stderr, "s_self_offset: %u, s_offset: %u, e_self_offset: %u, e_offset: %u, ovlp length: %lld, x_readLen: %d, y_readLen: %d\n", 
+    a[s_index].self_offset, a[s_index].offset, a[current_j].self_offset, a[current_j].offset,
+    get_chainLen(a[s_index].self_offset, a[current_j].self_offset, x_readLen, 
+                                     a[s_index].offset, a[current_j].offset, y_readLen), x_readLen, y_readLen);
+
+    if(indels != dp->indels[s_index])
+    {
+        fprintf(stderr, "indels: %lld, dp->indels[i]: %ld\n", indels, (long)dp->indels[s_index]);
+    }
+
+    if(selfLen != dp->self_length[s_index])
+    {
+        fprintf(stderr, "selfLen: %lld, dp->self_length[i]: %ld\n", selfLen, (long)dp->self_length[s_index]);
+    }
+    fprintf(stderr,"\n");
+}
+
 int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc, double bw_thres)
 {
 	int32_t i, tot_indel = 0, tot_len = 0;
@@ -348,7 +395,7 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc
 	if (i < n_a) return -1;
 	bw_pen = 1.0 / bw_thres;
 	dp->score[0] = a[0].good? min_sc : min_sc>>1;
-	dp->pre[0] = -1, dp->indels[0] = 0, dp->self_length[0] = 0;
+	dp->pre[0] = -1, dp->indels[0] = 0, dp->self_length[0] = 0, dp->occ[0] = 1;
 	for (i = 1; i < n_a; ++i) {
 		int32_t score, dg;
 		int32_t dx = (int32_t)a[i].offset - (int32_t)a[i-1].offset;
@@ -368,6 +415,7 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc
 		dp->pre[i] = i - 1;
 		dp->indels[i] = tot_indel;
 		dp->self_length[i] = tot_len;
+        dp->occ[i] = i + 1;
 	}
 	if (i < n_a) return -1;
 	return n_a;
@@ -391,6 +439,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
     resize_Chain_Data(dp, a_n);
 
 	ret = ha_chain_check(a, a_n, dp, min_score, band_width_threshold);
+
 	if (ret > 0) {
 		a_n = ret;
 		goto skip_dp;
@@ -435,7 +484,9 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
             ///min distance
             distance_min = distance_pos < distance_self_pos? distance_pos:distance_self_pos;
             score = distance_min < min_score? distance_min : min_score;
-			if (!a[j].good) score >>= 1;
+            ///need to be fixed in r305
+			///if (!a[j].good) score = (score >> 1) + (score & 1);
+            if (!a[j].good) score >>= 1;
 
             gap_rate = (double)((double)(total_indels)/(double)(total_self_length));
             ///if the gap rate > 0.06, score will be negative
@@ -444,7 +495,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
             score += dp->score[j];
 
 			///find a new max score
-			if (score > max_score) {
+			if (score > max_score) {///must use > instead of >=
 				max_score = score;
 				max_j = j;
 				max_indels = total_indels;
@@ -466,10 +517,12 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
         dp->pre[i] = max_j;
         dp->indels[i] = max_indels;
         dp->self_length[i] = max_self_length;
+        dp->occ[i] = 1;
+        if(max_j != -1) dp->occ[i] = dp->occ[max_j] + 1;
     }
 
     ///debug_chain(a, a_n, dp);
-
+    
 skip_dp:
 
     max_score = -1;
@@ -562,7 +615,7 @@ skip_dp:
     }
 }
 
-void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, 
+void calculate_overlap_region_by_chaining_back(Candidates_list* candidates, overlap_region_alloc* overlap_list, 
                                           uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end)
 {
     overlap_region tmp_region;
@@ -627,6 +680,71 @@ void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_r
     destory_fake_cigar(&(tmp_region.f_cigar));
 }
 
+
+void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx,
+                                          uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end, overlap_region* f_cigar)
+{
+    long long i = 0;
+    uint64_t current_ID;
+    uint64_t current_stand;
+
+    if (candidates->length == 0)
+    {
+        return;
+    }
+
+    long long sub_region_beg;
+    long long sub_region_end;
+
+    clear_fake_cigar(&((*f_cigar).f_cigar));
+
+    i = 0;
+    while (i < candidates->length)
+    {
+        chain_idx->a.n = 0;
+        current_ID = candidates->list[i].readID;
+        current_stand = candidates->list[i].strand;
+
+        ///reference read
+        (*f_cigar).x_id = readID;
+        (*f_cigar).x_pos_strand = current_stand;
+        ///query read
+        (*f_cigar).y_id = current_ID;
+        ///here the strand of query is always 0
+        (*f_cigar).y_pos_strand = 0;  
+
+        sub_region_beg = i;
+        sub_region_end = i;
+        i++;
+
+        while (i < candidates->length 
+        && 
+        current_ID == candidates->list[i].readID
+        &&
+        current_stand == candidates->list[i].strand)
+        {
+            sub_region_end = i;
+            i++;
+        }
+
+        if ((*f_cigar).x_id == (*f_cigar).y_id)
+        {
+            continue;
+        }
+
+		chain_DP(candidates->list + sub_region_beg,
+				sub_region_end - sub_region_beg + 1, &(candidates->chainDP), f_cigar, band_width_threshold,
+				25, Get_READ_LENGTH((*R_INF), (*f_cigar).x_id), Get_READ_LENGTH((*R_INF), (*f_cigar).y_id));
+
+        ///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1)
+        if ((*f_cigar).x_id != (*f_cigar).y_id)
+        {
+            append_inexact_overlap_region_alloc(overlap_list, f_cigar, R_INF, add_beg_end);
+        }
+    }
+}
+
+
 void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error,
                         int extra_begin, int extra_end, int error_threshold)
 {
@@ -704,6 +822,7 @@ void destory_Chain_Data(Chain_Data* x)
     free(x->pre);
     free(x->indels);
     free(x->self_length);
+    free(x->occ);
 	free(x->tmp);
 }
 
@@ -716,6 +835,7 @@ void resize_Chain_Data(Chain_Data* x, long long size)
 		REALLOC(x->pre, x->size);
 		REALLOC(x->indels, x->size);
 		REALLOC(x->self_length, x->size);
+        REALLOC(x->occ, x->size);
 		REALLOC(x->tmp, x->size);
 	}
 }

Hash_Table.h

@@ -119,6 +119,7 @@ typedef struct {
 	int64_t *pre;
 	int32_t *indels;
 	int32_t *self_length;
+    int32_t *occ;
 	int64_t *tmp; // MUST BE 64-bit integer
 	int64_t length;
 	int64_t size;
@@ -144,8 +145,8 @@ int extra_begin, int extra_end, int error_threshold);
 
 void overlap_region_sort_y_id(overlap_region *a, long long n);
 
-void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, 
+void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx,
-uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end);
+uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end, overlap_region* f_cigar);
 
 void init_fake_cigar(Fake_Cigar* x);
 void destory_fake_cigar(Fake_Cigar* x);

Makefile

@@ -1,10 +1,12 @@
 CXX=		g++
+CC=			gcc
 CXXFLAGS=	-g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall
+CFLAGS=		$(CXXFLAGS)
 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
 
@@ -13,12 +15,15 @@ ifneq ($(asan),)
 	LIBS+=-fsanitize=address
 endif
 
-.SUFFIXES:.cpp .o
+.SUFFIXES:.cpp .c .o
 .PHONY:all clean depend
 
 .cpp.o:
 		$(CXX) -c $(CXXFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@
 
+.c.o:
+		$(CC) -c $(CFLAGS) $(CPPFLAGS) $(INCLUDES) $< -o $@
+
 all:$(EXE)
 
 $(EXE):$(OBJS) main.o
@@ -36,8 +41,9 @@ Assembly.o: Assembly.h CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h
 Assembly.o: Hash_Table.h htab.h POA.h Correct.h Levenshtein_distance.h
 Assembly.o: kthread.h
 CommandLines.o: CommandLines.h ketopt.h
-Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h ksw2.h
+Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h
 Correct.o: kdq.h CommandLines.h Levenshtein_distance.h POA.h Assembly.h
+Correct.o: ksw2.h
 Hash_Table.o: Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h
 Hash_Table.o: CommandLines.h ksort.h
 Levenshtein_distance.o: Levenshtein_distance.h
@@ -65,4 +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
-ksw2_extz2_sse.o: ksw2.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, 
+void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int id, char* command);
-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

@@ -46,6 +46,7 @@ void destory_All_reads(All_reads* r)
 		if (r->read_sperate[i]) free(r->read_sperate[i]);
 		if (r->paf && r->paf[i].buffer) free(r->paf[i].buffer);
 		if (r->reverse_paf && r->reverse_paf[i].buffer) free(r->reverse_paf[i].buffer);
+		///if (r->pb_regions) kv_destroy(r->pb_regions[i].a);
 	}
 	free(r->paf);
 	free(r->reverse_paf);
@@ -55,6 +56,7 @@ void destory_All_reads(All_reads* r)
 	free(r->name_index);
 	free(r->read_length);
 	free(r->trio_flag);
+	///if (r->pb_regions) free(r->pb_regions);
 }
 
 void write_All_reads(All_reads* r, char* read_file_name)
@@ -192,6 +194,7 @@ int load_All_reads(All_reads* r, char* read_file_name)
 		r->second_round_cigar[i].lost_base_length = r->cigars[i].lost_base_length = 0;
 		r->second_round_cigar[i].lost_base = r->cigars[i].lost_base = NULL;
 	}
+	///r->pb_regions = NULL;
 
     free(index_name);    
     fclose(fp);
@@ -263,6 +266,8 @@ void malloc_All_reads(All_reads* r)
 	r->second_round_cigar = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
 	r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
 	r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
+	///r->pb_regions = (kvec_t_u64_warp*)malloc(r->total_reads*sizeof(kvec_t_u64_warp));
+
 	for (i = 0; i < (long long)r->total_reads; i++)
 	{
 		r->second_round_cigar[i].size = r->cigars[i].size = 0;
@@ -274,6 +279,7 @@ void malloc_All_reads(All_reads* r)
 		r->second_round_cigar[i].lost_base = r->cigars[i].lost_base = NULL;
 		init_ma_hit_t_alloc(&(r->paf[i]));
 		init_ma_hit_t_alloc(&(r->reverse_paf[i]));
+		///kv_init(r->pb_regions[i].a);
 	}
 
 	r->name = (char*)malloc(sizeof(char)*r->total_name_length);
@@ -332,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)
 {
@@ -509,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);
@@ -667,3 +684,54 @@ void reverse_complement(char* pattern, uint64_t length)
 		pattern[end] = RC_CHAR(pattern[end]);
 	}
 }
+
+
+void init_Debug_reads(Debug_reads* x, const char* file)
+{
+	int nameLen, i, bufLen = 1000;
+	if((uint64_t)(bufLen) < strlen(file) + 50) bufLen = strlen(file) + 50;
+	char* Name_Buffer = (char*)malloc(sizeof(char)*bufLen);
+	fprintf(stderr, "Queried debugging reads at: %s\n", file);
+
+	x->fp = fopen(file,"r");
+	x->query_num = 0;
+	
+	while(fgets(Name_Buffer, bufLen, x->fp))
+	{
+		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;
+	while(fgets(Name_Buffer, bufLen, x->fp))
+	{
+		nameLen = strlen(Name_Buffer) - 1;
+		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++;
+	}
+	
+    fclose(x->fp);
+
+	sprintf(Name_Buffer, "%s.debug.stdout", file);
+	x->fp = fopen(Name_Buffer,"w");
+	fprintf(stderr, "Print debugging information to: %s\n", Name_Buffer);
+	free(Name_Buffer);
+}
+
+void destory_Debug_reads(Debug_reads* x)
+{
+	uint64_t i;
+	for (i = 0; i < x->query_num; i++)
+	{
+		free(x->read_name[i]);
+		kv_destroy(x->candidate_count[i].a);
+	}
+	
+	free(x->read_name);
+	fclose(x->fp);
+}

Process_Read.h

@@ -22,7 +22,8 @@
 ///#define Get_READ(R_INF, ID) R_INF.read + (R_INF.index[ID]>>2) + ID
 #define Get_READ(R_INF, ID) (R_INF).read_sperate[(ID)]
 #define Get_NAME(R_INF, ID) ((R_INF).name + (R_INF).name_index[(ID)])
-
+#define CHECK_BY_NAME(R_INF, NAME, ID) (Get_NAME_LENGTH((R_INF),(ID))==strlen((NAME)) && \
+                                        memcmp((NAME), Get_NAME((R_INF), (ID)), Get_NAME_LENGTH((R_INF),(ID))) == 0)
 
 extern uint8_t seq_nt6_table[256];
 extern char bit_t_seq_table[256][4];
@@ -93,6 +94,7 @@ typedef struct
 	uint32_t new_length;
 } Compressed_Cigar_record;
 
+
 #define AMBIGU 0
 #define FATHER 1
 #define MOTHER 2
@@ -128,6 +130,8 @@ typedef struct
 
     ma_hit_t_alloc* paf;
     ma_hit_t_alloc* reverse_paf;
+
+    ///kvec_t_u64_warp* pb_regions;
 } All_reads;
 
 extern All_reads R_INF;
@@ -140,6 +144,15 @@ typedef struct
 	long long RID;
 } UC_Read;
 
+typedef struct
+{
+    char** read_name;
+    uint64_t query_num;
+    kvec_t_u64_warp* candidate_count;
+    FILE* fp;
+    pthread_mutex_t OutputMutex;
+} Debug_reads;
+
 void init_All_reads(All_reads* r);
 void malloc_All_reads(All_reads* r);
 void ha_insert_read_len(All_reads *r, int read_len, int name_len);
@@ -154,5 +167,7 @@ void write_All_reads(All_reads* r, char* read_file_name);
 int load_All_reads(All_reads* r, char* read_file_name);
 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

@@ -28,8 +28,8 @@ typedef struct {
 } seed1_t;
 
 struct ha_abuf_s {
-	uint64_t n_a, m_a;
+	uint64_t n_a, m_a;///number of anchors (seed positions)
-	uint32_t old_mz_m;
+	uint32_t old_mz_m;///number of seeds
 	ha_mz1_v mz;
 	seed1_t *seed;
 	anchor1_t *a;
@@ -50,17 +50,16 @@ uint64_t ha_abuf_mem(const ha_abuf_t *ab)
 	return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(ha_mz1_t) + sizeof(seed1_t)) + sizeof(ha_abuf_t);
 }
 
-static int ha_ov_type(const overlap_region *r, uint32_t len)
+int ha_ov_type(const overlap_region *r, uint32_t len)
 {
 	if (r->x_pos_s == 0 && r->x_pos_e == len - 1) return 2; // contained in a longer read
 	else if (r->x_pos_s > 0 && r->x_pos_e < len - 1) return 3; // containing a shorter read
 	else return r->x_pos_s == 0? 0 : 1;
 }
 
-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)
+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* dbg_ct)
 {
-	extern void *ha_flt_tab;
-	extern ha_pt_t *ha_idx;
 	uint32_t i, rlen;
 	uint64_t k, l;
 	double low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
@@ -74,7 +73,8 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	rlen = Get_READ_LENGTH(R_INF, rid); // read length
 
 	// get the list of anchors
-	ha_sketch(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);
+	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;
 		REALLOC(ab->seed, ab->old_mz_m);
@@ -93,6 +93,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	}
 	for (i = 0, k = 0; i < ab->mz.n; ++i) {
 		int j;
+		///z is one of the minimizer
 		ha_mz1_t *z = &ab->mz.a[i];
 		seed1_t *s = &ab->seed[i];
 		for (j = 0; j < s->n; ++j) {
@@ -116,6 +117,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 		}
 	}
 
+
 	// copy over to _cl_
 	if (ab->m_a >= (uint64_t)cl->size) {
 		cl->size = ab->m_a;
@@ -131,7 +133,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	}
 	cl->length = ab->n_a;
 
-	calculate_overlap_region_by_chaining(cl, overlap_list, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain);
+	calculate_overlap_region_by_chaining(cl, overlap_list, chain_idx, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain, f_cigar);
 
 	#if 0
 	if (overlap_list->length > 0) {
@@ -172,10 +174,130 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 		}
 	}
 
-	ks_introsort_or_xs(overlap_list->length, overlap_list->list);
+	///ks_introsort_or_xs(overlap_list->length, overlap_list->list);
 }
 
 
+void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf)
+{
+	uint64_t j = 0, inner_j = 0;
+	while (j < overlap_list->length && inner_j < paf->length)
+    {
+        if(overlap_list->list[j].y_id < paf->buffer[inner_j].tn)
+        {
+            j++;
+        }
+        else if(overlap_list->list[j].y_id > paf->buffer[inner_j].tn)
+        {
+            inner_j++;
+        }
+        else
+        {
+            if(overlap_list->list[j].y_pos_strand == paf->buffer[inner_j].rev)
+            {
+				overlap_list->list[j].is_match = 1;
+            }
+            j++;
+            inner_j++;
+        }
+    }
+}
+
+
+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, 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, dbg_ct);
+
+	if(ha_idx_hp)
+	{
+		uint32_t i, k, y_id, overlapLen, max_i;
+		int shared_seed;
+		overlap_region t;
+		overlap_region_sort_y_id(overlap_list->list, overlap_list->length);
+		ma_hit_sort_tn(paf->buffer, paf->length);
+		ma_hit_sort_tn(rev_paf->buffer, rev_paf->length);
+		lable_matched_ovlp(overlap_list, paf);
+		lable_matched_ovlp(overlap_list, rev_paf);
+
+		for (i = 0, k = 0; i < overlap_list->length; ++i) 
+		{
+			if(overlap_list->list[i].is_match == 1)
+			{
+				if(k != i)
+				{
+					t = overlap_list->list[k];
+					overlap_list->list[k] = overlap_list->list[i];
+					overlap_list->list[i] = t;
+					overlap_list->list[k].is_match = 0;
+				}
+				k++;
+			}
+		}
+		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, dbg_ct);	
+		
+		if(overlap_list->length + overlap_list_hp->length > overlap_list->size)
+		{
+			overlap_list->list = (overlap_region*)realloc(overlap_list->list, 
+				sizeof(overlap_region)*(overlap_list->length + overlap_list_hp->length));
+			memset(overlap_list->list + overlap_list->size, 0, sizeof(overlap_region)*
+				(overlap_list->length + overlap_list_hp->length - overlap_list->size));
+			overlap_list->size = overlap_list->length + overlap_list_hp->length;
+		}
+		
+		for (i = 0, k = overlap_list->length; i < overlap_list_hp->length; i++, k++)
+		{
+			t = overlap_list->list[k];
+			overlap_list->list[k] = overlap_list_hp->list[i];
+			overlap_list_hp->list[i] = t;
+		}
+		overlap_list->length = k;
+		
+		overlap_region_sort_y_id(overlap_list->list, overlap_list->length);
+
+		i = k = 0;
+		while (i < overlap_list->length)
+		{
+			y_id = overlap_list->list[i].y_id;
+			shared_seed = overlap_list->list[i].shared_seed;
+			overlapLen = overlap_list->list[i].overlapLen;
+			max_i = i;
+			i++;
+			while (i < overlap_list->length && overlap_list->list[i].y_id == y_id)
+			{
+				if((overlap_list->list[i].shared_seed > shared_seed) || 
+				  ((overlap_list->list[i].shared_seed == shared_seed) && (overlap_list->list[i].overlapLen <= overlapLen)))
+				{
+					y_id = overlap_list->list[i].y_id;
+					shared_seed = overlap_list->list[i].shared_seed;
+					overlapLen = overlap_list->list[i].overlapLen;
+					max_i = i;
+				}
+				i++;
+			}
+
+			if(k != max_i)
+			{
+				t = overlap_list->list[k];
+				overlap_list->list[k] = overlap_list->list[max_i];
+				overlap_list->list[max_i] = t;
+			}
+			k++;
+		}
+
+		overlap_list->length = k;
+	}
+	
+	ks_introsort_or_xs(overlap_list->length, overlap_list->list);
+}
 
 void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list)
 {

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

@@ -35,6 +35,9 @@ const unsigned char seq_nt4_table[256] = { // translate ACGT to 0123
 
 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 *
@@ -46,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)
@@ -71,7 +75,7 @@ typedef struct {
 	int n_shift, n_hashes;
 	uint8_t *b;
 } yak_bf_t;
-
+///in most cases, n_shift = 25, n_hashes = 4
 yak_bf_t *yak_bf_init(int n_shift, int n_hashes)
 {
 	yak_bf_t *b;
@@ -126,10 +130,12 @@ typedef struct {
 
 typedef struct {
 	int k, pre, n_hash, n_shift;
-	uint64_t tot;
+	uint64_t tot;  ///number of distinct k-mers
 	ha_ct1_t *h;
 } ha_ct_t;
 
+///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;
@@ -138,12 +144,15 @@ static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift)
 	CALLOC(h, 1);
 	h->k = k, h->pre = pre;
 	CALLOC(h->h, 1<<h->pre);
+	///i<h->pre = 4096
+	///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, 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)
-			h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash);
+			h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash); ///h->n_shift = 37, h->pre = 12, h->n_hash = 4 
 	}
 	return h;
 }
@@ -173,15 +182,23 @@ static int ha_ct_insert_list(ha_ct_t *h, int create_new, int n, const uint64_t *
 	int j, mask = (1<<h->pre) - 1, n_ins = 0;
 	ha_ct1_t *g;
 	if (n == 0) return 0;
+	///corresponding small hash index
 	g = &h->h[a[0]&mask];
 	for (j = 0; j < n; ++j) {
 		int ins = 1, absent;
+		///x is a 64-bit word, h->pre=12
+		///all elements at a have the same low 12 bits
+		///so low 12 bits are not useful
 		uint64_t x = a[j] >> h->pre;
 		khint_t k;
 		if ((a[j]&mask) != (a[0]&mask)) continue;
 		if (create_new) {
+			///for 0-th counting, g->b = NULL
 			if (g->b)
 				ins = (yak_bf_insert(g->b, x) == h->n_hash);
+			///for 0-th counting, g->b = NULL
+			///x = the high 52 bits of a[j] + low 12 bits 0
+			///the low 12 bits are used for counting
 			if (ins) {
 				k = yak_ct_put(g->h, x << YAK_COUNTER_BITS | (g->b? 1 : 0), &absent);
 				if (absent) ++n_ins;
@@ -219,6 +236,8 @@ static void worker_ct_hist(void *data, long i, int tid) // callback for kt_for()
 			++cnt[kh_key(g, k)&YAK_MAX_COUNT];
 }
 
+///YAK_N_COUNTS is also 4096
+///used for calculating k-mer histogram
 static void ha_ct_hist(const ha_ct_t *h, int64_t cnt[YAK_N_COUNTS], int n_thread)
 {
 	hist_aux_t a;
@@ -226,6 +245,7 @@ static void ha_ct_hist(const ha_ct_t *h, int64_t cnt[YAK_N_COUNTS], int n_thread
 	a.h = h;
 	memset(cnt, 0, YAK_N_COUNTS * sizeof(uint64_t));
 	CALLOC(a.cnt, n_thread);
+	///start 4096 threads
 	kt_for(n_thread, worker_ct_hist, &a, 1<<h->pre);
 	for (i = 0; i < YAK_N_COUNTS; ++i) cnt[i] = 0;
 	for (j = 0; j < n_thread; ++j)
@@ -265,6 +285,7 @@ static void ha_ct_shrink(ha_ct_t *h, int min, int max, int n_thread)
 	int i;
 	shrink_aux_t a;
 	a.h = h, a.min = min, a.max = max;
+	///still start 4096 threads
 	kt_for(n_thread, worker_ct_shrink, &a, 1<<h->pre);
 	for (i = 0, h->tot = 0; i < 1<<h->pre; ++i)
 		h->tot += kh_size(h->h[i].h);
@@ -306,6 +327,7 @@ static void worker_pt_gen(void *data, long i, int tid) // callback for kt_for()
 			int absent;
 			khint_t l;
 			l = yak_pt_put(b->h, kh_key(g, k) >> a->ct->pre << YAK_COUNTER_BITS, &absent);
+			///this should be the start index of kh_key's corresponding pos at ha_idxpos_t* a
 			kh_val(b->h, l) = b->n;
 			b->n += kh_key(g, k) & YAK_MAX_COUNT;
 		}
@@ -412,8 +434,12 @@ typedef struct {
 	ha_mz1_t *b;
 } ch_buf_t;
 
+///p = 12
 static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k-mer $y to a linear buffer
 {
+	///assign k-mer to one of the 4096 bins
+	///using low 12 bits for assigning
+	///so all elements at b have the same low 12 bits
 	int pre = y & ((1<<p) - 1);
 	ch_buf_t *b = &buf[pre];
 	if (b->n == b->m) {
@@ -425,6 +451,7 @@ static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k
 
 static inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y)
 {
+	///assign minimizer to one of 4096 bins by low 12 bits
 	int pre = y->x & ((1<<p) - 1);
 	ch_buf_t *b = &buf[pre];
 	if (b->n == b->m) {
@@ -434,13 +461,17 @@ static inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y)
 	b->b[b->n++] = *y;
 }
 
+///buf is the read block, k is the k-mer length, p = 12, len is the read length, seq is the read
 static void count_seq_buf(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf
 {
 	int i, l;
 	uint64_t x[4], mask = (1ULL<<k) - 1, shift = k - 1;
 	for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) {
 		int c = seq_nt4_table[(uint8_t)seq[i]];
+		///c = 00, 01, 10, 11
 		if (c < 4) { // not an "N" base
+			///x[0] & x[1] are the forward k-mer
+			///x[2] & x[3] are the reverse complementary k-mer
 			x[0] = (x[0] << 1 | (c&1))  & mask;
 			x[1] = (x[1] << 1 | (c>>1)) & mask;
 			x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1))  << shift;
@@ -483,12 +514,13 @@ KSEQ_INIT(gzFile, gzread)
 #define HAF_RS_WRITE_SEQ 0x8
 #define HAF_RS_READ      0x10
 #define HAF_CREATE_NEW   0x20
+#define HAF_SKIP_READ    0x40
 
 typedef struct { // global data structure for kt_pipeline()
 	const yak_copt_t *opt;
 	const void *flt_tab;
 	int flag, create_new, is_store;
-	uint64_t n_seq;
+	uint64_t n_seq; ///number of total reads
 	kseq_t *ks;
 	UC_Read ucr;
 	ha_ct_t *ct;
@@ -499,7 +531,8 @@ typedef struct { // global data structure for kt_pipeline()
 
 typedef struct { // data structure for each step in kt_pipeline()
 	pl_data_t *p;
-	uint64_t n_seq0;
+	uint64_t n_seq0;  ///the start index of current buffer block at R_INF
+	///sum_len = total bases, nk = number of k-mers
 	int n_seq, m_seq, sum_len, nk;
 	int *len;
 	char **seq;
@@ -514,15 +547,17 @@ static void worker_for_insert(void *data, long i, int tid) // callback for kt_fo
 	ch_buf_t *b = &s->buf[i];
 	if (s->p->pt)
 		b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b);
-	else
+	else///for 0-th count, go into here
 		b->n_ins += ha_ct_insert_list(s->p->ct, s->p->create_new, b->n, b->a);
 }
 
 static void worker_for_mz(void *data, long i, int tid)
 {
 	st_data_t *s = (st_data_t*)data;
+	///get the corresponding minimzer vector of this read
 	ha_mz1_v *b = &s->mz_buf[tid];
 	s->mz_buf[tid].n = 0;
+	///s->p->opt->w = 51, s->p->opt->k
 	ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab);
 	s->mz[i].n = s->mz[i].m = b->n;
 	MALLOC(s->mz[i].a, b->n);
@@ -540,6 +575,11 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		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;
@@ -559,12 +599,15 @@ 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);
 				}
 				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);
@@ -572,23 +615,25 @@ 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);
 					}
 				}
+				///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, 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;
 			}
@@ -596,18 +641,24 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		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]);
@@ -618,18 +669,21 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		} 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) {
+			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 {
+			} 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);
@@ -640,6 +694,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 			}
 			free(s->mz);
 		}
+		///just clean seq
 		free(s->seq); free(s->len);
 		s->seq = 0, s->len = 0;
 		return s;
@@ -647,7 +702,9 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		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);
@@ -666,8 +723,51 @@ 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
 	int read_rs = (rs && (flag & HAF_RS_READ));
 	pl_data_t pl;
 	gzFile fp = 0;
@@ -676,23 +776,29 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 	if (read_rs) {
 		pl.rs_in = rs;
 		init_UC_Read(&pl.ucr);
-	} else {
+	} else {///for 0-th counting, go into here
 		if ((fp = gzopen(fn, "r")) == 0) return 0;
 		pl.ks = kseq_init(fp);
 	}
+	///for 0-th counting, read all reads into pl.rs_out
 	if (rs && (flag & (HAF_RS_WRITE_LEN|HAF_RS_WRITE_SEQ)))
 		pl.rs_out = rs;
+	///for 0-th counting, flt_tab = NULL
+	///for 1-th counting, flt_tab = NULL
 	pl.flt_tab = flt_tab;
 	pl.opt = opt;
 	pl.flag = flag;
-	if (p0) {
+	if (p0) {///for 1-th counting, p0 = NULL
 		pl.pt = p0, pl.create_new = 0; // never create new elements in a position table
 		assert(p0->k == opt->k && p0->pre == opt->pre);
 	} else if (c0) {
 		pl.ct = c0, pl.create_new = !!(flag&HAF_CREATE_NEW);
 		assert(c0->k == opt->k && c0->pre == opt->pre);
-	} else {
+	} 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 = 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);
 	}
 	kt_pipeline(3, worker_count, &pl, 3);
@@ -713,6 +819,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
 	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);
@@ -721,10 +828,16 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
 	}
 	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)
@@ -771,23 +884,63 @@ 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 *
  *************************/
 
-void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov)
+void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode)
 {
 	yak_ft_t *flt_tab;
 	int64_t cnt[YAK_N_COUNTS];
-	int peak_hom, peak_het, cutoff;
+	int peak_hom, peak_het, cutoff = YAK_MAX_COUNT - 1, ex_flag = 0;
+	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, NULL, NULL, rs);
+	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs);
-	ha_ct_hist(h, cnt, asm_opt->thread_num);
+	if((asm_opt->flag & HA_F_VERBOSE_GFA))
-	peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
+	{
-	if (hom_cov) *hom_cov = peak_hom;
+		write_ct_index((void*)h, asm_opt->output_file_name);
-	if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
+		// load_ct_index(&ha_ct_table, asm_opt->output_file_name);
-	cutoff = (int)(peak_hom * asm_opt->high_factor);
+		// debug_ct_index((void*)h, ha_ct_table);
-	if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
+		// 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, 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
+		cutoff = (int)(peak_hom * asm_opt->high_factor);
+		if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
+	}
 	ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num);
 	flt_tab = gen_hh(h);
 	ha_ct_destroy(h);
@@ -796,36 +949,42 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov)
 	return (void*)flt_tab;
 }
 
-ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov, int *het_cov)
+ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov)
 {
 	int64_t cnt[YAK_N_COUNTS], tot_cnt;
 	int peak_hom, peak_het, i, extra_flag1, extra_flag2;
 	ha_ct_t *ct;
 	ha_pt_t *pt;
-	if (read_from_store) {
+	if (read_from_store) {///if reads have already been read
 		extra_flag1 = extra_flag2 = HAF_RS_READ;
-	} else if (rs->total_reads == 0) {
+	} else if (rs->total_reads == 0) {///if reads & length have not been scanned
 		extra_flag1 = HAF_RS_WRITE_LEN;
 		extra_flag2 = HAF_RS_WRITE_SEQ;
-	} else {
+	} else {///if length has been loaded but reads have not
 		extra_flag1 = HAF_RS_WRITE_SEQ;
 		extra_flag2 = HAF_RS_READ;
 	}
+	if(is_hp_mode) extra_flag1 |= HAF_SKIP_READ, extra_flag2 |= HAF_SKIP_READ;
+
 	ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs);
 	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
 			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);
+	///here ha_ct_shrink is mostly used to remove k-mer appearing only 1 time
 	if (flt_tab == 0) {
 		int cutoff = (int)(peak_hom * asm_opt->high_factor);
 		if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
+		if((extra_flag1 & HAF_SKIP_READ) && (extra_flag2 & HAF_SKIP_READ)) cutoff = YAK_MAX_COUNT - 1;
 		ha_ct_shrink(ct, 2, cutoff, asm_opt->thread_num);
 		for (i = 2, tot_cnt = 0; i <= cutoff; ++i) tot_cnt += cnt[i] * i;
 	} else {
+		///Note: here is just to remove minimizer appearing YAK_MAX_COUNT times
+		///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen
 		ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num);
 		for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i;
 	}
@@ -837,3 +996,239 @@ 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)pt->tot_pos);
 	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_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);
+    FILE* fp = fopen(gfa_name, "w");
+	if (!fp) {
+		free(gfa_name);
+        return 0;
+    }
+	yak_ft_t *ha_flt_tab = (yak_ft_t*)flt_tab;
+
+	if(ha_flt_tab)
+	{
+		fwrite("f", 1, 1, fp);
+		yak_ft_save(ha_flt_tab, fp);
+	}
+
+
+	if(ha_idx)
+	{
+		int i;
+		ha_pt1_t *g;
+		fwrite("h", 1, 1, fp);
+		fwrite(&ha_idx->k, sizeof(ha_idx->k), 1, fp);
+		fwrite(&ha_idx->pre, sizeof(ha_idx->pre), 1, fp);
+		fwrite(&ha_idx->tot, sizeof(ha_idx->tot), 1, fp);
+		fwrite(&ha_idx->tot_pos, sizeof(ha_idx->tot_pos), 1, fp);
+
+		for (i = 0; i < 1<<ha_idx->pre; ++i) 
+		{
+			g = &(ha_idx->h[i]);
+			yak_pt_save(g->h, fp);
+			fwrite(&g->n, sizeof(g->n), 1, fp);
+			fwrite(g->a, sizeof(ha_idxpos_t), g->n, fp);
+		}
+	}
+
+	fwrite(&opt->number_of_round, sizeof(opt->number_of_round), 1, fp);
+	fwrite(&opt->hom_cov, sizeof(opt->hom_cov), 1, fp);
+	fwrite(&opt->het_cov, sizeof(opt->het_cov), 1, fp);
+	fwrite(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp);
+
+
+	write_All_reads(r, gfa_name);
+
+	fprintf(stderr, "[M::%s] Index has been written.\n", __func__);
+    free(gfa_name);
+	fclose(fp);
+	return 1;
+}
+
+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);
+    FILE* fp = fopen(gfa_name, "r");
+	if (!fp) {
+		free(gfa_name);
+        return 0;
+    }
+
+	ha_pt_t *ha_idx = NULL;
+	char mode = 0;
+	int f_flag, absent, i;
+	double index_time, index_s_time, pos_time, pos_s_time;
+
+	
+	
+	f_flag += fread(&mode, 1, 1, fp);
+	if(mode == 'f')
+	{
+		index_time = yak_realtime();
+
+		yak_ft_load((yak_ft_t **)r_flt_tab, fp);
+		
+	    f_flag += fread(&mode, 1, 1, fp);
+
+		fprintf(stderr, "[M::%s::%.3f] ==> Loaded flt table\n", __func__, yak_realtime()-index_time);
+	}
+	///insert using multiple threads???
+	if(mode == 'h')
+	{
+		pos_time = index_time = 0;
+
+		CALLOC(ha_idx, 1);
+		ha_pt1_t *g;
+		f_flag += fread(&ha_idx->k, sizeof(ha_idx->k), 1, fp);
+		f_flag += fread(&ha_idx->pre, sizeof(ha_idx->pre), 1, fp);
+		f_flag += fread(&ha_idx->tot, sizeof(ha_idx->tot), 1, fp);
+		f_flag += fread(&ha_idx->tot_pos, sizeof(ha_idx->tot_pos), 1, fp);
+		CALLOC(ha_idx->h, 1<<ha_idx->pre);
+		for (i = 0; i < 1<<ha_idx->pre; ++i) 
+		{
+			index_s_time = yak_realtime();
+
+			g = &(ha_idx->h[i]);
+			yak_pt_load(&(g->h), fp);
+
+			index_time += yak_realtime() - index_s_time;
+			
+			pos_s_time = yak_realtime();
+
+			f_flag += fread(&g->n, sizeof(g->n), 1, fp);
+			MALLOC(g->a, g->n);
+			f_flag += fread(g->a, sizeof(ha_idxpos_t), g->n, fp);
+
+			pos_time += yak_realtime() - pos_s_time;
+		}
+		(*r_ha_idx) = ha_idx;
+
+		fprintf(stderr, "[M::%s::%.3f(index)/%.3f(pos)] ==> Loaded pos table\n", __func__, index_time, pos_time);
+	}
+
+	if(mode != 'h' && mode != 'f')
+	{
+		free(gfa_name);
+		fclose(fp);
+		return 0;
+	}
+
+
+	f_flag += fread(&absent, sizeof(absent), 1, fp);
+	if(absent != opt->number_of_round)
+	{
+		fprintf(stderr, "ERROR: different number of rounds!\n");
+		exit(1);
+	}
+
+	f_flag += fread(&opt->hom_cov, sizeof(opt->hom_cov), 1, fp);
+	f_flag += fread(&opt->het_cov, sizeof(opt->het_cov), 1, fp);
+	f_flag += fread(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp);
+
+
+	fclose(fp);
+	
+	if(!load_All_reads(r, gfa_name))
+	{
+		free(gfa_name);
+		return 0;
+	}
+
+
+	memset(r->trio_flag, AMBIGU, r->total_reads*sizeof(uint8_t));
+	r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
+    r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
+	for (i = 0; i < (long long)r->total_reads; i++)
+    {
+        init_ma_hit_t_alloc(&(r->paf[i]));
+        init_ma_hit_t_alloc(&(r->reverse_paf[i]));
+    }
+
+	fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__);
+
+	free(gfa_name);
+	return 1;
+}

htab.h

@@ -6,7 +6,9 @@
 #include "CommandLines.h"
 
 typedef struct {
-	uint64_t x;
+	uint64_t x; ///x is the hash key
+	///rid is the read id, pos is the end pos of this minimizer, rev is the direction
+	///span is the length of this k-mer. For non-HPC k-mer, span may not be equal to k
 	uint64_t rid:28, pos:27, rev:1, span:8;
 } ha_mz1_t;
 
@@ -25,15 +27,25 @@ typedef struct ha_abuf_s ha_abuf_t;
 extern const unsigned char seq_nt4_table[256];
 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);
+
+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);
 void ha_ft_destroy(void *h);
 
-ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov, int *het_cov);
+ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov);
 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_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);
 uint64_t ha_abuf_mem(const ha_abuf_t *ab);
@@ -48,7 +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);
-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
 {
@@ -76,6 +90,7 @@ static inline uint64_t yak_hash64_64(uint64_t key)
 
 static inline uint64_t yak_hash_long(uint64_t x[4])
 {
+	///compare forward k-mer and reverse complementary strand
 	int j = x[1] < x[3]? 0 : 1;
 	return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]);
 }

khashl.h

@@ -133,6 +133,30 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 			h->count = 0; \
 		} \
 	}
+#define __KHASHL_IMPL_S_L(SCOPE, HType, prefix, khkey_t) \
+	SCOPE khint_t prefix##_save(HType *h, FILE* fp) { \
+		if (!h) return 0; \
+		khint_t n_buckets = (h->keys? 1U<<h->bits : 0U); \
+		fwrite(&n_buckets, sizeof(n_buckets), 1, fp); \
+		fwrite(&h->bits, sizeof(h->bits), 1, fp); \
+		fwrite(&h->count, sizeof(h->count), 1, fp); \
+		fwrite(h->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \
+		fwrite(h->keys, sizeof(khkey_t), n_buckets, fp); \
+		return 1; \
+	} \
+	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); \
+		(*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); \
+		return 1; \
+	} \
 
 #define __KHASHL_IMPL_GET(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
 	SCOPE khint_t prefix##_getp(const HType *h, const khkey_t *key) { \
@@ -245,6 +269,7 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 #define KHASHL_INIT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
 	__KHASHL_TYPE(HType, khkey_t) \
 	__KHASHL_IMPL_BASIC(SCOPE, HType, prefix) \
+	__KHASHL_IMPL_S_L(SCOPE, HType, prefix, khkey_t) \
 	__KHASHL_IMPL_GET(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
 	__KHASHL_IMPL_RESIZE(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
 	__KHASHL_IMPL_PUT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \
@@ -264,6 +289,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 	KHASHL_INIT(KH_LOCAL, HType, prefix##_s, HType##_s_bucket_t, prefix##_s_hash, prefix##_s_eq) \
 	SCOPE HType *prefix##_init(void) { return prefix##_s_init(); } \
 	SCOPE void prefix##_destroy(HType *h) { prefix##_s_destroy(h); } \
+	SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_s_save(h, fp); } \
+	SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_s_load(h, fp); } \
 	SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_s_resize(h, new_n_buckets); } \
 	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_s_bucket_t t; t.key = key; return prefix##_s_getp(h, &t); } \
 	SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_s_del(h, k); } \
@@ -276,6 +303,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 	KHASHL_INIT(KH_LOCAL, HType, prefix##_m, HType##_m_bucket_t, prefix##_m_hash, prefix##_m_eq) \
 	SCOPE HType *prefix##_init(void) { return prefix##_m_init(); } \
 	SCOPE void prefix##_destroy(HType *h) { prefix##_m_destroy(h); } \
+	SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_m_save(h, fp); } \
+	SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_m_load(h, fp); } \
 	SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_m_resize(h, new_n_buckets); } \
 	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_m_bucket_t t; t.key = key; return prefix##_m_getp(h, &t); } \
 	SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_m_del(h, k); } \
@@ -287,6 +316,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 	KHASHL_INIT(KH_LOCAL, HType, prefix##_cs, HType##_cs_bucket_t, __kh_cached_hash, prefix##_cs_eq) \
 	SCOPE HType *prefix##_init(void) { return prefix##_cs_init(); } \
 	SCOPE void prefix##_destroy(HType *h) { prefix##_cs_destroy(h); } \
+	SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_cs_save(h, fp); } \
+	SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_cs_load(h, fp); } \
 	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_cs_bucket_t t; t.key = key; t.hash = __hash_fn(key); return prefix##_cs_getp(h, &t); } \
 	SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_cs_del(h, k); } \
 	SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_cs_bucket_t t; t.key = key, t.hash = __hash_fn(key); return prefix##_cs_putp(h, &t, absent); }
@@ -297,6 +328,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
 	KHASHL_INIT(KH_LOCAL, HType, prefix##_cm, HType##_cm_bucket_t, __kh_cached_hash, prefix##_cm_eq) \
 	SCOPE HType *prefix##_init(void) { return prefix##_cm_init(); } \
 	SCOPE void prefix##_destroy(HType *h) { prefix##_cm_destroy(h); } \
+	SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_cm_save(h, fp); } \
+	SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_cm_load(h, fp); } \
 	SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_cm_bucket_t t; t.key = key; t.hash = __hash_fn(key); return prefix##_cm_getp(h, &t); } \
 	SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_cm_del(h, k); } \
 	SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_cm_bucket_t t; t.key = key, t.hash = __hash_fn(key); return prefix##_cm_putp(h, &t, absent); }

sketch.cpp

@@ -37,7 +37,11 @@ static inline int tq_shift(tiny_queue_t *q)
  * @param p      minimizers
  */
 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)
-{
+{	///in default, w = 51, k = 51, is_hpc = 1
+	/**
+	 uint64_t x;
+	 uint64_t rid:28, pos:27, rev:1, span:8;
+	 **/
 	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};
 	int i, j, l, buf_pos, min_pos, kmer_span = 0;
@@ -45,8 +49,10 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 	tiny_queue_t tq;
 
 	assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
+	///sizeof(ha_mz1_t) = 16
 	memset(buf, 0xff, w * 16);
 	memset(&tq, 0, sizeof(tiny_queue_t));
+	///len/w is the evaluated minimizer numbers
 	kv_resize(ha_mz1_t, *p, p->n + len/w);
 
 	for (i = l = buf_pos = min_pos = 0; i < len; ++i) {
@@ -65,7 +71,11 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 				tq_push(&tq, skip_len);
 				kmer_span += skip_len;
 				if (tq.count > k) kmer_span -= tq_shift(&tq);
-			} else kmer_span = l + 1 < k? l + 1 : k;
+			} 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?
+			///kmer_span is used to calculate anchor pos on reverse complementary strand
+
 			kmer[0] = (kmer[0] << 1 | (c&1))  & mask;                  // forward k-mer
 			kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;
 			kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1))  << shift1; // reverse k-mer
@@ -80,6 +90,12 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 					info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;
 			}
 		} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
+
+
+		//for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i
+		//i is the real base iterator, while l is the HPC base iterator
+		//only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX)
+		//but even if l < k, infor is still stored into buf
 		buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
 		if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
 			for (j = buf_pos + 1; j < w; ++j)
@@ -87,15 +103,26 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 			for (j = 0; j < buf_pos; ++j)
 				if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
 		}
+		/**
+		 * There are three cases:
+		 * 1. info.x <= min.x, means info is a new minimizer
+		 * 2. info.x > min.x, info is not a new minimizer
+		 *    (1) buf_pos != min_pos, do nothing
+		 *    (2) buf_pos == min_pos, means current minimizer has moved outside the window
+		 * **/
+		///three cases: 1. 
 		if (info.x <= min.x) { // a new minimum; then write the old min
 			if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
 			min = info, min_pos = buf_pos;
 		} else if (buf_pos == min_pos) { // old min has moved outside the window
 			if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
+			///buf_pos == min_pos, means current minimizer has moved outside the window
+			///so for now we need to find a new minimizer at the current window (w k-mers)
 			for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers
 				if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer
 			for (j = 0; j <= buf_pos; ++j)
 				if (min.x >= buf[j].x) min = buf[j], min_pos = j;
+
 			if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
 				for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
 					if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
@@ -108,3 +135,127 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 	if (min.x != UINT64_MAX)
 		kv_push(ha_mz1_t, *p, min);
 }
+
+
+
+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)
+{	///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;
+	ha_mz1_t buf[256], min = dummy;
+	tiny_queue_t tq;
+	if(k_flag != NULL)
+	{
+		kv_resize(uint8_t, k_flag->a, (uint64_t)len);
+		k_flag->a.n = len; 
+		memset(k_flag->a.a, 0, k_flag->a.n);
+	} 
+
+
+	assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
+	///sizeof(ha_mz1_t) = 16
+	memset(buf, 0xff, w * 16);
+	memset(&tq, 0, sizeof(tiny_queue_t));
+	///len/w is the evaluated minimizer numbers
+	kv_resize(ha_mz1_t, *p, p->n + len/w);
+
+	for (i = l = buf_pos = min_pos = 0; i < len; ++i) {
+		int c = seq_nt4_table[(uint8_t)str[i]];
+		ha_mz1_t info = dummy;
+		if (c < 4) { // not an ambiguous base
+			int z;
+			if (is_hpc) {
+				int skip_len = 1;
+				if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) {
+					for (skip_len = 2; i + skip_len < len; ++skip_len)
+						if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c)
+							break;
+					i += skip_len - 1; // put $i at the end of the current homopolymer run
+				}
+				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
+			///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
+
+			kmer[0] = (kmer[0] << 1 | (c&1))  & mask;                  // forward k-mer
+			kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;
+			kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1))  << shift1; // reverse k-mer
+			kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1;
+			if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand
+			z = kmer[1] < kmer[3]? 0 : 1; // strand
+			++l;
+			if (l >= k && kmer_span < 256) {
+				uint64_t y;
+				y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
+
+				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;
+				if(k_flag != NULL) k_flag->a.a[i]++;
+				if(k_flag != NULL && filtered == 1) k_flag->a.a[i]++;
+			}
+		} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
+
+
+		//for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i
+		//i is the real base iterator, while l is the HPC base iterator
+		//only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX)
+		//but even if l < k, infor is still stored into buf
+		buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
+		if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
+			for (j = buf_pos + 1; j < w; ++j)
+				if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
+			for (j = 0; j < buf_pos; ++j)
+				if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
+		}
+		/**
+		 * There are three cases:
+		 * 1. info.x <= min.x, means info is a new minimizer
+		 * 2. info.x > min.x, info is not a new minimizer
+		 *    (1) buf_pos != min_pos, do nothing
+		 *    (2) buf_pos == min_pos, means current minimizer has moved outside the window
+		 * **/
+		///three cases: 1. 
+		if (info.x <= min.x) { // a new minimum; then write the old min
+			if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
+			min = info, min_pos = buf_pos;
+		} else if (buf_pos == min_pos) { // old min has moved outside the window
+			if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
+			///buf_pos == min_pos, means current minimizer has moved outside the window
+			///so for now we need to find a new minimizer at the current window (w k-mers)
+			for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers
+				if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer
+			for (j = 0; j <= buf_pos; ++j)
+				if (min.x >= buf[j].x) min = buf[j], min_pos = j;
+
+			if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
+				for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
+					if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
+				for (j = 0; j <= buf_pos; ++j)
+					if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
+			}
+		}
+		if (++buf_pos == w) buf_pos = 0;
+	}
+	if (min.x != UINT64_MAX)
+		kv_push(ha_mz1_t, *p, min);
+}