diff --git a/CommandLines.cpp b/CommandLines.cpp index 38a9a0a..821cb9a 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -234,6 +234,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->s_hap_cov = 3; asm_opt->ul_error_rate = 0.2/**0.15**/; asm_opt->ul_error_rate_low = 0.1; + asm_opt->ul_error_rate_hpc = 0.2; asm_opt->ul_ec_round = 3; asm_opt->is_dbg_het_cnt = 0; } diff --git a/CommandLines.h b/CommandLines.h index dee0ca4..1590730 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -122,7 +122,7 @@ typedef struct { int64_t hg_size; float kpt_rate; int64_t infor_cov, s_hap_cov; - double ul_error_rate, ul_error_rate_low; + double ul_error_rate, ul_error_rate_low, ul_error_rate_hpc; int32_t ul_ec_round; uint8_t is_dbg_het_cnt; } hifiasm_opt_t; diff --git a/Correct.cpp b/Correct.cpp index 6530a79..cb0d5d6 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include "Correct.h" #include "Levenshtein_distance.h" #include "Assembly.h" @@ -1007,6 +1008,11 @@ unsigned int* get_error, int* get_y_end, int* get_x_end, int* get_aligned_xLen) fill_subregion_ul(dumy->overlap_region, y_beg, o_len, y_pos_strand, uref, y_id, extra_begin, extra_end); + // if(y_id == 6) { + // fprintf(stderr, "-[M::%s::aln_dir->%d] qs->%lld, ts->%lld, thres->%d, aux_beg->%d, aux_end->%d, t_pri_l->%lld\n", + // __func__, alignment_strand, x_beg, y_beg, threshold, extra_begin, extra_end, o_len); + // } + char* x_string = g_read->seq + x_beg; char* y_string = dumy->overlap_region; @@ -1056,10 +1062,10 @@ uint32_t get_init_paras(All_reads* rref, const ul_idx_t *uref, overlap_region *z return 1; } -int64_t check_coverage_gap(const kvec_t_u64_warp* v_idx, uint64_t w_s, uint64_t w_e, int64_t block_s) +int64_t check_coverage_gap(uint64_t *v_idx, uint64_t w_s, uint64_t w_e, int64_t block_s) { - int64_t wid = w_s/block_s, a_n = (uint32_t)(v_idx->a.a[wid]), k; - uint64_t *a = v_idx->a.a + (v_idx->a.a[wid]>>32); + int64_t wid = w_s/block_s, a_n = (uint32_t)(v_idx[wid]), k; + uint64_t *a = v_idx + (v_idx[wid]>>32); for (k = 0; k < a_n; k++) { if(((a[k]>>32) == w_s) && (((uint32_t)(a[k])) == w_e)) return 1; } @@ -1079,6 +1085,10 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* rref, const ul_i assert(w_s == z->w_list.a[k].x_start && w_id < idx_e && k <= w_id); 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(z->y_id == 1) { + // fprintf(stderr, "+[M::%s] ws->%d, we->%d, tot_l->%ld, tot_e->%ld\n", + // __func__, z->w_list.a[k].x_start, z->w_list.a[k].x_end, tLen, tErr); + // } // 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++) { @@ -1088,13 +1098,19 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* rref, const ul_i ///check if there are some windows that cannot be algined by any overlaps/unitigs ///if no, it is likely that the UL read itself has issues if(uref && v_idx && z->is_match == 4) { - if(check_coverage_gap(v_idx, w_s, w_e, block_s)) { + if(check_coverage_gap(v_idx->a.a, w_s, w_e, block_s)) { tErr += THRESHOLD_MAX_SIZE; + // if(z->y_id == 1) { + // fprintf(stderr, "-[M::%s] ws->%ld, we->%ld, tot_l->%ld, tot_e->%ld\n", __func__, w_s, w_e, tLen, tErr); + // } continue; } } if(!get_init_paras(rref, uref, z, w_s, w_e, e_rate, block_s, &y_s, &ex_beg, &ex_end, &err_thre)) { tErr += x_len; + // if(z->y_id == 1) { + // fprintf(stderr, "-[M::%s] ws->%ld, we->%ld, tot_l->%ld, tot_e->%ld\n", __func__, w_s, w_e, tLen, tErr); + // } continue; } p_err_thre = err_thre; @@ -1110,11 +1126,11 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* rref, const ul_i y_beg_left = y_beg_right = -1; if(m == w_id+1) { ///if the previous window is mapped - y_beg_left = z->w_list.a[k].y_end + 1; + y_beg_left = z->w_list.a[k].y_end + 1;///incorrect } if(m+1 == idx_e && k+1 < aw) { ///if the next window is mapped - y_beg_right = z->w_list.a[k+1].y_start-x_len; + y_beg_right = z->w_list.a[k+1].y_start-x_len;///incorrect } if(y_beg_left == -1 && y_beg_right == -1) { @@ -1164,6 +1180,13 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* rref, const ul_i else if(aligned_xLen_right != 0) { tErr += r_error_right + (x_len - aligned_xLen_right); } + // if(z->y_id == 1) { + // fprintf(stderr, "*[M::%s] qs->%ld, ts->%ld, tb[0]->%ld, tb[1]->%ld, di[0]->%u, di[1]->%u, al[0]->%d, al[1]->%d, err_thre->%ld\n", __func__, + // w_s, y_s, y_beg_left, y_beg_right, r_error_left, r_error_right, aligned_xLen_left, aligned_xLen_right, err_thre); + // } + // if(z->y_id == 1) { + // fprintf(stderr, "-[M::%s] ws->%ld, we->%ld, tot_l->%ld, tot_e->%ld\n", __func__, w_s, w_e, tLen, tErr); + // } } idx_e = w_id; } @@ -1176,8 +1199,11 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* rref, const ul_i ///check if there are some windows that cannot be algined by any overlaps/unitigs ///if no, it is likely that the UL read itself has issues if(uref && v_idx && z->is_match == 4) { - if(check_coverage_gap(v_idx, w_s, w_e, block_s)) { + if(check_coverage_gap(v_idx->a.a, w_s, w_e, block_s)) { tErr += THRESHOLD_MAX_SIZE; + // if(z->y_id == 1) { + // fprintf(stderr, "-[M::%s] ws->%ld, we->%ld, tot_l->%ld, tot_e->%ld\n", __func__, w_s, w_e, tLen, tErr); + // } continue; } // else { @@ -1188,6 +1214,9 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* rref, const ul_i } if(!get_init_paras(rref, uref, z, w_s, w_e, e_rate, block_s, &y_s, &ex_beg, &ex_end, &err_thre)) { tErr += x_len; + // if(z->y_id == 1) { + // fprintf(stderr, "-[M::%s] ws->%ld, we->%ld, tot_l->%ld, tot_e->%ld\n", __func__, w_s, w_e, tLen, tErr); + // } continue; } p_err_thre = err_thre; @@ -1257,6 +1286,9 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* rref, const ul_i else if(aligned_xLen_right != 0) { tErr += r_error_right + (x_len - aligned_xLen_right); } + // if(z->y_id == 1) { + // fprintf(stderr, "-[M::%s] ws->%ld, we->%ld, tot_l->%ld, tot_e->%ld\n", __func__, w_s, w_e, tLen, tErr); + // } } } @@ -3505,11 +3537,12 @@ 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, uint32_t second_ck, double *e_rate) +uint32_t inline simi_pass(int64_t ol, int64_t aln_ol, uint32_t second_ck, double o_rate, double *e_rate) { if(aln_ol == 0 || ol == 0) return 0; if((!second_ck) && (!e_rate)) { - if((ol*OVERLAP_THRESHOLD_FILTER) <= aln_ol) return 1; + // if((ol*OVERLAP_THRESHOLD_FILTER) <= aln_ol) return 1; + if((ol*o_rate) <= aln_ol) return 1; } else if(e_rate) { if((ol*((double)(((double)1.0)-(*e_rate)))) <= aln_ol) return 1; } else if(second_ck) { @@ -3615,9 +3648,18 @@ char *qstr, char *tstr, char *tstr1, Correct_dumy* dumy, uint32_t rev, uint32_t 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) { + // int32_t dbg_e = ed_band_cal_global(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres); + // if(dbg_e != (int32_t)error && dbg_e <= 4) { + // fprintf(stderr, "\n[M::%s::] error::%u, ed_global::%d\n", __func__, error, dbg_e); + // fprintf(stderr, "[tstr] %.*s\n", t_end+1-r_ts, t_string+r_ts); + // fprintf(stderr, "[qstr] %.*s\n", (int32_t)ql, q_string); + // } + // assert(dbg_e <= (int32_t)error); + // assert(ed_band_cal_global(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres) == + // ed_band_cal_global_128bit(t_string+r_ts, t_end+1-r_ts, q_string, ql, thres)); + ///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)) { @@ -3807,7 +3849,7 @@ inline void refine_ed_aln(overlap_region_alloc* overlap_list, All_reads *rref, c for (j = 0; j < on; ++j) { // 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)) { + // dumy->overlap_region, dumy->overlap_region_fix, dumy, v_idx, block_s, e_rate, NULL, 1, &is_srt)) { // continue; // } z = &(overlap_list->list[j]); z->is_match = 0; is_srt = 1; @@ -3849,7 +3891,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, uref?1:0, NULL)) continue; + if(!simi_pass(ovl, mm_aln, uref?1:0, OVERLAP_THRESHOLD_FILTER, 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 @@ -3886,13 +3928,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, uref?1:0, NULL)) break; + if(!simi_pass(ovl, mm_aln, uref?1:0, OVERLAP_THRESHOLD_FILTER, 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, uref?1:0, NULL)) { + if(uref && simi_pass(ovl, z->align_length, uref?1:0, OVERLAP_THRESHOLD_FILTER, 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); @@ -3916,7 +3958,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, 0, &e_rate)) { + if(simi_pass(ovl, z->align_length, 0, OVERLAP_THRESHOLD_FILTER, &e_rate)) { a_nw = z->w_list.n; for (i = 0, is_srt = 1; i < a_nw; i++) { p = &(z->w_list.a[i]); @@ -3969,6 +4011,81 @@ inline void refine_ed_aln(overlap_region_alloc* overlap_list, All_reads *rref, c } } +double test_err_rate(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, UC_Read* g_read, +char* qstr, char *tstr, char *tstr_1, Correct_dumy* dumy, kvec_t_u64_warp* v_idx, int64_t block_s, double e_rate) +{ + int64_t ovl = z->x_pos_e+1-z->x_pos_s, a_nw = z->w_list.n, i; double error_rate; + window_list *p; int64_t y_id = z->y_id, y_strand = z->y_pos_strand; + if(!simi_pass(ovl, z->align_length, 0, OVERLAP_THRESHOLD_FILTER, &e_rate)) return DBL_MAX; + + for (i = 0; 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; + } + } + error_rate = non_trim_error_rate(z, rref, uref, v_idx, dumy, g_read, e_rate, block_s); + return error_rate; +} + +uint32_t align_ul_ed_post(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, char* qstr, char *tstr, char *tstr_1, +Correct_dumy* dumy, double e_rate, int64_t w_l, double ovlp_cut, void *km); +double gen_extend_err(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, uint64_t *v_idx, int64_t block_s, double ovlp_cut, double e_rate, double e_max, int64_t *r_e); +inline void refine_ed_aln_test(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) +{ + int64_t j, on, ovl; uint64_t k; double rr; overlap_region *z; + + overlap_list->mapped_overlaps_length = 0; on = overlap_list->length; + for (j = 0; j < on; ++j) { + z = &(overlap_list->list[j]); ovl = z->x_pos_e+1-z->x_pos_s; + if(!align_ul_ed_post(z, uref, NULL, g_read->seq, dumy->overlap_region, dumy->overlap_region_fix, + dumy, e_rate, block_s, OVERLAP_THRESHOLD_FILTER, NULL)) { + continue; + } + if(uref && simi_pass(ovl, z->align_length, uref?1:0, OVERLAP_THRESHOLD_FILTER, NULL)) { + z->is_match = 3; overlap_list->mapped_overlaps_length += z->align_length; + } + } + + if(uref && overlap_list->mapped_overlaps_length > 0) { + set_herror_win(overlap_list, dumy, v_idx, e_rate, g_read->length, block_s); + } + + double e_max = e_rate_final * 1.5; + overlap_list->mapped_overlaps_length = 0; on = overlap_list->length; + for (j = 0; j < on; j++) { + z = &(overlap_list->list[j]); ovl = z->x_pos_e + 1 - z->x_pos_s; + rr = gen_extend_err(z, uref, NULL, rref, g_read->seq, dumy->overlap_region, dumy->overlap_region_fix, + dumy, v_idx?v_idx->a.a:NULL, block_s, -1, e_rate, (e_max+0.000001), NULL); + z->is_match = 0;///must be here; + if (rr <= e_rate_final) { + for (k = 0; k < z->w_list.n; k++) { + if(z->w_list.a[k].clen) continue; + gen_backtrace_adv(&(z->w_list.a[k]), z, rref, NULL, uref, g_read->seq, dumy->overlap_region, dumy->overlap_region_fix, + dumy, z->y_pos_strand, z->y_id); + } + + 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); + } + } else if (rr <= e_max) { + z->is_match = 3; + } + } +} + + + inline void add_base_to_correct_read_directly(Correct_dumy* dumy, char base) { @@ -9386,18 +9503,23 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref Window_Pool w_inf; 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); // align_ul_ed(&(overlap_list->list[i]), uref, NULL, g_read->seq, dumy->overlap_region, max_ov_diff_ec, w_inf.window_length, km); } + **/ // 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); + **/ + refine_ed_aln_test(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 @@ -10616,9 +10738,9 @@ k_mer_hit* hit, int64_t n_hit) 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) +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 == 0) return 0; if(sc <= (*so)) return (*ho)-1; if(!rev) { while((*mo) < mn) { @@ -10710,10 +10832,8 @@ char *qstr, kvec_t_u32_warp* q_idx) } 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, + r = &(ol->list[i]); k = r->non_homopolymer_errors; + update_ov_track_hpc_0(&(r->f_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); } @@ -10721,10 +10841,8 @@ char *qstr, kvec_t_u32_warp* q_idx) } 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); + r = &(ol->list[i]); k = r->non_homopolymer_errors; + update_ov_track_0(&(r->f_cigar), r, apend_be, qlen, udb->ug->u.a[r->y_id].len, cl->list+k, cln-k); } } } @@ -10759,6 +10877,640 @@ uint64_t gen_hpc_str(const char *in, uint32_t in_l, UC_Read *z, uint64_t *in_hl) return hl; } +///ts do not have aux_beg, while te has +uint32_t push_wlst(const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, overlap_region* ol, + char* qstr, char *tstr, char *tstr_1, Correct_dumy* dumy, + int64_t qs, int64_t qe, int64_t ts, int64_t te, int64_t tl, + int64_t error, int64_t aux_beg, int64_t aux_end, int64_t thres, double e_rate, + int64_t block_s, uint32_t sec_check, double ovlp_cut, void *km) +{ + + window_list p, t, *a; int64_t w_e, w_s, ce = qs - 1, cs = ol->x_pos_s, toff, ovl, ualn, aln; + uint64_t a_n, k; + + p.x_start = qs; p.x_end = qe; p.y_start = ts; p.y_end = te; p.error = error; + p.extra_begin = aux_beg; p.extra_end = aux_end; + p.error_threshold = thres; p.cidx = p.clen = 0; + if(ol->w_list.n > 0) { //utilize the the end pos of pre-window in forward + w_e = ol->w_list.a[ol->w_list.n-1].x_end; + toff = ol->w_list.a[ol->w_list.n-1].y_end + 1 - ol->w_list.a[ol->w_list.n-1].extra_begin; + while ((w_e < ce) && (toff < tl)) { + w_s = w_e + 1; + get_win_id_by_s(ol, w_s, block_s, &w_e); + // x_start = w_s; x_end = w_e; + if(aln_wlst_adv(ol, rref, hpc_g, uref, qstr, tstr, tstr_1, dumy, + ol->y_pos_strand, ol->y_id, w_s, w_e, toff, block_s, e_rate, 0)) { + toff = ol->w_list.a[ol->w_list.n-1].y_end + 1 - ol->w_list.a[ol->w_list.n-1].extra_begin; + } else { + break; + } + } + cs = ol->w_list.a[ol->w_list.n-1].x_end + 1; + } + ///utilize the the start pos of next window in backward + a_n = ol->w_list.n; w_s = qs; + if(w_s > cs) { + gen_backtrace_adv(&p, ol, rref, hpc_g, uref, qstr, tstr, tstr_1, dumy, ol->y_pos_strand, ol->y_id); + p.y_end += p.extra_begin; + toff = p.y_start - 1; + while ((w_s > cs) && (toff > 0)) { + w_e = w_s - 1; + get_win_id_by_e(ol, w_e, block_s, &w_s); + // x_start = w_s; x_end = w_e; x_len = x_end + 1 - x_start; + if(aln_wlst_adv(ol, rref, hpc_g, uref, qstr, tstr, tstr_1, dumy, + ol->y_pos_strand, ol->y_id, w_s, w_e, toff+1-(w_e+1-w_s), block_s, e_rate, 1)) { + ///y_start has no shift, but y_end has shift + ol->w_list.a[ol->w_list.n-1].y_start -= ol->w_list.a[ol->w_list.n-1].extra_begin; + toff = ol->w_list.a[ol->w_list.n-1].y_start - 1; + } else { + break; + } + } + } + + ol->align_length += qe + 1 - qs; + ovl = ol->x_pos_e+1-ol->x_pos_s; ualn = (qe + 1 - ol->x_pos_s) - ol->align_length; aln = ovl-ualn; + if((!simi_pass(ovl, aln, 0, ovlp_cut, &e_rate)) && (!simi_pass(ovl, aln, sec_check, ovlp_cut, NULL))) { + kv_push(window_list, ol->w_list, p); + return 0; + } + + if(ol->w_list.n > a_n) { + a = ol->w_list.a + a_n; a_n = ol->w_list.n - a_n; toff = a_n; a_n >>=1; + for (k = 0; k < a_n; k++) { + t = a[k]; a[k] = a[toff-1-k]; a[toff-1-k] = t; + } + } + kv_push(window_list, ol->w_list, p); + return 1; +} + +uint32_t align_ul_ed_post(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, char* qstr, char *tstr, char *tstr_1, +Correct_dumy* dumy, double e_rate, int64_t w_l, double ovlp_cut, void *km) +{ + int64_t q_s, q_e, nw, k, q_l, t_tot_l, sec_check = (uref&&(!hpc_g))?1:0; + 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; z->align_length = 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); + // int32_t debug_t_end, debug_error; + // debug_t_end = ed_band_cal_semi(t_string, aln_l, q_string, q_l, thre, &debug_error); + // if((t_end != debug_t_end) || (t_end >= 0 && debug_t_end >= 0 && debug_error != (int32_t)error)) { + // fprintf(stderr, "[M::%s] debug_error->%d, error->%d\n", __func__, debug_error, error); + // } + + + if (error!=((unsigned int)-1)) { + ///t_s do not have aux_beg, while t_s + t_end (aka, te) has + if(!push_wlst(uref, hpc_g, NULL, z, qstr, tstr, tstr_1, dumy, q_s, q_e, t_s, t_s + t_end, + t_tot_l, error, aux_beg, aux_end, thre, e_rate, w_l, sec_check, ovlp_cut, km)) { + return 0; + } + // 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((!simi_pass(z->x_pos_e+1-z->x_pos_s, z->align_length, 0, ovlp_cut, &e_rate)) && + (!simi_pass(z->x_pos_e+1-z->x_pos_s, z->align_length, sec_check, ovlp_cut, NULL))) return 0; + return 1; +} + + +inline uint32_t ed_cut(const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, +char *qstr, char *tstr, 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 aln_dir, int64_t* r_err, int64_t* qoff, int64_t* toff, int64_t* aln_qlen) +{ + (*aln_qlen) = 0; (*r_err) = INT32_MAX; + if(qoff) (*qoff) = -1; if(toff) (*toff) = -1; + int64_t ql, aln_l, t_tot_l, aux_beg, aux_end, t_pri_l, thres; + char *q_string, *t_string; unsigned int error; int t_end, q_end; + + 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(id == 6) { + // fprintf(stderr, "-[M::%s::aln_dir->%u] qs->%ld, ts->%ld, thres->%ld, aux_beg->%ld, aux_end->%ld, t_pri_l->%ld\n", + // __func__, aln_dir, qs, t_s, thres, aux_beg, aux_end, t_pri_l); + // } + if(aln_dir == 0) { + Reserve_Banded_BPM_Extension(t_string, aln_l, q_string, ql, thres, &error, &t_end, &q_end); + } else { + Reserve_Banded_BPM_Extension_REV(t_string, aln_l, q_string, ql, thres, &error, &t_end, &q_end); + } + + if(t_end != -1 && q_end != -1) (*aln_qlen) = (aln_dir?(ql-q_end):(q_end+1)); + if(qoff) (*qoff) = q_end; if(toff) (*toff) = t_end; (*r_err) = error; + + if((*aln_qlen) == 0) return 0; + return 1; +} + +int64_t gen_extend_err_0(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, uint64_t *v_idx, int64_t block_s, double e_rate, +int64_t qs, int64_t qe, int64_t pk) +{ + int64_t tot_e = 0, ts, di[2], al[2], tb[2], an = z->w_list.n; double rr; + int64_t id = z->y_id, rev = z->y_pos_strand, ql = qe + 1 - qs; + ///check if there are some windows that cannot be algined by any overlaps/unitigs + ///if no, it is likely that the UL read itself has issues + if(uref && v_idx && z->is_match == 4) { + if(check_coverage_gap(v_idx, qs, qe, block_s)) { + tot_e += THRESHOLD_MAX_SIZE; return tot_e; + } + } + ts = (qs - z->x_pos_s) + z->y_pos_s; ts += y_start_offset(qs, &(z->f_cigar)); + + di[0] = di[1] = al[0] = al[1] = 0; tb[0] = tb[1] = -1; + if((pk > 0) && (qs == (z->w_list.a[pk].x_end + 1))) { + if(z->w_list.a[pk].clen == 0) {///do not have cigar + gen_backtrace_adv(&(z->w_list.a[pk]), z, rref, hpc_g, uref, qstr, tstr, tstr_1, dumy, rev, id); + } + tb[0] = z->w_list.a[pk].y_end + 1; + } + + if(((pk+1) < an) && ((qe+1) == (z->w_list.a[pk+1].x_start))) { + if(z->w_list.a[pk+1].clen == 0) {///do not have cigar + gen_backtrace_adv(&(z->w_list.a[pk+1]), z, rref, hpc_g, uref, qstr, tstr, tstr_1, dumy, rev, id); + } + tb[1] = z->w_list.a[pk+1].y_start-ql; + } + + if(tb[0] == -1 && tb[1] == -1) tb[0] = tb[1] = ts; + else if(tb[0] == -1 && tb[1] != -1) tb[0] = tb[1]; + else if(tb[1] == -1 && tb[0] != -1) tb[1] = tb[0]; + + if(tb[0] != -1) { + if(!ed_cut(uref, hpc_g, rref, qstr, tstr, dumy, rev, id, qs, qe, tb[0], block_s, e_rate, + 0, &(di[0]), NULL, NULL, &(al[0]))) { + di[0] = ql; al[0] = 0; + } + } + + if(tb[1] != -1) { + if(!ed_cut(uref, hpc_g, rref, qstr, tstr, dumy, rev, id, qs, qe, tb[1], block_s, e_rate, + 1, &(di[1]), NULL, NULL, &(al[1]))) { + di[1] = ql; al[1] = 0; + } + } + + if(al[0] && al[1]) {///matched in both sides + if((al[0] + al[1]) <= ql) { + tot_e += di[0] + di[1] + ql - (al[0] + al[1]); + } else { + rr = ((double)ql)/((double)(al[0] + al[1])); + tot_e += (di[0] + di[1])*rr; + } + } else if((!al[0]) && (!al[1])) {//failed + tot_e += ql; + } else if(al[0]) { + tot_e += di[0] + (ql - al[0]); + }else if(al[1]) { + tot_e += di[1] + (ql - al[1]); + } + // if(z->y_id == 6) { + // fprintf(stderr, "-[M::%s] qs->%ld, ts->%ld, tb[0]->%ld, tb[1]->%ld, di[0]->%ld, di[1]->%ld, al[0]->%ld, al[1]->%ld, block_s->%ld, e_rate->%f\n", __func__, + // qs, ts, tb[0], tb[1], di[0], di[1], al[0], al[1], block_s, e_rate); + // } + return tot_e; +} + +double gen_extend_err(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, uint64_t *v_idx, int64_t block_s, double ovlp_cut, double e_rate, double e_max, int64_t *r_e) +{ + int64_t ovl, k, ce, an = z->w_list.n, tot_l, tot_e, ws, we, ql; + ovl = z->x_pos_e+1-z->x_pos_s; if(r_e) (*r_e) = INT64_MAX; + if(!simi_pass(ovl, z->align_length, 0, ovlp_cut, &e_rate)) return DBL_MAX; + // nw = get_num_wins(z->x_pos_s, z->x_pos_e+1, block_s); + for (k = 0; k < an; k++) { + if(z->w_list.a[k].clen) z->w_list.a[k].y_end -= z->w_list.a[k].extra_begin; + } + + tot_l = tot_e = 0; + for (k = an-1, ce = z->x_pos_e; k >= 0; k--) { + // assert(k == 0 || z->w_list.a[k].x_end > z->w_list.a[k-1].x_start);//sorted + tot_l += z->w_list.a[k].x_end + 1 - z->w_list.a[k].x_start; + tot_e += z->w_list.a[k].error;///matched window + + we = z->w_list.a[k].x_end; + while (we < ce) { + ws = we+1; + get_win_id_by_s(z, ws, block_s, &we); + ql = we+1-ws; tot_l += ql; + tot_e += gen_extend_err_0(z, uref, hpc_g, rref, qstr, tstr, tstr_1, dumy, v_idx, block_s, e_rate, ws, we, k); + if((e_max > 0) && (tot_e > (ovl*e_max))) return DBL_MAX; + } + ce = z->w_list.a[k].x_start-1; + if((e_max > 0) && (tot_e > (ovl*e_max))) return DBL_MAX; + } + + if(ce >= ((int64_t)z->x_pos_s)) { + we = ((int64_t)z->x_pos_s)-1; + while (we < ce) { + ws = we+1; + get_win_id_by_s(z, ws, block_s, &we); + ql = we+1-ws; tot_l += ql; + tot_e += gen_extend_err_0(z, uref, hpc_g, rref, qstr, tstr, tstr_1, dumy, v_idx, block_s, e_rate, ws, we, k); + if((e_max > 0) && (tot_e > (ovl*e_max))) return DBL_MAX; + } + } + + assert(tot_l == ovl); if(r_e) (*r_e) = tot_e; + return (double)(tot_e)/(double)(tot_l); +} + + +void push_anchors(window_list *z, window_list_alloc *zidx, asg64_v *anchor, uint64_t *qhp, int64_t qhp_l, int64_t *qhp_k, uint32_t mcl) +{ + int64_t xi = 0, yi = 0, ci, cn = z->clen; uint8_t c = (uint8_t)-1; uint32_t cl = (uint32_t)-1; + for (ci = 0; ci < cn; ci++) { + get_cigar_cell(z, zidx, ci, &c, &cl); + if (c == 0) { //match + if(cl >= mcl) { + ; + ; + ; + ; + } + xi += cl; yi += cl; + } else if (c == 1) { + xi += cl; yi += cl; + } else if (c == 2) {///y has more bases than x + yi += cl; + } else if (c == 3) {///x has more bases than y + xi += cl; + } + } +} + +#define gen_hpc_max_len(x) ((x)+((x)>>1)+1) +///[off_s, off_e) +uint64_t extract_mm_hpc(char *in, int64_t len, int64_t off_s, int64_t off_e, int64_t w, uint64_t rev) +{ + int64_t i, o, l, trim, k, tl; uint64_t m, sf; uint8_t c; + ///forward + for (k = 1, trim = 0; k <= w; k++) { + m = 0; o = gen_hpc_max_len(k); sf = k<<1; + if(!rev) { + ///[off_s, off_e) + for (i = ((off_s>=o)?(off_s-o):(0)), l = 0; i < off_e; i++) { + c = seq_nt4_table[(uint8_t)in[i]]; + if((c < 4) && (((l >= k) && (((m>>sf)&3) == c)) || (l < k))) { + if(l < k) m = (m<<2) + c; + else sf = (sf?(sf):(k<<1))-2; + l++; + } else { + if(i > off_s) { + tl = i-off_s; + if((l >= o) && (trim < tl)) trim = tl; + l = -1; + break; + } + l = 0; sf = k<<1; + } + } + tl = i-off_s; + if((l!=-1) && (i > off_s) && (l >= o) && (trim < tl)) { + trim = tl; + if(trim >= (off_e-off_s)) break; + } + } else { + ///[off_s, off_e) + for (i = (((len-off_e)>=o)?(off_e+o):(len))-1, l = 0; i >= off_s; i--) { + c = seq_nt4_table[(uint8_t)in[i]]; + if((c < 4) && (((l >= k) && (((m>>sf)&3) == c)) || (l < k))) { + if(l < k) m = (m<<2) + c; + else sf = (sf?(sf):(k<<1))-2; + l++; + } else { + if(i+1 < off_e) { + tl = off_e-i-1; + if((l >= o) && (trim < tl)) trim = tl; + l = -1; + break; + } + l = 0; sf = k<<1; + } + } + tl = off_e-i-1; + if((l!=-1) && (i+1 < off_e) && (l >= o) && (trim < tl)) { + trim = tl; + if(trim >= (off_e-off_s)) break; + } + } + } + return trim; +} + +uint64_t trim_hpc(const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, char *tstr, int64_t hpc_max, +int64_t ql, int64_t tl, int64_t tid, int64_t trev, int64_t *rqs, int64_t *rqe, int64_t *rts, int64_t *rte) +{ + ///[qs, qe); [ts, te) + int64_t qs = *rqs, qe = *rqe, ts = *rts, te = *rte, trim[2], hl, aux_l, subl = qe-qs; char *ss; + if(hpc_max > 32) hpc_max = 32; trim[0] = trim[1] = 0; aux_l = gen_hpc_max_len(hpc_max); + + qs -= aux_l; if(qs < 0) qs = 0; + qe += aux_l; if(qe > ql) qe = ql; + ss = qstr + qs; + hl = extract_mm_hpc(ss, qe - qs, (*rqs)-qs, (*rqe)-qs, hpc_max, 0); + if(hl >= subl) return 0; trim[0] = hl; + hl = extract_mm_hpc(ss, qe - qs, (*rqs)-qs, (*rqe)-qs, hpc_max, 1); + if(hl >= subl) return 0; trim[1] = hl; + if(trim[0] + trim[1] >= subl) return 0; + + ts -= aux_l; if(ts < 0) ts = 0; + te += aux_l; if(te > tl) te = tl; + if(rref) { + fill_subregion(tstr, ts, te-ts, trev, rref, tid, 0, 0); ss = tstr; + } else { + ss = return_str_seq(tstr, ts, te-ts, trev, hpc_g, uref, tid, 0, 0); + } + hl = extract_mm_hpc(ss, te - ts, (*rts)-ts, (*rte)-ts, hpc_max, 0); + if(hl >= subl) return 0; if(hl > trim[0]) trim[0] = hl; + hl = extract_mm_hpc(ss, te - ts, (*rts)-ts, (*rte)-ts, hpc_max, 1); + if(hl >= subl) return 0; if(hl > trim[1]) trim[1] = hl; + if(trim[0] + trim[1] >= subl) return 0; + + (*rqs) += trim[0]; (*rts) += trim[0]; + (*rqe) -= trim[1]; (*rte) -= trim[1]; + return 1; +} + +#define cl_pushp(type, v, p) do { \ + if ((v).length == (v).size) { \ + (v).size = (v).size? (v).size<<1 : 2; \ + (v).list = (type*)realloc((v).list, sizeof(type) * (v).size); \ + } \ + *(p) = &(v).list[(v).length++]; \ + } while (0) + +///ai is the suffix of aj +int64_t inline traceback_sc(const k_mer_hit *ai, const k_mer_hit *aj) +{ + int64_t qsi = ai->self_offset-ai->cnt, qej = aj->self_offset; + int64_t tsi = ai->offset-ai->cnt, tej = aj->offset; + if(qsi >= qej && tsi >= tej) return ai->cnt; + return INT32_MIN; +} + +void split_long_anchors(Candidates_list *ac, int64_t block, int64_t block_n) +{ + int64_t i, m, an = ac->length; + if(block_n < 0) { + for (i = block_n = 0; i < an; i++) { + if(ac->list[i].cnt <= block) block_n++; + else block_n += (ac->list[i].cnt/block) + (((ac->list[i].cnt%block) > 0)?1:0); + } + } + if(block_n <= an) return; + if(block_n > ac->size) { + ac->size = block_n; REALLOC(ac->list, ac->size); + } + for (i = an-1, m = block_n-1; i >= 0; i--) { + if(ac->list[i].cnt <= block) { + ac->list[m--] = ac->list[i]; + } else { + while (ac->list[i].cnt > 0) { + ac->list[m] = ac->list[i]; + if(ac->list[i].cnt >= block) { + ac->list[m].cnt = block; + ac->list[i].cnt -= block; + ac->list[i].self_offset -= block; + ac->list[i].offset -= block; + } else { + ac->list[i].cnt = 0; + } + m--; + } + } + } + ac->length = block_n; assert(m == -1); +} + +int64_t gen_affine_traceback_dp(Candidates_list *ac, int64_t max_skip, int64_t max_iter, int64_t max_dis, int64_t block, int64_t block_n) +{ + if(ac->length < 1) return 0; + int64_t i, j, *p, *t, max_f, n_skip, max_j, end_j, st, max_ii, sc, max, tmp, msc_i, msc; + int32_t *f, cL; k_mer_hit* a = ac->list; int64_t a_n = ac->length; Chain_Data* dp; + + for (i = 1; i < a_n; ++i) { + sc = traceback_sc(&a[i], &a[i-1]); + if(sc == INT32_MIN) break; + } + if(i >= a_n) return a_n; + + split_long_anchors(ac, block, block_n); + a = ac->list; a_n = ac->length; dp = &(ac->chainDP); + + resize_Chain_Data(dp, a_n, NULL); t = dp->tmp; f = dp->score; p = dp->pre; + t[0] = 0; p[0] = -1; f[0] = a[0].cnt; + msc_i = msc = -1; i = 0; + + memset(t, 0, (a_n*sizeof((*t)))); + for (i = st = 0, max_ii = -1; i < a_n; ++i) { + max_f = a[i].cnt; 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 = traceback_sc(&a[i], &a[j]); + if(sc == INT32_MIN) break; + 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) && ((((int64_t)a[i].offset)-((int64_t)a[j].offset))<=max_dis); --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 = traceback_sc(&a[i], &a[max_ii]); + 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) { + msc = f[i]; msc_i = i; + } + } + + cL = 0; i = msc_i; + while (i >= 0) { + t[cL++] = i; i = p[i]; + } + for (i = 0; i < cL; i++) a[i] = a[t[cL-i-1]]; + return cL; +} + +uint64_t gen_affine_traceback(overlap_region *o, const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, char* qstr, char *tstr, uint64_t ql, +uint64_t *qhp, uint64_t qhp_l, Candidates_list *ac, uint32_t hpc_max, uint32_t min_ach, uint32_t block) +{ + if(o->w_list.n <= 0) return 0; + int64_t nw = o->w_list.n, snw, k, t, qi, ti, ci, p_qi, p_ti, cn, qs, qe, ts, te, tl; + window_list *z; uint8_t c; uint32_t cl, pcl, id = o->y_id, rev = o->y_pos_strand; k_mer_hit *p; + uint64_t hm, aocc = 0, bocc = 0; qi = ti = 0; p_qi = p_ti = INT32_MIN; pcl = 0; + clear_Candidates_list(ac); + if(o->w_list.n > (uint64_t)ac->size) { + ac->size = o->w_list.n; REALLOC(ac->list, ac->size); + } + + if(hpc_g) tl = hpc_len(*hpc_g, id); + else if(uref) tl = uref->ug->u.a[id].len; + else tl = Get_READ_LENGTH((*rref), id); + for (k = 0; k < nw; k++) { + z = &(o->w_list.a[k]); ci = 0; cn = z->clen; + qi = z->x_start; ti = z->y_start; + for (ci = 0; ci < cn; ci++) { + get_cigar_cell(z, &(o->w_list), ci, &c, &cl); + if (c == 0) { //match + if((p_qi == qi) && (p_ti == ti)) { + pcl += cl; + } else { + ///push + if(pcl > 0) { + hm = 0; qs = qi - pcl; qe = qi; ts = ti - pcl; te = ti; + if(pcl > min_ach) hm = trim_hpc(uref, hpc_g, rref, qstr, tstr, hpc_max, ql, tl, id, rev, &qs, &qe, &ts, &te); + if(hm) aocc++; + if(hm || aocc == 0) { + cl_pushp(k_mer_hit, *ac, &p); + p->readID = p->strand = !!hm; + p->cnt = qe - qs; p->self_offset = qe; p->offset = te; + if(p->cnt > min_ach) bocc++; + } + } + ///push + pcl = 0; + } + qi += cl; ti += cl; + p_qi = qi; p_ti = ti; + } else { + ///push + if(pcl > 0) { + hm = 0; qs = qi - pcl; qe = qi; ts = ti - pcl; te = ti; + if(pcl > min_ach) hm = trim_hpc(uref, hpc_g, rref, qstr, tstr, hpc_max, ql, tl, id, rev, &qs, &qe, &ts, &te); + if(hm) aocc++; + if(hm || aocc == 0) { + cl_pushp(k_mer_hit, *ac, &p); + p->readID = p->strand = !!hm; + p->cnt = qe - qs; p->self_offset = qe; p->offset = te; + if(p->cnt > min_ach) bocc++; + } + } + pcl = 0; + ///push + if (c == 1) { + qi += cl; ti += cl; + } if (c == 2) {///t has more bases than p + ti += cl; + } if (c == 3) {///p has more bases than t + qi += cl; + } + } + } + } + + ///push + if(pcl > 0) { + hm = 0; qs = qi - pcl; qe = qi; ts = ti - pcl; te = ti; + if(pcl > min_ach) hm = trim_hpc(uref, hpc_g, rref, qstr, tstr, hpc_max, ql, tl, id, rev, &qs, &qe, &ts, &te); + if(hm) aocc++; + if(hm || aocc == 0) { + cl_pushp(k_mer_hit, *ac, &p); + p->readID = p->strand = !!hm; + p->cnt = qe - qs; p->self_offset = qe; p->offset = te; + if(p->cnt > min_ach) bocc++; + } + } + ///push + + nw = ac->length; snw = -1; + if(aocc > 0) { + if(aocc < (uint64_t)ac->length) { + for (k = t = snw = 0; k < nw; k++) { + if(ac->list[k].readID) { + ac->list[t] = ac->list[k]; + if(ac->list[t].cnt <= block) snw++; + else snw += (ac->list[t].cnt/block) + (((ac->list[t].cnt%block) > 0)?1:0); + t++; + } + } + ac->length = t; + } + } else if (bocc > 0) { + if(bocc < (uint64_t)ac->length) { + for (k = t = snw = 0; k < nw; k++) { + if(ac->list[k].cnt > min_ach) { + ac->list[t] = ac->list[k]; + if(ac->list[t].cnt <= block) snw++; + else snw += (ac->list[t].cnt/block) + (((ac->list[t].cnt%block) > 0)?1:0); + t++; + } + } + ac->length = t; + } + } + nw = ac->length; + gen_affine_traceback_dp(ac, 25, 5000, 5000, block, snw); + return 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) { int64_t q_s, q_e, nw, k, q_l, t_tot_l; @@ -10796,104 +11548,150 @@ void align_ul_ed(overlap_region *z, const ul_idx_t *uref, hpc_t *hpc_g, char* qs 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) +double *e_rate_final, uint32_t sec_check, double ovlp_cut, 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 i, k, nw, a_nw, w_id, y_id, y_strand, real_y_start, x_start, x_end, x_len, ce, cs; 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; + uint64_t *w_idx = NULL, srt = 1; 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; + ovl = z->x_pos_e+1-z->x_pos_s; mm_ws = 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 + + for (i = a_nw-1, ce = z->x_pos_e; i >= 0; i--) { //utilize the the end pos of pre-window in forward + w_e = mm_we = z->w_list.a[i].x_end; 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; + while ((w_e < ce) && (total_y_start < y_readLen)) { 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; + mm_we = x_end; } else { break; } total_y_start = p->y_end + 1 - p->extra_begin; } + ce = z->w_list.a[i].x_start-1; + if(i == a_nw-1) {///only possiblity with the largest end pos + mm_aln = mm_we+1-mm_ws; + if(!simi_pass(ovl, mm_aln, sec_check, ovlp_cut, NULL)) break; + } } - 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; + + if(z->w_list.a[a_nw-1].x_end > z->w_list.a[z->w_list.n-1].x_end) { + srt = 0; if(is_sort) (*is_sort) = 0; + } + if(i >= 0) return 0; + if((!srt) && (z->w_list.n <= (nw*0.2))) {///if very few windows are mapped + radix_sort_window_list_xs_srt(z->w_list.a, z->w_list.a + z->w_list.n); + srt = 1; if(is_sort) (*is_sort) = 1; + } + if(!srt) {///need sort + 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; + } + ///deal with first window + mm_ws = z->x_pos_s; + if(w_idx[0] != (uint64_t)-1) { + w_s = z->w_list.a[w_idx[0]].x_start; mm_aln -= (w_s-mm_ws); + 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; - 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); + ///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, ovlp_cut, 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; + } else {//sorted + a_nw = z->w_list.n; mm_ws = z->x_pos_s; + for (i = 0, cs = z->x_pos_s; i < a_nw; i++) { + p = &(z->w_list.a[i]); + w_s = p->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]]); + if((w_s > cs) && (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); 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--) { + while ((w_s > cs) && (total_y_end > 0)) { 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; + // 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(!simi_pass(ovl, mm_aln, sec_check, ovlp_cut, NULL)) break; + mm_ws = cs = z->w_list.a[i].x_end+1; } - if(w_idx[i] != (uint64_t)-1) mm_ws = z->w_list.a[w_idx[i]].x_end+1; + if(a_nw < (int64_t)z->w_list.n) { + srt = 0; if(is_sort) (*is_sort) = 0; + } + if(i < a_nw) return 0; } - 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++) { + for (i = 0; i < a_nw; i++) { p = &(z->w_list.a[i]); ///check if the cigar of this window has been got if(p->clen == 0) { @@ -10902,7 +11700,6 @@ double *e_rate_final, uint32_t sec_check, int64_t *is_sort) } 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 @@ -10932,29 +11729,40 @@ double *e_rate_final, uint32_t sec_check, int64_t *is_sort) return 0; } +uint64_t col_errors(overlap_region *z) +{ + uint64_t i, e = 0; + for (i = 0; i < z->w_list.n; i++) e += z->w_list.a[i].error; + return e; +} -void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, char *qstr, + +void ul_lalign_hpc(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) + double e_rate, double eh_rate, int64_t wl, void *km) { - uint64_t i, qhl, bs, k; Window_Pool w; double err; overlap_region t; - + uint64_t i, qhl, bs, k, ovl, whl; Window_Pool w; double err; overlap_region t; overlap_region *z; + whl = MIN((((double)THRESHOLD_MAX_SIZE)/eh_rate), WINDOW); + ol->mapped_overlaps_length = 0; + if(ol->length <= 0) return; + ///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; + clear_Correct_dumy(dumy, ol, km); err = eh_rate; + init_Window_Pool(&w, qhl, whl, (int)(1.0/err)); + 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)) { + if(!align_ul_ed_post(&(ol->list[i]), uref, uref->hpc_g, qu->seq, tu->seq, tu->seq+bs, dumy, err, w.window_length, OVERLAP_THRESHOLD_FILTER_HPC, km)) { continue; } + // fprintf(stderr, "+++[M::%s] yid::%u, x::[%u, %u), y::[%u, %u), aln::%u, err::%lu\n", __func__, ol->list[i].y_id, + // ol->list[i].x_pos_s, ol->list[i].x_pos_e+1, ol->list[i].y_pos_s, ol->list[i].y_pos_e+1, + // ol->list[i].align_length, col_errors(&(ol->list[i]))); if(k != i) { t = ol->list[k]; ol->list[k] = ol->list[i]; @@ -10963,9 +11771,55 @@ void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur k++; } ol->length = k; + if(ol->length <= 0) return; - + ///base alignment + update_ol_track(ol, cl, NULL, uref, 1, ql, NULL, NULL); + resize_UC_Read(qu, ql); qu->length = ql; memcpy(qu->seq, qstr, ql); + clear_Correct_dumy(dumy, ol, km); err = e_rate; + init_Window_Pool(&w, ql, wl, (int)(1.0/err)); + bs = (w.window_length)+(THRESHOLD_MAX_SIZE<<1)+1; + resize_UC_Read(tu, bs<<1); + for (i = 0; i < ol->length; i++) { + z = &(ol->list[i]); ovl = z->x_pos_e+1-z->x_pos_s; + if(!align_ul_ed_post(z, uref, NULL, qu->seq, tu->seq, tu->seq+bs, dumy, err, w.window_length, -1, km)) { + continue; + } + // fprintf(stderr, "---[M::%s] yid::%u, x::[%u, %u), y::[%u, %u), aln::%u, err::%lu\n", __func__, ol->list[i].y_id, + // ol->list[i].x_pos_s, ol->list[i].x_pos_e+1, ol->list[i].y_pos_s, ol->list[i].y_pos_e+1, + // ol->list[i].align_length, col_errors(&(ol->list[i]))); + if(uref && simi_pass(ovl, z->align_length, uref?1:0, -1, NULL)) { + z->is_match = 3; ol->mapped_overlaps_length += z->align_length; + } + } + + if(uref && ol->mapped_overlaps_length > 0) { + set_herror_win(ol, dumy, v_idx, err, ql, w.window_length); + } + + double e_max = err*1.5, rr; + for (i = k = 0; i < ol->length; i++) { + z = &(ol->list[i]); ovl = z->x_pos_e + 1 - z->x_pos_s; + rr = gen_extend_err(z, uref, NULL, NULL, qu->seq, tu->seq, tu->seq+bs, + dumy, v_idx?v_idx->a.a:NULL, w.window_length, -1, err, (e_max+0.000001), NULL); + z->is_match = 0;///must be here; + if (rr <= err) { + for (k = 0; k < z->w_list.n; k++) { + if(z->w_list.a[k].clen) continue; + gen_backtrace_adv(&(z->w_list.a[k]), z, NULL, NULL, uref, qu->seq, tu->seq, tu->seq+bs, + dumy, z->y_pos_strand, z->y_id); + } + + ol->mapped_overlaps_length += ovl; k++; + z->is_match = 1; append_unmatched_wins(z, w.window_length); + calculate_ul_boundary_cigars(z, uref, dumy, qu, err, w.window_length); + } else if (rr <= e_max) { + z->is_match = 3; + } + } + + partition_ul_overlaps_advance(ol, uref, qu, tu, dumy, hap, 1, err, w.window_length, km); // recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read); // partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat); @@ -10998,4 +11852,73 @@ void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur (*fully_cov) = check_if_fully_covered(overlap_list, R_INF, g_read, dumy, g, abnormal); **/ +} + +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, + double e_rate, int64_t wl, uint64_t is_base, void *km) +{ + uint64_t i, bs, k, ovl; Window_Pool w; double err; overlap_region t; overlap_region *z; + ol->mapped_overlaps_length = 0; + if(ol->length <= 0) return; + + ///base alignment + clear_Correct_dumy(dumy, ol, km); err = e_rate; + init_Window_Pool(&w, ql, wl, (int)(1.0/err)); + bs = (w.window_length)+(THRESHOLD_MAX_SIZE<<1)+1; + resize_UC_Read(tu, bs<<1); + + if(is_base) { + resize_UC_Read(qu, ql); qu->length = ql; memcpy(qu->seq, qstr, ql); + for (i = 0; i < ol->length; i++) { + z = &(ol->list[i]); ovl = z->x_pos_e+1-z->x_pos_s; + if(!align_ul_ed_post(z, uref, NULL, qu->seq, tu->seq, tu->seq+bs, dumy, err, w.window_length, -1, km)) { + continue; + } + if(uref && simi_pass(ovl, z->align_length, uref?1:0, -1, NULL)) { + z->is_match = 3; ol->mapped_overlaps_length += z->align_length; + } + } + + if(uref && ol->mapped_overlaps_length > 0) { + set_herror_win(ol, dumy, v_idx, err, ql, w.window_length); + } + + double e_max = err*1.5, rr; int64_t re; + for (i = k = 0; i < ol->length; i++) { + z = &(ol->list[i]); ovl = z->x_pos_e + 1 - z->x_pos_s; + rr = gen_extend_err(z, uref, NULL, NULL, qu->seq, tu->seq, tu->seq+bs, + dumy, v_idx?v_idx->a.a:NULL, w.window_length, -1, err, (e_max+0.000001), &re); + z->is_match = 0;///must be here; + if (rr <= err) { + if(k != i) { + t = ol->list[k]; + ol->list[k] = ol->list[i]; + ol->list[i] = t; + } + ol->list[k].is_match = 1; ol->list[k].non_homopolymer_errors = re; + k++; + } + } + + ol->length = k; + // fprintf(stderr, "+[M::%s] on::%lu\n", __func__, ol->length); + if(ol->length <= 0) return; + } else { + // fprintf(stderr, "-[M::%s] on::%lu\n", __func__, ol->length); + if(ol->length <= 1) return; + for (i = 0; i < ol->length; i++) { + z = &(ol->list[i]); ovl = z->x_pos_e+1-z->x_pos_s; z->is_match = 1; + for (k = 0; k < z->w_list.n; k++) { + if(z->w_list.a[k].clen) continue; + gen_backtrace_adv(&(z->w_list.a[k]), z, NULL, NULL, uref, qu->seq, tu->seq, tu->seq+bs, + dumy, z->y_pos_strand, z->y_id); + } + ol->mapped_overlaps_length += ovl; + append_unmatched_wins(z, w.window_length); + calculate_ul_boundary_cigars(z, uref, dumy, qu, err, w.window_length); + } + partition_ul_overlaps_advance(ol, uref, qu, tu, dumy, hap, 1, err, w.window_length, km); + } } \ No newline at end of file diff --git a/Correct.h b/Correct.h index 801e85f..76cd241 100644 --- a/Correct.h +++ b/Correct.h @@ -1130,6 +1130,10 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref kvec_t_u64_warp* v_idx, window_list_alloc* win_ciagr_buf, int force_repeat, int is_consensus, int* fully_cov, int* abnormal, double max_ov_diff_ec, long long winLen, void *km); +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, + double e_rate, int64_t wl, uint64_t is_base, void *km); void lchain_align(overlap_region_alloc* overlap_list, const ul_idx_t *uref, UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, @@ -1267,6 +1271,7 @@ inline void push_cigar_cell(window_list_alloc *res, uint8_t c, uint32_t len) uint16_t p = c; p <<= 14; p += (uint16_t)len; kv_push(uint16_t, res->c, p); } +int64_t get_num_wins(int64_t s, int64_t e, int64_t block_s); #define FORWARD_KSW 0 #define BACKWARD_KSW 1 diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 6a0fa41..0da0d87 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -1441,7 +1441,7 @@ uint64_t lchain_dp_fciagr(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* 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) { + for (j = i - 1; (j >= st) && ((((int64_t)a[i].offset)-((int64_t)a[j].offset))<=max_dis); --j) { if (max < f[j]) { max = f[j], max_ii = j; } @@ -1581,7 +1581,7 @@ uint64_t lchain_dp(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, ov 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) { + for (j = i - 1; (j >= st) && ((((int64_t)a[i].offset)-((int64_t)a[j].offset))<=max_dis); --j) { if (max < f[j]) { max = f[j], max_ii = j; } diff --git a/Hash_Table.h b/Hash_Table.h index ecb503f..a5ec529 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -12,6 +12,7 @@ #define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25 #define THRESHOLD 15 #define OVERLAP_THRESHOLD_FILTER 0.9 +#define OVERLAP_THRESHOLD_FILTER_HPC 0.75 #define HIGH_HET_OVERLAP_THRESHOLD_FILTER 0.3 #define HIGH_HET_ERROR_RATE 0.08 #define THRESHOLD_MAX_SIZE 31 diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index 5e2a7bf..66133f7 100644 --- a/Levenshtein_distance.h +++ b/Levenshtein_distance.h @@ -8,9 +8,12 @@ #include #include #include +#include "kvec.h" +extern const unsigned char seq_nt4_table[256]; typedef uint64_t Word; typedef uint32_t Word_32; +typedef struct {size_t n, m; uint16_t *a; } asg16_v; inline void get_error(int t_length, int errthold, int init_err, Word VP, Word VN, unsigned int* return_err, int* back_site) @@ -208,6 +211,151 @@ unsigned int* return_err, int* return_p_end, int* return_t_end) return (*return_t_end); } +inline int Reserve_Banded_BPM_Extension_REV +(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, +unsigned int* return_err, int* return_p_end, int* return_t_end) +{ + (*return_err) = (unsigned int)-1; + (*return_p_end) = -1; + (*return_t_end) = -1; + + Word Peq[256]; + + unsigned int line_error = (unsigned int)-1; + int return_site; + int band_length = (errthold << 1) + 1; + int i = 0; + Word tmp_Peq_1 = (Word)1; + + Peq[(uint8_t)'A'] = (Word)0; + Peq[(uint8_t)'T'] = (Word)0; + Peq[(uint8_t)'G'] = (Word)0; + Peq[(uint8_t)'C'] = (Word)0; + + + Word Peq_A; + Word Peq_T; + Word Peq_C; + Word Peq_G; + + ///band_length = 2k + 1 + for (i = 0; i> 1; + VN = X&HP; + VP = HN | ~(X | HP); + + if (!(D0&err_mask)) + { + ++err; + if ((err - last_high)>errthold) + { + return (*return_t_end); + } + } + get_error(i + 1, errthold, err, VP, VN, &line_error, &return_site); + if(line_error != (unsigned int)-1) + { + (*return_t_end) = t_length-i-1; + (*return_p_end) = p_length-return_site-1; + (*return_err) = line_error; + } + + Peq[(uint8_t)'A'] = Peq[(uint8_t)'A'] >> 1; + Peq[(uint8_t)'C'] = Peq[(uint8_t)'C'] >> 1; + Peq[(uint8_t)'G'] = Peq[(uint8_t)'G'] >> 1; + Peq[(uint8_t)'T'] = Peq[(uint8_t)'T'] >> 1; + + + ++i; + ++i_bd; + Peq[(uint8_t)pattern[p_length-i_bd-1]] = Peq[(uint8_t)pattern[p_length-i_bd-1]] | Mask; + } + + + + + + X = Peq[(uint8_t)text[t_length-i-1]] | VN; + D0 = ((VP + (X&VP)) ^ VP) | X; + HN = VP&D0; + HP = VN | ~(VP | D0); + X = D0 >> 1; + VN = X&HP; + VP = HN | ~(X | HP); + if (!(D0&err_mask)) + { + ++err; + if ((err - last_high)>errthold) + { + return (*return_t_end); + } + } + ///i = t_length - 1 + get_error(i + 1, errthold, err, VP, VN, &line_error, &return_site); + if(line_error != (unsigned int)-1) + { + (*return_t_end) = t_length-i-1; + (*return_p_end) = p_length-return_site-1; + (*return_err) = line_error; + } + + return (*return_t_end); +} + inline void reverse_string(char* str, int strLen) { int i, Len; @@ -267,7 +415,494 @@ int* return_t_end, int* return_aligned_t_len) } } +///p_length might be samller than t_length + 2 * errthold +/// pattern is longer than text +inline int32_t ed_band_cal_semi(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t *re_err) +{ + (*re_err) = INT32_MAX; + Word Peq[5] = {0}, mm = (Word)1, VP = 0, VN = 0, X = 0, D0 = 0, HN = 0, HP = 0; + int32_t bd = (thre<<1)+1, i, err = 0, i_bd = (thre<<1), last_high = (thre<<1), tn0 = tn - 1; + int32_t cut = thre+last_high; + for (i = 0; i < bd; i++) { + Peq[seq_nt4_table[(uint8_t)pstr[i]]] |= mm; mm <<= 1; + } + ///should make Peq[4] = 0 if N is always an error + Peq[4] = 0; + i = 0; mm = ((Word)1 << (thre<<1));///for the incoming char/last char + + while (i < tn0) { + X = Peq[seq_nt4_table[(uint8_t)tstr[i]]] | VN; + + D0 = ((VP + (X&VP)) ^ VP) | X; + + HN = VP&D0; + HP = VN | ~(VP | D0); + + X = D0 >> 1; + VN = X&HP; + VP = HN | ~(X | HP); + + if (!(D0&(1ULL))) { + ++err; + if (err>cut) return -1; + } + + Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1; ///Peq[4] >>= 1; + + ++i; ++i_bd; + Peq[seq_nt4_table[(uint8_t)pstr[i_bd]]] |= mm; Peq[4] = 0; + } + + X = Peq[seq_nt4_table[(uint8_t)tstr[i]]] | VN; + D0 = ((VP + (X&VP)) ^ VP) | X; + HN = VP&D0; + HP = VN | ~(VP | D0); + X = D0 >> 1; + VN = X&HP; + VP = HN | ~(X | HP); + if (!(D0&(1ULL))) { + ++err; + if (err>cut) return -1; + } + + int32_t site = tn - 1, end = -1;///up bound + ///in most cases, ai = (thre<<1) + int32_t ai = pn - tn, uge = INT32_MAX; + if ((err <= thre) && (err<=(*re_err))) { + *re_err = err; end = site; + } + i = 0; + + while (i < ai) { + err += ((VP >> i)&(1ULL)); err -= ((VN >> i)&(1ULL)); ++i; + if ((err <= thre) && (err <= (*re_err))) { + *re_err = err; end = site + i; + } + if(i == thre) uge = err; + } + + if((uge<=thre) && (uge == (*re_err))) end = site + thre; + return end; +} + +inline void print_bit(Word z, int64_t w, const char *cmd) +{ + int64_t k;//, w = (sizeof(Word)<<3); + fprintf(stderr, "%s\t", cmd); + for (k = 0; k < w; k++) fprintf(stderr, "%llu", (z>>k)&(1ULL)); + fprintf(stderr, "\n"); +} + +inline void print_bits(Word *az, int64_t w, const char *cmd) +{ + int64_t k, m, s = (sizeof(*az)<<3), sw = (w/s) + (!!(w%s)), ks; + fprintf(stderr, "%s\t", cmd); + for (m = k = 0; m < sw && k < w; m++) { + for (ks = 0; ks < s && k < w; ks++, k++) fprintf(stderr, "%llu", (az[m]>>ks)&(1ULL)); + } + fprintf(stderr, "\n"); +} + +inline int32_t ed_band_cal_global(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre) +{ + if((pn > tn + thre) || (tn > pn + thre)) return INT32_MAX; + if((pn < thre + 1) || (tn < thre + 1)) return INT32_MAX; + Word Peq[5] = {0}, mm, VP = 0, VN = 0, X = 0, D0 = 0, HN = 0, HP = 0; + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd = thre+1, i_bd = thre; + // fprintf(stderr, "\n[M::%s::]\n", __func__); + for (i = 0, mm = (((Word)1)<> 1; + VN = X&HP; + VP = HN | ~(X | HP); + // fprintf(stderr, "\n[M::%s::i->%d]\n", __func__, i); + // print_bit(VN, (thre<<1)+1, "VN"); + // print_bit(VP, (thre<<1)+1, "VP"); + // print_bit(HN, (thre<<1)+1, "HN"); + // print_bit(HP, (thre<<1)+1, "HP"); + // print_bit(D0, (thre<<1)+1, "D0"); + + if (!(D0&(1ULL))) { + ++err; + if (err>cut) return INT32_MAX; + } + // fprintf(stderr, "[M::%s::i->%d] err->%d\n", __func__, i, err); + + Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1; ///Peq[4] >>= 1; + + ++i; ++i_bd; + if(i_bd < pn) { + Peq[seq_nt4_table[(uint8_t)pstr[i_bd]]] |= mm; Peq[4] = 0; + } + // if(i < pn) Peq[seq_nt4_table[(uint8_t)pstr[i]]] |= mm; + } + + X = Peq[seq_nt4_table[(uint8_t)tstr[i]]] | VN; + D0 = ((VP + (X&VP)) ^ VP) | X; + HN = VP&D0; + HP = VN | ~(VP | D0); + X = D0 >> 1; + VN = X&HP; + VP = HN | ~(X | HP); + // fprintf(stderr, "\n[M::%s::i->%d]\n", __func__, i); + // print_bit(VN, (thre<<1)+1, "VN"); + // print_bit(VP, (thre<<1)+1, "VP"); + // print_bit(HN, (thre<<1)+1, "HN"); + // print_bit(HP, (thre<<1)+1, "HP"); + // print_bit(D0, (thre<<1)+1, "D0"); + if (!(D0&(1ULL))) { + ++err; + if (err>cut) return INT32_MAX; + } + // fprintf(stderr, "[M::%s::i->%d] err->%d\n", __func__, i, err); + + int32_t site = tn - 1 - thre;///up bound + for (cut = pn - 1, i = 0; site < cut; site++, i++) { + // fprintf(stderr, "+[M::%s::site->%d] err->%d\n", __func__, site, err); + err += ((VP >> i)&(1ULL)); err -= ((VN >> i)&(1ULL)); + // fprintf(stderr, "-[M::%s::site->%d] err->%d\n", __func__, site, err); + } + + if (site == cut && err <= thre) return err; + return INT32_MAX; +} + +typedef uint64_t w_sig; +typedef struct {w_sig a[2];} w128_t; +#define bitw (6) +#define bitwbit (64) +#define bitz (63) +// typedef uint32_t w_sig; +// typedef struct {w_sig a[2];} w128_t; +// #define bitw (5) +// #define bitwbit (32) +// #define bitz (31) + +#define w128_bit(x, b) ((x).a[((b)>>bitw)]|=(((w_sig)1)<<((b)&bitz))) +#define w128_clear(x) ((x).a[0]=(x).a[1]=0) + +#define w128_self_not(x) ((x).a[0]=~(x).a[0], \ + (x).a[1]=~(x).a[1]) + +#define w128_self_or(x, y) ((x).a[0]|=(y).a[0], \ + (x).a[1]|=(y).a[1]) + +#define w128_or(r, x, y) ((r).a[0] = (x).a[0]|(y).a[0], \ + (r).a[1] = (x).a[1]|(y).a[1]) + +#define w128_and(r, x, y) ((r).a[0] = (x).a[0]&(y).a[0], \ + (r).a[1] = (x).a[1]&(y).a[1]) + +#define w128_self_xor(x, y) ((x).a[0]^=(y).a[0], \ + (x).a[1]^=(y).a[1]) + +// #define w128_self_lsft_l(x, l) ((x).a[1] = ((x).a[1]<<(l))|((x).a[0]>>(bitwbit-(l))), \ +// (x).a[0] <<= (l)) + +#define w128_self_lsft_1(x) ((x).a[1] = ((x).a[1]<<1)|((x).a[0]>>bitz), \ + (x).a[0] <<= 1) + +#define w128_self_rsft_1(x) ((x).a[0] = ((x).a[0]>>1)|((x).a[1]<>= 1) + +#define w128_self_add(x, y) ((x).a[0]+=(y).a[0], \ + (x).a[1]+=(y).a[1]+((x).a[0]<(y).a[0])) + +#define w128_set_bit_lsub(x, l) do { \ + (x).a[0] = (w_sig)-1, (x).a[1] = 0; \ + if((l) <= bitwbit) (x).a[0] = (((w_sig)1)<<(l))-1; \ + else (x).a[1] = (((w_sig)1)<<((l)-bitwbit))-1;\ + } while (0) \ + +#define ed_core_w128(RE) { \ + /**X = Peq[seq_nt4_table[(uint8_t)tstr[i]]] | VN;**/\ + c = seq_nt4_table[(uint8_t)tstr[i]]; w128_or(X, Peq[c], VN);\ + /**D0 = ((VP + (X&VP)) ^ VP) | X;**/\ + w128_and(D0, X, VP);\ + w128_self_add(D0, VP);\ + w128_self_xor(D0, VP);\ + w128_self_or(D0, X);\ + /**HN = VP&D0;**/\ + w128_and(HN, VP, D0);\ + /**HP = VN | ~(VP | D0);**/\ + w128_or(HP, VP, D0);\ + w128_self_not(HP);\ + w128_self_or(HP, VN);\ + /**X = D0 >> 1;**/\ + X = D0; w128_self_rsft_1(X);\ + /**VN = X&HP;**/\ + w128_and(VN, X, HP);\ + /**VP = HN | ~(X | HP);**/\ + w128_or(VP, X, HP);\ + w128_self_not(VP);\ + w128_self_or(VP, HN);\ + /**if (!(D0&(1ULL)))**/\ + if (!(D0.a[0]&(1ULL))) {\ + ++err;\ + if (err>cut) return RE;\ + }\ + /** Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1;**/\ + w128_self_rsft_1(Peq[0]); w128_self_rsft_1(Peq[1]);\ + w128_self_rsft_1(Peq[2]); w128_self_rsft_1(Peq[3]);\ + } + +inline int32_t ed_band_cal_global_128bit(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre) +{ + if((pn > tn + thre) || (tn > pn + thre)) return INT32_MAX; + if((pn < thre + 1) || (tn < thre + 1)) return INT32_MAX; + w128_t Peq[5], mm, VP, VN, X, D0, HN, HP; uint8_t c; + int32_t i, err, tn0 = tn - 1, cut = thre+(thre<<1), bd = thre+1, i_bd = thre; + w128_clear(Peq[0]); w128_clear(Peq[1]); w128_clear(Peq[2]); w128_clear(Peq[3]); w128_clear(Peq[4]); + + w128_clear(mm); w128_bit(mm, thre); ///mm = (((Word)1)<> 1; + X = D0; w128_self_rsft_1(X); + // VN = X&HP; + w128_and(VN, X, HP); + // VP = HN | ~(X | HP); + w128_or(VP, X, HP); + w128_self_not(VP); + w128_self_or(VP, HN); + + + //if (!(D0&(1ULL))) + if (!(D0.a[0]&(1ULL))) { + ++err; + if (err>cut) return INT32_MAX; + } + // fprintf(stderr, "[M::%s::i->%d] err->%d\n", __func__, i, err); + + // Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1; + w128_self_rsft_1(Peq[0]); w128_self_rsft_1(Peq[1]); + w128_self_rsft_1(Peq[2]); w128_self_rsft_1(Peq[3]); + + ++i; ++i_bd; + if(i_bd < pn) { + c = seq_nt4_table[(uint8_t)pstr[i_bd]]; + ///if(c < 4) Peq[c] |= mm; + if(c < 4) w128_self_or(Peq[c], mm); + } + } + + // X = Peq[seq_nt4_table[(uint8_t)tstr[i]]] | VN; + c = seq_nt4_table[(uint8_t)tstr[i]]; w128_or(X, Peq[c], VN); + // D0 = ((VP + (X&VP)) ^ VP) | X; + w128_and(D0, X, VP); + w128_self_add(D0, VP); + w128_self_xor(D0, VP); + w128_self_or(D0, X); + // HN = VP&D0; + w128_and(HN, VP, D0); + // HP = VN | ~(VP | D0); + w128_or(HP, VP, D0); + w128_self_not(HP); + w128_self_or(HP, VN); + // X = D0 >> 1; + X = D0; w128_self_rsft_1(X); + // VN = X&HP; + w128_and(VN, X, HP); + // VP = HN | ~(X | HP); + w128_or(VP, X, HP); + w128_self_not(VP); + w128_self_or(VP, HN); + + // if (!(D0&(1ULL))) { + if (!(D0.a[0]&(1ULL))) { + ++err; + if (err>cut) return INT32_MAX; + } + // fprintf(stderr, "[M::%s::i->%d] err->%d\n", __func__, i, err); + + int32_t site = tn - 1 - thre;///up bound + for (cut = pn - 1, i = 0; site < cut; site++, i++) { + // err += ((VP >> i)&(1ULL)); + err += VP.a[0]&(1ULL); w128_self_rsft_1(VP); + // err -= ((VN >> i)&(1ULL)); + err -= VN.a[0]&(1ULL); w128_self_rsft_1(VN); + } + + if (site == cut && err <= thre) return err; + return INT32_MAX; +} + +inline int32_t ed_band_cal_semi_128bit(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, int32_t *re_err) +{ + (*re_err) = INT32_MAX; + w128_t Peq[5], mm, VP, VN, X, D0, HN, HP; + //Peq[5] = {0}, mm = (Word)1, VP = 0, VN = 0, X, D0, HN, HP; + w128_clear(VP); w128_clear(VN); w128_clear(mm); w128_bit(mm, 0); + w128_clear(Peq[0]); w128_clear(Peq[1]); w128_clear(Peq[2]); w128_clear(Peq[3]); w128_clear(Peq[4]); + int32_t bd = (thre<<1)+1, i, err = 0, i_bd = (thre<<1), last_high = (thre<<1), tn0 = tn - 1; + int32_t cut = thre+last_high; uint8_t c; + + for (i = 0; i < bd; i++) { + w128_self_or(Peq[seq_nt4_table[(uint8_t)pstr[i]]], mm); w128_self_lsft_1(mm); + // Peq[seq_nt4_table[(uint8_t)pstr[i]]] |= mm; mm <<= 1; + } + ///should make Peq[4] = 0 if N is always an error + // Peq[4] = 0; + w128_clear(Peq[4]); + //mm = ((Word)1 << (thre<<1));///for the incoming char/last char + w128_clear(mm); w128_bit(mm, (thre<<1)); + + i = 0; + while (i < tn0) { + // X = Peq[seq_nt4_table[(uint8_t)tstr[i]]] | VN; + c = seq_nt4_table[(uint8_t)tstr[i]]; w128_or(X, Peq[c], VN); + + // D0 = ((VP + (X&VP)) ^ VP) | X; + w128_and(D0, X, VP); + w128_self_add(D0, VP); + w128_self_xor(D0, VP); + w128_self_or(D0, X); + + // HN = VP&D0; + w128_and(HN, VP, D0); + // HP = VN | ~(VP | D0); + w128_or(HP, VP, D0); + w128_self_not(HP); + w128_self_or(HP, VN); + + // X = D0 >> 1; + X = D0; w128_self_rsft_1(X); + // VN = X&HP; + w128_and(VN, X, HP); + // VP = HN | ~(X | HP); + w128_or(VP, X, HP); + w128_self_not(VP); + w128_self_or(VP, HN); + + // if (!(D0&(1ULL))) { + if (!(D0.a[0]&(1ULL))) { + ++err; + if (err>cut) return -1; + } + + // Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1; + w128_self_rsft_1(Peq[0]); w128_self_rsft_1(Peq[1]); + w128_self_rsft_1(Peq[2]); w128_self_rsft_1(Peq[3]); + + ++i; ++i_bd; + // Peq[seq_nt4_table[(uint8_t)pstr[i_bd]]] |= mm; Peq[4] = 0; + c = seq_nt4_table[(uint8_t)pstr[i_bd]]; + if(c < 4) w128_self_or(Peq[c], mm); + } + + // X = Peq[seq_nt4_table[(uint8_t)tstr[i]]] | VN; + c = seq_nt4_table[(uint8_t)tstr[i]]; w128_or(X, Peq[c], VN); + // D0 = ((VP + (X&VP)) ^ VP) | X; + w128_and(D0, X, VP); + w128_self_add(D0, VP); + w128_self_xor(D0, VP); + w128_self_or(D0, X); + // HN = VP&D0; + w128_and(HN, VP, D0); + // HP = VN | ~(VP | D0); + w128_or(HP, VP, D0); + w128_self_not(HP); + w128_self_or(HP, VN); + // X = D0 >> 1; + X = D0; w128_self_rsft_1(X); + // VN = X&HP; + w128_and(VN, X, HP); + // VP = HN | ~(X | HP); + w128_or(VP, X, HP); + w128_self_not(VP); + w128_self_or(VP, HN); + // if (!(D0&(1ULL))) { + if (!(D0.a[0]&(1ULL))) { + ++err; + if (err>cut) return -1; + } + + int32_t site = tn - 1, end = -1;///up bound + ///in most cases, ai = (thre<<1) + int32_t ai = pn - tn, uge = INT32_MAX; + if ((err <= thre) && (err<=(*re_err))) { + *re_err = err; end = site; + } + i = 0; + + while (i < ai) { + // err += ((VP >> i)&(1ULL)); + err += VP.a[0]&(1ULL); w128_self_rsft_1(VP); + // err -= ((VN >> i)&(1ULL)); + err -= VN.a[0]&(1ULL); w128_self_rsft_1(VN); + ++i; + if ((err <= thre) && (err <= (*re_err))) { + *re_err = err; end = site + i; + } + if(i == thre) uge = err; + } + + if((uge<=thre) && (uge == (*re_err))) end = site + thre; + return end; +} /** @@ -277,7 +912,9 @@ inline int Reserve_Banded_BPM (char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err) { (*return_err) = (unsigned int)-1; - + // int32_t rerr, rsite = ed_band_cal_semi_128bit(pattern, p_length, text, t_length, errthold, &rerr); + // if(rsite >= 0) (*return_err) = rerr; + // return rsite; Word Peq[256]; int band_length = (errthold << 1) + 1; @@ -888,7 +1525,6 @@ inline int Reserve_Banded_BPM_PATH return return_site; } - ////four patterns have the same p_length inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *pattern3, char *pattern4, int p_length, char *text, int t_length, int* return_sites, unsigned int* return_sites_error, unsigned short errthold, __m128i* Peq_SSE) @@ -1198,4 +1834,277 @@ inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *p } +#define EAC_M 0 +#define MIS_M 1 +#define MOR_YP 2 +#define MOR_XT 3 + +inline void push_trace(asg16_v *res, uint16_t c, uint32_t len) +{ + uint16_t p; c <<= 14; + while (len >= (0x3fff)) { + p = (c + (0x3fff)); kv_push(uint16_t, *res, p); len -= (0x3fff); + } + if(len) { + p = (c + len); kv_push(uint16_t, *res, p); + } +} + +// void move_trace_gap(uint16_t *trace, int32_t trace_n, int32_t trace_i, +// char *pstr, int32_t pi, char *tstr, int32_t ti, int32_t *err) +// { +// uint16_t c = trace[trace_i]>>14, l = (trace[trace_i]<<2)>>2; +// if(c != 3 && c != 2) return; +// trace_i--; +// if(c == 3) pi--; +// else if(c == 2) ti--; + + +// if() + +// } + +// void adjust_trace(uint16_t *trace, int32_t *trace_n, int32_t *p_beg, int32_t *p_end, int32_t *err, char *pstr, char *tstr) +// { +// if((*err) == 0) return; +// int32_t i, pi, ti; uint16_t c, l; +// for (i = 0; i < (*trace_n) && (trace[i]>>14) == 1; i++) { +// trace[i] <<= 2; trace[i] >>= 2; trace[i] += (((uint16_t)3)<<14); +// l = (trace[i]<<2)>>2; (*p_beg) += l; +// } +// for (i = (*trace_n) - 1; i >= 0 && (trace[i]>>14) == 1; i--) { +// trace[i] <<= 2; trace[i] >>= 2; trace[i] += (((uint16_t)3)<<14); +// l = (trace[i]<<2)>>2; (*p_end) -= l; +// } + +// i = 0; pi = (*p_beg); ti = 0; +// for (i = 0; i < (*trace_n); i++) { +// c = trace[i]>>14; l = (trace[i]<<2)>>2; +// if(c == 0 || c == 1) { +// pi += l; ti += l; +// } else if(c == 2) { +// // move_trace_gap(trace, *trace_n, i, pstr, pi, tstr, ti, err); +// pi += l; +// } else if(c == 3) { +// // move_trace_gap(trace, *trace_n, i, pstr, pi, tstr, ti, err); +// ti += l; +// } +// } +// } + +inline int32_t ungap_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t know_err, int32_t know_end, +int32_t *r_err, int32_t *r_beg, asg16_v *cigar, int32_t *cigar_l) +{ + int32_t cn = cigar->n, pk, tk, e, l; + (*r_err) = (*r_beg) = INT32_MAX; (*cigar_l) = 0; + if(know_err < 0 || know_end < 0) return -1; + if(know_err == 0) { + push_trace(cigar, EAC_M, tn); + (*r_err) = know_err; (*r_beg) = know_end + 1 - tn; (*cigar_l) = cigar->n - cn; + return know_end; + } + + pk = know_end+1-tn; tk = 0; e = 0; + for (l = 0; tk < tn; tk++, pk++) { + if(pstr[pk]!=tstr[tk]) { + e++; if(e > know_err) break; + if(tk > l) push_trace(cigar, EAC_M, tk-l); + push_trace(cigar, MIS_M, 1); l = tk + 1; + } + } + if(tk == tn) { + if(tk > l) push_trace(cigar, EAC_M, tk-l); + (*r_err) = know_err; (*r_beg) = know_end + 1 - tn; (*cigar_l) = cigar->n - cn; + return know_end; + } + + cigar->n = cn; + return -1; +} + + +// ///p_length might be samller than t_length + 2 * errthold +inline int32_t ed_band_cal_semi_trace(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, + int32_t know_err, int32_t know_end, int32_t *r_err, int32_t *r_beg, Word *buf, asg16_v *cigar, int32_t *cigar_l) { + int32_t cn = cigar->n; (*r_err) = (*r_beg) = INT32_MAX; (*cigar_l) = 0; + Word Peq[5] = {0}, mm = (Word)1, VP = 0, VN = 0, X = 0, D0 = 0, HN = 0, HP = 0, i_col, i_col_dux; + int32_t bd = (thre<<1)+1, i, err = 0, i_bd = (thre<<1), last_high = (thre<<1), tn0 = tn - 1; + int32_t cut = thre+last_high; ///kv_resize(uint16_t, *cigar, cigar->n+(uint32_t)know_err+2);//pre-alloc + + if(ungap_trace(pstr, pn, tstr, tn, know_err, know_end, r_err, r_beg, cigar, cigar_l) >= 0) { + return know_end; + } + + for (i = 0; i < bd; i++) { + Peq[seq_nt4_table[(uint8_t)pstr[i]]] |= mm; mm <<= 1; + } + Peq[4] = 0; + ///should make Peq[4] = 0 if N is always an error + i = i_col = 0; mm = ((Word)1 << (thre<<1));///for the incoming char/last char + + while (i < tn0) { + X = Peq[seq_nt4_table[(uint8_t)tstr[i]]] | VN; + + D0 = ((VP + (X&VP)) ^ VP) | X; + + HN = VP&D0; + HP = VN | ~(VP | D0); + + X = D0 >> 1; + VN = X&HP; + VP = HN | ~(X | HP); + + if (!(D0&(1ULL))) { + ++err; + if (err>cut) return -1; + } + + Peq[0] >>= 1; Peq[1] >>= 1; Peq[2] >>= 1; Peq[3] >>= 1; ///Peq[4] >>= 1; + + ++i; ++i_bd; + Peq[seq_nt4_table[(uint8_t)pstr[i_bd]]] |= mm; Peq[4] = 0; + + buf[i_col++] = D0; buf[i_col++] = VP; buf[i_col++] = VN; buf[i_col++] = HP; buf[i_col++] = HN; + } + + X = Peq[seq_nt4_table[(uint8_t)tstr[i]]] | VN; + D0 = ((VP + (X&VP)) ^ VP) | X; + HN = VP&D0; + HP = VN | ~(VP | D0); + X = D0 >> 1; + VN = X&HP; + VP = HN | ~(X | HP); + if (!(D0&(1ULL))) { + ++err; + if (err>cut) return -1; + } + + buf[i_col++] = D0; buf[i_col++] = VP; buf[i_col++] = VN; buf[i_col++] = HP; buf[i_col++] = HN; + + i_col_dux = i_col/tn; + + int32_t site = tn - 1, end = -1;///up bound + ///in most cases, ai = (thre<<1) + int32_t ai = pn - tn, uge = INT32_MAX; + if ((err <= thre) && (err<=(*r_err))) { + *r_err = err; end = site; + } + i = 0; + + while (i < ai) { + err += ((VP >> i)&(1ULL)); err -= ((VN >> i)&(1ULL)); ++i; + if ((err <= thre) && (err <= (*r_err))) { + *r_err = err; end = site + i; + } + if(i == thre) uge = err; + } + if((uge<=thre) && (uge == (*r_err))) end = site + thre; + if ((*r_err) > thre) return end; + + + ///need to correct pn here, since pn might be smaller than tn + 2* thre + pn = tn + (thre<<1); + int32_t beg = end, back_track_site = bd - (pn - end); + + Word v_value, h_value, delta_value, min_value, current_value; + ///Word direction; ///0 is match, 1 is mismatch, 2 is up, 3 is left + Word direction = 0, *ba, pd, pdl; ///0 is match, 1 is mismatch, 2 is up, 3 is left + i = tn; pd = (Word)-1; pdl = 0; + current_value = *r_err; + int32_t low_bound = bd - 1; + + while (i > 0) { + if (current_value == 0) break; + ba = buf + ((i*i_col_dux) - i_col_dux); + delta_value = current_value - ((~(ba[0]>>back_track_site))&(1ULL)); + + if (back_track_site == 0) { + ///HP + h_value = current_value - ((ba[3] >> back_track_site)&(1ULL)); + //HN + h_value = h_value + ((ba[4] >> back_track_site)&1ULL); + min_value = delta_value; direction = 0; + if (h_value < min_value) { + min_value = h_value; + direction = 3; + } + } else if (back_track_site == low_bound) { + v_value = current_value - ((ba[1]>>(back_track_site-1))&(1ULL)); + v_value = v_value + ((ba[2]>>(back_track_site-1))&(1ULL)); + + min_value = delta_value; direction = 0; + if (v_value < min_value) { + min_value = v_value; + direction = 2; + } + } + else { + h_value = current_value-((ba[3]>>back_track_site)&(1ULL)); + h_value = h_value+((ba[4]>>back_track_site)&(1ULL)); + + v_value = current_value - ((ba[1]>>(back_track_site-1))&(1ULL)); + v_value = v_value + ((ba[2]>>(back_track_site-1))&(1ULL)); + + min_value = delta_value; direction = 0; + if (v_value < min_value) { + min_value = v_value; + direction = 2; + } + + if (h_value < min_value) { + min_value = h_value; + direction = 3; + } + } + + + if (direction == 0) { + if (delta_value != current_value) { + direction = 1; + } + i--; beg--; + } + if (direction == 2) {///ru guo xiang shang yi dong, bing bu huan lie + back_track_site--; beg--; + } + else if (direction == 3) {///ru guo xiang zuo yi dong + i--; + back_track_site++; + } + + if(direction != pd) { + if(pdl > 0) push_trace(cigar, pd, pdl); + pd = direction; pdl = 1; + } else { + pdl++; + } + // path[path_length++] = direction; + current_value = min_value; + } + + + if (i > 0) { + direction = 0; beg -= i; + if(direction != pd) { + if(pdl > 0) push_trace(cigar, pd, pdl); + pd = direction; pdl = i; + } else { + pdl += i; + } + } + + if(pdl > 0) push_trace(cigar, pd, pdl); + if (direction != 3) beg++; + + uint16_t *trac = cigar->a + cn, tt; int32_t trac_n = cigar->n - cn; ai = trac_n>>1; + for (i = 0; i < ai; i++) { + tt = trac[i]; trac[i] = trac[trac_n-i-1]; trac[trac_n-i-1] = tt; + } + (*cigar_l) = cigar->n - cn; + + (*r_beg) = beg; (*cigar_l) = cigar->n - cn; + return end; +} + + #endif diff --git a/anchor.cpp b/anchor.cpp index 890ab5c..02c90d1 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -896,7 +896,7 @@ void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uin 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, srt = 0; overlap_region *r; + 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)); @@ -913,13 +913,11 @@ void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uin 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, 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); + r = &(ol->list[ol->length-1]); r->non_homopolymer_errors = m; + 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); } + if(gen_off) gen_fake_cigar(&(r->f_cigar), r, apend_be, cl->list+m, sm); m += sm; } } @@ -933,7 +931,7 @@ void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uin 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; + 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); @@ -959,7 +957,7 @@ 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; @@ -979,6 +977,7 @@ void lchain_gen(Candidates_list* cl, overlap_region_alloc* ol, uint32_t rid, uin } cl->length = m; } + **/ ks_introsort_or_xs(ol->length, ol->list); } diff --git a/inter.cpp b/inter.cpp index 4ccd729..0db5bee 100644 --- a/inter.cpp +++ b/inter.cpp @@ -335,7 +335,7 @@ void hc_gdpchain_destroy(gdpchain_t *b) } void init_mg_opt(mg_idxopt_t *opt, int is_HPC, int k, int w, int hap_n, int max_n_chain, double bw_thres, -double diff_ec_ul, double diff_ec_ul_low, int ec_ul_round) +double diff_ec_ul, double diff_ec_ul_low, double diff_ec_ul_hpc, int ec_ul_round) { opt->k = k; opt->w = w; @@ -364,6 +364,7 @@ double diff_ec_ul, double diff_ec_ul_low, int ec_ul_round) opt->bw_thres = bw_thres; opt->diff_ec_ul = diff_ec_ul; opt->diff_ec_ul_low = diff_ec_ul_low; + opt->diff_ec_ul_hpc = diff_ec_ul_hpc; opt->ec_ul_round = ec_ul_round; } @@ -4936,7 +4937,7 @@ void update_ul_vec_t_ug(const ul_idx_t *uref, ul_vec_t *rch, vec_mg_lchain_t *uc if(sp != (uint32_t)-1) l += ep - sp; l = (int64_t)rch->rlen - l; - // if(ulid == 1756) fprintf(stderr, "-ulid:%ld, l:%ld, rch->rlen:%u\n", ulid, l, rch->rlen); + // fprintf(stderr, "-ulid:%ld, l:%ld, rch->rlen:%u\n", ulid, l, rch->rlen); if(l == 0) { rch->dd = 1; } else if(l < ((int64_t)rch->rlen)*0.001) { @@ -5062,6 +5063,7 @@ int64_t debug_i, int64_t tid, void *km) ll->tk.n = ll->lo.n = 0; kv_ul_ov_t *idx = &(ll->lo); gl_chain_gen(olist, uref, idx, 0, hap, km);///no trans + // fprintf(stderr, "0-[M::%s] idx->n::%lu\n", __func__, (uint64_t)idx->n); if(idx->n == 0) return 0; // fprintf(stderr, "(beg0) [M::%s::tid:%ld] debug_i:%ld, qlen:%ld, # cis:%lu, # trans:%lu\n", __func__, tid, debug_i, qlen, (uint64_t)idx->n, o2); int64_t max_idx, occ = 0, f = 0; @@ -5076,6 +5078,7 @@ int64_t debug_i, int64_t tid, void *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**/); } } + // fprintf(stderr, "1-[M::%s] f::%ld\n", __func__, f); // fprintf(stderr, "(beg1) [M::%s] debug_i:%ld, qlen:%ld\n", __func__, debug_i, qlen); if(!f) { gl_chain_gen(olist, uref, idx, 0, hap, km);///no trans @@ -5099,6 +5102,376 @@ int64_t debug_i, int64_t tid, void *km) return 1; } +int64_t comput_err_partial_cigar(int64_t ol, overlap_region *z, int64_t *rk) +{ + int64_t k = 0, err = 0, e = z->x_pos_s+ol, wn = z->w_list.n; (*rk) = -1; + for (k = 0; k < wn; k++) { + if(z->w_list.a[k].x_start >= e) break; + if(z->w_list.a[k].y_end != -1) { + err += z->w_list.a[k].error; + } + } + k--; + if(k < 0) return 0; + if(z->w_list.a[k].y_end != -1) { + err -= z->w_list.a[k].error; + } + + if((int64_t)z->w_list.a[k].x_end+1 <= e) { + if(z->w_list.a[k].y_end != -1) { + err += z->w_list.a[k].error; + } + } else { + // assert(z->w_list.a[k].x_start < e); + if(z->w_list.a[k].y_end != -1) { + err += (((double)(e-z->w_list.a[k].x_start))/ + ((double)(z->w_list.a[k].x_end+1-z->w_list.a[k].x_start)))*z->w_list.a[k].error; + } + } + (*rk) = k; + return err; +} + +int64_t sum_w_err(window_list *a, int64_t n) +{ + int64_t k, err = 0; + for (k = 0; k < n; k++) { + if(a[k].y_end != -1) err += a[k].error; + } + return err; +} + +int64_t comput_sc_partial_cigar(int64_t sc, int64_t ol, double err_sc_r, overlap_region *z, int64_t *wi, int64_t *werr) +{ + int64_t k = wi?(*wi):0, wn = z->w_list.n, err = werr?(*werr):0, e = z->x_pos_s+ol; + if(ol == 0) return sc; + // int64_t pk, pe; + if((int64_t)(z->x_pos_e + 1 - z->x_pos_s) <= ol) return 0; + if(k == wn) { + k--; + if(z->w_list.a[k].y_end != -1) { + err -= z->w_list.a[k].error; + } + } + if(z->w_list.a[k].x_start >= e) { + if(z->w_list.a[k].y_end != -1) err += z->w_list.a[k].error; + for (;(k>=0) && (z->w_list.a[k].x_start>=e); k--) { + if(z->w_list.a[k].y_end != -1) { + err -= z->w_list.a[k].error; + } + } + } else { + for (;(kw_list.a[k].x_startw_list.a[k].y_end != -1) { + err += z->w_list.a[k].error; + } + } + k--; + } + // pk = (*wi); pe = (*werr); + if(k < 0) { + k = 0; err = 0; + if(wi) (*wi) = k; if(werr) (*werr) = err; + // assert(e <= z->w_list.a[0].x_start); + } else { + if(z->w_list.a[k].y_end != -1) { + err -= z->w_list.a[k].error; + } + if(wi) (*wi) = k; if(werr) (*werr) = err; + + // if(!(err >= 0 && k >= 0 && k < wn && z->w_list.a[k].x_start < e && z->w_list.a[k].x_end + 1 >= e)){ + // fprintf(stderr, "[M::%s] ol::%ld, e::%ld, z::[%u, %u], k::%ld, wn::%ld, w::[%d, %d], err::%ld\n", __func__, + // ol, e, z->x_pos_s, z->x_pos_e, k, wn, z->w_list.a[k].x_start, z->w_list.a[k].x_end, err); + // } + // assert(err >= 0 && k >= 0 && k < wn && z->w_list.a[k].x_start < e && + // (e <= z->w_list.a[k+1].x_start)); + + if((int64_t)z->w_list.a[k].x_end+1 <= e) { + if(z->w_list.a[k].y_end != -1) { + err += z->w_list.a[k].error; + } + } else { + // assert(z->w_list.a[k].x_start < e); + if(z->w_list.a[k].y_end != -1) { + err += (((double)(e-z->w_list.a[k].x_start))/ + ((double)(z->w_list.a[k].x_end+1-z->w_list.a[k].x_start)))*z->w_list.a[k].error; + } + } + } + + // int64_t dbg_k, dbg_e = comput_err_partial_cigar(ol, z, &dbg_k); + // if(err != dbg_e) { + // fprintf(stderr, "[M::%s] ol::%ld, e::%ld, z::[%u, %u], k::%ld, wn::%ld, w::[%d, %d], err::%ld, dbg_e::%ld, dbg_k::%ld, pe::%ld, pk::%ld, sum_pk_err::%ld, sum_k_err::%ld, werr::%ld\n", + // __func__, ol, e, z->x_pos_s, z->x_pos_e, k, wn, z->w_list.a[k].x_start, z->w_list.a[k].x_end, err, dbg_e, dbg_k, pe, pk, + // sum_w_err(z->w_list.a, pk), sum_w_err(z->w_list.a, k), *werr); + // } + // assert(err == dbg_e); + + ol -= (err*err_sc_r); sc -= ol; if(sc <= 0) sc = 1; + return sc; +} + +///mode: 0->ug; 1->read +int64_t ed_dp_c(overlap_region_alloc *o, kv_ul_ov_t *res, ul_ov_t *ex, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, +double diff_ec_ul, int64_t qlen, int64_t max_skip, uint64_t *srt, uint64_t *idx, uint64_t *track, double err_sc, +uint64_t mode, All_reads *ridx, ma_ug_t *ug) +{ + if(res->n == 0) return 0; + uint32_t li_v, lj_v, rev_n; + int64_t mm_ovlp, x, i, j, k, sc, csc, mm_sc, mm_idx, qo, qovl, share, minus_sc, pj, n_skip, wi, werr; + ul_ov_t *li = NULL, *lj = NULL, rev_t; + radix_sort_ul_ov_srt_qe(res->a, res->a + res->n); + for (i = 1, j = 0; i <= (int64_t)res->n; i++) { + if (i == (int64_t)res->n || res->a[i].qe != res->a[j].qe) { + if(i - j > 1) radix_sort_ul_ov_srt_qs(res->a+j, res->a+i); + j = i; + } + } + ///res->a[0].qe: min_qe; res->a[res->n-1].qs: max_qs + if(res->a[0].qe == qlen && res->a[res->n-1].qs == 0) {///all alignments are contained + for (i = 0; i < (int64_t)res->n; ++i) { + li = &(res->a[i]); assert(li->qs == 0 && li->qe == qlen); + csc = (li->qe-li->qs); minus_sc = (o->list[li->qn].non_homopolymer_errors*err_sc); + csc -= minus_sc; if(csc <= 0) csc = 1; mm_sc = csc; mm_idx = -1; + + if(mm_sc > ((int64_t)0x7fffffff)) mm_sc = ((int64_t)0x7fffffff); + track[i] = push_sc_pre(mm_sc, mm_idx); + srt[i] = track[i]>>32; srt[i] <<= 32; srt[i] |= i; + } + } else { + memset(idx, 0, (sizeof((*idx))*res->n)); + for (i = 0; i < (int64_t)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; + if(x < bw) x = bw; + 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 = (li->qe-li->qs); minus_sc = (o->list[li->qn].non_homopolymer_errors*err_sc); + csc -= minus_sc; if(csc <= 0) csc = 1; + mm_sc = csc; mm_idx = -1; n_skip = 0; wi = werr = 0; + for (j = x; j >= 0; --j) { // collect potential destination vertices + lj = &(res->a[j]); 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)) { + qovl = ((MIN(li->qe, lj->qe) > MAX(li->qs, lj->qs))? (MIN(li->qe, lj->qe) - MAX(li->qs, lj->qs)):0); + // fprintf(stderr, "[M::%s::] utg%.6dl->utg%.6dl, icsc::%ld, ierr::%u, ilen::%u, aln::%u, app_sc::%ld\n", + // __func__, (int32_t)li->tn+1, (int32_t)lj->tn+1, csc, o->list[li->qn].non_homopolymer_errors, + // li->qe - li->qs, o->list[li->qn].align_length, comput_sc_partial_cigar(csc, qovl, err_sc, &(o->list[li->qn]), &wi, &werr)); + sc = comput_sc_partial_cigar(csc, qovl, err_sc, &(o->list[li->qn]), &wi, &werr) + + pop_sc(track[j]); + if(sc > mm_sc) { + mm_sc = sc, mm_idx = j; + if (n_skip > 0) --n_skip; + } else if (idx[j] == (uint64_t)i) { + if (++n_skip > max_skip) + break; + } + pj = pop_pre(track[j]); + if(pj >= 0) idx[pj] = i; + } + } + if(mm_sc > ((int64_t)0x7fffffff)) mm_sc = ((int64_t)0x7fffffff); + track[i] = push_sc_pre(mm_sc, mm_idx); + srt[i] = track[i]>>32; srt[i] <<= 32; srt[i] |= i; + } + } + + int64_t n_v, n_u, n_v0; + radix_sort_gfa64(srt, srt+res->n); + for (k = (int64_t)res->n-1, n_v = n_u = 0; k >= 0; --k) { + n_v0 = n_v; + for (i = (uint32_t)srt[k]; i >= 0 && (track[i]&((uint64_t)0x80000000)) == 0;) { + ex[n_v++] = res->a[i]; track[i] |= ((uint64_t)0x80000000); + i = pop_pre(track[i]); + } + if(n_v0 == n_v) continue; + sc = (i<0?(pop_sc(srt[k])):(pop_sc(srt[k])-pop_sc(track[i]))); + if(sc < 0) { + n_v = n_v0; + continue; + } + idx[n_u++] = ((uint64_t)sc<<32)|(n_v-n_v0); + } + + for (k = 0, n_v = n_v0 = 0; k < n_u; k++) { + n_v0 = n_v; n_v += (uint32_t)idx[k]; + res->a[k].qn = idx[k]>>32;//score + res->a[k].ts = n_v0; res->a[k].te = n_v;///idx + + rev_n = ((uint32_t)idx[k])>>1; + ///we need to consider contained reads; so determining qs is not such easy + res->a[k].qs = (uint32_t)-1; res->a[k].qe = ex[n_v0].qe; + for (i = 0; i < rev_n; i++) { + rev_t = ex[n_v0+i]; ex[n_v0+i] = ex[n_v-i-1]; ex[n_v-i-1] = rev_t; + if(res->a[k].qs > ex[n_v0+i].qs) res->a[k].qs = ex[n_v0+i].qs; + if(res->a[k].qs > ex[n_v-i-1].qs) res->a[k].qs = ex[n_v-i-1].qs; + } + if(((uint32_t)idx[k])&1) { + if(res->a[k].qs > ex[n_v0+i].qs) res->a[k].qs = ex[n_v0+i].qs; + } + } + + res->n = n_u; + radix_sort_ul_ov_srt_qn(res->a, res->a + res->n);//sort by score + // fprintf(stderr, "---[M::%s] n_u:%ld, n_v:%ld\n", __func__, n_u, n_v); + return n_v; +} + + +void set_w_e(overlap_region *z, uint64_t *w_idx, int64_t wl, int64_t ql) +{ + int64_t wid, k, wn = z->w_list.n, ws, we; + for (k = 0; k < wn; k++) { + wid = z->w_list.a[k].x_start/wl; + ws = wid*wl; we = ws+wl; if(we > ql) we = ql; we--; + // fprintf(stderr, "[M::%s] ws::%ld, we::%ld, xs::%d, xe::%d, err::%d\n", __func__, + // ws, we, z->w_list.a[k].x_start, z->w_list.a[k].x_end, z->w_list.a[k].error); + if(ws == z->w_list.a[k].x_start && we == z->w_list.a[k].x_end && z->w_list.a[k].y_end != -1) { + if((w_idx[wid] == (uint64_t)-1) || (w_idx[wid] < (uint64_t)z->w_list.a[k].error)) { + w_idx[wid] = z->w_list.a[k].error; + } + } + } +} + +uint32_t ck_w_err(overlap_region *z, uint64_t *w_idx, int64_t wl, int64_t ql) +{ + int64_t wid, k, wn = z->w_list.n, ws, we, ol, e[2]; + ol = e[0] = e[1] = 0; + for (k = 0; k < wn; k++) { + wid = z->w_list.a[k].x_start/wl; + if(w_idx[wid] == (uint64_t)-1) continue; + ws = wid*wl; we = ws+wl; if(we > ql) we = ql; we--; + if(ws == z->w_list.a[k].x_start && we == z->w_list.a[k].x_end) { + ol += we+1-ws; e[0] += w_idx[wid]; + if(z->w_list.a[k].y_end != -1) e[1] += z->w_list.a[k].error; + else e[1] += THRESHOLD_MAX_SIZE + 1; + } + } + // 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]+(ol*0.01))) return 0; + if(e[1] > (e[0]+(e[0]*0.01))) return 0; + } + // if((e[1] > (e[0]+16)) && (e[1] > (e[0]+(ol*0.01)))) return 0; + return 1; +} + +int64_t filter_sec(overlap_region_alloc *ol, ul_ov_t *idx, int64_t idx_n, ul_ov_t *a, uint64_t *w_idx, uint64_t nw, uint64_t wl, uint64_t ql) +{ + if(idx_n <= 0) return 1; + int64_t on = ol->length, k, z, on_contain = 0, max_i = -1, max_k = -1, alt_occ = 0; overlap_region t; + memset(w_idx, -1, nw*sizeof((*w_idx))); + for (k = 0; k < on; k++) ol->list[k].is_match = 0; + for (k = 0; k < idx_n; k++) { + // fprintf(stderr, "[M::%s::pri_chain[%ld]] q_coord::[%u, %u), occ::%u\n", + // __func__, k, idx[k].qs, idx[k].qe, idx[k].te-idx[k].ts); + for (z = idx[k].ts; z < idx[k].te; z++) { + ol->list[a[z].qn].is_match = 2; + set_w_e(&(ol->list[a[z].qn]), w_idx, wl, ql); + // fprintf(stderr, "[M::%s::utg%.6dl]\n", __func__, (int32_t)a[z].tn+1); + } + on_contain += (((idx[k].te-idx[k].ts)==1)?1:0); + } + if(on_contain == idx_n) {///each primary chain only has one alignment + on_contain = 0; + } else { + on_contain = -on-1;///in this case, on_contain == z is always wrong + } + max_i = a[idx[idx_n-1].ts].qn; + for (k = z = 0; k < on; k++) { + if(!ol->list[k].is_match) ol->list[k].is_match = ck_w_err(&(ol->list[k]), w_idx, wl, ql); + if(!ol->list[k].is_match) continue; + if(z != k) { + t = ol->list[k]; + ol->list[k] = ol->list[z]; + ol->list[z] = t; + } + if(ol->list[z].x_pos_s == 0 && ol->list[z].x_pos_e == ql - 1) { + on_contain++; + if(max_i == k) max_k = z; + } + if(ol->list[z].is_match == 1) alt_occ++; + else ol->list[z].is_match = 1; + z++; + } + ol->length = z; + // fprintf(stderr, "+[M::%s] oln::%ld\n", __func__, ol->length); + if(on_contain == z) {///do not contribute to phase + k = max_k; z = 0; + if(z != k) { + t = ol->list[k]; + ol->list[k] = ol->list[z]; + ol->list[z] = t; + } + ol->length = 1; + } + // fprintf(stderr, "-[M::%s] oln::%ld\n", __func__, ol->length); + if(alt_occ == 0 || ol->length == 1) return 1;//if all alignments are primary or there is only one alignment + return 0; +} + +int64_t gl_chain_flter(overlap_region_alloc* olist, Correct_dumy* dumy, st_mt_t *sps, glchain_t *ll, const ul_idx_t *uref, double diff_ec_ul, int64_t wl, int64_t ql, const ug_opt_t *uopt, uint32_t *need_phase) +{ + (*need_phase) = 1; + uint64_t k, nw; ul_ov_t *p, *m; int64_t occ, i, ovlp, idx_n; + ll->tk.n = ll->lo.n = 0; + kv_ul_ov_t *idx = &(ll->lo); idx->n = 0; + kv_resize(ul_ov_t, *idx, olist->length); + for (k = 0; k < olist->length; k++) { + 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; + p->sec = olist->list[k].non_homopolymer_errors; + if(p->rev) { + p->ts = uref->ug->u.a[p->tn].len - (olist->list[k].y_pos_e+1); + p->te = uref->ug->u.a[p->tn].len - olist->list[k].y_pos_s; + } else { + p->ts = olist->list[k].y_pos_s; + p->te = olist->list[k].y_pos_e+1; + } + } + if(idx->n == 0) return 0; + kv_resize(uint64_t, ll->srt.a, idx->n); + kv_resize(uint64_t, *sps, idx->n); + kv_resize(ul_ov_t, ll->tk, idx->n); + + occ = ed_dp_c(olist, idx, ll->tk.a, uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, ql, 75, dumy->overlapID, ll->srt.a.a, sps->a, 1.25, 0, NULL, uref->ug); + if((!occ) || (!idx->n)) return 0; + idx_n = idx->n; p = &(idx->a[idx_n-1]); + // fprintf(stderr, "[M::%s] qs::%u, qe::%u, ql::%ld, occ::%u\n", __func__, p->qs, p->qe, ql, p->te - p->ts); + if(p->qe-p->qs <= (ql*0.25)) return 0;///primary chain is too short + i = idx_n-1; occ = p->te - p->ts; + if(p->qe-p->qs < ql && idx_n > 1) { + for (occ = 0; i >= 0; i--) { + p = &(idx->a[i]); + for (k = i + 1; k < idx->n; k++) { + m = &(idx->a[k]); + ovlp = ((MIN(m->qe, p->qe) > MAX(m->qs, p->qs))? (MIN(m->qe, p->qe) - MAX(m->qs, p->qs)):0); + if(((ovlp > ((m->qe-m->qs)*0.005)) || (ovlp > ((p->qe-p->qs)*0.015))) && ovlp > 32) break; + if((ovlp == (m->qe-m->qs)) || (ovlp == (p->qe-p->qs))) break; + } + if(k < idx->n) break; + occ += p->te - p->ts; + } + i++; + } + + // fprintf(stderr, "[M::%s] i::%ld, idx_n::%ld\n", __func__, i, ((int64_t)idx->n)); + if(occ == (int64_t)olist->length) return 1; + // if(i >= ((int64_t)idx->n)) return 0; + nw = get_num_wins(0, ql, wl); kv_resize(uint64_t, ll->srt.a, (uint64_t)nw); + if(filter_sec(olist, idx->a+i, idx->n-i, ll->tk.a, ll->srt.a.a, nw, wl, ql)) { + (*need_phase) = 0; + } + return 1; +} + uint64_t kv_ul_ov_t_statistics(kv_ul_ov_t *olist, uint64_t qn, int64_t *occ) { int64_t k, l = 0; @@ -5216,9 +5589,9 @@ 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; + uint32_t high_occ = 2, phase = 1; // uint64_t align = 0; - int fully_cov, abnormal; + // if(UL_INF.a[s->id+i].rlen != s->len[i]) { // fprintf(stderr, "[M::%s] rid:%ld, s->len:%lu, UL_INF->rlen:%u\n", __func__, s->id+i, s->len[i], UL_INF.a[s->id+i].rlen); // } @@ -5227,7 +5600,7 @@ 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!=47) return; + // if((s->id+i!=43) /**&& (s->id+i!=44) && (s->id+i!=948)**/) 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); @@ -5243,9 +5616,22 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // return; // b->num_correct_base += overlap_statistics(&b->olist, NULL, 0); - b->self_read.seq = s->seq[i]; b->self_read.length = s->len[i]; b->self_read.size = 0; - correct_ul_overlap(&b->olist, s->uu, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon, - &b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul, winLen, NULL); + // int fully_cov, abnormal; + // b->self_read.seq = s->seq[i]; b->self_read.length = s->len[i]; b->self_read.size = 0; + // correct_ul_overlap(&b->olist, s->uu, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon, + // &b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul, winLen, NULL); + // memset(&b->self_read, 0, sizeof(b->self_read)); + + ul_lalign(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, + &b->correct, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, 1, NULL); + + gl_chain_flter(&b->olist, &b->correct, &(s->sps[tid]), bl, s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, &phase); + + if(phase) { + ul_lalign(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read, + &b->correct, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, 0, NULL); + } + // exit(1); // uint64_t k; // for (k = 0; k < b->olist.length; k++) { @@ -5260,7 +5646,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // b->num_read_base += b->self_read.length; // b->num_correct_base += b->correct.corrected_base; // b->num_recorrect_base += b->round2.dumy.corrected_base; - memset(&b->self_read, 0, sizeof(b->self_read)); + if(UL_INF.a[s->id+i].dd) { free(s->seq[i]); s->seq[i] = NULL; b->num_correct_base++; } @@ -9150,6 +9536,10 @@ int rescall_ul_pipeline(uldat_t* sl, const enzyme *fn) // fprintf(stderr, "[M::%s::] ==> # bases: %lu; # corrected bases: %lu; # recorrected bases: %lu\n", // __func__, sl->num_bases, sl->num_corrected_bases, sl->num_recorrected_bases); // gen_ul_vec_rid_t(&UL_INF); + // for (i = 0; i < UL_INF.n; i++) { + // fprintf(stderr, "[M::%s] rid::%d, dd::%u\n", __func__, i, UL_INF.a[i].dd); + // } + return 1; } @@ -9666,7 +10056,7 @@ void ul_resolve(ma_ug_t *ug, const asg_t *rg, const ug_opt_t *uopt, int hap_n) { fprintf(stderr, "[M::%s::] ==> UL\n", __func__); mg_idxopt_t opt; - init_mg_opt(&opt, 0, 19, 10, hap_n, 0, 0, 0.05, asm_opt.ul_error_rate_low, asm_opt.ul_ec_round); + init_mg_opt(&opt, 0, 19, 10, hap_n, 0, 0, 0.05, asm_opt.ul_error_rate_low, asm_opt.ul_error_rate_hpc, asm_opt.ul_ec_round); int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name, NULL) : 0); if(exist == 0) uidx_build(ug, &opt); if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name, NULL); @@ -10655,7 +11045,7 @@ void ul_load(const ug_opt_t *uopt) int32_t cutoff; init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); cutoff = asm_opt.max_n_chain; - init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate, asm_opt.ul_error_rate_low, asm_opt.ul_ec_round); + init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate, asm_opt.ul_error_rate_low, asm_opt.ul_error_rate_hpc, asm_opt.ul_ec_round); init_uldat_t(&sl, NULL, NULL, &opt, CHUNK_SIZE, asm_opt.thread_num, uopt, NULL); if(!load_all_ul_t(&UL_INF, asm_opt.output_file_name, &R_INF, NULL)) { @@ -10685,7 +11075,7 @@ uint64_t ul_refine_alignment(const ug_opt_t *uopt, asg_t *sg) mg_idxopt_t opt; uldat_t sl; int32_t cutoff; init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); cutoff = asm_opt.max_n_chain; - init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate, asm_opt.ul_error_rate_low, asm_opt.ul_ec_round); + init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate, asm_opt.ul_error_rate_low, asm_opt.ul_error_rate_hpc, asm_opt.ul_ec_round); ul_idx_t *uu = gen_ul_idx_t(uopt, sg, 0, 0);///record contained reads; is_el = is_del = 0 init_uldat_t(&sl, NULL, NULL, &opt, CHUNK_SIZE, asm_opt.thread_num, uopt, uu); sl.rg = sg; if(work_ul_gchains(&sl)) { @@ -10738,7 +11128,7 @@ ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_c init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); cutoff = REA_ALIGN_CUTOFF; - init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate, asm_opt.ul_error_rate_low, asm_opt.ul_ec_round); + init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate, asm_opt.ul_error_rate_low, asm_opt.ul_error_rate_hpc, asm_opt.ul_ec_round); init_uldat_t(&sl, NULL, NULL, &opt, CHUNK_SIZE, asm_opt.thread_num, uopt, NULL); ma_ug_t *ug = gen_polished_ug(uopt, sg); // dd_ug(sg, ug, uopt->coverage_cut, uopt->sources, uopt->ruIndex, "UL.sa"); diff --git a/inter.h b/inter.h index 77ac3eb..8adc7c2 100644 --- a/inter.h +++ b/inter.h @@ -26,7 +26,7 @@ typedef struct { int min_gc_cnt, min_gc_score, sub_diff, best_n; float chn_pen_gap, mask_level, pri_ratio; ///base-alignment - double bw_thres, diff_ec_ul, diff_ec_ul_low; int max_n_chain, ec_ul_round; + double bw_thres, diff_ec_ul, diff_ec_ul_low, diff_ec_ul_hpc; int max_n_chain, ec_ul_round; } mg_idxopt_t; struct mg_tbuf_s { diff --git a/main.cpp b/main.cpp index 0518279..80a092a 100644 --- a/main.cpp +++ b/main.cpp @@ -12,6 +12,15 @@ int main(int argc, char *argv[]) yak_reset_realtime(); init_opt(&asm_opt); if (!CommandLine_process(argc, argv, &asm_opt)) return 0; + // fprintf(stderr, "[M::%s::] ed_global::%d, ed_global_128bit::%d\n", __func__, + // ed_band_cal_global((char *)"ACT", 3, (char *)"AAT", 3, 1), + // ed_band_cal_global_128bit((char *)"ACT", 3, (char *)"AAT", 3, 1)); + + // fprintf(stderr, "[M::%s::] ed_global::%d, ed_global_128bit::%d\n", __func__, + // ed_band_cal_global((char*)"ACTTTTTT", 8, (char*)"AATTTT", 6, 3), + // ed_band_cal_global_128bit((char*)"ACTTTTTT", 8, (char*)"AATTTT", 6, 3)); + // exit(1); + ret = ha_assemble(); destory_opt(&asm_opt); fprintf(stderr, "[M::%s] Version: %s\n", __func__, HA_VERSION);