From cc16626dfadbef7186b9debb37517df0821a7950 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Thu, 1 Sep 2022 23:16:08 -0400 Subject: [PATCH] debug cigar_offset --- Correct.cpp | 293 +++++++++++++++++++++++-------------------------- Hash_Table.cpp | 83 ++++++++++++-- Hash_Table.h | 2 +- Overlaps.h | 7 ++ anchor.cpp | 52 ++++++--- gchain_map.cpp | 2 +- inter.cpp | 179 +++++++++++++++++++++++++----- 7 files changed, 407 insertions(+), 211 deletions(-) diff --git a/Correct.cpp b/Correct.cpp index 2a62bad..f09da12 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -563,6 +563,11 @@ char* r_string, double max_ov_diff_ec, long long blockLen, long long max_error, continue; } + // if(overlap_list->list[currentID].y_id == 4) { + // fprintf(stderr, "[M::%s] q_s::%lld, t_s::%lld, t_pri_l::%lld, aux_beg::%d, aux_end::%d, aln_l::%lld\n", __func__, + // x_start, y_start, o_len, extra_begin, extra_end, Window_Len); + // } + fill_subregion_ul(dumy->overlap_region_group[groupLen], y_start, o_len, overlap_list->list[currentID].y_pos_strand, uref, overlap_list->list[currentID].y_id, extra_begin, extra_end); @@ -695,6 +700,11 @@ char* r_string, double max_ov_diff_ec, long long blockLen, long long max_error, &extra_begin, &extra_end, &y_start, &o_len)) { continue; } + + // if(overlap_list->list[currentID].y_id == 4) { + // fprintf(stderr, "[M::%s] q_s::%lld, t_s::%lld, t_pri_l::%lld, aux_beg::%d, aux_end::%d, aln_l::%lld\n", __func__, + // x_start, y_start, o_len, extra_begin, extra_end, Window_Len); + // } fill_subregion_ul(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, uref, overlap_list->list[currentID].y_id, extra_begin, extra_end); @@ -712,95 +722,120 @@ char* r_string, double max_ov_diff_ec, long long blockLen, long long max_error, } } -void verify_ul_ll_window(long long window_start, long long window_end, overlap_region_alloc* overlap_list,Correct_dumy* dumy, const ul_idx_t *uref, -char* r_string, double max_ov_diff_ec, long long blockLen, void *km) +int32_t init_waln(int64_t err, int64_t s, int64_t l, int64_t w_l, +int64_t* aux_beg, int64_t* aux_end, int64_t* r_s, int64_t* r_l) { - long long i, reverse_i = dumy->size - 1; - long long currentID, threshold = /**THRESHOLD_MAX_SIZE**/THRESHOLD; - long long x_start, y_start, o_len; - long long Window_Len = blockLen + (threshold << 1); - char* x_string = NULL; - char* y_string = NULL; - long long x_end, x_len; - int end_site; - unsigned int error; - int extra_begin; - int extra_end; + (*aux_beg) = (*aux_end) = (*r_s) = (*r_l) = -1; + ///since w_l == x_len + (err << 1) + if((s < 0) || (s >= l) || ((l-s+(2*err)+THRESHOLD_MAX_SIZE) < w_l)) return 0; + (*aux_beg) = (*aux_end) = 0; + ///s might be less than 0 + (*r_s) = s - err; + (*r_l) = l-(*r_s); if((*r_l) > w_l) (*r_l) = w_l; + (*aux_end) = w_l - (*r_l); - ///here are overlaps fully covered by blockLen - for (i = 0; i < (long long)dumy->length; i++) - { - extra_begin = extra_end = 0; - ///if the window has been fully covered, the interval at x is [window_start, window_end] - x_len = blockLen; - currentID = dumy->overlapID[i]; - x_start = window_start; - ///offset of y - y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s; - /****************************may have bugs********************************/ - y_start += y_start_offset(x_start, &(overlap_list->list[currentID].f_cigar)); - /****************************may have bugs********************************/ - - - if(!determine_overlap_region(threshold, y_start, overlap_list->list[currentID].y_id, Window_Len, uref->ug->u.a[overlap_list->list[currentID].y_id].len, - &extra_begin, &extra_end, &y_start, &o_len)) { - continue; - } - - fill_subregion_ul(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, - uref, overlap_list->list[currentID].y_id, extra_begin, extra_end); - - end_site = Reserve_Banded_BPM(dumy->overlap_region, Window_Len, x_string, blockLen, threshold, &error); - if (error!=(unsigned int)-1) { - overlap_list->list[currentID].align_length += x_len; - append_window_list(&overlap_list->list[currentID], window_start, window_end, y_start, y_start + end_site, (int)error, - extra_begin, extra_end, threshold, blockLen, km); - } - } - - ///here are overlaps partially covered by blockLen - for (i = 0; i < (long long)dumy->lengthNT; i++) - { - extra_begin = extra_end = 0; - currentID = dumy->overlapID[reverse_i--]; - x_start = MAX(window_start, (long long)overlap_list->list[currentID].x_pos_s); - x_end = MIN(window_end, (long long)overlap_list->list[currentID].x_pos_e); - - ///overlap length between [window_start, window_end] - x_len = x_end - x_start + 1; - threshold = x_len * max_ov_diff_ec; - /****************************may have bugs********************************/ - threshold = Adjust_Threshold(threshold, x_len); - if(threshold > THRESHOLD_MAX_SIZE) threshold = THRESHOLD_MAX_SIZE; - /****************************may have bugs********************************/ - - ///offset of y - y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s; - /****************************may have bugs********************************/ - y_start += y_start_offset(x_start, &(overlap_list->list[currentID].f_cigar)); - /****************************may have bugs********************************/ - - Window_Len = x_len + (threshold << 1); - - if(!determine_overlap_region(threshold, y_start, overlap_list->list[currentID].y_id, Window_Len, uref->ug->u.a[overlap_list->list[currentID].y_id].len, - &extra_begin, &extra_end, &y_start, &o_len)) { - continue; - } - - fill_subregion_ul(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, - uref, overlap_list->list[currentID].y_id, extra_begin, extra_end); - - x_string = r_string + x_start; - y_string = dumy->overlap_region; - - end_site = Reserve_Banded_BPM(y_string, Window_Len, x_string, x_len, threshold, &error); - - if (error!=(unsigned int)-1) { - overlap_list->list[currentID].align_length += x_len; - append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, y_start + end_site, (int)error, - extra_begin, extra_end, threshold, blockLen, km); - } + if ((*r_s) < 0) { + (*aux_beg) = -(*r_s); (*r_s) = 0; (*r_l) -= (*aux_beg); } + return 1; +} + +///[s, e) +int64_t get_num_wins(int64_t s, int64_t e, int64_t block_s) +{ + int64_t nl = e - ((s/block_s)*block_s), nw; + nw = (nl/block_s); if((nl%block_s)>0) nw++; + return nw; +} +///[w_s, w_e] +inline int64_t get_win_id_by_s(overlap_region *z, int64_t w_s, int64_t block_s, int64_t *w_e) +{ + int64_t n_s = ((z->x_pos_s/block_s)*block_s), wid = (w_s-n_s)/block_s; + if(w_e) { + (*w_e) = n_s + (wid+1)*block_s - 1; + if((*w_e) > z->x_pos_e) (*w_e) = z->x_pos_e; + } + return wid; +} + +///[w_s, w_e] +inline int64_t get_win_id_by_e(overlap_region *z, int64_t w_e, int64_t block_s, int64_t *w_s) +{ + int64_t n_s = ((z->x_pos_s/block_s)*block_s), wid = (w_e-n_s)/block_s; + if(w_s) { + (*w_s) = n_s + wid*block_s; + if((*w_s) < z->x_pos_s) (*w_s) = z->x_pos_s; + } + return wid; +} + +///[w_s, w_e] +inline void get_win_se_by_normalize_xs(overlap_region *z, int64_t norm_w_s, int64_t block_s, int64_t *w_s, int64_t *w_e) +{ + int64_t n_s = ((z->x_pos_s/block_s)*block_s), wid = (norm_w_s-n_s)/block_s; + if(w_s) { + (*w_s) = n_s + wid*block_s; + if((*w_s) < z->x_pos_s) (*w_s) = z->x_pos_s; + } + if(w_e) { + (*w_e) = n_s + (wid+1)*block_s - 1; + if((*w_e) > z->x_pos_e) (*w_e) = z->x_pos_e; + } +} + +void gen_str_seq(char *dst, int64_t s, int64_t pri_l, uint8_t rev, const ul_idx_t *uref, long long id, int64_t aux_beg, int64_t aux_end) +{ + // int64_t l = pri_l + aux_beg + aux_end; + memset(dst, 'N', aux_beg); + retrieve_u_seq(NULL, dst+aux_beg, &(uref->ug->u.a[id]), rev, s, pri_l, NULL); + memset(dst+aux_beg+pri_l, 'N', aux_end); +} + +void verify_ul_window_s(overlap_region *z, const ul_idx_t *uref, char* qstr, char *tstr, +double e_rate, int64_t w_l, int64_t e_max, void *km) +{ + int64_t q_s, q_e, nw, k, q_l; + int64_t aux_beg, aux_end, t_s, thre, aln_l, t_pri_l, t_end; + char *q_string, *t_string; unsigned int error; + z->w_list.n = 0; nw = get_num_wins(z->x_pos_s, z->x_pos_e+1, w_l); + get_win_se_by_normalize_xs(z, (z->x_pos_s/w_l)*w_l, w_l, &q_s, &q_e); + // q_s = z->x_pos_s; get_win_id_by_s(z, q_s, w_l, &q_e); + for (k = 0; k < nw; k++) { + aux_beg = aux_end = 0; q_l = 1 + q_e - q_s; + thre = q_l*e_rate; thre = Adjust_Threshold(thre, q_l); + if(thre > THRESHOLD_MAX_SIZE) thre = THRESHOLD_MAX_SIZE; + ///offset of y + t_s = (q_s - z->x_pos_s) + z->y_pos_s; + t_s += y_start_offset(q_s, &(z->f_cigar)); + + aln_l = q_l + (thre<<1); + // if(z->y_id == 115) { + // fprintf(stderr, "+[M::] q_s::%ld, t_s::%ld, t_pri_l::%ld, aux_beg::%ld, aux_end::%ld, aln_l::%ld, t_end::%ld, error::%u\n", + // q_s, t_s, t_pri_l, aux_beg, aux_end, aln_l, t_end, error); + // } + if(init_waln(thre, t_s, uref->ug->u.a[z->y_id].len, aln_l, &aux_beg, &aux_end, &t_s, &t_pri_l)) { + gen_str_seq(tstr, t_s, t_pri_l, z->y_pos_strand, uref, z->y_id, aux_beg, aux_end); + q_string = qstr+q_s; t_string = tstr; + t_end = Reserve_Banded_BPM(t_string, aln_l, q_string, q_l, thre, &error); + // if(z->y_id == 115) { + // fprintf(stderr, "-[M::] q_s::%ld, t_s::%ld, t_pri_l::%ld, aux_beg::%ld, aux_end::%ld, aln_l::%ld, t_end::%ld, error::%u, thre::%ld\n", + // q_s, t_s, t_pri_l, aux_beg, aux_end, aln_l, t_end, error, thre); + // } + if (error!=((unsigned int)-1)) { + z->align_length += q_l; + ///t_s do not have aux_beg, while t_s + t_end (aka, te) has + append_window_list(z, q_s, q_e, t_s, t_s + t_end, error, aux_beg, aux_end, thre, w_l, km); + } + } + q_s = q_e + 1; q_e = q_s + w_l - 1; + if(q_e >= (int64_t)z->x_pos_e) q_e = z->x_pos_e; + } + + // if(q_e != (int64_t)z->x_pos_e) { + // fprintf(stderr, "[M::%s] q_e::%ld, z->x_pos_s::%u, z->x_pos_e::%u, w_l::%ld, nw::%ld\n", __func__, + // q_e, z->x_pos_s, z->x_pos_e, w_l, nw); + // } + assert(q_e == (int64_t)z->x_pos_e); } ///error_rate should be 30% @@ -996,49 +1031,6 @@ unsigned int* get_error, int* get_y_end, int* get_x_end, int* get_aligned_xLen) } } -///[s, e) -int64_t get_num_wins(int64_t s, int64_t e, int64_t block_s) -{ - int64_t nl = e - ((s/block_s)*block_s), nw; - nw = (nl/block_s); if((nl%block_s)>0) nw++; - return nw; -} -///[w_s, w_e] -inline int64_t get_win_id_by_s(overlap_region *z, int64_t w_s, int64_t block_s, int64_t *w_e) -{ - int64_t n_s = ((z->x_pos_s/block_s)*block_s), wid = (w_s-n_s)/block_s; - if(w_e) { - (*w_e) = n_s + (wid+1)*block_s - 1; - if((*w_e) > z->x_pos_e) (*w_e) = z->x_pos_e; - } - return wid; -} - -///[w_s, w_e] -inline int64_t get_win_id_by_e(overlap_region *z, int64_t w_e, int64_t block_s, int64_t *w_s) -{ - int64_t n_s = ((z->x_pos_s/block_s)*block_s), wid = (w_e-n_s)/block_s; - if(w_s) { - (*w_s) = n_s + wid*block_s; - if((*w_s) < z->x_pos_s) (*w_s) = z->x_pos_s; - } - return wid; -} - -///[w_s, w_e] -inline void get_win_se_by_normalize_xs(overlap_region *z, int64_t norm_w_s, int64_t block_s, int64_t *w_s, int64_t *w_e) -{ - int64_t n_s = ((z->x_pos_s/block_s)*block_s), wid = (norm_w_s-n_s)/block_s; - if(w_s) { - (*w_s) = n_s + wid*block_s; - if((*w_s) < z->x_pos_s) (*w_s) = z->x_pos_s; - } - if(w_e) { - (*w_e) = n_s + (wid+1)*block_s - 1; - if((*w_e) > z->x_pos_e) (*w_e) = z->x_pos_e; - } -} - inline int64_t get_init_err_thres(int64_t len, double e_rate, int64_t block_s, int64_t block_err) { if(len >= block_s) return block_err; @@ -3590,6 +3582,10 @@ inline void refine_ed_aln(overlap_region_alloc* overlap_list, All_reads *rref, c overlap_list->mapped_overlaps_length = 0; on = overlap_list->length; for (j = 0; j < on; ++j) { z = &(overlap_list->list[j]); z->is_match = 0; is_srt = 1; + // /**if(z->y_id == 4)**/ { + // fprintf(stderr, "[M::%s::idx->%ld::] z->y_id::%u, z->w_list.n::%u, x_pos_s::%u, x_pos_e::%u, y_pos_s::%u, y_pos_e::%u\n", + // __func__, j, z->y_id, (uint32_t)z->w_list.n, z->x_pos_s, z->x_pos_e, z->y_pos_s, z->y_pos_e); + // } if(z->w_list.n == 0) continue;///no alignment nw = get_num_wins(z->x_pos_s, z->x_pos_e+1, block_s); a_nw = z->w_list.n; kv_resize(uint64_t, v_idx->a, (uint64_t)nw); @@ -3598,6 +3594,12 @@ inline void refine_ed_aln(overlap_region_alloc* overlap_list, All_reads *rref, c assert(z->w_list.a[i].y_end != -1); w_id = get_win_id_by_s(z, z->w_list.a[i].x_start, block_s, NULL); w_idx[w_id] = i; + // /**if(z->y_id == 4)**/ { + // fprintf(stderr, "[M::%s::i->%ld::] x_start::%d, x_end::%d, y_start::%d, y_end::%d, error::%d, error_threshold::%d, extra_begin::%d, extra_end::%d\n", + // __func__, i, z->w_list.a[i].x_start, z->w_list.a[i].x_end, + // z->w_list.a[i].y_start, z->w_list.a[i].y_end, z->w_list.a[i].error, + // z->w_list.a[i].error_threshold, z->w_list.a[i].extra_begin, z->w_list.a[i].extra_end); + // } } y_id = z->y_id; y_strand = z->y_pos_strand; @@ -9163,35 +9165,13 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref clear_Correct_dumy(dumy, overlap_list, km); - long long window_start, window_end; - Window_Pool w_inf; - - init_Window_Pool(&w_inf, g_read->length, /**WINDOW_UL**//**WINDOW_UL_H**/winLen, (int)(1.0/max_ov_diff_ec)); - int flag = 0; - - while(get_Window(&w_inf, &window_start, &window_end) && flag != -2) - { - dumy->length = 0; - dumy->lengthNT = 0; - flag = get_interval(window_start, window_end, overlap_list, dumy, w_inf.window_length); - - switch (flag) - { - case 1: ///no match here - break; - case 0: ///no match here - break; - case -2: ///if flag == -2, loop would be terminated - break; - } - - ///dumy->lengthNT represent how many overlaps that the length of them is not equal to WINDOW; may larger or less than WINDOW - ///dumy->length represent how many overlaps that the length of them is WINDOW - ///now the windows which are larger than WINDOW are verified one-by-one, to improve it, we can do it group-bygroup - verify_ul_window(window_start, window_end, overlap_list, dumy, uref, g_read->seq, max_ov_diff_ec, w_inf.window_length, /**THRESHOLD**/THRESHOLD_MAX_SIZE, km); - // verify_ul_ll_window(window_start, window_end, overlap_list, dumy, uref, g_read->seq, max_ov_diff_ec, w_inf.window_length, km); + init_Window_Pool(&w_inf, g_read->length, winLen, (int)(1.0/max_ov_diff_ec)); + uint64_t i; + for (i = 0; i < overlap_list->length; i++) { + verify_ul_window_s(&(overlap_list->list[i]), uref, g_read->seq, dumy->overlap_region, max_ov_diff_ec, w_inf.window_length, THRESHOLD_MAX_SIZE, km); } + // recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read); // partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat); @@ -9260,7 +9240,6 @@ void lchain_align(overlap_region_alloc* overlap_list, const ul_idx_t *uref, ///dumy->length represent how many overlaps that the length of them is WINDOW ///now the windows which are larger than WINDOW are verified one-by-one, to improve it, we can do it group-bygroup verify_ul_window(window_start, window_end, overlap_list, dumy, uref, g_read->seq, max_ov_diff_ec, w_inf.window_length, /**THRESHOLD**/THRESHOLD_MAX_SIZE, km); - // verify_ul_ll_window(window_start, window_end, overlap_list, dumy, uref, g_read->seq, max_ov_diff_ec, w_inf.window_length, km); } // recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read); diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 41156dc..7be0d7a 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -78,7 +78,54 @@ int get_fake_gap_shift(Fake_Cigar* x, int index) return result; } -int ovlp_chain_gen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int64_t yl, int64_t apend_be) +void gen_fake_cigar(Fake_Cigar* z, overlap_region *o, int64_t apend_be, k_mer_hit* hit, int64_t n_hit, int64_t print) +{ + int64_t k, dq, dr, dd, pdd; z->length = 0; + if(apend_be == 1) add_fake_cigar(z, o->x_pos_s, 0, NULL); + for (k = 0, pdd = INT32_MAX; k < n_hit; k++) { + dq = hit[k].self_offset - o->x_pos_s; + dr = hit[k].offset - o->y_pos_s; + dd = dr - dq; + if(print) { + fprintf(stderr, "[M::k->%lu] x::%u, y::%u, cnt::%u, dd::%ld, pdd::%ld, z->n::%u\n", + k, hit[k].self_offset, hit[k].offset, hit[k].cnt&(0xffu), dd, pdd, z->length); + } + if(dd != pdd) { + pdd = dd; + add_fake_cigar(z, hit[k].self_offset, pdd, NULL); + } + } + + if((apend_be == 1) && (get_fake_gap_pos(z, z->length-1)!=((int64_t)o->x_pos_e))) { + add_fake_cigar(z, o->x_pos_e, get_fake_gap_shift(z, z->length-1), NULL); + } +} + +void debug_cigar(Fake_Cigar* z, overlap_region *o, int64_t apend_be, k_mer_hit* hit, uint64_t n_hit) +{ + gen_fake_cigar(z, o, apend_be, hit, n_hit, 0); + if(!((z->length==o->f_cigar.length) && + (!memcmp(z->buffer, o->f_cigar.buffer, sizeof((*(o->f_cigar.buffer)))*o->f_cigar.length)))) { + uint64_t k; + fprintf(stderr, "\n[M::%s] z->n::%u, o->n::%u, rev::%u\n", __func__, z->length, o->f_cigar.length, o->y_pos_strand); + for (k = 0; k < z->length; k++) { + fprintf(stderr, "[z::k->%lu] pos::%d, off::%d\n", k, + get_fake_gap_pos(z, k), get_fake_gap_shift(z, k)); + } + for (k = 0; k < o->f_cigar.length; k++) { + fprintf(stderr, "[o::k->%lu] pos::%d, off::%d\n", k, + get_fake_gap_pos(&(o->f_cigar), k), get_fake_gap_shift(&(o->f_cigar), k)); + } + gen_fake_cigar(z, o, apend_be, hit, n_hit, 1); + // for (k = 0; k < o->f_cigar.length; k++) { + // fprintf(stderr, "[M::k->%lu] x::%u, y::%u, cnt::%u\n", + // k, hit[k].self_offset, hit[k].offset, hit[k].cnt&(0xffu)); + // } + + } +} + +int ovlp_chain_gen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int64_t yl, int64_t apend_be, k_mer_hit* hit, int64_t n_hit) { if (ol->length + 1 > ol->size) { uint64_t sl = ol->size; @@ -94,8 +141,7 @@ int ovlp_chain_gen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int6 ((ol->list[ol->length-1].shared_seed == t->shared_seed) && (ol->list[ol->length-1].overlapLen <= t->overlapLen))) { return 0; - } - else { + } else { ol->length--; } } @@ -126,7 +172,9 @@ int ovlp_chain_gen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int6 o->y_pos_strand = t->x_pos_strand; resize_fake_cigar(&(o->f_cigar), (t->f_cigar.length + 2), NULL); - if(apend_be == 1) add_fake_cigar(&(o->f_cigar), o->x_pos_s, 0, NULL); + if(apend_be == 1) { + add_fake_cigar(&(o->f_cigar), ((t->x_pos_strand)?(xl-t->x_pos_e-1):(t->x_pos_s)), 0, NULL); + } if (t->x_pos_strand == 1) { o->x_pos_e = xl-t->x_pos_s-1; o->x_pos_s = xl-t->x_pos_e-1; @@ -146,7 +194,7 @@ int ovlp_chain_gen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int6 dq = t->x_pos_e - t->x_pos_s; dr = t->y_pos_e - t->y_pos_s; - id = dr - dq; + id = dr - dq;///indel size from left pdd = INT32_MAX; fn = t->f_cigar.length; for (i = fn-1; i >= 0; i--) { @@ -159,10 +207,12 @@ int ovlp_chain_gen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int6 } if((apend_be == 1) && (get_fake_gap_pos(&(o->f_cigar), o->f_cigar.length-1) != ((int64_t)o->x_pos_e))) { - add_fake_cigar(&(o->f_cigar), o->x_pos_e, - get_fake_gap_shift(&(o->f_cigar), o->f_cigar.length-1), NULL); + add_fake_cigar(&(o->f_cigar), o->x_pos_e, get_fake_gap_shift(&(o->f_cigar), o->f_cigar.length-1), NULL); } - + + ///debug + debug_cigar(&(t->f_cigar), o, apend_be, hit, n_hit); + return 1; } @@ -1407,5 +1457,22 @@ uint64_t lchain_dp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, ov } res->overlapLen = get_chainLen(a[msc_i].self_offset, a[i].self_offset, xl, a[msc_i].offset, a[i].offset, yl); for (i = 0; i < cL; i++) des[i] = a[t[cL-i-1]]; + if(res->x_pos_strand) { + int64_t hcl = cL>>1; k_mer_hit kp; + for (i = 0; i < hcl; i++) { + j = cL-i-1; kp = des[i]; des[i] = des[j]; des[j] = kp; + + des[i].self_offset = xl-des[i].self_offset-1; + des[i].offset = yl-des[i].offset-1; + + des[j].self_offset = xl-des[j].self_offset-1; + des[j].offset = yl-des[j].offset-1; + } + if(cL&1) { + des[i].self_offset = xl-des[i].self_offset-1; + des[i].offset = yl-des[i].offset-1; + } + } + return cL; } \ No newline at end of file diff --git a/Hash_Table.h b/Hash_Table.h index c780b24..e96a4ab 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -190,7 +190,7 @@ long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* uint64_t lchain_dp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, overlap_region* res, int64_t max_skip, int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, int64_t xl, int64_t yl, int64_t quick_check); -int ovlp_chain_gen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int64_t yl, int64_t apend_be); +int ovlp_chain_gen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int64_t yl, int64_t apend_be, k_mer_hit* hit, int64_t n_hit); int append_utg_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, ma_utg_v *ua, int add_beg_end, void *km); diff --git a/Overlaps.h b/Overlaps.h index d0cff42..5a0ded8 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -249,8 +249,15 @@ typedef struct { uint64_t *idx; } cvert_t; +typedef struct { + ma_ug_t *hg; + size_t n, m; + uint64_t *a; +} hpc_t; + typedef struct { ma_ug_t *ug; + hpc_t *hpc_g; ucov_t *cc; ucov_t *cr; ul_contain *ct; diff --git a/anchor.cpp b/anchor.cpp index 75e11ed..27c035b 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -784,10 +784,18 @@ void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list) void minimizers_gen(ha_abufl_t *ab, char* rs, int64_t rl, uint64_t mz_w, uint64_t mz_k, Candidates_list *cl, kvec_t_u8_warp* k_flag, -void *ha_flt_tab, ha_pt_t *ha_idx, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t high_occ) +void *ha_flt_tab, ha_pt_t *ha_idx, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ) { - uint64_t i, k, l; int n, j; ha_mzl_t *z; seedl_t *s; - if(high_occ < 1) high_occ = 1; + // fprintf(stderr, "+[M::%s]\n", __func__); + uint64_t i, k, l, max_cnt = UINT32_MAX, min_cnt = 0; int n, j; ha_mzl_t *z; seedl_t *s; + if(high_occ) { + max_cnt = (*high_occ); + if(max_cnt < 2) max_cnt = 2; + } + if(low_occ) { + min_cnt = (*low_occ); + if(min_cnt < 2) min_cnt = 2; + } clear_Candidates_list(cl); ab->mz.n = 0, ab->n_a = 0; // get the list of anchors @@ -847,11 +855,12 @@ void *ha_flt_tab, ha_pt_t *ha_idx, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_ p->strand = ab->a[k].srt >> 32 & 1; p->offset = ab->a[k].other_off; p->self_offset = ab->a[k].self_off; - if((ab->a[k].cnt>>8) <= high_occ){ + if(((ab->a[k].cnt>>8) < max_cnt) && ((ab->a[k].cnt>>8) > min_cnt)){ p->cnt = 1; - } - else{ - p->cnt = 1 + (((ab->a[k].cnt>>8) + (high_occ<<1) - 1)/(high_occ<<1)); + } else if((ab->a[k].cnt>>8) <= min_cnt) { + p->cnt = 2; + } else{ + p->cnt = 1 + (((ab->a[k].cnt>>8) + (max_cnt<<1) - 1)/(max_cnt<<1)); p->cnt = pow(p->cnt, 1.1); } if(p->cnt > ((uint32_t)(0xffffffu))) p->cnt = 0xffffffu; @@ -860,10 +869,11 @@ void *ha_flt_tab, ha_pt_t *ha_idx, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_ cl->length = ab->n_a; } -void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint64_t rid, uint64_t rl, All_reads* rdb, - const ul_idx_t *udb, uint32_t beg_tail, 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) +void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uint64_t rl, All_reads* rdb, + const ul_idx_t *udb, uint32_t apend_be, overlap_region* tf, uint64_t max_n_chain, + int64_t max_skip, int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, int64_t quick_check, uint32_t gen_off) { + // fprintf(stderr, "+[M::%s]\n", __func__); uint64_t i, k, l, m, sm, cn = cl->length; clear_overlap_region_alloc(ol); clear_fake_cigar(&(tf->f_cigar)); @@ -877,10 +887,18 @@ void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint64_t rid, uin tf->x_pos_strand = cl->list[l].strand; tf->y_id = cl->list[l].readID; tf->y_pos_strand = 0;///always 0 + // fprintf(stderr, "+[M::%s] l::%lu, k::%lu\n", __func__, l, k); sm = lchain_dp(cl->list+l, k-l, cl->list+m, &(cl->chainDP), tf, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_rate, rl, rdb?Get_READ_LENGTH((*rdb), (*tf).y_id):udb->ug->u.a[(*tf).y_id].len, quick_check); - assert(sm > 0); - if(ovlp_chain_gen(ol, tf, rl, rdb?Get_READ_LENGTH((*rdb), (*tf).y_id):udb->ug->u.a[(*tf).y_id].len, beg_tail)) { + // assert(sm > 0); + // uint64_t z; + // for (z = 0; z < sm; z++) { + // fprintf(stderr, "+[M::z->%lu] x::%u, y::%u, cnt::%u\n", + // z, cl->list[m+z].self_offset, cl->list[m+z].offset, cl->list[m+z].cnt&(0xffu)); + // } + + // fprintf(stderr, "+[M::%s] l::%lu, k::%lu, sm::%lu, # cigar::%u, yid::%u, sc::%d\n", __func__, l, k, sm, tf->f_cigar.length, tf->y_id, tf->shared_seed); + if(ovlp_chain_gen(ol, tf, rl, rdb?Get_READ_LENGTH((*rdb), (*tf).y_id):udb->ug->u.a[(*tf).y_id].len, apend_be, cl->list+m, sm)) { m += sm; } } @@ -929,7 +947,7 @@ void set_lchain_dp_op(uint32_t is_accurate, uint32_t mz_k, int64_t *max_skip, in 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) = 0; (*max_skip) = 25; (*max_iter) = 5000; (*max_dis) = 5000; div = 0.1; pen_gap = 0.5f; pen_skip = 0.0005f; + (*quick_check) = 1; (*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 @@ -938,14 +956,14 @@ void set_lchain_dp_op(uint32_t is_accurate, uint32_t mz_k, int64_t *max_skip, in *chn_pen_skip = pen_skip * tmp; } -void ul_map_lchain(ha_abufl_t *ab, int64_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, - int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t high_occ, uint32_t is_accurate) +void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, + int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off) { extern void *ha_flt_tab; extern ha_pt_t *ha_idx; int64_t max_skip, max_iter, max_dis, quick_check; double chn_pen_gap, chn_pen_skip; set_lchain_dp_op(is_accurate, mz_k, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check); - minimizers_gen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, dbg_ct, sp, high_occ); - lchain_gen(cl, overlap_list, rid, rl, NULL, uref, keep_whole_chain, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check); + minimizers_gen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, dbg_ct, sp, high_occ, low_occ); + lchain_gen(cl, overlap_list, rid, rl, NULL, uref, apend_be, f_cigar, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off); ///no need to sort here, overlap_list has been sorted at lchain_gen } diff --git a/gchain_map.cpp b/gchain_map.cpp index 09f142b..eb4b25f 100644 --- a/gchain_map.cpp +++ b/gchain_map.cpp @@ -19,7 +19,7 @@ #include "gchain_map.h" KSEQ_INIT(gzFile, gzread) void ul_map_lchain(ha_abufl_t *ab, int64_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, - int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t high_occ, uint32_t is_accurate); + int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off); typedef struct { // global data structure for kt_pipeline() const void *ha_flt_tab; diff --git a/inter.cpp b/inter.cpp index 7effda3..8fff5b2 100644 --- a/inter.cpp +++ b/inter.cpp @@ -20,7 +20,8 @@ KSEQ_INIT(gzFile, gzread) void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t high_occ, void *km); - +void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, + int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off); #define MG_SEED_IGNORE (1ULL<<41) #define MG_SEED_TANDEM (1ULL<<42) #define MG_SEED_KEPT (1ULL<<43) @@ -5215,6 +5216,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call ha_ovec_buf_t *b = s->hab[tid]; glchain_t *bl = &(s->ll[tid]); int64_t /**rid = s->id+i,**/ winLen = MIN((((double)THRESHOLD_MAX_SIZE)/s->opt->diff_ec_ul), WINDOW); + uint32_t high_occ = 2; // uint64_t align = 0; int fully_cov, abnormal; // if(UL_INF.a[s->id+i].rlen != s->len[i]) { @@ -5225,15 +5227,17 @@ 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!=97) return; + // if(s->id+i!=47) 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); // if (memcmp(UL_INF.nid.a[s->id+i].a, "d0aab024-b3a7-40fb-83cc-22c3d6d951f8", UL_INF.nid.a[s->id+i].n-1)) return; // fprintf(stderr, "[M::%s::] ==> len: %lu\n", __func__, s->len[i]); - ha_get_ul_candidates_interface(b->abl, i, s->seq[i], s->len[i], s->opt->w, s->opt->k, s->uu, &b->olist, &b->olist_hp, &b->clist, s->opt->bw_thres, - s->opt->max_n_chain, 1, NULL/**&(b->k_flag)**/, &b->r_buf, &(b->tmp_region), NULL, &(b->sp), 1, NULL); - + // ha_get_ul_candidates_interface(b->abl, i, s->seq[i], s->len[i], s->opt->w, s->opt->k, s->uu, &b->olist, &b->olist_hp, &b->clist, s->opt->bw_thres, + // s->opt->max_n_chain, 1, NULL, &b->r_buf, &(b->tmp_region), NULL, &(b->sp), 1, NULL); + ul_map_lchain(b->abl, (uint32_t)-1, s->seq[i], s->len[i], s->opt->w, s->opt->k, s->uu, &b->olist, &b->clist, s->opt->bw_thres, + s->opt->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 0); + clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2); // return; @@ -9898,6 +9902,121 @@ ul_contain *ul_contain_gen(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc* src, int64_t return p; } +void gen_hpc_seq(const char *in, uint32_t in_len, ma_utg_t *ou) +{ + uint32_t k, l, m; memset(ou, 0, sizeof((*ou))); + for (l = 0, k = 1; k <= in_len; k++) { + if((k == in_len) || (in[k] != in[l]) || (seq_nt4_table[(uint8_t)in[l]] >= 4)) { + ou->len++; l = k; + } + } + + MALLOC(ou->s, ou->len); m = 0; + for (l = 0, k = 1; k <= in_len; k++) { + if((k == in_len) || (in[k] != in[l]) || (seq_nt4_table[(uint8_t)in[l]] >= 4)) { + ou->s[m++] = in[l]; + l = k; + } + } +} + +void gen_microsatellite(const char *in, uint32_t in_len, ma_utg_t *idx, hpc_t *res, uint32_t mcs_len) +{ + uint32_t k, i, l, c, o; char sk[256]; + idx->start = idx->end = res->n; + if(mcs_len > 256) mcs_len = 256; + for (k = 1; k <= mcs_len; k++) { + memset(sk, 'N', k); ///k->length of k-mer + o = k + (k>>1) + 1; + for (i = l = 0; i < in_len; i++) { + c = seq_nt4_table[(uint8_t)in[i]]; + if((c < 4) && (((l >= k) && (sk[l%k] == in[i])) || (l < k))) { + if(l < k) sk[l] = in[i]; + l++; + } else { + if(l >= o) { + kv_push(uint64_t, *res, (((uint64_t)(i-l))<<1)); + kv_push(uint64_t, *res, ((((uint64_t)i)<<1)|1)); + } + l = 0; + } + } + } + + radix_sort_gfa64(res->a+idx->start, res->a+res->n); + int64_t dp, old_dp; uint64_t st, en; + for (k = st = idx->start, dp = 0; k < res->n; k++) { + old_dp = dp; + if (res->a[k]&1) --dp;//qe + else ++dp; + if (old_dp < 1 && dp >= 1) {///qs + st = res->a[k]>>1; + } else if (old_dp >= 1 && dp < 1) {///qe + en = res->a[k]>>1; + res->a[idx->end++] = ((st<<32)|(en)); + } + } + res->n = idx->end; +} + +uint32_t hpc_l(char *s, int64_t hof, int64_t sof, int64_t scut) +{ + if(scut == sof) return hof; + if(scut < sof) sof = hof = 0; + int64_t l, k; + for (k=sof; (k>0)&&(k= 4)) { + hof++; l = k; + } + } + return hof; +} + +hpc_t *hpc_g_gen(ma_ug_t *ug) +{ + uint32_t k, i, ho, so, len; int32_t z; hpc_t *p; kvec_t(char) cc; asg_t *ng = asg_init(); + CALLOC(p, 1); CALLOC(p->hg, 1); kv_init(cc); + CALLOC(p->hg->u.a, ug->u.n); p->hg->u.n = p->hg->u.m = ug->u.n; + + for (k = 0; k < ug->u.n; k++) { + kv_resize(char, cc, ug->u.a[k].len); + retrieve_u_seq(NULL, cc.a, &(ug->u.a[k]), 0, 0, ug->u.a[k].len, NULL); + gen_hpc_seq(cc.a, ug->u.a[k].len, &(p->hg->u.a[k])); + gen_microsatellite(cc.a, ug->u.a[k].len, &(p->hg->u.a[k]), p, 6); + } + + ng->m_arc = ng->n_arc = ug->g->n_arc; CALLOC(ng->arc, ng->n_arc); + ng->m_seq = ng->n_seq = ug->g->n_seq; CALLOC(ng->seq, ng->n_seq); + for (k = 0; k < ng->n_seq; k++) { + ng->seq[k].del = ng->seq[k].c = 0; ng->seq[k].len = p->hg->u.a[k].len;///hpc len + } + memcpy(ng->arc, ug->g->arc, ng->n_arc*(sizeof((*(ng->arc))))); + for (k = 1, i = 0; k <= ng->n_arc; k++) { + if(k == ng->n_arc || (ng->arc[i].ul>>32) != (ng->arc[k].ul>>32)) { + //sorted by ol + len = ng->arc[i].ol; kv_resize(char, cc, len); + retrieve_u_seq(NULL, cc.a, &(ug->u.a[ng->arc[i].ul>>33]), ((ng->arc[i].ul>>32)&1)^1, 0, len, NULL); + // min_o = ng->arc[k-1].ol; max_o = ng->arc[i].ol; + len = ng->seq[ng->arc[i].ul>>33].len;///hpc len + for (z = k-1, ho = so = 0; z >= (int32_t)i; z--) { + ho = hpc_l(cc.a, ho, so, ng->arc[z].ol); + so = ng->arc[z].ol; + ng->arc[z].ol = ho; + ng->arc[z].ul>>=32; + ng->arc[z].ul<<=32; + ng->arc[z].ul += len - ng->arc[z].ol; + assert(z == (int32_t)k-1 || ng->arc[z].ol>=ng->arc[z+1].ol); + } + i = k; + } + } + asg_cleanup(ng); + p->hg->g = ng; + kv_destroy(cc); + return p; +} void debug_append_inexact_edges(ma_ug_t *ug, const ug_opt_t *uopt) { uint32_t n_asymm = 0, n_disconnect = 0, z, v, w, k, nv; asg_arc_t *av = NULL; @@ -10238,6 +10357,8 @@ ul_idx_t *gen_ul_idx(const ug_opt_t *uopt, ma_ug_t *ug, asg_t *sg) uu->ct = ul_contain_gen(ug, sg, src, min_ovlp, max_hang, 0, 1); uu->cc = gen_cov_track(ug, sg, uu->ct, src, min_ovlp, max_hang, 0, 1); uu->cr = gen_r_contain(ug, sg, src, R_INF.total_reads, min_ovlp, max_hang, asm_opt.thread_num, 0, 1); + uu->hpc_g = hpc_g_gen(ug); + return uu; } @@ -10300,35 +10421,39 @@ void destroy_ul_idx_t(ul_idx_t *uu) { if(!uu) return; if(uu->cc) { - if(uu->cc) { - free(uu->cc->idx); - free(uu->cc->interval.a); - free(uu->cc); - } + free(uu->cc->idx); + free(uu->cc->interval.a); + free(uu->cc); + } - if(uu->cr) { - free(uu->cr->idx); - free(uu->cr->interval.a); - free(uu->cr); - } + if(uu->cr) { + free(uu->cr->idx); + free(uu->cr->interval.a); + free(uu->cr); + } - if(uu->ct) { - free(uu->ct->idx.a); - free(uu->ct->rids.a); - free(uu->ct->is_c.a); - free(uu->ct); - } + if(uu->ct) { + free(uu->ct->idx.a); + free(uu->ct->rids.a); + free(uu->ct->is_c.a); + free(uu->ct); + } - if(uu->r_ug) { - free(uu->r_ug->idx); - free(uu->r_ug->p.a); - free(uu->r_ug); - } + if(uu->r_ug) { + free(uu->r_ug->idx); + free(uu->r_ug->p.a); + free(uu->r_ug); + } + + if(uu->hpc_g) { + free(uu->hpc_g->a); + ma_ug_destroy(uu->hpc_g->hg); + free(uu->hpc_g); + } // if(uu->ov) { // free(uu->ov->a); // free(uu->ov); // } - } ma_ug_destroy(uu->ug); // if(uu->nug) { // free(uu->nug->idx);