diff --git a/Process_Read.cpp b/Process_Read.cpp index 2b4d3f4..e461714 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -81,7 +81,6 @@ void destory_All_reads(All_reads* r) free(r->trio_flag); } - void write_All_reads(All_reads* r, char* read_file_name) { fprintf(stderr, "Writing reads to disk... \n"); @@ -112,9 +111,6 @@ void write_All_reads(All_reads* r, char* read_file_name) { fwrite(&zero, sizeof(zero), 1, fp); } - - - } fwrite(r->read_length, sizeof(uint64_t), r->total_reads, fp); @@ -132,8 +128,6 @@ void write_All_reads(All_reads* r, char* read_file_name) fprintf(stderr, "Reads has been written.\n"); } - - int load_All_reads(All_reads* r, char* read_file_name) { fprintf(stderr, "Loading reads from disk... \n"); @@ -164,7 +158,6 @@ int load_All_reads(All_reads* r, char* read_file_name) r->N_site = (uint64_t**)malloc(sizeof(uint64_t*)*r->total_reads); for (i = 0; i < r->total_reads; i++) { - f_flag += fread(&zero, sizeof(zero), 1, fp); if (zero) @@ -233,32 +226,33 @@ int load_All_reads(All_reads* r, char* read_file_name) return 1; } - - -inline void insert_read(All_reads* r, kstring_t* read, kstring_t* name) +void ha_insert_read_len(All_reads *r, int read_len, int name_len) { r->total_reads++; - r->total_reads_bases = r->total_reads_bases + read->l; - r->total_name_length = r->total_name_length + name->l; + r->total_reads_bases += (uint64_t)read_len; + r->total_name_length += (uint64_t)name_len; - ///must +1 - if (r->index_size < r->total_reads + 2) - { + // must +1 + if (r->index_size < r->total_reads + 2) { r->index_size = r->index_size * 2 + 2; - r->read_length = (uint64_t*)realloc(r->read_length,sizeof(uint64_t)*(r->index_size)); + 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->name_index = (uint64_t*)realloc(r->name_index, sizeof(uint64_t) * r->name_index_size); } - r->read_length[r->total_reads - 1] = read->l; - r->name_index[r->total_reads] = r->name_index[r->total_reads-1] + name->l; + r->read_length[r->total_reads - 1] = read_len; + r->name_index[r->total_reads] = r->name_index[r->total_reads - 1] + name_len; +} + +static inline void insert_read(All_reads* r, kstring_t* read, kstring_t* name) +{ + ha_insert_read_len(r, read->l, name->l); } void malloc_All_reads(All_reads* r) { - r->read_size = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads); - memcpy (r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads); + memcpy(r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads); r->read_sperate = (uint8_t**)malloc(sizeof(uint8_t*)*r->total_reads); long long i = 0; @@ -313,7 +307,6 @@ void init_aux_table() bit_t_seq_table_rc[i][2] = RC_CHAR(bit_t_seq_table[i][1]); bit_t_seq_table_rc[i][3] = RC_CHAR(bit_t_seq_table[i][0]); } - } } @@ -338,18 +331,12 @@ void init_UC_Read(UC_Read* r) bit_t_seq_table_rc[i][2] = RC_CHAR(bit_t_seq_table[i][1]); bit_t_seq_table_rc[i][3] = RC_CHAR(bit_t_seq_table[i][0]); } - } - } - -void recover_UC_Read_sub_region_begin_end -(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID, int extra_begin, int extra_end) +void recover_UC_Read_sub_region_begin_end(char* r, long long start_pos, long long length, uint8_t strand, + All_reads* R_INF, long long ID, int extra_begin, int extra_end) { - - - long long readLen = Get_READ_LENGTH((*R_INF), ID); uint8_t* src = Get_READ((*R_INF), ID); @@ -357,17 +344,11 @@ void recover_UC_Read_sub_region_begin_end long long copyLen; long long end_pos = start_pos + length - 1; - - - if (strand == 0) { - i = start_pos; copyLen = 0; - - long long initLen = start_pos % 4; if (initLen != 0) @@ -376,8 +357,6 @@ void recover_UC_Read_sub_region_begin_end copyLen = copyLen + 4 - initLen; i = i + copyLen; } - - while (copyLen < length) { memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4); @@ -385,7 +364,6 @@ void recover_UC_Read_sub_region_begin_end i = i + 4; } - if (R_INF->N_site[ID]) { for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) @@ -400,17 +378,12 @@ void recover_UC_Read_sub_region_begin_end } } } - - } else { - start_pos = readLen - start_pos - 1; end_pos = readLen - end_pos - 1; - - ///start_pos > end_pos i = start_pos; copyLen = 0; @@ -436,7 +409,6 @@ void recover_UC_Read_sub_region_begin_end for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) { - if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos) { r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N'; @@ -447,21 +419,11 @@ void recover_UC_Read_sub_region_begin_end } } } - - - } - } - - - void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID) { - - - long long readLen = Get_READ_LENGTH((*R_INF), ID); uint8_t* src = Get_READ((*R_INF), ID); @@ -471,7 +433,6 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, if (strand == 0) { - i = start_pos; copyLen = 0; @@ -483,8 +444,7 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, copyLen = copyLen + 4 - initLen; i = i + copyLen; } - - + while (copyLen < length) { memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4); @@ -492,7 +452,6 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, i = i + 4; } - if (R_INF->N_site[ID]) { for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) @@ -507,17 +466,12 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, } } } - - } else { - start_pos = readLen - start_pos - 1; end_pos = readLen - end_pos - 1; - - ///start_pos > end_pos i = start_pos; copyLen = 0; @@ -543,7 +497,6 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) { - if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos) { r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N'; @@ -554,11 +507,7 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, } } } - - - } - } @@ -623,7 +572,6 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID) index = index + 4; } - if (R_INF->N_site[ID]) { for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) @@ -631,11 +579,8 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID) r->seq[r->length - R_INF->N_site[ID][i] - 1] = 'N'; } } - } - - #define COMPRESS_BASE {c = seq_nt6_table[(uint8_t)src[i]];\ if (c >= 4)\ {\ @@ -647,7 +592,6 @@ 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) @@ -666,10 +610,8 @@ void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_l uint8_t tmp = 0; uint8_t c = 0; - while (i + 4 <= src_l) { - tmp = 0; COMPRESS_BASE; @@ -705,12 +647,8 @@ void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_l dest[dest_i] = tmp; dest_i++; } - } - - - void open_file(gz_files* nfps, char* name) { nfps->fp = gzopen(name, "r"); @@ -722,7 +660,6 @@ void open_file(gz_files* nfps, char* name) nfps->seq = kseq_init(nfps->fp); } - void close_file(gz_files* nfps) { kseq_destroy(nfps->seq); @@ -774,11 +711,9 @@ inline void exchage_kstring_t(kstring_t* a, kstring_t* b) int get_read(kseq_t *s, int adapterLen) { int l; - ///if ((l = kseq_read(seq)) >= 0) if ((l = read_item()) >= 0) { - exchage_kstring_t(&(fps.seq->comment), &s->comment); exchage_kstring_t(&(fps.seq->name), &s->name); exchage_kstring_t(&(fps.seq->qual), &s->qual); @@ -801,15 +736,12 @@ int get_read(kseq_t *s, int adapterLen) } } - return 1; } else { return 0; } - - } void init_R_buffer_block(R_buffer_block* curr_sub_block) @@ -823,6 +755,7 @@ void clear_R_buffer() RDB.all_read_end = 0; RDB.num = 0; } + void init_R_buffer(int thread_num) { RDB.all_read_end = 0; @@ -841,13 +774,11 @@ void init_R_buffer(int thread_num) } - void destory_R_buffer_block(R_buffer_block* curr_sub_block) { kseq_destroy(curr_sub_block->read); } - void destory_R_buffer() { int i = 0; @@ -858,7 +789,6 @@ void destory_R_buffer() } free(RDB.sub_block); - } @@ -868,12 +798,8 @@ inline void load_read_block(R_buffer_block* read_batch, int batch_read_size, int inner_i = 0; int file_flag = 1; - - - while (inner_iread[inner_i], adapterLen); if (file_flag == 1) @@ -881,7 +807,6 @@ inline void load_read_block(R_buffer_block* read_batch, int batch_read_size, read_batch->read[inner_i].ID = total_reads; total_reads++; - if (is_insert) { insert_read(&R_INF, &read_batch->read[inner_i].seq, @@ -903,14 +828,11 @@ inline void load_read_block(R_buffer_block* read_batch, int batch_read_size, *return_file_flag = file_flag; read_batch->num = inner_i; - } inline void push_R_block(R_buffer_block* tmp_sub_block) { - - ///only exchange pointers kseq_t *k1; k1 = RDB.sub_block[RDB.num].read; @@ -940,29 +862,20 @@ inline void pop_R_block(R_buffer_block* curr_sub_block) curr_sub_block->num = RDB.sub_block[RDB.num].num; RDB.sub_block[RDB.num].num = 0; - - - } - - void* input_reads_muti_threads(void* arg) { int is_insert = *((int*)arg); - total_reads = 0; - int file_flag = 1; R_buffer_block tmp_buf; init_R_buffer_block(&tmp_buf); - - while (1) { load_read_block(&tmp_buf, RDB.block_inner_size, &file_flag, is_insert, asm_opt.adapterLen); @@ -972,7 +885,6 @@ void* input_reads_muti_threads(void* arg) break; } - pthread_mutex_lock(&i_readinputMutex); while (IS_FULL(RDB)) { @@ -981,14 +893,12 @@ void* input_reads_muti_threads(void* arg) pthread_cond_wait(&i_readinputflushCond, &i_readinputMutex); } - push_R_block(&tmp_buf); pthread_cond_signal(&i_readinputstallCond); pthread_mutex_unlock(&i_readinputMutex); } - pthread_mutex_lock(&i_readinputMutex); RDB.all_read_end = 1; pthread_cond_signal(&i_readinputstallCond); //important @@ -999,34 +909,25 @@ void* input_reads_muti_threads(void* arg) fprintf(stderr, "Reads #: %lu\n", (unsigned long)total_reads); fprintf(stderr, "Bases #: %lu\n", (unsigned long)R_INF.total_reads_bases); - return NULL; } - - int get_reads_mul_thread(R_buffer_block* curr_sub_block) { - - pthread_mutex_lock(&i_readinputMutex); - while (IS_EMPTY(RDB) && RDB.all_read_end == 0) { - pthread_cond_signal(&i_readinputflushCond); pthread_cond_wait(&i_readinputstallCond, &i_readinputMutex); } - if (!IS_EMPTY(RDB)) { pop_R_block(curr_sub_block); pthread_cond_signal(&i_readinputflushCond); pthread_mutex_unlock(&i_readinputMutex); - return 1; } else @@ -1039,17 +940,8 @@ int get_reads_mul_thread(R_buffer_block* curr_sub_block) return 0; } - - } - - - - - - - void reverse_complement(char* pattern, uint64_t length) { uint64_t i = 0; @@ -1069,11 +961,8 @@ void reverse_complement(char* pattern, uint64_t length) { pattern[end] = RC_CHAR(pattern[end]); } - } - - typedef struct { char* tmp; long long tmpSize; diff --git a/Process_Read.h b/Process_Read.h index 23764e6..e359328 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -25,11 +25,9 @@ #define Get_NAME(R_INF, ID) R_INF.name + R_INF.name_index[(ID)] - KSEQ_INIT(gzFile, gzread) - extern uint8_t seq_nt6_table[256]; extern char bit_t_seq_table[256][4]; extern char bit_t_seq_table_rc[256][4]; @@ -43,7 +41,6 @@ extern char rc_Table[5]; void init_aux_table(); int get_read(kseq_t *s, int adapterLen); - typedef struct { uint64_t x_id; @@ -76,7 +73,6 @@ inline void init_PAF_alloc(PAF_alloc* list) list->list = (PAF*)malloc(sizeof(PAF)*list->size); } - inline void append_PAF_alloc(PAF_alloc* list, PAF* e) { if(list->length+1 > list->size) @@ -89,9 +85,6 @@ inline void append_PAF_alloc(PAF_alloc* list, PAF* e) list->length++; } - - - typedef struct { /**[0-1] bits are type:**/ @@ -104,7 +97,7 @@ typedef struct uint32_t lost_base_length; uint32_t lost_base_size; uint32_t new_length; -}Compressed_Cigar_record; +} Compressed_Cigar_record; #define AMBIGU 0 #define FATHER 1 @@ -112,13 +105,13 @@ typedef struct #define MIX_TRIO 3 #define NON_TRIO 4 #define DROP 5 + typedef struct { uint64_t** N_site; ///uint8_t* read; char* name; - uint8_t** read_sperate; uint64_t* read_length; uint64_t* read_size; @@ -129,7 +122,6 @@ typedef struct ///uint64_t* index; uint64_t index_size; - ///name start pos in char* name uint64_t* name_index; uint64_t name_index_size; @@ -143,7 +135,6 @@ typedef struct ma_hit_t_alloc* paf; ma_hit_t_alloc* reverse_paf; ma_sub_t* coverage_cut; - } All_reads; extern All_reads R_INF; @@ -154,10 +145,8 @@ typedef struct { kseq_t* read; long long num; - } R_buffer_block; - typedef struct { R_buffer_block* sub_block; @@ -167,7 +156,6 @@ typedef struct int all_read_end; } R_buffer; - typedef struct { char* seq; @@ -208,5 +196,6 @@ void clear_R_buffer(); void init_gz_files(hifiasm_opt_t* asm_opt); void destory_gz_files(); +void ha_insert_read_len(All_reads *r, int read_len, int name_len); #endif diff --git a/htab.cpp b/htab.cpp index 1a0f055..1de88d2 100644 --- a/htab.cpp +++ b/htab.cpp @@ -9,6 +9,7 @@ #include "kseq.h" #include "ksort.h" #include "htab.h" +#include "Process_Read.h" #include "CommandLines.h" #define YAK_COUNTER_BITS 12 @@ -308,7 +309,6 @@ static void worker_pt_gen(void *data, long i, int tid) // callback for kt_for() b->n += kh_key(g, k) & YAK_MAX_COUNT; } } -// fprintf(stderr, "X\t%ld\t%d\t%ld\n", i, kh_size(g), (long)b->n); yak_ct_destroy(g); a->ct->h[i].h = 0; CALLOC(b->a, b->n); @@ -464,12 +464,15 @@ static void count_seq_buf_HPC(ch_buf_t *buf, int k, int p, int len, const char * * K-mer counting * ******************/ -KSEQ_INIT(gzFile, gzread) +//KSEQ_INIT(gzFile, gzread) + +#define HAF_COUNT_EXACT 0x1 +#define HAF_COUNT_ALL 0x2 typedef struct { // global data structure for kt_pipeline() const yak_copt_t *opt; const void *flt_tab; - int create_new, is_store; + int flag, create_new, is_store; uint64_t batch_offset; uint64_t n_base; kseq_t *ks; @@ -608,7 +611,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip return 0; } -static ha_ct_t *yak_count(const char *fn, const yak_copt_t *opt, ha_pt_t *p0, ha_ct_t *c0, const void *flt_tab) +static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt_t *p0, ha_ct_t *c0, const void *flt_tab) { pl_data_t pl; gzFile fp; @@ -633,28 +636,28 @@ static ha_ct_t *yak_count(const char *fn, const yak_copt_t *opt, ha_pt_t *p0, ha return pl.ct; } -static ha_ct_t *yak_count_file(const yak_copt_t *opt, ha_pt_t *p0, int n_fn, char **fn, const void *flt_tab) +static ha_ct_t *yak_count_file(const yak_copt_t *opt, int flag, ha_pt_t *p0, int n_fn, char **fn, const void *flt_tab) { int i; ha_ct_t *h = 0; for (i = 0; i < n_fn; ++i) - h = yak_count(fn[i], opt, p0, h, flt_tab); + h = yak_count(opt, fn[i], flag, p0, h, flt_tab); if (h && opt->bf_shift > 0) ha_ct_destroy_bf(h); return h; } -ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, ha_pt_t *p0, int is_exact, int count_all, const void *flt_tab) +ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const void *flt_tab) { yak_copt_t opt; ha_ct_t *h; yak_copt_init(&opt); opt.k = asm_opt->k_mer_length; opt.is_HPC = !asm_opt->no_HPC; - opt.w = count_all? 1 : asm_opt->mz_win; - opt.bf_shift = is_exact? 0 : asm_opt->bf_shift; + opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win; + opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift; opt.n_thread = asm_opt->thread_num; - h = yak_count_file(&opt, p0, asm_opt->num_reads, asm_opt->read_file_names, flt_tab); + h = yak_count_file(&opt, flag, p0, asm_opt->num_reads, asm_opt->read_file_names, flt_tab); return h; } @@ -707,7 +710,7 @@ void *ha_gen_flt_tab(const hifiasm_opt_t *asm_opt) int64_t cnt[YAK_N_COUNTS]; int peak_hom, peak_het, cutoff; ha_ct_t *h; - h = ha_count(asm_opt, 0, 0, 1, 0); + h = ha_count(asm_opt, HAF_COUNT_ALL, NULL, NULL); ha_ct_hist(h, cnt, asm_opt->thread_num); peak_hom = yak_analyze_count(YAK_N_COUNTS, cnt, &peak_het); if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); @@ -727,7 +730,7 @@ void *ha_gen_mzidx(const hifiasm_opt_t *asm_opt, const void *flt_tab) int peak_hom, peak_het, i; ha_ct_t *ct; ha_pt_t *pt; - ct = ha_count(asm_opt, 0, 1, 0, flt_tab); + ct = ha_count(asm_opt, HAF_COUNT_EXACT, NULL, flt_tab); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__, yak_realtime(), yak_cputime() / yak_realtime(), (long)ct->tot); ha_ct_hist(ct, cnt, asm_opt->thread_num); @@ -737,7 +740,7 @@ void *ha_gen_mzidx(const hifiasm_opt_t *asm_opt, const void *flt_tab) ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num); for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i; pt = ha_pt_gen(ct, asm_opt->thread_num); - ha_count(asm_opt, pt, 1, 0, flt_tab); + ha_count(asm_opt, HAF_COUNT_EXACT, pt, flt_tab); assert((uint64_t)tot_cnt == pt->tot_pos); ha_pt_sort(pt, asm_opt->thread_num); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__,