diff --git a/Overlaps.cpp b/Overlaps.cpp index 53d566a..b9878c8 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -16280,6 +16280,18 @@ uint32_t cmp_untig_graph(ma_ug_t *src, ma_ug_t *dest) fprintf(stderr, "src->g->seq[%u].len: %u, dest->g->seq[%u].len: %u\n", v, src->g->seq[v].len, v, dest->g->seq[v].len); } + + if(src->u.a[v].n != dest->u.a[v].n) + { + fprintf(stderr, "src->u.a[%u].n: %u, dest->u.a[%u].n: %u\n", + v, src->u.a[v].n, v, dest->u.a[v].n); + } + + if(src->u.a[v].a[0] != dest->u.a[v].a[0] || + src->u.a[v].a[src->u.a[v].n-1] != dest->u.a[v].a[dest->u.a[v].n - 1]) + { + fprintf(stderr, "unequal beg/end node\n"); + } } @@ -28778,7 +28790,7 @@ void reset_bub(bubble_type* bub, ma_ug_t *ug, trans_chain* back_ug_chain, kvec_a destory_bubbles(bub); memset(bub, 0, sizeof(bubble_type)); - new_rtg_edges->a.n = 0; + if(new_rtg_edges) new_rtg_edges->a.n = 0; ///classify_untigs(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, max_hang, min_ovlp); identify_bubbles(ug, bub, back_ug_chain->ir_het, NULL); // update_bubble_chain(ug, bub, 0, 1); @@ -31398,19 +31410,23 @@ char *get_outfile_name(char* output_file_name) } void gen_ug_opt_t(ug_opt_t *opt, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int64_t max_hang, int64_t min_ovlp, -int64_t gap_fuzz, int64_t min_dp, uint64_t* readLen, ma_sub_t *coverage_cut, R_to_U* ruIndex) +int64_t gap_fuzz, int64_t min_dp, uint64_t* readLen, ma_sub_t *coverage_cut, R_to_U* ruIndex, long long tipsLen, +float tip_drop_ratio, long long stops_threshold, float chimeric_rate, float drop_ratio, bub_label_t* b_mask_t) { memset(opt, 0, sizeof((*opt))); opt->sources = sources; opt->reverse_sources = reverse_sources; opt->max_hang = max_hang; opt->min_ovlp = min_ovlp; opt->gap_fuzz = gap_fuzz; opt->min_dp = min_dp; opt->readLen = readLen; - opt->coverage_cut = coverage_cut; opt->ruIndex = ruIndex; + opt->coverage_cut = coverage_cut; opt->ruIndex = ruIndex; opt->tipsLen = tipsLen; + opt->tip_drop_ratio = tip_drop_ratio; opt->stops_threshold = stops_threshold; + opt->chimeric_rate = chimeric_rate; opt->drop_ratio = drop_ratio; opt->b_mask_t = b_mask_t; } void create_ul_info(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int64_t max_hang, int64_t min_ovlp, int64_t gap_fuzz, -int64_t min_dp, uint64_t* readLen, ma_sub_t *coverage_cut, R_to_U* ruIndex) +int64_t min_dp, uint64_t* readLen, ma_sub_t *coverage_cut, R_to_U* ruIndex, long long tipsLen, float tip_drop_ratio, long long stops_threshold, float chimeric_rate, float drop_ratio, bub_label_t* b_mask_t) { ug_opt_t opt; - gen_ug_opt_t(&opt, sources, reverse_sources, max_hang, min_ovlp, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex); + gen_ug_opt_t(&opt, sources, reverse_sources, max_hang, min_ovlp, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex, + tipsLen, tip_drop_ratio, stops_threshold, chimeric_rate, drop_ratio, b_mask_t); ul_load(&opt); } @@ -31480,13 +31496,16 @@ ma_sub_t **coverage_cut_ptr, int debug_g) { memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t)); } - + if(asm_opt.ar) init_all_ul_t(&UL_INF, &R_INF); // if (asm_opt.flag & HA_F_VERBOSE_GFA) { // write_debug_graph(NULL, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF); // debug_gfa:; // } ///should recover edges from sources by using UL alignments - if(asm_opt.ar) create_ul_info(sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex); + if(asm_opt.ar) { + create_ul_info(sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex, + (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, &b_mask_t); + } clean_weak_ma_hit_t(sources, reverse_sources, n_read, asm_opt.ar?UL_COV_THRES:(uint32_t)-1); sg = gen_init_sg(min_dp, n_read, mini_overlap_length, max_hang_length, gap_fuzz, sources, readLen, ruIndex, @@ -31519,14 +31538,16 @@ ma_sub_t **coverage_cut_ptr, int debug_g) free(unlean_name); } - gen_ug_opt_t(&uopt, sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex); + gen_ug_opt_t(&uopt, sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex, + (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, &b_mask_t); ul_clean_gfa(&uopt, sg, sources, reverse_sources, ruIndex, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.6, asm_opt.max_short_tip, &b_mask_t, !!asm_opt.ar, ha_opt_triobin(&asm_opt), UL_COV_THRES, o_file); if (asm_opt.flag & HA_F_VERBOSE_GFA) { write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF); debug_gfa:; - gen_ug_opt_t(&uopt, sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex); + gen_ug_opt_t(&uopt, sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex, + (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, &b_mask_t); } if(asm_opt.ar) ul_realignment_gfa(&uopt, sg); print_debug_gfa(sg, NULL, coverage_cut, "UL.debug", sources, ruIndex, max_hang_length, mini_overlap_length); diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 877ab69..6c2e746 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -9,6 +9,8 @@ #include "Correct.h" #include "inter.h" #include "Overlaps.h" +#include "hic.h" +#include "Purge_Dups.h" #define generic_key(x) (x) KRADIX_SORT_INIT(srt64, uint64_t, generic_key, 8) @@ -17,12 +19,174 @@ KRADIX_SORT_INIT(srt64, uint64_t, generic_key, 8) #define ASG_ET_TIP 1 #define ASG_ET_MULTI_OUT 2 #define ASG_ET_MULTI_NEI 3 +#define UL_TRAV_HERATE 0.2 +#define UL_TRAV_FT_RATE 0.8 typedef struct { asg_t *g; ma_hit_t_alloc *src; } sset_aux; + +typedef struct { + kvec_t(uint64_t) ref; + kvec_t(uint64_t) pat; + kvec_t(uint64_t) pat_cor; + kvec_t(uint8_t) g_flt; + kvec_t(uint8_t) m_dir; + kvec_t(int64_t) m_score; + uint64_t n, m; +} path_dp_t; + +typedef struct { + uint32_t p; // the optimal parent vertex + uint32_t d; // the shortest distance from the initial vertex + uint64_t c; // max count of positive reads + uint64_t m; // max count of negative reads + // uint32_t np; // max count of non-positive reads + uint32_t nc; // max count of reads, no matter positive or negative + uint32_t r:31, s:1; // r: the number of remaining incoming arc; s: state + //s: state, s=0, this edge has not been visited, otherwise, s=1 +} uinfo_t; + +typedef struct { + uint64_t c; // max count of positive reads + uint32_t i; +} uinfo_srt_t; + +#define uinfo_srt_t_c_key(p) ((p).c) +KRADIX_SORT_INIT(uinfo_srt_t_c, uinfo_srt_t, uinfo_srt_t_c_key, member_size(uinfo_srt_t, c)) + +typedef struct { + ///all information for each node + kvec_t(uinfo_t) a; + // kvec_t(uint32_t) u; + kvec_t(uint32_t) S; // set of vertices without parents, nodes with all incoming edges visited + kvec_t(uint32_t) T; // set of tips + kvec_t(uint32_t) b; // visited vertices + kvec_t(uint32_t) e; // visited edges/arcs + // kvec_t(uinfo_srt_t) srt; + kvec_t(uint8_t) us; + path_dp_t dp; +} ubuf_t; + +typedef struct{ + uint32_t bid, beg, occ; + uint32_t n_path, path_idx, path_occ; +}ul_sub_path_t; + + +typedef struct{ + ma_ug_t *buf_ug; + kvec_t(uint64_t) buf; +}ul_path_t; + +typedef struct{ + kvec_t(uint64_t) idx; + kvec_t(uint64_t) srt; + uint64_t ul_n; +}ul_path_srt_t; + +typedef struct{ + size_t n, m; + uint64_t *a; + uint32_t cn:31, is_cir:1; +}ul_str_t; + +typedef struct{ + kvec_t(uint64_t) idx; + kvec_t(uint64_t) occ; + kvec_t(ul_str_t) str; +}ul_str_idx_t; + +typedef struct{ + uint32_t v, pi, ai; + uint32_t k:31, is_gc:1; + int32_t dis; +} integer_seq_t; + +typedef struct{ + size_t n, m; + integer_seq_t *a; +} kv_integer_seq_t; + +typedef struct{ + uint32_t tk, vq; + uint64_t tn_rev_qk; +} integer_aln_t; + +#define integer_aln_t_vqk_key(x) ((x).tn_rev_qk) +KRADIX_SORT_INIT(integer_aln_t_srt, integer_aln_t, integer_aln_t_vqk_key, member_size(integer_aln_t, tn_rev_qk)) + + +typedef struct { + size_t n, m; + integer_aln_t *a; +} integer_aln_vec_t; + +typedef struct{ + uint32_t s, e, v; + uint64_t sc; + uint32_t q_sidx, q_eidx; + uint32_t t_sidx, t_eidx; +} ul_chain_t; + + + +typedef struct{ + uint64_t qidx_occ; + uint32_t tidx_occ; + uint32_t chain_id:31, is_rev:1; +} ul_snp_t; + +#define ul_snp_t_srt_key(x) ((x).qidx_occ) +KRADIX_SORT_INIT(ul_snp_t_srt, ul_snp_t, ul_snp_t_srt_key, member_size(ul_snp_t, qidx_occ)) + + +typedef struct { + kv_integer_seq_t q; + kv_integer_seq_t t; + integer_aln_vec_t b; + kvec_t(int64_t) f; + kvec_t(int64_t) p; + kvec_t(uint64_t) o; + kvec_t(uint64_t) u; + // kvec_t(uint64_t) srt; + // kvec_t(uint64_t) v; + // kvec_t(uint64_t) u; + // kvec_t(uint64_t) d; + kvec_t(ul_chain_t) sc; + kvec_t(ul_snp_t) snp; +}integer_t; + +typedef struct { + // ul_resolve_t *u; + integer_t *buf; + uint64_t n_thread; +}integer_ml_t; + +typedef struct{ + ma_ug_t *init_ug; + ma_ug_t *l0_ug; + ma_ug_t *l1_ug; + asg_t *sg; + bubble_type *bub; + all_ul_t *idx; + ul_path_t path; + uint8_t *r_het; + ubuf_t buf; + ul_str_idx_t pstr; + integer_ml_t str_b; + // ul_path_srt_t psrt; +}ul_resolve_t; + +void init_integer_ml_t(integer_ml_t *x, ul_resolve_t *u, uint64_t n_thread) +{ + memset(x, 0, sizeof((*x))); + x->n_thread = n_thread; ///x->u = u; + CALLOC(x->buf, n_thread); +} + int32_t if_sup_chimeric(ma_hit_t_alloc* src, uint64_t rLen, asg64_v *b, int if_exact); void print_edge(asg_arc_t *t, const char *cmd) @@ -1323,12 +1487,84 @@ void print_vw_edge(asg_t *sg, uint32_t v, uint32_t w, const char *cmd) if(i >= nv) fprintf(stderr, "[%s]\tno edges\n", cmd); } -void fill_containment_by_ul(asg_t *g, ma_hit_t_alloc *src) +int32_t gen_spec_edge(asg_t *rg, ug_opt_t *uopt, uint32_t v, uint32_t w, asg_arc_t *t) { - + uint32_t k, qn, tn; int32_t r; ma_hit_t_alloc *s = &(uopt->sources[v>>1]); asg_arc_t p; + for (k = 0; k < s->length; k++) { + qn = Get_qn(s->buffer[k]); tn = Get_tn(s->buffer[k]); + if(tn != (w>>1)) continue; + r = ma_hit2arc(&(s->buffer[k]), rg->seq[qn].len, rg->seq[tn].len, uopt->max_hang, asm_opt.max_hang_rate, uopt->min_ovlp, &p); + if(r < 0) continue; + if((p.ul>>32) != v || p.v != w) continue; + *t = p; t->ou = 0; + return 1; + } + return -1; } +void filter_sg_by_ug(asg_t *rg, ma_ug_t *ug, ug_opt_t *uopt) +{ + uint32_t i, m, v, w, nv, n_vx = ug->g->n_seq<<1, vx, wx; int32_t r; + asg_arc_t *av = NULL; ma_utg_t *u = NULL; asg_arc_t *p, t; + n_vx = rg->n_seq; rg->n_arc = 0; + for (v = 0; v < n_vx; v++) rg->seq[v].del = (!!1); + for (i = 0; i < ug->g->n_seq; ++i) { + ug->g->seq[i].c = PRIMARY_LABLE;; + if(ug->g->seq[i].del) continue; + u = &(ug->u.a[i]); + for (m = 0; m < u->n; m++) rg->seq[u->a[m]>>33].del = (!!0); + for (m = 0; (m + 1) < u->n; m++){ + v = u->a[m]>>32; w = u->a[m+1]>>32; + r = gen_spec_edge(rg, uopt, v, w, &t); + assert(r >= 0); p = asg_arc_pushp(rg); *p = t; + + r = gen_spec_edge(rg, uopt, w^1, v^1, &t); + assert(r >= 0); p = asg_arc_pushp(rg); *p = t; + } + + + v = i<<1; nv = asg_arc_n(ug->g, v); av = asg_arc_a(ug->g, v); + for (m = 0; m < nv; m++) { + if(av[m].del) continue; + w = av[m].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)); + + r = gen_spec_edge(rg, uopt, vx, wx, &t); + assert(r >= 0); p = asg_arc_pushp(rg); *p = t; + + // r = gen_spec_edge(rg, uopt, wx^1, vx^1, &t); + // assert(r >= 0); p = asg_arc_pushp(rg); *p = t; + } + + v = (i<<1)+1; nv = asg_arc_n(ug->g, v); av = asg_arc_a(ug->g, v); + for (m = 0; m < nv; m++) { + if(av[m].del) continue; + w = av[m].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)); + + r = gen_spec_edge(rg, uopt, vx, wx, &t); + assert(r >= 0); p = asg_arc_pushp(rg); *p = t; + + // r = gen_spec_edge(rg, uopt, wx^1, vx^1, &t); + // assert(r >= 0); p = asg_arc_pushp(rg); *p = t; + } + } + + free(rg->idx); + rg->idx = 0; + rg->is_srt = 0; + asg_cleanup(rg); + + + /*******************************for debug************************************/ + // ma_ug_t *dbg = ma_ug_gen(rg); + // for (i = 0; i < dbg->g->n_seq; ++i) dbg->g->seq[i].c = PRIMARY_LABLE; + // cmp_untig_graph(dbg, ug); + /*******************************for debug************************************/ +} void ul_clean_gfa(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, int64_t clean_round, double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, double ou_drop_rate, int64_t max_tip, bub_label_t *b_mask_t, int32_t is_ou, int32_t is_trio, uint32_t ou_thres, char *o_file) @@ -1425,14 +1661,1725 @@ double ou_drop_rate, int64_t max_tip, bub_label_t *b_mask_t, int32_t is_ou, int3 free(bu.a); } -void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg) + +bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het) +{ + kvec_asg_arc_t_warp new_rtg_edges; + kv_init(new_rtg_edges.a); + hap_cov_t *cov = NULL; + bubble_type *bub = NULL; + + asg_t *copy_sg = copy_read_graph(sg); + ma_ug_t *copy_ug = copy_untig_graph(ug); + + adjust_utg_by_primary(©_ug, copy_sg, TRIO_THRES, uopt->sources, uopt->reverse_sources, uopt->coverage_cut, + uopt->tipsLen, uopt->tip_drop_ratio, uopt->stops_threshold, uopt->ruIndex, uopt->chimeric_rate, uopt->drop_ratio, + uopt->max_hang, uopt->min_ovlp, &new_rtg_edges, &cov, uopt->b_mask_t, 0, 0); + ma_ug_destroy(copy_ug); copy_ug = NULL; + asg_destroy(copy_sg); copy_sg = NULL; + + CALLOC(bub, 1); (*ir_het) = cov->t_ch->ir_het; + cov->t_ch->ir_het = NULL; cov->is_r_het = NULL; + identify_bubbles(ug, bub, (*ir_het), NULL); + + kv_destroy(new_rtg_edges.a); destory_hap_cov_t(&cov); + return bub; +} + +void clear_path_dp_t(path_dp_t *x, asg_t *g) +{ + uint32_t n_vx = g->n_seq<<1; + x->ref.n = x->pat.n = x->pat_cor.n = 0; + kv_resize(uint8_t, x->g_flt, n_vx); x->g_flt.n = n_vx; + memset(x->g_flt.a, 0, sizeof(*(x->g_flt.a))*x->g_flt.n); +} + +void clear_ubuf_t(ubuf_t *x, asg_t *g, all_ul_t *ul_idx) +{ + uint32_t n_vx = g->n_seq<<1; + x->a.n = x->S.n = x->T.n = x->b.n = x->e.n = 0; + kv_resize(uinfo_t, x->a, n_vx); x->a.n = n_vx; memset(x->a.a, 0, sizeof(*(x->a.a))*x->a.n); + clear_path_dp_t(&(x->dp), g); +} + +uint64_t get_ul_read_weight(all_ul_t *ul, uint32_t *prg, uinfo_t *g_idx, ul_vec_t *p, uint32_t ii, uint32_t v, uint32_t w) +{ + assert((!p->bb.a[ii].base)&&(p->bb.a[ii].hid == (v>>1))&&(p->bb.a[ii].el)&&(p->bb.a[ii].pchain)); + if(p->bb.a[ii].aidx == (uint32_t)-1) return 0; ///not connected + uc_block_t *li = NULL, *lk = NULL; uint32_t li_v, lk_v; + li = &(p->bb.a[ii]); li_v = (((uint32_t)(li->hid))<<1)|((uint32_t)(li->rev)); //li_v^=1; + if(li_v != v) return 0; + + for (lk = &(p->bb.a[li->aidx]), w^=1; lk; ) { + lk_v = (((uint32_t)(lk->hid))<<1)|((uint32_t)(lk->rev)); //lk_v^=1; + while (lk_v != (w^1)) { + if(g_idx[w].p==(uint32_t)-1) break; + w = g_idx[w].p; + } + + if(lk_v == (w^1)) { + + } + + lk = ((lk->aidx==(uint32_t)-1)?NULL:&(p->bb.a[lk->aidx])); + + } + return 1; +} + +#define arc_first(g, v) ((g)->arc[(g)->idx[(v)]>>32]) +#define arc_cnt(g, v) ((uint32_t)(g)->idx[(v)]) + +uint64_t ulg_len_check(asg_t *g, uint32_t v, uc_block_t *p) +{ + // if(!(((v&1) && (p->ts==0)) || (((v&1)==0) && (p->te==g->seq[v>>1].len)))) { + // fprintf(stderr, "\n[M::%s::] v>>1:%u, v&1:%u, ts:%u, te:%u, tlen:%u, qs:%u, qe:%u, pidx:%u, aidx:%u, flag:%u\n", + // __func__, v>>1, v&1, p->ts, p->te, g->seq[v>>1].len, p->qs, p->qe, p->pidx, p->aidx, flag); + // } + // assert(((v&1) && (p->ts==0)) || (((v&1)==0) && (p->te==g->seq[v>>1].len))); + if(!(((v&1) && (p->ts==0)) || (((v&1)==0) && (p->te==g->seq[v>>1].len)))) return 0; + // if((v&1) && (p->ts!=0)) return 0; + // if(((v&1)==0) && (p->te!=g->seq[v>>1].len)) return 0; + if(p->ts==0 && p->te==g->seq[v>>1].len) return 1; + int32_t ol = arc_first(g, v).ol;///max len + int32_t tl = p->te - p->ts; + tl -= ol; + if(tl > 20000 || tl > (g->seq[v>>1].len*0.2)) return 1; + return 0; +} + +void update_l_coord(asg_t *g, uint64_t o_s, uint32_t v, uc_block_t *p, uint64_t *s, uint64_t *e) +{ + // if(!(((v&1) && (p->ts==0)) || (((v&1)==0) && (p->te==g->seq[v>>1].len)))) { + // fprintf(stderr, "\n[M::%s::] v>>1:%u, v&1:%u, ts:%u, te:%u, tlen:%u\n", + // __func__, v>>1, v&1, p->ts, p->te, g->seq[v>>1].len); + // } + // assert(((v&1) && (p->ts==0)) || (((v&1)==0) && (p->te==g->seq[v>>1].len))); + // assert((v&1) && (p->ts==0)); + // assert(((v&1)==0) && (p->te==g->seq[v>>1].len)); + if(((v&1)==0)) { + (*s) = o_s + p->ts; (*e) = o_s + p->te; + } else { + (*s) = o_s + g->seq[v>>1].len - p->te; + (*e) = o_s + g->seq[v>>1].len - p->ts; + } +} + +int64_t gen_ul_pat_seq(uc_block_t *a, path_dp_t *b, all_ul_t *ul, uint64_t idx, asg_t *g, uint32_t v, uint32_t is_backward, uint64_t *rl) +{ + uint32_t m = 0, pv, ai, pi; uint64_t *t, l = 0, tt, s, e; b->pat.n = 0; b->pat_cor.n = 0; + if(is_backward) { + for (ai = idx, l = 0; ai != (uint32_t)-1; ai = a[ai].pidx) { + pv = (((uint32_t)(a[ai].hid))<<1)|((uint32_t)(a[ai].rev)); pv^=1; + if((a[ai].pidx != (uint32_t)-1) && (!ulg_len_check(g, pv, &(a[ai])))) break; + if((b->pat.n > 0) && (!ulg_len_check(g, pv^1, &(a[ai])))) break; + if(pv == v) m++; + kv_pushp(uint64_t, b->pat, &t); + update_l_coord(g, l, pv, &(a[ai]), &s, &e); + kv_push(uint64_t, b->pat_cor, ((s<<32)|e)); + (*t) = pv; (*t) |= (l<<32); + + if(a[ai].pidx != (uint32_t)-1) l += a[ai].pdis; + else l += g->seq[pv>>1].len; + } + } else { + for (ai = idx, pi = 0; ai != (uint32_t)-1; ai = a[ai].aidx) { + pv = (((uint32_t)(a[ai].hid))<<1)|((uint32_t)(a[ai].rev)); + if((a[ai].aidx != (uint32_t)-1) && (!ulg_len_check(g, pv, &(a[ai])))) break; + if((pi > 0) && (!ulg_len_check(g, pv^1, &(a[ai])))) break; + l = ai; pi++; + } + + ai = l; + kv_resize(uint64_t, b->pat, pi); b->pat.n = pi; + kv_resize(uint64_t, b->pat_cor, pi); b->pat_cor.n = pi; + for (l = 0; ai != (uint32_t)-1 && ai >= idx; ai = a[ai].pidx) { + pv = (((uint32_t)(a[ai].hid))<<1)|((uint32_t)(a[ai].rev)); + pi--; if(pv == v) m++; + b->pat.a[pi] = pv; b->pat.a[pi] |= (l<<32); + update_l_coord(g, l, pv, &(a[ai]), &s, &e); + b->pat_cor.a[pi] = ((s<<32)|e); + if(a[ai].pidx != (uint32_t)-1 && a[ai].pidx >= idx) l += a[ai].pdis; + else l += g->seq[pv>>1].len; + } + assert(pi == 0); + + for (ai = 0; ai < b->pat.n; ai++) { + tt = (b->pat.a[ai]>>32) + g->seq[((uint32_t)b->pat.a[ai])>>1].len; + // if(l < tt) { + // fprintf(stderr, "\n[M::%s::] ai::%u, b->pat.n::%u, l::%lu, tt::%lu, beg::%lu\n", + // __func__, ai, (uint32_t)b->pat.n, l, tt, (b->pat.a[ai]>>32)); + // } + assert(l >= tt); + tt = l - tt; + b->pat.a[ai] = (uint32_t)b->pat.a[ai]; + b->pat.a[ai] += (tt<<32); + } + } + (*rl) = l; + + assert(m > 0); + return m; +} + +uint32_t get_arch_len(asg_t *g, uint32_t v, uint32_t w) +{ + + uint32_t i, an; asg_arc_t *av; + av = asg_arc_a(g, v); an = asg_arc_n(g, v); + for (i = 0; i < an; i++) { + if(av[i].del) continue; + if(av[i].v == w) return ((uint32_t)av[i].ul); + } + return (uint32_t)-1; +} + +void gen_ref_pat_seq(asg_t *g, ubuf_t *b, uint64_t rul, double diff_rate, uint32_t v, uint32_t w/**, uint32_t is_debug**/) +{ + int64_t ref_n = b->dp.ref.n, k; uint64_t x, l; uint32_t p, arc_l; + if(ref_n >= 2 && ((uint32_t)b->dp.ref.a[0]) == v && ((uint32_t)b->dp.ref.a[1]) == w) { + // if(is_debug) fprintf(stderr, "+[M::%s::] ref_n::%ld, rul::%lu\n", __func__, ref_n, rul); + if((ref_n) > 0 && ((b->dp.ref.a[ref_n-1]>>32) >= (rul*(1.0+diff_rate)))) { + for (k = ref_n-1; k >= 0; k--) { + if((b->dp.ref.a[k]>>32) < (rul*(1.0+diff_rate))) break; + b->dp.g_flt.a[(uint32_t)b->dp.ref.a[k]] = 0; + } + ref_n = k + 1; b->dp.ref.n = ref_n; + } else { + p = (uint32_t)b->dp.ref.a[ref_n-1]; p ^= 1; p = b->a.a[p].p; + for (l = b->dp.ref.a[ref_n-1]>>32; p != (uint32_t)-1; p = b->a.a[p].p) { + x = p^1; arc_l = get_arch_len(g, ((uint32_t)b->dp.ref.a[b->dp.ref.n-1]), x); + assert(arc_l != (uint32_t)-1); l += arc_l; + if(l >= (rul*(1.0+diff_rate))) break; + b->dp.g_flt.a[x] = 1; x += (l<<32); kv_push(uint64_t, b->dp.ref, x); + } + ref_n = b->dp.ref.n; + } + } else { + // if(is_debug) fprintf(stderr, "-[M::%s::] ref_n::%ld, rul::%lu\n", __func__, ref_n, rul); + b->dp.ref.n = 0; + for (k = 0; k < ref_n; k++) { + b->dp.g_flt.a[(uint32_t)b->dp.ref.a[k]] = 0; + } + + x = v; kv_push(uint64_t, b->dp.ref, x); b->dp.g_flt.a[x] = 1; + for (p = w^1, l = 0; p != (uint32_t)-1; p = b->a.a[p].p) { + x = p^1; arc_l = get_arch_len(g, ((uint32_t)b->dp.ref.a[b->dp.ref.n-1]), x); + assert(arc_l != (uint32_t)-1); l += arc_l; + if(l >= (rul*(1.0+diff_rate))) break; + b->dp.g_flt.a[x] = 1; x += (l<<32); kv_push(uint64_t, b->dp.ref, x); + } + ref_n = b->dp.ref.n; + } +} + +uint32_t quick_check(uc_block_t *a, asg_t *g, path_dp_t *b, double diff_rate, double filter_rate) +{ + uint64_t ref_l = 0, pat_l = 0, mm, min, max, l, k, s, e, match_s, match_e, unmatch_s, unmatch_e, match_l, unmatch_l; + if(b->ref.n == 0 || b->pat.n == 0) return 0; + assert(b->g_flt.a[(uint32_t)b->pat.a[0]]); + for (k = 1; k < b->pat.n; k++) {///k = 0, must be matched + if(b->g_flt.a[(uint32_t)b->pat.a[k]]) break; + } + if(k >= b->pat.n) return 0; + + ref_l = (b->ref.a[b->ref.n-1]>>32) + g->seq[((uint32_t)b->ref.a[b->ref.n-1])>>1].len; + pat_l = (b->pat.a[b->pat.n-1]>>32) + g->seq[((uint32_t)b->pat.a[b->pat.n-1])>>1].len; + mm = MIN(ref_l, pat_l); min = mm * (1.0 - diff_rate); max = mm * (1.0 + diff_rate); + + match_s = match_e = unmatch_s = unmatch_e = (uint64_t)-1; match_l = unmatch_l = 0; + for (k = 0; k < b->pat.n; k++) { + l = (b->pat.a[k]>>32) + g->seq[((uint32_t)b->pat.a[k])>>1].len; + s = b->pat_cor.a[k]>>32; e = (uint32_t)b->pat_cor.a[k]; + if(l <= min) { + if(unmatch_e == (uint64_t)-1 || s >= unmatch_e) { + if(unmatch_e != (uint64_t)-1) unmatch_l += unmatch_e - unmatch_s; + unmatch_s = s; unmatch_e = e; + } else { + if(e > unmatch_e) unmatch_e = e; + } + } + + if(l <= max) { + if(b->g_flt.a[(uint32_t)b->pat.a[k]]) { + if(match_e == (uint64_t)-1 || s >= match_e) { + if(match_e != (uint64_t)-1) match_l += match_e - match_s; + match_s = s; match_e = e; + } else { + if(e > match_e) match_e = e; + } + } + } + if(l > max) break; + } + + if(unmatch_e != (uint64_t)-1) unmatch_l += unmatch_e - unmatch_s; + if(match_e != (uint64_t)-1) match_l += match_e - match_s; + // fprintf(stderr, "[M::%s::] min::%lu, max::%lu, match_l::%lu, unmatch_l::%lu\n", __func__, + // min, max, match_l, unmatch_l); + if(match_l >= unmatch_l*filter_rate) return 1; + return 0; +} + +#define ul_dp_idx(dp, x, y) ((dp).m*(x)+(y)) +#define e_mdp 0 +#define ue_mdp 1 +#define lpat_dp 2 +#define lref_dp 3 +//0->match; 1->mismatch; 2->up (longer pat); 3->left (longer ref) +void init_ul_dp(asg_t *g, path_dp_t *dp) +{ + kv_resize(uint8_t, dp->m_dir, (dp->pat.n+1)*(dp->ref.n+1)); + kv_resize(int64_t, dp->m_score, (dp->pat.n+1)*(dp->ref.n+1)); + dp->n = dp->pat.n+1; dp->m = dp->ref.n+1; + /** + uint64_t k, s, e, sp, ep, sc; + dp->m_dir.a[0] = dp->m_score.a[0] = 0; ///[0, 0] + for (k = 1, sp = ep = (uint64_t)-1; k < dp->n; k++) {///pat; ul read + s = dp->pat.a[k-1]>>32; e = (dp->pat.a[k-1]>>32) + g->seq[(uint32_t)dp->pat.a[k-1]].len; sc = 0; + if(ep == (uint64_t)-1 || s >= ep) { + sc = e - s; sp = s; ep = e; + } else { + if(e > ep) sc = e - ep; + } + if(((uint32_t)dp->pat.a[k-1]) == ((uint32_t)dp->ref.a[0])) { + dp->m_score.a[ul_dp_idx(*dp, k, 0)] = dp->m_dir.a[ul_dp_idx(*dp, k, 0)] = 0; + } else { + dp->m_score.a[ul_dp_idx(*dp, k, 0)] = dp->m_score.a[ul_dp_idx(*dp, k-1, 0)] + sc; + dp->m_dir.a[ul_dp_idx(*dp, k, 0)] = lpat_dp; + } + } + + for (k = 1, sp = ep = (uint64_t)-1; k < dp->m; k++) {///ref; graph + s = dp->ref.a[k-1]>>32; e = (dp->ref.a[k-1]>>32) + g->seq[(uint32_t)dp->ref.a[k-1]].len; sc = 0; + if(ep == (uint64_t)-1 || s >= ep) { + sc = e - s; sp = s; ep = e; + } else { + if(e > ep) sc = e - ep; + } + dp->m_score.a[ul_dp_idx(*dp, 0, k)] = dp->m_score.a[ul_dp_idx(*dp, 0, k-1)] + sc; + dp->m_dir.a[ul_dp_idx(*dp, 0, k)] = lref_dp; + } + **/ + uint64_t k; int64_t min_weight = -1*((int64_t)(0xffffffff)); dp->m_dir.a[0] = dp->m_score.a[0] = 0; ///[0, 0] + + for (k = 1; k < dp->n; k++) {///pat; ul read + if(((uint32_t)dp->pat.a[k-1]) == ((uint32_t)dp->ref.a[0])) { + dp->m_score.a[ul_dp_idx(*dp, k, 0)] = dp->m_dir.a[ul_dp_idx(*dp, k, 0)] = 0; + } else { + dp->m_score.a[ul_dp_idx(*dp, k, 0)] = min_weight; + // dp->m_score.a[ul_dp_idx(*dp, k-1, 0)] - g->seq[((uint32_t)dp->pat.a[k-1])>>1].len; + dp->m_dir.a[ul_dp_idx(*dp, k, 0)] = lpat_dp; + } + } + + for (k = 1; k < dp->m; k++) {///ref; graph + dp->m_score.a[ul_dp_idx(*dp, 0, k)] = min_weight; + // dp->m_score.a[ul_dp_idx(*dp, 0, k-1)] - g->seq[((uint32_t)dp->ref.a[k-1])>>1].len; + dp->m_dir.a[ul_dp_idx(*dp, 0, k)] = lref_dp; + } +} + +uint64_t node_check(uint32_t ul_v, uint32_t g_v, uint32_t ul_v_len, uint32_t g_v_len, double diff_len, uint32_t ul_weight) +{ + if(ul_v != g_v) return 0; + uint32_t x = ((ul_v_len >= g_v_len)? (ul_v_len - g_v_len): (g_v_len - ul_v_len)); + if(x > g_v_len*diff_len) return 0; + if(g_v_len == 0) { + return ul_weight; + } else { + return (((double)(g_v_len-x))/((double)g_v_len))*ul_weight; + } +} + +void print_dp_matrix(path_dp_t *dp) +{ + uint64_t i, j; + for (i = 0; i < dp->n; i++) {//pat + for (j = 0; j < dp->m; j++) {//mat + fprintf(stderr, "%ld<%u>,\t", dp->m_score.a[ul_dp_idx(*dp, i, j)], dp->m_dir.a[ul_dp_idx(*dp, i, j)]); + } + fprintf(stderr, "\n"); + } + +} + +int64_t ul_dp0(bubble_type* bub, asg_t *g, path_dp_t *dp, double pass_thres/**, uint64_t is_debug**/) +{ + if(dp->pat.n < 2 || dp->ref.n < 2) return -1; + uint64_t i, j, d, w; int64_t sc0, sc1, sc2, sc, sc_i, sc_j; + init_ul_dp(g, dp); + // if(is_debug) { + // print_dp_matrix(dp); + // } + for (i = 1; i < dp->n; i++) {//pat + for (j = 1; j < dp->m; j++) {//mat + w = node_check(((uint32_t)dp->pat.a[i-1]), ((uint32_t)dp->ref.a[j-1]), + dp->pat.a[i-1]>>32, dp->ref.a[j-1]>>32, 0.04, g->seq[((uint32_t)dp->pat.a[i-1])>>1].len); + + if(w) { + // if(is_debug) fprintf(stderr, "\n[M::%s::] i::%lu, j::%lu, w::%lu\n", __func__, i, j, w); + dp->m_score.a[ul_dp_idx(*dp, i, j)] = dp->m_score.a[ul_dp_idx(*dp, i-1, j-1)] + w; + dp->m_dir.a[ul_dp_idx(*dp, i, j)] = e_mdp; + } else { + sc0 = dp->m_score.a[ul_dp_idx(*dp, i-1, j-1)] - + (int64_t)(MIN(g->seq[((uint32_t)dp->pat.a[i-1])>>1].len, g->seq[((uint32_t)dp->ref.a[j-1])>>1].len)); + sc1 = dp->m_score.a[ul_dp_idx(*dp, i-1, j)] - (int64_t)(g->seq[((uint32_t)dp->ref.a[j-1])>>1].len); + sc2 = dp->m_score.a[ul_dp_idx(*dp, i, j-1)] - (int64_t)(g->seq[((uint32_t)dp->pat.a[i-1])>>1].len); + + d = ue_mdp; sc = sc0; + if(sc < sc1) { + sc = sc1; d = lref_dp; + } + if(sc < sc2) { + sc = sc2; d = lpat_dp; + } + dp->m_score.a[ul_dp_idx(*dp, i, j)] = sc; + dp->m_dir.a[ul_dp_idx(*dp, i, j)] = d; + } + } + } + + sc = 0; sc_i = sc_j = -1; + for (j = 1, i = dp->n - 1; j < dp->m; j++) { + if(sc_i < 0 || sc < dp->m_score.a[ul_dp_idx(*dp, i, j)]) { + sc_i = i; sc_j = j; sc = dp->m_score.a[ul_dp_idx(*dp, i, j)]; + } + } + for (i = 1, j = dp->m - 1; i < dp->n; i++) { + if(sc_i < 0 || sc < dp->m_score.a[ul_dp_idx(*dp, i, j)]) { + sc_i = i; sc_j = j; sc = dp->m_score.a[ul_dp_idx(*dp, i, j)]; + } + } + + // if(is_debug) { + // fprintf(stderr, "[M::%s::] sc_i::%ld, sc_j::%ld, sc::%ld\n", __func__, sc_i, sc_j, sc); + // // print_dp_matrix(dp); + // } + + if (sc_i < 0 || sc_i == 1 || sc_j == 1) return 0; + int64_t match_all = 0, unmatch_all = 0, het_match_all = 0, het_unmatch_all = 0, het_occ = 0; + //pat_s = sc_i - 1; ref_s = ref_e = sc_j - 1; + for(i = sc_i, j = sc_j; i > 0 && j > 0;) { + d = dp->m_dir.a[ul_dp_idx(*dp, i, j)]; + //pat_s = i - 1; ref_s = j - 1; + // if(is_debug) fprintf(stderr, "[M::%s::] i::%lu, j::%lu, dir::%lu\n", __func__, i, j, d); + // if(pat_s == 0) continue; + if(d == e_mdp) { + // if(i != 1 || j != 1) {//skip (i == 1 && j == 1) + match_all += g->seq[((uint32_t)dp->pat.a[i-1])>>1].len; + if(!IF_HOM((((uint32_t)dp->pat.a[i-1])>>1), *bub)) { + het_match_all += g->seq[((uint32_t)dp->pat.a[i-1])>>1].len; + het_occ++; + } + i--; j--; + // if(is_debug) { + // fprintf(stderr, "[M::%s::] PAT::utg%.6dl(%c::len->%lu), REF::utg%.6dl(%c::len->%lu)\n", __func__, + // ((((uint32_t)dp->pat.a[i])>>1))+1, "+-"[((uint32_t)dp->pat.a[i])&1], dp->pat.a[i]>>32, + // ((((uint32_t)dp->ref.a[j])>>1))+1, "+-"[((uint32_t)dp->ref.a[j])&1], dp->ref.a[j]>>32); + // } + } + if(d == ue_mdp) { + unmatch_all += g->seq[((uint32_t)dp->pat.a[i-1])>>1].len; + if(!IF_HOM((((uint32_t)dp->pat.a[i-1])>>1), *bub)) { + het_unmatch_all += g->seq[((uint32_t)dp->pat.a[i-1])>>1].len; + } + i--; j--; + } + if(d == lpat_dp) { + j--; + } + if(d == lref_dp) { + unmatch_all += g->seq[((uint32_t)dp->pat.a[i-1])>>1].len; + if(!IF_HOM((((uint32_t)dp->pat.a[i-1])>>1), *bub)) { + het_unmatch_all += g->seq[((uint32_t)dp->pat.a[i-1])>>1].len; + } + i--; + } + } + // if(is_debug) { + // fprintf(stderr, "[M::%s::] het_occ::%ld, match_all::%ld, unmatch_all::%ld, het_match_all::%ld, het_unmatch_all::%ld\n", __func__, + // het_occ, match_all, unmatch_all, het_match_all, het_unmatch_all); + // } + + assert((j == 0) && (((uint32_t)dp->pat.a[i]) == ((uint32_t)dp->ref.a[j]))); + ///skip the beg node + match_all -= g->seq[((uint32_t)dp->pat.a[i])>>1].len; + if(!IF_HOM((((uint32_t)dp->pat.a[i])>>1), *bub)) { + het_match_all -= g->seq[((uint32_t)dp->pat.a[i])>>1].len; het_occ--; + } + + + if(het_occ < 1 || match_all == 0 || het_match_all == 0) return -1; + if(match_all <= (match_all + unmatch_all)*pass_thres) return -1; + if(het_match_all <= (het_match_all + het_unmatch_all)*pass_thres) return -1; + return (((double)het_match_all)/((double)(het_match_all + het_unmatch_all)))* + (((uint32_t)dp->pat_cor.a[i]) - (dp->pat_cor.a[i]>>32)); +} +uint64_t ul_dp(bubble_type* bub, all_ul_t *ul, ubuf_t *b, uint64_t *a, int64_t a_n, asg_t *g, uint32_t v, uint32_t w, uint32_t is_backward) +{ + int64_t i, m, a_i = -1; uc_block_t *p; uint32_t gv; uint64_t url; int64_t we; + for (i = m = 0; i < a_n; i++) { + p = &(ul->a[a[i]>>32].bb.a[(uint32_t)(a[i])]); + assert((!p->base)&&(p->hid == (v>>1))&&(p->el)/**&&(p->pchain)**/); + if(!p->pchain) continue; + gv = (((uint32_t)(p->hid))<<1)|((uint32_t)(p->rev)); + if(is_backward) { + gv ^= 1; + if(p->pidx == (uint32_t)-1) continue; + } else { + if(p->aidx == (uint32_t)-1) continue; + } + if(gv != v) continue; + // fprintf(stderr, "[M::%s::] i:%ld, a_n:%lu, v:%u, w:%u, ulid:%lu, uidx:%u, is_backward:%u\n", + // __func__, i, a_n, v, w, a[i]>>32, (uint32_t)(a[i]), is_backward); + // if(v == 6 && w == 5 && (a[i]>>32) == 95) { + // int64_t z; + // for (z = 0; z < ul->a[a[i]>>32].bb.n; z++) { + // fprintf(stderr, "[M::%s::vid->%u] qs:%u, qe:%u, qlen:%u, ts:%u, te:%u, tlen:%u\n", __func__, + // (((uint32_t)(ul->a[a[i]>>32].bb.a[z].hid))<<1)|((uint32_t)(ul->a[a[i]>>32].bb.a[z].rev)), + // ul->a[a[i]>>32].bb.a[z].qs, ul->a[a[i]>>32].bb.a[z].qe, ul->a[a[i]>>32].rlen, + // ul->a[a[i]>>32].bb.a[z].ts, ul->a[a[i]>>32].bb.a[z].te, g->seq[ul->a[a[i]>>32].bb.a[z].hid].len); + // } + // } + if(ulg_len_check(g, v, p) == 0) continue;///too short + m++; a_i = i; + } + + if(m == 0) return 0; + // uint32_t is_debug = (v == 231); + // if(is_debug) { + // fprintf(stderr, "\n+[M::%s::] m->%ld, ulid->%lu\n", __func__, m, a[a_i]>>32); + // } + + p = &(ul->a[a[a_i]>>32].bb.a[(uint32_t)(a[a_i])]);//longest one + m = gen_ul_pat_seq(ul->a[a[a_i]>>32].bb.a, &(b->dp), ul, (uint32_t)(a[a_i]), g, v, is_backward, &url); + assert(b->dp.pat.n > 0); + + // if(is_debug) { + // for (i = 0; i < (int64_t)b->dp.pat.n; i++) { + // fprintf(stderr, "pat[M::%s::] utg%.6dl(%c), d::%lu\n", __func__, + // (((uint32_t)b->dp.pat.a[i])>>1)+1, "+-"[((uint32_t)b->dp.pat.a[i])&1], b->dp.pat.a[i]>>32); + // } + // } + if(b->dp.pat.n == 1) return 0; + gen_ref_pat_seq(g, b, url, UL_TRAV_HERATE, v, w/**, is_debug**/); + + // if(is_debug) { + // for (i = 0; i < (int64_t)b->dp.ref.n; i++) { + // fprintf(stderr, "ref[M::%s::] utg%.6dl(%c), d::%lu\n", __func__, + // (((uint32_t)b->dp.ref.a[i])>>1)+1, "+-"[((uint32_t)b->dp.ref.a[i])&1], b->dp.ref.a[i]>>32); + // } + // } + + if(!quick_check(ul->a[a[a_i]>>32].bb.a, g, &(b->dp), UL_TRAV_HERATE, UL_TRAV_FT_RATE)) return 0; + we = ul_dp0(bub, g, &(b->dp), 0.9/**, is_debug**/); + if(we < 0) return 0; + return we; +} + +uint64_t get_eul_weight(bubble_type* bub, uint32_t v, uint32_t w, uinfo_t *g_idx, ubuf_t *b, all_ul_t *ul, asg_t *g) +{ + uint64_t *a, to = 0; int64_t a_n, l, k; + a = ul->ridx.occ.a + ul->ridx.idx.a[v>>1]; + a_n = ul->ridx.idx.a[(v>>1)+1] - ul->ridx.idx.a[v>>1]; + b->dp.pat.n = b->dp.pat_cor.n = b->dp.ref.n = 0; + for (l = 0, k = 1; k <= a_n; k++) { + if((k == a_n) || ((a[k]>>32) != (a[l]>>32))) { + // if((v == 106 && w == 103) || (v == 24 && w == 23)) { + // fprintf(stderr, "+[M::%s::] l->%ld, k->%ld, a_n->%ld\n", __func__, l, k, a_n); + // } + to += ul_dp(bub, ul, b, a + l, k - l, g, v, w, 1); + // if((v == 106 && w == 103) || (v == 24 && w == 23)) { + // fprintf(stderr, "-[M::%s::] l->%ld, k->%ld, a_n->%ld\n", __func__, l, k, a_n); + // } + to += ul_dp(bub, ul, b, a + l, k - l, g, v, w, 0); + // if((v == 106 && w == 103) || (v == 24 && w == 23)) { + // fprintf(stderr, "*[M::%s::] l->%ld, k->%ld, a_n->%ld\n", __func__, l, k, a_n); + // } + l = k; + } + } + + return to; +} + +void extract_paths(ubuf_t *b, ma_utg_v *gu, ul_path_t *res, uint32_t src, uint32_t dest) +{ + uint32_t i, v; uinfo_t *t; uint64_t m = 0, mi = (uint64_t)-1, avn; ///kv_resize(uinfo_srt_t, b->srt, b->b.n); b->srt.n = b->b.n; + for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices + t = &b->a.a[b->b.a[i]]; ///memset(t, 0, sizeof(*(t))); + //b->srt.a[i].c = ((uint64_t)-1) - t->c; b->srt.a[i].i = i; + if(m < t->c) { + m = t->c; mi = i; ///b->b.a[i]; + } + } + //radix_sort_uinfo_srt_t_c(b->srt.a, b->srt.a + b->srt.n); + // fprintf(stderr, "\n[M::%s::]->beg\n", __func__); + if(mi != (uint64_t)-1) { + ma_utg_t *p; kv_pushp(ma_utg_t, res->buf_ug->u, &p); + p->len = p->circ = p->n = 0; p->start = p->end = UINT32_MAX; p->a = NULL; p->s = NULL; + for(v = b->b.a[mi], avn = 0; v != (uint32_t)-1; v = b->a.a[v].p) { + if(b->a.a[v].c == 0 && avn == 0) break; + b->us.a[v>>1] = 1; + if(v != src && v != dest) kv_push(uint64_t, (*p), v); + // fprintf(stderr, "[M::%s::] utg%.6d%c(%c)\tC::%lu\n", __func__, (v>>1)+1, "lc"[gu->a[v>>1].circ], "+-"[v&1], b->a.a[v].c); + avn = b->a.a[v].c; + } + + if(p->n == 0) res->buf_ug->u.n = 0; + } +} + +uint32_t hc_simple_traversal(bubble_type* bub, asg_t *g, ma_utg_v *gu, ubuf_t *b, all_ul_t *ul, uint32_t src, uint32_t dest, ul_path_t *res) +{ + uint32_t v, nv, i, w, n_pending = 0, is_update = 0; uint64_t l, d, c, cc, nc, c_nc; asg_arc_t *av; uinfo_t *t; + if (g->seq[src>>1].del || g->seq[dest>>1].del) return 0; + clear_ubuf_t(b, g, ul); b->a.a[src].p = (uint32_t)-1; + kv_push(uint32_t, b->S, src); + while (b->S.n > 0) { + v = kv_pop(b->S); d = b->a.a[v].d; c = b->a.a[v].c; nc = b->a.a[v].nc; + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + for (i = 0; i < nv; ++i) { + if (av[i].del || g->seq[av[i].v>>1].del) continue; + w = av[i].v; l = (uint32_t)av[i].ul; t = &b->a.a[w]; + kv_push(uint32_t, b->e, ((g->idx[v]>>32)+i)); ///push the edge + // fprintf(stderr, "+[M::%s::] utg%.6d%c(%c:%u) -> utg%.6d%c(%c:%u)\n", __func__, + // (v>>1)+1, "lc"[gu->a[v>>1].circ], "+-"[v&1], v, + // (w>>1)+1, "lc"[gu->a[w>>1].circ], "+-"[w&1], w); + cc = c + get_eul_weight(bub, w^1, v^1, b->a.a, b, ul, g); + c_nc = nc + (gu?gu->a[w>>1].n:1); + // fprintf(stderr, "-[M::%s::] utg%.6d%c(%c:%u) -> utg%.6d%c(%c:%u)\n", __func__, + // (v>>1)+1, "lc"[gu->a[v>>1].circ], "+-"[v&1], v, + // (w>>1)+1, "lc"[gu->a[w>>1].circ], "+-"[w&1], w); + + if (t->s == 0) {///a new node + kv_push(uint32_t, b->b, w); // save it for revert + t->p = v; t->s = 1; t->d = d + l; + t->r = get_arcs(g, w^1, NULL, 0); + t->nc = c_nc; t->c = cc; + ++n_pending; + } else { + is_update = 0; + if(b->us.a[t->p>>1] == b->us.a[v>>1]) { + if(cc > t->c) is_update = 1; + if(cc == t->c && c_nc > t->nc) is_update = 1; + } else { + if(b->us.a[t->p>>1]) is_update = 1; + else is_update = 0; + } + + if(is_update) { + t->p = v; t->s = 1; t->d = d + l; t->nc = c_nc; t->c = cc; + } + } + + if (--(t->r) == 0) { + if(get_arcs(g, w, NULL, 0) > 0) kv_push(uint32_t, b->S, w); + --n_pending; + if(w == dest && n_pending == 0) goto pp_end; + } + } + } + pp_end: + extract_paths(b, gu, res, src, dest); + for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices + t = &b->a.a[b->b.a[i]]; memset(t, 0, sizeof(*(t))); + } + return 1; +} + + +// void clean_path_g(ul_path_t *p, asg_t *ref) +// { +// if(!(p->pg)) CALLOC(p->pg, 1); +// free(p->pg->idx); p->pg->idx = 0; p->pg->is_srt = 0; p->pg->is_symm = 0; +// REALLOC(p->pg->seq, ref->n_seq); p->pg->n_seq = ref->n_seq; p->pg->n_arc = 0; +// memcpy(p->pg->seq, ref->seq, p->pg->n_seq*sizeof(*(p->pg->seq))); +// asg_cleanup(p->pg); +// } + +void build_sub_graph(ul_resolve_t *uidx, ul_path_t *res) +{ + +} + +uint64_t get_chain_ul_cov(ul_resolve_t *uidx, uint64_t *a, uint64_t a_n, uint64_t bub_id, uint64_t inner_beg) +{ + if(a_n <= 0) return 0; + uint64_t w = 0, bid, rev; uint32_t rt, beg, sink; + uidx->path.buf.n = 0; + // kv_pushp(ul_sub_path_t, uidx->path, &p); memset(p, 0, sizeof((*p))); + // p->bid = bub_id; p->beg = inner_beg; p->occ = a_n; + + bid = a[0]>>33; rev = (a[0]>>32)&1; + get_bubbles(uidx->bub, bid, rev==0?&rt:NULL, rev==1?&rt:NULL, NULL, NULL, NULL); beg = rt; + + bid = a[a_n-1]>>33; rev = (a[a_n-1]>>32)&1; + get_bubbles(uidx->bub, bid, rev==1?&rt:NULL, rev==0?&rt:NULL, NULL, NULL, NULL); sink = rt^1; + + + + + + + + + // fprintf(stderr, "\n[M::%s::] # bubbles->%lu\n", __func__, a_n); + // uint64_t k + // for (k = 0; k < a_n; k++) { + // bid = a[k]>>33; rev = (a[k]>>32)&1; + // get_bubbles(uidx->bub, bid, rev==0?&rt:NULL, rev==1?&rt:NULL, NULL, NULL, NULL); + // // w += get_bub_ul_cov(a[k]>>33, ug, idx); + // fprintf(stderr, "[M::%s::] utg%.6d%c(%c)\n", __func__, (rt>>1)+1, "lc"[uidx->l1_ug->u.a[rt>>1].circ], "+-"[rt&1]); + // } + // if (a_n > 0) { + // bid = a[k-1]>>33; rev = (a[k-1]>>32)&1; + // get_bubbles(uidx->bub, bid, rev==1?&rt:NULL, rev==0?&rt:NULL, NULL, NULL, NULL); rt ^= 1; + // fprintf(stderr, "[M::%s::] utg%.6d%c(%c)\n", __func__, (rt>>1)+1, "lc"[uidx->l1_ug->u.a[rt>>1].circ], "+-"[rt&1]); + // } + // clean_path_g(&(uidx->path), uidx->l1_ug->g); + ma_utg_t *p; + kv_resize(uint8_t, uidx->buf.us, uidx->l1_ug->g->n_seq); uidx->buf.us.n = uidx->l1_ug->g->n_seq; memset(uidx->buf.us.a, 0, uidx->buf.us.n*sizeof(*(uidx->buf.us.a))); + hc_simple_traversal(uidx->bub, uidx->l1_ug->g, &(uidx->l1_ug->u), &(uidx->buf), uidx->idx, beg, sink, &(uidx->path)); + hc_simple_traversal(uidx->bub, uidx->l1_ug->g, &(uidx->l1_ug->u), &(uidx->buf), uidx->idx, beg, sink, &(uidx->path)); + + kv_pushp(ma_utg_t, uidx->path.buf_ug->u, &p); p->len = p->circ = p->n = 0; + p->start = p->end = UINT32_MAX; p->a = NULL; p->s = NULL; kv_push(uint64_t, *p, beg); + + kv_pushp(ma_utg_t, uidx->path.buf_ug->u, &p); p->len = p->circ = p->n = 0; + p->start = p->end = UINT32_MAX; p->a = NULL; p->s = NULL; kv_push(uint64_t, *p, sink); + + return w; +} + +void phrase_exact_chains(ul_resolve_t *uidx, uint64_t *a, int64_t a_n, asg_t *bg, uint64_t bub_id, uint64_t inner_beg) +{ + int64_t k, l; + for (k = l = 0; k < a_n; k++) { + if(k+1 >= a_n) continue; + if((arc_first(bg, a[k]>>32)).el) continue; + get_chain_ul_cov(uidx, a+l, k+1-l, bub_id, inner_beg + l); + l = k + 1; + } + if(l < a_n) get_chain_ul_cov(uidx, a+l, a_n-l, bub_id, inner_beg + l); +} +void resolve_dip_bub_chains(ul_resolve_t *uidx) +{ + bubble_type *bub = uidx->bub; + uint32_t i; int32_t k, l, un; ma_utg_t *u; + for (i = 0; i < bub->b_ug->u.n; i++) { + u = &(bub->b_ug->u.a[i]); un = u->n; + if(un <= 1) continue;///single bubble + for (l = -1, k = 0; k <= un; k++) { + if((k == un) || ((u->a[k]>>33) >= bub->f_bub)) {///not a complete bubble + // if(k < un) { + // fprintf(stderr, "++++++[M::%s::] l->%d, k->%d, bid->%lu, f_bub->%lu\n", + // __func__, l, k, (u->a[k]>>33), bub->f_bub); + // } + if(k - l > 1) {///at least a complete bubble + // fprintf(stderr, "\n[M::%s::] l->%d, k->%d\n", __func__, l, k); + phrase_exact_chains(uidx, u->a+l+1, k-l-1, bub->b_g, i, l+1); + } + l = k; + } + } + + } +} + +/** +static void fill_ul_path_srt_t(void *data, long i, int tid) // callback for kt_for() +{ + all_ul_t *idx = ((ul_resolve_t *)data)->idx; + ul_path_srt_t *idx_srt = &(((ul_resolve_t *)data)->psrt); + uc_block_t *a = NULL; uc_block_t *p; int64_t k, a_n, b_n, occ; uint64_t *b = NULL; + + a = idx->a[i].bb.a; a_n = idx->a[i].bb.n; + if(idx_srt->ul_n == 0) { + for (k = occ = 0; k < a_n; k++) { + p = &(a[k]); + if(p->base || (!p->el) || (!p->pchain)) continue; + occ++; + } + if(occ == 1) occ = 0;///UL read is too short + idx_srt->idx.a[i] = occ; + } else { + b_n = idx_srt->idx.a[i]>>33; b = idx_srt->srt.a + (uint32_t)idx_srt->idx.a[i]; + if(b_n == 0) return;//UL covers 0 or 1 node; not useful + for (k = occ = 0; k < a_n; k++) { + p = &(a[k]); + if(p->base || (!p->el) || (!p->pchain)) continue; + b[occ] = p->hid; b[occ] <<= 1; b[occ] |= p->rev; b[occ] <<= 32; b[occ] += (uint64_t)k; + occ++; + } + assert(occ == b_n); + radix_sort_srt64(b, b + b_n); + for (k = 1; k < b_n; k++) { + if((b[k]>>32) == (b[k-1]>>32)) break; + } + + if(k < b_n) idx_srt->idx.a[i] |= (uint64_t)(0x100000000); + } +} + + +void init_ul_path_srt_t(ul_resolve_t *p, all_ul_t *idx, ul_path_srt_t *idx_srt) +{ + idx_srt->ul_n = 0; + idx_srt->idx.n = idx_srt->idx.m = idx->n; CALLOC(idx_srt->idx.a, idx_srt->idx.n); + + kt_for(asm_opt.thread_num, fill_ul_path_srt_t, p, idx->n); + uint64_t l, k; + for (k = l = 0; k < idx_srt->idx.n; k++) { + idx_srt->idx.a[k] <<= 33; idx_srt->idx.a[k] += l; l += (idx_srt->idx.a[k]>>33); + } + MALLOC(idx_srt->srt.a, l); idx_srt->srt.n = idx_srt->srt.m = l; + + idx_srt->ul_n = idx->n; + kt_for(asm_opt.thread_num, fill_ul_path_srt_t, p, idx->n); +} +**/ + +uint64_t ug_occ_w(uint64_t is, uint64_t ie, ma_utg_t *u) +{ + if(is == 0 && ie == u->len) return u->n; + uint64_t l, i, us, ue, occ; + for (i = l = occ = 0; i < u->n; i++) { + us = l; ue = l + Get_READ_LENGTH(R_INF, (u->a[i]>>33)); + if(is <= us && ie >= ue) occ++; + if(us >= ie) break; + l += (uint32_t)u->a[i]; + } + return occ; +} + +// static void gen_ul_str_idx_t(void *data, long i, int tid) // callback for kt_for() +// { +// all_ul_t *idx = ((ul_resolve_t *)data)->idx; +// ma_ug_t *ug = ((ul_resolve_t *)data)->l1_ug; +// ul_str_t *str = &(((ul_resolve_t *)data)->pstr.str.a[i]); +// // uint64_t *str_idx = ((ul_resolve_t *)data)->pstr.idx.a; +// uc_block_t *a = NULL; uc_block_t *xk; uint64_t t; +// uint64_t k, a_n; + +// str->n = str->m = 0; str->a = NULL; +// a = idx->a[i].bb.a; a_n = idx->a[i].bb.n; +// kv_resize(uint64_t, *str, idx->a[i].bb.n); +// for (k = 0; k < a_n; k++) { +// xk = &(a[k]); +// if(xk->base || (!xk->el) || (!xk->pchain)) continue; +// t = k; t <<= 32; t += (xk->hid<<1); t += xk->rev; +// kv_push(uint64_t, *str, t); +// } +// str->cn = str->n; str->is_cir = 0; +// } + +void init_ul_str_idx_t(ul_resolve_t *p) +{ + uint64_t k, l, i, m, *a, a_n, x_n; uc_block_t *x_a; + all_ul_t *idx = p->idx; ul_str_idx_t *str = &(p->pstr); ul_str_t *z; uint32_t v0, v1; + MALLOC(str->str.a, idx->n); str->str.n = str->str.m = idx->n; + CALLOC(str->idx.a, p->l1_ug->u.n+1); str->idx.n = str->idx.m = p->l1_ug->u.n+1; + for (i = 0; i < str->str.n; i++) { + x_a = idx->a[i].bb.a; x_n = idx->a[i].bb.n; + memset(&(str->str.a[i]), 0, sizeof(str->str.a[i])); + for (k = 0; k < x_n; k++) { + if(x_a[k].base || (!x_a[k].el) || (!x_a[k].pchain)) continue; + l = k; l <<= 32; l += (x_a[k].hid<<1); l += x_a[k].rev; str->idx.a[x_a[k].hid]++; + kv_push(uint64_t, str->str.a[i], l); + } + str->str.a[i].cn = str->str.a[i].n; str->str.a[i].is_cir = 0; + } + + // kt_for(asm_opt.thread_num, gen_ul_str_idx_t, p, str->str.n); + for (k = l = 0; k < str->idx.n; k++) { + m = str->idx.a[k]; + str->idx.a[k] = l; + l += m; + } + + MALLOC(str->occ.a, l); str->occ.n = str->occ.m = l; + for (k = 0; k < p->l1_ug->u.n; k++) { + a = str->occ.a + str->idx.a[k]; + a_n = str->idx.a[k+1] - str->idx.a[k]; + if(a_n) a[a_n-1] = 0; + } + + for (k = 0; k < str->str.n; k++) { + z = &(str->str.a[k]); + for (i = 0; i < z->cn; i++) { + a = str->occ.a + str->idx.a[((uint32_t)z->a[i])>>1]; + a_n = str->idx.a[(((uint32_t)z->a[i])>>1)+1] - str->idx.a[((uint32_t)z->a[i])>>1]; + if(a_n) { + if(a[a_n-1] == a_n-1) { + m = a_n-1; a[a_n-1] = (k<<32)|i; + } else { + m = a[a_n-1]; a[a[a_n-1]++] = (k<<32)|i; + } + + v0 = (uint32_t)z->a[(uint32_t)a[m]]; + while (m > 0 && z->is_cir == 0) { + m--; + if((a[m]>>32) != k) break; + v1 = (uint32_t)z->a[(uint32_t)a[m]]; + if(v0 == v1) z->is_cir = 1; + } + } + } + } +} + +ul_resolve_t *init_ul_resolve_t(asg_t *sg, ma_ug_t *init_ug, bubble_type* bub, all_ul_t *idx, uint8_t *r_het) +{ + ul_resolve_t *p = NULL; CALLOC(p, 1); + p->sg = sg; p->init_ug = init_ug; p->bub = bub; p->idx = idx; p->r_het = r_het; + p->l1_ug = copy_untig_graph(p->init_ug); + init_integer_ml_t(&p->str_b, p, asm_opt.thread_num); + // init_ul_path_srt_t(p); + init_ul_str_idx_t(p); + return p; +} + +void print_bubble_gfa(FILE *fp, bubble_type *bub, const char* utg_pre, const char* bub_pre, const char* chain_pre) +{ + uint32_t i, k, m, *a, n, beg, sink, x; ma_utg_t *p; uint64_t occ; + ma_ug_t *b_ug = bub->b_ug; char name[32], bname[32]; uint8_t *f; CALLOC(f, bub->ug->u.n); + for (i = 0; i < b_ug->u.n; i++) { + p = &b_ug->u.a[i]; + if(p->n == 0) continue; + for (k = occ = 0; k < p->n; k++){ + x = p->a[k]>>33; + get_bubbles(bub, x, &beg, &sink, &a, &n, NULL); + + for (m = 0; m < n; m++) { + occ += bub->ug->u.a[a[m]>>1].n; f[a[m]>>1] = 1; + } + if(beg != (uint32_t)-1 && f[beg>>1] == 0) { + occ += bub->ug->u.a[beg>>1].n; f[beg>>1] = 1; + } + if(sink != (uint32_t)-1 && f[sink>>1] == 0) { + occ += bub->ug->u.a[sink>>1].n; f[sink>>1] = 1; + } + } + + sprintf(name, "%s%.6d%c", chain_pre, i + 1, "lc"[p->circ]); + fprintf(fp, "S\t%s\t*\tLN:i:%lu\n", name, occ); + for (k = 0; k < p->n; k++) { + x = p->a[k]>>33; + sprintf(bname, "%s%.6d", bub_pre, x + 1); + fprintf(fp, "B\t%s\t%c\tcid:i:%s\tsm:%c\n", bname, "+-"[(p->a[k]>>32)&1], name, "01"[xf_bub]); + + get_bubbles(bub, x, &beg, &sink, &a, &n, NULL); + if(beg != (uint32_t)-1) { + fprintf(fp, "U\t%s%.6d%c\t%c\tcid:i:%s\tbid:b:%s\thom:%c\n", + utg_pre, (beg>>1)+1, "lc"[bub->ug->u.a[(beg>>1)].circ], "+-"[beg&1], name, bname, "10"[IF_HOM((beg>>1), *bub)]); + } + + if(sink != (uint32_t)-1) { + fprintf(fp, "U\t%s%.6d%c\t%c\tcid:i:%s\tbid:s:%s\thom:%c\n", + utg_pre, (sink>>1)+1, "lc"[bub->ug->u.a[(sink>>1)].circ], "+-"[sink&1], name, bname, "10"[IF_HOM((sink>>1), *bub)]); + } + for (m = 0; m < n; m++) { + occ += bub->ug->u.a[a[m]>>1].n; + fprintf(fp, "U\t%s%.6d%c\t%c\tcid:i:%s\tbid:c:%s\thom:%c\n", + utg_pre, (a[m]>>1)+1, "lc"[bub->ug->u.a[(a[m]>>1)].circ], "+-"[a[m]&1], name, bname, "10"[IF_HOM((a[m]>>1), *bub)]); + } + } + } + + asg_arc_t* au = NULL; + uint32_t nu, u, v, j; + for (i = 0; i < b_ug->u.n; ++i) { + if(b_ug->u.a[i].m == 0) continue; + if(b_ug->u.a[i].circ) + { + fprintf(fp, "L\t%s%.6dc\t+\t%s%.6dc\t+\t%dM\tL1:i:%d\n", + chain_pre, i+1, chain_pre, i+1, 0, 0); + fprintf(fp, "L\t%s%.6dc\t-\t%s%.6dc\t-\t%dM\tL1:i:%d\n", + chain_pre, i+1, chain_pre, i+1, 0, 0); + } + u = i<<1; + au = asg_arc_a(b_ug->g, u); + nu = asg_arc_n(b_ug->g, u); + for (j = 0; j < nu; j++) + { + if(au[j].del) continue; + v = au[j].v; + fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\n", + chain_pre, (u>>1)+1, "lc"[b_ug->u.a[u>>1].circ], "+-"[u&1], + chain_pre, (v>>1)+1, "lc"[b_ug->u.a[v>>1].circ], "+-"[v&1], 0, 0); + } + + + u = (i<<1) + 1; + au = asg_arc_a(b_ug->g, u); + nu = asg_arc_n(b_ug->g, u); + for (j = 0; j < nu; j++) + { + if(au[j].del) continue; + v = au[j].v; + fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\n", + chain_pre, (u>>1)+1, "lc"[b_ug->u.a[u>>1].circ], "+-"[u&1], + chain_pre, (v>>1)+1, "lc"[b_ug->u.a[v>>1].circ], "+-"[v&1], 0, 0); + } + } + + for (i = 0; i < bub->ug->u.n; i++) { + if(f[i]) continue; + fprintf(fp, "U\t%s%.6d%c\t+\tcid:i:*\tbid:c:*\thom:%c\n", + utg_pre, i+1, "lc"[bub->ug->u.a[i].circ], "10"[IF_HOM(i, *bub)]); + } + free(f); +} + +void print_uls_ovlp(FILE *fp, all_ul_t *uls, const char* utg_pre, ma_ug_t *ug) +{ + ul_vec_t *p = NULL; nid_t *z = NULL; uc_block_t *m = NULL; uint64_t k, i; uint32_t a, la, occ; + kvec_t(uint8_t) f; kv_init(f); + for (k = 0; k < uls->n; k++) { + z = &(uls->nid.a[k]); + p = &(uls->a[k]); + for (i = 0; i < p->bb.n; i++) { + m = &(p->bb.a[i]); + if(m->base || (!m->el)) continue; + if(m->pchain) break; + } + if(i >= p->bb.n) continue; + kv_resize(uint8_t, f, p->bb.n); memset(f.a, 0, sizeof(*(f.a))*p->bb.n); + fprintf(fp, ">\t%.*s\n", (int32_t)z->n, z->a); + for (i = 0; i < p->bb.n; i++) { + if(f.a[i]) continue; + m = &(p->bb.a[i]); + if(m->base || (!m->el) || (!m->pchain)) continue; + for (a = i, la = i, occ = 0; a != (uint32_t)-1; a = p->bb.a[a].aidx) { + la = a; f.a[a] = 1; occ++; + } + fprintf(fp, "C\tLEN:%u\tS:%u\tE:%u\tOCC:%u\n", p->rlen, p->bb.a[i].qs, p->bb.a[la].qe, occ); + for (a = i; a != (uint32_t)-1; a = p->bb.a[a].aidx) { + m = &(p->bb.a[a]); + fprintf(fp, "%s%.6d%c[%c::", utg_pre, m->hid+1, "lc"[ug->u.a[(m->hid>>1)].circ], "+-"[m->rev]); + if(m->pdis == (uint32_t)-1) fprintf(fp, "*]\t"); + else fprintf(fp, "%u]\t", m->pdis); + } + fprintf(fp, "\n"); + } + } + + free(f.a); +} + +void print_debug_ul(const char* o_n, ma_ug_t *ug, asg_t *rg, const ma_sub_t *cov, ma_hit_t_alloc* src, R_to_U* ridx, +bubble_type *bub, all_ul_t *uls) +{ + char* gfa_name = (char*)malloc(strlen(o_n)+50); FILE *fn; + + if(ug && rg && cov && src && ridx) { + uint32_t k, m, v, w, an; asg_arc_t *p, *av; + for (k = 0; k < ug->g->n_arc; k++) { + p = &(ug->g->arc[k]); v = p->v^1; w = (p->ul>>32)^1; + av = asg_arc_a(ug->g, v); an = asg_arc_n(ug->g, v); + for (m = 0; m < an; m++) { + if(av[m].del == 0 && av[m].v == w) break; + } + assert(m < an && av[m].ou == p->ou); + } + sprintf(gfa_name, "%s.r_utg.noseq.gfa", o_n); fn = fopen(gfa_name, "w"); + ma_ug_print_simple(ug, rg, cov, src, ridx, "utg", fn); + fclose(fn); + } + + if(bub) { + sprintf(gfa_name, "%s.bub.noseq.gfa", o_n); fn = fopen(gfa_name, "w"); + print_bubble_gfa(fn, bub, "utg", "btg", "ctg"); + fclose(fn); + } + + if(uls && ug) { + sprintf(gfa_name, "%s.uls.ovlp", o_n); fn = fopen(gfa_name, "w"); + print_uls_ovlp(fn, uls, "utg", ug); + fclose(fn); + } + + free(gfa_name); + fprintf(stderr, "[M::%s::] done\n", __func__); + exit(1); +} + +int64_t normlize_gdis(ma_ug_t *ug, uc_block_t *i, uc_block_t *k) +{ + int64_t i_len, k_len; + ///i > k + if(!i->rev) { + i_len = i->qe + (ug->g->seq[i->hid].len - i->te); + } else { + i_len = i->qe + i->ts; + } + + if(!k->rev) { + k_len = k->qe + (ug->g->seq[k->hid].len - k->te); + } else { + k_len = k->qe + k->ts; + } + + assert(i_len >= k_len); + return i_len - k_len; +} + +int64_t normlize_gdis_exact(ma_ug_t *ug, uc_block_t *a, uint32_t i, uint32_t k) +{ + assert(i > k); + uint32_t li, lk, pk; uint64_t l; + for (li = i, l = 0; i != (uint32_t)-1 && i >= k; i = a[i].pidx) { + li = i; l += a[i].pdis; + } + i = li; + if(i == k) return l; + assert(i > k); pk = k; + for (lk = k; k != (uint32_t)-1 && k <= i; k = a[k].aidx) lk = k; + for (k = lk; k != (uint32_t)-1 && k != pk; k = a[k].pidx) l += a[k].pdis; + assert(k == pk); + k = lk; + assert(i > k); + return l + normlize_gdis(ug, &(a[i]), &(a[k])); +} + +uint32_t dis_check_integer_aln_t(all_ul_t *ul_idx, ul_str_idx_t *str_idx, ma_ug_t *ug, integer_aln_t *li, integer_aln_t *lk, +uint32_t qid, uint32_t tid, uint32_t is_rev, double diff_rate) +{ + assert(((uint32_t)li->tn_rev_qk) >= ((uint32_t)lk->tn_rev_qk)); + if(((uint32_t)li->tn_rev_qk) == ((uint32_t)lk->tn_rev_qk)) return 0; + if(li->tk <= lk->tk) return 0; + ul_str_t *q = &(str_idx->str.a[qid]), *t = &(str_idx->str.a[tid]); + uc_block_t *iq, *it, *kq, *kt; int64_t qlen, tlen, mm, dd, i_qk, k_qk, i_tk, k_tk; + i_qk = ((uint32_t)li->tn_rev_qk); k_qk = ((uint32_t)lk->tn_rev_qk); + if(!is_rev) { + i_tk = li->tk; k_tk = lk->tk; + } else { + i_tk = t->cn - lk->tk - 1; k_tk = t->cn - li->tk - 1; + } + iq = &(ul_idx->a[qid].bb.a[(q->a[i_qk]>>32)]); + it = &(ul_idx->a[tid].bb.a[(t->a[i_tk]>>32)]); + + kq = &(ul_idx->a[qid].bb.a[(q->a[k_qk]>>32)]); + kt = &(ul_idx->a[tid].bb.a[(t->a[k_tk]>>32)]); + qlen = normlize_gdis(ug, iq, kq); tlen = normlize_gdis(ug, it, kt); + if(qlen < tlen) { + mm = qlen; dd = tlen - qlen; + } else { + mm = tlen; dd = qlen - tlen; + } + + if(dd <= (mm*diff_rate)) return 1; + qlen = normlize_gdis_exact(ug, ul_idx->a[qid].bb.a, (q->a[i_qk]>>32), (q->a[k_qk]>>32)); + tlen = normlize_gdis_exact(ug, ul_idx->a[tid].bb.a, (t->a[i_tk]>>32), (t->a[k_tk]>>32)); + if(qlen < tlen) { + mm = qlen; dd = tlen - qlen; + } else { + mm = tlen; dd = qlen - tlen; + } + + if(dd <= (mm*diff_rate)) return 1; + + return 0; +} + +uint32_t integer_chain(uint32_t qid, integer_aln_t *a, int64_t a_n, int64_t offset, integer_t *buf, +ma_ug_t *ug, ul_str_idx_t *str_idx, all_ul_t *ul_idx, ul_chain_t *res) +{ + res->v = res->s = res->e = (uint32_t)-1; res->sc = (uint64_t)-1; + res->q_sidx = res->q_eidx = res->t_sidx = res->t_eidx = (uint32_t)-1; + if(a_n <= 0) return 0; + ///already sorted by qe + int64_t i, k, max_f, max_k, sc, csc, *p, *f, tf, ti, is_circle; + integer_aln_t *li, *lk; uint32_t tid = a[0].tn_rev_qk>>33; uint32_t is_rev = (a[0].tn_rev_qk>>32)&1; + for (i = 1, sc = 0; i < a_n; ++i) { + sc += ug->u.a[a[i].vq>>1].n; + if(a[i].tk <= a[i-1].tk) break;///== means there is a circle + if(((uint32_t)a[i].tn_rev_qk) <= ((uint32_t)a[i-1].tn_rev_qk)) break; + } + if(i >= a_n) { + sc += ug->u.a[a[0].vq>>1].n; + res->v = a[0].tn_rev_qk>>32; res->s = offset; res->e = offset + a_n; res->sc = sc; + return 1; + } + is_circle = 1; + if(str_idx->str.a[qid].is_cir == 0 && str_idx->str.a[tid].is_cir == 0) is_circle = 0; + + buf->p.n = buf->f.n = 0; + kv_resize(int64_t, buf->p, (uint64_t)a_n); p = buf->p.a; + kv_resize(int64_t, buf->f, (uint64_t)a_n); f = buf->f.a; + + tf = ti = -1; + for (i = 0; i < a_n; ++i) { + li = &(a[i]); csc = ug->u.a[li->vq>>1].n; + max_f = csc; max_k = -1; + for (k = i-1; k >= 0; --k) { + lk = &(a[k]); + if(lk->tk >= li->tk) continue; + if(is_circle && (!dis_check_integer_aln_t(ul_idx, str_idx, ug, li, lk, qid, tid, is_rev, 0.04))) continue; + sc = csc + f[k]; + if(sc > max_f) { + max_f = sc; max_k = k; + } + } + f[i] = max_f; p[i] = max_k; + if(tf < max_f) { + tf = max_f; ti = i; + } + } + + if(ti < 0) return 0; + for (i = ti, k = 0; i >= 0; i = p[i]) f[k++] = i; + for (i = sc = 0; k >= 0; k--) { + a[i] = a[f[k]]; sc += ug->u.a[a[i].vq>>1].n; i++; + } + res->v = a[0].tn_rev_qk>>32; res->s = offset; res->e = offset + i; res->sc = sc; + return 1; +} + +/** +void calculate_boundary_integer_length(all_ul_t *ul_idx, ul_str_t *str, int64_t qid, int64_t tid, int64_t qk, int64_t tk, +int64_t is_rev, int64_t is_prefix, int64_t is_suffix, int64_t *r_qoff, int64_t *r_toff) +{ + (*r_qoff) = (*r_toff) = -1; + int64_t qlen = ul_idx->a[qid].rlen, tlen = ul_idx->a[tid].rlen, qoff, toff, q_ext, t_ext; + ul_str_t *qstr = &(str[qid]), *tstr = &(str[tid]); + if(is_rev) tk = tstr->cn - tk; + uc_block_t *q_b = &(ul_idx->a[qid].bb.a[qstr->a[qk]>>32]); + uc_block_t *t_b = &(ul_idx->a[tid].bb.a[tstr->a[tk]>>32]); + if(is_prefix) { + qoff = q_b->qs; toff = (is_rev?(tlen-t_b->qe):(t_b->qs)); + } + if(is_suffix) { + qoff = qlen - q_b->qe; toff = (is_rev?(t_b->qs):(tlen-t_b->qe)); + } + + if(qoff <= toff) { + t_ext = qoff; q_ext = -1; + } else { + q_ext = toff; t_ext = -1; + } + + if(q_ext == -1) { + if(is_prefix) (*r_qoff) = 0; + if(is_suffix) (*r_qoff) = (int64_t)(qstr->cn); + } + +} + +void push_ul_snp_t(int64_t chain_id, ul_str_t *str, integer_t *buf, integer_aln_t *a, int64_t a_n, int64_t qid, int64_t tid, int64_t is_rev) +{ + if(a_n <= 0) return; + int64_t k, qk, tk, p_qk, p_tk; ul_snp_t *p; + ul_str_t *qstr = &(str[qid]), *tstr = &(str[tid]); + for (k = 0, p_qk = 0, p_tk = 0; k < a_n; k++) { + qk = (uint32_t)a[k].tn_rev_qk; tk = a[k].tk; + if(qk - p_qk > 0 || tk - p_tk > 0) { + if(p_qk == 0 && p_tk == 0) {///first window + ///qk == 0 || tk == 0 means we already reach the end + if(qk > 0 && tk > 0) { + ; + } + } else { + kv_pushp(ul_snp_t, buf->snp, &p); + p->chain_id = chain_id; p->is_rev = is_rev; + p->qidx_occ = p_qk; p->qidx_occ += (qk - p_qk); + p->tidx_occ = p_tk; p->tidx_occ += (tk - p_tk); + } + } + p_qk = qk + 1; p_tk = tk + 1; + } + + qk = qstr->cn; tk = tstr->cn; + if(qk - p_qk > 0 || tk - p_tk > 0) { + kv_pushp(ul_snp_t, buf->snp, &p); + p->chain_id = chain_id; p->is_rev = is_rev; + p->qidx_occ = p_qk; p->qidx_occ += (qk - p_qk); + p->tidx_occ = p_tk; p->tidx_occ += (tk - p_tk); + } +} + +void integer_variant_call(integer_t *buf, ul_snp_t *a, int64_t a_n) +{ + kv_resize(int64_t, buf->f, (uint64_t)a_n); + int64_t *f = buf->f.a; int64_t k, m; ul_snp_t *p; + memset(f, -1, sizeof((*f))*a_n); + for (k = 0; k < a_n; k++) { + if(f[k] != (uint64_t)-1) continue; + for (m = k+1, p = &(a[k]); m < a_n; m++) { + + } + } + + +} + +void integer_phase(ul_str_t *str, integer_t *buf, ul_chain_t *idx, int64_t idx_n, integer_aln_t *aln, int64_t qid) +{ + int64_t k; integer_aln_t *a; buf->snp.n = 0; + for (k = 0; k < idx_n; k++) { + push_ul_snp_t(k, str, buf, aln + idx[k].s, idx[k].e - idx[k].s, qid, idx[k].v>>1, idx[k].v&1); + } + + int64_t z, snp_n = buf->snp.n; + radix_sort_ul_snp_t_srt(buf->snp.a, buf->snp.a + buf->snp.n); + for (z = 0, k = 1; k <= snp_n; k++) { + if(k == snp_n || buf->snp.a[z].qidx_occ != buf->snp.a[k].qidx_occ) { + integer_variant_call(buf, buf->snp.a + z, k - z); + z = k; + } + } +} +**/ + +int64_t append_connective(integer_aln_t *aln, ul_chain_t *idx, int64_t str_i, uint64_t occ_thres, uint64_t *res) +{ + int64_t k, kl = idx->e - idx->s, str_k = -1, fp = 0; integer_aln_t *a = aln + idx->s; + assert(idx->sc <= (uint64_t)kl); + for (k = idx->sc; k < kl; k++) { + str_k = (uint32_t)a[k].tn_rev_qk; + if(str_k >= str_i) break; + } + idx->sc = k; + if(k >= kl || str_k != str_i) return 0; + + for (k -= 1; k >= 0; k--) { + str_k = (uint32_t)a[k].tn_rev_qk; + res[str_k]++; + if(res[str_k] == occ_thres) fp++; + } + return fp; +} +///occ_thres does not consider reference read itself; so the real coverage is (occ_thres+1) +int64_t integer_chain_dp(bubble_type *bub, integer_t *buf, ul_str_t *str, integer_aln_t *aln, ul_chain_t *idx, int64_t idx_n, ma_ug_t *ug, int64_t qid, uint64_t occ_thres) +{ + ul_str_t *qstr = &(str[qid]); int64_t k, q_n = qstr->cn, *f, *p, z, max_f, tf, tk, max_p, done_z, sc, csc; + kv_resize(int64_t, buf->f, qstr->cn); kv_resize(int64_t, buf->p, qstr->cn); + kv_resize(uint64_t, buf->o, qstr->cn); uint64_t *o; + f = buf->f.a; p = buf->p.a; o = buf->o.a; + + // radix_sort_ul_chain_t_srt(idx, idx + idx_n); + for (k = 0; k < idx_n; k++) idx[k].sc = 0; + + for (k = 0, tf = tk = -1; k < q_n; k++) { + csc = ug->u.a[((uint32_t)qstr->a[k])>>1].n; + max_p = -1; max_f = csc; + // if(!IF_HOM((((uint32_t)qstr->a[k])>>1), (*bub))) { + if(k > 0) { + memset(o, 0, sizeof((*o))*k); + for (z = done_z = 0; z < idx_n && done_z < k; z++) { + done_z += append_connective(aln, &(idx[z]), k, occ_thres, o); + } + assert(done_z <= k); + for (z = k - 1; z >= 0; z--) { + if(o[z] < occ_thres) continue; + sc = csc + f[z]; + if(sc > max_f) { + max_f = sc; max_p = z; + } + } + } + f[k] = max_f; p[k] = max_p; + if(tf < max_f) { + tf = max_f; tk = k; + } + + // fprintf(stderr, "[M::%s::k->%ld] f[k]->%ld, p[k]->%ld\n", __func__, k, f[k], p[k]); + } + + for (k = tk, done_z = 0; k >= 0; k = p[k]) done_z++; + for (k = tk, sc = done_z; k >= 0; k = p[k]) o[--done_z] = k; + return sc; +} + +void print_integer_seq(ma_ug_t *ug, ul_str_t *str, int64_t id, int64_t is_header) { uint64_t i; + if(is_header) { + fprintf(stderr,"[M::%s::tid->%ld] occ::%u\n", __func__, id, str[id].cn); + } + for (i = 0; i < str[id].cn; i++) { + fprintf(stderr, "utg%.6d%c(%c)\t", (((uint32_t)str[id].a[i])>>1)+1, + "lc"[ug->u.a[(((uint32_t)str[id].a[i])>>1)].circ], "+-"[(((uint32_t)str[id].a[i])&1)]); + } + fprintf(stderr,"\n"); +} + +int64_t utg_cover_read_occ_by_qs(ma_ug_t *ug, int64_t oqs, int64_t oqe, uc_block_t *x) +{ + assert(oqs >= (int64_t)x->qs && oqs <= (int64_t)x->qe); + assert(oqe >= (int64_t)x->qs && oqe <= (int64_t)x->qe); + assert(oqe > oqs); + int64_t s_off, e_off, k; + s_off = get_offset_adjust(oqs-x->qs, x->qe-x->qs, x->te-x->ts); + e_off = get_offset_adjust(x->qe-oqe, x->qe-x->qs, x->te-x->ts); + if(x->rev) { + k = s_off; s_off = e_off; e_off = k; + } + + return ug_occ_w(x->ts + s_off, x->te - e_off, &(ug->u.a[x->hid])); +} + +int64_t estimate_ul_len(ma_ug_t *ug, ul_vec_t *raw_ov, ul_str_t *str, int64_t idx, int64_t ext_len, uint64_t *cov_buf) +{ + if(ext_len == 0) return idx; + int64_t aim_s, aim_e, k = idx, qs, qe, ovlp_s, ovlp_e, cov_occ, str_n = str->cn; + if(ext_len < 0) { + aim_s = ((int64_t)(raw_ov->bb.a[str->a[idx]>>32].qs)) + ext_len; + if(aim_s < 0) aim_s = 0; + aim_e = raw_ov->bb.a[str->a[idx]>>32].qe; + for (k = idx - 1; k >= 0; k--) { + qs = raw_ov->bb.a[str->a[k]>>32].qs; qe = raw_ov->bb.a[str->a[k]>>32].qe; + // if(qe <= aim_s) break; + ovlp_s = MAX(aim_s, qs); ovlp_e = MIN(aim_e, qe); + if(ovlp_e <= ovlp_s) { + k++; + break; + } + cov_occ = utg_cover_read_occ_by_qs(ug, ovlp_s, ovlp_e, &(raw_ov->bb.a[str->a[k]>>32])); + if(cov_occ == 0) { + k++; + break; + } + cov_buf[k] = cov_occ; + } + if(k < 0) k++; + } else { + aim_s = raw_ov->bb.a[str->a[idx]>>32].qs; + aim_e = ((int64_t)(raw_ov->bb.a[str->a[idx]>>32].qe)) + ext_len; + if(aim_e > (int64_t)(raw_ov->rlen)) aim_e = raw_ov->rlen; + for (k = idx + 1; k < str_n; k++) { + qs = raw_ov->bb.a[str->a[k]>>32].qs; qe = raw_ov->bb.a[str->a[k]>>32].qe; + // if(qs >= aim_e) + ovlp_s = MAX(aim_s, qs); ovlp_e = MIN(aim_e, qe); + if(ovlp_e <= ovlp_s) { + k--; + break; + } + cov_occ = utg_cover_read_occ_by_qs(ug, ovlp_s, ovlp_e, &(raw_ov->bb.a[str->a[k]>>32])); + if(cov_occ == 0) { + k--; + break; + } + cov_buf[k] = cov_occ; + } + if(k >= str_n) k--; + } + return k; +} + +void print_integer_ovlps(ma_ug_t *ug, ul_str_t *str, integer_aln_t *aln, int64_t aln_occ, ul_chain_t *idx, int64_t idx_n, int64_t qid, int64_t consenus_occ) +{ + fprintf(stderr, "\n[M::%s::qid->%ld] qstr->cn::%u, aln_occ::%ld, idx_n::%ld, consenus_occ::%ld\n", + __func__, qid, str[qid].cn, aln_occ, idx_n, consenus_occ); + print_integer_seq(ug, str, qid, 0); + int64_t k, z, z_n, qk, tk, is_rev, tid; //uint64_t z; + for (k = 0; k < idx_n; k++) { + fprintf(stderr, "\n[M::%s::tid->%lu] rev->%lu, aln_n->%u\n", __func__, aln[idx[k].s].tn_rev_qk>>33, (aln[idx[k].s].tn_rev_qk>>32)&1, idx[k].e - idx[k].s); + print_integer_seq(ug, str, aln[idx[k].s].tn_rev_qk>>33, 0); + z = idx[k].s; z_n = idx[k].e; tid = aln[idx[k].s].tn_rev_qk>>33; + for (; z < z_n; z++) { + qk = (uint32_t)aln[z].tn_rev_qk; is_rev = ((aln[z].tn_rev_qk>>32)&1); + tk = ((is_rev == 0)? (aln[z].tk):(str[tid].cn - aln[z].tk - 1)); + fprintf(stderr, "[qk::%ld]utg%.6d%c(%c) <---> [tk::%ld]utg%.6d%c(%c)\n", + qk, (((uint32_t)str[qid].a[qk])>>1)+1, "lc"[ug->u.a[(((uint32_t)str[qid].a[qk])>>1)].circ], "+-"[(((uint32_t)str[qid].a[qk])&1)], + tk, (((uint32_t)str[tid].a[tk])>>1)+1, "lc"[ug->u.a[(((uint32_t)str[tid].a[tk])>>1)].circ], "+-"[(((uint32_t)str[tid].a[tk])&1)]); + } + + // for (z = idx[k].s; z < idx[k].e; z++) { + // aln[z].tn_rev_qk + // } + } +} + +int64_t integer_align_extention(all_ul_t *ul_idx, int64_t qid, ul_str_t *q_str, int64_t tid, ul_str_t *t_str, integer_aln_t *idx) +{ + ; +} + +int64_t refine_integer_ovlps(all_ul_t *ul_idx, bubble_type *bub, ma_ug_t *ug, ul_str_t *str, integer_aln_t *aln, ul_chain_t *idx, int64_t qid, integer_t *buf, +uint64_t *cns, uint64_t cns_occ) +{ + if(idx->e<=idx->s) return 0; + int64_t qk, tk, is_rev, tid, qoff, toff, qlen, tlen, ext, q_end, t_end, z; + ul_str_t *q_str, *t_str; uc_block_t *q_b, *t_b; integer_aln_t *x, *y; + tid = aln[idx->s].tn_rev_qk>>33; is_rev = ((aln[idx->s].tn_rev_qk>>32)&1); + q_str = &(str[qid]); t_str = &(str[tid]); qlen = ul_idx->a[qid].rlen; tlen = ul_idx->a[tid].rlen; + kv_resize(uint64_t, buf->u, buf->u.n + q_str->cn + t_str->cn); + uint64_t *q_cov_buf = buf->u.a + buf->u.n, *t_cov_buf = buf->u.a + buf->u.n + q_str->cn, k, rg_occ, cn_k; + memset(q_cov_buf, -1, sizeof((*q_cov_buf))*q_str->cn); memset(t_cov_buf, -1, sizeof((*t_cov_buf))*t_str->cn); + + //beg + x = &(aln[idx->s]); + qk = (uint32_t)(x->tn_rev_qk); tk = ((is_rev == 0)? (x->tk):(t_str->cn - x->tk - 1)); + + ///direction + if((qk > 0) && (x->tk > 0)) { + q_b = &(ul_idx->a[qid].bb.a[q_str->a[qk]>>32]); t_b = &(ul_idx->a[tid].bb.a[t_str->a[tk]>>32]); + ///direction + qoff = q_b->qs; toff = (is_rev?(tlen-t_b->qe):(t_b->qs)); + if (qoff <= toff) ext = qoff; + else ext = toff; + + q_end = estimate_ul_len(ug, &(ul_idx->a[qid]), q_str, qk, -ext, q_cov_buf); + t_end = estimate_ul_len(ug, &(ul_idx->a[tid]), t_str, tk, ((is_rev)?(ext):(-ext)), t_cov_buf); + // if(tid == 316) { + // fprintf(stderr, "[M::%s::tid->%ld] qoff::%ld, toff::%ld, ext::%ld, q_end::%ld, t_end::%ld\n", + // __func__, tid, qoff, toff, ext, q_end, t_end); + // } + idx->q_sidx = q_end; idx->t_sidx = ((is_rev == 0)? (t_end):(t_str->cn - t_end)); + } else { + idx->q_sidx = (uint32_t)(x->tn_rev_qk); idx->t_sidx = x->tk; + } + assert((idx->q_sidx <= ((uint32_t)(x->tn_rev_qk))) && (idx->t_sidx <= x->tk)); + + //end + x = &(aln[idx->e-1]); + qk = (uint32_t)(x->tn_rev_qk); tk = ((is_rev == 0)? (x->tk):(t_str->cn - x->tk - 1)); + // if(tid == 316) { + // fprintf(stderr, "[end-M::%s::tid->%ld] qk::%ld, tk::%ld, x->tk::%u, t_str->cn::%u\n", + // __func__, tid, qk, tk, x->tk, t_str->cn); + // } + ///direction + if((((uint32_t)qk + 1) < q_str->cn) && ((x->tk + 1) < t_str->cn)) { + q_b = &(ul_idx->a[qid].bb.a[q_str->a[qk]>>32]); t_b = &(ul_idx->a[tid].bb.a[t_str->a[tk]>>32]); + ///direction + qoff = qlen - q_b->qe; toff = (is_rev?(t_b->qs):(tlen-t_b->qe)); + if (qoff <= toff) ext = qoff; + else ext = toff; + q_end = estimate_ul_len(ug, &(ul_idx->a[qid]), q_str, qk, ext, q_cov_buf); + t_end = estimate_ul_len(ug, &(ul_idx->a[tid]), t_str, tk, ((is_rev)?(-ext):(ext)), t_cov_buf); + // if(tid == 316) { + // fprintf(stderr, "[end-M::%s::tid->%ld] qk::%ld, tk::%ld, qoff::%ld, toff::%ld, q_end::%ld, t_end::%ld, ext::%ld\n", + // __func__, tid, qk, tk, qoff, toff, q_end, t_end, ext); + // } + idx->q_eidx = q_end + 1; idx->t_eidx = ((is_rev == 0)? (t_end + 1):(t_str->cn - t_end)); + } else { + idx->q_eidx = (uint32_t)(x->tn_rev_qk)+1; idx->t_eidx = x->tk+1; + } + assert(idx->q_eidx > ((uint32_t)(x->tn_rev_qk)) && (idx->t_eidx > x->tk)); + assert((idx->q_eidx > idx->q_sidx) && (idx->t_eidx > idx->t_sidx)); + // if(tid == 316) { + // fprintf(stderr, "[M::%s::tid->%ld] idx->q_sidx::%u, idx->q_eidx::%u, idx->t_sidx::%u, idx->t_eidx::%u\n", + // __func__, tid, idx->q_sidx, idx->q_eidx, idx->t_sidx, idx->t_eidx); + // } + + ///mid + for (k = idx->s; k < idx->e; k++) {///go through all alignment pairs + x = &(aln[k]); y = ((k > idx->s)? (&(aln[k-1])):(NULL)); + qk = (uint32_t)(x->tn_rev_qk); + q_cov_buf[qk] = buf->u.a[qk]; + + tk = ((is_rev == 0)? (x->tk):(t_str->cn - x->tk - 1)); + t_b = &(ul_idx->a[tid].bb.a[t_str->a[tk]>>32]); + assert(((t_b->hid<<1)+t_b->rev)==((uint32_t)t_str->a[tk])); + t_cov_buf[tk] = ug_occ_w(t_b->ts, t_b->te, &(ug->u.a[t_b->hid])); + + q_cov_buf[qk] += ((uint64_t)(0x8000000000000000)); + t_cov_buf[tk] += ((uint64_t)(0x8000000000000000)); + + + if(!y) continue; + for (z = ((uint32_t)(y->tn_rev_qk)) + 1; z < qk; z++) { + q_cov_buf[z] = buf->u.a[z]; + } + + z = ((is_rev == 0)? (y->tk):(t_str->cn - x->tk - 1)) + 1; + tk = ((is_rev == 0)? (x->tk):(t_str->cn - y->tk - 1)); + assert(z <= tk); + for (; z < tk; z++) { + t_b = &(ul_idx->a[tid].bb.a[t_str->a[z]>>32]); + assert(((t_b->hid<<1)+t_b->rev)==((uint32_t)t_str->a[z])); + t_cov_buf[z] = ug_occ_w(t_b->ts, t_b->te, &(ug->u.a[t_b->hid])); + } + } + + if(cns) { + int64_t cns_cov_occ = 0, cns_het_occ = 0, match_cns_occ = 0, match_cns_het_occ = 0, hm; + rg_occ = idx->q_eidx; + for (k = idx->q_sidx, cn_k = 0; k < rg_occ; k++) { + assert((q_cov_buf[k]!=(uint64_t)-1) && (q_cov_buf[k] > 0)); + hm = ((q_cov_buf[k]<<1)>>1); + for (; cn_k < cns_occ; cn_k++) { + if(cns[cn_k] == k) { + cns_cov_occ += hm; + if(q_cov_buf[k]&((uint64_t)(0x8000000000000000))) match_cns_occ += hm; + if(!IF_HOM((((uint32_t)q_str->a[k])>>1), (*bub))) { + cns_het_occ += hm; + if(q_cov_buf[k]&((uint64_t)(0x8000000000000000))) match_cns_het_occ += hm; + } + } else if(cns[cn_k] > k) { + break; + } + } + } + + if(cns_cov_occ <= 0) return 0; + if((cns_het_occ > 0) && (match_cns_het_occ <= (cns_het_occ*0.5))) return 0; + if(match_cns_occ <= (cns_cov_occ*0.5)) return 0; + + k = idx->t_sidx; rg_occ = idx->t_eidx; + if(is_rev) { + k = t_str->cn - idx->t_eidx; rg_occ = t_str->cn - idx->t_sidx; + } + cns_cov_occ = match_cns_occ = 0; + assert(k < rg_occ); + for (; k < rg_occ; k++) { + assert((t_cov_buf[k]!=(uint64_t)-1) && (t_cov_buf[k] > 0)); + hm = ((t_cov_buf[k]<<1)>>1); + cns_cov_occ += hm; + if(t_cov_buf[k]&((uint64_t)(0x8000000000000000))) match_cns_occ += hm; + } + + if(cns_cov_occ <= 0 || match_cns_occ <= 0) return 0; + if(match_cns_occ <= (cns_cov_occ*0.5)) return 0; + } + + // print_integer_ovlps(ug, str, buf->b.a, buf->b.n, idx, 1, qid, buf->o.n); + // fprintf(stderr, "[M::%s::tid->%ld] idx->q_sidx::%u, idx->q_eidx::%u, idx->t_sidx::%u, idx->t_eidx::%u\n******************************************************\n", + // __func__, tid, idx->q_sidx, idx->q_eidx, idx->t_sidx, idx->t_eidx); + + return 1; +} + + + +void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid) +{ + if(qid != 267) return; + uint64_t k, z, m_het, m_het_occ, ref_occ, b_n, m; uint32_t vk, vz; integer_aln_t *p; ul_chain_t sc; + ul_str_idx_t *str_idx = &(uidx->pstr); ma_ug_t *ug = uidx->l1_ug; + uint64_t *hid_a, hid_n; uc_block_t *xi; + ul_str_t *str = &(str_idx->str.a[qid]); + if(str->cn < 2) return; ///directly filter out too short UL + kv_resize(uint64_t, buf->u, str->cn); buf->u.n = str->cn; + for (k = m_het = m_het_occ = ref_occ = 0; k < str->cn; k++) { + xi = &(uidx->idx->a[qid].bb.a[str->a[k]>>32]); + assert(((xi->hid<<1)+xi->rev)==((uint32_t)str->a[k])); + buf->u.a[k] = ug_occ_w(xi->ts, xi->te, &(ug->u.a[xi->hid])); + if(!IF_HOM((((uint32_t)str->a[k])>>1), (*uidx->bub))) { + m_het++; m_het_occ += buf->u.a[k]; + } + ref_occ += buf->u.a[k]; + } + if(m_het < 2 && m_het > 0) return;///if all matched unitigs are hom, is ok + + for (k = 0, buf->b.n = 0; k < str->cn; k++) { + vk = (uint32_t)str->a[k]; + hid_a = str_idx->occ.a + str_idx->idx.a[vk>>1]; + hid_n = str_idx->idx.a[(vk>>1)+1] - str_idx->idx.a[vk>>1]; + for (z = 0; z < hid_n; z++) { + if((hid_a[z]>>32) == qid) continue; + if(str_idx->str.a[hid_a[z]>>32].cn < 2) continue; + vz = (uint32_t)(str_idx->str.a[hid_a[z]>>32].a[(uint32_t)hid_a[z]]); + assert((vk>>1) == (vz>>1)); + kv_pushp(integer_aln_t, buf->b, &p); + p->vq = vk; p->tk = (uint32_t)hid_a[z]; + if((vk^vz)&1) p->tk = str_idx->str.a[hid_a[z]>>32].cn - p->tk - 1;///rev + p->tn_rev_qk = (hid_a[z]>>32); p->tn_rev_qk <<= 1; p->tn_rev_qk |= ((vk^vz)&1); + p->tn_rev_qk <<= 32; p->tn_rev_qk += k; + } + } + + radix_sort_integer_aln_t_srt(buf->b.a, buf->b.a + buf->b.n); + b_n = buf->b.n; buf->sc.n = 0; + for (k = 1, z = 0; k < b_n; k++) { + if(k == b_n || (buf->b.a[z].tn_rev_qk>>32) != (buf->b.a[k].tn_rev_qk>>32)) { + if(integer_chain(qid, buf->b.a + z, k - z, z, buf, ug, str_idx, uidx->idx, &sc) && sc.v != (uint32_t)-1) { + if((buf->sc.n > 0) && ((buf->sc.a[buf->sc.n-1].v>>1) == (sc.v>>1)) + && (buf->sc.a[buf->sc.n-1].sc < sc.sc)) { + buf->sc.n--; + } + kv_push(ul_chain_t, buf->sc, sc); + } + z = k; + } + } + + uint64_t *o, o_n, cns_het, cns_het_occ, ref_cns_occ; + o_n = integer_chain_dp(uidx->bub, buf, str_idx->str.a, buf->b.a, buf->sc.a, buf->sc.n, ug, qid, 2); + assert(o_n <= str->cn); + if(o_n <= 0 || o_n == str->cn) return; + + for (k = cns_het = cns_het_occ = ref_cns_occ = 0, o = buf->o.a; k < o_n; k++) { + if(!IF_HOM((((uint32_t)str->a[o[k]])>>1), (*uidx->bub))) { + cns_het++; cns_het_occ += buf->u.a[o[k]]; + } + ref_cns_occ += buf->u.a[o[k]]; + } + buf->o.n = o_n; + ///1. if the ref read only has hom unitigs, is fine + ///2. otherwise need to have consenus het untigs + if(m_het > 0 && cns_het <= 0) return; + if(m_het_occ > 0 && cns_het_occ <= (m_het_occ*0.25)) return; + if(ref_cns_occ <= (ref_occ*0.5)) return; + + for (k = m = 0; k < buf->sc.n; k++) { + // fprintf(stderr, "[M::%s::k->%lu] m::%lu\n", __func__, k, m); + // if(k == 72) { + // print_integer_ovlps(ug, str_idx->str.a, buf->b.a, buf->b.n, buf->sc.a+k, 1, qid, o_n); + // } + if(refine_integer_ovlps(uidx->idx, uidx->bub, ug, str_idx->str.a, buf->b.a, + &(buf->sc.a[k]), qid, buf, o, o_n)) { + buf->sc.a[m++] = buf->sc.a[k]; + } else { + print_integer_ovlps(ug, str_idx->str.a, buf->b.a, buf->b.n, buf->sc.a+k, 1, qid, o_n); + fprintf(stderr, "[M::%s::] idx->q_sidx::%u, idx->q_eidx::%u, idx->t_sidx::%u, idx->t_eidx::%u\n******************************************************\n", + __func__, buf->sc.a[k].q_sidx, buf->sc.a[k].q_eidx, buf->sc.a[k].t_sidx, buf->sc.a[k].t_eidx); + } + } + buf->sc.n = m; + // if(m != str->cn) print_integer_ovlps(uidx->l1_ug, str_idx->str.a, buf->b.a, buf->b.n, buf->sc.a, buf->sc.n, qid, m); + + + // radix_sort_ul_chain_t_srt(buf->sc.a, buf->sc.a + buf->sc.n); + + // integer_phase(str_idx->str.a, buf, buf->sc.a, buf->sc.n, buf->b.a, qid); + + // radix_sort_ul_chain_t_srt(buf->sc.a, buf->sc.a + buf->sc.n); +} + +static void worker_integer_correction(void *data, long i, int tid) // callback for kt_for() +{ + ul_resolve_t *uidx = (ul_resolve_t *)data; + integer_ml_t *sl = &(uidx->str_b); + integer_t *buf = &(sl->buf[tid]); + // uc_block_t *uls; uint64_t uls_n; uint32_t v; + // uint64_t *srt_a, srt_n, is_circle = ((uidx->psrt.idx.a[i]&((uint64_t)(0x100000000)))?1:0); + // srt_a = uidx->psrt.srt.a + (uint32_t)uidx->psrt.idx.a[i]; srt_n = uidx->psrt.idx.a[i]>>33; + // if(srt_n == 0) return; + // integer_candidate(uidx, srt_a, srt_n, is_circle, buf); + integer_candidate(uidx, buf, i); +} + + +void integer_correction(ul_resolve_t *uidx) +{ + integer_ml_t sl; + init_integer_ml_t(&sl, uidx, asm_opt.thread_num); +} + + +void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg) +{ + uint64_t i; uint8_t *r_het = NULL; bubble_type *bub = NULL; for (i = 0; i < sg->n_seq; ++i) { if(sg->seq[i].del) continue; sg->seq[i].c = PRIMARY_LABLE; } hic_clean(sg); - ul_realignment(uopt, sg); - // if(ul_refine_alignment(uopt, sg)) update_sg_uo(sg, src); + ma_ug_t *init_ug = ul_realignment(uopt, sg); + filter_sg_by_ug(sg, init_ug, uopt); + + bub = gen_bubble_chain(sg, init_ug, uopt, &r_het); + + ul_resolve_t *uidx = init_ul_resolve_t(sg, init_ug, bub, &UL_INF, r_het); + kt_for(asm_opt.thread_num, worker_integer_correction, uidx, uidx->idx->n); + + print_debug_ul("UL.debug", init_ug, sg, uopt->coverage_cut, uopt->sources, uopt->ruIndex, bub, &UL_INF); + + resolve_dip_bub_chains(uidx); + + // free(r_het); destory_bubbles(bub); free(bub); } \ No newline at end of file diff --git a/hic.cpp b/hic.cpp index 37b0ff3..cc06ae0 100644 --- a/hic.cpp +++ b/hic.cpp @@ -6416,15 +6416,15 @@ asg_arc_t *p, uint32_t check_het) if(x_1_b_id) id1 = (*x_1_b_id); if((x_0 != (uint32_t)-1) && (x_1 != (uint32_t)-1)) { - if(((x_0>>1) == (x_1>>1))) + if(((x_0>>1) == (x_1>>1)))///if we would like to find a edge bridging two nearby bubbles { if(x_0_b_id == NULL && x_1_b_id == NULL) { get_bub_id(bub, x_0>>1, &id0, &id1, check_het); } - get_bubbles(bub, id0, &beg_0, &sink_0, &a, &n, NULL); - get_bubbles(bub, id1, &beg_1, &sink_1, &a, &n, NULL); + get_bubbles(bub, id0, &beg_0, &sink_0, &a, &n, NULL);//first bubble + get_bubbles(bub, id1, &beg_1, &sink_1, &a, &n, NULL);//second bubble ori_0 = (uint64_t)-1; @@ -6716,10 +6716,11 @@ void detect_bub_graph(bubble_type* bub, asg_t *untig_sg) uint64_t pLen, rLEN, r_hetLen; ma_utg_t *u = NULL; asg_arc_t *t = NULL; - for (i = 0; i < ug->u.n; i++) + for (i = 0; i < ug->u.n; i++)///bubble chain graph; bubbles within the same chain have been merged { u = &(ug->u.a[i]); if(u->n == 0) continue; + ///u is a bubble chain for (k = pLen = rLEN = r_hetLen = beg_idx = 0, end_idx = -1; k < u->n; k++) { rId = u->a[k]>>33; @@ -6727,7 +6728,7 @@ void detect_bub_graph(bubble_type* bub, asg_t *untig_sg) get_bubbles(bub, rId, ori == 1?&root:&r_root, ori == 0?&root:&r_root, NULL, NULL, NULL); t = NULL; - if(k+1 < u->n) t = &(arc_first(bg, u->a[k]>>32)); + if(k+1 < u->n) t = &(arc_first(bg, u->a[k]>>32));//edge between two bubbles within the same chain pLen += bg->seq[rId].len;///path length in bubble if(end_idx < beg_idx) ///first bubble @@ -6739,7 +6740,7 @@ void detect_bub_graph(bubble_type* bub, asg_t *untig_sg) if(t) { - if(t->el == 0) + if(t->el == 0)///there is tangle between two bubbles { if(end_idx >= beg_idx) { @@ -6752,7 +6753,7 @@ void detect_bub_graph(bubble_type* bub, asg_t *untig_sg) pLen = rLEN = r_hetLen = 0; beg_idx = k + 1; end_idx = k; } - else + else///two bubbles directly connected with each other { pLen += t->ol; rLEN += t->ol; @@ -6784,7 +6785,7 @@ void get_bub_graph(ma_ug_t* ug, bubble_type* bub) uint32_t *pre = NULL; MALLOC(pre, n_vtx); uint32_t pre_id, adjecent, bub_occ; asg_t *bub_g = asg_init(); - for (v = 0; v < bub->f_bub; v++) + for (v = 0; v < bub->f_bub; v++)///all bubbles { uint64_t pathbase; uint32_t beg, sink; @@ -6793,11 +6794,12 @@ void get_bub_graph(ma_ug_t* ug, bubble_type* bub) bub_g->seq[v].c = PRIMARY_LABLE; } - //check all unitigs + //check all unitigs, instead of bubble nodes for (v = 0; v < n_vtx; ++v) { if(sg->seq[v>>1].del) continue; if(bub->b_s_idx.a[v>>1] == (uint64_t)-1) continue; ///if (v>>1) is not a beg or sink of bubbles + ///one node might be the beg/sink node of at most two bubbles bub_occ = connect_bub_occ(bub, v>>1, bub->check_het); if(bub_occ == 0) continue; if(bub_occ == 2) @@ -6812,6 +6814,7 @@ void get_bub_graph(ma_ug_t* ug, bubble_type* bub) } if(ma_2_bub_arc(bub, v, NULL, (uint32_t)-1, NULL, &t, bub->check_het) == 0) continue; + ///v is the beg/sink node of only one bubble get_shortest_path(v, &pq, sg, pre); for (k = 0; k < pq.dis.n; k++) { @@ -6833,7 +6836,7 @@ void get_bub_graph(ma_ug_t* ug, bubble_type* bub) if(adjecent == 0) { - if(ma_2_bub_arc(bub, v, NULL, k^1, NULL, &t, bub->check_het)) + if(ma_2_bub_arc(bub, v, NULL, k^1, NULL, &t, bub->check_het))///edges spanning tangles { t.el = 0; t.ol = pq.dis.a[k] + sg->seq[k>>1].len; p = asg_arc_pushp(bub_g); @@ -6963,7 +6966,7 @@ uint8_t* vis_flag, uint32_t vis_flag_n, kvec_t_u32_warp* stack, asg_t *bsg, asg_ radix_sort_u32(broken->a.a, broken->a.a + broken->a.n); for (i = n = 0, pre = (uint32_t)-1; i < broken->a.n; i++) { - if((broken->a.a[i]>>1) == (pre>>1)) continue; + if((broken->a.a[i]>>1) == (pre>>1)) continue;///skip same node like v and v^1 pre = broken->a.a[i]; broken->a.a[n] = pre; n++; @@ -8046,7 +8049,7 @@ void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub) occ_idx.a[occ_idx.n - k - 1] = tmp; } - for (k = 0; k < bub_ug->g->n_seq; k++) + for (k = 0; k < bub_ug->g->n_seq; k++)///start from the longest chain { v = ((uint32_t)(occ_idx.a[k]))<<1; if(is_used[v] == 0 && asg_arc_n(bub_ug->g, v) > 0) @@ -8192,9 +8195,9 @@ void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint ///uint64_t end_thres; uint8_t *bsg_idx = NULL; CALLOC(bsg_idx, n_vtx>>1); - for (i = 0; i < bub_ug->u.n; i++) + for (i = 0; i < bub_ug->u.n; i++)///label all unitigs within the bubble chains { - u = &(bub_ug->u.a[i]); + u = &(bub_ug->u.a[i]);///a bubble chain if(u->n == 0) continue; for (k_i = 0; k_i < u->n; k_i++) { @@ -8213,7 +8216,7 @@ void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint new_bub = bub->b_g->n_seq; for (i = 0; i < bub_ug->u.n; i++) { - u = &(bub_ug->u.a[i]); + u = &(bub_ug->u.a[i]);///bubble chain if(u->n == 0) continue; ///end_thres = calculate_chain_weight(u, bub, ug, &x); if(is_middle) @@ -8223,7 +8226,7 @@ void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint if(k_i+1 >= u->n) continue; ///note: must igore .del here, since bsg might be changed t = &(arc_first(bsg, u->a[k_i]>>32)); - if(t->el == 1) continue; + if(t->el == 1) continue;///if a[k_i] and a[k_i+1] are directly connected without any tangle involoved rId_0 = u->a[k_i]>>33; ori_0 = u->a[k_i]>>32&1; @@ -8233,7 +8236,7 @@ void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint ori_1 = (u->a[k_i+1]>>32&1)^1; get_bubbles(bub, rId_1, ori_1 == 1?&root_1:NULL, ori_1 == 0?&root_1:NULL, NULL, NULL, NULL); - broken.a.n = 0; + broken.a.n = 0;///just collect all nodes between root_0 and root_1 get_related_bub_nodes(&broken, bub, &pq, sg, pre, root_0, root_1, NULL); get_related_bub_nodes(&broken, bub, &pq, sg, pre, root_1, root_0, NULL); ///no need to cut the edge, we still have chance to flip by chain @@ -9353,8 +9356,8 @@ void build_bub_graph(ma_ug_t* ug, bubble_type* bub) // detect_bub_graph(bub, ug->g, 1); update_bubble_chain(ug, bub, 1, 0); ///print_bubble_chain(bub, "second round"); - update_bubble_chain(ug, bub, 0, 1); - resolve_bubble_chain_tangle(ug, bub); + update_bubble_chain(ug, bub, 0, 1);///resolve tangles within bubble chains + resolve_bubble_chain_tangle(ug, bub);///resolve tangles between bubble chains } void get_forward_distance(uint32_t src, uint32_t dest, asg_t *sg, hc_links* link, MT* M) diff --git a/inter.cpp b/inter.cpp index 2efb10d..30ad30d 100644 --- a/inter.cpp +++ b/inter.cpp @@ -8465,13 +8465,14 @@ void print_ovlp_src_bl_stat(all_ul_t *x, const ug_opt_t *uopt) __func__, tt[0]+tt[1]+tt[2]+tt[3], tt[1], tt[2], tt[3]); } -void gen_ul_vec_rid_t(all_ul_t *x) +void gen_ul_vec_rid_t(all_ul_t *x, All_reads *rdb, ma_ug_t *ug) { ul_vec_rid_t *ridx = &(x->ridx); - uint64_t k, i, l, m, *a, a_n; ul_vec_t *p = NULL; - ridx->idx.n = ridx->idx.m = R_INF.total_reads + 1; CALLOC(ridx->idx.a, ridx->idx.n); + uint64_t k, i, l, m, *a, a_n, idx_n; ul_vec_t *p = NULL; + idx_n = (rdb?rdb->total_reads:ug->u.n); + ridx->idx.n = ridx->idx.m = idx_n + 1; CALLOC(ridx->idx.a, ridx->idx.n); - for (k = 0; k < x->n; k++) { + for (k = 0; k < x->n; k++) {///each UL read p = &(x->a[k]); for (i = 0; i < p->bb.n; i++) { if(p->bb.a[i].base) continue; @@ -8486,7 +8487,7 @@ void gen_ul_vec_rid_t(all_ul_t *x) } ridx->occ.n = ridx->occ.m = l; MALLOC(ridx->occ.a, ridx->occ.n); - for (k = 0; k < R_INF.total_reads; k++) { + for (k = 0; k < idx_n; k++) { a = ridx->occ.a + ridx->idx.a[k]; a_n = ridx->idx.a[k+1] - ridx->idx.a[k]; if(a_n) a[a_n-1] = 0; @@ -8504,9 +8505,105 @@ void gen_ul_vec_rid_t(all_ul_t *x) } } } - } + +uint32_t ugl_cover_check(uint64_t is, uint64_t ie, ma_utg_t *u) +{ + if(is == 0 && ie == u->len) return 1; + uint64_t l, i, us, ue; + for (i = l = 0; i < u->n; i++) { + us = l; ue = l + Get_READ_LENGTH(R_INF, (u->a[i]>>33)); + if(is <= us && ie >= ue) return 1; + if(us >= ie) break; + l += (uint32_t)u->a[i]; + } + return 0; +} + +static void update_ug_arch_ul(void *data, long i, int tid) // callback for kt_for() +{ + const ma_ug_t *ug = (ma_ug_t *)data; + asg_arc_t *e = &(ug->g->arc[i]); e->ou = 0; + uint32_t v = e->ul>>32, w = e->v, k, uv, uw; + uint64_t *a, a_n; uc_block_t *p, *n; + a = UL_INF.ridx.occ.a + UL_INF.ridx.idx.a[v>>1]; + a_n = UL_INF.ridx.idx.a[(v>>1)+1] - UL_INF.ridx.idx.a[v>>1]; + for (k = 0; k < a_n; k++) { + p = &(UL_INF.a[a[k]>>32].bb.a[(uint32_t)(a[k])]); + if(p->base || (!p->el) || (!p->pchain)) continue; + uv = (((uint32_t)(p->hid))<<1)|((uint32_t)(p->rev)); + if((uv == v) && (p->aidx != (uint32_t)-1)) { + n = &(UL_INF.a[a[k]>>32].bb.a[p->aidx]); + // if(!((!n->base)&&(n->el)&&(n->pchain)&&(n->pidx==((uint32_t)(a[k]))))) { + // fprintf(stderr, "k::%u, n->base::%u, n->el::%u, n->pchain::%u, n->pidx::%u, p->aidx::%u\n", + // k, n->base, n->el, n->pchain, n->pidx, p->aidx); + // } + assert((!n->base)&&(n->el)&&(n->pchain)&&(n->pidx==((uint32_t)(a[k])))); + uw = (((uint32_t)(n->hid))<<1)|((uint32_t)(n->rev)); + if(uw == w) e->ou++; + } + + if(((uv^1) == v) && (p->pidx != (uint32_t)-1)) { + n = &(UL_INF.a[a[k]>>32].bb.a[p->pidx]); + assert((!n->base)&&(n->el)&&(n->pchain)&&(n->aidx==((uint32_t)(a[k])))); + uw = (((uint32_t)(n->hid))<<1)|((uint32_t)(n->rev)); uw ^= 1; + if(uw == w) e->ou++; + } + } +} + +static void filter_short_ulalignments(void *data, long i, int tid) // callback for kt_for() +{ + const ma_ug_t *ug = (ma_ug_t *)data; + uc_block_t *a = NULL; uc_block_t *p; int64_t k, a_n; uint32_t z, fz, lz, l, bz; + a = UL_INF.a[i].bb.a; a_n = UL_INF.a[i].bb.n; + for (k = a_n - 1; k >= 0; k--) { + p = &(a[k]); + if(p->base || (!p->el) || (!p->pchain)) continue; + if((p->pidx == (uint32_t)-1) && (p->aidx == (uint32_t)-1)) { + if(!ugl_cover_check(p->ts, p->te, &(ug->u.a[p->hid]))) p->pchain = 0; + continue; + } + if(p->pidx == (uint32_t)-1) continue; + if(ugl_cover_check(p->ts, p->te, &(ug->u.a[p->hid]))) continue; + for (z = p->pidx; z != (uint32_t)-1; z = a[z].pidx) { + if(ugl_cover_check(a[z].ts, a[z].te, &(ug->u.a[a[z].hid]))) break; + } + lz = z; fz = p->aidx; l = 0; if(fz != (uint32_t)-1) l = a[fz].pdis; + for (z = k; z != lz; z = bz) { + bz = a[z].pidx; l += a[z].pdis; + a[z].pidx = a[z].pdis = a[z].aidx = (uint32_t)-1; a[z].pchain = 0; + } + + if(fz != (uint32_t)-1 && lz != (uint32_t)-1) { + a[fz].pdis = l; a[fz].pidx = lz; a[lz].aidx = fz; + } else if(fz != (uint32_t)-1) { + a[fz].pdis = a[fz].pidx = (uint32_t)-1; + } else if(lz != (uint32_t)-1) { + a[lz].aidx = (uint32_t)-1; + } + } + + // for (k = a_n - 1; k >= 0; k--) { + // p = &(a[k]); + // if(p->base || (!p->el) || (!p->pchain)) continue; + // if(p->pidx != (uint32_t)-1) { + // assert(a[p->pidx].aidx == (uint32_t)k); + // } + // if(p->aidx != (uint32_t)-1) { + // assert(a[p->aidx].pidx == (uint32_t)k); + // } + // } +} + + +void filter_ul_ug(ma_ug_t *ug) +{ + kt_for(asm_opt.thread_num, filter_short_ulalignments, ug, UL_INF.n); +} + + int32_t find_ul_block_max_reverse(int32_t n, const uc_block_t *a, uint32_t x) { int32_t s = 0, e = n; @@ -8667,6 +8764,7 @@ static void update_ovlp_src(void *data, long i, int tid) // callback for kt_for( // fprintf(stderr, "--i->%d, a_n->%lu--\n", i, a_n); } + uint64_t* get_hifi2ul_list(all_ul_t *x, uint64_t hid, uint64_t* a_n) { (*a_n) = x->ridx.idx.a[hid+1] - x->ridx.idx.a[hid]; @@ -8710,7 +8808,7 @@ int scall_ul_pipeline(uldat_t* sl, const enzyme *fn) fprintf(stderr, "[M::%s::] ==> # reads: %lu, # bases: %lu\n", __func__, UL_INF.n, sl->total_base); fprintf(stderr, "[M::%s::] ==> # bases: %lu; # corrected bases: %lu; # recorrected bases: %lu\n", __func__, sl->num_bases, sl->num_corrected_bases, sl->num_recorrected_bases); - gen_ul_vec_rid_t(&UL_INF); + gen_ul_vec_rid_t(&UL_INF, &R_INF, NULL); return 1; } @@ -8718,13 +8816,7 @@ int scall_ul_pipeline(uldat_t* sl, const enzyme *fn) int rescall_ul_pipeline(uldat_t* sl, const enzyme *fn) { double index_time = yak_realtime(); - int32_t i; uint32_t k, rlen; - ///clean UL_INF - for (k = 0; k < UL_INF.n; k++) { - rlen = UL_INF.a[k].rlen; - free(UL_INF.a[k].bb.a); free(UL_INF.a[k].N_site.a); free(UL_INF.a[k].r_base.a); - memset(&(UL_INF.a[k]), 0, sizeof(UL_INF.a[k])); UL_INF.a[k].rlen = rlen; - } + int32_t i; for (i = 0; i < fn->n; i++){ gzFile fp; @@ -10009,6 +10101,7 @@ int32_t load_all_ul_t(all_ul_t *x, char* file_name, All_reads *hR, ma_ug_t *ug) return 0; } + destory_all_ul_t(x); memset(x, 0, sizeof(*x)); x->hR = hR; init_aux_table(); uint64_t k; ul_vec_t *p = NULL; @@ -10089,7 +10182,7 @@ uint64_t ul_refine_alignment(const ug_opt_t *uopt, asg_t *sg) init_uldat_t(&sl, NULL, NULL, &opt, CHUNK_SIZE, asm_opt.thread_num, uopt, uu); sl.rg = sg; if(work_ul_gchains(&sl)) { free(UL_INF.ridx.idx.a); free(UL_INF.ridx.occ.a); memset(&(UL_INF.ridx), 0, sizeof(UL_INF.ridx)); - gen_ul_vec_rid_t(&UL_INF); + gen_ul_vec_rid_t(&UL_INF, &R_INF, NULL); kt_for(sl.n_thread, update_ovlp_src, &sl, R_INF.total_reads); kt_for(sl.n_thread, update_ovlp_src_bl, &sl, R_INF.total_reads); destroy_ul_idx_t(uu); @@ -10114,6 +10207,18 @@ uint32_t dd_ug(asg_t *sg, ma_ug_t *ug, ma_sub_t* coverage_cut, ma_hit_t_alloc* s } +void clear_all_ul_t(all_ul_t *x) +{ + uint64_t k, rlen; + for (k = 0; k < x->n; k++) { + rlen = x->a[k].rlen; + free(x->a[k].bb.a); free(x->a[k].N_site.a); free(x->a[k].r_base.a); + memset(&(x->a[k]), 0, sizeof(x->a[k])); x->a[k].rlen = rlen; + } + free(x->ridx.idx.a); free(x->ridx.occ.a); memset(&(x->ridx), 0, sizeof((x->ridx))); +} + + ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg) { @@ -10131,11 +10236,14 @@ ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg) // dd_ug(sg, ug, uopt->coverage_cut, uopt->sources, uopt->ruIndex, "UL.sa"); // debug_sl_compress_base_disk_0(&sl, asm_opt.ar); // detect_outlier_len("ul_realignment"); + clear_all_ul_t(&UL_INF); if(!load_all_ul_t(&UL_INF, gfa_name, &R_INF, ug)) { gen_UL_reovlps(&sl, ug, sg, gfa_name, cutoff); write_all_ul_t(&UL_INF, gfa_name, ug); } - + filter_ul_ug(ug); + gen_ul_vec_rid_t(&UL_INF, NULL, ug); + kt_for(asm_opt.thread_num, update_ug_arch_ul, ug, ug->g->n_arc); // print_all_ul_t_stat(&UL_INF); // kt_for(sl.n_thread, update_ovlp_src, &sl, R_INF.total_reads); // kt_for(sl.n_thread, update_ovlp_src_bl, &sl, R_INF.total_reads); diff --git a/inter.h b/inter.h index 257c352..6e4d0c5 100644 --- a/inter.h +++ b/inter.h @@ -10,5 +10,6 @@ uint64_t ul_refine_alignment(const ug_opt_t *uopt, asg_t *sg); ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg); int32_t write_all_ul_t(all_ul_t *x, char* file_name, ma_ug_t *ug); int32_t load_all_ul_t(all_ul_t *x, char* file_name, All_reads *hR, ma_ug_t *ug); +uint32_t ugl_cover_check(uint64_t is, uint64_t ie, ma_utg_t *u); #endif