Release hifiasm-0.8 (r279)

Assembly.cpp

@@ -10,7 +10,6 @@
 #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);
 
@@ -1202,9 +1201,10 @@ int ha_assemble(void)
 		if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
 		ha_triobin(&asm_opt);
 	}
-	build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf, 
-			R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round, 
-			asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded);
+
+    build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf, 
+        R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round, 
+        asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded);
 	destory_All_reads(&R_INF);
 	return 0;
 }

CommandLines.cpp

@@ -19,7 +19,6 @@ 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 }
 };
 
@@ -51,6 +50,7 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "    -n INT        small removed unitig threshold [%d]\n", asm_opt->max_short_tip);
     fprintf(stderr, "    -x FLOAT      max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate);
     fprintf(stderr, "    -y FLOAT      min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate);
+    fprintf(stderr, "    -u            disable post join contigs step which may improve N50. Don't disable in default.\n");
     fprintf(stderr, "    --version     show version number\n");
     fprintf(stderr, "    -h            show help information\n");
 
@@ -362,7 +362,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 
     int c;
 
-    while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:w:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:1:2:3:4:l:s:O:", long_options)) >= 0) {
+    while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:w:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:1:2:3:4:l:s:O:eu", long_options)) >= 0) {
         if (c == 'h')
         {
             Print_H(asm_opt);
@@ -397,6 +397,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 'p') asm_opt->small_pop_bubble_size = atoll(opt.arg);
         else if (c == 'm') asm_opt->large_pop_bubble_size = atoll(opt.arg);
         else if (c == 'n') asm_opt->max_short_tip = atoll(opt.arg);
+        else if (c == 'e') asm_opt->flag |= HA_F_BAN_ASSEMBLY;
+        else if (c == 'u') asm_opt->flag |= HA_F_BAN_POST_JOIN;
 		else if (c == 301) asm_opt->flag |= HA_F_VERBOSE_GFA;
 		else if (c == 302) asm_opt->flag |= HA_F_WRITE_PAF;
 		else if (c == 303) asm_opt->flag |= HA_F_WRITE_EC;
@@ -405,7 +407,6 @@ 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-r257"
+#define HA_VERSION "0.8-r279"
 
 #define VERBOSE 0
 
@@ -15,7 +15,8 @@
 #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_F_BAN_POST_JOIN   0x100
+#define HA_F_BAN_ASSEMBLY    0x200
 
 #define HA_MIN_OV_DIFF       0.02 // min sequence divergence in an overlap
 

Hash_Table.cpp

@@ -4,7 +4,6 @@
 #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 overlap2.o sys.o
+			htab.o hist.o sketch.o anchor.o extract.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 utils.h
+Assembly.o: kthread.h
 CommandLines.o: CommandLines.h ketopt.h
 Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h
 Correct.o: kdq.h CommandLines.h Levenshtein_distance.h POA.h Assembly.h
 Hash_Table.o: Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h
-Hash_Table.o: CommandLines.h ksort.h utils.h
+Hash_Table.o: CommandLines.h ksort.h
 Levenshtein_distance.o: Levenshtein_distance.h
 Output.o: Output.h CommandLines.h
 Overlaps.o: Overlaps.h kvec.h kdq.h ksort.h Process_Read.h CommandLines.h
 Overlaps.o: Hash_Table.h htab.h Correct.h Levenshtein_distance.h POA.h
-Overlaps.o: Purge_Dups.h utils.h
+Overlaps.o: Purge_Dups.h
 POA.o: POA.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h
 POA.o: CommandLines.h Correct.h Levenshtein_distance.h
 Process_Read.o: Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
@@ -52,17 +52,16 @@ 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 utils.h
+Trio.o: CommandLines.h htab.h
 anchor.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
-anchor.o: ksort.h utils.h Hash_Table.h
+anchor.o: ksort.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
+extract.o: kseq.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 utils.h
+htab.o: kvec.h kdq.h CommandLines.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 utils.h
-overlap2.o: utils.h CommandLines.h Overlaps.h kvec.h kdq.h
+main.o: Levenshtein_distance.h htab.h
 sketch.o: kvec.h htab.h Process_Read.h Overlaps.h kdq.h CommandLines.h
-sys.o: utils.h
+sys.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h

Overlaps.cpp

@@ -9,7 +9,6 @@
 #include "Hash_Table.h"
 #include "Correct.h"
 #include "Purge_Dups.h"
-#include "utils.h"
 
 uint32_t debug_purge_dup = 0;
 
@@ -77,17 +76,21 @@ 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)
+    {
+        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);
+        }
 
-	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);
+        free(g->F_seq);
+    }
+    
+
+    free(g);
 }
 
 void asg_arc_sort(asg_t *g)
@@ -129,7 +132,8 @@ 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;
@@ -182,6 +186,7 @@ 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)
@@ -191,6 +196,7 @@ 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;
 }
 
@@ -270,15 +276,16 @@ 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
@@ -865,6 +872,51 @@ 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)
 {
@@ -975,6 +1027,108 @@ 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();
@@ -1540,7 +1694,9 @@ 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();
@@ -1952,6 +2108,71 @@ 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;
+}
 
 
 
@@ -4263,6 +4484,143 @@ 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)
 {
@@ -12893,16 +13251,19 @@ kvec_asg_arc_t_warp* new_rtg_edges)
 
     renew_utg(ug, read_g, new_rtg_edges);
 
-    rescue_missing_overlaps_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang,
-    min_ovlp, 0, 0, 1, NULL);
+    if (!(asm_opt.flag & HA_F_BAN_POST_JOIN))
+    {
+        rescue_missing_overlaps_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang,
+        min_ovlp, 0, 0, 1, NULL);
 
-    renew_utg(ug, read_g, new_rtg_edges);
+        renew_utg(ug, read_g, new_rtg_edges);
 
-    rescue_contained_reads_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang, 
-    min_ovlp, 0, 10, 0, 1, NULL, NULL);
-
-    renew_utg(ug, read_g, new_rtg_edges);
+        rescue_contained_reads_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang, 
+        min_ovlp, 0, 10, 0, 1, NULL, NULL);
 
+        renew_utg(ug, read_g, new_rtg_edges);
+    }
+    
     update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, 1, flag, drop_rate);
 
     update_hap_label(NULL, read_g);
@@ -21726,27 +22087,27 @@ kvec_asg_arc_t_warp* new_rtg_edges)
         renew_utg(ug, read_g, new_rtg_edges);
     }
     
-
-
-    rescue_missing_overlaps_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang,
-    min_ovlp, 0, 0, 1, NULL);
-    renew_utg(ug, read_g, new_rtg_edges);
-
-    rescue_contained_reads_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang, 
-    min_ovlp, 0, 10, 0, 1, NULL, NULL);
-    renew_utg(ug, read_g, new_rtg_edges);
-
-
-    if(asm_opt.purge_level_primary > 0)
+    if (!(asm_opt.flag & HA_F_BAN_POST_JOIN))
     {
-        just_contain = 0;
-        if(asm_opt.purge_level_primary == 1) just_contain = 1;
-
-        purge_dups(*ug, read_g, coverage_cut, reverse_sources, ruIndex, new_rtg_edges, 
-        asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, bubble_dist, 
-        drop_ratio, just_contain);
-        delete_useless_nodes(ug);
+        rescue_missing_overlaps_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang,
+        min_ovlp, 0, 0, 1, NULL);
         renew_utg(ug, read_g, new_rtg_edges);
+
+        rescue_contained_reads_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang, 
+        min_ovlp, 0, 10, 0, 1, NULL, NULL);
+        renew_utg(ug, read_g, new_rtg_edges);
+
+        if(asm_opt.purge_level_primary > 0)
+        {
+            just_contain = 0;
+            if(asm_opt.purge_level_primary == 1) just_contain = 1;
+
+            purge_dups(*ug, read_g, coverage_cut, reverse_sources, ruIndex, new_rtg_edges, 
+            asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, bubble_dist, 
+            drop_ratio, just_contain);
+            delete_useless_nodes(ug);
+            renew_utg(ug, read_g, new_rtg_edges);
+        }
     }
     
 
@@ -25728,12 +26089,13 @@ 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;
@@ -25748,7 +26110,7 @@ void clean_graph(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_so
     normalize_ma_hit_t_single_side_advance(sources, n_read);
     normalize_ma_hit_t_single_side_advance(reverse_sources, n_read);
     
-    // debug_info_of_specfic_read(">m64011_190830_220126/117834372/ccs", sources, reverse_sources, 
+    // debug_info_of_specfic_read("m64062_190803_042216/122882911/ccs", sources, reverse_sources, 
     // -1, "clean");
 
 
@@ -25761,7 +26123,6 @@ void clean_graph(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_so
         memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t));
     }
     
-
     ///print_binned_reads(sources, n_read, coverage_cut);
 
     clean_weak_ma_hit_t(sources, reverse_sources, n_read);
@@ -25774,8 +26135,6 @@ void clean_graph(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_so
     ma_hit_flt(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
 
     ///fix_binned_reads(sources, n_read, coverage_cut);
-
-
     ///just need to deal with trio here
     ma_hit_contained_advance(sources, n_read, coverage_cut, ruIndex, max_hang_length, mini_overlap_length);
 
@@ -25783,7 +26142,6 @@ void clean_graph(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_so
     asg_arc_del_trans(sg, gap_fuzz);
     asm_opt.coverage = get_coverage(sources, coverage_cut, n_read);
 
-
     if(VERBOSE >= 1)
     {
         char* unlean_name = (char*)malloc(strlen(output_file_name)+25);
@@ -25793,10 +26151,9 @@ void clean_graph(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_so
     }
 
     asg_cut_tip(sg, asm_opt.max_short_tip);
-
     
     ///drop_inexact_edegs_at_bubbles(sg, bubble_dist);
-    
+
     if(clean_round > 0)
     {
         double cut_step;
@@ -25812,6 +26169,7 @@ void clean_graph(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_so
         int i = 0;
         for (i = 0; i < clean_round; i++, drop_ratio += cut_step)
         {
+
             if(drop_ratio > max_ovlp_drop_ratio)
             {
                 drop_ratio = max_ovlp_drop_ratio;
@@ -25920,7 +26278,7 @@ void clean_graph(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_so
     debug_ma_hit_t(reverse_sources, coverage_cut, n_read, max_hang_length, 
     mini_overlap_length);
     **/
-
+ 
     ///note: don't apply asg_arc_del_too_short_overlaps() after this function!!!!
     rescue_contained_reads_aggressive(NULL, sg, sources, coverage_cut, ruIndex, max_hang_length, 
     mini_overlap_length, bubble_dist, 10, 1, 0, NULL, NULL);
@@ -25971,50 +26329,56 @@ void clean_graph(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_so
 	*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");
+    // 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;
 
-	///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;
+    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);
+    }
 
-	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");
+    if (!(asm_opt.flag & HA_F_BAN_ASSEMBLY))
+    {
+        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, 1);
-			asg_destroy(sg);
-			free(coverage_cut);
-			destory_R_to_U(&ruIndex);
-			return;
-		}
-	}
+        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);
+    }
 
-	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);
+    destory_R_to_U(&ruIndex);
 }

Overlaps.h

@@ -103,6 +103,7 @@ 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,10 +112,11 @@ 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 {
@@ -126,10 +128,6 @@ 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;
@@ -507,25 +505,6 @@ 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)
 {

README.md

@@ -4,18 +4,79 @@
 # Install hifiasm (requiring g++ and zlib)
 git clone https://github.com/chhylp123/hifiasm
 cd hifiasm && make
-# Assembly
-./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz
+
+# Run on test data (use -f0 for small datasets)
+wget https://github.com/chhylp123/hifiasm/releases/download/v0.7/chr11-2M.fa.gz
+./hifiasm -o test -t4 -f0 chr11-2M.fa.gz 2> test.log
+awk '/^S/{print ">"$1;print $2}' test.p_ctg.gfa > test.p_ctg.fa  # get primary contigs in FASTA
+
+# Assemble inbred/homozygous genomes (-l0 disables duplication purging)
+hifiasm -o CHM13.asm -t32 -l0 CHM13-HiFi.fa.gz 2> CHM13.asm.log
+# Assemble heterozygous with built-in duplication purging
+hifiasm -o HG002.asm -t32 HG002-file1.fq.gz HG002-file2.fq.gz
+
+# Trio binning assembly (requiring https://github.com/lh3/yak)
+yak count -b37 -t16 -o pat.yak <(cat pat_1.fq.gz pat_2.fq.gz) <(cat pat_1.fq.gz pat_2.fq.gz)
+yak count -b37 -t16 -o mat.yak <(cat mat_1.fq.gz mat_2.fq.gz) <(cat mat_1.fq.gz mat_2.fq.gz)
+hifiasm -o HG002.asm -t32 -1 pat.yak -2 mat.yak HG002-HiFi.fa.gz
 ```
 
 ## Introduction
 
-Hifiasm is a fast haplotype-resolved de novo assembler for PacBio
-Hifi reads. Unlike most existing assemblers, hifiasm starts from uncollapsed
-genome. Thus, it is able to keep the haplotype information as much as possible.
+Hifiasm is a fast haplotype-resolved de novo assembler for PacBio Hifi reads. 
+It can assemble a human genome in several hours and works with the California
+redwood genome, one of the most complex genomes sequenced so far. Hifiasm can
+produce primary/alternate assemblies of quality competitive with the best
+assemblers. It also introduces a new graph binning algorithm and achieves
+the best haplotype-resolved assembly given trio data.
 
-For non-trio assembly, the input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and its
-outputs consist of: 
+## Usage
+
+A typical hifiasm command line looks like:
+```sh
+hifiasm -o NA12878.asm -t 32 NA12878.fq.gz
+```
+where `NA12878.fq.gz` provides the input reads, `-t` sets the number of CPUs in
+use and `-o` specifies the prefix of output files. For this example, the
+primary contigs are written to `NA12878.asm.p_ctg.gfa` and alternate contigs to
+`NA12878.asm.a_ctg.gfa`. At the first run, hifiasm saves corrected reads and
+overlaps to disk as `NA12878.asm.*.bin`. It reuses the saved results to avoid
+the time-consuming all-vs-all overlap calculation next time. You may specify
+`-i` to ignore precomputed overlaps and redo overlapping from raw reads.
+
+Hifiasm purges haplotig duplications by default. For inbred or homozygous
+genomes, you may disable purging with option `-l0`. Old HiFi reads may contain
+short adapter sequences at the ends of reads. You can specify `-z20` to trim
+both ends of reads by 20bp. For small genomes, use `-f0` to disable the initial
+bloom filter which takes 16GB memory at the beginning. For genomes much larger
+than human, applying `-f38` or even `-f39` is preferred to save memory on k-mer
+counting.
+
+When parental short reads are available, hifiasm can generate a pair of
+haplotype-resolved assemblies with trio binning. To perform such assembly, you
+need to count k-mers first with [yak][yak] first and then do assembly:
+```sh
+yak count -k31 -b37 -t16 -o pat.yak paternal.fq.gz
+yak count -k31 -b37 -t16 -o mat.yak maternal.fq.gz
+hifiasm -o NA12878.asm -t 32 -1 pat.yak -2 mat.yak NA12878.fq.gz
+```
+Here `NA12878.asm.hap1.p_ctg.gfa` and `NA12878.asm.hap2.p_ctg.gfa` give the two
+haplotype assemblies. In the binning mode, hifiasm does not purge haplotig
+duplications by default. Because hifiasm reuses saved overlaps, you can
+generate both primary/alternate assemblies and trio binning assemblies with
+```sh
+hifiasm -o NA12878.asm -t 32 NA12878.fq.gz 2> NA12878.asm.pri.log
+hifiasm -o NA12878.asm -t 32 -1 pat.yak -2 mat.yak /dev/null 2> NA12878.asm.trio.log
+```
+The second command line will run much faster than the first. You can also dump
+error corrected in FASTA and/or overlaps in PAF with
+```sh
+hifiasm -o NA12878.asm -t 32 --write-paf --write-ec /dev/null
+```
+
+## Output files
+
+For non-trio assembly, hifiasm generates the following files:
 
 1. Haplotype-resolved raw [unitig][unitig] graph in [GFA][gfa] format
    (*prefix*.r\_utg.gfa). This graph keeps all haplotype information, including
@@ -28,7 +89,8 @@ outputs consist of:
 4. Alternate assembly contig graph (*prefix*.a\_ctg.gfa). This graph consists of all assemblies that
    are discarded in primary contig graph.
 
-For trio assembly, the input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and the paternal/maternal trio indexes generated by `yak count` (see https://github.com/lh3/yak). The outputs consist of:
+For trio assembly, hifiasm generates the following files:
+
 1. Haplotype-resolved raw [unitig][unitig] graph in [GFA][gfa] format
    (*prefix*.r\_utg.gfa). This graph keeps all haplotype information. 
 
@@ -38,77 +100,73 @@ For trio assembly, the input of hifiasm is the PacBio Hifi reads in fasta/fastq
 3. Phased maternal/haplotype2 contig graph (*prefix*.hap2.p\_ctg.gfa). This graph keeps the phased
    maternal/haplotype2 assembly.
 
+Hifiasm writes error corrected reads to the *prefix*.ec.bin binary file and
+writes overlaps to *prefix*.ovlp.source.bin and *prefix*.ovlp.reverse.bin.
 
+## Results
 
-In addition, hifiasm also outputs three binary files that save all overlap information (*prefix*.ec.bin, *prefix*.ovlp.reverse.bin, *prefix*.ovlp.source.bin). With these files, hifiasm can avoid the time-consuming all-to-all overlap calculation step, and do the assembly
-directly and quickly. This might be helpful when you want to get an optimized
-assembly by multiple rounds of experiments with different parameters.
+The following table shows the statistics of several hifiasm primary assemblies:
 
-Hifiasm is a standalone and lightweight assembler, which does not need external
-libraries (except zlib). For large genomes, it can generate high-quality
-assembly in a few hours. Hifiasm has been tested on human, butterfly, rice and drosophila. 
-In particular, hifiasm is able to assemble the 26.5Gb California redwood tree in a few days. 
-The results are as follows: 
-
-|<sub>Dataset<sub>|<sub>GSize<sub>|<sub>Cov<sub>|<sub>Asm options<sub>|<sub>CPU time<sub>|<sub>Wall time<sub>|<sub>RAM<sub>|<sub>[unitig][unitig]/[contig][unitig] N50<sup>[1]</sup><sub>|
+|<sub>Dataset<sub>|<sub>Size<sub>|<sub>Cov.<sub>|<sub>Asm options<sub>|<sub>CPU time<sub>|<sub>Wall time<sub>|<sub>RAM<sub>|<sub> N50<sub>|
 |:---------------|-----:|-----:|:---------------------|-------:|--------:|----:|----------------:|
-|<sub>[\[Redwood\]](https://downloads.pacbcloud.com/public/dataset/redwood2020/)<sub>|<sub>26.5Gb<sub>|<sub>x23<sub>|<sub>-k 40 -t 64 -r 2<sub>|<sub>7274h30m<sub>|<sub>141h30m<sub>|<sub>512G<sub>|<sub>1.7Mb/1.9Mb<sub>|
+|<sub>[Mouse (C57/BL6J)][mouse-data]</sub>|<sub>2.6Gb</sub> |<sub>&times;25</sub>|<sub>-t48 -l0</sub> |<sub>172.9h</sub> |<sub>4.8h</sub> |<sub>76G</sub> |<sub>21.1Mb</sub>|
+|<sub>[Maize (B73)][maize-data]</sub>     |<sub>2.2Gb</sub> |<sub>&times;22</sub>|<sub>-t48 -l0</sub> |<sub>203.2h</sub> |<sub>5.1h</sub> |<sub>68G</sub> |<sub>36.7Mb</sub>|
+|<sub>[Strawberry][strawberry-data]</sub> |<sub>0.8Gb</sub> |<sub>&times;36</sub>|<sub>-t48 -D10</sub>|<sub>152.7h</sub> |<sub>3.7h</sub> |<sub>91G</sub> |<sub>17.8Mb</sub>|
+|<sub>[Frog][frog-data]</sub>             |<sub>9.5Gb</sub> |<sub>&times;29</sub>|<sub>-t48</sub>     |<sub>2834.3h</sub>|<sub>69.0h</sub>|<sub>463G</sub>|<sub>9.3Mb</sub>|
+|<sub>[Redwood][redwood-data]</sub>       |<sub>35.6Gb</sub>|<sub>&times;28</sub>|<sub>-t80</sub>     |<sub>3890.3h</sub>|<sub>65.5h</sub>|<sub>699G</sub>|<sub>5.4Mb</sub>|
+|<sub>[Human (CHM13)][CHM13-data]</sub>   |<sub>3.1Gb</sub> |<sub>&times;32</sub>|<sub>-t48 -l0</sub> |<sub>310.7h</sub> |<sub>8.2h</sub> |<sub>114G</sub>|<sub>88.9Mb</sub>|
+|<sub>[Human (HG00733)][HG00733-data]</sub>|<sub>3.1Gb</sub>|<sub>&times;33</sub>|<sub>-t48</sub>     |<sub>269.1h</sub> |<sub>6.9h</sub> |<sub>135G</sub>|<sub>69.9Mb</sub>|
+|<sub>[Human (HG002)][NA24385-data]</sub> |<sub>3.1Gb</sub> |<sub>&times;36</sub>|<sub>-t48</sub>     |<sub>305.4h</sub> |<sub>7.7h</sub> |<sub>137G</sub>|<sub>98.7Mb</sub>|
 
-<sub>[1] unitig N50 is the N50 of assembly graph with haplotype information (i.e., bubbles), while the contig N50 is the N50 of haplotype collapsed assembly (i.e., without bubbles).<sub>
+[mouse-data]:      https://www.ncbi.nlm.nih.gov/sra/?term=SRR11606870
+[maize-data]:      https://www.ncbi.nlm.nih.gov/sra/?term=SRR11606869
+[strawberry-data]: https://www.ncbi.nlm.nih.gov/sra/?term=SRR11606867
+[frog-data]:       https://www.ncbi.nlm.nih.gov/sra?term=(SRR11606868)%20OR%20SRR12048570
+[redwood-data]:    https://www.ncbi.nlm.nih.gov/sra/?term=SRP251156
+[CHM13-data]:      https://www.ncbi.nlm.nih.gov/sra?term=(((SRR11292120)%20OR%20SRR11292121)%20OR%20SRR11292122)%20OR%20SRR11292123
 
-## Usage
+Hifiasm can assemble a 3.1Gb human genome in several hours or a ~30Gb hexaploid
+redwood genome in a few days on a single machine. For trio binning assembly:
 
-For Hifi reads assembly, a typical command line looks like:
+|<sub>Dataset<sub>|<sub>Cov.<sub>|<sub>CPU time<sub>|<sub>Elapsed time<sub>|<sub>RAM<sub>|<sub> N50<sub>|
+|:---------------|-----:|-------:|--------:|----:|----------------:|
+|<sub>[HG00733][HG00733-data], [\[father\]][HG00731-data], [\[mother\]][HG00732-data]</sub>|<sub>&times;33</sub>|<sub>269.1h</sub>|<sub>6.9h</sub>|<sub>135G</sub>|<sub>35.1Mb (paternal), 34.9Mb (maternal)</sub>|
+|<sub>[HG002][NA24385-data],   [\[father\]][NA24149-data], [\[mother\]][NA24143-data]</sup>|<sub>&times;36</sub>|<sub>305.4h</sub>|<sub>7.7h</sub>|<sub>137G</sub>|<sub>41.0Mb (paternal), 40.8Mb (maternal)</sub>|
+|<sub>[NA12878][NA12878-data], [\[father\]][NA12891-data], [\[mother\]][NA12892-data]</sub>|<sub>&times;30</sub>|<sub>180.8h</sub>|<sub>4.9h</sub>|<sub>123G</sub>|<sub>27.7Mb (paternal), 27.0Mb (maternal)</sub>|
 
-```sh
-./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz
+[HG00733-data]: https://www.ebi.ac.uk/ena/data/view/ERX3831682
+[HG00731-data]: https://www.ebi.ac.uk/ena/data/view/ERR3241754
+[HG00732-data]: https://www.ebi.ac.uk/ena/data/view/ERR3241755
+[NA24385-data]: https://www.ncbi.nlm.nih.gov/sra?term=(((SRR10382244)%20OR%20SRR10382245)%20OR%20SRR10382248)%20OR%20SRR10382249
+[NA24149-data]: https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/HG003_NA24149_father/NIST_HiSeq_HG003_Homogeneity-12389378/HG003Run01-13262252/
+[NA24143-data]: https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/AshkenazimTrio/HG004_NA24143_mother/NIST_HiSeq_HG004_Homogeneity-14572558/HG004Run01-15133132/
+[NA12878-data]: https://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/NA12878/PacBio_SequelII_CCS_11kb/
+[NA12891-data]: https://www.ebi.ac.uk/ena/data/view/ERR194160
+[NA12892-data]: https://www.ebi.ac.uk/ena/data/view/ERR194161
+
+Except NA12878, the assemblies above were produced by hifiasm v0.7 and can be
+downloaded at
+```txt
+ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm/submission/v0.7/
+```
+NA12878 was assembled with a more recent version of hifiasm and is available at
+```txt
+ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm/NA12878-r253/
 ```
 
-where `NA12878.fq.gz` is the input reads and `-o` specifies the output files.
-In this example, all output files can be found at `NA12878.asm.*`. `-t` specifies 
-the number of CPU threads. Note that at first run, hifiasm will save all overlaps 
-to disk, which can avoid the time-consuming all-to-all overlap calculation next time. 
-For hifiasm, once the overlap information has been obtained during the previous run 
-in advance, it is able to load all overlaps from disk and then directly do assembly. 
-If you want to ignore the pre-computed overlap information, please specify `-i`.
-
-Please note that some old Hifi reads may consist of short adapters. To improve
-the assembly quality, adapters should be removed by `-z` as follow:
-
-```sh
-./hifiasm -o butterfly.asm -t 42 -z 20 butterfly.fq.gz
-```
-
-In this example, hifiasm will remove 20 bases from both ends of each read.
-
-For trio assembly, first the trio indexes of paternal/maternal should be generated by 
-`yak count` (see https://github.com/lh3/yak):
-
-```sh
-./yak count -k31 -b37 -t16 -o mat.yak mat.fq.gz
-```
-```sh
-./yak count -k31 -b37 -t16 -o pat.yak pat.fq.gz
-```
-
-and then run hifiasm as follow:
-
-```sh
-./hifiasm -o NA12878.asm -t 32 -1 pat.yak -2 mat.yak NA12878_1.fq.gz NA12878_2.fq.gz
-```
 
 [unitig]: http://wgs-assembler.sourceforge.net/wiki/index.php/Celera_Assembler_Terminology
 [gfa]: https://github.com/pmelsted/GFA-spec/blob/master/GFA-spec.md
 [paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
+[yak]: https://github.com/lh3/yak
 
 ## Getting Help
 
-For detailed description of options, please see `man ./hifiasm.1`.
-The `-h` option of hifiasm also provides simple description of options. If you
-have further questions, please raise an issue at the issue page.
+For detailed description of options, please see `man ./hifiasm.1`. The `-h`
+option of hifiasm also provides brief description of options. If you have
+further questions, please raise an issue at the [issue
+page](https://github.com/chhylp123/hifiasm/issues).
 
-## Limitations and future works
+## Limitations
 
-1. The running time and memory usage should be further reduced.
-
-2. The N50 should be further improved. 
+1. Purging haplotig duplications may introduce misassemblies.

Trio.cpp

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

anchor.cpp

@@ -1,7 +1,6 @@
 #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,19 +3,21 @@
 #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);
-	MALLOC(dst, len + 1);
+	GFA_MALLOC(dst, len + 1);
 	memcpy(dst, src, len + 1);
 	return dst;
 }
@@ -23,7 +25,7 @@ char *gfa_strdup(const char *src)
 char *gfa_strndup(const char *src, size_t n)
 {
 	char *dst;
-	MALLOC(dst, n + 1);
+	GFA_MALLOC(dst, n + 1);
 	strncpy(dst, src, n);
 	dst[n] = 0;
 	return dst;
@@ -40,7 +42,7 @@ char **gv_read_list(const char *o, int *n_)
 			if (*p == ',' || *p == 0) {
 				if (n == m) {
 					m = m? m<<1 : 16;
-					REALLOC(s, m);
+					GFA_REALLOC(s, m);
 				}
 				s[n++] = gfa_strndup(q, p - q);
 				if (*p == 0) break;
@@ -61,7 +63,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;
-				REALLOC(s, m);
+				GFA_REALLOC(s, m);
 			}
 			s[n++] = gfa_strndup(str.s, p - str.s);
 		}
@@ -92,7 +94,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);
-		MALLOC(s, max_len + 1);
+		GFA_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));

hifiasm.1

@@ -1,4 +1,4 @@
-.TH hifiasm 1 "12 Apr 2020" "hifiasm-0.5.0" "Bioinformatics tools"
+.TH hifiasm 1 "27 June 2020" "hifiasm-0.8 (r279)" "Bioinformatics tools"
 
 .SH NAME
 .PP
@@ -322,5 +322,3 @@ For each graph, hifiasm also outputs a simplified version without sequences for
 the ease of visualization. Hifiasm keeps corrected reads and overlaps in three
 binary files such as it can regenerate assembly graphs from the binary files
 without redoing error correction.
-
-

htab.cpp

@@ -9,7 +9,6 @@
 #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,6 +38,13 @@ 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);
@@ -73,6 +80,18 @@ 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 "utils.h"
+#include "htab.h"
 
 int main(int argc, char *argv[])
 {

overlap2.cpp

@@ -1,407 +0,0 @@
-#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 "utils.h"
+#include "htab.h"
 
 int yak_verbose = 3;
 

utils.h

@@ -1,33 +0,0 @@
-#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