From b3b18ab1d08e17b954a2dafcb4234a8b11da918c Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Fri, 14 Mar 2025 01:40:21 -0400 Subject: [PATCH] r721 --- CommandLines.h | 2 +- Correct.cpp | 120 +++++++++++++++++++++++++++++++++------ Correct.h | 4 +- ecovlp.cpp | 149 ++++++++++++++++++++++++++++++++++++++++++++++--- 4 files changed, 246 insertions(+), 29 deletions(-) diff --git a/CommandLines.h b/CommandLines.h index 9725a49..511a5f5 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.25.0-r710" +#define HA_VERSION "0.25.0-r721" #define VERBOSE 0 diff --git a/Correct.cpp b/Correct.cpp index 4162d98..809e71a 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -17423,11 +17423,11 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u void hc_ovlp_base_direct(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, int64_t wl, All_reads *rref, char* qstr, UC_Read *tu, -bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, uint64_t rid) +bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, uint64_t rid, int64_t pre_mode) { - int64_t i, l, mode, q[2], t[2], qr, tr, is_done, zn; + int64_t i, l, mode, q[2], t[2], qr, tr, is_done, zn, si, ei; - if(z->non_homopolymer_errors == 0 && z->w_list.n) { + if((pre_mode < 0) && (z->non_homopolymer_errors == 0) && (z->w_list.n)) { zn = z->w_list.n; for (i = 1; i < zn; i++) { if((z->w_list.a[i].error == 0 && z->w_list.a[i-1].error == 0) && (z->w_list.a[i].x_start == z->w_list.a[i-1].x_end + 1) && @@ -17461,7 +17461,16 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u } } - for (l = -1, i = 0; i <= ch_n; i++) { + si = 0; ei = ch_n; + if(pre_mode == 0) { + si = 1; ei = ch_n - 1; + } else if(pre_mode == 1) { + si = 1; + } else if(pre_mode == 2) { + ei = ch_n - 1; + } + + for (l = si - 1, i = si; i <= ei; 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; @@ -17702,7 +17711,7 @@ All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate, int64_ // idx.qs, idx.qe, idx.ts, idx.te, idx.qn, idx.tn); // } // ovlp_base_aln(z, ch_a, ch_n, &idx, wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, (uint64_t)-1); - hc_ovlp_base_direct(z, ch_a, ch_n, wl, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, (uint64_t)-1); + hc_ovlp_base_direct(z, ch_a, ch_n, wl, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, (uint64_t)-1, mode); an = aux_o->w_list.n; q[0] = q[1] = t[0] = t[1] = 0; todo = 0; // if(z->x_id == 29033 && z->y_id == 21307) { // fprintf(stderr, "[M::%s]\tan::%ld\n", __func__, an); @@ -17820,11 +17829,17 @@ uint64_t gen_hc_fast_cigar0(overlap_region *z, Candidates_list *cl, uint64_t wl, aux_o->x_pos_s = z->x_pos_s; aux_o->x_pos_e = z->x_pos_e; aux_o->y_pos_s = z->y_pos_s; aux_o->y_pos_e = z->y_pos_e; - hc_ovlp_base_direct(z, ch_a, ch_n, wl, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, rid); + hc_ovlp_base_direct(z, ch_a, ch_n, wl, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, rid, -1); - int64_t aux_n = aux_o->w_list.n; for (i = 0; i < aux_n; i++) { + // if(z->y_id == 30129) { + // fprintf(stderr, "[aln::-i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, clen::%u, mode::%d\n", i, + // aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start, + // aux_o->w_list.a[i].x_start, aux_o->w_list.a[i].x_end+1, + // aux_o->w_list.a[i].y_start, aux_o->w_list.a[i].y_end+1, + // aux_o->w_list.a[i].error, aux_o->w_list.a[i].clen, aux_o->w_list.a[i].error_threshold); + // } if(!(is_ualn_win(aux_o->w_list.a[i]))) continue; // if((aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start) <= FORCE_CNS_L) { // fprintf(stderr, "[aln::-i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, clen::%u, mode::%d\n", i, @@ -19761,7 +19776,7 @@ inline uint64_t set_cgid(ul_ov_t *z, uint64_t *ia, uint64_t *iak, uint64_t ian, } ///need to print some examples for double check -uint64_t is_get_group(ul_ov_t *a, uint64_t *ga, uint64_t *ia, uint64_t gi, uint64_t tn) +uint64_t is_get_group(ul_ov_t *a/**, uint64_t an, uint64_t rid**/, uint64_t *ga, uint64_t *ia, uint64_t gi, uint64_t tn) { // fprintf(stderr, "+[M::%s] gi::%lu\n", __func__, gi); uint64_t k = ga[gi]>>32; @@ -19771,7 +19786,10 @@ uint64_t is_get_group(ul_ov_t *a, uint64_t *ga, uint64_t *ia, uint64_t gi, uint6 // fprintf(stderr, "-[M::%s] k::%lu\n", __func__, k); // } for (k = ga[gi]>>32; (k != ((uint64_t)-1)) && (a[k].tn != tn); k = ia[k]); - if((k != ((uint32_t)-1)) && (a[k].tn == tn)) return 0; + // if((k != ((uint32_t)-1)) && (k >= an)) { + // fprintf(stderr, "-[M::%s] rid::%lu, an::%lu, k::%lu\n", __func__, rid, an, k); + // } + if((k != ((uint64_t)-1)) && (a[k].tn == tn)) return 0; return 1; } @@ -19867,12 +19885,13 @@ int64_t rphase_lidel_cc(overlap_region_alloc* oa, ul_ov_t *a, int64_t an, double ol = (a[k].qe - a[k].qs) * len_st; if(ol < len_w) ol = len_w; if((a[k].qe - a[k].qs) < ol) continue; mk = -1; msw = -1; - for (z = k + 1; (z < an) && (a[z].qs < a[k].qe); z++) { + for (z = 0/**k + 1**/; (z < an) && (a[z].qs < a[k].qe); z++) { if(a[z].sec < c_sz) continue; if(a[z].ts == ((uint32_t)-1)) continue; + if(z == k) continue; sw = cal_lindel_dd(&(a[k]), &(a[z]), len_st, len_w, err_dif, c_sz); if(sw < 0) continue; - if((sw > msw) && (is_get_group(a, ca, ia, a[z].ts, a[k].tn))) { + if((sw > msw) && (is_get_group(a, /**an, rid,**/ ca, ia, a[z].ts, a[k].tn))) { mk = a[z].ts; msw = sw; } } @@ -22979,7 +22998,7 @@ int64_t return_t_chain(overlap_region *z, Candidates_list *cl) { int64_t i, cn = cl->length, scn; uint64_t pid; k_mer_hit *ca; - // if(z->y_id == 4378833 || z->y_id == 4378837) { + // if(z->y_id == 30129) { // fprintf(stderr, "\n-0-[M::%s]\tutg%.6ul\tx::[%u,\t%u)\t%c\tutg%.6ul\ty::[%u,\t%u)\n", // __func__, z->x_id+1, z->x_pos_s, z->x_pos_e+1, // "+-"[z->y_pos_strand], z->y_id+1, z->y_pos_s, z->y_pos_e+1); @@ -22998,14 +23017,14 @@ int64_t return_t_chain(overlap_region *z, Candidates_list *cl) scn = i - z->shared_seed; ca = cl->list+z->shared_seed; i = lchain_refine(ca, scn, ca, &(cl->chainDP), 50, 5000, 512, 16); cn = i; for (; i < scn; i++) ca[i].readID = ((uint32_t)(0x7fffffff)); + // if(z->y_id == 30129) fprintf(stderr, "\n-a-[M::%s]\tcn::%ld\n", __func__, cn); - - // if(z->y_id == 4378833 || z->y_id == 4378837) { + // if(z->y_id == 30129) { // fprintf(stderr, "\n-1-[M::%s]\tutg%.6ul\tx::[%u,\t%u)\t%c\tutg%.6ul\ty::[%u,\t%u)\n", // __func__, z->x_id+1, z->x_pos_s, z->x_pos_e+1, // "+-"[z->y_pos_strand], z->y_id+1, z->y_pos_s, z->y_pos_e+1); - // i = z->shared_seed; pid = cl->list[i].readID; + // i = z->shared_seed; pid = cl->list[i].readID; cn = cl->length; // for (; i < cn && cl->list[i].readID == pid && cl->list[i].readID != ((uint32_t)(0x7fffffff)); i++) { // fprintf(stderr, "i::%ld[M::%s]\treadID::%u\tself_offset::%u\toffset::%u\t%c\n", // i, __func__, cl->list[i].readID, cl->list[i].self_offset, cl->list[i].offset, @@ -25532,7 +25551,70 @@ uint32_t inline ff_tend(overlap_region *z, int64_t wn, int64_t dn, double dr, do return 0; } -void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop) +uint32_t inline ff_lunalign(overlap_region *z, double erate, double gap_rate, int64_t max_gap) +{ + if((z->w_list.n == 1) && (!(is_ualn_win(z->w_list.a[0]))) + && (z->w_list.a[0].x_start == ((int64_t)z->x_pos_s)) && (z->w_list.a[0].x_end == ((int64_t)z->x_pos_e)) + && (z->w_list.a[0].y_start == ((int64_t)z->y_pos_s)) && (z->w_list.a[0].y_end == ((int64_t)z->y_pos_e))) { + return 1; + } + + int64_t k, zwn = z->w_list.n, zq, zt, wq, wt, tot_e, tot_g, ql, tl; + + // fprintf(stderr, "[M::%s]\n", __func__); + + zq = z->x_pos_s; zt = z->y_pos_s; tot_e = tot_g = 0; + for (k = 0; k < zwn; k++) { + // fprintf(stderr, "[M::%s]\twk::%ld\tq::[%d, %d)\tt::[%d, %d)\terr::%d\n", __func__, + // k, z->w_list.a[k].x_start, z->w_list.a[k].x_end + 1, z->w_list.a[k].y_start, z->w_list.a[k].y_end + 1, z->w_list.a[k].error); + if(is_ualn_win(z->w_list.a[k])) continue; + wq = z->w_list.a[k].x_start; + wt = z->w_list.a[k].y_start; + + if(wq != zq) tot_g += ((wq>=zq)?(wq-zq):(zq-wq)); + if(wt != zt) tot_g += ((wt>=zt)?(wt-zt):(zt-wt)); + + zq = z->w_list.a[k].x_end + 1; + zt = z->w_list.a[k].y_end + 1; + tot_e += z->w_list.a[k].error; + + // fprintf(stderr, "[M::%s]\twk::%ld\tq::[%d, %d)\tt::[%d, %d)\terr::%d\n", __func__, + // k, z->w_list.a[k].x_start, z->w_list.a[k].x_end + 1, z->w_list.a[k].y_start, z->w_list.a[k].y_end + 1, z->w_list.a[k].error); + } + + wq = z->x_pos_e + 1; + wt = z->y_pos_e + 1; + if(wq != zq) tot_g += ((wq>=zq)?(wq-zq):(zq-wq)); + if(wt != zt) tot_g += ((wt>=zt)?(wt-zt):(zt-wt)); + + // fprintf(stderr, "[M::%s]\t%.*s(id::%u)\tq::[%u, %u)\t%.*s(id::%u)\tt::[%u, %u)\tre::%ld\trg::%ld\n", __func__, (int)Get_NAME_LENGTH(R_INF, z->x_id), Get_NAME(R_INF, z->x_id), z->x_id, z->x_pos_s, z->x_pos_e + 1, + // (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), z->y_id, z->y_pos_s, z->y_pos_e + 1, tot_e, tot_g); + + if(!tot_g) return 1; + + // fprintf(stderr, "-0-[M::%s]\n", __func__); + + if(tot_g > max_gap) return 0; + + // fprintf(stderr, "-1-[M::%s]\n", __func__); + + ql = z->x_pos_e + 1 - z->x_pos_s; + tl = z->y_pos_e + 1 - z->y_pos_s; + + if((tot_g > (ql*gap_rate)) || (tot_g > (tl*gap_rate))) return 0; + + // fprintf(stderr, "-2-[M::%s]\n", __func__); + + tot_e += tot_g; + + if((tot_e > (ql*erate)) || (tot_e > (tl*erate))) return 0; + + // fprintf(stderr, "-3-[M::%s]\n", __func__); + + return 1; +} + +void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max) { uint64_t i, bs, k, ql = qu->length; Window_Pool w; double err, e_max, rr; int64_t re; overlap_region t; overlap_region *z; //asg64_v iidx, buf, buf1; @@ -25566,6 +25648,8 @@ void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rre if(!gen_hc_fast_cigar(z, cl, rref, w.window_length, qu->seq, tu, exz, aux_o, e_rate, ql, rid, khit, &re)) continue; + if((align_gap_max >= 0) && (!ff_lunalign(z, err, align_gap_rate, align_gap_max))) continue; + if(chem_drop && ff_tend(z, 384, 2000, 0.1, (((e_rate*10)<0.36)?(e_rate*10):(0.36)), 128)) continue; @@ -25590,7 +25674,7 @@ void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rre } -void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop) +void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max) { uint64_t i, bs, k, ql = qu->length; Window_Pool w; double err, e_max, rr; int64_t re; overlap_region t; overlap_region *z; //asg64_v iidx, buf, buf1; @@ -25632,6 +25716,8 @@ void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads if(!gen_hc_fast_cigar(z, cl, rref, w.window_length, qu->seq, tu, exz, aux_o, e_rate, ql, rid, khit, &re)) continue; + if((align_gap_max >= 0) && (!ff_lunalign(z, err, align_gap_rate, align_gap_max))) continue; + if(chem_drop && ff_tend(z, 384, 2000, 0.1, (((e_rate*10)<0.36)?(e_rate*10):(0.36)), 128)) continue; // if(z->x_id == 3196 && z->y_id == 3199) fprintf(stderr, "-1-[M::%s] tid::%u\t%.*s\trr::%f\tre::%ld\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), rr, re); diff --git a/Correct.h b/Correct.h index bf1d6cb..2839e60 100644 --- a/Correct.h +++ b/Correct.h @@ -1391,8 +1391,8 @@ int64_t get_rid_backward_cigar_err(rtrace_iter *it, ul_ov_t *aln, kv_rtrace_t *t const ul_idx_t *uref, char* qstr, UC_Read *tu, overlap_region_alloc *ol, overlap_region *o, bit_extz_t *exz, double e_rate, int64_t qs); -void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop); -void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop); +void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max); +void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max); uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf); void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, Chain_Data *dp, asg8_v *q8, asg8_v *t8, uint8_t lindel); void set_exact_exz(bit_extz_t *exz, int64_t qs, int64_t qe, int64_t ts, int64_t te); diff --git a/ecovlp.cpp b/ecovlp.cpp index e3158fa..9b1fd45 100644 --- a/ecovlp.cpp +++ b/ecovlp.cpp @@ -65,6 +65,12 @@ typedef struct { uint8_t *cr; } ec_ovec_buf_t; +typedef struct { + ec_ovec_buf_t *p; + asg64_v idx; + ma_ug_t *ug; +} ec_polish_buf_t; + typedef struct { uint32_t n_thread, n_a, chunk_size, cn; FILE *fp; @@ -2801,7 +2807,7 @@ inline uint64_t exact_ec_check(char *qstr, uint64_t ql, char *tstr, uint64_t tl, return 0; } -void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in, uint8_t chem_drop) +void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max) { if(ol->length <= 0) return; @@ -2817,7 +2823,7 @@ void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads * } if(!(srt->n)) { - gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop); + gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max); } else { ///debug for memory // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); @@ -2852,9 +2858,10 @@ void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads * // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); if(on > nec) { - gen_hc_r_alin_nec(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop); + gen_hc_r_alin_nec(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max); } + // fprintf(stderr, "[M::%s] srt->n::%u, nec::%lu, on::%lu\n", __func__, (uint32_t)srt->n, nec, on); ///debug for memory // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); } @@ -3257,6 +3264,8 @@ static void worker_hap_ec(void *data, long i, int tid) ///for debug indel // if(i != 1238) return; // if(i != 2410) return; + // if(i != 4198005) return; + // if(i != 682) return; // debug_retrive_bqual(D, &b->v8t, i, 256); return; @@ -3276,7 +3285,7 @@ static void worker_hap_ec(void *data, long i, int tid) ///mz1_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); // if((asm_opt.is_ont) && (b->olist.length)) get_mz1(qu->seq, qu->length, RES_W, RES_K, 0, !(asm_opt.flag & HA_F_NO_HPC), b->ab, NULL, NULL, asm_opt.mz_sample_dist, NULL, NULL, NULL, -1, asm_opt.dp_min_len, -1, &(b->sp), asm_opt.mz_rewin, 0, NULL, 0); - gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT/**asm_opt.k_mer_length**/, 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont); + gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT/**asm_opt.k_mer_length**/, 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1)); ///for debug indel // prt_ovlp_sam(&b->olist, &b->ovlp_read, b->self_read.seq, b->self_read.length); @@ -3404,7 +3413,7 @@ static void worker_hap_ec_dbg_paf(void *data, long i, int tid) aux_o = fetch_aux_ovlp(&b->olist);///must be here // stderr_phase_ovlp(&b->olist); - gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, 1, &b->v16, &b->v64, &(R_INF.paf[i]), 0); + gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, 1, &b->v16, &b->v64, &(R_INF.paf[i]), 0, -1, -1); uint32_t k, m, tl; overlap_region *z; bit_extz_t ez; ma_hit_t *t; for (k = 0; k < b->olist.length; k++) { @@ -3796,6 +3805,42 @@ static void worker_hap_dc_ec(void *data, long i, int tid) refresh_ec_ovec_buf_t0(b, REFRESH_N); } +static void worker_update_dc_ec(void *data, long i, int tid) +{ + ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); + uint64_t k; ma_hit_t *z; + // fprintf(stderr, "-0-[M::%s-beg] rid->%ld\n", __func__, i); + // if (memcmp("m64012_190921_234837/139067658/ccs", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { + // fprintf(stderr, "-0-[M::%s-beg] rid->%ld\n", __func__, i); + // } else if (memcmp("m64012_190921_234837/28968323/ccs", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { + // fprintf(stderr, "-1-[M::%s-beg] rid->%ld\n", __func__, i); + // } else { + // return; + // } + // if(i != 2851) return; + + // if(scb.a[i].m < scc.a[i].n) { + // scb.a[i].m = scc.a[i].n; + // REALLOC(scb.a[i].a, scb.a[i].m); + // } + // scb.a[i].n = scc.a[i].n; + // memcpy(scb.a[i].a, scc.a[i].a, scc.a[i].n*sizeof((*(scb.a[i].a)))); + + + if(!(R_INF.paf[i].length)) return; + recover_UC_Read(&b->self_read, &R_INF, i); + for (k = 0; k < R_INF.paf[i].length; k++) { + z = &(R_INF.paf[i].buffer[k]); + if((z->el) && (quick_exact_match(z, &R_INF, &b->self_read, &b->ovlp_read, &scc))) { + z->el = 1; b->cnt[0]++; + } else { + z->el = 0; b->cnt[1]++; + } + } + + refresh_ec_ovec_buf_t0(b, REFRESH_N); +} + void flip_paf_rc(uint64_t rid, ma_hit_t_alloc *paf, All_reads *rref) { @@ -4280,6 +4325,52 @@ static void worker_hap_dc_ec_chemical_arc_mark(void *data, long i, int tid) refresh_ec_ovec_buf_t0(b, REFRESH_N); } +uint64_t get_candidate_rrs(ma_utg_t *u, uint64_t rz) +{ + // uint64_t rid, rs, re, rev, k, l[2], lr, ts, te; + // ma_hit_t_alloc *z = NULL; + // ma_hit_t *h = NULL; + + // rid = u->a[rz]>>33; rev = (u->a[rz]>>32)&1; + // rs = 0; re = (uint32_t)u->a[rz]; + // if(!rev) { + // rs = Get_READ_LENGTH(R_INF, rid) - ((uint32_t)u->a[rz]); + // re = Get_READ_LENGTH(R_INF, rid); + // } + + // for (k = lr = 0; k < rz; k++) lr += (uint32_t)u->a[k]; + // ts = lr; te = lr + (uint32_t)u->a[k]; + + // z = &(sources[rid]); + // for (k = 0; k < z->length; k++) { + // h = &(z->buffer[k]); + // if ((Get_qs(*h) <= rs) && (Get_qe(*h) >= re)) { + // l[0] = l[1] = 0; + // if(!(h->rev&rev)) { + // l[0] = Get_ts(*h); l[1] = Get_READ_LENGTH(R_INF, Get_tn(*h)) - Get_te(*h); + // } else { + // l[1] = Get_ts(*h); l[0] = Get_READ_LENGTH(R_INF, Get_tn(*h)) - Get_te(*h); + // } + + // } + // } + + + return 1; +} + +static void worker_ec_polish(void *data, long i, int tid) +{ + ec_ovec_buf_t0 *b = &(((ec_polish_buf_t*)data)->p->a[tid]); + uint64_t *idx = &(((ec_polish_buf_t*)data)->idx.a[i]); + + if(get_candidate_rrs(&(((ec_polish_buf_t*)data)->ug->u.a[(*idx)>>32]), (uint32_t)(*idx))) { + *idx = (uint64_t)-1; + } + + refresh_ec_ovec_buf_t0(b, REFRESH_N); +} + void gen_ovlst_paf(ma_hit_t_alloc *in_e, ma_hit_t_alloc *in_r, asg64_v *ou) { uint32_t n = 0, k; @@ -4482,7 +4573,7 @@ overlap_region* h_ec_lchain_re1(ha_abuf_t *ab, uint32_t rid, UC_Read *qu, UC_Rea // fprintf(stderr, "-0-[M::%s]\n", __func__); - gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf, 0); rs = qu->seq; + gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf, 0, -1, -1); rs = qu->seq; // fprintf(stderr, "-1-[M::%s]\n", __func__); @@ -5151,7 +5242,7 @@ overlap_region* h_ec_lchain_re3(ha_abuf_t *ab, uint32_t rid, UC_Read *qu, UC_Rea // fprintf(stderr, "-0-[M::%s]\n", __func__); - gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf, 0); rs = qu->seq; + gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf, 0, -1, -1); rs = qu->seq; // fprintf(stderr, "-1-[M::%s]\n", __func__); @@ -6001,6 +6092,22 @@ uint64_t cal_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, uint64 return num_correct; } +void cal_update_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a) +{ + double tt0 = yak_realtime_0(); + uint64_t k, num_ec_o = 0, num_nec_o = 0; + + for (k = 0; k < n_thre; ++k) b->a[k].cnt[0] = b->a[k].cnt[1] = 0; + + kt_for(n_thre, worker_update_dc_ec, b, n_a);///debug_for_fix + + for (k = 0; k < n_thre; ++k) { + num_ec_o += b->a[k].cnt[0]; num_nec_o += b->a[k].cnt[1]; + } + + fprintf(stderr, "[M::pec::%.3f] # exact o: %lu; # non-exact o: %lu\n", yak_realtime_0()-tt0, num_ec_o, num_nec_o); +} + void ha_print_ovlp_stat_1(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a) { @@ -6162,7 +6269,7 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u { // write_ec_reads("ec0.fa"); - fprintf(stderr, "[M::%s]\tn_thre::%lu, round::%lu, n_round::%lu, n_a::%lu, is_sv::%lu\n", __func__, n_thre, round, n_round, n_a, is_sv); + // fprintf(stderr, "[M::%s]\tn_thre::%lu, round::%lu, n_round::%lu, n_a::%lu, is_sv::%lu\n", __func__, n_thre, round, n_round, n_a, is_sv); ec_ovec_buf_t *b = NULL; uint64_t k, is_cr = (round&1); (*tot_b) = (*tot_e) = 0; @@ -6177,7 +6284,10 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u sl_ec_r(n_thre, n_a); } - if(is_sv) kt_for(n_thre, worker_hap_dc_ec, b, n_a);///update overlaps + cal_update_ec_multiple(b, n_thre, n_a);///update overlaps + + // if(is_sv) kt_for(n_thre, worker_hap_dc_ec, b, n_a);///update overlaps + if((!is_sv) || (is_sv && is_cr)) { kt_for(n_thre, worker_hap_post_rev, b, n_a); @@ -6379,4 +6489,25 @@ uint8_t* gen_chemical_arc_rf(uint64_t n_thre, uint64_t n_a) b->cr = NULL; destroy_ec_ovec_buf_t(b); return ra; +} + + +void gen_hc_polish(uint64_t n_thre, ma_ug_t *ug) +{ + ec_polish_buf_t b; memset(&b, 0, sizeof(b)); b.ug = ug; + uint64_t k, z, zn; + for (k = b.idx.n = 0; k < ug->u.n; k++) { + b.idx.n += ug->u.a[k].n; + } + MALLOC(b.idx.a, b.idx.n); + for (k = zn = 0; k < ug->u.n; k++) { + for (z = 0; z < ug->u.a[k].n; z++, zn++) { + b.idx.a[zn] = (k<<32)|z; + } + } + b.p = gen_ec_ovec_buf_t(n_thre); + + kt_for(n_thre, worker_ec_polish, &b, b.idx.n); + + destroy_ec_ovec_buf_t(b.p); free(b.idx.a); } \ No newline at end of file