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

NOT WORKING!!! backup only

Browse Source
This commit is contained in:
Heng Li
2020-03-24 00:44:56 -04:00
parent e6ef2fb56f
commit 59a9d62df9
8 changed files with 131 additions and 714 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+6 -72
View File
@@ -225,7 +225,6 @@ void* Build_hash_table(void* arg)
int avalible_k = 0; int avalible_k = 0;
while (file_flag != 0) while (file_flag != 0)
{ {
@@ -802,13 +801,11 @@ overlap_region_alloc* overlap_list, int flag)
} }
} }
return available_overlaps; return available_overlaps;
} }
void get_new_candidates(long long readID, UC_Read* g_read, overlap_region_alloc* overlap_list, k_mer_pos_list_alloc* array_list, void get_new_candidates(long long readID, UC_Read* g_read, overlap_region_alloc* overlap_list, k_mer_pos_list_alloc* array_list,
HeapSq* heap, Candidates_list* l, double band_width_threshold, int keep_whole_chain) Candidates_list* l, double band_width_threshold, int keep_whole_chain)
{ {
HPC_seq HPC_read; HPC_seq HPC_read;
Hash_code k_code; Hash_code k_code;
@@ -819,7 +816,6 @@ HeapSq* heap, Candidates_list* l, double band_width_threshold, int keep_whole_ch
uint64_t list_length; uint64_t list_length;
uint64_t sub_ID; uint64_t sub_ID;
clear_Heap(heap);
clear_Candidates_list(l); clear_Candidates_list(l);
clear_k_mer_pos_list_alloc(array_list); clear_k_mer_pos_list_alloc(array_list);
@@ -827,18 +823,15 @@ HeapSq* heap, Candidates_list* l, double band_width_threshold, int keep_whole_ch
recover_UC_Read(g_read, &R_INF, readID); recover_UC_Read(g_read, &R_INF, readID);
///forward strand
init_HPC_seq(&HPC_read, g_read->seq, g_read->length); init_HPC_seq(&HPC_read, g_read->seq, g_read->length);
init_Hash_code(&k_code); init_Hash_code(&k_code);
avalible_k = 0; avalible_k = 0;
while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6) while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
{ {
if(code < 4) if(code < 4)
{ {
k_mer_append(&k_code,code, asm_opt.k_mer_length); k_mer_append(&k_code, code, asm_opt.k_mer_length);
avalible_k++; avalible_k++;
if (avalible_k >= asm_opt.k_mer_length) if (avalible_k >= asm_opt.k_mer_length)
{ {
@@ -857,38 +850,8 @@ HeapSq* heap, Candidates_list* l, double band_width_threshold, int keep_whole_ch
} }
} }
// BIG CHANGES WILL GO HERE!!!
///reverse complement strand
reverse_complement(g_read->seq, g_read->length);
init_HPC_seq(&HPC_read, g_read->seq, g_read->length);
init_Hash_code(&k_code);
avalible_k = 0;
while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
{
if(code < 4)
{
k_mer_append(&k_code,code, asm_opt.k_mer_length);
avalible_k++;
if (avalible_k >= asm_opt.k_mer_length)
{
list_length = locate_Total_Pos_Table(&PCB, &k_code, &list, asm_opt.k_mer_length, &sub_ID);
if (list_length != 0)
{
append_k_mer_pos_list_alloc(array_list, list, list_length, end_pos, 1);
}
}
}
else
{
avalible_k = 0;
init_Hash_code(&k_code);
}
}
merge_k_mer_pos_list_alloc_heap_sort(array_list, l, heap);
calculate_overlap_region_by_chaining(l, overlap_list, readID, g_read->length, &R_INF, calculate_overlap_region_by_chaining(l, overlap_list, readID, g_read->length, &R_INF,
band_width_threshold, keep_whole_chain); band_width_threshold, keep_whole_chain);
} }
@@ -922,10 +885,6 @@ void* Overlap_calculate_heap_merge(void* arg)
overlap_region_alloc overlap_list; overlap_region_alloc overlap_list;
init_overlap_region_alloc(&overlap_list); init_overlap_region_alloc(&overlap_list);
HeapSq heap;
Init_Heap(&heap);
Correct_dumy correct; Correct_dumy correct;
init_Correct_dumy(&correct); init_Correct_dumy(&correct);
@@ -946,7 +905,7 @@ void* Overlap_calculate_heap_merge(void* arg)
for (i = thr_ID; i < (long long)R_INF.total_reads; i = i + asm_opt.thread_num) for (i = thr_ID; i < (long long)R_INF.total_reads; i = i + asm_opt.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, 1); get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.02, 1);
clear_Cigar_record(&current_cigar); clear_Cigar_record(&current_cigar);
clear_Round2_alignment(&second_round); clear_Round2_alignment(&second_round);
@@ -982,7 +941,6 @@ void* Overlap_calculate_heap_merge(void* arg)
destory_buffer_sub_block(&current_sub_buffer); destory_buffer_sub_block(&current_sub_buffer);
destory_Candidates_list(&l); destory_Candidates_list(&l);
destory_overlap_region_alloc(&overlap_list); destory_overlap_region_alloc(&overlap_list);
destory_Heap(&heap);
destory_k_mer_pos_list_alloc(&array_list); destory_k_mer_pos_list_alloc(&array_list);
destory_Graph(&POA_Graph); destory_Graph(&POA_Graph);
destory_Graph(&DAGCon); destory_Graph(&DAGCon);
@@ -1040,9 +998,6 @@ void* Output_related_reads(void* arg)
overlap_region_alloc overlap_list; overlap_region_alloc overlap_list;
init_overlap_region_alloc(&overlap_list); init_overlap_region_alloc(&overlap_list);
HeapSq heap;
Init_Heap(&heap);
Correct_dumy correct; Correct_dumy correct;
init_Correct_dumy(&correct); init_Correct_dumy(&correct);
@@ -1072,7 +1027,7 @@ void* Output_related_reads(void* arg)
memcmp(asm_opt.required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) memcmp(asm_opt.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, 1); get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.02, 1);
fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i),
Get_NAME((R_INF), i)); Get_NAME((R_INF), i));
@@ -1097,7 +1052,6 @@ void* Output_related_reads(void* arg)
destory_buffer_sub_block(&current_sub_buffer); destory_buffer_sub_block(&current_sub_buffer);
destory_Candidates_list(&l); destory_Candidates_list(&l);
destory_overlap_region_alloc(&overlap_list); destory_overlap_region_alloc(&overlap_list);
destory_Heap(&heap);
destory_k_mer_pos_list_alloc(&array_list); destory_k_mer_pos_list_alloc(&array_list);
destory_Graph(&POA_Graph); destory_Graph(&POA_Graph);
destory_Graph(&DAGCon); destory_Graph(&DAGCon);
@@ -1626,9 +1580,6 @@ void* Final_overlap_calculate_heap_merge(void* arg)
overlap_region_alloc overlap_list; overlap_region_alloc overlap_list;
init_overlap_region_alloc(&overlap_list); init_overlap_region_alloc(&overlap_list);
HeapSq heap;
Init_Heap(&heap);
Cigar_record_alloc cigarline; Cigar_record_alloc cigarline;
init_Cigar_record_alloc(&cigarline); init_Cigar_record_alloc(&cigarline);
@@ -1636,15 +1587,11 @@ void* Final_overlap_calculate_heap_merge(void* arg)
memset(c2n, 4, 256); memset(c2n, 4, 256);
c2n[(uint8_t)'A'] = c2n[(uint8_t)'a'] = 0; c2n[(uint8_t)'C'] = c2n[(uint8_t)'c'] = 1; c2n[(uint8_t)'A'] = c2n[(uint8_t)'a'] = 0; c2n[(uint8_t)'C'] = c2n[(uint8_t)'c'] = 1;
c2n[(uint8_t)'G'] = c2n[(uint8_t)'g'] = 2; c2n[(uint8_t)'T'] = c2n[(uint8_t)'t'] = 3; // build the encoding table c2n[(uint8_t)'G'] = c2n[(uint8_t)'g'] = 2; c2n[(uint8_t)'T'] = c2n[(uint8_t)'t'] = 3; // build the encoding table
for (i = thr_ID; i < R_INF.total_reads; i = i + asm_opt.thread_num) for (i = thr_ID; i < R_INF.total_reads; i = i + asm_opt.thread_num)
{ {
get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.001, 0); get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
/** /**
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,
@@ -1677,7 +1624,6 @@ void* Final_overlap_calculate_heap_merge(void* arg)
destory_Candidates_list(&l); destory_Candidates_list(&l);
destory_overlap_region_alloc(&overlap_list); destory_overlap_region_alloc(&overlap_list);
destory_Heap(&heap);
destory_k_mer_pos_list_alloc(&array_list); destory_k_mer_pos_list_alloc(&array_list);
destory_UC_Read(&g_read); destory_UC_Read(&g_read);
destory_UC_Read(&overlap_read); destory_UC_Read(&overlap_read);
@@ -1946,15 +1892,3 @@ void Correct_Reads(int last_round)
Correct_Reads(last_round - 1); Correct_Reads(last_round - 1);
} }
CommandLines.cpp
+1 -1
View File
@@ -63,7 +63,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->pat_index = NULL; asm_opt->pat_index = NULL;
asm_opt->mat_index = NULL; asm_opt->mat_index = NULL;
asm_opt->thread_num = 1; asm_opt->thread_num = 1;
asm_opt->k_mer_length = 40; asm_opt->k_mer_length = 39;
asm_opt->k_mer_min_freq = 3; asm_opt->k_mer_min_freq = 3;
asm_opt->k_mer_max_freq = 66; asm_opt->k_mer_max_freq = 66;
asm_opt->load_index_from_disk = 1; asm_opt->load_index_from_disk = 1;
Hash_Table.cpp
+29 -502
View File
@@ -18,141 +18,6 @@ void overlap_region_sort_y_id(overlap_region *a, long long n)
radix_sort_overlap_region_sort(a, a + n); radix_sort_overlap_region_sort(a, a + n);
} }
void Init_Heap(HeapSq* HBT)
{
HBT->MaxSize = 1000;
HBT->heap = (ElemType*)malloc(HBT->MaxSize*sizeof(ElemType));
HBT->index_i = (uint64_t*)malloc(HBT->MaxSize*sizeof(uint64_t));
HBT->len = 0;
}
void destory_Heap(HeapSq* HBT)
{
free(HBT->heap);
free(HBT->index_i);
}
void clear_Heap(HeapSq* HBT)
{
HBT->len = 0;
}
inline int cmp_ElemType(ElemType* x, ElemType* y)
{
if (x->node.readID < y->node.readID)
{
return 1;
}
else if (x->node.readID > y->node.readID)
{
return 2;
}
else
{ if (x->node.strand < y->node.strand)
{
return 1;
}
else if (x->node.strand > y->node.strand)
{
return 2;
}
else
{
if(x->node.offset < y->node.offset)
{
return 1;
}
else if(x->node.offset > y->node.offset)
{
return 2;
}
else
{
if (x->node.self_offset < y->node.self_offset)
{
return 1;
}
else if (x->node.self_offset > y->node.self_offset)
{
return 2;
}
else ///if both r_pos_x and self_offset are equal, offset must be equal
{
return 0;
}
}
}
}
}
inline void Insert_Heap(HeapSq* HBT, ElemType* x)
{
long long i, j;
if (HBT->len == HBT->MaxSize)
{
HBT->MaxSize = 2*HBT->MaxSize;
HBT->heap = (ElemType*)realloc(HBT->heap, HBT->MaxSize*sizeof(ElemType));
HBT->index_i = (uint64_t*)realloc(HBT->index_i,HBT->MaxSize*sizeof(uint64_t));
}
HBT->heap[HBT->len] = *x; //add element to tail
HBT->len++;
i = HBT->len - 1;
while (i != 0)
{
j = (i - 1) / 2;
///if (x >= HBT->heap[j])
///1: x<y
if (cmp_ElemType(x, &HBT->heap[j])!=1)
break;
HBT->heap[i] = HBT->heap[j];
i = j;
}
HBT->heap[i] = *x;
}
inline int DeleteHeap(HeapSq* HBT, ElemType* get)
{
ElemType temp, x;
int i, j;
if (HBT->len == 0)
{
return 0;
}
temp = HBT->heap[0];
HBT->len--;
if (HBT->len == 0)
{
*get = temp;
return 2;
}
x = HBT->heap[HBT->len];
i = 0;
j = 2 * i + 1;
while (j <= HBT->len - 1)
{
///if (j < HBT->len - 1 && HBT->heap[j] > HBT->heap[j+1])
if (j < HBT->len - 1 && cmp_ElemType(&HBT->heap[j], &HBT->heap[j + 1]) == 2)
j++;
///if (x <= HBT->heap[j])
if (cmp_ElemType(&x, &HBT->heap[j])!=2)
break;
HBT->heap[i] = HBT->heap[j];
i = j;
j = 2 * i + 1;
}
HBT->heap[i] = x;
*get = temp;
return 1;
}
void init_overlap_region_alloc(overlap_region_alloc* list) void init_overlap_region_alloc(overlap_region_alloc* list)
{ {
list->size = 1000; list->size = 1000;
@@ -219,10 +84,8 @@ int get_fake_gap_shift(Fake_Cigar* x, int index)
return result; return result;
} }
int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp,
All_reads* R_INF, int add_beg_end) All_reads* R_INF, int add_beg_end)
{ {
if (list->length + 1 > list->size) if (list->length + 1 > list->size)
@@ -264,8 +127,6 @@ All_reads* R_INF, int add_beg_end)
long long x_right_length = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_e - 1; 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; 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) if(x_right_length <= y_right_length)
{ {
tmp->x_pos_e = Get_READ_LENGTH((*R_INF), tmp->x_id) - 1; tmp->x_pos_e = Get_READ_LENGTH((*R_INF), tmp->x_id) - 1;
@@ -276,7 +137,6 @@ All_reads* R_INF, int add_beg_end)
tmp->x_pos_e = tmp->x_pos_e + y_right_length; tmp->x_pos_e = tmp->x_pos_e + y_right_length;
tmp->y_pos_e = Get_READ_LENGTH((*R_INF), tmp->y_id) - 1; tmp->y_pos_e = Get_READ_LENGTH((*R_INF), tmp->y_id) - 1;
} }
if (tmp->x_pos_strand == 1) if (tmp->x_pos_strand == 1)
{ {
@@ -290,9 +150,6 @@ All_reads* R_INF, int add_beg_end)
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_s = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_e - 1;
list->list[list->length].y_pos_strand = 1; list->list[list->length].y_pos_strand = 1;
resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2)); resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2));
if(add_beg_end == 1) if(add_beg_end == 1)
{ {
@@ -388,7 +245,6 @@ All_reads* R_INF, int add_beg_end)
return 1; return 1;
} }
void append_overlap_region_alloc_debug(overlap_region_alloc* list, overlap_region* tmp) void append_overlap_region_alloc_debug(overlap_region_alloc* list, overlap_region* tmp)
{ {
@@ -414,7 +270,6 @@ void append_overlap_region_alloc_debug(overlap_region_alloc* list, overlap_regio
list->length++; list->length++;
} }
int cmp_by_x_pos_s(const void * a, const void * b) 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) if ((*(overlap_region*)a).x_pos_s > (*(overlap_region*)b).x_pos_s)
@@ -474,7 +329,6 @@ int cmp_by_x_pos_e(const void * a, const void * b)
} }
} }
void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp) void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp)
{ {
long long i, j, current_j; long long i, j, current_j;
@@ -493,7 +347,7 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp)
if(j != -1) if(j != -1)
{ {
distance_self_pos = a[current_j].self_offset - a[j].self_offset; distance_self_pos = a[current_j].self_offset - a[j].self_offset;
distance_pos = a[current_j].offset - a[j].offset; distance_pos = ha_hit_get_offset(&a[current_j]) - ha_hit_get_offset(&a[j]);
distance_gap = distance_pos > distance_self_pos? distance_pos - distance_self_pos : distance_self_pos - distance_pos; distance_gap = distance_pos > distance_self_pos? distance_pos - distance_self_pos : distance_self_pos - distance_pos;
indels += distance_gap; indels += distance_gap;
@@ -517,7 +371,7 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp)
} }
long long get_chainLen(long long x_beg, long long x_end, long long xLen, long long get_chainLen(long long x_beg, long long x_end, long long xLen,
long long y_beg, long long y_end, long long yLen) long long y_beg, long long y_end, long long yLen)
{ {
if(x_beg <= y_beg) if(x_beg <= y_beg)
{ {
@@ -551,7 +405,7 @@ long long y_beg, long long y_end, long long yLen)
///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)
{ {
long long i, j; long long i, j;
long long self_pos, pos, max_j, max_i, max_score, score, n_skip; long long self_pos, pos, max_j, max_i, max_score, score, n_skip;
@@ -567,7 +421,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen)
// fill the score and backtrack arrays // fill the score and backtrack arrays
for (i = 0; i < a_n; ++i) for (i = 0; i < a_n; ++i)
{ {
pos = a[i].offset; pos = ha_hit_get_offset(&a[i]);
self_pos = a[i].self_offset; self_pos = a[i].self_offset;
max_j = -1; max_j = -1;
max_score = min_score; max_score = min_score;
@@ -579,7 +433,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen)
///may have a pre-cut condition for j ///may have a pre-cut condition for j
for (j = i - 1; j >= 0; --j) for (j = i - 1; j >= 0; --j)
{ {
distance_pos = pos - a[j].offset; distance_pos = pos - ha_hit_get_offset(&a[j]);
distance_self_pos = self_pos - a[j].self_offset; distance_self_pos = self_pos - a[j].self_offset;
///a has been sorted by a[].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 ///note for a, we do not have any two elements that have both equal offsets and self_offsets
@@ -636,12 +490,8 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen)
dp->self_length[i] = max_self_length; dp->self_length[i] = max_self_length;
} }
///debug_chain(a, a_n, dp); ///debug_chain(a, a_n, dp);
max_score = -1; max_score = -1;
max_i = -1; max_i = -1;
long long mini_xLen = x_readLen * 2 + 2, tmp_xLen; long long mini_xLen = x_readLen * 2 + 2, tmp_xLen;
@@ -652,12 +502,12 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen)
max_score = dp->score[i]; max_score = dp->score[i];
max_i = i; max_i = i;
mini_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen, mini_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen,
a[i].offset, a[i].offset, y_readLen); ha_hit_get_offset(&a[i]), ha_hit_get_offset(&a[i]), y_readLen);
} }
else if(dp->score[i] == max_score) else if(dp->score[i] == max_score)
{ {
tmp_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen, tmp_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen,
a[i].offset, a[i].offset, y_readLen); ha_hit_get_offset(&a[i]), ha_hit_get_offset(&a[i]), y_readLen);
if(tmp_xLen < mini_xLen) if(tmp_xLen < mini_xLen)
{ {
@@ -670,17 +520,16 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen)
} }
clear_fake_cigar(&(result->f_cigar)); clear_fake_cigar(&(result->f_cigar));
///not a has been sorted by offset, that means has been sorted by query offset ///not a has been sorted by offset, that means has been sorted by query offset
i = max_i; i = max_i;
result->x_pos_e = a[i].self_offset; result->x_pos_e = a[i].self_offset;
result->y_pos_e = a[i].offset; result->y_pos_e = ha_hit_get_offset(&a[i]);
result->shared_seed = max_score; result->shared_seed = max_score;
result->overlapLen = mini_xLen; result->overlapLen = mini_xLen;
distance_self_pos = result->x_pos_e - a[i].self_offset; distance_self_pos = result->x_pos_e - a[i].self_offset;
distance_pos = result->y_pos_e - a[i].offset; distance_pos = result->y_pos_e - ha_hit_get_offset(&a[i]);
long long pre_distance_gap = distance_pos - distance_self_pos; long long pre_distance_gap = distance_pos - distance_self_pos;
///record first site ///record first site
///the length of f_cigar should be at least 1 ///the length of f_cigar should be at least 1
@@ -692,7 +541,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen)
while (i >= 0) while (i >= 0)
{ {
distance_self_pos = result->x_pos_e - a[i].self_offset; distance_self_pos = result->x_pos_e - a[i].self_offset;
distance_pos = result->y_pos_e - a[i].offset; distance_pos = result->y_pos_e - ha_hit_get_offset(&a[i]);
distance_gap = distance_pos - distance_self_pos; distance_gap = distance_pos - distance_self_pos;
if(distance_gap != pre_distance_gap) if(distance_gap != pre_distance_gap)
{ {
@@ -703,7 +552,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen)
chainLen++; chainLen++;
result->x_pos_s = a[i].self_offset; result->x_pos_s = a[i].self_offset;
result->y_pos_s = a[i].offset; result->y_pos_s = ha_hit_get_offset(&a[i]);
i = dp->pre[i]; i = dp->pre[i];
} }
} }
@@ -713,7 +562,7 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen)
while (i >= 0) while (i >= 0)
{ {
distance_self_pos = result->x_pos_e - a[i].self_offset; distance_self_pos = result->x_pos_e - a[i].self_offset;
distance_pos = result->y_pos_e - a[i].offset; distance_pos = result->y_pos_e - ha_hit_get_offset(&a[i]);
distance_gap = distance_pos - distance_self_pos; distance_gap = distance_pos - distance_self_pos;
if(distance_gap == pre_distance_gap) if(distance_gap == pre_distance_gap)
{ {
@@ -728,16 +577,14 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen)
chainLen++; chainLen++;
result->x_pos_s = a[i].self_offset; result->x_pos_s = a[i].self_offset;
result->y_pos_s = a[i].offset; result->y_pos_s = ha_hit_get_offset(&a[i]);
i = dp->pre[i]; i = dp->pre[i];
} }
} }
} }
void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list,
uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end) uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end)
{ {
overlap_region tmp_region; overlap_region tmp_region;
long long i = 0; long long i = 0;
@@ -757,8 +604,8 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh
i = 0; i = 0;
while (i < candidates->length) while (i < candidates->length)
{ {
current_ID = candidates->list[i].readID; current_ID = ha_hit_get_readID(&candidates->list[i]);
current_stand = candidates->list[i].strand; current_stand = ha_hit_get_rev(&candidates->list[i]);
///reference read ///reference read
tmp_region.x_id = readID; tmp_region.x_id = readID;
@@ -776,9 +623,9 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh
while (i < candidates->length while (i < candidates->length
&& &&
current_ID == candidates->list[i].readID current_ID == ha_hit_get_readID(&candidates->list[i])
&& &&
current_stand == candidates->list[i].strand) current_stand == ha_hit_get_rev(&candidates->list[i]))
{ {
sub_region_end = i; sub_region_end = i;
i++; i++;
@@ -811,10 +658,8 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh
qsort(overlap_list->list, overlap_list->length, sizeof(overlap_region), cmp_by_x_pos_s); qsort(overlap_list->list, overlap_list->length, sizeof(overlap_region), cmp_by_x_pos_s);
} }
void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error, void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error,
int extra_begin, int extra_end, int error_threshold) int extra_begin, int extra_end, int error_threshold)
{ {
long long length = region->x_pos_e - region->x_pos_s + 1; long long length = region->x_pos_e - region->x_pos_s + 1;
@@ -843,8 +688,6 @@ int extra_begin, int extra_end, int error_threshold)
region->w_list_length++; region->w_list_length++;
} }
void init_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) void init_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list)
{ {
list->size = 1000; list->size = 1000;
@@ -864,9 +707,8 @@ void destory_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list)
} }
void append_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length, void append_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length,
uint64_t n_end_pos, uint8_t n_direction) uint64_t n_end_pos, uint8_t n_direction)
{ {
if (list->length + 1 > list->size) if (list->length + 1 > list->size)
{ {
list->size = list->size * 2; list->size = list->size * 2;
@@ -881,9 +723,6 @@ uint64_t n_end_pos, uint8_t n_direction)
list->length++; list->length++;
} }
int cmp_k_mer_pos_list(const void * a, const void * b) int cmp_k_mer_pos_list(const void * a, const void * b)
{ {
if ((*(k_mer_pos_list*)a).length > (*(k_mer_pos_list*)b).length) if ((*(k_mer_pos_list*)a).length > (*(k_mer_pos_list*)b).length)
@@ -898,18 +737,8 @@ int cmp_k_mer_pos_list(const void * a, const void * b)
{ {
return 0; return 0;
} }
} }
inline void append_pos_to_Candidates_list(Candidates_list* candidates, ElemType* x)
{
candidates->list[candidates->length] = x->node;
candidates->length++;
}
void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, uint64_t strand) void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, uint64_t strand)
{ {
uint64_t i; uint64_t i;
@@ -920,13 +749,13 @@ void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n
for (; j < candidates->length; j++) for (; j < candidates->length; j++)
{ {
if ( if (
n_list[i].offset == (uint64_t)candidates->list[j].offset n_list[i].offset == (uint64_t)ha_hit_get_offset(&candidates->list[j])
&& &&
n_list[i].readID == candidates->list[j].readID n_list[i].readID == ha_hit_get_readID(&candidates->list[j])
&& &&
end_pos == (uint64_t)candidates->list[j].self_offset end_pos == (uint64_t)candidates->list[j].self_offset
&& &&
strand == candidates->list[j].strand strand == ha_hit_get_rev(&candidates->list[j])
) )
{ {
break; break;
@@ -938,82 +767,8 @@ void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n
fprintf(stderr, "ERROR 4\n"); fprintf(stderr, "ERROR 4\n");
} }
} }
} }
void merge_k_mer_pos_list_alloc_heap_sort(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT)
{
clear_Heap(HBT);
uint64_t total_length = 0;
uint64_t i;
ElemType x, y;
//add the first element of each list to stack
for (i = 0; i < list->length; i++)
{
x.ID = i;
x.node.offset = list->list[i].list[0].offset;
x.node.readID = list->list[i].list[0].readID;
x.node.self_offset = list->list[i].end_pos;
x.node.strand = list->list[i].direction;
Insert_Heap(HBT, &x);
HBT->index_i[i] = 1;
total_length = total_length + list->list[i].length;
}
candidates->length = 0;
if(total_length > (uint64_t)candidates->size)
{
candidates->size = total_length;
candidates->list = (k_mer_hit*)realloc(candidates->list, sizeof(k_mer_hit)*candidates->size);
candidates->tmp = (k_mer_hit*)realloc(candidates->tmp, sizeof(k_mer_hit)*candidates->size);
}
uint64_t ID;
int flag;
///while (flag = DeleteHeap(HBT, &x))
while ((flag = DeleteHeap(HBT, &x)))
{
append_pos_to_Candidates_list(candidates, &x);
i = HBT->index_i[x.ID];
ID = x.ID;
if (flag == 2)
{
for (; i < list->list[x.ID].length; i++)
{
y.ID = ID;
y.node.offset = list->list[x.ID].list[i].offset;
y.node.readID = list->list[x.ID].list[i].readID;
y.node.self_offset = list->list[x.ID].end_pos;
y.node.strand = list->list[x.ID].direction;
append_pos_to_Candidates_list(candidates, &y);
}
break;
}
if (i < list->list[x.ID].length)
{
y.ID = ID;
y.node.offset = list->list[x.ID].list[i].offset;
y.node.readID = list->list[x.ID].list[i].readID;
y.node.self_offset = list->list[x.ID].end_pos;
y.node.strand = list->list[x.ID].direction;
Insert_Heap(HBT, &y);
HBT->index_i[ID]++;
}
}
}
void init_Count_Table(Count_Table** table) void init_Count_Table(Count_Table** table)
{ {
*table = ha_ct_init(); *table = ha_ct_init();
@@ -1233,127 +988,6 @@ void insert_H_peaks(H_peaks* h, long long index, long long value)
h->list[index] += value; h->list[index] += value;
} }
inline void RC_Hash_code(Hash_code* code, Hash_code* rc_code, int k)
{
rc_code->x[0] = 0;
rc_code->x[1] = 0;
int i;
for (i = 0; i < k; i++)
{
rc_code->x[0] = rc_code->x[0] << 1;
rc_code->x[1] = rc_code->x[1] << 1;
rc_code->x[0] |= (((uint64_t)((code->x[0] >> i) & 1))^((uint64_t)1));
rc_code->x[1] |= (((uint64_t)((code->x[1] >> i) & 1))^((uint64_t)1));
}
}
void get_peak_debug(Total_Count_Table* TCB, long long* min, long long* max)
{
int i;
Count_Table* h;
khint_t k;
long long c_count;
H_peaks LH;
LH.list = NULL;
LH.length = 0;
uint64_t sub_ID;
uint64_t sub_key;
Hash_code code, rc_code, debug_code;
char str[100];
char rc_str[100];
for (i = 0; i < TCB->size; i++)
{
h = TCB->sub_h[i];
for (k = 0; k != kh_end(h); ++k)
{
if (kh_exist(h, k)) // test if a bucket contains data
{
sub_ID = i;
sub_key = kh_key(h, k);
recover_hash_code(sub_ID, sub_key, &code, TCB->suffix_mode,
TCB->suffix_bits, asm_opt.k_mer_length);
RC_Hash_code(&code, &rc_code, asm_opt.k_mer_length);
RC_Hash_code(&rc_code, &debug_code, asm_opt.k_mer_length);
if(code.x[0] != debug_code.x[0] || code.x[1] != debug_code.x[1])
{
fprintf(stderr, "error\n");
}
Hashcode_to_string(&code, str, asm_opt.k_mer_length);
Hashcode_to_string(&rc_code, rc_str, asm_opt.k_mer_length);
reverse_complement(str, asm_opt.k_mer_length);
if(memcmp(str, rc_str, asm_opt.k_mer_length) != 0)
{
fprintf(stderr, "error\n");
int j;
for (j = 0; j < asm_opt.k_mer_length; j++)
{
fprintf(stderr, "%c",str[j]);
}
fprintf(stderr, "\n");
for (j = 0; j < asm_opt.k_mer_length; j++)
{
fprintf(stderr, "%c",rc_str[j]);
}
fprintf(stderr, "\n");
}
///get_Total_Count_Table(&TCB, &k_code, k_mer_length);
c_count = kh_val(h, k);
if(get_Total_Count_Table(TCB, &code, asm_opt.k_mer_length) != c_count)
{
fprintf(stderr, "error\n");
}
insert_H_peaks(&LH, c_count, c_count);
}
}
}
(*max) = -1;
(*min) = -1;
long long max_value = -1;
for (i = 0; i < (long long)LH.length; i++)
{
if(LH.list[i] >= max_value)
{
max_value = LH.list[i];
(*max) = i;
}
}
long long min_value = max_value;
for (i = 0; i < (long long)LH.length; i++)
{
if(LH.list[i] < min_value && LH.list[i] != 0)
{
min_value = LH.list[i];
(*min) = i;
}
}
for (i = 0; i < (long long)LH.length; i++)
{
///fprintf(stderr, "%d, %d\n", i, LH.list[i]);
fprintf(stderr, "%lld\n", LH.list[i]);
}
free(LH.list);
}
///1: a > b; -1: a < b; 0: a=b ///1: a > b; -1: a < b; 0: a=b
int cmp_Hash_code(Hash_code* a, Hash_code* b) int cmp_Hash_code(Hash_code* a, Hash_code* b)
{ {
@@ -1378,46 +1012,10 @@ int cmp_Hash_code(Hash_code* a, Hash_code* b)
return 0; return 0;
} }
int get_total_freq(Total_Count_Table* TCB, uint64_t sub_ID, uint64_t sub_key, long long* T_count) void get_total_freq(Total_Count_Table* TCB, uint64_t sub_ID, uint64_t sub_key, long long* T_count)
{ {
Hash_code code, rc_code; khint_t t = ha_ct_get(TCB->sub_h[sub_ID], sub_key);
long long count, rc_count; *T_count = kh_val(TCB->sub_h[sub_ID], t);
recover_hash_code(sub_ID, sub_key, &code, TCB->suffix_mode,
TCB->suffix_bits, asm_opt.k_mer_length);
RC_Hash_code(&code, &rc_code, asm_opt.k_mer_length);
count = get_Total_Count_Table(TCB, &code, asm_opt.k_mer_length);
rc_count = get_Total_Count_Table(TCB, &rc_code, asm_opt.k_mer_length);
(*T_count) = count + rc_count;
if(count == 0)
{
return 0;
}///count > 0 && rc_count == 0
else if(rc_count == 0)
{
return 1;
}///count > 0 && rc_count > 0
else
{
int flag = cmp_Hash_code(&code, &rc_code);
///code > rc_code
if(flag > 0)
{
return 1;
}///code < rc_code
else if(flag < 0)
{
return 0;
}
else
{
(*T_count) = (*T_count)/2;
return 1;
}
}
} }
void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long* up_boundary) void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long* up_boundary)
@@ -1432,8 +1030,6 @@ void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long*
uint64_t sub_ID; uint64_t sub_ID;
uint64_t sub_key; uint64_t sub_key;
for (i = 0; i < TCB->size; i++) for (i = 0; i < TCB->size; i++)
{ {
h = TCB->sub_h[i]; h = TCB->sub_h[i];
@@ -1444,10 +1040,8 @@ void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long*
sub_ID = i; sub_ID = i;
sub_key = kh_key(h, k); sub_key = kh_key(h, k);
if(get_total_freq(TCB, sub_ID, sub_key, &count)==1) get_total_freq(TCB, sub_ID, sub_key, &count);
{ insert_H_peaks(&LH, count, count);
insert_H_peaks(&LH, count, count);
}
} }
} }
} }
@@ -1492,9 +1086,6 @@ void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long*
(*up_boundary) = (*max) * 10; (*up_boundary) = (*max) * 10;
} }
long long min_value = max_value; long long min_value = max_value;
//// seed with freq 1 is useless //// seed with freq 1 is useless
for (i = 2; i < (long long)LH.length && i < (*max); i++) for (i = 2; i < (long long)LH.length && i < (*max); i++)
@@ -1509,8 +1100,6 @@ void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long*
free(LH.list); free(LH.list);
} }
void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k_mer_min_freq, int k_mer_max_freq) void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k_mer_min_freq, int k_mer_max_freq)
{ {
int i; int i;
@@ -1571,7 +1160,6 @@ void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k
///kh_val(PCB->sub_h[sub_ID], t)++; ///kh_val(PCB->sub_h[sub_ID], t)++;
fprintf(stderr, "ERROR\n"); fprintf(stderr, "ERROR\n");
} }
PCB->useful_k_mer++; PCB->useful_k_mer++;
PCB->total_occ = PCB->total_occ + kh_val(h, k); PCB->total_occ = PCB->total_occ + kh_val(h, k);
@@ -1580,7 +1168,6 @@ void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k
} }
} }
// fprintf(stdout, "useful_k_mer: %lld\n",PCB->useful_k_mer); // fprintf(stdout, "useful_k_mer: %lld\n",PCB->useful_k_mer);
// fprintf(stdout, "total_occ: %lld\n",PCB->total_occ); // fprintf(stdout, "total_occ: %lld\n",PCB->total_occ);
@@ -1615,13 +1202,8 @@ void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k
PCB->pos = (k_mer_pos*)malloc(sizeof(k_mer_pos)*PCB->total_occ); PCB->pos = (k_mer_pos*)malloc(sizeof(k_mer_pos)*PCB->total_occ);
memset(PCB->pos, 0, sizeof(k_mer_pos)*PCB->total_occ); memset(PCB->pos, 0, sizeof(k_mer_pos)*PCB->total_occ);
} }
int cmp_k_mer_pos(const void * a, const void * b) int cmp_k_mer_pos(const void * a, const void * b)
{ {
if ((*(k_mer_pos*)a).readID != (*(k_mer_pos*)b).readID) if ((*(k_mer_pos*)a).readID != (*(k_mer_pos*)b).readID)
@@ -1657,7 +1239,6 @@ void clear_Chain_Data(Chain_Data* x)
x->length = 0; x->length = 0;
} }
void destory_Chain_Data(Chain_Data* x) void destory_Chain_Data(Chain_Data* x)
{ {
free(x->score); free(x->score);
@@ -1666,7 +1247,6 @@ void destory_Chain_Data(Chain_Data* x)
free(x->self_length); free(x->self_length);
} }
void resize_Chain_Data(Chain_Data* x, long long size) void resize_Chain_Data(Chain_Data* x, long long size)
{ {
if(size > x->size) if(size > x->size)
@@ -1679,23 +1259,18 @@ void resize_Chain_Data(Chain_Data* x, long long size)
} }
} }
void init_Candidates_list(Candidates_list* l) void init_Candidates_list(Candidates_list* l)
{ {
l->length = 0; l->length = 0;
l->size = 0; l->size = 0;
l->list = NULL; l->list = NULL;
l->tmp = NULL; l->tmp = NULL;
l->foward_pos = 0;
l->rc_pos = 0;
init_Chain_Data(&(l->chainDP)); init_Chain_Data(&(l->chainDP));
} }
void clear_Candidates_list(Candidates_list* l) void clear_Candidates_list(Candidates_list* l)
{ {
l->length = 0; l->length = 0;
l->foward_pos = 0;
l->rc_pos = 0;
clear_Chain_Data(&(l->chainDP)); clear_Chain_Data(&(l->chainDP));
} }
@@ -1706,52 +1281,6 @@ void destory_Candidates_list(Candidates_list* l)
destory_Chain_Data(&(l->chainDP)); destory_Chain_Data(&(l->chainDP));
} }
int cmp_candidates_list(const void * a, const void * b)
{
long long r_pos_a, r_pos_b;
if((*((k_mer_hit*)a)).strand != (*((k_mer_hit*)b)).strand)
{
return (*((k_mer_hit*)a)).strand > (*((k_mer_hit*)b)).strand ? 1: -1;
}
else
{
if((*((k_mer_hit*)a)).readID != (*((k_mer_hit*)b)).readID)
{
return (*((k_mer_hit*)a)).readID > (*((k_mer_hit*)b)).readID ? 1: -1;
}
else
{
r_pos_a = (*((k_mer_hit*)a)).offset - (*((k_mer_hit*)a)).self_offset;
r_pos_b = (*((k_mer_hit*)b)).offset - (*((k_mer_hit*)b)).self_offset;
if(r_pos_a != r_pos_b)
{
return r_pos_a > r_pos_b ? 1: -1;
}
else
{
if((*((k_mer_hit*)a)).self_offset != (*((k_mer_hit*)b)).self_offset)
{
return (*((k_mer_hit*)a)).self_offset > (*((k_mer_hit*)b)).self_offset ? 1: -1;
}
else
{
return 0;
}
}
}
}
}
void init_fake_cigar(Fake_Cigar* x) void init_fake_cigar(Fake_Cigar* x)
{ {
x->buffer = NULL; x->buffer = NULL;
@@ -1847,7 +1376,5 @@ void resize_window_list_alloc(window_list_alloc* x, long long size)
{ {
x->buffer[i].error = -1; x->buffer[i].error = -1;
} }
x->length = 0; x->length = 0;
} }
Hash_Table.h
+54 -94
View File
@@ -35,8 +35,7 @@ KHASHL_MAP_INIT(static inline, Pos_Table, ha_pt, uint64_t, uint64_t, kh_hash_dum
typedef struct typedef struct
{ {
volatile int lock; volatile int lock;
} Hash_table_spin_lock;
}Hash_table_spin_lock;
typedef struct typedef struct
{ {
@@ -72,12 +71,11 @@ typedef struct
uint64_t length; uint64_t length;
} k_mer_pos_list_alloc; } k_mer_pos_list_alloc;
typedef struct typedef struct
{ {
int C_L[CIGAR_MAX_LENGTH]; int C_L[CIGAR_MAX_LENGTH];
char C_C[CIGAR_MAX_LENGTH]; char C_C[CIGAR_MAX_LENGTH];
int length; int length;
} CIGAR; } CIGAR;
typedef struct typedef struct
@@ -100,15 +98,14 @@ typedef struct
window_list* buffer; window_list* buffer;
long long length; long long length;
long long size; long long size;
}window_list_alloc; } window_list_alloc;
typedef struct typedef struct
{ {
uint64_t* buffer; uint64_t* buffer;
uint64_t length; uint64_t length;
uint64_t size; uint64_t size;
}Fake_Cigar; } Fake_Cigar;
typedef struct typedef struct
{ {
@@ -140,7 +137,6 @@ typedef struct
window_list_alloc boundary_cigars; window_list_alloc boundary_cigars;
} overlap_region; } overlap_region;
typedef struct typedef struct
{ {
overlap_region* list; overlap_region* list;
@@ -152,29 +148,24 @@ typedef struct
typedef struct typedef struct
{ {
///uint64_t offset; uint64_t opos;
long long offset; uint32_t self_offset; // offset on the target read
///uint64_t self_offset;
long long self_offset;
uint64_t readID;
uint8_t strand;
} k_mer_hit; } k_mer_hit;
static inline uint32_t ha_hit_get_readID(const k_mer_hit *h)
typedef struct
{ {
k_mer_hit node; return h->opos >> 33;
uint64_t ID; }
} ElemType;
static inline uint32_t ha_hit_get_rev(const k_mer_hit *h)
typedef struct
{ {
ElemType* heap; return h->opos >> 32 & 1;
uint64_t* index_i; }
int len;
int MaxSize; static inline uint32_t ha_hit_get_offset(const k_mer_hit *h)
} HeapSq; {
return (uint32_t)h->opos;
}
typedef struct typedef struct
{ {
@@ -192,8 +183,6 @@ typedef struct
k_mer_hit* tmp; k_mer_hit* tmp;
long long length; long long length;
long long size; long long size;
uint64_t foward_pos;
uint64_t rc_pos;
Chain_Data chainDP; Chain_Data chainDP;
} Candidates_list; } Candidates_list;
@@ -212,24 +201,9 @@ typedef struct
uint64_t* k_mer_index; uint64_t* k_mer_index;
} Total_Pos_Table; } Total_Pos_Table;
inline uint64_t mod_d(uint64_t h_key, uint64_t low_key, uint64_t d)
{
uint64_t result = (h_key >> 32) % d;
result = ((result << 32) + (h_key & (uint64_t)0xffffffff)) % d;
result = ((result << 32) + (low_key >> 32)) % d;
result = ((result << 32) + (low_key & (uint64_t)0xffffffff)) % d;
return result;
}
////suffix_bits = 64 in default ////suffix_bits = 64 in default
inline int recover_hash_code(uint64_t sub_ID, uint64_t sub_key, Hash_code* code, inline int recover_hash_code(uint64_t sub_ID, uint64_t sub_key, Hash_code* code,
uint64_t suffix_mode, int suffix_bits, int k) uint64_t suffix_mode, int suffix_bits, int k) // FIXME: not working right now
{ {
uint64_t h_key, low_key; uint64_t h_key, low_key;
h_key = low_key = 0; h_key = low_key = 0;
@@ -249,31 +223,39 @@ uint64_t suffix_mode, int suffix_bits, int k)
return 1; return 1;
} }
///inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, Total_Count_Table* TCB, Hash_code* code, int k) static inline int ha_code2rev(const Hash_code *code)
inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t suffix_mode, int suffix_bits,
Hash_code* code, int k)
{ {
uint64_t h_key, low_key; return code->x[1] < code->x[3]? 0 : 1;
///k might be 64，so it is unsafe
///low_key = code->x[0] | (code->x[1] << k);
low_key = code->x[0] | (code->x[1] << SAFE_SHIFT(k));
//k cannot be 0, so this shift is safe
h_key = code->x[1] >> (64 - k);
if(mod_d(h_key, low_key, MODE_VALUE) > 3)
{
return 0;
}
uint64_t sub_ID = (low_key >> SAFE_SHIFT(suffix_bits)) | (h_key << (64 - suffix_bits));
uint64_t sub_key = (low_key & suffix_mode);
*get_sub_ID = sub_ID;
*get_sub_key = sub_key;
return 1;
} }
static inline int ha_get_sub_table_short(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t suffix_mode, int suffix_bits, Hash_code* code, int k)
{ // for k < 32
int j = ha_code2rev(code);
uint64_t y = code->x[j<<1|1] << k | code->x[j<<1|0];
y = yak_hash64(y, (1ULL<<(k+k)) - 1);
if (y % MODE_VALUE > 3) return 0;
*get_sub_ID = y >> suffix_bits;
*get_sub_key = y & suffix_mode;
return 1;
}
static inline int ha_get_sub_table_long(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t suffix_mode, int suffix_bits, Hash_code* code, int k)
{ // for k > 32
int j = ha_code2rev(code);
uint64_t y = code->x[j<<1|1] << k | code->x[j<<1|0];
y = yak_hash64_64(y);
if (y % MODE_VALUE > 3) return 0;
int s = 64 - k;
uint64_t z = code->x[j<<1|1] >> s ^ y << (s + s) >> (s + s);
*get_sub_ID = y >> suffix_bits | z << (64 - suffix_bits);
*get_sub_key = y & suffix_mode;
return 1;
}
inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t suffix_mode, int suffix_bits, Hash_code* code, int k)
{ // not really working for k<=32
return ha_get_sub_table_long(get_sub_ID, get_sub_key, suffix_mode, suffix_bits, code, k);
}
inline int insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) inline int insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k)
{ {
@@ -286,7 +268,6 @@ inline int insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int
khint_t t; khint_t t;
int absent; int absent;
while (__sync_lock_test_and_set(&TCB->sub_h_lock[sub_ID].lock, 1)) while (__sync_lock_test_and_set(&TCB->sub_h_lock[sub_ID].lock, 1))
{ {
while (TCB->sub_h_lock[sub_ID].lock); while (TCB->sub_h_lock[sub_ID].lock);
@@ -310,7 +291,6 @@ inline int insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int
inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k)
{ {
uint64_t sub_ID, sub_key; uint64_t sub_ID, sub_key;
if(!get_sub_table(&sub_ID, &sub_key, TCB->suffix_mode, TCB->suffix_bits, code, k)) if(!get_sub_table(&sub_ID, &sub_key, TCB->suffix_mode, TCB->suffix_bits, code, k))
{ {
@@ -319,7 +299,6 @@ inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k)
khint_t t; khint_t t;
///query hash table，key is k
t = ha_ct_get(TCB->sub_h[sub_ID], sub_key); t = ha_ct_get(TCB->sub_h[sub_ID], sub_key);
if (t != kh_end(TCB->sub_h[sub_ID])) if (t != kh_end(TCB->sub_h[sub_ID]))
@@ -330,15 +309,10 @@ inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k)
{ {
return 0; return 0;
} }
} }
inline uint64_t get_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k, uint64_t* r_sub_ID) inline uint64_t get_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k, uint64_t* r_sub_ID)
{ {
uint64_t sub_ID, sub_key; uint64_t sub_ID, sub_key;
if(!get_sub_table(&sub_ID, &sub_key, PCB->suffix_mode, PCB->suffix_bits, code, k)) if(!get_sub_table(&sub_ID, &sub_key, PCB->suffix_mode, PCB->suffix_bits, code, k))
{ {
@@ -347,7 +321,6 @@ inline uint64_t get_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k
khint_t t; khint_t t;
///query hash table，key is k
t = ha_pt_get(PCB->sub_h[sub_ID], sub_key); t = ha_pt_get(PCB->sub_h[sub_ID], sub_key);
if (t != kh_end(PCB->sub_h[sub_ID])) if (t != kh_end(PCB->sub_h[sub_ID]))
@@ -359,8 +332,6 @@ inline uint64_t get_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k
{ {
return (uint64_t)-1; return (uint64_t)-1;
} }
} }
inline uint64_t count_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k) inline uint64_t count_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k)
@@ -405,7 +376,6 @@ inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, in
if (occ) if (occ)
{ {
while (__sync_lock_test_and_set(&PCB->sub_h_lock[sub_ID].lock, 1)) while (__sync_lock_test_and_set(&PCB->sub_h_lock[sub_ID].lock, 1))
{ {
while (PCB->sub_h_lock[sub_ID].lock); while (PCB->sub_h_lock[sub_ID].lock);
@@ -413,16 +383,16 @@ inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, in
if (list[0].offset + 1 < occ) if (list[0].offset + 1 < occ)
{ {
list[0].offset++; list[0].offset++; // if not the last k-mer, this field is reused to keep the number of inserted positions
list[list[0].offset].readID = readID; list[list[0].offset].readID = readID;
///list[list[0].offset].readID = readID|direction;
list[list[0].offset].offset = pos; list[list[0].offset].offset = pos;
list[list[0].offset].rev = ha_code2rev(code);
} }
else else // now comes to the last k-mer position; then save it to list[0]
{ {
list[0].readID = readID; list[0].readID = readID;
///list[0].readID = readID|direction;
list[0].offset = pos; list[0].offset = pos;
list[0].rev = ha_code2rev(code);
flag = 1; flag = 1;
} }
@@ -434,18 +404,14 @@ inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, in
qsort(list, occ, sizeof(k_mer_pos), cmp_k_mer_pos); qsort(list, occ, sizeof(k_mer_pos), cmp_k_mer_pos);
} }
return 1; return 1;
} }
else else
{ {
return 0; return 0;
} }
} }
void init_Total_Count_Table(int k, Total_Count_Table* TCB); void init_Total_Count_Table(int k, Total_Count_Table* TCB);
void init_Total_Pos_Table(Total_Pos_Table* TCB, Total_Count_Table* pre_TCB); void init_Total_Pos_Table(Total_Pos_Table* TCB, Total_Count_Table* pre_TCB);
void destory_Total_Count_Table(Total_Count_Table* TCB); void destory_Total_Count_Table(Total_Count_Table* TCB);
@@ -473,12 +439,6 @@ uint64_t n_end_pos, uint8_t n_direction);
void merge_k_mer_pos_list_alloc_heap_sort(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT);
void Init_Heap(HeapSq* HBT);
void destory_Heap(HeapSq* HBT);
void clear_Heap(HeapSq* HBT);
void init_overlap_region_alloc(overlap_region_alloc* list); void init_overlap_region_alloc(overlap_region_alloc* list);
void clear_overlap_region_alloc(overlap_region_alloc* list); void clear_overlap_region_alloc(overlap_region_alloc* list);
void destory_overlap_region_alloc(overlap_region_alloc* list); void destory_overlap_region_alloc(overlap_region_alloc* list);
Makefile
+1 -1
View File
@@ -50,7 +50,7 @@ POA.o: POA.h Hash_Table.h khashl.h kmer.h Process_Read.h kseq.h Overlaps.h
POA.o: kvec.h kdq.h CommandLines.h Correct.h Levenshtein_distance.h POA.o: kvec.h kdq.h CommandLines.h Correct.h Levenshtein_distance.h
Process_Read.o: Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h Process_Read.o: Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h
Trio.o: khashl.h kthread.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h Trio.o: khashl.h kthread.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h
Trio.o: CommandLines.h Trio.h Trio.o: CommandLines.h Trio.h kmer.h
kmer.o: kmer.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h kmer.o: kmer.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h
kthread.o: kthread.h kthread.o: kthread.h
main.o: CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h main.o: CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h
Trio.cpp
+2 -25
View File
@@ -8,6 +8,7 @@
#include "Process_Read.h" #include "Process_Read.h"
#include "Trio.h" #include "Trio.h"
#include "CommandLines.h" #include "CommandLines.h"
#include "kmer.h"
#define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr)))) #define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr))))
#define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr)))) #define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr))))
@@ -37,30 +38,6 @@ unsigned char seq_nt4_table[256] = { // translate ACGT to 0123
4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4
}; };
static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function
{
key = (~key + (key << 21)) & mask; // key = (key << 21) - key - 1;
key = key ^ key >> 24;
key = ((key + (key << 3)) + (key << 8)) & mask; // key * 265
key = key ^ key >> 14;
key = ((key + (key << 2)) + (key << 4)) & mask; // key * 21
key = key ^ key >> 28;
key = (key + (key << 31)) & mask;
return key;
}
static inline uint64_t yak_hash64_64(uint64_t key)
{
key = ~key + (key << 21);
key = key ^ key >> 24;
key = (key + (key << 3)) + (key << 8);
key = key ^ key >> 14;
key = (key + (key << 2)) + (key << 4);
key = key ^ key >> 28;
key = key + (key << 31);
return key;
}
static inline uint64_t yak_hash_long(uint64_t x[4]) static inline uint64_t yak_hash_long(uint64_t x[4])
{ {
int j = x[1] < x[3]? 0 : 1; int j = x[1] < x[3]? 0 : 1;
@@ -361,4 +338,4 @@ void trio_partition()
fprintf(stderr, "Trio binning has been done.\n"); fprintf(stderr, "Trio binning has been done.\n");
fprintf(stderr, "%-30s%18.2f\n\n", "Trio binning time:", Get_T() - start_time); fprintf(stderr, "%-30s%18.2f\n\n", "Trio binning time:", Get_T() - start_time);
} }
kmer.cpp
-8
View File
@@ -10,14 +10,6 @@ void init_HPC_seq(HPC_seq* seq, char* str, long long l)
seq->str = str; seq->str = str;
} }
void init_Hash_code(Hash_code* code)
{
code->x[0] = 0;
code->x[1] = 0;
}
void init_small_hash_table(small_hash_table* x) void init_small_hash_table(small_hash_table* x)
{ {
x->size = 0; x->size = 0;
kmer.h
+38 -11
View File
@@ -16,7 +16,7 @@
typedef struct typedef struct
{ {
//can represent at most 64-mer //can represent at most 64-mer
uint64_t x[2]; uint64_t x[4];
} Hash_code; } Hash_code;
typedef struct { typedef struct {
@@ -83,16 +83,45 @@ inline uint64_t get_HPC_code(HPC_seq* seq, uint64_t* end_pos)
} }
inline void k_mer_append(Hash_code* code, uint64_t c, int k) inline void init_Hash_code(Hash_code* code)
{ {
code->x[0] = code->x[1] = code->x[2] = code->x[3] = 0;
uint64_t mask = ALL >> (64 -k);
code->x[0] = ((code->x[0]<<1) | (c&1)) & mask;
code->x[1] = ((code->x[1]<<1) | (c>>1)) & mask;
} }
inline void Hashcode_to_string(Hash_code* code, char* str, int k) inline void k_mer_append(Hash_code* code, uint64_t c, int k)
{
uint64_t mask = ALL >> (64 - k), shift = k - 1;
code->x[0] = ((code->x[0]<<1) | (c&1)) & mask;
code->x[1] = ((code->x[1]<<1) | (c>>1)) & mask;
code->x[2] = code->x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift;
code->x[3] = code->x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift;
}
static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function
{
key = (~key + (key << 21)) & mask; // key = (key << 21) - key - 1;
key = key ^ key >> 24;
key = ((key + (key << 3)) + (key << 8)) & mask; // key * 265
key = key ^ key >> 14;
key = ((key + (key << 2)) + (key << 4)) & mask; // key * 21
key = key ^ key >> 28;
key = (key + (key << 31)) & mask;
return key;
}
static inline uint64_t yak_hash64_64(uint64_t key)
{
key = ~key + (key << 21);
key = key ^ key >> 24;
key = (key + (key << 3)) + (key << 8);
key = key ^ key >> 14;
key = (key + (key << 2)) + (key << 4);
key = key ^ key >> 28;
key = key + (key << 31);
return key;
}
inline void Hashcode_to_string(Hash_code* code, char* str, int k) // FIXME: not working
{ {
uint8_t c; uint8_t c;
int i; int i;
@@ -108,7 +137,5 @@ inline void Hashcode_to_string(Hash_code* code, char* str, int k)
} }
void init_HPC_seq(HPC_seq* seq, char* str, long long l); void init_HPC_seq(HPC_seq* seq, char* str, long long l);
void init_Hash_code(Hash_code* code);
#endif
#endif