r257: dropped contained unitigs

Assembly.cpp

@@ -10,6 +10,7 @@
 #include "Correct.h"
 #include "htab.h"
 #include "kthread.h"
+#include "utils.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);
 
@@ -1159,11 +1160,17 @@ void ha_overlap_final(void)
 
 int ha_assemble(void)
 {
+	extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
 	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 corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
+		if (asm_opt.extract_list) {
+			ha_extract_print_list(&R_INF, asm_opt.extract_iter, asm_opt.extract_list);
+			exit(0);
+		}
 		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_SKIP_TRIOBIN)) 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();
 	}

CommandLines.cpp

@@ -10,11 +10,16 @@
 hifiasm_opt_t asm_opt;
 
 static ko_longopt_t long_options[] = {
-	{ "version",      ko_no_argument,       300 },
+	{ "version",       ko_no_argument, 300 },
-	{ "dbg-gfa",      ko_no_argument,       301 },
+	{ "dbg-gfa",       ko_no_argument, 301 },
-	{ "write-paf",    ko_no_argument,       302 },
+	{ "write-paf",     ko_no_argument, 302 },
-	{ "write-ec",     ko_no_argument,       303 },
+	{ "write-ec",      ko_no_argument, 303 },
-	{ "skip-triobin", ko_no_argument,       304 },
+	{ "skip-triobin",  ko_no_argument, 304 },
+	{ "max-od-ec",     ko_no_argument, 305 },
+	{ "max-od-final",  ko_no_argument, 306 },
+	{ "ex-list",       ko_required_argument, 307 },
+	{ "ex-iter",       ko_required_argument, 308 },
+	{ "keep-contained",ko_no_argument, 309 },
 	{ 0, 0, 0 }
 };
 
@@ -83,7 +88,9 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->k_mer_length = 51;
 	asm_opt->mz_win = 51;
 	asm_opt->bf_shift = 37;
-	asm_opt->high_factor = 5.0f;
+	asm_opt->high_factor = 5.0;
+	asm_opt->max_ov_diff_ec = 0.04;
+	asm_opt->max_ov_diff_final = 0.03;
 	asm_opt->hom_cov = 20;
 	asm_opt->max_n_chain = 100;
     asm_opt->k_mer_min_freq = 3;
@@ -264,6 +271,15 @@ int check_option(hifiasm_opt_t* asm_opt)
         return 0;
     }
 
+	if (asm_opt->max_ov_diff_ec < asm_opt->max_ov_diff_final) {
+		fprintf(stderr, "[ERROR] max_ov_diff_ec shouldn't be smaller than max_ov_diff_final\n");
+		return 0;
+	}
+
+	if (asm_opt->max_ov_diff_ec < HA_MIN_OV_DIFF) {
+		fprintf(stderr, "[ERROR] max_ov_diff_ec shouldn't be smaller than %g\n", HA_MIN_OV_DIFF);
+		return 0;
+	}
 
     if(asm_opt->max_short_tip < 0)
     {
@@ -385,6 +401,11 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 		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 == 305) asm_opt->max_ov_diff_ec = atof(opt.arg);
+		else if (c == 306) asm_opt->max_ov_diff_final = atof(opt.arg);
+		else if (c == 307) asm_opt->extract_list = opt.arg;
+		else if (c == 308) asm_opt->extract_iter = atoi(opt.arg);
+		else if (c == 309) asm_opt->flag |= HA_F_KEEP_CONTAINED;
         else if (c == 'l')
         {   ///0: disable purge_dup; 1: purge containment; 2: purge overlap
             asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);

CommandLines.h

@@ -3,18 +3,21 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.5-dirty-r247"
+#define HA_VERSION "0.7-dirty-r257"
 
 #define VERBOSE 0
 
-#define HA_F_NO_HPC       0x1
+#define HA_F_NO_HPC          0x1
-#define HA_F_NO_KMER_FLT  0x2
+#define HA_F_NO_KMER_FLT     0x2
-#define HA_F_VERBOSE_GFA  0x4
+#define HA_F_VERBOSE_GFA     0x4
-#define HA_F_WRITE_EC     0x8
+#define HA_F_WRITE_EC        0x8
-#define HA_F_WRITE_PAF    0x10
+#define HA_F_WRITE_PAF       0x10
-#define HA_F_SKIP_TRIOBIN 0x20
+#define HA_F_SKIP_TRIOBIN    0x20
 #define HA_F_PURGE_CONTAIN   0x40
-#define HA_F_PURGE_JOIN   0x80
+#define HA_F_PURGE_JOIN      0x80
+#define HA_F_KEEP_CONTAINED  0x100
+
+#define HA_MIN_OV_DIFF       0.02 // min sequence divergence in an overlap
 
 typedef struct {
 	int flag;
@@ -24,11 +27,15 @@ typedef struct {
     char* required_read_name;
 	char *fn_bin_yak[2];
 	char *fn_bin_list[2];
+	char *extract_list;
+	int extract_iter;
     int thread_num;
     int k_mer_length;
 	int mz_win;
 	int bf_shift;
-	float high_factor; // coverage cutoff set to high_factor*hom_cov
+	double high_factor; // coverage cutoff set to high_factor*hom_cov
+	double max_ov_diff_ec;
+	double max_ov_diff_final;
 	int hom_cov;
 	int max_n_chain; // fall-back max number of chains to consider
     int k_mer_min_freq;

Correct.cpp

@@ -260,7 +260,7 @@ All_reads* R_INF)
     
 
     x_len = x_end - x_start + 1;
-    threshold = x_len * THRESHOLD_RATE;
+    threshold = x_len * asm_opt.max_ov_diff_ec;
     /****************************may have bugs********************************/
     threshold = Adjust_Threshold(threshold, x_len);
     /****************************may have bugs********************************/
@@ -483,7 +483,7 @@ char* r_string)
 
         ///overlap length between [window_start, window_end]
         x_len = x_end - x_start + 1;
-        threshold = x_len * THRESHOLD_RATE;
+        threshold = x_len * asm_opt.max_ov_diff_ec;
         /****************************may have bugs********************************/
         threshold = Adjust_Threshold(threshold, x_len);
         /****************************may have bugs********************************/
@@ -548,8 +548,8 @@ inline double trim_error_rate(overlap_region_alloc* overlap_list, long long ID)
             }
             else
             {
-                ///tError += (Adjust_Threshold(subWinLen*THRESHOLD_RATE, subWinLen) * 2);
+                ///tError += (Adjust_Threshold(subWinLen*asm_opt.max_ov_diff_ec, subWinLen) * 2);
-                tError += (Adjust_Threshold(subWinLen*THRESHOLD_RATE, subWinLen) * 3);
+                tError += Adjust_Threshold(subWinLen * asm_opt.max_ov_diff_ec, subWinLen) * 3;
             }
         }
     }
@@ -566,8 +566,8 @@ inline double trim_error_rate(overlap_region_alloc* overlap_list, long long ID)
             }
             else
             {
-                ///tError += (Adjust_Threshold(subWinLen*THRESHOLD_RATE, subWinLen) * 2);
+                ///tError += (Adjust_Threshold(subWinLen*asm_opt.max_ov_diff_ec, subWinLen) * 2);
-                tError += (Adjust_Threshold(subWinLen*THRESHOLD_RATE, subWinLen) * 3);
+                tError += Adjust_Threshold(subWinLen * asm_opt.max_ov_diff_ec, subWinLen) * 3;
             }
         }
     }
@@ -2118,7 +2118,7 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
                         x_end = overlap_list->list[j].w_list[i].x_end;
                         x_len = x_end - x_start + 1;
                         /****************************may have bugs********************************/
-                        ///threshold = x_len * THRESHOLD_RATE;
+                        ///threshold = x_len * asm_opt.max_ov_diff_ec;
                         threshold = overlap_list->list[j].w_list[i].error_threshold;
                         /****************************may have bugs********************************/
                         /****************************may have bugs********************************/
@@ -2415,7 +2415,7 @@ UC_Read* g_read)
         // }
 
 
-        threshold = xLen * THRESHOLD_RATE;
+        threshold = xLen * asm_opt.max_ov_diff_ec;
         threshold = Adjust_Threshold(threshold, xLen);
         threshold = double_error_threshold(threshold, xLen);
 
@@ -2931,7 +2931,7 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea
                         x_end = overlap_list->list[j].w_list[i].x_end;
                         x_len = x_end - x_start + 1;
                         /****************************may have bugs********************************/
-                        ///threshold = x_len * THRESHOLD_RATE;
+                        ///threshold = x_len * asm_opt.max_ov_diff_ec;
                         threshold = overlap_list->list[j].w_list[i].error_threshold;
                         /****************************may have bugs********************************/
                         /****************************may have bugs********************************/
@@ -2998,19 +2998,16 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea
                 }
             }
 
-            ///error_rate = trim_error_rate(overlap_list, j);
             error_rate = non_trim_error_rate(overlap_list, j, R_INF, dumy, g_read);
             
             
-            ///if(error_rate <= 0.015)
+            if (error_rate <= asm_opt.max_ov_diff_final)
-            ///if(error_rate <= 0.03)
-            if(error_rate <= FINAL_OVERLAP_ERROR_RATE)
             {
                 overlap_list->mapped_overlaps_length += overlap_length;
                 overlap_list->list[j].is_match = 1;
                 calculate_boundary_cigars(&(overlap_list->list[j]), R_INF, dumy, g_read);
             }
-            else if(error_rate <= 0.045)
+            else if (error_rate <= asm_opt.max_ov_diff_final * 1.5)
             {
                 overlap_list->list[j].is_match = 3;
             }
@@ -4552,7 +4549,7 @@ Cigar_record* current_cigar, long long uncorrected_window_start, Round2_alignmen
             x_start = corrected_window_start;
             x_end = corrected_window_end;
             x_len = x_end - x_start + 1;
-            threshold = x_len * THRESHOLD_RATE;
+            threshold = x_len * asm_opt.max_ov_diff_ec;
             /****************************may have bugs********************************/
             threshold = Adjust_Threshold(threshold, x_len);
             /****************************may have bugs********************************/
@@ -4741,7 +4738,7 @@ void generate_consensus(overlap_region_alloc* overlap_list, All_reads* R_INF,
 
 
     Window_Pool w_inf;
-    init_Window_Pool(&w_inf, g_read->length, WINDOW, TAIL_LENGTH);
+    init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec));
 
     int flag = 0;
     ///for last window
@@ -4885,7 +4882,7 @@ All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy, Graph* g, int* abnormal)
     (*abnormal) = 0;
     
     Window_Pool w_inf;
-    init_Window_Pool(&w_inf, g_read->length, WINDOW, TAIL_LENGTH);
+    init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec));
 
     int flag = 0;
     long long realLen = 0, tmpLen = 0;
@@ -6594,7 +6591,7 @@ Correct_dumy* dumy)
     {
         long long window_start, window_end;
         Window_Pool w_inf;
-        init_Window_Pool(&w_inf, read_length, WINDOW, TAIL_LENGTH);
+        init_Window_Pool(&w_inf, read_length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec));
         int flag = 0;
         long long realLen = 0, realLen_100 = 0;
         int to_recover = 0;
@@ -7046,7 +7043,7 @@ void partition_overlaps(overlap_region_alloc* overlap_list, All_reads* R_INF,
     long long num_availiable_win = 0;
     
     Window_Pool w_inf;
-    init_Window_Pool(&w_inf, g_read->length, WINDOW, TAIL_LENGTH);
+    init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec));
 
     int flag = 0;
     while(get_Window(&w_inf, &window_start, &window_end) && flag != -2)
@@ -7135,7 +7132,7 @@ void partition_overlaps_advance(overlap_region_alloc* overlap_list, All_reads* R
     long long num_availiable_win = 0;
     
     Window_Pool w_inf;
-    init_Window_Pool(&w_inf, g_read->length, WINDOW, TAIL_LENGTH);
+    init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec));
 
     int flag = 0;
     while(get_Window(&w_inf, &window_start, &window_end) && flag != -2)
@@ -7219,7 +7216,7 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
 
     Window_Pool w_inf;
 
-    init_Window_Pool(&w_inf, g_read->length, WINDOW, TAIL_LENGTH);
+    init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec));
 
     int flag = 0;
 

Correct.h

@@ -16,6 +16,7 @@
 #define INSERTION 2
 #define DELETION 3
 
+#define WINDOW_MAX_SIZE (WINDOW + (int)(1.0 / HA_MIN_OV_DIFF) + 3) // TODO: why 1/max_ov_diff?
 
 ///#define FLAG_THRE 0
 

Hash_Table.cpp

@@ -4,6 +4,7 @@
 #include <pthread.h>
 #include "Hash_Table.h"
 #include "ksort.h"
+#include "utils.h"
 pthread_mutex_t output_mutex;
 
 #define overlap_region_key(a) ((a).y_id)

Hash_Table.h

@@ -6,22 +6,13 @@
 #define MAX_SUFFIX_BITS 64
 #define MODE_VALUE 101
 
-///#define WINDOW 350
-///#define THRESHOLD  14
-
 #define WINDOW 375
-//#define WINDOW_BOUNDARY 150
 #define WINDOW_BOUNDARY 375
 ///for one side, the first or last WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY bases should not be corrected
 #define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25
 #define THRESHOLD  15
-#define THRESHOLD_RATE 0.04
-#define TAIL_LENGTH int(1/THRESHOLD_RATE)
-///#define OVERLAP_THRESHOLD 0.9
 #define OVERLAP_THRESHOLD_FILTER 0.9
-#define WINDOW_MAX_SIZE WINDOW + TAIL_LENGTH + 3
 #define THRESHOLD_MAX_SIZE  31
-#define FINAL_OVERLAP_ERROR_RATE 0.03
 
 #define GROUP_SIZE 4
 ///the max cigar likes 10M10D10M10D10M

Makefile

@@ -4,7 +4,7 @@ CPPFLAGS=
 INCLUDES=
 OBJS=		CommandLines.o Process_Read.o Assembly.o Hash_Table.o \
 			POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \
-			htab.o hist.o sketch.o anchor.o sys.o
+			htab.o hist.o sketch.o anchor.o extract.o overlap2.o sys.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm
 
@@ -34,32 +34,35 @@ 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: kthread.h
+Assembly.o: kthread.h utils.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 ksort.h
+Hash_Table.o: CommandLines.h ksort.h utils.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
+Overlaps.o: Purge_Dups.h utils.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: khashl.h kthread.h kseq.h Process_Read.h Overlaps.h kvec.h kdq.h
-Trio.o: CommandLines.h htab.h
+Trio.o: CommandLines.h htab.h utils.h
 anchor.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
-anchor.o: ksort.h Hash_Table.h
+anchor.o: ksort.h utils.h Hash_Table.h
+extract.o: Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h khashl.h
+extract.o: kseq.h utils.h
 hist.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
 htab.o: kthread.h khashl.h kseq.h ksort.h htab.h Process_Read.h Overlaps.h
-htab.o: kvec.h kdq.h CommandLines.h
+htab.o: kvec.h kdq.h CommandLines.h utils.h
 kthread.o: kthread.h
 main.o: CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h Assembly.h
-main.o: Levenshtein_distance.h htab.h
+main.o: Levenshtein_distance.h utils.h
+overlap2.o: utils.h CommandLines.h Overlaps.h kvec.h kdq.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
+sys.o: utils.h

Overlaps.cpp

@@ -9,6 +9,7 @@
 #include "Hash_Table.h"
 #include "Correct.h"
 #include "Purge_Dups.h"
+#include "utils.h"
 
 uint32_t debug_purge_dup = 0;
 
@@ -76,21 +77,17 @@ void asg_destroy(asg_t *g)
 {
 	if (g == 0) return;
 	free(g->seq); free(g->idx); free(g->arc); free(g->seq_vis);
+	free(g->contain); free(g->con_idx);
 
-    if(g->n_F_seq > 0 && g->F_seq)
+	if (g->n_F_seq > 0 && g->F_seq) {
-    {
+		uint32_t i = 0;
-        uint32_t i = 0;
+		for (i = 0; i < g->n_F_seq; ++i) {
-        for (i = 0; i < g->n_F_seq; i++)
+			if (g->F_seq[i].a) free(g->F_seq[i].a);
-        {
+			if (g->F_seq[i].s) free(g->F_seq[i].s);
-            if(g->F_seq[i].a) free(g->F_seq[i].a);
+		}
-            if(g->F_seq[i].s) free(g->F_seq[i].s);
+		free(g->F_seq);
-        }
+	}
-
+	free(g);
-        free(g->F_seq);
-    }
-    
-
-    free(g);
 }
 
 void asg_arc_sort(asg_t *g)
@@ -132,8 +129,7 @@ void remove_overlaps(ma_hit_t_alloc* source_paf, uint64_t* source_index, long lo
 }
 
 
-void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, 
+void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, uint64_t* source_index, long long listLen)
-uint64_t* source_index, long long listLen)
 {
     long long i;
     ma_hit_t ele;
@@ -186,7 +182,6 @@ uint64_t *asg_arc_index_core(size_t max_seq, size_t n, const asg_arc_t *a)
 	uint64_t *idx;
 	idx = (uint64_t*)calloc(max_seq * 2, 8);
 
-
     /**
 		 * ul: |____________31__________|__________1___________|______________32_____________|
 	                       qns            direction of overlap       length of this node (not overlap length)
@@ -196,7 +191,6 @@ uint64_t *asg_arc_index_core(size_t max_seq, size_t n, const asg_arc_t *a)
 		if (i == n || a[i-1].ul>>32 != a[i].ul>>32)
 			idx[a[i-1].ul>>32] = (uint64_t)last<<32 | (i - last), last = i;
 
-    
 	return idx;
 }
 
@@ -276,16 +270,15 @@ void asg_cleanup(asg_t *g)
 	                       qns            direction of overlap       length of this node (not overlap length)
 		**/
 		asg_arc_sort(g);
+		asg_con_sort(g);
 		g->is_srt = 1;
 	}
 	///index the overlaps in graph with query id
 	if (g->idx == 0) asg_arc_index(g);
+	if (g->con_idx == 0 && g->n_con > 0) asg_con_index(g);
 }
 
 
-
-
-
 // delete multi-arcs
 /**
  * remove edges like:   v has two out-edges to w
@@ -872,51 +865,6 @@ void drop_edges_by_trio(ma_hit_t_alloc* sources, long long num_sources)
 }
 
 
-
-ma_hit_t* get_specific_overlap_with_del(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, 
-uint32_t qn, uint32_t tn)
-{
-    if(coverage_cut[qn].del || coverage_cut[tn].del) return NULL;
-    ma_hit_t_alloc* x = &(sources[qn]);
-    uint32_t i;
-    for (i = 0; i < x->length; i++)
-    {
-        if(x->buffer[i].del) continue;
-        if(coverage_cut[Get_qn(x->buffer[i])].del) continue;
-        if(coverage_cut[Get_tn(x->buffer[i])].del) continue;
-
-        if(Get_tn(x->buffer[i])==tn 
-        &&
-           Get_qn(x->buffer[i])==qn)
-        {
-            return &(x->buffer[i]);
-        }        
-    }
-
-    return NULL;
-}
-
-
-
-void delete_single_edge(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, uint32_t qn, uint32_t tn)
-{
-    ma_hit_t* tmp = get_specific_overlap_with_del(sources, coverage_cut, qn, tn);
-    if(tmp != NULL) tmp->del = 1;
-}
-
-void delete_all_edges(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, uint32_t qn)
-{
-    ma_hit_t_alloc* x = &(sources[qn]);
-    uint32_t i;
-    for (i = 0; i < x->length; i++)
-    {
-        x->buffer[i].del = 1;
-        delete_single_edge(sources, coverage_cut, Get_tn(x->buffer[i]), Get_qn(x->buffer[i]));
-    }
-    coverage_cut[qn].del = 1;
-}
-
-
 uint32_t get_real_sources_length(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, 
 int max_hang, int min_ovlp, uint32_t query)
 {
@@ -1027,108 +975,6 @@ int max_hang, int min_ovlp, uint32_t qn)
     return keep_edge;
 }
 
-
-void ma_hit_contained_advance(ma_hit_t_alloc* sources, long long n_read, ma_sub_t *coverage_cut, 
-R_to_U* ruIndex, int max_hang, int min_ovlp)
-{
-    ///uint32_t qn_num = 0, no_fully_qn_num = 0, tn_num = 0, no_fully_tn_num = 0;
-    double startTime = Get_T();
-	int32_t r;
-	long long i, j, m;
-	asg_arc_t t;
-    ma_hit_t *h = NULL;
-    ma_sub_t *sq = NULL;
-    ma_sub_t *st = NULL;
-
-
-    for (i = 0; i < n_read; ++i) 
-    {
-        if(coverage_cut[i].del) continue;
-
-        for (j = 0; j < (long long)sources[i].length; j++)
-        {
-            h = &(sources[i].buffer[j]);
-            //check the corresponding two reads 
-		    sq = &(coverage_cut[Get_qn(*h)]);
-            st = &(coverage_cut[Get_tn(*h)]);
-            /****************************may have trio bugs********************************/
-            if(sq->del || st->del) continue;
-            if(h->del) continue;
-            /****************************may have trio bugs********************************/
-            r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
-            ///r could not be MA_HT_SHORT_OVLP or MA_HT_INT
-            if (r == MA_HT_QCONT) 
-            {
-                h->del = 1;
-                delete_single_edge(sources, coverage_cut, Get_tn(*h), Get_qn(*h));
-        
-                delete_all_edges(sources, coverage_cut, Get_qn(*h));
-                set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0);
-
-                // if(delete_all_edges_carefully(sources, coverage_cut, max_hang, min_ovlp, 
-                // Get_qn(*h))==0)
-                // {
-                //     set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0);
-                // }
-                // sq->del = 1;
-                // set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0);
-            }
-		    else if (r == MA_HT_TCONT) 
-            {
-                h->del = 1;
-                delete_single_edge(sources, coverage_cut, Get_tn(*h), Get_qn(*h));
-
-                delete_all_edges(sources, coverage_cut, Get_tn(*h));
-                set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0);
-
-                // if(delete_all_edges_carefully(sources, coverage_cut, max_hang, 
-                // min_ovlp, Get_tn(*h)) == 0)
-                // {
-                //     set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0);
-                //     no_fully_tn_num++;
-                // }
-                // st->del = 1;
-                // set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0);
-            }
-        }
-    }
-
-    transfor_R_to_U(ruIndex);
-
-
-    for (i = 0; i < n_read; ++i) 
-    {
-        m = 0;
-        for (j = 0; j < (long long)sources[i].length; j++)
-        {
-            ma_hit_t *h = &(sources[i].buffer[j]);
-            if(h->del) continue;
-            ///both the qn and tn have not been deleted
-            if(coverage_cut[Get_qn(*h)].del != 1 && coverage_cut[Get_tn(*h)].del != 1)
-            {
-                h->del = 0;
-                m++;
-            }
-            else
-            {
-                h->del = 1;
-            }
-        }
-
-        ///if sources[i].length == 0, that means all overlapped reads with read i are the contained reads
-        if(m == 0)
-        {
-            coverage_cut[i].del = 1;
-        }
-    }
-
-    if(VERBOSE >= 1)
-    {
-        fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime);
-    }
-}
-
-
 void ma_hit_flt(ma_hit_t_alloc* sources, long long n_read, ma_sub_t *coverage_cut, int max_hang, int min_ovlp)
 {
     double startTime = Get_T();
@@ -1694,9 +1540,7 @@ void print_overlaps(ma_hit_t_alloc* paf, long long rLen, long long interval_s, l
 
 }
 
-
+void detect_chimeric_reads(ma_hit_t_alloc* paf, long long n_read, uint64_t* readLen, ma_sub_t* coverage_cut, float shift_rate)
-void detect_chimeric_reads(ma_hit_t_alloc* paf, long long n_read, uint64_t* readLen, 
-ma_sub_t* coverage_cut, float shift_rate)
 {
     double startTime = Get_T();
     init_aux_table();
@@ -1725,10 +1569,9 @@ ma_sub_t* coverage_cut, float shift_rate)
         collect_contain(&(paf[i]), NULL, rLen, &max_left, &max_right, 0.1);
         ///collect_contain(&(paf[i]), &(rev_paf[i]), rLen, &max_left, &max_right, 0.1);
 
-        ////shift_rate should be (FINAL_OVERLAP_ERROR_RATE*2)
+        ////shift_rate should be (asm_opt.max_ov_diff_final*2)
         ///this read is a normal read
-        if(max_left.e > max_right.s && 
+        if (max_left.e > max_right.s && (max_left.e - max_right.s >= rLen * shift_rate))
-           (max_left.e - max_right.s >= rLen * shift_rate))
         {
             continue;
         }
@@ -2109,71 +1952,6 @@ static inline int asg_is_single_edge(const asg_t *g, uint32_t v, uint32_t start_
 }
 
 
-asg_t *ma_sg_gen(const ma_hit_t_alloc* sources, long long n_read, const ma_sub_t *coverage_cut, 
-int max_hang, int min_ovlp)
-{
-    double startTime = Get_T();
-	size_t i, j;
-	asg_t *g;
-	///just calloc
-	g = asg_init();
-
-	///add seq to graph, seq just save the length of each read
-	for (i = 0; i < (uint64_t)n_read; ++i) 
-    {
-        ///if a read has been deleted, should we still add them?
-		asg_seq_set(g, i, coverage_cut[i].e - coverage_cut[i].s, coverage_cut[i].del);
-        g->seq[i].c = coverage_cut[i].c;
-	}
-
-    g->seq_vis = (uint8_t*)calloc(g->n_seq*2, sizeof(uint8_t));
-
-    for (i = 0; i < (uint64_t)n_read; ++i) 
-    {
-        for (j = 0; j < sources[i].length; j++)
-        {
-            int r;
-		    asg_arc_t t, *p;
-		    const ma_hit_t *h = &(sources[i].buffer[j]);
-            if(h->del) continue;
-
-            //high coverage region [sub[qn].e, sub[qn].s) in query
-            int ql = coverage_cut[Get_qn(*h)].e - coverage_cut[Get_qn(*h)].s;
-            //high coverage region [sub[qn].e, sub[qn].s) in target
-            int tl = coverage_cut[Get_tn(*h)].e - coverage_cut[Get_tn(*h)].s;
-		    r = ma_hit2arc(h, ql, tl, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
-            /**
-            #define MA_HT_INT       (-1)
-            #define MA_HT_QCONT      (-2)
-            #define MA_HT_TCONT      (-3)
-            #define MA_HT_SHORT_OVLP (-4)
-            the short overlaps and the overlaps with contain reads have already been removed
-            here we should have overhang
-            so r should always >= 0
-            **/
-            if (r >= 0) 
-            {
-                ///push node?
-                p = asg_arc_pushp(g);
-                *p = t;
-		    }
-            else
-            {
-                fprintf(stderr, "error\n");
-            }
-        }
-    }
-
-    asg_cleanup(g);
-    g->r_seq = g->n_seq;
-
-    if(VERBOSE >= 1)
-    {
-        fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime);
-    }
-
-	return g;
-}
 
 
 
@@ -4485,143 +4263,6 @@ int asg_arc_del_cross_bubble(asg_t *g, long long max_dist)
     return n_reduced;
 }
 
-
-
-// transitive reduction; see Myers, 2005
-int asg_arc_del_trans(asg_t *g, int fuzz)
-{
-    double startTime = Get_T();
-
-	uint8_t *mark;
-	///n_vtx = number of seq * 2
-	///the reason is that each read has two direction (query->target, target->query)
-	uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0;
-	///at first, all nodes should be set to vacant
-	mark = (uint8_t*)calloc(n_vtx, 1);
-
-	/**v is the id+direction of a node, 
-     * the high 31-bit is the id, 
-     * and the lowest 1-bit is the direction
-     * (0 means query-to-target, 1 means target-to-query)**/
-	for (v = 0; v < n_vtx; ++v) {
-		///nv is the number of overlaps with v(qn+direction)
-		uint32_t L, i, nv = asg_arc_n(g, v);
-		///av is the array of v
-		asg_arc_t *av = asg_arc_a(g, v);
-        ///that means in this direction, read v is not overlapped with any other reads
-		if (nv == 0) continue; // no hits
-		
-        ///if the read itself has been removed
-		if (g->seq[v>>1].del) 
-        {
-			for (i = 0; i < nv; ++i) av[i].del = 1, ++n_reduced;
-			continue;
-		}
-
-
-
-        /**
-	********************************query-to-target overlap****************************
-	case 1: u = 0, rev = 0                           in the view of target: direction is 1 
-	query: CCCCCCCCTAATTAAAAT                        target: TAATTAAAATGGGGGG (use ex-target as query)
-	               ||||||||||         <--->                  ||||||||||
-	       target: TAATTAAAATGGGGGG           query: CCCCCCCCTAATTAAAAT (use ex-query as target)
-
-	case 2: u = 0, rev = 1                           in the view of target: direction is 0
-	query: CCCCCCCCTAATTAAAAT					     target: CCCCCCATTTTAATTA  (use ex-target as query)
-                   ||||||||||        <--->                         ||||||||||
-	       target: TAATTAAAATGGGGGG                         query: ATTTTAATTAGGGGGGGG  (use ex-query as target)
-	********************************query-to-target overlap****************************
-
-	********************************target-to-query overlap****************************
-	case 3: u = 1, rev = 0                           in the view of target: direction is 0
-			 query: AAATAATATCCCCCCGCG                target: GGGCCGGCAAATAATAT (use ex-target as query)
-					|||||||||          <--->                          |||||||||
-	target: GGGCCGGCAAATAATAT                                  query: AAATAATATCCCCCCGCG (use ex-query as target)
-
-	case 4: u = 1, rev = 1                          in the view of target: direction is 1
-	         query: AAATAATATCCCCCCGCG                        target: ATATTATTTGCCGGCCC (use ex-target as query)
-                    |||||||||             <--->                       |||||||||
-	target: GGGCCGGCAAATAATAT                          query: CGCGGGGGATATTATTT (use ex-query as target)
-	********************************target-to-query overlap****************************
-
-	p->ul: |____________31__________|__________1___________|______________32_____________|
-	                    qns            direction of overlap       length of this node (not overlap length)
-						                (in the view of query)
-	p->v : |___________31___________|__________1___________|
-				        tns             reverse direction of overlap
-						              (in the view of target)
-	p->ol: overlap length
-    **/
-
-
-		//all outnode of v should be set to "not reduce"
-		for (i = 0; i < nv; ++i) mark[av[i].v] = 1;
-
-		///length of node (not overlap length)
-		///av[nv-1] is longest out-dege
-		/**
-		 * v---------------
-		 *   w1---------------
-		 *      w2--------------
-		 *         w3--------------
-		 *            w4--------------
-		 *               w5-------------
-		 * for v, the longest out-edge is v->w5
-		 **/
-		L = asg_arc_len(av[nv-1]) + fuzz;
-
-
-		for (i = 0; i < nv; ++i) {
-			//w is an out-node of v
-			uint32_t w = av[i].v;
-			
-			uint32_t j, nw = asg_arc_n(g, w);
-			asg_arc_t *aw = asg_arc_a(g, w);
-			///if w has already been reduced
-			if (mark[av[i].v] != 1) continue;
-
-			for (j = 0; j < nw && asg_arc_len(aw[j]) + asg_arc_len(av[i]) <= L; ++j)
-				if (mark[aw[j].v]) mark[aw[j].v] = 2;
-		}
-		#if 0
-		for (i = 0; i < nv; ++i) {
-			uint32_t w = av[i].v;
-			uint32_t j, nw = asg_arc_n(g, w);
-			asg_arc_t *aw = asg_arc_a(g, w);
-			for (j = 0; j < nw && (j == 0 || asg_arc_len(aw[j]) < fuzz); ++j)
-				if (mark[aw[j].v]) mark[aw[j].v] = 2;
-		}
-		#endif
-        //remove edges
-		for (i = 0; i < nv; ++i) {
-			if (mark[av[i].v] == 2) av[i].del = 1, ++n_reduced;
-			mark[av[i].v] = 0;
-		}
-	}
-	free(mark);
-
-    if(VERBOSE >= 1)
-    {
-	    fprintf(stderr, "[M::%s] transitively reduced %d arcs\n", __func__, n_reduced);
-    }
-
-	if (n_reduced) {
-		asg_cleanup(g);
-		asg_symm(g);
-	}
-
-    if(VERBOSE >= 1)
-    {
-        fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime);
-    }
-	return n_reduced;
-}
-
-
-
-
-
 ///max_ext is 4
 int asg_cut_tip(asg_t *g, int max_ext)
 {
@@ -10435,7 +10076,7 @@ long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, uint32_t negativ
     ///the reason is that each read has two direction (query->target, target->query)
     uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, flag, is_hap, k, to_replace;
     long long ll, base_maxLen, base_maxPositive, base_minNegative, base_minNonPositive, base_best_i, all_covex, curPositive, curNonPositive, curNegative;
-    long long tmp, max_stop_nodeLen, max_stop_baseLen;
+    long long tmp, max_stop_nodeLen, max_stop_baseLen, cur_weight = 0, max_weight = 0;
     uint32_t non_positive_flag = (uint32_t)-1;
     if(positive_flag == FATHER) non_positive_flag = MOTHER;
     if(positive_flag == MOTHER) non_positive_flag = FATHER;
@@ -10501,37 +10142,22 @@ long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, uint32_t negativ
                 }
 
                 to_replace = 0;
-                if(curNegative < base_minNegative)
+                cur_weight = (long long)curPositive - ((long long)curNegative + (long long)curNonPositive);
+                max_weight = (long long)base_maxPositive - ((long long)base_minNegative + (long long)base_minNonPositive);
+                if(cur_weight > max_weight)
                 {
                     to_replace = 1;
                 }
-                else if(curNegative == base_minNegative)
+                else if (cur_weight == max_weight)
                 {
-                    if(curNonPositive < base_minNonPositive)
+                    if(ll > base_maxLen)
                     {
                         to_replace = 1;
                     }
-                    else if(curNonPositive == base_minNonPositive)
-                    {
-                        if(curPositive > base_maxPositive)
-                        {
-                            to_replace = 1;
-                        }
-                        else if(curPositive == base_maxPositive)
-                        {
-                            if(ll > base_maxLen)
-                            {
-                                to_replace = 1;
-                            }
-                        }   
-                    }
                 }
 
 
-
+                if(base_best_i == -1 || to_replace == 1)
-
-
-                if(to_replace == 1)
                 {
                     base_minNegative = curNegative;
                     base_minNonPositive = curNonPositive;
@@ -10610,6 +10236,7 @@ long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, uint32_t negativ
 }
 
 
+
 int asg_arc_cut_long_equal_tips_assembly_complex(asg_t *g, ma_hit_t_alloc* reverse_sources, 
 long long miniedgeLen, uint32_t stops_threshold, R_to_U* ruIndex)
 {
@@ -11378,7 +11005,7 @@ void print_untig(ma_ug_t *g, uint32_t uId, const char* info, uint32_t is_print_r
     asg_arc_t *aw = asg_arc_a(nsg, uId);
     uint32_t nw = asg_arc_n(nsg, uId);
 
-    fprintf(stderr, "\n%s: c = %u\n", info, nsg->seq[uId>>1].c);
+    fprintf(stderr, "\n%s: c = %u, del: %u\n", info, nsg->seq[uId>>1].c, nsg->seq[uId>>1].del);
     fprintf(stderr, "%s(%u): direction = 0...\n", info, uId>>1);
     for (i = 0; i < nw; i++)
     {
@@ -11939,22 +11566,26 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_thre
 }
 
 void renew_longest_tip_by_drop(asg_t *g, ma_ug_t *ug, asg_arc_t *av, uint32_t nv,
-long long* base_maxLen, long long* base_maxLen_i, uint32_t stops_threshold, buf_t* b)
+long long* base_maxLen, long long* base_maxLen_i, uint32_t stops_threshold, buf_t* b, uint32_t trio_flag)
 {
+    if(trio_flag != FATHER && trio_flag != MOTHER) return;
     Trio_counter max, cur;
     memset(&max, 0, sizeof(Trio_counter));
     memset(&cur, 0, sizeof(Trio_counter));
-    long long ll, tmp, max_stop_nodeLen, max_stop_baseLen;
+    long long ll, tmp, max_stop_nodeLen, max_stop_baseLen, max_weight = 0, cur_weight = 0;
     uint32_t convex, i, return_flag, best_tip_i, best_tip_len;
     
     b->b.n = 0;
     return_flag = get_unitig(g, ug, av[(*base_maxLen_i)].v, &convex, &tmp, &ll, 
     &max_stop_nodeLen, &max_stop_baseLen, stops_threshold, b);
-
     if(return_flag!=END_TIPS) return;
+
+
     get_trio_labs(b, ug, &max);
     ///means this unitig might be at current haplotype
-    if(max.drop_occ<(max.total*TRIO_DROP_THRES)) return;
+    ///if(max.drop_occ<(max.total*TRIO_DROP_THRES)) return;
+
+
     best_tip_i = (*base_maxLen_i);
     best_tip_len = (*base_maxLen);
 
@@ -11970,16 +11601,33 @@ long long* base_maxLen, long long* base_maxLen_i, uint32_t stops_threshold, buf_
         if(return_flag!=END_TIPS) return;
         ///this tip should be long enough
         if(ll<(TRIO_DROP_LENGTH_THRES*(*base_maxLen))) continue;
+
         get_trio_labs(b, ug, &cur);
-        ///this unitig is very likly at another haplotype, ignore it
+
-        if(cur.drop_occ>=(cur.total*TRIO_DROP_THRES)) continue;
+        if(trio_flag == FATHER)
-        if(cur.drop_occ<max.drop_occ)
+        {
+            max_weight = (long long)max.father_occ - (long long)max.drop_occ - (long long)max.mother_occ;
+            cur_weight = (long long)cur.father_occ - (long long)cur.drop_occ - (long long)cur.mother_occ;
+            ///this unitig is very likly at another haplotype, ignore it
+            if((cur.drop_occ+cur.mother_occ)>=(cur.total*TRIO_DROP_THRES)) continue;
+        }
+
+        if(trio_flag == MOTHER)
+        {
+            max_weight = (long long)max.mother_occ - (long long)max.drop_occ - (long long)max.father_occ;
+            cur_weight = (long long)cur.mother_occ - (long long)cur.drop_occ - (long long)cur.father_occ;
+            ///this unitig is very likly at another haplotype, ignore it
+            if((cur.drop_occ+cur.father_occ)>=(cur.total*TRIO_DROP_THRES)) continue;
+        }
+
+
+        if(cur_weight > max_weight)
         {
             max = cur;
             best_tip_i = i;
             best_tip_len = ll;
         }
-        else if(cur.drop_occ==max.drop_occ && ll>best_tip_len)
+        else if(cur_weight == max_weight && ll>best_tip_len)
         {
             max = cur;
             best_tip_i = i;
@@ -11992,7 +11640,7 @@ long long* base_maxLen, long long* base_maxLen_i, uint32_t stops_threshold, buf_
 }
 
 int asg_arc_cut_trio_long_equal_tips_assembly(asg_t *g, ma_ug_t *ug, asg_t *read_sg, 
-ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex)
+ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_t trio_flag)
 {
     double startTime = Get_T();
     uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, flag, is_hap, n_tips, return_flag, k;
@@ -12046,7 +11694,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex)
         if(n_arc == n_tips)
         {
             renew_longest_tip_by_drop(g, ug, av, nv, &base_maxLen, 
-            &base_maxLen_i, 1, &b);
+            &base_maxLen_i, 1, &b, trio_flag);
         }
 
         for (i = 0; i < nv; i++)
@@ -12625,12 +12273,11 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate)
     asg_t *g = ug->g;
     uint32_t is_first = 1;
 
-    
     redo:
+    ///if(trio_flag == MOTHER) fprintf(stderr, "(0) c: %u, del: %u, n: %u\n", ug->g->seq[28141].c, ug->g->seq[28141].del, ug->u.a[28141].n);
     asg_pop_bubble_primary_trio(ug, bubble_dist, trio_flag, DROP);
     untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, trio_flag, DROP); 
     magic_trio_phasing(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate);
-    
     ///drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex);
     
     /**********debug**********/
@@ -12654,21 +12301,21 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate)
         ///need consider tangles
         ///asg_pop_bubble_primary(g, bubble_dist);
         asg_pop_bubble_primary_trio(ug, bubble_dist, trio_flag, DROP);
-
-
-        
-
         /**********debug**********/
         if(just_bubble_pop == 0)
         {
             ///need consider tangles
             asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 
             2, tip_drop_ratio);
-            asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex);
+
+            asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag);
+            
             asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex,
             2, tip_drop_ratio, stops_threshold);
+
             asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 
             2, ruIndex, stops_threshold);
+
             detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate,
             ruIndex);
             ///need consider tangles
@@ -12686,13 +12333,10 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate)
         2);
     }
 
-
-    
     resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, trio_flag, drop_ratio);
     drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex);
     all_to_all_deduplicate(ug, trio_flag, trio_drop_rate, reverse_sources, ruIndex);
 
-
     if(is_first)
     {
         is_first = 0;
@@ -12735,7 +12379,7 @@ float drop_ratio)
             ///need consider tangles
             asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 
             2, tip_drop_ratio);
-            asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex);
+            asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1);
             asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex,
             2, tip_drop_ratio, stops_threshold);
             asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 
@@ -12801,10 +12445,10 @@ void set_drop_trio_flag(ma_ug_t *ug)
 
 
 void update_unitig_graph(ma_ug_t* ug, asg_t* read_g, ma_hit_t_alloc* reverse_sources, 
-R_to_U* ruIndex, uint8_t flag, float drop_rate)
+R_to_U* ruIndex, uint8_t is_double_check, uint8_t flag, float drop_rate)
 {
     asg_t* nsg = ug->g;
-    uint32_t v, n_vtx = nsg->n_seq, k, rId, flag_occ, non_flag_occ, n_reduce = 1;
+    uint32_t v, n_vtx = nsg->n_seq, k, rId, flag_occ, non_flag_occ, hap_label_occ, n_reduce = 1;
     ma_utg_t *u;
 
     drop_semi_circle(ug, nsg, read_g, reverse_sources, ruIndex);
@@ -12821,16 +12465,20 @@ R_to_U* ruIndex, uint8_t flag, float drop_rate)
             if((get_real_length(nsg, v<<1, NULL)!=0) 
                 && (get_real_length(nsg, ((v<<1)^1), NULL)!=0)) continue;
 
-            flag_occ = non_flag_occ = 0;
+            flag_occ = non_flag_occ = hap_label_occ = 0;
             for (k = 0; k < u->n; k++)
             {
                 rId = u->a[k]>>33;
+                if(read_g->seq[rId].c == HAP_LABLE) hap_label_occ++;
                 if(R_INF.trio_flag[rId] == AMBIGU) continue;
                 if(R_INF.trio_flag[rId] == DROP) continue;
                 if(R_INF.trio_flag[rId] == flag) flag_occ++;
                 if(R_INF.trio_flag[rId] != flag) non_flag_occ++;
             }
 
+            if(hap_label_occ == u->n) continue;
+            if(is_double_check && non_flag_occ < u->n*DOUBLE_CHECK_THRES) continue;
+
             if(non_flag_occ > ((non_flag_occ+flag_occ)*drop_rate))
             {
                 if(u->m != 0)
@@ -13168,7 +12816,40 @@ void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* re
 
 }
 
+void update_hap_label(ma_ug_t *ug, asg_t* read_g)
+{
+    uint32_t v, n_vtx, k;
+    uint64_t rId;
 
+    if(ug == NULL)
+    {
+        n_vtx = read_g->n_seq;
+        for (v = 0; v < n_vtx; ++v) 
+        {
+            if(read_g->seq[v].del) continue;
+            if(read_g->seq[v].c != HAP_LABLE) continue;
+            read_g->seq[v].c = PRIMARY_LABLE;
+        }
+        return;
+    }
+
+
+    asg_t* nsg = ug->g;
+    n_vtx = nsg->n_seq;
+    ma_utg_t* u = NULL;
+    for (v = 0; v < n_vtx; ++v) 
+    {
+        if(nsg->seq[v].del) continue;
+        if(nsg->seq[v].c != HAP_LABLE) continue;
+        u = &((ug)->u.a[v]);
+        if(u->m == 0) continue;
+        for (k = 0; k < u->n; k++)
+        {
+            rId = u->a[k]>>33;
+            read_g->seq[rId].c = HAP_LABLE;
+        }
+    }
+}
 
 
 void adjust_utg_by_trio(ma_ug_t **ug, asg_t* read_g, uint8_t flag, float drop_rate,
@@ -13180,14 +12861,15 @@ kvec_asg_arc_t_warp* new_rtg_edges)
     asg_t* nsg = (*ug)->g;
     uint32_t v, n_vtx = nsg->n_seq;
     
+    //if(flag == MOTHER) print_untig_by_read(*ug, "m54329U_190617_231905/65340614/ccs", (uint32_t)-1, sources, reverse_sources, "beg1");
+    //if(flag == MOTHER) print_untig_by_read(*ug, "m54329U_190827_173812/67108993/ccs", (uint32_t)-1, sources, reverse_sources, "beg2");
+    
     /**
     kvec_t_u32_warp new_rtg_nodes;
     kv_init(new_rtg_nodes.a);
     **/
-   
+    update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, 0, flag, drop_rate);
-    update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, flag, drop_rate);
     adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex);
-     
     nsg = (*ug)->g;
     n_vtx = nsg->n_seq;
     for (v = 0; v < n_vtx; ++v) 
@@ -13201,9 +12883,13 @@ kvec_asg_arc_t_warp* new_rtg_edges)
     tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, 
     chimeric_rate, 0, 0, drop_ratio, flag, drop_rate);
 
+    ///if(flag == MOTHER) fprintf(stderr, "(o.1) c: %u, del: %u, n: %u\n", (*ug)->g->seq[28141].c, (*ug)->g->seq[28141].del, (*ug)->u.a[28141].n);
+    
     delete_useless_nodes(ug);
 
-    update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, flag, drop_rate);
+    update_hap_label(*ug, read_g);
+
+    update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, 0, flag, drop_rate);
 
     renew_utg(ug, read_g, new_rtg_edges);
 
@@ -13217,7 +12903,9 @@ kvec_asg_arc_t_warp* new_rtg_edges)
 
     renew_utg(ug, read_g, new_rtg_edges);
 
-    update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, flag, drop_rate);
+    update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, 1, flag, drop_rate);
+
+    update_hap_label(NULL, read_g);
 
     renew_utg(ug, read_g, new_rtg_edges);
 
@@ -13543,6 +13231,7 @@ uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop)
 {   
 	uint32_t i, n_pending = 0, is_first = 1, cur_m, cur_c, cur_np, cur_nc, to_replace, n_tips, tip_end;
 	uint64_t n_pop = 0;
+    long long cur_weight = -1, max_weight = -1;
 	///if this node has been deleted
 	if (g->seq[v0>>1].del || g->seq[v0>>1].c == ALTER_LABLE) return 0; // already deleted
 	///if ((uint32_t)g->idx[v0] < 2) return 0; // no bubbles
@@ -13644,43 +13333,61 @@ uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop)
                     }
                     cur_nc = utg->u.a[(w>>1)].n;
                 }
-                ///select the way with less negative_flag, less non_positive_flag, more positive_flag, more distance
+                ///BUG: select the path with less negative_flag, less non_positive_flag, more positive_flag, more distance
+                ///FIXED: select the path with less (negative_flag+non_positive_flag), more positive_flag, more distance
                 to_replace = 0;
-                if(m + cur_m < t->m)
+
+                /****************************may have bugs********************************/
+                cur_weight = (long long)(c + cur_c) - ((long long)(m + cur_m) + (long long)(np + cur_np));
+                max_weight = (long long)t->c - ((long long)t->m + (long long)t->np);
+                if(cur_weight > max_weight)
                 {
                     to_replace = 1;
                 }
-                else if(m + cur_m == t->m)
+                else if(cur_weight == max_weight)
                 {
-                    if(np + cur_np < t->np)
+                    if(nc + cur_nc > t->nc)
                     {
                         to_replace = 1;
                     }
-                    else if(np + cur_np == t->np)
+                    else if(nc + cur_nc == t->nc)
                     {
-                        if(c + cur_c > t->c)
+                        if(d + l > t->d)
                         {
                             to_replace = 1;
                         }
-                        else if(c + cur_c == t->c)
+                    }
+                }
+                /****************************may have bugs********************************/
+
+                /**
+                if(((m + cur_m) + (np + cur_np)) < (t->m + t->np))
+                {
+                    to_replace = 1;
+                }
+                else if(((m + cur_m) + (np + cur_np)) == (t->m + t->np))
+                {
+                    if(c + cur_c > t->c)
+                    {
+                        to_replace = 1;
+                    }
+                    else if(c + cur_c == t->c)
+                    {
+                        if(nc + cur_nc > t->nc)
                         {
-                            if(nc + cur_nc > t->nc)
+                            to_replace = 1;
+                        }
+                        else if(nc + cur_nc == t->nc)
+                        {
+                            if(d + l > t->d)
                             {
                                 to_replace = 1;
                             }
-                            else if(nc + cur_nc == t->nc)
-                            {
-                                if(d + l > t->d)
-                                {
-                                    to_replace = 1;
-                                }
-                            }
-                            
                         }
                         
                     }
-                    
                 }
+                **/
                 
 
                 if(to_replace)
@@ -26021,13 +25728,12 @@ int max_hang, int min_ovlp)
 
 
 
-void clean_graph(
+void clean_graph(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
-int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
+		long long n_read, uint64_t* readLen, long long mini_overlap_length,
-long long n_read, uint64_t* readLen, long long mini_overlap_length, 
+		long long max_hang_length, long long clean_round, long long gap_fuzz,
-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,
-float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, 
+		long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr,
-long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr, 
+		ma_sub_t **coverage_cut_ptr, int debug_g)
-ma_sub_t **coverage_cut_ptr, int debug_g)
 {
 	ma_sub_t *coverage_cut = *coverage_cut_ptr;
 	asg_t *sg = *sg_ptr;
@@ -26062,7 +25768,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     ///ma_hit_sub is just use to init coverage_cut,
     ///it seems we do not need ma_hit_cut & ma_hit_flt
     ma_hit_sub(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut);
-    detect_chimeric_reads(sources, n_read, readLen, coverage_cut, FINAL_OVERLAP_ERROR_RATE*2);
+    detect_chimeric_reads(sources, n_read, readLen, coverage_cut, asm_opt.max_ov_diff_final * 2.0);
     ma_hit_cut(sources, n_read, readLen, mini_overlap_length, &coverage_cut);
     ///print_binned_reads(sources, n_read, coverage_cut);
     ma_hit_flt(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
@@ -26265,52 +25971,50 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
 	*sg_ptr = sg;
 }
 
-void build_string_graph_without_clean(
+void build_string_graph_without_clean(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
-int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
+		long long n_read, uint64_t* readLen, long long mini_overlap_length,
-long long n_read, uint64_t* readLen, long long mini_overlap_length, 
+		long long max_hang_length, long long clean_round, long long gap_fuzz,
-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,
-float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, 
+		long long bubble_dist, int read_graph, int write)
-long long bubble_dist, int read_graph, int write)
 {
-    R_to_U ruIndex;
+	R_to_U ruIndex;
-    init_R_to_U(&ruIndex, n_read);
+	init_R_to_U(&ruIndex, n_read);
-    asg_t *sg = NULL;
+	asg_t *sg = NULL;
-    ma_sub_t* coverage_cut = NULL;
+	ma_sub_t* coverage_cut = NULL;
 
-    // debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs", sources, reverse_sources, -1, "clean");
+	// debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs", sources, reverse_sources, -1, "clean");
 
-    ///actually min_thres = asm_opt.max_short_tip + 1 there are asm_opt.max_short_tip reads
+	///actually min_thres = asm_opt.max_short_tip + 1 there are asm_opt.max_short_tip reads
-    min_thres = asm_opt.max_short_tip + 1;
+	min_thres = asm_opt.max_short_tip + 1;
 
-    if (asm_opt.flag & HA_F_VERBOSE_GFA)
+	if (asm_opt.flag & HA_F_VERBOSE_GFA)
-    {
+	{
-        if(load_debug_graph(&sg, &sources, &coverage_cut, output_file_name, &reverse_sources, &ruIndex))
+		if (load_debug_graph(&sg, &sources, &coverage_cut, output_file_name, &reverse_sources, &ruIndex))
-        {
+		{
-            fprintf(stderr, "debug gfa has been loaded\n");
+			fprintf(stderr, "debug gfa has been loaded\n");
 
-            clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length, 
+			clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length,
-            max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 
+					max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
-            output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 1);
+					output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 1);
-            asg_destroy(sg);
+			asg_destroy(sg);
-            free(coverage_cut);
+			free(coverage_cut);
-            destory_R_to_U(&ruIndex);
+			destory_R_to_U(&ruIndex);
-            return;
+			return;
-        }
+		}
-    }
+	}
 
-    if (asm_opt.write_index_to_disk && write)
+	if (asm_opt.write_index_to_disk && write)
-    {
+	{
-        write_all_data_to_disk(sources, reverse_sources, 
+		write_all_data_to_disk(sources, reverse_sources, &R_INF, output_file_name);
-        &R_INF, output_file_name);
+	}
-    }
 
-    try_rescue_overlaps(sources, reverse_sources, n_read, 4); 
+	try_rescue_overlaps(sources, reverse_sources, n_read, 4);
 
-    clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length, 
+	clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length,
-    max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 
+			max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
-    output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 0);
+			output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 0);
 
-    asg_destroy(sg);
+	asg_destroy(sg);
-    free(coverage_cut);
+	free(coverage_cut);
-    destory_R_to_U(&ruIndex);
+	destory_R_to_U(&ruIndex);
 }

Overlaps.h

@@ -22,6 +22,7 @@
 #define ALTER_LABLE 1
 #define HAP_LABLE 2
 #define TRIO_THRES 0.9
+#define DOUBLE_CHECK_THRES 0.2
 // #define PRIMARY_LABLE 1
 // #define ALTER_LABLE 2
 // #define HAP_LABLE 4
@@ -102,7 +103,6 @@ typedef struct {
 	uint8_t no_l_indel;
 } asg_arc_t;
 
-
 typedef struct {
 	uint32_t len:31, circ:1; // len: length of the unitig; circ: circular if non-zero
 	uint32_t start, end; // start: starting vertex in the string graph; end: ending vertex
@@ -111,11 +111,10 @@ typedef struct {
 	char *s; // unitig sequence is not null
 } ma_utg_t;
 
-
-
 typedef struct {
 	uint32_t len:31, del:1;
 	uint8_t c;
+	uint8_t weak_contain;
 } asg_seq_t;
 
 typedef struct {
@@ -127,6 +126,10 @@ typedef struct {
 	asg_seq_t *seq;
 	uint64_t *idx;
 
+	uint32_t m_con, n_con;
+	uint64_t *contain;
+	uint64_t *con_idx;
+
 	uint8_t* seq_vis;
 
 	uint32_t n_F_seq;
@@ -504,6 +507,25 @@ void init_Edge_iter(asg_t* g, uint32_t v, asg_arc_t* new_edges, uint32_t new_edg
 int get_arc_t(Edge_iter* x, asg_arc_t* get);
 int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_flag, uint32_t negative_flag);
 
+/*************************************
+ * Routines modified for containment *
+ *************************************/
+
+void delete_single_edge(ma_hit_t_alloc *sources, const ma_sub_t *coverage_cut, uint32_t qn, uint32_t tn);
+void delete_all_edges(ma_hit_t_alloc *sources, ma_sub_t *coverage_cut, uint32_t qn);
+void ma_hit_contained_advance(ma_hit_t_alloc *sources, long long n_read, ma_sub_t *coverage_cut, R_to_U *ruIndex, int max_hang, int min_ovlp);
+void asg_con_sort(asg_t *g);
+void asg_con_index(asg_t *g);
+asg_t *ma_sg_gen(const ma_hit_t_alloc* sources, long long n_read, const ma_sub_t *coverage_cut, int max_hang, int min_ovlp);
+int asg_arc_del_trans(asg_t *g, int fuzz);
+int asg_drop_contained_utg(asg_t *g);
+
+#define asg_con_n(g, v) ((g)->contain? (uint32_t)(g)->con_idx[(v)] : 0)
+#define asg_con_a(g, v) ((g)->contain? &(g)->contain[(g)->con_idx[(v)]>>32] : 0)
+
+/*******************
+ * Other rountines *
+ *******************/
 
 inline int get_real_length(asg_t *g, uint32_t v, uint32_t* v_s)
 {

Trio.cpp

@@ -9,6 +9,7 @@
 #include "kseq.h"
 #include "Process_Read.h"
 #include "htab.h"
+#include "utils.h"
 #include "CommandLines.h"
 
 #define YAK_MAX_KMER     31

anchor.cpp

@@ -1,6 +1,7 @@
 #include <stdio.h>
 #include "htab.h"
 #include "ksort.h"
+#include "utils.h"
 #include "Hash_Table.h"
 
 #define HA_KMER_GOOD_RATIO 0.333
@@ -50,11 +51,18 @@ 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)
+{
+	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)
 {
 	extern void *ha_flt_tab;
 	extern ha_pt_t *ha_idx;
-	uint32_t i;
+	uint32_t i, rlen;
 	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);
@@ -64,6 +72,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
     clear_overlap_region_alloc(overlap_list);
 	recover_UC_Read(ucr, &R_INF, rid);
 	ab->mz.n = 0, ab->n_a = 0;
+	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);
@@ -126,30 +135,31 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	calculate_overlap_region_by_chaining(cl, overlap_list, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain);
 
 	#if 0
-	if (overlap_list->length > 2000) {
+	if (overlap_list->length > 0) {
-		fprintf(stderr, "B\t%ld\t%ld\t%ld\n", (long)rid, (long)overlap_list->length, (long)Get_READ_LENGTH(R_INF, rid));
+		fprintf(stderr, "B\t%ld\t%ld\t%d\n", (long)rid, (long)overlap_list->length, rlen);
 		for (int i = 0; i < (int)overlap_list->length; ++i) {
 			overlap_region *r = &overlap_list->list[i];
-			fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand],
+			fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand],
-					(int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed);
+					(int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed, ha_ov_type(r, rlen));
 		}
 	}
 	#endif
 
 	if ((int)overlap_list->length > max_n_chain) {
-		int32_t n[2], s[2];
+		int32_t w, n[4], s[4];
-		n[0] = n[1] = 0, s[0] = s[1] = 0;
+		n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0;
+		ks_introsort_or_ss(overlap_list->length, overlap_list->list);
 		for (i = 0; i < (uint32_t)overlap_list->length; ++i) {
 			const overlap_region *r = &overlap_list->list[i];
-			int dir = r->x_pos_s == 0? 0 : 1;
+			w = ha_ov_type(r, rlen);
-			++n[dir];
+			++n[w];
-			if ((int)n[dir] == max_n_chain) s[dir] = r->shared_seed;
+			if ((int)n[w] == max_n_chain) s[w] = r->shared_seed;
 		}
-		if (s[0] > 0 || s[1] > 0) {
+		if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) {
 			for (i = 0, k = 0; i < (uint32_t)overlap_list->length; ++i) {
 				overlap_region *r = &overlap_list->list[i];
-				int dir = r->x_pos_s == 0? 0 : 1;
+				w = ha_ov_type(r, rlen);
-				if (r->shared_seed > s[dir]) {
+				if (r->shared_seed >= s[w]) {
 					if ((uint32_t)k != i) {
 						overlap_region t;
 						t = overlap_list->list[k];

extract.cpp

@@ -0,0 +1,171 @@
+#include <zlib.h>
+#include <string.h>
+#include "Process_Read.h"
+#include "khashl.h"
+#include "kseq.h"
+#include "utils.h"
+
+typedef const char *cstr_t;
+KHASHL_CSET_INIT(KH_LOCAL, strset_t, ss, cstr_t, kh_hash_str, kh_eq_str)
+KHASHL_MAP_INIT(KH_LOCAL, hm64_t, h64, uint64_t, int, kh_hash_uint64, kh_eq_generic)
+KSTREAM_INIT(gzFile, gzread, 65536)
+
+char *gfa_strdup(const char *src)
+{
+	int32_t len;
+	char *dst;
+	len = strlen(src);
+	MALLOC(dst, len + 1);
+	memcpy(dst, src, len + 1);
+	return dst;
+}
+
+char *gfa_strndup(const char *src, size_t n)
+{
+	char *dst;
+	MALLOC(dst, n + 1);
+	strncpy(dst, src, n);
+	dst[n] = 0;
+	return dst;
+}
+
+char **gv_read_list(const char *o, int *n_)
+{
+	int n = 0, m = 0;
+	char **s = 0;
+	*n_ = 0;
+	if (*o != '@') {
+		const char *q = o, *p;
+		for (p = q;; ++p) {
+			if (*p == ',' || *p == 0) {
+				if (n == m) {
+					m = m? m<<1 : 16;
+					REALLOC(s, m);
+				}
+				s[n++] = gfa_strndup(q, p - q);
+				if (*p == 0) break;
+				q = p + 1;
+			}
+		}
+	} else {
+		gzFile fp;
+		kstream_t *ks;
+		kstring_t str = {0,0,0};
+		int dret;
+
+		fp = gzopen(o + 1, "r");
+		if (fp == 0) return 0;
+		ks = ks_init(fp);
+		while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) {
+			char *p;
+			for (p = str.s; *p && !isspace(*p); ++p);
+			if (n == m) {
+				m = m? m<<1 : 16;
+				REALLOC(s, m);
+			}
+			s[n++] = gfa_strndup(str.s, p - str.s);
+		}
+		ks_destroy(ks);
+		gzclose(fp);
+	}
+	if (s) s = (char**)realloc(s, n * sizeof(char*));
+	*n_ = n;
+	return s;
+}
+
+void ha_extract_print(const All_reads *rs, int n_rounds, int n, char **list)
+{
+	hm64_t *h;
+	khint_t k;
+	int i, absent, m, l;
+	uint64_t j;
+	const ma_hit_t_alloc *ov[2] = { rs->paf, rs->reverse_paf };
+	FILE *fp = stdout;
+
+	if (n > 0) {
+		int max_len = 0;
+		char *s = 0;
+		strset_t *ss;
+		ss = ss_init();
+		for (i = 0; i < n; ++i)
+			ss_put(ss, list[i], &absent);
+		for (j = 0; j < rs->total_reads; ++j)
+			if (max_len < (int)Get_NAME_LENGTH(*rs, j))
+				max_len = Get_NAME_LENGTH(*rs, j);
+		MALLOC(s, max_len + 1);
+		h = h64_init();
+		for (j = 0; j < rs->total_reads; ++j) {
+			strncpy(s, Get_NAME(*rs, j), Get_NAME_LENGTH(*rs, j));
+			s[Get_NAME_LENGTH(*rs, j)] = 0;
+			if (ss_get(ss, s) != kh_end(ss)) {
+				k = h64_put(h, j, &absent);
+				kh_val(h, k) = -1;
+			}
+		}
+		free(s);
+		ss_destroy(ss);
+	} else return;
+
+	for (m = 0; m < n_rounds; ++m) {
+		for (j = 0; j < rs->total_reads; ++j) {
+			for (l = 0; l < 2; ++l) {
+				const ma_hit_t_alloc *o = &ov[l][j];
+				for (i = 0; i < (int)o->length; ++i) {
+					uint64_t q = Get_qn(o->buffer[i]);
+					uint64_t t = Get_tn(o->buffer[i]);
+					int q_hit = 0, t_hit = 0;
+					k = h64_get(h, q);
+					q_hit = (k < kh_end(h) && kh_val(h, k) < m);
+					k = h64_get(h, t);
+					t_hit = (k < kh_end(h) && kh_val(h, k) < m);
+					if ((!q_hit && !t_hit) || (q_hit && t_hit)) continue;
+					if (!q_hit) {
+						k = h64_put(h, q, &absent);
+						if (absent) kh_val(h, k) = m;
+					}
+					if (!t_hit) {
+						k = h64_put(h, t, &absent);
+						if (absent) kh_val(h, k) = m;
+					}
+				}
+			}
+		}
+	}
+
+	for (j = 0; j < rs->total_reads; ++j) {
+		for (l = 0; l < 2; ++l) {
+			const ma_hit_t_alloc *o = &ov[l][j];
+			for (i = 0; i < (int)o->length; ++i) {
+				uint64_t q = Get_qn(o->buffer[i]);
+				uint64_t t = Get_tn(o->buffer[i]);
+				int q_hit = 0, t_hit = 0;
+				q_hit = (h64_get(h, q) < kh_end(h));
+				t_hit = (h64_get(h, t) < kh_end(h));
+				if (!q_hit && !t_hit) continue;
+				fwrite(Get_NAME(*rs, q), 1, Get_NAME_LENGTH(*rs, q), fp);
+				fwrite("\t", 1, 1, fp);
+				fprintf(fp, "%lu\t", (unsigned long)Get_READ_LENGTH(*rs, q));
+				fprintf(fp, "%d\t", Get_qs(o->buffer[i]));
+				fprintf(fp, "%d\t", Get_qe(o->buffer[i]));
+				fputs(o->buffer[i].rev? "-\t" : "+\t", fp);
+				fwrite(Get_NAME(*rs, t), 1, Get_NAME_LENGTH(*rs, t), fp);
+				fwrite("\t", 1, 1, fp);
+				fprintf(fp, "%lu\t", (unsigned long)Get_READ_LENGTH(*rs, t));
+				fprintf(fp, "%d\t", Get_ts(o->buffer[i]));
+				fprintf(fp, "%d\t%d\t%d\t%d\n", Get_te(o->buffer[i]), o->buffer[i].ml, o->buffer[i].bl, !l);
+			}
+		}
+	}
+
+	h64_destroy(h);
+}
+
+void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o)
+{
+	int i, n;
+	char **list;
+	list = gv_read_list(o, &n);
+	ha_extract_print(rs, n_rounds, n, list);
+	for (i = 0; i < n; ++i) free(list[i]);
+	free(list);
+}

htab.cpp

@@ -9,6 +9,7 @@
 #include "kseq.h"
 #include "ksort.h"
 #include "htab.h"
+#include "utils.h"
 
 #define YAK_COUNTER_BITS 12
 #define YAK_N_COUNTS     (1<<YAK_COUNTER_BITS)

htab.h

@@ -38,13 +38,6 @@ 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);
 
-double yak_cputime(void);
-void yak_reset_realtime(void);
-double yak_realtime(void);
-long yak_peakrss(void);
-double yak_peakrss_in_gb(void);
-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);
@@ -80,18 +73,6 @@ static inline uint64_t yak_hash_long(uint64_t x[4])
 	return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]);
 }
 
-#define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr))))
-#define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr))))
-#define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr))))
-
-#ifndef kroundup32
-#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
-#endif
-
-#ifndef kroundup64
-#define kroundup64(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, x|=(x)>>32, ++(x))
-#endif
-
 #ifndef klib_unused
 #if (defined __clang__ && __clang_major__ >= 3) || (defined __GNUC__ && __GNUC__ >= 3)
 #define klib_unused __attribute__ ((__unused__))

main.cpp

@@ -4,7 +4,7 @@
 #include "Process_Read.h"
 #include "Assembly.h"
 #include "Levenshtein_distance.h"
-#include "htab.h"
+#include "utils.h"
 
 int main(int argc, char *argv[])
 {

overlap2.cpp

@@ -0,0 +1,407 @@
+#include <assert.h>
+#include "utils.h"
+#include "CommandLines.h"
+#include "Overlaps.h"
+#include "Process_Read.h"
+
+/*******************************
+ * Dropping strong containment *
+ *******************************/
+
+static ma_hit_t *get_specific_overlap_with_del(ma_hit_t_alloc *sources, const ma_sub_t *coverage_cut, uint32_t qn, uint32_t tn)
+{
+	if (coverage_cut[qn].del || coverage_cut[tn].del) return NULL;
+	ma_hit_t_alloc *x = &sources[qn];
+	uint32_t i;
+	for (i = 0; i < x->length; i++) {
+		if (x->buffer[i].del) continue;
+		if (coverage_cut[Get_qn(x->buffer[i])].del) continue;
+		if (coverage_cut[Get_tn(x->buffer[i])].del) continue;
+		if (Get_tn(x->buffer[i]) == tn && Get_qn(x->buffer[i]) == qn)
+			return &x->buffer[i];
+	}
+	return NULL;
+}
+
+void delete_single_edge(ma_hit_t_alloc *sources, const ma_sub_t *coverage_cut, uint32_t qn, uint32_t tn)
+{
+	ma_hit_t *tmp = get_specific_overlap_with_del(sources, coverage_cut, qn, tn);
+	if (tmp != NULL) tmp->del = 1;
+}
+
+void delete_all_edges(ma_hit_t_alloc *sources, ma_sub_t *coverage_cut, uint32_t qn)
+{
+	ma_hit_t_alloc* x = &sources[qn];
+	uint32_t i;
+	for (i = 0; i < x->length; i++) {
+		x->buffer[i].del = 1;
+		delete_single_edge(sources, coverage_cut, Get_tn(x->buffer[i]), Get_qn(x->buffer[i]));
+	}
+	coverage_cut[qn].del = 1;
+}
+
+void ma_hit_contained_advance(ma_hit_t_alloc *sources, long long n_read, ma_sub_t *coverage_cut, R_to_U *ruIndex, int max_hang, int min_ovlp)
+{
+	int32_t r;
+	long long i, j, n_strong_contain = 0, n_weak_contain = 0;
+	asg_arc_t t;
+	ma_hit_t *h = NULL;
+	ma_sub_t *sq = NULL;
+	ma_sub_t *st = NULL;
+
+	for (i = 0; i < n_read; ++i) {
+		if (coverage_cut[i].del) continue;
+		for (j = 0; j < (long long)sources[i].length; j++) {
+			h = &sources[i].buffer[j];
+			//check the corresponding two reads 
+			sq = &coverage_cut[Get_qn(*h)];
+			st = &coverage_cut[Get_tn(*h)];
+			/****************************may have trio bugs********************************/
+			if (sq->del || st->del) continue;
+			if (h->del) continue;
+			/****************************may have trio bugs********************************/
+			r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
+			//assert(r != MA_HT_INT && r != MA_HT_SHORT_OVLP);
+			if (r == MA_HT_QCONT) {
+				if (h->ml || (asm_opt.flag & HA_F_KEEP_CONTAINED) == 0) {
+					h->del = 1;
+					delete_single_edge(sources, coverage_cut, Get_tn(*h), Get_qn(*h));
+					delete_all_edges(sources, coverage_cut, Get_qn(*h));
+					set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0);
+
+					// if (delete_all_edges_carefully(sources, coverage_cut, max_hang, min_ovlp, Get_qn(*h)) == 0)
+					//     set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0);
+					// sq->del = 1;
+					// set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0);
+				}
+				if (h->ml) ++n_strong_contain;
+				else ++n_weak_contain;
+			} else if (r == MA_HT_TCONT) {
+				if (h->ml || (asm_opt.flag & HA_F_KEEP_CONTAINED) == 0) {
+					h->del = 1;
+					delete_single_edge(sources, coverage_cut, Get_tn(*h), Get_qn(*h));
+					delete_all_edges(sources, coverage_cut, Get_tn(*h));
+					set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0);
+
+					// if (delete_all_edges_carefully(sources, coverage_cut, max_hang, min_ovlp, Get_tn(*h)) == 0)
+					//     set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0);
+					// st->del = 1;
+					// set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0);
+				}
+				if (h->ml) ++n_strong_contain;
+				else ++n_weak_contain;
+			}
+		}
+	}
+
+	transfor_R_to_U(ruIndex);
+
+	for (i = 0; i < n_read; ++i) {
+		int m = 0;
+		for (j = 0; j < (long long)sources[i].length; j++) {
+			ma_hit_t *h = &(sources[i].buffer[j]);
+			if (h->del) continue;
+			/// both the qn and tn have not been deleted
+			if (coverage_cut[Get_qn(*h)].del != 1 && coverage_cut[Get_tn(*h)].del != 1)
+				h->del = 0, ++m;
+			else h->del = 1;
+		}
+		/// sources[i].length == 0 means all overlapped reads with read i are the contained reads
+		if (m == 0) coverage_cut[i].del = 1;
+	}
+
+	fprintf(stderr, "[M::%s] %lld strong containments; %lld weak containments\n", __func__,
+			n_strong_contain, n_weak_contain);
+}
+
+/************************************
+ * Graph construction and reduction *
+ ************************************/
+
+static inline void asg_con_push(asg_t *g, uint32_t lower, uint32_t upper, int rev)
+{
+	if (g->n_con == g->m_con) {
+		g->m_con = g->m_con? g->m_con<<1 : 16;
+		REALLOC(g->contain, g->m_con);
+	}
+	g->contain[g->n_con++] = (uint64_t)lower << 32 | upper << 1 | (!!rev);
+}
+
+void asg_con_sort(asg_t *g)
+{
+	if (g->n_con > 1) radix_sort_ha64(g->contain, g->contain + g->n_con);
+}
+
+void asg_con_index(asg_t *g)
+{
+	uint32_t i, k;
+	if (g->n_con == 0 || g->contain == 0) return;
+	if (g->con_idx) free(g->con_idx);
+	CALLOC(g->con_idx, g->n_seq);
+	for (k = 0, i = 1; i < g->n_con; ++i)
+		if (g->contain[k] != g->contain[i])
+			g->contain[k++] = g->contain[i];
+	g->n_con = k;
+	for (i = 1, k = 0; i <= g->n_con; ++i)
+		if (i == g->n_con || g->contain[i-1]>>32 != g->contain[i]>>32)
+			g->con_idx[g->contain[i-1]>>32] = (uint64_t)k << 32 | (i - k), k = i;
+}
+
+asg_t *ma_sg_gen(const ma_hit_t_alloc* sources, long long n_read, const ma_sub_t *coverage_cut, int max_hang, int min_ovlp)
+{
+	size_t i, j;
+	asg_t *g;
+	g = asg_init();
+
+	// add seq to graph, seq just save the length of each read
+	for (i = 0; i < (uint64_t)n_read; ++i) {
+		///if a read has been deleted, should we still add them?
+		asg_seq_set(g, i, coverage_cut[i].e - coverage_cut[i].s, coverage_cut[i].del);
+		g->seq[i].c = coverage_cut[i].c;
+	}
+
+	g->seq_vis = (uint8_t*)calloc(g->n_seq*2, sizeof(uint8_t));
+
+	for (i = 0; i < (uint64_t)n_read; ++i) {
+		for (j = 0; j < sources[i].length; ++j) {
+			int r, ql, tl;
+			asg_arc_t t, *p;
+			const ma_hit_t *h = &sources[i].buffer[j];
+			uint32_t qn, tn;
+			if (h->del) continue;
+			qn = Get_qn(*h);
+			tn = Get_tn(*h);
+			ql = coverage_cut[qn].e - coverage_cut[qn].s;
+			tl = coverage_cut[tn].e - coverage_cut[tn].s;
+			r = ma_hit2arc(h, ql, tl, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
+			if (r >= 0) {
+				p = asg_arc_pushp(g);
+				*p = t;
+			} else if (r == MA_HT_QCONT) {
+				assert(h->ml == 0);
+				asg_con_push(g, h->qns>>32, h->tn, h->rev);
+			} else if (r == MA_HT_TCONT) {
+				assert(h->ml == 0);
+				asg_con_push(g, h->tn, h->qns>>32, h->rev);
+			}
+		}
+	}
+	asg_cleanup(g);
+	g->r_seq = g->n_seq;
+	return g;
+}
+
+typedef struct {
+	uint32_t len;
+	uint8_t mark; // can only be 0, 1 or 2
+} trinfo_t;
+
+// transitive reduction; see Myers, 2005
+int asg_arc_del_trans(asg_t *g, int fuzz)
+{
+	trinfo_t *info;
+	///n_vtx = number of seq * 2; the reason is that each read has two direction (query->target, target->query)
+	uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0;
+	///at first, all nodes should be set to vacant
+	CALLOC(info, n_vtx);
+
+	/**v is the id+direction of a node, 
+	 * the high 31-bit is the id, 
+	 * and the lowest 1-bit is the direction
+	 * (0 means query-to-target, 1 means target-to-query)**/
+	for (v = 0; v < n_vtx; ++v) {
+		///nv is the number of overlaps with v(qn+direction)
+		uint32_t L, i, nv = asg_arc_n(g, v);
+		///av is the array of v
+		asg_arc_t *av = asg_arc_a(g, v);
+		///that means in this direction, read v is not overlapped with any other reads
+		if (nv == 0) continue; // no hits
+
+		// if the read itself has been removed
+		if (g->seq[v>>1].del) {
+			for (i = 0; i < nv; ++i) av[i].del = 1, ++n_reduced;
+			continue;
+		}
+
+    /**
+	********************************query-to-target overlap****************************
+	case 1: u = 0, rev = 0                           in the view of target: direction is 1 
+	query: CCCCCCCCTAATTAAAAT                        target: TAATTAAAATGGGGGG (use ex-target as query)
+	               ||||||||||         <--->                  ||||||||||
+	       target: TAATTAAAATGGGGGG           query: CCCCCCCCTAATTAAAAT (use ex-query as target)
+
+	case 2: u = 0, rev = 1                           in the view of target: direction is 0
+	query: CCCCCCCCTAATTAAAAT					     target: CCCCCCATTTTAATTA  (use ex-target as query)
+                   ||||||||||        <--->                         ||||||||||
+	       target: TAATTAAAATGGGGGG                         query: ATTTTAATTAGGGGGGGG  (use ex-query as target)
+	********************************query-to-target overlap****************************
+
+	********************************target-to-query overlap****************************
+	case 3: u = 1, rev = 0                           in the view of target: direction is 0
+			 query: AAATAATATCCCCCCGCG                target: GGGCCGGCAAATAATAT (use ex-target as query)
+					|||||||||          <--->                          |||||||||
+	target: GGGCCGGCAAATAATAT                                  query: AAATAATATCCCCCCGCG (use ex-query as target)
+
+	case 4: u = 1, rev = 1                          in the view of target: direction is 1
+	         query: AAATAATATCCCCCCGCG                        target: ATATTATTTGCCGGCCC (use ex-target as query)
+                    |||||||||             <--->                       |||||||||
+	target: GGGCCGGCAAATAATAT                          query: CGCGGGGGATATTATTT (use ex-query as target)
+	********************************target-to-query overlap****************************
+
+	p->ul: |____________31__________|__________1___________|______________32_____________|
+	                    qns            direction of overlap       length of this node (not overlap length)
+						                (in the view of query)
+	p->v : |___________31___________|__________1___________|
+				        tns             reverse direction of overlap
+						              (in the view of target)
+	p->ol: overlap length
+    **/
+
+		// all outnode of v should be set to "not reduce"
+		for (i = 0; i < nv; ++i) {
+			uint32_t w = av[i].v;
+			info[w].mark = g->seq[w>>1].del? 2 : 1;
+			//if (asg_con_n(g, w>>1) > 0) info[w].mark = 2;
+			info[w].len = asg_arc_len(av[i]);
+		}
+
+		// remove contained reads
+		for (i = 0; i < nv; ++i) {
+			uint32_t j, nw, w = av[i].v;
+			uint64_t *aw;
+			if (info[w].mark != 1) continue;
+			nw = asg_con_n(g, w>>1);
+			if (nw == 0) continue;
+			aw = asg_con_a(g, w>>1);
+			for (j = 0; j < nw; ++j) {
+				uint32_t x = (uint32_t)aw[j];
+				if (w&1) x ^= 1;
+				if (info[x].mark == 1 && info[x].len <= info[w].len)
+					break;
+			}
+			if (j < nw) info[w].mark = 2;
+			//if (nw > 0) fprintf(stderr, "X\t%.*s\t%.*s\n", (int)Get_NAME_LENGTH(R_INF, w>>1), Get_NAME(R_INF, w>>1), (int)Get_NAME_LENGTH(R_INF, (uint32_t)aw[j]>>1), Get_NAME(R_INF, (uint32_t)aw[j]>>1));
+		}
+
+		// length of node (not overlap length)
+		// av[nv-1] is longest out-dege
+		/**
+		 * v---------------
+		 *   w1---------------
+		 *      w2--------------
+		 *         w3--------------
+		 *            w4--------------
+		 *               w5-------------
+		 * for v, the longest out-edge is v->w5
+		 **/
+		L = asg_arc_len(av[nv-1]) + fuzz;
+
+		for (i = 0; i < nv; ++i) {
+			uint32_t w = av[i].v;
+			uint32_t j, nw = asg_arc_n(g, w);
+			uint32_t is_con = (asg_con_n(g, w>>1) > 0);
+			asg_arc_t *aw = asg_arc_a(g, w);
+			if (info[w].mark != 1) continue;
+			for (j = 0; j < nw; ++j) {
+				uint32_t x, sum = asg_arc_len(aw[j]) + asg_arc_len(av[i]);
+				if (sum > L) break;
+				x = aw[j].v;
+				if (info[x].mark == 1 && sum < info[x].len + fuzz && sum + fuzz > info[x].len) {
+					if (!is_con || asg_con_n(g, x>>1) > 0)
+						info[x].mark = 2;
+				}
+			}
+		}
+		#if 0
+		for (i = 0; i < nv; ++i) {
+			uint32_t w = av[i].v;
+			uint32_t j, nw = asg_arc_n(g, w);
+			asg_arc_t *aw = asg_arc_a(g, w);
+			for (j = 0; j < nw && (j == 0 || asg_arc_len(aw[j]) < fuzz); ++j)
+				if (info[aw[j].v].mark) info[aw[j].v].mark = 2;
+		}
+		#endif
+
+		// remove edges
+		for (i = 0; i < nv; ++i) {
+			if (info[av[i].v].mark == 2) av[i].del = 1, ++n_reduced;
+			info[av[i].v].mark = 0;
+		}
+	}
+	free(info);
+
+	if (n_reduced) {
+		asg_cleanup(g);
+		asg_symm(g);
+		asg_drop_contained_utg(g);
+	}
+	fprintf(stderr, "[M::%s] transitively reduced %d arcs\n", __func__, n_reduced);
+	return n_reduced;
+}
+
+#define GFA_VT_MERGEABLE 0
+#define GFA_VT_TIP       1
+#define GFA_VT_MULTI_OUT 2
+#define GFA_VT_MULTI_IN  3
+
+static inline int32_t gfa_deg(const asg_t *g, uint32_t v, uint32_t *w)
+{
+	uint32_t i, nv, nv0, k;
+	const asg_arc_t *av;
+	if (w) *w = (uint32_t)-1;
+	if (g->seq[v>>1].del) return 0;
+	nv0 = k = asg_arc_n(g, v);
+	av = asg_arc_a(g, v);
+	for (i = nv = 0; i < nv0; ++i)
+		if (!av[i].del)
+			++nv, k = i;
+	if (w) *w = nv == 1? av[k].v : (uint32_t)-1;
+	return nv;
+}
+
+static inline int32_t gfa_vtype(const asg_t *g, uint32_t v, uint32_t *w_)
+{
+	int32_t nv, nw;
+	uint32_t w;
+	nv = gfa_deg(g, v, &w);
+	if (w_) *w_ = w;
+	if (nv == 0) return GFA_VT_TIP;
+	if (nv > 1) return GFA_VT_MULTI_OUT;
+	nw = gfa_deg(g, w^1, 0);
+	return nw == 1? GFA_VT_MERGEABLE : GFA_VT_MULTI_IN;
+}
+
+int asg_drop_contained_utg(asg_t *g)
+{
+	uint32_t n_vtx = g->n_seq * 2, v, cnt = 0;
+	if (g->contain == 0) return 0;
+	for (v = 0; v < n_vtx; ++v) {
+		int32_t vt, is_contained;
+		uint32_t w;
+		if (g->seq[v>>1].del) continue;
+		if (asg_con_n(g, v>>1) == 0) continue;
+		vt = gfa_vtype(g, v^1, &w);
+		if (vt == GFA_VT_MERGEABLE) continue;
+		w = v, is_contained = 1;
+		while (1) {
+			if (asg_con_n(g, w>>1) == 0) {
+				is_contained = 0;
+				break;
+			}
+			vt = gfa_vtype(g, w, &w);
+			if (vt != GFA_VT_MERGEABLE) break;
+		}
+		if (is_contained) {
+			w = v;
+			while (1) {
+				++cnt;
+				asg_seq_del(g, w>>1);
+				vt = gfa_vtype(g, w, &w);
+				if (vt != GFA_VT_MERGEABLE) break;
+			}
+		}
+	}
+	if (cnt > 0) asg_cleanup(g);
+	fprintf(stderr, "[M::%s] drop %d reads in contained unitigs\n", __func__, cnt);
+	return cnt;
+}

sys.cpp

@@ -1,6 +1,6 @@
 #include <sys/resource.h>
 #include <sys/time.h>
-#include "htab.h"
+#include "utils.h"
 
 int yak_verbose = 3;
 

utils.h

@@ -0,0 +1,33 @@
+#ifndef HA_UTILS_H
+#define HA_UTILS_H
+
+#include <stdint.h>
+
+#ifndef MALLOC
+#define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr))))
+#endif
+#ifndef CALLOC
+#define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr))))
+#endif
+#ifndef REALLOC
+#define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr))))
+#endif
+
+#ifndef kroundup32
+#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
+#endif
+
+#ifndef kroundup64
+#define kroundup64(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, x|=(x)>>32, ++(x))
+#endif
+
+void radix_sort_ha64(uint64_t *st, uint64_t *en);
+
+double yak_cputime(void);
+void yak_reset_realtime(void);
+double yak_realtime(void);
+long yak_peakrss(void);
+double yak_peakrss_in_gb(void);
+double yak_cpu_usage(void);
+
+#endif