r257: dropped contained unitigs

prepare for further changes
2026-09-16 05:07:55 +08:00 · 2020-06-08 15:43:21 -04:00 · 2020-06-06 21:55:41 -04:00
23 changed files with 767 additions and 2022 deletions
--- a/Assembly.cpp
+++ b/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);

@@ -476,7 +477,6 @@ static void worker_ovec(void *data, long i, int tid)
 	int fully_cov, abnormal;

 	ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1);
-    ///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 1);

 	clear_Cigar_record(&b->cigar1);
 	clear_Round2_alignment(&b->round2);
@@ -515,7 +515,6 @@ static void worker_ovec_related_reads(void *data, long i, int tid)
 			memcmp(asm_opt.required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
 	{
 		ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1);
-        ///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 1);

 		fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i));
 		recover_UC_Read(&b->self_read, &R_INF, i);
@@ -634,7 +633,7 @@ void Output_corrected_reads()

 void ha_overlap_and_correct(int round)
 {
-	int i, hom_cov, het_cov;
+	int i, hom_cov;
 	ha_ovec_buf_t **b;
 	ha_ecsave_buf_t *e;

@@ -642,7 +641,7 @@ void ha_overlap_and_correct(int round)
 	CALLOC(b, asm_opt.thread_num);
 	for (i = 0; i < asm_opt.thread_num; ++i)
 		b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1));
-	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF, &hom_cov, &het_cov); // build the index
+	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF, &hom_cov); // build the index
 	if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
 		ha_opt_update_cov(&asm_opt, hom_cov);
 	if (asm_opt.required_read_name)
@@ -959,7 +958,6 @@ static void worker_ov_final(void *data, long i, int tid)

 	//get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
 	ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.001, asm_opt.max_n_chain, 0);
-    ///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 0);

 	/**
 	  correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon,
@@ -1144,55 +1142,20 @@ long long readNum, long long rescue_threshold, float cluster_threshold)
    fprintf(stderr, "[M::%s] took %0.2fs, revise edges #: %lld\n\n", __func__, Get_T()-startTime, revises);
 }

-void hap_recalculate_peaks(char* output_file_name)
-{
-    destory_read_bin(&R_INF);
-    destory_ma_hit_t_alloc(R_INF.paf);
-    destory_ma_hit_t_alloc(R_INF.reverse_paf);
-
-    char* gfa_name = (char*)malloc(strlen(output_file_name)+25);
-	sprintf(gfa_name, "%s.ec", output_file_name);
-
-    int hom_cov, het_cov;
-    // construct hash table for high occurrence k-mers
-    if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) {
-        ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov);
-        ha_opt_update_cov(&asm_opt, hom_cov);
-    }
-    free(R_INF.read_length);
-    free(R_INF.name_index);
-
-    load_All_reads(&R_INF, gfa_name);
-
-    ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index
-    asm_opt.hom_cov = hom_cov;
-    asm_opt.het_cov = het_cov;
-    ha_pt_destroy(ha_idx);
-    ha_idx = 0;
-
-    destory_read_bin(&R_INF);
-    free(gfa_name);
-
-    load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name); 
-    fprintf(stderr, "M::%s has done.\n", __func__);   
-}
-
 void ha_overlap_final(void)
 {
-	int i, hom_cov, het_cov;
+	int i, hom_cov;
 	ha_ovec_buf_t **b;
 	CALLOC(b, asm_opt.thread_num);
 	for (i = 0; i < asm_opt.thread_num; ++i)
 		b[i] = ha_ovec_init(1, 1);
-	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index
+	ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov); // build the index
 	kt_for(asm_opt.thread_num, worker_ov_final, b, R_INF.total_reads);
 	ha_pt_destroy(ha_idx);
 	ha_idx = 0;
 	for (i = 0; i < asm_opt.thread_num; ++i)
 		ha_ovec_destroy(b[i]);
 	free(b);
-    asm_opt.hom_cov = hom_cov;
-    asm_opt.het_cov = het_cov;
 }

 int ha_assemble(void)
@@ -1207,10 +1170,9 @@ int ha_assemble(void)
 			exit(0);
 		}
 		if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN) && !(asm_opt.flag & HA_F_VERBOSE_GFA)) ha_triobin(&asm_opt);
-        ///if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt), ovlp_loaded = 2;
+        ///if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt);
 		if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();
 		if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
-        if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2;
 	}
 	if (!ovlp_loaded) {
 		// construct hash table for high occurrence k-mers
@@ -1240,11 +1202,9 @@ int ha_assemble(void)
 		if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
 		ha_triobin(&asm_opt);
 	}
-    if(ovlp_loaded == 2) ovlp_loaded = 0;
-
-    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;
 }
--- a/CommandLines.cpp
+++ b/CommandLines.cpp
@@ -1,8 +1,6 @@
 #include <zlib.h>
 #include <stdlib.h>
 #include <stdio.h>
-#include <stdint.h>
-#include <limits.h>
 #include <sys/time.h>
 #include "CommandLines.h"
 #include "ketopt.h"
@@ -21,8 +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 },
-    { "purge-cov",     ko_required_argument, 309 },
-    { "pri-range",     ko_required_argument, 310 },
+	{ "keep-contained",ko_no_argument, 309 },
 	{ 0, 0, 0 }
 };

@@ -37,46 +34,42 @@ void Print_H(hifiasm_opt_t* asm_opt)
 {
    fprintf(stderr, "Usage: hifiasm [options] <in_1.fq> <in_2.fq> <...>\n");
    fprintf(stderr, "Options:\n");
-	fprintf(stderr, "  Input/Output:\n");
-    fprintf(stderr, "    -o STR      prefix of output files [%s]\n", asm_opt->output_file_name);
-    fprintf(stderr, "    -i          ignore saved read correction and overlaps\n");
-    fprintf(stderr, "    -t INT      number of threads [%d]\n", asm_opt->thread_num);
-    fprintf(stderr, "    -z INT      length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
-    fprintf(stderr, "    --version   show version number\n");
-	fprintf(stderr, "  Overlap/Error correction:\n");
-    fprintf(stderr, "    -k INT      k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length);
-	fprintf(stderr, "    -w INT      minimizer window size [%d]\n", asm_opt->mz_win);
-	fprintf(stderr, "    -f INT      number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift);
-	fprintf(stderr, "    -D FLOAT    drop k-mers occurring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor);
-	fprintf(stderr, "    -N INT      consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
-    fprintf(stderr, "    -r INT      round of correction [%d]\n", asm_opt->number_of_round);
    fprintf(stderr, "  Assembly:\n");
-    fprintf(stderr, "    -a INT      round of assembly cleaning [%d]\n", asm_opt->clean_round);
-    fprintf(stderr, "    -m INT      pop bubbles of <INT in size in contig graphs [%lld]\n", asm_opt->large_pop_bubble_size);
-    fprintf(stderr, "    -p INT      pop bubbles of <INT in size in unitig graphs [%lld]\n", asm_opt->small_pop_bubble_size);
-    fprintf(stderr, "    -n INT      remove tip unitigs composed of <=INT reads [%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\n");
-//	fprintf(stderr, "    --pri-range INT1[,INT2]\n");
-//	fprintf(stderr, "                keep contigs with coverage in this range in p_ctg.gfa; -1 to disable [auto,inf]\n");
+    fprintf(stderr, "    -o FILE       prefix of output files [%s]\n", asm_opt->output_file_name);
+    fprintf(stderr, "    -t INT        number of threads [%d]\n", asm_opt->thread_num);
+    fprintf(stderr, "    -r INT        round of correction [%d]\n", asm_opt->number_of_round);
+    fprintf(stderr, "    -a INT        round of assembly cleaning [%d]\n", asm_opt->clean_round);
+    fprintf(stderr, "    -k INT        k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length);
+	fprintf(stderr, "    -w INT        minimizer window size [%d]\n", asm_opt->mz_win);
+	fprintf(stderr, "    -f INT        number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift);
+	fprintf(stderr, "    -D FLOAT      drop k-mers occuring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor);
+	fprintf(stderr, "    -N INT        consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
+    fprintf(stderr, "    -i            ignore saved overlaps in *.ovlp* files\n");
+    fprintf(stderr, "    -z INT        length of adapters that should be removed [%d]\n", asm_opt->adapterLen);
+    fprintf(stderr, "    -m INT        size of popped large bubbles for contig graph [%lld]\n", asm_opt->large_pop_bubble_size);
+    fprintf(stderr, "    -p INT        size of popped small bubbles for haplotype-resolved unitig graph [%lld]\n", asm_opt->small_pop_bubble_size);
+    fprintf(stderr, "    -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, "    --version     show version number\n");
+    fprintf(stderr, "    -h            show help information\n");

    fprintf(stderr, "  Trio-partition:\n");
-    fprintf(stderr, "    -1 FILE     hap1/paternal k-mer dump generated by \"yak count\" []\n");
-    fprintf(stderr, "    -2 FILE     hap2/maternal k-mer dump generated by \"yak count\" []\n");
-    fprintf(stderr, "    -c INT      lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt);
-    fprintf(stderr, "    -d INT      upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt);
-	fprintf(stderr, "    -3 FILE     list of hap1/paternal read names []\n");
-	fprintf(stderr, "    -4 FILE     list of hap2/maternal read names []\n");
+    fprintf(stderr, "    -1 FILE       hap1/paternal k-mer dump generated by \"yak count\" []\n");
+    fprintf(stderr, "    -2 FILE       hap2/maternal k-mer dump generated by \"yak count\" []\n");
+	fprintf(stderr, "    -3 FILE       list of hap1/paternal read names []\n");
+	fprintf(stderr, "    -4 FILE       list of hap2/maternal read names []\n");
+    fprintf(stderr, "    -c INT        lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt);
+    fprintf(stderr, "    -d INT        upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt);

    fprintf(stderr, "  Purge-dups:\n");
-    fprintf(stderr, "    -l INT      purge level. 0: no purging; 1: light; 2: aggressive [0 for trio; 2 for unzip]\n");
-    fprintf(stderr, "    -s FLOAT    similarity threshold for duplicate haplotigs [%g]\n", 
-                                     asm_opt->purge_simi_rate);
-    fprintf(stderr, "    -O INT      min number of overlapped reads for duplicate haplotigs [%d]\n", 
-                                     asm_opt->purge_overlap_len);
-    fprintf(stderr, "    --purge-cov INT\n");
-    fprintf(stderr, "                coverage upper bound of Purge-dups [auto]\n");
+    fprintf(stderr, "    -l INT        level of purge-dup. In default, [%d] for non-trio; [%d] for trio (see hifiasm.1 for details)\n", 
+                                       asm_opt->purge_level_primary, asm_opt->purge_level_trio);
+    fprintf(stderr, "    -s FLOAT      similarity threshold for duplicate haplotigs [%g]\n", 
+                                       asm_opt->purge_simi_rate);
+    fprintf(stderr, "    -O INT        min number of overlapped reads for duplicate haplotigs [%d]\n", 
+                                       asm_opt->purge_overlap_len);
+

    fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
    fprintf(stderr, "See `man ./hifiasm.1' for detailed description of these command-line options.\n");
@@ -99,7 +92,6 @@ void init_opt(hifiasm_opt_t* asm_opt)
 	asm_opt->max_ov_diff_ec = 0.04;
 	asm_opt->max_ov_diff_final = 0.03;
 	asm_opt->hom_cov = 20;
-    asm_opt->het_cov = -1024;
 	asm_opt->max_n_chain = 100;
    asm_opt->k_mer_min_freq = 3;
    asm_opt->k_mer_max_freq = 66;
@@ -124,9 +116,6 @@ void init_opt(hifiasm_opt_t* asm_opt)
    asm_opt->purge_level_trio = 0;
    asm_opt->purge_simi_rate = 0.75;
    asm_opt->purge_overlap_len = 1;
-    asm_opt->recover_atg_cov_min = -1024;
-    asm_opt->recover_atg_cov_max = INT_MAX;
-    asm_opt->hom_global_coverage = -1;
 }

 void destory_opt(hifiasm_opt_t* asm_opt)
@@ -310,11 +299,6 @@ int check_option(hifiasm_opt_t* asm_opt)
        return 0;
    }

-    if(asm_opt->hom_global_coverage < 0 && asm_opt->hom_global_coverage != -1)
-    {
-        fprintf(stderr, "[ERROR] purge duplication coverage threshold should be >= 0 (--purge-cov)\n");
-        return 0;
-    }


    if(asm_opt->fn_bin_yak[0] != NULL && check_file(asm_opt->fn_bin_yak[0], "YAK1") == 0) return 0;
@@ -378,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:eu", 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:", long_options)) >= 0) {
        if (c == 'h')
        {
            Print_H(asm_opt);
@@ -413,8 +397,6 @@ 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;
@@ -423,17 +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->hom_global_coverage = atoi(opt.arg);
-        else if (c == 310)
-        {
-            char* s = NULL;
-            asm_opt->recover_atg_cov_min = strtol(opt.arg, &s, 10);
-			if (*s == ',') asm_opt->recover_atg_cov_max = strtol(s + 1, &s, 10);
-            if(asm_opt->recover_atg_cov_min == -1 || asm_opt->recover_atg_cov_max == -1)
-            {
-                asm_opt->recover_atg_cov_min = asm_opt->recover_atg_cov_max = -1;
-            }
-        } 
+		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);
@@ -452,13 +424,11 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 		}
    }

-
    if (argc == opt.ind)
    {
        Print_H(asm_opt);
        return 0;
    }
-    ///fprintf(stderr, "max_ov_diff_ec: %f, max_ov_diff_final: %f\n", asm_opt->max_ov_diff_ec, asm_opt->max_ov_diff_final);

    get_queries(argc, argv, &opt, asm_opt);

--- a/CommandLines.h
+++ b/CommandLines.h
@@ -3,7 +3,7 @@

 #include <pthread.h>

-#define HA_VERSION "0.9-r289"
+#define HA_VERSION "0.7-dirty-r257"

 #define VERBOSE 0

@@ -15,8 +15,7 @@
 #define HA_F_SKIP_TRIOBIN    0x20
 #define HA_F_PURGE_CONTAIN   0x40
 #define HA_F_PURGE_JOIN      0x80
-#define HA_F_BAN_POST_JOIN   0x100
-#define HA_F_BAN_ASSEMBLY    0x200
+#define HA_F_KEEP_CONTAINED  0x100

 #define HA_MIN_OV_DIFF       0.02 // min sequence divergence in an overlap

@@ -38,7 +37,6 @@ typedef struct {
 	double max_ov_diff_ec;
 	double max_ov_diff_final;
 	int hom_cov;
-    int het_cov;
 	int max_n_chain; // fall-back max number of chains to consider
    int k_mer_min_freq;
    int k_mer_max_freq;
@@ -58,9 +56,6 @@ typedef struct {
    int purge_level_primary;
    int purge_level_trio;
    int purge_overlap_len;
-    int recover_atg_cov_min;
-    int recover_atg_cov_max;
-    int hom_global_coverage;

    float max_hang_rate;
    float min_drop_rate;
@@ -74,7 +69,6 @@ typedef struct {
    long long num_recorrected_bases;
 	long long mem_buf;
    long long coverage;
-    
 } hifiasm_opt_t;

 extern hifiasm_opt_t asm_opt;
--- a/Correct.h
+++ b/Correct.h
@@ -5,7 +5,6 @@
 #include "Levenshtein_distance.h"
 #include "POA.h"
 #include "Process_Read.h"
-#include "Correct.h"

 //#define CORRECT_THRESHOLD 0.70
 #define CORRECT_THRESHOLD 0.60
--- a/Hash_Table.cpp
+++ b/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)
--- a/21
+++ b/21
@@ -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
--- a/Overlaps.cpp
+++ b/Overlaps.cpp
--- a/Overlaps.h
+++ b/Overlaps.h
@@ -27,7 +27,6 @@
 // #define ALTER_LABLE 2
 // #define HAP_LABLE 4

-
 #define Get_qn(RECORD) ((uint32_t)((RECORD).qns>>32))
 #define Get_qs(RECORD) ((uint32_t)((RECORD).qns))
 #define Get_qe(RECORD) ((RECORD).qe)
@@ -104,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
@@ -113,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 {
@@ -129,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;
@@ -506,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)
 {
@@ -1047,8 +1067,8 @@ R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t is_
 uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges);
 void deduplicate(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig, 
 long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t resolve_tangle);
-void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, 
-ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
+void all_to_all_deduplicate(ma_ug_t* ug, uint8_t postive_flag, float drop_rate, 
+ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
 void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
 void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g,  ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
 ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges);
--- a/Process_Read.cpp
+++ b/Process_Read.cpp
@@ -98,9 +98,6 @@ void write_All_reads(All_reads* r, char* read_file_name)
 	fwrite(r->name, sizeof(char), r->total_name_length, fp);
 	fwrite(r->name_index, sizeof(uint64_t), r->name_index_size, fp);
 	fwrite(r->trio_flag, sizeof(uint8_t), r->total_reads, fp);
-	fwrite(&(asm_opt.hom_cov), sizeof(asm_opt.hom_cov), 1, fp);
-	fwrite(&(asm_opt.het_cov), sizeof(asm_opt.het_cov), 1, fp);
-
    free(index_name);    
 	fflush(fp);
    fclose(fp);
@@ -176,8 +173,6 @@ int load_All_reads(All_reads* r, char* read_file_name)
 	/****************************may have bugs********************************/
 	r->trio_flag = (uint8_t*)malloc(sizeof(uint8_t)*r->total_reads);
 	f_flag += fread(r->trio_flag, sizeof(uint8_t), r->total_reads, fp);
-	f_flag += fread(&(asm_opt.hom_cov), sizeof(asm_opt.hom_cov), 1, fp);
-    f_flag += fread(&(asm_opt.het_cov), sizeof(asm_opt.het_cov), 1, fp);
 	/****************************may have bugs********************************/

 	r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
@@ -200,35 +195,6 @@ int load_All_reads(All_reads* r, char* read_file_name)
 	return 1;
 }

-
-int destory_read_bin(All_reads* r)
-{
-
-	uint64_t i = 0;
-	for (i = 0; i < r->total_reads; i++)
-	{
-		if (r->N_site[i]) free(r->N_site[i]);
-		if (r->read_sperate[i]) free(r->read_sperate[i]);
-		if (r->cigars[i].record) free(r->cigars[i].record);
-		if (r->cigars[i].lost_base) free(r->cigars[i].lost_base);
-		if (r->second_round_cigar[i].record) free(r->second_round_cigar[i].record);
-		if (r->second_round_cigar[i].lost_base) free(r->second_round_cigar[i].lost_base);
-	}
-
-	free(r->N_site);
-	free(r->read_length);
-	free(r->read_size);
-	free(r->read_sperate);
-	free(r->name);
-	free(r->name_index);
-	free(r->trio_flag);
-	free(r->cigars);
-	free(r->second_round_cigar);
-	return 1;
-}
-
-
-
 void ha_insert_read_len(All_reads *r, int read_len, int name_len)
 {
 	r->total_reads++;
--- a/Process_Read.h
+++ b/Process_Read.h
@@ -153,6 +153,5 @@ void reverse_complement(char* pattern, uint64_t length);
 void write_All_reads(All_reads* r, char* read_file_name);
 int load_All_reads(All_reads* r, char* read_file_name);
 void destory_All_reads(All_reads* r);
-int destory_read_bin(All_reads* r);

 #endif
--- a/Purge_Dups.cpp
+++ b/Purge_Dups.cpp
@@ -6,9 +6,6 @@
 #include "Overlaps.h"
 #include "Correct.h"
 #include "kthread.h"
-#include "kdq.h"
-
-KDQ_INIT(uint64_t)

 #define Cal_Off(OFF) ((long long)((uint32_t)((OFF)>>32)) - (long long)((uint32_t)((OFF))))
 #define Get_xOff(OFF) ((long long)((uint32_t)((OFF)>>32)))
@@ -103,7 +100,6 @@ typedef struct {

    ma_ug_t *ug;
    asg_t *read_g;
-    ma_hit_t_alloc* sources;
    ma_hit_t_alloc* reverse_sources;
    R_to_U* ruIndex;
    ma_sub_t *coverage_cut; 
@@ -113,352 +109,9 @@ typedef struct {
    int min_ovlp;
    float chain_rate;
    hap_overlaps_list* all_ovlp;
-    long long cov_threshold;
 }hap_alignment_struct_pip;


-
-void print_peak_line(int c, int x, int exceed, int64_t cnt)
-{
-	int j;
-	if (c >= 0) fprintf(stderr, "[M::%s] %5d: ", __func__, c);
-	else fprintf(stderr, "[M::%s] %5s: ", __func__, "rest");
-	for (j = 0; j < x; ++j) fputc('*', stderr);
-	if (exceed) fputc('>', stderr);
-	fprintf(stderr, " %lld\n", (long long)cnt);
-}
-
-void print_peak(long long* cov_buf, long long cov_buf_length, long long max_i)
-{
-    long long i;
-    const long long hist_max = 100;
-    // print histogram
-	for (i = 0; i < cov_buf_length; ++i) 
-    {
-		long long x, exceed = 0;
-		x = (int)((double)hist_max * cov_buf[i] / cov_buf[max_i] + .499);
-		if (x > hist_max) exceed = 1, x = hist_max; // may happen if cnt[2] is higher
-		if (i > max_i && x == 0) break;
-		print_peak_line(i, x, exceed, cov_buf[i]);
-	}
-	{
-		long long x, exceed = 0;
-		long long rest = 0;
-		for (; i < cov_buf_length; ++i) rest += cov_buf[i];
-		x = (int)((double)hist_max * rest / cov_buf[max_i] + .499);
-		if (x > hist_max) exceed = 1, x = hist_max;
-		print_peak_line(-1, x, exceed, rest);
-	}
-}
-
-
-void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo_peak_cov, 
-long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* coverage_only)
-{
-    long long i, start, err_i, max_i, max2_i, max3_i, topo_peak_i, max, max2, max3, topo_peak, min;
-
-    i = start = err_i = max_i = max2_i = max3_i = topo_peak_i = -1;
-    max = max2 = max3 = topo_peak = min = -1;
-
-    ///cov_buf[0] is usually very large
-    for (i = 1; i < cov_buf_length; ++i)
-    {
-        if(cov_buf[i] > cov_buf[i-1]) break;
-    }
-    err_i = i - 1;
-    // find the global highest peak
-	max_i = err_i + 1, max = cov_buf[max_i];
-    for (i = max_i; i < cov_buf_length; ++i)
-    {
-        if (cov_buf[i] > max)
-        {
-            max = cov_buf[i];
-            max_i = i;
-        }
-    }
-
-    ///print_peak(cov_buf, cov_buf_length, max_i);
-
-    // look for smaller peak on the low end
-	max2 = -1; max2_i = -1;
-	for (i = max_i - 1; i > err_i; --i) 
-    {
-		///at first, it should be a peak
-		if (cov_buf[i] >= cov_buf[i-1] && cov_buf[i] >= cov_buf[i+1]) 
-        {
-			if (cov_buf[i] > max2)
-            {
-                max2 = cov_buf[i];
-                max2_i = i;
-            } 
-		}
-	}
-
-    ///fprintf(stderr, "***max2: %lld, max2_i: %lld\n", max2, max2_i);
-
-    if (max2_i != -1 && max2_i > err_i && max2_i < max_i) 
-    {
-		for (i = max2_i + 1, min = max; i < max_i; ++i)
-        {
-            if (cov_buf[i] < min) min = cov_buf[i];
-        }
-			
-		///if the second peak is not significant
-		if(max2 < max * 0.05 || min > max2 * 0.95) max2 = max2_i = -1;
-	}
-		
-
-    // look for smaller peak on the high end
-	max3 = -1; max3_i = -1;
-    // we'd better use i < cov_buf_length - 1, since cov_buf[cov_buf_length-1] may have problem
-	for (i = max_i + 1; i < cov_buf_length - 1; ++i) 
-    {
-		//at first, it should be a peak
-		if (cov_buf[i] >= cov_buf[i-1] && cov_buf[i] >= cov_buf[i+1]) 
-        {
-			if (cov_buf[i] > max3)
-            {
-                max3 = cov_buf[i], max3_i = i;
-            } 
-		}
-	}
-
-    ///fprintf(stderr, "***max3: %lld, max3_i: %lld\n", max3, max3_i);
-    
-    //if found a peak
-	if (max3 != -1 && max3_i > max_i) 
-    {
-		for (i = max_i + 1, min = max; i < max3_i; ++i)
-        {
-            if (cov_buf[i] < min) min = cov_buf[i];
-        }
-			
-		if (max3 < max * 0.05 || min > max3 * 0.95 || max3_i > max_i * 3) max3 = max3_i = -1;
-	}
-
-
-    
-
-    (*hom_peak) = (*het_peak) = -1;
-    if (topo_peak_cov && (*topo_peak_cov) < cov_buf_length)
-    {
-        topo_peak_i = (*topo_peak_cov);
-        topo_peak = cov_buf[topo_peak_i];
-        if (topo_peak <= max * 0.05) topo_peak_i = topo_peak = -1;
-    }
-    
-    long long k_mer_het, k_mer_hom, coverage_het, coverage_hom, alter_peak;
-    k_mer_het = k_mer_hom = coverage_het = coverage_hom = alter_peak = -1;
-
-    alter_peak = topo_peak_i;
-    k_mer_het = asm_opt.het_cov;
-    k_mer_hom = asm_opt.hom_cov;
-    if(max3_i > 0)
-    {
-        coverage_het = max_i;
-        coverage_hom = max3_i;
-    }
-    else
-    {
-        coverage_het = max2_i;
-        coverage_hom = max_i;
-    }
-
-    if(k_mer_het != -1)
-    {
-        (*het_peak) = k_mer_het;
-        (*hom_peak) = k_mer_hom;
-        return;
-    }
-    else if(coverage_het != -1)
-    {
-        (*het_peak) = coverage_het;
-        (*hom_peak) = coverage_hom;
-        return;
-    }
-    else if(k_mer_hom > coverage_hom*1.5)
-    {
-        (*het_peak) = coverage_hom;
-        (*hom_peak) = k_mer_hom;
-        return;
-    }
-    else if(alter_peak != -1)
-    {
-        ///if peak is het, coverage peak is more reliable 
-        if(coverage_hom >= alter_peak*0.8 && coverage_hom <= alter_peak*1.2) 
-        {
-            (*het_peak) = coverage_hom;
-            return;
-        }///if peak is homo, k-mer peak is more reliable 
-        else if(k_mer_hom >= alter_peak*0.8*2 && k_mer_hom <= alter_peak*1.2*2)
-        {
-            (*hom_peak) = k_mer_hom;
-            return;
-        }
-    }
-
-    (*k_mer_only) = k_mer_hom;
-    (*coverage_only) = coverage_hom;
-    
-
-    // fprintf(stderr, "max: %lld, max_i: %lld\n", max, max_i);
-    // fprintf(stderr, "max2: %lld, max2_i: %lld\n", max2, max2_i);
-    // fprintf(stderr, "max3: %lld, max3_i: %lld\n", max3, max3_i);
-    // fprintf(stderr, "[M::%s] Heterozygous k-mer peak: %d\n", __func__, asm_opt.het_cov);
-    // fprintf(stderr, "[M::%s] Homozygous k-mer peak: %d\n", __func__, asm_opt.hom_cov);
-    // fprintf(stderr, "[M::%s] Heterozygous coverage peak: %lld\n", __func__, (*het_peak));
-    // fprintf(stderr, "[M::%s] Homozygous coverage peak: %lld\n", __func__, (*hom_peak));
-    // fprintf(stderr, "[M::%s] Alter coverage peak: %lld\n", __func__, topo_peak_i);
-}
-
-
-
-
-long long get_alter_peak(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, uint64_t* position_index, 
-ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length)
-{
-    
-    ma_utg_t* u = NULL;
-    asg_t* nsg = ug->g;
-    uint64_t v, j, k, qn, n_vtx = nsg->n_seq, primary_bases = 0, alter_bases = 0;
-    uint32_t tn, is_Unitig;
-    long long* cov_buf = NULL;
-    ma_hit_t *h;
-    cov_buf = (long long*)calloc(cov_buf_length, sizeof(long long));
-    long long R_bases = 0, C_bases_primary = 0, C_bases_alter = 0, C_bases = 0;
-    memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);
-
-
-    for (v = 0; v < n_vtx; ++v) 
-    {
-        if(nsg->seq[v].del) continue;
-        if(nsg->seq[v].c == ALTER_LABLE) continue;
-        u = &(ug->u.a[v]);
-        if(u->m == 0) continue;
-        for (k = 0; k < u->n; k++)
-        {
-            qn = u->a[k]>>33;
-            position_index[qn] = 0;
-            R_bases = coverage_cut[qn].e - coverage_cut[qn].s;
-            primary_bases += R_bases;
-        }
-    }
-
-    for (qn = 0; qn < read_g->n_seq; qn++)
-    {
-        if(position_index[qn] == 0) continue;
-        if(read_g->seq[qn].del) continue;
-
-        C_bases = C_bases_primary = C_bases_alter = 0;
-        R_bases = coverage_cut[qn].e - coverage_cut[qn].s;
-        alter_bases += R_bases;
-        for (j = 0; j < (uint64_t)(sources[qn].length); j++)
-        {
-            h = &(sources[qn].buffer[j]);
-            if(h->el != 1) continue;
-            tn = Get_tn((*h));
-
-            if(read_g->seq[tn].del == 1)
-            {
-                ///get the id of read that contains it 
-                get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
-                if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue;
-            }
-
-            if(position_index[tn] == 0)
-            {
-                C_bases_primary += Get_qe((*h)) - Get_qs((*h));
-            }
-            else
-            {
-                C_bases_alter += Get_qe((*h)) - Get_qs((*h));
-            }
-        }
-
-        // if(qn == 1893151 || qn == 1929038)
-        // {
-        //     fprintf(stderr, "qn: %lu, C_bases_primary: %lld, C_bases_alter: %lld, C_bases: %lld\n",
-        //     qn, C_bases_primary, C_bases_alter, C_bases);
-        // }
-
-        C_bases = C_bases_primary + C_bases_alter;
-        if(C_bases_alter < C_bases * ALTER_COV_THRES) continue;
-
-        C_bases = C_bases/R_bases;
-        if(C_bases < 0 || C_bases >= cov_buf_length) continue;
-        cov_buf[C_bases]++;
-    }
-    
-    long long max_i = -1, max = -1;
-    for (j = 0; (long long)j < cov_buf_length; ++j)
-    {
-        if (cov_buf[j] > max)
-        {
-            max = cov_buf[j];
-            max_i = j;
-        }
-    }
-
-    ///fprintf(stderr, "alter max_i: %lld, max: %lld\n", max_i, max);
-    ///if(max_i < 5) max_i = max = -1;
-    if(alter_bases < primary_bases * REAL_ALTER_THRES) max_i = max = -1;
-
-    free(cov_buf);
-    memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);
-
-    return max_i;
-}
-
-long long get_read_coverage_thres(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, uint64_t* position_index, 
-ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, uint64_t n_read, long long cov_buf_length,
-long long* k_mer_only, long long* coverage_only)
-{
-    uint64_t i, j;
-    long long* cov_buf = NULL;
-    ma_hit_t *h;
-    cov_buf = (long long*)calloc(cov_buf_length, sizeof(long long));
-    long long R_bases = 0, C_bases = 0;
-    for (i = 0; i < n_read; ++i) 
-    {
-        C_bases = 0;
-        R_bases = coverage_cut[i].e - coverage_cut[i].s;
-        for (j = 0; j < (uint64_t)(sources[i].length); j++)
-        {
-            h = &(sources[i].buffer[j]);
-            if(h->el != 1) continue;
-            C_bases += Get_qe((*h)) - Get_qs((*h));
-        }
-        C_bases = C_bases/R_bases;
-        if(C_bases < 0 || C_bases >= cov_buf_length) continue;
-        cov_buf[C_bases]++;
-    }
-
-    long long alter_peak = -1, hom_peak = -1, het_peak = -1;
-    if(position_index)
-    {
-        alter_peak = get_alter_peak(ug, read_g, ruIndex, position_index, sources, coverage_cut, 
-        cov_buf_length);
-    }
-    
-    get_read_peak(cov_buf, cov_buf_length, alter_peak == -1? NULL: &alter_peak, &hom_peak, &het_peak,
-    k_mer_only, coverage_only);
-
-    free(cov_buf);
-
-    if(hom_peak != -1) return hom_peak*HOM_PEAK_RATE;
-    if(het_peak != -1) return het_peak*HET_PEAK_RATE;
-    return -1;
-}
-
-
-
-
-
-
-
-
-
-
 void init_hap_alignment_struct(hap_alignment_struct* x, uint32_t size)
 {
    x->vote_counting = (uint64_t*)malloc(sizeof(uint64_t)*size);
@@ -490,7 +143,7 @@ void destory_hap_alignment_struct(hap_alignment_struct* x)
 }

 void init_hap_alignment_struct_pip(hap_alignment_struct_pip* x, uint32_t num_threads, uint32_t n_seq,
-ma_ug_t *ug, asg_t *read_g,  ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, ma_sub_t *coverage_cut, 
+ma_ug_t *ug, asg_t *read_g,  ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, ma_sub_t *coverage_cut, 
 uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chain_rate, hap_overlaps_list* all_ovlp)
 {
    uint32_t i;
@@ -503,7 +156,6 @@ uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chai

    x->ug = ug;
    x->read_g = read_g;
-    x->sources = sources;
    x->reverse_sources = reverse_sources;
    x->ruIndex = ruIndex;
    x->coverage_cut = coverage_cut;
@@ -835,27 +487,6 @@ long long* r_yBeg, long long* r_yEnd)
    return XCY;
 }

-
-uint64_t get_pair_hap_coverage(uint64_t* readIDs, uint32_t Len, ma_hit_t_alloc* sources, ma_sub_t* coverage_cut)
-{
-    uint32_t m, n, qn;
-    ma_hit_t *h;
-    uint64_t R_bases = 0, C_bases = 0;
-
-    for (m = 0; m < Len; m++)
-    {
-        qn = readIDs[m]>>33;
-        R_bases += coverage_cut[qn].e - coverage_cut[qn].s;
-        for (n = 0; n < (uint64_t)(sources[qn].length); n++)
-        {
-            h = &(sources[qn].buffer[n]);
-            C_bases += Get_qe((*h)) - Get_qs((*h));
-        }
-    }
-
-    return C_bases/R_bases;
-}
-
 void get_pair_hap_similarity(uint64_t* readIDs, uint32_t Len, uint32_t target_uId, 
 ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match, double* Total)
 {
@@ -1827,10 +1458,11 @@ long long* r_y_pos_beg, long long* r_y_pos_end)


 uint32_t calculate_pair_hap_similarity_advance(hap_candidates* hap_can, 
-uint64_t* position_index, uint32_t xUid, uint32_t yUid, ma_utg_t* xReads, ma_utg_t* yReads,
-ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, ma_sub_t *coverage_cut, 
-float Hap_rate, int max_hang, int min_ovlp, uint64_t cov_threshold, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, 
-long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end)
+uint64_t* position_index, uint32_t xUid, uint32_t yUid, ma_utg_t* xReads, ma_utg_t* yReads, 
+ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, ma_sub_t *coverage_cut, 
+float Hap_rate, int max_hang, int min_ovlp, kvec_asg_arc_t_offset* u_buffer, 
+kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, long long* r_x_pos_beg, 
+long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end)
 {
    uint32_t max_count = 0, min_count = 0, flag;
    uint32_t xLen = xReads->n, xIndex;
@@ -1900,7 +1532,6 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon
    if(max_count > min_count*Hap_rate)
    {
        long long r_x_interval_beg, r_x_interval_end, r_y_interval_beg, r_y_interval_end;
-        uint64_t ploid_coverage = 0;

        ///for containment, don't need to do anything
        get_hap_alignment_boundary(xReads, yReads, flag, xLeftMatch, xLeftTotal, 
@@ -1934,15 +1565,6 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon
        if(hap_can->index_end == YCX && xReads->len > (yReads->len*2)) return NON_PLOID;
        if(hap_can->index_end == (uint32_t)-1) return NON_PLOID;

-        ploid_coverage = 0;
-        ploid_coverage += get_pair_hap_coverage(xReads->a+r_x_interval_beg, r_x_interval_end+1-r_x_interval_beg,
-        sources, coverage_cut);
-        ploid_coverage += get_pair_hap_coverage(yReads->a+r_y_interval_beg, r_y_interval_end+1-r_y_interval_beg, 
-        sources, coverage_cut);
-        ///fprintf(stderr, "ploid_coverage: %lu, cov_threshold: %lu\n", ploid_coverage, cov_threshold);
-
-        if(cov_threshold > 0 && ploid_coverage >= cov_threshold) return NON_PLOID;
-
        return PLOID;
    } 
 	return NON_PLOID;
@@ -2803,7 +2425,6 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid)
    hap_alignment_struct_pip* hap_buf = (hap_alignment_struct_pip*)_data;
    ma_ug_t *ug = hap_buf->ug;
    asg_t *read_g = hap_buf->read_g;
-    ma_hit_t_alloc* sources = hap_buf->sources;
    ma_hit_t_alloc* reverse_sources = hap_buf->reverse_sources;
    R_to_U* ruIndex = hap_buf->ruIndex;
    ma_sub_t *coverage_cut = hap_buf->coverage_cut;
@@ -2823,8 +2444,6 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid)
    kvec_t_i32_warp* prevIndex_vec = &(hap_buf->buf[tid].u_buffer_prevIndex);
    kvec_t_i32_warp* begIndex_vec = &(hap_buf->buf[tid].u_buffer_beg);
    kvec_t_u8_warp* flag_vec = &(hap_buf->buf[tid].u_buffer_flag);
-    uint64_t cov_threshold = hap_buf->cov_threshold;
-    if(hap_buf->cov_threshold < 0) cov_threshold = (uint64_t)-1;

    ma_utg_t *xReads = NULL, *yReads = NULL;
    ma_hit_t_alloc *xR = NULL;
@@ -2972,9 +2591,9 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid)
            if(u_can->a.a[k].weight < Get_match(hap_can)*Hap_rate) continue;

            if(calculate_pair_hap_similarity_advance(&(u_can->a.a[k]), position_index, xUid, yUid, 
-            xReads, yReads, sources, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate, max_hang, 
-            min_ovlp, cov_threshold, u_buffer, score_vc, prevIndex_vec, &r_x_pos_beg, &r_x_pos_end, 
-            &r_y_pos_beg, &r_y_pos_end)!=PLOID)
+            xReads, yReads, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate, max_hang, 
+            min_ovlp, u_buffer, score_vc, prevIndex_vec, &r_x_pos_beg, &r_x_pos_end, &r_y_pos_beg, 
+            &r_y_pos_end)!=PLOID)
            {
                continue;
            }
@@ -3968,211 +3587,9 @@ void print_all_purge_ovlp(ma_ug_t *ug, hap_overlaps_list* all_ovlp)
    }

 }
-
-inline int get_available_cnt(asg_t *g, uint32_t v, uint8_t* del, asg_arc_t* v_s)
-{
-    //v has direction
-    if(del && del[v>>1]) return 0;
-    uint32_t i, kv = 0;
-    asg_arc_t *av = asg_arc_a(g, v);
-    uint32_t nv = asg_arc_n(g, v);
-
-    for (i = 0, kv = 0; i < nv; i++)
-    {
-        if(!av[i].del)
-        {
-            if(del && del[av[i].v>>1]) continue;
-            if(v_s) v_s[kv] = av[i];
-            kv++;
-        }
-    }
-
-    return kv;
-}
-
-long long get_specific_contig_length(asg_t *g, uint8_t *del)
-{
-    asg_cleanup(g);
-    uint32_t v, n_vtx = g->n_seq * 2, q_occ;
-    uint8_t *mark = NULL;
-    ///is a queue
-	//kdq_t(uint64_t) *q;
-    ///each node has two directions
-    //q = kdq_init(uint64_t);
-
-
-    mark = (uint8_t*)calloc(n_vtx, 1);
-    
-    long long totalLen = 0;
-    for (v = 0; v < n_vtx; ++v) 
-    {
-        uint32_t w, x, l, start, end, len;
-        asg_arc_t arc;
-        if (g->seq[v>>1].del || mark[v]) continue;
-        if (get_available_cnt(g, v, del, NULL) == 0 && get_available_cnt(g, (v^1), del, NULL) != 0) continue;
-        if (del[v>>1]) continue;
-
-        mark[v] = 1;
-		//q->count = 0, start = v, end = v^1, len = 0;
-        q_occ =0, start = v, end = v^1, len = 0;
-		// forward
-		w = v;
-
-
-        while (1) 
-        {
-			/**
-			 * w----->x
-			 * w<-----x
-			 * that means the only suffix of w is x, and the only prefix of x is w
-			 **/
-			if (get_available_cnt(g, w, del, NULL) != 1) break;
-            get_available_cnt(g, w, del, &arc);
-			x = arc.v; // w->x
-			if (get_available_cnt(g, x^1, del, NULL) != 1) break;
-
-			/**
-			 * another direction of w would be marked as used (since w has been used)
-			**/
-			mark[x] = mark[w^1] = 1;
-			///l is the edge length, instead of overlap length
-            ///note: edge length is different with overlap length
-			///l = asg_arc_len(arc_first(g, w));
-            get_available_cnt(g, w, del, &arc);
-            l = ((uint32_t)((arc).ul));
-			//kdq_push(uint64_t, q, (uint64_t)w<<32 | l);
-            q_occ++;
-			end = x^1, len += l;
-			w = x;
-			if (x == v) break;
-		}
-
-
-        //if (start != (end^1) || kdq_size(q) == 0) { // linear unitig
-        if (start != (end^1) || q_occ == 0) { // linear unitig
-			///length of seq, instead of edge
-			l = g->seq[end>>1].len;
-			//kdq_push(uint64_t, q, (uint64_t)(end^1)<<32 | l);
-            q_occ++;
-			len += l;
-		} else { // circular unitig
-			start = end = UINT32_MAX;
-			goto add_unitig; // then it is not necessary to do the backward
-		}
-
-        // backward
-		x = v;
-		while (1) { // similar to forward but not the same
-			if (get_available_cnt(g, x^1, del, NULL) != 1) break;
-            get_available_cnt(g, x^1, del, &arc);
-            w = arc.v ^ 1;
-			if (get_available_cnt(g, w, del, NULL) != 1) break;
-			mark[x] = mark[w^1] = 1;
-			///l = asg_arc_len(arc_first(g, w));
-            get_available_cnt(g, w, del, &arc);
-            l = ((uint32_t)((arc).ul));
-			///w is the seq id + direction, l is the length of edge
-			///push element to the front of a queue
-			//kdq_unshift(uint64_t, q, (uint64_t)w<<32 | l);
-            q_occ++;
-
-			start = w, len += l;
-			x = w;
-		}
-
-
-        add_unitig:
-        if (start != UINT32_MAX) mark[start] = mark[end] = 1;
-        totalLen += len;
-    }
-    //kdq_destroy(uint64_t, q);
-    return totalLen;
-}
-
-
-void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen)
-{
-    uint8_t *del = (uint8_t *)malloc(sizeof(uint8_t)*g->n_seq);
-    uint32_t v, k;
-    ma_utg_t* u = NULL;
-    memset(del, 1, g->n_seq);
-    (*primaryLen) = (*alterLen) = 0;
-
-    for (v = 0; v < ug->g->n_seq; ++v) 
-    {
-        if(ug->g->seq[v].del) continue;
-        if(ug->g->seq[v].c == ALTER_LABLE) continue;
-        u = &(ug->u.a[v]);
-        if(u->m == 0) continue;
-        for (k = 0; k < u->n; k++)
-        {
-            del[u->a[k]>>33] = 0;
-        }
-    }
-    (*primaryLen) = get_specific_contig_length(g, del);
-
-
-    for (v = 0; v < g->n_seq; ++v) 
-    {
-        del[v] = 1 - del[v];
-    }
-
-    (*alterLen) = get_specific_contig_length(g, del);
-
-    free(del);
-}
-
-
-int if_ploid_sample(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, 
-ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, 
-hap_alignment_struct_pip* hap_buf, hap_overlaps_list* all_ovlp, hap_overlaps_list* back_all_ovlp, 
-uint32_t minLen, double purge_threshold)
-{
-    asg_t* nsg = ug->g;
-    uint64_t v, k, total_bases = 0, alter_bases = 0, primary_bases = 0, purge_bases = 0;
-    kt_for(asm_opt.thread_num, hap_alignment_advance_worker, hap_buf, nsg->n_seq);
-    filter_hap_overlaps_by_length(all_ovlp, minLen);
-    normalize_hap_overlaps_advance(all_ovlp, back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
-
-    get_contig_length(ug, read_g, &primary_bases, &alter_bases);
-    total_bases = primary_bases + alter_bases;
-    // fprintf(stderr, "primary_bases: %lu\n", primary_bases);
-    // fprintf(stderr, "alter_bases: %lu\n", alter_bases);
-    // fprintf(stderr, "total_bases: %lu\n", total_bases);
-
-
-    for (v = 0; v < all_ovlp->num; v++)
-    {
-        for (k = 0; k < all_ovlp->x[v].a.n; k++)
-        {
-            purge_bases += all_ovlp->x[v].a.a[k].x_end_pos - all_ovlp->x[v].a.a[k].x_beg_pos;
-        }
-    }
-    purge_bases = purge_bases/2;
-    ///fprintf(stderr, "purge_bases: %lu\n", purge_bases);
-    alter_bases = alter_bases + purge_bases;
-    ///fprintf(stderr, "new alter_bases: %lu\n", alter_bases);
-
-
-    for (v = 0; v < all_ovlp->num; v++)
-    {
-        all_ovlp->x[v].a.n = 0;
-    }
-
-    for (v = 0; v < back_all_ovlp->num; v++)
-    {
-        back_all_ovlp->x[v].a.n = 0;
-    }
-
-    if(alter_bases > total_bases * purge_threshold) return 1;
-    return 0;
-}
-
-
-void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, 
-ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, 
-uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, 
-uint32_t just_contain, uint32_t just_coverage)
+void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* reverse_sources, 
+R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t purege_minLen, int max_hang, 
+int min_ovlp, long long bubble_dist, float drop_ratio, uint32_t just_contain)
 {
    asg_t *purge_g = NULL;
    purge_g = asg_init();
@@ -4203,19 +3620,7 @@ uint32_t just_contain, uint32_t just_coverage)
    asg_arc_t* p = NULL;
    int r;
    hap_alignment_struct_pip hap_buf;
-    long long k_mer_only, coverage_only;

-    if(asm_opt.hom_global_coverage != -1)
-    {
-        hap_buf.cov_threshold = asm_opt.hom_global_coverage;
-    }
-    else
-    {
-        hap_buf.cov_threshold = get_read_coverage_thres(ug, read_g, ruIndex, position_index, 
-        sources, coverage_cut, read_g->n_seq, COV_COUNT, &k_mer_only, &coverage_only);
-    }
-    
-    
    for (v = 0; v < nsg->n_seq; v++)
    {
        uId = v;
@@ -4244,27 +3649,9 @@ uint32_t just_contain, uint32_t just_coverage)


    init_hap_alignment_struct_pip(&hap_buf, asm_opt.thread_num, nsg->n_seq, ug, read_g,  
-    sources, reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp, 
+    reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp, 
    0.05, &all_ovlp);

-    if(hap_buf.cov_threshold < 0)
-    {
-        if(if_ploid_sample(ug, read_g, ruIndex, sources, reverse_sources, coverage_cut, 
-        &hap_buf, &all_ovlp, &back_all_ovlp, purege_minLen, 0.333))
-        {
-            ///if peak is het, coverage peak is more reliable 
-            hap_buf.cov_threshold = coverage_only * HET_PEAK_RATE;
-        }
-        else
-        {
-            ///if peak is homo, k-mer peak is more reliable 
-            hap_buf.cov_threshold = k_mer_only * HOM_PEAK_RATE;
-        }
-    }
-    if(asm_opt.hom_global_coverage == -1) asm_opt.hom_global_coverage = hap_buf.cov_threshold;
-    fprintf(stderr, "[M::%s] purge duplication coverage threshold: %lld\n", __func__, hap_buf.cov_threshold);
-    if(just_coverage) goto end_coverage;
-
    ///kt_for(asm_opt.thread_num, hap_alignment_worker, &hap_buf, nsg->n_seq);
    kt_for(asm_opt.thread_num, hap_alignment_advance_worker, &hap_buf, nsg->n_seq);

@@ -4372,8 +3759,6 @@ uint32_t just_contain, uint32_t just_coverage)
        }
    }

-    end_coverage:
-
    uint32_t is_Unitig;
    for (v = 0; v < ruIndex->len; v++)
    {
--- a/Purge_Dups.h
+++ b/Purge_Dups.h
@@ -6,19 +6,13 @@
 #include "kdq.h"
 #include "Overlaps.h"
 #include "Hash_Table.h"
-#define COV_COUNT 1024
-#define HOM_PEAK_RATE 1.25
-#define HET_PEAK_RATE (HOM_PEAK_RATE*2)
-#define ALTER_COV_THRES 0.9
-#define REAL_ALTER_THRES 0.1

-void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, 
-ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, 
-uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, 
-uint32_t just_contain, uint32_t just_coverage);
+void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* reverse_sources, 
+R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t purege_minLen, int max_hang, 
+int min_ovlp, long long bubble_dist, float drop_ratio, uint32_t just_contain);
 void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge,
 uint32_t is_circle, uint64_t* rLen);
-void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen);
+
 void enable_debug_mode(uint32_t mode);

 #endif
--- a/README.md
+++ b/README.md
@@ -4,79 +4,18 @@
 # Install hifiasm (requiring g++ and zlib)
 git clone https://github.com/chhylp123/hifiasm
 cd hifiasm && make
-
-# Run on test data (use -f0 for small datasets)
-wget https://github.com/chhylp123/hifiasm/releases/download/v0.7/chr11-2M.fa.gz
-./hifiasm -o test -t4 -f0 chr11-2M.fa.gz 2> test.log
-awk '/^S/{print ">"$1;print $2}' test.p_ctg.gfa > test.p_ctg.fa  # get primary contigs in FASTA
-
-# 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
+# Assembly
+./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz
 ```

 ## Introduction

-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.
+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.

-## 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:
+For non-trio assembly, the input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and its
+outputs consist of: 

 1. Haplotype-resolved raw [unitig][unitig] graph in [GFA][gfa] format
   (*prefix*.r\_utg.gfa). This graph keeps all haplotype information, including
@@ -89,8 +28,7 @@ For non-trio assembly, hifiasm generates the following files:
 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, hifiasm generates the following files:
-
+For trio assembly, the input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and the paternal/maternal trio indexes generated by `yak count` (see https://github.com/lh3/yak). The outputs consist of:
 1. Haplotype-resolved raw [unitig][unitig] graph in [GFA][gfa] format
   (*prefix*.r\_utg.gfa). This graph keeps all haplotype information. 

@@ -100,73 +38,77 @@ For trio assembly, hifiasm generates the following files:
 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

-The following table shows the statistics of several hifiasm primary assemblies:
+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.

-|<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>|
+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>[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>[\[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>|

-[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
+<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>

-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:
+## Usage

-|<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>|
+For Hifi reads assembly, a typical command line looks like:

-[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/
+```sh
+./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz
 ```

+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 brief description of options. If you have
-further questions, please raise an issue at the [issue
-page](https://github.com/chhylp123/hifiasm/issues).
+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.

-## Limitations
+## Limitations and future works

-1. Purging haplotig duplications may introduce misassemblies.
+1. The running time and memory usage should be further reduced.
+
+2. The N50 should be further improved. 
--- a/Trio.cpp
+++ b/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
--- a/anchor.cpp
+++ b/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
--- a/extract.cpp
+++ b/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));
--- a/hifiasm.1
+++ b/hifiasm.1
@@ -1,4 +1,4 @@
-.TH hifiasm 1 "19 July 2020" "hifiasm-0.9 (r289)" "Bioinformatics tools"
+.TH hifiasm 1 "12 Apr 2020" "hifiasm-0.5.0" "Bioinformatics tools"

 .SH NAME
 .PP
@@ -192,13 +192,6 @@ and do the assembly directly and quickly.
 This might be helpful when users want to get an optimized assembly by multiple rounds of experiments
 with different parameters.

-.TP
-.BI --pri-range \ INT1[,INT2]
-Min and max coverage cutoff of primary contigs.
-Keep contigs with coverage in this range at p_ctg.gfa. 
-Inferred automatically in default.
-If INT2 is not specified, it is set to infinity. 
-Set -1 to disable 

 .SS Trio-partition options

@@ -259,11 +252,6 @@ Similarity threshold for duplicate haplotigs that should be purged [0.75].
 .BI -O \ FLOAT
 Min number of overlapped reads for duplicate haplotigs that should be purged [1].

-.TP
-.BI --purge-cov \ INT
-Coverage upper bound of Purge-dups, which is inferred automatically in default.
-If the coverage of a contig is higher than this bound, don't apply Purge-dups.
-
 .SS Debugging options

 .TP 10
@@ -334,3 +322,5 @@ 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.
+
+
--- a/htab.cpp
+++ b/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)
@@ -796,7 +797,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov)
 	return (void*)flt_tab;
 }

-ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov, int *het_cov)
+ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov)
 {
 	int64_t cnt[YAK_N_COUNTS], tot_cnt;
 	int peak_hom, peak_het, i, extra_flag1, extra_flag2;
@@ -818,7 +819,6 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 	fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]);
 	peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
 	if (hom_cov) *hom_cov = peak_hom;
-	if (het_cov) *het_cov = peak_het;
 	if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
 	if (flt_tab == 0) {
 		int cutoff = (int)(peak_hom * asm_opt->high_factor);
--- a/htab.h
+++ b/htab.h
@@ -30,7 +30,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov);
 int ha_ft_isflt(const void *hh, uint64_t y);
 void ha_ft_destroy(void *h);

-ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov, int *het_cov);
+ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov);
 void ha_pt_destroy(ha_pt_t *h);
 const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);

@@ -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__))
--- a/main.cpp
+++ b/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[])
 {
--- a/overlap2.cpp
+++ b/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;
+}
--- a/sys.cpp
+++ b/sys.cpp
@@ -1,6 +1,6 @@
 #include <sys/resource.h>
 #include <sys/time.h>
-#include "htab.h"
+#include "utils.h"

 int yak_verbose = 3;

--- a/utils.h
+++ b/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
Author	SHA1	Message	Date
Heng Li	36bfb02298	r257: dropped contained unitigs	2020-06-08 15:43:21 -04:00
Heng Li	afa93e4bbb	prepare for further changes	2020-06-06 21:55:41 -04:00