diff --git a/Correct.cpp b/Correct.cpp index b7c2842..b2054ff 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -3449,7 +3449,7 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea ///debug_scan_cigar(&(overlap_list->list[j])); ///only calculate cigar for high quality overlaps - if ((rref && (overlap_length*OVERLAP_THRESHOLD_FILTER <= z->align_length)) || + if ((rref && (overlap_length*OVERLAP_THRESHOLD_HIFI_FILTER <= z->align_length)) || (uref && (overlap_length*(1-e_rate) <= z->align_length))) { a_nw = z->w_list.n; // int64_t tt = 0; @@ -4178,7 +4178,7 @@ inline void refine_ed_aln(overlap_region_alloc* overlap_list, All_reads *rref, c } } mm_ws = z->x_pos_s; mm_aln = mm_we+1-mm_ws; - if(!simi_pass(ovl, mm_aln, uref?1:0, OVERLAP_THRESHOLD_FILTER, NULL)) continue; + if(!simi_pass(ovl, mm_aln, uref?1:0, OVERLAP_THRESHOLD_NOSI_FILTER, NULL)) continue; if(nw > 0 && w_idx[0] != (uint64_t)-1) mm_ws = z->w_list.a[w_idx[0]].x_end+1; for (i = 1; i < nw; i++) { //utilize the the start pos of next window in backward @@ -4215,13 +4215,13 @@ inline void refine_ed_aln(overlap_region_alloc* overlap_list, All_reads *rref, c } total_y_end = p->y_start - 1; } - if(!simi_pass(ovl, mm_aln, uref?1:0, OVERLAP_THRESHOLD_FILTER, NULL)) break; + if(!simi_pass(ovl, mm_aln, uref?1:0, OVERLAP_THRESHOLD_NOSI_FILTER, NULL)) break; } if(w_idx[i] != (uint64_t)-1) mm_ws = z->w_list.a[w_idx[i]].x_end+1; } if(i < nw) continue; - if(uref && simi_pass(ovl, z->align_length, uref?1:0, OVERLAP_THRESHOLD_FILTER, NULL)) { + if(uref && simi_pass(ovl, z->align_length, uref?1:0, OVERLAP_THRESHOLD_NOSI_FILTER, NULL)) { z->is_match = 3; overlap_list->mapped_overlaps_length += z->align_length; ///sort for set_herror_win if(!is_srt) radix_sort_window_list_xs_srt(z->w_list.a, z->w_list.a + z->w_list.n); @@ -4245,7 +4245,7 @@ inline void refine_ed_aln(overlap_region_alloc* overlap_list, All_reads *rref, c ///debug_scan_cigar(&(overlap_list->list[j])); ///only calculate cigar for high quality overlaps // int64_t tt = 0; - if(simi_pass(ovl, z->align_length, 0, OVERLAP_THRESHOLD_FILTER, &e_rate)) { + if(simi_pass(ovl, z->align_length, 0, OVERLAP_THRESHOLD_NOSI_FILTER, &e_rate)) { a_nw = z->w_list.n; for (i = 0, is_srt = 1; i < a_nw; i++) { p = &(z->w_list.a[i]); @@ -4298,26 +4298,6 @@ inline void refine_ed_aln(overlap_region_alloc* overlap_list, All_reads *rref, c } } -double test_err_rate(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, UC_Read* g_read, -char* qstr, char *tstr, char *tstr_1, Correct_dumy* dumy, kvec_t_u64_warp* v_idx, int64_t block_s, double e_rate) -{ - int64_t ovl = z->x_pos_e+1-z->x_pos_s, a_nw = z->w_list.n, i; double error_rate; - window_list *p; int64_t y_id = z->y_id, y_strand = z->y_pos_strand; - if(!simi_pass(ovl, z->align_length, 0, OVERLAP_THRESHOLD_FILTER, &e_rate)) return DBL_MAX; - - for (i = 0; i < a_nw; i++) { - p = &(z->w_list.a[i]); - ///check if the cigar of this window has been got - if(p->clen == 0) { - gen_backtrace_adv(p, z, rref, hpc_g, uref, qstr, tstr, tstr_1, dumy, y_strand, y_id); - assert(p->error != -1); - } else { - p->y_end -= p->extra_begin; - } - } - error_rate = non_trim_error_rate(z, rref, uref, v_idx, dumy, g_read, e_rate, block_s); - return error_rate; -} uint32_t align_ul_ed_post(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, char* qstr, char *tstr, char *tstr_1, Correct_dumy* dumy, double e_rate, int64_t w_l, double ovlp_cut, void *km); @@ -4332,10 +4312,10 @@ inline void refine_ed_aln_test(overlap_region_alloc* overlap_list, All_reads *rr for (j = 0; j < on; ++j) { z = &(overlap_list->list[j]); ovl = z->x_pos_e+1-z->x_pos_s; if(!align_ul_ed_post(z, uref, NULL, g_read->seq, dumy->overlap_region, dumy->overlap_region_fix, - dumy, e_rate, block_s, OVERLAP_THRESHOLD_FILTER, NULL)) { + dumy, e_rate, block_s, OVERLAP_THRESHOLD_NOSI_FILTER, NULL)) { continue; } - if(uref && simi_pass(ovl, z->align_length, uref?1:0, OVERLAP_THRESHOLD_FILTER, NULL)) { + if(uref && simi_pass(ovl, z->align_length, uref?1:0, OVERLAP_THRESHOLD_NOSI_FILTER, NULL)) { z->is_match = 3; overlap_list->mapped_overlaps_length += z->align_length; } } @@ -16503,6 +16483,9 @@ int64_t extract_subov(int64_t ts, int64_t te, overlap_region *o, double o_rate, if(!ovlp) continue; salnl += ovlp; } + // if(o->y_id == 700) { + // fprintf(stderr, "[M::%s] raw_t::[%ld, %ld), raw_q::[%ld, %ld), salnl::%ld\n", __func__, ts, te, qs, qe, salnl); + // } if((((qe-qs)*o_rate)<=salnl) && (salnl > 0)) return 1; else return 0; } @@ -18166,10 +18149,21 @@ uint64_t rid, rtrace_t *tc, ul_ov_t *res) // } } +void prt_aln_w(overlap_region *o) +{ + uint64_t k; + for (k = 0; k < o->w_list.n; k++) { + fprintf(stderr, "[M::%s::aln->%u] q::[%d, %d), t::[%d, %d), err::%d\n", + __func__, ((o->w_list.a[k].y_end != -1) && (!(is_ualn_win(o->w_list.a[k])))), + o->w_list.a[k].x_start, o->w_list.a[k].x_end+1, + o->w_list.a[k].y_start, o->w_list.a[k].y_end+1, o->w_list.a[k].error); + } +} + ///[ts, te) -> this is the reverse coordinates of t, not the original coordinates of t int64_t extract_subov_cigar(const ul_idx_t *uref, char* qstr, UC_Read *tu, bit_extz_t *exz, int64_t ts0, int64_t te0, overlap_region *o, double o_rate, int64_t *in_k, int64_t ql, -double e_rate, uint64_t rid, rtrace_t *trace, ul_ov_t *res) +double e_rate, uint64_t rid, uint64_t dbg_id, rtrace_t *trace, ul_ov_t *res) { int64_t rev = o->y_pos_strand, t[2], q[2], k = 0, wts, wte, ii[2]; int64_t wn = o->w_list.n, os, oe, ovlp, salnl = 0; @@ -18196,14 +18190,22 @@ double e_rate, uint64_t rid, rtrace_t *trace, ul_ov_t *res) if(k > ii[1]) ii[1] = k; salnl += ovlp; } + // if(dbg_id == 5525) { + // prt_aln_w(o); + // fprintf(stderr, "[M::%s::aln->%ld] ii::[%ld, %ld), q_aln0::[%d, %d), t0::[%ld, %ld)\n", + // __func__, salnl, ii[0], ii[1]+1, o->w_list.a[ii[0]].x_start, o->w_list.a[ii[0]].x_end+1, + // ts0, te0); + // } if((!salnl) || (ii[0] == INT32_MAX) || (ii[1] < 0)) return 0; if((ii[0] == ii[1]) && (is_ualn_win(o->w_list.a[ii[0]]))) return 0; //[ii[0], ii[1]] win_boundary_offset(o->w_list.a, o->w_list.n, ii[0], ts0, ql, &(q[0]), &(t[0])); win_boundary_offset(o->w_list.a, o->w_list.n, ii[1], te0-1, ql, &(q[1]), &(t[1])); q[1]++; t[1]++; - // fprintf(stderr, "\n[M::%s::aln->%ld] ii::[%ld, %ld), q::[%ld, %ld), t::[%ld, %ld)\n", - // __func__, salnl, ii[0], ii[1]+1, q[0], q[1], t[0], t[1]); + // if(dbg_id == 5525) { + // fprintf(stderr, "[M::%s::aln->%ld] ii::[%ld, %ld), q::[%ld, %ld), t::[%ld, %ld)\n", + // __func__, salnl, ii[0], ii[1]+1, q[0], q[1], t[0], t[1]); + // } if(salnl < ((t[1]-t[0])*o_rate)) return 0; gen_raln(uref, qstr, tu, o, exz, ii[0], ii[1]+1, ts0, te0, ql, o->y_id, rev, e_rate, rid, trace, res); @@ -18238,8 +18240,8 @@ uint64_t sid, uint64_t oid, kv_rtrace_t *trace, kv_ul_ov_t *res) if(e <= ts) continue; if(s >= te) break; t[0] = (rev?(u->len-e):(s)); t[1] = (rev?(u->len-s):(e)); - if(extract_subov_cigar(udb, qstr, tu, exz, t[0], t[1], o, o_rate, &k, ql, e_rate, sid, &tz, &z)) { - tz.oid = oid; + if(extract_subov_cigar(udb, qstr, tu, exz, t[0], t[1], o, o_rate, &k, ql, e_rate, sid, rid, &tz, &z)) { + tz.oid = oid; z.ts -= t[0]; z.te -= t[0]; z.rev = ((o->y_pos_strand == ((u->a[i]>>32)&1))?0:1); if(z.rev) { @@ -18257,8 +18259,9 @@ uint64_t sid, uint64_t oid, kv_rtrace_t *trace, kv_ul_ov_t *res) kv_push(rtrace_t, *trace, tz); // fprintf(stderr, "+[M::%s::rid->%lu::rev->%lu] utg_t::[%ld, %ld), ql::%ld\n", // __func__, rid, rev, t[0], t[1], ql); - // fprintf(stderr, "+[M::%s::rid->%u::%c] q::[%u, %u), ql::%ld, t::[%u, %u), tl::%lu, err::%u\n", - // __func__, z.tn, "+-"[z.rev], z.qs, z.qe, ql, z.ts, z.te, Get_READ_LENGTH(R_INF, z.tn), z.sec); + // fprintf(stderr, "+[M::%s::%.*s::%c] q::[%u, %u), ql::%ld, t::[%u, %u), tl::%lu, err::%u\n", + // __func__, (int)Get_NAME_LENGTH(R_INF, z.tn), Get_NAME(R_INF, z.tn), + // "+-"[z.rev], z.qs, z.qe, ql, z.ts, z.te, Get_READ_LENGTH(R_INF, z.tn), z.sec); } } } else { @@ -18267,10 +18270,10 @@ uint64_t sid, uint64_t oid, kv_rtrace_t *trace, kv_ul_ov_t *res) if(e <= ts) continue; if(s >= te) break; t[0] = (rev?(u->len-e):(s)); t[1] = (rev?(u->len-s):(e)); - if(extract_subov_cigar(udb, qstr, tu, exz, t[0], t[1], o, o_rate, &k, ql, e_rate, sid, &tz, &z)) { + if(extract_subov_cigar(udb, qstr, tu, exz, t[0], t[1], o, o_rate, &k, ql, e_rate, sid, rid, &tz, &z)) { tz.oid = oid; z.ts -= t[0]; z.te -= t[0]; - z.rev = ((o->y_pos_strand == ((u->a[i]>>32)&1))?0:1); + z.rev = ((o->y_pos_strand == (ct_a[i].x&1))?0:1); if(z.rev) { t[0] = z.ts; t[1] = z.te; z.ts = Get_READ_LENGTH(R_INF, rid) - t[1]; @@ -18286,8 +18289,9 @@ uint64_t sid, uint64_t oid, kv_rtrace_t *trace, kv_ul_ov_t *res) kv_push(rtrace_t, *trace, tz); // fprintf(stderr, "-[M::%s::rid->%lu::rev->%lu] utg_t::[%ld, %ld), ql::%ld\n", // __func__, rid, rev, t[0], t[1], ql); - // fprintf(stderr, "-[M::%s::rid->%u::%c] q::[%u, %u), ql::%ld, t::[%u, %u), tl::%lu, err::%u\n", - // __func__, z.tn, "+-"[z.rev], z.qs, z.qe, ql, z.ts, z.te, Get_READ_LENGTH(R_INF, z.tn), z.sec); + // fprintf(stderr, "-[M::%s::%.*s::%c] q::[%u, %u), ql::%ld, t::[%u, %u), tl::%lu, err::%u\n", + // __func__, (int)Get_NAME_LENGTH(R_INF, z.tn), Get_NAME(R_INF, z.tn), + // "+-"[z.rev], z.qs, z.qe, ql, z.ts, z.te, Get_READ_LENGTH(R_INF, z.tn), z.sec); } } } @@ -18366,9 +18370,13 @@ kv_rtrace_t *trace, uint64_t ql, uint64_t rid, uint64_t oid, uint64_t khit, over } else { ol = z->x_pos_e+1-z->x_pos_s; aln_ol = z->align_length; if(aln_ol <= min_ovlp) return 0; - // fprintf(stderr, "+[M::%s::aln_ol->%lu] utg%.6dl(%c), align::%u, q::[%u, %u), t::[%u, %u), ql::%lu\n", __func__, aln_ol, - // (int32_t)z->y_id + 1, "+-"[z->y_pos_strand], z->align_length, - // z->x_pos_s, z->x_pos_e+1, z->y_pos_s, z->y_pos_e+1, ql); + // if(z->y_id == 700) { + // fprintf(stderr, "+[M::%s::aln_ol->%lu::o_rate->%f] utg%.6dl(%c), align::%u, q::[%u, %u), t::[%u, %u), ql::%lu\n", + // __func__, aln_ol, o_rate, + // (int32_t)z->y_id + 1, "+-"[z->y_pos_strand], z->align_length, + // z->x_pos_s, z->x_pos_e+1, z->y_pos_s, z->y_pos_e+1, ql); + // print_aln_windows(z); + // } rln_0 = rln->n; if((ol*o_rate) <= aln_ol) { kv_pushp(ul_ov_t, *rln, &p); memset(p, 0, sizeof(*p)); @@ -18385,6 +18393,13 @@ kv_rtrace_t *trace, uint64_t ql, uint64_t rid, uint64_t oid, uint64_t khit, over if(rln->n <= rln_0) return 0; radix_sort_ul_ov_srt_qs1(rln->a+rln_0, rln->a+rln->n); + // if(z->y_id == 700) { + // for (k = 0; k < rln->n; k++) { + // p = &(rln->a[k]); + // fprintf(stderr, "+[M::%s::k->%lu] candidate_q::[%u, %u)\n", __func__, k, p->qs, p->qe); + // } + // } + for (k = mm = rln_0, p = NULL; k < rln->n; k++) { if((!p) || (rln->a[k].qs >= p->qe)) p = NULL; if(p) { @@ -18396,6 +18411,12 @@ kv_rtrace_t *trace, uint64_t ql, uint64_t rid, uint64_t oid, uint64_t khit, over } } rln->n = mm; + // if(z->y_id == 700) { + // for (k = 0; k < rln->n; k++) { + // p = &(rln->a[k]); + // fprintf(stderr, "-[M::%s::k->%lu] candidate_q::[%u, %u)\n", __func__, k, p->qs, p->qe); + // } + // } return_t_chain(z, cl); cigar_gen_by_chain_adv_local(z, cl, rln->a+rln_0, rln->n-rln_0, w_l, udb, NULL, NULL, qstr, tu, exz, aux_o, e_rate, ql, rid, khit); @@ -18468,7 +18489,7 @@ void ul_rid_lalign_adv(overlap_region_alloc* ol, Candidates_list *cl, const ul_i for (i = k = 0; i < ol->length; i++) { z = &(ol->list[i]); z->shared_seed = z->non_homopolymer_errors;///for index if(!ul_local_aln(z, cl, uref, qu->seq, tu, exz, err, w.window_length, - 1000, OVERLAP_THRESHOLD_FILTER, NULL, NULL, NULL, ql, sid, i, khit, NULL)) { + 1000, OVERLAP_THRESHOLD_NOSI_FILTER, NULL, NULL, NULL, ql, sid, i, khit, NULL)) { continue; } if(k != i) { @@ -18483,10 +18504,11 @@ void ul_rid_lalign_adv(overlap_region_alloc* ol, Candidates_list *cl, const ul_i } else { for (i = cln->n = trace->n = 0; i < ol->length; i++) { z = &(ol->list[i]); z->shared_seed = z->non_homopolymer_errors;///for index - // fprintf(stderr, "[M::%s] i::%ld, aln_l::%u, q::[%u, %u)\n", __func__, i, z->align_length, - // z->x_pos_s, z->x_pos_e+1); + // fprintf(stderr, "[M::%s::utg%.6dl(%c)] i::%ld, aln_l::%u, q::[%u, %u), ql::%u\n", __func__, + // (int32_t)z->y_id + 1, "+-"[z->y_pos_strand], i, z->align_length, + // z->x_pos_s, z->x_pos_e+1, z->x_pos_e+1-z->x_pos_s); ul_local_aln(z, cl, uref, qu->seq, tu, exz, err, w.window_length, 1000, - OVERLAP_THRESHOLD_FILTER, aln, cln, trace, ql, sid, i, khit, aux_o); + OVERLAP_THRESHOLD_NOSI_FILTER, aln, cln, trace, ql, sid, i, khit, aux_o); } ///contained reads diff --git a/Hash_Table.h b/Hash_Table.h index 57fceea..d0aebb3 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -11,7 +11,8 @@ ///for one side, the first or last WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY bases should not be corrected #define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25 #define THRESHOLD 15 -#define OVERLAP_THRESHOLD_FILTER 0.9 +#define OVERLAP_THRESHOLD_HIFI_FILTER 0.9 +#define OVERLAP_THRESHOLD_NOSI_FILTER 0.7 #define OVERLAP_THRESHOLD_FILTER_HPC 0.75 #define HIGH_HET_OVERLAP_THRESHOLD_FILTER 0.3 #define HIGH_HET_ERROR_RATE 0.08 diff --git a/Overlaps.cpp b/Overlaps.cpp index aa54243..cbd2c73 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -2695,6 +2695,43 @@ int max_hang, int min_ovlp) +asg_t *ma_sg_gen_ul(ma_hit_t_alloc* sources, int64_t n_read, const ma_sub_t *coverage_cut, +R_to_U* ruIndex, int64_t max_hang, int64_t min_ovlp, int64_t ul_occ) +{ + int64_t i, j, r; asg_arc_t t, *p; const ma_hit_t *h; + asg_t *g = asg_init(); + for (i = 0; i < n_read; ++i) { + asg_seq_set(g, i, coverage_cut[i].e - coverage_cut[i].s, coverage_cut[i].del); + g->seq[i].c = coverage_cut[i].c; + } + CALLOC(g->seq_vis, (g->n_seq<<1)); + // fprintf(stderr, "[M::%s::] n_read::%ld\n", __func__, n_read); + recover_contain_g(g, sources, ruIndex, max_hang, min_ovlp, ul_occ); + + for (i = 0; i < n_read; ++i) { + if(g->seq[i].del) continue; + for (j = 0; j < sources[i].length; j++) { + h = &(sources[i].buffer[j]); + if(h->del) continue; + r = ma_hit2arc(h, (coverage_cut[Get_qn(*h)].e-coverage_cut[Get_qn(*h)].s), + (coverage_cut[Get_tn(*h)].e-coverage_cut[Get_tn(*h)].s), max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + assert(r >= 0); + p = asg_arc_pushp(g); *p = t; + p->ou = ((h->bl>OU_MASK)?OU_MASK:h->bl); + // if(Get_qn(*h) == 10498 && Get_tn(*h) == 10505) { + // fprintf(stderr, "[M::%s::] qn::%u, tn::%u, p->ou::%u, h->bl::%u\n", __func__, + // Get_qn(*h), Get_tn(*h), p->ou, h->bl); + // } + } + } + + asg_cleanup(g); + g->r_seq = g->n_seq; + return g; +} + + + // pop bubbles from vertex v0; the graph MJUST BE symmetric: if u->v present, v'->u' must be present as well //note!!!!!!!! here we don't exculde the deleted edges @@ -5154,6 +5191,94 @@ int asg_arc_del_trans(asg_t *g, int fuzz) return n_reduced; } + +int asg_arc_del_trans_ul(asg_t *g, int fuzz) +{ + uint32_t v, n_vtx = g->n_seq<<1, n_reduced = 0, L, i, nv; + uint8_t *mark; asg_arc_t *av; CALLOC(mark, n_vtx); + uint32_t w, j, nw; asg_arc_t *aw; + + for (v = 0; v < n_vtx; ++v) { + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + if (nv == 0) continue; // no hits + if (g->seq[v>>1].del) { + for (i = 0; i < nv; ++i) av[i].del = 1, ++n_reduced; + continue; + } + /** + p->ul: |____________31__________|__________1___________|______________32_____________| + qns direction of overlap length of this node (not overlap length) + (in the view of query) + p->v : |___________31___________|__________1___________| + tns reverse direction of overlap + (in the view of target) + p->ol: overlap length + **/ + + + //all outnode of v should be set to "not reduce" + for (i = 0; i < nv; ++i) mark[av[i].v] = 1; + + ///length of node (not overlap length) + ///av[nv-1] is longest out-dege + /** + * v--------------- + * w1--------------- + * w2-------------- + * w3-------------- + * w4-------------- + * w5------------- + * for v, the longest out-edge is v->w5 + **/ + L = asg_arc_len(av[nv-1]) + fuzz; + + + for (i = 0; i < nv; ++i) { + //w is an out-node of v + w = av[i].v; + nw = asg_arc_n(g, w); aw = asg_arc_a(g, w); + ///if w has already been reduced + if (mark[av[i].v] != 1) continue; + + for (j = 0; j < nw && asg_arc_len(aw[j]) + asg_arc_len(av[i]) <= L; ++j) + if (mark[aw[j].v]) mark[aw[j].v] = 2; + } + + // for (i = 0; i < nv; ++i) { + // if((av[i].del) || (mark[av[i].v] != 1)) continue; + // w = av[i].v; + // nw = asg_arc_n(g, w); aw = asg_arc_a(g, w); + // for (j = 0; j < nw && asg_arc_len(aw[j]) + asg_arc_len(av[i]) <= L; ++j) { + // if(v == 20996 && w == 21011) { + // fprintf(stderr, "(0):v->%u, ou->%u\n", aw[j].v, aw[j].ou); + // } + + // if (mark[aw[j].v] == 2) { + // if((((uint32_t)av[i].ou) + ((uint32_t)aw[j].ou)) <= OU_MASK) { + // av[i].ou = av[i].ou + aw[j].ou; + // } else { + // av[i].ou = OU_MASK; + // } + // } + // } + // } + //remove edges + for (i = 0; i < nv; ++i) { + if (mark[av[i].v] == 2) { + av[i].del = 1, ++n_reduced; + } + mark[av[i].v] = 0; + } + } + free(mark); + asg_cleanup(g); + asg_symm(g); + // prt_specfic_sge(g, 10498, 10505, __func__); + // normalize_gou(g); + + return n_reduced; +} + ///max_ext is 4 int asg_cut_tip(asg_t *g, int max_ext) { @@ -31463,8 +31588,16 @@ ma_hit_t_alloc* src, uint64_t* readLen, R_to_U* ruIndex, bub_label_t *b_mask_t, ma_hit_cut(src, n_read, readLen, mini_overlap_length, cov); ma_hit_flt(src, n_read, *cov, max_hang_length, mini_overlap_length); ma_hit_contained_advance(src, n_read, *cov, ruIndex, max_hang_length, mini_overlap_length); - sg = ma_sg_gen(src, n_read, *cov, max_hang_length, mini_overlap_length); - asg_arc_del_trans(sg, gap_fuzz); + if(!ul) { + sg = ma_sg_gen(src, n_read, *cov, max_hang_length, mini_overlap_length); + asg_arc_del_trans(sg, gap_fuzz); + } else { + sg = ma_sg_gen_ul(src, n_read, *cov, ruIndex, max_hang_length, mini_overlap_length, UL_COV_THRES); + // prt_specfic_sge(sg, 10498, 10505, "--*--"); + asg_arc_del_trans_ul(sg, gap_fuzz); + // prt_specfic_sge(sg, 10498, 10505, "--#--"); + } + init_bub_label_t(b_mask_t, MIN(10, asm_opt.thread_num), sg->n_seq); asm_opt.coverage = get_coverage(src, *cov, n_read); return sg; @@ -31553,7 +31686,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) 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); + 0.6, asm_opt.max_short_tip, gap_fuzz, &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); diff --git a/Overlaps.h b/Overlaps.h index 30fa7c5..3f7996d 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -212,6 +212,20 @@ typedef struct { ma_utg_t* F_seq; } asg_t; +typedef struct { + ma_hit_t_alloc* src; + int64_t min_ovlp, max_hang, max_hang_rate, need_srt, gap_fuzz; + asg_t *g; + uint32_t *idx; + kvec_t(uint32_t) pi; + asg_arc_t *a; + size_t n, m; +} flex_asg_t; + +typedef struct { + uint32_t i[2]; +}flex_asg_e_retrive_t; + asg_t *asg_init(void); void asg_destroy(asg_t *g); void asg_arc_sort(asg_t *g); diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 43a5f3f..d81ee24 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -22,6 +22,7 @@ KRADIX_SORT_INIT(srt64, uint64_t, generic_key, 8) #define ASG_ET_MULTI_NEI 3 #define UL_TRAV_HERATE 0.2 #define UL_TRAV_FT_RATE 0.8 +#define is_contain_r(ri, z) (((z)<(ri).len)&&((ri).index[(z)]!=(uint32_t)(-1))&&(!((ri).index[(z)]>>31))) KDQ_INIT(uint64_t) @@ -161,6 +162,10 @@ KRADIX_SORT_INIT(ul2ul_srt, ul2ul_t, ul2ul_srt_key, member_size(ul2ul_t, hid)) typedef struct { asg_t *g; ma_hit_t_alloc *src; + R_to_U* ruIndex; + int64_t max_hang; + int64_t min_ovlp; + int64_t ul_occ; } sset_aux; @@ -446,6 +451,25 @@ uint32_t get_arcs(asg_t *g, uint32_t v, uint32_t* idx, uint32_t idx_n) return kv; } +#define flex_arcs0(res, fg, id) ((res) = ((!((id)&(((uint32_t)(0x80000000)))))?(&((fg).g->arc[(id)])):(&((fg).a[(id)])))); + +uint32_t get_flex_arcs(flex_asg_t *fg, uint32_t v, uint32_t* idx, uint32_t idx_n) +{ + asg_t *g = fg->g; uint32_t i, kv = 0; + uint32_t an = asg_arc_n(g, v), beg = g->idx[v]>>32; + for (i = 0, kv = 0; i < an; i++) { + if(g->arc[beg+i].del) continue; + if(idx && kvidx[v]; i != ((uint32_t)-1); i = fg->pi.a[i]) { + if(fg->a[i].del) continue; + if(idx && kvn_seq<<1, v, w, i, k, cnt = 0, nv, kv, pb, ou, mm_ou; + uint32_t n_vtx = g->n_seq<<1, v, w, i, k, cnt = 0, nv, kv, pb, ou, mm_ou, rr, is_u; asg_arc_t *av = NULL; uint64_t lw; if(in) b = in; else b = &tx; @@ -563,6 +587,37 @@ uint32_t asg_arc_cut_tips(asg_t *g, uint32_t max_ext, asg64_v *in, uint32_t is_o } b->n = pb; } + + if(ru && is_ou) { + for (v = b->n = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + av = asg_arc_a(g, v^1); nv = asg_arc_n(g, v^1); + for (i = kv = 0; i < nv; i++) { + if (av[i].del) continue; + kv++; break; + } + if(kv) continue; + + + get_R_to_U(ru, v>>1, &rr, &is_u); + if(rr == (uint32_t)-1 || is_u == 1) continue; + kv_push(uint64_t, *b, v); + + for (i = 0, w = v; i < max_ext; i++) { + if(asg_end(g, w^1, &lw, NULL)!=0) break; + w = (uint32_t)lw; + + get_R_to_U(ru, lw>>1, &rr, &is_u); + if(rr == (uint32_t)-1 || is_u == 1) break; + kv_push(uint64_t, *b, lw); + } + + for (i = 0; i < b->n; i++) { + asg_seq_del(g, ((uint32_t)b->a[i])>>1); + } + if(b->n) cnt++; + } + } /** @@ -596,6 +651,72 @@ uint32_t asg_arc_cut_tips(asg_t *g, uint32_t max_ext, asg64_v *in, uint32_t is_o return cnt; } +static void update_sg_contain(void *data, long i, int tid) +{ + sset_aux *sl = (sset_aux *)data; ma_hit_t *h, *z; asg_arc_t t; + ma_hit_t_alloc *src = sl->src; asg_t *g = sl->g; int64_t r, idx; + R_to_U *ridx = sl->ruIndex; uint32_t k, rr, qn, tn, is_u; + g->seq_vis[i] = 0; + if(!(g->seq[i].del)) return; + get_R_to_U(ridx, i, &rr, &is_u); + if(rr == (uint32_t)-1 || is_u == 1) return; + ma_hit_t_alloc *x = &(src[i]); + for (k = 0; k < x->length; k++) { + h = &(x->buffer[k]); + qn = Get_qn((*h)); tn = Get_tn((*h)); + if(h->bl < sl->ul_occ) continue; + if(g->seq[tn].del) { + get_R_to_U(ridx, tn, &rr, &is_u); + if(rr == (uint32_t)-1 || is_u == 1) continue; + } + r = ma_hit2arc(h, g->seq[qn].len, g->seq[tn].len, sl->max_hang, asm_opt.max_hang_rate, sl->min_ovlp, &t); + if(r < 0) continue; + idx = get_specific_overlap(&(src[tn]), tn, qn); + z = &(src[tn].buffer[idx]); + assert(z->bl == h->bl); + r = ma_hit2arc(z, g->seq[tn].len, g->seq[qn].len, sl->max_hang, asm_opt.max_hang_rate, sl->min_ovlp, &t); + if(r < 0) continue; + h->del = 0; if(!(g->seq[tn].del)) z->del = 0; + g->seq_vis[i] = 1; + } +} + +void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t max_hang, int64_t min_ovlp, int64_t ul_occ) +{ + sset_aux s; s.g = g; s.src = src; s.ruIndex = ruIndex; + s.max_hang = max_hang; s.min_ovlp = min_ovlp; s.ul_occ = ul_occ; + kt_for(asm_opt.thread_num, update_sg_contain, &s, g->n_seq); + uint32_t k; + for (k = 0; k < g->n_seq; k++) { + if(g->seq_vis[k]) g->seq[k].del = 0; + } + memset(g->seq_vis, 0, (sizeof(*(g->seq_vis))*(g->n_seq<<1))); +} + +static void normalize_gou0(void *data, long i, int tid) +{ + sset_aux *sl = (sset_aux *)data; + asg_t *g = sl->g; + asg_arc_t *e = &(g->arc[i]); + if(e->v > (e->ul>>32)) return; + uint32_t k, v = e->v^1, w = (e->ul>>32)^1, ou; + uint32_t nv = asg_arc_n(g, v); + asg_arc_t *av = asg_arc_a(g, v); + for (k = 0; k < nv; ++k) { + if (av[k].v == w) { + ou = MAX(av[k].ou, e->ou); + av[k].ou = e->ou = ou; + break; + } + } +} + +void normalize_gou(asg_t *g) +{ + sset_aux s; s.g = g; + kt_for(asm_opt.thread_num, normalize_gou0, &s, g->n_arc); +} + static void update_sg_uo_t(void *data, long i, int tid) { sset_aux *sl = (sset_aux *)data; @@ -776,7 +897,7 @@ void asg_arc_cut_chimeric(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t o if (cnt > 0) asg_cleanup(g); } -void asg_arc_cut_inexact(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max_ext, uint32_t is_ou, uint32_t is_trio/**, asg64_v *dbg**/) +void asg_arc_cut_inexact(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max_ext, uint32_t is_ou, uint32_t is_trio, float ou_rat/**, asg64_v *dbg**/) { asg64_v tx = {0,0,0}, *b = NULL; uint32_t v, w, i, k, n_vtx = g->n_seq<<1; @@ -845,7 +966,7 @@ void asg_arc_cut_inexact(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max if (kv < 1) continue; if (kv >= 2) { if (mm_ol >= ol_max) continue; - if (is_ou && mm_ou >= ou_max) continue; + if (is_ou && mm_ou > ou_max*ou_rat) continue; } for (i = kw = ol_max = ou_max = 0; i < nw; ++i) { @@ -861,7 +982,7 @@ void asg_arc_cut_inexact(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max if (kw < 1) continue; if (kw >= 2) { if (mm_ol >= ol_max) continue; - if (is_ou && mm_ou >= ou_max) continue; + if (is_ou && mm_ou > ou_max*ou_rat) continue; } if (kv <= 1 && kw <= 1) continue; @@ -1217,6 +1338,374 @@ uint32_t is_topo, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len) if (cnt > 0) asg_cleanup(g); } +asg_arc_t *iter_flex_asg(flex_asg_t *fg, flex_asg_e_retrive_t *rr, uint32_t v) +{ + asg_arc_t *av; uint32_t nv; asg_arc_t *z; + av = asg_arc_a(fg->g, v); nv = asg_arc_n(fg->g, v); + while(rr->i[0] < nv) { + return &(av[rr->i[0]++]); + } + + if(rr->i[0] >= nv && rr->i[0] != ((uint32_t)-1)) { + rr->i[0] = ((uint32_t)-1); rr->i[1] = fg->idx[v]; + } + while (rr->i[1] != ((uint32_t)-1)) { + z = &(fg->a[rr->i[1]]); rr->i[1] = fg->pi.a[rr->i[1]]; + return z; + } + return NULL; +} + +uint32_t detect_tip2(flex_asg_t *fg, uint32_t id, asg64_v *st) +{ + uint32_t v, w, kw; asg_arc_t *sv, *sw; + flex_asg_e_retrive_t rv, rw; + v = id<<1; rv.i[0] = 0; rv.i[1] = (uint32_t)-1; + while(1) { + sv = iter_flex_asg(fg, &rv, v); + if(!sv) break; + if(sv->del) continue; + w = sv->v^1; + if(fg->g->seq_vis[w]&128) continue; + if(fg->g->seq_vis[w>>1]&2) continue; + + rw.i[0] = 0; rw.i[1] = (uint32_t)-1; kw = 0; + while(1) { + sw = iter_flex_asg(fg, &rw, w); + if(!sw) break; + if(sw->del) continue; + if(fg->g->seq_vis[sw->v>>1]&2) continue; + kw++; + } + if(kw < 1) return 1; + fg->g->seq_vis[w] |= 128; kv_push(uint64_t, *st, w); + } + + v = (id<<1)+1; rv.i[0] = 0; rv.i[1] = (uint32_t)-1; + while(1) { + sv = iter_flex_asg(fg, &rv, v); + if(!sv) break; + if(sv->del) continue; + w = sv->v^1; + if(fg->g->seq_vis[w]&128) continue; + if(fg->g->seq_vis[w>>1]&2) continue; + + rw.i[0] = 0; rw.i[1] = (uint32_t)-1; kw = 0; + while(1) { + sw = iter_flex_asg(fg, &rw, w); + if(!sw) break; + if(sw->del) continue; + if(fg->g->seq_vis[sw->v>>1]&2) continue; + kw++; + } + if(kw < 1) return 1; + fg->g->seq_vis[w] |= 128; kv_push(uint64_t, *st, w); + } + + // av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + // for (i = 0; i < nv; ++i) { + // if(av[i].del) continue; + // w = av[i].v^1; + // if(g->seq_vis[w]&128) continue; + // if((g->seq_vis[w>>1]&3)==2) continue; + // aw = asg_arc_a(g, w); nw = asg_arc_n(g, w); + // for (k = kw = 0; k < nw && kw < 1; k++) { + // if(aw[k].del) continue; + // if((g->seq_vis[aw[k].v>>1]&3)==2) continue; + // kw++; + // } + // if(kw < 1) return 1; + // g->seq_vis[w] |= 128; kv_push(uint64_t, *st, w); + // } + + + // v = (id<<1)+1; + // av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + // for (i = 0; i < nv; ++i) { + // if(av[i].del) continue; + // w = av[i].v^1; + // if(g->seq_vis[w]&128) continue; + // if((g->seq_vis[w>>1]&3)==2) continue; + // aw = asg_arc_a(g, w); nw = asg_arc_n(g, w); + // for (k = kw = 0; k < nw && kw < 1; k++) { + // if(aw[k].del) continue; + // if((g->seq_vis[aw[k].v>>1]&3)==2) continue; + // kw++; + // } + // if(kw < 1) return 1; + // g->seq_vis[w] |= 128; kv_push(uint64_t, *st, w); + // } + + return 0; +} + +uint32_t trans_check(flex_asg_t *fg, asg_arc_t *z, asg64_v *b) +{ + flex_asg_e_retrive_t rv, rw; asg_arc_t *sv, *sw; + uint32_t v = z->ul>>32, w; + + rv.i[0] = 0; rv.i[1] = (uint32_t)-1; + while(1) { + sv = iter_flex_asg(fg, &rv, v); + if(!sv) break; + if(sv->del) continue; + if(sv->v == z->v) return 0; + w = sv->v; + // if((z->ul>>33) == 8340 && (z->v>>1) == 8352) { + // fprintf(stderr, "+[M::%s] v>>1::%u(%c), w>>1::%u(%c), sv->v::%u\n", + // __func__, (z->ul>>33), "+-"[(z->ul>>32)&1], (z->v>>1), "+-"[z->v&1], sv->v); + // } + + rw.i[0] = 0; rw.i[1] = (uint32_t)-1; + while (1) { + sw = iter_flex_asg(fg, &rw, w); + if(!sw) break; + if(sw->del) continue; + if(sw->v == z->v) return 0; + } + } + + + uint32_t bn = b->n; uint64_t vl, d, L = asg_arc_len((*z)) + fg->gap_fuzz; + kv_push(uint64_t, *b, v); + while (b->n > bn) { + vl = kv_pop(*b); v = (uint32_t)vl; vl >>= 32; + rv.i[0] = 0; rv.i[1] = (uint32_t)-1; + while(1) { + sv = iter_flex_asg(fg, &rv, v); + if(!sv) break; + if(sv->del) continue; + d = vl + asg_arc_len((*sv)); + if(d > L) continue; + if(sv->v == z->v) { + b->n = bn; + return 0; + } + d <<= 32; d |= sv->v; + kv_push(uint64_t, *b, d); + } + } + + b->n = bn; + return 1; +} + +void push_flex_asg_t(flex_asg_t *fg, asg_arc_t *z) +{ + asg_arc_t *av; uint32_t nv, k, v = z->ul>>32; + av = asg_arc_a(fg->g, v); nv = asg_arc_n(fg->g, v); + for (k = 0; k < nv; k++) { + if(av[k].del) { + av[k] = *z; + if(!(fg->need_srt)) { + if(k > 0 && av[k].ul < av[k-1].ul) fg->need_srt = 1; + if(k+1 < nv && av[k].ul > av[k+1].ul) fg->need_srt = 1; + } + return; + } + } + + k = fg->n; fg->need_srt = 1; + kv_push(asg_arc_t, *fg, *z); + kv_push(uint32_t, fg->pi, fg->idx[v]); fg->idx[v] = k; +} + +void append_notrans_e(flex_asg_t *fg, uint64_t *a, uint64_t a_n, asg64_v *b) +{ + ma_hit_t_alloc* src = fg->src; int32_t idx, r; + uint32_t i, k, v, w; ma_hit_t_alloc *z; asg_arc_t t0, t1; + for (i = 0; i < a_n; i++) { + v = a[i]; + z = &(src[v>>1]); + for (k = i + 1; k < a_n; k++) { + w = a[k]^1; + idx = get_specific_overlap(z, v>>1, w>>1); + if(idx < 0 || z->buffer[idx].del) continue; + r = ma_hit2arc(&(z->buffer[idx]), Get_READ_LENGTH(R_INF, v>>1), Get_READ_LENGTH(R_INF, w>>1), + fg->max_hang, fg->max_hang_rate, fg->min_ovlp, &t0); + if(r < 0) continue; + if((t0.ul>>32) != v || t0.v != w) continue; + t0.ou = ((z->buffer[idx].bl>OU_MASK)?OU_MASK:z->buffer[idx].bl); + // if((v>>1) == 8340 && (w>>1) == 8352) { + // fprintf(stderr, "[M::%s] v>>1::%u(%c), w>>1::%u(%c), t0.ou::%u\n", + // __func__, (v>>1), "+-"[v&1], (w>>1), "+-"[w&1], t0.ou); + // } + + + idx = get_specific_overlap(&(src[w>>1]), w>>1, v>>1); + if(idx < 0 || src[w>>1].buffer[idx].del) continue; + r = ma_hit2arc(&(src[w>>1].buffer[idx]), Get_READ_LENGTH(R_INF, w>>1), Get_READ_LENGTH(R_INF, v>>1), + fg->max_hang, fg->max_hang_rate, fg->min_ovlp, &t1); + if(r < 0) continue; + if((t1.ul>>32) != (w^1) || t1.v != (v^1)) continue; + t1.ou = ((src[w>>1].buffer[idx].bl>OU_MASK)?OU_MASK:src[w>>1].buffer[idx].bl); + // if((v>>1) == 8340 && (w>>1) == 8352) { + // fprintf(stderr, "[M::%s] v>>1::%u(%c), w>>1::%u(%c), t1.ou::%u\n", + // __func__, (v>>1), "+-"[v&1], (w>>1), "+-"[w&1], t1.ou); + // } + + if(trans_check(fg, &t0, b) && trans_check(fg, &t1, b)) { + push_flex_asg_t(fg, &t0); push_flex_asg_t(fg, &t1); + } + } + } +} + +uint32_t iter_contain_g(R_to_U* rI, flex_asg_t *fg, uint32_t v0, asg64_v *b, asg64_v *st) +{ + uint32_t v, w = (uint32_t)-1, x, i, kv, kw, ulen, cnt = 0, st_n, m, is_purge; + asg_arc_t *s; flex_asg_e_retrive_t rr; + b->n = st->n = ulen = 0; kv_push(uint64_t, *st, v0); + while (st->n) { + v = kv_pop(*st); kv = kw = (uint32_t)-1; ulen = 0; + if(fg->g->seq_vis[v>>1]) continue; + while ((!(fg->g->seq_vis[v>>1])) && (is_contain_r((*rI), (v>>1)))) { + kv = get_flex_arcs(fg, v, &w, 1); ulen++; + kv_push(uint64_t, *b, v); + if(kv == 1) { + flex_arcs0(s, (*fg), w); + w = s->v; + kw = get_flex_arcs(fg, w^1, NULL, 0); + if(kw == 1) v = w; + else break; + } else { + break; + } + } + // if((v0>>1) == 12321 || (v0>>1) == 12334) { + // fprintf(stderr, "0[M::%s] v0>>1::%u(%c), b->n::%u, ulen::%u\n", + // __func__, (v0>>1), "+-"[v0&1], (uint32_t)b->n, ulen); + // } + if((!ulen) || (fg->g->seq_vis[v>>1]) || (!(is_contain_r((*rI), (v>>1))))) { + for (i = b->n - ulen; i < b->n; i++) { + fg->g->seq_vis[b->a[i]>>1] = 1; b->a[i] |= ((uint64_t)0x100000000); + } + continue; + } + // fprintf(stderr, "1[M::%s] v0>>1::%u(%c), b->n::%u, ulen::%u\n", + // __func__, (v0>>1), "+-"[v0&1], (uint32_t)b->n, ulen); + + for (i = 0; i < b->n; i++) { + if(b->a[i]&(0x100000000)) continue; + fg->g->seq_vis[b->a[i]>>1] = 2; + } + // fprintf(stderr, "2[M::%s] v0>>1::%u(%c), b->n::%u, ulen::%u\n", + // __func__, (v0>>1), "+-"[v0&1], (uint32_t)b->n, ulen); + for (i = 0, st_n = st->n; i < b->n; i++) { + if(b->a[i]&(0x100000000)) continue; + if(detect_tip2(fg, b->a[i]>>1, st)) break; + } + if(i >= b->n) is_purge = 1; + else is_purge = 0; + // fprintf(stderr, "3[M::%s] v0>>1::%u(%c), b->n::%u, ulen::%u\n", + // __func__, (v0>>1), "+-"[v0&1], (uint32_t)b->n, ulen); + for (i = m = st_n; i < st->n; i++) { + if(fg->g->seq_vis[st->a[i]]&128) fg->g->seq_vis[st->a[i]] -= 128; + if(is_contain_r((*rI), (st->a[i]>>1))) continue; + // if((v0>>1) == 12321 || (v0>>1) == 12334) { + // fprintf(stderr, "bridge::[M::%s] v>>1::%lu(%c)\n", __func__, (st->a[i]>>1), "+-"[st->a[i]&1]); + // } + st->a[m++] = st->a[i]; + } + st->n = m; + // if((v0>>1) == 12321 || (v0>>1) == 12334) { + // fprintf(stderr, "4[M::%s] v0>>1::%u(%c), b->n::%u, i::%u\n", + // __func__, (v0>>1), "+-"[v0&1], (uint32_t)b->n, i); + // } + if(is_purge) { + for (i = 0; i < b->n; i++) { + x = ((uint32_t)b->a[i])>>1; + // if((v0>>1) == 12321 || (v0>>1) == 12334) { + // fprintf(stderr, "del::[M::%s] x::%u, seq_vis::%u\n", __func__, x, fg->g->seq_vis[x]); + // } + if(fg->g->seq_vis[x]&2) { + asg_seq_del(fg->g, x); cnt++; + } + fg->g->seq_vis[x] = 0; + } + append_notrans_e(fg, st->a + st_n, st->n-st_n, b); + st->n = st_n; + return cnt; + } + for (i = 0; i < b->n; i++) { + if(b->a[i]&(0x100000000)) continue; + fg->g->seq_vis[b->a[i]>>1] = 1; + } + st->n = st_n; + + rr.i[0] = 0; rr.i[1] = (uint32_t)-1; + while(1) { + s = iter_flex_asg(fg, &rr, v); + if(!s) break; + if(s->del) continue; + if(fg->g->seq_vis[s->v>>1]) continue; + if(!is_contain_r((*rI), (s->v>>1))) continue; + kv_push(uint64_t, *st, s->v); + } + } + for (i = 0; i < b->n; i++) fg->g->seq_vis[((uint32_t)b->a[i])>>1] = 0; + return cnt; +} + +void flex_asg_t_cleanup(flex_asg_t *fg) +{ + asg_arc_t *p; uint32_t i; + if(fg->n) { + for (i = 0; i < fg->n; i++) { + if(fg->a[i].del) continue; + p = asg_arc_pushp(fg->g); + *p = fg->a[i]; + } + free(fg->g->idx); + fg->g->idx = 0; + fg->g->is_srt = 0; + } else if(fg->need_srt){ + fg->g->is_srt = 0; + } + asg_cleanup(fg->g); +} + +void asg_arc_cut_contain(flex_asg_t *fg, asg64_v *in, asg64_v *in0, R_to_U* rI) +{ + // fprintf(stderr, "+[M::%s]\n", __func__); + asg64_v tx = {0,0,0}, tx0 = {0,0,0}, *b = NULL, *b0 = NULL; + uint32_t v, w = (uint32_t)-1, n_vtx = fg->g->n_seq<<1, cnt = 0; asg_arc_t *s; + b = in?in:&tx; b0 = in0?in0:&tx0; b->n = b0->n = 0; + fg->pi.n = fg->n = 0; memset(fg->idx, -1, (fg->g->n_seq<<1)*sizeof(*(fg->idx))); + fg->need_srt = 0; + + memset(fg->g->seq_vis, 0, sizeof(*(fg->g->seq_vis))*n_vtx); + for (v = 0; v < n_vtx; ++v) { + // if((v>>1) == 6236) { + // fprintf(stderr, "[M::%s] v>>1::%u(%c), del::%u, contain::%u, get_arcs(v)::%u, get_arcs(v^1)::%u\n", + // __func__, (v>>1), "+-"[v&1], g->seq[v>>1].del, is_contain_r((*rI), (v>>1)), + // get_arcs(g, v, NULL, 0), get_arcs(g, v^1, NULL, 0)); + // } + // fprintf(stderr, "+[M::%s] v>>1::%u(%c), del::%u, contain::%u, fg->n::%u, fg->need_srt::%u\n", + // __func__, (v>>1), "+-"[v&1], fg->g->seq[v>>1].del, is_contain_r((*rI), (v>>1)), (uint32_t)fg->n, (uint32_t)fg->need_srt); + if (fg->g->seq[v>>1].del) continue; + if(!is_contain_r((*rI), (v>>1))) continue; + if(get_flex_arcs(fg, v^1, &w, 1) == 1) { + flex_arcs0(s, (*fg), w); + if(get_flex_arcs(fg, s->v^1, NULL, 0) == 1) continue; + } + // if((v>>1) == 12321 || (v>>1) == 12334) { + // fprintf(stderr, "-[M::%s] v>>1::%u(%c), del::%u, contain::%u, fg->n::%u, fg->need_srt::%u\n", + // __func__, (v>>1), "+-"[v&1], fg->g->seq[v>>1].del, is_contain_r((*rI), (v>>1)), (uint32_t)fg->n, (uint32_t)fg->need_srt); + // } + + cnt += iter_contain_g(rI, fg, v, b, b0); + // if((v>>1) == 12321 || (v>>1) == 12334) { + // fprintf(stderr, "*[M::%s] v>>1::%u(%c), del::%u, contain::%u, fg->n::%u, fg->need_srt::%u\n", + // __func__, (v>>1), "+-"[v&1], fg->g->seq[v>>1].del, is_contain_r((*rI), (v>>1)), (uint32_t)fg->n, (uint32_t)fg->need_srt); + // } + } + // stats_sysm(g); + if(!in) free(tx.a); if(!in0) free(tx0.a); + if(cnt > 0) flex_asg_t_cleanup(fg); + // fprintf(stderr, "-[M::%s]\n", __func__); +} + uint32_t if_false_bub_links(uint32_t v, asg_t *g, buf_t *x, asg64_v *b, uint32_t bs, int32_t check_dist) { uint32_t i, mm = 1; @@ -1763,8 +2252,51 @@ void filter_sg_by_ug(asg_t *rg, ma_ug_t *ug, ug_opt_t *uopt) /*******************************for debug************************************/ } +void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd) +{ + uint32_t k, v, w, nv; asg_arc_t *av; + fprintf(stderr, "[M::%s::%s] src::%.*s(id::%u), dst::%.*s(id::%u)\n", __func__, cmd, + (int)Get_NAME_LENGTH(R_INF, src), Get_NAME(R_INF, src), src, + (int)Get_NAME_LENGTH(R_INF, dst), Get_NAME(R_INF, dst), dst); + v = src<<1; + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + for (k = 0; k < nv; ++k) { + if ((av[k].v>>1) == dst) { + w = av[k].v; + fprintf(stderr, "[M::%s::]\t%.*s(%c)\t%.*s(%c)\tou::%u\n", __func__, + (int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), "+-"[v&1], + (int)Get_NAME_LENGTH(R_INF, (w>>1)), Get_NAME(R_INF, (w>>1)), "+-"[w&1], av[k].ou); + } + } + + v = (src<<1)+1; + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + for (k = 0; k < nv; ++k) { + if ((av[k].v>>1) == dst) { + w = av[k].v; + fprintf(stderr, "[M::%s::]\t%.*s(%c)\t%.*s(%c)\tou::%u\n", __func__, + (int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), "+-"[v&1], + (int)Get_NAME_LENGTH(R_INF, (w>>1)), Get_NAME(R_INF, (w>>1)), "+-"[w&1], av[k].ou); + } + } +} + +flex_asg_t *init_flex_asg_t(asg_t *g, ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, int64_t max_hang_rate, int64_t gap_fuzz) +{ + flex_asg_t *z; CALLOC(z, 1); + z->g = g; z->src = src; z->min_ovlp = min_ovlp; z->gap_fuzz = gap_fuzz; + z->max_hang = max_hang; z->max_hang_rate = max_hang_rate; + MALLOC(z->idx, (z->g->n_seq<<1)); memset(z->idx, -1, (z->g->n_seq<<1)*sizeof(*(z->idx))); + return z; +} + +void des_flex_asg_t(flex_asg_t *z) +{ + free(z->idx); free(z->pi.a); free(z->a); +} + 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) +double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, int32_t is_ou, int32_t is_trio, uint32_t ou_thres, char *o_file) { #define HARD_OU_DROP 0.75 #define HARD_OL_DROP 0.6 @@ -1773,41 +2305,54 @@ double ou_drop_rate, int64_t max_tip, bub_label_t *b_mask_t, int32_t is_ou, int3 double step = (clean_round==1?max_ovlp_drop_ratio:((max_ovlp_drop_ratio-min_ovlp_drop_ratio)/(clean_round-1))); double drop = min_ovlp_drop_ratio; - int64_t i; asg64_v bu = {0,0,0}; uint32_t l_drop = 2000; - if(is_ou) update_sg_uo(sg, src); + int64_t i; asg64_v bu = {0,0,0}, ba = {0,0,0}; uint32_t l_drop = 2000; flex_asg_t *fg = NULL; + if(is_ou) fg = init_flex_asg_t(sg, uopt->sources, uopt->min_ovlp, uopt->max_hang, asm_opt.max_hang_rate, gap_fuzz); + // if(is_ou) update_sg_uo(sg, src);///do not do it here + // print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0); + // exit(1); // debug_info_of_specfic_node("m64012_190921_234837/111673711/ccs", sg, rI, "beg"); // debug_info_of_specfic_node("m64011_190830_220126/95028102/ccs", sg, rI, "beg"); - asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); // fprintf(stderr, "[M::%s] count_edges_v_w(sg, 49778, 49847)->%ld\n", __func__, count_edges_v_w(sg, 49778, 49847)); for (i = 0; i < clean_round; i++, drop += step) { if(drop > max_ovlp_drop_ratio) drop = max_ovlp_drop_ratio; // fprintf(stderr, "(0):i->%ld, drop->%f\n", i, drop); + // prt_specfic_sge(sg, 10531, 10519, "--0--"); + // print_vw_edge(sg, 34156, 34090, "0"); // stats_chimeric(sg, src, &bu); if(!is_ou) asg_iterative_semi_circ(sg, src, &bu, max_tip, 1); asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); asg_arc_cut_chimeric(sg, src, &bu, is_ou?ou_thres:(uint32_t)-1); - asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); + // prt_specfic_sge(sg, 10531, 10519, "--1--"); asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); - asg_arc_cut_inexact(sg, src, &bu, max_tip, is_ou, is_trio/**, NULL**//**&dbg**/); + asg_arc_cut_inexact(sg, src, &bu, max_tip, is_ou, is_trio, ou_drop_rate/**, NULL**//**&dbg**/); // debug_edges(&dbg, d, 2); - asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); + // prt_specfic_sge(sg, 10531, 10519, "--2--"); asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); asg_arc_cut_length(sg, &bu, max_tip, drop, ou_drop_rate, is_ou, is_trio, 1, NULL, NULL, NULL); - asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); + + // prt_specfic_sge(sg, 10531, 10519, "--3--"); + if(is_ou) asg_arc_cut_contain(fg, &bu, &ba, rI); asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); asg_arc_cut_bub_links(sg, &bu, HARD_OL_DROP, HARD_OL_SEC_DROP, HARD_OU_DROP, is_ou, asm_opt.large_pop_bubble_size, rev, rI, max_tip); - + // prt_specfic_sge(sg, 10531, 10519, "--4--"); + asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); asg_arc_cut_complex_bub_links(sg, &bu, HARD_OL_DROP, HARD_OU_DROP, is_ou, b_mask_t); - asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); + // prt_specfic_sge(sg, 10531, 10519, "--5--"); + if(is_ou) { if(ul_refine_alignment(uopt, sg)) update_sg_uo(sg, src); @@ -1820,20 +2365,19 @@ double ou_drop_rate, int64_t max_tip, bub_label_t *b_mask_t, int32_t is_ou, int3 asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); asg_cut_large_indel(sg, &bu, max_tip, HARD_OU_DROP, is_ou);///shoule we ignore ou here? - asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); ///asg_arc_del_triangular_directly might be unnecessary asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); asg_arc_cut_length(sg, &bu, max_tip, HARD_ORTHOLOGY_DROP/**min_ovlp_drop_ratio**/, ou_drop_rate, is_ou, 0/**is_trio**/, 0, rev, rI, NULL); - asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); asg_arc_cut_length(sg, &bu, max_tip, min_ovlp_drop_ratio, ou_drop_rate, is_ou, 0/**is_trio**/, 0, rev, rI, &l_drop); - asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL); if(!is_ou) asg_cut_semi_circ(sg, LIM_LEN, 1); - rescue_contained_reads_aggressive(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 10, 1, 0, NULL, NULL, b_mask_t); rescue_missing_overlaps_aggressive(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 1, 0, NULL, b_mask_t); rescue_missing_overlaps_backward(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 10, 1, 0, b_mask_t); @@ -1852,10 +2396,14 @@ double ou_drop_rate, int64_t max_tip, bub_label_t *b_mask_t, int32_t is_ou, int3 if(is_ou) { update_sg_uo(sg, src); } - + if(is_ou) { + des_flex_asg_t(fg); free(fg); + } + print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0); + exit(1); // print_node(sg, 17078); //print_node(sg, 8311); print_node(sg, 8294); - free(bu.a); + free(bu.a); free(ba.a); } diff --git a/gfa_ut.h b/gfa_ut.h index 5c3b6c6..a655e2a 100644 --- a/gfa_ut.h +++ b/gfa_ut.h @@ -3,11 +3,11 @@ #include "Overlaps.h" 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); -uint32_t asg_arc_cut_tips(asg_t *g, uint32_t max_ext, asg64_v *in, uint32_t is_ou); +double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, int32_t is_ou, int32_t is_trio, uint32_t ou_thres, char *o_file); +uint32_t asg_arc_cut_tips(asg_t *g, uint32_t max_ext, asg64_v *in, uint32_t is_ou, R_to_U *ru); void asg_iterative_semi_circ(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t normal_len, uint32_t pop_chimer, asg64_v *dbg); void asg_arc_cut_chimeric(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t ou_thres); -void asg_arc_cut_inexact(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max_ext, uint32_t is_ou, uint32_t is_trio); +void asg_arc_cut_inexact(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max_ext, uint32_t is_ou, uint32_t is_trio, float ou_rat/**, asg64_v *dbg**/); void asg_arc_cut_length(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio, uint32_t is_topo, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len); void asg_arc_cut_bub_links(asg_t *g, asg64_v *in, float len_rat, float sec_len_rat, float ou_rat, uint32_t is_ou, uint64_t check_dist, ma_hit_t_alloc *rev, R_to_U* rI, int32_t max_ext); @@ -16,4 +16,8 @@ uint32_t asg_cut_large_indel(asg_t *g, asg64_v *in, int32_t max_ext, float ou_ra uint32_t asg_cut_semi_circ(asg_t *g, uint32_t lim_len, uint32_t is_clean); void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio); +void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t max_hang, int64_t min_ovlp, int64_t ul_occ); +void normalize_gou(asg_t *g); +void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd); + #endif diff --git a/inter.cpp b/inter.cpp index d80496f..832d777 100644 --- a/inter.cpp +++ b/inter.cpp @@ -6583,26 +6583,36 @@ overlap_region *gen_aux_ovlp(overlap_region_alloc* ol) ///mode: 0->ug; 1->read -int64_t get_ecov_contain_adv(const ul_idx_t *uref, const ug_opt_t *uopt, uint32_t v, uint32_t w, int64_t bw, double diff_ec_ul, int64_t dq) +int64_t get_ecov_contain_adv(const ul_idx_t *uref, const ug_opt_t *uopt, uint32_t v, uint32_t w, int64_t bw, double diff_ec_ul, int64_t dq, int64_t *is_contain) { int64_t dt = -1, dif, mm; ma_hit_t_alloc* src = uopt->sources; int64_t min_ovlp = uopt->min_ovlp; int64_t max_hang = uopt->max_hang; uint64_t z, qn, tn, x = v>>1; int32_t r = 1; asg_arc_t e; + for (z = 0; z < src[x].length; z++) { qn = Get_qn(src[x].buffer[z]); tn = Get_tn(src[x].buffer[z]); if(tn != (w>>1)) continue; r = ma_hit2arc(&(src[x].buffer[z]), Get_READ_LENGTH(R_INF, qn), Get_READ_LENGTH(R_INF, tn), max_hang, asm_opt.max_hang_rate, min_ovlp, &e); if(r >= 0) { if((e.ul>>32) != v || e.v != w) continue; - dt = e.ol; break; - } else if(r == MA_HT_QCONT || r == MA_HT_TCONT) { - if(src[x].buffer[z].rev == ((uint32_t)(v^w))) { + dt = e.ol; if(is_contain) (*is_contain) = 0; + break; + } else if(r == MA_HT_TCONT) {///tn is contained in qn + // if(qn == 543 && tn == 548) { + // fprintf(stderr, "[M::%s]\t%.*s\t->%.*s\tr::%d\tsrc::%c\tqry::%c\n", __func__, (int)Get_NAME_LENGTH(R_INF, qn), + // Get_NAME(R_INF, qn), (int)Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn), r, + // "+-"[src[x].buffer[z].rev], "+-"[((uint32_t)(v^w))]); + // } + if((src[x].buffer[z].rev == (((uint32_t)(v^w))&1)) && (tn == (w>>1))) { dt = Get_qe(src[x].buffer[z]) - Get_qs(src[x].buffer[z]); if(dt < Get_te(src[x].buffer[z]) - Get_ts(src[x].buffer[z])) { dt = Get_te(src[x].buffer[z]) - Get_ts(src[x].buffer[z]); } + if(is_contain) (*is_contain) = 1; + // fprintf(stderr, "[M::%s]\t%.*s\t->%.*s\tdt::%ld\n", __func__, (int)Get_NAME_LENGTH(R_INF, qn), + // Get_NAME(R_INF, qn), (int)Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn), dt); break; } } @@ -6654,7 +6664,7 @@ int64_t trans_sc, All_reads *ridx, char* qstr, UC_Read *tu, int64_t rid, double if(lj->qe <= li->qs) break;//even this pair has a overlap, its length will be very small; just ignore if(lj->qs >= li->qs) continue;///no contain qo = infer_rovlp(li, lj, NULL, NULL, ridx, NULL); ///overlap length in query (UL read) - if(li_v != lj_v && get_ecov_contain_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo)) { + if(li_v != lj_v && get_ecov_contain_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, NULL)) { // fprintf(stderr, "[M::%s] j::%ld, jq::[%u, %u)\n", __func__, j, lj->qs, lj->qe); err = get_rid_backward_cigar_err(&tc, li, trace, NULL, uref, qstr, tu, ol, NULL, exz, e_rate, lj->qe); sc = f[j] + (li->qe - lj->qe) - (err*trans_sc); @@ -6683,7 +6693,7 @@ int64_t trans_sc, All_reads *ridx, char* qstr, UC_Read *tu, int64_t rid, double lj = &(res->a[max_ii]); lj_v = (lj->tn<<1)|lj->rev; if(lj->qe > li->qs && lj->qs < li->qs) { qo = infer_rovlp(li, lj, NULL, NULL, ridx, NULL); ///overlap length in query (UL read) - if(li_v != lj_v && get_ecov_contain_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo)) { + if(li_v != lj_v && get_ecov_contain_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, NULL)) { ///as max_ii < end_j, get_rid_backward_cigar_err still works // fprintf(stderr, "[M::%s] max_ii::%ld, max_ii::[%u, %u)\n", __func__, max_ii, lj->qs, lj->qe); err = get_rid_backward_cigar_err(&tc, li, trace, NULL, uref, qstr, tu, ol, NULL, exz, e_rate, lj->qe); @@ -6754,6 +6764,442 @@ int64_t trans_sc, All_reads *ridx, char* qstr, UC_Read *tu, int64_t rid, double return n_v; } + +void collapse_contain(ul_ov_t *a, int64_t a_n, int64_t i, int64_t *mm_idx, int64_t *mm_sc, int64_t *p, int32_t *s, int64_t min_s) +{ + if((*mm_idx) < 0) return; + +} + +#define rch_connect(x, i) ((((x)>>2)==(i))&&(((x)&2)!=3)) + +int64_t connect_detect(ul_ov_t *a, int64_t a_n, int64_t ai, int64_t aj, All_reads *ridx, int32_t *rch, +const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff, int64_t *p, int32_t *f, rtrace_iter *tc, +kv_rtrace_t *trace, char* qstr, UC_Read *tu, overlap_region_alloc* ol, bit_extz_t *exz, double e_rate, +int64_t trans_sc) +{ + ul_ov_t *li = &(a[ai]), *lj = &(a[aj]), *lk; + uint32_t li_v, lj_v, lk_v; int64_t qo, is_c, ak, afk, err, sc = INT32_MIN; + li_v = (li->tn<<1)|li->rev; lj_v = (lj->tn<<1)|lj->rev; + if(li_v == lj_v) return INT32_MIN; + //even this pair has a overlap, its length will be very small; just ignore + if(lj->qe <= li->qs) return INT32_MIN; + // if(lj->qs >= li->qs) continue;///no contain + + if((rch[aj]>>2) != ai) { + qo = infer_rovlp(li, lj, NULL, NULL, ridx, NULL); ///overlap length in query (UL read) + rch[aj] = (ai<<2); rch[aj] += 3; + // fprintf(stderr, "[j::%ld] (id::%u) %.*s\tqo::%ld\n", aj, lj->tn, + // (int)Get_NAME_LENGTH(R_INF, a[aj].tn), Get_NAME(R_INF, a[aj].tn), qo); + if(get_ecov_contain_adv(uref, uopt, li_v^1, lj_v^1, bw, diff, qo, &is_c)) { + rch[aj] = (ai<<2); rch[aj] += is_c; + } else { + // if(li->tn == 20171) { + // fprintf(stderr, "[j::%ld] %.*s\tconnect::0\n", aj, + // (int)Get_NAME_LENGTH(R_INF, a[aj].tn), Get_NAME(R_INF, a[aj].tn)); + // } + return INT32_MIN; + } + } + // if(li->tn == 20171) { + // fprintf(stderr, "[j::%ld] %.*s\tconnect::%u\n", aj, + // (int)Get_NAME_LENGTH(R_INF, a[aj].tn), Get_NAME(R_INF, a[aj].tn), rch_connect(rch[aj], ai)); + // } + + if(!rch_connect(rch[aj], ai)) return INT32_MIN; + is_c = rch[aj]&1; ak = afk = aj; + // fprintf(stderr, "+[j::%ld] %.*s\tis_c::%ld\n", aj, + // (int)Get_NAME_LENGTH(R_INF, a[aj].tn), Get_NAME(R_INF, a[aj].tn), is_c); + if(is_c) { + for (ak = p[aj]; ak >= 0; ak = p[ak]) { + if((rch[ak]>>2) != ai) { + rch[ak] = (ak<<2); rch[ak] += 3; + lk = &(a[ak]); lk_v = (lk->tn<<1)|lk->rev; + if(li_v == lk_v) break; + if(lk->qe <= li->qs) break; + qo = infer_rovlp(li, lk, NULL, NULL, ridx, NULL); ///overlap length in query (UL read) + if(!get_ecov_contain_adv(uref, uopt, li_v^1, lk_v^1, bw, diff, qo, &is_c)) break; + rch[ak] = (ai<<2); rch[ak] += is_c; + } + if(!rch_connect(rch[ak], ai)) break; + is_c = rch[ak]&1; if(is_c == 0) break; + } + + afk = ak; + if(ak >= 0) { + if(!rch_connect(rch[ak], ai)) return INT32_MIN;///go to a disconnected node + for (ak = p[ak]; ak >= 0 && a[afk].qe <= a[ak].qe + 256; ak = p[ak]) { + if((rch[ak]>>2) != ai) { + rch[ak] = (ak<<2); rch[ak] += 3; + lk = &(a[ak]); lk_v = (lk->tn<<1)|lk->rev; + if(li_v == lk_v) return INT32_MIN; + if(lk->qe <= li->qs) return INT32_MIN; + qo = infer_rovlp(li, lk, NULL, NULL, ridx, NULL); ///overlap length in query (UL read) + if(!get_ecov_contain_adv(uref, uopt, li_v^1, lk_v^1, bw, diff, qo, &is_c)) return INT32_MIN; + rch[ak] = (ai<<2); rch[ak] += is_c; + } + if(!rch_connect(rch[ak], ai)) return INT32_MIN; + } + } + } + // if(li->tn == 20171) { + // fprintf(stderr, "-[j::%ld] %.*s\tis_c::%ld\n", aj, + // (int)Get_NAME_LENGTH(R_INF, a[aj].tn), Get_NAME(R_INF, a[aj].tn), is_c); + // } + + if(afk >= 0) {///reach to one non-contained read + lj = &(a[aj]); + err = get_rid_backward_cigar_err(tc, li, trace, NULL, uref, qstr, tu, ol, NULL, exz, e_rate, lj->qe); + sc = f[aj] + (li->qe - lj->qe) - (err*trans_sc); + } else { + sc = li->qe - li->qs; sc -= (((int64_t)li->sec)*trans_sc); + } + return sc; +} + + +int64_t max_ovlp_src_contain(const ug_opt_t *uopt, uint32_t v) +{ + ma_hit_t_alloc* src = uopt->sources; + int64_t min_ovlp = uopt->min_ovlp, max_hang = uopt->max_hang; + uint32_t i, qn, tn, o = 0, x = v>>1, dt; asg_arc_t e; int32_t r = 1; + + for (i = 0; i < src[x].length; i++) { + qn = Get_qn(src[x].buffer[i]); tn = Get_tn(src[x].buffer[i]); + r = ma_hit2arc(&(src[x].buffer[i]), Get_READ_LENGTH(R_INF, qn), Get_READ_LENGTH(R_INF, tn), + max_hang, asm_opt.max_hang_rate, min_ovlp, &e); + // if(qn == 20171 && tn == 20172) { + // fprintf(stderr, "[r::%d]\t%.*s\t%c\tq::[%u, %u)\t%.*s\tt::[%u, %u)\n", r, + // (int)Get_NAME_LENGTH(R_INF, qn), Get_NAME(R_INF, qn), "+-"[src[x].buffer[i].rev], + // Get_qs(src[x].buffer[i]), Get_qe(src[x].buffer[i]), + // (int)Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn), + // Get_ts(src[x].buffer[i]), Get_te(src[x].buffer[i])); + // } + if(r >= 0) { + if((e.ul>>32) != v) continue; + if(o < e.ol) o = e.ol; + } else if(r == MA_HT_TCONT) {///tn is contained in qn + if(v&1) dt = Get_qe(src[x].buffer[i]); + else dt = Get_READ_LENGTH(R_INF, qn) - Get_qs(src[x].buffer[i]); + if(o < dt) o = dt; + } + } + + return o; +} + +int64_t flat_contain(All_reads *ridx, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, +double diff_ec_ul, int64_t qlen, int64_t max_skip, int64_t max_iter, int64_t max_dis, +ul_ov_t *a, int64_t a_n, int32_t *f, int32_t *c_n, int64_t *p, int64_t *t, ul_ov_t *idx) +{ + if(a_n <= 0) return 0; + int64_t mm_ovlp, x, i, j, st, max_ii, mm_sc, mm_n, mm_idx, n_skip, end_j, qo, sc, sn, is_c, cl; + uint32_t li_v, lj_v; ul_ov_t *li, *lj; int64_t max, max_n, tot_sc = INT32_MIN, tot_n = INT32_MIN, tot_i = -1; + for (i = 1, j = 0; i <= a_n; i++) { + if (i == a_n || a[i].qe != a[j].qe) { + if(i - j > 1) radix_sort_ul_ov_srt_qs(a+j, a+i); + j = i; + } + } + // fprintf(stderr, "\n[M::%s::] sc::%u\n", __func__, idx->qn); + memset(t, 0, (a_n*sizeof((*t)))); + for (i = st = 0, max_ii = -1; i < a_n; ++i) { + li = &(a[i]); li_v = (li->tn<<1)|li->rev; + // fprintf(stderr, "[i::%ld] (id::%u)%.*s\t%c\tq::[%u, %u)\tt::[%u, %u)\tc::%u\n", i, li->tn, + // (int)Get_NAME_LENGTH(R_INF, a[i].tn), Get_NAME(R_INF, a[i].tn), + // "+-"[a[i].rev], a[i].qs, a[i].qe, a[i].ts, a[i].te, !a[i].el); + mm_ovlp = max_ovlp_src(uopt, li_v^1); + x = (li->qs + mm_ovlp)*diff_ec_ul; + if(x < bw) x = bw; + x += li->qs + mm_ovlp; + if (x > qlen+1) x = qlen+1; + x = find_ul_ov_max(i, a, x+G_CHAIN_INDEL); + mm_sc = li->el; mm_n = 1; mm_idx = -1; n_skip = 0; end_j = -1; + if ((x-st) > max_iter) st = x-max_iter; + for (j = x; j >= st; --j) { // collect potential destination vertices + lj = &(a[j]); lj_v = (lj->tn<<1)|lj->rev; + if(lj->qe <= li->qs) break;//even this pair has a overlap, its length will be very small; just ignore + // if(lj->qs >= li->qs) continue;///no contain + qo = infer_rovlp(li, lj, NULL, NULL, ridx, NULL); ///overlap length in query (UL read) + if(li_v != lj_v && get_ecov_contain_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, &is_c)) { + if(is_c == 0) { + sc = f[j] + li->el; sn = c_n[j] + 1; + // if(li->tn == 20171) { + // fprintf(stderr, "[i::%ld] (id::%u)%.*s\t%c\tj::%ld\tsc::%ld\tsn::%ld\n", i, li->tn, + // (int)Get_NAME_LENGTH(R_INF, a[i].tn), Get_NAME(R_INF, a[i].tn), + // "+-"[a[i].rev], j, sc, sn); + // } + if((sc > mm_sc) || ((sc == mm_sc) && (sn > mm_n))) { + mm_sc = sc, mm_idx = j; mm_n = sn; + if (n_skip > 0) --n_skip; + } else if (t[j] == i) { + if (++n_skip > max_skip) + break; + } + if (p[j] >= 0) t[p[j]] = i; + } + } + } + + end_j = j; + if (max_ii < 0 || ((a[i].qe) > (a[max_ii].qe+max_dis))) {//too long + max = INT32_MIN; max_n = INT32_MIN; max_ii = -1; + for (j = i - 1; (j >= st) && (a[i].qe<=(max_dis+a[j].qe)); --j) { + if ((max < f[j]) || ((max == f[j]) && (max_n < c_n[j]))) { + max = f[j]; max_n = c_n[j]; max_ii = j; + } + } + } + + if (max_ii >= 0 && max_ii < end_j) {///just have a try with a[i]<->a[max_ii] + lj = &(a[max_ii]); lj_v = (lj->tn<<1)|lj->rev; + if(lj->qe > li->qs && lj->qs < li->qs) { + qo = infer_rovlp(li, lj, NULL, NULL, ridx, NULL); ///overlap length in query (UL read) + if(li_v != lj_v && get_ecov_contain_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, &is_c)) { + if(is_c == 0) { + sc = f[max_ii] + li->el; sn = c_n[max_ii] + 1; + if((sc > mm_sc) || ((sc == mm_sc) && (sn > mm_n))) { + mm_sc = sc; mm_idx = max_ii; mm_n = sn; + } + } + } + } + } + + f[i] = mm_sc; p[i] = mm_idx; c_n[i] = mm_n; + if ((max_ii < 0) || ((a[i].qe<=max_dis+a[max_ii].qe) && (f[max_ii]= 0) { + t[cl++] = i; i = p[i]; + } + idx->qs = a[t[cl-1]].qs; idx->qe = a[t[cl-1]].qe; + for (i = 0; i < cl; i++) { + a[i] = a[t[cl-i-1]]; + if(a[i].qs < idx->qs) idx->qs = a[i].qs; + if(a[i].qe > idx->qe) idx->qe = a[i].qe; + } + return cl; +} + +///li is the suffix +uint32_t if_qchain_cnn(const ul_idx_t *uref, const ug_opt_t *uopt, All_reads *ridx, int64_t bw, double diff, ul_ov_t *li, ul_ov_t *lj, int64_t *is_c) +{ + uint32_t li_v = (li->tn<<1)|li->rev, lj_v = (lj->tn<<1)|lj->rev; int64_t qo; + if((li_v == lj_v) || (lj->qe <= li->qs)) return 0; + qo = infer_rovlp(li, lj, NULL, NULL, ridx, NULL); ///overlap length in query (UL read) + if(get_ecov_contain_adv(uref, uopt, li_v^1, lj_v^1, bw, diff, qo, is_c)) return 1; + return 0; +} + +void propagate_transitive_reduction(const ul_idx_t *uref, const ug_opt_t *uopt, All_reads *ridx, int64_t bw, double diff, +ul_ov_t *a, int32_t a_n, int64_t ai, int32_t *rch, int32_t *f, int64_t *p, int32_t *c_n, int64_t *t, int64_t *mm_sc, +int64_t *mm_idx, int64_t *mm_n) +{ + if((*mm_idx) < 0) return; + int64_t mm_idx0 = (*mm_idx), j, k, is_c, sn; + for (j = mm_idx0 + 1; j < a_n; j++) { + t[j] = mm_idx0 - 1; + if(p[j] < 0) continue; + if((rch[j]>>2) != ai) { + rch[j] = (ai<<2); rch[j] += 3; + if(if_qchain_cnn(uref, uopt, ridx, bw, diff, &(a[ai]), &(a[j]), &is_c)) { + rch[j] = (ai<<2); rch[j] += is_c; + } + } + if(!rch_connect(rch[j], ai)) continue; + for (k = p[j]; k >= 0 && k > mm_idx0; k = p[k]) { + if(t[k] == mm_idx0) { + k = mm_idx0; break; + } else { + k = mm_idx0-1; break; + } + } + if (k != mm_idx0) continue; + t[j] = mm_idx0;//a[j] could reach mm_idx0; + if(p[j] != k) { + if(!(if_qchain_cnn(uref, uopt, ridx, bw, diff, &(a[j]), &(a[k]), &is_c))) continue; + } + sn = c_n[j] + 1; + if(sn >= (*mm_n)) {//must >= + (*mm_n) = sn; (*mm_idx) = j; + } + } + return; +} + +int64_t gl_rchain_lin_contain(overlap_region_alloc* ol, kv_ul_ov_t *res, ul_ov_t *ex, kv_rtrace_t *trace, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, +double diff_ec_ul, int64_t qlen, int64_t max_skip, int64_t max_iter, int64_t max_dis, Chain_Data* dp, bit_extz_t *exz, +int64_t trans_sc, All_reads *ridx, char* qstr, UC_Read *tu, int64_t rid, double e_rate, int64_t need_srt) +{ + if(res->n == 0) return 0; + uint32_t li_v, rev_n; int32_t *f, *c_n, *c_sc, *rch, *ssc; int64_t *p, *t, res_n = res->n, st, max_ii, max, max_n; + int64_t mm_ovlp, x, i, j, k, sc, csc, mm_sc, mm_idx, mm_n, sn, n_skip, end_j, plus; ul_ov_t *li, *lj, rev_t; rtrace_iter tc; + resize_Chain_Data(dp, res_n, NULL); + t = dp->tmp; f = dp->score; p = dp->pre; c_n = dp->occ; + c_sc = rch = dp->self_length; ssc = dp->indels; + if(need_srt) { + radix_sort_ul_ov_srt_qe(res->a, res->a + res_n); + for (i = 1, j = 0; i <= res_n; i++) { + res->a[i-1].qs = ((uint32_t)-1)-res->a[i-1].qs; + if (i == res_n || res->a[i].qe != res->a[j].qe) { + if(i - j > 1) radix_sort_ul_ov_srt_qs(res->a+j, res->a+i); + j = i; + } + } + } + + memset(t, 0, (res_n*sizeof((*t)))); + for (i = st = plus = 0, max_ii = -1; i < res_n; ++i) { + li = &(res->a[i]); li_v = (li->tn<<1)|li->rev; li->qs = ((uint32_t)-1)-li->qs; + rch[i] = INT32_MAX; ssc[i] = INT32_MIN; + mm_ovlp = max_ovlp_src_contain(uopt, li_v^1); + x = (li->qs + mm_ovlp)*diff_ec_ul; + if(x < bw) x = bw; + x += li->qs + mm_ovlp; + if (x > qlen+1) x = qlen+1; + // if(li->tn == 20171 || li->tn == 20209 || li->tn == 20204) { + // fprintf(stderr, "\n[i::%ld] (id::%u)%.*s\t%c\tq::[%u, %u)\tt::[%u, %u)\tc::%u\tmax_d::%ld\n", i, li->tn, + // (int)Get_NAME_LENGTH(R_INF, res->a[i].tn), Get_NAME(R_INF, res->a[i].tn), + // "+-"[res->a[i].rev], res->a[i].qs, res->a[i].qe, res->a[i].ts, res->a[i].te, + // !res->a[i].el, x+G_CHAIN_INDEL); + // } + x = find_ul_ov_max(i, res->a, x+G_CHAIN_INDEL); + // if(li->tn == 20171 || li->tn == 20209 || li->tn == 20204) { + // fprintf(stderr, "[i::%ld] (id::%u)%.*s\t%c\tq::[%u, %u)\tt::[%u, %u)\tc::%u\tmax_j::%ld\n", i, li->tn, + // (int)Get_NAME_LENGTH(R_INF, res->a[i].tn), Get_NAME(R_INF, res->a[i].tn), + // "+-"[res->a[i].rev], res->a[i].qs, res->a[i].qe, res->a[i].ts, res->a[i].te, + // !res->a[i].el, x); + // } + csc = li->qe - li->qs; csc -= (((int64_t)li->sec)*trans_sc); + mm_sc = csc; mm_idx = -1; mm_n = 1; + n_skip = 0; end_j = -1; tc.k = INT32_MAX; + if ((x-st) > max_iter) st = x-max_iter; + for (j = x; j >= st; --j) { // collect potential destination vertices + lj = &(res->a[j]); + if(lj->qe <= li->qs) break;//even this pair has a overlap, its length will be very small; just ignore + sc = connect_detect(res->a, res->n, i, j, ridx, rch, uref, uopt, bw, diff_ec_ul, p, f, &tc, + trace, qstr, tu, ol, exz, e_rate, trans_sc); + if(sc == INT32_MIN) continue; + sn = c_n[j] + 1; + // if(li->tn == 20209) { + // fprintf(stderr, "[i::%ld] (id::%u)%.*s\t%c\tj::%ld\tsc::%ld\tsn::%ld\n", i, li->tn, + // (int)Get_NAME_LENGTH(R_INF, res->a[i].tn), Get_NAME(R_INF, res->a[i].tn), + // "+-"[res->a[i].rev], j, sc, sn); + // } + if((sc > mm_sc) || ((sc == mm_sc) && (sn > mm_n))) { + mm_sc = sc, mm_idx = j; mm_n = sn; + if (n_skip > 0) --n_skip; + } else if (t[j] == i) { + if (++n_skip > max_skip) + break; + } + if (p[j] >= 0) t[p[j]] = i; + } + + end_j = j; + if (max_ii < 0 || (res->a[i].qe>(res->a[max_ii].qe+max_dis))) {//too long + max = INT32_MIN; max_n = INT32_MIN; max_ii = -1; + for (j = i - 1; (j >= st) && (res->a[i].qe<=(max_dis+res->a[j].qe)); --j) { + if ((max < f[j]) || ((max == f[j]) && (max_n < c_n[j]))) { + max = f[j]; max_n = c_n[j]; max_ii = j; + } + } + } + + if (max_ii >= 0 && max_ii < end_j) {///just have a try with a[i]<->a[max_ii] + lj = &(res->a[max_ii]); + if(lj->qe > li->qs/** && lj->qs < li->qs**/) { + ///as max_ii < end_j, get_rid_backward_cigar_err still works + sc = connect_detect(res->a, res->n, i, max_ii, ridx, rch, uref, uopt, bw, diff_ec_ul, p, f, &tc, + trace, qstr, tu, ol, exz, e_rate, trans_sc); + if(sc != INT32_MIN) { + sn = c_n[max_ii] + 1; + if((sc > mm_sc) || ((sc == mm_sc) && (sn > mm_n))) { + mm_sc = sc; mm_idx = max_ii; mm_n = sn; + } + } + } + } + if(mm_sc < 0) { + mm_sc = csc; mm_idx = -1; mm_n = 1; + } + + if(mm_idx >= 0) { + // propagate_transitive_reduction(res->a, x+1, i, rch, f, p, c_n, &mm_sc, &mm_idx, &mm_n); + propagate_transitive_reduction(uref, uopt, ridx, bw, diff_ec_ul, res->a, x+1, i, rch, f, p, c_n, t, &mm_sc, &mm_idx, &mm_n); + } + // collapse_contain(res->a, res_n, i, &mm_idx, &mm_sc, p, c_sc, end_j); + + f[i] = mm_sc; p[i] = mm_idx; c_n[i] = mm_n; + if(mm_idx < 0 || ssc[mm_idx] < mm_sc) ssc[i] = mm_sc; + else ssc[i] = ssc[mm_idx]; + + if ((max_ii < 0) || ((res->a[i].qe<=max_dis+res->a[max_ii].qe) && (f[max_ii]tn == 20171 || li->tn == 20209 || li->tn == 20204) { + // fprintf(stderr, "[i::%ld]\tf::%d\tp::%ld\n", i, f[i], p[i]); + // } + + } + + for (i = 0; i < res_n; ++i) {///make all f[] positive + ssc[i] -= plus; t[i] = ((uint64_t)ssc[i])<<32; t[i] += (i<<1); + } + + int64_t n_v, n_u, n_v0; + radix_sort_gfa64i(t, t + res_n); plus = 0; + for (k = res_n-1, n_v = n_u = 0; k >= 0; --k) { + n_v0 = n_v; + for (i = ((uint32_t)t[k])>>1; i >= 0 && (t[i]&1) == 0; ) { + ex[n_v++] = res->a[i]; t[i] |= 1; i = p[i]; + } + if(n_v0 == n_v) continue; + sc = (i<0?(t[k]>>32):((t[k]>>32)-f[i])); + c_n[n_u] = n_v-n_v0; c_sc[n_u] = sc; n_u++; if(sc < plus) plus = sc; + } + // fprintf(stderr, "---[M::%s] n_u:%ld, n_v:%ld\n", __func__, n_u, n_v); + for (k = 0, n_v = n_v0 = 0; k < n_u; k++) { + n_v0 = n_v; n_v += c_n[k]; + res->a[k].qn = c_sc[k]-plus;//score + res->a[k].ts = n_v0; res->a[k].te = n_v;///idx + // fprintf(stderr, "[M::%s] k:%ld, c_sc:%d\n", __func__, k, c_sc[k]); + + rev_n = c_n[k]>>1; + ///we need to consider contained reads; so determining qs is not such easy + // res->a[k].qs = (uint32_t)-1; res->a[k].qe = ex[n_v0].qe; + for (i = 0; i < rev_n; i++) { + rev_t = ex[n_v0+i]; ex[n_v0+i] = ex[n_v-i-1]; ex[n_v-i-1] = rev_t; + // if(res->a[k].qs > ex[n_v0+i].qs) res->a[k].qs = ex[n_v0+i].qs; + // if(res->a[k].qs > ex[n_v-i-1].qs) res->a[k].qs = ex[n_v-i-1].qs; + ex[n_v0+i].sec = ex[n_v-i-1].sec = SEC_MODE; + } + if(c_n[k]&1) { + // if(res->a[k].qs > ex[n_v0+i].qs) res->a[k].qs = ex[n_v0+i].qs; + ex[n_v0+i].sec = SEC_MODE; + } + // flat_contain(ex+n_v0, n_v-n_v0); + res->a[k].te = res->a[k].ts + flat_contain(ridx, uref, uopt, bw, diff_ec_ul, qlen, max_skip, max_iter, max_dis, + ex + res->a[k].ts, res->a[k].te - res->a[k].ts, f, ssc, p, t, &(res->a[k])); + } + res->n = n_u; + radix_sort_ul_ov_srt_qn(res->a, res->a + res->n);//sort by score + + // if(res->n > 0) { + // fprintf(stderr, "[M::%s::rid->%ld] qlen::%ld, q::[%u, %u), sc::%u\n", + // __func__, rid, qlen, res->a[res->n-1].qs, res->a[res->n-1].qe, res->a[res->n-1].qn); + // } + return n_v; +} + int64_t select_clean_chain(kv_ul_ov_t *idx, ul_ov_t *res_a, int64_t res_n, int64_t ulid_local, asg64_v *b64) { ul_ov_t kp, *m, *p, *idx_a = idx->a; uint64_t om, ovlp, min_sc, max_sc, ok, z; @@ -6864,6 +7310,20 @@ void prt_rid_raw_chain(kv_ul_ov_t *idx, int64_t rid, int64_t qlen) } +void prt_all_chain(kv_ul_ov_t *idx, ul_ov_t *a, int64_t ql) +{ + uint64_t i, k; + for (i = 0; i < idx->n; i++) { + fprintf(stderr, "\n[M::%s] q::[%u, %u), ql::%ld, sc::%u\n", __func__, idx->a[i].qs, idx->a[i].qe, ql, idx->a[i].qn); + for (k = idx->a[i].ts; k < idx->a[i].te; k++) { + fprintf(stderr, "%.*s\t%c\tq::[%u, %u)\tt::[%u, %u)\tc::%u\n", + (int)Get_NAME_LENGTH(R_INF, a[k].tn), Get_NAME(R_INF, a[k].tn), "+-"[a[k].rev], + a[k].qs, a[k].qe, a[k].ts, a[k].te, !a[k].el); + } + } + +} + void gen_rid_raw_chain(overlap_region_alloc* ol, glchain_t *ll, uint64_t cha_idx, Chain_Data* dp, const ul_idx_t *uref, double diff_ec_ul, int64_t qlen, const ug_opt_t *uopt, char* qstr, UC_Read *tu, bit_extz_t *exz, int64_t ulid_local, int64_t rid, ha_ovec_buf_t *bb) { @@ -6877,9 +7337,12 @@ int64_t rid, ha_ovec_buf_t *bb) memcpy(idx->a, res_a, res_n*sizeof(*(res->a))); // fprintf(stderr, "\n+[M::%s] rid::%ld, name::%.*s\n", __func__, rid, // (int32_t)UL_INF.nid.a[rid].n, UL_INF.nid.a[rid].a); - res_n = gl_rchain_lin(ol, idx, res_a, &(ll->tc), uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, qlen, UG_SKIP_N, UG_ITER_N, UG_DIS_N, dp, exz, tran_sc, &R_INF, qstr, tu, rid, diff_ec_ul, 1); + res_n = gl_rchain_lin_contain(ol, idx, res_a, &(ll->tc), uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, qlen, UG_SKIP_N, UG_ITER_N, UG_DIS_N, dp, exz, tran_sc, &R_INF, qstr, tu, rid, diff_ec_ul, 1); // fprintf(stderr, "-[M::%s] rid::%ld, name::%.*s\n", __func__, rid, // (int32_t)UL_INF.nid.a[rid].n, UL_INF.nid.a[rid].a); + + // prt_all_chain(idx, res_a, qlen); + copy_asg_arr(b64, ll->srt.a); res_n = select_clean_chain(idx, res_a, res_n, ulid_local, &b64); copy_asg_arr(ll->srt.a, b64); @@ -6887,7 +7350,11 @@ int64_t rid, ha_ovec_buf_t *bb) if((idx->n) && (idx->a[0].qe - idx->a[0].qs) >= (qlen*0.95)) { bb->num_read_base++; - } + } + // else { + // // idx->n = 1; + // prt_rid_raw_chain(idx, rid, qlen); + // } // prt_rid_raw_chain(idx, rid, qlen); // //debug @@ -6902,7 +7369,7 @@ static void worker_for_ul_scall_alignment(void *data, long i, int tid) // callba glchain_t *bl = &(s->ll[tid]); int64_t /**rid = s->id+i,**/ winLen = MIN((((double)THRESHOLD_MAX_SIZE)/s->opt->diff_ec_ul), WINDOW), cha_idx; uint32_t high_occ = 2; overlap_region *aux_o = NULL; - // if(s->id+i != 2555) return; + // if(s->id+i != 779) return; // fprintf(stderr, "\n[M::%s] rid::%ld, len::%lu, name::%.*s\n", __func__, s->id+i, s->len[i], // (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a); // if (memcmp(UL_INF.nid.a[s->id+i].a, "d0aab024-b3a7-40fb-83cc-22c3d6d951f8", UL_INF.nid.a[s->id+i].n-1)) return; @@ -12220,6 +12687,45 @@ void destroy_ul_idx_t(ul_idx_t *uu) free(uu); } +void print_raw_u2rgfa_seq(all_ul_t *aln, ul_idx_t *uu, uint32_t is_detail) +{ + uint64_t id, a_n, k, z; uc_block_t *a = NULL; + kvec_t(uint8_t) ff; kv_init(ff); + for (id = 0; id < aln->n; id++) { + a = aln->a[id].bb.a; a_n = aln->a[id].bb.n; + if(a_n == 0) continue; + fprintf(stderr,"\n%.*s\tid::%lu\trlen::%u", (int32_t)aln->nid.a[id].n, aln->nid.a[id].a, id, aln->a[id].rlen); + kv_resize(uint8_t, ff, a_n); memset(ff.a, 0, a_n*sizeof((*(ff.a)))); + if(is_detail) { + fprintf(stderr, "\n"); + for (k = 0; k < a_n; k++) { + if(ff.a[k]) continue; + for (z = k; z != (uint32_t)-1; z = a[z].aidx) { + fprintf(stderr, "%.*s\t%c\tq::[%u, %u)\tt::[%u, %u)\tid::%u\ttl::%lu\n", + (int)Get_NAME_LENGTH(R_INF, a[z].hid), Get_NAME(R_INF, a[z].hid), "+-"[a[z].rev], + a[z].qs, a[z].qe, a[z].ts, a[z].te, a[z].hid, Get_READ_LENGTH(R_INF, a[z].hid)); + assert(ff.a[z] == 0); + ff.a[z] = 1; + } + fprintf(stderr, "************\n"); + } + } else { + fprintf(stderr, "\t"); + for (k = 0; k < a_n; k++) { + if(ff.a[k]) continue; + for (z = k; z != (uint32_t)-1; z = a[z].aidx) { + fprintf(stderr, "%.*s\t", + (int)Get_NAME_LENGTH(R_INF, a[z].hid), Get_NAME(R_INF, a[z].hid)); + assert(ff.a[z] == 0); + ff.a[z] = 1; + } + fprintf(stderr, "\n"); + } + } + } + kv_destroy(ff); +} + void gen_UL_ovlps(uldat_t *sl, int32_t cutoff) { @@ -12229,6 +12735,7 @@ void gen_UL_ovlps(uldat_t *sl, int32_t cutoff) if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name, NULL); sl->ha_flt_tab = ha_flt_tab; sl->ha_idx = (ha_pt_t *)ha_idx; sl->uu = uu; ul_v_call(sl, asm_opt.ar); + // print_raw_u2rgfa_seq(&UL_INF, uu, 1); destroy_ul_idx_t(uu); ha_ft_destroy(ha_flt_tab); ha_pt_destroy(ha_idx); sl->ha_flt_tab = NULL; sl->ha_idx = NULL; sl->uu = NULL; } @@ -12383,8 +12890,6 @@ int32_t load_all_ul_t(all_ul_t *x, char* file_name, All_reads *hR, ma_ug_t *ug) return 1; } - - void ul_load(const ug_opt_t *uopt) { fprintf(stderr, "[M::%s::] ==> UL\n", __func__); @@ -12397,7 +12902,7 @@ void ul_load(const ug_opt_t *uopt) if(!load_all_ul_t(&UL_INF, asm_opt.output_file_name, &R_INF, NULL)) { gen_UL_ovlps(&sl, cutoff); - // write_all_ul_t(&UL_INF, asm_opt.output_file_name, NULL); + write_all_ul_t(&UL_INF, asm_opt.output_file_name, NULL); // exit(1); } // detect_outlier_len("ul_load");