From 2fc226826817cee600fcf3728c55ab8eb01d9831 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Tue, 27 Apr 2021 21:35:49 -0400 Subject: [PATCH 1/5] phasing improvement --- CommandLines.h | 2 +- Overlaps.cpp | 12 +- hic.cpp | 424 +++++++++++++++++++++++++++++++++++++++++++++---- rcut.cpp | 2 + 4 files changed, 403 insertions(+), 37 deletions(-) diff --git a/CommandLines.h b/CommandLines.h index 8e5e64d..b71d83d 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.15.1-r329" +#define HA_VERSION "0.15.1-r330" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index 37f5820..915c447 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -12764,12 +12764,12 @@ long long gap_fuzz, bub_label_t* b_mask_t) hic_analysis(ug, sg, cov?cov->t_ch:t_ch); - // char* gfa_name = (char*)malloc(strlen(output_file_name)+25); - // sprintf(gfa_name, "%s.d_utg.noseq.gfa", output_file_name); - // FILE* output_file = fopen(gfa_name, "w"); - // ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); - // fclose(output_file); - // free(gfa_name); + char* gfa_name = (char*)malloc(strlen(output_file_name)+25); + sprintf(gfa_name, "%s.d_utg.noseq.gfa", output_file_name); + FILE* output_file = fopen(gfa_name, "w"); + ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); + fclose(output_file); + free(gfa_name); if(cov) destory_hap_cov_t(&cov); diff --git a/hic.cpp b/hic.cpp index ef47805..75b0b99 100644 --- a/hic.cpp +++ b/hic.cpp @@ -4094,8 +4094,218 @@ void update_ug_by_tigs(asg_t *sg, hc_links *link) free(flag); kv_destroy(buf.a); } +/** +#define Get_Ovlp(xs, xe, ys, ye) ((a).occ) +///take care of circle +void idx_hc_links(kvec_pe_hit* hits, ha_ug_index* idx, bubble_type* bub); +void measure_disconnected_dis(ha_ug_index* idx, kvec_pe_hit *hits, hc_links *link, bubble_type *bub) +{ + uint32_t k, l, i, m, h_occ; + uint32_t qs[2], qe[2], ts[2], te[2], len; + uint32_t q_beg, q_end, t_beg, t_end, ovlpS[2], ovlpE[2]; + uint64_t shif = 64 - idx->uID_bits, qn, tn, u_dis; + pe_hit *h_a = NULL; + hc_linkeage *t = NULL; + u_trans_t *p = NULL; + kv_u_trans_t k_trans; + kv_init(k_trans); + if(hits->idx.n == 0) idx_hc_links(hits, idx, bub); -void measure_distance(const ma_ug_t* ug, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, kv_u_trans_t *ta) + for (qn = 0; qn < hits->idx.n; qn++) + { + if(IF_HOM(qn, *bub)) continue; + h_a = hits->a.a + (hits->idx.a[qn]>>32); + h_occ = (uint32_t)(hits->idx.a[qn]); + len = idx->ug->g->seq[qn].len; + qs[0] = 0; + qe[0] = (len&1? ((len>>1) + 1) : (len>>1)); + qs[1] = (len>>1); + qe[1] = len; + + for (k = 1, l = 0; k <= h_occ; ++k) ///same qn + { + if (k == h_occ || ((h_a[k].e<<1)>>shif) != ((h_a[l].e<<1)>>shif)) //same qn and tn + { + tn = ((h_a[l].e<<1)>>shif); + if(!IF_HOM(tn, *bub) && tn != qn) + { + t = &(link->a.a[qn]); + for (i = 0, u_dis = (uint64_t)-1; i < t->e.n; i++) + { + if(t->e.a[i].del || t->e.a[i].uID != tn) continue; + u_dis = (t->e.a[i].dis ==(uint64_t)-1? (uint64_t)-1 : t->e.a[i].dis>>3); + break; + } + + if(u_dis == (uint64_t)-1) + { + len = idx->ug->g->seq[tn].len; + ts[0] = 0; + te[0] = (len&1? ((len>>1) + 1) : (len>>1)); + ts[1] = (len>>1); + te[1] = len; + + for (i = l; i < k; i++) + { + ovlpS[0] = ovlpS[1] = 0; + ovlpE[0] = ovlpE[1] = 0; + interpr_hit(idx, h_a[i].s, h_a[i].len>>32, NULL, &s_beg, &s_end); + ovlpS[0] = ((MIN(i_tEcur, q->tEcur) > MAX(i_tScur, q->tScur))? + MIN(i_tEcur, q->tEcur) - MAX(i_tScur, q->tScur):0); + + ovlp = ((MIN(i_tEcur, q->tEcur) > MAX(i_tScur, q->tScur))? + MIN(i_tEcur, q->tEcur) - MAX(i_tScur, q->tScur):0); + + + interpr_hit(idx, h_a[i].e, (uint32_t)h_a[i].len, NULL, &e_beg, &e_end); + } + + + + + kv_pushp(u_trans_t, k_trans, &p); + p->qn = qn; p->tn = tn; p->occ = (k-l); + kv_pushp(u_trans_t, k_trans, &p); + p->qn = tn; p->tn = qn; p->occ = (k-l); + } + } + l = k; + } + } + } + + + radix_sort_u_trans_m(k_trans.a, k_trans.a + k_trans.n); + + for (k = 1, l = 0, m = 0; k <= k_trans.n; ++k) + { + if (k == k_trans.n || k_trans.a[l].qn != k_trans.a[k].qn || k_trans.a[l].tn != k_trans.a[k].tn) //same qn and tn + { + for (i = l, h_occ = 0; i < k; i++) + { + h_occ += k_trans.a[i].occ; + } + + k_trans.a[m] = k_trans.a[l]; + k_trans.a[m].occ = ((uint32_t)-1) - h_occ; + m++; + l = k; + } + } + k_trans.n = m; + + radix_sort_u_trans_occ(k_trans.a, k_trans.a + k_trans.n); + for (i = 0; i < k_trans.n; i++) + { + qn = k_trans.a[i].qn; + tn = k_trans.a[i].tn; + fprintf(stderr, "s-utg%.6lul<--->d-utg%.6lul(occ: %u)\n", qn + 1, tn + 1, + ((uint32_t)-1) - k_trans.a[i].occ); + } + + fprintf(stderr, "########hits########\n"); + char dir[2] = {'+', '-'}; + for (k = 0; k < hits->a.n; ++k) + { + fprintf(stderr, "%c\tutg%.6dl(len-%u)\t%lu\t%c\tutg%.6dl(len-%u)\t%lu\ti:%lu\n", + dir[hits->a.a[k].s>>63], (int)((hits->a.a[k].s<<1)>>shif)+1, + idx->ug->g->seq[((hits->a.a[k].s<<1)>>shif)].len, hits->a.a[k].s&idx->pos_mode, + dir[hits->a.a[k].e>>63], (int)((hits->a.a[k].e<<1)>>shif)+1, + idx->ug->g->seq[((hits->a.a[k].e<<1)>>shif)].len, hits->a.a[k].e&idx->pos_mode, + hits->a.a[k].id); + } + kv_destroy(k_trans); +} +**/ +void idx_hc_links(kvec_pe_hit* hits, ha_ug_index* idx, bubble_type* bub); +void filter_disconnect_edges(ha_ug_index* idx, kvec_pe_hit *hits, hc_links *link, bubble_type *bub, uint32_t thres, double rate) +{ + uint32_t k, l, i, m, h_occ, *occ = NULL; + uint64_t shif = 64 - idx->uID_bits, qn, tn, u_dis; + pe_hit *h_a = NULL; + hc_linkeage *t = NULL; + u_trans_t *p = NULL; + hc_edge *e = NULL; + kv_u_trans_t k_trans; + kv_init(k_trans); + + if(hits->idx.n == 0) idx_hc_links(hits, idx, bub); + CALLOC(occ, hits->idx.n); + for (qn = 0; qn < hits->idx.n; qn++) + { + if(IF_HOM(qn, *bub)) continue; + h_a = hits->a.a + (hits->idx.a[qn]>>32); + h_occ = (uint32_t)(hits->idx.a[qn]); + + for (k = 1, l = 0; k <= h_occ; ++k) ///same qn + { + if (k == h_occ || ((h_a[k].e<<1)>>shif) != ((h_a[l].e<<1)>>shif)) //same qn and tn + { + tn = ((h_a[l].e<<1)>>shif); + if(!IF_HOM(tn, *bub) && tn != qn) + { + t = &(link->a.a[qn]); + for (i = 0, u_dis = (uint64_t)-1; i < t->e.n; i++) + { + if(t->e.a[i].del || t->e.a[i].uID != tn) continue; + u_dis = (t->e.a[i].dis ==(uint64_t)-1? (uint64_t)-1 : t->e.a[i].dis>>3); + break; + } + + if(u_dis == (uint64_t)-1) + { + kv_pushp(u_trans_t, k_trans, &p); + p->qn = qn; p->tn = tn; p->occ = (k-l); + kv_pushp(u_trans_t, k_trans, &p); + p->qn = tn; p->tn = qn; p->occ = (k-l); + } + else + { + occ[qn] += (k-l); occ[tn] += (k-l); + } + } + l = k; + } + } + } + + + radix_sort_u_trans_m(k_trans.a, k_trans.a + k_trans.n); + + for (k = 1, l = 0, m = 0; k <= k_trans.n; ++k) + { + if (k == k_trans.n || k_trans.a[l].qn != k_trans.a[k].qn || k_trans.a[l].tn != k_trans.a[k].tn) //same qn and tn + { + if(k - l > 2) fprintf(stderr, "ERROR-3\n"); + for (i = l, h_occ = 0; i < k; i++) + { + h_occ += k_trans.a[i].occ; + } + + k_trans.a[m] = k_trans.a[l]; + // k_trans.a[m].occ = ((uint32_t)-1) - h_occ; + k_trans.a[m].occ = h_occ; + + if(h_occ > thres || h_occ >= (occ[k_trans.a[m].qn]*rate) || h_occ >= (occ[k_trans.a[m].tn]*rate)) + { + e = get_hc_edge(link, k_trans.a[m].qn, k_trans.a[m].tn, 0); + if(e->dis != (uint64_t)-1) fprintf(stderr, "ERROR-3-0\n"); + e->is_cc = 1; + + e = get_hc_edge(link, k_trans.a[m].tn, k_trans.a[m].qn, 0); + if(e->dis != (uint64_t)-1) fprintf(stderr, "ERROR-3-0\n"); + e->is_cc = 1; + } + m++; + l = k; + } + } + k_trans.n = m; + + kv_destroy(k_trans); free(occ); +} + +void measure_distance(ha_ug_index* idx, const ma_ug_t* ug, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, kv_u_trans_t *ta) { // double index_time = yak_realtime(); MT M; @@ -4143,9 +4353,20 @@ void measure_distance(const ma_ug_t* ug, kvec_pe_hit* hits, hc_links* link, bubb update_containment_distance(copy_sg, ta, link); // update_dis_connected_gfa(copy_sg, link, &M); asg_destroy(copy_sg); - - destory_MT(&M); + + for (i = 0; i < link->a.n; i++) + { + for (k = 0; k < link->a.a[i].e.n; k++) + { + link->a.a[i].e.a[k].is_cc = 0; + if(link->a.a[i].e.a[k].del) continue; + if(link->a.a[i].e.a[k].dis == (uint64_t)-1) continue; + link->a.a[i].e.a[k].is_cc = 1; + } + } + + // filter_disconnect_edges(idx, hits, link, bub, 1, 0.01); // fprintf(stderr, "[M::%s::%.3f] ==> Hi-C linkages have been counted\n", __func__, yak_realtime()-index_time); return; } @@ -5828,17 +6049,22 @@ G_partition* clean_bubbles(hc_links* link, bubble_type* bub, min_cut_t* m, const return x; } -uint64_t get_hic_distance(pe_hit* hit, hc_links* link, const ha_ug_index* idx) +uint64_t get_hic_distance(pe_hit* hit, hc_links* link, const ha_ug_index* idx, uint32_t *is_cc) { uint64_t s_uid, s_dir, e_uid, e_dir, u_dis, k; long long s_pos, e_pos; s_uid = ((hit->s<<1)>>(64 - idx->uID_bits)); s_pos = hit->s & idx->pos_mode; e_uid = ((hit->e<<1)>>(64 - idx->uID_bits)); e_pos = hit->e & idx->pos_mode; - if(s_uid == e_uid) return MAX(s_pos, e_pos) - MIN(s_pos, e_pos); + if(s_uid == e_uid) + { + if(is_cc) (*is_cc) = 1; + return MAX(s_pos, e_pos) - MIN(s_pos, e_pos); + } hc_linkeage* t = &(link->a.a[s_uid]); for (k = 0; k < t->e.n; k++) { if(t->e.a[k].del || t->e.a[k].uID != e_uid) continue; + if(is_cc) (*is_cc) = t->e.a[k].is_cc; s_dir = (!!(t->e.a[k].dis&(uint64_t)2)); e_dir = (!!(t->e.a[k].dis&(uint64_t)1)); u_dis = (t->e.a[k].dis ==(uint64_t)-1? (uint64_t)-1 : t->e.a[k].dis>>3); @@ -5890,27 +6116,35 @@ inline double get_trans_weight_advance(const ha_ug_index* idx, uint64_t x, trans { long double rate = 0; - if(x == (uint64_t)-1) x = dis->med; - if(x < dis->max) + ///if(x == (uint64_t)-1) x = dis->med; + if(x != (uint64_t)-1) { - uint64_t i; - for (i = 0; i < dis->n; i++) + if(x < dis->max) { - if(x < dis->a[i].end && x >= dis->a[i].beg) break; - } - if(i < dis->n) - { - rate = ((double)(dis->a[i].cnt_1))/((double)(dis->a[i].cnt_0 + dis->a[i].cnt_1)); + uint64_t i; + for (i = 0; i < dis->n; i++) + { + if(x < dis->a[i].end && x >= dis->a[i].beg) break; + } + if(i < dis->n) + { + rate = ((double)(dis->a[i].cnt_1))/((double)(dis->a[i].cnt_0 + dis->a[i].cnt_1)); + } + else + { + rate = get_trans(idx, x); + } } else { rate = get_trans(idx, x); - } + } } else { - rate = get_trans(idx, x); + rate = 0.2; } + if(rate < 0) rate = 0; rate += OFFSET_RATE; @@ -6010,7 +6244,7 @@ void weight_edges_advance(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, b if(IF_HOM(beg, *bub)) continue; if(IF_HOM(end, *bub)) continue; - t_d = get_hic_distance(&(hits->a.a[k]), link, idx); + t_d = get_hic_distance(&(hits->a.a[k]), link, idx, NULL); if(t_d == (uint64_t)-1) continue; e1 = get_hc_edge(link, beg, end, 0); @@ -9083,7 +9317,7 @@ H_partition* hap, int8_t *s, trans_idx* dis) if(IF_HOM(end, *bub)) continue; - t_d = get_hic_distance(&(hits->a.a[k]), link, idx); + t_d = get_hic_distance(&(hits->a.a[k]), link, idx, NULL); if(t_d == (uint64_t)-1) continue; if(beg == end) { @@ -9257,7 +9491,7 @@ H_partition* hap, int8_t *s, trans_idx* dis) if(IF_HOM(end, *bub)) continue; if(beg == end) continue; - t_d = get_hic_distance(&(hits->a.a[k]), link, idx); + t_d = get_hic_distance(&(hits->a.a[k]), link, idx, NULL); if(t_d == (uint64_t)-1) continue; kv_push(uint64_t, buf, t_d); } @@ -13748,7 +13982,7 @@ int alignment_worker_pipeline(sldat_t* sl, const enzyme *fn1, const enzyme *fn2) return 1; } -void debug_gfa_space(ma_ug_t* ug, trans_chain* t_ch, kv_u_trans_t *ref) +void debug_gfa_space(ha_ug_index* idx, ma_ug_t* ug, trans_chain* t_ch, kv_u_trans_t *ref) { bubble_type bub; memset(&bub, 0, sizeof(bubble_type)); @@ -13759,7 +13993,7 @@ void debug_gfa_space(ma_ug_t* ug, trans_chain* t_ch, kv_u_trans_t *ref) hc_links link; init_hc_links(&link, ug->g->n_seq, t_ch); - measure_distance(ug, NULL, &link, &bub, &(t_ch->k_trans)); + measure_distance(idx, ug, NULL, &link, &bub, &(t_ch->k_trans)); // uint32_t i, k; // for (i = 0; i < link.a.n; ++i) @@ -13858,6 +14092,7 @@ kv_u_trans_t *ta, trans_idx* dis) u_trans_t *e1 = NULL, *e2 = NULL; long double weight; u_trans_t *p = NULL; + uint32_t is_cc; for (i = 0, ta->idx.n = ta->n = 0; i < link->a.n; i++) { @@ -13886,8 +14121,10 @@ kv_u_trans_t *ta, trans_idx* dis) if(IF_HOM(beg, *bub)) continue; if(IF_HOM(end, *bub)) continue; - t_d = get_hic_distance(&(hits->a.a[k]), link, idx); - if(t_d == (uint64_t)-1) continue; + t_d = get_hic_distance(&(hits->a.a[k]), link, idx, &is_cc); + // if(t_d == (uint64_t)-1) continue; + // if(t_d == (uint64_t)-1 && is_cc == 0) continue; + // if(t_d == (uint64_t)-1 && !dis) continue; get_u_trans_spec(ta, beg, end, &e1, NULL); get_u_trans_spec(ta, end, beg, &e2, NULL); @@ -13956,7 +14193,7 @@ pe_hit *hits, uint32_t occ, uint32_t qid, uint32_t qs, uint32_t qe, uint32_t tid if(e_uid != tid) break; if(!(ts <= e_beg && te >= e_end)) continue; - t_d = get_hic_distance(&hits[k], link, idx); + t_d = get_hic_distance(&hits[k], link, idx, NULL); if(t_d == (uint64_t)-1) continue; weight = 1; @@ -13976,7 +14213,7 @@ pe_hit *hits, uint32_t occ, uint32_t qid, uint32_t qs, uint32_t qe, uint32_t tid if(e_uid != tid) break; if(!(ts <= e_beg && te >= e_end)) continue; - t_d = get_hic_distance(&hits[k], link, idx); + t_d = get_hic_distance(&hits[k], link, idx, NULL); if(t_d == (uint64_t)-1) continue; weight = 1; @@ -14151,7 +14388,7 @@ inline uint32_t get_trans_interval_weight(ha_ug_index* idx, hc_links* link, bubb trans_idx* dis, pe_hit *hit, uint32_t hit_n, uint32_t qn, uint32_t qs, uint32_t qe, uint32_t tn, uint32_t ts, uint32_t te, double *w_a) { - uint32_t i, s_uid, s_beg, s_end, e_uid, e_beg, e_end, found; + uint32_t i, s_uid, s_beg, s_end, e_uid, e_beg, e_end, found, is_cc; uint64_t t_d; double weight; (*w_a) = 0; found = 0; @@ -14165,8 +14402,10 @@ uint32_t ts, uint32_t te, double *w_a) if(e_uid != tn) continue; if(!(ts <= e_beg && te >= e_end)) continue; - t_d = get_hic_distance(&hit[i], link, idx); - if(t_d == (uint64_t)-1) continue; + t_d = get_hic_distance(&hit[i], link, idx, &is_cc); + // if(t_d == (uint64_t)-1) continue; + // if(t_d == (uint64_t)-1 && is_cc == 0) continue; + // if(t_d == (uint64_t)-1 && !dis) continue; weight = 1; if(dis) weight = get_trans_weight_advance(idx, t_d, dis); @@ -14265,7 +14504,10 @@ double merge_u_trans_list(u_trans_t* a, uint32_t a_n) weight += a[i].nw; } - w += (weight/(double)(a[l].occ)); + /*******************************for debug************************************/ + // w += (weight/(double)(a[l].occ)); + w += weight; + /*******************************for debug************************************/ l = k; } } @@ -14353,6 +14595,29 @@ void print_kv_weight(kv_u_trans_t *ta) } } } + +void print_debug_hc_links(ha_ug_index* idx, bubble_type* bub, hc_links* lk, kv_u_trans_t *ta, kvec_pe_hit* hits) +{ + uint64_t k, len = 0, shif = 64 - idx->uID_bits, beg, end; + hc_edge *e = NULL; + for (k = 0; k < hits->a.n; ++k) + { + beg = ((hits->a.a[k].s<<1)>>shif); + end = ((hits->a.a[k].e<<1)>>shif); + if(beg == end) continue; + if(IF_HOM(beg, *bub)) continue; + if(IF_HOM(end, *bub)) continue; + len += (hits->a.a[k].len>>32) + ((uint32_t)hits->a.a[k].len); + } + + fprintf(stderr, "# total Hi-C aligned bases: %lu\n", len); + for (k = 0; k < ta->n; k++) + { + e = get_hc_edge(lk, ta->a[k].qn, ta->a[k].tn, 0); + fprintf(stderr, "s-utg%.6ul\td-utg%.6ul\tD:%lu\tW:%f\n", + ta->a[k].qn+1, ta->a[k].tn+1, e->dis == (uint64_t)-1? (uint64_t)-1 : e->dis>>3, ta->a[k].nw); + } +} void renew_kv_u_trans(kv_u_trans_t *ta, hc_links *lk, kvec_pe_hit* hits, kv_u_trans_t *ref, ha_ug_index* idx, bubble_type* bub, int8_t *s, uint32_t ignore_dis) { @@ -14380,6 +14645,10 @@ ha_ug_index* idx, bubble_type* bub, int8_t *s, uint32_t ignore_dis) if(hits->idx.n == 0) idx_hc_links(hits, idx, bub); weight_kv_u_trans(idx, hits, lk, bub, ta, is_comples_weight == 1? &dis : NULL); + // if(bub->round_id == bub->n_round-1) + // { + // print_debug_hc_links(idx, bub, lk, ta, hits); + // } // adjust_weight_kv_u_trans(idx, hits, lk, bub, ta, ref, is_comples_weight == 1? &dis : NULL); adjust_weight_kv_u_trans_advance(idx, hits, lk, bub, ta, ref, is_comples_weight == 1? &dis : NULL); kv_destroy(dis); @@ -14429,6 +14698,94 @@ void verbose_het_stat(bubble_type *bub) fprintf(stderr, "[M::stat] # heterozygous bases: %lu; # homozygous bases: %lu\n", hetBase, homBase); } +void debug_output_disconnected_hits(ha_ug_index* idx, kvec_pe_hit *hits, hc_links *link, bubble_type *bub) +{ + uint32_t k, l, i, m, h_occ; + uint64_t shif = 64 - idx->uID_bits, qn, tn, u_dis; + pe_hit *h_a = NULL; + hc_linkeage *t = NULL; + u_trans_t *p = NULL; + kv_u_trans_t k_trans; + kv_init(k_trans); + + + for (qn = 0; qn < hits->idx.n; qn++) + { + if(IF_HOM(qn, *bub)) continue; + h_a = hits->a.a + (hits->idx.a[qn]>>32); + h_occ = (uint32_t)(hits->idx.a[qn]); + + for (k = 1, l = 0; k <= h_occ; ++k) ///same qn + { + if (k == h_occ || ((h_a[k].e<<1)>>shif) != ((h_a[l].e<<1)>>shif)) //same qn and tn + { + tn = ((h_a[l].e<<1)>>shif); + if(!IF_HOM(tn, *bub) && tn != qn) + { + t = &(link->a.a[qn]); + for (i = 0, u_dis = (uint64_t)-1; i < t->e.n; i++) + { + if(t->e.a[i].del || t->e.a[i].uID != tn) continue; + u_dis = (t->e.a[i].dis ==(uint64_t)-1? (uint64_t)-1 : t->e.a[i].dis>>3); + break; + } + + if(u_dis == (uint64_t)-1) + { + kv_pushp(u_trans_t, k_trans, &p); + p->qn = qn; p->tn = tn; p->occ = (k-l); + kv_pushp(u_trans_t, k_trans, &p); + p->qn = tn; p->tn = qn; p->occ = (k-l); + } + } + l = k; + } + } + } + + + radix_sort_u_trans_m(k_trans.a, k_trans.a + k_trans.n); + + for (k = 1, l = 0, m = 0; k <= k_trans.n; ++k) + { + if (k == k_trans.n || k_trans.a[l].qn != k_trans.a[k].qn || k_trans.a[l].tn != k_trans.a[k].tn) //same qn and tn + { + for (i = l, h_occ = 0; i < k; i++) + { + h_occ += k_trans.a[i].occ; + } + + k_trans.a[m] = k_trans.a[l]; + k_trans.a[m].occ = ((uint32_t)-1) - h_occ; + m++; + l = k; + } + } + k_trans.n = m; + + radix_sort_u_trans_occ(k_trans.a, k_trans.a + k_trans.n); + for (i = 0; i < k_trans.n; i++) + { + qn = k_trans.a[i].qn; + tn = k_trans.a[i].tn; + fprintf(stderr, "s-utg%.6lul<--->d-utg%.6lul(occ: %u)\n", qn + 1, tn + 1, + ((uint32_t)-1) - k_trans.a[i].occ); + } + + fprintf(stderr, "########hits########\n"); + char dir[2] = {'+', '-'}; + for (k = 0; k < hits->a.n; ++k) + { + fprintf(stderr, "%c\tutg%.6dl(len-%u)\t%lu\t%c\tutg%.6dl(len-%u)\t%lu\ti:%lu\n", + dir[hits->a.a[k].s>>63], (int)((hits->a.a[k].s<<1)>>shif)+1, + idx->ug->g->seq[((hits->a.a[k].s<<1)>>shif)].len, hits->a.a[k].s&idx->pos_mode, + dir[hits->a.a[k].e>>63], (int)((hits->a.a[k].e<<1)>>shif)+1, + idx->ug->g->seq[((hits->a.a[k].e<<1)>>shif)].len, hits->a.a[k].e&idx->pos_mode, + hits->a.a[k].id); + } + kv_destroy(k_trans); +} + int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) { double index_time = yak_realtime(); @@ -14470,7 +14827,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) identify_bubbles(idx->ug, &bub, idx->t_ch->is_r_het, &(idx->t_ch->k_trans)); if(bub.round_id == 0) { - measure_distance(idx->ug, &sl.hits, &link, &bub, &(idx->t_ch->k_trans)); + measure_distance(idx, idx->ug, &sl.hits, &link, &bub, &(idx->t_ch->k_trans)); } renew_kv_u_trans(&k_trans, &link, &sl.hits, &(idx->t_ch->k_trans), idx, &bub, s->s, 0); @@ -14481,6 +14838,13 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) (bub.round_id == 0? 1 : 0), s->s, 1, /**&bub**/NULL, &(idx->t_ch->k_trans)); /*******************************for debug************************************/ label_unitigs_sm(s->s, idx->ug); + + /*******************************for debug************************************/ + // if(bub.round_id == bub.n_round - 1) + // { + // debug_output_disconnected_hits(idx, &sl.hits, &link, &bub); + // } + /*******************************for debug************************************/ /** init_hic_advance((ha_ug_index*)sl.idx, &sl.hits, &link, &bub, &hap, 0); reset_H_partition(&hap, (bub.round_id == 0? 1 : 0)); diff --git a/rcut.cpp b/rcut.cpp index 8283e5c..e49fab8 100644 --- a/rcut.cpp +++ b/rcut.cpp @@ -2814,6 +2814,8 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u mc_opt_init(&opt, asm_opt.n_perturb, asm_opt.f_perturb, asm_opt.seed); mc_g_t *mg = init_mc_g_t(ug, read_g, s, renew_s); update_mc_edges(mg, ovlp, ta, t_ch, f_rate, is_sys); + + fprintf(stderr, "[M::%s:: # edges: %u]\n", __func__, (uint32_t)mg->e->ma.n); mb_solve_core(&opt, mg, ref, is_sys); ///debug_mc_g_t(mg); From 71e91f3fc2b6cbd6f69adf656c0d9f3446bb8eb6 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Wed, 28 Apr 2021 20:27:43 -0400 Subject: [PATCH 2/5] assgin disconnected parts --- CommandLines.h | 2 +- Overlaps.cpp | 11 ++++++----- hic.cpp | 4 +++- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/CommandLines.h b/CommandLines.h index b71d83d..04584f8 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.15.1-r330" +#define HA_VERSION "0.15.1-r331" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index 915c447..192d097 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -11298,8 +11298,9 @@ void refine_u_trans_t(u_trans_hit_t *q, kv_ca_buf_t* cb) - - if(q->tScur >= q->tEcur) fprintf(stderr, "ERROR5\n"); + ///might be equal + if(q->tScur > q->tEcur) fprintf(stderr, "ERROR5\n"); + // if(q->tScur >= q->tEcur) // { // fprintf(stderr, "\n###q->tScur: %u, s: %u, si: %u, q->tEcur: %u, e: %u, ei: %u\n", // q->tScur, s, si, q->tEcur, e, ei); @@ -11780,6 +11781,7 @@ ma_ug_t *ug, uint32_t flag, double score, const char* cmd) t_ch->c_buf.a[0].c_y_p = t_ch->c_buf.a[1].c_y_p - len; } + ///chain is [s, e) if(t_ch->c_buf.a[0].c_x_p != 0 && t_ch->c_buf.a[0].c_y_p != 0) fprintf(stderr, "ERROR1\n"); if(t_ch->c_buf.a[t_ch->c_buf.n-1].c_x_p!= pri_len && t_ch->c_buf.a[t_ch->c_buf.n-1].c_y_p!= aux_len) @@ -11817,6 +11819,8 @@ ma_ug_t *ug, uint32_t flag, double score, const char* cmd) for (i = bn; i < t_ch->k_t_b.n; i++) { kh = &(t_ch->k_t_b.a[i]); + if(kh->qEpre <= kh->qSpre) continue; + if(kh->tEpre <= kh->tSpre) continue; kv_pushp(u_trans_t, t_ch->k_trans, &kt); kt->f = flag; kt->rev = ((kh->qn ^ kh->tn) & 1); kt->del = 0; kt->qn = kh->qn>>1; kt->qs = kh->qSpre; kt->qe = kh->qEpre; @@ -11831,9 +11835,6 @@ ma_ug_t *ug, uint32_t flag, double score, const char* cmd) y_score = ((double)(kt->te-kt->ts)/(double)(t_ch->c_buf.a[t_ch->c_buf.n-1].c_y_p-t_ch->c_buf.a[0].c_y_p))*score; kt->nw = MIN(x_score, y_score); } - - // fprintf(stderr, "s-utg%.6ul\t%u\t%u\td-utg%.6ul\t%u\t%u\trev(%u)\n", - // kt->qn+1, kt->qs, kt->qe, kt->tn+1, kt->ts, kt->te, kt->rev); } diff --git a/hic.cpp b/hic.cpp index 75b0b99..1cc981f 100644 --- a/hic.cpp +++ b/hic.cpp @@ -14598,7 +14598,7 @@ void print_kv_weight(kv_u_trans_t *ta) void print_debug_hc_links(ha_ug_index* idx, bubble_type* bub, hc_links* lk, kv_u_trans_t *ta, kvec_pe_hit* hits) { - uint64_t k, len = 0, shif = 64 - idx->uID_bits, beg, end; + uint64_t k, len = 0, occ = 0, shif = 64 - idx->uID_bits, beg, end; hc_edge *e = NULL; for (k = 0; k < hits->a.n; ++k) { @@ -14608,9 +14608,11 @@ void print_debug_hc_links(ha_ug_index* idx, bubble_type* bub, hc_links* lk, kv_u if(IF_HOM(beg, *bub)) continue; if(IF_HOM(end, *bub)) continue; len += (hits->a.a[k].len>>32) + ((uint32_t)hits->a.a[k].len); + occ++; } fprintf(stderr, "# total Hi-C aligned bases: %lu\n", len); + fprintf(stderr, "# total Hi-C aligned pairs: %lu\n", occ); for (k = 0; k < ta->n; k++) { e = get_hc_edge(lk, ta->a[k].qn, ta->a[k].tn, 0); From 9b31d473790bf87568f1361f72ae317961cde04d Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Wed, 28 Apr 2021 23:10:51 -0400 Subject: [PATCH 3/5] exit(0) for --help --- main.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main.cpp b/main.cpp index ccc4ed1..0518279 100644 --- a/main.cpp +++ b/main.cpp @@ -11,7 +11,7 @@ int main(int argc, char *argv[]) int i, ret; yak_reset_realtime(); init_opt(&asm_opt); - if (!CommandLine_process(argc, argv, &asm_opt)) return 1; + if (!CommandLine_process(argc, argv, &asm_opt)) return 0; ret = ha_assemble(); destory_opt(&asm_opt); fprintf(stderr, "[M::%s] Version: %s\n", __func__, HA_VERSION); From 9c205c8271e855e681b000c89622ba91029e85c9 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Fri, 30 Apr 2021 22:44:48 -0400 Subject: [PATCH 4/5] resove tangle by hic --- CommandLines.h | 2 +- Overlaps.cpp | 13 +- hic.cpp | 917 +++++++++++++++++++++++++++++++++++++++++-------- 3 files changed, 781 insertions(+), 151 deletions(-) diff --git a/CommandLines.h b/CommandLines.h index 04584f8..52107ab 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.15.1-r331" +#define HA_VERSION "0.15.1-r333" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index 192d097..1d4d3da 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -12765,20 +12765,27 @@ long long gap_fuzz, bub_label_t* b_mask_t) hic_analysis(ug, sg, cov?cov->t_ch:t_ch); + + + char* gfa_name = (char*)malloc(strlen(output_file_name)+25); - sprintf(gfa_name, "%s.d_utg.noseq.gfa", output_file_name); + sprintf(gfa_name, "%s.clean_d_utg.noseq.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); fclose(output_file); free(gfa_name); - + + + + if(cov) destory_hap_cov_t(&cov); if(t_ch) destory_trans_chain(&t_ch); 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++) diff --git a/hic.cpp b/hic.cpp index 1cc981f..99fa1c6 100644 --- a/hic.cpp +++ b/hic.cpp @@ -41,6 +41,8 @@ KRADIX_SORT_INIT(u_trans_m, u_trans_t, u_trans_m_key, 8) #define u_trans_occ_key(a) ((a).occ) KRADIX_SORT_INIT(u_trans_occ, u_trans_t, u_trans_occ_key, member_size(u_trans_t, occ)) +#define is_hom_hit(a) ((a).id == (uint64_t)-1) + typedef struct{ kvec_t(char) name; kvec_t(uint64_t) name_Len; @@ -4094,129 +4096,6 @@ void update_ug_by_tigs(asg_t *sg, hc_links *link) free(flag); kv_destroy(buf.a); } -/** -#define Get_Ovlp(xs, xe, ys, ye) ((a).occ) -///take care of circle -void idx_hc_links(kvec_pe_hit* hits, ha_ug_index* idx, bubble_type* bub); -void measure_disconnected_dis(ha_ug_index* idx, kvec_pe_hit *hits, hc_links *link, bubble_type *bub) -{ - uint32_t k, l, i, m, h_occ; - uint32_t qs[2], qe[2], ts[2], te[2], len; - uint32_t q_beg, q_end, t_beg, t_end, ovlpS[2], ovlpE[2]; - uint64_t shif = 64 - idx->uID_bits, qn, tn, u_dis; - pe_hit *h_a = NULL; - hc_linkeage *t = NULL; - u_trans_t *p = NULL; - kv_u_trans_t k_trans; - kv_init(k_trans); - if(hits->idx.n == 0) idx_hc_links(hits, idx, bub); - - for (qn = 0; qn < hits->idx.n; qn++) - { - if(IF_HOM(qn, *bub)) continue; - h_a = hits->a.a + (hits->idx.a[qn]>>32); - h_occ = (uint32_t)(hits->idx.a[qn]); - len = idx->ug->g->seq[qn].len; - qs[0] = 0; - qe[0] = (len&1? ((len>>1) + 1) : (len>>1)); - qs[1] = (len>>1); - qe[1] = len; - - for (k = 1, l = 0; k <= h_occ; ++k) ///same qn - { - if (k == h_occ || ((h_a[k].e<<1)>>shif) != ((h_a[l].e<<1)>>shif)) //same qn and tn - { - tn = ((h_a[l].e<<1)>>shif); - if(!IF_HOM(tn, *bub) && tn != qn) - { - t = &(link->a.a[qn]); - for (i = 0, u_dis = (uint64_t)-1; i < t->e.n; i++) - { - if(t->e.a[i].del || t->e.a[i].uID != tn) continue; - u_dis = (t->e.a[i].dis ==(uint64_t)-1? (uint64_t)-1 : t->e.a[i].dis>>3); - break; - } - - if(u_dis == (uint64_t)-1) - { - len = idx->ug->g->seq[tn].len; - ts[0] = 0; - te[0] = (len&1? ((len>>1) + 1) : (len>>1)); - ts[1] = (len>>1); - te[1] = len; - - for (i = l; i < k; i++) - { - ovlpS[0] = ovlpS[1] = 0; - ovlpE[0] = ovlpE[1] = 0; - interpr_hit(idx, h_a[i].s, h_a[i].len>>32, NULL, &s_beg, &s_end); - ovlpS[0] = ((MIN(i_tEcur, q->tEcur) > MAX(i_tScur, q->tScur))? - MIN(i_tEcur, q->tEcur) - MAX(i_tScur, q->tScur):0); - - ovlp = ((MIN(i_tEcur, q->tEcur) > MAX(i_tScur, q->tScur))? - MIN(i_tEcur, q->tEcur) - MAX(i_tScur, q->tScur):0); - - - interpr_hit(idx, h_a[i].e, (uint32_t)h_a[i].len, NULL, &e_beg, &e_end); - } - - - - - kv_pushp(u_trans_t, k_trans, &p); - p->qn = qn; p->tn = tn; p->occ = (k-l); - kv_pushp(u_trans_t, k_trans, &p); - p->qn = tn; p->tn = qn; p->occ = (k-l); - } - } - l = k; - } - } - } - - - radix_sort_u_trans_m(k_trans.a, k_trans.a + k_trans.n); - - for (k = 1, l = 0, m = 0; k <= k_trans.n; ++k) - { - if (k == k_trans.n || k_trans.a[l].qn != k_trans.a[k].qn || k_trans.a[l].tn != k_trans.a[k].tn) //same qn and tn - { - for (i = l, h_occ = 0; i < k; i++) - { - h_occ += k_trans.a[i].occ; - } - - k_trans.a[m] = k_trans.a[l]; - k_trans.a[m].occ = ((uint32_t)-1) - h_occ; - m++; - l = k; - } - } - k_trans.n = m; - - radix_sort_u_trans_occ(k_trans.a, k_trans.a + k_trans.n); - for (i = 0; i < k_trans.n; i++) - { - qn = k_trans.a[i].qn; - tn = k_trans.a[i].tn; - fprintf(stderr, "s-utg%.6lul<--->d-utg%.6lul(occ: %u)\n", qn + 1, tn + 1, - ((uint32_t)-1) - k_trans.a[i].occ); - } - - fprintf(stderr, "########hits########\n"); - char dir[2] = {'+', '-'}; - for (k = 0; k < hits->a.n; ++k) - { - fprintf(stderr, "%c\tutg%.6dl(len-%u)\t%lu\t%c\tutg%.6dl(len-%u)\t%lu\ti:%lu\n", - dir[hits->a.a[k].s>>63], (int)((hits->a.a[k].s<<1)>>shif)+1, - idx->ug->g->seq[((hits->a.a[k].s<<1)>>shif)].len, hits->a.a[k].s&idx->pos_mode, - dir[hits->a.a[k].e>>63], (int)((hits->a.a[k].e<<1)>>shif)+1, - idx->ug->g->seq[((hits->a.a[k].e<<1)>>shif)].len, hits->a.a[k].e&idx->pos_mode, - hits->a.a[k].id); - } - kv_destroy(k_trans); -} -**/ void idx_hc_links(kvec_pe_hit* hits, ha_ug_index* idx, bubble_type* bub); void filter_disconnect_edges(ha_ug_index* idx, kvec_pe_hit *hits, hc_links *link, bubble_type *bub, uint32_t thres, double rate) { @@ -9315,6 +9194,7 @@ H_partition* hap, int8_t *s, trans_idx* dis) if(IF_HOM(beg, *bub)) continue; if(IF_HOM(end, *bub)) continue; + if(is_hom_hit(hits->a.a[k])) continue; t_d = get_hic_distance(&(hits->a.a[k]), link, idx, NULL); @@ -9489,6 +9369,7 @@ H_partition* hap, int8_t *s, trans_idx* dis) if(IF_HOM(beg, *bub)) continue; if(IF_HOM(end, *bub)) continue; + if(is_hom_hit(hits->a.a[k])) continue; if(beg == end) continue; t_d = get_hic_distance(&(hits->a.a[k]), link, idx, NULL); @@ -12688,11 +12569,10 @@ kvec_t_u32_warp* result, uint32_t e_flag, uint32_t ava_flag, kvec_t_u32_warp* re } root ^= 1; - ///fprintf(stderr, "root=utg%.6dl\n", (root>>1)+1); is_vis[root] = 0; stack->a.n = 0; kv_push(uint32_t, stack->a, root); - while (stack->a.n > 0) + while (stack->a.n > 0)///label all untigs not in any chain { stack->a.n--; cur = stack->a.a[stack->a.n]; @@ -12729,7 +12609,6 @@ kvec_t_u32_warp* result, uint32_t e_flag, uint32_t ava_flag, kvec_t_u32_warp* re root_source ^= 1; aim_1 = root_source>>1; - ///fprintf(stderr, "aim_0=utg%.6ul, aim_1=utg%.6ul\n", aim_0+1, aim_1+1); cur = root_source;///scan nodes that cannot be reached from root but can be reached from root_source ncur = asg_arc_n(ug->g, cur); @@ -12738,6 +12617,7 @@ kvec_t_u32_warp* result, uint32_t e_flag, uint32_t ava_flag, kvec_t_u32_warp* re { if(acur[k_i].del) continue; if(vis_flag[acur[k_i].v>>1] != 0) continue;///skip nodes that are already reachable + ///don't label any path that can reack other chains if_conflict_utg(acur[k_i].v, &aim_0, &aim_1, ug, vis_flag, is_vis, ava_flag, stack); } @@ -12810,7 +12690,7 @@ kvec_t_u32_warp* result, uint32_t e_flag, uint32_t ava_flag, kvec_t_u32_warp* re void clean_bubble_chain_by_hic(ma_ug_t* ug, kv_u_trans_t *ta, bubble_type* bub) { - // double index_time = yak_realtime(); + double index_time = yak_realtime(); ma_ug_t *bs_ug = bub->b_ug; uint32_t v, u, i, m, max_i, nv, rv, n_vx, root, flag_pri = 1, flag_aux = 2, flag_ava = 4, occ; double w, cutoff = 2; @@ -12913,7 +12793,7 @@ void clean_bubble_chain_by_hic(ma_ug_t* ug, kv_u_trans_t *ta, bubble_type* bub) av = asg_arc_a(bs_ug->g, v); nv = asg_arc_n(bs_ug->g, v); rv = get_real_length(bs_ug->g, v, NULL); - if(nv == rv) continue; + if(nv == rv) continue;///no edge drop if(rv != 1 || nv <= 1) continue; get_real_length(bs_ug->g, v, &u); u ^= 1; @@ -12943,7 +12823,6 @@ void clean_bubble_chain_by_hic(ma_ug_t* ug, kv_u_trans_t *ta, bubble_type* bub) set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, u^1, vis, flag_pri, NULL); for (; i < res_utg.a.n; i++) dedup[res_utg.a.a[i]>>1] |= 2; - for (i = m = 0; i < res_utg.a.n; i++) { if(dedup[res_utg.a.a[i]>>1] == 3) @@ -12955,9 +12834,6 @@ void clean_bubble_chain_by_hic(ma_ug_t* ug, kv_u_trans_t *ta, bubble_type* bub) } res_utg.a.n = m; - // fprintf(stderr, "res_utg.a.n: %u, m: %u, beg-utg%.6ul, sink-utg%.6ul\n", - // res_utg.a.n, m, (root_0>>1)+1, (root_1>>1)+1); - if(!IF_HOM(root_0>>1, *bub)) kv_push(uint32_t, res_utg.a, root_0); if(!IF_HOM(root_1>>1, *bub)) kv_push(uint32_t, res_utg.a, root_1); @@ -13004,6 +12880,604 @@ void clean_bubble_chain_by_hic(ma_ug_t* ug, kv_u_trans_t *ta, bubble_type* bub) free(vis); free(is_vis); free(is_used); free(dedup); free(b.b.a); free(e_w); free(e_occ); kv_destroy(stack.a); kv_destroy(result.a); kv_destroy(res_utg.a); kv_destroy(edges.a); ma_ug_destroy(back_bs_ug); + fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time); +} + +void clean_sub_tangle(bubble_type* bub, uint8_t* vis_flag, kvec_t_u32_warp *res_utg, +uint8_t *dedup, uint8_t *is_tangle, kvec_asg_arc_t_warp* edges, uint32_t flag, uint32_t s, +uint32_t e) +{ + uint32_t i, k_i, pv, v, occ, beg, sink, p_beg, p_sink, n, *a = NULL; + ma_utg_t *u = NULL; + + for (i = 0, pv = (uint32_t)-1; i < res_utg->a.n; i++) + { + if(is_tangle[res_utg->a.a[i]>>1] == 1) + { + if(pv == (uint32_t)-1) + { + pv = res_utg->a.a[i]>>1; + } + else if(pv != (res_utg->a.a[i]>>1)) + { + pv = (uint32_t)-1; + break; + } + } + else if(is_tangle[res_utg->a.a[i]>>1] == (uint8_t)-1) + { + pv = (uint32_t)-1; + break; + } + } + + if(pv == (uint32_t)-1) + { + for (i = 0; i < res_utg->a.n; i++) + { + is_tangle[res_utg->a.a[i]>>1] = (uint8_t)-1; + } + return; + } + + occ = 0; + u = &(bub->b_ug->u.a[s]); + for (i = 0, p_beg = p_sink = (uint32_t)-1; i < u->n; i++) + { + get_bubbles(bub, u->a[i]>>33, &beg, &sink, &a, &n, NULL); + + for (k_i = 0; k_i < n; k_i++) + { + occ += bub->ug->u.a[a[k_i]>>1].n; + } + if(beg != (uint32_t)-1 && (beg>>1) != (p_beg>>1) && (beg>>1) != (p_sink>>1)) + { + occ += bub->ug->u.a[beg>>1].n; + } + if(sink != (uint32_t)-1 && (sink>>1) != (p_beg>>1) && (sink>>1) != (p_sink>>1)) + { + occ += bub->ug->u.a[sink>>1].n; + } + p_beg = beg; p_sink = sink; + } + + u = &(bub->b_ug->u.a[e]); + for (i = 0, p_beg = p_sink = (uint32_t)-1; i < u->n; i++) + { + get_bubbles(bub, u->a[i]>>33, &beg, &sink, &a, &n, NULL); + + for (k_i = 0; k_i < n; k_i++) + { + occ += bub->ug->u.a[a[k_i]>>1].n; + } + if(beg != (uint32_t)-1 && (beg>>1) != (p_beg>>1) && (beg>>1) != (p_sink>>1)) + { + occ += bub->ug->u.a[beg>>1].n; + } + if(sink != (uint32_t)-1 && (sink>>1) != (p_beg>>1) && (sink>>1) != (p_sink>>1)) + { + occ += bub->ug->u.a[sink>>1].n; + } + p_beg = beg; p_sink = sink; + } + + if(occ <= bub->ug->u.a[pv].n*200) + { + for (i = 0; i < res_utg->a.n; i++) + { + is_tangle[res_utg->a.a[i]>>1] = (uint8_t)-1; + } + return; + } + + if(!edges) return; + + u = &(bub->b_ug->u.a[s]); + for (i = 0; i < u->n; i++) + { + get_bubbles(bub, u->a[i]>>33, &beg, &sink, &a, &n, NULL); + + for (k_i = 0; k_i < n; k_i++) + { + vis_flag[a[k_i]>>1] ^= flag; + } + if(beg != (uint32_t)-1) + { + vis_flag[beg>>1] ^= flag; + } + if(sink != (uint32_t)-1) + { + vis_flag[sink>>1] ^= flag; + } + } + u = &(bub->b_ug->u.a[e]); + for (i = 0; i < u->n; i++) + { + get_bubbles(bub, u->a[i]>>33, &beg, &sink, &a, &n, NULL); + + for (k_i = 0; k_i < n; k_i++) + { + vis_flag[a[k_i]>>1] ^= flag; + } + if(beg != (uint32_t)-1) + { + vis_flag[beg>>1] ^= flag; + } + if(sink != (uint32_t)-1) + { + vis_flag[sink>>1] ^= flag; + } + } + for (i = 0; i < res_utg->a.n; i++) + { + if(dedup[res_utg->a.a[i]>>1] == 3) + { + vis_flag[res_utg->a.a[i]>>1] ^= flag; + } + } + + + + asg_arc_t *av = NULL, *p = NULL; + uint32_t nv, c[2]; + + c[0] = c[1] = 0; + while (1) + { + c[0] = c[1] = 0; p = NULL; + v = pv<<1; + av = asg_arc_a(bub->ug->g, v); + nv = asg_arc_n(bub->ug->g, v); + for (i = 0; i < nv; i++) + { + if(av[i].del) continue; + c[(!!(vis_flag[av[i].v>>1]&flag))]++; + if(p == NULL || p->ol > av[i].ol) + { + p = &(av[i]); + } + } + + + + + v = (pv<<1) + 1; + av = asg_arc_a(bub->ug->g, v); + nv = asg_arc_n(bub->ug->g, v); + for (i = 0; i < nv; i++) + { + if(av[i].del) continue; + c[(!!(vis_flag[av[i].v>>1]&flag))]++; + if(p == NULL || p->ol > av[i].ol) + { + p = &(av[i]); + } + } + + if(p) + { + p->del = 1; + c[(!!(vis_flag[p->v>>1]&flag))]--; + if((p->v>>1) == pv) + { + asg_arc_del(bub->ug->g, (p->v)^1, (p->ul>>32)^1, 1); + c[0] = c[1] = 0; + v = pv<<1; + av = asg_arc_a(bub->ug->g, v); + nv = asg_arc_n(bub->ug->g, v); + for (i = 0; i < nv; i++) + { + if(av[i].del) continue; + c[(!!(vis_flag[av[i].v>>1]&flag))]++; + } + + v = (pv<<1) + 1; + av = asg_arc_a(bub->ug->g, v); + nv = asg_arc_n(bub->ug->g, v); + for (i = 0; i < nv; i++) + { + if(av[i].del) continue; + c[(!!(vis_flag[av[i].v>>1]&flag))]++; + } + } + + + } + if(c[0] == 0 || c[1] == 0) break; + } + + + + + v = pv<<1; + av = asg_arc_a(bub->ug->g, v); + nv = asg_arc_n(bub->ug->g, v); + for (i = 0; i < nv; i++) + { + if(!av[i].del) continue; + av[i].del = 0; + kv_push(asg_arc_t, edges->a, av[i]); + } + + + v = (pv<<1) + 1; + av = asg_arc_a(bub->ug->g, v); + nv = asg_arc_n(bub->ug->g, v); + for (i = 0; i < nv; i++) + { + if(!av[i].del) continue; + av[i].del = 0; + kv_push(asg_arc_t, edges->a, av[i]); + } + + is_tangle[pv] = (uint8_t)-1; + + + u = &(bub->b_ug->u.a[s]); + for (i = 0; i < u->n; i++) + { + get_bubbles(bub, u->a[i]>>33, &beg, &sink, &a, &n, NULL); + + for (k_i = 0; k_i < n; k_i++) + { + vis_flag[a[k_i]>>1] ^= flag; + } + if(beg != (uint32_t)-1) + { + vis_flag[beg>>1] ^= flag; + } + if(sink != (uint32_t)-1) + { + vis_flag[sink>>1] ^= flag; + } + } + u = &(bub->b_ug->u.a[e]); + for (i = 0; i < u->n; i++) + { + get_bubbles(bub, u->a[i]>>33, &beg, &sink, &a, &n, NULL); + + for (k_i = 0; k_i < n; k_i++) + { + vis_flag[a[k_i]>>1] ^= flag; + } + if(beg != (uint32_t)-1) + { + vis_flag[beg>>1] ^= flag; + } + if(sink != (uint32_t)-1) + { + vis_flag[sink>>1] ^= flag; + } + } + for (i = 0; i < res_utg->a.n; i++) + { + if(dedup[res_utg->a.a[i]>>1] == 3) + { + vis_flag[res_utg->a.a[i]>>1] ^= flag; + } + } + + +} + +void delete_sg_e_by_ug(asg_t* rg, ma_ug_t* ug, uint32_t v, uint32_t w) +{ + uint32_t vx, wx; + vx = (v&1?((ug->u.a[v>>1].a[0]>>32)^1):(ug->u.a[v>>1].a[ug->u.a[v>>1].n-1]>>32)); + wx = (w&1?((ug->u.a[w>>1].a[ug->u.a[w>>1].n-1]>>32)^1):(ug->u.a[w>>1].a[0]>>32)); + asg_arc_del(rg, vx, wx, 1); asg_arc_del(rg, wx^1, vx^1, 1); +} + +void resolve_bubble_chain_by_hic(ha_ug_index *idx, kv_u_trans_t *ta, bubble_type* bub) +{ + // double index_time = yak_realtime(); + ma_ug_t* ug = idx->ug; + ma_ug_t *bs_ug = bub->b_ug; + uint32_t v, u, i, max_i, nv, rv, n_vx, root, flag_pri = 1, flag_aux = 2, flag_ava = 4, occ; + double w, cutoff = 2; + uint32_t max_w_occ = 4; + asg_arc_t *av = NULL; + n_vx = bs_ug->g->n_seq << 1; + uint8_t *vis = NULL; CALLOC(vis, ug->g->n_seq<<1); + uint8_t *is_vis = NULL; CALLOC(is_vis, ug->g->n_seq<<1); + uint8_t *is_used = NULL; CALLOC(is_used, n_vx); + uint8_t *dedup = NULL; CALLOC(dedup, ug->g->n_seq<<1); + uint8_t *is_tangle = NULL; CALLOC(is_tangle, ug->g->n_seq); + buf_t b; memset(&b, 0, sizeof(buf_t)); + kvec_t_u32_warp stack, result, res_utg; + kv_init(stack.a); kv_init(result.a); kv_init(res_utg.a); + double *e_w = NULL; MALLOC(e_w, bs_ug->g->n_arc); + uint32_t *e_occ = NULL, *a_occ = NULL; CALLOC(e_occ, bs_ug->g->n_arc); + double *aw = NULL, max_w = 0; + kvec_asg_arc_t_warp edges; kv_init(edges.a); + ma_ug_t *back_bs_ug = copy_untig_graph(bs_ug); + + for (i = 0; i < bs_ug->g->n_arc; i++)///weight of bs_ug's edges + { + e_w[i] = -1; + } + + for (i = 0; i < bs_ug->g->n_seq; i++)///init all chain with flag_aux + { + set_b_utg_weight_flag(bub, &b, i<<1, vis, flag_aux, NULL); + } + + + for (v = 0; v < n_vx; v++) + { + av = asg_arc_a(bs_ug->g, v); + nv = asg_arc_n(bs_ug->g, v); + aw = (&e_w[bs_ug->g->idx[v]>>32]); + a_occ = (&e_occ[bs_ug->g->idx[v]>>32]); + if(nv <= 1 || get_real_length(bs_ug->g, v, NULL) <= 1) continue; + set_b_utg_weight_flag_xor(bub, bs_ug, &b, v^1, vis, flag_pri, NULL); + + for (i = 0; i < nv; i++) + { + if(av[i].del) continue; + w = get_chain_weight_hic(bub, bs_ug, &b, av[i].v, v, ta, vis, is_vis, ug, &stack, &result, flag_pri, flag_ava, NULL, &occ); + aw[i] = w; + a_occ[i] = occ; + } + + set_b_utg_weight_flag_xor(bub, bs_ug, &b, v^1, vis, flag_pri, NULL); + } + + + for (v = 0; v < n_vx; v++) + { + av = asg_arc_a(bs_ug->g, v); + nv = asg_arc_n(bs_ug->g, v); + aw = (&e_w[bs_ug->g->idx[v]>>32]); + a_occ = (&e_occ[bs_ug->g->idx[v]>>32]); + if(nv <= 1 || get_real_length(bs_ug->g, v, NULL) <= 1) continue; + + for (i = rv = 0, max_i = (uint32_t)-1; i < nv; i++) + { + if(av[i].del) continue; + if(max_i == (uint32_t)-1) + { + max_i = i; + max_w = aw[i]; + } + else if(max_w < aw[i]) + { + max_i = i; + max_w = aw[i]; + } + rv++; + } + + if(max_i == (uint32_t)-1) continue; + ///if(max_w <= max_w_cutoff) continue; //must be <= + if(a_occ[max_i] <= max_w_occ) continue; //must be <= + if(rv < 2) continue; + + for (i = 0; i < nv; i++) + { + if(av[i].del) continue; + if(i == max_i) continue; + ///if((av[i].v>>1) == (v>>1) && aw[i] <= max_w_cutoff) continue; ///might be not reasonable + if((av[i].v>>1) == (v>>1) && a_occ[i] <= max_w_occ) continue; ///might be not reasonable + if(aw[i]*cutoff < max_w && double_check_bub_branch(&av[i], bs_ug, e_w, e_occ, cutoff, max_w_occ)) + { + av[i].del = 1; asg_arc_del(bs_ug->g, (av[i].v)^1, (av[i].ul>>32)^1, 1); + } + } + } + + uint32_t rId_0, ori_0, rId_1, ori_1, root_0, root_1; + for (v = 0; v < n_vx; v++) + { + av = asg_arc_a(bs_ug->g, v); + nv = asg_arc_n(bs_ug->g, v); + rv = get_real_length(bs_ug->g, v, NULL); + if(nv == rv) continue;///no edge drop + if(rv != 1 || nv <= 1) continue; + get_real_length(bs_ug->g, v, &u); + u ^= 1; + if(get_real_length(bs_ug->g, u, NULL) != 1) continue; + drop_g_edges_by_utg(bub, bub->b_g, bs_ug, NULL, v, u); + if(is_used[v] || is_used[u]) continue; + + is_used[v] = is_used[u] = 1; + root = get_utg_end_from_btg(bub, bs_ug, v); + rId_0 = root>>1; + ori_0 = root&1; + get_bubbles(bub, rId_0, ori_0 == 1?&root_0:NULL, ori_0 == 0?&root_0:NULL, NULL, NULL, NULL); + + root = get_utg_end_from_btg(bub, bs_ug, u); + rId_1 = root>>1; + ori_1 = root&1; + get_bubbles(bub, rId_1, ori_1 == 1?&root_1:NULL, ori_1 == 0?&root_1:NULL, NULL, NULL, NULL); + + res_utg.a.n = 0; + set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, v^1, vis, flag_pri, NULL); + get_chain_weight_hic(bub, back_bs_ug, &b, u^1, v, ta, vis, is_vis, ug, &stack, &result, flag_pri, flag_ava, &res_utg, NULL); + set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, v^1, vis, flag_pri, NULL); + for (i = 0; i < res_utg.a.n; i++) dedup[res_utg.a.a[i]>>1] |= 1; + + set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, u^1, vis, flag_pri, NULL); + get_chain_weight_hic(bub, back_bs_ug, &b, v^1, u, ta, vis, is_vis, ug, &stack, &result, flag_pri, flag_ava, &res_utg, NULL); + set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, u^1, vis, flag_pri, NULL); + for (; i < res_utg.a.n; i++) dedup[res_utg.a.a[i]>>1] |= 2; + + for (i = 0; i < res_utg.a.n; i++) + { + if(dedup[res_utg.a.a[i]>>1] == 3) + { + is_tangle[res_utg.a.a[i]>>1] = 3; + } + else + { + if(is_tangle[res_utg.a.a[i]>>1] != 3) + { + is_tangle[res_utg.a.a[i]>>1] = 1; + } + } + dedup[res_utg.a.a[i]>>1] = 0; + } + } + + /*******************************for debug************************************/ + // for (i = 0; i < ug->g->n_seq; i++) + // { + // if(is_tangle[i] == 1) + // { + // fprintf(stderr, "*****tangle-utg%.6ul\n", i+1); + // } + // } + /*******************************for debug************************************/ + memset(is_used, 0, n_vx); + for (v = 0; v < n_vx; v++) + { + av = asg_arc_a(bs_ug->g, v); + nv = asg_arc_n(bs_ug->g, v); + rv = get_real_length(bs_ug->g, v, NULL); + if(nv == rv) continue;///no edge drop + if(rv != 1 || nv <= 1) continue; + get_real_length(bs_ug->g, v, &u); + u ^= 1; + if(get_real_length(bs_ug->g, u, NULL) != 1) continue; + drop_g_edges_by_utg(bub, bub->b_g, bs_ug, NULL, v, u); + if(is_used[v] || is_used[u]) continue; + + is_used[v] = is_used[u] = 1; + root = get_utg_end_from_btg(bub, bs_ug, v); + rId_0 = root>>1; + ori_0 = root&1; + get_bubbles(bub, rId_0, ori_0 == 1?&root_0:NULL, ori_0 == 0?&root_0:NULL, NULL, NULL, NULL); + + root = get_utg_end_from_btg(bub, bs_ug, u); + rId_1 = root>>1; + ori_1 = root&1; + get_bubbles(bub, rId_1, ori_1 == 1?&root_1:NULL, ori_1 == 0?&root_1:NULL, NULL, NULL, NULL); + + res_utg.a.n = 0; + set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, v^1, vis, flag_pri, NULL); + get_chain_weight_hic(bub, back_bs_ug, &b, u^1, v, ta, vis, is_vis, ug, &stack, &result, flag_pri, flag_ava, &res_utg, NULL); + set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, v^1, vis, flag_pri, NULL); + for (i = 0; i < res_utg.a.n; i++) dedup[res_utg.a.a[i]>>1] |= 1; + + set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, u^1, vis, flag_pri, NULL); + get_chain_weight_hic(bub, back_bs_ug, &b, v^1, u, ta, vis, is_vis, ug, &stack, &result, flag_pri, flag_ava, &res_utg, NULL); + set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, u^1, vis, flag_pri, NULL); + for (; i < res_utg.a.n; i++) dedup[res_utg.a.a[i]>>1] |= 2; + + + /*******************************for debug************************************/ + // fprintf(stderr, "\nres_utg.a.n: %u, beg-utg%.6ul, sink-utg%.6ul\n", + // (uint32_t)res_utg.a.n, (root_0>>1)+1, (root_1>>1)+1); + // for (i = 0; i < res_utg.a.n; i++) + // { + // if(dedup[res_utg.a.a[i]>>1] == 3) + // { + // fprintf(stderr, "share-utg%.6ul\n", (res_utg.a.a[i]>>1)+1); + // } + // } + // for (i = 0; i < res_utg.a.n; i++) + // { + // if(dedup[res_utg.a.a[i]>>1] == 1) + // { + // fprintf(stderr, "1-utg%.6ul\n", (res_utg.a.a[i]>>1)+1); + // } + // } + // for (i = 0; i < res_utg.a.n; i++) + // { + // if(dedup[res_utg.a.a[i]>>1] == 2) + // { + // fprintf(stderr, "2-utg%.6ul\n", (res_utg.a.a[i]>>1)+1); + // } + // } + /*******************************for debug************************************/ + + clean_sub_tangle(bub, vis, &res_utg, dedup, is_tangle, NULL, flag_pri, v>>1, u>>1); + + for (i = 0; i < res_utg.a.n; i++) + { + dedup[res_utg.a.a[i]>>1] = 0; + } + } + + + edges.a.n = 0; + memset(is_used, 0, n_vx); + for (v = 0; v < n_vx; v++) + { + av = asg_arc_a(bs_ug->g, v); + nv = asg_arc_n(bs_ug->g, v); + rv = get_real_length(bs_ug->g, v, NULL); + if(nv == rv) continue;///no edge drop + if(rv != 1 || nv <= 1) continue; + get_real_length(bs_ug->g, v, &u); + u ^= 1; + if(get_real_length(bs_ug->g, u, NULL) != 1) continue; + drop_g_edges_by_utg(bub, bub->b_g, bs_ug, NULL, v, u); + if(is_used[v] || is_used[u]) continue; + + is_used[v] = is_used[u] = 1; + root = get_utg_end_from_btg(bub, bs_ug, v); + rId_0 = root>>1; + ori_0 = root&1; + get_bubbles(bub, rId_0, ori_0 == 1?&root_0:NULL, ori_0 == 0?&root_0:NULL, NULL, NULL, NULL); + + root = get_utg_end_from_btg(bub, bs_ug, u); + rId_1 = root>>1; + ori_1 = root&1; + get_bubbles(bub, rId_1, ori_1 == 1?&root_1:NULL, ori_1 == 0?&root_1:NULL, NULL, NULL, NULL); + + res_utg.a.n = 0; + set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, v^1, vis, flag_pri, NULL); + get_chain_weight_hic(bub, back_bs_ug, &b, u^1, v, ta, vis, is_vis, ug, &stack, &result, flag_pri, flag_ava, &res_utg, NULL); + set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, v^1, vis, flag_pri, NULL); + for (i = 0; i < res_utg.a.n; i++) dedup[res_utg.a.a[i]>>1] |= 1; + + set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, u^1, vis, flag_pri, NULL); + get_chain_weight_hic(bub, back_bs_ug, &b, v^1, u, ta, vis, is_vis, ug, &stack, &result, flag_pri, flag_ava, &res_utg, NULL); + set_b_utg_weight_flag_xor(bub, back_bs_ug, &b, u^1, vis, flag_pri, NULL); + for (; i < res_utg.a.n; i++) dedup[res_utg.a.a[i]>>1] |= 2; + + clean_sub_tangle(bub, vis, &res_utg, dedup, is_tangle, &edges, flag_pri, v>>1, u>>1); + + for (i = 0; i < res_utg.a.n; i++) + { + dedup[res_utg.a.a[i]>>1] = 0; + } + } + + + /*******************************for debug************************************/ + // for (i = 0; i < ug->g->n_seq; i++) + // { + // if(is_tangle[i] == 1) + // { + // fprintf(stderr, "####tangle-utg%.6ul\n", i+1); + // } + // } + + for (i = 0; i < edges.a.n; i++) + { + fprintf(stderr, "s-utg%.6lul<------>d-utg%.6ul\n", (edges.a.a[i].ul>>33) + 1, (edges.a.a[i].v>>1) + 1); + + } + /*******************************for debug************************************/ + if(edges.a.n > 0) + { + for (i = 0; i < edges.a.n; i++) + { + asg_arc_del(bub->ug->g, edges.a.a[i].ul>>32, edges.a.a[i].v, 1); + asg_arc_del(bub->ug->g, (edges.a.a[i].v)^1, (edges.a.a[i].ul>>32)^1, 1); + delete_sg_e_by_ug(idx->read_g, idx->ug, edges.a.a[i].ul>>32, edges.a.a[i].v); + delete_sg_e_by_ug(idx->read_g, idx->ug, (edges.a.a[i].v)^1, (edges.a.a[i].ul>>32)^1); + } + asg_cleanup(bub->ug->g); + asg_cleanup(idx->read_g); + } + + + free(vis); free(is_vis); free(is_used); free(dedup); free(b.b.a); free(e_w); free(e_occ); free(is_tangle); + kv_destroy(stack.a); kv_destroy(result.a); kv_destroy(res_utg.a); kv_destroy(edges.a); + ma_ug_destroy(back_bs_ug); // fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time); } @@ -14041,8 +14515,8 @@ void idx_hc_links(kvec_pe_hit* hits, ha_ug_index* idx, bubble_type* bub) { qn = ((hits->a.a[l].s<<1)>>(64 - idx->uID_bits)); tn = ((hits->a.a[l].e<<1)>>(64 - idx->uID_bits)); - if(IF_HOM(qn, *bub)) continue; - if(IF_HOM(tn, *bub)) continue; + if(bub && IF_HOM(qn, *bub)) continue; + if(bub && IF_HOM(tn, *bub)) continue; hits->occ.a[qn]++; hits->occ.a[tn]++; } @@ -14120,6 +14594,7 @@ kv_u_trans_t *ta, trans_idx* dis) if(beg == end) continue; if(IF_HOM(beg, *bub)) continue; if(IF_HOM(end, *bub)) continue; + if(is_hom_hit(hits->a.a[k])) continue; t_d = get_hic_distance(&(hits->a.a[k]), link, idx, &is_cc); // if(t_d == (uint64_t)-1) continue; @@ -14394,6 +14869,7 @@ uint32_t ts, uint32_t te, double *w_a) (*w_a) = 0; found = 0; for (i = 0; i < hit_n; i++)///all hits already have the same qn and tn { + if(is_hom_hit(hit[i])) continue; interpr_hit(idx, hit[i].s, hit[i].len>>32, &s_uid, &s_beg, &s_end); if(s_uid != qn) continue; if(!(qs <= s_beg && qe >= s_end)) continue; @@ -14688,6 +15164,80 @@ void destory_ps_t(ps_t **s) free((*s)); } +uint32_t is_hom_map(uint64_t x, uint32_t len, mc_interval_t *p, uint32_t *p_idx, ha_ug_index* idx) +{ + mc_interval_t *a = NULL; + uint32_t uid, qs, qe, as, ae, occ, k; + uint64_t oLen; + interpr_hit(idx, x, len, &uid, &qs, &qe); + + a = p + p_idx[uid]; + occ = p_idx[uid+1] - p_idx[uid]; + for (k = 0; k < occ; k++) + { + as = a[k].bS; + ae = a[k].bE; + oLen = ((MIN(qe, ae) >= MAX(qs, as))? MIN(qe, ae) - MAX(qs, as) + 1 : 0); + if(oLen == 0) continue; + if(oLen > len*0.2) return 1; + } + return 0; +} + +void update_hits(ha_ug_index* idx, kvec_pe_hit* hits, uint8_t *r_het) +{ + ma_ug_t *ug = idx->ug; + asg_t *rg = idx->read_g; + uint32_t k, v, l, offset, l_pos; + asg_t* nsg = ug->g; + ma_utg_t *u = NULL; + mc_interval_t *t = NULL; + kvec_t(mc_interval_t) p; kv_init(p); + kvec_t(uint32_t) p_idx; kv_init(p_idx); + + kv_push(uint32_t, p_idx, 0); + for (v = 0; v < nsg->n_seq; v++) + { + u = &(ug->u.a[v]); + for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k) + { + if (k == u->n || r_het[u->a[k]>>33] != r_het[u->a[l]>>33]) + { + if(r_het[u->a[l]>>33] == N_HET)///only keep hom suregions + { + kv_pushp(mc_interval_t, p, &t); + t->uID = v; + t->hs = r_het[u->a[l]>>33]; + + t->bS = l_pos; + t->bE = offset + rg->seq[u->a[k-1]>>33].len - 1; + + t->nS = l; + t->nE = k - 1; + } + l = k; + l_pos = offset + (uint32_t)u->a[k-1]; + } + offset += (uint32_t)u->a[k-1]; + } + kv_push(uint32_t, p_idx, p.n); + } + + + for (k = 0; k < hits->a.n; ++k) + { + hits->a.a[k].id = (uint64_t)-1; + if(is_hom_map(hits->a.a[k].s, hits->a.a[k].len>>32, p.a, p_idx.a, idx) || + is_hom_map(hits->a.a[k].e, (uint32_t)hits->a.a[k].len, p.a, p_idx.a, idx)) + { + continue; + } + hits->a.a[k].id = 0; + } + + kv_destroy(p); kv_destroy(p_idx); +} + void verbose_het_stat(bubble_type *bub) { uint64_t i, hetBase = 0, homBase = 0; @@ -14724,15 +15274,15 @@ void debug_output_disconnected_hits(ha_ug_index* idx, kvec_pe_hit *hits, hc_link tn = ((h_a[l].e<<1)>>shif); if(!IF_HOM(tn, *bub) && tn != qn) { - t = &(link->a.a[qn]); - for (i = 0, u_dis = (uint64_t)-1; i < t->e.n; i++) - { - if(t->e.a[i].del || t->e.a[i].uID != tn) continue; - u_dis = (t->e.a[i].dis ==(uint64_t)-1? (uint64_t)-1 : t->e.a[i].dis>>3); - break; - } + // t = &(link->a.a[qn]); + // for (i = 0, u_dis = (uint64_t)-1; i < t->e.n; i++) + // { + // if(t->e.a[i].del || t->e.a[i].uID != tn) continue; + // u_dis = (t->e.a[i].dis ==(uint64_t)-1? (uint64_t)-1 : t->e.a[i].dis>>3); + // break; + // } - if(u_dis == (uint64_t)-1) + // if(u_dis == (uint64_t)-1) { kv_pushp(u_trans_t, k_trans, &p); p->qn = qn; p->tn = tn; p->occ = (k-l); @@ -14770,8 +15320,17 @@ void debug_output_disconnected_hits(ha_ug_index* idx, kvec_pe_hit *hits, hc_link { qn = k_trans.a[i].qn; tn = k_trans.a[i].tn; - fprintf(stderr, "s-utg%.6lul<--->d-utg%.6lul(occ: %u)\n", qn + 1, tn + 1, - ((uint32_t)-1) - k_trans.a[i].occ); + + t = &(link->a.a[qn]); + for (k = 0, u_dis = (uint64_t)-1; k < t->e.n; k++) + { + if(t->e.a[k].del || t->e.a[k].uID != tn) continue; + u_dis = (t->e.a[k].dis ==(uint64_t)-1? (uint64_t)-1 : t->e.a[k].dis>>3); + break; + } + + fprintf(stderr, "s-utg%.6lul<--->d-utg%.6lul(occ: %u):(dis-%lu)\n", qn + 1, tn + 1, + ((uint32_t)-1) - k_trans.a[i].occ, u_dis); } fprintf(stderr, "########hits########\n"); @@ -14788,6 +15347,69 @@ void debug_output_disconnected_hits(ha_ug_index* idx, kvec_pe_hit *hits, hc_link kv_destroy(k_trans); } +void resolve_tangles_hic(ha_ug_index *idx, bubble_type *bub, kvec_pe_hit *hits, kv_u_trans_t *ta) +{ + uint32_t i, k, l, m, h_occ; + uint64_t shif = 64 - idx->uID_bits, qn, tn; + pe_hit *h_a = NULL; + u_trans_t *p = NULL; + + identify_bubbles(idx->ug, bub, idx->t_ch->is_r_het, &(idx->t_ch->k_trans)); + + ta->idx.n = ta->n = 0; + if(hits->idx.n == 0) idx_hc_links(hits, idx, NULL); + + for (qn = 0; qn < hits->idx.n; qn++) + { + h_a = hits->a.a + (hits->idx.a[qn]>>32); + h_occ = (uint32_t)(hits->idx.a[qn]); + + for (k = 1, l = 0; k <= h_occ; ++k) ///same qn + { + if (k == h_occ || ((h_a[k].e<<1)>>shif) != ((h_a[l].e<<1)>>shif)) //same qn and tn + { + tn = ((h_a[l].e<<1)>>shif); + if(tn != qn) + { + kv_pushp(u_trans_t, *ta, &p); + p->qn = qn; p->tn = tn; p->occ = (k-l); + kv_pushp(u_trans_t, *ta, &p); + p->qn = tn; p->tn = qn; p->occ = (k-l); + } + l = k; + } + } + } + + radix_sort_u_trans_m(ta->a, ta->a + ta->n); + + for (k = 1, l = 0, m = 0; k <= ta->n; ++k) + { + if (k == ta->n || ta->a[l].qn != ta->a[k].qn || ta->a[l].tn != ta->a[k].tn) //same qn and tn + { + for (i = l, h_occ = 0; i < k; i++) + { + h_occ += ta->a[i].occ; + } + + qn = ta->a[l].qn; + tn = ta->a[l].tn; + ta->a[m] = ta->a[l]; + ta->a[m].occ = h_occ; + ta->a[m].nw = (h_occ*SCALL)/(MIN(hits->occ.a[qn], hits->occ.a[tn])); + m++; + l = k; + } + } + ta->n = m; + kt_u_trans_t_idx(ta, idx->ug->g->n_seq); + + + resolve_bubble_chain_by_hic(idx, ta, bub); + + ta->idx.n = ta->n = 0; +} + int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) { double index_time = yak_realtime(); @@ -14808,7 +15430,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) alignment_worker_pipeline(&sl, fn1, fn2); write_hc_hits(&sl.hits, asm_opt.output_file_name); } - + // update_hits(idx, &sl.hits, idx->t_ch->is_r_het); ///debug_hc_hits_v14(&sl.hits, asm_opt.output_file_name, sl.idx); ////dedup_hits(&(sl.hits), sl.idx); ///write_hc_hits_v14(&sl.hits, asm_opt.output_file_name); @@ -14826,15 +15448,16 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) bub.round_id = 0; bub.n_round = asm_opt.n_weight; for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++) { - identify_bubbles(idx->ug, &bub, idx->t_ch->is_r_het, &(idx->t_ch->k_trans)); + // identify_bubbles(idx->ug, &bub, idx->t_ch->is_r_het, &(idx->t_ch->k_trans)); if(bub.round_id == 0) { + resolve_tangles_hic(idx, &bub, &sl.hits, &k_trans); measure_distance(idx, idx->ug, &sl.hits, &link, &bub, &(idx->t_ch->k_trans)); } renew_kv_u_trans(&k_trans, &link, &sl.hits, &(idx->t_ch->k_trans), idx, &bub, s->s, 0); // if(bub.round_id == 0) init_phase(idx, &k_trans, &bub, s); - update_trans_g(idx, &k_trans, &bub); + // 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, /**&bub**/NULL, &(idx->t_ch->k_trans)); From 24e5453781db3682cbee9dbaff92090ba71b7c7d Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Sun, 2 May 2021 14:16:36 -0400 Subject: [PATCH 5/5] bug fix --- CommandLines.cpp | 1 - CommandLines.h | 3 ++- Overlaps.cpp | 29 ++++++-------------- Purge_Dups.cpp | 70 +++++++++++++++++++++++++++++++++++++++++++----- Purge_Dups.h | 1 + hic.cpp | 15 +++++++++++ 6 files changed, 90 insertions(+), 29 deletions(-) diff --git a/CommandLines.cpp b/CommandLines.cpp index 85c2f3d..beac8f2 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -175,7 +175,6 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->purge_level_trio = 0; asm_opt->purge_simi_rate_l2 = 0.75; asm_opt->purge_simi_rate_l3 = 0.55; - ///asm_opt->purge_simi_rate_hic = 0.85; asm_opt->purge_overlap_len = 1; ///asm_opt->purge_overlap_len_hic = 50; asm_opt->recover_atg_cov_min = -1024; diff --git a/CommandLines.h b/CommandLines.h index 52107ab..718051e 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.15.1-r333" +#define HA_VERSION "0.15.1-r334" #define VERBOSE 0 @@ -86,6 +86,7 @@ typedef struct { float purge_simi_rate_l2; float purge_simi_rate_l3; float purge_simi_thres; + ///float purge_simi_rate_hic; long long small_pop_bubble_size; diff --git a/Overlaps.cpp b/Overlaps.cpp index 1d4d3da..0c03a76 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -11200,7 +11200,7 @@ uint32_t only_len) break; } } - if(k == nv) fprintf(stderr, "ERROR\n"); + if(k == nv) fprintf(stderr, "ERROR-set_utg_offset\n"); } p_v = v; len += l; @@ -11272,8 +11272,8 @@ void refine_u_trans_t(u_trans_hit_t *q, kv_ca_buf_t* cb) if(si != cb->n && ei != cb->n) break; } - if(si == 0 || ei == 0) fprintf(stderr, "ERROR\n"); - if(si >= cb->n || ei >= cb->n) fprintf(stderr, "ERROR\n"); + if(si == 0 || ei == 0) fprintf(stderr, "ERROR-si-ei-0\n"); + if(si >= cb->n || ei >= cb->n) fprintf(stderr, "ERROR-si-ei-1\n"); si--; ei--; if(s < cb->a[si].c_x_p || ((si + 1) < cb->n && s >= cb->a[si + 1].c_x_p)) @@ -11452,7 +11452,7 @@ uint32_t get_u_trans_hit(u_trans_hit_idx *t, u_trans_hit_t *hit) break; } } - if(k == nv) fprintf(stderr, "ERROR\n"); + if(k == nv) fprintf(stderr, "ERROR-nv\n"); } ///[t->cBeg, t->cEnd) @@ -12762,12 +12762,6 @@ long long gap_fuzz, bub_label_t* b_mask_t) if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_trans_chain(cov->t_ch, output_file_name); } - hic_analysis(ug, sg, cov?cov->t_ch:t_ch); - - - - - char* gfa_name = (char*)malloc(strlen(output_file_name)+25); sprintf(gfa_name, "%s.clean_d_utg.noseq.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); @@ -12775,9 +12769,7 @@ long long gap_fuzz, bub_label_t* b_mask_t) fclose(output_file); free(gfa_name); - - - + hic_analysis(ug, sg, cov?cov->t_ch:t_ch); if(cov) destory_hap_cov_t(&cov); if(t_ch) destory_trans_chain(&t_ch); @@ -16593,13 +16585,13 @@ void chain_origin_trans_uid_s_bubble(buf_t *pri, buf_t* aux, uint32_t beg, uint3 for (i = 0; i < nv; ++i) { if(av[i].del) continue; - if(av[i].v == pri_v) priEnd = pri_len - av[i].ol - 1; - if(av[i].v == aux_v) auxEnd = aux_len - av[i].ol - 1; + 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); } if(priBeg == (uint32_t)-1 || priEnd == (uint32_t)-1 || auxBeg == (uint32_t)-1 || auxEnd == (uint32_t)-1) { - fprintf(stderr, "ERROR\n"); + fprintf(stderr, "ERROR-s_bubble\n"); } cov->u_buffer.a.n = cov->tailIndex.a.n = 0; @@ -27630,11 +27622,6 @@ bub_label_t* b_mask_t) ma_ug_t *ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); - // FILE* output_file = fopen("straw-debug.noseq.gfa", "w"); - // ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); - // fclose(output_file); - - hap_cov_t *cov = NULL; asg_t *copy_sg = copy_read_graph(sg); ma_ug_t *copy_ug = copy_untig_graph(ug); diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index 410c55b..2e81c19 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -11,6 +11,7 @@ #include "rcut.h" KDQ_INIT(uint64_t) +KSORT_INIT_GENERIC(uint64_t) uint8_t debug_enable = 0; @@ -2423,7 +2424,7 @@ uint64_t* position_index, uint32_t xUid, uint32_t yUid, ma_utg_t* xReads, ma_utg ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, ma_sub_t *coverage_cut, float Hap_rate, int is_local, int max_hang, int min_ovlp, uint64_t cov_threshold, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, hap_cov_t *cov, long long* r_x_pos_beg, long long* r_x_pos_end, -long long* r_y_pos_beg, long long* r_y_pos_end) +long long* r_y_pos_beg, long long* r_y_pos_end, float *sim) { uint32_t max_count = 0, min_count = 0, flag; uint32_t xLen = xReads->n, xIndex; @@ -2537,6 +2538,7 @@ long long* r_y_pos_beg, long long* r_y_pos_end) get_pair_hap_similarity_by_base(xReads, read_g, yUid, reverse_sources, ruIndex, *r_x_pos_beg, *r_x_pos_end, &xLeftMatch, &xLeftTotal); + (*sim) = ((double)xLeftMatch)/((double)xLeftTotal); if(xLeftMatch == 0 || xLeftTotal == 0 || xLeftMatch <= xLeftTotal*Hap_rate) { return NON_PLOID; @@ -2842,6 +2844,63 @@ int filter_secondary_chain(long long max_score, long long cur_score, double rate return 1; } + +void filter_secondary_ovlp(kvec_hap_overlaps *x, kvec_t_u64_warp *a, float sim_flt, float ovlp_flt) +{ + if(sim_flt == 0 || ovlp_flt == 0 || x->a.n == 0) return; + #define f_ovlp(s_0, e_0, s_1, e_1) ((MIN((e_0), (e_1)) > MAX((s_0), (s_1)))? MIN((e_0), (e_1)) - MAX((s_0), (s_1)):0) + uint32_t i, m, k; + uint64_t t, ovlp; + hap_overlaps *p = NULL; + a->a.n = 0; + for (i = 0; i < x->a.n; i++) + { + if(x->a.a[i].s < sim_flt) continue; + t = x->a.a[i].x_beg_pos; t<<=32; t |= x->a.a[i].x_end_pos; + kv_push(uint64_t, a->a, t); + } + + if(a->a.n == 0) return; + ks_introsort_uint64_t(a->a.n, a->a.a); + + for (i = m = 1; i < a->a.n; ++i) + { + t = a->a.a[m-1]; + ovlp = f_ovlp(t>>32, (uint32_t)t, a->a.a[i]>>32, (uint32_t)a->a.a[i]); + if(ovlp == 0) + { + a->a.a[m] = a->a.a[i]; + m++; + } + else + { + t = MIN(a->a.a[m-1]>>32, a->a.a[i]>>32); + t<<=32; + t |= MAX((uint32_t)a->a.a[m-1], (uint32_t)a->a.a[i]); + a->a.a[m-1] = t; + } + } + a->a.n = m; + + for (i = m = 0; i < x->a.n; i++) + { + p = &(x->a.a[i]); + if(p->s < sim_flt) + { + for (k = ovlp = 0; k < a->a.n; k++) + { + ovlp += f_ovlp(p->x_beg_pos, p->x_end_pos, a->a.a[k]>>32, (uint32_t)a->a.a[k]); + if(ovlp >= ovlp_flt*(p->x_end_pos-p->x_beg_pos)) break; + } + if(k < a->a.n) continue; + if(ovlp >= ovlp_flt*(p->x_end_pos-p->x_beg_pos)) continue; + } + x->a.a[m] = x->a.a[i]; + m++; + } + x->a.n = m; +} + static void hap_alignment_advance_worker(void *_data, long eid, int tid) { hap_alignment_struct_pip* hap_buf = (hap_alignment_struct_pip*)_data; @@ -2852,7 +2911,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) R_to_U* ruIndex = hap_buf->ruIndex; ma_sub_t *coverage_cut = hap_buf->coverage_cut; uint64_t* position_index = hap_buf->position_index; - float Hap_rate = hap_buf->Hap_rate; + float Hap_rate = hap_buf->Hap_rate/**MIN(hap_buf->Hap_rate, 0.2)**/, sim; int max_hang = hap_buf->max_hang; int min_ovlp = hap_buf->min_ovlp; float chain_rate = hap_buf->chain_rate; @@ -3017,7 +3076,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) if(calculate_pair_hap_similarity_advance(&(u_can->a.a[k]), position_index, xUid, yUid, xReads, yReads, sources, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate, (asm_opt.purge_level_primary<=2? 0:1), max_hang, min_ovlp, cov_threshold, u_buffer, - score_vc, prevIndex_vec, cov, &r_x_pos_beg, &r_x_pos_end, &r_y_pos_beg, &r_y_pos_end)!=PLOID) + score_vc, prevIndex_vec, cov, &r_x_pos_beg, &r_x_pos_end, &r_y_pos_beg, &r_y_pos_end, &sim)!=PLOID) { continue; } @@ -3050,6 +3109,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) hap_align.xUid = xUid; hap_align.yUid = yUid; hap_align.status = SELF_EXIST; + hap_align.s = sim; kv_push(hap_overlaps, all_ovlp->x[hap_align.xUid].a, hap_align); } /** @@ -3091,7 +3151,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) all_ovlp->x[xUid].a.n = m + 1; } } - + // filter_secondary_ovlp(&all_ovlp->x[xUid], u_vecs, hap_buf->Hap_rate, 0.7); } int get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn) @@ -4513,7 +4573,6 @@ void print_all_purge_ovlp(ma_ug_t *ug, hap_overlaps_list* all_ovlp, const char* for (v = 0; v < all_ovlp->num; v++) { uId = v; - ///if(uId != 96 && uId != 272) continue; for (i = 0; i < all_ovlp->x[uId].a.n; i++) { print_hap_paf(ug, &(all_ovlp->x[uId].a.a[i])); @@ -5242,7 +5301,6 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle if(asm_opt.polyploidy <= 2) { mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL); - ///pt_solve(&all_ovlp, cov->t_ch, ug, read_g, 0.8, R_INF.trio_flag); } if(collect_p_trans && collect_p_trans_f == 1) diff --git a/Purge_Dups.h b/Purge_Dups.h index 337f590..5c4f273 100644 --- a/Purge_Dups.h +++ b/Purge_Dups.h @@ -52,6 +52,7 @@ typedef struct { uint32_t yUid; uint32_t weight; long long score; + float s; }hap_overlaps; typedef struct { diff --git a/hic.cpp b/hic.cpp index 99fa1c6..097e7cc 100644 --- a/hic.cpp +++ b/hic.cpp @@ -15410,6 +15410,20 @@ void resolve_tangles_hic(ha_ug_index *idx, bubble_type *bub, kvec_pe_hit *hits, ta->idx.n = ta->n = 0; } + +void print_kv_u_trans_t(kv_u_trans_t *ta) +{ + uint32_t i; + u_trans_t *p = NULL; + for (i = 0; i < ta->n; i++) + { + p = &(ta->a[i]); + fprintf(stderr, "q-utg%.6ul\tqs(%u)\tqe(%u)\tt-utg%.6ul\tts(%u)\tte(%u)\trev(%u)\tw(%f)\tf(%u)\n", + p->qn+1, p->qs, p->qe, p->tn+1, p->ts, p->te, p->rev, p->nw, p->f); + } + fprintf(stderr, "[M::%s::] \n", __func__); +} + int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) { double index_time = yak_realtime(); @@ -15434,6 +15448,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) ///debug_hc_hits_v14(&sl.hits, asm_opt.output_file_name, sl.idx); ////dedup_hits(&(sl.hits), sl.idx); ///write_hc_hits_v14(&sl.hits, asm_opt.output_file_name); + // print_kv_u_trans_t(&(idx->t_ch->k_trans)); hc_links link; init_hc_links(&link, idx->ug->g->n_seq, idx->t_ch);