diff --git a/CommandLines.h b/CommandLines.h index cfe6287..b0b583d 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.17.1-r428" +#define HA_VERSION "0.17.2-r432" #define VERBOSE 0 diff --git a/Correct.cpp b/Correct.cpp index dfe257d..b930f07 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -12989,19 +12989,24 @@ int64_t ql, uint64_t rid, ul_ov_t *des) -int64_t push_adp_k_hits(Candidates_list *cl, int64_t cln, uint64_t qs, uint64_t qe, uint64_t ts, uint64_t te, uint64_t readID, int64_t ci) +int64_t push_adp_k_hits(Candidates_list *cl, int64_t cln, uint64_t qs, uint64_t qe, uint64_t ts, uint64_t te, uint64_t readID, int64_t ci, int64_t dbgid) { - int64_t k = ci; k_mer_hit *p = NULL, *m = NULL; + int64_t k = ci, m = -1; k_mer_hit *p = NULL; if(cl->length > cln) { - m = &(cl->list[cl->length-1]); assert(is_alnw(*m)); + m = cl->length-1; assert(is_alnw(cl->list[m])); } for (; k < cln && cl->list[k].readID == readID && cl->list[k].self_offset < qs; k++) { p = &(cl->list[k]); // fprintf(stderr, "[M::%s::] p->q::%ld, p->t::%ld\n", __func__, p->self_offset, p->offset); if(p->offset >= ts) continue; - if((m) && ((m->offset >= p->offset) || (m->self_offset >= p->self_offset))) continue; + if((m>=0) && ((cl->list[m].offset >= p->offset) || (cl->list[m].self_offset >= p->self_offset))) { + continue; + } kv_pushp_cl(k_mer_hit, (*cl), &p); *p = cl->list[k]; p->readID = p->cnt; - // fprintf(stderr, "[M::%s::+]\n", __func__); + // if(dbgid == 109111 || dbgid == 75436) { + // fprintf(stderr, "[M::%s::k->%ld] self_offset::%u, offset::%u, cl->length::%lld, cln::%ld, p->readID::%u\n", + // __func__, k, cl->list[k].self_offset, cl->list[k].offset, cl->length, cln, p->readID); + // } } if(qs == (uint64_t)-1 || qe == (uint64_t)-1) return k; // if(!(k >= cln || cl->list[k].self_offset >= qe)) { @@ -13028,7 +13033,7 @@ int64_t push_adp_k_hits(Candidates_list *cl, int64_t cln, uint64_t qs, uint64_t int64_t gen_weight_khits0(uint32_t qs, uint32_t qe, k_mer_hit *a, int64_t an, int64_t k, uint64_t dp)///[qs, qe) { for (; k >= 0 && a[k].self_offset >= qs; k--); - for (((k>=0)?k:0); k < an && a[k].self_offset < qe; k++) { + for (k = ((k>=0)?k:0); k < an && a[k].self_offset < qe; k++) { if((a[k].self_offset >= qs) && (a[k].self_offset < qe) && (!is_alnw(a[k]))) { a[k].readID = dp; } @@ -13055,7 +13060,25 @@ void gen_weight_khits(asg64_v* idx, k_mer_hit *a, int64_t an) } } -int64_t gen_cns_chain(overlap_region *z, Candidates_list *cl, asg64_v* iidx, int64_t max_lgap, double sgap_rate, int64_t need_filter_khit) +void prt_khit(Candidates_list *cl, overlap_region_alloc* ol, overlap_region *z, uint64_t utg_id, const char *cmd) +{ + uint64_t k, cid; int64_t i; + if(!z) { + for (k = 0; k < ol->length; k++) { + z = &(ol->list[k]); + if(z->y_id == utg_id) break; + } + } + if((z) || (z->y_id == utg_id)) { + fprintf(stderr, "[M::%s::cmd->%s] ******\n", __func__, (char *)cmd); + for (i = z->shared_seed, cid = cl->list[i].readID; i < cl->length && cl->list[i].readID == cid; i++) { + fprintf(stderr, "[M::%s::i->%lu] qoff::%u, toff::%u, cid::%lu, cl->length::%lld\n", __func__, i, + cl->list[i].self_offset, cl->list[i].offset, cid, cl->length); + } + } +} + +int64_t gen_cns_chain(overlap_region_alloc* ol, overlap_region *z, Candidates_list *cl, asg64_v* iidx, int64_t max_lgap, double sgap_rate, int64_t need_filter_khit) { int64_t k, wn = z->w_list.n, aln_n, qs, qe, ts, te, ci, id, rcn = cl->length, kn; k_mer_hit *ka; if(wn <= 0) return 0; @@ -13068,29 +13091,41 @@ int64_t gen_cns_chain(overlap_region *z, Candidates_list *cl, asg64_v* iidx, int qe = z->w_list.a[k].x_end+1; te = z->w_list.a[k].y_end+1; } else { if(qs != -1) { - // fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld, cn::%lld\n", - // __func__, qs, qe, ts, te, ci, cl->length); - ci = push_adp_k_hits(cl, rcn, qs, qe, ts, te, id, ci); + // if(z->y_id == 109111 || z->y_id == 75436) { + // fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld, cn::%lld\n", + // __func__, qs, qe, ts, te, ci, cl->length); + // } + ci = push_adp_k_hits(cl, rcn, qs, qe, ts, te, id, ci, z->y_id); aln_n++;//[qs, qe); [ts, te) } qs = qe = ts = te = -1; } } if(qs != -1) { - // fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld, cn::%lld\n", - // __func__, qs, qe, ts, te, ci, cl->length); - ci = push_adp_k_hits(cl, rcn, qs, qe, ts, te, id, ci); + // if(z->y_id == 109111 || z->y_id == 75436) { + // fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld, cn::%lld\n", + // __func__, qs, qe, ts, te, ci, cl->length); + // } + ci = push_adp_k_hits(cl, rcn, qs, qe, ts, te, id, ci, z->y_id); aln_n++;//[qs, qe); [ts, te) } - // fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld, cn::%lld\n", - // __func__, qs, qe, ts, te, ci, cl->length); - push_adp_k_hits(cl, rcn, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, id, ci); + // if(z->y_id == 109111 || z->y_id == 75436) { + // fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld, cn::%lld, rcn::%ld\n", + // __func__, qs, qe, ts, te, ci, cl->length, rcn); + // } + push_adp_k_hits(cl, rcn, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, id, ci, z->y_id); ka = cl->list+rcn; kn = cl->length-rcn; // if(z->y_id == 66) { // fprintf(stderr, "[M::%s::] kn::%ld, q_pos::%u, t_pos::%u\n", __func__, kn, ka[kn-1].self_offset, ka[kn-1].offset); // } + // if(z->y_id == 109111 || z->y_id == 75436) { + // prt_khit(cl, ol, NULL, 109111, "sa"); + // } if(iidx) gen_weight_khits(iidx, ka, kn); + // if(z->y_id == 109111 || z->y_id == 75436) { + // prt_khit(cl, ol, NULL, 109111, "sb"); + // } if(need_filter_khit) { kn = lchain_dp_trace(ka, kn, max_lgap, sgap_rate, SGAP); cl->length = rcn + kn; } @@ -13313,6 +13348,7 @@ void count_k_hits_filter(overlap_region_alloc* ol, Candidates_list *cl, asg64_v* } + void count_k_hits_adv(All_reads *rref, const ul_idx_t *uref, char* qstr, UC_Read *buf, overlap_region_alloc* ol, Candidates_list *cl, asg64_v* ii, Chain_Data *dp, double sgap_rate, uint64_t khit, uint64_t basec) @@ -13326,7 +13362,7 @@ uint64_t khit, uint64_t basec) refine_khits(&(ol->list[k]), cl, dp, sgap_rate); } radix_sort_bc64(ii->a, ii->a+ii->n); - + // prt_khit(cl, ol, NULL, 109111, "a"); for (k = m = 0; k < on; k++) { z = &(ol->list[(uint32_t)ii->a[k]]); i = z->shared_seed; pid = cl->list[i].readID; z->shared_seed = m; @@ -13336,15 +13372,20 @@ uint64_t khit, uint64_t basec) cl->list[m].cnt = 0; t = cl->list[m].self_offset; t <<= 32; t |= m; kv_push(uint64_t, (*ii), t); + // if(z->y_id == 109111) { + // fprintf(stderr, "[M::%s::m->%ld] qoff::%u, toff::%u\n", __func__, m, + // cl->list[m].self_offset, cl->list[m].offset); + // } m++; } } cl->length = cn = m; srt = ii->a+on; srt_n = ii->n-on; + // prt_khit(cl, ol, NULL, 109111, "b"); radix_sort_bc64(srt, srt+srt_n); if(basec) { resize_UC_Read(buf, (khit<<1)); str0 = buf->seq; str1 = buf->seq + khit; } - + // prt_khit(cl, ol, NULL, 109111, "c"); for (k = 1, l = 0; k <= srt_n; k++) { if(k == cn || (srt[l]>>32) != (srt[k]>>32)) { ff = k - l; @@ -13360,6 +13401,7 @@ uint64_t khit, uint64_t basec) l = k; } } + // prt_khit(cl, ol, NULL, 109111, "d"); } @@ -14436,6 +14478,15 @@ ul_ov_t *res, int64_t wl, int64_t ql, int64_t tl, int64_t ch_i) // __func__, ch_a[iend-1].self_offset, ch_a[iend-1].offset); // } + // if(z->y_id == 109111) { + // int64_t di; + // fprintf(stderr, "[M::%s] ch_n::%ld, s::%lu, e::%lu, ibeg::%ld, iend::%ld\n", __func__, ch_n, s, e, ibeg, iend); + // for (di = 0; di < ch_n; di++) { + // fprintf(stderr, "[M::%s::di->%ld] qoff::%u, toff::%u\n", __func__, di, + // ch_a[di].self_offset, ch_a[di].offset); + // } + // } + if(ibeg >= 0) { res->qs = ch_a[ibeg].self_offset; res->ts = ch_a[ibeg].offset; @@ -14537,10 +14588,15 @@ overlap_region *aux_o, double e_rate) } else { is_done = hc_aln_exz_adv(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, ql, e_rate, MAX_SIN_L, MAX_SIN_E, FORCE_SIN_L, -1, aux_o); } + if(!is_done) { push_unmap_alnw(aux_o, q[0], q[1]-1, t[0], t[1]-1, mode); // chain_win_aln(z, dp, cl, q[0], q[1], t[0], t[1], ql, tl, wl, exz); } + // if(aux_o->y_id == 109111) { + // fprintf(stderr, "<[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld, is_done::%ld\n", + // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode, is_done); + // } } return 0; } @@ -14584,20 +14640,26 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u if((mode == 0) && is_alnw(ch_a[l]) && is_alnw(ch_a[i]) && (ch_a[l].strand == 0) && (ch_a[i].strand == 1)) { is_done = 1; - // fprintf(stderr, "*[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n", - // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); + // if(aux_o->y_id == 109111) { + // fprintf(stderr, "*[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n", + // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); + // } // push_replace_alnw(aux_o, q[0], q[1]-1, t[0], t[1]-1, mode);///no need this push_replace_alnw_adv(z, wl, aux_o, q[0], q[1]-1, t[0], t[1]-1, mode); } else if(mode != 3) { if(mode == 1 || mode == 2) adjust_ext_offset(&(q[0]), &(q[1]), &(t[0]), &(t[1]), ql, tl, 0, mode); - // fprintf(stderr, "#[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n", - // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); + // if(aux_o->y_id == 109111) { + // fprintf(stderr, "#[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n", + // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); + // } is_done = hc_aln_exz_adv(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, ql, e_rate, MAX_CNS_L, MAX_CNS_E, FORCE_CNS_L, -1, aux_o); } if(!is_done) {///postprocess - // fprintf(stderr, ">[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n", - // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); + // if(aux_o->y_id == 109111) { + // fprintf(stderr, ">[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n", + // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); + // } is_done = ovlp_base_aln_all(z, ch_a, ch_n, l, i, uref, hpc_g, rref, qstr, tu, ov, ql, tl, wl, exz, aux_o, e_rate); } // if(rid == (uint64_t)-1) { @@ -15034,26 +15096,45 @@ UC_Read *tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, for (i = 0; i < on; i++) { - // fprintf(stderr, "[M::%s::i->%ld] ovq::[%u, %u), ovt::[%u, %u), hits::[%d, %d)\n", __func__, i, - // ov->qs, ov->qe, ov->ts, ov->te, - // (ov->qn!=((uint32_t)-1))?(int32_t)ov->qn:-1, (int32_t)ov->tn); + // if(aux_o->y_id == 109111) { + // fprintf(stderr, "[M::%s::i->%ld] ovq::[%u, %u), ovt::[%u, %u), hits::[%d, %d)\n", __func__, i, + // ov->qs, ov->qe, ov->ts, ov->te, + // (ov->qn!=((uint32_t)-1))?(int32_t)ov->qn:-1, (int32_t)ov->tn); + // } assert((i<=0)||(ov[i].qs>ov[i-1].qe)); ovlp_base_aln(z, ch_a, ch_n, &(ov[i]), wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, rid); } int64_t aux_n = aux_o->w_list.n; + ///for debug + // if(aux_o->y_id == 109111) { + // for (i = 0; i < ((int64_t)aux_o->w_list.n); i++) { + // fprintf(stderr, "-0-[aln::i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, clen::%u, extra_end::%d, mode::%d\n", i, + // aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start, + // aux_o->w_list.a[i].x_start, aux_o->w_list.a[i].x_end+1, + // aux_o->w_list.a[i].y_start, aux_o->w_list.a[i].y_end+1, + // aux_o->w_list.a[i].error, aux_o->w_list.a[i].clen, + // aux_o->w_list.a[i].extra_end, aux_o->w_list.a[i].error_threshold); + // } + // } + for (i = 0; i < aux_n; i++) { if(!(is_ualn_win(aux_o->w_list.a[i]))) continue; - // if((aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start) <= FORCE_CNS_L) { - // fprintf(stderr, "[aln::i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, clen::%u, mode::%d\n", i, - // aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start, - // aux_o->w_list.a[i].x_start, aux_o->w_list.a[i].x_end+1, - // aux_o->w_list.a[i].y_start, aux_o->w_list.a[i].y_end+1, - // aux_o->w_list.a[i].error, aux_o->w_list.a[i].clen, aux_o->w_list.a[i].error_threshold); - // } //will overwrite ch_a; does not matter rechain_aln(z, cl, aux_o, i, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, tl, h_khit, rid); } + + ///for debug + // if(aux_o->y_id == 109111) { + // for (i = 0; i < ((int64_t)aux_o->w_list.n); i++) { + // fprintf(stderr, "-2-[aln::i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, clen::%u, mode::%d\n", i, + // aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start, + // aux_o->w_list.a[i].x_start, aux_o->w_list.a[i].x_end+1, + // aux_o->w_list.a[i].y_start, aux_o->w_list.a[i].y_end+1, + // aux_o->w_list.a[i].error, aux_o->w_list.a[i].clen, aux_o->w_list.a[i].error_threshold); + // } + // } + if(((int64_t)aux_o->w_list.n) > aux_n) { for (i = m = 0; i < ((int64_t)aux_o->w_list.n); i++) { if((i < aux_n) && (is_ualn_win(aux_o->w_list.a[i]))) continue; @@ -15779,7 +15860,7 @@ void prt_overlap_region_phase_stat(overlap_region *z) } } -void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t *uopt, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, asg64_v* buf1) +void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t *uopt, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, asg64_v* buf1, int64_t ulid) { int64_t on = ol->length, k, i, zwn, q[2], t[2], w[2]; uint64_t m; overlap_region *z; ul_ov_t *cp; @@ -15787,6 +15868,12 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t kv_resize(ul_ov_t, *c_idx, ol->length); for (k = idx->n = c_idx->n = 0; k < on; k++) { z = &(ol->list[k]); zwn = z->w_list.n; + // if(ulid == 35437 && z->y_id == 109111) { + // fprintf(stderr, "+0+[M::%s::utg%.6dl::%c] ulid::%ld, all::%u, non-best::%ld, best::%u\n", __func__, + // (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ulid, z->overlapLen, zwn, z->align_length); + // prt_overlap_region_stat(&(ol->list[k])); + // prt_overlap_region_phase_stat(&(ol->list[k])); + // } // z->align_length = z->overlapLen = z->x_pos_e+1-z->x_pos_s; z->align_length = 0; z->overlapLen = (uint32_t)-1; z->non_homopolymer_errors = 0; @@ -15832,6 +15919,12 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t ovlp_cur_xoff(*cp) = z->w_list.a[w[0]].x_start; ///cur xpos ovlp_cur_coff(*cp) = 0; ///cur cigar off in cur window } + // if(ulid == 35437 && z->y_id == 109111) { + // fprintf(stderr, "+1+[M::%s::utg%.6dl::%c] ulid::%ld, all::%u, non-best::%ld, best::%u\n", __func__, + // (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ulid, z->overlapLen, zwn, z->align_length); + // prt_overlap_region_stat(&(ol->list[k])); + // prt_overlap_region_phase_stat(&(ol->list[k])); + // } } radix_sort_bc64(idx->a, idx->a+idx->n); gen_gov_idx(ol, uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, buf1); @@ -15869,6 +15962,15 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t for (k = 0; k < on; k++) { z = &(ol->list[k]); + // if(ulid == 35437 && z->y_id == 109111) { + // fprintf(stderr, "+2+[M::%s::utg%.6dl::%c] ulid::%ld, all::%u, non-best::%ld, best::%u\n", __func__, + // (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ulid, z->overlapLen, zwn, z->align_length); + // prt_overlap_region_stat(&(ol->list[k])); + // prt_overlap_region_phase_stat(&(ol->list[k])); + // } + + + z->overlapLen = z->x_pos_e+1-z->x_pos_s; z->non_homopolymer_errors = 0; zwn = 0; for (i = m = 0; i < z->align_length; i++) { @@ -15879,7 +15981,14 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t } m++; } - z->w_list.n = m; assert(zwn <= z->overlapLen); + z->w_list.n = m; + // if(ulid == 35437 && z->y_id == 109111/**zwn > z->overlapLen**/) { + // fprintf(stderr, "+3+[M::%s::utg%.6dl::%c] ulid::%ld, all::%u, non-best::%ld, best::%u\n", __func__, + // (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ulid, z->overlapLen, zwn, z->align_length); + // prt_overlap_region_stat(&(ol->list[k])); + // prt_overlap_region_phase_stat(&(ol->list[k])); + // } + assert(zwn <= z->overlapLen); z->align_length = z->overlapLen - zwn; // fprintf(stderr, "[M::%s::utg%.6dl::%c] all::%u, non-best::%ld, best::%u\n", __func__, // (int32_t)z->y_id+1, "+-"[z->y_pos_strand], z->overlapLen, zwn, z->align_length); @@ -15899,14 +16008,17 @@ int64_t max_lgap) int64_t k, l, ch_n, a_n = aln->n; uint64_t pqn, pk; overlap_region *z; //k_mer_hit *ch_a; // count_k_hits(rref, uref, qstr, tu, ol, cl, buf, khit, base_chekc_k_hit); count_k_hits_adv(rref, uref, qstr, tu, ol, cl, buf, &(cl->chainDP), e_rate, khit, base_chekc_k_hit); + // prt_khit(cl, ol, NULL, 109111, "e"); for (k = 1, l = 0, pqn = 0; k <= a_n; k++) { if(k == a_n || aln->a[l].qn != aln->a[k].qn) { z = &(ol->list[aln->a[l].qn]); assert(z->align_length == l); // fprintf(stderr, "[M::%s::] oid::[%lu, %u)\n", __func__, pqn, aln->a[l].qn); for (pk = pqn; pk < aln->a[l].qn; pk++) ol->list[pk].w_list.n = 0; pqn = aln->a[l].qn+1; - - ch_n = gen_cns_chain(z, cl, iidx, max_lgap, e_rate, 0); + // fprintf(stderr, "[M::%s::utg%.6dl::%c]\n", __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand]); + // prt_khit(cl, ol, NULL, 109111, "f"); + ch_n = gen_cns_chain(ol, z, cl, iidx, max_lgap, e_rate, 0); + // prt_khit(cl, ol, NULL, 109111, "g"); if(ch_n) { cigar_gen_by_chain_adv(z, cl, cl->length, ch_n, aln->a+l, k-l, wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, rid, khit); } @@ -16233,7 +16345,7 @@ void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur copy_asg_arr(hap->snp_srt, iidx); copy_asg_arr(v_idx->a, buf); copy_asg_arr((*stb), buf1); copy_asg_arr(iidx, hap->snp_srt); copy_asg_arr(buf, v_idx->a); copy_asg_arr(buf1, (*stb)); - region_phase(ol, uref, uopt, aln, &iidx, &buf, &buf1); + region_phase(ol, uref, uopt, aln, &iidx, &buf, &buf1, sid); copy_asg_arr(hap->snp_srt, iidx); copy_asg_arr(v_idx->a, buf); copy_asg_arr((*stb), buf1); } } diff --git a/Overlaps.cpp b/Overlaps.cpp index 6073176..d048900 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -930,6 +930,7 @@ void normalize_ma_hit_t_single_side_advance(ma_hit_t_alloc* sources, long long n } } + typedef struct { kvec_t_u64_warp *buf; ma_hit_t_alloc* src; diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 783f5e5..8e4ecd9 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -441,7 +441,7 @@ 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)])))); +#define flex_arcs0(res, fg, id) ((res) = ((!((id)&(((uint32_t)(0x80000000)))))?(&((fg).g->arc[(id)])):(&((fg).a[(id)-(((uint32_t)(0x80000000)))])))); uint32_t get_flex_arcs(flex_asg_t *fg, uint32_t v, uint32_t* idx, uint32_t idx_n) { @@ -683,6 +683,38 @@ void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t m memset(g->seq_vis, 0, (sizeof(*(g->seq_vis))*(g->n_seq<<1))); } +// static void update_norm_arc(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; int64_t r, idx; +// uint32_t k, rr, qn, tn, is_u; +// 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 normalize_ma_hit_t_mul(ma_hit_t_alloc *src, uint32_t n_src) +// { +// sset_aux s; s.src = src; +// kt_for(asm_opt.thread_num, update_norm_arc, &s, n_src); +// } + static void normalize_gou0(void *data, long i, int tid) { sset_aux *sl = (sset_aux *)data; @@ -9982,6 +10014,7 @@ int usg_naive_topocut_aux_sec(usg_t *g, uint32_t v0, int max_ext) kv++; } if(kv != 1) break; + if(v == v0) return 0;///circle, it is ok to remove it } v = v0^1; @@ -10007,6 +10040,7 @@ int usg_naive_topocut_aux_sec(usg_t *g, uint32_t v0, int max_ext) kv++; } if(kv != 1) break; + if(v == (v0^1)) return 0;///circle, it is ok to remove it } n_ext -= g->a[v0>>1].occ; @@ -10021,12 +10055,14 @@ int usg_naive_topocut_aux_sec(usg_t *g, uint32_t v0, int max_ext) void usg_arc_cut_length(usg_t *g, asg64_v *in_0, asg64_v *in_1, int32_t max_ext, float len_rat, uint32_t is_trio, uint32_t is_topo, uint32_t *max_drop_len) { - // fprintf(stderr, "+[M::%s::] max_ext::%d, len_rat::%f\n", __func__, max_ext, len_rat); + // if(len_rat > 0.7) { + // fprintf(stderr, "+[M::%s::] max_ext::%d, len_rat::%f\n", __func__, max_ext, len_rat); + // } asg64_v tx = {0,0,0}, tz = {0,0,0}, *b = NULL, *ub = NULL; uint32_t i, k, v, w, n_vtx = g->n<<1, nv, nw, kv, kw, /**trioF = (uint32_t)-1, ntrioF = (uint32_t)-1,**/ ol_max, ou_max, to_del, cnt = 0, mm_ol; usg_arc_t *av, *aw, *ve, *we; uint64_t x, kocc[2], ou; uint8_t *f; CALLOC(f, g->n); b = ((in_0)?(in_0):(&tx)); ub = ((in_1)?(in_1):(&tz)); - + for (v = 0, b->n = ub->n = 0; v < n_vtx; ++v) { if (g->a[v>>1].del) continue; av = usg_arc_a(g, v); nv = usg_arc_n(g, v); @@ -10050,6 +10086,10 @@ uint32_t is_topo, uint32_t *max_drop_len) kv_push(uint64_t, *ub, ((((uint64_t)(v))<<32)|((uint64_t)(i)))); } } + // if(len_rat > 0.7) { + // fprintf(stderr, "[M::%s::] max_ext::%d, len_rat::%f, b->n::%u, ub->n::%u\n", + // __func__, max_ext, len_rat, (uint32_t)b->n, (uint32_t)ub->n); + // } radix_sort_srt64(b->a, b->a + b->n); for (k = 0; k < b->n; k++) { @@ -10099,6 +10139,10 @@ uint32_t is_topo, uint32_t *max_drop_len) } if (kv <= 1 && kw <= 1) continue; + // if(len_rat > 0.7) { + // fprintf(stderr, "0[M::%s::] v>>1::%u(%c), w>>1::%u(%c), kv::%u, kw::%u\n", + // __func__, v>>1, "+-"[v&1], w>>1, "+-"[w&1], kv, kw); + // } to_del = 1; if(is_topo) { @@ -10111,6 +10155,10 @@ uint32_t is_topo, uint32_t *max_drop_len) if (usg_naive_topocut_aux(g, v^1, max_ext, f, b, ub) < max_ext) to_del = 1; } } + // if(len_rat > 0.7) { + // fprintf(stderr, "1[M::%s::] v>>1::%u(%c), w>>1::%u(%c), kv::%u, kw::%u, to_del::%u\n", + // __func__, v>>1, "+-"[v&1], w>>1, "+-"[w&1], kv, kw, to_del); + // } if (to_del) { ve->del = we->del = 1; @@ -10130,6 +10178,10 @@ uint32_t is_topo, uint32_t *max_drop_len) } } + // if(len_rat > 0.7) { + // fprintf(stderr, "2[M::%s::] v>>1::%u(%c), w>>1::%u(%c), kv::%u, kw::%u, to_del::%u\n", + // __func__, v>>1, "+-"[v&1], w>>1, "+-"[w&1], kv, kw, to_del); + // } } if(in_0) free(tx.a); if(in_1) free(tz.a); @@ -10836,7 +10888,7 @@ void integer_realign_g(ul_resolve_t *uidx, usg_t *ng, uinfo_srt_warp_t *seq, uin for (t = 0; t < zt->n; t++) { j = l + t - zt->n; vj = (zt->a[t]<<1)|(seq->a[z-1].v&1); - if(check_hybrid_connect(ng, seq_id, vj, z-1, vi, z)) { + if(check_hybrid_connect(ng, seq_id, vj, z-1, vi, z)) {///seq_id:: integer contig id sc = csc + (((uint32_t)-1) - (track[j]>>32)); if(sc > mm_sc) { mm_sc = sc; mm_idx = j; @@ -10915,7 +10967,9 @@ void integer_realign_g(ul_resolve_t *uidx, usg_t *ng, uinfo_srt_warp_t *seq, uin } buf->u.n = n_u; buf->res_dump.n = l; for (k = 0; k < buf->u.n; k++) { - t = (uint32_t)-1; t <<= 32; t |= (((uint32_t)buf->u.a[k])-(buf->u.a[k]>>32)); + // t = (uint32_t)-1; t <<= 32; + t = seq_id; t <<= 32; t |= ((uint64_t)0x8000000000000000); + t |= (((uint32_t)buf->u.a[k])-(buf->u.a[k]>>32)); kv_push(uint64_t, buf->res_dump, t); n_v0 = buf->u.a[k]>>32; n_v = (uint32_t)buf->u.a[k]; for (z = n_v0; z < n_v; z++) { @@ -10951,7 +11005,7 @@ uint64_t get_ug_occ_v(uint32_t i_ug_occ) return v; } - +///this function might be wrong uint32_t usg_unique_arcs_cluster(asg64_v *b64, uint64_t a_n, uint64_t *idx, uint64_t *integ_seq) { uint64_t bn = b64->n, k, v; @@ -10959,17 +11013,22 @@ uint32_t usg_unique_arcs_cluster(asg64_v *b64, uint64_t a_n, uint64_t *idx, uint memset(b64->a + bn, -1, sizeof((*(b64->a)))*a_n); uint64_t *cidx = b64->a + bn, s, e, i, zs, ze, z; + ///b64.a[0, a_n]:: all resolvable paths with unique beg && end + ///b64.a[a_n, bn]:: (raw unitig/non-unqiue node id)|(resolvable path id) + ///idx:: the idx for b64.a[a_n, bn] for (i = 0; i < a_n; i++) {///available interval with beg/end with unique arcs s = b64->a[i]>>32; e = (uint32_t)(b64->a[i]); assert(e > s); if(cidx[i] != (uint64_t)-1) continue; for (k = s + 1; k < e; k++) {///note: here is [s, e] - v = integ_seq[k]; + v = integ_seq[k];///v is the raw unitig id zs = (idx[v]<<1)>>33; ze = (uint32_t)idx[v]; for (z = zs; z < ze; z++) { // if((b64->a[z]>>32)!=v) { // fprintf(stderr, "[M::%s::] v::%lu, b64->a[z]::%lu, zs::%lu, ze::%lu, a_n::%lu\n", // __func__, v, b64->a[z]>>32, zs, ze, a_n); // } + ///(b64->a[z]>>32):: raw unitig id + ///(uint32_t)b64->a[z]:: available interval id assert((b64->a[z]>>32)==v); cidx[((uint32_t)b64->a[z])] = (i<<32)|((uint32_t)b64->a[z]); } @@ -11001,7 +11060,8 @@ uint32_t usg_unique_arcs_cluster(asg64_v *b64, uint64_t a_n, uint64_t *idx, uint if(k == a_n || (cidx[i]>>32) != (cidx[k]>>32)) { for (z = i; z < k; z++) { assert(cidx[z] != (uint64_t)-1); - b64->a[b64->n++] = (v<<32)|((uint32_t)cidx[z]);///cluest integer seqs-> (cluster id)|(integer seq id) + ///cluest integer seqs-> (cluster id)|(available interval id) + b64->a[b64->n++] = (v<<32)|((uint32_t)cidx[z]); } i = k; v++; } @@ -11011,6 +11071,104 @@ uint32_t usg_unique_arcs_cluster(asg64_v *b64, uint64_t a_n, uint64_t *idx, uint return v;///how many cluster } +void iter_unique_arcs(asg64_v *buf, asg64_v *b64, uint64_t a_n, uint64_t *idx, uint64_t *integ_seq, uint64_t *cidx, uint64_t i0) +{ + uint64_t s, e, x, k, v, zs, ze, z; + buf->n = 0; + kv_push(uint64_t, *buf, i0); + while(buf->n) { + x = kv_pop(*buf); + if(cidx[x] != (uint64_t)-1) continue; + cidx[x] = (i0<<32)|(x); + s = b64->a[x]>>32; e = (uint32_t)(b64->a[x]); assert(e > s); + for (k = s + 1; k < e; k++) {///note: here is [s, e] + v = integ_seq[k];///v is the raw unitig id + zs = (idx[v]<<1)>>33; ze = (uint32_t)idx[v]; + for (z = zs; z < ze; z++) { + ///(b64->a[z]>>32):: raw unitig id + ///(uint32_t)b64->a[z]:: available interval id + assert((b64->a[z]>>32)==v); + if(cidx[((uint32_t)b64->a[z])] != (uint64_t)-1) { + assert((cidx[((uint32_t)b64->a[z])]>>32)==i0); + continue; + } + // cidx[((uint32_t)b64->a[z])] = (i0<<32)|((uint32_t)b64->a[z]); + kv_push(uint64_t, *buf, ((uint32_t)b64->a[z])); + } + + v = integ_seq[k]^1; + zs = (idx[v]<<1)>>33; ze = (uint32_t)idx[v]; + for (z = zs; z < ze; z++) { + ///(b64->a[z]>>32):: raw unitig id + ///(uint32_t)b64->a[z]:: available interval id + assert((b64->a[z]>>32)==v); + if(cidx[((uint32_t)b64->a[z])] != (uint64_t)-1) { + assert((cidx[((uint32_t)b64->a[z])]>>32)==i0); + continue; + } + // cidx[((uint32_t)b64->a[z])] = (i0<<32)|((uint32_t)b64->a[z]); + kv_push(uint64_t, *buf, ((uint32_t)b64->a[z])); + } + } + + v = integ_seq[s]; + zs = (idx[v]<<1)>>33; ze = (uint32_t)idx[v]; + for (z = zs; z < ze; z++) { + assert((b64->a[z]>>32)==v); + if(cidx[((uint32_t)b64->a[z])] != (uint64_t)-1) { + assert((cidx[((uint32_t)b64->a[z])]>>32)==i0); + continue; + } + // cidx[((uint32_t)b64->a[z])] = (i0<<32)|((uint32_t)b64->a[z]); + kv_push(uint64_t, *buf, ((uint32_t)b64->a[z])); + } + + v = integ_seq[e]^1; + zs = (idx[v]<<1)>>33; ze = (uint32_t)idx[v]; + for (z = zs; z < ze; z++) { + assert((b64->a[z]>>32)==v); + if(cidx[((uint32_t)b64->a[z])] != (uint64_t)-1) { + assert((cidx[((uint32_t)b64->a[z])]>>32)==i0); + continue; + } + // cidx[((uint32_t)b64->a[z])] = (i0<<32)|((uint32_t)b64->a[z]); + kv_push(uint64_t, *buf, ((uint32_t)b64->a[z])); + } + } +} + +uint32_t usg_unique_arcs_cluster_adv(asg64_v *b64, uint64_t a_n, uint64_t *idx, uint64_t *integ_seq, asg64_v *buf) +{ + uint64_t bn = b64->n, k, z, v; buf->n = 0; + kv_resize(uint64_t, *b64, b64->n + a_n); b64->n += a_n; + memset(b64->a + bn, -1, sizeof((*(b64->a)))*a_n); + uint64_t *cidx = b64->a + bn, i; + + ///b64.a[0, a_n]:: all resolvable paths with unique beg && end + ///b64.a[a_n, bn]:: (raw unitig/non-unqiue node id)|(resolvable path id) + ///idx:: the idx for b64.a[a_n, bn] + for (i = 0; i < a_n; i++) {///available interval with beg/end with unique arcs + if(cidx[i] != (uint64_t)-1) continue; + iter_unique_arcs(buf, b64, a_n, idx, integ_seq, cidx, i); + } + + radix_sort_srt64(cidx, cidx + a_n); b64->n = a_n; + for (i = 0, k = 1, v = 0; k <= a_n; k++) { + if(k == a_n || (cidx[i]>>32) != (cidx[k]>>32)) { + for (z = i; z < k; z++) { + assert(cidx[z] != (uint64_t)-1); + ///cluest integer seqs-> (cluster id)|(available interval id) + b64->a[b64->n++] = (v<<32)|((uint32_t)cidx[z]); + } + i = k; v++; + } + } + + buf->n = 0; + assert(b64->n == (a_n<<1)); + return v;///how many cluster +} + uint32_t ava_pass_unique_bridge(uint64_t *idx, uint64_t *integer_seq, uint64_t s, uint64_t e) { uint64_t k; @@ -11030,14 +11188,15 @@ uint32_t ava_pass_unique_bridge(uint64_t *idx, uint64_t *integer_seq, uint64_t s uint32_t ava_pass_unique_bridge_tips(usg_t *g, asg64_v *b64, uint64_t g_s, uint64_t g_e, uint64_t *integer_seq, uint8_t *f, uint64_t max_ext) { uint64_t i, k, z, s, e, v, nv, bn = b64->n, kv, n_ext = 0; usg_arc_t *av; - for (i = g_s; i < g_e; i++) { + for (i = g_s; i < g_e; i++) {///available intervals within the same cluster s = b64->a[((uint32_t)b64->a[i])]>>32; e = ((uint32_t)b64->a[((uint32_t)b64->a[i])]); assert(s < e); for (k = s + 1; k < e; k++) {///note: here is [s, e] f[integer_seq[k]] = f[integer_seq[k]^1] = 1; } f[integer_seq[s]] = f[integer_seq[e]^1] = 1; } - + + ///collect nodes within raw unitig graph that are linked by the clusters but not in the cluster for (i = g_s; i < g_e; i++) { s = b64->a[((uint32_t)b64->a[i])]>>32; e = ((uint32_t)b64->a[((uint32_t)b64->a[i])]); assert(s < e); for (k = s + 1; k < e; k++) {///note: here is [s, e] @@ -11233,7 +11392,7 @@ void update_usg_t_threading_0(usg_t *ng, uint64_t *a, uint64_t a_n, uint32_t *oc void update_usg_t_threading(ul_resolve_t *uidx, usg_t *ng, uint64_t *arcs, uint64_t *arcs_g, uint64_t arcs_gn, uint64_t *integ_seq, uint32_t *occ, asg64_v *b) { uint64_t i, k, s, e, nvtx = ng->n<<1; memset(occ, 0, sizeof((*occ))*nvtx); - for (i = 0; i < arcs_gn; i++) { + for (i = 0; i < arcs_gn; i++) {///set cluster s = arcs[((uint32_t)arcs_g[i])]>>32; e = ((uint32_t)arcs[((uint32_t)arcs_g[i])]); assert(s < e); for (k = s + 1; k < e; k++) {///note: here is [s, e] occ[integ_seq[k]]++; occ[integ_seq[k]^1]++; @@ -11525,10 +11684,43 @@ ma_ug_t *ma_ug_hybrid_gen(usg_t *g) +void prt_thread_info(uint64_t *interval, uint64_t interval_n, uint64_t *cluster, uint64_t cluster_n, +const char *nn) +{ + char* gfa_name = NULL; MALLOC(gfa_name, strlen(nn)+70); + sprintf(gfa_name, "%s.thread_info.log", nn); + FILE* fp = fopen(gfa_name, "w"); free(gfa_name); + if (!fp) return; + uint64_t k; + if(interval) { + for (k = 0; k < interval_n; k++) fprintf(fp,"it_val::[%lu, %u)\n", interval[k]>>32, (uint32_t)interval[k]); + } + if(cluster) { + for (k = 0; k < cluster_n; k++) fprintf(fp,"cluster::%lu\tit_id::%u)\n", cluster[k]>>32, (uint32_t)cluster[k]); + } + fclose(fp); +} + +void prt_intg_info(uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, const char *nn) +{ + char* gfa_name = NULL; MALLOC(gfa_name, strlen(nn)+70); + sprintf(gfa_name, "%s.intg_info.log", nn); + FILE* fp = fopen(gfa_name, "w"); free(gfa_name); + if (!fp) return; + uint64_t k, i, s, e; + for (i = 0; i < int_idx_n; i++) {///scan all integer contigs + s = int_idx[i]>>32; e = s + ((uint32_t)int_idx[i]); + assert(e > s + 1);//the length is at least 2 + fprintf(fp,"idx::[%lu, %lu)\n", s, e); + for (k = s; k < e; k++) { + fprintf(fp,"%lu\n", int_a[k]); + } + } + fclose(fp); +} - - +void prt_usg_t(ul_resolve_t *uidx, usg_t *ng, const char *cmd); uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, uint32_t max_ext) { @@ -11548,6 +11740,8 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint } ma_ug_t *un_g = ma_ug_hybrid_gen(ng); ma_utg_t *u; + // print_debug_gfa(uidx->sg, un_g, uidx->uopt->coverage_cut, "iig0", uidx->uopt->sources, + // uidx->uopt->ruIndex, uidx->uopt->max_hang, uidx->uopt->min_ovlp, 0, 0, 0); int32_t ui, un; for (k = 0; k < un_g->u.n; k++) {///all unitigs of raw utg u = &(un_g->u.a[k]); @@ -11570,7 +11764,7 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint ma_ug_destroy(un_g); // fprintf(stderr, ">>>>>>[M::%s::] int_idx[0]::%lu, int_idx[1]::%lu\n", __func__, int_idx[0], int_idx[1]); - + // prt_intg_info(int_idx, int_idx_n, int_a, "intg"); for (i = b64.n = ua_n = a_n = 0; i < int_idx_n; i++) {///scan all integer contigs s = int_idx[i]>>32; e = s + ((uint32_t)int_idx[i]); assert(e > s + 1);//the length is at least 2 @@ -11603,12 +11797,13 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint } } assert(a_n + ua_n == b64.n); + // prt_thread_info(b64.a, a_n+ua_n, NULL, 0, "tt_minus"); // fprintf(stderr, "**0**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); uint64_t *i_idx, n_clus; CALLOC(i_idx, ng->n<<1); radix_sort_srt64(b64.a, b64.a + b64.n);///keeps the coordinates within int_a[] - + // prt_thread_info(b64.a, a_n, NULL, 0, "tt0"); /*********debugging*********/ // for (i = 0; i < a_n; i++) { ///available intervals // s = b64.a[i]>>32; e = (uint32_t)b64.a[i]; @@ -11641,6 +11836,7 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint i_idx[int_a[s]] |= ((uint64_t)0x8000000000000000); i_idx[int_a[e]^1] |= ((uint64_t)0x8000000000000000); } + // prt_thread_info(b64.a, a_n, NULL, 0, "tt1"); ///unavailable intervals are useless b64.n = a_n; assert(ua_n == 0); for (i = 0; i < a_n; i++) { ///available intervals @@ -11648,7 +11844,7 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint assert(!(b64.a[i]&((uint64_t)0x8000000000000000))); for (k = s + 1; k < e; k++) {///note: here is [s, e]; s && e are unique, but [s+1, e-1] are not unique kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[k]]++; - (*pz) = int_a[k]; (*pz) <<= 32; (*pz) |= i;///(raw unitig node id)|(integer contig id) + (*pz) = int_a[k]; (*pz) <<= 32; (*pz) |= i;///(raw unitig/non-unqiue node id)|(integer contig id) kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[k]^1]++; (*pz) = int_a[k]^1; (*pz) <<= 32; (*pz) |= i; @@ -11659,23 +11855,31 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[e]^1]++; (*pz) = int_a[e]^1; (*pz) <<= 32; (*pz) |= i; } + // prt_thread_info(b64.a, a_n, NULL, 0, "tt2"); // fprintf(stderr, "**1**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); ///index - radix_sort_srt64(b64.a + a_n, b64.a + b64.n); + radix_sort_srt64(b64.a + a_n, b64.a + b64.n);///(raw unitig node id)|(integer contig id) for (k = a_n + 1, i = a_n; k <= b64.n; k++) { if(k == b64.n || (b64.a[k]>>32) != (b64.a[i]>>32)) { - i_idx[b64.a[i]>>32] |= (((uint64_t)i)<<32)|((uint64_t)k); + i_idx[b64.a[i]>>32] |= (((uint64_t)i)<<32)|((uint64_t)k);///b64.a[i]>>32 appear once (unique ends)/multipe times i = k; } } + // prt_thread_info(b64.a, a_n, NULL, 0, "tt3"); // fprintf(stderr, "**2**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); - n_clus = usg_unique_arcs_cluster(&b64, a_n, i_idx, int_a); + ///b64.a[0, a_n]:: all resolvable paths with unique beg && end + ///b64.a[a_n, b64.n]:: (raw unitig/non-unqiue node id)|(resolvable path id) + // n_clus = usg_unique_arcs_cluster(&b64, a_n, i_idx, int_a);///this function might be wrong + n_clus = usg_unique_arcs_cluster_adv(&b64, a_n, i_idx, int_a, &ub64); + // fprintf(stderr, "**3**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus); + // prt_thread_info(b64.a, a_n, NULL, 0, "tt4"); assert(b64.n == (a_n<<1)); for (k = a_n + 1, i = a_n, mm = a_n; k <= b64.n; k++) { if(k == b64.n || (b64.a[k]>>32) != (b64.a[i]>>32)) { for (z = i; z < k; z++) {///all intger seqs within the same cluster s = b64.a[((uint32_t)b64.a[z])]>>32; e = ((uint32_t)b64.a[((uint32_t)b64.a[z])]); assert(e > s); + ///[s, e]:: available interval if(!ava_pass_unique_bridge(i_idx, int_a, s, e)) break; } if(z >= k) {///all arcs in this cluster is fine -> each of arch is reliable @@ -11685,6 +11889,7 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint } } assert(n_clus == 0); + // prt_thread_info(b64.a, a_n, NULL, 0, "tt5"); // fprintf(stderr, "**4**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus); b64.n = mm; n_clus = 0; for (k = a_n + 1, i = a_n, mm = a_n; k <= b64.n; k++) { @@ -11696,9 +11901,10 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint i = k; } } + // prt_thread_info(b64.a, a_n, NULL, 0, "tt6"); b64.n = mm; // fprintf(stderr, "**5**[M::%s::] a_n::%lu, n_clus::%lu, ng->n::%lu\n", __func__, a_n, n_clus, ng->n); - + // prt_thread_info(b64.a, a_n, b64.a + a_n, b64.n - a_n, "thred"); if(n_clus > 0) { update_usg_t_threading(uidx, ng, b64.a, b64.a + a_n, b64.n - a_n, int_a, ng_occ, &ub64); } @@ -11711,24 +11917,68 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint void u2g_hybrid_detan(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, asg64_v *in, asg64_v *ib) { - uint64_t k, i, x; asg64_v tx = {0,0,0}, tb = {0,0,0}, *ob = NULL, *ub = NULL; + uint64_t k, i, x, m, *tmp, sn; asg64_v tx = {0,0,0}, tb = {0,0,0}, *ob = NULL, *ub = NULL; ob = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); ob->n = ub->n = 0; for (k = 0; k < uidx->str_b.n_thread; k++) { uidx->str_b.buf[k].res_dump.n = uidx->str_b.buf[k].u.n = uidx->str_b.buf[k].o.n = 0; } kt_for(uidx->str_b.n_thread, worker_integer_realign_g, uidx, uidx->uovl.i_ug->u.n); - for (k = ob->n = ub->n = 0; k < uidx->str_b.n_thread; k++) { + + for (k = ob->n = ub->n = m = 0; k < uidx->str_b.n_thread; k++) { for (i = 0; i < uidx->str_b.buf[k].res_dump.n; i++) { - if((uidx->str_b.buf[k].res_dump.a[i]>>32)==((uint32_t)-1)) { - x = ob->n; x <<= 32; x |= ((uint32_t)uidx->str_b.buf[k].res_dump.a[i]); - kv_push(uint64_t, *ub, x); + x = uidx->str_b.buf[k].res_dump.a[i]; + kv_push(uint64_t, *ob, x);//aln details + if(x&((uint64_t)0x8000000000000000)) { + x -= ((uint64_t)0x8000000000000000); x >>= 32; x <<= 32;//seq_id + x |= ob->n;//offset + kv_push(uint64_t, *ub, x);///idx:: seq_id|offset_in_ob } else { - kv_push(uint64_t, *ob, uidx->str_b.buf[k].res_dump.a[i]); + m++; } - } + } } + radix_sort_srt64(ub->a, ub->a + ub->n); + kv_resize(uint64_t, *ub, ub->n+m); tmp = ub->a + ub->n; m = 0; + for (k = 0; k < ub->n; k++) { + sn = ((uint32_t)(ob->a[((uint32_t)ub->a[k])-1])); + memcpy(tmp + m, ob->a + ((uint32_t)ub->a[k]), sn*sizeof((*tmp))); + ub->a[k] = m; ub->a[k] <<= 32; ub->a[k] |= sn;//offset_in_ob|occ + m += sn; + } + assert(m <= ob->n); + memcpy(ob->a, tmp, m*sizeof((*tmp))); ob->n = m; + + + // for (k = 0, p = NULL; k < ub->n; k++) { + // p = &(ob->a[((uint32_t)ub->a[k])-1]); + // x = ub->a[k]<<32;//offset + // x |= ((uint32_t)(*p));///occ + // ub->a[k] = x; + // (*p) >>= 32; (*p) <<= 32; (*p) |= k; + // } + // for (k = m = 0; k < ob->n; k++) { + // if(ob->a[k]&((uint64_t)0x8000000000000000)) { + // x = m; x <<= 32; x |= ((uint32_t)ub->a[(uint32_t)ob->a[k]]); + // ub->a[(uint32_t)ob->a[k]] = x; + // } else { + // ob->a[m++] = ob->a[k]; + // } + // } + // ob->n = m; + + // for (k = ob->n = ub->n = 0; k < uidx->str_b.n_thread; k++) { + // for (i = 0; i < uidx->str_b.buf[k].res_dump.n; i++) { + // if((uidx->str_b.buf[k].res_dump.a[i]>>32)==((uint32_t)-1)) { + // x = ob->n; x <<= 32; x |= ((uint32_t)uidx->str_b.buf[k].res_dump.a[i]); + // kv_push(uint64_t, *ub, x);///idx:: offset_in_ob|occ + // } else { + // kv_push(uint64_t, *ob, uidx->str_b.buf[k].res_dump.a[i]);//aln details + // } + // } + // } + /*********debugging*********/ fprintf(stderr, "[M::%s::] # iug::%u, # gchain::%u\n", __func__, (uint32_t)uidx->uovl.i_ug->u.n, (uint32_t)ub->n); // for (k = 0; k < ub->n; k++) { @@ -11747,6 +11997,7 @@ void u2g_hybrid_detan(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, asg64_v * // u_ug = gen_unique_g(uidx, ng, ng_occ, ub->a, ub->n, ob->a);//ma_ug_hybrid_gen(ng); ///debug debug_sysm_usg_t(ng, __func__); + // prt_usg_t(uidx, ng, "ng4"); if(gen_unique_g_adv(uidx, ng, ub->a, ub->n, ob->a, max_ext)) { // usg_arc_t *z = get_usg_arc(ng, 2, 576), *q = get_usg_arc(ng, 577, 3); // fprintf(stderr, "xxxx0xxx[M::%s::] p->del::%u, q->del::%u\n", @@ -11762,6 +12013,7 @@ void u2g_hybrid_detan(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, asg64_v * debug_sysm_usg_t(ng, __func__); // fprintf(stderr, "-[M::%s::] ng->n::%u\n", __func__, (uint32_t)ng->n); } + // prt_usg_t(uidx, ng, "ng_dbg"); if(!in) free(tx.a); if(!ib) free(tb.a); } @@ -11825,31 +12077,65 @@ void u2g_hybrid_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, usg_t *ng, asg64_v * ((ulopt->max_ovlp_drop_ratio-ulopt->min_ovlp_drop_ratio)/(ulopt->clean_round-1))); double drop = ulopt->min_ovlp_drop_ratio; ///CALLOC(iug->g->seq_vis, iug->g->n_seq*2); fprintf(stderr, "\n[M::%s::] Starting hybrid clean, mm_tip::%ld\n", __func__, mm_tip); - - + // prt_usg_t(uidx, ng, "ng0"); usg_arc_cut_tips(ng, mm_tip, 0, b); + // prt_usg_t(uidx, ng, "ng1"); + // char sb[1000]; for (ss = 1; ss <= 1/**6**/; ss++) { mm_tip = ulopt->max_tip_hifi*ss; + fprintf(stderr, "\n[M::%s::] ss::%ld, mm_tip::%ld, ulopt->clean_round::%ld\n", + __func__, ss, mm_tip, ulopt->clean_round); for (i = 0, drop = ulopt->min_ovlp_drop_ratio; i < ulopt->clean_round; i++, drop += step) { if(drop > ulopt->max_ovlp_drop_ratio) drop = ulopt->max_ovlp_drop_ratio; + // fprintf(stderr, "-0-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop); + // sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_a", ss, i, drop); + // prt_usg_t(uidx, ng, sb); usg_arc_cut_length(ng, b, ub, mm_tip>>1, drop, ulopt->is_trio, 1, NULL); + // fprintf(stderr, "-1-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop); + // sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_b", ss, i, drop); + // prt_usg_t(uidx, ng, sb); usg_arc_cut_tips(ng, mm_tip, 0, b); + // fprintf(stderr, "-2-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop); + // sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_c", ss, i, drop); + // prt_usg_t(uidx, ng, sb); usg_arc_cut_length(ng, b, ub, mm_tip, drop, ulopt->is_trio, 1, NULL); + // fprintf(stderr, "-3-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop); + // sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_d", ss, i, drop); + // prt_usg_t(uidx, ng, sb); usg_arc_cut_tips(ng, mm_tip, 1, b); + // sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_e", ss, i, drop); + // prt_usg_t(uidx, ng, sb); + // fprintf(stderr, "-4-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop); } drop = 1; + // fprintf(stderr, "-0-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop); + // sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_a", ss, i, drop); + // prt_usg_t(uidx, ng, sb); usg_arc_cut_length(ng, b, ub, mm_tip>>1, drop, ulopt->is_trio, 1, NULL); + // fprintf(stderr, "-1-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop); + // sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_b", ss, i, drop); + // prt_usg_t(uidx, ng, sb); usg_arc_cut_tips(ng, mm_tip, 0, b); + // fprintf(stderr, "-2-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop); + // sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_c", ss, i, drop); + // prt_usg_t(uidx, ng, sb); usg_arc_cut_length(ng, b, ub, mm_tip, drop, ulopt->is_trio, 1, NULL); + // fprintf(stderr, "-3-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop); + // sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_d", ss, i, drop); + // prt_usg_t(uidx, ng, sb); usg_arc_cut_tips(ng, mm_tip, 1, b); + // sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_e", ss, i, drop); + // prt_usg_t(uidx, ng, sb); + // fprintf(stderr, "-4-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop); } + // prt_usg_t(uidx, ng, "ng2"); ///debug debug_sysm_usg_t(ng, __func__); /******for debug******/ - // prt_usg_t(uidx, ng, "ng1"); + prt_usg_t(uidx, ng, "ng_dbg"); /******for debug******/ // u2g_hybrid_extend(ng, NULL, b, ub); @@ -11928,9 +12214,9 @@ ma_ug_t *gen_hybrid_ug(ul_resolve_t *uidx, usg_t *ng) ug->g->seq[i].c = PRIMARY_LABLE; u = &(ug->u.a[i]); if(u->m == 0) continue; - fprintf(stderr, "+[M::%s::] i::%u\n", __func__, i); + // fprintf(stderr, "+[M::%s::] i::%u\n", __func__, i); merge_hybrid_utg_content(u, uidx->l1_ug, uidx->sg, ng, &e); - fprintf(stderr, "-[M::%s::] i::%u\n", __func__, i); + // fprintf(stderr, "-[M::%s::] i::%u\n", __func__, i); ug->g->seq[i].len = u->len; } kv_destroy(e.a); @@ -12690,6 +12976,7 @@ double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_ print_raw_uls_seq(uidx, asm_opt.output_file_name); ul_re_correct(uidx, 3); init_ulg_opt_t(&uu, uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_tip<<1, b_mask_t, is_trio); + print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug0", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1); /**ul2ul_idx_t *u2o = **/gen_ul2ul(uidx, uopt, &uu, 0); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-3"); diff --git a/inter.cpp b/inter.cpp index 759e006..746c52e 100644 --- a/inter.cpp +++ b/inter.cpp @@ -5184,7 +5184,7 @@ const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, uint6 if(lj->qs >= li->qs) return INT32_MIN; uint32_t li_v = (li->tn<<1)|li->rev, lj_v = (lj->tn<<1)|lj->rev; int64_t qo = infer_rovlp(li, lj, NULL, NULL, ridx, ug), trans_l = 0, sec_err = 0, sc; ///overlap length in query (UL read) - if(li_v != lj_v && get_ecov_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, mode, NULL)) { + if(/**li_v != lj_v &&**/ get_ecov_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, mode, NULL)) { trans_l = get_overlap_region_sub_err(&(ol[li->qn]), tc, lj->qe, &sec_err); // int64_t trans_l_debug, sec_err_debug; @@ -5445,7 +5445,7 @@ int64_t *t, int64_t n_ext) if(lj->qs >= li->qs) continue; set_ul_ov_t_by_mg_lchain_t(&uj, lj); qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug); ///overlap length in query (UL read) - if(li->v!=lj->v && get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, ng_diff_thre, qo, 0, NULL)) { + if(/**li->v!=lj->v &&**/get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, ng_diff_thre, qo, 0, NULL)) { sc = csc + f[j]; if(sc > mm_sc) { mm_sc = sc, mm_idx = j; @@ -5473,7 +5473,7 @@ int64_t *t, int64_t n_ext) if(lj->qe+G_CHAIN_INDEL > li->qs && lj->qs < li->qs) { set_ul_ov_t_by_mg_lchain_t(&uj, lj); qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug);///overlap length in query (UL read) - if(li->v!=lj->v && get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, ng_diff_thre, qo, 0, NULL)) { + if(/**li->v!=lj->v &&**/get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, ng_diff_thre, qo, 0, NULL)) { sc = csc + f[max_ii]; if(sc > mm_sc) { mm_sc = sc; mm_idx = max_ii; @@ -5651,7 +5651,7 @@ st_mt_t *bf, Chain_Data* dp, int64_t max_skip, int64_t max_iter, int64_t max_dis if((!is_f) && (lj->qe+G_CHAIN_INDEL > li->qs)) { set_ul_ov_t_by_mg_lchain_t(&uj, lj); qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug); - if(li->v!=lj->v && get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, N_GCHAIN_RATE, qo, 0, NULL)) { + if(/**li->v!=lj->v &&**/ get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, N_GCHAIN_RATE, qo, 0, NULL)) { is_f = 1; if(n_skip > 0) n_skip--; if(n_skip < (max_skip>>1)) n_skip= (max_skip>>1); } @@ -5882,7 +5882,7 @@ int64_t need_srt) if((!is_f) && (lj->qe+G_CHAIN_INDEL > li->qs)) { set_ul_ov_t_by_mg_lchain_t(&uj, lj); qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug); - if(li->v!=lj->v && get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, N_GCHAIN_RATE, qo, 0, NULL)) { + if(/**li->v!=lj->v &&**/ get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, N_GCHAIN_RATE, qo, 0, NULL)) { is_f = 1; if(n_skip > 0) n_skip--; if(n_skip < (max_skip>>1)) n_skip= (max_skip>>1); } @@ -6095,9 +6095,12 @@ const asg_t *g, st_mt_t *dst_done, vec_sp_node_t *out, vec_mg_pathv_t *res, vec_ if(res->a[i].v == (uint32_t)-1) { gchains->a[i+gchains->n] = a[res->a[i].pre + g_item->rs]; gchains->a[i+gchains->n].dist_pre = res->a[i].d; - - // fprintf(stderr, "+[M::%s::]\tutg%.6dl(%c)\n", __func__, - // (int32_t)(gchains->a[i+gchains->n].v>>1)+1, "+-"[gchains->a[i+gchains->n].v&1]); + // if(ulid == 14714) { + // fprintf(stderr, "+[M::%s::]\tutg%.6dl(%c)\tq::[%d, %d)\tt::[%d, %d)\n", __func__, + // (int32_t)(gchains->a[i+gchains->n].v>>1)+1, "+-"[gchains->a[i+gchains->n].v&1], + // gchains->a[i+gchains->n].qs, gchains->a[i+gchains->n].qe, + // gchains->a[i+gchains->n].rs, gchains->a[i+gchains->n].re); + // } } else { gchains->a[i+gchains->n].v = res->a[i].v; gchains->a[i+gchains->n].off = -1; @@ -6108,6 +6111,12 @@ const asg_t *g, st_mt_t *dst_done, vec_sp_node_t *out, vec_mg_pathv_t *res, vec_ // fprintf(stderr, "aaaaaaa, ulid->%ld\n", ulid); // fprintf(stderr, "-[M::%s::]\tutg%.6dl(%c)\n", __func__, // (int32_t)(gchains->a[i+gchains->n].v>>1)+1, "+-"[gchains->a[i+gchains->n].v&1]); + // if(ulid == 14714) { + // fprintf(stderr, "-[M::%s::]\tutg%.6dl(%c)\tq::[%d, %d)\tt::[%d, %d)\n", __func__, + // (int32_t)(gchains->a[i+gchains->n].v>>1)+1, "+-"[gchains->a[i+gchains->n].v&1], + // gchains->a[i+gchains->n].qs, gchains->a[i+gchains->n].qe, + // gchains->a[i+gchains->n].rs, gchains->a[i+gchains->n].re); + // } } } g_item->cnt = res->n; @@ -8803,6 +8812,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // if(s->id+i!=3046/** && s->id+i!=3111**/) return; // if((s->id+i!=871) && (s->id+i!=963) && (s->id+i!=980)) return; // if(s->id+i!=963) return; + // if(s->id+i != 35437) 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); @@ -9603,7 +9613,7 @@ utg_rid_dt *get_r_ug_region(utg_rid_t *idx, uint64_t *n, uint64_t rid) return (*n)?idx->p.a + idx->idx[rid]:NULL; } -void rov2uov(uint64_t rid, const ul_idx_t *uref, utg_rid_dt *ru_map, uc_block_t *rovlp, ul_ov_t *res, uint32_t adjust_rev) +uint64_t rov2uov(uint64_t rid, const ul_idx_t *uref, utg_rid_dt *ru_map, uc_block_t *rovlp, ul_ov_t *res, uint32_t adjust_rev, int64_t ulid) { uint64_t ori = ru_map->u&1, ts, te; if(!ori) { @@ -9612,15 +9622,27 @@ void rov2uov(uint64_t rid, const ul_idx_t *uref, utg_rid_dt *ru_map, uc_block_t ts = uref->r_ug->rg->seq[rid].len - rovlp->te; te = uref->r_ug->rg->seq[rid].len - rovlp->ts; } + // if(ulid == 14714) { + // fprintf(stderr, "[M::%s::]\tori::%lu\trovlp->rev::%u\tro_t::[%u, %u)\tt::[%lu, %lu)\toff::%u\n", __func__, + // ori, rovlp->rev, rovlp->ts, rovlp->te, ts, te, ru_map->off); + // } ts += ru_map->off; te += ru_map->off; - memset(res, 0, sizeof(*res)); - res->qn = 0; res->qs = rovlp->qs; res->qe = rovlp->qe; - res->tn = ru_map->u>>1; res->ts = ts; res->te = te; - res->el = rovlp->el; res->rev = (rovlp->rev == ori?0:1); - if(adjust_rev && res->rev) {///for linear chaining - res->ts = uref->ug->g->seq[res->tn].len - te; - res->te = uref->ug->g->seq[res->tn].len - ts; - } + if(ts >= 0 && te <= uref->ug->g->seq[ru_map->u>>1].len) { + memset(res, 0, sizeof(*res)); + res->qn = 0; res->qs = rovlp->qs; res->qe = rovlp->qe; + res->tn = ru_map->u>>1; res->ts = ts; res->te = te; + res->el = rovlp->el; res->rev = (rovlp->rev == ori?0:1); + if(adjust_rev && res->rev) {///for linear chaining + res->ts = uref->ug->g->seq[res->tn].len - te; + res->te = uref->ug->g->seq[res->tn].len - ts; + } + return 1; + } + return 0; + // if(ulid == 14714) { + // fprintf(stderr, "[M::%s::]\tulen::%u\trlen::%u\tro_t::[%u, %u)\tt::[%lu, %lu)\toff::%u\n", __func__, + // uref->ug->g->seq[res->tn].len, uref->r_ug->rg->seq[rid].len, rovlp->ts, rovlp->te, ts, te, ru_map->off); + // } } void print_ul_ov_t(ul_ov_t *xs, const char* cmd) @@ -9629,9 +9651,9 @@ void print_ul_ov_t(ul_ov_t *xs, const char* cmd) "+-"[xs->rev], (int)Get_NAME_LENGTH(R_INF, ((xs->tn<<1)>>1)), Get_NAME(R_INF, ((xs->tn<<1)>>1)), xs->ts, xs->te); } -void gl_rg2ug_gen(ul_vec_t *r_cl, kv_ul_ov_t *u_cl, const ul_idx_t *uref, uint64_t is_el, uint64_t n_pchain) +void gl_rg2ug_gen(ul_vec_t *r_cl, kv_ul_ov_t *u_cl, const ul_idx_t *uref, uint64_t is_el, uint64_t n_pchain, int64_t ulid) { - uint64_t k, a_k, a_n; uc_block_t *z; utg_rid_dt *a; ul_ov_t *p; + uint64_t k, a_k, a_n; uc_block_t *z; utg_rid_dt *a; ul_ov_t p; u_cl->n = 0; for (k = 0; k < r_cl->bb.n; k++) { z = &(r_cl->bb.a[k]); @@ -9641,14 +9663,22 @@ void gl_rg2ug_gen(ul_vec_t *r_cl, kv_ul_ov_t *u_cl, const ul_idx_t *uref, uint64 a = get_r_ug_region(uref->r_ug, &a_n, z->hid); if(!a) continue; for (a_k = 0; a_k < a_n; a_k++) { - kv_pushp(ul_ov_t, *u_cl, &p); - // fprintf(stderr, "\n+[M::%s::] %u\t%u\t%c\t%.*s(%u)\t%u\t%u\n", __func__, z->qs, z->qe, "+-"[z->rev], - // (int)Get_NAME_LENGTH(R_INF, z->hid), Get_NAME(R_INF, z->hid), (uint32_t)Get_READ_LENGTH(R_INF, z->hid), z->ts, z->te); - // fprintf(stderr, "*[M::%s::] utg%.6d%c(%u)\t%c\t%u\n", __func__, - // (int32_t)(a[a_k].u>>1)+1, "lc"[uref->ug->u.a[a[a_k].u>>1].circ], uref->ug->u.a[a[a_k].u>>1].len, - // "+-"[a[a_k].u&1], a[a_k].off); - rov2uov(z->hid, uref, &(a[a_k]), z, p, 1); - p->el = 1; p->tn <<= 1; p->tn |= p->rev; p->qn = k/**uref->r_ug->idx[z->hid] + a_k**/;//for linear chain + // if(ulid == 14714) { + // fprintf(stderr, "\n+[M::%s::]\tq::[%u, %u)\t%c\t%.*s(%u)\tt::[%u, %u)\n", __func__, z->qs, z->qe, "+-"[z->rev], + // (int)Get_NAME_LENGTH(R_INF, z->hid), Get_NAME(R_INF, z->hid), (uint32_t)Get_READ_LENGTH(R_INF, z->hid), + // z->ts, z->te); + // fprintf(stderr, "*[M::%s::] utg%.6d%c(%u)\t%c\toff::%u\tpos::%u\n", __func__, + // (int32_t)(a[a_k].u>>1)+1, "lc"[uref->ug->u.a[a[a_k].u>>1].circ], uref->ug->u.a[a[a_k].u>>1].len, + // "+-"[a[a_k].u&1], a[a_k].off, a[a_k].pos); + // } + if(!rov2uov(z->hid, uref, &(a[a_k]), z, &p, 1, ulid)) continue; + p.el = 1; p.tn <<= 1; p.tn |= p.rev; p.qn = k/**uref->r_ug->idx[z->hid] + a_k**/;//for linear chain + kv_push(ul_ov_t, *u_cl, p); + // if(ulid == 14714) { + // fprintf(stderr, "[M::%s::]\tutg%.6dl(%c)\tq::[%d, %d)\tt::[%d, %d)\ttlen::%u\n", __func__, + // (int32_t)(p->tn>>1)+1, "+-"[p->tn&1], p->qs, p->qe, + // p->ts, p->te, uref->ug->g->seq[p->tn>>1].len); + // } // if(k == 2) { // fprintf(stderr, "[M::%s::] p->ts:%u, p->te:%u, z->ts:%u, z->te:%u, a[a_k].off:%u\n", __func__, p->ts, p->te, z->ts, z->te, a[a_k].off); // } @@ -10142,7 +10172,7 @@ void extend_end_coord(mg_lchain_t *li, ul_ov_t *ui, const int64_t qlen, const in } } -void dump_linear_chain(ma_ug_t *ug, kv_ul_ov_t *lidx, vec_mg_lchain_t *res, int64_t qlen) +void dump_linear_chain(ma_ug_t *ug, kv_ul_ov_t *lidx, vec_mg_lchain_t *res, int64_t qlen, int64_t ulid) { uint64_t i; int64_t iqs, iqe, its, ite; mg_lchain_t *p; kv_resize(mg_lchain_t, *res, lidx->n); @@ -10154,12 +10184,21 @@ void dump_linear_chain(ma_ug_t *ug, kv_ul_ov_t *lidx, vec_mg_lchain_t *res, int6 p->off = i; p->score = lidx->a[i].sec; p->qs = lidx->a[i].qs; p->qe = lidx->a[i].qe; p->rs = lidx->a[i].ts; p->re = lidx->a[i].te; + // if(ulid == 14714) { + // fprintf(stderr, "+[M::%s::]\tutg%.6dl(%c)\tq::[%d, %d)\tqlen::%ld\tt::[%d, %d)\ttlen::%u\n", __func__, + // (int32_t)(p->v>>1)+1, "+-"[p->v&1], p->qs, p->qe, qlen, + // p->rs, p->re, ug->g->seq[p->v>>1].len); + // } extend_end_coord(p, NULL, qlen, ug->g->seq[p->v>>1].len, &iqs, &iqe, &its, &ite); p->qs = iqs; p->qe = iqe; p->rs = its; p->re = ite; // if(!ugl_cover_check(p->rs, p->re, &(ug->u.a[p->v>>1]))) res->n--; // fprintf(stderr, "chain_id:%d\t%u\t%u\t%c\tutg%.6dl(%u)\t%u\t%u\n", // res->a[k].off, res->a[k].qs, res->a[k].qe, "+-"[res->a[k].v&1], (int32_t)(res->a[k].v>>1)+1, // g->seq[res->a[k].v>>1].len, res->a[k].rs, res->a[k].re); + // if(ulid == 14714) { + // fprintf(stderr, "-[M::%s::]\tutg%.6dl(%c)\tq::[%d, %d)\tt::[%d, %d)\n", __func__, + // (int32_t)(p->v>>1)+1, "+-"[p->v&1], p->qs, p->qe, p->rs, p->re); + // } } } @@ -11670,7 +11709,7 @@ mg_lchain_t *uo, kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn) } void gl_ug2rg_gen(const asg_t *rg, ul_vec_t *rch, ma_ug_t *ug, mg_lchain_t *uo, vec_mg_lchain_t *res, int64_t tOff, -kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn) +int64_t ulid, kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn) { // fprintf(stderr, "\n[M::%s::] uo->qs:%d, uo->qe:%d\n", __func__, uo->qs, uo->qe); ///uo is a unitig alignment @@ -11684,6 +11723,10 @@ kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn) if(p.s >= re) break; assert(extract_rovlp_by_ug(&p, uo, res, tOff)); + // if(!extract_rovlp_by_ug(&p, uo, res, tOff)) { + // fprintf(stderr, "[M::%s::ulid::%ld]u_rs->%lu, u_re->%lu, p.s->%u, p.e->%u\n", __func__, ulid, rs, re, p.s, p.e); + // exit(1); + // } res->a[res->n-1].score = uo->v; res->a[res->n-1].cnt = i; // fprintf(stderr, "[M::%s::]u_rs->%lu, u_re->%lu, p.s->%u, p.e->%u\n", __func__, rs, re, p.s, p.e); @@ -11972,7 +12015,7 @@ int64_t flat_rovlp_chain(ma_ug_t *ug, mg_lchain_t *x, int64_t x_n, Chain_Data* d } void gen_rovlp_chain_by_ul(const asg_t *rg, ul_vec_t *rch, const ul_idx_t *uref, kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn, mg_lchain_t *a, int64_t a_n, vec_mg_lchain_t *res, Chain_Data* dp, -int64_t dp_max_skip, int64_t dp_max_iter, int64_t dp_max_dis) +int64_t dp_max_skip, int64_t dp_max_iter, int64_t dp_max_dis, int64_t ulid) { if(a_n == 0) return; int64_t k, l, res_n0 = res->n, tt = 0; ma_ug_t *ug = uref->ug; @@ -11981,7 +12024,7 @@ int64_t dp_max_skip, int64_t dp_max_iter, int64_t dp_max_dis) for (k = 0, l = ug->g->seq[a[0].v>>1].len; k < a_n; k++) { // fprintf(stderr, ">k->%ld, ls->%ld, le->%ld, rev->%c\n", k, l - ug->g->seq[a[k].v>>1].len, l, "+-"[a[k].v&1]); l -= ug->g->seq[a[k].v>>1].len; - gl_ug2rg_gen(rg, rch, ug, &(a[k]), res, l, raw_idx, raw_chn); + gl_ug2rg_gen(rg, rch, ug, &(a[k]), res, l, ulid, raw_idx, raw_chn); l += ug->g->seq[a[k].v>>1].len + a[k].dist_pre; } mg_lchain_t *x = res->a + res_n0; int64_t x_n = res->n - res_n0, fn; @@ -12296,10 +12339,12 @@ int64_t dp_max_skip, int64_t dp_max_iter, int64_t dp_max_dis) int64_t k, ucn = uc->n; mg_lchain_t *ix; for (k = 0, swap->n = 0; k < ucn; k += ix->cnt + 1) { ix = &(uc->a[k]); assert(ix->v == (uint32_t)-1); - // fprintf(stderr, "\n[M::%s::ucn->%ld, k->%ld, kcnt->%d]\n", __func__, ucn, k, ix->cnt); - // print_debug_gchain(uref, uc->a + k + 1, ix->cnt, rch); + // if(ulid == 14714) { + // fprintf(stderr, "\n[M::%s::ucn->%ld, k->%ld, kcnt->%d]\n", __func__, ucn, k, ix->cnt); + // print_debug_gchain(uref, uc->a + k + 1, ix->cnt, rch); + // } gen_rovlp_chain_by_ul(rg, rch, uref, raw_idx, raw_chn, uc->a + k + 1, ix->cnt, swap, dp, - dp_max_skip, dp_max_iter, dp_max_dis); + dp_max_skip, dp_max_iter, dp_max_dis, ulid); } ///up to now, given a in swap @@ -12345,7 +12390,7 @@ int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid, Chain_Data* dp, c // if(ulid != 86660) return 0; kv_ul_ov_t *idx = &(ll->lo), *init = &(ll->tk); int64_t max_idx; idx->n = init->n = 0; - gl_rg2ug_gen(rch, idx, uref, 1, 2); + gl_rg2ug_gen(rch, idx, uref, 1, 2, ulid); if(idx->n == 0) return 0; ///generate linear chains gen_linear_chains(idx, init, uref, uopt, bw, diff_ec_ul, rch->rlen, dp); @@ -12358,7 +12403,7 @@ int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid, Chain_Data* dp, c // fprintf(stderr, "\n++[M::%s::%.*s(id:%ld), len:%u] idx->n:%lu\n", __func__, UL_INF.nid.a[ulid].n, UL_INF.nid.a[ulid].a, // ulid, rch->rlen, (uint64_t)idx->n); - dump_linear_chain(uref->ug, idx, &(gdp->l), rch->rlen); + dump_linear_chain(uref->ug, idx, &(gdp->l), rch->rlen, ulid); if(gdp->l.n == 0) return 0; // fprintf(stderr, "\n+++[M::%s::id->%ld, len->%u] idx->n:%lu\n", __func__, ulid, rch->rlen, (uint64_t)idx->n); // kv_resize(uint64_t, ll->srt.a, idx->n); kv_resize(uint64_t, hap->snp_srt, idx->n); kv_resize(uint64_t, gdp->v, idx->n); @@ -14645,7 +14690,7 @@ ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_c ///for debug interval if(!load_all_ul_t(&UL_INF, gfa_name, &R_INF, ug)) { gen_UL_reovlps(&sl, ug, sg, gfa_name, cutoff); - // exit(1); + exit(1); write_all_ul_t(&UL_INF, gfa_name, ug); } else if(double_check_cache){ if(drenew_UL_reovlps(&sl, ug, sg, gfa_name, cutoff)) { diff --git a/inter.h b/inter.h index b361c1d..1684581 100644 --- a/inter.h +++ b/inter.h @@ -13,7 +13,7 @@ #define G_CHAIN_GAP 0.1 #define UG_SKIP 5 #define RG_SKIP 25 -#define UG_SKIP_GRAPH_N 50 +#define UG_SKIP_GRAPH_N 72 #define UG_SKIP_N 100 #define UG_ITER_N 5000 #define UG_DIS_N 50000