From 4105dd360e6bc2f217d9e1cb18555f7eda5c4ef2 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Tue, 11 May 2021 13:01:06 -0400 Subject: [PATCH] correction --- Overlaps.cpp | 3 +- Overlaps.h | 7 +- hic.cpp | 54 +++- hic.h | 1 + horder.cpp | 686 ++++++++++++++++++++++++++++++++++++++++++++++++--- horder.h | 1 + 6 files changed, 711 insertions(+), 41 deletions(-) diff --git a/Overlaps.cpp b/Overlaps.cpp index e72746d..549b746 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -29580,8 +29580,7 @@ void clean_sg_by_utg(asg_t *sg, ma_ug_t *ug) } } } - - + asg_cleanup(sg); /*******************************for debug************************************/ diff --git a/Overlaps.h b/Overlaps.h index deb8482..c9fb08c 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -4,6 +4,7 @@ #include #include "kvec.h" #include "kdq.h" +#include "ksort.h" ///#define MIN_OVERLAP_LEN 2000 ///#define MIN_OVERLAP_LEN 500 @@ -320,11 +321,10 @@ static inline asg_arc_t *asg_arc_pushp(asg_t *g) // set asg_arc_t::del for v->w static inline void asg_arc_del(asg_t *g, uint32_t v, uint32_t w, int del) { - uint32_t i, nv = asg_arc_n(g, v)/**, found = 0**/; + uint32_t i, nv = asg_arc_n(g, v); asg_arc_t *av = asg_arc_a(g, v); for (i = 0; i < nv; ++i) - if (av[i].v == w) av[i].del = !!del/**, found = 1**/; - /**if(found == 0) fprintf(stderr, "ERROR\n");**/ + if (av[i].v == w) av[i].del = !!del; } // set asg_arc_t::del and asg_seq_t::del to 1 for sequence s and all its associated arcs @@ -1217,6 +1217,7 @@ ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t); +uint32_t cmp_untig_graph(ma_ug_t *src, ma_ug_t *dest); #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/hic.cpp b/hic.cpp index 776e718..d7f8840 100644 --- a/hic.cpp +++ b/hic.cpp @@ -15151,6 +15151,7 @@ ps_t* init_ps_t(uint64_t seed, uint64_t n) ps_t *s = NULL; CALLOC(s, 1); s->xs = seed; CALLOC(s->s, n); + s->n = n; return s; } @@ -15424,6 +15425,39 @@ void print_kv_u_trans_t(kv_u_trans_t *ta) fprintf(stderr, "[M::%s::] \n", __func__); } +void write_ps_t(ps_t *s, const char *fn) +{ + char *buf = (char*)calloc(strlen(fn) + 25, 1); + sprintf(buf, "%s.hic.pst.bin", fn); + FILE* fp = fopen(buf, "w"); + + fwrite(&(s->xs), sizeof(s->xs), 1, fp); + fwrite(&(s->n), sizeof(s->n), 1, fp); + fwrite(s->s, sizeof(int8_t), s->n, fp); + + fclose(fp); + free(buf); +} + +int load_ps_t(ps_t **s, const char *fn) +{ + uint64_t flag = 0; + char *buf = (char*)calloc(strlen(fn) + 25, 1); + sprintf(buf, "%s.hic.pst.bin", fn); + FILE* fp = NULL; + fp = fopen(buf, "r"); + if(!fp) return 0; + CALLOC(*s, 1); + flag += fread(&((*s)->xs), sizeof((*s)->xs), 1, fp); + flag += fread(&((*s)->n), sizeof((*s)->n), 1, fp); + MALLOC((*s)->s, (*s)->n); + flag += fread((*s)->s, sizeof(int8_t), (*s)->n, fp); + + fclose(fp); + free(buf); + return 1; +} + int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_opt_t *opt) { double index_time = yak_realtime(); @@ -15456,20 +15490,24 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o bubble_type bub; kv_u_trans_t k_trans; kv_init(k_trans); kv_init(k_trans.idx); - ps_t *s = init_ps_t(11, idx->ug->g->n_seq); + ps_t *s = NULL; mb_nodes_t u; kv_init(u.bid); kv_init(u.idx); kv_init(u.u); memset(&bub, 0, sizeof(bubble_type)); bub.round_id = 0; bub.n_round = asm_opt.n_weight; + + resolve_tangles_hic(idx, &bub, &sl.hits, &k_trans); + measure_distance(idx, idx->ug, &sl.hits, &link, &bub, &(idx->t_ch->k_trans)); + if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_ps_t(&s, asm_opt.output_file_name)) + { + bub.round_id = bub.n_round; + label_unitigs_sm(s->s, idx->ug); + goto skip_flipping; + } + s = init_ps_t(11, idx->ug->g->n_seq); 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)); - 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); @@ -15493,7 +15531,9 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o label_unitigs(&(hap.g_p), idx->ug); **/ } + write_ps_t(s, asm_opt.output_file_name); + skip_flipping: verbose_het_stat(&bub); horder_t *ho = init_horder_t(&sl.hits, idx->uID_bits, idx->pos_mode, idx->read_g, idx->ug, &bub, opt); diff --git a/hic.h b/hic.h index 6a197b0..8f14dfc 100644 --- a/hic.h +++ b/hic.h @@ -53,6 +53,7 @@ typedef struct { typedef struct { int8_t *s; uint64_t xs; + uint64_t n; } ps_t; typedef struct { diff --git a/horder.cpp b/horder.cpp index b1939da..ddcce27 100644 --- a/horder.cpp +++ b/horder.cpp @@ -21,10 +21,30 @@ KDQ_INIT(uint64_t) KRADIX_SORT_INIT(pe_hit_idx_hn1, pe_hit, pe_hit_an1_idx_key, member_size(pe_hit, s)) #define pe_hit_an2_idx_key(x) ((x).e<<1) KRADIX_SORT_INIT(pe_hit_idx_hn2, pe_hit, pe_hit_an2_idx_key, member_size(pe_hit, e)) +#define generic_key(x) (x) +KRADIX_SORT_INIT(ho64, uint64_t, generic_key, 8) +#define get_hit_srev(x, k) ((x).a.a[(k)].s>>63) +#define get_hit_slen(x, k) ((x).a.a[(k)].len>>32) #define get_hit_suid(x, k) (((x).a.a[(k)].s<<1)>>(64 - (x).uID_bits)) #define get_hit_spos(x, k) ((x).a.a[(k)].s & (x).pos_mode) +#define get_hit_spos_e(x, k) (get_hit_srev((x),(k))?\ + ((get_hit_spos((x),(k))+1>=get_hit_slen((x),(k)))?\ + (get_hit_spos((x),(k))+1-get_hit_slen((x),(k))):0)\ + :(get_hit_spos((x),(k))+get_hit_slen((x),(k))-1)) + +#define get_hit_erev(x, k) ((x).a.a[(k)].e>>63) +#define get_hit_elen(x, k) ((uint32_t)((x).a.a[(k)].len)) #define get_hit_euid(x, k) (((x).a.a[(k)].e<<1)>>(64 - (x).uID_bits)) #define get_hit_epos(x, k) ((x).a.a[(k)].e & (x).pos_mode) +#define get_hit_epos_e(x, k) (get_hit_erev((x),(k))?\ + ((get_hit_epos((x),(k))+1>=get_hit_elen((x),(k)))?\ + (get_hit_epos((x),(k))+1-get_hit_elen((x),(k))):0)\ + :(get_hit_epos((x),(k))+get_hit_elen((x),(k))-1)) +#define OVL(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) +#define BREAK_THRES 1000000 +#define BREAK_CUTOFF 0.1 +#define BREAK_BOUNDARY 0.01 +#define GAP_LEN 100 typedef struct { uint64_t ruid; @@ -37,10 +57,47 @@ typedef struct { kvec_t(uint64_t) idx; } u_hits_t; +typedef struct { + uint64_t s, e, dp; +} h_cov_t; +typedef struct { + h_cov_t *a; + size_t n, m; +} h_covs; + +#define h_cov_s_key(x) ((x).s) +KRADIX_SORT_INIT(h_cov_s, h_cov_t, h_cov_s_key, member_size(h_cov_t, s)) #define hit_aux_ruid_key(x) ((x).ruid) KRADIX_SORT_INIT(hit_aux_ruid, hit_aux_t, hit_aux_ruid_key, member_size(hit_aux_t, ruid)) +void print_N50(ma_ug_t* ug) +{ + kvec_t(uint64_t) b; kv_init(b); + uint64_t i, s, len; + for (i = s = 0; i < ug->g->n_seq; ++i) + { + kv_push(uint64_t, b, ug->u.a[i].len); + s += ug->u.a[i].len; + } + len = s; + fprintf(stderr, "[M::%s::] Genome Size: %lu, # Contigs: %u\n", __func__, len, ug->g->n_seq); + + radix_sort_ho64(b.a, b.a+b.n); + i = b.n; s = 0; + while (i > 0) + { + i--; + s += b.a[i]; + if(s >= (len>>1)) + { + fprintf(stderr, "[M::%s::] N50: %lu\n", __func__, b.a[i]); + break; + } + } + + kv_destroy(b); +} void get_r_hits(kvec_pe_hit *u_hits, kvec_pe_hit *r_hits, asg_t* r_g, ma_ug_t* ug, bubble_type* bub, uint64_t uID_bits, uint64_t pos_mode) { @@ -164,6 +221,26 @@ uint64_t get_corresp_usite(uint64_t rid, uint64_t rpos, uint64_t rev, uint64_t r return a_n; } +void idx_hits(kvec_pe_hit* hits, uint64_t n) +{ + uint64_t k, l; + kv_resize(uint64_t, hits->idx, n); + hits->idx.n = n; + memset(hits->idx.a, 0, hits->idx.n*sizeof(uint64_t)); + + radix_sort_pe_hit_idx_hn1(hits->a.a, hits->a.a + hits->a.n); + for (k = 1, l = 0; k <= hits->a.n; ++k) + { + if (k == hits->a.n || (get_hit_suid(*hits, k) != get_hit_suid(*hits, l))) + { + if (k - l > 1) radix_sort_pe_hit_idx_hn2(hits->a.a + l, hits->a.a + k); + + hits->idx.a[get_hit_suid(*hits, l)] = (uint64_t)l << 32 | (k - l); + l = k; + } + } +} + void update_u_hits(kvec_pe_hit *u_hits, kvec_pe_hit *r_hits, ma_ug_t* ug, asg_t* r_g) { u_hits_t x; memset(&x, 0, sizeof(x)); @@ -172,17 +249,24 @@ void update_u_hits(kvec_pe_hit *u_hits, kvec_pe_hit *r_hits, ma_ug_t* ug, asg_t* pe_hit *t = NULL; uint64_t v, i, l, k, offset, occ_1, occ_2, *a_1, *a_2, i_1, i_2; - for (v = 0; v < ug->g->n_seq; v++) + for (v = 0; v < ug->u.n; v++) { u = &(ug->u.a[v]); for (i = offset = 0; i < u->n; i++) { - kv_pushp(hit_aux_t, x, &p); - p->ruid = u->a[i]>>32; - p->ruid <<= 32; - p->ruid |= v; - p->off = offset; - offset += (uint32_t)u->a[i]; + if(u->a[i] != (uint64_t)-1) + { + kv_pushp(hit_aux_t, x, &p); + p->ruid = u->a[i]>>32; + p->ruid <<= 32; + p->ruid |= v; + p->off = offset; + offset += (uint32_t)u->a[i]; + } + else + { + offset += GAP_LEN; + } } } @@ -227,8 +311,8 @@ void update_u_hits(kvec_pe_hit *u_hits, kvec_pe_hit *r_hits, ma_ug_t* ug, asg_t* } } - free(x.a); free(x.idx.a); kv_destroy(buf.a); + idx_hits(u_hits, ug->u.n); } ma_ug_t* get_trio_unitig_graph(asg_t *sg, uint8_t flag, ug_opt_t *opt) @@ -248,12 +332,549 @@ ma_ug_t* get_trio_unitig_graph(asg_t *sg, uint8_t flag, ug_opt_t *opt) return ug; } +static inline void asg_arc_unique_del(asg_t *g, uint32_t v, uint32_t w, int del) +{ + uint32_t i, nv = asg_arc_n(g, v); + asg_arc_t *av = asg_arc_a(g, v); + for (i = 0; i < nv; ++i) + { + if (av[i].v == w) + { + av[i].del = !!del; + break; + } + } +} +void horder_clean_sg_by_utg(asg_t *sg, ma_ug_t *ug) +{ + uint32_t i, v, n_vx, w, k, m, nv, vx, wx; + asg_arc_t *av = NULL; + ma_utg_t *u = NULL; + + n_vx = sg->n_seq<<1; + for (v = 0; v < n_vx; v++) + { + nv = asg_arc_n(sg, v); + av = asg_arc_a(sg, v); + for (m = 0; m < nv; m++) av[m].del = (!!1); + sg->seq[v>>1].del = (!!1); + } + + for (i = 0; i < ug->g->n_seq; ++i) + { + if(ug->g->seq[i].del) continue; + u = &(ug->u.a[i]); + for (k = 0; (k + 1) < u->n; k++) + { + v = u->a[k]>>32; w = u->a[k+1]>>32; + + asg_arc_unique_del(sg, v, w, 0); + asg_arc_unique_del(sg, w^1, v^1, 0); + } + for (k = 0; k < u->n; k++) + { + sg->seq[u->a[k]>>33].del = (!!0); + } + + v = i<<1; + nv = asg_arc_n(ug->g, v); av = asg_arc_a(ug->g, v); + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + w = av[k].v; + + 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_unique_del(sg, vx, wx, 0); asg_arc_unique_del(sg, wx^1, vx^1, 0); + } + + v = (i<<1)+1; + nv = asg_arc_n(ug->g, v); av = asg_arc_a(ug->g, v); + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + w = av[k].v; + + 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_unique_del(sg, vx, wx, 0); asg_arc_unique_del(sg, wx^1, vx^1, 0); + } + } + asg_cleanup(sg); + + /*******************************for debug************************************/ + // ma_ug_t *dbg = ma_ug_gen(sg); + // print_N50(dbg); + // print_N50(ug); + + // uint8_t *end = NULL; CALLOC(end, sg->n_seq<<1); + // for (i = 0; i < dbg->g->n_seq; ++i) + // { + // u = &(dbg->u.a[i]); + // if(u->n == 0) continue; + // end[(u->a[0]>>32)^1] = 1; + // end[u->a[u->n-1]>>32] = 2; + // } + + // for (i = 0; i < ug->g->n_seq; ++i) + // { + // u = &(ug->u.a[i]); + // if(u->n == 0) continue; + // for (k = 1; (k + 1) < u->n; k++) + // { + // if(end[(u->a[k]>>32)]) + // { + // fprintf(stderr, "(1) node-%lu, v-%lu, w-%lu, sg(v).n: %u\n", + // u->a[k]>>33, (u->a[k]>>32), (u->a[k+1]>>32), asg_arc_n(sg, (u->a[k]>>32))); + // } + + // if(end[(u->a[k]>>32)^1]) + // { + // fprintf(stderr, "(2) node-%lu, v-%lu, w-%lu, sg(v).n: %u\n", + // u->a[k]>>33, (u->a[k]>>32)^1, (u->a[k-1]>>32)^1, asg_arc_n(sg, (u->a[k]>>32)^1)); + // } + // } + // } + // free(end); + // ma_ug_destroy(dbg); + /*******************************for debug************************************/ +} + +uint64_t get_hic_cov_interval(uint64_t *b, uint64_t b_n, int64_t min_dp, int64_t *boundS, int64_t *boundE, +h_covs *res) +{ + if(res) res->n = 0; + if(min_dp == 0 || b_n == 0) return (uint64_t)-1; + uint64_t i, len = 0; + int64_t dp, old_dp, start = 0, bs = b[0]>>1, be = b[b_n-1]>>1, olen; + h_cov_t *p = NULL; + if(boundS) bs = (*boundS); + if(boundE) be = (*boundE); + for (i = 0, dp = 0, start = 0; i < b_n; ++i) + { + old_dp = dp; + ///if a[j] is qe + if (b[i]&1) --dp; + else ++dp; + + if (old_dp < min_dp && dp >= min_dp) ///old_dp < dp, b.a[j] is qs + { + ///case 2, a[j] is qs + start = b[i]>>1; + } + else if (old_dp >= min_dp && dp < min_dp) ///old_dp > min_dp, b.a[j] is qe + { + olen = OVL(start, (int64_t)(b[i]>>1), bs, be); + if(olen == 0) continue; + if(res) + { + kv_pushp(h_cov_t, *res, &p); + p->s = MAX(start, bs); + p->e = MIN((int64_t)(b[i]>>1), be); + p->dp = old_dp; + } + len += olen; + } + } + return len; +} + +void get_hic_breakpoint(uint64_t *b, uint64_t b_n, int64_t cutoff, h_covs *res, +int64_t cov_s_pos, int64_t cov_e_pos, uint64_t *s, uint64_t *e) +{ + uint64_t i; + (*s) = (*e) = (uint64_t)-1; + res->n = 0; + get_hic_cov_interval(b, b_n, cutoff, &cov_s_pos, &cov_e_pos, res); + if(res->n == 0) return; + int64_t max = -1, max_cur = 0; + int64_t max_s_idx, max_e_idx, cur_s_idx; + for (i = 0; i < res->n; i++)//all intervals have cov >= cutoff + { + if(i > 0 && (res->a[i].s - res->a[i-1].e) > 0)//cov < cutoff + { + max_cur += (res->a[i].s - res->a[i-1].e);//at least positive + if(max < max_cur) + { + max_s_idx = (max < 0? res->a[i-1].e:cur_s_idx); + max_e_idx = res->a[i].s; + max = max_cur; + cur_s_idx = max_s_idx; + } + } + + //cov >= cutoff + max_cur -= (res->a[i].e - res->a[i].s); + if(max_cur < 0) + { + max_cur = 0; + cur_s_idx = res->a[i].e; + } + } + + if(max > 0) + { + (*s) = max_s_idx; + (*e) = max_e_idx; + } +} + +void get_consensus_break(h_covs *res, h_covs *tmp) +{ + uint64_t i, k, n, m, max_cut = 0; + h_cov_t *p = NULL; + tmp->n = 0; + if(res->n == 0) return; + n = res->n; + for (i = 0; i < n; i++) + { + res->a[i].s <<= 1; + if(max_cut < res->a[i].dp) max_cut = res->a[i].dp; + kv_pushp(h_cov_t, *res, &p); + *p = res->a[i]; + p->s = (res->a[i].e<<1)|1; + } + radix_sort_h_cov_s(res->a, res->a+res->n); + int64_t dp, old_dp, start = 0, max_dp; + p = NULL; max_dp = -1; tmp->n = 0; + for (i = 0, dp = 0, start = 0; i < res->n; ++i) + { + old_dp = dp; + ///if a[j] is qe + if (res->a[i].s&1) --dp; + else ++dp; + + if (old_dp < dp) ///old_dp < dp, b.a[j] is qs + { + ///case 2, a[j] is qs + start = res->a[i].s>>1; + } + else if (old_dp > dp) ///old_dp > min_dp, b.a[j] is qe + { + if(max_dp < old_dp) + { + max_dp = old_dp; + tmp->n = 0; + kv_pushp(h_cov_t, *tmp, &p); + p->s = start; p->e = res->a[i].s>>1; p->dp = max_cut; + } + else if(max_dp == old_dp) + { + kv_pushp(h_cov_t, *tmp, &p); + p->s = start; p->e = res->a[i].s>>1; p->dp = max_cut; + } + } + } + if(tmp->n == 0) fprintf(stderr, "ERROR-break-0\n"); + if(tmp->n == 1) return; + + for (i = m = 0; i < res->n; ++i) + { + if(res->a[i].s&1) continue; + res->a[m] = res->a[i]; + res->a[m].s >>= 1; + m++; + } + res->n = m; + if(res->n != n) fprintf(stderr, "ERROR-break-1\n"); + for (i = 0; i < res->n; ++i) + { + for (k = 0; k < tmp->n; k++) + { + if(OVL(res->a[i].s, res->a[i].e, tmp->a[k].s, tmp->a[k].e) == 0) continue; + tmp->a[k].dp = MIN(tmp->a[k].dp, res->a[i].dp); + if(max_cut > tmp->a[k].dp) max_cut = tmp->a[k].dp; + } + } + + for (k = m = 0; k < tmp->n; k++) + { + if(max_cut != tmp->a[k].dp) continue; + tmp->a[m] = tmp->a[k]; + m++; + } + tmp->n = m; + if(tmp->n == 0) fprintf(stderr, "ERROR-break-2\n"); +} + +int64_t update_r_break(uint64_t rs, uint64_t re, h_covs *hits) +{ + uint64_t i, hs, he; + int64_t dp, old_dp, max_dp = 0; + for (i = 0, dp = 0; i < hits->n; ++i) + { + if(hits->a[i].s&1) + { + hs = hits->a[i].e; + he = hits->a[i].s>>1; + } + else + { + hs = hits->a[i].s>>1; + he = hits->a[i].e; + } + if(hs <= rs && he >= re) + { + old_dp = dp; + ///if a[j] is qe + if (hits->a[i].s&1) --dp; + else ++dp; + ///hits->a[i].s is qe + if (old_dp > dp && max_dp < old_dp) + { + max_dp = old_dp; + } + } + } + + return max_dp; +} + +void get_read_breaks(ma_utg_t *u, asg_t* r_g, h_covs *cov, h_covs *hit_tmp, +kvec_pe_hit *hits, uint64_t sidx, uint64_t eidx, uint64_t ulen, uint64_t *idx, uint64_t *rdp) +{ + (*idx) = (*rdp) = (uint64_t)-1; + uint64_t i, k, offset, beg, end, n = cov->n, min_ovlp, o, mi, dp, min_dp; + uint64_t p0s, p0e, p1s, p1e, span_s, span_e; + h_cov_t *p = NULL, *a = NULL; + for (i = offset = 0; i < u->n; i++) + { + end = offset + ((u->a[i] != (uint64_t)-1?r_g->seq[u->a[i]>>33].len:GAP_LEN)); + offset += (u->a[i] != (uint64_t)-1? (uint32_t)u->a[i]:GAP_LEN); + beg = offset; + + if(u->a[i] == (uint64_t)-1) beg -= GAP_LEN; + if(end <= beg && i + 1 < u->n && u->a[i+1] == (uint64_t)-1) + { + end = beg + GAP_LEN; + } + + for (k = 0; k < n; k++) + { + if(OVL(beg, end, cov->a[k].s, cov->a[k].e) == 0) continue; + kv_pushp(h_cov_t, *cov, &p); + p->s = beg; p->e = end; p->dp = i; + break; + } + } + a = cov->a + n; + n = cov->n - n; + if(n == 0) fprintf(stderr, "ERROR-r-break\n"); + hit_tmp->n = 0; + for (i = sidx; i < eidx; i++)///keep all hic hits that contain interval we want + { + if(get_hit_suid(*hits, i) != get_hit_euid(*hits, i)) continue; + p0s = get_hit_spos(*hits, i); + p0e = get_hit_spos_e(*hits, i); + p1s = get_hit_epos(*hits, i); + p1e = get_hit_epos_e(*hits, i); + + span_s = MIN(MIN(p0s, p0e), MIN(p1s, p1e)); + span_s = MIN(span_s, ulen-1); + span_e = MAX(MAX(p0s, p0e), MAX(p1s, p1e)); + span_e = MIN(span_e, ulen-1) + 1; + + //if(span_e - span_s <= ulen*BREAK_CUTOFF)//need it or not? + { + for (k = 0; k < n; k++) + { + if(span_s <= a[k].s && span_e >= a[k].e) break; + } + if(k >= n) continue; + kv_pushp(h_cov_t, *hit_tmp, &p); + p->s = (span_s<<1); p->e = span_e; + kv_pushp(h_cov_t, *hit_tmp, &p); + p->s = ((span_e<<1)|1); p->e = span_s; + } + } + + radix_sort_h_cov_s(hit_tmp->a, hit_tmp->a+hit_tmp->n); + for (k = 0, min_dp = (uint64_t)-1; k < n; k++) + { + //a[k].s, a[k].e + dp = update_r_break(a[k].s, a[k].e, hit_tmp); + a[k].dp = (uint32_t)a[k].dp; + a[k].dp += (dp << 32); + if(dp < min_dp) min_dp = dp; + } + + min_ovlp = mi = (uint64_t)-1; + for (k = 0; k < n; k++) + { + if((a[k].dp>>32) != min_dp) continue; + i = (uint32_t)a[k].dp; + o = 0; + if(u->a[i] != (uint64_t)-1) + { + o = r_g->seq[u->a[i]>>33].len - (uint32_t)u->a[i]; + } + if(o < min_ovlp) + { + min_ovlp = o; + mi = i; + } + } + if(mi != (uint64_t)-1) (*idx) = mi, (*rdp) = min_dp; +} + +void debug_sub_cov(kvec_pe_hit *hits, uint64_t sidx, uint64_t eidx, uint64_t ulen, ma_utg_t *u, asg_t* r_g, +uint64_t rid, uint64_t i_cnt) +{ + uint64_t i, offset, beg, end, rs, re; + uint64_t p0s, p0e, p1s, p1e, span_s, span_e, cnt = 0; + rs = re = (uint64_t)-1; + for (i = offset = 0; i < u->n; i++) + { + end = offset + ((u->a[i] != (uint64_t)-1?r_g->seq[u->a[i]>>33].len:GAP_LEN)); + offset += (u->a[i] != (uint64_t)-1? (uint32_t)u->a[i]:GAP_LEN); + beg = offset; + + if(u->a[i] == (uint64_t)-1) beg -= GAP_LEN; + if(end <= beg && i + 1 < u->n && u->a[i+1] == (uint64_t)-1) + { + end = beg + GAP_LEN; + } + if(rid == i) + { + rs = beg; + re = end; + } + } + + for (i = sidx; i < eidx; i++)///keep all hic hits that contain interval we want + { + if(get_hit_suid(*hits, i) != get_hit_euid(*hits, i)) continue; + p0s = get_hit_spos(*hits, i); + p0e = get_hit_spos_e(*hits, i); + p1s = get_hit_epos(*hits, i); + p1e = get_hit_epos_e(*hits, i); + + span_s = MIN(MIN(p0s, p0e), MIN(p1s, p1e)); + span_s = MIN(span_s, ulen-1); + span_e = MAX(MAX(p0s, p0e), MAX(p1s, p1e)); + span_e = MIN(span_e, ulen-1) + 1; + + //if(span_e - span_s <= ulen*BREAK_CUTOFF)//need it or not? + { + if(span_s <= rs && span_e >= re) cnt++; + } + } + + // if(cnt != i_cnt) fprintf(stderr, "cnt-%lu, i_cnt-%lu\n", cnt, i_cnt); + fprintf(stderr, "******cnt-%lu, i_cnt-%lu\n", cnt, i_cnt); +} + +void break_utg_horder(horder_t *h, h_covs *b_points) +{ + uint64_t i, rid, uid; + for (i = 0; i < b_points->n; i++) + { + rid = b_points->a[i].s; + uid = b_points->a[i].e; + } + +} + +void break_contig_init(horder_t *h, uint64_t cutoff_s, uint64_t cutoff_e) +{ + uint64_t k, l, i, p0s, p0e, p1s, p1e, ulen, cov_hic, cov_utg, cov_ava, span_s, span_e, cutoff, bs, be, dp; + kvec_t(uint64_t) b; kv_init(b); + h_covs cov_buf; kv_init(cov_buf); + h_covs res; kv_init(res); + h_covs b_points; kv_init(b_points); + h_cov_t *p = NULL; + ma_ug_t *ug = h->ug; + kvec_pe_hit *hits = &(h->u_hits); + b_points.n = 0; + for (k = 1, l = 0; k <= hits->a.n; ++k) + { + if (k == hits->a.n || (get_hit_suid(*hits, k) != get_hit_suid(*hits, l))) + { + ulen = ug->u.a[get_hit_suid(*hits, l)].len; + b.n = 0; cov_hic = cov_utg = 0; + if(ulen >= BREAK_THRES) + { + for (i = l; i < k; i++) + { + if(get_hit_suid(*hits, i) != get_hit_euid(*hits, i)) continue; + p0s = get_hit_spos(*hits, i); + p0e = get_hit_spos_e(*hits, i); + p1s = get_hit_epos(*hits, i); + p1e = get_hit_epos_e(*hits, i); + + span_s = MIN(MIN(p0s, p0e), MIN(p1s, p1e)); + span_s = MIN(span_s, ulen-1); + span_e = MAX(MAX(p0s, p0e), MAX(p1s, p1e)); + span_e = MIN(span_e, ulen-1) + 1; + //if(span_e - span_s <= ulen*BREAK_CUTOFF)//need it or not? + { + kv_push(uint64_t, b, (span_s<<1)); + kv_push(uint64_t, b, (span_e<<1)|1); + cov_hic += (span_e - span_s); + } + } + + radix_sort_ho64(b.a, b.a+b.n); + cov_utg = get_hic_cov_interval(b.a, b.n, 1, NULL, NULL, NULL); + cov_ava = (cov_utg? cov_hic/cov_utg:0); + ///if cov_ava == 0, do nothing or break? + + fprintf(stderr, "\n[M::%s::] utg%.6lul, ulen: %lu, # hic hits: %lu, map cov: %lu, utg cov: %lu, average: %lu\n", + __func__, get_hit_suid(*hits, l)+1, ulen, (uint64_t)(b.n>>1), cov_hic, cov_utg, cov_ava); + + + res.n = 0; + for (i = cutoff_s; i <= cutoff_e; i++) + { + if(i == 0) continue; + cutoff = cov_ava/i; + if(cutoff == 0) continue; + get_hic_breakpoint(b.a, b.n, cutoff, &cov_buf, ulen*BREAK_BOUNDARY, ulen - ulen*BREAK_BOUNDARY, &bs, &be); + if(bs != (uint64_t)-1 && be != (uint64_t)-1) + { + kv_pushp(h_cov_t, res, &p); + p->s = bs; p->e = be; p->dp = cutoff; + fprintf(stderr, "cutoff: %lu, bs: %lu, be: %lu\n", cutoff, bs, be); + } + } + + if(res.n > 0) + { + get_consensus_break(&res, &cov_buf); + for (i = 0; i < cov_buf.n; i++) + { + fprintf(stderr, "consensus_break-s: %lu, e: %lu\n", cov_buf.a[i].s, cov_buf.a[i].e); + } + get_read_breaks(&(ug->u.a[get_hit_suid(*hits, l)]), h->r_g, &cov_buf, + &res, hits, l, k, ulen, &bs, &dp); + if(bs == (uint64_t)-1) fprintf(stderr, "ERROR-read\n"); + + kv_pushp(h_cov_t, b_points, &p); + p->s = bs; p->e = get_hit_suid(*hits, l); p->dp = dp; + fprintf(stderr, "consensus_break-rid: %lu, cov: %lu\n", bs, dp); + // debug_sub_cov(hits, l, k, ulen, &(ug->u.a[get_hit_suid(*hits, l)]), h->r_g, bs, dp); + } + + } + l = k; + } + } + + + + kv_destroy(b); + kv_destroy(cov_buf); + kv_destroy(res); + kv_destroy(b_points); +} + + void generate_haplotypes(horder_t *h, ug_opt_t *opt) { uint64_t i, off; ma_ug_t *ug_1 = NULL, *ug_2 = NULL; - asg_cleanup(h->r_g); - asg_arc_del_trans(h->r_g, asm_opt.gap_fuzz);///must + // asg_cleanup(h->r_g); + // asg_arc_del_trans(h->r_g, asm_opt.gap_fuzz);///must ug_1 = get_trio_unitig_graph(h->r_g, FATHER, opt); ug_2 = get_trio_unitig_graph(h->r_g, MOTHER, opt); @@ -296,28 +917,30 @@ void generate_haplotypes(horder_t *h, ug_opt_t *opt) h->ug->g->seq[i].c = (i < off? FATHER:MOTHER); } /*******************************for debug************************************/ - if(h->ug->g->n_seq != h->ug->u.n) - { - fprintf(stderr, "ERROR-non-equal-length\n"); - } - fprintf(stderr, "h->ug->u.n-%u, off-%lu, ug_2->u.n-%u\n", (uint32_t)h->ug->u.n, off, (uint32_t)ug_2->u.n); - for (i = 0; i < h->ug->g->n_seq; i++) - { - if(h->ug->g->seq[i].len != h->ug->u.a[i].len) - { - fprintf(stderr, "****ERROR-1-%lu: g->seq[i].len-%u, u.a[i].len-%u\n", - i, (uint32_t)h->ug->g->seq[i].len, (uint32_t)h->ug->u.a[i].len); - } - pu = i < off? (&ug_1->u.a[i]) : (&ug_2->u.a[i - off]); - if(h->ug->g->seq[i].len != pu->len) - { - fprintf(stderr, "ERROR-2\n"); - } - } + // if(h->ug->g->n_seq != h->ug->u.n) + // { + // fprintf(stderr, "ERROR-non-equal-length\n"); + // } + // fprintf(stderr, "h->ug->u.n-%u, off-%lu, ug_2->u.n-%u\n", (uint32_t)h->ug->u.n, off, (uint32_t)ug_2->u.n); + // for (i = 0; i < h->ug->g->n_seq; i++) + // { + // if(h->ug->g->seq[i].len != h->ug->u.a[i].len) + // { + // fprintf(stderr, "****ERROR-1-%lu: g->seq[i].len-%u, u.a[i].len-%u\n", + // i, (uint32_t)h->ug->g->seq[i].len, (uint32_t)h->ug->u.a[i].len); + // } + // pu = i < off? (&ug_1->u.a[i]) : (&ug_2->u.a[i - off]); + // if(h->ug->g->seq[i].len != pu->len) + // { + // fprintf(stderr, "ERROR-2\n"); + // } + // } /*******************************for debug************************************/ ug_1 = NULL; ma_ug_destroy(ug_2); + asg_destroy(h->ug->g); + h->ug->g = NULL; } horder_t *init_horder_t(kvec_pe_hit *i_hits, uint64_t i_hits_uid_bits, uint64_t i_hits_pos_mode, @@ -326,8 +949,11 @@ asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt) horder_t *h = NULL; CALLOC(h, 1); get_r_hits(i_hits, &(h->r_hits), i_rg, i_ug, bub, i_hits_uid_bits, i_hits_pos_mode); h->r_g = copy_read_graph(i_rg); + horder_clean_sg_by_utg(h->r_g, i_ug); generate_haplotypes(h, opt); - update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); + print_N50(h->ug); + update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); + break_contig_init(h, 10, 20); return h; } @@ -342,6 +968,8 @@ void destory_horder_t(horder_t **h) kv_destroy((*h)->u_hits.idx); kv_destroy((*h)->u_hits.occ); + kv_destroy((*h)->hp); + ma_ug_destroy((*h)->ug); asg_destroy((*h)->r_g); free((*h)); diff --git a/horder.h b/horder.h index f71738b..22bf244 100644 --- a/horder.h +++ b/horder.h @@ -6,6 +6,7 @@ typedef struct { kvec_pe_hit r_hits, u_hits; ma_ug_t* ug; asg_t* r_g; + kvec_t(uint8_t) hp; }horder_t; horder_t *init_horder_t(kvec_pe_hit *i_hits, uint64_t i_hits_uid_bits, uint64_t i_hits_pos_mode,