Reduce memory requirement

This commit is contained in:
Haoyu Cheng
2019-05-16 00:22:20 -04:00
parent 45355a1103
commit 3c5799aa1b
12 changed files with 320 additions and 39 deletions
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
+99 -13
View File
@@ -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
+3
View File
@@ -1,6 +1,9 @@
#ifndef __ASSEMBLY__
#define __ASSEMBLY__
#define FORWARD 0
#define REVERSE_COMPLEMENT (0x8000000000000000)
void Counting_multiple_thr();
+3
View File
@@ -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;
}
+157 -1
View File
@@ -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()
{
+44 -3
View File
@@ -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);
+1
View File
@@ -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;
}
+13 -22
View File
@@ -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
{