508 —> futhur graph cleaning

CommandLines.cpp

@@ -53,6 +53,7 @@ static ko_longopt_t long_options[] = {
     { "low-het",     ko_no_argument, 339},
     { "s-base",     ko_required_argument, 340},
     { "bin-only",     ko_no_argument, 341},
+    { "ul-round",     ko_required_argument, 342},
 	{ 0, 0, 0 }
 };
 
@@ -262,7 +263,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->is_topo_trans = 1;
     asm_opt->is_bub_trans = 1;
     asm_opt->bin_only = 0;
-    asm_opt->ul_clean_round = 2;
+    asm_opt->ul_clean_round = 1;
 }
 
 void destory_enzyme(enzyme* f)
@@ -795,6 +796,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
             if(asm_opt->trans_base_rate_sec < 0) asm_opt->is_base_trans = 0;
         } 
         else if (c == 341) asm_opt->bin_only = 1;
+        else if (c == 342) asm_opt->ul_clean_round = atol(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);
         }

CommandLines.h

@@ -4,7 +4,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.18.6-r502"
+#define HA_VERSION "0.18.6-r509"
 
 #define VERBOSE 0
 

Overlaps.cpp

@@ -64,6 +64,7 @@ KSORT_INIT_GENERIC(uint32_t)
 
 void reduce_hamming_error_adv(ma_ug_t *iug, asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, 
 int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub);
+void print_vw_edge(asg_t *sg, uint32_t vid, uint32_t wid, const char *cmd);
 
 typedef struct {
     uint32_t d, tot, ma, p;
@@ -10101,13 +10102,29 @@ uint64_t *n_utg)
     return C_bases/R_bases;
 }
 
+uint32_t cal_circle_ov(const ma_ug_t *ug, uint32_t uid, uint32_t rev, asg_t *sg)
+{
+    uint32_t v, w, vx, wx, k;
+    v = w = (uid<<1)+(!!rev);
+    vx = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32));
+    wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32));
+    v = vx; w = wx;
 
+    if(sg) {
+        asg_arc_t *av; uint32_t nv;
+        av = asg_arc_a(sg, vx); nv = asg_arc_n(sg, vx);
+        for (k = 0; k < nv; k++) {
+            if(av[k].v == wx) return av[k].ol;
+        }
+    }
+    return 0;
+}
 
 void ma_ug_print2(const ma_ug_t *ug, All_reads *RNF, asg_t* read_g, const ma_sub_t *coverage_cut, 
 ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FILE *fp)
 {
     uint8_t* primary_flag = read_g?(uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t)):NULL;
-	uint32_t i, j, l, pc = read_g && coverage_cut && sources && ruIndex?1:0;
+	uint32_t i, j, l, pc = read_g && coverage_cut && sources && ruIndex?1:0, co;
 	char name[32];
 	for (i = 0; i < ug->u.n; ++i) { // the Segment lines in GFA
 		ma_utg_t *p = &ug->u.a[i];
@@ -10154,36 +10171,38 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FIL
         uint32_t nu, u, v;
         for (i = 0; i < ug->u.n; ++i) {
             if(ug->u.a[i].m == 0) continue;
-            if(ug->u.a[i].circ)
+            if(ug->u.a[i].circ) {
-            {
+                co = cal_circle_ov(ug, i, 0, read_g);
-                fprintf(fp, "L\t%s%.6dc\t+\t%s%.6dc\t+\t%dM\tL1:i:%d\n", 
+                fprintf(fp, "L\t%s%.6dc\t+\t%s%.6dc\t+\t%dM\tL1:i:%d\tL2:i:%u\n", 
-                prefix, i+1, prefix, i+1, 0, ug->u.a[i].len);
+                prefix, i+1, prefix, i+1, co, ((ug->u.a[i].len>=co)?(ug->u.a[i].len-co):0), 0);
-                fprintf(fp, "L\t%s%.6dc\t-\t%s%.6dc\t-\t%dM\tL1:i:%d\n", 
+                co = cal_circle_ov(ug, i, 1, read_g);
-                prefix, i+1, prefix, i+1, 0, ug->u.a[i].len);
+                fprintf(fp, "L\t%s%.6dc\t-\t%s%.6dc\t-\t%dM\tL1:i:%d\tL2:i:%u\n", 
-            } 
+                prefix, i+1, prefix, i+1, co, ((ug->u.a[i].len>=co)?(ug->u.a[i].len-co):0), 0);
-            u = i<<1;
+            } else {
-            au = asg_arc_a(ug->g, u);
+                u = i<<1;
-            nu = asg_arc_n(ug->g, u);
+                au = asg_arc_a(ug->g, u);
-            for (j = 0; j < nu; j++)
+                nu = asg_arc_n(ug->g, u);
-            {
+                for (j = 0; j < nu; j++)
-                if(au[j].del) continue;
+                {
-                v = au[j].v;
+                    if(au[j].del) continue;
-                fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\tL2:i:%u\n", 
+                    v = au[j].v;
-                prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1],
+                    fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\tL2:i:%u\n", 
-                prefix,	(v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), 0/**au[j].ou**/);
+                    prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1],
-            }
+                    prefix,	(v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), 0/**au[j].ou**/);
+                }
 
 
-            u = (i<<1) + 1;
+                u = (i<<1) + 1;
-            au = asg_arc_a(ug->g, u);
+                au = asg_arc_a(ug->g, u);
-            nu = asg_arc_n(ug->g, u);
+                nu = asg_arc_n(ug->g, u);
-            for (j = 0; j < nu; j++)
+                for (j = 0; j < nu; j++)
-            {
+                {
-                if(au[j].del) continue;
+                    if(au[j].del) continue;
-                v = au[j].v;
+                    v = au[j].v;
-                fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\tL2:i:%u\n", 
+                    fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\tL2:i:%u\n", 
-                prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1],
+                    prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1],
-                prefix,	(v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), 0/**au[j].ou**/);
+                    prefix,	(v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), 0/**au[j].ou**/);
+                }
             } 
         }
     }
@@ -13494,12 +13513,16 @@ void hic_clean(asg_t* read_g)
 
 void hic_clean_adv(asg_t *sg, ug_opt_t *uopt)
 {
-    uint32_t i, k, m, z, v, w; ma_utg_t *u = NULL; uint32_t *ba, bn, n_vtx, beg, end, n0, n1; 
+    uint32_t i, k, m, z, v, w, mk; ma_utg_t *u = NULL; uint32_t *ba, bn, n_vtx, beg, end, n0, n1; ma_utg_t *mz = NULL;
     ma_ug_t *ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); n_vtx = ug->g->n_seq<<1; double bub_rate = 0.1;
     uint8_t *bf = NULL; bubble_type *bub = gen_bubble_chain(sg, ug, uopt, &bf); 
     uint64_t tLen, vocc, socc, pocc; buf_t b; memset(&b, 0, sizeof(buf_t)); CALLOC(b.a, n_vtx);
     REALLOC(bf, n_vtx); memset(bf, 0, sizeof((*bf))*n_vtx);
     kvec_t(uint64_t) buf; kv_init(buf); n0 = n1 = 0;
+    for (i = 0; i < ug->g->n_seq; ++i) {
+        if(ug->g->seq[i].del) continue;
+        ug->g->seq[i].c = PRIMARY_LABLE;
+    }
 
     for (i = 0; i < bub->b_ug->u.n; i++) {
         u = &(bub->b_ug->u.a[i]);
@@ -13568,6 +13591,20 @@ void hic_clean_adv(asg_t *sg, ug_opt_t *uopt)
                         if((pocc+socc) >= (vocc*bub_rate)) continue;
                         pocc += socc;
                         asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0, NULL);
+                        for (z = 0; z < b.b.n; z++) {
+                            if(b.b.a[z]==v || b.b.a[z]==b.S.a[0]) continue;
+                            // socc += ug->u.a[b.b.a[i]>>1].n;
+                            if(ug->g->seq[b.b.a[z]>>1].del) continue;
+                            if(ug->g->seq[b.b.a[z]>>1].c != ALTER_LABLE) continue;
+                            mz = &(ug->u.a[b.b.a[z]>>1]);
+                            if(mz->m == 0) continue;
+                            for (mk = 0; mk < mz->n; mk++) asg_seq_del(sg, mz->a[mk]>>33);
+                            asg_seq_del(ug->g, b.b.a[z]>>1);
+                            if(ug->u.a[b.b.a[z]>>1].m) {
+                                ug->u.a[b.b.a[z]>>1].m = ug->u.a[b.b.a[z]>>1].n = 0;
+                                free(ug->u.a[b.b.a[z]>>1].a); ug->u.a[b.b.a[z]>>1].a = NULL;
+                            }
+                        }
                         // fprintf(stderr, "+utg%.6dl\tutg%.6dl\n", (int32_t)(v>>1)+1, (int32_t)(b.S.a[0]>>1)+1);
                         n0++;
                     }
@@ -13579,6 +13616,22 @@ void hic_clean_adv(asg_t *sg, ug_opt_t *uopt)
         }
     }
 
+    for (v = 0; v < ug->g->n_seq; ++v) {
+        if(ug->g->seq[v].del) continue;
+        if(ug->g->seq[v].c != ALTER_LABLE) continue;
+        mz = &(ug->u.a[v]);
+        if(mz->m == 0) continue;
+        for (k = 0; k < mz->n; k++) asg_seq_del(sg, mz->a[k]>>33);
+        asg_seq_del(ug->g, v);
+        if(ug->u.a[v].m) {
+            ug->u.a[v].m = ug->u.a[v].n = 0;
+            free(ug->u.a[v].a); ug->u.a[v].a = NULL;
+        }
+    }
+    asg_cleanup(ug->g);
+
+
+
     tLen = get_bub_pop_max_dist_advance(ug->g, &b);
     for (v = buf.n = 0; v < n_vtx; ++v) {
         if(ug->g->seq[v>>1].del) continue;
@@ -13624,14 +13677,43 @@ void hic_clean_adv(asg_t *sg, ug_opt_t *uopt)
                     vocc += ug->u.a[b.b.a[k]>>1].n;
                 }
                 if((socc) >= (vocc*bub_rate)) continue;
+
                 asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0, NULL);
+                for (i = 0; i < b.b.n; i++) {
+                    if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue;
+                    // socc += ug->u.a[b.b.a[i]>>1].n;
+                    if(ug->g->seq[b.b.a[i]>>1].del) continue;
+                    if(ug->g->seq[b.b.a[i]>>1].c != ALTER_LABLE) continue;
+                    mz = &(ug->u.a[b.b.a[i]>>1]);
+                    if(mz->m == 0) continue;
+                    for (mk = 0; mk < mz->n; mk++) asg_seq_del(sg, mz->a[mk]>>33);
+                    asg_seq_del(ug->g, b.b.a[i]>>1);
+                    if(ug->u.a[b.b.a[i]>>1].m) {
+                        ug->u.a[b.b.a[i]>>1].m = ug->u.a[b.b.a[i]>>1].n = 0;
+                        free(ug->u.a[b.b.a[i]>>1].a); ug->u.a[b.b.a[i]>>1].a = NULL;
+                    }
+                }
                 // fprintf(stderr, "-utg%.6dl\tutg%.6dl\n", (int32_t)(v>>1)+1, (int32_t)(b.S.a[0]>>1)+1);
                 n1++;
             }
         }
     }
-    filter_sg_by_ug(sg, ug, uopt);
+    // filter_sg_by_ug(sg, ug, uopt);
+    for (v = 0; v < ug->g->n_seq; ++v) {
+        if(ug->g->seq[v].del) continue;
+        if(ug->g->seq[v].c != ALTER_LABLE) continue;
+        mz = &(ug->u.a[v]);
+        if(mz->m == 0) continue;
+        for (k = 0; k < mz->n; k++) asg_seq_del(sg, mz->a[k]>>33);
+        asg_seq_del(ug->g, v);
+        if(ug->u.a[v].m) {
+            ug->u.a[v].m = ug->u.a[v].n = 0;
+            free(ug->u.a[v].a); ug->u.a[v].a = NULL;
+        }
+    }
+    asg_cleanup(ug->g);
 
+    asg_cleanup(sg);
     free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); 
     ma_ug_destroy(ug); free(bf); kv_destroy(buf);
     destory_bubbles(bub); free(bub);
@@ -15627,7 +15709,6 @@ bub_label_t* b_mask_t)
     ma_ug_t *ug = NULL;
     ug = ma_ug_gen_primary(sg, PRIMARY_LABLE);
 
-
     hap_cov_t *cov = NULL;
     asg_t *copy_sg = copy_read_graph(sg);
     ma_ug_t *copy_ug = copy_untig_graph(ug);
@@ -27603,7 +27684,7 @@ R_to_U* ruIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* new_rtg_edges,
 
 void output_contig_graph_primary_pre(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, 
 ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, uint64_t bubble_dist, long long tipsLen, 
-R_to_U* ruIndex, int max_hang, int min_ovlp)
+R_to_U* ruIndex, int max_hang, int min_ovlp, const ug_opt_t *uopt)
 {
     kvec_asg_arc_t_warp new_rtg_edges;
     kv_init(new_rtg_edges.a);
@@ -27656,6 +27737,10 @@ R_to_U* ruIndex, int max_hang, int min_ovlp)
         fclose(output_file);
     }
 
+    ///for debug
+    // graph_ovlp_binning(ug, sg, uopt);
+    // gen_hpc_re_t(ug);
+
     free(gfa_name);
     ma_ug_destroy(ug);
     kv_destroy(new_rtg_edges.a);
@@ -33584,10 +33669,10 @@ ma_hit_t_alloc* src, uint64_t* readLen, R_to_U* ruIndex, bub_label_t *b_mask_t,
 }
 
 void renew_g(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, long long *n_read, 
-uint64_t **readLen, ma_sub_t **coverage_cut, R_to_U *ruIndex, asg_t **sg, 
+uint64_t **readLen, ma_sub_t **coverage_cut, R_to_U *ruIndex, asg_t **sg, int64_t mini_overlap_length, 
-int64_t mini_overlap_length, int64_t max_hang_length,
+int64_t max_hang_length, ug_opt_t *uopt, int64_t clean_round, double min_ovlp_drop_ratio, 
-ug_opt_t *uopt, int64_t clean_round, double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, 
+double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, 
-int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, const char *bin_file, uint64_t free_uld)
+char *o_file, const char *bin_file, uint64_t free_uld, uint64_t is_bridg, uint64_t deep_clean)
 {
     ul_renew_t nopt; memset(&nopt, 0, sizeof(nopt));
     nopt.src = sources; nopt.r_src = reverse_sources; nopt.ruIndex = ruIndex;
@@ -33595,7 +33680,7 @@ int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, const ch
     nopt.cov = coverage_cut; nopt.b_mask_t = b_mask_t;
     nopt.max_hang = max_hang_length; nopt.mini_ovlp = mini_overlap_length;
     ul_realignment_gfa(uopt, *sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 
-        asm_opt.max_short_tip, asm_opt.max_short_ul_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file, &nopt, bin_file, free_uld);
+        asm_opt.max_short_tip, asm_opt.max_short_ul_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file, &nopt, bin_file, free_uld, is_bridg, deep_clean);
     // ma_ug_t *iug = ul_realignment_gfa(uopt, *sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 
     //                                                                 asm_opt.max_short_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file);
     // asg_t *ng = gen_ng(iug, *sg, uopt, coverage_cut, ruIndex, 256);
@@ -33618,7 +33703,7 @@ bub_label_t *b_mask_t, uint32_t is_trio, int32_t ul_aln_round, char *o_file, con
     int32_t k, strl = strlen(bin_file)+1, kt, cl, sl; char *id = NULL;
     renew_g(src, rev_src, n_read, readLen, cov, ruIndex, sg, mini_overlap_length, max_hang_length,
         uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, asm_opt.max_short_tip, b_mask_t, 
-        is_trio, o_file, bin_file, (ul_aln_round<=1)?1:0);
+        is_trio, o_file, bin_file, (ul_aln_round<=1)?1:0, 1, 1);
     gen_ug_opt_t(uopt, *src, *rev_src, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, *readLen, 
             *cov, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, b_mask_t);
     ug_ext_gfa(uopt, *sg, ug_ext_len);
@@ -33634,7 +33719,7 @@ bub_label_t *b_mask_t, uint32_t is_trio, int32_t ul_aln_round, char *o_file, con
         sprintf(id, "%s%d", bin_file, k);
         renew_g(src, rev_src, n_read, readLen, cov, ruIndex, sg, mini_overlap_length, max_hang_length,
             uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, asm_opt.max_short_tip, b_mask_t, 
-            is_trio, o_file, id, ((k+1)==ul_aln_round)?1:0);
+            is_trio, o_file, id, ((k+1)==ul_aln_round)?1:0, 0, 0);
         gen_ug_opt_t(uopt, *src, *rev_src, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, *readLen, 
             *cov, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, b_mask_t);
         ug_ext_gfa(uopt, *sg, ug_ext_len);
@@ -33902,7 +33987,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
     **/
 
     output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources, 
-        asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length);
+        asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length, &uopt);
     /**
     if (asm_opt.flag & HA_F_VERBOSE_GFA)
     {

Overlaps.h

@@ -67,6 +67,22 @@ typedef struct {
 	size_t n, m;
 } kv_ul_ov_t;
 
+typedef struct {
+	uint32_t tn, rn, el;
+	uint32_t qs, qe, ts, te;
+	uint8_t dir:5, pe:1, full:1, rev:1;
+} emask_t;
+
+typedef struct {
+	emask_t *a; 
+	size_t n, m;
+} kv_emask_t;
+
+typedef struct {
+	kv_emask_t *a; 
+	uint32_t n;
+} idx_emask_t;
+
 typedef struct {
 	///off: start idx in mg128_t * a[];
 	///cnt: how many eles in this chain
@@ -118,8 +134,6 @@ void ma_hit_sort_qns(ma_hit_t *a, long long n);
 int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, 
 char* output_file_name);
 
-
-
 typedef struct {
 	uint32_t s:31, del:1, e;
 	uint8_t c;
@@ -239,6 +253,7 @@ void print_gfa(asg_t *g);
 typedef struct { size_t n, m; uint64_t *a; } asg64_v;
 typedef struct { size_t n, m; uint32_t *a; } asg32_v;
 typedef struct { size_t n, m; ma_utg_t *a;} ma_utg_v;
+typedef struct { asg64_v idx; kv_ul_ov_t srt;} mask_ul_ov_t;
 
 typedef struct {
 	ma_utg_v u;
@@ -294,9 +309,33 @@ typedef struct {
 	hmap_t *mm;
 } hpc_t;
 
+typedef struct {
+	uint32_t s, e;
+	uint8_t k;
+} hpc_ss_t;
+
+typedef struct {
+	uint32_t s, e;
+} hpc_idx_t;
+
+typedef struct {
+	size_t n, m; 
+	hpc_ss_t *a;
+	hpc_idx_t *idx;
+	uint64_t idx_n;
+} hpc_re_t;
+
 #define hpc_len(x, id) ((x).hg->u.a[(id)].len>>1)
 #define hpc_str(x, id, rev) (((x).hg->u.a[(id)].s)+((rev)?((x).hg->u.a[(id)].len>>1):(0)))
 
+typedef struct {
+	uint32_t n;
+	uint8_t *a; 
+} bit_mask_t;
+
+#define set_bit_mask_t(x, i) ((x).a[(i)>>3]|=(((uint8_t)1)<<((i)&((uint32_t)7))))
+#define get_bit_mask_t(x, i) ((x).a[(i)>>3]&(((uint8_t)1)<<((i)&((uint32_t)7))))
+
 typedef struct {
 	ma_ug_t *ug;
 	hpc_t *hpc_g;
@@ -1162,6 +1201,8 @@ void set_reverse_overlap(ma_hit_t* dest, ma_hit_t* source);
 // int* b_low_cov, int* b_high_cov, double m_rate);
 void ma_hit_contained_advance(ma_hit_t_alloc* sources, long long n_read, ma_sub_t *coverage_cut, 
 R_to_U* ruIndex, int max_hang, int min_ovlp);
+void hic_clean_adv(asg_t *sg, ug_opt_t *uopt);
+void update_ug_ou(ma_ug_t *ug, asg_t *sg);
 
 #define JUNK_COV 5
 #define DISCARD_RATE 0.8

Process_Read.cpp

@@ -764,6 +764,11 @@ void destory_all_ul_t(all_ul_t *x) {
 	for (i = 0; i < x->nid.n; i++) free(x->nid.a[i].a);
 	free(x->nid.a);
 	free(x->ridx.idx.a); free(x->ridx.occ.a);
+
+	// if(x->mm) {
+	// 	for (i = 0; i < x->mm->n; i++) free(x->mm->a[i].a);
+	// 	free(x->mm->a); free(x->mm); x->mm = NULL;
+	// }
 }
 
 void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn, uint64_t nn_offset)

Process_Read.h

@@ -198,6 +198,7 @@ typedef struct
     ul_vec_t *a;
     size_t n, m;
     All_reads *hR;
+    // idx_emask_t *mm;
     // uint32_t mm;
 } all_ul_t;
 

anchor.cpp

@@ -1806,3 +1806,23 @@ int64_t ug_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint6
 		return 0;
 	}
 }
+
+int64_t ug_map_lchain_simple(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres,
+								 int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, 
+								 uint32_t gen_off, double mcopy_rate, uint32_t mcopy_khit_cut, uint32_t is_hpc, ha_mzl_t *res, uint64_t res_n, ha_mzl_t *idx, uint64_t idx_n, uint64_t mzl_cutoff, uint64_t chain_cutoff)
+{
+	int64_t max_skip, max_iter, max_dis, quick_check; double chn_pen_gap, chn_pen_skip;
+	set_lchain_dp_op(is_accurate, mz_k, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check);
+	if((!overlap_list) || (!cl)) {
+		ab->mz.n = 0;
+		mz2_ha_sketch(rs, rl, mz_w, mz_k, rid, is_hpc, &ab->mz, NULL, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, NULL);
+		if(res) memcpy(res, ab->mz.a, ab->mz.n * (sizeof((*(ab->mz.a)))));
+		return ab->mz.n;
+	} else {
+		sp->n = 0;
+		minimizers_qgen_input(ab, rid, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, NULL, uref, dbg_ct, sp, high_occ, low_occ, res, res_n, idx, idx_n, mzl_cutoff, chain_cutoff, NULL);
+		// lchain_qgen_mcopy_input(cl, overlap_list, rid, rl, NULL, uref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, bw_thres_sec, quick_check, gen_off, mcopy_rate, mcopy_khit_cut, sp);
+		lchain_qgen_mcopy(cl, overlap_list, rid, rl, NULL, uref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off, mcopy_rate, chain_cutoff, mcopy_khit_cut, sp);
+		return 0;
+	}
+}

gfa_ut.cpp

@@ -380,6 +380,25 @@ typedef struct{
     // ul_path_srt_t psrt;
 }ul_resolve_t;
 
+typedef struct{
+    ug_opt_t *uopt;
+    asg_t *sg;
+    ma_ug_t *ug;
+    buf_t b;
+    uint32_t *idx, *bid;
+    uint32_t idx_n, bid_n;
+    uint8_t is_ou;
+    uint8_t is_trio;
+    int64_t max_ext;
+    uint64_t tlen;
+    double len_rat;
+    double ou_rat;
+    int64_t min_ou;
+}ug_clean_t;
+
+void deep_graph_clean(ug_opt_t *uopt, asg_t *sg, uint8_t is_ou, uint8_t is_trio, int64_t max_ext, 
+double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, double ou_rat, int64_t min_ou, int64_t clean_round, int64_t long_tip);
+
 void init_integer_ml_t(integer_ml_t *x, ul_resolve_t *u, uint64_t n_thread)
 {
     memset(x, 0, sizeof((*x))); 
@@ -1368,6 +1387,44 @@ uint32_t is_topo, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len)
     if (cnt > 0) asg_cleanup(g);
 }
 
+uint8_t get_ug_tip_trio_infor(ma_ug_t *ug, uint32_t begNode)
+{
+    uint32_t v = begNode, w, k, s;
+    uint32_t kv;
+    uint32_t eLen = 0, uLen = 0;
+    uint32_t father_occ = 0, mother_occ = 0, ambigious_occ = 0;
+    ma_utg_t *u;
+
+    while (1) {
+        kv = get_real_length(ug->g, v, NULL);
+        u = &(ug->u.a[v>>1]); eLen += u->n; 
+
+        for (k = 0; k < u->n; k++) {
+            s = u->a[k]>>33;
+            if(R_INF.trio_flag[s]==FATHER) father_occ++;
+            else if(R_INF.trio_flag[s]==MOTHER) mother_occ++;
+            else if((R_INF.trio_flag[s]==AMBIGU) || (R_INF.trio_flag[s]==DROP)) ambigious_occ++;
+        }
+        if(kv!=1) break;
+        ///kv must be 1 here
+        kv = get_real_length(ug->g, v, &w);
+        if(get_real_length(ug->g, w^1, NULL)!=1) break;
+        v = w;
+        if(v == begNode) break;
+    }
+
+    uLen = eLen;///how many nodes
+    eLen = father_occ + mother_occ;///haplotype-sepcific nodes
+    if(eLen == 0) return AMBIGU;
+    if(father_occ >= mother_occ) {
+        if((father_occ > TRIO_THRES*eLen) && (father_occ >= DOUBLE_CHECK_THRES*uLen)) return FATHER;
+    } else {
+        if((mother_occ > TRIO_THRES*eLen) && (mother_occ >= DOUBLE_CHECK_THRES*uLen)) return MOTHER;
+    }
+    return AMBIGU;
+}
+
+
 asg_arc_t *iter_flex_asg(flex_asg_t *fg, flex_asg_e_retrive_t *rr, uint32_t v)
 {
     asg_arc_t *av; uint32_t nv; asg_arc_t *z;
@@ -2255,31 +2312,65 @@ void print_node(asg_t *sg, uint32_t src)
     }
 }
 
-void print_vw_edge(asg_t *sg, uint32_t v, uint32_t w, const char *cmd)
+void print_vw_edge(asg_t *sg, uint32_t vid, uint32_t wid, const char *cmd)
 {
-    asg_arc_t *av; uint32_t nv, i;
+    asg_arc_t *av; uint32_t nv, i, sid, eid;
-    av = asg_arc_a(sg, v); nv = asg_arc_n(sg, v);
+    if(vid >= sg->n_seq || wid >= sg->n_seq) return;
+
+    sid = vid; eid = wid;
+    av = asg_arc_a(sg, (sid<<1)); nv = asg_arc_n(sg, (sid<<1));
     for (i = 0; i < nv; i++) {
-        if(av[i].v == w) {
+        if((av[i].v>>1) == eid) {
             fprintf(stderr, "[%s]\t%.*s(%c)<id:%lu>\t%.*s(%c)<id:%u>\tol:%u\tou:%u\tdel:%u\n", cmd, 
                 (int32_t)Get_NAME_LENGTH(R_INF, (av[i].ul>>33)), Get_NAME(R_INF, (av[i].ul>>33)), "+-"[(av[i].ul>>32)&1], av[i].ul>>33, 
                 (int32_t)Get_NAME_LENGTH(R_INF, (av[i].v>>1)), Get_NAME(R_INF, (av[i].v>>1)), "+-"[av[i].v&1], av[i].v>>1, av[i].ol, av[i].ou, av[i].del);
             break;
         }
     }
-    if(i >= nv) fprintf(stderr, "[%s]\tno edges\n", cmd);
+    av = asg_arc_a(sg, ((sid<<1)+1)); nv = asg_arc_n(sg, ((sid<<1)+1));
+    for (i = 0; i < nv; i++) {
+        if((av[i].v>>1) == eid) {
+            fprintf(stderr, "[%s]\t%.*s(%c)<id:%lu>\t%.*s(%c)<id:%u>\tol:%u\tou:%u\tdel:%u\n", cmd, 
+                (int32_t)Get_NAME_LENGTH(R_INF, (av[i].ul>>33)), Get_NAME(R_INF, (av[i].ul>>33)), "+-"[(av[i].ul>>32)&1], av[i].ul>>33, 
+                (int32_t)Get_NAME_LENGTH(R_INF, (av[i].v>>1)), Get_NAME(R_INF, (av[i].v>>1)), "+-"[av[i].v&1], av[i].v>>1, av[i].ol, av[i].ou, av[i].del);
+            break;
+        }
+    }
+
+    sid = wid; eid = vid;
+    av = asg_arc_a(sg, (sid<<1)); nv = asg_arc_n(sg, (sid<<1));
+    for (i = 0; i < nv; i++) {
+        if((av[i].v>>1) == eid) {
+            fprintf(stderr, "[%s]\t%.*s(%c)<id:%lu>\t%.*s(%c)<id:%u>\tol:%u\tou:%u\tdel:%u\n", cmd, 
+                (int32_t)Get_NAME_LENGTH(R_INF, (av[i].ul>>33)), Get_NAME(R_INF, (av[i].ul>>33)), "+-"[(av[i].ul>>32)&1], av[i].ul>>33, 
+                (int32_t)Get_NAME_LENGTH(R_INF, (av[i].v>>1)), Get_NAME(R_INF, (av[i].v>>1)), "+-"[av[i].v&1], av[i].v>>1, av[i].ol, av[i].ou, av[i].del);
+            break;
+        }
+    }
+    av = asg_arc_a(sg, ((sid<<1)+1)); nv = asg_arc_n(sg, ((sid<<1)+1));
+    for (i = 0; i < nv; i++) {
+        if((av[i].v>>1) == eid) {
+            fprintf(stderr, "[%s]\t%.*s(%c)<id:%lu>\t%.*s(%c)<id:%u>\tol:%u\tou:%u\tdel:%u\n", cmd, 
+                (int32_t)Get_NAME_LENGTH(R_INF, (av[i].ul>>33)), Get_NAME(R_INF, (av[i].ul>>33)), "+-"[(av[i].ul>>32)&1], av[i].ul>>33, 
+                (int32_t)Get_NAME_LENGTH(R_INF, (av[i].v>>1)), Get_NAME(R_INF, (av[i].v>>1)), "+-"[av[i].v&1], av[i].v>>1, av[i].ol, av[i].ou, av[i].del);
+            break;
+        }
+    }
+    // if(i >= nv) fprintf(stderr, "[%s]\tno edges\n", cmd);
 }
 
-void prt_spec_edge(asg_t *rg, ma_hit_t_alloc *src, uint32_t tot_rid, uint32_t vid, uint32_t wid, const char *cmd)
+void prt_spec_edge(asg_t *rg, ma_hit_t_alloc *src, uint32_t tot_rid, uint32_t vid, uint32_t wid, ug_opt_t *uopt, const char *cmd)
 {
     if(vid >= tot_rid) return;
-    uint32_t k, qn, tn; ma_hit_t_alloc *s = &(src[vid]); 
+    uint32_t k, qn, tn; int32_t r; asg_arc_t p; ma_hit_t_alloc *s = &(src[vid]); 
     for (k = 0; k < s->length; k++) {
         qn = Get_qn(s->buffer[k]); tn = Get_tn(s->buffer[k]);
         if(qn == vid && tn == wid) {
-            fprintf(stderr, "[M::%s]%s\tqn::%u\tq::[%u,\t%u)\t%c\ttn::%u\tt::[%u,\t%u)\n", 
+            r = ma_hit2arc(&(s->buffer[k]), rg->seq[qn].len, rg->seq[tn].len, uopt->max_hang, asm_opt.max_hang_rate, uopt->min_ovlp, &p);
-                    __func__, cmd, qn, Get_qs(s->buffer[k]), Get_qe(s->buffer[k]), 
+            fprintf(stderr, "[M::%s::]%s\tqn::%u(%c)\tq::[%u,\t%u)\t%c\ttn::%u(%c)\tt::[%u,\t%u)\n", 
-                    "+-"[s->buffer[k].rev], tn, Get_ts(s->buffer[k]), Get_te(s->buffer[k]));
+                    __func__, cmd, qn, (r>=0)?("+-"[(p.ul>>32)&1]):('*'), Get_qs(s->buffer[k]), Get_qe(s->buffer[k]), 
+                    "+-"[s->buffer[k].rev], 
+                    tn, (r>=0)?("+-"[p.v&1]):('*'), Get_ts(s->buffer[k]), Get_te(s->buffer[k]));
         }
     }
 }
@@ -2348,12 +2439,27 @@ void filter_sg_by_ug(asg_t *rg, ma_ug_t *ug, ug_opt_t *uopt)
             // r = gen_spec_edge(rg, uopt, wx^1, vx^1, &t);
             // assert(r >= 0); p = asg_arc_pushp(rg); *p = t;
         }
+
+        if(u->circ) {
+            v = w = i<<1;
+            vx = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32));
+            wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32));
+            r = gen_spec_edge(rg, uopt, vx, wx, &t);
+            assert(r >= 0); p = asg_arc_pushp(rg); *p = t;
+
+            v = w = (i<<1)^1;
+            vx = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32));
+            wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32));
+            r = gen_spec_edge(rg, uopt, vx, wx, &t);
+            assert(r >= 0); p = asg_arc_pushp(rg); *p = t;
+        }
     }
 
     free(rg->idx);
     rg->idx = 0;
     rg->is_srt = 0;
     asg_cleanup(rg);
+    asg_symm(rg);
 
 
     /*******************************for debug************************************/
@@ -5823,10 +5929,6 @@ void update_raw_integer_seq(poa_g_t *pg, ma_ug_t *ug, uint32_t *cns_seq, uint32_
 
 void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_t is_hom)
 {
-    // if(qid != 3165) return;
-    // if(qid != 17165) return;
-    // if(qid != 24100) return;///circle
-    // if(qid != 27512) return;
     uint64_t k, z, m_het, m_het_occ, ref_occ, b_n, m; uint32_t vk, vz, is_circle = 0; integer_aln_t *p; ul_chain_t sc;
     ul_str_idx_t *str_idx = &(uidx->pstr); ma_ug_t *ug = uidx->l1_ug; 
     uint64_t *hid_a, hid_n; uc_block_t *xi;
@@ -5839,7 +5941,7 @@ void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_
         buf->u.a[k] = ug_occ_w(xi->ts, xi->te, &(ug->u.a[xi->hid]));
         assert(buf->u.a[k] > 0);
         if((!is_hom) && (!IF_HOM((((uint32_t)str->a[k])>>1), (*uidx->bub)))) {
-            m_het++; m_het_occ += buf->u.a[k];
+            m_het++; m_het_occ += buf->u.a[k];///m_het_occ: how many het HiFi reads
         }
         ref_occ += buf->u.a[k];
         // fprintf(stderr, "[M::%s::k->%lu] buf->u.a[k]->%lu, ts->%u, te->%u, pchain->%u\n", __func__, k, buf->u.a[k], xi->ts, xi->te, xi->pchain);
@@ -14344,6 +14446,11 @@ asg_t *renew_ng(usg_t *eg, ma_ug_t *rug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **c
     memset(ext.idx_a.a, -1, sizeof(*(ext.idx_a.a))*ext.idx_a.n);//map
     ug1 = convert_usg_t(eg, rug);
 
+    // fprintf(stderr, "-0-[M::%s]\n", __func__);
+    // print_vw_edge(sg, 681356, 73351, __func__);
+    // prt_spec_edge(sg, R_INF.paf, R_INF.total_reads, 681356, 73351, uopt, "src");
+    // // prt_spec_edge(sg, R_INF.reverse_paf, R_INF.total_reads, 681356, 73351, uopt, "rsrc");
+
     for (i = 0; i < ug1->u.n; ++i) {
         if(ug1->g->seq[i].del) continue;
         u = &(ug1->u.a[i]);        
@@ -14375,6 +14482,7 @@ asg_t *renew_ng(usg_t *eg, ma_ug_t *rug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **c
     ng->r_seq = ng->n_seq;
     assert(ng->n_seq == ext.idx_a.n);
 
+
     for (i = 0, nocc = ng->r_seq, sa.n = 0; i < eg->n; ++i) {
         if(ug1->g->seq[i].del) continue;
         ///there shouldn't any scaffolding within the nodes
@@ -14435,11 +14543,9 @@ asg_t *renew_ng(usg_t *eg, ma_ug_t *rug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **c
     CALLOC(ng->seq_vis, (ng->n_seq<<1));
     ///# scaffolding nodes
     if(sa.n > 0) fill_scaffolds(sa.a, sa.n, rug, uopt->min_ovlp, bin_file);
-
     realloc_rdb_adv(&(R_INF), &UL_INF, cov, ruI, ext.idx_a.a, ext.idx_a.n, scaffold_len, (char *)"scaf", ng, uopt, rug, sa.a);     
     update_paf(R_INF.paf, R_INF.reverse_paf, ext.idx_a.a, ext.idx_a.n, sg->n_seq, ng);
 
-
     for (i = 0, nocc = ng->r_seq; i < eg->n; ++i) {
         if(ug1->g->seq[i].del) continue;
         u = &(ug1->u.a[i]);
@@ -14542,6 +14648,7 @@ asg_t *renew_ng(usg_t *eg, ma_ug_t *rug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **c
 
     asg_cleanup(ng); ng->r_seq = ng->n_seq;
     free(ext.cnt.a); free(ext.idx_a.a); free(sa.a); ma_ug_destroy(ug1);
+
     fprintf(stderr, "[M::%s] nocc::%lu, ng->n_seq::%u, sg->n_seq::%u\n", 
     __func__, nocc, (uint32_t)ng->n_seq, (uint32_t)sg->n_seq);
     assert(nocc == ng->n_seq);
@@ -15350,7 +15457,7 @@ ma_ug_t *gen_hybrid_ug(ul_resolve_t *uidx, usg_t *ng)
 
 void renew_ul2_utg(ul_resolve_t *uidx);
 
-void u2g_threading(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint64_t cov_cutoff, asg64_v *b, asg64_v *ub)
+void u2g_threading(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint64_t cov_cutoff, uint64_t is_bridg, asg64_v *b, asg64_v *ub)
 {
     renew_ul2_utg(uidx);
     ul2ul_idx_t *idx = &(uidx->uovl); ma_ug_t *iug = idx->i_ug; ma_ug_t *raw = uidx->l1_ug;
@@ -15425,7 +15532,7 @@ void u2g_threading(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint64_t cov_cutoff, as
                         if(p->ou < (uint64_t)tm) p->ou = tm;
                         p = get_usg_arc(ng, seq->a[i+1].v^1, seq->a[i].v^1);
                         if(p->ou < (uint64_t)tm) p->ou = tm;
-                    } else {
+                    } else if(is_bridg) {
                         v = seq->a[i].v; w = seq->a[i+1].v; p = NULL;
                         vx = (v&1?((raw->u.a[v>>1].a[0]>>32)^1):(raw->u.a[v>>1].a[raw->u.a[v>>1].n-1]>>32));
                         wx = (w&1?((raw->u.a[w>>1].a[raw->u.a[w>>1].n-1]>>32)^1):(raw->u.a[w>>1].a[0]>>32));
@@ -15480,7 +15587,7 @@ void u2g_threading(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint64_t cov_cutoff, as
     // fprintf(stderr, "-[M::%s::] idx->hybrid_ug->g->n_seq::%u\n", __func__, (uint32_t)idx->hybrid_ug->g->n_seq);
 }
 
-void u2g_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint32_t keep_raw_utg)
+void u2g_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint32_t keep_raw_utg, uint64_t is_bridg)
 {
     ul2ul_idx_t *idx = &(uidx->uovl); asg64_v bu = {0,0,0}, uu = {0,0,0}; int64_t max_tip_hifi0 = ulopt->max_tip_hifi;
     ma_ug_t *iug = idx->i_ug; int64_t i, mm_tip = ulopt->max_tip; uint64_t cnt = 1, topo_level, ss = 0;
@@ -15527,7 +15634,7 @@ void u2g_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint32_t keep_raw_utg)
     }
     ulopt->max_tip_hifi = max_tip_hifi0;
 
-    u2g_threading(uidx, ulopt, 3, &bu, &uu);
+    u2g_threading(uidx, ulopt, 3, is_bridg, &bu, &uu);
 
 
     // fill_u2g(uidx, &bu, &uu);
@@ -15995,7 +16102,7 @@ void renew_ul2_utg(ul_resolve_t *uidx)
 }
 
 
-ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt, uint32_t keep_raw_utg)
+ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt, uint32_t keep_raw_utg, uint64_t is_bridg)
 {
     uint64_t k, m;
     ma_ug_t *ug = uidx->l1_ug; all_ul_t *uls = uidx->idx; 
@@ -16039,7 +16146,7 @@ ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt, uin
     renew_u2g_bg(uidx);
     
     // output_integer_graph(uidx, z->i_ug, asm_opt.output_file_name);
-    u2g_clean(uidx, ulopt, keep_raw_utg);
+    u2g_clean(uidx, ulopt, keep_raw_utg, is_bridg);
     // renew_ul2_utg(uidx);
 
     // output_integer_graph(uidx, z->i_ug, asm_opt.output_file_name, 0);
@@ -16763,7 +16870,7 @@ R_to_U* ruIndex, int max_hang, int min_ovlp)
 
 void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, 
 double max_ovlp_drop_ratio, int64_t max_tip, int64_t max_ul_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file,
-ul_renew_t *ropt, const char *bin_file, uint64_t free_uld)
+ul_renew_t *ropt, const char *bin_file, uint64_t free_uld, uint64_t is_bridg, uint64_t deep_clean)
 {
     uint64_t i, bn = 0, idn = 0; uint32_t *bl = NULL; uint8_t *r_het = NULL; bubble_type *bub = NULL; ulg_opt_t uu;
     for (i = 0; i < sg->n_seq; ++i) {
@@ -16772,10 +16879,25 @@ ul_renew_t *ropt, const char *bin_file, uint64_t free_uld)
     }
     // fprintf(stderr, "0[M::%s]\n", __func__);
     hic_clean(sg);
+    if(deep_clean) {
+        // print_debug_gfa(sg, NULL, uopt->coverage_cut, "bclean", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
+        deep_graph_clean(uopt, sg, 1, is_trio, asm_opt.max_short_tip, asm_opt.min_drop_rate, 
+        MIN(asm_opt.max_drop_rate, 0.7), 0.75, 1, asm_opt.clean_round, 20);
+        // print_debug_gfa(sg, NULL, uopt->coverage_cut, "aclean", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
+    }
+    
     // fprintf(stderr, "1[M::%s]\n", __func__);
     ma_ug_t *init_ug = ul_realignment(uopt, sg, 0, bin_file);
     // fprintf(stderr, "2[M::%s]\n", __func__);
     // exit(1);
+
+    // char* gfa_name = NULL; MALLOC(gfa_name, strlen(o_file)+strlen(bin_file)+50);
+    // sprintf(gfa_name, "%s.%s", o_file, bin_file);
+    // print_debug_gfa(sg, init_ug, uopt->coverage_cut, gfa_name, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
+    // // print_debug_gfa(sg, init_ug, uopt->coverage_cut, gfa_name, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1);
+    // free(gfa_name);
+    
+    
     filter_sg_by_ug(sg, init_ug, uopt);
     // fprintf(stderr, "-0-[M::%s]\tUL_INF.a[25].rlen::%u\n", __func__, UL_INF.a[25].rlen);
     // print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
@@ -16791,12 +16913,15 @@ ul_renew_t *ropt, const char *bin_file, uint64_t free_uld)
     // fprintf(stderr, "5[M::%s]\n", __func__);
     // print_ul_alignment(init_ug, &UL_INF, 47072, "after-2");
     // exit(1);
-    // print_raw_uls_seq(uidx, asm_opt.output_file_name);
+    // if(free_uld) {
-    // print_raw_uls_aln(uidx, asm_opt.output_file_name);
+    //     print_raw_uls_seq(uidx, asm_opt.output_file_name);
+    //     print_raw_uls_aln(uidx, asm_opt.output_file_name);
+    // }
+    
     ul_re_correct(uidx, 3); 
     init_ulg_opt_t(&uu, uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_ul_tip, b_mask_t, is_trio);
     // print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug0", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1);
-    /**ul2ul_idx_t *u2o = **/gen_ul2ul(uidx, uopt, &uu, 0);
+    /**ul2ul_idx_t *u2o = **/gen_ul2ul(uidx, uopt, &uu, 0, is_bridg);
     // print_ul_alignment(init_ug, &UL_INF, 47072, "after-3");
 
     // print_debug_ul("UL.debug", init_ug, sg, uopt->coverage_cut, uopt->sources, uopt->ruIndex, bub, &UL_INF);
@@ -16804,11 +16929,12 @@ ul_renew_t *ropt, const char *bin_file, uint64_t free_uld)
     // resolve_dip_bub_chains(uidx);
 
     // free(r_het); destory_bubbles(bub); free(bub);
-    // uidx->uovl.hybrid_ug = gen_hybrid_ug(uidx, uidx->uovl.h_usg);
+    // if(free_uld) {
-    // print_debug_gfa(sg, uidx->uovl.hybrid_ug, uopt->coverage_cut, "hybrid_ug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1);
+    //     uidx->uovl.hybrid_ug = gen_hybrid_ug(uidx, uidx->uovl.h_usg);
-    // print_debug_gfa(sg, uidx->uovl.hybrid_ug, uopt->coverage_cut, "hybrid_ug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
+    //     // print_debug_gfa(sg, uidx->uovl.hybrid_ug, uopt->coverage_cut, "hybrid_ug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1);
-    // print_debug_gfa(sg, init_ug, uopt->coverage_cut, bin_file, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
+    //     print_debug_gfa(sg, uidx->uovl.hybrid_ug, uopt->coverage_cut, "hybrid_ug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
-    
+    //     // print_debug_gfa(sg, init_ug, uopt->coverage_cut, bin_file, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
+    // }
     // if(is_trio) gen_ul_trio_graph(uopt, uidx, o_file);    
     // exit(0);
     // return uidx->uovl.hybrid_ug;
@@ -16840,3 +16966,503 @@ ul_renew_t *ropt, const char *bin_file, uint64_t free_uld)
     // print_raw_uls_aln(uidx, asm_opt.output_file_name);
     // exit(0);
 }
+
+ug_clean_t *init_ug_clean_t(ug_opt_t *uopt, asg_t *sg, uint8_t is_ou, uint8_t is_trio, int64_t max_ext, double len_rat, double ou_rat, int64_t min_ou)
+{
+    uint64_t i;
+    ug_clean_t *sl; CALLOC(sl, 1);
+    ma_ug_t *ug = ma_ug_gen(sg);
+    for (i = 0; i < ug->g->n_seq; ++i) {
+        if(ug->g->seq[i].del) continue;
+        ug->g->seq[i].c = PRIMARY_LABLE;
+    }
+    MALLOC(sl->idx, (ug->g->n_seq<<1)); memset(sl->idx, -1, sizeof((*(sl->idx)))*(ug->g->n_seq<<1));
+    MALLOC(sl->bid, (ug->g->n_seq<<1)); memset(sl->bid, -1, sizeof((*(sl->bid)))*(ug->g->n_seq<<1));
+    sl->uopt = uopt; sl->sg = sg; sl->ug = ug; sl->idx_n = sl->bid_n = ug->g->n_seq<<1;
+    memset(&(sl->b), 0, sizeof(sl->b)); CALLOC(sl->b.a, (ug->g->n_seq<<1));
+    sl->is_ou = is_ou;
+    sl->is_trio = is_trio;
+    sl->max_ext = max_ext;
+    sl->len_rat = len_rat;
+    sl->ou_rat = ou_rat;
+    sl->min_ou = min_ou;
+    if(is_ou) update_ug_ou(sl->ug, sl->sg);
+    return sl;
+}
+
+void destroy_ug_clean_t(ug_clean_t *sl)
+{
+    free(sl->idx); free(sl->bid); ma_ug_destroy(sl->ug);
+    free(sl->b.a); free(sl->b.S.a); free(sl->b.T.a); free(sl->b.b.a); free(sl->b.e.a); 
+}
+
+void update_ug_clean_t(ug_clean_t *sl)
+{
+    uint32_t v, k, i, n_vtx = sl->ug->g->n_seq<<1;
+    sl->tlen = get_bub_pop_max_dist_advance(sl->ug->g, &(sl->b));
+    for (k = 0; k < sl->ug->g->n_seq; ++k) {
+        sl->idx[k<<1] = sl->idx[(k<<1)+1] = sl->bid[k<<1] = sl->bid[(k<<1)+1] = (uint32_t)-1;
+        if(sl->ug->g->seq[k].del) continue;
+        sl->ug->g->seq[k].c = PRIMARY_LABLE;
+    }
+
+    for (v = 0; v < n_vtx; ++v) {
+        if(sl->ug->g->seq[v>>1].del) continue;
+        if(asg_arc_n(sl->ug->g, v) < 2) continue;
+        if(sl->bid[v] != ((uint32_t)-1)) continue;
+        if(asg_bub_pop1_primary_trio(sl->ug->g, NULL, v, sl->tlen, &(sl->b), (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) {
+            //beg is v, end is b.S.a[0]
+            //note b.b include end, does not include beg
+            for (i = 0; i < sl->b.b.n; i++) {
+                if(sl->b.b.a[i]==v || sl->b.b.a[i]==sl->b.S.a[0]) continue;
+                sl->bid[sl->b.b.a[i]] = sl->bid[sl->b.b.a[i]^1] = 1;
+            }
+            sl->bid[v] = 2; sl->bid[sl->b.S.a[0]^1] = 3;
+        }
+    }
+
+    for (v = 0; v < n_vtx; ++v) {
+        if(sl->bid[v] !=2) continue;
+        if(asg_bub_pop1_primary_trio(sl->ug->g, NULL, v, sl->tlen, &(sl->b), (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) {
+            //note b.b include end, does not include beg
+            for (i = 0; i < sl->b.b.n; i++) {
+                if(sl->b.b.a[i]==v || sl->b.b.a[i]==sl->b.S.a[0]) continue;
+                sl->idx[sl->b.b.a[i]] = v;
+                sl->idx[sl->b.b.a[i]^1] = sl->b.S.a[0]^1;
+            }
+            sl->idx[v] = v; sl->idx[sl->b.S.a[0]^1] = sl->b.S.a[0]^1;
+        }
+    }
+    memcpy(sl->bid, sl->idx, sizeof((*(sl->bid)))*n_vtx);
+    // memset(sl->idx, -1, sizeof((*(sl->idx)))*n_vtx);
+}
+
+uint32_t get_best_het_node(ma_ug_t *ug, uint32_t v, uint32_t *bid, uint8_t is_ou, uint8_t is_trio, double len_rat, double ou_rat, uint32_t min_ou)
+{
+    uint32_t w = (uint32_t)-1, k, nv, kv, ol_max, ol_k;
+    asg_arc_t *av; uint32_t trioF = (uint32_t)-1, ntrioF = (uint32_t)-1;
+    if(ug->g->seq[v>>1].del) return (uint32_t)-1;
+    av = asg_arc_a(ug->g, v); nv = asg_arc_n(ug->g, v);
+    for (k = kv = 0; k < nv; k++) {
+        if(av[k].del) continue;
+        ///could not connect to the beg/sink node
+        if((av[k].v>>1)==(bid[v]>>1) || (av[k].v>>1)==(bid[v^1]>>1)) return (uint32_t)-1;
+        kv++; w = av[k].v;
+    }
+    if(kv < 2) return w;///kv == 0, return (uint32_t)-1; kv == 1, return node; 
+    // if(v == 351) fprintf(stderr, "-0-[M::%s]\tutg%.6ul(%c)\tv::%u\tkv::%u\n", __func__, (v>>1)+1, "+-"[v&1], v, kv);
+
+    if(is_trio) {
+        trioF = get_ug_tip_trio_infor(ug, v^1);
+        ntrioF = (trioF==FATHER? MOTHER : (trioF==MOTHER? FATHER : (uint32_t)-1));
+    }
+
+    ol_max = 0; ol_k = (uint32_t)-1;
+    for (k = 0; k < nv; ++k) {
+        if(av[k].del) continue;
+        // if(is_trio && get_ug_tip_trio_infor(ug, av[k].v) == ntrioF) continue;
+        if(ol_max < av[k].ol) ol_max = av[k].ol, ol_k = k;
+    }
+    // if(v == 351) fprintf(stderr, "-1-[M::%s]\tutg%.6ul(%c)\tv::%u\tol_k::%u\n", __func__, (v>>1)+1, "+-"[v&1], v, ol_k);
+    if(ol_k == (uint32_t)-1) return (uint32_t)-1;
+    for (k = 0; k < nv; ++k) {
+        if(av[k].del || k == ol_k) continue;
+        // if(is_trio && get_ug_tip_trio_infor(ug, av[k].v) == ntrioF) continue;
+        if(av[k].ol > av[ol_k].ol*len_rat) return (uint32_t)-1;
+        if((is_ou) && (av[k].ou > min_ou) && ((av[k].ou) > (av[ol_k].ou*ou_rat))) return (uint32_t)-1;
+    }
+    if(is_trio && get_ug_tip_trio_infor(ug, av[ol_k].v) == ntrioF) return (uint32_t)-1;
+    // if(v == 351) fprintf(stderr, "-2-[M::%s]\tutg%.6ul(%c)\tv::%u\tav[ol_k].v::%u\n", __func__, (v>>1)+1, "+-"[v&1], v, av[ol_k].v);
+    return av[ol_k].v;
+}
+
+static void cal_bub_best(void *data, long i, int tid)
+{
+    ug_clean_t *sl = (ug_clean_t *)data;
+    ma_ug_t *ug = sl->ug; uint32_t v, w;
+    sl->idx[i<<1] = sl->idx[(i<<1)+1] = (uint32_t)-1;
+    if(ug->g->seq[i].del) return;
+    if((sl->bid[i<<1] == (uint32_t)-1) || (sl->bid[(i<<1)+1] == (uint32_t)-1)) return;///not within a bubble
+    
+    v = i<<1;
+    if(sl->bid[v] != v) {///not the beg/sink node
+        w = get_best_het_node(ug, v, sl->bid, sl->is_ou, sl->is_trio, sl->len_rat, sl->ou_rat, sl->min_ou);
+        // if(v == 352) fprintf(stderr, "-v-[M::%s]\tutg%.6ul(%c)\tv::%u\tw::%u\n", __func__, (v>>1)+1, "+-"[v&1], v, w);
+        if((w != (uint32_t)-1) && ((v^1) == get_best_het_node(ug, w^1, sl->bid, sl->is_ou, sl->is_trio, sl->len_rat, sl->ou_rat, sl->min_ou))) {
+            sl->idx[v] = w;
+        }
+    }
+
+    v = (i<<1)+1;
+    if(sl->bid[v] != v) {///not the beg/sink node
+        w = get_best_het_node(ug, v, sl->bid, sl->is_ou, sl->is_trio, sl->len_rat, sl->ou_rat, sl->min_ou);
+        if((w != (uint32_t)-1) && ((v^1) == get_best_het_node(ug, w^1, sl->bid, sl->is_ou, sl->is_trio, sl->len_rat, sl->ou_rat, sl->min_ou))) {
+            sl->idx[v] = w;
+        }
+    }
+}
+
+uint32_t usg_topocut_aux_unambi1(asg_t *g, uint32_t v)
+{
+	asg_arc_t *av = asg_arc_a(g, v);
+	uint32_t i, nv = asg_arc_n(g, v);
+	uint32_t k = nv, kv;
+	for (i = 0, kv = 0; i < nv; ++i)
+		if (!av[i].del) ++kv, k = i;
+	if (kv != 1) return (uint32_t)-1;
+	return av[k].v;
+}
+
+int usg_topocut_aux_del(ma_ug_t *ug, uint32_t v, int max_ext, asg64_v *b)
+{
+	int32_t n_ext; 
+	for (n_ext = 1; n_ext < max_ext && v != (uint32_t)-1; ++n_ext) {
+		if (usg_topocut_aux_unambi1(ug->g, v^1) == (uint32_t)-1) {
+			--n_ext;
+			break;
+		}
+        kv_push(uint64_t, *b, v);
+		v = usg_topocut_aux_unambi1(ug->g, v);
+	}
+	return n_ext;
+}
+
+#define is_best_arc(sl, v, w) (((sl).idx[(v)]==(w))&&((sl).idx[((w)^1)]==((v)^1)))
+
+inline void asg_arc_del_by_ug(asg_t *sg, ma_ug_t *ug, uint32_t uv, uint32_t uw, uint32_t del)
+{
+    uint32_t rv, rw, i, nv; asg_arc_t *av;
+    if(uv&1) rv = ug->u.a[uv>>1].start^1;
+    else rv = ug->u.a[uv>>1].end^1;
+
+    if(uw&1) rw = ug->u.a[uw>>1].end;
+    else rw = ug->u.a[uw>>1].start;
+
+    av = asg_arc_a(sg, rv);
+    nv = asg_arc_n(sg, rv);
+    for (i = 0; i < nv; i++) {
+        if(av[i].v == rw) break;
+    }
+    assert(i < nv);
+    av[i].del = del;
+}
+
+void cal_bub_best_by_len(ug_clean_t *sl, asg64_v *in, uint32_t max_ext, uint32_t is_trio, uint32_t is_ou, double len_rat, double ou_rat, uint32_t min_ou)
+{
+    // fprintf(stderr, "[M::%s]\tStart\n", __func__);
+    ma_ug_t *ug = sl->ug; asg_t *g = sl->ug->g; uint32_t ol_max, ou_max;
+    uint32_t v, w, n_vtx = (g->n_seq<<1), nv, nw, i, k, z, kv, kw, bb, to_del, bn, tip;
+    asg64_v tx = {0,0,0}, *b = NULL; asg_arc_t *av, *aw, *ve, *we; ma_utg_t *u;
+    uint32_t trioF = (uint32_t)-1, ntrioF = (uint32_t)-1, mm_ol, mm_ou, cnt = 0, del_v, del_w;
+
+    if(in) b = in;
+    else b = &tx;
+    b->n = 0;
+
+    for (v = 0; v < n_vtx; ++v) {
+        if (g->seq[v>>1].del) continue;
+        av = asg_arc_a(g, v); 
+        nv = asg_arc_n(g, v);
+        if (nv < 2) continue;
+
+        for (i = kv = bb = 0; i < nv; ++i) {
+            if(av[i].del) continue;
+            if(is_best_arc((*sl), v, av[i].v)) bb++;
+            kv++;
+        }
+        if((kv < 2) || (!bb) || (kv<=bb)) continue;///it is impossible that kv <= bv
+
+        for (i = 0; i < nv; ++i) {
+            if(av[i].del) continue;
+            if(is_best_arc((*sl), v, av[i].v)) continue;
+            kv_push(uint64_t, *b, (((uint64_t)av[i].ol)<<32) | ((uint64_t)(av-g->arc+i)));   
+        }
+    }
+
+    radix_sort_srt64(b->a, b->a + b->n);
+    for (k = 0; k < b->n; k++) {
+        if(g->arc[(uint32_t)b->a[k]].del) continue;
+
+        v = g->arc[(uint32_t)b->a[k]].ul>>32; w = g->arc[(uint32_t)b->a[k]].v^1;
+        if(g->seq[v>>1].del || g->seq[w>>1].del) continue;
+        nv = asg_arc_n(g, v); nw = asg_arc_n(g, w);
+        av = asg_arc_a(g, v); aw = asg_arc_a(g, w);
+        if(nv<=1 && nw <= 1) continue;
+
+        if(is_trio) {
+            if(get_arcs(g, v, NULL, 0)<=1 && get_arcs(g, w, NULL, 0)<=1) continue;///speedup
+            trioF = get_ug_tip_trio_infor(ug, v^1);
+            ntrioF = (trioF==FATHER? MOTHER : (trioF==MOTHER? FATHER : (uint32_t)-1));
+        }
+
+        ve = &(g->arc[(uint32_t)b->a[k]]);
+        for (i = 0; i < nw; ++i) {
+            if (aw[i].v == (v^1)) {
+                we = &(aw[i]);
+                break;
+            }
+        }
+        ///mm_ol and mm_ou are used to make edge with long indel more easy to be cutted
+        mm_ol = MIN(ve->ol, we->ol); mm_ou = MIN(ve->ou, we->ou);
+
+        for (i = kv = ol_max = ou_max = 0; i < nv; ++i) {
+            if(av[i].del) continue;
+            kv++; 
+            if(is_trio && get_ug_tip_trio_infor(ug, av[i].v) == ntrioF) continue;
+            if(ol_max < av[i].ol) ol_max = av[i].ol; 
+            if(ou_max < av[i].ou) ou_max = av[i].ou;
+        }
+        if (kv < 1) continue;
+        if (kv >= 2) {
+            if (mm_ol > ol_max*len_rat) continue;
+            if ((is_ou) && (mm_ou > min_ou) && (mm_ou > (ou_max*ou_rat))) continue;
+        }
+        
+
+        for (i = kw = ol_max = ou_max = 0; i < nw; ++i) {
+            if(aw[i].del) continue;
+            kw++; 
+            if(is_trio && get_ug_tip_trio_infor(ug, aw[i].v) == ntrioF) continue;
+            if(ol_max < aw[i].ol) ol_max = aw[i].ol; 
+            if(ou_max < aw[i].ou) ou_max = aw[i].ou;
+        }
+        if (kw < 1) continue;
+        if (kw >= 2) {
+            if (mm_ol > ol_max*len_rat) continue;
+            if ((is_ou) && (mm_ou > min_ou) && (mm_ou > (ou_max*ou_rat))) continue;
+        }
+
+        if (kv <= 1 && kw <= 1) continue;
+
+        to_del = 0; del_v = del_w = (uint32_t)-1; tip = 0;
+        if (kv > 1 && kw > 1) {
+            to_del = 1;
+        } else if (kw == 1) {
+            tip = asg_topocut_aux(g, w^1, max_ext);
+            if (tip < max_ext) {
+                to_del = 1; del_w = w^1;
+            }
+        } else if (kv == 1) {
+            tip = asg_topocut_aux(g, v^1, max_ext);
+            if (tip < max_ext) {
+                to_del = 1; del_v = v^1;
+            }
+        }
+
+        if (to_del) {
+            bn = b->n;
+            if(del_v != ((uint32_t)-1)) usg_topocut_aux_del(ug, del_v, max_ext, b);
+            if(del_w != ((uint32_t)-1)) usg_topocut_aux_del(ug, del_w, max_ext, b);
+            for (i = bn; i < b->n; i++) {
+                u = &(ug->u.a[b->a[i]>>1]);
+                if(u->m == 0) continue;
+                for (z = 0; z < u->n; z++) asg_seq_del(sl->sg, u->a[z]>>33);
+                asg_seq_del(ug->g, (b->a[i]>>1));
+                if(u->m) {
+                    u->m = u->n = 0; free(u->a); u->a = NULL;
+                }
+            }
+            // assert(tip == (b->n-bn));
+            b->n = bn;
+
+            ve->del = we->del = 1, ++cnt;
+            asg_arc_del_by_ug(sl->sg, ug, ve->ul>>32, ve->v, 1);
+            asg_arc_del_by_ug(sl->sg, ug, we->ul>>32, we->v, 1);
+            
+        }
+    }
+
+    if(!in) free(tx.a);
+    if (cnt > 0) {
+        asg_cleanup(g); asg_cleanup(sl->sg);
+    }
+    // fprintf(stderr, "[M::%s]\tEnd\n", __func__);
+}
+
+/**
+void cal_bub_best_by_topo(ug_clean_t *sl, asg64_v *in, uint32_t max_ext, uint32_t is_trio, uint32_t is_ou, double len_rat, double ou_rat, uint32_t min_ou, uint32_t long_tip)
+{
+    // fprintf(stderr, "[M::%s]\tStart\n", __func__);
+    ma_ug_t *ug = sl->ug; asg_t *g = sl->ug->g; uint32_t ol_max, ou_max;
+    uint32_t v, w, n_vtx = (g->n_seq<<1), nv, nw, i, k, z, kv, kw, to_del, bn, tip, avi, awi;
+    asg64_v tx = {0,0,0}, *b = NULL; asg_arc_t *av, *aw, *ve, *we; ma_utg_t *u;
+    uint32_t trioF = (uint32_t)-1, mm_ol, mm_ou, cnt = 0, del_v, del_w;
+
+    if(in) b = in;
+    else b = &tx;
+    b->n = 0;
+
+    for (v = 0; v < n_vtx; ++v) {
+        if (g->seq[v>>1].del) continue;
+        if((sl->bid[v] == (uint32_t)-1) || (sl->bid[v^1] == (uint32_t)-1)) continue;///not within a bubble
+
+        av = asg_arc_a(g, v); 
+        nv = asg_arc_n(g, v);
+        if (nv < 2) continue;
+
+        for (i = kv = 0; i < nv; ++i) {
+            if(av[i].del) continue;
+            ///could not connect to the beg/sink node
+            if((av[k].v>>1)==(sl->bid[v]>>1) || (av[k].v>>1)==(sl->bid[v^1]>>1)) break;
+            kv++;
+        }
+        if(kv < 2 || i < nv) continue;///it is impossible that kv <= bv
+
+        for (i = 0; i < nv; ++i) {
+            if(av[i].del) continue;
+            kv_push(uint64_t, *b, (((uint64_t)av[i].ol)<<32) | ((uint64_t)(av-g->arc+i)));   
+        }
+    }
+
+    radix_sort_srt64(b->a, b->a + b->n);
+    for (k = 0; k < b->n; k++) {
+        if(g->arc[(uint32_t)b->a[k]].del) continue;
+
+        v = g->arc[(uint32_t)b->a[k]].ul>>32; w = g->arc[(uint32_t)b->a[k]].v^1;
+        if(g->seq[v>>1].del || g->seq[w>>1].del) continue;
+        nv = asg_arc_n(g, v); nw = asg_arc_n(g, w);
+        av = asg_arc_a(g, v); aw = asg_arc_a(g, w);
+        if(nv < 2 || nw < 2) continue;
+        kv = get_arcs(g, v, NULL, 0); kw = get_arcs(g, w, NULL, 0);
+        if(kv < 2 || kw < 2) continue;
+
+        if(is_trio) {
+            trioF = get_ug_tip_trio_infor(ug, v^1);
+            if(trioF == FATHER || trioF == MOTHER) {
+                if(get_ug_tip_trio_infor(ug, w^1) == trioF) continue;
+            }
+        }
+
+        avi = awi = (uint32_t)-1;
+
+        avi = ((uint32_t)b->a[k]) - ((uint64_t)(av-g->arc));
+        ve = &(g->arc[(uint32_t)b->a[k]]);
+
+        for (i = 0; i < nw; ++i) {
+            if (aw[i].v == (v^1)) {
+                we = &(aw[i]); awi = i;
+                break;
+            }
+        }
+
+        for (i = kv = ol_max = ou_max = 0; i < nv; ++i) {
+            if(av[i].del) continue;
+        }
+
+        ///mm_ol and mm_ou are used to make edge with long indel more easy to be cutted
+        mm_ol = MIN(ve->ol, we->ol); mm_ou = MIN(ve->ou, we->ou);
+
+        for (i = kv = ol_max = ou_max = 0; i < nv; ++i) {
+            if(av[i].del) continue;
+            kv++; 
+            if(is_trio && get_ug_tip_trio_infor(ug, av[i].v) == ntrioF) continue;
+            if(ol_max < av[i].ol) ol_max = av[i].ol; 
+            if(ou_max < av[i].ou) ou_max = av[i].ou;
+        }
+        if (kv < 1) continue;
+        if (kv >= 2) {
+            if (mm_ol > ol_max*len_rat) continue;
+            if ((is_ou) && (mm_ou > min_ou) && (mm_ou > (ou_max*ou_rat))) continue;
+        }
+        
+
+        for (i = kw = ol_max = ou_max = 0; i < nw; ++i) {
+            if(aw[i].del) continue;
+            kw++; 
+            if(is_trio && get_ug_tip_trio_infor(ug, aw[i].v) == ntrioF) continue;
+            if(ol_max < aw[i].ol) ol_max = aw[i].ol; 
+            if(ou_max < aw[i].ou) ou_max = aw[i].ou;
+        }
+        if (kw < 1) continue;
+        if (kw >= 2) {
+            if (mm_ol > ol_max*len_rat) continue;
+            if ((is_ou) && (mm_ou > min_ou) && (mm_ou > (ou_max*ou_rat))) continue;
+        }
+
+        if (kv <= 1 && kw <= 1) continue;
+
+        to_del = 0; del_v = del_w = (uint32_t)-1; tip = 0;
+        if (kv > 1 && kw > 1) {
+            to_del = 1;
+        } else if (kw == 1) {
+            tip = asg_topocut_aux(g, w^1, max_ext);
+            if (tip < max_ext) {
+                to_del = 1; del_w = w^1;
+            }
+        } else if (kv == 1) {
+            tip = asg_topocut_aux(g, v^1, max_ext);
+            if (tip < max_ext) {
+                to_del = 1; del_v = v^1;
+            }
+        }
+
+        if (to_del) {
+            bn = b->n;
+            if(del_v != ((uint32_t)-1)) usg_topocut_aux_del(ug, del_v, max_ext, b);
+            if(del_w != ((uint32_t)-1)) usg_topocut_aux_del(ug, del_w, max_ext, b);
+            for (i = bn; i < b->n; i++) {
+                u = &(ug->u.a[b->a[i]>>1]);
+                if(u->m == 0) continue;
+                for (z = 0; z < u->n; z++) asg_seq_del(sl->sg, u->a[z]>>33);
+                asg_seq_del(ug->g, (b->a[i]>>1));
+                if(u->m) {
+                    u->m = u->n = 0; free(u->a); u->a = NULL;
+                }
+            }
+            // assert(tip == (b->n-bn));
+            b->n = bn;
+
+            ve->del = we->del = 1, ++cnt;
+            asg_arc_del_by_ug(sl->sg, ug, ve->ul>>32, ve->v, 1);
+            asg_arc_del_by_ug(sl->sg, ug, we->ul>>32, we->v, 1);
+            
+        }
+    }
+
+    if(!in) free(tx.a);
+    if (cnt > 0) {
+        asg_cleanup(g); asg_cleanup(sl->sg);
+    }
+    // fprintf(stderr, "[M::%s]\tEnd\n", __func__);
+}
+**/
+
+
+
+void deep_graph_clean(ug_opt_t *uopt, asg_t *sg, uint8_t is_ou, uint8_t is_trio, int64_t max_ext, 
+double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, double ou_rat, int64_t min_ou, int64_t clean_round, int64_t long_tip)
+{
+    ug_clean_t *sl; asg64_v b; double step; int64_t i;
+    kv_init(b); hic_clean_adv(sg, uopt);
+    step = (clean_round==1?max_ovlp_drop_ratio:((max_ovlp_drop_ratio-min_ovlp_drop_ratio)/(clean_round-1)));
+    sl = init_ug_clean_t(uopt, sg, is_ou, is_trio, max_ext, min_ovlp_drop_ratio, ou_rat, min_ou);
+
+    // print_debug_gfa(sg, NULL, uopt->coverage_cut, "debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
+    for (i = 0, sl->len_rat = min_ovlp_drop_ratio; i < clean_round; i++, sl->len_rat += step) {
+        if(sl->len_rat > max_ovlp_drop_ratio) sl->len_rat = max_ovlp_drop_ratio;
+        update_ug_clean_t(sl);
+        kt_for(asm_opt.thread_num, cal_bub_best, sl, sl->ug->g->n_seq);
+        cal_bub_best_by_len(sl, &b, max_ext, is_trio, sl->is_ou, sl->len_rat, sl->ou_rat, min_ou);
+    }
+
+    if(is_ou) {
+        uint64_t k; sl->is_ou = 0;
+        for (k = 0; k < sl->ug->g->n_arc; k++) sl->ug->g->arc[k].ou = 0;
+        max_ovlp_drop_ratio = 0.6;
+        if(max_ovlp_drop_ratio > min_ovlp_drop_ratio) {
+            step = (clean_round==1?max_ovlp_drop_ratio:((max_ovlp_drop_ratio-min_ovlp_drop_ratio)/(clean_round-1)));
+            for (i = 0, sl->len_rat = min_ovlp_drop_ratio; i < clean_round; i++, sl->len_rat += step) {
+                if(sl->len_rat > max_ovlp_drop_ratio) sl->len_rat = max_ovlp_drop_ratio;
+                update_ug_clean_t(sl);
+                kt_for(asm_opt.thread_num, cal_bub_best, sl, sl->ug->g->n_seq);
+                cal_bub_best_by_len(sl, &b, max_ext, is_trio, sl->is_ou, sl->len_rat, sl->ou_rat, min_ou);
+            }
+        }
+    }
+
+    // update_ug_clean_t(sl);
+    // cal_bub_best_by_topo(sl, &b, max_ext, is_trio, sl->is_ou, sl->len_rat, sl->ou_rat, min_ou, long_tip);
+
+    kv_destroy(b); destroy_ug_clean_t(sl); free(sl);
+    hic_clean_adv(sg, uopt);
+}

gfa_ut.h

@@ -27,7 +27,7 @@ void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float o
 uint32_t asg_cut_large_indel(asg_t *g, asg64_v *in, int32_t max_ext, float ou_rat, uint32_t is_ou);
 uint32_t asg_cut_semi_circ(asg_t *g, uint32_t lim_len, uint32_t is_clean);
 void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, 
-double max_ovlp_drop_ratio, int64_t max_tip, int64_t max_ul_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, ul_renew_t *ropt, const char *bin_file, uint64_t free_uld);
+double max_ovlp_drop_ratio, int64_t max_tip, int64_t max_ul_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, ul_renew_t *ropt, const char *bin_file, uint64_t free_uld, uint64_t is_bridg, uint64_t deep_clean);
 void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t max_hang, int64_t min_ovlp, int64_t ul_occ);
 void normalize_gou(asg_t *g);
 void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd);

horder.cpp

@@ -842,6 +842,13 @@ void horder_clean_sg_by_utg(asg_t *sg, ma_ug_t *ug)
             wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32));
             asg_arc_unique_del(sg, vx, wx, 0); asg_arc_unique_del(sg, wx^1, vx^1, 0);
         }
+
+        if(u->circ) {
+            v = w = i<<1;
+            vx = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32));
+            wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32));
+            asg_arc_unique_del(sg, vx, wx, 0); asg_arc_unique_del(sg, wx^1, vx^1, 0);
+        }         
     }
     asg_cleanup(sg);
 

inter.h

@@ -118,5 +118,11 @@ void update_ug_arch_ul_mul(ma_ug_t *ug);
 void print_ul_alignment(ma_ug_t *ug, all_ul_t *aln, uint32_t id, const char* cmd);
 void clear_all_ul_t(all_ul_t *x);
 void trans_base_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res, bubble_type *bub);
+hpc_re_t *gen_hpc_re_t(ma_ug_t *ug);
+idx_emask_t* graph_ovlp_binning(ma_ug_t *ug, asg_t *sg, const ug_opt_t *uopt);
+uint32_t gen_src_shared_interval_simple(uint32_t src, ma_ug_t *ug, kv_ul_ov_t *res);
+uint64_t check_ul_ov_t_consist(ul_ov_t *x, ul_ov_t *y, int64_t ql, int64_t tl, double diff);
+uint32_t infer_se(uint32_t qs, uint32_t qe, uint32_t ts, uint32_t te, uint32_t rev, 
+uint32_t rqs, uint32_t rqe, uint32_t *rts, uint32_t *rte);
 
 #endif