diff --git a/Overlaps.cpp b/Overlaps.cpp index 8b9aa7a..da3000b 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -11587,14 +11587,7 @@ uint32_t tn, kv_u_trans_hit_t* ktb, uint32_t bn) } } -typedef struct {///[cBeg, cEnd) - uint32_t u_i, r_i, len, s_pos_cur, s_pre_v, s_pre_w, p_v, p_idx, p_uId, cBeg, cEnd; - ///buf_t* x; - uint32_t *a, an; - ma_ug_t *ug; - asg_t *read_sg; - trans_chain* t_ch; -} u_trans_hit_idx; + void reset_u_trans_hit_idx(u_trans_hit_idx *t, uint32_t* i_x_a, uint32_t i_x_n, ma_ug_t *i_ug, asg_t *i_read_sg, trans_chain* i_t_ch, uint32_t i_cBeg, uint32_t i_cEnd) @@ -17889,37 +17882,6 @@ hap_cov_t *cov, uint32_t is_update_chain, uint32_t keep_d) } } /****************************may have bugs********************************/ - - /** - if(((m + cur_m) + (np + cur_np)) < (t->m + t->np)) - { - to_replace = 1; - } - else if(((m + cur_m) + (np + cur_np)) == (t->m + t->np)) - { - if(c + cur_c > t->c) - { - to_replace = 1; - } - else if(c + cur_c == t->c) - { - if(nc + cur_nc > t->nc) - { - to_replace = 1; - } - else if(nc + cur_nc == t->nc) - { - if(d + l > t->d) - { - to_replace = 1; - } - } - - } - } - **/ - - if(to_replace) { t->p = v; @@ -29942,6 +29904,178 @@ void flat_bubbles(asg_t *sg, uint8_t* r_het) ma_ug_destroy(ug); free(bs_flag); } +uint64_t get_primary_path_len(asg_t *sg, ma_ug_t *ug, uint32_t v0, buf_t *b) +{ + uint32_t v, u, k; + uint64_t len; + ma_utg_t *p = NULL; + + v = b->S.a[0]; + len = 0; + while (1) + { + u = b->a[v].p; // u->v + if(u == v0) break; + + p = &(ug->u.a[u>>1]); + + for (k = 0; k < p->n; k++) + { + len += sg->seq[p->a[k]>>33].len; + } + + v = u; + } + + return len; +} + +void flat_bubbles_advance(asg_t *sg, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint64_t het_thres) +{ + fprintf(stderr, "het_thres-%lu\n", het_thres); + ma_ug_t *ug = NULL; + ug = ma_ug_gen(sg); + ma_utg_t *u = NULL; + ma_hit_t *h; + uint32_t n_vtx = ug->g->n_seq<<1, v, i, k, m, rId, tn, is_Unitig, n_pop = 0; + uint64_t C_bases, R_bases; + buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); + uint8_t* bs_flag = (uint8_t*)calloc(n_vtx, 1); + uint8_t* r_flag = (uint8_t*)calloc(sg->n_seq, sizeof(uint8_t)); + uint64_t tLen = get_bub_pop_max_dist_advance(ug->g, &b); + + for (v = 0; v < ug->g->n_seq; ++v) + { + if(ug->g->seq[v].del) continue; + ug->g->seq[v].c = PRIMARY_LABLE; + EvaluateLen(ug->u, v) = ug->u.a[v].n; + } + + n_pop = 1; ///round = 0; + while(n_pop > 0) + { + for (v = n_pop = 0; v < n_vtx; ++v) + { + if(ug->g->seq[v>>1].del) continue; + if(asg_arc_n(ug->g, v) < 2) continue; + if(get_real_length(ug->g, v, NULL) < 2) continue; + if(bs_flag[v] == 1) continue; + if(bs_flag[v] == 0) bs_flag[v] = 1; + + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0)) + { + bs_flag[v] = 2; bs_flag[b.S.a[0]^1] = 2; + //note b.b include end, does not include beg + for (i = 0; i < b.b.n; i++) + { + u = &(ug->u.a[b.b.a[i]>>1]); + for (k = 0; k < u->n; k++) r_flag[u->a[k]>>33] = 1; + } + + u = &(ug->u.a[v>>1]); + for (k = 0; k < u->n; k++) r_flag[u->a[k]>>33] = 1; + + u = &(ug->u.a[b.S.a[0]>>1]); + for (k = 0; k < u->n; k++) r_flag[u->a[k]>>33] = 1; + + + + + + + C_bases = 0; + for (i = 0; i < b.b.n; i++) + { + if((b.b.a[i]>>1) == (v>>1) || (b.b.a[i]>>1) == (b.S.a[0]>>1)) continue; + u = &(ug->u.a[b.b.a[i]>>1]); + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + // R_bases += sg->seq[rId].len; + for (m = 0; m < (uint64_t)(sources[rId].length); m++) + { + h = &(sources[rId].buffer[m]); + ///if(h->el != 1) continue; + tn = Get_tn((*h)); + if(sg->seq[tn].del == 1) + { + ///get the id of read that contains it + get_R_to_U(ruIndex, tn, &tn, &is_Unitig); + if(tn == (uint32_t)-1 || is_Unitig == 1 || sg->seq[tn].del == 1) continue; + } + if(sg->seq[tn].del == 1) continue; + if(r_flag[tn] == 0) continue; + C_bases += (Get_qe((*h)) - Get_qs((*h))); + } + } + } + + R_bases = get_primary_path_len(sg, ug, v, &b); + + if((C_bases/R_bases) <= het_thres) + { + fprintf(stderr, "s-utg%.6ul\te-utg%.6ul\tC_bases:%lu\tR_bases:%lu\n", (v>>1)+1, (b.S.a[0]>>1)+1, C_bases, R_bases); + asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0); + n_pop++; + } + + + for (i = 0; i < b.b.n; i++) + { + u = &(ug->u.a[b.b.a[i]>>1]); + for (k = 0; k < u->n; k++) r_flag[u->a[k]>>33] = 0; + } + + u = &(ug->u.a[v>>1]); + for (k = 0; k < u->n; k++) r_flag[u->a[k]>>33] = 0; + + u = &(ug->u.a[b.S.a[0]>>1]); + for (k = 0; k < u->n; k++) r_flag[u->a[k]>>33] = 0; + } + } + ///round++; + } + + /*******************************for debug************************************/ + // kvec_t(uint64_t) occ_sort; kv_init(occ_sort); + // for (v = n_pop = 0; v < ug->g->n_seq; ++v) + // { + // if(ug->g->seq[v].del) continue; + // if(ug->g->seq[v].c != ALTER_LABLE) continue; + // u = &(ug->u.a[v]); + + // kv_push(uint64_t, occ_sort, (uint64_t)((uint32_t)(-1) - (uint32_t)(u->n)) << 32 | (v)); + // } + // radix_sort_arch64(occ_sort.a, occ_sort.a + occ_sort.n); + // for (i = 0; i < occ_sort.n; ++i) + // { + // fprintf(stderr, "-utg%.6ul, n=%u\n", ((uint32_t)occ_sort.a[i])+1, + // (uint32_t)(-1) - (uint32_t)(occ_sort.a[i]>>32)); + // } + // kv_destroy(occ_sort); + /*******************************for debug************************************/ + + clean_sg_by_utg(sg, ug); + + free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); + ma_ug_destroy(ug); free(bs_flag); free(r_flag); +} + + +void flat_soma_v(asg_t *sg, ma_hit_t_alloc* sources, R_to_U* ruIndex) +{ + uint64_t dip_thre_max; + if(asm_opt.hom_global_coverage_set) + { + dip_thre_max = asm_opt.hom_global_coverage; + } + else + { + dip_thre_max = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE); + } + flat_bubbles_advance(sg, sources, ruIndex, (((double)(dip_thre_max)*1.15)/asm_opt.polyploidy)); +} + char *get_outfile_name(char* output_file_name) { char *buf = NULL; @@ -30162,7 +30296,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g) ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz, &b_mask_t); output_unitig_graph(sg, coverage_cut, o_file, sources, ruIndex, max_hang_length, mini_overlap_length); - flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL; + // flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL; + flat_soma_v(sg, sources, ruIndex); free(ruIndex->is_het); ruIndex->is_het = NULL; output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources, asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length); diff --git a/Overlaps.h b/Overlaps.h index 4b870b9..4a8a961 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1228,6 +1228,21 @@ void ma_ug_print(const ma_ug_t *ug, asg_t* read_g, const ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, const char* prefix, FILE *fp); void ma_ug_print_simple(const ma_ug_t *ug, asg_t* read_g, const ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, const char* prefix, FILE *fp); +trans_chain* init_trans_chain(ma_ug_t *ug, uint64_t r_num); +void destory_trans_chain(trans_chain **x); + +typedef struct {///[cBeg, cEnd) + uint32_t u_i, r_i, len, s_pos_cur, s_pre_v, s_pre_w, p_v, p_idx, p_uId, cBeg, cEnd; + ///buf_t* x; + uint32_t *a, an; + ma_ug_t *ug; + asg_t *read_sg; + trans_chain* t_ch; +} u_trans_hit_idx; +void reset_u_trans_hit_idx(u_trans_hit_idx *t, uint32_t* i_x_a, uint32_t i_x_n, ma_ug_t *i_ug, +asg_t *i_read_sg, trans_chain* i_t_ch, uint32_t i_cBeg, uint32_t i_cEnd); +uint32_t get_u_trans_hit(u_trans_hit_idx *t, u_trans_hit_t *hit); + #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/hic.cpp b/hic.cpp index aad0c24..d9763c1 100644 --- a/hic.cpp +++ b/hic.cpp @@ -15536,7 +15536,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o 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, 3); + // horder_t *ho = init_horder_t(&sl.hits, idx->uID_bits, idx->pos_mode, idx->read_g, idx->ug, &bub, &(idx->t_ch->k_trans), opt, 3); ///print_hc_links(&link, 0, &hap); // print_kv_u_trans(&k_trans, &link, s->s); @@ -15549,7 +15549,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o ///print_debug_bubble_graph(&bub, idx->ug, asm_opt.output_file_name); // print_bubble_chain(&bub); // destory_contig_partition(&hap); - destory_horder_t(&ho); + // destory_horder_t(&ho); kv_destroy(back_hc_edge.a); kv_destroy(sl.hits.a); kv_destroy(sl.hits.idx); diff --git a/horder.cpp b/horder.cpp index 76f1486..5de9413 100644 --- a/horder.cpp +++ b/horder.cpp @@ -26,6 +26,7 @@ KRADIX_SORT_INIT(ho64, uint64_t, generic_key, 8) #define osg_arc_key(a) ((a).u) KRADIX_SORT_INIT(osg, osg_arc_t, osg_arc_key, member_size(osg_arc_t, u)) + #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)) @@ -61,7 +62,7 @@ typedef struct { } u_hits_t; typedef struct { - uint64_t e; + uint64_t e, d; double w; } hw_aux_t; @@ -70,6 +71,24 @@ typedef struct { size_t n, m; } h_w_t; +#define hw_e_key(x) ((x).e) +KRADIX_SORT_INIT(hw_e, hw_aux_t, hw_e_key, member_size(hw_aux_t, e)) + +#define hw_d_key(x) ((x).d) +KRADIX_SORT_INIT(hw_d, hw_aux_t, hw_d_key, member_size(hw_aux_t, d)) + +#define hw_ew_key(x) ((uint32_t)((x).e)) +KRADIX_SORT_INIT(hw_ew, hw_aux_t, hw_ew_key, member_size(hw_aux_t, e)) + +#define hw_dw_key(x) ((uint32_t)((x).d)) +KRADIX_SORT_INIT(hw_dw, hw_aux_t, hw_dw_key, member_size(hw_aux_t, d)) + +typedef struct { + kvec_t(uint64_t) pos; + uint64_t *a; + size_t n, m; +} dens_idx_t; + typedef struct { uint64_t s, e, dp; } h_cov_t; @@ -109,6 +128,18 @@ KRADIX_SORT_INIT(h_cov_e, h_cov_t, h_cov_e_key, member_size(h_cov_t, e)) #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)) +typedef struct { + kv_u_trans_t *ref; + trans_chain* idx; +} trans_col_t; + +typedef struct { + uint32_t Spre, Epre, Scur, Ecur, uCur, uPre;///[qSp, qEp) && [qSn, qEn] +} u_hit_t; + +#define u_hit_t_key(x) ((x).uPre) +KRADIX_SORT_INIT(u_hit, u_hit_t, u_hit_t_key, member_size(u_hit_t, uPre)) + void print_N50(ma_ug_t* ug) { kvec_t(uint64_t) b; kv_init(b); @@ -175,6 +206,21 @@ void print_N50_layout(ma_ug_t* ug, sc_lay_t* sl) kv_destroy(b); } +trans_col_t *init_trans_col(ma_ug_t *ug, uint64_t r_num, kv_u_trans_t *ref) +{ + trans_col_t *p = NULL; + CALLOC(p, 1); + p->ref = ref; + p->idx = init_trans_chain(ug, r_num); + return p; +} + +void destory_trans_col(trans_col_t **p) +{ + destory_trans_chain(&((*p)->idx)); + free(*p); +} + 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) { uint64_t k, l, i, r_i, offset, rid, rev, rBeg, rEnd, ubits, p_mode, upos, rpos, update; @@ -1470,12 +1516,264 @@ double get_max_weight(uint32_t u, uint32_t v, osg_t *g) return max; } -void update_scg(horder_t *h) +dens_idx_t *build_interval_idx(kvec_pe_hit *hits, ma_ug_t *ug) { - uint64_t i, k, l, p0s, p0e, p1s, p1e, span_s, span_e, suid, euid, v, w, slen, elen, *ep = NULL, dis; + uint64_t i, k, l, p0s, p0e, p1s, p1e, suid, euid, slen, elen, span_s, span_e; + dens_idx_t *idx = NULL; + CALLOC(idx, 1); + + for (i = 0; i < hits->a.n; i++) + { + if(!hits->a.a[i].id) continue; + suid = get_hit_suid(*hits, i); + euid = get_hit_euid(*hits, i); + slen = ug->u.a[suid].len; + elen = ug->u.a[euid].len; + + p0s = get_hit_spos(*hits, i); + p0e = get_hit_spos_e(*hits, i); + span_s = MIN(p0s, p0e); + span_s = MIN(span_s, slen-1); + span_e = MAX(p0s, p0e); + span_e = MIN(span_e, slen-1); + span_s = ((span_s+span_e)>>1); + kv_push(uint64_t, idx->pos, (suid<<32)|span_s); + + + p1s = get_hit_epos(*hits, i); + p1e = get_hit_epos_e(*hits, i); + span_s = MIN(p1s, p1e); + span_s = MIN(span_s, elen-1); + span_e = MAX(p1s, p1e); + span_e = MIN(span_e, elen-1); + span_s = ((span_s+span_e)>>1); + kv_push(uint64_t, idx->pos, (euid<<32)|span_s); + } + radix_sort_ho64(idx->pos.a, idx->pos.a + idx->pos.n); + idx->n = idx->m = ug->u.n; + CALLOC(idx->a, idx->n); + for (k = 1, l = 0; k <= idx->pos.n; ++k) + { + if (k == idx->pos.n || ((idx->pos.a[k]>>32) != (idx->pos.a[l]>>32))) + { + idx->a[idx->pos.a[l]>>32] = (uint64_t)l << 32 | (k - l); + l = k; + } + } + return idx; +} + +uint64_t get_vw_hits_num(h_w_t *e, uint64_t x, uint64_t y) +{ + uint64_t i, occ; + hw_aux_t *ep = NULL; + for (i = occ = 0; i < e->n; i++) + { + ep = &(e->a[i]); + if(((ep->e>>33) == x && (((uint32_t)(ep->e))>>1) == y) || + ((ep->e>>33) == y && (((uint32_t)(ep->e))>>1) == x)) + { + occ++; + } + } + return occ; +} + +void update_h_w(h_w_t *e, dens_idx_t *idx, double *max_div) +{ + uint64_t k, l, ii, pi, pos, ori, uid, *id = NULL, idn; + if(max_div) (*max_div) = 0; + radix_sort_hw_e(e->a, e->a+e->n); + for (k = 1, l = 0; k <= e->n; ++k) + { + if (k == e->n || (e->a[k].e>>32) != (e->a[l].e>>32)) ///same uid + { + if(k - l > 1) radix_sort_hw_d(e->a+l, e->a+k); + ori = (e->a[l].e>>32)&1; + uid = e->a[l].e>>33; + id = idx->pos.a + (idx->a[uid]>>32); + idn = (uint32_t)(idx->a[uid]); + ii = 0; + for (pi = l; pi < k; pi++) + { + pos = e->a[pi].d>>32;/// + while (ii < idn) + { + if(((uint32_t)id[ii]) == pos) + { + if(ori) + { + break; + } + else + { + while (ii < idn && (((uint32_t)id[ii]) == pos)) + { + ii++; + } + ii--; + break; + } + } + ii++; + } + if(ii >= idn) fprintf(stderr, "ERROR-1\n"); + e->a[pi].w += (ori? idn-ii: ii+1); + if(max_div) (*max_div) = MAX((*max_div), e->a[pi].w); + } + l = k; + } + } + + radix_sort_hw_ew(e->a, e->a+e->n); + for (k = 1, l = 0; k <= e->n; ++k) + { + if (k == e->n || ((uint32_t)(e->a[k].e)) != ((uint32_t)(e->a[l].e))) ///same uid + { + if(k - l > 1) radix_sort_hw_dw(e->a+l, e->a+k); + ori = e->a[l].e&1; + uid = (((uint32_t)e->a[l].e)>>1); + id = idx->pos.a + (idx->a[uid]>>32); + idn = (uint32_t)(idx->a[uid]); + ii = 0; + for (pi = l; pi < k; pi++) + { + pos = (uint32_t)(e->a[pi].d);/// + while (ii < idn) + { + if(((uint32_t)id[ii]) == pos) + { + if(ori) + { + break; + } + else + { + while (ii < idn && (((uint32_t)id[ii]) == pos)) + { + ii++; + } + ii--; + break; + } + } + ii++; + } + if(ii >= idn) fprintf(stderr, "ERROR-2\n"); + e->a[pi].w += (ori? idn-ii: ii+1); + if(max_div) (*max_div) = MAX((*max_div), e->a[pi].w); + } + l = k; + } + } + // radix_sort_hw_e(e->a, e->a+e->n); + if(max_div) (*max_div) *= 2;///different with slsa2 +} + +void print_specfic_hic_hits(kvec_pe_hit *hits, uint64_t v, uint64_t w, ma_ug_t *ug) +{ + uint64_t i, suid, euid, slen, elen, p0s, p0e, p1s, p1e, span_s, span_e, sd, ed; + for (i = 0; i < hits->a.n; i++) + { + if(!hits->a.a[i].id) continue; + suid = get_hit_suid(*hits, i); + euid = get_hit_euid(*hits, i); + if((suid == v && euid == w) || (suid == w && euid == v)) + { + slen = ug->u.a[suid].len; + elen = ug->u.a[euid].len; + + p0s = get_hit_spos(*hits, i); + p0e = get_hit_spos_e(*hits, i); + span_s = MIN(p0s, p0e); + span_s = MIN(span_s, slen-1); + span_e = MAX(p0s, p0e); + span_e = MIN(span_e, slen-1); + span_s = ((span_s+span_e)>>1); + sd = span_s; + + + p1s = get_hit_epos(*hits, i); + p1e = get_hit_epos_e(*hits, i); + span_s = MIN(p1s, p1e); + span_s = MIN(span_s, elen-1); + span_e = MAX(p1s, p1e); + span_e = MIN(span_e, elen-1); + span_s = ((span_s+span_e)>>1); + ed = span_s; + + fprintf(stderr, "u-stg%.6lul\t%lu\tv-stg%.6lul\t%lu\n", suid+1, sd, euid + 1, ed); + } + } +} + +kv_u_trans_t *get_update_trans_idx(horder_t *h, trans_col_t *t_idx) +{ + uint64_t i, k, l; + uint32_t uid/**, q_n, t_n**/; + ma_ug_t *ug = h->ug; + asg_t *rg = h->r_g; + kv_u_trans_t *idx = NULL; CALLOC(idx, 1); + u_trans_hit_idx iter; + u_trans_hit_t hit; + kvec_t(u_hit_t) b; kv_init(b); + uint64_t *b_idx = NULL; CALLOC(b_idx, t_idx->ref->idx.n); + /**u_hit_t *q, *t;**/ + u_hit_t *p = NULL; + for (i = 0; i < ug->u.n; i++) + { + uid = i<<1; + reset_u_trans_hit_idx(&iter, &uid, 1, ug, rg, t_idx->idx, 0, ug->u.a[i].len); + while(get_u_trans_hit(&iter, &hit))///get [qScur, qEcur), [qSpre, qEpre) + { + kv_pushp(u_hit_t, b, &p); + p->Scur = hit.qScur; + p->Ecur = hit.qEcur; + p->Spre = hit.qSpre; + p->Epre = hit.qEpre; + p->uPre = hit.qn; + p->uCur = uid; + } + } + + radix_sort_u_hit(b.a, b.a + b.n); + for (k = 1, l = 0; k <= b.n; ++k) + { + if (k == b.n || (b.a[k].uPre>>1) != (b.a[l].uPre>>1)) + { + b_idx[b.a[l].uPre>>1] = (uint64_t)l << 32 | (k - l); + l = k; + } + } + + /** + for (i = 0; i < t_idx->ref->idx.n; i++) + { + q_n = (uint32_t)b_idx[i]; + q = b.a + (b_idx[i]>>32); + if(!q_n) continue; + } + **/ + + // for (i = 0; i < b.n; i++) + // { + + // } + + + free(b_idx); + kv_destroy(b); + return idx; +} + +void update_scg(horder_t *h, trans_col_t *t_idx) +{ + uint64_t i, k, l, p0s, p0e, p1s, p1e, span_s, span_e, suid, euid, v, w, slen, elen, sd, ed; uint64_t t_hits = 0, a_hits = 0; - double div, max_div; - kvec_t(uint64_t) e; kv_init(e); + double max_div, we; + h_w_t e; kv_init(e); + dens_idx_t *idx = NULL; + hw_aux_t *ep = NULL; ma_ug_t *ug = h->ug; kvec_pe_hit *hits = &(h->u_hits); osg_arc_t *p = NULL; @@ -1488,9 +1786,11 @@ void update_scg(horder_t *h) h->sg.g->seq[i].ez[0] = ug->u.a[i].len>>1; h->sg.g->seq[i].ez[1] = ug->u.a[i].len - (ug->u.a[i].len>>1); } - - for (i = 0, max_div = 1, e.n = 0; i < hits->a.n; i++) + idx = build_interval_idx(hits, ug); + + for (i = 0, e.n = 0; i < hits->a.n; i++) { + if(!hits->a.a[i].id) continue; suid = get_hit_suid(*hits, i); euid = get_hit_euid(*hits, i); if(suid == euid) continue; @@ -1504,6 +1804,7 @@ void update_scg(horder_t *h) span_e = MAX(p0s, p0e); span_e = MIN(span_e, slen-1); span_s = ((span_s+span_e)>>1); + sd = span_s; v = suid << 1; if(span_s > (slen>>1)) v++; @@ -1515,32 +1816,45 @@ void update_scg(horder_t *h) span_e = MAX(p1s, p1e); span_e = MIN(span_e, elen-1); span_s = ((span_s+span_e)>>1); + ed = span_s; w = euid << 1; if(span_s > (elen>>1)) w++; t_hits++; + + kv_pushp(hw_aux_t, e, &ep); + ep->w = 0; + ep->e = (v<<32)|w; + ep->d = (sd<<32)|ed; - if(hits->a.a[i].id) + if(v > w) { - kv_pushp(uint64_t, e, &ep); - (*ep) = (v<<32)|w; - if(v > w) (*ep) = (w<<32)|v; - // (*ep) <<= 1; (*ep) |= ((uint64_t)(!!hits->a.a[i].id)); - div = h->sg.g->seq[v>>1].ez[v&1] + h->sg.g->seq[w>>1].ez[w&1]; - max_div = MAX(max_div, div); - a_hits++; - } + ep->e = (w<<32)|v; + ep->d = (ed<<32)|sd; + } + + // div = h->sg.g->seq[v>>1].ez[v&1] + h->sg.g->seq[w>>1].ez[w&1]; + // max_div = MAX(max_div, div); + a_hits++; } - max_div *= 2;///different with slsa2 - radix_sort_ho64(e.a, e.a+e.n); + // fprintf(stderr, "sa-0-sa: occ-%lu\n", get_vw_hits_num(&e, 12, 690)); + // print_specfic_hic_hits(hits, 12, 690, ug); + + + + update_h_w(&e, idx, &max_div); + + // fprintf(stderr, "sa-1-sa: occ-%lu\n", get_vw_hits_num(&e, 676, 738)); + + radix_sort_hw_e(e.a, e.a+e.n); for (k = 1, l = 0; k <= e.n; ++k) { - if (k == e.n || e.a[k] != e.a[l]) + if (k == e.n || e.a[k].e != e.a[l].e) { for (i = 0; i < h->avoid.n; i++) { - v = e.a[l]; - w = e.a[l]<<32; w |= (e.a[l]>>32); + v = e.a[l].e; + w = e.a[l].e<<32; w |= (e.a[l].e>>32); if(h->avoid.a[i] == v || h->avoid.a[i] == w) { break; @@ -1549,29 +1863,34 @@ void update_scg(horder_t *h) if(i >= h->avoid.n) { - div = h->sg.g->seq[e.a[l]>>33].ez[(e.a[l]>>32)&1] + - h->sg.g->seq[((uint32_t)e.a[l])>>1].ez[e.a[l]&1]; + for (i = l, we = 0; i < k; i++) + { + if(e.a[i].w == 0) fprintf(stderr, "ERROR-3\n"); + we += (max_div/e.a[i].w); + } + p = osg_arc_pushp(h->sg.g); p->u = p->v = p->occ = p->del = p->w = p->nw = 0; - p->u = e.a[l]>>32; p->v = (uint32_t)e.a[l]; - p->occ = k - l; - if(div != 0) p->w = (double)(k - l)*(max_div/div); + p->u = e.a[l].e>>32; p->v = (uint32_t)e.a[l].e; + p->occ = k - l; p->w = we; + // if(div != 0) p->w = (double)(k - l)*(max_div/div); p = osg_arc_pushp(h->sg.g); p->u = p->v = p->occ = p->del = p->w = p->nw = 0; - p->u = (uint32_t)e.a[l]; p->v = e.a[l]>>32; - p->occ = k - l; - if(div != 0) p->w = (double)(k - l)*(max_div/div); + p->u = (uint32_t)e.a[l].e; p->v = e.a[l].e>>32; + p->occ = k - l; p->w = we; + - h->sg.g->seq[e.a[l]>>33].mw[(e.a[l]>>32)&1] - = MAX(h->sg.g->seq[e.a[l]>>33].mw[(e.a[l]>>32)&1], p->w); - h->sg.g->seq[((uint32_t)e.a[l])>>1].mw[e.a[l]&1] - = MAX(h->sg.g->seq[((uint32_t)e.a[l])>>1].mw[e.a[l]&1], p->w); + h->sg.g->seq[e.a[l].e>>33].mw[(e.a[l].e>>32)&1] + = MAX(h->sg.g->seq[e.a[l].e>>33].mw[(e.a[l].e>>32)&1], p->w); + h->sg.g->seq[((uint32_t)e.a[l].e)>>1].mw[e.a[l].e&1] + = MAX(h->sg.g->seq[((uint32_t)e.a[l].e)>>1].mw[e.a[l].e&1], p->w); } l = k; } } + // fprintf(stderr, "sa-2-sa: occ-%lu\n", get_vw_hits_num(&e, 676, 738)); osg_cleanup(h->sg.g); double bestAlt; uint64_t eg_edges = 0; @@ -1593,6 +1912,7 @@ void update_scg(horder_t *h) __func__, h->sg.g->n_seq, h->sg.g->n_arc, eg_edges, t_hits, a_hits); /*******************************for debug************************************/ + /** for (i = 0; i < h->sg.g->n_arc; i++) { p = &(h->sg.g->arc[i]); @@ -1601,6 +1921,8 @@ void update_scg(horder_t *h) (p->v>>1)+1, "+-"[p->v&1], h->sg.g->seq[p->v>>1].ez[p->v&1], p->occ, p->w, p->nw); } fprintf(stderr, "sbsbsbsb\n\n\n\n\n\n"); + **/ + // uint32_t u, nv, f; // osg_arc_t *av = NULL; // for (k = 0; k < h->sg.g->n_arc; k++) @@ -1648,7 +1970,10 @@ void update_scg(horder_t *h) // } /*******************************for debug************************************/ - kv_destroy(e); + kv_destroy(e); + free(idx->pos.a); + free(idx->a); + free(idx); } int cmp_arc_nw(const void * a, const void * b) @@ -2129,7 +2454,7 @@ uint32_t get_sl_occ(sc_lay_t *sl) void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis) { - uint32_t k, m, max_utg, max_sc; + uint32_t k, m/**, max_utg, max_sc**/; lay_t *p = NULL; uint8_t *sgv = NULL; MALLOC(sgv, sl->n); double *w = NULL; MALLOC(w, sl->n); @@ -2494,7 +2819,7 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp) kv_destroy(new_rtg_edges.a); } -void scaffold_hap(horder_t *h, ug_opt_t *opt, uint32_t round, char *output_file_name, uint8_t flag) +void scaffold_hap(horder_t *h, ug_opt_t *opt, trans_col_t *t_idx, uint32_t round, char *output_file_name, uint8_t flag) { uint32_t i; kv_destroy(h->u_hits.a); @@ -2518,7 +2843,7 @@ void scaffold_hap(horder_t *h, ug_opt_t *opt, uint32_t round, char *output_file_ for (i = 0; i < round; i++) { update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); - update_scg(h); + update_scg(h, t_idx); layout_scg(h, 1.001, 19); renew_scaffold(h); } @@ -2545,19 +2870,21 @@ void output_hic_rtg(ma_ug_t *ug, asg_t *rg, ug_opt_t *opt, char* output_file_nam } horder_t *init_horder_t(kvec_pe_hit *i_hits, uint64_t i_hits_uid_bits, uint64_t i_hits_pos_mode, -asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt, uint32_t round) +asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, kv_u_trans_t *ref, ug_opt_t *opt, uint32_t round) { uint32_t i; + trans_col_t *t_idx = NULL; 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); + t_idx = init_trans_col(i_ug, h->r_g->n_seq, ref); // output_hic_rtg(i_ug, h->r_g, opt, asm_opt.output_file_name); // reduce_hamming_error(h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz); /** - scaffold_hap(h, opt, round, asm_opt.output_file_name, FATHER); - scaffold_hap(h, opt, round, asm_opt.output_file_name, MOTHER); + scaffold_hap(h, t_idx, opt, round, asm_opt.output_file_name, FATHER); + scaffold_hap(h, t_idx, opt, round, asm_opt.output_file_name, MOTHER); **/ @@ -2569,7 +2896,7 @@ asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt, uint32_t round) for (i = 0; i < round; i++) { update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); - update_scg(h); + update_scg(h, t_idx); layout_scg(h, 1.001, 19); renew_scaffold(h); } @@ -2580,6 +2907,7 @@ asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt, uint32_t round) opt->stops_threshold, opt->ruIndex, opt->chimeric_rate, opt->drop_ratio, opt->max_hang, opt->min_ovlp); + destory_trans_col(&t_idx); exit(1); return h; } diff --git a/horder.h b/horder.h index 8e10fde..3d98e39 100644 --- a/horder.h +++ b/horder.h @@ -37,6 +37,6 @@ typedef struct { }horder_t; horder_t *init_horder_t(kvec_pe_hit *i_hits, uint64_t i_hits_uid_bits, uint64_t i_hits_pos_mode, -asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt, uint32_t round); +asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, kv_u_trans_t *ref, ug_opt_t *opt, uint32_t round); void destory_horder_t(horder_t **h); #endif