diff --git a/CommandLines.cpp b/CommandLines.cpp index 5a658f5..c4dccf0 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -136,7 +136,9 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->purge_level_primary = 2; asm_opt->purge_level_trio = 0; asm_opt->purge_simi_rate = 0.75; + asm_opt->purge_simi_rate_hic = 0.85; asm_opt->purge_overlap_len = 1; + asm_opt->purge_overlap_len_hic = 50; asm_opt->recover_atg_cov_min = -1024; asm_opt->recover_atg_cov_max = INT_MAX; asm_opt->hom_global_coverage = -1; diff --git a/CommandLines.h b/CommandLines.h index ba03f03..a431de5 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -67,6 +67,7 @@ typedef struct { int purge_level_primary; int purge_level_trio; int purge_overlap_len; + int purge_overlap_len_hic; int recover_atg_cov_min; int recover_atg_cov_max; int hom_global_coverage; @@ -76,6 +77,7 @@ typedef struct { float min_drop_rate; float max_drop_rate; float purge_simi_rate; + float purge_simi_rate_hic; long long small_pop_bubble_size; long long large_pop_bubble_size; diff --git a/Overlaps.cpp b/Overlaps.cpp index b5cb80e..796f77e 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -178,6 +178,7 @@ void ma_ug_destroy(ma_ug_t *ug) } free(ug->u.a); asg_destroy(ug->g); + kv_destroy(ug->occ); free(ug); } @@ -10966,7 +10967,7 @@ void set_pre_uid(buf_t* b, hc_links* link, ma_ug_t *ug) { u = &(ug->u.a[b->b.a[i]>>1]); if(u->m == 0) continue; - for (k = 0; k < u->m; k++) + for (k = 0; k < u->n; k++) { rId = u->a[k]>>33; if(link->u_idx[rId] == (uint32_t)-1) continue; @@ -10980,7 +10981,7 @@ void set_pre_uid(buf_t* b, hc_links* link, ma_ug_t *ug) } -void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug) +void collect_reverse_unitigs_back(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug) { uint32_t k, m; uint64_t d = (uint64_t)-1; @@ -10999,6 +11000,221 @@ void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug } +void collect_reverse_unitigs_back(buf_t* b_0, uint32_t b_0_uLen, buf_t* b_1, uint32_t b_1_uLen, hc_links* link, ma_ug_t *ug) +{ + uint32_t b_0_i, b_0_i_l, b_0_k, b_1_i, b_1_i_l, b_1_k, pre_0, pre_1, rId_0, rId_1; + uint64_t d = (uint64_t)-1; + ma_utg_t* u_b_0 = NULL; + ma_utg_t* u_b_1 = NULL; + if(b_0->b.n == 0 || b_1->b.n == 0) return; + if(b_0_uLen == 0 || b_1_uLen == 0) return; + + for (b_0_i = b_0_i_l = 0, pre_0 = (uint32_t)-1; b_0_i < b_0->b.n; b_0_i++) + { + u_b_0 = &(ug->u.a[b_0->b.a[b_0_i]>>1]); + if(u_b_0->n == 0) continue; + for (b_0_k = 0; b_0_k < u_b_0->n; b_0_k++) + { + rId_0 = u_b_0->a[b_0_k]>>33; + if(link->u_idx[rId_0] == (uint32_t)-1) continue; + if(pre_0 == link->u_idx[rId_0]) continue; + pre_0 = link->u_idx[rId_0]; + b_0_i_l++; + if(b_0_i_l > b_0_uLen) return; + + for (b_1_i = b_1_i_l = 0, pre_1 = (uint32_t)-1; b_1_i < b_1->b.n; b_1_i++) + { + u_b_1 = &(ug->u.a[b_1->b.a[b_1_i]>>1]); + if(u_b_1->n == 0) continue; + for (b_1_k = 0; b_1_k < u_b_1->n; b_1_k++) + { + rId_1 = u_b_1->a[b_1_k]>>33; + if(link->u_idx[rId_1] == (uint32_t)-1) continue; + if(pre_1 == link->u_idx[rId_1]) continue; + pre_1 = link->u_idx[rId_1]; + b_1_i_l++; + if(b_1_i_l > b_1_uLen) goto b_1_i_end; + + push_hc_edge(&(link->a.a[pre_0]), pre_1, 1, 1, &d); + push_hc_edge(&(link->a.a[pre_1]), pre_0, 1, 1, &d); + } + } + + b_1_i_end:; + } + } +} + +inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t* len_thre, uint64_t* occ) +{ + if(len_thre && occ)(*occ) = (uint64_t)-1; + uint32_t ori, uid, v, nv, l, k, idx; + uint32_t *a = b->b.a, a_n = b->b.n; + uint32_t u_i, r_i, len, p_v; + asg_arc_t *av = NULL; + ma_utg_t* u = NULL; + for (u_i = r_i = len = idx = 0, p_v = (uint32_t)-1; u_i < a_n; u_i++) + { + uid = a[u_i] >> 1; + ori = a[u_i] & 1; + u = &(ug->u.a[uid]); + if(u->n == 0) continue; + if(ori == 1) + { + for (r_i = 0; r_i < u->n; r_i++, idx++) + { + v = ((uint64_t)((u->a[u->n - r_i - 1])^(uint64_t)(0x100000000)))>>32; + ///w = ((uint64_t)((u->a[u->n - x->r_i - 2])^(uint64_t)(0x100000000)))>>32; + if(p_v == (uint32_t)-1) + { + p_v = v; + continue; + } + + av = asg_arc_a(read_sg, p_v); + nv = asg_arc_n(read_sg, p_v); + l = 0; + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + if(av[k].v == v) + { + l = asg_arc_len(av[k]); + break; + } + } + if(k == nv) fprintf(stderr, "ERROR\n"); + len += l; + if(len_thre && occ && len >= (*len_thre) && (*occ) != (uint64_t)-1) + { + (*occ) = idx - 1; + return len; + } + p_v = v; + } + } + else + { + for (r_i = 0; r_i < u->n; r_i++, idx++) + { + v = ((uint64_t)(u->a[r_i]))>>32; + ///w = ((uint64_t)(u->a[x->r_i + 1]))>>32; + if(p_v == (uint32_t)-1) + { + p_v = v; + continue; + } + + + av = asg_arc_a(read_sg, p_v); + nv = asg_arc_n(read_sg, p_v); + l = 0; + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + if(av[k].v == v) + { + l = asg_arc_len(av[k]); + break; + } + } + if(k == nv) fprintf(stderr, "ERROR\n"); + len += l; + if(len_thre && occ && len >= (*len_thre) && (*occ) != (uint64_t)-1) + { + (*occ) = idx - 1; + return len; + } + p_v = v; + } + } + } + + if(p_v != (uint32_t)-1) + { + len += read_sg->seq[p_v>>1].len; + if(len_thre && occ && len >= (*len_thre) && (*occ) != (uint64_t)-1) + { + (*occ) = idx - 1; + return len; + } + } + + return len; +} + +void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug, asg_t *read_sg) +{ + uint32_t b_0_i, b_0_k, b_1_i, b_1_k, pre_0, pre_1, rId_0, rId_1, ori_0, ori_1; + uint64_t d = (uint64_t)-1, len_0, len_1, thre_0, thre_1; + ma_utg_t* u_b_0 = NULL; + ma_utg_t* u_b_1 = NULL; + if(b_0->b.n == 0 || b_1->b.n == 0) return; + + len_0 = get_utg_len(b_0, ug, read_sg, NULL, NULL); + len_1 = get_utg_len(b_1, ug, read_sg, NULL, NULL); + + len_0 = MIN(len_0, len_1); + get_utg_len(b_0, ug, read_sg, &len_0, &thre_0); + get_utg_len(b_1, ug, read_sg, &len_0, &thre_1); + + + for (b_0_i = len_0 = 0, pre_0 = (uint32_t)-1; b_0_i < b_0->b.n; b_0_i++) + { + ori_0 = b_0->b.a[b_0_i]&1; + u_b_0 = &(ug->u.a[b_0->b.a[b_0_i]>>1]); + if(u_b_0->n == 0) continue; + for (b_0_k = 0; b_0_k < u_b_0->n; b_0_k++) + { + len_0++; + if(len_0 > thre_0) return; + if(ori_0 == 1) + { + rId_0 = u_b_0->a[u_b_0->n - b_0_k - 1]>>33; + } + else + { + rId_0 = u_b_0->a[b_0_k]>>33; + } + + if(link->u_idx[rId_0] == (uint32_t)-1) continue; + if(pre_0 == link->u_idx[rId_0]) continue; + pre_0 = link->u_idx[rId_0]; + + for (b_1_i = len_1 = 0, pre_1 = (uint32_t)-1; b_1_i < b_1->b.n; b_1_i++) + { + ori_1 = b_1->b.a[b_1_i]&1; + u_b_1 = &(ug->u.a[b_1->b.a[b_1_i]>>1]); + if(u_b_1->n == 0) continue; + for (b_1_k = 0; b_1_k < u_b_1->n; b_1_k++) + { + len_1++; + if(len_1 > thre_1) goto b_1_i_end; + if(ori_1 == 1) + { + rId_1 = u_b_1->a[u_b_1->n - b_1_k - 1]>>33; + } + else + { + rId_1 = u_b_1->a[b_1_k]>>33; + } + if(link->u_idx[rId_1] == (uint32_t)-1) continue; + if(pre_1 == link->u_idx[rId_1]) continue; + pre_1 = link->u_idx[rId_1]; + + + push_hc_edge(&(link->a.a[pre_0]), pre_1, 1, 1, &d); + push_hc_edge(&(link->a.a[pre_1]), pre_0, 1, 1, &d); + } + } + + b_1_i_end:; + } + } +} + + + int untig_asg_arc_simple_large_bubbles_trio(ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, uint32_t negative_flag, hc_links* link) { @@ -11128,7 +11344,7 @@ long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, uint32_t negativ asg_seq_drop(g, buffer.b.a[k]>>1); } - if(link) collect_reverse_unitigs(&b_0, &b_1, link, ug); + if(link) collect_reverse_unitigs(&b_0, &b_1, link, ug, read_sg); is_hap++; } @@ -11775,20 +11991,59 @@ 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 init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num) +{ + kv_malloc(link->a, ug_num); link->a.n = ug_num; + kv_malloc(link->enzymes, ug_num); link->enzymes.n = ug_num; + uint64_t i; + for (i = 0; i < link->a.n; i++) + { + kv_init(link->a.a[i].e); + kv_init(link->a.a[i].f); + } + MALLOC(link->u_idx, r_num); + memset(link->u_idx, -1, r_num*sizeof(uint32_t)); +} + +void destory_hc_links(hc_links* link) +{ + uint64_t i; + for (i = 0; i < link->a.n; i++) + { + kv_destroy(link->a.a[i].e); + kv_destroy(link->a.a[i].f); + } + kv_destroy(link->a); + free(link->u_idx); +} + + void output_hic_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, -ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, int max_hang, -int min_ovlp) +ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist, +long long tipsLen, float tip_drop_ratio, long long stops_threshold, +R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp) { kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); + ma_ug_t *ug = NULL, *copy = NULL; + ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); - ma_ug_t *ug = NULL; - ug = ma_ug_gen(sg); + + hc_links link; + init_hc_links(&link, ug->g->n_seq, R_INF.total_reads); + copy = copy_untig_graph(ug); + asm_opt.purge_overlap_len = asm_opt.purge_overlap_len_hic; + asm_opt.purge_simi_rate = asm_opt.purge_simi_rate_hic; + adjust_utg_by_primary(©, sg, TRIO_THRES, sources, reverse_sources, coverage_cut, + bubble_dist, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, + max_hang, min_ovlp, &new_rtg_edges, &link); + ma_ug_destroy(copy); + + + new_rtg_edges.a.n = 0; ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp); - - - /** fprintf(stderr, "Writing raw unitig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+25); @@ -11813,11 +12068,90 @@ int min_ovlp) **/ classify_untigs(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, &new_rtg_edges, max_hang, min_ovlp); - hic_analysis(ug); + hic_analysis(ug, sg, &link); + destory_hc_links(&link); ma_ug_destroy(ug); kv_destroy(new_rtg_edges.a); } +ma_ug_t* merge_utg(ma_ug_t **dest, ma_ug_t **src) +{ + asg_t *g_d = (*dest)->g, *g_s = (*src)->g; + uint64_t occ_d = g_d->n_seq, occ_s = g_s->n_seq, i; + asg_arc_t *p = NULL; + g_d->is_srt = g_d->is_symm = 0; + + for (i = 0; i < occ_s; i++) + { + asg_seq_set(g_d, i+occ_d, g_s->seq[i].len, g_s->seq[i].del); + g_d->seq[i+occ_d].c = g_s->seq[i].c; + } + + g_d->seq_vis = (uint8_t*)realloc(g_d->seq_vis, g_d->n_seq*2*sizeof(uint8_t)); + + + for (i = 0; i < g_s->n_arc; i++) + { + p = asg_arc_pushp(g_d); + (*p) = g_s->arc[i]; + p->ul += (occ_d<<33); + p->v += (occ_d<<1); + } + + asg_cleanup(g_d); + g_d->r_seq = g_d->n_seq; + + if(g_s->n_F_seq > 0 && g_s->F_seq) + { + uint64_t n_F_seq = g_d->n_F_seq + g_s->n_F_seq; + g_d->F_seq = (ma_utg_t*)realloc(g_d->F_seq, n_F_seq*sizeof(ma_utg_t)); + memcpy(g_d->F_seq + g_d->n_F_seq, g_s->F_seq, g_s->n_F_seq*sizeof(ma_utg_t)); + g_d->n_F_seq = n_F_seq; + free(g_s->F_seq); + g_s->F_seq = NULL; + g_s->n_F_seq = 0; + } + + ma_utg_v *u_d = &((*dest)->u), *u_s = &((*src)->u); + if(u_s->n > 0) + { + uint64_t n = u_d->n + u_s->n; + u_d->a = (ma_utg_t*)realloc(u_d->a, n*sizeof(ma_utg_t)); + memcpy(u_d->a + u_d->n, u_s->a, u_s->n*sizeof(ma_utg_t)); + u_d->n = u_d->m = n; + free(u_s->a); + u_s->a = NULL; + u_s->n = u_s->m = 0; + } + + kv_push(uint64_t, (*dest)->occ, occ_d); + kv_push(uint64_t, (*dest)->occ, occ_s); + + ma_ug_destroy(*src); + return (*dest); +} + +void benchmark_hic_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, +ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist, +long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, +float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp) +{ + ma_ug_t *ug_1 = output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, + reverse_sources, bubble_dist, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, + chimeric_rate, drop_ratio, max_hang, min_ovlp, 1); + + ma_ug_t *ug_2 = output_trio_unitig_graph(sg, coverage_cut, output_file_name, MOTHER, sources, + reverse_sources, bubble_dist, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, + chimeric_rate, drop_ratio, max_hang, min_ovlp, 1); + fprintf(stderr, "ug_1->u.n: %u, ug_2->u.n: %u\n", (uint32_t)ug_1->u.n, (uint32_t)ug_2->u.n); + ma_ug_t *ug = merge_utg(&ug_1, &ug_2); + fprintf(stderr, "ug->u.n: %u\n", (uint32_t)ug->u.n); + + hic_benchmark(ug, sg); + + ma_ug_destroy(ug); +} + void merge_unitig_content(ma_utg_t* collection, ma_ug_t* ug, asg_t* read_g, kvec_asg_arc_t_warp* edge) { if(collection->m == 0) return; @@ -12467,7 +12801,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hc_links* link) // #define NON_PLOID 1 if(flag == NON_PLOID) operation = CUT; - for (k = 0; k < b.b.n; k++) { g->seq[b.b.a[k]>>1].c = ALTER_LABLE; @@ -12486,9 +12819,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hc_links* link) } } - - if(link && operation != CUT) collect_reverse_unitigs(&b_0, &b_1, link, ug); - + if(link && operation != CUT) collect_reverse_unitigs(&b_0, &b_1, link, ug, read_sg); } } } @@ -12612,7 +12943,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_thre } } - if(link && operation != CUT) collect_reverse_unitigs(&b_0, &b_1, link, ug); + if(link && operation != CUT) collect_reverse_unitigs(&b_0, &b_1, link, ug, read_sg); break; } @@ -12799,7 +13130,7 @@ hc_links* link) asg_seq_drop(g, b.b.a[k]>>1); } - if(link) collect_reverse_unitigs(&b_0, &b_1, link, ug); + if(link) collect_reverse_unitigs(&b_0, &b_1, link, ug, read_sg); is_hap++; } @@ -13125,7 +13456,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_ asg_seq_drop(g, b.b.a[k]>>1); } - if(link) collect_reverse_unitigs(&b_0, &b_1, link, ug); + if(link) collect_reverse_unitigs(&b_0, &b_1, link, ug, read_sg); ///lable the primary one b_0.b.n = 0; @@ -13583,12 +13914,12 @@ float drop_ratio, hc_links* link) if(just_bubble_pop == 0) { ///need consider tangles - asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, link); + asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, link); asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, link); asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, link); asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, link); detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex); - + if(round != T_ROUND) { unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, @@ -13597,21 +13928,17 @@ float drop_ratio, hc_links* link) } cur_cons = get_graph_statistic(g); } - untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, DROP, link); + untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, DROP, link); if(just_bubble_pop == 0) { cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 2); } - - resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, (uint32_t)-1, drop_ratio); drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex); - unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, reverse_sources, 0, 1); - if(round > 0) { if(round != T_ROUND) @@ -14382,15 +14709,16 @@ int debug_untig_length(ma_ug_t *g, uint32_t tipsLen, const char* name) -void output_trio_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, +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, long long bubble_dist, 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) +float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench) { char* gfa_name = (char*)malloc(strlen(output_file_name)+100); sprintf(gfa_name, "%s.%s.p_ctg.gfa", output_file_name, (flag==FATHER?"hap1":"hap2")); fprintf(stderr, "Writing %s to disk... \n", gfa_name); - FILE* output_file = fopen(gfa_name, "w"); + FILE* output_file = NULL; + if(is_bench == 0) output_file = fopen(gfa_name, "w"); ma_ug_t *ug = NULL; ug = ma_ug_gen(sg); @@ -14406,6 +14734,12 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp) ///debug_untig_length(ug, tipsLen, gfa_name); ///print_untig_by_read(ug, "m64011_190901_095311/125831121/ccs", 2310925, "end"); ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp); + if(is_bench) + { + free(gfa_name); + kv_destroy(new_rtg_edges.a); + return ug; + } ma_ug_print(ug, &R_INF, sg, coverage_cut, sources, ruIndex, (flag==FATHER?"h1tg":"h2tg"), output_file); fclose(output_file); @@ -14425,6 +14759,7 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp) free(gfa_name); ma_ug_destroy(ug); kv_destroy(new_rtg_edges.a); + return NULL; } @@ -14668,12 +15003,50 @@ static void asg_bub_backtrack_primary(asg_t *g, uint32_t v0, buf_t *b) } +// in a resolved bubble, mark unused vertices and arcs as "reduced" +uint64_t asg_bub_backtrack_primary_length(asg_t *g, uint32_t v0, buf_t *b) +{ + uint32_t i, v, u, nv; + uint64_t len = 0; + ///b->S.a[0] is the sink of this bubble + asg_arc_t *av = NULL; + + ///v is the sink of this bubble + v = b->S.a[0]; + len = 0; + ///recover node + while (1) + { + u = b->a[v].p; // u->v + if(u == v0) break; + if(v == b->S.a[0]) + { + len += g->seq[u>>1].len; + } + else + { + nv = asg_arc_n(g, u); + av = asg_arc_a(g, u); + for (i = 0; i < nv; ++i) + { + if(av[i].del) continue; + if(av[i].v == v) break; + } + ///if(i == nv) fprintf(stderr, "ERROR\n"); + len += (uint32_t)av[i].ul; + } + v = u; + } + + return len; +} + // pop bubbles from vertex v0; the graph MJUST BE symmetric: if u->v present, v'->u' must be present as well uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, int max_dist, buf_t *b, -uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop) +uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len) { uint32_t i, n_pending = 0, is_first = 1, cur_m, cur_c, cur_np, cur_nc, to_replace, n_tips, tip_end; uint64_t n_pop = 0; @@ -14897,6 +15270,7 @@ uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop) if(is_pop) asg_bub_backtrack_primary(g, v0, b); + if(path_base_len) (*path_base_len) = asg_bub_backtrack_primary_length(g, v0, b); n_pop = 1; pop_reset: @@ -14929,7 +15303,7 @@ int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_fla for (i = 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); + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL); } 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); @@ -20443,13 +20817,13 @@ uint32_t positive_flag, uint32_t negative_flag) v = beg; if((!g->seq[v>>1].del)&&(g->seq[v>>1].c!=ALTER_LABLE)&&get_real_length(g, v, NULL)>=2) { - n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1); + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL); } v = end^1; if((!g->seq[v>>1].del)&&(g->seq[v>>1].c!=ALTER_LABLE)&&get_real_length(g, v, NULL)>=2) { - n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1); + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL); } @@ -20464,7 +20838,7 @@ uint32_t positive_flag, uint32_t negative_flag) { v = v|k; if(get_real_length(g, v, NULL)<=1) continue; - n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1); + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL); } } @@ -23124,7 +23498,32 @@ R_to_U* ruIndex) } -void append_utg(ma_ug_t* ptg, ma_ug_t* atg) +void reset_reverse_unitigs(hc_links* link, ma_utg_t *u) +{ + uint32_t k = 0, i = 0, rId, pre = (uint32_t)-1; + hc_edge *e = NULL; + if(u->n == 0 || u->m == 0) return; + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + if(link->u_idx[rId] == (uint32_t)-1) continue; + if(pre == link->u_idx[rId]) continue; + pre = link->u_idx[rId]; + if(link->a.a[pre].f.n == 0) continue; + + for (i = 0; i < link->a.a[pre].f.n; i++) + { + if(link->a.a[pre].f.a[i].del) continue; + link->a.a[pre].f.a[i].del = 1; + e = get_hc_edge(link, link->a.a[pre].f.a[i].uID, pre, 1); + if(e == NULL) continue; + e->del = 1; + } + } +} + + +void append_utg(ma_ug_t* ptg, ma_ug_t* atg, hc_links* link) { uint64_t num_nodes = 0; asg_t* nsg = atg->g; @@ -23149,6 +23548,7 @@ void append_utg(ma_ug_t* ptg, ma_ug_t* atg) for (v = 0; v < atg->g->n_seq; ++v) { if(atg->g->seq[v].del || atg->u.a[v].m == 0) continue; + if(link) reset_reverse_unitigs(link, &(atg->u.a[v])); p = &(ptg->u.a[ptg->u.n]); p->len = atg->u.a[v].len; @@ -23203,7 +23603,7 @@ void print_utg_coverage(ma_ug_t *ug, ma_sub_t* coverage_cut, uint32_t v, ma_hit_ } void recover_utg_by_coverage(ma_ug_t **ptg, asg_t* read_g, ma_sub_t* coverage_cut, -ma_hit_t_alloc* sources, R_to_U* ruIndex) +ma_hit_t_alloc* sources, R_to_U* ruIndex, hc_links* link) { if(asm_opt.recover_atg_cov_min == -1) return; if(asm_opt.recover_atg_cov_max == -1) return; @@ -23290,7 +23690,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex) { asg_cleanup(nsg); asg_symm(nsg); - append_utg(*ptg, atg); + append_utg(*ptg, atg, link); n_vtx = read_g->n_seq; for (v = 0; v < n_vtx; v++) @@ -23361,6 +23761,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link) } drop_semi_circle((*ug), nsg, read_g, reverse_sources, ruIndex); + asg_cleanup(nsg); adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex); @@ -23392,13 +23793,11 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link) 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, min_ovlp, 0, 0, 1, NULL); renew_utg(ug, read_g, new_rtg_edges); - rescue_contained_reads_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang, min_ovlp, 0, 10, 0, 1, NULL, NULL); renew_utg(ug, read_g, new_rtg_edges); @@ -23411,6 +23810,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link) purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, bubble_dist, drop_ratio, just_contain, 0, link); + delete_useless_nodes(ug); renew_utg(ug, read_g, new_rtg_edges); } @@ -23473,7 +23873,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link) __func__, asm_opt.recover_atg_cov_min); } - recover_utg_by_coverage(ug, read_g, coverage_cut, sources, ruIndex); + recover_utg_by_coverage(ug, read_g, coverage_cut, sources, ruIndex, link); } @@ -23530,6 +23930,7 @@ long long tipsLen, R_to_U* ruIndex, int max_hang, int min_ovlp) kv_destroy(new_rtg_edges.a); } + void output_contig_graph_primary(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, @@ -24591,7 +24992,7 @@ void lable_all_bubbles(asg_t *r_g, long long bubble_dist) ///if this is a bubble ///if(asg_bub_finder_with_del_advance(r_g, v, bubble_dist, &b) == 1) - if(asg_bub_pop1_primary_trio(r_g, NULL, v, bubble_dist, &b, (uint32_t)-1, (uint32_t)-1, 0)) + if(asg_bub_pop1_primary_trio(r_g, NULL, v, bubble_dist, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL)) { //beg is v, end is b.S.a[0] //note b.b include end, does not include beg @@ -27635,7 +28036,15 @@ ma_sub_t **coverage_cut_ptr, int debug_g) /*******************************for debug***************************************/ } - if (ha_opt_triobin(&asm_opt)) + if (ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt)) + { + char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); + sprintf(buf, "%s.hic.bench", output_file_name); + benchmark_hic_graph(sg, coverage_cut, buf, sources, reverse_sources, bubble_dist, + (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length); + free(buf); + } + else if (ha_opt_triobin(&asm_opt)) { char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); sprintf(buf, "%s.dip", output_file_name); @@ -27644,16 +28053,18 @@ ma_sub_t **coverage_cut_ptr, int debug_g) output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, reverse_sources, bubble_dist, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, - 0.05, 0.9, max_hang_length, mini_overlap_length); + 0.05, 0.9, max_hang_length, mini_overlap_length, 0); output_trio_unitig_graph(sg, coverage_cut, output_file_name, MOTHER, sources, reverse_sources, bubble_dist, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, - 0.05, 0.9, max_hang_length, mini_overlap_length); + 0.05, 0.9, max_hang_length, mini_overlap_length, 0); } else if(ha_opt_hic(&asm_opt)) { char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); sprintf(buf, "%s.hic", output_file_name); - output_hic_graph(sg, coverage_cut, buf, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length);; + output_hic_graph(sg, coverage_cut, buf, sources, reverse_sources, bubble_dist, + (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, + mini_overlap_length); free(buf); } else diff --git a/Overlaps.h b/Overlaps.h index 1fd741f..8419147 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -155,6 +155,7 @@ typedef struct { size_t n, m; ma_utg_t *a; } ma_utg_v; typedef struct { ma_utg_v u; asg_t *g; + kvec_t(uint64_t) occ; } ma_ug_t; typedef struct { @@ -743,13 +744,12 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold) stops_threshold, b_0) == LOOP) { return UNAVAILABLE; - } + } if(get_unitig(nsg, ug, v_1, &vEnd, &ELen_1, &tmp, &max_stop_nodeLen, &max_stop_baseLen, stops_threshold, b_1) == LOOP) { return UNAVAILABLE; } - if(ELen_0<=min_edge_length || ELen_1<=min_edge_length) return UNAVAILABLE; rIdContig b_max, b_min; @@ -770,7 +770,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold) uint32_t max_count = 0, min_count = 0; ma_utg_t *node_min = NULL, *node_max = NULL; - if(ug != NULL) { /*****************************label all unitigs****************************************/ @@ -786,7 +785,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold) } /*****************************label all unitigs****************************************/ - ///each unitig for (b_min.untigI = 0; b_min.untigI < b_min.b_0->b.n; b_min.untigI++) { @@ -820,7 +818,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold) } } } - /*****************************label all unitigs****************************************/ for (b_max.untigI = 0; b_max.untigI < b_max.b_0->b.n; b_max.untigI++) { @@ -833,7 +830,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold) } } /*****************************label all unitigs****************************************/ - } else { @@ -896,8 +892,6 @@ R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold) return NON_PLOID; } - - inline uint32_t check_different_haps_naive(asg_t *nsg, ma_ug_t *ug, asg_t *read_sg, uint32_t v_0, uint32_t v_1, ma_hit_t_alloc* reverse_sources, buf_t* b_0, buf_t* b_1, R_to_U* ruIndex, uint32_t min_edge_length, uint32_t stops_threshold) @@ -1061,7 +1055,7 @@ uint32_t is_bubble_check, uint32_t is_primary_check); uint32_t get_edge_from_source(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t query, uint32_t target, asg_arc_t* t); uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, int max_dist, buf_t *b, -uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop); +uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len); int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio); @@ -1089,7 +1083,8 @@ typedef struct{ hc_edge *a; }hc_edge_warp; - +void init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num); +void destory_hc_links(hc_links* link); void clean_primary_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen, @@ -1100,8 +1095,12 @@ ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link); -void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug); - +void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug, asg_t *read_sg); +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, long long bubble_dist, +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, int is_bench); #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index c7d2048..5818c4b 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -3292,7 +3292,7 @@ int asg_pop_bubble_purge_graph(asg_t *purge_g, int max_dist) for (i = 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(purge_g, NULL, v, max_dist, &b, (uint32_t)-1, DROP, 1); + n_pop += asg_bub_pop1_primary_trio(purge_g, NULL, v, max_dist, &b, (uint32_t)-1, DROP, 1, NULL); } free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); if (n_pop) asg_cleanup(purge_g); @@ -3928,23 +3928,23 @@ kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edg -void collect_reverse_unitig_pair(hc_links* link, ma_ug_t *ug, uint32_t b_0, uint32_t b_1) +void collect_reverse_unitig_pair(hc_links* link, ma_ug_t *ug, hap_overlaps* t) { - uint32_t i = 0, k = 0, rId_0, rId_1, pre_0 , pre_1; + uint32_t i = 0, k = 0, rId_0, rId_1, pre_0, pre_1, b_0 = t->xUid, b_1 = t->yUid; uint64_t d = (uint64_t)-1; ma_utg_t* u_b_0 = &(ug->u.a[b_0]); ma_utg_t* u_b_1 = &(ug->u.a[b_1]); if(u_b_0->n == 0) return; if(u_b_1->n == 0) return; - for (i = 0, pre_0 = (uint32_t)-1; i < u_b_0->n; i++) + for (i = t->x_beg_id, pre_0 = (uint32_t)-1; i < t->x_end_id; i++) { rId_0 = u_b_0->a[i]>>33; if(link->u_idx[rId_0] == (uint32_t)-1) continue; if(pre_0 == link->u_idx[rId_0]) continue; pre_0 = link->u_idx[rId_0]; - for (k = 0, pre_1 = (uint32_t)-1; k < u_b_1->n; k++) + for (k = t->y_beg_id, pre_1 = (uint32_t)-1; k < t->y_end_id; k++) { rId_1 = u_b_1->a[k]>>33; if(link->u_idx[rId_1] == (uint32_t)-1) continue; @@ -3958,13 +3958,16 @@ void collect_reverse_unitig_pair(hc_links* link, ma_ug_t *ug, uint32_t b_0, uint } -void collect_reverse_unitigs_purge(buf_t* b_0, hc_links* link, ma_ug_t *ug) +void collect_reverse_unitigs_purge(buf_t* b_0, hc_links* link, ma_ug_t *ug, hap_overlaps_list* all_ovlp) { if(b_0->b.n <= 1) return; uint32_t k; + int index = 0; for (k = 0; k < b_0->b.n - 1; k++) { - collect_reverse_unitig_pair(link, ug, b_0->b.a[k]>>1, b_0->b.a[k+1]>>1); + index = get_specific_hap_overlap(&(all_ovlp->x[b_0->b.a[k]>>1]), b_0->b.a[k]>>1, b_0->b.a[k+1]>>1); + if(index == -1) continue; + collect_reverse_unitig_pair(link, ug, &(all_ovlp->x[b_0->b.a[k]>>1].a.a[index])); } } @@ -3994,7 +3997,7 @@ hc_links* link) continue; } - if(link) collect_reverse_unitigs_purge(&b_0, link, ug); + if(link) collect_reverse_unitigs_purge(&b_0, link, ug, all_ovlp); purge_merge(purge_g, ug, all_ovlp, &b_0, ruIndex, reverse_sources, coverage_cut, read_g, position_index, u_buffer, tailIndex, prevIndex,max_hang, min_ovlp, edge, visit); } @@ -4347,7 +4350,7 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link) 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 = 1; all_ovlp.x[uId].a.a[i].status = DELETE; - if(link) collect_reverse_unitig_pair(link, ug, all_ovlp.x[uId].a.a[i].xUid, all_ovlp.x[uId].a.a[i].yUid); + if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i])); } if(all_ovlp.x[uId].a.a[i].type == XCY) @@ -4356,7 +4359,7 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link) 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 = 1; all_ovlp.x[uId].a.a[i].status = DELETE; - if(link) collect_reverse_unitig_pair(link, ug, all_ovlp.x[uId].a.a[i].xUid, all_ovlp.x[uId].a.a[i].yUid); + if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i])); } ///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i])); } diff --git a/hic.cpp b/hic.cpp index 461faf0..2f7b737 100644 --- a/hic.cpp +++ b/hic.cpp @@ -1,5 +1,6 @@ #define __STDC_LIMIT_MACROS #include "float.h" +#include #include "hic.h" #include "htab.h" #include "assert.h" @@ -13,8 +14,6 @@ #include "kdq.h" KSEQ_INIT(gzFile, gzread) KDQ_INIT(uint64_t) -#define kdq_clear(q) ((q)->count = (q)->front = 0) -#define kv_malloc(v, s) ((v).n = 0, (v).m = (s), MALLOC((v).a, (s))) #define HIC_COUNTER_BITS 12 #define HIC_MAX_COUNT ((1<index) free(bub->index); kv_destroy(bub->list); kv_destroy(bub->num); + kv_destroy(bub->pathLen); } -inline void get_bubbles(bubble_type* bub, uint64_t id, uint32_t* beg, uint32_t* sink, uint32_t** a, uint32_t* n) +inline void get_bubbles(bubble_type* bub, uint64_t id, uint32_t* beg, uint32_t* sink, uint32_t** a, uint32_t* n, uint64_t* pathBase) { (*a) = bub->list.a + bub->num.a[id] + 2; (*n) = bub->num.a[id+1] - bub->num.a[id] - 2; (*beg) = bub->list.a[bub->num.a[id]]; (*sink) = bub->list.a[bub->num.a[id] + 1]; + if(pathBase) (*pathBase) = bub->pathLen.a[id]; } void identify_bubbles(ma_ug_t* ug, bubble_type* bub) @@ -1402,7 +1420,8 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub) if (!ug->g->is_symm) asg_symm(ug->g); memset(bub, 0, sizeof(bubble_type)); uint32_t v, n_vtx = ug->g->n_seq * 2, tLen, i, mode = (((uint32_t)-1)<<2); - ///bub->ug = ug; + uint64_t pathLen; + bub->ug = ug; CALLOC(bub->index, n_vtx); for (i = 0; i < ug->g->n_seq; i++) { @@ -1412,7 +1431,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub) ug->g->seq[i].c = 0; } } - kv_init(bub->list); kv_init(bub->num); + kv_init(bub->list); kv_init(bub->num); kv_init(bub->pathLen); buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); for (i = 0, tLen = 1; i < ug->u.n; i++) tLen += ug->u.a[i].len; @@ -1421,7 +1440,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub) if(ug->g->seq[v>>1].del) continue; if(asg_arc_n(ug->g, v) < 2) continue; if((bub->index[v]&(uint32_t)3) != 0) continue; - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0)) + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL)) { //beg is v, end is b.S.a[0] //note b.b include end, does not include beg @@ -1441,8 +1460,9 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub) { if((bub->index[v]&(uint32_t)3) !=2) continue; kv_push(uint32_t, bub->num, bub->list.n); - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0)) + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen)) { + kv_push(uint64_t, bub->pathLen, pathLen); //beg is v, end is b.S.a[0] kv_push(uint32_t, bub->list, v); kv_push(uint32_t, bub->list, b.S.a[0]^1); @@ -1481,7 +1501,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub) uint32_t beg, sink, n, *a; for (i = 0; i < bub->num.n-1; i++) { - get_bubbles(bub, i, &beg, &sink, &a, &n); + get_bubbles(bub, i, &beg, &sink, &a, &n, NULL); for (v = 0; v < n; v++) { bub->index[(a[v]>>1)] = i; @@ -1510,7 +1530,7 @@ void print_bubbles(ma_ug_t* ug, bubble_type* bub, kvec_pe_hit* hits, hc_links* l fprintf(stderr, "[M::%s] # unitigs: %lu, # bases: %lu\n", __func__, bub->ug->u.n, tLen); for (i = 0, tLen = 0, t_utg = 0; i < bub->num.n-1; i++) { - get_bubbles(bub, i, &beg, &sink, &a, &n); + get_bubbles(bub, i, &beg, &sink, &a, &n, NULL); t_utg += n; for (k = 0; k < n; k++) { @@ -1634,19 +1654,20 @@ void print_bubbles(ma_ug_t* ug, bubble_type* bub, kvec_pe_hit* hits, hc_links* l free(flag); - // fprintf(stderr, "************bubble utgs************\n"); - // for (i = 0, tLen = 0, t_utg = 0; i < bub->num.n-1; i++) - // { - // get_bubbles(bub, i, &beg, &sink, &a, &n); - // t_utg += n; - // fprintf(stderr, "(%lu)\tbeg:utg%.6u\tsink:utg%.6u\n", i, (beg>>1)+1, (sink>>1)+1); - // for (k = 0; k < n; k++) - // { - // tLen +=bub->ug->u.a[(a[k]>>1)].len; - // fprintf(stderr, "utg%.6u,", (a[k]>>1)+1); - // } - // fprintf(stderr, "\n"); - // } + fprintf(stderr, "************bubble utgs************\n"); + uint64_t pathLen; + for (i = 0, tLen = 0, t_utg = 0; i < bub->num.n-1; i++) + { + get_bubbles(bub, i, &beg, &sink, &a, &n, &pathLen); + t_utg += n; + fprintf(stderr, "(%lu)\tbeg:utg%.6u\tsink:utg%.6u\tpathLen:%lu\n", i, (beg>>1)+1, (sink>>1)+1, pathLen); + for (k = 0; k < n; k++) + { + tLen +=bub->ug->u.a[(a[k]>>1)].len; + fprintf(stderr, "utg%.6u,", (a[k]>>1)+1); + } + fprintf(stderr, "\n"); + } // fprintf(stderr, "************het utgs************\n"); // for (i = 0; i < ug->g->n_seq; i++) @@ -1659,31 +1680,6 @@ void print_bubbles(ma_ug_t* ug, bubble_type* bub, kvec_pe_hit* hits, hc_links* l -void init_hc_links(hc_links* link, uint64_t ug_num) -{ - kv_malloc(link->a, ug_num); link->a.n = ug_num; - kv_malloc(link->enzymes, ug_num); link->enzymes.n = ug_num; - uint64_t i; - for (i = 0; i < link->a.n; i++) - { - kv_init(link->a.a[i].e); - kv_init(link->a.a[i].f); - } - MALLOC(link->u_idx, ug_num); - memset(link->u_idx, -1, ug_num*sizeof(uint32_t)); -} - -void destory_hc_links(hc_links* link) -{ - uint64_t i; - for (i = 0; i < link->a.n; i++) - { - kv_destroy(link->a.a[i].e); - kv_destroy(link->a.a[i].f); - } - kv_destroy(link->a); - free(link->u_idx); -} void push_hc_edge(hc_linkeage* x, uint64_t uID, int weight, int dir, uint64_t* d) { @@ -2011,7 +2007,7 @@ uint64_t get_LCA_bubble(uint32_t x, uint64_t xLen, uint32_t y, uint64_t yLen, ui uint64_t u, d = (uint64_t)-1, tmp; uint8_t rev; uint32_t root[2], a_n, *a; - get_bubbles(bub, bub->index[x>>1], &root[0], &root[1], &a, &a_n); + get_bubbles(bub, bub->index[x>>1], &root[0], &root[1], &a, &a_n, NULL); root[0] ^= 1; root[1] ^= 1; if(root[0] > root[1]) { @@ -2310,7 +2306,7 @@ static void worker_for_dis(void *data, long i, int tid) ///might be wrong if(bub->index[i] < bub->num.n && bub->index[u] < bub->num.n - && t->e.a[k].dis != (uint64_t)-1) + && bub->index[i] != bub->index[u] && t->e.a[k].dis != (uint64_t)-1) { continue; } @@ -2398,8 +2394,8 @@ void collect_hc_links(const ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, if(bub->index[end] > bub->num.n) continue; t_d = (uint64_t)-1; - push_hc_edge(&(link->a.a[beg]), end, 1, 0, &t_d); - push_hc_edge(&(link->a.a[end]), beg, 1, 0, &t_d); + push_hc_edge(&(link->a.a[beg]), end, 0, 0, &t_d); + push_hc_edge(&(link->a.a[end]), beg, 0, 0, &t_d); } uint32_t n_vtx = idx->ug->g->n_seq<<1, v; @@ -2500,13 +2496,57 @@ void set_reverse_links(uint32_t* bub, uint32_t n, kvec_t_u32_warp* reach, uint32 void collect_hc_reverse_links(hc_links* link, ma_ug_t* ug, bubble_type* bub) { - uint64_t i, k, d = 0; + uint64_t i, j, k, d = 0, m, pre; uint32_t beg, sink, n, v, *a = NULL; kvec_t_u32_warp stack, result; + hc_edge *e = NULL; kv_init(stack.a); kv_init(result.a); for (i = 0; i < bub->num.n-1; i++) { - get_bubbles(bub, i, &beg, &sink, &a, &n); + get_bubbles(bub, i, &beg, &sink, &a, &n, NULL); + for (k = 0; k < n; k++) + { + v = a[k]>>1; + for (j = 0; j < link->a.a[v].f.n; j++) + { + if(link->a.a[v].f.a[j].del) continue; + e = get_hc_edge(link, link->a.a[v].f.a[j].uID, v, 1); + if(e == NULL) fprintf(stderr, "ERROR\n"); + e->del = 1; + } + link->a.a[v].f.n = 0; + } + + v = beg>>1; + if(bub->index[v] > bub->num.n) + { + for (j = 0; j < link->a.a[v].f.n; j++) + { + if(link->a.a[v].f.a[j].del) continue; + e = get_hc_edge(link, link->a.a[v].f.a[j].uID, v, 1); + if(e == NULL) fprintf(stderr, "ERROR\n"); + e->del = 1; + } + link->a.a[v].f.n = 0; + } + + v = sink>>1; + if(bub->index[v] > bub->num.n) + { + for (j = 0; j < link->a.a[v].f.n; j++) + { + if(link->a.a[v].f.a[j].del) continue; + e = get_hc_edge(link, link->a.a[v].f.a[j].uID, v, 1); + if(e == NULL) fprintf(stderr, "ERROR\n"); + e->del = 1; + } + link->a.a[v].f.n = 0; + } + } + + for (i = 0; i < bub->num.n-1; i++) + { + get_bubbles(bub, i, &beg, &sink, &a, &n, NULL); if(n == 2) { push_hc_edge(&(link->a.a[a[0]>>1]), a[1]>>1, 1, 1, &d); @@ -2522,6 +2562,50 @@ void collect_hc_reverse_links(hc_links* link, ma_ug_t* ug, bubble_type* bub) } } kv_destroy(stack.a); kv_destroy(result.a); + + for (i = 0; i < link->a.n; i++) + { + for (k = m = 0; k < link->a.a[i].f.n; k++) + { + if(link->a.a[i].f.a[k].del) continue; + link->a.a[i].f.a[m] = link->a.a[i].f.a[k]; + m++; + } + link->a.a[i].f.n = m; + radix_sort_hc_edge_u(link->a.a[i].f.a, link->a.a[i].f.a + link->a.a[i].f.n); + + for (k = m = 0, pre = (uint64_t)-1; k < link->a.a[i].f.n; k++) + { + if(link->a.a[i].f.a[k].del) continue; + if(link->a.a[i].f.a[k].uID == pre) + { + if(link->a.a[i].f.a[k].dis == 0) link->a.a[i].f.a[m-1].dis = 0; + continue; + } + + pre = link->a.a[i].f.a[k].uID; + link->a.a[i].f.a[m] = link->a.a[i].f.a[k]; + m++; + } + link->a.a[i].f.n = m; + radix_sort_hc_edge_d(link->a.a[i].f.a, link->a.a[i].f.a + link->a.a[i].f.n); + } + + + + + + // hc_edge *e = NULL; + // for (i = 0; i < link->a.n; i++) + // { + // for (k = 0; k < link->a.a[i].f.n; k++) + // { + // if(link->a.a[i].f.a[k].del) continue; + // e = get_hc_edge(link, link->a.a[i].f.a[k].uID, i, 1); + // if(e == NULL) fprintf(stderr, "ERROR\n"); + // } + // } + } void write_hc_links(hc_links* link, kvec_pe_hit* hits, const char *fn) @@ -2552,6 +2636,20 @@ void write_hc_links(hc_links* link, kvec_pe_hit* hits, const char *fn) free(buf); } + +void write_hc_hits(kvec_pe_hit* hits, const char *fn) +{ + char *buf = (char*)calloc(strlen(fn) + 25, 1); + sprintf(buf, "%s.hic.lk", fn); + FILE* fp = fopen(buf, "w"); + + fwrite(&hits->a.n, sizeof(hits->a.n), 1, fp); + fwrite(hits->a.a, sizeof(pe_hit), hits->a.n, fp); + + fclose(fp); + free(buf); +} + int load_hc_links(hc_links* link, kvec_pe_hit* hits, const char *fn) { uint64_t k, flag = 0; @@ -2594,34 +2692,60 @@ int load_hc_links(hc_links* link, kvec_pe_hit* hits, const char *fn) return 1; } -void print_hc_links(hc_links* link) +int load_hc_hits(kvec_pe_hit* hits, const char *fn) +{ + uint64_t flag = 0; + char *buf = (char*)calloc(strlen(fn) + 25, 1); + sprintf(buf, "%s.hic.lk", fn); + + FILE* fp = NULL; + fp = fopen(buf, "r"); + if(!fp) return 0; + + kv_init(hits->a); + flag += fread(&hits->a.n, sizeof(hits->a.n), 1, fp); + hits->a.m = hits->a.n; MALLOC(hits->a.a, hits->a.n); + flag += fread(hits->a.a, sizeof(pe_hit), hits->a.n, fp); + + fclose(fp); + free(buf); + fprintf(stderr, "[M::%s::] ==> Hi-C linkages have been loaded\n", __func__); + return 1; +} + +void print_hc_links(hc_links* link, int dir) { uint64_t i, k; - for (i = 0; i < link->a.n; ++i) - { - for (k = 0; k < link->a.a[i].e.n; k++) - { - if(link->a.a[i].e.a[k].del) continue; - fprintf(stderr, "s-utg%.6d(%c)\td-utg%.6d(%c)\t%lu\t%c\te\n", - (int)(i+1), "01"[!!(link->a.a[i].e.a[k].dis&(uint64_t)2)], - (int)(link->a.a[i].e.a[k].uID+1), "01"[!!(link->a.a[i].e.a[k].dis&(uint64_t)1)], - link->a.a[i].e.a[k].dis == (uint64_t)-1? (uint64_t)-1 : link->a.a[i].e.a[k].dis>>3, - "fb"[!!(link->a.a[i].e.a[k].dis&(uint64_t)4)]); + if(dir == 0) + { + for (i = 0; i < link->a.n; ++i) + { + for (k = 0; k < link->a.a[i].e.n; k++) + { + if(link->a.a[i].e.a[k].del) continue; + fprintf(stderr, "s-utg%.6d(%c)\td-utg%.6d(%c)\t%lu\t%c\t%f\te\n", + (int)(i+1), "01"[!!(link->a.a[i].e.a[k].dis&(uint64_t)2)], + (int)(link->a.a[i].e.a[k].uID+1), "01"[!!(link->a.a[i].e.a[k].dis&(uint64_t)1)], + link->a.a[i].e.a[k].dis == (uint64_t)-1? (uint64_t)-1 : link->a.a[i].e.a[k].dis>>3, + "fb"[!!(link->a.a[i].e.a[k].dis&(uint64_t)4)], link->a.a[i].e.a[k].weight); + } } } + - - // for (i = 0; i < link->a.n; ++i) - // { - // for (k = 0; k < link->a.a[i].f.n; k++) - // { - // if(link->a.a[i].f.a[k].del) continue; - // fprintf(stderr, "utg%.6d\tutg%.6d\t%f\t-", - // (int)(i+1), (int)(link->a.a[i].f.a[k].uID+1), link->a.a[i].f.a[k].weight); - // if(link->a.a[i].f.a[k].weight > 2) fprintf(stderr,"\tcomplex"); - // fprintf(stderr,"\n"); - // } - // } + if(dir == 1) + { + for (i = 0; i < link->a.n; ++i) + { + for (k = 0; k < link->a.a[i].f.n; k++) + { + if(link->a.a[i].f.a[k].del) continue; + fprintf(stderr, "s-utg%.6d\td-utg%.6d\t%lu\te\n", + (int)(i+1), (int)(link->a.a[i].f.a[k].uID+1), link->a.a[i].f.a[k].dis); + } + } + } + } void normalize_hc_links(hc_links* link) @@ -2651,6 +2775,7 @@ hc_edge* get_rGraph_edge(min_cut_t* x, uint64_t src, uint64_t dest) return NULL; } + void init_min_cut_t(min_cut_t* x, hc_links* link, const bubble_type* bub, const ma_ug_t *ug) { uint64_t utg_num = link->a.n, i, k, u, v; @@ -3237,7 +3362,7 @@ uint64_t select_bmer(uint32_t src, uint64_t k, const bubble_type* bub, min_cut_t b_mer_d = d; if(bub->index[u] < bub->num.n && x->bmerVis.a[u] == 0) { - get_bubbles((bubble_type*)bub, bub->index[u], &beg, &sink, &a, &n); + get_bubbles((bubble_type*)bub, bub->index[u], &beg, &sink, &a, &n, NULL); for (j = 0; j < n; j++) x->bmerVis.a[(a[j]>>1)] = 1; } //must be here @@ -3265,7 +3390,7 @@ uint64_t select_bmer(uint32_t src, uint64_t k, const bubble_type* bub, min_cut_t b_mer_d = d; if(bub->index[u] < bub->num.n && x->bmerVis.a[u] == 0) { - get_bubbles((bubble_type*)bub, bub->index[u], &beg, &sink, &a, &n); + get_bubbles((bubble_type*)bub, bub->index[u], &beg, &sink, &a, &n, NULL); for (j = 0; j < n; j++) x->bmerVis.a[(a[j]>>1)] = 1; } //must be here @@ -3303,7 +3428,7 @@ uint32_t bub_only, uint32_t bub_extend) { if(bub_extend && x->bmerVis.a[v>>1] == 0) { - get_bubbles((bubble_type*)bub, bub->index[v>>1], &beg, &sink, &a, &n); + get_bubbles((bubble_type*)bub, bub->index[v>>1], &beg, &sink, &a, &n, NULL); for (j = 0; j < n; j++) x->bmerVis.a[(a[j]>>1)] = 1; } x->bmerVis.a[v>>1] = 1; @@ -3335,7 +3460,7 @@ void get_bmer_unitgs(min_cut_t* x, const bubble_type* bub, uint64_t k, uint64_t { uint32_t beg, sink, n, *a; if(bub->index[src] >= bub->num.n) return; - get_bubbles((bubble_type*)bub, bub->index[src], &beg, &sink, &a, &n); + get_bubbles((bubble_type*)bub, bub->index[src], &beg, &sink, &a, &n, NULL); memset(x->bmerVis.a, 0, x->bmerVis.n); ///select_bmer(src, k, bub, x, 1); select_bmer_distance(beg^1, k, bub, x, 1, 1); @@ -3667,7 +3792,7 @@ min_cut_t* m, hc_links* link, G_partition* x) kv_pushp(partition_warp, *x, &res); memset(flag, 0, ug->g->n_seq); uint32_t beg, sink, n, *a, i, k; - get_bubbles(bub, bid, &beg, &sink, &a, &n); + get_bubbles(bub, bid, &beg, &sink, &a, &n, NULL); res->full_bub = 0; trace_phase_path((ma_ug_t *)ug, beg, sink, b, m, flag, HAP1_LAB); trace_phase_path((ma_ug_t *)ug, beg, sink, b, m, flag, HAP2_LAB); @@ -3796,7 +3921,7 @@ void print_phased_bubble(G_partition* x, bubble_type* bub, uint32_t utg_n) } uint32_t n, *a; - get_bubbles(bub, bubID, &beg, &sink, &a, &n); + get_bubbles(bub, bubID, &beg, &sink, &a, &n, NULL); if(n > 2) fprintf(stderr, "complex\n"); } @@ -3845,9 +3970,9 @@ void debug_hc_links(ha_ug_index* idx, hc_links* link, sldat_t* sl, bubble_type* { ///print_hits(idx, &sl->hits, fn1); destory_hc_links(link); - init_hc_links(link, sl->idx->ug->g->n_seq); + init_hc_links(link, sl->idx->ug->g->n_seq, R_INF.total_reads); collect_hc_links(sl->idx, &sl->hits, link, bub); - print_hc_links(link); + print_hc_links(link, 0); } uint64_t get_hic_distance(pe_hit* hit, hc_links* link, const ha_ug_index* idx) @@ -3864,6 +3989,10 @@ uint64_t get_hic_distance(pe_hit* hit, hc_links* link, const ha_ug_index* idx) s_dir = (!!(t->e.a[k].dis&(uint64_t)2)); e_dir = (!!(t->e.a[k].dis&(uint64_t)1)); u_dis = (t->e.a[k].dis ==(uint64_t)-1? (uint64_t)-1 : t->e.a[k].dis>>3); + // if(s_uid == 24684 && s_pos == 124953 && e_uid == 16950 && e_pos == 93039) + // { + // fprintf(stderr, "*****************s_dir: %lu, e_dir: %lu, u_dis: %lu\n", s_dir, e_dir, u_dis); + // } if(u_dis == (uint64_t)-1) return (uint64_t)-1; if(s_dir == 1) s_pos = (long long)idx->ug->g->seq[s_uid].len - s_pos - 1; if(e_dir == 1) e_pos = (long long)idx->ug->g->seq[e_uid].len - e_pos - 1; @@ -3899,12 +4028,183 @@ hc_edge* get_hc_edge(hc_links* link, uint64_t src, uint64_t dest, uint64_t dir) return NULL; } -void init_hic_p(const ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub) +inline double get_trans(const ha_ug_index* idx, uint64_t x) { - uint64_t k, beg, end, t_d; - kvec_t(uint64_t) buf; - kv_init(buf); + return idx->a*(x/idx->frac) + idx->b; +} + +inline double get_trans_weight(const ha_ug_index* idx, uint64_t x) +{ + #define OFFSET_RATE 0.0000001 + #define OFFSET_RATE_THRES 16.118095551 + long double rate = get_trans(idx, x); + if(rate < 0) rate = 0; + rate += OFFSET_RATE; + if(rate > 0.5) rate = 0.5; + + double w = log((1-rate)/rate); + if(w < 0) w = 0; + if(w > OFFSET_RATE_THRES) w = OFFSET_RATE_THRES; + return w; +} + +void LeastSquare(uint64_t* vec, uint64_t len, ha_ug_index* idx, uint64_t med) +{ + #define SCAL_RATE 1000 + long double t1=0, t2=0, t3=0, t4=0, x, y, thres; + uint64_t i, len_convince, m; + + + for (i = 0; i < len; i += 4) + { + if(vec[i+1] > med) break; + x = ((double)(vec[i] + vec[i+1]))/2; + y = ((double)(vec[i+3]))/((double)(vec[i+2] + vec[i+3])); + + t1 += x*x; + t2 += x; + t3 += x*y; + t4 += y; + } + len_convince = i; + + if(t2 > t4) + { + idx->frac = t2/t4; + if(idx->frac > SCAL_RATE) idx->frac = idx->frac / SCAL_RATE; + } + + t1 /= (idx->frac*idx->frac); + t2 /= idx->frac; + t3 /= idx->frac; + idx->a = idx->b = 0; + if((t1*(len_convince>>2) - t2*t2) != 0) + { + idx->a = (t3*(len_convince>>2) - t2*t4) / (t1*(len_convince>>2) - t2*t2); + } + if((t1*(len_convince>>2) - t2*t2) != 0) + { + idx->b = (t1*t4 - t2*t3) / (t1*(len_convince>>2) - t2*t2); + } + + + if(len > 0) + { + vec[len - 3] = vec[len - 4] + (vec[1] - vec[0]); + } + if(len_convince >= len) return; + + thres = get_trans(idx, vec[len_convince] + vec[len_convince+1]); + fprintf(stderr, "len_convince: %lu, len: %lu, t1: %f, t2: %f, t3: %f, t4: %f, idx->a: %f, idx->b: %f, thres: %f\n", + len_convince, len, (double)t1, (double)t2, (double)t3, (double)t4, (double)idx->a, (double)idx->b, (double)thres); + + + + for (i = m = 0; i < len; i += 4) + { + x = ((double)(vec[i] + vec[i+1]))/2; + y = ((double)(vec[i+3]))/((double)(vec[i+2] + vec[i+3])); + if(vec[i+1] > med && y < thres) continue; + + t1 += x*x; + t2 += x; + t3 += x*y; + t4 += y; + m++; + } + + if(t2 > t4) + { + idx->frac = t2/t4; + if(idx->frac > SCAL_RATE) idx->frac = idx->frac / SCAL_RATE; + } + + len = m; + t1 /= (idx->frac*idx->frac); + t2 /= idx->frac; + t3 /= idx->frac; + ///fprintf(stderr, "len: %lu, t1: %f, t2: %f, t3: %f, t4: %f\n", len, (double)t1, (double)t2, (double)t3, (double)t4); + if((t1*(len>>2) - t2*t2) != 0) + { + idx->a = (t3*(len>>2) - t2*t4) / (t1*(len>>2) - t2*t2); + } + if((t1*(len>>2) - t2*t2) != 0) + { + idx->b = (t1*t4 - t2*t3) / (t1*(len>>2) - t2*t2); + } + fprintf(stderr, "len: %lu, t1: %f, t2: %f, t3: %f, t4: %f, idx->a: %f, idx->b: %f\n", + len, (double)t1, (double)t2, (double)t3, (double)t4, (double)idx->a, (double)idx->b); +} + + + +void weight_edges(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub) +{ + uint64_t k, i, shif = 64 - idx->uID_bits, beg, end, t_d; + hc_edge *e1 = NULL, *e2 = NULL; + double weight; + + for (i = 0; i < link->a.n; i++) + { + for (k = 0; k < link->a.a[i].e.n; k++) + { + if(link->a.a[i].e.a[k].del) continue; + link->a.a[i].e.a[k].weight = 0; + } + } + + for (k = 0; k < hits->a.n; ++k) + { + beg = ((hits->a.a[k].s<<1)>>shif); + end = ((hits->a.a[k].e<<1)>>shif); + + if(beg == end) continue; + if(bub->index[beg] > bub->num.n) continue; + if(bub->index[end] > bub->num.n) continue; + + t_d = get_hic_distance(&(hits->a.a[k]), link, idx); + if(t_d == (uint64_t)-1) continue; + + e1 = get_hc_edge(link, beg, end, 0); + e2 = get_hc_edge(link, end, beg, 0); + if(e1 == NULL || e2 == NULL) continue; + weight = get_trans_weight(idx, t_d); + + e1->weight += weight; + e2->weight += weight; + } +} + +void init_hic_p(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub) +{ + uint64_t k, i, m, beg, end, t_d, r_idx, f_idx, b_size, med = (uint64_t)-1, uID; + uint32_t b_beg, b_end, n, *a, b_cnt; + kvec_t(uint64_t) buf, buf_idx; + kv_init(buf); + kv_init(buf_idx); + + buf.n = 0; + for (i = 0; i < bub->num.n-1; i++) + { + get_bubbles(bub, i, &b_beg, &b_end, &a, &n, &b_size); + for (k = b_cnt = 0; k < n; k++) + { + b_cnt +=bub->ug->u.a[(a[k]>>1)].n; + } + ///too small + if(b_cnt <= 3) continue; + kv_push(uint64_t, buf, b_size); + } + radix_sort_hc64(buf.a, buf.a+buf.n); + med = buf.a[(uint64_t)(buf.n*0.75)]; + // if((buf.n&1) == 1) med = buf.a[buf.n>>1]; + // if((buf.n&1) == 0) med = (buf.a[buf.n>>1] + buf.a[(buf.n>>1)-1])/2; + + + + + buf.n = 0; for (k = 0; k < hits->a.n; ++k) { beg = ((hits->a.a[k].s<<1)>>(64 - idx->uID_bits)); @@ -3913,33 +4213,685 @@ void init_hic_p(const ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubbl if(bub->index[beg] > bub->num.n) continue; if(bub->index[end] > bub->num.n) continue; + + if(beg == end) { t_d = get_hic_distance(&(hits->a.a[k]), link, idx); + if(t_d == (uint64_t)-1) continue; + t_d = (t_d << 1); kv_push(uint64_t, buf, t_d); - continue; - } - - if(get_hc_edge(link, beg, end, 1)) - { - t_d = get_hic_distance(&(hits->a.a[k]), link, idx); - t_d = (t_d << 1) + 1; - kv_push(uint64_t, buf, t_d); } + else + { + if(get_hc_edge(link, beg, end, 1)) + { + t_d = get_hic_distance(&(hits->a.a[k]), link, idx); + if(t_d == (uint64_t)-1) continue; + t_d = (t_d << 1) + 1; + kv_push(uint64_t, buf, t_d); + } + } } ///might have bias, we may not use right linkage larger than trans rc linkage radix_sort_hc64(buf.a, buf.a+buf.n); + for (k = 0, r_idx = f_idx = (uint64_t)-1; k < buf.n; k++) + { + if((buf.a[k]&1) == 0) r_idx = k; + if((buf.a[k]&1) == 1) f_idx = k; + + ///fprintf(stderr, "%lu\t%lu\n", buf.a[k]>>1, buf.a[k]&1); + } + buf.n = MIN(r_idx, f_idx); + + for (k = 0; k < buf.n; k++) + { + if((buf.a[k]&1) == 1) + { + kv_push(uint64_t, buf_idx, buf.a[k]>>1); + } + } + + uint64_t cutoff = buf_idx.n * 0.9, t = buf_idx.n * 0.005, pre, step; + for (k = cutoff - t, pre = buf_idx.a[cutoff - t - 1], t_d = 0; k < cutoff + t; k++) + { + t_d += (buf_idx.a[k] - pre); + pre = buf_idx.a[k]; + } + step = (t_d/(t * 2))*100; + + buf_idx.n = 0; + uint64_t step_s = 0, step_e = step, cnt[2], k_end; + if(buf.n>0) step_s = buf.a[0]>>1, step_e = (buf.a[0]>>1) + step; + for (k = cnt[0] = cnt[1] = 0; k < buf.n; k++) + { + if((buf.a[k]>>1) < step_e && (buf.a[k]>>1) >= step_s) + { + cnt[buf.a[k]&1]++; + } + + if((buf.a[k]>>1) >= step_e) + { + while (!((buf.a[k]>>1) < step_e && (buf.a[k]>>1) >= step_s)) + { + // fprintf(stderr, "i: %u, step_s: %lu, step_e: %lu, cnt[0]: %lu, cnt[1]: %lu, rate: %f\n", + // (uint32_t)(buf_idx.n>>2), step_s, step_e, cnt[0], cnt[1], ((double)cnt[1])/(double)(cnt[0] + cnt[1])); + kv_push(uint64_t, buf_idx, step_s); + kv_push(uint64_t, buf_idx, step_e); + kv_push(uint64_t, buf_idx, cnt[0]); + kv_push(uint64_t, buf_idx, cnt[1]); + step_s += step; + step_e += step; + cnt[0] = cnt[1] = 0; + } + } + } + + if(cnt[0] > 0 || cnt[1] > 0) + { + kv_push(uint64_t, buf_idx, step_s); + kv_push(uint64_t, buf_idx, step_e); + kv_push(uint64_t, buf_idx, cnt[0]); + kv_push(uint64_t, buf_idx, cnt[1]); + } + + uint64_t smooth_step = 20, k_i, cnt_0; + for (k = 0; k+smooth_step < (buf_idx.n>>2); k++) + { + for (k_i = cnt_0 = 0; k_i < smooth_step; k_i++) + { + if(buf_idx.a[((k+k_i)<<2)+2] == 0 || + buf_idx.a[((k+k_i)<<2)+3] == 0) + { + cnt_0++; + } + } + + if(cnt_0 >= smooth_step * 0.3) + { + break; + } + } + + for (k_end = k; k < (buf_idx.n>>2); k++) + { + buf_idx.a[(k_end<<2)+1] = buf_idx.a[(k<<2)+1]; + buf_idx.a[(k_end<<2)+2] += buf_idx.a[(k<<2)+2]; + buf_idx.a[(k_end<<2)+3] += buf_idx.a[(k<<2)+3]; + } + + buf_idx.n = (k_end+1)<<2; + if(k_end == 0) buf_idx.n = 0; + + for (k = i = 0; k < buf_idx.n; k += 4) + { + if(buf_idx.a[k+2] == 0 && buf_idx.a[k+3] == 0) continue; + buf_idx.a[i] = buf_idx.a[k]; + buf_idx.a[i+1] = buf_idx.a[k+1]; + buf_idx.a[i+2] = buf_idx.a[k+2]; + buf_idx.a[i+3] = buf_idx.a[k+3]; + i += 4; + } + + buf_idx.n = i; + + // for (k = 0; k < buf_idx.n; k += 4) + // { + // fprintf(stderr, "step_s: %lu, step_e: %lu, rate: %f\n", buf_idx.a[k], buf_idx.a[k+1], + // (double)(buf_idx.a[k+3])/(double)(buf_idx.a[k+2] + buf_idx.a[k+3])); + // } + ///idx->step = step; + LeastSquare(buf_idx.a, buf_idx.n, idx, med); + + fprintf(stderr, "idx->a: %f, idx->b: %f, idx->frac: %f, med: %lu\n", + (double)idx->a, (double)idx->b, (double)idx->frac, med); + + + hc_edge *e = NULL; + for (i = 0; i < link->a.n; i++) + { + for (k = 0; k < link->a.a[i].f.n; k++) + { + if(link->a.a[i].f.a[k].del) continue; + if(link->a.a[i].f.a[k].dis != 0) continue; + uID = link->a.a[i].f.a[k].uID; + e = get_hc_edge(link, i, uID, 0); + if(e) e->del = 1; + e = get_hc_edge(link, uID, i, 0); + if(e) e->del = 1; + } + } + + for (i = 0; i < link->a.n; i++) + { + for (k = m = 0; k < link->a.a[i].e.n; k++) + { + if(link->a.a[i].e.a[k].del) continue; + if(link->a.a[i].e.a[k].dis == (uint64_t)-1) continue; + + link->a.a[i].e.a[m] = link->a.a[i].e.a[k]; + link->a.a[i].e.a[m].weight = 0; + m++; + } + link->a.a[i].e.n = m; + } + + weight_edges(idx, hits, link, bub); + + kv_destroy(buf); + kv_destroy(buf_idx); +} + +#define is_hap_set(i, Hap) (!!((Hap).hap[i]&((Hap).m[0]|(Hap).m[1]))) + + +double get_path_weight(uint32_t query, uint32_t v0, uint32_t root, bub_p_t_warp *b, hc_links* x) +{ + if(v0 == root) return 0; + uint32_t v, u; + hc_edge *p = NULL; + double weight = 0; + v = v0; + do { + u = b->a[v].p; // u->v + p = get_hc_edge(x, query>>1, v>>1, 0); + if(p) weight += p->weight; + v = u; + } while (v != root); + + return weight; +} + +void get_related_weight(uint32_t x, H_partition* hap, double* w0, double* w1) +{ + (*w0) = (*w1) = 0; + if(x >= hap->link->a.n) return; + uint32_t i, a_n = hap->link->a.a[x].e.n; + hc_edge* a = hap->link->a.a[x].e.a; + for (i = 0; i < a_n; i++) + { + if(a[i].del) continue; + if((hap->hap[a[i].uID] & hap->m[0])) (*w0)+= a[i].weight; + if((hap->hap[a[i].uID] & hap->m[1])) (*w1)+= a[i].weight; + } + return; +} + +void set_path_hap(bub_p_t_warp *b, uint32_t root, H_partition* hap) +{ + uint32_t v, u, label; + double w0 = 0, w1 = 0, cur_w0, cur_w1; + ///v is the sink of this bubble + v = b->S.a[0]; + do { + u = b->a[v].p; // u->v + if(v != b->S.a[0]) + { + get_related_weight(v>>1, hap, &cur_w0, &cur_w1); + w0 += cur_w0; w1 += cur_w1; + } + v = u; + } while (v != root); + + if(w0 >= w1) + { + label = hap->label | hap->m[0]; + } + else + { + label = hap->label | hap->m[1]; + } + + v = b->S.a[0]; + do { + u = b->a[v].p; // u->v + if(v != b->S.a[0]) hap->hap[v>>1] |= label; + v = u; + } while (v != root); +} + + +uint64_t get_phase_path(ma_ug_t *ug, uint32_t s, uint32_t d, bub_p_t_warp *b, H_partition* hap) +{ + asg_t *g = ug->g; + if(g->seq[s>>1].del) return 0; // already deleted + if(get_real_length(g, s, NULL)<2) return 0; + uint32_t i, n_pending, is_first, to_replace, cur_nc, cur_uc, cur_ac, n_tips, tip_end, n_pop; + double cur_nh, cur_w0, cur_w1, cur_rate, max_rate, cur_weight, max_weight; + ///S saves nodes with all incoming edges visited + b->S.n = b->T.n = b->b.n = b->e.n = 0; + ///for each node, b->a saves all related information + b->a[s].d = b->a[s].nc = b->a[s].ac = b->a[s].uc = 0; b->a[s].nh = b->a[s].w[0] = b->a[s].w[1] = 0; + ///b->S is the nodes with all incoming edges visited + kv_push(uint32_t, b->S, s); + n_pop = n_tips = n_pending = 0; + tip_end = (uint32_t)-1; + is_first = 1; + + do { + ///v is a node that all incoming edges have been visited + ///d is the distance from v0 to v + uint32_t v = kv_pop(b->S); + uint32_t d = b->a[v].d, nc = b->a[v].nc, uc = b->a[v].uc, ac = b->a[v].ac; + double nh = b->a[v].nh; + double nw_0 = b->a[v].w[0], nw_1 = b->a[v].w[1]; + + uint32_t nv = asg_arc_n(g, v); + asg_arc_t *av = asg_arc_a(g, v); + for (i = 0; i < nv; ++i) { + uint32_t w = av[i].v, l = (uint32_t)av[i].ul; // v->w with length l, not overlap length + bub_p_t *t = &b->a[w]; + //got a circle + if ((w>>1) == (s>>1)) goto pop_reset; + //important when poping at long untig graph + if(is_first) l = 0; + if (av[i].del) continue; + ///push the edge + kv_push(uint32_t, b->e, (g->idx[v]>>32) + i); + + if (t->s == 0) + { // this vertex has never been visited + kv_push(uint32_t, b->b, w); // save it for revert + ///t->p is the parent node of + ///t->s = 1 means w has been visited + ///d is len(v0->v), l is len(v->w), so t->d is len(v0->w) + t->p = v, t->s = 1, t->d = d + l, t->nc = nc + ug->u.a[(w>>1)].n; + t->r = get_real_length(g, w^1, NULL); + /**need fix**/ + t->nh = nh + get_path_weight(w, v, s, b, hap->link); + + get_related_weight(w>>1, hap, &(t->w[0]), &(t->w[1])); + t->w[0] += nw_0; t->w[1] += nw_1; + + t->ac = ac + ((!is_hap_set(w>>1, *hap))?ug->u.a[(w>>1)].n : 0); + t->uc = uc + ((is_hap_set(w>>1, *hap))?ug->u.a[(w>>1)].n : 0); + + ++n_pending; + } + else { + to_replace = 0; + + cur_nc = nc + ug->u.a[(w>>1)].n; + /**need fix**/ + cur_nh = nh + get_path_weight(w, v, s, b, hap->link); + get_related_weight(w>>1, hap, &cur_w0, &cur_w1); + cur_w0 += nw_0; cur_w1 += nw_1; + cur_weight = cur_nh + MAX(cur_w0, cur_w1) - MIN(cur_w0, cur_w1); + max_weight = t->nh + MAX(t->w[0], t->w[1]) - MIN(t->w[0], t->w[1]); + + cur_ac = ac + ((!is_hap_set(w>>1, *hap))?ug->u.a[(w>>1)].n : 0);; + cur_uc = uc + ((is_hap_set(w>>1, *hap))?ug->u.a[(w>>1)].n : 0); + cur_rate = ((double)(cur_ac)/(double)(cur_ac+cur_uc)); + max_rate = ((double)(t->ac)/(double)(t->ac+t->uc)); + + if(cur_rate > max_rate) + { + to_replace = 1; + } + else if(cur_rate == max_rate) + { + ///if(cur_nh > t->nh) + if(cur_weight > max_weight) + { + to_replace = 1; + } + else if(cur_weight == max_weight)///(cur_nh == t->nh) + { + if(cur_nc > t->nc) + { + to_replace = 1; + } + else if(cur_nc == t->nc) + { + if(d + l > t->d) + { + to_replace = 1; + } + } + } + } + + + if(to_replace) + { + t->p = v; + t->nc = cur_nc; + t->nh = cur_nh; + t->ac = cur_ac; + t->uc = cur_uc; + t->w[0] = cur_w0; + t->w[1] = cur_w1; + } + + + if (d + l < t->d) t->d = d + l; // update dist + } + + if (--(t->r) == 0) { + uint32_t x = get_real_length(g, w, NULL); + if(x > 0) + { + kv_push(uint32_t, b->S, w); + } + else + { + ///at most one tip + if(n_tips != 0) goto pop_reset; + n_tips++; + tip_end = w; + } + --n_pending; + } + } + is_first = 0; + + + if(n_tips == 1) + { + if(tip_end != (uint32_t)-1 && n_pending == 0 && b->S.n == 0) + { + ///sink is b.S.a[0] + kv_push(uint32_t, b->S, tip_end); + break; + } + else + { + goto pop_reset; + } + } + + if (i < nv || b->S.n == 0) goto pop_reset; + }while (b->S.n > 1 || n_pending); + + + n_pop = 1; + /**need fix**/ + set_path_hap(b, s, hap); + pop_reset: + + for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices + bub_p_t *t = &b->a[b->b.a[i]]; + t->p = t->d = t->nc = t->ac = t->uc = t->r = t->s = 0; + t->nh = t->w[0] = t->w[1] = 0; + } + + return n_pop; +} + +uint32_t get_available_com(H_partition* hap, bubble_type* bub, ma_ug_t *ug, uint32_t check_self, uint32_t check_others) +{ + hc_links* link = hap->link; + uint32_t beg, sink, n, *a, i, j, k, uID, max_bub_i, max_non_bub_i, max_i, is_ava; + double w, max_bub_w, max_non_bub_w; + max_i = (uint32_t)-1; + + for (i = 0, max_bub_w = -1, max_bub_i = (uint32_t)-1; i < bub->num.n-1; i++) + { + get_bubbles(bub, i, &beg, &sink, &a, &n, NULL); + + for (j = 0, w = 0, is_ava = 0; j < n; j++) + { + uID = a[j]>>1; + if(check_self && is_hap_set(uID, *hap)) break; + for (k = 0; k < link->a.a[uID].e.n; k++) + { + if(link->a.a[uID].e.a[k].del) continue; + if(check_others && (!is_hap_set(link->a.a[uID].e.a[k].uID, *hap))) continue; + w += link->a.a[uID].e.a[k].weight; + is_ava = 1; + } + } + if(j != n) continue; + if(is_ava == 0) continue; + + if(w > max_bub_w) + { + max_bub_w = w; + max_bub_i = i; + } + } + + for (i = 0, max_non_bub_w = -1, max_non_bub_i = (uint32_t)-1; i < ug->u.n; i++) + { + if(bub->index[i] == bub->num.n) + { + uID = i; + if(check_self && is_hap_set(uID, *hap)) continue; + for (k = 0, w = 0, is_ava = 0; k < link->a.a[uID].e.n; k++) + { + if(link->a.a[uID].e.a[k].del) continue; + if(check_others && (!is_hap_set(link->a.a[uID].e.a[k].uID, *hap))) continue; + w += link->a.a[uID].e.a[k].weight; + is_ava = 1; + } + if(is_ava == 0) continue; + + if(w > max_non_bub_w) + { + max_non_bub_w = w; + max_non_bub_i = i; + } + } + } + + if(max_bub_i != (uint32_t)-1 && max_non_bub_i != (uint32_t)-1) + { + if(max_non_bub_w > max_bub_w) + { + max_i = (max_non_bub_i << 1) + 1; + w = max_non_bub_w; + } + else + { + max_i = (max_bub_i << 1); + w = max_bub_w; + } + } + else if(max_bub_i != (uint32_t)-1) + { + max_i = (max_bub_i << 1); + w = max_bub_w; + } + else if(max_non_bub_i != (uint32_t)-1) + { + max_i = (max_non_bub_i << 1) + 1; + w = max_non_bub_w; + } + + if(max_i == (uint32_t)-1) + { + for (i = 0, max_non_bub_w = -1, max_non_bub_i = (uint32_t)-1; i < ug->u.n; i++) + { + if(bub->index[i] < bub->num.n) + { + uID = i; + if(check_self && is_hap_set(uID, *hap)) continue; + for (k = 0, w = 0, is_ava = 0; k < link->a.a[uID].e.n; k++) + { + if(link->a.a[uID].e.a[k].del) continue; + if(check_others && (!is_hap_set(link->a.a[uID].e.a[k].uID, *hap))) continue; + w += link->a.a[uID].e.a[k].weight; + is_ava = 1; + } + if(is_ava == 0) continue; + + if(w > max_non_bub_w) + { + max_non_bub_w = w; + max_non_bub_i = i; + } + } + } + + if(max_non_bub_i != (uint32_t)-1) + { + max_i = (max_non_bub_i << 1) + 1; + w = max_non_bub_w; + } + } + + // if(max_i == (uint32_t)-1) + // { + // fprintf(stderr, "-Cannot find!\n"); + // } + // else if(max_i & 1) + // { + // fprintf(stderr, "-utg-%uth, phasing ID: %u, w: %f, max_bub_i: %u, max_bub_w: %f, max_non_bub_i: %u, max_non_bub_w: %f\n", + // max_i>>1, hap->label>>3, w, max_bub_i, max_bub_w, max_non_bub_i, max_non_bub_w); + // } + // else + // { + // fprintf(stderr, "-bubble-%uth, phasing ID: %u, w: %f, max_bub_i: %u, max_bub_w: %f, max_non_bub_i: %u, max_non_bub_w: %f\n", + // max_i>>1, hap->label>>3, w, max_bub_i, max_bub_w, max_non_bub_i, max_non_bub_w); + // } + + return max_i; +} + +uint32_t get_unset_com(H_partition* hap, bubble_type* bub, ma_ug_t *ug) +{ + uint32_t max_i = get_available_com(hap, bub, ug, 1, 1); + + if(max_i == (uint32_t)-1) + { + max_i = get_available_com(hap, bub, ug, 1, 0); + if(max_i != (uint32_t)-1) hap->label += hap->label_add; + } + + return max_i; +} + +void phase_com(H_partition* hap, ma_ug_t *ug, bub_p_t_warp* b, bubble_type* bub, uint32_t bid) +{ + + if((bid & 1) == 0) ///bubble + { + uint32_t beg = (uint32_t)-1, sink = (uint32_t)-1, n, *a; + get_bubbles(bub, bid>>1, &beg, &sink, &a, &n, NULL); + fprintf(stderr, "bubble-%uth, beg: %u, sink: %u, phasing ID: %u\n", bid>>1, beg>>1, sink>>1, hap->label>>3); + get_phase_path(ug, beg, sink, b, hap); + get_phase_path(ug, beg, sink, b, hap); + } + else + { + double cur_w0, cur_w1; + get_related_weight(bid>>1, hap, &cur_w0, &cur_w1); + fprintf(stderr, "utg-%uth, phasing ID: %u\n", bid>>1, hap->label>>3); + if(cur_w0 >= cur_w1) + { + hap->hap[bid>>1] |= (hap->label | hap->m[0]); + } + else + { + hap->hap[bid>>1] |= (hap->label | hap->m[1]); + } + } +} + +inline uint32_t get_phase_status(H_partition* hap, uint32_t uID) +{ + int d = -2; + if(hap->hap[uID] & hap->m[0]) d = 1; + if(hap->hap[uID] & hap->m[1]) d = -1; + if(hap->hap[uID] & hap->m[2]) d = 0; + return d; +} + +uint32_t init_contig_partition(H_partition* hap, ha_ug_index* idx, bubble_type* bub) +{ + hc_links* link = idx->link; + ma_ug_t *ug = idx->ug; + bub_p_t_warp b; + memset(&b, 0, sizeof(bub_p_t_warp)); + CALLOC(b.a, ug->g->n_seq*2); + uint32_t i, k, nv = ug->g->n_seq * 2; + int s_d, o_d; + for (i = 0; i < nv; i++) + { + b.a[i].w[0] = b.a[i].w[1] = b.a[i].nh = 0; + b.a[i].p =b.a[i].d = b.a[i].nc = b.a[i].uc = b.a[i].ac = b.a[i].r = b.a[i].s = 0; + } + + + hap->n = ug->u.n; + MALLOC(hap->hap, hap->n); + memset(hap->hap, 0, hap->n*sizeof(uint32_t)); + MALLOC(hap->lock, hap->n); + memset(hap->lock, 0, hap->n); + MALLOC(hap->weight, hap->n); + hap->m[0] = 1; hap->m[1] = 2; hap->m[2] = 4; + hap->link = link; + hap->label = 0; + hap->label_add = 8; + + uint32_t max_i = get_available_com(hap, bub, ug, 0, 0); + + if(max_i == (uint32_t)-1) return 0; + + hap->label = 0; + while (1) + { + max_i = get_unset_com(hap, bub, ug); + if(max_i == (uint32_t)-1) break; + phase_com(hap, ug, &b, bub, max_i); + } + + free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); + + s_d = 0; + for (i = 0; i < hap->n; i++) + { + if(bub->index[i] > bub->num.n) continue; + s_d = get_phase_status(hap, i); + ///if(s_d == -2 && link->a.a[i].e.n > 0) fprintf(stderr, "ERROR\n"); + if(s_d == -2) continue; + hap->weight[i] = 0; + + for (k = 0; k < link->a.a[i].e.n; k++) + { + if(link->a.a[i].e.a[k].del) continue; + o_d = get_phase_status(hap, link->a.a[i].e.a[k].uID); + ///if(o_d == -2) fprintf(stderr, "ERROR\n"); + hap->weight[i] += (s_d*o_d*link->a.a[i].e.a[k].weight); + } + } + + return 1; +} + +uint32_t get_max_unitig(H_partition* hap, hc_links* link, ma_ug_t *ug, bubble_type* bub) +{ + uint32_t i; + for (i = 0; i < hap->n; i++) + { + + } + + return 0; +} + +uint32_t phasing_improevment(H_partition* hap, ha_ug_index* idx, bubble_type* bub) +{ + ; + + return 0; +} + +void destory_contig_partition(H_partition* hap) +{ + free(hap->lock); + free(hap->hap); + free(hap->weight); } int hic_short_align(const char *fn1, const char *fn2, ha_ug_index* idx) { double index_time = yak_realtime(); - hc_links link; sldat_t sl; gzFile fp1, fp2; if ((fp1 = gzopen(fn1, "r")) == 0) return 0; @@ -3956,11 +4908,9 @@ int hic_short_align(const char *fn1, const char *fn2, ha_ug_index* idx) bubble_type bub; identify_bubbles(idx->ug, &bub); - init_hc_links(&link, sl.idx->ug->g->n_seq); - collect_hc_reverse_links(link, idx->ug, bub); + ///print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, idx->link, idx); - - if(!load_hc_links(&link, &sl.hits, asm_opt.output_file_name)) + if(!load_hc_hits(&sl.hits, asm_opt.output_file_name)) { /*******************************for debug************************************/ // load_reads(&R1, fn1); @@ -3975,36 +4925,41 @@ int hic_short_align(const char *fn1, const char *fn2, ha_ug_index* idx) // print_hits(idx, &sl.hits, fn1); /*******************************for debug************************************/ dedup_hits(&sl.hits); - - collect_hc_links(sl.idx, &sl.hits, &link, &bub); - write_hc_links(&link, &sl.hits, asm_opt.output_file_name); + write_hc_hits(&sl.hits, asm_opt.output_file_name); } + collect_hc_links(sl.idx, &sl.hits, idx->link, &bub); + collect_hc_reverse_links(idx->link, idx->ug, &bub); + + ///debug_hc_links(idx, idx->link, &sl, &bub, fn1); + + init_hic_p((ha_ug_index*)sl.idx, &sl.hits, idx->link, &bub); + ///print_hc_links(idx->link, 0); + H_partition hap; + init_contig_partition(&hap, idx, &bub); + destory_contig_partition(&hap); - debug_hc_links(idx, &link, &sl, &bub, fn1); - init_hic_p(sl.idx, &sl.hits, &link, &bub); return 1; /*******************************for debug************************************/ // destory_reads(&R1); // destory_reads(&R2); /*******************************for debug************************************/ - print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, &link, idx); - collect_hc_reverse_links(&link, idx->ug, &bub); - normalize_hc_links(&link); + print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, idx->link, idx); + collect_hc_reverse_links(idx->link, idx->ug, &bub); + normalize_hc_links(idx->link); /*******************************for debug************************************/ ///print_hc_links(&link); /*******************************for debug************************************/ - min_cut_t* cut = clean_hap(&link, &bub, idx->ug); + min_cut_t* cut = clean_hap(idx->link, &bub, idx->ug); ///print_bubbles(idx->ug, &bub, NULL, &link, idx); - G_partition* gp = clean_bubbles(&link, &bub, cut, idx->ug); + G_partition* gp = clean_bubbles(idx->link, &bub, cut, idx->ug); ///print_hc_links(&link); destory_min_cut_t(cut); free(cut); destory_G_partition(gp); free(gp); kv_destroy(sl.hits.a); - destory_hc_links(&link); destory_bubbles(&bub); kseq_destroy(sl.ks1); kseq_destroy(sl.ks2); @@ -4015,16 +4970,343 @@ int hic_short_align(const char *fn1, const char *fn2, ha_ug_index* idx) return 1; } -void hic_analysis(ma_ug_t *ug) + +void hic_analysis(ma_ug_t *ug, asg_t* read_g, hc_links* link) { - ug_index = NULL; int exist = load_hc_pt_index(&ug_index, asm_opt.output_file_name); if(exist == 0) ug_index = build_unitig_index(ug, asm_opt.hic_mer_length); if(exist == 0) write_hc_pt_index(ug_index, asm_opt.output_file_name); ug_index->ug = ug; + ug_index->read_g = read_g; + ug_index->link = link; ///test_unitig_index(ug_index, ug); hic_short_align(asm_opt.hic_files[0], asm_opt.hic_files[1], ug_index); destory_hc_pt_index(ug_index); +} + +typedef struct{ + //[uID_start, uID_end) + uint64_t uID_start; + uint64_t uID_end; + uint64_t u_n; + uint64_t r_n; + uint64_t* r_idx; +} bench_utg; + +typedef struct{ + uint64_t s, e; +}homo_interval; + +typedef struct{ + kvec_t(bench_utg) ug_idx; + uint64_t uID_bits; + uint64_t pos_mode; + hc_links link; + kvec_t(homo_interval) regions; +}bench_idx; + +uint64_t* set_bench_idx(ma_ug_t *ug, asg_t* read_g, uint64_t uID_start, uint64_t uID_end, uint64_t uID_bits, uint64_t r_n) +{ + uint64_t *idx = (uint64_t*)malloc(sizeof(uint64_t)*r_n), i, k; + memset(idx, -1, sizeof(uint64_t)*r_n); + uint64_t rId, ori, start, l; + ma_utg_t *u = NULL; + for (i = uID_start; i < uID_end; i++) + { + u = &(ug->u.a[i]); + if(u->n == 0) continue; + for (k = l = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + ori = u->a[k]>>32&1; + start = l; + l += (uint32_t)u->a[k]; + if(idx[rId] != (uint64_t)-1) + { + idx[rId] = (uint64_t)-1; + } + else + { + idx[rId] = (ori<<63) + ((i<<(64-uID_bits))>>1) + start; + if(ori) idx[rId] = idx[rId] + read_g->seq[rId].len - 1; + } + } + } + + return idx; +} + +void get_r_utg_bench(uint64_t index, bench_idx* idx, ma_ug_t *ug) +{ + + bench_utg* a_list = idx->ug_idx.a; + uint64_t a_n = idx->ug_idx.n; + bench_utg *x = &(a_list[index]), *y = NULL; + uint64_t i, k, t, rev, x_uid, y_uid, y_pos, x_pos, d; + uint64_t rId, ori; + ma_utg_t *u = NULL; + for (i = x->uID_start; i < x->uID_end; i++) + { + u = &(ug->u.a[i]); + x_uid = i; + if(u->n == 0) continue; + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + ori = u->a[k]>>32&1; + if(x->r_idx[rId] == (uint64_t)-1) continue; + x_pos = x->r_idx[rId] & idx->pos_mode; + + for (t = 0; t < a_n; t++) + { + if(t == index) continue; + y = &(a_list[t]); + if(y->r_idx[rId] == (uint64_t)-1) continue; + rev = 0; + if((y->r_idx[rId]>>63) != ori) rev = 1; + y_uid = (y->r_idx[rId]<<1)>>(64 - idx->uID_bits); + y_pos = y->r_idx[rId] & idx->pos_mode; + if(rev) y_pos = ug->u.a[y_uid].len - y_pos - 1; + ///if(ori) x_pos = ug->u.a[x_uid].len - x_pos - 1, y_pos = ug->u.a[y_uid].len - y_pos - 1; + d = MAX(x_pos, y_pos) - MIN(x_pos, y_pos); + d = (d<<2) + (rev<<1); + if(y_pos > x_pos) d = d + 1; + push_hc_edge(&(idx->link.a.a[x_uid]), y_uid, 1, 0, &d); + if(x_pos != y_pos) d = d ^ 1; + push_hc_edge(&(idx->link.a.a[y_uid]), x_uid, 1, 0, &d); + } + } + } +} + +void hap_ID(bench_idx* idx, uint64_t ID, uint64_t* hapID, uint64_t* uID) +{ + uint64_t i; + (*hapID) = (*uID) = (uint64_t)-1; + for (i = 0; i < idx->ug_idx.n; i++) + { + if(ID >= idx->ug_idx.a[i].uID_start && ID < idx->ug_idx.a[i].uID_end) + { + (*hapID) = i; + (*uID) = ID - idx->ug_idx.a[i].uID_start; + return; + } + } + return; +} + +void print_bench_idx(bench_idx* idx, ma_ug_t *ug) +{ + uint64_t i, k, s_uID, s_hapID, d_uID, d_hapID; + long long x[2] = {1, -1}; + for (i = 0; i < idx->link.a.n; i++) + { + for (k = 0; k < idx->link.a.a[i].e.n; k++) + { + if(idx->link.a.a[i].e.a[k].del) continue; + hap_ID(idx, i, &s_hapID, &s_uID); + hap_ID(idx, idx->link.a.a[i].e.a[k].uID, &d_hapID, &d_uID); + fprintf(stderr, "s-hap%lu-utg%.6d\td-hap%lu-utg%.6d\t%c\t%lld\n", + s_hapID, (int)(s_uID+1), d_hapID, (int)(d_uID+1), + "+-"[!!(idx->link.a.a[i].e.a[k].dis&(uint64_t)2)], + ((long long)(idx->link.a.a[i].e.a[k].dis>>2))*x[idx->link.a.a[i].e.a[k].dis&(uint64_t)1]); + } + } + + +} + +uint64_t get_hic_distance_bench(pe_hit* hit, hc_links* link, bench_idx* idx, ma_ug_t *ug, uint64_t* is_trans) +{ + (*is_trans) = (uint64_t)-1; + uint64_t s_uid, e_uid; + long long s_pos, e_pos; + s_uid = ((hit->s<<1)>>(64 - idx->uID_bits)); s_pos = hit->s & idx->pos_mode; + e_uid = ((hit->e<<1)>>(64 - idx->uID_bits)); e_pos = hit->e & idx->pos_mode; + if(s_uid == e_uid) + { + (*is_trans) = 0; + return MAX(s_pos, e_pos) - MIN(s_pos, e_pos); + } + + uint64_t s_i, e_i, k, ori; + for (s_i = 0; s_i < idx->ug_idx.n; s_i++) + { + if(s_uid >= idx->ug_idx.a[s_i].uID_start && s_uid < idx->ug_idx.a[s_i].uID_end) break; + } + for (e_i = 0; e_i < idx->ug_idx.n; e_i++) + { + if(e_uid >= idx->ug_idx.a[e_i].uID_start && e_uid < idx->ug_idx.a[e_i].uID_end) break; + } + if(s_i == idx->ug_idx.n || e_i == idx->ug_idx.n) return (uint64_t)-1; + if(s_i == e_i) + { + (*is_trans) = 0; + return (uint64_t)-1; + } + + (*is_trans) = 1; + hc_linkeage* t = &(link->a.a[s_uid]); + long long m_x[2] = {1, -1}, dis; + for (k = 0; k < t->e.n; k++) + { + if(t->e.a[k].del || t->e.a[k].uID != e_uid) continue; + ori = !!(t->e.a[k].dis & (uint64_t)2); + dis = (long long)(t->e.a[k].dis>>2) * m_x[t->e.a[k].dis & (uint64_t)1]; + if(ori) e_pos = ug->u.a[e_uid].len - e_pos - 1; + e_pos = e_pos + dis; + return MAX(s_pos, e_pos) - MIN(s_pos, e_pos); + } + + return (uint64_t)-1; +} + +void init_bench_idx(bench_idx* idx, asg_t* read_g, ma_ug_t *ug) +{ + uint64_t i, occ; + kv_init(idx->ug_idx); + kv_init(idx->regions); + kv_malloc(idx->ug_idx, ug->occ.n); idx->ug_idx.n = ug->occ.n; + for (idx->uID_bits = 1; (uint64_t)(1<uID_bits)<(uint64_t)ug->u.n; idx->uID_bits++); + idx->pos_mode = ((uint64_t)-1)>>(idx->uID_bits+1); + for (i = occ = 0; i < ug->occ.n; i++) + { + idx->ug_idx.a[i].uID_start = occ; + occ += ug->occ.a[i]; + idx->ug_idx.a[i].uID_end = occ; + idx->ug_idx.a[i].u_n = ug->occ.a[i]; + + idx->ug_idx.a[i].r_n = read_g->n_seq; + idx->ug_idx.a[i].r_idx + = set_bench_idx(ug, read_g, idx->ug_idx.a[i].uID_start, idx->ug_idx.a[i].uID_end, + idx->uID_bits, idx->ug_idx.a[i].r_n); + } + + init_hc_links(&(idx->link), ug->u.n, ug->g->n_seq); + + for (i = 0; i < idx->ug_idx.n; i++) + { + get_r_utg_bench(i, idx, ug); + } +} + +void evaluate_bench_idx(bench_idx* idx, kvec_pe_hit* hits, ma_ug_t *ug) +{ + uint64_t k, distance, is_trans, trans[2]; + kvec_t(uint64_t) buf; + kv_init(buf); + for (k = trans[0] = trans[1] = 0; k < hits->a.n; ++k) + { + distance = get_hic_distance_bench(&(hits->a.a[k]), &(idx->link), idx, ug, &is_trans); + if(is_trans != (uint64_t)-1) trans[is_trans]++; + if(distance == (uint64_t)-1 || is_trans == (uint64_t)-1) continue; + distance = (distance << 1) + is_trans; + kv_push(uint64_t, buf, distance); + } + + radix_sort_hc64(buf.a, buf.a+buf.n); + + for (k = 0; k < buf.n; k++) + { + fprintf(stderr, "%lu\t%lu\n", buf.a[k]>>1, buf.a[k]&1); + } + + // uint64_t up_dis = buf.a[(uint64_t)(buf.n*0.99)]>>1, step = 1000; + // uint64_t step_s = 0, step_e = step, cnt[2]; + // for (k = cnt[0] = cnt[1] = 0; k < buf.n; k++) + // { + // if(step_s > up_dis) step_e = (buf.a[buf.n-1]>>1) + 1; + // if((buf.a[k]>>1) < step_e && (buf.a[k]>>1) >= step_s) + // { + // cnt[buf.a[k]&1]++; + // } + // if((buf.a[k]>>1) >= step_e) + // { + // while (!((buf.a[k]>>1) < step_e && (buf.a[k]>>1) >= step_s)) + // { + // fprintf(stderr, "i: %lu, step_s: %lu, step_e: %lu, cnt[0]: %lu, cnt[1]: %lu, rate: %f\n", + // step_s/step, step_s, step_e, cnt[0], cnt[1], ((double)cnt[1])/(double)(cnt[0] + cnt[1])); + // step_s += step; + // step_e += step; + // cnt[0] = cnt[1] = 0; + // } + // } + // } + + // if(cnt[0] > 0 || cnt[1] > 0) + // { + // fprintf(stderr, "i: %lu, step_s: %lu, step_e: %lu, cnt[0]: %lu, cnt[1]: %lu, rate: %f\n", + // step_s/step, step_s, step_e, cnt[0], cnt[1], ((double)cnt[1])/(double)(cnt[0] + cnt[1])); + // } + + kv_destroy(buf); +} + +void destory_bench_idx(bench_idx* idx) +{ + uint64_t i; + for (i = 0; i < idx->ug_idx.n; i++) + { + free(idx->ug_idx.a[i].r_idx); + } + kv_destroy(idx->ug_idx); + kv_destroy(idx->regions); + destory_hc_links(&(idx->link)); +} + + +int hic_short_align_bench(const char *fn1, const char *fn2, const char *output_file_name, ha_ug_index* idx) +{ + double index_time = yak_realtime(); + sldat_t sl; + gzFile fp1, fp2; + if ((fp1 = gzopen(fn1, "r")) == 0) return 0; + if ((fp2 = gzopen(fn2, "r")) == 0) return 0; + sl.ks1 = kseq_init(fp1); + sl.ks2 = kseq_init(fp2); + sl.idx = idx; + sl.chunk_size = 20000000; + sl.n_thread = asm_opt.thread_num; + sl.total_base = sl.total_pair = 0; + idx->max_cnt = 5; + kv_init(sl.hits.a); + fprintf(stderr, "u.n: %d, uID_bits: %lu, pos_bits: %lu\n", (uint32_t)idx->ug->u.n, idx->uID_bits, idx->pos_bits); + + if(!load_hc_hits(&sl.hits, output_file_name)) + { + kt_pipeline(3, worker_pipeline, &sl, 3); + dedup_hits(&sl.hits); + write_hc_hits(&sl.hits, output_file_name); + } + bench_idx bench; + init_bench_idx(&bench, idx->read_g, idx->ug); + ///print_bench_idx(&bench, idx->ug); + evaluate_bench_idx(&bench, &sl.hits, idx->ug); + + destory_bench_idx(&bench); + kv_destroy(sl.hits.a); + kseq_destroy(sl.ks1); + kseq_destroy(sl.ks2); + gzclose(fp1); + gzclose(fp2); + fprintf(stderr, "[M::%s::%.3f] processed %lu pairs; %lu bases\n", __func__, yak_realtime()-index_time, sl.total_pair, sl.total_base); + return 1; +} + +void hic_benchmark(ma_ug_t *ug, asg_t* read_g) +{ + char *output_file_name = (char*)calloc(strlen(asm_opt.output_file_name) + 25, 1); + sprintf(output_file_name, "%s.bench", asm_opt.output_file_name); + ug_index = NULL; + int exist = load_hc_pt_index(&ug_index, output_file_name); + if(exist == 0) ug_index = build_unitig_index(ug, asm_opt.hic_mer_length); + if(exist == 0) write_hc_pt_index(ug_index, output_file_name); + ug_index->ug = ug; + ug_index->read_g = read_g; + + hic_short_align_bench(asm_opt.hic_files[0], asm_opt.hic_files[1], output_file_name, ug_index); + + free(output_file_name); } \ No newline at end of file diff --git a/hic.h b/hic.h index 12daefb..98662ef 100644 --- a/hic.h +++ b/hic.h @@ -3,7 +3,12 @@ #include #include "Overlaps.h" +#define kdq_clear(q) ((q)->count = (q)->front = 0) +#define kv_malloc(v, s) ((v).n = 0, (v).m = (s), MALLOC((v).a, (s))) + +hc_edge* get_hc_edge(hc_links* link, uint64_t src, uint64_t dest, uint64_t dir); void push_hc_edge(hc_linkeage* x, uint64_t uID, int weight, int dir, uint64_t* d); -void hic_analysis(ma_ug_t *ug); +void hic_analysis(ma_ug_t *ug, asg_t* read_g, hc_links* link); +void hic_benchmark(ma_ug_t *ug, asg_t* read_g); #endif