diff --git a/Assembly.cpp b/Assembly.cpp index 7d47d6e..6d828b3 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -13,7 +13,10 @@ Total_Count_Table TCB; Total_Pos_Table PCB; All_reads R_INF; - +pthread_mutex_t statistics; +long long total_matched_overlap_0 = 0; +long long total_matched_overlap_1 = 0; +long long complete_threads = 0; void* Perform_Counting(void* arg) @@ -679,6 +682,14 @@ void* Overlap_calculate(void* arg) void* Overlap_calculate_heap_merge(void* arg) { + /************需要注释掉**********/ + // long long debug_overlap = 0; + // long long filtered_debug_overlap = 0; + /************需要注释掉**********/ + + long long matched_overlap_0 = 0; + long long matched_overlap_1 = 0; + long long j; int thr_ID = *((int*)arg); uint64_t POA_i; @@ -721,6 +732,7 @@ void* Overlap_calculate_heap_merge(void* arg) for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num) { + clear_Heap(&heap); clear_Candidates_list(&l); @@ -749,10 +761,14 @@ void* Overlap_calculate_heap_merge(void* arg) if (list_length != 0) { - append_k_mer_pos_list_alloc(&array_list, list, list_length, end_pos, 0); + /************需要注释掉**********/ + // debug_overlap = debug_overlap + list_length; + /************需要注释掉**********/ + ///append_k_mer_pos_list_alloc(&array_list, list, list_length, end_pos, 0); + append_k_mer_pos_list_alloc_prefilter(&array_list, list, list_length, end_pos, 0, &g_read, &R_INF, &correct); } ///merge_Candidates_list(&l, list, list_length, end_pos, 0); - //merge_Candidates_list_version(&debug_l, list, list_length, end_pos, 0); + ///merge_Candidates_list_version(&debug_l, list, list_length, end_pos, 0); } } else @@ -805,7 +821,11 @@ void* Overlap_calculate_heap_merge(void* arg) list_length = locate_Total_Pos_Table(&PCB, &k_code, &list, k_mer_length, &sub_ID); if (list_length != 0) { - append_k_mer_pos_list_alloc(&array_list, list, list_length, end_pos, 1); + /************需要注释掉**********/ + // debug_overlap = debug_overlap + list_length; + /************需要注释掉**********/ + ///append_k_mer_pos_list_alloc(&array_list, list, list_length, end_pos, 1); + append_k_mer_pos_list_alloc_prefilter(&array_list, list, list_length, end_pos, 1, &g_read, &R_INF, &correct); } ///merge_Candidates_list(&l, list, list_length, end_pos, 1); //merge_Candidates_list_version(&debug_l, list, list_length, end_pos, 1); @@ -820,6 +840,16 @@ void* Overlap_calculate_heap_merge(void* arg) ///HPC_base++; } + + + + /************需要注释掉**********/ + // for (int ijk = 0; ijk < array_list.length; ijk++) + // { + // filtered_debug_overlap = filtered_debug_overlap + array_list.list[ijk].length; + // } + /************需要注释掉**********/ + ///merge_k_mer_pos_list_alloc(&array_list, &l); merge_k_mer_pos_list_alloc_heap_sort(&array_list, &l, &heap); /** @@ -840,6 +870,22 @@ void* Overlap_calculate_heap_merge(void* arg) correct_overlap(&overlap_list, &R_INF, &g_read, &correct); + for (j = 0; j < overlap_list.length; j++) + { + long long Len_x = overlap_list.list[j].x_pos_e - overlap_list.list[j].x_pos_s + 1; + + if (Len_x * 0.6 <= overlap_list.list[j].align_length) + { + if (overlap_list.list[j].y_pos_strand == 0) + { + matched_overlap_0++; + } + else + { + matched_overlap_1++; + } + } + } /** POA_i = 0; @@ -918,6 +964,16 @@ void* Overlap_calculate_heap_merge(void* arg) } + + /************需要注释掉**********/ + // fprintf(stderr, "debug_overlap: %llu\n", debug_overlap); + // fprintf(stderr, "filtered_debug_overlap: %llu\n", filtered_debug_overlap); + /************需要注释掉**********/ + + + + + /** fprintf(stderr, "candidate_overlap_reads: %llu\n", candidate_overlap_reads); fprintf(stderr, "total_shared_seed: %llu\n", total_shared_seed); @@ -927,12 +983,25 @@ void* Overlap_calculate_heap_merge(void* arg) //destory_Candidates_list(&debug_l); destory_Heap(&heap); - destory_k_mer_pos_list_alloc(&array_list); + ///destory_k_mer_pos_list_alloc(&array_list); + destory_k_mer_pos_list_alloc_prefilter(&array_list); destory_Graph(&POA_Graph); destory_UC_Read(&g_read); destory_Correct_dumy(&correct); + + + pthread_mutex_lock(&statistics); + total_matched_overlap_0 += matched_overlap_0; + total_matched_overlap_1 += matched_overlap_1; + complete_threads++; + if(complete_threads == thread_num) + { + fprintf(stderr, "total_matched_overlap_0: %llu\n", total_matched_overlap_0); + fprintf(stderr, "total_matched_overlap_1: %llu\n", total_matched_overlap_1); + } + pthread_mutex_unlock(&statistics); } diff --git a/Correct.cpp b/Correct.cpp index 6f8565d..6f820a8 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -12,10 +12,7 @@ long long T_total_mis=0; pthread_mutex_t debug_statistics ; -#define MAX(x, y) ((x >= y)?x:y) -#define MIN(x, y) ((x <= y)?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 + @@ -392,6 +389,7 @@ char* r_string) int groupLen = 0; int return_sites[GROUP_SIZE]; unsigned int return_sites_error[GROUP_SIZE]; + uint64_t overlapID[GROUP_SIZE]; ///这些是整个window被完全覆盖的 for (i = 0; i < dumy->length; i++) @@ -402,6 +400,13 @@ char* r_string) x_start = window_start; ///y上的相对位置 y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s; + + // /************需要注释掉**********/ + // if(x_start < overlap_list->list[currentID].x_pos_s) + // { + // fprintf(stderr, "ERROR\n"); + // } + // /************需要注释掉**********/ ///y上的起始 y_start = y_start - THRESHOLD; if (y_start < 0) @@ -412,6 +417,13 @@ char* r_string) currentIDLen = Get_READ_LENGTH((*R_INF), overlap_list->list[currentID].y_id); ///不能超过y的剩余长度 o_len = MIN(Window_Len, currentIDLen - y_start); + + // /************需要注释掉**********/ + // if(o_len < x_len) + // { + // fprintf(stderr, "ERROR\n"); + // } + // /************需要注释掉**********/ ///这个长度足够,可以用来一起比 @@ -420,6 +432,8 @@ char* r_string) recover_UC_Read_sub_region(dumy->overlap_region_group[groupLen], y_start, o_len, overlap_list->list[currentID].y_pos_strand, R_INF, overlap_list->list[currentID].y_id); + + overlapID[groupLen] = currentID; x_string = r_string + x_start; @@ -432,30 +446,29 @@ char* r_string) return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE); groupLen = 0; + if (return_sites_error[0]!=(unsigned int)-1) + { + overlap_list->list[overlapID[0]].align_length += x_len; + } + + if (return_sites_error[1]!=(unsigned int)-1) + { + overlap_list->list[overlapID[1]].align_length += x_len; + } + + if (return_sites_error[2]!=(unsigned int)-1) + { + overlap_list->list[overlapID[2]].align_length += x_len; + } + + if (return_sites_error[3]!=(unsigned int)-1) + { + overlap_list->list[overlapID[3]].align_length += x_len; + } - // end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], o_len, x_string, x_len, THRESHOLD, &error); - // if (error!=return_sites_error[0]) - // { - // fprintf(stderr, "error\n"); - // } - // end_site = Reserve_Banded_BPM(dumy->overlap_region_group[1], o_len, x_string, x_len, THRESHOLD, &error); - // if (error!=return_sites_error[1]) - // { - // fprintf(stderr, "error\n"); - // } - // end_site = Reserve_Banded_BPM(dumy->overlap_region_group[2], o_len, x_string, x_len, THRESHOLD, &error); - // if (error!=return_sites_error[2]) - // { - // fprintf(stderr, "error\n"); - // } - // end_site = Reserve_Banded_BPM(dumy->overlap_region_group[3], o_len, x_string, x_len, THRESHOLD, &error); - // if (error!=return_sites_error[3]) - // { - // fprintf(stderr, "error\n"); - // } - // /************需要注释掉**********/ + /************需要注释掉**********/ // for (size_t ijk = 0; ijk < GROUP_SIZE; ijk++) // { // if (return_sites_error[ijk]!=(unsigned int)-1) @@ -469,7 +482,7 @@ char* r_string) // test_edit_distance_by_edlib(x_string, dumy->overlap_region_group[ijk], x_len, o_len, THRESHOLD, // return_sites_error[ijk], &total_mis); // } - // /************需要注释掉********* */ + /************需要注释掉********* */ } @@ -486,8 +499,14 @@ char* r_string) ///o_len < Window_Len, 说明y的长度不够,需要在y后面补N memset (y_string + o_len, 0, Window_Len - o_len); end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, THRESHOLD, &error); + + + if (error!=(unsigned int)-1) + { + overlap_list->list[currentID].align_length += x_len; + } - // /************需要注释掉**********/ + /************需要注释掉**********/ // if (error!=(unsigned int)-1) // { // total_match++; @@ -497,7 +516,7 @@ char* r_string) // total_unmatch++; // } // test_edit_distance_by_edlib(x_string, y_string, x_len, o_len, THRESHOLD, error, &total_mis); - // /************需要注释掉********* */ + /************需要注释掉********* */ } } @@ -506,7 +525,13 @@ char* r_string) if (groupLen == 1) { end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], o_len, x_string, x_len, THRESHOLD, &error); - // /************需要注释掉**********/ + + if (error!=(unsigned int)-1) + { + overlap_list->list[overlapID[0]].align_length += x_len; + } + + /************需要注释掉**********/ // if (error!=(unsigned int)-1) // { // total_match++; @@ -516,14 +541,22 @@ char* r_string) // total_unmatch++; // } // test_edit_distance_by_edlib(x_string, dumy->overlap_region_group[0], x_len, o_len, THRESHOLD, error, &total_mis); - // /************需要注释掉********* */ + /************需要注释掉********* */ } 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], o_len, x_string, x_len, return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE); - // /************需要注释掉**********/ + + for (i = 0; i < groupLen; i++) + { + if (return_sites_error[i]!=(unsigned int)-1) + { + overlap_list->list[overlapID[i]].align_length += x_len; + } + } + /************需要注释掉**********/ // for (size_t ijk = 0; ijk < groupLen; ijk++) // { // if (return_sites_error[ijk]!=(unsigned int)-1) @@ -537,7 +570,7 @@ char* r_string) // test_edit_distance_by_edlib(x_string, dumy->overlap_region_group[ijk], x_len, o_len, THRESHOLD, // return_sites_error[ijk], &total_mis); // } - // /************需要注释掉**********/ + /************需要注释掉**********/ groupLen = 0; } @@ -554,13 +587,23 @@ char* r_string) ///这个是和当前窗口重叠的长度 x_len = x_end - x_start + 1; - if (x_len <= 0) - { - fprintf(stderr, "ERROR\n"); - } + // /************需要注释掉**********/ + // if (x_len <= 0) + // { + // fprintf(stderr, "ERROR\n"); + // } + // /************需要注释掉**********/ ///y上的相对位置 y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s; + + // /************需要注释掉**********/ + // if(x_start < overlap_list->list[currentID].x_pos_s) + // { + // fprintf(stderr, "ERROR\n"); + // } + // /************需要注释掉**********/ + ///y上的起始 y_start = y_start - THRESHOLD; if (y_start < 0) @@ -575,6 +618,15 @@ char* r_string) ///不能超过y的剩余长度 Window_Len = x_len + (THRESHOLD << 1); o_len = MIN(Window_Len, currentIDLen - y_start); + + + // /************需要注释掉**********/ + // if(o_len < x_len) + // { + // fprintf(stderr, "ERROR\n"); + // } + // /************需要注释掉**********/ + recover_UC_Read_sub_region(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, R_INF, overlap_list->list[currentID].y_id); @@ -590,8 +642,13 @@ char* r_string) } end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, THRESHOLD, &error); + if (error!=(unsigned int)-1) + { + overlap_list->list[currentID].align_length += x_len; + } - // /************需要注释掉********* */ + + /************需要注释掉********* */ // if (error!=(unsigned int)-1) // { // total_match++; @@ -601,19 +658,19 @@ char* r_string) // total_unmatch++; // } // test_edit_distance_by_edlib(x_string, y_string, x_len, o_len, THRESHOLD, error, &total_mis); - // /************需要注释掉********* */ + /************需要注释掉********* */ } - // /************需要注释掉********* */ + /************需要注释掉********* */ // pthread_mutex_lock(&debug_statistics); // T_total_match = T_total_match + total_match; // T_total_unmatch = T_total_unmatch + total_unmatch; // T_total_mis = T_total_mis + total_mis; // pthread_mutex_unlock(&debug_statistics); - // /************需要注释掉********* */ + /************需要注释掉********* */ } @@ -670,12 +727,12 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Re - // /************需要注释掉********* */ + /************需要注释掉********* */ // pthread_mutex_lock(&debug_statistics); // fprintf(stderr, "total_match: %u, total_unmatch: %u, total_mis: %u\n", // T_total_match, T_total_unmatch, T_total_mis); // pthread_mutex_unlock(&debug_statistics); - // /************需要注释掉********* */ + /************需要注释掉********* */ } @@ -710,4 +767,357 @@ void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list) list->overlapID = (uint64_t*)realloc(list->overlapID, list->size*sizeof(uint64_t)); } -} \ No newline at end of file +} + + +void pre_filter_by_nearby_single(k_mer_pos* new_n_list, k_mer_pos* old_n_list, uint64_t n_length, uint64_t n_end_pos, UC_Read* g_read, +All_reads* R_INF, Correct_dumy* dumy, uint64_t* new_n_length) +{ + (*new_n_length) = 0; + ///这种就是0,new_n_list也不需要有数据 + if (n_length == 0) + { + return; + } + + char* x_string = NULL; + char* y_string = NULL; + + + long long x_offset = (long long)(n_end_pos) - WINDOW + 1; + if (x_offset < 0) + { + x_offset = 0; + } + + ///x_length是x上待验证区间的实际长度 + ///x_offset是x上待验证区间的起始位置 + ///如果是向前取待验证区间,那么x_length至少大于等于40 (compressed k-mer长度,不是实际长度) + long long x_length = n_end_pos - x_offset + 1; + ///这种也不要过滤了,直接把old_n_list全部赋值过来就好了 + if (x_length < (WINDOW/2)) + { + (*new_n_length) = n_length; + memcpy(new_n_list, old_n_list, sizeof(k_mer_pos)*n_length); + return; + } + long long Window_Len = x_length + (THRESHOLD << 1); + x_string = g_read->seq + x_offset; + + + long long y_offset; + long long y_length; + + long long i = 0; + + long long y_read_length; + int end_site; + unsigned int error; + + for (i = 0; i < n_length; i++) + { + ////old_n_list[i].offset是y上k-mer的结束位置 + ///n_end_pos是x上k-mer的结束位置 + ///x_length是x上区间长度 + y_offset = (long long)(old_n_list[i].offset) - x_length; + ///这种情况下弃疗 + if (y_offset < 0) + { + new_n_list[(*new_n_length)].readID = old_n_list[i].readID; + new_n_list[(*new_n_length)].offset = old_n_list[i].offset; + (*new_n_length)++; + continue; + } + + y_offset = y_offset - THRESHOLD; + ///还能抢救 + if (y_offset < 0) + { + y_offset = 0; + } + + ///y的read的总长度 + y_read_length = Get_READ_LENGTH((*R_INF), old_n_list[i].readID); + y_length = MIN(Window_Len, y_read_length - y_offset); + + ///如果y的长度比x长度还小,那就直接弃疗了 + if (y_length < x_length) + { + new_n_list[(*new_n_length)].readID = old_n_list[i].readID; + new_n_list[(*new_n_length)].offset = old_n_list[i].offset; + (*new_n_length)++; + continue; + } + + ///y的方向都是0,因为索引里都是0 + recover_UC_Read_sub_region(dumy->overlap_region, y_offset, y_length, 0, R_INF, old_n_list[i].readID); + y_string = dumy->overlap_region; + + memset (y_string + y_length, 0, Window_Len - y_length); + end_site = Reserve_Banded_BPM(y_string, y_length, x_string, x_length, THRESHOLD, &error); + + if (error!=(unsigned int)-1) + { + new_n_list[(*new_n_length)].readID = old_n_list[i].readID; + new_n_list[(*new_n_length)].offset = old_n_list[i].offset; + (*new_n_length)++; + } + } +} + + + +void pre_filter_by_nearby(k_mer_pos* new_n_list, k_mer_pos* old_n_list, uint64_t n_length, uint64_t n_end_pos, UC_Read* g_read, +All_reads* R_INF, Correct_dumy* dumy, uint64_t* new_n_length) +{ + (*new_n_length) = 0; + ///这种就是0,new_n_list也不需要有数据 + if (n_length == 0) + { + return; + } + + char* x_string = NULL; + char* y_string = NULL; + + + long long x_offset = (long long)(n_end_pos) - WINDOW + 1; + if (x_offset < 0) + { + x_offset = 0; + } + + ///x_length是x上待验证区间的实际长度 + ///x_offset是x上待验证区间的起始位置 + ///如果是向前取待验证区间,那么x_length至少大于等于40 (compressed k-mer长度,不是实际长度) + long long x_length = n_end_pos - x_offset + 1; + ///这种也不要过滤了,直接把old_n_list全部赋值过来就好了 + if (x_length < (WINDOW/2)) + { + (*new_n_length) = n_length; + memcpy(new_n_list, old_n_list, sizeof(k_mer_pos)*n_length); + return; + } + long long Window_Len = x_length + (THRESHOLD << 1); + x_string = g_read->seq + x_offset; + + + long long y_offset; + long long y_length; + + long long i = 0; + + long long y_read_length; + int end_site; + unsigned int error; + + int groupLen = 0; + int return_sites[GROUP_SIZE]; + unsigned int return_sites_error[GROUP_SIZE]; + uint64_t readID[GROUP_SIZE]; + uint64_t offset[GROUP_SIZE]; + + for (i = 0; i < n_length; i++) + { + ////old_n_list[i].offset是y上k-mer的结束位置 + ///n_end_pos是x上k-mer的结束位置 + ///x_length是x上区间长度 + y_offset = (long long)(old_n_list[i].offset) - x_length; + ///这种情况下弃疗 + if (y_offset < 0) + { + new_n_list[(*new_n_length)].readID = old_n_list[i].readID; + new_n_list[(*new_n_length)].offset = old_n_list[i].offset; + (*new_n_length)++; + continue; + } + + y_offset = y_offset - THRESHOLD; + ///还能抢救 + if (y_offset < 0) + { + y_offset = 0; + } + + ///y的read的总长度 + y_read_length = Get_READ_LENGTH((*R_INF), old_n_list[i].readID); + y_length = MIN(Window_Len, y_read_length - y_offset); + + ///如果y的长度比x长度还小,那就直接弃疗了 + if (y_length < x_length) + { + new_n_list[(*new_n_length)].readID = old_n_list[i].readID; + new_n_list[(*new_n_length)].offset = old_n_list[i].offset; + (*new_n_length)++; + continue; + } + + if(y_length == Window_Len) + { + ///y的方向都是0,因为索引里都是0 + recover_UC_Read_sub_region(dumy->overlap_region_group[groupLen], y_offset, y_length, 0, + R_INF, old_n_list[i].readID); + readID[groupLen] = old_n_list[i].readID; + offset[groupLen] = old_n_list[i].offset; + + + 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], y_length, x_string, x_length, + return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE); + groupLen = 0; + + if (return_sites_error[0]!=(unsigned int)-1) + { + new_n_list[(*new_n_length)].readID = readID[0]; + new_n_list[(*new_n_length)].offset = offset[0]; + (*new_n_length)++; + } + + if (return_sites_error[1]!=(unsigned int)-1) + { + new_n_list[(*new_n_length)].readID = readID[1]; + new_n_list[(*new_n_length)].offset = offset[1]; + (*new_n_length)++; + } + + if (return_sites_error[2]!=(unsigned int)-1) + { + new_n_list[(*new_n_length)].readID = readID[2]; + new_n_list[(*new_n_length)].offset = offset[2]; + (*new_n_length)++; + } + + if (return_sites_error[3]!=(unsigned int)-1) + { + new_n_list[(*new_n_length)].readID = readID[3]; + new_n_list[(*new_n_length)].offset = offset[3]; + (*new_n_length)++; + } + + } + + } + else + { + ///y的方向都是0,因为索引里都是0 + recover_UC_Read_sub_region(dumy->overlap_region, y_offset, y_length, 0, R_INF, old_n_list[i].readID); + y_string = dumy->overlap_region; + + memset (y_string + y_length, 0, Window_Len - y_length); + end_site = Reserve_Banded_BPM(y_string, y_length, x_string, x_length, THRESHOLD, &error); + + if (error!=(unsigned int)-1) + { + new_n_list[(*new_n_length)].readID = old_n_list[i].readID; + new_n_list[(*new_n_length)].offset = old_n_list[i].offset; + (*new_n_length)++; + } + } + + + + } + + + + if (groupLen == 1) + { + end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], Window_Len, x_string, x_length, THRESHOLD, &error); + + if (error!=(unsigned int)-1) + { + new_n_list[(*new_n_length)].readID = readID[0]; + new_n_list[(*new_n_length)].offset = offset[0]; + (*new_n_length)++; + } + + } + else + { + 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, x_length, + return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE); + + + for (i = 0; i < groupLen; i++) + { + if (return_sites_error[i]!=(unsigned int)-1) + { + new_n_list[(*new_n_length)].readID = readID[i]; + new_n_list[(*new_n_length)].offset = offset[i]; + (*new_n_length)++; + } + } + + groupLen = 0; + + + } + + ///要排序.... + if ((*new_n_length)>1) + { + qsort(new_n_list, (*new_n_length), sizeof(k_mer_pos), cmp_k_mer_pos); + } +} + + + +/**********************for prefilter************************ */ + +void destory_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list) +{ + long long i = 0; + for (i = 0; i < list->size; i++) + { + if (list->list[i].size != 0) + { + free(list->list[i].list); + } + } + + free(list->list); +} + + +void append_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length, +uint64_t n_end_pos, uint8_t n_direction, UC_Read* g_read, All_reads* R_INF, Correct_dumy* dumy) +{ + + if (list->length + 1 > list->size) + { + list->size = list->size * 2; + list->list = (k_mer_pos_list*)realloc(list->list, sizeof(k_mer_pos_list)*list->size); + ///新分配空间要初始化 + memset(list->list + (list->size/2), 0, sizeof(k_mer_pos_list)*(list->size/2)); + } + + if (list->list[list->length].size < n_length) + { + list->list[list->length].size = n_length; + list->list[list->length].list = (k_mer_pos*)realloc(list->list[list->length].list, + sizeof(k_mer_pos)*list->list[list->length].size); + } + + + ///list->list[list->length].list = n_list; + ///memcpy(list->list[list->length].list, n_list, sizeof(k_mer_pos)*n_length); + pre_filter_by_nearby(list->list[list->length].list, n_list, n_length, n_end_pos, + g_read, R_INF, dumy, &n_length); + + + if(n_length > 0) + { + list->list[list->length].length = n_length; + list->list[list->length].direction = n_direction; + list->list[list->length].end_pos = n_end_pos; + + list->length++; + } +} + + +/**********************for prefilter************************ */ diff --git a/Correct.h b/Correct.h index 225a4f6..7274977 100644 --- a/Correct.h +++ b/Correct.h @@ -8,6 +8,11 @@ #define THRESHOLD 15 #define GROUP_SIZE 4 +#define MAX(x, y) ((x >= y)?x:y) +#define MIN(x, y) ((x <= y)?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 + typedef struct { uint64_t* overlapID; @@ -26,6 +31,15 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Re 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 pre_filter_by_nearby(k_mer_pos* new_n_list, k_mer_pos* old_n_list, uint64_t n_length, uint64_t n_end_pos, UC_Read* g_read, +All_reads* R_INF, Correct_dumy* dumy, uint64_t* new_n_length); +void pre_filter_by_nearby_single(k_mer_pos* new_n_list, k_mer_pos* old_n_list, uint64_t n_length, uint64_t n_end_pos, UC_Read* g_read, +All_reads* R_INF, Correct_dumy* dumy, uint64_t* new_n_length); +/**********************for prefilter************************ */ +void destory_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list); +void append_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length, +uint64_t n_end_pos, uint8_t n_direction, UC_Read* g_read, All_reads* R_INF, Correct_dumy* dumy); +/**********************for prefilter************************ */ #endif \ No newline at end of file diff --git a/Hash_Table.cpp b/Hash_Table.cpp index b13bfc2..23e6eb6 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -317,6 +317,7 @@ void append_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp list->list[list->length].shared_seed = tmp->shared_seed; + list->list[list->length].align_length = 0; list->length++; } @@ -661,12 +662,8 @@ void init_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) { list->size = 1000; list->length = 0; - list->list = (k_mer_pos_list*)malloc(sizeof(k_mer_pos_list)*list->size); -} - -void destory_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) -{ - free(list->list); + //list->list = (k_mer_pos_list*)malloc(sizeof(k_mer_pos_list)*list->size); + list->list = (k_mer_pos_list*)calloc(list->size, sizeof(k_mer_pos_list)); } void clear_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) @@ -674,6 +671,10 @@ void clear_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) list->length = 0; } +void destory_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) +{ + free(list->list); +} void append_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length, uint64_t n_end_pos, uint8_t n_direction) @@ -693,6 +694,9 @@ uint64_t n_end_pos, uint8_t n_direction) list->length++; } + + + int cmp_k_mer_pos_list(const void * a, const void * b) { if ((*(k_mer_pos_list*)a).length > (*(k_mer_pos_list*)b).length) @@ -716,18 +720,7 @@ inline void append_pos_to_Candidates_list(Candidates_list* candidates, ElemType* candidates->length++; } -void pre_filter_by_nearby(k_mer_pos* n_list, uint64_t n_length, uint64_t n_end_pos, uint8_t n_direction, UC_Read* g_read) -{ - char* x_string = NULL; - char* y_string = NULL; - - long long i = 0; - for (i = 0; i < n_length; i++) - { - /* code */ - } - -} + void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, uint64_t strand) diff --git a/Hash_Table.h b/Hash_Table.h index bb0c68e..802061f 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -41,6 +41,7 @@ typedef struct { k_mer_pos* list; uint64_t length; + uint64_t size; uint8_t direction; uint64_t end_pos; } k_mer_pos_list; @@ -65,6 +66,7 @@ typedef struct uint64_t y_pos_strand; uint64_t shared_seed; + uint64_t align_length; } overlap_region; @@ -387,6 +389,9 @@ void destory_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list); void clear_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list); void append_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length, uint64_t n_end_pos, uint8_t n_direction); + + + void merge_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, Candidates_list* candidates); void merge_k_mer_pos_list_alloc_heap_sort(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT); void merge_k_mer_pos_list_alloc_heap_sort_advance(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT); @@ -443,6 +448,9 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF); + + +