From 06556e4a9b6899248c595901f8c5e058a5d404cb Mon Sep 17 00:00:00 2001 From: Haoyu Cheng Date: Tue, 17 Dec 2019 16:58:20 -0500 Subject: [PATCH] backup without bundaries processing --- Assembly.cpp | 699 +++++++++++++------------- Correct.cpp | 1306 ++++++++++++++++++++++++++++++++++++++++++++++-- Correct.h | 29 ++ Hash_Table.cpp | 521 ++++++++++++++++++- Hash_Table.h | 46 +- Overlaps.cpp | 2 +- kmer.h | 18 +- 7 files changed, 2227 insertions(+), 394 deletions(-) diff --git a/Assembly.cpp b/Assembly.cpp index addab11..e6c56b3 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -578,15 +578,15 @@ void Build_hash_table_multiple_thr() destory_kseq(); destory_R_buffer(); ///destory_Total_Count_Table(&TCB); - if (write_index_to_disk) - { - write_Total_Pos_Table(&PCB, read_file_name); - ///destory_Total_Pos_Table(&PCB); - ///load_Total_Pos_Table(&PCB, read_file_name); - write_All_reads(&R_INF, read_file_name); - ///destory_All_reads(&R_INF); - ///load_All_reads(&R_INF, read_file_name); - } + // if (write_index_to_disk) + // { + // write_Total_Pos_Table(&PCB, read_file_name); + // ///destory_Total_Pos_Table(&PCB); + // ///load_Total_Pos_Table(&PCB, read_file_name); + // write_All_reads(&R_INF, read_file_name); + // ///destory_All_reads(&R_INF); + // ///load_All_reads(&R_INF, read_file_name); + // } } ///destory_Total_Count_Table(&TCB); @@ -1157,29 +1157,44 @@ inline void push_cigar(Compressed_Cigar_record* records, long long ID, Cigar_rec } -void push_overlaps(ma_hit_t_alloc* paf, overlap_region_alloc* overlap_list, int flag) +void push_overlaps(ma_hit_t_alloc* paf, overlap_region_alloc* overlap_list, int flag, +All_reads* R_INF, int if_reverse) { - long long i = 0; + long long i = 0, xLen, yLen; ma_hit_t tmp; clear_ma_hit_t_alloc(paf); for (i = 0; i < overlap_list->length; i++) { if (overlap_list->list[i].is_match == flag) { + xLen = Get_READ_LENGTH((*R_INF), overlap_list->list[i].x_id); + yLen = Get_READ_LENGTH((*R_INF), overlap_list->list[i].y_id); + tmp.qns = overlap_list->list[i].x_id; tmp.qns = tmp.qns << 32; - tmp.qns = tmp.qns | (uint64_t)(overlap_list->list[i].x_pos_s); - - tmp.qe = overlap_list->list[i].x_pos_e; - tmp.tn = overlap_list->list[i].y_id; - tmp.ts = overlap_list->list[i].y_pos_s; - tmp.te = overlap_list->list[i].y_pos_e; - + + if(if_reverse != 0) + { + tmp.qns = tmp.qns | (uint64_t)(xLen - overlap_list->list[i].x_pos_s - 1); + tmp.qe = xLen - overlap_list->list[i].x_pos_e - 1; + tmp.ts = yLen - overlap_list->list[i].y_pos_s - 1; + tmp.te = yLen - overlap_list->list[i].y_pos_e - 1; + } + else + { + tmp.qns = tmp.qns | (uint64_t)(overlap_list->list[i].x_pos_s); + tmp.qe = overlap_list->list[i].x_pos_e; + tmp.ts = overlap_list->list[i].y_pos_s; + tmp.te = overlap_list->list[i].y_pos_e; + } + + ///for overlap_list, the x_strand of all overlaps are 0, so the tmp.rev is the same as the y_strand tmp.rev = overlap_list->list[i].y_pos_strand; - tmp.bl = R_INF.read_length[overlap_list->list[i].y_id]; + ///tmp.bl = R_INF.read_length[overlap_list->list[i].y_id]; + tmp.bl = Get_READ_LENGTH((*R_INF), overlap_list->list[i].y_id); tmp.ml = overlap_list->list[i].strong; tmp.no_l_indel = overlap_list->list[i].without_large_indel; @@ -1238,7 +1253,7 @@ long long xBeg, long long xEnd, long long yBeg, long long yEnd) } long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, -overlap_region_alloc* overlap_list, UC_Read* x_read, UC_Read* y_read, int flag, int test_exact) +overlap_region_alloc* overlap_list, int flag) { long long i = 0; long long available_overlaps = 0; @@ -1286,23 +1301,7 @@ overlap_region_alloc* overlap_list, UC_Read* x_read, UC_Read* y_read, int flag, tmp.ml = overlap_list->list[i].strong; tmp.no_l_indel = overlap_list->list[i].without_large_indel; - - - if(test_exact == 1) - { - if(overlap_list->list[i].y_pos_strand == 0) - { - recover_UC_Read(y_read, &R_INF, overlap_list->list[i].y_id); - } - else - { - recover_UC_Read_RC(y_read, &R_INF, overlap_list->list[i].y_id); - } - - tmp.el = if_exact_match(x_read->seq, x_read->length, y_read->seq, y_read->length, - overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e, - overlap_list->list[i].y_pos_s, overlap_list->list[i].y_pos_e); - } + tmp.el = overlap_list->list[i].shared_seed; add_ma_hit_t_alloc(paf, &tmp); } @@ -1736,7 +1735,7 @@ HeapSq* heap, Candidates_list* l, small_hash_table* forward, small_hash_table* r } void get_new_candidates(long long readID, UC_Read* g_read, overlap_region_alloc* overlap_list, k_mer_pos_list_alloc* array_list, -HeapSq* heap, Candidates_list* l, double band_width_threshold) +HeapSq* heap, Candidates_list* l, double band_width_threshold, int keep_whole_chain) { HPC_seq HPC_read; Hash_code k_code; @@ -1821,7 +1820,8 @@ HeapSq* heap, Candidates_list* l, double band_width_threshold) ///以x_pos_e，即结束位置为主元排序 ///calculate_overlap_region(l, overlap_list, readID, g_read->length, &R_INF); - calculate_overlap_region_by_chaining(l, overlap_list, readID, g_read->length, &R_INF, band_width_threshold); + calculate_overlap_region_by_chaining(l, overlap_list, readID, g_read->length, &R_INF, + band_width_threshold, keep_whole_chain); } @@ -1910,7 +1910,6 @@ void* Overlap_calculate_heap_merge(void* arg) init_small_hash_table(&reverse); - uint8_t c2n[256]; memset(c2n, 4, 256); c2n['A'] = c2n['a'] = 0; c2n['C'] = c2n['c'] = 1; @@ -1920,7 +1919,7 @@ void* Overlap_calculate_heap_merge(void* arg) for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num) { ///get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); - get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02); + get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02, 1); clear_Cigar_record(¤t_cigar); clear_Round2_alignment(&second_round); @@ -1955,8 +1954,8 @@ void* Overlap_calculate_heap_merge(void* arg) - push_overlaps(&(R_INF.paf[i]), &overlap_list, 1); - push_overlaps(&(R_INF.reverse_paf[i]), &overlap_list, 2); + push_overlaps(&(R_INF.paf[i]), &overlap_list, 1, &R_INF, roundID%2); + push_overlaps(&(R_INF.reverse_paf[i]), &overlap_list, 2, &R_INF, roundID%2); } @@ -1971,17 +1970,11 @@ void* Overlap_calculate_heap_merge(void* arg) - /** - fprintf(stderr, "candidate_overlap_reads: %llu\n", candidate_overlap_reads); - fprintf(stderr, "total_shared_seed: %llu\n", total_shared_seed); - **/ destory_Candidates_list(&l); destory_overlap_region_alloc(&overlap_list); - //destory_Candidates_list(&debug_l); destory_Heap(&heap); destory_k_mer_pos_list_alloc(&array_list); - ///destory_k_mer_pos_list_alloc_prefilter(&array_list); destory_Graph(&POA_Graph); destory_Graph(&DAGCon); @@ -2121,7 +2114,7 @@ void* Output_related_reads(void* arg) memcmp(required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { ////get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); - get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02); + get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02, 1); fprintf(stderr, ">%.*s\n", Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i)); @@ -2893,21 +2886,303 @@ void debug_print_overlap(char* y_name, overlap_region_alloc* overlap_list, All_r } } +int debug_diff(int a, int b) +{ + if(DIFF(a, b)!= abs(a-b)) + { + fprintf(stderr, "sbsbsbsb\n"); + } +} + +void update_overlaps(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf, +UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact) +{ + + long long inner_j = 0; + long long j = 0; + long long x_overlapLen, y_overlapLen; + while (j < overlap_list->length && inner_j < paf->length) + { + if(overlap_list->list[j].y_id < paf->buffer[inner_j].tn) + { + j++; + } + else if(overlap_list->list[j].y_id > paf->buffer[inner_j].tn) + { + inner_j++; + } + else + { + if(overlap_list->list[j].y_pos_strand == paf->buffer[inner_j].rev) + { + x_overlapLen = Get_qe(paf->buffer[inner_j]) - Get_qs(paf->buffer[inner_j]) + 1; + y_overlapLen = Get_te(paf->buffer[inner_j]) - Get_ts(paf->buffer[inner_j]) + 1; + if(x_overlapLen < y_overlapLen) x_overlapLen = y_overlapLen; + x_overlapLen = x_overlapLen * 0.1; + + // debug_diff(overlap_list->list[j].x_pos_s, Get_qs(paf->buffer[inner_j])); + // debug_diff(overlap_list->list[j].x_pos_e, Get_qe(paf->buffer[inner_j])); + // debug_diff(overlap_list->list[j].y_pos_s, Get_ts(paf->buffer[inner_j])); + // debug_diff(overlap_list->list[j].y_pos_e, Get_te(paf->buffer[inner_j])); + + ///fprintf(stderr, "hehe\n"); + + // if( + // ((DIFF(overlap_list->list[j].x_pos_s, Get_qs(paf->buffer[inner_j])) < x_overlapLen) + // && (DIFF(overlap_list->list[j].x_pos_e, Get_qe(paf->buffer[inner_j])) < x_overlapLen)) + // || + // ((DIFF(overlap_list->list[j].y_pos_s, Get_ts(paf->buffer[inner_j])) < x_overlapLen) + // && (DIFF(overlap_list->list[j].y_pos_e, Get_te(paf->buffer[inner_j])) < x_overlapLen))) + if( + ((DIFF(overlap_list->list[j].x_pos_s, Get_qs(paf->buffer[inner_j])) < x_overlapLen) + && (DIFF(overlap_list->list[j].x_pos_e, Get_qe(paf->buffer[inner_j])) < x_overlapLen)) + || + ((DIFF(overlap_list->list[j].y_pos_s, Get_ts(paf->buffer[inner_j])) < x_overlapLen) + && (DIFF(overlap_list->list[j].y_pos_e, Get_te(paf->buffer[inner_j])) < x_overlapLen)) + ) + { + overlap_list->list[j].is_match = is_match; + overlap_list->list[j].strong = paf->buffer[inner_j].ml; + overlap_list->list[j].without_large_indel = paf->buffer[inner_j].no_l_indel; + if(is_exact == 1) + { + if(overlap_list->list[j].y_pos_strand == 0) + { + recover_UC_Read(overlap_read, &R_INF, overlap_list->list[j].y_id); + } + else + { + recover_UC_Read_RC(overlap_read, &R_INF, overlap_list->list[j].y_id); + } + if(if_exact_match(g_read->seq, g_read->length, overlap_read->seq, overlap_read->length, + overlap_list->list[j].x_pos_s, overlap_list->list[j].x_pos_e, + overlap_list->list[j].y_pos_s, overlap_list->list[j].y_pos_e)) + { + overlap_list->list[j].shared_seed = 1; + } + else + { + overlap_list->list[j].shared_seed = 0; + } + } + } + else + { + overlap_list->list[j].is_match = 3; + } + } + else + { + overlap_list->list[j].is_match = 3; + } + + j++; + inner_j++; + } + } +} + +void update_exact_overlaps(overlap_region_alloc* overlap_list, UC_Read* g_read, UC_Read* overlap_read) +{ + long long j; + for (j = 0; j < overlap_list->length; j++) + { + if (overlap_list->list[j].is_match != 1) + { + if(overlap_list->list[j].y_pos_strand == 0) + { + recover_UC_Read(overlap_read, &R_INF, overlap_list->list[j].y_id); + } + else + { + recover_UC_Read_RC(overlap_read, &R_INF, overlap_list->list[j].y_id); + } + + if(if_exact_match(g_read->seq, g_read->length, overlap_read->seq, overlap_read->length, + overlap_list->list[j].x_pos_s, overlap_list->list[j].x_pos_e, + overlap_list->list[j].y_pos_s, overlap_list->list[j].y_pos_e)) + { + overlap_list->list[j].is_match = 1; + overlap_list->list[j].strong = 0; + overlap_list->list[j].without_large_indel = 1; + overlap_list->list[j].shared_seed = 1; + } + } + } +} + + +void statistic(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long readNum) +{ + + long long forward, reverse, strong, weak, exact, no_l_indel; + no_l_indel = forward = reverse = exact = strong = weak = 0; + long long i, j; + + for (i = 0; i < readNum; i++) + { + forward += paf[i].length; + reverse += rev_paf[i].length; + for (j = 0; j < paf[i].length; j++) + { + if(paf[i].buffer[j].el == 1) exact++; + if(paf[i].buffer[j].ml == 1) strong++; + if(paf[i].buffer[j].ml == 0) weak++; + if(paf[i].buffer[j].no_l_indel == 1) no_l_indel++; + } + } + + + fprintf(stdout, "****************statistic for overlaps****************\n"); + fprintf(stdout, "overlaps #: %lld\n", forward); + fprintf(stdout, "strong overlaps #: %lld\n", strong); + fprintf(stdout, "weak overlaps #: %lld\n", weak); + fprintf(stdout, "exact overlaps #: %lld\n", exact); + fprintf(stdout, "inexact overlaps #: %lld\n", forward - exact); + fprintf(stdout, "overlaps without large indels#: %lld\n", no_l_indel); + fprintf(stdout, "reverse overlaps #: %lld\n", reverse); + fprintf(stdout, "****************statistic for overlaps****************\n"); +} + + +void fill_chain(Fake_Cigar* chain, char* x_string, char* y_string, long long xBeg, long long yBeg, +long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n) +{ + long long i, xOffset, yOffset, xRegionLen, yRegionLen, /**bandLen,**/ maxXpos, maxYpos, mapScore, zdroped; + float band_rate = 0.08; + int endbouns; + if(chain->length <= 0) return; + + kvec_t(uint8_t) x_num; + kvec_t(uint8_t) y_num; + kv_init(x_num); + kv_init(y_num); + + ///deal with region 0 backward + i = 0; + endbouns = 0; + + xOffset = get_fake_gap_pos(chain, 0); + xOffset = xOffset - 1; + yOffset = (xOffset - xBeg) + yBeg + get_fake_gap_shift(chain, 0); + if(xOffset >= 0 && yOffset >= 0) + { + xRegionLen = xOffset + 1; + yRegionLen = yOffset + 1; + //note here cannot use DIFF(xRegionLen, yRegionLen) + // bandLen = (MIN(xRegionLen, yRegionLen))*band_rate; + // if(bandLen == 0) bandLen = MIN(xRegionLen, yRegionLen); + + ///do alignment backward + kv_resize(uint8_t, x_num, xRegionLen); + kv_resize(uint8_t, y_num, yRegionLen); + ///text is x, query is y + afine_gap_alignment(x_string, x_num.a, xRegionLen, y_string, y_num.a, 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, &mapScore, &zdroped); + fprintf(stderr, "* xOffset: %d, yOffset: %d, xRegionLen: %d, yRegionLen: %d, bandLen: %d, maxXpos: %d, maxYpos: %d, zdroped: %d\n", + xOffset, yOffset, xRegionLen, yRegionLen, BAND_KSW, maxXpos, maxYpos, zdroped); + } + + ///align forward + for (i = 0; i < chain->length; i++) + { + // xOffset = get_fake_gap_pos(chain, i); + // yOffset = xOffset + get_fake_gap_shift(chain, i); + xOffset = get_fake_gap_pos(chain, i); + yOffset = (xOffset - xBeg) + yBeg + get_fake_gap_shift(chain, i); + ///last region + if(i == chain->length - 1) + { + endbouns = 0; + xRegionLen = x_readLen - xOffset; + yRegionLen = y_readLen - yOffset; + //note here cannot use DIFF(xRegionLen, yRegionLen) + // bandLen = (MIN(xRegionLen, yRegionLen))*band_rate; + // if(bandLen == 0) bandLen = MIN(xRegionLen, yRegionLen); + } + else + { + ///higher endbouns for middle regions + endbouns = MATCH_SCORE_KSW; + xRegionLen = get_fake_gap_pos(chain, i+1) - xOffset; + yRegionLen = (get_fake_gap_pos(chain, i+1) + get_fake_gap_shift(chain, i+1)) - + (get_fake_gap_pos(chain, i) + get_fake_gap_shift(chain, i)); + + // bandLen = MAX((MIN(xRegionLen, yRegionLen))*band_rate, DIFF(xRegionLen, yRegionLen)); + // if(bandLen == 0) bandLen = MIN(xRegionLen, yRegionLen); + } + + + ///do alignment forward + kv_resize(uint8_t, x_num, xRegionLen); + kv_resize(uint8_t, y_num, yRegionLen); + ///text is x, query is y + afine_gap_alignment(x_string+xOffset, x_num.a, xRegionLen, y_string+yOffset, y_num.a, 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, &mapScore, &zdroped); + fprintf(stderr, "# xOffset: %d, yOffset: %d, xRegionLen: %d, yRegionLen: %d, bandLen: %d, maxXpos: %d, maxYpos: %d, zdroped: %d\n", + xOffset, yOffset, xRegionLen, yRegionLen, BAND_KSW, maxXpos, maxYpos, zdroped); + } + + + kv_destroy(x_num); + kv_destroy(y_num); +} +void Final_phasing(overlap_region_alloc* overlap_list, Cigar_record_alloc* cigarline, +UC_Read* g_read, UC_Read* overlap_read, uint8_t* c2n) +{ + long long i, xLen, yLen, yStrand; + char* x_string; + char* y_string; + Cigar_record* cigar; + resize_Cigar_record_alloc(cigarline, overlap_list->length); + + + + for (i = 0; i < overlap_list->length; i++) + { + if(overlap_list->list[i].is_match == 1 || + overlap_list->list[i].is_match == 2 || + overlap_list->list[i].is_match == 3) + { + xLen = overlap_list->list[i].x_pos_e - overlap_list->list[i].x_pos_s + 1; + yLen = overlap_list->list[i].y_pos_e - overlap_list->list[i].y_pos_s + 1; + yStrand = overlap_list->list[i].y_pos_strand; + cigar = &(cigarline->buffer[i]); + ///has already been matched exactly + if(overlap_list->list[i].is_match == 1 && overlap_list->list[i].shared_seed == 1) + { + add_cigar_record(g_read->seq + overlap_list->list[i].x_pos_s, xLen, cigar, 0); + } + else + { + if(yStrand == 0) + { + recover_UC_Read(overlap_read, &R_INF, overlap_list->list[i].y_id); + } + else + { + recover_UC_Read_RC(overlap_read, &R_INF, overlap_list->list[i].y_id); + } + x_string = g_read->seq; + y_string = overlap_read->seq; + + fill_chain(&(overlap_list->list[i].f_cigar), x_string, y_string, + overlap_list->list[i].x_pos_s, overlap_list->list[i].y_pos_s, + Get_READ_LENGTH(R_INF, overlap_list->list[i].x_id), + Get_READ_LENGTH(R_INF, overlap_list->list[i].y_id), cigar, c2n); + } + + } + } + +} + void* Final_overlap_calculate_heap_merge(void* arg) { - long long matched_overlap_0 = 0; - long long matched_overlap_1 = 0; - long long potiental_matched_overlap_0 = 0; - long long potiental_matched_overlap_1 = 0; - long long num_correct_base = 0; - long long num_read_base = 0; - long long num_second_correct_base = 0; - long long j, inner_j; - int thr_ID = *((int*)arg); - uint64_t POA_i; long long i = 0; - int avalible_k = 0; UC_Read g_read; init_UC_Read(&g_read); @@ -2915,27 +3190,8 @@ void* Final_overlap_calculate_heap_merge(void* arg) UC_Read overlap_read; init_UC_Read(&overlap_read); - HPC_seq HPC_read; - Hash_code k_code; - uint64_t code; - uint64_t end_pos; - k_mer_pos* list; - uint64_t list_length; - uint64_t sub_ID; - long long total_shared_seed = 0; - long long candidate_overlap_reads = 0; - - - Candidates_list l; - //Candidates_list debug_l; - Graph POA_Graph; - Graph DAGCon; - init_Graph(&DAGCon); - init_Graph(&POA_Graph); - - + Candidates_list l; init_Candidates_list(&l); - //init_Candidates_list(&debug_l); k_mer_pos_list_alloc array_list; init_k_mer_pos_list_alloc(&array_list); @@ -2944,63 +3200,24 @@ void* Final_overlap_calculate_heap_merge(void* arg) init_overlap_region_alloc(&overlap_list); HeapSq heap; - Init_Heap(&heap); - Correct_dumy correct; - init_Correct_dumy(&correct); + Cigar_record_alloc cigarline; + init_Cigar_record_alloc(&cigarline); - - Output_buffer_sub_block current_sub_buffer; - - init_buffer_sub_block(¤t_sub_buffer); - - Cigar_record current_cigar; - init_Cigar_record(¤t_cigar); - - haplotype_evdience_alloc hap; - InitHaplotypeEvdience(&hap); - - - Round2_alignment second_round; - init_Round2_alignment(&second_round); - - small_hash_table forward, reverse; - init_small_hash_table(&forward); - init_small_hash_table(&reverse); - - - long long pre_r_overlaps = 0; - long long cur_r_overlaps = 0; - long long pre_overlaps = 0; - long long cur_overlaps = 0; - if(thr_ID == 0) - { - for (i = 0; i < R_INF.total_reads; i++) - { - pre_r_overlaps += R_INF.reverse_paf[i].length; - pre_overlaps += R_INF.paf[i].length; - } - } - - // if(thr_ID == 0) - // { - // debug_info_of_specfic_read("m64016_190918_162737/130811282/ccs", - // R_INF.paf, R_INF.reverse_paf, -1, "xxxx"); - - // debug_info_of_specfic_read("m64016_190918_162737/179635219/ccs", - // R_INF.paf, R_INF.reverse_paf, -1, "xxxx"); - // } + uint8_t c2n[256]; + memset(c2n, 4, 256); + c2n['A'] = c2n['a'] = 0; c2n['C'] = c2n['c'] = 1; + c2n['G'] = c2n['g'] = 2; c2n['T'] = c2n['t'] = 3; // build the encoding table + + for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num) { - - - ////get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); - get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.001); + get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.001, 0); /** correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon, &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1, @@ -3012,148 +3229,23 @@ void* Final_overlap_calculate_heap_merge(void* arg) overlap_region_sort_y_id(overlap_list.list, overlap_list.length); ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length); ma_hit_sort_tn(R_INF.reverse_paf[i].buffer, R_INF.reverse_paf[i].length); + reverse_complement(g_read.seq, g_read.length); - // if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i), - // Get_NAME_LENGTH((R_INF), i)) == 0) - // { - // fprintf(stderr, "\n1\n"); - // debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "first"); - // } + update_overlaps(&overlap_list, &(R_INF.paf[i]), &g_read, &overlap_read, 1, 1); + update_overlaps(&overlap_list, &(R_INF.reverse_paf[i]), &g_read, &overlap_read, 2, 0); + ///recover missing exact overlaps + update_exact_overlaps(&overlap_list, &g_read, &overlap_read); - // if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i), - // Get_NAME_LENGTH((R_INF), i)) == 0) - // { - // fprintf(stderr, "\n1\n"); - // debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "first"); - // } + ///Final_phasing(&overlap_list, &cigarline, &g_read, &overlap_read, c2n); + push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, + &overlap_list, 1); + push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, + &overlap_list, 2); + - - overlap_list.mapped_overlaps_length = 0; - inner_j = 0; - j = 0; - while (j < overlap_list.length && inner_j < R_INF.paf[i].length) - { - if(overlap_list.list[j].y_id < R_INF.paf[i].buffer[inner_j].tn) - { - j++; - } - else if(overlap_list.list[j].y_id > R_INF.paf[i].buffer[inner_j].tn) - { - inner_j++; - } - else - { - if(overlap_list.list[j].y_pos_strand == R_INF.paf[i].buffer[inner_j].rev) - { - - overlap_list.list[j].is_match = 1; - overlap_list.list[j].strong = R_INF.paf[i].buffer[inner_j].ml; - overlap_list.list[j].without_large_indel = R_INF.paf[i].buffer[inner_j].no_l_indel; - overlap_list.mapped_overlaps_length++; - - if(overlap_list.list[j].strong == 1) - { - matched_overlap_1++; - } - else if(overlap_list.list[j].strong == 0) - { - matched_overlap_0++; - } - else - { - fprintf(stderr, "error\n"); - } - - - if(overlap_list.list[j].without_large_indel == 0) - { - num_second_correct_base++; - } - - - - } - j++; - inner_j++; - } - } - - - - - inner_j = 0; - j = 0; - while (j < overlap_list.length && inner_j < R_INF.reverse_paf[i].length) - { - if(overlap_list.list[j].y_id < R_INF.reverse_paf[i].buffer[inner_j].tn) - { - j++; - } - else if(overlap_list.list[j].y_id > R_INF.reverse_paf[i].buffer[inner_j].tn) - { - inner_j++; - } - else - { - if(overlap_list.list[j].y_pos_strand == R_INF.reverse_paf[i].buffer[inner_j].rev) - { - - overlap_list.list[j].is_match = 2; - overlap_list.list[j].strong = 0; - overlap_list.list[j].without_large_indel = 1; - } - j++; - inner_j++; - } - } - - - ///recover missing exact overlaps - reverse_complement(g_read.seq, g_read.length); - for (j = 0; j < overlap_list.length; j++) - { - if (overlap_list.list[j].is_match != 1) - { - if(overlap_list.list[j].y_pos_strand == 0) - { - recover_UC_Read(&overlap_read, &R_INF, overlap_list.list[j].y_id); - } - else - { - recover_UC_Read_RC(&overlap_read, &R_INF, overlap_list.list[j].y_id); - } - - if(if_exact_match(g_read.seq, g_read.length, overlap_read.seq, overlap_read.length, - overlap_list.list[j].x_pos_s, overlap_list.list[j].x_pos_e, - overlap_list.list[j].y_pos_s, overlap_list.list[j].y_pos_e)) - { - overlap_list.list[j].is_match = 1; - overlap_list.list[j].strong = 0; - overlap_list.list[j].without_large_indel = 1; - overlap_list.mapped_overlaps_length++; - potiental_matched_overlap_0++; - } - - } - } - - - if(R_INF.paf[i].is_fully_corrected) - { - potiental_matched_overlap_1++; - } - - num_correct_base += - push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, - &overlap_list, &g_read, &overlap_read, 1, 1); - - - push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, - &overlap_list, &g_read, &overlap_read, 2, 0); - } @@ -3161,70 +3253,20 @@ void* Final_overlap_calculate_heap_merge(void* arg) finish_output_buffer(); - destory_buffer_sub_block(¤t_sub_buffer); destory_Candidates_list(&l); destory_overlap_region_alloc(&overlap_list); destory_Heap(&heap); - destory_k_mer_pos_list_alloc(&array_list); - destory_Graph(&POA_Graph); - destory_Graph(&DAGCon); + destory_k_mer_pos_list_alloc(&array_list); destory_UC_Read(&g_read); destory_UC_Read(&overlap_read); - destory_Cigar_record(¤t_cigar); - destory_Correct_dumy(&correct); - destoryHaplotypeEvdience(&hap); - destory_Round2_alignment(&second_round); - destory_small_hash_table(&forward); - destory_small_hash_table(&reverse); - - + destory_Cigar_record_alloc(&cigarline); + pthread_mutex_lock(&statistics); - total_matched_overlap_0 += matched_overlap_0; - total_matched_overlap_1 += matched_overlap_1; - total_potiental_matched_overlap_0 += potiental_matched_overlap_0; - total_potiental_matched_overlap_1 += potiental_matched_overlap_1; - total_num_correct_base += num_correct_base; - total_second_num_correct_base += num_second_correct_base; - - complete_threads++; if(complete_threads == thread_num) { - - - fprintf(stderr, "overlaps with large indels: %llu\n", total_second_num_correct_base); - fprintf(stderr, "weak overlaps: %llu\n", total_matched_overlap_0); - fprintf(stderr, "strong overlaps: %llu\n", total_matched_overlap_1); - fprintf(stderr, "recover weak overlaps: %llu\n", total_potiental_matched_overlap_0); - fprintf(stderr, "final available overlaps: %llu\n", total_num_correct_base); - - fprintf(stderr, "fully corrected reads: %llu\n", total_potiental_matched_overlap_1); - - - - for (i = 0; i < R_INF.total_reads; i++) - { - cur_r_overlaps += R_INF.reverse_paf[i].length; - cur_overlaps += R_INF.paf[i].length; - } - fprintf(stderr, "pre_r_overlaps: %d, cur_r_overlaps: %d\n", - pre_r_overlaps, cur_r_overlaps); - - fprintf(stderr, "pre_overlaps: %d, cur_overlaps: %d\n", - pre_overlaps, cur_overlaps); - - - - - - // debug_info_of_specfic_read("m64016_190918_162737/130811282/ccs", - // R_INF.paf, R_INF.reverse_paf, -1, "yyyy"); - - // debug_info_of_specfic_read("m64016_190918_162737/179635219/ccs", - // R_INF.paf, R_INF.reverse_paf, -1, "yyyy"); - - + statistic(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads); } pthread_mutex_unlock(&statistics); free(arg); @@ -3363,24 +3405,9 @@ void Correct_Reads(int last_round) roundID = number_of_round - last_round; fprintf(stdout, "Error correction: start the %d-th round ...\n", roundID); - ///only the first round correction can load index from disk - if (roundID == 0 && load_index_from_disk && load_pre_cauculated_index()) - { - ; - } - else - { - Counting_multiple_thr(); - - Build_hash_table_multiple_thr(); - } - - - fprintf(stdout, "Total pos in hash tabe: %d\n", PCB.total_occ); - fprintf(stdout, "k_mer_min_freq in hashtable: %d\n", k_mer_min_freq); - fprintf(stdout, "k_mer_max_freq in hashtable: %d\n", k_mer_max_freq); - + Counting_multiple_thr(); + Build_hash_table_multiple_thr(); ///verify_Position_hash_table(); Overlap_calculate_multipe_thr(); diff --git a/Correct.cpp b/Correct.cpp index bd2d824..d4bb5bc 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -22,7 +22,7 @@ const uint8_t *c, int sc_mch, int sc_mis, int gapo, int gape, int bandLen, int z int* max_q_pos, int* max_t, int* score) { int i, a = sc_mch, b = sc_mis < 0? sc_mis : -sc_mis; // a>0 and b<0 - int8_t mat[25] = { a,b,b,b,0, b,a,b,b,0, b,b,a,b,0, b,b,b,a,0, 0,0,0,0,0 }; + int8_t mat[25] = {a,b,b,b,0, b,a,b,b,0, b,b,a,b,0, b,b,b,a,0, 0,0,0,0,0}; uint8_t *ts, *qs; ksw_extz_t ez; @@ -43,6 +43,42 @@ int* max_q_pos, int* max_t, int* score) } +void afine_gap_alignment(const char *tseq, uint8_t* tnum, const int tl, +const char *qseq, uint8_t* qnum, const int ql, const uint8_t *c2n, const int strand, +int sc_mch, int sc_mis, int gapo, int gape, int bandLen, int zdrop, int end_bonus, +long long* max_t_pos, long long* max_q_pos, long long* score, long long* droped) +{ + (*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 + int8_t mat[25] = {a,b,b,b,0, b,a,b,b,0, b,b,a,b,0, b,b,b,a,0, 0,0,0,0,0}; + ksw_extz_t ez; + memset(&ez, 0, sizeof(ksw_extz_t)); + + if(strand == FORWARD_KSW) + { + for (i = 0; i < tl; ++i) tnum[i] = c2n[(uint8_t)tseq[i]]; // encode to 0/1/2/3 + for (i = 0; i < ql; ++i) qnum[i] = c2n[(uint8_t)qseq[i]]; + } + else if(strand == BACKWARD_KSW) + { + for (i = 0; i < tl; ++i) tnum[i] = c2n[(uint8_t)tseq[tl - i - 1]]; // encode to 0/1/2/3 + for (i = 0; i < ql; ++i) qnum[i] = c2n[(uint8_t)qseq[ql - i - 1]]; + } + + ksw_extz2_sse(0, ql, qnum, tl, tnum, 5, mat, gapo, gape, bandLen, zdrop, end_bonus, 0, &ez); + + (*score) = ez.max; + (*max_t_pos) = ez.max_t; + (*max_q_pos) = ez.max_q; + (*droped) = ez.zdropped; + + /** + for (i = 0; i < ez.n_cigar; ++i) // print CIGAR + printf("%d%c", ez.cigar[i]>>4, "MID"[ez.cigar[i]&0xf]); + putchar('\n'); + **/ + free(ez.cigar); +} void clear_Round2_alignment(Round2_alignment* h) @@ -1574,6 +1610,239 @@ int verify_cigar(char* x, int x_len, char* y, int y_len, CIGAR* cigar, int error } + +int scan_cigar(CIGAR* cigar, int* get_error, int scanXLen, int direction) +{ + (*get_error) = -1; + if(cigar->length == 1 && cigar->C_C[0] == 0) + { + (*get_error) = 0; + return 1; + } + + + int x_i, y_i, cigar_i; + x_i = 0; + y_i = 0; + + int operation; + int operationLen; + int i; + int cigar_error = 0; + + ///0 is match, 1 is mismatch, 2 is up, 3 is left + ///2: there are more bases at y, 3: there are more bases at x + if(direction == 0) + { + cigar_i = 0; + while (cigar_i < cigar->length) + { + operation = cigar->C_C[cigar_i]; + operationLen = cigar->C_L[cigar_i]; + cigar_i++; + + if (operation == 0) + { + x_i = x_i + operationLen; + y_i = y_i + operationLen; + if(x_i >= scanXLen) + { + (*get_error) = cigar_error; + return 1; + } + } + else if (operation == 1) + { + for (i = 0; i < operationLen; i++) + { + x_i++; + y_i++; + cigar_error++; + if(x_i >= scanXLen) + { + (*get_error) = cigar_error; + return 1; + } + } + }///2是x缺字符（y多字符） + else if (operation == 2) + { + cigar_error += operationLen; + y_i += operationLen; + }///3是y缺字符（x多字符） + else if (operation == 3) + { + for (i = 0; i < operationLen; i++) + { + x_i++; + cigar_error++; + if(x_i >= scanXLen) + { + (*get_error) = cigar_error; + return 1; + } + } + } + } + } + else + { + cigar_i = cigar->length - 1; + while (cigar_i >= 0) + { + operation = cigar->C_C[cigar_i]; + operationLen = cigar->C_L[cigar_i]; + cigar_i--; + if (operation == 0) + { + x_i = x_i + operationLen; + y_i = y_i + operationLen; + if(x_i >= scanXLen) + { + (*get_error) = cigar_error; + return 1; + } + } + else if (operation == 1) + { + for (i = 0; i < operationLen; i++) + { + x_i++; + y_i++; + cigar_error++; + if(x_i >= scanXLen) + { + (*get_error) = cigar_error; + return 1; + } + } + }///2是x缺字符（y多字符） + else if (operation == 2) + { + cigar_error += operationLen; + y_i += operationLen; + }///3是y缺字符（x多字符） + else if (operation == 3) + { + for (i = 0; i < operationLen; i++) + { + x_i++; + cigar_error++; + if(x_i >= scanXLen) + { + (*get_error) = cigar_error; + return 1; + } + } + } + } + } + + (*get_error) = cigar_error; + return 0; +} + +///[scanXbeg, scanXend] +int scan_cigar_interval(CIGAR* cigar, int* get_error, int scanXbeg, int scanXend) +{ + (*get_error) = -1; + if(cigar->length == 1 && cigar->C_C[0] == 0) + { + (*get_error) = 0; + return 1; + } + + + + int x_i, y_i, cigar_i; + x_i = 0; + y_i = 0; + + int operation; + int operationLen; + int i; + int cigar_error = 0; + + ///0 is match, 1 is mismatch, 2 is up, 3 is left + ///2: there are more bases at y, 3: there are more bases at x + cigar_i = 0; + while (cigar_i < cigar->length) + { + operation = cigar->C_C[cigar_i]; + operationLen = cigar->C_L[cigar_i]; + cigar_i++; + + if (operation == 0) + { + for (i = 0; i < operationLen; i++) + { + if(x_i == scanXbeg) + { + cigar_error = 0; + } + + x_i++; + y_i++; + + if(x_i == scanXend + 1) + { + (*get_error) = cigar_error; + return 1; + } + } + } + else if (operation == 1) + { + for (i = 0; i < operationLen; i++) + { + if(x_i == scanXbeg) + { + cigar_error = 0; + } + + x_i++; + y_i++; + cigar_error++; + + if(x_i == scanXend + 1) + { + (*get_error) = cigar_error; + return 1; + } + } + }///2是x缺字符（y多字符） + else if (operation == 2) + { + cigar_error += operationLen; + y_i += operationLen; + }///3是y缺字符（x多字符） + else if (operation == 3) + { + for (i = 0; i < operationLen; i++) + { + if(x_i == scanXbeg) + { + cigar_error = 0; + } + + + x_i++; + cigar_error++; + + if(x_i == scanXend + 1) + { + (*get_error) = cigar_error; + return 1; + } + } + } + } + + + (*get_error) = cigar_error; + return 0; +} + inline int move_gap_greedy(char* path, int path_i, int path_length, char* x, int x_i, char* y, int y_i, unsigned int* new_error) { if(path[path_i] < 2) @@ -3513,32 +3782,6 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN ///note!!! need notification end_site = Reserve_Banded_BPM(y_string, Window_Len, x_string, x_len, threshold, &error); - - - // unsigned int debug_error; - // int return_y_end, return_x_end; - // Reserve_Banded_BPM_Extension(y_string, Window_Len, x_string, x_len, threshold, - // &debug_error, &return_y_end, &return_x_end); - // if(return_x_end != x_len - 1 && end_site != -1) - // { - // fprintf(stderr, "error\n"); - // } - - // if(return_x_end == x_len - 1 && end_site == -1) - // { - // fprintf(stderr, "error\n"); - // } - - // if(return_x_end == x_len - 1 && end_site != -1) - // { - // if(return_y_end != end_site || debug_error != error) - // { - // fprintf(stderr, "error\n"); - // } - // } - - - ///error等于-1说明没匹配 if (error!=(unsigned int)-1) { @@ -3885,6 +4128,788 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN } } + +void debug_scan_cigar(overlap_region* sub_list) +{ + long long i; + int f_err, b_err, fLen, bLen, xLen; + for (i = 0; i < sub_list->w_list_length; i++) + { + if(sub_list->w_list[i].y_end == -1 || sub_list->w_list[i].cigar.length == -1) + { + continue; + } + xLen = sub_list->w_list[i].x_end - sub_list->w_list[i].x_start + 1; + + scan_cigar(&(sub_list->w_list[i].cigar), &b_err, + xLen, 1); + scan_cigar(&(sub_list->w_list[i].cigar), &f_err, + xLen, 0); + + if(b_err != sub_list->w_list[i].error || f_err != sub_list->w_list[i].error) + { + fprintf(stderr, "sbsbsbsbsb1\n"); + } + + scan_cigar(&(sub_list->w_list[i].cigar), &b_err, + WINDOW, 1); + scan_cigar(&(sub_list->w_list[i].cigar), &f_err, + WINDOW, 0); + + if(b_err != sub_list->w_list[i].error || f_err != sub_list->w_list[i].error) + { + fprintf(stderr, "sbsbsbsbsb2\n"); + } + + scan_cigar_interval(&(sub_list->w_list[i].cigar), &b_err, 0, xLen-1); + if(b_err != sub_list->w_list[i].error) + { + fprintf(stderr, "sbsbsbsbsb3\n"); + } + + fLen = xLen / 3; + scan_cigar(&(sub_list->w_list[i].cigar), &f_err, fLen, 0); + scan_cigar_interval(&(sub_list->w_list[i].cigar), &b_err, 0, fLen-1); + if(f_err != b_err) + { + fprintf(stderr, "sbsbsbsbsb4\n"); + } + + fLen = xLen / 3; + scan_cigar(&(sub_list->w_list[i].cigar), &f_err, fLen, 1); + scan_cigar_interval(&(sub_list->w_list[i].cigar), &b_err, xLen-fLen, xLen-1); + if(f_err != b_err) + { + fprintf(stderr, "\nsbsbsbsbsb5\n"); + fprintf(stderr, "b_err: %d, f_err: %d\n",b_err, f_err); + long long j; + for (j = 0; j < sub_list->w_list[i].cigar.length; j++) + { + fprintf(stderr, "len: %d, opera: %d\n", + sub_list->w_list[i].cigar.C_L[j], sub_list->w_list[i].cigar.C_C[j]); + } + } + + // bLen = xLen / 3; + // fLen = xLen - bLen; + + // scan_cigar(&(sub_list->w_list[i].cigar), &b_err, + // bLen, 1); + // scan_cigar(&(sub_list->w_list[i].cigar), &f_err, + // fLen, 0); + + // if(b_err + f_err != sub_list->w_list[i].error) + // { + // fprintf(stderr, "\nsub_list->w_list[i].error: %d, bLen: %d, b_err: %d, fLen: %d, f_err: %d\n", + // sub_list->w_list[i].error, bLen, b_err, fLen, f_err); + // long long j; + // for (j = 0; j < sub_list->w_list[i].cigar.length; j++) + // { + // fprintf(stderr, "len: %d, opera: %d\n", + // sub_list->w_list[i].cigar.C_L[j], sub_list->w_list[i].cigar.C_C[j]); + // } + + // } + } +} + +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); + int x_id = sub_list->x_id; + int y_id = sub_list->y_id; + int y_strand = sub_list->y_pos_strand; + long long x_readLen = Get_READ_LENGTH((*R_INF), x_id); + long long y_readLen = Get_READ_LENGTH((*R_INF), y_id); + long long i, y_distance; + int f_err, b_err, m_error; + int scanLen = 10; + long long boundaryLen = 200; + long long single_sideLen = boundaryLen/2; + long long useless_side = single_sideLen/2; + int alpha = 1; + long long y_start, x_start; + long long y_end, 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; + for (i = 0; i + 1 < sub_list->w_list_length; i++) + { + ///if two windows are not aligned + 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; + } + + 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_cigar(&(sub_list->w_list[i].cigar), &b_err, scanLen, 1); + 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; + } + 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; + } + else + { + sub_list->boundary_cigars.buffer[i].error = -1; + sub_list->boundary_cigars.buffer[i].y_end = -1; + continue; + } + + + 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 + leftLen = MIN(MIN((x_start - sub_list->w_list[i].x_start), y_start), single_sideLen); + rightLen = MIN(MIN((sub_list->w_list[i+1].x_end + 1 - x_start), y_readLen - y_start), + single_sideLen); + + + xLen = leftLen + rightLen; + x_start = x_start - leftLen; + x_end = x_start + xLen - 1; + y_start = y_start - leftLen; + + ///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; + } + + + threshold = xLen * THRESHOLD_RATE; + threshold = Adjust_Threshold(threshold, xLen); + threshold = double_error_threshold(threshold, xLen); + + yLen = xLen + (threshold << 1); + if(!determine_overlap_region(threshold, y_start, y_id, yLen, R_INF, + &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(dumy->overlap_region, y_start, o_len, y_strand, + R_INF, 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); + + 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); + + + ///y_distance can be less than 0, or larger than 0 + 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; + } + + + scan_cigar_interval(&(sub_list->boundary_cigars.buffer[i].cigar), &m_error, + useless_side, xLen-useless_side-1); + scan_cigar(&(sub_list->w_list[i].cigar), &b_err, leftLen-useless_side, 1); + scan_cigar(&(sub_list->w_list[i+1].cigar), &f_err, rightLen-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; + + + } + else + { + sub_list->boundary_cigars.buffer[i].error = -1; + sub_list->boundary_cigars.buffer[i].y_end = -1; + continue; + } + + } + +} + +inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_reads* R_INF, + UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read) +{ + long long j, k, i; + long long Len_x; + int threshold; + long long y_len; + long long currentIDLen; + long long matches; + 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 y_end; + 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; + + + ///j负责遍历整个overlap list + for (j = 0; j < overlap_list->length; j++) + { + + y_id = overlap_list->list[j].y_id; + y_strand = overlap_list->list[j].y_pos_strand; + y_readLen = Get_READ_LENGTH((*R_INF), y_id); + + //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 < 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_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, R_INF, + &extra_begin, &extra_end, &y_start, &o_len)) + { + break; + } + + if(o_len + threshold < x_len) + { + break; + } + + fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, + R_INF, 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); + + ///error等于-1说明没匹配 + 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 < 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(dumy->overlap_region, y_start, o_len, y_strand, + R_INF, 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_boundary(x_string, x_len, threshold, y_start, real_y_start, + end_site, extra_begin, extra_end, y_id, Window_Len, R_INF, 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); + + // if(real_y_start < extra_begin || end_site >= Window_Len - extra_end) + // { + // fprintf(stderr, "\nreal_y_start: %d, extra_begin: %d\n", + // real_y_start, extra_begin); + + // fprintf(stderr, "end_site: %d, Window_Len: %d, extra_end: %d\n", + // end_site, Window_Len, extra_end); + // } + + + ///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; + ///k遍历匹配window左侧所有不匹配的window + ///如果i匹配，则k从i-1开始 + ///直到第一个匹配的window结束 + ///因为i!=0，所以k的大小不用担心 + 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_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, R_INF, + &extra_begin, &extra_end, &y_start, &o_len)) + { + break; + } + + if(o_len + threshold < x_len) + { + break; + } + + fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, + R_INF, 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_boundary(x_string, x_len, threshold, y_start, real_y_start, end_site, + extra_begin, extra_end, y_id, Window_Len, R_INF, 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); + + // if(real_y_start < extra_begin || end_site >= Window_Len - extra_end) + // { + // fprintf(stderr, "\nreal_y_start: %d, extra_begin: %d\n", + // real_y_start, extra_begin); + + // fprintf(stderr, "end_site: %d, Window_Len: %d, extra_end: %d\n", + // end_site, Window_Len, extra_end); + // } + + ///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->mapped_overlaps_length = 0; + + int pre_threshold; + long long tLen, tError; + double error_rate; + for (j = 0; j < overlap_list->length; j++) + { + y_id = overlap_list->list[j].y_id; + y_strand = overlap_list->list[j].y_pos_strand; + y_readLen = Get_READ_LENGTH((*R_INF), y_id); + 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 * OVERLAP_THRESHOLD_FILTER <= overlap_list->list[j].align_length) + { + + for (i = 0; i < 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) + { + ///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 * THRESHOLD_RATE; + 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(dumy->overlap_region, y_start, o_len, y_strand, + R_INF, 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 != overlap_list->list[j].w_list[i].error) + // { + // fprintf(stderr, "error\n"); + // } + + if (error != (unsigned int)-1) + { + if (end_site == Window_Len - 1 || real_y_start == 0) + { + + if(fix_boundary(x_string, x_len, threshold, y_start, real_y_start, end_site, + extra_begin, extra_end, y_id, Window_Len, R_INF, 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); + + + // if(real_y_start < extra_begin || end_site >= Window_Len - extra_end) + // { + // fprintf(stderr, "\nreal_y_start: %d, extra_begin: %d\n", + // real_y_start, extra_begin); + + // fprintf(stderr, "end_site: %d, Window_Len: %d, extra_end: %d\n", + // end_site, Window_Len, extra_end); + // } + + ///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 = trim_error_rate(overlap_list, j); + error_rate = non_trim_error_rate(overlap_list, j, R_INF, dumy, g_read); + + + ///if(error_rate <= 0.015) + if(error_rate <= 0.03) + { + overlap_list->mapped_overlaps_length += overlap_length; + overlap_list->list[j].is_match = 1; + calculate_boundary_cigars(&(overlap_list->list[j]), R_INF, dumy, g_read); + } + else if(error_rate <= 0.045) + { + overlap_list->list[j].is_match = 3; + } + } + } + + + /** + for (j = 0; j < overlap_list->length; j++) + { + y_id = overlap_list->list[j].y_id; + y_strand = overlap_list->list[j].y_pos_strand; + if(overlap_list->list[j].is_match == 1) + { + for (i = 0; i < overlap_list->list[j].w_list_length; i++) + { + if(overlap_list->list[j].w_list[i].y_end != -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; + + x_string = g_read->seq + x_start; + + y_start = overlap_list->list[j].w_list[i].y_start; + y_end = overlap_list->list[j].w_list[i].y_end; + y_len = y_end - y_start + 1; + + recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, R_INF, y_id); + y_string = dumy->overlap_region; + + + if(verify_cigar(x_string, x_len, y_string, y_len, &overlap_list->list[j].w_list[i].cigar, + overlap_list->list[j].w_list[i].error)) + { + fprintf(stderr, "j: %d, i: %d, y_id: %d, y_start: %d, y_end: %d\n", j, i, y_id, y_start, y_end); + } + } + + } + + + + for (i = 0; i < overlap_list->list[j].boundary_cigars.length; i++) + { + if(overlap_list->list[j].boundary_cigars.buffer[i].y_end != -1) + { + x_start = overlap_list->list[j].boundary_cigars.buffer[i].x_start; + x_end = overlap_list->list[j].boundary_cigars.buffer[i].x_end; + x_len = x_end - x_start + 1; + x_string = g_read->seq + x_start; + + y_start = overlap_list->list[j].boundary_cigars.buffer[i].y_start; + y_end = overlap_list->list[j].boundary_cigars.buffer[i].y_end; + y_len = y_end - y_start + 1; + + recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, + R_INF, y_id); + y_string = dumy->overlap_region; + + + if(verify_cigar(x_string, x_len, y_string, y_len, + &overlap_list->list[j].boundary_cigars.buffer[i].cigar, + overlap_list->list[j].boundary_cigars.buffer[i].error)) + { + fprintf(stderr, "j: %d, i: %d, y_id: %d, y_start: %d, y_end: %d\n", j, i, y_id, y_start, y_end); + } + } + } + } + } + **/ + +} + + + + + /** inline void adjust_alignment_windows(overlap_region* overlap, All_reads* R_INF) { @@ -7671,6 +8696,126 @@ CIGAR* cigar, haplotype_evdience_alloc* hap) } +///[xBeg, xEnd] +void markSNP_detail(CIGAR* cigar_record, uint8_t* flag, +long long xBeg, long long xEnd, long long flag_offset) +{ + 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; + + while (cigar_i < cigar_record->length) + { + operation = cigar_record->C_C[cigar_i]; + operationLen = cigar_record->C_L[cigar_i]; + if(x_i > xEnd) + { + break; + } + + + ///match + if (operation == 0) + { + x_i += operationLen; + y_i += operationLen; + } + else if(operation == 1) ///mismatch + { + for (i = 0; i < operationLen; i++) + { + ///flag_offset + if(flag[x_i - flag_offset] < 127 && x_i >= xBeg && x_i <= xEnd) + { + flag[x_i - flag_offset]++; + } + + x_i++; + y_i++; + } + }///insertion + else if (operation == 2) + { + y_i += operationLen; + } + else if (operation == 3) + { + x_i += operationLen; + } + + cigar_i++; + } +} + +void markSNP_advance( +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) +{ + long long x_total_end = x_total_start + x_length - 1; + ///mismatches based on the offset of x + long long inner_offset = x_total_start - window_offset; + CIGAR* cigar_record; + long long useless_side, xleftLen, xrightLen, x_interval_beg, x_interval_end; + long long current_cigar_beg, current_cigar_end; + + ///for current_cigar, [current_cigar_beg, current_cigar_end] + current_cigar_beg = 0; + current_cigar_end = x_length - 1; + + if(beg_cigar!=NULL && beg_cigar->y_end!=-1) + { + useless_side = beg_cigar->error_threshold; + cigar_record = &(beg_cigar->cigar); + + xleftLen = x_total_start - beg_cigar->x_start; + xrightLen = beg_cigar->x_end - x_total_start + 1; + ///actually xleftLen could be no larger than useless_side + ///but such window has already been filtered out in calculate_boundary_cigars + if(xleftLen > useless_side && xrightLen > useless_side) + { + x_interval_beg = xleftLen; + x_interval_end = x_interval_beg + (xrightLen - useless_side) - 1; + ///this is the offset of the current cigar + current_cigar_beg = xrightLen - useless_side; + markSNP_detail(cigar_record, hap->flag + inner_offset, x_interval_beg, + x_interval_end, x_interval_beg); + } + } + + if(end_cigar!=NULL && end_cigar->y_end!=-1) + { + useless_side = end_cigar->error_threshold; + cigar_record = &(end_cigar->cigar); + xleftLen = x_total_end - end_cigar->x_start; + xrightLen = end_cigar->x_end - x_total_end + 1; + + ///actually xrightLen could be no larger than useless_side + ///but such window has already been filtered out in calculate_boundary_cigars + if(xleftLen > useless_side && xrightLen > useless_side) + { + x_interval_end = xleftLen; + x_interval_beg = x_interval_end - (xleftLen + 1 - useless_side) + 1; + current_cigar_end = (x_length - 1) - (xleftLen + 1 - useless_side); + markSNP_detail(cigar_record, hap->flag + end_cigar->x_start - window_offset, + x_interval_beg, x_interval_end, 0); + } + } + + + + + // markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, 0, x_length/2); + // markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, x_length/2+1, x_length - 1); + ///markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, 0, x_length - 1, 0); + markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, current_cigar_beg, + current_cigar_end, 0); +} + + void addSNPtohaplotype( @@ -7781,6 +8926,8 @@ CIGAR* cigar, haplotype_evdience_alloc* hap, int snp_threshold) + + void cluster(char* r_string, long long window_start, long long window_end, overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplotype_evdience_alloc* hap) { @@ -7907,17 +9054,26 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplot } -void cluster_advance(char* r_string, long long window_start, long long window_end, -overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplotype_evdience_alloc* hap) +void get_related_cigars(window_list_alloc* boundary_cigars, long long id, window_list** beg_cigar, +window_list** end_cigar) { + (*beg_cigar) = &(boundary_cigars->buffer[id*2]); + (*end_cigar) = &(boundary_cigars->buffer[id*2+1]); +} + +void cluster_advance(char* r_string, long long window_start, long long window_end, +overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, +haplotype_evdience_alloc* hap) +{ + 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; - long long x_length; + long long x_start, x_end, x_length; char* x_string; char* y_string; long long i; - long long y_start, y_length; + long long y_start, y_end, y_length; long long overlapID, windowID; long long startNodeID, endNodeID, currentNodeID; @@ -7926,6 +9082,8 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplot int snp_threshold; snp_threshold = 1; + ///resize_window_list_alloc(boundary_cigars, dumy->length * 2); + ///all overlaps related to the current window [window_start, window_end] ///first mark all snp pos for (i = 0; i < dumy->length; i++) @@ -7949,14 +9107,30 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplot 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; + x_end = overlap_list->list[overlapID].w_list[windowID].x_end; + 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; + y_end = overlap_list->list[overlapID].w_list[windowID].y_end; + + beg_cigar = end_cigar = NULL; + + if(windowID >= 1) + { + beg_cigar = &(overlap_list->list[overlapID].boundary_cigars.buffer[windowID-1]); + } + + if(windowID < overlap_list->list[overlapID].w_list_length - 1) + { + end_cigar = &(overlap_list->list[overlapID].boundary_cigars.buffer[windowID]); + } - markSNP(window_start, x_start, x_length, y_start, y_length, - &(overlap_list->list[overlapID].w_list[windowID].cigar), hap); + + markSNP_advance(window_start, x_start, x_length, y_start, y_length, + &(overlap_list->list[overlapID].w_list[windowID]), beg_cigar, end_cigar, hap); } @@ -8024,6 +9198,17 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplot y_string = dumy->overlap_region; + beg_cigar = end_cigar = NULL; + if(windowID >= 1) + { + beg_cigar = &(overlap_list->list[overlapID].boundary_cigars.buffer[windowID-1]); + } + if(windowID < overlap_list->list[overlapID].w_list_length - 1) + { + end_cigar = &(overlap_list->list[overlapID].boundary_cigars.buffer[windowID]); + } + + addSNPtohaplotype(window_start, overlapID, x_string, x_start, x_length, y_string, y_start, y_length, &(overlap_list->list[overlapID].w_list[windowID].cigar), hap, snp_threshold); @@ -11760,12 +12945,11 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, } recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read); + ///recalcate_window_advance(overlap_list, R_INF, g_read, dumy, overlap_read); - - ////partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat); - partition_overlaps_advance(overlap_list, R_INF, g_read, dumy, hap, force_repeat); - + partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat); + ///partition_overlaps_advance(overlap_list, R_INF, g_read, dumy, hap, force_repeat); // print_overlap("m64016_190918_162737/174131552/ccs", // overlap_list->list[0].x_id, overlap_list, R_INF, 1); @@ -12225,3 +13409,49 @@ uint64_t n_end_pos, uint8_t n_direction, UC_Read* g_read, All_reads* R_INF, Corr /**********************for prefilter************************ */ + + +void init_Cigar_record_alloc(Cigar_record_alloc* x) +{ + x->length = 0; + x->size = 0; + x->buffer = NULL; +} + +void resize_Cigar_record_alloc(Cigar_record_alloc* x, long long new_size) +{ + long long i; + if(new_size > x->size) + { + x->buffer = (Cigar_record*)realloc(x->buffer, new_size*sizeof(Cigar_record)); + for (i = 0; i < x->size; i++) + { + clear_Cigar_record(&(x->buffer[i])); + } + for (; i < new_size; i++) + { + init_Cigar_record(&(x->buffer[i])); + clear_Cigar_record(&(x->buffer[i])); + } + + x->size = new_size; + } + else + { + for (i = 0; i < new_size; i++) + { + clear_Cigar_record(&(x->buffer[i])); + } + } + + x->length = 0; +} +void destory_Cigar_record_alloc(Cigar_record_alloc* x) +{ + long long i; + for (i = 0; i < x->size; i++) + { + destory_Cigar_record(&(x->buffer[i])); + } + free(x->buffer); +} diff --git a/Correct.h b/Correct.h index 13b0630..95b068c 100644 --- a/Correct.h +++ b/Correct.h @@ -16,10 +16,13 @@ #define INSERTION 2 #define DELETION 3 +#define MIN(x,y) ((x)<=(y)?(x):(y)) + ///#define FLAG_THRE 0 #define MAX(x, y) ((x >= y)?x:y) #define MIN(x, y) ((x <= y)?x:y) +#define DIFF(x, y) ((MAX((x), (y))) - (MIN((x), (y)))) #define OVERLAP(x_start, x_end, y_start, y_end) (MIN(x_end, y_end) - MAX(x_start, y_start) + 1) ///#define OVERLAP(x_start, x_end, y_start, y_end) MIN(x_end, y_end) - MAX(x_start, y_start) + 1 @@ -131,6 +134,14 @@ typedef struct }Cigar_record; +typedef struct +{ + long long length; + long long size; + Cigar_record* buffer; +}Cigar_record_alloc; + + typedef struct { ////the position of snp in read itself @@ -1436,4 +1447,22 @@ void append_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list, k_mer_pos uint64_t n_end_pos, uint8_t n_direction, UC_Read* g_read, All_reads* R_INF, Correct_dumy* dumy); /**********************for prefilter************************ */ + +void init_Cigar_record_alloc(Cigar_record_alloc* x); +void resize_Cigar_record_alloc(Cigar_record_alloc* x, long long new_size); +void destory_Cigar_record_alloc(Cigar_record_alloc* x); + +void afine_gap_alignment(const char *tseq, uint8_t* tnum, const int tl, +const char *qseq, uint8_t* qnum, const int ql, const uint8_t *c2n, const int strand, +int sc_mch, int sc_mis, int gapo, int gape, int bandLen, int zdrop, int end_bonus, +long long* max_t_pos, long long* max_q_pos, long long* score, long long* droped); + +#define FORWARD_KSW 0 +#define BACKWARD_KSW 1 +#define MATCH_SCORE_KSW 2 +#define MISMATCH_SCORE_KSW 4 +#define GAP_OPEN_KSW 4 +#define GAP_EXT_KSW 2 +#define Z_DROP_KSW 400 +#define BAND_KSW 50 #endif \ No newline at end of file diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 1b40777..8026afa 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -291,6 +291,7 @@ void init_overlap_region_alloc(overlap_region_alloc* list) for (i = 0; i < list->size; i++) { init_fake_cigar(&(list->list[i].f_cigar)); + init_window_list_alloc(&(list->list[i].boundary_cigars)); } } void clear_overlap_region_alloc(overlap_region_alloc* list) @@ -302,6 +303,7 @@ void clear_overlap_region_alloc(overlap_region_alloc* list) { list->list[i].w_list_length = 0; clear_fake_cigar(&(list->list[i].f_cigar)); + clear_window_list_alloc(&(list->list[i].boundary_cigars)); } } @@ -315,6 +317,7 @@ void destory_overlap_region_alloc(overlap_region_alloc* list) free(list->list[i].w_list); } destory_fake_cigar(&(list->list[i].f_cigar)); + destory_window_list_alloc(&(list->list[i].boundary_cigars)); } free(list->list); } @@ -851,7 +854,8 @@ int append_inexact_overlap_region_alloc_back(overlap_region_alloc* list, overlap } -int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, All_reads* R_INF) +int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, +All_reads* R_INF, int add_beg_end) { if (list->length + 1 > list->size) @@ -923,9 +927,16 @@ 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)); - add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0); + if(add_beg_end == 1) + { + add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0); + } + long long distance_gap; - long long pre_distance_gap = 0; + /****************************may have bugs********************************/ + ///long long pre_distance_gap = 0; + long long pre_distance_gap = 0xfffffffffffffff; + /****************************may have bugs********************************/ long long i = 0; for (i = 0; i < tmp->f_cigar.length; i++) { @@ -939,8 +950,8 @@ int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_regi } } - if(get_fake_gap_pos(&(list->list[list->length].f_cigar), - list->list[list->length].f_cigar.length - 1) != list->list[list->length].x_pos_e) + if(add_beg_end == 1 && get_fake_gap_pos(&(list->list[list->length].f_cigar), + list->list[list->length].f_cigar.length - 1) != list->list[list->length].x_pos_e) { add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_e, @@ -996,11 +1007,18 @@ 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)); - add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0); + if(add_beg_end == 1) + { + add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0); + } + long long distance_self_pos = tmp->x_pos_e - tmp->x_pos_s; long long distance_pos = tmp->y_pos_e - tmp->y_pos_s; long long init_distance_gap = distance_pos - distance_self_pos; - long long pre_distance_gap = init_distance_gap; + /****************************may have bugs********************************/ + ///long long pre_distance_gap = init_distance_gap; + long long pre_distance_gap = 0xfffffffffffffff; + /****************************may have bugs********************************/ long long distance_gap; long long i = 0; for (i = tmp->f_cigar.length - 1; i >= 0; i--) @@ -1015,7 +1033,7 @@ int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_regi } } - if(get_fake_gap_pos(&(list->list[list->length].f_cigar), + if(add_beg_end == 1 && get_fake_gap_pos(&(list->list[list->length].f_cigar), list->list[list->length].f_cigar.length - 1) != list->list[list->length].x_pos_e) { add_fake_cigar(&(list->list[list->length].f_cigar), @@ -1059,6 +1077,19 @@ int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_regi /******************************for debug********************************/ } + // if(list->list[list->length].f_cigar.length < 3 && tmp->f_cigar.length != 1) + // { + // fprintf(stderr, "\n original cigar:\n"); + // print_fake_gap(&tmp->f_cigar); + // fprintf(stderr, "new cigar:\n"); + // print_fake_gap(&list->list[list->length].f_cigar); + // fprintf(stderr, "xs: %d, xe: %d, strand: %d, xLen: %d\n", + // list->list[list->length].x_pos_s, + // list->list[list->length].x_pos_e, + // tmp->x_pos_strand, + // Get_READ_LENGTH((*R_INF), tmp->x_id)); + // } + list->list[list->length].shared_seed = tmp->shared_seed; list->list[list->length].align_length = 0; @@ -1778,7 +1809,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) clear_fake_cigar(&(result->f_cigar)); - + ///not a has been sorted by offset, that means has been sorted by query offset i = max_i; result->x_pos_e = a[i].self_offset; result->y_pos_e = a[i].offset; @@ -1790,6 +1821,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) long long pre_distance_gap = distance_pos - distance_self_pos; ///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); long long chainLen = 0; if(result->x_pos_strand == 1) @@ -1938,7 +1970,7 @@ All_reads* R_INF) void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, -uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold) +uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end) { overlap_region tmp_region; uint64_t i = 0; @@ -1967,12 +1999,13 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh current_ID = candidates->list[i].readID; current_stand = candidates->list[i].strand; - ///这个是查询read的信息 + ///reference read tmp_region.x_id = readID; tmp_region.x_pos_strand = current_stand; - ///这个是被查询的read的信息 + ///query read tmp_region.y_id = current_ID; - tmp_region.y_pos_strand = 0; ///永远是0 + ///here the strand of query is always 0 + tmp_region.y_pos_strand = 0; @@ -2007,7 +2040,7 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh ///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1) if (tmp_region.x_id != tmp_region.y_id) { - append_inexact_overlap_region_alloc(overlap_list, &tmp_region, R_INF); + append_inexact_overlap_region_alloc(overlap_list, &tmp_region, R_INF, add_beg_end); ///append_inexact_overlap_region_alloc_back(overlap_list, &tmp_region, R_INF); } } @@ -2158,7 +2191,7 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF) ///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1) if (tmp_region.x_id != tmp_region.y_id) { - append_inexact_overlap_region_alloc(overlap_list, &tmp_region, R_INF); + append_inexact_overlap_region_alloc(overlap_list, &tmp_region, R_INF, 1); } @@ -2937,8 +2970,310 @@ int load_Total_Pos_Table(Total_Pos_Table* TCB, char* read_file_name) return 1; } +typedef struct +{ + long long* list; + uint64_t length; +} H_peaks; -void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k_mer_min_freq, int k_mer_max_freq) +void insert_H_peaks(H_peaks* h, long long index, long long value) +{ + if(h->length <= index) + { + long long newLen = index + 1; + h->list = (long long*)realloc(h->list, newLen*sizeof(long long)); + memset(h->list + h->length, 0, sizeof(long long) * (newLen - h->length)); + h->length = newLen; + } + + h->list[index] += value; +} + +inline void RC_Hash_code(Hash_code* code, Hash_code* rc_code, int k) +{ + rc_code->x[0] = 0; + rc_code->x[1] = 0; + int i; + for (i = 0; i < k; i++) + { + rc_code->x[0] = rc_code->x[0] << 1; + rc_code->x[1] = rc_code->x[1] << 1; + rc_code->x[0] |= (((uint64_t)((code->x[0] >> i) & 1))^((uint64_t)1)); + rc_code->x[1] |= (((uint64_t)((code->x[1] >> i) & 1))^((uint64_t)1)); + } +} + + + +void get_peak_debug(Total_Count_Table* TCB, long long* min, long long* max) +{ + int i; + Count_Table* h; + khint_t k; + long long c_count; + H_peaks LH; + LH.list = NULL; + LH.length = 0; + uint64_t sub_ID; + uint64_t sub_key; + Hash_code code, rc_code, debug_code; + char str[100]; + char rc_str[100]; + + + + for (i = 0; i < TCB->size; i++) + { + h = TCB->sub_h[i]; + for (k = kh_begin(h); k != kh_end(h); ++k) + { + if (kh_exist(h, k)) // test if a bucket contains data + { + sub_ID = i; + sub_key = kh_key(h, k); + + recover_hash_code(sub_ID, sub_key, &code, TCB->suffix_mode, + TCB->suffix_bits, k_mer_length); + RC_Hash_code(&code, &rc_code, k_mer_length); + RC_Hash_code(&rc_code, &debug_code, k_mer_length); + if(code.x[0] != debug_code.x[0] || code.x[1] != debug_code.x[1]) + { + fprintf(stderr, "sbsbsb\n"); + } + + Hashcode_to_string(&code, str, k_mer_length); + Hashcode_to_string(&rc_code, rc_str, k_mer_length); + reverse_complement(str, k_mer_length); + if(memcmp(str, rc_str, k_mer_length) != 0) + { + fprintf(stderr, "hehehehe\n"); + int j; + for (j = 0; j < k_mer_length; j++) + { + fprintf(stderr, "%c",str[j]); + } + fprintf(stderr, "\n"); + + for (j = 0; j < k_mer_length; j++) + { + fprintf(stderr, "%c",rc_str[j]); + } + fprintf(stderr, "\n"); + + } + + + ///get_Total_Count_Table(&TCB, &k_code, k_mer_length); + + c_count = kh_value(h, k); + + if(get_Total_Count_Table(TCB, &code, k_mer_length) != c_count) + { + fprintf(stderr, "sbsbsb\n"); + } + + insert_H_peaks(&LH, c_count, c_count); + } + } + } + + (*max) = -1; + (*min) = -1; + long long max_value = -1; + for (i = 0; i < LH.length; i++) + { + if(LH.list[i] >= max_value) + { + max_value = LH.list[i]; + (*max) = i; + } + } + + long long min_value = max_value; + for (i = 0; i < LH.length; i++) + { + if(LH.list[i] < min_value && LH.list[i] != 0) + { + min_value = LH.list[i]; + (*min) = i; + } + } + + for (i = 0; i < LH.length; i++) + { + ///fprintf(stderr, "%d, %d\n", i, LH.list[i]); + fprintf(stderr, "%d\n", LH.list[i]); + } + + + + free(LH.list); +} + +///1: a > b; -1: a < b; 0: a=b +int cmp_Hash_code(Hash_code* a, Hash_code* b) +{ + if(a->x[1] > b->x[1]) + { + return 1; + } + if(a->x[1] < b->x[1]) + { + return -1; + } + ///a->x[1] == b->x[1] + if(a->x[0] > b->x[0]) + { + return 1; + } + if(a->x[0] < b->x[0]) + { + return -1; + } + + return 0; +} + +int get_total_freq(Total_Count_Table* TCB, uint64_t sub_ID, uint64_t sub_key, long long* T_count) +{ + Hash_code code, rc_code; + long long count, rc_count; + + recover_hash_code(sub_ID, sub_key, &code, TCB->suffix_mode, + TCB->suffix_bits, k_mer_length); + RC_Hash_code(&code, &rc_code, k_mer_length); + + count = get_Total_Count_Table(TCB, &code, k_mer_length); + rc_count = get_Total_Count_Table(TCB, &rc_code, k_mer_length); + (*T_count) = count + rc_count; + + if(count == 0) + { + return 0; + }///count > 0 && rc_count == 0 + else if(rc_count == 0) + { + return 1; + }///count > 0 && rc_count > 0 + else + { + int flag = cmp_Hash_code(&code, &rc_code); + + ///code > rc_code + if(flag > 0) + { + return 1; + }///code < rc_code + else if(flag < 0) + { + return 0; + } + else + { + (*T_count) = (*T_count)/2; + return 1; + } + } +} + +void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long* up_boundary) +{ + int i; + Count_Table* h; + khint_t k; + long long count; + H_peaks LH; + LH.list = NULL; + LH.length = 0; + uint64_t sub_ID; + uint64_t sub_key; + + + + for (i = 0; i < TCB->size; i++) + { + h = TCB->sub_h[i]; + for (k = kh_begin(h); k != kh_end(h); ++k) + { + if (kh_exist(h, k)) // test if a bucket contains data + { + sub_ID = i; + sub_key = kh_key(h, k); + + if(get_total_freq(TCB, sub_ID, sub_key, &count)==1) + { + insert_H_peaks(&LH, count, count); + } + } + } + } + + (*max) = -1; + (*min) = -1; + long long max_value = -1; + //// seed with freq 1 is useless + for (i = 2; i < LH.length; i++) + { + if(LH.list[i] >= max_value) + { + max_value = LH.list[i]; + (*max) = i; + } + } + + long long opt = 4; + (*up_boundary) = -1; + for (i = (*max) + opt; i < LH.length; i++) + { + if(LH.list[i] > LH.list[i-opt]) + { + long long j = i-opt; + for (; j < i; j++) + { + if(LH.list[j] < LH.list[j+1]) + { + (*up_boundary) = j; + goto end_opt; + } + } + + (*up_boundary) = i; + goto end_opt; + } + } + + end_opt: + if((*up_boundary) == -1 || (*up_boundary) > (*max) * 10) + { + (*up_boundary) = (*max) * 10; + } + + + + + long long min_value = max_value; + //// seed with freq 1 is useless + for (i = 2; i < LH.length && i < (*max); i++) + { + if(LH.list[i] < min_value && LH.list[i] != 0) + { + min_value = LH.list[i]; + (*min) = i; + } + } + + // for (i = 0; i < LH.length; i++) + // { + // ///fprintf(stderr, "%d, %d\n", i, LH.list[i]); + // fprintf(stderr, "%d\n", LH.list[i]); + // } + // fflush(stderr); + + + free(LH.list); +} + +void Traverse_Counting_Table_back(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k_mer_min_freq, int k_mer_max_freq) { int i; Count_Table* h; @@ -2948,6 +3283,12 @@ void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k PCB->useful_k_mer = 0; PCB->total_occ = 0; + long long freq_min, freq_max, freq_up; + ///get_peak_debug(TCB, &freq_min, &freq_max); + get_peak(TCB, &freq_min, &freq_max, &freq_up); + fprintf(stderr, "freq_min: %d, freq_max: %d, freq_up: %d\n", + freq_min, freq_max, freq_up); + ///init_Total_Pos_Table(PCB, TCB); khint_t t; ///这就是个迭代器 @@ -3025,6 +3366,115 @@ void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k } +void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k_mer_min_freq, int k_mer_max_freq) +{ + int i; + Count_Table* h; + khint_t k; + uint64_t sub_key; + uint64_t sub_ID; + PCB->useful_k_mer = 0; + PCB->total_occ = 0; + + long long freq_min, max, freq_up; + ///get_peak_debug(TCB, &freq_min, &freq_max); + get_peak(TCB, &freq_min, &max, &freq_up); + fprintf(stdout, "freq_min: %d, freq_max: %d, freq_up:%d\n", + freq_min, max, freq_up); + if(freq_min < k_mer_min_freq) + { + k_mer_min_freq = freq_min; + } + if(freq_up > k_mer_max_freq) + { + k_mer_max_freq = freq_up; + } + + fprintf(stdout, "k_mer_min_freq: %d, k_mer_max_freq: %d\n", + k_mer_min_freq, k_mer_max_freq); + + + khint_t t; ///这就是个迭代器 + int absent; + long long count; + + /******************************************** + hash_table(key) ----> PCB->k_mer_index ------> PCB->pos + ********************************************/ + for (i = 0; i < TCB->size; i++) + { + h = TCB->sub_h[i]; + for (k = kh_begin(h); k != kh_end(h); ++k) + { + if (kh_exist(h, k)) // test if a bucket contains data + { + sub_ID = i; + sub_key = kh_key(h, k); + get_total_freq(TCB, sub_ID, sub_key, &count); + + ///只有符合频率范围要求的k-mer，才会被加入到pos table中 + if (count>=k_mer_min_freq && count<=k_mer_max_freq) + { + t = kh_put(POS64, PCB->sub_h[sub_ID], sub_key, &absent); + + if (absent) + { + ///kh_value(PCB->sub_h[sub_ID], t) = useful_k_mer + total_occ; + kh_value(PCB->sub_h[sub_ID], t) = PCB->useful_k_mer; + } + else ///哈希表中已有的元素 + { + ///kh_value(PCB->sub_h[sub_ID], t)++; + fprintf(stderr, "ERROR\n"); + } + + + PCB->useful_k_mer++; + PCB->total_occ = PCB->total_occ + kh_value(h, k); + } + } + } + } + + + fprintf(stdout, "useful_k_mer: %lld\n",PCB->useful_k_mer); + fprintf(stdout, "total_occ: %lld\n",PCB->total_occ); + + PCB->k_mer_index = (uint64_t*)malloc(sizeof(uint64_t)*(PCB->useful_k_mer+1)); + + PCB->k_mer_index[0] = 0; + + PCB->total_occ = 0; + PCB->useful_k_mer = 0; + + for (i = 0; i < TCB->size; i++) + { + h = TCB->sub_h[i]; + for (k = kh_begin(h); k != kh_end(h); ++k) + { + if (kh_exist(h, k)) // test if a bucket contains data + { + sub_ID = i; + sub_key = kh_key(h, k); + get_total_freq(TCB, sub_ID, sub_key, &count); + + ///if (kh_value(h, k)>=k_mer_min_freq && kh_value(h, k)<=k_mer_max_freq) + if (count>=k_mer_min_freq && count<=k_mer_max_freq) + { + PCB->useful_k_mer++; + PCB->total_occ = PCB->total_occ + kh_value(h, k); + PCB->k_mer_index[PCB->useful_k_mer] = PCB->total_occ; + } + } + } + } + + PCB->pos = (k_mer_pos*)malloc(sizeof(k_mer_pos)*PCB->total_occ); + memset(PCB->pos, 0, sizeof(k_mer_pos)*PCB->total_occ); + + +} + @@ -3802,5 +4252,44 @@ void resize_fake_cigar(Fake_Cigar* x, long long size) x->buffer = (uint64_t*)realloc(x->buffer, sizeof(uint64_t) * x->size); } + x->length = 0; +} + + +void init_window_list_alloc(window_list_alloc* x) +{ + x->buffer = NULL; + x->length = 0; + x->size = 0; +} + +void clear_window_list_alloc(window_list_alloc* x) +{ + x->length = 0; +} + +void destory_window_list_alloc(window_list_alloc* x) +{ + if(x->size != 0) + { + free((x->buffer)); + } +} + +void resize_window_list_alloc(window_list_alloc* x, long long size) +{ + if(size > x->size) + { + x->size = size; + x->buffer = (window_list*)realloc(x->buffer, sizeof(window_list) * x->size); + } + + long long i; + for (i = 0; i < x->size; i++) + { + x->buffer[i].error = -1; + } + + x->length = 0; } \ No newline at end of file diff --git a/Hash_Table.h b/Hash_Table.h index 8b48caf..175084b 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -100,6 +100,14 @@ typedef struct } window_list; +typedef struct +{ + window_list* buffer; + long long length; + long long size; +}window_list_alloc; + + typedef struct { uint64_t* buffer; @@ -133,6 +141,8 @@ typedef struct uint64_t w_list_length; int8_t strong; Fake_Cigar f_cigar; + + window_list_alloc boundary_cigars; } overlap_region; @@ -250,7 +260,25 @@ inline int if_k_mer_available(Hash_code* code, int k) return 1; } +////suffix_bits = 64 in default +inline int recover_hash_code(uint64_t sub_ID, uint64_t sub_key, Hash_code* code, +uint64_t suffix_mode, int suffix_bits, int k) +{ + uint64_t h_key, low_key; + h_key = low_key = 0; + low_key = sub_ID << SAFE_SHIFT(suffix_bits); + low_key = low_key | sub_key; + + h_key = sub_ID >> (64 - suffix_bits); + + code->x[0] = code->x[1] = 0; + uint64_t mask = ALL >> (64 - k); + code->x[0] = low_key & mask; + + code->x[1] = h_key << (64 - k); + code->x[1] = code->x[1] | (low_key >> SAFE_SHIFT(k)); +} ///inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, Total_Count_Table* TCB, Hash_code* code, int k) inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t suffix_mode, int suffix_bits, @@ -279,7 +307,13 @@ Hash_code* code, int k) *get_sub_ID = sub_ID; *get_sub_key = sub_key; - return 1; + + // Hash_code de_code; + // recover_hash_code(sub_ID, sub_key, &de_code, suffix_mode, suffix_bits, k); + ///if(de_code.x[0] != (*code).x[0] || de_code.x[1] != (*code).x[1]) fprintf(stderr, "hehe\n"); + ///if(de_code.x[0] == (*code).x[0] || de_code.x[1] == (*code).x[1]) fprintf(stderr, "hehe\n"); + + return 1; } @@ -514,7 +548,7 @@ void calculate_inexact_overlap_region(Candidates_list* candidates, overlap_regio uint64_t readID, uint64_t readLength, All_reads* R_INF); void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, -uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold); +uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end); @@ -703,4 +737,12 @@ void append_overlap_region_alloc_from_existing(overlap_region_alloc* list, overl int cmp_by_x_pos_s(const void * a, const void * b); void resize_Chain_Data(Chain_Data* x, long long size); + + +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); + + #endif \ No newline at end of file diff --git a/Overlaps.cpp b/Overlaps.cpp index 0abeeff..7458830 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -9890,7 +9890,7 @@ char* output_file_name, long long bubble_dist, int read_graph, int write) // debug_info_of_specfic_node("m64016_190918_162737/141297762/ccs", sg); out: - output_tips(sg, &R_INF); + ///output_tips(sg, &R_INF); output_unitig_graph(sg, coverage_cut, output_file_name, n_read); diff --git a/kmer.h b/kmer.h index b49be00..b3f735d 100644 --- a/kmer.h +++ b/kmer.h @@ -4,8 +4,9 @@ ///#define ALL (0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff) #define ALL (0xffffffffffffffff) - +/****************************may have bugs********************************/ #define SAFE_SHIFT(k) k & ((k < 64)?ALL:0) +/****************************may have bugs********************************/ @@ -94,6 +95,21 @@ inline void k_mer_append(Hash_code* code, uint64_t c, int k) code->x[1] = ((code->x[1]<<1) | (c>>1)) & mask; } +inline void Hashcode_to_string(Hash_code* code, char* str, int k) +{ + uint8_t c; + int i; + for (i = 0; i < k; i++) + { + c = (code->x[1] >> (k - i - 1)) & ((uint64_t)1); + c = c << 1; + c = c | ((code->x[0] >> (k - i - 1)) & ((uint64_t)1)); + + str[i] = s_H[c]; + } + +} + void init_HPC_seq(HPC_seq* seq, char* str, long long l); void init_Hash_code(Hash_code* code);