From 9ea4f6917f5e165d01ce1dea59523ed84d2e0c90 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Thu, 6 May 2021 07:59:51 -0400 Subject: [PATCH] resolve centromere --- Overlaps.cpp | 60 +++++++++++-------------- Overlaps.h | 2 + Purge_Dups.cpp | 46 ++++++++++++++++--- hic.cpp | 119 ++++++++++++++++++++++++++++++++++++++----------- hic.h | 3 ++ 5 files changed, 163 insertions(+), 67 deletions(-) diff --git a/Overlaps.cpp b/Overlaps.cpp index 0c03a76..552462b 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -8994,8 +8994,8 @@ UC_Read* r_read, UC_Read* q_read, int max_hang, int min_ovlp, kvec_asg_arc_t_war // generate unitig sequences -int ma_ug_seq(ma_ug_t *g, asg_t *read_g, All_reads *RNF, ma_sub_t *coverage_cut, -ma_hit_t_alloc* sources, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, kvec_asg_arc_t_warp *E) +int ma_ug_seq(ma_ug_t *g, asg_t *read_g, ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, +kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, kvec_asg_arc_t_warp *E, uint32_t is_polish) { UC_Read g_read; init_UC_Read(&g_read); @@ -9013,8 +9013,12 @@ ma_hit_t_alloc* sources, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, for (i = 0; i < g->u.n; ++i) { ma_utg_t *u = &g->u.a[i]; if(u->m == 0) continue; - polish_unitig(u, read_g, sources, coverage_cut, edge, max_hang, min_ovlp, E); - polish_unitig_advance(u, read_g, RNF, sources, coverage_cut, edge, &g_read, &tmp, max_hang, min_ovlp, E); + if(is_polish) + { + polish_unitig(u, read_g, sources, coverage_cut, edge, max_hang, min_ovlp, E); + polish_unitig_advance(u, read_g, &R_INF, sources, coverage_cut, edge, &g_read, &tmp, max_hang, min_ovlp, E); + } + g->g->seq[i].len = u->len; uint32_t l = 0; @@ -9031,12 +9035,12 @@ ma_hit_t_alloc* sources, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, if(eLen == 0) continue; if(rId < read_g->r_seq) { - recover_UC_Read(&g_read, RNF, rId); + recover_UC_Read(&g_read, &R_INF, rId); } else { recover_fake_read(&g_read, &tmp, &(read_g->F_seq[rId-read_g->r_seq]), - RNF, coverage_cut); + &R_INF, coverage_cut); } readS = g_read.seq + coverage_cut[rId].s; @@ -12153,7 +12157,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) ma_ug_t *ug = NULL; ug = ma_ug_gen(sg); - ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0); + ma_ug_seq(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 1); fprintf(stderr, "Writing raw unitig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+25); @@ -12493,7 +12497,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp, kvec_asg_a 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); + ma_ug_seq(ug, sg, coverage_cut, sources, new_rtg_edges, max_hang, min_ovlp, 0, 1); char* gfa_name = (char*)malloc(strlen(output_file_name)+35); @@ -12723,7 +12727,7 @@ long long gap_fuzz, bub_label_t* b_mask_t) ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); new_rtg_edges.a.n = 0; - ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, &d_edges);///polish + ma_ug_seq(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, &d_edges, 1);///polish hap_cov_t *cov = NULL; @@ -12763,7 +12767,7 @@ long long gap_fuzz, bub_label_t* b_mask_t) } char* gfa_name = (char*)malloc(strlen(output_file_name)+25); - sprintf(gfa_name, "%s.clean_d_utg.noseq.gfa", output_file_name); + sprintf(gfa_name, "%s.pre.clean_d_utg.noseq.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); fclose(output_file); @@ -12774,6 +12778,14 @@ long long gap_fuzz, bub_label_t* b_mask_t) if(cov) destory_hap_cov_t(&cov); if(t_ch) destory_trans_chain(&t_ch); + + gfa_name = (char*)malloc(strlen(output_file_name)+25); + sprintf(gfa_name, "%s.after.clean_d_utg.noseq.gfa", output_file_name); + output_file = fopen(gfa_name, "w"); + ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); + fclose(output_file); + free(gfa_name); + ma_ug_destroy(ug); kv_destroy(new_rtg_edges.a); @@ -13752,7 +13764,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) } } - ma_ug_seq(ug, read_g, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0); + ma_ug_seq(ug, read_g, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 1); fprintf(stderr, "Writing raw unitig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+25); @@ -16160,7 +16172,7 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench, ///debug_utg_graph(ug, sg, 0, 0); ///debug_untig_length(ug, tipsLen, gfa_name); ///print_untig_by_read(ug, "m64011_190901_095311/125831121/ccs", 2310925, "end"); - ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0); + ma_ug_seq(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 1); if(is_bench) { free(gfa_name); @@ -18831,24 +18843,6 @@ long long asg_arc_del_simple_circle_untig(ma_hit_t_alloc* sources, ma_sub_t* cov } - -uint32_t* build_unitig_index(asg_t *sg, ma_ug_t *ug) -{ - if(sg == NULL && ug == NULL) return NULL; - uint32_t* index = (uint32_t*)calloc(sg->n_seq, sizeof(uint32_t)); - uint32_t i, j; - for (i = 0; i < ug->u.n; ++i) - { - ma_utg_t *u = &(ug->u.a[i]); - for (j = 0; j < u->n; j++) - { - index[(u->a[j]>>33)] = i; - } - } - return index; -} - - int double_check_tangle(uint32_t vBeg, uint32_t vEnd, uint32_t* u_vecs, uint32_t n, asg_t *nsg) { uint32_t v, nv, k, i, j; @@ -24049,7 +24043,7 @@ R_to_U* ruIndex, int max_hang, int min_ovlp) } - ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0); + ma_ug_seq(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 1); fprintf(stderr, "Writing processed unitig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+35); @@ -24100,7 +24094,7 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov 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); + ma_ug_seq(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 1); fprintf(stderr, "Writing primary contig GFA to disk... \n"); @@ -24144,7 +24138,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) // break_ug_contig(&ug, sg, &R_INF, coverage_cut, sources, ruIndex, &new_rtg_edges, max_hang, min_ovlp, asm_opt.b_low_cov); // } - ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0); + ma_ug_seq(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 1); fprintf(stderr, "Writing alternate contig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+35); diff --git a/Overlaps.h b/Overlaps.h index 761b26e..51b4d75 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1192,6 +1192,8 @@ ma_ug_t *ug, uint32_t flag, double overall_score, const char* cmd); int asg_arc_del_trans(asg_t *g, int fuzz); void kt_u_trans_t_idx(kv_u_trans_t *ta, uint32_t n); uint32_t get_u_trans_spec(kv_u_trans_t *ta, uint32_t qn, uint32_t tn, u_trans_t **r_a, uint32_t *occ); +int ma_ug_seq(ma_ug_t *g, asg_t *read_g, ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, +kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, kvec_asg_arc_t_warp *E, uint32_t is_polish); #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index 2e81c19..582d6c2 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -2242,10 +2242,10 @@ uint32_t* xBeg, uint32_t* xEnd, uint32_t* yBeg, uint32_t* yEnd) KRADIX_SORT_INIT(i32, int32_t, generic_key, sizeof(int32_t)) long long get_chain_score(ma_utg_t *xReads, asg_t *read_g, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* idx, -ma_hit_t_alloc* reverse_sources, long long xBegPos, long long xEndPos) +ma_hit_t_alloc* reverse_sources, long long xBegPos, long long xEndPos, uint32_t xUid, uint32_t yUid) { long long offset, r_beg, r_end, inp_beg, inp_end, hap_beg, hap_end, inp_match, hap_match, ovlp; - uint64_t i, k, rId; + uint64_t i, k, rId, all, found; idx->a.n = 0; for (i = k = 0; i < tailIndex->a.n; i++) @@ -2331,8 +2331,19 @@ ma_hit_t_alloc* reverse_sources, long long xBegPos, long long xEndPos) // if(inp_match > hap_match) fprintf(stderr, "ERROR3\n"); // fprintf(stderr, "tailIndex->a.n: %u, xReads->n: %u, total_match: %lld, hap_match: %lld, inp_match: %lld\n", // tailIndex->a.n, xReads->n, (xEndPos - xBegPos + 1), hap_match, inp_match); - - return ((double)(inp_match)*CHAIN_MATCH) - ((double)(hap_match-inp_match)*CHAIN_UNMATCH); + double base_w = 0, k_w = 1; + base_w = ((double)(inp_match)*CHAIN_MATCH) - ((double)(hap_match-inp_match)*CHAIN_UNMATCH); + if(base_w <= 0) return base_w; ///hard filter + + if(!count_unique_k_mers(xReads->s + xBegPos, xEndPos+1-xBegPos, xUid, yUid, &all, &found)) + { + return base_w; + } + + if(all > 0) k_w += ((double)(found)/(double)(all)); + + return (base_w*k_w)/2; + // return ((double)(inp_match)*CHAIN_MATCH) - ((double)(hap_match-inp_match)*CHAIN_UNMATCH); } @@ -2378,7 +2389,7 @@ kvec_t_i32_warp* prevIndex, long long* r_x_pos_beg, long long* r_x_pos_end, long hap_can->weight = xLeftMatch; hap_can->index_beg = xLeftTotal; hap_can->score = get_chain_score(xReads, read_g, u_buffer, tailIndex, prevIndex, reverse_sources, - (*r_x_pos_beg), (*r_x_pos_end)); + (*r_x_pos_beg), (*r_x_pos_end), xUid, yUid); if(hap_can->score <= 0) return (uint32_t)-1; return hap_can->index_end; } @@ -2541,7 +2552,14 @@ long long* r_y_pos_beg, long long* r_y_pos_end, float *sim) (*sim) = ((double)xLeftMatch)/((double)xLeftTotal); if(xLeftMatch == 0 || xLeftTotal == 0 || xLeftMatch <= xLeftTotal*Hap_rate) { - return NON_PLOID; + uint64_t all, found; + if(!count_unique_k_mers(xReads->s + (*r_x_pos_beg), (*r_x_pos_end)+1-(*r_x_pos_beg), + xUid, yUid, &all, &found)) + { + return NON_PLOID; + } + (*sim) = MAX((*sim), (((double)found)/((double)all))); + if(found == 0 || all == 0 || found <= all*Hap_rate) return NON_PLOID; } return PLOID; @@ -5263,7 +5281,12 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle } } - + if(just_coverage == 0) + { + ma_ug_seq(ug, read_g, coverage_cut, sources, edge, max_hang, min_ovlp, 0, 0); + } + init_ug_idx(ug, asm_opt.k_mer_length, asm_opt.polyploidy, 2, !just_coverage); + init_hap_alignment_struct_pip(&hap_buf, asm_opt.thread_num, nsg->n_seq, ug, read_g, sources, reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp, 0.1, &all_ovlp, cov); @@ -5392,4 +5415,13 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle else free(position_index); destory_hap_alignment_struct_pip(&hap_buf); destory_p_g_t(&pg); + if(just_coverage == 0) + { + des_ug_idx(); + for (i = 0; i < ug->u.n; i++) + { + free(ug->u.a[i].s); + ug->u.a[i].s = NULL; + } + } } diff --git a/hic.cpp b/hic.cpp index 097e7cc..ac582e2 100644 --- a/hic.cpp +++ b/hic.cpp @@ -153,7 +153,7 @@ typedef struct { uint64_t pre; uint64_t tot; uint64_t tot_pos; - uint64_t up_bound; + uint64_t up_bound, low_bound; hc_pt1_t* idx_buf; long double a, b, frac, max_d; } ha_ug_index; @@ -310,7 +310,8 @@ void print_debug_bubble_graph(bubble_type* bub, ma_ug_t* ug, const char *fn); void build_bub_graph(ma_ug_t* ug, bubble_type* bub); -void init_ha_ug_index_opt(ha_ug_index* idx, ma_ug_t *ug, int k, pldat_t* p) +void init_ha_ug_index_opt(ha_ug_index* idx, ma_ug_t *ug, int k, pldat_t* p, uint64_t up_occ, +uint64_t low_occ, uint64_t thread_num) { uint64_t i, n; for (idx->uID_bits=1; (uint64_t)(1<uID_bits)<(uint64_t)ug->u.n; idx->uID_bits++); @@ -324,7 +325,8 @@ void init_ha_ug_index_opt(ha_ug_index* idx, ma_ug_t *ug, int k, pldat_t* p) idx->tot = 1 << idx->pre; idx->tot_pos = 0; ///idx->up_bound = 1; - idx->up_bound = asm_opt.hap_occ; + idx->up_bound = up_occ; + idx->low_bound = low_occ; CALLOC(idx->idx_buf, idx->tot); for (i = 0; i < idx->tot; i++) { @@ -344,7 +346,7 @@ void init_ha_ug_index_opt(ha_ug_index* idx, ma_ug_t *ug, int k, pldat_t* p) kv_init(p->cnt[i].a); kv_init(p->buf[i].a); } - p->n_thread = asm_opt.thread_num; + p->n_thread = thread_num; } inline uint64_t get_k_direction(uint64_t x[4]) @@ -463,16 +465,17 @@ void test_unitig_index(ha_ug_index* idx, ma_ug_t *ug) fprintf(stderr, "[M::%s::%.3f] ==> Test has been passed\n", __func__, yak_realtime()-index_time); } -void hc_pt_t_gen_single(hc_pt1_t* pt, uint64_t* up_bound) +void hc_pt_t_gen_single(hc_pt1_t* pt, uint64_t* up_bound, uint64_t* low_bound) { khint_t k; uint64_t c; - if(up_bound) + if(up_bound || low_bound) { for (k = 0; k != kh_end(pt->h); ++k) { if (kh_exist(pt->h, k)) { - if(kh_val(pt->h, k) > (*up_bound)) + if((up_bound && kh_val(pt->h, k) > (*up_bound)) || + (low_bound && kh_val(pt->h, k) < (*low_bound))) { kh_val(pt->h, k) = 0; kh_key(pt->h, k) = (kh_key(pt->h, k)&HIC_KEY_MODE)| @@ -614,7 +617,7 @@ void hc_pt_t_gen(ha_ug_index* idx, pldat_t* pl) uint64_t i; for (i = 0; i < idx->tot; i++) { - hc_pt_t_gen_single(&(idx->idx_buf[i]), &(idx->up_bound)); + hc_pt_t_gen_single(&(idx->idx_buf[i]), &(idx->up_bound), &(idx->low_bound)); } } else @@ -785,12 +788,12 @@ void parallel_count_hc_pt1(pldat_t* pl) } } -ha_ug_index* build_unitig_index(ma_ug_t *ug, int k) +ha_ug_index* build_unitig_index(ma_ug_t *ug, int k, uint64_t up_occ, uint64_t low_occ, uint64_t thread_num) { ha_ug_index* idx = NULL; CALLOC(idx, 1); pldat_t pl; pl.h = idx; pl.is_cnt = 1; double index_time = yak_realtime(), beg_time; - init_ha_ug_index_opt(idx, ug, k, &pl); + init_ha_ug_index_opt(idx, ug, k, &pl, up_occ, low_occ, thread_num); beg_time = yak_realtime(); pl.is_cnt = 1; @@ -15250,7 +15253,7 @@ void verbose_het_stat(bubble_type *bub) fprintf(stderr, "[M::stat] # heterozygous bases: %lu; # homozygous bases: %lu\n", hetBase, homBase); } -void debug_output_disconnected_hits(ha_ug_index* idx, kvec_pe_hit *hits, hc_links *link, bubble_type *bub) +void debug_output_disconnected_hits(ha_ug_index* idx, kv_u_trans_t *ta, kvec_pe_hit *hits, hc_links *link, bubble_type *bub, int8_t *s) { uint32_t k, l, i, m, h_occ; uint64_t shif = 64 - idx->uID_bits, qn, tn, u_dis; @@ -15328,22 +15331,26 @@ void debug_output_disconnected_hits(ha_ug_index* idx, kvec_pe_hit *hits, hc_link u_dis = (t->e.a[k].dis ==(uint64_t)-1? (uint64_t)-1 : t->e.a[k].dis>>3); break; } + get_u_trans_spec(ta, qn, tn, &p, NULL); - fprintf(stderr, "s-utg%.6lul<--->d-utg%.6lul(occ: %u):(dis-%lu)\n", qn + 1, tn + 1, - ((uint32_t)-1) - k_trans.a[i].occ, u_dis); + fprintf(stderr, "s-utg%.6lul[hap-%d]<--->d-utg%.6lul[hap-%d](occ: %u, ", + qn + 1, s[qn], tn + 1, s[tn], ((uint32_t)-1) - k_trans.a[i].occ); + if(p) fprintf(stderr, "weight: %f", p->nw); + else fprintf(stderr, "weight: NA"); + fprintf(stderr, "):(dis-%lu)\n", u_dis); } - fprintf(stderr, "########hits########\n"); - char dir[2] = {'+', '-'}; - for (k = 0; k < hits->a.n; ++k) - { - fprintf(stderr, "%c\tutg%.6dl(len-%u)\t%lu\t%c\tutg%.6dl(len-%u)\t%lu\ti:%lu\n", - dir[hits->a.a[k].s>>63], (int)((hits->a.a[k].s<<1)>>shif)+1, - idx->ug->g->seq[((hits->a.a[k].s<<1)>>shif)].len, hits->a.a[k].s&idx->pos_mode, - dir[hits->a.a[k].e>>63], (int)((hits->a.a[k].e<<1)>>shif)+1, - idx->ug->g->seq[((hits->a.a[k].e<<1)>>shif)].len, hits->a.a[k].e&idx->pos_mode, - hits->a.a[k].id); - } + // fprintf(stderr, "########hits########\n"); + // char dir[2] = {'+', '-'}; + // for (k = 0; k < hits->a.n; ++k) + // { + // fprintf(stderr, "%c\tutg%.6dl(len-%u)\t%lu\t%c\tutg%.6dl(len-%u)\t%lu\ti:%lu\n", + // dir[hits->a.a[k].s>>63], (int)((hits->a.a[k].s<<1)>>shif)+1, + // idx->ug->g->seq[((hits->a.a[k].s<<1)>>shif)].len, hits->a.a[k].s&idx->pos_mode, + // dir[hits->a.a[k].e>>63], (int)((hits->a.a[k].e<<1)>>shif)+1, + // idx->ug->g->seq[((hits->a.a[k].e<<1)>>shif)].len, hits->a.a[k].e&idx->pos_mode, + // hits->a.a[k].id); + // } kv_destroy(k_trans); } @@ -15482,7 +15489,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) /*******************************for debug************************************/ // if(bub.round_id == bub.n_round - 1) // { - // debug_output_disconnected_hits(idx, &sl.hits, &link, &bub); + // debug_output_disconnected_hits(idx, &k_trans, &sl.hits, &link, &bub, s->s); // } /*******************************for debug************************************/ /** @@ -15540,7 +15547,7 @@ void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch) ug_index = NULL; int exist = (asm_opt.load_index_from_disk? load_hc_pt_index(&ug_index, asm_opt.output_file_name) : 0); - if(exist == 0) ug_index = build_unitig_index(ug, asm_opt.hic_mer_length); + if(exist == 0) ug_index = build_unitig_index(ug, asm_opt.hic_mer_length, asm_opt.hap_occ, 0, asm_opt.thread_num); if(exist == 0) write_hc_pt_index(ug_index, asm_opt.output_file_name); ug_index->ug = ug; ug_index->read_g = read_g; @@ -15551,6 +15558,64 @@ void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch) destory_hc_pt_index(ug_index); } +void init_ug_idx(ma_ug_t *ug, uint64_t k, uint64_t up_bound, uint64_t low_bound, uint64_t build_idx) +{ + ug_index = NULL; + if(build_idx) + { + ug_index = build_unitig_index(ug, k, up_bound, low_bound, asm_opt.thread_num); + } +} + +void des_ug_idx() +{ + destory_hc_pt_index(ug_index); +} + +uint64_t count_unique_k_mers(char *r, uint64_t len, uint64_t query, uint64_t target, uint64_t *all, uint64_t *found) +{ + if(!ug_index) return 0; + uint64_t i, j, l = 0, skip, *pos_list = NULL, cnt, uID, is_q, k_mer = ug_index->k; + uint64_t x[4], mask = (1ULL<>1)) & mask; + x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; + x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; + if (++l >= k_mer) + { + hash = hc_hash_long(x, &skip, k_mer); + if(skip == (uint64_t)-1) continue; + cnt = get_hc_pt1_count((ha_ug_index*)ug_index, hash, &pos_list); + if(cnt <= 0) continue; + + for (j = 0, is_q = 0; j < cnt; j++) + { + uID = (pos_list[j] << 1) >> (64 - ug_index->uID_bits); + if(query == uID) is_q = 1; + if(target == uID) break; + } + + if(is_q) + { + (*found)++; + if(j < cnt) (*all)++; + } + } + + } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart + } + + return 1; +} + + typedef struct{ //[uID_start, uID_end) uint64_t uID_start; @@ -15958,7 +16023,7 @@ void hic_benchmark(ma_ug_t *ug, asg_t* read_g) sprintf(output_file_name, "%s.bench", asm_opt.output_file_name); ug_index = NULL; int exist = load_hc_pt_index(&ug_index, output_file_name); - if(exist == 0) ug_index = build_unitig_index(ug, asm_opt.hic_mer_length); + if(exist == 0) ug_index = build_unitig_index(ug, asm_opt.hic_mer_length, asm_opt.hap_occ, 0, asm_opt.thread_num); if(exist == 0) write_hc_pt_index(ug_index, output_file_name); ug_index->ug = ug; ug_index->read_g = read_g; diff --git a/hic.h b/hic.h index fe1d166..dd3e71b 100644 --- a/hic.h +++ b/hic.h @@ -74,5 +74,8 @@ void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, u void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint32_t is_end); void set_b_utg_weight_flag(bubble_type* bub, buf_t* b, uint32_t v, uint8_t* vis_flag, uint32_t flag, uint32_t* occ); void debug_gfa_space(ma_ug_t* ug, hap_cov_t *cov); +void init_ug_idx(ma_ug_t *ug, uint64_t k, uint64_t up_bound, uint64_t low_bound, uint64_t build_idx); +void des_ug_idx(); +uint64_t count_unique_k_mers(char *r, uint64_t len, uint64_t query, uint64_t target, uint64_t *all, uint64_t *found); #endif