From 786f8bdff024a7290274b402fb077e7cbc5f0d40 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Fri, 2 Dec 2022 09:40:25 -0500 Subject: [PATCH] r465 --- CommandLines.h | 2 +- Overlaps.cpp | 167 +++++++++++++- gfa_ut.cpp | 21 +- gfa_ut.h | 3 + hic.cpp | 201 ++++++++++++++++- rcut.cpp | 600 ++++++++++++++++++++++++++++++++++++++++++++++++- 6 files changed, 967 insertions(+), 27 deletions(-) diff --git a/CommandLines.h b/CommandLines.h index bf89041..3760230 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.17.7-r461" +#define HA_VERSION "0.18.0-r465" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index 14a010f..e231d6d 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -7718,7 +7718,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen) kv_push(uint32_t, b_r, w); aw = asg_arc_a(g, w); - min_edge = (u_int32_t)-1; + min_edge = (uint32_t)-1; for (t = 0; t < nw; t++) { if(aw[t].del) continue; @@ -13463,6 +13463,153 @@ void hic_clean(asg_t* read_g) kv_destroy(ax); } + +void hic_clean_adv(asg_t *sg, ug_opt_t *uopt) +{ + uint32_t i, k, m, z, v, w; ma_utg_t *u = NULL; uint32_t *ba, bn, n_vtx, beg, end, n0, n1; + ma_ug_t *ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); n_vtx = ug->g->n_seq<<1; double bub_rate = 0.1; + uint8_t *bf = NULL; bubble_type *bub = gen_bubble_chain(sg, ug, uopt, &bf); + uint64_t tLen, vocc, socc, pocc; buf_t b; memset(&b, 0, sizeof(buf_t)); CALLOC(b.a, n_vtx); + REALLOC(bf, n_vtx); memset(bf, 0, sizeof((*bf))*n_vtx); + kvec_t(uint64_t) buf; kv_init(buf); n0 = n1 = 0; + + for (i = 0; i < bub->b_ug->u.n; i++) { + u = &(bub->b_ug->u.a[i]); + if(u->n == 0) continue; + for (k = 0; k < u->n; k++) {///bubble chain + get_bubbles(bub, u->a[k]>>33, &beg, &end, &ba, &bn, NULL);///bubble + for (m = vocc = tLen = 0; m < bn; m++) { + bf[ba[m]] = bf[ba[m]^1] = 1; + tLen += ug->u.a[ba[m]>>1].len; + if(IF_HOM((ba[m]>>1), *bub)) continue; + if(ug->g->seq[ba[m]>>1].del) continue; + vocc += ug->u.a[ba[m]>>1].n; + } + if(beg != (uint32_t)-1) tLen += ug->u.a[beg>>1].len; + if(end != (uint32_t)-1) tLen += ug->u.a[end>>1].len; + + if(vocc) { + for (m = buf.n = 0; m < bn; m++) { + v = ba[m]; + if(ug->g->seq[v>>1].del) continue; + if(asg_arc_n(ug->g, v) < 2) continue; + if(get_real_length(ug->g, v, NULL) < 2) continue; + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) { + //beg is v, end is b.S.a[0] + //note b.b include end, does not include beg + for (z = socc = 0; z < b.b.n; z++) { + if(b.b.a[z]==v || b.b.a[z]==b.S.a[0]) continue; + socc += ug->u.a[b.b.a[z]>>1].n; + if((!bf[b.b.a[z]])&&(!bf[b.b.a[z]^1])) break; + } + if(z < b.b.n) continue; + kv_push(uint64_t, buf, ((socc<<32)|v)); + } + + v ^= 1; + if(asg_arc_n(ug->g, v) < 2) continue; + if(get_real_length(ug->g, v, NULL) < 2) continue; + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) { + //beg is v, end is b.S.a[0] + //note b.b include end, does not include beg + for (z = socc = 0; z < b.b.n; z++) { + if(b.b.a[z]==v || b.b.a[z]==b.S.a[0]) continue; + socc += ug->u.a[b.b.a[z]>>1].n; + if((!bf[b.b.a[z]])&&(!bf[b.b.a[z]^1])) break; + } + if(z < b.b.n) continue; + kv_push(uint64_t, buf, ((socc<<32)|v)); + } + } + + radix_sort_arch64(buf.a, buf.a + buf.n); + for (m = pocc = 0; m < buf.n; m++) { + v = (uint32_t)buf.a[m]; + if(ug->g->seq[v>>1].del) continue; + if(asg_arc_n(ug->g, v) < 2) continue; + if(get_real_length(ug->g, v, NULL) < 2) continue; + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) { + //beg is v, end is b.S.a[0] + //note b.b include end, does not include beg + for (z = socc = 0; z < b.b.n; z++) { + if(b.b.a[z]==v || b.b.a[z]==b.S.a[0]) continue; + socc += ug->u.a[b.b.a[z]>>1].n; + if((!bf[b.b.a[z]])&&(!bf[b.b.a[z]^1])) break; + } + if(z < b.b.n) continue; + if((pocc+socc) >= (vocc*bub_rate)) continue; + pocc += socc; + asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0, NULL); + // fprintf(stderr, "+utg%.6dl\tutg%.6dl\n", (int32_t)(v>>1)+1, (int32_t)(b.S.a[0]>>1)+1); + n0++; + } + } + } + for (m = 0; m < bn; m++) { + bf[ba[m]] = bf[ba[m]^1] = 0; + } + } + } + + tLen = get_bub_pop_max_dist_advance(ug->g, &b); + for (v = buf.n = 0; v < n_vtx; ++v) { + if(ug->g->seq[v>>1].del) continue; + if(asg_arc_n(ug->g, v) < 2) continue; + if(get_real_length(ug->g, v, NULL) < 2) continue; + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) { + //beg is v, end is b.S.a[0] + //note b.b include end, does not include beg + for (i = socc = 0; i < b.b.n; i++) { + if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue; + socc += ug->u.a[b.b.a[i]>>1].n; + } + if(socc <= 16) kv_push(uint64_t, buf, ((socc<<32)|v)); + } + } + + uint32_t convex; long long ll, tmp, max_stop_nodeLen, max_stop_baseLen; + radix_sort_arch64(buf.a, buf.a + buf.n); + for (m = 0; m < buf.n; m++) { + v = (uint32_t)buf.a[m]; + if(ug->g->seq[v>>1].del) continue; + if(asg_arc_n(ug->g, v) < 2) continue; + if(get_real_length(ug->g, v, NULL) < 2) continue; + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) { + w = b.S.a[0]^1; + for (i = socc = 0; i < b.b.n; i++) { + if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue; + socc += ug->u.a[b.b.a[i]>>1].n; + } + if(socc <= 16) { + b.b.n = 0; + get_unitig(ug->g, NULL, v^1, &convex, &ll, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 1, &b); + for (k = vocc = 0; k < b.b.n; k++) { + if(IF_HOM((b.b.a[k]>>1), *bub)) break; + vocc += ug->u.a[b.b.a[k]>>1].n; + } + if((socc) >= (vocc*bub_rate)) continue; + + b.b.n = 0; + get_unitig(ug->g, NULL, w^1, &convex, &ll, &tmp, &max_stop_nodeLen, &max_stop_baseLen, 1, &b); + for (k = vocc = 0; k < b.b.n; k++) { + if(IF_HOM((b.b.a[k]>>1), *bub)) break; + vocc += ug->u.a[b.b.a[k]>>1].n; + } + if((socc) >= (vocc*bub_rate)) continue; + asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0, NULL); + // fprintf(stderr, "-utg%.6dl\tutg%.6dl\n", (int32_t)(v>>1)+1, (int32_t)(b.S.a[0]>>1)+1); + n1++; + } + } + } + filter_sg_by_ug(sg, ug, uopt); + + free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); + ma_ug_destroy(ug); free(bf); kv_destroy(buf); + destory_bubbles(bub); free(bub); + // fprintf(stderr, "[M::%s::] # type0::%u, # type1::%u\n", __func__, n0, n1); +} + void update_dump_trio(uint8_t* trio_flag, uint32_t rn, uint8_t *rf, ma_ug_t *ug) { uint32_t k, i, x; @@ -18387,7 +18534,7 @@ void chain_origin_trans_uid_s_bubble(buf_t *pri, buf_t* aux, uint32_t beg, uint3 if(av[i].v == pri_v) priEnd = ((pri_len > av[i].ol)? (pri_len - av[i].ol - 1) : 0); if(av[i].v == aux_v) auxEnd = ((aux_len > av[i].ol)? (aux_len - av[i].ol - 1) : 0); } - + ///[priBeg, priEnd) && [auxBeg, auxEnd) if(priBeg == (uint32_t)-1 || priEnd == (uint32_t)-1 || auxBeg == (uint32_t)-1 || auxEnd == (uint32_t)-1) { fprintf(stderr, "ERROR-s_bubble\n"); @@ -18508,7 +18655,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha for (k = 0; k < p->n; k++) { rId = p->a[k]>>33; - cov->cov[rId] += (uCov * cov->read_g->seq[rId].len); + cov->cov[rId] += (uCov * cov->read_g->seq[rId].len);///this the average coverage of the whole bubble } } v = u; @@ -18517,6 +18664,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha if(t_ch) { + ///is this requirement appropriate? if(get_real_length(ug->g, v0, NULL) == 2 && get_real_length(ug->g, b->S.a[0]^1, NULL) == 2) { long long tmp, max_stop_nodeLen, max_stop_baseLen, bch_occ[2]; @@ -18530,16 +18678,17 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha &max_stop_nodeLen, &max_stop_baseLen, 1, NULL); get_unitig(ug->g, NULL, bch[1], &convex[1], &bch_occ[1], &tmp, &max_stop_nodeLen, &max_stop_baseLen, 1, NULL); + ///if this is a simple bubble if(((bch_occ[0] + bch_occ[1] + 1) == (uint32_t)b->b.n) && get_real_length(ug->g, convex[0], NULL) == 1 && get_real_length(ug->g, convex[1], NULL) == 1) { get_real_length(ug->g, convex[0], &convex[0]); get_real_length(ug->g, convex[1], &convex[1]); - if(convex[0] == b->S.a[0] && convex[1] == b->S.a[0]) + if(convex[0] == b->S.a[0] && convex[1] == b->S.a[0])///double check if it is a simple bubble { t_ch->b_buf_0.b.n = 0; get_unitig(ug->g, NULL, bch[0], &convex[0], &bch_occ[0], &tmp, &max_stop_nodeLen, &max_stop_baseLen, 1, &(t_ch->b_buf_0)); - for (i = 0; i < t_ch->b_buf_0.b.n; ++i) + for (i = 0; i < t_ch->b_buf_0.b.n; ++i)///retrive one side of the bubble { uId = t_ch->b_buf_0.b.a[i]>>1; p = &(ug->u.a[uId]); @@ -18553,7 +18702,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha t_ch->b_buf_1.b.n = 0; get_unitig(ug->g, NULL, bch[1], &convex[1], &bch_occ[1], &tmp, &max_stop_nodeLen, &max_stop_baseLen, 1, &(t_ch->b_buf_1)); - for (i = 0; i < t_ch->b_buf_1.b.n; ++i) + for (i = 0; i < t_ch->b_buf_1.b.n; ++i)///retrive another side of the bubble { uId = t_ch->b_buf_1.b.a[i]>>1; p = &(ug->u.a[uId]); @@ -18564,7 +18713,7 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha t_ch->ir_het[(ori == 1?((p->a[p->n-k-1])>>33):(p->a[k]>>33))] |= P_HET; } } - + ///generate read-to-read overlaps chain_origin_trans_uid_s_bubble(&(t_ch->b_buf_0), &(t_ch->b_buf_1), v0, b->S.a[0]^1, ug, cov); return; @@ -31963,6 +32112,10 @@ ma_sub_t **coverage_cut_ptr, int debug_g) // flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL; flat_soma_v(sg, sources, ruIndex); **/ + if(!(ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt))) { + // output_unitig_graph(sg, coverage_cut, "pre_clean", sources, ruIndex, max_hang_length, mini_overlap_length); + hic_clean_adv(sg, &uopt); + } output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources, asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length); diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 3944763..c475c77 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -13127,15 +13127,17 @@ asg64_v *b0, asg64_v *b1, double cutoff) for (z = 0; z < zn && b0->a[z] == b1->a[z]; z++); ///z: first raw unitig that is different between two paths assert(z > 0); w0 = w1 = (uint64_t)-1; if(z < b0->n) get_integer_seq_ovlps(uidx, b0->a, b0->n, z - 1, 0, NULL, &w0); + else return 0;///b0 is contained if(z < b1->n) get_integer_seq_ovlps(uidx, b1->a, b1->n, z - 1, 0, NULL, &w1); - // if(is_debug) { + else continue;///b1 is contained + // if(((v>>1) == 8722) || ((v>>1) == 56768)) { // prt_sub_integer_path(b0, it, w0, w1, zn, z); // prt_sub_integer_path(b1, it, w0, w1, zn, z); // } if(w0 == (uint64_t)-1) w0 = 0; if(w1 == (uint64_t)-1) w1 = 0; - if(b0->n == zn) return 0;///b0 is shorter - if(b1->n == zn) continue;///b1 is shorter + // if(b0->n == zn) return 0;///b0 is shorter + // if(b1->n == zn) continue;///b1 is shorter if((min_w0 == (uint64_t)-1) || (z == zn) || (min_w0 > w0) || (min_w0 == w0 && min_w1 < w1)) { min_w0 = w0; min_w1 = w1; } @@ -13163,9 +13165,17 @@ uint64_t *ridx, asg64_v *res) v = int_a[k]; if((!f[v])&&(!f[v^1])) continue; if((pi != (uint64_t)-1) && (f[v^1])) { - // fprintf(stderr, "+[M::%s::] utg%.6dl(%c), f[v^1]::%u ,v^1::%lu\n", - // __func__, (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1], f[v^1], v^1); + // if(((v>>1) == 8722) || ((v>>1) == 56768)) { + // fprintf(stderr, "+[M::%s::ii[%lu, %lu)] utg%.6dl(%c), f[v^1]::%u, v^1::%lu, putg%.6dl(%c), f[pv]::%u, pv::%lu\n", + // __func__, s, e, (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1], f[v^1], v^1, + // (int32_t)(int_a[pi]>>1)+1, "+-"[int_a[pi]&1], f[int_a[pi]], int_a[pi]); + // } if(is_best_path(uidx, ng, int_idx, int_a, s, e, k, v^1, ridx_a, ridx, b0, b1, 0.51)) { + // if(((v>>1) == 8722) || ((v>>1) == 56768)) { + // fprintf(stderr, "-[M::%s::ii[%lu, %lu)] utg%.6dl(%c), f[v^1]::%u, v^1::%lu, putg%.6dl(%c), f[pv]::%u, pv::%lu\n", + // __func__, s, e, (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1], f[v^1], v^1, + // (int32_t)(int_a[pi]>>1)+1, "+-"[int_a[pi]&1], f[int_a[pi]], int_a[pi]); + // } // fprintf(stderr, "[M::%s::] utg%.6dl(%c)->utg%.6dl(%c)\n", __func__, // (int32_t)(int_a[pi]>>1)+1, "+-"[int_a[pi]&1], (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1]); pz = (res->n > res_n)? &(res->a[res->n-1]):(NULL); @@ -16361,4 +16371,5 @@ ul_renew_t *ropt) ma_hit_contained_advance((*(ropt->src)), (*(ropt->n_read)), (*(ropt->cov)), ropt->ruIndex, ropt->max_hang, ropt->mini_ovlp); post_rescue(uopt, (*(ropt->sg)), (*(ropt->src)), (*(ropt->r_src)), ropt->ruIndex, ropt->b_mask_t, 0); // print_raw_uls_aln(uidx, asm_opt.output_file_name); + // exit(0); } \ No newline at end of file diff --git a/gfa_ut.h b/gfa_ut.h index a2dfb5e..8fc2274 100644 --- a/gfa_ut.h +++ b/gfa_ut.h @@ -1,6 +1,7 @@ #ifndef __GFA_UT__ #define __GFA_UT__ #include "Overlaps.h" +#include "hic.h" typedef struct { asg_t *g; @@ -31,5 +32,7 @@ void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd); asg_t *gen_ng(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **cov, R_to_U *ruI, uint64_t scaffold_len); void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover); // void print_raw_u2rgfa_seq(all_ul_t *aln, R_to_U* rI, uint32_t is_detail); +bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het); +void filter_sg_by_ug(asg_t *rg, ma_ug_t *ug, ug_opt_t *uopt); #endif diff --git a/hic.cpp b/hic.cpp index 33bb101..ee3e875 100644 --- a/hic.cpp +++ b/hic.cpp @@ -15,6 +15,7 @@ #include "kseq.h" // FASTA/Q parser #include "kdq.h" #include "horder.h" +#include "gfa_ut.h" KSEQ_INIT(gzFile, gzread) KDQ_INIT(uint64_t) @@ -16481,15 +16482,14 @@ void optimize_u_trans(kv_u_trans_t *ovlp, kvec_pe_hit* hits, ha_ug_index* idx) x = &(ovlp->a[i]); if(x->qn > x->tn) continue; if(x->f != RC_2 || x->del) continue; - occ = get_oe_occ(x->qn, x->tn, hits, idx) + get_oe_occ(x->tn, x->qn, hits, idx); + occ = get_oe_occ(x->qn, x->tn, hits, idx) + get_oe_occ(x->tn, x->qn, hits, idx);///how many UL bridging qn and tn kv_pushp(u_trans_t, k_trans, &p); (*p) = (*x); p->nw = (x->nw*(1-(((double)(occ<<1))/((double)(hits->occ.a[x->qn]+hits->occ.a[x->tn]))))); if(p->nw < 0) fprintf(stderr, "ERROR-nw\n"); if(p->nw == 0) p->nw = x->nw*0.005; if(p->nw == 0) { k_trans.n--; - } - else { + } else { kv_pushp(u_trans_t, k_trans, &p); (*p) = k_trans.a[k_trans.n-2]; p->qn = k_trans.a[k_trans.n-2].tn; p->qs = k_trans.a[k_trans.n-2].ts; p->qe = k_trans.a[k_trans.n-2].te; @@ -16540,6 +16540,179 @@ ha_ug_index* idx, uint64_t step, uint64_t total) } +void prt_hits_noid(ha_ug_index* idx, ma_ug_t* ug, kvec_pe_hit* hits, FILE *fn) +{ + uint64_t k, shif = 64 - idx->uID_bits; + char dir[2] = {'+', '-'}; + for (k = 0; k < hits->a.n; ++k) { + fprintf(fn, "r-%lu-th\t%c\trs-utg%.6d%c\t%lu\t%c\tre-utg%.6d%c\t%lu\n", + hits->a.a[k].id, + dir[hits->a.a[k].s>>63], (int)((hits->a.a[k].s<<1)>>shif)+1, + "lc"[ug->u.a[((hits->a.a[k].s<<1)>>shif)].circ], hits->a.a[k].s&idx->pos_mode, + dir[hits->a.a[k].e>>63], (int)((hits->a.a[k].e<<1)>>shif)+1, + "lc"[ug->u.a[((hits->a.a[k].e<<1)>>shif)].circ], hits->a.a[k].e&idx->pos_mode); + } +} + +void prt_utg_trans(kv_u_trans_t *ta, ma_ug_t* ug, FILE *fn) +{ + uint32_t i; + u_trans_t *p = NULL; + for (i = 0; i < ta->n; i++) { + p = &(ta->a[i]); + fprintf(fn, "utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\tw(%f)\tf(%u)\n", + p->qn+1, "lc"[ug->u.a[p->qn].circ], ug->u.a[p->qn].len, p->qs, p->qe, "+-"[p->rev], + p->tn+1, "lc"[ug->u.a[p->tn].circ], ug->u.a[p->tn].len, p->ts, p->te, p->nw, p->f); + } +} + + +void prt_kv_u_trans(kv_u_trans_t *ta, hc_links* lk, int8_t *s, FILE *fn) +{ + uint32_t i; + u_trans_t *p = NULL; + hc_edge *e = NULL; + + for (i = 0; i < ta->n; i++) { + p = &(ta->a[i]); + e = get_hc_edge(lk, p->qn, p->tn, 0); + fprintf(fn, "s-utg%.6ul\tS(%d)\td-utg%.6ul\tS(%d)\trev(%u)\td(%lld)\ttw(%f)\n", + p->qn+1, s[p->qn], p->tn+1, s[p->tn], p->rev, + (e == NULL || e->dis == (uint64_t)-1)? -1 : (long long)(e->dis>>3), p->nw); + } +} + + +void prt_bubble_gfa_adv(FILE *fp, bubble_type *bub, const char* utg_pre, const char* bub_pre, const char* chain_pre) +{ + uint32_t i, k, m, *a, n, beg, sink, x; ma_utg_t *p; uint64_t occ; + ma_ug_t *b_ug = bub->b_ug; char name[32], bname[32]; uint8_t *f; CALLOC(f, bub->ug->u.n); + for (i = 0; i < b_ug->u.n; i++) { + p = &b_ug->u.a[i]; + if(p->n == 0) continue; + for (k = occ = 0; k < p->n; k++){ + x = p->a[k]>>33; + get_bubbles(bub, x, &beg, &sink, &a, &n, NULL); + + for (m = 0; m < n; m++) { + occ += bub->ug->u.a[a[m]>>1].n; f[a[m]>>1] = 1; + } + if(beg != (uint32_t)-1 && f[beg>>1] == 0) { + occ += bub->ug->u.a[beg>>1].n; f[beg>>1] = 1; + } + if(sink != (uint32_t)-1 && f[sink>>1] == 0) { + occ += bub->ug->u.a[sink>>1].n; f[sink>>1] = 1; + } + } + + sprintf(name, "%s%.6d%c", chain_pre, i + 1, "lc"[p->circ]); + fprintf(fp, "S\t%s\t*\tLN:i:%lu\n", name, occ); + for (k = 0; k < p->n; k++) { + x = p->a[k]>>33; + sprintf(bname, "%s%.6d", bub_pre, x + 1); + fprintf(fp, "B\t%s\t%c\tcid:i:%s\tsm:%c\n", bname, "+-"[(p->a[k]>>32)&1], name, "01"[xf_bub]); + + get_bubbles(bub, x, &beg, &sink, &a, &n, NULL); + if(beg != (uint32_t)-1) { + fprintf(fp, "U\t%s%.6d%c\t%c\tcid:i:%s\tbid:b:%s\thom:%c\n", + utg_pre, (beg>>1)+1, "lc"[bub->ug->u.a[(beg>>1)].circ], "+-"[beg&1], name, bname, "10"[IF_HOM((beg>>1), *bub)]); + } + + if(sink != (uint32_t)-1) { + fprintf(fp, "U\t%s%.6d%c\t%c\tcid:i:%s\tbid:s:%s\thom:%c\n", + utg_pre, (sink>>1)+1, "lc"[bub->ug->u.a[(sink>>1)].circ], "+-"[sink&1], name, bname, "10"[IF_HOM((sink>>1), *bub)]); + } + for (m = 0; m < n; m++) { + occ += bub->ug->u.a[a[m]>>1].n; + fprintf(fp, "U\t%s%.6d%c\t%c\tcid:i:%s\tbid:c:%s\thom:%c\n", + utg_pre, (a[m]>>1)+1, "lc"[bub->ug->u.a[(a[m]>>1)].circ], "+-"[a[m]&1], name, bname, "10"[IF_HOM((a[m]>>1), *bub)]); + } + } + } + + asg_arc_t* au = NULL; + uint32_t nu, u, v, j; + for (i = 0; i < b_ug->u.n; ++i) { + if(b_ug->u.a[i].m == 0) continue; + if(b_ug->u.a[i].circ) + { + fprintf(fp, "L\t%s%.6dc\t+\t%s%.6dc\t+\t%dM\tL1:i:%d\n", + chain_pre, i+1, chain_pre, i+1, 0, 0); + fprintf(fp, "L\t%s%.6dc\t-\t%s%.6dc\t-\t%dM\tL1:i:%d\n", + chain_pre, i+1, chain_pre, i+1, 0, 0); + } + u = i<<1; + au = asg_arc_a(b_ug->g, u); + nu = asg_arc_n(b_ug->g, u); + for (j = 0; j < nu; j++) + { + if(au[j].del) continue; + v = au[j].v; + fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\n", + chain_pre, (u>>1)+1, "lc"[b_ug->u.a[u>>1].circ], "+-"[u&1], + chain_pre, (v>>1)+1, "lc"[b_ug->u.a[v>>1].circ], "+-"[v&1], 0, 0); + } + + + u = (i<<1) + 1; + au = asg_arc_a(b_ug->g, u); + nu = asg_arc_n(b_ug->g, u); + for (j = 0; j < nu; j++) + { + if(au[j].del) continue; + v = au[j].v; + fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\n", + chain_pre, (u>>1)+1, "lc"[b_ug->u.a[u>>1].circ], "+-"[u&1], + chain_pre, (v>>1)+1, "lc"[b_ug->u.a[v>>1].circ], "+-"[v&1], 0, 0); + } + } + + for (i = 0; i < bub->ug->u.n; i++) { + if(f[i]) continue; + fprintf(fp, "U\t%s%.6d%c\t+\tcid:i:*\tbid:c:*\thom:%c\n", + utg_pre, i+1, "lc"[bub->ug->u.a[i].circ], "10"[IF_HOM(i, *bub)]); + } + free(f); +} + + +void prt_debug_hic(const char* o_n, ma_ug_t* ug, ha_ug_index* idx, ug_opt_t *opt, kvec_pe_hit* hits, +kv_u_trans_t *utg_trans, kv_u_trans_t *p_arcs, hc_links* lk, int8_t *s, bubble_type *bub) +{ + char* gfa_name = (char*)malloc(strlen(o_n)+100); FILE *fn = NULL; + + sprintf(gfa_name, "%s.hic.dbg", o_n); + print_debug_gfa(idx->read_g, idx->ug, opt->coverage_cut, gfa_name, opt->sources, opt->ruIndex, + opt->max_hang, opt->min_ovlp, 0, 0, 0); + + if(hits) { + sprintf(gfa_name, "%s.hic.hits.log", o_n); fn = fopen(gfa_name, "w"); + prt_hits_noid(idx, ug, hits, fn); + fclose(fn); + } + + if(utg_trans) { + sprintf(gfa_name, "%s.hic.utg.trans.log", o_n); fn = fopen(gfa_name, "w"); + prt_utg_trans(utg_trans, ug, fn); + fclose(fn); + } + + if(p_arcs) { + sprintf(gfa_name, "%s.hic.parcs.log", o_n); fn = fopen(gfa_name, "w"); + prt_kv_u_trans(p_arcs, lk, s, fn); + fclose(fn); + } + + if(bub) { + sprintf(gfa_name, "%s.bub.noseq.gfa", o_n); fn = fopen(gfa_name, "w"); + prt_bubble_gfa_adv(fn, bub, "utg", "btg", "ctg"); + fclose(fn); + } + + free(gfa_name); + fprintf(stderr, "[M::%s::] done\n", __func__); +} + int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_opt_t *opt, kvec_pe_hit **rhits) { double index_time = yak_realtime(); @@ -16608,11 +16781,15 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o // if(bub.round_id == 0) init_phase(idx, &k_trans, &bub, s); // update_trans_g(idx, &k_trans, &bub); /*******************************for debug************************************/ - // mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, - // (bub.round_id == 0? 1 : 0), s->s, 1, (asm_opt.ar)?(&bub):(NULL), &(idx->t_ch->k_trans), 0, - // (((bub.round_id+1) == bub.n_round)?1:0)); mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, - (bub.round_id == 0? 1 : 0), s->s, 1, NULL, &(idx->t_ch->k_trans), 0, 0); + (bub.round_id == 0? 1 : 0), s->s, 1, &bub, + &(idx->t_ch->k_trans), 0, 0/**(((bub.round_id+1) == bub.n_round)?1:0)**/); + // mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, + // (bub.round_id == 0? 1 : 0), s->s, 1, NULL, &(idx->t_ch->k_trans), 0, 0); + // if((bub.round_id+1) == bub.n_round) { + // prt_debug_hic(asm_opt.output_file_name, idx->ug, idx, opt, &sl.hits, + // &(idx->t_ch->k_trans), &k_trans, &link, s->s, &(bub)); + // } /*******************************for debug************************************/ label_unitigs_sm(s->s, NULL, idx->ug); @@ -16642,15 +16819,15 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o // horder_t *ho = init_horder_t(&sl.hits, idx->uID_bits, idx->pos_mode, idx->read_g, idx->ug, &bub, &(idx->t_ch->k_trans), opt, 3); - ///print_hc_links(&link, 0, &hap); + // print_hc_links(&link, 0, &hap); + // print_hits_simp(idx, &sl.hits); + // print_kv_u_trans_t(&(idx->t_ch->k_trans)); // print_kv_u_trans(&k_trans, &link, s->s); - // print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, NULL/**idx->link**/, idx); // print_hits(idx, &sl.hits, fn1, fn2); - - - ///print_debug_bubble_graph(&bub, idx->ug, asm_opt.output_file_name); + // print_debug_bubble_graph(&bub, idx->ug, asm_opt.output_file_name); // print_bubble_chain(&bub); + // destory_contig_partition(&hap); // destory_horder_t(&ho); kv_destroy(sl.hits.a); kv_destroy(sl.hits.idx); kv_destroy(sl.hits.occ); diff --git a/rcut.cpp b/rcut.cpp index 0746c73..0b716a8 100644 --- a/rcut.cpp +++ b/rcut.cpp @@ -1,6 +1,8 @@ #define __STDC_LIMIT_MACROS #include #include +#include +#include "assert.h" #include "rcut.h" #include "Purge_Dups.h" #include "Correct.h" @@ -94,6 +96,28 @@ typedef struct{ bits_p *vis; }mc_bp_t; +typedef struct{ + kvec_t(uint32_t) nn; + kvec_t(uint64_t) ng; +} nn_clus_t; + +typedef struct{ + bits_p vis; + t_w_t w; + uint32_t off, occ; +}clus_flip_aux; + +typedef struct{ + nn_clus_t cc; + bubble_type* bub; + const mc_opt_t *opt; + mc_g_t *mg; + uint8_t *lock, lock_max, dbg; + uint32_t n, n_thread; + clus_flip_aux *aux; + mc_svaux_t *baux; +} mc_clus_t; + typedef struct { uint64_t x; // RNG uint32_t cc_off, cc_size; @@ -1685,6 +1709,7 @@ static t_w_t mc_optimize_local(const mc_opt_t *opt, const mc_match_t *ma, mc_sva if (b->z[k].z[0] == b->z[k].z[1]) continue; s = b->z[k].z[0] > b->z[k].z[1]? -1 : 1; if (b->s[k] != s) { + // fprintf(stderr, "utg%.6dl, s[k]::%d, s::%d\n", (int32_t)(k)+1, b->s[k], s); mc_set_spin(ma, b, k, s);///no need to change the score of k itself ++n_flip; } @@ -2209,6 +2234,517 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3 return n_iter; } +#define clus_a(b, id) (((b)).cc.nn.a+(((b)).cc.ng.a[(id)]>>32)) +#define clus_n(b, id) (((uint32_t)((((b)).cc.ng.a[(id)])))) + + +void reorder_bub(const mc_match_t *ma, uint32_t *a, uint32_t a_n, uint8_t *ff, uint8_t lf, uint8_t rf, uint8_t cf, +uint8_t cuf, double *sc_l, double *sc_r, double *sc_m, uint32_t *res) +{ + uint32_t k, o, j, n, t, z; mc_edge_t *e; + double mm_l, mm_r, m_inner, mm; int64_t mmlk, mmrk, rev, m_inner_k, mmk; + for (k = 0; k < a_n; k++) ff[a[k]] = cf; + m_inner = 0; m_inner_k = -1; mm_l = mm_r = 0; mmlk = mmrk = -1; + for (k = 0, rev = -1; k < a_n; k++) { + sc_l[k] = sc_r[k] = sc_m[k] = 0; + o = ma->idx.a[a[k]] >> 32; n = (uint32_t)ma->idx.a[a[k]]; + for (j = 0; j < n; ++j) { + e = &ma->ma.a[o + j]; + t = ma_y(*e); + if(ff[t] == lf) sc_l[k] += fabs(e->w); + if(ff[t] == rf) sc_r[k] += fabs(e->w); + if(ff[t] == cf) sc_m[k] += fabs(e->w); + } + if(sc_l[k] > 0) { + if(sc_l[k] > mm_l) { + mm_l = sc_l[k]; mmlk = k; + } + } else if(sc_r[k] > 0) { + if(sc_r[k] > mm_r) { + mm_r = sc_r[k]; mmrk = k; + } + } else if(sc_m[k] > 0 && sc_m[k] > m_inner) { + m_inner = sc_m[k]; m_inner_k = k; + } + } + if(mmlk == -1 && mmrk == -1 && m_inner_k == -1) return; + if(mmlk != -1) { + mm = mm_l; mmk = mmlk; rev = 0; + } else if(mmrk != -1) { + mm = mm_r; mmk = mmrk; rev = 1; + } else { + mm = m_inner; mmk = m_inner_k; rev = 0; + } + + double *sc = (!rev)?sc_l:sc_r; + // for (k = 0; k < a_n; k++) { + // // sc_m[k] = 0; + // fprintf(stderr, "+[M::%s] a_n::%u, a[%u]::%u\n", __func__, a_n, k, a[k]); + // } + + z = 0; res[z++] = a[mmk]; ff[a[mmk]] = cuf; + for (; z < a_n; ) { + for (k = 0, mm = 0, mmk = -1; k < a_n; k++) { + if(ff[a[k]] == cuf) continue; + o = ma->idx.a[a[k]] >> 32; n = (uint32_t)ma->idx.a[a[k]]; + for (j = 0; j < n; ++j) { + e = &ma->ma.a[o + j]; + t = ma_y(*e); + if(ff[t] == cuf) sc[k] += fabs(e->w); + // if(ff[t] == cuf) sc_m[k] += fabs(e->w); + } + + if(sc[k] > 0) { + if(sc[k] > mm) { + mm = sc[k]; mmk = k; + } + } + } + if(mmk == -1) break; + // fprintf(stderr, "-[M::%s] mmk::%ld, a[%ld]::%u, z::%u\n", __func__, mmk, mmk, a[mmk], z); + res[z++] = a[mmk]; ff[a[mmk]] = cuf; + } + // fprintf(stderr, "[M::%s] a_n::%u, z::%u\n", __func__, a_n, k, z); + if(z < a_n) { + for (k = 0; k < a_n; k++) { + if(ff[a[k]] == cuf) continue; + res[z++] = a[k]; ff[a[k]] = cuf; + } + } + assert(z == a_n); + if(!rev) { + for (k = 0; k < a_n; k++) { + a[k] = res[k]; + // fprintf(stderr, "[M::%s] a_n::%u, a[k]::%u, res[k]::%u\n", __func__, a_n, a[k], res[k]); + ff[a[k]] = lf; + } + } else { + for (k = 0; k < a_n; k++) { + a[k] = res[a_n-k-1]; ff[a[k]] = lf; + } + } +} + +void prt_bub(uint32_t *a, uint32_t a_n, const char *cmd) +{ + uint32_t k; + fprintf(stderr, "%s\n", cmd); + for (k = 0; k < a_n; k++) { + fprintf(stderr, "utg%.6dl\t", (int32_t)a[k]+1); + } + fprintf(stderr, "\n"); + +} + +void renew_mc_clus_t(mc_clus_t *bc, uint32_t *a, uint32_t a_n) +{ + if(!bc) return; + // fprintf(stderr, "[M::%s] a_n::%u\n", __func__, a_n); + uint32_t k, i, *ba, bn, m, cocc, iin, bub_occ = 0, bbn = 0; uint64_t *p; ma_utg_t *u = NULL; + kvec_t(double) sc_l; kvec_t(double) sc_r; kvec_t(double) sc_m; kvec_t(uint32_t) tmp; + kv_init(sc_l); kv_init(sc_r); kv_init(sc_m); kv_init(tmp); + + bc->cc.ng.n = bc->cc.nn.n = 0; + memset(bc->lock, 0, sizeof((*(bc->lock)))*bc->n); + for (k = 0; k < a_n; k++) bc->lock[a[k]] = 1; + + kv_resize(uint32_t, bc->cc.nn, a_n); + // for (i = 0; i < bc->bub->chain_weight.n; i++) { + // if(bc->bub->chain_weight.a[i].del) continue; + // u = &(bc->bub->b_ug->u.a[bc->bub->chain_weight.a[i].id]);///list of bubbles + for (i = 0; i < bc->bub->b_ug->u.n; i++) { + u = &(bc->bub->b_ug->u.a[i]); + if(u->n == 0) continue; + bub_occ += u->n; + // fprintf(stderr, "[M::%s] i::%u, u->n::%u\n", __func__, i, (uint32_t)u->n); + for (k = cocc = 0, iin = bc->cc.nn.n; k < u->n; k++) { + get_bubbles(bc->bub, u->a[k]>>33, NULL, NULL, &ba, &bn, NULL); + kv_pushp(uint64_t, bc->cc.ng, &p); bbn += bn; + *p = bc->cc.nn.n;///a bubble + for (m = 0; m < bn; m++) { + if(!(bc->lock[ba[m]>>1])) continue; + kv_push(uint32_t, bc->cc.nn, (ba[m]>>1)); + bc->lock[ba[m]>>1] = 2; + } + if(bc->cc.nn.n <= (*p)) {///no node in this bubble + bc->cc.ng.n--; + continue; + } + *p <<= 32; *p |= (bc->cc.nn.n-((*p)>>32)); cocc++; + } + //split chains + if(cocc > 0) {//cocc: # of bubbles in this chain + for (k = bc->cc.ng.n - cocc; k < bc->cc.ng.n; k++) { + kv_resize(double, sc_l, clus_n((*bc), k)); + kv_resize(double, sc_r, clus_n((*bc), k)); + kv_resize(double, sc_m, clus_n((*bc), k)); + kv_resize(uint32_t, tmp, clus_n((*bc), k)); + // prt_bub(clus_a((*bc), k), clus_n((*bc), k), "-0-"); + reorder_bub(bc->mg->e, clus_a((*bc), k), clus_n((*bc), k), bc->lock, 3, 2, 4, 5, + sc_l.a, sc_r.a, sc_m.a, tmp.a); + // prt_bub(clus_a((*bc), k), clus_n((*bc), k), "-1-"); + } + for (k = iin; k < bc->cc.nn.n; k++) bc->lock[bc->cc.nn.a[k]] = 1;//reset + kv_pushp(uint64_t, bc->cc.ng, &p); *p = (uint64_t)-1; + kv_push(uint32_t, bc->cc.nn, ((uint32_t)-1)); ///split + } + } + + for (k = 0; k < a_n; k++) bc->lock[a[k]] = 0; + kv_destroy(sc_l); kv_destroy(sc_r); kv_destroy(sc_m); kv_destroy(tmp); + // fprintf(stderr, "[M::%s] a_n::%u, bc->cc.nn.n::%u, bc->cc.ng.n::%u, bbn::%u, bub_occ::%u, bc->bub->b_ug->u.n::%u\n", __func__, + // a_n, (uint32_t)bc->cc.nn.n, (uint32_t)bc->cc.ng.n, bbn, bub_occ, (uint32_t)bc->bub->b_ug->u.n); +} + +void clean_clus_flip_aux(clus_flip_aux *z) +{ + z->w = -1; z->occ = z->off = (uint32_t)-1; + memset(z->vis.a, 0, sizeof(*(z->vis.a))*z->vis.n); +} + +#define is_set_bits_p(v, i) (((v).a[((i)>>3)]>>(i&7))&1) +#define set_bits_p(v, i) (((v).a[((i)>>3)])|=(((uint8_t)1)<<(i&7))); + +t_w_t clus_weight(mc_svaux_t *b_aux, mc_match_t *ma, bits_p *vis, uint32_t uid) +{ + mc_edge_t *o = NULL; + uint32_t n, i, t; + t_w_t w = ((t_w_t)(b_aux->s[uid])) * (b_aux->z[uid].z[0] - b_aux->z[uid].z[1]) * 2; + t_w_t w_off = 0; + o = pt_a(*ma, uid); + n = pt_n(*ma, uid); + for (i = 0; i < n; ++i) { + t = ma_y(o[i]); + if(!(is_set_bits_p((*vis), t))) continue; + // if(vis[t] == 0) continue; + if(t == uid) continue; + w_off += (b_aux->s[uid]*b_aux->s[t]*o[i].w); + } + return w - (w_off*4);//2 for self; 4 for both directions +} + +void cal_clus_sc0(mc_match_t *ma, mc_svaux_t *b_aux, mc_clus_t* bc, uint8_t *lock, uint8_t lock_max, +uint32_t *a, uint32_t a_n, uint32_t id, clus_flip_aux *r, uint32_t tid) +{ + uint32_t i, max_occ = (uint32_t)-1;//, len, mm0 = (uint32_t)-1, mm1 = 0; + bits_p *vis = &(r->vis); t_w_t w = 0, max_w = -1; + + memset(vis->a, 0, sizeof(*(vis->a))*vis->n); + // if(bc->dbg) { + // if(a[id] == 487 || a[id] == 47) { + // fprintf(stderr, "[M::%s::] id::%u, a[id]::%u, s::%d\n", __func__, id, a[id], b_aux->s[a[id]]); + // } + // } + for (i = id; i < a_n && a[i] != (uint32_t)-1; i++) { + if(lock[a[i]] >= lock_max) continue; + if(is_set_bits_p((*vis), a[i])) continue; + w += clus_weight(b_aux, ma, vis, a[i]); + // if(bc->dbg) { + // if(a[id] == 487 || a[id] == 47) { + // fprintf(stderr, "[M::%s::id->%u] a[%u]::%u, s::%d, lock::%u, lock_max::%u, is_set::%u, w::%f\n", + // __func__, id, i, a[i], b_aux->s[a[i]], lock[a[i]], lock_max, is_set_bits_p((*vis), a[i]), w); + // } + // } + set_bits_p((*vis), a[i]); + // mm1 = a[i]; if(a[i] < mm0) mm0 = a[i]; + ///update max_w + if(max_w < w) { + max_w = w; max_occ = i + 1 - id; + } + } + // if(mm0 != (uint32_t)-1 && mm1 != (uint32_t)-1) { + // len = MIN(((mm1>>3)+1), vis->n) - (mm0>>3); + // memset(vis->a+(mm0>>3), 0, sizeof(*(vis->a))*len); + // } + for (i = id; i < a_n && a[i] != (uint32_t)-1; i++) vis->a[a[i]>>3] = 0;///reset + if(max_w <= 0.000001 || max_occ == (uint32_t)-1) return; + if((max_w > r->w) || (max_w == r->w && id < r->off)) { + r->w = max_w; r->off = id; r->occ = max_occ; + } + // if(bc->dbg) { + // if(a[id] == 487 || a[id] == 47) { + // fprintf(stderr, "[M::%s::] id::%u, a[id]::%u, w::%f, off::%u, occ::%u\n", __func__, id, a[id], r->w, r->off, r->occ); + // } + // } +} + +static void worker_cal_clus_sc(void *data, long i, int tid) // callback for kt_for() +{ + mc_clus_t *bc = (mc_clus_t *)data; + uint32_t *a = bc->cc.nn.a, a_n = bc->cc.nn.n; + if(a[i] == (uint32_t)-1) return; + cal_clus_sc0(bc->mg->e, bc->baux, bc, bc->lock, bc->lock_max, a, a_n, i, &(bc->aux[tid]), tid); +} + + +uint32_t gen_best_clus(mc_clus_t *bc, uint32_t *off, uint32_t *occ, double *rw) +{ + uint32_t i; (*off) = (*occ) = (uint32_t)-1; (*rw) = -1; + for (i = 0; i < bc->n_thread; i++) clean_clus_flip_aux(&(bc->aux[i])); + + kt_for(bc->n_thread, worker_cal_clus_sc, bc, bc->cc.nn.n); + + for (i = 0; i < bc->n_thread; i++) { + if(bc->aux[i].off == (uint32_t)-1) continue; + if(bc->aux[i].w < 0) continue; + + if((bc->aux[i].w > (*rw)) || (bc->aux[i].w == (*rw) && bc->aux[i].off < (*off))) { + (*off) = bc->aux[i].off; (*occ) = bc->aux[i].occ; (*rw) = bc->aux[i].w; + } + } + if((*off) != (uint32_t)-1) return 1; + return 0; +} + +double flip_chain(const mc_match_t *ma, mc_svaux_t *b, uint32_t *a, uint32_t a_n, uint32_t off, uint32_t occ, +bits_p *vis, uint8_t *lock, uint8_t lock_max) +{ + uint32_t k, kn = off + occ; + for (k = off; k < kn; k++) { + vis->a[a[k]>>3] = 0; + // if(occ == 4 && a[off] == 11279) dbg = 1; + } + for (k = off; k < kn; k++) { + if(lock[a[k]] >= lock_max) continue; + if(is_set_bits_p((*vis), a[k])) continue; + set_bits_p((*vis), a[k]); lock[a[k]]++; + // if(dbg) fprintf(stderr, "[M::%s::] a[%u]::%u, s::%d\n", __func__, k, a[k], b->s[a[k]]); + mc_set_spin(ma, b, a[k], -b->s[a[k]]); + } + return mc_score(ma, b); +} + + +void test_flip_sc(const mc_match_t *ma, mc_svaux_t *b, uint32_t *a, uint32_t a_n, bits_p *vis) +{ + mc_edge_t *o = NULL; t_w_t w = 0, w_off = 0; + t_w_t sc_new = mc_score(ma, b), sc; + uint32_t n, i, t, k, uid, z; + for (k = 0; k < a_n; k++) mc_set_spin(ma, b, a[k], -b->s[a[k]]); + for (k = 0; k < a_n; k++) { + uid = a[k]; + fprintf(stderr, "+[M::%s::] uid::%u, z[0]::%f, z[1]::%f\n", __func__, uid, b->z[uid].z[0], b->z[uid].z[1]); + w += ((t_w_t)(b->s[uid])) * (b->z[uid].z[0] - b->z[uid].z[1]) * 2; + o = pt_a(*ma, uid); + n = pt_n(*ma, uid); + for (i = 0; i < n; ++i) { + t = ma_y(o[i]); + for (z = 0; z < k; z++) { + if(a[z] == t) break; + } + if(z >= k) continue; + // if(!(is_set_bits_p((*vis), t))) continue; + // if(vis[t] == 0) continue; + if(t == uid) continue; + fprintf(stderr, "+[M::%s::] uid::%u, t::%u, w::%f\n", __func__, uid, t, o[i].w); + w_off += (b->s[uid]*b->s[t]*o[i].w); + } + } + w -= (w_off*4); + sc = mc_score(ma, b); + fprintf(stderr, "+[M::%s::] sc::%f, sc_new::%f, w::%f\n", __func__, sc, sc_new, w); + + w = w_off = 0; memset(vis->a, 0, sizeof(*(vis->a))*vis->n); + for (k = 0; k < a_n; k++) { + uid = a[k]; + fprintf(stderr, "-[M::%s::] uid::%u, z[0]::%f, z[1]::%f\n", __func__, uid, b->z[uid].z[0], b->z[uid].z[1]); + w += ((t_w_t)(b->s[uid])) * (b->z[uid].z[0] - b->z[uid].z[1]) * 2; + o = pt_a(*ma, uid); + n = pt_n(*ma, uid); + for (i = 0; i < n; ++i) { + t = ma_y(o[i]); + // for (z = 0; z < k; z++) { + // if(a[z] == t) break; + // } + // if(z >= k) continue; + if(!(is_set_bits_p((*vis), t))) continue; + // if(vis[t] == 0) continue; + if(t == uid) continue; + fprintf(stderr, "-[M::%s::] uid::%u, t::%u, w::%f\n", __func__, uid, t, o[i].w); + w_off += (b->s[uid]*b->s[t]*o[i].w); + } + set_bits_p((*vis), uid); + } + w -= (w_off*4); + sc = mc_score(ma, b); + fprintf(stderr, "-[M::%s::] sc::%f, sc_new::%f, w::%f\n", __func__, sc, sc_new, w); +} + +t_w_t mc_solve_clus(mc_clus_t *bc) +{ + uint32_t off, occ/**, r = 0**/; double rw, sc, sc_opt = mc_score(bc->mg->e, bc->baux); + memset(bc->lock, 0, bc->n*sizeof(*(bc->lock))); + while (gen_best_clus(bc, &off, &occ, &rw)) { + sc = flip_chain(bc->mg->e, bc->baux, bc->cc.nn.a, bc->cc.nn.n, off, occ, &(bc->aux[0].vis), bc->lock, bc->lock_max); + // if(r%10000) { + // fprintf(stderr, "[M::%s::] rw::%f, sc_opt::%f, sc::%f, off::%u, occ::%u, r::%u\n", + // __func__, rw, sc_opt, sc, off, occ, r); + // } + if(sc < sc_opt) { + fprintf(stderr, "\nwrong::[M::%s::] rw::%f, sc_opt::%f, sc::%f, off::%u, occ::%u\n", + __func__, rw, sc_opt, sc, off, occ); + // if(occ == 2) { + // uint32_t k; + // for (k = off; k < off + occ; k++) { + // fprintf(stderr, "[M::%s::] a[%u]::%u, s::%d\n", __func__, k, bc->cc.nn.a[k], bc->baux->s[bc->cc.nn.a[k]]); + // } + // test_flip_sc(bc->mg->e, bc->baux, bc->cc.nn.a+off, occ, &(bc->aux[0].vis)); + // } + } + sc_opt = sc; //r++; + } + return sc_opt; +} + +t_w_t mc_clus_cc(mc_clus_t *bc) +{ + // fprintf(stderr, "+[M::%s::]\n", __func__); + // double index_time = yak_realtime(); + // uint32_t r = 1; + t_w_t sc_opt, sc; + // fprintf(stderr, "-[M::%s::]\n", __func__); + mc_reset_z(bc->mg->e, bc->baux); + // fprintf(stderr, "*[M::%s::]\n", __func__); + sc_opt = mc_score(bc->mg->e, bc->baux); + // fprintf(stderr, "[M::%s::] sc_opt: %f\n", __func__, sc_opt); + while (1) { + sc = mc_solve_clus(bc); + // fprintf(stderr, "[M::%s::# round: %u] sc_opt: %f, sc: %f\n", __func__, r, sc_opt, sc); + if(sc <= (sc_opt+0.0000001)) break; + sc_opt = sc; //r++; + } + // fprintf(stderr, "[M::%s::%.3f] ==> round %u\n", __func__, yak_realtime()-index_time, r); + return sc; +} + +uint32_t mc_solve_cc_adv(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint32_t cc_off, uint32_t cc_size, mc_clus_t *bc) +{ + uint32_t j, k, n_iter = 0, flush = opt->max_iter * 50, n_skip, n_skip_flush = opt->n_perturb/16; + t_w_t sc_opt = -(1<<30), sc;///problem-w + b->cc_off = cc_off, b->cc_size = cc_size; + if (b->cc_size < 2) return 0; + sc_opt = mc_init_spin(mg->e, b); + // print_sc(opt, mg, b, sc_opt, n_iter); + if (b->cc_size == 2) return 0; + for (j = 0; j < b->cc_size; ++j) {///backup s and z in s_opt and z_opt + b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]]; ///hap status of each unitig + b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; ///z[0]: positive weight; z[1]: positive weight + } + renew_mc_clus_t(bc, b->cc_node, b->cc_size); + // fprintf(stderr, "\ncc_size: %u, cc_off: %u\n", b->cc_size, b->cc_off); + // print_sc(opt, mg->e, b, sc_opt, n_iter); + sc = mc_optimize_local(opt, mg->e, b, &n_iter); + if (sc > sc_opt) { + for (j = 0; j < b->cc_size; ++j) { + b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]]; + b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; + } + sc_opt = sc; + } else { + for (j = 0; j < b->cc_size; ++j) { + b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]]; + b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]]; + } + } + // print_mc_node(mg->e, b, 36880); + // print_sc(opt, mg, b, sc_opt, n_iter); + // mc_reset_z_debug(mg->e, b); + // print_sc(opt, mg->e, b, sc_opt, n_iter); + // fprintf(stderr, "\ncc_size: %u, cc_off: %u\n", b->cc_size, b->cc_off); + + if(bc) { + sc = mc_clus_cc(bc); + if (sc > sc_opt) { + for (j = 0; j < b->cc_size; ++j) { + b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]]; + b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; + } + sc_opt = sc; + } else { + for (j = 0; j < b->cc_size; ++j) { + b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]]; + b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]]; + } + } + } + + for (k = n_skip = 0; k < (uint32_t)opt->n_perturb; ++k) { + if (k&1) mc_perturb(opt, mg->e, b); + else mc_perturb_node(opt, mg->e, b, 3); + sc = mc_optimize_local(opt, mg->e, b, &n_iter); + // if((k%256) == 0) fprintf(stderr, "(%u) sc_opt::%f, sc::%f, flush::%u, n_iter::%u\n", k, sc_opt, sc, flush, n_iter); + if (sc > sc_opt) { + for (j = 0; j < b->cc_size; ++j) { + b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]]; + b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; + } + sc_opt = sc; n_skip = 0; + // print_sc(opt, mg, b, sc_opt, n_iter); + } else { + for (j = 0; j < b->cc_size; ++j) { + b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]]; + b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]]; + } + n_skip++; + } + if(n_skip >= n_skip_flush && bc) { + sc = mc_clus_cc(bc); + if (sc > sc_opt) { + for (j = 0; j < b->cc_size; ++j) { + b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]]; + b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; + } + sc_opt = sc; + } else { + for (j = 0; j < b->cc_size; ++j) { + b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]]; + b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]]; + } + } + n_skip = 0; + } + + if((n_iter%flush) == 0) { + mc_reset_z(mg->e, b); + sc = mc_score(mg->e, b); + + for (j = 0; j < b->cc_size; ++j) { + b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; + } + sc_opt = sc; + } + } + if(bc) { + sc = mc_clus_cc(bc); + if (sc > sc_opt) { + for (j = 0; j < b->cc_size; ++j) { + b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]]; + b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]]; + } + sc_opt = sc; + } else { + for (j = 0; j < b->cc_size; ++j) { + b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]]; + b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]]; + } + } + } + + // if(bc) { + // bc->dbg = 1; + // mc_clus_cc(bc); + // bc->dbg = 0; + // } + + for (j = 0; j < b->cc_size; ++j) + { + b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]]; + b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]]; + } + // exit(1); + return n_iter; +} + void reset_mb_g_t_z(mb_g_t *mbg); uint32_t mb_solve_cc(const mc_opt_t *opt, mb_g_t *mbg, mb_svaux_t *b, uint32_t cc_off, uint32_t cc_size) @@ -2689,6 +3225,63 @@ void mc_solve_core(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub) fprintf(stderr, "[M::%s::%.3f] ==> Partition\n", __func__, yak_realtime()-index_time); } +mc_clus_t *init_mc_clus_t(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub, uint32_t n_thread, mc_svaux_t *b, uint32_t flp_max) +{ + if((!bub)) return NULL; + mc_clus_t *p; CALLOC(p, 1); + p->bub = bub; p->opt = opt; p->mg = mg; p->n = bub->ug->g->n_seq; + CALLOC(p->lock, p->n); + if(n_thread > 64) n_thread = 64; p->n_thread = n_thread; + CALLOC(p->aux, p->n_thread); + uint32_t k, ss = (p->n>>3)+(!!(p->n&7)); + for (k = 0; k < p->n_thread; k++) { + kv_resize(uint8_t, p->aux[k].vis, ss); p->aux[k].vis.n = ss; + } + p->baux = b; p->lock_max = ((flp_max<=255)?flp_max:255); if(p->lock_max < 1) p->lock_max = 1; + return p; +} + +void mc_solve_core_adv(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub) +{ + double index_time = yak_realtime(); + uint32_t st, i; + mc_svaux_t *b; mc_clus_t *bc; + // mc_bp_t *bp = NULL; + mc_g_cc(mg->e); + b = mc_svaux_init(mg, opt->seed); + bc = init_mc_clus_t(opt, mg, bub, asm_opt.thread_num, b, 16); + // bc = gen_mc_clus_t(mg->e, b, bub, ref, asm_opt.thread_num); + // if(bub) bp = mc_bp_t_init(mg->e, b, bub, asm_opt.thread_num); + /*******************************for debug************************************/ + // if(bp) mc_init_spin_all(opt, mg, NULL, b); + // if(bp) mc_solve_bp(bp); + /*******************************for debug************************************/ + if(VERBOSE_CUT) + { + fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a)); + } + + for (st = 0, i = 1; i <= mg->e->n_seq; ++i) { + if (i == mg->e->n_seq || mg->e->cc[st]>>32 != mg->e->cc[i]>>32) { + mc_solve_cc_adv(opt, mg, b, st, i - st, bc); + st = i; + } + } + + if(VERBOSE_CUT) + { + fprintf(stderr, "##############end-[---M::%s::score->%f] ==> Partition\n", __func__, mc_score_all(mg->e, b)); + mc_status_all(mg->e, mg->s.a); + } + + + // if(bp) mc_solve_bp(bp); + ///mc_write_info(g, b); + mc_svaux_destroy(b); + // if(bp) destroy_mc_bp_t(&bp); + fprintf(stderr, "[M::%s::%.3f] ==> Partition\n", __func__, yak_realtime()-index_time); +} + void set_p_flag(mc_g_t *mg, uint32_t uID, uint8_t* trio_flag, trans_chain* t_ch, int8_t s) { uint32_t i; @@ -2908,6 +3501,7 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u { if(is_dump) { dump_debug_phasing(MC_NAME, ta, ug, read_g, f_rate, renew_s, s, is_sys, bub, ref); + // bub = NULL; } mc_opt_t opt; @@ -2920,7 +3514,8 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u mb_solve_core(&opt, mg, ref, is_sys); ///debug_mc_g_t(mg); // if(renew_s == 0) write_mc_g_t(&opt, mg, MC_NAME); - mc_solve_core(&opt, mg, bub); + // mc_solve_core(&opt, mg, bub); + mc_solve_core_adv(&opt, mg, bub); if((asm_opt.flag & HA_F_PARTITION) && t_ch) { @@ -4005,9 +4600,10 @@ void quick_debug_phasing(const char* fn) mc_solve(NULL, NULL, ta, ug, read_g, f_rate, NULL, renew_s, s, is_sys, bub, ref, 0, 0); + // mc_solve_core_adv(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub, kv_u_trans_t *ref) for (k = 0; k < ug->g->n_seq; k++) { - fprintf(stderr, "utg%.6dl(len::%u)\n", (int32_t)(k)+1, ug->g->seq[k].len); + fprintf(stderr, "utg%.6dl(len::%u), s[k]::%d\n", (int32_t)(k)+1, ug->g->seq[k].len, s[k]); }