update r351

Added citation

CommandLines.cpp

@@ -23,13 +23,15 @@ static ko_longopt_t long_options[] = {
 	{ "ex-iter",       ko_required_argument, 308 },
     { "purge-cov",     ko_required_argument, 309 },
     { "pri-range",     ko_required_argument, 310 },
-    { "high-het",      ko_no_argument, 311 },
+    ///{ "high-het",      ko_no_argument, 311 },
     { "lowQ",          ko_required_argument, 312 },
 	{ "min-hist-cnt",  ko_required_argument, 313 },
     { "h1",      ko_required_argument, 314 },
     { "h2",      ko_required_argument, 315 },
     { "enzyme",      ko_required_argument, 316 },
     { "b-cov",      ko_required_argument, 317 },
+    { "h-cov",      ko_required_argument, 318 },
+    { "m-rate",      ko_required_argument, 319 },
 	{ 0, 0, 0 }
 };
 
@@ -68,7 +70,12 @@ void Print_H(hifiasm_opt_t* asm_opt)
     fprintf(stderr, "    --lowQ      INT\n");
     fprintf(stderr, "                output contig regions with >=INT%% inconsistency in BED format; 0 to disable [%d]\n", asm_opt->bed_inconsist_rate);
     fprintf(stderr, "    --b-cov     INT\n");
-    fprintf(stderr, "                break contigs at breakpoints with coverage drop at <INT-fold coverage [%d]\n", asm_opt->break_cov);
+    fprintf(stderr, "                break contigs at positions with <INT-fold coverage; work with '--m-rate'; 0 to disable [%d]\n", asm_opt->b_low_cov);
+    fprintf(stderr, "    --h-cov     INT\n");
+    fprintf(stderr, "                break contigs at positions with >INT-fold coverage; work with '--m-rate'; -1 to disable [%d]\n", asm_opt->b_high_cov);
+    fprintf(stderr, "    --m-rate    FLOAT\n");
+    fprintf(stderr, "                break contigs at positions with <=FLOAT*coverage exact overlaps;\n");
+    fprintf(stderr, "                only work with '--b-cov' or '--h-cov'[%.2f]\n", asm_opt->m_rate);
     
 //	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");
@@ -82,14 +89,14 @@ void Print_H(hifiasm_opt_t* asm_opt)
 	fprintf(stderr, "    -4 FILE     list of hap2/maternal read names []\n");
 
     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, "    -l INT      purge level. 0: no purging; 1: light; 2/3: aggressive [0 for trio; 2 for unzip]\n");
+    fprintf(stderr, "    -s FLOAT    similarity threshold for duplicate haplotigs [%g for -l1/-l2, %g for -l3]\n", 
+                                     asm_opt->purge_simi_rate_l2, asm_opt->purge_simi_rate_l3);
     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, "    --high-het  enable this mode for high heterozygosity sample [experimental, not stable]\n");
+    ///fprintf(stderr, "    --high-het  enable this mode for high heterozygosity sample [experimental, not stable]\n");
     
     fprintf(stderr, "  Hi-C-partition [experimental, not stable]:\n");
     fprintf(stderr, "    --h1 FILEs   file names of Hi-C R1  [r1_1.fq,r1_2.fq,...]\n");
@@ -142,10 +149,11 @@ void init_opt(hifiasm_opt_t* asm_opt)
     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_simi_rate_hic = 0.85;
+    asm_opt->purge_simi_rate_l2 = 0.75;
+    asm_opt->purge_simi_rate_l3 = 0.55;
+    ///asm_opt->purge_simi_rate_hic = 0.85;
     asm_opt->purge_overlap_len = 1;
-    asm_opt->purge_overlap_len_hic = 50;
+    ///asm_opt->purge_overlap_len_hic = 50;
     asm_opt->recover_atg_cov_min = -1024;
     asm_opt->recover_atg_cov_max = INT_MAX;
     asm_opt->hom_global_coverage = -1;
@@ -153,7 +161,10 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->hic_inconsist_rate = 30;
     ///asm_opt->bub_mer_length = 3;
     asm_opt->bub_mer_length = 1000000;
-    asm_opt->break_cov = 0;
+    asm_opt->b_low_cov = 0;
+    asm_opt->b_high_cov = -1;
+    asm_opt->m_rate = 0.75;
+    asm_opt->hap_occ = 1;
 }
 
 void destory_enzyme(enzyme* f)
@@ -350,9 +361,9 @@ int check_option(hifiasm_opt_t* asm_opt)
         return 0;
     }
 
-    if(asm_opt->purge_level_primary < 0 || asm_opt->purge_level_primary > 2)
+    if(asm_opt->purge_level_primary < 0 || asm_opt->purge_level_primary > 3)
     {
-        fprintf(stderr, "[ERROR] the level of purge-dup should be [0, 2] (-l)\n");
+        fprintf(stderr, "[ERROR] the level of purge-dup should be [0, 3] (-l)\n");
         return 0;
     }
 
@@ -419,6 +430,36 @@ int check_option(hifiasm_opt_t* asm_opt)
         return 0;
     }
 
+    if(asm_opt->b_low_cov < 0)
+    {
+        fprintf(stderr, "[ERROR] must >= 0 (--b-cov)\n");
+        return 0;
+    }
+
+    if(asm_opt->b_high_cov != -1 && asm_opt->b_high_cov < 0)
+    {
+        fprintf(stderr, "[ERROR] must >= 0 (--h-cov)\n");
+        return 0;
+    }
+
+    if(asm_opt->m_rate < 0)
+    {
+        fprintf(stderr, "[ERROR] must >= 0 (--m-rate)\n");
+        return 0;
+    }
+
+    if(asm_opt->b_high_cov != -1 && asm_opt->b_high_cov <= asm_opt->b_low_cov)
+    {
+        fprintf(stderr, "[ERROR] [--h-cov] must >= [--b-cov]\n");
+        return 0;
+    }
+
+    if(asm_opt->purge_simi_thres < 0)
+    {
+        fprintf(stderr, "[ERROR] [-s] must >= 0\n");
+        return 0;
+    }
+
     // fprintf(stderr, "input file num: %d\n", asm_opt->num_reads);
     // fprintf(stderr, "output file: %s\n", asm_opt->output_file_name);
     // fprintf(stderr, "number of threads: %d\n", asm_opt->thread_num);
@@ -436,7 +477,7 @@ int check_option(hifiasm_opt_t* asm_opt)
     // 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_simi_thres: %f\n", asm_opt->purge_simi_thres);
     // fprintf(stderr, "purge_overlap_len: %d\n", asm_opt->purge_overlap_len);
 
     return 1;
@@ -586,18 +627,20 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
                 asm_opt->recover_atg_cov_min = asm_opt->recover_atg_cov_max = -1;
             }
         }
-        else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET;
+        ///else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET;
         else if (c == 312) asm_opt->bed_inconsist_rate = atoi(opt.arg);
         else if (c == 313) asm_opt->min_hist_kmer_cnt = atoi(opt.arg);
         else if (c == 314) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[0]), 0);
         else if (c == 315) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[1]), 0);
         else if (c == 316) get_hic_enzymes(opt.arg, &(asm_opt->hic_enzymes), 1);
-        else if (c == 317) asm_opt->break_cov = atoi(opt.arg);
+        else if (c == 317) asm_opt->b_low_cov = atoi(opt.arg);
+        else if (c == 318) asm_opt->b_high_cov = atoi(opt.arg);
+        else if (c == 319) asm_opt->m_rate = atof(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 == 's') asm_opt->purge_simi_rate_l2 = asm_opt->purge_simi_rate_l3 = atof(opt.arg);
         else if (c == 'O') asm_opt->purge_overlap_len = atoll(opt.arg);
         else if (c == ':') 
         {
@@ -611,6 +654,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
 		}
     }
 
+    if(asm_opt->purge_level_primary > 2) asm_opt->purge_simi_thres = asm_opt->purge_simi_rate_l3;
+    else asm_opt->purge_simi_thres = asm_opt->purge_simi_rate_l2;
+
     
     if (argc == opt.ind)
     {

CommandLines.h

@@ -3,7 +3,7 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.14-r309"
+#define HA_VERSION "0.14.2-r315"
 
 #define VERBOSE 0
 
@@ -49,7 +49,9 @@ typedef struct {
 	double max_ov_diff_final;
 	int hom_cov;
     int het_cov;
-    int break_cov;
+    int b_low_cov;
+    int b_high_cov;
+    double m_rate;
 	int max_n_chain; // fall-back max number of chains to consider
 	int min_hist_kmer_cnt;
     int load_index_from_disk;
@@ -68,7 +70,7 @@ typedef struct {
     int purge_level_primary;
     int purge_level_trio;
     int purge_overlap_len;
-    int purge_overlap_len_hic;
+    ///int purge_overlap_len_hic;
     int recover_atg_cov_min;
     int recover_atg_cov_max;
     int hom_global_coverage;
@@ -78,8 +80,10 @@ typedef struct {
     float max_hang_rate;
     float min_drop_rate;
     float max_drop_rate;
-    float purge_simi_rate;
-    float purge_simi_rate_hic;
+    float purge_simi_rate_l2;
+    float purge_simi_rate_l3;
+    float purge_simi_thres;
+    ///float purge_simi_rate_hic;
 
     long long small_pop_bubble_size;
     long long large_pop_bubble_size;
@@ -88,6 +92,7 @@ typedef struct {
     long long num_recorrected_bases;
 	long long mem_buf;
     long long coverage;
+    int hap_occ;
     
 } hifiasm_opt_t;
 

Makefile

@@ -43,7 +43,7 @@ Assembly.o: kthread.h
 CommandLines.o: CommandLines.h ketopt.h
 Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h
 Correct.o: kdq.h CommandLines.h Levenshtein_distance.h POA.h Assembly.h
-Correct.o: ksw2.h
+Correct.o: ksw2.h ksort.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
 Levenshtein_distance.o: Levenshtein_distance.h

Overlaps.h

@@ -477,8 +477,7 @@ void destory_R_to_U(R_to_U* x);
 void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig, uint8_t* flag);
 void get_R_to_U(R_to_U* x, uint32_t rID, uint32_t* uID, uint32_t* is_Unitig);
 void transfor_R_to_U(R_to_U* x);
-void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, int require_equal_nv, int test_tangle);
-int asg_pop_bubble_primary(asg_t *g, int max_dist);
+void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, kvec_asg_arc_t_warp* edge, int require_equal_nv, int test_tangle);
 long long asg_arc_del_simple_circle_untig(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *g, long long circleLen, int is_drop);
 
 typedef struct {
@@ -493,9 +492,21 @@ typedef struct {
 	uint32_t new_edges_i;
 } Edge_iter;
 
+typedef struct {
+	asg_arc_t x;
+    uint64_t Off;
+    uint64_t weight;
+}asg_arc_t_offset;
+
+typedef struct {
+	kvec_t(asg_arc_t_offset) a;
+	uint64_t i;
+}kvec_asg_arc_t_offset;
+
+
+
 void init_Edge_iter(asg_t* g, uint32_t v, asg_arc_t* new_edges, uint32_t new_edges_n, Edge_iter* x);
 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);
 
 
 inline int get_real_length(asg_t *g, uint32_t v, uint32_t* v_s)
@@ -540,55 +551,6 @@ inline uint32_t check_tip(asg_t *sg, uint32_t begNode, uint32_t* endNode, buf_t*
     }
 }
 
-inline uint32_t get_unitig_back(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode, 
-long long* nodeLen, long long* baseLen, buf_t* b)
-{
-	ma_utg_v* u = NULL; 
-    uint32_t v = begNode, w, k;
-    uint32_t kv;
-	(*nodeLen) = (*baseLen) = 0;
-    (*endNode) = (uint32_t)-1;
-	if(ug!=NULL) u = &(ug->u);
-
-    while (1)
-    {
-        kv = get_real_length(sg, v, NULL);
-        (*endNode) = v;
-		if(u == NULL) 
-		{
-			(*nodeLen)++;
-		}
-		else
-		{
-			(*nodeLen) += EvaluateLen((*u), v>>1);
-		}
-        if(b) kv_push(uint32_t, b->b, v);
-		///means reach the end of a unitig
-		if(kv!=1) (*baseLen) += sg->seq[v>>1].len;
-		if(kv==0) return END_TIPS;
-        if(kv>1) return MUL_OUTPUT;
-        ///kv must be 1 here
-        kv = get_real_length(sg, v, &w);
-		///means reach the end of a unitig
-        if(get_real_length(sg, w^1, NULL)!=1)
-		{
-			(*baseLen) += sg->seq[v>>1].len;
-			return MUL_INPUT;
-		}
-
-		for (k = 0; k < asg_arc_n(sg, v); k++)
-        {
-			if(asg_arc_a(sg, v)[k].del) continue;
-			///here is just one undeleted edge
-			(*baseLen) += asg_arc_len(asg_arc_a(sg, v)[k]);
-			break;
-		}
-
-        v = w;
-        if(v == begNode) return LOOP;
-    }
-}
-
 inline uint32_t get_unitig(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode, 
 long long* nodeLen, long long* baseLen, long long* max_stop_nodeLen, long long* max_stop_baseLen, 
 uint32_t stops_threshold, buf_t* b)
@@ -1019,31 +981,49 @@ typedef struct {
 	uint32_t total;
 } Trio_counter;
 
+typedef struct {
+	uint32_t p; // the optimal parent vertex
+	uint32_t d; // the shortest distance from the initial vertex
+	uint32_t r:31, s:1; // r: the number of remaining incoming arc; s: state
+} binfo_s_t;
+
+typedef struct {
+	///all information for each node
+	binfo_s_t *a;
+	kvec_t(uint32_t) S; // set of vertices without parents, nodes with all incoming edges visited
+	kvec_t(uint32_t) b; // visited vertices
+	kvec_t(uint32_t) e; // visited edges/arcs
+} buf_s_t;
+
+typedef struct{
+    buf_s_t *b;
+    uint32_t n_thres, n_reads;
+    asg_t *g;
+    uint32_t check_cross;
+	uint64_t bub_dist;
+} bub_label_t;
+
 void resolve_tangles(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 trio_flag, 
 float drop_ratio);
-void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
+void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, bub_label_t* b_mask_t);
 void rescue_contained_reads_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
-R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t chainLenThres, uint32_t is_bubble_check, 
-uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* new_rtg_nodes);
+R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t chainLenThres, uint32_t is_bubble_check, 
+uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* new_rtg_nodes, bub_label_t* b_mask_t);
 void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
-R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t is_bubble_check, 
-uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges);
+R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t is_bubble_check, uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t);
 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, float double_check_rate);
 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);
+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, bub_label_t* b_mask_t);
 void get_unitig_trio_flag(ma_utg_t* nsu, uint32_t flag, uint32_t* require, uint32_t* non_require, uint32_t* ambigious);
 void rescue_missing_overlaps_backward(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
-R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t backward_steps, 
-uint32_t is_bubble_check, uint32_t is_primary_check);
+R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t backward_steps, uint32_t is_bubble_check, uint32_t is_primary_check, bub_label_t* b_mask_t);
 uint32_t get_edge_from_source(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, 
 R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t query, uint32_t target, asg_arc_t* t);
-uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, int max_dist, buf_t *b, 
-uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes);
 int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio);
-
+void asg_bub_backtrack_primary(asg_t *g, uint32_t v0, buf_t *b);
 
 typedef struct{
     double weight;
@@ -1071,35 +1051,62 @@ typedef struct{
 typedef struct{
     kvec_t(hc_linkeage) a;
     kvec_t(uint64_t) enzymes;
-	kvec_t(bed_in) bed;
-    uint32_t* u_idx;
-	uint64_t r_num;
 } hc_links;
 
+
+typedef struct{
+	kvec_t(uint32_t) uIDs;
+	kvec_t(uint32_t) iDXs;
+	kvec_t(uint32_t) rescue_hom;
+	uint32_t* u_idx;
+	uint8_t* is_het;
+	uint32_t r_num, u_num;
+	uint32_t chain_num;
+	uint32_t l0_chain, l1_chain;
+	kvec_t(bed_in) bed;
+}trans_chain;
+
+typedef struct {
+    uint32_t n;
+    uint32_t* cov;
+    uint64_t* pos_idx;
+    ma_hit_t_alloc* reverse_sources;
+    ma_sub_t *coverage_cut;
+    R_to_U* ruIndex;
+	asg_t *read_g;
+    int max_hang;
+    int min_ovlp;
+    kvec_asg_arc_t_offset u_buffer; 
+    kvec_t_i32_warp tailIndex;
+    kvec_t_i32_warp prevIndex;
+	///hc_links* link;
+	trans_chain* t_ch;
+}hap_cov_t;
+
 typedef struct{
     ///kvec_t(hc_edge) a;
     size_t n, m; 
     hc_edge *a;
 }hc_edge_warp;
 
-void init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num);
+void init_hc_links(hc_links* link, uint64_t ug_num, trans_chain* t_ch);
 void destory_hc_links(hc_links* link);
-void clean_primary_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources,
-long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
-R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen, 
-uint32_t miniBiGraph, float chimeric_rate, int is_final_clean, int just_bubble_pop, 
-float drop_ratio, hc_links* link);
+uint64_t get_bub_pop_max_dist(asg_t *g, buf_t *b);
+uint64_t get_bub_pop_max_dist_advance(asg_t *g, buf_t *b);
+int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov);
+uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, uint64_t max_dist, buf_t *b, uint32_t positive_flag, 
+uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes, hap_cov_t *cov);
+
 void adjust_utg_by_primary(ma_ug_t **ug, asg_t* read_g, float drop_rate,
 ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, 
-long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
+long long tipsLen, float tip_drop_ratio, long long stops_threshold, 
 R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp,
-kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link);
-void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug, asg_t *read_sg);
+kvec_asg_arc_t_warp* new_rtg_edges, hap_cov_t *i_cov, bub_label_t* b_mask_t);
 ma_ug_t* copy_untig_graph(ma_ug_t *src);
 ma_ug_t* output_trio_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, 
-uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist, 
+uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
 long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, 
-float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench);
+float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench, bub_label_t* b_mask_t);
 asg_t* copy_read_graph(asg_t *src);
 ma_ug_t *ma_ug_gen(asg_t *g);
 void ma_ug_destroy(ma_ug_t *ug);
@@ -1112,6 +1119,9 @@ inline int inter_interval(int a_s, int a_e, int b_s, int b_e, int* i_s, int* i_e
     return 1;
 }
 
+
+
+
 #define JUNK_COV 5
 #define DISCARD_RATE 0.8
 

Purge_Dups.h

@@ -11,14 +11,20 @@
 #define HET_PEAK_RATE (HOM_PEAK_RATE*2)
 #define ALTER_COV_THRES 0.9
 #define REAL_ALTER_THRES 0.1
+#define CHAIN_FILTER_RATE 0.7
 
 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, hc_links* link);
+uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain, 
+uint32_t just_coverage, hap_cov_t *cov);
 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);
+hap_cov_t* init_hap_cov_t(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* sources, R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, 
+ma_sub_t *coverage_cut, int max_hang, int min_ovlp, uint32_t is_collect_trans);
+void destory_hap_cov_t(hap_cov_t **x);
+void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd);
+
 
 #endif

README.md

@@ -192,4 +192,8 @@ page](https://github.com/chhylp123/hifiasm/issues).
 
 ## Limitations
 
-1. Purging haplotig duplications may introduce misassemblies.
+1. Purging haplotig duplications may introduce misassemblies.
+
+## Citation
+
+Cheng, H., Concepcion, G.T., Feng, X., Zhang, H., Li H. Haplotype-resolved de novo assembly using phased assembly graphs with hifiasm. Nat Methods 18, 170–175 (2021). https://doi.org/10.1038/s41592-020-01056-5

hic.h

@@ -11,7 +11,7 @@
 
 hc_edge* get_hc_edge(hc_links* link, uint64_t src, uint64_t dest, uint64_t dir);
 void push_hc_edge(hc_linkeage* x, uint64_t uID, double weight, int dir, uint64_t* d);
-void hic_analysis(ma_ug_t *ug, asg_t* read_g, hc_links* link);
+void hic_analysis(ma_ug_t *ug, asg_t* read_g, hap_cov_t *cov);
 void hic_benchmark(ma_ug_t *ug, asg_t* read_g);
 
 typedef struct {
@@ -62,7 +62,7 @@ void get_bubbles(bubble_type* bub, uint64_t id, uint32_t* beg, uint32_t* sink, u
 int load_hc_links(hc_links* link, const char *fn);
 void write_hc_links(hc_links* link, const char *fn);
 void destory_bubbles(bubble_type* bub);
-void identify_bubbles(ma_ug_t* ug, bubble_type* bub, hc_links* link);
+void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *het_flag);
 void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub);
 uint32_t connect_bub_occ(bubble_type* bub, uint32_t root_id, uint32_t check_het);
 void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, uint32_t check_het);

hifiasm.1

@@ -1,4 +1,4 @@
-.TH hifiasm 1 "19 July 2020" "hifiasm-0.9 (r289)" "Bioinformatics tools"
+.TH hifiasm 1 "20 Mar 2021" "hifiasm-0.14.2 (r315)" "Bioinformatics tools"
 
 .SH NAME
 .PP
@@ -212,6 +212,32 @@ with suffix
 .B lowQ.bed
 [70]. Set 0 to disable. 
 
+
+.TP
+.BI --b-cov \ INT
+Break contigs at potential misassemblies with <INT-fold coverage [0].
+Work with 
+.B --m-rate.
+Set 0 to disable. 
+
+.TP
+.BI --h-cov \ INT
+Break contigs at potential misassemblies with >INT-fold coverage [-1].
+Work with 
+.B --m-rate.
+Set -1 to disable. 
+
+.TP
+.BI --m-rate \ FLOAT
+Break contigs with <=FLOAT*coverage exact overlaps [0.75].
+Only work with 
+.B --b-cov
+and
+.B --h-cov.
+
+
+
+
 .SS Trio-partition options
 
 .TP 10
@@ -260,12 +286,14 @@ times in the other sample.
 .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.
+2 to purge all types of haplotigs, 3 to purge all types of haplotigs in most aggressive way 
+for high heterozygosity sample. 
+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].
+Similarity threshold for duplicate haplotigs that should be purged [0.75 for -l1/-l2, 0.55 for -l3].
 
 .TP
 .BI -O \ FLOAT
@@ -276,11 +304,6 @@ Min number of overlapped reads for duplicate haplotigs that should be purged [1]
 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.
 
-.TP
-.BI --high-het \ INT
-Enable this mode for high heterozygosity sample, which will increase running time. 
-For ordinary samples, no need to enable this mode [experimental, not stable].
-
 
 .SS Debugging options
 
@@ -289,6 +312,17 @@ For ordinary samples, no need to enable this mode [experimental, not stable].
 Write additional files to speed up the debugging of graph cleaning.
 
 
+.SS Hi-C-partition options [experimental, not stable]
+
+.TP
+.BI --h1 \ FILEs 
+File names of input Hi-C R1 [r1_1.fq,r1_2.fq,...]
+
+.TP
+.BI --h2 \ FILEs 
+File names of input Hi-C R2 [r2_1.fq,r2_2.fq,...]
+
+
 .SH OUTPUTS
 
 .PP
@@ -347,6 +381,21 @@ phased maternal/haplotype2 contig graph. This graph keeps the phased
 maternal/haplotype2 assembly.
 .RE
 
+.PP
+With Hi-C partition, hifiasm outputs the assembly graphs like trio partition,
+but with additional prefix
+.B [hic].
+In this mode, hifiasm keeps Hi-C alignment results and Hi-C index in two bin 
+files: 
+.B *hic.lk.bin
+and 
+.B *hic.tlb.bin.
+Rerunning hifiasm with different Hi-C reads needs to delete these bin files.
+.RE
+
+
+
+
 .PP
 For each graph, hifiasm also outputs a simplified version without sequences for
 the ease of visualization. Hifiasm keeps corrected reads and overlaps in three