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

fix missed ug when bubble popping

Browse Source
This commit is contained in:
chhylp123
2021-05-19 10:16:23 -04:00
parent b12e008416
commit 5c1680e998
2 changed files with 6 additions and 107 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Overlaps.cpp
+4 -105
View File
@@ -12712,104 +12712,6 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp, kvec_asg_a
/*******************************for debug************************************/
}
void hic_clean(asg_t* read_g)
{
uint32_t n_vtx, v, u;
uint64_t i, k, k_i, tLen, v_occ, u_occ, utg_occ;
double bub_rate = 0.1;
ma_ug_t *ug = NULL;
ug = ma_ug_gen_primary(read_g, PRIMARY_LABLE);
n_vtx = ug->g->n_seq * 2;
buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t));
///for (i = 0, tLen = 1; i < ug->u.n; i++) tLen += ug->u.a[i].len;
tLen = get_bub_pop_max_dist_advance(ug->g, &b);
uint8_t* bs_flag = (uint8_t*)calloc(n_vtx, 1);
kvec_t(uint32_t) ax;
kv_init(ax);
for (v = 0; v < ug->g->n_seq; ++v)
{
if(ug->g->seq[v].del) continue;
ug->g->seq[v].c = PRIMARY_LABLE;
EvaluateLen(ug->u, v) = ug->u.a[v].n;
}
for (v = 0; v < n_vtx; ++v)
{
if(ug->g->seq[v>>1].del) continue;
if(asg_arc_n(ug->g, v) < 2) continue;
if(bs_flag[v] != 0) continue;
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0))
{
//beg is v, end is b.S.a[0]
//note b.b include end, does not include beg
for (i = 0; i < b.b.n; i++)
{
if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue;
bs_flag[b.b.a[i]] = bs_flag[b.b.a[i]^1] = 1;
}
bs_flag[v] = 2; bs_flag[b.S.a[0]^1] = 3;
}
}
for (v = 0; v < n_vtx; ++v)
{
if(bs_flag[v] !=2) continue;
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0))
{
//note b.b include end, does not include beg
for (i = v_occ = ax.n = 0; i < b.b.n; i++)
{
if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue;
v_occ += ug->u.a[b.b.a[i]>>1].n;
kv_push(uint32_t, ax, b.b.a[i]>>1);
}
for (i = 0; i < ax.n; i++)
{
for (k = 0; k < 2; k++)
{
u = (ax.a[i]<<1) + k;
if(asg_arc_n(ug->g, u) < 2) continue;
if(asg_bub_pop1_primary_trio(ug->g, NULL, u, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0))
{
for (k_i = u_occ = utg_occ = 0; k_i < b.b.n; k_i++)
{
if(b.b.a[k_i]==u || b.b.a[k_i]==b.S.a[0]) continue;
u_occ += ug->u.a[b.b.a[k_i]>>1].n;
utg_occ++;
}
if(u_occ >= v_occ*bub_rate) continue;
if(u_occ > 3) continue;
if(utg_occ > 2) continue;
asg_bub_pop1_primary_trio(ug->g, NULL, u, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0);
}
}
}
}
}
ma_utg_t* m = NULL;
for (v = 0; v < ug->g->n_seq; ++v)
{
if(ug->g->seq[v].del) continue;
if(ug->g->seq[v].c != ALTER_LABLE) continue;
m = &(ug->u.a[v]);
if(m->m == 0) continue;
for (k = 0; k < m->n; k++)
{
asg_seq_del(read_g, m->a[k]>>33);
}
}
asg_cleanup(read_g);
free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); free(bs_flag);
ma_ug_destroy(ug);
kv_destroy(ax);
}
void write_trans_chain(trans_chain* t_ch, const char *fn)
{
char *buf = (char*)calloc(strlen(fn) + 25, 1);
@@ -12899,7 +12801,6 @@ long long tipsLen, float tip_drop_ratio, long long stops_threshold,
R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp,
long long gap_fuzz, bub_label_t* b_mask_t)
{
hic_clean(sg);
ug_opt_t opt; memset(&opt, 0, sizeof(opt));
opt.coverage_cut = coverage_cut;
opt.sources = sources;
@@ -13675,8 +13576,6 @@ long long tipsLen, float tip_drop_ratio, long long stops_threshold,
R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp,
bub_label_t* b_mask_t)
{
hic_clean(sg);
kvec_asg_arc_t_warp new_rtg_edges;
kv_init(new_rtg_edges.a);
ma_ug_t *ug = NULL;
@@ -29782,7 +29681,7 @@ void flat_bubbles(asg_t *sg, uint8_t* r_het)
if(bs_flag[v] == 1) continue;
if(bs_flag[v] == 0) bs_flag[v] = 1;
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0))
if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0))
{
//note b.b include end, does not include beg
for (i = path = 0; i < b.b.n; i++)
@@ -29871,7 +29770,7 @@ void flat_bubbles(asg_t *sg, uint8_t* r_het)
if(is_het_b > path && is_het_s > path && (is_het_b+is_het_s)>(path<<2))
{
asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0);
asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0);
n_pop++;
}
}
@@ -29962,7 +29861,7 @@ void flat_bubbles_advance(asg_t *sg, ma_hit_t_alloc* sources, R_to_U* ruIndex, u
if(bs_flag[v] == 1) continue;
if(bs_flag[v] == 0) bs_flag[v] = 1;
if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0))
if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0))
{
bs_flag[v] = 2; bs_flag[b.S.a[0]^1] = 2;
//note b.b include end, does not include beg
@@ -30015,7 +29914,7 @@ void flat_bubbles_advance(asg_t *sg, ma_hit_t_alloc* sources, R_to_U* ruIndex, u
if((C_bases/R_bases) <= het_thres)
{
fprintf(stderr, "s-utg%.6ul\te-utg%.6ul\tC_bases:%lu\tR_bases:%lu\n", (v>>1)+1, (b.S.a[0]>>1)+1, C_bases, R_bases);
asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0);
asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL, 0, 0);
n_pop++;
}