update README

Probably has no practical effect

Assembly.cpp

@@ -8,13 +8,12 @@
 #include "Hash_Table.h"
 #include "POA.h"
 #include "Correct.h"
-#include "Output.h"
 #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);
 
 All_reads R_INF;
-pthread_mutex_t statistics;
 
 void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length, char* new_read, int* new_length)
 {
@@ -251,7 +250,12 @@ inline void push_cigar(Compressed_Cigar_record* records, long long ID, Cigar_rec
 void push_overlaps(ma_hit_t_alloc* paf, overlap_region_alloc* overlap_list, int flag, All_reads* R_INF, int if_reverse)
 {
     long long i = 0, xLen, yLen;
+	int32_t size = 0;
     ma_hit_t tmp;
+	for (i = 0; i < (long long)overlap_list->length; ++i)
+		if (overlap_list->list[i].is_match == flag)
+			++size;
+	resize_ma_hit_t_alloc(paf, size);
     clear_ma_hit_t_alloc(paf);
     for (i = 0; i < (long long)overlap_list->length; i++)
     {
@@ -279,7 +283,6 @@ void push_overlaps(ma_hit_t_alloc* paf, overlap_region_alloc* overlap_list, int
                 tmp.te = overlap_list->list[i].y_pos_e;
             }
 
-
             ///for overlap_list, the x_strand of all overlaps are 0, so the tmp.rev is the same as the y_strand
             tmp.rev = overlap_list->list[i].y_pos_strand;
 
@@ -323,7 +326,7 @@ long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_l
     long long i = 0;
     long long available_overlaps = 0;
     ma_hit_t tmp;
-    clear_ma_hit_t_alloc(paf);
+    clear_ma_hit_t_alloc(paf); // paf has been preallocated, so we don't need preallocation
     for (i = 0; i < (long long)overlap_list->length; i++)
     {
         if (overlap_list->list[i].is_match == flag)
@@ -375,13 +378,14 @@ long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_l
 }
 
 typedef struct {
-	int is_final;
+	int is_final, save_ov;
 	// chaining and overlapping related buffers
 	UC_Read self_read, ovlp_read;
 	Candidates_list clist;
 	overlap_region_alloc olist;
 	ha_abuf_t *ab;
 	// error correction related buffers
+	int64_t num_read_base, num_correct_base, num_recorrect_base;
 	Cigar_record cigar1;
 	Graph POA_Graph;
 	Graph DAGCon;
@@ -390,11 +394,11 @@ typedef struct {
 	Round2_alignment round2;
 } ha_ovec_buf_t;
 
-ha_ovec_buf_t *ha_ovec_init(int is_final)
+ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov)
 {
 	ha_ovec_buf_t *b;
 	CALLOC(b, 1);
-	b->is_final = !!is_final;
+	b->is_final = !!is_final, b->save_ov = !!save_ov;
 	init_UC_Read(&b->self_read);
 	init_UC_Read(&b->ovlp_read);
 	init_Candidates_list(&b->clist);
@@ -466,204 +470,141 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b)
 	return mem;
 }
 
-void* Overlap_calculate_heap_merge(void* arg)
+static void worker_ovec(void *data, long i, int tid)
 {
-	long long num_read_base = 0;
-	long long num_correct_base = 0;
-	long long num_recorrect_base = 0;
-	long long mem_buf;
+	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
 	int fully_cov, abnormal;
 
-	int thr_ID = *((int*)arg);
-	long long i = 0;
-	ha_ovec_buf_t *b;
+	ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1);
 
-	b = ha_ovec_init(0);
-	for (i = thr_ID; i < (long long)R_INF.total_reads; i = i + asm_opt.thread_num) {
-		//get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.02, 1);
+	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;
+
+	if (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);
+	}
+}
+
+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);
 
-		clear_Cigar_record(&b->cigar1);
-		clear_Round2_alignment(&b->round2);
+		fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i));
+		recover_UC_Read(&b->self_read, &R_INF, i);
+		fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq);
 
-		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);
-
-		num_read_base += b->self_read.length;
-		num_correct_base += b->correct.corrected_base;
-		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;
-
-		push_overlaps(&(R_INF.paf[i]), &b->olist, 1, &R_INF, asm_opt.roundID%2);
-		push_overlaps(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, asm_opt.roundID%2);
-	}
-	finish_output_buffer();
-	mem_buf = ha_ovec_mem(b);
-	ha_ovec_destroy(b);
-
-	pthread_mutex_lock(&statistics);
-	asm_opt.num_bases += num_read_base;
-	asm_opt.num_corrected_bases += num_correct_base;
-	asm_opt.num_recorrected_bases += num_recorrect_base;
-	asm_opt.mem_buf += mem_buf;
-	pthread_mutex_unlock(&statistics);
-	return NULL;
-}
-
-void* Output_related_reads(void* arg)
-{
-	int thr_ID = *((int*)arg);
-	long long i = 0;
-	ha_ovec_buf_t *b;
-
-	long long required_read_name_length = strlen(asm_opt.required_read_name);
-	b = ha_ovec_init(0);
-	for (i = thr_ID; i < (long long)R_INF.total_reads; i = i + asm_opt.thread_num) {
-		if (required_read_name_length == (long long)Get_NAME_LENGTH((R_INF),i)
-				&&
-				memcmp(asm_opt.required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
-		{
-			//get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.02, 1);
-			ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1);
-
-			fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i));
-			recover_UC_Read(&b->self_read, &R_INF, i);
+		for (k = 0; k < b->olist.length; k++) {
+			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);
 			fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq);
-
-			uint64_t k;
-			for (k = 0; k < b->olist.length; k++) {
-				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);
-				fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq);
-			}
 		}
 	}
-	finish_output_buffer();
-	ha_ovec_destroy(b);
-	return NULL;
 }
 
-inline long long get_N_occ(char* seq, long long length)
+static inline long long get_N_occ(char* seq, long long length)
 {
-    long long N_occ = 0;
-    long long j;
-    for (j = 0; j < length; j++)
-    {
-        if(seq_nt6_table[(uint8_t)seq[j]] >= 4)
-        {
-            N_occ++;
-        }
-    }
-    return N_occ;
+	long long j, N_occ = 0;
+	for (j = 0; j < length; j++)
+		if(seq_nt6_table[(uint8_t)seq[j]] >= 4)
+			N_occ++;
+	return N_occ;
 }
 
+typedef struct {
+	UC_Read g_read;
+	int first_round_read_size;
+	int second_round_read_size;
+	char *first_round_read;
+	char *second_round_read;
+} ha_ecsave_buf_t;
 
-void* Save_corrected_reads(void* arg)
+static void worker_ec_save(void *data, long i, int tid)
 {
-    int thr_ID = *((int*)arg);
-    long long i;
-    UC_Read g_read;
-    init_UC_Read(&g_read);
+	ha_ecsave_buf_t *e = (ha_ecsave_buf_t*)data + tid;
 
-    int first_round_read_size = 10000;
-    char* first_round_read = (char*)malloc(first_round_read_size);
+	Cigar_record cigar;
+	int first_round_read_length;
+	int second_round_read_length;
+	uint64_t N_occ;
 
-    int second_round_read_size = 10000;
-    char* second_round_read = (char*)malloc(second_round_read_size);
+	char *new_read;
+	int new_read_length;
 
-    Cigar_record cigar;
-    int first_round_read_length;
-    int second_round_read_length;
-    uint64_t N_occ;
+	recover_UC_Read(&e->g_read, &R_INF, i);
 
-    char* new_read;
-    int new_read_length;
+	// round 1
+	if ((long long)R_INF.cigars[i].new_length > e->first_round_read_size) {
+		e->first_round_read_size = R_INF.cigars[i].new_length;
+		REALLOC(e->first_round_read, e->first_round_read_size);
+	}
 
-    for (i = thr_ID; i < (long long)R_INF.total_reads; i = i + asm_opt.thread_num)
-    {
-        recover_UC_Read(&g_read, &R_INF, i);
+	cigar.length = R_INF.cigars[i].length;
+	cigar.lost_base_length = R_INF.cigars[i].lost_base_length;
+	cigar.record = R_INF.cigars[i].record;
+	cigar.lost_base = R_INF.cigars[i].lost_base;
 
-        /********************************1 round******************************/
-        if((long long)R_INF.cigars[i].new_length > first_round_read_size)
-        {
-            first_round_read_size = R_INF.cigars[i].new_length;
-            first_round_read = (char*)realloc(first_round_read, first_round_read_size);
-        }
+	get_corrected_read_from_cigar(&cigar, e->g_read.seq, e->g_read.length, e->first_round_read, &first_round_read_length);
 
-        cigar.length = R_INF.cigars[i].length;
-        cigar.lost_base_length = R_INF.cigars[i].lost_base_length;
-        cigar.record = R_INF.cigars[i].record;
-        cigar.lost_base = R_INF.cigars[i].lost_base;
+	// round 2
+	if ((long long)R_INF.second_round_cigar[i].new_length > e->second_round_read_size) {
+		e->second_round_read_size = R_INF.second_round_cigar[i].new_length;
+		REALLOC(e->second_round_read, e->second_round_read_size);
+	}
+	cigar.length = R_INF.second_round_cigar[i].length;
+	cigar.lost_base_length = R_INF.second_round_cigar[i].lost_base_length;
+	cigar.record = R_INF.second_round_cigar[i].record;
+	cigar.lost_base = R_INF.second_round_cigar[i].lost_base;
 
-        get_corrected_read_from_cigar(&cigar, g_read.seq, g_read.length, first_round_read, &first_round_read_length);
+	get_corrected_read_from_cigar(&cigar, e->first_round_read, first_round_read_length, e->second_round_read, &second_round_read_length);
 
-        /********************************1 round******************************/
+	new_read = e->second_round_read;
+	new_read_length = second_round_read_length;
 
-        /********************************2 round******************************/
-        if((long long)R_INF.second_round_cigar[i].new_length > second_round_read_size)
-        {
-            second_round_read_size = R_INF.second_round_cigar[i].new_length;
-            second_round_read = (char*)realloc(second_round_read, second_round_read_size);
-        }
-        cigar.length = R_INF.second_round_cigar[i].length;
-        cigar.lost_base_length = R_INF.second_round_cigar[i].lost_base_length;
-        cigar.record = R_INF.second_round_cigar[i].record;
-        cigar.lost_base = R_INF.second_round_cigar[i].lost_base;
-        get_corrected_read_from_cigar(&cigar, first_round_read, first_round_read_length, 
-        second_round_read, &second_round_read_length);
+	if (asm_opt.roundID != asm_opt.number_of_round - 1)
+	{
+		///need modification
+		reverse_complement(new_read, new_read_length);
+	}
+	else if(asm_opt.number_of_round % 2 == 0)
+	{
+		///need modification
+		reverse_complement(new_read, new_read_length);
+	}
 
-        /********************************2 round******************************/
+	N_occ = get_N_occ(new_read, new_read_length);
 
-
-
-
-        new_read = second_round_read;
-        new_read_length = second_round_read_length;
-
-
-        if (asm_opt.roundID != asm_opt.number_of_round - 1)
-        {
-            ///need modification
-            reverse_complement(new_read, new_read_length);
-        }
-        else if(asm_opt.number_of_round % 2 == 0)
-        {
-            ///need modification
-            reverse_complement(new_read, new_read_length);
-        }
-
-
-        N_occ = get_N_occ(new_read, new_read_length);
-
-
-        if((long long)R_INF.read_size[i] < new_read_length)
-        {
-            R_INF.read_size[i] = new_read_length;
-            R_INF.read_sperate[i] = (uint8_t*)realloc(R_INF.read_sperate[i], R_INF.read_size[i]/4+1);
-        }
-
-        R_INF.read_length[i] = new_read_length;
-
-
-        ha_compress_base(Get_READ(R_INF, i),
-            new_read, new_read_length, 
-            &R_INF.N_site[i], N_occ);
-    }
-
-    destory_UC_Read(&g_read);
-    free(first_round_read);
-    free(second_round_read);
-
-    return NULL;
+	if ((long long)R_INF.read_size[i] < new_read_length) {
+		R_INF.read_size[i] = new_read_length;
+		REALLOC(R_INF.read_sperate[i], R_INF.read_size[i]/4+1);
+	}
+	R_INF.read_length[i] = new_read_length;
+	ha_compress_base(Get_READ(R_INF, i), new_read, new_read_length, &R_INF.N_site[i], N_occ);
 }
 
 void Output_corrected_reads()
@@ -691,33 +632,51 @@ void Output_corrected_reads()
 
 void ha_overlap_and_correct(int round)
 {
-	int i, *args;
-	pthread_t *_r_threads;
-	MALLOC(_r_threads, asm_opt.thread_num);
-	args = (int*)alloca(sizeof(int) * asm_opt.thread_num);
+	int i, hom_cov;
+	ha_ovec_buf_t **b;
+	ha_ecsave_buf_t *e;
 
-	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF); // build the index
-	for (i = 0; i < asm_opt.thread_num; i++) {
-		args[i] = i;
-		if (!asm_opt.required_read_name)
-			pthread_create(_r_threads + i, NULL, Overlap_calculate_heap_merge, (void*)&args[i]);
-		else
-			pthread_create(_r_threads + i, NULL, Output_related_reads, (void*)&args[i]);
-	}
-	for (i = 0; i < asm_opt.thread_num; i++)
-		pthread_join(_r_threads[i], NULL);
+	// 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); // build the index
+	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);
 	ha_pt_destroy(ha_idx);
 	ha_idx = 0;
 
+	// collect statistics
+	for (i = 0; i < asm_opt.thread_num; ++i) {
+		asm_opt.num_bases += b[i]->num_read_base;
+		asm_opt.num_corrected_bases += b[i]->num_correct_base;
+		asm_opt.num_recorrected_bases += b[i]->num_recorrect_base;
+		asm_opt.mem_buf += ha_ovec_mem(b[i]);
+		ha_ovec_destroy(b[i]);
+	}
+	free(b);
+
 	if (asm_opt.required_read_name) exit(0); // for debugging only
 
-	for (i = 0; i < asm_opt.thread_num; i++) {
-		args[i] = i;
-		pthread_create(_r_threads + i, NULL, Save_corrected_reads, (void*)&args[i]);
+	// save corrected reads to R_INF
+	CALLOC(e, asm_opt.thread_num);
+	for (i = 0; i < asm_opt.thread_num; ++i) {
+		init_UC_Read(&e[i].g_read);
+		e[i].first_round_read_size = e[i].second_round_read_size = 50000;
+		CALLOC(e[i].first_round_read, e[i].first_round_read_size);
+		CALLOC(e[i].second_round_read, e[i].second_round_read_size);
 	}
-	for (i = 0; i < asm_opt.thread_num; i++)
-		pthread_join(_r_threads[i], NULL);
-	free(_r_threads);
+	kt_for(asm_opt.thread_num, worker_ec_save, e, R_INF.total_reads);
+	for (i = 0; i < asm_opt.thread_num; ++i) {
+		destory_UC_Read(&e[i].g_read);
+		free(e[i].first_round_read);
+		free(e[i].second_round_read);
+	}
+	free(e);
 }
 
 void update_overlaps(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf, 
@@ -988,47 +947,36 @@ UC_Read* g_read, UC_Read* overlap_read, uint8_t* c2n)
 
 }
 
-void* Final_overlap_calculate_heap_merge(void* arg)
+static void worker_ov_final(void *data, long i, int tid)
 {
-    int thr_ID = *((int*)arg);
-    uint64_t i = 0;
-	ha_ovec_buf_t *b;
-
-    uint8_t c2n[256];
+	ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid];
+    uint8_t c2n[256]; // this may be moved to ha_ovec_buf_t, but it should be fast to populate anyway
     memset(c2n, 4, 256);
     c2n[(uint8_t)'A'] = c2n[(uint8_t)'a'] = 0; c2n[(uint8_t)'C'] = c2n[(uint8_t)'c'] = 1;
 	c2n[(uint8_t)'G'] = c2n[(uint8_t)'g'] = 2; c2n[(uint8_t)'T'] = c2n[(uint8_t)'t'] = 3; // build the encoding table
 
-	b = ha_ovec_init(1);
-    for (i = thr_ID; i < R_INF.total_reads; i = i + asm_opt.thread_num)
-    {
-        //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);
+	//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);
 
-        /**
-        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);
-        **/
+	/**
+	  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);
-        ma_hit_sort_tn(R_INF.reverse_paf[i].buffer, R_INF.reverse_paf[i].length);
+	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);
+	ma_hit_sort_tn(R_INF.reverse_paf[i].buffer, R_INF.reverse_paf[i].length);
 
-        update_overlaps(&b->olist, &(R_INF.paf[i]), &b->self_read, &b->ovlp_read, 1, 1);
-        update_overlaps(&b->olist, &(R_INF.reverse_paf[i]), &b->self_read, &b->ovlp_read, 2, 0);
-        ///recover missing exact overlaps 
-        update_exact_overlaps(&b->olist, &b->self_read, &b->ovlp_read);
+	update_overlaps(&b->olist, &(R_INF.paf[i]), &b->self_read, &b->ovlp_read, 1, 1);
+	update_overlaps(&b->olist, &(R_INF.reverse_paf[i]), &b->self_read, &b->ovlp_read, 2, 0);
+	///recover missing exact overlaps
+	update_exact_overlaps(&b->olist, &b->self_read, &b->ovlp_read);
 
-        ///Final_phasing(&overlap_list, &cigarline, &g_read, &overlap_read, c2n);
-
-        push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, &b->olist, 1);
-        push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
-    }
-    finish_output_buffer();
-	ha_ovec_destroy(b);
-    return NULL;
+	///Final_phasing(&overlap_list, &cigarline, &g_read, &overlap_read, c2n);
+	push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, &b->olist, 1);
+	push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2);
 }
 
 void Output_PAF()
@@ -1195,40 +1143,40 @@ long long readNum, long long rescue_threshold, float cluster_threshold)
 
 void ha_overlap_final(void)
 {
-	int i, *args;
-	pthread_t *_r_threads;
-
-	MALLOC(_r_threads, asm_opt.thread_num);
-	args = (int*)alloca(sizeof(int) * asm_opt.thread_num);
-
-	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF);
-	for (i = 0; i < asm_opt.thread_num; i++) {
-		args[i] = i;
-		pthread_create(_r_threads + i, NULL, Final_overlap_calculate_heap_merge, (void*)&args[i]);
-	}
-	for (i = 0; i < asm_opt.thread_num; i++)
-		pthread_join(_r_threads[i], NULL);
-	free(_r_threads);
+	int i, hom_cov;
+	ha_ovec_buf_t **b;
+	CALLOC(b, asm_opt.thread_num);
+	for (i = 0; i < asm_opt.thread_num; ++i)
+		b[i] = ha_ovec_init(1, 1);
+	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov); // build the index
+	kt_for(asm_opt.thread_num, worker_ov_final, b, R_INF.total_reads);
 	ha_pt_destroy(ha_idx);
 	ha_idx = 0;
-	///rescue_edges(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads, 4, 0.985);
+	for (i = 0; i < asm_opt.thread_num; ++i)
+		ha_ovec_destroy(b[i]);
+	free(b);
 }
 
 int ha_assemble(void)
 {
-	int r, ovlp_loaded = 0;
+	int r, hom_cov = -1, ovlp_loaded = 0;
 	if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
 		ovlp_loaded = 1;
-		fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
+		fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
+		if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN) && !(asm_opt.flag & HA_F_VERBOSE_GFA)) ha_triobin(&asm_opt);
+		if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();
+		if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
 	}
 	if (!ovlp_loaded) {
 		// construct hash table for high occurrence k-mers
-		if (!asm_opt.no_kmer_flt)
-			ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF);
+		if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
+			ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov);
+			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) {
-			clear_opt(&asm_opt, r); // this update asm_opt.roundID and a few other fields
+			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(),
 					yak_cpu_usage(), yak_peakrss_in_gb(), r + 1);
@@ -1236,17 +1184,16 @@ int ha_assemble(void)
 					asm_opt.num_bases, asm_opt.num_corrected_bases, asm_opt.num_recorrected_bases);
 			fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0);
 		}
-		Output_corrected_reads();
-		fprintf(stderr, "[M::%s::%.3f*%.2f] ==> written corrected reads to disk\n", __func__, yak_realtime(), yak_cpu_usage());
+		if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();
 		// overlap between corrected reads
-		clear_opt(&asm_opt, asm_opt.number_of_round);
+		ha_opt_reset_to_round(&asm_opt, asm_opt.number_of_round);
 		ha_overlap_final();
 		fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(),
 				yak_cpu_usage(), yak_peakrss_in_gb());
 		ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads);
 		ha_ft_destroy(ha_flt_tab);
-		Output_PAF();
-		trio_partition();
+		if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
+		ha_triobin(&asm_opt);
 	}
 	build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf, 
 			R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round, 

CommandLines.cpp

@@ -1,14 +1,23 @@
-#include "CommandLines.h"
 #include <zlib.h>
 #include <stdlib.h>
 #include <stdio.h>
-#include "ketopt.h"
 #include <sys/time.h>
+#include "CommandLines.h"
+#include "ketopt.h"
 
 #define DEFAULT_OUTPUT "hifiasm.asm"
 
 hifiasm_opt_t asm_opt;
 
+static ko_longopt_t long_options[] = {
+	{ "version",      ko_no_argument,       300 },
+	{ "dbg-gfa",      ko_no_argument,       301 },
+	{ "write-paf",    ko_no_argument,       302 },
+	{ "write-ec",     ko_no_argument,       303 },
+	{ "skip-triobin", ko_no_argument,       304 },
+	{ 0, 0, 0 }
+};
+
 double Get_T(void)
 {
   struct timeval t;
@@ -27,26 +36,24 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "    -a INT        round of assembly cleaning [%d]\n", asm_opt->clean_round);
     fprintf(stderr, "    -k INT        k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length);
 	fprintf(stderr, "    -w INT        minimizer window size [%d]\n", asm_opt->mz_win);
-	fprintf(stderr, "    -f INT        number of bits for bloom filter [%d]\n", asm_opt->bf_shift);
+	fprintf(stderr, "    -f INT        number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift);
 	fprintf(stderr, "    -D FLOAT      drop k-mers occuring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor);
-	fprintf(stderr, "    -N INT        consider up to INT overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
-    ///fprintf(stderr, "    -l            load all overlaps from disk, can avoid overlap calculation [%d]\n", asm_opt->load_index_from_disk);
-    ///fprintf(stderr, "    -i            ignore saved overlaps in *.ovlp*.bin files\n");
+	fprintf(stderr, "    -N INT        consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
     fprintf(stderr, "    -i            ignore saved overlaps in *.ovlp* files\n");
     fprintf(stderr, "    -z INT        length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
-    fprintf(stderr, "    -m INT        size of popped large bubbles for contig graph [%lld]\n", 
-    asm_opt->large_pop_bubble_size);
-    fprintf(stderr, "    -p INT        size of popped small bubbles for haplotype-resolved unitig graph [%lld]\n", 
-    asm_opt->small_pop_bubble_size);
+    fprintf(stderr, "    -m INT        size of popped large bubbles for contig graph [%lld]\n", asm_opt->large_pop_bubble_size);
+    fprintf(stderr, "    -p INT        size of popped small bubbles for haplotype-resolved unitig graph [%lld]\n", asm_opt->small_pop_bubble_size);
     fprintf(stderr, "    -n INT        small removed unitig threshold [%d]\n", asm_opt->max_short_tip);
     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, "    -v            show version number\n");
+    fprintf(stderr, "    --version     show version number\n");
     fprintf(stderr, "    -h            show help information\n");
 
     fprintf(stderr, "  Trio-partition:\n");
-    fprintf(stderr, "    -P FILE       paternal trio index generated by \"yak count\" [NULL]\n");
-    fprintf(stderr, "    -M FILE       maternal trio index generated by \"yak count\" [NULL]\n");
+    fprintf(stderr, "    -1 FILE       hap1/paternal k-mer dump generated by \"yak count\" []\n");
+    fprintf(stderr, "    -2 FILE       hap2/maternal k-mer dump generated by \"yak count\" []\n");
+	fprintf(stderr, "    -3 FILE       list of hap1/paternal read names []\n");
+	fprintf(stderr, "    -4 FILE       list of hap2/maternal read names []\n");
     fprintf(stderr, "    -c INT        lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt);
     fprintf(stderr, "    -d INT        upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt);
 
@@ -56,21 +63,20 @@ void Print_H(hifiasm_opt_t* asm_opt)
 
 void init_opt(hifiasm_opt_t* asm_opt)
 {
+	memset(asm_opt, 0, sizeof(hifiasm_opt_t));
+	asm_opt->flag = 0;
     asm_opt->coverage = -1;
     asm_opt->num_reads = 0;
     asm_opt->read_file_names = NULL;
     asm_opt->output_file_name = (char*)(DEFAULT_OUTPUT);
     asm_opt->required_read_name = NULL;
-    asm_opt->pat_index = NULL;
-    asm_opt->mat_index = NULL;
     asm_opt->thread_num = 1;
     asm_opt->k_mer_length = 51;
 	asm_opt->mz_win = 51;
 	asm_opt->bf_shift = 37;
 	asm_opt->high_factor = 5.0f;
-	asm_opt->no_HPC = 0;
-	asm_opt->no_kmer_flt = 0;
-	asm_opt->max_n_chain = 400;
+	asm_opt->hom_cov = 20;
+	asm_opt->max_n_chain = 100;
     asm_opt->k_mer_min_freq = 3;
     asm_opt->k_mer_max_freq = 66;
     asm_opt->load_index_from_disk = 1;
@@ -100,7 +106,7 @@ void destory_opt(hifiasm_opt_t* asm_opt)
     }
 }
 
-void clear_opt(hifiasm_opt_t* asm_opt, int round)
+void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round)
 {
     asm_opt->num_bases = 0;
     asm_opt->num_corrected_bases = 0;
@@ -109,7 +115,16 @@ void clear_opt(hifiasm_opt_t* asm_opt, int round)
     asm_opt->roundID = round;
 }
 
-int check_file(char* name, const char* opt)
+void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov)
+{
+	int max_n_chain = (int)(hom_cov * opt->high_factor + .499);
+	opt->hom_cov = hom_cov;
+	if (opt->max_n_chain < max_n_chain)
+		opt->max_n_chain = max_n_chain;
+	fprintf(stderr, "[M::%s] updated max_n_chain to %d\n", __func__, opt->max_n_chain);
+}
+
+static int check_file(char* name, const char* opt)
 {
     if(!name)
     {
@@ -244,8 +259,10 @@ int check_option(hifiasm_opt_t* asm_opt)
     }
 
 
-    if(asm_opt->pat_index != NULL && check_file(asm_opt->pat_index, "P") == 0) return 0;
-    if(asm_opt->mat_index != NULL && check_file(asm_opt->mat_index, "M") == 0) 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;
 
     // fprintf(stderr, "input file num: %d\n", asm_opt->num_reads);
     // fprintf(stderr, "output file: %s\n", asm_opt->output_file_name);
@@ -260,8 +277,6 @@ int check_option(hifiasm_opt_t* asm_opt)
     // fprintf(stderr, "size of popped large bubbles: %lld\n", asm_opt->large_pop_bubble_size);
     // fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip);
     // fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip);
-    // fprintf(stderr, "pat_index: %s\n", asm_opt->pat_index);
-    // fprintf(stderr, "mat_index: %s\n", asm_opt->mat_index);
     // fprintf(stderr, "min_cnt: %d\n", asm_opt->min_cnt);
     // fprintf(stderr, "mid_cnt: %d\n", asm_opt->mid_cnt);
 
@@ -301,13 +316,13 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 
     int c;
 
-    while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:lw:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:", 0)) >= 0) {
+    while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:w:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:1:2:3:4:", long_options)) >= 0) {
         if (c == 'h')
         {
             Print_H(asm_opt);
             return 0;
         } 
-        else if (c == 'v')
+        else if (c == 'v' || c == 300)
         {
 			puts(HA_VERSION);
             return 0;
@@ -318,23 +333,28 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 'r') asm_opt->number_of_round = atoi(opt.arg);
         else if (c == 'k') asm_opt->k_mer_length = atoi(opt.arg);
         else if (c == 'i') asm_opt->load_index_from_disk = 0; 
-        else if (c == 'l') asm_opt->load_index_from_disk = 1; 
         else if (c == 'w') asm_opt->mz_win = atoi(opt.arg);
 		else if (c == 'D') asm_opt->high_factor = atof(opt.arg);
-		else if (c == 'F') asm_opt->no_kmer_flt = 1;
+		else if (c == 'F') asm_opt->flag |= HA_F_NO_KMER_FLT;
 		else if (c == 'N') asm_opt->max_n_chain = atoi(opt.arg);
         else if (c == 'a') asm_opt->clean_round = atoi(opt.arg); 
         else if (c == 'z') asm_opt->adapterLen = atoi(opt.arg);
         else if (c == 'b') asm_opt->required_read_name = opt.arg;
         else if (c == 'c') asm_opt->min_cnt = atoi(opt.arg);
         else if (c == 'd') asm_opt->mid_cnt = atoi(opt.arg);
-        else if (c == 'P') asm_opt->pat_index = opt.arg;
-        else if (c == 'M') asm_opt->mat_index = opt.arg;
+        else if (c == '1' || c == 'P') asm_opt->fn_bin_yak[0] = opt.arg; // -P/-M reserved for backward compatibility
+        else if (c == '2' || c == 'M') asm_opt->fn_bin_yak[1] = opt.arg;
+        else if (c == '3') asm_opt->fn_bin_list[0] = opt.arg;
+        else if (c == '4') asm_opt->fn_bin_list[1] = opt.arg;
         else if (c == 'x') asm_opt->max_drop_rate = atof(opt.arg);
         else if (c == 'y') asm_opt->min_drop_rate = atof(opt.arg);
         else if (c == 'p') asm_opt->small_pop_bubble_size = atoll(opt.arg);
         else if (c == 'm') asm_opt->large_pop_bubble_size = atoll(opt.arg);
         else if (c == 'n') asm_opt->max_short_tip = atoll(opt.arg);
+		else if (c == 301) asm_opt->flag |= HA_F_VERBOSE_GFA;
+		else if (c == 302) asm_opt->flag |= HA_F_WRITE_PAF;
+		else if (c == 303) asm_opt->flag |= HA_F_WRITE_EC;
+		else if (c == 304) asm_opt->flag |= HA_F_SKIP_TRIOBIN;
         else if (c == ':') 
         {
 			fprintf(stderr, "[ERROR] missing option argument in \"%s\"\n", argv[opt.i - 1]);

CommandLines.h

@@ -3,26 +3,32 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.3.0-dirty-r192"
+#define HA_VERSION "0.4.0"
 
 #define VERBOSE 0
-#define VERBOSE_GFA 0
+
+#define HA_F_NO_HPC       0x1
+#define HA_F_NO_KMER_FLT  0x2
+#define HA_F_VERBOSE_GFA  0x4
+#define HA_F_WRITE_EC     0x8
+#define HA_F_WRITE_PAF    0x10
+#define HA_F_SKIP_TRIOBIN 0x20
 
 typedef struct {
+	int flag;
     int num_reads;
     char** read_file_names;
     char* output_file_name;
     char* required_read_name;
-    char* pat_index;
-    char* mat_index;
+	char *fn_bin_yak[2];
+	char *fn_bin_list[2];
     int thread_num;
     int k_mer_length;
 	int mz_win;
 	int bf_shift;
-	float high_factor;
-	int no_HPC;
-	int no_kmer_flt;
-	int max_n_chain;
+	float high_factor; // coverage cutoff set to high_factor*hom_cov
+	int hom_cov;
+	int max_n_chain; // fall-back max number of chains to consider
     int k_mer_min_freq;
     int k_mer_max_freq;
     int load_index_from_disk;
@@ -56,8 +62,14 @@ extern hifiasm_opt_t asm_opt;
 
 void init_opt(hifiasm_opt_t* asm_opt);
 void destory_opt(hifiasm_opt_t* asm_opt);
-void clear_opt(hifiasm_opt_t* asm_opt, int round);
+void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round);
+void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov);
 int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt);
 double Get_T(void);
 
+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]));
+}
+
 #endif

Hash_Table.cpp

@@ -3,9 +3,6 @@
 #include <stdint.h>
 #include <pthread.h>
 #include "Hash_Table.h"
-#include "Process_Read.h"
-#include "Correct.h"
-#include "CommandLines.h"
 #include "ksort.h"
 pthread_mutex_t output_mutex;
 
@@ -347,7 +344,7 @@ static int32_t ha_kmer_hit_lis(int32_t n, const k_mer_hit *a, int32_t *b, int32_
 		int32_t lo = 1, hi = L, newL;
 		while (lo <= hi) {
 			int32_t mid = (lo + hi + 1) >> 1;
-			if (a[M[mid]].offset < a[i].offset) lo = mid + 1;
+			if (a[M[mid]].self_offset < a[i].self_offset) lo = mid + 1;
 			else hi = mid - 1;
 		}
 		newL = lo, P[i] = M[newL - 1], M[newL] = i;
@@ -364,16 +361,19 @@ int32_t ha_chain_lis_core(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min
 	int32_t *tmp = (int32_t*)dp->tmp;
 	int32_t i, m, *b = tmp, *M = tmp + n_a;
 	int32_t tot_indel = 0, tot_len = 0;
+	double bw_pen;
 	if (n_a < 2) return -1;
 	for (i = 1; i < n_a; ++i)
-		if (a[i-1].offset >= a[i].offset)
+		if (a[i-1].self_offset >= a[i].self_offset)
 			break;
 	if (i == n_a) {
 		for (i = 0; i < n_a; ++i)
 			b[i] = i;
 		m = n_a;
 	} else m = ha_kmer_hit_lis(n_a, a, b, M);
-	dp->score[0] = 0, dp->pre[0] = -1, dp->indels[0] = 0, dp->self_length[0] = 0;
+	bw_pen = 1.0 / bw_thres;
+	dp->score[0] = a[b[0]].good? min_sc : min_sc>>1;
+	dp->pre[0] = -1, dp->indels[0] = 0, dp->self_length[0] = 0;
 	for (i = 1; i < m; ++i) {
 		int32_t j0 = b[i-1], j1 = b[i], score, dg;
 		int32_t dx = (int32_t)a[j1].offset - (int32_t)a[j0].offset;
@@ -382,12 +382,13 @@ int32_t ha_chain_lis_core(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min
 		double gap_rate;
 		tot_indel += dd;
 		tot_len += dy;
-		if (tot_indel > tot_len * bw_thres)
-			break;
+		if (tot_indel > tot_len * bw_thres) break;
 		dg = dx < dy? dx : dy;
+		if (dd > THRESHOLD_MAX_SIZE && dd > dg * bw_thres) break;
 		score = dg < min_sc? dg : min_sc;
+		if (!a[j1].good) score >>= 1;
 		gap_rate = (double)tot_indel / tot_len;
-		score -= (int)(gap_rate * score * bw_thres);
+		score -= (int)(gap_rate * score * bw_pen);
 		dp->score[i] = dp->score[i-1] + score;
 		dp->pre[i] = i - 1;
 		dp->indels[i] = tot_indel;
@@ -431,7 +432,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
         pos = a[i].offset;
         self_pos = a[i].self_offset;
         max_j = -1;
-        max_score = min_score;
+        max_score = a[i].good? min_score : min_score>>1;
         max_indels = 0;
         max_self_length = 0;
 
@@ -460,6 +461,7 @@ 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;
 
             gap_rate = (double)((double)(total_indels)/(double)(total_self_length));
             ///if the gap rate > 0.06, score will be negative
@@ -525,7 +527,7 @@ skip_dp:
     }
 
     clear_fake_cigar(&(result->f_cigar));
-    ///not a has been sorted by offset, that means has been sorted by query offset
+    ///note a has been sorted by offset, that means has been sorted by query offset
     i = max_i;
     result->x_pos_e = a[i].self_offset;
     result->y_pos_e = a[i].offset;

Hash_Table.h

@@ -92,7 +92,7 @@ typedef struct
     uint32_t y_pos_strand;
 
     uint32_t overlapLen;
-    uint32_t shared_seed;
+    int32_t shared_seed;
     uint32_t align_length;
     uint8_t is_match;
     uint8_t without_large_indel;
@@ -117,7 +117,7 @@ typedef struct
 
 typedef struct
 {
-	uint32_t readID:31, strand:1;
+	uint32_t readID:30, strand:1, good:1;
 	uint32_t offset, self_offset;
 } k_mer_hit;
 
@@ -193,5 +193,6 @@ void init_window_list_alloc(window_list_alloc* x);
 void clear_window_list_alloc(window_list_alloc* x);
 void destory_window_list_alloc(window_list_alloc* x);
 void resize_window_list_alloc(window_list_alloc* x, long long size);
+void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, double band_width_threshold, int max_skip, int x_readLen, int y_readLen);
 
 #endif

Makefile

@@ -2,8 +2,8 @@ CXX=		g++
 CXXFLAGS=	-g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall
 CPPFLAGS=
 INCLUDES=
-OBJS=		Output.o CommandLines.o Process_Read.o Assembly.o Hash_Table.o \
-			POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o \
+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 sys.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm
@@ -34,19 +34,23 @@ depend:
 
 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: Output.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
 Correct.o: kdq.h CommandLines.h Levenshtein_distance.h POA.h Assembly.h
 Hash_Table.o: Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h
-Hash_Table.o: CommandLines.h Correct.h Levenshtein_distance.h POA.h ksort.h
+Hash_Table.o: CommandLines.h ksort.h
 Levenshtein_distance.o: Levenshtein_distance.h
 Output.o: Output.h CommandLines.h
 Overlaps.o: Overlaps.h kvec.h kdq.h ksort.h Process_Read.h CommandLines.h
 Overlaps.o: Hash_Table.h htab.h Correct.h Levenshtein_distance.h POA.h
+Overlaps.o: Purge_Dups.h
 POA.o: POA.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h
 POA.o: CommandLines.h Correct.h Levenshtein_distance.h
 Process_Read.o: Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
+Purge_Dups.o: ksort.h Purge_Dups.h kvec.h kdq.h Overlaps.h Hash_Table.h
+Purge_Dups.o: htab.h Process_Read.h CommandLines.h Correct.h
+Purge_Dups.o: Levenshtein_distance.h POA.h kthread.h
 Trio.o: khashl.h kthread.h Process_Read.h Overlaps.h kvec.h kdq.h
 Trio.o: CommandLines.h htab.h
 anchor.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h

Overlaps.h

@@ -68,8 +68,7 @@ typedef struct {
 
 void init_ma_hit_t_alloc(ma_hit_t_alloc* x);
 void clear_ma_hit_t_alloc(ma_hit_t_alloc* x);
-void resize_ma_hit_t_alloc(ma_hit_t_alloc* x, uint64_t size);
-void destory_all_ma_hit_t_alloc(ma_hit_t_alloc* x, uint64_t n_read);
+void resize_ma_hit_t_alloc(ma_hit_t_alloc* x, uint32_t size);
 void destory_ma_hit_t_alloc(ma_hit_t_alloc* x);
 void add_ma_hit_t_alloc(ma_hit_t_alloc* x, ma_hit_t* element);
 void ma_hit_sort_tn(ma_hit_t *a, long long n);
@@ -134,6 +133,16 @@ typedef struct {
 	ma_utg_t* F_seq;
 } asg_t;
 
+asg_t *asg_init(void);
+void asg_destroy(asg_t *g);
+void asg_arc_sort(asg_t *g);
+void asg_seq_set(asg_t *g, int sid, int len, int del);
+void asg_arc_index(asg_t *g);
+void asg_cleanup(asg_t *g);
+void asg_symm(asg_t *g);
+void print_gfa(asg_t *g);
+
+
 typedef struct { size_t n, m; uint64_t *a; } asg64_v;
 
 
@@ -355,11 +364,21 @@ typedef struct {
 	uint32_t pre_n_seq, seqID; 
 } C_graph;
 
+typedef struct {
+	kvec_t(uint8_t) a;
+	uint32_t i;
+} kvec_t_u8_warp;
+
 typedef struct {
 	kvec_t(uint32_t) a;
 	uint32_t i;
 } kvec_t_u32_warp;
 
+typedef struct {
+	kvec_t(int32_t) a;
+	uint32_t i;
+} kvec_t_i32_warp;
+
 typedef struct {
 	kvec_t(uint64_t) a;
 	uint64_t i;
@@ -370,6 +389,7 @@ typedef struct {
 	uint64_t i;
 }kvec_asg_arc_t_warp;
 
+void sort_kvec_t_u64_warp(kvec_t_u64_warp* u_vecs, uint32_t is_descend);
 
 
 typedef struct {
@@ -774,7 +794,10 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
         	{
 				qn = node_min->a[b_min.readI]>>33;
 
-				if(reverse_sources[qn].length>=0) min_count++;
+				/************************BUG: don't forget****************************/
+				if(reverse_sources[qn].length > 0) min_count++;
+				///if(reverse_sources[qn].length >= 0) min_count++;
+				/************************BUG: don't forget****************************/
 				for (j = 0; j < (long long)reverse_sources[qn].length; j++)
 				{
 					tn = Get_tn(reverse_sources[qn].buffer[j]);
@@ -823,7 +846,10 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 		{
 			qn = (b_min.b_0->b.a[b_min.untigI]>>1);
 
-			if(reverse_sources[qn].length>=0) min_count++;
+			/************************BUG: don't forget****************************/
+			if(reverse_sources[qn].length > 0) min_count++;
+			///if(reverse_sources[qn].length >= 0) min_count++;
+			/************************BUG: don't forget****************************/
 
 			for (j = 0; j < (long long)reverse_sources[qn].length; j++)
         	{
@@ -867,6 +893,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 }
 
 
+
 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)
@@ -920,7 +947,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
         	{
 				qn = node_min->a[b_min.readI]>>33;
 
-				if(reverse_sources[qn].length>=0) min_count++;
+				if(reverse_sources[qn].length > 0) min_count++;
 				for (j = 0; j < (long long)reverse_sources[qn].length; j++)
 				{
 					tn = Get_tn(reverse_sources[qn].buffer[j]);
@@ -957,7 +984,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold)
 		{
 			qn = (b_min.b_0->b.a[b_min.untigI]>>1);
 
-			if(reverse_sources[qn].length>=0) min_count++;
+			if(reverse_sources[qn].length > 0) min_count++;
 
 			for (j = 0; j < (long long)reverse_sources[qn].length; j++)
         	{
@@ -1017,7 +1044,7 @@ void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_allo
 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);
+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, uint8_t postive_flag, float drop_rate, 
 ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
 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);
@@ -1029,5 +1056,11 @@ R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t bac
 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);
+int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio);
+
+#define JUNK_COV 5
+#define DISCARD_RATE 0.8
 
 #endif

Process_Read.cpp

@@ -42,11 +42,10 @@ void destory_All_reads(All_reads* r)
 {
 	uint64_t i = 0;
 	for (i = 0; i < r->total_reads; i++) {
-		if (r->N_site[i] != NULL)
-			free(r->N_site[i]);
-		free(r->read_sperate[i]);
-		if (r->paf) free(r->paf[i].buffer);
-		if (r->reverse_paf) free(r->reverse_paf[i].buffer);
+		if (r->N_site[i]) free(r->N_site[i]);
+		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);
 	}
 	free(r->paf);
 	free(r->reverse_paf);

Purge_Dups.h

@@ -0,0 +1,19 @@
+#ifndef __PURGEDUPS__
+#define __PURGEDUPS__
+#include <stdio.h>
+#include <stdint.h>
+#include "kvec.h"
+#include "kdq.h"
+#include "Overlaps.h"
+#include "Hash_Table.h"
+
+void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* reverse_sources, 
+R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t bi_graph_Len, uint32_t long_hap_overlap, 
+float lable_match_rate, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio,
+uint32_t just_contain);
+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 enable_debug_mode(uint32_t mode);
+
+#endif

README.md

@@ -21,53 +21,40 @@ outputs consist of:
    (*prefix*.r\_utg.gfa). This graph keeps all haplotype information, including
    somatic mutations and recurrent sequencing errors.
 2. Haplotype-resolved processed unitig graph without small bubbles
-   (*prefix*.p\_utg.gfa). This is usually the preferred output for highly
-   heterozygous genomes.
-3. Primary assembly [contig][unitig] graph (*prefix*.p\_ctg.gfa). This is the
-   preferred output for inbred strains or human. For highly heterozygous
-   genomes, this graph may represent multiple haplotypes. We plan to change
-   this to represent one set of haplotypes.
-4. Alternate assembly contig graph (*prefix*.a\_ctg.gfa).
-5. Haplotype-aware error corrected reads in fasta format (*prefix*.ec.fa).
-6. All-to-all overlaps in the [PAF][paf] format (*prefix*.ovlp.paf).
+   (*prefix*.p\_utg.gfa). Small bubbles might be caused by somatic mutations or noise in data, 
+   which are not the real haplotype information.
+3. Primary assembly [contig][unitig] graph (*prefix*.p\_ctg.gfa). This graph collapses different
+   haplotypes.
+4. Alternate assembly contig graph (*prefix*.a\_ctg.gfa). This graph consists of all assemblies that
+   are discarded in primary contig graph.
 
 For trio assembly, the input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and the paternal/maternal trio indexes generated by `yak count` (see https://github.com/lh3/yak). The outputs consist of:
 1. Haplotype-resolved raw [unitig][unitig] graph in [GFA][gfa] format
    (*prefix*.r\_utg.gfa). This graph keeps all haplotype information. 
 
-2. Phased maternal unitig graph (*prefix*.m.r\_utg.gfa).
-This graph keeps the phased maternal assembly.
+2. Phased paternal/haplotype1 contig graph (*prefix*.hap1.p\_ctg.gfa). This graph keeps the phased
+   paternal/haplotype1 assembly.
 
-3. Phased paternal unitig graph (*prefix*.p.r\_utg.gfa).
-This graph keeps the phased paternal assembly.
+3. Phased maternal/haplotype2 contig graph (*prefix*.hap2.p\_ctg.gfa). This graph keeps the phased
+   maternal/haplotype2 assembly.
 
 
-In addition, hifiasm also outputs three binary files that save all overlap inforamtion
-(hifiasm.asm.ovlp, hifiasm.asm.ovlp.source, hifiasm.asm.ovlp.reverse in default). With these files, hifiasm can avoid the time-consuming all-to-all overlap calculation step, and do the assembly
+
+In addition, hifiasm also outputs three binary files that save all overlap information (*prefix*.ec.bin, *prefix*.ovlp.reverse.bin, *prefix*.ovlp.source.bin). With these files, hifiasm can avoid the time-consuming all-to-all overlap calculation step, and do the assembly
 directly and quickly. This might be helpful when you want to get an optimized
 assembly by multiple rounds of experiments with different parameters.
 
 Hifiasm is a standalone and lightweight assembler, which does not need external
 libraries (except zlib). For large genomes, it can generate high-quality
-assembly in a few hours. Hifiasm has been tested on the following datasets:
+assembly in a few hours. Hifiasm has been tested on human, butterfly, rice and drosophila. 
+In particular, hifiasm is able to assemble the 26.5Gb California redwood tree in a few days. 
+The results are as follows: 
 
 |<sub>Dataset<sub>|<sub>GSize<sub>|<sub>Cov<sub>|<sub>Asm options<sub>|<sub>CPU time<sub>|<sub>Wall time<sub>|<sub>RAM<sub>|<sub>[unitig][unitig]/[contig][unitig] N50<sup>[1]</sup><sub>|
 |:---------------|-----:|-----:|:---------------------|-------:|--------:|----:|----------------:|
-|<sub>[Human NA12878]<sub>|<sub>3Gb<sub>|<sub>x28<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>200h<sub>|    <sub>5h32m<sub>|<sub>114G<sub>|<sub>93.5Kb/28.2Mb<sub>|
-|<sub>[Human HG002]<sub>|<sub>3Gb<sub>|<sub>x43<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>405h10m<sub>|<sub>12h7m<sub>|<sub>146G<sub>|<sub>320kb/46.0Mb<sub>|
-|<sub>[Human CHM13]<sub>|<sub>3Gb<sub>|<sub>x27<sub>|<sub>-k 40 -t 42 -r 2<sub>|<sub>157h28m<sub>|<sub>5h10m<sub>|<sub>85.8G<sub>|<sub>NA<sup>[2]</sup>/41.4Mb<sub>|
-|<sub>[Butterfly]<sub>|<sub>358Mb<sub>|<sub>x35<sub>|<sub>-k 40 -t 42 -r 2 -z 20<sub>|<sub>17h6m<sub>|<sub>36m<sub>|<sub>16G<sub>|<sub>7.5Mb/NA<sup>[3]</sup><sub>|
 |<sub>[\[Redwood\]](https://downloads.pacbcloud.com/public/dataset/redwood2020/)<sub>|<sub>26.5Gb<sub>|<sub>x23<sub>|<sub>-k 40 -t 64 -r 2<sub>|<sub>7274h30m<sub>|<sub>141h30m<sub>|<sub>512G<sub>|<sub>1.7Mb/1.9Mb<sub>|
 
-<sub>[1] unitig N50 is the N50 of assembly graph with haplotype information (i.e., bubbles), while the contig N50 is the N50 of haplotype collapsed assembly (i.e., without bubbles).
-[2] CHM13 is a homozygous sample, so that unitig N50 makes no sense.
-[3] Butterfly has high heterozygous rate, so that most chromosomes have been fully separated into two haplotypes. In this case, contig N50 makes no sense.<sub>
-
-Note that different species need different assembly graphs. For homozygous genomes (i.e., Human CHM13), the primary assembly contig graph is the best choice. 
-For species with high heterozygous rate (i.e., Butterfly), different haplotypes can be fully separated. It is important to remove small bubbles from the haplotype-resolved unitig graph. The
-reason is that some small bubbles are caused by somatic mutations or noise in data, which are not
-the real haplotype information. In this case, haplotype-resolved processed unitig graph
-without small bubbles should be better. For ordinary human genome (i.e., Human NA12878 and HG002), different haplotypes cannot be fully separated due to the low heterozygous rate. There are many small bubbles including haplotype information, which cannot be simply removed. Thus, it is necessary to use the haplotype-resolved raw unitig graph. **Hifiasm will generate a universal haplotype-resolved contig graph for all species in the near future.**
+<sub>[1] unitig N50 is the N50 of assembly graph with haplotype information (i.e., bubbles), while the contig N50 is the N50 of haplotype collapsed assembly (i.e., without bubbles).<sub>
 
 ## Usage
 
@@ -78,14 +65,12 @@ For Hifi reads assembly, a typical command line looks like:
 ```
 
 where `NA12878.fq.gz` is the input reads and `-o` specifies the output files.
-In this example, all output files can be found at `NA12878.asm.*`. `-k`, `-t`
-and `-r` specify the length of k-mer, the number of CPU threads, and the number
-of correction rounds, respectively. Note that at first run, hifiasm will save
-all overlaps to disk, which can avoid the time-consuming all-to-all overlap
-calculation next time. For hifiasm, once the overlap information has been
-obtained during the previous run in advance, it is able to load all overlaps
-from disk and then directly do assembly. If you want to ignore the pre-computed
-overlap information, please specify `-i`.
+In this example, all output files can be found at `NA12878.asm.*`. `-t` specifies 
+the number of CPU threads. Note that at first run, hifiasm will save all overlaps 
+to disk, which can avoid the time-consuming all-to-all overlap calculation next time. 
+For hifiasm, once the overlap information has been obtained during the previous run 
+in advance, it is able to load all overlaps from disk and then directly do assembly. 
+If you want to ignore the pre-computed overlap information, please specify `-i`.
 
 Please note that some old Hifi reads may consist of short adapters. To improve
 the assembly quality, adapters should be removed by `-z` as follow:
@@ -109,7 +94,7 @@ For trio assembly, first the trio indexes of paternal/maternal should be generat
 and then run hifiasm as follow:
 
 ```sh
-./hifiasm -o NA12878.asm -t 32 -P pat.yak -M mat.yak NA12878_1.fq.gz NA12878_2.fq.gz
+./hifiasm -o NA12878.asm -t 32 -1 pat.yak -2 mat.yak NA12878_1.fq.gz NA12878_2.fq.gz
 ```
 
 [unitig]: http://wgs-assembler.sourceforge.net/wiki/index.php/Celera_Assembler_Terminology
@@ -124,13 +109,6 @@ have further questions, please raise an issue at the issue page.
 
 ## Limitations and future works
 
-1. For genome with low heterozygous rate, hifiasm only outputs
-   haplotype-resolved assembly graph, instead of the phased chromosome-level
-   assembly (will support such output in future).
+1. The running time and memory usage should be further reduced.
 
-2. For different species, hifiasm outputs different assembly graphs, which are not easy to use.
-   Hifiasm will generate a universal haplotype-resolved contig graph for all species in future.
-
-3. The running time and memory usage should be further reduced.
-
-4. The N50 should be further improved. 
+2. The N50 should be further improved. 

Trio.cpp

@@ -3,8 +3,10 @@
 #include <stdarg.h>
 #include <string.h>
 #include <assert.h>
+#include <zlib.h>
 #include "khashl.h" // hash table
 #include "kthread.h"
+#include "kseq.h"
 #include "Process_Read.h"
 #include "htab.h"
 #include "CommandLines.h"
@@ -24,6 +26,11 @@
 #define yak_ch_hash(a) ((a)>>YAK_COUNTER_BITS)
 KHASHL_SET_INIT(static klib_unused, yak_ht_t, yak_ht, uint64_t, yak_ch_hash, yak_ch_eq)
 
+typedef const char *ha_cstr_t;
+KHASHL_MAP_INIT(static klib_unused, cstr_ht_t, cstr_ht, ha_cstr_t, int64_t, kh_hash_str, kh_eq_str)
+
+KSTREAM_INIT(gzFile, gzread, 65536)
+
 typedef struct {
 	struct yak_ht_t *h;
 } yak_ch1_t;
@@ -248,26 +255,18 @@ static void tb_worker(void *_data, long k, int tid)
 	if(type == 'm') aux->seq->trio_flag[k] = MOTHER;
 }
 
-void trio_partition(void)
+static void ha_triobin_yak(const hifiasm_opt_t *opt)
 {
-    if (asm_opt.pat_index == NULL || asm_opt.mat_index == NULL) {
-		memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t));
-		return;
-	}
-
-	double start_time = Get_T();
-	fprintf(stderr, "Start trio binning ...... \n");
-
     yak_ch_t *ch;
     int i /**, min_cnt = 2, mid_cnt = 5**/;
     tb_shared_t aux;
     memset(&aux, 0, sizeof(tb_shared_t));
-	aux.n_threads = asm_opt.thread_num, aux.print_diff = 0;
+	aux.n_threads = opt->thread_num, aux.print_diff = 0;
 	aux.ratio_thres = 0.33;
 	aux.seq = &R_INF;
 
-    ch = yak_ch_restore_core(0,  asm_opt.pat_index, YAK_LOAD_TRIOBIN1, asm_opt.min_cnt, asm_opt.mid_cnt);
-	ch = yak_ch_restore_core(ch, asm_opt.mat_index, YAK_LOAD_TRIOBIN2, asm_opt.min_cnt, asm_opt.mid_cnt);
+    ch = yak_ch_restore_core(0,  opt->fn_bin_yak[0], YAK_LOAD_TRIOBIN1, opt->min_cnt, opt->mid_cnt);
+	ch = yak_ch_restore_core(ch, opt->fn_bin_yak[1], YAK_LOAD_TRIOBIN2, opt->min_cnt, opt->mid_cnt);
 
     aux.k = ch->k;
     aux.ch = ch;
@@ -286,6 +285,72 @@ void trio_partition(void)
 	free(aux.bseq);
     yak_ch_destroy(ch);
 
-	fprintf(stderr, "Trio binning has been done.\n");
-	fprintf(stderr, "%-30s%18.2f\n\n", "Trio binning time:", Get_T() - start_time); 
+	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> partitioned reads using yak dumps\n", __func__, yak_realtime(), yak_cpu_usage());
+}
+
+static int ha_triobin_set_list(const cstr_ht_t *h, const char *fn, int flag)
+{
+	gzFile fp;
+	kstream_t *ks;
+	kstring_t str = {0,0,0};
+	int dret;
+	int64_t n_tot = 0, n_bin = 0;
+	fp = gzopen(fn, "r");
+	if (fp == 0) {
+		fprintf(stderr, "ERROR: failed to open file '%s'\n", fn);
+		return -1;
+	}
+	ks = ks_init(fp);
+	while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) {
+		char *p;
+		khint_t k;
+		++n_tot;
+		for (p = str.s; *p; ++p)
+			if (*p == '\t' || *p == ' ')
+				*p = 0;
+		k = cstr_ht_get(h, str.s);
+		if (k != kh_end(h)) {
+			R_INF.trio_flag[kh_val(h, k)] = flag;
+			++n_bin;
+		}
+	}
+	free(str.s);
+	ks_destroy(ks);
+	gzclose(fp);
+	fprintf(stderr, "[M::%s::%.3f*%.2f] flagged %ld reads, out of %ld lines in file '%s'\n",
+			__func__, yak_realtime(), yak_cpu_usage(), (long)n_bin, (long)n_tot, fn);
+	return 0;
+}
+
+static void ha_triobin_list(const hifiasm_opt_t *opt)
+{
+	int64_t i;
+	khint_t k;
+	cstr_ht_t *h;
+	assert(R_INF.total_reads < (uint32_t)-1);
+	h = cstr_ht_init();
+	for (i = 0; i < (int64_t)R_INF.total_reads; ++i) {
+		int absent;
+		char *str = (char*)calloc(Get_NAME_LENGTH(R_INF, i) + 1, 1);
+		strncpy(str, Get_NAME(R_INF, i), Get_NAME_LENGTH(R_INF, i));
+		k = cstr_ht_put(h, str, &absent);
+		if (absent) kh_val(h, k) = i;
+	}
+	fprintf(stderr, "[M::%s::%.3f*%.2f] created the hash table for read names\n", __func__, yak_realtime(), yak_cpu_usage());
+	ha_triobin_set_list(h, opt->fn_bin_list[0], FATHER);
+	ha_triobin_set_list(h, opt->fn_bin_list[1], MOTHER);
+	for (k = 0; k < kh_end(h); ++k)
+		if (kh_exist(h, k))
+			free((char*)kh_key(h, k));
+	cstr_ht_destroy(h);
+	fprintf(stderr, "[M::%s::%.3f*%.2f] ==> partitioned reads with external lists\n", __func__, yak_realtime(), yak_cpu_usage());
+}
+
+void ha_triobin(const hifiasm_opt_t *opt)
+{
+	memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads * sizeof(uint8_t));
+	if (opt->fn_bin_list[0] && opt->fn_bin_list[1])
+		ha_triobin_list(opt);
+	if (opt->fn_bin_yak[0] && opt->fn_bin_yak[1])
+		ha_triobin_yak(opt);
 }

anchor.cpp

@@ -3,9 +3,11 @@
 #include "ksort.h"
 #include "Hash_Table.h"
 
+#define HA_KMER_GOOD_RATIO 0.333
+
 typedef struct { // this struct is not strictly necessary; we can use k_mer_pos instead, with modifications
 	uint64_t srt;
-	uint32_t self_off;
+	uint32_t self_off:31, good:1;
 	uint32_t other_off;
 } anchor1_t;
 
@@ -21,7 +23,7 @@ KSORT_INIT(or_xs, overlap_region, oreg_xs_lt)
 KSORT_INIT(or_ss, overlap_region, oreg_ss_lt)
 
 typedef struct {
-	int n;
+	int n, good;
 	const ha_idxpos_t *a;
 } seed1_t;
 
@@ -54,6 +56,8 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	extern ha_pt_t *ha_idx;
 	uint32_t i;
 	uint64_t k, l;
+	double low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
+	double high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
 
 	// prepare
     clear_Candidates_list(cl);
@@ -62,14 +66,17 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	ab->mz.n = 0, ab->n_a = 0;
 
 	// get the list of anchors
-	ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !asm_opt.no_HPC, &ab->mz, ha_flt_tab);
+	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);
 	if (ab->mz.m > ab->old_mz_m) {
 		ab->old_mz_m = ab->mz.m;
 		REALLOC(ab->seed, ab->old_mz_m);
 	}
 	for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) {
-		ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &ab->seed[i].n);
-		ab->n_a += ab->seed[i].n;
+		int n;
+		ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &n);
+		ab->seed[i].n = n;
+		ab->seed[i].good = (n > low_occ && n < high_occ);
+		ab->n_a += n;
 	}
 	if (ab->n_a > ab->m_a) {
 		ab->m_a = ab->n_a;
@@ -86,6 +93,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 			uint8_t rev = z->rev == y->rev? 0 : 1;
 			an->other_off = y->pos;
 			an->self_off = rev? ucr->length - 1 - (z->pos + 1 - z->span) : z->pos;
+			an->good = s->good;
 			an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off;
 		}
 	}
@@ -111,6 +119,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 		p->strand = ab->a[k].srt >> 32 & 1;
 		p->offset = ab->a[k].other_off;
 		p->self_offset = ab->a[k].self_off;
+		p->good = ab->a[k].good;
 	}
 	cl->length = ab->n_a;
 
@@ -128,7 +137,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	#endif
 
 	if ((int)overlap_list->length > max_n_chain) {
-		uint32_t n[2], s[2];
+		int32_t n[2], s[2];
 		n[0] = n[1] = 0, s[0] = s[1] = 0;
 		for (i = 0; i < (uint32_t)overlap_list->length; ++i) {
 			const overlap_region *r = &overlap_list->list[i];

hifiasm.1

@@ -1,38 +1,60 @@
-.TH hifiasm 1 "22 Mar 2020" "hifiasm-0.3.0" "Bioinformatics tools"
+.TH hifiasm 1 "12 Apr 2020" "hifiasm-0.4.0" "Bioinformatics tools"
 
 .SH NAME
 .PP
 hifiasm - haplotype-resolved de novo assembler for PacBio Hifi reads.
 
 .SH SYNOPSIS
-.PP
-hifiasm
+
+* Assemble HiFi reads:
+.RS 4
+.B hifiasm
 .RB [ -o
-.IR outPrefix ]
+.IR prefix ]
 .RB [ -t
-.IR numThres ]
-.RB [ -r
-.IR roundCorrection ]
-.RB [ -a
-.IR roundGraphClean ]
+.IR nThreads ]
+.RB [ -z
+.IR endTrimLen ]
+.R [options]
+.I input1.fq
+.RI [ input2.fq
+.R [...]]
+.RE
+
+* Trio binning assembly with yak dumps:
+.RS 4
+.B yak count
+.B -o
+.I paternal.yak
+.B -b37
+.RB [ -t
+.IR nThreads ]
 .RB [ -k
 .IR kmerLen ]
-.RB [ -z
-.IR adapterLen ]
-.RB [ -m
-.IR maxLargeBubbles ]
-.RB [ -p
-.IR maxSmallBubbles ]
-.RB [ -n
-.IR maxSmallUnitig ]
-.RB [ -x
-.IR maxDropRatio ]
-.RB [ -y
-.IR minDropRatio ]
-.RB [ -i ]
-.RB [ -v ]
-.RB [ -h ]
-.I <in_1.fq> <in_2.fq> <...>
+.I paternal.fq.gz
+.br
+.B yak count
+.B -o
+.I maternal.yak
+.B -b37
+.RB [ -t
+.IR nThreads ]
+.RB [ -k
+.IR kmerLen ]
+.I maternal.fq.gz
+.br
+.B hifiasm
+.RB [ -o
+.IR prefix ]
+.RB [ -t
+.IR nThreads ]
+.R [options]
+.B -1
+.I paternal.yak
+.B -2
+.I maternal.yak
+.I child.hifi.fq.gz
+.RE
 
 .SH DESCRIPTION
 .PP
@@ -49,49 +71,59 @@ outputs consist of multiple types of assembly graph in GFA format.
 
 .TP 10
 .BI -o \ FILE
-Prefix of output files [hifiasm.asm]. The outputs of hifiasm include error corrected
-reads in fasta format, all-to-all overlaps in paf format, and four types of assembly
-graph in GFA format. For detailed description of all assembly graphs, please see
-.I 'Outputs'
+Prefix of output files [hifiasm.asm]. For detailed description of all assembly
+graphs, please see the
+.B OUTPUTS
 section of this man-page. 
 
 .TP 10
 .BI -t \ INT
 Number of CPU threads used by hifiasm [1]. 
 
+.TP
+.BI -h
+Show help information.
 
-.TP 10
-.BI -v
+.TP
+.BI --version
 Show version number. 
 
-.TP 10
-.BI -h
-Show help information. 
 
 .SS Error correction options
 
 .TP 10
 .BI -k \ INT
-K-mer length [40]. This option must be less than 64.
+K-mer length [51]. This option must be less than 64.
 
-.TP 10
+.TP
+.BI -w \ INT
+Minimizer window size [51].
+
+.TP
+.BI -f \ INT
+Number of bits for bloom filter; 0 to disable [37]. This bloom filter is used
+to filter out singleton k-mers when counting all k-mers. It takes
+.RI 2^( INT -3)
+bytes of memory. A proper setting saves memory. 37 is recommended for human
+assembly.
+
+.TP
 .BI -r \ INT
 Rounds of haplotype-aware error corrections [2]. This option affects all outputs of hifiasm.
 
 .SS Assembly options
 
-.TP 10
+.TP
 .BI -a \ INT
 Rounds of assembly graph cleaning [4]. This option is used with
-.I [-x maxDropRatio]
+.B -x
 and
-.I [-y minDropRatio].
+.BR -y .
 Note that unlike
-.I [-r],
+.BR -r ,
 this option does not affect error corrected reads and all-to-all overlaps.
 
-
-.TP 10
+.TP
 .BI -z \ INT
 Length of adapters that should be removed [0]. This option remove
 .I INT
@@ -100,40 +132,35 @@ Some old Hifi reads may consist of
 short adapters (e.g., 20bp adapter at one end). For such data, trimming short adapters would 
 significantly improve the assembly quality.
 
-
-.TP 10
+.TP
 .BI -m \ INT
-Maximal probing distance for bubble popping when generating primary/alternate assembly
+Maximal probing distance for bubble popping when generating primary/alternate
 contig graphs [10000000]. Bubbles longer than
 .I INT
-bases will not be popped. For detailed description of these graphs, please see
-.I 'Outputs'
+bases will not be popped. For detailed description of these graphs, please see the
+.B OUTPUTS
 section of this man-page. 
 
-
-.TP 10
+.TP
 .BI -p \ INT
 Maximal probing distance for bubble popping when generating haplotype-resolved processed unitig graph
 without small bubbles [100000]. Bubbles longer than
 .I INT
 bases will not be popped. Small bubbles might be caused by somatic mutations or noise in data, which
-are not the real haplotype information. For detailed description of this graph, please see
-.I 'Outputs'
+are not the real haplotype information. For detailed description of this graph, please see the
+.B OUTPUTS
 section of this man-page. 
 
-
-.TP 10
+.TP
 .BI -n \ INT
 A unitig is considered small if it is composed of less than 
 .I INT
 reads [3]. Hifiasm may try to remove small unitigs at various steps.
 
-
-
-.TP 10
+.TP
 .BI -x \ FLOAT, -y \ FLOAT
 Max and min overlap drop ratio [0.8, 0.2]. This option is used with
-.I [-r roundCorrection].
+.BR -r .
 Given a node
 .I N
 in the assembly graph, let max(N)
@@ -143,19 +170,21 @@ Hifiasm iteratively drops overlaps of
 .I N
 if their length / max(N)
 are below a threshold controlled by
-.I [-x maxDropRatio]
+.B -x
 and
-.I [-y minDropRatio].
+.BR -y .
 Hifiasm applies
-.I [-r roundCorrection]
+.B -r
 rounds of short overlap removal with an increasing threshold between
-.I [-x maxDropRatio]
+.B -x
 and
-.I [-y minDropRatio].
+.BR -y .
 
-.TP 10
+.TP
 .BI -i
-Ignore saved overlaps in [*.ovlp*] files.
+Ignore error corrected reads and overlaps saved in
+.IR prefix .*.bin
+files.
 Apart from assembly graphs, hifiasm also outputs three binary files
 that save all overlap information during assembly step.
 With these files, hifiasm can avoid the time-consuming all-to-all overlap calculation step,
@@ -167,139 +196,115 @@ with different parameters.
 .SS Trio-partition options
 
 .TP 10
-.BI -P \ FILE
-Paternal trio index. This index should be generated by
-.I [yak count]
-with the paternal short reads. For details of yak, please see
-.I [https://github.com/lh3/yak]
+.BI -1 \ FILE
+K-mer dump generated by
+.B yak count
+from the paternal/haplotype1 reads []
 
+.TP
+.BI -2 \ FILE
+K-mer dump generated by
+.B yak count
+from the maternal/haplotype2 reads []
 
-.TP 10
-.BI -M \ FILE
-Maternal trio index. This index should be generated by
-.I [yak count]
-with the maternal short reads. For details of yak, please see
-.I [https://github.com/lh3/yak]
+.TP
+.BI -3 \ FILE
+List of paternal/haplotype1 read names []
 
-.TP 10
+.TP
+.BI -4 \ FILE
+List of maternal/haplotype2 read names []
+
+.TP
 .BI -c \ INT
 Lower bound of the binned k-mer's frequency [2]. When doing trio binning, 
 a k-mer is said to be differentiating if it occurs >=
-.I [-d upper_bound] 
+.B -d
 times in one sample 
 but occurs <
-.I [-c lower_bound]  
+.B -c
 times in the other sample.
 
-.TP 10
+.TP
 .BI -d \ INT
 Upper bound of the binned k-mer's frequency [5]. When doing trio binning, 
 a k-mer is said to be differentiating if it occurs >=
-.I [-d upper_bound] 
+.B -d
 times in one sample 
 but occurs <
-.I [-c lower_bound]  
+.B -c
 times in the other sample.
 
 
+.SS Debugging options
 
-
-
-.SH EXAMPLES
-
-.TP
-.BR ./hifiasm " " \-o " " NA12878.asm " " \-t " " 32 " " NA12878_1.fq.gz " " NA12878_2.fq.gz
-In this example, hifiasm will be run with 32 CPU threads. The input read files are [NA12878_1.fq.gz]
-and [NA12878_2.fq.gz],
-while all output files can be found at [NA12878.asm.*].
-
-.TP
-.BR ./hifiasm " " \-o " " butterfly.asm " " \-t " " 32 " " \-z " " 20 " " butterfly.fq.gz
-In this example, hifiasm will be run with 32 CPU threads. The input read file is [butterfly.fq.gz],
-while all output files can be found at [butterfly.asm.*].
-With
-.I [-z 20],
-hifiasm will remove 20 bases from both ends of each read.
-
-.SH EXAMPLES FRO TRIO
-.TP
-.BR ./yak " " count " " \-k31 " " \-b37 " " \-t16 " " \-o " " mat.yak " " mat.fq.gz
-Build maternal trio index from mat.fq.gz.
-
-.TP
-.BR ./yak " " count " " \-k31 " " \-b37 " " \-t16 " " \-o " " pat.yak " " pat.fq.gz
-Build paternal trio index from pat.fq.gz.
-
-.TP
-.BR ./hifiasm " " \-o " " NA12878.asm " " \-t " " 32 " " \-P " " pat.yak " " \-M " " mat.yak " " NA12878_1.fq.gz " " NA12878_2.fq.gz
-In this example, hifiasm will do trio assembly with 32 CPU threads. The paternal assembly can be found at [NA12878.asm.p.r_utg.gfa],
-and the maternal assembly can be found at [NA12878.asm.m.r_utg.gfa].
-
+.TP 10
+.B --dbg-gfa
+Write additional files to speed up the debugging of graph cleaning
 
 
 .SH OUTPUTS
 
+.PP
+Without trio partition options
+.B -1
+and
+.BR -2 ,
+hifiasm generates the following assembly graphs in the GFA format:
+
+.RS 2
+.TP 2
+*
+.IR prefix .r_utg.gfa:
+haplotype-resolved raw unitig graph. This graph keeps all haplotype information.
+
+.TP
+*
+.IR prefix .p_utg.gfa:
+haplotype-resolved processed unitig graph without small bubbles.  Small bubbles
+might be caused by somatic mutations or noise in data, which are not the real
+haplotype information.  The size of popped small bubbles should be specified by
+.BR -p .
+
+.TP
+*
+.IR prefix .p_ctg.gfa:
+assembly graph of primary contigs. This graph collapses different haplotypes.
+
+.TP
+*
+.IR prefix .a_ctg.gfa:
+assembly graph of alternate contigs. This graph consists of all assemblies that
+are discarded in primary contig graph.
+
+.RE
 
 .PP
-Consider the prefix of output files has been specified by
-.I [-o outPrefix].
-During the error correction step, hifiasm outputs the following two files:
+With trio partition, hifiasm outputs the following assembly graphs:
 
-.IP
-1. Haplotype-aware error corrected reads in fasta format [outPrefix.ec.fa].
+.RS 2
+.TP 2
+*
+.IR prefix .dip.r_utg.gfa:
+haplotype-resolved raw unitig graph. This graph keeps all haplotype information.
 
-2. All-to-all overlaps in paf format [outPrefix.ovlp.paf].
+.TP
+*
+.IR prefix .hap1.p_ctg.gfa:
+phased paternal/haplotype1 contig graph. This graph keeps the phased
+paternal/haplotype1 assembly.
+
+.TP
+*
+.IR prefix .hap2.p_ctg.gfa:
+phased maternal/haplotype2 contig graph. This graph keeps the phased
+maternal/haplotype2 assembly.
+.RE
 
 .PP
-During the non-trio assembly step, hifiasm outputs the following four assembly graphs in GFA format:
+For each graph, hifiasm also outputs a simplified version without sequences for
+the ease of visualization. Hifiasm keeps corrected reads and overlaps in three
+binary files such as it can regenerate assembly graphs from the binary files
+without redoing error correction.
 
 
-.IP
-1. Haplotype-resolved raw unitig graph [outPrefix.r_utg.gfa].
-This graph keeps all haplotype information.
-
-
-2. Haplotype-resolved processed unitig graph without small bubbles [outPrefix.p_utg.gfa].
-Small bubbles might be caused by somatic mutations or noise in data, which are not the real haplotype information.
-The size of popped small bubbles should be specified by
-.I [-p maxSmallBubbles].
-
-
-3. Primary assembly contig graph [outPrefix.p_ctg.gfa].
-This graph collapses different haplotypes.
-
-4. Alternate assembly contig graph [outPrefix.a_ctg.gfa].
-This graph consists of all assemblies that are discarded in primary assembly contig graph.
-
-
-.PP
-If you have trio information, hifiasm outputs the following three assembly graphs in GFA format:
-
-.IP
-1. Phased maternal unitig graph [outPrefix.m.r_utg.gfa].
-This graph keeps the phased maternal assembly.
-
-2. Phased paternal unitig graph [outPrefix.p.r_utg.gfa].
-This graph keeps the phased paternal assembly.
-
-3. Haplotype-resolved raw unitig graph [outPrefix.r_utg.gfa].
-This graph keeps all haplotype information.
-
-
-.PP
-For each graph, hifiasm also outputs a simplified version without sequences. These simplified
-graphs can be easily visualized.
-
-.PP
-Note that different species need different assembly graphs. For homozygous genomes,
-the primary assembly contig graph is the best choice. 
-For species with high heterozygous rate, different haplotypes can be fully separated.
-It is important to remove small bubbles from the haplotype-resolved unitig graph. The
-reason is that some small bubbles are caused by somatic mutations or noise in data,
-which are not the real haplotype information. In this case, haplotype-resolved processed
-unitig graph without small bubbles should be better.
-For ordinary human genome, different haplotypes cannot be fully separated due to the low
-heterozygous rate. There are many small bubbles including haplotype information,
-which cannot be simply removed. Thus, it is necessary to use the haplotype-resolved raw
-unitig graph.
-

htab.cpp

@@ -482,6 +482,7 @@ KSEQ_INIT(gzFile, gzread)
 #define HAF_RS_WRITE_LEN 0x4
 #define HAF_RS_WRITE_SEQ 0x8
 #define HAF_RS_READ      0x10
+#define HAF_CREATE_NEW   0x20
 
 typedef struct { // global data structure for kt_pipeline()
 	const yak_copt_t *opt;
@@ -665,12 +666,13 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 	return 0;
 }
 
-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)
+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)
 {
 	int read_rs = (rs && (flag & HAF_RS_READ));
 	pl_data_t pl;
 	gzFile fp = 0;
 	memset(&pl, 0, sizeof(pl_data_t));
+	pl.n_seq = *n_seq;
 	if (read_rs) {
 		pl.rs_in = rs;
 		init_UC_Read(&pl.ucr);
@@ -684,13 +686,13 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 	pl.opt = opt;
 	pl.flag = flag;
 	if (p0) {
-		pl.pt = p0, pl.create_new = 0;
+		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 = 0;
+		pl.ct = c0, pl.create_new = !!(flag&HAF_CREATE_NEW);
 		assert(c0->k == opt->k && c0->pre == opt->pre);
 	} else {
-		pl.create_new = 1;
+		pl.create_new = 1; // alware create new elements if the count table is empty
 		pl.ct = ha_ct_init(opt->k, opt->pre, opt->bf_n_hash, opt->bf_shift);
 	}
 	kt_pipeline(3, worker_count, &pl, 3);
@@ -700,12 +702,14 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 		kseq_destroy(pl.ks);
 		gzclose(fp);
 	}
+	*n_seq = pl.n_seq;
 	return pl.ct;
 }
 
 ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const void *flt_tab, All_reads *rs)
 {
 	int i;
+	int64_t n_seq = 0;
 	yak_copt_t opt;
 	ha_ct_t *h = 0;
 	assert(!(flag & HAF_RS_WRITE_LEN) || !(flag & HAF_RS_WRITE_SEQ)); // not both
@@ -717,12 +721,12 @@ 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;
-	opt.is_HPC = !asm_opt->no_HPC;
+	opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC);
 	opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win;
 	opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift;
 	opt.n_thread = asm_opt->thread_num;
 	for (i = 0; i < asm_opt->num_reads; ++i)
-		h = yak_count(&opt, asm_opt->read_file_names[i], flag, p0, h, flt_tab, rs);
+		h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, &n_seq);
 	if (h && opt.bf_shift > 0)
 		ha_ct_destroy_bf(h);
 	return h;
@@ -771,7 +775,7 @@ void ha_ft_destroy(void *h)
  * High-level interfaces *
  *************************/
 
-void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs)
+void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov)
 {
 	yak_ft_t *flt_tab;
 	int64_t cnt[YAK_N_COUNTS];
@@ -780,6 +784,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs)
 	h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN, NULL, NULL, rs);
 	ha_ct_hist(h, cnt, asm_opt->thread_num);
 	peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
+	if (hom_cov) *hom_cov = peak_hom;
 	if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
 	cutoff = (int)(peak_hom * asm_opt->high_factor);
 	if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
@@ -791,7 +796,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs)
 	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)
+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)
 {
 	int64_t cnt[YAK_N_COUNTS], tot_cnt;
 	int peak_hom, peak_het, i, extra_flag1, extra_flag2;
@@ -812,6 +817,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 	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);
+	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);
 	if (flt_tab == 0) {
 		int cutoff = (int)(peak_hom * asm_opt->high_factor);

htab.h

@@ -26,11 +26,11 @@ extern const unsigned char seq_nt4_table[256];
 extern void *ha_flt_tab;
 extern ha_pt_t *ha_idx;
 
-void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs);
+void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov);
 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);
+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);
 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);
 
@@ -45,7 +45,7 @@ long yak_peakrss(void);
 double yak_peakrss_in_gb(void);
 double yak_cpu_usage(void);
 
-void trio_partition(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);

ketopt.h

@@ -17,7 +17,7 @@ typedef struct {
 } ketopt_t;
 
 typedef struct {
-	char *name;
+	const char *name;
 	int has_arg;
 	int val;
 } ko_longopt_t;