diff --git a/CommandLines.h b/CommandLines.h index 4bf61be..bb7eb72 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.25.0-r877" +#define HA_VERSION "0.25.1-r878" #define VERBOSE 0 diff --git a/Correct.cpp b/Correct.cpp index d00367e..e76fe68 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -35989,6 +35989,192 @@ void gen_hc_r_alin_adv_adp_smp(gen_hc_aln_t *ez, uint32_t *a_cu, uint32_t *a_ci, +void gen_hc_r_alin_adv_adp_smp_0(gen_hc_aln_t *ez, uint8_t set_match) +{ + uint64_t i, bs, k, w, ql = ez->qu->length, *wsrt = NULL, wsrt_n = 0, spn0 = 0, tot_b = 0; double err, e_max, align_gap_max, align_gap_rate, chem_drop; overlap_region t, *z; asg64_v *sp = ez->srt; + ez->ol->mapped_overlaps_length = 0; uint32_t *ocn = ez->v32->a, *osc = ez->v32->a + ez->ol->length; + if(ez->ol->length <= 0) return; + + ///base alignment + bs = (MAX(ez->wl[0], ez->wl[1]))+(THRESHOLD_MAX_SIZE<<1)+1; + resize_UC_Read(ez->tu, bs<<1); spn0 = sp->n; + + wsrt = mmp_chn_select(ez->ol, ez->cl, sp, ez->ocw, ocn, osc, ql, &wsrt_n, set_match, ez->max_n_chain, ez->max_n_chain_f, ez->chain_cutoff, ez->ave_cov_min, 0.333333, 16, 16); + + for (i = 0; i < wsrt_n; i++) { + z = &(ez->ol->list[(uint32_t)wsrt[i]]); + + if(z->y_id < ez->t_cut) { + err = ez->e_rate[1]; w = ez->wl[1]; align_gap_max = ez->align_gap_max[1]; align_gap_rate = ez->align_gap_rate[1]; chem_drop = ez->chem_drop[1]; + } else { + err = ez->e_rate[0]; w = ez->wl[0]; align_gap_max = ez->align_gap_max[0]; align_gap_rate = ez->align_gap_rate[0]; chem_drop = ez->chem_drop[0]; + } + e_max = err * 1.5; + + + if((z->is_match == 0) && (!gen_hc_r_alin_flt_1_smp(z, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, align_gap_rate, align_gap_max, NULL/**ez->hpz->a**/, ez->buf, 0, &tot_b))) { + continue; + } + + z->is_match = 1; z->strong = z->without_large_indel = 0; + } + + for (i = k = 0; i < ez->ol->length; i++) {///primary chain + z = &(ez->ol->list[i]); + if(z->is_match == 0) continue; + if(k != i) { + t = ez->ol->list[k]; + ez->ol->list[k] = ez->ol->list[i]; + ez->ol->list[i] = t; + } + k++; + } + + ez->ol->length = k; sp->n = spn0; +} + +void batch_simd_ck(gen_hc_aln_t *ez, uint64_t *wsrt, uint64_t wsrt_n, uint64_t ql, double err, uint64_t w, double align_gap_max, double align_gap_rate, double chem_drop) +{ + if(!wsrt_n) return; + uint8_t rr; uint64_t fi[AVX_GS]; int32_t aux_beg[AVX_GS], aux_end[AVX_GS], t_s[AVX_GS], t_pri_l[AVX_GS]; + uint64_t i, k, zk, in0, s, e, q[2], os, oe; overlap_region *z; + + radix_sort_bc64(wsrt, wsrt + wsrt_n); + for (k = 1, i = 0; k <= wsrt_n; k++) { + if (k == wsrt_n || (wsrt[k]>>32) != (wsrt[i]>>32)) { + if(k - i > 1) { + for (zk = i; zk < k; zk++) { + wsrt[zk] = (((uint64_t)(ez->ol->list[(uint32_t)wsrt[zk]].x_pos_e + 1))<<32)|((uint32_t)wsrt[zk]); + } + radix_sort_bc64(wsrt + i, wsrt + k); + } + i = k; + } + } + + ez->v32->n = (ez->ol->length<<1); in0 = ez->v32->n; + i = 0; s = 0; e = w; e = ((e<=ql)?e:ql); rr = 0; + for (; s < ql; ) {///[s, e) + if(rr) { + for (k = zk = in0; k < ez->v32->n; k++) { + if((wsrt[ez->v32->a[k]]>>32) == UINT32_MAX) continue;///passed + z = &(ez->ol->list[(uint32_t)wsrt[ez->v32->a[k]]]); + q[0] = z->x_pos_s; + q[1] = z->x_pos_e + 1; + ///[s, e) && [q[0], q[1]) + os = MAX(q[0], s); oe = MIN(q[1], e); + if(oe > os) ez->v32->a[zk++] = ez->v32->a[k]; + } + ez->v32->n = zk; + } + + for (; i < wsrt_n; ++i) { + z = &(ez->ol->list[(uint32_t)wsrt[i]]); + q[0] = z->x_pos_s; + q[1] = z->x_pos_e + 1; + if(q[0] >= e) break; + if(z->is_match == 0) { + os = MAX(q[0], s); oe = MIN(q[1], e); + if(oe > os) kv_push(uint32_t, *(ez->v32), i); + } + } + + rr = hc_aln_simd(ez->ol->list, wsrt, ez->v32->a + in0, ez->v32->n - in0, s, e, w, err, ez->rref, ez->qu->seq, ez->tu, fi, aux_beg, aux_end, t_s, t_pri_l, ez->exz, OVERLAP_THRESHOLD_HIFI_FILTER, 0, NULL/**&tot_b**/); + s += w; e += w; e = ((e<=ql)?e:ql); + } +} + +void gen_hc_r_alin_adv_adp_smp_1(gen_hc_aln_t *ez, uint8_t set_match) +{ + uint64_t i, bs, k, w, ql = ez->qu->length, nol_1, *wsrt = NULL, wsrt_n = 0, spn0 = 0; double err, e_max, align_gap_max, align_gap_rate, chem_drop; overlap_region t, *z; asg64_v *sp = ez->srt; + ez->ol->mapped_overlaps_length = 0; uint32_t *ocn = ez->v32->a, *osc = ez->v32->a + ez->ol->length; + if(ez->ol->length <= 0) return; + + ///base alignment + bs = (MAX(ez->wl[0], ez->wl[1]))+(THRESHOLD_MAX_SIZE<<1)+1; + resize_UC_Read(ez->tu, bs<<1); spn0 = sp->n; + + wsrt = mmp_chn_select(ez->ol, ez->cl, sp, ez->ocw, ocn, osc, ql, &wsrt_n, set_match, ez->max_n_chain, ez->max_n_chain_f, ez->chain_cutoff, ez->ave_cov_min, 0.333333, 16, 16); + + + for (i = nol_1 = 0; i < wsrt_n; i++) { + z = &(ez->ol->list[(uint32_t)wsrt[i]]); + if(z->is_match == 0) { + z->w_list.n = 0; z->align_length = 0; + } + wsrt[i] = ((((uint64_t)z->x_pos_s))<<32)|((uint32_t)wsrt[i]); + + if(z->y_id < ez->t_cut) { + if(nol_1 != i) { + bs = wsrt[i]; wsrt[i] = wsrt[nol_1]; wsrt[nol_1] = bs; + } + nol_1++; + } + } + + batch_simd_ck(ez, wsrt, nol_1, ql, ez->e_rate[1], ez->wl[1], ez->align_gap_max[1], ez->align_gap_rate[1], ez->chem_drop[1]); + batch_simd_ck(ez, wsrt + nol_1, wsrt_n - nol_1, ql, ez->e_rate[0], ez->wl[0], ez->align_gap_max[0], ez->align_gap_rate[0], ez->chem_drop[0]); + ocn = ez->v32->a; osc = ez->v32->a + ez->ol->length; + + for (i = 0; i < wsrt_n; i++) { + if((wsrt[i]>>32) == UINT32_MAX) { + continue; + } + z = &(ez->ol->list[(uint32_t)wsrt[i]]); + + if(z->y_id < ez->t_cut) { + err = ez->e_rate[1]; w = ez->wl[1]; align_gap_max = ez->align_gap_max[1]; align_gap_rate = ez->align_gap_rate[1]; chem_drop = ez->chem_drop[1]; + } else { + err = ez->e_rate[0]; w = ez->wl[0]; align_gap_max = ez->align_gap_max[0]; align_gap_rate = ez->align_gap_rate[0]; chem_drop = ez->chem_drop[0]; + } + e_max = err * 1.5; + + if(z->is_match == 0) { + if((!pass_qovlp(z->x_pos_e+1-z->x_pos_s, z->align_length, OVERLAP_THRESHOLD_HIFI_FILTER)) || (!gen_hc_r_alin_flt_1_smp(z, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, + align_gap_rate, align_gap_max, NULL/**ez->hpz->a**/, ez->buf, 1, NULL/**&tot_b**/))) { + continue; + } + } + // if((z->is_match == 0) && (!gen_hc_r_alin_flt_1_smp(z, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, align_gap_rate, align_gap_max, NULL/**ez->hpz->a**/, ez->buf, 0, &tot_b))) { + // continue; + // } + + z->is_match = 1; z->strong = z->without_large_indel = 0; + } + + /** + // uint64_t kqs, kqe, kts, kte, kerr; + for (i = 0; i < wsrt_n; i++) { + if((wsrt[i]>>32) == UINT32_MAX) { + continue; + } + z = &(ol->list[(uint32_t)wsrt[i]]); + + if(z->is_match == 0) { + if((!pass_qovlp(z->x_pos_e+1-z->x_pos_s, z->align_length, OVERLAP_THRESHOLD_HIFI_FILTER)) || (!gen_hc_r_alin_flt_1_smp(z, cl, rref, qu, tu, exz, aux_o, err, e_max, + e_rate, w.window_length, ql, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, 1, &tot_b))) { + continue; + } + } + + z->is_match = 1; z->strong = z->without_large_indel = 0; + } + **/ + + for (i = k = 0; i < ez->ol->length; i++) {///primary chain + z = &(ez->ol->list[i]); + if(z->is_match == 0) continue; + if(k != i) { + t = ez->ol->list[k]; + ez->ol->list[k] = ez->ol->list[i]; + ez->ol->list[i] = t; + } + k++; + } + + ez->ol->length = k; sp->n = spn0; +} + uint64_t 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 sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *kp, uint8_t *hpf) { diff --git a/Correct.h b/Correct.h index 266e0a2..59cbda9 100644 --- a/Correct.h +++ b/Correct.h @@ -1442,6 +1442,8 @@ uint64_t gen_hc_r_alin_adp_mmp_1(overlap_region_alloc* ol, Candidates_list *cl, uint64_t gen_hc_r_alin_adp_mmp_0(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, asg32_v *v32, 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); 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); +void gen_hc_r_alin_adv_adp_smp_0(gen_hc_aln_t *ez, uint8_t set_match); +void gen_hc_r_alin_adv_adp_smp_1(gen_hc_aln_t *ez, uint8_t set_match); void pp_chn_a(overlap_region *z, Candidates_list *cl, uint8_t is_raw); uint64_t 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 sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *kp, uint8_t *hpf); void gen_hc_r_alin_flt(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, overlap_region *aux_b, 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, diff --git a/ecovlp.cpp b/ecovlp.cpp index 6f80ce6..50213fe 100644 --- a/ecovlp.cpp +++ b/ecovlp.cpp @@ -3686,6 +3686,73 @@ void gen_hc_r_alin_ea_adv_flt(gen_hc_aln_t *ez) } +void gen_hc_r_alin_ea_adv_flt_mmp(gen_hc_aln_t *ez) +{ + if(ez->ol->length <= 0) return; + + uint64_t i, k, m, *ei, en, *oi, on, tid, trev, nec; overlap_region *z; ma_hit_t *p; + ez->v32->n = ez->ol->length<<1; kv_resize(uint32_t, *(ez->v32), ez->v32->n); + for (i = 0; i < ez->ol->length; i++) { + ez->v32->a[i] = ez->ol->list[i].align_length; + ez->v32->a[i+ez->ol->length] = ez->ol->list[i].shared_seed; + ez->ol->list[i].align_length = 0; + } + + // uint64_t k, i, m, *ei, en, *oi, on, tid, trev, nec; overlap_region *z; ma_hit_t *p; + for (k = ez->srt->n = 0; k < ez->in->length; k++) { + if(ez->in->buffer[k].el) { + m = ez->in->buffer[k].tn; m <<= 1; m |= ez->in->buffer[k].rev; + m <<= 32; m |= k; kv_push(uint64_t, (*(ez->srt)), m); + } + } + + if(!(ez->srt->n)) { + // gen_hc_r_alin_adv_adp_smp(ez, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 1); + gen_hc_r_alin_adv_adp_smp_1(ez, 1); + } 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); + kv_resize(uint64_t, *(ez->srt), (ez->srt->n + ez->ol->length)); + ei = ez->srt->a; en = ez->srt->n; oi = ez->srt->a + ez->srt->n; on = ez->ol->length; + for (k = 0; k < on; k++) { + z = &(ez->ol->list[k]); z->is_match = z->strong = z->without_large_indel = 0; + oi[k] = z->y_id; oi[k] <<= 1; oi[k] |= z->y_pos_strand; + oi[k] <<= 32; oi[k] |= k; + } + + radix_sort_ec64(ei, ei + en); radix_sort_ec64(oi, oi + on); + for (k = i = nec = 0; k < on; k++) { + z = &(ez->ol->list[(uint32_t)oi[k]]); tid = z->y_id; trev = z->y_pos_strand; + for (; (i < en) && ((ei[i]>>32) < ((tid<<1)|trev)); i++); + if((i < en) && ((ei[i]>>32) == ((tid<<1)|trev))) { + p = &(ez->in->buffer[(uint32_t)ei[i]]); + if((z->x_pos_s == ((uint32_t)p->qns)) && (z->x_pos_e + 1 == p->qe) && + (z->y_pos_s == p->ts) && (z->y_pos_e + 1 == p->te)) { + resize_UC_Read(ez->tu, p->te - p->ts); recover_UC_Read_sub_region(ez->tu->seq, p->ts, p->te - p->ts, trev, ez->rref, tid); + if(exact_ec_check(ez->qu->seq, ez->qu->length, ez->tu->seq, p->te - p->ts, ((uint32_t)p->qns), p->qe, 0, p->te - p->ts)) { + z->is_match = 1; z->shared_seed = z->non_homopolymer_errors;///for index + z->non_homopolymer_errors = 0; z->strong = z->without_large_indel = 0; + set_exact_exz(ez->exz, z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1); push_alnw(z, ez->exz); + nec++; + } + } + } + } + ///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); + + if(on > nec) { + // gen_hc_r_alin_adv_adp_smp(ez, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 0); + gen_hc_r_alin_adv_adp_smp_1(ez, 0); + } + // 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); + } +} + + + void prt_ovlp_sam_0(char *cm, FILE *fp, char *ref_id, int32_t ref_id_n, char *qry_id, int32_t qry_id_n, char *qry_seq, uint64_t qry_seq_n, uint64_t rs, uint64_t re, uint64_t qs, uint64_t qe, uint64_t flag, uint64_t err0, bit_extz_t *ez) { uint64_t ci = 0, err1 = 0; uint16_t c; uint32_t cl, cl0 = 0; char c0 = (char)-1; @@ -5129,14 +5196,15 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid) aux_o = fetch_aux_ovlp(&b->olist, NULL/**&rse_o**/);///must be here copy_asg_arr(buf0, b->sp); - init_gen_hc_aln_t(&ez, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, rse_o, &b->v8q, + init_gen_hc_aln_t(&ez, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, rse_o, /**&b->v8q**/NULL, e_l, e_h, (qid < R_INF.tqn)?(WINDOW_OHC):(WINDOW_HC), WINDOW_OHC, i, E_KHIT, 1, &b->v16, &b->v64, &(R_INF.paf[i]), ((qid < R_INF.tqn)?(1):(0)), 1, (qid < R_INF.tqn)?(0.006):(-1), 0.006, (qid < R_INF.tqn)?(64):(-1), 64, (qid < R_INF.tqn)?(512):(0), (qid < R_INF.tqn)?(6):(0), (qid < R_INF.tqn)?(1.5):(-1), (qid < R_INF.tqn)?(0.1):(-1), NULL, &(b->v32), &buf0, b->ab, (asm_opt.max_n_chain>0)?(asm_opt.max_n_chain):(1), ((asm_opt.max_n_chain*HC_MF_R)>0)?(asm_opt.max_n_chain*HC_MF_R):1, asm_opt.chn_occ, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), qw, R_INF.tqn, asm_opt.hom_cov); // gen_hc_r_alin_ea_adv(&ez); - gen_hc_r_alin_ea_adv_flt(&ez); + // gen_hc_r_alin_ea_adv_flt(&ez); + gen_hc_r_alin_ea_adv_flt_mmp(&ez); copy_asg_arr(b->sp, buf0); // fprintf(stderr, "-c-[M::%s] rid::%ld\n", __func__, i);