From 39c19618e8e17558da4898f3b5f4876f00cd4d0e Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Sun, 18 Apr 2021 21:57:16 -0400 Subject: [PATCH] backup hic --- Overlaps.cpp | 63 ++++++++++++++++++++++++++++++---------------------- hic.cpp | 51 +++++++++++++++++++++++++++++++++++++----- 2 files changed, 81 insertions(+), 33 deletions(-) diff --git a/Overlaps.cpp b/Overlaps.cpp index 413ec9e..542a888 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -12578,7 +12578,7 @@ trans_chain* load_hc_hits(const char *fn) fclose(fp); free(buf); - // fprintf(stderr, "[M::%s::] ==> Hi-C cov have been loaded\n", __func__); + fprintf(stderr, "[M::%s::] ==> Hi-C cov have been loaded\n", __func__); return t_ch; } @@ -12634,14 +12634,14 @@ bub_label_t* b_mask_t) hic_analysis(ug, sg, cov?cov->t_ch:t_ch); - /** - char* gfa_name = (char*)malloc(strlen(output_file_name)+25); - sprintf(gfa_name, "%s.d_utg.noseq.gfa", output_file_name); - FILE* output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); - fclose(output_file); - free(gfa_name); - **/ + + // char* gfa_name = (char*)malloc(strlen(output_file_name)+25); + // sprintf(gfa_name, "%s.d_utg.noseq.gfa", output_file_name); + // FILE* output_file = fopen(gfa_name, "w"); + // ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); + // fclose(output_file); + // free(gfa_name); + if(cov) destory_hap_cov_t(&cov); if(t_ch) destory_trans_chain(&t_ch); @@ -16210,6 +16210,8 @@ void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sou fprintf(stderr, "bin files have been written.\n"); } +int load_debug_graph(asg_t** sg, ma_hit_t_alloc** sources, ma_sub_t** coverage_cut, +char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex); int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, char* output_file_name) { char* gfa_name = (char*)malloc(strlen(output_file_name)+25); @@ -16218,6 +16220,14 @@ int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_s free(gfa_name); return 0; } + if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_debug_graph(NULL, NULL, NULL, output_file_name, NULL, NULL)) + { + (*sources) = NULL; + (*reverse_sources) = NULL; + free(gfa_name); + return 1; + } + sprintf(gfa_name, "%s.ovlp.source", output_file_name); if (!load_ma_hit_ts(sources, gfa_name)) { free(gfa_name); @@ -22858,8 +22868,7 @@ int write_ruIndex(R_to_U* ruIndex, char* read_file_name) fwrite(&ruIndex->len, sizeof(ruIndex->len), 1, fp); fwrite(ruIndex->index, sizeof(ruIndex->index[0]), ruIndex->len, fp); fwrite(R_INF.trio_flag, sizeof(R_INF.trio_flag[0]), ruIndex->len, fp); - fwrite(ruIndex->is_het, 1, ruIndex->len, fp); - + // fwrite(ruIndex->is_het, 1, ruIndex->len, fp); free(index_name); fflush(fp); fclose(fp); @@ -22884,8 +22893,8 @@ int load_ruIndex(R_to_U* ruIndex, char* read_file_name) R_INF.trio_flag = (uint8_t*)malloc(sizeof(uint8_t)*(ruIndex)->len); f_flag += fread(R_INF.trio_flag, sizeof(R_INF.trio_flag[0]), (ruIndex)->len, fp); - CALLOC(ruIndex->is_het, ruIndex->len); - f_flag += fread(ruIndex->is_het, 1, ruIndex->len, fp); + // CALLOC(ruIndex->is_het, ruIndex->len); + // f_flag += fread(ruIndex->is_het, 1, ruIndex->len, fp); free(index_name); fflush(fp); @@ -23099,7 +23108,11 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex) fp = fopen(gfa_name, "r"); if(!fp) return 0; sprintf(gfa_name, "%s.all.debug.asg_t.bin", output_file_name); fp = fopen(gfa_name, "r"); if(!fp) return 0; - + if((sg == NULL) || (sources == NULL) || (coverage_cut == NULL) || (reverse_sources == NULL) || + (ruIndex == NULL)) + { + return 1; + } if((*sources)!=NULL) { @@ -29791,7 +29804,7 @@ float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr, ma_sub_t **coverage_cut_ptr, int debug_g) { - char *o_file = NULL; + char *o_file = get_outfile_name(output_file_name); ma_sub_t *coverage_cut = *coverage_cut_ptr; asg_t *sg = *sg_ptr; bub_label_t b_mask_t; @@ -29982,24 +29995,20 @@ ma_sub_t **coverage_cut_ptr, int debug_g) ruIndex->is_het = rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz, &b_mask_t); - if (asm_opt.flag & HA_F_VERBOSE_GFA) - { - /*******************************for debug***************************************/ - write_debug_graph(sg, sources, coverage_cut, output_file_name, n_read, reverse_sources, ruIndex); - debug_gfa:; - /*******************************for debug***************************************/ - set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, - max_hang_length, mini_overlap_length); - } - - - o_file = get_outfile_name(output_file_name); + output_unitig_graph(sg, coverage_cut, o_file, sources, ruIndex, max_hang_length, mini_overlap_length); flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL; output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources, asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length); + if (asm_opt.flag & HA_F_VERBOSE_GFA) + { + write_debug_graph(sg, sources, coverage_cut, output_file_name, n_read, reverse_sources, ruIndex); + debug_gfa:; + set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length); + } + if (ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt)) { if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION; diff --git a/hic.cpp b/hic.cpp index a5afba0..bebdb69 100644 --- a/hic.cpp +++ b/hic.cpp @@ -13514,6 +13514,39 @@ void idx_hc_links(kvec_pe_hit* hits, ha_ug_index* idx, bubble_type* bub) } } +inline uint32_t trans_checking_pass(bubble_type* bub, kv_u_trans_t *ref, uint32_t x, uint32_t y) +{ + if(u_trans_n(*ref, x) == 0 || u_trans_n(*ref, y) == 0) return 0; + u_trans_t *a = NULL; + uint32_t n, i, f[2], qn, tn; + + qn = x; tn = y; + a = u_trans_a(*ref, qn); n = u_trans_n(*ref, qn); + for (i = 0, f[0] = f[1] = 0; i < n; i++) + { + if(a[i].del) continue; + if(a[i].f == RC_2) continue; + if(IF_HOM(a[i].tn, *bub)) continue; + f[(a[i].tn == tn && a[i].f != RC_2)]++; + } + if(f[0] != 0) return 0; + if(f[1] == 0) return 0; + + qn = y; tn = x; + a = u_trans_a(*ref, qn); n = u_trans_n(*ref, qn); + for (i = 0, f[0] = f[1] = 0; i < n; i++) + { + if(a[i].del) continue; + if(IF_HOM(a[i].tn, *bub)) continue; + f[(a[i].tn == tn && a[i].f != RC_2)]++; + } + if(f[0] != 0) return 0; + if(f[1] == 0) return 0; + + fprintf(stderr, "M::%s::s-utg%.6ul<----->d-utg%.6ul\n", __func__, x+1, y+1); + return 1; +} + void weight_kv_u_trans(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, kv_u_trans_t *ta, trans_idx* dis) { @@ -13653,7 +13686,7 @@ pe_hit *hits, uint32_t occ, uint32_t qid, uint32_t qs, uint32_t qe, uint32_t tid return w; } -double get_hits_weight(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, trans_idx* dis, +double get_hits_weight(ha_ug_index* idx, bubble_type* bub, kvec_pe_hit* hits, hc_links* link, trans_idx* dis, u_trans_t *t_a, uint32_t t_n, uint32_t qid, kv_u_trans_t *ta_idx) { /****************************may have bugs********************************/ @@ -13666,6 +13699,8 @@ u_trans_t *t_a, uint32_t t_n, uint32_t qid, kv_u_trans_t *ta_idx) double w, i_w; for (k = 0, w = 0; k < t_n; k++) { + if(IF_HOM(qid, *bub)) continue; + if(IF_HOM(t_a[k].tn, *bub)) continue; /****************************may have bugs********************************/ if(qid != t_a[k].tn) { @@ -13705,6 +13740,7 @@ u_trans_t *t_a, uint32_t t_n, uint32_t qid, kv_u_trans_t *ta_idx) return w; } + void adjust_weight_kv_u_trans(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, kv_u_trans_t *ta, kv_u_trans_t *ref, trans_idx* dis) { @@ -13723,12 +13759,15 @@ kv_u_trans_t *ta, kv_u_trans_t *ref, trans_idx* dis) { ta->a[m].nw /= (double)(MIN(hits->occ.a[ta->a[m].qn], hits->occ.a[ta->a[m].tn])); } + /*******************************for debug************************************/ + trans_checking_pass(bub, ref, ta->a[i].qn, ta->a[i].tn); + // if(trans_checking_pass(bub, ref, ta->a[i].qn, ta->a[i].tn)) ta->a[m].nw = 0; + /*******************************for debug************************************/ m++; } ta->n = m; - // fprintf(stderr, "+++++ta->n=%u\n", (uint32_t)ta->n); for (k = 0; k < ta->idx.n; k++)///all nodes { if(IF_HOM(k, *bub)) continue; @@ -13743,10 +13782,10 @@ kv_u_trans_t *ta, kv_u_trans_t *ref, trans_idx* dis) if(IF_HOM(a[i].qn, *bub)) continue; if(IF_HOM(a[i].tn, *bub)) continue; ///(qn, tn^) - a[i].nw += get_hits_weight(idx, hits, link, dis, + a[i].nw += get_hits_weight(idx, bub, hits, link, dis, u_trans_a(*ref, a[i].tn), u_trans_n(*ref, a[i].tn), a[i].qn, ta); ///(qn^1, tn) - a[i].nw += get_hits_weight(idx, hits, link, dis, + a[i].nw += get_hits_weight(idx, bub, hits, link, dis, u_trans_a(*ref, a[i].qn), u_trans_n(*ref, a[i].qn), a[i].tn, ta); } @@ -13759,10 +13798,10 @@ kv_u_trans_t *ta, kv_u_trans_t *ref, trans_idx* dis) if(vis[tn]) continue; if(qn == tn) continue; ///(qn, tn^) - w += get_hits_weight(idx, hits, link, dis, + w += get_hits_weight(idx, bub, hits, link, dis, u_trans_a(*ref, tn), u_trans_n(*ref, tn), qn, ta); ///(qn^1, tn) - w += get_hits_weight(idx, hits, link, dis, + w += get_hits_weight(idx, bub, hits, link, dis, u_trans_a(*ref, qn), u_trans_n(*ref, qn), tn, ta); if(w == 0) continue; kv_pushp(u_trans_t, *ta, &p);