From 86b7dd424ae0674f7b0ab099b1299cb9a6ca4f6b Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Wed, 22 Mar 2023 12:23:53 -0400 Subject: [PATCH] avoid misassemblies; better polyploidy graph --- CommandLines.cpp | 2 +- CommandLines.h | 2 +- Correct.cpp | 85 +++- Overlaps.cpp | 1001 +++++++++++++++++++++++++++++++++++++++++++++- Overlaps.h | 23 ++ anchor.cpp | 18 +- gfa_ut.cpp | 41 +- gfa_ut.h | 2 +- hic.cpp | 279 ++++++++++++- hic.h | 2 +- inter.cpp | 554 +++++++++++++++++++++++-- inter.h | 1 + 12 files changed, 1931 insertions(+), 79 deletions(-) diff --git a/CommandLines.cpp b/CommandLines.cpp index d9fa080..fda2142 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -285,7 +285,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->min_path_drop_rate = 0.2; asm_opt->max_path_drop_rate = 0.6; asm_opt->hifi_pst_join = 1; - asm_opt->ul_pst_join = 0; + asm_opt->ul_pst_join = 1; } void destory_enzyme(enzyme* f) diff --git a/CommandLines.h b/CommandLines.h index d172968..a60c268 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.19.2-r560" +#define HA_VERSION "0.19.2-r572" #define VERBOSE 0 diff --git a/Correct.cpp b/Correct.cpp index 216131b..34c4450 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -15216,8 +15216,10 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u else ibeg = ov->qn; iend = ov->tn; assert(iend>=ibeg+1); - // fprintf(stderr, "\n***[M::%s::rid->%lu] utg%.6dl(%c), s::%u, e::%u, z::[%u, %u), ibeg::%ld, iend::%ld, ch_n::%ld\n", - // __func__, rid, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ov->qs, ov->qe, z->x_pos_s, z->x_pos_e+1, ibeg, iend, ch_n); + // if(z->x_id == 57 && z->y_id == 2175) { + // fprintf(stderr, "\n***[M::%s::rid->%lu] utg%.6dl(%c), s::%u, e::%u, z::[%u, %u), ibeg::%ld, iend::%ld, ch_n::%ld\n", + // __func__, rid, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ov->qs, ov->qe, z->x_pos_s, z->x_pos_e+1, ibeg, iend, ch_n); + // } for (l = ibeg, i = ibeg + 1; i <= iend; i++) { q[0] = q[1] = t[0] = t[1] = mode = -1; is_done = 0; if(l >= 0) { @@ -15243,14 +15245,14 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u } if(mode == 1 || mode == 2) adjust_ext_offset(&(q[0]), &(q[1]), &(t[0]), &(t[1]), ql, tl, 0, mode); - // if(aux_o->x_id == 29033 && aux_o->y_id == 21307) { + // if(z->x_id == 57 && z->y_id == 2175) { // fprintf(stderr, "#[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n", // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); // } is_done = hc_aln_exz_adv(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, ql, e_rate, MAX_SIN_L, MAX_SIN_E, FORCE_SIN_L, -1, aux_o); - // if(aux_o->x_id == 29033 && aux_o->y_id == 21307) { + // if(z->x_id == 57 && z->y_id == 2175) { // fprintf(stderr, "-is_done::%ld[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n", // is_done, __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); // } @@ -19123,11 +19125,40 @@ uint64_t gen_sub_ov_adv(const ul_idx_t *udb, overlap_region* o, double o_rate, u int64_t return_t_chain(overlap_region *z, Candidates_list *cl) { int64_t i, cn = cl->length, scn; uint64_t pid; k_mer_hit *ca; + + // if(z->x_id == 57 && z->y_id == 2175) { + // fprintf(stderr, "\n-0-[M::%s]\tutg%.6ul\tx::[%u,\t%u)\t%c\tutg%.6ul\ty::[%u,\t%u)\n", + // __func__, z->x_id+1, z->x_pos_s, z->x_pos_e+1, + // "+-"[z->y_pos_strand], z->y_id+1, z->y_pos_s, z->y_pos_e+1); + // i = z->shared_seed; pid = cl->list[i].readID; + // for (; i < cn && cl->list[i].readID == pid && cl->list[i].readID != ((uint32_t)(0x7fffffff)); i++) { + // fprintf(stderr, "i::%ld[M::%s]\treadID::%u\tself_offset::%u\toffset::%u\t%c\n", + // i, __func__, cl->list[i].readID, cl->list[i].self_offset, cl->list[i].offset, + // "+-"[cl->list[i].strand]); + // } + // } + + + i = z->shared_seed; pid = cl->list[i].readID; for (; i < cn && cl->list[i].readID == pid && cl->list[i].readID != ((uint32_t)(0x7fffffff)); i++); scn = i - z->shared_seed; ca = cl->list+z->shared_seed; i = lchain_refine(ca, scn, ca, &(cl->chainDP), 50, 5000, 512, 16); cn = i; for (; i < scn; i++) ca[i].readID = ((uint32_t)(0x7fffffff)); + + + + // if(z->x_id == 57 && z->y_id == 2175) { + // fprintf(stderr, "\n-1-[M::%s]\tutg%.6ul\tx::[%u,\t%u)\t%c\tutg%.6ul\ty::[%u,\t%u)\n", + // __func__, z->x_id+1, z->x_pos_s, z->x_pos_e+1, + // "+-"[z->y_pos_strand], z->y_id+1, z->y_pos_s, z->y_pos_e+1); + // i = z->shared_seed; pid = cl->list[i].readID; + // for (; i < cn && cl->list[i].readID == pid && cl->list[i].readID != ((uint32_t)(0x7fffffff)); i++) { + // fprintf(stderr, "i::%ld[M::%s]\treadID::%u\tself_offset::%u\toffset::%u\t%c\n", + // i, __func__, cl->list[i].readID, cl->list[i].self_offset, cl->list[i].offset, + // "+-"[cl->list[i].strand]); + // } + // } return cn; } @@ -21017,6 +21048,20 @@ bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t tl = udb->ug->u.a[id].len; for (i = ch_idx; i < cl->length && cl->list[i].readID == cl->list[ch_idx].readID; i++); ch_n = i-ch_idx; + // if(z->x_id == 57 && z->y_id == 2175) { + // fprintf(stderr, "\n-*-[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\tch_idx::%ld\tch_n::%ld\n", + // __func__, + // z->x_id+1, "lc"[udb->ug->u.a[z->x_id].circ], udb->ug->u.a[z->x_id].len, z->x_pos_s, z->x_pos_e+1, + // "+-"[z->y_pos_strand], + // z->y_id+1, "lc"[udb->ug->u.a[z->y_id].circ], udb->ug->u.a[z->y_id].len, z->y_pos_s, z->y_pos_e+1, + // ch_idx, ch_n); + // for (i = ch_idx; i < cl->length && cl->list[i].readID == cl->list[ch_idx].readID; i++) { + // fprintf(stderr, "i::%ld[M::%s]\treadID::%u\tself_offset::%u\toffset::%u\t%c\n", + // i, __func__, cl->list[i].readID, cl->list[i].self_offset, cl->list[i].offset, + // "+-"[cl->list[i].strand]); + // } + // } + // on = fusion_chain_ovlp(z, ch_a, ch_n, ov, on, wl, ql, tl); aux_o->w_list.n = aux_o->w_list.c.n = 0; aux_o->y_id = z->y_id; aux_o->y_pos_strand = z->y_pos_strand; @@ -21033,8 +21078,8 @@ bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t int64_t aux_n = aux_o->w_list.n; - // if(z->x_id == 77960 && z->y_id == 27340) { - // fprintf(stderr, "\n[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\n", + // if(z->x_id == 57 && z->y_id == 2175) { + // fprintf(stderr, "\n-0-[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\n", // __func__, // z->x_id+1, "lc"[udb->ug->u.a[z->x_id].circ], udb->ug->u.a[z->x_id].len, z->x_pos_s, z->x_pos_e+1, // "+-"[z->y_pos_strand], @@ -21057,8 +21102,8 @@ bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t rechain_aln(z, cl, aux_o, i, w_l, udb, NULL, NULL, qstr, tu, exz, e_rate, ql, tl, khit, rid); } - // if(z->x_id == 77960 && z->y_id == 27340) { - // fprintf(stderr, "\n[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\n", + // if(z->x_id == 57 && z->y_id == 2175) { + // fprintf(stderr, "\n-1-[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\n", // __func__, // z->x_id+1, "lc"[udb->ug->u.a[z->x_id].circ], udb->ug->u.a[z->x_id].len, z->x_pos_s, z->x_pos_e+1, // "+-"[z->y_pos_strand], @@ -21085,6 +21130,7 @@ bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t ///update z by aux_o update_overlap_region(z, aux_o, ql, tl); + assert((z->x_pos_e>=z->x_pos_s) && (z->y_pos_e>=z->y_pos_s)); int64_t zwn = z->w_list.n, zerr = 0, zlen = z->x_pos_e+1-z->x_pos_s; for (i = 0; i < zwn; i++) { @@ -21097,8 +21143,8 @@ bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t } z->non_homopolymer_errors = zerr; - // if(z->x_id == 77960 && z->y_id == 27340) { - // fprintf(stderr, "\n[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\n", + // if(z->x_id == 57 && z->y_id == 2175) { + // fprintf(stderr, "\n-2-[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\n", // __func__, // z->x_id+1, "lc"[udb->ug->u.a[z->x_id].circ], udb->ug->u.a[z->x_id].len, z->x_pos_s, z->x_pos_e+1, // "+-"[z->y_pos_strand], @@ -21114,6 +21160,15 @@ bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t // fprintf(stderr, "\n"); // } + // if(z->x_id == 57 && z->y_id == 2175) { + // fprintf(stderr, "\n-3-[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\tzerr::%ld\tzlen::%ld\te_rate::%f\n", + // __func__, + // z->x_id+1, "lc"[udb->ug->u.a[z->x_id].circ], udb->ug->u.a[z->x_id].len, z->x_pos_s, z->x_pos_e+1, + // "+-"[z->y_pos_strand], + // z->y_id+1, "lc"[udb->ug->u.a[z->y_id].circ], udb->ug->u.a[z->y_id].len, z->y_pos_s, z->y_pos_e+1, + // zerr, zlen, e_rate); + // } + if(zerr >= zlen) return 0; if(zerr <= 0 && zlen > 0) return 1; if(zerr > (zlen*e_rate)) return 0; @@ -21319,7 +21374,7 @@ void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur rr = gen_extend_err_exz(z, uref, NULL, NULL, in/**qu->seq**/, tu->seq, exz, NULL/**v_idx?v_idx->a.a:NULL**/, w.window_length, -1, err, (e_max+0.000001), &re); z->is_match = 0;///must be here; - // if(z->x_id == 77960 && z->y_id == 27340) { + // if(z->x_id == 57 && z->y_id == 2175) { // fprintf(stderr, "+utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\talign_length::%u\terr::%f\trr::%f\twn::%u\tk::%lu\ti::%lu\n", // z->x_id+1, "lc"[uref->ug->u.a[z->x_id].circ], uref->ug->u.a[z->x_id].len, // z->x_pos_s, z->x_pos_e+1, "+-"[z->y_pos_strand], @@ -21339,7 +21394,7 @@ void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur } ol->length = k; - // if(sid == 77960) { + // if(sid == 57) { // fprintf(stderr, "+utg%.6ld%c\tol->length::%lu\n", sid+1, "lc"[uref->ug->u.a[sid].circ], ol->length); // } if(ol->length <= 0) return; @@ -21351,7 +21406,7 @@ void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur } for (i = k = 0; i < ol->length; i++) { z = &(ol->list[i]); - // if(z->x_id == 77960 && z->y_id == 27340) { + // if(z->x_id == 57 && z->y_id == 2175) { // fprintf(stderr, "-1-utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\talign_length::%u\terr::%f\twn::%u\tk::%lu\ti::%lu\n", // z->x_id+1, "lc"[uref->ug->u.a[z->x_id].circ], uref->ug->u.a[z->x_id].len, // z->x_pos_s, z->x_pos_e+1, "+-"[z->y_pos_strand], @@ -21364,7 +21419,7 @@ void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur continue; } - // if(z->x_id == 77960 && z->y_id == 27340) { + // if(z->x_id == 57 && z->y_id == 2175) { // fprintf(stderr, "-2-utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\talign_length::%u\terr::%f\twn::%u\tk::%lu\ti::%lu\n", // z->x_id+1, "lc"[uref->ug->u.a[z->x_id].circ], uref->ug->u.a[z->x_id].len, // z->x_pos_s, z->x_pos_e+1, "+-"[z->y_pos_strand], @@ -21382,7 +21437,7 @@ void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur } ol->length = k; - // if(sid == 77960) { + // if(sid == 57) { // fprintf(stderr, "-utg%.6ld%c\tol->length::%lu\n", sid+1, "lc"[uref->ug->u.a[sid].circ], ol->length); // } } diff --git a/Overlaps.cpp b/Overlaps.cpp index 1b9d22c..0372f9d 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -14121,8 +14121,19 @@ static void worker_for_trans_clean_re(void *data, long i, int tid) // callback f z = cz->qe; cz->qe = cz->te; cz->te = z; } + // if(id == 394 || id == 1698 || id == 84) { + // fprintf(stderr, "+[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\n", __func__, + // cz->qn+1, "lc"[s->ug->u.a[cz->qn].circ], s->ug->u.a[cz->qn].len, cz->qs, cz->qe, "+-"[cz->rev], + // cz->tn+1, "lc"[s->ug->u.a[cz->tn].circ], s->ug->u.a[cz->tn].len, cz->ts, cz->te); + // } + if((cz->nw >= 0) && ((!trans_ovlp_connect1(cz, s->ug)) - || (!ovlp_rocc(cz, s->ug, s->rg, s->small_ov_rate, s->ov_cutoff)))) { + || (!ovlp_rocc(cz, s->ug, s->rg, s->small_ov_rate, s->ov_cutoff)))) { + // if(id == 394 || id == 1698 || id == 84) { + // fprintf(stderr, "-[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\n", __func__, + // cz->qn+1, "lc"[s->ug->u.a[cz->qn].circ], s->ug->u.a[cz->qn].len, cz->qs, cz->qe, "+-"[cz->rev], + // cz->tn+1, "lc"[s->ug->u.a[cz->tn].circ], s->ug->u.a[cz->tn].len, cz->ts, cz->te); + // } res->n--; } } @@ -14396,6 +14407,90 @@ static void worker_for_trans_sec_cut(void *data, long i, int tid) // callback fo } } + +uint64_t trans_sec_cut0(kv_u_trans_t *ta, asg64_v *srt, uint32_t id, double sec_rate, uint64_t bd, ma_ug_t *ug) +{ + u_trans_t *a = NULL, *mz, *cz, t, *ca; + uint32_t a_n, r_n, c_n, k, l, z, i, m, wm, wl, len, b_n = srt->n; + uint32_t ovq, os, oe; uint64_t *ss, *hs, *bu, *wu, occ = 0; + + a = u_trans_a(*ta, id); a_n = u_trans_n(*ta, id); + for (k = r_n = 0; k < a_n; k++) {///RC_0/RC_1 to a[0, r_n) + if(a[k].del) continue; + if(k != r_n) { + t = a[k]; a[k] = a[r_n]; a[r_n] = t; + } + r_n++; + } + if(r_n <= 1) return occ; + qsort(a, r_n, sizeof((*a)), cmp_u_trans_weight); + + kv_resize(uint64_t, *srt, (b_n+(r_n<<2))); + for (k = 0; k < r_n; k++) { + srt->a[srt->n] = a[k].qs; + srt->a[srt->n] <<= 32; + srt->a[srt->n] |= k; + srt->n++; + } + + ss = srt->a + b_n; hs = ss + r_n; bu = hs + r_n; wu = bu + r_n; + radix_sort_arch64(ss, ss + r_n); + for (k = 0; k < r_n; k++) { + ss[k] = (uint32_t)ss[k]; hs[ss[k]] = k; + } + + for (k = 0; k < r_n; k++) { + cz = &(a[k]); len = (cz->qe-cz->qs)*sec_rate; m = wm = 0; + for (z = l = wl = 0; z < r_n; z++) { + if(z == hs[k]) { + assert(ss[z] == k); + continue; + } + if(ss[z] > k) continue;///smaller nw than a[k] + mz = &(a[ss[z]]); + if(mz->qs >= cz->qe) break; + if(mz->del) continue;///has been deleted + // if((mz->f != RC_0) && (mz->f != RC_1) && (cz->nw > (mz->nw*0.95))) continue; + os = MAX(cz->qs, mz->qs); oe = MIN(cz->qe, mz->qe); + ovq = ((oe > os)? (oe - os):0); + if(!ovq) continue; + if(is_connect_arc(cz, mz, ug, 0.04) || is_connect_arc(mz, cz, ug, 0.04)) continue; + + if((m > 0) && (((uint32_t)bu[m-1]) >= os)) { + if(oe > ((uint32_t)bu[m-1])) { + l += (oe - ((uint32_t)bu[m-1])); + bu[m-1] += (oe - ((uint32_t)bu[m-1])); + } + } else { + l += oe - os; + bu[m] = os; bu[m] <<= 32; bu[m] |= oe; m++; + } + + if((wm > 0) && (((uint32_t)wu[wm-1]) >= mz->qs)) { + if(mz->qe > ((uint32_t)wu[wm-1])) { + wl += (mz->qe - ((uint32_t)wu[wm-1])); + wu[wm-1] += (mz->qe - ((uint32_t)wu[wm-1])); + } + } else { + wl += mz->qe - mz->qs; + wu[wm] = mz->qs; wu[wm] <<= 32; wu[wm] |= mz->qe; wm++; + } + + if((l >= bd) && ((l>=len) || (l>=(wl*sec_rate)))) { + cz->del = 1; occ++; + ca = u_trans_a(*ta, cz->tn); c_n = u_trans_n(*ta, cz->tn); + for (i = 0; i < c_n; i++) { + if(ca[i].tn == cz->qn) ca[i].del = 1; + } + break; + } + } + } + + srt->n = b_n; + return occ; +} + 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); @@ -14617,6 +14712,646 @@ void refine_hic_trans(ug_opt_t *opt, kv_u_trans_t *ta, asg_t *sg, ma_ug_t *ug) // dbg_prt_utg_trans(ta, ug, "after"); } +void set_rset(ma_ug_t *ug, uint32_t uid, uint8_t *ff, uint8_t s) +{ + if(ug->g->seq[uid].del) return; + ma_utg_t* u; uint32_t k, v, nv, w, z; + asg_arc_t *av = NULL; + + u = &(ug->u.a[uid]); + for (k = 0; k < u->n; k++) ff[u->a[k]>>33] = s; + + v = uid<<1; + nv = asg_arc_n(ug->g, v); + av = asg_arc_a(ug->g, v); + for (k = 0; k < nv; k++) { + w = av[k].v; + if(av[k].del) continue; + if(ug->g->seq[w>>1].del) continue; + u = &(ug->u.a[w>>1]); + for (z = 0; z < u->n; z++) ff[u->a[z]>>33] = s; + } + + v = (uid<<1)+1; + nv = asg_arc_n(ug->g, v); + av = asg_arc_a(ug->g, v); + for (k = 0; k < nv; k++) { + w = av[k].v; + if(av[k].del) continue; + if(ug->g->seq[w>>1].del) continue; + u = &(ug->u.a[w>>1]); + for (z = 0; z < u->n; z++) ff[u->a[z]>>33] = s; + } +} + +uint32_t check_nc_status(asg_t *sg, ma_hit_t_alloc *src, uint8_t *ff, uint32_t id) +{ + uint32_t i, tn; ma_hit_t *h; + for (i = 0; i < src[id].length; i++) { + h = &(src[id].buffer[i]); tn = Get_tn((*h)); + if(sg->seq[tn].del) continue; + if(!ff[tn]) continue; + if((Get_qs((*h)) == 0) && (Get_qe((*h)) == sg->seq[id].len)) return 1; + } + return 0; +} + +uint64_t *gen_cov_rg(ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc *src) +{ + uint32_t i, k, rid, j, tn; uint8_t *ff; + uint64_t C_bases, *cc; ma_hit_t *h; ma_utg_t *u; + CALLOC(cc, sg->n_seq); CALLOC(ff, sg->n_seq); + + for (i = 0; i < ug->g->n_seq; i++) { + if(ug->g->seq[i].del) continue; + set_rset(ug, i, ff, 1); + + u = &(ug->u.a[i]); + for (k = 0; k < u->n; k++) { + C_bases = 0; rid = u->a[k]>>33; + for (j = 0; j < src[rid].length; j++) { + h = &(src[rid].buffer[j]); + tn = Get_tn((*h)); + if((!sg->seq[tn].del) && (!ff[tn])) continue; + if((sg->seq[tn].del) && (!check_nc_status(sg, src, ff, tn))) continue; + C_bases += (Get_qe((*h)) - Get_qs((*h))); + } + if(C_bases > cc[rid]) cc[rid] = C_bases; + } + + set_rset(ug, i, ff, 0); + } + free(ff); + + return cc; +} + +void fill_cov_arr(ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc *src, uint8_t *ff, uint64_t uid, uint64_t *res) +{ + ma_hit_t *h; ma_utg_t *u; + uint64_t C_bases, k, j, rid, tn; + set_rset(ug, uid, ff, 1); + + u = &(ug->u.a[uid]); + for (k = 0; k < u->n; k++) { + C_bases = 0; rid = u->a[k]>>33; + for (j = 0; j < src[rid].length; j++) { + h = &(src[rid].buffer[j]); + tn = Get_tn((*h)); + if((!sg->seq[tn].del) && (!ff[tn])) continue; + if((sg->seq[tn].del) && (!check_nc_status(sg, src, ff, tn))) continue; + C_bases += (Get_qe((*h)) - Get_qs((*h))); + } + res[k] = C_bases; + } + + set_rset(ug, uid, ff, 0); +} + +uint64_t infer_mmhap_copy(ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc *src, uint8_t *ff, uint64_t uid, uint64_t het_cov, uint64_t n_hap) +{ + ma_hit_t *h; ma_utg_t *u; + uint64_t C_bases, R_bases, cc, dd, k, j, rid, tn, md, mk; + set_rset(ug, uid, ff, 1); + + u = &(ug->u.a[uid]); + for (k = C_bases = R_bases = 0; k < u->n; k++) { + rid = u->a[k]>>33; R_bases += sg->seq[rid].len; + for (j = 0; j < src[rid].length; j++) { + h = &(src[rid].buffer[j]); + tn = Get_tn((*h)); + if((!sg->seq[tn].del) && (!ff[tn])) continue; + if((sg->seq[tn].del) && (!check_nc_status(sg, src, ff, tn))) continue; + C_bases += (Get_qe((*h)) - Get_qs((*h))); + } + // res[k] = C_bases; + } + + set_rset(ug, uid, ff, 0); + + if(R_bases) C_bases /= R_bases; + else C_bases = 0; + + md = mk = (uint64_t)-1; + for (k = 1; k <= n_hap; k++) { + cc = het_cov*k; + dd = ((C_bases>=cc)?(C_bases-cc):(cc-C_bases)); + if(dd <= md) { + md = dd; mk = k; + } + } + if(mk == ((uint64_t)-1)) mk = n_hap; + return mk; +} + +ug_rid_cov_t* gen_ug_rid_cov_t(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc *src) +{ + uint64_t k; uint8_t *ff; CALLOC(ff, rg->n_seq); + ug_rid_cov_t *p; CALLOC(p, 1); + p->rg = rg; p->ug = ug; + MALLOC(p->idx, ug->g->n_seq); + p->cov.n = p->cov.m = 0; p->cov.a = NULL; + for (k = 0; k < p->ug->u.n; k++) { + p->idx[k] = p->cov.n; p->cov.n += p->ug->u.a[k].n; + } + p->cov.m = p->cov.n; MALLOC(p->cov.a, p->cov.n); + for (k = 0; k < p->ug->u.n; k++) { + fill_cov_arr(ug, rg, src, ff, k, p->cov.a+p->idx[k]); + } + free(ff); + + uint64_t hom_cov, het_cov, hom_cut; + if(asm_opt.hom_global_coverage_set) { + hom_cov = asm_opt.hom_global_coverage; + } else { + hom_cov = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE); + } + het_cov = hom_cov/asm_opt.polyploidy; + hom_cut = hom_cov + (het_cov*(0.5+(((double)asm_opt.polyploidy)*0.05))); + p->hom_max = hom_cut; + p->hom_min = (het_cov*(asm_opt.polyploidy-1)) + (het_cov*(0.5+(((double)(asm_opt.polyploidy-1))*0.05))); + p->hom_cov = hom_cov; p->het_cov = het_cov; + return p; +} + +void destory_ug_rid_cov_t(ug_rid_cov_t *p) +{ + if(!p) return; + free(p->idx); free(p->cov.a); +} + +uint64_t get_ovlp_cov(ma_utg_t *in, uint64_t *cc, asg_t *sg, int64_t s, int64_t e, int64_t *itk, int64_t *itl, uint64_t reflen) +{ + int64_t k, l, n = in->n; uint64_t rid, o, ol; int64_t qs, qe, os, oe; + k = (*itk); l = (*itl); o = ol = 0; + if(k >= n) { + k = n-1; l = ((int64_t)(in->len))-((int64_t)(sg->seq[k].len)); + } + if(k < 0 || l < 0) { + k = 0; l = 0; + } + + + for (; k > 0; k--) { + rid = in->a[k]>>33; + qs = l; qe = l + sg->seq[rid].len; + if(qe <= s) break; + l -= (int64_t)((uint32_t)in->a[k]); + } + + for (; k < n; k++) { + rid = in->a[k]>>33; + qs = l; qe = l + sg->seq[rid].len; + if(qs >= e) break; + l += (uint32_t)in->a[k]; + if(qe <= s) continue; + + os = MAX(qs, s); oe = MIN(qe, e); + if(oe <= os) continue; + assert(oe > os); + o += (((double)(oe-os))/((double)(sg->seq[rid].len)))*(cc[rid]); + ol += oe - os; + } + + (*itk) = k; (*itl) = l; + o = ((ol)?(o/ol):(0)); + return o*reflen; +} + +uint64_t get_ovlp_cov_1(ma_utg_t *in, uint64_t *idx, asg_t *sg, int64_t s, int64_t e, int64_t *itk, int64_t *itl, uint64_t reflen) +{ + int64_t k, l, n = in->n; uint64_t rid, o, ol; int64_t qs, qe, os, oe; + k = (*itk); l = (*itl); o = ol = 0; + if(k >= n) { + k = n-1; l = ((int64_t)(in->len))-((int64_t)(sg->seq[k].len)); + } + if(k < 0 || l < 0) { + k = 0; l = 0; + } + + + for (; k > 0; k--) { + rid = in->a[k]>>33; + qs = l; qe = l + sg->seq[rid].len; + if(qe <= s) break; + l -= (int64_t)((uint32_t)in->a[k]); + } + + for (; k < n; k++) { + rid = in->a[k]>>33; + qs = l; qe = l + sg->seq[rid].len; + if(qs >= e) break; + l += (uint32_t)in->a[k]; + if(qe <= s) continue; + + os = MAX(qs, s); oe = MIN(qe, e); + if(oe <= os) continue; + assert(oe > os); + o += (((double)(oe-os))/((double)(sg->seq[rid].len)))*(idx[k]); + ol += oe - os; + } + + (*itk) = k; (*itl) = l; + o = ((ol)?(o/ol):(0)); + return o*reflen; +} + +uint32_t append_cov_line_ug_rid_cov_t(uint64_t uid, uint64_t *qcc, u_trans_t *p, ug_rid_cov_t *idx, uint64_t hom_cut, double cut_rate) +{ + uint32_t k, rid; uint64_t qs, qe, oqs, oqe, ovq, ots, ote, no, tot, ava; + ma_utg_t *qu, *tu; uint64_t s_shift, e_shift; int64_t itk, itl, l; + qu = &(idx->ug->u.a[p->qn]); tu = &(idx->ug->u.a[p->tn]); + if(qu->n <= 0 || tu->n <= 0) return 0; + if(p->rev) { + itk = tu->n-1; itl = ((int64_t)(tu->len))-((int64_t)(idx->rg->seq[itk].len)); + } else { + itk = 0; itl = 0; + } + // if(uid == 7929) { + // fprintf(stderr, "*****[M::%s]\tutg%.6lu%c\tqu->n::%u\ttu->n::%u\n", __func__, + // uid+1, "lc"[idx->ug->u.a[uid].circ], (uint32_t)qu->n, (uint32_t)tu->n); + // } + for (k = l = 0, tot = ava = 0; k < qu->n; k++) { + rid = qu->a[k]>>33; + qs = l; qe = l + idx->rg->seq[rid].len; + l += (uint32_t)qu->a[k]; + if(qe <= p->qs) continue; + if(qs >= p->qe) break; + ///qe > p->qs && qs < p->qe && qe > qs + + oqs = MAX(qs, p->qs); oqe = MIN(qe, p->qe); + ovq = ((oqe > oqs)? (oqe - oqs):0); + + // if(!(ovq > 0)) { + // fprintf(stderr, "[M::%s::] qn::%u, tn::%u, q::[%lu, %lu), p->q::[%u, %u)\n", + // __func__, p->qn, p->tn, qs, qe, p->qs, p->qe); + // } + if(ovq <= 0) continue; + assert(ovq > 0); + if((hom_cut != ((uint64_t)-1)) && (cut_rate >= 0) && (ovq <= (idx->rg->seq[rid].len*0.55))) { + continue; + } + + s_shift = get_offset_adjust(oqs-p->qs, p->qe-p->qs, p->te-p->ts); + e_shift = get_offset_adjust(p->qe-oqe, p->qe-p->qs, p->te-p->ts); + if(p->rev) { + ovq = s_shift; s_shift = e_shift; e_shift = ovq; + } + ots = p->ts + s_shift; ote = p->te-e_shift; + if(ote <= ots) continue; + no = get_ovlp_cov_1(tu, idx->cov.a + idx->idx[p->tn], idx->rg, ots, ote, &itk, &itl, idx->rg->seq[rid].len); + // if(uid == 7929) { + // fprintf(stderr, "[%u]\trid::%u\tno::%lu\tqcc[k]::%lu\trlen::%u\n", + // k, rid, no, qcc[k], idx->rg->seq[rid].len); + // } + // to = qcc[k] + no; + // qcc[k] = to; + if((hom_cut == ((uint64_t)-1)) || (cut_rate < 0)) { + qcc[k] += no; + } else { + tot++; + if((qcc[k] + no) <= (hom_cut*((uint64_t)idx->rg->seq[rid].len))) ava++; + } + } + + // if(uid == 7929) { + // fprintf(stderr, "[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\tis_reliable::%u\tw::%f\n", __func__, + // p->qn+1, "lc"[idx->ug->u.a[p->qn].circ], idx->ug->u.a[p->qn].len, p->qs, p->qe, "+-"[p->rev], + // p->tn+1, "lc"[idx->ug->u.a[p->tn].circ], idx->ug->u.a[p->tn].len, p->ts, p->te, ((p->f == RC_0) || (p->f == RC_1))?1:0, p->nw); + // fprintf(stderr, "[M::%s]\tutg%.6lu%c\thom_cut::%lu\tcut_rate::%f\tava::%lu\ttot::%lu\n", __func__, + // uid+1, "lc"[idx->ug->u.a[uid].circ], hom_cut, cut_rate, ava, tot); + // } + + if((hom_cut != ((uint64_t)-1)) && (cut_rate >= 0) && (ava >= (tot*cut_rate))) { + return 1; + } + return 0; +} + +uint32_t append_cov_line(uint64_t uid, uint64_t *qcc, uint64_t *cc, u_trans_t *p, ma_ug_t *ug, asg_t *sg, uint64_t hom_cut, double cut_rate) +{ + uint32_t k, rid; uint64_t qs, qe, oqs, oqe, ovq, ots, ote, no, tot, ava; + ma_utg_t *qu, *tu; uint64_t s_shift, e_shift; int64_t itk, itl, l; + qu = &(ug->u.a[p->qn]); tu = &(ug->u.a[p->tn]); + if(qu->n <= 0 || tu->n <= 0) return 0; + if(p->rev) { + itk = tu->n-1; itl = ((int64_t)(tu->len))-((int64_t)(sg->seq[itk].len)); + } else { + itk = 0; itl = 0; + } + // if(uid == 315 || uid == 1055) { + // fprintf(stderr, "*****[M::%s]\tutg%.6lu%c\tqu->n::%u\ttu->n::%u\n", __func__, + // uid+1, "lc"[ug->u.a[uid].circ], (uint32_t)qu->n, (uint32_t)tu->n); + // } + for (k = l = 0, tot = ava = 0; k < qu->n; k++) { + rid = qu->a[k]>>33; + qs = l; qe = l + sg->seq[rid].len; + l += (uint32_t)qu->a[k]; + if(qe <= p->qs) continue; + if(qs >= p->qe) break; + ///qe > p->qs && qs < p->qe && qe > qs + + oqs = MAX(qs, p->qs); oqe = MIN(qe, p->qe); + ovq = ((oqe > oqs)? (oqe - oqs):0); + + // if(!(ovq > 0)) { + // fprintf(stderr, "[M::%s::] qn::%u, tn::%u, q::[%lu, %lu), p->q::[%u, %u)\n", + // __func__, p->qn, p->tn, qs, qe, p->qs, p->qe); + // } + if(ovq <= 0) continue; + assert(ovq > 0); + if((hom_cut != ((uint64_t)-1)) && (cut_rate >= 0) && (ovq <= (sg->seq[rid].len*0.55))) { + continue; + } + + s_shift = get_offset_adjust(oqs-p->qs, p->qe-p->qs, p->te-p->ts); + e_shift = get_offset_adjust(p->qe-oqe, p->qe-p->qs, p->te-p->ts); + if(p->rev) { + ovq = s_shift; s_shift = e_shift; e_shift = ovq; + } + ots = p->ts + s_shift; ote = p->te-e_shift; + if(ote <= ots) continue; + no = get_ovlp_cov(tu, cc, sg, ots, ote, &itk, &itl, sg->seq[rid].len); + // if(uid == 315 || uid == 1055) { + // fprintf(stderr, "[%u]\trid::%u\tno::%lu\tqcc[k]::%lu\trlen::%u\n", k, rid, no, qcc[k], sg->seq[rid].len); + // } + // to = qcc[k] + no; + // qcc[k] = to; + if((hom_cut == ((uint64_t)-1)) || (cut_rate < 0)) { + qcc[k] += no; + } else { + tot++; + if((qcc[k] + no) <= (hom_cut*((uint64_t)sg->seq[rid].len))) ava++; + } + } + + // if(uid == 315 || uid == 1055) { + // fprintf(stderr, "[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\tis_reliable::%u\tw::%f\n", __func__, + // 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->f == RC_0) || (p->f == RC_1))?1:0, p->nw); + // fprintf(stderr, "[M::%s]\tutg%.6lu%c\thom_cut::%lu\tcut_rate::%f\tava::%lu\ttot::%lu\n", __func__, + // uid+1, "lc"[ug->u.a[uid].circ], hom_cut, cut_rate, ava, tot); + // } + + if((hom_cut != ((uint64_t)-1)) && (cut_rate >= 0) && (ava >= (tot*cut_rate))) { + return 1; + } + return 0; +} + +void purge_ovlp_cov(uint32_t id, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *sg, asg64_v *b64, uint64_t *cc, uint64_t hom_cut) +{ + ma_utg_t *u = &(ug->u.a[id]); uint32_t k, n, cn; + u_trans_t *a = NULL, t; + kv_resize(uint64_t, *b64, u->n); + for (k = 0; k < u->n; k++) b64->a[k] = cc[u->a[k]>>33]; + a = u_trans_a(*ta, id); n = u_trans_n(*ta, id); + // if(id == 315 || id == 1055) { + // fprintf(stderr, "\n[M::%s]\tutg%.6u%c\tn::%u\n", __func__, id+1, "lc"[ug->u.a[id].circ], n); + // } + for (k = cn = 0; k < n; k++) {///RC_0/RC_1 to a[0, r_n) + // if(id == 315 || id == 1055) { + // fprintf(stderr, "-0-[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\tis_reliable::%u\tw::%f\n", __func__, + // a[k].qn+1, "lc"[ug->u.a[a[k].qn].circ], ug->u.a[a[k].qn].len, a[k].qs, a[k].qe, "+-"[a[k].rev], + // a[k].tn+1, "lc"[ug->u.a[a[k].tn].circ], ug->u.a[a[k].tn].len, a[k].ts, a[k].te, + // ((a[k].f == RC_0) || (a[k].f == RC_1))?1:0, a[k].nw); + // } + if((a[k].f == RC_0) || (a[k].f == RC_1)) { + if(k != cn) { + t = a[k]; a[k] = a[cn]; a[cn] = t; + } + cn++; + } + } + // if(id == 315 || id == 1055) { + // fprintf(stderr, "[M::%s]\tutg%.6u%c\tn::%u\tcn::%u\n", __func__, + // id+1, "lc"[ug->u.a[id].circ], n, cn); + // } + if(cn >= n) return; + if(n-cn > 1) qsort(a+cn, n-cn, sizeof((*a)), cmp_u_trans_weight); + for (k = 0; k < cn; k++) {///calculate coverage for + append_cov_line(id, b64->a, cc, &(a[k]), ug, sg, ((uint64_t)-1), -1); + } + + for (k = cn; k < n; k++) { + if(append_cov_line(id, b64->a, cc, &(a[k]), ug, sg, hom_cut, 0.51)) { + append_cov_line(id, b64->a, cc, &(a[k]), ug, sg, ((uint64_t)-1), -1); + continue; + } + a[k].del = 1; + } +} + +void trans_sec_cut_filter_mmhap(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc* src) +{ + uint64_t k, hom_cov, het_cov, hom_cut; + asg64_v b64; kv_init(b64); + uint64_t *cc = gen_cov_rg(ug, sg, src); + if(asm_opt.hom_global_coverage_set) { + hom_cov = asm_opt.hom_global_coverage; + } else { + hom_cov = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE); + } + het_cov = hom_cov/asm_opt.polyploidy; + hom_cut = hom_cov + (het_cov*(0.5+(((double)asm_opt.polyploidy)*0.05))); + + fprintf(stderr, "+[M::%s]\thom_cov::%lu\thet_cov::%lu\thom_cut::%lu\n", __func__, + hom_cov, het_cov, hom_cut); + + for (k = 0; k < ug->g->n_seq; k++) { + purge_ovlp_cov(k, ta, ug, sg, &b64, cc, hom_cut); + } + free(cc); kv_destroy(b64); +} + + +void purge_ovlp_cov_adv(uint32_t id, kv_u_trans_t *ta, asg64_v *b64, ug_rid_cov_t *cc, uint64_t hom_cut) +{ + ma_utg_t *u = &(cc->ug->u.a[id]); uint32_t k, n, cn; + u_trans_t *a = NULL, t; + kv_resize(uint64_t, *b64, u->n); + memcpy(b64->a, cc->cov.a+cc->idx[id], sizeof((*(b64->a)))*cc->ug->u.a[id].n); + a = u_trans_a(*ta, id); n = u_trans_n(*ta, id); + // if(id == 7929) { + // fprintf(stderr, "\n[M::%s]\tutg%.6u%c\tn::%u\n", __func__, id+1, "lc"[cc->ug->u.a[id].circ], n); + // } + for (k = cn = 0; k < n; k++) {///RC_0/RC_1 to a[0, r_n) + // if(id == 7929) { + // fprintf(stderr, "-0-[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\tis_reliable::%u\tw::%f\n", __func__, + // a[k].qn+1, "lc"[cc->ug->u.a[a[k].qn].circ], cc->ug->u.a[a[k].qn].len, a[k].qs, a[k].qe, "+-"[a[k].rev], + // a[k].tn+1, "lc"[cc->ug->u.a[a[k].tn].circ], cc->ug->u.a[a[k].tn].len, a[k].ts, a[k].te, + // ((a[k].f == RC_0) || (a[k].f == RC_1))?1:0, a[k].nw); + // } + if((a[k].f == RC_0) || (a[k].f == RC_1)) { + if(k != cn) { + t = a[k]; a[k] = a[cn]; a[cn] = t; + } + cn++; + } + } + // if(id == 7929) { + // fprintf(stderr, "[M::%s]\tutg%.6u%c\tn::%u\tcn::%u\n", __func__, + // id+1, "lc"[cc->ug->u.a[id].circ], n, cn); + // } + if(cn >= n) return; + if(n-cn > 1) qsort(a+cn, n-cn, sizeof((*a)), cmp_u_trans_weight); + for (k = 0; k < cn; k++) {///calculate coverage for + append_cov_line_ug_rid_cov_t(id, b64->a, &(a[k]), cc, ((uint64_t)-1), -1); + } + + for (k = cn; k < n; k++) { + if(append_cov_line_ug_rid_cov_t(id, b64->a, &(a[k]), cc, hom_cut, 0.51)) { + append_cov_line_ug_rid_cov_t(id, b64->a, &(a[k]), cc, ((uint64_t)-1), -1); + continue; + } + a[k].del = 1; + } +} + +void trans_sec_cut_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc* src) +{ + uint64_t k; + asg64_v b64; kv_init(b64); + ug_rid_cov_t *cc = gen_ug_rid_cov_t(ug, sg, src); + + fprintf(stderr, "+[M::%s]\thom_cov::%lu\thet_cov::%lu\thom_cut::%lu\n", __func__, + cc->hom_cov, cc->het_cov, cc->hom_max); + + for (k = 0; k < ug->g->n_seq; k++) { + purge_ovlp_cov_adv(k, ta, &b64, cc, cc->hom_max); + } + + destory_ug_rid_cov_t(cc); free(cc); kv_destroy(b64); +} + +static void worker_for_trans_sec_simple_cut(void *data, long i, int tid) // callback for kt_for() +{ + u_trans_clean_t *s = (u_trans_clean_t*)data; + kv_u_trans_t *ta = s->ta; u_trans_t *a, *za; + uint32_t id = i, n, k, z, zn; + + a = u_trans_a(*ta, id); n = u_trans_n(*ta, id); + for (k = 0; k < n; k++) { + if((!a[k].del)) continue; + za = u_trans_a(*ta, a[k].tn); zn = u_trans_n(*ta, a[k].tn); + for (z = 0; z < zn; z++) { + if(za[z].tn == id) break; + } + assert(z < zn); + if(za[z].del) { ///both a[k] and za[z] have been deleted + if(a[k].qn < a[k].tn) a[k].occ = za[z].occ = ((uint32_t)-1); + continue; + } + ///a[k].del == 1 && za[z].del == 0 + if(za[z].qe-za[z].qs >= 3000000) { + a[k].occ = za[z].occ = 0; + } else { + a[k].occ = za[z].occ = ((uint32_t)-1); + } + } +} + + +void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* src) +{ + u_trans_clean_t sl; uint64_t k, i, l, st, occ; ha_mzl_t *tz; + ha_mzl_v srt_a; kv_u_trans_t *bl; u_trans_t *z; + memset(&sl, 0, sizeof(sl)); kv_init(srt_a); + + kt_u_trans_t_idx(ta, ug->g->n_seq); + sl.n_thread = asm_opt.thread_num; sl.ug = ug; sl.rg = read_g; sl.ta = ta; + CALLOC(sl.res, sl.n_thread); + sl.ov_cutoff = 3; sl.small_ov_rate = 0.8; + // dbg_prt_utg_trans(ta, ug, "pre"); + kt_for(sl.n_thread, worker_for_trans_clean_re, &sl, sl.ta->idx.n); + + for (i = srt_a.n = occ = 0; i < sl.n_thread; i++) { + bl = &(sl.res[i]); + if(!(bl->n)) continue; + for (k = 1, l = 0; k <= bl->n; k++) { + if(k == bl->n || bl->a[k].qn != bl->a[l].qn) { + if(k > l) { + kv_pushp(ha_mzl_t, srt_a, &tz); + tz->x = bl->a[l].qn; tz->x <<= 32; tz->x |= i; + tz->rid = l>>32; tz->pos = (uint32_t)l; + occ += (k - l); + } + l = k; + } + } + } + // fprintf(stderr, "[M::%s::] occ::%lu \n", __func__, occ); + assert(srt_a.n <= sl.ug->u.n); + radix_sort_ha_mzl_t_srt1(srt_a.a, srt_a.a + srt_a.n); + occ <<= 1; ta->n = 0; kv_resize(u_trans_t, *ta, occ); + for (i = 0; i < srt_a.n; i++) { + tz = &(srt_a.a[i]); + bl = &(sl.res[(uint32_t)(tz->x)]); + k = tz->rid; k <<= 32; k += tz->pos; + assert(bl->a[k].qn == (tz->x>>32)); + for (; (k < bl->n) && (bl->a[k].qn == (tz->x>>32)); k++) { + if(bl->a[k].qn == bl->a[k].tn) continue; + if(bl->a[k].qe-bl->a[k].qs <= 0) continue; + if(bl->a[k].te-bl->a[k].ts <= 0) continue; + kv_pushp(u_trans_t, *ta, &z); + (*z) = bl->a[k]; z->occ = 0; if(z->f == RC_3) z->f = RC_2; + // if(z->ts >= z->te || z->qs >= z->qe) { + // fprintf(stderr, "+[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\n", __func__, + // z->qn+1, "lc"[ug->u.a[z->qn].circ], ug->u.a[z->qn].len, z->qs, z->qe, "+-"[z->rev], + // z->tn+1, "lc"[ug->u.a[z->tn].circ], ug->u.a[z->tn].len, z->ts, z->te); + // } + + kv_pushp(u_trans_t, *ta, &z); + (*z) = bl->a[k]; z->occ = 0; if(z->f == RC_3) z->f = RC_2; + z->qn = bl->a[k].tn; z->qs = bl->a[k].ts; z->qe = bl->a[k].te; + z->tn = bl->a[k].qn; z->ts = bl->a[k].qs; z->te = bl->a[k].qe; + + // if(z->ts >= z->te || z->qs >= z->qe) { + // fprintf(stderr, "-[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\n", __func__, + // z->qn+1, "lc"[ug->u.a[z->qn].circ], ug->u.a[z->qn].len, z->qs, z->qe, "+-"[z->rev], + // z->tn+1, "lc"[ug->u.a[z->tn].circ], ug->u.a[z->tn].len, z->ts, z->te); + // } + } + } + + for (k = 0; k < sl.n_thread; k++) free(sl.res[k].a); + free(sl.res); kv_destroy(srt_a); + kt_u_trans_t_idx(ta, ug->g->n_seq); + // dbg_prt_utg_trans(ta, ug, "after"); + trans_sec_cut_filter_mmhap_adv(ta, ug, read_g, src); + kt_for(sl.n_thread, worker_for_trans_sec_simple_cut, &sl, sl.ta->idx.n); + + // CALLOC(sl.srt, sl.n_thread); sl.sec_rate = 0.5; + // kt_for(sl.n_thread, worker_for_trans_sec_cut, &sl, sl.ta->idx.n); + // sl.cu = gen_u_trans_cluster(ta); + // kt_for(sl.n_thread, worker_for_trans_sec_cut, &sl, sl.cu->idx.n); + // free(sl.cu->idx.a); free(sl.cu->z.a); free(sl.cu); + // for (k = 0; k < sl.n_thread; k++) free(sl.srt[k].a); free(sl.srt); + for (k = st = 0; k < ta->n; k++) { + if(ta->a[k].occ == ((uint32_t)-1)) continue; + ta->a[st] = ta->a[k]; ta->a[st].del = 0; st++; + } + ta->n = st; + for (st = 0, i = 1; i <= ta->n; ++i) { + if (i == ta->n || ta->a[i].qn != ta->a[st].qn) { + ta->idx.a[ta->a[st].qn] = (((uint64_t)st)<<32)|(i-st); + st = i; + } + } + // dbg_prt_utg_extra_trans(ta, ug, asm_opt.output_file_name); +} + +void refine_hic_trans_mmhap(ug_opt_t *opt, kv_u_trans_t *ta, asg_t *sg, ma_ug_t *ug) +{ + filter_u_trans(ta, asm_opt.is_bub_trans, asm_opt.is_topo_trans, asm_opt.is_read_trans, asm_opt.is_base_trans); + if(asm_opt.is_base_trans) { + trans_base_mmhap_infer(ug, sg, opt, ta); + } + // dbg_prt_utg_trans(ta, ug, "pre"); + // clean_u_trans_t_idx_filter_mmhap_adv(ta, ug, sg, opt->sources); + // dbg_prt_utg_trans(ta, ug, "after"); +} + void output_contig_graph_alternative(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp); void output_hic_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, @@ -14691,7 +15426,7 @@ long long gap_fuzz, bub_label_t* b_mask_t, ug_opt_t *opt) - hic_analysis(ug, sg, cov?cov->t_ch:t_ch, opt, 0, asm_opt.scffold?&rhits:NULL); + hic_analysis(ug, sg, cov?cov->t_ch:t_ch, opt, NULL, asm_opt.scffold?&rhits:NULL); if(!rhits && cov) destory_hap_cov_t(&cov); @@ -14747,6 +15482,256 @@ long long gap_fuzz, bub_label_t* b_mask_t, ug_opt_t *opt) } } + +uint64_t dump_trans_ovlp(trans_chain* t_ch, const char *fn, ma_ug_t *ug) +{ + char* gfa_name = NULL; MALLOC(gfa_name, strlen(fn)+50); + sprintf(gfa_name, "%s.trans.bin", fn); + FILE* fp = fopen(gfa_name, "w"); free(gfa_name); + if (!fp) return 0; + if(ug) write_dbug(ug, fp); + + fwrite(&t_ch->r_num, sizeof(t_ch->r_num), 1, fp); + fwrite(t_ch->ir_het, sizeof(uint8_t), t_ch->r_num, fp); + uint64_t i; + fwrite(&t_ch->bed.n, sizeof(t_ch->bed.n), 1, fp); + for (i = 0; i < t_ch->bed.n; i++) { + fwrite(&t_ch->bed.a[i].n, sizeof(t_ch->bed.a[i].n), 1, fp); + fwrite(t_ch->bed.a[i].a, sizeof(bed_interval), t_ch->bed.a[i].n, fp); + } + + fwrite(&t_ch->k_trans.n, sizeof(t_ch->k_trans.n), 1, fp); + fwrite(t_ch->k_trans.a, sizeof(u_trans_t), t_ch->k_trans.n, fp); + + fwrite(&t_ch->k_trans.idx.n, sizeof(t_ch->k_trans.idx.n), 1, fp); + fwrite(t_ch->k_trans.idx.a, sizeof(uint64_t), t_ch->k_trans.idx.n, fp); + + fclose(fp); + fprintf(stderr, "[M::%s] Dump trans overlaps\n", __func__); + return 1; +} + +trans_chain* load_trans_ovlp(const char *fn, ma_ug_t *ug) +{ + char* gfa_name = NULL; MALLOC(gfa_name, strlen(fn)+50); + sprintf(gfa_name, "%s.trans.bin", fn); + FILE* fp = fopen(gfa_name, "r"); free(gfa_name); + if (!fp) return NULL; + if(ug && (!test_dbug(ug, fp))) { + fprintf(stderr, "[M::%s] Renew trans overlaps\n", __func__); + fclose(fp); + return NULL; + } + + uint64_t flag = 0; + trans_chain *t_ch = NULL; + CALLOC(t_ch, 1); + + flag += fread(&t_ch->r_num, sizeof(t_ch->r_num), 1, fp); + MALLOC(t_ch->ir_het, t_ch->r_num); + flag += fread(t_ch->ir_het, sizeof(uint8_t), t_ch->r_num, fp); + + uint64_t i; + flag += fread(&t_ch->bed.n, sizeof(t_ch->bed.n), 1, fp); + MALLOC(t_ch->bed.a, t_ch->bed.n); t_ch->bed.m = t_ch->bed.n; + for (i = 0; i < t_ch->bed.n; i++) { + flag += fread(&t_ch->bed.a[i].n, sizeof(t_ch->bed.a[i].n), 1, fp); + MALLOC(t_ch->bed.a[i].a, t_ch->bed.a[i].n); t_ch->bed.a[i].m = t_ch->bed.a[i].n; + flag += fread(t_ch->bed.a[i].a, sizeof(bed_interval), t_ch->bed.a[i].n, fp); + } + + flag += fread(&t_ch->k_trans.n, sizeof(t_ch->k_trans.n), 1, fp); + MALLOC(t_ch->k_trans.a, t_ch->k_trans.n); t_ch->k_trans.m = t_ch->k_trans.n; + flag += fread(t_ch->k_trans.a, sizeof(u_trans_t), t_ch->k_trans.n, fp); + + flag += fread(&t_ch->k_trans.idx.n, sizeof(t_ch->k_trans.idx.n), 1, fp); + MALLOC(t_ch->k_trans.idx.a, t_ch->k_trans.idx.n); t_ch->k_trans.idx.m = t_ch->k_trans.idx.n; + flag += fread(t_ch->k_trans.idx.a, sizeof(uint64_t), t_ch->k_trans.idx.n, fp); + + fclose(fp); + fprintf(stderr, "[M::%s::] ==> Trans overlaps have been loaded\n", __func__); + return t_ch; +} + +void update_trio_mmhap(uint32_t hapid, ma_ug_t *mm_ug, mmhap_t *rh, asg_t *sg, uint32_t n_hap) +{ + uint64_t k, z, m, rid; ma_utg_t *u; + for (k = 0; k < sg->n_seq; k++) { + if(R_INF.trio_flag[k] == DROP) continue; + R_INF.trio_flag[k] = AMBIGU; + } + + for (k = 0; k < mm_ug->u.n; k++) { + if(rh->h.a[k].m == n_hap || rh->h.a[k].n == 0) { + m = AMBIGU; + } else { + for (z = 0, m = MOTHER; z < rh->h.a[k].n; z++) { + if(rh->a.a[rh->h.a[k].a+z] == hapid) {m = FATHER; break;} + } + } + + u = &(mm_ug->u.a[k]); + for (z = 0; z < u->n; z++) { + rid = u->a[z]>>33; + if(R_INF.trio_flag[rid] == DROP) continue; + R_INF.trio_flag[rid] = m; + } + fprintf(stderr, "utg%.6lu%c\tm::%lu\n", k+1, "lc"[mm_ug->u.a[k].circ], m); + } +} + +void dbg_prt_trio_mmhap_label(ma_ug_t *ug, mmhap_t *rh, const char *o_n) +{ + char* gfa_name = (char*)malloc(strlen(o_n)+100); + FILE *fn = NULL; sprintf(gfa_name, "%s.mmhap.binning.log", o_n); + uint32_t i, z; fn = fopen(gfa_name, "w"); + for (i = 0; i < ug->g->n_seq; i++) { + fprintf(fn, "utg%.6u%c\tm::%u\tn::%u", i+1, "lc"[ug->u.a[i].circ], rh->h.a[i].m, rh->h.a[i].n); + for (z = 0; z < rh->h.a[i].n; z++) { + fprintf(fn, "\t%u", rh->a.a[rh->h.a[i].a+z]); + } + fprintf(fn, "\n"); + } + fclose(fn); free(gfa_name); fprintf(stderr, "[M::%s::] done\n", __func__); +} + +void output_trio_mmhap(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, long long gap_fuzz, bub_label_t* b_mask_t, ug_opt_t *opt, +ma_ug_t *mm_ug, mmhap_t *rh, uint32_t n_hap) +{ + uint32_t i; char *fp = NULL; MALLOC(fp, 100); + memset(R_INF.trio_flag, 0, sizeof((*(R_INF.trio_flag)))*sg->n_seq); + for (i = 0; i < n_hap; i++){ + sprintf(fp, "hap%u", i+1); + update_trio_mmhap(i, mm_ug, rh, sg, n_hap); + output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, reverse_sources, tipsLen, tip_drop_ratio, + stops_threshold, ruIndex, chimeric_rate, drop_ratio, max_hang, min_ovlp, 0, b_mask_t, fp, NULL, NULL); + } + free(fp); +} + +void output_hic_graph_mmhap(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, long long gap_fuzz, bub_label_t* b_mask_t, ug_opt_t *opt) +{ + hic_clean(sg); + mmhap_t *rh = NULL; + + kvec_asg_arc_t_warp new_rtg_edges, d_edges; + kv_init(new_rtg_edges.a); kv_init(d_edges.a); + ma_ug_t *ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); + + new_rtg_edges.a.n = 0; + ma_ug_seq(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, &d_edges, 1);///polish + + hap_cov_t *cov = NULL; + trans_chain* t_ch = NULL; + t_ch = load_trans_ovlp(output_file_name, ug); + // if((asm_opt.flag & HA_F_VERBOSE_GFA)) t_ch = load_hc_trans(output_file_name); + + if(!t_ch) { + new_rtg_edges.a.n = 0; + asg_t *copy_sg = copy_read_graph(sg); + ma_ug_t *copy_ug = copy_untig_graph(ug); + ///asm_opt.purge_overlap_len = asm_opt.purge_overlap_len_hic; + ///asm_opt.purge_simi_thres = asm_opt.purge_simi_rate_hic; + adjust_utg_by_primary(©_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, + tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, + max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 0); + print_utg(copy_ug, copy_sg, coverage_cut, output_file_name, sources, ruIndex, max_hang, + min_ovlp, &new_rtg_edges); + + if(asm_opt.is_alt) { + output_contig_graph_alternative(copy_sg, coverage_cut, output_file_name, sources, ruIndex, max_hang, + min_ovlp); + } + ma_ug_destroy(copy_ug); + asg_destroy(copy_sg); + // clean_u_trans_t_idx(&(cov->t_ch->k_trans), ug, sg); + + + new_rtg_edges.a.n = 0; + ma_ug_print_bed(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, + max_hang, min_ovlp, asm_opt.hic_inconsist_rate, NULL, NULL, cov->t_ch); + + refine_hic_trans_mmhap(opt, &(cov->t_ch->k_trans), sg, ug); + dump_trans_ovlp(cov->t_ch, output_file_name, ug); + // if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_trans_chain(cov->t_ch, output_file_name); + } + + dbg_prt_utg_trans(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, "pre"); + clean_u_trans_t_idx_filter_mmhap_adv(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, sg, opt->sources); + dbg_prt_utg_trans(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, "after"); + + // refine_hic_trans_mmhap(opt, &(cov?cov->t_ch->k_trans:t_ch->k_trans), sg, ug); + ///for debug + // dbg_prt_utg_trans(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, "dd"); + // char* gfa_name = (char*)malloc(strlen(output_file_name)+50); FILE* output_file = NULL; + + // sprintf(gfa_name, "%s.pre.clean_d_utg.noseq.gfa", output_file_name); + // output_file = fopen(gfa_name, "w"); + // ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file); + // fclose(output_file); + + // sprintf(gfa_name, "%s.pre.clean_d_utg.gfa", output_file_name); + // output_file = fopen(gfa_name, "w"); + // ma_ug_print(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, &rh, NULL); + + + if(cov) destory_hap_cov_t(&cov); + if(t_ch) destory_trans_chain(&t_ch); + + + // char* gfa_name = (char*)malloc(strlen(output_file_name)+50); + // sprintf(gfa_name, "%s.after.clean_d_utg.noseq.gfa", output_file_name); + // FILE* output_file = fopen(gfa_name, "w"); + // ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file); + // fclose(output_file); + // free(gfa_name); + + // ma_ug_destroy(ug); + kv_destroy(new_rtg_edges.a); + + + asg_arc_t* av = NULL; + uint32_t v, w, k, i, nv; + for (i = 0; i < d_edges.a.n; i++) + { + v = d_edges.a.a[i].ul>>32; + w = d_edges.a.a[i].v; + av = asg_arc_a(sg, v); + nv = asg_arc_n(sg, v); + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + if(av[k].v == w) + { + av[k].del = 1; + break; + } + } + } + kv_destroy(d_edges.a); + asg_cleanup(sg); + + dbg_prt_trio_mmhap_label(ug, rh, output_file_name); + + output_trio_mmhap(sg, coverage_cut, output_file_name, sources, reverse_sources, tipsLen, tip_drop_ratio, + stops_threshold, ruIndex, chimeric_rate, drop_ratio, max_hang, min_ovlp, gap_fuzz, b_mask_t, opt, ug, rh, asm_opt.polyploidy); + + ma_ug_destroy(ug); + free(rh->a.a); free(rh->h.a); free(rh); +} + void print_debug_gfa(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, ma_hit_t_alloc* sources, R_to_U* ruIndex) { @@ -14851,7 +15836,7 @@ long long gap_fuzz, bub_label_t* b_mask_t) if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_trans_chain(t_ch, output_file_name); } - hic_analysis(ug, sg, t_ch, &opt, 0, NULL); + hic_analysis(ug, sg, t_ch, &opt, NULL, NULL); destory_trans_chain(&t_ch); ma_ug_destroy(ug); asg_cleanup(sg); @@ -34056,8 +35041,14 @@ ma_sub_t **coverage_cut_ptr, int debug_g) else if(ha_opt_hic(&asm_opt)) { if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION; - output_hic_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, gap_fuzz, &b_mask_t, &uopt); + if(asm_opt.polyploidy <= 2) { + output_hic_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, gap_fuzz, &b_mask_t, &uopt); + } else { + output_hic_graph_mmhap(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, gap_fuzz, &b_mask_t, &uopt); + } + // output_hic_graph_polyploid(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, gap_fuzz, &b_mask_t); } diff --git a/Overlaps.h b/Overlaps.h index 5d0107f..65d116e 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1211,4 +1211,27 @@ int asg_arc_del_trans_ul(asg_t *g, int fuzz); #define JUNK_COV 5 #define DISCARD_RATE 0.8 +typedef struct { + uint32_t n, m, a; +} mmhap_status_t; + +typedef struct { + kvec_t(mmhap_status_t) h; + kvec_t(uint32_t) a; +} mmhap_t; + +typedef struct { // global data structure for kt_pipeline() + ma_ug_t *ug; + asg_t *rg; + uint64_t *idx; + asg64_v cov; + uint64_t hom_min, hom_max, hom_cov, het_cov; +} ug_rid_cov_t; + +ug_rid_cov_t* gen_ug_rid_cov_t(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc *src); +void destory_ug_rid_cov_t(ug_rid_cov_t *p); +uint32_t append_cov_line_ug_rid_cov_t(uint64_t uid, uint64_t *qcc, u_trans_t *p, ug_rid_cov_t *idx, uint64_t hom_cut, double cut_rate); +uint64_t infer_mmhap_copy(ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc *src, uint8_t *ff, uint64_t uid, uint64_t het_cov, uint64_t n_hap); +uint64_t trans_sec_cut0(kv_u_trans_t *ta, asg64_v *srt, uint32_t id, double sec_rate, uint64_t bd, ma_ug_t *ug); + #endif diff --git a/anchor.cpp b/anchor.cpp index 2e79b39..24b553a 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -1620,7 +1620,7 @@ inline uint64_t special_lchain(Candidates_list* cl, overlap_region_alloc* ol, ui double chn_pen_skip, double bw_rate, int64_t quick_check, uint32_t gen_off, double mcopy_rate, uint32_t mcopy_khit_cut, st_mt_t *sp, uint64_t *si, uint64_t m, uint64_t l, uint64_t k) { - uint64_t z, s, e, yid, ol0 = ol->length, ol1; + uint64_t z, s, e, yid, ol0 = ol->length, ol1, m0 = m; if(l >= k) return m; yid = cl->list[l].readID; for (z = l; z < k && cl->list[z].strand == cl->list[l].strand; z++); @@ -1678,7 +1678,21 @@ inline uint64_t special_lchain(Candidates_list* cl, overlap_region_alloc* ol, ui } s++; } - ol->length = s; + + if(s != ol->length) { + ol->length = s; + for (z = ol0; z < ol->length; z++) { + s = ol->list[z].non_homopolymer_errors; + pi = cl->list[s].readID; + ol->list[z].non_homopolymer_errors = m0; + for (; s < m && cl->list[s].readID == pi; s++) { + cl->list[m0] = cl->list[s]; + cl->list[m0].readID = z; m0++; + } + } + m = m0; + } + } return m; } diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 01defe8..89e2860 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -1262,7 +1262,7 @@ uint32_t minLen, asg64_v *t) } void asg_arc_cut_length(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio, -uint32_t is_topo, uint32_t min_diff, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len) +uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len) { asg64_v tx = {0,0,0}, *b = NULL; uint32_t i, k, v, w, n_vtx = g->n_seq<<1, nv, nw, kv, kw, trioF = (uint32_t)-1, ntrioF = (uint32_t)-1, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou; @@ -1341,7 +1341,7 @@ uint32_t is_topo, uint32_t min_diff, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t * if (kv < 1) continue; if (kv >= 2) { if (mm_ol > ol_max*len_rat) continue; - if (is_ou && mm_ou > ou_max*ou_rat) continue; + if (is_ou && mm_ou > ou_max*ou_rat && mm_ou > min_ou) continue; if ((mm_ol + min_diff) > ol_max) continue; } @@ -1356,7 +1356,7 @@ uint32_t is_topo, uint32_t min_diff, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t * if (kw < 1) continue; if (kw >= 2) { if (mm_ol > ol_max*len_rat) continue; - if (is_ou && mm_ou > ou_max*ou_rat) continue; + if (is_ou && mm_ou > ou_max*ou_rat && mm_ou > min_ou) continue; if ((mm_ol + min_diff) > ol_max) continue; } @@ -2798,7 +2798,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i // prt_specfic_sge(sg, 10531, 10519, "--2--"); asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); - asg_arc_cut_length(sg, &bu, max_tip, drop, ou_drop_rate, is_ou, is_trio, 1, min_diff, NULL, NULL, NULL); + asg_arc_cut_length(sg, &bu, max_tip, drop, ou_drop_rate, is_ou, is_trio, 1, min_diff, 1, NULL, NULL, NULL); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); // prt_specfic_sge(sg, 10531, 10519, "--3--"); @@ -2851,13 +2851,18 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i if(!is_ou) { ///asg_arc_del_triangular_directly might be unnecessary asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); - asg_arc_cut_length(sg, &bu, max_tip, HARD_ORTHOLOGY_DROP/**min_ovlp_drop_ratio**/, ou_drop_rate, is_ou, 0/**is_trio**/, is_ou?1:0, min_diff, rev, rI, NULL); + asg_arc_cut_length(sg, &bu, max_tip, HARD_ORTHOLOGY_DROP/**min_ovlp_drop_ratio**/, ou_drop_rate, is_ou, 0/**is_trio**/, is_ou?1:0, min_diff, 1, rev, rI, NULL); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); // print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty7.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0); asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); - asg_arc_cut_length(sg, &bu, max_tip, min_ovlp_drop_ratio, ou_drop_rate, is_ou, 0/**is_trio**/, is_ou?1:0, min_diff, rev, rI, &l_drop); + asg_arc_cut_length(sg, &bu, max_tip, min_ovlp_drop_ratio, ou_drop_rate, is_ou, 0/**is_trio**/, is_ou?1:0, min_diff, 1, rev, rI, &l_drop); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); + } else { + min_diff = step_diff; l_drop = 6000; + asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); + asg_arc_cut_length(sg, &bu, max_tip, 0.3, 0.9, is_ou, is_trio, 1, min_diff, 8, NULL, NULL, &l_drop); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); } @@ -11923,16 +11928,18 @@ uint32_t ava_pass_unique_bridge(uint64_t *idx, uint64_t *integer_seq, uint64_t s return 1; } -uint32_t ava_pass_unique_bridge_tips(usg_t *g, asg64_v *b64, uint64_t g_s, uint64_t g_e, uint64_t *integer_seq, uint8_t *f, uint64_t max_ext) +uint32_t ava_pass_unique_bridge_tips(usg_t *g, asg64_v *b64, uint64_t g_s, uint64_t g_e, uint64_t *integer_seq, uint8_t *f, uint64_t max_ext, double max_ext_rate, uint64_t ext_up) { - uint64_t i, k, z, s, e, v, nv, bn = b64->n, kv, n_ext = 0; usg_arc_t *av; + uint64_t i, k, z, s, e, v, nv, bn = b64->n, kv, n_ext = 0, nkeep = 0, ncut = 0; usg_arc_t *av; for (i = g_s; i < g_e; i++) {///available intervals within the same cluster s = b64->a[((uint32_t)b64->a[i])]>>32; e = ((uint32_t)b64->a[((uint32_t)b64->a[i])]); assert(s < e); for (k = s + 1; k < e; k++) {///note: here is [s, e] - f[integer_seq[k]] = f[integer_seq[k]^1] = 1; + f[integer_seq[k]] = f[integer_seq[k]^1] = 1; + nkeep += g->a[integer_seq[k]>>1].occ; } f[integer_seq[s]] = f[integer_seq[e]^1] = 1; } + nkeep = nkeep*max_ext_rate; ///collect nodes within raw unitig graph that are linked by the clusters but not in the cluster for (i = g_s; i < g_e; i++) { @@ -11964,7 +11971,10 @@ uint32_t ava_pass_unique_bridge_tips(usg_t *g, asg64_v *b64, uint64_t g_s, uint6 } } - while (b64->n > bn && n_ext < max_ext) { + ncut = MAX(nkeep, max_ext); + if(ncut > ext_up) ncut = ext_up; + if(ncut < max_ext) ncut = max_ext; + while (b64->n > bn && n_ext < ncut) { v = b64->a[--b64->n]; if(f[v]) continue; av = usg_arc_a(g, v^1); nv = usg_arc_n(g, v^1); for (i = kv = 0; i < nv && kv < 1; i++) { @@ -12014,7 +12024,7 @@ uint32_t ava_pass_unique_bridge_tips(usg_t *g, asg64_v *b64, uint64_t g_s, uint6 } } } - + if((n_ext < ncut) && (n_ext > max_ext - 1)) n_ext = max_ext - 1; return n_ext; } @@ -12954,7 +12964,7 @@ uint8_t *ff, uint32_t *ng_occ, uint32_t max_ext) n_clus = 0; for (k = a_n + 1, i = a_n, mm = a_n; k <= b64->n; k++) { if(k == b64->n || (b64->a[k]>>32) != (b64->a[i]>>32)) { - tip_l = ava_pass_unique_bridge_tips(ng, b64, i, k, int_a, ff, max_ext); + tip_l = ava_pass_unique_bridge_tips(ng, b64, i, k, int_a, ff, max_ext, 0.03, 16); if(tip_l < max_ext) {///no long tip if(/**(!tip_l) || (**/ava_pass_unique_bridge_cov(uidx, ng, b64, i, k, int_a, tip_l)) { for (z = i; z < k; z++) { @@ -15934,11 +15944,14 @@ void u2g_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint32_t keep_raw_utg, uint // char sb[1000]; // output_integer_graph(uidx, iug, "ig_h0", 0); - // fprintf(stderr, "\n[M::%s::] max_path_drop_ratio::%f, min_path_drop_ratio::%f\n", __func__, ulopt->max_path_drop_ratio, ulopt->min_path_drop_ratio); + // fprintf(stderr, "\n[M::%s::] max_path_drop_ratio::%f, min_path_drop_ratio::%f, max_tip_hifi::%ld, max_tip::%ld\n", + // __func__, ulopt->max_path_drop_ratio, ulopt->min_path_drop_ratio, + // ulopt->max_tip_hifi, ulopt->max_tip); for (ss = 0; ss < 2; ss++) { for (i = 0, drop = ulopt->min_path_drop_ratio; i < ulopt->clean_round; i++, drop += step) { if(drop > ulopt->max_path_drop_ratio) drop = ulopt->max_path_drop_ratio; - // fprintf(stderr, "[M::%s::] Starting round-%ld, drop::%f\n", __func__, i, drop); + // fprintf(stderr, "[M::%s::] Starting round-%ld, drop::%f, max_tip_hifi::%ld\n", + // __func__, i, drop, ulopt->max_tip_hifi); cnt = 1; topo_level = 2; mm_tip = ulopt->max_tip; while (cnt) { cnt = 0; diff --git a/gfa_ut.h b/gfa_ut.h index 32c5675..c39e0d8 100644 --- a/gfa_ut.h +++ b/gfa_ut.h @@ -22,7 +22,7 @@ void asg_iterative_semi_circ(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_ void asg_arc_cut_chimeric(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t ou_thres); void asg_arc_cut_inexact(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max_ext, uint32_t is_ou, uint32_t is_trio, uint32_t min_diff, float ou_rat/**, asg64_v *dbg**/); void asg_arc_cut_length(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio, -uint32_t is_topo, uint32_t min_diff, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len); +uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len); void asg_arc_cut_bub_links(asg_t *g, asg64_v *in, float len_rat, float sec_len_rat, float ou_rat, uint32_t is_ou, uint64_t check_dist, ma_hit_t_alloc *rev, R_to_U* rI, int32_t max_ext); void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float ou_rat, uint32_t is_ou, bub_label_t *b_mask_t); uint32_t asg_cut_large_indel(asg_t *g, asg64_v *in, int32_t max_ext, float ou_rat, uint32_t is_ou, uint32_t min_diff); diff --git a/hic.cpp b/hic.cpp index 59ec0c1..ac666cd 100644 --- a/hic.cpp +++ b/hic.cpp @@ -17430,8 +17430,280 @@ int hic_short_align_poy(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, return 1; } +mmhap_t* gen_mmhap_t(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc *src) +{ + uint64_t k, z, hom_cov, het_cov, s, *bs = NULL; uint8_t *ff; mmhap_t *p; + if(asm_opt.hom_global_coverage_set) { + hom_cov = asm_opt.hom_global_coverage; + } else { + hom_cov = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE); + } + het_cov = hom_cov/asm_opt.polyploidy; + CALLOC(ff, rg->n_seq); CALLOC(p, 1); CALLOC(bs, asm_opt.polyploidy+1); + p->h.n = p->h.m = ug->u.n; CALLOC(p->h.a, p->h.n); + for (k = p->a.n = s = 0; k < ug->u.n; k++) { + p->h.a[k].a = p->a.n; + p->h.a[k].n = 0; + p->h.a[k].m = infer_mmhap_copy(ug, rg, src, ff, k, het_cov, asm_opt.polyploidy); + if(p->h.a[k].m == (uint64_t)asm_opt.polyploidy) { + p->h.a[k].n = p->h.a[k].m; s = 1; + } + p->a.n += p->h.a[k].m; bs[p->h.a[k].m] += ug->g->seq[k].len; + } + free(ff); + p->a.m = p->a.n; MALLOC(p->a.a, p->a.n); memset(p->a.a, -1, sizeof((*(p->a.a)))*p->a.n); + if(s) { + for (k = p->a.n = 0; k < ug->u.n; k++) { + if(p->h.a[k].n == p->h.a[k].m) { + for (z = 0; z < p->h.a[k].m; z++) p->a.a[p->h.a[k].a+z] = z; + } + } + } -void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, uint32_t is_poy, kvec_pe_hit **rhits) + for (k = 1; k <= (uint64_t)asm_opt.polyploidy; k++) { + fprintf(stderr, "[M::stat] # %lu-copy bases: %lu\n", k, bs[k]); + } + free(bs); + return p; +} + +bubble_type *gen_mmhap_bub(ma_ug_t* ug, uint8_t *r_het_flag, kv_u_trans_t *ref, mmhap_t *hh) +{ + uint64_t k; + bubble_type *p; CALLOC(p, 1); p->n_round = asm_opt.n_weight; p->round_id = 0; + identify_bubbles(ug, p, r_het_flag, ref); + for (k = 0; k < ug->g->n_seq; k++) { + if(IF_BUB(k, (*p))) continue; + if(IF_HOM(k, (*p))) { + if(hh->h.a[k].n < hh->h.a[k].m) p->index[k] = p->num.n; + } else if(IF_HET(k, (*p))) { + if(hh->h.a[k].n >= hh->h.a[k].m) p->index[k] = (uint32_t)-1; + } + } + return p; +} + +void purge_phase_0(ha_ug_index* idx, hc_links *link, bubble_type *bub, ps_t *s, kvec_pe_hit* hits, kv_u_trans_t *k_trans, uint8_t *del, mmhap_t *hh) +{ + k_trans->idx.n = k_trans->n = 0; hits->idx.n = 0; + for (bub->round_id = 0; bub->round_id < bub->n_round; bub->round_id++) { + renew_kv_u_trans(k_trans, link, hits, &(idx->t_ch->k_trans), idx, bub, s->s, NULL, 1/**0**/); + mc_solve(NULL, NULL, k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, + (bub->round_id==0?1:0), s->s, 1, bub, &(idx->t_ch->k_trans), 0, + /**(((bub.round_id+1) == bub.n_round)?1:0)**/0); + /*******************************for debug************************************/ + // label_unitigs_sm(s->s, NULL, idx->ug); + } + + uint64_t l[2], k, len; int64_t p; ma_ug_t *ug = idx->ug; + for (k = l[0] = l[1] = 0; k < ug->g->n_seq; k++) { + if((del[k] == 1) || (s->s[k] == 0)) continue; + if(s->s[k] > 0) l[0] += ug->g->seq[k].len; + else l[1] += ug->g->seq[k].len; + } + + p = (l[0]>=l[1])?1:-1; + for (k = len = 0; k < ug->g->n_seq; k++) { + if((del[k] == 1) || (s->s[k] == 0)) { + if((del[k] == 2) || + ((hh->h.a[k].n == hh->h.a[k].m) && (hh->h.a[k].m == ((uint64_t)asm_opt.polyploidy)))) { + len += ug->g->seq[k].len; + } + continue; + } + if(s->s[k] == p) { + del[k] = 2; len += ug->g->seq[k].len; + } else { + del[k] = 1; bub->index[k] = (uint32_t)-1; + } + s->s[k] = 0;///reset + } + fprintf(stderr, "[M::%s::stat] # remaining bases: %lu\n", __func__, len); +} + +void exchange_kv_u_trans_t(kv_u_trans_t *a, kv_u_trans_t *b) +{ + uint64_t k, *ua; u_trans_t *u; + k = a->n; a->n = b->n; b->n = k; + k = a->m; a->m = b->m; b->m = k; + u = a->a; a->a = b->a; b->a = u; + + k = a->idx.n; a->idx.n = b->idx.n; b->idx.n = k; + k = a->idx.m; a->idx.m = b->idx.m; b->idx.m = k; + ua = a->idx.a; a->idx.a = b->idx.a; b->idx.a = ua; +} + +uint64_t cal_ave_ovlp(u_trans_t *a, uint64_t an, double top) +{ + if(!an) return 0; + uint64_t k, len, cut, occ, tot; + for (k = len = 0; k < an; k++) { + len += a[k].qe - a[k].qs; + } + cut = len - (len*top); + for (k = occ = tot = 0; k < an; k++) { + if(a[k].qe - a[k].qs < cut) continue; + occ++; tot += a[k].qe - a[k].qs; + } + if(!occ) { + occ = an; tot = len; + } + return tot/occ; +} + + +void clean_trans_ovlp(bubble_type *bub, kv_u_trans_t *des, kv_u_trans_t *src, asg64_v *srt, ma_ug_t *ug) +{ + uint64_t k, st, i, m = 0, ncut; + for (k = des->n = 0; k < src->n; k++) { + if(IF_HOM(src->a[k].qn, (*bub))) continue; + if(IF_HOM(src->a[k].tn, (*bub))) continue; + kv_push(u_trans_t, *des, src->a[k]); + } + + kv_resize(uint64_t, des->idx, src->idx.n); des->idx.n = src->idx.n; + memset(des->idx.a, 0, des->idx.n*sizeof((*(des->idx.a)))); + for (st = 0, i = 1; i <= des->n; ++i) { + if (i == des->n || des->a[i].qn != des->a[st].qn) { + des->idx.a[des->a[st].qn] = (((uint64_t)st)<<32)|(i-st); m++; + st = i; + } + } + + if(srt && m) { + ncut = 0; kv_resize(uint64_t, *srt, m); + for (k = srt->n = 0; k < des->idx.n; k++) { + if(!((uint32_t)(des->idx.a[k]))) continue; + m = cal_ave_ovlp(des->a+(des->idx.a[k]>>32), ((uint32_t)(des->idx.a[k])), 0.9); + m = ((uint64_t)-1)-m; m <<= 32; m += k; + kv_push(uint64_t, *srt, m); + } + radix_sort_b64(srt->a, srt->a + srt->n); + for (k = 0; k < srt->n; k++) { + ncut += trans_sec_cut0(des, srt, (uint32_t)(srt->a[k]), 0.2, 256, ug); + } + + if(ncut) {///renew idx + for (i = k = 0; i < des->n; i++) { + if(des->a[i].del) continue; + des->a[k++] = des->a[i]; + } + des->n = k; + + memset(des->idx.a, 0, des->idx.n*sizeof((*(des->idx.a)))); + for (st = 0, i = 1; i <= des->n; ++i) { + if (i == des->n || des->a[i].qn != des->a[st].qn) { + des->idx.a[des->a[st].qn] = (((uint64_t)st)<<32)|(i-st); m++; + st = i; + } + } + } + } +} + +void purge_phase(ha_ug_index* idx, mmhap_t *hh, uint64_t hapid, uint64_t max_round, hc_links *link, kvec_pe_hit* hits, +bubble_type *bub, ps_t *s, kv_u_trans_t *k_trans, uint8_t *del, uint32_t *bidx, kv_u_trans_t *ref, asg64_v *srt) +{ + uint64_t k, len; uint32_t *bm; + ma_ug_t *ug = idx->ug; + if(max_round <= 0) { + for (k = 0; k < ug->g->n_seq; k++) { + if(hh->h.a[k].n >= hh->h.a[k].m) continue; + hh->a.a[hh->h.a[k].a+hh->h.a[k].n] = hapid; + hh->h.a[k].n++; + } + return; + } + + for (k = 0; k < ug->g->n_seq; k++) { + del[k] = 0; s->s[k] = 0; + if(hh->h.a[k].n >= hh->h.a[k].m) {///all haplotypes have been set + bub->index[k] = (uint32_t)-1; del[k] = 1; + } else { + if(IF_HOM(k, (*bub))) bub->index[k] = bub->num.n; + } + bidx[k] = bub->index[k]; + } + bm = bub->index; bub->index = bidx; bidx = bm; + + kv_resize(u_trans_t, *ref, idx->t_ch->k_trans.n); ref->n = idx->t_ch->k_trans.n; + memcpy(ref->a, idx->t_ch->k_trans.a, ref->n*sizeof((*(ref->a)))); + kv_resize(uint64_t, ref->idx, idx->t_ch->k_trans.idx.n); ref->idx.n = idx->t_ch->k_trans.idx.n; + memcpy(ref->idx.a, idx->t_ch->k_trans.idx.a, ref->idx.n*sizeof((*(ref->idx.a)))); + exchange_kv_u_trans_t(ref, &(idx->t_ch->k_trans)); + + for (k = 0; k < max_round; k++) { + clean_trans_ovlp(bub, &(idx->t_ch->k_trans), ref, ((k+1)g->n_seq; k++) { + if(del[k] != 2) { + if((hh->h.a[k].n == hh->h.a[k].m) && (hh->h.a[k].m == ((uint64_t)asm_opt.polyploidy))) { + len += ug->g->seq[k].len; + } + continue; + } + assert(hh->h.a[k].n < hh->h.a[k].m); + hh->a.a[hh->h.a[k].a+hh->h.a[k].n] = hapid; + hh->h.a[k].n++; len += ug->g->seq[k].len; + } + fprintf(stderr, "[M::%s::stat] # hap%lu bases: %lu\n", __func__, hapid+1, len); + + bm = bub->index; bub->index = bidx; bidx = bm; + exchange_kv_u_trans_t(ref, &(idx->t_ch->k_trans)); +} + +int hic_short_align_mmhap(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_opt_t *opt, kvec_pe_hit **rhits, mmhap_t **rh) +{ + if(rh) (*rh) = NULL; + sldat_t sl; + sl.idx = idx; + sl.t_ch = idx->t_ch; + sl.chunk_size = 20000000; + sl.n_thread = asm_opt.thread_num; + sl.total_base = sl.total_pair = 0; + idx->hap_cnt = asm_opt.hap_occ; + kv_init(sl.hits.a); kv_init(sl.hits.idx); kv_init(sl.hits.occ); + + + if(!load_hc_hits(&sl.hits, idx->ug, asm_opt.output_file_name)) { + alignment_worker_pipeline(&sl, fn1, fn2); + write_hc_hits(&sl.hits, idx->ug, asm_opt.output_file_name); + } + sl.hits.uID_bits = idx->uID_bits; sl.hits.pos_mode = idx->pos_mode; + + if(sl.hits.idx.n == 0) idx_hc_links(&(sl.hits), idx, NULL); + + mmhap_t *hh = gen_mmhap_t(idx->ug, idx->read_g, opt->sources); + bubble_type *bub = gen_mmhap_bub(idx->ug, idx->t_ch->ir_het, &(idx->t_ch->k_trans), hh); + hc_links link; init_hc_links(&link, idx->ug->g->n_seq, idx->t_ch); + measure_distance(idx, idx->ug, &sl.hits, &link, bub, &(idx->t_ch->k_trans)); + kv_u_trans_t k_trans; kv_init(k_trans); kv_init(k_trans.idx); + ps_t *s = init_ps_t(11, idx->ug->g->n_seq); ///H_partition hap; + uint64_t k, n_hap = asm_opt.polyploidy; asg64_v srt; kv_init(srt); + uint8_t *ff; CALLOC(ff, idx->ug->g->n_seq); + uint32_t *bidx; MALLOC(bidx, idx->ug->g->n_seq); + kv_u_trans_t r_trans_buf; kv_init(r_trans_buf); kv_init(r_trans_buf.idx); + + for (k = 0; k < n_hap; k++) { + purge_phase(idx, hh, k, ((n_hap>k)?(n_hap-k-1):(0)), &link, &sl.hits, bub, s, &k_trans, ff, bidx, &r_trans_buf, &srt); + } + if(rh) (*rh) = hh; + + // if(rhits) (*rhits) = get_r_hits_order(&sl.hits, idx->uID_bits, idx->pos_mode, idx->read_g, idx->ug, &bub); + + kv_destroy(sl.hits.a); kv_destroy(sl.hits.idx); kv_destroy(sl.hits.occ); + destory_hc_links(&link); + kv_destroy(k_trans); kv_destroy(k_trans.idx); + destory_ps_t(&s); destory_bubbles(bub); free(bub); + kv_destroy(srt); free(ff); free(bidx); + kv_destroy(r_trans_buf); kv_destroy(r_trans_buf.idx); + return 1; +} + + +void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, mmhap_t **rh, kvec_pe_hit **rhits) { ug_index = NULL; int exist = (asm_opt.load_index_from_disk? @@ -17442,8 +17714,9 @@ void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, ug_index->read_g = read_g; ug_index->t_ch = t_ch; ///test_unitig_index(ug_index, ug); - if(!is_poy) hic_short_align(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt, rhits); - else hic_short_align_poy(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt); + if(!rh) hic_short_align(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt, rhits); + else hic_short_align_mmhap(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt, rhits, rh); + // else hic_short_align_poy(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt); destory_hc_pt_index(ug_index);free(ug_index); diff --git a/hic.h b/hic.h index 321c044..61760de 100644 --- a/hic.h +++ b/hic.h @@ -113,7 +113,7 @@ pdq* pqw, uint32_t* path_w, buf_t *resw, asg_t *sg, uint8_t *dest, uint8_t df, u long long *dis); void set_utg_by_dis(uint32_t v, pdq* pq, asg_t *g, kvec_t_u32_warp *res, uint32_t dis); void dedup_hits(kvec_pe_hit* hits, uint64_t is_dup); -void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, uint32_t is_poy, kvec_pe_hit **rhits); +void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, mmhap_t **rh, kvec_pe_hit **rhits); spg_t *hic_pre_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, kvec_pe_hit **rhits); void prt_bubble_gfa_adv(FILE *fp, bubble_type *bub, const char* utg_pre, const char* bub_pre, const char* chain_pre); void bp_solve(ug_opt_t *opt, kv_u_trans_t *ref, ma_ug_t *ug, asg_t *sg, bubble_type *bub, double cis_rate); diff --git a/inter.cpp b/inter.cpp index eb49c7d..2cf8218 100644 --- a/inter.cpp +++ b/inter.cpp @@ -64,6 +64,8 @@ int64_t ug_map_lchain_simple(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl // #define GBIN_L 15000 #define GBIN_L 256 +#define FREE_BATCH 16 + #define generic_key(x) (x) KRADIX_SORT_INIT(gfa64, uint64_t, generic_key, 8) @@ -417,6 +419,8 @@ typedef struct { // global data structure for kt_pipeline() int32_t is_HPC, bw, max_gap, chn_pen_gap, n_thread, is_cnt, is_ovlp, mini_cut, chain_cut, keep_unsymm_arc; ul_idx_t udb; kv_u_trans_t *filter; + uint32_t *free_cnt; + ug_rid_cov_t *ccov; } ug_trans_t; typedef struct { // global data structure for kt_pipeline() @@ -441,6 +445,7 @@ typedef struct { // global data structure for kt_pipeline() // bit_mask_t *bm; } ug_bin_t; + void hc_glchain_destroy(glchain_t *b) { if (!b) return; @@ -9573,7 +9578,8 @@ static void worker_for_trans_ovlp(void *data, long i, int tid) // callback for k } -void filter_by_reliable_ovlp_adv(uint32_t id, kv_u_trans_t *idx, st_mt_t *sp, overlap_region_alloc* ol, const ul_idx_t *udb, double sec_rate, uint64_t avoid_dup_aln, uint64_t *occ1) +void filter_by_reliable_ovlp_adv(uint32_t id, kv_u_trans_t *idx, st_mt_t *sp, overlap_region_alloc* ol, const ul_idx_t *udb, double sec_rate, uint64_t avoid_dup_aln, +uint64_t dedup_by_reliable_ovlp, uint64_t *occ1) { (*occ1) = 0; u_trans_t *a; uint64_t n, k, l, s, e, s0, e0, z, ov, rr, r1, spn; overlap_region *m, t; @@ -9629,25 +9635,28 @@ void filter_by_reliable_ovlp_adv(uint32_t id, kv_u_trans_t *idx, st_mt_t *sp, ov } } - for (k = sp->n = 0; k < n; k++) { - if(a[k].del) continue; - if(a[k].f == RC_0 || a[k].f == RC_1) { - kv_push(uint64_t, *sp, (((uint64_t)a[k].qs)<<32)|((uint64_t)a[k].qe)); - } - } - if(sp->n > 1) { - radix_sort_gfa64(sp->a, sp->a + sp->n); - for (k = z = 0; k < sp->n; k++) { - s = sp->a[k]>>32; e = (uint32_t)sp->a[k]; - if(z > 0 && s <= ((uint32_t)sp->a[z-1])) { - if(e > ((uint32_t)sp->a[z-1])) { - sp->a[z-1] >>= 32; sp->a[z-1] <<= 32; sp->a[z-1] |= e; - } - } else { - sp->a[z++] = sp->a[k]; + sp->n = 0; + if(dedup_by_reliable_ovlp) { + for (k = sp->n = 0; k < n; k++) { + if(a[k].del) continue; + if(a[k].f == RC_0 || a[k].f == RC_1) { + kv_push(uint64_t, *sp, (((uint64_t)a[k].qs)<<32)|((uint64_t)a[k].qe)); } } - sp->n = z; + if(sp->n > 1) { + radix_sort_gfa64(sp->a, sp->a + sp->n); + for (k = z = 0; k < sp->n; k++) { + s = sp->a[k]>>32; e = (uint32_t)sp->a[k]; + if(z > 0 && s <= ((uint32_t)sp->a[z-1])) { + if(e > ((uint32_t)sp->a[z-1])) { + sp->a[z-1] >>= 32; sp->a[z-1] <<= 32; sp->a[z-1] |= e; + } + } else { + sp->a[z++] = sp->a[k]; + } + } + sp->n = z; + } } for (k = rr = r1 = 0, spn = sp->n; k < ol->length; k++) { @@ -9659,7 +9668,7 @@ void filter_by_reliable_ovlp_adv(uint32_t id, kv_u_trans_t *idx, st_mt_t *sp, ov // ol->list[k].y_id+1, "lc"[udb->ug->u.a[ol->list[k].y_id].circ], udb->ug->u.a[ol->list[k].y_id].len, // ol->list[k].y_pos_s, ol->list[k].y_pos_e+1, ol->list[k].x_pos_strand); // } - if(m->x_pos_strand == 0) { + if((dedup_by_reliable_ovlp) && (m->x_pos_strand == 0)) { s = m->x_pos_s; e = m->x_pos_e + 1; l = 0; for (z = 0; z < spn; z++) { s0 = sp->a[z]>>32; e0 = (uint32_t)sp->a[z]; @@ -9743,8 +9752,8 @@ uint64_t split_ug_lalign(uint64_t ol_h, overlap_region_alloc* ol, double errh, d overlap_region *aux_o, int64_t sid, uint64_t khit, uint64_t chain_cut, void *km) { uint64_t ol_l = ol->length - ol_h, on0, k, m; int64_t wl; double erate; overlap_region t; - // if(sid == 160) fprintf(stderr, "[M::%s] errh::%f, errl::%f\n", __func__, errh, errl); - // if(sid == 160) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "st"); + // if(sid == 57) fprintf(stderr, "[M::%s] errh::%f, errl::%f\n", __func__, errh, errl); + // if(sid == 57) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "st"); if(ol_h) { erate = errh; on0 = ol->length; wl = MIN((((double)THRESHOLD_MAX_SIZE)/erate), WINDOW); @@ -9760,7 +9769,7 @@ uint64_t split_ug_lalign(uint64_t ol_h, overlap_region_alloc* ol, double errh, d } ol_h = ol->length; ol->length = m; ol_l = ol->length - ol_h; } - // if(sid == 160) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "mi"); + // if(sid == 57) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "mi"); if(ol_l) { erate = errl; on0 = ol->length; wl = MIN((((double)THRESHOLD_MAX_SIZE)/erate), WINDOW); @@ -9774,10 +9783,10 @@ uint64_t split_ug_lalign(uint64_t ol_h, overlap_region_alloc* ol, double errh, d m++; } } - // if(sid == 7) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "sw"); + // if(sid == 57) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "sw"); ol->length = ol_l; ug_lalign(ol, cl, uref, uopt, qstr, ql, qu, tu, dumy, exz, aux_o, erate, wl, sid, khit, chain_cut, km); - // if(sid == 7) { + // if(sid == 57) { // fprintf(stderr, "[M::%s::]\ton0::%lu\tol_l::%lu\tol_h::%lu\tol->length::%lu\n", __func__, // on0, ol_l, ol_h, ol->length); // prt_split_ovs(ol, uref->ug, ol_h, ol_l, "u0"); @@ -9791,7 +9800,7 @@ uint64_t split_ug_lalign(uint64_t ol_h, overlap_region_alloc* ol, double errh, d m++; } ol_l = ol->length; ol->length = m; ol_h = ol->length - ol_l; - // if(sid == 7) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "u1"); + // if(sid == 57) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "u1"); if(ol_l) {///swap ol_h and ol_l for (k = ol_l, m = 0; k < ol->length; k++) { if(k != m) { @@ -9803,7 +9812,7 @@ uint64_t split_ug_lalign(uint64_t ol_h, overlap_region_alloc* ol, double errh, d } } } - // if(sid == 160) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "ed"); + // if(sid == 57) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "ed"); return ol_h; } @@ -9870,6 +9879,50 @@ uint32_t test_het_aln(ma_ug_t *ug, uint64_t rid, u_trans_t *a, uint64_t a_n, st_ return 0; } +uint32_t is_mmhom_node(uint64_t *ca, ma_utg_t *u, asg_t *sg, uint64_t cov_bd, double cut_rate) +{ + if(cut_rate < 0) cut_rate = 0; if(cut_rate > 1.0) cut_rate = 1.0; + uint64_t k, a, na, a_cut = u->n*cut_rate, na_cut = u->n*(1.0-cut_rate); + for (k = a = na = 0; k < u->n; k++) { + if(ca[k] > (cov_bd*((uint64_t)sg->seq[u->a[k]>>33].len))) { + a++; if((a) && (a>=a_cut)) return 1; + } else { + na++; if((na) && (na>=na_cut)) return 0; + } + } + + if((a) && (a>=a_cut)) return 1; + return 0; +} + +uint32_t test_het_aln_mmhap(uint64_t uid, ug_rid_cov_t *ccov, u_trans_t *a, uint64_t a_n, overlap_region_alloc* ol, st_mt_t *sp) +{ + uint64_t k; u_trans_t p; + if(a_n == 0 && ol->length == 0) return 0; + kv_resize(uint64_t, *sp, ccov->ug->u.a[uid].n); + memcpy(sp->a, ccov->cov.a+ccov->idx[uid], sizeof((*(sp->a)))*ccov->ug->u.a[uid].n); + + for (k = 0; k < a_n; k++) { + if(a[k].del) continue; + append_cov_line_ug_rid_cov_t(uid, sp->a, &(a[k]), ccov, ((uint64_t)-1), -1); + } + for (k = 0; k < ol->length; k++) { + p.qn = uid; p.tn = ol->list[k].y_id; + p.rev = ol->list[k].y_pos_strand; p.f = RC_3; p.nw = 0; + p.qs = ol->list[k].x_pos_s; p.qe = ol->list[k].x_pos_e+1; + if(p.rev) { + p.ts = ccov->ug->u.a[p.tn].len - (ol->list[k].y_pos_e+1); + p.te = ccov->ug->u.a[p.tn].len - ol->list[k].y_pos_s; + } else { + p.ts = ol->list[k].y_pos_s; + p.te = ol->list[k].y_pos_e+1; + } + append_cov_line_ug_rid_cov_t(uid, sp->a, &p, ccov, ((uint64_t)-1), -1); + } + + return is_mmhom_node(sp->a, &(ccov->ug->u.a[uid]), ccov->rg, ccov->hom_min, 0.8); +} + void push_ul_ov_t(ul_idx_t *udb, u_trans_t *a, uint64_t a_n, uint64_t rid, st_mt_t *sp, overlap_region_alloc* ol, uint64_t len, uint64_t is_arc_filter, double max_err, kv_ul_ov_t *res) { uint64_t cnt, z, k, l, m, spn; ul_ov_t *p; @@ -9921,11 +9974,10 @@ void push_ul_ov_t(ul_idx_t *udb, u_trans_t *a, uint64_t a_n, uint64_t rid, st_mt // fprintf(stderr, ">0<[M::%s] utg%.6u%c -> utg%.6u%c\n", __func__, // p->qn+1, "lc"[s->udb.ug->u.a[p->qn].circ], // p->tn+1, "lc"[s->udb.ug->u.a[p->tn].circ]); - // if(i == 5) - // { - // fprintf(stderr, "***utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\ti::%ld\n", - // p->qn+1, "lc"[s->ug->u.a[p->qn].circ], s->ug->u.a[p->qn].len, p->qs, p->qe, "+-"[p->rev], - // p->tn+1, "lc"[s->ug->u.a[p->tn].circ], s->ug->u.a[p->tn].len, p->ts, p->te, i); + // if(p->ts >= p->te || p->qs >= p->qe) { + // fprintf(stderr, "+[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\n", __func__, + // p->qn+1, "lc"[udb->ug->u.a[p->qn].circ], udb->ug->u.a[p->qn].len, p->qs, p->qe, "+-"[p->rev], + // p->tn+1, "lc"[udb->ug->u.a[p->tn].circ], udb->ug->u.a[p->tn].len, p->ts, p->te); // } if((is_arc_filter) && (!trans_ovlp_connect(p, udb->ug))) res->n--; // fprintf(stderr, ">1<[M::%s] utg%.6u%c -> utg%.6u%c\n", __func__, @@ -9948,6 +10000,11 @@ void push_ul_ov_t(ul_idx_t *udb, u_trans_t *a, uint64_t a_n, uint64_t rid, st_mt p->qs = a[sp->a[k]].qs; p->qe = a[sp->a[k]].qe; p->ts = a[sp->a[k]].ts; p->te = a[sp->a[k]].te; p->sec = (p->qe-p->qs)*max_err; + // if(p->ts >= p->te || p->qs >= p->qe) { + // fprintf(stderr, "-[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\n", __func__, + // p->qn+1, "lc"[udb->ug->u.a[p->qn].circ], udb->ug->u.a[p->qn].len, p->qs, p->qe, "+-"[p->rev], + // p->tn+1, "lc"[udb->ug->u.a[p->tn].circ], udb->ug->u.a[p->tn].len, p->ts, p->te); + // } } // if(is_sec_filter) { @@ -10038,7 +10095,7 @@ uint64_t gen_trans_adaptive_aln(ug_trans_t *s, uint64_t rid, ha_ovec_buf_t *b, k // fprintf(stderr, "-2-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); // } - filter_by_reliable_ovlp_adv(rid, s->filter, &(b->sp), &b->olist, &(s->udb), s->sec_cutoff, 1, &ol_h); + filter_by_reliable_ovlp_adv(rid, s->filter, &(b->sp), &b->olist, &(s->udb), s->sec_cutoff, 1, 1, &ol_h); clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2); if(!fi) ol_h = 0; @@ -10064,6 +10121,47 @@ uint64_t gen_trans_adaptive_aln(ug_trans_t *s, uint64_t rid, ha_ovec_buf_t *b, k return pass_aln; } +void clear_count_buf(ug_trans_t *s, uint32_t tid, uint32_t free_count) +{ + // fprintf(stderr, "[M::%s]\n", __func__); + ha_ovec_buf_t *b = s->hab[tid]; + destory_fake_cigar(&(b->tmp_region.f_cigar)); + free(b->tmp_region.w_list.a); free(b->tmp_region.w_list.c.a); + memset(&(b->tmp_region), 0, sizeof(b->tmp_region)); + init_fake_cigar(&(b->tmp_region.f_cigar)); + memset(&(b->tmp_region.w_list), 0, sizeof(b->tmp_region.w_list)); + CALLOC(b->tmp_region.w_list.a, 1); b->tmp_region.w_list.n = b->tmp_region.w_list.m = 1; + + ha_abufl_destroy(b->abl); b->abl = ha_abufl_init(); + + kv_destroy(b->sp); memset(&(b->sp), 0, sizeof((b->sp))); + if(free_count) return; + + destory_Candidates_list(&b->clist); + memset((&(b->clist)), 0, sizeof(b->clist)); + init_Candidates_list(&b->clist); + + destory_overlap_region_alloc(&b->olist); + memset((&(b->olist)), 0, sizeof(b->olist)); + init_overlap_region_alloc(&b->olist); + + destory_UC_Read(&b->self_read); + memset((&(b->self_read)), 0, sizeof(b->self_read)); + init_UC_Read(&b->self_read); + + destory_UC_Read(&b->ovlp_read); + memset((&(b->ovlp_read)), 0, sizeof(b->ovlp_read)); + init_UC_Read(&b->ovlp_read); + + destory_Correct_dumy(&b->correct); + memset((&(b->correct)), 0, sizeof(b->correct)); + init_Correct_dumy(&b->correct); + + destroy_bit_extz_t(&(b->exz)); + memset((&(b->exz)), 0, sizeof(b->exz)); + init_bit_extz_t(&(b->exz), 31); +} + static void worker_for_trans_ovlp_adv(void *data, long i, int tid) // callback for kt_for() { ug_trans_t *s = (ug_trans_t*)data; @@ -10085,6 +10183,10 @@ static void worker_for_trans_ovlp_adv(void *data, long i, int tid) // callback f s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3, s->is_HPC, s->idx_a.a + s->idx_n.a[i], 0, NULL, 0, s->mini_cut, s->chain_cut, NULL); assert(cnt == ((s->idx_n.a[i+1]-s->idx_n.a[i]))); } + if(s->free_cnt[tid] >= FREE_BATCH) { + clear_count_buf(s, tid, 1); s->free_cnt[tid] = 0; + } + s->free_cnt[tid]++; return; } @@ -10096,6 +10198,240 @@ static void worker_for_trans_ovlp_adv(void *data, long i, int tid) // callback f if(!gen_trans_adaptive_aln(s, i, b, bl, seq, len, s->filter, s->diff_ec_ul, s->diff_ec_ul_double, s->bw_thres, s->bw_thres_double)) { gen_trans_adaptive_aln(s, i, b, bl, seq, len, NULL, s->diff_ec_ul_double, s->diff_ec_ul_double, s->bw_thres_double, s->bw_thres_double); } + + if(s->free_cnt[tid] >= FREE_BATCH) { + clear_count_buf(s, tid, 0); s->free_cnt[tid] = 0; + } + s->free_cnt[tid]++; +} + +uint64_t *gen_reliable_cov_arr(uint32_t id, kv_u_trans_t *idx, ug_rid_cov_t *ccov, st_mt_t *sp) +{ + u_trans_t *a; uint64_t n, k; + a = u_trans_a(*idx, id); n = u_trans_n(*idx, id); + kv_resize(uint64_t, *sp, ccov->ug->u.a[id].n); + memcpy(sp->a, ccov->cov.a+ccov->idx[id], sizeof((*(sp->a)))*ccov->ug->u.a[id].n); + for (k = 0; k < n; k++) { + if(a[k].del) continue; + if(a[k].f == RC_0 || a[k].f == RC_1) { + append_cov_line_ug_rid_cov_t(id, sp->a, &(a[k]), ccov, ((uint64_t)-1), -1); + } + } + return sp->a; +} + +uint64_t is_above_cov(uint64_t uid, overlap_region *o, uint64_t *fc, ug_rid_cov_t *ccov, double sec_rate) +{ + u_trans_t p; + p.qn = uid; p.tn = o->y_id; + p.rev = o->y_pos_strand; p.f = RC_3; p.nw = 0; + p.qs = o->x_pos_s; p.qe = o->x_pos_e+1; + if(p.rev) { + p.ts = ccov->ug->u.a[p.tn].len - (o->y_pos_e+1); + p.te = ccov->ug->u.a[p.tn].len - o->y_pos_s; + } else { + p.ts = o->y_pos_s; + p.te = o->y_pos_e+1; + } + + if(append_cov_line_ug_rid_cov_t(uid, fc, &p, ccov, ccov->hom_max, sec_rate)) return 0; + return 1; +} + +void filter_by_reliable_ovlp_mmhap_adv(uint32_t id, kv_u_trans_t *idx, st_mt_t *sp, overlap_region_alloc* ol, const ul_idx_t *udb, double sec_rate, uint64_t avoid_dup_aln, +uint64_t dedup_by_reliable_ovlp, ug_rid_cov_t *ccov, uint64_t *occ1) +{ + (*occ1) = 0; + u_trans_t *a; uint64_t n, k, l, z, rr, r1, *fc; overlap_region *m, t; + a = u_trans_a(*idx, id); n = u_trans_n(*idx, id); + if(avoid_dup_aln) { + kv_resize(uint64_t, *sp, (ol->length)+n); + for (k = sp->n = 0; k < n; k++) { + if(a[k].del) continue; + if(a[k].f == RC_0 || a[k].f == RC_1) { + z = a[k].tn; z <<= 1; z |= a[k].rev; z <<= 32; + kv_push(uint64_t, *sp, z); + } + } + if(sp->n > 0) { + for (k = 0; k < ol->length; k++) { + z = ol->list[k].y_id; z <<= 1; z |= ol->list[k].y_pos_strand; + z <<= 32; z |= k; z |= ((uint64_t)0x80000000); + kv_push(uint64_t, *sp, z); + } + + radix_sort_gfa64(sp->a, sp->a + sp->n); + for (k = 1, l = 0, rr = 0; k <= sp->n; k++) { + if(k == sp->n || (sp->a[l]>>32)!=(sp->a[k]>>32)) { + if((k - l > 1) && (!(sp->a[l]&((uint64_t)0x80000000)))) {///overlap within bck + for (z = l; z < k; z++) { + if(sp->a[z]&((uint64_t)0x80000000)) { + ol->list[(uint32_t)(sp->a[z]-((uint64_t)0x80000000))].y_id = ((uint32_t)-1); + rr++; + } + } + } + l = k; + } + } + + if(rr > 0) { + for (k = rr = 0; k < ol->length; k++) { + if(ol->list[k].y_id == ((uint32_t)-1)) continue; + if(rr != k) { + t = ol->list[rr]; + ol->list[rr] = ol->list[k]; + ol->list[k] = t; + } + rr++; + } + ol->length = rr; + } + } + } + + sp->n = 0; fc = NULL; + if(dedup_by_reliable_ovlp) { + for (k = 0; k < n; k++) { + if(a[k].del) continue; + if(a[k].f == RC_0 || a[k].f == RC_1) break; + } + if(k < n) fc = gen_reliable_cov_arr(id, idx, ccov, sp); + } + + for (k = rr = r1 = 0; k < ol->length; k++) { + m = &(ol->list[k]); + if((dedup_by_reliable_ovlp) && (m->x_pos_strand == 0) && (fc)) { + if(is_above_cov(id, m, fc, ccov, sec_rate)) continue; + } + + if(rr != k) { + t = ol->list[k]; + ol->list[k] = ol->list[rr]; + ol->list[rr] = t; + } + if(ol->list[rr].x_pos_strand) { + ol->list[rr].x_pos_strand = 0; + if(r1 != rr) { + t = ol->list[r1]; + ol->list[r1] = ol->list[rr]; + ol->list[rr] = t; + } + r1++; + } + rr++; + } + // if(id == 1576) { + // fprintf(stderr, "[M::%s] utg%.6ul, ol->length0::%lu, ol->length::%lu\n", __func__, id+1, ol->length, rr); + // } + ol->length = rr; (*occ1) = r1; +} + + +uint64_t gen_trans_adaptive_mmhap_aln(ug_trans_t *s, uint64_t rid, ha_ovec_buf_t *b, kv_ul_ov_t *bl, char *seq, uint64_t len, kv_u_trans_t *fi, double err_low, double err_high, double bw_low, double bw_high) +{ + uint64_t cnt = ((s->idx_n.a[rid+1]-s->idx_n.a[rid])), ol_h = 0, pass_aln = 0; + uint32_t high_occ = asm_opt.polyploidy + 1; overlap_region *aux_o = NULL; + ///note: high_occ is different + ug_map_lchain(b->abl, rid, seq, len, s->w, s->k, &(s->udb), &b->olist, &b->clist, bw_low, bw_high, + s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3, + s->is_HPC, s->idx_a.a + s->idx_n.a[rid], cnt, s->srt_a.a, s->srt_a.n, s->mini_cut, s->chain_cut, fi); + // if(rid == 57) { + // fprintf(stderr, "-1-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); + // } + ///remove candidate chains that have been calculated + if(!fi) backward_dedup_ol(rid, bl, &(b->sp), &b->olist);///it is ok + // if(rid == 57) { + // fprintf(stderr, "-2-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); + // } + filter_by_reliable_ovlp_mmhap_adv(rid, s->filter, &(b->sp), &b->olist, &(s->udb), s->sec_cutoff, 1, 1, s->ccov, &ol_h); + clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2); + if(!fi) ol_h = 0; + + // if(rid == 57) { + // fprintf(stderr, "-3-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); + // } + + ol_h = split_ug_lalign(ol_h, &b->olist, err_high, err_low, + &b->clist, &(s->udb), s->uopt, seq, len, &b->self_read, &b->ovlp_read, + &b->correct, &b->exz, aux_o, rid, s->k, s->chain_cut, NULL); + + // if(rid == 57) { + // fprintf(stderr, "-4-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); + // } + + aux_o = gen_aux_ovlp(&b->olist);///must be here + + // if(rid == 57) { + // fprintf(stderr, "-5-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); + // } + + ol_h = split_ug_lalign(ol_h, &b->olist, err_high, err_low, + &b->clist, &(s->udb), s->uopt, seq, len, &b->self_read, &b->ovlp_read, + &b->correct, &b->exz, aux_o, rid, s->k, s->chain_cut, NULL); + + // if(rid == 57) { + // fprintf(stderr, "-6-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); + // } + + if(fi) {///first round + pass_aln = test_het_aln_mmhap(rid, s->ccov, u_trans_a((*(s->filter)), rid), u_trans_n((*(s->filter)), rid), &b->olist, &(b->sp)); + push_ul_ov_t(&(s->udb), u_trans_a((*(s->filter)), rid), u_trans_n((*(s->filter)), rid), rid, &(b->sp), &b->olist, len, pass_aln, err_high, bl); + // fprintf(stderr, "-1-[M::%s] utg%.6lu%c, rid::%lu, pass_aln::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, pass_aln); + } else {///second round + push_ul_ov_t(&(s->udb), NULL, 0, rid, &(b->sp), &b->olist, len, 0, err_high, bl); + remove_trans_ovlp_connect(s->udb.ug, rid, bl); + } + return pass_aln; +} + +static void worker_for_trans_ovlp_mmhap_adv(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]; kv_ul_ov_t *bl = &(s->ll[tid].tk); + uint32_t high_occ = asm_opt.polyploidy + 1; uint64_t cnt; + char *seq = s->ug->u.a[i].s; int64_t len = s->ug->u.a[i].len; + if((!s->is_ovlp) && (s->is_cnt)) s->idx_n.a[i] = 0; + if(s->ug->g->seq[i].del) return; + if(is_mmhom_node(s->ccov->cov.a+s->ccov->idx[i], &(s->ug->u.a[i]), s->ccov->rg, s->ccov->hom_min, 0.9)) return; + // 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; + + if(!s->is_ovlp) { + if(s->is_cnt) { + s->idx_n.a[i] = ug_map_lchain(b->abl, i, seq, len, s->w, s->k, &(s->udb), NULL, NULL, s->bw_thres, s->bw_thres_double, + s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3, s->is_HPC, NULL, 0, NULL, 0, s->mini_cut, s->chain_cut, NULL); + } else { + cnt = ug_map_lchain(b->abl, i, seq, len, s->w, s->k, &(s->udb), NULL, NULL, s->bw_thres, s->bw_thres_double, + s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3, s->is_HPC, s->idx_a.a + s->idx_n.a[i], 0, NULL, 0, s->mini_cut, s->chain_cut, NULL); + assert(cnt == ((s->idx_n.a[i+1]-s->idx_n.a[i]))); + } + if(s->free_cnt[tid] >= FREE_BATCH) { + clear_count_buf(s, tid, 1); s->free_cnt[tid] = 0; + } + s->free_cnt[tid]++; + return; + } + + // if(i == 58) { + // fprintf(stderr, "\n-1-[M::%s] utg%.6u%c, rid::%ld, is_ovlp::%d, is_cnt::%d, len::%ld, str::%u\n", + // __func__, (uint32_t)i+1, "lc"[s->ug->u.a[i].circ], i, s->is_ovlp, s->is_cnt, len, (uint32_t)(!!seq)); + // } + + if(!gen_trans_adaptive_mmhap_aln(s, i, b, bl, seq, len, s->filter, s->diff_ec_ul, s->diff_ec_ul_double, s->bw_thres, s->bw_thres_double)) { + gen_trans_adaptive_mmhap_aln(s, i, b, bl, seq, len, NULL, s->diff_ec_ul_double, s->diff_ec_ul_double, s->bw_thres_double, s->bw_thres_double); + } + if(s->free_cnt[tid] >= FREE_BATCH) { + clear_count_buf(s, tid, 0); s->free_cnt[tid] = 0; + } + s->free_cnt[tid]++; } @@ -19562,7 +19898,7 @@ void clear_all_ul_t(all_ul_t *x) void init_ug_trans_t(ug_trans_t *opt, ug_opt_t *uopt, int32_t is_HPC, int32_t k, int32_t w, int32_t max_n_chain, double bw_thres, double diff_ec_ul, double bw_thres_double, double diff_ec_ul_double, double sec_cutoff, int32_t n_thread, -int32_t mini_cut, int32_t chain_cut, int32_t keep_unsymm_arc, ma_ug_t *ug, asg_t *sg, bubble_type *bub) +int32_t mini_cut, int32_t chain_cut, int32_t keep_unsymm_arc, ma_ug_t *ug, asg_t *sg, bubble_type *bub, uint8_t gen_bub) { int64_t i; uint8_t *bf = NULL; memset(opt, 0, sizeof((*opt))); @@ -19586,6 +19922,7 @@ int32_t mini_cut, int32_t chain_cut, int32_t keep_unsymm_arc, ma_ug_t *ug, asg_t opt->rg = sg; opt->n_thread = ((n_thread>=1)?n_thread:1); + CALLOC(opt->free_cnt, opt->n_thread); CALLOC(opt->hab, opt->n_thread); CALLOC(opt->ll, opt->n_thread); for (i = 0; i < opt->n_thread; ++i) { @@ -19598,7 +19935,7 @@ int32_t mini_cut, int32_t chain_cut, int32_t keep_unsymm_arc, ma_ug_t *ug, asg_t opt->udb.ug = ug; if(bub) { opt->bub = bub; - } else { + } else if(gen_bub) { opt->bub = gen_bubble_chain(sg, ug, uopt, &bf, 0); free(bf); } } @@ -19727,6 +20064,7 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res) clean_u_trans_t_idx_adv(res, p->ug, p->rg); p->filter = res; // fprintf(stderr, "[M::%s::] ==> 0\n", __func__); p->is_cnt = 1; p->is_ovlp = 0; + memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread); kt_for(p->n_thread, worker_for_trans_ovlp_adv, p, p->ug->u.n); for (i = l = 0; i < p->ug->u.n; i++) { occ = p->idx_n.a[i]; p->idx_n.a[i] = l; l += occ; @@ -19736,6 +20074,7 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res) p->idx_a.n = p->idx_a.m = l; MALLOC(p->idx_a.a, p->idx_a.n); // fprintf(stderr, "[M::%s::] ==> 1\n", __func__); p->is_cnt = 0; p->is_ovlp = 0; + memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread); kt_for(p->n_thread, worker_for_trans_ovlp_adv, p, p->ug->u.n); p->srt_a.n = p->srt_a.m = p->idx_a.n; MALLOC(p->srt_a.a, p->srt_a.n); // fprintf(stderr, "[M::%s::] p->idx_a.n::%lu \n", __func__, (uint64_t)p->idx_a.n); @@ -19770,13 +20109,14 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res) // fprintf(stderr, "[M::%s::] ==> 2\n", __func__); p->is_cnt = 0; p->is_ovlp = 1; + memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread); kt_for(p->n_thread, worker_for_trans_ovlp_adv, p, p->ug->u.n); // fprintf(stderr, "[M::%s::] ==> 3\n", __func__); for (i = 0; (int64_t)i < p->n_thread; i++) { ha_ovec_destroy(p->hab[i]); free(p->ll[i].lo.a); free(p->ll[i].srt.a.a); free(p->ll[i].tc.a); } - free(p->idx_a.a); free(p->idx_n.a); free(p->hab); + free(p->idx_a.a); free(p->idx_n.a); free(p->hab); free(p->free_cnt); // destory_bubbles(p->bub); free(p->bub); // fprintf(stderr, "[M::%s::] ==> 4\n", __func__); @@ -19838,13 +20178,144 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res) fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time); } + +void gen_trans_base_count_mmhap_comp(ug_trans_t *p, kv_u_trans_t *res) +{ + double index_time = yak_realtime(); + // ha_flt_tab = NULL; + uint64_t i, k, l, occ, m, cc; kv_ul_ov_t *bl = NULL; + u_trans_t *z; ha_mzl_t *tz; double ww; + p->ccov = gen_ug_rid_cov_t(p->ug, p->rg, p->uopt->sources); + clean_u_trans_t_idx_adv(res, p->ug, p->rg); p->filter = res; + // fprintf(stderr, "[M::%s::] ==> 0\n", __func__); + p->is_cnt = 1; p->is_ovlp = 0; + memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread); + kt_for(p->n_thread, worker_for_trans_ovlp_mmhap_adv, p, p->ug->u.n); + for (i = l = 0; i < p->ug->u.n; i++) { + occ = p->idx_n.a[i]; p->idx_n.a[i] = l; l += occ; + } + // fprintf(stderr, "[M::%s::] i::%lu, l::%lu\n", __func__, i, l); + p->idx_n.a[i] = l; + p->idx_a.n = p->idx_a.m = l; MALLOC(p->idx_a.a, p->idx_a.n); + // fprintf(stderr, "[M::%s::] ==> 1\n", __func__); + p->is_cnt = 0; p->is_ovlp = 0; + memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread); + kt_for(p->n_thread, worker_for_trans_ovlp_mmhap_adv, p, p->ug->u.n); + p->srt_a.n = p->srt_a.m = p->idx_a.n; MALLOC(p->srt_a.a, p->srt_a.n); + // fprintf(stderr, "[M::%s::] p->idx_a.n::%lu \n", __func__, (uint64_t)p->idx_a.n); + // memcpy(p->srt_a.a, p->idx_a.a, p->srt_a.n*sizeof((*(p->srt_a.a)))); + for (i = 0; i < p->srt_a.n; i++) { + p->srt_a.a[i] = p->idx_a.a[i]; + p->srt_a.a[i].pos = (uint32_t)i; + p->srt_a.a[i].rid = i>>32; + } + radix_sort_ha_mzl_t_srt(p->srt_a.a, p->srt_a.a + p->srt_a.n); + kvec_t(uint64_t) cut; kv_init(cut); + for (k = 1, l = 0; k <= p->srt_a.n; k++) { + if(k == p->srt_a.n || p->srt_a.a[l].x != p->srt_a.a[k].x) { + for (i = l; i < k; i++) { + m = p->srt_a.a[i].rid; m <<= 32; m |= p->srt_a.a[i].pos; + assert(p->srt_a.a[i].x == p->idx_a.a[m].x); + p->srt_a.a[i] = p->idx_a.a[m]; p->idx_a.a[m].x = i; + } + kv_push(uint64_t, cut, (k - l)); + l = k; + } + } + + if(cut.n > 0) { + radix_sort_gfa64(cut.a, cut.a + cut.n); + m = cut.n * 0.0002; cc = cut.a[cut.n-1] + 1; + if(m > 0 && m <= cut.n) cc = cut.a[cut.n-m] + 1; + if(cc < (uint64_t)p->mini_cut) p->mini_cut = cc; + } + kv_destroy(cut); + // fprintf(stderr, "[M::%s::] p->mini_cut::%d \n", __func__, p->mini_cut); + + // fprintf(stderr, "[M::%s::] ==> 2\n", __func__); + p->is_cnt = 0; p->is_ovlp = 1; + memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread); + kt_for(p->n_thread, worker_for_trans_ovlp_mmhap_adv, p, p->ug->u.n); + // fprintf(stderr, "[M::%s::] ==> 3\n", __func__); + for (i = 0; (int64_t)i < p->n_thread; i++) { + ha_ovec_destroy(p->hab[i]); + free(p->ll[i].lo.a); free(p->ll[i].srt.a.a); free(p->ll[i].tc.a); + } + free(p->idx_a.a); free(p->idx_n.a); free(p->hab); free(p->free_cnt); + destory_ug_rid_cov_t(p->ccov); free(p->ccov); + // destory_bubbles(p->bub); free(p->bub); + + // fprintf(stderr, "[M::%s::] ==> 4\n", __func__); + ///make results consistent + kv_resize(ha_mzl_t, p->srt_a, p->ug->u.n); p->srt_a.n = p->ug->u.n; + for (i = 0; i < p->srt_a.n; i++) { + tz = &(p->srt_a.a[i]); + tz->x = (uint64_t)-1; tz->rev = 0; + tz->pos = tz->rid = tz->span = 0; + } + // memset(p->srt_a.a, 0, sizeof((*(p->srt_a.a)))*p->srt_a.n); + for (i = 0, occ = res->n; (int64_t)i < p->n_thread; i++) { + bl = &(p->ll[i].tk); + if(!(bl->n)) continue; + for (k = 1, l = 0; k <= bl->n; k++) { + if(k == bl->n || bl->a[k].qn != bl->a[l].qn) { + if(k > l) { + tz = &(p->srt_a.a[bl->a[l].qn]); + tz->x = bl->a[l].qn; tz->x <<= 32; tz->x |= i; + tz->rid = l>>32; tz->pos = (uint32_t)l; tz->rev = 1; + occ += (k - l); + } + l = k; + } + } + } + // fprintf(stderr, "[M::%s::] ==> 5\n", __func__); + // kt_for(p->n_thread, worker_for_sysm_trans_ovlp, p, p->ug->u.n);///not correct + // assert(p->srt_a.n <= p->ug->u.n); + // radix_sort_ha_mzl_t_srt(p->srt_a.a, p->srt_a.a + p->srt_a.n); + kv_resize(u_trans_t, *res, occ); + for (i = 0; i < p->srt_a.n; i++) { + tz = &(p->srt_a.a[i]); + if(!(tz->rev)) continue; + bl = &(p->ll[(uint32_t)(tz->x)].tk); + k = tz->rid; k <<= 32; k += tz->pos; + assert(bl->a[k].qn == (tz->x>>32)); + for (; (k < bl->n) && (bl->a[k].qn == (tz->x>>32)); k++) { + if(bl->a[k].qn == bl->a[k].tn) continue; + ww = cal_trans_ov_w(&(bl->a[k])); + if(ww <= 0) continue; + + kv_pushp(u_trans_t, *res, &z); + z->f = RC_3; z->rev = bl->a[k].rev; z->del = 0; + z->qn = bl->a[k].qn; z->qs = bl->a[k].qs; z->qe = bl->a[k].qe; + z->tn = bl->a[k].tn; z->ts = bl->a[k].ts; z->te = bl->a[k].te; + z->nw = ww; + // if(z->ts >= z->te || z->qs >= z->qe) { + // fprintf(stderr, "[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\n", __func__, + // z->qn+1, "lc"[p->ug->u.a[z->qn].circ], p->ug->u.a[z->qn].len, z->qs, z->qe, "+-"[z->rev], + // z->tn+1, "lc"[p->ug->u.a[z->tn].circ], p->ug->u.a[z->tn].len, z->ts, z->te); + // } + // if(z->qn == 56 || z->qn == 160 || z->tn == 56 || z->tn == 160) { + // fprintf(stderr, ">>>utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\tnw::%f\n", + // z->qn+1, "lc"[p->ug->u.a[z->qn].circ], p->ug->u.a[z->qn].len, z->qs, z->qe, "+-"[z->rev], + // z->tn+1, "lc"[p->ug->u.a[z->tn].circ], p->ug->u.a[z->tn].len, z->ts, z->te, z->nw); + // } + // if(z->nw <= 0) res->n--; + } + } + // fprintf(stderr, "[M::%s::] ==> 6\n", __func__); + for (i = 0; (int64_t)i < p->n_thread; i++) free(p->ll[i].tk.a); + free(p->srt_a.a); free(p->ll); + fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time); +} + void trans_base_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res, bubble_type *bub) { ug_trans_t sl; init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); init_ug_trans_t(&sl, uopt, 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, 1.0-asm_opt.trans_base_rate_sec, 1.0-asm_opt.trans_base_rate_sec, - 0.85, asm_opt.thread_num, 512, 3, 1, ug, sg, bub); + 0.85, asm_opt.thread_num, 512, 3, 1, ug, sg, bub, 1); // gen_trans_base_count(&sl, res); gen_trans_base_count_comp(&sl, res); if(!bub) { @@ -19852,6 +20323,17 @@ void trans_base_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res, } } +void trans_base_mmhap_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res) +{ + ug_trans_t sl; + init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); + init_ug_trans_t(&sl, uopt, 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, 1.0-asm_opt.trans_base_rate_sec, 1.0-asm_opt.trans_base_rate_sec, + 0.85, asm_opt.thread_num, 512, 3, 1, ug, sg, NULL, 0); + // gen_trans_base_count(&sl, res); + gen_trans_base_count_mmhap_comp(&sl, res); +} + void init_ug_bin_t(ug_bin_t *sl, const ug_opt_t *uopt, int32_t is_HPC, int32_t k, int32_t w, int32_t max_n_chain, double bw_thres, double diff_ov, double diff_bin, uint64_t max_diff, uint64_t min_bin_len, int32_t n_thread, int32_t mini_cut, int32_t chain_cut, int32_t keep_unsymm_arc, ma_ug_t *ug, asg_t *sg) diff --git a/inter.h b/inter.h index e5b74b5..648157a 100644 --- a/inter.h +++ b/inter.h @@ -126,5 +126,6 @@ uint32_t infer_se(uint32_t qs, uint32_t qe, uint32_t ts, uint32_t te, uint32_t r uint32_t rqs, uint32_t rqe, uint32_t *rts, uint32_t *rte); uint32_t clean_contain_g(const ug_opt_t *uopt, asg_t *sg, uint32_t push_trans); void dedup_contain_g(const ug_opt_t *uopt, asg_t *sg); +void trans_base_mmhap_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res); #endif