diff --git a/CommandLines.h b/CommandLines.h index 84fe77f..a3a9893 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.25.0-r857" +#define HA_VERSION "0.25.0-r861" #define VERBOSE 0 diff --git a/Correct.cpp b/Correct.cpp index 59c210e..0ce20c7 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -33084,6 +33084,120 @@ void gen_hc_aln_small_chn_smp(overlap_region_alloc* ol, Candidates_list *cl, All } +void gen_hc_aln_small_chn_smp_ff_ec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, uint64_t ql, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, + uint64_t *wcut, uint64_t wcut_n, uint64_t ocw, uint64_t *srt_a, uint64_t i_srt_n, uint32_t *ocn, uint32_t *osc, uint64_t chain_cutoff, uint64_t ch_occ, uint64_t ch_sc) +{ + // if(!((i_srt_n > 0) && (i_srt_n <= ol->length))) { + // fprintf(stderr, "-1-[M::%s] (*i_srt_n)::%lu, ol->length::%lu\n", __func__, + // i_srt_n, ol->length); + // } + + assert((i_srt_n > 0) && (i_srt_n <= ol->length)); + // uint64_t ch_occ = 4, ch_sc = 16; + uint64_t k, m, s, e, ak, *srt_b = srt_a + i_srt_n, srt_ba = 0, bk, ncut = chain_cutoff*ch_occ; + for (k = s = bk = ak = 0; k < i_srt_n; k++) { + e = (uint32_t)srt_a[k]; + assert(s <= e); + for (m = s; m < e; m++){ + // if(!(ol->list[m].non_homopolymer_errors != UINT32_MAX)) { + // fprintf(stderr, "[M::%s::oid->%lu]\tqid::%u\ttid::%u\t%.*s\tq::[%u,%u)\tt::[%u,%u)\tsc::%u\tcn::%u\tis_match::%u\terr::%u\n", __func__, m, ol->list[m].x_id, + // ol->list[m].y_id, (int)Get_NAME_LENGTH(R_INF, ol->list[m].y_id), Get_NAME(R_INF, ol->list[m].y_id), + // ol->list[m].x_pos_s, ol->list[m].x_pos_e + 1, ol->list[m].y_pos_s, ol->list[m].y_pos_e + 1, + // osc[m], ocn[m], ol->list[m].is_match, ol->list[m].non_homopolymer_errors); + // } + assert(ol->list[m].non_homopolymer_errors != UINT32_MAX); + if(ocn[m] < ncut) continue; + srt_b[bk] = ol->list[m].x_pos_s; + srt_b[bk] <<= 32; srt_b[bk] |= m; bk++; + } + srt_ba += e - s; + s = e + 1; + if((srt_a[k]>>32) == ((uint32_t)-1)) { + m = e; + srt_a[ak] = ol->list[m].x_pos_s; + srt_a[ak] <<= 32; srt_a[ak] |= m; ak++; + } + } + + e = ol->length; + assert(s <= e); + for (m = s; m < e; m++){ + // if(!(ol->list[m].non_homopolymer_errors != UINT32_MAX)) { + // fprintf(stderr, "[M::%s::oid->%lu]\tqid::%u\ttid::%u\t%.*s\tq::[%u,%u)\tt::[%u,%u)\tsc::%u\tcn::%u\tis_match::%u\terr::%u\n", __func__, m, ol->list[m].x_id, + // ol->list[m].y_id, (int)Get_NAME_LENGTH(R_INF, ol->list[m].y_id), Get_NAME(R_INF, ol->list[m].y_id), + // ol->list[m].x_pos_s, ol->list[m].x_pos_e + 1, ol->list[m].y_pos_s, ol->list[m].y_pos_e + 1, + // osc[m], ocn[m], ol->list[m].is_match, ol->list[m].non_homopolymer_errors); + // } + assert(ol->list[m].non_homopolymer_errors != UINT32_MAX); + if(ocn[m] < ncut) continue; + srt_b[bk] = ol->list[m].x_pos_s; + srt_b[bk] <<= 32; srt_b[bk] |= m; bk++; + } + srt_ba += e - s; + + assert(srt_ba + i_srt_n == ol->length); + assert(ak > 0); + + radix_sort_bc64(srt_a, srt_a+ak); + radix_sort_bc64(srt_b, srt_b+bk); + + // fprintf(stderr, "-1-[M::%s] srt_ba::%lu, (*i_srt_n)::%lu, ol->length::%lu, ak::%lu, bk::%lu, ncut::%lu\n", __func__, + // srt_ba, i_srt_n, ol->length, ak, bk, ncut); + + uint64_t zs, ze, ob, zsc, zcn, rs, re, os, oe, oi, rr, kn, cs, ce, cn = cl->length; uint8_t f; overlap_region *zm, *rm; + for (m = 0; m < ak; m++) { + zm = &(ol->list[(uint32_t)srt_a[m]]); + zs = zm->x_pos_s; ze = zm->x_pos_e + 1; + ob = (ze - zs)*0.95; if(ob < 16) ob = 16; + zsc = osc[(uint32_t)srt_a[m]]*ch_sc; + zcn = ocn[(uint32_t)srt_a[m]]*ch_occ; + + // fprintf(stderr, "\n[M::%s::zm::oid->%u]\ttid::%u\t%.*s\tq::[%u,%u)\tt::[%u,%u)\tsc::%u\tcn::%u\tis_match::%u\terr::%u\n", __func__, (uint32_t)srt_a[m], + // zm->y_id, (int)Get_NAME_LENGTH(R_INF, zm->y_id), Get_NAME(R_INF, zm->y_id), + // zm->x_pos_s, zm->x_pos_e + 1, zm->y_pos_s, zm->y_pos_e + 1, + // osc[(uint32_t)srt_a[m]], ocn[(uint32_t)srt_a[m]], zm->is_match, zm->non_homopolymer_errors); + + for (k = f = 0; (k < bk) && (ze > ol->list[(uint32_t)srt_b[k]].x_pos_s); k++) { + rm = &(ol->list[(uint32_t)srt_b[k]]); + // fprintf(stderr, "[M::%s::rm::oid->%u]\ttid::%u\t%.*s\tq::[%u,%u)\tt::[%u,%u)\tsc::%u\tcn::%u\tis_match::%u\terr::%u\n", __func__, (uint32_t)srt_b[k], + // rm->y_id, (int)Get_NAME_LENGTH(R_INF, rm->y_id), Get_NAME(R_INF, rm->y_id), + // rm->x_pos_s, rm->x_pos_e + 1, rm->y_pos_s, rm->y_pos_e + 1, + // osc[(uint32_t)srt_b[k]], ocn[(uint32_t)srt_b[k]], rm->is_match, rm->non_homopolymer_errors); + if(osc[(uint32_t)srt_b[k]] < zsc) continue; + if(ocn[(uint32_t)srt_b[k]] < zcn) continue; + rs = rm->x_pos_s; re = rm->x_pos_e + 1; + os = ((rs>=zs)?rs:zs); oe = ((re<=ze)?re:ze); + if((oe > os) && (oe - os) >= ob) { + oi = rm->overlapLen; rr = cl->list[oi].readID; kn = 0;///idx has been moved to rm->overlapLen + assert(cl->list[oi].readID == rr); + for (; (oi < cn) && (cl->list[oi].readID == rr) && (kn < zcn); oi++) { + ce = cl->list[oi].self_offset; cs = ce - (cl->list[oi].cnt&(0xffu)); + if((cs >= os) && (ce <= oe)) kn++; + // fprintf(stderr, "[M::%s::rm::oid->%u]\tcs::%lu\tce::%lu\tos::%lu\toe::%lu\tkn::%lu\n", __func__, (uint32_t)srt_b[k], cs, ce, os, oe, kn); + } + if(kn >= zcn) { + f = 1; break; + } + } + } + // fprintf(stderr, "[M::%s::zm::oid->%u]\ttid::%u\t%.*s\tq::[%u,%u)\tt::[%u,%u)\tsc::%u\tcn::%u\tis_match::%u\terr::%u\tf::%u\n", __func__, (uint32_t)srt_a[m], + // zm->y_id, (int)Get_NAME_LENGTH(R_INF, zm->y_id), Get_NAME(R_INF, zm->y_id), + // zm->x_pos_s, zm->x_pos_e + 1, zm->y_pos_s, zm->y_pos_e + 1, + // osc[(uint32_t)srt_a[m]], ocn[(uint32_t)srt_a[m]], zm->is_match, zm->non_homopolymer_errors, f); + if(f) continue; + + zm->non_homopolymer_errors = UINT32_MAX - 1; + if((zm->x_pos_e-zm->x_pos_s) != (zm->y_pos_e-zm->y_pos_s)) continue; + resize_UC_Read(tu, zm->y_pos_e+1-zm->y_pos_s); + recover_UC_Read_sub_region(tu->seq, zm->y_pos_s, zm->y_pos_e+1-zm->y_pos_s, zm->y_pos_strand, rref, zm->y_id); + if(!exact_ec_check(qu->seq, qu->length, tu->seq, zm->y_pos_e+1-zm->y_pos_s, zm->x_pos_s, zm->x_pos_e+1, 0, zm->y_pos_e+1-zm->y_pos_s)) continue; + zm->is_match = 1; zm->strong = 0; zm->without_large_indel = 1; zm->shared_seed = 1; zm->non_homopolymer_errors = 0; + + if(wcut) update_mm_wins(zm, wcut, wcut_n, ocw, ql, 0, 0.15, 16, 0, 1); + } + +} + void gen_hc_aln_small_chn_smp_adv(gen_hc_aln_t *ez, uint64_t ql, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, uint64_t *wcut, uint64_t wcut_n, uint64_t *tot_b, uint64_t *srt_a, uint64_t i_srt_n, uint64_t ch_occ, uint64_t ch_sc) @@ -33593,6 +33707,176 @@ uint64_t gen_hc_r_alin_adp_smp(overlap_region_alloc* ol, Candidates_list *cl, Al return tot_b; } +///need to consider coverage, this information is missing right now (currently only use numbers) +uint64_t gen_hc_r_alin_adp_smp_ff_ec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, int64_t rid, asg64_v *sp, uint64_t ocw, uint32_t *ocn, uint32_t *osc, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min) +{ + uint64_t i, k, ql = qu->length, *wcut = NULL, wcut_n = 0, *wsrt = NULL, wsrt_n = 0, lch = 0, wsrt_n0, spn0 = 0, tz, tot_b = 0, t_cov0, z_cov, focv = 0, fcov_0 = 0; + overlap_region *z, t; uint8_t fc = ((ol->length > max_n_chain)?(1):(0)), wf, ff; int32_t s[4]; s[0] = s[1] = s[2] = s[3] = 0; uint32_t scn[4]; + scn[0] = scn[1] = scn[2] = scn[3] = 0; + + ol->mapped_overlaps_length = 0; + if(ol->length <= 0 || ql <= 0) return tot_b; + + // prt_chain_cluster(ol, cl, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 0, NULL); + + spn0 = sp->n; sp->n += ol->length + ol->length; + kv_resize(uint64_t, (*sp), sp->n); + wsrt = sp->a + spn0; + if(fc) { + wcut = infer_chn_bar_0(ol, max_n_chain, ocn, osc, ql, ocw, s, sp, &wcut_n, &focv); + if ((s[0] <= 0) && (s[1] <= 0) && (s[2] <= 0) && (s[3] <= 0)) fc = 0; + wsrt = sp->a + spn0; + + focv = focv/ql; if(focv < 0) focv = 1; focv *= 1.05; + if(focv < ave_cov_min) focv = ave_cov_min; + if(focv > (max_n_chain<<1)) focv = max_n_chain<<1; + fcov_0 = max_n_chain>>1; + } + + s[0] = s[1] = s[2] = s[3] = 0;///reset it + for (i = wsrt_n = 0; i < ol->length; i++) {///primary chain + z = &(ol->list[i]); + + if(z->is_match == 0) {///the one we need to verify + // z->shared_seed = z->non_homopolymer_errors; + z->shared_seed = 0;/// the chain offset has been kept in z->overlapLen + z->non_homopolymer_errors = UINT32_MAX;///for index + // fprintf(stderr, "[M::%s]\trid::%u(%c)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\n", __func__, + // z->y_id, "+-"[z->y_pos_strand], (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), + // z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1); + } + + wf = ha_ov_type(z, ql); ff = 1; + + if(fc) { + if(((scn[wf] <= max_n_chain) || ((scn[wf] <= (max_n_chain<<1)) && (((int64_t)osc[i]) == s[wf])) || (z->is_match != 0))) { + ff = update_mm_wins(z, wcut, wcut_n, ocw, ql, ((scn[wf]>max_n_chain_f) && (z->is_match == 0))?1:0, 0.15, 16, focv, 0); + } else { + ff = 0; + } + } + + + if((!ff) && (z->is_match == 0)) {///fitered out due to coverage + wsrt[wsrt_n++] = (((uint64_t)osc[i])<<32)|(i); + continue; + } + + scn[wf]++; + if (scn[wf] == max_n_chain) { + s[wf] = osc[i]; + } + + if((ocn[i] < chain_cutoff) && (z->is_match == 0)) {///fitered out due to no enough minimizers + lch = 1; + wsrt[wsrt_n++] = (((uint64_t)-1)<<32)|(i); + continue; + } + + if(z->is_match == ((uint8_t)-1)) continue;///has been fitered out outside this function + + if(z->is_match == 0) { + z->non_homopolymer_errors = UINT32_MAX - 1; + if((z->x_pos_e-z->x_pos_s) != (z->y_pos_e-z->y_pos_s)) continue; + resize_UC_Read(tu, z->y_pos_e+1-z->y_pos_s); + recover_UC_Read_sub_region(tu->seq, z->y_pos_s, z->y_pos_e+1-z->y_pos_s, z->y_pos_strand, rref, z->y_id); + if(!exact_ec_check(qu->seq, qu->length, tu->seq, z->y_pos_e+1-z->y_pos_s, z->x_pos_s, z->x_pos_e+1, 0, z->y_pos_e+1-z->y_pos_s)) continue; + z->is_match = 1; z->strong = 0; z->without_large_indel = 1; z->shared_seed = 1; z->non_homopolymer_errors = 0; + } + + if(fc) update_mm_wins(z, wcut, wcut_n, ocw, ql, 0, 0.15, 16, 0, 1); + } + + // fprintf(stderr, "[M::%s] wsrt_n::%lu, fc::%u, ez->ol->length::%lu, ez->max_n_chain::%lu, n_cu::%lu\n", __func__, wsrt_n, fc, ol->length, max_n_chain, n_cu); + // for (i = 0; i < wsrt_n; i++) { + // z = &(ol->list[(uint32_t)wsrt[i]]); + // fprintf(stderr, "[M::%s::i->%lu]\tqid::%u\ttid::%u\t%.*s\toid::%u\tsc::%lu\n", __func__, i, + // z->x_id, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), + // (uint32_t)wsrt[i], wsrt[i]>>32); + // } + + // fprintf(stderr, "-1-[M::%s] tot_b::%lu\n", __func__, tot_b); + + if(fc && wsrt_n) { + t_cov0 = cal_mm_wins_cov(wcut, wcut_n, ocw, ql, fcov_0); + + wsrt_n0 = wsrt_n; wsrt_n = 0; lch = 0; + for (tz = z_cov = 0; (tz < wsrt_n0) && (z_cov <= t_cov0); tz++) { + i = (uint32_t)wsrt[tz]; z = &(ol->list[i]); + assert(z->non_homopolymer_errors == UINT32_MAX); + assert(z->is_match == 0); + + + if((wsrt[tz]>>32) == ((uint32_t)-1)) {///fitered out due to no enough minimizers + lch = 1; wsrt[wsrt_n++] = (((uint64_t)-1)<<32)|(i); + continue; + } + + + if(update_mm_wins(z, wcut, wcut_n, ocw, ql, 1, 0.15, 16, fcov_0, 0) == 0) { + wsrt[wsrt_n++] = (((uint64_t)osc[i])<<32)|(i); + continue; + } + + if(ocn[i] < chain_cutoff) {///fitered out due to no enough minimizers + lch = 1; wsrt[wsrt_n++] = (((uint64_t)-1)<<32)|(i); + continue; + } + + z->non_homopolymer_errors = UINT32_MAX - 1; + if((z->x_pos_e-z->x_pos_s) != (z->y_pos_e-z->y_pos_s)) continue; + resize_UC_Read(tu, z->y_pos_e+1-z->y_pos_s); z_cov += z->y_pos_e+1-z->y_pos_s; + recover_UC_Read_sub_region(tu->seq, z->y_pos_s, z->y_pos_e+1-z->y_pos_s, z->y_pos_strand, rref, z->y_id); + if(!exact_ec_check(qu->seq, qu->length, tu->seq, z->y_pos_e+1-z->y_pos_s, z->x_pos_s, z->x_pos_e+1, 0, z->y_pos_e+1-z->y_pos_s)) continue; + z->is_match = 1; z->strong = 0; z->without_large_indel = 1; z->shared_seed = 1; z->non_homopolymer_errors = 0; + + if(fc) update_mm_wins(z, wcut, wcut_n, ocw, ql, 0, 0.15, 16, 0, 1); + } + + for (; tz < wsrt_n0; tz++) { + i = (uint32_t)wsrt[tz]; z = &(ol->list[i]); + assert(z->non_homopolymer_errors == UINT32_MAX); + assert(z->is_match == 0); + + if((wsrt[tz]>>32) == ((uint32_t)-1)) {///fitered out due to no enough minimizers + lch = 1; wsrt[wsrt_n++] = (((uint64_t)-1)<<32)|(i); + continue; + } + wsrt[wsrt_n++] = (((uint64_t)osc[i])<<32)|(i); + } + + tot_b += z_cov; + } + + + if(lch) { + gen_hc_aln_small_chn_smp_ff_ec(ol, cl, rref, ql, qu, tu, exz, wcut, wcut_n, ocw, wsrt, wsrt_n, ocn, osc, chain_cutoff, /**3, 3**/16, 16); + } + + + for (i = k = 0; i < ol->length; i++) {///primary chain + z = &(ol->list[i]); + if(z->is_match == 1 || z->is_match == 2) { + if(k != i) { + t = ol->list[k]; + ol->list[k] = ol->list[i]; + ol->list[i] = t; + } + k++; + } + } + // print_mm_wins_all(wcut, wcut_n, ocw, ql, 0); + // print_mm_wins_all(wcut, wcut_n, ocw, ql, 1); + + ol->length = k; + // prt_chain_cluster(ol, cl, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 1, NULL); + // fprintf(stderr, "-[M::%s]\trid::%ld\ttot_b::%lu\tql::%lu\tmax_n_chain::%ld\ttot_b_cov:::%ld\n", __func__, rid, tot_b, ql, max_n_chain, tot_b/ql); + // exit(1); + sp->n = spn0; + // if(ol->length <= 0) return tot_b; + return tot_b; +} + void gen_hc_r_alin_adv(gen_hc_aln_t *ez) diff --git a/Correct.h b/Correct.h index fa1fb5f..83268bf 100644 --- a/Correct.h +++ b/Correct.h @@ -1432,7 +1432,7 @@ 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); - +uint64_t gen_hc_r_alin_adp_smp_ff_ec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, int64_t rid, asg64_v *sp, uint64_t ocw, uint32_t *ocn, uint32_t *osc, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min); uint64_t gen_hc_r_alin_adp_smp(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 sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *sp, uint64_t ocw, uint8_t *hpf, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, uint8_t set_match, uint8_t is_dedup); void gen_hc_r_alin_adv_adp_smp(gen_hc_aln_t *ez, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, uint8_t set_match); @@ -1452,7 +1452,12 @@ void push_alnw(overlap_region *aux_o, bit_extz_t *exz); void cal_exz_global(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez); void get_wqual(uint64_t zid, uint64_t zpos, uint64_t zrev, asg8_v *v, uint8_t *va, uint64_t scw, uint64_t *tqual, uint64_t *wqual); void gen_reseed_re(overlap_region_alloc *ol, Candidates_list *cl, overlap_region *aux_o, overlap_region *rse_o, All_reads *rref, UC_Read* qu, UC_Read *tu, bit_extz_t *exz, kv_ul_ov_t *c_idx, asg64_v *idx, asg64_v *res, int64_t bd, int64_t mzw, int64_t kl, int64_t rid, double err_h, double err_l, asg16_v *b16, uint64_t tqn, uint8_t *hpf); - +inline uint64_t exact_ec_check(char *qstr, uint64_t ql, char *tstr, uint64_t tl, int64_t qs, int64_t qe, int64_t ts, int64_t te) +{ + if(qe - qs != te - ts) return 0; + if(memcmp(qstr + qs, tstr + ts, qe - qs) == 0) return 1; + return 0; +} #define ovlp_id(x) ((x).tn) #define ovlp_min_wid(x) ((x).ts) diff --git a/ecovlp.cpp b/ecovlp.cpp index 204cf31..c402145 100644 --- a/ecovlp.cpp +++ b/ecovlp.cpp @@ -2928,12 +2928,6 @@ void check_well_cal(asg16_v *sc, asg64_v *idx, uint8_t *f_ec, uint8_t *abnormal, } } -inline uint64_t exact_ec_check(char *qstr, uint64_t ql, char *tstr, uint64_t tl, int64_t qs, int64_t qe, int64_t ts, int64_t te) -{ - if(qe - qs != te - ts) return 0; - if(memcmp(qstr + qs, tstr + ts, qe - qs) == 0) return 1; - return 0; -} void gen_ff_hpc(asg8_v *fi, char *a, int64_t an, int64_t hpc_rr, int64_t hpc_cutoff) { @@ -4909,9 +4903,435 @@ static void worker_hap_dc_ec_gen(void *data, long i, int tid) refresh_ec_ovec_buf_t0(b, REFRESH_N); } +inline uint64_t dd_ec_check(All_reads *rref, overlap_region_alloc *ol, uint64_t *oidx, uint64_t oidx_n, uint64_t *eidx, uint64_t eidx_n, uint64_t *e_des, uint32_t *osc, double sc_bd, ma_hit_t_alloc *in0, ma_hit_t_alloc *in1, char* rs, uint64_t rl, UC_Read *tu, double orate, uint8_t *push_oidx) +{ + overlap_region *z; uint64_t ok, ek, tid, trev, is_match, diff, strong, without_large_indel, shared_seed; int64_t om; ma_hit_t *p = NULL; uint64_t aq[2], at[2], bq[2], bt[2], ovlp, os, oe; + uint64_t is_match_0, diff_0, strong_0, without_large_indel_0, shared_seed_0; int64_t mm_sc = INT32_MIN, sc, mm_ok; uint64_t mm_ol = ((uint64_t)-1), db_eidx = 0; uint8_t push_eidx = 0; uint64_t e_des_n = 0; + + + if(oidx_n && eidx_n) { + is_match_0 = diff_0 = strong_0 = without_large_indel_0 = shared_seed_0 = (uint64_t)-1; mm_ok = -1; db_eidx = 0; + for (ok = 0; ok < oidx_n; ok++) { + z = &(ol->list[(uint32_t)oidx[ok]]); tid = z->y_id; trev = z->y_pos_strand; + z->non_homopolymer_errors = 0; + z->shared_seed = 0; z->is_match = 0;///may have bugs? + // fprintf(stderr, "[M::%s]\trid::%u(%c)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\n", __func__, + // z->y_id, "+-"[z->y_pos_strand], (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), + // z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1); + for (ek = 0; ek < eidx_n; ek++) { + om = 1; ///z->shared_seed = 0; z->is_match = 0; + if(eidx[ek]&1) { + p = &(in1->buffer[((uint32_t)eidx[ek])>>1]); is_match = 2; + } else { + p = &(in0->buffer[((uint32_t)eidx[ek])>>1]); is_match = 1; + if(p->el) db_eidx++; + } + diff = 0; strong = p->ml; without_large_indel = p->no_l_indel; shared_seed = 0;///indicate exact match + ///z->strong = p->ml; z->without_large_indel = p->no_l_indel; + + + aq[0] = (uint32_t)p->qns; aq[1] = p->qe; + at[0] = p->ts; at[1] = p->te; + + bq[0] = z->x_pos_s; bq[1] = z->x_pos_e + 1; + bt[0] = z->y_pos_s; bt[1] = z->y_pos_e + 1; + // fprintf(stderr, "[a]\tq::[%lu,%lu)\tt::[%lu,%lu)\n", aq[0], aq[1], at[0], at[1]); + // fprintf(stderr, "[b]\tq::[%lu,%lu)\tt::[%lu,%lu)\n", bq[0], bq[1], bt[0], bt[1]); + + os = MAX(aq[0], bq[0]); oe = MIN(aq[1], bq[1]); + ovlp = ((oe>os)? (oe-os):0); + if(!((ovlp) && (ovlp >= ((aq[1] - aq[0])*orate)) && ((ovlp >= ((bq[1] - bq[0])*orate))))) om = 0; + diff += ((aq[1] - aq[0]) - ovlp) + ((bq[1] - bq[0]) - ovlp); + // fprintf(stderr, "[q]\tqos::[%lu,%lu)\tovlp::%lu\tdiff::%lu\n", os, oe, ovlp, diff); + // z->non_homopolymer_errors += ((aq[1] - aq[0]) - ovlp); + + os = MAX(at[0], bt[0]); oe = MIN(at[1], bt[1]); + ovlp = ((oe>os)? (oe-os):0); + if(!((ovlp) && (ovlp >= ((at[1] - at[0])*orate)) && ((ovlp >= ((bt[1] - bt[0])*orate))))) om = 0; + diff += ((at[1] - at[0]) - ovlp) + ((bt[1] - bt[0]) - ovlp); + // fprintf(stderr, "[t]\ttos::[%lu,%lu)\tovlp::%lu\tdiff::%lu\n", os, oe, ovlp, diff); + // z->non_homopolymer_errors += ((at[1] - at[0]) - ovlp); + // fprintf(stderr, "[M::%s]\trid::%u(%c)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\ttype::%lu\tel::%u\tdiff::%lu\n", __func__, p->tn, "+-"[p->rev], (int)Get_NAME_LENGTH(R_INF, p->tn), Get_NAME(R_INF, p->tn), + // (uint32_t)p->qns, p->qe, p->ts, p->te, is_match, ((p->el) && (is_match == 1))?1:0, diff); + + if(om) { + // if(is_match == 1 && p->el == 1) p->el = 0; + if((bq[1] - bq[0]) == (bt[1] - bt[0])) { + resize_UC_Read(tu, bt[1] - bt[0]); + recover_UC_Read_sub_region(tu->seq, bt[0], bt[1] - bt[0], trev, rref, tid); + if(exact_ec_check(rs, rl, tu->seq, bt[1] - bt[0], bq[0], bq[1], 0, bt[1] - bt[0])) { + if(is_match == 2) { + strong = 0; without_large_indel = 1; + // z->strong = 0; z->without_large_indel = 1; + } + is_match = 1; shared_seed = 1;///z->shared_seed = 1; + } + } + // z->is_match = is_match; + } else { + is_match = 0; + } + + if(om) { + sc = diff; sc = - sc; sc *= sc_bd; sc += osc[(uint32_t)oidx[ok]]; + // fprintf(stderr, "[M::%s]\tsc::%ld\tdiff::%lu\tosc[]::%u\n", __func__, sc, diff, osc[(uint32_t)oidx[ok]]); + if((sc > mm_sc) || ((sc == mm_sc) && (mm_ol < (z->x_pos_e+1-z->x_pos_s)))) { + mm_sc = sc; mm_ol = z->x_pos_e+1-z->x_pos_s; ///mm_k = s; + is_match_0 = is_match; + diff_0 = diff; + strong_0 = strong; + without_large_indel_0 = without_large_indel; + shared_seed_0 = shared_seed; + mm_ok = ok; + } + } + } + } + + if(mm_ok != -1) {///found an intersection + z = &(ol->list[(uint32_t)oidx[mm_ok]]); + z->is_match = is_match_0; + z->non_homopolymer_errors = diff_0; + z->strong = strong_0; + z->without_large_indel = without_large_indel_0; + z->shared_seed = shared_seed_0; + // fprintf(stderr, "mm_ok->[M::%s]\trid::%u(%c)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\n", __func__, + // z->y_id, "+-"[z->y_pos_strand], (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), + // z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1); + if(oidx_n > 1) { + for (ok = 0; ok < oidx_n; ok++) { + if(ok == ((uint64_t)mm_ok)) continue; + z = &(ol->list[(uint32_t)oidx[ok]]); + z->is_match = (uint8_t)-1;///skipped + } + } + + if(db_eidx && eidx_n > 1 && e_des) { + z = &(ol->list[(uint32_t)oidx[mm_ok]]); + for (ek = 0; ek < eidx_n; ek++) { + if(eidx[ek]&1) { + p = &(in1->buffer[((uint32_t)eidx[ek])>>1]); is_match = 2; + } else { + p = &(in0->buffer[((uint32_t)eidx[ek])>>1]); is_match = 1; + } + + if(is_match == 1 && p->el) { + aq[0] = (uint32_t)p->qns; aq[1] = p->qe; + at[0] = p->ts; at[1] = p->te; + + bq[0] = z->x_pos_s; bq[1] = z->x_pos_e + 1; + bt[0] = z->y_pos_s; bt[1] = z->y_pos_e + 1; + + os = MAX(aq[0], bq[0]); oe = MIN(aq[1], bq[1]); + ovlp = ((oe>os)? (oe-os):0); + if(!((ovlp) && (ovlp >= ((aq[1] - aq[0])*orate)) && ((ovlp >= ((bq[1] - bq[0])*orate))))) { + e_des[e_des_n++] = eidx[ek]; + continue; + } + + os = MAX(at[0], bt[0]); oe = MIN(at[1], bt[1]); + ovlp = ((oe>os)? (oe-os):0); + if(!((ovlp) && (ovlp >= ((at[1] - at[0])*orate)) && ((ovlp >= ((bt[1] - bt[0])*orate))))) { + e_des[e_des_n++] = eidx[ek]; + continue; + } + } + } + } + } else {///no intersection, need to double check both sides + push_eidx = 1; (*push_oidx) = 1; + } + } else if(oidx_n) {///in this case, if ol->list is non-exact matched, directly discard them; so no need to to consider the score of chain + for (ok = 0; ok < oidx_n; ok++) { + z = &(ol->list[(uint32_t)oidx[ok]]); + z->non_homopolymer_errors = 0; + z->shared_seed = 0; z->is_match = 0;///may have bugs? + } + (*push_oidx) = 1; + } else { + push_eidx = 1; + } + + + if(e_des && eidx && push_eidx) { + for (ek = 0; ek < eidx_n; ek++) { + if(eidx[ek]&1) { + p = &(in1->buffer[((uint32_t)eidx[ek])>>1]); is_match = 2; + } else { + p = &(in0->buffer[((uint32_t)eidx[ek])>>1]); is_match = 1; + } + if(is_match == 1 && p->el) { + e_des[e_des_n++] = eidx[ek]; + } + } + } + + return e_des_n; +} + +inline void push_overlap_region_by_ma_hit_t_0(overlap_region_alloc *ol, ma_hit_t *p, uint32_t rid, bit_extz_t *exz) +{ + overlap_region *z = NULL; + kv_pushp_ol(overlap_region, (*ol), &z); + clear_fake_cigar(&(z->f_cigar)); + clear_window_list_alloc(&(z->w_list)); + clear_window_list_alloc(&(z->boundary_cigars)); + + z->x_id = rid; z->y_id = p->tn; + z->x_pos_strand = 0; z->y_pos_strand = p->rev; + z->x_pos_s = (uint32_t)p->qns; + z->x_pos_e = p->qe - 1; + z->y_pos_s = p->ts; + z->y_pos_e = p->te - 1; + + z->align_length = z->overlapLen = z->x_pos_e + 1 - z->x_pos_s; + z->non_homopolymer_errors = 0; + + z->is_match = 1; z->shared_seed = 1; + z->strong = p->ml; z->without_large_indel = p->no_l_indel; + + set_exact_exz(exz, z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1); push_alnw(z, exz); +} + +void h_ec_lchain_fast_new_adv(uint32_t rid, UC_Read *qu, UC_Read *tu, All_reads *rref, overlap_region_alloc *ol, Candidates_list *cl, bit_extz_t *exz, + ma_hit_t_alloc *in0, ma_hit_t_alloc *in1, double orate, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, int64_t wl, asg64_v *srt, uint64_t ocw, asg32_v *v32) +{ + if(ol->length <= 0) return; + uint32_t *a_cu = NULL, *a_ci = NULL, *ocn = NULL, *osc = NULL; uint64_t *idx_cu = NULL, n_cu = 0; + uint64_t ok, o_l, ek, e_l, *osrt_a, osrt_n, *esrt_a, esrt_n, esrt_m, oid, eid; v32->n = 0; + uint64_t m, *ei, en, *oi, on, m0; overlap_region *z; ma_hit_t *oa = NULL, *p = NULL; char* rs = qu->seq; uint64_t rl = qu->length; uint8_t push_oidx = 0; + ggen_chain_clus_0(ol, v32, NULL, &a_cu, &a_ci, &ocn, &osc, &idx_cu, &n_cu); + + srt->n = 0; kv_resize(uint64_t, *srt, (in0->length + in1->length)); + oa = in0->buffer; on = in0->length; m0 = 0; + for (ok = 0; ok < on; ok++) { + m = oa[ok].tn; m <<= 1; m |= oa[ok].rev; + m <<= 32; m |= (ok<<1); m |= m0; + kv_push(uint64_t, (*srt), m); + } + oa = in1->buffer; on = in1->length; m0 = 1; + for (ok = 0; ok < on; ok++) { + m = oa[ok].tn; m <<= 1; m |= oa[ok].rev; + m <<= 32; m |= (ok<<1); m |= m0; + kv_push(uint64_t, (*srt), m); + } + + ///for on-the-fly chains + kv_resize(uint64_t, *srt, (srt->n + ol->length)); + ei = srt->a; en = srt->n; oi = srt->a + srt->n; on = ol->length; + srt->n += ol->length; + for (ok = 0; ok < on; ok++) { + z = &(ol->list[ok]); z->is_match = z->strong = z->without_large_indel = z->shared_seed = 0; + z->overlapLen = z->non_homopolymer_errors; ///just for chain filtering, keeping the start pos of this chain + oi[ok] = z->y_id; oi[ok] <<= 1; oi[ok] |= z->y_pos_strand; + oi[ok] <<= 32; oi[ok] |= ok; + // fprintf(stderr, "[M::%s::ok->%lu]\trid::%u(%c)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\n", __func__, ok, z->y_id, "+-"[z->y_pos_strand], (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), + // z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1); + } + + radix_sort_ec64(ei, ei + en); radix_sort_ec64(oi, oi + on); + + // for (ok = 0; ok < on; ok++) { + // z = &(ol->list[(uint32_t)oi[ok]]); + // fprintf(stderr, "[M::%s::ok->%lu]\trid::%u(%c)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\n", __func__, ok, z->y_id, "+-"[z->y_pos_strand], (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), + // z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1); + // } + // for (ek = 0; ek < en; ek++) { + // if(ei[ek]&1) { + // p = &(in1->buffer[((uint32_t)ei[ek])>>1]); m0 = 2; + // } else { + // p = &(in0->buffer[((uint32_t)ei[ek])>>1]); m0 = 1; + // } + + // fprintf(stderr, "[M::%s::ek->%lu]\trid::%u(%c)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\ttype::%lu\tel::%u\n", __func__, ek, p->tn, "+-"[p->rev], (int)Get_NAME_LENGTH(R_INF, p->tn), Get_NAME(R_INF, p->tn), + // (uint32_t)p->qns, p->qe, p->ts, p->te, m0, ((p->el) && (m0 == 1))?1:0); + // } + + + + ok = 1; o_l = 0; ek = 1; e_l = 0; oid = eid = (uint64_t)-1; + osrt_a = NULL; osrt_n = 0; esrt_a = NULL; esrt_n = 0; + esrt_m = 0; ///keep overlaps in in0 that are exactly the match but not covered by ol + while ((ok <= on) || (ek <= en)) { + // uint64_t o_l0 = o_l, e_l0 = e_l, dbg_k;; + if(!osrt_a) { + for (osrt_a = oi + o_l, osrt_n = 0; ok <= on; ok++) { + if (ok == on || (oi[ok]>>32) != (oi[o_l]>>32)) { + osrt_n = ok - o_l; oid = oi[o_l]>>32; + // fprintf(stderr, "+[M::%s::okl]\to_l::%lu\tok::%lu\ton::%lu\n", __func__, o_l, ok, on); + // fprintf(stderr, "-[M::%s::okl]\to_l::%lu\tok::%lu\ton::%lu\toid::%lu\n", __func__, o_l, ok, on, oid); + o_l = ok; ok++; + break; + } + } + } + if(!esrt_a) { + for (esrt_a = ei + e_l, esrt_n = 0; ek <= en; ek++) { + if (ek == en || (ei[ek]>>32) != (ei[e_l]>>32)) { + esrt_n = ek - e_l; eid = ei[e_l]>>32; + // fprintf(stderr, "+[M::%s::ekl]\te_l::%lu\tek::%lu\ten::%lu\n", __func__, e_l, ek, en); + // fprintf(stderr, "-[M::%s::ekl]\te_l::%lu\tek::%lu\ten::%lu\teid::%lu\n", __func__, e_l, ek, en, eid); + e_l = ek; ek++; + break; + } + } + } + + oid = eid = (uint64_t)-1; + if(osrt_n) oid = osrt_a[0]>>32; + if(esrt_n) eid = esrt_a[0]>>32; + + // if(osrt_n) { + // for (dbg_k = 0; dbg_k < osrt_n; dbg_k++) { + // z = &(ol->list[(uint32_t)osrt_a[dbg_k]]); + // fprintf(stderr, "[M::%s::ok->%lu]\trid::%u(%c)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\n", __func__, dbg_k + o_l0, z->y_id, "+-"[z->y_pos_strand], (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), + // z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1); + // } + // } + // if(esrt_n) { + // for (dbg_k = 0; dbg_k < esrt_n; dbg_k++) { + // if(esrt_a[dbg_k]&1) { + // p = &(in1->buffer[((uint32_t)esrt_a[dbg_k])>>1]); m0 = 2; + // } else { + // p = &(in0->buffer[((uint32_t)esrt_a[dbg_k])>>1]); m0 = 1; + // } + // fprintf(stderr, "[M::%s::ek->%lu]\trid::%u(%c)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\ttype::%lu\tel::%u\n", __func__, dbg_k + e_l0, p->tn, "+-"[p->rev], (int)Get_NAME_LENGTH(R_INF, p->tn), Get_NAME(R_INF, p->tn), + // (uint32_t)p->qns, p->qe, p->ts, p->te, m0, ((p->el) && (m0 == 1))?1:0); + // } + // } + // fprintf(stderr, "[M::%s]\toid::%lu(%c)\toidx::[%lu(%d),%lu)\to_l::%lu\teid::%lu(%c)\teidx::[%lu(%d),%lu)\te_l::%lu\n\n", __func__, oid>>1, "+-"[oid&1], o_l0, (int32_t)(osrt_a?osrt_a-oi:-1), o_l0 + osrt_n, o_l, + // eid>>1, "+-"[eid&1], e_l0, (int32_t)(esrt_a?esrt_a-ei:-1), e_l0 + esrt_n, e_l); + + + if(osrt_n && esrt_n) { + if(oid == eid) { + esrt_m += dd_ec_check(rref, ol, osrt_a, osrt_n, esrt_a, esrt_n, ei + esrt_m, osc, 8, in0, in1, rs, rl, tu, orate, &push_oidx); + osrt_a = NULL; osrt_n = 0; esrt_a = NULL; esrt_n = 0; + } else if(oid < eid) { + dd_ec_check(rref, ol, osrt_a, osrt_n, NULL, 0, NULL, osc, 8, in0, in1, rs, rl, tu, orate, &push_oidx); + osrt_a = NULL; osrt_n = 0; + } else { + esrt_m += dd_ec_check(rref, ol, NULL, 0, esrt_a, esrt_n, ei + esrt_m, osc, 8, in0, in1, rs, rl, tu, orate, &push_oidx); + esrt_a = NULL; esrt_n = 0; + } + } else if(osrt_n) {///ek > en + dd_ec_check(rref, ol, osrt_a, osrt_n, NULL, 0, NULL, osc, 8, in0, in1, rs, rl, tu, orate, &push_oidx); + osrt_a = NULL; osrt_n = 0; + } else if(esrt_n) {///ok > on + esrt_m += dd_ec_check(rref, ol, NULL, 0, esrt_a, esrt_n, ei + esrt_m, osc, 8, in0, in1, rs, rl, tu, orate, &push_oidx); + esrt_a = NULL; esrt_n = 0; + } + // if(esrt_m) fprintf(stderr, "[M::%s] add-esrt_m::%lu\n", __func__, esrt_m); + } + + // fprintf(stderr, "[M::%s] rid::%u, esrt_m::%lu, en::%lu\n", __func__, rid, esrt_m, en); + + + if(push_oidx) {////go for verify + gen_hc_r_alin_adp_smp_ff_ec(ol, cl, rref, qu, tu, exz, rid, srt, ocw, ocn, osc, max_n_chain>0?max_n_chain:1, max_n_chain_f>0?max_n_chain_f:1, chain_cutoff, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1)); + ei = srt->a; en = srt->n - ol->length; + oi = srt->a + srt->n - ol->length; on = ol->length; + } + + if(esrt_m) {///exact matches but not intersected with ol + for (ek = 0; ek < esrt_m; ek++) { + // if(!((ei[ek]&1) == 0)) { + // fprintf(stderr, "[M::%s] rid::%u, esrt_m::%lu, en::%lu, ek::%lu\n", __func__, rid, esrt_m, en, ek); + // } + assert((ei[ek]&1) == 0); + p = &(in0->buffer[((uint32_t)ei[ek])>>1]); + assert(p->el == 1); + push_overlap_region_by_ma_hit_t_0(ol, p, rid, exz); + } + } + + overlap_region_sort_y_id(ol->list, ol->length); + + if(ol->length > 1) {///for duplicated chains + uint64_t mm_k, s, k, l; int64_t mm_sc, sc; overlap_region t; + for (k = 1, l = m = 0; k <= ol->length; k++) { + if(k == ol->length || ol->list[k].y_id != ol->list[l].y_id) { + mm_k = l; + if(k - l > 1) { + for (s = l, mm_sc = INT32_MIN, mm_k = ((uint64_t)-1); s < k; s++) { + z = &(ol->list[s]); + sc = z->non_homopolymer_errors; sc = - sc; + if((sc > mm_sc) || ((sc == mm_sc) && ((ol->list[mm_k].x_pos_e+1-ol->list[mm_k].x_pos_s) < (z->x_pos_e+1-z->x_pos_s)))) { + mm_sc = sc; mm_k = s; + } + } + } + if(mm_k != ((uint64_t)-1)) { + if(mm_k != m) { + t = ol->list[mm_k]; + ol->list[mm_k] = ol->list[m]; + ol->list[m] = t; + } + m++; + } + l = k; + } + } + ol->length = m; + } +} + + static void worker_hap_dc_ec_gen_new_idx(void *data, long i, int tid) { + ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); + uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO); + uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; uint32_t qlen = 0, qw = 0; + // if(i != 4318) return; + + recover_UC_Read(&b->self_read, &R_INF, i); qlen = b->self_read.length; + qw = ((qlen < (COV_W_AC<<1))?(qlen>>1):(COV_W_AC)); if(!qw) qw = 1; + + // h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, /**0.02**/0.001, asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32, COV_W, 0); + h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, /**0.02**/0.001, asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32, COV_W, 1); + + + // overlap_region_sort_y_id(b->olist.list, b->olist.length); + // ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length); + // ma_hit_sort_tn(R_INF.reverse_paf[i].buffer, R_INF.reverse_paf[i].length); + + // R_INF.paf[i].is_fully_corrected = is_well_cal(&b->v64, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), b->self_read.length, 4); + + // R_INF.paf[i].is_abnormal = abnormal; + // R_INF.trio_flag[i] = AMBIGU; + + // overlap_region_sort_y_id(b->olist.list, b->olist.length); + // h_ec_lchain_fast_new(b->ab, i, &b->self_read, &b->ovlp_read, &R_INF, &b->olist, &b->clist, &b->exz, &b->v16, &b->v64, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), 0.866666); + h_ec_lchain_fast_new_adv(i, &b->self_read, &b->ovlp_read, &R_INF, &b->olist, &b->clist, &b->exz, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), 0.866666, + asm_opt.max_n_chain, asm_opt.max_n_chain*HC_MF_R, asm_opt.chn_occ, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), &b->v64, qw, &b->v32); + + if((asm_opt.is_ont) && (is_uncorrected_read(&b->olist, &b->v64, qlen, 1600))) { + b->olist.length = 0; + // fprintf(stderr, "[M::%s] rid::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + } + + push_ff_ovlp(&(R_INF.paf[i]), &b->olist, 1, &R_INF, b->cnt); + push_ff_ovlp(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, b->cnt); + + /** + copy_asg_arr(buf0, b->sp); + if(gen_hap_dc_cov(&(b->v64), &buf0, &(R_INF.paf[i]), &R_INF, WINDOW_HC_FAST, 4, 0.7, &b->self_read, &b->ovlp_read, &scc, i)) { + scc.f[i] = 1; b->num_read_base++; + // fprintf(stderr, "-mm-[M::%s]\tqn::%u::%.*s\n", __func__, (uint32_t)(i), (int)Get_NAME_LENGTH(R_INF, i), Get_NAME((R_INF), i)); + } else { + scc.f[i] = 0; b->num_correct_base++; + // fprintf(stderr, "-um-[M::%s]\tqn::%u::%.*s\n", __func__, (uint32_t)(i), (int)Get_NAME_LENGTH(R_INF, i), Get_NAME((R_INF), i)); + } + copy_asg_arr(b->sp, buf0); + **/ + refresh_ec_ovec_buf_t0(b, REFRESH_N); +// exit(1); +} + + void worker_hap_dc_ec_gen_new_idx_back(void *data, long i, int tid) +{ + ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO); uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; uint32_t qlen = 0;