diff --git a/Correct.cpp b/Correct.cpp index 0b73d00..70f0053 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -13,6 +13,8 @@ #include "kalloc.h" #include "htab.h" #include "Overlaps.h" +#define A_L 16 +#define ext_w 6 #define generic_key(x) (x) KRADIX_SORT_INIT(b32, uint32_t, generic_key, 4) @@ -1291,7 +1293,7 @@ void append_unmatched_wins(overlap_region *z, int64_t block_s) z->w_list.a[m].cidx = z->w_list.a[m].clen = 0; z->w_list.a[m].y_start = z->w_list.a[m].y_end = -1; z->w_list.a[m].error = z->w_list.a[m].error_threshold = -1; - z->w_list.a[m].extra_begin = z->w_list.a[m].extra_end = -1; + z->w_list.a[m].extra_begin = z->w_list.a[m].extra_end = 0; z->w_list.a[m].x_start = w_e + 1; wn_id = get_win_id_by_s(z, z->w_list.a[m].x_start, block_s, &w_e); z->w_list.a[m].x_end = w_e; @@ -12651,8 +12653,93 @@ void adjust_ext_offset(int64_t *qs, int64_t *qe, int64_t *ts, int64_t *te, int64 if((*te) > tl) (*te) = tl; } +void adjust_specific_ext_offset(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, bit_extz_t *exz, char* qstr, UC_Read *tu, +int64_t *qs0, int64_t *qe0, int64_t *ts0, int64_t *te0, int64_t ql, int64_t tl, int64_t wl, int64_t *mode) +{ + ///[qs, qe) + if((*mode) != 1 && (*mode) != 2) return; + int64_t qs = *qs0, qe = *qe0, ts = *ts0, te = *te0, k, update, wid, we, ws, sl, gq, gt, gg, min_g = INT32_MAX, min_id = -1, gl; + adjust_ext_offset(&qs, &qe, &ts, &te, ql, tl, 0, *mode); + if((*mode) == 1) {///forward extension + //qs0 and ts0 are fixed; te0 = -1, qe0 is unreliable + if(qe > (*qe0)) { + we = (*qe0); we/=wl; we *= wl; we +=wl; we--; ///next window + for (k = 0; we < qe && k < ext_w; we+=wl) {//[ws, we]; [qs, qe) + wid = get_win_id_by_e(z, we, wl, NULL); + if(z->w_list.a[wid].y_end == -1) continue;//unmapped + if(z->w_list.a[wid].x_end < (*qs0)) continue; + if(z->w_list.a[wid].y_end < (*ts0)) continue; + sl = z->w_list.a[wid].x_end+1-z->w_list.a[wid].x_start; + if(z->w_list.a[wid].error > (sl/A_L)) continue; + gq = z->w_list.a[wid].x_end-(*qs0); + gt = z->w_list.a[wid].y_end-(*ts0); + gl = MIN(gq, gt); gl/=16; if(gl <= 0) gl = 1; + gg = (gq>=gt)?(gq-gt):(gt-gq); gg /= gl; + update = 0; + if(min_g>gg) { + update = 1; + } else if((min_g==gg)&&(z->w_list.a[min_id].error>z->w_list.a[wid].error)) { + update = 1; + } + if(update) { + min_g = gg; min_id = wid; + } + k++; + } + } + if(min_id >= 0) { + (*qe0) = z->w_list.a[min_id].x_end+1; + if(((*qe0)+wl) >= qe) { + (*qe0) = qe; (*te0) = te; ///still extension + return; + } + (*te0) = z->w_list.a[min_id].y_end+1; (*mode) = 0;///global + } else { + (*qe0) = qe; (*te0) = te; ///still extension + } + } else if((*mode) == 2) {///backward extension + //qe0 and te0 are fixed; ts0 = -1, qs0 is unreliable + if(qs < (*qs0)) { + ws = (*qs0)-1; ws/=wl; ws*=wl; + for (k = 0; ws >= qs && k < ext_w; ws-=wl) {//[ws, we]; [qs, qe) + wid = get_win_id_by_s(z, ws, wl, NULL); + if(z->w_list.a[wid].y_end == -1) continue;//unmapped + if(z->w_list.a[wid].x_start >= (*qe0)) continue; + if(z->w_list.a[wid].y_start >= (*te0)) continue; + sl = z->w_list.a[wid].x_end+1-z->w_list.a[wid].x_start; + if(z->w_list.a[wid].error > (sl/A_L)) continue; + gq = (*qe0) - z->w_list.a[wid].x_start; + gt = (*te0) - z->w_list.a[wid].y_start; + gl = MIN(gq, gt); gl/=16; if(gl <= 0) gl = 1; + gg = (gq>=gt)?(gq-gt):(gt-gq); gg /= gl; + update = 0; + if(min_g>gg) { + update = 1; + } else if((min_g==gg)&&(z->w_list.a[min_id].error>z->w_list.a[wid].error)) { + update = 1; + } + if(update) { + min_g = gg; min_id = wid; + } + k++; + } + } + if(min_id >= 0) { + (*qs0) = z->w_list.a[min_id].x_start; + if((qs+wl) >= (*qs0)) { + (*qs0) = qs; (*ts0) = ts; ///still extension + return; + } + (*ts0) = z->w_list.a[min_id].y_start; (*mode) = 0;///global + } else { + (*qs0) = qs; (*ts0) = ts; ///still extension + } + } +} + int64_t cal_exz_infi(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, bit_extz_t *exz, char* qstr, UC_Read *tu, int64_t qs, int64_t qe, int64_t ts, int64_t te, int64_t thre, int64_t q_tot_l, int64_t mode) { + clear_align(*exz); int64_t aux_beg = 0, bd = (((thre)<<1)+1), ql, tl, t_tot_l = -1; int32_t nword = ((bd>>bitw)+(!!(bd&bitz))); char *q_string, *t_string; int32_t rev = z->y_pos_strand, id = z->y_id; ql = qe - qs; if(hpc_g) t_tot_l = hpc_len(*hpc_g, id); @@ -12674,7 +12761,6 @@ int64_t cal_exz_infi(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_ } else { t_string = return_str_seq_exz(tu->seq, ts, tl, rev, hpc_g, uref, id); } - clear_align(*exz); // fprintf(stderr, ", nword::%d", nword); // if(ql < 0) fprintf(stderr, "qs::%ld, qe::%ld\n", qs, qe); // return 0; @@ -12719,15 +12805,29 @@ int64_t cal_exz_infi(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_ return 0; } -void hc_aln_exz(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, +void hc_aln_exz(overlap_region *z, Candidates_list *cl, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, int64_t qs, int64_t qe, int64_t ts, int64_t te, int64_t estimate_err, int64_t mode, int64_t wl, bit_extz_t *exz, int64_t q_tot, double e_rate) { - int64_t thre, ql = qe - qs, thre0; + clear_align(*exz); + int64_t thre, ql = qe - qs, thre0, t_tot, id = z->y_id; if(((ts == -1) && (te == -1))) mode = 3;///set to semi-global + + if(mode == 1 || mode == 2) { + if(hpc_g) t_tot = hpc_len(*hpc_g, id); + else if(uref) t_tot = uref->ug->u.a[id].len; + else t_tot = Get_READ_LENGTH((*rref), id); + // fprintf(stderr, "\n+[M::%s::ql->%ld::tl->%ld::mode->%ld] q::[%ld, %ld), t::[%ld, %ld), zx[%d, %d), zy[%d, %d)\n", + // __func__, q_tot, t_tot, mode, qs, qe, ts, te, z->x_pos_s, z->x_pos_e+1, z->y_pos_s, z->y_pos_e+1); + ///find a aligned window >= qe to convert extension to global; or do extension directly + adjust_specific_ext_offset(z, uref, hpc_g, rref, exz, qstr, tu, &qs, &qe, &ts, &te, q_tot, t_tot, wl, &mode); + // fprintf(stderr, "-[M::%s::ql->%ld::tl->%ld::mode->%ld] q::[%ld, %ld), t::[%ld, %ld), zx[%d, %d), zy[%d, %d)\n", + // __func__, q_tot, t_tot, mode, qs, qe, ts, te, z->x_pos_s, z->x_pos_e+1, z->y_pos_s, z->y_pos_e+1); + } + // fprintf(stderr, "[M::%s::ql::%ld] qs::%ld, qe::%ld, ts::%ld, te::%ld, mode::%ld, estimate_err::%ld, e_rate::%f", // __func__, ql, qs, qe, ts, te, mode, estimate_err, e_rate); - + if(ql <= MAX_L && (estimate_err*1.2) <= MAX_E) { thre = scale_ed_thre(estimate_err); if(thre > ql) thre = ql; if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, q_tot, mode)) { @@ -12763,31 +12863,36 @@ int64_t mode, int64_t wl, bit_extz_t *exz, int64_t q_tot, double e_rate) } } // fprintf(stderr, ", err::%d, thre::%d\n", INT32_MAX, exz->thre); + // anchor_aln(z, cl, uref, hpc_g, rref, qstr, tu, qs, qe, ts, te, thre, mode, wl, exz, q_tot, A_L); } -void sub_ciagar_gen(overlap_region *z, uint64_t s, uint64_t e, uint64_t wl, +void sub_ciagar_gen(overlap_region *z, Candidates_list *cl, uint64_t s, uint64_t e, uint64_t wl, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate, int64_t ql, uint64_t rid) { - uint64_t qs, qe, sid, eid, k, l, m, tot_e, c_e; int64_t q[2], t[2], mode; + uint64_t qs, qe, sid, eid, k, l, m, tot_e, c_e; int64_t q[2], t[2], o[2], mode, kocc; qs = (s/wl)*wl; if(qs < z->x_pos_s) qs = z->x_pos_s; if(qs > z->x_pos_e) return; qe = (e/wl)*wl; if(qe < e) qe += wl; if(qe > z->x_pos_e+1) qe = z->x_pos_e+1; if(qe <= 0) return; sid = get_win_id_by_s(z, qs, wl, NULL); eid = get_win_id_by_e(z, qe-1, wl, NULL) + 1;///must qe-1 instead of qe!!!!!! if(sid >= eid) return; - // fprintf(stderr, "***[M::%s] s::%lu, e::%lu, n_qs::%lu, n_qe::%lu, z::[%u, %u), sid::%lu, eid::%lu, w_list.n::%lu\n", - // __func__, s, e, qs, qe, z->x_pos_s, z->x_pos_e+1, sid, eid, (uint64_t)z->w_list.n); + // fprintf(stderr, "\n***[M::%s::rid->%lu] s::%lu, e::%lu, n_qs::%lu, n_qe::%lu, z::[%u, %u), sid::%lu, eid::%lu, w_list.n::%lu\n", + // __func__, rid, s, e, qs, qe, z->x_pos_s, z->x_pos_e+1, sid, eid, (uint64_t)z->w_list.n); for (k = sid+1, l = sid; k <= eid; k++) {//[sid, eid) if(k == eid || z->w_list.a[k].extra_end < 0) { if(k - l > 1 || z->w_list.a[l].extra_end >= 0) { - q[0] = q[1] = t[0] = t[1] = -1; mode = -1; tot_e = 0; + q[0] = q[1] = t[0] = t[1] = -1; + o[0] = o[1] = -1; mode = -1; tot_e = 0; if(z->w_list.a[l].extra_end < 0) { q[0] = z->w_list.a[l].x_end+1; if(z->w_list.a[l].y_end != -1) { t[0] = z->w_list.a[l].y_end+1; } + if(z->w_list.a[l].extra_end != INT16_MIN) { + o[0] = -z->w_list.a[l].extra_end; + } } else {///first window q[0] = qs; if(z->w_list.a[l].y_end != -1) { @@ -12809,7 +12914,10 @@ int64_t ql, uint64_t rid) // } assert(z->w_list.a[k-1].y_end != -1); t[1] = z->w_list.a[k-1].y_end+1; - } + } + if(z->w_list.a[k].extra_end != INT16_MIN) { + o[1] = -z->w_list.a[k].extra_end; + } } else {///last window q[1] = qe; } @@ -12837,7 +12945,8 @@ int64_t ql, uint64_t rid) // fprintf(stderr, "[M::%s::ql::%lu] qs::%lu, qe::%lu, ts::%lu, te::%lu, mode::%ld, tot_e::%lu\n", // __func__, q[1]-q[0], q[0], q[1], t[0], t[1], mode, tot_e); // } - hc_aln_exz(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], tot_e, mode, wl, exz, ql, e_rate); + kocc = MAX(o[0], o[1]); if(kocc < 0) kocc = 1; + hc_aln_exz(z, cl, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], tot_e, mode, wl, exz, ql, e_rate); } l = k; @@ -12845,37 +12954,149 @@ int64_t ql, uint64_t rid) } } -void cigar_gen(overlap_region *z, ul_ov_t *ov, uint64_t on, uint64_t qn, uint64_t wl, -const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate, int64_t ql, uint64_t rid) + +uint64_t cigar_gen(overlap_region *z, Candidates_list *cl, ul_ov_t *ov, uint64_t on, uint64_t qn, uint64_t wl, +const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate, +int64_t ql, uint64_t rid, ul_ov_t *des) { - uint64_t i, qs = (uint64_t)-1, qe = (uint64_t)-1; - for (i = 0; i < on && ov[i].qn == qn; i++) { - assert((i<=0)||(ov[i].qs >= ov[i-1].qe)); - if(ov[i].qs <= qe && qe != (uint64_t)-1) { - qe = ov[i].qe; - } else { - if(qs != (uint64_t)-1) sub_ciagar_gen(z, qs, qe, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); - qs = ov[i].qs; qe = ov[i].qe; + if(on <= 0) return 0; + uint64_t i; + // des[0] = ov[0]; + // for (i = m = 1; i < on; i++) { + // fusion_merge(); + // } + + + for (i = 0; i < on; i++) { + assert((i<=0)||(ov[i].qs > ov[i-1].qe)); + sub_ciagar_gen(z, cl, ov[i].qs, ov[i].qe, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); + } + return on; + // uint64_t i, qs = (uint64_t)-1, qe = (uint64_t)-1; + // for (i = 0; i < on && ov[i].qn == qn; i++) { + // assert((i<=0)||(ov[i].qs >= ov[i-1].qe)); + // if(ov[i].qs <= qe && qe != (uint64_t)-1) { + // qe = ov[i].qe; + // } else { + // if(qs != (uint64_t)-1) sub_ciagar_gen(z, cl, qs, qe, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); + // qs = ov[i].qs; qe = ov[i].qe; + // } + // } + // if(qs != (uint64_t)-1) sub_ciagar_gen(z, cl, qs, qe, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); +} + + + +int64_t push_adp_k_hits(Candidates_list *cl, uint64_t qs, uint64_t qe, uint64_t ts, uint64_t te, uint64_t readID, int64_t ci) +{ + int64_t k = ci, cln = cl->length; k_mer_hit *p = NULL, *m = NULL; + if(cln > 0) { + m = &(cl->list[cln-1]); if(!is_alnw(*m)) m = NULL; + } + for (; k < cln && cl->list[k].readID == readID && cl->list[k].self_offset < qs; k++) { + p = &(cl->list[k]); + if(p->offset >= ts) continue; + if((m) && ((m->offset >= p->offset) || (m->self_offset >= p->self_offset))) continue; + kv_pushp_cl(k_mer_hit, (*cl), &p); *p = cl->list[k]; p->readID = ((uint32_t)(0x7fffffff)); + } + + if(qs == (uint64_t)-1 || qe == (uint64_t)-1) return k; + ///push qs, ts + kv_pushp_cl(k_mer_hit, (*cl), &p); + p->readID = ((uint32_t)(0x7fffffff)); + p->cnt = (uint32_t)-1; p->strand = 0; + p->self_offset = qs; p->offset = ts; + + for (; k < cln && cl->list[k].readID == readID && cl->list[k].self_offset < qe; k++); + + ///push qe, te + kv_pushp_cl(k_mer_hit, (*cl), &p); + p->readID = ((uint32_t)(0x7fffffff)); + p->cnt = (uint32_t)-1; p->strand = 1; + p->self_offset = qe-1; p->offset = te-1; + return k; +} + +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++) { + if((a[k].self_offset >= qs) && (a[k].self_offset < qe) && (!is_alnw(a[k]))) { + a[k].readID = dp; } } - if(qs != (uint64_t)-1) sub_ciagar_gen(z, qs, qe, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); + return k; } -void ul_gap_filling(overlap_region_alloc* ol, kv_ul_ov_t *aln, uint64_t wl, -const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, -double e_rate, int64_t ql, uint64_t rid) +void gen_weight_khits(asg64_v* idx, k_mer_hit *a, int64_t an) { - int64_t i, on = ol->length; - for (i = 0; i < on; i++) { - if(ol->list[i].align_length == (uint32_t)-1) continue; - cigar_gen(&(ol->list[i]), aln->a+ol->list[i].align_length, aln->n-ol->list[i].align_length, i, wl, - uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); + int64_t i, idx_n = idx->n, dp, beg, end, k; + for (i = k = 0, dp = 0, beg = 0, end = -1; i < idx_n; ++i) {///[beg, end) but coordinates in idx is [, ] + ///if idx->a.a[] is qe + if ((idx->a[i]>>32)&1) { + --dp; end = (idx->a[i]>>33)+1; + }else { + //meet a new overlap; the overlaps are pushed by the x_pos_s + ++dp; end = (idx->a[i]>>33); + } + + if(end > beg) k = gen_weight_khits0(beg, end, a, an, k, dp); + beg = end; } } +int64_t gen_cns_chain(overlap_region *z, Candidates_list *cl, asg64_v* iidx, int64_t max_lgap, double sgap_rate) +{ + 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; + qs = qe = ts = te = -1; ci = z->shared_seed; id = cl->list[ci].readID; + for (k = aln_n = 0; k < wn; k++) { + if(z->w_list.a[k].y_end == -1) continue; + if(z->w_list.a[k].extra_end < 0) {///anchor + if(qs == -1) { + qs = z->w_list.a[k].x_start; ts = z->w_list.a[k-1].y_end+1; + } + qe = z->w_list.a[k].x_end+1; te = z->w_list.a[k].y_end+1; + } else { + if(qs != -1) { + fprintf(stderr, "[M::%s::] q::[%ld, %ld)\n", __func__, qs, qe); + ci = push_adp_k_hits(cl, qs, qe, ts, te, id, ci); + aln_n++;//[qs, qe); [ts, te) + } + qs = qe = ts = te = -1; + } + } + if(qs != -1) { + fprintf(stderr, "[M::%s::] q::[%ld, %ld)\n", __func__, qs, qe); + ci = push_adp_k_hits(cl, qs, qe, ts, te, id, ci); + aln_n++;//[qs, qe); [ts, te) + } + + push_adp_k_hits(cl, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, id, ci); + + ka = cl->list+rcn; kn = cl->length-rcn; + if(iidx) gen_weight_khits(iidx, ka, kn); + kn = lchain_dp_trace(ka, kn, max_lgap, sgap_rate, SGAP); cl->length = rcn + kn; + cl->length = rcn; + return kn; +} + + +uint64_t cigar_gen_cns(overlap_region *z, Candidates_list *cl, ul_ov_t *ov, uint64_t on, uint64_t qn, uint64_t wl, +const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate, +int64_t ql, uint64_t rid, asg64_v* iidx, ul_ov_t *des) +{ + if(on <= 0 || iidx->n <= 0) return 0; + uint64_t i; + for (i = 0; i < on; i++) { + assert((i<=0)||(ov[i].qs > ov[i-1].qe)); + sub_ciagar_gen(z, cl, ov[i].qs, ov[i].qe, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); + } + return on; +} ///[ys, ye) -uint64_t get_win_aln(overlap_region *z, uint64_t wid, int64_t *ys, int64_t *ye, int64_t *err) +uint64_t inline get_win_aln(overlap_region *z, uint64_t wid, int64_t *ys, int64_t *ye, int64_t *err) { (*err) = -2; if((wid > 0) && (z->w_list.a[wid].y_end != -1) && (z->w_list.a[wid-1].y_end != -1) && (z->w_list.a[wid].y_end > z->w_list.a[wid-1].y_end)) { @@ -12895,12 +13116,397 @@ char* retrive_str_piece_exz(All_reads *rref, const ul_idx_t *uref, char *buf, in return buf; } +uint64_t k_hits_bcheck(All_reads *rref, const ul_idx_t *uref, overlap_region_alloc* ol, Candidates_list *cl, +uint64_t khit, uint64_t *a, uint64_t a_n, char *str0, char *str1) +{ + if(a_n < 2) return 1; + uint64_t k, e; char *ref0, *ref1; overlap_region *z; + + e = cl->list[(uint32_t)a[0]].offset; + if(e >= khit) e-=khit; + else return 0; + z = &(ol->list[cl->list[(uint32_t)a[0]].readID]); + ref0 = retrive_str_piece_exz(rref, uref, str0, e, khit, z->y_pos_strand, z->y_id); + fprintf(stderr, "\n[M::%s] str0::%.*s\n", __func__, (int32_t)khit, ref0); + + for (k = 1; k < a_n; k++) { + e = cl->list[(uint32_t)a[k]].offset; + if(e >= khit) e-=khit; + else return 0; + z = &(ol->list[cl->list[(uint32_t)a[k]].readID]); + ref1 = retrive_str_piece_exz(rref, uref, str1, e, khit, z->y_pos_strand, z->y_id); + fprintf(stderr, "[M::%s] str1::%.*s\n", __func__, (int32_t)khit, ref1); + if(memcmp(ref0, ref1, khit)) return 0; + } + return 1; +} + +void count_k_hits(All_reads *rref, const ul_idx_t *uref, UC_Read *buf, +overlap_region_alloc* ol, Candidates_list *cl, asg64_v* idx, uint64_t khit, uint64_t basec) +{ + int64_t k, l, on = ol->length, m = 0, i, cn = cl->length; + overlap_region *z; uint64_t t, ff, pid; char *str0, *str1; + for (k = 0, idx->n = 0; k < on; k++) { + z = &(ol->list[k]); i = z->shared_seed; pid = cl->list[i].readID; z->shared_seed = m; + fprintf(stderr, "[M::%s::] utg%.6dl(%c)\n", __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand]); + for (; i < cn && cl->list[i].readID == pid; i++) { + cl->list[m] = cl->list[i]; + cl->list[m].readID = k; + cl->list[m].cnt = 0; + t = cl->list[m].self_offset; t <<= 32; t |= m; + kv_push(uint64_t, (*idx), t); + m++; + } + } + cl->length = cn = m; radix_sort_bc64(idx->a, idx->a+idx->n); + resize_UC_Read(buf, (khit<<1)); str0 = buf->seq; str1 = buf->seq + khit; + + for (k = 1, l = 0; k <= (int64_t)idx->n; k++) { + if(k == cn || (idx->a[l]>>32) != (idx->a[k]>>32)) { + ff = k - l; + if(basec && ff > 1) { + if(!k_hits_bcheck(rref, uref, ol, cl, khit, idx->a+l, k-l, str0, str1)) ff = 1; + } + for (i = l; i < k; i++) { + cl->list[(uint32_t)idx->a[i]].cnt = ff; + fprintf(stderr, "[M::%s::] pos::%u, cnt::%lu\n", __func__, + cl->list[(uint32_t)idx->a[i]].self_offset, ff); + } + + l = k; + } + } + + +} + + +void tuning_ext_offset(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, int64_t ql, int64_t tl, int64_t wl, +int64_t fusion_k_len, int64_t fusion_win_occ, int64_t *ch_s0, int64_t *ch_e0, int64_t *qs0, int64_t *qe0, int64_t *ts0, int64_t *te0, int64_t *mode) +{ + ///[qs, qe) + if((*mode) != 1 && (*mode) != 2) return; + int64_t qs = *qs0, qe = *qe0, ts = *ts0, te = *te0, ch_s = *ch_s0, ch_e = *ch_e0; + int64_t ke = (*qe0)+fusion_k_len, ks=(*qs0)-fusion_k_len, p[3], we, ws, k, wid; + int64_t sl, gq, gt, gl, gg, update, min_g = INT32_MAX, min_id = -1; + adjust_ext_offset(&qs, &qe, &ts, &te, ql, tl, 0, *mode); p[0] = p[1] = p[2] = -1; + if((*mode) == 1) {///forward extension + //qs0 and ts0 are fixed; te0 = -1, qe0 is unreliable + if(qe > (*qe0)) { + ///find k-mer hit + for (ch_e = ch_s; (ch_e < ch_n) && (ch_a[ch_e].self_offset <= ke); ch_e++) { + if(ch_a[ch_e].self_offset<(*qe0)) continue; + if((ch_a[ch_e].offset<=(*ts0))||(ch_a[ch_e].self_offset<=(*qs0))) continue;//not co-linear + if(is_pri_aln(ch_a[ch_e])) { + p[2] = ch_e; break; + } else if(ch_a[ch_e].cnt > 1 && p[1] == -1) { + p[1] = ch_e; + } else if(p[0] == -1) { + p[0] = ch_e; + } + } + if(p[2] != -1) p[0] = p[2]; + else if(p[1] != -1) p[0] = p[1]; + if(p[0] != -1) { + (*mode) = 0;///global + (*qe0) = ch_a[p[0]].self_offset; (*te0) = ch_a[p[0]].offset; (*ch_e0) = p[0]; + return; + } + + ///find aligned window + we = (*qe0); we/=wl; we *= wl; we +=wl; we--; ///next window + for (k = 0; we < qe && k < fusion_win_occ; we+=wl) {//[ws, we]; [qs, qe) + wid = get_win_id_by_e(z, we, wl, NULL); + if(z->w_list.a[wid].y_end == -1) continue;//unmapped + if(z->w_list.a[wid].x_end < (*qs0)) continue; + if(z->w_list.a[wid].y_end < (*ts0)) continue; + sl = z->w_list.a[wid].x_end+1-z->w_list.a[wid].x_start; + if(z->w_list.a[wid].error > (sl/A_L)) continue; + gq = z->w_list.a[wid].x_end-(*qs0); + gt = z->w_list.a[wid].y_end-(*ts0); + gl = MIN(gq, gt); gl/=16; if(gl <= 0) gl = 1; + gg = (gq>=gt)?(gq-gt):(gt-gq); gg /= gl; + update = 0; + if(min_g>gg) { + update = 1; + } else if((min_g==gg)&&(z->w_list.a[min_id].error>z->w_list.a[wid].error)) { + update = 1; + } + if(update) { + min_g = gg; min_id = wid; + } + k++; + } + } + + if(min_id != -1) { + (*qe0) = z->w_list.a[min_id].x_end+1; + (*te0) = z->w_list.a[min_id].y_end+1; + (*mode) = 0;///global + } else { + (*qe0) = qe; (*te0) = te;//extension + } + } else if((*mode) == 2) {///backward extension + //qe0 and te0 are fixed; ts0 = -1, qs0 is unreliable + if(qs < (*qs0)) { + ///find k-mer hit + for (ch_s = ch_e; (ch_s > 0) && (ch_a[ch_s].self_offset >= ks); ch_s--) { + if(ch_a[ch_s].self_offset>(*qs0)) continue; + if((ch_a[ch_s].offset>=(*te0))||(ch_a[ch_s].self_offset>=(*qe0))) continue;//not co-linear + if(is_pri_aln(ch_a[ch_s])) { + p[2] = ch_s; break; + } else if(ch_a[ch_s].cnt > 1 && p[1] == -1) { + p[1] = ch_s; + } else if(p[0] == -1) { + p[0] = ch_s; + } + } + if(p[2] != -1) p[0] = p[2]; + else if(p[1] != -1) p[0] = p[1]; + if(p[0] != -1) { + (*mode) = 0;///global + (*qs0) = ch_a[p[0]].self_offset; (*ts0) = ch_a[p[0]].offset; (*ch_s0) = p[0]; + return; + } + + ///find aligned window + ws = (*qs0)-1; ws/=wl; ws*=wl; + for (k = 0; ws >= qs && k < fusion_win_occ; ws-=wl) {//[ws, we]; [qs, qe) + wid = get_win_id_by_s(z, ws, wl, NULL); + if(z->w_list.a[wid].y_end == -1) continue;//unmapped + if(z->w_list.a[wid].x_start >= (*qe0)) continue; + if(z->w_list.a[wid].y_start >= (*te0)) continue; + sl = z->w_list.a[wid].x_end+1-z->w_list.a[wid].x_start; + if(z->w_list.a[wid].error > (sl/A_L)) continue; + gq = (*qe0) - z->w_list.a[wid].x_start; + gt = (*te0) - z->w_list.a[wid].y_start; + gl = MIN(gq, gt); gl/=16; if(gl <= 0) gl = 1; + gg = (gq>=gt)?(gq-gt):(gt-gq); gg /= gl; + update = 0; + if(min_g>gg) { + update = 1; + } else if((min_g==gg)&&(z->w_list.a[min_id].error>z->w_list.a[wid].error)) { + update = 1; + } + if(update) { + min_g = gg; min_id = wid; + } + k++; + } + } + + if(min_id != -1) { + (*qs0) = z->w_list.a[min_id].x_start; + (*ts0) = z->w_list.a[min_id].y_start; + (*mode) = 0;///global + } else { + (*qs0) = qs; (*ts0) = ts; ///still extension + } + } + + if((*ch_e0) == -1) { + for ((*ch_e0)=(*ch_s0);((*ch_e0)=(*qs0))&&(ch_a[(*ch_e0)].self_offset<(*qe0)); (*ch_e0)++); + } + + if((*ch_s0) == -1) { + for ((*ch_s0)=(*ch_e0);((*ch_s0)>=0)&&(ch_a[(*ch_s0)].self_offset>=(*qs0))&&(ch_a[(*ch_s0)].self_offset<(*qe0)); (*ch_s0)--); + (*ch_s0)++; + } +} + +int64_t cal_estimate_err(overlap_region *z, int64_t wl, int64_t qs, int64_t qe) +{ + int64_t k, ws, we, wid, os, oe, ovlp, tot; + ws = qs/wl; ws *= wl; wid = get_win_id_by_s(z, ws, wl, NULL); + for (k=wid, tot=0; wsw_list.a[k].y_end == -1) continue; + // assert(z->w_list.a[k].x_start==ws); + we = z->w_list.a[k].x_end+1; + os = MAX(qs, ws); oe = MIN(qe, we); + ovlp = ((oe>os)? (oe-os):0); + assert(ovlp); + if(ovlp == (we-ws)) { + tot += z->w_list.a[k].error; + } else { + tot += ((double)z->w_list.a[k].error)*((double)ovlp)/((double)(we-ws)); + } + } + return tot; +} + +void hc_aln_exz_adv(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, const ul_idx_t *uref, +hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, int64_t qs, int64_t qe, int64_t ts, +int64_t te, int64_t mode, int64_t wl, bit_extz_t *exz, int64_t q_tot, double e_rate) +{ + clear_align(*exz); + if(((ts == -1) && (te == -1))) mode = 3;///set to semi-global + int64_t estimate_err = cal_estimate_err(z, wl, qs, qe); + int64_t thre, ql = qe - qs, thre0; + + fprintf(stderr, "[M::%s::ql::%ld] qs::%ld, qe::%ld, ts::%ld, te::%ld, mode::%ld, estimate_err::%ld, e_rate::%f", + __func__, ql, qs, qe, ts, te, mode, estimate_err, e_rate); + if(ql <= MAX_L && (estimate_err*1.2) <= MAX_E) { + thre = scale_ed_thre(estimate_err); if(thre > ql) thre = ql; + if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, q_tot, mode)) { + fprintf(stderr, ", err::%d, thre::%d, scale::%ld(+)\n", exz->err, exz->thre, thre); + return; + } + + thre0 = thre; thre = ql*e_rate; + thre = scale_ed_thre(thre); if(thre > ql) thre = ql; + if(thre > thre0) { + if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, q_tot, mode)) { + fprintf(stderr, ", err::%d, thre::%d, scale::%ld(-)\n", exz->err, exz->thre, thre); + return; + } + } + + thre0 = thre; thre <<= 1; + thre = scale_ed_thre(thre); if(thre > ql) thre = ql; + if(thre > thre0) { + if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, q_tot, mode)) { + fprintf(stderr, ", err::%d, thre::%d, scale::%ld(-)\n", exz->err, exz->thre, thre); + return; + } + } + + thre0 = thre; thre = ql*0.51; + thre = scale_ed_thre(thre); if(thre > ql) thre = ql; + if(thre > thre0) { + if(cal_exz_infi(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, thre, q_tot, mode)) { + fprintf(stderr, ", err::%d, thre::%d, scale::%ld(*)\n", exz->err, exz->thre, thre); + return; + } + } + } + fprintf(stderr, ", err::%d, thre::%d\n", INT32_MAX, exz->thre); +} + + +int64_t sub_base_aln(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, uint64_t pre_e, +uint64_t s, uint64_t e, int64_t wl, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, +char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate, int64_t ql, int64_t tl, +int64_t ch_i, uint64_t rid) +{ + int64_t i = ch_i, l, ibeg, iend, mode, q[2], t[2], ch_s, ch_e; + for (; i >= 0 && ch_a[i].self_offset >= s; i--); + for (((i>=0)?i:0), ibeg=iend=-1 ; i < ch_n && ch_a[i].self_offset < e; i++) { + if((ch_a[i].self_offset >= s) && (ch_a[i].self_offset < e)) { + if(ibeg < 0) ibeg = i; + iend = i+1; + } + } + fprintf(stderr, "\n***[M::%s::rid->%lu] s::%lu, e::%lu, z::[%u, %u)\n", + __func__, rid, s, e, z->x_pos_s, z->x_pos_e+1); + ch_i = i; + if(ibeg > 0 && iend > 0 && iend > ibeg) {///find some anchors[ibeg, iend) + for (i = l = ibeg; i <= iend; i++) { + if(i == iend || is_pri_aln(ch_a[i])) { + q[0] = q[1] = t[0] = t[1] = -1; mode = ch_s = ch_e = -1; + if(l < i && l < iend && is_pri_aln(ch_a[l])) { + q[0] = ch_a[l].self_offset; t[0] = ch_a[l].offset; ch_s = l; + } else { + q[0] = s; + } + + if(i < iend && is_pri_aln(ch_a[i])) { + q[1] = ch_a[i].self_offset; t[1] = ch_a[i].offset; ch_e = i; + } else { + q[1] = e; + } + + if((t[0] != -1) && (t[1] != -1)) { + mode = 0;//global + } else if((t[0] != -1) && (t[1] == -1)) { + mode = 1;///forward extension + } else if((t[0] == -1) && (t[1] != -1)) { + mode = 2;///backward extension + } else { + mode = 3;//semi-global + } + + if(mode == 1 || mode == 2) { + tuning_ext_offset(z, ch_a, ch_n, ql, tl, wl, 256, 4, &ch_s, &ch_e, &q[0], &q[1], &t[0], &t[1], &mode); + } + + hc_aln_exz_adv(z, ch_a+ch_s, ch_e-ch_s, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, ql, e_rate); + + l = i; + } + } + } else {//totoally no anchor; probably semi-global + + } + + return ch_i; +} + +void cigar_gen_by_chain(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, ul_ov_t *ov, int64_t on, uint64_t wl, +const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate, +int64_t ql, uint64_t rid) +{ + if(on <= 0) return; + int64_t i, wn = z->w_list.n, ch_i, tl, id = z->y_id; uint64_t pe = (uint64_t)-1; + if(hpc_g) tl = hpc_len(*hpc_g, id); + else if(uref) tl = uref->ug->u.a[id].len; + else tl = Get_READ_LENGTH((*rref), id); + for (i = 0; i < wn; i++) z->w_list.a[i].clen = 0;///clean cigar + + for (i = ch_i = 0; i < on; i++) { + assert((i<=0)||(ov[i].qs > ov[i-1].qe)); + ch_i = sub_base_aln(z, ch_a, ch_n, pe, ov[i].qs, ov[i].qe, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, tl, ch_i, rid); + pe = ov[i].qe; + } +} + +void ul_gap_filling_adv(overlap_region_alloc* ol, Candidates_list *cl, kv_ul_ov_t *aln, uint64_t wl, +const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, +asg64_v* buf, asg64_v* iidx, double e_rate, int64_t ql, uint64_t rid, int64_t khit, int64_t base_chekc_k_hit, +int64_t max_lgap) +{ + int64_t k, l, ch_n, a_n = aln->n; overlap_region *z; k_mer_hit *ch_a; + count_k_hits(rref, uref, tu, ol, cl, buf, khit, base_chekc_k_hit); + for (k = 1, l = 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); + ch_n = gen_cns_chain(z, cl, iidx, max_lgap, e_rate); + if(ch_n) { + ch_a = cl->list + cl->length; + cigar_gen_by_chain(z, ch_a, ch_n, aln->a+l, k-l, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); + // m = fusion_coordinates(z, ch_a, ch_n, aln->a+l, k-l); + } + // m += cigar_gen_cns(z, cl, aln->a+l, k-l, i, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid, aln->a+m); + l = k; + } + } +} + +inline uint32_t ovlp_win_check(overlap_region *z, uint32_t id0, uint32_t id1, int64_t max_lgap, double small_bw_rate, int64_t min_small_bw) +{ + if(id0 == (uint32_t)-1 || id1 == (uint32_t)-1) return 1; + int64_t qs0, qe0, ts0, te0, qs1, qe1, ts1, te1, err, dd, dm, dq, dr; + if(!get_win_aln(z, id0, &ts0, &te0, &err)) return 1; + qs0 = z->w_list.a[id0].x_start; qe0 = z->w_list.a[id0].x_end+1; + if(!get_win_aln(z, id1, &ts1, &te1, &err)) return 1; + qs1 = z->w_list.a[id1].x_start; qe1 = z->w_list.a[id1].x_end+1; + + if(qs1 < qs0 || qe1 < qe0) return 0; + if(ts1 < ts0 || te1 < te0) return 0; + dq = qe1 - qs0; dr = te1 - ts0; dd = dq>=dr? ((dq)-(dr)): ((dr)-(dq)); + if(ts1 < te0) {//has overlap in y + dm = dq>=dr?dr:dq; + if((dd > (dm*small_bw_rate)) && (dd > min_small_bw)) return 0; + } else { + if(dd > max_lgap) return 0; + } + return 1; +} + uint64_t dp_commen_sketch(kv_ul_ov_t *aln, overlap_region_alloc* ol, uint64_t *id_a, int64_t id_n, -uint64_t *win_a, int64_t win_n, uint64_t *dp, int64_t n_skip, int64_t wl) +uint64_t *win_a, int64_t win_n, uint64_t *dp, int64_t n_skip, int64_t wl, int64_t cov, int64_t max_lgap, double sgap_rate, int64_t sgap) { if(!win_n) return 0; - int64_t k, ws, ws0, i, m, wid, wid0, s, e, err, s0, e0, err0, sc, max_sc, p, long_sc, long_idx; - overlap_region *z; s = e = err = s0 = e0 = err0 = -1; + int64_t k, ws, ws0, i, m, wid, wid0, sc, max_sc, p, long_sc, long_idx; overlap_region *z; if(n_skip < win_n) { long_sc = long_idx = -1; for (k = 0; k < n_skip; k++) { @@ -12918,10 +13524,7 @@ uint64_t *win_a, int64_t win_n, uint64_t *dp, int64_t n_skip, int64_t wl) z = &(ol->list[aln->a[(uint32_t)id_a[i]].qn]); wid = get_win_id_by_s(z, ws, wl, NULL); wid0 = get_win_id_by_s(z, ws0, wl, NULL); - - get_win_aln(z, wid, &s, &e, &err); - get_win_aln(z, wid0, &s0, &e0, &err0); - if(e0 > s) break; + if(!ovlp_win_check(z, wid0, wid, max_lgap, sgap_rate, sgap)) break; } if((i >= id_n) && ((sc + ((uint32_t)dp[m])) > max_sc)) { max_sc = (sc + ((uint32_t)dp[m])); p = m; @@ -12946,21 +13549,23 @@ uint64_t *win_a, int64_t win_n, uint64_t *dp, int64_t n_skip, int64_t wl) win_n = m; } - for (k = 0; k < win_n; k++) { + for (k = 0, cov = -cov; k < win_n; k++) { ws = win_a[k]; for (i = 0; i < id_n; i++) { z = &(ol->list[aln->a[(uint32_t)id_a[i]].qn]); wid = get_win_id_by_s(z, ws, wl, NULL); - z->w_list.a[wid].extra_end = -1; - get_win_aln(z, wid, &s, &e, &err); - if(z->align_length < (uint32_t)e) z->align_length = e;///for conliner + if(cov > INT16_MIN) z->w_list.a[wid].extra_end = cov; + else z->w_list.a[wid].extra_end = INT16_MIN; + z->align_length = wid;///for conliner } } return win_n; } + uint64_t gen_commen_sketch(All_reads *rref, const ul_idx_t *uref, overlap_region_alloc* ol, uint64_t *id_a, uint64_t id_n, uint64_t s, uint64_t e, uint64_t ql, uint64_t wl, -uint64_t *buf, uint64_t dp, char *str0, char *str1, kv_ul_ov_t *aln, asg64_v *trace)///[s, e) +uint64_t *buf, uint64_t dp, char *str0, char *str1, kv_ul_ov_t *aln, asg64_v *trace, +int64_t max_lgap, double sgap_rate, int64_t sgap)///[s, e) { if(!id_n) return id_n; uint64_t i, m, k, rm_n = 0, buf_n = 0, qs, qe, wid, co, occ; char *qstring, *tstring; @@ -12972,7 +13577,8 @@ uint64_t *buf, uint64_t dp, char *str0, char *str1, kv_ul_ov_t *aln, asg64_v *tr //idx_a[] is sorted by aln[].qs for (k = 0; k < id_n; k++) { if(aln->a[id_a[k]].qs<=qs && aln->a[id_a[k]].qe>=qe) { - buf[buf_n++] = id_a[k]; + buf[buf_n] = ol->list[aln->a[id_a[k]].qn].align_length; buf[buf_n] <<= 32; buf[buf_n] |= id_a[k]; + buf_n++; } if(aln->a[id_a[k]].qe < e) rm_n++; } @@ -12987,13 +13593,13 @@ uint64_t *buf, uint64_t dp, char *str0, char *str1, kv_ul_ov_t *aln, asg64_v *tr z = &(ol->list[aln->a[(uint32_t)buf[0]].qn]); wid = get_win_id_by_s(z, ws, wl, NULL); if(!get_win_aln(z, wid, &(r_y[0]), &(r_y[1]), &r_err)) continue; - if(r_y[0] < z->align_length) continue;///not co-linear + if(!ovlp_win_check(z, z->align_length, wid, max_lgap, sgap_rate, sgap)) continue;///not co-linear qstring = tstring = NULL; for (i = 1; i < buf_n; i++) { z = &(ol->list[aln->a[(uint32_t)buf[i]].qn]); wid = get_win_id_by_s(z, ws, wl, NULL); if(!get_win_aln(z, wid, &(p_y[0]), &(p_y[1]), &p_err)) break; - if(p_y[0] < z->align_length) continue;///not co-linear + if(!ovlp_win_check(z, z->align_length, wid, max_lgap, sgap_rate, sgap)) continue;///not co-linear if(((r_y[1]-r_y[0]) != (p_y[1]-p_y[0])) || (r_err != p_err)) break; if(r_err == 0) continue; @@ -13002,7 +13608,7 @@ uint64_t *buf, uint64_t dp, char *str0, char *str1, kv_ul_ov_t *aln, asg64_v *tr ol->list[aln->a[(uint32_t)buf[0]].qn].y_pos_strand, ol->list[aln->a[(uint32_t)buf[0]].qn].y_id); } tstring = retrive_str_piece_exz(rref, uref, str1, p_y[0], p_y[1]-p_y[0], z->y_pos_strand, z->y_id); - if(memcmp(str0, str1, (p_y[1]-p_y[0]))) { + if(memcmp(qstring, tstring, (p_y[1]-p_y[0]))) { // fprintf(stderr, "[M::%s::] qs::%ld, qe::%lu, ts::%ld, te::%lu, err::%ld, rts::%ld, rte::%lu, err::%ld\n", // __func__, ws, we, p_y[0], p_y[1], r_err, r_y[0], r_y[1], p_err); // fprintf(stderr, "str0::%.*s\n", (int32_t)(r_y[1]-r_y[0]), str0); @@ -13011,23 +13617,21 @@ uint64_t *buf, uint64_t dp, char *str0, char *str1, kv_ul_ov_t *aln, asg64_v *tr } } if(i < buf_n) continue; - for (i = 0, p_y[0] = p_y[1] = -1; i < buf_n; i++) { - z = &(ol->list[aln->a[(uint32_t)buf[i]].qn]); - wid = get_win_id_by_s(z, ws, wl, NULL); - - if(co) { - get_win_aln(z, wid, &(p_y[0]), &(p_y[1]), &p_err); - if((buf[i]>>32) > (uint64_t)p_y[0]) co = 0; - m = p_y[1]; m <<= 32; m |= (uint32_t)buf[i]; buf[i] = m; + if(co) { + for (i = 0; i < buf_n; i++) { + z = &(ol->list[aln->a[(uint32_t)buf[i]].qn]); + wid = get_win_id_by_s(z, ws, wl, NULL); + co = ovlp_win_check(z, buf[i]>>32, wid, max_lgap, sgap_rate, sgap); + m = wid; m <<= 32; m |= (uint32_t)buf[i]; buf[i] = m; } + if(co) occ++; } - if(co) occ++; kv_push(uint64_t, *trace, ws); } if(trace->n) { if(!co) kv_resize(uint64_t, *trace, trace->n<<1); - trace->n = dp_commen_sketch(aln, ol, buf, buf_n, trace->a, trace->n, trace->a + trace->n, occ, wl); + trace->n = dp_commen_sketch(aln, ol, buf, buf_n, trace->a, trace->n, trace->a + trace->n, occ, wl, dp, max_lgap, sgap_rate, sgap); } } if(rm_n) { @@ -13041,8 +13645,10 @@ uint64_t *buf, uint64_t dp, char *str0, char *str1, kv_ul_ov_t *aln, asg64_v *tr } void update_sketch_trace(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t *uopt, -All_reads *rref, UC_Read* tu, asg64_v* idx, asg64_v *b0, asg64_v *b1, int64_t ql, int64_t wl, kv_ul_ov_t *aln, uint64_t rid) +All_reads *rref, UC_Read* tu, asg64_v* idx, asg64_v *b0, asg64_v *b1, int64_t ql, int64_t wl, +kv_ul_ov_t *aln, uint64_t rid, int64_t max_lgap, double sgap_rate) { + idx->n = 0; if(!aln->n) return; uint64_t i, k, own, srt_n; int64_t dp, old_dp, beg, end; overlap_region *z; for (i = 0; i < ol->length; i++) { @@ -13056,14 +13662,14 @@ All_reads *rref, UC_Read* tu, asg64_v* idx, asg64_v *b0, asg64_v *b1, int64_t ql kv_resize(uint64_t, *idx, (aln->n<<1)); kv_resize(uint64_t, *b0, aln->n); for (i = srt_n = 0; i < aln->n; i++) { // if(i == 0 || aln->a[i].qn != aln->a[i-1].qn) ol->list[aln->a[i].qn].align_length = i; - ol->list[aln->a[i].qn].align_length = 0;///for co-linear + ol->list[aln->a[i].qn].align_length = (uint32_t)-1;///for co-linear idx->a[srt_n] = aln->a[i].qs<<1; idx->a[srt_n] <<= 32; idx->a[srt_n] += i; srt_n++; idx->a[srt_n] = ((aln->a[i].qe-1)<<1)+1; idx->a[srt_n] <<= 32; idx->a[srt_n] += i; srt_n++; aln->a[i].el = 1; } radix_sort_bc64(idx->a, idx->a+srt_n); idx->n = srt_n; resize_UC_Read(tu, (wl<<1)); - for (i = 0, dp = 0, beg = 0, end = -1; i < srt_n; ++i) {///[beg, end] + for (i = 0, dp = 0, beg = 0, end = -1; i < srt_n; ++i) {///[beg, end), but the idx saves [qs, qe] old_dp = dp; ///if idx->a.a[] is qe if ((idx->a[i]>>32)&1) { @@ -13076,7 +13682,8 @@ All_reads *rref, UC_Read* tu, asg64_v* idx, asg64_v *b0, asg64_v *b1, int64_t ql if((end > beg) && (end - beg > wl) && (old_dp >= 2) ) { idx->n = srt_n + - gen_commen_sketch(rref, uref, ol, idx->a+srt_n, idx->n-srt_n, beg, end, ql, wl, b0->a, old_dp, tu->seq, tu->seq+wl, aln, b1); + gen_commen_sketch(rref, uref, ol, idx->a+srt_n, idx->n-srt_n, beg, end, ql, wl, b0->a, + old_dp, tu->seq, tu->seq+wl, aln, b1, max_lgap, sgap_rate, SGAP); } beg = end; } @@ -13084,7 +13691,7 @@ All_reads *rref, UC_Read* tu, asg64_v* idx, asg64_v *b0, asg64_v *b1, int64_t ql for (i = srt_n = 0; i < aln->n; i++) { if(i == 0 || aln->a[i].qn != aln->a[i-1].qn) ol->list[aln->a[i].qn].align_length = i; } - + idx->n = srt_n; return; } @@ -13092,10 +13699,10 @@ All_reads *rref, UC_Read* tu, asg64_v* idx, asg64_v *b0, asg64_v *b1, int64_t ql void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, char *qstr, uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, bit_extz_t *exz, haplotype_evdience_alloc* hap, kvec_t_u64_warp* v_idx, - double e_rate, int64_t wl, kv_ul_ov_t *aln, int64_t sid, void *km) + double e_rate, int64_t wl, kv_ul_ov_t *aln, int64_t sid, uint64_t khit, void *km) { uint64_t i, bs, k, ovl/**, on**/; Window_Pool w; double err; - /**int64_t sc;**/ overlap_region t; overlap_region *z; + /**int64_t sc;**/ overlap_region t; overlap_region *z; asg64_v iidx, buf; ol->mapped_overlaps_length = 0; if(ol->length <= 0) return; @@ -13108,7 +13715,8 @@ void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur if(!aln) { resize_UC_Read(qu, ql); qu->length = ql; memcpy(qu->seq, qstr, ql); for (i = 0; i < ol->length; i++) { - z = &(ol->list[i]); ovl = z->x_pos_e+1-z->x_pos_s; + z = &(ol->list[i]); ovl = z->x_pos_e+1-z->x_pos_s; + z->shared_seed = z->non_homopolymer_errors;///for index if(!align_ul_ed_post_extz(z, uref, NULL, qu->seq, tu->seq, exz, err, w.window_length, -1, km)) { continue; } @@ -13143,7 +13751,10 @@ void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur } else { // fprintf(stderr, "-[M::%s] on::%lu\n", __func__, ol->length); if(ol->length <= 1) return; - ul_gap_filling(ol, aln, wl, uref, NULL, NULL, qu->seq, tu, exz, err, ql, sid); + ///coordinates for all intervals with cov > 1 + copy_asg_arr(iidx, hap->snp_srt); copy_asg_arr(buf, v_idx->a); + ul_gap_filling_adv(ol, cl, aln, wl, uref, NULL, NULL, qu->seq, tu, exz, &buf, &iidx, err, ql, sid, khit, 1, MAX_LGAP(ql)); + copy_asg_arr(hap->snp_srt, iidx); copy_asg_arr(v_idx->a, buf); // for (i = 0; (i < ol->length) && (ol->list[i].is_match == 1); i++); on = i; // if(on <= 1) return; // kv_resize(uint64_t, v_idx->a, (on<<1)); v_idx->a.n = on; diff --git a/Correct.h b/Correct.h index c8cef1e..da8c7f9 100644 --- a/Correct.h +++ b/Correct.h @@ -19,6 +19,8 @@ #define DELETION 3 #define ERROR_RATE 1.25 #define UL_TOPN 50 +#define SGAP 16 +#define MAX_LGAP(ql) ((((ql)*0.33)>2000)?((ql)*0.33):2000) #define WINDOW_MAX_SIZE (WINDOW + (int)(1.0 / HA_MIN_OV_DIFF) + 3) // TODO: why 1/max_ov_diff? @@ -1135,7 +1137,7 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, char *qstr, uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, bit_extz_t *exz, haplotype_evdience_alloc* hap, kvec_t_u64_warp* v_idx, - double e_rate, int64_t wl, kv_ul_ov_t *aln, int64_t sid, void *km); + double e_rate, int64_t wl, kv_ul_ov_t *aln, int64_t sid, uint64_t hpc_k, void *km); void ul_lalign_old_ed(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, char *qstr, uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, @@ -1325,7 +1327,8 @@ void inline resize_UC_Read(UC_Read *z, int64_t s) } void update_sketch_trace(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t *uopt, -All_reads *rref, UC_Read* tu, asg64_v* idx, asg64_v *b0, asg64_v *b1, int64_t ql, int64_t wl, kv_ul_ov_t *aln, uint64_t rid); +All_reads *rref, UC_Read* tu, asg64_v* idx, asg64_v *b0, asg64_v *b1, int64_t ql, int64_t wl, +kv_ul_ov_t *aln, uint64_t rid, int64_t max_lgap, double sgap_rate); #define copy_asg_arr(des, src) ((des).a = (src).a, (des).n = (src).n, (des).m = (src).m) diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 0da0d87..2b6c1e2 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -1622,4 +1622,82 @@ uint64_t lchain_dp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, ov res->overlapLen = get_chainLen(res->x_pos_s, res->x_pos_e, xl, res->y_pos_s, res->y_pos_e, yl); for (i = 0; i < cL; i++) des[i] = a[t[cL-i-1]]; return cL; +} + + +inline int64_t hit_long_gap(k_mer_hit *a, k_mer_hit *b, int64_t max_lgap, double small_bw_rate, int64_t min_small_bw) +{ + int64_t dq, dr, dd, dm; + dq = b->self_offset-a->self_offset; + dr = b->offset-a->offset; + dd = dq>=dr? ((dq)-(dr)): ((dr)-(dq)); + if(max_lgap>=0) { + if(dd <= max_lgap) return 1; + return 0; + } else { + dm = dq>=dr?dr:dq; + if((dd > (dm*small_bw_rate)) && (dd > min_small_bw)) return 0; + return 1; + } +} + +int64_t filter_bad_seed_dp(k_mer_hit *sk, k_mer_hit *ek, k_mer_hit* a, int64_t a_n, int64_t max_lgap, double small_bw_rate, int64_t min_small_bw) +{ + int64_t k = 0; k_mer_hit *z; double bw_r; int64_t bw, mmgap, occ = 0; + bw_r = small_bw_rate; bw = min_small_bw; mmgap = max_lgap; + if(sk) { + for (k = 0; k < a_n; k++) { + z = &(a[k]); + if(z->cnt < z->readID) continue; + mmgap = max_lgap; + if(z->cnt <= 1) mmgap = -1; + if((sk && (!hit_long_gap(sk, z, mmgap, bw_r, bw))) || + (ek && (!hit_long_gap(z, ek, mmgap, bw_r, bw)))) { + z->offset = z->self_offset = (uint32_t)-1; occ++; + } else { + break; + } + } + } + + if(ek && k < a_n) { + for (k = a_n-1; k >= 0; k--) { + z = &(a[k]); + if(z->cnt < z->readID) continue; + mmgap = max_lgap; + if(z->cnt <= 1) mmgap = -1; + if((sk && (!hit_long_gap(sk, z, mmgap, bw_r, bw))) || + (ek && (!hit_long_gap(z, ek, mmgap, bw_r, bw)))) { + z->offset = z->self_offset = (uint32_t)-1; occ++; + } else { + break; + } + } + } + + return occ; +} + +uint64_t lchain_dp_trace(k_mer_hit* a, int64_t a_n, int64_t max_lgap, double sgap_rate, int64_t sgap) +{ + if(a_n <= 0) return 0; + int64_t i, st, occ = 0; + + for (i = 1, st = 0; i <= a_n; ++i) { + if((i == a_n) || (is_alnw(a[i]))) {///[st, i) + if(i > st) { + occ += filter_bad_seed_dp((st>0)?&(a[st-1]):NULL, (i= (a).readID)) + +uint64_t lchain_dp_trace(k_mer_hit* a, int64_t a_n, int64_t max_lgap, double sgap_rate, int64_t sgap); + #endif diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index 7938777..d286ef2 100644 --- a/Levenshtein_distance.h +++ b/Levenshtein_distance.h @@ -543,7 +543,7 @@ inline uint32_t pop_trace(asg16_v *res, uint32_t i, uint16_t *c, uint32_t *len) ///511 -> 32 64-bits #define MAX_E 1023 -#define MAX_L 8192 +#define MAX_L 5120 typedef uint64_t w_sig; #define bitw (6) diff --git a/anchor.cpp b/anchor.cpp index 02c90d1..502a1f3 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -891,6 +891,32 @@ void inline reverse_k_mer_hit(k_mer_hit *a, uint64_t a_n, uint64_t xl, uint64_t } } + +void inline reset_k_mer_hit(k_mer_hit *a, uint64_t a_n, uint64_t xl, uint64_t yl, uint64_t rev, uint64_t *nid) +{ + uint64_t z, han = a_n>>1; k_mer_hit *ai, *aj, ka; + if(rev) { + for (z = 0; z < han; z++) { + ai = &(a[z]); aj = &(a[a_n-z-1]); + ka = (*ai); (*ai) = (*aj); (*aj) = ka; + + ai->self_offset = xl-ai->self_offset-1; + ai->offset = yl-ai->offset-1; + + aj->self_offset = xl-aj->self_offset-1; + aj->offset = yl-aj->offset-1; + if(nid) ai->readID = aj->readID = (*nid); + } + if(a_n&1) { + a[z].self_offset = xl-a[z].self_offset-1; + a[z].offset = yl-a[z].offset-1; + if(nid) a[z].readID = (*nid); + } + } else if(nid) { + for (z = 0; z < a_n; z++) a[z].readID = (*nid); + } +} + void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uint64_t rl, All_reads* rdb, const ul_idx_t *udb, uint32_t apend_be, overlap_region* tf, uint64_t max_n_chain, int64_t max_skip, int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, int64_t quick_check, uint32_t gen_off) @@ -913,10 +939,11 @@ void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uin rl, rdb?Get_READ_LENGTH((*rdb), (*tf).y_id):udb->ug->u.a[(*tf).y_id].len, quick_check); // assert(sm > 0); if(ovlp_chain_gen(ol, tf, rl, rdb?Get_READ_LENGTH((*rdb), (*tf).y_id):udb->ug->u.a[(*tf).y_id].len, apend_be, cl->list+m, sm)) { - r = &(ol->list[ol->length-1]); r->non_homopolymer_errors = m; - if(r->y_pos_strand) { - reverse_k_mer_hit(cl->list+m, sm, rl, rdb?Get_READ_LENGTH((*rdb), r->y_id):udb->ug->u.a[r->y_id].len); - } + r = &(ol->list[ol->length-1]); r->non_homopolymer_errors = m; + // if(r->y_pos_strand) { + // reverse_k_mer_hit(cl->list+m, sm, rl, rdb?Get_READ_LENGTH((*rdb), r->y_id):udb->ug->u.a[r->y_id].len); + // } + reset_k_mer_hit(cl->list+m, sm, rl, rdb?Get_READ_LENGTH((*rdb), r->y_id):udb->ug->u.a[r->y_id].len, r->y_pos_strand, &(ol->length)); if(gen_off) gen_fake_cigar(&(r->f_cigar), r, apend_be, cl->list+m, sm); m += sm; } diff --git a/inter.cpp b/inter.cpp index 6a3172c..1830dc6 100644 --- a/inter.cpp +++ b/inter.cpp @@ -6050,7 +6050,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // if(s->id+i!=41927 && s->id+i!=47072 && s->id+i!=67641 && s->id+i!=90305 && s->id+i!=698342 && s->id+i!=329421) { // return; // } - // if((s->id+i!=230) /**&& (s->id+i!=44) && (s->id+i!=948)**/) return; + if((s->id+i!=1) /**&& (s->id+i!=44) && (s->id+i!=948)**/) 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); @@ -6073,7 +6073,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // memset(&b->self_read, 0, sizeof(b->self_read)); ul_lalign(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, - &b->correct, &b->exz, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, NULL, s->id+i, NULL); + &b->correct, &b->exz, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, NULL, s->id+i, s->opt->k, NULL); // ul_lalign_old_ed(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, // &b->correct, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, 1, NULL); ton = b->olist.length;//all alignments pass similary check @@ -6086,13 +6086,13 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call filter_topN(&b->olist, &(bl->lo), s->len[i], winLen, UL_TOPN, bl); // update_shared_intervals(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b->r_buf, &(s->sps[tid]), s->len[i], winLen, &(bl->lo), s->id+i); - copy_asg_arr(b0, b->r_buf.a); copy_asg_arr(b1, s->sps[tid]); copy_asg_arr(b2, b->hap.snp_srt); + copy_asg_arr(b0, b->hap.snp_srt); copy_asg_arr(b1, s->sps[tid]); copy_asg_arr(b2, b->r_buf.a); // update_shared_intervals(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b0, &b1, &b2, s->len[i], winLen, &(bl->lo), s->id+i); - update_sketch_trace(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b0, &b1, &b2, s->len[i], winLen, &(bl->lo), s->id+i); - copy_asg_arr(b->r_buf.a, b0); copy_asg_arr(s->sps[tid], b1); copy_asg_arr(b->hap.snp_srt, b2); + update_sketch_trace(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b0, &b1, &b2, s->len[i], winLen, &(bl->lo), s->id+i, MAX_LGAP(s->len[i]), s->opt->diff_ec_ul); + copy_asg_arr(b->hap.snp_srt, b0); copy_asg_arr(s->sps[tid], b1); copy_asg_arr(b->r_buf.a, b2); ul_lalign(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, - &b->correct, &b->exz, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, &(bl->lo), s->id+i, NULL); + &b->correct, &b->exz, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, &(bl->lo), s->id+i, s->opt->k, NULL); // ul_lalign_old_ed(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, // &b->correct, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, 0, NULL); }