diff --git a/CommandLines.h b/CommandLines.h index e8166af..faa534e 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.16.4-r4000" +#define HA_VERSION "0.16.4-r402" #define VERBOSE 0 diff --git a/Correct.cpp b/Correct.cpp index 375773c..1ea4c1f 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -487,6 +487,7 @@ char* r_string) threshold = x_len * asm_opt.max_ov_diff_ec; /****************************may have bugs********************************/ threshold = Adjust_Threshold(threshold, x_len); + if(threshold > THRESHOLD_MAX_SIZE) threshold = THRESHOLD_MAX_SIZE; /****************************may have bugs********************************/ ///offset of y @@ -1038,23 +1039,25 @@ inline void get_win_se_by_normalize_xs(overlap_region *z, int64_t norm_w_s, int6 } } -inline int64_t get_init_err_thres(int64_t len, double e_rate, int64_t block_s) +inline int64_t get_init_err_thres(int64_t len, double e_rate, int64_t block_s, int64_t block_err) { - if(len >= block_s) return THRESHOLD; + if(len >= block_s) return block_err; int64_t thres = len * e_rate; - return Adjust_Threshold(thres, len); + thres = Adjust_Threshold(thres, len); + if(thres > THRESHOLD_MAX_SIZE) thres = THRESHOLD_MAX_SIZE; + return thres; } -uint32_t get_init_paras(All_reads* R_INF, overlap_region *z, int64_t x_s, int64_t x_e, double e_rate, int64_t block_s, +uint32_t get_init_paras(All_reads* rref, const ul_idx_t *uref, overlap_region *z, int64_t x_s, int64_t x_e, double e_rate, int64_t block_s, int64_t *r_ys, int64_t *r_ex_beg, int64_t *r_ex_end, int64_t *r_err_thre) { int e, ex_beg, ex_end; long long y_s, o_len, Window_Len; - e = get_init_err_thres(x_e+1-x_s, e_rate, block_s); + e = get_init_err_thres(x_e+1-x_s, e_rate, block_s, rref?THRESHOLD:THRESHOLD_MAX_SIZE); y_s = (x_s-z->x_pos_s) + z->y_pos_s; y_s += y_start_offset(x_s, &(z->f_cigar)); Window_Len = (x_e+1-x_s) + (e<<1); - if(!determine_overlap_region(e, y_s, z->y_id, Window_Len, Get_READ_LENGTH((*R_INF), z->y_id), + if(!determine_overlap_region(e, y_s, z->y_id, Window_Len, (rref?(Get_READ_LENGTH((*rref), z->y_id)):(uref->ug->u.a[z->y_id].len)), &ex_beg, &ex_end, &y_s, &o_len)) { return 0; } @@ -1072,13 +1075,14 @@ int64_t check_coverage_gap(const kvec_t_u64_warp* v_idx, uint64_t w_s, uint64_t return 0; } -inline double non_trim_error_rate(overlap_region *z, All_reads* R_INF, const ul_idx_t *uref, const kvec_t_u64_warp* v_idx, Correct_dumy* dumy, UC_Read* g_read, double e_rate, int64_t block_s) +inline double non_trim_error_rate(overlap_region *z, All_reads* rref, const ul_idx_t *uref, const kvec_t_u64_warp* v_idx, Correct_dumy* dumy, UC_Read* g_read, double e_rate, int64_t block_s) { int64_t nw, aw = z->w_list.n, k, m, w_id, wn_id, w_s, w_e, idx_e, tErr = 0, tLen = 0, y_s, ex_beg, ex_end, err_thre, p_err_thre; int64_t x_len, Window_Len, y_beg_left, y_beg_right; unsigned int r_error_left, r_error_right; int32_t r_x_end_left, r_y_end_left, aligned_xLen_left, r_x_end_right, r_y_end_right, aligned_xLen_right; nw = get_num_wins(z->x_pos_s, z->x_pos_e+1, block_s); assert(nw >= aw && aw > 0); + for (k = aw-1, idx_e = nw; k >= 0; k--) { w_id = get_win_id_by_e(z, z->w_list.a[k].x_end, block_s, &w_s); assert(w_s == z->w_list.a[k].x_start && w_id < idx_e && k <= w_id); @@ -1095,14 +1099,19 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* R_INF, const ul_ if(check_coverage_gap(v_idx, w_s, w_e, block_s)) { tErr += THRESHOLD_MAX_SIZE; continue; - } + } } - if(!get_init_paras(R_INF, z, w_s, w_e, e_rate, block_s, &y_s, &ex_beg, &ex_end, &err_thre)) { + 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; continue; } p_err_thre = err_thre; - err_thre = double_error_threshold(err_thre, x_len); + + if(rref) { + err_thre = double_error_threshold(err_thre, x_len); + } else { + err_thre = double_ul_error_threshold(err_thre, x_len); + } Window_Len = x_len + (err_thre << 1); r_error_left = r_error_right = 0; aligned_xLen_left = aligned_xLen_right = 0; @@ -1126,8 +1135,8 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* R_INF, const ul_ if(y_beg_right == -1 && y_beg_left != -1) y_beg_right = y_beg_left; if(y_beg_left != -1) { - if(R_INF) { - verify_sub_window(R_INF, dumy, g_read, w_s, x_len, y_beg_left, Window_Len, + if(rref) { + verify_sub_window(rref, dumy, g_read, w_s, x_len, y_beg_left, Window_Len, z->y_id, z->y_pos_strand, err_thre, 0, &r_error_left, &r_y_end_left, &r_x_end_left, &aligned_xLen_left); } else { verify_ul_sub_window(uref, dumy, g_read, w_s, x_len, y_beg_left, Window_Len, @@ -1136,8 +1145,8 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* R_INF, const ul_ } if(y_beg_right != -1) { - if(R_INF) { - verify_sub_window(R_INF, dumy, g_read, w_s, x_len, y_beg_right, Window_Len, z->y_id, z->y_pos_strand, + if(rref) { + verify_sub_window(rref, dumy, g_read, w_s, x_len, y_beg_right, Window_Len, z->y_id, z->y_pos_strand, err_thre, 1, &r_error_right, &r_y_end_right, &r_x_end_right, &aligned_xLen_right); } else { verify_ul_sub_window(uref, dumy, g_read, w_s, x_len, y_beg_right, Window_Len, z->y_id, z->y_pos_strand, @@ -1172,12 +1181,29 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* R_INF, const ul_ w_s = w_e + 1; wn_id = get_win_id_by_s(z, w_s, block_s, &w_e); assert(wn_id == m); x_len = w_e + 1 - w_s; tLen += x_len; - if(!get_init_paras(R_INF, z, w_s, w_e, e_rate, block_s, &y_s, &ex_beg, &ex_end, &err_thre)) { + ///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)) { + tErr += THRESHOLD_MAX_SIZE; + continue; + } + // else { + // if(z->y_id == 575) { + // fprintf(stderr, "---[M::%s::] z::y_id->%u, w_s->%ld, w_e->%ld\n", __func__, z->y_id, w_s, w_e); + // } + // } + } + 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; continue; } p_err_thre = err_thre; - err_thre = double_error_threshold(err_thre, x_len); + if(rref) { + err_thre = double_error_threshold(err_thre, x_len); + } else { + err_thre = double_ul_error_threshold(err_thre, x_len); + } Window_Len = x_len + (err_thre << 1); r_error_left = r_error_right = 0; aligned_xLen_left = aligned_xLen_right = 0; @@ -1202,8 +1228,8 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* R_INF, const ul_ if(y_beg_right == -1 && y_beg_left != -1) y_beg_right = y_beg_left; if(y_beg_left != -1) { - if(R_INF) { - verify_sub_window(R_INF, dumy, g_read, w_s, x_len, y_beg_left, Window_Len, + if(rref) { + verify_sub_window(rref, dumy, g_read, w_s, x_len, y_beg_left, Window_Len, z->y_id, z->y_pos_strand, err_thre, 0, &r_error_left, &r_y_end_left, &r_x_end_left, &aligned_xLen_left); } else { verify_ul_sub_window(uref, dumy, g_read, w_s, x_len, y_beg_left, Window_Len, @@ -1212,8 +1238,8 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* R_INF, const ul_ } if(y_beg_right != -1) { - if(R_INF) { - verify_sub_window(R_INF, dumy, g_read, w_s, x_len, y_beg_right, Window_Len, z->y_id, z->y_pos_strand, + if(rref) { + verify_sub_window(rref, dumy, g_read, w_s, x_len, y_beg_right, Window_Len, z->y_id, z->y_pos_strand, err_thre, 1, &r_error_right, &r_y_end_right, &r_x_end_right, &aligned_xLen_right); } else { verify_ul_sub_window(uref, dumy, g_read, w_s, x_len, y_beg_right, Window_Len, z->y_id, z->y_pos_strand, @@ -2853,8 +2879,14 @@ int64_t get_adjust_winid(overlap_region *z, int64_t win_beg, int64_t win_len) win_id = (win_beg-((z->x_pos_s/win_len)*win_len))/win_len; if((uint64_t)win_id < z->w_list.n && z->w_list.a[win_id].x_start == win_beg) return win_id; if(z->w_list.n == 0) return -1; - assert(z->w_list.a[win_id].x_start > win_beg); - for (k = win_id - 1; k >= 0; k++) { + // if(z->w_list.a[win_id].x_start <= win_beg) { + // fprintf(stderr, "z->w_list.n::%u, z->w_list.a[%ld].x_start::%d, win_beg::%ld\n", + // (uint32_t)z->w_list.n, win_id, z->w_list.a[win_id].x_start, win_beg); + // } + if((uint64_t)win_id > z->w_list.n) win_id = z->w_list.n; + // assert((z->w_list.a[win_id].x_start > win_beg); + for (k = win_id - 1; k >= 0; k--) { + // if(k < 0 || k >= (int64_t)z->w_list.n) fprintf(stderr, "win_id::%ld, k::%ld, z->w_list.n::%ld\n", win_id, k, (int64_t)z->w_list.n); if(z->w_list.a[k].x_start == win_beg) return k; if(z->w_list.a[k].x_start < win_beg) return -1; } @@ -2896,17 +2928,7 @@ void set_herror_win(overlap_region_alloc* ovlp, Correct_dumy* du, kvec_t_u64_war if(ovlp->list[cID].is_match!=3 && ovlp->list[cID].is_match!=4) continue; ovlp->list[cID].is_match = 4; ovlp->list[cID].align_length += window_end + 1 - window_start; - // ovlp->list[cID].non_homopolymer_errors++; mm++; - // w_list_id = (window_start - (ovlp->list[cID].x_pos_s / w_inf.window_length)* - // w_inf.window_length)/w_inf.window_length; - // if (ovlp->list[cID].w_list[w_list_id].y_end == -1) { - // ovlp->list[cID].w_list[w_list_id].y_end = -2; - // ovlp->list[cID].w_list[w_list_id].error = THRESHOLD_MAX_SIZE; - // ovlp->list[cID].align_length += ovlp->list[cID].w_list[w_list_id].x_end + 1 - // - ovlp->list[cID].w_list[w_list_id].x_start; - // ovlp->list[cID].is_match = 4; - // } } if(mm > 0) { kv_push(uint64_t, v_idx->a, (((uint64_t)window_start)<<32)|((uint64_t)window_end)); @@ -2918,7 +2940,8 @@ void set_herror_win(overlap_region_alloc* ovlp, Correct_dumy* du, kvec_t_u64_war for (i = du->size-du->lengthNT, mLen = du->size-du->lengthNT, fc = 0; i < (int64_t)du->size; i++) { cID = (uint32_t)du->overlapID[i]; if(ovlp->list[cID].is_match!=3 && ovlp->list[cID].is_match!=4) continue; - w_list_id = get_adjust_winid(&(ovlp->list[cID]), window_start, w_inf.window_length); + get_win_se_by_normalize_xs(&(ovlp->list[cID]), window_start, blockLen, &ws, &we); + w_list_id = get_adjust_winid(&(ovlp->list[cID]), ws, blockLen); if (w_list_id >= 0) {///matched cID = w_list_id; cID <<= 32; cID += (uint32_t)du->overlapID[i]; du->overlapID[i] = cID; if(mLen != i) { @@ -2950,21 +2973,15 @@ void set_herror_win(overlap_region_alloc* ovlp, Correct_dumy* du, kvec_t_u64_war if(k >= mLen) {///no matched window can cover the unmatched window ovlp->list[cID].is_match = 4; ovlp->list[cID].align_length += we + 1 - ws; - // ovlp->list[cID].non_homopolymer_errors++; kv_push(uint64_t, v_idx->a, (((uint64_t)ws)<<32)|((uint64_t)we)); v_idx->a.a[idx_i-1]++; - // ovlp->list[cID].w_list[w_list_id].y_end = -2; - // ovlp->list[cID].w_list[w_list_id].error = THRESHOLD_MAX_SIZE; - // ovlp->list[cID].align_length += ovlp->list[cID].w_list[w_list_id].x_end + 1 - // - ovlp->list[cID].w_list[w_list_id].x_start; - // ovlp->list[cID].is_match = 4; } } } } -inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_reads* R_INF, const ul_idx_t *uref, +inline void recalcate_window_advance(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) { long long j, k, i; @@ -2993,9 +3010,10 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea overlap_region *z; window_list *p = NULL; + overlap_list->mapped_overlaps_length = 0; for (j = 0; j < (long long)overlap_list->length; j++) { - z = &(overlap_list->list[j]); z->is_match = 0; + z = &(overlap_list->list[j]); z->is_match = 0; is_srt = 1; if(z->w_list.n == 0) continue;///no alignment nw = get_num_wins(z->x_pos_s, z->x_pos_e+1, block_s); a_nw = z->w_list.n; kv_resize(uint64_t, v_idx->a, (uint64_t)nw); memset(v_idx->a.a, -1, sizeof((*v_idx->a.a))*nw); w_idx = v_idx->a.a; @@ -3010,7 +3028,7 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea // } y_id = z->y_id; y_strand = z->y_pos_strand; - y_readLen = (R_INF?(Get_READ_LENGTH((*R_INF), y_id)):(uref->ug->u.a[y_id].len)); + y_readLen = (rref?(Get_READ_LENGTH((*rref), y_id)):(uref->ug->u.a[y_id].len)); 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); // if(z->w_list.a[i].x_end != w_e) { @@ -3031,22 +3049,22 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea ///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(R_INF) { - threshold = double_error_threshold(get_init_err_thres(x_len, e_rate, block_s), x_len); + if(rref) { + threshold = double_error_threshold(get_init_err_thres(x_len, e_rate, block_s, THRESHOLD), x_len); } else { - threshold = double_ul_error_threshold(get_init_err_thres(x_len, e_rate, block_s), x_len); + threshold = double_ul_error_threshold(get_init_err_thres(x_len, e_rate, block_s, THRESHOLD_MAX_SIZE), x_len); } Window_Len = x_len + (threshold << 1); - if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, (R_INF?(Get_READ_LENGTH((*R_INF), y_id)):(uref->ug->u.a[y_id].len)), + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, (rref?(Get_READ_LENGTH((*rref), y_id)):(uref->ug->u.a[y_id].len)), &extra_begin, &extra_end, &y_start, &o_len)) { break; } if(o_len + threshold < x_len) break; - if(R_INF) { - fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, R_INF, y_id, extra_begin, extra_end); + if(rref) { + fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, rref, y_id, extra_begin, extra_end); } else { fill_subregion_ul(dumy->overlap_region, y_start, o_len, y_strand, uref, y_id, extra_begin, extra_end); } @@ -3067,6 +3085,8 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea p->cidx = p->clen = 0; z->align_length += x_len; w_idx[k] = z->w_list.n - 1; + + if(is_srt && z->w_list.n > 1 && p->x_start < z->w_list.a[z->w_list.n-2].x_start) is_srt = 0; } else { break; @@ -3097,8 +3117,8 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea ///y_start is the real y_start y_start = p->y_start; extra_begin = p->extra_begin; extra_end = p->extra_end; o_len = Window_Len - extra_end - extra_begin; - if(R_INF) { - fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, R_INF, y_id, extra_begin, extra_end); + if(rref) { + fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, rref, y_id, extra_begin, extra_end); } else { fill_subregion_ul(dumy->overlap_region, y_start, o_len, y_strand, uref, y_id, extra_begin, extra_end); } @@ -3114,9 +3134,9 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea ///this condition is always wrong ///in best case, real_y_start = threshold, end_site = Window_Len - threshold - 1 if (end_site == Window_Len - 1 || real_y_start == 0) { - if(R_INF) { + if(rref) { if(fix_boundary(x_string, x_len, threshold, y_start, real_y_start, - end_site, extra_begin, extra_end, y_id, Window_Len, R_INF, dumy, + end_site, extra_begin, extra_end, y_id, Window_Len, rref, dumy, y_strand, error, &y_start, &real_y_start, &end_site, &extra_begin, &extra_end, &error)) { p->error = error; p->extra_begin = extra_begin; p->extra_end = extra_end; @@ -3158,25 +3178,25 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea ///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(R_INF) { - threshold = double_error_threshold(get_init_err_thres(x_len, e_rate, block_s), x_len); + if(rref) { + threshold = double_error_threshold(get_init_err_thres(x_len, e_rate, block_s, THRESHOLD), x_len); } else { - threshold = double_ul_error_threshold(get_init_err_thres(x_len, e_rate, block_s), x_len); + threshold = double_ul_error_threshold(get_init_err_thres(x_len, e_rate, block_s, THRESHOLD_MAX_SIZE), x_len); } Window_Len = x_len + (threshold << 1); if(total_y_end <= 0) break; ///y_start might be less than 0 y_start = total_y_end - x_len + 1; - if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, (R_INF?(Get_READ_LENGTH((*R_INF), y_id)):(uref->ug->u.a[y_id].len)), + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, (rref?(Get_READ_LENGTH((*rref), y_id)):(uref->ug->u.a[y_id].len)), &extra_begin, &extra_end, &y_start, &o_len)) { break; } if(o_len + threshold < x_len) break; - if(R_INF) { - fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, R_INF, y_id, extra_begin, extra_end); + if(rref) { + fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, rref, y_id, extra_begin, extra_end); } else { fill_subregion_ul(dumy->overlap_region, y_start, o_len, y_strand, uref, y_id, extra_begin, extra_end); } @@ -3191,9 +3211,9 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea ///this condition is always wrong ///in best case, real_y_start = threshold, end_site = Window_Len - threshold - 1 if (end_site == Window_Len - 1 || real_y_start == 0) { - if(R_INF) { + if(rref) { fix_boundary(x_string, x_len, threshold, y_start, real_y_start, - end_site, extra_begin, extra_end, y_id, Window_Len, R_INF, dumy, + end_site, extra_begin, extra_end, y_id, Window_Len, rref, dumy, y_strand, error, &y_start, &real_y_start, &end_site, &extra_begin, &extra_end, &error); } else { @@ -3216,6 +3236,8 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea p->extra_end = extra_end; p->error_threshold = threshold; z->align_length += x_len; w_idx[k] = z->w_list.n - 1; + + if(is_srt && z->w_list.n > 1 && p->x_start < z->w_list.a[z->w_list.n-2].x_start) is_srt = 0; } else { break; @@ -3235,10 +3257,21 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea z->is_match = 0; if((((z->x_pos_e + 1 - z->x_pos_s)*MIN_UL_ALIN_RATE) <= z->align_length) && (z->align_length >= MIN_UL_ALIN_LEN)){ 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); } } } + // fprintf(stderr, "+++[M::%s::idx->%d::y_id->%u] z::align_length->%u, e_threshold->%f\n", + // __func__, 27, overlap_list->list[27].y_id, overlap_list->list[27].align_length, e_rate); + + // fprintf(stderr, "+++[M::%s::idx->%d::y_id->%u] z::align_length->%u, e_threshold->%f\n", + // __func__, 45, overlap_list->list[45].y_id, overlap_list->list[45].align_length, e_rate); + + // fprintf(stderr, "+++[M::%s::idx->%d::y_id->%u] z::align_length->%u, e_threshold->%f\n", + // __func__, 277, overlap_list->list[277].y_id, overlap_list->list[277].align_length, e_rate); + if(uref && overlap_list->mapped_overlaps_length > 0) { set_herror_win(overlap_list, dumy, v_idx, e_rate, g_read->length, block_s); } @@ -3246,12 +3279,13 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea overlap_list->mapped_overlaps_length = 0; for (j = 0; j < (long long)overlap_list->length; j++) { z = &(overlap_list->list[j]); - y_id = z->y_id; y_strand = z->y_pos_strand; y_readLen = Get_READ_LENGTH((*R_INF), y_id); + y_id = z->y_id; y_strand = z->y_pos_strand; + y_readLen = (rref?(Get_READ_LENGTH((*rref), y_id)):(uref->ug->u.a[y_id].len)); overlap_length = z->x_pos_e + 1 - z->x_pos_s; //z->is_match = 0; - + ///debug_scan_cigar(&(overlap_list->list[j])); ///only calculate cigar for high quality overlaps - if ((R_INF && (overlap_length*OVERLAP_THRESHOLD_FILTER <= z->align_length)) || + if ((rref && (overlap_length*OVERLAP_THRESHOLD_FILTER <= z->align_length)) || (uref && (overlap_length*(1-e_rate) <= z->align_length))) { a_nw = z->w_list.n; for (i = 0, is_srt = 1; i < a_nw; i++) { @@ -3275,8 +3309,8 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea ///for the window with cigar, y_start has already reduced extra_begin y_start = p->y_start; extra_begin = p->extra_begin; extra_end = p->extra_end; o_len = Window_Len - extra_end - extra_begin; - if(R_INF) { - fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, R_INF, y_id, extra_begin, extra_end); + if(rref) { + fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, rref, y_id, extra_begin, extra_end); } else { fill_subregion_ul(dumy->overlap_region, y_start, o_len, y_strand, uref, y_id, extra_begin, extra_end); } @@ -3290,9 +3324,9 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea { if (end_site == Window_Len - 1 || real_y_start == 0) { - if(R_INF) { + if(rref) { if(fix_boundary(x_string, x_len, threshold, y_start, real_y_start, end_site, - extra_begin, extra_end, y_id, Window_Len, R_INF, dumy, y_strand, error, + extra_begin, extra_end, y_id, Window_Len, rref, dumy, y_strand, error, &y_start, &real_y_start, &end_site, &extra_begin, &extra_end, &error)) { p->error = error; p->extra_begin = extra_begin; p->extra_end = extra_end; } @@ -3323,7 +3357,7 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea } if(!is_srt) radix_sort_window_list_xs_srt(z->w_list.a, z->w_list.a + z->w_list.n); - error_rate = non_trim_error_rate(z, R_INF, uref, v_idx, dumy, g_read, e_rate, block_s); + error_rate = non_trim_error_rate(z, rref, uref, v_idx, dumy, g_read, e_rate, block_s); z->is_match = 0; if (error_rate <= e_rate_final/**asm_opt.max_ov_diff_final**/) { @@ -3333,8 +3367,8 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea // fprintf(stderr, "3-[M::%s] j::%lld, nw::%ld, a_nw::%ld, z->x_pos_s::%u, z->x_pos_e::%u, z->y_pos_s::%u, z->y_pos_e::%u, w_idx[0]::%lu, w_idx[0]->cidx::%u, w_idx[0]->clen::%u, w_idx[0]->cigar[0]:%u\n", __func__, // j, nw, a_nw, z->x_pos_s, z->x_pos_e, z->y_pos_s, z->y_pos_e, w_idx[0], z->w_list.a[w_idx[0]].cidx, z->w_list.a[w_idx[0]].clen, z->w_list.c.a[z->w_list.a[w_idx[0]].cidx]); // } - if(R_INF) { - calculate_boundary_cigars(z, R_INF, dumy, g_read, e_rate); + 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); } @@ -3354,8 +3388,8 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea } else if (error_rate <= /**asm_opt.max_ov_diff_final**/e_rate_final * 1.5) { z->is_match = 3; } - // fprintf(stderr, "[M::%s::idx->%ld::is_match->%u] z::x_pos_s->%u, z::x_pos_e->%u, error_rate->%f, e_threshold->%f\n", - // __func__, j, z->is_match, z->x_pos_s, z->x_pos_e, error_rate, e_rate); + // fprintf(stderr, "[M::%s::idx->%lld::is_match->%u] z::y_id->%u, z::x_pos_s->%u, z::x_pos_e->%u, error_rate->%f, e_threshold->%f\n", + // __func__, j, z->is_match, z->y_id, z->x_pos_s, z->x_pos_e, error_rate, e_rate); } else {///it impossible to be matched z->is_match = 0; // fprintf(stderr, "[M::%s::idx->%ld::is_match->%u] z::x_pos_s->%u, z::x_pos_e->%u, error_rate->-1, e_threshold->%f\n", @@ -3364,10 +3398,11 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea } - ///debug_window_cigar(overlap_list, g_read, dumy, R_INF, 1, 1); + ///debug_window_cigar(overlap_list, g_read, dumy, rref, 1, 1); } + inline void add_base_to_correct_read_directly(Correct_dumy* dumy, char base) { @@ -8372,9 +8407,12 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref // 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); + // fprintf(stderr, "[M::%s-beg] occ[0]->%lu, occ[1]->%lu, occ[2]->%lu, occ[3]->%lu\n", __func__, + // ovlp_occ(overlap_list, 0), ovlp_occ(overlap_list, 1), ovlp_occ(overlap_list, 2), ovlp_occ(overlap_list, 3)); ///after this function, overlap_list is sorted by x_pos_e; used for g_chain partition_ul_overlaps_advance(overlap_list, uref, g_read, overlap_read, dumy, hap, force_repeat, max_ov_diff_ec, w_inf.window_length, km); - + // fprintf(stderr, "[M::%s-end] occ[0]->%lu, occ[1]->%lu, occ[2]->%lu, occ[3]->%lu\n", __func__, + // ovlp_occ(overlap_list, 0), ovlp_occ(overlap_list, 1), ovlp_occ(overlap_list, 2), ovlp_occ(overlap_list, 3)); // debug_phasing_status(overlap_list, uref->ug, 0, hap, g_read, 20, 1176); // debug_phasing_status(overlap_list, uref->ug, 0, hap, g_read, 20, 1167); // debug_phasing_status(overlap_list, uref->ug, 0, hap, g_read, 20, 1170); diff --git a/Overlaps.cpp b/Overlaps.cpp index 8e0ae45..53d566a 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -24825,12 +24825,16 @@ char* output_file_name, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, all_ul write_debug_ma_hit_ts(sources, R_INF.total_reads, gfa_name); sprintf(gfa_name, "%s.all.debug.reverse", output_file_name); write_debug_ma_hit_ts(reverse_sources, R_INF.total_reads, gfa_name); - sprintf(gfa_name, "%s.all.debug.coverage_cut", output_file_name); - write_coverage_cut(coverage_cut, gfa_name, R_INF.total_reads); + if(coverage_cut) { + sprintf(gfa_name, "%s.all.debug.coverage_cut", output_file_name); + write_coverage_cut(coverage_cut, gfa_name, R_INF.total_reads); + } sprintf(gfa_name, "%s.all.debug.ruIndex", output_file_name); write_ruIndex(ruIndex, gfa_name); - sprintf(gfa_name, "%s.all.debug.asg_t", output_file_name); - write_asg_t(sg, gfa_name); + if(sg) { + sprintf(gfa_name, "%s.all.debug.asg_t", output_file_name); + write_asg_t(sg, gfa_name); + } if(ul_r_inf) { sprintf(gfa_name, "%s.all.debug.ul.rinfor", output_file_name); write_all_ul_t(ul_r_inf, gfa_name, NULL); @@ -24850,19 +24854,22 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_u fp = fopen(gfa_name, "r"); if(!fp) return 0; fclose(fp); sprintf(gfa_name, "%s.all.debug.reverse.bin", output_file_name); fp = fopen(gfa_name, "r"); if(!fp) return 0; fclose(fp); - sprintf(gfa_name, "%s.all.debug.coverage_cut.bin", output_file_name); - fp = fopen(gfa_name, "r"); if(!fp) return 0; fclose(fp); + if(coverage_cut) { + sprintf(gfa_name, "%s.all.debug.coverage_cut.bin", output_file_name); + fp = fopen(gfa_name, "r"); if(!fp) return 0; fclose(fp); + } sprintf(gfa_name, "%s.all.debug.ruIndex.bin", output_file_name); fp = fopen(gfa_name, "r"); if(!fp) return 0; fclose(fp); - sprintf(gfa_name, "%s.all.debug.asg_t.bin", output_file_name); - fp = fopen(gfa_name, "r"); if(!fp) return 0; fclose(fp); + if(sg) { + sprintf(gfa_name, "%s.all.debug.asg_t.bin", output_file_name); + fp = fopen(gfa_name, "r"); if(!fp) return 0; fclose(fp); + } if(ul_r_inf) { sprintf(gfa_name, "%s.all.debug.ul.rinfor.ul.ovlp.bin", output_file_name); fp = fopen(gfa_name, "r"); if(!fp) return 0; fclose(fp); } - if((sg == NULL) || (sources == NULL) || (coverage_cut == NULL) || (reverse_sources == NULL) || - (ruIndex == NULL)) + if((sources == NULL) || (reverse_sources == NULL) || (ruIndex == NULL)) { return 1; } @@ -24877,7 +24884,7 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_u destory_ma_hit_t_alloc((*reverse_sources)); } - if((*coverage_cut)!=NULL) + if(coverage_cut && (*coverage_cut)!=NULL) { free((*coverage_cut)); } @@ -24887,7 +24894,7 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_u destory_R_to_U((ruIndex)); } - if((*sg)!=NULL) + if(sg && (*sg)!=NULL) { asg_destroy(*sg); } @@ -24906,10 +24913,12 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_u return 0; } - sprintf(gfa_name, "%s.all.debug.coverage_cut", output_file_name); - if(!load_coverage_cut(coverage_cut, gfa_name)) - { - return 0; + if(coverage_cut) { + sprintf(gfa_name, "%s.all.debug.coverage_cut", output_file_name); + if(!load_coverage_cut(coverage_cut, gfa_name)) + { + return 0; + } } sprintf(gfa_name, "%s.all.debug.ruIndex", output_file_name); @@ -24918,11 +24927,14 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_u return 0; } - sprintf(gfa_name, "%s.all.debug.asg_t", output_file_name); - if(!load_asg_t(sg, gfa_name)) - { - return 0; + if(sg) { + sprintf(gfa_name, "%s.all.debug.asg_t", output_file_name); + if(!load_asg_t(sg, gfa_name)) + { + return 0; + } } + if(ul_r_inf) { sprintf(gfa_name, "%s.all.debug.ul.rinfor", output_file_name); @@ -31446,7 +31458,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) if(debug_g) { - init_bub_label_t(&b_mask_t, MIN(10, asm_opt.thread_num), sg->n_seq); + init_bub_label_t(&b_mask_t, MIN(10, asm_opt.thread_num), n_read); goto debug_gfa; } ///just for debug @@ -31468,6 +31480,11 @@ ma_sub_t **coverage_cut_ptr, int debug_g) { memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t)); } + + // if (asm_opt.flag & HA_F_VERBOSE_GFA) { + // write_debug_graph(NULL, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF); + // debug_gfa:; + // } ///should recover edges from sources by using UL alignments if(asm_opt.ar) create_ul_info(sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex); @@ -31510,7 +31527,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g) write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF); debug_gfa:; gen_ug_opt_t(&uopt, sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex); - // if(asm_opt.ar) create_ul_info(sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex); } if(asm_opt.ar) ul_realignment_gfa(&uopt, sg); print_debug_gfa(sg, NULL, coverage_cut, "UL.debug", sources, ruIndex, max_hang_length, mini_overlap_length); @@ -31722,7 +31738,7 @@ long long bubble_dist, int read_graph, int write) min_thres = asm_opt.max_short_tip + 1; if (asm_opt.flag & HA_F_VERBOSE_GFA) { - if(load_debug_graph(&sg, &sources, &coverage_cut, output_file_name, &reverse_sources, &ruIndex, &UL_INF)) + if(load_debug_graph(/**NULL**/&sg, &sources, /**NULL**/&coverage_cut, output_file_name, &reverse_sources, &ruIndex, &UL_INF)) { fprintf(stderr, "debug gfa has been loaded\n"); diff --git a/inter.cpp b/inter.cpp index 3df8ef4..cc1bcc5 100644 --- a/inter.cpp +++ b/inter.cpp @@ -2697,6 +2697,7 @@ int64_t gen_contain_chain(const ul_idx_t *uref, utg_ct_t *p, overlap_region* o, /**x->tn = p->x>>1;**/x->tn = (uint32_t)(0x80000000); x->tn |= (p->x>>1); x->ts = q_s; x->te = q_e; x->el = 1;x->sec = 0; x->rev = ((o->y_pos_strand == (p->x&1))?0:1); + // if(((x->tn<<1)>>1) == 23113) fprintf(stderr, "x->tn:%u, o->y_id:%u\n", (x->tn<<1)>>1, o->y_id); // if(x->qn == 0 /**&& ((x->tn<<1)>>1) == 302**/) { // /**if(o->x_id == 0 && (o->y_id == 46 || o->y_id == 48))**/ { // // fprintf(stderr, "\nUL[%u]\t%u\t%u\t%c\tUTG[%u]\t%u\t%u\n", o->x_id, o->x_pos_s, o->x_pos_e, @@ -4016,33 +4017,37 @@ void *km) ul_ov_t *m = &(idx->a[idx->n-1]); //largest chain // fprintf(stderr, "[M::%s] m->score:%u, m->qs:%u, m->qe:%u, chain_n:%u\n", __func__, m->qn, m->qs, m->qe, m->te-m->ts); if((m->qe-m->qs) <= (qlen*cov_rate)) return 0; - if(check_trans_rate(a+m->ts, m->te-m->ts, trans_thres)) return 1; - if(olist && hap && uref) { - int64_t idx_n = idx->n, z, i, het_n, resc_tk = 0, f = 0; - uint64_t si; ma_utg_t *u = NULL; - for (z = m->ts; z < m->te; z++) { - if(a[z].el) { - kv_push_km(km, ul_ov_t, *idx, a[z]); - } else { - i = a[z].qn; si = 0; - het_n = update_ava_het_site(hap, i, &si, NULL, 1); - assert(het_n > 0 && olist->list[i].is_match == 2); - u = &(uref->ug->u.a[olist->list[i].y_id]); - if(u->n > 1) { - resc_tk += rescue_trans_ul_chains(uref, &(olist->list[i]), hap->list+si, het_n, u, - idx, diff_ec_ul, winLen, 0, NULL, km); + if(trans_thres >= 0) { + if(check_trans_rate(a+m->ts, m->te-m->ts, trans_thres)) return 1; + if(olist && hap && uref) { + int64_t idx_n = idx->n, z, i, het_n, resc_tk = 0, f = 0; + uint64_t si; ma_utg_t *u = NULL; + for (z = m->ts; z < m->te; z++) { + if(a[z].el) { + kv_push_km(km, ul_ov_t, *idx, a[z]); + } else { + i = a[z].qn; si = 0; + het_n = update_ava_het_site(hap, i, &si, NULL, 1); + assert(het_n > 0 && olist->list[i].is_match == 2); + u = &(uref->ug->u.a[olist->list[i].y_id]); + if(u->n > 1) { + resc_tk += rescue_trans_ul_chains(uref, &(olist->list[i]), hap->list+si, het_n, u, + idx, diff_ec_ul, winLen, 0, NULL, km); + } } } - } - if(resc_tk) { - radix_sort_ul_ov_srt_qe(idx->a+idx_n, idx->a+idx->n); - f = check_trans_rate(idx->a+idx_n, idx->n-idx_n, trans_thres); + if(resc_tk) { + radix_sort_ul_ov_srt_qe(idx->a+idx_n, idx->a+idx->n); + f = check_trans_rate(idx->a+idx_n, idx->n-idx_n, trans_thres); + } + idx->n = idx_n; + return f; } - idx->n = idx_n; - return f; + return 0; + } else { + return 1; } - return 0; } void dump_chain(kv_ul_ov_t *des, ul_ov_t *src, ul_ov_t *chain, void *km) @@ -4353,7 +4358,7 @@ int64_t debug_i, void *km) ///the first two rounds could reuse dumy->overlapID. But for the last round, dumy->overlapID is not long enough occ = gl_chain_advance(idx, ll->tk.a+ll->tk.n, uref, uopt, G_CHAIN_BW, diff_ec_ul, qlen, UG_SKIP, dumy->overlapID, ll->srt.a.a, hap->snp_srt.a, G_CHAIN_TRANS_WEIGHT, 0, NULL, uref->ug, debug_i, km); if(occ) { - if(ff_chain(idx, qlen, P_CHAIN_COV, G_CHAIN_TRANS_RATE, ll->tk.a+ll->tk.n, NULL, NULL, NULL, diff_ec_ul, winLen, km)) { + if(ff_chain(idx, qlen, P_CHAIN_COV, -1/**G_CHAIN_TRANS_RATE**/, ll->tk.a+ll->tk.n, NULL, NULL, NULL, diff_ec_ul, winLen, km)) { f = 1; //dump_chain(idx, ll->tk.a+ll->tk.n, &(idx->a[idx->n-1]), km); for (k = idx->a[idx->n-1].ts; k < idx->a[idx->n-1].te; k++) { olist->list[ll->tk.a[ll->tk.n+k].qn].x_pos_strand = 1; @@ -5011,7 +5016,8 @@ void debug_ul_vec_t_chain(void *km, const asg_t *g, ul_vec_t *rch, st_mt_t *dst_ } } -int64_t gl_chain_refine_advance_combine(mg_tbuf_t *b, ul_vec_t *rch, overlap_region_alloc* olist, Correct_dumy* dumy, haplotype_evdience_alloc *hap, st_mt_t *sps, glchain_t *ll, gdpchain_t *gdp, const ul_idx_t *uref, double diff_ec_ul, int64_t winLen, int64_t qlen, const ug_opt_t *uopt, + +int64_t gl_chain_refine_advance_combine_with_trans(mg_tbuf_t *b, ul_vec_t *rch, overlap_region_alloc* olist, Correct_dumy* dumy, haplotype_evdience_alloc *hap, st_mt_t *sps, glchain_t *ll, gdpchain_t *gdp, const ul_idx_t *uref, double diff_ec_ul, int64_t winLen, int64_t qlen, const ug_opt_t *uopt, int64_t debug_i, int64_t tid, void *km) { ll->tk.n = ll->lo.n = 0; @@ -5027,7 +5033,7 @@ int64_t debug_i, int64_t tid, void *km) ///chain exact U-matches occ = gl_chain_advance(idx, ll->tk.a, uref, uopt, G_CHAIN_BW, /**diff_ec_ul**/N_GCHAIN_RATE, qlen, UG_SKIP, dumy->overlapID, ll->srt.a.a, hap->snp_srt.a, G_CHAIN_TRANS_WEIGHT, 0, NULL, uref->ug, debug_i, km); if(occ) { - if(ff_chain(idx, qlen, 0.99/**P_CHAIN_COV**/, G_CHAIN_TRANS_RATE, ll->tk.a, NULL, NULL, NULL, diff_ec_ul, winLen, km)) { + if(ff_chain(idx, qlen, 0.99/**P_CHAIN_COV**/, -1/**G_CHAIN_TRANS_RATE**/, ll->tk.a, NULL, NULL, NULL, diff_ec_ul, winLen, km)) { f = l2g_res_chain(uref->ug, ll->tk.a+idx->a[idx->n-1].ts, idx->a[idx->n-1].te-idx->a[idx->n-1].ts, &(gdp->swap), -1/**N_GCHAIN_RATE**/); } else if(o2) {///means there are trans overlaps gl_chain_gen(olist, uref, idx, 1, hap, km); @@ -5063,6 +5069,49 @@ int64_t debug_i, int64_t tid, void *km) return 1; } +int64_t gl_chain_refine_advance_combine(mg_tbuf_t *b, ul_vec_t *rch, overlap_region_alloc* olist, Correct_dumy* dumy, haplotype_evdience_alloc *hap, st_mt_t *sps, glchain_t *ll, gdpchain_t *gdp, const ul_idx_t *uref, double diff_ec_ul, int64_t winLen, int64_t qlen, const ug_opt_t *uopt, +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 + 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; + + kv_resize_km(km, uint64_t, ll->srt.a, idx->n); + kv_resize_km(km, uint64_t, hap->snp_srt, idx->n); + kv_resize_km(km, ul_ov_t, ll->tk, idx->n); + ///chain exact U-matches + occ = gl_chain_advance(idx, ll->tk.a, uref, uopt, G_CHAIN_BW, /**diff_ec_ul**/N_GCHAIN_RATE, qlen, UG_SKIP, dumy->overlapID, ll->srt.a.a, hap->snp_srt.a, G_CHAIN_TRANS_WEIGHT, 0, NULL, uref->ug, debug_i, km); + if(occ) { + if(ff_chain(idx, qlen, 0.99/**P_CHAIN_COV**/, -1/**G_CHAIN_TRANS_RATE**/, ll->tk.a, NULL, NULL, NULL, diff_ec_ul, winLen, km)) { + f = l2g_res_chain(uref->ug, ll->tk.a+idx->a[idx->n-1].ts, idx->a[idx->n-1].te-idx->a[idx->n-1].ts, &(gdp->swap), -1/**N_GCHAIN_RATE**/); + } + } + // 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 + l2g_chain(uref, idx, &(gdp->l)); ll->tk.n = 0; + ///buffer + kv_resize(uint64_t, ll->srt.a, gdp->l.n); kv_resize(uint64_t, hap->snp_srt, gdp->l.n); + kv_resize(uint64_t, gdp->v, gdp->l.n); kv_resize(int64_t, gdp->f, gdp->l.n); + max_idx = hc_gchain1_dp(b->km, uref, uref->ug, &(gdp->l), &(gdp->swap), &(gdp->dst), &(gdp->out), &(gdp->path), rch->rlen, + uopt, G_CHAIN_BW, diff_ec_ul, -1, ll->srt.a.a, sps, gdp->f.a, hap->snp_srt.a, gdp->v.a); + if(max_idx >= 0 && gen_max_gchain_adv(b->km, uref, debug_i, sps, &(gdp->l), &(ll->tk), NULL, rch->rlen, P_CHAIN_COV, 0.3/**P_FRAGEMENT_PRIMARY_CHAIN_COV**/, + 0.1/**P_FRAGEMENT_PRIMARY_SECOND_COV**/, PRIMARY_UL_CHAIN_MIN, uref->ug->g, &(gdp->dst_done), &(gdp->out), &(gdp->path), ll->srt.a.a, &(gdp->swap))) { + // print_raw_chains(&(gdp->swap), debug_i); + // f = check_trans_rate_gap(&(gdp->swap), &(ll->tk), olist, hap, uref, diff_ec_ul, winLen, G_CHAIN_TRANS_RATE); + f = 1; + } + } + // if(debug_i == 1756) fprintf(stderr, "[M::%s] ulid:%ld, qlen:%ld, f:%ld\n", __func__, debug_i, qlen, f); + if(f) update_ul_vec_t_ug(uref, rch, &(gdp->swap), debug_i); + // debug_ul_vec_t_chain(km, uref->ug->g, rch, &(gdp->dst_done), &(gdp->out)); + // fprintf(stderr, "(beg3) [M::%s::tid:%ld] debug_i:%ld, qlen:%ld\n", __func__, tid, debug_i, qlen); + 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; @@ -5121,7 +5170,7 @@ static void worker_for_ul_scall_alignment(void *data, long i, int tid) // callba uint64_t align = 0; int fully_cov, abnormal; void *km = s->buf?(s->buf[tid]?s->buf[tid]->km:NULL):NULL; - // if(s->id+i!=601) return; + // if(s->id+i!=3196) return; // fprintf(stderr, "[M::%s] rid:%ld\n", __func__, s->id+i); // if (memcmp(UL_INF.nid.a[s->id+i].a, "d0aab024-b3a7-40fb-83cc-22c3d6d951f8", UL_INF.nid.a[s->id+i].n-1)) return; // fprintf(stderr, "[M::%s::] ==> len: %lu\n", __func__, s->len[i]); @@ -5156,7 +5205,7 @@ static void worker_for_ul_scall_alignment(void *data, long i, int tid) // callba free(s->seq[i]); s->seq[i] = NULL; } b->num_correct_base += align; - + // exit(1); // uint64_t k; // b->num_read_base += overlap_statistics(&b->olist, NULL, NULL, 1); // for (k = 0; k < bl->tk.n; k++) { @@ -5182,9 +5231,12 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call int64_t /**rid = s->id+i,**/ winLen = MIN((((double)THRESHOLD_MAX_SIZE)/s->opt->diff_ec_ul), WINDOW); // 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); + } assert(UL_INF.a[s->id+i].rlen == s->len[i]); // void *km = s->buf?(s->buf[tid]?s->buf[tid]->km:NULL):NULL; - // if(s->id+i!=873) return; + // if(s->id+i!=3373) return; // fprintf(stderr, "\n[M::%s] rid:%ld, len:%lu\n", __func__, s->id+i, s->len[i]); // if (memcmp(UL_INF.nid.a[s->id+i].a, "d0aab024-b3a7-40fb-83cc-22c3d6d951f8", UL_INF.nid.a[s->id+i].n-1)) return; // fprintf(stderr, "[M::%s::] ==> len: %lu\n", __func__, s->len[i]); @@ -5218,16 +5270,16 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call free(s->seq[i]); s->seq[i] = NULL; b->num_correct_base++; } s->hab[tid]->num_read_base++; - - int64_t mem[6], mem_hab[6]; - if(get_utepdat_t_mem_tid(s, tid, mem, mem_hab)>((int64_t)5*(int64_t)1073741824)) { - fprintf(stderr, "[M::%s::tid->%d::rid->%ld] buffer[0]: %.3fGB(%.3fGB::%.3fGB::%.3fGB::%.3fGB::%.3fGB), buffer[1]: %.3fGB, buffer[2]: %.3fGB, buffer[3]: %.3fGB, buffer[4]: %.3fGB, buffer[5]: %.3fGB\n", - __func__, tid, i, mem[0]/1073741824.0, - mem_hab[0]/1073741824.0, mem_hab[1]/1073741824.0, mem_hab[2]/1073741824.0, - mem_hab[3]/1073741824.0, mem_hab[4]/1073741824.0, - mem[1]/1073741824.0, mem[2]/1073741824.0, - mem[3]/1073741824.0, mem[4]/1073741824.0, mem[5]/1073741824.0); - } + // fprintf(stderr, "[M::%s] rid:%ld, dd:%u\n", __func__, s->id+i, UL_INF.a[s->id+i].dd); + // int64_t mem[6], mem_hab[6]; + // if(get_utepdat_t_mem_tid(s, tid, mem, mem_hab)>((int64_t)5*(int64_t)1073741824)) { + // fprintf(stderr, "[M::%s::tid->%d::rid->%ld] buffer[0]: %.3fGB(%.3fGB::%.3fGB::%.3fGB::%.3fGB::%.3fGB), buffer[1]: %.3fGB, buffer[2]: %.3fGB, buffer[3]: %.3fGB, buffer[4]: %.3fGB, buffer[5]: %.3fGB\n", + // __func__, tid, i, mem[0]/1073741824.0, + // mem_hab[0]/1073741824.0, mem_hab[1]/1073741824.0, mem_hab[2]/1073741824.0, + // mem_hab[3]/1073741824.0, mem_hab[4]/1073741824.0, + // mem[1]/1073741824.0, mem[2]/1073741824.0, + // mem[3]/1073741824.0, mem[4]/1073741824.0, mem[5]/1073741824.0); + // } // align = kv_ul_ov_t_statistics(&(bl->tk), i, &(b->num_recorrect_base)); // if(align == s->len[i]) { @@ -5525,11 +5577,11 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac free(s->ll[i].tk.a); } for (i = 0; i < (uint64_t)s->n; ++i) { - rid = s->id + i; - if(UL_INF.n > rid && UL_INF.a[rid].rlen != s->len[i]) { + rid = s->id + i; + if((UL_INF.n <= rid) || (UL_INF.n > rid && UL_INF.a[rid].rlen != s->len[i])) { append_ul_t(&UL_INF, &rid, NULL, 0, s->seq[i], s->len[i], NULL, 0, P_CHAIN_COV, s->uopt, 0); - } - free(s->seq[i]); + } + free(s->seq[i]); } fprintf(stderr, "[M::%s::dump_done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); /** @@ -5585,6 +5637,9 @@ static void *worker_ul_rescall_pipeline(void *data, int step, void *in) // callb s->sum_len += l; memcpy(s->seq[s->n], p->ks->seq.s, l); // fprintf(stderr, "s->n->%d, l->%lu\n", s->n, l); + // if(s->id + s->n == 13706) { + // fprintf(stderr, "+++++rid->%lu, l->%lu, %.*s\n", + // s->id + s->n, l, (int32_t)p->ks->name.l, p->ks->name.s);} s->len[s->n++] = l; if (s->sum_len >= p->chunk_size) break; } @@ -5630,6 +5685,7 @@ static void *worker_ul_rescall_pipeline(void *data, int step, void *in) // callb // if(s->seq[i] == NULL) fprintf(stderr, "[M::%s::]rid->%ld, len->%lu\n", __func__, rid, s->len[i]); write_compress_base_disk(p->ucr_s->fp, rid, s->seq[i], s->len[i], &(p->ucr_s->u)); } + // if(UL_INF.a[rid].dd) fprintf(stderr, "rid->%ld\n", rid); free(s->seq[i]); } fprintf(stderr, "[M::%s::dump_done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n); @@ -8241,6 +8297,17 @@ static void worker_for_ul_gchains_alignment(void *data, long i, int tid) // gl_chain_refine_advance(&b->olist, &b->correct, &b->hap, bl, s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, s->id+i, km); } +void detect_outlier_len(const char* cmd) +{ + uint64_t k; + for (k = 0; k < UL_INF.n; k++) { + if(UL_INF.a[k].rlen == 0) { + fprintf(stderr, "[%s] rid->%lu, rlen->%u, %.*s\n", + cmd, k, UL_INF.a[k].rlen, (int32_t)UL_INF.nid.a[k].n, UL_INF.nid.a[k].a); + } + } +} + uint64_t work_ul_gchains(uldat_t *sl) { utepdat_t s; uint64_t i; memset(&s, 0, sizeof(s)); @@ -8252,8 +8319,12 @@ uint64_t work_ul_gchains(uldat_t *sl) s.hab[i] = ha_ovec_init(0, 0, 1); s.buf[i] = mg_tbuf_init(); } + // detect_outlier_len("+++work_ul_gchains"); + kt_for(sl->n_thread, worker_for_ul_gchains_alignment, &s, UL_INF.n); + // detect_outlier_len("---work_ul_gchains"); + for (i = 0; i < sl->n_thread; ++i) { s.sum_len += s.hab[i]->num_read_base; s.n += s.hab[i]->num_correct_base; ha_ovec_destroy(s.hab[i]); mg_tbuf_destroy(s.buf[i]); hc_glchain_destroy(&(s.ll[i])); @@ -8283,9 +8354,9 @@ void print_ul_ovlps(all_ul_t *x, int32_t prt_ovlp) fprintf(stderr, "B\t%.*s\t%u\t%u\t%u\n", (int32_t)z->n, z->a, p->rlen, (m->qs+((m->hid>>15)&FLANK_M)), (m->qe-(m->hid&FLANK_M))); } else { - fprintf(stderr, "A\t%.*s\t%u\t%u\t%u\t%c\t%.*s\t%u\t%u\t%u\n", - (int32_t)z->n, z->a, p->rlen, m->qs, m->qe, "+-"[m->rev], - (int32_t)Get_NAME_LENGTH(R_INF, m->hid), Get_NAME(R_INF, m->hid), + fprintf(stderr, "A\t%.*s(%lu)\t%u\t%u\t%u\t%c\t%.*s(%u)\t%u\t%u\t%u\n", + (int32_t)z->n, z->a, k, p->rlen, m->qs, m->qe, "+-"[m->rev], + (int32_t)Get_NAME_LENGTH(R_INF, m->hid), Get_NAME(R_INF, m->hid), m->hid, (uint32_t)Get_READ_LENGTH(R_INF, m->hid), m->ts, m->te); if(m->el) cov_occ++; } @@ -8345,7 +8416,12 @@ void print_ovlp_src_bl_stat(all_ul_t *x, const ug_opt_t *uopt) __func__, R_INF.total_reads, tc, ta); uint64_t tt[4] = {0}; - for (k = 0; k < x->n; k++) tt[x->a[k].dd]++; + for (k = 0; k < x->n; k++) { + tt[x->a[k].dd]++; + // if(x->a[k].dd == 1) { + // fprintf(stderr, "(%lu) %.*s\n", k, (int32_t)x->nid.a[k].n, x->nid.a[k].a); + // } + } fprintf(stderr, "[M::%s::] ==> # passed UL reads:%lu, # fully corrected UL reads:%lu, # almost fully corrected UL reads:%lu, # UL reads have primary chains:%lu\n", __func__, tt[0]+tt[1]+tt[2]+tt[3], tt[1], tt[2], tt[3]); @@ -9931,6 +10007,8 @@ int32_t load_all_ul_t(all_ul_t *x, char* file_name, All_reads *hR, ma_ug_t *ug) return 1; } + + void ul_load(const ug_opt_t *uopt) { fprintf(stderr, "[M::%s::] ==> UL\n", __func__); @@ -9945,7 +10023,7 @@ void ul_load(const ug_opt_t *uopt) gen_UL_ovlps(&sl, cutoff); write_all_ul_t(&UL_INF, asm_opt.output_file_name, NULL); } - + // detect_outlier_len("ul_load"); // print_all_ul_t_stat(&UL_INF); // fprintf(stderr, "**1**\n"); kt_for(sl.n_thread, update_ovlp_src, &sl, R_INF.total_reads); @@ -9954,7 +10032,9 @@ void ul_load(const ug_opt_t *uopt) // fprintf(stderr, "**3**\n"); print_ovlp_src_bl_stat(&UL_INF, sl.uopt); - // print_ul_ovlps(&UL_INF, 0); print_ul_ovlps(&UL_INF, 1); + // exit(1); + // print_ul_ovlps(&UL_INF, 0); + // print_ul_ovlps(&UL_INF, 1); // destory_all_ul_t(&UL_INF); } @@ -10012,7 +10092,7 @@ ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg) ma_ug_t *ug = gen_polished_ug(uopt, sg); // dd_ug(sg, ug, uopt->coverage_cut, uopt->sources, uopt->ruIndex, "UL.sa"); // debug_sl_compress_base_disk_0(&sl, asm_opt.ar); - + // detect_outlier_len("ul_realignment"); if(!load_all_ul_t(&UL_INF, gfa_name, &R_INF, ug)) { gen_UL_reovlps(&sl, ug, sg, gfa_name, cutoff); write_all_ul_t(&UL_INF, gfa_name, ug);