From a3b72457f0c643747a670d2a6157a4aaee93b4e8 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Tue, 6 Dec 2022 09:52:00 -0500 Subject: [PATCH] cleaner assembly graph --- CommandLines.cpp | 3 + CommandLines.h | 5 +- Overlaps.cpp | 25 +++++- Purge_Dups.cpp | 2 +- anchor.cpp | 13 +++ gfa_ut.cpp | 223 ++++++++++++++++++++++++++++++++++++++++------- hic.cpp | 3 +- inter.cpp | 190 +++++++++++++++++++++++++++++++++++++++- 8 files changed, 426 insertions(+), 38 deletions(-) diff --git a/CommandLines.cpp b/CommandLines.cpp index dd7f771..c9e6eda 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -172,8 +172,10 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->k_mer_length = 51; asm_opt->hic_mer_length = 31; asm_opt->ul_mer_length = 19; + asm_opt->trans_mer_length = 31; asm_opt->mz_win = 51; asm_opt->ul_mz_win = 19; + asm_opt->trans_win = 31; asm_opt->mz_rewin = 1000; asm_opt->ul_mz_rewin = 360; asm_opt->mz_sample_dist = 500; @@ -208,6 +210,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->purge_level_trio = 0; asm_opt->purge_simi_rate_l2 = 0.75; asm_opt->purge_simi_rate_l3 = 0.55; + asm_opt->trans_base_rate = 0.8; asm_opt->purge_overlap_len = 1; ///asm_opt->purge_overlap_len_hic = 50; asm_opt->recover_atg_cov_min = -1024; diff --git a/CommandLines.h b/CommandLines.h index 3760230..eb062e8 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.18.0-r465" +#define HA_VERSION "0.18.1-r466" #define VERBOSE 0 @@ -48,9 +48,11 @@ typedef struct { int k_mer_length; int hic_mer_length; int ul_mer_length; + int trans_mer_length; int bub_mer_length; int mz_win; int ul_mz_win; + int trans_win; int mz_rewin; int ul_mz_rewin; int mz_sample_dist; @@ -97,6 +99,7 @@ typedef struct { float purge_simi_rate_l2; float purge_simi_rate_l3; float purge_simi_thres; + float trans_base_rate; ///float purge_simi_rate_hic; diff --git a/Overlaps.cpp b/Overlaps.cpp index e231d6d..f86158a 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -13823,6 +13823,8 @@ 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); @@ -13830,7 +13832,11 @@ long long gap_fuzz, bub_label_t* b_mask_t) // ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file); // fclose(output_file); // free(gfa_name); + // exit(1); ///for debug + + + hic_analysis(ug, sg, cov?cov->t_ch:t_ch, &opt, 0, asm_opt.scffold?&rhits:NULL); @@ -14876,12 +14882,29 @@ void filter_u_trans_t(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, uint32_t thr kt_u_trans_t_idx(ta, ug->g->n_seq); } + +void dbg_prt_utg_trans(kv_u_trans_t *ta, ma_ug_t *ug, const char *o_n) +{ + char* gfa_name = (char*)malloc(strlen(o_n)+100); + FILE *fn = NULL; sprintf(gfa_name, "%s.tran.dbg.ovlp.log", o_n); + u_trans_t *p = NULL; uint32_t i; fn = fopen(gfa_name, "w"); + for (i = 0; i < ta->n; i++) { + p = &(ta->a[i]); + fprintf(fn, "utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\tw(%f)\tf(%u)\n", + p->qn+1, "lc"[ug->u.a[p->qn].circ], ug->u.a[p->qn].len, p->qs, p->qe, "+-"[p->rev], + p->tn+1, "lc"[ug->u.a[p->tn].circ], ug->u.a[p->tn].len, p->ts, p->te, p->nw, p->f); + } + fclose(fn); free(gfa_name); fprintf(stderr, "[M::%s::] done\n", __func__); +} + void clean_u_trans_t_idx(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g) { + // dbg_prt_utg_trans(ta, ug, "pre"); kt_u_trans_t_idx(ta, ug->g->n_seq); kt_u_trans_t_symm(ta, ug); filter_u_trans_t(ta, ug, read_g, 3); ///debug_u_trans_t(ta); + // dbg_prt_utg_trans(ta, ug, "after"); } @@ -32112,7 +32135,7 @@ 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); **/ - if(!(ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt))) { + 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); } diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index 42c91af..2f81441 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -1788,7 +1788,7 @@ void chain_trans_ovlp(hap_cov_t *cover, utg_trans_t *o, ma_ug_t *ug, asg_t *read } p_v = v; len += l; - if(len >= targetBaseLen) + if(len >= targetBaseLen)///might be not right for centromeres { xOcc = aOcc; break; diff --git a/anchor.cpp b/anchor.cpp index 3cafdd1..e497b82 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -1278,3 +1278,16 @@ void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t ///no need to sort here, overlap_list has been sorted at lchain_gen lchain_qgen_mcopy(cl, overlap_list, rid, rl, NULL, uref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off); } + +void ug_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, + int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, uint32_t is_hpc) +{ + int64_t max_skip, max_iter, max_dis, quick_check; double chn_pen_gap, chn_pen_skip; + set_lchain_dp_op(is_accurate, mz_k, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check); + if((!overlap_list) || (!cl)) { + ab->mz.n = 0; + mz2_ha_sketch(rs, rl, mz_w, mz_k, 0, is_hpc, &ab->mz, NULL, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, NULL); + } else { + + } +} diff --git a/gfa_ut.cpp b/gfa_ut.cpp index c475c77..9f997d8 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -11510,18 +11510,20 @@ void debug_sysm_usg_t(usg_t *ng, const char *cmd) } } -void update_usg_t_threading_0(usg_t *ng, uint64_t *a, uint64_t a_n, uint32_t *occ, asg64_v *b) +void update_usg_t_threading_0(ul_resolve_t *uidx, usg_t *ng, uint64_t *a, uint64_t a_n, uint32_t *occ, asg64_v *b) { uint64_t k, i, *p, nid, nnid, v, w; b->n = 0; usg_seq_t *s; - 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, "utg%.6dl(%c)(hom::%u)\t", (int32_t)(a[k]>>1)+1, "+-"[a[k]&1], (IF_HOM((a[k]>>1), (*(uidx->bub))))); // } // fprintf(stderr, "\n"); + 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]]--; + 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) {///copy a node @@ -11627,16 +11629,29 @@ void update_usg_t_threading(ul_resolve_t *uidx, usg_t *ng, uint64_t *arcs, uint6 { uint64_t i, k, s, e, nvtx = ng->n<<1; memset(occ, 0, sizeof((*occ))*nvtx); for (i = 0; i < arcs_gn; i++) {///set cluster + // prt = 0; s = arcs[((uint32_t)arcs_g[i])]>>32; e = ((uint32_t)arcs[((uint32_t)arcs_g[i])]); assert(s < e); for (k = s + 1; k < e; k++) {///note: here is [s, e] occ[integ_seq[k]]++; occ[integ_seq[k]^1]++; + // if((integ_seq[k]>>1) == 2736) prt = 1; } occ[integ_seq[s]]++; occ[integ_seq[e]^1]++; + // if((integ_seq[s]>>1) == 2736) prt = 1; + // if(((integ_seq[e]^1)>>1) == 2736) prt = 1; // if(occ[integ_seq[s]] > 1 || occ[integ_seq[e]^1] > 1) { // fprintf(stderr, "s::utg%.6dl(%c)(occ::%u), e::utg%.6dl(%c)(occ::%u)\n", // (int32_t)(integ_seq[s]>>1)+1, "+-"[integ_seq[s]&1], occ[integ_seq[s]], // (int32_t)((integ_seq[e]^1)>>1)+1, "+-"[((integ_seq[e]^1)&1)], occ[integ_seq[e]^1]); // } + // if(prt) { + // fprintf(stderr, "[M::%s::] a_n::%lu\n", __func__, e + 1 - s); + // for (k = s; k <= e; k++) { + // fprintf(stderr, "utg%.6dl(%c)(hom::%u)\t", + // (int32_t)(integ_seq[k]>>1)+1, "+-"[integ_seq[k]&1], + // (IF_HOM((integ_seq[k]>>1), (*(uidx->bub))))); + // } + // fprintf(stderr, "\n"); + // } } for (i = 0; i < nvtx; i++) { @@ -11644,7 +11659,7 @@ void update_usg_t_threading(ul_resolve_t *uidx, usg_t *ng, uint64_t *arcs, uint6 } for (i = 0; i < arcs_gn; i++) {///arcs_g[]>>32 is the group id; (uint32_t)arcs_g[] s = arcs[((uint32_t)arcs_g[i])]>>32; e = ((uint32_t)arcs[((uint32_t)arcs_g[i])]); - update_usg_t_threading_0(ng, integ_seq + s, e + 1 - s, occ, b); + update_usg_t_threading_0(uidx, ng, integ_seq + s, e + 1 - s, occ, b); } } @@ -12245,39 +12260,68 @@ uint64_t get_arc_support(ul_resolve_t *uidx, uint64_t v, uint64_t w) // } // } +uint64_t get_arc_support_chain(ul_resolve_t *uidx, uint64_t *a, uint64_t a_n, usg_t *ng) +{ + if(a_n <= 0) return 0; + uint64_t *hid_a, hid_n; ul_str_idx_t *str_idx = &(uidx->pstr); + uint64_t z, vz, occ = 0, v, w, ai; ul_str_t *str; int64_t s, s_n; + v = (ng->a[a[0]>>1].mm<<1)|(a[0]&1); + hid_a = str_idx->occ.a + str_idx->idx.a[v>>1]; + hid_n = str_idx->idx.a[(v>>1)+1] - str_idx->idx.a[v>>1]; + + for (z = occ = 0; z < hid_n; z++) { + str = &(str_idx->str.a[hid_a[z]>>32]); s_n = str->cn; + if(s_n < 2) continue; + vz = (uint32_t)(str->a[(uint32_t)hid_a[z]]); + assert((v>>1) == (vz>>1)); + + if(v == vz) { + for(s = ((uint32_t)hid_a[z]) + 1, ai = 1; s < s_n && ai < a_n; s++, ai++) { + w = (ng->a[a[ai]>>1].mm<<1)|(a[ai]&1); + if(((uint32_t)(str->a[s])) != w) break; + } + if(ai >= a_n) occ++; + } else { + for(s = ((int32_t)((uint32_t)hid_a[z]))-1, ai = 1; s >= 0 && ai < a_n; s--, ai++) { + w = (ng->a[a[ai]>>1].mm<<1)|(a[ai]&1); + if(((uint32_t)(str->a[s])) != (w^1)) break; + } + if(ai >= a_n) occ++; + } + } + return occ; +} + static void worker_unique_bridge_check_s(void *data, long i, int tid) // callback for kt_for() { 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 *integer_seq = uaux->integer_seq, s, e, vk, wk, nse, k; uint64_t *gidx = uaux->gidx, *interval = uaux->interval_idx; 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++) { + for (k = s, vk = wk = (uint32_t)-1; k <= e; k++) { if(k > s) { - 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)); + // 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)); + nse = get_arc_support_chain(uaux->uidx, integer_seq+k-1, 2, ng); 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], - // ((int32_t)(integer_seq[k]>>1))+1, "+-"[integer_seq[k]&1], nse); gidx[i] |= ((uint64_t)0x8000000000000000); return; } } if(IF_HOM((integer_seq[k]>>1), *bub)) continue; - 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); + // w = (ng->a[integer_seq[k]>>1].mm<<1)|(integer_seq[k]&1); + wk = k; + if(vk != (uint32_t)-1) { + // nse = get_arc_support(uaux->uidx, v, w); + nse = get_arc_support_chain(uaux->uidx, integer_seq+vk, wk+1-vk, ng); if(nse < unique_bridge_occ) { - // fprintf(stderr, "-utg%.6dl(%c)\tutg%.6dl(%c)\tnse::%lu\n", - // ((int32_t)(v>>1))+1, "+-"[v&1], - // ((int32_t)(w>>1))+1, "+-"[w&1], nse); gidx[i] |= ((uint64_t)0x8000000000000000); return; } } - v = w; + vk = wk; } } @@ -12382,7 +12426,7 @@ uint8_t *ff, uint32_t *ng_occ, uint32_t max_ext) for (i = 0, ub64->n = 0; i < res.n; i++) { s = res.a[i]>>32; e = ((uint32_t)res.a[i]); - update_usg_t_threading_0(ng, int_a + s, e + 1 - s, ng_occ, ub64); + update_usg_t_threading_0(uidx, ng, int_a + s, e + 1 - s, ng_occ, ub64); } } kv_destroy(res); @@ -13011,7 +13055,8 @@ void renew_usg_t_bub(ul_resolve_t *uidx, usg_t *ng, uint32_t *id_map, uint8_t *f // fprintf(stderr, "[M::%s::k->%u] utg%.6dl(%c), ug->u.a[k].n::%u, rocc_cut::%u, is_hom::%u\n", __func__, k, // i+1, "+-"[0], (uint32_t)ug->u.a[k].n, rocc_cut, IF_HOM(k, (*(uidx->bub)))); // } - if((ug->u.a[k].n >= rocc_cut) && (!(IF_HOM(k, (*(uidx->bub)))))) { + if((ug->u.a[k].n >= rocc_cut) && + ((!(IF_HOM(k, (*(uidx->bub))))) || (asm_opt.purge_level_primary == 0))) { if(id_map[i]&p[0]) { ff[i<<1] = 2; // fprintf(stderr, "[M::%s::] utg%.6dl(%c), ug->u.a[k].n::%u, rocc_cut::%u\n", __func__, @@ -13110,8 +13155,8 @@ asg64_v *b0, asg64_v *b1, double cutoff) // uint64_t is_debug = 0; - // if(((v>>1) == 10600)/** && (b0->n > 1 && (b0->a[b0->n-1]>>1) == 34944)**/) { - // is_debug = 1; + // if((((v>>1) == 2736) && (v&1))/** && (b0->n > 1 && (b0->a[b0->n-1]>>1) == 34944)**/) { + // // is_debug = 1; // fprintf(stderr, "[M::%s::] utg%.6dl(%c)\n", __func__, // (int32_t)(v>>1)+1, "+-"[v&1]); // prt_sub_integer_path(b0, it, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1); @@ -13130,12 +13175,13 @@ asg64_v *b0, asg64_v *b1, double cutoff) else return 0;///b0 is contained if(z < b1->n) get_integer_seq_ovlps(uidx, b1->a, b1->n, z - 1, 0, NULL, &w1); else continue;///b1 is contained - // if(((v>>1) == 8722) || ((v>>1) == 56768)) { + // if(((v>>1) == 2736) && (v&1)) { // prt_sub_integer_path(b0, it, w0, w1, zn, z); // prt_sub_integer_path(b1, it, w0, w1, zn, z); // } if(w0 == (uint64_t)-1) w0 = 0; if(w1 == (uint64_t)-1) w1 = 0; + if((w0 <= w1) || (w1 > (w0*cutoff))) return 0; // if(b0->n == zn) return 0;///b0 is shorter // if(b1->n == zn) continue;///b1 is shorter if((min_w0 == (uint64_t)-1) || (z == zn) || (min_w0 > w0) || (min_w0 == w0 && min_w1 < w1)) { @@ -13143,6 +13189,10 @@ asg64_v *b0, asg64_v *b1, double cutoff) } is_contain = 0; } + + // if(((v>>1) == 2736) && (v&1)) { + // fprintf(stderr, "[M::%s::] min_w0::%lu, min_w1::%lu\n\n", __func__, min_w0, min_w1); + // } // fprintf(stderr, "[M::%s::] utg%.6dl(%c)->utg%.6dl(%c)\n", __func__, // (int32_t)(v>>1)+1, "+-"[v&1], (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1]); if(alt_n == 0) return 1; @@ -13165,13 +13215,13 @@ uint64_t *ridx, asg64_v *res) v = int_a[k]; if((!f[v])&&(!f[v^1])) continue; if((pi != (uint64_t)-1) && (f[v^1])) { - // if(((v>>1) == 8722) || ((v>>1) == 56768)) { + // if(((v>>1) == 2736)) { // fprintf(stderr, "+[M::%s::ii[%lu, %lu)] utg%.6dl(%c), f[v^1]::%u, v^1::%lu, putg%.6dl(%c), f[pv]::%u, pv::%lu\n", // __func__, s, e, (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1], f[v^1], v^1, // (int32_t)(int_a[pi]>>1)+1, "+-"[int_a[pi]&1], f[int_a[pi]], int_a[pi]); // } if(is_best_path(uidx, ng, int_idx, int_a, s, e, k, v^1, ridx_a, ridx, b0, b1, 0.51)) { - // if(((v>>1) == 8722) || ((v>>1) == 56768)) { + // if(((v>>1) == 2736)) { // fprintf(stderr, "-[M::%s::ii[%lu, %lu)] utg%.6dl(%c), f[v^1]::%u, v^1::%lu, putg%.6dl(%c), f[pv]::%u, pv::%lu\n", // __func__, s, e, (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1], f[v^1], v^1, // (int32_t)(int_a[pi]>>1)+1, "+-"[int_a[pi]&1], f[int_a[pi]], int_a[pi]); @@ -13199,6 +13249,81 @@ uint64_t *ridx, asg64_v *res) } } + +uint64_t gen_sub_integer_path_ff(uint64_t *int_a, int64_t s, int64_t e, int64_t it, uint64_t v, uint8_t *f, asg64_v *res) +{ + int64_t k; res->n = 0; assert((int_a[it]>>1) == (v>>1)); + if(int_a[it] == v) {///forward + kv_resize(uint64_t, *res, (uint64_t)(e - it)); + for (k = it; k < e; k++) { + res->a[res->n++] = int_a[k]; + if(res->n > 1 && f[res->a[res->n-1]^1]) return res->n; + } + } else {//reverse + kv_resize(uint64_t, *res, (uint64_t)(it - s)); + for (k = it; k >= s; k--) { + res->a[res->n++] = int_a[k]^1; + if(res->n > 1 && f[res->a[res->n-1]^1]) return res->n; + } + } + res->n = 0; return 0; +} + +uint64_t is_best_pair(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t *int_a, +uint64_t s, uint64_t e, uint64_t it, uint64_t v, uint64_t *ridx_a, uint64_t *ridx, +asg64_v *b0, asg64_v *b1, uint8_t *f, double cutoff) +{ + uint64_t *arc_a, arc_n, k, is, ie, w0 = 0, w1 = 0, sup_cut = 3; + if(!gen_sub_integer_path_ff(int_a, s, e, it, v, f, b0)) return 0; + assert(b0->n >= 2); assert(f[b0->a[0]] && f[b0->a[b0->n-1]^1]); + w0 = get_arc_support_chain(uidx, b0->a, b0->n, ng); + if(w0 < sup_cut) return 0; + + arc_a = ridx_a + (ridx[v>>1]>>32); + arc_n = (uint32_t)ridx[v>>1]; + for (k = 0; k < arc_n; k++) { + if(it == ((uint32_t)arc_a[k])) continue; + is = int_idx[arc_a[k]>>32]>>32; + ie = is + ((uint32_t)(int_idx[arc_a[k]>>32])); + if(!gen_sub_integer_path_ff(int_a, is, ie, ((uint32_t)arc_a[k]), v, f, b1)) continue; + assert(b1->n >= 2); assert(f[b1->a[0]] && f[b1->a[b1->n-1]^1]); + if((b0->n == b1->n) && (memcmp(b0->a, b1->a, sizeof((*(b0->a)))*b0->n) == 0)) { + if(it < ((uint32_t)arc_a[k])) continue; + else return 0;///only keep 1 equal interval + } + w1 = get_arc_support_chain(uidx, b1->a, b1->n, ng); + // fprintf(stderr, "[M::%s::] w0::%lu, w1::%lu\n", __func__, w0, w1); + if((w0 <= w1) || (w1 > (w0*cutoff))) return 0; + } + return 1; +} + +void get_best_pair(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t *int_a, +uint8_t *f, uint64_t s, uint64_t e, asg64_v *b0, asg64_v *b1, uint64_t *ridx_a, +uint64_t *ridx, asg64_v *res) +{ + uint64_t pi, v, k, *pz; + b0->n = b1->n = 0; + for (k = s, pi = (uint64_t)-1; k < e; k++) { + v = int_a[k]; + if((f[v^1])) { + if(pi != (uint64_t)-1) { + if(is_best_pair(uidx, ng, int_idx, int_a, s, e, k, v^1, ridx_a, ridx, b0, b1, f, 0.51)) { + kv_pushp(uint64_t, *res, &pz); + *pz = pi; (*pz) <<= 32; (*pz) |= k; + } + } + pi = (uint64_t)-1; + } + if(f[v]) { + pi = (uint64_t)-1; + if(is_best_pair(uidx, ng, int_idx, int_a, s, e, k, v, ridx_a, ridx, b0, b1, f, 0.51)) { + pi = k;//keep the shortest pi<->k + } + } + } +} + static void worker_ul_aln_path(void *data, long i, int tid) // callback for kt_for() { unique_bridge_check_t *u_aux = (unique_bridge_check_t*)data; @@ -13226,6 +13351,33 @@ static void worker_ul_aln_path(void *data, long i, int tid) // callback for kt_f buf->res_dump.a = res.a; buf->res_dump.n = res.n; buf->res_dump.m = res.m; } +static void worker_ul_aln_pair(void *data, long i, int tid) // callback for kt_for() +{ + unique_bridge_check_t *u_aux = (unique_bridge_check_t*)data; + ul_resolve_t *uidx = u_aux->uidx; + integer_t *buf = &(uidx->str_b.buf[tid]); + uint64_t s, e; + // uint64_t *x = &(uidx->uovl.iug_tra->a[i]); + // asg_arc_t *ve = &(uidx->uovl.i_ug->g->arc[*x]); + + asg64_v b_v, b_r, res; + b_v.a = buf->u.a; b_v.n = buf->u.n; b_v.m = buf->u.m; + b_r.a = buf->o.a; b_r.n = buf->o.n; b_r.m = buf->o.m; + res.a = buf->res_dump.a; res.n = buf->res_dump.n; res.m = buf->res_dump.m; + + b_v.n = b_r.n = 0; + s = u_aux->int_idx[i]>>32; ///the i-th integer contig/path + e = s + ((uint32_t)(u_aux->int_idx[i])); + assert(e > s + 1);//the length is at least 2 + + get_best_pair(uidx, u_aux->ng, u_aux->int_idx, u_aux->int_a, u_aux->f, s, e, &b_v, &b_r, u_aux->ridx_a, u_aux->ridx, &res); + // get_ul_arc_supports(uidx, ve, &b_v, &b_r, 1, &w_v, &w_r); + + buf->u.a = b_v.a; buf->u.n = b_v.n; buf->u.m = b_v.m; + buf->o.a = b_r.a; buf->o.n = b_r.n; buf->o.m = b_r.m; + buf->res_dump.a = res.a; buf->res_dump.n = res.n; buf->res_dump.m = res.m; +} + uint32_t ug_ext_0(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, uint32_t max_ext, uint8_t *ff, uint32_t *ng_occ, uint64_t *i_idx, asg64_v *b64, asg64_v *ub64, uint64_t a_n) { @@ -13294,7 +13446,8 @@ 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) +uint8_t **ff, uint32_t **ng_occ, uint64_t **i_idx, asg64_v *b64, asg64_v *ub64, uint32_t rocc_cut, +uint32_t thread_path) { // 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; @@ -13313,7 +13466,13 @@ uint8_t **ff, uint32_t **ng_occ, uint64_t **i_idx, asg64_v *b64, asg64_v *ub64, for (k = 0; k < uidx->str_b.n_thread; k++) { uidx->str_b.buf[k].res_dump.n = uidx->str_b.buf[k].u.n = uidx->str_b.buf[k].o.n = 0; } - kt_for(uidx->str_b.n_thread, worker_ul_aln_path, &u_aux, u_aux.int_idx_n); + + if(thread_path) { + kt_for(uidx->str_b.n_thread, worker_ul_aln_path, &u_aux, u_aux.int_idx_n); + } else { + kt_for(uidx->str_b.n_thread, worker_ul_aln_pair, &u_aux, u_aux.int_idx_n); + } + for (k = b64->n = a_n = 0; k < uidx->str_b.n_thread; k++) { a_n += uidx->str_b.buf[k].res_dump.n; kv_resize(uint64_t, *b64, a_n); @@ -14229,9 +14388,11 @@ void u2g_hybrid_detan_iter(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint // prt_usg_t(uidx, ng, "ng0"); for (k = 0; k < clean_round; k++) { ncut += 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, 48); - ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 16); + ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 48, 1); + // prt_usg_t(uidx, ng, "ng_python"); + ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 16, 1); + ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 0, 0); ///renew bubble for ug_ext_strict n_vtx = ng->n<<1; REALLOC(ff, n_vtx); REALLOC(ng_occ, n_vtx); REALLOC(i_idx, n_vtx); diff --git a/hic.cpp b/hic.cpp index ee3e875..36caaba 100644 --- a/hic.cpp +++ b/hic.cpp @@ -285,7 +285,6 @@ KRADIX_SORT_INIT(k_trans_qs, u_trans_t, k_trans_qs_key, member_size(u_trans_t, q #define get_hit_euid(x, k) (((x).a.a[(k)].e<<1)>>(64 - (x).uID_bits)) #define get_hit_epos(x, k) ((x).a.a[(k)].e & (x).pos_mode) - typedef struct { kvec_t(kvec_t_u64_warp) matrix; uint64_t uID_shift, dis_mode; @@ -17551,4 +17550,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 ce3f22f..57af6af 100644 --- a/inter.cpp +++ b/inter.cpp @@ -25,6 +25,9 @@ void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint6 int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t high_occ, void *km); void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off); +void ug_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, + int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, uint32_t is_hpc); + #define MG_SEED_IGNORE (1ULL<<41) #define MG_SEED_TANDEM (1ULL<<42) #define MG_SEED_KEPT (1ULL<<43) @@ -248,8 +251,6 @@ KRADIX_SORT_INIT(eg_srt_x, eg_srt_t, eg_srt_x_key, member_size(eg_srt_t, x)) #define eg_srt_d_key(p) ((p).d) KRADIX_SORT_INIT(eg_srt_d, eg_srt_t, eg_srt_d_key, member_size(eg_srt_t, d)) - - // shortest path typedef struct { // input @@ -374,6 +375,28 @@ typedef struct { // data structure for each step in kt_pipeline() } utepdat_t; +typedef struct { // global data structure for kt_pipeline() + uint64_t x, i; +} ha_mzl_srt_t; + +typedef struct { size_t n, m; ha_mzl_srt_t *a; } ha_mzl_srt_v; + +typedef struct { // global data structure for kt_pipeline() + ma_ug_t *ug; + asg_t *rg; + ha_ovec_buf_t **hab; + glchain_t *ll; + int32_t w, k; + ///base-alignment + double bw_thres, diff_ec_ul; + int32_t max_n_chain; + ha_mzl_v idx_a; + asg64_v idx_n; + ha_mzl_srt_v srt_a; + int32_t is_HPC, bw, max_gap, chn_pen_gap, n_thread, is_cnt; + ul_idx_t udb; +} ug_trans_t; + void hc_glchain_destroy(glchain_t *b) { if (!b) return; @@ -9184,6 +9207,126 @@ uint32_t ck_ul_alignment(ul_vec_t *x) } +static void worker_for_trans_ovlp(void *data, long i, int tid) // callback for kt_for() +{ + /** + ug_trans_t *s = (ug_trans_t*)data; + ha_ovec_buf_t *b = s->hab[tid]; + glchain_t *bl = &(s->ll[tid]); + int64_t winLen = MIN((((double)THRESHOLD_MAX_SIZE)/s->diff_ec_ul), WINDOW), ton = 0; + uint32_t high_occ = 2, phase = 1, k; + asg64_v b0, b1, b2; window_list p; memset(&p, 0, sizeof(p)); + overlap_region *aux_o = NULL; + char *seq = s->ug->u.a[i].s; + int64_t len = s->ug->u.a[i].len; + **/ + // uint64_t align = 0; + + // if(UL_INF.a[s->id+i].rlen != s->len[i]) { + // fprintf(stderr, "[M::%s] rid:%ld, s->len:%lu, UL_INF->rlen:%u\n", __func__, s->id+i, s->len[i], UL_INF.a[s->id+i].rlen); + // } + // void *km = s->buf?(s->buf[tid]?s->buf[tid]->km:NULL):NULL; + // if(s->id+i!=3046) return; + // if((s->id+i!=871) && (s->id+i!=963) && (s->id+i!=980)) return; + // if(s->id+i!=944) return; + // if(s->id+i != 35437) return; + // if((s->id+i != 7086) && (s->id+i != 51705) && (s->id+i != 266022) && (s->id+i != 353608) + // && (s->id+i != 399416) && (s->id+i != 403014) && (s->id+i != 420915) && (s->id+i != 603855) + // && (s->id+i != 680134) && (s->id+i != 766261) && (s->id+i != 794527)) { + // return; + // } + // char *as = NULL; + // asprintf(&as, "\n[M::%s] rid::%ld, len::%lu, name::%.*s\n", __func__, s->id+i, s->len[i], (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a); + // push_vlog(&(overall_zdbg->a[s->id+i]), as); free(as); as = NULL; + // fprintf(stderr, "\n[M::%s] rid::%ld, len::%lu, name::%.*s\n", __func__, s->id+i, s->len[i], + // (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a); + // fprintf(stderr, ">%.*s\n%.*s\n", (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a, + // (int32_t)s->len[i], s->seq[i]); + + // if (memcmp(UL_INF.nid.a[s->id+i].a, "d0aab024-b3a7-40fb-83cc-22c3d6d951f8", UL_INF.nid.a[s->id+i].n-1)) return; + // fprintf(stderr, "[M::%s::] ==> len: %lu\n", __func__, s->len[i]); + // ha_get_ul_candidates_interface(b->abl, i, s->seq[i], s->len[i], s->opt->w, s->opt->k, s->uu, &b->olist, &b->olist_hp, &b->clist, s->opt->bw_thres, + // s->opt->max_n_chain, 1, NULL, &b->r_buf, &(b->tmp_region), NULL, &(b->sp), 1, NULL); + /** + ug_map_lchain(b->abl, (uint32_t)-1, seq, len, s->w, s->k, &(s->udb), NULL, NULL, s->bw_thres, + s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, s->is_HPC); + **/ + /** + clear_Cigar_record(&b->cigar1); + clear_Round2_alignment(&b->round2); + // return; + // b->num_correct_base += overlap_statistics(&b->olist, NULL, 0); + + // int fully_cov, abnormal; + // b->self_read.seq = s->seq[i]; b->self_read.length = s->len[i]; b->self_read.size = 0; + // correct_ul_overlap(&b->olist, s->uu, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon, + // &b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul, winLen, NULL); + // memset(&b->self_read, 0, sizeof(b->self_read)); + + ul_lalign(&b->olist, &b->clist, s->uu, s->uopt, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, + &b->correct, &b->exz, &b->hap, &b->r_buf, aux_o, s->opt->diff_ec_ul, winLen, NULL, s->id+i, s->opt->k, &(s->sps[tid]), NULL); + // ul_lalign_old_ed(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, + // &b->correct, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, 1, NULL); + ton = b->olist.length;//all alignments pass similary check + aux_o = gen_aux_ovlp(&b->olist);///must be here + gl_chain_flter(&b->olist, &b->correct, &(s->sps[tid]), bl, s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, &phase); + + // fprintf(stderr, "[M::%s] rid::%ld, len::%lu, name::%.*s, phase::%u\n", __func__, s->id+i, s->len[i], + // (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a, phase); + if(phase && gen_shared_intervals(&b->olist, s->uu, s->uopt, winLen, &b->r_buf, &(bl->lo))) { + filter_topN(&b->olist, &(bl->lo), s->len[i], winLen, UL_TOPN, bl); + // update_shared_intervals(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b->r_buf, &(s->sps[tid]), s->len[i], winLen, &(bl->lo), s->id+i); + + copy_asg_arr(b0, b->hap.snp_srt); copy_asg_arr(b1, s->sps[tid]); copy_asg_arr(b2, b->r_buf.a); + // update_shared_intervals(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b0, &b1, &b2, s->len[i], winLen, &(bl->lo), s->id+i); + update_sketch_trace(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b0, &b1, &b2, s->len[i], winLen, &(bl->lo), s->id+i, MAX_LGAP(s->len[i]), s->opt->diff_ec_ul); + copy_asg_arr(b->hap.snp_srt, b0); copy_asg_arr(s->sps[tid], b1); copy_asg_arr(b->r_buf.a, b2); + + ul_lalign(&b->olist, &b->clist, s->uu, s->uopt, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, + &b->correct, &b->exz, &b->hap, &b->r_buf, aux_o, s->opt->diff_ec_ul, winLen, &(bl->lo), s->id+i, s->opt->k, &(s->sps[tid]), NULL); + // ul_lalign_old_ed(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, + // &b->correct, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, 0, NULL); + ///recover alignments + for (k = b->olist.length; k < ton; k++) { + b->olist.list[k].w_list.n = 0; + p.x_start = b->olist.list[k].x_pos_s; + p.x_end = b->olist.list[k].x_pos_e+1; + p.clen = b->olist.list[k].non_homopolymer_errors; + kv_push(window_list, b->olist.list[k].w_list, p); + b->olist.list[k].align_length = 0; + b->olist.list[k].overlapLen = b->olist.list[k].x_pos_e+1-b->olist.list[k].x_pos_s; + } + b->olist.length = ton; + } else { + for (k = 0; k < b->olist.length; k++) { + b->olist.list[k].w_list.n = 0; + p.x_start = b->olist.list[k].x_pos_s; + p.x_end = b->olist.list[k].x_pos_e+1; + p.clen = 0; + kv_push(window_list, b->olist.list[k].w_list, p); + b->olist.list[k].align_length = b->olist.list[k].overlapLen = + b->olist.list[k].x_pos_e+1-b->olist.list[k].x_pos_s; + b->olist.list[k].non_homopolymer_errors = 0; + } + } + + // prt_overlap_region_alloc_ol(&(b->olist), s->id+i); + + + gl_chain(s->buf[tid], &(UL_INF.a[s->id+i]), &b->olist, &(b->clist.chainDP), &b->hap, &(s->sps[tid]), bl, &(s->gdp[tid]), s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, s->id+i, tid, NULL); + + + 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++; + **/ +} + + static void worker_for_ul_recorrect_alignment(void *data, long i, int tid) // callback for kt_for() { utepdat_t *s = (utepdat_t*)data; @@ -14891,6 +15034,49 @@ void clear_all_ul_t(all_ul_t *x) free(x->ridx.idx.a); free(x->ridx.occ.a); memset(&(x->ridx), 0, sizeof((x->ridx))); } +void init_ug_trans_t(ug_trans_t *opt, int32_t is_HPC, int32_t k, int32_t w, int32_t max_n_chain, +double bw_thres, double diff_ec_ul, int32_t n_thread, ma_ug_t *ug, asg_t *sg) +{ + int64_t i; + memset(opt, 0, sizeof((*opt))); + opt->k = k; + opt->w = w; + opt->is_HPC = is_HPC; + opt->max_n_chain = max_n_chain; + opt->bw_thres = bw_thres; + opt->diff_ec_ul = diff_ec_ul; + + opt->bw = 10000;///2000 in minigraph + opt->max_gap = 500000;///5000 in minigraph + opt->ug = ug; + opt->rg = sg; + opt->n_thread = ((n_thread>=1)?n_thread:1); + + CALLOC(opt->hab, opt->n_thread); + CALLOC(opt->ll, opt->n_thread); + for (i = 0; i < opt->n_thread; ++i) { + opt->hab[i] = ha_ovec_init(0, 0, 1); + } + + opt->idx_n.n = opt->idx_n.m = ug->u.n; + MALLOC(opt->idx_n.a, opt->idx_n.n); + + opt->udb.ug = ug; +} + +void trans_base_count(ug_trans_t *p) +{ + p->is_cnt = 1; ha_flt_tab = NULL; + kt_for(p->n_thread, worker_for_trans_ovlp, p, p->ug->u.n); + +} + +void trans_base_infer(ma_ug_t *ug, asg_t *sg) +{ + ug_trans_t sl; + init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); + init_ug_trans_t(&sl, 0, asm_opt.trans_mer_length, asm_opt.trans_win, asm_opt.max_n_chain, 1.0-asm_opt.trans_base_rate, 1.0-asm_opt.trans_base_rate, asm_opt.thread_num, ug, sg); +} ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_cache)