diff --git a/Assembly.cpp b/Assembly.cpp index bc719ea..07c82b4 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1775,6 +1775,7 @@ int ha_assemble(void) { // debug_mc_g_t(MC_NAME); // debug_mc_gg_t(MC_NAME, 0, 0); + // quick_debug_phasing(MC_NAME); extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o); int r, hom_cov = -1, ovlp_loaded = 0; if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) { diff --git a/CommandLines.cpp b/CommandLines.cpp index 821cb9a..dd7f771 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -143,6 +143,12 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " --l-msjoin INT\n"); fprintf(stderr, " detect misjoined unitigs of >=INT in size; 0 to disable [%lu]\n", asm_opt->misjoin_len); + fprintf(stderr, " Ultra-Long-integration (beta):\n"); + fprintf(stderr, " --ul FILEs file names of Ultra-Long reads [r1.fq,r2.fq,...]\n"); + fprintf(stderr, " --ul-rate FLOAT\n"); + fprintf(stderr, " similarity threshold for UL-to-HiFi alignment [%.3g]\n", asm_opt->ul_error_rate); + + fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n"); fprintf(stderr, "See `https://hifiasm.readthedocs.io/en/latest/' or `man ./hifiasm.1' for complete documentation.\n"); } diff --git a/CommandLines.h b/CommandLines.h index 1fd6948..cf237b0 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.17.2-r433" +#define HA_VERSION "0.17.3-r439" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index c1911cb..687a4fa 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -817,7 +817,7 @@ long long get_specific_overlap(ma_hit_t_alloc* x, uint32_t qn, uint32_t tn) -inline void set_reverse_overlap(ma_hit_t* dest, ma_hit_t* source) +void set_reverse_overlap(ma_hit_t* dest, ma_hit_t* source) { dest->qns = Get_tn(*source); dest->qns = dest->qns << 32; @@ -924,7 +924,7 @@ void normalize_ma_hit_t_single_side_advance(ma_hit_t_alloc* sources, long long n } } - // if(VERBOSE >= 1) + if(VERBOSE >= 1) { fprintf(stderr, "[M::%s] takes %0.2fs\n\n", __func__, Get_T()-startTime); } @@ -10140,7 +10140,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FIL v = au[j].v; fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\tL2:i:%u\n", prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], - prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), au[j].ou); + prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), 0/**au[j].ou**/); } @@ -10153,7 +10153,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FIL v = au[j].v; fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\tL2:i:%u\n", prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], - prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), au[j].ou); + prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j]), 0/**au[j].ou**/); } } } @@ -13674,14 +13674,14 @@ long long gap_fuzz, bub_label_t* b_mask_t) if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_trans_chain(cov->t_ch, output_file_name); } - + ///for debug // char* gfa_name = (char*)malloc(strlen(output_file_name)+50); // sprintf(gfa_name, "%s.pre.clean_d_utg.noseq.gfa", output_file_name); // FILE* output_file = fopen(gfa_name, "w"); // ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file); // fclose(output_file); // free(gfa_name); - + ///for debug hic_analysis(ug, sg, cov?cov->t_ch:t_ch, &opt, 0, asm_opt.scffold?&rhits:NULL); @@ -16641,52 +16641,79 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov) redo: ///print_untig((ug), 61955, "i-0:", 0); + + // fprintf(stderr, "[M::%s] 0\n", __func__); asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov, NULL, 1); + // fprintf(stderr, "[M::%s] 1\n", __func__); magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate); + // fprintf(stderr, "[M::%s] 2\n", __func__); /**********debug**********/ if(just_bubble_pop == 0) { cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2); } + // fprintf(stderr, "[M::%s] 3\n", __func__); /**********debug**********/ long long pre_cons = get_graph_statistic(g); long long cur_cons = 0; while(pre_cons != cur_cons) { + // fprintf(stderr, "[M::%s] 4\n", __func__); pre_cons = get_graph_statistic(g); + // fprintf(stderr, "[M::%s] 5\n", __func__); ///need consider tangles asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov, NULL, 1); + // fprintf(stderr, "[M::%s] 6\n", __func__); /**********debug**********/ 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, trio_flag, cov, NULL); + // fprintf(stderr, "[M::%s] 7\n", __func__); // if(trio_flag == MOTHER) print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, cov, NULL); + // fprintf(stderr, "[M::%s] 8\n", __func__); asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL, trio_flag); + // fprintf(stderr, "[M::%s] 9\n", __func__); asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL, trio_flag); + // fprintf(stderr, "[M::%s] 10\n", __func__); detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL, cov->is_r_het); + // fprintf(stderr, "[M::%s] 11\n", __func__); ///need consider tangles ///note we need both the read graph and the untig graph } /**********debug**********/ cur_cons = get_graph_statistic(g); + // fprintf(stderr, "[M::%s] 12\n", __func__); } if(just_bubble_pop == 0) { + // fprintf(stderr, "[M::%s] 13\n", __func__); cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2); + // fprintf(stderr, "[M::%s] 14\n", __func__); } // print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); + // fprintf(stderr, "[M::%s] 15\n", __func__); magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate); + // fprintf(stderr, "[M::%s] 16\n", __func__); + + + // print_debug_gfa(read_g, ug, coverage_cut, "resolve_tangles", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len, 0, 0, 0); + // exit(1); + ///bug here resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, cov->is_r_het, trio_flag, drop_ratio); + // fprintf(stderr, "[M::%s] 17\n", __func__); drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het); + // fprintf(stderr, "[M::%s] 18\n", __func__); all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, cov->is_r_het, DOUBLE_CHECK_THRES, asm_opt.trio_flag_occ_thres); + // fprintf(stderr, "[M::%s] 19\n", __func__); // if(trio_flag == MOTHER) print_untig_by_read(ug, "m54329U_190827_173812/30214441/ccs", (uint32_t)-1, NULL, NULL, "bf-16"); if(is_first) { is_first = 0; unitig_arc_del_short_diploid_by_length(ug->g, drop_ratio); + // fprintf(stderr, "[M::%s] 20\n", __func__); goto redo; } } @@ -17624,12 +17651,13 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) __func__, asm_opt.recover_atg_cov_min); } - + // fprintf(stderr, "[M::%s] 0\n", __func__); adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t, cov->is_r_het); + // fprintf(stderr, "[M::%s] 1\n", __func__); ///primary_flag = get_utg_attributes(*ug, read_g, coverage_cut, sources, ruIndex); update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, cov->is_r_het, 0, DOUBLE_CHECK_THRES, flag, drop_rate); - + // fprintf(stderr, "[M::%s] 2\n", __func__); nsg = (*ug)->g; n_vtx = nsg->n_seq; for (v = 0; v < n_vtx; ++v) @@ -17638,36 +17666,44 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) nsg->seq[v].c = PRIMARY_LABLE; EvaluateLen((*ug)->u, v) = (*ug)->u.a[v].n; } + // fprintf(stderr, "[M::%s] 3\n", __func__); clean_trio_untig_graph(*ug, read_g, coverage_cut, sources, reverse_sources, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, chimeric_rate, 0, 0, drop_ratio, flag, drop_rate, cov); - + // fprintf(stderr, "[M::%s] 4\n", __func__); ///delete_useless_nodes(ug); delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex); + // fprintf(stderr, "[M::%s] 5\n", __func__); update_hap_label(*ug, read_g); + // fprintf(stderr, "[M::%s] 6\n", __func__); update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, cov->is_r_het, 0, DOUBLE_CHECK_THRES, flag, drop_rate); + // fprintf(stderr, "[M::%s] 7\n", __func__); force_trio_clean((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, flag, 0.55, 0.01, 5); - + // fprintf(stderr, "[M::%s] 8\n", __func__); ///if(flag == MOTHER) print_debug_gfa(read_g, *ug, coverage_cut, "debug_trio_1", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); renew_utg(ug, read_g, new_rtg_edges); - + // fprintf(stderr, "[M::%s] 9\n", __func__); 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, 1, NULL, b_mask_t); + // fprintf(stderr, "[M::%s] 10\n", __func__); renew_utg(ug, read_g, new_rtg_edges); + // fprintf(stderr, "[M::%s] 11\n", __func__); rescue_contained_reads_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang, min_ovlp, 10, 0, 1, NULL, NULL, b_mask_t); + // fprintf(stderr, "[M::%s] 12\n", __func__); renew_utg(ug, read_g, new_rtg_edges); + // fprintf(stderr, "[M::%s] 13\n", __func__); } ///if(flag == MOTHER) print_untig_by_read(*ug, "m64043_200627_000137/124716590/ccs", 2789716, NULL, NULL, "beg"); @@ -17675,13 +17711,17 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, cov->is_r_het, 1, FINAL_DOUBLE_CHECK_THRES, flag, drop_rate); + // fprintf(stderr, "[M::%s] 14\n", __func__); update_hap_label(NULL, read_g); + // fprintf(stderr, "[M::%s] 15\n", __func__); renew_utg(ug, read_g, new_rtg_edges); + // fprintf(stderr, "[M::%s] 16\n", __func__); ///delete_useless_nodes(ug); delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex); + // fprintf(stderr, "[M::%s] 17\n", __func__); if(asm_opt.purge_level_trio == 1) @@ -17692,13 +17732,15 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) ///delete_useless_nodes(ug); delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex); } - + // fprintf(stderr, "[M::%s] 18\n", __func__); set_drop_trio_flag(*ug); + // fprintf(stderr, "[M::%s] 19\n", __func__); destory_hap_cov_t(&cov); - + // fprintf(stderr, "[M::%s] 20\n", __func__); // purge_dump(*ug); renew_utg(ug, read_g, new_rtg_edges); + // fprintf(stderr, "[M::%s] 21\n", __func__); } @@ -17736,28 +17778,23 @@ char *f_prefix, uint8_t *kpt_buf, kvec_asg_arc_t_warp *r_edges) kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); - adjust_utg_by_trio(&ug, sg, flag, TRIO_THRES, sources, reverse_sources, coverage_cut, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, max_hang, min_ovlp, &new_rtg_edges, b_mask_t); - if(asm_opt.b_low_cov > 0) { break_ug_contig(&ug, sg, &R_INF, coverage_cut, sources, ruIndex, &new_rtg_edges, max_hang, min_ovlp, &asm_opt.b_low_cov, NULL, asm_opt.m_rate); } - if(asm_opt.b_high_cov > 0) { break_ug_contig(&ug, sg, &R_INF, coverage_cut, sources, ruIndex, &new_rtg_edges, max_hang, min_ovlp, NULL, &asm_opt.b_high_cov, asm_opt.m_rate); } - if(kpt_buf) { update_dump_trio(R_INF.trio_flag, sg->n_seq, kpt_buf, ug); } - if(is_bench) { free(gfa_name); @@ -17775,7 +17812,6 @@ char *f_prefix, uint8_t *kpt_buf, kvec_asg_arc_t_warp *r_edges) ///debug_untig_length(ug, tipsLen, gfa_name); ///print_untig_by_read(ug, "m64011_190901_095311/125831121/ccs", 2310925, "end"); ma_ug_seq(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 1); - ma_ug_print(ug, sg, coverage_cut, sources, ruIndex, (flag==FATHER?"h1tg":"h2tg"), output_file); fclose(output_file); @@ -21075,7 +21111,7 @@ uint32_t type) kvec_t(uint32_t) u_vecs; uint32_t i = 0, maxEvaluateLen = 0, maxBaseLen = 0, v, totalEvaluateLen = 0; u_vecs.a = a, u_vecs.n = n; - if(u_vecs.n > 0) + if(u_vecs.n > 0) ///u_vecs does not contain startID && endId { v_x = &(ug->u.a[u_vecs.a[0]>>1]); asg_seq_del(nsg, u_vecs.a[0]>>1); @@ -23476,7 +23512,6 @@ buf_t* bb, uint8_t* visit, uint32_t trio_flag, long long totalNodeLen, long long } n_reduce = drop_useless_edges(ug, bb, visit, beg, end, nodes, nodes_n, totalNodeLen); } - // if(pop_bubble_at_tangle(ug, 10000000, nodes, nodes_n, beg, end, trio_flag, DROP)!=0) // { // fprintf(stderr, "false bubble popping: beg: %u, end: %u\n", beg>>1, end>>1); @@ -23497,19 +23532,18 @@ buf_t* bb, uint8_t* visit, uint32_t trio_flag, long long totalNodeLen, long long v = w; if(v == beg) break; } - return is_found; } uint32_t cut_edges_progressive(ma_ug_t *ug, kvec_t_u64_warp* edges, float drop_ratio, uint64_t* nodes, uint64_t nodes_n, uint32_t beg, uint32_t end, buf_t* bb, uint8_t* visit, -uint32_t trio_flag, long long totalNodeLen, long long totalBaseLen) +uint32_t trio_flag, long long totalNodeLen, long long totalBaseLen, uint32_t max_arc_n) { uint32_t k, v, i, kv, w, nv, ov_max, ban, is_found = 0; asg_arc_t *a = NULL, *av = NULL; asg_t* nsg = ug->g; radix_sort_arch64(edges->a.a, edges->a.a + edges->a.n); - for (k = 0; k < edges->a.n; k++) + for (k = 0; k < edges->a.n && k < max_arc_n; k++) { a = &nsg->arc[(uint32_t)edges->a.a[k]]; if(a->del) continue; @@ -23586,20 +23620,17 @@ float drop_ratio) // if(debug_is_circle==0) fprintf(stderr, "Not circle: beg: %u, end: %u\n", beg>>1, end>>1); // if(debug_is_circle==1) fprintf(stderr, "Circle: beg: %u, end: %u\n", beg>>1, end>>1); /**************************debug**************************/ - is_found = process_tangles(ug, nodes, nodes_n, beg, end, bb, visit, trio_flag, totalNodeLen, totalBaseLen); - ///if(is_found == 1) fprintf(stderr, "***Found: beg>>1: %u, end>>1: %u\n", beg>>1, end>>1); if(is_found == 0) { is_found = cut_edges_progressive(ug, edges, drop_ratio, nodes, nodes_n, - beg, end, bb, visit, trio_flag, totalNodeLen, totalBaseLen); + beg, end, bb, visit, trio_flag, totalNodeLen, totalBaseLen, 48); ///if(is_found == 1) fprintf(stderr, "***Cutting Found: beg>>1: %u, end>>1: %u\n", beg>>1, end>>1); } - if(is_found == 1) { @@ -23631,10 +23662,8 @@ float drop_ratio) is_found = 0; } } - recover_edges(nsg, edges, nodes, nodes_n, beg, end, beg_c, end_c, 1-is_found); - } void resolve_tangles(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig, @@ -23674,9 +23703,8 @@ uint32_t trio_flag, float drop_ratio) uint8_t* visit = NULL; visit = (uint8_t*)malloc(sizeof(uint8_t) * nsg->n_seq); - uint32_t n_reduce, flag; + uint32_t n_reduce, flag, dbg_round = 0; n_vtx = nsg->n_seq * 2; - while (1) { n_reduce = 0; @@ -23699,8 +23727,8 @@ uint32_t trio_flag, float drop_ratio) } } if(n_reduce == 0) break; + dbg_round++; } - asg_cleanup(nsg); asg_symm(nsg); @@ -23720,7 +23748,6 @@ uint32_t trio_flag, float drop_ratio) } } - n_vtx = nsg->n_seq; for (v = 0; v < n_vtx; ++v) { @@ -23740,7 +23767,6 @@ uint32_t trio_flag, float drop_ratio) } nsu->n = m; } - n_vtx = nsg->n_seq; for (i = 0; i < n_vtx; ++i) { @@ -23764,7 +23790,6 @@ uint32_t trio_flag, float drop_ratio) if(get_real_length(nsg, w^1, NULL) != 1) continue; if(IsMerge(ug->u, w>>1) != 0) continue; end = w; - ///we have three types of merged nodes ///1) CONVEX_M: one direction has two out-nodes, another direction has one out-node ///2) UNROLL_E: one direction has one out-node, another direction doesn't has out-node @@ -23774,13 +23799,11 @@ uint32_t trio_flag, float drop_ratio) unroll_tangle(src, ug, nsu->a, nsu->n, beg, end, &e_vecs, trio_flag, &b_0, visit, drop_ratio); } - for (v = 0; v < ruIndex->len; v++) { get_R_to_U(ruIndex, v, &uId, &is_Unitig); if(is_Unitig == 1) ruIndex->index[v] = (uint32_t)-1; } - kv_destroy(u_vecs.a); kv_destroy(e_vecs.a); @@ -23799,7 +23822,6 @@ uint32_t trio_flag, float drop_ratio) EvaluateLen(src->u, v) = src->u.a[v].n; } ///print_untig_by_read(src, "m64076_200203_181219/82511682/ccs", 2429597, NULL, NULL, "end-1"); - } @@ -24842,18 +24864,20 @@ int load_ruIndex(R_to_U* ruIndex, char* read_file_name) int f_flag = 0; f_flag += fread(&(ruIndex)->len, sizeof((ruIndex)->len), 1, fp); (ruIndex)->index = (uint32_t*)malloc(sizeof(uint32_t)*(ruIndex)->len); - f_flag += fread((ruIndex)->index, sizeof((ruIndex)->index[0]), (ruIndex)->len, fp); + f_flag += fread((ruIndex)->index, sizeof((*((ruIndex)->index))), (ruIndex)->len, fp); if(!(asm_opt.ar)) { 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); + f_flag += fread(R_INF.trio_flag, sizeof((*(R_INF.trio_flag))), (ruIndex)->len, fp); + } else { + fseek(fp, sizeof((*(R_INF.trio_flag)))*(ruIndex)->len, SEEK_CUR); } // CALLOC(ruIndex->is_het, ruIndex->len); // f_flag += fread(ruIndex->is_het, 1, ruIndex->len, fp); - // f_flag += fread(&(asm_opt.hom_global_coverage_set), sizeof(asm_opt.hom_global_coverage_set), 1, fp); - // f_flag += fread(&(asm_opt.hom_global_coverage), sizeof(asm_opt.hom_global_coverage), 1, fp); + f_flag += fread(&(asm_opt.hom_global_coverage_set), sizeof(asm_opt.hom_global_coverage_set), 1, fp); + f_flag += fread(&(asm_opt.hom_global_coverage), sizeof(asm_opt.hom_global_coverage), 1, fp); free(index_name); fflush(fp); @@ -29291,9 +29315,8 @@ int max_hang, int min_ovlp, bubble_type* bub, long long gap_fuzz) void rescue_bubble_by_chain(asg_t *sg, ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, -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, uint32_t chainLenThres, long long gap_fuzz, -bub_label_t* b_mask_t) +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, uint32_t chainLenThres, long long gap_fuzz, bub_label_t* b_mask_t, long long no_trio_recover) { kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); @@ -29327,7 +29350,7 @@ bub_label_t* b_mask_t) beg_idx = bub.f_bub; occ = bub.b_bub + bub.b_end_bub + bub.tangle_bub; rescue_bubbles_by_missing_ovlp_backward(ug, sg, sources, coverage_cut, ruIndex, max_hang, min_ovlp, chainLenThres, beg_idx, occ, &bub, b_mask_t); - if(ha_opt_triobin(&asm_opt)) + if((!no_trio_recover) && (ha_opt_triobin(&asm_opt))) { ma_ug_destroy(ug); ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); reset_bub(&bub, ug, cov->t_ch, &new_rtg_edges); @@ -31675,6 +31698,25 @@ ma_hit_t_alloc* src, uint64_t* readLen, R_to_U* ruIndex, bub_label_t *b_mask_t, return sg; } +void renew_g(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, long long *n_read, +uint64_t **readLen, ma_sub_t **coverage_cut, R_to_U *ruIndex, asg_t **sg, +int64_t mini_overlap_length, int64_t max_hang_length, +ug_opt_t *uopt, int64_t clean_round, double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, +int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file) +{ + ma_ug_t *iug = ul_realignment_gfa(uopt, *sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, + asm_opt.max_short_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file); + asg_t *ng = gen_ng(iug, *sg, uopt, coverage_cut, ruIndex, 256); + ma_ug_destroy(iug); asg_destroy(*sg); + (*sources) = R_INF.paf; + (*reverse_sources) = R_INF.reverse_paf; + (*n_read) = R_INF.total_reads; + (*readLen) = R_INF.read_length; + (*sg) = ng; + ma_hit_contained_advance(*sources, *n_read, *coverage_cut, ruIndex, max_hang_length, mini_overlap_length); + post_rescue(uopt, *sg, (*sources), (*reverse_sources), ruIndex, b_mask_t, 0); +} + void clean_graph( int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long n_read, uint64_t* readLen, long long mini_overlap_length, @@ -31770,11 +31812,17 @@ ma_sub_t **coverage_cut_ptr, int debug_g) debug_gfa:; gen_ug_opt_t(&uopt, sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, &b_mask_t); - set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length); + // set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length); } if(asm_opt.ar) { - ul_realignment_gfa(&uopt, sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, - asm_opt.max_short_tip, &b_mask_t, ha_opt_triobin(&asm_opt), o_file); + renew_g(&sources, &reverse_sources, &n_read, &readLen, &coverage_cut, ruIndex, &sg, mini_overlap_length, max_hang_length, + &uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, asm_opt.max_short_tip, &b_mask_t, + ha_opt_triobin(&asm_opt), o_file); + + // ma_ug_t *iug = ul_realignment_gfa(&uopt, sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, + // asm_opt.max_short_tip, &b_mask_t, ha_opt_triobin(&asm_opt), o_file); + // gen_ng(iug, sg, &uopt, &coverage_cut, ruIndex, 100); + // exit(1); } // print_debug_gfa(sg, NULL, coverage_cut, "UL.debug", sources, ruIndex, max_hang_length, mini_overlap_length, 0, 0, 0); /** diff --git a/Overlaps.h b/Overlaps.h index 3a05afa..2cc2188 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -32,7 +32,7 @@ // #define PRIMARY_LABLE 1 // #define ALTER_LABLE 2 // #define HAP_LABLE 4 - +#define HA_RE_UL_ID "re" #define Get_qn(RECORD) ((uint32_t)((RECORD).qns>>32)) #define Get_qs(RECORD) ((uint32_t)((RECORD).qns)) @@ -880,7 +880,7 @@ ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, int m void rescue_bubble_by_chain(asg_t *sg, ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 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, uint32_t chainLenThres, long long gap_fuzz, -bub_label_t* b_mask_t); +bub_label_t* b_mask_t, long long no_trio_recover); typedef struct{ double weight; @@ -1140,7 +1140,8 @@ void hic_clean(asg_t* read_g); int64_t count_edges_v_w(asg_t *g, uint32_t v, uint32_t w); void renew_utg(ma_ug_t **ug, asg_t* read_g, kvec_asg_arc_t_warp* edge); void merge_unitig_content(ma_utg_t* collection, ma_ug_t* ug, asg_t* read_g, kvec_asg_arc_t_warp* edge); - +void reset_bub_label_t(bub_label_t* x, asg_t *g, uint64_t bub_dist, uint32_t check_cross); +void set_reverse_overlap(ma_hit_t* dest, ma_hit_t* source); // void break_ug_contig(ma_ug_t **ug, asg_t *read_g, All_reads *RNF, ma_sub_t *coverage_cut, // ma_hit_t_alloc* sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, // int* b_low_cov, int* b_high_cov, double m_rate); diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index 0ef3be0..42c91af 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -5623,7 +5623,7 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle if(asm_opt.polyploidy <= 2) { - mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL, 1); + mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL, 1, 0); } if(collect_p_trans && collect_p_trans_f == 1) diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 6b6b4fe..d13a2ef 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -28,6 +28,14 @@ KDQ_INIT(uint64_t) typedef struct { size_t n, m; char *a; } asgc8_v; +typedef struct { + asg64_v cnt; + asg64_v idx_a; + uint64_t idx_n; + asg_t *ext; + uint64_t a_n; +} asg_ext_t; + typedef struct { uint32_t v, uid, off; } usg_arc_mm_t; @@ -2426,6 +2434,20 @@ void print_raw_u2rgfa_seq(all_ul_t *aln, R_to_U* rI, uint32_t is_detail) } +void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover) +{ + rescue_contained_reads_aggressive(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 10, 1, 0, NULL, NULL, b_mask_t); + rescue_missing_overlaps_aggressive(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 1, 0, NULL, b_mask_t); + rescue_missing_overlaps_backward(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 10, 1, 0, b_mask_t); + // rescue_wrong_overlaps_to_unitigs(NULL, sg, sources, reverse_sources, coverage_cut, ruIndex, + // max_hang_length, mini_overlap_length, bubble_dist, NULL); + // rescue_no_coverage_aggressive(sg, sources, reverse_sources, &coverage_cut, ruIndex, max_hang_length, + // mini_overlap_length, bubble_dist, 10); + set_hom_global_coverage(&asm_opt, sg, uopt->coverage_cut, src, rev, rI, uopt->max_hang, uopt->min_ovlp); + rescue_bubble_by_chain(sg, uopt->coverage_cut, src, rev, (asm_opt.max_short_tip*2), 0.15, 3, rI, 0.05, 0.9, uopt->max_hang, + uopt->min_ovlp, 10, uopt->gap_fuzz, b_mask_t, no_trio_recover); +} + void ul_clean_gfa(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, int64_t clean_round, double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, int32_t is_ou, int32_t is_trio, uint32_t ou_thres, char *o_file) { @@ -2522,7 +2544,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); if(!is_ou) asg_cut_semi_circ(sg, LIM_LEN, 1); - + /** rescue_contained_reads_aggressive(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 10, 1, 0, NULL, NULL, b_mask_t); rescue_missing_overlaps_aggressive(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 1, 0, NULL, b_mask_t); rescue_missing_overlaps_backward(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 10, 1, 0, b_mask_t); @@ -2532,6 +2554,9 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i // mini_overlap_length, bubble_dist, 10); set_hom_global_coverage(&asm_opt, sg, uopt->coverage_cut, src, rev, rI, uopt->max_hang, uopt->min_ovlp); rescue_bubble_by_chain(sg, uopt->coverage_cut, src, rev, (asm_opt.max_short_tip*2), 0.15, 3, rI, 0.05, 0.9, uopt->max_hang, uopt->min_ovlp, 10, uopt->gap_fuzz, b_mask_t); + **/ + post_rescue(uopt, sg, src, rev, rI, b_mask_t, is_ou); + output_unitig_graph(sg, uopt->coverage_cut, o_file, src, rI, uopt->max_hang, uopt->min_ovlp); // flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL; flat_soma_v(sg, src, rI); @@ -3381,7 +3406,7 @@ void init_ul_str_idx_t(ul_resolve_t *p) for (k = 0; k < x_n; k++) { if(x_a[k].base || (!x_a[k].el) || (!x_a[k].pchain)) continue; l = k; l <<= 32; l += (x_a[k].hid<<1); l += x_a[k].rev; str->idx.a[x_a[k].hid]++; - kv_push(uint64_t, str->str.a[i], l); + kv_push(uint64_t, str->str.a[i], l);///idx|hid|rev } str->str.a[i].cn = str->str.a[i].n; str->str.a[i].is_cir = 0; } @@ -3757,7 +3782,7 @@ ma_ug_t *ug, ul_str_idx_t *str_idx, all_ul_t *ul_idx, ul_chain_t *res) } } // if(tid == 269 || tid == 276 || tid == 277 || tid == 278) fprintf(stderr, "tid->%u, i->%ld, an->%ld\n", tid, i, a_n); - if(i >= a_n) { + if(i >= a_n) {//the whole chain is co-inear sc += a[0].sc; res->v = a[0].tn_rev_qk>>32; res->s = offset; res->e = offset + a_n; res->sc = sc; return 1; @@ -5737,6 +5762,7 @@ void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_ b_n = buf->b.n; buf->sc.n = 0; for (k = 1, z = 0; k <= b_n; k++) { if(k == b_n || (buf->b.a[z].tn_rev_qk>>32) != (buf->b.a[k].tn_rev_qk>>32)) { + ///get the chain for if(integer_chain(qid, buf->b.a + z, k - z, z, buf, ug, str_idx, uidx->idx, &sc) && sc.v != (uint32_t)-1) { if((buf->sc.n > 0) && ((buf->sc.a[buf->sc.n-1].v>>1) == (sc.v>>1))) { if(buf->sc.a[buf->sc.n-1].sc < sc.sc) { @@ -9224,7 +9250,7 @@ asg64_v *b, asg64_v *ub, uint32_t *r_w_c, uint32_t *r_w_m) uint64_t ulg_pop_bubble(ul_resolve_t *uidx, ma_ug_t *ug, uint64_t* i_max_dist, uint32_t max_ext, uint32_t max_ext_hifi, uint32_t skip_hom, asg64_v *in, asg64_v *ib) { - fprintf(stderr, "[M::%s::] Starting...\n", __func__); + // fprintf(stderr, "[M::%s::] Starting...\n", __func__); asg_t *g = ug->g; asg64_v tx = {0,0,0}, tb = {0,0,0}, *ob = NULL, *ub = NULL; uint32_t v, w, n_vtx = g->n_seq<<1, n_arc, nv, i, wc[2], wm[2], mm_c, mm_m, mm_v; @@ -9276,7 +9302,7 @@ uint64_t ulg_pop_bubble(ul_resolve_t *uidx, ma_ug_t *ug, uint64_t* i_max_dist, u free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); if(n_pop) asg_cleanup(g); if(!in) free(tx.a); if(!ib) free(tb.a); - fprintf(stderr, "[M::%s::] Done...\n", __func__); + // fprintf(stderr, "[M::%s::] Done...\n", __func__); return n_pop; } @@ -9640,7 +9666,7 @@ uint32_t skip_hom, uint32_t *max_drop_len, asg64_v *in, asg64_v *ib) if(!in) free(tx.a); if(!ib) free(tb.a); if (cnt > 0) asg_cleanup(g); - fprintf(stderr, "[M::%s::] cnt::%u\n", __func__, cnt); + // fprintf(stderr, "[M::%s::] cnt::%u\n", __func__, cnt); return cnt; } @@ -11326,9 +11352,15 @@ void update_usg_t_threading_0(usg_t *ng, uint64_t *a, uint64_t a_n, uint32_t *oc assert(a_n > 1); assert(occ[a[0]] == 1); kv_pushp(uint64_t, *b, &p); (*p) = a[0]; (*p) <<= 32; (*p) |= a[0]; occ[a[0]]--; + // fprintf(stderr, "[M::%s::] a_n::%lu\n", __func__, a_n); + // for (k = 0; k < a_n; k++) { + // fprintf(stderr, "utg%.6dl(%c)\t", (int32_t)(a[k]>>1)+1, "+-"[a[k]&1]); + // } + // fprintf(stderr, "\n"); + for (k = 1; k + 1 < a_n; k++) { assert(occ[a[k]] == occ[a[k]^1]); assert(occ[a[k]] > 0); - if(occ[a[k]] > 1) { + if(occ[a[k]] > 1) {///copy a node nid = a[k]>>1; nnid = ng->n; kv_push(uint32_t, ng->mp.a[ng->a[nid].mm], nnid); s = push_usg_t_node(ng, nnid); @@ -11352,8 +11384,8 @@ void update_usg_t_threading_0(usg_t *ng, uint64_t *a, uint64_t a_n, uint32_t *oc for (k = 0; k < b->n; k += 2) { v = b->a[k]; w = b->a[k+1]; - if(((v>>32) == ((uint32_t)v)) && ((w>>32) == ((uint32_t)w))) continue; - remap_gen_arcs(ng, v>>32, (w>>32)^1, ((uint32_t)v), ((uint32_t)w)^1); + if(((v>>32) == ((uint32_t)v)) && ((w>>32) == ((uint32_t)w))) continue;//no need arc + remap_gen_arcs(ng, v>>32, (w>>32)^1, ((uint32_t)v), ((uint32_t)w)^1);///v>>32: old id; ((uint32_t)v): new id remap_gen_arcs(ng, w>>32, (v>>32)^1, ((uint32_t)w), ((uint32_t)v)^1); } @@ -11373,6 +11405,18 @@ void update_usg_t_threading_0(usg_t *ng, uint64_t *a, uint64_t a_n, uint32_t *oc } else if(av[i].del == 0) { av[i].del = 1; z = get_usg_arc(ng, av[i].v^1, (av[i].ul>>32)^1); + // if(!(z && (!z->del))) { + // fprintf(stderr, "[M::%s::] z::%u, z->del::%u\n", __func__, z?1:0, z?z->del:1); + // fprintf(stderr, "[M::%s::arc] utg%.6dl(%c)->utg%.6dl(%c)\n", + // __func__, (int32_t)(av[i].ul>>33)+1, "+-"[(av[i].ul>>32)&1], + // (int32_t)(av[i].v>>1)+1, "+-"[(av[i].v)&1]); + // fprintf(stderr, "[M::%s::new] utg%.6dl(%c)->utg%.6dl(%c)\n", + // __func__, (int32_t)(v1>>1)+1, "+-"[v1&1], + // (int32_t)(w1>>1)+1, "+-"[w1&1]); + // fprintf(stderr, "[M::%s::old] utg%.6dl(%c)->utg%.6dl(%c)\n", + // __func__, (int32_t)(v0>>1)+1, "+-"[v0&1], + // (int32_t)(w0>>1)+1, "+-"[w0&1]); + // } assert(z && (!z->del)); z->del = 1; } } @@ -12042,14 +12086,15 @@ static void worker_unique_bridge_check_s(void *data, long i, int tid) // callbac unique_bridge_check_t *uaux = (unique_bridge_check_t *)data; uint64_t *integer_seq = uaux->integer_seq, s, e, v, w, nse, k; uint64_t *gidx = uaux->gidx, *interval = uaux->interval_idx; - bubble_type *bub = uaux->uidx->bub; + bubble_type *bub = uaux->uidx->bub; usg_t *ng = uaux->ng; s = interval[((uint32_t)gidx[i])]>>32; e = ((uint32_t)interval[((uint32_t)gidx[i])]); assert(s < e); for (k = s, v = w = (uint32_t)-1; k <= e; k++) { if(k > s) { - nse = get_arc_support(uaux->uidx, integer_seq[k-1], integer_seq[k]); + nse = get_arc_support(uaux->uidx, (ng->a[integer_seq[k-1]>>1].mm<<1)|(integer_seq[k-1]&1), + (ng->a[integer_seq[k]>>1].mm<<1)|(integer_seq[k]&1)); if(nse < unique_bridge_occ) { // fprintf(stderr, "+utg%.6dl(%c)\tutg%.6dl(%c)\tnse::%lu\n", // ((int32_t)(integer_seq[k-1]>>1))+1, "+-"[integer_seq[k-1]&1], @@ -12058,7 +12103,7 @@ static void worker_unique_bridge_check_s(void *data, long i, int tid) // callbac } } if(IF_HOM((integer_seq[k]>>1), *bub)) continue; - w = integer_seq[k]; + w = (ng->a[integer_seq[k]>>1].mm<<1)|(integer_seq[k]&1); if(v != (uint32_t)-1) { nse = get_arc_support(uaux->uidx, v, w); if(nse < unique_bridge_occ) { @@ -12079,7 +12124,7 @@ uint32_t ava_pass_unique_bridge_cov(ul_resolve_t *uidx, usg_t *g, asg64_v *b64, uaux.uidx = uidx; uaux.integer_seq = integer_seq; uaux.gidx = b64->a + g_s; uaux.gidx_n = g_e - g_s; - uaux.interval_idx = b64->a; + uaux.interval_idx = b64->a; uaux.ng = g; kt_for(uidx->str_b.n_thread, worker_unique_bridge_check_s, (&uaux), g_e-g_s);///seq->n > 1 for (i = g_s; i < g_e; i++) {///available intervals within the same cluster @@ -12789,6 +12834,8 @@ void renew_usg_t_bub(ul_resolve_t *uidx, usg_t *ng, uint32_t *id_map, uint8_t *f kv_destroy(e.a); destory_bubbles(uidx->bub); free(uidx->bub); CALLOC(uidx->bub, 1); + // fprintf(stderr, "[M::%s] homozygous read coverage threshold: %d\n", __func__, asm_opt.hom_global_coverage_set? + // asm_opt.hom_global_coverage:(int)(((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE))); // identify_bubbles(ug, uidx->bub, uidx->r_het, NULL); identify_bubbles_recal(uidx->sg, ug, uidx->bub, uidx->r_het, uidx->uopt->sources, uidx->uopt->ruIndex, NULL); // fprintf(stderr, "0[M::%s::] f[51]::%u\n", __func__, ff[51]); @@ -13070,7 +13117,7 @@ void debug_prt_renew_aln(usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint6 uint32_t ug_ext_free(asg64_v *ob, asg64_v *ub, ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint8_t **ff, uint32_t **ng_occ, uint64_t **i_idx, asg64_v *b64, asg64_v *ub64, uint32_t rocc_cut) { - fprintf(stderr, "[M::%s::] rocc_cut::%u\n", __func__, rocc_cut); + // fprintf(stderr, "[M::%s::] rocc_cut::%u\n", __func__, rocc_cut); uint32_t k, n_vtx = ng->n<<1, a_n; unique_bridge_check_t u_aux; REALLOC((*ff), n_vtx); REALLOC((*ng_occ), n_vtx); REALLOC((*i_idx), n_vtx); renew_usg_t_bub(uidx, ng, *ng_occ, *ff, rocc_cut); @@ -13105,7 +13152,7 @@ uint8_t **ff, uint32_t **ng_occ, uint64_t **i_idx, asg64_v *b64, asg64_v *ub64, uint32_t ug_ext_strict(asg64_v *ob, asg64_v *ub, ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint8_t **ff, uint32_t **ng_occ, uint64_t **i_idx, asg64_v *b64, asg64_v *ub64) { - fprintf(stderr, "[M::%s::]\n", __func__); + // fprintf(stderr, "[M::%s::]\n", __func__); uint32_t k, i, n_vtx = ng->n<<1, a_n; REALLOC((*ff), n_vtx); REALLOC((*ng_occ), n_vtx); REALLOC((*i_idx), n_vtx); u2g_hybrid_aln(uidx, ng, ob, ub); b64->n = ub64->n = 0; @@ -13173,6 +13220,405 @@ uint8_t **ff, uint32_t **ng_occ, uint64_t **i_idx, asg64_v *b64, asg64_v *ub64) // if(a_n) usg_cleanup(ng); return a_n; } +/** +typedef struct { + uint32_t nid, ulid; + uint32_t raw_sid, raw_eid; + uint32_t raw_sof, raw_eof; +} usc_t; + +typedef struct { + usc_t *a; + size_t n, m; +} usc_vec_t; + +typedef struct { + // ul_resolve_t *uidx; + usc_t *a; + uint64_t a_n; + ma_ug_t *ug; +} scaf_mul_t; + +int64_t reload_uovl(all_ul_t *x, char* file_name) +{ + char* gfa_name = NULL; MALLOC(gfa_name, strlen(file_name)+100); + sprintf(gfa_name, "%s.%s.ul.ovlp.bin", file_name, HA_RE_UL_ID); + FILE* fp = fopen(gfa_name, "r"); free(gfa_name); + if (!fp) return 0; + uint64_t k, tn; size_t kn; ul_vec_t *p = NULL; + + fseek(fp, sizeof(x->nid.n)*1, SEEK_CUR); + for (k = 0; k < x->nid.n; k++) { + fread(&tn, sizeof(tn), 1, fp); + fseek(fp, sizeof((*(x->nid.a->a)))*tn, SEEK_CUR); + } + + // free(x->ridx.idx.a); x->ridx.idx.a = NULL; + fread(&kn, sizeof(kn), 1, fp); + fseek(fp, sizeof((*(x->ridx.idx.a)))*kn, SEEK_CUR); + + // free(x->ridx.occ.a); x->ridx.occ.a = NULL; + fread(&kn, sizeof(kn), 1, fp); + fseek(fp, sizeof((*(x->ridx.occ.a)))*kn, SEEK_CUR); + + fread(&x->n, sizeof(x->n), 1, fp); + if(x->n > x->m) kv_resize(ul_vec_t, *x, x->n); + for (k = 0; k < x->n; k++) { + p = &(x->a[k]); + fread(&p->dd, sizeof(p->dd), 1, fp); + fread(&p->rlen, sizeof(p->rlen), 1, fp); + + fread(&p->r_base.n, sizeof(p->r_base.n), 1, fp); + if(p->r_base.n > p->r_base.m) kv_resize(uint8_t, p->r_base, p->r_base.n); + fread(p->r_base.a, sizeof((*(p->r_base.a))), p->r_base.n, fp); + + fread(&p->bb.n, sizeof(p->bb.n), 1, fp); + if(p->bb.n > p->bb.m) kv_resize(uc_block_t, p->bb, p->bb.n); + fread(p->bb.a, sizeof((*(p->bb.a))), p->bb.n, fp); + + fread(&p->N_site.n, sizeof(p->N_site.n), 1, fp); + if(p->N_site.n > p->N_site.m) kv_resize(uint32_t, p->N_site, p->N_site.n); + fread(p->N_site.a, sizeof((*(p->N_site.a))), p->N_site.n, fp); + } + fclose(fp); +} + +static void gen_scaffold_id(void *data, long i, int tid) // callback for kt_for() +{ + scaf_mul_t *ss = (scaf_mul_t *)data; + ma_ug_t *ug = ss->ug; + usc_t *z = &(ss->a[i]); + uint32_t v = z->raw_sid, w = z->raw_eid, uv, uw, rev, id; + uint64_t *a, a_n, k; int64_t ql, mmql; + uc_block_t *p, *n, *m0, *m1; + p = n = m0 = m1 = NULL; + a = UL_INF.ridx.occ.a + UL_INF.ridx.idx.a[v>>1]; + a_n = UL_INF.ridx.idx.a[(v>>1)+1] - UL_INF.ridx.idx.a[v>>1]; + for (k = 0, mmql = INT32_MAX; k < a_n; k++) { + p = &(UL_INF.a[a[k]>>32].bb.a[(uint32_t)(a[k])]); assert(p->hid == (v>>1)); + if(UL_INF.a[a[k]>>32].dd != 3) continue;///dd = 3 is saved + if(p->base || (!p->el) || (!p->pchain)) continue; + uv = (((uint32_t)(p->hid))<<1)|((uint32_t)(p->rev)); + + if((uv == v) && (p->aidx != (uint32_t)-1)) { + n = &(UL_INF.a[a[k]>>32].bb.a[p->aidx]); + assert((!n->base)&&(n->el)&&(n->pchain)&&(n->pidx==((uint32_t)(a[k])))); + uw = (((uint32_t)(n->hid))<<1)|((uint32_t)(n->rev)); + if(uw == w) { + ql = ((int64_t)n->qs) - ((int64_t)p->qe); + if(ql < INT32_MAX) { + mmql = ql; rev = 0; m0 = p; m1 = n; id = a[k]>>32; + } + } + } + + if(((uv^1) == v) && (p->pidx != (uint32_t)-1)) { + n = &(UL_INF.a[a[k]>>32].bb.a[p->pidx]); + assert((!n->base)&&(n->el)&&(n->pchain)&&(n->aidx==((uint32_t)(a[k])))); + uw = (((uint32_t)(n->hid))<<1)|((uint32_t)(n->rev)); uw ^= 1; + if(uw == w) { + ql = ((int64_t)p->qs) - ((int64_t)n->qe); + if(ql < INT32_MAX) { + mmql = ql; rev = 1; m0 = n; m1 = p; id = a[k]>>32; + } + } + } + } + + if(m0 && m1) { + z->ulid = (id<<1)|rev; + z->raw_sof = m0->qe; + z->raw_eof = m1->qs; + } +} + +#define B4Lg(x) (((x)>>2)+(((x)&3)?1:0)) +uint32_t load_scaf_base(all_ul_t *x, char* file_name) +{ + char *gfa_name = (char*)malloc(strlen(file_name)+50); + sprintf(gfa_name, "%s.uidx.ucr.bin", file_name); + FILE *fp = fopen(gfa_name, "r"); free(gfa_name); + if (!fp) return 0; + uint64_t rid; uint32_t len; ul_vec_t ss, *z; memset(&ss, 0, sizeof(ss)); + while(1) { + fread(&rid, sizeof(rid), 1, fp); + if(feof(fp)) break; + z = &(x->a[rid]); + fread(&len, sizeof(len), 1, fp); assert(z->rlen == len); + fread(&(ss.N_site.n), sizeof(ss.N_site.n), 1, fp); + kv_resize(uint32_t, ss.N_site, ss.N_site.n); + fread(ss.N_site.a, sizeof((*(ss.N_site.a))), ss.N_site.n, fp); + ss.r_base.n = B4Lg(len); kv_resize(uint8_t, ss.r_base, ss.r_base.n); + fread(ss.r_base.a, sizeof((*(ss.r_base.a))), ss.r_base.n, fp); + if(z->dd != 4) continue; + + + kv_resize(uint32_t, z->N_site, ss.N_site.n); z->N_site.n = ss.N_site.n; + memcpy(z->N_site.a, ss.N_site.a, sizeof((*(z->N_site.a)))*z->N_site.n); + + kv_resize(uint8_t, z->r_base, ss.r_base.n); z->r_base.n = ss.r_base.n; + memcpy(z->r_base.a, ss.r_base.a, sizeof((*(z->r_base.a)))*z->r_base.n); + } + // load_compress_base_disk(fp, &rid, des.a, &ulen, &(sl->ucr_s->u)); + + fclose(fp); free(ss.N_site.a); free(ss.r_base.a); +} + +int64_t fill_scaffolds(usc_t *a, uint64_t a_n, ma_ug_t *raw_g) +{ + scaf_mul_t ss; uint32_t k, occ = 0; + reload_uovl(&UL_INF, asm_opt.output_file_name); + filter_ul_ug(raw_g); + free(UL_INF.ridx.idx.a); UL_INF.ridx.idx.n = UL_INF.ridx.idx.m = 0; + free(UL_INF.ridx.occ.a); UL_INF.ridx.occ.n = UL_INF.ridx.occ.m = 0; + gen_ul_vec_rid_t(&UL_INF, NULL, raw_g); + ss.a = a; ss.a_n = a_n; ss.ug = raw_g; + kt_for(asm_opt.thread_num, gen_scaffold_id, &ss, a_n); + for (k = 0; k < a_n; k++) { + if(a[k].ulid == ((uint32_t)-1)) continue; + if(a[k].raw_eof <= a[k].raw_sof) continue; + UL_INF.a[a[k].ulid>>1].dd = 4;///load + occ++; + } + + if(occ) load_scaf_base(&UL_INF, asm_opt.output_file_name); +} + +uint64_t reset_scaf_node_uinfo_srt_t(usc_t *z, int64_t min_arc_len, int64_t scaf_len, int64_t *nlen) +{ + (*nlen) = scaf_len + (min_arc_len<<1); + if(z->ulid == ((uint32_t)-1)) return 0;///a scaffold node + if(z->raw_eof + min_arc_len <= z->raw_sof) {///has an overlap longer than min_arc_len + (*nlen) = scaf_len; + return 1;///a scaffold node; no need this node, could directly use existing nodes + } + int64_t s, e; + s = z->raw_sof; e = z->raw_eof; + (*nlen) = (min_arc_len<<1) + (e - s); + return 2; +} + + +int64_t push_scaf_bases(uint8_t *des, All_reads* rdb, usc_t *z, int64_t min_arc_len, int64_t scaf_len, int64_t nlen, UC_Read *g_read) +{ + int64_t nlen0, ff; + ff = reset_scaf_node_uinfo_srt_t(z, min_arc_len, scaf_len, &nlen0); assert(nlen0 == nlen); + if(ff == 1) {///has an overlap longer than min_arc_len between z->raw_sid and z->raw_eid + memset(des, 0, sizeof((*(des))*(nlen/4+1))); + return ff; + } + if(z->raw_eof >= z->raw_sof) { + + } + return ff; +} + +void realloc_rdb_adv(All_reads* rdb, ma_sub_t **cov, R_to_U *ruI, uint64_t *rmap, uint64_t rid_n, +uint64_t scaf_len, char *scaf_id, asg_t *ng, ug_opt_t *uopt, usc_t *a) +{ + uint64_t i, rid_n0 = rdb->total_reads, tname, cname; usc_t *z; + uint64_t scaf_id_len = strlen(scaf_id); char *des, *src; int64_t nlen, ff; + UC_Read g_read; init_UC_Read(&g_read); + rdb->total_reads = rid_n; + tname = rdb->name_index[rid_n0]; + fprintf(stderr, "+[M::%s] rid_n0::%lu, rid_n::%lu\n", __func__, rid_n0, rid_n); + ///for read bases + REALLOC(rdb->N_site, rdb->total_reads); + REALLOC(rdb->read_length, rdb->total_reads); + REALLOC(rdb->read_size, rdb->total_reads); + REALLOC(rdb->read_sperate, rdb->total_reads); + REALLOC(rdb->trio_flag, rdb->total_reads); + REALLOC(rdb->name_index, rdb->total_reads+1);///total_reads+1 + REALLOC((*cov), rdb->total_reads); + for (i = rid_n0; i < rdb->total_reads; i++) { + // fprintf(stderr, "[M::%s] i::%lu\n", __func__, i); + rdb->N_site[i] = NULL; + // rdb->read_length[i] = scaf_len; + // rdb->read_size[i] = scaf_len; + rdb->trio_flag[i] = AMBIGU; + (*cov)[i].c = (*cov)[i].del = 0; + // (*cov)[i].s = 0; (*cov)[i].e = scaf_len; + cname = scaf_id_len; + if(rmap[i] != ((uint64_t)-1)) {///not a scaffold node + if(rdb->N_site[rmap[i]] != NULL) { + MALLOC(rdb->N_site[i], rdb->N_site[rmap[i]][0]+1); + memcpy(rdb->N_site[i], rdb->N_site[rmap[i]], + sizeof((*(rdb->N_site[i])))*(rdb->N_site[rmap[i]][0]+1)); + } + rdb->read_length[i] = rdb->read_length[rmap[i]]; + rdb->read_size[i] = rdb->read_length[rmap[i]]; + (*cov)[i].s = (*cov)[rmap[i]].s; + (*cov)[i].e = (*cov)[rmap[i]].e; + rdb->trio_flag[i] = rdb->trio_flag[rmap[i]]; + (*cov)[i] = (*cov)[rmap[i]]; + cname = Get_NAME_LENGTH((*rdb), (rmap[i])); + } else { + z = &(a[(i-ng->r_seq)]); assert(z->nid == i); + reset_scaf_node_uinfo_srt_t(z, uopt->min_ovlp, scaf_len, &nlen); + rdb->read_length[i] = nlen; + rdb->read_size[i] = nlen; + (*cov)[i].s = 0; (*cov)[i].e = nlen; + ng->seq[i].len = nlen;///update length + } + rdb->name_index[i] = tname; tname += cname; + rdb->total_reads_bases += rdb->read_length[i]; + + MALLOC(rdb->read_sperate[i], (rdb->read_length[i]/4+1)); + if(rmap[i] != ((uint64_t)-1)) {///not a scaffold node + memcpy(rdb->read_sperate[i], rdb->read_sperate[rmap[i]], + sizeof((*(rdb->read_sperate[i])))*(rdb->read_length[i]/4+1)); + } else { + ///set to A + z = &(a[(i-ng->r_seq)]); assert(z->nid == i); + ff = push_scaf_bases(rdb->read_sperate[i], rdb, z, uopt->min_ovlp, scaf_len, rdb->read_length[i], &g_read); + if(ff == 2) ng->seq[i].del = 1; + // memset(rdb->read_sperate[i], 0, sizeof((*(rdb->read_sperate[i])))*(rdb->read_length[i]/4+1)); + } + } + + rdb->index_size = rdb->total_reads; + rdb->name_index[i] = tname; + rdb->total_name_length = tname; + rdb->name_index_size = rdb->total_reads+1; + REALLOC(rdb->name, tname); + for (i = rid_n0; i < rdb->total_reads; i++) { + des = Get_NAME((*rdb), i); src = scaf_id; cname = scaf_id_len; + if(rmap[i] != ((uint64_t)-1)) { + src = Get_NAME((*rdb), rmap[i]); cname = Get_NAME_LENGTH((*rdb), (rmap[i])); + } + memcpy(des, src, sizeof((*(des)))*cname); + } + + REALLOC(rdb->paf, rdb->total_reads); + memset(rdb->paf+rid_n0, 0, (rdb->total_reads-rid_n0)*sizeof((*rdb->paf))); + REALLOC(rdb->reverse_paf, rdb->total_reads); + memset(rdb->reverse_paf+rid_n0, 0, (rdb->total_reads-rid_n0)*sizeof((*rdb->reverse_paf))); + + ruI->len = rdb->total_reads; + REALLOC(ruI->index, ruI->len); + memset(ruI->index, -1, sizeof((*(ruI->index)))*(ruI->len)); + + reset_bub_label_t(uopt->b_mask_t, ng, 0, 0); + uopt->coverage_cut = (*cov); + uopt->reverse_sources = rdb->reverse_paf; + uopt->sources = rdb->paf; + destory_UC_Read(&g_read); + fprintf(stderr, "-[M::%s] rid_n0::%lu, rid_n::%lu\n", __func__, rid_n0, rid_n); +} + +asg_t *renew_ng(usg_t *tg, ma_ug_t *raw_g, asg_t *sg, ug_opt_t *uopt, ma_sub_t **cov, R_to_U *ruI, uint64_t scaffold_len) +{ + ma_utg_t *u; uint64_t i, v, w, m, h, z, raw_v, raw_w, nocc, nv, vx, wx; int32_t r; + asg_arc_t t, *p; usg_arc_t *av = NULL; uint64_t slen = scaffold_len + (uopt->min_ovlp*2); + asg_ext_t ext; memset(&ext, 0, sizeof(ext)); ext.ext = asg_init(); asg_t *ng = ext.ext; + usc_vec_t sa; memset(&sa, 0, sizeof(sa)); usc_t *psa; + ext.a_n = sg->n_seq; + ext.cnt.n = ext.cnt.m = ext.a_n; + CALLOC(ext.cnt.a, ext.cnt.n);///count + + ext.idx_a.n = ext.idx_a.m = sg->n_seq; + MALLOC(ext.idx_a.a, ext.idx_a.n); + memset(ext.idx_a.a, -1, sizeof(*(ext.idx_a.a))*ext.idx_a.n);//map + + for (i = 0; i < tg->n; ++i) { + u = &(raw_g->u.a[tg->a[i].mm]); + for (m = 0; m < u->n; m++) { + v = u->a[m]>>32; + ext.cnt.a[v>>1]++; + if(ext.cnt.a[v>>1] == 1) { + h = v>>1; ext.idx_a.a[h] = v>>1; + } else { + h = ext.idx_a.n; + kv_push(uint64_t, ext.idx_a, (v>>1)); + } + + z = (h<<1)|(v&1); z <<= 32; z |= ((uint32_t)u->a[m]); + u->a[m] = z; + } + } + + for (i = 0; i < sg->n_seq; ++i) { + asg_seq_set(ng, i, sg->seq[i].len, ext.idx_a.a[i]==((uint64_t)-1)?1:0); + ng->seq[i].c = 0; + if(ext.idx_a.a[i]!=((uint64_t)-1)) assert(ext.idx_a.a[i] == i); + ext.idx_a.a[i] = i;///set for delted read + } + for (; i < ext.idx_a.n; i++) { + asg_seq_set(ng, i, sg->seq[ext.idx_a.a[i]].len, 0); + ng->seq[i].c = 0; assert(!(ng->seq[ext.idx_a.a[i]].del)); + } + ng->r_seq = ng->n_seq; + assert(ng->n_seq == ext.idx_a.n); + + + for (i = 0, nocc = ng->r_seq, sa.n = 0; i < tg->n; ++i) { + ///there shouldn't any scaffolding within the nodes + v = i<<1; nv = usg_arc_n(tg, v); av = usg_arc_a(tg, v); + for (m = 0; m < nv; m++) { + if(av[m].del) continue; + w = av[m].v; + vx = (v&1?((raw_g->u.a[tg->a[v>>1].mm].a[0]>>32)^1): + (raw_g->u.a[tg->a[v>>1].mm].a[raw_g->u.a[tg->a[v>>1].mm].n-1]>>32)); + wx = (w&1?((raw_g->u.a[tg->a[w>>1].mm].a[raw_g->u.a[tg->a[w>>1].mm].n-1]>>32)^1): + (raw_g->u.a[tg->a[w>>1].mm].a[0]>>32)); + raw_v = (ext.idx_a.a[vx>>1]<<1)|(vx&1); + raw_w = (ext.idx_a.a[wx>>1]<<1)|(wx&1); + assert((raw_v>>1) < sg->n_seq); + assert((raw_w>>1) < sg->n_seq); + if(vx > wx) continue; + if((vx == wx) && (vx&1)) continue; + if(gen_spec_edge(sg, uopt, raw_v, raw_w, &t) < 0) { + asg_seq_set(ng, nocc, slen, 0); + ng->seq[nocc].c = 0; + kv_push(uint64_t, ext.idx_a, ((uint64_t)-1)); + + kv_pushp(usc_t, sa, &psa); + psa->ulid = (uint32_t)-1; psa->nid = nocc; + psa->raw_sid = (tg->a[v>>1].mm<<1)|(v&1); + psa->raw_eid = (tg->a[w>>1].mm<<1)|(w&1); + psa->raw_sof = psa->raw_eof = (uint32_t)-1; + nocc++;///a scaffold node + } + } + + v = (i<<1)+1; nv = usg_arc_n(tg, v); av = usg_arc_a(tg, v); + for (m = 0; m < nv; m++) { + if(av[m].del) continue; + w = av[m].v; + vx = (v&1?((raw_g->u.a[tg->a[v>>1].mm].a[0]>>32)^1): + (raw_g->u.a[tg->a[v>>1].mm].a[raw_g->u.a[tg->a[v>>1].mm].n-1]>>32)); + wx = (w&1?((raw_g->u.a[tg->a[w>>1].mm].a[raw_g->u.a[tg->a[w>>1].mm].n-1]>>32)^1): + (raw_g->u.a[tg->a[w>>1].mm].a[0]>>32)); + raw_v = (ext.idx_a.a[vx>>1]<<1)|(vx&1); + raw_w = (ext.idx_a.a[wx>>1]<<1)|(wx&1); + assert((raw_v>>1) < sg->n_seq); + assert((raw_w>>1) < sg->n_seq); + if(vx > wx) continue; + if((vx == wx) && (vx&1)) continue; + if(gen_spec_edge(sg, uopt, raw_v, raw_w, &t) < 0) { + asg_seq_set(ng, nocc, slen, 0); + ng->seq[nocc].c = 0; + kv_push(uint64_t, ext.idx_a, ((uint64_t)-1)); + + kv_pushp(usc_t, sa, &psa); + psa->ulid = (uint32_t)-1; psa->nid = nocc; + psa->raw_sid = (tg->a[v>>1].mm<<1)|(v&1); + psa->raw_eid = (tg->a[w>>1].mm<<1)|(w&1); + psa->raw_sof = psa->raw_eof = (uint32_t)-1; + nocc++;///a scaffold node + } + } + } + + assert(ng->n_seq == ext.idx_a.n); + CALLOC(ng->seq_vis, (ng->n_seq<<1)); + ///# scaffolding nodes + if(sa.n > 0) fill_scaffolds(sa.a, sa.n, raw_g); + realloc_rdb_adv(&(R_INF), cov, ruI, ext.idx_a.a, ext.idx_a.n, slen, (char *)"scaf", ng, uopt, sa.a); + + free(sa.a); +} +**/ void u2g_hybrid_detan_iter(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint32_t clean_round, asg64_v *in, asg64_v *ib) { @@ -13192,7 +13638,7 @@ void u2g_hybrid_detan_iter(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint } // // ug_ext_strict(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64); // ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 50); - // // prt_usg_t(uidx, ng, "ng0"); + // prt_usg_t(uidx, ng, "ng.db"); // ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 10); // ncut += ug_ext_strict(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64); if(ncut) { @@ -13334,7 +13780,7 @@ void u2g_hybrid_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, usg_t *ng, asg64_v * double step = (ulopt->clean_round==1?ulopt->max_ovlp_drop_ratio: ((ulopt->max_ovlp_drop_ratio-ulopt->min_ovlp_drop_ratio)/(ulopt->clean_round-1))); double drop = ulopt->min_ovlp_drop_ratio; ///CALLOC(iug->g->seq_vis, iug->g->n_seq*2); - fprintf(stderr, "\n[M::%s::] Starting hybrid clean, mm_tip::%ld\n", __func__, mm_tip); + // fprintf(stderr, "\n[M::%s::] Starting hybrid clean, mm_tip::%ld\n", __func__, mm_tip); // prt_usg_t(uidx, ng, "ng0"); usg_arc_cut_tips(ng, mm_tip, 0, b); // prt_usg_t(uidx, ng, "ng1"); @@ -13465,9 +13911,9 @@ void merge_hybrid_utg_content(ma_utg_t* cc, ma_ug_t* raw, asg_t* rg, usg_t *ng, ma_ug_t *gen_hybrid_ug(ul_resolve_t *uidx, usg_t *ng) { - fprintf(stderr, "[M::%s::] ng->n::%u\n", __func__, (uint32_t)ng->n); + // fprintf(stderr, "[M::%s::] ng->n::%u\n", __func__, (uint32_t)ng->n); ma_ug_t *ug = ma_ug_hybrid_gen(ng); - fprintf(stderr, "[M::%s::] ug->g->n_seq::%u\n", __func__, (uint32_t)ug->g->n_seq); + // fprintf(stderr, "[M::%s::] ug->g->n_seq::%u\n", __func__, (uint32_t)ug->g->n_seq); uint32_t i; ma_utg_t *u; kvec_asg_arc_t_warp e; kv_init(e.a); e.i = 0; for (i = 0; i < ug->u.n; i++) { ug->g->seq[i].c = PRIMARY_LABLE; @@ -13490,7 +13936,7 @@ void u2g_threading(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint64_t cov_cutoff, as ul2ul_idx_t *idx = &(uidx->uovl); ma_ug_t *iug = idx->i_ug; ma_ug_t *raw = uidx->l1_ug; uint64_t k, z, i, t_s, t_e; uinfo_srt_warp_t *seq; usg_t *ng; CALLOC(ng, 1); usg_seq_t *s; usg_arc_warp *sv; usg_arc_t *p; - asg_arc_t *av; uint32_t nv, v, w; int64_t tt, tl, tm; + asg_arc_t *av; uint32_t nv, v, w, vx, wx; int64_t tt, tl, tm; ng->mp.n = ng->mp.m = raw->g->n_seq; MALLOC(ng->mp.a, ng->mp.n); for (k = 0; k < raw->g->n_seq; k++) { @@ -13560,19 +14006,26 @@ void u2g_threading(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint64_t cov_cutoff, as p = get_usg_arc(ng, seq->a[i+1].v^1, seq->a[i].v^1); if(p->ou < (uint64_t)tm) p->ou = tm; } else { - v = seq->a[i].v; w = seq->a[i+1].v; - kv_pushp(usg_arc_t, (ng->a[v>>1].arc[v&1]), &p); - p->del = 0; p->ou = tm; p->v = w; p->ol = 0; p->idx = 0; - p->ul = (((uint64_t)v)<<32)|(raw->g->seq[v>>1].len); + v = seq->a[i].v; w = seq->a[i+1].v; p = NULL; + vx = (v&1?((raw->u.a[v>>1].a[0]>>32)^1):(raw->u.a[v>>1].a[raw->u.a[v>>1].n-1]>>32)); + wx = (w&1?((raw->u.a[w>>1].a[raw->u.a[w>>1].n-1]>>32)^1):(raw->u.a[w>>1].a[0]>>32)); + ///(v>>1) != (w>>1) or v == w, need to handle (v>>1) == (w>>1) && (v&1) != (w&1) later + if(((v^w) != 1) && (uidx->sg->seq[vx>>1].len > uidx->uopt->min_ovlp) && + (uidx->sg->seq[wx>>1].len > uidx->uopt->min_ovlp)) { + kv_pushp(usg_arc_t, (ng->a[v>>1].arc[v&1]), &p); + p->del = 0; p->ou = tm; p->v = w; p->ol = 0; p->idx = 0; + p->ul = (((uint64_t)v)<<32)|(raw->g->seq[v>>1].len); - v = seq->a[i+1].v^1; w = seq->a[i].v^1; - kv_pushp(usg_arc_t, (ng->a[v>>1].arc[v&1]), &p); - p->del = 0; p->ou = tm; p->v = w; p->ol = 0; p->idx = 0; - p->ul = (((uint64_t)v)<<32)|(raw->g->seq[v>>1].len); + v = seq->a[i+1].v^1; w = seq->a[i].v^1; + kv_pushp(usg_arc_t, (ng->a[v>>1].arc[v&1]), &p); + p->del = 0; p->ou = tm; p->v = w; p->ol = 0; p->idx = 0; + p->ul = (((uint64_t)v)<<32)|(raw->g->seq[v>>1].len); + } + } + if(p) { + pushp_usg_arc_mm(ng, seq->a[i].v, seq->a[i+1].v, k, i); + pushp_usg_arc_mm(ng, seq->a[i+1].v^1, seq->a[i].v^1, k, i); } - - pushp_usg_arc_mm(ng, seq->a[i].v, seq->a[i+1].v, k, i); - pushp_usg_arc_mm(ng, seq->a[i+1].v^1, seq->a[i].v^1, k, i); } } t_s = ub->a[z] + 1; @@ -13643,7 +14096,7 @@ void u2g_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint32_t keep_raw_utg) cnt += ulg_arc_cut_supports(uidx, iug, mm_tip, ulopt->max_tip_hifi, drop, ulopt->is_trio, topo_level, 1, NULL, keep_raw_utg, &bu, &uu); cnt += ulg_arc_cut_tips(uidx, iug, mm_tip, ulopt->max_tip_hifi, 0, &bu, &uu); } - fprintf(stderr, "[M::%s::] Done round-%ld, drop::%f\n", __func__, i, drop); + // fprintf(stderr, "[M::%s::] Done round-%ld, drop::%f\n", __func__, i, drop); } ulg_pop_bubble(uidx, iug, NULL, ((int64_t)0x7fffffff), ulopt->max_tip_hifi, 1, &bu, &uu); @@ -14149,8 +14602,8 @@ ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt, uin kt_for(uidx->str_b.n_thread, worker_integert_debug_sym, uidx, z->tot);///all ul + ug print_integert_ovlp_stat(z); - print_uls_seq(uidx, asm_opt.output_file_name); - print_uls_ovs(uidx, asm_opt.output_file_name); + // print_uls_seq(uidx, asm_opt.output_file_name); + // print_uls_ovs(uidx, asm_opt.output_file_name); z->i_g = integer_sg_gen(uidx, uopt->min_ovlp); asg_arc_del_trans(z->i_g, uopt->gap_fuzz); @@ -14164,7 +14617,7 @@ ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt, uin u2g_clean(uidx, ulopt, keep_raw_utg); // renew_ul2_utg(uidx); - output_integer_graph(uidx, z->i_ug, asm_opt.output_file_name, 0); + // output_integer_graph(uidx, z->i_ug, asm_opt.output_file_name, 0); return z; } @@ -14265,13 +14718,436 @@ ma_ug_t* output_trio_unitig_graph_ul(ug_opt_t *uopt, ul_resolve_t *uidx, char* o return NULL; } + + +void realloc_rdb(All_reads* rdb, ma_sub_t **cov, R_to_U *ruI, uint64_t *rmap, uint64_t rid_n, uint64_t scaf_len, char *scaf_id, +asg_t *ng, ug_opt_t *uopt) +{ + uint64_t i, rid_n0 = rdb->total_reads, tname, cname; + uint64_t scaf_id_len = strlen(scaf_id); char *des, *src; + rdb->total_reads = rid_n; + tname = rdb->name_index[rid_n0]; + fprintf(stderr, "+[M::%s] rid_n0::%lu, rid_n::%lu\n", __func__, rid_n0, rid_n); + ///for read bases + REALLOC(rdb->N_site, rdb->total_reads); + REALLOC(rdb->read_length, rdb->total_reads); + REALLOC(rdb->read_size, rdb->total_reads); + REALLOC(rdb->read_sperate, rdb->total_reads); + REALLOC(rdb->trio_flag, rdb->total_reads); + REALLOC(rdb->name_index, rdb->total_reads+1);///total_reads+1 + REALLOC((*cov), rdb->total_reads); + for (i = rid_n0; i < rdb->total_reads; i++) { + // fprintf(stderr, "[M::%s] i::%lu\n", __func__, i); + rdb->N_site[i] = NULL; + rdb->read_length[i] = scaf_len; + rdb->read_size[i] = scaf_len; + rdb->trio_flag[i] = AMBIGU; + (*cov)[i].c = (*cov)[i].del = 0; + (*cov)[i].s = 0; (*cov)[i].e = scaf_len; + cname = scaf_id_len; + if(rmap[i] != ((uint64_t)-1)) {///not a scaffold node + if(rdb->N_site[rmap[i]] != NULL) { + MALLOC(rdb->N_site[i], rdb->N_site[rmap[i]][0]+1); + memcpy(rdb->N_site[i], rdb->N_site[rmap[i]], + sizeof((*(rdb->N_site[i])))*(rdb->N_site[rmap[i]][0]+1)); + } + rdb->read_length[i] = rdb->read_length[rmap[i]]; + rdb->read_size[i] = rdb->read_length[rmap[i]]; + rdb->trio_flag[i] = rdb->trio_flag[rmap[i]]; + (*cov)[i] = (*cov)[rmap[i]]; + cname = Get_NAME_LENGTH((*rdb), (rmap[i])); + } + rdb->name_index[i] = tname; tname += cname; + rdb->total_reads_bases += rdb->read_length[i]; + + MALLOC(rdb->read_sperate[i], (rdb->read_length[i]/4+1)); + if(rmap[i] != ((uint64_t)-1)) {///not a scaffold node + memcpy(rdb->read_sperate[i], rdb->read_sperate[rmap[i]], + sizeof((*(rdb->read_sperate[i])))*(rdb->read_length[i]/4+1)); + } else { + ///set to A + memset(rdb->read_sperate[i], 0, sizeof((*(rdb->read_sperate[i])))*(rdb->read_length[i]/4+1)); + } + } + + rdb->index_size = rdb->total_reads; + rdb->name_index[i] = tname; + rdb->total_name_length = tname; + rdb->name_index_size = rdb->total_reads+1; + REALLOC(rdb->name, tname); + for (i = rid_n0; i < rdb->total_reads; i++) { + des = Get_NAME((*rdb), i); src = scaf_id; cname = scaf_id_len; + if(rmap[i] != ((uint64_t)-1)) { + src = Get_NAME((*rdb), rmap[i]); cname = Get_NAME_LENGTH((*rdb), (rmap[i])); + } + memcpy(des, src, sizeof((*(des)))*cname); + } + + REALLOC(rdb->paf, rdb->total_reads); + memset(rdb->paf+rid_n0, 0, (rdb->total_reads-rid_n0)*sizeof((*rdb->paf))); + REALLOC(rdb->reverse_paf, rdb->total_reads); + memset(rdb->reverse_paf+rid_n0, 0, (rdb->total_reads-rid_n0)*sizeof((*rdb->reverse_paf))); + + ruI->len = rdb->total_reads; + REALLOC(ruI->index, ruI->len); + memset(ruI->index, -1, sizeof((*(ruI->index)))*(ruI->len)); + + reset_bub_label_t(uopt->b_mask_t, ng, 0, 0); + uopt->coverage_cut = (*cov); + uopt->reverse_sources = rdb->reverse_paf; + uopt->sources = rdb->paf; + fprintf(stderr, "-[M::%s] rid_n0::%lu, rid_n::%lu\n", __func__, rid_n0, rid_n); +} + +inline void update_qtn(ma_hit_t *z, uint64_t qn, uint64_t tn) +{ + z->qns <<= 32; z->qns >>= 32; z->qns |= (qn<<32); z->tn = tn; +} + +void renew_paf0(ma_hit_t_alloc *paf, uint64_t *a0, uint64_t a0n, uint64_t *a1, uint64_t a1n) +{ + if(a0n <= 1 && a1n <= 1) return; + int64_t idx; uint64_t qn = ((uint32_t)a0[0]), tn = ((uint32_t)a1[0]), i, k; + ma_hit_t e01, e10; + idx = get_specific_overlap(&(paf[qn]), qn, tn); + if(idx < 0) return; + e01 = paf[qn].buffer[idx]; + idx = get_specific_overlap(&(paf[tn]), tn, qn); + e10 = paf[tn].buffer[idx]; + + for (i = 0; i < a0n; i++) { + qn = ((uint32_t)a0[i]); + for (k = 0; k < a1n; k++) { + if(i == 0 && k == 0) continue; + tn = ((uint32_t)a1[k]); + update_qtn(&e01, qn, tn); add_ma_hit_t_alloc(&(paf[qn]), &e01); + update_qtn(&e10, tn, qn); add_ma_hit_t_alloc(&(paf[tn]), &e10); + } + } +} + +void renew_paf1(ma_hit_t_alloc *paf, uint64_t *a, uint64_t an, uint64_t len) +{ + uint64_t i, k, qn, tn; ma_hit_t arc; + arc.qns = 0; arc.qe = len; + arc.tn = 0; arc.ts = 0; arc.te = len; + arc.rev = arc.el = arc.ml = arc.no_l_indel = arc.bl = 0; + + for (i = 0; i < an; i++) { + qn = ((uint32_t)a[i]); + for (k = 0; k < an; k++) { + tn = ((uint32_t)a[k]); + update_qtn(&arc, qn, tn); + add_ma_hit_t_alloc(&(paf[qn]), &arc); + } + } +} + +void update_paf(ma_hit_t_alloc *src, ma_hit_t_alloc *r_src, uint64_t *rmap, uint64_t rid_n, uint64_t pre_gn, asg_t *ng) +{ + asg64_v clus; uint64_t k, l, dn = 0, j, qn, tn; + uint64_t *idx; ma_hit_t_alloc *z; + uint64_t *a0, *a1, a0n, a1n; + CALLOC(idx, pre_gn); + kv_init(clus); kv_resize(uint64_t, clus, rid_n); + for (k = 0; k < rid_n; k++) { + if(rmap[k] == ((uint64_t)-1)) continue;///scaffold + clus.a[clus.n++] = (rmap[k]<<32)|k; + } + radix_sort_srt64(clus.a, clus.a+clus.n); + for (k = 1, l = 0; k <= clus.n; k++) { + if(k == clus.n || (clus.a[k]>>32) != (clus.a[l]>>32)) { + idx[(clus.a[l]>>32)] = (l<<32)|k; + l = k; dn++; + } + } + assert(dn == pre_gn); + fprintf(stderr, "+[M::%s] dn::%lu\n", __func__, dn); + + for (k = 0; k < dn; k++) { + a0 = clus.a + (idx[k]>>32); a0n = ((uint32_t)idx[k]) - (idx[k]>>32); + assert(a0n > 0 && ((uint32_t)a0[0]) == k); + + z = &(src[k]); + for (j = 0; j < z->length; j++) { + qn = Get_qn(z->buffer[j]); + tn = Get_tn(z->buffer[j]); + if(tn >= dn) continue; + if(qn > tn) continue; + a1 = clus.a + (idx[tn]>>32); + a1n = ((uint32_t)idx[tn]) - (idx[tn]>>32); + // fprintf(stderr, "+[M::%s] qn::%lu, tn::%lu, a0n::%lu, a1n::%lu\n", __func__, qn, tn, a0n, a1n); + assert(a1n > 0 && ((uint32_t)a1[0]) == tn); + renew_paf0(src, a0, a0n, a1, a1n); + } + + z = &(r_src[k]); + for (j = 0; j < z->length; j++) { + qn = Get_qn(z->buffer[j]); + tn = Get_tn(z->buffer[j]); + if(tn >= dn) continue; + if(qn > tn) continue; + a1 = clus.a + (idx[tn]>>32); + a1n = ((uint32_t)idx[tn]) - (idx[tn]>>32); + // fprintf(stderr, "-[M::%s] qn::%lu, tn::%lu, a0n::%lu, a1n::%lu\n", __func__, qn, tn, a0n, a1n); + assert(a1n > 0 && ((uint32_t)a1[0]) == tn); + renew_paf0(r_src, a0, a0n, a1, a1n); + } + + + + // for (l = k; l < dn; l++) { + // a1 = clus.a + (idx[l]>>32); a1n = ((uint32_t)idx[l]) - (idx[l]>>32); + // fprintf(stderr, "[M::%s] k::%lu, l::%lu, a0n::%lu, a1n::%lu\n", __func__, k, l, a0n, a1n); + // assert(a1n > 0 && ((uint32_t)a1[0]) == l); + // renew_paf0(src, a0, a0n, a1, a1n); + // renew_paf0(r_src, a0, a0n, a1, a1n); + // } + if(a0n > 1) renew_paf1(r_src, a0, a0n, ng->seq[k].len); + } + free(idx); free(clus.a); +} + +void push_scaff_node(ma_hit_t_alloc *src, uint64_t v, uint64_t w, uint64_t ol, asg_t *ng) +{ + uint64_t vl = ng->seq[v>>1].len, wl = ng->seq[w>>1].len; + ma_hit_t arc, arc1; + + arc.qns = (v>>1)<<32; + if(!(v&1)) { + arc.qns += vl - ol; arc.qe = vl; + } else { + arc.qns += 0; arc.qe = ol; + } + arc.tn = w>>1; + if(!(w&1)) { + arc.ts = 0; arc.te = ol; + } else { + arc.ts = wl - ol; arc.te = wl; + } + + arc.rev = (v^w)&1; arc.el = 0; + arc.ml = arc.no_l_indel = arc.bl = 0; + + add_ma_hit_t_alloc(&(src[v>>1]), &arc); + set_reverse_overlap(&arc1, &arc); + add_ma_hit_t_alloc(&(src[w>>1]), &arc1); +} + +asg_t *gen_ng(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **cov, R_to_U *ruI, +uint64_t scaffold_len) +{ + ma_utg_t *u; uint64_t i, v, w, m, h, z, raw_v, raw_w, nocc, nv, vx, wx; int32_t r; + asg_arc_t t, *p; asg_arc_t *av = NULL; uint64_t slen = scaffold_len + (uopt->min_ovlp*2); + asg_ext_t ext; memset(&ext, 0, sizeof(ext)); ext.ext = asg_init(); asg_t *ng = ext.ext; + ext.a_n = sg->n_seq; + ext.cnt.n = ext.cnt.m = ext.a_n; + CALLOC(ext.cnt.a, ext.cnt.n);///count + + ext.idx_a.n = ext.idx_a.m = sg->n_seq; + MALLOC(ext.idx_a.a, ext.idx_a.n); + memset(ext.idx_a.a, -1, sizeof(*(ext.idx_a.a))*ext.idx_a.n);//map + + for (i = 0; i < ug->g->n_seq; ++i) { + u = &(ug->u.a[i]); + for (m = 0; m < u->n; m++) { + v = u->a[m]>>32; + ext.cnt.a[v>>1]++; + if(ext.cnt.a[v>>1] == 1) { + h = v>>1; ext.idx_a.a[h] = v>>1; + } else { + h = ext.idx_a.n; + kv_push(uint64_t, ext.idx_a, (v>>1)); + } + + z = (h<<1)|(v&1); z <<= 32; z |= ((uint32_t)u->a[m]); + u->a[m] = z; + } + } + + for (i = 0; i < sg->n_seq; ++i) { + asg_seq_set(ng, i, sg->seq[i].len, ext.idx_a.a[i]==((uint64_t)-1)?1:0); + ng->seq[i].c = 0; + if(ext.idx_a.a[i]!=((uint64_t)-1)) assert(ext.idx_a.a[i] == i); + ext.idx_a.a[i] = i;///set for delted read + } + for (; i < ext.idx_a.n; i++) { + asg_seq_set(ng, i, sg->seq[ext.idx_a.a[i]].len, 0); + ng->seq[i].c = 0; assert(!(ng->seq[ext.idx_a.a[i]].del)); + } + ng->r_seq = ng->n_seq; + // fprintf(stderr, "+[M::%s] ng->n_seq::%u, ext.idx_a.n::%u\n", + // __func__, (uint32_t)ng->n_seq, (uint32_t)ext.idx_a.n); + assert(ng->n_seq == ext.idx_a.n); + + for (i = 0, nocc = ng->r_seq; i < ug->g->n_seq; ++i) { + u = &(ug->u.a[i]); + for (m = 1; m < u->n; m++) { + v = u->a[m-1]>>32; w = u->a[m]>>32; + raw_v = (ext.idx_a.a[v>>1]<<1)|(v&1); + raw_w = (ext.idx_a.a[w>>1]<<1)|(w&1); + assert((raw_v>>1) < sg->n_seq); + assert((raw_w>>1) < sg->n_seq); + if(gen_spec_edge(sg, uopt, raw_v, raw_w, &t) < 0) { + asg_seq_set(ng, nocc, slen, 0); + ng->seq[nocc].c = 0; + kv_push(uint64_t, ext.idx_a, ((uint64_t)-1)); + nocc++;///a scaffold node + } + } + + v = i<<1; nv = asg_arc_n(ug->g, v); av = asg_arc_a(ug->g, v); + for (m = 0; m < nv; m++) { + if(av[m].del) continue; + w = av[m].v; + vx = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32)); + wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32)); + raw_v = (ext.idx_a.a[vx>>1]<<1)|(vx&1); + raw_w = (ext.idx_a.a[wx>>1]<<1)|(wx&1); + assert((raw_v>>1) < sg->n_seq); + assert((raw_w>>1) < sg->n_seq); + if(vx > wx) continue; + if((vx == wx) && (vx&1)) continue; + if(gen_spec_edge(sg, uopt, raw_v, raw_w, &t) < 0) { + asg_seq_set(ng, nocc, slen, 0); + ng->seq[nocc].c = 0; + kv_push(uint64_t, ext.idx_a, ((uint64_t)-1)); + nocc++;///a scaffold node + } + } + + v = (i<<1)+1; nv = asg_arc_n(ug->g, v); av = asg_arc_a(ug->g, v); + for (m = 0; m < nv; m++) { + if(av[m].del) continue; + w = av[m].v; + vx = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32)); + wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32)); + raw_v = (ext.idx_a.a[vx>>1]<<1)|(vx&1); + raw_w = (ext.idx_a.a[wx>>1]<<1)|(wx&1); + assert((raw_v>>1) < sg->n_seq); + assert((raw_w>>1) < sg->n_seq); + if(vx > wx) continue; + if((vx == wx) && (vx&1)) continue; + if(gen_spec_edge(sg, uopt, raw_v, raw_w, &t) < 0) { + asg_seq_set(ng, nocc, slen, 0); + ng->seq[nocc].c = 0; + kv_push(uint64_t, ext.idx_a, ((uint64_t)-1)); + nocc++;///a scaffold node + } + } + } + // fprintf(stderr, "+[M::%s] ng->n_seq::%u, ext.idx_a.n::%u\n", + // __func__, (uint32_t)ng->n_seq, (uint32_t)ext.idx_a.n); + assert(ng->n_seq == ext.idx_a.n); + CALLOC(ng->seq_vis, (ng->n_seq<<1)); + realloc_rdb(&(R_INF), cov, ruI, ext.idx_a.a, ext.idx_a.n, slen, (char *)"scaf", ng, uopt); + update_paf(R_INF.paf, R_INF.reverse_paf, ext.idx_a.a, ext.idx_a.n, sg->n_seq, ng); + + + for (i = 0, nocc = ng->r_seq; i < ug->g->n_seq; ++i) { + u = &(ug->u.a[i]); + for (m = 1; m < u->n; m++) { + // fprintf(stderr, "[M::%s] i::%lu, m::%lu\n", __func__, i, m); + v = u->a[m-1]>>32; w = u->a[m]>>32; + r = gen_spec_edge(ng, uopt, v, w, &t); + if(r >= 0) { + p = asg_arc_pushp(ng); *p = t; + r = gen_spec_edge(ng, uopt, w^1, v^1, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + } else { + z = nocc<<1; nocc++; + push_scaff_node(R_INF.paf, v, z, uopt->min_ovlp, ng); + push_scaff_node(R_INF.paf, z, w, uopt->min_ovlp, ng); + + r = gen_spec_edge(ng, uopt, v, z, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + r = gen_spec_edge(ng, uopt, z^1, v^1, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + + r = gen_spec_edge(ng, uopt, z, w, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + r = gen_spec_edge(ng, uopt, w^1, z^1, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + } + } + + v = i<<1; nv = asg_arc_n(ug->g, v); av = asg_arc_a(ug->g, v); + for (m = 0; m < nv; m++) { + if(av[m].del) continue; + w = av[m].v; + vx = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32)); + wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32)); + if(vx > wx) continue; + if((vx == wx) && (vx&1)) continue;///it is possible + + r = gen_spec_edge(ng, uopt, vx, wx, &t); + if(r >= 0) { + p = asg_arc_pushp(ng); *p = t; + r = gen_spec_edge(ng, uopt, wx^1, vx^1, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + } else { + z = nocc<<1; nocc++; + push_scaff_node(R_INF.paf, vx, z, uopt->min_ovlp, ng); + push_scaff_node(R_INF.paf, z, wx, uopt->min_ovlp, ng); + + r = gen_spec_edge(ng, uopt, vx, z, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + r = gen_spec_edge(ng, uopt, z^1, vx^1, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + + r = gen_spec_edge(ng, uopt, z, wx, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + r = gen_spec_edge(ng, uopt, wx^1, z^1, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + } + } + + v = (i<<1)+1; nv = asg_arc_n(ug->g, v); av = asg_arc_a(ug->g, v); + for (m = 0; m < nv; m++) { + if(av[m].del) continue; + w = av[m].v; + vx = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32)); + wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32)); + if(vx > wx) continue; + if((vx == wx) && (vx&1)) continue;///it is possible + + r = gen_spec_edge(ng, uopt, vx, wx, &t); + if(r >= 0) { + p = asg_arc_pushp(ng); *p = t; + r = gen_spec_edge(ng, uopt, wx^1, vx^1, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + } else { + z = nocc<<1; nocc++; + push_scaff_node(R_INF.paf, vx, z, uopt->min_ovlp, ng); + push_scaff_node(R_INF.paf, z, wx, uopt->min_ovlp, ng); + + r = gen_spec_edge(ng, uopt, vx, z, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + r = gen_spec_edge(ng, uopt, z^1, vx^1, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + + r = gen_spec_edge(ng, uopt, z, wx, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + r = gen_spec_edge(ng, uopt, wx^1, z^1, &t); + assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + } + } + } + asg_cleanup(ng); ng->r_seq = ng->n_seq; + free(ext.cnt.a); free(ext.idx_a.a); + fprintf(stderr, "[M::%s] nocc::%lu, ng->n_seq::%u, sg->n_seq::%u\n", + __func__, nocc, (uint32_t)ng->n_seq, (uint32_t)sg->n_seq); + assert(nocc == ng->n_seq); + return ng; +} + void gen_ul_trio_graph(ug_opt_t *uopt, ul_resolve_t *uidx, char *o_file) { output_trio_unitig_graph_ul(uopt, uidx, o_file, FATHER); output_trio_unitig_graph_ul(uopt, uidx, o_file, MOTHER); } -void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, +ma_ug_t *ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file) { uint64_t i; uint8_t *r_het = NULL; bubble_type *bub = NULL; ulg_opt_t uu; @@ -14302,12 +15178,11 @@ double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_ // resolve_dip_bub_chains(uidx); // free(r_het); destory_bubbles(bub); free(bub); - print_debug_gfa(sg, uidx->uovl.hybrid_ug, uopt->coverage_cut, "hybrid_ug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1); - print_debug_gfa(sg, uidx->uovl.hybrid_ug, uopt->coverage_cut, "hybrid_ug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0); - print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0); + // print_debug_gfa(sg, uidx->uovl.hybrid_ug, uopt->coverage_cut, "hybrid_ug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1); + // print_debug_gfa(sg, uidx->uovl.hybrid_ug, uopt->coverage_cut, "hybrid_ug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0); + // print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0); - if(is_trio) gen_ul_trio_graph(uopt, uidx, o_file); - - - exit(0); + // if(is_trio) gen_ul_trio_graph(uopt, uidx, o_file); + // exit(0); + return uidx->uovl.hybrid_ug; } \ No newline at end of file diff --git a/gfa_ut.h b/gfa_ut.h index 761cc71..6a14b43 100644 --- a/gfa_ut.h +++ b/gfa_ut.h @@ -23,11 +23,13 @@ void asg_arc_cut_bub_links(asg_t *g, asg64_v *in, float len_rat, float sec_len_r void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float ou_rat, uint32_t is_ou, bub_label_t *b_mask_t); uint32_t asg_cut_large_indel(asg_t *g, asg64_v *in, int32_t max_ext, float ou_rat, uint32_t is_ou); uint32_t asg_cut_semi_circ(asg_t *g, uint32_t lim_len, uint32_t is_clean); -void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, +ma_ug_t *ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file); void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t max_hang, int64_t min_ovlp, int64_t ul_occ); void normalize_gou(asg_t *g); void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd); +asg_t *gen_ng(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **cov, R_to_U *ruI, uint64_t scaffold_len); +void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover); // void print_raw_u2rgfa_seq(all_ul_t *aln, R_to_U* rI, uint32_t is_detail); #endif diff --git a/hic.cpp b/hic.cpp index a1d51ea..33bb101 100644 --- a/hic.cpp +++ b/hic.cpp @@ -16497,7 +16497,7 @@ void optimize_u_trans(kv_u_trans_t *ovlp, kvec_pe_hit* hits, ha_ug_index* idx) } } kt_u_trans_t_idx(&k_trans, idx->ug->g->n_seq); - mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL, 1); + mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL, 1, 0); for (i = m = 0; i < ovlp->n; i++){ x = &(ovlp->a[i]); if(x->del) continue; @@ -16524,7 +16524,7 @@ ha_ug_index* idx, uint64_t test_block_flip, uint64_t n_perturb) __func__, test_block_flip, n_perturb); renew_kv_u_trans(k_trans, link, &sl->hits, &(idx->t_ch->k_trans), idx, bub, s->s, NULL, 0); mc_solve(NULL, NULL, k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, - (bub->round_id == 0? 1 : 0), s->s, 1, NULL, test_block_flip?&(idx->t_ch->k_trans):0, 0); + (bub->round_id == 0? 1 : 0), s->s, 1, NULL, test_block_flip?&(idx->t_ch->k_trans):0, 0, 0); } @@ -16596,6 +16596,9 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o // label_unitigs_sm(s->s, NULL, idx->ug); // goto skip_flipping; // } + // print_debug_gfa(idx->read_g, idx->ug, opt->coverage_cut, "hic.phasing", opt->sources, opt->ruIndex, + // opt->max_hang, opt->min_ovlp, 0, 0, 0); + s = init_ps_t(11, idx->ug->g->n_seq); // debug_round_test(s, 11, &bub, &k_trans, &link, &sl, idx, 1000, 10000); for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++) @@ -16605,8 +16608,11 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o // if(bub.round_id == 0) init_phase(idx, &k_trans, &bub, s); // update_trans_g(idx, &k_trans, &bub); /*******************************for debug************************************/ + // mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, + // (bub.round_id == 0? 1 : 0), s->s, 1, (asm_opt.ar)?(&bub):(NULL), &(idx->t_ch->k_trans), 0, + // (((bub.round_id+1) == bub.n_round)?1:0)); mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, - (bub.round_id == 0? 1 : 0), s->s, 1, /**&bub**/NULL, &(idx->t_ch->k_trans), 0); + (bub.round_id == 0? 1 : 0), s->s, 1, NULL, &(idx->t_ch->k_trans), 0, 0); /*******************************for debug************************************/ label_unitigs_sm(s->s, NULL, idx->ug); @@ -16639,9 +16645,8 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o ///print_hc_links(&link, 0, &hap); // print_kv_u_trans(&k_trans, &link, s->s); - - ///print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, idx->link, idx); - ///print_hits(idx, &sl.hits, fn1); + // print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, NULL/**idx->link**/, idx); + // print_hits(idx, &sl.hits, fn1, fn2); ///print_debug_bubble_graph(&bub, idx->ug, asm_opt.output_file_name); @@ -17369,4 +17374,4 @@ void hic_benchmark(ma_ug_t *ug, asg_t* read_g) hic_short_align_bench(asm_opt.hic_reads[0], asm_opt.hic_reads[1], output_file_name, ug_index); free(output_file_name); -} \ No newline at end of file +} diff --git a/inter.cpp b/inter.cpp index baa7638..1a6b28e 100644 --- a/inter.cpp +++ b/inter.cpp @@ -4876,7 +4876,8 @@ void update_ul_vec_t_ug(const ul_idx_t *uref, ul_vec_t *rch, vec_mg_lchain_t *uc rch->bb.a[rch->bb.a[k].pidx].aidx = k; } // fprintf(stderr, "+ulid->%ld\n", ulid); - uint32_t sp = (uint32_t)-1, ep = (uint32_t)-1; k = l/**a_n - 1**/;///start from the max chain + uint32_t sp = (uint32_t)-1, ep = (uint32_t)-1, ch_n = 0; + k = l/**a_n - 1**/;///start from the max chain for (l = 0; k >= 0; ) { if(sp == (uint32_t)-1 || rch->bb.a[k].qe <= sp) { if(sp != (uint32_t)-1) l += ep - sp; @@ -4886,6 +4887,7 @@ void update_ul_vec_t_ug(const ul_idx_t *uref, ul_vec_t *rch, vec_mg_lchain_t *uc } if(rch->bb.a[k].pidx == (uint32_t)-1) k = -1; else k = rch->bb.a[k].pidx; + ch_n++; } rch->dd = 0; @@ -4897,8 +4899,9 @@ void update_ul_vec_t_ug(const ul_idx_t *uref, ul_vec_t *rch, vec_mg_lchain_t *uc rch->dd = 1; } else if(l < ((int64_t)rch->rlen)*0.001) { rch->dd = 2; + } else if(ch_n < rch->bb.n) {///multiple chain, might be useful for the scaffolding + rch->dd = 3; } - // fprintf(stderr, "[M::%s::] rch->dd::%u, rch->bb.n::%u\n", // __func__, rch->dd, (uint32_t)rch->bb.n); } @@ -8975,9 +8978,12 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // b->num_correct_base += b->correct.corrected_base; // b->num_recorrect_base += b->round2.dumy.corrected_base; - if(UL_INF.a[s->id+i].dd) { - free(s->seq[i]); s->seq[i] = NULL; b->num_correct_base++; + if(UL_INF.a[s->id+i].dd == 1 || UL_INF.a[s->id+i].dd == 2) { + b->num_correct_base++; } + if(UL_INF.a[s->id+i].dd != 3) { + free(s->seq[i]); s->seq[i] = NULL; + } s->hab[tid]->num_read_base++; // fprintf(stderr, "[M::%s] rid:%ld, dd:%u\n", __func__, s->id+i, UL_INF.a[s->id+i].dd); // int64_t mem[6], mem_hab[6]; @@ -9364,9 +9370,9 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac // s->hab[i] = ha_ovec_buf_init(NULL, 0, 0, 1); s->hab[i] = ha_ovec_init(0, 0, 1); } - fprintf(stderr, "[M::%s::Start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); + // fprintf(stderr, "[M::%s::Start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); kt_for(p->n_thread, worker_for_ul_scall_alignment, s, s->n); - fprintf(stderr, "[M::%s::Done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); + // fprintf(stderr, "[M::%s::Done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); ///debug /** uint64_t i; @@ -9408,7 +9414,7 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac p->num_bases += s->num_bases; p->num_corrected_bases += s->num_corrected_bases; p->num_recorrected_bases += s->num_recorrected_bases; - fprintf(stderr, "[M::%s::dump_start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); + // fprintf(stderr, "[M::%s::dump_start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); for (i = 0; i < p->n_thread; ++i) { push_uc_block_t(s->uopt, &(s->ll[i].tk), s->seq, s->len, s->id); free(s->ll[i].tk.a); @@ -9420,7 +9426,7 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac } free(s->seq[i]); } - fprintf(stderr, "[M::%s::dump_done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); + // fprintf(stderr, "[M::%s::dump_done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); /** for (i = 0; i < (uint64_t)s->n; ++i) { ///debug @@ -9495,10 +9501,10 @@ static void *worker_ul_rescall_pipeline(void *data, int step, void *in) // callb for (i = 0; i < p->n_thread; ++i) { s->hab[i] = ha_ovec_init(0, 0, 1); s->buf[i] = mg_tbuf_init(); } - fprintf(stderr, "[M::%s::Start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); + // fprintf(stderr, "[M::%s::Start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); kt_for(p->n_thread, worker_for_ul_rescall_alignment, s, s->n); - fprintf(stderr, "[M::%s::Done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); - get_utepdat_t_mem(s, 1); + // fprintf(stderr, "[M::%s::Done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); + // get_utepdat_t_mem(s, 1); for (i = 0; i < p->n_thread; ++i) { p->num_bases += s->hab[i]->num_read_base; @@ -9514,10 +9520,10 @@ static void *worker_ul_rescall_pipeline(void *data, int step, void *in) // callb return s; } else if (step == 2) { // step 3: dump utepdat_t *s = (utepdat_t*)in; int64_t i, rid; - fprintf(stderr, "[M::%s::dump_start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); + // fprintf(stderr, "[M::%s::dump_start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); for (i = 0; i < s->n; ++i) { rid = s->id + i; - if(UL_INF.a[rid].dd == 0 && p->ucr_s && p->ucr_s->flag == 1) { + if(UL_INF.a[rid].dd == 3 && p->ucr_s && p->ucr_s->flag == 1) {///for the scaffolding assert(s->seq[i]); // if(s->seq[i] == NULL) fprintf(stderr, "[M::%s::]rid->%ld, len->%lu\n", __func__, rid, s->len[i]); ///for debug interval @@ -9526,7 +9532,7 @@ static void *worker_ul_rescall_pipeline(void *data, int step, void *in) // callb // if(UL_INF.a[rid].dd) fprintf(stderr, "rid->%ld\n", rid); free(s->seq[i]); } - fprintf(stderr, "[M::%s::dump_done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); + // fprintf(stderr, "[M::%s::dump_done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); free(s->len); free(s->seq); free(s); } return 0; @@ -14752,7 +14758,7 @@ ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_c mg_idxopt_t opt; uldat_t sl; int32_t cutoff; char* gfa_name = NULL; MALLOC(gfa_name, strlen(asm_opt.output_file_name)+50); - sprintf(gfa_name, "%s.re", asm_opt.output_file_name); + sprintf(gfa_name, "%s.%s", asm_opt.output_file_name, HA_RE_UL_ID); init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); cutoff = REA_ALIGN_CUTOFF; diff --git a/rcut.cpp b/rcut.cpp index b3b49f7..0746c73 100644 --- a/rcut.cpp +++ b/rcut.cpp @@ -634,7 +634,7 @@ mb_g_t *init_mb_g_t(mc_match_t* e, kv_u_trans_t *ref, uint32_t is_sys) for (m = 0; m < n; m++) { tn = ma_y(o[m]); - tb = p->u->idx.a[tn]>>1; + tb = p->u->idx.a[tn]>>1;///tn is the unitig id; tb is the block id if(tb == qb) { continue; @@ -1901,7 +1901,7 @@ uint8_t *lock, bits_p *vis, uint64_t id, mc_bp_res* r) reset_mc_bp_iter(bp, &i, id); memset(vis->a, 0, vis->n); max_bid = max_uid = (uint32_t)-1; - f_bid = i.bid; f_uid = i.uid; + f_bid = i.bid; f_uid = i.uid;///f_bid::bubble id, f_uid::unitig id while (1) { @@ -2902,14 +2902,20 @@ void print_hap_s(int8_t *s, uint32_t sn) } } -void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov) +void dump_debug_phasing(const char* fn, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, +uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref); +void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov, int is_dump) { + if(is_dump) { + dump_debug_phasing(MC_NAME, ta, ug, read_g, f_rate, renew_s, s, is_sys, bub, ref); + } + mc_opt_t opt; mc_opt_init(&opt, asm_opt.n_perturb, asm_opt.f_perturb, asm_opt.seed); mc_g_t *mg = init_mc_g_t(ug, read_g, s, renew_s); update_mc_edges(mg, ovlp, ta, t_ch, f_rate, is_sys); - // fprintf(stderr, "[M::%s:: # edges: %u]\n", __func__, (uint32_t)mg->e->ma.n); + fprintf(stderr, "[M::%s:: # edges: %u]\n", __func__, (uint32_t)mg->e->ma.n); mb_solve_core(&opt, mg, ref, is_sys); ///debug_mc_g_t(mg); @@ -3788,4 +3794,222 @@ void mc_solve_general(kv_u_trans_t *ta, uint32_t un, kv_gg_status *s, uint16_t h destory_mc_gg_t(&mg); if(write_dump) write_mc_gg_dump(ta, un, s, hapN, MC_NAME); -} \ No newline at end of file +} + +void dump_kv_u_trans_t(kv_u_trans_t *z, FILE *fp) +{ + fwrite(&(z->n), sizeof(z->n), 1, fp); + fwrite(z->a, sizeof((*(z->a))), z->n, fp); + fwrite(&z->idx.n, sizeof(z->idx.n), 1, fp); + fwrite(z->idx.a, sizeof((*(z->idx.a))), z->idx.n, fp); +} + +void dump_asg_t(asg_t *g, FILE *fp) +{ + uint32_t tmp, Len; + tmp = g->n_arc; + fwrite(&tmp, sizeof(tmp), 1, fp); + + tmp = g->is_srt; + fwrite(&tmp, sizeof(tmp), 1, fp); + + tmp = g->n_seq; + fwrite(&tmp, sizeof(tmp), 1, fp); + + tmp = g->is_symm; + fwrite(&tmp, sizeof(tmp), 1, fp); + + tmp = g->r_seq; + fwrite(&tmp, sizeof(tmp), 1, fp); + + + // Len = g->n_seq*2; + // fwrite(g->seq_vis, sizeof((*(g->seq_vis))), Len, fp); + + Len = g->n_seq*2; + fwrite(g->idx, sizeof((*(g->idx))), Len, fp); + fwrite(g->arc, sizeof((*(g->arc))), g->n_arc, fp); + fwrite(g->seq, sizeof((*(g->seq))), g->n_seq, fp); +} + +void dump_ma_ug_t(ma_ug_t *ug, FILE *fp) +{ + ma_utg_t *u = NULL; uint32_t t, i; + fwrite(&(ug->u.n), sizeof(ug->u.n), 1, fp); + for (i = 0; i < ug->u.n; i++) { + u = &(ug->u.a[i]); + t = u->len; + fwrite(&t, sizeof(t), 1, fp); + t = u->circ; + fwrite(&t, sizeof(t), 1, fp); + fwrite(&(u->start), sizeof(u->start), 1, fp); + fwrite(&(u->end), sizeof(u->end), 1, fp); + fwrite(&(u->n), sizeof(u->n), 1, fp); + fwrite(u->a, sizeof(uint64_t), u->n, fp); + } + dump_asg_t(ug->g, fp); +} + +void dump_bubble_type(bubble_type* bub, FILE *fp) +{ + fwrite(&(bub->chain_weight.n), sizeof(bub->chain_weight.n), 1, fp); + fwrite(bub->chain_weight.a, sizeof((*(bub->chain_weight.a))), bub->chain_weight.n, fp); + + fwrite(&(bub->list.n), sizeof(bub->list.n), 1, fp); + fwrite(bub->list.a, sizeof((*(bub->list.a))), bub->list.n, fp); + + fwrite(&(bub->num.n), sizeof(bub->num.n), 1, fp); + fwrite(bub->num.a, sizeof((*(bub->num.a))), bub->num.n, fp); + + fwrite(&(bub->pathLen.n), sizeof(bub->pathLen.n), 1, fp); + fwrite(bub->pathLen.a, sizeof((*(bub->pathLen.a))), bub->pathLen.n, fp); + + dump_ma_ug_t(bub->b_ug, fp); + dump_asg_t(bub->b_g, fp); + dump_ma_ug_t(bub->ug, fp); +} + +void dump_debug_phasing(const char* fn, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, +uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref) +{ + fprintf(stderr, "\n[M::%s]\n", __func__); + char *buf = (char*)calloc(strlen(fn) + 50, 1); + sprintf(buf, "%s.hic.dbg.dump.bin", fn); + FILE *fp = fopen(buf, "w"); + + dump_kv_u_trans_t(ta, fp);///ta + dump_ma_ug_t(ug, fp);///ug + dump_asg_t(read_g, fp);///read_g + fwrite(&f_rate, sizeof(f_rate), 1, fp);///f_rate + fwrite(&renew_s, sizeof(renew_s), 1, fp);///renew_s + fwrite(s, sizeof((*s)), ug->g->n_seq, fp);///s + fwrite(&is_sys, sizeof(is_sys), 1, fp);///is_sys + dump_bubble_type(bub, fp);///bub + dump_kv_u_trans_t(ref, fp);///ta + + fclose(fp); + free(buf); +} + +void load_kv_u_trans_t(kv_u_trans_t *z, FILE *fp) +{ + fread(&(z->n), sizeof(z->n), 1, fp); + z->m = z->n; MALLOC(z->a, z->n); + fread(z->a, sizeof((*(z->a))), z->n, fp); + + fread(&z->idx.n, sizeof(z->idx.n), 1, fp); + z->idx.m = z->idx.n; MALLOC(z->idx.a, z->idx.n); + fread(z->idx.a, sizeof((*(z->idx.a))), z->idx.n, fp); +} + +void load_asg_t(asg_t *g, FILE *fp) +{ + uint32_t tmp, Len; + fread(&tmp, sizeof(tmp), 1, fp); + g->n_arc = g->m_arc = tmp; + + fread(&tmp, sizeof(tmp), 1, fp); + g->is_srt = tmp; + + fread(&tmp, sizeof(tmp), 1, fp); + g->n_seq = g->m_seq = tmp; + + fread(&tmp, sizeof(tmp), 1, fp); + g->is_symm = tmp; + + fread(&tmp, sizeof(tmp), 1, fp); + g->r_seq = tmp; + + // Len = g->n_seq*2; + // fwrite(g->seq_vis, sizeof((*(g->seq_vis))), Len, fp); + + Len = g->n_seq*2; + CALLOC(g->idx, Len); fread(g->idx, sizeof((*(g->idx))), Len, fp); + CALLOC(g->arc, g->n_arc); fread(g->arc, sizeof((*(g->arc))), g->n_arc, fp); + CALLOC(g->seq, g->n_seq); fread(g->seq, sizeof((*(g->seq))), g->n_seq, fp); +} + +void load_ma_ug_t(ma_ug_t *ug, FILE *fp) +{ + ma_utg_t *u = NULL; uint32_t t, i; + fread(&(ug->u.n), sizeof(ug->u.n), 1, fp); + ug->u.m = ug->u.n; CALLOC(ug->u.a, ug->u.n); + for (i = 0; i < ug->u.n; i++) { + u = &(ug->u.a[i]); + fread(&t, sizeof(t), 1, fp); u->len = t; + fread(&t, sizeof(t), 1, fp); u->circ = t; + fread(&(u->start), sizeof(u->start), 1, fp); + fread(&(u->end), sizeof(u->end), 1, fp); + fread(&(u->n), sizeof(u->n), 1, fp); + u->m = u->n; MALLOC(u->a, u->n); + fread(u->a, sizeof((*(u->a))), u->n, fp); + } + CALLOC(ug->g, 1); + load_asg_t(ug->g, fp); +} + +void load_bubble_type(bubble_type* bub, FILE *fp) +{ + fread(&(bub->chain_weight.n), sizeof(bub->chain_weight.n), 1, fp); + bub->chain_weight.m = bub->chain_weight.n; CALLOC(bub->chain_weight.a, bub->chain_weight.n); + fread(bub->chain_weight.a, sizeof((*(bub->chain_weight.a))), bub->chain_weight.n, fp); + + fread(&(bub->list.n), sizeof(bub->list.n), 1, fp); + bub->list.m = bub->list.n; CALLOC(bub->list.a, bub->list.n); + fread(bub->list.a, sizeof((*(bub->list.a))), bub->list.n, fp); + + fread(&(bub->num.n), sizeof(bub->num.n), 1, fp); + bub->num.m = bub->num.n; CALLOC(bub->num.a, bub->num.n); + fread(bub->num.a, sizeof((*(bub->num.a))), bub->num.n, fp); + + fread(&(bub->pathLen.n), sizeof(bub->pathLen.n), 1, fp); + bub->pathLen.m = bub->pathLen.n; CALLOC(bub->pathLen.a, bub->pathLen.n); + fread(bub->pathLen.a, sizeof((*(bub->pathLen.a))), bub->pathLen.n, fp); + + CALLOC(bub->b_ug, 1); load_ma_ug_t(bub->b_ug, fp); + CALLOC(bub->b_g, 1); load_asg_t(bub->b_g, fp); + CALLOC(bub->ug, 1); load_ma_ug_t(bub->ug, fp); +} + +void load_debug_phasing(const char* fn, kv_u_trans_t **ta, ma_ug_t **ug, asg_t **read_g, double *f_rate, +uint32_t *renew_s, int8_t **s, uint32_t *is_sys, bubble_type **bub, kv_u_trans_t **ref) +{ + fprintf(stderr, "\n[M::%s]\n", __func__); + char *buf = (char*)calloc(strlen(fn) + 50, 1); + sprintf(buf, "%s.hic.dbg.dump.bin", fn); + FILE *fp = fopen(buf, "r"); + + CALLOC((*ta), 1); load_kv_u_trans_t(*ta, fp);///ta + CALLOC((*ug), 1); load_ma_ug_t(*ug, fp);///ug + CALLOC((*read_g), 1); load_asg_t((*read_g), fp);///read_g + + fread(f_rate, sizeof((*f_rate)), 1, fp);///f_rate + fread(renew_s, sizeof((*renew_s)), 1, fp);///renew_s + CALLOC((*s), (*ug)->g->n_seq); fread(*s, sizeof((*(*s))), (*ug)->g->n_seq, fp);///s + fread(is_sys, sizeof((*is_sys)), 1, fp);///is_sys + + CALLOC((*bub), 1); load_bubble_type(*bub, fp);///bub + CALLOC((*ref), 1); load_kv_u_trans_t(*ref, fp);///ta + + fclose(fp); + free(buf); +} + + +void quick_debug_phasing(const char* fn) +{ + kv_u_trans_t *ta; ma_ug_t *ug; asg_t *read_g; double f_rate; + uint32_t renew_s; int8_t *s; uint32_t is_sys, k; bubble_type *bub; kv_u_trans_t *ref; + load_debug_phasing(fn, &ta, &ug, &read_g, &f_rate, &renew_s, &s, &is_sys, &bub, &ref); + + + + mc_solve(NULL, NULL, ta, ug, read_g, f_rate, NULL, renew_s, s, is_sys, bub, ref, 0, 0); + + for (k = 0; k < ug->g->n_seq; k++) { + fprintf(stderr, "utg%.6dl(len::%u)\n", (int32_t)(k)+1, ug->g->seq[k].len); + } + + + exit(1); +} diff --git a/rcut.h b/rcut.h index 6510cb8..e2baf47 100644 --- a/rcut.h +++ b/rcut.h @@ -118,11 +118,12 @@ static inline double kr_drand_r(uint64_t *x) return u.d - 1.0; } -void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov); +void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov, int is_dump); void debug_mc_g_t(const char* name); void mc_solve_general(kv_u_trans_t *ta, uint32_t un, kv_gg_status *s, uint16_t hapN, uint16_t update_ta, uint16_t write_dump); kv_gg_status *init_mc_gg_status(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint64_t t_cov, uint16_t hapN); void destory_mc_gg_t(mc_gg_t **p); void debug_mc_gg_t(const char* fn, uint32_t update_ta, uint32_t convert_mc_g_t); +void quick_debug_phasing(const char* fn); #endif \ No newline at end of file