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);
 

CommandLines.cpp

@@ -19,6 +19,7 @@ static ko_longopt_t long_options[] = {
 	{ "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 }
 };
 
@@ -404,6 +405,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 		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,7 +3,7 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.7-dirty-r255"
+#define HA_VERSION "0.7-dirty-r257"
 
 #define VERBOSE 0
 
@@ -15,6 +15,7 @@
 #define HA_F_SKIP_TRIOBIN    0x20
 #define HA_F_PURGE_CONTAIN   0x40
 #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
 

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)

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 extract.o sys.o
+			htab.o hist.o sketch.o anchor.o extract.o overlap2.o sys.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm
 
@@ -34,17 +34,17 @@ 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
@@ -52,16 +52,17 @@ 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 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
+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);
-     
-    if(g->n_F_seq > 0 && g->F_seq)
-    {
-        uint32_t i = 0;
-        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);
-        }
+	free(g->contain); free(g->con_idx);
 
-        free(g->F_seq);
-    }
-    
-
-    free(g);
+	if (g->n_F_seq > 0 && g->F_seq) {
+		uint32_t i = 0;
+		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);
+		}
+		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, 
-uint64_t* source_index, long long listLen)
+void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, 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();
@@ -2108,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;
-}
 
 
 
@@ -4484,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)
 {
@@ -26086,13 +25728,12 @@ int max_hang, int min_ovlp)
 
 
 
-void clean_graph(
-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 max_hang_length, long long clean_round, long long gap_fuzz,
-float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, 
-long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr, 
-ma_sub_t **coverage_cut_ptr, int debug_g)
+void clean_graph(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 max_hang_length, long long clean_round, long long gap_fuzz,
+		float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name,
+		long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr,
+		ma_sub_t **coverage_cut_ptr, int debug_g)
 {
 	ma_sub_t *coverage_cut = *coverage_cut_ptr;
 	asg_t *sg = *sg_ptr;
@@ -26330,52 +25971,50 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
 	*sg_ptr = sg;
 }
 
-void build_string_graph_without_clean(
-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 max_hang_length, long long clean_round, long long gap_fuzz,
-float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, 
-long long bubble_dist, int read_graph, int write)
+void build_string_graph_without_clean(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 max_hang_length, long long clean_round, long long gap_fuzz,
+		float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name,
+		long long bubble_dist, int read_graph, int write)
 {
-    R_to_U ruIndex;
-    init_R_to_U(&ruIndex, n_read);
-    asg_t *sg = NULL;
-    ma_sub_t* coverage_cut = NULL;
+	R_to_U ruIndex;
+	init_R_to_U(&ruIndex, n_read);
+	asg_t *sg = NULL;
+	ma_sub_t* coverage_cut = NULL;
 
-    // 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
-    min_thres = asm_opt.max_short_tip + 1;
+	// debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs", sources, reverse_sources, -1, "clean");
 
-    if (asm_opt.flag & HA_F_VERBOSE_GFA)
-    {
-        if(load_debug_graph(&sg, &sources, &coverage_cut, output_file_name, &reverse_sources, &ruIndex))
-        {
-            fprintf(stderr, "debug gfa has been loaded\n");
-            
-            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, 
-            output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 1);
-            asg_destroy(sg);
-            free(coverage_cut);
-            destory_R_to_U(&ruIndex);
-            return;
-        }
-    }
-    
-    if (asm_opt.write_index_to_disk && write)
-    {
-        write_all_data_to_disk(sources, reverse_sources, 
-        &R_INF, output_file_name);
-    }
+	///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;
 
-    try_rescue_overlaps(sources, reverse_sources, n_read, 4); 
+	if (asm_opt.flag & HA_F_VERBOSE_GFA)
+	{
+		if (load_debug_graph(&sg, &sources, &coverage_cut, output_file_name, &reverse_sources, &ruIndex))
+		{
+			fprintf(stderr, "debug gfa has been loaded\n");
 
-    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, 
-    output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 0);
-    
-    asg_destroy(sg);
-    free(coverage_cut);
-    destory_R_to_U(&ruIndex);
+			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,
+					output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 1);
+			asg_destroy(sg);
+			free(coverage_cut);
+			destory_R_to_U(&ruIndex);
+			return;
+		}
+	}
+
+	if (asm_opt.write_index_to_disk && write)
+	{
+		write_all_data_to_disk(sources, reverse_sources, &R_INF, output_file_name);
+	}
+
+	try_rescue_overlaps(sources, reverse_sources, n_read, 4);
+
+	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,
+			output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 0);
+
+	asg_destroy(sg);
+	free(coverage_cut);
+	destory_R_to_U(&ruIndex);
 }

Overlaps.h

@@ -103,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
@@ -112,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 {
@@ -128,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;
@@ -505,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

extract.cpp

@@ -3,21 +3,19 @@
 #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)
 
-#define GFA_MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr))))
-#define GFA_REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr))))
-
 char *gfa_strdup(const char *src)
 {
 	int32_t len;
 	char *dst;
 	len = strlen(src);
-	GFA_MALLOC(dst, len + 1);
+	MALLOC(dst, len + 1);
 	memcpy(dst, src, len + 1);
 	return dst;
 }
@@ -25,7 +23,7 @@ char *gfa_strdup(const char *src)
 char *gfa_strndup(const char *src, size_t n)
 {
 	char *dst;
-	GFA_MALLOC(dst, n + 1);
+	MALLOC(dst, n + 1);
 	strncpy(dst, src, n);
 	dst[n] = 0;
 	return dst;
@@ -42,7 +40,7 @@ char **gv_read_list(const char *o, int *n_)
 			if (*p == ',' || *p == 0) {
 				if (n == m) {
 					m = m? m<<1 : 16;
-					GFA_REALLOC(s, m);
+					REALLOC(s, m);
 				}
 				s[n++] = gfa_strndup(q, p - q);
 				if (*p == 0) break;
@@ -63,7 +61,7 @@ char **gv_read_list(const char *o, int *n_)
 			for (p = str.s; *p && !isspace(*p); ++p);
 			if (n == m) {
 				m = m? m<<1 : 16;
-				GFA_REALLOC(s, m);
+				REALLOC(s, m);
 			}
 			s[n++] = gfa_strndup(str.s, p - str.s);
 		}
@@ -94,7 +92,7 @@ void ha_extract_print(const All_reads *rs, int n_rounds, int n, char **list)
 		for (j = 0; j < rs->total_reads; ++j)
 			if (max_len < (int)Get_NAME_LENGTH(*rs, j))
 				max_len = Get_NAME_LENGTH(*rs, j);
-		GFA_MALLOC(s, max_len + 1);
+		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));

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