From ee51a6dddc8b9886f3d448c91af3f85728dafb4c Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Sun, 4 Sep 2022 00:18:51 -0400 Subject: [PATCH] seperate alignment --- Correct.cpp | 702 +++++++++++++++++++++++++++++++++++++++++++++++-- Hash_Table.cpp | 173 +++++++++++- Hash_Table.h | 1 + Overlaps.h | 10 + anchor.cpp | 68 ++++- inter.cpp | 32 ++- 6 files changed, 925 insertions(+), 61 deletions(-) diff --git a/Correct.cpp b/Correct.cpp index f09da12..6530a79 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -722,8 +722,7 @@ char* r_string, double max_ov_diff_ec, long long blockLen, long long max_error, } } -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) +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) { (*aux_beg) = (*aux_end) = (*r_s) = (*r_l) = -1; ///since w_l == x_len + (err << 1) @@ -1081,6 +1080,7 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* rref, const ul_i tLen += z->w_list.a[k].x_end + 1 - z->w_list.a[k].x_start; tErr += z->w_list.a[k].error;///matched window // if(k != w_id) z->w_list.a[w_id] = z->w_list.a[k]; + ///from mapped window w_list.a[k] to the following unmapped windows for (m = w_id+1, w_e = z->w_list.a[k].x_end; m < idx_e; m++) { w_s = w_e + 1; wn_id = get_win_id_by_s(z, w_s, block_s, &w_e); @@ -1126,7 +1126,7 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* rref, const ul_i if(y_beg_left == -1 && y_beg_right != -1) y_beg_left = y_beg_right; if(y_beg_right == -1 && y_beg_left != -1) y_beg_right = y_beg_left; - if(y_beg_left != -1) { + if(y_beg_left != -1) {///note: this function will change tstr/qstr if(rref) { verify_sub_window(rref, dumy, g_read, w_s, x_len, y_beg_left, Window_Len, z->y_id, z->y_pos_strand, err_thre, 0, &r_error_left, &r_y_end_left, &r_x_end_left, &aligned_xLen_left); @@ -2307,6 +2307,96 @@ int* r_extra_begin, int* r_extra_end, unsigned int* r_error) return 0; } +inline char *return_str_seq(char *buf, int64_t s, int64_t pri_l, uint8_t rev, hpc_t *hpc_g, const ul_idx_t *uref, int64_t id, int64_t aux_beg, int64_t aux_end) +{ + if(!hpc_g) { + memset(buf, 'N', aux_beg); + retrieve_u_seq(NULL, buf+aux_beg, &(uref->ug->u.a[id]), rev, s, pri_l, NULL); + memset(buf+aux_beg+pri_l, 'N', aux_end); + return buf; + } else { + char *z = hpc_str(*hpc_g, id, rev); + if((aux_beg == 0) && (aux_end == 0)) { + return z+s; + } else { + memset(buf, 'N', aux_beg); + memcpy(buf+aux_beg, z+s, pri_l); + memset(buf+aux_beg+pri_l, 'N', aux_end); + return buf; + } + } +} + +///cannot use tstr in-place +inline int recal_boundary(char* qstr, char* tstr1, int64_t ql, int64_t thres, +int64_t global_ts0, int64_t local_ts0, int64_t local_te0, +int64_t aux_beg0, int64_t aux_end0, unsigned int err0, +int64_t tid, int64_t aln_l, uint32_t rev, +Correct_dumy* dumy, All_reads* rref, hpc_t *hpc_g, const ul_idx_t *uref, +int64_t* global_ts1, int* local_ts1, int* local_te1, +int64_t* aux_beg1, int64_t* aux_end1, unsigned int* err1) +{ + int64_t ts, t_tot_l, aux_beg, aux_end, t_pri_l, t_end; + char *q_string = qstr, *t_string; unsigned int error = (unsigned int)-1; + int r_ts = 0, path_length = 0; + if(hpc_g) t_tot_l = hpc_len(*hpc_g, tid); + else if(uref) t_tot_l = uref->ug->u.a[tid].len; + else t_tot_l = Get_READ_LENGTH((*rref), tid); + + if(local_ts0 == 0) {//left boundary + if(aux_beg0 > 0) return 0;///shift to the left cannot get a new start pos + ts = global_ts0; + } else if((local_te0 + 1) == aln_l) {//right boundary + if(aux_end0 > 0) return 0;///shift to the right cannot get a new start pos + ts = global_ts0 + local_te0 - ql + 1; + } else { + return 0; + } + if(!init_waln(thres, ts, t_tot_l, aln_l, &aux_beg, &aux_end, &ts, &t_pri_l)) return 0; + if(ts == global_ts0) return 0;//unchanged, make no sense + + if(rref) { + fill_subregion(tstr1, ts, t_pri_l, rev, rref, tid, aux_beg, aux_end); t_string = tstr1; + } else { + t_string = return_str_seq(tstr1, ts, t_pri_l, rev, hpc_g, uref, tid, aux_beg, aux_end); + } + + t_end = Reserve_Banded_BPM_PATH(t_string, aln_l, q_string, ql, thres, &error, &r_ts, + &path_length, dumy->matrix_bit, dumy->path_fix, -1, -1); + + if (error != (unsigned int)-1 && error < err0) { + (*global_ts1) = ts; + (*local_ts1) = r_ts; + (*local_te1) = t_end; + (*aux_beg1) = aux_beg; + (*aux_end1) = aux_end; + (*err1) = error; + + dumy->path_length = path_length; + memcpy(dumy->path, dumy->path_fix, path_length); + // memcpy(tstr0, t_string, aln_l); + return 1; + } + return 0; +} + +inline char *update_des_str(char *des, int64_t s, int64_t pri_l, uint8_t rev, All_reads *rref, hpc_t *hpc_g, + const ul_idx_t *uref, int64_t id, int64_t aux_beg, int64_t aux_end, char *src) +{ + if(src) { + // memcpy(des, src, (pri_l+aux_beg+aux_end)); + // return des; + return src; + } else { + if(rref) { + fill_subregion(des, s, pri_l, rev, rref, id, aux_beg, aux_end); + return des; + } else { + return return_str_seq(des, s, pri_l, rev, hpc_g, uref, id, aux_beg, aux_end); + } + } +} + /** void debug_scan_cigar(overlap_region* sub_list) { @@ -3415,20 +3505,27 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea ///debug_window_cigar(overlap_list, g_read, dumy, rref, 1, 1); } -uint32_t inline simi_pass(int64_t ol, int64_t aln_ol, All_reads *rref, const ul_idx_t *uref, double *e_rate) +uint32_t inline simi_pass(int64_t ol, int64_t aln_ol, uint32_t second_ck, double *e_rate) { if(aln_ol == 0 || ol == 0) return 0; - - if(rref) { + if((!second_ck) && (!e_rate)) { if((ol*OVERLAP_THRESHOLD_FILTER) <= aln_ol) return 1; - } else if(uref) { - if(e_rate) { - if((ol*((double)(((double)1.0)-(*e_rate)))) <= aln_ol) return 1; - } else { - if(((ol*MIN_UL_ALIN_RATE) <= aln_ol) && (aln_ol >= MIN_UL_ALIN_LEN)) return 1; - } + } else if(e_rate) { + if((ol*((double)(((double)1.0)-(*e_rate)))) <= aln_ol) return 1; + } else if(second_ck) { + if(((ol*MIN_UL_ALIN_RATE) <= aln_ol) && (aln_ol >= MIN_UL_ALIN_LEN)) return 1; } + // if(rref) { + // if((ol*OVERLAP_THRESHOLD_FILTER) <= aln_ol) return 1; + // } else if(uref) { + // if(e_rate) { + // if((ol*((double)(((double)1.0)-(*e_rate)))) <= aln_ol) return 1; + // } else { + // if(((ol*MIN_UL_ALIN_RATE) <= aln_ol) && (aln_ol >= MIN_UL_ALIN_LEN)) return 1; + // } + // } + return 0; } @@ -3493,6 +3590,131 @@ int32_t y_strand, int32_t y_id) return 0; } +inline uint32_t gen_backtrace_adv(window_list *p, overlap_region *z, All_reads *rref, hpc_t *hpc_g, const ul_idx_t *uref, +char *qstr, char *tstr, char *tstr1, Correct_dumy* dumy, uint32_t rev, uint32_t id) +{ + int64_t qs, qe, ql, aln_l, t_pri_l, thres, ts; + int r_ts = 0, t_end; int64_t aux_beg, aux_end; + char *q_string, *t_string; unsigned int error; + ///there is no problem for x + qs = p->x_start; qe = p->x_end; ql = qe + 1 - qs; + thres = p->error_threshold; aln_l = ql + (thres<<1); + + ///y_start is the real y_start + ///for the window with cigar, y_start has already reduced extra_begin + ts = p->y_start; aux_beg = p->extra_begin; aux_end = p->extra_end; + t_pri_l = aln_l - aux_beg - aux_end; + + q_string = qstr + qs; + if(rref) { + fill_subregion(tstr, ts, t_pri_l, rev, rref, id, aux_beg, aux_end); t_string = tstr; + } else { + t_string = return_str_seq(tstr, ts, t_pri_l, rev, hpc_g, uref, id, aux_beg, aux_end); + } + + t_end = Reserve_Banded_BPM_PATH(t_string, aln_l, q_string, ql, thres, &error, &r_ts, + &(dumy->path_length), dumy->matrix_bit, dumy->path, p->error, p->y_end - ts); + + + // assert(error != (unsigned int)-1); + if(error != (unsigned int)-1) { + ///this condition is always wrong + ///in best case, r_ts = threshold, t_end = aln_l - thres - 1 + if (((t_end+1) == aln_l) || (r_ts == 0)) { + if(recal_boundary(q_string, tstr1, ql, thres, ts, r_ts, t_end, + aux_beg, aux_end, error, id, aln_l, rev, dumy, rref, hpc_g, uref, + &ts, &r_ts, &t_end, &aux_beg, &aux_end, &error)) { + p->error = error; p->extra_begin = aux_beg; p->extra_end = aux_end; + t_string = update_des_str(tstr, ts, aln_l-aux_beg-aux_end, rev, rref, hpc_g, uref, + id, aux_beg, aux_end, hpc_g?NULL:tstr1); + } + } + + generate_cigar(dumy->path, dumy->path_length, p, &(z->w_list), &r_ts, &t_end, &error, q_string, ql, t_string); + + p->y_start = ts + r_ts - aux_beg; + p->y_end = ts + t_end - aux_beg; + p->error = error; + return 1; + } + p->error = -1; + return 0; +} + +inline uint32_t aln_wlst_adv(overlap_region *z, All_reads *rref, hpc_t *hpc_g, +const ul_idx_t *uref, char *qstr, char *tstr, char *tstr1, Correct_dumy* dumy, +uint32_t rev, uint32_t id, int64_t qs, int64_t qe, int64_t t_s, int64_t block_s, +double e_rate, uint32_t is_cigar) +{ + int64_t ql, aln_l, t_tot_l; window_list *p = NULL; int r_ts = 0, t_end; + int64_t aux_beg, aux_end, t_pri_l; + int64_t thres; char *q_string, *t_string; unsigned int error; + ql = qe + 1 - qs; + ///there are two potiential reasons for unmatched window: + ///1. this window has a large number of differences + ///2. DP does not start from the right offset + if(rref) { + thres = double_error_threshold(get_init_err_thres(ql, e_rate, block_s, THRESHOLD), ql); + } else { + thres = double_ul_error_threshold(get_init_err_thres(ql, e_rate, block_s, THRESHOLD_MAX_SIZE), ql); + } + aln_l = ql + (thres << 1); + if(hpc_g) t_tot_l = hpc_len(*hpc_g, id); + else if(uref) t_tot_l = uref->ug->u.a[id].len; + else t_tot_l = Get_READ_LENGTH((*rref), id); + + if(!init_waln(thres, t_s, t_tot_l, aln_l, &aux_beg, &aux_end, &t_s, &t_pri_l)) return 0; + if(t_pri_l + thres < ql) return 0; + + q_string = qstr + qs; + if(rref) { + fill_subregion(tstr, t_s, t_pri_l, rev, rref, id, aux_beg, aux_end); t_string = tstr; + } else { + t_string = return_str_seq(tstr, t_s, t_pri_l, rev, hpc_g, uref, id, aux_beg, aux_end); + } + + if(is_cigar) { + ///note!!! need notification + t_end = Reserve_Banded_BPM_PATH(t_string, aln_l, q_string, ql, thres, &error, &r_ts, + &(dumy->path_length), dumy->matrix_bit, dumy->path, -1, -1); + } else { + ///note!!! need notification + t_end = Reserve_Banded_BPM(t_string, aln_l, q_string, ql, thres, &error); + } + if(error!=(unsigned int)-1) { + if(is_cigar) { + ///this condition is always wrong + ///in best case, r_ts = threshold, t_end = aln_l - thres - 1 + if (((t_end+1) == aln_l) || (r_ts == 0)) { + if(recal_boundary(q_string, tstr1, ql, thres, t_s, r_ts, t_end, + aux_beg, aux_end, error, id, aln_l, rev, dumy, rref, hpc_g, uref, + &t_s, &r_ts, &t_end, &aux_beg, &aux_end, &error)) { + t_string = update_des_str(tstr, t_s, aln_l-aux_beg-aux_end, rev, rref, hpc_g, uref, + id, aux_beg, aux_end, hpc_g?NULL:tstr1); + } + } + } + + kv_pushp(window_list, z->w_list, &p); + p->x_start = qs; p->x_end = qe; ///must set x_start/x_end here + if(is_cigar) { + generate_cigar(dumy->path, dumy->path_length, p, &(z->w_list), &r_ts, &t_end, &error, q_string, ql, t_string); + } else { + p->cidx = p->clen = 0; + } + p->y_start = t_s + r_ts;///difference + p->y_end = t_s + t_end; + p->error = error; + p->extra_begin = aux_beg; + p->extra_end = aux_end; + p->error_threshold = thres; + z->align_length += ql; + + return 1; + } + return 0; +} + inline uint32_t aln_wlst(overlap_region *z, All_reads *rref, const ul_idx_t *uref, UC_Read* g_read, Correct_dumy* dumy, int32_t y_strand, int32_t y_id, int64_t x_start, int64_t x_end, long long y_start, int64_t block_s, double e_rate, int32_t is_cigar) { @@ -3570,7 +3792,9 @@ int32_t y_strand, int32_t y_id, int64_t x_start, int64_t x_end, long long y_star } return 0; } - +uint64_t realign_ed(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, +char *tstr, char *tstr_1, Correct_dumy* dumy, kvec_t_u64_warp* v_idx, int64_t block_s, double e_rate, +double *e_rate_final, uint32_t sec_check, int64_t *is_sort); inline void refine_ed_aln(overlap_region_alloc* overlap_list, All_reads *rref, const ul_idx_t *uref, UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, kvec_t_u64_warp* v_idx, int64_t block_s, double e_rate, double e_rate_final) { @@ -3581,11 +3805,12 @@ 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); + // z = &(overlap_list->list[j]); ovl = z->x_pos_e+1-z->x_pos_s; + // if(!realign_ed(z, uref, NULL, rref, g_read->seq, + // dumy->overlap_region, dumy->overlap_region_fix, dumy, v_idx, block_s, e_rate, 1, &is_srt)) { + // continue; // } + z = &(overlap_list->list[j]); z->is_match = 0; is_srt = 1; 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); @@ -3594,12 +3819,6 @@ 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; @@ -3630,7 +3849,7 @@ inline void refine_ed_aln(overlap_region_alloc* overlap_list, All_reads *rref, c } } mm_ws = z->x_pos_s; mm_aln = mm_we+1-mm_ws; - if(!simi_pass(ovl, mm_aln, rref, uref, NULL)) continue; + if(!simi_pass(ovl, mm_aln, uref?1:0, NULL)) continue; if(nw > 0 && w_idx[0] != (uint64_t)-1) mm_ws = z->w_list.a[w_idx[0]].x_end+1; for (i = 1; i < nw; i++) { //utilize the the start pos of next window in backward @@ -3667,13 +3886,13 @@ inline void refine_ed_aln(overlap_region_alloc* overlap_list, All_reads *rref, c } total_y_end = p->y_start - 1; } - if(!simi_pass(ovl, mm_aln, rref, uref, NULL)) break; + if(!simi_pass(ovl, mm_aln, uref?1:0, NULL)) break; } if(w_idx[i] != (uint64_t)-1) mm_ws = z->w_list.a[w_idx[i]].x_end+1; } if(i < nw) continue; - if(uref && simi_pass(ovl, z->align_length, NULL, uref, NULL)) { + if(uref && simi_pass(ovl, z->align_length, uref?1:0, NULL)) { z->is_match = 3; overlap_list->mapped_overlaps_length += z->align_length; ///sort for set_herror_win if(!is_srt) radix_sort_window_list_xs_srt(z->w_list.a, z->w_list.a + z->w_list.n); @@ -3697,7 +3916,7 @@ inline void refine_ed_aln(overlap_region_alloc* overlap_list, All_reads *rref, c ///debug_scan_cigar(&(overlap_list->list[j])); ///only calculate cigar for high quality overlaps // int64_t tt = 0; - if(simi_pass(ovl, z->align_length, rref, uref, &e_rate)) { + if(simi_pass(ovl, z->align_length, 0, &e_rate)) { a_nw = z->w_list.n; for (i = 0, is_srt = 1; i < a_nw; i++) { p = &(z->w_list.a[i]); @@ -9153,7 +9372,7 @@ void print_ovlp_occ_stat(overlap_region_alloc* overlap_list, uint32_t xlen, uint overlap_list->list[k].x_pos_s, overlap_list->list[k].x_pos_e+1); } } - +void align_ul_ed(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, char* qstr, char *tstr, double e_rate, int64_t w_l, void *km); void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref, UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g, Graph* DAGCon, Cigar_record* current_cigar, @@ -9170,6 +9389,7 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref 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); + // align_ul_ed(&(overlap_list->list[i]), uref, NULL, g_read->seq, dumy->overlap_region, max_ov_diff_ec, w_inf.window_length, km); } @@ -10354,4 +10574,428 @@ void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_I recalcate_high_het_overlap(overlap_list, R_INF, g_read, dumy, overlap_read); } -**/ \ No newline at end of file +**/ + +uint64_t update_ov_track_0(Fake_Cigar* z, overlap_region *o, int64_t apend_be, int64_t xl, int64_t yl, +k_mer_hit* hit, int64_t n_hit) +{ + int64_t k, dq, dr, dd, pdd = INT32_MAX, xr, yr; + z->length = 0; + if(hit[0].readID != o->y_id || hit[0].strand != o->y_pos_strand) return 0; + ///update o->s + o->x_pos_s = hit[0].self_offset; o->y_pos_s = hit[0].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; + } + + if(apend_be == 1) add_fake_cigar(z, o->x_pos_s, 0, NULL); + for (k = 0; (k < n_hit) && (hit[k].readID == o->y_id) && (hit[k].strand == o->y_pos_strand); k++) { + dq = hit[k].self_offset - o->x_pos_s; + dr = hit[k].offset - o->y_pos_s; + dd = dr - dq; + if(dd != pdd) { + pdd = dd; + add_fake_cigar(z, hit[k].self_offset, pdd, NULL); + } + } + ///update o->s + o->x_pos_e = hit[k-1].self_offset; o->y_pos_e = hit[k-1].offset; + 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; + } + + 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); + } + + return k; +} + +int64_t iter_hpc(uint8_t *m, int64_t mn, int64_t *mo, int64_t rev, +int64_t *so, int64_t *ho, int64_t sc) +{ + if(sc <= (*so)) return (*ho)-1; + if(!rev) { + while((*mo) < mn) { + for (; (*mo) < mn && m[(*mo)] == 255; (*mo)++) { + (*so) += m[(*mo)]; + } + (*so) += m[(*mo)]; (*mo)++; (*ho)++; + if(sc <= (*so)) return (*ho)-1; + } + } else { + while((*mo) < mn) { + (*so) += m[mn-(*mo)-1]; (*mo)++; (*ho)++; + for (; (*mo) < mn && m[mn-(*mo)-1] == 255; (*mo)++) { + (*so) += m[mn-(*mo)-1]; + } + if(sc <= (*so)) return (*ho)-1; + } + } + + return -1; +} + +uint64_t update_ov_track_hpc_0(Fake_Cigar* z, overlap_region *o, int64_t apend_be, int64_t xhl, +uint32_t *x_idx, int64_t y_idx_map_l, uint8_t *y_idx_map, hpc_t *hpc_g, k_mer_hit* hit, int64_t n_hit) +{ + int64_t k, dq, dr, dd, pdd = INT32_MAX, xr, yr, yhl, mo = 0, so = 0, ho = 0, x1, y1; + z->length = 0; + if(hit[0].readID != o->y_id || hit[0].strand != o->y_pos_strand) return 0; + yhl = hpc_len(*hpc_g, o->y_id); + + x1 = x_idx[hit[0].self_offset]; + y1 = iter_hpc(y_idx_map, y_idx_map_l, &mo, o->y_pos_strand, &so, &ho, hit[0].offset); + assert(y1 >= 0); + o->x_pos_s = x1; o->y_pos_s = y1; + + ///update o->s + 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; + } + + + if(apend_be == 1) add_fake_cigar(z, o->x_pos_s, 0, NULL); + for (k = 0; (k < n_hit) && (hit[k].readID == o->y_id) && (hit[k].strand == o->y_pos_strand); k++) { + x1 = x_idx[hit[k].self_offset]; + y1 = iter_hpc(y_idx_map, y_idx_map_l, &mo, o->y_pos_strand, &so, &ho, hit[k].offset); + assert(y1 >= 0); + dq = x1 - o->x_pos_s; + dr = y1 - o->y_pos_s; + dd = dr - dq; + if(dd != pdd) { + pdd = dd; + add_fake_cigar(z, x1, pdd, NULL); + } + } + ///update o->s + x1 = x_idx[hit[k-1].self_offset]; + y1 = iter_hpc(y_idx_map, y_idx_map_l, &mo, o->y_pos_strand, &so, &ho, hit[k-1].offset); + assert(y1 >= 0); + o->x_pos_e = x1; o->y_pos_e = y1; + xr = xhl-o->x_pos_e-1; yr = yhl-o->y_pos_e-1; + if(xr <= yr) { + o->x_pos_e = xhl-1; o->y_pos_e += xr; + } else { + o->y_pos_e = yhl-1; o->x_pos_e += yr; + } + + 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); + } + + return k; +} + +///(char *qstr, kvec_t_u64_warp* q_idx) -> only used for hpc +uint64_t update_ol_track(overlap_region_alloc* ol, Candidates_list *cl, hpc_t *hpc_g, const ul_idx_t *udb, uint32_t apend_be, uint64_t qlen, +char *qstr, kvec_t_u32_warp* q_idx) +{ + uint64_t cln = cl->length, i, k, l, m = 0; overlap_region *r; + if(hpc_g) { + if(ol->length) { + q_idx->a.n = 0; kv_resize(uint32_t, q_idx->a, qlen); m = 0; + for (l = 0, k = 1; k <= qlen; k++) { + if((k == qlen) || (qstr[k] != qstr[l]) || (seq_nt4_table[(uint8_t)qstr[l]] >= 4)) { + for (i = l; i < k; i++) q_idx->a.a[i] = m; + l = k; m++; + } + } + + for (i = k = 0; i < ol->length; ++i) { + r = &(ol->list[i]); + for (;(klist[k].readID!=r->y_id)||(cl->list[k].strand!=r->y_pos_strand)); k++); + // assert(kf_cigar), r, apend_be, m, q_idx->a.a, + (uint32_t)hpc_g->mm->idx[r->y_id], hpc_g->mm->a + (hpc_g->mm->idx[r->y_id]>>32), + hpc_g, cl->list+k, cln-k); + } + } + } else { + if(ol->length) { + for (i = k = 0; i < ol->length; ++i) { + r = &(ol->list[i]); + for (;(klist[k].readID!=r->y_id)||(cl->list[k].strand!=r->y_pos_strand)); k++); + // assert(kf_cigar), r, apend_be, qlen, udb->ug->u.a[r->y_id].len, cl->list+k, cln-k); + } + } + } + return m;///hpc length +} + +void inline resize_UC_Read(UC_Read *z, int64_t s) +{ + if(z->size < s) { + REALLOC(z->seq, s); z->size = s; + } +} + +uint64_t gen_hpc_str(const char *in, uint32_t in_l, UC_Read *z, uint64_t *in_hl) +{ + uint64_t hl, k, l; + if(in_hl) { + hl = (*in_hl); + } else { + for (l = hl = 0, k = 1; k <= in_l; k++) { + if((k == in_l) || (in[k] != in[l]) || (seq_nt4_table[(uint8_t)in[l]] >= 4)) { + hl++; l = k; + } + } + } + resize_UC_Read(z, hl); z->length = 0; + for (l = 0, k = 1; k <= in_l; k++) { + if((k == in_l) || (in[k] != in[l]) || (seq_nt4_table[(uint8_t)in[l]] >= 4)) { + z->seq[z->length++] = in[l]; l = k; + } + } + return hl; +} + +void align_ul_ed(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, char* qstr, char *tstr, double e_rate, int64_t w_l, void *km) +{ + int64_t q_s, q_e, nw, k, q_l, t_tot_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; z->is_match = 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); + 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); t_tot_l = hpc_g?hpc_len(*hpc_g, z->y_id):uref->ug->u.a[z->y_id].len; + if(init_waln(thre, t_s, t_tot_l, aln_l, &aux_beg, &aux_end, &t_s, &t_pri_l)) { + q_string = qstr+q_s; + t_string = return_str_seq(tstr, t_s, t_pri_l, z->y_pos_strand, hpc_g, uref, z->y_id, aux_beg, aux_end); + + t_end = Reserve_Banded_BPM(t_string, aln_l, q_string, q_l, thre, &error); + 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; + } + + assert(q_e == (int64_t)z->x_pos_e); +} + +uint64_t realign_ed(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, +char *tstr, char *tstr_1, Correct_dumy* dumy, kvec_t_u64_warp* v_idx, int64_t block_s, double e_rate, +double *e_rate_final, uint32_t sec_check, int64_t *is_sort) +{ + int64_t i, k, nw, a_nw, w_id, y_id, y_strand, real_y_start, x_start, x_end, x_len; + int64_t w_s, w_e, mm_we, mm_ws, mm_aln, ovl, y_readLen, total_y_start, total_y_end; + uint64_t *w_idx; window_list *p = NULL; if(sec_check && (!uref)) sec_check = 0; + z->is_match = 0; if(is_sort) (*is_sort) = 1; + if(z->w_list.n == 0) return 0; + 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); w_idx = v_idx->a.a; + memset(v_idx->a.a, -1, sizeof((*v_idx->a.a))*nw); + for (i = 0; i < a_nw; i++) { ///w_idx[] == (uint64_t) if unmatched + 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; + } + y_id = z->y_id; y_strand = z->y_pos_strand; + ovl = z->x_pos_e+1-z->x_pos_s; mm_we = z->x_pos_s; mm_aln = 0; + + if(hpc_g) y_readLen = hpc_len(*hpc_g, y_id); + else if(uref) y_readLen = uref->ug->u.a[y_id].len; + else y_readLen = Get_READ_LENGTH((*rref), y_id); + + for (i = a_nw-1; i >= 0; i--) { //utilize the the end pos of pre-window in forward + w_id = get_win_id_by_s(z, z->w_list.a[i].x_start, block_s, &w_e); + assert(z->w_list.a[i].x_end == w_e); + if(w_e > mm_we) mm_we = w_e; + ///in most cases, extra_begin = 0 + total_y_start = z->w_list.a[i].y_end + 1 - z->w_list.a[i].extra_begin; + for (k = w_id + 1; k < nw && total_y_start < y_readLen; k++) { + if(w_idx[k] != (uint64_t)-1) break; + w_s = w_e + 1; + w_id = get_win_id_by_s(z, w_s, block_s, &w_e); + assert(w_id == k); + x_start = w_s; x_end = w_e; + if(aln_wlst_adv(z, rref, hpc_g, uref, qstr, tstr, tstr_1, dumy, + y_strand, y_id, x_start, x_end, total_y_start, block_s, e_rate, 0)) { + p = &(z->w_list.a[z->w_list.n-1]); + w_idx[k] = z->w_list.n - 1; + if(x_end > mm_we) mm_we = x_end; + if(is_sort && (*is_sort) && z->w_list.n > 1 && p->x_start < z->w_list.a[z->w_list.n-2].x_start) (*is_sort) = 0; + } else { + break; + } + total_y_start = p->y_end + 1 - p->extra_begin; + } + } + mm_ws = z->x_pos_s; mm_aln = mm_we+1-mm_ws; + if(!simi_pass(ovl, mm_aln, sec_check, NULL)) return 0; + if(nw > 0 && w_idx[0] != (uint64_t)-1) mm_ws = z->w_list.a[w_idx[0]].x_end+1; + + for (i = 1; i < nw; i++) { //utilize the the start pos of next window in backward + ///find the first matched window, which should not be the first window + ///the pre-window of this matched window must be unmatched + if(w_idx[i] != (uint64_t)-1 && w_idx[i-1] == (uint64_t)-1) { + w_s = z->w_list.a[w_idx[i]].x_start; mm_aln -= (w_s-mm_ws); + ///check if the start pos of this matched window has been calculated + if(z->w_list.a[w_idx[i]].clen == 0) { + p = &(z->w_list.a[w_idx[i]]); + gen_backtrace_adv(p, z, rref, hpc_g, uref, qstr, tstr, tstr_1, dumy, y_strand, y_id); + assert(p->error != -1); + p->y_end += p->extra_begin; + } + real_y_start = p->y_start; + + ///the end pos for pre window is real_y_start - 1 + total_y_end = real_y_start - 1; + ///find the unmatched window on the left of current matched window + ///k starts from i - 1 + for (k = i - 1; k >= 0 && w_idx[k] == (uint64_t)-1 && total_y_end > 0; k--) { + w_e = w_s - 1; + w_id = get_win_id_by_e(z, w_e, block_s, &w_s); + assert(w_id == k); + x_start = w_s; x_end = w_e; x_len = x_end + 1 - x_start; + if(aln_wlst_adv(z, rref, hpc_g, uref, qstr, tstr, tstr_1, dumy, + y_strand, y_id, x_start, x_end, total_y_end+1-x_len, block_s, e_rate, 1)) { + p = &(z->w_list.a[z->w_list.n-1]); + p->y_start -= p->extra_begin; ///y_start has no shift, but y_end has shift + w_idx[k] = z->w_list.n - 1; + mm_aln += x_len; + if(is_sort && (*is_sort) && z->w_list.n > 1 && p->x_start < z->w_list.a[z->w_list.n-2].x_start) (*is_sort) = 0; + } else { + break; + } + total_y_end = p->y_start - 1; + } + if(!simi_pass(ovl, mm_aln, sec_check, NULL)) break; + } + if(w_idx[i] != (uint64_t)-1) mm_ws = z->w_list.a[w_idx[i]].x_end+1; + } + + if(i < nw) return 0; + + if(e_rate_final) { + /** + uint64_t srt; + if(simi_pass(ovl, z->align_length, 0, &e_rate)) { + a_nw = z->w_list.n; + for (i = 0, srt = 1; i < a_nw; i++) { + p = &(z->w_list.a[i]); + ///check if the cigar of this window has been got + if(p->clen == 0) { + gen_backtrace_adv(p, z, rref, hpc_g, uref, qstr, tstr, tstr_1, dumy, y_strand, y_id); + assert(p->error != -1); + } else { + p->y_end -= p->extra_begin; + } + if(srt && i > 0 && p->x_start < z->w_list.a[i-1].x_start) srt = 0; + } + if(!srt) radix_sort_window_list_xs_srt(z->w_list.a, z->w_list.a + z->w_list.n); + ///note: this function will change tstr/qstr + error_rate = non_trim_error_rate(z, rref, uref, v_idx, dumy, g_read, e_rate, block_s); + z->is_match = 0;///must be here; + + if (error_rate <= e_rate_final) { + overlap_list->mapped_overlaps_length += ovl; + z->is_match = 1; append_unmatched_wins(z, block_s); + if(rref) { + calculate_boundary_cigars(z, rref, dumy, g_read, e_rate); + } else { + calculate_ul_boundary_cigars(z, uref, dumy, g_read, e_rate, block_s); + } + // assert(get_num_wins(z->x_pos_s, z->x_pos_e+1, block_s)==(int64_t)z->w_list.n); + // assert((int64_t)z->x_pos_s==z->w_list.a[0].x_start && + // (int64_t)z->x_pos_e==z->w_list.a[z->w_list.n-1].x_end); + } else if (error_rate <= e_rate_final * 1.5) { + z->is_match = 3; + } + } + **/ + return 1; + } else { + return 1; + } + return 0; +} + + +void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, char *qstr, + uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, + haplotype_evdience_alloc* hap, kvec_t_u64_warp* v_idx, kvec_t_u32_warp* q_idx, + double e_rate, double eh_rate, int64_t wl, int64_t whl, void *km) +{ + uint64_t i, qhl, bs, k; Window_Pool w; double err; overlap_region t; + + ///hpc alignment + ///init hpc seq + qhl = update_ol_track(ol, cl, uref->hpc_g, uref, 1, ql, qstr, q_idx); + gen_hpc_str(qstr, ql, qu, &qhl); + ///verify hpc seq + clear_Correct_dumy(dumy, ol, km); + init_Window_Pool(&w, ql, whl, (int)(1.0/err)); + err = eh_rate; bs = (w.window_length)+(THRESHOLD_MAX_SIZE<<1)+1; + resize_UC_Read(tu, bs<<1); + for (i = k = 0; i < ol->length; i++) { + align_ul_ed(&(ol->list[i]), uref, uref->hpc_g, qu->seq, tu->seq, err, w.window_length, km); + if(!realign_ed(&(ol->list[i]), uref, uref->hpc_g, NULL, qu->seq, + tu->seq, tu->seq+bs, dumy, v_idx, w.window_length, err, &err, 0, NULL)) { + continue; + } + if(k != i) { + t = ol->list[k]; + ol->list[k] = ol->list[i]; + ol->list[i] = t; + } + k++; + } + ol->length = k; + + + + + // recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read); + // partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat); + // recalcate_window_ul_advance(overlap_list, uref, g_read, dumy, overlap_read, max_ov_diff_ec, w_inf.window_length, km); + // recalcate_window_advance(overlap_list, NULL, uref, g_read, dumy, overlap_read, v_idx, w_inf.window_length, max_ov_diff_ec, max_ov_diff_ec); + /** + refine_ed_aln(overlap_list, NULL, uref, g_read, dumy, overlap_read, v_idx, w_inf.window_length, max_ov_diff_ec, max_ov_diff_ec); + **/ + // fprintf(stderr, "[M::%s-beg] occ[0]->%lu, occ[1]->%lu, occ[2]->%lu, occ[3]->%lu\n", __func__, + // ovlp_occ(overlap_list, 0), ovlp_occ(overlap_list, 1), ovlp_occ(overlap_list, 2), ovlp_occ(overlap_list, 3)); + ///after this function, overlap_list is sorted by x_pos_e; used for g_chain + /** + partition_ul_overlaps_advance(overlap_list, uref, g_read, overlap_read, dumy, hap, force_repeat, max_ov_diff_ec, w_inf.window_length, km); + **/ + // print_ovlp_occ_stat(overlap_list, g_read->length, 1); + // print_ovlp_occ_stat(overlap_list, g_read->length, 2); + // fprintf(stderr, "[M::%s-end] occ[0]->%lu, occ[1]->%lu, occ[2]->%lu, occ[3]->%lu\n", __func__, + // ovlp_occ(overlap_list, 0), ovlp_occ(overlap_list, 1), ovlp_occ(overlap_list, 2), ovlp_occ(overlap_list, 3)); + // debug_phasing_status(overlap_list, uref->ug, 0, hap, g_read, 20, 1176); + // debug_phasing_status(overlap_list, uref->ug, 0, hap, g_read, 20, 1167); + // debug_phasing_status(overlap_list, uref->ug, 0, hap, g_read, 20, 1170); + /** + + + if(is_consensus) + { + generate_consensus(overlap_list, R_INF, g_read, dumy, g, DAGCon, current_cigar, second_round); + } + + + (*fully_cov) = check_if_fully_covered(overlap_list, R_INF, g_read, dumy, g, abnormal); + **/ +} \ No newline at end of file diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 7be0d7a..6a0fa41 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -78,7 +78,7 @@ int get_fake_gap_shift(Fake_Cigar* x, int index) return result; } -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) +void gen_fake_cigar(Fake_Cigar* z, overlap_region *o, int64_t apend_be, k_mer_hit* hit, int64_t n_hit) { int64_t k, dq, dr, dd, pdd; z->length = 0; if(apend_be == 1) add_fake_cigar(z, o->x_pos_s, 0, NULL); @@ -86,10 +86,10 @@ void gen_fake_cigar(Fake_Cigar* z, overlap_region *o, int64_t apend_be, k_mer_hi 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(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); @@ -103,7 +103,7 @@ void gen_fake_cigar(Fake_Cigar* z, overlap_region *o, int64_t apend_be, k_mer_hi 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); + gen_fake_cigar(z, o, apend_be, hit, n_hit); 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; @@ -116,7 +116,7 @@ void debug_cigar(Fake_Cigar* z, overlap_region *o, int64_t apend_be, k_mer_hit* 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); + // 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)); @@ -125,6 +125,7 @@ void debug_cigar(Fake_Cigar* z, overlap_region *o, int64_t apend_be, k_mer_hit* } } +///for backup 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) { @@ -146,6 +147,65 @@ int ovlp_chain_gen(overlap_region_alloc* ol, overlap_region* t, int64_t xl, int6 } } + int64_t xr, yr; + if(t->x_pos_s <= t->y_pos_s) { + t->y_pos_s -= t->x_pos_s; t->x_pos_s = 0; + } else { + t->x_pos_s -= t->y_pos_s; t->y_pos_s = 0; + } + + xr = xl-t->x_pos_e-1; yr = yl-t->y_pos_e-1; + if(xr <= yr) { + t->x_pos_e = xl-1; t->y_pos_e += xr; + } else { + t->y_pos_e = yl-1; t->x_pos_e += yr; + } + + overlap_region *o = &(ol->list[ol->length++]); + o->shared_seed = t->shared_seed; + o->align_length = 0; + o->is_match = 0; + o->non_homopolymer_errors = 0; + o->strong = 0; + o->x_id = t->x_id; + o->y_id = t->y_id; + o->x_pos_strand = 0;///always 0 + o->y_pos_strand = t->x_pos_strand; + + 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; + o->y_pos_e = yl-t->y_pos_s-1; o->y_pos_s = yl-t->y_pos_e-1; + } else { + o->x_pos_e = t->x_pos_e; o->x_pos_s = t->x_pos_s; + o->y_pos_e = t->y_pos_e; o->y_pos_s = t->y_pos_s; + } + ///debug + // debug_cigar(&(t->f_cigar), o, apend_be, hit, n_hit); + + return 1; +} + +int ovlp_chain_gen_fcigar(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; + ol->size = ol->length + 1; + kroundup64(ol->size); + REALLOC(ol->list, ol->size); + /// need to set new space to be 0 + memset(ol->list + sl, 0, sizeof(overlap_region)*(ol->size - sl)); + } + + if ((ol->length!=0) && (ol->list[ol->length-1].y_id==t->y_id)) { + if((ol->list[ol->length-1].shared_seed > t->shared_seed) || + ((ol->list[ol->length-1].shared_seed == t->shared_seed) && + (ol->list[ol->length-1].overlapLen <= t->overlapLen))) { + return 0; + } else { + ol->length--; + } + } + int64_t xr, yr, dd, pdd, dq, dr, id, i, fn; if(t->x_pos_s <= t->y_pos_s) { t->y_pos_s -= t->x_pos_s; t->x_pos_s = 0; @@ -1259,7 +1319,7 @@ int32_t lchain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, double bw_thres) dr = (int32_t)a[n_a-1].offset - (int32_t)a[0].offset; dd = ((dq>=dr)? (dq-dr): (dr-dq));//gap dg = ((dq>=dr)? (dr): (dq));///len - if (dd > (dg*bw_thres)) return -1; + if (dg == 0 || dd > (dg*bw_thres)) return -1; } for (i = 1; i < n_a; ++i) {///a[] is sorted by offset, instead of self_offset; but offset might be equal @@ -1275,6 +1335,7 @@ int32_t lchain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, double bw_thres) dr = (int32_t)a[i].offset - (int32_t)a[i-1].offset; dd = ((dq>=dr)? (dq-dr): (dr-dq));//gap dg = ((dq>=dr)? (dr): (dq));///len + if(dg == 0) break; tot_g += dd; if (dd > THRESHOLD_MAX_SIZE && dd > (dg*bw_thres)) break; @@ -1318,9 +1379,9 @@ inline int32_t comput_sc_ch(const k_mer_hit *ai, const k_mer_hit *aj, double bw_ ///ai is the suffix of aj int32_t dq, dr, dd, dg, q_span, sc; dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset); - if(dq < 0) return INT32_MIN; + if(dq <= 0) return INT32_MIN; dr = (int64_t)(ai->offset) - (int64_t)(aj->offset); - if(dr < 0) return INT32_MIN; + if(dr <= 0) return INT32_MIN; dd = dr > dq? dr - dq : dq - dr;//gap if((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, sl, ol))) return INT32_MIN; dg = dr < dq? dr : dq;//len @@ -1337,8 +1398,8 @@ inline int32_t comput_sc_ch(const k_mer_hit *ai, const k_mer_hit *aj, double bw_ } return sc; } - -uint64_t lchain_dp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, overlap_region* res, +///for backuo +uint64_t lchain_dp_fciagr(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) { @@ -1455,7 +1516,7 @@ uint64_t lchain_dp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, ov msc_i = i; i = p[i]; } } - res->overlapLen = get_chainLen(a[msc_i].self_offset, a[i].self_offset, xl, a[msc_i].offset, a[i].offset, yl); + 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]]; if(res->x_pos_strand) { int64_t hcl = cL>>1; k_mer_hit kp; @@ -1474,5 +1535,91 @@ uint64_t lchain_dp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, ov } } + return cL; +} + + +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) +{ + int64_t *p, *t, max_f, n_skip, st, max_j, end_j, sc, msc, msc_i, bw, max_ii, ovl, movl; + int32_t *f, max, tmp; int64_t i, j, ret, cL = 0; + resize_Chain_Data(dp, a_n, NULL); + 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(quick_check) { + ret = lchain_check(a, a_n, dp, bw_rate); + if (ret > 0) { + a_n = ret; msc_i = a_n-1; msc = f[msc_i]; + goto skip_ldp; + } + } + + memset(t, 0, (a_n*sizeof((*t)))); + for (i = st = 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; + + 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 || ((int64_t)a[i].offset) - ((int64_t)a[max_ii].offset) > max_dis) { + max = INT32_MIN; max_ii = -1; + for (j = i - 1; j >= st; --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) || (((((int64_t)a[i].offset)-((int64_t)a[max_ii].offset))<=max_dis) && (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; + } + } + } + + skip_ldp: + ///a[] has been sorted by offset + i = msc_i; + res->x_pos_s = res->x_pos_e = a[i].self_offset; + res->y_pos_s = res->y_pos_e = a[i].offset; + res->shared_seed = msc; + + cL = 0; + while (i >= 0) { + t[cL++] = i; msc_i = i; i = p[i]; + } + + 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]]; return cL; } \ No newline at end of file diff --git a/Hash_Table.h b/Hash_Table.h index e96a4ab..ecb503f 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -191,6 +191,7 @@ uint64_t lchain_dp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, ov 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, k_mer_hit* hit, int64_t n_hit); +void gen_fake_cigar(Fake_Cigar* z, overlap_region *o, 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 5a0ded8..7f60f13 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -249,12 +249,22 @@ typedef struct { uint64_t *idx; } cvert_t; +typedef struct { + size_t n, m; + uint8_t *a; + uint64_t *idx; +} hmap_t; + typedef struct { ma_ug_t *hg; size_t n, m; uint64_t *a; + hmap_t *mm; } hpc_t; +#define hpc_len(x, id) ((x).hg->u.a[(id)].len>>1) +#define hpc_str(x, id, rev) (((x).hg->u.a[(id)].s)+((rev)?((x).hg->u.a[(id)].len>>1):(0))) + typedef struct { ma_ug_t *ug; hpc_t *hpc_g; diff --git a/anchor.cpp b/anchor.cpp index 27c035b..890ab5c 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -26,6 +26,9 @@ KSORT_INIT(or_xs, overlap_region, oreg_xs_lt) #define oreg_ss_lt(a, b) ((a).shared_seed > (b).shared_seed) // in the decending order KSORT_INIT(or_ss, overlap_region, oreg_ss_lt) +#define oreg_id_lt(a, b) ((a).y_id < (b).y_id) +KSORT_INIT(or_id, overlap_region, oreg_id_lt) + typedef struct { int n; const ha_idxpos_t *a; @@ -869,16 +872,34 @@ 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 inline reverse_k_mer_hit(k_mer_hit *a, uint64_t a_n, uint64_t xl, uint64_t yl) +{ + uint64_t z, han = a_n>>1; k_mer_hit *ai, *aj, ka; + for (z = 0; z < han; z++) { + ai = &(a[z]); aj = &(a[a_n-z-1]); + ka = (*ai); (*ai) = (*aj); (*aj) = ka; + + ai->self_offset = xl-ai->self_offset-1; + ai->offset = yl-ai->offset-1; + + aj->self_offset = xl-aj->self_offset-1; + aj->offset = yl-aj->offset-1; + } + if(a_n&1) { + a[z].self_offset = xl-a[z].self_offset-1; + a[z].offset = yl-a[z].offset-1; + } +} + 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; + uint64_t i, k, l, m, sm, cn = cl->length, srt = 0; overlap_region *r; clear_overlap_region_alloc(ol); clear_fake_cigar(&(tf->f_cigar)); - ///calculate_overlap_region_by_chaining(cl, overlap_list, chain_idx, rid, rl, NULL, uref, bw_thres, keep_whole_chain, f_cigar); for (l = 0, k = 1, m = 0; k <= cn; k++) { if((k == cn) || (cl->list[k].readID != cl->list[l].readID) || (cl->list[k].strand != cl->list[l].strand)) { @@ -891,14 +912,14 @@ void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uin 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); - // 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)) { + if(gen_off) { + r = &(ol->list[ol->length-1]); + if(r->y_pos_strand) { + reverse_k_mer_hit(cl->list+m, sm, rl, rdb?Get_READ_LENGTH((*rdb), r->y_id):udb->ug->u.a[r->y_id].len); + } + gen_fake_cigar(&(r->f_cigar), r, apend_be, cl->list+m, sm); + } m += sm; } } @@ -908,13 +929,13 @@ void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uin cl->length = m; - + 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); + ks_introsort_or_ss(ol->length, ol->list); srt = 1; for (i = 0; i < ol->length; ++i) { - const overlap_region *r = &(ol->list[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; @@ -922,7 +943,7 @@ void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uin 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) { - overlap_region *r = &(ol->list[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))) { @@ -938,6 +959,27 @@ void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uin ol->length = k; } } + + if(!gen_off) { + if(srt) ks_introsort_or_id(ol->length, ol->list); + uint64_t cln = cl->length; + for (i = k = m = 0; i < ol->length; ++i) { + r = &(ol->list[i]); + for (;(klist[k].readID!=r->y_id)||(cl->list[k].strand!=r->y_pos_strand)); k++); + // assert(klist[k].readID!=r->y_id)||(cl->list[k].strand!=r->y_pos_strand)) break; + if(m != k) cl->list[m] = cl->list[k]; + } + // assert(m - cn > 0); + if(r->y_pos_strand) { + reverse_k_mer_hit(cl->list+cn, m-cn, rl, rdb?Get_READ_LENGTH((*rdb), r->y_id):udb->ug->u.a[r->y_id].len); + } + // gen_fake_cigar(&(r->f_cigar), r, apend_be, cl->list+cn, m-cn); + } + cl->length = m; + } + ks_introsort_or_xs(ol->length, ol->list); } diff --git a/inter.cpp b/inter.cpp index 8fff5b2..4ccd729 100644 --- a/inter.cpp +++ b/inter.cpp @@ -5236,7 +5236,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // 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); + s->opt->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1); clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2); @@ -9902,22 +9902,38 @@ 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) +void gen_hpc_seq(const char *in, uint32_t in_len, ma_utg_t *ou, hmap_t *mm, uint32_t hpc_id) { - uint32_t k, l, m; memset(ou, 0, sizeof((*ou))); + uint32_t k, l, m; int64_t n; 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; } } + kv_resize(uint8_t, *mm, mm->n+ou->len); + mm->idx[hpc_id] = ((uint64_t)mm->n)<<32; + ou->len <<= 1; 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]; + n = k-l; + while(n >= 0) {///even if n == 0, need to keep it + if(n < 255) { + kv_push(uint8_t, *mm, n); break; + } else { + kv_push(uint8_t, *mm, 255); n -= 255; + } + } l = k; } } + l = m; mm->idx[hpc_id] += (mm->n - (mm->idx[hpc_id]>>32)); + for (k = 0; k < l; k++) ou->s[m++] = RC_CHAR(ou->s[l-k-1]); + + + } void gen_microsatellite(const char *in, uint32_t in_len, ma_utg_t *idx, hpc_t *res, uint32_t mcs_len) @@ -9977,20 +9993,21 @@ uint32_t hpc_l(char *s, int64_t hof, int64_t sof, int64_t scut) 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, 1); CALLOC(p->hg, 1); kv_init(cc); CALLOC(p->mm, 1); CALLOC(p->hg->u.a, ug->u.n); p->hg->u.n = p->hg->u.m = ug->u.n; + CALLOC(p->mm->idx, 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_hpc_seq(cc.a, ug->u.a[k].len, &(p->hg->u.a[k]), p->mm, 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 + ng->seq[k].del = ng->seq[k].c = 0; ng->seq[k].len = p->hg->u.a[k].len>>1;///hpc len } memcpy(ng->arc, ug->g->arc, ng->n_arc*(sizeof((*(ng->arc))))); for (k = 1, i = 0; k <= ng->n_arc; k++) { @@ -10448,6 +10465,9 @@ void destroy_ul_idx_t(ul_idx_t *uu) if(uu->hpc_g) { free(uu->hpc_g->a); ma_ug_destroy(uu->hpc_g->hg); + free(uu->hpc_g->mm->a); + free(uu->hpc_g->mm->idx); + free(uu->hpc_g->mm); free(uu->hpc_g); } // if(uu->ov) {