diff --git a/Assembly.cpp b/Assembly.cpp index 42bf649..36d5ded 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -468,6 +468,7 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov, int is_ug) kv_init(b->k_flag.a); kv_init(b->sp); init_bit_extz_t(&(b->exz), 31); + init_bit_extz_t(&(b->exz1), 31); if(!is_ug) b->ab = ha_abuf_init(); else b->abl = ha_abufl_init(); if (!b->is_final) { @@ -497,6 +498,7 @@ void ha_ovec_destroy(ha_ovec_buf_t *b) kv_destroy(b->k_flag.a); kv_destroy(b->sp); destroy_bit_extz_t(&(b->exz)); + destroy_bit_extz_t(&(b->exz1)); if (!b->is_final) { destory_Cigar_record(&b->cigar1); destory_Graph(&b->POA_Graph); diff --git a/Assembly.h b/Assembly.h index 68ab3ff..7bb38c5 100644 --- a/Assembly.h +++ b/Assembly.h @@ -40,7 +40,7 @@ typedef struct { overlap_region tmp_region; ma_utg_v *ua; st_mt_t sp; - bit_extz_t exz; + bit_extz_t exz, exz1; } ha_ovec_buf_t; int ha_assemble(void); diff --git a/Correct.cpp b/Correct.cpp index 70f0053..85889b3 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -32,6 +32,8 @@ KRADIX_SORT_INIT(overlap_region_dp_srt, overlap_region, overlap_region_dp_key, m #define window_list_xs_key(x) ((x).x_start) KRADIX_SORT_INIT(window_list_xs_srt, window_list, window_list_xs_key, member_size(window_list, x_start)) +#define uov_qs_key(p) ((p).qs) +KRADIX_SORT_INIT(uov_srt_qs, ul_ov_t, uov_qs_key, member_size(ul_ov_t, qs)) int ha_ov_type(const overlap_region *r, uint32_t len); @@ -12600,11 +12602,11 @@ void ul_lalign_old_ed(overlap_region_alloc* ol, Candidates_list *cl, const ul_id } } -inline uint64_t scale_ed_thre(uint32_t err) +inline uint64_t scale_ed_thre(uint32_t err, uint32_t max_err) { uint64_t bd = (err<<1)+1, w; w = (bd>>bitw); w <<= bitw; if(w < bd) w += bitwbit; - err = (w-1)>>1; if(err > MAX_E) err = MAX_E; + err = (w-1)>>1; if(err > max_err) err = max_err; return err; } @@ -12807,7 +12809,7 @@ int64_t cal_exz_infi(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_ 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 mode, int64_t wl, bit_extz_t *exz, int64_t q_tot, double e_rate, int64_t maxl, int64_t maxe) { clear_align(*exz); int64_t thre, ql = qe - qs, thre0, t_tot, id = z->y_id; @@ -12828,15 +12830,15 @@ int64_t mode, int64_t wl, bit_extz_t *exz, int64_t q_tot, double e_rate) // 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(ql <= maxl && (estimate_err*1.2) <= maxe) { + thre = scale_ed_thre(estimate_err, maxe); 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; + thre = scale_ed_thre(thre, maxe); 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); @@ -12845,7 +12847,7 @@ int64_t mode, int64_t wl, bit_extz_t *exz, int64_t q_tot, double e_rate) } thre0 = thre; thre <<= 1; - thre = scale_ed_thre(thre); if(thre > ql) thre = ql; + thre = scale_ed_thre(thre, maxe); 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); @@ -12854,7 +12856,7 @@ int64_t mode, int64_t wl, bit_extz_t *exz, int64_t q_tot, double e_rate) } thre0 = thre; thre = ql*0.51; - thre = scale_ed_thre(thre); if(thre > ql) thre = ql; + thre = scale_ed_thre(thre, maxe); 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); @@ -12946,7 +12948,7 @@ int64_t ql, uint64_t rid) // __func__, q[1]-q[0], q[0], q[1], t[0], t[1], mode, tot_e); // } 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); + 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, MAX_CNS_L, MAX_CNS_E); } l = k; @@ -12987,26 +12989,32 @@ int64_t ql, uint64_t rid, ul_ov_t *des) -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 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 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; + int64_t k = ci; k_mer_hit *p = NULL, *m = NULL; + if(cl->length > cln) { + m = &(cl->list[cl->length-1]); assert(is_alnw(*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; - kv_pushp_cl(k_mer_hit, (*cl), &p); *p = cl->list[k]; p->readID = ((uint32_t)(0x7fffffff)); + kv_pushp_cl(k_mer_hit, (*cl), &p); *p = cl->list[k]; p->readID = p->cnt; + // fprintf(stderr, "[M::%s::+]\n", __func__); } - if(qs == (uint64_t)-1 || qe == (uint64_t)-1) return k; + // if(!(k >= cln || cl->list[k].self_offset >= qe)) { + // fprintf(stderr, "[M::%s::] k::%ld, self_offset::%u, offset::%u\n", + // __func__, k, cl->list[k].self_offset, cl->list[k].offset); + // } + assert(k >= cln || cl->list[k].readID != readID || cl->list[k].self_offset >= qe); ///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 @@ -13030,53 +13038,62 @@ int64_t gen_weight_khits0(uint32_t qs, uint32_t qe, k_mer_hit *a, int64_t an, in void gen_weight_khits(asg64_v* idx, k_mer_hit *a, int64_t an) { - 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 [, ] + int64_t i, idx_n = idx->n, dp, old_dp, beg, end, k; + // fprintf(stderr, "[M::%s::] idx->n::%ld\n", __func__, (int64_t)idx->n); + for (i = k = 0, dp = old_dp = 0, beg = 0, end = -1; i < idx_n; ++i) {///[beg, end) but coordinates in idx is [, ] ///if idx->a.a[] is qe + old_dp = dp; 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); + // fprintf(stderr, "[M::%s::] beg::%ld, end::%ld, old_dp::%ld\n", __func__, beg, end, old_dp); + if(end > beg) k = gen_weight_khits0(beg, end, a, an, k, old_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 gen_cns_chain(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; 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(z->w_list.a[k].extra_end < 0 && z->w_list.a[k].y_end != -1) {///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); + // 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); 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); + // 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); 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); + // 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); 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(iidx) gen_weight_khits(iidx, ka, kn); - kn = lchain_dp_trace(ka, kn, max_lgap, sgap_rate, SGAP); cl->length = rcn + kn; + if(need_filter_khit) { + kn = lchain_dp_trace(ka, kn, max_lgap, sgap_rate, SGAP); cl->length = rcn + kn; + } cl->length = rcn; return kn; } @@ -13117,7 +13134,7 @@ char* retrive_str_piece_exz(All_reads *rref, const ul_idx_t *uref, char *buf, in } 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) +uint64_t khit, uint64_t *a, uint64_t a_n, char* qstr, char *str0, char *str1) { if(a_n < 2) return 1; uint64_t k, e; char *ref0, *ref1; overlap_region *z; @@ -13127,7 +13144,10 @@ uint64_t khit, uint64_t *a, uint64_t a_n, char *str0, char *str1) 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); + // fprintf(stderr, "\n[M::%s::] qstr::%.*s\n", __func__, + // (int32_t)khit, qstr+cl->list[(uint32_t)a[0]].self_offset-khit); + // fprintf(stderr, "[M::%s::qoff->%u::toff->%u::%c] tstr0::%.*s\n", __func__, + // cl->list[(uint32_t)a[0]].self_offset, cl->list[(uint32_t)a[0]].offset, "+-"[z->y_pos_strand], (int32_t)khit, ref0); for (k = 1; k < a_n; k++) { e = cl->list[(uint32_t)a[k]].offset; @@ -13135,42 +13155,224 @@ uint64_t khit, uint64_t *a, uint64_t a_n, char *str0, char *str1) 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); + // fprintf(stderr, "[M::%s::qoff->%u::toff->%u::%c] tstr1::%.*s\n", __func__, + // cl->list[(uint32_t)a[k]].self_offset, cl->list[(uint32_t)a[k]].offset, "+-"[z->y_pos_strand], (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) +inline int64_t khit_long_gap(k_mer_hit *a, k_mer_hit *b, double small_bw_rate, int64_t min_small_bw) { - 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++) { + 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)); + + dm = dq>=dr?dr:dq; + if((dd > (dm*small_bw_rate)) && (dd > min_small_bw)) return 0; + return 1; +} + +int64_t filter_bad_khits(k_mer_hit *sk, k_mer_hit *ek, k_mer_hit* a, int64_t a_n, double small_bw_rate, int64_t min_small_bw) +{ + int64_t k = 0; k_mer_hit *z; double bw_r; int64_t bw, occ = 0; + bw_r = small_bw_rate; bw = min_small_bw; + if(sk) { + for (k = 0; k < a_n; k++) { + z = &(a[k]); + if((sk && (!khit_long_gap(sk, z, bw_r, bw))) || + (ek && (!khit_long_gap(z, ek, 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; + if((sk && (!khit_long_gap(sk, z, bw_r, bw))) || + (ek && (!khit_long_gap(z, ek, bw_r, bw)))) { + z->offset = z->self_offset = (uint32_t)-1; occ++; + } else { + break; + } + } + } + + if(occ) { + for (k = occ = 0; k < a_n; k++) { + if(a[k].offset == (uint32_t)-1) continue; + a[occ++] = a[k]; + } + a_n = occ; + } + return a_n; +} + +int64_t filter_adp_k_hits(k_mer_hit *ka, int64_t kn, Chain_Data *dp, uint64_t qs, uint64_t qe, uint64_t ts, uint64_t te, uint64_t readID, int64_t ci, int64_t *cmi, k_mer_hit *m, +double small_bw_rate, int64_t min_small_bw) +{ + int64_t k = ci, cmi0 = (*cmi); k_mer_hit *p = NULL, n; + for (; k < kn && ka[k].readID == readID && ka[k].self_offset < qs; k++) { + p = &(ka[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; + ka[(*cmi)++] = *p; + } + n.self_offset = qs; n.offset = ts; p = &n; + if(qs == (uint64_t)-1 || qe == (uint64_t)-1) p = NULL; + // fprintf(stderr, "[M::%s::] cmi0::%ld, cmi::%ld\n", __func__, cmi0, (*cmi)); back = (*cmi); + if((*cmi) > cmi0) { + (*cmi) = cmi0 + filter_bad_khits(m, p, ka+cmi0, (*cmi)-cmi0, small_bw_rate, min_small_bw); + if((*cmi) > cmi0) { + (*cmi) = cmi0 + lchain_refine(ka+cmi0, (*cmi)-cmi0, ka+cmi0, dp, 50, 5000, 512, 16); + } + } + // if(back != (*cmi)) { + // fprintf(stderr, "sbsbsbsb[M::%s::] cmi0::%ld, cmi::%ld\n", __func__, cmi0, (*cmi)); + // } + + for (; k < kn && ka[k].readID == readID && ka[k].self_offset < qe; k++); + return k; +} + + +void refine_khits(overlap_region *z, Candidates_list *cl, Chain_Data *dp, double sgap_rate) +{ + int64_t k, wn = z->w_list.n, aln_n, qs, qe, ts, te, ci, cmi, id, mm_gap = 64; + k_mer_hit p; p.self_offset = (uint32_t)-1; p.offset = (uint32_t)-1; + if(wn <= 0) return; + qs = qe = ts = te = -1; ci = cmi = z->shared_seed; id = cl->list[ci].readID; + for (k = aln_n = 0; k < wn; k++) { + if(z->w_list.a[k].extra_end < 0 && z->w_list.a[k].y_end != -1) {///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, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld\n", + // __func__, qs, qe, ts, te, ci); + ci = filter_adp_k_hits(cl->list, cl->length, dp, qs, qe, ts, te, id, ci, &cmi, + aln_n?&p:NULL, sgap_rate, mm_gap); + p.self_offset = qe - 1; p.offset = te -1; + 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\n", + // __func__, qs, qe, ts, te, ci); + ci = filter_adp_k_hits(cl->list, cl->length, dp, qs, qe, ts, te, id, ci, &cmi, + aln_n?&p:NULL, sgap_rate, mm_gap); + p.self_offset = qe - 1; p.offset = te -1; + aln_n++;//[qs, qe); [ts, te) + } + // fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld\n", + // __func__, qs, qe, ts, te, ci); + ci = filter_adp_k_hits(cl->list, cl->length, dp, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, id, ci, &cmi, + aln_n?&p:NULL, sgap_rate, mm_gap); + for (k = cmi; k < ci; k++) cl->list[k].readID = ((uint32_t)(0x7fffffff)); +} + + +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) +{ + int64_t k, l, on = ol->length, m = 0, i, cn = cl->length, srt_n; + overlap_region *z; uint64_t t, ff, pid, *srt; char *str0 = NULL, *str1 = NULL; + kv_resize(uint64_t, *ii, ol->length); + for (k = 0, ii->n = 0; k < on; k++, ii->n++) { + ii->a[ii->n] = ol->list[k].shared_seed; + ii->a[ii->n] <<= 32; ii->a[ii->n] |= k; + refine_khits(&(ol->list[k]), cl, dp, sgap_rate); + } + radix_sort_bc64(ii->a, ii->a+ii->n); + + 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; + for (; i < cn && cl->list[i].readID == pid && cl->list[i].readID != ((uint32_t)(0x7fffffff)); i++) { cl->list[m] = cl->list[i]; - cl->list[m].readID = k; + cl->list[m].readID = (uint32_t)ii->a[k]; cl->list[m].cnt = 0; t = cl->list[m].self_offset; t <<= 32; t |= m; - kv_push(uint64_t, (*idx), t); + kv_push(uint64_t, (*ii), 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; + cl->length = cn = m; srt = ii->a+on; srt_n = ii->n-on; + radix_sort_bc64(srt, srt+srt_n); + if(basec) { + 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)) { + for (k = 1, l = 0; k <= srt_n; k++) { + if(k == cn || (srt[l]>>32) != (srt[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; + if(!k_hits_bcheck(rref, uref, ol, cl, khit, srt+l, k-l, qstr, 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); + cl->list[(uint32_t)srt[i]].cnt = ff; + // fprintf(stderr, "[M::%s::] pos::%u, cnt::%lu\n", __func__, + // cl->list[(uint32_t)srt[i]].self_offset, ff); + } + + l = k; + } + } + + +} + +void count_k_hits(All_reads *rref, const ul_idx_t *uref, char* qstr, UC_Read *buf, +overlap_region_alloc* ol, Candidates_list *cl, asg64_v* ii, uint64_t khit, uint64_t basec) +{ + int64_t k, l, on = ol->length, m = 0, i, cn = cl->length, srt_n; + overlap_region *z; uint64_t t, ff, pid, *srt; char *str0 = NULL, *str1 = NULL; + kv_resize(uint64_t, *ii, ol->length); + for (k = 0, ii->n = 0; k < on; k++, ii->n++) { + ii->a[ii->n] = ol->list[k].shared_seed; + ii->a[ii->n] <<= 32; ii->a[ii->n] |= k; + } + radix_sort_bc64(ii->a, ii->a+ii->n); + + 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; + for (; i < cn && cl->list[i].readID == pid; i++) { + cl->list[m] = cl->list[i]; + cl->list[m].readID = (uint32_t)ii->a[k]; + cl->list[m].cnt = 0; + t = cl->list[m].self_offset; t <<= 32; t |= m; + kv_push(uint64_t, (*ii), t); + m++; + } + } + cl->length = cn = m; srt = ii->a+on; srt_n = ii->n-on; + radix_sort_bc64(srt, srt+srt_n); + if(basec) { + resize_UC_Read(buf, (khit<<1)); str0 = buf->seq; str1 = buf->seq + khit; + } + + for (k = 1, l = 0; k <= srt_n; k++) { + if(k == cn || (srt[l]>>32) != (srt[k]>>32)) { + ff = k - l; + if(basec && ff > 1) { + if(!k_hits_bcheck(rref, uref, ol, cl, khit, srt+l, k-l, qstr, str0, str1)) ff = 1; + } + for (i = l; i < k; i++) { + cl->list[(uint32_t)srt[i]].cnt = ff; + // fprintf(stderr, "[M::%s::] pos::%u, cnt::%lu\n", __func__, + // cl->list[(uint32_t)srt[i]].self_offset, ff); } l = k; @@ -13187,14 +13389,18 @@ int64_t fusion_k_len, int64_t fusion_win_occ, int64_t *ch_s0, int64_t *ch_e0, in ///[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 ke, ks, 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++) { + ke = ch_a[ch_s].self_offset; + for (ch_e = ch_s; (ch_e < ch_n) && (ch_a[ch_e].self_offset < (*qe0)); ch_e++) { + ke = ch_a[ch_e].self_offset; + } + for (ke += fusion_k_len; (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])) { @@ -13215,7 +13421,7 @@ int64_t fusion_k_len, int64_t fusion_win_occ, int64_t *ch_s0, int64_t *ch_e0, in ///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) + for (k = 0; we < qe && we <= z->x_pos_e && 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; @@ -13250,7 +13456,11 @@ int64_t fusion_k_len, int64_t fusion_win_occ, int64_t *ch_s0, int64_t *ch_e0, in //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--) { + ks = ch_a[ch_e].self_offset; + for (ch_s = ch_e; (ch_s >= 0) && (ch_a[ch_s].self_offset > (*qs0)); ch_s--) { + ks = ch_a[ch_s].self_offset; + } + for (ks -= fusion_k_len; (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])) { @@ -13261,6 +13471,8 @@ int64_t fusion_k_len, int64_t fusion_win_occ, int64_t *ch_s0, int64_t *ch_e0, in p[0] = ch_s; } } + // fprintf(stderr, "-[M::%s::] utg%.6dl(%c), p[0]::%ld, p[1]::%ld, p[2]::%ld\n", + // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], p[0], p[1], p[2]); if(p[2] != -1) p[0] = p[2]; else if(p[1] != -1) p[0] = p[1]; if(p[0] != -1) { @@ -13271,7 +13483,7 @@ int64_t fusion_k_len, int64_t fusion_win_occ, int64_t *ch_s0, int64_t *ch_e0, in ///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) + for (k = 0; ws >= qs && ws >= z->x_pos_s && 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; @@ -13305,18 +13517,22 @@ int64_t fusion_k_len, int64_t fusion_win_occ, int64_t *ch_s0, int64_t *ch_e0, in } if((*ch_e0) == -1) { - for ((*ch_e0)=(*ch_s0);((*ch_e0)=(*qs0))&&(ch_a[(*ch_e0)].self_offset<(*qe0)); (*ch_e0)++); + 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)--); + 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)++; } + ///boundry + for (;((*ch_e0)=(*qs0))&&(ch_a[(*ch_e0)].self_offset<=(*qe0)); (*ch_e0)++); } 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; + if(qs < z->x_pos_s) qs = z->x_pos_s; + if(qe > z->x_pos_e+1) qe = z->x_pos_e+1; 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; @@ -13324,6 +13540,10 @@ int64_t cal_estimate_err(overlap_region *z, int64_t wl, int64_t qs, int64_t qe) we = z->w_list.a[k].x_end+1; os = MAX(qs, ws); oe = MIN(qe, we); ovlp = ((oe>os)? (oe-os):0); + // if(!ovlp) { + // fprintf(stderr, "\n[M::%s::] utg%.6dl(%c), q::[%ld, %ld), w::[%ld, %ld), z::::[%d, %d)\n", + // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], qs, qe, ws, we, z->x_pos_s, z->x_pos_e+1); + // } assert(ovlp); if(ovlp == (we-ws)) { tot += z->w_list.a[k].error; @@ -13334,73 +13554,367 @@ int64_t cal_estimate_err(overlap_region *z, int64_t wl, int64_t qs, int64_t qe) 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) +#define UC_Read_resize(v, s) do {\ + if ((v).size<(s)) {REALLOC((v).seq,(s));(v).size=(s);}\ + } while (0) + +///[s, e); [ps, pe) +inline char *retrieve_str_seq_exz(UC_Read *tu, int64_t s, int64_t l, +int64_t ps, int64_t pl, uint8_t rev, const ul_idx_t *uref, hpc_t *hpc_g, +All_reads *rref, int64_t id) +{ + if(!hpc_g) { + char *str; int64_t ss = s, sl = l; tu->length = l; + UC_Read_resize(*tu, sl); str = tu->seq; + if(s == ps) { + if(l <= pl) return tu->seq; + str = tu->seq + pl; ss = ps + pl; sl = l - pl; + } + if(uref) { + retrieve_u_seq(NULL, str, &(uref->ug->u.a[id]), rev, ss, sl, NULL); + } else if(rref) { + recover_UC_Read_sub_region(str, ss, sl, rev, rref, id); + } + return str; + } else { + return hpc_str(*hpc_g, id, rev) + s; + } +} + +void cal_exz_global(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez) +{ + int32_t bd, nword; + bd = (((thre)<<1)+1); nword = ((bd>>bitw)+(!!(bd&bitz))); + + if(nword <= 1) { + ed_band_cal_global_64_w_trace(pstr, pn, tstr, tn, thre, ez); + } else if(nword == 2) { + ed_band_cal_global_128_w_trace(pstr, pn, tstr, tn, thre, ez); + } else { + ed_band_cal_global_infi_w_trace(pstr, pn, tstr, tn, thre, &nword, ez); + } +} + +void cal_exz_extension_0(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez) +{ + int32_t bd, nword; + bd = (((thre)<<1)+1); nword = ((bd>>bitw)+(!!(bd&bitz))); + + if(nword <= 1) { + ed_band_cal_extension_64_0_w_trace(pstr, pn, tstr, tn, thre, ez); + } else if(nword == 2) { + ed_band_cal_extension_128_0_w_trace(pstr, pn, tstr, tn, thre, ez); + } else { + ed_band_cal_extension_infi_0_w_trace(pstr, pn, tstr, tn, thre, &nword, ez); + } +} + +void cal_exz_extension_1(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez) +{ + int32_t bd, nword; + bd = (((thre)<<1)+1); nword = ((bd>>bitw)+(!!(bd&bitz))); + + if(nword <= 1) { + ed_band_cal_extension_64_1_w_trace(pstr, pn, tstr, tn, thre, ez); + } else if(nword == 2) { + ed_band_cal_extension_128_1_w_trace(pstr, pn, tstr, tn, thre, ez); + } else { + ed_band_cal_extension_infi_1_w_trace(pstr, pn, tstr, tn, thre, &nword, ez); + } +} + +void cal_exz_semi(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t aux_beg, bit_extz_t *ez) +{ + int32_t bd, nword; + bd = (((thre)<<1)+1); nword = ((bd>>bitw)+(!!(bd&bitz))); + + if(nword <= 1) { + ed_band_cal_semi_64_w_absent_diag_trace(pstr, pn, tstr, tn, thre, aux_beg, ez); + } else if(nword == 2) { + ed_band_cal_semi_128_w_absent_diag_trace(pstr, pn, tstr, tn, thre, aux_beg, ez); + } else { + ed_band_cal_semi_infi_w_absent_diag_trace(pstr, pn, tstr, tn, thre, aux_beg, &nword, ez); + } +} + +int64_t cal_exz_infi_adv(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 *pts, int64_t *pte, int64_t thre, int64_t *pthre, int64_t q_tot_l, int64_t mode) { 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; - } - } + int64_t aux_beg = 0, ql, tl, t_tot_l = -1, dd; + char *q_string, *t_string; int32_t rev = z->y_pos_strand, id = z->y_id; + ql = qe - qs; tl = te - ts; dd = MAX(ql, tl); + if(hpc_g) t_tot_l = hpc_len(*hpc_g, id); + else if(uref) t_tot_l = uref->ug->u.a[id].len; + else t_tot_l = Get_READ_LENGTH((*rref), id); + + if(mode == 3) { + update_semi_coord(uref, hpc_g, rref, z, qs, qe, ((thre>dd)?dd:thre), &ts, &te, &aux_beg); + } else if(mode == 1 || mode == 2) { + adjust_ext_offset(&qs, &qe, &ts, &te, q_tot_l, t_tot_l, ((thre>dd)?dd:thre), mode); } - fprintf(stderr, ", err::%d, thre::%d\n", INT32_MAX, exz->thre); + + if((qe > qs) && (te > ts) && (ts != -1) && (te != -1)) { + ql = qe - qs; tl = te - ts; + dd = MAX(ql, tl); + if(thre > dd) thre = dd; + if(thre <= (*pthre)) return 0; + (*pthre) = thre; + + q_string = qstr + qs; + t_string = retrieve_str_seq_exz(tu, ts, tl, (*pts), (*pte)-(*pts), rev, uref, hpc_g, rref, id); + (*pts) = ts; (*pte) = te; + + if(mode == 0) { //global + cal_exz_global(t_string, tl, q_string, ql, thre, exz); + } else if(mode == 1) {///forward extension + cal_exz_extension_0(t_string, tl, q_string, ql, thre, exz); + } else if(mode == 2) {///backward extension + cal_exz_extension_1(t_string, tl, q_string, ql, thre, exz); + } else if(mode == 3) {//semi-global + cal_exz_semi(t_string, tl, q_string, ql, thre, aux_beg, exz); + } + + if(is_align(*exz)) { + // if(exz->err < 0) { + // fprintf(stderr, "\n[M::%s::ql::%ld] qs::%ld, qe::%ld, ts::%ld, te::%ld, mode::%ld, err::%d, thre::%d\n", + // __func__, ql, qs, qe, ts, te, mode, exz->err, exz->thre); + // fprintf(stderr, "[M::%s::] pstr::%.*s\n", __func__, (int32_t)tu->length, tu->seq); + // fprintf(stderr, "[M::%s::] tstr::%.*s\n", __func__, (int32_t)(qe-qs), qstr+qs); + // } + return 1; + } + return 0; + } + return 0; } -int64_t sub_base_aln(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, uint64_t pre_e, +int64_t cal_exact_exz(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 *pts, int64_t *pte, int64_t q_tot_l, int64_t mode) +{ + clear_align(*exz); exz->thre = 0; + int64_t ql, tl, t_tot_l = -1; + 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); + else if(uref) t_tot_l = uref->ug->u.a[id].len; + else t_tot_l = Get_READ_LENGTH((*rref), id); + + if(mode == 3) {//semi + ts = (qs - z->x_pos_s) + z->y_pos_s; ts += y_start_offset(qs, &(z->f_cigar)); + te = ts + ql; + } else if(mode == 1) {///forward extension + te = ts + ql; + } else if(mode == 2) {///backward extension + ts = te - ql; + } + if(ts < 0) ts = 0; + if(ts > t_tot_l) ts = t_tot_l; + if(te > t_tot_l) te = t_tot_l; + ql = qe - qs; tl = te - ts; + if(ql != tl) return 0; + + q_string = qstr + qs; + t_string = retrieve_str_seq_exz(tu, ts, tl, (*pts), (*pte)-(*pts), rev, uref, hpc_g, rref, id); + (*pts) = ts; (*pte) = te; + + if(memcmp(q_string, t_string, ql)) return 0; + exz->err = 0; push_trace(&(exz->cigar), 0, ql); + exz->pl = tl; exz->ps = 0; exz->pe = tl; + exz->tl = ql; exz->ts = 0; exz->te = ql; + return 1; +} + +int64_t hc_aln_exz_adv(overlap_region *z, 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, int64_t maxl, int64_t maxe, int64_t force_l, +int64_t estimate_err) +{ + clear_align(*exz); exz->thre = 0; + if(((ts == -1) && (te == -1))) mode = 3;///set to semi-global + int64_t thre, ql = qe - qs, thre0, pts = -1, pte = -1, pthre = -1; + if(ql <= 0) return 0; + if(estimate_err < 0) { + if(ql > wl) estimate_err = cal_estimate_err(z, wl, qs, qe); + else estimate_err = ql*e_rate; + } + + + // fprintf(stderr, "[M::%s::ql::%ld] qs::[%ld, %ld), ts::[%ld, %ld), mode::%ld, est_err::%ld, e_rate::%f, maxe::%ld", + // __func__, ql, qs, qe, ts, te, mode, estimate_err, e_rate, maxe); + if(ql <= 16) { + if(cal_exact_exz(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, &pts, &pte, q_tot, mode)) { + // fprintf(stderr, ", err::%d, thre::%d, scale::0(+)\n", exz->err, exz->thre); + return 1; + } + } + + if(ql <= maxl && (estimate_err>>1) <= maxe) { + thre = scale_ed_thre(estimate_err, maxe); + if(cal_exz_infi_adv(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, &pts, &pte, thre, &pthre, q_tot, mode)) { + // fprintf(stderr, ", err::%d, thre::%d, scale::%ld(+)\n", exz->err, exz->thre, thre); + return 1; + } + + thre0 = thre; thre = ql*e_rate; thre = scale_ed_thre(thre, maxe); + if(thre > thre0) { + if(cal_exz_infi_adv(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, &pts, &pte, thre, &pthre, q_tot, mode)) { + // fprintf(stderr, ", err::%d, thre::%d, scale::%ld(-)\n", exz->err, exz->thre, thre); + return 1; + } + } + + thre0 = thre; thre <<= 1; thre = scale_ed_thre(thre, maxe); + if(thre > thre0) { + if(cal_exz_infi_adv(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, &pts, &pte, thre, &pthre, q_tot, mode)) { + // fprintf(stderr, ", err::%d, thre::%d, scale::%ld(-)\n", exz->err, exz->thre, thre); + return 1; + } + } + + thre0 = thre; thre = ql*0.51; thre = scale_ed_thre(thre, maxe); + if(thre > thre0) { + if(cal_exz_infi_adv(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, &pts, &pte, thre, &pthre, q_tot, mode)) { + // fprintf(stderr, ", err::%d, thre::%d, scale::%ld(*)\n", exz->err, exz->thre, thre); + return 1; + } + } + + if(ql <= force_l) { + thre = maxe; + if(cal_exz_infi_adv(z, uref, hpc_g, rref, exz, qstr, tu, qs, qe, ts, te, &pts, &pte, thre, &pthre, q_tot, mode)) { + // fprintf(stderr, ", err::%d, thre::%d, scale::%ld(*)\n", exz->err, exz->thre, thre); + return 1; + } + } + } + // fprintf(stderr, ", err::%d, thre::%d\n", INT32_MAX, exz->thre); + // if(mode == 0) { + // fprintf(stderr, "[M::%s::] pstr::%.*s\n", __func__, (int32_t)tu->length, tu->seq); + // fprintf(stderr, "[M::%s::] tstr::%.*s\n", __func__, (int32_t)(qe-qs), qstr+qs); + // } + return 0; + +} + +void prt_k_mer_hit(k_mer_hit *ch_a, int64_t ch_n) +{ + int64_t k; + for (k = 0; k < ch_n; k++) { + fprintf(stderr, "[M::%s::k->%ld] q_pos::%u, t_pos::%u, cnt::%u, cov::%u\n", + __func__, k, ch_a[k].self_offset, ch_a[k].offset, ch_a[k].cnt, ch_a[k].readID); + } + +} + +void debug_iter_k_mer_hit(k_mer_hit *ch_a, int64_t ch_n, uint64_t s, uint64_t e, int64_t ibeg, int64_t iend) +{ + int64_t i, beg = -1, end = -1; + for (i = 0; i < ch_n; i++) { + if((ch_a[i].self_offset >= s) && (ch_a[i].self_offset < e)) { + if(beg == -1) beg = i; + end = i+1; + } + } + assert(ibeg==beg && iend==end); +} + +int64_t chain_aln(overlap_region *z, Chain_Data *dp, 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, +int64_t min_chain_aln, uint64_t rid) +{ + int64_t *m, mn = 0, k; int64_t q[2], t[2], is_chain_aln = 1; + if(mode == 0) { + ///wrong + // assert((ch_n >= 2) && (ch_a[0].self_offset == qs) && (ch_a[0].offset == ts) && + // (ch_a[ch_n-1].self_offset == qe) && (ch_a[ch_n-1].offset == te)); + //does not work with CNS alignment; it also could not work here + if(ch_n == 2 && ((qe-qs)+128) < min_chain_aln) is_chain_aln = 0; + } else if(mode == 1) { + // assert((ch_n >= 1) && (ch_a[0].self_offset == qs) && (ch_a[0].offset == ts)); + //does not work with CNS alignment; it also could not work here + if(ch_n == 1 && ((qe-qs)+128) < min_chain_aln) is_chain_aln = 0; + } else if(mode == 2) { + // assert((ch_n >= 1) && (ch_a[ch_n-1].self_offset == qe) && (ch_a[ch_n-1].offset == te)); + //does not work with CNS alignment; it also could not work here + if(ch_n == 1 && ((qe-qs)+128) < min_chain_aln) is_chain_aln = 0; + } + + if(is_chain_aln) { + mn = lchain_refine(ch_a, ch_n, NULL, dp, 50, 5000, 512, 16); m = dp->tmp; + q[0] = qs; t[0] = ts; + for (k = 0; k < mn; k++) { + q[1] = ch_a[m[k]].self_offset; t[1] = ch_a[m[k]].offset; + if(q[1] > q[0]) { + 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(!hc_aln_exz_adv(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, q_tot, e_rate, MAX_SIN_L, MAX_SIN_E, FORCE_SIN_L, -1)) { + + } + } + q[0] = q[1]; t[0] = t[1]; + } + q[1] = qe; t[1] = te; + if(q[1] > q[0]) { + 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(!hc_aln_exz_adv(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, q_tot, e_rate, MAX_SIN_L, MAX_SIN_E, FORCE_SIN_L, -1)) { + + } + } + + } else { + + } + + // { + // q[0] = qs; t[0] = ts; q[1] = qe; t[1] = te; + // if(q[1] > q[0]) { + // hc_aln_exz_adv(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, q_tot, e_rate, MAX_SIN_L, MAX_SIN_E); + // } + // } + return 0; +} + +int64_t sub_base_aln(overlap_region *z, Chain_Data *dp, 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(i<0) i = 0; ibeg=iend=-1; + for (; 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); + // debug_iter_k_mer_hit(ch_a, ch_n, s, e, ibeg, iend); + // fprintf(stderr, "***[M::%s::rid->%lu] utg%.6dl(%c), s::%lu, e::%lu, z::[%u, %u), ibeg::%ld, iend::%ld\n", + // __func__, rid, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], s, e, z->x_pos_s, z->x_pos_e+1, ibeg, iend); ch_i = i; - if(ibeg > 0 && iend > 0 && iend > ibeg) {///find some anchors[ibeg, iend) + if(ibeg >= 0 && iend >= 0 && iend > ibeg) {///find some anchors[ibeg, iend) for (i = l = ibeg; i <= iend; i++) { + // fprintf(stderr, "\n[M::%s::i->%ld] q::%u, t::%u, cnt::%u, readID::%u\n", __func__, + // i, ch_a[i].self_offset, ch_a[i].offset, ch_a[i].cnt, ch_a[i].readID); 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])) { @@ -13425,23 +13939,36 @@ int64_t ch_i, uint64_t rid) 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); + if((mode == 0) && is_alnw(ch_a[l]) && is_alnw(ch_a[i]) + && (ch_a[l].strand == 0) && (ch_a[i].strand == 1)) { + ; + } else { + // 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(mode == 1 || mode == 2) { + tuning_ext_offset(z, ch_a, ch_n, ql, tl, wl, MAX_SIN_L, 4, &ch_s, &ch_e, &q[0], &q[1], &t[0], &t[1], &mode); + } else if(ch_e >= 0) {//global + ch_e++; + } + if(!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)) { + if(ch_s < 0) ch_s = ibeg; if(ch_e < 0) ch_e = iend; + assert(ch_e > ch_s); + // chain_aln(z, dp, 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, MAX_CNS_L, rid); + } } - - 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 - + // fprintf(stderr, "\n***[M::%s::rid->%lu] utg%.6dl(%c), s::%lu, e::%lu, z::[%u, %u), ibeg::%ld, iend::%ld\n", + // __func__, rid, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], s, e, z->x_pos_s, z->x_pos_e+1, ibeg, iend); } 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, +void cigar_gen_by_chain(overlap_region *z, Chain_Data *dp, 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) { @@ -13450,29 +13977,247 @@ int64_t ql, uint64_t rid) 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 = 0; i < wn; i++) z->w_list.a[i].clen = 0;///clean cigar + // if(on > 1) { + // fprintf(stderr, "[M::%s::] rid::%lu, on::%ld\n", __func__, rid, on); + // } + // if(z->y_id == 126) prt_k_mer_hit(ch_a, ch_n); 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); + ch_i = sub_base_aln(z, dp, 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; } } + +int64_t adjust_base_coordinates(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, +ul_ov_t *res, int64_t wl, int64_t ql, int64_t tl, int64_t ch_i) +{ + res->ts = res->te = (uint32_t)-1; + res->sec = 3; res->qn = res->tn = (uint32_t)-1; //semi-global + if(ch_n == 0) return ch_i; + int64_t i = ch_i, ibeg, iend; uint64_t s = res->qs, e = res->qe; + for (; i >= 0 && ch_a[i].self_offset >= s; i--); if(i < 0) i = 0; + + for (; i < ch_n && ch_a[i].self_offset < s; i++); + if((i > 0) && ((ch_a[i].self_offset > s) || (i >= ch_n))) i--; ibeg = i; ///if ch_a[i].self_offset == s, do nothing + for (; i < ch_n && ch_a[i].self_offset < e; i++); iend = i; + ch_i = i;///ch_i must be here + ///ibeg might be < 0, iend might be == ch_n + ///1. [s, e) contain anchors + ///2. anchors contain [s, e) + + if(ibeg >= 0) { + res->qs = ch_a[ibeg].self_offset; + res->ts = ch_a[ibeg].offset; + res->qn = ibeg; + } else {//extension to left + res->qs = 0; res->qn = (uint32_t)-1; + } + + if(iend < ch_n) { + res->qe = ch_a[iend].self_offset; + res->te = ch_a[iend].offset; + res->tn = iend; + } else {//extension to right + res->qe = ql; res->tn = ch_n; + } + + return ch_i; +} + +inline int64_t translate_double_mode(uint64_t double_mode, uint64_t is_backward) +{ + if(double_mode == 0) return 0; + if(double_mode == 4) return 3; + if(double_mode == 1 || double_mode == 2) return double_mode; +} + +int64_t fusion_chain_ovlp(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, ul_ov_t *ov, int64_t on, uint64_t wl, int64_t ql, int64_t tl) +{ + int64_t i, srt, ch_i, m, os, oe, ovlp; ul_ov_t *p; + for (i = ch_i = 0, srt = 1; i < on; i++) { + ov[i].sec = 6;///do not know the aln type + ch_i = adjust_base_coordinates(z, ch_a, ch_n, &(ov[i]), wl, ql, tl, ch_i); + if(i > 0 && ov[i].qs < ov[i-1].qe) srt = 0; + } + if(on <= 1) return on; + + if(!srt) radix_sort_uov_srt_qs(ov, ov+on); + for (i = m = 1; i < on; i++) { + p = &(ov[m-1]); + os = MAX(p->qs, ov[i].qs); + oe = MIN(p->qe, ov[i].qe); + ovlp = oe - os; + if(ovlp >= 0) {//merge + p->qe = MAX(p->qe, ov[i].qe); + p->te = MAX(p->te, ov[i].te); + p->tn = MAX(p->tn, ov[i].tn); + } else {//new + ov[m++] = ov[i]; + } + } + + on = m; + return on; +} + +// void chain_win_aln(overlap_region *z, Chain_Data *dp, Candidates_list *cl, int64_t qs, int64_t qe, +// int64_t ts, int64_t te, int64_t ql, int64_t tl, int64_t wl, int64_t mode, bit_extz_t *exz) +// { + +// } + +int64_t ovlp_base_aln_all(overlap_region *z, Chain_Data *dp, Candidates_list *cl, int64_t ch_tot_beg, +int64_t ch_tot_n, int64_t soff, int64_t eoff, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, +char* qstr, UC_Read *tu, ul_ov_t *ov, int64_t ql, int64_t tl, int64_t wl, bit_extz_t *exz, double e_rate) +{ + k_mer_hit *ch_a = cl->list + ch_tot_beg; + int64_t ch_n = ch_tot_n, ibeg, iend, i, l, mode, q[2], t[2], is_done; + ibeg = soff; iend = eoff; + for (l = ibeg, i = ibeg + 1; i <= iend; i++) { + l = i - 1; + q[0] = q[1] = t[0] = t[1] = mode = -1; is_done = 0; + if(l >= 0) { + q[0] = ch_a[l].self_offset; t[0] = ch_a[l].offset; + } else { + q[0] = ov->qs; + } + + if(i < ch_n) { + q[1] = ch_a[i].self_offset; t[1] = ch_a[i].offset; + } else { + q[1] = ov->qe; + } + + 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;///no primary hit within [ibeg, iend] + } + assert(mode != 3); + if(mode == 1 || mode == 2) adjust_ext_offset(&(q[0]), &(q[1]), &(t[0]), &(t[1]), ql, tl, 0, mode); + ///at cns chain, the base alignment fails; there is no anchor between soff and eoff + if((eoff-soff<=1) && (((q[1]-q[0])>>1) < MAX_CNS_E)) { + is_done = 0; + } 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); + } + if(!is_done) { + // chain_win_aln(z, dp, cl, q[0], q[1], t[0], t[1], ql, tl, wl, exz); + ch_a = cl->list + ch_tot_beg;//update ch_a + } + } + return 0; +} + +void ovlp_base_aln(overlap_region *z, Chain_Data *dp, Candidates_list *cl, int64_t ch_beg, int64_t ch_n, +ul_ov_t *ov, int64_t wl, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, +bit_extz_t *exz, bit_extz_t *exz1, double e_rate, int64_t ql, int64_t tl, uint64_t rid) +{ + int64_t ibeg, iend, i, l, mode, q[2], t[2], is_done; + k_mer_hit *ch_a = cl->list + ch_beg; + if(ov->qn == ((uint32_t)-1)) ibeg = -1; + else ibeg = ov->qn; + iend = ov->tn; + assert(iend>=ibeg+1); + exz1->err = 0; exz1->thre = 0; + exz1->ps = exz1->pe = (uint32_t)-1; + + for (l = ibeg, i = ibeg + 1; i <= iend; i++) { + if(i == iend || is_pri_aln(ch_a[i])) { + q[0] = q[1] = t[0] = t[1] = mode = -1; is_done = 0; + if(l >= 0) { + q[0] = ch_a[l].self_offset; t[0] = ch_a[l].offset; + } else { + q[0] = ov->qs; + } + + if(i < ch_n) { + q[1] = ch_a[i].self_offset; t[1] = ch_a[i].offset; + } else { + q[1] = ov->qe; + } + + 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;///no primary hit within [ibeg, iend] + } + + 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; + } else if(mode != 3) { + if(mode == 1 || mode == 2) adjust_ext_offset(&(q[0]), &(q[1]), &(t[0]), &(t[1]), ql, tl, 0, 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); + } + + if(!is_done) {///postprocess + is_done = ovlp_base_aln_all(z, dp, cl, ch_beg, ch_n, l, i, uref, hpc_g, rref, qstr, tu, ov, ql, tl, wl, exz, e_rate); + ch_a = cl->list + ch_beg;///update ch_a + } + l = i; + } + } + +} + +void cigar_gen_by_chain_adv(overlap_region *z, Candidates_list *cl, int64_t ch_beg, 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, bit_extz_t *exz1, double e_rate, int64_t ql, uint64_t rid) +{ + if(on <= 0) return; + int64_t i, tl, id = z->y_id; k_mer_hit *ch_a; Chain_Data *dp; + 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); + dp = &(cl->chainDP); ch_a = cl->list + ch_beg; + on = fusion_chain_ovlp(z, ch_a, ch_n, ov, on, wl, ql, tl); + + for (i = 0; i < on; i++) { + // ch_a = cl->list + ch_beg; + ovlp_base_aln(z, dp, cl, ch_beg, ch_n, &(ov[i]), wl, uref, hpc_g, rref, qstr, tu, exz, exz1, e_rate, ql, tl, rid); + } + + // for (i = 0; i < wn; i++) z->w_list.a[i].clen = 0;///clean cigar + // if(on > 1) { + // fprintf(stderr, "[M::%s::] rid::%lu, on::%ld\n", __func__, rid, on); + // } + // if(z->y_id == 126) prt_k_mer_hit(ch_a, ch_n); + // for (i = ch_i = 0; i < on; i++) { + // assert((i<=0)||(ov[i].qs > ov[i-1].qe)); + // ov[i].sec = 16;///do not know the aln type + // ch_i = sub_base_aln(z, dp, 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, +const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, bit_extz_t *exz1, 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); + int64_t k, l, ch_n, a_n = aln->n; overlap_region *z; + // 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); 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); + z = &(ol->list[aln->a[l].qn]); assert(z->align_length == l); + ch_n = gen_cns_chain(z, cl, iidx, max_lgap, e_rate, 0); 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); + // ch_a = cl->list + cl->length; + // cigar_gen_by_chain(z, &(cl->chainDP), ch_a, ch_n, aln->a+l, k-l, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, rid); + cigar_gen_by_chain_adv(z, cl, cl->length, ch_n, aln->a+l, k-l, wl, uref, hpc_g, rref, qstr, tu, exz, exz1, 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); @@ -13493,7 +14238,7 @@ inline uint32_t ovlp_win_check(overlap_region *z, uint32_t id0, uint32_t id1, in 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 + if((ts1 < te0) && (qe0 == qs1)) {//has overlap in y dm = dq>=dr?dr:dq; if((dd > (dm*small_bw_rate)) && (dd > min_small_bw)) return 0; } else { @@ -13503,10 +14248,15 @@ inline uint32_t ovlp_win_check(overlap_region *z, uint32_t id0, uint32_t id1, in } 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, int64_t cov, int64_t max_lgap, double sgap_rate, int64_t sgap) +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, sc, max_sc, p, long_sc, long_idx; overlap_region *z; + // fprintf(stderr, "\n[M::%s::] n_skip::%ld, win_n::%ld\n", __func__, n_skip, win_n); + // for (k = 0; k < win_n; k++) { + // fprintf(stderr, "[M::%s::] win_a[%ld]::%lu\n", __func__, k, win_a[k]); + // } if(n_skip < win_n) { long_sc = long_idx = -1; for (k = 0; k < n_skip; k++) { @@ -13535,7 +14285,7 @@ uint64_t *win_a, int64_t win_n, uint64_t *dp, int64_t n_skip, int64_t wl, int64_ long_sc = max_sc; long_idx = k; } } - + // fprintf(stderr, "[M::%s::] long_idx::%lu, long_sc::%lu\n", __func__, long_idx, long_sc); k = long_idx; while (k >= 0) { win_a[k] |= ((uint64_t)0x8000000000000000); @@ -13544,9 +14294,10 @@ uint64_t *win_a, int64_t win_n, uint64_t *dp, int64_t n_skip, int64_t wl, int64_ for (k = 0, m = 0; k < win_n; k++) { if(!(win_a[k]&((uint64_t)0x8000000000000000))) continue; - win_a[m++] = win_a[k]; + win_a[m++] = (win_a[k]<<1)>>1; } win_n = m; + // fprintf(stderr, "[M::%s::] win_n::%lu\n", __func__, win_n); } for (k = 0, cov = -cov; k < win_n; k++) { @@ -13554,6 +14305,10 @@ uint64_t *win_a, int64_t win_n, uint64_t *dp, int64_t n_skip, int64_t wl, int64_ 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); + // if(wid >= z->w_list.n) { + // fprintf(stderr, "[M::%s::] ws::%ld, wl::%ld, wid::%ld, wn::%ld\n", + // __func__, ws, wl, wid, ((int64_t)z->w_list.n)); + // } 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 @@ -13681,6 +14436,7 @@ kv_ul_ov_t *aln, uint64_t rid, int64_t max_lgap, double sgap_rate) } if((end > beg) && (end - beg > wl) && (old_dp >= 2) ) { + // fprintf(stderr, "+++[M::%s::] beg::%ld, end::%ld\n", __func__, beg, end); 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, max_lgap, sgap_rate, SGAP); @@ -13688,17 +14444,18 @@ kv_ul_ov_t *aln, uint64_t rid, int64_t max_lgap, double sgap_rate) beg = end; } - for (i = srt_n = 0; i < aln->n; i++) { + for (i = 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; + // fprintf(stderr, "+++[M::%s::] idx->n::%ld\n", __func__, idx->n); return; } 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, + bit_extz_t *exz1, haplotype_evdience_alloc* hap, kvec_t_u64_warp* v_idx, 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; @@ -13753,7 +14510,8 @@ void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur if(ol->length <= 1) return; ///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)); + // fprintf(stderr, "\n[M::%s] iidx_n::%ld\n", __func__, (int64_t)iidx.n); + ul_gap_filling_adv(ol, cl, aln, wl, uref, NULL, NULL, qu->seq, tu, exz, exz1, &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; diff --git a/Correct.h b/Correct.h index da8c7f9..ac23396 100644 --- a/Correct.h +++ b/Correct.h @@ -20,7 +20,7 @@ #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 MAX_LGAP(ql) ((((ql)*0.2)<256)?((ql)*0.2):256) #define WINDOW_MAX_SIZE (WINDOW + (int)(1.0 / HA_MIN_OV_DIFF) + 3) // TODO: why 1/max_ov_diff? @@ -1136,7 +1136,7 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref double max_ov_diff_ec, long long winLen, void *km); 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, + bit_extz_t *exz1, haplotype_evdience_alloc* hap, kvec_t_u64_warp* v_idx, 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, diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 2b6c1e2..ec4d4d5 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -185,6 +185,60 @@ int ovlp_chain_gen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int6 return 1; } +int ovlp_chain_qgen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int64_t yl, int64_t apend_be, k_mer_hit* hit, int64_t n_hit) +{ + if (ol->length + 1 > ol->size) { + uint64_t sl = ol->size; + ol->size = ol->length + 1; + kroundup64(ol->size); + REALLOC(ol->list, ol->size); + /// need to set new space to be 0 + memset(ol->list + sl, 0, sizeof(overlap_region)*(ol->size - sl)); + } + + if ((ol->length!=0) && (ol->list[ol->length-1].y_id==t->y_id)) { + if((ol->list[ol->length-1].shared_seed > t->shared_seed) || + ((ol->list[ol->length-1].shared_seed == t->shared_seed) && + (ol->list[ol->length-1].overlapLen <= t->overlapLen))) { + return 0; + } else { + ol->length--; + } + } + + int64_t xr, yr; + if(t->x_pos_s <= t->y_pos_s) { + t->y_pos_s -= t->x_pos_s; t->x_pos_s = 0; + } else { + t->x_pos_s -= t->y_pos_s; t->y_pos_s = 0; + } + + xr = xl-t->x_pos_e-1; yr = yl-t->y_pos_e-1; + if(xr <= yr) { + t->x_pos_e = xl-1; t->y_pos_e += xr; + } else { + t->y_pos_e = yl-1; t->x_pos_e += yr; + } + + overlap_region *o = &(ol->list[ol->length++]); + o->shared_seed = t->shared_seed; + o->align_length = 0; + o->is_match = 0; + o->non_homopolymer_errors = 0; + o->strong = 0; + o->x_id = t->x_id; + o->y_id = t->y_id; + o->x_pos_strand = 0;///always 0 + o->y_pos_strand = t->x_pos_strand; + + o->x_pos_e = t->x_pos_e; o->x_pos_s = t->x_pos_s; + o->y_pos_e = t->y_pos_e; o->y_pos_s = t->y_pos_s; + ///debug + // debug_cigar(&(t->f_cigar), o, apend_be, hit, n_hit); + + return 1; +} + int ovlp_chain_gen_fcigar(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int64_t yl, int64_t apend_be, k_mer_hit* hit, int64_t n_hit) { if (ol->length + 1 > ol->size) { @@ -1359,6 +1413,58 @@ int32_t lchain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, double bw_thres) return n_a; } + +int32_t lchain_qcheck(k_mer_hit *a, int32_t n_a, Chain_Data *dp, double bw_thres) +{ + int32_t i, tot_g = 0, sc, dg, dq, dr, dd, span; + double bw_pen; + if (n_a == 0) return -1; + if (n_a > 1) { + if ((a[0].self_offset >= a[n_a-1].self_offset)||(a[0].offset >= a[n_a-1].offset)) return -1; + dq = (int32_t)a[n_a-1].self_offset - (int32_t)a[0].self_offset; + dr = (int32_t)a[n_a-1].offset - (int32_t)a[0].offset; + dd = ((dq>=dr)? (dq-dr): (dr-dq));//gap + dg = ((dq>=dr)? (dr): (dq));///len + if (dg == 0 || dd > (dg*bw_thres)) return -1; + } + + for (i = 1; i < n_a; ++i) {///a[] is sorted by self_offset + if(a[i-1].self_offset >= a[i].self_offset) break; + if(a[i-1].offset >= a[i].offset) break; + } + if (i < n_a) return -1; + + bw_pen = 1.0 / bw_thres; + dp->score[0] = normal_w((a[0].cnt&(0xffu)), (a[0].cnt>>8)); dp->pre[0] = -1; + for (i = 1; i < n_a; ++i) { + dq = (int32_t)a[i].self_offset - (int32_t)a[i-1].self_offset; + dr = (int32_t)a[i].offset - (int32_t)a[i-1].offset; + dd = ((dq>=dr)? (dq-dr): (dr-dq));//gap + dg = ((dq>=dr)? (dr): (dq));///len + if(dg == 0) break; + + tot_g += dd; + if (dd > THRESHOLD_MAX_SIZE && dd > (dg*bw_thres)) break; + span = a[i].cnt&(0xffu); + sc = dg < span? dg : span; + sc = normal_w(sc, ((int32_t)(a[i].cnt>>8))); + sc -= (int32_t)((((double)dd)/((double)dg))*bw_pen*((double)sc));///bw_pen is 20 for HiFi + + dp->score[i] = dp->score[i-1] + sc; + dp->pre[i] = i - 1; + } + if (i < n_a) return -1; + + if(n_a > 1) { + dq = (int32_t)a[n_a-1].self_offset - (int32_t)a[0].self_offset; + dr = (int32_t)a[n_a-1].offset - (int32_t)a[0].offset; + dg = ((dq>=dr)? (dr): (dq));///len + dd = tot_g;///gap + if (dd > (dg*bw_thres)) return -1; + } + return n_a; +} + inline int32_t cal_bw(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, int64_t sf_l, int64_t ot_l) { ///ai is the suffix of aj @@ -1624,6 +1730,180 @@ uint64_t lchain_dp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, ov return cL; } +uint64_t lchain_qdp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, overlap_region* res, + 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 xl, int64_t yl, int64_t quick_check) +{ + int64_t *p, *t, max_f, n_skip, st, max_j, end_j, sc, msc, msc_i, bw, max_ii, ovl, movl; + int32_t *f, max, tmp; int64_t i, j, ret, cL = 0; + resize_Chain_Data(dp, a_n, NULL); + t = dp->tmp; f = dp->score; p = dp->pre; + bw = ((xl < yl)?xl:yl); bw *= bw_rate; + msc = msc_i = -1; movl = INT32_MAX; + + if(quick_check) { + ret = lchain_qcheck(a, a_n, dp, bw_rate); + if (ret > 0) { + a_n = ret; msc_i = a_n-1; msc = f[msc_i]; + goto skip_ldp; + } + } + + memset(t, 0, (a_n*sizeof((*t)))); + for (i = st = 0, max_ii = -1; i < a_n; ++i) { + max_f = a[i].cnt&(0xffu); + n_skip = 0; max_j = end_j = -1; + if ((i-st) > max_iter) st = i-max_iter; + + for (j = i - 1; j >= st; --j) { + sc = comput_sc_ch(&a[i], &a[j], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl); + if (sc == INT32_MIN) continue; + sc += f[j]; + if (sc > max_f) { + max_f = sc, max_j = j; + if (n_skip > 0) --n_skip; + } else if (t[j] == (int32_t)i) { + if (++n_skip > max_skip) + break; + } + if (p[j] >= 0) t[p[j]] = i; + } + end_j = j; + + if (max_ii < 0 || ((int64_t)a[i].self_offset) - ((int64_t)a[max_ii].self_offset) > max_dis) { + max = INT32_MIN; max_ii = -1; + for (j = i - 1; (j >= st) && ((((int64_t)a[i].self_offset)-((int64_t)a[j].self_offset))<=max_dis); --j) { + if (max < f[j]) { + max = f[j], max_ii = j; + } + } + } + + if (max_ii >= 0 && max_ii < end_j) {///just have a try with a[i]<->a[max_ii] + tmp = comput_sc_ch(&a[i], &a[max_ii], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl); + if (tmp != INT32_MIN && max_f < tmp + f[max_ii]) + max_f = tmp + f[max_ii], max_j = max_ii; + } + f[i] = max_f; p[i] = max_j; + if ((max_ii < 0) || (((((int64_t)a[i].self_offset)-((int64_t)a[max_ii].self_offset))<=max_dis) && (f[max_ii]= msc) { + ovl = get_chainLen(a[i].self_offset, a[i].self_offset, xl, a[i].offset, a[i].offset, yl); + if(f[i] > msc || ovl < movl) { + msc = f[i]; msc_i = i; movl = ovl; + } + } + } + + skip_ldp: + ///a[] has been sorted by self_offset + i = msc_i; + res->x_pos_s = res->x_pos_e = a[i].self_offset; + res->y_pos_s = res->y_pos_e = a[i].offset; + res->shared_seed = msc; + + cL = 0; + while (i >= 0) { + t[cL++] = i; msc_i = i; i = p[i]; + } + + res->x_pos_s = a[t[cL-1]].self_offset; + res->y_pos_s = a[t[cL-1]].offset; + 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; +} + + +uint64_t lchain_refine(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, + int64_t max_skip, int64_t max_iter, int64_t max_dis, int64_t long_gap) +{ + if(a_n <= 0) return 0; + int64_t *p, *t, max_f, n_skip, st, max_j, sc, msc, msc_i, dq, dr, dd; + int32_t *f; int64_t i, j, cL = 0; + resize_Chain_Data(dp, a_n, NULL); + t = dp->tmp; f = dp->score; p = dp->pre; msc = msc_i = -1; + + for (i = 1, f[0] = 0, p[0] = -1, msc_i = a_n - 1; i < a_n; i++) { + j = i-1; + dq = (int64_t)(a[i].self_offset) - (int64_t)(a[j].self_offset); + dr = (int64_t)(a[i].offset) - (int64_t)(a[j].offset); + dd = dr > dq? dr - dq : dq - dr;//gap + if(dd <= long_gap || dq > max_dis) { + p[i] = i - 1; f[i] = i; + } else { + break; + } + } + if(i >= a_n) goto ss_kip; + + + memset(t, 0, (a_n*sizeof((*t)))); + f[0] = 0; p[0] = -1; + + for (i = 1, st = 0; i < a_n; ++i) { + max_f = INT32_MIN; n_skip = 0; max_j = -1; + if ((i-st) > max_iter) st = i-max_iter; + ///i-1 + j = i - 1; + dq = (int64_t)(a[i].self_offset) - (int64_t)(a[j].self_offset); + dr = (int64_t)(a[i].offset) - (int64_t)(a[j].offset); + dd = dr > dq? dr - dq : dq - dr;//gap + if(dd <= long_gap) dd = 0; + sc = f[j] - dd; + if (sc > max_f) { + max_f = sc, max_j = j; + if (n_skip > 0) --n_skip; + } else if (t[j] == (int32_t)i) { + if (++n_skip > max_skip) + break; + } + if (p[j] >= 0) t[p[j]] = i; + ///[st, i-2] + for (--j; (j >= st) && (a[i].self_offset <= (max_dis + a[j].self_offset)); --j) { + dq = (int64_t)(a[i].self_offset) - (int64_t)(a[j].self_offset); + dr = (int64_t)(a[i].offset) - (int64_t)(a[j].offset); + dd = dr > dq? dr - dq : dq - dr;//gap + if(dd <= long_gap) dd = 0; + sc = f[j] - dd; + if (sc > max_f) { + max_f = sc, max_j = j; + if (n_skip > 0) --n_skip; + } else if (t[j] == (int32_t)i) { + if (++n_skip > max_skip) + break; + } + if (p[j] >= 0) t[p[j]] = i; + } + + f[i] = max_f; p[i] = max_j; + } + + i = a_n-1; msc = f[i]; msc_i = i; + for (j = i-1; (j >= 0) && (a[i].self_offset <= (max_dis + a[j].self_offset)); --j) { + if(msc < f[j] && p[j] >= 0) {///hold at least two hits in th final chain + msc = f[j]; msc_i = j; + } + } + + ss_kip: + ///a[] has been sorted by self_offset + i = msc_i; + cL = 0; + while (i >= 0) { + t[cL++] = i; i = p[i]; + } + + n_skip = cL>>1; + for (i = 0; i < n_skip; i++) { + msc_i = t[i]; t[i] = t[cL-i-1]; t[cL-i-1] = msc_i; + } + if(des) { + for (i = 0; i < cL; i++) des[i] = a[t[i]]; + } + 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) { @@ -1680,10 +1960,13 @@ int64_t filter_bad_seed_dp(k_mer_hit *sk, k_mer_hit *ek, k_mer_hit* a, int64_t a uint64_t lchain_dp_trace(k_mer_hit* a, int64_t a_n, int64_t max_lgap, double sgap_rate, int64_t sgap) { + // fprintf(stderr, "[M::%s::] a_n::%ld\n", __func__, a_n); if(a_n <= 0) return 0; int64_t i, st, occ = 0; for (i = 1, st = 0; i <= a_n; ++i) { + // fprintf(stderr, "[M::%s::i->%ld] q::%u, t::%u, cnt::%u, readID::%u\n", __func__, + // i-1, a[i-1].self_offset, a[i-1].offset, a[i-1].cnt, a[i-1].readID); 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); - +uint64_t lchain_qdp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, overlap_region* res, + 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 xl, int64_t yl, int64_t quick_check); +int ovlp_chain_qgen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int64_t yl, int64_t apend_be, k_mer_hit* hit, int64_t n_hit); +uint64_t lchain_refine(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, + int64_t max_skip, int64_t max_iter, int64_t max_dis, int64_t long_gap); #endif diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index d286ef2..706d33d 100644 --- a/Levenshtein_distance.h +++ b/Levenshtein_distance.h @@ -542,8 +542,14 @@ inline uint32_t pop_trace(asg16_v *res, uint32_t i, uint16_t *c, uint32_t *len) // #define MAX_L 2500 ///511 -> 32 64-bits -#define MAX_E 1023 -#define MAX_L 5120 +#define MAX_CNS_E 1023 +#define MAX_CNS_L 3072 +#define FORCE_CNS_L 256 + + +#define MAX_SIN_E 2047 +#define MAX_SIN_L 10000 +#define FORCE_SIN_L 512 typedef uint64_t w_sig; #define bitw (6) @@ -1226,7 +1232,7 @@ inline void ed_band_cal_extension_##sf##_0_w(char *pstr, int32_t pn, char *tstr, {\ poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/\ if(k >= 0) {\ - if(k == 0) {\ + if(tmp_e == INT32_MAX) {\ tmp_e = err;\ for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {\ bd = (k>>bitw); k_bd = (k&bitz);\ @@ -1316,7 +1322,7 @@ inline void ed_band_cal_extension_##sf##_1_w(char *pstr, int32_t pn, char *tstr, {\ poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/\ if(k >= 0) {\ - if(k == 0) {\ + if(tmp_e == INT32_MAX) {\ tmp_e = err;\ for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {\ bd = (k>>bitw); k_bd = (k&bitz);\ @@ -1572,7 +1578,7 @@ inline void ed_band_cal_extension_##sf##_0_w_trace(char *pstr, int32_t pn, char if(!done) {\ poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/\ if(k >= 0) {\ - if(k == 0) {\ + if(tmp_e == INT32_MAX) {\ tmp_e = err;\ for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {\ bd = (k>>bitw); k_bd = (k&bitz);\ @@ -1690,7 +1696,7 @@ inline void ed_band_cal_extension_##sf##_1_w_trace(char *pstr, int32_t pn, char if(!done) {\ poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/\ if(k >= 0) {\ - if(k == 0) {\ + if(tmp_e == INT32_MAX) {\ tmp_e = err;\ for ((k) = 0; (poff) < (pe); (poff)++, (k)++) {\ bd = (k>>bitw); k_bd = (k&bitz);\ @@ -2157,7 +2163,7 @@ inline void ed_band_cal_extension_infi_0_w(char *pstr, int32_t pn, char *tstr, i { poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ if(k >= 0) { - if(k == 0) { + if(tmp_e == INT32_MAX) { tmp_e = err; for ((k) = 0; (poff) < (pe); (poff)++, (k)++) { bd = (k>>bitw); k_bd = (k&bitz); @@ -2251,7 +2257,7 @@ inline void ed_band_cal_extension_infi_1_w(char *pstr, int32_t pn, char *tstr, i { poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ if(k >= 0) { - if(k == 0) { + if(tmp_e == INT32_MAX) { tmp_e = err; for ((k) = 0; (poff) < (pe); (poff)++, (k)++) { bd = (k>>bitw); k_bd = (k&bitz); @@ -2540,7 +2546,7 @@ inline void ed_band_cal_extension_infi_0_w_trace(char *pstr, int32_t pn, char *t if(!done) {//diff poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ if(k >= 0) { - if(k == 0) { + if(tmp_e == INT32_MAX) { tmp_e = err; for ((k) = 0; (poff) < (pe); (poff)++, (k)++) { bd = (k>>bitw); k_bd = (k&bitz); @@ -2664,7 +2670,7 @@ inline void ed_band_cal_extension_infi_1_w_trace(char *pstr, int32_t pn, char *t if(!done) {//diff poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ if(k >= 0) { - if(k == 0) { + if(tmp_e == INT32_MAX) { tmp_e = err; for ((k) = 0; (poff) < (pe); (poff)++, (k)++) { bd = (k>>bitw); k_bd = (k&bitz); @@ -3038,7 +3044,7 @@ inline void ed_band_cal_extension_64_0_w(char *pstr, int32_t pn, char *tstr, int { poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ if(k >= 0) { - if(k == 0) { + if(tmp_e == INT32_MAX) { tmp_e = err; for (k = 0; poff < pe; poff++, (k)++) { tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL)); @@ -3111,7 +3117,7 @@ inline void ed_band_cal_extension_64_1_w(char *pstr, int32_t pn, char *tstr, int { poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ if(k >= 0) { - if(k == 0) { + if(tmp_e == INT32_MAX) { tmp_e = err; for (k = 0; poff < pe; poff++, (k)++) { tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL)); @@ -3337,7 +3343,7 @@ inline void ed_band_cal_extension_64_0_w_trace(char *pstr, int32_t pn, char *tst if(!done) {//diff poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ if(k >= 0) { - if(k == 0) { + if(tmp_e == INT32_MAX) { tmp_e = err; for (k = 0; poff < pe; poff++, (k)++) { tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL)); @@ -3351,9 +3357,10 @@ inline void ed_band_cal_extension_64_0_w_trace(char *pstr, int32_t pn, char *tst if(tmp_e <= (*ez).thre && tmp_e < (*ez).err) { (*ez).err = tmp_e; (*ez).pe = pe; (*ez).te = i; } + // fprintf(stderr, "i::%d, (*ez).err::%d, (*ez).pe::%d, (*ez).te::%d, thre::%d, err::%d, tmp_e::%d, poff::%d, k::%d, pe::%d\n", + // i, (*ez).err, (*ez).pe, (*ez).te, thre, err, tmp_e, i-thre, i+thre-pe, pe); } } - // fprintf(stderr, "i::%d, tmp_e::%d, err::%d\n", i, tmp_e, err); Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1; ++i; ++i_bd; @@ -3446,7 +3453,7 @@ inline void ed_band_cal_extension_64_1_w_trace(char *pstr, int32_t pn, char *tst if(!done) {//diff poff = i-thre; k = i+thre-pe;/**poff:[i-thre, i+thre]**/ if(k >= 0) { - if(k == 0) { + if(tmp_e == INT32_MAX) { tmp_e = err; for (k = 0; poff < pe; poff++, (k)++) { tmp_e += ((VP>>k)&(1ULL)); tmp_e -= ((VN>>k)&(1ULL)); diff --git a/anchor.cpp b/anchor.cpp index 502a1f3..79ee111 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -17,8 +17,11 @@ typedef struct { // this struct is not strictly necessary; we can use k_mer_pos #define an_key1(a) ((a).srt) #define an_key2(a) ((a).self_off) +#define an_key3(a) ((a).other_off) KRADIX_SORT_INIT(ha_an1, anchor1_t, an_key1, 8) KRADIX_SORT_INIT(ha_an2, anchor1_t, an_key2, 4) +KRADIX_SORT_INIT(ha_an3, anchor1_t, an_key3, 4) + #define oreg_xs_lt(a, b) (((uint64_t)(a).x_pos_s<<32|(a).x_pos_e) < ((uint64_t)(b).x_pos_s<<32|(b).x_pos_e)) KSORT_INIT(or_xs, overlap_region, oreg_xs_lt) @@ -872,6 +875,102 @@ void *ha_flt_tab, ha_pt_t *ha_idx, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_ cl->length = ab->n_a; } + +void minimizers_qgen(ha_abufl_t *ab, char* rs, int64_t rl, uint64_t mz_w, uint64_t mz_k, Candidates_list *cl, kvec_t_u8_warp* k_flag, +void *ha_flt_tab, ha_pt_t *ha_idx, All_reads* rdb, const ul_idx_t *udb, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, +uint32_t *low_occ) +{ + // fprintf(stderr, "+[M::%s]\n", __func__); + uint64_t i, k, l, max_cnt = UINT32_MAX, min_cnt = 0; int n, j; ha_mzl_t *z; seedl_t *s; + if(high_occ) { + max_cnt = (*high_occ); + if(max_cnt < 2) max_cnt = 2; + } + if(low_occ) { + min_cnt = (*low_occ); + if(min_cnt < 2) min_cnt = 2; + } + clear_Candidates_list(cl); ab->mz.n = 0, ab->n_a = 0; + + // get the list of anchors + mz2_ha_sketch(rs, rl, mz_w, mz_k, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, NULL); + + // minimizer of queried read + if (ab->mz.m > ab->old_mz_m) { + ab->old_mz_m = ab->mz.m; + REALLOC(ab->seed, ab->old_mz_m); + } + + for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) { + ab->seed[i].a = ha_ptl_get(ha_idx, ab->mz.a[i].x, &n); + ab->seed[i].n = n; + ab->n_a += n; + } + + if (ab->n_a > ab->m_a) { + ab->m_a = ab->n_a; + REALLOC(ab->a, ab->m_a); + } + + for (i = 0, k = 0; i < ab->mz.n; ++i) { + ///z is one of the minimizer + z = &ab->mz.a[i]; s = &ab->seed[i]; + for (j = 0; j < s->n; ++j) { + const ha_idxposl_t *y = &s->a[j]; + anchor1_t *an = &ab->a[k++]; + uint8_t rev = z->rev == y->rev? 0 : 1; + an->other_off = rev?((uint32_t)-1)-1-(y->pos+1-y->span):y->pos; + an->self_off = z->pos; + ///an->cnt: cnt<<8|span + an->cnt = s->n; if(an->cnt > ((uint32_t)(0xffffffu))) an->cnt = 0xffffffu; + an->cnt <<= 8; an->cnt |= ((z->span <= ((uint32_t)(0xffu)))?z->span:((uint32_t)(0xffu))); + an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->self_off; + } + } + + // copy over to _cl_ + if (ab->m_a >= (uint64_t)cl->size) { + cl->size = ab->m_a; + REALLOC(cl->list, cl->size); + } + + k_mer_hit *p; uint64_t tid = (uint64_t)-1, tl = (uint64_t)-1; + radix_sort_ha_an1(ab->a, ab->a + ab->n_a); + for (k = 1, l = 0; k <= ab->n_a; ++k) { + if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) { + if (k-l>1) radix_sort_ha_an3(ab->a+l, ab->a+k); + if((ab->a[l].srt>>33)!=tid) { + tid = ab->a[l].srt>>33; + tl = rdb?Get_READ_LENGTH((*rdb), tid):udb->ug->u.a[tid].len; + } + for (i = l; i < k; i++) { + p = &cl->list[i]; + p->readID = ab->a[i].srt>>33; + p->strand = (ab->a[i].srt>>32)&1; + if(!(p->strand)) { + p->offset = ab->a[i].other_off; + } else { + p->offset = ((uint32_t)-1)-ab->a[i].other_off; + p->offset = tl-p->offset; + } + p->self_offset = ab->a[i].self_off; + if(((ab->a[i].cnt>>8) < max_cnt) && ((ab->a[i].cnt>>8) > min_cnt)){ + p->cnt = 1; + } else if((ab->a[i].cnt>>8) <= min_cnt) { + p->cnt = 2; + } else{ + p->cnt = 1 + (((ab->a[i].cnt>>8) + (max_cnt<<1) - 1)/(max_cnt<<1)); + p->cnt = pow(p->cnt, 1.1); + } + if(p->cnt > ((uint32_t)(0xffffffu))) p->cnt = 0xffffffu; + p->cnt <<= 8; p->cnt |= (((uint32_t)(0xffu))&(ab->a[i].cnt)); + } + l = k; + } + } + cl->length = ab->n_a; +} + void inline reverse_k_mer_hit(k_mer_hit *a, uint64_t a_n, uint64_t xl, uint64_t yl) { uint64_t z, han = a_n>>1; k_mer_hit *ai, *aj, ka; @@ -1009,6 +1108,81 @@ void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uin ks_introsort_or_xs(ol->length, ol->list); } +void lchain_qgen(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) +{ + // fprintf(stderr, "+[M::%s]\n", __func__); + uint64_t i, k, l, m, sm, cn = cl->length; overlap_region *r; ///srt = 0 + clear_overlap_region_alloc(ol); + clear_fake_cigar(&(tf->f_cigar)); + + for (l = 0, k = 1, m = 0; k <= cn; k++) { + if((k == cn) || (cl->list[k].readID != cl->list[l].readID) + || (cl->list[k].strand != cl->list[l].strand)) { + if(cl->list[l].readID != rid) { + tf->x_id = rid; + tf->x_pos_strand = cl->list[l].strand; + tf->y_id = cl->list[l].readID; + tf->y_pos_strand = 0;///always 0 + // fprintf(stderr, "+[M::%s] l::%lu, k::%lu\n", __func__, l, k); + sm = lchain_qdp(cl->list+l, k-l, cl->list+m, &(cl->chainDP), tf, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate, + 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_qgen(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(tf->y_id == 66 || tf->y_id == 66) { + // fprintf(stderr, "\n[M::%s::] utg%.6dl(%c), i::%lu\n", + // __func__, (int32_t)tf->y_id+1, "+-"[tf->x_pos_strand], m); + // } + // 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)); + for (i = 0; i < sm; i++) { + cl->list[m+i].readID = ol->length; + // if(tf->y_id == 126) fprintf(stderr, "[M::%s::qoff->%u::toff->%u]\n", __func__, cl->list[m+i].self_offset, cl->list[m+i].offset); + } + if(gen_off) gen_fake_cigar(&(r->f_cigar), r, apend_be, cl->list+m, sm); + m += sm; + } + } + l = k; + } + } + cl->length = m; + + + k = ol->length; + if (ol->length > max_n_chain) { + int32_t w, n[4], s[4]; overlap_region t; + n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0; + ks_introsort_or_ss(ol->length, ol->list); ///srt = 1; + for (i = 0; i < ol->length; ++i) { + r = &(ol->list[i]); + w = ha_ov_type(r, rl); + ++n[w]; + if (((uint64_t)n[w]) == max_n_chain) s[w] = r->shared_seed; + } + if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) { + // n[0] = n[1] = n[2] = n[3] = 0; + for (i = 0, k = 0; i < ol->length; ++i) { + r = &(ol->list[i]); + w = ha_ov_type(r, rl); + // ++n[w]; + // if (((int)n[w] <= max_n_chain) || (r->shared_seed >= s[w] && s[w] >= (asm_opt.k_mer_length<<1))) { + if (r->shared_seed >= s[w]) { + if (k != i) { + t = ol->list[k]; + ol->list[k] = ol->list[i]; + ol->list[i] = t; + } + ++k; + } + } + ol->length = k; + } + } + ks_introsort_or_xs(ol->length, ol->list); +} + void set_lchain_dp_op(uint32_t is_accurate, uint32_t mz_k, int64_t *max_skip, int64_t *max_iter, int64_t *max_dis, double *chn_pen_gap, double *chn_pen_skip, int64_t *quick_check) { double div, pen_gap, pen_skip, tmp; @@ -1031,7 +1205,9 @@ void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t extern ha_pt_t *ha_idx; int64_t max_skip, max_iter, max_dis, quick_check; double chn_pen_gap, chn_pen_skip; set_lchain_dp_op(is_accurate, mz_k, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check); - minimizers_gen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, dbg_ct, sp, high_occ, low_occ); - lchain_gen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off); + // minimizers_gen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, dbg_ct, sp, high_occ, low_occ); + minimizers_qgen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, NULL, uref, dbg_ct, sp, high_occ, low_occ); + // lchain_gen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off); + lchain_qgen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off); ///no need to sort here, overlap_list has been sorted at lchain_gen } diff --git a/inter.cpp b/inter.cpp index 1830dc6..13fa368 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!=1) /**&& (s->id+i!=44) && (s->id+i!=948)**/) return; + // if((s->id+i!=291) /**&& (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, s->opt->k, NULL); + &b->correct, &b->exz, &b->exz1, &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 @@ -6092,7 +6092,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call 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, s->opt->k, NULL); + &b->correct, &b->exz, &b->exz1, &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); } diff --git a/main.cpp b/main.cpp index 19fa21d..d58a8e1 100644 --- a/main.cpp +++ b/main.cpp @@ -12,16 +12,18 @@ int main(int argc, char *argv[]) yak_reset_realtime(); init_opt(&asm_opt); if (!CommandLine_process(argc, argv, &asm_opt)) return 0; - /** - bit_extz_t exz, exz64; init_bit_extz_t(&exz, 2); init_bit_extz_t(&exz64, 2); - char *pstr = "AAATATT", *tsrt = "AATATTTT"; int32_t thre = 1; - ed_band_cal_semi_64_w_absent_diag((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, 0, &exz); - fprintf(stderr, "\n[M::%s::] exz.err::%d, exz.ps::%d, exz.pe::%d, exz.ts::%d, exz.te::%d\n", __func__, - exz.err, exz.ps, exz.pe, exz.ts, exz.te); - ed_band_cal_extension_64_0_w((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, &exz); - fprintf(stderr, "\n[M::%s::] exz.err::%d, exz.ps::%d, exz.pe::%d, exz.ts::%d, exz.te::%d\n", __func__, - exz.err, exz.ps, exz.pe, exz.ts, exz.te); + // bit_extz_t exz, exz64; init_bit_extz_t(&exz, 2); init_bit_extz_t(&exz64, 2); + + // char *pstr = "GACCCAG", *tsrt = "GTTGTTAATTCCAT"; int32_t thre = 14; clear_align(exz); clear_align(exz64); + // ed_band_cal_extension_64_0_w_trace((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, &exz); + // // // ed_band_cal_semi_64_w_absent_diag((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, 0, &exz); + // fprintf(stderr, "\n[M::%s::] exz.err::%d, exz.ps::%d, exz.pe::%d, exz.ts::%d, exz.te::%d\n", __func__, + // exz.err, exz.ps, exz.pe, exz.ts, exz.te); + // cigar_check((char*)pstr, (char*)tsrt, &(exz)); + // ed_band_cal_extension_64_0_w((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, &exz); + // fprintf(stderr, "\n[M::%s::] exz.err::%d, exz.ps::%d, exz.pe::%d, exz.ts::%d, exz.te::%d\n", __func__, + // exz.err, exz.ps, exz.pe, exz.ts, exz.te); // ed_band_cal_semi_infi_w((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, NULL, &exz); // ed_band_cal_semi_64_w((char*)pstr, strlen(pstr), (char*)tsrt, strlen(tsrt), thre, &exz64); // fprintf(stderr, "\n[M::%s::] exz.err::%d, exz64.err::%d, exz.ps::%d, exz64.ps::%d, exz.pe::%d, exz64.pe::%d, exz.ts::%d, exz64.ts::%d, exz.te::%d, exz64.te::%d\n", __func__, @@ -47,8 +49,7 @@ int main(int argc, char *argv[]) // ed_band_cal_extension_128bit((char *)"AAGTTTA", 7, (char *)"CCTTTTTT", 8, 4, &exz); // ed_band_cal_extension_128bit((char *)"AA", 2, (char *)"ACTTTTTT", 8, 1, &exz); // fprintf(stderr, "ed_extension::%d, pe::%d, te::%d\n", exz.err, exz.pe, exz.te); - exit(1); - **/ + // exit(1);