From 0ce5a7ba5bb1d2151f179cb21bab4bb12bdd6799 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Mon, 17 Jan 2022 18:29:48 -0500 Subject: [PATCH] integration/print het sites --- Assembly.cpp | 65 ++ Assembly.h | 1 + CommandLines.cpp | 6 + CommandLines.h | 4 +- Correct.cpp | 1876 +++++++++++++++++++++++++++++++++++++++++----- Correct.h | 76 +- Hash_Table.cpp | 194 +++-- Hash_Table.h | 26 +- Overlaps.cpp | 14 +- Overlaps.h | 23 +- Process_Read.cpp | 225 +++++- Process_Read.h | 7 +- anchor.cpp | 45 +- htab.cpp | 4 +- htab.h | 9 +- inter.cpp | 623 ++++++++++++++- kalloc.h | 30 + sketch.cpp | 41 +- 18 files changed, 2884 insertions(+), 385 deletions(-) diff --git a/Assembly.cpp b/Assembly.cpp index d012c5b..d20bf6f 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -11,6 +11,7 @@ #include "htab.h" #include "kthread.h" #include "rcut.h" +#include "kalloc.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, kvec_t_u64_warp* dbg_ct, st_mt_t *sp); @@ -21,6 +22,7 @@ void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list); All_reads R_INF; Debug_reads R_INF_FLAG; all_ul_t UL_INF; +uint32_t *het_cnt = NULL; void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length, char* new_read, int* new_length) { @@ -418,6 +420,36 @@ long long push_final_overlaps_increment(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev return available_overlaps; } +ha_ovec_buf_t *ha_ovec_buf_init(void *km, int is_final, int save_ov, int is_ug) +{ + ha_ovec_buf_t *b; + KCALLOC(km, b, 1); + b->is_final = !!is_final, b->save_ov = !!save_ov; + init_UC_Read(&b->self_read);//set 0 + init_UC_Read(&b->ovlp_read);//set 0 + init_Candidates_list(&b->clist);//set 0 + memset(&b->olist, 0, sizeof(overlap_region_alloc)); + memset(&b->olist_hp, 0, sizeof(overlap_region_alloc)); + // init_overlap_region_alloc(&b->olist); + // init_overlap_region_alloc(&b->olist_hp); + init_fake_cigar(&(b->tmp_region.f_cigar));//set 0 + kv_init(b->b_buf.a);//set 0 + kv_init(b->r_buf.a);//set 0 + kv_init(b->k_flag.a);//set 0 + kv_init(b->sp);//set 0 + if(!is_ug) b->ab = ha_abuf_init_buf(km); + else b->abl = ha_abufl_init_buf(km); + if (!b->is_final) { + init_Cigar_record_buf(&b->cigar1, km); + // init_Graph(&b->POA_Graph); + // init_Graph(&b->DAGCon); + init_Correct_dumy_buf(&b->correct, km);//set 0 + InitHaplotypeEvdience_buf(&b->hap, km); + init_Round2_alignment_buf(&b->round2, km); + } + return b; +} + ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov, int is_ug) { ha_ovec_buf_t *b; @@ -518,6 +550,17 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b) return mem; } +uint32_t get_het_cnt(haplotype_evdience_alloc *hap) +{ + uint32_t i, cnt; + for (i = cnt = 0; i < hap->snp_stat.n; i++) { + if(hap->snp_stat.a[i].score == 1 && (!(hap->snp_stat.a[i].occ_0 < 2 || hap->snp_stat.a[i].occ_1 < 2))) { + cnt++; + } + } + return cnt; +} + static void worker_ovec(void *data, long i, int tid) { ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid]; @@ -559,6 +602,8 @@ static void worker_ovec(void *data, long i, int tid) 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); } + + if(het_cnt) het_cnt[i] = get_het_cnt(&b->hap); } @@ -814,6 +859,20 @@ void rescue_hp_reads(ha_ovec_buf_t **b) } +void print_het_cnt_log(uint32_t *het_cnt) +{ + if(!het_cnt) return; + char* gfa_name = (char*)malloc(strlen(asm_opt.output_file_name)+35); + sprintf(gfa_name, "%s.het_cnt.log", asm_opt.output_file_name); + FILE* output_file = fopen(gfa_name, "w"); + fprintf(stderr, "[M::%s::] ==> print cnt of het sites to %s...\n", __func__, gfa_name); + free(gfa_name); + uint64_t i; + for (i = 0; i < R_INF.total_reads; i++){ + fprintf(output_file, ">%.*s\t%u\n", (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), het_cnt[i]); + } + fclose(output_file); +} void ha_overlap_and_correct(int round) @@ -838,6 +897,8 @@ void ha_overlap_and_correct(int round) ///debug_adapter(&asm_opt, &R_INF); if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated ha_opt_update_cov(&asm_opt, hom_cov); + het_cnt = NULL; + if(round == asm_opt.number_of_round-1 && asm_opt.is_dbg_het_cnt) CALLOC(het_cnt, R_INF.total_reads); if (asm_opt.required_read_name) kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads); else @@ -847,6 +908,10 @@ void ha_overlap_and_correct(int round) ha_pt_destroy(ha_idx); ha_idx = NULL; + if(het_cnt) { + print_het_cnt_log(het_cnt); free(het_cnt); het_cnt = NULL; + } + // collect statistics for (i = 0; i < asm_opt.thread_num; ++i) { asm_opt.num_bases += b[i]->num_read_base; diff --git a/Assembly.h b/Assembly.h index 7f37338..11df5f2 100644 --- a/Assembly.h +++ b/Assembly.h @@ -45,6 +45,7 @@ typedef struct { int ha_assemble(void); void ug_idx_build(ma_ug_t *ug, int hap_n); ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov, int is_ug); +ha_ovec_buf_t *ha_ovec_buf_init(void *km, int is_final, int save_ov, int is_ug); void ha_ovec_destroy(ha_ovec_buf_t *b); #endif diff --git a/CommandLines.cpp b/CommandLines.cpp index 6311edb..0211ff5 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -47,6 +47,8 @@ static ko_longopt_t long_options[] = { { "ul", ko_required_argument, 333}, { "unskew", ko_no_argument, 334}, { "kpt-rate", ko_required_argument, 335}, + { "ul-rate", ko_required_argument, 336}, + { "dbg-het-cnt", ko_no_argument, 337}, { 0, 0, 0 } }; @@ -230,6 +232,8 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->kpt_rate = -1; asm_opt->infor_cov = 3; asm_opt->s_hap_cov = 3; + asm_opt->ul_error_rate = 0.15; + asm_opt->is_dbg_het_cnt = 0; } void destory_enzyme(enzyme* f) @@ -745,6 +749,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 333) get_hic_enzymes(opt.arg, &(asm_opt->ar), 0); else if (c == 334) asm_opt->flag |= HA_F_USKEW; else if (c == 335) asm_opt->kpt_rate = atof(opt.arg); + else if (c == 336) asm_opt->ul_error_rate = atof(opt.arg); + else if (c == 337) asm_opt->is_dbg_het_cnt = 1; else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg); diff --git a/CommandLines.h b/CommandLines.h index 0dc877f..e0f14e6 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.16.2-r382" +#define HA_VERSION "0.16.2-r383" #define VERBOSE 0 @@ -122,6 +122,8 @@ typedef struct { int64_t hg_size; float kpt_rate; int64_t infor_cov, s_hap_cov; + double ul_error_rate; + uint8_t is_dbg_het_cnt; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Correct.cpp b/Correct.cpp index 1f9b9a9..a8c876e 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -6,8 +6,10 @@ #include "Levenshtein_distance.h" #include "Assembly.h" #include "CommandLines.h" -#include "ksw2.h" +// #include "ksw2.h" #include "ksort.h" +#include "kalloc.h" +#include "htab.h" #define generic_key(x) (x) KRADIX_SORT_INIT(b32, uint32_t, generic_key, 4) @@ -19,6 +21,9 @@ KRADIX_SORT_INIT(haplotype_evdience_srt, haplotype_evdience, haplotype_evdience_ #define haplotype_evdience_id_key(x) ((x).overlapID) KRADIX_SORT_INIT(haplotype_evdience_id_srt, haplotype_evdience, haplotype_evdience_id_key, member_size(haplotype_evdience, overlapID)) +#define overlap_region_dp_key(x) ((x).x_pos_e) +KRADIX_SORT_INIT(overlap_region_dp_srt, overlap_region, overlap_region_dp_key, member_size(overlap_region, x_pos_e)) + int ha_ov_type(const overlap_region *r, uint32_t len); @@ -40,6 +45,14 @@ void init_Round2_alignment(Round2_alignment* h) h->obtained_cigar_length = 0; } +void init_Round2_alignment_buf(Round2_alignment* h, void *km) +{ + init_Correct_dumy_buf(&(h->dumy), km); + init_Cigar_record_buf(&(h->cigar), km); + init_Cigar_record_buf(&(h->tmp_cigar), km); + h->obtained_cigar_length = 0; +} + void destory_Round2_alignment(Round2_alignment* h) { destory_Correct_dumy(&(h->dumy)); @@ -49,11 +62,11 @@ void destory_Round2_alignment(Round2_alignment* h) -inline int get_interval(long long window_start, long long window_end, overlap_region_alloc* overlap_list, Correct_dumy* dumy) +inline int get_interval(long long window_start, long long window_end, overlap_region_alloc* overlap_list, Correct_dumy* dumy, long long blockLen) { - uint64_t i; + uint64_t i, fud = 0; long long Len; - + if(window_start == 0) dumy->start_i = 0; for (i = dumy->start_i; i < overlap_list->length; i++) { ///this interval is smaller than all overlaps @@ -84,13 +97,15 @@ inline int get_interval(long long window_start, long long window_end, overlap_re dumy->length = 0; dumy->lengthNT = 0; + fud = 0; for (; i < overlap_list->length; i++) { if((Len = OVERLAP(window_start, window_end, (long long)overlap_list->list[i].x_pos_s, (long long)overlap_list->list[i].x_pos_e)) > 0) { ///sometimes the length of window > WINDOW, but overlap length == WINDOW - if (Len == WINDOW && window_end - window_start + 1 == WINDOW) + // if (Len == WINDOW && window_end - window_start + 1 == WINDOW) + if (Len == blockLen && window_end - window_start + 1 == blockLen) { dumy->overlapID[dumy->length] = i; dumy->length++; @@ -100,6 +115,7 @@ inline int get_interval(long long window_start, long long window_end, overlap_re dumy->lengthNT++; dumy->overlapID[dumy->size - dumy->lengthNT] = i; } + if(fud == 0) fud = 1, dumy->start_i = i; } if((long long)overlap_list->list[i].x_pos_s > window_end) @@ -121,10 +137,9 @@ inline int get_interval(long long window_start, long long window_end, overlap_re inline int get_available_interval(long long window_start, long long window_end, overlap_region_alloc* overlap_list, Correct_dumy* dumy) { - uint64_t i; + uint64_t i, fud = 0; long long Len; - - + if(window_start == 0) dumy->start_i = 0; for (i = dumy->start_i; i < overlap_list->length; i++) { ///this interval is smaller than all overlaps @@ -154,7 +169,7 @@ inline int get_available_interval(long long window_start, long long window_end, dumy->length = 0; dumy->lengthNT = 0; - + fud = 0; long long fake_length = 0; @@ -173,6 +188,7 @@ inline int get_available_interval(long long window_start, long long window_end, dumy->overlapID[dumy->length] = i; dumy->length++; } + if(fud == 0) fud = 1, dumy->start_i = i; } if((long long)overlap_list->list[i].x_pos_s > window_end) @@ -218,7 +234,15 @@ int extra_begin, int extra_end) memset(r+extra_begin+length, 'N', extra_end); } -int determine_overlap_region(int threshold, long long y_start, long long y_ID, long long Window_Len, All_reads* R_INF, +void fill_subregion_ul(char* r, long long start_pos, long long length, uint8_t strand, const ul_idx_t *uref, long long ID, +int extra_begin, int extra_end) +{ + retrieve_u_seq(NULL, r+extra_begin, &(uref->ug->u.a[ID]), strand, start_pos, length, NULL); + memset(r, 'N', extra_begin); + memset(r+extra_begin+length, 'N', extra_end); +} + +int determine_overlap_region(int threshold, long long y_start, long long y_ID, long long Window_Len, /**All_reads* R_INF**/long long y_len, int* r_extra_begin, int* r_extra_end, long long* r_y_start, long long* r_y_length) { int extra_begin; @@ -227,7 +251,8 @@ int* r_extra_begin, int* r_extra_end, long long* r_y_start, long long* r_y_lengt long long o_len; ///the length of y - currentIDLen = Get_READ_LENGTH((*R_INF), y_ID); + // currentIDLen = Get_READ_LENGTH((*R_INF), y_ID); + currentIDLen = y_len; ///since Window_Len == x_len + (threshold << 1) if(y_start < 0 || currentIDLen <= y_start || @@ -282,7 +307,7 @@ All_reads* R_INF) Window_Len = x_len + (threshold << 1); - if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF, + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, Get_READ_LENGTH((*R_INF), y_id), &extra_begin, &extra_end, &y_start, &o_len)) { return 0; @@ -350,11 +375,11 @@ char* r_string) /****************************may have bugs********************************/ - if(!determine_overlap_region(THRESHOLD, y_start, overlap_list->list[currentID].y_id, Window_Len, R_INF, + if(!determine_overlap_region(THRESHOLD, y_start, overlap_list->list[currentID].y_id, Window_Len, Get_READ_LENGTH((*R_INF), overlap_list->list[currentID].y_id), &extra_begin, &extra_end, &y_start, &o_len)) { append_window_list(&overlap_list->list[currentID], window_start, window_end, - -1, -1, -1, -1, -1, -1); + -1, -1, -1, -1, -1, -1, WINDOW, NULL); continue; } @@ -384,12 +409,12 @@ char* r_string) append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], y_startGroup[0] + return_sites[0], (int)return_sites_error[0], - y_extra_begin[0], y_extra_end[0], error_threshold[0]); + y_extra_begin[0], y_extra_end[0], error_threshold[0], WINDOW, NULL); } else { append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], -1, -1, - y_extra_begin[0], y_extra_end[0], error_threshold[0]); + y_extra_begin[0], y_extra_end[0], error_threshold[0], WINDOW, NULL); } @@ -399,12 +424,12 @@ char* r_string) append_window_list(&overlap_list->list[overlapID[1]], window_start, window_end, y_startGroup[1], y_startGroup[1] + return_sites[1], (int)return_sites_error[1], - y_extra_begin[1], y_extra_end[1], error_threshold[1]); + y_extra_begin[1], y_extra_end[1], error_threshold[1], WINDOW, NULL); } else { append_window_list(&overlap_list->list[overlapID[1]], window_start, window_end, y_startGroup[1], -1, -1, - y_extra_begin[1], y_extra_end[1], error_threshold[1]); + y_extra_begin[1], y_extra_end[1], error_threshold[1], WINDOW, NULL); } @@ -414,12 +439,12 @@ char* r_string) append_window_list(&overlap_list->list[overlapID[2]], window_start, window_end, y_startGroup[2], y_startGroup[2] + return_sites[2], (int)return_sites_error[2], - y_extra_begin[2], y_extra_end[2], error_threshold[2]); + y_extra_begin[2], y_extra_end[2], error_threshold[2], WINDOW, NULL); } else { append_window_list(&overlap_list->list[overlapID[2]], window_start, window_end, y_startGroup[2], -1, -1, - y_extra_begin[2], y_extra_end[2], error_threshold[2]); + y_extra_begin[2], y_extra_end[2], error_threshold[2], WINDOW, NULL); } @@ -429,12 +454,12 @@ char* r_string) append_window_list(&overlap_list->list[overlapID[3]], window_start, window_end, y_startGroup[3], y_startGroup[3] + return_sites[3], (int)return_sites_error[3], - y_extra_begin[3], y_extra_end[3], error_threshold[3]); + y_extra_begin[3], y_extra_end[3], error_threshold[3], WINDOW, NULL); } else { append_window_list(&overlap_list->list[overlapID[3]], window_start, window_end, y_startGroup[3], -1, -1, - y_extra_begin[3], y_extra_end[3], error_threshold[3]); + y_extra_begin[3], y_extra_end[3], error_threshold[3], WINDOW, NULL); } } } @@ -449,12 +474,12 @@ char* r_string) append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], y_startGroup[0] + end_site, (int)error, - y_extra_begin[0], y_extra_end[0], error_threshold[0]); + y_extra_begin[0], y_extra_end[0], error_threshold[0], WINDOW, NULL); } else { append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], -1, -1, - y_extra_begin[0], y_extra_end[0], error_threshold[0]); + y_extra_begin[0], y_extra_end[0], error_threshold[0], WINDOW, NULL); } } else if (groupLen > 1) @@ -470,12 +495,12 @@ char* r_string) overlap_list->list[overlapID[i]].align_length += x_len; append_window_list(&overlap_list->list[overlapID[i]], window_start, window_end, y_startGroup[i], y_startGroup[i] + return_sites[i], (int)return_sites_error[i], - y_extra_begin[i], y_extra_end[i], error_threshold[i]); + y_extra_begin[i], y_extra_end[i], error_threshold[i], WINDOW, NULL); } else { append_window_list(&overlap_list->list[overlapID[i]], window_start, window_end, y_startGroup[i], -1, -1, - y_extra_begin[i], y_extra_end[i], error_threshold[i]); + y_extra_begin[i], y_extra_end[i], error_threshold[i], WINDOW, NULL); } } @@ -509,11 +534,11 @@ char* r_string) Window_Len = x_len + (threshold << 1); - if(!determine_overlap_region(threshold, y_start, overlap_list->list[currentID].y_id, Window_Len, R_INF, + if(!determine_overlap_region(threshold, y_start, overlap_list->list[currentID].y_id, Window_Len, Get_READ_LENGTH((*R_INF), overlap_list->list[currentID].y_id), &extra_begin, &extra_end, &y_start, &o_len)) { append_window_list(&overlap_list->list[currentID], x_start, x_end, - -1, -1, -1, -1, -1, -1); + -1, -1, -1, -1, -1, -1, WINDOW, NULL); continue; } @@ -529,24 +554,24 @@ char* r_string) { overlap_list->list[currentID].align_length += x_len; append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, y_start + end_site, (int)error, - extra_begin, extra_end, threshold); + extra_begin, extra_end, threshold, WINDOW, NULL); } else { append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, -1, -1, - extra_begin, extra_end, threshold); + extra_begin, extra_end, threshold, WINDOW, NULL); } } } -void verify_ul_window(long long window_start, long long window_end, overlap_region_alloc* overlap_list,Correct_dumy* dumy, const ma_ug_t *uref, -char* r_string) +void verify_ul_window(long long window_start, long long window_end, overlap_region_alloc* overlap_list,Correct_dumy* dumy, const ul_idx_t *uref, +char* r_string, double max_ov_diff_ec, long long blockLen, long long max_error, void *km) { long long i; long long currentID; long long x_start, y_start, o_len; - long long Window_Len = WINDOW + (THRESHOLD << 1); + long long Window_Len = blockLen + (max_error << 1); char* x_string = NULL; char* y_string = NULL; long long x_end, x_len; @@ -554,7 +579,7 @@ char* r_string) unsigned int error; int groupLen = 0; int return_sites[GROUP_SIZE]; - unsigned int return_sites_error[GROUP_SIZE]; + unsigned int return_sites_error[GROUP_SIZE]; uint64_t overlapID[GROUP_SIZE]; uint64_t y_startGroup[GROUP_SIZE]; int y_extra_begin[GROUP_SIZE]; @@ -563,12 +588,12 @@ char* r_string) int extra_begin; int extra_end; - ///here are overlaps fully covered by WINDOW + ///here are overlaps fully covered by blockLen for (i = 0; i < (long long)dumy->length; i++) { extra_begin = extra_end = 0; ///if the window has been fully covered, the interval at x is [window_start, window_end] - x_len = WINDOW; + x_len = blockLen; currentID = dumy->overlapID[i]; x_start = window_start; ///offset of y @@ -578,33 +603,45 @@ char* r_string) /****************************may have bugs********************************/ - if(!determine_overlap_region(THRESHOLD, y_start, overlap_list->list[currentID].y_id, Window_Len, R_INF, + if(!determine_overlap_region(max_error, y_start, overlap_list->list[currentID].y_id, Window_Len, uref->ug->u.a[overlap_list->list[currentID].y_id].len, &extra_begin, &extra_end, &y_start, &o_len)) { append_window_list(&overlap_list->list[currentID], window_start, window_end, - -1, -1, -1, -1, -1, -1); + -1, -1, -1, -1, -1, -1, blockLen, km); continue; } - fill_subregion(dumy->overlap_region_group[groupLen], y_start, o_len, overlap_list->list[currentID].y_pos_strand, - R_INF, overlap_list->list[currentID].y_id, extra_begin, extra_end); + fill_subregion_ul(dumy->overlap_region_group[groupLen], y_start, o_len, overlap_list->list[currentID].y_pos_strand, + uref, overlap_list->list[currentID].y_id, extra_begin, extra_end); y_extra_begin[groupLen] = extra_begin; y_extra_end[groupLen] = extra_end; overlapID[groupLen] = currentID; y_startGroup[groupLen] = y_start; - error_threshold[groupLen] = THRESHOLD; + error_threshold[groupLen] = max_error; x_string = r_string + x_start; groupLen++; if (groupLen == GROUP_SIZE) { - Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], - dumy->overlap_region_group[2], dumy->overlap_region_group[3], Window_Len, x_string, WINDOW, - return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE); - groupLen = 0; + // Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], + // dumy->overlap_region_group[2], dumy->overlap_region_group[3], Window_Len, x_string, blockLen, + // return_sites, return_sites_error, max_error, dumy->Peq_SSE); + return_sites[0] = + Reserve_Banded_BPM(dumy->overlap_region_group[0], Window_Len, x_string, blockLen, max_error, &return_sites_error[0]); + return_sites[1] = + Reserve_Banded_BPM(dumy->overlap_region_group[1], Window_Len, x_string, blockLen, max_error, &return_sites_error[1]); + return_sites[2] = + Reserve_Banded_BPM(dumy->overlap_region_group[2], Window_Len, x_string, blockLen, max_error, &return_sites_error[2]); + return_sites[3] = + Reserve_Banded_BPM(dumy->overlap_region_group[3], Window_Len, x_string, blockLen, max_error, &return_sites_error[3]); + + + + + groupLen = 0; if (return_sites_error[0]!=(unsigned int)-1) { @@ -612,12 +649,12 @@ char* r_string) append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], y_startGroup[0] + return_sites[0], (int)return_sites_error[0], - y_extra_begin[0], y_extra_end[0], error_threshold[0]); + y_extra_begin[0], y_extra_end[0], error_threshold[0], blockLen, km); } else { append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], -1, -1, - y_extra_begin[0], y_extra_end[0], error_threshold[0]); + y_extra_begin[0], y_extra_end[0], error_threshold[0], blockLen, km); } @@ -627,12 +664,12 @@ char* r_string) append_window_list(&overlap_list->list[overlapID[1]], window_start, window_end, y_startGroup[1], y_startGroup[1] + return_sites[1], (int)return_sites_error[1], - y_extra_begin[1], y_extra_end[1], error_threshold[1]); + y_extra_begin[1], y_extra_end[1], error_threshold[1], blockLen, km); } else { append_window_list(&overlap_list->list[overlapID[1]], window_start, window_end, y_startGroup[1], -1, -1, - y_extra_begin[1], y_extra_end[1], error_threshold[1]); + y_extra_begin[1], y_extra_end[1], error_threshold[1], blockLen, km); } @@ -642,12 +679,12 @@ char* r_string) append_window_list(&overlap_list->list[overlapID[2]], window_start, window_end, y_startGroup[2], y_startGroup[2] + return_sites[2], (int)return_sites_error[2], - y_extra_begin[2], y_extra_end[2], error_threshold[2]); + y_extra_begin[2], y_extra_end[2], error_threshold[2], blockLen, km); } else { append_window_list(&overlap_list->list[overlapID[2]], window_start, window_end, y_startGroup[2], -1, -1, - y_extra_begin[2], y_extra_end[2], error_threshold[2]); + y_extra_begin[2], y_extra_end[2], error_threshold[2], blockLen, km); } @@ -657,19 +694,19 @@ char* r_string) append_window_list(&overlap_list->list[overlapID[3]], window_start, window_end, y_startGroup[3], y_startGroup[3] + return_sites[3], (int)return_sites_error[3], - y_extra_begin[3], y_extra_end[3], error_threshold[3]); + y_extra_begin[3], y_extra_end[3], error_threshold[3], blockLen, km); } else { append_window_list(&overlap_list->list[overlapID[3]], window_start, window_end, y_startGroup[3], -1, -1, - y_extra_begin[3], y_extra_end[3], error_threshold[3]); + y_extra_begin[3], y_extra_end[3], error_threshold[3], blockLen, km); } } } if (groupLen == 1) { - end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], Window_Len, x_string, WINDOW, THRESHOLD, &error); + end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], Window_Len, x_string, blockLen, max_error, &error); if (error!=(unsigned int)-1) { @@ -677,33 +714,35 @@ char* r_string) append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], y_startGroup[0] + end_site, (int)error, - y_extra_begin[0], y_extra_end[0], error_threshold[0]); + y_extra_begin[0], y_extra_end[0], error_threshold[0], blockLen, km); } else { append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], -1, -1, - y_extra_begin[0], y_extra_end[0], error_threshold[0]); + y_extra_begin[0], y_extra_end[0], error_threshold[0], blockLen, km); } } else if (groupLen > 1) { - Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], - dumy->overlap_region_group[2], dumy->overlap_region_group[3], Window_Len, x_string, WINDOW, - return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE); + // Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], + // dumy->overlap_region_group[2], dumy->overlap_region_group[3], Window_Len, x_string, blockLen, + // return_sites, return_sites_error, max_error, dumy->Peq_SSE); for (i = 0; i < groupLen; i++) { + return_sites[i] = Reserve_Banded_BPM(dumy->overlap_region_group[i], Window_Len, x_string, blockLen, max_error, &return_sites_error[i]); + if (return_sites_error[i]!=(unsigned int)-1) { overlap_list->list[overlapID[i]].align_length += x_len; append_window_list(&overlap_list->list[overlapID[i]], window_start, window_end, y_startGroup[i], y_startGroup[i] + return_sites[i], (int)return_sites_error[i], - y_extra_begin[i], y_extra_end[i], error_threshold[i]); + y_extra_begin[i], y_extra_end[i], error_threshold[i], blockLen, km); } else { append_window_list(&overlap_list->list[overlapID[i]], window_start, window_end, y_startGroup[i], -1, -1, - y_extra_begin[i], y_extra_end[i], error_threshold[i]); + y_extra_begin[i], y_extra_end[i], error_threshold[i], blockLen, km); } } @@ -714,7 +753,7 @@ char* r_string) long long reverse_i = dumy->size - 1; int threshold; - ///here are overlaps partially covered by WINDOW + ///here are overlaps partially covered by blockLen for (i = 0; i < (long long)dumy->lengthNT; i++) { extra_begin = extra_end = 0; @@ -724,9 +763,10 @@ char* r_string) ///overlap length between [window_start, window_end] x_len = x_end - x_start + 1; - threshold = x_len * asm_opt.max_ov_diff_ec; + threshold = x_len * max_ov_diff_ec; /****************************may have bugs********************************/ threshold = Adjust_Threshold(threshold, x_len); + if(threshold > THRESHOLD_MAX_SIZE) threshold = THRESHOLD_MAX_SIZE; /****************************may have bugs********************************/ ///offset of y @@ -737,16 +777,16 @@ char* r_string) Window_Len = x_len + (threshold << 1); - if(!determine_overlap_region(threshold, y_start, overlap_list->list[currentID].y_id, Window_Len, R_INF, + if(!determine_overlap_region(threshold, y_start, overlap_list->list[currentID].y_id, Window_Len, uref->ug->u.a[overlap_list->list[currentID].y_id].len, &extra_begin, &extra_end, &y_start, &o_len)) { append_window_list(&overlap_list->list[currentID], x_start, x_end, - -1, -1, -1, -1, -1, -1); + -1, -1, -1, -1, -1, -1, blockLen, km); continue; } - fill_subregion(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, - R_INF, overlap_list->list[currentID].y_id, extra_begin, extra_end); + fill_subregion_ul(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, + uref, overlap_list->list[currentID].y_id, extra_begin, extra_end); x_string = r_string + x_start; y_string = dumy->overlap_region; @@ -757,12 +797,121 @@ char* r_string) { overlap_list->list[currentID].align_length += x_len; append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, y_start + end_site, (int)error, - extra_begin, extra_end, threshold); + extra_begin, extra_end, threshold, blockLen, km); } else { append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, -1, -1, - extra_begin, extra_end, threshold); + extra_begin, extra_end, threshold, blockLen, km); + } + } +} + +void verify_ul_ll_window(long long window_start, long long window_end, overlap_region_alloc* overlap_list,Correct_dumy* dumy, const ul_idx_t *uref, +char* r_string, double max_ov_diff_ec, long long blockLen, void *km) +{ + long long i, reverse_i = dumy->size - 1; + long long currentID, threshold = /**THRESHOLD_MAX_SIZE**/THRESHOLD; + long long x_start, y_start, o_len; + long long Window_Len = blockLen + (threshold << 1); + char* x_string = NULL; + char* y_string = NULL; + long long x_end, x_len; + int end_site; + unsigned int error; + int extra_begin; + int extra_end; + + ///here are overlaps fully covered by blockLen + for (i = 0; i < (long long)dumy->length; i++) + { + extra_begin = extra_end = 0; + ///if the window has been fully covered, the interval at x is [window_start, window_end] + x_len = blockLen; + currentID = dumy->overlapID[i]; + x_start = window_start; + ///offset of y + y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s; + /****************************may have bugs********************************/ + y_start += y_start_offset(x_start, &(overlap_list->list[currentID].f_cigar)); + /****************************may have bugs********************************/ + + + if(!determine_overlap_region(threshold, y_start, overlap_list->list[currentID].y_id, Window_Len, uref->ug->u.a[overlap_list->list[currentID].y_id].len, + &extra_begin, &extra_end, &y_start, &o_len)) + { + append_window_list(&overlap_list->list[currentID], window_start, window_end, + -1, -1, -1, -1, -1, -1, blockLen, km); + continue; + } + + fill_subregion_ul(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, + uref, overlap_list->list[currentID].y_id, extra_begin, extra_end); + + end_site = Reserve_Banded_BPM(dumy->overlap_region, Window_Len, x_string, blockLen, threshold, &error); + if (error!=(unsigned int)-1) + { + overlap_list->list[currentID].align_length += x_len; + append_window_list(&overlap_list->list[currentID], window_start, window_end, y_start, y_start + end_site, (int)error, + extra_begin, extra_end, threshold, blockLen, km); + } + else + { + append_window_list(&overlap_list->list[currentID], window_start, window_end, y_start, -1, -1, + extra_begin, extra_end, threshold, blockLen, km); + } + } + + ///here are overlaps partially covered by blockLen + for (i = 0; i < (long long)dumy->lengthNT; i++) + { + extra_begin = extra_end = 0; + currentID = dumy->overlapID[reverse_i--]; + x_start = MAX(window_start, (long long)overlap_list->list[currentID].x_pos_s); + x_end = MIN(window_end, (long long)overlap_list->list[currentID].x_pos_e); + + ///overlap length between [window_start, window_end] + x_len = x_end - x_start + 1; + threshold = x_len * max_ov_diff_ec; + /****************************may have bugs********************************/ + threshold = Adjust_Threshold(threshold, x_len); + if(threshold > THRESHOLD_MAX_SIZE) threshold = THRESHOLD_MAX_SIZE; + /****************************may have bugs********************************/ + + ///offset of y + y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s; + /****************************may have bugs********************************/ + y_start += y_start_offset(x_start, &(overlap_list->list[currentID].f_cigar)); + /****************************may have bugs********************************/ + + Window_Len = x_len + (threshold << 1); + + if(!determine_overlap_region(threshold, y_start, overlap_list->list[currentID].y_id, Window_Len, uref->ug->u.a[overlap_list->list[currentID].y_id].len, + &extra_begin, &extra_end, &y_start, &o_len)) + { + append_window_list(&overlap_list->list[currentID], x_start, x_end, + -1, -1, -1, -1, -1, -1, blockLen, km); + continue; + } + + fill_subregion_ul(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, + uref, overlap_list->list[currentID].y_id, extra_begin, extra_end); + + x_string = r_string + x_start; + y_string = dumy->overlap_region; + + end_site = Reserve_Banded_BPM(y_string, Window_Len, x_string, x_len, threshold, &error); + + if (error!=(unsigned int)-1) + { + overlap_list->list[currentID].align_length += x_len; + append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, y_start + end_site, (int)error, + extra_begin, extra_end, threshold, blockLen, km); + } + else + { + append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, -1, -1, + extra_begin, extra_end, threshold, blockLen, km); } } } @@ -860,7 +1009,7 @@ UC_Read* g_read) SubWindowLen = SubLen + (SubThreshold << 1); - if(determine_overlap_region(SubThreshold, y_start, y_id, SubWindowLen, R_INF, + if(determine_overlap_region(SubThreshold, y_start, y_id, SubWindowLen, Get_READ_LENGTH((*R_INF), y_id), &extra_begin, &extra_end, &y_start, &o_len) == 0) { T_error = T_error + (x_end - x_start + 1) * error_rate * 1.5; @@ -915,6 +1064,16 @@ inline int double_error_threshold(int pre_threshold, int x_len) return threshold; } +inline int double_ul_error_threshold(int pre_threshold, int x_len) +{ + + pre_threshold = Adjust_Threshold(pre_threshold, x_len); + int threshold = THRESHOLD_UL_MAX * x_len; + if(threshold < pre_threshold) threshold = pre_threshold; + if(threshold > THRESHOLD_MAX_SIZE) threshold = THRESHOLD_MAX_SIZE; + return threshold; +} + inline int verify_sub_window(All_reads* R_INF, Correct_dumy* dumy, UC_Read* g_read, long long x_beg, long long xLen, long long y_beg, long long yLen, uint64_t y_id, @@ -930,7 +1089,7 @@ unsigned int* get_error, int* get_y_end, int* get_x_end, int* get_aligned_xLen) long long o_len; unsigned int r_error; if(!determine_overlap_region(threshold, y_beg, y_id, yLen, - R_INF, &extra_begin, &extra_end, &y_beg, &o_len)) + Get_READ_LENGTH((*R_INF), y_id), &extra_begin, &extra_end, &y_beg, &o_len)) { return 0; } @@ -961,6 +1120,51 @@ unsigned int* get_error, int* get_y_end, int* get_x_end, int* get_aligned_xLen) } } +inline int verify_ul_sub_window(const ul_idx_t *uref, Correct_dumy* dumy, UC_Read* g_read, +long long x_beg, long long xLen, long long y_beg, long long yLen, uint64_t y_id, +uint64_t y_pos_strand, int threshold, int alignment_strand, +unsigned int* get_error, int* get_y_end, int* get_x_end, int* get_aligned_xLen) +{ + (*get_aligned_xLen) = 0; + (*get_y_end) = -1; + (*get_x_end) = -1; + (*get_error) = (unsigned int)-1; + + int extra_begin, extra_end, r_x_end, r_y_end, aligned_xLen; + long long o_len; + unsigned int r_error; + if(!determine_overlap_region(threshold, y_beg, y_id, yLen, + uref->ug->u.a[y_id].len, &extra_begin, &extra_end, &y_beg, &o_len)) + { + return 0; + } + + fill_subregion_ul(dumy->overlap_region, y_beg, o_len, y_pos_strand, uref, + y_id, extra_begin, extra_end); + + char* x_string = g_read->seq + x_beg; + char* y_string = dumy->overlap_region; + + aligned_xLen = 0; + + alignment_extension(y_string, yLen, x_string, xLen, threshold, + alignment_strand, &r_error, &r_y_end, &r_x_end, &aligned_xLen); + + (*get_error) = r_error; + (*get_y_end) = r_y_end; + (*get_x_end) = r_x_end; + (*get_aligned_xLen) = aligned_xLen; + + if(aligned_xLen == 0) + { + return 0; + } + else + { + return 1; + } +} + inline double non_trim_error_rate(overlap_region_alloc* overlap_list, long long ID, All_reads* R_INF, Correct_dumy* dumy, UC_Read* g_read) { @@ -1084,6 +1288,130 @@ All_reads* R_INF, Correct_dumy* dumy, UC_Read* g_read) return error_rate; } + +inline double non_trim_ul_error_rate(overlap_region_alloc* overlap_list, long long ID, +const ul_idx_t *uref, Correct_dumy* dumy, UC_Read* g_read) +{ + long long tLen, tError,i, subWinLen, subWinNum; + + tLen = 0; + tError = 0; + + subWinNum = overlap_list->list[ID].w_list_length; + + + for (i = 0; i < subWinNum; i++) + { + subWinLen = overlap_list->list[ID].w_list[i].x_end - overlap_list->list[ID].w_list[i].x_start + 1; + tLen += subWinLen; + + if(overlap_list->list[ID].w_list[i].y_end != -1) + { + tError += overlap_list->list[ID].w_list[i].error; + } + else + { + int x_len = subWinLen; + int threshold = double_ul_error_threshold(overlap_list->list[ID].w_list[i].error_threshold, x_len); + int Window_Len = x_len + (threshold << 1); + unsigned int r_error_left = 0; + int r_x_end_left, r_y_end_left, aligned_xLen_left; + unsigned int r_error_right = 0; + int r_x_end_right, r_y_end_right, aligned_xLen_right; + long long y_beg_left, y_beg_right; + + aligned_xLen_left = aligned_xLen_right = 0; + y_beg_left = y_beg_right = -1; + + if(overlap_list->list[ID].w_list[i].y_start == -1) + { + tError += x_len; + continue; + } + + ///if the previous window is mapped + if(i > 0 && overlap_list->list[ID].w_list[i - 1].y_end != -1) + { + y_beg_left = overlap_list->list[ID].w_list[i - 1].y_end + 1; + } + + ///if the next window is mapped + if(i < (long long)(overlap_list->list[ID].w_list_length - 1) && overlap_list->list[ID].w_list[i + 1].y_end != -1) + { + y_beg_right = 1 + overlap_list->list[ID].w_list[i + 1].y_start - 1 - x_len; + } + + + if(y_beg_left == -1 && y_beg_right == -1) + { + y_beg_left = overlap_list->list[ID].w_list[i].y_start; + if(overlap_list->list[ID].w_list[i].extra_begin >= 0) + { + y_beg_left = y_beg_left + overlap_list->list[ID].w_list[i].error_threshold - + overlap_list->list[ID].w_list[i].extra_begin; + } + y_beg_right = y_beg_left; + } + + if(y_beg_left == -1 && y_beg_right != -1) + { + y_beg_left = y_beg_right; + } + + if(y_beg_right == -1 && y_beg_left != -1) + { + y_beg_right = y_beg_left; + } + + + if(y_beg_left != -1) + { + verify_ul_sub_window(uref, dumy, g_read, overlap_list->list[ID].w_list[i].x_start, + x_len, y_beg_left, Window_Len, overlap_list->list[ID].y_id, overlap_list->list[ID].y_pos_strand, + threshold, 0, &r_error_left, &r_y_end_left, &r_x_end_left, &aligned_xLen_left); + } + + if(y_beg_right != -1) + { + verify_ul_sub_window(uref, dumy, g_read, overlap_list->list[ID].w_list[i].x_start, + x_len, y_beg_right, Window_Len, overlap_list->list[ID].y_id, overlap_list->list[ID].y_pos_strand, + threshold, 1, &r_error_right, &r_y_end_right, &r_x_end_right, &aligned_xLen_right); + } + + ///aligned in both direction + if(aligned_xLen_left != 0 && aligned_xLen_right != 0) + { + if(aligned_xLen_left + aligned_xLen_right <= x_len) + { + tError = tError + r_error_left + r_error_right + + (x_len - aligned_xLen_left - aligned_xLen_right); + } + else + { + float E_rate = (float)(x_len)/(float)(aligned_xLen_left + aligned_xLen_right); + tError = tError + (r_error_left + r_error_right)*E_rate; + } + }///not aligned in both direction + else if(aligned_xLen_left == 0 && aligned_xLen_right == 0) + { + tError += x_len; + }///only aligned in left + else if(aligned_xLen_left != 0) + { + tError = tError + r_error_left + (x_len - aligned_xLen_left); + }///only aligned in right + else if(aligned_xLen_right != 0) + { + tError = tError + r_error_right + (x_len - aligned_xLen_right); + } + } + } + + double error_rate = (double)(tError)/(double)(tLen); + + return error_rate; +} + int calculate_hpm_errors(char* x, int x_len, char* y, int y_len, CIGAR* cigar, int error) { int x_i, y_i, cigar_i; @@ -1912,6 +2240,126 @@ int verify_cigar_2(char* x, int x_len, char* y, int y_len, Cigar_record* cigar, } +inline int fix_ul_boundary(char* x_string, long long x_len, int threshold, +long long total_y_start, long long local_y_start, long long local_y_end, +long long old_extra_begin, long long old_extra_end, +long long y_ID, long long Window_Len, const ul_idx_t *uref, +Correct_dumy* dumy, int y_strand, unsigned int old_error, +long long* r_total_y_start, int* r_start_site, int* r_end_site, +int* r_extra_begin, int* r_extra_end, unsigned int* r_error) +{ + + + int new_extra_begin, new_extra_end; + long long new_y_start, new_y_length; + int new_end_site, new_start_site; + unsigned int new_error; + char* y_string; + + + int path_length; + + ///if the start pos at the left boundary + if(local_y_start == 0) + { + total_y_start = total_y_start + local_y_start; + ///if local_y_start == 0 and old_extra_begin != 0 + ///this means total_y_start == 0, so shift to the left cannot get a new start pos + if(old_extra_begin != 0) + { + return 0; + } + + ///if the begining of alignment is 0, we should try to shift the window to find a better result + ///shift to the left by threshold-1 bases + if(!determine_overlap_region(threshold, total_y_start, y_ID, Window_Len, uref->ug->u.a[y_ID].len, + &new_extra_begin, &new_extra_end, &new_y_start, &new_y_length)) + { + return 0; + } + + ///if new_y_start is equal to total_y_start, recalculate makes no sense + if(new_y_start == total_y_start) + { + return 0; + } + + fill_subregion_ul(dumy->overlap_region_fix, new_y_start, new_y_length, y_strand, uref, y_ID, + new_extra_begin, new_extra_end); + + y_string = dumy->overlap_region_fix; + + new_end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &new_error, &new_start_site, + &path_length, dumy->matrix_bit, dumy->path_fix, -1, -1); + + if (new_error != (unsigned int)-1 && new_error < old_error) + { + (*r_total_y_start) = new_y_start; + (*r_start_site) = new_start_site; + (*r_end_site) = new_end_site; + (*r_extra_begin) = new_extra_begin; + (*r_extra_end) = new_extra_end; + (*r_error) = new_error; + + dumy->path_length = path_length; + memcpy(dumy->path, dumy->path_fix, path_length); + memcpy(dumy->overlap_region, dumy->overlap_region_fix, Window_Len); + return 1; + } + } + else if(local_y_end == Window_Len - 1) + { + ///if local_y_end == Window_Len - 1 and old_extra_end > 0 + ///this means local_y_end is the end of the y + ///so shit to the right makes no sense + if(old_extra_end != 0) + { + return 0; + } + long long total_y_end = total_y_start + local_y_end; + + total_y_start = total_y_end - x_len + 1; + + if(!determine_overlap_region(threshold, total_y_start, y_ID, Window_Len, uref->ug->u.a[y_ID].len, + &new_extra_begin, &new_extra_end, &new_y_start, &new_y_length)) + { + return 0; + } + + if(new_y_start == total_y_end - local_y_end) + { + return 0; + } + + fill_subregion_ul(dumy->overlap_region_fix, new_y_start, new_y_length, y_strand, uref, y_ID, + new_extra_begin, new_extra_end); + + y_string = dumy->overlap_region_fix; + + new_end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &new_error, &new_start_site, + &path_length, dumy->matrix_bit, dumy->path_fix, -1, -1); + + if (new_error != (unsigned int)-1 && new_error < old_error) + { + (*r_total_y_start) = new_y_start; + (*r_start_site) = new_start_site; + (*r_end_site) = new_end_site; + (*r_extra_begin) = new_extra_begin; + (*r_extra_end) = new_extra_end; + (*r_error) = new_error; + + dumy->path_length = path_length; + memcpy(dumy->path, dumy->path_fix, path_length); + memcpy(dumy->overlap_region, dumy->overlap_region_fix, Window_Len); + return 1; + } + + + } + return 0; +} + + inline int fix_boundary(char* x_string, long long x_len, int threshold, long long total_y_start, long long local_y_start, long long local_y_end, long long old_extra_begin, long long old_extra_end, @@ -1944,7 +2392,7 @@ int* r_extra_begin, int* r_extra_end, unsigned int* r_error) ///if the begining of alignment is 0, we should try to shift the window to find a better result ///shift to the left by threshold-1 bases - if(!determine_overlap_region(threshold, total_y_start, y_ID, Window_Len, R_INF, + if(!determine_overlap_region(threshold, total_y_start, y_ID, Window_Len, Get_READ_LENGTH((*R_INF), y_ID), &new_extra_begin, &new_extra_end, &new_y_start, &new_y_length)) { return 0; @@ -1992,7 +2440,7 @@ int* r_extra_begin, int* r_extra_end, unsigned int* r_error) total_y_start = total_y_end - x_len + 1; - if(!determine_overlap_region(threshold, total_y_start, y_ID, Window_Len, R_INF, + if(!determine_overlap_region(threshold, total_y_start, y_ID, Window_Len, Get_READ_LENGTH((*R_INF), y_ID), &new_extra_begin, &new_extra_end, &new_y_start, &new_y_length)) { return 0; @@ -2102,7 +2550,7 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN y_start = total_y_start; Window_Len = x_len + (threshold << 1); - if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF, + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, Get_READ_LENGTH((*R_INF), y_id), &extra_begin, &extra_end, &y_start, &o_len)) { break; @@ -2266,7 +2714,7 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN ///y_start might be less than 0 y_start = total_y_end - x_len + 1; - if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF, + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, Get_READ_LENGTH((*R_INF), y_id), &extra_begin, &extra_end, &y_start, &o_len)) { break; @@ -2549,7 +2997,7 @@ void debug_scan_cigar(overlap_region* sub_list) void calculate_boundary_cigars(overlap_region* sub_list, All_reads* R_INF, Correct_dumy* dumy, UC_Read* g_read) { - resize_window_list_alloc(&(sub_list->boundary_cigars), sub_list->w_list_length - 1); + resize_window_list_alloc(&(sub_list->boundary_cigars), sub_list->w_list_length - 1, NULL); int y_id = sub_list->y_id; int y_strand = sub_list->y_pos_strand; long long y_readLen = Get_READ_LENGTH((*R_INF), y_id); @@ -2663,7 +3111,7 @@ UC_Read* g_read) threshold = double_error_threshold(threshold, xLen); yLen = xLen + (threshold << 1); - if(!determine_overlap_region(threshold, y_start, y_id, yLen, R_INF, + if(!determine_overlap_region(threshold, y_start, y_id, yLen, Get_READ_LENGTH((*R_INF), y_id), &extra_begin, &extra_end, &y_start, &o_len)) { sub_list->boundary_cigars.buffer[i].error = -1; @@ -2769,6 +3217,232 @@ UC_Read* g_read) } +void calculate_ul_boundary_cigars(overlap_region* sub_list, const ul_idx_t *uref, Correct_dumy* dumy, +UC_Read* g_read, double max_ov_diff_ec, long long blockLen, void *km) +{ + resize_window_list_alloc(&(sub_list->boundary_cigars), sub_list->w_list_length - 1, km); + int y_id = sub_list->y_id; + int y_strand = sub_list->y_pos_strand; + long long y_readLen = uref->ug->u.a[y_id].len; + long long i, y_distance; + int f_err, b_err, m_error; + int scanLen = 10; + long long boundaryLen = WINDOW_UL_BOUND_RATE*blockLen; + boundaryLen >>= 2; boundaryLen <<= 2; + if(boundaryLen < WINDOW_UL_BOUND) boundaryLen = WINDOW_UL_BOUND; + long long single_sideLen = boundaryLen/2; + long long force_useless_side = single_sideLen/2; + long long L_useless_side, R_useless_side; + int alpha = 1; + long long y_start, x_start; + long long x_end; + long long yLen, xLen; + long long leftLen, rightLen; + long long threshold; + int extra_begin, extra_end; + long long o_len; + char* x_string; + char* y_string; + int end_site; + unsigned int error; + int real_y_start; + sub_list->boundary_cigars.length = sub_list->w_list_length - 1; + ///the (i)-th boundary between the (i)-th window and the (i+1)-th window + ///that means it includes (the tail of (i)-th window) and (the header of (i+1)-th window) + ///note the (i)-th boundary is calculated at the (i)-th window + for (i = 0; i + 1 < (long long)sub_list->w_list_length; i++) + { + ///if both of the two windows are not aligned + ///it is not necessary to calculate the boundary + ///if(sub_list->w_list[i].y_end == -1 && sub_list->w_list[i+1].y_end == -1) + if(sub_list->w_list[i].y_end == -1 || sub_list->w_list[i+1].y_end == -1) + { + sub_list->boundary_cigars.buffer[i].error = -1; + sub_list->boundary_cigars.buffer[i].y_end = -1; + continue; + } + + ///note if w_list[i+1].y_start or sub_list->w_list[i].y_end is -1 + ///y_distance might have some problems at the last of this function + ///we need to deal with it carefully + y_distance = sub_list->w_list[i+1].y_start - sub_list->w_list[i].y_end - 1; + + ///if two windows are aligned + if(sub_list->w_list[i].y_end != -1 && sub_list->w_list[i+1].y_end != -1 && y_distance == 0) + { + ///scan backward + scan_cigar(&(sub_list->w_list[i].cigar), &b_err, scanLen, 1); + ///scan forward + scan_cigar(&(sub_list->w_list[i+1].cigar), &f_err, scanLen, 0); + if(b_err == 0 && f_err == 0) + { + sub_list->boundary_cigars.buffer[i].error = -2; + sub_list->boundary_cigars.buffer[i].y_end = -1; + continue; + } + } + + + ///y_distance can be less than 0, or larger than 0 + if(sub_list->w_list[i].y_end != -1) + { + y_start = sub_list->w_list[i].y_end; + x_start = sub_list->w_list[i].x_end; + }///if the (i)-th window is not matched, have a look at the (i+1)-th window + else if(sub_list->w_list[i+1].y_end != -1) + { + y_start = sub_list->w_list[i+1].y_start; + x_start = sub_list->w_list[i+1].x_start; + }///if both of these two windows are not matched, directly skip + else + { + sub_list->boundary_cigars.buffer[i].error = -1; + sub_list->boundary_cigars.buffer[i].y_end = -1; + continue; + } + + ///it seems we don't need to record x_start and y_start + sub_list->boundary_cigars.buffer[i].extra_begin = x_start; + sub_list->boundary_cigars.buffer[i].extra_end = y_start; + + ///leftLen and rightLen are used for x + ///x should be at [sub_list->w_list[i].x_start, sub_list->w_list[i+1].x_end] + ///y shouldn't have limitation + ///note that the x_start and x_end should not be -1 in any case + ///up to now, x_start and y_start are not -1 + ///leftLen does not include x_start itself, rightLen does + ///gnerally speaking, rightLen should be always larger than leftLen + leftLen = MIN(MIN((x_start - (long long)sub_list->w_list[i].x_start), y_start), single_sideLen); + rightLen = MIN(MIN(((long long)sub_list->w_list[i+1].x_end + 1 - x_start), y_readLen - y_start), + single_sideLen); + + ///xLen should be the sum length of two windows + xLen = leftLen + rightLen; + x_start = x_start - leftLen; + x_end = x_start + xLen - 1; + y_start = y_start - leftLen; + + ///if we don't have enough leftLen and rightLen + // if(leftLen <= useless_side || rightLen <= useless_side) + // { + // sub_list->boundary_cigars.buffer[i].error = -1; + // sub_list->boundary_cigars.buffer[i].y_end = -1; + // continue; + // } + + + threshold = xLen * max_ov_diff_ec; + threshold = Adjust_Threshold(threshold, xLen); + threshold = double_ul_error_threshold(threshold, xLen); + + yLen = xLen + (threshold << 1); + if(!determine_overlap_region(threshold, y_start, y_id, yLen, uref->ug->u.a[y_id].len, + &extra_begin, &extra_end, &y_start, &o_len)) + { + sub_list->boundary_cigars.buffer[i].error = -1; + sub_list->boundary_cigars.buffer[i].y_end = -1; + continue; + } + + if(o_len < xLen) + { + sub_list->boundary_cigars.buffer[i].error = -1; + sub_list->boundary_cigars.buffer[i].y_end = -1; + continue; + } + + fill_subregion_ul(dumy->overlap_region, y_start, o_len, y_strand, + uref, y_id, extra_begin, extra_end); + + x_string = g_read->seq + x_start; + y_string = dumy->overlap_region; + + end_site = Reserve_Banded_BPM_PATH(y_string, yLen, x_string, xLen, threshold, &error, + &real_y_start, &(dumy->path_length), dumy->matrix_bit, dumy->path, -1, -1); + + ///means this window is matched + if (error!=(unsigned int)-1) + { + sub_list->boundary_cigars.buffer[i].x_start = x_start; + sub_list->boundary_cigars.buffer[i].x_end = x_end; + + generate_cigar(dumy->path, dumy->path_length, &(sub_list->boundary_cigars.buffer[i]), + &real_y_start, &end_site, &error, x_string, xLen, y_string); + ///should not adjust cigar here, adjust cigar may cause problem + ///that is not what we want + + ///y_distance can be less than 0, or larger than 0 + ///please if one of the two windows is not matched, + ///y_distance may have potential problems + if(y_distance < 0) y_distance = y_distance * (-1); + ///leftLen, rightLen + // if(leftLen <= useless_side || rightLen <= useless_side) + // { + // sub_list->boundary_cigars.buffer[i].error = -1; + // sub_list->boundary_cigars.buffer[i].y_end = -1; + // continue; + // } + L_useless_side = R_useless_side = force_useless_side; + + ///first window + if((i == 0) && (x_start == (long long)sub_list->w_list[0].x_start)) + { + L_useless_side = 0; + } + ///last window + if((i == (long long)(sub_list->w_list_length - 2)) && + (x_end == (long long)(sub_list->w_list[sub_list->w_list_length - 1].x_end))) + { + R_useless_side = 0; + } + + if(leftLen <= L_useless_side || rightLen <= R_useless_side) + { + sub_list->boundary_cigars.buffer[i].error = -1; + sub_list->boundary_cigars.buffer[i].y_end = -1; + continue; + } + + + + ///up to now, if we require (i)-th window and (i+1)-th window are matched + ///boundary_cigars.buffer[i].cigar, w_list[i].cigar and w_list[i+1].cigar are avaiable + ///get the error excluding the first and the last useless_side bases + scan_cigar_interval(&(sub_list->boundary_cigars.buffer[i].cigar), &m_error, + L_useless_side, xLen-R_useless_side-1); + scan_cigar(&(sub_list->w_list[i].cigar), &b_err, leftLen-L_useless_side, 1); + scan_cigar(&(sub_list->w_list[i+1].cigar), &f_err, rightLen-R_useless_side, 0); + + if(f_err + b_err + y_distance + alpha < m_error) + { + sub_list->boundary_cigars.buffer[i].error = -1; + sub_list->boundary_cigars.buffer[i].y_end = -1; + continue; + } + + sub_list->boundary_cigars.buffer[i].error = error; + sub_list->boundary_cigars.buffer[i].y_start = y_start + real_y_start - extra_begin; + sub_list->boundary_cigars.buffer[i].y_end = y_start + end_site - extra_begin; + + sub_list->boundary_cigars.buffer[i].x_start = x_start; + sub_list->boundary_cigars.buffer[i].x_end = x_end; + ///sub_list->boundary_cigars.buffer[i].error_threshold = useless_side; + sub_list->boundary_cigars.buffer[i].extra_begin = L_useless_side; + sub_list->boundary_cigars.buffer[i].extra_end = R_useless_side; + } + else + { + sub_list->boundary_cigars.buffer[i].error = -1; + sub_list->boundary_cigars.buffer[i].y_end = -1; + continue; + } + + } + +} + + + void debug_window_cigar(overlap_region_alloc* overlap_list, UC_Read* g_read, Correct_dumy* dumy, All_reads* R_INF, int test_window, int test_boundary) { @@ -2933,7 +3607,7 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea y_start = total_y_start; Window_Len = x_len + (threshold << 1); - if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF, + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, Get_READ_LENGTH((*R_INF), y_id), &extra_begin, &extra_end, &y_start, &o_len)) { break; @@ -3087,7 +3761,7 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea ///y_start might be less than 0 y_start = total_y_end - x_len + 1; - if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF, + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, Get_READ_LENGTH((*R_INF), y_id), &extra_begin, &extra_end, &y_start, &o_len)) { break; @@ -3261,6 +3935,507 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea } + +void set_herror_win(overlap_region_alloc* ovlp, Correct_dumy* dumy, double max_ov_diff_ec, long long rLen, long long blockLen) +{ + Window_Pool w_inf; + init_Window_Pool(&w_inf, rLen, blockLen, (int)(1.0/max_ov_diff_ec)); + int flag = 0; long long window_start, window_end, cID, w_list_id, i, k, mLen, fc, fw; + while(get_Window(&w_inf, &window_start, &window_end) && flag != -2) { + dumy->length = 0; + dumy->lengthNT = 0; + flag = get_interval(window_start, window_end, ovlp, dumy, w_inf.window_length); + switch (flag) { + case 1: ///no match here + break; + case 0: ///no match here + break; + case -2: ///if flag == -2, loop would be terminated + break; + } + + + for (i = 0; i < (long long)dumy->length; i++) { + cID = dumy->overlapID[i]; + if(ovlp->list[cID].is_match!=3 && ovlp->list[cID].is_match!=4) continue; + w_list_id = (window_start - (ovlp->list[cID].x_pos_s / w_inf.window_length)* + w_inf.window_length)/w_inf.window_length; + if (ovlp->list[cID].w_list[w_list_id].y_end != -1) break; + } + if(i < (long long)dumy->length) continue; + for (i = 0; i < (long long)dumy->length; i++) { + cID = dumy->overlapID[i]; + if(ovlp->list[cID].is_match!=3 && ovlp->list[cID].is_match!=4) continue; + w_list_id = (window_start - (ovlp->list[cID].x_pos_s / w_inf.window_length)* + w_inf.window_length)/w_inf.window_length; + if (ovlp->list[cID].w_list[w_list_id].y_end == -1) { + ovlp->list[cID].w_list[w_list_id].y_end = -2; + ovlp->list[cID].w_list[w_list_id].error = THRESHOLD_MAX_SIZE; + ovlp->list[cID].align_length += ovlp->list[cID].w_list[w_list_id].x_end + 1 + - ovlp->list[cID].w_list[w_list_id].x_start; + ovlp->list[cID].is_match = 4; + } + } + + + ///shorter than blockLen + for (i = dumy->size - dumy->lengthNT, mLen = dumy->size - dumy->lengthNT; i < (long long)dumy->size; i++) { + cID = dumy->overlapID[i]; + if(ovlp->list[cID].is_match!=3 && ovlp->list[cID].is_match!=4) continue; + w_list_id = (window_start - (ovlp->list[cID].x_pos_s / w_inf.window_length)* + w_inf.window_length)/w_inf.window_length; + if (ovlp->list[cID].w_list[w_list_id].y_end != -1) { + if(mLen != i) { + k = dumy->overlapID[i]; + dumy->overlapID[i] = dumy->overlapID[mLen]; + dumy->overlapID[mLen] = k; + } + mLen++; + } + } + if(mLen == (long long)dumy->size) continue; + for (i = mLen; i < (long long)dumy->size; i++){ + cID = dumy->overlapID[i]; + if(ovlp->list[cID].is_match!=3 && ovlp->list[cID].is_match!=4) continue; + w_list_id = (window_start - (ovlp->list[cID].x_pos_s / w_inf.window_length)* + w_inf.window_length)/w_inf.window_length; + if (ovlp->list[cID].w_list[w_list_id].y_end != -1) fprintf(stderr, "ERROR\n"); + for (k = dumy->size - dumy->lengthNT; k < mLen; k++) { + fc = dumy->overlapID[k]; + fw = (window_start - (ovlp->list[fc].x_pos_s / w_inf.window_length)* + w_inf.window_length)/w_inf.window_length; + if (ovlp->list[fc].w_list[fw].y_end == -1 || (ovlp->list[fc].is_match!=3 && ovlp->list[fc].is_match!=4)) fprintf(stderr, "ERROR\n"); + if(ovlp->list[fc].w_list[fw].x_start<=ovlp->list[cID].w_list[w_list_id].x_start && + ovlp->list[fc].w_list[fw].x_end>=ovlp->list[cID].w_list[w_list_id].x_end) { + break; + } + } + + if(k >= mLen) { + ovlp->list[cID].w_list[w_list_id].y_end = -2; + ovlp->list[cID].w_list[w_list_id].error = THRESHOLD_MAX_SIZE; + ovlp->list[cID].align_length += ovlp->list[cID].w_list[w_list_id].x_end + 1 + - ovlp->list[cID].w_list[w_list_id].x_start; + ovlp->list[cID].is_match = 4; + } + } + } +} + +inline void recalcate_window_ul_advance(overlap_region_alloc* overlap_list, const ul_idx_t *uref, + UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, double max_ov_diff_ec, long long blockLen, void *km) +{ + long long j, k, i; + int threshold; + long long y_id; + int y_strand; + long long y_readLen; + long long x_start; + long long x_end; + long long x_len; + long long total_y_start; + long long total_y_end; + long long y_start; + long long Window_Len; + char* x_string; + char* y_string; + int end_site; + unsigned int error; + int real_y_start; + long long overlap_length; + int extra_begin, extra_end; + long long o_len; + + overlap_list->mapped_overlaps_length = 0; + for (j = 0; j < (long long)overlap_list->length; j++) + { + + y_id = overlap_list->list[j].y_id; + y_strand = overlap_list->list[j].y_pos_strand; + y_readLen = uref->ug->u.a[y_id].len; + + //i corresponding to each window of a overlap + //utilize the the end pos of pre-window in backwards + for (i = overlap_list->list[j].w_list_length - 1; i >= 0; i--) + { + ///the first matched window + if(overlap_list->list[j].w_list[i].y_end != -1) + { + ///note!!! need notification + ///this is the actual end postion in ystring + total_y_start = overlap_list->list[j].w_list[i].y_end + - overlap_list->list[j].w_list[i].extra_begin + 1; + + ///k corresponding to all unmatched windows at the right side of overlap_list->list[j].w_list[i] + ///so k starts from i + 1, and end to the first matched window + for (k = i + 1; k < (long long)overlap_list->list[j].w_list_length && overlap_list->list[j].w_list[k].y_end == -1; k++) + { + extra_begin = extra_end = 0; + + ///if y_start > y_readLen, direct terminate + if (total_y_start >= y_readLen) + { + break; + } + + ///there is no problem for x + x_start = overlap_list->list[j].w_list[k].x_start; + x_end = overlap_list->list[j].w_list[k].x_end; + x_len = x_end - x_start + 1; + ///there are two potiential reasons for unmatched window: + ///1. this window has a large number of differences + ///2. DP does not start from the right offset + threshold = double_ul_error_threshold(overlap_list->list[j].w_list[k].error_threshold, x_len); + + y_start = total_y_start; + Window_Len = x_len + (threshold << 1); + + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, uref->ug->u.a[y_id].len, + &extra_begin, &extra_end, &y_start, &o_len)) + { + break; + } + + if(o_len + threshold < x_len) + { + break; + } + + fill_subregion_ul(dumy->overlap_region, y_start, o_len, y_strand, uref, y_id, extra_begin, extra_end); + + x_string = g_read->seq + x_start; + y_string = dumy->overlap_region; + + ///note!!! need notification + end_site = Reserve_Banded_BPM(y_string, Window_Len, x_string, x_len, threshold, &error); + + ///if error==-1, unmatched + if (error!=(unsigned int)-1) + { + overlap_list->list[j].w_list[k].cigar.length = -1; + overlap_list->list[j].w_list[k].y_start = y_start; + overlap_list->list[j].w_list[k].y_end = y_start + end_site; + overlap_list->list[j].w_list[k].error = (int)error; + ///note!!! need notification + overlap_list->list[j].w_list[k].extra_begin = extra_begin; + overlap_list->list[j].w_list[k].extra_end = extra_end; + overlap_list->list[j].w_list[k].error_threshold = threshold; + + overlap_list->list[j].align_length += x_len; + } + else + { + break; + } + + ///note!!! need notification + total_y_start = y_start + end_site - extra_begin + 1; + } + + } + + } + + + + //i corresponding to each window of a overlap + //utilize the the start pos of next window in forward + for (i = 0; i < (long long)overlap_list->list[j].w_list_length; i++) + { + ///find the first matched window, which should not be the first window + ///the pre-window of this matched window must be unmatched + if(overlap_list->list[j].w_list[i].y_end != -1 && i != 0 && overlap_list->list[j].w_list[i - 1].y_end == -1) + { + ///check if the start pos of this matched window has been calculated + if(overlap_list->list[j].w_list[i].cigar.length == -1) + { + ///there is no problem for x + x_start = overlap_list->list[j].w_list[i].x_start; + x_end = overlap_list->list[j].w_list[i].x_end; + x_len = x_end - x_start + 1; + /****************************may have bugs********************************/ + threshold = overlap_list->list[j].w_list[i].error_threshold; + /****************************may have bugs********************************/ + /****************************may have bugs********************************/ + ///should not adjust threshold, since this window can be matched by the old threshold + ///threshold = Adjust_Threshold(threshold, x_len); + /****************************may have bugs********************************/ + Window_Len = x_len + (threshold << 1); + + + ///y_start is the real y_start + y_start = overlap_list->list[j].w_list[i].y_start; + extra_begin = overlap_list->list[j].w_list[i].extra_begin; + extra_end = overlap_list->list[j].w_list[i].extra_end; + o_len = Window_Len - extra_end - extra_begin; + fill_subregion_ul(dumy->overlap_region, y_start, o_len, y_strand, uref, y_id, extra_begin, extra_end); + x_string = g_read->seq + x_start; + y_string = dumy->overlap_region; + + ///note!!! need notification + end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start, + &(dumy->path_length), dumy->matrix_bit, dumy->path, + overlap_list->list[j].w_list[i].error, overlap_list->list[j].w_list[i].y_end - y_start); + + + ///y_start has already been calculated + if (error != (unsigned int)-1) + { + ///this condition is always wrong + ///in best case, real_y_start = threshold, end_site = Window_Len - threshold - 1 + if (end_site == Window_Len - 1 || real_y_start == 0) + { + if(fix_ul_boundary(x_string, x_len, threshold, y_start, real_y_start, + end_site, extra_begin, extra_end, y_id, Window_Len, uref, dumy, + y_strand, error, &y_start, &real_y_start, &end_site, &extra_begin, + &extra_end, &error)) + { + overlap_list->list[j].w_list[i].error = error; + overlap_list->list[j].w_list[i].extra_begin = extra_begin; + overlap_list->list[j].w_list[i].extra_end = extra_end; + } + } + + generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[i]), + &real_y_start, &end_site, &error, x_string, x_len, y_string); + + ///note!!! need notification + real_y_start = y_start + real_y_start - extra_begin; + overlap_list->list[j].w_list[i].y_start = real_y_start; + ///I forget why don't reduce the extra_begin for y_end + ///it seems extra_begin will be reduced at the end of this function + overlap_list->list[j].w_list[i].y_end = y_start + end_site; + overlap_list->list[j].w_list[i].error = error; + } + else + { + fprintf(stderr, "error\n"); + } + } + else + { + real_y_start = overlap_list->list[j].w_list[i].y_start; + } + + + ///the end pos for pre window is real_y_start - 1 + total_y_end = real_y_start - 1; + ///find the unmatched window on the left of current matched window + ///k starts from i - 1 + for (k = i - 1; k >= 0 && overlap_list->list[j].w_list[k].y_end == -1; k--) + { + ///there is no problem in x + x_start = overlap_list->list[j].w_list[k].x_start; + x_end = overlap_list->list[j].w_list[k].x_end; + x_len = x_end - x_start + 1; + ///there are two potiential reasons for unmatched window: + ///1. this window has a large number of differences + ///2. DP does not start from the right offset + threshold = double_ul_error_threshold(overlap_list->list[j].w_list[k].error_threshold, x_len); + + Window_Len = x_len + (threshold << 1); + + if(total_y_end <= 0) + { + break; + } + + ///y_start might be less than 0 + y_start = total_y_end - x_len + 1; + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, uref->ug->u.a[y_id].len, + &extra_begin, &extra_end, &y_start, &o_len)) + { + break; + } + + if(o_len + threshold < x_len) + { + break; + } + + fill_subregion_ul(dumy->overlap_region, y_start, o_len, y_strand, + uref, y_id, extra_begin, extra_end); + x_string = g_read->seq + x_start; + y_string = dumy->overlap_region; + + ///note!!! need notification + end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start, + &(dumy->path_length), dumy->matrix_bit, dumy->path, -1, -1); + + if (error!=(unsigned int)-1) + { + ///this condition is always wrong + ///in best case, real_y_start = threshold, end_site = Window_Len - threshold - 1 + if (end_site == Window_Len - 1 || real_y_start == 0) + { + fix_ul_boundary(x_string, x_len, threshold, y_start, real_y_start, end_site, + extra_begin, extra_end, y_id, Window_Len, uref, dumy, y_strand, error, + &y_start, &real_y_start, &end_site, + &extra_begin, &extra_end, &error); + } + + generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[k]), + &real_y_start, &end_site, &error, x_string, x_len, y_string); + + ///y_start has no shift, but y_end has shift + overlap_list->list[j].w_list[k].y_start = y_start + real_y_start - extra_begin; + overlap_list->list[j].w_list[k].y_end = y_start + end_site; + overlap_list->list[j].w_list[k].error = error; + overlap_list->list[j].align_length += x_len; + overlap_list->list[j].w_list[k].extra_begin = extra_begin; + overlap_list->list[j].w_list[k].extra_end = extra_end; + overlap_list->list[j].w_list[k].error_threshold = threshold; + } + else + { + break; + } + + total_y_end = y_start + real_y_start - 1 - extra_begin; + } + } + } + + overlap_list->list[j].is_match = 0; + if((((overlap_list->list[j].x_pos_e + 1 - overlap_list->list[j].x_pos_s)*MIN_UL_ALIN_RATE) + <= overlap_list->list[j].align_length) && (overlap_list->list[j].align_length >= MIN_UL_ALIN_LEN)){ + overlap_list->list[j].is_match = 3; + overlap_list->mapped_overlaps_length += overlap_list->list[j].align_length; + } + } + + if(overlap_list->mapped_overlaps_length > 0) set_herror_win(overlap_list, dumy, max_ov_diff_ec, g_read->length, blockLen); + + overlap_list->mapped_overlaps_length = 0; + + double error_rate; + for (j = 0; j < (long long)overlap_list->length; j++) + { + y_id = overlap_list->list[j].y_id; + y_strand = overlap_list->list[j].y_pos_strand; + y_readLen = uref->ug->u.a[y_id].len; + overlap_length = overlap_list->list[j].x_pos_e - overlap_list->list[j].x_pos_s + 1; + // overlap_list->list[j].is_match = 0; + + ///debug_scan_cigar(&(overlap_list->list[j])); + + ///only calculate cigar for high quality overlaps + if (overlap_length * (1-max_ov_diff_ec) <= overlap_list->list[j].align_length) + { + + for (i = 0; i < (long long)overlap_list->list[j].w_list_length; i++) + { + ///first we need to check if this window is matched + if(overlap_list->list[j].w_list[i].y_end != -1 && overlap_list->list[j].w_list[i].y_end != -2) + { + ///second check if the cigar of this window has been got + if(overlap_list->list[j].w_list[i].cigar.length == -1) + { + ///there is no problem for x + x_start = overlap_list->list[j].w_list[i].x_start; + x_end = overlap_list->list[j].w_list[i].x_end; + x_len = x_end - x_start + 1; + /****************************may have bugs********************************/ + ///threshold = x_len * asm_opt.max_ov_diff_ec; + threshold = overlap_list->list[j].w_list[i].error_threshold; + /****************************may have bugs********************************/ + /****************************may have bugs********************************/ + ///should not adjust threshold, since this window can be matched by the old threshold + ///threshold = Adjust_Threshold(threshold, x_len); + /****************************may have bugs********************************/ + Window_Len = x_len + (threshold << 1); + + + ///y_start is the real y_start + ///for the window with cigar, y_start has already reduced extra_begin + y_start = overlap_list->list[j].w_list[i].y_start; + extra_begin = overlap_list->list[j].w_list[i].extra_begin; + extra_end = overlap_list->list[j].w_list[i].extra_end; + o_len = Window_Len - extra_end - extra_begin; + fill_subregion_ul(dumy->overlap_region, y_start, o_len, y_strand, + uref, y_id, extra_begin, extra_end); + x_string = g_read->seq + x_start; + y_string = dumy->overlap_region; + + + ///note!!! need notification + end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start, + &(dumy->path_length), dumy->matrix_bit, dumy->path, + overlap_list->list[j].w_list[i].error, overlap_list->list[j].w_list[i].y_end - y_start); + + if (error != (unsigned int)-1) + { + if (end_site == Window_Len - 1 || real_y_start == 0) + { + + if(fix_ul_boundary(x_string, x_len, threshold, y_start, real_y_start, end_site, + extra_begin, extra_end, y_id, Window_Len, uref, dumy, y_strand, error, + &y_start, &real_y_start, &end_site, + &extra_begin, &extra_end, &error)) + { + overlap_list->list[j].w_list[i].error = error; + overlap_list->list[j].w_list[i].extra_begin = extra_begin; + overlap_list->list[j].w_list[i].extra_end = extra_end; + } + + } + + generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[i]), + &real_y_start, &end_site, &error, x_string, x_len, y_string); + + ///note!!! need notification + real_y_start = y_start + real_y_start - extra_begin; + overlap_list->list[j].w_list[i].y_start = real_y_start; + overlap_list->list[j].w_list[i].y_end = y_start + end_site - extra_begin; + overlap_list->list[j].w_list[i].error = error; + } + else + { + fprintf(stderr, "error\n"); + } + } + else + { + overlap_list->list[j].w_list[i].y_end -= overlap_list->list[j].w_list[i].extra_begin; + } + + + } + } + + error_rate = non_trim_ul_error_rate(overlap_list, j, uref, dumy, g_read); + + if(overlap_list->list[j].is_match == 4) { + for (i = 0; i < (long long)overlap_list->list[j].w_list_length; i++) { + if(overlap_list->list[j].w_list[i].y_end == -2) overlap_list->list[j].w_list[i].y_end = -1; + } + } + overlap_list->list[j].is_match = 0; + + if (error_rate <= max_ov_diff_ec) + { + overlap_list->mapped_overlaps_length += overlap_length; + overlap_list->list[j].is_match = 1; + calculate_ul_boundary_cigars(&(overlap_list->list[j]), uref, dumy, g_read, max_ov_diff_ec, blockLen, km); + } + else if (error_rate <= max_ov_diff_ec * 1.5) + { + overlap_list->list[j].is_match = 3; + } + } else { + if(overlap_list->list[j].is_match == 4) { + for (i = 0; i < (long long)overlap_list->list[j].w_list_length; i++) { + if(overlap_list->list[j].w_list[i].y_end == -2) overlap_list->list[j].w_list[i].y_end = -1; + } + } + overlap_list->list[j].is_match = 0; + } + } + + ///debug_window_cigar(overlap_list, g_read, dumy, R_INF, 1, 1); +} + + inline void add_base_to_correct_read_directly(Correct_dumy* dumy, char base) { @@ -4808,7 +5983,7 @@ Cigar_record* current_cigar, long long uncorrected_window_start, Round2_alignmen Window_Len = x_len + (threshold << 1); error =(unsigned int)-1; - if(determine_overlap_region(threshold, y_start, overlap_list->list[overlapID].y_id, Window_Len, R_INF, + if(determine_overlap_region(threshold, y_start, overlap_list->list[overlapID].y_id, Window_Len, Get_READ_LENGTH((*R_INF), overlap_list->list[overlapID].y_id), &extra_begin, &extra_end, &y_start, &o_len)) { fill_subregion(dumy->overlap_region, y_start, o_len, overlap_list->list[overlapID].y_pos_strand, @@ -4855,7 +6030,7 @@ Cigar_record* current_cigar, long long uncorrected_window_start, Round2_alignmen } error =(unsigned int)-1; - if(determine_overlap_region(threshold, y_start, overlap_list->list[overlapID].y_id, Window_Len, R_INF, + if(determine_overlap_region(threshold, y_start, overlap_list->list[overlapID].y_id, Window_Len, Get_READ_LENGTH((*R_INF), overlap_list->list[overlapID].y_id), &extra_begin, &extra_end, &y_start, &o_len)) { fill_subregion(dumy->overlap_region, y_start, o_len, overlap_list->list[overlapID].y_pos_strand, @@ -5050,11 +6225,11 @@ void generate_consensus(overlap_region_alloc* overlap_list, All_reads* R_INF, inline int get_available_fully_covered_interval(long long window_start, long long window_end, overlap_region_alloc* overlap_list, Correct_dumy* dumy, long long* real_length, long long* real_length_100) { - long long i; + long long i, fud = 0; long long Len; long long overlap_length; - + if(window_start == 0) dumy->start_i = 0; for (i = dumy->start_i; i < (long long)overlap_list->length; i++) { if (window_end < (long long)overlap_list->list[i].x_pos_s) @@ -5082,7 +6257,7 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, long long* real_length, long long fake_length = 0; overlap_length = window_end - window_start + 1; - (*real_length) = 0; + (*real_length) = 0; fud = 0; for (; i < (long long)overlap_list->length; i++) { @@ -5099,6 +6274,7 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, long long* real_length, { (*real_length_100)++; } + if(fud == 0) fud = 1, dumy->start_i = i; } if((long long)overlap_list->list[i].x_pos_s > window_end) @@ -5235,13 +6411,14 @@ CIGAR* cigar, haplotype_evdience_alloc* hap) ///mark SNPs at [xBeg, xEnd], note we need to deal with flag_offset carefully void markSNP_detail(CIGAR* cigar_record, uint8_t* flag, -long long xBeg, long long xEnd, long long flag_offset) +long long xBeg, long long xEnd, long long flag_offset, const ul_idx_t *uref, long long y_total_start, int y_strand, int yid) { if(xBeg > xEnd) return; int operation; int operationLen; long long x_i, y_i, cigar_i, i; + int64_t pi = 0, cc = 0; i = cigar_i = x_i = y_i = 0; while (cigar_i < cigar_record->length) @@ -5267,7 +6444,14 @@ long long xBeg, long long xEnd, long long flag_offset) ///if(flag[x_i - flag_offset] < 127 && x_i >= xBeg && x_i <= xEnd) if(x_i >= xBeg && x_i <= xEnd && flag[x_i - flag_offset] < 127)///Fix-attention { - flag[x_i - flag_offset]++; + if(uref) { + cc = retrieve_u_cov(uref, yid, y_strand, y_total_start + y_i, y_strand, &pi); + cc += flag[x_i - flag_offset]; + flag[x_i - flag_offset] = (cc <= 127?cc:127); + // if(cc <= 127) flag[x_i - flag_offset] = cc; + } else { + flag[x_i - flag_offset]++; + } } x_i++; @@ -5295,7 +6479,7 @@ long long window_offset, long long x_total_start, long long x_length, long long y_total_start, long long y_length, window_list* current_cigar, window_list* beg_cigar, window_list* end_cigar, -haplotype_evdience_alloc* hap) +haplotype_evdience_alloc* hap, const ul_idx_t *uref, int strand, int yid) { long long x_total_end = x_total_start + x_length - 1; ///mismatches based on the offset of x @@ -5335,8 +6519,8 @@ haplotype_evdience_alloc* hap) ///current_cigar_beg = xrightLen - useless_side; current_cigar_beg = xrightLen - R_useless_side; - markSNP_detail(cigar_record, hap->flag + inner_offset, x_interval_beg, - x_interval_end, x_interval_beg); + markSNP_detail(cigar_record, hap->flag + inner_offset, x_interval_beg, x_interval_end, + x_interval_beg, uref, beg_cigar->y_start, strand, yid); } } @@ -5371,13 +6555,13 @@ haplotype_evdience_alloc* hap) ///current_cigar_end = (x_length - 1) - ((xleftLen + 1) - useless_side); current_cigar_end = (x_length - 1) - ((xleftLen + 1) - L_useless_side); - markSNP_detail(cigar_record, hap->flag + end_cigar->x_start - window_offset, - x_interval_beg, x_interval_end, 0); + markSNP_detail(cigar_record, hap->flag + end_cigar->x_start - window_offset, x_interval_beg, + x_interval_end, 0, uref, end_cigar->y_start, strand, yid); } } markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, current_cigar_beg, - current_cigar_end, 0); + current_cigar_end, 0, uref, current_cigar->y_start, strand, yid); } @@ -5422,7 +6606,7 @@ CIGAR* cigar, haplotype_evdience_alloc* hap, int snp_threshold) ev.site = x_total_start + x_i; ev.overlapSite = y_total_start + y_i; ev.type = 0; - addHaplotypeEvdience(hap, &ev); + addHaplotypeEvdience(hap, &ev, NULL); } @@ -5444,7 +6628,7 @@ CIGAR* cigar, haplotype_evdience_alloc* hap, int snp_threshold) ev.site = x_total_start + x_i; ev.overlapSite = y_total_start + y_i; ev.type = 1; - addHaplotypeEvdience(hap, &ev); + addHaplotypeEvdience(hap, &ev, NULL); } inner_offset++; @@ -5468,7 +6652,7 @@ CIGAR* cigar, haplotype_evdience_alloc* hap, int snp_threshold) ev.site = x_total_start + x_i; ev.overlapSite = y_total_start + y_i; ev.type = 2; - addHaplotypeEvdience(hap, &ev); + addHaplotypeEvdience(hap, &ev, NULL); } inner_offset++; @@ -5486,13 +6670,14 @@ CIGAR* cigar, haplotype_evdience_alloc* hap, int snp_threshold) void addSNPtohaplotype_details(CIGAR* cigar_record, uint8_t* flag, char* x_string, char* y_string, long long x_total_start, long long y_total_start, long long xBeg, long long xEnd, int overlapID, long long flag_offset, -haplotype_evdience_alloc* hap, long long snp_threshold) +haplotype_evdience_alloc* hap, long long snp_threshold, const ul_idx_t *uref, int y_strand, int yid, void *km) { if(xBeg > xEnd) return; int operation; int operationLen; long long x_i, y_i, cigar_i, i; i = cigar_i = x_i = y_i = 0; + int64_t pi = 0; haplotype_evdience ev; @@ -5524,7 +6709,8 @@ haplotype_evdience_alloc* hap, long long snp_threshold) ev.site = x_total_start + x_i; ev.overlapSite = y_total_start + y_i; ev.type = 0; - addHaplotypeEvdience(hap, &ev); + ev.cov = uref?retrieve_u_cov(uref, yid, y_strand, y_total_start + y_i, y_strand, &pi):1; + addHaplotypeEvdience(hap, &ev, km); } ///inner_offset++; x_i++; @@ -5546,7 +6732,8 @@ haplotype_evdience_alloc* hap, long long snp_threshold) ev.site = x_total_start + x_i; ev.overlapSite = y_total_start + y_i; ev.type = 1; - addHaplotypeEvdience(hap, &ev); + ev.cov = uref?retrieve_u_cov(uref, yid, y_strand, y_total_start + y_i, y_strand, &pi):1; + addHaplotypeEvdience(hap, &ev, km); } ///inner_offset++; @@ -5574,7 +6761,8 @@ haplotype_evdience_alloc* hap, long long snp_threshold) ev.site = x_total_start + x_i; ev.overlapSite = y_total_start + y_i; ev.type = 2; - addHaplotypeEvdience(hap, &ev); + ev.cov = uref?retrieve_u_cov(uref, yid, y_strand, y_total_start + y_i, y_strand, &pi):1; + addHaplotypeEvdience(hap, &ev, km); } ///inner_offset++; @@ -5592,7 +6780,7 @@ long long window_offset, int overlapID, long long x_total_start, long long x_length, long long y_total_start, long long y_length, window_list* current_cigar, window_list* beg_cigar, window_list* end_cigar, -haplotype_evdience_alloc* hap, int snp_threshold, char* x_T_string, char* y_T_string) +haplotype_evdience_alloc* hap, int snp_threshold, char* x_T_string, char* y_T_string, const ul_idx_t *uref, int strand, int yid, void *km) { long long x_total_end = x_total_start + x_length - 1; long long inner_offset = x_total_start - window_offset; @@ -5636,7 +6824,7 @@ haplotype_evdience_alloc* hap, int snp_threshold, char* x_T_string, char* y_T_st addSNPtohaplotype_details(cigar_record, hap->flag + inner_offset, x_T_string + beg_cigar->x_start, y_T_string + beg_cigar->y_start, beg_cigar->x_start, beg_cigar->y_start, x_interval_beg, x_interval_end, - overlapID, x_interval_beg, hap, snp_threshold); + overlapID, x_interval_beg, hap, snp_threshold, uref, strand, yid, km); } } @@ -5677,7 +6865,7 @@ haplotype_evdience_alloc* hap, int snp_threshold, char* x_T_string, char* y_T_st addSNPtohaplotype_details(cigar_record, hap->flag + end_cigar->x_start - window_offset, x_T_string + end_cigar->x_start, y_T_string + end_cigar->y_start, end_cigar->x_start, end_cigar->y_start, x_interval_beg, x_interval_end, - overlapID, 0, hap, snp_threshold); + overlapID, 0, hap, snp_threshold, uref, strand, yid, km); } } @@ -5686,7 +6874,7 @@ haplotype_evdience_alloc* hap, int snp_threshold, char* x_T_string, char* y_T_st addSNPtohaplotype_details(&(current_cigar->cigar), hap->flag + inner_offset, x_T_string + current_cigar->x_start, y_T_string + current_cigar->y_start, current_cigar->x_start, current_cigar->y_start, current_cigar_beg, - current_cigar_end, overlapID, 0, hap, snp_threshold); + current_cigar_end, overlapID, 0, hap, snp_threshold, uref, strand, yid, km); } @@ -5912,7 +7100,8 @@ haplotype_evdience_alloc* hap, UC_Read* overlap_read, int snp_threshold) markSNP_advance(window_start, x_start, x_length, y_start, y_length, - &(overlap_list->list[overlapID].w_list[windowID]), beg_cigar, end_cigar, hap); + &(overlap_list->list[overlapID].w_list[windowID]), beg_cigar, end_cigar, hap, NULL, + overlap_list->list[overlapID].y_pos_strand, overlap_list->list[overlapID].y_id); } @@ -6001,7 +7190,8 @@ haplotype_evdience_alloc* hap, UC_Read* overlap_read, int snp_threshold) addSNPtohaplotype_advance(window_start, overlapID, x_start, x_length, y_start, y_length, &(overlap_list->list[overlapID].w_list[windowID]), - beg_cigar, end_cigar, hap, snp_threshold, x_string, y_string); + beg_cigar, end_cigar, hap, snp_threshold, x_string, y_string, NULL, + overlap_list->list[overlapID].y_pos_strand, overlap_list->list[overlapID].y_id, NULL); } RsetInitHaplotypeEvdienceFlag(hap, first_snp, last_snp + 1 - first_snp); @@ -6009,6 +7199,160 @@ haplotype_evdience_alloc* hap, UC_Read* overlap_read, int snp_threshold) if(hap->length - ll > 1 && lr > 1) radix_sort_haplotype_evdience_srt(hap->list+ll, hap->list + hap->length); } +void cluster_ul_advance(char* r_string, long long window_start, long long window_end, +overlap_region_alloc* overlap_list, Correct_dumy* dumy, const ul_idx_t *uref, +haplotype_evdience_alloc* hap, UC_Read* overlap_read, int snp_threshold, long long blockLen, void *km) +{ + window_list* beg_cigar; + window_list* end_cigar; + ///window_start, window_end, and useful_length correspond to x, instead of y + long long useful_length = window_end - window_start + 1; + long long x_start, x_length, ll = hap->length, lr; + char* x_string; + char* y_string; + long long i; + long long y_start, y_length; + long long overlapID, windowID; + + long long correct_x_pos_s; + + + ///all overlaps related to the current window [window_start, window_end] + ///first mark all snp pos + for (i = 0; i < (long long)dumy->length; i++) + { + ///overlap id, instead of the window id or the y id + overlapID = dumy->overlapID[i]; + + ///overlap_list->list[overlapID].x_pos_s is the begining of the whole overlap + correct_x_pos_s = (overlap_list->list[overlapID].x_pos_s / blockLen) * blockLen; + ///window_start is the begining of this window in the whole x_read + windowID = (window_start - correct_x_pos_s) / blockLen; + + ///skip if this window is not matched + if (overlap_list->list[overlapID].w_list[windowID].y_end == -1) + { + continue; + } + + ///both x_start and y_start are the offsets of the whole x_read and y_read + ///instead of the offsets of window + x_start = overlap_list->list[overlapID].w_list[windowID].x_start; + x_length = overlap_list->list[overlapID].w_list[windowID].x_end + - overlap_list->list[overlapID].w_list[windowID].x_start + 1; + + y_start = overlap_list->list[overlapID].w_list[windowID].y_start; + y_length = overlap_list->list[overlapID].w_list[windowID].y_end + - overlap_list->list[overlapID].w_list[windowID].y_start + 1; + + beg_cigar = end_cigar = NULL; + + if(windowID >= 1) + { + beg_cigar = &(overlap_list->list[overlapID].boundary_cigars.buffer[windowID-1]); + } + + if(windowID < (long long)(overlap_list->list[overlapID].w_list_length - 1)) + { + end_cigar = &(overlap_list->list[overlapID].boundary_cigars.buffer[windowID]); + } + + + markSNP_advance(window_start, x_start, x_length, y_start, y_length, + &(overlap_list->list[overlapID].w_list[windowID]), beg_cigar, end_cigar, hap, uref, + overlap_list->list[overlapID].y_pos_strand, overlap_list->list[overlapID].y_id); + } + + + /****************************may have bugs********************************/ + long long last_snp = -1; + long long first_snp = -1; + for (i = 0, lr = 0; i < useful_length; i++) + { + if(hap->flag[i] != 0) + { + last_snp = i; + if(first_snp == -1) + { + first_snp = i; + } + } + ///for a real snp, the coverage should be at least 2 + if(hap->flag[i] > snp_threshold) + { + // hap->snp++; + hap->nn_snp++; + lr++; + } + } + ///if there are any >0 elements, both first_snp and last_snp should be != -1 + if(first_snp == -1 || last_snp == -1) + { + first_snp = 0; + last_snp = -1; + } + /****************************may have bugs********************************/ + + + + ///add the information related to snp to haplotype_evdience_alloc + for (i = 0; i < (long long)dumy->length; i++) + { + ///overlap ID, instead of the window ID + overlapID = dumy->overlapID[i]; + + ///overlap_list->list[overlapID].x_pos_s is the begining of the whole overlap + correct_x_pos_s = (overlap_list->list[overlapID].x_pos_s / blockLen) * blockLen; + ///window_start is the begining of this window in the whole x_read + windowID = (window_start - correct_x_pos_s) / blockLen; + + ///skip if this window is not matched + if (overlap_list->list[overlapID].w_list[windowID].y_end == -1) + { + continue; + } + + ///both x_start and y_start are the offsets of the whole x_read and y_read + ///instead of the offsets of window + x_start = overlap_list->list[overlapID].w_list[windowID].x_start; + x_length = overlap_list->list[overlapID].w_list[windowID].x_end + - overlap_list->list[overlapID].w_list[windowID].x_start + 1; + + y_start = overlap_list->list[overlapID].w_list[windowID].y_start; + y_length = overlap_list->list[overlapID].w_list[windowID].y_end + - overlap_list->list[overlapID].w_list[windowID].y_start + 1; + + retrieve_u_seq(overlap_read, NULL, &uref->ug->u.a[overlap_list->list[overlapID].y_id], + overlap_list->list[overlapID].y_pos_strand, 0, -1, km); + + + x_string = r_string; + y_string = overlap_read->seq; + + + beg_cigar = end_cigar = NULL; + if(windowID >= 1) + { + beg_cigar = &(overlap_list->list[overlapID].boundary_cigars.buffer[windowID-1]); + } + if(windowID < (long long)(overlap_list->list[overlapID].w_list_length - 1)) + { + end_cigar = &(overlap_list->list[overlapID].boundary_cigars.buffer[windowID]); + } + + + addSNPtohaplotype_advance(window_start, overlapID, x_start, x_length, + y_start, y_length, &(overlap_list->list[overlapID].w_list[windowID]), + beg_cigar, end_cigar, hap, snp_threshold, x_string, y_string, uref, + overlap_list->list[overlapID].y_pos_strand, overlap_list->list[overlapID].y_id, km); + } + + RsetInitHaplotypeEvdienceFlag(hap, first_snp, last_snp + 1 - first_snp); + + if(hap->length - ll > 1 && lr > 1) radix_sort_haplotype_evdience_srt(hap->list+ll, hap->list + hap->length); +} + + int cmp_snp_stats(const void * a, const void * b) { if ((*(SnpStats*)a).score != (*(SnpStats*)b).score) @@ -6778,8 +8122,7 @@ void process_repeat_snps(haplotype_evdience_alloc* hap, overlap_region_alloc* ov } -void lable_large_indels(overlap_region_alloc* overlap_list, All_reads* R_INF, long long read_length, -Correct_dumy* dumy) +void lable_large_indels(overlap_region_alloc* overlap_list, long long read_length, Correct_dumy* dumy, double max_ov_diff_ec) { long long i, j; long long cigar_i, operation, operationLen; @@ -6837,7 +8180,7 @@ Correct_dumy* dumy) { long long window_start, window_end; Window_Pool w_inf; - init_Window_Pool(&w_inf, read_length, WINDOW, (int)(1.0/asm_opt.max_ov_diff_ec)); + init_Window_Pool(&w_inf, read_length, WINDOW, (int)(1.0/max_ov_diff_ec)); int flag = 0; long long realLen = 0, realLen_100 = 0; int to_recover = 0; @@ -6881,20 +8224,23 @@ Correct_dumy* dumy) } - for (i = 0; i < (long long)(overlap_list->length); i++) + for (i = is_delete = 0; i < (long long)(overlap_list->length); i++) { if (overlap_list->list[i].is_match == 1) { overlap_list->list[i].without_large_indel = 1; + is_delete++; } if (overlap_list->list[i].is_match == 100) { overlap_list->list[i].is_match = 1; overlap_list->list[i].without_large_indel = 0; + is_delete++; } } + // if(is_delete) radix_sort_overlap_region_dp_srt(overlap_list->list, overlap_list->list+overlap_list->length); } int debug_print_snp_stat(char* name, haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, All_reads* R_INF) @@ -7280,7 +8626,7 @@ int force_repeat) } -void generate_haplotypes_naive_advance(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list) +void generate_haplotypes_naive_advance(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, void *km) { if(hap->length == 0) return; uint64_t k, l, i, o, *a, ii; @@ -7291,14 +8637,16 @@ void generate_haplotypes_naive_advance(haplotype_evdience_alloc* hap, overlap_re for (k = 1, l = 0; k <= hap->length; ++k) { if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { for (i = l, o = 0; i < k; i++) { - if(hap->list[i].type!=1) continue; + if(hap->list[i].type!=1) continue;///mismatch s = &(hap->snp_stat.a[hap->list[i].overlapSite]); if(s->occ_0 < 2 || s->occ_1 < 2) continue; - if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) o++; + if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) o++;///allels must be real } if(o > 0) { o = ((uint32_t)-1) - o; - o <<= 32; o += l; kv_push(uint64_t, hap->snp_srt, o); + o <<= 32; o += l; + if(!km) kv_push(uint64_t, hap->snp_srt, o); + else kv_push_km(km, uint64_t, hap->snp_srt, o); } l = k; } @@ -7333,11 +8681,9 @@ void generate_haplotypes_naive_advance(haplotype_evdience_alloc* hap, overlap_re for (z = hap->list[i].overlapSite; z >= 0; z--) { t = &(hap->snp_stat.a[z]); if(s->site!=t->site) break; - if(t->occ_0 < 2) fprintf(stderr, "WRONG-CORRECTION\n"); - // if(t->occ_0 == asm_opt.s_hap_cov && t->occ_1 >= asm_opt.infor_cov){ - // fprintf(stderr, "t->id=%u, t->site=%u\n", t->id, t->site); - // } - t->occ_0--; + // t->occ_0--; + t->occ_0 -= hap->list[i].cov; + if(t->occ_0 < 1) fprintf(stderr, "WRONG-CORRECTION\n"); } } } @@ -7373,7 +8719,9 @@ void generate_haplotypes_naive_advance(haplotype_evdience_alloc* hap, overlap_re s = &(hap->snp_stat.a[hap->list[i].overlapSite]); if(s->occ_0 < 2 || s->occ_1 < 2) continue; if(s->score == 1) continue; - o++; kv_push(uint64_t, hap->snp_srt, hap->list[i].overlapSite); + o++; + if(!km) kv_push(uint64_t, hap->snp_srt, hap->list[i].overlapSite); + else kv_push_km(km, uint64_t, hap->snp_srt, hap->list[i].overlapSite); } hap->snp_srt.n -= o; if(o >= 2) {///there are at least two variants at one read @@ -7484,7 +8832,7 @@ void partition_overlaps(overlap_region_alloc* overlap_list, All_reads* R_INF, ///debug_hap_information(overlap_list, R_INF, g_read, hap, dumy); - SetSnpMatrix(hap, &(hap->nn_snp), &(overlap_list->length), 1); + SetSnpMatrix(hap, &(hap->nn_snp), &(overlap_list->length), 1, NULL); uint64_t pre_site = (uint64_t)-1; @@ -7521,7 +8869,7 @@ void partition_overlaps(overlap_region_alloc* overlap_list, All_reads* R_INF, generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat); ///generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat); - lable_large_indels(overlap_list, R_INF, g_read->length, dumy); + lable_large_indels(overlap_list, g_read->length, dumy, asm_opt.max_ov_diff_ec); ///debug_snp_matrix(hap); @@ -7571,7 +8919,7 @@ void partition_overlaps_advance_back(overlap_region_alloc* overlap_list, All_rea ///very time-consuming ///Fix-attention ---> able to be sorted locally // qsort(hap->list, hap->length, sizeof(haplotype_evdience), cmp_haplotype_evdience); - SetSnpMatrix(hap, &(hap->nn_snp), &(overlap_list->length), 1); + SetSnpMatrix(hap, &(hap->nn_snp), &(overlap_list->length), 1, NULL); uint64_t pre_site = (uint64_t)-1; uint64_t num_of_snps = 0; @@ -7606,14 +8954,15 @@ void partition_overlaps_advance_back(overlap_region_alloc* overlap_list, All_rea generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat); ///generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat); - lable_large_indels(overlap_list, R_INF, g_read->length, dumy); + lable_large_indels(overlap_list, g_read->length, dumy, asm_opt.max_ov_diff_ec); } -inline void insert_snp_vv(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, char misBase, UC_Read* g_read) +inline void insert_snp_vv(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, char misBase, UC_Read* g_read, void *km) { if(a_n == 0) return; SnpStats *p = NULL; uint64_t /**nn = 0,**/ i; - kv_pushp(SnpStats, h->snp_stat, &p); + if(!km) kv_pushp(SnpStats, h->snp_stat, &p); + else kv_pushp_km(km, SnpStats, h->snp_stat, &p); p->id = h->snp_stat.n-1; p->occ_0 = 1; p->occ_1 = 0; @@ -7621,45 +8970,34 @@ inline void insert_snp_vv(haplotype_evdience_alloc* h, haplotype_evdience* a, ui p->overlap_num = 0; p->site = a[0].site; p->is_homopolymer = if_is_homopolymer_strict(p->site, g_read->seq, g_read->length); - // nn = h->snp_stat.n*h->overlap; - // if(nn > h->snp_matrix_size) { - // REALLOC(h->snp_matrix, nn); - // memset(h->snp_matrix+h->snp_matrix_size, -1, nn-h->snp_matrix_size); - // h->snp_matrix_size = nn; - // } - // int8_t* vector = Get_SNP_Vector((*h), p->id); for (i = 0; i < a_n; i++) { if(a[i].type == 0) { - // vector[a[i].overlapID] = 0; a[i].overlapSite = p->id; - h->snp_stat.a[p->id].occ_0++; + h->snp_stat.a[p->id].occ_0 += a[i].cov; } else if(a[i].type == 1 && a[i].misBase == misBase) { - // vector[a[i].overlapID] = 1; a[i].overlapSite = p->id; - h->snp_stat.a[p->id].occ_1++; + h->snp_stat.a[p->id].occ_1 += a[i].cov; } else { - // vector[a[i].overlapID] = 2; - // a[i].overlapSite = (uint32_t)-1; - h->snp_stat.a[p->id].occ_2++; + h->snp_stat.a[p->id].occ_2 += a[i].cov; } - h->snp_stat.a[p->id].overlap_num++; + h->snp_stat.a[p->id].overlap_num += a[i].cov; } h->snp_stat.a[p->id].score = -1; } -int insert_snp_ee(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, haplotype_evdience* u_a, UC_Read* g_read) +int insert_snp_ee(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, haplotype_evdience* u_a, UC_Read* g_read, void *km) { uint64_t i, m, occ_0, occ_1[5], occ_2, diff; occ_0 = occ_2 = diff = 0; memset(occ_1, 0, sizeof(uint64_t)*5); for (i = 0; i < a_n; i++) { if(a[i].type == 0){ - occ_0++; + occ_0 += a[i].cov; }else if(a[i].type == 1){ - occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]]++; - diff++; + occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]] += a[i].cov; + diff += a[i].cov; } // else if(a[i].type == 2){ // occ_2++; @@ -7676,7 +9014,7 @@ int insert_snp_ee(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a if(occ_0 == 0 || diff <= 1) return 0; for (i = m = 0; i < 4; i++) { if(occ_1[i] >= 2) { - insert_snp_vv(h, a, a_n, s_H[i], g_read); + insert_snp_vv(h, a, a_n, s_H[i], g_read, km); m++; } } @@ -7733,20 +9071,93 @@ void partition_overlaps_advance(overlap_region_alloc* overlap_list, All_reads* R ///very time-consuming ///Fix-attention ---> able to be sorted locally // qsort(hap->list, hap->length, sizeof(haplotype_evdience), cmp_haplotype_evdience); - SetSnpMatrix(hap, &(hap->nn_snp), &(overlap_list->length), 0); + SetSnpMatrix(hap, &(hap->nn_snp), &(overlap_list->length), 0, NULL); for (k = 1, l = m = 0; k <= hap->length; ++k) { if (k == hap->length || hap->list[k].site != hap->list[l].site) { - m += insert_snp_ee(hap, hap->list+l, k-l, hap->list+m, g_read); + m += insert_snp_ee(hap, hap->list+l, k-l, hap->list+m, g_read, NULL); l = k; } } hap->length = m; - generate_haplotypes_naive_advance(hap, overlap_list); + generate_haplotypes_naive_advance(hap, overlap_list, NULL); // generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat); // generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat); - lable_large_indels(overlap_list, R_INF, g_read->length, dumy); + lable_large_indels(overlap_list, g_read->length, dumy, asm_opt.max_ov_diff_ec); +} + + +void debug_phasing_snp_site_status(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, overlap_region_alloc* olist) +{ + uint64_t i; + for (i = 0; i < a_n; i++) { + if(a[i].site == 76046) { + fprintf(stderr, "[M::utg%.6dl::]x_site->%u, y_site->%u, type->%u, cov->%u, misBase->%c\n", + (int)olist->list[a[i].overlapID].y_id+1, a[i].site, + a[i].overlapSite, a[i].type, a[i].cov, a[i].misBase); + } + } +} + + +void partition_ul_overlaps_advance(overlap_region_alloc* overlap_list, const ul_idx_t *uref, + UC_Read* g_read, UC_Read* overlap_read, Correct_dumy* dumy, haplotype_evdience_alloc* hap, + int force_repeat, double max_ov_diff_ec, long long blockLen, void *km) +{ + ResizeInitHaplotypeEvdience(hap); + + uint64_t k, l, m; + long long window_start, window_end; + long long num_availiable_win = 0; + + Window_Pool w_inf; + init_Window_Pool(&w_inf, g_read->length, blockLen, (int)(1.0/max_ov_diff_ec)); + + int flag = 0; + while(get_Window(&w_inf, &window_start, &window_end) && flag != -2) + { + dumy->length = 0; + dumy->lengthNT = 0; + ///return overlaps that is overlaped with [window_start, window_end] + flag = get_available_interval(window_start, window_end, overlap_list, dumy); + switch (flag) + { + case 1: ///found matched overlaps + break; + case 0: ///do not find any matched overlaps + break; + case -2: ///do not find any matched overlaps, and the next window also cannot match + break; + } + + num_availiable_win = num_availiable_win + dumy->length; + + ///need to deal with + cluster_ul_advance(g_read->seq, window_start, window_end, overlap_list, dumy, uref, hap, overlap_read, 1, w_inf.window_length, km); + } + + + + ///very time-consuming + ///Fix-attention ---> able to be sorted locally + // qsort(hap->list, hap->length, sizeof(haplotype_evdience), cmp_haplotype_evdience); + SetSnpMatrix(hap, &(hap->nn_snp), &(overlap_list->length), 0, km); + for (k = 1, l = m = 0; k <= hap->length; ++k) { + if (k == hap->length || hap->list[k].site != hap->list[l].site) { + // debug_phasing_snp_site_status(hap, hap->list+l, k-l, overlap_list); + m += insert_snp_ee(hap, hap->list+l, k-l, hap->list+m, g_read, km); + l = k; + } + } + hap->length = m; + + generate_haplotypes_naive_advance(hap, overlap_list, km); + + // generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat); + // generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat); + + lable_large_indels(overlap_list, g_read->length, dumy, max_ov_diff_ec); } @@ -7926,7 +9337,7 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, haplotype_evdience_alloc* hap, Round2_alignment* second_round, int force_repeat, int is_consensus, int* fully_cov, int* abnormal) { - clear_Correct_dumy(dumy, overlap_list); + clear_Correct_dumy(dumy, overlap_list, NULL); long long window_start, window_end; @@ -7940,7 +9351,7 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, { dumy->length = 0; dumy->lengthNT = 0; - flag = get_interval(window_start, window_end, overlap_list, dumy); + flag = get_interval(window_start, window_end, overlap_list, dumy, w_inf.window_length); switch (flag) { @@ -7958,7 +9369,7 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, ///now the windows which are larger than WINDOW are verified one-by-one, to improve it, we can do it group-bygroup verify_window(window_start, window_end, overlap_list, dumy, R_INF, g_read->seq); } - + // fprintf(stderr, "###dumy->start_i:%lu, overlap_list->length:%lu\n\n", dumy->start_i, overlap_list->length); // recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read); // partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat); @@ -7974,29 +9385,79 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, (*fully_cov) = check_if_fully_covered(overlap_list, R_INF, g_read, dumy, g, abnormal); } +void debug_phasing_status(overlap_region_alloc *olist, ma_ug_t *ug, uint64_t print_w_list, +haplotype_evdience_alloc* hap, UC_Read* g_read, int64_t flanking, uint64_t yid) +{ + uint64_t i, k, l, ii; + int64_t t; + SnpStats *s = NULL; + + for (i = 0; i < olist->length; i++) { + if(olist->list[i].y_id != yid) continue; + fprintf(stderr, "\n[M::utg%.6d%c::is_match->%u] rev->%u, x->[%u, %u), y->[%u, %u)\n", + (int)olist->list[i].y_id+1, "lc"[ug->u.a[olist->list[i].y_id].circ], olist->list[i].is_match, + olist->list[i].y_pos_strand, olist->list[i].x_pos_s, olist->list[i].x_pos_e+1, olist->list[i].y_pos_s, olist->list[i].y_pos_e+1); + if(print_w_list) { + for (k = 0; k < olist->list[i].w_list_length; k++) { + if(olist->list[i].w_list[k].y_end != -1) { + fprintf(stderr, "x->[%lu, %lu), y->[%d, %d), e->%d\n", + olist->list[i].w_list[k].x_start, olist->list[i].w_list[k].x_end+1, + olist->list[i].w_list[k].y_start, olist->list[i].w_list[k].y_end+1, + olist->list[i].w_list[k].error); + } else { + fprintf(stderr, "x->[-1, -1), y->[-1, -1), e->-1\n"); + } + } + } + } + -void correct_ul_overlap(overlap_region_alloc* overlap_list, const ma_ug_t *uref, + for (k = 1, l = 0; k <= hap->length; ++k) { + if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { + ii = hap->list[l].overlapID; + if(olist->list[ii].y_id != yid) { + l = k; + continue; + } + for (i = l; i < k; i++) { + if(hap->list[i].type!=1) continue; + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->score == 1 && (!(s->occ_0 < 2 || s->occ_1 < 2))) { + fprintf(stderr, "s->site:%u, s->occ_0:%u, s->occ_1:%u, s->occ_2:%u\n", s->site, s->occ_0, s->occ_1, s->occ_2); + for (t = s->site>=flanking?s->site-flanking:0; tlength && t<=s->site+flanking; t++){ + if(t == s->site) fprintf(stderr,"["); + fprintf(stderr,"%c", g_read->seq[t]); + if(t == s->site) fprintf(stderr,"]"); + } + fprintf(stderr,"\n"); + } + } + l = k; + } + } +} + +void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref, UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g, Graph* DAGCon, Cigar_record* current_cigar, haplotype_evdience_alloc* hap, Round2_alignment* second_round, - int force_repeat, int is_consensus, int* fully_cov, int* abnormal, double max_ov_diff_ec) + int force_repeat, int is_consensus, int* fully_cov, int* abnormal, double max_ov_diff_ec, void *km) { - clear_Correct_dumy(dumy, overlap_list); + clear_Correct_dumy(dumy, overlap_list, km); long long window_start, window_end; Window_Pool w_inf; - init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/max_ov_diff_ec)); - + init_Window_Pool(&w_inf, g_read->length, /**WINDOW_UL**//**WINDOW_UL_H**/MIN((((double)THRESHOLD_MAX_SIZE)/max_ov_diff_ec),WINDOW), (int)(1.0/max_ov_diff_ec)); int flag = 0; while(get_Window(&w_inf, &window_start, &window_end) && flag != -2) { dumy->length = 0; dumy->lengthNT = 0; - flag = get_interval(window_start, window_end, overlap_list, dumy); + flag = get_interval(window_start, window_end, overlap_list, dumy, w_inf.window_length); switch (flag) { @@ -8008,18 +9469,24 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ma_ug_t *uref, break; } - ///dumy->lengthNT represent how many overlaps that the length of them is not equal to WINDOW; may larger or less than WINDOW ///dumy->length represent how many overlaps that the length of them is WINDOW ///now the windows which are larger than WINDOW are verified one-by-one, to improve it, we can do it group-bygroup - verify_ul_window(window_start, window_end, overlap_list, dumy, uref, g_read->seq); + verify_ul_window(window_start, window_end, overlap_list, dumy, uref, g_read->seq, max_ov_diff_ec, w_inf.window_length, /**THRESHOLD**/THRESHOLD_MAX_SIZE, km); + // verify_ul_ll_window(window_start, window_end, overlap_list, dumy, uref, g_read->seq, max_ov_diff_ec, w_inf.window_length, km); } - - /** + // recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read); // partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat); - recalcate_window_advance(overlap_list, R_INF, g_read, dumy, overlap_read); - partition_overlaps_advance(overlap_list, R_INF, g_read, overlap_read, dumy, hap, force_repeat); + recalcate_window_ul_advance(overlap_list, uref, g_read, dumy, overlap_read, max_ov_diff_ec, w_inf.window_length, km); + ///after this function, overlap_list is sorted by x_pos_e; used for g_chain + partition_ul_overlaps_advance(overlap_list, uref, g_read, overlap_read, dumy, hap, force_repeat, max_ov_diff_ec, w_inf.window_length, km); + + // debug_phasing_status(overlap_list, uref->ug, 0, hap, g_read, 20, 1176); + // debug_phasing_status(overlap_list, uref->ug, 0, hap, g_read, 20, 1167); + // debug_phasing_status(overlap_list, uref->ug, 0, hap, g_read, 20, 1170); + /** + if(is_consensus) { @@ -8049,6 +9516,12 @@ void init_Cigar_record(Cigar_record* dummy) dummy->current_operation = 127; } +void init_Cigar_record_buf(Cigar_record* dummy, void *km) +{ + memset(dummy, 0, sizeof(*dummy)); + dummy->current_operation = 127; +} + void destory_Cigar_record(Cigar_record* dummy) { @@ -8067,6 +9540,16 @@ void clear_Cigar_record(Cigar_record* dummy) } +void init_Correct_dumy_buf(Correct_dumy* list, void *km) +{ + memset(list, 0, sizeof(Correct_dumy)); + int i; + for (i = 0; i < 256; i++){ + list->Peq_SSE[i] = _mm_setzero_si128(); + } +} + + void init_Correct_dumy(Correct_dumy* list) { list->size = 0; @@ -8094,16 +9577,33 @@ void destory_Correct_dumy(Correct_dumy* list) free(list->corrected_read); } -void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list) +void destory_Correct_dumy_buf(void *km, Correct_dumy* list, int is_z) +{ + kfree(km, list->overlapID); + kfree(km, list->corrected_read); + if(is_z) { + list->size = 0; + list->length = 0; + list->lengthNT = 0; + list->start_i = 0; + list->overlapID = NULL; + list->corrected_read_size = 0; + list->corrected_read_length = 0; + list->corrected_read = NULL; + list->corrected_base = 0; + } +} + +void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list, void *km) { list->length = 0; list->lengthNT = 0; list->start_i = 0; - if (list->size < overlap_list->length) - { + if (list->size < overlap_list->length){ list->size = overlap_list->length; - list->overlapID = (uint64_t*)realloc(list->overlapID, list->size*sizeof(uint64_t)); + if(!km) REALLOC(list->overlapID, list->size); + else KREALLOC(km, list->overlapID, list->size); } list->last_boundary_length = 0; @@ -8264,7 +9764,8 @@ long long* extention_score, long long* q_boundary_score, long long* q_boundary_t long long* t_boundary_score, long long* t_boundary_q_coordinate, long long* droped, int mode) { - /**************for ksw2**************/ + /** + // for ksw2 (*max_t_pos) = (*max_q_pos) = -1; int i, a = sc_mch, b = sc_mis < 0? sc_mis : -sc_mis; // a>0 and b<0 @@ -8297,19 +9798,21 @@ long long* droped, int mode) (*droped) = ez.zdropped; free(ez.cigar); - /** - for (i = 0; i < ez.n_cigar; ++i) // print CIGAR - printf("%d%c", ez.cigar[i]>>4, "MID"[ez.cigar[i]&0xf]); - putchar('\n'); - **/ - /**************for ksw2**************/ + + // for (i = 0; i < ez.n_cigar; ++i) // print CIGAR + // printf("%d%c", ez.cigar[i]>>4, "MID"[ez.cigar[i]&0xf]); + // putchar('\n'); + + // for ksw2 + **/ } int fill_chain_by_affine_gap_debug(Fake_Cigar* chain, char* x_string, char* y_string, overlap_region* ovc, long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n, uint8_t* x_num, uint8_t* y_num, long long* minus_score_thres, long long* final_scores) { - long long i, xOffset, yOffset, xRegionLen, yRegionLen, /**bandLen,**/ maxXpos, maxYpos, zdroped; + /** + long long i, xOffset, yOffset, xRegionLen, yRegionLen, maxXpos, maxYpos, zdroped; long long mapGlobalScore, mapExtentScore; long long xBuoundaryScore, xBuoundaryYcoordinate, yBuoundaryScore, yBuoundaryXcoordinate; ///float band_rate = 0.08; @@ -8348,7 +9851,7 @@ long long* minus_score_thres, long long* final_scores) mode = KSW_EZ_SCORE_ONLY; afine_gap_alignment(x_string, x_num, xRegionLen, y_string, y_num, yRegionLen, c2n, BACKWARD_KSW, MATCH_SCORE_KSW, MISMATCH_SCORE_KSW, GAP_OPEN_KSW, GAP_EXT_KSW, - /**bandLen,**/BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, + BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, &mapExtentScore, &xBuoundaryScore, &xBuoundaryYcoordinate, &yBuoundaryScore, &yBuoundaryXcoordinate, &zdroped, mode); @@ -8434,7 +9937,7 @@ long long* minus_score_thres, long long* final_scores) ///text is x, query is y afine_gap_alignment(x_string+xOffset, x_num, xRegionLen, y_string+yOffset, y_num, yRegionLen, c2n, FORWARD_KSW, MATCH_SCORE_KSW, MISMATCH_SCORE_KSW, GAP_OPEN_KSW, GAP_EXT_KSW, - /**bandLen,**/BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, + BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, &mapExtentScore, &xBuoundaryScore, &xBuoundaryYcoordinate, &yBuoundaryScore, &yBuoundaryXcoordinate, &zdroped, mode); // fprintf(stderr, "# xOffset: %lld, yOffset: %lld, xRegionLen: %lld, yRegionLen: %lld, bandLen: %lld, maxXpos: %lld, maxYpos: %lld, zdroped: %lld\n", @@ -8513,6 +10016,7 @@ long long* minus_score_thres, long long* final_scores) // i, get_fake_gap_pos(chain, i), get_fake_gap_shift(chain, i)); // } } + **/ return 1; } @@ -8523,7 +10027,8 @@ int fill_chain_by_affine_gap(Fake_Cigar* chain, char* x_string, char* y_string, long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n, uint8_t* x_num, uint8_t* y_num, long long* minus_score_thres, long long* final_scores) { - long long i, xOffset, yOffset, xRegionLen, yRegionLen, /**bandLen,**/ maxXpos, maxYpos, zdroped; + /** + long long i, xOffset, yOffset, xRegionLen, yRegionLen, maxXpos, maxYpos, zdroped; long long mapGlobalScore, mapExtentScore; long long xBuoundaryScore, xBuoundaryYcoordinate, yBuoundaryScore, yBuoundaryXcoordinate; ///float band_rate = 0.08; @@ -8576,7 +10081,7 @@ long long* minus_score_thres, long long* final_scores) ///text is x, query is y afine_gap_alignment(x_string+xOffset, x_num, xRegionLen, y_string+yOffset, y_num, yRegionLen, c2n, FORWARD_KSW, MATCH_SCORE_KSW, MISMATCH_SCORE_KSW, GAP_OPEN_KSW, GAP_EXT_KSW, - /**bandLen,**/BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, + BAND_KSW, Z_DROP_KSW, endbouns, &maxXpos, &maxYpos, &mapGlobalScore, &mapExtentScore, &xBuoundaryScore, &xBuoundaryYcoordinate, &yBuoundaryScore, &yBuoundaryXcoordinate, &zdroped, mode); @@ -8603,6 +10108,7 @@ long long* minus_score_thres, long long* final_scores) } (*final_scores) = sum_score; + **/ return 1; } @@ -8630,7 +10136,7 @@ uint8_t* y_num, uint64_t EstimateXOlen, uint64_t EstimateYOlen) x_string = g_read->seq; y_string = overlap_read->seq; - long long sum; + long long sum = 0; if(fill_chain_by_affine_gap(&(ovc->f_cigar), x_string, y_string, ovc, Get_READ_LENGTH(R_INF, ovc->x_id), Get_READ_LENGTH(R_INF, ovc->y_id), NULL, seq_nt6_table, x_num, y_num, &minus_score_thres, &sum) == 0) @@ -8714,7 +10220,7 @@ void recalcate_high_het_overlap(overlap_region_alloc* overlap_list, All_reads* R y_start = total_y_start; Window_Len = x_len + (threshold << 1); - if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF, + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, Get_READ_LENGTH((*R_INF), y_id), &extra_begin, &extra_end, &y_start, &o_len)) { break; @@ -8868,7 +10374,7 @@ void recalcate_high_het_overlap(overlap_region_alloc* overlap_list, All_reads* R ///y_start might be less than 0 y_start = total_y_end - x_len + 1; - if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF, + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, Get_READ_LENGTH((*R_INF), y_id), &extra_begin, &extra_end, &y_start, &o_len)) { break; @@ -9057,7 +10563,7 @@ void recalcate_high_het_overlap(overlap_region_alloc* overlap_list, All_reads* R void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read) { - clear_Correct_dumy(dumy, overlap_list); + clear_Correct_dumy(dumy, overlap_list, NULL); long long window_start, window_end; @@ -9071,7 +10577,7 @@ void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_I { dumy->length = 0; dumy->lengthNT = 0; - flag = get_interval(window_start, window_end, overlap_list, dumy); + flag = get_interval(window_start, window_end, overlap_list, dumy, w_inf.window_length); switch (flag) { diff --git a/Correct.h b/Correct.h index f15c420..9d3c34a 100644 --- a/Correct.h +++ b/Correct.h @@ -6,6 +6,7 @@ #include "POA.h" #include "Process_Read.h" #include "Correct.h" +#include "kalloc.h" //#define CORRECT_THRESHOLD 0.70 #define CORRECT_THRESHOLD 0.60 @@ -132,6 +133,7 @@ typedef struct uint32_t overlapSite; ///there are several types: 0: equal to read 1: not equal to read, but it is a mismatch 2: is a gap uint8_t type; + uint32_t cov; ///misbase char misBase; }haplotype_evdience; @@ -790,10 +792,11 @@ inline int calculate_score(int new_occ_0, int new_occ_1) return consensus; } -inline void SetSnpMatrix(haplotype_evdience_alloc* h, uint32_t *nn_snp, uint64_t *overlap_num, int32_t set_matrix) +inline void SetSnpMatrix(haplotype_evdience_alloc* h, uint32_t *nn_snp, uint64_t *overlap_num, int32_t set_matrix, void *km) { if(nn_snp && overlap_num) { - kv_resize(SnpStats, h->snp_stat, *nn_snp); + if(!km) kv_resize(SnpStats, h->snp_stat, *nn_snp); + else kv_resize_km(km, SnpStats, h->snp_stat, *nn_snp); h->snp_stat.n = 0; h->overlap = *overlap_num; h->core_snp = 0; } @@ -803,13 +806,15 @@ inline void SetSnpMatrix(haplotype_evdience_alloc* h, uint32_t *nn_snp, uint64_t uint64_t new_size = n_snp* n_ovlp; if(h->snp_matrix_size < new_size) { h->snp_matrix_size = new_size; - REALLOC(h->snp_matrix, h->snp_matrix_size); + if(!km) REALLOC(h->snp_matrix, h->snp_matrix_size); + else KREALLOC(km, h->snp_matrix, h->snp_matrix_size); } memset(h->snp_matrix, -1, n_snp * n_ovlp); if(h->r_snp_size < n_ovlp) { h->r_snp_size = n_ovlp; - REALLOC(h->r_snp, h->r_snp_size); + if(!km) REALLOC(h->r_snp, h->r_snp_size); + else KREALLOC(km, h->r_snp, h->r_snp_size); } } } @@ -923,6 +928,18 @@ inline void init_DP_matrix(DP_matrix* dp, uint32_t snp_num) } +inline void InitHaplotypeEvdience_buf(haplotype_evdience_alloc* h, void *km) +{ + memset(h, 0, sizeof(haplotype_evdience_alloc)); + /****************************may have bugs********************************/ + memset(h->flag, 0, WINDOW_MAX_SIZE * sizeof(uint8_t)); + /****************************may have bugs********************************/ + // init_SNP_IDs(&(h->dp.SNP_IDs)); + memset(&(h->dp.SNP_IDs), 0, sizeof(h->dp.SNP_IDs)); + h->dp.SNP_IDs.max_snp_id = -1; +} + + inline void InitHaplotypeEvdience(haplotype_evdience_alloc* h) { @@ -969,7 +986,6 @@ inline void InitHaplotypeEvdience(haplotype_evdience_alloc* h) h->dp.max_buffer = NULL; init_SNP_IDs(&(h->dp.SNP_IDs)); - } inline void StarSubListHaplotypeEvdience(haplotype_evdience_alloc* h) @@ -1001,6 +1017,29 @@ inline void destoryHaplotypeEvdience(haplotype_evdience_alloc* h) destory_SNP_IDs(&(h->dp.SNP_IDs)); } +inline void destoryHaplotypeEvdience_buf(void *km, haplotype_evdience_alloc* h, int is_z) +{ + kfree(km, h->list); + kfree(km, h->snp_stat.a); + kfree(km, h->snp_srt.a); + kfree(km, h->snp_matrix); + kfree(km, h->r_snp); + kfree(km, h->dp.backtrack); + kfree(km, h->dp.max); + kfree(km, h->dp.max_for_sort); + kfree(km, h->dp.visit); + kfree(km, h->dp.backtrack_length); + kfree(km, h->dp.buffer); + kfree(km, h->dp.max_buffer); + kfree(km, h->dp.SNP_IDs.buffer); + kfree(km, h->dp.SNP_IDs.IDs); + if(is_z) { + memset(h, 0, sizeof(*h)); + memset(&(h->dp.SNP_IDs), 0, sizeof(h->dp.SNP_IDs)); + h->dp.SNP_IDs.max_snp_id = -1; + } +} + inline void ResizeInitHaplotypeEvdience(haplotype_evdience_alloc* h) { // h->snp = 0; @@ -1021,17 +1060,14 @@ inline void RsetInitHaplotypeEvdienceFlag(haplotype_evdience_alloc* h, long long /****************************may have bugs********************************/ } -inline void addHaplotypeEvdience(haplotype_evdience_alloc* h, haplotype_evdience* ev) +inline void addHaplotypeEvdience(haplotype_evdience_alloc* h, haplotype_evdience* ev, void *km) { - uint32_t new_length = h->length + 1; - if(new_length > h->size) - { - h->size = h->size * 2; - if(h->size < new_length) - { - h->size = new_length; - } - h->list = (haplotype_evdience*)realloc(h->list, sizeof(haplotype_evdience)*h->size); + if(h->length + 1 > h->size){ + h->size = h->length + 1; + kroundup32(h->size); + if(!km) REALLOC(h->list, h->size); + else KREALLOC(km, h->list, h->size); + // h->list = (haplotype_evdience*)realloc(h->list, sizeof(haplotype_evdience)*h->size); } h->list[h->length] = (*ev); @@ -1080,6 +1116,7 @@ typedef struct } Round2_alignment; +void init_Round2_alignment_buf(Round2_alignment* h, void *km); void init_Round2_alignment(Round2_alignment* h); void destory_Round2_alignment(Round2_alignment* h); void clear_Round2_alignment(Round2_alignment* h); @@ -1091,13 +1128,16 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, Cigar_record* current_cigar, haplotype_evdience_alloc* hap, Round2_alignment* second_round, int force_repeat, int is_consensus, int* fully_cov, int* abnormal); +void init_Correct_dumy_buf(Correct_dumy* list, void *km); void init_Correct_dumy(Correct_dumy* list); void destory_Correct_dumy(Correct_dumy* list); -void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list); +void destory_Correct_dumy_buf(void *km, Correct_dumy* list, int is_z); +void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list, void *km); void clear_Correct_dumy_pure(Correct_dumy* list); void get_seq_from_Graph(Graph* backbone, Graph* DAGCon, Correct_dumy* dumy, Cigar_record* current_cigar, char* self_string, char* r_string, long long r_string_length, long long r_string_site); void init_Cigar_record(Cigar_record* dummy); +void init_Cigar_record_buf(Cigar_record* dummy, void *km); void destory_Cigar_record(Cigar_record* dummy); void clear_Cigar_record(Cigar_record* dummy); void add_new_cell_to_cigar_record(Cigar_record* dummy, uint32_t len, uint32_t type); @@ -1106,11 +1146,11 @@ void add_new_cell_to_cigar_record_with_different_base(Cigar_record* dummy, uint3 void add_existing_cell_to_cigar_record_with_different_base(Cigar_record* dummy, uint32_t len, uint32_t type, char* seq); -void correct_ul_overlap(overlap_region_alloc* overlap_list, const ma_ug_t *uref, +void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref, UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g, Graph* DAGCon, Cigar_record* current_cigar, haplotype_evdience_alloc* hap, Round2_alignment* second_round, - int force_repeat, int is_consensus, int* fully_cov, int* abnormal, double max_ov_diff_ec); + int force_repeat, int is_consensus, int* fully_cov, int* abnormal, double max_ov_diff_ec, void *km); /*** type: diff --git a/Hash_Table.cpp b/Hash_Table.cpp index c1b7dbd..17b3dc5 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -4,6 +4,7 @@ #include #include "Hash_Table.h" #include "ksort.h" +#include "kalloc.h" pthread_mutex_t output_mutex; #define overlap_region_key(a) ((a).y_id) @@ -58,6 +59,18 @@ void destory_overlap_region_alloc(overlap_region_alloc* list) free(list->list); } +void destory_overlap_region_alloc_buf(void *km, overlap_region_alloc* list, int is_z) +{ + uint64_t i = 0; + for (i = 0; i < list->size; i++) { + if(list->list[i].w_list_size>0) kfree(km, list->list[i].w_list); + if(list->list[i].f_cigar.size>0) kfree(km, list->list[i].f_cigar.buffer); + if(list->list[i].boundary_cigars.size>0) kfree(km, list->list[i].boundary_cigars.buffer); + } + kfree(km, list->list); + if(is_z) memset(list, 0, sizeof(*list)); +} + int get_fake_gap_pos(Fake_Cigar* x, int index) { return (x->buffer[index]>>32); @@ -83,15 +96,21 @@ int get_fake_gap_shift(Fake_Cigar* x, int index) } int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, - All_reads* R_INF, int add_beg_end) + long long xLen, long long yLen, int add_beg_end, void *km) { if (list->length + 1 > list->size) { - list->size = list->size * 2; - list->list = (overlap_region*)realloc(list->list, sizeof(overlap_region)*list->size); + uint64_t sl = list->size; + list->size = list->length + 1; + kroundup64(list->size); + if(!km) { + REALLOC(list->list, list->size); + } else { + KREALLOC(km, list->list, list->size); + } /// need to set new space to be 0 - memset(list->list + (list->size/2), 0, sizeof(overlap_region)*(list->size/2)); + memset(list->list + sl, 0, sizeof(overlap_region)*(list->size - sl)); } if (list->length!=0 && list->list[list->length - 1].y_id==tmp->y_id) @@ -122,36 +141,36 @@ int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_regi } - long long x_right_length = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_e - 1; - long long y_right_length = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_e - 1; + long long x_right_length = xLen - tmp->x_pos_e - 1; + long long y_right_length = yLen - tmp->y_pos_e - 1; if(x_right_length <= y_right_length) { - tmp->x_pos_e = Get_READ_LENGTH((*R_INF), tmp->x_id) - 1; + tmp->x_pos_e = xLen - 1; tmp->y_pos_e = tmp->y_pos_e + x_right_length; } else { tmp->x_pos_e = tmp->x_pos_e + y_right_length; - tmp->y_pos_e = Get_READ_LENGTH((*R_INF), tmp->y_id) - 1; + tmp->y_pos_e = yLen - 1; } if (tmp->x_pos_strand == 1) { list->list[list->length].x_id = tmp->x_id; - list->list[list->length].x_pos_e = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_s - 1; - list->list[list->length].x_pos_s = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_e - 1; + list->list[list->length].x_pos_e = xLen - tmp->x_pos_s - 1; + list->list[list->length].x_pos_s = xLen - tmp->x_pos_e - 1; list->list[list->length].x_pos_strand = 0; list->list[list->length].y_id = tmp->y_id; - list->list[list->length].y_pos_e = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_s - 1; - list->list[list->length].y_pos_s = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_e - 1; + list->list[list->length].y_pos_e = yLen - tmp->y_pos_s - 1; + list->list[list->length].y_pos_s = yLen - tmp->y_pos_e - 1; list->list[list->length].y_pos_strand = 1; - resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2)); + resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2), km); if(add_beg_end == 1) { - add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0); + add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0, km); } long long distance_gap; @@ -166,9 +185,8 @@ int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_regi if(distance_gap != pre_distance_gap) { pre_distance_gap = distance_gap; - add_fake_cigar(&(list->list[list->length].f_cigar), - Get_READ_LENGTH((*R_INF), tmp->x_id) - get_fake_gap_pos(&(tmp->f_cigar), i) - 1, - pre_distance_gap); + add_fake_cigar(&(list->list[list->length].f_cigar), xLen - get_fake_gap_pos(&(tmp->f_cigar), i) - 1, + pre_distance_gap, km); } } @@ -178,7 +196,7 @@ int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_regi add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_e, get_fake_gap_shift(&(list->list[list->length].f_cigar), - list->list[list->length].f_cigar.length - 1)); + list->list[list->length].f_cigar.length - 1), km); } } else @@ -195,10 +213,10 @@ int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_regi - resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2)); + resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2), km); if(add_beg_end == 1) { - add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0); + add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0, km); } long long distance_self_pos = tmp->x_pos_e - tmp->x_pos_s; @@ -218,7 +236,7 @@ int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_regi pre_distance_gap = distance_gap; add_fake_cigar(&(list->list[list->length].f_cigar), - get_fake_gap_pos(&(tmp->f_cigar), i), init_distance_gap - pre_distance_gap); + get_fake_gap_pos(&(tmp->f_cigar), i), init_distance_gap - pre_distance_gap, km); } } @@ -228,7 +246,7 @@ int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_regi add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_e, get_fake_gap_shift(&(list->list[list->length].f_cigar), - list->list[list->length].f_cigar.length - 1)); + list->list[list->length].f_cigar.length - 1), km); } } @@ -245,15 +263,20 @@ int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_regi int append_utg_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, - ma_utg_v *ua, int add_beg_end) + ma_utg_v *ua, int add_beg_end, void *km) { - if (list->length + 1 > list->size) { - list->size = list->size * 2; - list->list = (overlap_region*)realloc(list->list, sizeof(overlap_region)*list->size); + uint64_t sl = list->size; + list->size = list->length + 1; + kroundup64(list->size); + if(!km) { + REALLOC(list->list, list->size); + } else { + KREALLOC(km, list->list, list->size); + } /// need to set new space to be 0 - memset(list->list + (list->size/2), 0, sizeof(overlap_region)*(list->size/2)); + memset(list->list + sl, 0, sizeof(overlap_region)*(list->size - sl)); } if (list->length!=0 && list->list[list->length - 1].y_id==tmp->y_id) @@ -310,10 +333,10 @@ int append_utg_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_ list->list[list->length].y_pos_s = ua->a[tmp->y_id].len - tmp->y_pos_e - 1; list->list[list->length].y_pos_strand = 1; - resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2)); + resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2), km); if(add_beg_end == 1) { - add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0); + add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0, km); } long long distance_gap; @@ -330,7 +353,7 @@ int append_utg_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_ pre_distance_gap = distance_gap; add_fake_cigar(&(list->list[list->length].f_cigar), ua->a[tmp->x_id].len - get_fake_gap_pos(&(tmp->f_cigar), i) - 1, - pre_distance_gap); + pre_distance_gap, km); } } @@ -340,7 +363,7 @@ int append_utg_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_ add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_e, get_fake_gap_shift(&(list->list[list->length].f_cigar), - list->list[list->length].f_cigar.length - 1)); + list->list[list->length].f_cigar.length - 1), km); } } else @@ -357,10 +380,10 @@ int append_utg_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_ - resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2)); + resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2), km); if(add_beg_end == 1) { - add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0); + add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0, km); } long long distance_self_pos = tmp->x_pos_e - tmp->x_pos_s; @@ -380,7 +403,7 @@ int append_utg_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_ pre_distance_gap = distance_gap; add_fake_cigar(&(list->list[list->length].f_cigar), - get_fake_gap_pos(&(tmp->f_cigar), i), init_distance_gap - pre_distance_gap); + get_fake_gap_pos(&(tmp->f_cigar), i), init_distance_gap - pre_distance_gap, km); } } @@ -390,7 +413,7 @@ int append_utg_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_ add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_e, get_fake_gap_shift(&(list->list[list->length].f_cigar), - list->list[list->length].f_cigar.length - 1)); + list->list[list->length].f_cigar.length - 1), km); } } @@ -627,7 +650,7 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc ///double band_width_threshold = 0.05; long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, - double band_width_threshold, int max_skip, int x_readLen, int y_readLen) + double band_width_threshold, int max_skip, int x_readLen, int y_readLen, void *km) { long long i, j; long long self_pos, pos, max_j, max_i, max_score, score; @@ -640,7 +663,7 @@ long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* long long total_indels, total_self_length; int32_t ret; - resize_Chain_Data(dp, a_n); + resize_Chain_Data(dp, a_n, km); ret = ha_chain_check(a, a_n, dp, min_score, band_width_threshold); if (ret > 0) { @@ -771,7 +794,7 @@ skip_dp: ///record first site ///the length of f_cigar should be at least 1 ///record the offset of reference - add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap); + add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap, km); long long chainLen = 0; if(result->x_pos_strand == 1) { @@ -784,7 +807,7 @@ skip_dp: { pre_distance_gap = distance_gap; ///record this site - add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap); + add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap, km); } chainLen++; @@ -803,12 +826,12 @@ skip_dp: if(distance_gap == pre_distance_gap) { result->f_cigar.length--; - add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap); + add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap, km); } else { pre_distance_gap = distance_gap; - add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap); + add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap, km); } chainLen++; @@ -823,16 +846,13 @@ skip_dp: void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx, - uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end, overlap_region* f_cigar) + uint64_t readID, uint64_t readLength, All_reads* R_INF, const ul_idx_t *uref, double band_width_threshold, int add_beg_end, overlap_region* f_cigar, void *km) { long long i = 0; uint64_t current_ID; uint64_t current_stand; - if (candidates->length == 0) - { - return; - } + if (candidates->length == 0) return; long long sub_region_beg; long long sub_region_end; @@ -875,32 +895,34 @@ void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_r chain_DP(candidates->list + sub_region_beg, sub_region_end - sub_region_beg + 1, &(candidates->chainDP), f_cigar, band_width_threshold, - 25, Get_READ_LENGTH((*R_INF), (*f_cigar).x_id), Get_READ_LENGTH((*R_INF), (*f_cigar).y_id)); + 25, /**Get_READ_LENGTH((*R_INF), (*f_cigar).x_id)**/readLength, + R_INF?Get_READ_LENGTH((*R_INF), (*f_cigar).y_id):uref->ug->u.a[(*f_cigar).y_id].len, km); ///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1) if ((*f_cigar).x_id != (*f_cigar).y_id) { - append_inexact_overlap_region_alloc(overlap_list, f_cigar, R_INF, add_beg_end); + append_inexact_overlap_region_alloc(overlap_list, f_cigar, readLength, R_INF?Get_READ_LENGTH((*R_INF), (*f_cigar).y_id):uref->ug->u.a[(*f_cigar).y_id].len, add_beg_end, km); } } } void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error, - int extra_begin, int extra_end, int error_threshold) + int extra_begin, int extra_end, int error_threshold, int blockLen, void *km) { long long length = region->x_pos_e - region->x_pos_s + 1; ///the length of window may large or small than WINDOW /****************************may have bugs********************************/ - uint64_t num_windows = length / WINDOW + 4; + uint64_t num_windows = length / blockLen + 4; /****************************may have bugs********************************/ ///w_list_length has alredy set to be 0 at clear_overlap_region_alloc if (num_windows > region->w_list_size) { region->w_list_size = num_windows; - region->w_list = (window_list*)realloc(region->w_list, region->w_list_size*sizeof(window_list)); + if(!km) REALLOC(region->w_list, region->w_list_size); + else KREALLOC(km, region->w_list, region->w_list_size); } @@ -967,17 +989,26 @@ void destory_Chain_Data(Chain_Data* x) free(x->tmp); } -void resize_Chain_Data(Chain_Data* x, long long size) +void resize_Chain_Data(Chain_Data* x, long long size, void *km) { if (size + 1 > x->size) { x->size = size + 1; kroundup64(x->size); - REALLOC(x->score, x->size); - REALLOC(x->pre, x->size); - REALLOC(x->indels, x->size); - REALLOC(x->self_length, x->size); - REALLOC(x->occ, x->size); - REALLOC(x->tmp, x->size); + if(!km) { + REALLOC(x->score, x->size); + REALLOC(x->pre, x->size); + REALLOC(x->indels, x->size); + REALLOC(x->self_length, x->size); + REALLOC(x->occ, x->size); + REALLOC(x->tmp, x->size); + } else { + KREALLOC(km, x->score, x->size); + KREALLOC(km, x->pre, x->size); + KREALLOC(km, x->indels, x->size); + KREALLOC(km, x->self_length, x->size); + KREALLOC(km, x->occ, x->size); + KREALLOC(km, x->tmp, x->size); + } } } @@ -1001,6 +1032,18 @@ void destory_Candidates_list(Candidates_list* l) destory_Chain_Data(&(l->chainDP)); } +void destory_Candidates_list_buf(void *km, Candidates_list* l, int is_z) +{ + kfree(km, l->list); + kfree(km, l->chainDP.score); + kfree(km, l->chainDP.pre); + kfree(km, l->chainDP.indels); + kfree(km, l->chainDP.self_length); + kfree(km, l->chainDP.occ); + kfree(km, l->chainDP.tmp); + if(is_z) memset(l, 0, sizeof(*l)); +} + void init_fake_cigar(Fake_Cigar* x) { x->buffer = NULL; @@ -1021,12 +1064,18 @@ void clear_fake_cigar(Fake_Cigar* x) x->length = 0; } -void add_fake_cigar(Fake_Cigar* x, uint32_t gap_site, int32_t gap_shift) +void add_fake_cigar(Fake_Cigar* x, uint32_t gap_site, int32_t gap_shift, void *km) { if(x->length + 1 > x->size) { - x->size = (x->length + 1) * 2; - x->buffer = (uint64_t*)realloc(x->buffer, sizeof(uint64_t) * x->size); + x->size = x->length + 1; + kroundup32(x->size); + if(!km) { + REALLOC(x->buffer, x->size); + } else { + KREALLOC(km, x->buffer, x->size); + } + // x->buffer = (uint64_t*)realloc(x->buffer, sizeof(uint64_t) * x->size); } x->buffer[x->length] = gap_site; @@ -1051,11 +1100,17 @@ void add_fake_cigar(Fake_Cigar* x, uint32_t gap_site, int32_t gap_shift) } -void resize_fake_cigar(Fake_Cigar* x, uint64_t size) +void resize_fake_cigar(Fake_Cigar* x, uint64_t size, void *km) { if (size > x->size) { x->size = size; - REALLOC(x->buffer, x->size); + if(!km) { + REALLOC(x->buffer, x->size); + } + else { + KREALLOC(km, x->buffer, x->size); + } + } x->length = 0; } @@ -1081,18 +1136,15 @@ void destory_window_list_alloc(window_list_alloc* x) } } -void resize_window_list_alloc(window_list_alloc* x, long long size) +void resize_window_list_alloc(window_list_alloc* x, long long size, void *km) { - if(size > x->size) - { + if(size > x->size){ x->size = size; - x->buffer = (window_list*)realloc(x->buffer, sizeof(window_list) * x->size); + if(!km) REALLOC(x->buffer, x->size); + else KREALLOC(km, x->buffer, x->size); } long long i; - for (i = 0; i < x->size; i++) - { - x->buffer[i].error = -1; - } + for (i = 0; i < x->size; i++) x->buffer[i].error = -1; x->length = 0; } diff --git a/Hash_Table.h b/Hash_Table.h index 2948309..40a541b 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -15,6 +15,14 @@ #define HIGH_HET_OVERLAP_THRESHOLD_FILTER 0.3 #define HIGH_HET_ERROR_RATE 0.08 #define THRESHOLD_MAX_SIZE 31 +#define THRESHOLD_UL_MAX 0.2 +#define WINDOW_UL 75 +#define WINDOW_UL_H 200 +// #define WINDOW_UL_H 150 +#define MIN_UL_ALIN_RATE 0.5 +#define MIN_UL_ALIN_LEN (WINDOW_UL*6) +#define WINDOW_UL_BOUND 48 +#define WINDOW_UL_BOUND_RATE 0.55 #define GROUP_SIZE 4 ///the max cigar likes 10M10D10M10D10M @@ -136,23 +144,25 @@ typedef struct void init_Candidates_list(Candidates_list* l); void clear_Candidates_list(Candidates_list* l); void destory_Candidates_list(Candidates_list* l); +void destory_Candidates_list_buf(void *km, Candidates_list* l, int is_z); void init_overlap_region_alloc(overlap_region_alloc* list); void clear_overlap_region_alloc(overlap_region_alloc* list); void destory_overlap_region_alloc(overlap_region_alloc* list); +void destory_overlap_region_alloc_buf(void *km, overlap_region_alloc* list, int is_z); void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error, -int extra_begin, int extra_end, int error_threshold); +int extra_begin, int extra_end, int error_threshold, int blockLen, void *km); void overlap_region_sort_y_id(overlap_region *a, long long n); void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx, -uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end, overlap_region* f_cigar); +uint64_t readID, uint64_t readLength, All_reads* R_INF, const ul_idx_t *uref, double band_width_threshold, int add_beg_end, overlap_region* f_cigar, void *km); void init_fake_cigar(Fake_Cigar* x); void destory_fake_cigar(Fake_Cigar* x); void clear_fake_cigar(Fake_Cigar* x); -void add_fake_cigar(Fake_Cigar* x, uint32_t gap_site, int32_t gap_shift); -void resize_fake_cigar(Fake_Cigar* x, uint64_t size); +void add_fake_cigar(Fake_Cigar* x, uint32_t gap_site, int32_t gap_shift, void *km); +void resize_fake_cigar(Fake_Cigar* x, uint64_t size, void *km); int get_fake_gap_pos(Fake_Cigar* x, int index); int get_fake_gap_shift(Fake_Cigar* x, int index); @@ -182,12 +192,12 @@ static inline long long y_start_offset(long long x_start, Fake_Cigar* o) return get_fake_gap_shift(o, i - 1); } -void resize_Chain_Data(Chain_Data* x, long long size); +void resize_Chain_Data(Chain_Data* x, long long size, void *km); void init_window_list_alloc(window_list_alloc* x); void clear_window_list_alloc(window_list_alloc* x); void destory_window_list_alloc(window_list_alloc* x); -void resize_window_list_alloc(window_list_alloc* x, long long size); -long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, double band_width_threshold, int max_skip, int x_readLen, int y_readLen); +void resize_window_list_alloc(window_list_alloc* x, long long size, void *km); +long long chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, double band_width_threshold, int max_skip, int x_readLen, int y_readLen, void *km); int append_utg_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, - ma_utg_v *ua, int add_beg_end); + ma_utg_v *ua, int add_beg_end, void *km); #endif diff --git a/Overlaps.cpp b/Overlaps.cpp index 20f0bd2..014e887 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -9757,19 +9757,17 @@ uint64_t *n_utg) void ma_ug_print2(const ma_ug_t *ug, All_reads *RNF, asg_t* read_g, const ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FILE *fp) { - uint8_t* primary_flag = (uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t)); - uint32_t i, j, l; + uint8_t* primary_flag = read_g?(uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t)):NULL; + uint32_t i, j, l, pc = read_g && coverage_cut && sources && ruIndex?1:0; char name[32]; for (i = 0; i < ug->u.n; ++i) { // the Segment lines in GFA ma_utg_t *p = &ug->u.a[i]; if(p->m == 0) continue; sprintf(name, "%s%.6d%c", prefix, i + 1, "lc"[p->circ]); if (print_seq) fprintf(fp, "S\t%s\t%s\tLN:i:%d\trd:i:%u\n", name, p->s? p->s : "*", p->len, - get_ug_coverage(p, read_g, coverage_cut, sources, ruIndex, primary_flag)); + pc?get_ug_coverage(p, read_g, coverage_cut, sources, ruIndex, primary_flag):0); else fprintf(fp, "S\t%s\t*\tLN:i:%d\trd:i:%u\n", name, p->len, - get_ug_coverage(p, read_g, coverage_cut, sources, ruIndex, primary_flag)); - // if (print_seq) fprintf(fp, "S\t%s\t%s\tLN:i:%d\n", name, p->s? p->s : "*", p->len); - // else fprintf(fp, "S\t%s\t*\tLN:i:%d\n", name, p->len); + pc?get_ug_coverage(p, read_g, coverage_cut, sources, ruIndex, primary_flag):0); for (j = l = 0; j < p->n; j++) { if(p->a[j] != (uint64_t)-1) @@ -9779,13 +9777,13 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FIL { fprintf(fp, "A\t%s\t%d\t%c\t%.*s\t%d\t%d\tid:i:%d\tHG:A:%c\n", name, l, "+-"[p->a[j]>>32&1], (int)Get_NAME_LENGTH((*RNF), x), Get_NAME((*RNF), x), - coverage_cut[x].s, coverage_cut[x].e, x, + coverage_cut?coverage_cut[x].s:0, coverage_cut?coverage_cut[x].e:(int)Get_READ_LENGTH((*RNF), x), x, "apmaaa"[((RNF->trio_flag[x]!=FATHER && RNF->trio_flag[x]!=MOTHER)?AMBIGU:RNF->trio_flag[x])]); } else { fprintf(fp, "A\t%s\t%d\t%c\t%s\t%d\t%d\tid:i:%d\tHG:A:%c\n", name, l, "+-"[p->a[j]>>32&1], - "FAKE", coverage_cut[x].s, coverage_cut[x].e, x, '*'); + "FAKE", coverage_cut?coverage_cut[x].s:0, coverage_cut?coverage_cut[x].e:(int)Get_READ_LENGTH((*RNF), x), x, '*'); } } else diff --git a/Overlaps.h b/Overlaps.h index 2bf7342..6cd9deb 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -57,9 +57,14 @@ typedef struct { uint32_t qn, qs, qe; uint32_t tn, ts, te; - uint8_t sec:6, el:1, rev:1; + uint32_t sec:30, el:1, rev:1; } ul_ov_t; +typedef struct { + ul_ov_t *a; + size_t n, m; +} kv_ul_ov_t; + ///query is the read itself typedef struct { @@ -207,6 +212,22 @@ typedef struct { uint32_t utg:31, ori:1, start, len; } utg_intv_t; +typedef struct { + uint32_t *idx; + kvec_t(uint64_t) interval; +} ucov_t; + +typedef struct { + kvec_t(uint64_t) idx; + kvec_t(uint32_t) rids; +} ul_contain; + +typedef struct { + ma_ug_t *ug; + ucov_t *cc; + ul_contain *ct; + // kv_ul_ov_t *ov; +} ul_idx_t; #define MA_HT_INT (-1) #define MA_HT_QCONT (-2) diff --git a/Process_Read.cpp b/Process_Read.cpp index c3a325a..06222aa 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -5,6 +5,7 @@ #include "Process_Read.h" #include "htab.h" #include "Correct.h" +#include "kalloc.h" uint8_t seq_nt6_table[256] = { 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, @@ -1059,7 +1060,7 @@ void retrieve_ul_t(UC_Read* i_r, char *i_s, all_ul_t *ref, uint64_t ID, uint8_t } -void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_t s, int64_t l) +void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_t s, int64_t l, void *km) { if(u->m == 0 || u->n == 0) return; if(l < 0) l = u->len; @@ -1070,7 +1071,9 @@ void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_ i_r->length = l; i_r->RID = 0; if(i_r->length > i_r->size) { i_r->size = i_r->length; - i_r->seq = (char*)realloc(i_r->seq,sizeof(char)*(i_r->size)); + if(!km) REALLOC(i_r->seq, i_r->size); + else KREALLOC(km, i_r->seq, i_r->size); + // i_r->seq = (char*)realloc(i_r->seq,sizeof(char)*(i_r->size)); } r = i_r->seq; } @@ -1108,6 +1111,151 @@ void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_ } } +uint32_t retrieve_u_cov(const ul_idx_t *ul, uint64_t id, uint8_t strand, uint64_t pos, uint8_t dir, int64_t *pi) +{ + uint64_t *a = ul->cc->interval.a + ul->cc->idx[id], cc = 0, ff = 0; + int64_t a_n = ul->cc->idx[id+1]-ul->cc->idx[id], k = 0, cc_i = pi? *pi:0; + if(a_n == 0) return 0; + if(cc_i + 1 >= a_n || cc_i < 0) cc_i = 0; + if(strand) pos = ul->ug->u.a[id].len - pos - 1; + if(dir == 0) { + for (k = cc_i; k + 1 < a_n; k++) { + if(pos>=(a[k]>>32) && pos<(a[k+1]>>32)) { + cc = (uint32_t)a[k]; + ff = 1; + break; + } + } + + if(ff == 0) { + for (k = 0; k < cc_i; k++) { + if(pos>=(a[k]>>32) && pos<(a[k+1]>>32)) { + cc = (uint32_t)a[k]; + ff = 1; + break; + } + } + } + } else { + for (k = cc_i; k >= 0; k--) { + if(pos>=(a[k]>>32) && pos<(a[k+1]>>32)) { + cc = (uint32_t)a[k]; + ff = 1; + break; + } + } + + if(ff == 0) { + for (k = cc_i+1; k + 1 < a_n; k++) { + if(pos>=(a[k]>>32) && pos<(a[k+1]>>32)) { + cc = (uint32_t)a[k]; + ff = 1; + break; + } + } + } + } + + if(pi) *pi = ff?k:0; + return cc; +} + +uint64_t retrieve_u_cov_region(const ul_idx_t *ul, uint64_t id, uint8_t strand, uint64_t s, uint64_t e, int64_t *pi) +{ + uint64_t *a = ul->cc->interval.a + ul->cc->idx[id], cc = 0, o = 0, tk, ts, te, tcc = 0; + int64_t a_n = ul->cc->idx[id+1]-ul->cc->idx[id], k = 0, cc_i = pi? *pi:0; + if(a_n == 0) return 0; + if(cc_i + 1 >= a_n || cc_i < 0) cc_i = 0; + if(strand) { + tk = s; + s = ul->ug->u.a[id].len - e; + e = ul->ug->u.a[id].len - tk; + } + // fprintf(stderr,"\nul->ug->u.a[id].len:%u, fs:%lu, fe:%lu\n", ul->ug->u.a[id].len, a[k]>>32, a[k+1]>>32); + k = cc_i; tk = s; + if(tk < (a[k]>>32)) { + for (; k >= 0; k--) { + if(tk>=(a[k]>>32) && tk<(a[k+1]>>32)) break; + } + } else if(tk >= (a[k+1]>>32)) { + for (; k + 1 < a_n; k++) { + if(tk>=(a[k]>>32) && tk<(a[k+1]>>32)) break; + } + } + if(pi) *pi = k; + + + + for (; k + 1 < a_n; k++) { + ts = a[k]>>32; te = a[k+1]>>32; cc = (uint32_t)a[k]; + o = (MIN(e, te) > MAX(s, ts))?(MIN(e, te)-MAX(s, ts)):0; + tcc += o*cc; + // fprintf(stderr, ">>k:%ld, s:%lu, e:%lu, ts:%lu, te:%lu, o:%lu, cc:%lu\n", k, s, e, ts, te, o, cc); + if(e>=(a[k]>>32) && e<(a[k+1]>>32)) break; + } + + + return tcc; +} + + +uint32_t produce_u_cov(ul_idx_t *ul, uint64_t id, uint8_t strand, uint64_t pos, ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, int64_t gap_fuzz, +uint8_t *sset, kvec_t_u64_warp *buf) +{ + uint64_t k, l, i, z, s = 0, e = 0, qn, tn, qs, qe; + ma_utg_t *u = NULL; + int64_t dp, r; + asg_arc_t t; + + if(strand == (uint8_t)-1 || pos == (uint64_t)-1) { + u = &(ul->ug->u.a[id]); + buf->a.n = 0; kv_resize(uint64_t, buf->a, u->n*2); + for (k = l = 0; k < u->n; k++) { + kv_push(uint64_t, buf->a, l<<1); + kv_push(uint64_t, buf->a, ((l + Get_READ_LENGTH(R_INF, u->a[k]>>33))<<1)|1); + + + i = u->a[k]>>33;///rid + for (z = 0; z < src[i].length; z++) { + if(!src[i].buffer[z].el) continue; + qn = Get_qn(src[i].buffer[z]); tn = Get_tn(src[i].buffer[z]); + if(sset[tn]) continue; + if((Get_qe(src[i].buffer[z]) - Get_qs(src[i].buffer[z])) < min_ovlp) continue; + if((Get_te(src[i].buffer[z]) - Get_ts(src[i].buffer[z])) < min_ovlp) continue; + r = ma_hit2arc(&(src[i].buffer[z]), Get_READ_LENGTH(R_INF, qn), Get_READ_LENGTH(R_INF, tn), + max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + if(r != MA_HT_TCONT) continue;///tn is contained + if(((u->a[k]>>32)&1) == 0) { + qs = Get_qs(src[i].buffer[z]); qe = Get_qe(src[i].buffer[z]); + } else { + qs = (Get_READ_LENGTH(R_INF, i)) - Get_qe(src[i].buffer[z]); + qe = (Get_READ_LENGTH(R_INF, i)) - Get_qs(src[i].buffer[z]); + } + kv_push(uint64_t, buf->a, (l+qs)<<1); + kv_push(uint64_t, buf->a, ((l+qe)<<1)|1); + } + + l += (uint32_t)u->a[k]; + } + sort_kvec_t_u64_warp(buf, 0); + return 0; + } + + if(strand) pos = ul->ug->u.a[id].len - pos - 1; + for (k = 0, dp = 0, s = e = 0; k < buf->a.n; k++) { + e = buf->a.a[k]>>1; + // fprintf(stderr, "[M::%s::k:%lu] [s, e)->[%lu, %lu), dp->%ld\n", __func__, k, s, e, dp); + if(pos >= s && pos < e) break; + s = buf->a.a[k]>>1; + if (buf->a.a[k]&1) --dp; + else ++dp; + } + + return dp; +} + + + void produce_u_seq(char* r, ma_utg_t *u, UC_Read *buf) { if(u->m == 0 || u->n == 0) return; @@ -1147,9 +1295,9 @@ void produce_u_seq(char* r, ma_utg_t *u, UC_Read *buf) } } -void debug_retrieve_rc_sub(all_ul_t *ref, const All_reads *R_INF, ma_utg_v *u, uint32_t n_step) +void debug_retrieve_rc_sub(const ug_opt_t *uopt, all_ul_t *ref, const All_reads *R_INF, ul_idx_t *ul, uint32_t n_step) { - uint64_t i, step, s, e, occ; + uint64_t i, step, s, e, occ, qc, rc; UC_Read f, r; init_UC_Read(&f); init_UC_Read(&r); kvec_t(char) ss; kv_init(ss); @@ -1224,31 +1372,84 @@ void debug_retrieve_rc_sub(all_ul_t *ref, const All_reads *R_INF, ma_utg_v *u, u fprintf(stderr, "[M::%s::# checking: %lu] ==> All_reads\n", __func__, occ); } - if(u) { - for (i = 0, occ = 0; i < u->n; i++) { - kv_resize(char, ss, u->a[i].len); ss.n = u->a[i].len; - retrieve_u_seq(&f, NULL, &(u->a[i]), 0, 0, -1); - produce_u_seq(ss.a, &(u->a[i]), &r); + if(ul) { + + uint8_t *sset = NULL; CALLOC(sset, R_INF->total_reads); + for (i = 0; i < ul->ug->u.n; i++) { + for (s = 0; s < ul->ug->u.a[i].n; s++){ + sset[ul->ug->u.a[i].a[s]>>33] = 1; + } + } + + kvec_t_u64_warp buf; memset(&buf, 0, sizeof(buf)); int64_t pi = 0; + // produce_u_cov(ul, 0, (uint8_t)-1, (uint64_t)-1, &buf); + // produce_u_cov(ul, 0, 0, ul->ug->u.a[0].len>>1, &buf); + + for (i = 0, occ = 0; i < ul->ug->u.n; i++) { + kv_resize(char, ss, ul->ug->u.a[i].len); ss.n = ul->ug->u.a[i].len; + retrieve_u_seq(&f, NULL, &(ul->ug->u.a[i]), 0, 0, -1, NULL); + produce_u_seq(ss.a, &(ul->ug->u.a[i]), &r); if(memcmp(ss.a, f.seq, ss.n)) fprintf(stderr, "4-Wrong whole reverse-read, id: %lu\n", i); - retrieve_u_seq(&r, NULL, &(u->a[i]), 1, 0, -1); + retrieve_u_seq(&r, NULL, &(ul->ug->u.a[i]), 1, 0, -1, NULL); memcpy(ss.a, r.seq, ss.n); reverse_complement(ss.a, ss.n); if(memcmp(ss.a, f.seq, ss.n)) fprintf(stderr, "3-Wrong whole reverse-read, id: %lu\n", i); + produce_u_cov(ul, i, (uint8_t)-1, (uint64_t)-1, uopt->sources, uopt->min_ovlp, uopt->max_hang, uopt->gap_fuzz, sset, &buf); step = ss.n/n_step; if(step <= 0) step = 1; for (s = 0; s < ss.n; s += step) { e = MIN(s+step, ss.n); - retrieve_u_seq(NULL, ss.a, &(u->a[i]), 0, s, e-s); + retrieve_u_seq(NULL, ss.a, &(ul->ug->u.a[i]), 0, s, e-s, NULL); if(memcmp(ss.a, f.seq + s, e - s)) fprintf(stderr, "3-Wrong sub forward-read, id: %lu, [%lu, %lu)\n", i, s, e); - retrieve_u_seq(NULL, ss.a, &(u->a[i]), 1, s, e-s); + retrieve_u_seq(NULL, ss.a, &(ul->ug->u.a[i]), 1, s, e-s, NULL); if(memcmp(ss.a, r.seq + s, e - s)) fprintf(stderr, "3-Wrong sub reverse-read, id: %lu, [%lu, %lu)\n", i, s, e); + + /**if(ul->ug->u.a[i].n > 1)**/ { + rc = produce_u_cov(ul, i, 0, s, uopt->sources, uopt->min_ovlp, uopt->max_hang, uopt->gap_fuzz, sset, &buf); + + qc = retrieve_u_cov(ul, i, 0, s, 0, &pi); + if(rc != qc) fprintf(stderr, "4-Wrong coverage, id: %lu, rc:%lu, qc:%lu, pos:%lu\n", i, rc, qc, s); + pi++; + qc = retrieve_u_cov(ul, i, 0, s, 0, &pi); + if(rc != qc) fprintf(stderr, "4-Wrong coverage, id: %lu, rc:%lu, qc:%lu, pos:%lu\n", i, rc, qc, s); + pi--; + qc = retrieve_u_cov(ul, i, 0, s, 0, &pi); + if(rc != qc) fprintf(stderr, "4-Wrong coverage, id: %lu, rc:%lu, qc:%lu, pos:%lu\n", i, rc, qc, s); + pi++; + qc = retrieve_u_cov(ul, i, 0, s, 1, &pi); + if(rc != qc) fprintf(stderr, "4-Wrong coverage, id: %lu, rc:%lu, qc:%lu, pos:%lu\n", i, rc, qc, s); + pi--; + qc = retrieve_u_cov(ul, i, 0, s, 1, &pi); + if(rc != qc) fprintf(stderr, "4-Wrong coverage, id: %lu, rc:%lu, qc:%lu, pos:%lu\n", i, rc, qc, s); + + + + + rc = produce_u_cov(ul, i, 1, s, uopt->sources, uopt->min_ovlp, uopt->max_hang, uopt->gap_fuzz, sset, &buf); + qc = retrieve_u_cov(ul, i, 1, s, 0, &pi); + if(rc != qc) fprintf(stderr, "4-Wrong coverage, id: %lu, rc:%lu, qc:%lu, pos:%lu\n", i, rc, qc, s); + pi++; + qc = retrieve_u_cov(ul, i, 1, s, 0, &pi); + if(rc != qc) fprintf(stderr, "4-Wrong coverage, id: %lu, rc:%lu, qc:%lu, pos:%lu\n", i, rc, qc, s); + pi--; + qc = retrieve_u_cov(ul, i, 1, s, 0, &pi); + if(rc != qc) fprintf(stderr, "4-Wrong coverage, id: %lu, rc:%lu, qc:%lu, pos:%lu\n", i, rc, qc, s); + pi++; + qc = retrieve_u_cov(ul, i, 1, s, 1, &pi); + if(rc != qc) fprintf(stderr, "4-Wrong coverage, id: %lu, rc:%lu, qc:%lu, pos:%lu\n", i, rc, qc, s); + pi--; + qc = retrieve_u_cov(ul, i, 1, s, 1, &pi); + if(rc != qc) fprintf(stderr, "4-Wrong coverage, id: %lu, rc:%lu, qc:%lu, pos:%lu\n", i, rc, qc, s); + } occ++; } } + + kv_destroy(buf.a); free(sset); fprintf(stderr, "[M::%s::# checking: %lu] ==> ma_utg_v\n", __func__, occ); } diff --git a/Process_Read.h b/Process_Read.h index 0c46d32..30c76c3 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -193,6 +193,7 @@ typedef struct } all_ul_t; extern all_ul_t UL_INF; +// extern uint32_t *het_cnt; void init_All_reads(All_reads* r); void malloc_All_reads(All_reads* r); @@ -216,7 +217,9 @@ void init_all_ul_t(all_ul_t *x, All_reads *hR); void destory_all_ul_t(all_ul_t *x); void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on); void retrieve_ul_t(UC_Read* i_r, char *i_s, all_ul_t *ref, uint64_t ID, uint8_t strand, int64_t s, int64_t l); -void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_t s, int64_t l); -void debug_retrieve_rc_sub(all_ul_t *ref, const All_reads *R_INF, ma_utg_v *u, uint32_t n_step); +void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_t s, int64_t l, void *km); +void debug_retrieve_rc_sub(const ug_opt_t *uopt, all_ul_t *ref, const All_reads *R_INF, ul_idx_t *ul, uint32_t n_step); +uint32_t retrieve_u_cov(const ul_idx_t *ul, uint64_t id, uint8_t strand, uint64_t pos, uint8_t dir, int64_t *pi); +uint64_t retrieve_u_cov_region(const ul_idx_t *ul, uint64_t id, uint8_t strand, uint64_t s, uint64_t e, int64_t *pi); #endif diff --git a/anchor.cpp b/anchor.cpp index d852e53..7bc73c8 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -3,6 +3,7 @@ #include "htab.h" #include "ksort.h" #include "Hash_Table.h" +#include "kalloc.h" #define HA_KMER_GOOD_RATIO 0.333 @@ -51,7 +52,10 @@ struct ha_abufl_s { }; #define HA_ABUF_INIT(HType, MZType, SDType, sf) \ +HType *sf##_init_buf(void *km){HType *b = NULL; KCALLOC((km), b, 1); return b;}\ HType *sf##_init(void){return (HType*)calloc(1, sizeof(HType));}\ +void sf##_free_buf(void *km, HType *ab, int is_z){if(ab){kfree(km, ab->seed); kfree(km, ab->a); kfree(km, ab->mz.a); if((is_z)){memset(ab, 0, sizeof(*ab));}}}\ +void sf##_destroy_buf(void *km, HType *ab){if(ab){kfree(km, ab->seed); kfree(km, ab->a); kfree(km, ab->mz.a); kfree(km, ab);}}\ void sf##_destroy(HType *ab){if(ab){free(ab->seed); free(ab->a); free(ab->mz.a); free(ab);}}\ uint64_t sf##_mem(const HType *ab){\ return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(MZType) + sizeof(SDType)) + sizeof(HType);\ @@ -84,7 +88,7 @@ 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 - mz1_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, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0); + mz1_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, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, NULL); // minimizer of queried read if (ab->mz.m > ab->old_mz_m) { ab->old_mz_m = ab->mz.m; @@ -152,7 +156,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg } cl->length = ab->n_a; - calculate_overlap_region_by_chaining(cl, overlap_list, chain_idx, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain, f_cigar); + calculate_overlap_region_by_chaining(cl, overlap_list, chain_idx, rid, ucr->length, &R_INF, NULL, bw_thres, keep_whole_chain, f_cigar, NULL); #if 0 if (overlap_list->length > 0) { @@ -199,13 +203,12 @@ 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); } -void ha_get_new_ul_candidates(ha_abufl_t *ab, int64_t rid, char* rs, int64_t rl, uint64_t mz_w, uint64_t mz_k, 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* dbg_ct, st_mt_t *sp) +void ha_get_new_ul_candidates(ha_abufl_t *ab, int64_t rid, char* rs, int64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, 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* dbg_ct, st_mt_t *sp, void *km) { uint32_t i; uint64_t k, l; uint32_t high_occ = asm_opt.hom_cov >= 1?asm_opt.hom_cov:1; - // prepare clear_Candidates_list(cl); @@ -213,13 +216,14 @@ void ha_get_new_ul_candidates(ha_abufl_t *ab, int64_t rid, char* rs, int64_t rl, ab->mz.n = 0, ab->n_a = 0; // get the list of anchors - mz2_ha_sketch(rs, rl, mz_w, mz_k, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0); + mz2_ha_sketch(rs, rl, mz_w, mz_k, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, km); // 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); + KREALLOC(km, ab->seed, ab->old_mz_m); } + for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) { int n; ab->seed[i].a = ha_ptl_get(ha_idx, ab->mz.a[i].x, &n); @@ -228,9 +232,9 @@ void ha_get_new_ul_candidates(ha_abufl_t *ab, int64_t rid, char* rs, int64_t rl, } if (ab->n_a > ab->m_a) { ab->m_a = ab->n_a; - kroundup64(ab->m_a); - REALLOC(ab->a, ab->m_a); + KREALLOC(km, ab->a, ab->m_a); } + for (i = 0, k = 0; i < ab->mz.n; ++i) { int j; ///z is one of the minimizer @@ -261,8 +265,9 @@ void ha_get_new_ul_candidates(ha_abufl_t *ab, int64_t rid, char* rs, int64_t rl, // copy over to _cl_ if (ab->m_a >= (uint64_t)cl->size) { cl->size = ab->m_a; - REALLOC(cl->list, cl->size); + KREALLOC(km, cl->list, cl->size); } + for (k = 0; k < ab->n_a; ++k) { k_mer_hit *p = &cl->list[k]; p->readID = ab->a[k].srt >> 33; @@ -279,7 +284,7 @@ void ha_get_new_ul_candidates(ha_abufl_t *ab, int64_t rid, char* rs, int64_t rl, } cl->length = ab->n_a; - calculate_overlap_region_by_chaining(cl, overlap_list, chain_idx, rid, rl, &R_INF, bw_thres, keep_whole_chain, f_cigar); + calculate_overlap_region_by_chaining(cl, overlap_list, chain_idx, rid, rl, NULL, uref, bw_thres, keep_whole_chain, f_cigar, km); #if 0 if (overlap_list->length > 0) { @@ -381,13 +386,13 @@ void calculate_ug_chaining(Candidates_list* candidates, overlap_region_alloc* ov chain_len = chain_DP(candidates->list + sub_region_beg, sub_region_end - sub_region_beg + 1, &(candidates->chainDP), f_cigar, band_width_threshold, - 50, ua->a[(*f_cigar).x_id].len, ua->a[(*f_cigar).y_id].len); + 50, ua->a[(*f_cigar).x_id].len, ua->a[(*f_cigar).y_id].len, NULL); // if ((*f_cigar).x_id != (*f_cigar).y_id) if ((*f_cigar).x_id != (*f_cigar).y_id && chain_len > mz_occ*mz_rate) { - append_utg_inexact_overlap_region_alloc(overlap_list, f_cigar, ua, add_beg_end); + append_utg_inexact_overlap_region_alloc(overlap_list, f_cigar, ua, add_beg_end, NULL); } } } @@ -406,7 +411,7 @@ void ha_get_inter_candidates(ha_abufl_t *ab, uint64_t id, char* r, uint64_t rlen // get the list of anchors mz2_ha_sketch(r, rlen, rw, rk, 0, is_hpc, &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, - NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 1); + NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 1, NULL); // minimizer of queried read if (ab->mz.m > ab->old_mz_m) { @@ -467,7 +472,7 @@ void ha_get_inter_candidates(ha_abufl_t *ab, uint64_t id, char* r, uint64_t rlen } cl->length = ab->n_a; - calculate_overlap_region_by_chaining(cl, ol, chain_idx, id, rlen, &R_INF, bw_thres, keep_whole_chain, f_cigar); + calculate_overlap_region_by_chaining(cl, ol, chain_idx, id, rlen, /**&R_INF**/NULL, NULL, bw_thres, keep_whole_chain, f_cigar, NULL); #if 0 if (ol->length > 0) { @@ -754,13 +759,17 @@ void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overl } -void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, 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, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp) +void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, 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, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, void *km) { extern void *ha_flt_tab; extern ha_pt_t *ha_idx; - ha_get_new_ul_candidates(ab, rid, rs, rl, mz_w, mz_k, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar, dbg_ct, sp); + ha_get_new_ul_candidates(ab, rid, rs, rl, mz_w, mz_k, uref, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar, dbg_ct, sp, km); + if(km) { + ha_abufl_free_buf(km, ab, 1); + destory_Candidates_list_buf(km, cl, 1); + } ks_introsort_or_xs(overlap_list->length, overlap_list->list); } diff --git a/htab.cpp b/htab.cpp index a207807..f447002 100644 --- a/htab.cpp +++ b/htab.cpp @@ -674,7 +674,7 @@ static void sf##_worker_for_mz(void *data, long i, int tid)\ VType *b = &s->mz_buf[tid];\ s->mz_buf[tid].n = 0;\ sf##_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, asm_opt.mz_sample_dist, 0, 0, \ - (s->p->pt&&(s->p->flag&HAF_COUNT_REFINE))?s->p->pt:NULL, s->p->opt->min_rcnt, asm_opt.dp_min_len, asm_opt.dp_e, &(s->mt[tid]), asm_opt.mz_rewin, s->uq);\ + (s->p->pt&&(s->p->flag&HAF_COUNT_REFINE))?s->p->pt:NULL, s->p->opt->min_rcnt, asm_opt.dp_min_len, asm_opt.dp_e, &(s->mt[tid]), asm_opt.mz_rewin, s->uq, NULL);\ s->mz[i].n = s->mz[i].m = b->n;\ MALLOC(s->mz[i].a, b->n);\ MEMCPY(s->mz[i].a, b->a, b->n);\ @@ -733,7 +733,7 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for }\ MALLOC(s->seq[s->n_seq], u->len);\ if(u->s) memcpy(s->seq[s->n_seq], u->s, u->len);\ - else retrieve_u_seq(NULL, s->seq[s->n_seq], u, 0, 0, -1);\ + else retrieve_u_seq(NULL, s->seq[s->n_seq], u, 0, 0, -1, NULL);\ s->len[s->n_seq++] = u->len;\ ++p->n_seq;\ s->sum_len += u->len;\ diff --git a/htab.h b/htab.h index ec9a22c..35dbdf6 100644 --- a/htab.h +++ b/htab.h @@ -94,6 +94,11 @@ 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_buf(void *km); +ha_abufl_t *ha_abufl_init_buf(void *km); +void ha_abuf_destroy_buf(void *km, ha_abuf_t *ab); +void ha_abufl_destroy_buf(void *km, ha_abufl_t *ab); +void ha_abufl_free_buf(void *km, ha_abufl_t *ab, int is_z); 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); @@ -111,8 +116,8 @@ double yak_cpu_usage(void); void ha_triobin(const hifiasm_opt_t *opt); uint32_t *ha_polybin_list(const hifiasm_opt_t *opt); -void mz1_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique); -void mz2_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mzl_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique); +void mz1_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique, void *km); +void mz2_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mzl_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique, void *km); int ha_analyze_count(int n_cnt, int start_cnt, int m_peak_hom, const int64_t *cnt, int *peak_het); int adj_m_peak_hom(int m_peak_hom, int max_i, int max2_i, int max3_i, int *peak_het); void print_hist_lines(int n_cnt, int start_cnt, const int64_t *cnt); diff --git a/inter.cpp b/inter.cpp index d5b4b73..64eee48 100644 --- a/inter.cpp +++ b/inter.cpp @@ -18,8 +18,9 @@ #include "Assembly.h" KSEQ_INIT(gzFile, gzread) -void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, 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, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp); +void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, 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, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, void *km); +#define G_CHAIN_BW 128 #define MG_SEED_IGNORE (1ULL<<41) #define MG_SEED_TANDEM (1ULL<<42) @@ -39,6 +40,15 @@ void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint6 #define generic_key(x) (x) KRADIX_SORT_INIT(gfa64, uint64_t, generic_key, 8) +#define ul_ov_srt_qe_key(p) ((p).qe) +KRADIX_SORT_INIT(ul_ov_srt_qe, ul_ov_t, ul_ov_srt_qe_key, member_size(ul_ov_t, qe)) + + +#define ul_ov_srt_qs_key(p) ((p).qs) +KRADIX_SORT_INIT(ul_ov_srt_qs, ul_ov_t, ul_ov_srt_qs_key, member_size(ul_ov_t, qs)) + + + struct mg_tbuf_s { void *km; int frag_gap; @@ -160,6 +170,7 @@ typedef struct { // global data structure for kt_pipeline() const ma_ug_t *ug; const asg_t *rg; const ug_opt_t *uopt; + const ul_idx_t *uu; kseq_t *ks; int64_t chunk_size; uint64_t n_thread; @@ -283,6 +294,11 @@ KSORT_INIT(sp, sp_node_p, sp_node_lt) KHASH_MAP_INIT_INT(sp, sp_topk_t) KHASH_MAP_INIT_INT(sp2, uint64_t) +typedef struct { + kv_ul_ov_t lo; + kv_ul_ov_t tk; + kvec_t_u64_warp srt; +}glchain_t; typedef struct { // data structure for each step in kt_pipeline() const mg_idxopt_t *opt; @@ -291,6 +307,7 @@ typedef struct { // data structure for each step in kt_pipeline() const ma_ug_t *ug; const asg_t *rg; const ug_opt_t *uopt; + const ul_idx_t *uu; int n, m, sum_len; uint64_t *len, id; char **seq; @@ -299,6 +316,7 @@ typedef struct { // data structure for each step in kt_pipeline() mg_gchains_t **gcs;///useless mg_tbuf_t **buf;///useless ha_ovec_buf_t **hab; + glchain_t *ll; uint64_t num_bases, num_corrected_bases, num_recorrected_bases; } utepdat_t; @@ -1399,7 +1417,7 @@ int32_t mg_gchain1_dp(void *km, const ma_ug_t *ug, const asg_t *rg, int32_t *n_l // extend_lchain(lc, n_lc, qlen, ug); KMALLOC(km, a, n_lc); - ///n_lc how many linear chains + ///n_lc how many linear chains; just filter some linear chains for (i = n_ext = 0; i < n_lc; ++i) { // a[] is a view of frag[]; for sorting mg_lchain_t *r = &lc[i]; gc_frag_t *ai = &a[i]; @@ -2050,7 +2068,7 @@ st_mt_t *sp, mg_tbuf_t *b, int32_t w, int32_t k, int32_t hpc, int32_t mz_sd, int hash = __ac_Wang_hash(hash); mz->n = 0; - mz2_ha_sketch(qseq, qlen, w, k, 0, hpc, mz, ha_flt_tab, mz_sd, NULL, NULL, NULL, -1, -1, -1, sp, mz_rewin, 1); + mz2_ha_sketch(qseq, qlen, w, k, 0, hpc, mz, ha_flt_tab, mz_sd, NULL, NULL, NULL, -1, -1, -1, sp, mz_rewin, 1, NULL); ///a[]->y: weight(8)seg_id(8)flag(8)span(8)pos(32);--->query ///a[]->x: rid(31)rev(1)rpos(33);--->reference a = collect_seed_hits(b->km, opt, 1/**opt->hap_n**/, ha_flt_tab, ha_idx, ug, mz, &n_a, &rep_len, &n_mini_pos, &mini_pos); @@ -2126,28 +2144,349 @@ static void worker_for_ul_alignment(void *data, long i, int tid) // callback for s->opt->is_HPC, asm_opt.mz_sample_dist, asm_opt.mz_rewin, s->opt, s->uopt, &(s->gcs[i])); } +uint32_t overlap_statistics(overlap_region_alloc* olist, ma_ug_t *ug, int64_t *tt, uint8_t mm) +{ + uint32_t k, sp = (uint32_t)-1, ep = (uint32_t)-1, l = 0; + for (k = 0; k < olist->length; k++) { + /** + if(b->olist.list[k].y_id != 38) continue; + **/ + /** + + for (z = 0, te = ta = tua = 0; z < b->olist.list[k].w_list_length; z++) { + if(b->olist.list[k].w_list[z].y_end != -1) { + te += b->olist.list[k].w_list[z].error; + ta += b->olist.list[k].w_list[z].x_end + 1 - b->olist.list[k].w_list[z].x_start; + fprintf(stderr, "x->[%lu, %lu), y->[%d, %d), e->%d\n", b->olist.list[k].w_list[z].x_start, b->olist.list[k].w_list[z].x_end+1, + b->olist.list[k].w_list[z].y_start, b->olist.list[k].w_list[z].y_end+1, b->olist.list[k].w_list[z].error); + } + else { + tua += b->olist.list[k].w_list[z].x_end + 1 - b->olist.list[k].w_list[z].x_start; + } + } + fprintf(stderr, "[M::utg%.6d%c::is_match:%u] x->[%u, %u); y->[%u, %u), ualigned->%u, e_rate->%f\n", b->olist.list[k].y_id+1, "lc"[s->ug->u.a[b->olist.list[k].y_id].circ], + b->olist.list[k].is_match == 1, b->olist.list[k].x_pos_s, b->olist.list[k].x_pos_e+1, b->olist.list[k].y_pos_s, b->olist.list[k].y_pos_e+1, tua, (float)te/(float)ta); + **/ + if(tt){ + uint32_t z; + for (z = 0; z < olist->list[k].w_list_length; z++) { + if(olist->list[k].w_list[z].y_end != -1) { + if(tt) *tt += olist->list[k].w_list[z].x_end+1-olist->list[k].w_list[z].x_start; + } + } + } + if(olist->list[k].is_match == mm) { + if(sp == (uint32_t)-1 || ep < olist->list[k].x_pos_s) { + if(sp != (uint32_t)-1) l += ep + 1 - sp; + sp = olist->list[k].x_pos_s; + ep = olist->list[k].x_pos_e; + } else { + ep = MAX(ep, olist->list[k].x_pos_e); + } + if(ug) { + fprintf(stderr, "[M::utg%.6d%c::is_match->%u] rev->%u, x->[%u, %u), y->[%u, %u)\n", (int)olist->list[k].y_id+1, "lc"[ug->u.a[olist->list[k].y_id].circ], olist->list[k].is_match, + olist->list[k].y_pos_strand, olist->list[k].x_pos_s, olist->list[k].x_pos_e+1, olist->list[k].y_pos_s, olist->list[k].y_pos_e+1); + } + } + } + if(sp != (uint32_t)-1) l += ep + 1 - sp; + return l; +} +/** +void replace_ul(overlap_region_alloc* olist, Correct_dumy* dumy, haplotype_evdience_alloc* hap, const ul_idx_t *uu) +{ + int64_t k, z, n = 0, c_qs, c_qe, c_ts, c_te, c_rev, p_qs, p_qe, p_te, p_ts, p_rev; + uint64_t *sc = NULL, *track = NULL; + overlap_region *c = NULL, *p = NULL; + dumy->length = 0; + for (k = 0; k < olist->length; k++) {///has already sorted by x_pos_e + if(olist->list[k].is_match!=1) continue; + dumy->overlapID[dumy->length] = (uint64_t)-1; + dumy->overlapID[dumy->length] <<= 32; + dumy->overlapID[dumy->length] |= k; + dumy->length++; + } + + kv_resize(uint64_t, hap->snp_srt, dumy->length); + hap->snp_srt.n = dumy->length; + memset(hap->snp_srt.a, 0, hap->snp_srt.n*sizeof(uint64_t)); + + sc = dumy->overlapID; track = hap->snp_srt.a; n = dumy->length; + for (k = 0; k < n; k++) { + c = &(olist->list[(uint32_t)track[k]]); + for (z = k-1; z >= 0; z--) { + p = &(olist->list[(uint32_t)track[z]]); + } + } +} +**/ + + +void gl_chain_gen(overlap_region_alloc* olist, const ul_idx_t *uref, kv_ul_ov_t *res, void *km) +{ + uint64_t k; ul_ov_t *p = NULL; + res->n = 0; + for (k = 0; k < olist->length; k++) { + if(olist->list[k].is_match!=1) continue; + kv_pushp_km(km, ul_ov_t, *res, &p); + p->qn = olist->list[k].x_id; p->qs = olist->list[k].x_pos_s; p->qe = olist->list[k].x_pos_e+1; + p->tn = olist->list[k].y_id; p->el = 1; p->sec = 0; p->rev = olist->list[k].y_pos_strand; + if(p->rev) { + p->ts = uref->ug->u.a[p->tn].len - (olist->list[k].y_pos_e+1); + p->te = uref->ug->u.a[p->tn].len - olist->list[k].y_pos_s; + } else { + p->ts = olist->list[k].y_pos_s; + p->te = olist->list[k].y_pos_e+1; + } + } +} + +int32_t find_ul_ov_max(int32_t n, const ul_ov_t *a, uint32_t x) +{ + int32_t s = 0, e = n; + if (n == 0) return -1; + if (a[n-1].qe < x) return n - 1; + if (a[0].qe >= x) return -1; + while (e > s) { // TODO: finish this block + int32_t m = s + (e - s) / 2; + if (a[m].qe >= x) e = m; + else s = m + 1; + } + assert(s == e); + return s; +} + + + +int64_t get_ecov(const ul_idx_t *uref, ul_ov_t *lv, ul_ov_t *lw, int64_t qlen, int64_t bw, double diff_ec_ul) +{ + int64_t dis_q = lv->qe - lw->qe, dis_t = 0, dif, mm; + uint32_t i, v = ((lv->tn<<1)|lv->rev)^1, w = ((lw->tn<<1)|lw->rev)^1; + const asg_t *g = uref->ug->g; + uint32_t nv = asg_arc_n(g, v); + asg_arc_t *av = asg_arc_a(g, v); + for (i = 0; i < nv; i++) { + if(av[i].del || av[i].v != w) continue; + dis_t = ((uint32_t)av[i].ul); + dis_t -= (lv->rev?lv->ts:g->seq[v>>1].len-lv->te); + break; + } + + dif = (dis_q>dis_t? dis_q-dis_t:dis_t-dis_q); + mm = MAX(dis_q, dis_t); mm *= diff_ec_ul; if(mm < bw) mm = bw; + // if((v>>1) == 1163 && (w>>1) == 1168) fprintf(stderr, ">>>>>>dis_q:%ld, dis_t:%ld, dif:%ld, mm:%ld\n", dis_q, dis_t, dif, mm); + if(dif <= mm) return 1; + return 0; +} + +int64_t gl_exact_chain(kv_ul_ov_t *res, kv_ul_ov_t *ex, const ul_idx_t *uref, int64_t bw, double diff_ec_ul, +int64_t qlen, uint64_t *srt, uint64_t *idx, uint64_t *track, void *km) +{ + // fprintf(stderr, "*****************\n"); + uint32_t li_v, lj_v; + int64_t mm_ovlp, x, i, j, k, sc, csc, mm_sc, mm_idx; + ul_ov_t *li = NULL, *lj = NULL; + const asg_t *g = uref->ug->g; + radix_sort_ul_ov_srt_qe(res->a, res->a + res->n); + for (i = 0; i < (int64_t)res->n; ++i) { + li = &(res->a[i]); li_v = (li->tn<<1)|li->rev; + mm_ovlp = max_ovlp(g, li_v^1); + x = (li->qs + mm_ovlp)*diff_ec_ul; + if(x < bw) x = bw; + x += li->qs + mm_ovlp; + if (x > qlen+1) x = qlen+1; + x = find_ul_ov_max(i, res->a, x); + csc = retrieve_u_cov_region(uref, li->tn, 0, li->ts, li->te, NULL); + mm_sc = csc; mm_idx = -1; + // fprintf(stderr, "---i:%ld, csc:%ld, li->tn:%u, li->ts:%u, li->te:%u\n", i, csc, li->tn, li->ts, li->te); + for (j = x; j >= 0; --j) { // collect potential destination vertices + lj = &(res->a[j]); lj_v = (lj->tn<<1)|lj->rev; + // if(lj->qs >= li->qs) continue; // lj is contained in li on the query coordinate + if(li_v != lj_v && get_ecov(uref, li, lj, qlen, bw, diff_ec_ul)) { + sc = csc + (track[j]>>32); + if(sc > mm_sc) mm_sc = sc, mm_idx = j; + } + } + // 4294967295L + track[i] = mm_sc; track[i] <<= 32; + track[i] |= (mm_idx>=0?mm_idx:((uint64_t)0x7FFFFFFF)); + srt[i] = mm_sc; srt[i] <<= 32; srt[i] |= i; + // fprintf(stderr, "+++i:%ld, mm_idx:%ld, mm_sc:%ld\n", i, mm_idx, mm_sc); + // fprintf(stderr, "[M::utg%.6d%c] qs->%u; qe->%u\n\n", li->tn+1, "lc"[uref->ug->u.a[li->tn].circ], li->qs, li->qe); + } + + int64_t n_v, n_u, n_v0; + radix_sort_gfa64(srt, srt+res->n); ex->n = res->n; + for (k = (int64_t)res->n-1, n_v = n_u = 0; k >= 0; --k) { + // fprintf(stderr, "\nk:%ld\n", k); + n_v0 = n_v; + for (i = (uint32_t)srt[k]; i >= 0 && (track[i]&((uint64_t)0x80000000)) == 0;) { + ex->a[n_v++] = res->a[i]; track[i] |= ((uint64_t)0x80000000); + // fprintf(stderr, "+i:%ld, ", i); + // fprintf(stderr, "[M::utg%.6d%c] qs->%u; qe->%u\n", res->a[i].tn+1, "lc"[uref->ug->u.a[res->a[i].tn].circ], res->a[i].qs, res->a[i].qe); + if((track[i]&((uint64_t)0x7FFFFFFF)) == ((uint64_t)0x7FFFFFFF)) i = -1; + else i = track[i]&((uint64_t)0x7FFFFFFF); + // if(i>=(int64_t)res->n) fprintf(stderr, "ERROR->i:%ld, res->n:%d, n_v:%ld, qlen:%ld\n", i, (int32_t)res->n, n_v, qlen); + // fprintf(stderr, "next_i:%ld\n", i); + // i = (olist->list[i].y_id == (uint32_t)-1?-1:olist->list[i].y_id); + // fprintf(stderr, "-i:%ld\n", i); + } + if(n_v0 == n_v) continue; + ///keep the whole score; do not cut score like minigraph + // sc = (i<0?(srt[k]>>32):((srt[k]>>32)-olist->list[i].x_id)); + sc = srt[k]>>32; + idx[n_u++] = ((uint64_t)sc<<32)|(n_v-n_v0); + } + // if(n_v != (int64_t)res->n) { + // fprintf(stderr, "\nERROR->n_v:%ld, res->n:%d, qlen:%ld\n", n_v, (int32_t)res->n, qlen); + // for (k = 0; k < (int64_t)res->n; k++) { + // fprintf(stderr, "(%ld)srt-sc:%lu, srt-i:%u\n", k, srt[k]>>32, (uint32_t)srt[k]); + // } + + // for (k = 0; k < (int64_t)res->n; k++) { + // fprintf(stderr, "(%ld)track-sc:%lu, track-pi:%lu\n", k, track[k]>>32, track[k]&((uint64_t)0x7FFFFFFF)); + // } + // } + + for (k = 0, n_v = n_v0 = 0; k < n_u; k++) { + n_v0 = n_v; n_v += (uint32_t)idx[k]; + res->a[k].qn = idx[k]>>32; + res->a[k].ts = n_v0; res->a[k].te = n_v; + res->a[k].qs = ex->a[n_v-1].qs; + res->a[k].qe = ex->a[n_v0].qe; + } + res->n = n_u; + return res->n; +} + +uint64_t get_het_site(haplotype_evdience_alloc *hap, uint32_t oid) +{ + uint64_t k, l, i, occ = 0; SnpStats *s = NULL; + for (k = 1, l = 0; k <= hap->length; ++k) { + if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { + if(hap->list[l].overlapID != oid) { + l = k; + continue; + } + for (i = l; i < k; i++) { + if(hap->list[i].type!=1) continue; + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->score == 1 && (!(s->occ_0 < 2 || s->occ_1 < 2))) { + occ++; + } + } + l = k; + } + } + + return (occ&((uint64_t)0x3FFFFFFF)); +} + +int64_t gl_chain_refine(overlap_region_alloc* olist, Correct_dumy* dumy, haplotype_evdience_alloc *hap, glchain_t *ll, const ul_idx_t *uref, double diff_ec_ul, int64_t qlen, void *km) +{ + ll->tk.n = ll->lo.n = 0; + kv_ul_ov_t *idx = &(ll->lo); + gl_chain_gen(olist, uref, idx, km); + if(idx->n == 0) return 0; + kv_resize_km(km, ul_ov_t, ll->tk, idx->n); + kv_resize_km(km, uint64_t, ll->srt.a, idx->n); + kv_resize_km(km, uint64_t, hap->snp_srt, idx->n); + if(gl_exact_chain(idx, &(ll->tk), uref, G_CHAIN_BW, diff_ec_ul, qlen, dumy->overlapID, ll->srt.a.a, hap->snp_srt.a, km)) { + kv_ul_ov_t *chains = &(ll->tk); ul_ov_t *p = NULL; uint64_t k, z, ff, s, e, sft = 50; + radix_sort_ul_ov_srt_qs(idx->a, idx->a + idx->n); + for (k = 0; k < olist->length; k++) { + if(olist->list[k].is_match!=2) continue; + s = olist->list[k].x_pos_s; e = olist->list[k].x_pos_e+1; + for (z = ff = 0; z < idx->n; z++) { + if((s+sft) >= idx->a[z].qs && e <= (idx->a[z].qe+sft)) { + ff = 1; + break; + } + if(idx->a[z].qs >= (e+sft)) break; + } + + if(ff) continue; + kv_pushp_km(km, ul_ov_t, *chains, &p); + p->qn = olist->list[k].x_id; p->qs = olist->list[k].x_pos_s; p->qe = olist->list[k].x_pos_e+1; + p->tn = olist->list[k].y_id; p->el = 1; p->sec = get_het_site(hap, k); + p->rev = olist->list[k].y_pos_strand; + if(p->rev) { + p->ts = uref->ug->u.a[p->tn].len - (olist->list[k].y_pos_e+1); + p->te = uref->ug->u.a[p->tn].len - olist->list[k].y_pos_s; + } else { + p->ts = olist->list[k].y_pos_s; + p->te = olist->list[k].y_pos_e+1; + } + } + } + return 1; +} + +void rescue_contain_reads(overlap_region_alloc* olist) +{ + +} + + static void worker_for_ul_scall_alignment(void *data, long i, int tid) // callback for kt_for() { utepdat_t *s = (utepdat_t*)data; ha_ovec_buf_t *b = s->hab[tid]; + glchain_t *bl = &(s->ll[tid]); int64_t rid = s->id+i; int fully_cov, abnormal; + void *km = s->buf?(s->buf[tid]?s->buf[tid]->km:NULL):NULL; + + + // if (memcmp(UL_INF.nid.a[s->id+i].a, "d0aab024-b3a7-40fb-83cc-22c3d6d951f8", UL_INF.nid.a[s->id+i].n-1)) return; + // fprintf(stderr, "[M::%s::] ==> len: %lu\n", __func__, s->len[i]); + ha_get_ul_candidates_interface(b->abl, rid, s->seq[i], s->len[i], s->opt->w, s->opt->k, s->uu, &b->olist, &b->olist_hp, &b->clist, s->opt->bw_thres, + s->opt->max_n_chain, 1, &(b->k_flag), &b->r_buf, &(b->tmp_region), NULL, &(b->sp), km); - ha_get_ul_candidates_interface(b->abl, rid, s->seq[i], s->len[i], s->opt->w, s->opt->k, &b->olist, &b->olist_hp, &b->clist, s->opt->bw_thres, - s->opt->max_n_chain, 1, &(b->k_flag), &b->r_buf, &(b->tmp_region), NULL, &(b->sp)); - // mg_map_frag(s->ha_flt_tab, s->ha_idx, s->ug, s->rg, s->id+i, s->len[i], s->seq[i], &(s->mzs[tid]), &(s->sps[tid]), s->buf[tid], s->opt->w, s->opt->k, - // s->opt->is_HPC, asm_opt.mz_sample_dist, asm_opt.mz_rewin, s->opt, s->uopt, &(s->gcs[i])); clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2); + // return; + // b->num_correct_base += overlap_statistics(&b->olist, NULL, 0); b->self_read.seq = s->seq[i]; b->self_read.length = s->len[i]; b->self_read.size = 0; - correct_ul_overlap(&b->olist, s->ug, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon, - &b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul); + correct_ul_overlap(&b->olist, s->uu, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon, + &b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul, km); - 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; + // uint64_t k; + // for (k = 0; k < b->olist.length; k++) { + // if(b->olist.list[k].is_match == 1) b->num_correct_base += b->olist.list[k].x_pos_e+1-b->olist.list[k].x_pos_s; + // if(b->olist.list[k].is_match == 2) b->num_recorrect_base += b->olist.list[k].x_pos_e+1-b->olist.list[k].x_pos_s; + // } + + + + gl_chain_refine(&b->olist, &b->correct, &b->hap, bl, s->uu, s->opt->diff_ec_ul, s->len[i], km); + // return; + // 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; memset(&b->self_read, 0, sizeof(b->self_read)); + + uint64_t k; + b->num_correct_base += overlap_statistics(&b->olist, NULL, NULL, 1); + // for (k = 0; k < bl->tk.n; k++) { + // if(bl->tk.a[k].sec == 0) b->num_correct_base += bl->tk.a[k].qe - bl->tk.a[k].qs; + // if(bl->tk.a[k].sec > 0) b->num_recorrect_base += bl->tk.a[k].qe - bl->tk.a[k].qs; + // } + for (k = 0; k < bl->lo.n; k++) { + b->num_read_base += bl->lo.a[k].qe - bl->lo.a[k].qs; + } + + // uint32_t l1 = overlap_statistics(&b->olist, s->uu->ug, 1), l2 = overlap_statistics(&b->olist, s->uu->ug, 2); + // + // if(l1 == 0 && l2 > 0) fprintf(stderr, "[M::%s::%lu::no_match]\n", UL_INF.nid.a[s->id+i].a, s->len[i]); + // fprintf(stderr, "[M::%s::%lu::] l1->%u; l2->%u\n", UL_INF.nid.a[s->id+i].a, s->len[i], l1, l2); + if(km) { + destory_overlap_region_alloc_buf(km, &b->olist, 1); + destory_Correct_dumy_buf(km, &b->correct, 1); + destoryHaplotypeEvdience_buf(km, &b->hap, 1); + } + } void dump_gaf(mg_gres_a *hits, const mg_gchains_t *gs, uint32_t only_p) @@ -2303,7 +2642,7 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac utepdat_t *s; CALLOC(s, 1); s->ha_flt_tab = p->ha_flt_tab; s->ha_idx = p->ha_idx; s->id = p->total_pair; - s->opt = p->opt; s->ug = p->ug; s->uopt = p->uopt; s->rg = p->rg; + s->opt = p->opt; s->uu = p->uu; s->uopt = p->uopt; s->rg = p->rg; while ((ret = kseq_read(p->ks)) >= 0) { if (p->ks->seq.l < (uint64_t)p->opt->k) continue; @@ -2330,7 +2669,15 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac uint64_t i; CALLOC(s->hab, p->n_thread); - for (i = 0; i < p->n_thread; ++i) s->hab[i] = ha_ovec_init(0, 0, 1); + CALLOC(s->ll, p->n_thread); + // CALLOC(s->buf, p->n_thread); + for (i = 0; i < p->n_thread; ++i) { + // s->buf[i] = mg_tbuf_init(); + // s->hab[i] = ha_ovec_buf_init(s->buf[i]->km, 0, 0, 1); + // s->buf[i] = NULL; + // s->hab[i] = ha_ovec_buf_init(NULL, 0, 0, 1); + s->hab[i] = ha_ovec_init(0, 0, 1); + } kt_for(p->n_thread, worker_for_ul_scall_alignment, s, s->n); ///debug /** @@ -2358,10 +2705,12 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac s->num_bases += s->hab[i]->num_read_base; s->num_corrected_bases += s->hab[i]->num_correct_base; s->num_recorrected_bases += s->hab[i]->num_recorrect_base; - ha_ovec_destroy(s->hab[i]); + // mg_tbuf_destroy(s->buf[i]); + ha_ovec_destroy(s->hab[i]); + free(s->ll[i].lo.a); free(s->ll[i].tk.a); free(s->ll[i].srt.a.a); } - free(s->hab); - // free(s->buf); free(s->mzs); free(s->sps); + free(s->hab); free(s->ll); // free(s->buf); + //free(s->mzs); free(s->sps); return s; } else if (step == 2) { // step 3: dump @@ -2914,19 +3263,20 @@ void ul_resolve(ma_ug_t *ug, const asg_t *rg, const ug_opt_t *uopt, int hap_n) uidx_destory(); } -int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const enzyme *fn, void *ha_flt_tab, ha_pt_t *ha_idx, ma_ug_t *ug) +int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const enzyme *fn, void *ha_flt_tab, ha_pt_t *ha_idx, ul_idx_t *uu) { uldat_t sl; memset(&sl, 0, sizeof(sl)); sl.ha_flt_tab = ha_flt_tab; sl.ha_idx = ha_idx; sl.opt = opt; - sl.chunk_size = 100000000; + sl.chunk_size = 500000000; sl.n_thread = asm_opt.thread_num; - sl.ug = ug; + sl.uu = uu; sl.uopt = uopt; scall_ul_pipeline(&sl, fn); + // UL_INF; // print_ul_rs(&UL_INF); - // debug_retrieve_rc_sub(&UL_INF, &R_INF, &(ug->u), 100); + // debug_retrieve_rc_sub(uopt, &UL_INF, &R_INF, (ul_idx_t *)sl.uu, 100); // if(!load_ul_hits(&sl.hits, &sl.nn, asm_opt.output_file_name)) { // scall_ul_pipeline(&sl, fn); // write_ul_hits(&sl.hits, &sl.nn, asm_opt.output_file_name); @@ -2935,10 +3285,171 @@ int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const enzyme *fn, void *ha return 1; } -ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, int64_t gap_fuzz) +void print_dedup_HiFis_seq(ma_ug_t *ug) { - uint64_t i, k, qn, tn, n_read = R_INF.total_reads; + uint64_t i; + ma_utg_t *p = NULL; + for (i = 0; i < ug->u.n; i++) { + p = &(ug->u.a[i]); + CALLOC(p->s, p->len+1); + retrieve_u_seq(NULL, p->s, p, 0, 0, -1, NULL); + p->s[p->len] = '\0'; + } + + FILE* output_file = fopen("dedup_HiFis_seq.gfa", "w"); + ma_ug_print(ug, NULL, NULL, NULL, NULL, "utg", output_file); + fclose(output_file); + + output_file = fopen("dedup_HiFis_seq.noseq.gfa", "w"); + ma_ug_print_simple(ug, NULL, NULL, NULL, NULL, "utg", output_file); + fclose(output_file); + exit(1); +} + +void push_coverage_track(ucov_t *cc, uint64_t uid, ma_utg_t *u, asg_t *rg, ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, int64_t gap_fuzz) +{ + uint64_t k, l, i, z, dp, qn, tn, qs, qe, ori; + int32_t r; asg_arc_t t; + cc->idx[uid] = cc->interval.n; + for (k = l = 0; k < u->n; k++) { + kv_push(uint64_t, cc->interval, l<<1); + kv_push(uint64_t, cc->interval, ((l + Get_READ_LENGTH(R_INF, u->a[k]>>33))<<1)|1); + i = u->a[k]>>33;///rid + for (z = 0; z < src[i].length; z++) { + if(!src[i].buffer[z].el) continue; + qn = Get_qn(src[i].buffer[z]); tn = Get_tn(src[i].buffer[z]); + if(!rg->seq[tn].del) continue; + if((Get_qe(src[i].buffer[z]) - Get_qs(src[i].buffer[z])) < min_ovlp) continue; + if((Get_te(src[i].buffer[z]) - Get_ts(src[i].buffer[z])) < min_ovlp) continue; + r = ma_hit2arc(&(src[i].buffer[z]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + if(r != MA_HT_TCONT) continue;///tn is contained + ori = (u->a[k]>>32)&1; + if(ori == 0) { + qs = Get_qs(src[i].buffer[z]); qe = Get_qe(src[i].buffer[z]); + } else { + qs = (Get_READ_LENGTH(R_INF, i)) - Get_qe(src[i].buffer[z]); + qe = (Get_READ_LENGTH(R_INF, i)) - Get_qs(src[i].buffer[z]); + } + kv_push(uint64_t, cc->interval, (l+qs)<<1); + kv_push(uint64_t, cc->interval, ((l+qe)<<1)|1); + } + l += (uint32_t)u->a[k]; + } + cc->idx[uid+1] = cc->interval.n; + + + radix_sort_gfa64(cc->interval.a+cc->idx[uid], cc->interval.a+cc->interval.n); + for (k = cc->idx[uid], dp = 0; k < cc->interval.n; ++k) { + ///if a[j] is qe + if (cc->interval.a[k]&1) --dp; + else ++dp; + l = cc->interval.a[k]>>1; l <<= 32; l += dp; + cc->interval.a[k] = l; + } +} + + +uint32_t check_if_fully_contain(uint32_t sid, uint32_t lid, uint32_t ori, uint8_t *rset, asg_t *rg, +ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, int64_t gap_fuzz) +{ + uint32_t rid, k, qn, tn, ff = 1; int32_t r; asg_arc_t t; + return 1; + + rid = lid; + for (k = 0; k < src[rid].length; k++) { + if(!src[rid].buffer[k].el) continue; + qn = Get_qn(src[rid].buffer[k]); tn = Get_tn(src[rid].buffer[k]); + if(rg->seq[qn].del || rg->seq[tn].del) continue; + if((Get_qe(src[rid].buffer[k]) - Get_qs(src[rid].buffer[k])) < min_ovlp) continue; + if((Get_te(src[rid].buffer[k]) - Get_ts(src[rid].buffer[k])) < min_ovlp) continue; + r = ma_hit2arc(&(src[rid].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + if(r < 0) continue; + rset[t.v] = ((t.ul>>32)&1)+1; + } + + + + rid = sid; + for (k = 0; k < src[rid].length; k++) { + if(!src[rid].buffer[k].el) continue; + qn = Get_qn(src[rid].buffer[k]); tn = Get_tn(src[rid].buffer[k]); + if(rg->seq[qn].del || rg->seq[tn].del) continue; + if((Get_qe(src[rid].buffer[k]) - Get_qs(src[rid].buffer[k])) < min_ovlp) continue; + if((Get_te(src[rid].buffer[k]) - Get_ts(src[rid].buffer[k])) < min_ovlp) continue; + r = ma_hit2arc(&(src[rid].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + if(r < 0) continue; + if(rset[t.v] != ((((t.ul>>32)&1)^ori)+1)) { + ff = 0; + break; + } + } + + + rid = lid; + for (k = 0; k < src[rid].length; k++) { + if(!src[rid].buffer[k].el) continue; + qn = Get_qn(src[rid].buffer[k]); tn = Get_tn(src[rid].buffer[k]); + if(rg->seq[qn].del || rg->seq[tn].del) continue; + if((Get_qe(src[rid].buffer[k]) - Get_qs(src[rid].buffer[k])) < min_ovlp) continue; + if((Get_te(src[rid].buffer[k]) - Get_ts(src[rid].buffer[k])) < min_ovlp) continue; + r = ma_hit2arc(&(src[rid].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + if(r < 0) continue; + rset[t.v] = 0; + } + + return ff; +} + +ul_contain *ul_contain_gen(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, int64_t gap_fuzz) +{ + uint64_t k, l, i, z, t, qn, tn, ori, qs, qe; + ul_contain *p = NULL; ma_utg_t *u = NULL; + int32_t r; asg_arc_t e; + p->idx.n = p->idx.m = ug->u.n; CALLOC(p->idx.a, p->idx.n); + for (t = 0; t < ug->u.n; t++) { + u = &(ug->u.a[t]); + p->idx.a[t] = p->rids.n; p->idx.a[t] <<= 32; + + for (k = l = 0; k < u->n; k++) { + i = u->a[k]>>33;///rid + + for (z = 0; z < src[i].length; z++) { + if(!src[i].buffer[z].el) continue; + qn = Get_qn(src[i].buffer[z]); tn = Get_tn(src[i].buffer[z]); + if(!rg->seq[tn].del) continue; + if((Get_qe(src[i].buffer[z]) - Get_qs(src[i].buffer[z])) < min_ovlp) continue; + if((Get_te(src[i].buffer[z]) - Get_ts(src[i].buffer[z])) < min_ovlp) continue; + r = ma_hit2arc(&(src[i].buffer[z]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &e); + if(r != MA_HT_TCONT) continue;///tn is contained + + ori = (u->a[k]>>32)&1; + if(ori == 0) { + qs = Get_qs(src[i].buffer[z]); qe = Get_qe(src[i].buffer[z]); + } else { + qs = (Get_READ_LENGTH(R_INF, i)) - Get_qe(src[i].buffer[z]); + qe = (Get_READ_LENGTH(R_INF, i)) - Get_qs(src[i].buffer[z]); + } + qs += l; qe += l; + + kv_push(uint32_t, p->rids, ((tn<<1)|src[i].buffer[z].rev)); + kv_push(uint32_t, p->rids, qs); + kv_push(uint32_t, p->rids, qe); + } + + l += (uint32_t)u->a[k]; + } + p->idx.a[t] |= (p->rids.n - (p->idx.a[t]>>32)); + } + + + return p; +} + +ul_idx_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, int64_t gap_fuzz) +{ + uint64_t i, k, qn, tn, m, n_read = R_INF.total_reads, cc_num = 0; int32_t r; asg_arc_t t, *p = NULL; + uint8_t *rset = NULL; CALLOC(rset, n_read<<1); asg_t *rg = asg_init(); rg->m_seq = rg->n_seq = n_read; MALLOC(rg->seq, rg->m_seq); @@ -2953,17 +3464,18 @@ ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, in if((Get_qe(src[i].buffer[k]) - Get_qs(src[i].buffer[k])) < min_ovlp) continue; if((Get_te(src[i].buffer[k]) - Get_ts(src[i].buffer[k])) < min_ovlp) continue; r = ma_hit2arc(&(src[i].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); - if (r == MA_HT_QCONT) { + if (r == MA_HT_QCONT/** && check_if_fully_contain(qn, tn, src[i].buffer[k].rev, rset, rg, src, min_ovlp, max_hang, gap_fuzz)**/) { rg->seq[qn].del = 1; - } else if(r == MA_HT_TCONT) { + } else if(r == MA_HT_TCONT/** && check_if_fully_contain(tn, qn, src[i].buffer[k].rev, rset, rg, src, min_ovlp, max_hang, gap_fuzz)**/) { rg->seq[tn].del = 1; } if(rg->seq[i].del) break; } } + for (i = 0; i < n_read; i++) { - if(rg->seq[i].del) continue; + if(rg->seq[i].del) {cc_num++; continue;} for (k = 0; k < src[i].length; k++) { if(!src[i].buffer[k].el) continue; qn = Get_qn(src[i].buffer[k]); tn = Get_tn(src[i].buffer[k]); @@ -2980,29 +3492,66 @@ ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, in asg_cleanup(rg); asg_symm(rg); asg_arc_del_trans(rg, gap_fuzz); - ma_ug_t *ug = NULL; + ma_ug_t *ug = NULL; ug = ma_ug_gen(rg); - asg_destroy(rg); + + ul_idx_t *uu = NULL; CALLOC(uu, 1); + uu->ug = ug; CALLOC(uu->cc, 1); + MALLOC(uu->cc->idx, ug->u.n+1); kv_init(uu->cc->interval); - for (i = k = 0; i < ug->u.n; i++) k += ug->u.a[i].len; - fprintf(stderr, "[M::%s::] # unitigs: %lu, # bases: %lu\n", __func__, (uint64_t)ug->u.n, k); - return ug; + for (i = k = m = 0; i < ug->u.n; i++) { + k += ug->u.a[i].len; + push_coverage_track(uu->cc, i, &(ug->u.a[i]), rg, src, min_ovlp, max_hang, gap_fuzz); + } + + // uu->ct = ul_contain_gen(ug, rg, src, min_ovlp, max_hang, gap_fuzz); + + // uu->ov = compress_dedup_HiFis(ug, src); + asg_destroy(rg); free(rset); + + fprintf(stderr, "[M::%s::] # unitigs: %lu, # bases: %lu, # edges: %lu, # cc_num: %lu\n", __func__, (uint64_t)ug->u.n, k, (uint64_t)ug->g->n_arc, cc_num); + // print_dedup_HiFis_seq(ug); + return uu; +} + +void destroy_ul_idx_t(ul_idx_t *uu) +{ + if(!uu) return; + if(uu->cc) { + if(uu->cc) { + free(uu->cc->idx); + free(uu->cc->interval.a); + free(uu->cc); + } + + if(uu->ct) { + free(uu->ct->idx.a); + free(uu->ct->rids.a); + free(uu->ct); + } + // if(uu->ov) { + // free(uu->ov->a); + // free(uu->ov); + // } + } + ma_ug_destroy(uu->ug); + free(uu); } void ul_load(const ug_opt_t *uopt) { fprintf(stderr, "[M::%s::] ==> UL\n", __func__); mg_idxopt_t opt; - ma_ug_t *ug = dedup_HiFis(uopt->sources, uopt->min_ovlp, uopt->max_hang, uopt->gap_fuzz); + ul_idx_t *uu = dedup_HiFis(uopt->sources, uopt->min_ovlp, uopt->max_hang, uopt->gap_fuzz); int cutoff; init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); cutoff = asm_opt.max_n_chain; - init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, 0.05, 0.05); + init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate); int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name) : 0); - if(exist == 0) uidx_l_build(ug, &opt, cutoff); + if(exist == 0) uidx_l_build(uu->ug, &opt, cutoff); if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name); - ul_v_call(&opt, uopt, asm_opt.ar, ha_flt_tab, ha_idx, ug); - ma_ug_destroy(ug); destory_all_ul_t(&UL_INF); + ul_v_call(&opt, uopt, asm_opt.ar, ha_flt_tab, ha_idx, uu); + destroy_ul_idx_t(uu); destory_all_ul_t(&UL_INF); } \ No newline at end of file diff --git a/kalloc.h b/kalloc.h index 5ef432a..b52fda2 100644 --- a/kalloc.h +++ b/kalloc.h @@ -34,6 +34,36 @@ void km_stat(const void *_km, km_stat_t *s); KREALLOC((km), (a), (m)); \ } while (0) +#define kv_resize_km(km, type, v, s) do { \ + if ((v).m < (s)) { \ + (v).m = (s); \ + kv_roundup32((v).m); \ + KREALLOC((km), (v).a, (v).m); \ + } \ + } while (0) + +#define kv_copy_km(km, type, v1, v0) do { \ + if ((v1).m < (v0).n) kv_resize_km((km), type, v1, (v0).n); \ + (v1).n = (v0).n; \ + memcpy((v1).a, (v0).a, sizeof(type) * (v0).n); \ + } while (0) \ + +#define kv_push_km(km, type, v, x) do { \ + if ((v).n == (v).m) { \ + (v).m = (v).m? (v).m<<1 : 2; \ + KREALLOC((km), (v).a, (v).m); \ + } \ + (v).a[(v).n++] = (x); \ + } while (0) + +#define kv_pushp_km(km, type, v, p) do { \ + if ((v).n == (v).m) { \ + (v).m = (v).m? (v).m<<1 : 2; \ + KREALLOC((km), (v).a, (v).m); \ + } \ + *(p) = &(v).a[(v).n++]; \ + } while (0) + #ifndef klib_unused #if (defined __clang__ && __clang_major__ >= 3) || (defined __GNUC__ && __GNUC__ >= 3) #define klib_unused __attribute__ ((__unused__)) diff --git a/sketch.cpp b/sketch.cpp index df6f60e..5ac0e74 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -6,6 +6,7 @@ #include "htab.h" #include "ksort.h" #include "Correct.h" +#include "kalloc.h" #define MAX_HIGH_OCC 8 // TODO: don't hard code if we need to tune this parameter #define MAX_MAX_HIGH_OCC 16 @@ -327,11 +328,11 @@ static void sf##_select_mz_h(VType *p, st_mt_t *mt, int len, int sample_dist, in p->a[n++] = p->a[i];\ p->n = n;\ }\ -void sf##_refine_select(VType *mz, int32_t sidx, int32_t eidx, int32_t sn, int32_t min_freq, st_mt_t *mm, int32_t *rsi, int32_t *rei)\ +void sf##_refine_select(VType *mz, int32_t sidx, int32_t eidx, int32_t sn, int32_t min_freq, st_mt_t *mm, int32_t *rsi, int32_t *rei, void *km)\ {\ int32_t n = sn, m = eidx + 1 - sidx, i, k, t, mk=-1;\ uint64_t ix, kx, ks;\ - kv_resize(uint64_t, *mm, mm->n+n*m);\ + kv_resize_km(km, uint64_t, *mm, mm->n+n*m);\ HType *ma = mz->a + sidx;\ uint64_t *mmt = mm->a + mm->n;\ /**fprintf(stderr, "[M::%s::] ==> +n: %d, m: %d, sn: %d, sidx: %d, eidx: %d\n", __func__, n, m, sn, sidx, eidx);**/\ @@ -372,12 +373,12 @@ void sf##_refine_select(VType *mz, int32_t sidx, int32_t eidx, int32_t sn, int32 if(rsi) (*rsi) = ix + sidx;\ if(rei) (*rei) = kx + sidx;\ }\ -void sf##_refine_sketch(VType *p, ha_pt_t *pt, int32_t rlen, int32_t dp_min_len, float er, int32_t min_freq, st_mt_t *mt)\ +void sf##_refine_sketch(VType *p, ha_pt_t *pt, int32_t rlen, int32_t dp_min_len, float er, int32_t min_freq, st_mt_t *mt, void *km)\ {\ /**fprintf(stderr, "[M::%s::] ==> #########10#########, rlen: %d\n", __func__, rlen);**/\ int32_t i, n = p->n, bd, len = MIN(rlen, dp_min_len), sublen, cnt, ei, li, ri;\ int32_t sn = len*er + 1;\ - kv_resize(uint64_t, *mt, (int64_t)p->n);\ + kv_resize_km(km, uint64_t, *mt, (int64_t)p->n);\ mt->n = p->n; memset(mt->a, 0, sizeof(uint64_t)*p->n);\ for (i = 0; i < n; i++) p->a[i].rid = ha_pt_cnt(pt, p->a[i].x);\ for (i = cnt = 0, bd = -1, ei = -1; i < n; i++){\ @@ -391,7 +392,7 @@ void sf##_refine_sketch(VType *p, ha_pt_t *pt, int32_t rlen, int32_t dp_min_len, }\ }\ /**fprintf(stderr, "[M::%s::] ==> +cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);**/\ - if(cnt >= sn) sf##_refine_select(p, 0, ei, sn, min_freq, mt, NULL, &li);\ + if(cnt >= sn) sf##_refine_select(p, 0, ei, sn, min_freq, mt, NULL, &li, km);\ else{\ li = i-1;\ for (i = 0; i <= li; i++) mt->a[i] = 1;\ @@ -408,7 +409,7 @@ void sf##_refine_sketch(VType *p, ha_pt_t *pt, int32_t rlen, int32_t dp_min_len, }\ }\ /**fprintf(stderr, "[M::%s::] ==> -cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);**/\ - if(cnt >= sn) sf##_refine_select(p, ei, n-1, sn, min_freq, mt, &ri, NULL);\ + if(cnt >= sn) sf##_refine_select(p, ei, n-1, sn, min_freq, mt, &ri, NULL, km);\ else {\ ri = i+1;\ for (i = ri; i <= n-1; i++) mt->a[i] = 1;\ @@ -425,7 +426,7 @@ void sf##_refine_sketch(VType *p, ha_pt_t *pt, int32_t rlen, int32_t dp_min_len, if(cnt >= sn) break;\ }\ }\ - if(cnt >= sn) sf##_refine_select(p, li, ri, sn, min_freq, mt, NULL, NULL);\ + if(cnt >= sn) sf##_refine_select(p, li, ri, sn, min_freq, mt, NULL, NULL, km);\ else for (i = li; i <= ri; i++) mt->a[i] = 1;\ }\ }\ @@ -450,7 +451,7 @@ void sf##_refine_sketch(VType *p, ha_pt_t *pt, int32_t rlen, int32_t dp_min_len, * @param is_hpc homopolymer-compressed or not\ * @param p minimizers\ */\ -void sf##_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, VType *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique)\ +void sf##_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, VType *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique, void *km)\ { /**in default, w = 51, k = 51, is_hpc = 1**/\ extern void *ha_ct_table;\ static const HType dummy = { UINT64_MAX, (((uint64_t)1)< 0 && (int64_t)(len) < (int64_t)((((uint64_t)1)< 0 && w < 256) && (k > 0 && k <= 63));\ if (dbg_ct != NULL) dbg_ct->a.n = 0;\ if (k_flag != NULL) {\ - kv_resize(uint8_t, k_flag->a, (uint64_t)len);\ + kv_resize_km(km, uint8_t, k_flag->a, (uint64_t)len);\ k_flag->a.n = len;\ memset(k_flag->a.a, 0, k_flag->a.n);\ }\ memset(buf, 0xff, w * sizeof(HType));\ memset(&tq, 0, sizeof(tiny_queue_t));\ /**len/w is the evaluated minimizer numbers**/\ - kv_resize(HType, *p, p->n + len/w);\ - kv_resize(uint64_t, *mt, (int64_t)p->m); mt->n = p->n;\ + kv_resize_km(km, HType, *p, p->n + len/w);\ + kv_resize_km(km, uint64_t, *mt, (int64_t)p->m); mt->n = p->n;\ for (i = l = tl = buf_pos = min_pos = 0; i < len; ++i) {\ int c = seq_nt4_table[(uint8_t)str[i]];\ HType info = dummy;\ @@ -511,7 +512,7 @@ void sf##_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is filtered = (cnt == 0);\ cnt = (cnt == 1? 0:cnt);\ }\ - if (dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i));\ + if (dbg_ct != NULL) kv_push_km(km, uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i));\ if (!filtered) info.x = y, info.rid = cnt, info.pos = i, info.rev = z, info.span = kmer_span; /** initially ha_mz1_t::rid keeps the k-mer count**/\ if (k_flag != NULL) k_flag->a.a[i]++;\ if (k_flag != NULL && filtered > 0) k_flag->a.a[i]++;\ @@ -522,12 +523,12 @@ void sf##_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is 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 (sf##_mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos){\ - kv_push(HType, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);\ + kv_push_km(km, HType, *p, buf[j]); kv_push_km(km, uint64_t, *mt, buf_p[j]);\ }\ }\ for (j = 0; j < buf_pos; ++j){\ if (sf##_mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos){\ - kv_push(HType, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);\ + kv_push_km(km, HType, *p, buf[j]); kv_push_km(km, uint64_t, *mt, buf_p[j]);\ }\ }\ }\ @@ -541,12 +542,12 @@ void sf##_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is /**three cases: 1.**/\ if (sf##_mzcmp(&min, &info) >= 0) { /**a new minimum; then write the old min**/\ if (l >= w + k && min.x != UINT64_MAX){\ - kv_push(HType, *p, min); kv_push(uint64_t, *mt, min_s);\ + kv_push_km(km, HType, *p, min); kv_push_km(km, uint64_t, *mt, min_s);\ }\ min = info, min_pos = buf_pos, min_s = buf_p[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(HType, *p, min); kv_push(uint64_t, *mt, min_s);\ + kv_push_km(km, HType, *p, min); kv_push_km(km, uint64_t, *mt, min_s);\ }\ /**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)**/\ @@ -557,22 +558,22 @@ void sf##_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is 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 (sf##_mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos){\ - kv_push(HType, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);\ + kv_push_km(km, HType, *p, buf[j]); kv_push_km(km, uint64_t, *mt, buf_p[j]);\ }\ for (j = 0; j <= buf_pos; ++j)\ if (sf##_mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos){\ - kv_push(HType, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);\ + kv_push_km(km, HType, *p, buf[j]); kv_push_km(km, uint64_t, *mt, buf_p[j]);\ }\ }\ }\ if (++buf_pos == w) buf_pos = 0;\ }\ if (min.x != UINT64_MAX){\ - kv_push(HType, *p, min); kv_push(uint64_t, *mt, min_s);\ + kv_push_km(km, HType, *p, min); kv_push_km(km, uint64_t, *mt, min_s);\ }\ /**debug_pl(str, len, w, k, is_hpc, p, hf, mt);**/\ if (sample_dist > w) sf##_select_mz_h(p, mt, len, sample_dist, ws, k, tl);\ - if (dp_min_len > 0 && pt && mt) sf##_refine_sketch(p, pt, len, dp_min_len, dp_e, min_freq, mt);\ + if (dp_min_len > 0 && pt && mt) sf##_refine_sketch(p, pt, len, dp_min_len, dp_e, min_freq, mt, km);\ for (i = 0; i < (int)p->n; ++i) /**populate .rid as this was keeping counts**/\ p->a[i].rid = rid;\ }