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

debug cigar_offset

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