diff --git a/.vscode/ipch/2549b0c830b4882d/POA.ipch b/.vscode/ipch/2549b0c830b4882d/POA.ipch deleted file mode 100644 index a84cd9b..0000000 Binary files a/.vscode/ipch/2549b0c830b4882d/POA.ipch and /dev/null differ diff --git a/.vscode/ipch/2549b0c830b4882d/mmap_address.bin b/.vscode/ipch/2549b0c830b4882d/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/2549b0c830b4882d/mmap_address.bin and /dev/null differ diff --git a/.vscode/ipch/491621f0e7654f5/mmap_address.bin b/.vscode/ipch/491621f0e7654f5/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/491621f0e7654f5/mmap_address.bin and /dev/null differ diff --git a/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch b/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch deleted file mode 100644 index a33202d..0000000 Binary files a/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch and /dev/null differ diff --git a/.vscode/ipch/4eee59091dc293a3/mmap_address.bin b/.vscode/ipch/4eee59091dc293a3/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/4eee59091dc293a3/mmap_address.bin and /dev/null differ diff --git a/.vscode/ipch/756c62c8d5e4771d/COMMANDLINES.ipch b/.vscode/ipch/756c62c8d5e4771d/COMMANDLINES.ipch deleted file mode 100644 index 895d32c..0000000 Binary files a/.vscode/ipch/756c62c8d5e4771d/COMMANDLINES.ipch and /dev/null differ diff --git a/.vscode/ipch/756c62c8d5e4771d/mmap_address.bin b/.vscode/ipch/756c62c8d5e4771d/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/756c62c8d5e4771d/mmap_address.bin and /dev/null differ diff --git a/.vscode/ipch/856c69690c514632/KMER.ipch b/.vscode/ipch/856c69690c514632/KMER.ipch deleted file mode 100644 index 23c0c45..0000000 Binary files a/.vscode/ipch/856c69690c514632/KMER.ipch and /dev/null differ diff --git a/.vscode/ipch/856c69690c514632/mmap_address.bin b/.vscode/ipch/856c69690c514632/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/856c69690c514632/mmap_address.bin and /dev/null differ diff --git a/.vscode/ipch/920fff44ce64d7c2/mmap_address.bin b/.vscode/ipch/920fff44ce64d7c2/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/920fff44ce64d7c2/mmap_address.bin and /dev/null differ diff --git a/.vscode/ipch/a64815cd9ef3d0a4/mmap_address.bin b/.vscode/ipch/a64815cd9ef3d0a4/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/a64815cd9ef3d0a4/mmap_address.bin and /dev/null differ diff --git a/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch b/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch deleted file mode 100644 index b11e08a..0000000 Binary files a/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch and /dev/null differ diff --git a/.vscode/ipch/b49c1e8b711f08aa/mmap_address.bin b/.vscode/ipch/b49c1e8b711f08aa/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/b49c1e8b711f08aa/mmap_address.bin and /dev/null differ diff --git a/.vscode/ipch/b5c6b274f2611ee7/KMER.ipch b/.vscode/ipch/b5c6b274f2611ee7/KMER.ipch deleted file mode 100644 index 9221230..0000000 Binary files a/.vscode/ipch/b5c6b274f2611ee7/KMER.ipch and /dev/null differ diff --git a/.vscode/ipch/b5c6b274f2611ee7/mmap_address.bin b/.vscode/ipch/b5c6b274f2611ee7/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/b5c6b274f2611ee7/mmap_address.bin and /dev/null differ diff --git a/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch b/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch deleted file mode 100644 index fcf8a67..0000000 Binary files a/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch and /dev/null differ diff --git a/.vscode/ipch/c0e71cfe49f0fe81/mmap_address.bin b/.vscode/ipch/c0e71cfe49f0fe81/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/c0e71cfe49f0fe81/mmap_address.bin and /dev/null differ diff --git a/.vscode/ipch/e64bc513afbf8204/mmap_address.bin b/.vscode/ipch/e64bc513afbf8204/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/e64bc513afbf8204/mmap_address.bin and /dev/null differ diff --git a/.vscode/ipch/efd6dc00fdc3a0f/HASH_TABLE.ipch b/.vscode/ipch/efd6dc00fdc3a0f/HASH_TABLE.ipch deleted file mode 100644 index 5e161a1..0000000 Binary files a/.vscode/ipch/efd6dc00fdc3a0f/HASH_TABLE.ipch and /dev/null differ diff --git a/.vscode/ipch/efd6dc00fdc3a0f/mmap_address.bin b/.vscode/ipch/efd6dc00fdc3a0f/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/efd6dc00fdc3a0f/mmap_address.bin and /dev/null differ diff --git a/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch b/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch deleted file mode 100644 index bb86d9c..0000000 Binary files a/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch and /dev/null differ diff --git a/.vscode/ipch/f2b98741f94c3e10/mmap_address.bin b/.vscode/ipch/f2b98741f94c3e10/mmap_address.bin deleted file mode 100644 index 862b842..0000000 Binary files a/.vscode/ipch/f2b98741f94c3e10/mmap_address.bin and /dev/null differ diff --git a/.vscode/settings.json b/.vscode/settings.json index bf562e6..18a7993 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,3 +1,8 @@ { - "C_Cpp.errorSquiggles": "Enabled" + "C_Cpp.errorSquiggles": "Enabled", + "files.associations": { + "limits": "cpp", + "random": "cpp", + "functional": "cpp" + } } \ No newline at end of file diff --git a/Assembly.cpp b/Assembly.cpp index 128ac7b..e0cc482 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -7,6 +7,7 @@ #include "kmer.h" #include "Hash_Table.h" #include "POA.h" +#include "Correct.h" Total_Count_Table TCB; Total_Pos_Table PCB; @@ -693,6 +694,9 @@ void* Overlap_calculate_heap_merge(void* arg) 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; @@ -712,16 +716,27 @@ void* Overlap_calculate_heap_merge(void* arg) Init_Heap(&heap); + Correct_dumy correct; + init_Correct_dumy(&correct); + for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num) { /** - if (thr_ID == 0 && i % 1000 == 0) + if (i < 4) { - fprintf(stderr, "i: %llu\n", i); + continue; } + + if (i > 4) + { + break; + } + fprintf(stderr, "i: %u\n", i); + fflush(stderr); **/ + clear_Heap(&heap); clear_Candidates_list(&l); ///clear_Candidates_list(&debug_l); @@ -833,8 +848,48 @@ void* Overlap_calculate_heap_merge(void* arg) } **/ + ///以x_pos_e，即结束位置为主元排序 calculate_overlap_region(&l, &overlap_list, i, g_read.length, &R_INF); + + correct_overlap(&overlap_list, &R_INF, &g_read, &correct); + + + + /** + POA_i = 0; + fprintf(stderr, "\n\n**************\ni: %u\n", i); + + for (POA_i = 0; POA_i < overlap_list.length; POA_i++) + { + fprintf(stderr, "x_id: %u, y_id: %u\n", overlap_list.list[POA_i].x_id, overlap_list.list[POA_i].y_id); + fprintf(stderr, "x_strand: %u, y_strand: %u\n", overlap_list.list[POA_i].x_pos_strand, overlap_list.list[POA_i].y_pos_strand); + fprintf(stderr, "x_pos_s: %u\n", overlap_list.list[POA_i].x_pos_s); + + } + **/ + + /** + POA_i = 0; + + for (POA_i = 1; POA_i < overlap_list.length; POA_i++) + { + if(overlap_list.list[POA_i].x_pos_s < overlap_list.list[POA_i - 1].x_pos_s) + { + fprintf(stderr, "1 sbsbsbsbs\n"); + } + else if(overlap_list.list[POA_i].x_pos_s == overlap_list.list[POA_i - 1].x_pos_s && + overlap_list.list[POA_i].x_pos_e > overlap_list.list[POA_i - 1].x_pos_e) + { + fprintf(stderr, "2 sbsbsbsbs\n"); + } + + + } + **/ + + + ///merge_k_mer_pos_list_alloc_heap_sort_advance(&array_list, &l, &heap); @@ -843,10 +898,12 @@ void* Overlap_calculate_heap_merge(void* arg) ///debug_merge_result(&l, &debug_l); + /** clear_Graph(&POA_Graph); Perform_POA(&POA_Graph, &overlap_list, &R_INF, &g_read); + **/ ///fprintf(stderr, "i: %u\n", i); @@ -875,7 +932,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); @@ -885,6 +946,8 @@ void* Overlap_calculate_heap_merge(void* arg) destory_Graph(&POA_Graph); destory_UC_Read(&g_read); + + destory_Correct_dumy(&correct); } diff --git a/Correct.cpp b/Correct.cpp new file mode 100644 index 0000000..c8970c6 --- /dev/null +++ b/Correct.cpp @@ -0,0 +1,685 @@ +#include +#include +#include +#include "Correct.h" +#include "Levenshtein_distance.h" +#include "edlib.h" + +long long T_total_match=0; +long long T_total_unmatch=0; +long long T_total_mis=0; + + +#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 + + + +///y_length > x_length +unsigned int edit_distance_normal_test_banded(char* y, int y_length, char* x, int x_length, int error_cut, int matrix[1000][1000] ) +{ memset(matrix, 0, sizeof(matrix)); + + int i, j; + for (i = 0; i <= x_length; i++) + { + matrix[i][0] = i; + } + + int digonal, up, left; + unsigned int min; + + ///一列列算的 + for (i = 0; i < x_length; i++) + { + for (j = 0; j < y_length; j++) + { + ///matrix[i + 1][j + 1] + digonal = matrix[i][j] + (x[i] != y[j]); + up = matrix[i + 1][j] + 1; + left = matrix[i][j + 1] + 1; + min = digonal; + if (up < min) + { + min = up; + } + + if (left< min) + { + min = left; + } + + matrix[i + 1][j + 1] = min; + } + } + + min = (unsigned int)-1; + for (j = x_length; j <= y_length; j++) + { + if (matrix[i][j] < min) + { + min = matrix[i][j]; + } + } + + + return min <= error_cut?min:(unsigned int)(-1); +} + +void verify_get_interval(long long window_start, long long window_end, overlap_region_alloc* overlap_list,Correct_dumy* dumy) +{ + long long i; + long long match_length = 0; + long long match_lengthNT = 0; + long long Len; + + for (i = 0; i < overlap_list->length; i++) + { + if((Len = OVERLAP(window_start, window_end, overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e)) > 0) + { + if (Len == WINDOW) + { + match_length++; + long long j; + for (j = 0; j < dumy->length; j++) + { + if (i==dumy->overlapID[j]) + { + break; + } + } + + if (j >= dumy->length) + { + fprintf(stderr, "+ERROR interval\n"); + + fprintf(stderr, "i: %u, window_start: %u, window_end: %u, x_pos_s: %u, x_pos_e: %u\n", + i, window_start, window_end, overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e); + } + } + else + { + match_lengthNT++; + long long j; + for (j = 0; j < dumy->lengthNT; j++) + { + if (i==dumy->overlapID[dumy->size - j - 1]) + { + break; + } + } + + if (j >= dumy->lengthNT) + { + fprintf(stderr, "-ERROR interval\n"); + + fprintf(stderr, "i: %u, window_start: %u, window_end: %u, x_pos_s: %u, x_pos_e: %u, dumy->lengthNT: %u\n", + i, window_start, window_end, overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e, dumy->lengthNT); + } + } + + + + + + } + } + + if (match_length != dumy->length || match_lengthNT != dumy->lengthNT) + { + fprintf(stderr, "****************ERROR interval length*******************\n"); + fprintf(stderr, "match_length: %u\n", match_length); + fprintf(stderr, "dumy->length: %u\n", dumy->length); + fprintf(stderr, "window_start: %u, window_end: %u\n", window_start, window_end); + } + +} + + +inline int get_interval_back(long long window_start, long long window_end, overlap_region_alloc* overlap_list, Correct_dumy* dumy) +{ + long long i; + int flag = 0; + + for (i = dumy->start_i; i < overlap_list->length; i++) + { + ///只会发生在这个interval比list里所有元素都小的情况 + ///这种情况下一个interval需要从0开始 + if (window_start < overlap_list->list[i].x_pos_s) + { + dumy->start_i = 0; + dumy->length = 0; + return -1; + } + else if(window_start >= overlap_list->list[i].x_pos_s && window_start <= overlap_list->list[i].x_pos_e) + { + dumy->start_i = i; + break; + } + } + + ///只会发生在这个window比list里所有元素都大的情况 + ///这种情况下一个window也无需遍历了 + if (i >= overlap_list->length) + { + dumy->start_i = overlap_list->length; + dumy->length = 0; + return -2; + } + + ///走到这里的时候，至少window_start的要求是满足了 + dumy->length = 0; + + for (; i < overlap_list->length; i++) + { + if(overlap_list->list[i].x_pos_s <= window_start && overlap_list->list[i].x_pos_e >= window_end) + { + dumy->overlapID[dumy->length] = i; + dumy->length++; + } + else if(overlap_list->list[i].x_pos_s > window_start) + { + break; + } + } + + if ( dumy->length == 0) + { + return 0; + } + else + { + return 1; + } +} + + + +inline int get_interval(long long window_start, long long window_end, overlap_region_alloc* overlap_list, Correct_dumy* dumy) +{ + long long i; + int flag = 0; + long long Begin, End, Len; + + + for (i = dumy->start_i; i < overlap_list->length; i++) + { + ///只会发生在这个interval比list里所有元素都小的情况 + ///这种情况下一个interval需要从0开始 + if (window_end < overlap_list->list[i].x_pos_s) + { + dumy->start_i = 0; + dumy->length = 0; + return 0; + } + else ///只要window_end >= overlap_list->list[i].x_pos_s，就有可能重叠 + { + dumy->start_i = i; + break; + } + } + + ///只会发生在这个window比list里所有元素都大的情况 + ///这种情况下一个window也无需遍历了 + if (i >= overlap_list->length) + { + dumy->start_i = overlap_list->length; + dumy->length = 0; + return -2; + } + + dumy->length = 0; + dumy->lengthNT = 0; + + for (; i < overlap_list->length; i++) + { + + if((Len = OVERLAP(window_start, window_end, overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e)) > 0) + { + if (Len == WINDOW) + { + dumy->overlapID[dumy->length] = i; + dumy->length++; + } + else + { + dumy->lengthNT++; + dumy->overlapID[dumy->size - dumy->lengthNT] = i; + } + } + + if(overlap_list->list[i].x_pos_s > window_end) + { + break; + } + } + + if ( dumy->length + dumy->lengthNT == 0) + { + return 0; + } + else + { + return 1; + } +} + +///Len = OVERLAP(window_start, window_end, overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e)) + +void print_string(char* s, int l) +{ + for (size_t i = 0; i < l; i++) + { + fprintf(stderr, "%c", s[i]); + } + + fprintf(stderr, "\n"); + +} + +void test_edit_distance_by_edlib(char* x_string, char* y_string, long long x_len, +long long o_len, int threashold, int error, long long* total_mis) +{ + EdlibAlignResult result = edlibAlign(x_string, x_len, y_string, o_len, + edlibNewAlignConfig(threashold, EDLIB_MODE_HW, EDLIB_TASK_PATH, NULL, 0)); + + if (result.status == EDLIB_STATUS_OK) { + if (result.editDistance != error) + ///if (result.editDistance != error && result.endLocations[0] > x_len) + { + + (*total_mis)++; + + + if ((int)error != -1 && result.editDistance==-1) + { + fprintf(stderr, "ERROR1\n"); + } + + if ((int)error < result.editDistance==-1 && + (int)error != -1 && result.editDistance!=-1) + { + fprintf(stderr, "ERROR2\n"); + } + + int up_length = o_len - result.endLocations[0] - 1; + int left_length = result.endLocations[0] - x_len; + + + char* cigar = edlibAlignmentToCigar(result.alignment, result.alignmentLength, EDLIB_CIGAR_STANDARD); + int cigar_length = strlen(cigar); + + /** + for (int i = cigar_length - 1; i >= 0; i--) + { + switch (cigar[i]) + { + case 'I': + + break; + + default: + break; + } + } + **/ + + ///fprintf(stderr,"%s\n", cigar); + free(cigar); + + /** + fprintf(stderr, "****\ni: %u, edlib: %d, alignmentLength: %d, startLocations: %d, endLocations: %d\n", + i, result.editDistance, result.alignmentLength, result.startLocations[0], result.endLocations[0]); + char* cigar = edlibAlignmentToCigar(result.alignment, result.alignmentLength, EDLIB_CIGAR_STANDARD); + fprintf(stderr,"%s\n", cigar); + free(cigar); + print_string(x_string, x_len); + print_string(y_string, o_len); + + fprintf(stderr, "BPM: %d\n", error); + **/ + + } + } + edlibFreeAlignResult(result); +} + +void verify_window(long long window_start, long long window_end, overlap_region_alloc* overlap_list,Correct_dumy* dumy, All_reads* R_INF, +char* r_string) +{ + + long long i; + long long currentID, currentIDLen; + long long x_start, y_start, o_len; + long long Window_Len = WINDOW + (THRESHOLD << 1); + char* x_string = NULL; + char* y_string = NULL; + long long x_end, x_len; + int end_site; + unsigned int error; + long long total_match=0; + long long total_unmatch=0; + long long total_mis=0; + + + ///这些是整个window被完全覆盖的 + for (i = 0; i < dumy->length; i++) + { + ///整个window被覆盖的话，read本身上的区间就是[window_start, window_end] + x_len = WINDOW; + currentID = dumy->overlapID[i]; + x_start = window_start; + ///y上的相对位置 + y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s; + ///y上的起始 + y_start = y_start - THRESHOLD; + if (y_start < 0) + { + y_start = 0; + } + ///当前y的长度 + currentIDLen = Get_READ_LENGTH((*R_INF), overlap_list->list[currentID].y_id); + ///不能超过y的剩余长度 + o_len = MIN(Window_Len, currentIDLen - y_start); + + 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); + + x_string = r_string + x_start; + y_string = dumy->overlap_region; + + + + + ///这个长度足够，可以用来一起比 + if (Window_Len == o_len) + { + + end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, THRESHOLD, &error); + + + if (error!=(unsigned int)-1) + { + total_match++; + } + else + { + total_unmatch++; + } + + + test_edit_distance_by_edlib(x_string, y_string, x_len, + o_len, THRESHOLD, error, &total_mis); + + } + else ///这个不够，只能单个比 ///不够的地方要置N + { + /** + end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, THRESHOLD, &error); + + + if (error!=-1) + { + total_match++; + } + else + { + total_unmatch++; + } + **/ + } + } + + long long reverse_i = dumy->size - 1; + + ///这些是整个window被部分覆盖的 + for (i = 0; i < dumy->lengthNT; i++) + { + currentID = dumy->overlapID[reverse_i--]; + x_start = MAX(window_start, overlap_list->list[currentID].x_pos_s); + x_end = MIN(window_end, overlap_list->list[currentID].x_pos_e); + ///这个是和当前窗口重叠的长度 + x_len = x_end - x_start + 1; + + 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; + ///y上的起始 + y_start = y_start - THRESHOLD; + if (y_start < 0) + { + y_start = 0; + } + + ///当前y的长度 + currentIDLen = Get_READ_LENGTH((*R_INF), overlap_list->list[currentID].y_id); + + + ///不能超过y的剩余长度 + Window_Len = x_len + (THRESHOLD << 1); + o_len = MIN(Window_Len, currentIDLen - y_start); + + 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); + + x_string = r_string + x_start; + y_string = dumy->overlap_region; + + ///不够的地方要置N + + end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, THRESHOLD, &error); + + /** + if (error!=-1) + { + total_match++; + } + else + { + total_unmatch++; + } + **/ + } + + /** + if (total_match!=0) + { + + fprintf(stderr, "total_match: %u\n", total_match); + fprintf(stderr, "total_unmatch: %u\n", total_unmatch); + fprintf(stderr, "total_mis: %u\n", total_mis); + + } + **/ + + T_total_match = T_total_match + total_match; + T_total_unmatch = T_total_unmatch + total_unmatch; + T_total_mis = T_total_mis + total_mis; + + + + + + +} + +void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy) +{ + reverse_complement(g_read->seq, g_read->length); + + clear_Correct_dumy(dumy, overlap_list); + + long long window_num = (g_read->length + WINDOW - 1) / WINDOW; + + long long i; + + long long window_start, window_end; + + window_start = 0; + window_end = WINDOW - 1; + + + + /** + UC_Read debug_read; + init_UC_Read(&debug_read); + + for (i = 0; i < overlap_list->length; i++) + { + + if (overlap_list->list[i].y_pos_strand) + { + recover_UC_Read_RC(&debug_read, R_INF, overlap_list->list[i].y_id); + } + else + { + recover_UC_Read(&debug_read, R_INF, overlap_list->list[i].y_id); + } + + + long long x_length = overlap_list->list[i].x_pos_e - overlap_list->list[i].x_pos_s + 1; + long long y_length = overlap_list->list[i].y_pos_e - overlap_list->list[i].y_pos_s + 1; + + EdlibAlignResult result = edlibAlign(g_read->seq+overlap_list->list[i].x_pos_s, + x_length, + debug_read.seq + overlap_list->list[i].y_pos_s, + y_length, + edlibNewAlignConfig(-1, EDLIB_MODE_NW, EDLIB_TASK_DISTANCE, NULL, 0)); + + + if (result.status == EDLIB_STATUS_OK) { + if (result.editDistance>0 && result.editDistance < x_length*0.02) + { + fprintf(stderr, "inner_i: %u, x_length: %u, y_length: %u, editDistance: %u\n", + i, x_length, y_length, result.editDistance); + fprintf(stderr, "x_pos_s: %u, x_pos_e: %u\n", + overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e); + fprintf(stderr, "y_pos_s: %u, y_pos_e: %u, y_pos_strand: %u\n", + overlap_list->list[i].y_pos_s, overlap_list->list[i].y_pos_e, overlap_list->list[i].y_pos_strand); + + + EdlibAlignResult n_result = edlibAlign(g_read->seq+overlap_list->list[i].x_pos_s, + 350, + debug_read.seq + overlap_list->list[i].y_pos_s - 15, + 380, + edlibNewAlignConfig(-1, EDLIB_MODE_HW, EDLIB_TASK_DISTANCE, NULL, 0)); + + fprintf(stderr, "n_editDistance: %u\n", + n_result.editDistance); + + + + + edlibFreeAlignResult(n_result); + } + + + } + edlibFreeAlignResult(result); + + + } + + destory_UC_Read(&debug_read); + **/ + ///return; + + + + + + + + + + + + + + + + + + + int flag; + + for (i = 0; i < window_num; i++) + { + + dumy->length = 0; + dumy->lengthNT = 0; + flag = get_interval(window_start, window_end, overlap_list, dumy); + + switch (flag) + { + case 1: ///找到匹配 + break; + case 0: ///没找到匹配 + break; + case -2: ///下一个window也不会存在匹配, 直接跳出 + i = window_num; + break; + } + + if(dumy->length + dumy->lengthNT>overlap_list->length) + { + fprintf(stderr, "error length\n"); + } + + ///verify_get_interval(window_start, window_end, overlap_list, dumy); + + verify_window(window_start, window_end, overlap_list, dumy, R_INF, g_read->seq); + + window_start = window_start + WINDOW; + window_end = window_end + WINDOW; + if (window_end >= g_read->length) + { + window_end = g_read->length - 1; + } + + ///break; + + } + + + fprintf(stderr, "total_match: %u, total_unmatch: %u, total_mis: %u\n", + T_total_match, T_total_unmatch, T_total_mis); + + + + +} + + +void init_Correct_dumy(Correct_dumy* list) +{ + list->size = 0; + list->length = 0; + list->lengthNT = 0; + list->start_i = 0; + list->overlapID = NULL; +} + +void destory_Correct_dumy(Correct_dumy* list) +{ + free(list->overlapID); +} + +void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list) +{ + list->length = 0; + list->lengthNT = 0; + list->start_i = 0; + + if (list->size < overlap_list->length) + { + list->size = overlap_list->length; + list->overlapID = (uint64_t*)realloc(list->overlapID, list->size*sizeof(uint64_t)); + } + +} \ No newline at end of file diff --git a/Correct.h b/Correct.h new file mode 100644 index 0000000..3990ac3 --- /dev/null +++ b/Correct.h @@ -0,0 +1,27 @@ +#ifndef __CORRECT__ +#define __CORRECT__ +#include +#include "Hash_Table.h" + +#define WINDOW 350 +#define THRESHOLD 15 + +typedef struct +{ + uint64_t* overlapID; + uint64_t length; + uint64_t lengthNT; + uint64_t size; + uint64_t start_i; + char overlap_region[WINDOW + THRESHOLD*2 + 10]; +} Correct_dumy; + + +void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy); + +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); + + +#endif \ No newline at end of file diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 9ddae36..02fda01 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -483,6 +483,64 @@ void print_overlap_region(Candidates_list* candidates, overlap_region_alloc* ove } +int cmp_by_x_pos_s(const void * a, const void * b) +{ + if ((*(overlap_region*)a).x_pos_s > (*(overlap_region*)b).x_pos_s) + { + return 1; + } + else if ((*(overlap_region*)a).x_pos_s < (*(overlap_region*)b).x_pos_s) + { + return -1; + } + else + { + + if ((*(overlap_region*)a).x_pos_e > (*(overlap_region*)b).x_pos_e) + { + return 1; + } + else if ((*(overlap_region*)a).x_pos_e < (*(overlap_region*)b).x_pos_e) + { + return -1; + } + else + { + return 0; + } + + } +} + + +int cmp_by_x_pos_e(const void * a, const void * b) +{ + if ((*(overlap_region*)a).x_pos_e > (*(overlap_region*)b).x_pos_e) + { + return 1; + } + else if ((*(overlap_region*)a).x_pos_e < (*(overlap_region*)b).x_pos_e) + { + return -1; + } + else + { + + if ((*(overlap_region*)a).x_pos_s > (*(overlap_region*)b).x_pos_s) + { + return 1; + } + else if ((*(overlap_region*)a).x_pos_s < (*(overlap_region*)b).x_pos_s) + { + return -1; + } + else + { + return 0; + } + + } +} ///r->length = Get_READ_LENGTH((*R_INF), ID); @@ -500,6 +558,7 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF) long long constant_distance = 5; double error_rate = 0.05; + if (candidates->length == 0) { return; @@ -556,14 +615,19 @@ 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_overlap_region_alloc(overlap_list, &tmp_region, R_INF); } } + ///以x_pos_e，即结束位置为主元排序 + ///qsort(overlap_list->list, overlap_list->length, sizeof(overlap_region), cmp_by_x_pos_e); + qsort(overlap_list->list, overlap_list->length, sizeof(overlap_region), cmp_by_x_pos_s); ///debug_overlap_region(candidates, overlap_list, readID); ///print_overlap_region(candidates, overlap_list, R_INF); diff --git a/Levenshtein_distance.cpp b/Levenshtein_distance.cpp new file mode 100644 index 0000000..4622c9d --- /dev/null +++ b/Levenshtein_distance.cpp @@ -0,0 +1 @@ +#include "Levenshtein_distance.h" diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h new file mode 100644 index 0000000..cd8d239 --- /dev/null +++ b/Levenshtein_distance.h @@ -0,0 +1,1050 @@ +#ifndef __LEVENSHTEIN__ +#define __LEVENSHTEIN__ +#include +#include "emmintrin.h" +#include "nmmintrin.h" +#include "smmintrin.h" +#include +#include +#include +#include + +typedef uint64_t Word; + + + + +inline void output_bit_myers(Word x, int length) +{ + int i = 0; + while (i < length) + { + fprintf(stderr, "%u", (x >> i) & ((Word)1)); + i++; + } + fprintf(stderr, "\n"); +} + + +inline void prase_vertical(Word VP, Word VN, int length, int matrix[1000][1000], int i) +{ + + i++; + int k = 0; + int j = i; + int diff; + while (k < length) + { + + int x_p = (VP >> k) & ((Word)1); + int x_n = (VN >> k) & ((Word)1); + if (x_p == 1 && x_n == 1) + { + fprintf(stderr, "error\n"); + } + + + + if (x_p == 1) + { + diff = 1; + ///fprintf(stderr, "[+1]"); + } + + if (x_n == 1) + { + diff = -1; + ///fprintf(stderr, "[-1]"); + } + + if (x_p == 0 && x_n == 0) + { + diff = 0; + ///fprintf(stderr, "[+0]"); + } + j++; + if (matrix[i][j] - matrix[i][j - 1] != diff) + { + fprintf(stderr, "*************V(k): %u\n", k); + ///return; + } + k++; + } + ///fprintf(stderr, "\n"); +} + + + + +inline void prase_D0(Word D0, int length, int matrix[1000][1000], int i) +{ + + i++; + int k = 0; + int j = i; + int diff; + while (k < length) + { + + int diff = (D0 >> k) & ((Word)1); + + if(diff==matrix[i][j] - matrix[i-1][j-1]) + { + fprintf(stderr, "*************D(k): %u\n", k); + fprintf(stderr, "diff: %u, matrix[i][j]: %u, matrix[i-1][j-1]: %u\n", diff, matrix[i][j], matrix[i-1][j-1]); + ///return; + } + j++; + k++; + } + ///fprintf(stderr, "\n"); +} + +inline void prase_H(Word HP, Word HN, int length, int matrix[1000][1000], int i) +{ + + i++; + int k = 0; + int j = i; + int diff; + while (k < length) + { + + int x_p = (HP >> k) & ((Word)1); + int x_n = (HN >> k) & ((Word)1); + if (x_p == 1 && x_n == 1) + { + fprintf(stderr, "error\n"); + } + if (x_p == 1) + { + diff = 1; + ///fprintf(stderr, "[+1]"); + } + if (x_n == 1) + { + diff = -1; + ///fprintf(stderr, "[-1]"); + } + if (x_p == 0 && x_n == 0) + { + diff = 0; + ///fprintf(stderr, "[+0]"); + } + + if(diff!=matrix[i][j] - matrix[i-1][j]) + { + fprintf(stderr, "*************H(k): %u\n", k); + ///return; + } + + j++; + k++; + } + ///fprintf(stderr, "\n"); +} + + +/** + pattern是长的那个，是y + p_length是长的那个的长度, p_length实际没用 + text是短的那个，是x + t_length是短的那个的长度 + errthold是阈值 + return_err是编辑距离 + 返回值是结束位置 + **/ +inline int Reserve_Banded_BPM_debug +(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err, int matrix[1000][1000]) +{ + (*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= 6) + if (i >= 0) + { + /** + fprintf(stderr, "VP:\n"); + output_bit_myers(VP, band_length); + + fprintf(stderr, "VN:\n"); + output_bit_myers(VN, band_length); + + fprintf(stderr, "HP:\n"); + output_bit_myers(HP, band_length); + + fprintf(stderr, "HN:\n"); + output_bit_myers(HN, band_length); + + fprintf(stderr, "D0:\n"); + output_bit_myers(D0, band_length); + + fprintf(stderr, "text[i]: %c\n", text[i]); + + fprintf(stderr, "Peq[text[i]]:\n"); + output_bit_myers(Peq[text[i]], band_length); + + for (size_t j = i; j < i + band_length; j++) + { + fprintf(stderr, "%c", pattern[j]); + } + fprintf(stderr, "\n"); + + fprintf(stderr, "Previous begin.\n"); + prase_vertical(VP, VN, band_length, matrix, i-1); + prase_D0(D0, band_length, matrix, i-1); + prase_H(HP, HN, band_length, matrix, i-1); + fprintf(stderr, "Previous test done.\n"); + **/ + + X = Peq[text[i]] | VN; + /** + fprintf(stderr, "#X:\n"); + output_bit_myers(X, band_length); + **/ + + D0 = ((VP + (X&VP)) ^ VP) | X; + /** + fprintf(stderr, "#(X&VP):\n"); + output_bit_myers((X&VP), band_length); + + fprintf(stderr, "#(VP + (X&VP)):\n"); + output_bit_myers((VP + (X&VP)), band_length); + + fprintf(stderr, "#((VP + (X&VP)) ^ VP):\n"); + output_bit_myers(((VP + (X&VP)) ^ VP), band_length); + + fprintf(stderr, "#D0:\n"); + output_bit_myers(D0, band_length); + **/ + HN = VP&D0; + HP = VN | ~(VP | D0); + + X = D0 >> 1; + VN = X&HP; + VP = HN | ~(X | HP); + } + else + { + ///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0] + 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); + } + + + + + + + + /** + for (size_t j = i + 1; j <= i + band_length + 1; j++) + { + fprintf(stderr, "[%u]", matrix[i + 1][j]); + } + fprintf(stderr, "\n"); + **/ + + prase_vertical(VP, VN, band_length, matrix, i); + prase_D0(D0, band_length, matrix, i); + prase_H(HP, HN, band_length, matrix, i); + /** + fprintf(stderr, "VP:\n"); + output_bit_myers(VP, band_length); + + fprintf(stderr, "VN:\n"); + output_bit_myers(VN, band_length); + **/ + + + 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']; + } + + + + + + 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; + } + + prase_vertical(VP, VN, band_length, matrix, i); + prase_D0(D0, band_length, matrix, i); + prase_H(HP, HN, band_length, matrix, i); + + + fprintf(stderr, "err: %d, matrix[][]: %d\n", err, matrix[i+1][i+1]); + fprintf(stderr, "VP:\n"); + output_bit_myers(VP, band_length); + + fprintf(stderr, "VN:\n"); + output_bit_myers(VN, band_length); + + ////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; + + fprintf(stderr, "*i: %u, err: %d\n", i, err); + + 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; + + fprintf(stderr, "*i: %u, err: %d\n", i, err); + + 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; + } + + return return_site; + +} + + +/** + pattern是长的那个，是y + p_length是长的那个的长度, p_length实际没用 + text是短的那个，是x + t_length是短的那个的长度 + errthold是阈值 + return_err是编辑距离 + 返回值是结束位置 + **/ +inline int Reserve_Banded_BPM +(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err) +{ + (*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']; + } + + + + + + 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; + } + + + ////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; + } + + return return_site; + +} + + + + + + + + +inline int BS_Reserve_Banded_BPM +(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, unsigned int* return_err) +{ + (*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 + ///ÕâÊÇ°ÑpatternµÄÇ°2k + 1¸ö×Ö·ûÔ¤´¦Àí + ///pattern[0]¶ÔÓ¦Peq[0] + ///pattern[2k]¶ÔÓ¦Peq[2k] + for (i = 0; i> 1; + VN = X&HP; + VP = HN | ~(X | HP); + ///Èç¹ûб¶Ô½ÇÏß·½ÏòÆ¥ÅäÔòD0ÊÇ1 + ///Èç¹û²»Æ¥ÅäÔòD0ÊÇ0 + ///Õâ¸öÒâ˼ÊÇÈç¹û×îÉÏÃæÄÇÌõ¶Ô½ÇÏßÉϵÄб¶Ô½ÇÏß·½Ïò·¢ÉúÎóÅä,ÔòÖ´ÐÐÄÚ²¿³ÌÐò + /// + if (!(D0&err_mask)) + { + ++err; + + ///¼´Ê¹È«²¿µÝ¼õ£¬Ò²¾Í¼õ2k + if ((err - last_high)>errthold) + return -1; + } + + ///pattern[0]ÔÚPeq[2k], ¶øpattern[2k]ÔÚPeq[0] + //ÓÒÒÆʵ¼ÊÉÏÊÇ°Ñpattern[0]ÒƵôÁË + Peq['A'] = Peq['A'] >> 1; + Peq['C'] = Peq['C'] >> 1; + Peq['G'] = Peq['G'] >> 1; + Peq['T'] = Peq['T'] >> 1; + + + ++i; + ++i_bd; + ///ÕâÊÇ°ÑеÄpattern[2k]¼Ó½øÀ´, ÕâòËÆÊǼӵ½Peq[2k]ÉÏÁË + Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask; + + + Peq['T'] = Peq['T'] | Peq['C']; + } + + + + + ///fprintf(stderr, "sucess(1)\n"); + + + ///Õâ¸öÑ­»·ÄóöÀ´ÊÇΪÁË·ÀÖ¹ÄÚ´æй¶ + ///ÆäʵҲ¾ÍÊÇÑ­»·ÀïµÄ×îºóÒ»ÐÐÓï¾ä°É + ///ÍêÈ«¿ÉÒÔ°ÑpatternÔö´óһλ + ///²»¹ýÕâÑùÒ²ºÃ£¬¿ÉÒÔ¼õÉÙ¼ÆË㿪Ïú + 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; + } + + + + + + ////fprintf(stderr, "sucess(2)\n"); + + /// last_high = 2k + /// site = (SEQ_LENGTH + 2k) - 2k -1 + /// site = SEQ_LENGTH - 1 + ///´ËʱÕâ¸ösiteòËÆÊÇ×îÉÏÃæÄÇÌõ¶Ô½ÇÏßµÄλÖà + ///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; + } + + return return_site; + +} + + + +inline int Reserve_Banded_BPM_new(char *pattern,int p_length,char *text,int t_length, + unsigned short errthold,unsigned short band_down,unsigned short band_below,unsigned short band_length,int* return_err, int thread_id) +{ + + Word Peq[128]; + char Peq_index[4]= {'A','C','G','T'}; + int symbol = 0; + int r; + 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; + + for (r =0; r>1; + VN=X&HP; + VP=HN|~(X|HP); + if(!(D0&err_mask)) + { + ++err; + if((err-last_high)>errthold) + 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; + } + + + ///这个循环拿出来是为了防止内存泄露 + + 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; + } + + + int site=p_length-last_high-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; + } + + + } + return return_site; + +} + + +#endif \ No newline at end of file diff --git a/Makefile b/Makefile index fc0fe01..c4fa7f3 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ CC=g++ CFLAGS = -w -c -msse4.2 -mpopcnt -fomit-frame-pointer -Winline -O3 -lz LDFLAGS = -lm -lz -lpthread -O3 -mpopcnt -msse4.2 -lz -w -SOURCES = main.cpp CommandLines.cpp Process_Read.cpp Assembly.cpp kmer.cpp Hash_Table.cpp POA.cpp +SOURCES = main.cpp CommandLines.cpp Process_Read.cpp Assembly.cpp kmer.cpp Hash_Table.cpp POA.cpp Correct.cpp Levenshtein_distance.cpp edlib.cpp OBJECTS = $(SOURCES:.c=.o) EXECUTABLE = ccs_assembly diff --git a/Process_Read.cpp b/Process_Read.cpp index d3d4e03..c95d19c 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -237,6 +237,111 @@ void init_UC_Read(UC_Read* r) } +void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID) +{ + + + + long long readLen = Get_READ_LENGTH((*R_INF), ID); + uint8_t* src = Get_READ((*R_INF), ID); + + long long i; + long long copyLen; + long long end_pos = start_pos + length - 1; + + if (strand == 0) + { + + i = start_pos; + copyLen = 0; + + long long initLen = start_pos % 4; + + if (initLen != 0) + { + memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen); + copyLen = copyLen + 4 - initLen; + i = i + copyLen; + } + + + while (copyLen < length) + { + memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4); + copyLen = copyLen + 4; + i = i + 4; + } + + + if (R_INF->N_site[ID]) + { + for (i = 1; i <= R_INF->N_site[ID][0]; i++) + { + if (R_INF->N_site[ID][i] >= start_pos && R_INF->N_site[ID][i] <= end_pos) + { + r[R_INF->N_site[ID][i] - start_pos] = 'N'; + } + else if(R_INF->N_site[ID][i] > end_pos) + { + break; + } + } + } + + + } + else + { + + start_pos = readLen - start_pos - 1; + end_pos = readLen - end_pos - 1; + + + + ///start_pos > end_pos + i = start_pos; + copyLen = 0; + long long initLen = (start_pos + 1) % 4; + + if (initLen != 0) + { + memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen); + copyLen = copyLen + initLen; + i = i - initLen; + } + + while (copyLen < length) + { + memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4); + copyLen = copyLen + 4; + i = i - 4; + } + + if (R_INF->N_site[ID]) + { + long long offset = readLen - start_pos - 1; + + for (i = 1; i <= R_INF->N_site[ID][0]; i++) + { + + if (R_INF->N_site[ID][i] >= end_pos && R_INF->N_site[ID][i] <= start_pos) + { + r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N'; + } + else if(R_INF->N_site[ID][i] > start_pos) + { + break; + } + } + } + + + + } + +} + + void recover_UC_Read(UC_Read* r, All_reads* R_INF, uint64_t ID) { r->length = Get_READ_LENGTH((*R_INF), ID); diff --git a/Process_Read.h b/Process_Read.h index 9623601..ae8357a 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -111,6 +111,7 @@ void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_l void init_UC_Read(UC_Read* r); void recover_UC_Read(UC_Read* r, All_reads* R_INF, uint64_t ID); void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID); +void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID); void destory_UC_Read(UC_Read* r); void reverse_complement(char* pattern, uint64_t length); void write_All_reads(All_reads* r, char* read_file_name); diff --git a/edlib.cpp b/edlib.cpp new file mode 100644 index 0000000..1d75a84 --- /dev/null +++ b/edlib.cpp @@ -0,0 +1,1461 @@ +#include "edlib.h" + +#include +#include +#include +#include +#include +#include + +using namespace std; + +typedef uint64_t Word; +static const int WORD_SIZE = sizeof(Word) * 8; // Size of Word in bits +static const Word WORD_1 = (Word)1; +static const Word HIGH_BIT_MASK = WORD_1 << (WORD_SIZE - 1); // 100..00 +static const int MAX_UCHAR = 255; + +// Data needed to find alignment. +struct AlignmentData { + Word* Ps; + Word* Ms; + int* scores; + int* firstBlocks; + int* lastBlocks; + + AlignmentData(int maxNumBlocks, int targetLength) { + // We build a complete table and mark first and last block for each column + // (because algorithm is banded so only part of each columns is used). + // TODO: do not build a whole table, but just enough blocks for each column. + Ps = new Word[maxNumBlocks * targetLength]; + Ms = new Word[maxNumBlocks * targetLength]; + scores = new int[maxNumBlocks * targetLength]; + firstBlocks = new int[targetLength]; + lastBlocks = new int[targetLength]; + } + + ~AlignmentData() { + delete[] Ps; + delete[] Ms; + delete[] scores; + delete[] firstBlocks; + delete[] lastBlocks; + } +}; + +struct Block { + Word P; // Pvin + Word M; // Mvin + int score; // score of last cell in block; + + Block() {} + Block(Word P, Word M, int score) :P(P), M(M), score(score) {} +}; + + +/** + * Defines equality relation on alphabet characters. + * By default each character is always equal only to itself, but you can also provide additional equalities. + */ +class EqualityDefinition { +private: + bool matrix[MAX_UCHAR + 1][MAX_UCHAR + 1]; +public: + EqualityDefinition(const string& alphabet, + const EdlibEqualityPair* additionalEqualities = NULL, + const int additionalEqualitiesLength = 0) { + for (int i = 0; i < (int) alphabet.size(); i++) { + for (int j = 0; j < (int) alphabet.size(); j++) { + matrix[i][j] = (i == j); + } + } + if (additionalEqualities != NULL) { + for (int i = 0; i < additionalEqualitiesLength; i++) { + size_t firstTransformed = alphabet.find(additionalEqualities[i].first); + size_t secondTransformed = alphabet.find(additionalEqualities[i].second); + if (firstTransformed != string::npos && secondTransformed != string::npos) { + matrix[firstTransformed][secondTransformed] = matrix[secondTransformed][firstTransformed] = true; + } + } + } + } + + /** + * @param a Element from transformed sequence. + * @param b Element from transformed sequence. + * @return True if a and b are defined as equal, false otherwise. + */ + bool areEqual(unsigned char a, unsigned char b) const { + return matrix[a][b]; + } +}; + +static int myersCalcEditDistanceSemiGlobal(const Word* Peq, int W, int maxNumBlocks, + int queryLength, + const unsigned char* target, int targetLength, + int k, EdlibAlignMode mode, + int* bestScore_, int** positions_, int* numPositions_); + +static int myersCalcEditDistanceNW(const Word* Peq, int W, int maxNumBlocks, + int queryLength, + const unsigned char* target, int targetLength, + int k, int* bestScore_, + int* position_, bool findAlignment, + AlignmentData** alignData, int targetStopPosition); + + +static int obtainAlignment( + const unsigned char* query, const unsigned char* rQuery, int queryLength, + const unsigned char* target, const unsigned char* rTarget, int targetLength, + const EqualityDefinition& equalityDefinition, int alphabetLength, int bestScore, + unsigned char** alignment, int* alignmentLength); + +static int obtainAlignmentHirschberg( + const unsigned char* query, const unsigned char* rQuery, int queryLength, + const unsigned char* target, const unsigned char* rTarget, int targetLength, + const EqualityDefinition& equalityDefinition, int alphabetLength, int bestScore, + unsigned char** alignment, int* alignmentLength); + +static int obtainAlignmentTraceback(int queryLength, int targetLength, + int bestScore, const AlignmentData* alignData, + unsigned char** alignment, int* alignmentLength); + +static string transformSequences(const char* queryOriginal, int queryLength, + const char* targetOriginal, int targetLength, + unsigned char** queryTransformed, + unsigned char** targetTransformed); + +static inline int ceilDiv(int x, int y); + +static inline unsigned char* createReverseCopy(const unsigned char* seq, int length); + +static inline Word* buildPeq(const int alphabetLength, + const unsigned char* query, + const int queryLength, + const EqualityDefinition& equalityDefinition); + + +/** + * Main edlib method. + */ +extern "C" EdlibAlignResult edlibAlign(const char* const queryOriginal, const int queryLength, + const char* const targetOriginal, const int targetLength, + const EdlibAlignConfig config) { + EdlibAlignResult result; + result.status = EDLIB_STATUS_OK; + result.editDistance = -1; + result.endLocations = result.startLocations = NULL; + result.numLocations = 0; + result.alignment = NULL; + result.alignmentLength = 0; + result.alphabetLength = 0; + + /*------------ TRANSFORM SEQUENCES AND RECOGNIZE ALPHABET -----------*/ + unsigned char* query, * target; + string alphabet = transformSequences(queryOriginal, queryLength, targetOriginal, targetLength, + &query, &target); + result.alphabetLength = (int) alphabet.size(); + /*-------------------------------------------------------*/ + + /*--------------------- INITIALIZATION ------------------*/ + int maxNumBlocks = ceilDiv(queryLength, WORD_SIZE); // bmax in Myers + int W = maxNumBlocks * WORD_SIZE - queryLength; // number of redundant cells in last level blocks + EqualityDefinition equalityDefinition(alphabet, config.additionalEqualities, config.additionalEqualitiesLength); + Word* Peq = buildPeq((int) alphabet.size(), query, queryLength, equalityDefinition); + /*-------------------------------------------------------*/ + + /*------------------ MAIN CALCULATION -------------------*/ + // TODO: Store alignment data only after k is determined? That could make things faster. + int positionNW; // Used only when mode is NW. + AlignmentData* alignData = NULL; + bool dynamicK = false; + int k = config.k; + if (k < 0) { // If valid k is not given, auto-adjust k until solution is found. + dynamicK = true; + k = WORD_SIZE; // Gives better results than smaller k. + } + + do { + if (config.mode == EDLIB_MODE_HW || config.mode == EDLIB_MODE_SHW) { + myersCalcEditDistanceSemiGlobal(Peq, W, maxNumBlocks, + queryLength, target, targetLength, + k, config.mode, &(result.editDistance), + &(result.endLocations), &(result.numLocations)); + } else { // mode == EDLIB_MODE_NW + myersCalcEditDistanceNW(Peq, W, maxNumBlocks, + queryLength, target, targetLength, + k, &(result.editDistance), &positionNW, + false, &alignData, -1); + } + k *= 2; + } while(dynamicK && result.editDistance == -1); + + if (result.editDistance >= 0) { // If there is solution. + // If NW mode, set end location explicitly. + if (config.mode == EDLIB_MODE_NW) { + result.endLocations = (int *) malloc(sizeof(int) * 1); + result.endLocations[0] = targetLength - 1; + result.numLocations = 1; + } + + // Find starting locations. + if (config.task == EDLIB_TASK_LOC || config.task == EDLIB_TASK_PATH) { + result.startLocations = (int*) malloc(result.numLocations * sizeof(int)); + if (config.mode == EDLIB_MODE_HW) { // If HW, I need to calculate start locations. + const unsigned char* rTarget = createReverseCopy(target, targetLength); + const unsigned char* rQuery = createReverseCopy(query, queryLength); + // Peq for reversed query. + Word* rPeq = buildPeq((int) alphabet.size(), rQuery, queryLength, equalityDefinition); + for (int i = 0; i < result.numLocations; i++) { + int endLocation = result.endLocations[i]; + if (endLocation == -1) { + // NOTE: Sometimes one of optimal solutions is that query starts before target, like this: + // AAGG <- target + // CCTT <- query + // It will never be only optimal solution and it does not happen often, however it is + // possible and in that case end location will be -1. What should we do with that? + // Should we just skip reporting such end location, although it is a solution? + // If we do report it, what is the start location? -4? -1? Nothing? + // TODO: Figure this out. This has to do in general with how we think about start + // and end locations. + // Also, we have alignment later relying on this locations to limit the space of it's + // search -> how can it do it right if these locations are negative or incorrect? + result.startLocations[i] = 0; // I put 0 for now, but it does not make much sense. + } else { + int bestScoreSHW, numPositionsSHW; + int* positionsSHW; + myersCalcEditDistanceSemiGlobal( + rPeq, W, maxNumBlocks, + queryLength, rTarget + targetLength - endLocation - 1, endLocation + 1, + result.editDistance, EDLIB_MODE_SHW, + &bestScoreSHW, &positionsSHW, &numPositionsSHW); + // Taking last location as start ensures that alignment will not start with insertions + // if it can start with mismatches instead. + result.startLocations[i] = endLocation - positionsSHW[numPositionsSHW - 1]; + free(positionsSHW); + } + } + delete[] rTarget; + delete[] rQuery; + delete[] rPeq; + } else { // If mode is SHW or NW + for (int i = 0; i < result.numLocations; i++) { + result.startLocations[i] = 0; + } + } + } + + // Find alignment -> all comes down to finding alignment for NW. + // Currently we return alignment only for first pair of locations. + if (config.task == EDLIB_TASK_PATH) { + int alnStartLocation = result.startLocations[0]; + int alnEndLocation = result.endLocations[0]; + const unsigned char* alnTarget = target + alnStartLocation; + const int alnTargetLength = alnEndLocation - alnStartLocation + 1; + const unsigned char* rAlnTarget = createReverseCopy(alnTarget, alnTargetLength); + const unsigned char* rQuery = createReverseCopy(query, queryLength); + obtainAlignment(query, rQuery, queryLength, + alnTarget, rAlnTarget, alnTargetLength, + equalityDefinition, (int) alphabet.size(), result.editDistance, + &(result.alignment), &(result.alignmentLength)); + delete[] rAlnTarget; + delete[] rQuery; + } + } + /*-------------------------------------------------------*/ + + //--- Free memory ---// + delete[] Peq; + free(query); + free(target); + if (alignData) delete alignData; + //-------------------// + + return result; +} + +extern "C" char* edlibAlignmentToCigar(const unsigned char* const alignment, const int alignmentLength, + const EdlibCigarFormat cigarFormat) { + if (cigarFormat != EDLIB_CIGAR_EXTENDED && cigarFormat != EDLIB_CIGAR_STANDARD) { + return 0; + } + + // Maps move code from alignment to char in cigar. + // 0 1 2 3 + char moveCodeToChar[] = {'=', 'I', 'D', 'X'}; + if (cigarFormat == EDLIB_CIGAR_STANDARD) { + moveCodeToChar[0] = moveCodeToChar[3] = 'M'; + } + + vector* cigar = new vector(); + char lastMove = 0; // Char of last move. 0 if there was no previous move. + int numOfSameMoves = 0; + for (int i = 0; i <= alignmentLength; i++) { + // if new sequence of same moves started + if (i == alignmentLength || (moveCodeToChar[alignment[i]] != lastMove && lastMove != 0)) { + // Write number of moves to cigar string. + int numDigits = 0; + for (; numOfSameMoves; numOfSameMoves /= 10) { + cigar->push_back('0' + numOfSameMoves % 10); + numDigits++; + } + reverse(cigar->end() - numDigits, cigar->end()); + // Write code of move to cigar string. + cigar->push_back(lastMove); + // If not at the end, start new sequence of moves. + if (i < alignmentLength) { + // Check if alignment has valid values. + if (alignment[i] > 3) { + delete cigar; + return 0; + } + numOfSameMoves = 0; + } + } + if (i < alignmentLength) { + lastMove = moveCodeToChar[alignment[i]]; + numOfSameMoves++; + } + } + cigar->push_back(0); // Null character termination. + char* cigar_ = (char*) malloc(cigar->size() * sizeof(char)); + memcpy(cigar_, &(*cigar)[0], cigar->size() * sizeof(char)); + delete cigar; + + return cigar_; +} + +/** + * Build Peq table for given query and alphabet. + * Peq is table of dimensions alphabetLength+1 x maxNumBlocks. + * Bit i of Peq[s * maxNumBlocks + b] is 1 if i-th symbol from block b of query equals symbol s, otherwise it is 0. + * NOTICE: free returned array with delete[]! + */ +static inline Word* buildPeq(const int alphabetLength, + const unsigned char* const query, + const int queryLength, + const EqualityDefinition& equalityDefinition) { + int maxNumBlocks = ceilDiv(queryLength, WORD_SIZE); + // table of dimensions alphabetLength+1 x maxNumBlocks. Last symbol is wildcard. + Word* Peq = new Word[(alphabetLength + 1) * maxNumBlocks]; + + // Build Peq (1 is match, 0 is mismatch). NOTE: last column is wildcard(symbol that matches anything) with just 1s + for (unsigned char symbol = 0; symbol <= alphabetLength; symbol++) { + for (int b = 0; b < maxNumBlocks; b++) { + if (symbol < alphabetLength) { + Peq[symbol * maxNumBlocks + b] = 0; + for (int r = (b+1) * WORD_SIZE - 1; r >= b * WORD_SIZE; r--) { + Peq[symbol * maxNumBlocks + b] <<= 1; + // NOTE: We pretend like query is padded at the end with W wildcard symbols + if (r >= queryLength || equalityDefinition.areEqual(query[r], symbol)) + Peq[symbol * maxNumBlocks + b] += 1; + } + } else { // Last symbol is wildcard, so it is all 1s + Peq[symbol * maxNumBlocks + b] = (Word)-1; + } + } + } + + return Peq; +} + + +/** + * Returns new sequence that is reverse of given sequence. + * Free returned array with delete[]. + */ +static inline unsigned char* createReverseCopy(const unsigned char* const seq, const int length) { + unsigned char* rSeq = new unsigned char[length]; + for (int i = 0; i < length; i++) { + rSeq[i] = seq[length - i - 1]; + } + return rSeq; +} + +/** + * Corresponds to Advance_Block function from Myers. + * Calculates one word(block), which is part of a column. + * Highest bit of word (one most to the left) is most bottom cell of block from column. + * Pv[i] and Mv[i] define vin of cell[i]: vin = cell[i] - cell[i-1]. + * @param [in] Pv Bitset, Pv[i] == 1 if vin is +1, otherwise Pv[i] == 0. + * @param [in] Mv Bitset, Mv[i] == 1 if vin is -1, otherwise Mv[i] == 0. + * @param [in] Eq Bitset, Eq[i] == 1 if match, 0 if mismatch. + * @param [in] hin Will be +1, 0 or -1. + * @param [out] PvOut Bitset, PvOut[i] == 1 if vout is +1, otherwise PvOut[i] == 0. + * @param [out] MvOut Bitset, MvOut[i] == 1 if vout is -1, otherwise MvOut[i] == 0. + * @param [out] hout Will be +1, 0 or -1. + */ +static inline int calculateBlock(Word Pv, Word Mv, Word Eq, const int hin, + Word &PvOut, Word &MvOut) { + // hin can be 1, -1 or 0. + // 1 -> 00...01 + // 0 -> 00...00 + // -1 -> 11...11 (2-complement) + + Word hinIsNeg = (Word)(hin >> 2) & WORD_1; // 00...001 if hin is -1, 00...000 if 0 or 1 + + Word Xv = Eq | Mv; + // This is instruction below written using 'if': if (hin < 0) Eq |= (Word)1; + Eq |= hinIsNeg; + Word Xh = (((Eq & Pv) + Pv) ^ Pv) | Eq; + + Word Ph = Mv | ~(Xh | Pv); + Word Mh = Pv & Xh; + + int hout = 0; + // This is instruction below written using 'if': if (Ph & HIGH_BIT_MASK) hout = 1; + hout = (Ph & HIGH_BIT_MASK) >> (WORD_SIZE - 1); + // This is instruction below written using 'if': if (Mh & HIGH_BIT_MASK) hout = -1; + hout -= (Mh & HIGH_BIT_MASK) >> (WORD_SIZE - 1); + + Ph <<= 1; + Mh <<= 1; + + // This is instruction below written using 'if': if (hin < 0) Mh |= (Word)1; + Mh |= hinIsNeg; + // This is instruction below written using 'if': if (hin > 0) Ph |= (Word)1; + Ph |= (Word)((hin + 1) >> 1); + + PvOut = Mh | ~(Xv | Ph); + MvOut = Ph & Xv; + + return hout; +} + +/** + * Does ceiling division x / y. + * Note: x and y must be non-negative and x + y must not overflow. + */ +static inline int ceilDiv(const int x, const int y) { + return x % y ? x / y + 1 : x / y; +} + +static inline int min(const int x, const int y) { + return x < y ? x : y; +} + +static inline int max(const int x, const int y) { + return x > y ? x : y; +} + + +/** + * @param [in] block + * @return Values of cells in block, starting with bottom cell in block. + */ +static inline vector getBlockCellValues(const Block block) { + vector scores(WORD_SIZE); + int score = block.score; + Word mask = HIGH_BIT_MASK; + for (int i = 0; i < WORD_SIZE - 1; i++) { + scores[i] = score; + if (block.P & mask) score--; + if (block.M & mask) score++; + mask >>= 1; + } + scores[WORD_SIZE - 1] = score; + return scores; +} + +/** + * Writes values of cells in block into given array, starting with first/top cell. + * @param [in] block + * @param [out] dest Array into which cell values are written. Must have size of at least WORD_SIZE. + */ +static inline void readBlock(const Block block, int* const dest) { + int score = block.score; + Word mask = HIGH_BIT_MASK; + for (int i = 0; i < WORD_SIZE - 1; i++) { + dest[WORD_SIZE - 1 - i] = score; + if (block.P & mask) score--; + if (block.M & mask) score++; + mask >>= 1; + } + dest[0] = score; +} + +/** + * Writes values of cells in block into given array, starting with last/bottom cell. + * @param [in] block + * @param [out] dest Array into which cell values are written. Must have size of at least WORD_SIZE. + */ +static inline void readBlockReverse(const Block block, int* const dest) { + int score = block.score; + Word mask = HIGH_BIT_MASK; + for (int i = 0; i < WORD_SIZE - 1; i++) { + dest[i] = score; + if (block.P & mask) score--; + if (block.M & mask) score++; + mask >>= 1; + } + dest[WORD_SIZE - 1] = score; +} + +/** + * @param [in] block + * @param [in] k + * @return True if all cells in block have value larger than k, otherwise false. + */ +static inline bool allBlockCellsLarger(const Block block, const int k) { + vector scores = getBlockCellValues(block); + for (int i = 0; i < WORD_SIZE; i++) { + if (scores[i] <= k) return false; + } + return true; +} + + +/** + * Uses Myers' bit-vector algorithm to find edit distance for one of semi-global alignment methods. + * @param [in] Peq Query profile. + * @param [in] W Size of padding in last block. + * TODO: Calculate this directly from query, instead of passing it. + * @param [in] maxNumBlocks Number of blocks needed to cover the whole query. + * TODO: Calculate this directly from query, instead of passing it. + * @param [in] queryLength + * @param [in] target + * @param [in] targetLength + * @param [in] k + * @param [in] mode EDLIB_MODE_HW or EDLIB_MODE_SHW + * @param [out] bestScore_ Edit distance. + * @param [out] positions_ Array of 0-indexed positions in target at which best score was found. + Make sure to free this array with free(). + * @param [out] numPositions_ Number of positions in the positions_ array. + * @return Status. + */ +static int myersCalcEditDistanceSemiGlobal( + const Word* const Peq, const int W, const int maxNumBlocks, + const int queryLength, + const unsigned char* const target, const int targetLength, + int k, const EdlibAlignMode mode, + int* const bestScore_, int** const positions_, int* const numPositions_) { + *positions_ = NULL; + *numPositions_ = 0; + + // firstBlock is 0-based index of first block in Ukkonen band. + // lastBlock is 0-based index of last block in Ukkonen band. + int firstBlock = 0; + int lastBlock = min(ceilDiv(k + 1, WORD_SIZE), maxNumBlocks) - 1; // y in Myers + Block *bl; // Current block + + Block* blocks = new Block[maxNumBlocks]; + + // For HW, solution will never be larger then queryLength. + if (mode == EDLIB_MODE_HW) { + k = min(queryLength, k); + } + + // Each STRONG_REDUCE_NUM column is reduced in more expensive way. + // This gives speed up of about 2 times for small k. + const int STRONG_REDUCE_NUM = 2048; + + // Initialize P, M and score + bl = blocks; + for (int b = 0; b <= lastBlock; b++) { + bl->score = (b + 1) * WORD_SIZE; + bl->P = (Word)-1; // All 1s + bl->M = (Word)0; + bl++; + } + + int bestScore = -1; + vector positions; // TODO: Maybe put this on heap? + const int startHout = mode == EDLIB_MODE_HW ? 0 : 1; // If 0 then gap before query is not penalized; + const unsigned char* targetChar = target; + for (int c = 0; c < targetLength; c++) { // for each column + const Word* Peq_c = Peq + (*targetChar) * maxNumBlocks; + + //----------------------- Calculate column -------------------------// + int hout = startHout; + bl = blocks + firstBlock; + Peq_c += firstBlock; + for (int b = firstBlock; b <= lastBlock; b++) { + hout = calculateBlock(bl->P, bl->M, *Peq_c, hout, bl->P, bl->M); + bl->score += hout; + bl++; Peq_c++; + } + bl--; Peq_c--; + //------------------------------------------------------------------// + + //---------- Adjust number of blocks according to Ukkonen ----------// + if ((lastBlock < maxNumBlocks - 1) && (bl->score - hout <= k) // bl is pointing to last block + && ((*(Peq_c + 1) & WORD_1) || hout < 0)) { // Peq_c is pointing to last block + // If score of left block is not too big, calculate one more block + lastBlock++; bl++; Peq_c++; + bl->P = (Word)-1; // All 1s + bl->M = (Word)0; + bl->score = (bl - 1)->score - hout + WORD_SIZE + calculateBlock(bl->P, bl->M, *Peq_c, hout, bl->P, bl->M); + } else { + while (lastBlock >= firstBlock && bl->score >= k + WORD_SIZE) { + lastBlock--; bl--; Peq_c--; + } + } + + // Every some columns, do some expensive but also more efficient block reducing. + // This is important! + // + // Reduce the band by decreasing last block if possible. + if (c % STRONG_REDUCE_NUM == 0) { + while (lastBlock >= 0 && lastBlock >= firstBlock && allBlockCellsLarger(*bl, k)) { + lastBlock--; bl--; Peq_c--; + } + } + // For HW, even if all cells are > k, there still may be solution in next + // column because starting conditions at upper boundary are 0. + // That means that first block is always candidate for solution, + // and we can never end calculation before last column. + if (mode == EDLIB_MODE_HW && lastBlock == -1) { + lastBlock++; bl++; Peq_c++; + } + + // Reduce band by increasing first block if possible. Not applicable to HW. + if (mode != EDLIB_MODE_HW) { + while (firstBlock <= lastBlock && blocks[firstBlock].score >= k + WORD_SIZE) { + firstBlock++; + } + if (c % STRONG_REDUCE_NUM == 0) { // Do strong reduction every some blocks + while (firstBlock <= lastBlock && allBlockCellsLarger(blocks[firstBlock], k)) { + firstBlock++; + } + } + } + + // If band stops to exist finish + if (lastBlock < firstBlock) { + *bestScore_ = bestScore; + if (bestScore != -1) { + *positions_ = (int *) malloc(sizeof(int) * (int) positions.size()); + *numPositions_ = (int) positions.size(); + copy(positions.begin(), positions.end(), *positions_); + } + delete[] blocks; + return EDLIB_STATUS_OK; + } + //------------------------------------------------------------------// + + //------------------------- Update best score ----------------------// + if (lastBlock == maxNumBlocks - 1) { + int colScore = bl->score; + if (colScore <= k) { // Scores > k dont have correct values (so we cannot use them), but are certainly > k. + // NOTE: Score that I find in column c is actually score from column c-W + if (bestScore == -1 || colScore <= bestScore) { + if (colScore != bestScore) { + positions.clear(); + bestScore = colScore; + // Change k so we will look only for equal or better + // scores then the best found so far. + k = bestScore; + } + positions.push_back(c - W); + } + } + } + //------------------------------------------------------------------// + + targetChar++; + } + + + // Obtain results for last W columns from last column. + if (lastBlock == maxNumBlocks - 1) { + vector blockScores = getBlockCellValues(*bl); + for (int i = 0; i < W; i++) { + int colScore = blockScores[i + 1]; + if (colScore <= k && (bestScore == -1 || colScore <= bestScore)) { + if (colScore != bestScore) { + positions.clear(); + k = bestScore = colScore; + } + positions.push_back(targetLength - W + i); + } + } + } + + *bestScore_ = bestScore; + if (bestScore != -1) { + *positions_ = (int *) malloc(sizeof(int) * (int) positions.size()); + *numPositions_ = (int) positions.size(); + copy(positions.begin(), positions.end(), *positions_); + } + + delete[] blocks; + return EDLIB_STATUS_OK; +} + + +/** + * Uses Myers' bit-vector algorithm to find edit distance for global(NW) alignment method. + * @param [in] Peq Query profile. + * @param [in] W Size of padding in last block. + * TODO: Calculate this directly from query, instead of passing it. + * @param [in] maxNumBlocks Number of blocks needed to cover the whole query. + * TODO: Calculate this directly from query, instead of passing it. + * @param [in] queryLength + * @param [in] target + * @param [in] targetLength + * @param [in] k + * @param [out] bestScore_ Edit distance. + * @param [out] position_ 0-indexed position in target at which best score was found. + * @param [in] findAlignment If true, whole matrix is remembered and alignment data is returned. + * Quadratic amount of memory is consumed. + * @param [out] alignData Data needed for alignment traceback (for reconstruction of alignment). + * Set only if findAlignment is set to true, otherwise it is NULL. + * Make sure to free this array using delete[]. + * @param [out] targetStopPosition If set to -1, whole calculation is performed normally, as expected. + * If set to p, calculation is performed up to position p in target (inclusive) + * and column p is returned as the only column in alignData. + * @return Status. + */ +static int myersCalcEditDistanceNW(const Word* const Peq, const int W, const int maxNumBlocks, + const int queryLength, + const unsigned char* const target, const int targetLength, + int k, int* const bestScore_, + int* const position_, const bool findAlignment, + AlignmentData** const alignData, const int targetStopPosition) { + if (targetStopPosition > -1 && findAlignment) { + // They can not be both set at the same time! + return EDLIB_STATUS_ERROR; + } + + // Each STRONG_REDUCE_NUM column is reduced in more expensive way. + const int STRONG_REDUCE_NUM = 2048; // TODO: Choose this number dinamically (based on query and target lengths?), so it does not affect speed of computation + + if (k < abs(targetLength - queryLength)) { + *bestScore_ = *position_ = -1; + return EDLIB_STATUS_OK; + } + + k = min(k, max(queryLength, targetLength)); // Upper bound for k + + // firstBlock is 0-based index of first block in Ukkonen band. + // lastBlock is 0-based index of last block in Ukkonen band. + int firstBlock = 0; + // This is optimal now, by my formula. + int lastBlock = min(maxNumBlocks, ceilDiv(min(k, (k + queryLength - targetLength) / 2) + 1, WORD_SIZE)) - 1; + Block* bl; // Current block + + Block* blocks = new Block[maxNumBlocks]; + + // Initialize P, M and score + bl = blocks; + for (int b = 0; b <= lastBlock; b++) { + bl->score = (b + 1) * WORD_SIZE; + bl->P = (Word)-1; // All 1s + bl->M = (Word)0; + bl++; + } + + // If we want to find alignment, we have to store needed data. + if (findAlignment) + *alignData = new AlignmentData(maxNumBlocks, targetLength); + else if (targetStopPosition > -1) + *alignData = new AlignmentData(maxNumBlocks, 1); + else + *alignData = NULL; + + const unsigned char* targetChar = target; + for (int c = 0; c < targetLength; c++) { // for each column + const Word* Peq_c = Peq + *targetChar * maxNumBlocks; + + //----------------------- Calculate column -------------------------// + int hout = 1; + bl = blocks + firstBlock; + for (int b = firstBlock; b <= lastBlock; b++) { + hout = calculateBlock(bl->P, bl->M, Peq_c[b], hout, bl->P, bl->M); + bl->score += hout; + bl++; + } + bl--; + //------------------------------------------------------------------// + // bl now points to last block + + // Update k. I do it only on end of column because it would slow calculation too much otherwise. + // NOTICE: I add W when in last block because it is actually result from W cells to the left and W cells up. + k = min(k, bl->score + + max(targetLength - c - 1, queryLength - ((1 + lastBlock) * WORD_SIZE - 1) - 1) + + (lastBlock == maxNumBlocks - 1 ? W : 0)); + + //---------- Adjust number of blocks according to Ukkonen ----------// + //--- Adjust last block ---// + // If block is not beneath band, calculate next block. Only next because others are certainly beneath band. + if (lastBlock + 1 < maxNumBlocks + && !(//score[lastBlock] >= k + WORD_SIZE || // NOTICE: this condition could be satisfied if above block also! + ((lastBlock + 1) * WORD_SIZE - 1 + > k - bl->score + 2 * WORD_SIZE - 2 - targetLength + c + queryLength))) { + lastBlock++; bl++; + bl->P = (Word)-1; // All 1s + bl->M = (Word)0; + int newHout = calculateBlock(bl->P, bl->M, Peq_c[lastBlock], hout, bl->P, bl->M); + bl->score = (bl - 1)->score - hout + WORD_SIZE + newHout; + hout = newHout; + } + + // While block is out of band, move one block up. + // NOTE: Condition used here is more loose than the one from the article, since I simplified the max() part of it. + // I could consider adding that max part, for optimal performance. + while (lastBlock >= firstBlock + && (bl->score >= k + WORD_SIZE + || ((lastBlock + 1) * WORD_SIZE - 1 > + // TODO: Does not work if do not put +1! Why??? + k - bl->score + 2 * WORD_SIZE - 2 - targetLength + c + queryLength + 1))) { + lastBlock--; bl--; + } + //-------------------------// + + //--- Adjust first block ---// + // While outside of band, advance block + while (firstBlock <= lastBlock + && (blocks[firstBlock].score >= k + WORD_SIZE + || ((firstBlock + 1) * WORD_SIZE - 1 < + blocks[firstBlock].score - k - targetLength + queryLength + c))) { + firstBlock++; + } + //--------------------------/ + + + // TODO: consider if this part is useful, it does not seem to help much + if (c % STRONG_REDUCE_NUM == 0) { // Every some columns do more expensive but more efficient reduction + while (lastBlock >= firstBlock) { + // If all cells outside of band, remove block + vector scores = getBlockCellValues(*bl); + int numCells = lastBlock == maxNumBlocks - 1 ? WORD_SIZE - W : WORD_SIZE; + int r = lastBlock * WORD_SIZE + numCells - 1; + bool reduce = true; + for (int i = WORD_SIZE - numCells; i < WORD_SIZE; i++) { + // TODO: Does not work if do not put +1! Why??? + if (scores[i] <= k && r <= k - scores[i] - targetLength + c + queryLength + 1) { + reduce = false; + break; + } + r--; + } + if (!reduce) break; + lastBlock--; bl--; + } + + while (firstBlock <= lastBlock) { + // If all cells outside of band, remove block + vector scores = getBlockCellValues(blocks[firstBlock]); + int numCells = firstBlock == maxNumBlocks - 1 ? WORD_SIZE - W : WORD_SIZE; + int r = firstBlock * WORD_SIZE + numCells - 1; + bool reduce = true; + for (int i = WORD_SIZE - numCells; i < WORD_SIZE; i++) { + if (scores[i] <= k && r >= scores[i] - k - targetLength + c + queryLength) { + reduce = false; + break; + } + r--; + } + if (!reduce) break; + firstBlock++; + } + } + + + // If band stops to exist finish + if (lastBlock < firstBlock) { + *bestScore_ = *position_ = -1; + delete[] blocks; + return EDLIB_STATUS_OK; + } + //------------------------------------------------------------------// + + + //---- Save column so it can be used for reconstruction ----// + if (findAlignment && c < targetLength) { + bl = blocks + firstBlock; + for (int b = firstBlock; b <= lastBlock; b++) { + (*alignData)->Ps[maxNumBlocks * c + b] = bl->P; + (*alignData)->Ms[maxNumBlocks * c + b] = bl->M; + (*alignData)->scores[maxNumBlocks * c + b] = bl->score; + (*alignData)->firstBlocks[c] = firstBlock; + (*alignData)->lastBlocks[c] = lastBlock; + bl++; + } + } + //----------------------------------------------------------// + //---- If this is stop column, save it and finish ----// + if (c == targetStopPosition) { + for (int b = firstBlock; b <= lastBlock; b++) { + (*alignData)->Ps[b] = (blocks + b)->P; + (*alignData)->Ms[b] = (blocks + b)->M; + (*alignData)->scores[b] = (blocks + b)->score; + (*alignData)->firstBlocks[0] = firstBlock; + (*alignData)->lastBlocks[0] = lastBlock; + } + *bestScore_ = -1; + *position_ = targetStopPosition; + delete[] blocks; + return EDLIB_STATUS_OK; + } + //----------------------------------------------------// + + targetChar++; + } + + if (lastBlock == maxNumBlocks - 1) { // If last block of last column was calculated + // Obtain best score from block -> it is complicated because query is padded with W cells + int bestScore = getBlockCellValues(blocks[lastBlock])[W]; + if (bestScore <= k) { + *bestScore_ = bestScore; + *position_ = targetLength - 1; + delete[] blocks; + return EDLIB_STATUS_OK; + } + } + + *bestScore_ = *position_ = -1; + delete[] blocks; + return EDLIB_STATUS_OK; +} + + +/** + * Finds one possible alignment that gives optimal score by moving back through the dynamic programming matrix, + * that is stored in alignData. Consumes large amount of memory: O(queryLength * targetLength). + * @param [in] queryLength Normal length, without W. + * @param [in] targetLength Normal length, without W. + * @param [in] bestScore Best score. + * @param [in] alignData Data obtained during finding best score that is useful for finding alignment. + * @param [out] alignment Alignment. + * @param [out] alignmentLength Length of alignment. + * @return Status code. + */ +static int obtainAlignmentTraceback(const int queryLength, const int targetLength, + const int bestScore, const AlignmentData* const alignData, + unsigned char** const alignment, int* const alignmentLength) { + const int maxNumBlocks = ceilDiv(queryLength, WORD_SIZE); + const int W = maxNumBlocks * WORD_SIZE - queryLength; + + *alignment = (unsigned char*) malloc((queryLength + targetLength - 1) * sizeof(unsigned char)); + *alignmentLength = 0; + int c = targetLength - 1; // index of column + int b = maxNumBlocks - 1; // index of block in column + int currScore = bestScore; // Score of current cell + int lScore = -1; // Score of left cell + int uScore = -1; // Score of upper cell + int ulScore = -1; // Score of upper left cell + Word currP = alignData->Ps[c * maxNumBlocks + b]; // P of current block + Word currM = alignData->Ms[c * maxNumBlocks + b]; // M of current block + // True if block to left exists and is in band + bool thereIsLeftBlock = c > 0 && b >= alignData->firstBlocks[c-1] && b <= alignData->lastBlocks[c-1]; + // We set initial values of lP and lM to 0 only to avoid compiler warnings, they should not affect the + // calculation as both lP and lM should be initialized at some moment later (but compiler can not + // detect it since this initialization is guaranteed by "business" logic). + Word lP = 0, lM = 0; + if (thereIsLeftBlock) { + lP = alignData->Ps[(c - 1) * maxNumBlocks + b]; // P of block to the left + lM = alignData->Ms[(c - 1) * maxNumBlocks + b]; // M of block to the left + } + currP <<= W; + currM <<= W; + int blockPos = WORD_SIZE - W - 1; // 0 based index of current cell in blockPos + + // TODO(martin): refactor this whole piece of code. There are too many if-else statements, + // it is too easy for a bug to hide and to hard to effectively cover all the edge-cases. + // We need better separation of logic and responsibilities. + while (true) { + if (c == 0) { + thereIsLeftBlock = true; + lScore = b * WORD_SIZE + blockPos + 1; + ulScore = lScore - 1; + } + + // TODO: improvement: calculate only those cells that are needed, + // for example if I calculate upper cell and can move up, + // there is no need to calculate left and upper left cell + //---------- Calculate scores ---------// + if (lScore == -1 && thereIsLeftBlock) { + lScore = alignData->scores[(c - 1) * maxNumBlocks + b]; // score of block to the left + for (int i = 0; i < WORD_SIZE - blockPos - 1; i++) { + if (lP & HIGH_BIT_MASK) lScore--; + if (lM & HIGH_BIT_MASK) lScore++; + lP <<= 1; + lM <<= 1; + } + } + if (ulScore == -1) { + if (lScore != -1) { + ulScore = lScore; + if (lP & HIGH_BIT_MASK) ulScore--; + if (lM & HIGH_BIT_MASK) ulScore++; + } + else if (c > 0 && b-1 >= alignData->firstBlocks[c-1] && b-1 <= alignData->lastBlocks[c-1]) { + // This is the case when upper left cell is last cell in block, + // and block to left is not in band so lScore is -1. + ulScore = alignData->scores[(c - 1) * maxNumBlocks + b - 1]; + } + } + if (uScore == -1) { + uScore = currScore; + if (currP & HIGH_BIT_MASK) uScore--; + if (currM & HIGH_BIT_MASK) uScore++; + currP <<= 1; + currM <<= 1; + } + //-------------------------------------// + + // TODO: should I check if there is upper block? + + //-------------- Move --------------// + // Move up - insertion to target - deletion from query + if (uScore != -1 && uScore + 1 == currScore) { + currScore = uScore; + lScore = ulScore; + uScore = ulScore = -1; + if (blockPos == 0) { // If entering new (upper) block + if (b == 0) { // If there are no cells above (only boundary cells) + (*alignment)[(*alignmentLength)++] = EDLIB_EDOP_INSERT; // Move up + for (int i = 0; i < c + 1; i++) // Move left until end + (*alignment)[(*alignmentLength)++] = EDLIB_EDOP_DELETE; + break; + } else { + blockPos = WORD_SIZE - 1; + b--; + currP = alignData->Ps[c * maxNumBlocks + b]; + currM = alignData->Ms[c * maxNumBlocks + b]; + if (c > 0 && b >= alignData->firstBlocks[c-1] && b <= alignData->lastBlocks[c-1]) { + thereIsLeftBlock = true; + lP = alignData->Ps[(c - 1) * maxNumBlocks + b]; // TODO: improve this, too many operations + lM = alignData->Ms[(c - 1) * maxNumBlocks + b]; + } else { + thereIsLeftBlock = false; + // TODO(martin): There may not be left block, but there can be left boundary - do we + // handle this correctly then? Are l and ul score set correctly? I should check that / refactor this. + } + } + } else { + blockPos--; + lP <<= 1; + lM <<= 1; + } + // Mark move + (*alignment)[(*alignmentLength)++] = EDLIB_EDOP_INSERT; + } + // Move left - deletion from target - insertion to query + else if (lScore != -1 && lScore + 1 == currScore) { + currScore = lScore; + uScore = ulScore; + lScore = ulScore = -1; + c--; + if (c == -1) { // If there are no cells to the left (only boundary cells) + (*alignment)[(*alignmentLength)++] = EDLIB_EDOP_DELETE; // Move left + int numUp = b * WORD_SIZE + blockPos + 1; + for (int i = 0; i < numUp; i++) // Move up until end + (*alignment)[(*alignmentLength)++] = EDLIB_EDOP_INSERT; + break; + } + currP = lP; + currM = lM; + if (c > 0 && b >= alignData->firstBlocks[c-1] && b <= alignData->lastBlocks[c-1]) { + thereIsLeftBlock = true; + lP = alignData->Ps[(c - 1) * maxNumBlocks + b]; + lM = alignData->Ms[(c - 1) * maxNumBlocks + b]; + } else { + if (c == 0) { // If there are no cells to the left (only boundary cells) + thereIsLeftBlock = true; + lScore = b * WORD_SIZE + blockPos + 1; + ulScore = lScore - 1; + } else { + thereIsLeftBlock = false; + } + } + // Mark move + (*alignment)[(*alignmentLength)++] = EDLIB_EDOP_DELETE; + } + // Move up left - (mis)match + else if (ulScore != -1) { + unsigned char moveCode = ulScore == currScore ? EDLIB_EDOP_MATCH : EDLIB_EDOP_MISMATCH; + currScore = ulScore; + uScore = lScore = ulScore = -1; + c--; + if (c == -1) { // If there are no cells to the left (only boundary cells) + (*alignment)[(*alignmentLength)++] = moveCode; // Move left + int numUp = b * WORD_SIZE + blockPos; + for (int i = 0; i < numUp; i++) // Move up until end + (*alignment)[(*alignmentLength)++] = EDLIB_EDOP_INSERT; + break; + } + if (blockPos == 0) { // If entering upper left block + if (b == 0) { // If there are no more cells above (only boundary cells) + (*alignment)[(*alignmentLength)++] = moveCode; // Move up left + for (int i = 0; i < c + 1; i++) // Move left until end + (*alignment)[(*alignmentLength)++] = EDLIB_EDOP_DELETE; + break; + } + blockPos = WORD_SIZE - 1; + b--; + currP = alignData->Ps[c * maxNumBlocks + b]; + currM = alignData->Ms[c * maxNumBlocks + b]; + } else { // If entering left block + blockPos--; + currP = lP; + currM = lM; + currP <<= 1; + currM <<= 1; + } + // Set new left block + if (c > 0 && b >= alignData->firstBlocks[c-1] && b <= alignData->lastBlocks[c-1]) { + thereIsLeftBlock = true; + lP = alignData->Ps[(c - 1) * maxNumBlocks + b]; + lM = alignData->Ms[(c - 1) * maxNumBlocks + b]; + } else { + if (c == 0) { // If there are no cells to the left (only boundary cells) + thereIsLeftBlock = true; + lScore = b * WORD_SIZE + blockPos + 1; + ulScore = lScore - 1; + } else { + thereIsLeftBlock = false; + } + } + // Mark move + (*alignment)[(*alignmentLength)++] = moveCode; + } else { + // Reached end - finished! + break; + } + //----------------------------------// + } + + *alignment = (unsigned char*) realloc(*alignment, (*alignmentLength) * sizeof(unsigned char)); + reverse(*alignment, *alignment + (*alignmentLength)); + return EDLIB_STATUS_OK; +} + + +/** + * Finds one possible alignment that gives optimal score (bestScore). + * It will split problem into smaller problems using Hirschberg's algorithm and when they are small enough, + * it will solve them using traceback algorithm. + * @param [in] query + * @param [in] rQuery Reversed query. + * @param [in] queryLength + * @param [in] target + * @param [in] rTarget Reversed target. + * @param [in] targetLength + * @param [in] equalityDefinition + * @param [in] alphabetLength + * @param [in] bestScore Best(optimal) score. + * @param [out] alignment Sequence of edit operations that make target equal to query. + * @param [out] alignmentLength Length of alignment. + * @return Status code. + */ +static int obtainAlignment( + const unsigned char* const query, const unsigned char* const rQuery, const int queryLength, + const unsigned char* const target, const unsigned char* const rTarget, const int targetLength, + const EqualityDefinition& equalityDefinition, const int alphabetLength, const int bestScore, + unsigned char** const alignment, int* const alignmentLength) { + + // Handle special case when one of sequences has length of 0. + if (queryLength == 0 || targetLength == 0) { + *alignmentLength = targetLength + queryLength; + *alignment = (unsigned char*) malloc((*alignmentLength) * sizeof(unsigned char)); + for (int i = 0; i < *alignmentLength; i++) { + (*alignment)[i] = queryLength == 0 ? EDLIB_EDOP_DELETE : EDLIB_EDOP_INSERT; + } + return EDLIB_STATUS_OK; + } + + const int maxNumBlocks = ceilDiv(queryLength, WORD_SIZE); + const int W = maxNumBlocks * WORD_SIZE - queryLength; + int statusCode; + + // TODO: think about reducing number of memory allocations in alignment functions, probably + // by sharing some memory that is allocated only once. That refers to: Peq, columns in Hirschberg, + // and it could also be done for alignments - we could have one big array for alignment that would be + // sparsely populated by each of steps in recursion, and at the end we would just consolidate those results. + + // If estimated memory consumption for traceback algorithm is smaller than 1MB use it, + // otherwise use Hirschberg's algorithm. By running few tests I choose boundary of 1MB as optimal. + long long alignmentDataSize = (long long) (2 * sizeof(Word) + sizeof(int)) * maxNumBlocks * targetLength + + (long long) 2 * sizeof(int) * targetLength; + if (alignmentDataSize < 1024 * 1024) { + int score_, endLocation_; // Used only to call function. + AlignmentData* alignData = NULL; + Word* Peq = buildPeq(alphabetLength, query, queryLength, equalityDefinition); + myersCalcEditDistanceNW(Peq, W, maxNumBlocks, + queryLength, + target, targetLength, + bestScore, + &score_, &endLocation_, true, &alignData, -1); + //assert(score_ == bestScore); + //assert(endLocation_ == targetLength - 1); + + statusCode = obtainAlignmentTraceback(queryLength, targetLength, + bestScore, alignData, alignment, alignmentLength); + delete alignData; + delete[] Peq; + } else { + statusCode = obtainAlignmentHirschberg(query, rQuery, queryLength, + target, rTarget, targetLength, + equalityDefinition, alphabetLength, bestScore, + alignment, alignmentLength); + } + return statusCode; +} + + +/** + * Finds one possible alignment that gives optimal score (bestScore). + * Uses Hirschberg's algorithm to split problem into two sub-problems, solve them and combine them together. + * @param [in] query + * @param [in] rQuery Reversed query. + * @param [in] queryLength + * @param [in] target + * @param [in] rTarget Reversed target. + * @param [in] targetLength + * @param [in] alphabetLength + * @param [in] bestScore Best(optimal) score. + * @param [out] alignment Sequence of edit operations that make target equal to query. + * @param [out] alignmentLength Length of alignment. + * @return Status code. + */ +static int obtainAlignmentHirschberg( + const unsigned char* const query, const unsigned char* const rQuery, const int queryLength, + const unsigned char* const target, const unsigned char* const rTarget, const int targetLength, + const EqualityDefinition& equalityDefinition, const int alphabetLength, const int bestScore, + unsigned char** const alignment, int* const alignmentLength) { + + const int maxNumBlocks = ceilDiv(queryLength, WORD_SIZE); + const int W = maxNumBlocks * WORD_SIZE - queryLength; + + Word* Peq = buildPeq(alphabetLength, query, queryLength, equalityDefinition); + Word* rPeq = buildPeq(alphabetLength, rQuery, queryLength, equalityDefinition); + + // Used only to call functions. + int score_, endLocation_; + + // Divide dynamic matrix into two halfs, left and right. + const int leftHalfWidth = targetLength / 2; + const int rightHalfWidth = targetLength - leftHalfWidth; + + // Calculate left half. + AlignmentData* alignDataLeftHalf = NULL; + int leftHalfCalcStatus = myersCalcEditDistanceNW( + Peq, W, maxNumBlocks, queryLength, target, targetLength, bestScore, + &score_, &endLocation_, false, &alignDataLeftHalf, leftHalfWidth - 1); + + // Calculate right half. + AlignmentData* alignDataRightHalf = NULL; + int rightHalfCalcStatus = myersCalcEditDistanceNW( + rPeq, W, maxNumBlocks, queryLength, rTarget, targetLength, bestScore, + &score_, &endLocation_, false, &alignDataRightHalf, rightHalfWidth - 1); + + delete[] Peq; + delete[] rPeq; + + if (leftHalfCalcStatus == EDLIB_STATUS_ERROR || rightHalfCalcStatus == EDLIB_STATUS_ERROR) { + if (alignDataLeftHalf) delete alignDataLeftHalf; + if (alignDataRightHalf) delete alignDataRightHalf; + return EDLIB_STATUS_ERROR; + } + + // Unwrap the left half. + int firstBlockIdxLeft = alignDataLeftHalf->firstBlocks[0]; + int lastBlockIdxLeft = alignDataLeftHalf->lastBlocks[0]; + // TODO: avoid this allocation by using some shared array? + // scoresLeft contains scores from left column, starting with scoresLeftStartIdx row (query index) + // and ending with scoresLeftEndIdx row (0-indexed). + int scoresLeftLength = (lastBlockIdxLeft - firstBlockIdxLeft + 1) * WORD_SIZE; + int* scoresLeft = new int[scoresLeftLength]; + for (int blockIdx = firstBlockIdxLeft; blockIdx <= lastBlockIdxLeft; blockIdx++) { + Block block(alignDataLeftHalf->Ps[blockIdx], alignDataLeftHalf->Ms[blockIdx], + alignDataLeftHalf->scores[blockIdx]); + readBlock(block, scoresLeft + (blockIdx - firstBlockIdxLeft) * WORD_SIZE); + } + int scoresLeftStartIdx = firstBlockIdxLeft * WORD_SIZE; + // If last block contains padding, shorten the length of scores for the length of padding. + if (lastBlockIdxLeft == maxNumBlocks - 1) { + scoresLeftLength -= W; + } + + // Unwrap the right half (I also reverse it while unwraping). + int firstBlockIdxRight = alignDataRightHalf->firstBlocks[0]; + int lastBlockIdxRight = alignDataRightHalf->lastBlocks[0]; + int scoresRightLength = (lastBlockIdxRight - firstBlockIdxRight + 1) * WORD_SIZE; + int* scoresRight = new int[scoresRightLength]; + int* scoresRightOriginalStart = scoresRight; + for (int blockIdx = firstBlockIdxRight; blockIdx <= lastBlockIdxRight; blockIdx++) { + Block block(alignDataRightHalf->Ps[blockIdx], alignDataRightHalf->Ms[blockIdx], + alignDataRightHalf->scores[blockIdx]); + readBlockReverse(block, scoresRight + (lastBlockIdxRight - blockIdx) * WORD_SIZE); + } + int scoresRightStartIdx = queryLength - (lastBlockIdxRight + 1) * WORD_SIZE; + // If there is padding at the beginning of scoresRight (that can happen because of reversing that we do), + // move pointer forward to remove the padding (that is why we remember originalStart). + if (scoresRightStartIdx < 0) { + //assert(scoresRightStartIdx == -1 * W); + scoresRight += W; + scoresRightStartIdx += W; + scoresRightLength -= W; + } + + delete alignDataLeftHalf; + delete alignDataRightHalf; + + //--------------------- Find the best move ----------------// + // Find the query/row index of cell in left column which together with its lower right neighbour + // from right column gives the best score (when summed). We also have to consider boundary cells + // (those cells at -1 indexes). + // x| + // -+- + // |x + int queryIdxLeftStart = max(scoresLeftStartIdx, scoresRightStartIdx - 1); + int queryIdxLeftEnd = min(scoresLeftStartIdx + scoresLeftLength - 1, + scoresRightStartIdx + scoresRightLength - 2); + int leftScore, rightScore; + int queryIdxLeftAlignment; // Query/row index of cell in left column where alignment is passing through. + bool queryIdxLeftAlignmentFound = false; + for (int queryIdx = queryIdxLeftStart; queryIdx <= queryIdxLeftEnd; queryIdx++) { + leftScore = scoresLeft[queryIdx - scoresLeftStartIdx]; + rightScore = scoresRight[queryIdx + 1 - scoresRightStartIdx]; + if (leftScore + rightScore == bestScore) { + queryIdxLeftAlignment = queryIdx; + queryIdxLeftAlignmentFound = true; + break; + } + } + // Check boundary cells. + if (!queryIdxLeftAlignmentFound && scoresLeftStartIdx == 0 && scoresRightStartIdx == 0) { + leftScore = leftHalfWidth; + rightScore = scoresRight[0]; + if (leftScore + rightScore == bestScore) { + queryIdxLeftAlignment = -1; + queryIdxLeftAlignmentFound = true; + } + } + if (!queryIdxLeftAlignmentFound && scoresLeftStartIdx + scoresLeftLength == queryLength + && scoresRightStartIdx + scoresRightLength == queryLength) { + leftScore = scoresLeft[scoresLeftLength - 1]; + rightScore = rightHalfWidth; + if (leftScore + rightScore == bestScore) { + queryIdxLeftAlignment = queryLength - 1; + queryIdxLeftAlignmentFound = true; + } + } + + delete[] scoresLeft; + delete[] scoresRightOriginalStart; + + if (queryIdxLeftAlignmentFound == false) { + // If there was no move that is part of optimal alignment, then there is no such alignment + // or given bestScore is not correct! + return EDLIB_STATUS_ERROR; + } + //----------------------------------------------------------// + + // Calculate alignments for upper half of left half (upper left - ul) + // and lower half of right half (lower right - lr). + const int ulHeight = queryIdxLeftAlignment + 1; + const int lrHeight = queryLength - ulHeight; + const int ulWidth = leftHalfWidth; + const int lrWidth = rightHalfWidth; + unsigned char* ulAlignment = NULL; int ulAlignmentLength; + int ulStatusCode = obtainAlignment(query, rQuery + lrHeight, ulHeight, + target, rTarget + lrWidth, ulWidth, + equalityDefinition, alphabetLength, leftScore, + &ulAlignment, &ulAlignmentLength); + unsigned char* lrAlignment = NULL; int lrAlignmentLength; + int lrStatusCode = obtainAlignment(query + ulHeight, rQuery, lrHeight, + target + ulWidth, rTarget, lrWidth, + equalityDefinition, alphabetLength, rightScore, + &lrAlignment, &lrAlignmentLength); + if (ulStatusCode == EDLIB_STATUS_ERROR || lrStatusCode == EDLIB_STATUS_ERROR) { + if (ulAlignment) free(ulAlignment); + if (lrAlignment) free(lrAlignment); + return EDLIB_STATUS_ERROR; + } + + // Build alignment by concatenating upper left alignment with lower right alignment. + *alignmentLength = ulAlignmentLength + lrAlignmentLength; + *alignment = (unsigned char*) malloc((*alignmentLength) * sizeof(unsigned char)); + memcpy(*alignment, ulAlignment, ulAlignmentLength); + memcpy(*alignment + ulAlignmentLength, lrAlignment, lrAlignmentLength); + + free(ulAlignment); + free(lrAlignment); + return EDLIB_STATUS_OK; +} + + +/** + * Takes char query and char target, recognizes alphabet and transforms them into unsigned char sequences + * where elements in sequences are not any more letters of alphabet, but their index in alphabet. + * Most of internal edlib functions expect such transformed sequences. + * This function will allocate queryTransformed and targetTransformed, so make sure to free them when done. + * Example: + * Original sequences: "ACT" and "CGT". + * Alphabet would be recognized as "ACTG". Alphabet length = 4. + * Transformed sequences: [0, 1, 2] and [1, 3, 2]. + * @param [in] queryOriginal + * @param [in] queryLength + * @param [in] targetOriginal + * @param [in] targetLength + * @param [out] queryTransformed It will contain values in range [0, alphabet length - 1]. + * @param [out] targetTransformed It will contain values in range [0, alphabet length - 1]. + * @return Alphabet as a string of unique characters, where index of each character is its value in transformed + * sequences. + */ +static string transformSequences(const char* const queryOriginal, const int queryLength, + const char* const targetOriginal, const int targetLength, + unsigned char** const queryTransformed, + unsigned char** const targetTransformed) { + // Alphabet is constructed from letters that are present in sequences. + // Each letter is assigned an ordinal number, starting from 0 up to alphabetLength - 1, + // and new query and target are created in which letters are replaced with their ordinal numbers. + // This query and target are used in all the calculations later. + *queryTransformed = (unsigned char *) malloc(sizeof(unsigned char) * queryLength); + *targetTransformed = (unsigned char *) malloc(sizeof(unsigned char) * targetLength); + + string alphabet = ""; + + // Alphabet information, it is constructed on fly while transforming sequences. + // letterIdx[c] is index of letter c in alphabet. + unsigned char letterIdx[MAX_UCHAR + 1]; + bool inAlphabet[MAX_UCHAR + 1]; // inAlphabet[c] is true if c is in alphabet + for (int i = 0; i < MAX_UCHAR + 1; i++) inAlphabet[i] = false; + + for (int i = 0; i < queryLength; i++) { + unsigned char c = static_cast(queryOriginal[i]); + if (!inAlphabet[c]) { + inAlphabet[c] = true; + letterIdx[c] = (unsigned char) alphabet.size(); + alphabet += queryOriginal[i]; + } + (*queryTransformed)[i] = letterIdx[c]; + } + for (int i = 0; i < targetLength; i++) { + unsigned char c = static_cast(targetOriginal[i]); + if (!inAlphabet[c]) { + inAlphabet[c] = true; + letterIdx[c] = (unsigned char) alphabet.size(); + alphabet += targetOriginal[i]; + } + (*targetTransformed)[i] = letterIdx[c]; + } + + return alphabet; +} + + +extern "C" EdlibAlignConfig edlibNewAlignConfig(int k, EdlibAlignMode mode, EdlibAlignTask task, + EdlibEqualityPair* additionalEqualities, + int additionalEqualitiesLength) { + EdlibAlignConfig config; + config.k = k; + config.mode = mode; + config.task = task; + config.additionalEqualities = additionalEqualities; + config.additionalEqualitiesLength = additionalEqualitiesLength; + return config; +} + +extern "C" EdlibAlignConfig edlibDefaultAlignConfig(void) { + return edlibNewAlignConfig(-1, EDLIB_MODE_NW, EDLIB_TASK_DISTANCE, NULL, 0); +} + +extern "C" void edlibFreeAlignResult(EdlibAlignResult result) { + if (result.endLocations) free(result.endLocations); + if (result.startLocations) free(result.startLocations); + if (result.alignment) free(result.alignment); +} diff --git a/edlib.h b/edlib.h new file mode 100644 index 0000000..4ed7c3b --- /dev/null +++ b/edlib.h @@ -0,0 +1,258 @@ +#ifndef EDLIB_H +#define EDLIB_H + +/** + * @file + * @author Martin Sosic + * @brief Main header file, containing all public functions and structures. + */ + +#ifdef __cplusplus +extern "C" { +#endif + +// Status codes +#define EDLIB_STATUS_OK 0 +#define EDLIB_STATUS_ERROR 1 + + /** + * Alignment methods - how should Edlib treat gaps before and after query? + */ + typedef enum { + /** + * Global method. This is the standard method. + * Useful when you want to find out how similar is first sequence to second sequence. + */ + EDLIB_MODE_NW, + /** + * Prefix method. Similar to global method, but with a small twist - gap at query end is not penalized. + * What that means is that deleting elements from the end of second sequence is "free"! + * For example, if we had "AACT" and "AACTGGC", edit distance would be 0, because removing "GGC" from the end + * of second sequence is "free" and does not count into total edit distance. This method is appropriate + * when you want to find out how well first sequence fits at the beginning of second sequence. + */ + EDLIB_MODE_SHW, + /** + * Infix method. Similar as prefix method, but with one more twist - gaps at query end and start are + * not penalized. What that means is that deleting elements from the start and end of second sequence is "free"! + * For example, if we had ACT and CGACTGAC, edit distance would be 0, because removing CG from the start + * and GAC from the end of second sequence is "free" and does not count into total edit distance. + * This method is appropriate when you want to find out how well first sequence fits at any part of + * second sequence. + * For example, if your second sequence was a long text and your first sequence was a sentence from that text, + * but slightly scrambled, you could use this method to discover how scrambled it is and where it fits in + * that text. In bioinformatics, this method is appropriate for aligning read to a sequence. + */ + EDLIB_MODE_HW + } EdlibAlignMode; + + /** + * Alignment tasks - what do you want Edlib to do? + */ + typedef enum { + EDLIB_TASK_DISTANCE, //!< Find edit distance and end locations. + EDLIB_TASK_LOC, //!< Find edit distance, end locations and start locations. + EDLIB_TASK_PATH //!< Find edit distance, end locations and start locations and alignment path. + } EdlibAlignTask; + + /** + * Describes cigar format. + * @see http://samtools.github.io/hts-specs/SAMv1.pdf + * @see http://drive5.com/usearch/manual/cigar.html + */ + typedef enum { + EDLIB_CIGAR_STANDARD, //!< Match: 'M', Insertion: 'I', Deletion: 'D', Mismatch: 'M'. + EDLIB_CIGAR_EXTENDED //!< Match: '=', Insertion: 'I', Deletion: 'D', Mismatch: 'X'. + } EdlibCigarFormat; + +// Edit operations. +#define EDLIB_EDOP_MATCH 0 //!< Match. +#define EDLIB_EDOP_INSERT 1 //!< Insertion to target = deletion from query. +#define EDLIB_EDOP_DELETE 2 //!< Deletion from target = insertion to query. +#define EDLIB_EDOP_MISMATCH 3 //!< Mismatch. + + /** + * @brief Defines two given characters as equal. + */ + typedef struct { + char first; + char second; + } EdlibEqualityPair; + + /** + * @brief Configuration object for edlibAlign() function. + */ + typedef struct { + /** + * Set k to non-negative value to tell edlib that edit distance is not larger than k. + * Smaller k can significantly improve speed of computation. + * If edit distance is larger than k, edlib will set edit distance to -1. + * Set k to negative value and edlib will internally auto-adjust k until score is found. + */ + int k; + + /** + * Alignment method. + * EDLIB_MODE_NW: global (Needleman-Wunsch) + * EDLIB_MODE_SHW: prefix. Gap after query is not penalized. + * EDLIB_MODE_HW: infix. Gaps before and after query are not penalized. + */ + EdlibAlignMode mode; + + /** + * Alignment task - tells Edlib what to calculate. Less to calculate, faster it is. + * EDLIB_TASK_DISTANCE - find edit distance and end locations of optimal alignment paths in target. + * EDLIB_TASK_LOC - find edit distance and start and end locations of optimal alignment paths in target. + * EDLIB_TASK_PATH - find edit distance, alignment path (and start and end locations of it in target). + */ + EdlibAlignTask task; + + /** + * List of pairs of characters, where each pair defines two characters as equal. + * This way you can extend edlib's definition of equality (which is that each character is equal only + * to itself). + * This can be useful if you have some wildcard characters that should match multiple other characters, + * or e.g. if you want edlib to be case insensitive. + * Can be set to NULL if there are none. + */ + EdlibEqualityPair* additionalEqualities; + + /** + * Number of additional equalities, which is non-negative number. + * 0 if there are none. + */ + int additionalEqualitiesLength; + } EdlibAlignConfig; + + /** + * Helper method for easy construction of configuration object. + * @return Configuration object filled with given parameters. + */ + EdlibAlignConfig edlibNewAlignConfig(int k, EdlibAlignMode mode, EdlibAlignTask task, + EdlibEqualityPair* additionalEqualities, + int additionalEqualitiesLength); + + /** + * @return Default configuration object, with following defaults: + * k = -1, mode = EDLIB_MODE_NW, task = EDLIB_TASK_DISTANCE, no additional equalities. + */ + EdlibAlignConfig edlibDefaultAlignConfig(void); + + + /** + * Container for results of alignment done by edlibAlign() function. + */ + typedef struct { + /** + * EDLIB_STATUS_OK or EDLIB_STATUS_ERROR. If error, all other fields will have undefined values. + */ + int status; + + /** + * -1 if k is non-negative and edit distance is larger than k. + */ + int editDistance; + + /** + * Array of zero-based positions in target where optimal alignment paths end. + * If gap after query is penalized, gap counts as part of query (NW), otherwise not. + * Set to NULL if edit distance is larger than k. + * If you do not free whole result object using edlibFreeAlignResult(), do not forget to use free(). + */ + int* endLocations; + + /** + * Array of zero-based positions in target where optimal alignment paths start, + * they correspond to endLocations. + * If gap before query is penalized, gap counts as part of query (NW), otherwise not. + * Set to NULL if not calculated or if edit distance is larger than k. + * If you do not free whole result object using edlibFreeAlignResult(), do not forget to use free(). + */ + int* startLocations; + + /** + * Number of end (and start) locations. + */ + int numLocations; + + /** + * Alignment is found for first pair of start and end locations. + * Set to NULL if not calculated. + * Alignment is sequence of numbers: 0, 1, 2, 3. + * 0 stands for match. + * 1 stands for insertion to target. + * 2 stands for insertion to query. + * 3 stands for mismatch. + * Alignment aligns query to target from begining of query till end of query. + * If gaps are not penalized, they are not in alignment. + * If you do not free whole result object using edlibFreeAlignResult(), do not forget to use free(). + */ + unsigned char* alignment; + + /** + * Length of alignment. + */ + int alignmentLength; + + /** + * Number of different characters in query and target together. + */ + int alphabetLength; + } EdlibAlignResult; + + /** + * Frees memory in EdlibAlignResult that was allocated by edlib. + * If you do not use it, make sure to free needed members manually using free(). + */ + void edlibFreeAlignResult(EdlibAlignResult result); + + + /** + * Aligns two sequences (query and target) using edit distance (levenshtein distance). + * Through config parameter, this function supports different alignment methods (global, prefix, infix), + * as well as different modes of search (tasks). + * It always returns edit distance and end locations of optimal alignment in target. + * It optionally returns start locations of optimal alignment in target and alignment path, + * if you choose appropriate tasks. + * @param [in] query First sequence. + * @param [in] queryLength Number of characters in first sequence. + * @param [in] target Second sequence. + * @param [in] targetLength Number of characters in second sequence. + * @param [in] config Additional alignment parameters, like alignment method and wanted results. + * @return Result of alignment, which can contain edit distance, start and end locations and alignment path. + * Make sure to clean up the object using edlibFreeAlignResult() or by manually freeing needed members. + */ + EdlibAlignResult edlibAlign(const char* query, int queryLength, + const char* target, int targetLength, + const EdlibAlignConfig config); + + + /** + * Builds cigar string from given alignment sequence. + * @param [in] alignment Alignment sequence. + * 0 stands for match. + * 1 stands for insertion to target. + * 2 stands for insertion to query. + * 3 stands for mismatch. + * @param [in] alignmentLength + * @param [in] cigarFormat Cigar will be returned in specified format. + * @return Cigar string. + * I stands for insertion. + * D stands for deletion. + * X stands for mismatch. (used only in extended format) + * = stands for match. (used only in extended format) + * M stands for (mis)match. (used only in standard format) + * String is null terminated. + * Needed memory is allocated and given pointer is set to it. + * Do not forget to free it later using free()! + */ + char* edlibAlignmentToCigar(const unsigned char* alignment, int alignmentLength, + EdlibCigarFormat cigarFormat); + + + +#ifdef __cplusplus +} +#endif + +#endif // EDLIB_H diff --git a/main.cpp b/main.cpp index 1b55d49..9f8e916 100644 --- a/main.cpp +++ b/main.cpp @@ -3,7 +3,8 @@ #include "CommandLines.h" #include "Process_Read.h" #include "Assembly.h" - +#include "Levenshtein_distance.h" +#include "edlib.h" /********************************for debug***************************************/ ///使用这个函数的时候，必须把Counting_multiple_thr()里的destory_Total_Count_Table(&TCB)注释掉 void debug_Counting() @@ -15,6 +16,182 @@ void debug_Counting() } +int matrix[1000][1000] = {0}; +///y_length > x_length +int edit_distance_normal(char* y, int y_length, char* x, int x_length) +{ memset(matrix, 0, sizeof(matrix)); + + int i, j; + for (i = 0; i <= x_length; i++) + { + matrix[i][0] = i; + } + + int digonal, up, left, min; + + ///一列列算的 + for (i = 0; i < x_length; i++) + { + for (j = 0; j < y_length; j++) + { + ///matrix[i + 1][j + 1] + digonal = matrix[i][j] + (x[i] != y[j]); + up = matrix[i + 1][j] + 1; + left = matrix[i][j + 1] + 1; + min = digonal; + if (up < min) + { + min = up; + } + + if (left< min) + { + min = left; + } + + matrix[i + 1][j + 1] = min; + } + } + + min = 999999; + for (j = 0; j <= y_length; j++) + { + if (matrix[i][j] < min) + { + min = matrix[i][j]; + } + } + + + return min; +} + + +///y_length > x_length +int edit_distance_normal_banded(char* y, int y_length, char* x, int x_length, int error) +{ memset(matrix, 0, sizeof(matrix)); + + int i, j; + for (i = 0; i <= x_length; i++) + { + for (j = 0; j <= y_length; j++) + { + matrix[i][j] = 1000000; + } + } + + for (i = 0; i <= x_length; i++) + { + matrix[i][0] = i; + } + + for (i = 0; i <= y_length; i++) + { + matrix[0][i] = 0; + } + + int banded_length = error*2 + 1; + + int digonal, up, left, min; + + + for (i = 0; i < x_length; i++) + { + ///for (j = 0; j < y_length; j++) + for (j = i; j < banded_length + i; j++) + { + ///matrix[i + 1][j + 1] + digonal = matrix[i][j] + (x[i] != y[j]); + up = matrix[i + 1][j] + 1; + left = matrix[i][j + 1] + 1; + min = digonal; + if (up < min) + { + min = up; + } + + if (left< min) + { + min = left; + } + + matrix[i + 1][j + 1] = min; + } + } + min = 999999; + for (j = 0; j <= y_length; j++) + ///for (j = i; j < banded_length + i; j++) + { + if (matrix[i][j] < min) + { + min = matrix[i][j]; + } + } + + + return min; +} + +void debug_edit_distance() +{ + /** + char* x = "TTCCATACGATTCCATTCAATTCGAGACCATTCTATTCCTGTCCATTCCTTGTGGTTCGATTCCATTTCACTCTAGTCCATTCCATTCCATTCAATTCCATTCGACTCTATTCCGTTCCACTCAATTCCATTCCATTCGATTCCATTTTTTTCGAGAACCTTCCATTACACTCCCTTCCATTCCAGTGCATTCCATTCCAGTCTCTTCAGTTCGATTCCATTCCATTCGTTTCGATTCCTTTCCATTCCAGCCCATTCCATTCCATTCCATTCCTTTCCTTTCCGTTTCATTAGATTCCATTGCATTCGATTCCATTCAAATCAATTCCGTTCTATTCAATTTGATTCAT"; + char* y = "CCATACGATTCCATTCAATTCGAGACCATTCTATTCCTGTCCATTCCTTGTGGTTCGATTCCATTTCACTCTAGTCCATTCCATTCCATTCAATTCCATTCGACTCTATTCCGTTCCATTCAATTCCATTCCATTCGATTCCATTTTTTTCGAGAACCTTCCATTACACTCCCTTCCATTCCAGTGCATTCCATTCCAGTCTCTTCACTTCGATTCCATTCCATTCGTTTCGATTCCTTTCCATTCCAGCCCATTCCATTCCATTCCATTCCTTTCCTTTCCGTTTCATTAGATTCCATTGCATTCCATTCCATTCAATTCAATTCCGTGCTATTCAATTTGATTCATTTCCATTTAATTCCATTCCATTAGATTCCATT"; + **/ + unsigned short toold = 7; + char* x = "TTCCATACGATTCCATTCAATTCGAGACCATTCTATTCCT"; + char* y = "CCATACGATTCCATTCAATTCGAGACCATTCTATTCCTGTCCATTCCTTGTGGT"; + fprintf(stderr, "x_length: %u\n", strlen(x)); + fprintf(stderr, "y_length: %u\n", strlen(y)); + + + EdlibAlignResult result = edlibAlign(x, strlen(x), y, strlen(y), + edlibNewAlignConfig(toold, EDLIB_MODE_HW, EDLIB_TASK_PATH, NULL, 0)); + + if (result.status == EDLIB_STATUS_OK) { + + fprintf(stderr, "****\nedlib: %d, alignmentLength: %d, startLocations: %d, endLocations: %d\n", + result.editDistance, result.alignmentLength, result.startLocations[0], result.endLocations[0]); + char* cigar = edlibAlignmentToCigar(result.alignment, result.alignmentLength, EDLIB_CIGAR_STANDARD); + fprintf(stderr,"%s\n", cigar); + free(cigar); + } + edlibFreeAlignResult(result); + + + + unsigned int error; + int end_site = Reserve_Banded_BPM(y, strlen(y), x, strlen(x), toold, &error); + + fprintf(stderr, "BPM: error: %u, end_site: %u\n", error, end_site); + + + unsigned short band_length=(toold+1)*3-1-1-toold; + unsigned short band_down=toold-1; + unsigned short band_blew=2*(toold+1)-1-1; + + + + int return_err = 99999; + ///注意pattern/text和band_down/band_blew是反的 + Reserve_Banded_BPM_new(y, strlen(y), x, strlen(x), + toold,band_blew,band_down,band_length, &return_err, 0); + + fprintf(stderr, "new BPM: error: %u\n", return_err); + + + return_err = edit_distance_normal(y, strlen(y), x, strlen(x)); + fprintf(stderr, "edit_distance_normal: error: %u\n", return_err); + + return_err = edit_distance_normal_banded(y, strlen(y), x, strlen(x), toold); + fprintf(stderr, "edit_distance_normal_banded: error: %u\n", return_err); + + end_site = Reserve_Banded_BPM_debug(y, strlen(y), x, strlen(x), toold, &error, matrix); + + fprintf(stderr, "BPM debug: error: %u, end_site: %u\n", error, end_site); +} + + int main(int argc, char *argv[]) { @@ -22,6 +199,13 @@ int main(int argc, char *argv[]) return 1; + /** + debug_edit_distance(); + + return 1; + **/ + + fprintf(stdout, "k-mer length: %d\n",k_mer_length); if (load_index_from_disk && load_pre_cauculated_index())