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

Merge pull request #421 from chhylp123/hifiasm_dev_debug

Browse Source 
Hifiasm dev debug
This commit is contained in:
chhylp123
2023-03-13 11:14:38 -04:00
committed by GitHub
parent f69166ee3b f67faa9c34
commit 7bb366688e
8 changed files with 367 additions and 28 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.cpp
+13
View File
@@ -58,6 +58,9 @@ static ko_longopt_t long_options[] = {
{ "prt-raw", ko_no_argument, 343},
{ "integer-correct", ko_required_argument, 344},
{ "dbg-ovec", ko_no_argument, 345},
{ "path-max", ko_required_argument, 346},
{ "path-min", ko_required_argument, 347},
// { "path-round", ko_required_argument, 348},
{ 0, 0, 0 }
};
@@ -160,6 +163,12 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " error rate of Ultra-Long reads [%.3g]\n", asm_opt->ul_error_rate);
fprintf(stderr, " --ul-tip INT\n");
fprintf(stderr, " remove tip unitigs composed of <=INT reads for the UL assembly [%d]\n", asm_opt->max_short_ul_tip);
fprintf(stderr, " --path-max FLOAT\n");
fprintf(stderr, " max path drop ratio [%.2g]; higher number may make the assembly cleaner\n", asm_opt->max_path_drop_rate);
fprintf(stderr, " but may lead to more misassemblies\n");
fprintf(stderr, " --path-min FLOAT\n");
fprintf(stderr, " min path drop ratio [%.2g]; higher number may make the assembly cleaner\n", asm_opt->min_path_drop_rate);
fprintf(stderr, " but may lead to more misassemblies\n");
// fprintf(stderr, " --low-het enable it for genomes with very low het heterozygosity rate (<0.0001%%)\n");
fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n");
@@ -271,6 +280,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->prt_dbg_gfa = 0;
asm_opt->integer_correct_round = 0;
asm_opt->dbg_ovec_cal = 0;
asm_opt->min_path_drop_rate = 0.2;
asm_opt->max_path_drop_rate = 0.6;
}
void destory_enzyme(enzyme* f)
@@ -807,6 +818,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 343) asm_opt->prt_dbg_gfa = 1;
else if (c == 344) asm_opt->integer_correct_round = atol(opt.arg);
else if (c == 345) asm_opt->dbg_ovec_cal = 1;
else if (c == 346) asm_opt->max_path_drop_rate = atof(opt.arg);
else if (c == 347) asm_opt->min_path_drop_rate = atof(opt.arg);
else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap
asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);
}
CommandLines.h
+4 -1
View File
@@ -5,7 +5,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.19.0-r555"
#define HA_VERSION "0.19.1-r559"
#define VERBOSE 0
@@ -104,6 +104,9 @@ typedef struct {
float purge_simi_thres;
float trans_base_rate;
float trans_base_rate_sec;
float min_path_drop_rate;
float max_path_drop_rate;
// uint64_t path_clean_round;
///float purge_simi_rate_hic;
Overlaps.cpp
+12 -3
View File
@@ -13528,7 +13528,7 @@ void hic_clean_adv(asg_t *sg, ug_opt_t *uopt)
{
uint32_t i, k, m, z, v, w, mk; ma_utg_t *u = NULL; uint32_t *ba, bn, n_vtx, beg, end, n0, n1; ma_utg_t *mz = NULL;
ma_ug_t *ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); n_vtx = ug->g->n_seq<<1; double bub_rate = 0.1;
uint8_t *bf = NULL; bubble_type *bub = gen_bubble_chain(sg, ug, uopt, &bf);
uint8_t *bf = NULL; bubble_type *bub = gen_bubble_chain(sg, ug, uopt, &bf, ((asm_opt.polyploidy>2)?1:0));
uint64_t tLen, vocc, socc, pocc; buf_t b; memset(&b, 0, sizeof(buf_t)); CALLOC(b.a, n_vtx);
REALLOC(bf, n_vtx); memset(bf, 0, sizeof((*bf))*n_vtx);
kvec_t(uint64_t) buf; kv_init(buf); n0 = n1 = 0;
@@ -15605,7 +15605,7 @@ void clean_u_trans_t_idx(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g)
void gen_bp_phasing(ug_opt_t *opt, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *sg)
{
uint8_t *bf = NULL;
bubble_type *bub = gen_bubble_chain(sg, ug, opt, &bf); free(bf);
bubble_type *bub = gen_bubble_chain(sg, ug, opt, &bf, 0); free(bf);
filter_u_trans(ta, asm_opt.is_bub_trans, asm_opt.is_topo_trans, asm_opt.is_read_trans, asm_opt.is_base_trans);
// dbg_prt_utg_trans(ta, ug, "pre");
if(asm_opt.is_base_trans) {
@@ -33734,7 +33734,7 @@ bub_label_t *b_mask_t, uint32_t is_trio, int32_t ul_aln_round, char *o_file, con
int32_t k, strl = strlen(bin_file)+1, kt, cl, sl; char *id = NULL;
renew_g(src, rev_src, n_read, readLen, cov, ruIndex, sg, mini_overlap_length, max_hang_length,
uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, asm_opt.max_short_tip, b_mask_t,
is_trio, o_file, bin_file, (ul_aln_round<=1)?1:0, 1, /**((is_trio)?(0):(1))**/1);
is_trio, o_file, bin_file, (ul_aln_round<=1)?1:0, 1, /**((is_trio)?(0):(1))**/((asm_opt.polyploidy<=2)?1:0));
gen_ug_opt_t(uopt, *src, *rev_src, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, *readLen,
*cov, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, b_mask_t);
ug_ext_gfa(uopt, *sg, ug_ext_len);
@@ -34019,6 +34019,15 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources,
asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length, &uopt);
// if(asm_opt.ar) {
// char *op_file = get_outfile_name(o_file);
// ul_clean_gfa(&uopt, sg, sources, reverse_sources, ruIndex, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
// 0.6, asm_opt.max_short_tip, gap_fuzz, &b_mask_t, 0/**!!asm_opt.ar**/, ha_opt_triobin(&asm_opt), UL_COV_THRES, op_file);
// output_contig_graph_primary_pre(sg, coverage_cut, op_file, sources, reverse_sources,
// asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length, &uopt);
// free(op_file);
// }
/**
if (asm_opt.flag & HA_F_VERBOSE_GFA)
{
gfa_ut.cpp
+19 -7
View File
@@ -2995,7 +2995,7 @@ void ug_ext_gfa(ug_opt_t *uopt, asg_t *sg, uint32_t max_len)
}
bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het)
bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het, uint8_t avoid_het)
{
kvec_asg_arc_t_warp new_rtg_edges;
kv_init(new_rtg_edges.a);
@@ -3013,7 +3013,14 @@ bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **
// fprintf(stderr, "1[M::%s]\n", __func__);
CALLOC(bub, 1); (*ir_het) = cov->t_ch->ir_het;
cov->t_ch->ir_het = NULL; cov->is_r_het = NULL;
identify_bubbles(ug, bub, (*ir_het), NULL);
if(!avoid_het) {
identify_bubbles(ug, bub, (*ir_het), NULL);
} else {
uint8_t *pr_het; CALLOC(pr_het, sg->n_seq);
identify_bubbles_recal_poy(sg, ug, bub, pr_het, uopt->sources, uopt->ruIndex, NULL);
free(pr_het);
}
// fprintf(stderr, "2[M::%s]\n", __func__);
kv_destroy(new_rtg_edges.a); destory_hap_cov_t(&cov);
if(asm_opt.purge_level_primary == 0) {///all nodes are het
@@ -13621,7 +13628,12 @@ void renew_usg_t_bub(ul_resolve_t *uidx, usg_t *ng, uint32_t *id_map, uint8_t *f
// fprintf(stderr, "[M::%s] homozygous read coverage threshold: %d\n", __func__, asm_opt.hom_global_coverage_set?
// asm_opt.hom_global_coverage:(int)(((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE)));
// identify_bubbles(ug, uidx->bub, uidx->r_het, NULL);
identify_bubbles_recal(uidx->sg, ug, uidx->bub, uidx->r_het, uidx->uopt->sources, uidx->uopt->ruIndex, NULL);
if(asm_opt.polyploidy <= 2) {
identify_bubbles_recal(uidx->sg, ug, uidx->bub, uidx->r_het, uidx->uopt->sources, uidx->uopt->ruIndex, NULL);
} else {
identify_bubbles_recal_poy(uidx->sg, ug, uidx->bub, uidx->r_het, uidx->uopt->sources, uidx->uopt->ruIndex, NULL);
}
// fprintf(stderr, "0[M::%s::] f[51]::%u\n", __func__, ff[51]);
for (i = 0; i < ng->n; i++) {
if(id_map[i] != (uint32_t)-1) {
@@ -15922,11 +15934,11 @@ void u2g_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint32_t keep_raw_utg, uint
// char sb[1000];
// output_integer_graph(uidx, iug, "ig_h0", 0);
// fprintf(stderr, "\n[M::%s::] max_path_drop_ratio::%f, min_path_drop_ratio::%f\n", __func__, ulopt->max_path_drop_ratio, ulopt->min_path_drop_ratio);
for (ss = 0; ss < 2; ss++) {
for (i = 0, drop = ulopt->min_path_drop_ratio; i < ulopt->clean_round; i++, drop += step) {
if(drop > ulopt->max_path_drop_ratio) drop = ulopt->max_path_drop_ratio;
// fprintf(stderr, "\n[M::%s::] Starting round-%ld, drop::%f\n", __func__, i, drop);
// fprintf(stderr, "[M::%s::] Starting round-%ld, drop::%f\n", __func__, i, drop);
cnt = 1; topo_level = 2; mm_tip = ulopt->max_tip;
while (cnt) {
cnt = 0;
@@ -17245,7 +17257,7 @@ ul_renew_t *ropt, const char *bin_file, uint64_t free_uld, uint64_t is_bridg, ui
// fprintf(stderr, "-0-[M::%s]\tUL_INF.a[25].rlen::%u\n", __func__, UL_INF.a[25].rlen);
// print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0);
// print_ul_alignment(init_ug, &UL_INF, 47072, "after-0");
bub = gen_bubble_chain(sg, init_ug, uopt, &r_het);
bub = gen_bubble_chain(sg, init_ug, uopt, &r_het, ((asm_opt.polyploidy>2)?1:0));
// fprintf(stderr, "4[M::%s]\n", __func__);
// print_ul_alignment(init_ug, &UL_INF, 47072, "after-1");
ul_resolve_t *uidx = init_ul_resolve_t(sg, init_ug, bub, &UL_INF, uopt, r_het);
@@ -17262,7 +17274,7 @@ ul_renew_t *ropt, const char *bin_file, uint64_t free_uld, uint64_t is_bridg, ui
// }
ul_re_correct(uidx, asm_opt.integer_correct_round/**3**/);
init_ulg_opt_t(&uu, uopt, clean_round, 0.2, 0.6, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_ul_tip, b_mask_t, is_trio);
init_ulg_opt_t(&uu, uopt, clean_round, asm_opt.min_path_drop_rate, asm_opt.max_path_drop_rate, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_ul_tip, b_mask_t, is_trio);
// print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug0", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1);
/**ul2ul_idx_t *u2o = **/gen_ul2ul(uidx, uopt, &uu, 0, is_bridg);
// print_ul_alignment(init_ug, &UL_INF, 47072, "after-3");
gfa_ut.h
+1 -1
View File
@@ -35,7 +35,7 @@ void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd);
asg_t *gen_ng(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **cov, R_to_U *ruI, uint64_t scaffold_len);
void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover);
// void print_raw_u2rgfa_seq(all_ul_t *aln, R_to_U* rI, uint32_t is_detail);
bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het);
bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het, uint8_t avoid_het);
void filter_sg_by_ug(asg_t *rg, ma_ug_t *ug, ug_opt_t *uopt);
void ug_ext_gfa(ug_opt_t *uopt, asg_t *sg, uint32_t max_len);
void update_sg_uo(asg_t *g, ma_hit_t_alloc *src);
hic.cpp
+315 -15
View File
@@ -2520,7 +2520,8 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t
// fprintf(stderr, "-bub->index[18759]: %u, bub->num.n: %u\n", (uint32_t)bub->index[18759], bub->num.n);
}
uint32_t get_unitig_het_fly(ma_ug_t* ug, uint32_t uid, asg_t* sg, int64_t het_cov_thres,
uint32_t get_unitig_het_fly(ma_ug_t* ug, uint32_t uid, asg_t* sg, /**int64_t het_cov_thres,**/
int64_t m_het_cov_thres, int64_t m_hom_cov_thres,
ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag, uint32_t m_het_occ, uint32_t m_het_label,
uint32_t p_het_label, uint32_t n_het_label)
{
@@ -2596,8 +2597,10 @@ uint32_t p_het_label, uint32_t n_het_label)
// fprintf(stderr, "[M::%s::uid->%u] u->n::%u, C_bases::%ld, R_bases::%ld, het_cov_thres::%ld, m_het_occ::%u\n",
// __func__, uid, (uint32_t)u->n, C_bases, R_bases, het_cov_thres, m_het_occ);
// }
if((cov <= (het_cov_thres*1.333333)) && (u->n >= m_het_occ)) return m_het_label; ///must het
if((cov >= (het_cov_thres*1.6))) return n_het_label; ///hom
// if((cov <= (het_cov_thres*1.333333)) && (u->n >= m_het_occ)) return m_het_label; ///must het
// if((cov >= (het_cov_thres*1.6))) return n_het_label; ///hom
if((cov <= m_het_cov_thres) && (u->n >= m_het_occ)) return m_het_label; ///must het
if(cov >= m_hom_cov_thres) return n_het_label; ///hom
return p_het_label; ///potential het
}
@@ -2608,9 +2611,20 @@ kv_u_trans_t *ref)
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;
uint64_t pathLen, hom_cov, het_cov, m_het_cov, m_hom_cov;
bub->ug = ug;
bub->b_bub = bub->b_end_bub = bub->tangle_bub = bub->cross_bub = bub->mess_bub = 0;
if(asm_opt.hom_global_coverage_set) {
hom_cov = asm_opt.hom_global_coverage;
} else {
hom_cov = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE);
}
het_cov = hom_cov/asm_opt.polyploidy;
m_het_cov = hom_cov - het_cov + (het_cov*0.333333);
m_hom_cov = hom_cov - het_cov + (het_cov*0.6);
// fprintf(stderr, "hom_cov::%lu, het_cov::%lu, m_het_cov::%lu, m_hom_cov::%lu\n", hom_cov, het_cov, m_het_cov, m_hom_cov);
if(bub->round_id == 0)
{
buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t));
@@ -2697,20 +2711,11 @@ kv_u_trans_t *ref)
kv_push(uint32_t, bub->num, bub->list.n);
free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
bub->f_bub = bub->num.n - 1; ///bub->s_bub = bub->num.n - 1;
uint64_t dip_thre_max;
memset(r_het_flag, 0, sizeof((*r_het_flag))*sg->n_seq);
if(asm_opt.hom_global_coverage_set) {
dip_thre_max = asm_opt.hom_global_coverage;
} else {
dip_thre_max = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE);
}
// dip_thre_max = (double)(dip_thre_max) - (((double)(dip_thre_max)*0.5)/asm_opt.polyploidy);
dip_thre_max = (double)(dip_thre_max) - ((double)(dip_thre_max)/asm_opt.polyploidy);
for (i = 0; i < ug->g->n_seq; i++) {
bub->index[i] = get_unitig_het_fly(ug, i, sg, dip_thre_max, sources, ruIndex,
r_het_flag, 20, M_het(*bub), P_het(*bub), (uint32_t)-1);
bub->index[i] = get_unitig_het_fly(ug, i, sg, m_het_cov, m_hom_cov,
sources, ruIndex, r_het_flag, 20, M_het(*bub), P_het(*bub), (uint32_t)-1);
}
for (i = 0; i < bub->f_bub; i++)
@@ -2782,6 +2787,301 @@ kv_u_trans_t *ref)
// fprintf(stderr, "-bub->index[18759]: %u, bub->num.n: %u\n", (uint32_t)bub->index[18759], bub->num.n);
}
void reset_inner_bub_het_poy(asg_t* sg, ma_ug_t* ug, bubble_type* bub, uint32_t bid, uint64_t tLen, buf_t *b,
uint64_t m_het_cov, uint64_t m_hom_cov, uint8_t *r_het_flag, ma_hit_t_alloc* sources, R_to_U* ruIndex)
{
uint32_t beg, sink, *ba, bn, m, v, z, socc;
get_bubbles(bub, bid, &beg, &sink, &ba, &bn, NULL);
for (m = 0; m < bn; m++) {
v = ba[m];
if(IF_HOM((v>>1), *bub)) continue;
if(ug->g->seq[v>>1].del) continue;
if(bub->index[v>>1] == (uint32_t)-1) continue;
if(get_unitig_het_fly(ug, v>>1, sg, m_het_cov, m_hom_cov, sources, ruIndex, r_het_flag, 20, 1, 0, (uint32_t)-1) == 1) {
continue;
}
if(asg_arc_n(ug->g, v) < 2) continue;
if(get_real_length(ug->g, v, NULL) < 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, NULL)) {
//beg is v, end is b.S.a[0]
//note b.b include end, does not include beg
for (z = socc = 0; z < b->b.n; z++) {
if(b->b.a[z]==v || b->b.a[z]==b->S.a[0]) continue;
socc += ug->u.a[b->b.a[z]>>1].n;
}
if((socc > 10) && (socc > (ug->u.a[v>>1].n*3))) {
bub->index[v>>1] = (uint32_t)-1;
}
}
v ^= 1;
if(bub->index[v>>1] == (uint32_t)-1) continue;
if(asg_arc_n(ug->g, v) < 2) continue;
if(get_real_length(ug->g, v, NULL) < 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, NULL)) {
//beg is v, end is b.S.a[0]
//note b.b include end, does not include beg
for (z = socc = 0; z < b->b.n; z++) {
if(b->b.a[z]==v || b->b.a[z]==b->S.a[0]) continue;
socc += ug->u.a[b->b.a[z]>>1].n;
}
if((socc > 10) && (socc > (ug->u.a[v>>1].n*3))) {
bub->index[v>>1] = (uint32_t)-1;
}
}
}
}
uint32_t deter_unitig_het_fly(ma_ug_t* ug, uint32_t uid, asg_t* sg, int64_t het_cov_thres,
ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag)
{
ma_utg_t *u = &(ug->u.a[uid]);
uint32_t k, i, j, rId, nv, tn, is_Unitig;
asg_arc_t *av = NULL; ma_hit_t *h;
int64_t R_bases = 0, C_bases = 0, cov;
///set
u = &(ug->u.a[uid]);
for (k = 0; k < u->n; k++) {
rId = u->a[k]>>33;
r_flag[rId] = 1;
}
for (i = 0; i < 2; i++) {
nv = asg_arc_n(ug->g, (uid<<1)+i);
av = asg_arc_a(ug->g, (uid<<1)+i);
for (j = 0; j < nv; j++) {
u = &(ug->u.a[av[j].v>>1]);
for (k = 0; k < u->n; k++) {
rId = u->a[k]>>33;
r_flag[rId] = 2;
}
}
}
u = &(ug->u.a[uid]);
for (k = 0; k < u->n; k++) {
if(u->a[k] == (uint64_t)-1) continue;
rId = u->a[k]>>33;
R_bases += sg->seq[rId].len;
for (j = 0; j < (uint64_t)(sources[rId].length); j++) {
h = &(sources[rId].buffer[j]);
if(h->el != 1) continue;
tn = Get_tn((*h));
if(sg->seq[tn].del == 1) {
///get the id of read that contains it
get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
if(tn == (uint32_t)-1 || is_Unitig == 1 || sg->seq[tn].del == 1) continue;
}
if(sg->seq[tn].del == 1) continue;
if(r_flag[tn] == 0) continue;
if(r_flag[tn] == 1) {
C_bases += (Get_qe((*h)) - Get_qs((*h)));
}
if(r_flag[tn] == 2) {
C_bases += ((Get_qe((*h)) - Get_qs((*h)))/2);
}
}
}
///reset
u = &(ug->u.a[uid]);
for (k = 0; k < u->n; k++) {
rId = u->a[k]>>33;
r_flag[rId] = 0;
}
for (i = 0; i < 2; i++) {
nv = asg_arc_n(ug->g, (uid<<1)+i);
av = asg_arc_a(ug->g, (uid<<1)+i);
for (j = 0; j < nv; j++) {
u = &(ug->u.a[av[j].v>>1]);
for (k = 0; k < u->n; k++) {
rId = u->a[k]>>33;
r_flag[rId] = 0;
}
}
}
u = &(ug->u.a[uid]); cov = 0;
if(R_bases > 0) cov = C_bases/R_bases;
if(cov <= het_cov_thres) return 1;
return 0;
}
void identify_bubbles_recal_poy(asg_t* sg, ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, ma_hit_t_alloc* sources, R_to_U* ruIndex,
kv_u_trans_t *ref)
{
asg_cleanup(ug->g);
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, hom_cov, het_cov, m_het_cov, m_hom_cov;
bub->ug = ug;
bub->b_bub = bub->b_end_bub = bub->tangle_bub = bub->cross_bub = bub->mess_bub = 0;
if(asm_opt.hom_global_coverage_set) {
hom_cov = asm_opt.hom_global_coverage;
} else {
hom_cov = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE);
}
het_cov = hom_cov/asm_opt.polyploidy;
m_het_cov = het_cov + (het_cov*0.333333);///hom_cov - het_cov + (het_cov*0.333333);
m_hom_cov = het_cov + (het_cov*0.6);///hom_cov - het_cov + (het_cov*0.6);
// fprintf(stderr, "hom_cov::%lu, het_cov::%lu, m_het_cov::%lu, m_hom_cov::%lu\n", hom_cov, het_cov, m_het_cov, m_hom_cov);
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);
bub->b_s_idx.n = ug->g->n_seq;
memset(bub->b_s_idx.a, -1, bub->b_s_idx.n * sizeof(uint64_t));
CALLOC(bub->index, n_vtx);
for (i = 0; i < ug->g->n_seq; i++) ug->g->seq[i].c = 0;
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((bub->index[v]&(uint32_t)3) != 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, NULL))
{
//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;
bub->index[b.b.a[i]] &= mode; bub->index[b.b.a[i]] += 1;
bub->index[b.b.a[i]^1] &= mode; bub->index[b.b.a[i]^1] += 1;
}
bub->index[v] &= mode; bub->index[v] += 2;
bub->index[b.S.a[0]^1] &= mode; bub->index[b.S.a[0]^1] += 3;
}
}
kvec_t_u32_warp stack, result;
kv_init(stack.a); kv_init(result.a);
for (v = 0; v < n_vtx; ++v)
{
if((bub->index[v]&(uint32_t)3) !=2) continue;
if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL, 0, 0, NULL))
{
//note b.b include end, does not include beg
i = b.b.n + 1;
if(b.b.n == 2 || b.b.n == 3 || b.b.n == 5)
{
for (i = 0; i < b.b.n; i++)
{
if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue;
dfs_bubble(ug->g, &stack, &result, b.b.a[i]>>1, v>>1, b.S.a[0]>>1);
if((result.a.n + 3) != b.b.n && (result.a.n + 2) != b.b.n) break;
}
}
if(i == b.b.n)
{
kv_push(uint32_t, bub->num, v);
}
else
{
kv_push(uint32_t, bub->num, v + (1<<31));
}
}
}
kv_destroy(stack.a); kv_destroy(result.a);
radix_sort_u32(bub->num.a, bub->num.a + bub->num.n);
bub->s_bub = 0;
for (k = 0; k < bub->num.n; k++)
{
if((bub->num.a[k]>>31) == 0) bub->s_bub++;
v = (bub->num.a[k]<<1)>>1;
bub->num.a[k] = bub->list.n;
if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL, 0, 0, NULL))
{
kv_push(uint64_t, bub->pathLen, pathLen);
//beg is v, end is b.S.a[0]
kv_push(uint32_t, bub->list, v);
kv_push(uint32_t, bub->list, b.S.a[0]^1);
//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;
kv_push(uint32_t, bub->list, b.b.a[i]);
}
}
}
kv_push(uint32_t, bub->num, bub->list.n);
// free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
bub->f_bub = bub->num.n - 1; ///bub->s_bub = bub->num.n - 1;
memset(r_het_flag, 0, sizeof((*r_het_flag))*sg->n_seq);
for (i = 0; i < ug->g->n_seq; i++) {
bub->index[i] = get_unitig_het_fly(ug, i, sg, m_het_cov, m_hom_cov,
sources, ruIndex, r_het_flag, 20, M_het(*bub), P_het(*bub), (uint32_t)-1);
}
for (i = 0; i < bub->f_bub; i++) {
get_bubbles(bub, i, &beg, &sink, &a, &n, &pathLen);
for (v = n_occ = 0; v < n; v++) {
bub->index[(a[v]>>1)] = i;
n_occ += ug->u.a[a[v]>>1].n;
}
if(n_occ > 3) {
if(bub->index[(beg>>1)] != M_het(*bub)) bub->index[(beg>>1)] = (uint32_t)-1;
if(bub->index[(sink>>1)] != M_het(*bub)) bub->index[(sink>>1)] = (uint32_t)-1;
}
v = beg>>1;
if(bub->b_s_idx.a[v] == (uint64_t)-1) {
bub->b_s_idx.a[v] <<= 32;
bub->b_s_idx.a[v] |= i;
} else if((bub->b_s_idx.a[v] & 0xffffffff00000000) == 0xffffffff00000000) {
bub->b_s_idx.a[v] <<= 32;
bub->b_s_idx.a[v] |= i;
}
v = sink>>1;
if(bub->b_s_idx.a[v] == (uint64_t)-1) {
bub->b_s_idx.a[v] <<= 32;
bub->b_s_idx.a[v] |= i;
} else if((bub->b_s_idx.a[v] & 0xffffffff00000000) == 0xffffffff00000000) {
bub->b_s_idx.a[v] <<= 32;
bub->b_s_idx.a[v] |= i;
}
}
for (i = 0; i < ug->g->n_seq; i++) {
if(bub->index[i] == M_het(*bub)) bub->index[i] = P_het(*bub);
}
bub->b_g = NULL;
bub->b_ug = NULL;
build_bub_graph(ug, bub);
///make het nodes to be hom
ma_utg_t *u = NULL;
for (i = 0; i < bub->b_ug->u.n; i++) {
u = &(bub->b_ug->u.a[i]);
if(u->n == 0) continue;
for (k = 0; k < u->n; k++) {
reset_inner_bub_het_poy(sg, ug, bub, u->a[k]>>33, tLen, &b, m_het_cov, m_hom_cov, r_het_flag, sources, ruIndex);
}
}
free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a);
for (i = 0; i < ug->g->n_seq; i++) {
if(IF_HOM(i, *bub)) continue;
if(deter_unitig_het_fly(ug, i, sg, het_cov*1.15, sources, ruIndex, r_het_flag)) continue;
bub->index[i] = (uint32_t)-1;
}
// fprintf(stderr, "-bub->index[18759]: %u, bub->num.n: %u\n", (uint32_t)bub->index[18759], bub->num.n);
}
void print_bubbles(ma_ug_t* ug, bubble_type* bub, kvec_pe_hit* hits, hc_links* link, ha_ug_index* idx)
{
uint64_t tLen, t_utg, i, k;
hic.h
+2
View File
@@ -95,6 +95,8 @@ void destory_bubbles(bubble_type* bub);
void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_trans_t *ref);
void identify_bubbles_recal(asg_t* sg, ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, ma_hit_t_alloc* sources, R_to_U* ruIndex,
kv_u_trans_t *ref);
void identify_bubbles_recal_poy(asg_t* sg, ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, ma_hit_t_alloc* sources, R_to_U* ruIndex,
kv_u_trans_t *ref);
void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub);
uint32_t connect_bub_occ(bubble_type* bub, uint32_t root_id, uint32_t check_het);
void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, uint32_t check_het);
inter.cpp
+1 -1
View File
@@ -19599,7 +19599,7 @@ int32_t mini_cut, int32_t chain_cut, int32_t keep_unsymm_arc, ma_ug_t *ug, asg_t
if(bub) {
opt->bub = bub;
} else {
opt->bub = gen_bubble_chain(sg, ug, uopt, &bf); free(bf);
opt->bub = gen_bubble_chain(sg, ug, uopt, &bf, 0); free(bf);
}
}