From 8a0ec3d5750785d327e5d34c6b6fd34f25dd3605 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Sun, 12 Dec 2021 20:59:48 -0500 Subject: [PATCH] change ug extraction --- Overlaps.cpp | 3 +- Process_Read.cpp | 15 ++++++- htab.cpp | 109 +++++++++++++++++++++++++++++++++-------------- htab.h | 2 + inter.cpp | 74 +++++++++++++++++++++++++++++++- 5 files changed, 167 insertions(+), 36 deletions(-) diff --git a/Overlaps.cpp b/Overlaps.cpp index e5832bf..5980665 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -31117,7 +31117,7 @@ float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr, ma_sub_t **coverage_cut_ptr, int debug_g) { - // if(asm_opt.ar) create_ul_info(); + if(asm_opt.ar) create_ul_info(); @@ -31400,6 +31400,7 @@ long long bubble_dist, int read_graph, int write) init_R_to_U(&ruIndex, n_read); asg_t *sg = NULL; ma_sub_t* coverage_cut = NULL; + init_aux_table(); ///actually min_thres = asm_opt.max_short_tip + 1 there are asm_opt.max_short_tip reads min_thres = asm_opt.max_short_tip + 1; diff --git a/Process_Read.cpp b/Process_Read.cpp index cf70f84..d9dcde3 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -881,7 +881,6 @@ void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, } } - void retrieve_ul_t(UC_Read* r, all_ul_t *ref, uint64_t ID, uint8_t strand) { ul_vec_t *p = &(ref->a[ID]); uc_block_t *b = NULL; @@ -931,4 +930,18 @@ void retrieve_ul_t(UC_Read* r, all_ul_t *ref, uint64_t ID, uint8_t strand) { } for (k = 0; k < (int64_t)p->N_site.n; k++) r->seq[r->length - p->N_site.a[k] - 1] = 'N'; } +} + +void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand) +{ + if(i_r) { + + } + + + if (r->length + 4 > r->size) + { + r->size = r->length + 4; + r->seq = (char*)realloc(r->seq,sizeof(char)*(r->size)); + } } \ No newline at end of file diff --git a/htab.cpp b/htab.cpp index 16d7fee..97b3d18 100644 --- a/htab.cpp +++ b/htab.cpp @@ -8,6 +8,7 @@ #include "kseq.h" #include "ksort.h" #include "htab.h" +#include "Process_Read.h" #define YAK_COUNTER_BITS 12 #define YAK_N_COUNTS (1<p = p;\ s->n_seq0 = p->n_seq;\ + s->uq = p->opt->uq;\ if (p->rs_in && (p->flag & HAF_RS_READ)) {\ while (p->n_seq < p->rs_in->total_reads) {\ if ((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU) {\ @@ -721,7 +723,7 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for break;\ }\ } else if(p->us_in) {\ - ma_utg_t *u; s->uq = 1;\ + ma_utg_t *u;\ while (p->n_seq < p->us_in->n) {\ u = &(p->us_in->a[p->n_seq]);\ if (s->n_seq == s->m_seq) {\ @@ -730,7 +732,8 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for REALLOC(s->seq, s->m_seq);\ }\ MALLOC(s->seq[s->n_seq], u->len);\ - memcpy(s->seq[s->n_seq], u->s, u->len);\ + if(u->s) memcpy(s->seq[s->n_seq], u->s, u->len);\ + else memcpy(s->seq[s->n_seq], u->s, u->len);\ s->len[s->n_seq++] = u->len;\ ++p->n_seq;\ s->sum_len += u->len;\ @@ -905,7 +908,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt pl_data_t pl; gzFile fp = 0; memset(&pl, 0, sizeof(pl_data_t)); - pl.n_seq = *n_seq; + pl.n_seq = *n_seq; if(ug_rs) { pl.us_in = us; } else if (read_rs) { @@ -950,7 +953,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt return pl.ct; } -ha_ct_t *ha_count(const hifiasm_opt_t *asm_o, int flag, int HPC, int k, int w, ha_pt_t *p0, const void *flt_tab, All_reads *rs, ma_utg_v *us, int keep_adapter, int *low_freq) +ha_ct_t *ha_count(const hifiasm_opt_t *asm_o, int flag, int HPC, int k, int w, ha_pt_t *p0, const void *flt_tab, All_reads *rs, ma_utg_v *us, int keep_adapter, int *low_freq, int unique_only) { int i; int64_t n_seq = 0; @@ -976,6 +979,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_o, int flag, int HPC, int k, int w, h opt.n_thread = asm_o->thread_num; opt.adaLen = (keep_adapter? asm_o->adapterLen : 0); opt.min_rcnt = (low_freq?*low_freq:-1); + opt.uq = (unique_only?1:0); ///asm_opt->num_reads is the number of fastq files for (i = n_bs = 0; i < (us?1:asm_o->num_reads); ++i){ h = yak_count(&opt, asm_o->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq); @@ -1062,39 +1066,35 @@ void debug_ct_index(void* q_ct_idx, void* r_ct_idx) /************************* * High-level interfaces * *************************/ +void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff) +{ + yak_ft_t *flt_tab; + ha_ct_t *h; + ///HAF_COUNT_EXACT ---> no bf; HAF_COUNT_ALL ---> no minimizer + h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ, !(asm_opt->flag&HA_F_NO_HPC), k, w, NULL, NULL, NULL, us, 0, NULL, 0); + + // cutoff = (int)(asm_opt->hom_cov * asm_opt->high_factor); + if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1; + ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num); + flt_tab = gen_hh(h, asm_opt->max_kmer_cnt); + ha_ct_destroy(h); + fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> filtered out %ld k-mers occurring %d or more times\n", __func__, + yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb(), (long)kh_size(flt_tab), cutoff); + return (void*)flt_tab; +} void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq) { yak_ft_t *flt_tab; ha_ct_t *h; ///HAF_COUNT_EXACT ---> no bf; HAF_COUNT_ALL ---> no minimizer - h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, is_HPC, k, w, NULL, NULL, NULL, us, 0, NULL); + h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, is_HPC, k, w, NULL, NULL, NULL, us, 0, NULL, 0); ha_ct_shrink(h, min_freq, max_freq>YAK_MAX_COUNT-1?YAK_MAX_COUNT-1:max_freq, asm_opt->thread_num); flt_tab = gen_hh(h, YAK_MAX_COUNT); ha_ct_destroy(h); return (void*)flt_tab; } - -ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int is_HPC, int k, int w, int min_freq) -{ - ha_ct_t *ct; - ha_pt_t *pt; - ///HAF_COUNT_EXACT: no bf - ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, NULL, flt_tab, NULL, us, 0, NULL); - fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__, - yak_realtime(), yak_cpu_usage(), (long)ct->tot); - ///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen - ha_ct_shrink(ct, min_freq, YAK_MAX_COUNT - 1, asm_opt->thread_num); - - pt = ha_pt_gen(ct, asm_opt->thread_num, 1); - ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, pt, flt_tab, NULL, us, 0, NULL); - //ha_pt_sort(pt, asm_opt->thread_num); - fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__, - yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos); - return pt; -} - void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode) { yak_ft_t *flt_tab; @@ -1102,7 +1102,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i int peak_hom, peak_het, cutoff = YAK_MAX_COUNT - 1, ex_flag = 0; if(is_hp_mode) ex_flag = HAF_RS_READ|HAF_SKIP_READ; ha_ct_t *h; - h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, NULL, rs, NULL, 1, NULL); + h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, NULL, rs, NULL, 1, NULL, 0); if((asm_opt->flag & HA_F_VERBOSE_GFA)) { write_ct_index((void*)h, asm_opt->output_file_name); @@ -1130,16 +1130,61 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i return (void*)flt_tab; } +ha_pt_t *ha_pt_ul_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int k, int w, int cutoff) +{ + ha_ct_t *ct; + ha_pt_t *pt; + ///HAF_COUNT_EXACT: no bf + ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, !(asm_opt->flag&HA_F_NO_HPC), k, w, NULL, flt_tab, NULL, us, 0, NULL, 0); + fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__, + yak_realtime(), yak_cpu_usage(), (long)ct->tot); + ///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen + if (flt_tab == 0) { + if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1; + ha_ct_shrink(ct, 2, cutoff, asm_opt->thread_num); + } else { + ///Note: here is just to remove minimizer appearing YAK_MAX_COUNT times + ///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen + ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num); + } + + pt = ha_pt_gen(ct, asm_opt->thread_num, 1); + ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, !(asm_opt->flag&HA_F_NO_HPC), k, w, pt, flt_tab, NULL, us, 0, NULL, 0); + //ha_pt_sort(pt, asm_opt->thread_num); + fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__, + yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos); + return pt; +} + +ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int is_HPC, int k, int w, int min_freq) +{ + ha_ct_t *ct; + ha_pt_t *pt; + ///HAF_COUNT_EXACT: no bf + ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, NULL, flt_tab, NULL, us, 0, NULL, 1); + fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__, + yak_realtime(), yak_cpu_usage(), (long)ct->tot); + ///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen + ha_ct_shrink(ct, min_freq, YAK_MAX_COUNT - 1, asm_opt->thread_num); + + pt = ha_pt_gen(ct, asm_opt->thread_num, 1); + ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, pt, flt_tab, NULL, us, 0, NULL, 1); + //ha_pt_sort(pt, asm_opt->thread_num); + fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__, + yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos); + return pt; +} + ha_pt_t *ha_pt_gen_dp(const hifiasm_opt_t *asm_opt, ha_ct_t *ct, int flag, int n_thread, const void *flt_tab, All_reads *rs, int peak_hom, int peak_het) { int low_freq = mz_low_b(peak_hom, peak_het); ha_pt_t *pt = ha_pt_gen_count(ct, n_thread); ///key = cnt, val = 0 - ha_count(asm_opt, HAF_COUNT_EXACT|HAF_COUNT_REFINE|flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, &low_freq); + ha_count(asm_opt, HAF_COUNT_EXACT|HAF_COUNT_REFINE|flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, &low_freq, 0); uint64_t occ = ha_pt_shrink(pt, n_thread); if(flag&HAF_RS_WRITE_LEN) flag -= HAF_RS_WRITE_LEN; if(flag&HAF_RS_WRITE_SEQ) flag -= HAF_RS_WRITE_SEQ; flag |= HAF_RS_READ; pt->tot_pos = 0; - ha_count(asm_opt, HAF_COUNT_EXACT|flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, NULL); + ha_count(asm_opt, HAF_COUNT_EXACT|flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, NULL, 0); // fprintf(stderr, "[M::%s::] counted %lu distinct minimizer k-mers\n", __func__, pt->tot); // fprintf(stderr, "[M::%s::] collected %lu minimizers\n\n\n", __func__, pt->tot_pos); assert(occ == pt->tot_pos); @@ -1163,7 +1208,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f } if(is_hp_mode) extra_flag1 |= HAF_SKIP_READ, extra_flag2 |= HAF_SKIP_READ; - ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, flt_tab, rs, NULL, 1, NULL); + ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, flt_tab, rs, NULL, 1, NULL, 0); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__, yak_realtime(), yak_cpu_usage(), (long)ct->tot); ha_ct_hist(ct, cnt, asm_opt->thread_num); @@ -1189,7 +1234,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f { fprintf(stderr, "[M::%s::] counting in normal mode\n", __func__); pt = ha_pt_gen(ct, asm_opt->thread_num, 0); - ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, NULL); + ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, NULL, 0); assert((uint64_t)tot_cnt == pt->tot_pos); } else diff --git a/htab.h b/htab.h index 41549ff..ec9a22c 100644 --- a/htab.h +++ b/htab.h @@ -72,11 +72,13 @@ extern void *ha_flt_tab_hp; extern ha_pt_t *ha_idx_hp; extern void *ha_ct_table; +void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff); void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq); void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode); int32_t ha_ft_cnt(const void *hh, uint64_t y); void ha_ft_destroy(void *h); +ha_pt_t *ha_pt_ul_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int k, int w, int cutoff); ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int is_HPC, int k, int w, int min_freq); ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov); void ha_pt_destroy(ha_pt_t *h); diff --git a/inter.cpp b/inter.cpp index f723f8a..e3ab5ea 100644 --- a/inter.cpp +++ b/inter.cpp @@ -322,6 +322,13 @@ void init_mg_opt(mg_idxopt_t *opt, int is_HPC, int k, int w, int hap_n) opt->pri_ratio = 0.8f; } +void uidx_l_build(ma_ug_t *ug, mg_idxopt_t *opt, int cutoff) +{ + ha_flt_tab = ha_ft_ul_gen(&asm_opt, &(ug->u), opt->k, opt->w, cutoff); + ha_idx = ha_pt_ul_gen(&asm_opt, ha_flt_tab, &(ug->u), opt->k, opt->w, cutoff); + fprintf(stderr, "[M::%s] Index has been built.\n", __func__); +} + void uidx_build(ma_ug_t *ug, mg_idxopt_t *opt) { int flag = asm_opt.flag; @@ -2345,7 +2352,6 @@ int scall_ul_pipeline(uldat_t* sl, const enzyme *fn) { double index_time = yak_realtime(); int i; - init_aux_table(); init_all_ul_t(&UL_INF, &R_INF); for (i = 0; i < fn->n; i++){ @@ -2885,10 +2891,74 @@ int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const asg_t *rg, const enz return 1; } +ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, uint64_t n_read, int64_t readLen) +{ + uint64_t i, k, qn, tn; + int32_t r; asg_arc_t t, *p = NULL; + asg_t *rg = asg_init(); + + rg->m_seq = n_read; MALLOC(rg->seq, rg->m_seq); + for (i = 0; i < n_read; ++i) rg->seq[i].len = Get_READ_LENGTH(R_INF, i); + + for (i = 0; i < n_read; i++) { + if(rg->seq[i].del) continue; + for (k = 0; k < src[i].length; k++) { + if(!src[i].buffer[k].el) continue; + qn = Get_qn(src[i].buffer[k]); tn = Get_tn(src[i].buffer[k]); + if(rg->seq[qn].del || rg->seq[tn].del) continue; + r = ma_hit2arc(&(src[i].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + if (r == MA_HT_QCONT) { + rg->seq[qn].del = 1; + } else if(r == MA_HT_TCONT) { + rg->seq[tn].del = 1; + } + if(rg->seq[i].del) break; + } + } + + for (i = 0; i < n_read; i++) { + if(rg->seq[i].del) continue; + for (k = 0; k < src[i].length; k++) { + if(!src[i].buffer[k].el) continue; + qn = Get_qn(src[i].buffer[k]); tn = Get_tn(src[i].buffer[k]); + if(rg->seq[qn].del || rg->seq[tn].del) continue; + r = ma_hit2arc(&(src[i].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + if (r >= 0) { + p = asg_arc_pushp(rg); + *p = t; + } + } + } + + asg_cleanup(rg); asg_symm(rg); + + + + + ma_ug_t *ug = NULL; + ug = ma_ug_gen(rg); + + + + + + + asg_destroy(rg); ma_ug_destroy(ug); +} + void ul_load(const ug_opt_t *uopt) { fprintf(stderr, "[M::%s::] ==> UL\n", __func__); mg_idxopt_t opt; - init_mg_opt(&opt, 0, 19, 10, 4095); + ma_ug_t *ug = NULL; + int cutoff = asm_opt.hom_cov * asm_opt.high_factor; + init_aux_table(); + init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff); + int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name) : 0); + if(exist == 0) uidx_l_build(ug, &opt, cutoff); + if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name); + + + ul_v_call(&opt, uopt, /**rg**/NULL, asm_opt.ar, ha_flt_tab, ha_idx, /**ug**/NULL); } \ No newline at end of file