xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-10-10 03:30:56 +08:00
Code Issues Packages Projects Releases Wiki Activity

for break points at chr11

Browse Source
This commit is contained in:
Haoyu Cheng
2019-12-12 00:31:24 -05:00
parent 09603432e5
commit f1cde019a9
5 changed files with 1426 additions and 186 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+7 -73
View File
@@ -1919,8 +1919,9 @@ void* Overlap_calculate_heap_merge(void* arg)
for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num) for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
{ {
get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); ///get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5);
get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02);
clear_Cigar_record(&current_cigar); clear_Cigar_record(&current_cigar);
clear_Round2_alignment(&second_round); clear_Round2_alignment(&second_round);
@@ -2119,8 +2120,8 @@ void* Output_related_reads(void* arg)
&& &&
memcmp(required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) memcmp(required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
{ {
get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); ////get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5);
get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02);
fprintf(stderr, ">%.*s\n", Get_NAME_LENGTH((R_INF), i), fprintf(stderr, ">%.*s\n", Get_NAME_LENGTH((R_INF), i),
Get_NAME((R_INF), i)); Get_NAME((R_INF), i));
@@ -2998,7 +2999,8 @@ void* Final_overlap_calculate_heap_merge(void* arg)
get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); ////get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5);
get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.001);
/** /**
correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon, correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon,
&matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1, &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1,
@@ -3081,22 +3083,6 @@ void* Final_overlap_calculate_heap_merge(void* arg)
// if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i),
// Get_NAME_LENGTH((R_INF), i)) == 0)
// {
// fprintf(stderr, "\n2\n");
// debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "second");
// }
// if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i),
// Get_NAME_LENGTH((R_INF), i)) == 0)
// {
// fprintf(stderr, "\n2\n");
// debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "second");
// }
inner_j = 0; inner_j = 0;
j = 0; j = 0;
@@ -3125,22 +3111,6 @@ void* Final_overlap_calculate_heap_merge(void* arg)
} }
// if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i),
// Get_NAME_LENGTH((R_INF), i)) == 0)
// {
// fprintf(stderr, "\n3\n");
// debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "third");
// }
// if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i),
// Get_NAME_LENGTH((R_INF), i)) == 0)
// {
// fprintf(stderr, "\n3\n");
// debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "third");
// }
///recover missing exact overlaps ///recover missing exact overlaps
reverse_complement(g_read.seq, g_read.length); reverse_complement(g_read.seq, g_read.length);
for (j = 0; j < overlap_list.length; j++) for (j = 0; j < overlap_list.length; j++)
@@ -3171,23 +3141,6 @@ void* Final_overlap_calculate_heap_merge(void* arg)
} }
// if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i),
// Get_NAME_LENGTH((R_INF), i)) == 0)
// {
// fprintf(stderr, "\n4\n");
// debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "fourth");
// }
// if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i),
// Get_NAME_LENGTH((R_INF), i)) == 0)
// {
// fprintf(stderr, "\n4\n");
// debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "fourth");
// }
if(R_INF.paf[i].is_fully_corrected) if(R_INF.paf[i].is_fully_corrected)
{ {
potiental_matched_overlap_1++; potiental_matched_overlap_1++;
@@ -3201,25 +3154,6 @@ void* Final_overlap_calculate_heap_merge(void* arg)
push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf,
&overlap_list, &g_read, &overlap_read, 2, 0); &overlap_list, &g_read, &overlap_read, 2, 0);
// if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i),
// Get_NAME_LENGTH((R_INF), i)) == 0)
// {
// fprintf(stderr, "\n5\n");
// debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "last");
// debug_print_overlap(NULL, &overlap_list, &R_INF, "SET");
// }
// if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i),
// Get_NAME_LENGTH((R_INF), i)) == 0)
// {
// fprintf(stderr, "\n5\n");
// debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "last");
// debug_print_overlap(NULL, &overlap_list, &R_INF, "SET");
// }
} }
Correct.cpp
+649 -33
View File
@@ -537,19 +537,12 @@ int* r_extra_begin, int* r_extra_end, long long* r_y_start, long long* r_y_lengt
currentIDLen = Get_READ_LENGTH((*R_INF), y_ID); currentIDLen = Get_READ_LENGTH((*R_INF), y_ID);
///since Window_Len == x_len + (threshold << 1) ///since Window_Len == x_len + (threshold << 1)
if(currentIDLen <= y_start || if(y_start < 0 || currentIDLen <= y_start ||
currentIDLen - y_start + 2 * threshold + THRESHOLD_MAX_SIZE < Window_Len) currentIDLen - y_start + 2 * threshold + THRESHOLD_MAX_SIZE < Window_Len)
{ {
return 0; return 0;
} }
/**
if(currentIDLen <= y_start)
{
return 0;
}
**/
extra_begin = extra_end = 0; extra_begin = extra_end = 0;
///y maybe less than 0 ///y maybe less than 0
y_start = y_start - threshold; y_start = y_start - threshold;
@@ -1054,6 +1047,71 @@ UC_Read* g_read)
return T_error; return T_error;
} }
inline int double_error_threshold(int pre_threshold, int x_len)
{
pre_threshold = Adjust_Threshold(pre_threshold, x_len);
int threshold = pre_threshold * 2;
///may have some bugs
if(x_len >= 300 && threshold < THRESHOLD_MAX_SIZE)
{
threshold = THRESHOLD_MAX_SIZE;
}
if(threshold > THRESHOLD_MAX_SIZE)
{
threshold = THRESHOLD_MAX_SIZE;
}
return threshold;
}
inline int verify_sub_window(All_reads* R_INF, Correct_dumy* dumy, UC_Read* g_read,
long long x_beg, long long xLen, long long y_beg, long long yLen, uint64_t y_id,
uint64_t y_pos_strand, int threshold, int alignment_strand,
unsigned int* get_error, int* get_y_end, int* get_x_end, int* get_aligned_xLen)
{
(*get_aligned_xLen) = 0;
(*get_y_end) = -1;
(*get_x_end) = -1;
(*get_error) = (unsigned int)-1;
int extra_begin, extra_end, r_x_end, r_y_end, aligned_xLen;
long long o_len;
unsigned int r_error;
if(!determine_overlap_region(threshold, y_beg, y_id, yLen,
R_INF, &extra_begin, &extra_end, &y_beg, &o_len))
{
return 0;
}
fill_subregion(dumy->overlap_region, y_beg, o_len, y_pos_strand, R_INF,
y_id, extra_begin, extra_end);
char* x_string = g_read->seq + x_beg;
char* y_string = dumy->overlap_region;
aligned_xLen = 0;
alignment_extension(y_string, yLen, x_string, xLen, threshold,
alignment_strand, &r_error, &r_y_end, &r_x_end, &aligned_xLen);
(*get_error) = r_error;
(*get_y_end) = r_y_end;
(*get_x_end) = r_x_end;
(*get_aligned_xLen) = aligned_xLen;
if(aligned_xLen == 0)
{
return 0;
}
else
{
return 1;
}
}
inline double non_trim_error_rate(overlap_region_alloc* overlap_list, long long ID, inline double non_trim_error_rate(overlap_region_alloc* overlap_list, long long ID,
All_reads* R_INF, Correct_dumy* dumy, UC_Read* g_read) All_reads* R_INF, Correct_dumy* dumy, UC_Read* g_read)
{ {
@@ -1076,13 +1134,96 @@ All_reads* R_INF, Correct_dumy* dumy, UC_Read* g_read)
} }
else else
{ {
tError += int x_len = subWinLen;
get_high_error(overlap_list->list[ID].w_list[i].x_start, int threshold = double_error_threshold(overlap_list->list[ID].w_list[i].error_threshold, x_len);
overlap_list->list[ID].w_list[i].x_end, overlap_list->list[ID].w_list[i].y_start, int Window_Len = x_len + (threshold << 1);
overlap_list->list[ID].w_list[i].y_end, overlap_list->list[ID].y_id, unsigned int r_error_left;
overlap_list->list[ID].y_pos_strand, int r_x_end_left, r_y_end_left, aligned_xLen_left;
overlap_list->list[ID].w_list[i].error_threshold - overlap_list->list[ID].w_list[i].extra_begin, unsigned int r_error_right;
3, 0.3, R_INF, dumy, g_read); int r_x_end_right, r_y_end_right, aligned_xLen_right;
long long y_beg_left, y_beg_right;
aligned_xLen_left = aligned_xLen_right = 0;
y_beg_left = y_beg_right = -1;
if(overlap_list->list[ID].w_list[i].y_start == -1)
{
tError += x_len;
continue;
}
///if the previous window is mapped
if(i > 0 && overlap_list->list[ID].w_list[i - 1].y_end != -1)
{
y_beg_left = overlap_list->list[ID].w_list[i - 1].y_end + 1;
}
///if the next window is mapped
if(i < overlap_list->list[ID].w_list_length - 1 && overlap_list->list[ID].w_list[i + 1].y_end != -1)
{
y_beg_right = 1 + overlap_list->list[ID].w_list[i + 1].y_start - 1 - x_len;
}
///cannot apply this strategy to the leftmost and the rightmost window
if(y_beg_left == -1 && y_beg_right == -1 && i > 0 && i < subWinNum - 1)
{
y_beg_left = overlap_list->list[ID].w_list[i].y_start;
if(overlap_list->list[ID].w_list[i].extra_begin >= 0)
{
y_beg_left = y_beg_left + overlap_list->list[ID].w_list[i].error_threshold -
overlap_list->list[ID].w_list[i].extra_begin;
}
y_beg_right = y_beg_left;
}
if(y_beg_left != -1)
{
verify_sub_window(R_INF, dumy, g_read, overlap_list->list[ID].w_list[i].x_start,
x_len, y_beg_left, Window_Len, overlap_list->list[ID].y_id, overlap_list->list[ID].y_pos_strand,
threshold, 0, &r_error_left, &r_y_end_left, &r_x_end_left, &aligned_xLen_left);
}
if(y_beg_right != -1)
{
verify_sub_window(R_INF, dumy, g_read, overlap_list->list[ID].w_list[i].x_start,
x_len, y_beg_right, Window_Len, overlap_list->list[ID].y_id, overlap_list->list[ID].y_pos_strand,
threshold, 1, &r_error_right, &r_y_end_right, &r_x_end_right, &aligned_xLen_right);
}
///here one of aligned_xLen_left and aligned_xLen_right must be 0
if((i == 0 || i == subWinNum - 1) &&
(x_len - (aligned_xLen_left + aligned_xLen_right) > 10))
{
return 1.0;
}
///aligned in both direction
if(aligned_xLen_left != 0 && aligned_xLen_right != 0)
{
if(aligned_xLen_left + aligned_xLen_right <= x_len)
{
tError = tError + r_error_left + r_error_right +
(x_len - aligned_xLen_left - aligned_xLen_right);
}
else
{
float E_rate = (float)(x_len)/(float)(aligned_xLen_left + aligned_xLen_right);
tError = tError + (r_error_left + r_error_right)*E_rate;
}
}///not aligned in both direction
else if(aligned_xLen_left == 0 && aligned_xLen_right == 0)
{
tError += x_len;
}///only aligned in left
else if(aligned_xLen_left != 0)
{
tError = tError + r_error_left + (x_len - aligned_xLen_left);
}///only aligned in right
else if(aligned_xLen_right != 0)
{
tError = tError + r_error_right + (x_len - aligned_xLen_right);
}
} }
} }
@@ -1912,6 +2053,7 @@ int* r_extra_begin, int* r_extra_end, unsigned int* r_error)
int path_length; int path_length;
///if the start pos at the left boundary
if(local_y_start == 0) if(local_y_start == 0)
{ {
total_y_start = total_y_start + local_y_start; total_y_start = total_y_start + local_y_start;
@@ -2012,24 +2154,7 @@ int* r_extra_begin, int* r_extra_end, unsigned int* r_error)
} }
inline int double_error_threshold(int pre_threshold, int x_len)
{
pre_threshold = Adjust_Threshold(pre_threshold, x_len);
int threshold = pre_threshold * 2;
///may have some bugs
if(x_len >= 300 && threshold < THRESHOLD_MAX_SIZE)
{
threshold = THRESHOLD_MAX_SIZE;
}
if(threshold > THRESHOLD_MAX_SIZE)
{
threshold = THRESHOLD_MAX_SIZE;
}
return threshold;
}
inline void recalcate_window_back(overlap_region_alloc* overlap_list, All_reads* R_INF, inline void recalcate_window_back(overlap_region_alloc* overlap_list, All_reads* R_INF,
@@ -2716,7 +2841,7 @@ inline void recalcate_window_back(overlap_region_alloc* overlap_list, All_reads*
} }
inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_INF, inline void recalcate_window_simple(overlap_region_alloc* overlap_list, All_reads* R_INF,
UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read) UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read)
{ {
long long j, k, i; long long j, k, i;
@@ -3275,6 +3400,497 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN
inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_INF,
UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read)
{
long long j, k, i;
long long Len_x;
int threshold;
long long y_len;
long long currentIDLen;
long long matches;
long long y_id;
int y_strand;
long long y_readLen;
long long x_start;
long long x_end;
long long x_len;
long long total_y_start;
long long total_y_end;
long long y_start;
long long y_end;
long long Window_Len;
char* x_string;
char* y_string;
int end_site;
unsigned int error;
int real_y_start;
long long overlap_length;
int extra_begin, extra_end;
long long o_len;
///j负责遍历整个overlap list
for (j = 0; j < overlap_list->length; j++)
{
y_id = overlap_list->list[j].y_id;
y_strand = overlap_list->list[j].y_pos_strand;
y_readLen = Get_READ_LENGTH((*R_INF), y_id);
//i corresponding to each window of a overlap
//utilize the the end pos of pre-window in backwards
for (i = overlap_list->list[j].w_list_length - 1; i >= 0; i--)
{
///the first matched window
if(overlap_list->list[j].w_list[i].y_end != -1)
{
///note!!! need notification
///this is the actual end postion in ystring
total_y_start = overlap_list->list[j].w_list[i].y_end
- overlap_list->list[j].w_list[i].extra_begin + 1;
///k corresponding to all unmatched windows at the right side of overlap_list->list[j].w_list[i]
///so k starts from i + 1, and end to the first matched window
for (k = i + 1; k < overlap_list->list[j].w_list_length && overlap_list->list[j].w_list[k].y_end == -1; k++)
{
extra_begin = extra_end = 0;
///if y_start > y_readLen, direct terminate
if (total_y_start >= y_readLen)
{
break;
}
///there is no problem for x
x_start = overlap_list->list[j].w_list[k].x_start;
x_end = overlap_list->list[j].w_list[k].x_end;
x_len = x_end - x_start + 1;
///there are two potiential reasons for unmatched window:
///1. this window has a large number of differences
///2. DP does not start from the right offset
threshold = double_error_threshold(overlap_list->list[j].w_list[k].error_threshold, x_len);
y_start = total_y_start;
Window_Len = x_len + (threshold << 1);
if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF,
&extra_begin, &extra_end, &y_start, &o_len))
{
break;
}
if(o_len + threshold < x_len)
{
break;
}
fill_subregion(dumy->overlap_region, y_start, o_len, y_strand,
R_INF, y_id, extra_begin, extra_end);
x_string = g_read->seq + x_start;
y_string = dumy->overlap_region;
///note!!! need notification
end_site = Reserve_Banded_BPM(y_string, Window_Len, x_string, x_len, threshold, &error);
// unsigned int debug_error;
// int return_y_end, return_x_end;
// Reserve_Banded_BPM_Extension(y_string, Window_Len, x_string, x_len, threshold,
// &debug_error, &return_y_end, &return_x_end);
// if(return_x_end != x_len - 1 && end_site != -1)
// {
// fprintf(stderr, "error\n");
// }
// if(return_x_end == x_len - 1 && end_site == -1)
// {
// fprintf(stderr, "error\n");
// }
// if(return_x_end == x_len - 1 && end_site != -1)
// {
// if(return_y_end != end_site || debug_error != error)
// {
// fprintf(stderr, "error\n");
// }
// }
///error等于-1说明没匹配
if (error!=(unsigned int)-1)
{
overlap_list->list[j].w_list[k].cigar.length = -1;
overlap_list->list[j].w_list[k].y_start = y_start;
overlap_list->list[j].w_list[k].y_end = y_start + end_site;
overlap_list->list[j].w_list[k].error = (int)error;
///note!!! need notification
overlap_list->list[j].w_list[k].extra_begin = extra_begin;
overlap_list->list[j].w_list[k].extra_end = extra_end;
overlap_list->list[j].w_list[k].error_threshold = threshold;
overlap_list->list[j].align_length += x_len;
}
else
{
break;
}
///note!!! need notification
total_y_start = y_start + end_site - extra_begin + 1;
}
}
}
//i corresponding to each window of a overlap
//utilize the the start pos of next window in forward
for (i = 0; i < overlap_list->list[j].w_list_length; i++)
{
///find the first matched window, which should not be the first window
///the pre-window of this matched window must be unmatched
if(overlap_list->list[j].w_list[i].y_end != -1 && i != 0 && overlap_list->list[j].w_list[i - 1].y_end == -1)
{
///check if the start pos of this matched window has been calculated
if(overlap_list->list[j].w_list[i].cigar.length == -1)
{
///there is no problem for x
x_start = overlap_list->list[j].w_list[i].x_start;
x_end = overlap_list->list[j].w_list[i].x_end;
x_len = x_end - x_start + 1;
/****************************may have bugs********************************/
threshold = overlap_list->list[j].w_list[i].error_threshold;
/****************************may have bugs********************************/
/****************************may have bugs********************************/
///should not adjust threshold, since this window can be matched by the old threshold
///threshold = Adjust_Threshold(threshold, x_len);
/****************************may have bugs********************************/
Window_Len = x_len + (threshold << 1);
///y_start is the real y_start
y_start = overlap_list->list[j].w_list[i].y_start;
extra_begin = overlap_list->list[j].w_list[i].extra_begin;
extra_end = overlap_list->list[j].w_list[i].extra_end;
o_len = Window_Len - extra_end - extra_begin;
fill_subregion(dumy->overlap_region, y_start, o_len, y_strand,
R_INF, y_id, extra_begin, extra_end);
x_string = g_read->seq + x_start;
y_string = dumy->overlap_region;
///note!!! need notification
end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start,
&(dumy->path_length), dumy->matrix_bit, dumy->path,
overlap_list->list[j].w_list[i].error, overlap_list->list[j].w_list[i].y_end - y_start);
///y_start has already been calculated
if (error != (unsigned int)-1)
{
///this condition is always wrong
///in best case, real_y_start = threshold, end_site = Window_Len - threshold - 1
if (end_site == Window_Len - 1 || real_y_start == 0)
{
if(fix_boundary(x_string, x_len, threshold, y_start, real_y_start,
end_site, extra_begin, extra_end, y_id, Window_Len, R_INF, dumy,
y_strand, error, &y_start, &real_y_start, &end_site, &extra_begin,
&extra_end, &error))
{
overlap_list->list[j].w_list[i].error = error;
overlap_list->list[j].w_list[i].extra_begin = extra_begin;
overlap_list->list[j].w_list[i].extra_end = extra_end;
}
}
generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[i]),
&real_y_start, &end_site, &error, x_string, x_len, y_string);
// if(real_y_start < extra_begin || end_site >= Window_Len - extra_end)
// {
// fprintf(stderr, "\nreal_y_start: %d, extra_begin: %d\n",
// real_y_start, extra_begin);
// fprintf(stderr, "end_site: %d, Window_Len: %d, extra_end: %d\n",
// end_site, Window_Len, extra_end);
// }
///note!!! need notification
real_y_start = y_start + real_y_start - extra_begin;
overlap_list->list[j].w_list[i].y_start = real_y_start;
///I forget why don't reduce the extra_begin for y_end
///it seems extra_begin will be reduced at the end of this function
overlap_list->list[j].w_list[i].y_end = y_start + end_site;
overlap_list->list[j].w_list[i].error = error;
}
else
{
fprintf(stderr, "error\n");
}
}
else
{
real_y_start = overlap_list->list[j].w_list[i].y_start;
}
///the end pos for pre window is real_y_start - 1
total_y_end = real_y_start - 1;
///k遍历匹配window左侧所有不匹配的window
///如果i匹配，则k从i-1开始
///直到第一个匹配的window结束
///因为i!=0，所以k的大小不用担心
for (k = i - 1; k >= 0 && overlap_list->list[j].w_list[k].y_end == -1; k--)
{
///there is no problem in x
x_start = overlap_list->list[j].w_list[k].x_start;
x_end = overlap_list->list[j].w_list[k].x_end;
x_len = x_end - x_start + 1;
///there are two potiential reasons for unmatched window:
///1. this window has a large number of differences
///2. DP does not start from the right offset
threshold = double_error_threshold(overlap_list->list[j].w_list[k].error_threshold, x_len);
Window_Len = x_len + (threshold << 1);
if(total_y_end <= 0)
{
break;
}
///y_start might be less than 0
y_start = total_y_end - x_len + 1;
if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF,
&extra_begin, &extra_end, &y_start, &o_len))
{
break;
}
if(o_len + threshold < x_len)
{
break;
}
fill_subregion(dumy->overlap_region, y_start, o_len, y_strand,
R_INF, y_id, extra_begin, extra_end);
x_string = g_read->seq + x_start;
y_string = dumy->overlap_region;
///note!!! need notification
end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start,
&(dumy->path_length), dumy->matrix_bit, dumy->path, -1, -1);
if (error!=(unsigned int)-1)
{
///this condition is always wrong
///in best case, real_y_start = threshold, end_site = Window_Len - threshold - 1
if (end_site == Window_Len - 1 || real_y_start == 0)
{
fix_boundary(x_string, x_len, threshold, y_start, real_y_start, end_site,
extra_begin, extra_end, y_id, Window_Len, R_INF, dumy, y_strand, error,
&y_start, &real_y_start, &end_site,
&extra_begin, &extra_end, &error);
}
generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[k]),
&real_y_start, &end_site, &error, x_string, x_len, y_string);
// if(real_y_start < extra_begin || end_site >= Window_Len - extra_end)
// {
// fprintf(stderr, "\nreal_y_start: %d, extra_begin: %d\n",
// real_y_start, extra_begin);
// fprintf(stderr, "end_site: %d, Window_Len: %d, extra_end: %d\n",
// end_site, Window_Len, extra_end);
// }
///y_start has no shift, but y_end has shift
overlap_list->list[j].w_list[k].y_start = y_start + real_y_start - extra_begin;
overlap_list->list[j].w_list[k].y_end = y_start + end_site;
overlap_list->list[j].w_list[k].error = error;
overlap_list->list[j].align_length += x_len;
overlap_list->list[j].w_list[k].extra_begin = extra_begin;
overlap_list->list[j].w_list[k].extra_end = extra_end;
overlap_list->list[j].w_list[k].error_threshold = threshold;
}
else
{
break;
}
total_y_end = y_start + real_y_start - 1 - extra_begin;
}
}
}
}
overlap_list->mapped_overlaps_length = 0;
int pre_threshold;
long long tLen, tError;
double error_rate;
for (j = 0; j < overlap_list->length; j++)
{
y_id = overlap_list->list[j].y_id;
y_strand = overlap_list->list[j].y_pos_strand;
y_readLen = Get_READ_LENGTH((*R_INF), y_id);
overlap_length = overlap_list->list[j].x_pos_e - overlap_list->list[j].x_pos_s + 1;
overlap_list->list[j].is_match = 0;
///only calculate cigar for high quality overlaps
if (overlap_length * OVERLAP_THRESHOLD_FILTER <= overlap_list->list[j].align_length)
{
for (i = 0; i < overlap_list->list[j].w_list_length; i++)
{
///first we need to check if this window is matched
if(overlap_list->list[j].w_list[i].y_end != -1)
{
///second check if the cigar of this window has been got
if(overlap_list->list[j].w_list[i].cigar.length == -1)
{
///there is no problem for x
x_start = overlap_list->list[j].w_list[i].x_start;
x_end = overlap_list->list[j].w_list[i].x_end;
x_len = x_end - x_start + 1;
/****************************may have bugs********************************/
///threshold = x_len * THRESHOLD_RATE;
threshold = overlap_list->list[j].w_list[i].error_threshold;
/****************************may have bugs********************************/
/****************************may have bugs********************************/
///should not adjust threshold, since this window can be matched by the old threshold
///threshold = Adjust_Threshold(threshold, x_len);
/****************************may have bugs********************************/
Window_Len = x_len + (threshold << 1);
///y_start is the real y_start
///for the window with cigar, y_start has already reduced extra_begin
y_start = overlap_list->list[j].w_list[i].y_start;
extra_begin = overlap_list->list[j].w_list[i].extra_begin;
extra_end = overlap_list->list[j].w_list[i].extra_end;
o_len = Window_Len - extra_end - extra_begin;
fill_subregion(dumy->overlap_region, y_start, o_len, y_strand,
R_INF, y_id, extra_begin, extra_end);
x_string = g_read->seq + x_start;
y_string = dumy->overlap_region;
///note!!! need notification
end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start,
&(dumy->path_length), dumy->matrix_bit, dumy->path,
overlap_list->list[j].w_list[i].error, overlap_list->list[j].w_list[i].y_end - y_start);
// if(error != overlap_list->list[j].w_list[i].error)
// {
// fprintf(stderr, "error\n");
// }
if (error != (unsigned int)-1)
{
if (end_site == Window_Len - 1 || real_y_start == 0)
{
if(fix_boundary(x_string, x_len, threshold, y_start, real_y_start, end_site,
extra_begin, extra_end, y_id, Window_Len, R_INF, dumy, y_strand, error,
&y_start, &real_y_start, &end_site,
&extra_begin, &extra_end, &error))
{
overlap_list->list[j].w_list[i].error = error;
overlap_list->list[j].w_list[i].extra_begin = extra_begin;
overlap_list->list[j].w_list[i].extra_end = extra_end;
}
}
generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[i]),
&real_y_start, &end_site, &error, x_string, x_len, y_string);
// if(real_y_start < extra_begin || end_site >= Window_Len - extra_end)
// {
// fprintf(stderr, "\nreal_y_start: %d, extra_begin: %d\n",
// real_y_start, extra_begin);
// fprintf(stderr, "end_site: %d, Window_Len: %d, extra_end: %d\n",
// end_site, Window_Len, extra_end);
// }
///note!!! need notification
real_y_start = y_start + real_y_start - extra_begin;
overlap_list->list[j].w_list[i].y_start = real_y_start;
overlap_list->list[j].w_list[i].y_end = y_start + end_site - extra_begin;
overlap_list->list[j].w_list[i].error = error;
}
else
{
fprintf(stderr, "error\n");
}
}
else
{
overlap_list->list[j].w_list[i].y_end -= overlap_list->list[j].w_list[i].extra_begin;
}
}
}
///error_rate = trim_error_rate(overlap_list, j);
error_rate = non_trim_error_rate(overlap_list, j, R_INF, dumy, g_read);
///if(error_rate <= 0.015)
if(error_rate <= 0.03)
{
overlap_list->mapped_overlaps_length += overlap_length;
overlap_list->list[j].is_match = 1;
}
else if(error_rate <= 0.045)
{
overlap_list->list[j].is_match = 3;
}
}
}
}
/**
inline void adjust_alignment_windows(overlap_region* overlap, All_reads* R_INF)
{
long long y_id = overlap->y_id;
long long y_strand = overlap->y_pos_strand;
long long y_readLen = Get_READ_LENGTH((*R_INF), y_id);
long long overlap_length = overlap->x_pos_e - overlap->x_pos_s + 1;
long long i;
///note here we start from i = 1, instead of i = 0
for (i = 1; i < overlap->w_list_length; i++)
{
///that means we have both start pos and end pos
if(overlap->w_list[i].y_end != -1 && overlap->w_list[i].cigar.length != -1)
{
;
}
}
}
**/
inline void add_base_to_correct_read_directly(Correct_dumy* dumy, char base) inline void add_base_to_correct_read_directly(Correct_dumy* dumy, char base)
{ {
Hash_Table.cpp
+381 -79
View File
@@ -632,6 +632,225 @@ void append_overlap_region_alloc_from_existing(overlap_region_alloc* list, overl
} }
int append_inexact_overlap_region_alloc_back(overlap_region_alloc* list, overlap_region* tmp, All_reads* R_INF)
{
if (list->length + 1 > list->size)
{
list->size = list->size * 2;
list->list = (overlap_region*)realloc(list->list, sizeof(overlap_region)*list->size);
///新分配空间要初始化
memset(list->list + (list->size/2), 0, sizeof(overlap_region)*(list->size/2));
}
if (list->length!=0 &&
list->list[list->length - 1].y_id==tmp->y_id
)
{
if(list->list[list->length - 1].shared_seed >= tmp->shared_seed)
{
return 0;
}
else
{
list->length--;
}
}
if(tmp->x_pos_s <= tmp->y_pos_s)
{
tmp->y_pos_s = tmp->y_pos_s - tmp->x_pos_s;
tmp->x_pos_s = 0;
}
else
{
tmp->x_pos_s = tmp->x_pos_s - tmp->y_pos_s;
tmp->y_pos_s = 0;
}
long long x_right_length = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_e - 1;
long long y_right_length = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_e - 1;
if(x_right_length <= y_right_length)
{
tmp->x_pos_e = Get_READ_LENGTH((*R_INF), tmp->x_id) - 1;
tmp->y_pos_e = tmp->y_pos_e + x_right_length;
}
else
{
tmp->x_pos_e = tmp->x_pos_e + y_right_length;
tmp->y_pos_e = Get_READ_LENGTH((*R_INF), tmp->y_id) - 1;
}
if (tmp->x_pos_strand == 1)
{
list->list[list->length].x_id = tmp->x_id;
list->list[list->length].x_pos_e = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_s - 1;
list->list[list->length].x_pos_s = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_e - 1;
list->list[list->length].x_pos_strand = 0;
list->list[list->length].y_id = tmp->y_id;
list->list[list->length].y_pos_e = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_s - 1;
list->list[list->length].y_pos_s = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_e - 1;
list->list[list->length].y_pos_strand = 1;
resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2));
add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0);
long long distance_gap;
long long pre_distance_gap = 0;
long long i = 0;
for (i = 0; i < tmp->f_cigar.length; i++)
{
distance_gap = get_fake_gap_shift(&(tmp->f_cigar), i);
if(distance_gap != pre_distance_gap)
{
pre_distance_gap = distance_gap;
add_fake_cigar(&(list->list[list->length].f_cigar),
Get_READ_LENGTH((*R_INF), tmp->x_id) - get_fake_gap_pos(&(tmp->f_cigar), i) - 1,
pre_distance_gap);
}
}
if(get_fake_gap_pos(&(list->list[list->length].f_cigar),
list->list[list->length].f_cigar.length - 1) != list->list[list->length].x_pos_e)
{
add_fake_cigar(&(list->list[list->length].f_cigar),
list->list[list->length].x_pos_e,
get_fake_gap_shift(&(list->list[list->length].f_cigar),
list->list[list->length].f_cigar.length - 1));
}
/******************************for debug********************************/
// long long distance_self_pos = tmp->x_pos_e - tmp->x_pos_s;
// long long distance_pos = tmp->y_pos_e - tmp->y_pos_s;
// distance_gap = distance_pos - distance_self_pos;
// if(distance_gap !=
// get_fake_gap_shift(&(list->list[list->length].f_cigar),
// list->list[list->length].f_cigar.length - 1))
// {
// fprintf(stderr, "error\n");
// }
// long long distance_self_pos = tmp->x_pos_e - tmp->x_pos_s;
// long long distance_pos = tmp->y_pos_e - tmp->y_pos_s;
// distance_gap = distance_pos - distance_self_pos;
// fprintf(stderr, "\nx_s: %d, x_e: %d, y_s: %d, y_e: %d, distance_gap: %d, xLen: %d\n",
// list->list[list->length].x_pos_s, list->list[list->length].x_pos_e,
// list->list[list->length].y_pos_s, list->list[list->length].y_pos_e, distance_gap,
// Get_READ_LENGTH((*R_INF), tmp->x_id));
// for (i = 0; i < list->list[list->length].f_cigar.length; i++)
// {
// fprintf(stderr, "##i: %d, gap_pos_in_x: %d, gap_shift: %d\n",
// i, get_fake_gap_pos(&(list->list[list->length].f_cigar), i),
// get_fake_gap_shift(&(list->list[list->length].f_cigar), i));
// }
// for (i = 0; i < tmp->f_cigar.length; i++)
// {
// fprintf(stderr, "**i: %d, gap_pos_in_x: %d, gap_shift: %d\n",
// i, get_fake_gap_pos(&(tmp->f_cigar), i),
// get_fake_gap_shift(&(tmp->f_cigar), i));
// }
/******************************for debug********************************/
}
else
{
list->list[list->length].x_id = tmp->x_id;
list->list[list->length].x_pos_e = tmp->x_pos_e;
list->list[list->length].x_pos_s = tmp->x_pos_s;
list->list[list->length].x_pos_strand = tmp->x_pos_strand;
list->list[list->length].y_id = tmp->y_id;
list->list[list->length].y_pos_e = tmp->y_pos_e;
list->list[list->length].y_pos_s = tmp->y_pos_s;
list->list[list->length].y_pos_strand = tmp->y_pos_strand;
resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2));
add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0);
long long distance_self_pos = tmp->x_pos_e - tmp->x_pos_s;
long long distance_pos = tmp->y_pos_e - tmp->y_pos_s;
long long init_distance_gap = distance_pos - distance_self_pos;
long long pre_distance_gap = init_distance_gap;
long long distance_gap;
long long i = 0;
for (i = tmp->f_cigar.length - 1; i >= 0; i--)
{
distance_gap = get_fake_gap_shift(&(tmp->f_cigar), i);
if(distance_gap != pre_distance_gap)
{
pre_distance_gap = distance_gap;
add_fake_cigar(&(list->list[list->length].f_cigar),
get_fake_gap_pos(&(tmp->f_cigar), i), init_distance_gap - pre_distance_gap);
}
}
if(get_fake_gap_pos(&(list->list[list->length].f_cigar),
list->list[list->length].f_cigar.length - 1) != list->list[list->length].x_pos_e)
{
add_fake_cigar(&(list->list[list->length].f_cigar),
list->list[list->length].x_pos_e,
get_fake_gap_shift(&(list->list[list->length].f_cigar),
list->list[list->length].f_cigar.length - 1));
}
/******************************for debug********************************/
// distance_self_pos = tmp->x_pos_e - tmp->x_pos_s;
// distance_pos = tmp->y_pos_e - tmp->y_pos_s;
// distance_gap = distance_pos - distance_self_pos;
// if(distance_gap !=
// get_fake_gap_shift(&(list->list[list->length].f_cigar),
// list->list[list->length].f_cigar.length - 1))
// {
// fprintf(stderr, "error\n");
// }
// distance_self_pos = tmp->x_pos_e - tmp->x_pos_s;
// distance_pos = tmp->y_pos_e - tmp->y_pos_s;
// distance_gap = distance_pos - distance_self_pos;
// fprintf(stderr, "\nx_s: %d, x_e: %d, y_s: %d, y_e: %d, distance_gap: %d\n",
// list->list[list->length].x_pos_s, list->list[list->length].x_pos_e,
// list->list[list->length].y_pos_s, list->list[list->length].y_pos_e, distance_gap);
// for (i = 0; i < list->list[list->length].f_cigar.length; i++)
// {
// fprintf(stderr, "##i: %d, gap_pos_in_x: %d, gap_shift: %d\n",
// i, get_fake_gap_pos(&(list->list[list->length].f_cigar), i),
// get_fake_gap_shift(&(list->list[list->length].f_cigar), i));
// }
// for (i = 0; i < tmp->f_cigar.length; i++)
// {
// fprintf(stderr, "**i: %d, gap_pos_in_x: %d, gap_shift: %d\n",
// i, get_fake_gap_pos(&(tmp->f_cigar), i),
// get_fake_gap_shift(&(tmp->f_cigar), i));
// }
/******************************for debug********************************/
}
list->list[list->length].shared_seed = tmp->shared_seed;
list->list[list->length].align_length = 0;
list->list[list->length].is_match = 0;
list->list[list->length].non_homopolymer_errors = 0;
list->list[list->length].strong = 0;
list->length++;
return 1;
}
int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, All_reads* R_INF) int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, All_reads* R_INF)
{ {
@@ -1278,6 +1497,164 @@ long long y_beg, long long y_end, long long yLen)
return x_end - x_beg + 1; return x_end - x_beg + 1;
} }
void chain_DP_back(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result,
double band_width_threshold)
{
long long i, j;
long long self_pos, pos, max_j, max_i, max_score, score, n_skip;
long long distance_pos, distance_self_pos, distance_gap, log_distance_gap, distance_min;
///double band_width_threshold = 0.05;
double band_width_penalty = 1 / band_width_threshold;
long long min_score = k_mer_length;
long long max_indels, max_self_length;
double gap_rate;
long long total_indels, total_self_length;
resize_Chain_Data(dp, a_n);
// fill the score and backtrack arrays
for (i = 0; i < a_n; ++i)
{
pos = a[i].offset;
self_pos = a[i].self_offset;
max_j = -1;
max_score = min_score;
n_skip = 0;
max_indels = 0;
max_self_length = 0;
///may have a pre-cut condition for j
for (j = i - 1; j >= 0; --j)
{
distance_pos = pos - a[j].offset;
distance_self_pos = self_pos - a[j].self_offset;
///a has been sorted by a[].offset
///note for a, we do not have any two elements that have both equal offsets and self_offsets
///but there maybe two elements that have equal offsets or equal self_offsets
if(distance_pos == 0 || distance_self_pos <= 0)
{
continue;
}
distance_gap = distance_pos > distance_self_pos? distance_pos - distance_self_pos : distance_self_pos - distance_pos;
total_indels = dp->indels[j] + distance_gap;
total_self_length = dp->self_length[j] + distance_self_pos;
if(total_indels > band_width_threshold * total_self_length)
{
continue;
}
distance_min = distance_pos < distance_self_pos? distance_pos:distance_self_pos;
score = distance_min < min_score? distance_min : min_score;
/**
log_distance_gap = distance_gap? ilog2_32(distance_gap) : 0;
score -= (long long)(distance_gap * 0.01 * min_score) + (log_distance_gap/2);
**/
gap_rate = (double)((double)(total_indels)/(double)(total_self_length));
///if the gap rate > 0.05, score will be negative
score -= (long long)(gap_rate * score * band_width_penalty);
score += dp->score[j];
if(score > max_score)
{
max_score = score;
max_j = j;
max_indels = total_indels;
max_self_length = total_self_length;
if (n_skip > 0)
{
n_skip--;
}
}
}
dp->score[i] = max_score;
dp->pre[i] = max_j;
dp->indels[i] = max_indels;
dp->self_length[i] = max_self_length;
}
///debug_chain(a, a_n, dp);
max_score = -1;
max_i = -1;
for (i = 0; i < a_n; ++i)
{
if(dp->score[i] > max_score)
{
max_score = dp->score[i];
max_i = i;
}
}
clear_fake_cigar(&(result->f_cigar));
i = max_i;
result->x_pos_e = a[i].self_offset;
result->y_pos_e = a[i].offset;
result->shared_seed = max_score;
distance_self_pos = result->x_pos_e - a[i].self_offset;
distance_pos = result->y_pos_e - a[i].offset;
long long pre_distance_gap = distance_pos - distance_self_pos;
///record first site
///the length of f_cigar should be at least 1
add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap);
long long chainLen = 0;
if(result->x_pos_strand == 1)
{
while (i >= 0)
{
distance_self_pos = result->x_pos_e - a[i].self_offset;
distance_pos = result->y_pos_e - a[i].offset;
distance_gap = distance_pos - distance_self_pos;
if(distance_gap != pre_distance_gap)
{
pre_distance_gap = distance_gap;
///record this site
add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap);
}
chainLen++;
result->x_pos_s = a[i].self_offset;
result->y_pos_s = a[i].offset;
i = dp->pre[i];
}
}
else
{
while (i >= 0)
{
distance_self_pos = result->x_pos_e - a[i].self_offset;
distance_pos = result->y_pos_e - a[i].offset;
distance_gap = distance_pos - distance_self_pos;
if(distance_gap == pre_distance_gap)
{
result->f_cigar.length--;
add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap);
}
else
{
pre_distance_gap = distance_gap;
add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap);
}
chainLen++;
result->x_pos_s = a[i].self_offset;
result->y_pos_s = a[i].offset;
i = dp->pre[i];
}
}
}
///double band_width_threshold = 0.05; ///double band_width_threshold = 0.05;
void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result,
double band_width_threshold, int max_skip, int x_readLen, int y_readLen) double band_width_threshold, int max_skip, int x_readLen, int y_readLen)
@@ -1618,95 +1995,20 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh
continue; continue;
} }
// if(
// ((memcmp("m64016_190918_162737/130811282/ccs",
// Get_NAME((*R_INF), tmp_region.x_id), Get_NAME_LENGTH((*R_INF), tmp_region.x_id)) == 0)
// &&
// (memcmp("m64016_190918_162737/179635219/ccs",
// Get_NAME((*R_INF), tmp_region.y_id), Get_NAME_LENGTH((*R_INF), tmp_region.y_id)) == 0))
// ||
// ((memcmp("m64016_190918_162737/179635219/ccs",
// Get_NAME((*R_INF), tmp_region.x_id), Get_NAME_LENGTH((*R_INF), tmp_region.x_id)) == 0)
// &&
// (memcmp("m64016_190918_162737/130811282/ccs",
// Get_NAME((*R_INF), tmp_region.y_id), Get_NAME_LENGTH((*R_INF), tmp_region.y_id)) == 0)))
// {
// fprintf(stderr, "****************x_name: %.*s****************\n",
// Get_NAME_LENGTH((*R_INF), tmp_region.x_id), Get_NAME((*R_INF), tmp_region.x_id));
// fprintf(stderr, "****************y_name: %.*s****************\n",
// Get_NAME_LENGTH((*R_INF), tmp_region.y_id), Get_NAME((*R_INF), tmp_region.y_id));
// k_mer_hit* k_list = candidates->list + sub_region_beg;
// long long k_listLen = sub_region_end - sub_region_beg + 1;
// long long k = 0;
// fprintf(stderr, "k_listLen: %d\n", k_listLen);
// for (k = 0; k < k_listLen; k++)
// {
// fprintf(stderr, "k: %d, readID: %d, strand: %d, offset: %d, self_offset: %d\n",
// k, k_list[k].readID, k_list[k].strand, k_list[k].offset, k_list[k].self_offset);
// }
// }
chain_DP(candidates->list + sub_region_beg, chain_DP(candidates->list + sub_region_beg,
sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold, sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold,
50, Get_READ_LENGTH((*R_INF), tmp_region.x_id), Get_READ_LENGTH((*R_INF), tmp_region.y_id)); 50, Get_READ_LENGTH((*R_INF), tmp_region.x_id), Get_READ_LENGTH((*R_INF), tmp_region.y_id));
// chain_DP_back(candidates->list + sub_region_beg,
// sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold);
// if(
// ((memcmp("m64016_190918_162737/130811282/ccs",
// Get_NAME((*R_INF), tmp_region.x_id), Get_NAME_LENGTH((*R_INF), tmp_region.x_id)) == 0)
// &&
// (memcmp("m64016_190918_162737/179635219/ccs",
// Get_NAME((*R_INF), tmp_region.y_id), Get_NAME_LENGTH((*R_INF), tmp_region.y_id)) == 0))
// ||
// ((memcmp("m64016_190918_162737/179635219/ccs",
// Get_NAME((*R_INF), tmp_region.x_id), Get_NAME_LENGTH((*R_INF), tmp_region.x_id)) == 0)
// &&
// (memcmp("m64016_190918_162737/130811282/ccs",
// Get_NAME((*R_INF), tmp_region.y_id), Get_NAME_LENGTH((*R_INF), tmp_region.y_id)) == 0)))
// {
// fprintf(stderr, "****************x_name: %.*s****************\n",
// Get_NAME_LENGTH((*R_INF), tmp_region.x_id), Get_NAME((*R_INF), tmp_region.x_id));
// fprintf(stderr, "****************y_name: %.*s****************\n",
// Get_NAME_LENGTH((*R_INF), tmp_region.y_id), Get_NAME((*R_INF), tmp_region.y_id));
// k_mer_hit* k_list = candidates->list + sub_region_beg;
// long long k_listLen = sub_region_end - sub_region_beg + 1;
// long long k = 0;
// fprintf(stderr, "k_listLen: %d\n", k_listLen);
// for (k = 0; k < k_listLen; k++)
// {
// fprintf(stderr, "k: %d, score: %d, pre: %d, indels: %d, self_length: %d\n",
// k, candidates->chainDP.score[k], candidates->chainDP.pre[k],
// candidates->chainDP.indels[k], candidates->chainDP.self_length[k]);
// }
// }
///自己和自己重叠的要排除 ///自己和自己重叠的要排除
///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1) ///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1)
if (tmp_region.x_id != tmp_region.y_id) if (tmp_region.x_id != tmp_region.y_id)
{ {
///append_overlap_region_alloc(overlap_list, &tmp_region, R_INF);
/**
if(append_inexact_overlap_region_alloc(overlap_list, &tmp_region, R_INF))
{
debug_seed_offset(candidates->list + sub_region_beg,
sub_region_end - sub_region_beg + 1,
&(candidates->chainDP), &(overlap_list->list[overlap_list->length - 1]), R_INF);
}
**/
append_inexact_overlap_region_alloc(overlap_list, &tmp_region, R_INF); append_inexact_overlap_region_alloc(overlap_list, &tmp_region, R_INF);
///append_inexact_overlap_region_alloc_back(overlap_list, &tmp_region, R_INF);
} }
} }
Levenshtein_distance.h
+388
View File
@@ -13,6 +13,394 @@ typedef uint64_t Word;
typedef uint32_t Word_32; typedef uint32_t Word_32;
typedef struct
{
/**
* vec[num_words-1], vec[num_words-2], ..., vec[1], vec[0]
* num_bits-1, num_bits-2, num_bits-3, ..., 8, 7, 6, 5, 4, 3, 2, 1, 0
**/
Word* vec;
Word num_words;
Word num_bits;
///for vec[num_words-1]
Word last_bit_mode;
Word last_bit_shift;
} bit_vectors;
inline void init_bit_vector(bit_vectors* b, Word num_bits)
{
b->num_bits = num_bits;
b->num_words = ((b->num_bits)>>6);
b->last_bit_mode = (Word)-1;
b->last_bit_shift = b->num_bits&63;
if(b->last_bit_shift != 0)
{
b->num_words++;
b->last_bit_mode = b->last_bit_mode >> (64 - b->last_bit_shift);
}
b->vec = (Word*)calloc(b->num_words, sizeof(Word));
}
inline void destory_bit_vector(bit_vectors* b, Word num_bits)
{
free(b->vec);
}
///&
inline int AND_bit_vector(bit_vectors* x, bit_vectors* y, bit_vectors* dest)
{
if(x->num_bits != y->num_bits || x->num_bits != dest->num_bits) return 0;
Word i;
for (i = 0; i < x->num_words; i++)
{
dest->vec[i] = x->vec[i] & y->vec[i];
}
if(dest->num_words > 0)
{
dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode;
}
return 1;
}
///|
inline int OR_bit_vector(bit_vectors* x, bit_vectors* y, bit_vectors* dest)
{
if(x->num_bits != y->num_bits || x->num_bits != dest->num_bits) return 0;
Word i;
for (i = 0; i < x->num_words; i++)
{
dest->vec[i] = x->vec[i] | y->vec[i];
}
if(dest->num_words > 0)
{
dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode;
}
return 1;
}
///^
inline int XOR_bit_vector(bit_vectors* x, bit_vectors* y, bit_vectors* dest)
{
if(x->num_bits != y->num_bits || x->num_bits != dest->num_bits) return 0;
Word i;
for (i = 0; i < x->num_words; i++)
{
dest->vec[i] = x->vec[i] ^ y->vec[i];
}
if(dest->num_words > 0)
{
dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode;
}
return 1;
}
///~
inline int NOT_bit_vector(bit_vectors* source, bit_vectors* dest)
{
if(source->num_bits != dest->num_bits) return 0;
Word i;
for (i = 0; i < source->num_words; i++)
{
dest->vec[i] = ~(source->vec[i]);
}
if(dest->num_words > 0)
{
dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode;
}
return 1;
}
///<< 1
inline int L_shift_1_bit_vector(bit_vectors* source, bit_vectors* dest)
{
/**
* vec[num_words-1], vec[num_words-2], ..., vec[1], vec[0]
* num_bits-1, num_bits-2, num_bits-3, ..., 8, 7, 6, 5, 4, 3, 2, 1, 0
**/
if(source->num_bits != dest->num_bits || source->num_words < 1) return 0;
Word i;
for (i = source->num_words - 1; i >= 1; i--)
{
dest->vec[i] = (source->vec[i])<<1;
dest->vec[i] = dest->vec[i] | ((source->vec[i-1])>>63);
}
dest->vec[0] = (source->vec[0])<<1;
dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode;
return 1;
}
inline void get_error(int t_length, int errthold, int init_err, Word VP, Word VN,
unsigned int* return_err, int* back_site)
{
(*return_err) = (unsigned int)-1;
int site = t_length - 1;
int return_site = -1;
///p_length大部分情况下应该是t_length + 2 * errthold，这是i要小于last_high = 2 * errthold
///也就是p_length - t_length
///那么当p_length < t_length + 2 * errthold, available_i也应该是这个值
///int available_i = p_length - t_length;
int available_i = 2 * errthold;
if ((init_err <= errthold) && (init_err <= (*return_err)))
{
(*return_err) = init_err;
return_site = site;
}
int i = 0;
unsigned int ungap_error = (unsigned int)-1;
while (i < available_i)
{
init_err = init_err + ((VP >> i)&(Word)1);
init_err = init_err - ((VN >> i)&(Word)1);
++i;
if ((init_err <= errthold) && (init_err <= *return_err))
{
*return_err = init_err;
return_site = site + i;
}
/****************************may have bugs********************************/
if(i == errthold)
{
ungap_error = init_err;
}
/****************************may have bugs********************************/
}
/****************************may have bugs********************************/
if((ungap_error<=errthold) && (ungap_error == (*return_err)))
{
return_site = site + errthold;
}
/****************************may have bugs********************************/
(*back_site) = return_site;
}
inline int Reserve_Banded_BPM_Extension
(char *pattern, int p_length, char *text, int t_length, unsigned short errthold,
unsigned int* return_err, int* return_p_end, int* return_t_end)
{
(*return_err) = (unsigned int)-1;
(*return_p_end) = -1;
(*return_t_end) = -1;
Word Peq[256];
unsigned int line_error = (unsigned int)-1;
int return_site;
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<band_length; i++)
{
Peq[pattern[i]] = Peq[pattern[i]] | tmp_Peq_1;
tmp_Peq_1 = tmp_Peq_1 << 1;
}
Peq_A = Peq['A'];
Peq_C = Peq['C'];
Peq_T = Peq['T'];
Peq_G = Peq['G'];
memset(Peq, 0, sizeof(Word)* 256);
Peq['A'] = Peq_A;
Peq['C'] = Peq_C;
Peq['T'] = Peq_T;
Peq['G'] = Peq_G;
Word Mask = ((Word)1 << (errthold << 1));
Word VP = 0;
Word VN = 0;
Word X = 0;
Word D0 = 0;
Word HN = 0;
Word HP = 0;
i = 0;
int err = 0;
Word err_mask = (Word)1;
int i_bd = (errthold << 1);
int last_high = (errthold << 1);
int t_length_1 = t_length - 1;
while (i<t_length_1)
{
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 (*return_t_end);
}
}
get_error(i + 1, errthold, err, VP, VN, &line_error, &return_site);
if(line_error != (unsigned int)-1)
{
(*return_t_end) = i;
(*return_p_end) = return_site;
(*return_err) = line_error;
}
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 (*return_t_end);
}
}
///i = t_length - 1
get_error(i + 1, errthold, err, VP, VN, &line_error, &return_site);
if(line_error != (unsigned int)-1)
{
(*return_t_end) = i;
(*return_p_end) = return_site;
(*return_err) = line_error;
}
return (*return_t_end);
}
inline void reverse_string(char* str, int strLen)
{
int i, Len;
char k;
Len = strLen / 2;
for (i = 0; i < Len; i++)
{
k = str[i];
str[i] = str[strLen - i - 1];
str[strLen - i - 1] = k;
}
}
inline int alignment_extension(char *pattern, int p_length, char *text, int t_length,
unsigned short errthold, int direction, unsigned int* return_err, int* return_p_end,
int* return_t_end, int* return_aligned_t_len)
{
(*return_aligned_t_len) = 0;
if(direction == 0)
{
Reserve_Banded_BPM_Extension(pattern, p_length, text, t_length, errthold, return_err,
return_p_end, return_t_end);
if((*return_p_end) != -1 && (*return_t_end) != -1)
{
(*return_aligned_t_len) = (*return_t_end) + 1;
return 1;
}
else
{
return -1;
}
}
else
{
reverse_string(pattern, p_length);
reverse_string(text, t_length);
Reserve_Banded_BPM_Extension(pattern, p_length, text, t_length, errthold, return_err,
return_p_end, return_t_end);
reverse_string(pattern, p_length);
reverse_string(text, t_length);
if((*return_p_end) != -1 && (*return_t_end) != -1)
{
(*return_aligned_t_len) = (*return_t_end) + 1;
(*return_p_end) = p_length - (*return_p_end);
(*return_t_end) = t_length - (*return_t_end);
return 1;
}
else
{
return -1;
}
}
}
Overlaps.cpp
+1 -1
View File
@@ -9890,7 +9890,7 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
// debug_info_of_specfic_node("m64016_190918_162737/141297762/ccs", sg); // debug_info_of_specfic_node("m64016_190918_162737/141297762/ccs", sg);
out: out:
///output_tips(sg, &R_INF); output_tips(sg, &R_INF);
output_unitig_graph(sg, coverage_cut, output_file_name, n_read); output_unitig_graph(sg, coverage_cut, output_file_name, n_read);