xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-09-23 05:58:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

prefilter

Browse Source
This commit is contained in:
Haoyu Cheng
2019-07-09 15:11:04 -04:00
parent 1892b4d2db
commit 4eddb02dbd
5 changed files with 560 additions and 66 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Correct.cpp
+453 -43
View File
@@ -12,10 +12,7 @@ long long T_total_mis=0;
pthread_mutex_t debug_statistics ;
#define MAX(x, y) ((x >= y)?x:y)
#define MIN(x, y) ((x <= y)?x:y)
#define OVERLAP(x_start, x_end, y_start, y_end) (MIN(x_end, y_end) - MAX(x_start, y_start) + 1)
///#define OVERLAP(x_start, x_end, y_start, y_end) MIN(x_end, y_end) - MAX(x_start, y_start) + 1
@@ -392,6 +389,7 @@ char* r_string)
int groupLen = 0;
int return_sites[GROUP_SIZE];
unsigned int return_sites_error[GROUP_SIZE];
uint64_t overlapID[GROUP_SIZE];
///这些是整个window被完全覆盖的
for (i = 0; i < dumy->length; i++)
@@ -402,6 +400,13 @@ char* r_string)
x_start = window_start;
///y上的相对位置
y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s;
// /************需要注释掉**********/
// if(x_start < overlap_list->list[currentID].x_pos_s)
// {
// fprintf(stderr, "ERROR\n");
// }
// /************需要注释掉**********/
///y上的起始
y_start = y_start - THRESHOLD;
if (y_start < 0)
@@ -412,6 +417,13 @@ char* r_string)
currentIDLen = Get_READ_LENGTH((*R_INF), overlap_list->list[currentID].y_id);
///不能超过y的剩余长度
o_len = MIN(Window_Len, currentIDLen - y_start);
// /************需要注释掉**********/
// if(o_len < x_len)
// {
// fprintf(stderr, "ERROR\n");
// }
// /************需要注释掉**********/
///这个长度足够，可以用来一起比
@@ -420,6 +432,8 @@ char* r_string)
recover_UC_Read_sub_region(dumy->overlap_region_group[groupLen], y_start, o_len, overlap_list->list[currentID].y_pos_strand,
R_INF, overlap_list->list[currentID].y_id);
overlapID[groupLen] = currentID;
x_string = r_string + x_start;
@@ -432,30 +446,29 @@ char* r_string)
return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE);
groupLen = 0;
if (return_sites_error[0]!=(unsigned int)-1)
{
overlap_list->list[overlapID[0]].align_length += x_len;
}
if (return_sites_error[1]!=(unsigned int)-1)
{
overlap_list->list[overlapID[1]].align_length += x_len;
}
if (return_sites_error[2]!=(unsigned int)-1)
{
overlap_list->list[overlapID[2]].align_length += x_len;
}
if (return_sites_error[3]!=(unsigned int)-1)
{
overlap_list->list[overlapID[3]].align_length += x_len;
}
// end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], o_len, x_string, x_len, THRESHOLD, &error);
// if (error!=return_sites_error[0])
// {
// fprintf(stderr, "error\n");
// }
// end_site = Reserve_Banded_BPM(dumy->overlap_region_group[1], o_len, x_string, x_len, THRESHOLD, &error);
// if (error!=return_sites_error[1])
// {
// fprintf(stderr, "error\n");
// }
// end_site = Reserve_Banded_BPM(dumy->overlap_region_group[2], o_len, x_string, x_len, THRESHOLD, &error);
// if (error!=return_sites_error[2])
// {
// fprintf(stderr, "error\n");
// }
// end_site = Reserve_Banded_BPM(dumy->overlap_region_group[3], o_len, x_string, x_len, THRESHOLD, &error);
// if (error!=return_sites_error[3])
// {
// fprintf(stderr, "error\n");
// }
// /************需要注释掉**********/
/************需要注释掉**********/
// for (size_t ijk = 0; ijk < GROUP_SIZE; ijk++)
// {
// if (return_sites_error[ijk]!=(unsigned int)-1)
@@ -469,7 +482,7 @@ char* r_string)
// test_edit_distance_by_edlib(x_string, dumy->overlap_region_group[ijk], x_len, o_len, THRESHOLD,
// return_sites_error[ijk], &total_mis);
// }
// /************需要注释掉********* */
/************需要注释掉********* */
}
@@ -486,8 +499,14 @@ char* r_string)
///o_len < Window_Len, 说明y的长度不够，需要在y后面补N
memset (y_string + o_len, 0, Window_Len - o_len);
end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, THRESHOLD, &error);
if (error!=(unsigned int)-1)
{
overlap_list->list[currentID].align_length += x_len;
}
// /************需要注释掉**********/
/************需要注释掉**********/
// if (error!=(unsigned int)-1)
// {
// total_match++;
@@ -497,7 +516,7 @@ char* r_string)
// total_unmatch++;
// }
// test_edit_distance_by_edlib(x_string, y_string, x_len, o_len, THRESHOLD, error, &total_mis);
// /************需要注释掉********* */
/************需要注释掉********* */
}
}
@@ -506,7 +525,13 @@ char* r_string)
if (groupLen == 1)
{
end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], o_len, x_string, x_len, THRESHOLD, &error);
// /************需要注释掉**********/
if (error!=(unsigned int)-1)
{
overlap_list->list[overlapID[0]].align_length += x_len;
}
/************需要注释掉**********/
// if (error!=(unsigned int)-1)
// {
// total_match++;
@@ -516,14 +541,22 @@ char* r_string)
// total_unmatch++;
// }
// test_edit_distance_by_edlib(x_string, dumy->overlap_region_group[0], x_len, o_len, THRESHOLD, error, &total_mis);
// /************需要注释掉********* */
/************需要注释掉********* */
}
else if (groupLen > 1)
{
Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1],
dumy->overlap_region_group[2], dumy->overlap_region_group[3], o_len, x_string, x_len,
return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE);
// /************需要注释掉**********/
for (i = 0; i < groupLen; i++)
{
if (return_sites_error[i]!=(unsigned int)-1)
{
overlap_list->list[overlapID[i]].align_length += x_len;
}
}
/************需要注释掉**********/
// for (size_t ijk = 0; ijk < groupLen; ijk++)
// {
// if (return_sites_error[ijk]!=(unsigned int)-1)
@@ -537,7 +570,7 @@ char* r_string)
// test_edit_distance_by_edlib(x_string, dumy->overlap_region_group[ijk], x_len, o_len, THRESHOLD,
// return_sites_error[ijk], &total_mis);
// }
// /************需要注释掉**********/
/************需要注释掉**********/
groupLen = 0;
}
@@ -554,13 +587,23 @@ char* r_string)
///这个是和当前窗口重叠的长度
x_len = x_end - x_start + 1;
if (x_len <= 0)
{
fprintf(stderr, "ERROR\n");
}
// /************需要注释掉**********/
// if (x_len <= 0)
// {
// fprintf(stderr, "ERROR\n");
// }
// /************需要注释掉**********/
///y上的相对位置
y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s;
// /************需要注释掉**********/
// if(x_start < overlap_list->list[currentID].x_pos_s)
// {
// fprintf(stderr, "ERROR\n");
// }
// /************需要注释掉**********/
///y上的起始
y_start = y_start - THRESHOLD;
if (y_start < 0)
@@ -575,6 +618,15 @@ char* r_string)
///不能超过y的剩余长度
Window_Len = x_len + (THRESHOLD << 1);
o_len = MIN(Window_Len, currentIDLen - y_start);
// /************需要注释掉**********/
// if(o_len < x_len)
// {
// fprintf(stderr, "ERROR\n");
// }
// /************需要注释掉**********/
recover_UC_Read_sub_region(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand,
R_INF, overlap_list->list[currentID].y_id);
@@ -590,8 +642,13 @@ char* r_string)
}
end_site = Reserve_Banded_BPM(y_string, o_len, x_string, x_len, THRESHOLD, &error);
if (error!=(unsigned int)-1)
{
overlap_list->list[currentID].align_length += x_len;
}
// /************需要注释掉********* */
/************需要注释掉********* */
// if (error!=(unsigned int)-1)
// {
// total_match++;
@@ -601,19 +658,19 @@ char* r_string)
// total_unmatch++;
// }
// test_edit_distance_by_edlib(x_string, y_string, x_len, o_len, THRESHOLD, error, &total_mis);
// /************需要注释掉********* */
/************需要注释掉********* */
}
// /************需要注释掉********* */
/************需要注释掉********* */
// pthread_mutex_lock(&debug_statistics);
// T_total_match = T_total_match + total_match;
// T_total_unmatch = T_total_unmatch + total_unmatch;
// T_total_mis = T_total_mis + total_mis;
// pthread_mutex_unlock(&debug_statistics);
// /************需要注释掉********* */
/************需要注释掉********* */
}
@@ -670,12 +727,12 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Re
// /************需要注释掉********* */
/************需要注释掉********* */
// pthread_mutex_lock(&debug_statistics);
// fprintf(stderr, "total_match: %u, total_unmatch: %u, total_mis: %u\n",
// T_total_match, T_total_unmatch, T_total_mis);
// pthread_mutex_unlock(&debug_statistics);
// /************需要注释掉********* */
/************需要注释掉********* */
}
@@ -710,4 +767,357 @@ void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list)
list->overlapID = (uint64_t*)realloc(list->overlapID, list->size*sizeof(uint64_t));
}
}
}
void pre_filter_by_nearby_single(k_mer_pos* new_n_list, k_mer_pos* old_n_list, uint64_t n_length, uint64_t n_end_pos, UC_Read* g_read,
All_reads* R_INF, Correct_dumy* dumy, uint64_t* new_n_length)
{
(*new_n_length) = 0;
///这种就是0，new_n_list也不需要有数据
if (n_length == 0)
{
return;
}
char* x_string = NULL;
char* y_string = NULL;
long long x_offset = (long long)(n_end_pos) - WINDOW + 1;
if (x_offset < 0)
{
x_offset = 0;
}
///x_length是x上待验证区间的实际长度
///x_offset是x上待验证区间的起始位置
///如果是向前取待验证区间，那么x_length至少大于等于40 (compressed k-mer长度，不是实际长度)
long long x_length = n_end_pos - x_offset + 1;
///这种也不要过滤了，直接把old_n_list全部赋值过来就好了
if (x_length < (WINDOW/2))
{
(*new_n_length) = n_length;
memcpy(new_n_list, old_n_list, sizeof(k_mer_pos)*n_length);
return;
}
long long Window_Len = x_length + (THRESHOLD << 1);
x_string = g_read->seq + x_offset;
long long y_offset;
long long y_length;
long long i = 0;
long long y_read_length;
int end_site;
unsigned int error;
for (i = 0; i < n_length; i++)
{
////old_n_list[i].offset是y上k-mer的结束位置
///n_end_pos是x上k-mer的结束位置
///x_length是x上区间长度
y_offset = (long long)(old_n_list[i].offset) - x_length;
///这种情况下弃疗
if (y_offset < 0)
{
new_n_list[(*new_n_length)].readID = old_n_list[i].readID;
new_n_list[(*new_n_length)].offset = old_n_list[i].offset;
(*new_n_length)++;
continue;
}
y_offset = y_offset - THRESHOLD;
///还能抢救
if (y_offset < 0)
{
y_offset = 0;
}
///y的read的总长度
y_read_length = Get_READ_LENGTH((*R_INF), old_n_list[i].readID);
y_length = MIN(Window_Len, y_read_length - y_offset);
///如果y的长度比x长度还小，那就直接弃疗了
if (y_length < x_length)
{
new_n_list[(*new_n_length)].readID = old_n_list[i].readID;
new_n_list[(*new_n_length)].offset = old_n_list[i].offset;
(*new_n_length)++;
continue;
}
///y的方向都是0，因为索引里都是0
recover_UC_Read_sub_region(dumy->overlap_region, y_offset, y_length, 0, R_INF, old_n_list[i].readID);
y_string = dumy->overlap_region;
memset (y_string + y_length, 0, Window_Len - y_length);
end_site = Reserve_Banded_BPM(y_string, y_length, x_string, x_length, THRESHOLD, &error);
if (error!=(unsigned int)-1)
{
new_n_list[(*new_n_length)].readID = old_n_list[i].readID;
new_n_list[(*new_n_length)].offset = old_n_list[i].offset;
(*new_n_length)++;
}
}
}
void pre_filter_by_nearby(k_mer_pos* new_n_list, k_mer_pos* old_n_list, uint64_t n_length, uint64_t n_end_pos, UC_Read* g_read,
All_reads* R_INF, Correct_dumy* dumy, uint64_t* new_n_length)
{
(*new_n_length) = 0;
///这种就是0，new_n_list也不需要有数据
if (n_length == 0)
{
return;
}
char* x_string = NULL;
char* y_string = NULL;
long long x_offset = (long long)(n_end_pos) - WINDOW + 1;
if (x_offset < 0)
{
x_offset = 0;
}
///x_length是x上待验证区间的实际长度
///x_offset是x上待验证区间的起始位置
///如果是向前取待验证区间，那么x_length至少大于等于40 (compressed k-mer长度，不是实际长度)
long long x_length = n_end_pos - x_offset + 1;
///这种也不要过滤了，直接把old_n_list全部赋值过来就好了
if (x_length < (WINDOW/2))
{
(*new_n_length) = n_length;
memcpy(new_n_list, old_n_list, sizeof(k_mer_pos)*n_length);
return;
}
long long Window_Len = x_length + (THRESHOLD << 1);
x_string = g_read->seq + x_offset;
long long y_offset;
long long y_length;
long long i = 0;
long long y_read_length;
int end_site;
unsigned int error;
int groupLen = 0;
int return_sites[GROUP_SIZE];
unsigned int return_sites_error[GROUP_SIZE];
uint64_t readID[GROUP_SIZE];
uint64_t offset[GROUP_SIZE];
for (i = 0; i < n_length; i++)
{
////old_n_list[i].offset是y上k-mer的结束位置
///n_end_pos是x上k-mer的结束位置
///x_length是x上区间长度
y_offset = (long long)(old_n_list[i].offset) - x_length;
///这种情况下弃疗
if (y_offset < 0)
{
new_n_list[(*new_n_length)].readID = old_n_list[i].readID;
new_n_list[(*new_n_length)].offset = old_n_list[i].offset;
(*new_n_length)++;
continue;
}
y_offset = y_offset - THRESHOLD;
///还能抢救
if (y_offset < 0)
{
y_offset = 0;
}
///y的read的总长度
y_read_length = Get_READ_LENGTH((*R_INF), old_n_list[i].readID);
y_length = MIN(Window_Len, y_read_length - y_offset);
///如果y的长度比x长度还小，那就直接弃疗了
if (y_length < x_length)
{
new_n_list[(*new_n_length)].readID = old_n_list[i].readID;
new_n_list[(*new_n_length)].offset = old_n_list[i].offset;
(*new_n_length)++;
continue;
}
if(y_length == Window_Len)
{
///y的方向都是0，因为索引里都是0
recover_UC_Read_sub_region(dumy->overlap_region_group[groupLen], y_offset, y_length, 0,
R_INF, old_n_list[i].readID);
readID[groupLen] = old_n_list[i].readID;
offset[groupLen] = old_n_list[i].offset;
groupLen++;
if (groupLen == GROUP_SIZE)
{
Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1],
dumy->overlap_region_group[2], dumy->overlap_region_group[3], y_length, x_string, x_length,
return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE);
groupLen = 0;
if (return_sites_error[0]!=(unsigned int)-1)
{
new_n_list[(*new_n_length)].readID = readID[0];
new_n_list[(*new_n_length)].offset = offset[0];
(*new_n_length)++;
}
if (return_sites_error[1]!=(unsigned int)-1)
{
new_n_list[(*new_n_length)].readID = readID[1];
new_n_list[(*new_n_length)].offset = offset[1];
(*new_n_length)++;
}
if (return_sites_error[2]!=(unsigned int)-1)
{
new_n_list[(*new_n_length)].readID = readID[2];
new_n_list[(*new_n_length)].offset = offset[2];
(*new_n_length)++;
}
if (return_sites_error[3]!=(unsigned int)-1)
{
new_n_list[(*new_n_length)].readID = readID[3];
new_n_list[(*new_n_length)].offset = offset[3];
(*new_n_length)++;
}
}
}
else
{
///y的方向都是0，因为索引里都是0
recover_UC_Read_sub_region(dumy->overlap_region, y_offset, y_length, 0, R_INF, old_n_list[i].readID);
y_string = dumy->overlap_region;
memset (y_string + y_length, 0, Window_Len - y_length);
end_site = Reserve_Banded_BPM(y_string, y_length, x_string, x_length, THRESHOLD, &error);
if (error!=(unsigned int)-1)
{
new_n_list[(*new_n_length)].readID = old_n_list[i].readID;
new_n_list[(*new_n_length)].offset = old_n_list[i].offset;
(*new_n_length)++;
}
}
}
if (groupLen == 1)
{
end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], Window_Len, x_string, x_length, THRESHOLD, &error);
if (error!=(unsigned int)-1)
{
new_n_list[(*new_n_length)].readID = readID[0];
new_n_list[(*new_n_length)].offset = offset[0];
(*new_n_length)++;
}
}
else
{
Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1],
dumy->overlap_region_group[2], dumy->overlap_region_group[3], Window_Len, x_string, x_length,
return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE);
for (i = 0; i < groupLen; i++)
{
if (return_sites_error[i]!=(unsigned int)-1)
{
new_n_list[(*new_n_length)].readID = readID[i];
new_n_list[(*new_n_length)].offset = offset[i];
(*new_n_length)++;
}
}
groupLen = 0;
}
///要排序....
if ((*new_n_length)>1)
{
qsort(new_n_list, (*new_n_length), sizeof(k_mer_pos), cmp_k_mer_pos);
}
}
/**********************for prefilter************************ */
void destory_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list)
{
long long i = 0;
for (i = 0; i < list->size; i++)
{
if (list->list[i].size != 0)
{
free(list->list[i].list);
}
}
free(list->list);
}
void append_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length,
uint64_t n_end_pos, uint8_t n_direction, UC_Read* g_read, All_reads* R_INF, Correct_dumy* dumy)
{
if (list->length + 1 > list->size)
{
list->size = list->size * 2;
list->list = (k_mer_pos_list*)realloc(list->list, sizeof(k_mer_pos_list)*list->size);
///新分配空间要初始化
memset(list->list + (list->size/2), 0, sizeof(k_mer_pos_list)*(list->size/2));
}
if (list->list[list->length].size < n_length)
{
list->list[list->length].size = n_length;
list->list[list->length].list = (k_mer_pos*)realloc(list->list[list->length].list,
sizeof(k_mer_pos)*list->list[list->length].size);
}
///list->list[list->length].list = n_list;
///memcpy(list->list[list->length].list, n_list, sizeof(k_mer_pos)*n_length);
pre_filter_by_nearby(list->list[list->length].list, n_list, n_length, n_end_pos,
g_read, R_INF, dumy, &n_length);
if(n_length > 0)
{
list->list[list->length].length = n_length;
list->list[list->length].direction = n_direction;
list->list[list->length].end_pos = n_end_pos;
list->length++;
}
}
/**********************for prefilter************************ */