xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-10-11 12:10:55 +08:00
Code Issues Packages Projects Releases Wiki Activity

small bugs fixed

Browse Source
This commit is contained in:
chhylp123
2022-05-14 22:07:55 -04:00
parent 832e43fe3b
commit 74ad15846e
4 changed files with 281 additions and 147 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.h
+1 -1
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.16.4-r4000"
#define HA_VERSION "0.16.4-r402"
#define VERBOSE 0
Correct.cpp
+108 -70
View File
@@ -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);
@@ -1097,12 +1101,17 @@ inline double non_trim_error_rate(overlap_region *z, All_reads* R_INF, const ul_
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;
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;
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);
Overlaps.cpp
+23 -7
View File
@@ -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);
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);
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);
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);
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,11 +24913,13 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_u
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);
if(!load_ruIndex(ruIndex, gfa_name))
@@ -24918,11 +24927,14 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_u
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");
inter.cpp
+104 -24
View File
@@ -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,6 +4017,7 @@ 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(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;
@@ -4043,6 +4045,9 @@ void *km)
return f;
}
return 0;
} else {
return 1;
}
}
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]) {
@@ -5526,7 +5578,7 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac
}
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]) {
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]);
@@ -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);