From 24a19d797635074c31ada39e5f62cdbedba09d3c Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Fri, 26 Mar 2021 18:07:30 -0400 Subject: [PATCH] more accurate purging --- Overlaps.cpp | 234 +++++++++++++++++++++-------------- Overlaps.h | 2 +- Purge_Dups.cpp | 327 +++++++++++++++++++++++++++++++++++++++++++------ hic.cpp | 21 ++-- 4 files changed, 446 insertions(+), 138 deletions(-) diff --git a/Overlaps.cpp b/Overlaps.cpp index 0be3fcb..d5f9ee8 100644 --- a/Overlaps.cpp +++ b/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; - asm_opt.hom_global_coverage = -1; + int is_set = ((asm_opt.hom_global_coverage != -1)? 1 : 0); - 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; - asm_opt.hom_global_coverage = tmp_cov; - ///fprintf(stderr, "dip_thre_max: %lu\n", dip_thre_max); - + dip_thre_max = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE); + } + 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(©_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); - else cov = i_cov; + ma_ug_t *cov_ug = NULL; + 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(collect_p_trans) + if(asm_opt.purge_level_primary > 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, - 0, 0, cov, collect_p_trans); + if(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(©_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; diff --git a/Overlaps.h b/Overlaps.h index ddd621e..346d80a 100644 --- a/Overlaps.h +++ b/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, diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index c6cdc1a..35acc48 100644 --- a/Purge_Dups.cpp +++ b/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; - purge_g = asg_init(); + p_g_t *pg = NULL; 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(purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].c == ALTER_LABLE|| - purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].del|| - purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].c == ALTER_LABLE|| - purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].del) + if(pg->pg_h_lev->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|| + pg->pg_h_lev->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) { 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_symm(purge_g); + asg_cleanup(pg->pg_h_lev); + 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); } - diff --git a/hic.cpp b/hic.cpp index 1583178..cd6dc47 100644 --- a/hic.cpp +++ b/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) - { - bub->index[(beg>>1)] = (uint32_t)-1; - bub->index[(sink>>1)] = (uint32_t)-1; - } + // 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[(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); }