Release hifiasm-0.8 (r279)

Assembly.cpp

@@ -1159,11 +1159,17 @@ void ha_overlap_final(void)
 
 int ha_assemble(void)
 {
+	extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
 	int r, hom_cov = -1, ovlp_loaded = 0;
 	if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
 		ovlp_loaded = 1;
 		fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
+		if (asm_opt.extract_list) {
+			ha_extract_print_list(&R_INF, asm_opt.extract_iter, asm_opt.extract_list);
+			exit(0);
+		}
 		if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN) && !(asm_opt.flag & HA_F_VERBOSE_GFA)) ha_triobin(&asm_opt);
+        ///if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt);
 		if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();
 		if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
 	}
@@ -1195,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

@@ -10,11 +10,15 @@
 hifiasm_opt_t asm_opt;
 
 static ko_longopt_t long_options[] = {
-	{ "version",      ko_no_argument,       300 },
-	{ "dbg-gfa",      ko_no_argument,       301 },
-	{ "write-paf",    ko_no_argument,       302 },
-	{ "write-ec",     ko_no_argument,       303 },
-	{ "skip-triobin", ko_no_argument,       304 },
+	{ "version",       ko_no_argument, 300 },
+	{ "dbg-gfa",       ko_no_argument, 301 },
+	{ "write-paf",     ko_no_argument, 302 },
+	{ "write-ec",      ko_no_argument, 303 },
+	{ "skip-triobin",  ko_no_argument, 304 },
+	{ "max-od-ec",     ko_no_argument, 305 },
+	{ "max-od-final",  ko_no_argument, 306 },
+	{ "ex-list",       ko_required_argument, 307 },
+	{ "ex-iter",       ko_required_argument, 308 },
 	{ 0, 0, 0 }
 };
 
@@ -46,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");
 
@@ -57,6 +62,15 @@ void Print_H(hifiasm_opt_t* asm_opt)
     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        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");
 }
@@ -74,7 +88,9 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->k_mer_length = 51;
 	asm_opt->mz_win = 51;
 	asm_opt->bf_shift = 37;
-	asm_opt->high_factor = 5.0f;
+	asm_opt->high_factor = 5.0;
+	asm_opt->max_ov_diff_ec = 0.04;
+	asm_opt->max_ov_diff_final = 0.03;
 	asm_opt->hom_cov = 20;
 	asm_opt->max_n_chain = 100;
     asm_opt->k_mer_min_freq = 3;
@@ -96,6 +112,10 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->max_short_tip = 3;
     asm_opt->min_cnt = 2;
     asm_opt->mid_cnt = 5;
+    asm_opt->purge_level_primary = 2;
+    asm_opt->purge_level_trio = 0;
+    asm_opt->purge_simi_rate = 0.75;
+    asm_opt->purge_overlap_len = 1;
 }
 
 void destory_opt(hifiasm_opt_t* asm_opt)
@@ -251,6 +271,15 @@ int check_option(hifiasm_opt_t* asm_opt)
         return 0;
     }
 
+	if (asm_opt->max_ov_diff_ec < asm_opt->max_ov_diff_final) {
+		fprintf(stderr, "[ERROR] max_ov_diff_ec shouldn't be smaller than max_ov_diff_final\n");
+		return 0;
+	}
+
+	if (asm_opt->max_ov_diff_ec < HA_MIN_OV_DIFF) {
+		fprintf(stderr, "[ERROR] max_ov_diff_ec shouldn't be smaller than %g\n", HA_MIN_OV_DIFF);
+		return 0;
+	}
 
     if(asm_opt->max_short_tip < 0)
     {
@@ -258,6 +287,19 @@ int check_option(hifiasm_opt_t* asm_opt)
         return 0;
     }
 
+    if(asm_opt->purge_level_primary < 0 || asm_opt->purge_level_primary > 2)
+    {
+        fprintf(stderr, "[ERROR] the level of purge-dup should be [0, 2] (-l)\n");
+        return 0;
+    }
+
+    if(ha_opt_triobin(asm_opt) && ((asm_opt->purge_level_trio < 0 || asm_opt->purge_level_trio > 1)))
+    {
+        fprintf(stderr, "[ERROR] the level of purge-dup for trio should be [0, 1] (-l)\n");
+        return 0;
+    }
+
+
 
     if(asm_opt->fn_bin_yak[0] != NULL && check_file(asm_opt->fn_bin_yak[0], "YAK1") == 0) return 0;
     if(asm_opt->fn_bin_yak[1] != NULL && check_file(asm_opt->fn_bin_yak[1], "YAK2") == 0) return 0;
@@ -279,6 +321,10 @@ int check_option(hifiasm_opt_t* asm_opt)
     // fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip);
     // fprintf(stderr, "min_cnt: %d\n", asm_opt->min_cnt);
     // fprintf(stderr, "mid_cnt: %d\n", asm_opt->mid_cnt);
+    // fprintf(stderr, "purge_level_primary: %d\n", asm_opt->purge_level_primary);
+    // fprintf(stderr, "purge_level_trio: %d\n", asm_opt->purge_level_trio);
+    // fprintf(stderr, "purge_simi_rate: %f\n", asm_opt->purge_simi_rate);
+    // fprintf(stderr, "purge_overlap_len: %d\n", asm_opt->purge_overlap_len);
 
     return 1;
 }
@@ -316,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:", 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);
@@ -351,23 +397,35 @@ 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;
 		else if (c == 304) asm_opt->flag |= HA_F_SKIP_TRIOBIN;
+		else if (c == 305) asm_opt->max_ov_diff_ec = atof(opt.arg);
+		else if (c == 306) asm_opt->max_ov_diff_final = atof(opt.arg);
+		else if (c == 307) asm_opt->extract_list = opt.arg;
+		else if (c == 308) asm_opt->extract_iter = atoi(opt.arg);
+        else if (c == 'l')
+        {   ///0: disable purge_dup; 1: purge containment; 2: purge overlap
+            asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);
+        }
+        else if (c == 's') asm_opt->purge_simi_rate = atof(opt.arg);
+        else if (c == 'O') asm_opt->purge_overlap_len = atoll(opt.arg);
         else if (c == ':') 
         {
 			fprintf(stderr, "[ERROR] missing option argument in \"%s\"\n", argv[opt.i - 1]);
-			return 0;
+			return 1;
 		} 
         else if (c == '?') 
         {
 			fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[opt.i - 1]);
-			return 0;
+			return 1;
 		}
     }
 
-    if (argc == 1)
+    if (argc == opt.ind)
     {
         Print_H(asm_opt);
         return 0;

CommandLines.h

@@ -3,16 +3,22 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.5.0"
+#define HA_VERSION "0.8-r279"
 
 #define VERBOSE 0
 
-#define HA_F_NO_HPC       0x1
-#define HA_F_NO_KMER_FLT  0x2
-#define HA_F_VERBOSE_GFA  0x4
-#define HA_F_WRITE_EC     0x8
-#define HA_F_WRITE_PAF    0x10
-#define HA_F_SKIP_TRIOBIN 0x20
+#define HA_F_NO_HPC          0x1
+#define HA_F_NO_KMER_FLT     0x2
+#define HA_F_VERBOSE_GFA     0x4
+#define HA_F_WRITE_EC        0x8
+#define HA_F_WRITE_PAF       0x10
+#define HA_F_SKIP_TRIOBIN    0x20
+#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_MIN_OV_DIFF       0.02 // min sequence divergence in an overlap
 
 typedef struct {
 	int flag;
@@ -22,11 +28,15 @@ typedef struct {
     char* required_read_name;
 	char *fn_bin_yak[2];
 	char *fn_bin_list[2];
+	char *extract_list;
+	int extract_iter;
     int thread_num;
     int k_mer_length;
 	int mz_win;
 	int bf_shift;
-	float high_factor; // coverage cutoff set to high_factor*hom_cov
+	double high_factor; // coverage cutoff set to high_factor*hom_cov
+	double max_ov_diff_ec;
+	double max_ov_diff_final;
 	int hom_cov;
 	int max_n_chain; // fall-back max number of chains to consider
     int k_mer_min_freq;
@@ -44,10 +54,14 @@ typedef struct {
     int max_short_tip;
     int min_cnt;
     int mid_cnt;
+    int purge_level_primary;
+    int purge_level_trio;
+    int purge_overlap_len;
 
     float max_hang_rate;
     float min_drop_rate;
     float max_drop_rate;
+    float purge_simi_rate;
 
     long long small_pop_bubble_size;
     long long large_pop_bubble_size;

Correct.cpp

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

Correct.h

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

Hash_Table.h

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

Makefile

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

Overlaps.cpp

@@ -1725,10 +1725,9 @@ ma_sub_t* coverage_cut, float shift_rate)
         collect_contain(&(paf[i]), NULL, rLen, &max_left, &max_right, 0.1);
         ///collect_contain(&(paf[i]), &(rev_paf[i]), rLen, &max_left, &max_right, 0.1);
 
-        ////shift_rate should be (FINAL_OVERLAP_ERROR_RATE*2)
+        ////shift_rate should be (asm_opt.max_ov_diff_final*2)
         ///this read is a normal read
-        if(max_left.e > max_right.s && 
-           (max_left.e - max_right.s >= rLen * shift_rate))
+        if (max_left.e > max_right.s && (max_left.e - max_right.s >= rLen * shift_rate))
         {
             continue;
         }
@@ -10435,7 +10434,7 @@ long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, uint32_t negativ
     ///the reason is that each read has two direction (query->target, target->query)
     uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, flag, is_hap, k, to_replace;
     long long ll, base_maxLen, base_maxPositive, base_minNegative, base_minNonPositive, base_best_i, all_covex, curPositive, curNonPositive, curNegative;
-    long long tmp, max_stop_nodeLen, max_stop_baseLen;
+    long long tmp, max_stop_nodeLen, max_stop_baseLen, cur_weight = 0, max_weight = 0;
     uint32_t non_positive_flag = (uint32_t)-1;
     if(positive_flag == FATHER) non_positive_flag = MOTHER;
     if(positive_flag == MOTHER) non_positive_flag = FATHER;
@@ -10501,37 +10500,22 @@ long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, uint32_t negativ
                 }
 
                 to_replace = 0;
-                if(curNegative < base_minNegative)
+                cur_weight = (long long)curPositive - ((long long)curNegative + (long long)curNonPositive);
+                max_weight = (long long)base_maxPositive - ((long long)base_minNegative + (long long)base_minNonPositive);
+                if(cur_weight > max_weight)
                 {
                     to_replace = 1;
                 }
-                else if(curNegative == base_minNegative)
+                else if (cur_weight == max_weight)
                 {
-                    if(curNonPositive < base_minNonPositive)
+                    if(ll > base_maxLen)
                     {
                         to_replace = 1;
                     }
-                    else if(curNonPositive == base_minNonPositive)
-                    {
-                        if(curPositive > base_maxPositive)
-                        {
-                            to_replace = 1;
-                        }
-                        else if(curPositive == base_maxPositive)
-                        {
-                            if(ll > base_maxLen)
-                            {
-                                to_replace = 1;
-                            }
-                        }   
-                    }
                 }
 
 
-
-
-
-                if(to_replace == 1)
+                if(base_best_i == -1 || to_replace == 1)
                 {
                     base_minNegative = curNegative;
                     base_minNonPositive = curNonPositive;
@@ -10610,6 +10594,7 @@ long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, uint32_t negativ
 }
 
 
+
 int asg_arc_cut_long_equal_tips_assembly_complex(asg_t *g, ma_hit_t_alloc* reverse_sources, 
 long long miniedgeLen, uint32_t stops_threshold, R_to_U* ruIndex)
 {
@@ -11378,7 +11363,7 @@ void print_untig(ma_ug_t *g, uint32_t uId, const char* info, uint32_t is_print_r
     asg_arc_t *aw = asg_arc_a(nsg, uId);
     uint32_t nw = asg_arc_n(nsg, uId);
 
-    fprintf(stderr, "\n%s: c = %u\n", info, nsg->seq[uId>>1].c);
+    fprintf(stderr, "\n%s: c = %u, del: %u\n", info, nsg->seq[uId>>1].c, nsg->seq[uId>>1].del);
     fprintf(stderr, "%s(%u): direction = 0...\n", info, uId>>1);
     for (i = 0; i < nw; i++)
     {
@@ -11939,22 +11924,26 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_thre
 }
 
 void renew_longest_tip_by_drop(asg_t *g, ma_ug_t *ug, asg_arc_t *av, uint32_t nv,
-long long* base_maxLen, long long* base_maxLen_i, uint32_t stops_threshold, buf_t* b)
+long long* base_maxLen, long long* base_maxLen_i, uint32_t stops_threshold, buf_t* b, uint32_t trio_flag)
 {
+    if(trio_flag != FATHER && trio_flag != MOTHER) return;
     Trio_counter max, cur;
     memset(&max, 0, sizeof(Trio_counter));
     memset(&cur, 0, sizeof(Trio_counter));
-    long long ll, tmp, max_stop_nodeLen, max_stop_baseLen;
+    long long ll, tmp, max_stop_nodeLen, max_stop_baseLen, max_weight = 0, cur_weight = 0;
     uint32_t convex, i, return_flag, best_tip_i, best_tip_len;
-
+    
     b->b.n = 0;
     return_flag = get_unitig(g, ug, av[(*base_maxLen_i)].v, &convex, &tmp, &ll, 
     &max_stop_nodeLen, &max_stop_baseLen, stops_threshold, b);
-
     if(return_flag!=END_TIPS) return;
+
+
     get_trio_labs(b, ug, &max);
     ///means this unitig might be at current haplotype
-    if(max.drop_occ<(max.total*TRIO_DROP_THRES)) return;
+    ///if(max.drop_occ<(max.total*TRIO_DROP_THRES)) return;
+
+
     best_tip_i = (*base_maxLen_i);
     best_tip_len = (*base_maxLen);
 
@@ -11970,16 +11959,33 @@ long long* base_maxLen, long long* base_maxLen_i, uint32_t stops_threshold, buf_
         if(return_flag!=END_TIPS) return;
         ///this tip should be long enough
         if(ll<(TRIO_DROP_LENGTH_THRES*(*base_maxLen))) continue;
+
         get_trio_labs(b, ug, &cur);
-        ///this unitig is very likly at another haplotype, ignore it
-        if(cur.drop_occ>=(cur.total*TRIO_DROP_THRES)) continue;
-        if(cur.drop_occ<max.drop_occ)
+
+        if(trio_flag == FATHER)
+        {
+            max_weight = (long long)max.father_occ - (long long)max.drop_occ - (long long)max.mother_occ;
+            cur_weight = (long long)cur.father_occ - (long long)cur.drop_occ - (long long)cur.mother_occ;
+            ///this unitig is very likly at another haplotype, ignore it
+            if((cur.drop_occ+cur.mother_occ)>=(cur.total*TRIO_DROP_THRES)) continue;
+        }
+
+        if(trio_flag == MOTHER)
+        {
+            max_weight = (long long)max.mother_occ - (long long)max.drop_occ - (long long)max.father_occ;
+            cur_weight = (long long)cur.mother_occ - (long long)cur.drop_occ - (long long)cur.father_occ;
+            ///this unitig is very likly at another haplotype, ignore it
+            if((cur.drop_occ+cur.father_occ)>=(cur.total*TRIO_DROP_THRES)) continue;
+        }
+
+
+        if(cur_weight > max_weight)
         {
             max = cur;
             best_tip_i = i;
             best_tip_len = ll;
         }
-        else if(cur.drop_occ==max.drop_occ && ll>best_tip_len)
+        else if(cur_weight == max_weight && ll>best_tip_len)
         {
             max = cur;
             best_tip_i = i;
@@ -11992,7 +11998,7 @@ long long* base_maxLen, long long* base_maxLen_i, uint32_t stops_threshold, buf_
 }
 
 int asg_arc_cut_trio_long_equal_tips_assembly(asg_t *g, ma_ug_t *ug, asg_t *read_sg, 
-ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex)
+ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_t trio_flag)
 {
     double startTime = Get_T();
     uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, flag, is_hap, n_tips, return_flag, k;
@@ -12046,7 +12052,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex)
         if(n_arc == n_tips)
         {
             renew_longest_tip_by_drop(g, ug, av, nv, &base_maxLen, 
-            &base_maxLen_i, 1, &b);
+            &base_maxLen_i, 1, &b, trio_flag);
         }
 
         for (i = 0; i < nv; i++)
@@ -12625,12 +12631,11 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate)
     asg_t *g = ug->g;
     uint32_t is_first = 1;
 
-    
     redo:
+    ///if(trio_flag == MOTHER) fprintf(stderr, "(0) c: %u, del: %u, n: %u\n", ug->g->seq[28141].c, ug->g->seq[28141].del, ug->u.a[28141].n);
     asg_pop_bubble_primary_trio(ug, bubble_dist, trio_flag, DROP);
     untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, trio_flag, DROP); 
     magic_trio_phasing(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate);
-    
     ///drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex);
     
     /**********debug**********/
@@ -12654,21 +12659,21 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate)
         ///need consider tangles
         ///asg_pop_bubble_primary(g, bubble_dist);
         asg_pop_bubble_primary_trio(ug, bubble_dist, trio_flag, DROP);
-
-
-        
-
         /**********debug**********/
         if(just_bubble_pop == 0)
         {
             ///need consider tangles
             asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 
             2, tip_drop_ratio);
-            asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex);
+
+            asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag);
+            
             asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex,
             2, tip_drop_ratio, stops_threshold);
+
             asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 
             2, ruIndex, stops_threshold);
+
             detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate,
             ruIndex);
             ///need consider tangles
@@ -12686,13 +12691,10 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate)
         2);
     }
 
-
-    
     resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, trio_flag, drop_ratio);
     drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex);
     all_to_all_deduplicate(ug, trio_flag, trio_drop_rate, reverse_sources, ruIndex);
 
-
     if(is_first)
     {
         is_first = 0;
@@ -12735,7 +12737,7 @@ float drop_ratio)
             ///need consider tangles
             asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 
             2, tip_drop_ratio);
-            asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex);
+            asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1);
             asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex,
             2, tip_drop_ratio, stops_threshold);
             asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 
@@ -12801,10 +12803,10 @@ void set_drop_trio_flag(ma_ug_t *ug)
 
 
 void update_unitig_graph(ma_ug_t* ug, asg_t* read_g, ma_hit_t_alloc* reverse_sources, 
-R_to_U* ruIndex, uint8_t flag, float drop_rate)
+R_to_U* ruIndex, uint8_t is_double_check, uint8_t flag, float drop_rate)
 {
     asg_t* nsg = ug->g;
-    uint32_t v, n_vtx = nsg->n_seq, k, rId, flag_occ, non_flag_occ, n_reduce = 1;
+    uint32_t v, n_vtx = nsg->n_seq, k, rId, flag_occ, non_flag_occ, hap_label_occ, n_reduce = 1;
     ma_utg_t *u;
 
     drop_semi_circle(ug, nsg, read_g, reverse_sources, ruIndex);
@@ -12821,16 +12823,20 @@ R_to_U* ruIndex, uint8_t flag, float drop_rate)
             if((get_real_length(nsg, v<<1, NULL)!=0) 
                 && (get_real_length(nsg, ((v<<1)^1), NULL)!=0)) continue;
 
-            flag_occ = non_flag_occ = 0;
+            flag_occ = non_flag_occ = hap_label_occ = 0;
             for (k = 0; k < u->n; k++)
             {
                 rId = u->a[k]>>33;
+                if(read_g->seq[rId].c == HAP_LABLE) hap_label_occ++;
                 if(R_INF.trio_flag[rId] == AMBIGU) continue;
                 if(R_INF.trio_flag[rId] == DROP) continue;
                 if(R_INF.trio_flag[rId] == flag) flag_occ++;
                 if(R_INF.trio_flag[rId] != flag) non_flag_occ++;
             }
 
+            if(hap_label_occ == u->n) continue;
+            if(is_double_check && non_flag_occ < u->n*DOUBLE_CHECK_THRES) continue;
+
             if(non_flag_occ > ((non_flag_occ+flag_occ)*drop_rate))
             {
                 if(u->m != 0)
@@ -13168,7 +13174,40 @@ void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* re
 
 }
 
+void update_hap_label(ma_ug_t *ug, asg_t* read_g)
+{
+    uint32_t v, n_vtx, k;
+    uint64_t rId;
 
+    if(ug == NULL)
+    {
+        n_vtx = read_g->n_seq;
+        for (v = 0; v < n_vtx; ++v) 
+        {
+            if(read_g->seq[v].del) continue;
+            if(read_g->seq[v].c != HAP_LABLE) continue;
+            read_g->seq[v].c = PRIMARY_LABLE;
+        }
+        return;
+    }
+
+
+    asg_t* nsg = ug->g;
+    n_vtx = nsg->n_seq;
+    ma_utg_t* u = NULL;
+    for (v = 0; v < n_vtx; ++v) 
+    {
+        if(nsg->seq[v].del) continue;
+        if(nsg->seq[v].c != HAP_LABLE) continue;
+        u = &((ug)->u.a[v]);
+        if(u->m == 0) continue;
+        for (k = 0; k < u->n; k++)
+        {
+            rId = u->a[k]>>33;
+            read_g->seq[rId].c = HAP_LABLE;
+        }
+    }
+}
 
 
 void adjust_utg_by_trio(ma_ug_t **ug, asg_t* read_g, uint8_t flag, float drop_rate,
@@ -13180,14 +13219,15 @@ kvec_asg_arc_t_warp* new_rtg_edges)
     asg_t* nsg = (*ug)->g;
     uint32_t v, n_vtx = nsg->n_seq;
     
+    //if(flag == MOTHER) print_untig_by_read(*ug, "m54329U_190617_231905/65340614/ccs", (uint32_t)-1, sources, reverse_sources, "beg1");
+    //if(flag == MOTHER) print_untig_by_read(*ug, "m54329U_190827_173812/67108993/ccs", (uint32_t)-1, sources, reverse_sources, "beg2");
+    
     /**
     kvec_t_u32_warp new_rtg_nodes;
     kv_init(new_rtg_nodes.a);
     **/
-   
-    update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, flag, drop_rate);
+    update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, 0, flag, drop_rate);
     adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex);
-     
     nsg = (*ug)->g;
     n_vtx = nsg->n_seq;
     for (v = 0; v < n_vtx; ++v) 
@@ -13201,31 +13241,44 @@ kvec_asg_arc_t_warp* new_rtg_edges)
     tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, 
     chimeric_rate, 0, 0, drop_ratio, flag, drop_rate);
 
+    ///if(flag == MOTHER) fprintf(stderr, "(o.1) c: %u, del: %u, n: %u\n", (*ug)->g->seq[28141].c, (*ug)->g->seq[28141].del, (*ug)->u.a[28141].n);
+    
     delete_useless_nodes(ug);
 
-    update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, flag, drop_rate);
+    update_hap_label(*ug, read_g);
+
+    update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, 0, flag, drop_rate);
 
     renew_utg(ug, read_g, new_rtg_edges);
 
-    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);
+        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);
+        renew_utg(ug, read_g, new_rtg_edges);
+    }
+    
+    update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, 1, flag, drop_rate);
 
-    update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, flag, drop_rate);
+    update_hap_label(NULL, read_g);
 
     renew_utg(ug, read_g, new_rtg_edges);
 
     delete_useless_nodes(ug);
-    purge_dups(*ug, read_g, coverage_cut, reverse_sources, ruIndex, new_rtg_edges, 0.75, 50, 50, 0.5, max_hang,
-    min_ovlp, bubble_dist, drop_ratio, 1);
-    delete_useless_nodes(ug);
 
+    if(asm_opt.purge_level_trio == 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, 1);
+		delete_useless_nodes(ug);
+	}
 
     set_drop_trio_flag(*ug);
     
@@ -13272,7 +13325,7 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp)
     ///print_untig_by_read(ug, "m64011_190830_220126/117834372/ccs", 865264, sources, reverse_sources, "beg");
     adjust_utg_by_trio(&ug, sg, flag, TRIO_THRES, sources, reverse_sources, coverage_cut, 
     bubble_dist, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, 
-    max_hang, min_ovlp, &new_rtg_edges);
+    max_hang, min_ovlp, &new_rtg_edges);    
     ///debug_utg_graph(ug, sg, 0, 0);
     ///debug_untig_length(ug, tipsLen, gfa_name);
     ///print_untig_by_read(ug, "m64011_190901_095311/125831121/ccs", 2310925, "end");
@@ -13539,6 +13592,7 @@ uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop)
 {   
 	uint32_t i, n_pending = 0, is_first = 1, cur_m, cur_c, cur_np, cur_nc, to_replace, n_tips, tip_end;
 	uint64_t n_pop = 0;
+    long long cur_weight = -1, max_weight = -1;
 	///if this node has been deleted
 	if (g->seq[v0>>1].del || g->seq[v0>>1].c == ALTER_LABLE) return 0; // already deleted
 	///if ((uint32_t)g->idx[v0] < 2) return 0; // no bubbles
@@ -13640,43 +13694,61 @@ uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop)
                     }
                     cur_nc = utg->u.a[(w>>1)].n;
                 }
-                ///select the way with less negative_flag, less non_positive_flag, more positive_flag, more distance
+                ///BUG: select the path with less negative_flag, less non_positive_flag, more positive_flag, more distance
+                ///FIXED: select the path with less (negative_flag+non_positive_flag), more positive_flag, more distance
                 to_replace = 0;
-                if(m + cur_m < t->m)
+
+                /****************************may have bugs********************************/
+                cur_weight = (long long)(c + cur_c) - ((long long)(m + cur_m) + (long long)(np + cur_np));
+                max_weight = (long long)t->c - ((long long)t->m + (long long)t->np);
+                if(cur_weight > max_weight)
                 {
                     to_replace = 1;
                 }
-                else if(m + cur_m == t->m)
+                else if(cur_weight == max_weight)
                 {
-                    if(np + cur_np < t->np)
+                    if(nc + cur_nc > t->nc)
                     {
                         to_replace = 1;
                     }
-                    else if(np + cur_np == t->np)
+                    else if(nc + cur_nc == t->nc)
                     {
-                        if(c + cur_c > t->c)
+                        if(d + l > t->d)
                         {
                             to_replace = 1;
                         }
-                        else if(c + cur_c == t->c)
+                    }
+                }
+                /****************************may have bugs********************************/
+
+                /**
+                if(((m + cur_m) + (np + cur_np)) < (t->m + t->np))
+                {
+                    to_replace = 1;
+                }
+                else if(((m + cur_m) + (np + cur_np)) == (t->m + t->np))
+                {
+                    if(c + cur_c > t->c)
+                    {
+                        to_replace = 1;
+                    }
+                    else if(c + cur_c == t->c)
+                    {
+                        if(nc + cur_nc > t->nc)
                         {
-                            if(nc + cur_nc > t->nc)
+                            to_replace = 1;
+                        }
+                        else if(nc + cur_nc == t->nc)
+                        {
+                            if(d + l > t->d)
                             {
                                 to_replace = 1;
                             }
-                            else if(nc + cur_nc == t->nc)
-                            {
-                                if(d + l > t->d)
-                                {
-                                    to_replace = 1;
-                                }
-                            }
-                            
                         }
                         
                     }
-                    
                 }
+                **/
                 
 
                 if(to_replace)
@@ -21973,7 +22045,7 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov
 kvec_asg_arc_t_warp* new_rtg_edges)
 {
     asg_t* nsg = (*ug)->g;
-    uint32_t v, n_vtx = nsg->n_seq, k, rId;
+    uint32_t v, n_vtx = nsg->n_seq, k, rId, just_contain;
     ma_utg_t* u = NULL;
     
     
@@ -21999,42 +22071,46 @@ kvec_asg_arc_t_warp* new_rtg_edges)
     tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, 
     chimeric_rate, 0, 0, drop_ratio);
 
+    ///delete_useless_nodes(ug);    
     delete_useless_nodes(ug);
+    renew_utg(ug, read_g, new_rtg_edges);
 
-    ///deduplicate_advance(*ug, read_g, coverage_cut, sources, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, 1);
+    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);
+    }
     
-    delete_useless_nodes(ug);
+    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);
 
+        if(asm_opt.purge_level_primary > 0)
+        {
+            just_contain = 0;
+            if(asm_opt.purge_level_primary == 1) just_contain = 1;
 
-    ///enable_debug_mode(1);
-    purge_dups(*ug, read_g, coverage_cut, reverse_sources, ruIndex, new_rtg_edges, 0.75, 50, 50, 0.5, max_hang,
-    min_ovlp, bubble_dist, drop_ratio, 0);
-    delete_useless_nodes(ug);
-    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);
-
-
-    ///debug_purge_dup = 1;
-    ///deduplicate_advance(*ug, read_g, coverage_cut, sources, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, 0);
-    ///delete_useless_nodes(ug);
-    ///enable_debug_mode(0);
-    purge_dups(*ug, read_g, coverage_cut, reverse_sources, ruIndex, new_rtg_edges, 0.75, 50, 50, 0.5, max_hang,
-    min_ovlp, bubble_dist, drop_ratio, 0);
-    delete_useless_nodes(ug);
-    renew_utg(ug, read_g, new_rtg_edges);
+            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);
+        }
+    }
     
 
-
     n_vtx = read_g->n_seq;
     for (v = 0; v < n_vtx; v++)
     {
@@ -26034,7 +26110,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     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");
 
 
@@ -26047,21 +26123,18 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
         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);
     ///ma_hit_sub is just use to init coverage_cut,
     ///it seems we do not need ma_hit_cut & ma_hit_flt
     ma_hit_sub(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut);
-    detect_chimeric_reads(sources, n_read, readLen, coverage_cut, FINAL_OVERLAP_ERROR_RATE*2);
+    detect_chimeric_reads(sources, n_read, readLen, coverage_cut, asm_opt.max_ov_diff_final * 2.0);
     ma_hit_cut(sources, n_read, readLen, mini_overlap_length, &coverage_cut);
     ///print_binned_reads(sources, n_read, coverage_cut);
     ma_hit_flt(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
 
     ///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);
 
@@ -26069,7 +26142,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     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);
@@ -26079,10 +26151,9 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     }
 
     asg_cut_tip(sg, asm_opt.max_short_tip);
-
     
     ///drop_inexact_edegs_at_bubbles(sg, bubble_dist);
-    
+
     if(clean_round > 0)
     {
         double cut_step;
@@ -26098,6 +26169,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
         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;
@@ -26206,7 +26278,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     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);
@@ -26296,13 +26368,17 @@ long long bubble_dist, int read_graph, int write)
         &R_INF, output_file_name);
     }
 
-    try_rescue_overlaps(sources, reverse_sources, n_read, 4); 
+    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, 0);
+        
+        asg_destroy(sg);
+        free(coverage_cut);
+    }
 
-    clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length, 
-    max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 
-    output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 0);
-    
-    asg_destroy(sg);
-    free(coverage_cut);
     destory_R_to_U(&ruIndex);
 }

Overlaps.h

@@ -22,6 +22,7 @@
 #define ALTER_LABLE 1
 #define HAP_LABLE 2
 #define TRIO_THRES 0.9
+#define DOUBLE_CHECK_THRES 0.2
 // #define PRIMARY_LABLE 1
 // #define ALTER_LABLE 2
 // #define HAP_LABLE 4

Purge_Dups.cpp

@@ -2777,6 +2777,25 @@ ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex)
     }
 }
 
+void filter_hap_overlaps_by_length(hap_overlaps_list* all_ovlp, uint32_t minLen)
+{
+    hap_overlaps *x = NULL;
+    uint32_t v, i, m, uId;    
+
+    for (v = 0; v < all_ovlp->num; v++)
+    {
+        uId = v;
+        m = 0;
+        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
+        {
+            x = &(all_ovlp->x[uId].a.a[i]);
+            if(x->x_end_id - x->x_beg_id < minLen) continue;
+            all_ovlp->x[uId].a.a[m] = (*x);
+            m++;
+        }
+        all_ovlp->x[uId].a.n = m;
+    }
+}
 
 void debug_hap_overlaps(hap_overlaps_list* all_ovlp, hap_overlaps_list* back_all_ovlp)
 {
@@ -3569,9 +3588,8 @@ void print_all_purge_ovlp(ma_ug_t *ug, hap_overlaps_list* all_ovlp)
 
 }
 void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* reverse_sources, 
-R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t bi_graph_Len, uint32_t long_hap_overlap, 
-float lable_match_rate, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio,
-uint32_t just_contain)
+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();
@@ -3651,6 +3669,9 @@ uint32_t just_contain)
     //     0.05, &all_ovlp);
     // }
 
+
+    filter_hap_overlaps_by_length(&all_ovlp, purege_minLen);
+
     ///normalize_hap_overlaps(&all_ovlp, &back_all_ovlp);
     normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
     ///debug_hap_overlaps(&all_ovlp, &back_all_ovlp);

Purge_Dups.h

@@ -8,9 +8,8 @@
 #include "Hash_Table.h"
 
 void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* reverse_sources, 
-R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t bi_graph_Len, uint32_t long_hap_overlap, 
-float lable_match_rate, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio,
-uint32_t just_contain);
+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);
 

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.

anchor.cpp

@@ -50,11 +50,18 @@ uint64_t ha_abuf_mem(const ha_abuf_t *ab)
 	return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(ha_mz1_t) + sizeof(seed1_t)) + sizeof(ha_abuf_t);
 }
 
+static int ha_ov_type(const overlap_region *r, uint32_t len)
+{
+	if (r->x_pos_s == 0 && r->x_pos_e == len - 1) return 2; // contained in a longer read
+	else if (r->x_pos_s > 0 && r->x_pos_e < len - 1) return 3; // containing a shorter read
+	else return r->x_pos_s == 0? 0 : 1;
+}
+
 void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain)
 {
 	extern void *ha_flt_tab;
 	extern ha_pt_t *ha_idx;
-	uint32_t i;
+	uint32_t i, rlen;
 	uint64_t k, l;
 	double low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
 	double high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
@@ -64,6 +71,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
     clear_overlap_region_alloc(overlap_list);
 	recover_UC_Read(ucr, &R_INF, rid);
 	ab->mz.n = 0, ab->n_a = 0;
+	rlen = Get_READ_LENGTH(R_INF, rid); // read length
 
 	// get the list of anchors
 	ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab);
@@ -126,30 +134,31 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	calculate_overlap_region_by_chaining(cl, overlap_list, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain);
 
 	#if 0
-	if (overlap_list->length > 2000) {
-		fprintf(stderr, "B\t%ld\t%ld\t%ld\n", (long)rid, (long)overlap_list->length, (long)Get_READ_LENGTH(R_INF, rid));
+	if (overlap_list->length > 0) {
+		fprintf(stderr, "B\t%ld\t%ld\t%d\n", (long)rid, (long)overlap_list->length, rlen);
 		for (int i = 0; i < (int)overlap_list->length; ++i) {
 			overlap_region *r = &overlap_list->list[i];
-			fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand],
-					(int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed);
+			fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand],
+					(int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed, ha_ov_type(r, rlen));
 		}
 	}
 	#endif
 
 	if ((int)overlap_list->length > max_n_chain) {
-		int32_t n[2], s[2];
-		n[0] = n[1] = 0, s[0] = s[1] = 0;
+		int32_t w, n[4], s[4];
+		n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0;
+		ks_introsort_or_ss(overlap_list->length, overlap_list->list);
 		for (i = 0; i < (uint32_t)overlap_list->length; ++i) {
 			const overlap_region *r = &overlap_list->list[i];
-			int dir = r->x_pos_s == 0? 0 : 1;
-			++n[dir];
-			if ((int)n[dir] == max_n_chain) s[dir] = r->shared_seed;
+			w = ha_ov_type(r, rlen);
+			++n[w];
+			if ((int)n[w] == max_n_chain) s[w] = r->shared_seed;
 		}
-		if (s[0] > 0 || s[1] > 0) {
+		if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) {
 			for (i = 0, k = 0; i < (uint32_t)overlap_list->length; ++i) {
 				overlap_region *r = &overlap_list->list[i];
-				int dir = r->x_pos_s == 0? 0 : 1;
-				if (r->shared_seed > s[dir]) {
+				w = ha_ov_type(r, rlen);
+				if (r->shared_seed >= s[w]) {
 					if ((uint32_t)k != i) {
 						overlap_region t;
 						t = overlap_list->list[k];

extract.cpp

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

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
@@ -236,11 +236,27 @@ but occurs <
 times in the other sample.
 
 
+.SS Purge-dups options
+
+.TP 10
+.BI -l \ INT
+Level of purge-dup. 0 to disable purge-dup, 1 to only purge contained haplotigs, 
+2 to purge all types of haplotigs. In default, [2] for non-trio assembly, [0] for trio assembly.
+For trio assembly, only level 0 and level 1 are allowed.
+
+.TP
+.BI -s \ FLOAT
+Similarity threshold for duplicate haplotigs that should be purged [0.75].
+
+.TP
+.BI -O \ FLOAT
+Min number of overlapped reads for duplicate haplotigs that should be purged [1].
+
 .SS Debugging options
 
 .TP 10
 .B --dbg-gfa
-Write additional files to speed up the debugging of graph cleaning
+Write additional files to speed up the debugging of graph cleaning.
 
 
 .SH OUTPUTS
@@ -306,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.
-
-