From 89fcb35a64588665f1b49bce7a93ff1074edd0fb Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Sun, 11 Oct 2020 14:32:33 -0400 Subject: [PATCH] for debugging --- Assembly.cpp | 30 +++++++---- CommandLines.h | 2 +- Overlaps.cpp | 34 ++++++++---- Process_Read.cpp | 3 ++ Process_Read.h | 1 + anchor.cpp | 10 ++-- htab.cpp | 136 ++++++++++++++++++++++++++++++++++++++++++----- htab.h | 11 ++-- sketch.cpp | 6 ++- 9 files changed, 192 insertions(+), 41 deletions(-) diff --git a/Assembly.cpp b/Assembly.cpp index f22c685..a87345c 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -12,7 +12,7 @@ #include "kthread.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); +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); void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list); All_reads R_INF; @@ -561,7 +561,7 @@ static void worker_ovec(void *data, long i, int tid) int fully_cov, abnormal; ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, - 0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region)); + 0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL); clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2); @@ -606,7 +606,6 @@ static void worker_ovec_related_reads(void *data, long i, int tid) uint64_t k, queryNameLen; for (k = 0; k < R_INF_FLAG.query_num; k++) { - queryNameLen = strlen(R_INF_FLAG.read_name[k]); if (queryNameLen != Get_NAME_LENGTH((R_INF),i)) continue; if (memcmp(R_INF_FLAG.read_name[k], Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) @@ -617,10 +616,10 @@ static void worker_ovec_related_reads(void *data, long i, int tid) if(k < R_INF_FLAG.query_num) { - int fully_cov, abnormal; + int fully_cov, abnormal, q_idx = k; ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, - 0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region)); + 0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), &(R_INF_FLAG.candidate_count[q_idx])); clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2); @@ -649,6 +648,18 @@ static void worker_ovec_related_reads(void *data, long i, int tid) fprintf(R_INF_FLAG.fp, "\n>%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i)); fprintf(R_INF_FLAG.fp, "%d-th round, len: %lu, hom_cov: %d, max_n_chain: %d\n", asm_opt.number_of_round, Get_READ_LENGTH(R_INF, i), asm_opt.hom_cov, asm_opt.max_n_chain); + + fprintf(R_INF_FLAG.fp, "***************************k-mer counts (%d)***************************\n", (int)(R_INF_FLAG.candidate_count[q_idx].a.n)); + + sort_kvec_t_u64_warp(&(R_INF_FLAG.candidate_count[q_idx]), 0); + for (k = 0; k < R_INF_FLAG.candidate_count[q_idx].a.n; k++) + { + fprintf(R_INF_FLAG.fp, "[%lu] Count(%u): %lu, filtered: %lu\n", k, + (uint32_t)R_INF_FLAG.candidate_count[q_idx].a.a[k], R_INF_FLAG.candidate_count[q_idx].a.a[k]>>33, + (R_INF_FLAG.candidate_count[q_idx].a.a[k]>>32)&(uint64_t)1); + } + + fprintf(R_INF_FLAG.fp, "***************************forward ovlp***************************\n"); for (k = 0; k < b->olist.length; k++) { @@ -692,6 +703,7 @@ static void worker_ovec_related_reads(void *data, long i, int tid) // } fprintf(R_INF_FLAG.fp, "R_INF.trio_flag[%ld]: %u\n", i, R_INF.trio_flag[i]); + pthread_mutex_unlock(&R_INF_FLAG.OutputMutex); } @@ -867,7 +879,7 @@ void ha_overlap_and_correct(int round) else kt_for(asm_opt.thread_num, worker_ovec, b, R_INF.total_reads); - if (r_out) write_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name); + if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name); ha_pt_destroy(ha_idx); ha_idx = NULL; @@ -1298,7 +1310,7 @@ static void worker_ov_final(void *data, long i, int tid) //get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0); ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 0.001, - asm_opt.max_n_chain, 0, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region)); + asm_opt.max_n_chain, 0, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL); overlap_region_sort_y_id(b->olist.list, b->olist.length); ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length); @@ -1364,7 +1376,7 @@ static void worker_ov_final_high_het(void *data, long i, int tid) ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid]; ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, HIGH_HET_ERROR_RATE, - asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region)); + asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL); overlap_region_sort_y_id(b->olist.list, b->olist.length); ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length); @@ -1631,7 +1643,7 @@ int ha_assemble(void) } if (!ovlp_loaded) { ha_flt_tab = ha_idx = NULL; - if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name); + if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name), load_ct_index(&ha_ct_table, asm_opt.output_file_name); // construct hash table for high occurrence k-mers if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL) diff --git a/CommandLines.h b/CommandLines.h index c431096..10a9f4b 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.12-r304" +#define HA_VERSION "0.12-r305" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index 20a3aca..ca1c72d 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -2107,20 +2107,30 @@ static inline int asg_is_single_edge(const asg_t *g, uint32_t v, uint32_t start_ return nv; } -void debug_info_of_specfic_node(const char* name, asg_t *g, char* command) +void debug_info_of_specfic_node(char* name, asg_t *g, R_to_U* ruIndex, char* command) { fprintf(stderr, "\n\n\n"); - uint32_t v, n_vtx = g->n_seq * 2, queryLen = strlen(name), flag = 0; + uint32_t v, n_vtx = g->n_seq * 2, queryLen = strlen(name), flag = 0, contain_rId, is_Unitig; for (v = 0; v < n_vtx; ++v) { if(queryLen == Get_NAME_LENGTH(R_INF, (v>>1)) && memcmp(name, Get_NAME(R_INF, (v>>1)), Get_NAME_LENGTH(R_INF, (v>>1))) == 0) { if(flag == 0) fprintf(stderr, "\nafter %s\n", command); - fprintf(stderr, "****************graph ref_read: %.*s, dir: %u****************\n", - (int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), v&1); + fprintf(stderr, "****************graph ref_read: %.*s, id: %u, dir: %u****************\n", + (int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), v>>1, v&1); if(g->seq[v>>1].del) { - fprintf(stderr, "read has already been deleted.\n"); + get_R_to_U(ruIndex, (v>>1), &contain_rId, &is_Unitig); + if(contain_rId != (uint32_t)-1 && is_Unitig != 1) + { + fprintf(stderr, "read is deleted as a contained read by: %.*s\n contain_rId: %u, del: %u\n", + (int)Get_NAME_LENGTH(R_INF, contain_rId), Get_NAME(R_INF, contain_rId), contain_rId, g->seq[contain_rId].del); + } + else + { + fprintf(stderr, "read has already been deleted.\n"); + } + return; } @@ -9757,12 +9767,13 @@ ma_hit_t_alloc* reverse_sources, int id, const char* command) for (j = 0; j < reverse_sources[i].length; j++) { tn = Get_tn(reverse_sources[i].buffer[j]); - fprintf(stderr, "target: %.*s, qs: %u, qe: %u, ts: %u, te: %u, del: %u\n", + fprintf(stderr, "target: %.*s, qs: %u, qe: %u, ts: %u, te: %u, rev: %u, del: %u\n", (int)Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn), Get_qs(reverse_sources[i].buffer[j]), Get_qe(reverse_sources[i].buffer[j]), Get_ts(reverse_sources[i].buffer[j]), Get_te(reverse_sources[i].buffer[j]), + reverse_sources[i].buffer[j].rev, (uint32_t)sources[i].buffer[j].del); } @@ -26919,6 +26930,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) ///it seems we do not need ma_hit_cut & ma_hit_flt ma_hit_sub(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut); detect_chimeric_reads(sources, n_read, readLen, coverage_cut, asm_opt.max_ov_diff_final * 2.0); + 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); @@ -26926,8 +26938,10 @@ ma_sub_t **coverage_cut_ptr, int debug_g) ///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); - sg = ma_sg_gen(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length); + + ///debug_info_of_specfic_node((char*)"m64062_190804_172951/130483063/ccs", sg, ruIndex, (char*)"sbsbsb"); + asg_arc_del_trans(sg, gap_fuzz); asm_opt.coverage = get_coverage(sources, coverage_cut, n_read); @@ -27154,9 +27168,11 @@ long long bubble_dist, int read_graph, int write) &R_INF, output_file_name); } + debug_info_of_specfic_read("m64062_190803_042216/120916077/ccs", sources, reverse_sources, -1, "beg"); + debug_info_of_specfic_read("m64062_190806_063919/71436446/ccs", sources, reverse_sources, -1, "beg"); + debug_info_of_specfic_read("m64062_190803_042216/82117654/ccs", sources, reverse_sources, -1, "beg"); + - ///debug_info_of_specfic_read("m64062_190803_042216/15205346/ccs", sources, reverse_sources, -1, "beg"); - if (!(asm_opt.flag & HA_F_BAN_ASSEMBLY)) { try_rescue_overlaps(sources, reverse_sources, n_read, 4); diff --git a/Process_Read.cpp b/Process_Read.cpp index 316a5d6..cb23a98 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -690,6 +690,7 @@ void init_Debug_reads(Debug_reads* x, const char* file) x->query_num++; } x->read_name = (char**)malloc(sizeof(char*)*x->query_num); + x->candidate_count = (kvec_t_u64_warp*)malloc(sizeof(kvec_t_u64_warp)*x->query_num); fseek(x->fp, 0, SEEK_SET); i = 0; @@ -699,6 +700,7 @@ void init_Debug_reads(Debug_reads* x, const char* file) x->read_name[i] = (char*)malloc(sizeof(char)*(nameLen+1)); memcpy(x->read_name[i], Name_Buffer, sizeof(char)*nameLen); x->read_name[i][nameLen] = '\0'; + kv_init(x->candidate_count[i].a); i++; } @@ -716,6 +718,7 @@ void destory_Debug_reads(Debug_reads* x) for (i = 0; i < x->query_num; i++) { free(x->read_name[i]); + kv_destroy(x->candidate_count[i].a); } free(x->read_name); diff --git a/Process_Read.h b/Process_Read.h index d259bd4..01cd42c 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -148,6 +148,7 @@ typedef struct { char** read_name; uint64_t query_num; + kvec_t_u64_warp* candidate_count; FILE* fp; pthread_mutex_t OutputMutex; } Debug_reads; diff --git a/anchor.cpp b/anchor.cpp index 7e0bd55..50a644f 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -58,7 +58,7 @@ int ha_ov_type(const overlap_region *r, uint32_t len) } void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, 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, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar) +kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct) { uint32_t i, rlen; uint64_t k, l; @@ -73,7 +73,7 @@ kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f rlen = Get_READ_LENGTH(R_INF, rid); // read length // get the list of anchors - ha_sketch_query(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag); + ha_sketch_query(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag, dbg_ct); // minimizer of queried read if (ab->mz.m > ab->old_mz_m) { ab->old_mz_m = ab->mz.m; @@ -205,14 +205,14 @@ void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf) 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) +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) { extern void *ha_flt_tab; extern ha_pt_t *ha_idx; extern void *ha_flt_tab_hp; extern ha_pt_t *ha_idx_hp; - ha_get_new_candidates(ab, rid, ucr, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar); + ha_get_new_candidates(ab, rid, ucr, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar, dbg_ct); if(ha_idx_hp) { @@ -242,7 +242,7 @@ int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* overlap_list->length = k; - ha_get_new_candidates(ab, rid, ucr, overlap_list_hp, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab_hp, ha_idx_hp, f_cigar); + ha_get_new_candidates(ab, rid, ucr, overlap_list_hp, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab_hp, ha_idx_hp, f_cigar, dbg_ct); if(overlap_list->length + overlap_list_hp->length > overlap_list->size) { diff --git a/htab.cpp b/htab.cpp index 0e1ec4f..dcbf9c9 100644 --- a/htab.cpp +++ b/htab.cpp @@ -37,6 +37,7 @@ void *ha_flt_tab; ha_pt_t *ha_idx; void *ha_flt_tab_hp; ha_pt_t *ha_idx_hp; +void *ha_ct_table; /*************************** * Yak specific parameters * @@ -132,8 +133,8 @@ typedef struct { ha_ct1_t *h; } ha_ct_t; -///for 0-th counting, k = 51, pre = 12, n_hash = 4, n_shift = 0 -///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 37 +///for 0-th counting, k = 51, pre = 12, n_hash = 4, n_shift = 37 +///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0 static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift) { ha_ct_t *h; @@ -146,8 +147,7 @@ static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift) ///it seems there is a large hash table h, consisting 4096 small hash tables for (i = 0; i < 1<pre; ++i) h->h[i].h = yak_ct_init(); - ///for 0-th counting, don't enter here - ///seems used for minimzer + ///for 0-th counting, enter here; used for bloom filter if (n_hash > 0 && n_shift > h->pre) { h->n_hash = n_hash, h->n_shift = n_shift; for (i = 0; i < 1<pre; ++i) @@ -749,10 +749,10 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt } else if (c0) { pl.ct = c0, pl.create_new = !!(flag&HAF_CREATE_NEW); assert(c0->k == opt->k && c0->pre == opt->pre); - } else {///for 0-th counting and 1-th counting, go into here + } else {///for ft-th counting and 1-th counting, go into here pl.create_new = 1; // alware create new elements if the count table is empty - ///for 0-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0 - ///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 37 + ///for 0-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 37 + ///for 1-th counting, opt.k = 51, opt->pre = 12, opt->bf_n_hash = 4, opt.bf_shift = 0 ///building a large hash table consisting of 4096 small hash tables pl.ct = ha_ct_init(opt->k, opt->pre, opt->bf_n_hash, opt->bf_shift); } @@ -785,11 +785,10 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi opt.k = asm_opt->k_mer_length; ///always 0 opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC); - ///for 0-th counting, shoud be 1 - ///for 1-th counting, shoud be 51 + ///for ft-counting, shoud be 1 opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win; - ///for 0-th counting, shoud be 0 - ///for 1-th counting, shoud be 37 + ///for ft-counting, shoud be 37 + ///for ha_pt_gen, shoud be 0 opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift; opt.n_thread = asm_opt->thread_num; ///asm_opt->num_reads is the number of fastq files @@ -839,6 +838,32 @@ void ha_ft_destroy(void *h) if (h) yak_ft_destroy((yak_ft_t*)h); } + +void debug_ct_index(void* q_ct_idx, void* r_ct_idx) +{ + ha_ct_t* ct_idx = (ha_ct_t*)q_ct_idx; + yak_ct_t *g = NULL; + uint64_t i; + khint_t k; + for (i = 0; (int)i < 1<pre; i++) + { + g = ct_idx->h[i].h; + for (k = 0; k < kh_end(g); ++k) + { + if (kh_exist(g, k)) + { + int c = kh_key(g, k) & YAK_MAX_COUNT; + uint64_t hash = ((kh_key(g, k) >> ct_idx->pre)<pre) | i; + int q = query_ct_index(r_ct_idx, hash); + if(q!=c) + { + fprintf(stderr, "ERROR:c: %d, q: %d\n", c, q); + } + } + } + } +} + /************************* * High-level interfaces * *************************/ @@ -851,6 +876,15 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i 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, NULL, NULL, rs); + if((asm_opt->flag & HA_F_VERBOSE_GFA)) + { + write_ct_index((void*)h, asm_opt->output_file_name); + // load_ct_index(&ha_ct_table, asm_opt->output_file_name); + // debug_ct_index((void*)h, ha_ct_table); + // debug_ct_index(ha_ct_table, (void*)h); + // ha_ct_destroy((ha_ct_t *)ha_ct_table); + } + if(!(ex_flag & HAF_SKIP_READ)) { ha_ct_hist(h, cnt, asm_opt->thread_num); @@ -917,9 +951,85 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f return pt; } +int query_ct_index(void* ct_idx, uint64_t hash) +{ + ha_ct1_t *g = &(((ha_ct_t*)ct_idx)->h[hash & ((1ULL<<((ha_ct_t*)ct_idx)->pre) - 1)]); + khint_t k; + k = yak_ct_get(g->h, hash); + if (k == kh_end(g->h)) return 0; + return kh_key(g->h, k)&YAK_MAX_COUNT; +} -int write_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name) +int write_ct_index(void *i_ct_idx, char* file_name) +{ + char* gfa_name = (char*)malloc(strlen(file_name)+25); + sprintf(gfa_name, "%s.ct_flt", file_name); + FILE* fp = fopen(gfa_name, "w"); + if (!fp) { + free(gfa_name); + return 0; + } + ha_ct_t* ct_idx = (ha_ct_t*)i_ct_idx; + int i; + ha_ct1_t *g; + fwrite(&ct_idx->k, sizeof(ct_idx->k), 1, fp); + fwrite(&ct_idx->pre, sizeof(ct_idx->pre), 1, fp); + fwrite(&ct_idx->n_hash, sizeof(ct_idx->n_hash), 1, fp); + fwrite(&ct_idx->n_shift, sizeof(ct_idx->n_shift), 1, fp); + fwrite(&ct_idx->tot, sizeof(ct_idx->tot), 1, fp); + for (i = 0; i < 1<pre; i++) + { + g = &(ct_idx->h[i]); + yak_ct_save(g->h, fp); + } + + + fprintf(stderr, "[M::%s] Index has been written.\n", __func__); + free(gfa_name); + fclose(fp); + return 1; +} + +int load_ct_index(void **i_ct_idx, char* file_name) +{ + char* gfa_name = (char*)malloc(strlen(file_name)+25); + sprintf(gfa_name, "%s.ct_flt", file_name); + FILE* fp = fopen(gfa_name, "r"); + if (!fp) { + free(gfa_name); + return 0; + } + ha_ct_t** ct_idx = (ha_ct_t**)i_ct_idx; + double index_time = 0; + int i; + ha_ct_t *h = 0; + ha_ct1_t *g; + CALLOC(h, 1); + + fread(&h->k, sizeof(h->k), 1, fp); + fread(&h->pre, sizeof(h->pre), 1, fp); + fread(&h->n_hash, sizeof(h->n_hash), 1, fp); + fread(&h->n_shift, sizeof(h->n_shift), 1, fp); + fread(&h->tot, sizeof(h->tot), 1, fp); + CALLOC(h->h, 1<pre); + + + index_time = yak_realtime(); + for (i = 0; i < 1<pre; ++i) + { + g = &(h->h[i]); + yak_ct_load(&(g->h), fp); + } + + (*ct_idx) = h; + fprintf(stderr, "[M::%s::%.3f] ==> Loaded count table\n", __func__, yak_realtime() - index_time); + fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__); + free(gfa_name); + return 1; +} + +int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name) { char* gfa_name = (char*)malloc(strlen(file_name)+25); sprintf(gfa_name, "%s.pt_flt", file_name); @@ -970,7 +1080,7 @@ int write_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt return 1; } -int load_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name) +int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name) { char* gfa_name = (char*)malloc(strlen(file_name)+25); sprintf(gfa_name, "%s.pt_flt", file_name); diff --git a/htab.h b/htab.h index 42696f6..c6feba6 100644 --- a/htab.h +++ b/htab.h @@ -29,6 +29,8 @@ extern void *ha_flt_tab; extern ha_pt_t *ha_idx; extern void *ha_flt_tab_hp; extern ha_pt_t *ha_idx_hp; +extern void *ha_ct_table; + void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode); int ha_ft_isflt(const void *hh, uint64_t y); @@ -38,8 +40,11 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f void ha_pt_destroy(ha_pt_t *h); const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n); -int write_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name); -int load_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name); +int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name); +int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name); +int write_ct_index(void *ct_idx, char* file_name); +int load_ct_index(void **ct_idx, char* file_name); +int query_ct_index(void* ct_idx, uint64_t hash); ha_abuf_t *ha_abuf_init(void); void ha_abuf_destroy(ha_abuf_t *ab); @@ -55,7 +60,7 @@ double yak_cpu_usage(void); void ha_triobin(const hifiasm_opt_t *opt); void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf); -void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag); +void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct); int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het); static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function diff --git a/sketch.cpp b/sketch.cpp index 046bdf0..60694d6 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -139,12 +139,15 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, -kvec_t_u8_warp* k_flag) +kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct) { ///in default, w = 51, k = 51, is_hpc = 1 /** uint64_t x; uint64_t rid:28, pos:27, rev:1, span:8; **/ + extern void *ha_ct_table; + if(dbg_ct != NULL) dbg_ct->a.n = 0; + static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 }; uint64_t shift1 = k - 1, mask = (1ULL<a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i)); ///if (hf == 0 || ha_ft_isflt(hf, y) == 0) if(filtered == 0) info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;