From 1d0562ebd183911f91359aa888dbb46d1b8b9902 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Sun, 9 Oct 2022 16:59:54 -0400 Subject: [PATCH] done second alignment --- Correct.cpp | 46 +++++++++- Hash_Table.cpp | 244 ++++++++++++++++++++++++++++++++++++++++++++++++- Hash_Table.h | 13 +++ anchor.cpp | 73 ++++++++++++++- inter.cpp | 75 ++++++++++----- 5 files changed, 421 insertions(+), 30 deletions(-) diff --git a/Correct.cpp b/Correct.cpp index 647d7fe..ae0638e 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -14973,6 +14973,8 @@ uint64_t query_gen_gov_idx(asg64_v *ovidx, uint64_t v, uint64_t w) ///[s, e) uint64_t gen_region_phase(overlap_region* ol, uint64_t *id_a, uint64_t id_n, uint64_t s, uint64_t e, uint64_t dp, ul_ov_t *c_idx, uint64_t *buf, asg64_v *ovidx) { + #define id_mm ((uint64_t)0x7fffffffffffffff) + #define id_set ((uint64_t)0x8000000000000000) if(!id_n) return id_n; uint64_t k, m, mn, q[2], buf_n, rm_n, i; int64_t err, msc, msc_k, msc_n; overlap_region *z; ul_ov_t *p; @@ -15031,14 +15033,47 @@ uint64_t gen_region_phase(overlap_region* ol, uint64_t *id_a, uint64_t id_n, uin buf[k] = ovlp_id((c_idx[(uint32_t)buf[k]])); // if(k < mn) fprintf(stderr, "d::bst::[M::%s::utg%.6dl]\n", __func__, (int32_t)ol[buf[k]].y_id+1); } + + for (k = 0; k < mn; k++) { + if(buf[k]&id_set) continue; + for (i = 0; i < k; i++) { + if(query_gen_gov_idx(ovidx, buf[k]&id_mm, buf[i]&id_mm)) break; + } + if(i < k) { + if(!(buf[k]&id_set)) { + ol[buf[k]&id_mm].overlapLen -= e - s; + ol[buf[k]&id_mm].align_length -= e - s; + } + buf[k] |= id_set; + + if(!(buf[i]&id_set)) { + ol[buf[i]&id_mm].overlapLen -= e - s; + ol[buf[i]&id_mm].align_length -= e - s; + } + buf[i] |= id_set; + } + } + + + for (k = mn; k < buf_n; k++) { - z = &(ol[buf[k]]); + z = &(ol[buf[k]&id_mm]); z->align_length -= e - s; for (i = 0; i < mn; i++) { - if(query_gen_gov_idx(ovidx, buf[k], buf[i])) break; + if(query_gen_gov_idx(ovidx, buf[k]&id_mm, buf[i]&id_mm)) break; } if(i < mn) { z->align_length += e - s; + if(!(buf[k]&id_set)) { + ol[buf[k]&id_mm].overlapLen -= e - s; + ol[buf[k]&id_mm].align_length -= e - s; + } + buf[k] |= id_set; + if(!(buf[i]&id_set)) { + ol[buf[i]&id_mm].overlapLen -= e - s; + ol[buf[i]&id_mm].align_length -= e - s; + } + buf[i] |= id_set; // fprintf(stderr, "i::bst::[M::%s::utg%.6dl]\n", __func__, (int32_t)ol[buf[k]].y_id+1); } } @@ -15328,7 +15363,12 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t } ///hap->length - + // for (k = 0; k < on; k++) { + // z = &(ol->list[k]); + // if(z->align_length == z->overlapLen) {///prefer alignments without any trans hit + // z->align_length = z->overlapLen = z->x_pos_e+1-z->x_pos_s; + // } + // } } void ul_gap_filling_adv(overlap_region_alloc* ol, Candidates_list *cl, kv_ul_ov_t *aln, uint64_t wl, diff --git a/Hash_Table.cpp b/Hash_Table.cpp index aaf47e7..06f9d1b 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -10,6 +10,12 @@ pthread_mutex_t output_mutex; #define overlap_region_key(a) ((a).y_id) KRADIX_SORT_INIT(overlap_region_sort, overlap_region, overlap_region_key, member_size(overlap_region, y_id)) +#define generic_key(x) (x) +KRADIX_SORT_INIT(hc64i, int64_t, generic_key, 8) + +#define oreg_sss_lt(a, b) ((a).shared_seed > (b).shared_seed) // in the decending order +KSORT_INIT(or_sss, overlap_region, oreg_sss_lt) + #define normal_w(x, y) ((x)>=(y)?(x)/(y):1) void overlap_region_sort_y_id(overlap_region *a, long long n) @@ -1765,7 +1771,10 @@ uint64_t lchain_qdp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, o 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(a_n && a[0].readID == 0) { + // fprintf(stderr, "---[M::%s::utg%.6dl::%c]\n", + // __func__, (int32_t)a[0].readID+1, "+-"[a[0].strand]); + // } if(quick_check) { ret = lchain_qcheck(a, a_n, dp, bw_rate); if (ret > 0) { @@ -1819,6 +1828,11 @@ uint64_t lchain_qdp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, o msc = f[i]; msc_i = i; movl = ovl; } } + // if(a_n && a[0].readID == 0) { + // fprintf(stderr, "i::%ld[M::%s::utg%.6dl::%c] x::%u, y::%u, st::%ld, max_ii::%ld, f[i]::%d, p[i]::%ld, msc_i::%ld, msc::%ld, movl::%ld\n", + // i, __func__, (int32_t)a[i].readID+1, "+-"[a[i].strand], + // a[i].self_offset, a[i].offset, st, max_ii, f[i], p[i], msc_i, msc, movl); + // } } skip_ldp: @@ -1836,11 +1850,237 @@ uint64_t lchain_qdp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, o 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]]; + for (i = 0; i < cL; i++) { + des[i] = a[t[cL-i-1]]; + // if(a_n && a[0].readID == 0) { + // fprintf(stderr, "i::%ld[M::%s::utg%.6dl::%c] x::%u, y::%u, cL::%ld\n", + // i, __func__, (int32_t)des[i].readID+1, "+-"[des[i].strand], des[i].self_offset, des[i].offset, cL); + // } + } return cL; } +#define kv_pushp_ol(type, v, p) do { \ + if ((v).length == (v).size) { \ + (v).list = (type*)realloc((v).list, sizeof(type)*((v).size?((v).size<<1):(2))); \ + memset((v).list+(v).size, 0, sizeof(overlap_region)*(((v).size?((v).size<<1):2)-(v).size));\ + (v).size = (v).size?((v).size<<1):(2); \ + } \ + *(p) = &((v).list[(v).length++]); \ + } while (0) + +void push_ovlp_chain_qgen(overlap_region* o, uint32_t xid, int64_t xl, int64_t yl, int64_t sc, +k_mer_hit *beg, k_mer_hit *end) +{ + int64_t xr, yr; + o->x_id = xid; o->y_id = beg->readID; + o->x_pos_strand = 0; o->y_pos_strand = beg->strand; + o->x_pos_s = beg->self_offset; o->y_pos_s = beg->offset; + o->x_pos_e = end->self_offset; o->y_pos_e = end->offset; + + if(o->x_pos_s <= o->y_pos_s) { + o->y_pos_s -= o->x_pos_s; o->x_pos_s = 0; + } else { + o->x_pos_s -= o->y_pos_s; o->y_pos_s = 0; + } + + xr = xl-o->x_pos_e-1; yr = yl-o->y_pos_e-1; + if(xr <= yr) { + o->x_pos_e = xl-1; o->y_pos_e += xr; + } else { + o->y_pos_e = yl-1; o->x_pos_e += yr; + } + + o->shared_seed = sc; + o->align_length = 0; + o->is_match = 0; + o->non_homopolymer_errors = 0; + o->strong = 0; + o->overlapLen = 0; +} + +int64_t filter_non_ovlp_chains(overlap_region *a, int64_t a_n, int64_t *n_v) +{ + int64_t k, i, n_mchain, omx, omy, opx, opy, ovx, ovy, os, oe; overlap_region *m, *p, t; + for (k = n_mchain = (*n_v) = 0; k < a_n; k++) { + m = &(a[k]); + omx = m->x_pos_e + 1 - m->x_pos_s; + omy = m->y_pos_e + 1 - m->y_pos_s; + for (i = 0; i < n_mchain; i++) { + p = &(a[i]); + opx = p->x_pos_e + 1 - p->x_pos_s; + opy = p->y_pos_e + 1 - p->y_pos_s; + + os = ((m->x_pos_s>=p->x_pos_s)? m->x_pos_s:p->x_pos_s); + oe = ((m->x_pos_e<=p->x_pos_e)? m->x_pos_e:p->x_pos_e) + 1; + ovx = oe>os?oe-os:0; + if((ovx > omx*0.1) || (ovx > opx*0.1)) break; + + os = ((m->y_pos_s>=p->y_pos_s)? m->y_pos_s:p->y_pos_s); + oe = ((m->y_pos_e<=p->y_pos_e)? m->y_pos_e:p->y_pos_e) + 1; + ovy = oe>os?oe-os:0; + if((ovy > omy*0.1) || (ovy > opy*0.1)) break; + } + if(i < n_mchain) continue; + + if (n_mchain != k) { + t = a[k]; a[k] = a[n_mchain]; a[n_mchain] = t; + } + (*n_v) += a[n_mchain].align_length; + n_mchain++; + } + return n_mchain; +} + +uint64_t lchain_qdp_mcopy(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64_t des_idx, + Chain_Data* dp, overlap_region_alloc* res, int64_t max_skip, int64_t max_iter, + int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, + uint32_t xid, int64_t xl, int64_t yl, int64_t quick_check, uint32_t apend_be, + int64_t gen_cigar) +{ + int64_t *p, *t, max_f, n_skip, st, max_j, end_j, sc, msc, msc_i, bw, max_ii, ovl, movl, plus = 0, min_sc, ch_n; + int32_t *f, max, tmp, *ii; int64_t i, k, j, cL = 0; k_mer_hit* a; k_mer_hit* des; k_mer_hit *swap; overlap_region *z; + resize_Chain_Data(dp, a_n, NULL); + t = dp->tmp; f = dp->score; p = dp->pre; ii = dp->occ; + bw = ((xl < yl)?xl:yl); bw *= bw_rate; + msc = msc_i = INT32_MIN; movl = INT32_MAX; ch_n = 1; + a = cl->list + a_idx; des = cl->list + des_idx; + // if(a_n && a[0].readID == 0) { + // fprintf(stderr, "---[M::%s::utg%.6dl::%c]\n", + // __func__, (int32_t)a[0].readID+1, "+-"[a[0].strand]); + // } + + memset(t, 0, (a_n*sizeof((*t)))); + for (i = st = plus = 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; + while (a[i].strand != a[st].strand) ++st; + + 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) || (a[i].self_offset>a[max_ii].self_offset+max_dis) || (a[i].strand!=a[max_ii].strand)) { + max = INT32_MIN; max_ii = -1; + for (j=i-1; (j>=st) && (a[i].self_offset<=max_dis+a[j].self_offset)&&(a[i].strand==a[j].strand); --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) || ((a[i].self_offset<=max_dis+a[max_ii].self_offset)&&(a[i].strand==a[max_ii].strand)&&(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; + } + } + if(f[i] < plus) plus = f[i]; + // if(a_n && a[0].readID == 0) { + // fprintf(stderr, "i::%ld[M::%s::utg%.6dl::%c] x::%u, y::%u, st::%ld, max_ii::%ld, f[i]::%d, p[i]::%ld, msc_i::%ld, msc::%ld, movl::%ld\n", + // i, __func__, (int32_t)a[i].readID+1, "+-"[a[i].strand], + // a[i].self_offset, a[i].offset, st, max_ii, f[i], p[i], msc_i, msc, movl); + // } + } + if((movl < xl) && (movl < yl)) { + msc -= plus; min_sc = msc*0.2; + for (i = ch_n = 0; i < a_n; ++i) {///make all f[] positive + f[i] -= plus; if(i >= ch_n) t[i] = 0; + if(f[i] >= min_sc) { + t[ch_n] = ((uint64_t)f[i])<<32; t[ch_n] += (i<<1); ch_n++; + } + } + + int64_t n_v, n_v0, ni, n_u, n_u0 = res->length; + radix_sort_hc64i(t, t + ch_n); + for (k = ch_n-1, n_v = n_u = 0; k >= 0; --k) { + n_v0 = n_v; + for (i = ((uint32_t)t[k])>>1; i >= 0 && (t[i]&1) == 0; ) { + ii[n_v++] = i; t[i] |= 1; i = p[i]; + } + // if(a_n && a[0].readID == 0) { + // fprintf(stderr, "init_sc::%ld[M::%s::k->%ld] n_v0::%ld, n_v::%ld\n", t[k]>>32, __func__, k, n_v0, n_v); + // } + if(n_v0 == n_v) continue; + sc = (i<0?(t[k]>>32):((t[k]>>32)-f[i])); + if(sc >= min_sc) { + // if(a_n && a[0].readID == 0) { + // fprintf(stderr, "sc::%ld[M::%s::] n_v0::%ld, n_v::%ld, k::%ld, ch_n::%ld, msc::%ld, min_sc::%ld\n", + // sc, __func__, n_v0, n_v, k, ch_n, msc, min_sc); + // } + kv_pushp_ol(overlap_region, (*res), &z); + push_ovlp_chain_qgen(z, xid, xl, yl, sc+plus, &(a[ii[n_v-1]]), &(a[ii[n_v0]])); + z->align_length = n_v-n_v0; z->x_id = n_v0; + n_u++; + } else { + n_v = n_v0; + } + } + + ks_introsort_or_sss(n_u, res->list + n_u0); + res->length = n_u0 + filter_non_ovlp_chains(res->list + n_u0, n_u, &n_v); + n_u = res->length; + + kv_resize_cl(k_mer_hit, (*cl), (n_v+cl->length)); + a = cl->list + a_idx; des = cl->list + des_idx; swap = cl->list + cl->length; + for (k = n_u0, i = n_v0 = n_v = 0; k < n_u; k++) { + z = &(res->list[k]); + z->non_homopolymer_errors = des_idx + i; + n_v0 = z->x_id; ni = z->align_length; + // if(n_u > 1) { + // fprintf(stderr, "\nk::%ld[M::%s::utg%.6dl::%c] ni::%ld\n", + // i, __func__, (int32_t)des[i].readID+1, "+-"[des[i].strand], ni); + // } + for (j = 0; j < ni; j++, i++) { + ///k0 + (ni - j - 1) + swap[i] = a[ii[n_v0 + (ni- j - 1)]]; + swap[i].readID = k; + // if(a_n && a[0].readID == 0) { + // fprintf(stderr, "i::%ld[M::%s::utg%.6dl::%c] x::%u, y::%u, ni::%ld\n", + // i, __func__, (int32_t)swap[i].readID+1, "+-"[swap[i].strand], + // swap[i].self_offset, swap[i].offset, ni); + // } + } + z->x_id = xid; + if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, swap+i-ni, ni); + z->align_length = 0; + } + memcpy(des, swap, i*sizeof((*swap))); //assert(i == ch_n); + return i; + } + ///a[] has been sorted by self_offset + i = msc_i; cL = 0; + while (i >= 0) {t[cL++] = i; i = p[i];} + kv_pushp_ol(overlap_region, (*res), &z); + push_ovlp_chain_qgen(z, xid, xl, yl, msc, &(a[t[cL-1]]), &(a[t[0]])); + for (i = 0; i < cL; i++) {des[i] = a[t[cL-i-1]]; des[i].readID = res->length-1;} + z->non_homopolymer_errors = des_idx; + if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, des, cL); + return cL; +} + #define rev_khit(an, xl, yl) do { \ diff --git a/Hash_Table.h b/Hash_Table.h index 928f860..57fceea 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -204,6 +204,14 @@ int append_utg_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_ *(p) = &((v).list[(v).length++]); \ } while (0) +#define kv_resize_cl(type, v, s) do { \ + if ((v).size < (s)) { \ + (v).size = (s); \ + kv_roundup32((v).size); \ + (v).list = (type*)realloc((v).list, sizeof(type) * (v).size); \ + } \ + } while (0) + #define is_alnw(a) (((a).readID) == ((uint32_t)(0x7fffffff))) #define is_pri_aln(a) ((((a).readID) == ((uint32_t)(0x7fffffff)))||((a).cnt >= (a).readID)) @@ -220,5 +228,10 @@ uint64_t lchain_qdp_fix(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_s int64_t left_fix, int64_t right_fix); uint64_t lchain_simple(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, int64_t max_skip, int64_t max_iter); +uint64_t lchain_qdp_mcopy(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64_t des_idx, + Chain_Data* dp, overlap_region_alloc* res, int64_t max_skip, int64_t max_iter, + int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, + uint32_t xid, int64_t xl, int64_t yl, int64_t quick_check, uint32_t apend_be, + int64_t gen_cigar); #endif diff --git a/anchor.cpp b/anchor.cpp index 79ee111..3cafdd1 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -1112,7 +1112,6 @@ void lchain_qgen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, ui 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)); @@ -1148,6 +1147,12 @@ void lchain_qgen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, ui } } cl->length = m; + // fprintf(stderr, "+[M::%s] cn::%lu, m::%lu, ol->length::%lu\n", __func__, cn, m, ol->length); + // for (k = 0; k < ol->length; k++) { + // fprintf(stderr, "---[M::%s::utg%.6dl] q[%d, %d), t[%d, %d), khit_off::%u\n", __func__, + // (int32_t)ol->list[k].y_id+1, ol->list[k].x_pos_s, ol->list[k].x_pos_e+1, + // ol->list[k].y_pos_s, ol->list[k].y_pos_e+1, ol->list[k].non_homopolymer_errors); + // } k = ol->length; @@ -1183,13 +1188,74 @@ void lchain_qgen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, ui ks_introsort_or_xs(ol->length, ol->list); } + +void lchain_qgen_mcopy(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, 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, cn = cl->length, yid; overlap_region *r; ///srt = 0 + clear_overlap_region_alloc(ol); + + for (l = 0, k = 1, m = 0; k <= cn; k++) { + if((k == cn) || (cl->list[k].readID != cl->list[l].readID)) { + if(cl->list[l].readID != rid) { + yid = cl->list[l].readID; + m += lchain_qdp_mcopy(cl, l, k-l, m, &(cl->chainDP), ol, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate, + rid, rl, rdb?Get_READ_LENGTH((*rdb), yid):udb->ug->u.a[yid].len, quick_check, apend_be, gen_off); + } + l = k; + } + } + cl->length = m; + + // for (k = 0; k < ol->length; k++) { + // fprintf(stderr, "---[M::%s::utg%.6dl] q[%d, %d), t[%d, %d), khit_off::%u\n", __func__, + // (int32_t)ol->list[k].y_id+1, ol->list[k].x_pos_s, ol->list[k].x_pos_e+1, + // ol->list[k].y_pos_s, ol->list[k].y_pos_e+1, ol->list[k].non_homopolymer_errors); + // } + + 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; if(is_accurate) { (*quick_check) = 1; (*max_skip) = 25; (*max_iter) = 5000; (*max_dis) = 5000; div = 0.01; pen_gap = 0.5f; pen_skip = 0.0005f; } else { - (*quick_check) = 1; (*max_skip) = 25; (*max_iter) = 5000; (*max_dis) = 5000; div = 0.1; pen_gap = 0.5f; pen_skip = 0.0005f; + (*quick_check) = 0; (*max_skip) = 25; (*max_iter) = 5000; (*max_dis) = 5000; div = 0.1; pen_gap = 0.5f; pen_skip = 0.0005f; } tmp = expf(-div * (double)mz_k);///0.60049557881 -> HiFi; 0.18268352405 -> ont *chn_pen_gap = pen_gap * tmp;///0.300247789405 -> HiFi; 0.091341762025 -> ont @@ -1208,6 +1274,7 @@ void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t // 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); + // 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 + lchain_qgen_mcopy(cl, overlap_list, rid, rl, NULL, uref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off); } diff --git a/inter.cpp b/inter.cpp index 7e02c1a..b486d1c 100644 --- a/inter.cpp +++ b/inter.cpp @@ -4438,6 +4438,32 @@ int64_t l2g_res_chain(ma_ug_t *ug, ul_ov_t *a, uint64_t a_n, vec_mg_lchain_t *gc return 0; } + +int64_t l2g_res_chain_sc(ma_ug_t *ug, ul_ov_t *a, uint64_t a_n, vec_mg_lchain_t *gchains) +{ + // fprintf(stderr, "[M::%s::] a_n::%lu\n", __func__, a_n); + if(a_n <= 0) return 0; + uint64_t k, m; int64_t l; a_n++; asg_t *g = ug->g; + gchains->n = 0; kv_resize(mg_lchain_t, *gchains, a_n); gchains->n = a_n; + memset(&(gchains->a[0]), 0, sizeof(gchains->a[0])); + gchains->a[0].cnt = a_n - 1; gchains->a[0].v = (uint32_t)-1; + for (k = 1, m = 0, l = 0; k < a_n; k++, m++) { + memset(&(gchains->a[k]), 0, sizeof(gchains->a[k])); + gchains->a[k].v = (a[m].tn<<1)|(a[m].rev); gchains->a[k].dist_pre = -1; + gchains->a[k].off = a[m].qn; gchains->a[k].score = a[m].sec; + gchains->a[k].qs = a[m].qs; gchains->a[k].qe = a[m].qe; + gchains->a[k].rs = a[m].ts; gchains->a[k].re = a[m].te; + if(k > 1) { + gchains->a[k-1].dist_pre = g_adjacent_dis(g, gchains->a[k].v^1, gchains->a[k-1].v^1); + assert(gchains->a[k-1].dist_pre >= 0); + l += g->seq[gchains->a[k-1].v>>1].len + gchains->a[k-1].dist_pre; + } + // fprintf(stderr, "[M::%s::k->%lu] utg%.6dl(%c)\n", __func__, k, (int32_t)(gchains->a[k].v>>1)+1, "+-"[gchains->a[k].v&1]); + } + + return 1; +} + int64_t check_elen_gchain(ul_ov_t *a, int64_t a_n, float trans_thres) { uint32_t sp_e, ep_e, ts, te, tl = 0, el = 0, iel = 0; @@ -5046,7 +5072,8 @@ int64_t debug_i, int64_t tid, void *km) return 1; } -#define aln_sc(a, w) (((int64_t)((a).sec))-((int64_t)(((a).qe-(a).qs-(a).sec)*(w)))) +// #define aln_sc(a, w) (((int64_t)((a).sec))-((int64_t)(((a).qe-(a).qs-(a).sec)*(w)))) +#define aln_sc(a, w) (((int64_t)((a).align_length))-((int64_t)(((a).overlapLen-(a).align_length)*(w)))) void gen_gl_aln(overlap_region_alloc* olist, const ul_idx_t *uref, kv_ul_ov_t *res) { @@ -5057,7 +5084,7 @@ void gen_gl_aln(overlap_region_alloc* olist, const ul_idx_t *uref, kv_ul_ov_t *r // (int32_t)olist->list[k].y_id+1, olist->list[k].x_pos_s, olist->list[k].x_pos_e+1, // olist->list[k].y_pos_s, olist->list[k].y_pos_e+1, // olist->list[k].overlapLen, olist->list[k].align_length); - p = &(res->a[res->n++]); + p = &(res->a[res->n++]); assert(olist->list[k].overlapLen >= olist->list[k].align_length); p->qn = k; p->qs = olist->list[k].x_pos_s; p->qe = olist->list[k].x_pos_e+1; p->tn = olist->list[k].y_id; p->sec = olist->list[k].align_length; p->rev = olist->list[k].y_pos_strand; p->el = (olist->list[k].is_match==1?1:0); @@ -5078,8 +5105,7 @@ void gen_gg_aln(overlap_region_alloc* olist, const ul_idx_t *uref, int64_t trans for (k = 0; k < olist->length; k++) { p = &(res->a[res->n++]); memset(p, 0, sizeof((*p))); p->v = ((olist->list[k].y_id<<1)|(olist->list[k].y_pos_strand)); p->off = k; - p->score = (((int64_t)(olist->list[k].align_length)) - -((int64_t)((olist->list[k].overlapLen-olist->list[k].align_length)*(trans_sc)))); + p->score = aln_sc((olist->list[k]), (trans_sc)); p->qs = olist->list[k].x_pos_s; p->qe = olist->list[k].x_pos_e+1; if((p->v&1)) { p->rs = uref->ug->u.a[p->v>>1].len - (olist->list[k].y_pos_e+1); @@ -5116,6 +5142,7 @@ uint64_t mode, All_reads *ridx, ma_ug_t *ug, int64_t need_srt) memset(t, 0, (res_n*sizeof((*t)))); for (i = st = plus = 0, max_ii = -1; i < res_n; ++i) { + li = &(res->a[i]); li_v = (li->tn<<1)|li->rev; mm_ovlp = mode?max_ovlp_src(uopt, li_v^1):max_ovlp(uref->ug->g, li_v^1); x = (li->qs + mm_ovlp)*diff_ec_ul; @@ -5123,9 +5150,8 @@ uint64_t mode, All_reads *ridx, ma_ug_t *ug, int64_t need_srt) x += li->qs + mm_ovlp; if (x > qlen+1) x = qlen+1; x = find_ul_ov_max(i, res->a, x+G_CHAIN_INDEL); - csc = aln_sc((*li), trans_sc); + csc = aln_sc(ol[(*li).qn], trans_sc); mm_sc = csc; mm_idx = -1; - n_skip = 0; end_j = -1; if ((x-st) > max_iter) st = x-max_iter; for (j = x; j >= st; --j) { // collect potential destination vertices @@ -5145,7 +5171,7 @@ uint64_t mode, All_reads *ridx, ma_ug_t *ug, int64_t need_srt) if (p[j] >= 0) t[p[j]] = i; } } - + end_j = j; if (max_ii < 0 || (res->a[i].qe>(res->a[max_ii].qe+max_dis))) {//too long max = INT32_MIN; max_ii = -1; @@ -5155,16 +5181,16 @@ uint64_t mode, All_reads *ridx, ma_ug_t *ug, int64_t need_srt) } } } - + if (max_ii >= 0 && max_ii < end_j) {///just have a try with a[i]<->a[max_ii] lj = &(res->a[max_ii]); lj_v = (lj->tn<<1)|lj->rev; - if(lj->qe+G_CHAIN_INDEL <= li->qs) break;//even this pair has a overlap, its length will be very small; just ignore - if(lj->qs >= li->qs) continue; - qo = infer_rovlp(li, lj, NULL, NULL, ridx, ug); ///overlap length in query (UL read) - if(li_v != lj_v && get_ecov_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, mode, &share)) { - sc = csc + f[j]; - if(sc > mm_sc) { - mm_sc = sc; mm_idx = max_ii; + if(lj->qe+G_CHAIN_INDEL > li->qs && lj->qs < li->qs) { + qo = infer_rovlp(li, lj, NULL, NULL, ridx, ug); ///overlap length in query (UL read) + if(li_v != lj_v && get_ecov_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, mode, &share)) { + sc = csc + f[j]; + if(sc > mm_sc) { + mm_sc = sc; mm_idx = max_ii; + } } } } @@ -5172,13 +5198,12 @@ uint64_t mode, All_reads *ridx, ma_ug_t *ug, int64_t need_srt) mm_sc = csc; mm_idx = -1; } f[i] = mm_sc; p[i] = mm_idx; - if ((max_ii < 0) || ((res->a[i].qe<=max_dis+res->a[max_ii].qe) && (f[max_ii]tn+1, csc, f[i], p[i], li->qs, li->qe); + // fprintf(stderr, "-5-[M::%s::utg%.6dl] i::%ld, res_n::%ld, csc::%ld, f[i]::%d, p[i]::%ld, q::[%u, %u)\n", + // __func__, (int32_t)li->tn+1, i, res_n, csc, f[i], p[i], li->qs, li->qe); } for (i = 0; i < res_n; ++i) {///make all f[] positive f[i] -= plus; t[i] = ((uint64_t)f[i])<<32; t[i] += (i<<1); @@ -5641,7 +5666,8 @@ int64_t qlen, const ug_opt_t *uopt, int64_t debug_i, int64_t tid, void *km) if(occ) { if(ff_chain(idx, qlen, 0.99/**P_CHAIN_COV**/, -1/**G_CHAIN_TRANS_RATE**/, ll->tk.a, NULL, NULL, NULL, diff_ec_ul, winLen, km)) { - f = l2g_res_chain(uref->ug, ll->tk.a+idx->a[idx->n-1].ts, idx->a[idx->n-1].te-idx->a[idx->n-1].ts, &(gdp->swap), -1/**N_GCHAIN_RATE**/); + f = l2g_res_chain_sc(uref->ug, ll->tk.a+idx->a[idx->n-1].ts, idx->a[idx->n-1].te-idx->a[idx->n-1].ts, &(gdp->swap)); + // f = l2g_res_chain(uref->ug, ll->tk.a+idx->a[idx->n-1].ts, idx->a[idx->n-1].te-idx->a[idx->n-1].ts, &(gdp->swap), -1/**N_GCHAIN_RATE**/); } } // fprintf(stderr, "1-[M::%s] f::%ld\n", __func__, f); @@ -5974,9 +6000,9 @@ uint32_t ck_w_err(overlap_region *z, uint64_t *w_idx, int64_t wl, int64_t ql) } // fprintf(stderr, "[M::%s::utg%.6dl] x::[%u, %u), ol::%ld, e[0]::%ld, e[1]::%ld\n", // __func__, (int32_t)z->y_id+1, z->x_pos_s, z->x_pos_e+1, ol, e[0], e[1]); - if(e[1] > (e[0]+32)) { + if(e[1] > (e[0]+64)) { if(e[1] > (e[0]+(ol*0.01))) return 0; - if(e[1] > (e[0]+(e[0]*0.01))) return 0; + if(e[1] > (e[0]+(e[0]*0.03))) return 0; } // if((e[1] > (e[0]+16)) && (e[1] > (e[0]+(ol*0.01)))) return 0; return 1; @@ -6043,6 +6069,10 @@ void regen_ul_ov_t_lst(const ul_idx_t *uref, overlap_region_alloc* olist, kv_ul_ uint64_t k; ul_ov_t *p; idx->n = 0; kv_resize(ul_ov_t, *idx, olist->length); for (k = 0; k < olist->length; k++) { + // fprintf(stderr, "---[M::%s::utg%.6dl] q[%d, %d), t[%d, %d), tot::%u, cis::%u\n", __func__, + // (int32_t)olist->list[k].y_id+1, olist->list[k].x_pos_s, olist->list[k].x_pos_e+1, + // olist->list[k].y_pos_s, olist->list[k].y_pos_e+1, + // olist->list[k].overlapLen, olist->list[k].align_length); p = &(idx->a[idx->n++]); p->qn = k; p->qs = olist->list[k].x_pos_s; p->qe = olist->list[k].x_pos_e+1; p->tn = olist->list[k].y_id; p->el = 1; p->rev = olist->list[k].y_pos_strand; @@ -6581,7 +6611,8 @@ 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!=319) /**&& (s->id+i!=44) && (s->id+i!=948)**/) return; + // if((s->id+i!=871) && (s->id+i!=963) && (s->id+i!=980)) return; + // if(s->id+i!=963) 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);