more accurate purging

Overlaps.cpp

@@ -11717,15 +11717,20 @@ kvec_asg_arc_t_warp* new_rtg_edges, int max_hang, int min_ovlp)
     uint8_t* primary_flag = (uint8_t*)calloc(sg->n_seq, sizeof(uint8_t));
     hap_cov_t *cov = init_hap_cov_t(ug, sg, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, 0);
 
-    int tmp_cov = asm_opt.hom_global_coverage;
+    int is_set = ((asm_opt.hom_global_coverage != -1)? 1 : 0);
-    asm_opt.hom_global_coverage = -1;
 
-    purge_dups(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, 
+    if(is_set == 0)
+    {
+        purge_dups(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, 
         asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, 0, 0, 1, cov, 0);
-    dip_thre_max = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE)*0.70;
+        dip_thre_max = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE);
-    asm_opt.hom_global_coverage = tmp_cov;
+    }
-    ///fprintf(stderr, "dip_thre_max: %lu\n", dip_thre_max);
+    else
-    
+    {
+        dip_thre_max = asm_opt.hom_global_coverage;
+    }
+    dip_thre_max *= 0.70;
+            
     for (i = 0; i < ug->g->n_seq; i++)
     {
         dip_thres = dip_thre_max;
@@ -11749,6 +11754,96 @@ kvec_asg_arc_t_warp* new_rtg_edges, int max_hang, int min_ovlp)
     ///fprintf(stderr, "[M::%s] diploid coverage threshold: %lu\n", __func__, dip_thres);
 }
 
+void update_hc_links_by_trans_chain(trans_chain* t_ch)
+{
+    if(t_ch == NULL) return;
+    ///fprintf(stderr, "sbsbsbsbsbsb1sbsbsbsbsbsb, l0_chain: %u, chain_num: %u\n", t_ch->l0_chain, t_ch->chain_num);
+    uint32_t i, k, *x = NULL, x_occ, *y = NULL, y_occ, x_k, y_k;
+    memset(t_ch->is_het, 0, t_ch->u_num);
+    for (i = 0; i < t_ch->l0_chain; i++)
+    {
+        x_occ = y_occ = 0;
+        get_chain_trans(t_ch, i, &x, &x_occ, &y, &y_occ);
+
+        for (k = x_k = 0; k < x_occ; k++)
+        {
+            if(x[k] == (uint32_t)-1) continue;
+            x_k++;
+        }
+        for (k = y_k = 0; k < y_occ; k++)
+        {
+            if(y[k] == (uint32_t)-1) continue;
+            y_k++;
+        }
+
+        if(x_k == 0 || y_k == 0)
+        {
+            for (k = 0; k < x_occ; k++) x[k] = (uint32_t)-1;
+            for (k = 0; k < y_occ; k++) y[k] = (uint32_t)-1;
+        }
+        else
+        {
+            for (k = 0; k < x_occ; k++)
+            {
+                if(x[k] == (uint32_t)-1) continue;
+                t_ch->is_het[x[k]>>1] = 1;
+            }
+            for (k = 0; k < y_occ; k++)
+            {
+                if(y[k] == (uint32_t)-1) continue;
+                t_ch->is_het[y[k]>>1] = 1;
+            }
+        }
+    }
+
+    /*******************************for debug************************************/
+    // for (i = 0; i < t_ch->l0_chain; i++)
+    // {
+    //     x_occ = y_occ = 0;
+    //     get_chain_trans(t_ch, i, &x, &x_occ, &y, &y_occ);
+    //     fprintf(stderr, "\n******x******");
+    //     for (k = 0; k < x_occ; k++)
+    //     {
+    //         fprintf(stderr, "utg%.6ul\t%u\n", (x[k]>>1)+1, x[k]&1);
+    //     }
+    //     fprintf(stderr, "******y******");
+    //     for (k = 0; k < y_occ; k++)
+    //     {
+    //         fprintf(stderr, "utg%.6ul\t%u\n", (y[k]>>1)+1, y[k]&1);
+    //     }
+    // }
+    // exit(0);
+    /*******************************for debug************************************/
+}
+
+
+void check_utg_het(ma_ug_t *ug, asg_t *sg, 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* new_rtg_edges, int max_hang, int min_ovlp, hap_cov_t *cov)
+{
+    ///set by cov->t_ch
+    update_hc_links_by_trans_chain(cov->t_ch);
+    bubble_type bub;
+    memset(&bub, 0, sizeof(bubble_type));
+    ///set by coverage
+    classify_untigs(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, max_hang, min_ovlp);
+    ///set by bubble
+    identify_bubbles(ug, &bub, cov->t_ch->is_het);
+    uint32_t i;
+    for (i = 0; i < ug->g->n_seq; i++)
+    {
+        if(IF_HOM(i, bub))
+        {
+            cov->t_ch->is_het[i] = 0;
+        }
+        else
+        {
+            cov->t_ch->is_het[i] = 1;
+        }
+    }
+    destory_bubbles(&bub);
+}
+
 
 trans_chain* init_trans_chain(ma_ug_t *ug, uint64_t r_num)
 {
@@ -11988,16 +12083,18 @@ bub_label_t* b_mask_t)
 
 
 
-    hap_cov_t *cov = init_hap_cov_t(ug, sg, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, 1);
+    hap_cov_t *cov = NULL;
     asg_t *copy_sg = copy_read_graph(sg);
     ma_ug_t *copy_ug = copy_untig_graph(ug);    
     ///asm_opt.purge_overlap_len = asm_opt.purge_overlap_len_hic;
     ///asm_opt.purge_simi_thres = asm_opt.purge_simi_rate_hic;
     adjust_utg_by_primary(&copy_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, 
     tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, 
-    max_hang, min_ovlp, &new_rtg_edges, cov, b_mask_t, 1);
+    max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1);
     ma_ug_destroy(copy_ug);
     asg_destroy(copy_sg);
+
+
     ma_ug_print_bed(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, 
             max_hang, min_ovlp, asm_opt.hic_inconsist_rate, NULL, NULL, cov);
 
@@ -14515,7 +14612,10 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
 {
     asg_t* nsg = (*ug)->g;
     uint32_t v, n_vtx = nsg->n_seq;    
-    hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, 0);
+    ma_ug_t *cov_ug = NULL;
+    if(asm_opt.purge_level_trio > 0) cov_ug = copy_untig_graph(*ug);
+    hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources, 
+                                                    coverage_cut, max_hang, min_ovlp, asm_opt.purge_level_trio>0?1:0);
 
     purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, 
         asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 
@@ -14555,8 +14655,12 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
     clean_trio_untig_graph(*ug, read_g, coverage_cut, sources, reverse_sources, tipsLen, 
     tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, chimeric_rate, 0, 0, drop_ratio, flag, drop_rate, cov);
 
-    ///if(flag == MOTHER) fprintf(stderr, "(o.1) c: %u, del: %u, n: %u\n", (*ug)->g->seq[28141].c, (*ug)->g->seq[28141].del, (*ug)->u.a[28141].n);
+    if(cov_ug)
-    
+    {
+        check_utg_het(cov_ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, 
+        new_rtg_edges, max_hang, min_ovlp, cov);
+    }
+
     ///delete_useless_nodes(ug);
     delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex);
 
@@ -14567,6 +14671,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
 
     renew_utg(ug, read_g, new_rtg_edges);
 
+
     if (!(asm_opt.flag & HA_F_BAN_POST_JOIN))
     {
         rescue_missing_overlaps_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang,
@@ -22299,86 +22404,23 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, trans_chain* t_ch)
     ma_ug_destroy(atg);
 }
 
-void update_hc_links_by_trans_chain(trans_chain* t_ch)
-{
-    ///fprintf(stderr, "sbsbsbsbsbsb1sbsbsbsbsbsb, l0_chain: %u, chain_num: %u\n", t_ch->l0_chain, t_ch->chain_num);
-    uint32_t i, k, *x = NULL, x_occ, *y = NULL, y_occ, x_k, y_k;
-    memset(t_ch->is_het, 0, t_ch->u_num);
-    for (i = 0; i < t_ch->l0_chain; i++)
-    {
-        x_occ = y_occ = 0;
-        get_chain_trans(t_ch, i, &x, &x_occ, &y, &y_occ);
-
-        for (k = x_k = 0; k < x_occ; k++)
-        {
-            if(x[k] == (uint32_t)-1) continue;
-            x_k++;
-        }
-        for (k = y_k = 0; k < y_occ; k++)
-        {
-            if(y[k] == (uint32_t)-1) continue;
-            y_k++;
-        }
-
-        if(x_k == 0 || y_k == 0)
-        {
-            for (k = 0; k < x_occ; k++) x[k] = (uint32_t)-1;
-            for (k = 0; k < y_occ; k++) y[k] = (uint32_t)-1;
-        }
-        else
-        {
-            for (k = 0; k < x_occ; k++)
-            {
-                if(x[k] == (uint32_t)-1) continue;
-                t_ch->is_het[x[k]>>1] = 1;
-            }
-            for (k = 0; k < y_occ; k++)
-            {
-                if(y[k] == (uint32_t)-1) continue;
-                t_ch->is_het[y[k]>>1] = 1;
-            }
-        }
-    }
-
-    /*******************************for debug************************************/
-    // for (i = 0; i < t_ch->l0_chain; i++)
-    // {
-    //     x_occ = y_occ = 0;
-    //     get_chain_trans(t_ch, i, &x, &x_occ, &y, &y_occ);
-    //     fprintf(stderr, "\n******x******");
-    //     for (k = 0; k < x_occ; k++)
-    //     {
-    //         fprintf(stderr, "utg%.6ul\t%u\n", (x[k]>>1)+1, x[k]&1);
-    //     }
-    //     fprintf(stderr, "******y******");
-    //     for (k = 0; k < y_occ; k++)
-    //     {
-    //         fprintf(stderr, "utg%.6ul\t%u\n", (y[k]>>1)+1, y[k]&1);
-    //     }
-    // }
-    // exit(0);
-    /*******************************for debug************************************/
-}
-
 
 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 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, hap_cov_t *i_cov, bub_label_t* b_mask_t, 
+kvec_asg_arc_t_warp* new_rtg_edges, hap_cov_t **i_cov, bub_label_t* b_mask_t, 
 uint32_t collect_p_trans)
 {
     asg_t* nsg = (*ug)->g;
     uint32_t v, n_vtx = nsg->n_seq, k, rId, just_contain;
     ma_utg_t* u = NULL;
     hap_cov_t *cov = NULL;
-    if(i_cov == NULL) cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, 0);
+    ma_ug_t *cov_ug = NULL;
-    else cov = i_cov;
+    if(asm_opt.purge_level_primary > 0) cov_ug = copy_untig_graph(*ug);
+    cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources, 
+                                                coverage_cut, max_hang, min_ovlp, asm_opt.purge_level_primary>0?1:0);
 
-    ///print_utg_coverage(*ug, coverage_cut, 440, sources);
-    ///exit(0);    
-    // drop_semi_circle((*ug), nsg, read_g, reverse_sources, ruIndex);
-    // asg_cleanup(nsg);
     adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t);
     
     nsg = (*ug)->g;
@@ -22394,7 +22436,12 @@ uint32_t collect_p_trans)
     delete_useless_nodes(ug);
     renew_utg(ug, read_g, new_rtg_edges);
 
-    if(i_cov) update_hc_links_by_trans_chain(cov->t_ch);
+    if(cov_ug)
+    {
+        check_utg_het(cov_ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, 
+        new_rtg_edges, max_hang, min_ovlp, cov);
+    }
+    
     if(i_cov) goto skip_purge;
 
     if(asm_opt.purge_level_primary > 0)
@@ -22490,16 +22537,25 @@ uint32_t collect_p_trans)
     recover_utg_by_coverage(ug, read_g, coverage_cut, sources, ruIndex, cov->t_ch);
     if(i_cov)
     {
-        update_hc_links_by_trans_chain(cov->t_ch);
+        if(asm_opt.purge_level_primary > 0)
-        if(collect_p_trans)
         {
-            purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, 
+            if(collect_p_trans)
-            asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 
+            {
-            0, 0, cov, collect_p_trans);
+                // purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, 
+                // asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 
+                // 0, 0, cov, collect_p_trans);
+            }
         }
+        
+        update_hc_links_by_trans_chain(cov->t_ch);
+        (*i_cov) = cov;
+    } 
+    else
+    {
+        destory_hap_cov_t(&cov);
     } 
     
-    else destory_hap_cov_t(&cov);
+    ma_ug_destroy(cov_ug);
 }
 
 
@@ -26514,12 +26570,12 @@ bub_label_t* b_mask_t)
     ma_ug_t *ug = NULL;
     ug = ma_ug_gen_primary(sg, PRIMARY_LABLE);
 
-    hap_cov_t *cov = init_hap_cov_t(ug, sg, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, 1);
+    hap_cov_t *cov = NULL;
     asg_t *copy_sg = copy_read_graph(sg);
     ma_ug_t *copy_ug = copy_untig_graph(ug);    
     adjust_utg_by_primary(&copy_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, 
     tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, 
-    max_hang, min_ovlp, &new_rtg_edges, cov, b_mask_t, 0);
+    max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 0);
     ma_ug_destroy(copy_ug); copy_ug = NULL;
     asg_destroy(copy_sg); copy_sg = NULL;
 

Overlaps.h

@@ -1104,7 +1104,7 @@ 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 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, hap_cov_t *i_cov, bub_label_t* b_mask_t, uint32_t collect_p_trans);
+kvec_asg_arc_t_warp* new_rtg_edges, hap_cov_t **i_cov, bub_label_t* b_mask_t, uint32_t collect_p_trans);
 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, 

Purge_Dups.cpp

@@ -114,14 +114,20 @@ typedef struct {
     uint32_t baseBeg, baseEnd;
     uint32_t nodeBeg, nodeEnd;
     uint32_t h_lev_idx;
-    uint32_t b_ug_id;
+    uint32_t b_ug_id, c_ug_id;
 }p_node_t;
 
+typedef struct {
+    uint32_t beg;
+    uint32_t occ;
+}p_g_in_t;
+
 typedef struct {
     ma_ug_t *ug;
     kvec_t(p_node_t) pg_het_node;
     asg_t *pg_het;
     asg_t *pg_h_lev;
+    kvec_t(p_g_in_t) pg_h_lev_idx;
 }p_g_t;
 
 void print_peak_line(int c, int x, int exceed, int64_t cnt)
@@ -3269,6 +3275,151 @@ ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex)
     }
 }
 
+void get_p_nodes(p_g_t *pg, p_node_t **x, uint32_t *x_occ, uint32_t id)
+{
+    if(x) (*x) = pg->pg_het_node.a + pg->pg_h_lev_idx.a[id].beg;
+    if(x_occ) (*x_occ) = pg->pg_h_lev_idx.a[id].occ;
+}
+
+void normalize_hap_overlaps_advance_by_p_g_t(hap_overlaps_list* all_ovlp, hap_overlaps_list* back_all_ovlp,
+ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, p_g_t *pg, hap_cov_t *cov, 
+double filter_rate)
+{
+    hap_overlaps *x = NULL, *y = NULL;
+    uint32_t v, i, uId, qn, tn;
+    uint32_t types[4];
+    
+
+    types[X2Y] = Y2X; types[Y2X] = X2Y; types[XCY] = YCX; types[YCX] = XCY;
+    int index;
+    uint32_t k, qs, qe, ts, te, occ, as, ae, ovlp, hetLen, homLen;
+    p_node_t *a = NULL;
+
+    /*******************************for debug************************************/
+    // for (v = 0; v < cov->t_ch->u_num; v++)
+    // {
+    //     fprintf(stderr, "utg%.6ul-is_het=%u\n", v+1, cov->t_ch->is_het[v]);
+    // }
+    /*******************************for debug************************************/
+
+
+
+
+
+
+    for (v = 0; v < all_ovlp->num; v++)
+    {
+        uId = v;
+        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
+        {
+            /*****************qn*****************/
+            qn = all_ovlp->x[uId].a.a[i].xUid;
+            qs = all_ovlp->x[uId].a.a[i].x_beg_pos;
+            qe = all_ovlp->x[uId].a.a[i].x_end_pos - 1;
+            get_p_nodes(pg, &a, &occ, qn);
+            for (k = 0, hetLen = 0, homLen = 0; k < occ; k++)
+            {
+                as = a[k].baseBeg;
+                ae = a[k].baseEnd;
+                ovlp = ((MIN(qe, ae) >= MAX(qs, as))? MIN(qe, ae) - MAX(qs, as) + 1 : 0);
+                if(homLen + hetLen > 0 && ovlp == 0) break;
+                if(ovlp == 0) continue;
+                if(cov->t_ch->is_het[a[k].b_ug_id>>1] == 0)
+                {
+                    homLen += ovlp;
+                } 
+                else
+                {
+                    hetLen += ovlp;
+                }                
+            }
+
+            if(hetLen <= ((hetLen + homLen) * filter_rate))
+            {
+                all_ovlp->x[uId].a.a[i].status = DELETE;
+                continue;
+            }
+            /*****************qn*****************/
+
+
+            /*****************tn*****************/
+            tn = all_ovlp->x[uId].a.a[i].yUid;
+            ts = all_ovlp->x[uId].a.a[i].y_beg_pos;
+            te = all_ovlp->x[uId].a.a[i].y_end_pos - 1;
+            get_p_nodes(pg, &a, &occ, tn);
+            for (k = 0, hetLen = 0, homLen = 0; k < occ; k++)
+            {
+                as = a[k].baseBeg;
+                ae = a[k].baseEnd;
+                ovlp = ((MIN(te, ae) >= MAX(ts, as))? MIN(te, ae) - MAX(ts, as) + 1 : 0);
+                if(homLen + hetLen > 0 && ovlp == 0) break;
+                if(ovlp == 0) continue;
+                if(cov->t_ch->is_het[a[k].b_ug_id>>1] == 0)
+                {
+                    homLen += ovlp;
+                } 
+                else
+                {
+                    hetLen += ovlp;
+                }                
+            }
+
+            if(hetLen <= ((hetLen + homLen) * filter_rate))
+            {
+                all_ovlp->x[uId].a.a[i].status = DELETE;
+                continue;
+            }
+            /*****************tn*****************/
+        }
+    }
+
+    for (v = 0; v < all_ovlp->num; v++)
+    {
+        uId = v;
+        k = 0;
+        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
+        {
+            if(all_ovlp->x[uId].a.a[i].status == DELETE) continue;
+            all_ovlp->x[uId].a.a[k] = all_ovlp->x[uId].a.a[i];
+            k++;
+        }
+        all_ovlp->x[uId].a.n = k;
+    }
+
+
+    for (v = 0; v < all_ovlp->num; v++)
+    {
+        uId = v;
+        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
+        {
+            qn = all_ovlp->x[uId].a.a[i].xUid;
+            tn = all_ovlp->x[uId].a.a[i].yUid;
+            x = &(all_ovlp->x[uId].a.a[i]);
+            index = get_specific_hap_overlap(&(all_ovlp->x[tn]), tn, qn);
+            if(index != -1)
+            {
+                y = &(all_ovlp->x[tn].a.a[index]);
+                if(x->rev == y->rev && types[x->type]==y->type) continue;
+                if(x->score >= y->score)
+                {
+                    kv_push(hap_overlaps, back_all_ovlp->x[tn].a, (*y));
+                    set_reverse_hap_overlap(y, x, types);
+                }
+                else
+                {
+                    kv_push(hap_overlaps, back_all_ovlp->x[qn].a, (*x));
+                    set_reverse_hap_overlap(x, y, types);
+                }
+            }
+            else
+            {
+                kv_pushp(hap_overlaps, all_ovlp->x[tn].a, &y);
+                set_reverse_hap_overlap(y, x, types);
+            }
+        }
+    }
+}
+
 void filter_hap_overlaps_by_length(hap_overlaps_list* all_ovlp, uint32_t minLen)
 {
     if(minLen == 0) return;
@@ -4863,7 +5014,7 @@ int max_hang, int min_ovlp, float drop_ratio)
 
 }
 
-void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, 
+void purge_dups_back(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, float drop_ratio, uint32_t just_contain, 
 uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
@@ -5036,18 +5187,82 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
     destory_hap_alignment_struct_pip(&hap_buf);
 }
 
+void debug_p_g_t(p_g_t* pg, asg_t *read_g)
+{
+    fprintf(stderr, "----------[M::%s]----------\n", __func__);
+    uint32_t i, offset, v, sid, eid, spos, epos, puid, occ;
+    p_node_t *t = NULL;
+    ma_utg_t *u = NULL;
+    p_node_t *a = NULL;
+
+    for (v = 0; v < pg->ug->u.n; v++)
+    {
+        ///fprintf(stderr, "\nu->n: %u, uid: %u\n", (uint32_t)(pg->ug->u.a[v].n), v);
+        get_p_nodes(pg, &a, &occ, v);
+        for (i = 0; i < occ; i++)
+        {
+            if(a[i].c_ug_id != v) fprintf(stderr, "sbsbsbsbsb\n");
+            ///fprintf(stderr, "sid: %u, eid: %u\n", a[i].nodeBeg, a[i].nodeEnd);
+        }
+    }
+    
+
+    for (v = 0, puid = (uint32_t)-1; v < pg->pg_het_node.n; v++)
+    {
+        t = &(pg->pg_het_node.a[v]);
+        sid = t->nodeBeg;
+        eid = t->nodeEnd;
+        spos = t->baseBeg;
+        epos = t->baseEnd;
+        // fprintf(stderr, "sid: %u, eid: %u, spos: %u, epos: %u, t->b_ug_id: %u\n", 
+        //                                                 sid, eid, spos, epos, (uint32_t)t->b_ug_id);
+        // fprintf(stderr, "pg->pg_het_node.n: %u\n", (uint32_t)pg->pg_het_node.n);
+        if(puid == t->b_ug_id) fprintf(stderr, "ERROR-(-1)\n");
+        if(t->b_ug_id == (uint32_t)-1) fprintf(stderr, "ERROR-(-2)\n");
+        puid = t->b_ug_id;
+
+        u = &(pg->ug->u.a[t->c_ug_id]);
+        ///fprintf(stderr, "u->n: %u, sid: %u, eid: %u\n", (uint32_t)u->n, sid, eid);
+        for (i = offset = 0; i < u->n; i++)
+        {
+            if(i == sid)
+            {
+                if(spos != offset)
+                {
+                    fprintf(stderr, "ERROR-1\n");
+                }
+            }
+
+            if(i == eid)
+            {
+                if(epos != (offset+read_g->seq[u->a[i]>>33].len - 1))
+                {
+                    fprintf(stderr, "ERROR-2, real end: %u\n", 
+                                    (uint32_t)(offset+read_g->seq[u->a[i]>>33].len - 1));
+                }
+            }
+            offset += (uint32_t)u->a[i];
+        }
+    }
+    
+    
+}
+
 
 p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
 {
-    uint32_t v, uId, k_uId, l_uid, k, l, offset, l_pos;
+    uint32_t v, uId, k_uId, l_uid, k, l, offset, l_pos, g_beg_idx, occ, ovlp, tLen, zLen;
     p_g_t *pg = NULL; CALLOC(pg, 1);
     pg->ug = ug;
     asg_t* nsg = pg->ug->g;
     ma_utg_t *u = NULL;
     p_node_t *t = NULL, *z = NULL;
+    asg_arc_t *e = NULL;
+    p_g_in_t *x = NULL;
     pg->pg_het = asg_init();
     pg->pg_h_lev = asg_init();
     kv_init(pg->pg_het_node);
+    kv_init(pg->pg_h_lev_idx);
 
     for (v = 0; v < nsg->n_seq; v++)
     {
@@ -5069,6 +5284,8 @@ p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
         if(nsg->seq[uId].del || nsg->seq[uId].c == ALTER_LABLE) continue;
 
         u = &(ug->u.a[uId]);
+        g_beg_idx = pg->pg_het_node.n;
+        ///fprintf(stderr, "\n+v: %u, pg->pg_het_node.n: %u\n", v, (uint32_t)pg->pg_het_node.n);
         for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k) 
         {   
             l_uid = k_uId = (uint32_t)-1;
@@ -5077,15 +5294,17 @@ p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
 
             if (k == u->n || k_uId != l_uid)
             {
+                ///fprintf(stderr, "l: %u, k: %u, u->n: %u, l_uid: %u\n", l, k, (uint32_t)u->n, l_uid);
                 if(l_uid != (uint32_t)-1)
                 {
                     kv_pushp(p_node_t, pg->pg_het_node, &t);
+                    t->c_ug_id = uId;
                     t->b_ug_id = l_uid;
                     t->baseBeg = l_pos;
                     t->baseEnd = offset + read_g->seq[u->a[k-1]>>33].len - 1;
                     t->nodeBeg = l;
                     t->nodeEnd = k - 1;
-                    if(pg->pg_het_node.n > 1)
+                    if(pg->pg_het_node.n > 1 && pg->pg_het_node.n >= (g_beg_idx + 2))
                     {
                         z = &(pg->pg_het_node.a[pg->pg_het_node.n - 2]);
                         if(t->b_ug_id == z->b_ug_id)
@@ -5093,9 +5312,10 @@ p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
                             z->baseEnd = t->baseEnd;
                             z->nodeEnd = t->nodeEnd;
                             t = z;
+                            pg->pg_het_node.n--;
                         }
-                        pg->pg_het_node.n--;
                     }
+                    ///if(t->b_ug_id == (uint32_t)-1) fprintf(stderr, "xxxx\n");
                     asg_seq_set(pg->pg_het, pg->pg_het_node.n-1, t->baseEnd+1-t->baseBeg, 0);
                 }
                 l = k;
@@ -5103,20 +5323,59 @@ p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g)
             }
             offset += (uint32_t)u->a[k-1];
         }
+
+        occ = pg->pg_het_node.n - g_beg_idx;
+        kv_pushp(p_g_in_t, pg->pg_h_lev_idx, &x);
+        x->beg = g_beg_idx; x->occ = occ;
+        ///fprintf(stderr, "-v: %u, pg->pg_het_node.n: %u\n", v, (uint32_t)pg->pg_het_node.n);
+        
+        if(occ > 1)
+        {
+            for (k = g_beg_idx; (k + 1) < pg->pg_het_node.n; ++k) 
+            {
+                t = &(pg->pg_het_node.a[k]); tLen = t->baseEnd + 1 - t->baseBeg;
+                z = &(pg->pg_het_node.a[k+1]); zLen = z->baseEnd + 1 - z->baseBeg;
+
+                ovlp = ((MIN(t->baseEnd, z->baseEnd) >= MAX(t->baseBeg, z->baseBeg))? 
+                                MIN(t->baseEnd, z->baseEnd) - MAX(t->baseBeg, z->baseBeg) + 1 : 0);
+
+                e = asg_arc_pushp(pg->pg_het);
+                e->ol = ovlp;
+                e->ul = (k<<1); e->ul <<= 32; e->ul += (tLen - ovlp);
+                e->v = ((k+1)<<1); e->del = 0; e->el = e->no_l_indel = e->strong = 1;
+
+                e = asg_arc_pushp(pg->pg_het);
+                e->ol = ovlp;
+                e->ul = ((k+1)<<1)+1; e->ul <<= 32; e->ul += (zLen - ovlp);
+                e->v = (k<<1)+1; e->del = 0; e->el = e->no_l_indel = e->strong = 1;
+            }
+        }
     }
-
+    asg_cleanup(pg->pg_het);
-
+    ///debug_p_g_t(pg, read_g);
 
     return pg;
 }
 
-void purge_dups_advance(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, 
+void destory_p_g_t(p_g_t **pg)
+{
+    if(pg && (*pg))
+    {
+        kv_destroy((*pg)->pg_het_node);
+        kv_destroy((*pg)->pg_h_lev_idx);
+        asg_destroy((*pg)->pg_het);
+        asg_destroy((*pg)->pg_h_lev);
+        free((*pg)); 
+        (*pg) = NULL;
+    }
+}
+
+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, float drop_ratio, uint32_t just_contain, 
 uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
 {
-    asg_t *purge_g = NULL;
+    p_g_t *pg = NULL;
-    purge_g = asg_init();
     asg_t* nsg = ug->g;
     uint32_t v, rId, uId, i, offset;
     ma_utg_t* reads = NULL;
@@ -5149,12 +5408,6 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
     for (v = 0; v < nsg->n_seq; v++)
     {
         uId = v;
-        if(nsg->seq[uId].del || nsg->seq[uId].c == ALTER_LABLE)
-        {
-            asg_seq_set(purge_g, uId, 0, 1);
-            purge_g->seq[uId].c = ALTER_LABLE;
-            continue;
-        } 
         reads = &(ug->u.a[uId]);
         for (i = 0, offset = 0; i < reads->n; i++)
         {
@@ -5167,9 +5420,6 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
 
             offset += (uint32_t)reads->a[i];
         }
-
-        asg_seq_set(purge_g, uId, offset, 0);
-        purge_g->seq[uId].c = PRIMARY_LABLE;
     }
 
 
@@ -5196,30 +5446,34 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
     if(just_coverage) goto end_coverage;
 
     kt_for(asm_opt.thread_num, hap_alignment_advance_worker, &hap_buf, nsg->n_seq);
-
+    
     ///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp);
     filter_hap_overlaps_by_length(&all_ovlp, purege_minLen);
 
+    pg = init_p_g_t(ug, cov, read_g);
     ///normalize_hap_overlaps(&all_ovlp, &back_all_ovlp);
-    normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
+
+    normalize_hap_overlaps_advance_by_p_g_t(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, 
+    ruIndex, pg, cov, 0.8);
+
+    ///normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
     ///debug_hap_overlaps(&all_ovlp, &back_all_ovlp);
 
-    remove_contained_haplotig(&all_ovlp, ug, nsg, purge_g, cov);
+    remove_contained_haplotig(&all_ovlp, ug, nsg, pg->pg_h_lev, cov);
     
     if(just_contain == 0)
     {
         for (v = 0; v < all_ovlp.num; v++)
         {
             uId = v;
-            if(purge_g->seq[uId].del || purge_g->seq[uId].c == ALTER_LABLE) continue;
+            if(pg->pg_h_lev->seq[uId].del || pg->pg_h_lev->seq[uId].c == ALTER_LABLE) continue;
             for (i = 0; i < all_ovlp.x[uId].a.n; i++)
             {
                 if(all_ovlp.x[uId].a.a[i].status == DELETE) continue;
-                ///if(all_ovlp.x[uId].a.a[i].type == )
+                if(pg->pg_h_lev->seq[all_ovlp.x[uId].a.a[i].xUid].c == ALTER_LABLE||
-                if(purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].c == ALTER_LABLE||
+                   pg->pg_h_lev->seq[all_ovlp.x[uId].a.a[i].xUid].del||
-                   purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].del||
+                   pg->pg_h_lev->seq[all_ovlp.x[uId].a.a[i].yUid].c == ALTER_LABLE||
-                   purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].c == ALTER_LABLE||
+                   pg->pg_h_lev->seq[all_ovlp.x[uId].a.a[i].yUid].del)
-                   purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].del)
                 {
                     continue;
                 }
@@ -5237,20 +5491,20 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
                 // }
                 
                 if(r < 0) continue;
-                p = asg_arc_pushp(purge_g);
+                p = asg_arc_pushp(pg->pg_h_lev);
                 *p = t;
             }
         }
 
-        asg_cleanup(purge_g);
+        asg_cleanup(pg->pg_h_lev);
-		asg_symm(purge_g);
+		asg_symm(pg->pg_h_lev);
         ///may need to do transitive reduction
-        clean_purge_graph(purge_g, drop_ratio, 1);
+        clean_purge_graph(pg->pg_h_lev, drop_ratio, 1);
 
         // if(debug_enable) print_purge_gfa(ug, purge_g);
         // if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp);
 
-        link_unitigs(purge_g, ug, &all_ovlp, ruIndex, reverse_sources, coverage_cut, read_g, position_index, 
+        link_unitigs(pg->pg_h_lev, ug, &all_ovlp, ruIndex, reverse_sources, coverage_cut, read_g, position_index, 
         &(hap_buf.buf[0].u_buffer), &(hap_buf.buf[0].u_buffer_tailIndex), &(hap_buf.buf[0].u_buffer_prevIndex),
         max_hang, min_ovlp, edge, hap_buf.buf[0].visit, cov);
     }
@@ -5258,28 +5512,25 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
     for (v = 0; v < all_ovlp.num; v++)
     {
         uId = v;
-        if(purge_g->seq[uId].c == ALTER_LABLE)
+        if(pg->pg_h_lev->seq[uId].c == ALTER_LABLE)
         {
             ug->g->seq[uId].c = ALTER_LABLE;
         }
     }
 
     end_coverage:
-
     uint32_t is_Unitig;
     for (v = 0; v < ruIndex->len; v++)
     {
         get_R_to_U(ruIndex, v, &uId, &is_Unitig);
         if(is_Unitig == 1) ruIndex->index[v] = (uint32_t)-1;
     }
-
     asg_cleanup(nsg);
     destory_hap_overlaps_list(&all_ovlp);
     destory_hap_overlaps_list(&back_all_ovlp);
-    asg_destroy(purge_g);
     if(cov) memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);
     else free(position_index);
     
     destory_hap_alignment_struct_pip(&hap_buf);
+    destory_p_g_t(&pg);
 }
-

hic.cpp

@@ -2678,17 +2678,17 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *het_flag)
 
         for (i = 0; i < ug->g->n_seq; i++)
         {
-            if((bub->index[i]>>2) == 0)
+            if((bub->index[i]>>2) == 0) 
             {
-                bub->index[i] = (uint32_t)-1; 
+                bub->index[i] = (uint32_t)-1; ///hom
             }
             else
             {
-                if((bub->index[i]>>2) == 1)
+                if((bub->index[i]>>2) == 1) 
                 {
                     bub->index[i] = P_het(*bub); ///potential het
                 }
-                else
+                else 
                 {
                     bub->index[i] = M_het(*bub); ///must het
                 }
@@ -2705,11 +2705,11 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *het_flag)
                 n_occ += ug->u.a[a[v]>>1].n;
             }
 
-            if((pathLen*2) >= ug->g->seq[beg>>1].len && (pathLen*2) >= ug->g->seq[sink>>1].len)
+            // if((pathLen*2) >= ug->g->seq[beg>>1].len && (pathLen*2) >= ug->g->seq[sink>>1].len)
-            {
+            // {
-                bub->index[(beg>>1)] = (uint32_t)-1;
+            //     bub->index[(beg>>1)] = (uint32_t)-1;
-                bub->index[(sink>>1)] = (uint32_t)-1;
+            //     bub->index[(sink>>1)] = (uint32_t)-1;
-            }
+            // }
 
             if(n_occ > 3)
             {
@@ -2749,7 +2749,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *het_flag)
             if(bub->index[i] == M_het(*bub)) bub->index[i] = P_het(*bub);
             if(bub->index[i] > P_het(*bub))
             {
-                if(het_flag && het_flag[i] == 1)
+                if(het_flag && het_flag[i])
                 {
                     bub->index[i] = P_het(*bub);
                 }
@@ -2769,6 +2769,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *het_flag)
     }
     bub->b_g = NULL;
     bub->b_ug = NULL;
+
     build_bub_graph(ug, bub);
 }