diff --git a/Overlaps.cpp b/Overlaps.cpp index 6a3a2bd..0456f52 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -11304,9 +11304,13 @@ void collect_trans_cov(const char* cmd, buf_t* pri, buf_t* aux, ma_ug_t *ug, asg chain_trans_ovlp(cov, ug, read_sg, pri, len_aux, &thre_pri); if(thre_pri > 0) { - // fprintf(stderr, "\n%s, thre_pri: %u\n", cmd, thre_pri); + /*******************************for debug************************************/ + // fprintf(stderr, "\n%s, thre_pri: %u, len_aux: %lu\n", cmd, thre_pri, len_aux); // print_buf_t(ug, pri, "pri"); // print_buf_t(ug, aux, "aux"); + /*******************************for debug************************************/ + + for (i = uCov = 0, p_uId = (uint32_t)-1; i < aux->b.n; i++) { @@ -12111,25 +12115,32 @@ bub_label_t* b_mask_t) 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; + fprintf(stderr, "sb0sb\n"); 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); + fprintf(stderr, "sb1sb\n"); ma_ug_destroy(copy_ug); asg_destroy(copy_sg); - + fprintf(stderr, "sb2sb\n"); 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); + fprintf(stderr, "sb3sb\n"); - - + fprintf(stderr, "sb4sb\n"); new_rtg_edges.a.n = 0; ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp); ///classify_untigs(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, &new_rtg_edges, max_hang, min_ovlp); + fprintf(stderr, "sb5sb\n"); + hic_analysis(ug, sg, cov); + + fprintf(stderr, "sb6sb\n"); + destory_hap_cov_t(&cov); ma_ug_destroy(ug); kv_destroy(new_rtg_edges.a); @@ -14107,7 +14118,7 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov) #define T_ROUND 2 asg_t *g = ug->g; int round = T_ROUND; - + redo: ///print_graph_statistic(g, "beg"); ///print_debug_gfa(read_g, ug, coverage_cut, "debug_trans_ovlp_hg002", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); @@ -16899,8 +16910,9 @@ uint64_t get_s_bub_pop_max_dist_advance(asg_t *g, buf_s_t *b) int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov, uint32_t is_update_chain) { asg_t *g = ug->g; - uint32_t v, n_vtx = g->n_seq * 2; + uint32_t v, n_vtx = g->n_seq * 2, n_arc, nv, i; uint64_t n_pop = 0, max_dist; + asg_arc_t *av = NULL; buf_t b; if (!g->is_symm) asg_symm(g); memset(&b, 0, sizeof(buf_t)); @@ -16908,27 +16920,57 @@ int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t posi b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); if(i_max_dist) max_dist = (*i_max_dist); else max_dist = get_bub_pop_max_dist_advance(g, &b); + uint8_t* bs_flag = NULL; CALLOC(bs_flag, n_vtx); + if(max_dist > 0) { - //traverse all node with two directions - for (v = 0; v < n_vtx; ++v) { - uint32_t i, n_arc = 0, nv = asg_arc_n(g, v); - asg_arc_t *av = asg_arc_a(g, v); + for (v = 0; v < n_vtx; ++v) + { + if(bs_flag[v] != 0) continue; + nv = asg_arc_n(g, v); + av = asg_arc_a(g, v); ///some node could be deleted if (nv < 2 || g->seq[v>>1].del || g->seq[v>>1].c == ALTER_LABLE) continue; ///some edges could be deleted - for (i = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs + for (i = n_arc = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs + if (!av[i].del) ++n_arc; + if (n_arc < 2) continue; + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, max_dist, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0)) + { + //beg is v, end is b.S.a[0] + //note b.b include end, does not include beg + for (i = 0; i < b.b.n; i++) + { + if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue; + bs_flag[b.b.a[i]] = bs_flag[b.b.a[i]^1] = 1; + } + bs_flag[v] = 2; bs_flag[b.S.a[0]^1] = 3; + } + } + + //traverse all node with two directions + for (v = 0; v < n_vtx; ++v) { + if(bs_flag[v] !=2) continue; + nv = asg_arc_n(g, v); + av = asg_arc_a(g, v); + ///some node could be deleted + if (nv < 2 || g->seq[v>>1].del || g->seq[v>>1].c == ALTER_LABLE) continue; + ///some edges could be deleted + for (i = n_arc = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs if (!av[i].del) ++n_arc; if (n_arc > 1) n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, cov, is_update_chain); } - free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); - if (n_pop) asg_cleanup(g); + if(VERBOSE >= 1) { fprintf(stderr, "[M::%s] popped %lu bubbles\n", __func__, (unsigned long)n_pop); } } + + free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); + if (n_pop) asg_cleanup(g); + free(bs_flag); return n_pop; } @@ -22675,6 +22717,7 @@ uint32_t collect_p_trans) EvaluateLen((*ug)->u, v) = (*ug)->u.a[v].n; } + clean_primary_untig_graph(*ug, read_g, sources, reverse_sources, coverage_cut, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, chimeric_rate, 0, 0, drop_ratio, cov); delete_useless_nodes(ug); @@ -22684,6 +22727,7 @@ uint32_t collect_p_trans) if(asm_opt.purge_level_primary > 0) { + ///print_debug_gfa(read_g, *ug, coverage_cut, "debug_purge", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); just_contain = 0; if(asm_opt.purge_level_primary == 1) just_contain = 1; purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, @@ -22701,7 +22745,7 @@ uint32_t collect_p_trans) rescue_contained_reads_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang, min_ovlp, 10, 0, 1, NULL, NULL, b_mask_t); renew_utg(ug, read_g, new_rtg_edges); - + // if(asm_opt.purge_level_primary > 0) // { // just_contain = 0; diff --git a/Overlaps.h b/Overlaps.h index c4e0608..4381018 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1133,7 +1133,7 @@ inline uint32_t get_origin_uid(uint32_t v, trans_chain* t_ch) if(t_ch->u_idx[v>>1] == (uint32_t)-1) return (uint32_t)-1; return ((t_ch->u_idx[v>>1]>>1)<<1) + ((t_ch->u_idx[v>>1]^v)&1); } - +void get_chain_trans(trans_chain* t_ch, uint32_t id, uint32_t** x, uint32_t* x_occ, uint32_t** y, uint32_t* y_occ); #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index 0a1bd02..85c6538 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -1718,7 +1718,8 @@ void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads uint64_t tmp; asg_arc_t_offset t_offset; u_buffer->a.n = 0; - (*xEnd) = (uint32_t)-1; + ///(*xEnd) = (uint32_t)-1; + (*xEnd) = 0; uint32_t u_i, r_i, k, j, m, len, p_v, *a = xReads->b.a, uid, ori, l, aOcc, nv, xOcc = (uint32_t)-1; ma_utg_t* u = NULL; asg_arc_t *av = NULL; @@ -2389,6 +2390,7 @@ kvec_t_i32_warp* prevIndex, long long* r_x_pos_beg, long long* r_x_pos_end, long hap_can->index_beg = xLeftTotal; hap_can->score = get_chain_score(xReads, read_g, u_buffer, tailIndex, prevIndex, reverse_sources, (*r_x_pos_beg), (*r_x_pos_end)); + if(hap_can->score <= 0) return (uint32_t)-1; return hap_can->index_end; } @@ -3260,13 +3262,6 @@ double filter_rate) 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++) { @@ -3289,7 +3284,7 @@ double filter_rate) { homLen += ovlp; } - else if(asm_opt.polyploidy <= 2 && (a[k].h_status&S_HET)) + else if(asm_opt.polyploidy <= 2 && (a[k].h_status&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].h_status&S_HET)) { homLen += ovlp; } @@ -3323,7 +3318,7 @@ double filter_rate) { homLen += ovlp; } - else if(asm_opt.polyploidy <= 2 && (a[k].h_status&S_HET)) + else if(asm_opt.polyploidy <= 2 && (a[k].h_status&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].h_status&S_HET)) { homLen += ovlp; } @@ -3382,24 +3377,11 @@ double filter_rate) } else { - x->status = DELETE; + kv_pushp(hap_overlaps, all_ovlp->x[tn].a, &y); + set_reverse_hap_overlap(y, x, types); } } } - - - 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; - } } void filter_hap_overlaps_by_length(hap_overlaps_list* all_ovlp, uint32_t minLen) @@ -4388,6 +4370,166 @@ void collect_reverse_unitigs_purge(buf_t* b_0, hc_links* link, ma_ug_t *ug, hap_ } **/ +void print_het_ovlp(p_g_t *pg, ma_ug_t *ug, hap_overlaps_list* ha, double filter_rate) +{ + uint32_t v, i, k, n_vtx = pg->pg_h_lev->n_seq * 2, nv, qn, qs, qe, tn, ts, te, as, ae, occ, ovlp, hetLen, homLen; + asg_arc_t *av = NULL; + hap_overlaps *x = NULL; + p_node_t *a = NULL; + int index; + for (v = 0; v < n_vtx; v++) + { + av = asg_arc_a(pg->pg_h_lev, v); + nv = asg_arc_n(pg->pg_h_lev, v); + for (i = 0; i < nv; i++) + { + if(av[i].del) continue; + index = get_specific_hap_overlap(&(ha->x[av[i].ul>>33]), av[i].ul>>33, av[i].v>>1); + if(index == -1) fprintf(stderr, "ERROR\n"); + x = &(ha->x[av[i].ul>>33].a.a[index]); + + qn = x->xUid; + qs = x->x_beg_pos; + qe = x->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(a[k].h_status == N_HET) + { + homLen += ovlp; + } + else if(asm_opt.polyploidy <= 2 && (a[k].h_status&P_HET)) + { + homLen += ovlp; + } + else + { + hetLen += ovlp; + } + } + + if(hetLen <= ((hetLen + homLen) * filter_rate)) + { + ///all_ovlp->x[uId].a.a[i].status = DELETE; + fprintf(stderr, "********XY********\n"); + print_hap_paf(ug, x); + } + + + tn = x->yUid; + ts = x->y_beg_pos; + te = x->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(a[k].h_status == N_HET) + { + homLen += ovlp; + } + else if(asm_opt.polyploidy <= 2 && (a[k].h_status&P_HET)) + { + homLen += ovlp; + } + else + { + hetLen += ovlp; + } + } + + if(hetLen <= ((hetLen + homLen) * filter_rate)) + { + ///all_ovlp->x[uId].a.a[i].status = DELETE; + fprintf(stderr, "********YX********\n"); + print_hap_paf(ug, x); + } + } + } + + for (v = 0; v < ha->num; v++) + { + for (i = 0; i < ha->x[v].a.n; i++) + { + if(ha->x[v].a.a[i].status == DELETE) + { + x = &(ha->x[v].a.a[i]); + + qn = x->xUid; + qs = x->x_beg_pos; + qe = x->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(a[k].h_status == N_HET) + { + homLen += ovlp; + } + else if(asm_opt.polyploidy <= 2 && (a[k].h_status&P_HET)) + { + homLen += ovlp; + } + else + { + hetLen += ovlp; + } + } + + if(hetLen <= ((hetLen + homLen) * filter_rate)) + { + fprintf(stderr, "********C(X)********hetLen-%u, homLen-%u\n", hetLen, homLen); + print_hap_paf(ug, x); + } + + + tn = x->yUid; + ts = x->y_beg_pos; + te = x->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(a[k].h_status == N_HET) + { + homLen += ovlp; + } + else if(asm_opt.polyploidy <= 2 && (a[k].h_status&P_HET)) + { + homLen += ovlp; + } + else + { + hetLen += ovlp; + } + } + + if(hetLen <= ((hetLen + homLen) * filter_rate)) + { + fprintf(stderr, "********C(Y)********hetLen-%u, homLen-%u\n", hetLen, homLen); + print_hap_paf(ug, x); + } + } + } + } +} + void link_unitigs(asg_t *purge_g, ma_ug_t *ug, hap_overlaps_list* all_ovlp, R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, @@ -5392,6 +5534,56 @@ int max_hang, int min_ovlp, float drop_ratio, p_g_t *pg) } +void collect_purge_trans_cov(ma_ug_t *ug, hap_overlaps_list* ha, trans_chain* t_ch) +{ + uint32_t v, i, k, e, s, o, c_uId, p_uId, x_occ, y_occ; + ma_utg_t *q = NULL; + hap_overlaps *x = NULL; + for (v = 0; v < ha->num; v++) + { + for (i = 0; i < ha->x[v].a.n; i++) + { + if(x->yUid < x->xUid) continue; + x = &(ha->x[v].a.a[i]); + + q = &(ug->u.a[x->xUid]); s = x->x_beg_id; e = x->x_end_id; o = 0; + for (k = s, p_uId = (uint32_t)-1; k < e; k++) + { + c_uId = get_origin_uid((o == 1?((q->a[e-k-1]^(uint64_t)(0x100000000))>>32):(q->a[k]>>32)), t_ch); + if(c_uId == (uint32_t)-1 || p_uId == c_uId) continue; + p_uId = c_uId; + kv_push(uint32_t, t_ch->uIDs, c_uId); + } + kv_push(uint32_t, t_ch->iDXs, t_ch->uIDs.n); + + + q = &(ug->u.a[x->yUid]); s = x->y_beg_id; e = x->y_end_id; o = x->rev; + for (k = s, p_uId = (uint32_t)-1; k < e; k++) + { + c_uId = get_origin_uid((o == 1?((q->a[e-k-1]^(uint64_t)(0x100000000))>>32):(q->a[k]>>32)), t_ch); + if(c_uId == (uint32_t)-1 || p_uId == c_uId) continue; + p_uId = c_uId; + kv_push(uint32_t, t_ch->uIDs, c_uId); + } + kv_push(uint32_t, t_ch->iDXs, t_ch->uIDs.n); + + + x_occ = y_occ = 0; + get_chain_trans(t_ch, t_ch->chain_num, NULL, &x_occ, NULL, &y_occ); + if(x_occ == 0 || y_occ == 0) + { + t_ch->uIDs.n -= (x_occ + y_occ); + t_ch->iDXs.n -= 2; + } + else + { + t_ch->chain_num++; + t_ch->l1_chain++; + } + } + } +} + 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, @@ -5472,12 +5664,17 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans) ///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp); filter_hap_overlaps_by_length(&all_ovlp, purege_minLen); + normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex); + if(asm_opt.polyploidy <= 2) pt_solve(&all_ovlp, cov->t_ch, ug, read_g, 0.8, R_INF.trio_flag); - if(collect_p_trans) goto end_coverage; - + if(collect_p_trans) + { + collect_purge_trans_cov(ug, &all_ovlp, cov->t_ch); + goto end_coverage; + } + pg = init_p_g_t(ug, cov, read_g); - ///normalize_hap_overlaps(&all_ovlp, &back_all_ovlp); - + normalize_hap_overlaps_advance_by_p_g_t(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex, pg, cov, 0.8); @@ -5509,6 +5706,12 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans) r = get_hap_arch(&(all_ovlp.x[uId].a.a[i]), ug->u.a[all_ovlp.x[uId].a.a[i].xUid].len, ug->u.a[all_ovlp.x[uId].a.a[i].yUid].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + // if(all_ovlp.x[uId].a.a[i].xUid == 118 && all_ovlp.x[uId].a.a[i].yUid == 82) + // { + // fprintf(stderr, "r: %d\n", r); + // print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i])); + // } + if(r < 0) continue; p = asg_arc_pushp(pg->pg_h_lev); *p = t; @@ -5522,6 +5725,9 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans) // if(debug_enable) print_purge_gfa(ug, purge_g); // if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp); + /*******************************for debug************************************/ + // print_het_ovlp(pg, ug, &all_ovlp, 0.8); + /*******************************for debug************************************/ 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), diff --git a/Purge_Dups.h b/Purge_Dups.h index 18763a5..a191468 100644 --- a/Purge_Dups.h +++ b/Purge_Dups.h @@ -65,5 +65,6 @@ void destory_hap_cov_t(hap_cov_t **x); void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd); int get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn); void set_reverse_hap_overlap(hap_overlaps* dest, hap_overlaps* source, uint32_t* types); +void print_hap_paf(ma_ug_t *ug, hap_overlaps* ovlp); #endif \ No newline at end of file diff --git a/hic.cpp b/hic.cpp index 18c1afb..c7438db 100644 --- a/hic.cpp +++ b/hic.cpp @@ -13364,6 +13364,11 @@ const char* aln) } **/ +void reduce_trans_chain(hc_links *link) +{ + ; +} + int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) { double index_time = yak_realtime(); @@ -13416,28 +13421,28 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) bub.round_id = 0; bub.n_round = 2; for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++) { - fprintf(stderr, "0********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); + ///fprintf(stderr, "0********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); identify_bubbles(idx->ug, &bub, idx->cov->t_ch->is_r_het); - fprintf(stderr, "1********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); + ///fprintf(stderr, "1********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); if(bub.round_id == 0) { - fprintf(stderr, "2********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); + ///fprintf(stderr, "2********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); collect_hc_links(sl.idx, &sl.hits, &link, &bub, &M); - fprintf(stderr, "3********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); + ///fprintf(stderr, "3********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); collect_hc_reverse_links(&link, idx->ug, &bub); - fprintf(stderr, "4********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); + ///fprintf(stderr, "4********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); } init_hic_p((ha_ug_index*)sl.idx, &sl.hits, &link, &bub, &back_hc_edge, &M, &hap, 0); - fprintf(stderr, "5********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); + ///fprintf(stderr, "5********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); ///init_hic_p_new((ha_ug_index*)sl.idx, &sl.hits, idx->link, &bub, &back_hc_edge, &M); reset_H_partition(&hap, (bub.round_id == 0? 1 : 0)); - fprintf(stderr, "6********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); + ///fprintf(stderr, "6********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); init_contig_partition(&hap, idx, &bub, &link); - fprintf(stderr, "7********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); + ///fprintf(stderr, "7********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); phasing_improvement(&hap, &(hap.g_p), idx, &bub, &link); - fprintf(stderr, "8********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); + ///fprintf(stderr, "8********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); label_unitigs(&(hap.g_p), idx->ug); - fprintf(stderr, "9********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); + ///fprintf(stderr, "9********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); ///print_hc_links(idx->link, 0, &hap); } diff --git a/partig.cpp b/partig.cpp index ab26f79..974486c 100644 --- a/partig.cpp +++ b/partig.cpp @@ -721,6 +721,7 @@ int8_t *s) } } + void clean_ovlp(ma_ug_t *ug, asg_t *read_g, hap_overlaps_list* ha, pt_g_t *pg, int8_t* s) { uint32_t v, i, k, qn, tn, types[4]; @@ -751,11 +752,35 @@ void clean_ovlp(ma_ug_t *ug, asg_t *read_g, hap_overlaps_list* ha, pt_g_t *pg, i } } + + // for (v = 0; v < ha->num; v++) + // { + // for (i = k = 0; i < ha->x[v].a.n; i++) + // { + // if(ha->x[v].a.a[i].status != FLIP) + // { + // if(s[ha->x[v].a.a[i].xUid]*s[ha->x[v].a.a[i].yUid] == -1) + // { + // fprintf(stderr, "\ns[0]=%d, s[1]=%d\n", s[ha->x[v].a.a[i].xUid], s[ha->x[v].a.a[i].yUid]); + // print_hap_paf(ug, &(ha->x[v].a.a[i])); + // } + + // continue; + // } + + + // ha->x[v].a.a[k] = ha->x[v].a.a[i]; + // k++; + // } + // ha->x[v].a.n = k; + // } + for (v = 0; v < ha->num; v++) { for (i = k = 0; i < ha->x[v].a.n; i++) { if(ha->x[v].a.a[i].status != FLIP) continue; + ha->x[v].a.a[k] = ha->x[v].a.a[i]; k++; } @@ -775,28 +800,31 @@ void pt_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, ma_ug_t *ug, asg_t *re pt_cc(pg->e); s = pt_solve_core(&opt, pg->e); + if(asm_opt.flag & HA_F_PARTITION) + { + MALLOC(buf, pg->e->ma.n); // FIXME: this is over-allocation for convenience + for (i = 0; i < pg->e->n_seq; ++i) { + uint64_t z[2]; + uint32_t o = pg->e->idx.a[i] >> 32; + uint32_t n = (uint32_t)pg->e->idx.a[i], j; - MALLOC(buf, pg->e->ma.n); // FIXME: this is over-allocation for convenience - for (i = 0; i < pg->e->n_seq; ++i) { - uint64_t z[2]; - uint32_t o = pg->e->idx.a[i] >> 32; - uint32_t n = (uint32_t)pg->e->idx.a[i], j; + set_trio_flag(ug, read_g, i, trio_flag, t_ch, ovlp, pg->e, s[i]); - set_trio_flag(ug, read_g, i, trio_flag, t_ch, ovlp, pg->e, s[i]); - - pg->info.a[i].s = s[i]; - for (j = 0; j < n; ++j) { - const pt_match1_t *m = &pg->e->ma.a[o + j]; - buf[j] = (uint64_t)((uint32_t)-1 - m->w) << 32 | (o + j); + pg->info.a[i].s = s[i]; + for (j = 0; j < n; ++j) { + const pt_match1_t *m = &pg->e->ma.a[o + j]; + buf[j] = (uint64_t)((uint32_t)-1 - m->w) << 32 | (o + j); + } + radix_sort_gfa64(buf, buf + n); + for (j = 0, z[0] = z[1] = 0; j < n; ++j) { + const pt_match1_t *m = &pg->e->ma.a[(uint32_t)buf[j]]; + if (s[m->sid[1]] > 0) z[0] += m->w; + else if (s[m->sid[1]] < 0) z[1] += m->w; + } + pg->info.a[i].m[0] = z[0], pg->info.a[i].m[1] = z[1]; } - radix_sort_gfa64(buf, buf + n); - for (j = 0, z[0] = z[1] = 0; j < n; ++j) { - const pt_match1_t *m = &pg->e->ma.a[(uint32_t)buf[j]]; - if (s[m->sid[1]] > 0) z[0] += m->w; - else if (s[m->sid[1]] < 0) z[1] += m->w; - } - pg->info.a[i].m[0] = z[0], pg->info.a[i].m[1] = z[1]; } + clean_ovlp(ug, read_g, ovlp, pg, s);