diff --git a/Assembly.cpp b/Assembly.cpp index 6d828b3..c7e94e1 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -16,6 +16,9 @@ All_reads R_INF; pthread_mutex_t statistics; long long total_matched_overlap_0 = 0; long long total_matched_overlap_1 = 0; +long long total_potiental_matched_overlap_0 = 0; +long long total_potiental_matched_overlap_1 = 0; + long long complete_threads = 0; @@ -689,6 +692,8 @@ void* 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 j; int thr_ID = *((int*)arg); @@ -698,6 +703,10 @@ void* Overlap_calculate_heap_merge(void* arg) UC_Read g_read; init_UC_Read(&g_read); + + UC_Read overlap_read; + init_UC_Read(&overlap_read); + HPC_seq HPC_read; Hash_code k_code; uint64_t code; @@ -852,29 +861,30 @@ void* Overlap_calculate_heap_merge(void* arg) ///merge_k_mer_pos_list_alloc(&array_list, &l); merge_k_mer_pos_list_alloc_heap_sort(&array_list, &l, &heap); - /** - if (array_list.length < 3) - { - merge_k_mer_pos_list_alloc(&array_list, &l); - } - else - { - merge_k_mer_pos_list_alloc_heap_sort_advance(&array_list, &l, &heap); - } - **/ + + // if (array_list.length < 3) + // { + // merge_k_mer_pos_list_alloc(&array_list, &l); + // } + // else + // { + // merge_k_mer_pos_list_alloc_heap_sort_advance(&array_list, &l, &heap); + // } + ///以x_pos_e，即结束位置为主元排序 calculate_overlap_region(&l, &overlap_list, i, g_read.length, &R_INF); - correct_overlap(&overlap_list, &R_INF, &g_read, &correct); - - + correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, + &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1); + + /** 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 (Len_x * 0.9 <= overlap_list.list[j].align_length) { if (overlap_list.list[j].y_pos_strand == 0) { @@ -886,6 +896,7 @@ void* Overlap_calculate_heap_merge(void* arg) } } } + **/ /** POA_i = 0; @@ -988,6 +999,8 @@ void* Overlap_calculate_heap_merge(void* arg) destory_Graph(&POA_Graph); destory_UC_Read(&g_read); + destory_UC_Read(&overlap_read); + destory_Correct_dumy(&correct); @@ -995,11 +1008,15 @@ void* Overlap_calculate_heap_merge(void* arg) 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; 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); + fprintf(stderr, "total_potiental_matched_overlap_0: %llu\n", total_potiental_matched_overlap_0); + fprintf(stderr, "total_potiental_matched_overlap_1: %llu\n", total_potiental_matched_overlap_1); } pthread_mutex_unlock(&statistics); } diff --git a/Correct.cpp b/Correct.cpp index 6f820a8..311ef7d 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -674,7 +674,11 @@ char* r_string) } -void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy) + +void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, + UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, + long long* matched_overlap_0, long long* matched_overlap_1, + long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1) { reverse_complement(g_read->seq, g_read->length); @@ -733,6 +737,86 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Re // T_total_match, T_total_unmatch, T_total_mis); // pthread_mutex_unlock(&debug_statistics); /************需要注释掉********* */ + /** + long long j; + long long Len_x; + int threshold; + long long y_start; + long long Len_y; + long long currentIDLen; + for (j = 0; j < overlap_list->length; j++) + { + 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) + { + threshold = Len_x * 0.04; + + y_start = overlap_list->list[j].y_pos_s - threshold; + if (y_start < 0) + { + y_start = 0; + } + + + ///当前y的长度 + currentIDLen = Get_READ_LENGTH((*R_INF), overlap_list->list[j].y_id); + ///不能超过y的剩余长度 + Len_y = MIN(Len_x + 2 * threshold, currentIDLen - y_start); + + + if (overlap_list->list[j].y_pos_strand == 0) + { + recover_UC_Read(overlap_read, R_INF, overlap_list->list[j].y_id); + (*matched_overlap_0)++; + } + else + { + recover_UC_Read_RC(overlap_read, R_INF, overlap_list->list[j].y_id); + (*matched_overlap_1)++; + } + + + + + EdlibAlignResult result = edlibAlign(g_read->seq + overlap_list->list[j].x_pos_s, Len_x, + overlap_read->seq + y_start, Len_y, + edlibNewAlignConfig(threshold, EDLIB_MODE_HW, EDLIB_TASK_PATH, NULL, 0)); + + + if (result.status == EDLIB_STATUS_OK && result.editDistance != -1) { + char* cigar = edlibAlignmentToCigar(result.alignment, result.alignmentLength, EDLIB_CIGAR_STANDARD); + int cigar_length = strlen(cigar); + free(cigar); + + + if (overlap_list->list[j].y_pos_strand == 0) + { + (*potiental_matched_overlap_0)++; + } + else + { + (*potiental_matched_overlap_1)++; + } + + // if (overlap_list->list[j].shared_seed == 1 && Len_x >= 1000) + // { + // (*matched_overlap_0)++; + // } + + // if (overlap_list->list[j].shared_seed == 2 && Len_x >= 1000) + // { + // (*matched_overlap_1)++; + // } + + + } + + edlibFreeAlignResult(result); + } + } + **/ + } diff --git a/Correct.h b/Correct.h index 7274977..1611afb 100644 --- a/Correct.h +++ b/Correct.h @@ -26,8 +26,10 @@ typedef struct } Correct_dumy; -void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy); - +void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, + UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, + long long* matched_overlap_0, long long* matched_overlap_1, + long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1); 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); diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index 2bdc58a..93734bc 100644 --- a/Levenshtein_distance.h +++ b/Levenshtein_distance.h @@ -228,6 +228,386 @@ inline int Reserve_Banded_BPM } + + + +inline int Reserve_Banded_BPM_PATH +(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, + unsigned int* return_err, int* return_start_site, int* return_path_length, Word* matrix_bit, char* path) +{ + (*return_err) = (unsigned int)-1; + + Word Peq[256]; + + int band_length = (errthold << 1) + 1; + int i = 0; + Word tmp_Peq_1 = (Word)1; + + Peq['A'] = (Word)0; + Peq['T'] = (Word)0; + Peq['G'] = (Word)0; + Peq['C'] = (Word)0; + + + Word Peq_A; + Word Peq_T; + Word Peq_C; + Word Peq_G; + + ///band_length = 2k + 1 + for (i = 0; i> 1; + VN = X&HP; + VP = HN | ~(X | HP); + + if (!(D0&err_mask)) + { + ++err; + + ///¼´Ê¹È«²¿µÝ¼õ£¬Ò²¾Í¼õ2k + if ((err - last_high)>errthold) + { + ///fprintf(stderr, "0 ######, i: %u\n", i); + return -1; + } + + } + + + Peq['A'] = Peq['A'] >> 1; + Peq['C'] = Peq['C'] >> 1; + Peq['G'] = Peq['G'] >> 1; + Peq['T'] = Peq['T'] >> 1; + + + ++i; + ++i_bd; + Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask; + + + ///Peq['T'] = Peq['T'] | Peq['C']; + + column_start = i << 3; + matrix_bit[column_start] = D0; + matrix_bit[column_start + 1] = VP; + matrix_bit[column_start + 2] = VN; + matrix_bit[column_start + 3] = HP; + matrix_bit[column_start + 4] = HN; + } + + + + + + X = Peq[text[i]] | VN; + D0 = ((VP + (X&VP)) ^ VP) | X; + HN = VP&D0; + HP = VN | ~(VP | D0); + X = D0 >> 1; + VN = X&HP; + VP = HN | ~(X | HP); + if (!(D0&err_mask)) + { + ++err; + if ((err - last_high)>errthold) + return -1; + } + + + column_start = (i + 1) << 3; + matrix_bit[column_start] = D0; + matrix_bit[column_start + 1] = VP; + matrix_bit[column_start + 2] = VN; + matrix_bit[column_start + 3] = HP; + matrix_bit[column_start + 4] = HN; + + + ////fprintf(stderr, "sucess(2)\n"); + + /// last_high = 2k + /// site = (SEQ_LENGTH + 2k) - 2k -1 + /// site = SEQ_LENGTH - 1 + ///int site = p_length - last_high - 1; + int site = t_length - 1; + int return_site = -1; + if ((err <= errthold) && (err<=*return_err)) + { + *return_err = err; + return_site = site; + } + int i_last = i; + i = 0; + + + + + while (i> i)&(Word)1); + err = err - ((VN >> i)&(Word)1); + ++i; + + if ((err <= errthold) && (err <= *return_err)) + { + *return_err = err; + return_site = site + i; + } + } + + + unsigned int ungap_err; + ungap_err = err; + + + while (i> i)&(Word)1); + err = err - ((VN >> i)&(Word)1); + ++i; + + if ((err <= errthold) && (err<=*return_err)) + { + *return_err = err; + return_site = site + i; + } + + + + + } + + + if ((ungap_err <= errthold) && (ungap_err == *return_err)) + { + return_site = site + errthold; + } + + if ((*return_err) == (unsigned int)-1) + { + return return_site; + } + + + + ///end_site是正确的 + int end_site = return_site; + int start_site = end_site; + ///这个是各个bit-vector里面，end_site对应bit所在的位置 + int back_track_site = band_length - (p_length - end_site); + + Word v_value, h_value, delta_value, min_value, current_value; + Word direction, is_mismatch; ///0 is match, 1 is mismatch, 2 is up, 3 is left + + ///代表pattern到哪了，就是短的那个到哪了 + i = t_length; + int path_length = 0; + + ///到0就结束了，后面的路径可以直接match + current_value = *return_err; + + + int low_bound = band_length - 1; + + + while (i>0) + { + if (current_value == 0) + { + break; + } + + column_start = i << 3; + + delta_value = current_value - + ((~(matrix_bit[column_start] >> back_track_site))&err_mask); + + + if (back_track_site == 0) + { + ///HP + h_value = current_value - ((matrix_bit[column_start + 3] >> back_track_site)&err_mask); + //HN + h_value = h_value + ((matrix_bit[column_start + 4] >> back_track_site)&err_mask); + + + min_value = delta_value; + direction = 0; + + if (h_value < min_value) + { + min_value = h_value; + direction = 3; + } + + } + else if (back_track_site == low_bound) + { + v_value = current_value - ((matrix_bit[column_start + 1] >> (back_track_site - 1))&err_mask); + v_value = v_value + ((matrix_bit[column_start + 2] >> (back_track_site - 1))&err_mask); + + min_value = delta_value; + direction = 0; + + if (v_value < min_value) + { + min_value = v_value; + direction = 2; + } + + } + else + { + + h_value = current_value - ((matrix_bit[column_start + 3] >> back_track_site)&(Word)1); + + + h_value = h_value + ((matrix_bit[column_start + 4] >> back_track_site)&(Word)1); + + + v_value = current_value - ((matrix_bit[column_start + 1] >> (back_track_site - 1))&err_mask); + v_value = v_value + ((matrix_bit[column_start + 2] >> (back_track_site - 1))&err_mask); + + + min_value = delta_value; + direction = 0; + + if (v_value < min_value) + { + min_value = v_value; + direction = 2; + } + + + if (h_value < min_value) + { + min_value = h_value; + direction = 3; + } + + } + + + if (direction == 0) + { + + if (delta_value != current_value) + { + direction = 1; + } + + i--; + + start_site--; + + } + if (direction == 2)///ru guo xiang shang yi dong, bing bu huan lie + { + back_track_site--; + start_site--; + } + else if (direction == 3)///ru guo xiang zuo yi dong + { + i--; + back_track_site++; + } + + + path[path_length++] = direction; + + + current_value = min_value; + + } + + + if (i > 0) + { + memset(path + path_length, 0, i); + start_site = start_site - i; + direction = 0; + path_length = path_length + i; + } + + if (direction != 3) + { + start_site++; + } + (*return_start_site) = start_site; + (*return_path_length) = path_length; + + + + + return return_site; + +} + inline int Reserve_Banded_BPM_4_SSE_only(char *pattern1, char *pattern2, char *pattern3, char *pattern4, int p_length, char *text, int t_length, int* return_sites, unsigned int* return_sites_error, unsigned short errthold, __m128i* Peq_SSE)