bug fix for bubble scanning

This commit is contained in:
chhylp123
2021-04-23 20:40:43 -04:00
parent a5b29b00bd
commit 7b07e355b7
8 changed files with 183 additions and 97 deletions
+2
View File
@@ -10,6 +10,7 @@
#include "Correct.h"
#include "htab.h"
#include "kthread.h"
#include "rcut.h"
void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, 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, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct);
@@ -1605,6 +1606,7 @@ void ha_overlap_final(void)
int ha_assemble(void)
{
// debug_mc_g_t(MC_NAME);
extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o);
int r, hom_cov = -1, ovlp_loaded = 0;
if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
+40 -53
View File
@@ -12736,7 +12736,7 @@ bub_label_t* b_mask_t)
///asm_opt.purge_simi_thres = asm_opt.purge_simi_rate_hic;
adjust_utg_by_primary(&copy_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut,
tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio,
max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1);
max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 1/**0**/);
print_utg(copy_ug, copy_sg, coverage_cut, output_file_name, sources, ruIndex, max_hang,
min_ovlp, &new_rtg_edges);
@@ -13466,14 +13466,12 @@ bub_label_t* b_mask_t)
hap_cov_t *cov = NULL;
asg_t *copy_sg = copy_read_graph(sg);
ma_ug_t *copy_ug = copy_untig_graph(ug);
ma_ug_t *copy_ug = copy_untig_graph(ug);
adjust_utg_by_primary(&copy_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut,
tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio,
max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1);
max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 1);
print_utg(copy_ug, copy_sg, coverage_cut, output_file_name, sources, ruIndex, max_hang,
min_ovlp, &new_rtg_edges);
ma_ug_destroy(copy_ug);
asg_destroy(copy_sg);
@@ -15318,7 +15316,6 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov)
{
pre_cons = get_graph_statistic(g);
asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, 1);
if(just_bubble_pop == 0)
{
///need consider tangles
@@ -16088,7 +16085,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t)
{
purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges,
asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 1, 0,
cov, 0);
cov, 0, 0);
///delete_useless_nodes(ug);
delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex);
}
@@ -16388,6 +16385,7 @@ void dfs_trans_chain_bub(asg_t *g, hap_cov_t *cov, uint32_t v, uint32_t beg, uin
{
b->b.n--;
cur = b->b.a[b->b.n];
if(flag[cur>>1] == 0 && (cur>>1) != (v>>1)) continue;
flag[cur>>1] = 0;
ncur = asg_arc_n(g, cur);
@@ -16396,6 +16394,7 @@ void dfs_trans_chain_bub(asg_t *g, hap_cov_t *cov, uint32_t v, uint32_t beg, uin
{
if(acur[i].del) continue;
if((acur[i].v>>1) == beg || (acur[i].v>>1) == sink) continue;
if(flag[acur[i].v>>1] == 0) continue;
kv_push(uint32_t, b->b, acur[i].v);
}
}
@@ -16406,6 +16405,7 @@ void dfs_trans_chain_bub(asg_t *g, hap_cov_t *cov, uint32_t v, uint32_t beg, uin
{
b->b.n--;
cur = b->b.a[b->b.n];
if(flag[cur>>1] == 0 && (cur>>1) != (v>>1)) continue;
flag[cur>>1] = 0;
ncur = asg_arc_n(g, cur);
@@ -16414,6 +16414,7 @@ void dfs_trans_chain_bub(asg_t *g, hap_cov_t *cov, uint32_t v, uint32_t beg, uin
{
if(acur[i].del) continue;
if((acur[i].v>>1) == beg || (acur[i].v>>1) == sink) continue;
if(flag[acur[i].v>>1] == 0) continue;
kv_push(uint32_t, b->b, acur[i].v);
}
}
@@ -16641,14 +16642,14 @@ void chain_origin_trans_uid_c_bubble(uint32_t query, buf_t *target, buf_t *idx,
// in a resolved bubble, mark unused vertices and arcs as "reduced"
static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, hap_cov_t *cov, uint32_t is_update_chain)
{
uint32_t i, k, k_i, v, u, uLen = 0, uCov = 0, uId, rId, ori;
uint32_t i, k, k_i, v, u, uLen = 0, uCov = 0, uId, rId, ori;
ma_utg_t* p = NULL;
trans_chain* t_ch = (is_update_chain?cov->t_ch:NULL);
///b->S.a[0] is the sink of this bubble
///assert(b->S.n == 1);
///first remove all nodes in this bubble
for (i = 0; i < b->b.n; ++i)
///assert(b->S.n == 1);
///first remove all nodes in this bubble
for (i = 0; i < b->b.n; ++i)
{
uId = b->b.a[i]>>1;
if(uId == (b->S.a[0]>>1)) continue;
@@ -16664,10 +16665,10 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
///v is the sink of this bubble
v = b->S.a[0];
///recover node
do {
u = b->a[v].p; // u->v
v = b->S.a[0];
///recover node
do {
u = b->a[v].p; // u->v
if(v != b->S.a[0])
{
uId = v>>1;
@@ -16680,16 +16681,16 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
uLen += cov->read_g->seq[rId].len;
}
}
v = u;
} while (v != v0);
v = u;
} while (v != v0);
uCov = (uLen == 0? 0 : uCov / uLen);
///v is the sink of this bubble
v = b->S.a[0];
///recover node
do {
u = b->a[v].p; // u->v
v = b->S.a[0];
///recover node
do {
u = b->a[v].p; // u->v
if(v != b->S.a[0])
{
uId = v>>1;
@@ -16701,8 +16702,8 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
cov->cov[rId] += (uCov * cov->read_g->seq[rId].len);
}
}
v = u;
} while (v != v0);
v = u;
} while (v != v0);
if(t_ch)
@@ -16763,7 +16764,6 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
}
}
topologicalSortUtil(ug->g, cov, v0, b->S.a[0]);
///if(cov->t_ch->topo_res.n != b->b.n - 1) fprintf(stderr, "ERROR-4\n");
if(cov->t_ch->topo_res.n == 0) return;
@@ -16773,8 +16773,10 @@ static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, ha
if(uId == (b->S.a[0]>>1)) continue;
dfs_trans_chain_bub(ug->g, cov, uId, v0>>1, b->S.a[0]>>1);
if(cov->t_ch->b_buf_0.b.n == 0) continue;
chain_origin_trans_uid_c_bubble(cov->t_ch->topo_res.a[i], &(t_ch->b_buf_0), b, ug, cov);
/***********************x***********************/
uId = cov->t_ch->topo_res.a[i]>>1;
p = &(ug->u.a[uId]);
@@ -17528,7 +17530,6 @@ hap_cov_t *cov, uint32_t is_update_chain, uint32_t keep_d)
/****************************may have bugs********************************/
///if(keep_d != 0) debug_asg_bub_pop1_primary_trio(g, utg, v0, max_dist, b, positive_flag, negative_flag, 1);
/****************************may have bugs********************************/
if(cov && utg) asg_bub_backtrack_primary_cov(utg, v0, b, cov, is_update_chain);
if(is_pop) asg_bub_backtrack_primary(g, v0, b);
if(path_base_len || path_nodes) asg_bub_backtrack_primary_length(g, utg, v0, b, path_base_len, path_nodes);
@@ -23678,7 +23679,7 @@ ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut
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, hap_cov_t **i_cov, bub_label_t* b_mask_t,
uint32_t collect_p_trans)
uint32_t collect_p_trans, uint32_t collect_p_trans_f)
{
asg_t* nsg = (*ug)->g;
uint32_t v, n_vtx = nsg->n_seq, k, rId, just_contain;
@@ -23686,9 +23687,8 @@ uint32_t collect_p_trans)
hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources,
coverage_cut, max_hang, min_ovlp, (asm_opt.purge_level_primary>0||i_cov)?1:0);
if(cov->t_ch) set_r_het_flag(*ug, read_g, coverage_cut, sources, ruIndex, cov->t_ch);
adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex, b_mask_t);
nsg = (*ug)->g;
n_vtx = nsg->n_seq;
for (v = 0; v < n_vtx; ++v)
@@ -23698,14 +23698,11 @@ uint32_t collect_p_trans)
EvaluateLen((*ug)->u, v) = (*ug)->u.a[v].n;
}
clean_primary_untig_graph(*ug, read_g, sources, reverse_sources, coverage_cut, tipsLen, tip_drop_ratio,
stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, chimeric_rate, 0, 0, drop_ratio, cov);
delete_useless_nodes(ug);
renew_utg(ug, read_g, new_rtg_edges);
if(i_cov && collect_p_trans == 0) goto skip_purge;
if(asm_opt.purge_level_primary > 0)
{
///print_debug_gfa(read_g, *ug, coverage_cut, "debug_purge", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
@@ -23713,7 +23710,7 @@ uint32_t collect_p_trans)
if(asm_opt.purge_level_primary == 1) just_contain = 1;
purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges,
asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio,
just_contain, 0, cov, !!(cov->t_ch&&collect_p_trans));
just_contain, 0, cov, !!(cov->t_ch&&collect_p_trans), collect_p_trans_f);
delete_useless_nodes(ug);
renew_utg(ug, read_g, new_rtg_edges);
}
@@ -23855,26 +23852,21 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov
kvec_asg_arc_t_warp new_rtg_edges;
kv_init(new_rtg_edges.a);
adjust_utg_by_primary(&ug, sg, TRIO_THRES, sources, reverse_sources, coverage_cut,
tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio,
max_hang, min_ovlp, &new_rtg_edges, NULL, b_mask_t, 0);
max_hang, min_ovlp, &new_rtg_edges, NULL, b_mask_t, 0, 0);
if(asm_opt.b_low_cov > 0)
{
break_ug_contig(&ug, sg, &R_INF, coverage_cut, sources, ruIndex, &new_rtg_edges, max_hang, min_ovlp,
&asm_opt.b_low_cov, NULL, asm_opt.m_rate);
}
if(asm_opt.b_high_cov > 0)
{
break_ug_contig(&ug, sg, &R_INF, coverage_cut, sources, ruIndex, &new_rtg_edges, max_hang, min_ovlp,
NULL, &asm_opt.b_high_cov, asm_opt.m_rate);
}
ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0);
fprintf(stderr, "Writing primary contig GFA to disk... \n");
@@ -27096,8 +27088,8 @@ void reset_bub(bubble_type* bub, ma_ug_t *ug, trans_chain* back_ug_chain, kvec_a
new_rtg_edges->a.n = 0;
///classify_untigs(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, max_hang, min_ovlp);
identify_bubbles(ug, bub, back_ug_chain->is_r_het, NULL);
update_bubble_chain(ug, bub, 0, 1);
resolve_bubble_chain_tangle(ug, bub);
// update_bubble_chain(ug, bub, 0, 1);
// resolve_bubble_chain_tangle(ug, bub);
// fprintf(stderr, "bub.f_bub: %lu, bub.b_bub: %lu, bub.b_end_bub: %lu, bub.tangle_bub: %lu, bub.cross_bub: %lu\n",
// bub->f_bub, bub->b_bub, bub->b_end_bub, bub->tangle_bub, bub->cross_bub);
}
@@ -27396,16 +27388,20 @@ bub_label_t* b_mask_t)
ma_ug_t *ug = NULL;
ug = ma_ug_gen_primary(sg, PRIMARY_LABLE);
// FILE* output_file = fopen("straw-debug.noseq.gfa", "w");
// ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file);
// fclose(output_file);
hap_cov_t *cov = NULL;
asg_t *copy_sg = copy_read_graph(sg);
ma_ug_t *copy_ug = copy_untig_graph(ug);
adjust_utg_by_primary(&copy_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut,
tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio,
max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 0);
max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 0, 0);
ma_ug_destroy(copy_ug); copy_ug = NULL;
asg_destroy(copy_sg); copy_sg = NULL;
uint32_t beg_idx, occ;
bubble_type bub;
memset(&bub, 0, sizeof(bubble_type));
@@ -29276,7 +29272,7 @@ ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, int m
hap_cov_t *cov = init_hap_cov_t(ug, sg, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, 0);
purge_dups(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, NULL,
opt->purge_simi_thres, opt->purge_overlap_len, max_hang, min_ovlp, 0, 0, 1, cov, 0);
opt->purge_simi_thres, opt->purge_overlap_len, max_hang, min_ovlp, 0, 0, 1, cov, 0, 0);
destory_hap_cov_t(&cov);
ma_ug_destroy(ug);
@@ -29562,7 +29558,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
init_bub_label_t(&b_mask_t, MIN(10, asm_opt.thread_num), sg->n_seq);
goto debug_gfa;
}
///just for debug
renew_graph_init(sources, reverse_sources, sg, coverage_cut, ruIndex, n_read);
@@ -29571,7 +29566,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
///normalize_ma_hit_t_single_side(sources, n_read);
normalize_ma_hit_t_single_side_advance(sources, n_read);
normalize_ma_hit_t_single_side_advance(reverse_sources, n_read);
if (ha_opt_triobin(&asm_opt))
{
drop_edges_by_trio(sources, n_read);
@@ -29580,7 +29574,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
{
memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t));
}
///print_binned_reads(sources, n_read, coverage_cut);
clean_weak_ma_hit_t(sources, reverse_sources, n_read);
@@ -29592,7 +29585,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
ma_hit_cut(sources, n_read, readLen, mini_overlap_length, &coverage_cut);
///print_binned_reads(sources, n_read, coverage_cut);
ma_hit_flt(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length);
///fix_binned_reads(sources, n_read, coverage_cut);
///just need to deal with trio here
ma_hit_contained_advance(sources, n_read, coverage_cut, ruIndex, max_hang_length, mini_overlap_length);
@@ -29601,7 +29593,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
///debug_info_of_specfic_node((char*)"m64043_200504_050026/93784180/ccs", sg, ruIndex, (char*)"sbsbsb");
init_bub_label_t(&b_mask_t, MIN(10, asm_opt.thread_num), sg->n_seq);
asg_arc_del_trans(sg, gap_fuzz);
asm_opt.coverage = get_coverage(sources, coverage_cut, n_read);
@@ -29617,7 +29608,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
asg_cut_tip(sg, asm_opt.max_short_tip);
///debug_info_of_specfic_node("m64043_200505_112554/8849050/ccs", sg, "inner_1");
///drop_inexact_edegs_at_bubbles(sg, bubble_dist);
if(clean_round > 0)
{
@@ -29696,7 +29686,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
asg_cut_tip(sg, asm_opt.max_short_tip);
}
}
if(VERBOSE >= 1)
{
fprintf(stderr, "\n\n**********final clean**********\n");
@@ -29715,7 +29704,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
asg_arc_del_orthology_multiple_way(sg, reverse_sources, 0.4, asm_opt.max_short_tip, ruIndex);
asg_cut_tip(sg, asm_opt.max_short_tip);
@@ -29725,7 +29713,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
asg_cut_tip(sg, asm_opt.max_short_tip);
asg_arc_del_simple_circle_untig(sources, coverage_cut, sg, 100, 0);
///note: don't apply asg_arc_del_too_short_overlaps() after this function!!!!
rescue_contained_reads_aggressive(NULL, sg, sources, coverage_cut, ruIndex, max_hang_length,
mini_overlap_length, 10, 1, 0, NULL, NULL, &b_mask_t);
@@ -29738,13 +29725,13 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
// rescue_no_coverage_aggressive(sg, sources, reverse_sources, &coverage_cut, ruIndex, max_hang_length,
// mini_overlap_length, bubble_dist, 10);
set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex,
max_hang_length, mini_overlap_length);
ruIndex->is_het = rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3,
ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz, &b_mask_t);
output_unitig_graph(sg, coverage_cut, o_file, sources, ruIndex, max_hang_length, mini_overlap_length);
flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL;
+1 -1
View File
@@ -1160,7 +1160,7 @@ 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 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, hap_cov_t **i_cov, bub_label_t* b_mask_t, uint32_t collect_p_trans);
kvec_asg_arc_t_warp* new_rtg_edges, hap_cov_t **i_cov, bub_label_t* b_mask_t, uint32_t collect_p_trans, uint32_t collect_p_trans_f);
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,
+7 -3
View File
@@ -5156,7 +5156,7 @@ void collect_purge_trans_cov(ma_ug_t *ug, hap_overlaps_list* ha, hap_cov_t *cov,
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, float drop_ratio, uint32_t just_contain,
uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t collect_p_trans_f)
{
p_g_t *pg = NULL;
asg_t* nsg = ug->g;
@@ -5235,16 +5235,20 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans)
normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex);
if(collect_p_trans && collect_p_trans_f == 0)
{
collect_purge_trans_cov(ug, &all_ovlp, cov, position_index);
}
if(asm_opt.polyploidy <= 2)
{
mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL);
///pt_solve(&all_ovlp, cov->t_ch, ug, read_g, 0.8, R_INF.trio_flag);
}
if(collect_p_trans)
if(collect_p_trans && collect_p_trans_f == 1)
{
collect_purge_trans_cov(ug, &all_ovlp, cov, position_index);
///goto end_coverage;
}
pg = init_p_g_t(ug, cov, read_g);
+1 -1
View File
@@ -66,7 +66,7 @@ typedef struct {
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, float drop_ratio, uint32_t just_contain,
uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans);
uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t collect_p_trans_f);
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);
+1 -1
View File
@@ -2267,7 +2267,7 @@ void dfs_bubble(asg_t *g, kvec_t_u32_warp* stack, kvec_t_u32_warp* result, uint3
uint32_t get_unitig_het_arb(ma_ug_t* ug, uint32_t uid, uint8_t *r_het_flag, kv_u_trans_t *ref, uint32_t m_het_occ,
uint32_t m_het_label, uint32_t p_het_label, uint32_t n_het_label)
{
if(u_trans_n(*ref, uid) > 0) return m_het_label;
if(ref && u_trans_n(*ref, uid) > 0) return m_het_label;
ma_utg_t *u = &(ug->u.a[uid]);
uint32_t k, rId;
uint32_t het_occ, hom_occ;
+129 -37
View File
@@ -28,6 +28,8 @@ uint8_t bit_filed[8] = {1, 2, 4, 8, 16, 32, 64, 128};
typedef struct {
int32_t max_iter;
int32_t n_perturb;
int32_t n_b_perturb;
int32_t n_s_perturb;
double f_perturb;
uint64_t seed;
} mc_opt_t;
@@ -111,13 +113,12 @@ typedef struct {
void mc_opt_init(mc_opt_t *opt, int32_t n_perturb, double f_perturb, uint64_t seed)
{
memset(opt, 0, sizeof(mc_opt_t));
// opt->n_perturb = 50000;
opt->n_perturb = n_perturb;
// opt->f_perturb = 0.1;
opt->f_perturb = f_perturb;
opt->max_iter = 1000;
// opt->seed = 11;
opt->seed = seed;
opt->n_s_perturb = n_perturb;
opt->n_b_perturb = n_perturb*0.5;
}
void mc_merge_dup(mc_g_t *mg) // MUST BE sorted
@@ -1774,6 +1775,7 @@ static void mc_perturb_node(const mc_opt_t *opt, const mc_match_t *ma, mc_svaux_
{
uint32_t i, k, n_bfs = 0;
k = (uint32_t)(kr_drand_r(&b->x) * b->cc_size + .499);
if(k >= b->cc_size) k = b->cc_size - 1;
k = (uint32_t)ma->cc[b->cc_off + k];///node id
n_bfs = mc_bfs(ma, b, k, bfs_round, (int32_t)(b->cc_size * opt->f_perturb));
for (i = 0; i < n_bfs; ++i)
@@ -1808,6 +1810,7 @@ static void mb_perturb_node(const mc_opt_t *opt, mb_g_t *mbg, mb_svaux_t *b, int
{
uint32_t i, k, n_bfs = 0;
k = (uint32_t)(kr_drand_r(&b->x) * b->cc_size + .499);
if(k >= b->cc_size) k = b->cc_size - 1;
k = (uint32_t)mbg->e->cc[b->cc_off + k];///node id
n_bfs = mb_bfs(mbg->e, b, k, bfs_round, (int32_t)(b->cc_size * opt->f_perturb));
for (i = 0; i < n_bfs; ++i)
@@ -2052,11 +2055,51 @@ void mc_solve_bp(mc_bp_t *bp)
fprintf(stderr, "[M::%s::%.3f] ==> round %u\n", __func__, yak_realtime()-index_time, r);
}
void print_sc(const mc_opt_t *opt, const mc_match_t *ma, mc_svaux_t *b, t_w_t sc_opt, uint32_t n_iter)
t_w_t mc_score_all_advance(const mc_match_t *ma, int8_t *s)
{
t_w_t w = mc_score(ma, b);
if(w != sc_opt) fprintf(stderr, "ERROR\n");
fprintf(stderr, "# iter: %u, sc_opt: %f, sc-local: %f, sc-global: %f\n", n_iter, sc_opt, w, mc_score_all(ma, b));
uint32_t k;
t_w_t z[2], zt = 0;
for (k = 0; k < ma->n_seq; ++k)
{
uint32_t o = ma->idx.a[k] >> 32;
uint32_t j, n = (uint32_t)ma->idx.a[k];
z[0] = z[1] = 0;
for (j = 0; j < n; ++j) {
const mc_edge_t *e = &ma->ma.a[o + j];
uint32_t t = ma_y(*e);
if (s[t] > 0) z[0] += e->w;
else if (s[t] < 0) z[1] += e->w;
}
zt += -((t_w_t)(s[k])) * (z[0] - z[1]);
}
return zt;
}
t_w_t mb_score_all_advance(const mc_match_t *ma, mb_g_t *mbg)
{
uint32_t k;
t_w_t z[2], zt = 0;
for (k = 0; k < ma->n_seq; ++k)
{
uint32_t o = ma->idx.a[k] >> 32;
uint32_t j, n = (uint32_t)ma->idx.a[k];
z[0] = z[1] = 0;
for (j = 0; j < n; ++j) {
const mc_edge_t *e = &ma->ma.a[o + j];
uint32_t t = ma_y(*e);
if (mbg->u->u.a[mbg->u->idx.a[t]>>1].s[mbg->u->idx.a[t]&1] > 0) z[0] += e->w;
else if (mbg->u->u.a[mbg->u->idx.a[t]>>1].s[mbg->u->idx.a[t]&1] < 0) z[1] += e->w;
}
zt += -((t_w_t)(mbg->u->u.a[mbg->u->idx.a[k]>>1].s[mbg->u->idx.a[k]&1])) * (z[0] - z[1]);
}
return zt;
}
void print_sc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, t_w_t sc_opt, uint32_t n_iter)
{
t_w_t w = mc_score(mg->e, b);
// if(w != sc_opt) fprintf(stderr, "ERROR\n");
fprintf(stderr, "# iter: %u, sc_opt: %f, sc-local: %f, sc-global: %f\n", n_iter, sc_opt, w, mc_score_all_advance(mg->e, mg->s.a));
}
uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint32_t cc_off, uint32_t cc_size)
@@ -2065,8 +2108,8 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3
t_w_t sc_opt = -(1<<30), sc;///problem-w
b->cc_off = cc_off, b->cc_size = cc_size;
if (b->cc_size < 2) return 0;
// print_sc(opt, mg->e, b, sc_opt, (uint32_t)-1);
sc_opt = mc_init_spin(mg->e, b);
// print_sc(opt, mg, b, sc_opt, n_iter);
if (b->cc_size == 2) return 0;
for (j = 0; j < b->cc_size; ++j) {///backup s and z in s_opt and z_opt
b->s_opt[b->cc_node[j]] = b->s[b->cc_node[j]]; ///hap status of each unitig
@@ -2086,6 +2129,7 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3
b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]];
}
}
// print_sc(opt, mg, b, sc_opt, n_iter);
// mc_reset_z_debug(mg->e, b);
// print_sc(opt, mg->e, b, sc_opt, n_iter);
// fprintf(stderr, "\ncc_size: %u, cc_off: %u\n", b->cc_size, b->cc_off);
@@ -2100,6 +2144,7 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3
b->z_opt[b->cc_node[j]] = b->z[b->cc_node[j]];
}
sc_opt = sc;
// print_sc(opt, mg, b, sc_opt, n_iter);
} else {
for (j = 0; j < b->cc_size; ++j) {
b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]];
@@ -2117,11 +2162,14 @@ uint32_t mc_solve_cc(const mc_opt_t *opt, const mc_g_t *mg, mc_svaux_t *b, uint3
}
sc_opt = sc;
}
// print_sc(opt, mg->e, b, sc_opt, n_iter);
}
for (j = 0; j < b->cc_size; ++j)
{
b->s[b->cc_node[j]] = b->s_opt[b->cc_node[j]];
b->z[b->cc_node[j]] = b->z_opt[b->cc_node[j]];
}
return n_iter;
}
@@ -2291,6 +2339,7 @@ void mc_init_spin_all(const mc_opt_t *opt, mc_g_t *mg, mb_g_t *mbg, mc_svaux_t *
}
if(!mbg) return;
// memcpy(b->s_opt, b->s, sizeof(int8_t)*mg->e->n_seq);
///adjust by block
kvec_t(uint32_t) s; kv_init(s);
@@ -2385,6 +2434,7 @@ void mc_init_spin_all(const mc_opt_t *opt, mc_g_t *mg, mb_g_t *mbg, mc_svaux_t *
/*******************************for debug************************************/
// debug_mbg(mbg, mg, b);
/*******************************for debug************************************/
// memcpy(b->s, b->s_opt, sizeof(int8_t)*mg->e->n_seq);
}
@@ -2448,6 +2498,9 @@ void mb_g_cc(mb_g_t *mbg)
void mc_set_by_mbg(mc_g_t *mg, mb_g_t *mbg)
{
// t_w_t z0 = mc_score_all_advance(mg->e, mg->s.a);
// t_w_t z1 = mb_score_all_advance(mg->e, mbg);
// if(z0 >= z1) return;
uint32_t i, k, qn, *a[2], a_n[2];
int8_t s[2];
for (i = 0; i < mbg->u->u.n; i++)
@@ -2502,27 +2555,9 @@ void print_mb_g_blcok(mb_g_t *mbg)
}
}
t_w_t mc_score_all_advance(const mc_match_t *ma, int8_t *s)
{
uint32_t k;
t_w_t z[2], zt = 0;
for (k = 0; k < ma->n_seq; ++k)
{
uint32_t o = ma->idx.a[k] >> 32;
uint32_t j, n = (uint32_t)ma->idx.a[k];
z[0] = z[1] = 0;
for (j = 0; j < n; ++j) {
const mc_edge_t *e = &ma->ma.a[o + j];
uint32_t t = ma_y(*e);
if (s[t] > 0) z[0] += e->w;
else if (s[t] < 0) z[1] += e->w;
}
zt += -((t_w_t)(s[k])) * (z[0] - z[1]);
}
return zt;
}
void mb_solve_core(const mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t is_sys)
void mb_solve_core(mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t is_sys)
{
if(!ref) return;
double index_time = yak_realtime();
@@ -2530,11 +2565,11 @@ void mb_solve_core(const mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t
mb_g_t *mbg = init_mb_g_t(mg, ref, is_sys);
mb_svaux_t *bb;
/**************************init**************************/
fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
mc_svaux_t *b;
mc_g_cc(mg->e);
b = mc_svaux_init(mg, opt->seed);
mc_init_spin_all(opt, mg, mbg, b);
mc_svaux_destroy(b);
free(mg->e->cc);
mg->e->cc = NULL;
@@ -2543,18 +2578,21 @@ void mb_solve_core(const mc_opt_t *opt, mc_g_t *mg, kv_u_trans_t *ref, uint32_t
bb = mb_svaux_init(mbg, opt->seed);
/*******************************for debug************************************/
// print_mb_g_blcok(mbg);
fprintf(stderr, "*********before-[M::%s::mc_score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
fprintf(stderr, "*********before-[M::%s::mb_score->%f] ==> Partition\n", __func__, mb_score_all_advance(mg->e, mbg));
/*******************************for debug************************************/
fprintf(stderr, "\n\n\n\n\n*************beg-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
opt->n_perturb = opt->n_b_perturb;
for (st = 0, i = 1; i <= mbg->e->n_seq; ++i) {
if (i == mbg->e->n_seq || mbg->e->cc[st]>>32 != mbg->e->cc[i]>>32) {
mb_solve_cc(opt, mbg, bb, st, i - st);
st = i;
}
}
opt->n_perturb = opt->n_s_perturb - opt->n_b_perturb;
/*******************************for debug************************************/
// debug_mb_solve_core(mbg);
fprintf(stderr, "*********after-[M::%s::mc_score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
fprintf(stderr, "*********after-[M::%s::mb_score->%f] ==> Partition\n", __func__, mb_score_all_advance(mg->e, mbg));
/*******************************for debug************************************/
mc_set_by_mbg(mg, mbg);
fprintf(stderr, "##############end-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
@@ -2583,7 +2621,7 @@ void mc_solve_core(const mc_opt_t *opt, mc_g_t *mg, bubble_type* bub)
st = i;
}
}
fprintf(stderr, "##############end-[M::%s::score->%f] ==> Partition\n", __func__, mc_score_all_advance(mg->e, mg->s.a));
fprintf(stderr, "##############end-[---M::%s::score->%f] ==> Partition\n", __func__, mc_score_all(mg->e, b));
if(bp) mc_solve_bp(bp);
///mc_write_info(g, b);
mc_svaux_destroy(b);
@@ -2694,6 +2732,60 @@ void p_nodes(mc_g_t *mg, trans_chain* t_ch, uint8_t* trio_flag)
}
}
void write_mc_g_t(mc_opt_t *opt, mc_g_t *mg, const char *name)
{
FILE* fp = fopen(name, "w");
fwrite(opt, sizeof(mc_opt_t), 1, fp);
fwrite(&(mg->s.n), sizeof(mg->s.n), 1, fp);
fwrite(mg->s.a, sizeof(mc_node_t), mg->s.n, fp);
fwrite(&(mg->e->n_seq), sizeof(mg->e->n_seq), 1, fp);
fwrite(&(mg->e->ma.n), sizeof(mg->e->ma.n), 1, fp);
fwrite(mg->e->ma.a, sizeof(mc_edge_t), mg->e->ma.n, fp);
fwrite(&(mg->e->idx.n), sizeof(mg->e->idx.n), 1, fp);
fwrite(mg->e->idx.a, sizeof(uint64_t), mg->e->idx.n, fp);
fclose(fp);
}
mc_g_t* load_mc_g_t(mc_opt_t *opt, const char *name)
{
FILE* fp = NULL;
fp = fopen(name, "r");
if(!fp) return NULL;
uint64_t flag = 0;
mc_g_t *mg = NULL; CALLOC(mg, 1);
kv_init(mg->s); CALLOC(mg->e, 1);
flag += fread(opt, sizeof(mc_opt_t), 1, fp);
flag += fread(&(mg->s.n), sizeof(mg->s.n), 1, fp);
mg->s.m = mg->s.n; MALLOC(mg->s.a, mg->s.n);
flag += fread(mg->s.a, sizeof(mc_node_t), mg->s.n, fp);
flag += fread(&(mg->e->n_seq), sizeof(mg->e->n_seq), 1, fp);
flag += fread(&(mg->e->ma.n), sizeof(mg->e->ma.n), 1, fp);
mg->e->ma.m = mg->e->ma.n; MALLOC(mg->e->ma.a, mg->e->ma.n);
flag += fread(mg->e->ma.a, sizeof(mc_edge_t), mg->e->ma.n, fp);
flag += fread(&(mg->e->idx.n), sizeof(mg->e->idx.n), 1, fp);
mg->e->idx.m = mg->e->idx.n; MALLOC(mg->e->idx.a, mg->e->idx.n);
flag += fread(mg->e->idx.a, sizeof(uint64_t), mg->e->idx.n, fp);
fclose(fp);
return mg;
}
void debug_mc_g_t(const char* name)
{
mc_opt_t opt;
mc_g_t *mg = load_mc_g_t(&opt, name);
mc_solve_core(&opt, mg, NULL);
destory_mc_g_t(&mg);
exit(1);
}
void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref)
{
@@ -2704,6 +2796,7 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u
mb_solve_core(&opt, mg, ref, is_sys);
///debug_mc_g_t(mg);
if(renew_s == 0) write_mc_g_t(&opt, mg, MC_NAME);
mc_solve_core(&opt, mg, bub);
if((asm_opt.flag & HA_F_PARTITION) && t_ch)
@@ -2713,6 +2806,5 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u
if(ovlp) clean_ovlp_by_mc(mg, ovlp);
destory_mc_g_t(&mg);
destory_mc_g_t(&mg);
}
+2 -1
View File
@@ -22,7 +22,7 @@ typedef struct {
#define w_t double
#define t_w_t double
#define w_cast(x) ((t_w_t)((x)))
#define MC_NAME "debug_mc.bin"
typedef struct {
uint64_t x; ///(uint64_t)nid1 << 32 | nid2;
@@ -89,4 +89,5 @@ static inline double kr_drand_r(uint64_t *x)
}
void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref);
void debug_mc_g_t(const char* name);
#endif