diff --git a/CommandLines.cpp b/CommandLines.cpp index 3dac9cf..37978f6 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -262,6 +262,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->is_topo_trans = 1; asm_opt->is_bub_trans = 1; asm_opt->bin_only = 0; + asm_opt->ul_clean_round = 2; } void destory_enzyme(enzyme* f) diff --git a/CommandLines.h b/CommandLines.h index 6f6eee5..5756596 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.18.5-r500" +#define HA_VERSION "0.18.6-r502" #define VERBOSE 0 @@ -137,6 +137,7 @@ typedef struct { uint8_t is_topo_trans; uint8_t is_bub_trans; uint8_t bin_only; + int32_t ul_clean_round; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Overlaps.cpp b/Overlaps.cpp index 30838a6..cb845e3 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -18754,10 +18754,11 @@ void filter_set_kug(uint8_t* trio_flag, asg_t *rg, uint8_t *rf, kvec_asg_arc_t_w void output_trio_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, -ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, -long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, -float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench, bub_label_t* b_mask_t) +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, int is_bench, long long gap_fuzz, ug_opt_t *opt, bub_label_t* b_mask_t) { + reduce_hamming_error_adv(NULL, sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz, opt->ruIndex, NULL); uint8_t *rf = NULL; if(asm_opt.kpt_rate > 0) CALLOC(rf, sg->n_seq); @@ -31231,7 +31232,7 @@ int max_hang, int min_ovlp, uint32_t chainLenThres, long long gap_fuzz, bub_labe reset_bub(&bub, ug, cov->t_ch, &new_rtg_edges); 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((!no_trio_recover) && (ha_opt_triobin(&asm_opt))) { ma_ug_destroy(ug); ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); @@ -31239,6 +31240,7 @@ int max_hang, int min_ovlp, uint32_t chainLenThres, long long gap_fuzz, bub_labe // 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, ruIndex, &bub); } + **/ destory_bubbles(&bub); destory_hap_cov_t(&cov); @@ -33585,7 +33587,7 @@ void renew_g(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, long lo 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) +int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, const char *bin_file, uint64_t free_uld) { ul_renew_t nopt; memset(&nopt, 0, sizeof(nopt)); nopt.src = sources; nopt.r_src = reverse_sources; nopt.ruIndex = ruIndex; @@ -33593,7 +33595,7 @@ int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file) nopt.cov = coverage_cut; nopt.b_mask_t = b_mask_t; nopt.max_hang = max_hang_length; nopt.mini_ovlp = mini_overlap_length; ul_realignment_gfa(uopt, *sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, - asm_opt.max_short_tip, asm_opt.max_short_ul_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file, &nopt);; + asm_opt.max_short_tip, asm_opt.max_short_ul_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file, &nopt, bin_file, free_uld); // 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); @@ -33607,6 +33609,40 @@ int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file) // post_rescue(uopt, *sg, (*sources), (*reverse_sources), ruIndex, b_mask_t, 0); } +void gradually_renew_g(ma_hit_t_alloc **src, ma_hit_t_alloc **rev_src, long long *n_read, +uint64_t **readLen, ma_sub_t **cov, 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, int64_t gap_fuzz, int64_t min_dp, +bub_label_t *b_mask_t, uint32_t is_trio, int32_t ul_aln_round, char *o_file, const char *bin_file) +{ + int32_t k, strl = strlen(bin_file)+1, kt, cl, sl; char *id = NULL; + renew_g(src, rev_src, n_read, readLen, cov, 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, + is_trio, o_file, bin_file, (ul_aln_round<=1)?1:0); + gen_ug_opt_t(uopt, *src, *rev_src, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, *readLen, + *cov, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, b_mask_t); + ug_ext_gfa(uopt, *sg, ug_ext_len); + /**if(!ha_opt_triobin(&asm_opt))**/ hic_clean_adv(*sg, uopt); + + cl = strl+1; MALLOC(id, cl); + for (k = 1; k < ul_aln_round; k++) { + for(kt = k, sl = strl+1; kt > 0; kt/=10) sl++; + if(cl < sl) { + cl = sl; REALLOC(id, cl); + } + + sprintf(id, "%s%d", bin_file, k); + renew_g(src, rev_src, n_read, readLen, cov, 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, + is_trio, o_file, id, ((k+1)==ul_aln_round)?1:0); + gen_ug_opt_t(uopt, *src, *rev_src, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, *readLen, + *cov, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, b_mask_t); + ug_ext_gfa(uopt, *sg, ug_ext_len); + /**if(!ha_opt_triobin(&asm_opt))**/ hic_clean_adv(*sg, uopt); + } + free(id); +} + 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, @@ -33711,6 +33747,20 @@ ma_sub_t **coverage_cut_ptr, int debug_g) // set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length); } + if(asm_opt.ar) { + gradually_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, gap_fuzz, min_dp, &b_mask_t, + ha_opt_triobin(&asm_opt), asm_opt.ul_clean_round, o_file, "re"); + } else { + ug_ext_gfa(&uopt, sg, ug_ext_len); + if(!ha_opt_triobin(&asm_opt)) { + // output_unitig_graph(sg, coverage_cut, "pre_clean", sources, ruIndex, max_hang_length, mini_overlap_length); + hic_clean_adv(sg, &uopt); + } + } + + /** if(asm_opt.ar) { 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, @@ -33718,11 +33768,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g) ///make sure uopt has been updated; not necessary 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); - // 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); /** asg_cut_tip(sg, asm_opt.max_short_tip); @@ -33853,12 +33900,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g) // flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL; flat_soma_v(sg, sources, ruIndex); **/ - ug_ext_gfa(&uopt, sg, ug_ext_len); - - if(!ha_opt_triobin(&asm_opt)) { - // output_unitig_graph(sg, coverage_cut, "pre_clean", sources, ruIndex, max_hang_length, mini_overlap_length); - hic_clean_adv(sg, &uopt); - } output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources, asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length); @@ -33885,7 +33926,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) { if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION; output_trio_graph(sg, coverage_cut, o_file, sources, reverse_sources, (asm_opt.max_short_tip*2), - 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 0, &b_mask_t); + 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 0, gap_fuzz, &uopt, &b_mask_t); } else if(ha_opt_hic(&asm_opt)) { diff --git a/Overlaps.h b/Overlaps.h index 8d49dad..22b66c1 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -32,7 +32,6 @@ // #define PRIMARY_LABLE 1 // #define ALTER_LABLE 2 // #define HAP_LABLE 4 -#define HA_RE_UL_ID "re" #define ug_ext_len 75000 #define Get_qn(RECORD) ((uint32_t)((RECORD).qns>>32)) diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 720f8ed..641f3ee 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -2270,6 +2270,20 @@ void print_vw_edge(asg_t *sg, uint32_t v, uint32_t w, const char *cmd) if(i >= nv) fprintf(stderr, "[%s]\tno edges\n", cmd); } +void prt_spec_edge(asg_t *rg, ma_hit_t_alloc *src, uint32_t tot_rid, uint32_t vid, uint32_t wid, const char *cmd) +{ + if(vid >= tot_rid) return; + uint32_t k, qn, tn; ma_hit_t_alloc *s = &(src[vid]); + for (k = 0; k < s->length; k++) { + qn = Get_qn(s->buffer[k]); tn = Get_tn(s->buffer[k]); + if(qn == vid && tn == wid) { + fprintf(stderr, "[M::%s]%s\tqn::%u\tq::[%u,\t%u)\t%c\ttn::%u\tt::[%u,\t%u)\n", + __func__, cmd, qn, Get_qs(s->buffer[k]), Get_qe(s->buffer[k]), + "+-"[s->buffer[k].rev], tn, Get_ts(s->buffer[k]), Get_te(s->buffer[k])); + } + } +} + int32_t gen_spec_edge(asg_t *rg, ug_opt_t *uopt, uint32_t v, uint32_t w, asg_arc_t *t) { uint32_t k, qn, tn; int32_t r; ma_hit_t_alloc *s = &(uopt->sources[v>>1]); asg_arc_t p; @@ -13717,10 +13731,10 @@ typedef struct { #define B4Lg(x) (((x)>>2)+(((x)&3)?1:0)) -uint32_t load_scaf_base(all_ul_t *x, char* file_name) +uint32_t load_scaf_base(all_ul_t *x, char* file_name, const char *bin_file) { char *gfa_name = (char*)malloc(strlen(file_name)+50); - sprintf(gfa_name, "%s.%s.uidx.ucr.bin", file_name, HA_RE_UL_ID); + sprintf(gfa_name, "%s.%s.uidx.ucr.bin", file_name, bin_file); // fprintf(stderr, "[M::%s] open %s...\n", __func__, gfa_name); FILE *fp = fopen(gfa_name, "r"); free(gfa_name); if (!fp) return 0; @@ -13986,7 +14000,7 @@ uint64_t rid_n, uint64_t scaf_len, char *scaf_id, asg_t *ng, ug_opt_t *uopt, ma_ 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); + sprintf(gfa_name, "%s.%s.ul.ovlp.bin", file_name, "re"); FILE* fp = fopen(gfa_name, "r"); free(gfa_name); if (!fp) return 0; uint64_t k; size_t kn; ul_vec_t *p = NULL; @@ -14128,10 +14142,10 @@ static void gen_scaffold_id(void *data, long i, int tid) // callback for kt_for( } } -void reload_uu(all_ul_t *x, ma_ug_t *raw_g, char* file_name) +void reload_uu(all_ul_t *x, ma_ug_t *raw_g, char* file_name, const char *bin_file) { char* gfa_name = NULL; MALLOC(gfa_name, strlen(asm_opt.output_file_name)+50); - sprintf(gfa_name, "%s.%s", asm_opt.output_file_name, HA_RE_UL_ID); + sprintf(gfa_name, "%s.%s", asm_opt.output_file_name, bin_file); clear_all_ul_t(x); load_all_ul_t(x, gfa_name, &R_INF, raw_g); // filter_ul_ug(raw_g);//no filter since we would like to use all alignments @@ -14139,10 +14153,10 @@ void reload_uu(all_ul_t *x, ma_ug_t *raw_g, char* file_name) free(gfa_name); } -void fill_scaffolds(usc_t *a, uint64_t a_n, ma_ug_t *raw_g, int64_t min_arc_len) +void fill_scaffolds(usc_t *a, uint64_t a_n, ma_ug_t *raw_g, int64_t min_arc_len, const char *bin_file) { scaf_mul_t ss; uint32_t k, occ = 0; - reload_uu(&UL_INF, raw_g, asm_opt.output_file_name); + reload_uu(&UL_INF, raw_g, asm_opt.output_file_name, bin_file); // 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; @@ -14164,7 +14178,7 @@ void fill_scaffolds(usc_t *a, uint64_t a_n, ma_ug_t *raw_g, int64_t min_arc_len) } // fprintf(stderr, "[M::%s] occ::%u\n", __func__, occ); - if(occ) load_scaf_base(&UL_INF, asm_opt.output_file_name); + if(occ) load_scaf_base(&UL_INF, asm_opt.output_file_name, bin_file); } 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); @@ -14239,15 +14253,22 @@ void push_direct_scaff_node(usc_t *psa, asg_t *ng, ug_opt_t *uopt, uint64_t vx, if(dif > 16) r = -1; } } - // if((((vx>>1) == 3313901) && ((wx>>1) == 3555473)) || (((vx>>1) == 3555473) && ((wx>>1) == 3313901))) { - // fprintf(stderr, "[M::%s::]\t%.*s(%c)(sv>>1::%lu)\t%.*s(%c)(sw>>1::%lu)\tr::%d\n", __func__, - // (int)Get_NAME_LENGTH(R_INF, (vx>>1)), Get_NAME(R_INF, (vx>>1)), "+-"[vx&1], vx>>1, - // (int)Get_NAME_LENGTH(R_INF, (wx>>1)), Get_NAME(R_INF, (wx>>1)), "+-"[wx&1], wx>>1, r); - // } if(r >= 0) { - p = asg_arc_pushp(ng); *p = t; r = gen_spec_rc_edge(ng, uopt, wx^1, vx^1, &t, &rc1); - assert(r >= 0); p = asg_arc_pushp(ng); *p = t; + if(uol > t.ol) dif = uol - t.ol; + else dif = t.ol - uol; + if((dif > (uol*0.05)) || (dif > (t.ol*0.05))) { + if(dif > 16) r = -1; + } + } + + if(r >= 0) { + r = gen_spec_rc_edge(ng, uopt, vx, wx, &t, &rc0); assert(r >= 0); + p = asg_arc_pushp(ng); *p = t; + + r = gen_spec_rc_edge(ng, uopt, wx^1, vx^1, &t, &rc1); assert(r >= 0); + p = asg_arc_pushp(ng); *p = t; + add_ma_hit_t_alloc(&(R_INF.paf[Get_qn((*rc0))]), rc0); add_ma_hit_t_alloc(&(R_INF.paf[Get_qn((*rc1))]), rc1); } else { @@ -14307,7 +14328,7 @@ void print_debug_scaffold_nodes(ug_opt_t *uopt, usc_t *a, uint32_t a_n, asg_t *n destory_UC_Read(&r0); destory_UC_Read(&r1); } -asg_t *renew_ng(usg_t *eg, ma_ug_t *rug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **cov, R_to_U *ruI, uint64_t scaffold_len) +asg_t *renew_ng(usg_t *eg, ma_ug_t *rug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **cov, R_to_U *ruI, uint64_t scaffold_len, const char *bin_file) { init_aux_table(); ma_utg_t *u; uint64_t i, v, w, m, h, z, raw_v, raw_w, nocc, nv, vx, wx; int32_t r; @@ -14323,7 +14344,6 @@ asg_t *renew_ng(usg_t *eg, ma_ug_t *rug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **c memset(ext.idx_a.a, -1, sizeof(*(ext.idx_a.a))*ext.idx_a.n);//map ug1 = convert_usg_t(eg, rug); - for (i = 0; i < ug1->u.n; ++i) { if(ug1->g->seq[i].del) continue; u = &(ug1->u.a[i]); @@ -14355,7 +14375,6 @@ asg_t *renew_ng(usg_t *eg, ma_ug_t *rug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **c 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 < eg->n; ++i) { if(ug1->g->seq[i].del) continue; ///there shouldn't any scaffolding within the nodes @@ -14415,13 +14434,12 @@ asg_t *renew_ng(usg_t *eg, ma_ug_t *rug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **c 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, rug, uopt->min_ovlp); - realloc_rdb_adv(&(R_INF), &UL_INF, cov, ruI, ext.idx_a.a, ext.idx_a.n, scaffold_len, (char *)"scaf", ng, uopt, rug, sa.a); + if(sa.n > 0) fill_scaffolds(sa.a, sa.n, rug, uopt->min_ovlp, bin_file); + + realloc_rdb_adv(&(R_INF), &UL_INF, cov, ruI, ext.idx_a.a, ext.idx_a.n, scaffold_len, (char *)"scaf", ng, uopt, rug, sa.a); 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 < eg->n; ++i) { if(ug1->g->seq[i].del) continue; u = &(ug1->u.a[i]); @@ -14449,6 +14467,13 @@ asg_t *renew_ng(usg_t *eg, ma_ug_t *rug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **c if(r >= 0) { p = asg_arc_pushp(ng); *p = t; r = gen_spec_edge(ng, uopt, wx^1, vx^1, &t); + // if(!(r >= 0)) { + // fprintf(stderr, "[M::%s]\tvx>>1::%lu(%c)\twx>>1::%lu(%c)\n", + // __func__, vx>>1, "+-"[vx&1], wx>>1, "+-"[wx&1]); + // fprintf(stderr, "[M::%s]\t%.*s(%c)->%.*s(%c)\n", + // __func__, (int)Get_NAME_LENGTH(R_INF, (vx>>1)), Get_NAME(R_INF, (vx>>1)), "+-"[vx&1], + // (int)Get_NAME_LENGTH(R_INF, (wx>>1)), Get_NAME(R_INF, (wx>>1)), "+-"[wx&1]); + // } assert(r >= 0); p = asg_arc_pushp(ng); *p = t; } else { psa = &(sa.a[nocc-ng->r_seq]); z = nocc<<1; @@ -14514,7 +14539,6 @@ asg_t *renew_ng(usg_t *eg, ma_ug_t *rug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **c } } } - // print_debug_scaffold_nodes(uopt, sa.a, sa.n, ng); asg_cleanup(ng); ng->r_seq = ng->n_seq; free(ext.cnt.a); free(ext.idx_a.a); free(sa.a); ma_ug_destroy(ug1); @@ -16010,7 +16034,7 @@ ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt, uin asg_arc_del_trans(z->i_g, uopt->gap_fuzz); z->i_ug = ma_ug_gen(z->i_g); - CALLOC(z->i_ug->g->seq_vis, z->i_ug->g->n_seq*2); + // CALLOC(z->i_ug->g->seq_vis, z->i_ug->g->n_seq*2); renew_u2g_cov(uidx); renew_u2g_bg(uidx); @@ -16205,11 +16229,66 @@ 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) +inline uint64_t dup_paf_check(ma_hit_t_alloc *x, ma_hit_t *p) +{ + int64_t i; ma_hit_t *z; + for (i = 0; i < x->length; i++) { + z = &(x->buffer[i]); + if((z->qns == p->qns) && (z->tn == p->tn) && (z->qe == p->qe) && (z->ts == p->ts) && (z->te == p->te) && + (z->cc == p->cc) && (z->ml == p->ml) && (z->rev == p->rev) && (z->bl == p->bl) && (z->del == p->del) && + (z->el == p->el) && (z->no_l_indel == p->no_l_indel)) { + return 0; + } + } + return 1; +} + +void renew_paf0(ma_hit_t_alloc *paf, uint64_t *a0, uint64_t a0n, uint64_t *a1, uint64_t a1n, asg64_v *buf) { 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; + uint64_t qn = ((uint32_t)a0[0]), tn = ((uint32_t)a1[0]), i, k, sf = 0; + ma_hit_t e01, e10; uint64_t qi, ti, *qa, *ta, qlen, tlen; buf->n = 0; ma_hit_t_alloc *qo, *to; + if(qn == tn) sf = 1; + qo = &(paf[qn]); qi = qo->length; + for (i = qlen = 0, qa = NULL; i < qi; i++) { + if(qo->buffer[i].tn == tn) kv_push(uint64_t, *buf, i); + } + qlen = buf->n; + + to = &(paf[tn]); ti = to->length; + for (i = tlen = 0, ta = NULL; i < ti; i++) { + if(to->buffer[i].tn == qn) kv_push(uint64_t, *buf, i); + } + tlen = buf->n - qlen; + assert(qlen && tlen); + qa = buf->a; ta = buf->a + qlen; + + for (qi = 0; qi < qlen; qi++) { + e01 = qo->buffer[qa[qi]]; + for (ti = 0; ti < tlen; ti++) { + e10 = to->buffer[ta[ti]]; + + 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); + if((!sf) || (dup_paf_check(&(paf[qn]), &e01))) { + add_ma_hit_t_alloc(&(paf[qn]), &e01); + } + + update_qtn(&e10, tn, qn); + if((!sf) || (dup_paf_check(&(paf[tn]), &e10))) { + add_ma_hit_t_alloc(&(paf[tn]), &e10); + } + } + } + } + } + + /** idx = get_specific_overlap(&(paf[qn]), qn, tn); if(idx < 0) return; e01 = paf[qn].buffer[idx]; @@ -16221,10 +16300,15 @@ void renew_paf0(ma_hit_t_alloc *paf, uint64_t *a0, uint64_t a0n, uint64_t *a1, u for (k = 0; k < a1n; k++) { if(i == 0 && k == 0) continue; tn = ((uint32_t)a1[k]); + if((qn == 5619628 && tn == 5619629) || (tn == 5619628 && qn == 5619629)) { + fprintf(stderr, "[M::%s]\tqn::%lu(qg::%u)\ttn::%lu(tg::%u)\n", __func__, + qn, ((uint32_t)a0[0]), tn, ((uint32_t)a1[0])); + } 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) @@ -16246,10 +16330,10 @@ void renew_paf1(ma_hit_t_alloc *paf, uint64_t *a, uint64_t an, uint64_t len) 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; + asg64_v clus, bb; 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); + CALLOC(idx, pre_gn); kv_init(bb); 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 @@ -16268,7 +16352,12 @@ void update_paf(ma_hit_t_alloc *src, ma_hit_t_alloc *r_src, uint64_t *rmap, uint 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); - + // if(k == 4265018) { + // fprintf(stderr, "\n[M::%s]\tgroup::%lu\n", __func__, k); + // for (j = 0; j < a0n; j++) { + // fprintf(stderr, "[M::%s]\tnid::%uu\n", __func__, (uint32_t)a0[j]); + // } + // } z = &(src[k]); for (j = 0; j < z->length; j++) { qn = Get_qn(z->buffer[j]); @@ -16279,7 +16368,7 @@ void update_paf(ma_hit_t_alloc *src, ma_hit_t_alloc *r_src, uint64_t *rmap, uint 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); + renew_paf0(src, a0, a0n, a1, a1n, &bb); } z = &(r_src[k]); @@ -16292,7 +16381,7 @@ void update_paf(ma_hit_t_alloc *src, ma_hit_t_alloc *r_src, uint64_t *rmap, uint 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); + renew_paf0(r_src, a0, a0n, a1, a1n, &bb); } @@ -16306,14 +16395,13 @@ void update_paf(ma_hit_t_alloc *src, ma_hit_t_alloc *r_src, uint64_t *rmap, uint // } if(a0n > 1) renew_paf1(r_src, a0, a0n, ng->seq[k].len); } - free(idx); free(clus.a); + free(idx); free(clus.a); free(bb.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; @@ -16348,7 +16436,8 @@ uint64_t scaffold_len) 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 - + fprintf(stderr, "\n+[M::%s] 0\n", __func__); + for (i = 0; i < ug->g->n_seq; ++i) { u = &(ug->u.a[i]); for (m = 0; m < u->n; m++) { @@ -16380,6 +16469,7 @@ uint64_t scaffold_len) // 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); + fprintf(stderr, "\n+[M::%s] 1\n", __func__); for (i = 0, nocc = ng->r_seq; i < ug->g->n_seq; ++i) { u = &(ug->u.a[i]); @@ -16437,13 +16527,16 @@ uint64_t scaffold_len) } } } + fprintf(stderr, "\n+[M::%s] 2\n", __func__); // 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); + fprintf(stderr, "\n+[M::%s] 3\n", __func__); + update_paf(R_INF.paf, R_INF.reverse_paf, ext.idx_a.a, ext.idx_a.n, sg->n_seq, ng); + fprintf(stderr, "\n+[M::%s] 4\n", __func__); for (i = 0, nocc = ng->r_seq; i < ug->g->n_seq; ++i) { u = &(ug->u.a[i]); @@ -16534,6 +16627,7 @@ uint64_t scaffold_len) } } } + fprintf(stderr, "\n+[M::%s] 5\n", __func__); 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", @@ -16669,9 +16763,9 @@ R_to_U* ruIndex, int max_hang, int min_ovlp) void 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, int64_t max_ul_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, -ul_renew_t *ropt) +ul_renew_t *ropt, const char *bin_file, uint64_t free_uld) { - uint64_t i; uint8_t *r_het = NULL; bubble_type *bub = NULL; ulg_opt_t uu; + uint64_t i, bn = 0, idn = 0; uint32_t *bl = NULL; uint8_t *r_het = NULL; bubble_type *bub = NULL; ulg_opt_t uu; for (i = 0; i < sg->n_seq; ++i) { if(sg->seq[i].del) continue; sg->seq[i].c = PRIMARY_LABLE; @@ -16679,17 +16773,21 @@ ul_renew_t *ropt) // fprintf(stderr, "0[M::%s]\n", __func__); hic_clean(sg); // fprintf(stderr, "1[M::%s]\n", __func__); - ma_ug_t *init_ug = ul_realignment(uopt, sg, 0); + ma_ug_t *init_ug = ul_realignment(uopt, sg, 0, bin_file); // fprintf(stderr, "2[M::%s]\n", __func__); // exit(1); filter_sg_by_ug(sg, init_ug, uopt); - // fprintf(stderr, "3[M::%s]\n", __func__); + // fprintf(stderr, "-0-[M::%s]\tUL_INF.a[25].rlen::%u\n", __func__, UL_INF.a[25].rlen); // 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_ul_alignment(init_ug, &UL_INF, 47072, "after-0"); bub = gen_bubble_chain(sg, init_ug, uopt, &r_het); // fprintf(stderr, "4[M::%s]\n", __func__); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-1"); ul_resolve_t *uidx = init_ul_resolve_t(sg, init_ug, bub, &UL_INF, uopt, r_het); + if(!free_uld) {///backup + bn = UL_INF.n; idn = UL_INF.nid.n; MALLOC(bl, bn); + for (i = 0; i < bn; i++) bl[i] = UL_INF.a[i].rlen; + } // fprintf(stderr, "5[M::%s]\n", __func__); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-2"); // exit(1); @@ -16709,16 +16807,30 @@ ul_renew_t *ropt) // uidx->uovl.hybrid_ug = gen_hybrid_ug(uidx, uidx->uovl.h_usg); // 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, init_ug, uopt->coverage_cut, bin_file, 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); // return uidx->uovl.hybrid_ug; - asg_t *ng = renew_ng(uidx->uovl.h_usg, init_ug, sg, uopt, ropt->cov, ropt->ruIndex, 16); - + asg_t *ng = renew_ng(uidx->uovl.h_usg, init_ug, sg, uopt, ropt->cov, ropt->ruIndex, 16, bin_file); destory_bubbles(uidx->bub); free(uidx->bub); ma_ug_destroy(init_ug); free(r_het); - destroy_ul_resolve_t(uidx); destory_all_ul_t(&UL_INF); + destroy_ul_resolve_t(uidx); + if(free_uld) { + destory_all_ul_t(&UL_INF); memset((&UL_INF), 0, sizeof(UL_INF)); + } else { + for (i = 0; i < bn; i++) UL_INF.a[i].rlen = bl[i]; + for (i = bn; i < UL_INF.n; i++) { + free(UL_INF.a[i].N_site.a); free(UL_INF.a[i].r_base.a); + free(UL_INF.a[i].bb.a); memset(&(UL_INF.a[i]), 0, sizeof(UL_INF.a[i])); + } + UL_INF.n = bn; + + for (i = idn; i < UL_INF.nid.n; i++) { + free(UL_INF.nid.a[i].a); memset(&(UL_INF.nid.a[i]), 0, sizeof(UL_INF.nid.a[i])); + } + UL_INF.nid.n = idn; free(bl); + } asg_destroy((*(ropt->sg))); (*(ropt->sg)) = ng; (*(ropt->src)) = R_INF.paf; (*(ropt->r_src)) = R_INF.reverse_paf; diff --git a/gfa_ut.h b/gfa_ut.h index c309526..87d513c 100644 --- a/gfa_ut.h +++ b/gfa_ut.h @@ -27,7 +27,7 @@ void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float o 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, -double max_ovlp_drop_ratio, int64_t max_tip, int64_t max_ul_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, ul_renew_t *ropt); +double max_ovlp_drop_ratio, int64_t max_tip, int64_t max_ul_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, ul_renew_t *ropt, const char *bin_file, uint64_t free_uld); 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); diff --git a/inter.cpp b/inter.cpp index 5d11de8..d49b41d 100644 --- a/inter.cpp +++ b/inter.cpp @@ -16125,13 +16125,13 @@ void trans_base_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res, } -ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_cache) +ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_cache, const char *bin_file) { fprintf(stderr, "[M::%s::] ==> starting UL\n", __func__); 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.%s", asm_opt.output_file_name, HA_RE_UL_ID); + sprintf(gfa_name, "%s.%s", asm_opt.output_file_name, bin_file); init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); cutoff = REA_ALIGN_CUTOFF; @@ -16147,9 +16147,13 @@ ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_c gen_UL_reovlps(&sl, ug, sg, gfa_name, cutoff); // exit(1); write_all_ul_t(&UL_INF, gfa_name, ug); - } else if(double_check_cache){ - if(drenew_UL_reovlps(&sl, ug, sg, gfa_name, cutoff)) { - write_all_ul_t(&UL_INF, gfa_name, ug); + } else{ + free(UL_INF.ridx.idx.a); free(UL_INF.ridx.occ.a); + memset(&(UL_INF.ridx), 0, sizeof((UL_INF.ridx))); + if(double_check_cache){ + if(drenew_UL_reovlps(&sl, ug, sg, gfa_name, cutoff)) { + write_all_ul_t(&UL_INF, gfa_name, ug); + } } } diff --git a/inter.h b/inter.h index 0a6f961..c25729f 100644 --- a/inter.h +++ b/inter.h @@ -108,7 +108,7 @@ void ul_resolve(ma_ug_t *ug, const asg_t *rg, const ug_opt_t *uopt, int hap_n); void ul_load(const ug_opt_t *uopt); uint64_t* get_hifi2ul_list(all_ul_t *x, uint64_t hid, uint64_t* a_n); uint64_t ul_refine_alignment(const ug_opt_t *uopt, asg_t *sg); -ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_cache); +ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_cache, const char *bin_file); int32_t write_all_ul_t(all_ul_t *x, char* file_name, ma_ug_t *ug); int32_t load_all_ul_t(all_ul_t *x, char* file_name, All_reads *hR, ma_ug_t *ug); uint32_t ugl_cover_check(uint64_t is, uint64_t ie, ma_utg_t *u);