From 3c0c201cdf55ecdb7eab39e2fd5259e34e71c081 Mon Sep 17 00:00:00 2001 From: Haoyu Cheng Date: Mon, 26 Aug 2019 02:22:20 -0400 Subject: [PATCH] improve speed --- .vscode/settings.json | 5 +- Assembly.cpp | 86 +++++++++++++++++++++++++++++++++-- CommandLines.cpp | 12 +++-- Correct.cpp | 14 ++++++ Correct.h | 2 +- Hash_Table.cpp | 10 +++- Hash_Table.h | 2 + POA.cpp | 4 +- POA.h | 21 +++++++-- Process_Read.cpp | 103 +++++++++++++++++++++++++++++++++++------- Process_Read.h | 21 +++++++-- main.cpp | 8 +++- 12 files changed, 249 insertions(+), 39 deletions(-) diff --git a/.vscode/settings.json b/.vscode/settings.json index edb48a6..44c3105 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -7,6 +7,9 @@ "deque": "cpp", "vector": "cpp", "initializer_list": "cpp", - "*.tcc": "cpp" + "*.tcc": "cpp", + "bitset": "cpp", + "algorithm": "cpp", + "hashtable": "cpp" } } \ No newline at end of file diff --git a/Assembly.cpp b/Assembly.cpp index 0a2f486..2429e52 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -217,14 +217,17 @@ void* Build_hash_table(void* arg) ///HPC_base++; } + + ///load read compress_base(Get_READ(R_INF, curr_sub_block.read[i].ID), curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l, &R_INF.N_site[curr_sub_block.read[i].ID], HPC_read.N_occ); - + memcpy(R_INF.name+R_INF.name_index[curr_sub_block.read[i].ID], curr_sub_block.read[i].name.s, curr_sub_block.read[i].name.l); + /** ///reverse complement strand @@ -342,6 +345,7 @@ void Counting_multiple_thr() fprintf(stdout, "%-30s%18.2f\n\n", "Counting time:", Get_T() - start_time); + free(is_insert); } @@ -421,8 +425,10 @@ void Build_hash_table_multiple_thr() ///destory_All_reads(&R_INF); ///load_All_reads(&R_INF, read_file_name); } - + destory_Total_Count_Table(&TCB); + + free(is_insert); } @@ -717,6 +723,69 @@ Output_buffer_sub_block* current_sub_buffer) } +int get_required_read(const char *required_name, long long RID, All_reads* R_INF) +{ + + + + int required_name_length = strlen(required_name); + char* debug_name = Get_NAME((*R_INF), RID); + int debug_name_length = Get_NAME_LENGTH((*R_INF), RID); + int i; + + + if (required_name_length == debug_name_length) + { + for (i = 0; i < debug_name_length; i++) + { + if (required_name[i] != debug_name[i]) + { + break; + } + } + + if (i == debug_name_length) + { + fprintf(stderr, "required_name: %s\n", required_name); + + return 1; + ///fprintf(stderr, "read_length: %d\n", R_INF->g_read->length); + + /** + int aviable_overlap_name = 0; + for (i = 0; i < overlap_list->length; i++) + { + long long Len_x = overlap_list->list[i].x_pos_e - overlap_list->list[i].x_pos_s + 1; + + if (Len_x * OVERLAP_THRESHOLD <= overlap_list->list[i].align_length) + { + fprintf(stderr, "a_i: %d\n", aviable_overlap_name); + fprintf(stderr, "x_pos_s: %d, x_pos_e: %d\n", + overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e); + aviable_overlap_name++; + + debug_name = Get_NAME((*R_INF), overlap_list->list[i].y_id); + debug_name_length = Get_NAME_LENGTH((*R_INF), overlap_list->list[i].y_id); + + int j = 0; + for (j = 0; j < debug_name_length; j++) + { + fprintf(stderr, "%c", debug_name[j]); + } + fprintf(stderr, "\n"); + + } + } + **/ + } + + } + + + return 0; + +} + void* Overlap_calculate_heap_merge(void* arg) { @@ -781,10 +850,16 @@ void* Overlap_calculate_heap_merge(void* arg) init_buffer_sub_block(¤t_sub_buffer); + + for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num) { - - + /** + if(!get_required_read("m54238_180916_191625\/43253934\/ccs", i, &R_INF)) + { + continue; + } + **/ clear_Heap(&heap); clear_Candidates_list(&l); @@ -1069,6 +1144,9 @@ void* Overlap_calculate_heap_merge(void* arg) void Overlap_calculate_multipe_thr() { + + fprintf(stderr, "k_mer_min_freq: %lld, CORRECT_THRESHOLD: %f\n", k_mer_min_freq, CORRECT_THRESHOLD); + double start_time = Get_T(); init_output_buffer(thread_num); diff --git a/CommandLines.cpp b/CommandLines.cpp index a90da2d..0a50fee 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -10,8 +10,8 @@ char* output_file_name = NULL; int thread_num = 1; int k_mer_length = 40; //int k_mer_min_freq = 9; -//int k_mer_min_freq = 3; -int k_mer_min_freq = 2; +int k_mer_min_freq = 3; +//int k_mer_min_freq = 2; int k_mer_max_freq = 66; int load_index_from_disk = 0; int write_index_to_disk = 0; @@ -36,21 +36,25 @@ int CommandLine_process (int argc, char *argv[]) { static ko_longopt_t longopts[] = { - { "help", ko_no_argument, 100}, + { "help", ko_no_argument, 100}, { "seq", ko_required_argument, 101}, { "output", ko_required_argument, 102}, { "thread", ko_required_argument, 103}, + { "k_mer_min_freq", ko_required_argument, 104}, + { "k_mer_max_freq", ko_required_argument, 105}, { NULL, 0, 0 } }; ketopt_t opt = KETOPT_INIT; int i, c; - while ((c = ketopt(&opt, argc, argv, 1, "ht:o:q:k:lw", longopts)) >= 0) { + while ((c = ketopt(&opt, argc, argv, 1, "ht:o:q:k:lwm:n:", longopts)) >= 0) { if (c == 100 || c == 'h') Print_H(); else if (c == 103 || c == 't') thread_num = atoi(opt.arg); else if (c == 102 || c == 'o') output_file_name = opt.arg; else if (c == 101 || c == 'q') read_file_name = opt.arg; + else if (c == 104 || c == 'n') k_mer_min_freq = atoi(opt.arg); + else if (c == 105 || c == 'm') k_mer_max_freq = atoi(opt.arg); else if (c == 'k') k_mer_length = atoi(opt.arg); else if (c == 'l') load_index_from_disk = 1; else if (c == 'w') write_index_to_disk = 1; diff --git a/Correct.cpp b/Correct.cpp index daeb78b..09705ee 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -2604,6 +2604,20 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, generate_consensus(overlap_list, R_INF, g_read, dumy, g); + + + + + + + + + + + + + + /** ///fprintf(stderr, "length: %lld, corrected_base: %lld\n", g_read->length, dumy->corrected_base); diff --git a/Correct.h b/Correct.h index f727669..e147e32 100644 --- a/Correct.h +++ b/Correct.h @@ -5,7 +5,7 @@ #include "Levenshtein_distance.h" #include "POA.h" -#define CORRECT_THRESHOLD 0.65 +#define CORRECT_THRESHOLD 0.7 #define MIN_COVERAGE_THRESHOLD 4 #define CORRECT_INDEL_LENGTH 2 #define MISMATCH 1 diff --git a/Hash_Table.cpp b/Hash_Table.cpp index eacdf2a..61eef07 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -4,6 +4,7 @@ #include "Hash_Table.h" #include "Process_Read.h" #include "Correct.h" +#include "CommandLines.h" #include pthread_mutex_t output_mutex; @@ -1304,6 +1305,8 @@ void write_Total_Pos_Table(Total_Pos_Table* TCB, char* read_file_name) char* index_name = (char*)malloc(strlen(read_file_name)+5); sprintf(index_name, "%s.idx", read_file_name); FILE* fp = fopen(index_name, "w"); + fwrite(&k_mer_min_freq, sizeof(k_mer_min_freq), 1, fp); + fwrite(&k_mer_max_freq, sizeof(k_mer_max_freq), 1, fp); fwrite(&TCB->prefix_bits, sizeof(TCB->prefix_bits), 1, fp); fwrite(&TCB->suffix_bits, sizeof(TCB->suffix_bits), 1, fp); fwrite(&TCB->suffix_mode, sizeof(TCB->suffix_mode), 1, fp); @@ -1337,7 +1340,9 @@ int load_Total_Pos_Table(Total_Pos_Table* TCB, char* read_file_name) return 0; } - + + fread(&k_mer_min_freq, sizeof(k_mer_min_freq), 1, fp); + fread(&k_mer_max_freq, sizeof(k_mer_max_freq), 1, fp); fread(&TCB->prefix_bits, sizeof(TCB->prefix_bits), 1, fp); fread(&TCB->suffix_bits, sizeof(TCB->suffix_bits), 1, fp); fread(&TCB->suffix_mode, sizeof(TCB->suffix_mode), 1, fp); @@ -1403,6 +1408,9 @@ void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k khint_t t; ///这就是个迭代器 int absent; + /******************************************** + hash_table(key) ----> PCB->k_mer_index ------> PCB->pos + ********************************************/ for (i = 0; i < TCB->size; i++) { h = TCB->sub_h[i]; diff --git a/Hash_Table.h b/Hash_Table.h index 52b2fdf..e1fbb61 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -44,6 +44,7 @@ typedef struct Hash_table_spin_lock* sub_h_lock; int prefix_bits; int suffix_bits; + ///number of subtable int size; uint64_t suffix_mode; uint64_t non_unique_k_mer; @@ -163,6 +164,7 @@ typedef struct Hash_table_spin_lock* sub_h_lock; int prefix_bits; int suffix_bits; + ///number of subtable int size; uint64_t suffix_mode; k_mer_pos* pos; diff --git a/POA.cpp b/POA.cpp index 3f1819f..b7ed636 100644 --- a/POA.cpp +++ b/POA.cpp @@ -283,7 +283,7 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri else if (operation == 2) { ///cigar的起始和结尾不可能是2，所以这里-1没问题 - if (operationLen <= CORRECT_INDEL_LENGTH) + ///if (operationLen <= CORRECT_INDEL_LENGTH) { add_insertionEdge_weight(backbone, currentNodeID, y_string + y_i, operationLen); backbone->g_nodes.list[currentNodeID].num_insertions++; @@ -295,7 +295,7 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri ///3是y缺字符（x多字符），也就是backbone多字符 ///这个相当于在backbone对应字符处变成了‘——’ ///因此可以用mismatch类似的方法处理 - if (operationLen <= CORRECT_INDEL_LENGTH) + ///if (operationLen <= CORRECT_INDEL_LENGTH) { ///add_deletion_to_backbone(backbone, ¤tNodeID, operationLen); ///在编辑距离中，前面是个插入，后面是个删除，这种情况是不存在的 diff --git a/POA.h b/POA.h index 0e2e39e..822fbcd 100644 --- a/POA.h +++ b/POA.h @@ -203,9 +203,10 @@ inline void add_deletionEdge_weight(Graph* g, long long alignNodeID, long long d add_single_deletionEdge_weight(g, alignNodeID, alignNodeID + 2, 0); add_single_deletionEdge_weight(g, alignNodeID + 1, alignNodeID + 2, 0); } - else + else if (deletion_length > 2) { - fprintf(stderr, "too long deletion!\n"); + ///fprintf(stderr, "too long deletion!\n"); + add_single_deletionEdge_weight(g, alignNodeID, alignNodeID + deletion_length, 0); } } @@ -385,9 +386,21 @@ inline void add_insertionEdge_weight(Graph* g, long long alignNodeID, char* inse /**********************两个字符******************* */ } - else + else if (insert_length > 2) { - fprintf(stderr, "too long insertion\n"); + ////fprintf(stderr, "too long insertion\n"); + /*************************大于2个字符************************** */ + + edgeID = get_insertion_Edges(g, edge, insert_length, insert); + if (edgeID != -1) + { + ///这条路均只有一个出度 + edge->list[edgeID].weight++; + } + else + { + create_insertion_Edges(g, alignNodeID, insert_length, insert); + } } diff --git a/Process_Read.cpp b/Process_Read.cpp index c95d19c..c08f9e6 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -24,9 +24,15 @@ pthread_mutex_t i_doneMutex; void init_All_reads(All_reads* r) { r->index_size = READ_INIT_NUMBER; - r->index = (uint64_t*)malloc(sizeof(uint64_t)*r->index_size); - r->index[0] = 0; - r->read = NULL; + /**********should remove**********/ + ///r->index = (uint64_t*)malloc(sizeof(uint64_t)*r->index_size); + ///r->index[0] = 0; + ///r->read = NULL; + /**********should remove**********/ + r->read_length = (uint64_t*)malloc(sizeof(uint64_t)*r->index_size); + r->read_sperate = NULL; + + r->N_site = NULL; r->total_reads_bases = 0; @@ -50,11 +56,18 @@ void destory_All_reads(All_reads* r) { free(r->N_site[i]); } + free(r->read_sperate[i]); } free(r->N_site); - free(r->read); + free(r->read_sperate); + + + + ///free(r->read); free(r->name); free(r->name_index); + free(r->read_length); + } @@ -94,9 +107,22 @@ void write_All_reads(All_reads* r, char* read_file_name) } - fwrite(r->read, sizeof(uint8_t), (r->total_reads_bases/4 + r->total_reads + 5), fp); + /**********should remove**********/ + ///fwrite(r->index, sizeof(uint64_t), r->index_size, fp); + /**********should remove**********/ + fwrite(r->read_length, sizeof(uint64_t), r->total_reads, fp); + + /**********should remove**********/ + ///fwrite(r->read, sizeof(uint8_t), (r->total_reads_bases/4 + r->total_reads + 5), fp); + /**********should remove**********/ + for (i = 0; i < r->total_reads; i++) + { + fwrite(r->read_sperate[i], sizeof(uint8_t), r->read_length[i]/4+1, fp); + } + + + fwrite(r->name, sizeof(char), r->total_name_length, fp); - fwrite(r->index, sizeof(uint64_t), r->index_size, fp); fwrite(r->name_index, sizeof(uint64_t), r->name_index_size, fp); @@ -152,15 +178,29 @@ int load_All_reads(All_reads* r, char* read_file_name) } - r->read = (uint8_t*)malloc(sizeof(uint8_t)*(r->total_reads_bases/4 + r->total_reads + 5)); - fread(r->read, sizeof(uint8_t), (r->total_reads_bases/4 + r->total_reads + 5), fp); + + /**********should remove**********/ + ///r->index = (uint64_t*)malloc(sizeof(uint64_t)*r->index_size); + ///fread(r->index, sizeof(uint64_t), r->index_size, fp); + /**********should remove**********/ + r->read_length = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads); + fread(r->read_length, sizeof(uint64_t), r->total_reads, fp); + + /**********should remove**********/ + ///r->read = (uint8_t*)malloc(sizeof(uint8_t)*(r->total_reads_bases/4 + r->total_reads + 5)); + ///fread(r->read, sizeof(uint8_t), (r->total_reads_bases/4 + r->total_reads + 5), fp); + /**********should remove**********/ + r->read_sperate = (uint8_t**)malloc(sizeof(uint8_t*)*r->total_reads); + for (i = 0; i < r->total_reads; i++) + { + r->read_sperate[i] = (uint8_t*)malloc(sizeof(uint8_t)*(r->read_length[i]/4+1)); + fread(r->read_sperate[i], sizeof(uint8_t), r->read_length[i]/4+1, fp); + } + r->name = (char*)malloc(sizeof(char)*r->total_name_length); fread(r->name, sizeof(char), r->total_name_length, fp); - r->index = (uint64_t*)malloc(sizeof(uint64_t)*r->index_size); - fread(r->index, sizeof(uint64_t), r->index_size, fp); - r->name_index = (uint64_t*)malloc(sizeof(uint64_t)*r->name_index_size); fread(r->name_index, sizeof(uint64_t), r->name_index_size, fp); @@ -184,12 +224,20 @@ inline void insert_read(All_reads* r, kstring_t* read, kstring_t* name) if (r->index_size < r->total_reads + 2) { r->index_size = r->index_size * 2 + 2; - r->index = (uint64_t*)realloc(r->index,sizeof(uint64_t)*(r->index_size)); + /**********should remove**********/ + ///r->index = (uint64_t*)realloc(r->index,sizeof(uint64_t)*(r->index_size)); + /**********should remove**********/ + r->read_length = (uint64_t*)realloc(r->read_length,sizeof(uint64_t)*(r->index_size)); r->name_index_size = r->name_index_size * 2 + 2; r->name_index = (uint64_t*)realloc(r->name_index,sizeof(uint64_t)*(r->name_index_size)); } - r->index[r->total_reads] = r->index[r->total_reads-1] + read->l; + /**********should remove**********/ + ///r->index[r->total_reads] = r->index[r->total_reads-1] + read->l; + /**********should remove**********/ + r->read_length[r->total_reads - 1] = read->l; + + //r->index[r->total_reads] = r->index[r->total_reads-1] + read->l/4 + 1; r->name_index[r->total_reads] = r->name_index[r->total_reads-1] + name->l; @@ -197,11 +245,24 @@ inline void insert_read(All_reads* r, kstring_t* read, kstring_t* name) void malloc_All_reads(All_reads* r) { + ///必须加r->total_reads - r->read = (uint8_t*)malloc(sizeof(uint8_t)*(r->total_reads_bases/4 + r->total_reads + 5)); + /**********should remove**********/ + ///r->read = (uint8_t*)malloc(sizeof(uint8_t)*(r->total_reads_bases/4 + r->total_reads + 5)); + /**********should remove**********/ + r->read_sperate = (uint8_t**)malloc(sizeof(uint8_t*)*r->total_reads); + long long i = 0; + for (i = 0; i < r->total_reads; i++) + { + r->read_sperate[i] = (uint8_t*)malloc(sizeof(uint8_t)*(r->read_length[i]/4+1)); + } + + + + r->name = (char*)malloc(sizeof(char)*r->total_name_length); r->N_site = (uint64_t**)calloc(r->total_reads, sizeof(uint64_t*)); - + } void destory_UC_Read(UC_Read* r) @@ -369,6 +430,8 @@ void recover_UC_Read(UC_Read* r, All_reads* R_INF, uint64_t ID) r->seq[R_INF->N_site[ID][i]] = 'N'; } } + + r->RID = ID; } @@ -426,7 +489,8 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID) void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ) { - + ///N_site_lis saves the pos of all Ns in this read + ///N_site_lis[0] is the number of Ns if (N_site_occ) { (*N_site_lis) = (uint64_t*)malloc(sizeof(uint64_t)*(N_site_occ + 1)); @@ -442,7 +506,10 @@ void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_l uint64_t dest_i = 0; uint8_t tmp = 0; uint8_t c = 0; - + /** + fprintf(stderr, "src_l: %lld\n", src_l); + fflush(stderr); + **/ while (i + 4 <= src_l) @@ -597,6 +664,8 @@ inline void load_read_block(R_buffer_block* read_batch, int batch_read_size, read_batch->read[inner_i].ID = total_reads; total_reads++; + ///fprintf(stderr, "is_insert: %d\n", is_insert); + if (is_insert) { insert_read(&R_INF, &read_batch->read[inner_i].seq, &read_batch->read[inner_i].name); diff --git a/Process_Read.h b/Process_Read.h index ae8357a..4981e04 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -14,9 +14,11 @@ #define IS_FULL(buffer) ((buffer.num >= buffer.size)?1:0) #define IS_EMPTY(buffer) ((buffer.num == 0)?1:0) -#define Get_READ_LENGTH(R_INF, ID) (R_INF.index[ID+1] - R_INF.index[ID]) +///#define Get_READ_LENGTH(R_INF, ID) (R_INF.index[ID+1] - R_INF.index[ID]) +#define Get_READ_LENGTH(R_INF, ID) R_INF.read_length[ID] #define Get_NAME_LENGTH(R_INF, ID) (R_INF.name_index[ID+1] - R_INF.name_index[ID]) -#define Get_READ(R_INF, ID) R_INF.read + (R_INF.index[ID]>>2) + ID +///#define Get_READ(R_INF, ID) R_INF.read + (R_INF.index[ID]>>2) + ID +#define Get_READ(R_INF, ID) R_INF.read_sperate[ID] #define Get_NAME(R_INF, ID) R_INF.name + R_INF.name_index[ID] @@ -61,10 +63,20 @@ int get_read(kseq_t *s); typedef struct { uint64_t** N_site; - uint8_t* read; + ///uint8_t* read; char* name; - uint64_t* index; + + + uint8_t** read_sperate; + uint64_t* read_length; + + ///seq start pos in uint8_t* read + ///do not need it + ///uint64_t* index; uint64_t index_size; + + + ///name start pos in char* name uint64_t* name_index; uint64_t name_index_size; uint64_t total_reads; @@ -100,6 +112,7 @@ typedef struct char* seq; long long length; long long size; + long long RID; } UC_Read; void init_R_buffer(int thread_num); diff --git a/main.cpp b/main.cpp index 8d2bb8b..11ea66f 100644 --- a/main.cpp +++ b/main.cpp @@ -207,7 +207,8 @@ int main(int argc, char *argv[]) return 1; **/ - + fprintf(stdout, "defined k_mer_min_freq by user: %d\n", k_mer_min_freq); + fprintf(stdout, "defined k_mer_max_freq by user: %d\n", k_mer_max_freq); fprintf(stdout, "k-mer length: %d\n",k_mer_length); @@ -222,6 +223,11 @@ int main(int argc, char *argv[]) Build_hash_table_multiple_thr(); } + + fprintf(stdout, "k_mer_min_freq in hashtable: %d\n", k_mer_min_freq); + fprintf(stdout, "k_mer_max_freq in hashtable: %d\n", k_mer_max_freq); + + ///verify_Position_hash_table();