diff --git a/Assembly.cpp b/Assembly.cpp index 07c82b4..f24d9c4 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1517,31 +1517,6 @@ void Output_PAF() fprintf(stderr, "PAF has been written.\n"); } -void Output_yak_binning() -{ - fprintf(stderr, "Writing binning to disk ...... \n"); - char* paf_name = (char*)malloc(strlen(asm_opt.output_file_name)+50); - sprintf(paf_name, "%s.hap1.bin.log", asm_opt.output_file_name); - FILE* oh1 = fopen(paf_name, "w"); - sprintf(paf_name, "%s.hap2.bin.log", asm_opt.output_file_name); - FILE* oh2 = fopen(paf_name, "w"); - uint64_t i; - - for (i = 0; i < R_INF.total_reads; i++) { - if(R_INF.trio_flag[i]==FATHER) { - fprintf(oh1, "%.*s\n", (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); - } - if(R_INF.trio_flag[i]==MOTHER) { - fprintf(oh2, "%.*s\n", (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); - } - } - - free(paf_name); - fclose(oh1); fclose(oh2); - fprintf(stderr, "Binning has been written.\n"); -} - - int check_cluster(uint64_t* list, long long listLen, ma_hit_t_alloc* paf, float threshold) { long long i, k; @@ -1785,13 +1760,9 @@ int ha_assemble(void) ha_extract_print_list(&R_INF, asm_opt.extract_iter, asm_opt.extract_list); exit(0); } - // if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN) && !(asm_opt.flag & HA_F_VERBOSE_GFA)) ha_triobin(&asm_opt); - if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt); - // if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt), ovlp_loaded = 2; if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads(); if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF(); if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2; - if (asm_opt.fn_bin_yak[0] && asm_opt.fn_bin_yak[1]) Output_yak_binning(); } if (!ovlp_loaded) { ha_flt_tab = ha_idx = NULL; @@ -1826,6 +1797,7 @@ int ha_assemble(void) ha_triobin(&asm_opt); } if(ovlp_loaded == 2) ovlp_loaded = 0; + ha_opt_update_cov_min(&asm_opt, asm_opt.hom_cov, MIN_N_CHAIN); build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf, R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round, diff --git a/CommandLines.cpp b/CommandLines.cpp index 5cebb2f..3dac9cf 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -52,6 +52,7 @@ static ko_longopt_t long_options[] = { { "ul-tip", ko_required_argument, 338}, { "low-het", ko_no_argument, 339}, { "s-base", ko_required_argument, 340}, + { "bin-only", ko_no_argument, 341}, { 0, 0, 0 } }; @@ -193,7 +194,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->max_ov_diff_final = 0.03; asm_opt->hom_cov = 20; asm_opt->het_cov = -1024; - asm_opt->max_n_chain = 100; + asm_opt->max_n_chain = MIN_N_CHAIN; asm_opt->min_hist_kmer_cnt = 5; asm_opt->load_index_from_disk = 1; asm_opt->write_index_to_disk = 1; @@ -260,6 +261,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->is_read_trans = 1; asm_opt->is_topo_trans = 1; asm_opt->is_bub_trans = 1; + asm_opt->bin_only = 0; } void destory_enzyme(enzyme* f) @@ -304,6 +306,14 @@ void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov) fprintf(stderr, "[M::%s] updated max_n_chain to %d\n", __func__, opt->max_n_chain); } +void ha_opt_update_cov_min(hifiasm_opt_t *opt, int hom_cov, int min_chain) +{ + int max_n_chain = (int)(hom_cov * opt->high_factor + .499); + opt->hom_cov = hom_cov; opt->max_n_chain = max_n_chain; + if(opt->max_n_chain < min_chain) opt->max_n_chain = min_chain; + fprintf(stderr, "[M::%s] updated max_n_chain to %d\n", __func__, opt->max_n_chain); +} + static int check_file(char* name, const char* opt) { if(!name) @@ -782,8 +792,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 340) { asm_opt->trans_base_rate_sec = atof(opt.arg); if(asm_opt->trans_base_rate_sec < 0) asm_opt->is_base_trans = 0; - } else if (c == 'l') - { ///0: disable purge_dup; 1: purge containment; 2: purge overlap + } + else if (c == 341) asm_opt->bin_only = 1; + else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg); } else if (c == 's') asm_opt->purge_simi_rate_l2 = asm_opt->purge_simi_rate_l3 = atof(opt.arg); diff --git a/CommandLines.h b/CommandLines.h index a61bff9..c2348aa 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.18.4-r496" +#define HA_VERSION "0.18.5-r499" #define VERBOSE 0 @@ -21,9 +21,10 @@ #define HA_F_HIGH_HET 0x400 #define HA_F_PARTITION 0x800 #define HA_F_FAST 0x1000 -#define HA_F_USKEW 0x2000 +#define HA_F_USKEW 0x2000 #define HA_MIN_OV_DIFF 0.02 // min sequence divergence in an overlap +#define MIN_N_CHAIN 100 typedef struct{ int *l, n; @@ -135,6 +136,7 @@ typedef struct { uint8_t is_read_trans; uint8_t is_topo_trans; uint8_t is_bub_trans; + uint8_t bin_only; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; @@ -143,6 +145,7 @@ void init_opt(hifiasm_opt_t* asm_opt); void destory_opt(hifiasm_opt_t* asm_opt); void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round); void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov); +void ha_opt_update_cov_min(hifiasm_opt_t *opt, int hom_cov, int min_chain); int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt); double Get_T(void); diff --git a/Overlaps.cpp b/Overlaps.cpp index 11cae2f..0e57aaa 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -63,7 +63,7 @@ KRADIX_SORT_INIT(ha_mzl_t_srt1, ha_mzl_t, ha_mzl_t_key, member_size(ha_mzl_t, x) KSORT_INIT_GENERIC(uint32_t) void reduce_hamming_error_adv(ma_ug_t *iug, asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, -int max_hang, int min_ovlp, long long gap_fuzz, bubble_type* bub); +int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub); typedef struct { uint32_t d, tot, ma, p; @@ -14636,7 +14636,7 @@ long long gap_fuzz, bub_label_t* b_mask_t, ug_opt_t *opt) asg_cleanup(sg); // reduce_hamming_error(sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz); - reduce_hamming_error_adv(NULL, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, NULL); + reduce_hamming_error_adv(NULL, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, opt->ruIndex, NULL); ug_fa = output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang, min_ovlp, rhits?1:0, b_mask_t, NULL, NULL, NULL); @@ -14762,7 +14762,7 @@ long long gap_fuzz, bub_label_t* b_mask_t) asg_cleanup(sg); // reduce_hamming_error(sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz); - reduce_hamming_error_adv(NULL, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, NULL); + reduce_hamming_error_adv(NULL, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, opt.ruIndex, NULL); output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang, min_ovlp, 0, b_mask_t, NULL, NULL, NULL); @@ -18998,12 +18998,82 @@ void write_debug_ma_hit_ts(ma_hit_t_alloc* x, long long n_read, char* read_file_ fprintf(stderr, "ma_hit_ts has been written.\n"); } +void write_yak_binning(char *ou, char *fn_bin_yak1, char *fn_bin_yak2) +{ + fprintf(stderr, "Writing binning to disk ...... \n"); + char* paf_name = (char*)malloc(strlen(ou)+50); + sprintf(paf_name, "%s.hap1.phase.bin", ou); + FILE* oh1 = fopen(paf_name, "w"); + sprintf(paf_name, "%s.hap2.phase.bin", ou); + FILE* oh2 = fopen(paf_name, "w"); + uint64_t i, s1, s2; + + for (i = s1 = s2 = 0; i < R_INF.total_reads; i++) { + if(R_INF.trio_flag[i]==FATHER) { + if(s1) { + fprintf(oh1, "%.*s\n", (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + } else { + fprintf(oh1, "%.*s\t%s\n", (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), fn_bin_yak1); + } + s1 = 1; + } + if(R_INF.trio_flag[i]==MOTHER) { + if(s2) { + fprintf(oh2, "%.*s\n", (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + } else { + fprintf(oh2, "%.*s\t%s\n", (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), fn_bin_yak2); + } + s2 = 1; + } + } + free(paf_name); + fclose(oh1); fclose(oh2); + fprintf(stderr, "Binning has been written.\n"); +} + +uint32_t load_yak_binning(hifiasm_opt_t *opt, char *ou)///asm_opt +{ + char* paf_name = (char*)malloc(strlen(ou)+50); + uint32_t len[2] = {0}; char *lst0, *lst1, *yak0, *yak1; + lst0 = lst1 = yak0 = yak1 = NULL; + + sprintf(paf_name, "%s.hap1.phase.bin", ou); len[0] = strlen(paf_name)+1; + if(!test_yak_binning(paf_name, opt->fn_bin_yak[0])) { + free(paf_name); return 0; + } + + sprintf(paf_name, "%s.hap2.phase.bin", ou); len[1] = strlen(paf_name)+1; + if(!test_yak_binning(paf_name, opt->fn_bin_yak[1])) { + free(paf_name); return 0; + } + free(paf_name); paf_name = NULL; + + lst0 = opt->fn_bin_list[0]; MALLOC(opt->fn_bin_list[0], len[0]); + sprintf(opt->fn_bin_list[0], "%s.hap1.phase.bin", ou); + + lst1 = opt->fn_bin_list[1]; MALLOC(opt->fn_bin_list[1], len[1]); + sprintf(opt->fn_bin_list[1], "%s.hap2.phase.bin", ou); + + memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads * sizeof(uint8_t)); + yak0 = opt->fn_bin_yak[0]; opt->fn_bin_yak[0] = NULL; + yak1 = opt->fn_bin_yak[1]; opt->fn_bin_yak[1] = NULL; + ha_triobin(opt); + opt->fn_bin_yak[0] = yak0; opt->fn_bin_yak[1] = yak1; + free(opt->fn_bin_list[0]); opt->fn_bin_list[0] = lst0; + free(opt->fn_bin_list[1]); opt->fn_bin_list[1] = lst1; + return 1; +} + void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, All_reads *RNF, char* output_file_name) { char* gfa_name = (char*)malloc(strlen(output_file_name)+25); sprintf(gfa_name, "%s.ec", output_file_name); write_All_reads(RNF, gfa_name); + if((ha_opt_triobin(&asm_opt)) && (asm_opt.fn_bin_yak[0] && asm_opt.fn_bin_yak[1])) { + write_yak_binning(asm_opt.output_file_name, asm_opt.fn_bin_yak[0], asm_opt.fn_bin_yak[1]); + } + sprintf(gfa_name, "%s.ovlp.source", output_file_name); write_ma_hit_ts(sources, RNF->total_reads, gfa_name); @@ -19012,6 +19082,7 @@ void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sou free(gfa_name); fprintf(stderr, "bin files have been written.\n"); + if(asm_opt.bin_only) exit(1); } int load_debug_graph(asg_t** sg, ma_hit_t_alloc** sources, ma_sub_t** coverage_cut, @@ -19024,6 +19095,15 @@ int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_s free(gfa_name); return 0; } + + if(ha_opt_triobin(&asm_opt)) { + if(!((asm_opt.fn_bin_yak[0]) && (asm_opt.fn_bin_yak[1]) && + (load_yak_binning(&asm_opt, asm_opt.output_file_name)))) { + ha_triobin(&asm_opt); + write_yak_binning(asm_opt.output_file_name, asm_opt.fn_bin_yak[0], asm_opt.fn_bin_yak[1]); + } + } + if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_debug_graph(NULL, NULL, NULL, output_file_name, NULL, NULL, NULL)) { (*sources) = NULL; @@ -20045,7 +20125,7 @@ buf_t *b, uint64_t tLen, uint64_t vis_f, asg_t *res, asg64_v *sv) for (k_j = 0; k_j < nsu->n; k_j++) { for (k_v = 0; k_v < 2; k_v++) { v = ((nsu->a[k_j]>>33)<<1) + k_v; - if(vis_r_flag[v] != vis_f) break; + if(vis_r_flag[v] != vis_f) continue;; x = &(src[v>>1]); za = asg_arc_a(sg, v); zn = asg_arc_n(sg, v); @@ -20100,7 +20180,10 @@ int asg_arc_del_trans_aux(asg_t *g, asg_t *aux, uint8_t *mark, int fuzz) if(kv == 0) continue; ///av[nv-1] is longest out-dege - L = MAX(asg_arc_len(av0[nv0-1]), asg_arc_len(av1[nv1-1])) + fuzz; + L = 0; + if(nv0) L = asg_arc_len(av0[nv0-1]); + if(nv1 && L < asg_arc_len(av1[nv1-1])) L = asg_arc_len(av1[nv1-1]); + L += fuzz; for (i = 0; i < nv0; ++i) { //w is an out-node of v w = av0[i].v; if (mark[w] != 1) continue; ///w has already been reduced @@ -20128,16 +20211,27 @@ int asg_arc_del_trans_aux(asg_t *g, asg_t *aux, uint8_t *mark, int fuzz) return n_reduced; } - +/** +#define ba_fetch(pa, ca, v0, v) (((v)==(v0))?(pa)[(v)]:(ca)[(v)]) uint64_t rd_hm_bub(asg_t *g, asg_t *ref, uint32_t v0, uint64_t max_dist, buf_t *b) { uint32_t i0, i1, n_pending = 0, is_first = 1, n_tips, tip_end; uint64_t n_pop = 0; - uint32_t v, w, d, nv0, nv1, l, x, i; asg_arc_t *av0, *av1; binfo_t *t; + uint32_t v, w, d, nv0, nv1, l, x, i; asg_arc_t *av0, *av1; binfo_t *t, *pa, *ca; if (g->seq[v0>>1].del) return 0; // already deleted nv0 = g?get_real_length(g, v0, NULL):0; nv1 = ref?get_real_length(ref, v0, NULL):0; - if((nv0+nv1)<2) return 0; - b->a[v0].c = b->a[v0].d = b->a[v0].m = b->a[v0].nc = b->a[v0].np = 0; + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, "[M::%s] v0>>1::%u(v0&1::%u), nv0::%u, nv1::%u\n", __func__, v0>>1, v0&1, + // nv0, nv1); + // } + if((nv0+nv1) < 2) return 0; + pa = b->a; ca = b->a + (g->n_seq<<1); + // if(v0 >= (ref->n_seq<<1)) { + // fprintf(stderr, "[M::%s] v0::%u, ref->n_seq::%u\n", __func__, v0, ref->n_seq); + // } + b->S.n = b->T.n = b->b.n = b->e.n = 0; + v = v0; t = &(ba_fetch(pa, ca, v0, v)); + t->c = t->d = t->m = t->nc = t->np = 0; ///b->S is the nodes with all incoming edges visited kv_push(uint32_t, b->S, v0); n_tips = 0; tip_end = (uint32_t)-1; @@ -20145,7 +20239,7 @@ uint64_t rd_hm_bub(asg_t *g, asg_t *ref, uint32_t v0, uint64_t max_dist, buf_t * do { ///v is a node that all incoming edges have been visited ///d is the distance from v0 to v - v = kv_pop(b->S); d = b->a[v].d; + v = kv_pop(b->S); d = (ba_fetch(pa, ca, v0, v)).d; nv0 = 0; av0 = NULL; nv1 = 0; av1 = NULL; i0 = i1 = 0; if(g) { nv0 = asg_arc_n(g, v); av0 = asg_arc_a(g, v); @@ -20153,10 +20247,17 @@ uint64_t rd_hm_bub(asg_t *g, asg_t *ref, uint32_t v0, uint64_t max_dist, buf_t * if(ref) { nv1 = asg_arc_n(ref, v); av1 = asg_arc_a(ref, v); } + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, "[M::%s] v0>>1::%u(v0&1::%u), v>>1::%u(v&1::%u)\n", __func__, v0>>1, v0&1, v>>1, v&1); + // } ///all out-edges of v for (i0 = 0; i0 < nv0; ++i0) { // loop through v's neighbors if (av0[i0].del) continue; - w = av0[i0].v; l = (uint32_t)av0[i0].ul; t = &b->a[w]; + w = av0[i0].v; l = (uint32_t)av0[i0].ul; t = &((ba_fetch(pa, ca, v0, w))); + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, "[M::%s] v0>>1::%u(v0&1::%u), w>>1::%u(w&1::%u)\n", + // __func__, v0>>1, v0&1, w>>1, w&1); + // } if ((w>>1) == (v0>>1)) goto pop_rd_hm_bub; if(is_first) l = 0; if (d + l > max_dist) break; // too far @@ -20170,12 +20271,27 @@ uint64_t rd_hm_bub(asg_t *g, asg_t *ref, uint32_t v0, uint64_t max_dist, buf_t * ///the shortest path if (d + l < t->d) t->d = d + l; // update dist } - + // if((v0 == 4386) && (!ref) && ((w>>1) == 2597)) { + // fprintf(stderr, "[M::%s] w>>1::%u(w&1::%u), t->r::%u\n", + // __func__, w>>1, w&1, t->r); + // } //if all incoming edges of w have visited //push it to b->S if (--(t->r) == 0) { x = (g?get_real_length(g, w, NULL):0)+(ref?get_real_length(ref, w, NULL):0); + // if((v0 == 4386) && (!ref) && ((w>>1) == 2597)) { + // fprintf(stderr, "[M::%s] w>>1::%u(w&1::%u), t->r::%u, x::%u, b->S.n::%u\n", + // __func__, w>>1, w&1, t->r, x, (uint32_t)b->S.n); + // } if(x > 0) { + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, "+[M::%s] v0>>1::%u(v0&1::%u), push_w>>1::%u(w&1::%u)\n", + // __func__, v0>>1, v0&1, w>>1, w&1); + // uint32_t m; + // for (m = 0; m < b->S.n; m++) { + // fprintf(stderr, "+[M::%s] m>>1::%u(m&1::%u)\n", __func__, b->S.a[m]>>1, b->S.a[m]&1); + // } + // } kv_push(uint32_t, b->S, w); } else { ///at most one tip @@ -20188,7 +20304,7 @@ uint64_t rd_hm_bub(asg_t *g, asg_t *ref, uint32_t v0, uint64_t max_dist, buf_t * if (i0 >= nv0) { for (i1 = 0; i1 < nv1; ++i1) { // loop through v's neighbors if (av1[i1].del) continue; - w = av1[i1].v; l = (uint32_t)av1[i1].ul; t = &b->a[w]; + w = av1[i1].v; l = (uint32_t)av1[i1].ul; t = &((ba_fetch(pa, ca, v0, w))); if ((w>>1) == (v0>>1)) goto pop_rd_hm_bub; if(is_first) l = 0; if (d + l > max_dist) break; // too far @@ -20208,6 +20324,10 @@ uint64_t rd_hm_bub(asg_t *g, asg_t *ref, uint32_t v0, uint64_t max_dist, buf_t * if (--(t->r) == 0) { x = (g?get_real_length(g, w, NULL):0)+(ref?get_real_length(ref, w, NULL):0); if(x > 0) { + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, "-[M::%s] v0>>1::%u(v0&1::%u), push_w>>1::%u(w&1::%u)\n", + // __func__, v0>>1, v0&1, w>>1, w&1); + // } kv_push(uint32_t, b->S, w); } else { ///at most one tip @@ -20223,19 +20343,27 @@ uint64_t rd_hm_bub(asg_t *g, asg_t *ref, uint32_t v0, uint64_t max_dist, buf_t * //if found a tip if(n_tips == 1) { if(tip_end != (uint32_t)-1 && n_pending == 0 && b->S.n == 0) { + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, ">[M::%s] v0>>1::%u(v0&1::%u), push_w>>1::%u(w&1::%u)\n", + // __func__, v0>>1, v0&1, tip_end>>1, tip_end&1); + // } kv_push(uint32_t, b->S, tip_end); break; } else { goto pop_rd_hm_bub; } } + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, "[M::%s] v0>>1::%u(v0&1::%u), b->S.n::%u, n_pending::%u\n", + // __func__, v0>>1, v0&1, (uint32_t)b->S.n, n_pending); + // } ///if i < nv, that means (d + l > max_dist) if (i0 < nv0 || i1 < nv1 || b->S.n == 0) goto pop_rd_hm_bub; } while (b->S.n > 1 || n_pending); n_pop = 1; pop_rd_hm_bub: for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices - t = &b->a[b->b.a[i]]; + t = &((ba_fetch(pa, ca, v0, b->b.a[i]))); t->s = t->c = t->d = t->m = t->nc = t->np = 0; } return n_pop; @@ -20279,23 +20407,25 @@ uint64_t rd_hm_drop(asg_t *g, asg_t *ref, uint32_t v0, uint32_t v1, double cutof uint32_t i1, ncut = 0; uint32_t v, w, nv1, i; asg_arc_t *av1; if (g->seq[v0>>1].del) return 0; // already deleted + b->S.n = 0; binfo_t *pa, *ca; + pa = b->a; ca = b->a + (g->n_seq<<1); ///b->S is the nodes with all incoming edges visited kv_push(uint32_t, b->S, v0); while(b->S.n) { v = kv_pop(b->S); - if(b->a[v].s) continue; - b->a[v].s = 1; + if((ba_fetch(pa, ca, v0, v)).s) continue; + (ba_fetch(pa, ca, v0, v)).s = 1; kv_push(uint32_t, b->b, v); // save it for revert nv1 = asg_arc_n(ref, v); av1 = asg_arc_a(ref, v); for (i1 = 0; i1 < nv1; ++i1) { // loop through v's neighbors if (av1[i1].del) continue; w = av1[i1].v; - if(b->a[w].s || w == v1) continue; + if((ba_fetch(pa, ca, v0, w)).s || w == v1) continue; kv_push(uint32_t, b->S, w); } } for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices - v = b->b.a[i]; b->a[b->b.a[i]].s = 0; + v = b->b.a[i]; (ba_fetch(pa, ca, v0, b->b.a[i])).s = 0; if(v == v0 || v == v1) continue; ncut += rd_hm_drop0(g, ref, v, cutoff); ncut += rd_hm_drop0(g, ref, v^1, cutoff); @@ -20311,9 +20441,20 @@ static void rd_hamming_symm(void *data, long i, int tid) // callback for kt_for( uint32_t st = s->rr->a[i]>>32, ed = (uint32_t)(s->rr->a[i]), p, k, ncut; double step = 0.2, cuttoff; uint64_t max_dist = s->max_dist; p = rd_hm_bub(s->g, s->ref, st, max_dist, b); - if(p) return; + if(p) { + // if(!(b->S.a[0] == (ed^1))) { + // fprintf(stderr, "[M::%s] st>>1::%u(st&1::%u), ed>>1::%u(ed&1::%u), S[0]>>1::%u(S[0]&1::%u), max_dist::%lu\n", + // __func__, st>>1, st&1, ed>>1, ed&1, b->S.a[0]>>1, b->S.a[0]&1, max_dist); + // } + assert(b->S.a[0] == (ed^1)); + return; + } ///recalculate max_dist p = rd_hm_bub(s->ref, NULL, st, max_dist, b); + // if(!p) { + // fprintf(stderr, "[M::%s] st>>1::%u(st&1::%u), ed>>1::%u(ed&1::%u), max_dist::%lu\n", + // __func__, st>>1, st&1, ed>>1, ed&1, max_dist); + // } assert(p); assert(b->S.a[0] == (ed^1)); for (k = max_dist = 0; k < b->b.n; ++k) { if(b->b.a[k]==st || b->b.a[k]==b->S.a[0]) continue; @@ -20322,21 +20463,30 @@ static void rd_hamming_symm(void *data, long i, int tid) // callback for kt_for( max_dist += s->ref->seq[st>>1].len; max_dist += s->ref->seq[b->S.a[0]>>1].len; p = rd_hm_bub(s->g, s->ref, st, max_dist, b); - if(p) return; + if(p) { + assert(b->S.a[0] == (ed^1)); + return; + } for (cuttoff = step; cuttoff < 1.0; cuttoff += step) { ncut = rd_hm_drop(s->g, s->ref, st, ed^1, cuttoff, b); p = rd_hm_bub(s->g, s->ref, st, max_dist, b); - if(p) return; + if(p) { + assert(b->S.a[0] == (ed^1)); + return; + } if(!ncut) break; } rd_hm_drop(s->g, s->ref, st, ed^1, 1024, b); p = rd_hm_bub(s->g, s->ref, st, max_dist, b); - assert(p); + if(p) { + assert(b->S.a[0] == (ed^1)); + return; + } } void reduce_hamming_error_adv(ma_ug_t *iug, asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, -int max_hang, int min_ovlp, long long gap_fuzz, bubble_type* bub) +int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub) { double index_time = yak_realtime(); ma_ug_t *ug = NULL; rd_hamming_t aux_t; memset((&aux_t), 0, sizeof(aux_t)); @@ -20397,7 +20547,7 @@ int max_hang, int min_ovlp, long long gap_fuzz, bubble_type* bub) } } } - free(vis_flag); free(bs_flag); if(!iug) ma_ug_destroy(ug); + free(bs_flag); if(!iug) ma_ug_destroy(ug); if(sv.n > 0) { ig->n_seq = ig->m_seq = sg->n_seq; @@ -20405,8 +20555,10 @@ int max_hang, int min_ovlp, long long gap_fuzz, bubble_type* bub) memcpy(ig->seq, sg->seq, (sizeof((*(ig->seq)))*ig->n_seq)); asg_cleanup(ig); asg_arc_del_trans_aux(ig, sg, vis_flag, gap_fuzz); aux_t.n_thread = asm_opt.thread_num; CALLOC(aux_t.a, aux_t.n_thread); + REALLOC(b.a, (ig->n_seq<<2)); memset(b.a, 0, sizeof((*(b.a)))*(ig->n_seq<<2)); for (i = 0; i < aux_t.n_thread; i++) aux_t.a[i].a = b.a; - aux_t.g = ig; aux_t.ref = sg; aux_t.rr = &sv; + aux_t.g = ig; aux_t.ref = sg; aux_t.rr = &sv; aux_t.max_dist = max_dist; + // print_debug_gfa(ug, sg, coverage_cut, "debug_hamming", sources, ru); kt_for(aux_t.n_thread, rd_hamming_symm, &aux_t, aux_t.rr->n);///all ul + ug for (i = 0; i < aux_t.n_thread; i++) { free(aux_t.a[i].S.a); free(aux_t.a[i].T.a); @@ -20414,7 +20566,7 @@ int max_hang, int min_ovlp, long long gap_fuzz, bubble_type* bub) } free(aux_t.a); } - free(sv.a); + free(sv.a); free(vis_flag); for (i = n_pop = 0; i < ig->n_arc; i++) { if(ig->arc[i].del) continue; @@ -20432,8 +20584,421 @@ int max_hang, int min_ovlp, long long gap_fuzz, bubble_type* bub) fprintf(stderr, "[M::%s::%.3f] # inserted edges: %u, # fixed bubbles: %u\n", __func__, yak_realtime() - index_time, sg->n_arc - n_arc_0, fix_bub); } +**/ +uint64_t rd_hm_bub(asg_t *g, asg_t *ref, uint32_t v0, uint64_t max_dist, buf_t *b) +{ + uint32_t i0, i1, n_pending = 0, is_first = 1, n_tips, tip_end; uint64_t n_pop = 0; + uint32_t v, w, d, nv0, nv1, l, x, i; asg_arc_t *av0, *av1; binfo_t *t; + if (g->seq[v0>>1].del) return 0; // already deleted + nv0 = g?get_real_length(g, v0, NULL):0; + nv1 = ref?get_real_length(ref, v0, NULL):0; + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, "[M::%s] v0>>1::%u(v0&1::%u), nv0::%u, nv1::%u\n", __func__, v0>>1, v0&1, + // nv0, nv1); + // } + if((nv0+nv1)<2) return 0; + // if(v0 >= (ref->n_seq<<1)) { + // fprintf(stderr, "[M::%s] v0::%u, ref->n_seq::%u\n", __func__, v0, ref->n_seq); + // } + b->S.n = b->T.n = b->b.n = b->e.n = 0; + b->a[v0].c = b->a[v0].d = b->a[v0].m = b->a[v0].nc = b->a[v0].np = 0; + ///b->S is the nodes with all incoming edges visited + kv_push(uint32_t, b->S, v0); + n_tips = 0; tip_end = (uint32_t)-1; + do { + ///v is a node that all incoming edges have been visited + ///d is the distance from v0 to v + v = kv_pop(b->S); d = b->a[v].d; + nv0 = 0; av0 = NULL; nv1 = 0; av1 = NULL; i0 = i1 = 0; + if(g) { + nv0 = asg_arc_n(g, v); av0 = asg_arc_a(g, v); + } + if(ref) { + nv1 = asg_arc_n(ref, v); av1 = asg_arc_a(ref, v); + } + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, "[M::%s] v0>>1::%u(v0&1::%u), v>>1::%u(v&1::%u)\n", __func__, v0>>1, v0&1, v>>1, v&1); + // } + ///all out-edges of v + for (i0 = 0; i0 < nv0; ++i0) { // loop through v's neighbors + if (av0[i0].del) continue; + w = av0[i0].v; l = (uint32_t)av0[i0].ul; t = &b->a[w]; + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, "[M::%s] v0>>1::%u(v0&1::%u), w>>1::%u(w&1::%u)\n", + // __func__, v0>>1, v0&1, w>>1, w&1); + // } + if ((w>>1) == (v0>>1)) goto pop_rd_hm_bub; + if(is_first) l = 0; + if (d + l > max_dist) break; // too far + ///unvisited node + if (t->s == 0) { // this vertex has never been visited + kv_push(uint32_t, b->b, w); // save it for revert + t->p = v, t->s = 1, t->d = d + l; + t->r = (g?get_real_length(g, w^1, NULL):0)+(ref?get_real_length(ref, w^1, NULL):0); + ++n_pending; + } else { // visited before + ///the shortest path + if (d + l < t->d) t->d = d + l; // update dist + } + // if((v0 == 4386) && (!ref) && ((w>>1) == 2597)) { + // fprintf(stderr, "[M::%s] w>>1::%u(w&1::%u), t->r::%u\n", + // __func__, w>>1, w&1, t->r); + // } + //if all incoming edges of w have visited + //push it to b->S + if (--(t->r) == 0) { + x = (g?get_real_length(g, w, NULL):0)+(ref?get_real_length(ref, w, NULL):0); + // if((v0 == 4386) && (!ref) && ((w>>1) == 2597)) { + // fprintf(stderr, "[M::%s] w>>1::%u(w&1::%u), t->r::%u, x::%u, b->S.n::%u\n", + // __func__, w>>1, w&1, t->r, x, (uint32_t)b->S.n); + // } + if(x > 0) { + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, "+[M::%s] v0>>1::%u(v0&1::%u), push_w>>1::%u(w&1::%u)\n", + // __func__, v0>>1, v0&1, w>>1, w&1); + // uint32_t m; + // for (m = 0; m < b->S.n; m++) { + // fprintf(stderr, "+[M::%s] m>>1::%u(m&1::%u)\n", __func__, b->S.a[m]>>1, b->S.a[m]&1); + // } + // } + kv_push(uint32_t, b->S, w); + } else { + ///at most one tip + if(n_tips != 0) goto pop_rd_hm_bub; + n_tips++; tip_end = w; + } + --n_pending; + } + } + if (i0 >= nv0) { + for (i1 = 0; i1 < nv1; ++i1) { // loop through v's neighbors + if (av1[i1].del) continue; + w = av1[i1].v; l = (uint32_t)av1[i1].ul; t = &b->a[w]; + if ((w>>1) == (v0>>1)) goto pop_rd_hm_bub; + if(is_first) l = 0; + if (d + l > max_dist) break; // too far + ///unvisited node + if (t->s == 0) { // this vertex has never been visited + kv_push(uint32_t, b->b, w); // save it for revert + t->p = v, t->s = 1, t->d = d + l; + t->r = (g?get_real_length(g, w^1, NULL):0)+(ref?get_real_length(ref, w^1, NULL):0); + ++n_pending; + } else { // visited before + ///the shortest path + if (d + l < t->d) t->d = d + l; // update dist + } + + //if all incoming edges of w have visited + //push it to b->S + if (--(t->r) == 0) { + x = (g?get_real_length(g, w, NULL):0)+(ref?get_real_length(ref, w, NULL):0); + if(x > 0) { + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, "-[M::%s] v0>>1::%u(v0&1::%u), push_w>>1::%u(w&1::%u)\n", + // __func__, v0>>1, v0&1, w>>1, w&1); + // } + kv_push(uint32_t, b->S, w); + } else { + ///at most one tip + if(n_tips != 0) goto pop_rd_hm_bub; + n_tips++; tip_end = w; + } + --n_pending; + } + } + } + + is_first = 0; + //if found a tip + if(n_tips == 1) { + if(tip_end != (uint32_t)-1 && n_pending == 0 && b->S.n == 0) { + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, ">[M::%s] v0>>1::%u(v0&1::%u), push_w>>1::%u(w&1::%u)\n", + // __func__, v0>>1, v0&1, tip_end>>1, tip_end&1); + // } + kv_push(uint32_t, b->S, tip_end); + break; + } else { + goto pop_rd_hm_bub; + } + } + // if((v0 == 4386) && (!ref)) { + // fprintf(stderr, "[M::%s] v0>>1::%u(v0&1::%u), b->S.n::%u, n_pending::%u\n", + // __func__, v0>>1, v0&1, (uint32_t)b->S.n, n_pending); + // } + ///if i < nv, that means (d + l > max_dist) + if (i0 < nv0 || i1 < nv1 || b->S.n == 0) goto pop_rd_hm_bub; + } while (b->S.n > 1 || n_pending); + n_pop = 1; + pop_rd_hm_bub: + for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices + t = &b->a[b->b.a[i]]; + t->s = t->c = t->d = t->m = t->nc = t->np = 0; + } + return n_pop; +} + +uint64_t rd_hm_drop0(asg_t *g, asg_t *ref, uint32_t v, double cutoff) +{ + uint32_t nv0, nv1, mol = 0, i0, i1, ncut = 0; asg_arc_t *av0, *av1; + nv0 = asg_arc_n(g, v); av0 = asg_arc_a(g, v); + nv1 = asg_arc_n(ref, v); av1 = asg_arc_a(ref, v); + if(cutoff < 1) { + for (i0 = 0; i0 < nv0; ++i0) { // loop through v's neighbors + if (av0[i0].del) continue; + if(mol < av0[i0].ol) mol = av0[i0].ol; + } + for (i1 = 0; i1 < nv1; ++i1) { // loop through v's neighbors + if (av1[i1].del) continue; + if(mol < av1[i1].ol) mol = av1[i1].ol; + } + if(mol > 0) { + for (i0 = 0; i0 < nv0; ++i0) { // loop through v's neighbors + if (av0[i0].del) continue; + if(av0[i0].ol < (mol*cutoff)) { + av0[i0].del = 1; asg_arc_del(g, av0[i0].v^1, (av0[i0].ul>>32)^1, 1); + ncut++; + } + } + } + } else { + for (i0 = 0; i0 < nv0; ++i0) { // loop through v's neighbors + if (av0[i0].del) continue; + av0[i0].del = 1; asg_arc_del(g, av0[i0].v^1, (av0[i0].ul>>32)^1, 1); + ncut++; + } + } + return ncut; +} + +uint64_t rd_hm_drop(asg_t *g, asg_t *ref, uint32_t v0, uint32_t v1, double cutoff, buf_t *b) +{ + uint32_t i1, ncut = 0; + uint32_t v, w, nv1, i; asg_arc_t *av1; + if (g->seq[v0>>1].del) return 0; // already deleted + ///b->S is the nodes with all incoming edges visited + b->S.n = 0; + kv_push(uint32_t, b->S, v0); + while(b->S.n) { + v = kv_pop(b->S); + if(b->a[v].s) continue; + b->a[v].s = 1; + kv_push(uint32_t, b->b, v); // save it for revert + nv1 = asg_arc_n(ref, v); av1 = asg_arc_a(ref, v); + for (i1 = 0; i1 < nv1; ++i1) { // loop through v's neighbors + if (av1[i1].del) continue; + w = av1[i1].v; + if(b->a[w].s || w == v1) continue; + kv_push(uint32_t, b->S, w); + } + } + for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices + v = b->b.a[i]; b->a[b->b.a[i]].s = 0; + if(v == v0 || v == v1) continue; + ncut += rd_hm_drop0(g, ref, v, cutoff); + ncut += rd_hm_drop0(g, ref, v^1, cutoff); + } + ncut += rd_hm_drop0(g, ref, v0, cutoff); + ncut += rd_hm_drop0(g, ref, v1^1, cutoff); + return ncut; +} + +void rd_hamming_symm(void *data, long i, int tid) // callback for kt_for() +{ + rd_hamming_t *s = (rd_hamming_t *)data; buf_t *b = &(s->a[tid]); + uint32_t st = s->rr->a[i]>>32, ed = (uint32_t)(s->rr->a[i]), p, k, ncut; + double step = 0.2, cuttoff; uint64_t max_dist = s->max_dist; + p = rd_hm_bub(s->g, s->ref, st, max_dist, b); + if(p) { + assert(b->S.a[0] == (ed^1)); + return; + } + ///recalculate max_dist + p = rd_hm_bub(s->ref, NULL, st, max_dist, b); + // if(!p) { + // fprintf(stderr, "[M::%s] st>>1::%u(st&1::%u), ed>>1::%u(ed&1::%u), max_dist::%lu\n", + // __func__, st>>1, st&1, ed>>1, ed&1, max_dist); + // } + assert(p); assert(b->S.a[0] == (ed^1)); + for (k = max_dist = 0; k < b->b.n; ++k) { + if(b->b.a[k]==st || b->b.a[k]==b->S.a[0]) continue; + max_dist += s->ref->seq[b->b.a[k]>>1].len; + } + max_dist += s->ref->seq[st>>1].len; + max_dist += s->ref->seq[b->S.a[0]>>1].len; + p = rd_hm_bub(s->g, s->ref, st, max_dist, b); + if(p) { + assert(b->S.a[0] == (ed^1)); + return; + } + + for (cuttoff = step; cuttoff < 1.0; cuttoff += step) { + ncut = rd_hm_drop(s->g, s->ref, st, ed^1, cuttoff, b); + p = rd_hm_bub(s->g, s->ref, st, max_dist, b); + if(p) { + assert(b->S.a[0] == (ed^1)); + return; + } + if(!ncut) break; + } + rd_hm_drop(s->g, s->ref, st, ed^1, 1024, b); + p = rd_hm_bub(s->g, s->ref, st, max_dist, b); + if(p) { + assert(b->S.a[0] == (ed^1)); + return; + } +} + +void rd_hamming_symm_simple(rd_hamming_t *s, uint32_t st, uint32_t ed) // callback for kt_for() +{ + buf_t *b = &(s->a[0]); double step = 0.2, cuttoff; uint64_t max_dist = s->max_dist; + // uint32_t st = s->rr->a[i]>>32, ed = (uint32_t)(s->rr->a[i]); + uint32_t p, k, ncut; + p = rd_hm_bub(s->g, s->ref, st, max_dist, b); + if(p) { + assert(b->S.a[0] == (ed^1)); + return; + } + ///recalculate max_dist + p = rd_hm_bub(s->ref, NULL, st, max_dist, b); + // if(!p) { + // fprintf(stderr, "[M::%s] st>>1::%u(st&1::%u), ed>>1::%u(ed&1::%u), max_dist::%lu\n", + // __func__, st>>1, st&1, ed>>1, ed&1, max_dist); + // } + assert(p); assert(b->S.a[0] == (ed^1)); + for (k = max_dist = 0; k < b->b.n; ++k) { + if(b->b.a[k]==st || b->b.a[k]==b->S.a[0]) continue; + max_dist += s->ref->seq[b->b.a[k]>>1].len; + } + max_dist += s->ref->seq[st>>1].len; + max_dist += s->ref->seq[b->S.a[0]>>1].len; + p = rd_hm_bub(s->g, s->ref, st, max_dist, b); + if(p) { + assert(b->S.a[0] == (ed^1)); + return; + } + + for (cuttoff = step; cuttoff < 1.0; cuttoff += step) { + ncut = rd_hm_drop(s->g, s->ref, st, ed^1, cuttoff, b); + p = rd_hm_bub(s->g, s->ref, st, max_dist, b); + if(p) { + assert(b->S.a[0] == (ed^1)); + return; + } + if(!ncut) break; + } + rd_hm_drop(s->g, s->ref, st, ed^1, 1024, b); + p = rd_hm_bub(s->g, s->ref, st, max_dist, b); + if(p) { + assert(b->S.a[0] == (ed^1)); + return; + } +} + +void reduce_hamming_error_adv(ma_ug_t *iug, asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, +int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub) +{ + double index_time = yak_realtime(); + ma_ug_t *ug = NULL; rd_hamming_t aux_t; memset((&aux_t), 0, sizeof(aux_t)); + ug = (iug)?(iug):(ma_ug_gen_primary(sg, PRIMARY_LABLE)); + uint8_t* vis_flag = NULL; CALLOC(vis_flag, sg->n_seq*2); + uint32_t fix_bub = 0; asg_t *g = ug->g; + uint32_t v, n_vtx = g->n_seq * 2, n_arc, n_arc_0 = sg->n_arc, nv, i; + uint64_t n_pop = 0, max_dist; asg_arc_t *p; + asg_arc_t *av; asg_t *ig = asg_init(); asg64_v sv; kv_init(sv); + buf_t b; memset(&b, 0, sizeof(buf_t)); + b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); + uint8_t* bs_flag = NULL; CALLOC(bs_flag, n_vtx); + for (i = 0; i < ug->g->n_seq; i++) ug->g->seq[i].c = 0; + max_dist = get_bub_pop_max_dist_advance(g, &b); + + if(max_dist > 0) { + if(bub) { + for (i = 0; i < bub->f_bub; i++) { + get_bubbles(bub, i, &v, NULL, NULL, NULL, NULL); + fix_bub += gen_switch_phasing(sg, ug, v, sources, coverage_cut, max_hang, min_ovlp, + R_INF.trio_flag, vis_flag, &b, max_dist, 1, ig, &sv); + } + } else { + for (v = 0; v < n_vtx; ++v) { + if(bs_flag[v] != 0) continue; + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + ///some node could be deleted + if (nv < 2 || g->seq[v>>1].del) continue; + ///some edges could be deleted + for (i = n_arc = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs + if (!av[i].del) ++n_arc; + if (n_arc < 2) continue; + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, max_dist, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) { + //beg is v, end is b.S.a[0] + //note b.b include end, does not include beg + for (i = 0; i < b.b.n; i++) { + if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue; + bs_flag[b.b.a[i]] = bs_flag[b.b.a[i]^1] = 1; + } + bs_flag[v] = 2; bs_flag[b.S.a[0]^1] = 3; + } + } + + //traverse all node with two directions + for (v = 0; v < n_vtx; ++v) { + if(bs_flag[v] !=2) continue; + nv = asg_arc_n(g, v); + av = asg_arc_a(g, v); + ///some node could be deleted + if (nv < 2 || g->seq[v>>1].del) continue; + ///some edges could be deleted + for (i = n_arc = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs + if (!av[i].del) ++n_arc; + if (n_arc > 1) { + fix_bub += gen_switch_phasing(sg, ug, v, sources, coverage_cut, max_hang, min_ovlp, + R_INF.trio_flag, vis_flag, &b, max_dist, 1, ig, &sv); + } + } + } + } + free(bs_flag); if(!iug) ma_ug_destroy(ug); + + if(sv.n > 0) { + ig->n_seq = ig->m_seq = sg->n_seq; + MALLOC(ig->seq, ig->n_seq); + memcpy(ig->seq, sg->seq, (sizeof((*(ig->seq)))*ig->n_seq)); + asg_cleanup(ig); asg_arc_del_trans_aux(ig, sg, vis_flag, gap_fuzz); + aux_t.n_thread = 1/**asm_opt.thread_num**/; CALLOC(aux_t.a, aux_t.n_thread); + REALLOC(b.a, (ig->n_seq<<1)); memset(b.a, 0, sizeof((*(b.a)))*(ig->n_seq<<1)); + for (i = 0; i < aux_t.n_thread; i++) aux_t.a[i].a = b.a; + aux_t.g = ig; aux_t.ref = sg; aux_t.rr = &sv; aux_t.max_dist = max_dist; + // print_debug_gfa(ug, sg, coverage_cut, "debug_hamming", sources, ru); + // kt_for(aux_t.n_thread, rd_hamming_symm, &aux_t, aux_t.rr->n); + for (i = 0; i < aux_t.rr->n; i++) { + rd_hamming_symm_simple(&aux_t, aux_t.rr->a[i]>>32, (uint32_t)(aux_t.rr->a[i])); + } + for (i = 0; i < aux_t.n_thread; i++) { + free(aux_t.a[i].S.a); free(aux_t.a[i].T.a); + free(aux_t.a[i].b.a); free(aux_t.a[i].e.a); + } + free(aux_t.a); + } + free(sv.a); free(vis_flag); + + for (i = n_pop = 0; i < ig->n_arc; i++) { + if(ig->arc[i].del) continue; + p = asg_arc_pushp(sg); *p = (ig->arc[i]); n_pop++; + } + if(n_pop) { + free(sg->idx); + sg->idx = 0; + sg->is_srt = 0; + asg_cleanup(sg); + asg_symm(sg); + asg_arc_del_trans(sg, gap_fuzz); + } + free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); asg_destroy(ig); + fprintf(stderr, "[M::%s::%.3f] # inserted edges: %u, # fixed bubbles: %u\n", + __func__, yak_realtime() - index_time, sg->n_arc - n_arc_0, fix_bub); +} void reduce_hamming_error(asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, int max_hang, int min_ovlp, long long gap_fuzz) @@ -30672,7 +31237,7 @@ int max_hang, int min_ovlp, uint32_t chainLenThres, long long gap_fuzz, bub_labe ma_ug_destroy(ug); ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); reset_bub(&bub, ug, cov->t_ch, &new_rtg_edges); // rescue_missing_hap_ovlp(ug, sg, sources, coverage_cut, max_hang, min_ovlp, &bub, gap_fuzz); - reduce_hamming_error_adv(ug, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, &bub); + reduce_hamming_error_adv(ug, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, ruIndex, &bub); } destory_bubbles(&bub); diff --git a/Trio.cpp b/Trio.cpp index f03db1b..878b5c2 100644 --- a/Trio.cpp +++ b/Trio.cpp @@ -349,6 +349,27 @@ static void ha_triobin_list(const hifiasm_opt_t *opt) fprintf(stderr, "[M::%s::%.3f*%.2f] ==> partitioned reads with external lists\n", __func__, yak_realtime(), yak_cpu_usage()); } +uint32_t test_yak_binning(char* fn, char *cmd) +{ + gzFile fp; kstream_t *ks; kstring_t str = {0,0,0}; + int dret, eq = 0; fp = gzopen(fn, "r"); + if (fp == 0) return eq; + ks = ks_init(fp); + if(ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) { + uint64_t sl = strlen(str.s), z; + for (z = 0; (z < sl) && (str.s[z]!='\t'); z++); + if(((z+1)r_g, opt, asm_opt.output_file_name); // reduce_hamming_error(h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz); - reduce_hamming_error_adv(NULL, h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz, NULL); + reduce_hamming_error_adv(NULL, h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz, opt->ruIndex, NULL); /** scaffold_hap(h, t_idx, opt, round, asm_opt.output_file_name, FATHER); scaffold_hap(h, t_idx, opt, round, asm_opt.output_file_name, MOTHER); diff --git a/htab.cpp b/htab.cpp index d63f098..e147153 100644 --- a/htab.cpp +++ b/htab.cpp @@ -1546,6 +1546,7 @@ int uidx_load(void **r_flt_tab, ha_pt_t **r_ha_idx, char* file_name, ma_ug_t *ug free(gfa_name); fclose(fp); return 0; } + // fprintf(stderr, "[M::%s]\t%s\tftell::%ld\n", __func__, file_name, ftell(fp)); ha_pt_t *ha_idx = NULL; char mode = 0; diff --git a/htab.h b/htab.h index b013ac4..afe44b7 100644 --- a/htab.h +++ b/htab.h @@ -114,6 +114,7 @@ double yak_peakrss_in_gb(void); double yak_cpu_usage(void); void ha_triobin(const hifiasm_opt_t *opt); +uint32_t test_yak_binning(char* fn, char *cmd); uint32_t *ha_polybin_list(const hifiasm_opt_t *opt); void mz1_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique, void *km); diff --git a/inter.cpp b/inter.cpp index e44d57a..a5c6b55 100644 --- a/inter.cpp +++ b/inter.cpp @@ -10449,7 +10449,7 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac } // 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; @@ -10580,7 +10580,7 @@ static void *worker_ul_rescall_pipeline(void *data, int step, void *in) // callb } // 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); + // 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) { @@ -15594,6 +15594,7 @@ void gen_UL_ovlps(uldat_t *sl, int32_t cutoff) ul_idx_t *uu = dedup_HiFis(sl->uopt, 1, 0); int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name, NULL) : 0); if(exist == 0) uidx_l_build(uu->ug, (mg_idxopt_t *)sl->opt, cutoff); + // print_debug_gfa(sl->uopt, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name, NULL); sl->ha_flt_tab = ha_flt_tab; sl->ha_idx = (ha_pt_t *)ha_idx; sl->uu = uu; ul_v_call(sl, asm_opt.ar); @@ -15699,6 +15700,7 @@ int32_t write_all_ul_t(all_ul_t *x, char* file_name, ma_ug_t *ug) fprintf(stderr, "[M::%s] Index has been written.\n", __func__); fclose(fp); + if(asm_opt.bin_only) exit(1); return 1; }