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

Merge pull request #368 from chhylp123/hifiasm_dev_debug

Browse Source 
fixed het/hom node bugs
This commit is contained in:
chhylp123
2022-12-08 03:07:43 -05:00
committed by GitHub
parent 6743dcf204 a674110b23
commit 34948a093c
6 changed files with 67 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.cpp
+10 -2
View File
@@ -49,6 +49,8 @@ static ko_longopt_t long_options[] = {
{ "kpt-rate", ko_required_argument, 335},
{ "ul-rate", ko_required_argument, 336},
{ "dbg-het-cnt", ko_no_argument, 337},
{ "ul-tip", ko_required_argument, 338},
{ "low-het", ko_no_argument, 339},
{ 0, 0, 0 }
};
@@ -146,8 +148,10 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " Ultra-Long-integration (beta):\n");
fprintf(stderr, " --ul FILEs file names of Ultra-Long reads [r1.fq,r2.fq,...]\n");
fprintf(stderr, " --ul-rate FLOAT\n");
fprintf(stderr, " similarity threshold for UL-to-HiFi alignment [%.3g]\n", asm_opt->ul_error_rate);
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, " --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");
fprintf(stderr, "See `https://hifiasm.readthedocs.io/en/latest/' or `man ./hifiasm.1' for complete documentation.\n");
@@ -204,6 +208,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->min_overlap_Len = 50;
asm_opt->min_overlap_coverage = 0;
asm_opt->max_short_tip = 3;
asm_opt->max_short_ul_tip = 6;
asm_opt->min_cnt = 2;
asm_opt->mid_cnt = 5;
asm_opt->purge_level_primary = 3;
@@ -246,6 +251,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->ul_error_rate_hpc = 0.2;
asm_opt->ul_ec_round = 3;
asm_opt->is_dbg_het_cnt = 0;
asm_opt->is_low_het_ul = 0;
}
void destory_enzyme(enzyme* f)
@@ -763,6 +769,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 335) asm_opt->kpt_rate = atof(opt.arg);
else if (c == 336) asm_opt->ul_error_rate = atof(opt.arg);
else if (c == 337) asm_opt->is_dbg_het_cnt = 1;
else if (c == 338) asm_opt->max_short_ul_tip = atol(opt.arg);
else if (c == 339) asm_opt->is_low_het_ul = 1;
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
+3 -1
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.18.1-r466"
#define HA_VERSION "0.18.2-r467"
#define VERBOSE 0
@@ -79,6 +79,7 @@ typedef struct {
int min_overlap_Len;
int min_overlap_coverage;
int max_short_tip;
int max_short_ul_tip;
int min_cnt;
int mid_cnt;
int purge_level_primary;
@@ -128,6 +129,7 @@ typedef struct {
double ul_error_rate, ul_error_rate_low, ul_error_rate_hpc;
int32_t ul_ec_round;
uint8_t is_dbg_het_cnt;
uint8_t is_low_het_ul;
} hifiasm_opt_t;
extern hifiasm_opt_t asm_opt;
Overlaps.cpp
+13 -13
View File
@@ -3594,7 +3594,7 @@ uint32_t startNode, uint32_t endNode, int max_dist, buf_t* bub)
{
Ns_first[0] = av[0].v; Ns_first[1] = av[1].v;
f1 = asg_bub_end_finder_with_del_advance(g, Ns_first, 2, max_dist,
bub, 0, (u_int32_t)-1, &convex1);
bub, 0, (uint32_t)-1, &convex1);
l1 = min_thres + 10;
}
@@ -3602,7 +3602,7 @@ uint32_t startNode, uint32_t endNode, int max_dist, buf_t* bub)
{
Ns_second[0] = aw[0].v; Ns_second[1] = aw[1].v;
f2 = asg_bub_end_finder_with_del_advance(g, Ns_second, 2, max_dist,
bub, 0, (u_int32_t)-1, &convex2);
bub, 0, (uint32_t)-1, &convex2);
l2 = min_thres + 10;
}
@@ -4034,7 +4034,7 @@ uint32_t startNode, uint32_t endNode, int max_dist, buf_t* bub)
{
Ns_first[0] = Nodes1[0]^1; Ns_first[1] = Nodes2[1]^1;
f1 = asg_bub_end_finder_with_del_advance(g, Ns_first, 2, max_dist,
bub, 0, (u_int32_t)-1, &convex1);
bub, 0, (uint32_t)-1, &convex1);
l1 = min_thres + 10;
}
@@ -4042,7 +4042,7 @@ uint32_t startNode, uint32_t endNode, int max_dist, buf_t* bub)
{
Ns_second[0] = Nodes2[0]^1; Ns_second[1] = Nodes1[1]^1;
f2 = asg_bub_end_finder_with_del_advance(g, Ns_second, 2, max_dist,
bub, 0, (u_int32_t)-1, &convex2);
bub, 0, (uint32_t)-1, &convex2);
l2 = min_thres + 10;
}
@@ -6896,7 +6896,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex)
kv_push(uint32_t, b_r, w);
aw = asg_arc_a(g, w);
min_edge = (u_int32_t)-1;
min_edge = (uint32_t)-1;
for (t = 0; t < nw; t++)
{
if(aw[t].del) continue;
@@ -6996,7 +6996,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex)
}
//check the length
if(to_del_l == 0 && asg_bub_end_finder_with_del_advance(g,
b_f.a, b_f.n, max_dist, &bub, 0, (u_int32_t)-1, &sink)==1)
b_f.a, b_f.n, max_dist, &bub, 0, (uint32_t)-1, &sink)==1)
{
to_del_l = 1;
}
@@ -7143,7 +7143,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex)
kv_push(uint32_t, b_r, w);
aw = asg_arc_a(g, w);
min_edge = (u_int32_t)-1;
min_edge = (uint32_t)-1;
for (t = 0; t < nw; t++)
{
if(aw[t].del) continue;
@@ -7187,7 +7187,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex)
kv_push(uint32_t, b_r, w);
aw = asg_arc_a(g, w);
min_edge = (u_int32_t)-1;
min_edge = (uint32_t)-1;
for (t = 0; t < nw; t++)
{
if(aw[t].del) continue;
@@ -7288,7 +7288,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex)
}
//check the length
if(to_del_l == 0 && asg_bub_end_finder_with_del_advance(g,
b_f.a, b_f.n, max_dist, &bub, 0, (u_int32_t)-1, &sink)==1)
b_f.a, b_f.n, max_dist, &bub, 0, (uint32_t)-1, &sink)==1)
{
to_del_l = 1;
}
@@ -7431,7 +7431,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex)
kv_push(uint32_t, b_r, w);
aw = asg_arc_a(g, w);
min_edge = (u_int32_t)-1;
min_edge = (uint32_t)-1;
for (t = 0; t < nw; t++)
{
if(aw[t].del) continue;
@@ -7475,7 +7475,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex)
kv_push(uint32_t, b_r, w);
aw = asg_arc_a(g, w);
min_edge = (u_int32_t)-1;
min_edge = (uint32_t)-1;
for (t = 0; t < nw; t++)
{
if(aw[t].del) continue;
@@ -7575,7 +7575,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex)
}
//check the length
if(to_del_l == 0 && asg_bub_end_finder_with_del_advance(g,
b_f.a, b_f.n, max_dist, &bub, 0, (u_int32_t)-1, &sink)==1)
b_f.a, b_f.n, max_dist, &bub, 0, (uint32_t)-1, &sink)==1)
{
to_del_l = 1;
}
@@ -31884,7 +31884,7 @@ int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file)
nopt.cov = coverage_cut; nopt.b_mask_t = b_mask_t;
nopt.max_hang = max_hang_length; nopt.mini_ovlp = mini_overlap_length;
ul_realignment_gfa(uopt, *sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
asm_opt.max_short_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file, &nopt);;
asm_opt.max_short_tip, asm_opt.max_short_ul_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file, &nopt);;
// ma_ug_t *iug = ul_realignment_gfa(uopt, *sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
// asm_opt.max_short_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file);
// asg_t *ng = gen_ng(iug, *sg, uopt, coverage_cut, ruIndex, 256);
gfa_ut.cpp
+37 -4
View File
@@ -14385,9 +14385,12 @@ void u2g_hybrid_detan_iter(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint
asg64_v b64, ub64; kv_init(b64); kv_init(ub64);
asg64_v tx = {0,0,0}, tb = {0,0,0}, *ob = NULL, *ub = NULL;
ob = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); ob->n = ub->n = 0;
// fprintf(stderr, "\n[M::%s::] asm_opt.is_low_het_ul::%u, max_ext::%u\n",
// __func__, asm_opt.is_low_het_ul, max_ext);
// prt_usg_t(uidx, ng, "ng0");
for (k = 0; k < clean_round; k++) {
ncut += ug_ext_strict(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64);
// if(asm_opt.is_low_het_ul) break;
ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 48, 1);
// prt_usg_t(uidx, ng, "ng_python");
ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 16, 1);
@@ -14403,7 +14406,7 @@ void u2g_hybrid_detan_iter(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint
// prt_usg_t(uidx, ng, "ng.db");
// ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 10);
// ncut += ug_ext_strict(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64);
// prt_usg_t(uidx, ng, "ng1");
// prt_usg_t(uidx, ng, "ng.db");
if(ncut) {
usg_cleanup(ng); usg_arc_cut_tips(ng, max_ext, 1, ub);
}
@@ -16479,8 +16482,38 @@ void destroy_ul_resolve_t(ul_resolve_t *uidx)
free(uidx);
}
static void clear_ma_hit_t_alloc(void *data, long i, int tid)
{
ma_hit_t_alloc *src = (ma_hit_t_alloc *)data;
ma_hit_t_alloc *z = &(src[i]); uint32_t k;
for (k = 0; k < z->length; k++) z->buffer[k].del = 0;
}
static void reset_ma_sub_t(void *data, long i, int tid)
{
sset_aux *s = (sset_aux *)data;
if(s->g && s->g->seq[i].del) s->cov[i].del = 1;
else s->cov[i].del = 0;
}
void renew_R_to_U(asg_t *ng, ma_hit_t_alloc* src, ma_hit_t_alloc* r_src, int64_t n_read, ma_sub_t *coverage_cut,
R_to_U* ruIndex, int max_hang, int min_ovlp)
{
sset_aux s; memset(&s, 0, sizeof(s));
kt_for(asm_opt.thread_num, clear_ma_hit_t_alloc, src, n_read);
kt_for(asm_opt.thread_num, clear_ma_hit_t_alloc, r_src, n_read);
s.g = NULL; s.cov = coverage_cut;
kt_for(asm_opt.thread_num, reset_ma_sub_t, &s, n_read);
ma_hit_contained_advance(src, n_read, coverage_cut, ruIndex, max_hang, min_ovlp);
// s.g = ng; s.cov = coverage_cut;
// kt_for(asm_opt.thread_num, reset_ma_sub_t, &s, n_read);
}
void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio,
double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file,
double max_ovlp_drop_ratio, int64_t max_tip, int64_t max_ul_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file,
ul_renew_t *ropt)
{
uint64_t i; uint8_t *r_het = NULL; bubble_type *bub = NULL; ulg_opt_t uu;
@@ -16502,7 +16535,7 @@ ul_renew_t *ropt)
// print_raw_uls_seq(uidx, asm_opt.output_file_name);
// print_raw_uls_aln(uidx, asm_opt.output_file_name);
ul_re_correct(uidx, 3);
init_ulg_opt_t(&uu, uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_tip<<1, b_mask_t, is_trio);
init_ulg_opt_t(&uu, uopt, clean_round, 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);
// print_ul_alignment(init_ug, &UL_INF, 47072, "after-3");
@@ -16529,7 +16562,7 @@ ul_renew_t *ropt)
asg_destroy((*(ropt->sg))); (*(ropt->sg)) = ng;
(*(ropt->src)) = R_INF.paf; (*(ropt->r_src)) = R_INF.reverse_paf;
(*(ropt->n_read)) = R_INF.total_reads; (*(ropt->readLen)) = R_INF.read_length;
ma_hit_contained_advance((*(ropt->src)), (*(ropt->n_read)), (*(ropt->cov)), ropt->ruIndex, ropt->max_hang, ropt->mini_ovlp);
renew_R_to_U(ng, (*(ropt->src)), (*(ropt->r_src)), (*(ropt->n_read)), (*(ropt->cov)), ropt->ruIndex, ropt->max_hang, ropt->mini_ovlp);
post_rescue(uopt, (*(ropt->sg)), (*(ropt->src)), (*(ropt->r_src)), ropt->ruIndex, ropt->b_mask_t, 0);
// print_raw_uls_aln(uidx, asm_opt.output_file_name);
// exit(0);
gfa_ut.h
+2 -1
View File
@@ -6,6 +6,7 @@
typedef struct {
asg_t *g;
ma_hit_t_alloc *src;
ma_sub_t *cov;
R_to_U* ruIndex;
int64_t max_hang;
int64_t min_ovlp;
@@ -25,7 +26,7 @@ void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float o
uint32_t asg_cut_large_indel(asg_t *g, asg64_v *in, int32_t max_ext, float ou_rat, uint32_t is_ou);
uint32_t asg_cut_semi_circ(asg_t *g, uint32_t lim_len, uint32_t is_clean);
void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio,
double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, ul_renew_t *ropt);
double max_ovlp_drop_ratio, int64_t max_tip, int64_t max_ul_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, ul_renew_t *ropt);
void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t max_hang, int64_t min_ovlp, int64_t ul_occ);
void normalize_gou(asg_t *g);
void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd);
hic.h
+2
View File
@@ -110,4 +110,6 @@ void set_utg_by_dis(uint32_t v, pdq* pq, asg_t *g, kvec_t_u32_warp *res, uint32_
void dedup_hits(kvec_pe_hit* hits, uint64_t is_dup);
void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, uint32_t is_poy, kvec_pe_hit **rhits);
spg_t *hic_pre_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, kvec_pe_hit **rhits);
void prt_bubble_gfa_adv(FILE *fp, bubble_type *bub, const char* utg_pre, const char* bub_pre, const char* chain_pre);
#endif