purge_dups for high het

CommandLines.cpp

@@ -89,7 +89,7 @@ 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, "    -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]\n", 
                                      asm_opt->purge_simi_rate);
     fprintf(stderr, "    -O INT      min number of overlapped reads for duplicate haplotigs [%d]\n", 
@@ -360,9 +360,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;
     }
 

CommandLines.h

@@ -3,7 +3,7 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.14-r312"
+#define HA_VERSION "0.14-r313"
 
 #define VERBOSE 0
 

Overlaps.h

@@ -478,7 +478,6 @@ void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig, uint8
 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, kvec_asg_arc_t_warp* edge, int require_equal_nv, int test_tangle);
-int asg_pop_bubble_primary(asg_t *g, int max_dist);
 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,35 @@ 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;
+
+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;
+}hap_cov_t;
+
 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);
+int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov);
 
 
 inline int get_real_length(asg_t *g, uint32_t v, uint32_t* v_s)
@@ -1040,10 +1065,10 @@ R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t bac
 uint32_t is_bubble_check, uint32_t is_primary_check);
 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);
+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, hap_cov_t *cov);
 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;
@@ -1084,11 +1109,6 @@ typedef struct{
 
 void init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num);
 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);
 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, 
@@ -1112,6 +1132,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 just_contain, uint32_t just_coverage, hc_links* link, 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);
+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

hifiasm.1

@@ -286,7 +286,8 @@ 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. 
+In default, [2] for non-trio assembly, [0] for trio assembly.
 For trio assembly, only level 0 and level 1 are allowed.
 
 .TP