diff --git a/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch b/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch index f38974f..1b643f6 100644 Binary files a/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch and b/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch differ diff --git a/.vscode/ipch/856c69690c514632/KMER.ipch b/.vscode/ipch/856c69690c514632/KMER.ipch index 9114f47..f25e0b1 100644 Binary files a/.vscode/ipch/856c69690c514632/KMER.ipch and b/.vscode/ipch/856c69690c514632/KMER.ipch differ diff --git a/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch b/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch index 99f78ac..1b88f27 100644 Binary files a/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch and b/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch differ diff --git a/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch b/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch index f8248ca..a9c0d8c 100644 Binary files a/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch and b/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch differ diff --git a/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch b/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch index b61119d..238f6ea 100644 Binary files a/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch and b/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch differ diff --git a/Assembly.cpp b/Assembly.cpp index c4a693b..f0f441e 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -48,6 +48,7 @@ void* Perform_Counting(void* arg) for (i = 0; i < curr_sub_block.num; i++) { + ///forward strand init_HPC_seq(&HPC_read, curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l); init_Hash_code(&k_code); @@ -79,7 +80,44 @@ void* Perform_Counting(void* arg) } } - + + + /** + ///reverse complement strand + reverse_complement(curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l); + init_HPC_seq(&HPC_read, curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l); + init_Hash_code(&k_code); + + avalible_k = 0; + + while ((code = get_HPC_code(&HPC_read)) != 6) + { + if(code < 4) + { + k_mer_append(&k_code,code,k_mer_length); + avalible_k++; + if (avalible_k>=k_mer_length) + { + ///插入 + if(insert_Total_Count_Table(&TCB, &k_code, k_mer_length)) + { + select_k_mer_number++; + } + + k_mer_number++; + + } + + } + else + { + avalible_k = 0; + init_Hash_code(&k_code); + } + + } + **/ + } @@ -117,7 +155,6 @@ void* Build_hash_table(void* arg) int avalible_k = 0; - char* dest; while (file_flag != 0) { @@ -127,15 +164,57 @@ void* Build_hash_table(void* arg) for (i = 0; i < curr_sub_block.num; i++) { + + ///forward strand + init_HPC_seq(&HPC_read, curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l); + init_Hash_code(&k_code); - dest = R_INF.read+R_INF.index[curr_sub_block.read[i].ID]; - memcpy(dest, curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l); + avalible_k = 0; + + HPC_base = 0; + + while ((code = get_HPC_code(&HPC_read)) != 6) + { + if(code < 4) + { + k_mer_append(&k_code,code,k_mer_length); + avalible_k++; + if (avalible_k>=k_mer_length) + { + + ///选取的k-mer满足两个要求 + ///1. hash(k-mer) % 101 <= 3 + ///2. occ(k-mer)要满足范围 + ///TCB表中的元素仅满足第一个要求,而PCB表中的元素满足两个要求 + ///所以如果当前k-mer在PCB表中存在,则他的位置一定要加入到候选位置中去 + ///insert_Total_Pos_Table(&PCB, &k_code, k_mer_length, curr_sub_block.read[i].ID, HPC_base - k_mer_length + 1, FORWARD); + insert_Total_Pos_Table(&PCB, &k_code, k_mer_length, + curr_sub_block.read[i].ID, HPC_base - k_mer_length + 1); + } + + } + else + { + avalible_k = 0; + init_Hash_code(&k_code); + } + + HPC_base++; + } - dest = R_INF.name+R_INF.name_index[curr_sub_block.read[i].ID]; - memcpy(dest, curr_sub_block.read[i].name.s, curr_sub_block.read[i].name.l); - - + ///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 + reverse_complement(curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l); + init_HPC_seq(&HPC_read, curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l); init_Hash_code(&k_code); @@ -158,7 +237,7 @@ void* Build_hash_table(void* arg) ///TCB表中的元素仅满足第一个要求,而PCB表中的元素满足两个要求 ///所以如果当前k-mer在PCB表中存在,则他的位置一定要加入到候选位置中去 insert_Total_Pos_Table(&PCB, &k_code, k_mer_length, - curr_sub_block.read[i].ID, HPC_base - k_mer_length + 1); + curr_sub_block.read[i].ID, HPC_base - k_mer_length + 1, REVERSE_COMPLEMENT); } } @@ -170,8 +249,7 @@ void* Build_hash_table(void* arg) HPC_base++; } - - + **/ } @@ -269,6 +347,7 @@ void Build_hash_table_multiple_thr() Traverse_Counting_Table(&TCB, &PCB, k_mer_min_freq, k_mer_max_freq); fprintf(stdout, "%-30s%18.2f\n\n", "Traverse time:", Get_T() - T_start_time); + fflush(stdout); @@ -408,6 +487,9 @@ void verify_Position_hash_table() k_mer_pos* list; uint64_t sub_ID; + UC_Read g_read; + init_UC_Read(&g_read); + fprintf(stdout, "Start Verifying Position Table...\n"); while (get_read(seq)) @@ -419,9 +501,13 @@ void verify_Position_hash_table() fprintf(stderr, "seq error\n"); } - if(memcmp(seq->seq.s, Get_READ(R_INF, read_number), seq->seq.l)) + recover_UC_Read(&g_read, &R_INF, read_number); + + + if(memcmp(seq->seq.s, g_read.seq, seq->seq.l)) { - fprintf(stderr, "seq error\n"); + fprintf(stderr, "\nseq error ID: %llu, length: %llu\n",read_number, seq->seq.l); + } if (seq->name.l diff --git a/Assembly.h b/Assembly.h index 38fe873..eae139c 100644 --- a/Assembly.h +++ b/Assembly.h @@ -1,6 +1,9 @@ #ifndef __ASSEMBLY__ #define __ASSEMBLY__ +#define FORWARD 0 +#define REVERSE_COMPLEMENT (0x8000000000000000) + void Counting_multiple_thr(); diff --git a/Hash_Table.h b/Hash_Table.h index 3cac6a4..a3074f7 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -234,6 +234,7 @@ inline uint64_t locate_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, k_ int cmp_k_mer_pos(const void * a, const void * b); +//inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k, uint64_t readID, uint64_t pos, uint64_t direction) inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k, uint64_t readID, uint64_t pos) { k_mer_pos* list; @@ -253,11 +254,13 @@ inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, in { list[0].offset++; list[list[0].offset].readID = readID; + ///list[list[0].offset].readID = readID|direction; list[list[0].offset].offset = pos; } else { list[0].readID = readID; + ///list[0].readID = readID|direction; list[0].offset = pos; flag = 1; } diff --git a/Process_Read.cpp b/Process_Read.cpp index c773616..e8ae813 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -27,6 +27,7 @@ void init_All_reads(All_reads* r) r->index = (uint64_t*)malloc(sizeof(uint64_t)*r->index_size); r->index[0] = 0; r->read = NULL; + r->N_site = NULL; r->total_reads_bases = 0; @@ -57,14 +58,141 @@ inline void insert_read(All_reads* r, kstring_t* read, kstring_t* name) 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; + //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; } void malloc_All_reads(All_reads* r) { - r->read = (char*)malloc(sizeof(char)*r->total_reads_bases); + ///必须加r->total_reads + r->read = (uint8_t*)malloc(sizeof(uint8_t)*(r->total_reads_bases/4 + r->total_reads + 5)); r->name = (char*)malloc(sizeof(char)*r->total_name_length); + r->N_site = (uint64_t**)malloc(sizeof(uint64_t*)*r->total_reads); +} + +void init_UC_Read(UC_Read* r) +{ + r->length = 0; + r->size = 0; + r->seq = NULL; + if (bit_t_seq_table[0][0] == 0) + { + uint64_t i = 0; + + for (i = 0; i < 256; i++) + { + bit_t_seq_table[i][0] = s_H[((i >> 6)&(uint64_t)3)]; + bit_t_seq_table[i][1] = s_H[((i >> 4)&(uint64_t)3)]; + bit_t_seq_table[i][2] = s_H[((i >> 2)&(uint64_t)3)]; + bit_t_seq_table[i][3] = s_H[(i&(uint64_t)3)]; + } + + } + +} + + +void recover_UC_Read(UC_Read* r, All_reads* R_INF, uint64_t ID) +{ + r->length = Get_READ_LENGTH((*R_INF), ID); + uint8_t* src = Get_READ((*R_INF), ID); + + if (r->length + 4 > r->size) + { + r->size = r->length + 4; + r->seq = (char*)realloc(r->seq,sizeof(char)*(r->size)); + } + + uint64_t i = 0; + + while (i < r->length) + { + memcpy(r->seq+i, bit_t_seq_table[src[i>>2]], 4); + i = i + 4; + } + + + if (R_INF->N_site[ID]) + { + for (i = 1; i <= R_INF->N_site[ID][0]; i++) + { + r->seq[R_INF->N_site[ID][i]] = 'N'; + } + } + +} + + + +#define COMPRESS_BASE {c = seq_nt6_table[src[i]];\ + if (c >= 4)\ + {\ + c = 0;\ + (*N_site_lis)[N_site_i] = i;\ + N_site_i++;\ + }\ + i++;}\ + +void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ) +{ + + + if (N_site_occ) + { + (*N_site_lis) = (uint64_t*)malloc(sizeof(uint64_t)*(N_site_occ + 1)); + (*N_site_lis)[0] = N_site_occ; + } + else + { + (*N_site_lis) = NULL; + } + + uint64_t i = 0; + uint64_t N_site_i = 1; + uint64_t dest_i = 0; + uint8_t tmp = 0; + uint8_t c = 0; + + + + while (i + 4 <= src_l) + { + tmp = 0; + + COMPRESS_BASE; + tmp = tmp | (c<<6); + + COMPRESS_BASE; + tmp = tmp | (c<<4); + + COMPRESS_BASE; + tmp = tmp | (c<<2); + + COMPRESS_BASE; + tmp = tmp | c; + + dest[dest_i] = tmp; + dest_i++; + } + + //最多还剩3个字符 + uint64_t shift = 6; + if (i < src_l) + { + tmp = 0; + + while (i < src_l) + { + COMPRESS_BASE; + tmp = tmp | (c << shift); + shift = shift -2; + } + + dest[dest_i] = tmp; + dest_i++; + } + } void init_kseq(char* file) @@ -348,6 +476,34 @@ int get_reads_mul_thread(R_buffer_block* curr_sub_block) + + + + + +void reverse_complement(char* pattern, uint64_t length) +{ + int i = 0; + uint64_t end = length / 2; + char k; + uint64_t index; + + for (i = 0; i < end; i++) + { + index = length - i - 1; + k = pattern[index]; + pattern[index] = RC_CHAR(pattern[i]); + pattern[i] = RC_CHAR(k); + } + + if(length&(uint64_t)1) + { + pattern[end] = RC_CHAR(pattern[end]); + } + +} + + void Counting_block() { diff --git a/Process_Read.h b/Process_Read.h index a02d75b..9a537bb 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -16,13 +16,41 @@ #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_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]) -#define Get_NAME(R_INF, ID) (R_INF.name + R_INF.name_index[ID]) +#define Get_READ(R_INF, ID) R_INF.read + (R_INF.index[ID]>>2) + ID +#define Get_NAME(R_INF, ID) R_INF.name + R_INF.name_index[ID] KSEQ_INIT(gzFile, gzread) + + +static uint8_t seq_nt6_table[256] = { + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 0, 5, 1, 5, 5, 5, 2, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 0, 5, 1, 5, 5, 5, 2, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5 +}; + +static char bit_t_seq_table[256][4] = {0}; +static char s_H[4] = {'A', 'C', 'G', 'T'}; +static char rc_Table[4] = {'T', 'G', 'C', 'A'}; + +#define RC_CHAR(x) rc_Table[seq_nt6_table[(uint8_t)x]] + + void init_kseq(char* file); void destory_kseq(); int get_read(kseq_t *s); @@ -31,7 +59,8 @@ int get_read(kseq_t *s); typedef struct { - char* read; + uint64_t** N_site; + uint8_t* read; char* name; uint64_t* index; uint64_t index_size; @@ -64,11 +93,23 @@ typedef struct int all_read_end; } R_buffer; + +typedef struct +{ + char* seq; + long long length; + long long size; +} UC_Read; + void init_R_buffer(int thread_num); void init_All_reads(All_reads* r); void* input_reads_muti_threads(void*); void init_R_buffer_block(R_buffer_block* curr_sub_block); int get_reads_mul_thread(R_buffer_block* curr_sub_block); +void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ); +void init_UC_Read(UC_Read* r); +void recover_UC_Read(UC_Read* r, All_reads* R_INF, uint64_t ID); +void reverse_complement(char* pattern, uint64_t length); void Counting_block(); void destory_R_buffer_block(R_buffer_block* curr_sub_block); diff --git a/kmer.cpp b/kmer.cpp index ec3b380..e9e9647 100644 --- a/kmer.cpp +++ b/kmer.cpp @@ -6,6 +6,7 @@ void init_HPC_seq(HPC_seq* seq, char* str, long long l) { seq->i = 0; seq->l = l; + seq->N_occ = 0; seq->str = str; } diff --git a/kmer.h b/kmer.h index 55e35a3..9b1e5c7 100644 --- a/kmer.h +++ b/kmer.h @@ -2,29 +2,13 @@ #define __KMER__ #include "Process_Read.h" -#define ALL (0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff) +///#define ALL (0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff) +#define ALL (0xffffffffffffffff) #define SAFE_SHIFT(k) k & ((k < 64)?ALL:0) -static unsigned char seq_nt6_table[256] = { - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 0, 5, 1, 5, 5, 5, 2, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 0, 5, 1, 5, 5, 5, 2, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, - 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5 -}; + @@ -44,6 +28,7 @@ typedef struct char* str; long long l; long long i; + long long N_occ; } HPC_seq; @@ -53,17 +38,23 @@ inline uint64_t get_HPC_code(HPC_seq* seq) if(seq->i < seq ->l) { - char code = seq->str[seq->i]; + uint8_t code = seq_nt6_table[(uint8_t)seq->str[seq->i]]; for (; seq->i < seq->l; seq->i++) { - if (seq->str[seq->i] != code) + ///统计N的个数 + if (seq_nt6_table[(uint8_t)seq->str[seq->i]] >= 4) + { + seq->N_occ++; + } + + if (seq_nt6_table[(uint8_t)seq->str[seq->i]] != code) { break; } } - return (uint64_t)seq_nt6_table[(uint8_t)code]; + return (uint64_t)code; } else {