From 03cf0d475c4174202229b5a1524906da53f828e0 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Sun, 23 Aug 2020 08:46:33 -0400 Subject: [PATCH] save bfore multiple overlaps --- Assembly.cpp | 234 +++++++++++++++++++++++++++++++------- Assembly.h | 5 + CommandLines.cpp | 2 +- CommandLines.h | 2 +- Correct.cpp | 104 ++++++++++++++++- Correct.h | 1 + Hash_Table.cpp | 21 ++++ Overlaps.cpp | 184 +++++++++++++++++++++--------- Process_Read.cpp | 54 +++++++++ Process_Read.h | 14 ++- anchor.cpp | 148 ++++++++++++++++++++++-- htab.cpp | 286 +++++++++++++++++++++++++++++++++++++++++++---- htab.h | 15 ++- khashl.h | 32 ++++++ sketch.cpp | 149 +++++++++++++++++++++++- 15 files changed, 1120 insertions(+), 131 deletions(-) diff --git a/Assembly.cpp b/Assembly.cpp index 2286562..8edfe4a 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -11,10 +11,12 @@ #include "htab.h" #include "kthread.h" -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); +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, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf); void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list); All_reads R_INF; +Debug_reads R_INF_FLAG; void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length, char* new_read, int* new_length) { @@ -442,6 +444,7 @@ typedef struct { UC_Read self_read, ovlp_read; Candidates_list clist; overlap_region_alloc olist; + overlap_region_alloc olist_hp; ha_abuf_t *ab; // error correction related buffers int64_t num_read_base, num_correct_base, num_recorrect_base; @@ -451,6 +454,9 @@ typedef struct { Correct_dumy correct; haplotype_evdience_alloc hap; Round2_alignment round2; + kvec_t_u32_warp b_buf; + kvec_t_u64_warp r_buf; + kvec_t_u8_warp k_flag; } ha_ovec_buf_t; ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov) @@ -462,6 +468,10 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov) init_UC_Read(&b->ovlp_read); init_Candidates_list(&b->clist); init_overlap_region_alloc(&b->olist); + init_overlap_region_alloc(&b->olist_hp); + kv_init(b->b_buf.a); + kv_init(b->r_buf.a); + kv_init(b->k_flag.a); b->ab = ha_abuf_init(); if (!b->is_final) { init_Cigar_record(&b->cigar1); @@ -480,7 +490,11 @@ void ha_ovec_destroy(ha_ovec_buf_t *b) destory_UC_Read(&b->ovlp_read); destory_Candidates_list(&b->clist); destory_overlap_region_alloc(&b->olist); + destory_overlap_region_alloc(&b->olist_hp); ha_abuf_destroy(b->ab); + kv_destroy(b->b_buf.a); + kv_destroy(b->r_buf.a); + kv_destroy(b->k_flag.a); if (!b->is_final) { destory_Cigar_record(&b->cigar1); destory_Graph(&b->POA_Graph); @@ -510,6 +524,7 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b) { int64_t i, mem = 0, mem_clist, mem_olist; mem_clist = b->clist.size * sizeof(k_mer_hit) + b->clist.chainDP.size * 7 * 4; + mem_olist = b->olist.size * sizeof(overlap_region); for (i = 0; i < (int64_t)b->olist.size; ++i) { const overlap_region *r = &b->olist.list[i]; @@ -517,6 +532,14 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b) mem_olist += r->f_cigar.size * 8; mem_olist += r->boundary_cigars.size * sizeof(window_list); } + mem_olist += b->olist_hp.size * sizeof(overlap_region); + for (i = 0; i < (int64_t)b->olist_hp.size; ++i) { + const overlap_region *r = &b->olist_hp.list[i]; + mem_olist += r->w_list_size * sizeof(window_list); + mem_olist += r->f_cigar.size * 8; + mem_olist += r->boundary_cigars.size * sizeof(window_list); + } + mem = ha_abuf_mem(b->ab) + mem_clist + mem_olist; if (!b->is_final) { mem += sizeof(Cigar_record) + b->cigar1.lost_base_size + b->cigar1.size * 4; @@ -534,8 +557,8 @@ static void worker_ovec(void *data, long i, int tid) ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid]; int fully_cov, abnormal; - ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1); - ///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 1); + 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), &(R_INF.paf[i]), &(R_INF.reverse_paf[i])); clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2); @@ -557,7 +580,13 @@ static void worker_ovec(void *data, long i, int tid) } R_INF.paf[i].is_abnormal = abnormal; - if (b->save_ov) { + R_INF.trio_flag[i] = AMBIGU; + if(ha_idx_hp == NULL) + { + R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &b->b_buf, &b->r_buf, &(b->k_flag), RESEED_HP_RATE, NULL); + } + + if (R_INF.trio_flag[i] != AMBIGU || b->save_ov) { int is_rev = (asm_opt.number_of_round % 2 == 0); push_overlaps(&(R_INF.paf[i]), &b->olist, 1, &R_INF, is_rev); push_overlaps(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, is_rev); @@ -568,25 +597,97 @@ static void worker_ovec(void *data, long i, int tid) static void worker_ovec_related_reads(void *data, long i, int tid) { ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid]; - int required_read_name_length = strlen(asm_opt.required_read_name); - uint64_t k; - if (required_read_name_length == (int)Get_NAME_LENGTH((R_INF),i) - && - memcmp(asm_opt.required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) - { - ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1); - ///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 1); - fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i)); - recover_UC_Read(&b->self_read, &R_INF, i); - fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq); + 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) + { + break; + } + } - for (k = 0; k < b->olist.length; k++) { - fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id)); - recover_UC_Read(&b->self_read, &R_INF, b->olist.list[k].y_id); - fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq); + if(k < R_INF_FLAG.query_num) + { + 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), &(R_INF.paf[i]), &(R_INF.reverse_paf[i])); + + clear_Cigar_record(&b->cigar1); + clear_Round2_alignment(&b->round2); + + correct_overlap(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon, + &b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal); + + b->num_read_base += b->self_read.length; + b->num_correct_base += b->correct.corrected_base; + b->num_recorrect_base += b->round2.dumy.corrected_base; + + push_cigar(R_INF.cigars, i, &b->cigar1); + push_cigar(R_INF.second_round_cigar, i, &b->round2.cigar); + + R_INF.paf[i].is_fully_corrected = 0; + if (fully_cov) { + if (get_cigar_errors(&b->cigar1) == 0 && get_cigar_errors(&b->round2.cigar) == 0) + R_INF.paf[i].is_fully_corrected = 1; + } + R_INF.paf[i].is_abnormal = abnormal; + + + + pthread_mutex_lock(&R_INF_FLAG.OutputMutex); + + 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, "***************************forward ovlp***************************\n"); + for (k = 0; k < b->olist.length; k++) + { + if(b->olist.list[k].is_match != 1) continue; + fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id)); + fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n", + b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e, + b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel, + Get_READ_LENGTH(R_INF, b->olist.list[k].y_id)); } - } + + fprintf(R_INF_FLAG.fp, "***************************reverse ovlp***************************\n"); + for (k = 0; k < b->olist.length; k++) + { + if(b->olist.list[k].is_match != 2) continue; + fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id)); + fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n", + b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e, + b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel, + Get_READ_LENGTH(R_INF, b->olist.list[k].y_id)); + } + + fprintf(R_INF_FLAG.fp, "***************************unmatched ovlp***************************\n"); + for (k = 0; k < b->olist.length; k++) + { + if(b->olist.list[k].is_match == 1) continue; + if(b->olist.list[k].is_match == 2) continue; + fprintf(R_INF_FLAG.fp, "%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id)); + fprintf(R_INF_FLAG.fp, "qs: %u, qe: %u, ts: %u, te: %u, rev: %u, strong: %u, no_l_indel: %u, len: %lu\n", + b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e, + b->olist.list[k].y_pos_strand, b->olist.list[k].strong, b->olist.list[k].without_large_indel, + Get_READ_LENGTH(R_INF, b->olist.list[k].y_id)); + } + + R_INF.trio_flag[i] = AMBIGU; + if(ha_idx_hp == NULL) + { + R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &b->b_buf, &b->r_buf, &(b->k_flag), RESEED_HP_RATE, R_INF_FLAG.fp); + } + + 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); + } } @@ -693,25 +794,75 @@ void Output_corrected_reads() fclose(output_file); } +void debug_print_pob_regions() +{ + uint64_t i, total = 0; + for (i = 0; i < R_INF.total_reads; i++) + { + if(R_INF.trio_flag[i]!=AMBIGU) + { + total++; + fprintf(stderr, "(%lu) %.*s\n", i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + } + } + fprintf(stderr, "total hp reads: %lu, R_INF.total_reads: %lu\n", total, R_INF.total_reads); + exit(1); +} + +void rescue_hp_reads(ha_ovec_buf_t **b) +{ + int hom_cov, het_cov; + ha_flt_tab_hp = ha_idx_hp = NULL; + if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) { + ha_flt_tab_hp = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 1); + } + ha_idx_hp = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 1, &R_INF, &hom_cov, &het_cov); + + + if (asm_opt.required_read_name) + kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads); + else + kt_for(asm_opt.thread_num, worker_ovec, b, R_INF.total_reads); + + + + + ha_ft_destroy(ha_flt_tab_hp); ha_flt_tab_hp = NULL; + ha_pt_destroy(ha_idx_hp); ha_idx_hp = NULL; +} + + + + void ha_overlap_and_correct(int round) { - int i, hom_cov, het_cov; + int i, hom_cov, het_cov, r_out = 0; ha_ovec_buf_t **b; ha_ecsave_buf_t *e; + ha_flt_tab_hp = ha_idx_hp = NULL; + if((ha_idx == NULL)&&(asm_opt.flag & HA_F_VERBOSE_GFA)&&(round == asm_opt.number_of_round - 1)) + { + r_out = 1; + } + + if(asm_opt.required_read_name) init_Debug_reads(&R_INF_FLAG, asm_opt.required_read_name); // for debugging only // overlap and correct reads CALLOC(b, asm_opt.thread_num); for (i = 0; i < asm_opt.thread_num; ++i) b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1)); - ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF, &hom_cov, &het_cov); // build the index + if(ha_idx) hom_cov = asm_opt.hom_cov; + if(ha_idx == NULL) ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, 0, &R_INF, &hom_cov, &het_cov); // build the index if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated ha_opt_update_cov(&asm_opt, hom_cov); if (asm_opt.required_read_name) kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads); 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); ha_pt_destroy(ha_idx); - ha_idx = 0; + ha_idx = NULL; // collect statistics for (i = 0; i < asm_opt.thread_num; ++i) { @@ -723,8 +874,8 @@ void ha_overlap_and_correct(int round) } free(b); - if (asm_opt.required_read_name) exit(0); // for debugging only - + if (asm_opt.required_read_name) destory_Debug_reads(&R_INF_FLAG), exit(0); // for debugging only + // save corrected reads to R_INF CALLOC(e, asm_opt.thread_num); for (i = 0; i < asm_opt.thread_num; ++i) { @@ -740,6 +891,7 @@ void ha_overlap_and_correct(int round) free(e[i].second_round_read); } free(e); + debug_print_pob_regions(); } @@ -1138,15 +1290,8 @@ static void worker_ov_final(void *data, long i, int tid) ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid]; //get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0); - ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.001, asm_opt.max_n_chain, 0); - ///ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.08, asm_opt.max_n_chain, 0); - - /** - correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon, - &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1, - ¤t_cigar, &hap, &second_round, 0, 0); - push_final_overlaps(&(R_INF.paf[i]), &overlap_list); - **/ + 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), &(R_INF.paf[i]), &(R_INF.reverse_paf[i])); 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); @@ -1211,7 +1356,8 @@ 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_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, HIGH_HET_ERROR_RATE, asm_opt.max_n_chain, 1); + 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), &(R_INF.paf[i]), &(R_INF.reverse_paf[i])); 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); @@ -1410,7 +1556,7 @@ void hap_recalculate_peaks(char* output_file_name) int hom_cov, het_cov; // construct hash table for high occurrence k-mers if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) { - ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov); + ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0); ha_opt_update_cov(&asm_opt, hom_cov); } free(R_INF.read_length); @@ -1418,7 +1564,7 @@ void hap_recalculate_peaks(char* output_file_name) load_All_reads(&R_INF, gfa_name); - ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index + ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index asm_opt.hom_cov = hom_cov; asm_opt.het_cov = het_cov; ha_pt_destroy(ha_idx); @@ -1435,10 +1581,12 @@ void ha_overlap_final(void) { int i, hom_cov, het_cov; ha_ovec_buf_t **b; + ha_flt_tab_hp = ha_idx_hp = NULL; + CALLOC(b, asm_opt.thread_num); for (i = 0; i < asm_opt.thread_num; ++i) b[i] = ha_ovec_init(asm_opt.flag & HA_F_HIGH_HET, 1);///b[i] = ha_ovec_init(1, 1); - ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov, &het_cov); // build the index + ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index if(asm_opt.flag & HA_F_HIGH_HET) { kt_for(asm_opt.thread_num, worker_ov_final_high_het, b, R_INF.total_reads); @@ -1475,14 +1623,18 @@ int ha_assemble(void) if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2; } 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); + // construct hash table for high occurrence k-mers - if (!(asm_opt.flag & HA_F_NO_KMER_FLT)) { - ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov); + if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL) + { + ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0); ha_opt_update_cov(&asm_opt, hom_cov); } // error correction assert(asm_opt.number_of_round > 0); - for (r = 0; r < asm_opt.number_of_round; ++r) { + for (r = ha_idx?asm_opt.number_of_round-1:0; r < asm_opt.number_of_round; ++r) { ha_opt_reset_to_round(&asm_opt, r); // this update asm_opt.roundID and a few other fields ha_overlap_and_correct(r); fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(), diff --git a/Assembly.h b/Assembly.h index 339d5e1..552004a 100644 --- a/Assembly.h +++ b/Assembly.h @@ -8,6 +8,11 @@ #define Get_Cigar_Type(RECORD) (RECORD&3) #define Get_Cigar_Length(RECORD) (RECORD>>2) +#define RESEED_DP 4 +#define RESEED_PEAK_RATE 0.15 +#define RESEED_LEN 2000 +#define RESEED_HP_RATE 0.9 + int ha_assemble(void); #endif diff --git a/CommandLines.cpp b/CommandLines.cpp index 271b712..221fb0e 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -324,7 +324,7 @@ int check_option(hifiasm_opt_t* asm_opt) if(asm_opt->fn_bin_yak[1] != NULL && check_file(asm_opt->fn_bin_yak[1], "YAK2") == 0) return 0; if(asm_opt->fn_bin_list[0] != NULL && check_file(asm_opt->fn_bin_list[0], "LIST1") == 0) return 0; if(asm_opt->fn_bin_list[1] != NULL && check_file(asm_opt->fn_bin_list[1], "LIST2") == 0) return 0; - + if(asm_opt->required_read_name != NULL && check_file(asm_opt->required_read_name, "b") == 0) return 0; // fprintf(stderr, "input file num: %d\n", asm_opt->num_reads); // fprintf(stderr, "output file: %s\n", asm_opt->output_file_name); // fprintf(stderr, "number of threads: %d\n", asm_opt->thread_num); diff --git a/CommandLines.h b/CommandLines.h index 95759bd..2d13ba3 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.11-r302" +#define HA_VERSION "0.11-r303" #define VERBOSE 0 diff --git a/Correct.cpp b/Correct.cpp index d18ffa7..0389134 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -7,8 +7,10 @@ #include "Assembly.h" #include "CommandLines.h" #include "ksw2.h" +#include "ksort.h" - +#define generic_key(x) (x) +KRADIX_SORT_INIT(b32, uint32_t, generic_key, 4) void clear_Round2_alignment(Round2_alignment* h) @@ -7202,7 +7204,107 @@ void partition_overlaps_advance(overlap_region_alloc* overlap_list, All_reads* R } +void collect_no_cov_regions(overlap_region_alloc* overlap_list, All_reads* R_INF, +kvec_t_u32_warp* b, kvec_t_u64_warp* r, int min_dp, int min_len) +{ + b->a.n = r->a.n = 0; + ///if(overlap_list->length == 0) return; + long long i = 0, xLen = Get_READ_LENGTH((*R_INF), overlap_list->list[0].x_id); + uint32_t qs, qe; + uint64_t tmp; + int dp, old_dp, s_start = 0, s_end = 0; + ///at least 1 + if(min_len < 1) min_len = 1; + + for (i = 0; i < (long long)overlap_list->length; i++) + { + if (overlap_list->list[i].is_match != 1 && overlap_list->list[i].is_match != 2) continue; + + qs = overlap_list->list[i].x_pos_s; + qe = overlap_list->list[i].x_pos_e + 1; + kv_push(uint32_t, b->a, qs<<1); + kv_push(uint32_t, b->a, qe<<1|1); + } + + + ///we can identify the qs and qe by the 0-th bit + radix_sort_b32(b->a.a, b->a.a + b->a.n); + + for (i = 0, dp = 0; i < (long long)b->a.n; ++i) + { + old_dp = dp; + //if a[j] is qe + if (b->a.a[i]&1) --dp; + else ++dp; + /** + min_dp is the coverage drop threshold + there are two cases: + 1. old_dp = dp + 1 (b.a[j] is qe); 2. old_dp = dp - 1 (b.a[j] is qs); + **/ + if (old_dp < min_dp && dp >= min_dp) ///old_dp < dp, b.a[j] is qs + { + ///case 2, a[j] is qs + s_end = b->a.a[i]>>1; + ///at least 1 + if(s_end-s_start >= min_len) + { + tmp = s_start; tmp = tmp << 32; tmp = tmp | (uint64_t)(s_end-1); + kv_push(uint64_t, r->a, tmp); + } + } + else if (old_dp >= min_dp && dp < min_dp) ///old_dp > min_dp, b.a[j] is qe + { + s_start = b->a.a[i]>>1; + } + } + + + if(s_start < xLen && xLen-s_start >= min_len) + { + s_end = xLen; + tmp = s_start; tmp = tmp << 32; tmp = tmp | (uint64_t)(s_end-1); + kv_push(uint64_t, r->a, tmp); + } +} + +int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u32_warp* b, kvec_t_u64_warp* r, kvec_t_u8_warp* k_flag, float hp_rate, FILE* fp) +{ + int i, k, qs, qe, ava_k_mer = 0, hp_k_mer = 0; + int min_dp = RESEED_DP; + if(asm_opt.hom_cov > 0) min_dp = asm_opt.hom_cov * RESEED_PEAK_RATE; + if(min_dp > RESEED_DP) min_dp = RESEED_DP; + collect_no_cov_regions(olist, R_INF, b, r, min_dp, RESEED_LEN); + + for (i = 0; i < (int)r->a.n; i++) + { + ///[qs, qe] + qs = r->a.a[i]>>32; + qe = (r->a.a[i]<<32)>>32; + + for (k = qs; k <= qe; k++) + { + if(k_flag->a.a[k] > 1) ava_k_mer++; + if(k_flag->a.a[k] > 2) hp_k_mer++; + } + + if(fp) fprintf(fp, "qs: %d, qe: %d, ava_k_mer: %d, hp_k_mer: %d\n", qs, qe, ava_k_mer, hp_k_mer); + } + + if(fp) fprintf(fp, "ava_k_mer: %d, hp_k_mer: %d, hp_rate: %f\n", ava_k_mer, hp_k_mer, hp_rate); + + if(fp) + { + for (k = 0; k < (int)k_flag->a.n; k++) + { + if(k_flag->a.a[k] > 0) fprintf(fp, "(%d) %u\n", k, k_flag->a.a[k]); + } + } + + if(hp_k_mer > ava_k_mer*hp_rate) return 1; ///must use '>' instead of '>=' + r->a.n = 0; + return 0; +} void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, diff --git a/Correct.h b/Correct.h index 48af7f9..8bb4735 100644 --- a/Correct.h +++ b/Correct.h @@ -1173,6 +1173,7 @@ void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_I UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read); long long get_affine_gap_score(overlap_region* ovc, UC_Read* g_read, UC_Read* overlap_read, uint8_t* x_num, uint8_t* y_num, uint64_t EstimateXOlen, uint64_t EstimateYOlen); +int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u32_warp* b, kvec_t_u64_warp* r, kvec_t_u8_warp* k_flag, float hp_rate, FILE* fp); #define FORWARD_KSW 0 #define BACKWARD_KSW 1 diff --git a/Hash_Table.cpp b/Hash_Table.cpp index b8a3318..fdda1e0 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -391,6 +391,16 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul resize_Chain_Data(dp, a_n); ret = ha_chain_check(a, a_n, dp, min_score, band_width_threshold); + /***************************************debug**************************************/ + if(a_n > 0 && Get_NAME_LENGTH((R_INF),a[0].readID)==strlen("m64062_190803_042216/128778853/ccs")) + { + if (memcmp("m64062_190803_042216/128778853/ccs", Get_NAME((R_INF), a[0].readID), + Get_NAME_LENGTH((R_INF), a[0].readID)) == 0) + { + fprintf(stderr, "ret: %d\n", ret); + } + } + /***************************************debug**************************************/ if (ret > 0) { a_n = ret; goto skip_dp; @@ -500,6 +510,17 @@ skip_dp: } + /***************************************debug**************************************/ + if(a_n > 0 && Get_NAME_LENGTH((R_INF),a[0].readID)==strlen("m64062_190803_042216/128778853/ccs")) + { + if (memcmp("m64062_190803_042216/128778853/ccs", Get_NAME((R_INF), a[0].readID), + Get_NAME_LENGTH((R_INF), a[0].readID)) == 0) + { + fprintf(stderr, "max_i: %lld\n", max_i); + } + } + /***************************************debug**************************************/ + clear_fake_cigar(&(result->f_cigar)); ///note a has been sorted by offset, that means has been sorted by query offset i = max_i; diff --git a/Overlaps.cpp b/Overlaps.cpp index b9b542d..20a3aca 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -2107,6 +2107,37 @@ 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) +{ + fprintf(stderr, "\n\n\n"); + uint32_t v, n_vtx = g->n_seq * 2, queryLen = strlen(name), flag = 0; + 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); + if(g->seq[v>>1].del) + { + fprintf(stderr, "read has already been deleted.\n"); + return; + } + + asg_arc_t *av = asg_arc_a(g, v); + uint32_t i, nv = asg_arc_n(g, v); + for (i = 0; i < nv; ++i) + { + fprintf(stderr, "target: %.*s, el: %u, strong: %u, ol: %u, del: %u\n", + (int)Get_NAME_LENGTH(R_INF, (av[i].v>>1)), + Get_NAME(R_INF, (av[i].v>>1)), + av[i].el, av[i].strong, av[i].ol, av[i].del); + } + flag = 1; + } + } +} + asg_t *ma_sg_gen(const ma_hit_t_alloc* sources, long long n_read, const ma_sub_t *coverage_cut, int max_hang, int min_ovlp) @@ -8206,13 +8237,96 @@ int asg_arc_del_short_diploi_by_suspect_edge(asg_t *g, int max_ext) return n_cut; } +int if_potential_false_node(ma_hit_t_alloc* paf, uint64_t rLen, ma_sub_t* max_left, ma_sub_t* max_right, int if_set) +{ + max_left->s = max_right->s = rLen; + max_left->e = max_right->e = 0; + + long long j; + uint32_t qs, qe; + for (j = 0; j < paf->length; j++) + { + if(paf->buffer[j].del) continue; + if(paf->buffer[j].el != 1) continue; + + qs = Get_qs(paf->buffer[j]); + qe = Get_qe(paf->buffer[j]); + + ///overlaps from left side + if(qs == 0) + { + if(qs < max_left->s) max_left->s = qs; + if(qe > max_left->e) max_left->e = qe; + } + + ///overlaps from right side + if(qe == rLen) + { + if(qs < max_right->s) max_right->s = qs; + if(qe > max_right->e) max_right->e = qe; + } + + ///note: if (qs == 0 && qe == rLen) + ///this overlap would be added to both b_left and b_right + ///that is what we want + } + + if (max_left->e > max_right->s) return 0; + + long long new_left_e, new_right_s; + new_left_e = max_left->e; + new_right_s = max_right->s; + + for (j = 0; j < paf->length; j++) + { + if(paf->buffer[j].del) continue; + if(paf->buffer[j].el != 1) continue; + + qs = Get_qs(paf->buffer[j]); + qe = Get_qe(paf->buffer[j]); + ///check contained overlaps + if(qs != 0 && qe != rLen) + { + ///[qs, qe), [max_left.s, max_left.e) + if(qs < max_left->e && qe > max_left->e) + { + ///if(qe > max_left->e) max_left->e = qe; + if(qe > max_left->e && qe > new_left_e) new_left_e = qe; + } + + ///[qs, qe), [max_right.s, max_right.e) + if(qs < max_right->s && qe > max_right->s) + { + ///if(qs < max_right->s) max_right->s = qs; + if(qs < max_right->s && qs < new_right_s) new_right_s = qs; + } + + + + if(if_set) + { + max_left->e = new_left_e; + max_right->s = new_right_s; + } + } + } + + max_left->e = new_left_e; + max_right->s = new_right_s; + + if (max_left->e > max_right->s) return 0; + + return 1; +} + //reomve edge between two chromesomes //this node must be a single read -int asg_arc_del_false_node(asg_t *g, int max_ext) +int asg_arc_del_false_node(asg_t *g, ma_hit_t_alloc* sources, int max_ext) { double startTime = Get_T(); kvec_t(uint64_t) b; memset(&b, 0, sizeof(b)); + ma_sub_t max_left, max_right; uint32_t v, n_vtx = g->n_seq * 2; long long n_cut = 0; @@ -8242,6 +8356,11 @@ int asg_arc_del_false_node(asg_t *g, int max_ext) continue; } + if(if_potential_false_node(&(sources[v>>1]), g->seq[v>>1].len, &max_left, &max_right, 1)==0) + { + continue; + } + asg_arc_t *av = asg_arc_a(g, v); kv_push(uint64_t, b, (uint64_t)((uint64_t)av[0].ol << 32 | (av - g->arc))); } @@ -9582,35 +9701,6 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source -void debug_info_of_specfic_node(const char* name, asg_t *g, char* command) -{ - fprintf(stderr, "\n\n\n"); - uint32_t v, n_vtx = g->n_seq * 2, queryLen = strlen(name); - 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) - { - fprintf(stderr, "\nafter %s\n****************graph ref_read: %.*s, dir: %u****************\n", - command, (int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), v&1); - if(g->seq[v>>1].del) - { - fprintf(stderr, "read has already been deleted.\n"); - return; - } - - asg_arc_t *av = asg_arc_a(g, v); - uint32_t i, nv = asg_arc_n(g, v); - for (i = 0; i < nv; ++i) - { - fprintf(stderr, "target: %.*s, el: %u, strong: %u, ol: %u, del: %u\n", - (int)Get_NAME_LENGTH(R_INF, (av[i].v>>1)), - Get_NAME(R_INF, (av[i].v>>1)), - av[i].el, av[i].strong, av[i].ol, av[i].del); - } - return; - } - } -} void debug_info_of_specfic_read(const char* name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int id, const char* command) @@ -9636,8 +9726,8 @@ ma_hit_t_alloc* reverse_sources, int id, const char* command) { fprintf(stderr, "\n\n\nafter %s\n", command); - fprintf(stderr, "****************ma_hit_t (%lld)ref_read: %.*s****************\n", - i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + fprintf(stderr, "****************ma_hit_t (%lld)ref_read: %.*s, len: %lu****************\n", + i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), Get_READ_LENGTH(R_INF, i)); fprintf(stderr, "sources Len: %d, is_fully_corrected: %d\n", @@ -12954,12 +13044,12 @@ float drop_ratio) #define T_ROUND 2 asg_t *g = ug->g; int round = T_ROUND; - + redo: asg_pop_bubble_primary_trio(ug, bubble_dist, (uint32_t)-1, DROP); untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, DROP); - + if(just_bubble_pop == 0) { cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, @@ -13013,7 +13103,7 @@ float drop_ratio) { if(round != T_ROUND) { - unitig_arc_del_short_diploid_by_length(ug->g, drop_ratio); + unitig_arc_del_short_diploid_by_length(ug->g, drop_ratio); } round--; goto redo; @@ -22764,7 +22854,7 @@ kvec_asg_arc_t_warp* new_rtg_edges) delete_useless_nodes(ug); renew_utg(ug, read_g, new_rtg_edges); } - + if (!(asm_opt.flag & HA_F_BAN_POST_JOIN)) { rescue_missing_overlaps_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang, @@ -26792,9 +26882,6 @@ int max_hang, int min_ovlp) } } - - - void clean_graph( int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long n_read, uint64_t* readLen, long long mini_overlap_length, @@ -26815,7 +26902,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g) ///normalize_ma_hit_t_single_side(sources, n_read); normalize_ma_hit_t_single_side_advance(sources, n_read); normalize_ma_hit_t_single_side_advance(reverse_sources, n_read); - if (ha_opt_triobin(&asm_opt)) { @@ -26840,11 +26926,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); - - ///debug_info_of_specfic_read("m54329U_190827_173812/166332272/ccs", sources, reverse_sources, -1, "clean"); - + sg = ma_sg_gen(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length); asg_arc_del_trans(sg, gap_fuzz); + asm_opt.coverage = get_coverage(sources, coverage_cut, n_read); if(VERBOSE >= 1) @@ -26856,9 +26941,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g) } asg_cut_tip(sg, asm_opt.max_short_tip); - + ///debug_info_of_specfic_node("m64062_190803_042216/15205346/ccs", sg, "inner_1"); ///drop_inexact_edegs_at_bubbles(sg, bubble_dist); - ///debug_info_of_specfic_node("m54329U_190827_173812/166332272/ccs", sg, "beg"); if(clean_round > 0) { @@ -26896,7 +26980,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) asg_arc_identify_simple_bubbles_multi(sg, 1); //reomve edge between two chromesomes //this node must be a single read - asg_arc_del_false_node(sg, asm_opt.max_short_tip); + asg_arc_del_false_node(sg, sources, asm_opt.max_short_tip); asg_cut_tip(sg, asm_opt.max_short_tip); /****************************may have bugs********************************/ /****************************may have bugs********************************/ @@ -27044,10 +27128,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; - - - // debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs", sources, reverse_sources, -1, "clean"); - + ///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; @@ -27073,6 +27154,9 @@ long long bubble_dist, int read_graph, int write) &R_INF, output_file_name); } + + ///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 635b7c8..316a5d6 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -46,6 +46,7 @@ void destory_All_reads(All_reads* r) if (r->read_sperate[i]) free(r->read_sperate[i]); if (r->paf && r->paf[i].buffer) free(r->paf[i].buffer); if (r->reverse_paf && r->reverse_paf[i].buffer) free(r->reverse_paf[i].buffer); + ///if (r->pb_regions) kv_destroy(r->pb_regions[i].a); } free(r->paf); free(r->reverse_paf); @@ -55,6 +56,7 @@ void destory_All_reads(All_reads* r) free(r->name_index); free(r->read_length); free(r->trio_flag); + ///if (r->pb_regions) free(r->pb_regions); } void write_All_reads(All_reads* r, char* read_file_name) @@ -192,6 +194,7 @@ int load_All_reads(All_reads* r, char* read_file_name) r->second_round_cigar[i].lost_base_length = r->cigars[i].lost_base_length = 0; r->second_round_cigar[i].lost_base = r->cigars[i].lost_base = NULL; } + ///r->pb_regions = NULL; free(index_name); fclose(fp); @@ -263,6 +266,8 @@ void malloc_All_reads(All_reads* r) r->second_round_cigar = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads); r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads); r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads); + ///r->pb_regions = (kvec_t_u64_warp*)malloc(r->total_reads*sizeof(kvec_t_u64_warp)); + for (i = 0; i < (long long)r->total_reads; i++) { r->second_round_cigar[i].size = r->cigars[i].size = 0; @@ -274,6 +279,7 @@ void malloc_All_reads(All_reads* r) r->second_round_cigar[i].lost_base = r->cigars[i].lost_base = NULL; init_ma_hit_t_alloc(&(r->paf[i])); init_ma_hit_t_alloc(&(r->reverse_paf[i])); + ///kv_init(r->pb_regions[i].a); } r->name = (char*)malloc(sizeof(char)*r->total_name_length); @@ -667,3 +673,51 @@ void reverse_complement(char* pattern, uint64_t length) pattern[end] = RC_CHAR(pattern[end]); } } + + +void init_Debug_reads(Debug_reads* x, const char* file) +{ + int nameLen, i, bufLen = 1000; + if((uint64_t)(bufLen) < strlen(file) + 50) bufLen = strlen(file) + 50; + char* Name_Buffer = (char*)malloc(sizeof(char)*bufLen); + fprintf(stderr, "Queried debugging reads at: %s\n", file); + + x->fp = fopen(file,"r"); + x->query_num = 0; + + while(fgets(Name_Buffer, bufLen, x->fp)) + { + x->query_num++; + } + x->read_name = (char**)malloc(sizeof(char*)*x->query_num); + fseek(x->fp, 0, SEEK_SET); + + i = 0; + while(fgets(Name_Buffer, bufLen, x->fp)) + { + nameLen = strlen(Name_Buffer) - 1; + 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'; + i++; + } + + fclose(x->fp); + + sprintf(Name_Buffer, "%s.debug.stdout", file); + x->fp = fopen(Name_Buffer,"w"); + fprintf(stderr, "Print debugging information to: %s\n", Name_Buffer); + free(Name_Buffer); +} + +void destory_Debug_reads(Debug_reads* x) +{ + uint64_t i; + for (i = 0; i < x->query_num; i++) + { + free(x->read_name[i]); + } + + free(x->read_name); + fclose(x->fp); +} \ No newline at end of file diff --git a/Process_Read.h b/Process_Read.h index 5af0c3d..fc29101 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -93,6 +93,7 @@ typedef struct uint32_t new_length; } Compressed_Cigar_record; + #define AMBIGU 0 #define FATHER 1 #define MOTHER 2 @@ -128,6 +129,8 @@ typedef struct ma_hit_t_alloc* paf; ma_hit_t_alloc* reverse_paf; + + ///kvec_t_u64_warp* pb_regions; } All_reads; extern All_reads R_INF; @@ -140,6 +143,14 @@ typedef struct long long RID; } UC_Read; +typedef struct +{ + char** read_name; + uint64_t query_num; + FILE* fp; + pthread_mutex_t OutputMutex; +} Debug_reads; + void init_All_reads(All_reads* r); void malloc_All_reads(All_reads* r); void ha_insert_read_len(All_reads *r, int read_len, int name_len); @@ -154,5 +165,6 @@ void write_All_reads(All_reads* r, char* read_file_name); int load_All_reads(All_reads* r, char* read_file_name); void destory_All_reads(All_reads* r); int destory_read_bin(All_reads* r); - +void init_Debug_reads(Debug_reads* x, const char* file); +void destory_Debug_reads(Debug_reads* x); #endif diff --git a/anchor.cpp b/anchor.cpp index f8aa2c1..a1f9ed3 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -28,8 +28,8 @@ typedef struct { } seed1_t; struct ha_abuf_s { - uint64_t n_a, m_a; - uint32_t old_mz_m; + uint64_t n_a, m_a;///number of anchors (seed positions) + uint32_t old_mz_m;///number of seeds ha_mz1_v mz; seed1_t *seed; anchor1_t *a; @@ -57,10 +57,9 @@ static int ha_ov_type(const overlap_region *r, uint32_t len) else return r->x_pos_s == 0? 0 : 1; } -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) +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, +void *ha_flt_tab, ha_pt_t *ha_idx) { - extern void *ha_flt_tab; - extern ha_pt_t *ha_idx; uint32_t i, rlen; uint64_t k, l; double low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; @@ -74,7 +73,8 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg rlen = Get_READ_LENGTH(R_INF, rid); // read length // get the list of anchors - ha_sketch(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); + 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); + // minimizer of queried read if (ab->mz.m > ab->old_mz_m) { ab->old_mz_m = ab->mz.m; REALLOC(ab->seed, ab->old_mz_m); @@ -93,6 +93,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg } for (i = 0, k = 0; i < ab->mz.n; ++i) { int j; + ///z is one of the minimizer ha_mz1_t *z = &ab->mz.a[i]; seed1_t *s = &ab->seed[i]; for (j = 0; j < s->n; ++j) { @@ -116,6 +117,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg } } + // copy over to _cl_ if (ab->m_a >= (uint64_t)cl->size) { cl->size = ab->m_a; @@ -128,6 +130,18 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg p->offset = ab->a[k].other_off; p->self_offset = ab->a[k].self_off; p->good = ab->a[k].good; + + /***************************************debug**************************************/ + if(Get_NAME_LENGTH((R_INF),p->readID)==strlen("m64062_190803_042216/128778853/ccs")) + { + if (memcmp("m64062_190803_042216/128778853/ccs", Get_NAME((R_INF), p->readID), + Get_NAME_LENGTH((R_INF), p->readID)) == 0) + { + fprintf(stderr, "(%lu) readID: %u, strand: %u, offset: %u, self_offset: %u\n", + k, p->readID, p->strand, p->offset, p->self_offset); + } + } + /***************************************debug**************************************/ } cl->length = ab->n_a; @@ -172,10 +186,130 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg } } - ks_introsort_or_xs(overlap_list->length, overlap_list->list); + ///ks_introsort_or_xs(overlap_list->length, overlap_list->list); } +void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf) +{ + uint64_t j = 0, inner_j = 0; + while (j < overlap_list->length && inner_j < paf->length) + { + if(overlap_list->list[j].y_id < paf->buffer[inner_j].tn) + { + j++; + } + else if(overlap_list->list[j].y_id > paf->buffer[inner_j].tn) + { + inner_j++; + } + else + { + if(overlap_list->list[j].y_pos_strand == paf->buffer[inner_j].rev) + { + overlap_list->list[j].is_match = 1; + } + j++; + inner_j++; + } + } +} + + +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, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf) +{ + 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, ha_flt_tab, ha_idx); + + if(ha_idx_hp) + { + uint32_t i, k, y_id, overlapLen, max_i; + int shared_seed; + overlap_region t; + overlap_region_sort_y_id(overlap_list->list, overlap_list->length); + ma_hit_sort_tn(paf->buffer, paf->length); + ma_hit_sort_tn(rev_paf->buffer, rev_paf->length); + lable_matched_ovlp(overlap_list, paf); + lable_matched_ovlp(overlap_list, rev_paf); + + for (i = 0, k = 0; i < overlap_list->length; ++i) + { + if(overlap_list->list[i].is_match == 1) + { + if(k != i) + { + t = overlap_list->list[k]; + overlap_list->list[k] = overlap_list->list[i]; + overlap_list->list[i] = t; + overlap_list->list[k].is_match = 0; + } + k++; + } + } + 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, ha_flt_tab_hp, ha_idx_hp); + + if(overlap_list->length + overlap_list_hp->length > overlap_list->size) + { + overlap_list->list = (overlap_region*)realloc(overlap_list->list, + sizeof(overlap_region)*(overlap_list->length + overlap_list_hp->length)); + memset(overlap_list->list + overlap_list->size, 0, sizeof(overlap_region)* + (overlap_list->length + overlap_list_hp->length - overlap_list->size)); + overlap_list->size = overlap_list->length + overlap_list_hp->length; + } + + for (i = 0, k = overlap_list->length; i < overlap_list_hp->length; i++, k++) + { + t = overlap_list->list[k]; + overlap_list->list[k] = overlap_list_hp->list[i]; + overlap_list_hp->list[i] = t; + } + overlap_list->length = k; + + overlap_region_sort_y_id(overlap_list->list, overlap_list->length); + + i = k = 0; + while (i < overlap_list->length) + { + y_id = overlap_list->list[i].y_id; + shared_seed = overlap_list->list[i].shared_seed; + overlapLen = overlap_list->list[i].overlapLen; + max_i = i; + i++; + while (i < overlap_list->length && overlap_list->list[i].y_id == y_id) + { + if((overlap_list->list[i].shared_seed > shared_seed) || + ((overlap_list->list[i].shared_seed == shared_seed) && (overlap_list->list[i].overlapLen <= overlapLen))) + { + y_id = overlap_list->list[i].y_id; + shared_seed = overlap_list->list[i].shared_seed; + overlapLen = overlap_list->list[i].overlapLen; + max_i = i; + } + i++; + } + + if(k != max_i) + { + t = overlap_list->list[k]; + overlap_list->list[k] = overlap_list->list[max_i]; + overlap_list->list[max_i] = t; + } + k++; + } + + overlap_list->length = k; + } + + ks_introsort_or_xs(overlap_list->length, overlap_list->list); +} void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list) { diff --git a/htab.cpp b/htab.cpp index 8442efb..0e1ec4f 100644 --- a/htab.cpp +++ b/htab.cpp @@ -35,6 +35,8 @@ const unsigned char seq_nt4_table[256] = { // translate ACGT to 0123 void *ha_flt_tab; ha_pt_t *ha_idx; +void *ha_flt_tab_hp; +ha_pt_t *ha_idx_hp; /*************************** * Yak specific parameters * @@ -71,7 +73,7 @@ typedef struct { int n_shift, n_hashes; uint8_t *b; } yak_bf_t; - +///in most cases, n_shift = 25, n_hashes = 4 yak_bf_t *yak_bf_init(int n_shift, int n_hashes) { yak_bf_t *b; @@ -126,10 +128,12 @@ typedef struct { typedef struct { int k, pre, n_hash, n_shift; - uint64_t tot; + uint64_t tot; ///number of distinct k-mers 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 static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift) { ha_ct_t *h; @@ -138,12 +142,16 @@ static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift) CALLOC(h, 1); h->k = k, h->pre = pre; CALLOC(h->h, 1<pre); + ///ipre = 4096 + ///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 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) - h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash); + h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash); ///h->n_shift = 37, h->pre = 12, h->n_hash = 4 } return h; } @@ -173,15 +181,23 @@ static int ha_ct_insert_list(ha_ct_t *h, int create_new, int n, const uint64_t * int j, mask = (1<pre) - 1, n_ins = 0; ha_ct1_t *g; if (n == 0) return 0; + ///corresponding small hash index g = &h->h[a[0]&mask]; for (j = 0; j < n; ++j) { int ins = 1, absent; + ///x is a 64-bit word, h->pre=12 + ///all elements at a have the same low 12 bits + ///so low 12 bits are not useful uint64_t x = a[j] >> h->pre; khint_t k; if ((a[j]&mask) != (a[0]&mask)) continue; if (create_new) { + ///for 0-th counting, g->b = NULL if (g->b) ins = (yak_bf_insert(g->b, x) == h->n_hash); + ///for 0-th counting, g->b = NULL + ///x = the high 52 bits of a[j] + low 12 bits 0 + ///the low 12 bits are used for counting if (ins) { k = yak_ct_put(g->h, x << YAK_COUNTER_BITS | (g->b? 1 : 0), &absent); if (absent) ++n_ins; @@ -219,6 +235,8 @@ static void worker_ct_hist(void *data, long i, int tid) // callback for kt_for() ++cnt[kh_key(g, k)&YAK_MAX_COUNT]; } +///YAK_N_COUNTS is also 4096 +///used for calculating k-mer histogram static void ha_ct_hist(const ha_ct_t *h, int64_t cnt[YAK_N_COUNTS], int n_thread) { hist_aux_t a; @@ -226,6 +244,7 @@ static void ha_ct_hist(const ha_ct_t *h, int64_t cnt[YAK_N_COUNTS], int n_thread a.h = h; memset(cnt, 0, YAK_N_COUNTS * sizeof(uint64_t)); CALLOC(a.cnt, n_thread); + ///start 4096 threads kt_for(n_thread, worker_ct_hist, &a, 1<pre); for (i = 0; i < YAK_N_COUNTS; ++i) cnt[i] = 0; for (j = 0; j < n_thread; ++j) @@ -265,6 +284,7 @@ static void ha_ct_shrink(ha_ct_t *h, int min, int max, int n_thread) int i; shrink_aux_t a; a.h = h, a.min = min, a.max = max; + ///still start 4096 threads kt_for(n_thread, worker_ct_shrink, &a, 1<pre); for (i = 0, h->tot = 0; i < 1<pre; ++i) h->tot += kh_size(h->h[i].h); @@ -306,6 +326,7 @@ static void worker_pt_gen(void *data, long i, int tid) // callback for kt_for() int absent; khint_t l; l = yak_pt_put(b->h, kh_key(g, k) >> a->ct->pre << YAK_COUNTER_BITS, &absent); + ///this should be the start index of kh_key's corresponding pos at ha_idxpos_t* a kh_val(b->h, l) = b->n; b->n += kh_key(g, k) & YAK_MAX_COUNT; } @@ -412,8 +433,12 @@ typedef struct { ha_mz1_t *b; } ch_buf_t; +///p = 12 static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k-mer $y to a linear buffer { + ///assign k-mer to one of the 4096 bins + ///using low 12 bits for assigning + ///so all elements at b have the same low 12 bits int pre = y & ((1<n == b->m) { @@ -425,6 +450,7 @@ static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k static inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y) { + ///assign minimizer to one of 4096 bins by low 12 bits int pre = y->x & ((1<n == b->m) { @@ -434,13 +460,17 @@ static inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y) b->b[b->n++] = *y; } +///buf is the read block, k is the k-mer length, p = 12, len is the read length, seq is the read static void count_seq_buf(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf { int i, l; uint64_t x[4], mask = (1ULL<>1)) & mask; x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; @@ -483,12 +513,13 @@ KSEQ_INIT(gzFile, gzread) #define HAF_RS_WRITE_SEQ 0x8 #define HAF_RS_READ 0x10 #define HAF_CREATE_NEW 0x20 +#define HAF_SKIP_READ 0x40 typedef struct { // global data structure for kt_pipeline() const yak_copt_t *opt; const void *flt_tab; int flag, create_new, is_store; - uint64_t n_seq; + uint64_t n_seq; ///number of total reads kseq_t *ks; UC_Read ucr; ha_ct_t *ct; @@ -499,7 +530,8 @@ typedef struct { // global data structure for kt_pipeline() typedef struct { // data structure for each step in kt_pipeline() pl_data_t *p; - uint64_t n_seq0; + uint64_t n_seq0; ///the start index of current buffer block at R_INF + ///sum_len = total bases, nk = number of k-mers int n_seq, m_seq, sum_len, nk; int *len; char **seq; @@ -514,15 +546,17 @@ static void worker_for_insert(void *data, long i, int tid) // callback for kt_fo ch_buf_t *b = &s->buf[i]; if (s->p->pt) b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b); - else + else///for 0-th count, go into here b->n_ins += ha_ct_insert_list(s->p->ct, s->p->create_new, b->n, b->a); } static void worker_for_mz(void *data, long i, int tid) { st_data_t *s = (st_data_t*)data; + ///get the corresponding minimzer vector of this read ha_mz1_v *b = &s->mz_buf[tid]; s->mz_buf[tid].n = 0; + ///s->p->opt->w = 51, s->p->opt->k ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab); s->mz[i].n = s->mz[i].m = b->n; MALLOC(s->mz[i].a, b->n); @@ -540,6 +574,11 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip s->n_seq0 = p->n_seq; 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) + { + ++p->n_seq; + continue; + } int l; recover_UC_Read(&p->ucr, p->rs_in, p->n_seq); l = p->ucr.length; @@ -565,6 +604,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip exit(1); } if (p->rs_out) { + ///for 0-th count, just insert read length to R_INF, instead of read if (p->flag & HAF_RS_WRITE_LEN) { assert(p->n_seq == p->rs_out->total_reads); ha_insert_read_len(p->rs_out, l, p->ks->name.l); @@ -578,6 +618,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l); } } + ///for 0-th count, insert both seq and length to local block if (s->n_seq == s->m_seq) { s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1); REALLOC(s->len, s->m_seq); @@ -589,6 +630,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip ++p->n_seq; s->sum_len += l; s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0; + ///p->opt->chunk_size is the block max size if (s->sum_len >= p->opt->chunk_size) break; } @@ -596,18 +638,24 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip if (s->sum_len == 0) free(s); else return s; } else if (step == 1) { // step 2: extract k-mers + ///s is the block of reads st_data_t *s = (st_data_t*)in; + ///for 0-th counting, n_pre = 4096 int i, n_pre = 1<opt->pre, m; // allocate the k-mer buffer CALLOC(s->buf, n_pre); m = (int)(s->nk * 1.2 / n_pre) + 1; + //pre-allocate memory for each of 4096 buffer for (i = 0; i < n_pre; ++i) { s->buf[i].m = m; + ///for 0-th counting, p->pt = NULL if (p->pt) MALLOC(s->buf[i].b, m); else MALLOC(s->buf[i].a, m); } // fill the buffer + ///for 0-th counting, p->opt->w == 1 if (p->opt->w == 1) { // enumerate all k-mers + ///scan all reads for (i = 0; i < s->n_seq; ++i) { if (p->opt->is_HPC) count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]); @@ -618,18 +666,21 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip } else { // minimizers only uint32_t j; // compute minimizers + // s->n_seq is how many reads at this buffer + // s->mz && s->mz_buf are lists of minimzer vectors CALLOC(s->mz, s->n_seq); CALLOC(s->mz_buf, p->opt->n_thread); + ///calculate minimzers for each read, each read corresponds to one thread kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq); for (i = 0; i < p->opt->n_thread; ++i) free(s->mz_buf[i].a); free(s->mz_buf); // insert minimizers - if (p->pt) { + if (p->pt) {///insert whole minimizer for (i = 0; i < s->n_seq; ++i) for (j = 0; j < s->mz[i].n; ++j) pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]); - } else { + } else {///just insert the hash key of minimizer for (i = 0; i < s->n_seq; ++i) for (j = 0; j < s->mz[i].n; ++j) ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x); @@ -640,6 +691,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip } free(s->mz); } + ///just clean seq free(s->seq); free(s->len); s->seq = 0, s->len = 0; return s; @@ -647,7 +699,9 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip st_data_t *s = (st_data_t*)in; int i, n = 1<opt->pre; uint64_t n_ins = 0; + ///for 0-th counting, p->pt = NULL kt_for(p->opt->n_thread, worker_for_insert, s, n); + ///n_ins is number of distinct k-mers for (i = 0; i < n; ++i) { n_ins += s->buf[i].n_ins; if (p->pt) free(s->buf[i].b); @@ -668,6 +722,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt_t *p0, ha_ct_t *c0, const void *flt_tab, All_reads *rs, int64_t *n_seq) { + ///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN|HAF_CREATE_NEW int read_rs = (rs && (flag & HAF_RS_READ)); pl_data_t pl; gzFile fp = 0; @@ -676,23 +731,29 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt if (read_rs) { pl.rs_in = rs; init_UC_Read(&pl.ucr); - } else { + } else {///for 0-th counting, go into here if ((fp = gzopen(fn, "r")) == 0) return 0; pl.ks = kseq_init(fp); } + ///for 0-th counting, read all reads into pl.rs_out if (rs && (flag & (HAF_RS_WRITE_LEN|HAF_RS_WRITE_SEQ))) pl.rs_out = rs; + ///for 0-th counting, flt_tab = NULL + ///for 1-th counting, flt_tab = NULL pl.flt_tab = flt_tab; pl.opt = opt; pl.flag = flag; - if (p0) { + if (p0) {///for 1-th counting, p0 = NULL pl.pt = p0, pl.create_new = 0; // never create new elements in a position table assert(p0->k == opt->k && p0->pre == opt->pre); } else if (c0) { pl.ct = c0, pl.create_new = !!(flag&HAF_CREATE_NEW); assert(c0->k == opt->k && c0->pre == opt->pre); - } else { + } else {///for 0-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 + ///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); } kt_pipeline(3, worker_count, &pl, 3); @@ -713,6 +774,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi yak_copt_t opt; ha_ct_t *h = 0; assert(!(flag & HAF_RS_WRITE_LEN) || !(flag & HAF_RS_WRITE_SEQ)); // not both + ///for 0-th counting, flag = HAF_COUNT_ALL|HAF_RS_WRITE_LEN if (rs) { if (flag & HAF_RS_WRITE_LEN) init_All_reads(rs); @@ -721,10 +783,16 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi } yak_copt_init(&opt); 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 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 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 for (i = 0; i < asm_opt->num_reads; ++i) h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, &n_seq); if (h && opt.bf_shift > 0) @@ -775,19 +843,24 @@ void ha_ft_destroy(void *h) * High-level interfaces * *************************/ -void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov) +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; int64_t cnt[YAK_N_COUNTS]; - int peak_hom, peak_het, cutoff; + 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, NULL, NULL, rs); - ha_ct_hist(h, cnt, asm_opt->thread_num); - peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het); - if (hom_cov) *hom_cov = peak_hom; - if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); - cutoff = (int)(peak_hom * asm_opt->high_factor); - if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1; + h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs); + if(!(ex_flag & HAF_SKIP_READ)) + { + ha_ct_hist(h, cnt, asm_opt->thread_num); + peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het); + if (hom_cov) *hom_cov = peak_hom; + if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); + ///in default, asm_opt->high_factor = 5.0 + cutoff = (int)(peak_hom * 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); ha_ct_destroy(h); @@ -796,21 +869,23 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov) return (void*)flt_tab; } -ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov, int *het_cov) +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) { int64_t cnt[YAK_N_COUNTS], tot_cnt; int peak_hom, peak_het, i, extra_flag1, extra_flag2; ha_ct_t *ct; ha_pt_t *pt; - if (read_from_store) { + if (read_from_store) {///if reads have already been read extra_flag1 = extra_flag2 = HAF_RS_READ; - } else if (rs->total_reads == 0) { + } else if (rs->total_reads == 0) {///if reads & length have not been scanned extra_flag1 = HAF_RS_WRITE_LEN; extra_flag2 = HAF_RS_WRITE_SEQ; - } else { + } else {///if length has been loaded but reads have not extra_flag1 = HAF_RS_WRITE_SEQ; extra_flag2 = HAF_RS_READ; } + if(is_hp_mode) extra_flag1 |= HAF_SKIP_READ, extra_flag2 |= HAF_SKIP_READ; + ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__, yak_realtime(), yak_cpu_usage(), (long)ct->tot); @@ -820,12 +895,16 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f if (hom_cov) *hom_cov = peak_hom; if (het_cov) *het_cov = peak_het; if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); + ///here ha_ct_shrink is mostly used to remove k-mer appearing only 1 time if (flt_tab == 0) { int cutoff = (int)(peak_hom * asm_opt->high_factor); if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1; + if((extra_flag1 & HAF_SKIP_READ) && (extra_flag2 & HAF_SKIP_READ)) cutoff = YAK_MAX_COUNT - 1; ha_ct_shrink(ct, 2, cutoff, asm_opt->thread_num); for (i = 2, tot_cnt = 0; i <= cutoff; ++i) tot_cnt += cnt[i] * i; } 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); for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i; } @@ -837,3 +916,162 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos); return pt; } + + + +int write_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); + FILE* fp = fopen(gfa_name, "w"); + if (!fp) { + free(gfa_name); + return 0; + } + yak_ft_t *ha_flt_tab = (yak_ft_t*)flt_tab; + + if(ha_flt_tab) + { + fwrite("f", 1, 1, fp); + yak_ft_save(ha_flt_tab, fp); + } + + + if(ha_idx) + { + int i; + ha_pt1_t *g; + fwrite("h", 1, 1, fp); + fwrite(&ha_idx->k, sizeof(ha_idx->k), 1, fp); + fwrite(&ha_idx->pre, sizeof(ha_idx->pre), 1, fp); + fwrite(&ha_idx->tot, sizeof(ha_idx->tot), 1, fp); + fwrite(&ha_idx->tot_pos, sizeof(ha_idx->tot_pos), 1, fp); + + for (i = 0; i < 1<pre; ++i) + { + g = &(ha_idx->h[i]); + yak_pt_save(g->h, fp); + fwrite(&g->n, sizeof(g->n), 1, fp); + fwrite(g->a, sizeof(ha_idxpos_t), g->n, fp); + } + } + + fwrite(&opt->number_of_round, sizeof(opt->number_of_round), 1, fp); + fwrite(&opt->hom_cov, sizeof(opt->hom_cov), 1, fp); + fwrite(&opt->het_cov, sizeof(opt->het_cov), 1, fp); + fwrite(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp); + + + write_All_reads(r, gfa_name); + + fprintf(stderr, "[M::%s] Index has been written.\n", __func__); + free(gfa_name); + fclose(fp); + 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) +{ + char* gfa_name = (char*)malloc(strlen(file_name)+25); + sprintf(gfa_name, "%s.pt_flt", file_name); + FILE* fp = fopen(gfa_name, "r"); + if (!fp) { + free(gfa_name); + return 0; + } + + ha_pt_t *ha_idx = NULL; + char mode = 0; + int f_flag, absent, i; + double index_time, index_s_time, pos_time, pos_s_time; + + + + f_flag += fread(&mode, 1, 1, fp); + if(mode == 'f') + { + index_time = yak_realtime(); + + yak_ft_load((yak_ft_t **)r_flt_tab, fp); + + f_flag += fread(&mode, 1, 1, fp); + + fprintf(stderr, "[M::%s::%.3f] ==> Loaded flt table\n", __func__, yak_realtime()-index_time); + } + ///insert using multiple threads??? + if(mode == 'h') + { + pos_time = index_time = 0; + + CALLOC(ha_idx, 1); + ha_pt1_t *g; + f_flag += fread(&ha_idx->k, sizeof(ha_idx->k), 1, fp); + f_flag += fread(&ha_idx->pre, sizeof(ha_idx->pre), 1, fp); + f_flag += fread(&ha_idx->tot, sizeof(ha_idx->tot), 1, fp); + f_flag += fread(&ha_idx->tot_pos, sizeof(ha_idx->tot_pos), 1, fp); + CALLOC(ha_idx->h, 1<pre); + for (i = 0; i < 1<pre; ++i) + { + index_s_time = yak_realtime(); + + g = &(ha_idx->h[i]); + yak_pt_load(&(g->h), fp); + + index_time += yak_realtime() - index_s_time; + + pos_s_time = yak_realtime(); + + f_flag += fread(&g->n, sizeof(g->n), 1, fp); + MALLOC(g->a, g->n); + f_flag += fread(g->a, sizeof(ha_idxpos_t), g->n, fp); + + pos_time += yak_realtime() - pos_s_time; + } + (*r_ha_idx) = ha_idx; + + fprintf(stderr, "[M::%s::%.3f(index)/%.3f(pos)] ==> Loaded pos table\n", __func__, index_time, pos_time); + } + + if(mode != 'h' && mode != 'f') + { + free(gfa_name); + fclose(fp); + return 0; + } + + + f_flag += fread(&absent, sizeof(absent), 1, fp); + if(absent != opt->number_of_round) + { + fprintf(stderr, "ERROR: different number of rounds!\n"); + exit(1); + } + + f_flag += fread(&opt->hom_cov, sizeof(opt->hom_cov), 1, fp); + f_flag += fread(&opt->het_cov, sizeof(opt->het_cov), 1, fp); + f_flag += fread(&opt->max_n_chain, sizeof(opt->max_n_chain), 1, fp); + + + fclose(fp); + + if(!load_All_reads(r, gfa_name)) + { + free(gfa_name); + return 0; + } + + + memset(r->trio_flag, AMBIGU, r->total_reads*sizeof(uint8_t)); + r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads); + r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads); + for (i = 0; i < (long long)r->total_reads; i++) + { + init_ma_hit_t_alloc(&(r->paf[i])); + init_ma_hit_t_alloc(&(r->reverse_paf[i])); + } + + fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__); + + free(gfa_name); + return 1; +} \ No newline at end of file diff --git a/htab.h b/htab.h index 2655be8..42696f6 100644 --- a/htab.h +++ b/htab.h @@ -6,7 +6,9 @@ #include "CommandLines.h" typedef struct { - uint64_t x; + uint64_t x; ///x is the hash key + ///rid is the read id, pos is the end pos of this minimizer, rev is the direction + ///span is the length of this k-mer. For non-HPC k-mer, span may not be equal to k uint64_t rid:28, pos:27, rev:1, span:8; } ha_mz1_t; @@ -25,15 +27,20 @@ typedef struct ha_abuf_s ha_abuf_t; extern const unsigned char seq_nt4_table[256]; extern void *ha_flt_tab; extern ha_pt_t *ha_idx; +extern void *ha_flt_tab_hp; +extern ha_pt_t *ha_idx_hp; -void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov); +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); void ha_ft_destroy(void *h); -ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov, int *het_cov); +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); 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); + ha_abuf_t *ha_abuf_init(void); void ha_abuf_destroy(ha_abuf_t *ab); uint64_t ha_abuf_mem(const ha_abuf_t *ab); @@ -48,6 +55,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); 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 @@ -76,6 +84,7 @@ static inline uint64_t yak_hash64_64(uint64_t key) static inline uint64_t yak_hash_long(uint64_t x[4]) { + ///compare forward k-mer and reverse complementary strand int j = x[1] < x[3]? 0 : 1; return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]); } diff --git a/khashl.h b/khashl.h index 9381850..4a8f2dc 100644 --- a/khashl.h +++ b/khashl.h @@ -133,6 +133,29 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26 h->count = 0; \ } \ } +#define __KHASHL_IMPL_S_L(SCOPE, HType, prefix, khkey_t) \ + SCOPE khint_t prefix##_save(HType *h, FILE* fp) { \ + if (!h) return 0; \ + khint_t n_buckets = (h->keys? 1U<bits : 0U); \ + fwrite(&n_buckets, sizeof(n_buckets), 1, fp); \ + fwrite(&h->bits, sizeof(h->bits), 1, fp); \ + fwrite(&h->count, sizeof(h->count), 1, fp); \ + fwrite(h->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \ + fwrite(h->keys, sizeof(khkey_t), n_buckets, fp); \ + return 1; \ + } \ + SCOPE khint_t prefix##_load(HType **h, FILE* fp) { \ + (*h) = prefix##_init(); \ + khint_t n_buckets; \ + fread(&n_buckets, sizeof(n_buckets), 1, fp); \ + fread(&(*h)->bits, sizeof((*h)->bits), 1, fp); \ + fread(&(*h)->count, sizeof((*h)->count), 1, fp); \ + (*h)->used = (khint32_t*)kmalloc(__kh_fsize(n_buckets) * sizeof(khint32_t)); \ + (*h)->keys = (khkey_t*)kmalloc(n_buckets * sizeof(khkey_t)); \ + fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \ + fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); \ + return 1; \ + } \ #define __KHASHL_IMPL_GET(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \ SCOPE khint_t prefix##_getp(const HType *h, const khkey_t *key) { \ @@ -245,6 +268,7 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26 #define KHASHL_INIT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \ __KHASHL_TYPE(HType, khkey_t) \ __KHASHL_IMPL_BASIC(SCOPE, HType, prefix) \ + __KHASHL_IMPL_S_L(SCOPE, HType, prefix, khkey_t) \ __KHASHL_IMPL_GET(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \ __KHASHL_IMPL_RESIZE(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \ __KHASHL_IMPL_PUT(SCOPE, HType, prefix, khkey_t, __hash_fn, __hash_eq) \ @@ -264,6 +288,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26 KHASHL_INIT(KH_LOCAL, HType, prefix##_s, HType##_s_bucket_t, prefix##_s_hash, prefix##_s_eq) \ SCOPE HType *prefix##_init(void) { return prefix##_s_init(); } \ SCOPE void prefix##_destroy(HType *h) { prefix##_s_destroy(h); } \ + SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_s_save(h, fp); } \ + SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_s_load(h, fp); } \ SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_s_resize(h, new_n_buckets); } \ SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_s_bucket_t t; t.key = key; return prefix##_s_getp(h, &t); } \ SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_s_del(h, k); } \ @@ -276,6 +302,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26 KHASHL_INIT(KH_LOCAL, HType, prefix##_m, HType##_m_bucket_t, prefix##_m_hash, prefix##_m_eq) \ SCOPE HType *prefix##_init(void) { return prefix##_m_init(); } \ SCOPE void prefix##_destroy(HType *h) { prefix##_m_destroy(h); } \ + SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_m_save(h, fp); } \ + SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_m_load(h, fp); } \ SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_m_resize(h, new_n_buckets); } \ SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_m_bucket_t t; t.key = key; return prefix##_m_getp(h, &t); } \ SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_m_del(h, k); } \ @@ -287,6 +315,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26 KHASHL_INIT(KH_LOCAL, HType, prefix##_cs, HType##_cs_bucket_t, __kh_cached_hash, prefix##_cs_eq) \ SCOPE HType *prefix##_init(void) { return prefix##_cs_init(); } \ SCOPE void prefix##_destroy(HType *h) { prefix##_cs_destroy(h); } \ + SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_cs_save(h, fp); } \ + SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_cs_load(h, fp); } \ SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_cs_bucket_t t; t.key = key; t.hash = __hash_fn(key); return prefix##_cs_getp(h, &t); } \ SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_cs_del(h, k); } \ SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_cs_bucket_t t; t.key = key, t.hash = __hash_fn(key); return prefix##_cs_putp(h, &t, absent); } @@ -297,6 +327,8 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26 KHASHL_INIT(KH_LOCAL, HType, prefix##_cm, HType##_cm_bucket_t, __kh_cached_hash, prefix##_cm_eq) \ SCOPE HType *prefix##_init(void) { return prefix##_cm_init(); } \ SCOPE void prefix##_destroy(HType *h) { prefix##_cm_destroy(h); } \ + SCOPE khint_t prefix##_save(HType *h, FILE* fp) { return prefix##_cm_save(h, fp); } \ + SCOPE khint_t prefix##_load(HType **h, FILE* fp) { return prefix##_cm_load(h, fp); } \ SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_cm_bucket_t t; t.key = key; t.hash = __hash_fn(key); return prefix##_cm_getp(h, &t); } \ SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_cm_del(h, k); } \ SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_cm_bucket_t t; t.key = key, t.hash = __hash_fn(key); return prefix##_cm_putp(h, &t, absent); } diff --git a/sketch.cpp b/sketch.cpp index c3546ed..046bdf0 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -37,7 +37,11 @@ static inline int tq_shift(tiny_queue_t *q) * @param p minimizers */ 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) -{ +{ ///in default, w = 51, k = 51, is_hpc = 1 + /** + uint64_t x; + uint64_t rid:28, pos:27, rev:1, span:8; + **/ static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 }; uint64_t shift1 = k - 1, mask = (1ULL< 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63)); + ///sizeof(ha_mz1_t) = 16 memset(buf, 0xff, w * 16); memset(&tq, 0, sizeof(tiny_queue_t)); + ///len/w is the evaluated minimizer numbers kv_resize(ha_mz1_t, *p, p->n + len/w); for (i = l = buf_pos = min_pos = 0; i < len; ++i) { @@ -65,7 +71,11 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, tq_push(&tq, skip_len); kmer_span += skip_len; if (tq.count > k) kmer_span -= tq_shift(&tq); - } else kmer_span = l + 1 < k? l + 1 : k; + } else kmer_span = l + 1 < k? l + 1 : k; + ///kmer_span should be used for HPC k-mer + ///so for non-HPC k-mer, kmer_span should be k in any case? + ///kmer_span is used to calculate anchor pos on reverse complementary strand + kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer kmer[1] = (kmer[1] << 1 | (c>>1)) & mask; kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer @@ -80,6 +90,12 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span; } } else l = 0, tq.count = tq.front = 0, kmer_span = 0; + + + //for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i + //i is the real base iterator, while l is the HPC base iterator + //only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX) + //but even if l < k, infor is still stored into buf buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet for (j = buf_pos + 1; j < w; ++j) @@ -87,15 +103,26 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, for (j = 0; j < buf_pos; ++j) if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]); } + /** + * There are three cases: + * 1. info.x <= min.x, means info is a new minimizer + * 2. info.x > min.x, info is not a new minimizer + * (1) buf_pos != min_pos, do nothing + * (2) buf_pos == min_pos, means current minimizer has moved outside the window + * **/ + ///three cases: 1. if (info.x <= min.x) { // a new minimum; then write the old min if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); min = info, min_pos = buf_pos; } else if (buf_pos == min_pos) { // old min has moved outside the window if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); + ///buf_pos == min_pos, means current minimizer has moved outside the window + ///so for now we need to find a new minimizer at the current window (w k-mers) for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer for (j = 0; j <= buf_pos; ++j) if (min.x >= buf[j].x) min = buf[j], min_pos = j; + if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]); @@ -108,3 +135,121 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, if (min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); } + + + +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) +{ ///in default, w = 51, k = 51, is_hpc = 1 + /** + uint64_t x; + uint64_t rid:28, pos:27, rev:1, span:8; + **/ + static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 }; + uint64_t shift1 = k - 1, mask = (1ULL<a, (uint64_t)len); + k_flag->a.n = len; + memset(k_flag->a.a, 0, k_flag->a.n); + } + + + assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63)); + ///sizeof(ha_mz1_t) = 16 + memset(buf, 0xff, w * 16); + memset(&tq, 0, sizeof(tiny_queue_t)); + ///len/w is the evaluated minimizer numbers + kv_resize(ha_mz1_t, *p, p->n + len/w); + + for (i = l = buf_pos = min_pos = 0; i < len; ++i) { + int c = seq_nt4_table[(uint8_t)str[i]]; + ha_mz1_t info = dummy; + if (c < 4) { // not an ambiguous base + int z; + if (is_hpc) { + int skip_len = 1; + if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) { + for (skip_len = 2; i + skip_len < len; ++skip_len) + if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c) + break; + i += skip_len - 1; // put $i at the end of the current homopolymer run + } + tq_push(&tq, skip_len); + kmer_span += skip_len; + if (tq.count > k) kmer_span -= tq_shift(&tq); + } else kmer_span = l + 1 < k? l + 1 : k; + ///kmer_span should be used for HPC k-mer + ///so for non-HPC k-mer, kmer_span should be k in any case? + ///kmer_span is used to calculate anchor pos on reverse complementary strand + + if(k_flag != NULL) k_flag->a.a[i] = 1;///lable all useful base, which are not ignored by HPC + + kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer + kmer[1] = (kmer[1] << 1 | (c>>1)) & mask; + kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer + kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1; + if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand + z = kmer[1] < kmer[3]? 0 : 1; // strand + ++l; + if (l >= k && kmer_span < 256) { + uint64_t y; + y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]); + + filtered = 0; + if(hf != 0) filtered = ha_ft_isflt(hf, y); + ///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; + if(k_flag != NULL) k_flag->a.a[i]++; + if(k_flag != NULL && filtered == 1) k_flag->a.a[i]++; + } + } else l = 0, tq.count = tq.front = 0, kmer_span = 0; + + + //for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i + //i is the real base iterator, while l is the HPC base iterator + //only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX) + //but even if l < k, infor is still stored into buf + buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below + if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet + for (j = buf_pos + 1; j < w; ++j) + if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]); + for (j = 0; j < buf_pos; ++j) + if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]); + } + /** + * There are three cases: + * 1. info.x <= min.x, means info is a new minimizer + * 2. info.x > min.x, info is not a new minimizer + * (1) buf_pos != min_pos, do nothing + * (2) buf_pos == min_pos, means current minimizer has moved outside the window + * **/ + ///three cases: 1. + if (info.x <= min.x) { // a new minimum; then write the old min + if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); + min = info, min_pos = buf_pos; + } else if (buf_pos == min_pos) { // old min has moved outside the window + if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); + ///buf_pos == min_pos, means current minimizer has moved outside the window + ///so for now we need to find a new minimizer at the current window (w k-mers) + for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers + if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer + for (j = 0; j <= buf_pos; ++j) + if (min.x >= buf[j].x) min = buf[j], min_pos = j; + + if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers + for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted + if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]); + for (j = 0; j <= buf_pos; ++j) + if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]); + } + } + if (++buf_pos == w) buf_pos = 0; + } + if (min.x != UINT64_MAX) + kv_push(ha_mz1_t, *p, min); +} \ No newline at end of file