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

cuttoff before debugging

Browse Source
This commit is contained in:
chhylp123
2022-11-05 09:03:23 -04:00
parent e66c003b65
commit 612c1ab712
8 changed files with 1115 additions and 93 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.17.0-r427"
#define HA_VERSION "0.17.1-r428"
#define VERBOSE 0
Correct.cpp
+80
View File
@@ -15213,6 +15213,75 @@ UC_Read *tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql,
#define ovlp_cur_xoff(x) ((x).qs)
#define ovlp_cur_coff(x) ((x).qe)
int64_t retrieve_cigar_err_debug(bit_extz_t *ez, int64_t s, int64_t e)
{
if(!ez->cigar.n) return 0;
int64_t err = 0, xk = ez->ts; int64_t ws, we, os, oe, ovlp;
uint32_t ck = 0, cl; uint16_t op;
os = MAX(s, ez->ts); oe = MIN(e, ez->te+1);
ovlp = ((oe>os)? (oe-os):0);
if(!ovlp) {
fprintf(stderr, "[M::%s::] s::%ld, e::%ld, ez->ts::%u, ez->te::%u\n",
__func__, s, e, ez->ts, ez->te);
}
assert(ovlp);
//some cigar will span s or e
while (ck < ez->cigar.n && xk < e) {//[s, e)
ws = xk;
ck = pop_trace(&(ez->cigar), ck, &op, &cl);
if(op!=2) xk += cl;
we = xk;
os = MAX(s, ws); oe = MIN(e, we);
ovlp = ((oe>os)? (oe-os):0);
if((op==2) && (ws>=s) && (ws<e)) {
ovlp = cl;
}
if((!ovlp) || (!op)) continue;
err += ovlp;
}
return err;
}
///[s, e)
int64_t extract_sub_cigar_err_debug(overlap_region *z, int64_t s, int64_t e)
{
int64_t wk = 0, wn = z->w_list.n, ws, we, os, oe, ovlp;
window_list *m; int64_t xl, yl, werr, err; bit_extz_t ez;
for (wk = err = 0; wk < wn; wk++) {
m = &(z->w_list.a[wk]);
ws = m->x_start; we = m->x_end+1;
if(ws >= e) break;
os = MAX(s, ws); oe = MIN(e, we);
ovlp = ((oe>os)? (oe-os):0);
if(ovlp) {
xl = m->x_end+1-m->x_start;
yl = m->y_end+1-m->y_start;
if((is_ualn_win((*m))) || (is_est_aln((*m)))) {
if(is_ualn_win((*m))) { //unmapped
werr = gen_err_unaligned(xl, yl);
} else {
werr = m->error;//shared window
}
if(ovlp < xl) {
werr = (((double)ovlp)/((double)xl))*((double)werr);
}
//skip the whole window
err += werr;
} else {
if(ovlp == xl) {
//skip the whole window
err += m->error;
} else {
set_bit_extz_t(ez, (*z), wk);
err += retrieve_cigar_err_debug(&ez, os, oe);
}
}
}
}
return err;
}
int64_t retrieve_cigar_err(bit_extz_t *ez, int64_t s, int64_t e, int64_t *xk, int64_t *ck)
{
@@ -15243,6 +15312,13 @@ int64_t retrieve_cigar_err(bit_extz_t *ez, int64_t s, int64_t e, int64_t *xk, in
if((!ovlp) || (!op)) continue;
err += ovlp;
}
// int64_t debug_err = retrieve_cigar_err_debug(ez, s, e);
// if(!(err == debug_err)) {
// fprintf(stderr, "[M::%s::] err::%ld, debug_err::%ld, s::%ld, e::%ld, ez->ts::%u, ez->te::%u\n",
// __func__, err, debug_err, s, e, ez->ts, ez->te);
// }
// assert(err == debug_err);
return err;
}
///[s, e)
@@ -15401,6 +15477,10 @@ uint64_t gen_region_phase(overlap_region* ol, uint64_t *id_a, uint64_t id_n, uin
// fprintf(stderr, "+++[M::%s::utg%.6dl] wid::%u, xoff::%u, coff::%u\n", __func__,
// (int32_t)ol[ovlp_id(*p)].y_id+1, ovlp_cur_wid(*p), ovlp_cur_xoff(*p), ovlp_cur_coff(*p));
err = extract_sub_cigar_err(z, s, e, p);
// int64_t debug_err = extract_sub_cigar_err_debug(z, s, e);
// assert(err == debug_err);
// fprintf(stderr, "[M::%s::] err::%ld, debug_err::%ld\n", __func__, err, debug_err);
// fprintf(stderr, "---[M::%s::utg%.6dl] wid::%u, xoff::%u, coff::%u, err::%ld\n", __func__,
// (int32_t)ol[ovlp_id(*p)].y_id+1, ovlp_cur_wid(*p), ovlp_cur_xoff(*p), ovlp_cur_coff(*p), err);
assert(err >= 0);
Correct.h
+1 -1
View File
@@ -17,7 +17,7 @@
#define MISMATCH 1
#define INSERTION 2
#define DELETION 3
#define ERROR_RATE 1.25
#define ERROR_RATE 1.5
#define UL_TOPN 50
#define SGAP 16
#define MAX_LGAP(ql) ((((ql)*0.2)<256)?((ql)*0.2):256)
Overlaps.cpp
-1
View File
@@ -10238,7 +10238,6 @@ static void update_weak_by_contain(void *data, long i, int tid)
}
}
void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long num_sources, uint32_t ou_thres)
{
double startTime = Get_T();
Process_Read.cpp
+12 -1
View File
@@ -1111,7 +1111,18 @@ void determine_chain_distance(ul_ov_t *o, int64_t on, ul_vec_t *p, ma_hit_t_allo
}
}
// if(t>=x->length) {
// fprintf(stderr, "m->%ld, i->%ld, rid->%ld, on->%ld\n", m, i, rid, on);
// fprintf(stderr, "m->%ld(id::%u), i->%ld(id::%u), rid->%ld, on->%ld\n",
// m, (o[m].tn<<1)>>1, i, (o[i].tn<<1)>>1, rid, on);
// // fprintf(stderr, "[i::%ld]id::%u\t\tq::[%u, %u)\tt::[%u, %u)\n",
// // i, o[i].tn, o[i].qs, o[i].qe, o[i].ts, o[i].te);
// // fprintf(stderr, "[m::%ld]id::%u\t\tq::[%u, %u)\tt::[%u, %u)\n",
// // m, o[m].tn, o[m].qs, o[m].qe, o[m].ts, o[m].te);
// fprintf(stderr, "[i::%ld]id::%u\t%.*s\t%c\tq::[%u, %u)\tt::[%u, %u)\n",
// i, (o[i].tn<<1)>>1, (int)Get_NAME_LENGTH(R_INF, ((o[i].tn<<1)>>1)), Get_NAME(R_INF, ((o[i].tn<<1)>>1)),
// "+-"[o[i].rev], o[i].qs, o[i].qe, o[i].ts, o[i].te);
// fprintf(stderr, "[m::%ld]id::%u\t%.*s\t%c\tq::[%u, %u)\tt::[%u, %u)\n",
// m, (o[m].tn<<1)>>1, (int)Get_NAME_LENGTH(R_INF, ((o[m].tn<<1)>>1)), Get_NAME(R_INF, ((o[m].tn<<1)>>1)),
// "+-"[o[m].rev], o[m].qs, o[m].qe, o[m].ts, o[m].te);
// exit(1);
// }
assert(t<x->length);
gfa_ut.cpp
+57 -2
View File
@@ -2352,6 +2352,48 @@ void des_flex_asg_t(flex_asg_t *z)
free(z->idx); free(z->pi.a); free(z->a);
}
void print_raw_u2rgfa_seq(all_ul_t *aln, R_to_U* rI, uint32_t is_detail)
{
uint64_t id, a_n, k, z; uc_block_t *a = NULL;
kvec_t(uint8_t) ff; kv_init(ff);
for (id = 0; id < aln->n; id++) {
a = aln->a[id].bb.a; a_n = aln->a[id].bb.n;
if(a_n == 0) continue;
fprintf(stderr,"\n%.*s\tid::%lu\trlen::%u", (int32_t)aln->nid.a[id].n, aln->nid.a[id].a, id, aln->a[id].rlen);
kv_resize(uint8_t, ff, a_n); memset(ff.a, 0, a_n*sizeof((*(ff.a))));
if(is_detail) {
fprintf(stderr, "\n");
for (k = 0; k < a_n; k++) {
if(ff.a[k]) continue;
for (z = k; z != (uint32_t)-1; z = a[z].aidx) {
fprintf(stderr, "%.*s\t%c\tq::[%u, %u)\tt::[%u, %u)\tid::%u\ttl::%lu\tc::%u\n",
(int)Get_NAME_LENGTH(R_INF, a[z].hid), Get_NAME(R_INF, a[z].hid), "+-"[a[z].rev],
a[z].qs, a[z].qe, a[z].ts, a[z].te, a[z].hid, Get_READ_LENGTH(R_INF, a[z].hid),
rI?is_contain_r((*rI), a[z].hid):0);
assert(ff.a[z] == 0);
ff.a[z] = 1;
}
fprintf(stderr, "************\n");
}
} else {
fprintf(stderr, "\t");
for (k = 0; k < a_n; k++) {
if(ff.a[k]) continue;
for (z = k; z != (uint32_t)-1; z = a[z].aidx) {
fprintf(stderr, "%.*s\t",
(int)Get_NAME_LENGTH(R_INF, a[z].hid), Get_NAME(R_INF, a[z].hid));
assert(ff.a[z] == 0);
ff.a[z] = 1;
}
fprintf(stderr, "\n");
}
}
}
kv_destroy(ff);
}
void ul_clean_gfa(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, int64_t clean_round, double min_ovlp_drop_ratio, double max_ovlp_drop_ratio,
double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, int32_t is_ou, int32_t is_trio, uint32_t ou_thres, char *o_file)
{
@@ -2367,7 +2409,20 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
// if(is_ou) update_sg_uo(sg, src);///do not do it here
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
// exit(1);
// print_raw_u2rgfa_seq(&UL_INF, rI, 1);
// exit(1);
// fprintf(stderr, "%.*s\tid::%u\tis_c::%u\n",
// (int)Get_NAME_LENGTH(R_INF, 10785), Get_NAME(R_INF, 10785), 10785, is_contain_r((*rI), 10785));
// fprintf(stderr, "%.*s\tid::%u\tis_c::%u\n",
// (int)Get_NAME_LENGTH(R_INF, 10790), Get_NAME(R_INF, 10790), 10790, is_contain_r((*rI), 10790));
// fprintf(stderr, "%.*s\tid::%u\tis_c::%u\n",
// (int)Get_NAME_LENGTH(R_INF, 10805), Get_NAME(R_INF, 10805), 10805, is_contain_r((*rI), 10805));
// fprintf(stderr, "%.*s\tid::%u\tis_c::%u\n",
// (int)Get_NAME_LENGTH(R_INF, 10809), Get_NAME(R_INF, 10809), 10809, is_contain_r((*rI), 10809));
// fprintf(stderr, "%.*s\tid::%u\tis_c::%u\n",
// (int)Get_NAME_LENGTH(R_INF, 10819), Get_NAME(R_INF, 10819), 10819, is_contain_r((*rI), 10819));
// debug_info_of_specfic_node("m64012_190921_234837/111673711/ccs", sg, rI, "beg");
// debug_info_of_specfic_node("m64011_190830_220126/95028102/ccs", sg, rI, "beg");
@@ -11794,7 +11849,7 @@ void u2g_hybrid_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, usg_t *ng, asg64_v *
debug_sysm_usg_t(ng, __func__);
/******for debug******/
prt_usg_t(uidx, ng, "ng1");
// prt_usg_t(uidx, ng, "ng1");
/******for debug******/
// u2g_hybrid_extend(ng, NULL, b, ub);
@@ -11994,7 +12049,7 @@ void u2g_threading(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint64_t cov_cutoff, as
idx->h_usg = ng;
/******for debug******/
prt_usg_t(uidx, ng, "ng0");
// prt_usg_t(uidx, ng, "ng0");
/******for debug******/
gfa_ut.h
+1
View File
@@ -28,5 +28,6 @@ double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_
void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t max_hang, int64_t min_ovlp, int64_t ul_occ);
void normalize_gou(asg_t *g);
void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd);
// void print_raw_u2rgfa_seq(all_ul_t *aln, R_to_U* rI, uint32_t is_detail);
#endif
inter.cpp
+963 -87
View File
File diff suppressed because it is too large Load Diff