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

fix distance bug

Browse Source
This commit is contained in:
chhylp123
2021-04-13 01:18:49 -04:00
parent 36afbce9bc
commit 67c7218264
6 changed files with 677 additions and 1007 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Overlaps.cpp
+3 -2
View File
@@ -13222,7 +13222,7 @@ void clean_u_trans_t_idx(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g)
kt_u_trans_t_idx(ta, ug->g->n_seq);
kt_u_trans_t_symm(ta, ug);
filter_u_trans_t(ta, ug, read_g, 3);
debug_u_trans_t(ta);
///debug_u_trans_t(ta);
}
@@ -13256,8 +13256,9 @@ bub_label_t* b_mask_t)
clean_u_trans_t_idx(&(cov->t_ch->k_trans), ug, sg);
// print_untig_by_read(copy_ug, "m64011_190830_220126/88867583/ccs", 603738, NULL, NULL, "sb");
set_trio_flag_by_cov(ug, sg, cov);
///debug_gfa_space(ug, cov);
set_trio_flag_by_cov(ug, sg, cov);
// print_r_het(cov, R_INF.trio_flag, "out-1");
Overlaps.h
+1
View File
@@ -1180,6 +1180,7 @@ void chain_origin_trans_uid_by_distance(hap_cov_t *cov, asg_t *read_sg,
uint32_t *pri_a, uint32_t pri_n, uint32_t pri_beg, uint64_t *i_pri_len,
uint32_t *aux_a, uint32_t aux_n, uint32_t aux_beg, uint64_t *i_aux_len,
ma_ug_t *ug, uint32_t flag, double overall_score, const char* cmd);
int asg_arc_del_trans(asg_t *g, int fuzz);
#define JUNK_COV 5
#define DISCARD_RATE 0.8
Purge_Dups.cpp
-313
View File
@@ -4846,179 +4846,6 @@ void remove_contained_haplotig(hap_overlaps_list* all_ovlp, ma_ug_t *ug, asg_t*
// }
}
void purge_dups_back(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density,
uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain,
uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
{
asg_t *purge_g = NULL;
purge_g = asg_init();
asg_t* nsg = ug->g;
uint32_t v, rId, uId, i, offset;
ma_utg_t* reads = NULL;
uint64_t* position_index = NULL;
if(cov) position_index = cov->pos_idx;
else position_index = (uint64_t*)malloc(sizeof(uint64_t)*read_g->n_seq);
memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);
hap_overlaps_list all_ovlp;
init_hap_overlaps_list(&all_ovlp, nsg->n_seq);
hap_overlaps_list back_all_ovlp;
init_hap_overlaps_list(&back_all_ovlp, nsg->n_seq);
///uint32_t junk_cov, hap_cov, dip_cov, junk_occ, repeat_occ, single_cov;
asg_arc_t t, *p = NULL;
int r;
hap_alignment_struct_pip hap_buf;
long long k_mer_only, coverage_only;
if(asm_opt.hom_global_coverage != -1)
{
hap_buf.cov_threshold = asm_opt.hom_global_coverage;
}
else
{
hap_buf.cov_threshold = get_read_coverage_thres(ug, read_g, ruIndex, position_index,
sources, coverage_cut, read_g->n_seq, COV_COUNT, &k_mer_only, &coverage_only);
}
for (v = 0; v < nsg->n_seq; v++)
{
uId = v;
if(nsg->seq[uId].del || nsg->seq[uId].c == ALTER_LABLE)
{
asg_seq_set(purge_g, uId, 0, 1);
purge_g->seq[uId].c = ALTER_LABLE;
continue;
}
reads = &(ug->u.a[uId]);
for (i = 0, offset = 0; i < reads->n; i++)
{
rId = reads->a[i]>>33;
set_R_to_U(ruIndex, rId, uId, 1, &(read_g->seq[rId].c));
position_index[rId] = offset;
position_index[rId] = position_index[rId] << 32;
position_index[rId] = position_index[rId] | (uint64_t)i;
offset += (uint32_t)reads->a[i];
}
asg_seq_set(purge_g, uId, offset, 0);
purge_g->seq[uId].c = PRIMARY_LABLE;
}
init_hap_alignment_struct_pip(&hap_buf, asm_opt.thread_num, nsg->n_seq, ug, read_g,
sources, reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp,
0.1, &all_ovlp, cov);
if(hap_buf.cov_threshold < 0)
{
if(if_ploid_sample(ug, read_g, ruIndex, sources, reverse_sources, coverage_cut,
&hap_buf, &all_ovlp, &back_all_ovlp, purege_minLen, 0.333))
{
///if peak is het, coverage peak is more reliable
hap_buf.cov_threshold = coverage_only * HET_PEAK_RATE;
}
else
{
///if peak is homo, k-mer peak is more reliable
hap_buf.cov_threshold = k_mer_only * HOM_PEAK_RATE;
}
}
if(asm_opt.hom_global_coverage == -1) asm_opt.hom_global_coverage = hap_buf.cov_threshold;
fprintf(stderr, "[M::%s] purge duplication coverage threshold: %lld\n", __func__, hap_buf.cov_threshold);
if(just_coverage) goto end_coverage;
kt_for(asm_opt.thread_num, hap_alignment_advance_worker, &hap_buf, nsg->n_seq);
///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp);
filter_hap_overlaps_by_length(&all_ovlp, purege_minLen);
///normalize_hap_overlaps(&all_ovlp, &back_all_ovlp);
normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
///debug_hap_overlaps(&all_ovlp, &back_all_ovlp);
remove_contained_haplotig(&all_ovlp, ug, nsg, purge_g, cov);
if(just_contain == 0)
{
for (v = 0; v < all_ovlp.num; v++)
{
uId = v;
if(purge_g->seq[uId].del || purge_g->seq[uId].c == ALTER_LABLE) continue;
for (i = 0; i < all_ovlp.x[uId].a.n; i++)
{
if(all_ovlp.x[uId].a.a[i].status == DELETE) continue;
///if(all_ovlp.x[uId].a.a[i].type == )
if(purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].c == ALTER_LABLE||
purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].del||
purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].c == ALTER_LABLE||
purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].del)
{
continue;
}
///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
r = get_hap_arch(&(all_ovlp.x[uId].a.a[i]), ug->u.a[all_ovlp.x[uId].a.a[i].xUid].len,
ug->u.a[all_ovlp.x[uId].a.a[i].yUid].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
// if(all_ovlp.x[uId].a.a[i].xUid == 118 && all_ovlp.x[uId].a.a[i].yUid == 82)
// {
// fprintf(stderr, "r: %d\n", r);
// print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
// }
if(r < 0) continue;
p = asg_arc_pushp(purge_g);
*p = t;
}
}
asg_cleanup(purge_g);
asg_symm(purge_g);
///may need to do transitive reduction
clean_purge_graph(purge_g, drop_ratio, 1);
// if(debug_enable) print_purge_gfa(ug, purge_g);
// if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp);
link_unitigs(purge_g, ug, &all_ovlp, ruIndex, reverse_sources, coverage_cut, read_g, position_index,
&(hap_buf.buf[0].u_buffer), &(hap_buf.buf[0].u_buffer_tailIndex), &(hap_buf.buf[0].u_buffer_prevIndex),
max_hang, min_ovlp, edge, hap_buf.buf[0].visit, cov);
}
for (v = 0; v < all_ovlp.num; v++)
{
uId = v;
if(purge_g->seq[uId].c == ALTER_LABLE)
{
ug->g->seq[uId].c = ALTER_LABLE;
}
}
end_coverage:
uint32_t is_Unitig;
for (v = 0; v < ruIndex->len; v++)
{
get_R_to_U(ruIndex, v, &uId, &is_Unitig);
if(is_Unitig == 1) ruIndex->index[v] = (uint32_t)-1;
}
asg_cleanup(nsg);
destory_hap_overlaps_list(&all_ovlp);
destory_hap_overlaps_list(&back_all_ovlp);
asg_destroy(purge_g);
if(cov) memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);
else free(position_index);
destory_hap_alignment_struct_pip(&hap_buf);
}
void debug_p_g_t(p_g_t* pg, hap_cov_t *cov, asg_t *read_g)
{
fprintf(stderr, "----------[M::%s]----------\n", __func__);
@@ -5229,146 +5056,6 @@ void destory_p_g_t(p_g_t **pg)
}
void partition_contigs(hap_overlaps_list* all_ovlp, ma_ug_t *ug, hap_cov_t *cov, double keep_rate,
int max_hang, int min_ovlp, float drop_ratio, p_g_t *pg)
{
int r, index;
uint32_t v, i, uId, m;
hap_overlaps *p = NULL;
asg_arc_t t, *p_t = NULL;
for (v = 0; v < all_ovlp->num; v++)
{
uId = v;
///has been removed as contained
if(pg->pg_h_lev->seq[uId].del || pg->pg_h_lev->seq[uId].c == ALTER_LABLE) continue;
for (i = 0; i < all_ovlp->x[uId].a.n; i++)
{
if(all_ovlp->x[uId].a.a[i].type == YCX) continue;
if(all_ovlp->x[uId].a.a[i].type == XCY) continue;
/****************************may have bugs********************************/
if(all_ovlp->x[uId].a.a[i].score <= 0) continue;
/****************************may have bugs********************************/
if(pg->pg_h_lev->seq[all_ovlp->x[uId].a.a[i].xUid].c == ALTER_LABLE||
pg->pg_h_lev->seq[all_ovlp->x[uId].a.a[i].xUid].del||
pg->pg_h_lev->seq[all_ovlp->x[uId].a.a[i].yUid].c == ALTER_LABLE||
pg->pg_h_lev->seq[all_ovlp->x[uId].a.a[i].yUid].del)
{
continue;
}
///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
r = get_hap_arch(&(all_ovlp->x[uId].a.a[i]), ug->u.a[all_ovlp->x[uId].a.a[i].xUid].len,
ug->u.a[all_ovlp->x[uId].a.a[i].yUid].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
if(r < 0) continue;
p_t = asg_arc_pushp(pg->pg_h_lev);
*p_t = t;
}
}
asg_cleanup(pg->pg_h_lev);
asg_symm(pg->pg_h_lev);
clean_purge_graph(pg->pg_h_lev, keep_rate, 0);
asg_arc_t *av = NULL;
uint32_t n_vtx = pg->pg_h_lev->n_seq<<1, nv, a, b;
for (v = 0; v < n_vtx; v++)
{
av = asg_arc_a(pg->pg_h_lev, v);
nv = asg_arc_n(pg->pg_h_lev, v);
for (i = 0; i < nv; ++i)
{
if(av[i].del) continue;
a = av[i].ul>>33;
b = av[i].v>>1;
index = get_specific_hap_overlap(&(all_ovlp->x[a]), a, b);
p = &(all_ovlp->x[a].a.a[index]);
p->status = MIXED;
index = get_specific_hap_overlap(&(all_ovlp->x[b]), b, a);
p = &(all_ovlp->x[b].a.a[index]);
p->status = MIXED;
}
}
for (v = 0; v < all_ovlp->num; v++)
{
uId = v;
for (i = m = 0; i < all_ovlp->x[uId].a.n; i++)
{
p = (&all_ovlp->x[uId].a.a[i]);
/****************************may have bugs********************************/
if(p->score <= 0) continue;
/****************************may have bugs********************************/
if(p->type == YCX || p->type == XCY || p->status == MIXED)
{
all_ovlp->x[uId].a.a[m] = (*p);
all_ovlp->x[uId].a.a[m].status = SELF_EXIST;
m++;
}
}
all_ovlp->x[uId].a.n = m;
}
for (v = 0; v < all_ovlp->num; v++)
{
uId = v; p = NULL;
if(all_ovlp->x[uId].a.n == 0) continue;
for (i = 0; i < all_ovlp->x[uId].a.n; i++)
{
if(p == NULL || p->score < all_ovlp->x[uId].a.a[i].score)
{
p = &(all_ovlp->x[uId].a.a[i]);
}
}
if(!p) continue;
for (i = m = 0; i < all_ovlp->x[uId].a.n; i++)
{
if(all_ovlp->x[uId].a.a[i].type == YCX || all_ovlp->x[uId].a.a[i].type == XCY)
{
if(filter_secondary_chain(p->score, all_ovlp->x[uId].a.a[i].score, keep_rate))
{
all_ovlp->x[uId].a.a[m] = all_ovlp->x[uId].a.a[i];
m++;
}
}
else
{
all_ovlp->x[uId].a.a[m] = all_ovlp->x[uId].a.a[i];
m++;
}
}
all_ovlp->x[uId].a.n = m;
}
for (v = 0; v < all_ovlp->num; v++)
{
uId = v;
for (i = m = 0; i < all_ovlp->x[uId].a.n; i++)
{
p = (&all_ovlp->x[uId].a.a[i]);
index = get_specific_hap_overlap(&(all_ovlp->x[p->yUid]), p->yUid, p->xUid);
if(index != -1)
{
all_ovlp->x[uId].a.a[m] = (*p);
m++;
}
}
all_ovlp->x[uId].a.n = m;
}
}
void chain_origin_trans_uid_by_purge(hap_overlaps *x, ma_ug_t *ug, hap_cov_t *cov, uint64_t* position_index)
{
uint32_t pri_uid, aux_uid, r_x, r_y;
Purge_Dups.h
+3
View File
@@ -85,5 +85,8 @@ uint32_t classify_hap_overlap(long long xBeg, long long xEnd, long long xLen,
long long yBeg, long long yEnd, long long yLen, long long* r_xBeg, long long* r_xEnd,
long long* r_yBeg, long long* r_yEnd);
int cmp_hap_alignment_chaining(const void * a, const void * b);
uint32_t classify_hap_overlap(long long xBeg, long long xEnd, long long xLen,
long long yBeg, long long yEnd, long long yLen, long long* r_xBeg, long long* r_xEnd,
long long* r_yBeg, long long* r_yEnd);
#endif
hic.cpp
+669 -692
View File
File diff suppressed because it is too large Load Diff
hic.h
+1
View File
@@ -69,5 +69,6 @@ void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, u
void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint32_t is_end);
void set_b_utg_weight_flag(bubble_type* bub, buf_t* b, uint32_t v, uint8_t* vis_flag, uint32_t flag, uint32_t* occ);
uint32_t get_unitig_het_arb(ma_utg_t* u, uint8_t *r_het_flag, uint32_t m_het_label, uint32_t p_het_label, uint32_t n_het_label);
void debug_gfa_space(ma_ug_t* ug, hap_cov_t *cov);
#endif