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

bubble label

Browse Source
This commit is contained in:
chhylp123
2021-03-14 03:25:19 -04:00
parent 4c2ce6fc6b
commit 8e4b98f0a6
5 changed files with 1161 additions and 406 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Overlaps.cpp
+1097 -318
View File
File diff suppressed because it is too large Load Diff
Overlaps.h
+40 -65
View File
@@ -551,55 +551,6 @@ inline uint32_t check_tip(asg_t *sg, uint32_t begNode, uint32_t* endNode, buf_t*
}
}
inline uint32_t get_unitig_back(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode,
long long* nodeLen, long long* baseLen, buf_t* b)
{
ma_utg_v* u = NULL;
uint32_t v = begNode, w, k;
uint32_t kv;
(*nodeLen) = (*baseLen) = 0;
(*endNode) = (uint32_t)-1;
if(ug!=NULL) u = &(ug->u);
while (1)
{
kv = get_real_length(sg, v, NULL);
(*endNode) = v;
if(u == NULL)
{
(*nodeLen)++;
}
else
{
(*nodeLen) += EvaluateLen((*u), v>>1);
}
if(b) kv_push(uint32_t, b->b, v);
///means reach the end of a unitig
if(kv!=1) (*baseLen) += sg->seq[v>>1].len;
if(kv==0) return END_TIPS;
if(kv>1) return MUL_OUTPUT;
///kv must be 1 here
kv = get_real_length(sg, v, &w);
///means reach the end of a unitig
if(get_real_length(sg, w^1, NULL)!=1)
{
(*baseLen) += sg->seq[v>>1].len;
return MUL_INPUT;
}
for (k = 0; k < asg_arc_n(sg, v); k++)
{
if(asg_arc_a(sg, v)[k].del) continue;
///here is just one undeleted edge
(*baseLen) += asg_arc_len(asg_arc_a(sg, v)[k]);
break;
}
v = w;
if(v == begNode) return LOOP;
}
}
inline uint32_t get_unitig(asg_t *sg, ma_ug_t *ug, uint32_t begNode, uint32_t* endNode,
long long* nodeLen, long long* baseLen, long long* max_stop_nodeLen, long long* max_stop_baseLen,
uint32_t stops_threshold, buf_t* b)
@@ -1030,25 +981,47 @@ typedef struct {
uint32_t total;
} Trio_counter;
typedef struct {
uint32_t p; // the optimal parent vertex
uint32_t d; // the shortest distance from the initial vertex
uint32_t r:31, s:1; // r: the number of remaining incoming arc; s: state
} binfo_s_t;
typedef struct {
///all information for each node
binfo_s_t *a;
kvec_t(uint32_t) S; // set of vertices without parents, nodes with all incoming edges visited
kvec_t(uint32_t) b; // visited vertices
kvec_t(uint32_t) e; // visited edges/arcs
} buf_s_t;
typedef struct{
buf_s_t *b;
uint32_t n_thres, n_reads;
asg_t *g;
uint32_t check_cross;
uint64_t bub_dist;
} bub_label_t;
void resolve_tangles(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig,
long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t trio_flag,
float drop_ratio);
void adjust_utg_advance(asg_t *sg, ma_ug_t *ug, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
void rescue_contained_reads_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t chainLenThres, uint32_t is_bubble_check,
uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* new_rtg_nodes);
uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* new_rtg_nodes, bub_label_t* b_mask_t);
void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t is_bubble_check,
uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges);
uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t);
void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut,
ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate);
void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex);
void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges);
ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges, bub_label_t* b_mask_t);
void get_unitig_trio_flag(ma_utg_t* nsu, uint32_t flag, uint32_t* require, uint32_t* non_require, uint32_t* ambigious);
void rescue_missing_overlaps_backward(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t backward_steps,
uint32_t is_bubble_check, uint32_t is_primary_check);
uint32_t is_bubble_check, uint32_t is_primary_check, bub_label_t* b_mask_t);
uint32_t get_edge_from_source(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t query, uint32_t target, asg_arc_t* t);
int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio);
@@ -1085,12 +1058,15 @@ typedef struct{
uint64_t r_num;
} hc_links;
typedef struct{
kvec_t(uint32_t) uIDs;
kvec_t(uint32_t) idx;
uint32_t chain_num;
kvec_t(uint32_t) iDXs;
kvec_t(uint32_t) rescue_hom;
uint32_t* u_idx;
uint64_t r_num;
uint32_t r_num;
uint32_t chain_num;
uint32_t l0_chain, l1_chain;
}trans_chain;
typedef struct {
@@ -1106,8 +1082,8 @@ typedef struct {
kvec_asg_arc_t_offset u_buffer;
kvec_t_i32_warp tailIndex;
kvec_t_i32_warp prevIndex;
hc_links* link;
///trans_chain t_ch;
///hc_links* link;
trans_chain* t_ch;
}hap_cov_t;
typedef struct{
@@ -1118,28 +1094,27 @@ typedef struct{
void init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num);
void destory_hc_links(hc_links* link);
int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov);
uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, int max_dist, buf_t *b, uint32_t positive_flag,
uint64_t get_bub_pop_max_dist(asg_t *g, buf_t *b);
uint64_t get_bub_pop_max_dist_advance(asg_t *g, buf_t *b);
int asg_pop_bubble_primary_trio(ma_ug_t *ug, uint64_t* i_max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov);
uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, uint64_t max_dist, buf_t *b, uint32_t positive_flag,
uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes, hap_cov_t *cov);
void adjust_utg_by_primary(ma_ug_t **ug, asg_t* read_g, float drop_rate,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut,
long long bubble_dist, 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,
kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link);
kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link, bub_label_t* b_mask_t);
void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug, asg_t *read_sg);
ma_ug_t* copy_untig_graph(ma_ug_t *src);
ma_ug_t* output_trio_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name,
uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist,
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, int is_bench);
float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench, bub_label_t* b_mask_t);
asg_t* copy_read_graph(asg_t *src);
ma_ug_t *ma_ug_gen(asg_t *g);
void ma_ug_destroy(ma_ug_t *ug);
inline int inter_interval(int a_s, int a_e, int b_s, int b_e, int* i_s, int* i_e)
{
if(a_s > b_e || b_s > a_e) return 0;
Purge_Dups.cpp
+19 -18
View File
@@ -2575,11 +2575,12 @@ long long* r_y_pos_beg, long long* r_y_pos_end)
void print_hap_paf(ma_ug_t *ug, hap_overlaps* ovlp)
{
fprintf(stderr, "utg%.6d%c\t%u(%u)\t%u(%u)\t%u(%u)\t%c\tutg%.6d%c\t%u(%u)\t%u(%u)\t%u(%u)\t%u\t%u\n",
fprintf(stderr, "utg%.6d%c\t%u(%u)\t%u(%u)\t%u(%u)\t%c\tutg%.6d%c\t%u(%u)\t%u(%u)\t%u(%u)\t%u\t%u\t%lld(%u)\n",
ovlp->xUid+1, "lc"[ug->u.a[ovlp->xUid].circ], ug->u.a[ovlp->xUid].len, ug->u.a[ovlp->xUid].n,
ovlp->x_beg_pos, ovlp->x_beg_id, ovlp->x_end_pos, ovlp->x_end_id, "+-"[ovlp->rev],
ovlp->yUid+1, "lc"[ug->u.a[ovlp->yUid].circ], ug->u.a[ovlp->yUid].len, ug->u.a[ovlp->yUid].n,
ovlp->y_beg_pos, ovlp->y_beg_id, ovlp->y_end_pos, ovlp->y_end_id, ovlp->type, (uint32_t)ovlp->weight);
ovlp->y_beg_pos, ovlp->y_beg_id, ovlp->y_end_pos, ovlp->y_end_id, ovlp->type, ovlp->weight,
ovlp->score, ovlp->status);
}
inline long long get_max_index(asg_arc_t_offset* x, int32_t* Scores, uint8_t* Flag, long long n,
@@ -3148,6 +3149,7 @@ void set_reverse_hap_overlap(hap_overlaps* dest, hap_overlaps* source, uint32_t*
dest->x_end_id = source->y_end_id;
dest->y_beg_id = source->x_beg_id;
dest->y_end_id = source->x_end_id;
dest->score = source->score;
}
/**
@@ -3409,13 +3411,13 @@ uint64_t asg_bub_pop1_purge_graph(asg_t *g, uint32_t v0, int max_dist, buf_t *b)
///assert(nv > 0);
///all out-edges of v
for (i = 0; i < nv; ++i) { // loop through v's neighbors
///if this edge has been deleted
if (av[i].del) continue;
uint32_t w = av[i].v; // v->w with length l
binfo_t *t = &b->a[w];
///that means there is a circle, directly terminate the whole bubble poping
///if (w == v0) goto pop_reset;
if ((w>>1) == (v0>>1)) goto pop_reset;
///if this edge has been deleted
if (av[i].del) continue;
c_s = decode_score((uint32_t)av[i].ul, av[i].ol);
///push the edge
///high 32-bit of g->idx[v] is the start point of v's edges
@@ -3499,7 +3501,7 @@ pop_reset:
// pop bubbles
int asg_pop_bubble_purge_graph(asg_t *purge_g, int max_dist)
int asg_pop_bubble_purge_graph(asg_t *purge_g)
{
uint32_t v, n_vtx = purge_g->n_seq * 2;
uint64_t n_pop = 0;
@@ -3641,13 +3643,13 @@ int purge_g_arc_del_short_diploid_by_score(asg_t *g, float drop_ratio)
}
void clean_purge_graph(asg_t *purge_g, int max_dist, float drop_ratio)
void clean_purge_graph(asg_t *purge_g, float drop_ratio)
{
uint64_t operation = 1;
while (operation > 0)
{
operation = 0;
operation += asg_pop_bubble_purge_graph(purge_g, max_dist);
operation += asg_pop_bubble_purge_graph(purge_g);
operation += purge_g_arc_del_short_diploid_by_score(purge_g, drop_ratio);
}
@@ -4122,7 +4124,7 @@ hap_cov_t *cov)
}
purge_g->seq[w>>1].c = ALTER_LABLE;
collect_trans_purge_joint_cov(cov, ug, x);
if(cov) collect_trans_purge_joint_cov(cov, ug, x);
// if(buffer.n > 1)
// {
@@ -4263,20 +4265,21 @@ hap_cov_t *cov)
continue;
}
if(cov->link) collect_reverse_unitigs_purge(&b_0, cov->link, ug, all_ovlp);
///if(cov->link) collect_reverse_unitigs_purge(&b_0, cov->link, ug, all_ovlp);
purge_merge(purge_g, ug, all_ovlp, &b_0, ruIndex, reverse_sources, coverage_cut,
read_g, position_index, u_buffer, tailIndex, prevIndex,max_hang, min_ovlp, edge, visit, cov);
}
free(b_0.b.a);
}
void print_all_purge_ovlp(ma_ug_t *ug, hap_overlaps_list* all_ovlp)
void print_all_purge_ovlp(ma_ug_t *ug, hap_overlaps_list* all_ovlp, const char* cmd)
{
fprintf(stderr, "\n%s--->ug->u.n: %u\n", cmd, (uint32_t)ug->u.n);
uint32_t v, uId, i;
for (v = 0; v < all_ovlp->num; v++)
{
uId = v;
if(uId != 96 && uId != 272) continue;
///if(uId != 96 && uId != 272) continue;
for (i = 0; i < all_ovlp->x[uId].a.n; i++)
{
print_hap_paf(ug, &(all_ovlp->x[uId].a.a[i]));
@@ -4571,7 +4574,7 @@ void remove_contained_haplotig(hap_overlaps_list* all_ovlp, ma_ug_t *ug, asg_t*
purge_g->seq[xUid].del = 1;
all_ovlp->x[uId].a.a[i].status = DELETE;
if(cov->link) collect_reverse_unitig_pair(cov->link, ug, &(all_ovlp->x[uId].a.a[i]));
///if(cov->link) collect_reverse_unitig_pair(cov->link, ug, &(all_ovlp->x[uId].a.a[i]));
collect_trans_purge_cov(cov, ug, &(all_ovlp->x[uId].a.a[i]), 0);
}
@@ -4612,8 +4615,8 @@ void remove_contained_haplotig(hap_overlaps_list* all_ovlp, ma_ug_t *ug, asg_t*
void purge_dups(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, long long bubble_dist, float drop_ratio,
uint32_t just_contain, uint32_t just_coverage, hap_cov_t *cov)
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)
{
asg_t *purge_g = NULL;
purge_g = asg_init();
@@ -4697,17 +4700,15 @@ uint32_t just_contain, uint32_t just_coverage, hap_cov_t *cov)
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++)
@@ -4747,7 +4748,7 @@ uint32_t just_contain, uint32_t just_coverage, hap_cov_t *cov)
asg_cleanup(purge_g);
asg_symm(purge_g);
///may need to do transitive reduction
clean_purge_graph(purge_g, bubble_dist, drop_ratio);
clean_purge_graph(purge_g, drop_ratio);
// if(debug_enable) print_purge_gfa(ug, purge_g);
// if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp);
Purge_Dups.h
+3 -3
View File
@@ -15,14 +15,14 @@
void purge_dups(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, long long bubble_dist, float drop_ratio,
uint32_t just_contain, uint32_t just_coverage, hap_cov_t *cov);
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);
void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge,
uint32_t is_circle, uint64_t* rLen);
void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen);
void enable_debug_mode(uint32_t mode);
hap_cov_t* init_hap_cov_t(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* sources, R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources,
ma_sub_t *coverage_cut, int max_hang, int min_ovlp, hc_links* link);
ma_sub_t *coverage_cut, int max_hang, int min_ovlp, uint32_t is_collect_trans);
void destory_hap_cov_t(hap_cov_t **x);
void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd);
hic.cpp
+2 -2
View File
@@ -2566,14 +2566,14 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, hc_links* link)
if (!ug->g->is_symm) asg_symm(ug->g);
uint32_t v, n_vtx = ug->g->n_seq * 2, i, k, mode = (((uint32_t)-1)<<2);
uint32_t beg, sink, n, *a, n_occ;
uint64_t pathLen, tLen;
uint64_t pathLen;
bub->ug = ug;
for (i = 0, tLen = 1; i < ug->u.n; i++) tLen += ug->u.a[i].len;
bub->b_bub = bub->b_end_bub = bub->tangle_bub = bub->cross_bub = bub->mess_bub = 0;
if(bub->round_id == 0)
{
buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t));
uint64_t tLen = get_bub_pop_max_dist_advance(ug->g, &b);
kv_init(bub->list); kv_init(bub->num); kv_init(bub->pathLen);
kv_init(bub->b_s_idx); kv_malloc(bub->b_s_idx, ug->g->n_seq);
bub->b_ug = NULL; kv_init(bub->chain_weight);