diff --git a/Assembly.cpp b/Assembly.cpp index f0a3aa6..4afb606 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -265,7 +265,7 @@ void* Build_hash_table(void* arg) } ///load read - compress_base(Get_READ(R_INF, curr_sub_block.read[i].ID), + ha_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); @@ -1170,7 +1170,7 @@ void* Save_corrected_reads(void* arg) R_INF.read_length[i] = new_read_length; - compress_base(Get_READ(R_INF, i), + ha_compress_base(Get_READ(R_INF, i), new_read, new_read_length, &R_INF.N_site[i], N_occ); } diff --git a/Assembly.h b/Assembly.h index ef7caba..557aa8a 100644 --- a/Assembly.h +++ b/Assembly.h @@ -8,9 +8,6 @@ #define Get_Cigar_Type(RECORD) (RECORD&3) #define Get_Cigar_Length(RECORD) (RECORD>>2) -void *ha_gen_flt_tab(const hifiasm_opt_t *asm_opt); -void *ha_gen_mzidx(const hifiasm_opt_t *asm_opt, const void *flt_tab); - void Counting_multiple_thr(); void Build_hash_table_multiple_thr(); void Overlap_calculate_multipe_thr(); diff --git a/Makefile b/Makefile index 6925f23..c5b22d0 100644 --- a/Makefile +++ b/Makefile @@ -52,10 +52,11 @@ POA.o: kvec.h kdq.h CommandLines.h htab.h Correct.h Levenshtein_distance.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: CommandLines.h htab.h -hist.o: htab.h -htab.o: kthread.h khashl.h kseq.h htab.h CommandLines.h +hist.o: htab.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h +htab.o: kthread.h khashl.h kseq.h ksort.h htab.h Process_Read.h Overlaps.h +htab.o: kvec.h kdq.h CommandLines.h kthread.o: kthread.h main.o: CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h main.o: Assembly.h Levenshtein_distance.h htab.h -sketch.o: kvec.h htab.h -sys.o: htab.h +sketch.o: kvec.h htab.h Process_Read.h kseq.h Overlaps.h kdq.h CommandLines.h +sys.o: htab.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h diff --git a/Process_Read.cpp b/Process_Read.cpp index e461714..c24654f 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -590,7 +590,7 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID) }\ i++;}\ -void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ) +void ha_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 diff --git a/Process_Read.h b/Process_Read.h index e359328..ef6cbe9 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -18,11 +18,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.read_length[(ID)] -#define Get_NAME_LENGTH(R_INF, ID) (R_INF.name_index[(ID)+1] - R_INF.name_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_sperate[(ID)] -#define Get_NAME(R_INF, ID) R_INF.name + R_INF.name_index[(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)]) KSEQ_INIT(gzFile, gzread) @@ -139,8 +139,6 @@ typedef struct extern All_reads R_INF; -void malloc_All_reads(All_reads* r); - typedef struct { kseq_t* read; @@ -175,10 +173,12 @@ typedef struct void init_R_buffer(int thread_num); void init_All_reads(All_reads* r); +void malloc_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 ha_insert_read_len(All_reads *r, int read_len, int name_len); +void ha_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 recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID); @@ -196,6 +196,4 @@ 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 1de88d2..859d623 100644 --- a/htab.cpp +++ b/htab.cpp @@ -9,8 +9,6 @@ #include "kseq.h" #include "ksort.h" #include "htab.h" -#include "Process_Read.h" -#include "CommandLines.h" #define YAK_COUNTER_BITS 12 #define YAK_N_COUNTS (1<mz_buf[tid]; s->mz_buf[tid].n = 0; - ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->p->batch_offset + i, s->p->opt->is_HPC, b, s->p->flt_tab); + ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab); s->mz[i].n = s->mz[i].m = b->n; MALLOC(s->mz[i].a, b->n); memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t)); @@ -519,25 +522,42 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip st_data_t *s; CALLOC(s, 1); s->p = p; - while ((ret = kseq_read(p->ks)) >= 0) { - int l = p->ks->seq.l; - if (p->batch_offset + s->n_seq >= (1<<28) - 1) { - fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", (1<<28) - 1); - exit(1); + s->n_seq0 = p->n_seq; + if (p->rs_in && (p->flag & HAF_RS_READ)) { + } else { + while ((ret = kseq_read(p->ks)) >= 0) { + int l = p->ks->seq.l; + if (p->n_seq >= 1<<28) { + fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28); + exit(1); + } + if (p->rs_out) { + if (p->flag & HAF_RS_WRITE_LEN) { + assert(p->n_seq == p->rs_out->total_reads); + ha_insert_read_len(p->rs_out, l, p->ks->name.l); + } else if (p->flag & HAF_RS_WRITE_SEQ) { + int i, n_N; + for (i = n_N = 0; i < l; ++i) // count number of ambiguous bases + if (seq_nt4_table[(uint8_t)p->ks->seq.s[i]] >= 4) + ++n_N; + assert(l == (int)p->rs_out->read_length[p->n_seq]); + ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s, l, &p->rs_out->N_site[p->n_seq], n_N); + } + } + if (s->n_seq == s->m_seq) { + s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1); + REALLOC(s->len, s->m_seq); + REALLOC(s->seq, s->m_seq); + } + MALLOC(s->seq[s->n_seq], l); + memcpy(s->seq[s->n_seq], p->ks->seq.s, l); + s->len[s->n_seq++] = l; + ++p->n_seq; + s->sum_len += l; + s->nk += l - p->opt->k + 1; + if (s->sum_len >= p->opt->chunk_size) + break; } - p->n_base += l; - if (s->n_seq == s->m_seq) { - s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1); - REALLOC(s->len, s->m_seq); - REALLOC(s->seq, s->m_seq); - } - MALLOC(s->seq[s->n_seq], l); - memcpy(s->seq[s->n_seq], p->ks->seq.s, l); - s->len[s->n_seq++] = l; - s->sum_len += l; - s->nk += l - p->opt->k + 1; - if (s->sum_len >= p->opt->chunk_size) - break; } if (s->sum_len == 0) free(s); else return s; @@ -601,17 +621,16 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip } if (p->ct) p->ct->tot += n_ins; if (p->pt) p->pt->tot_pos += n_ins; - p->batch_offset += s->n_seq; free(s->buf); fprintf(stderr, "[M::%s::%.3f*%.2f] processed %ld sequences; %ld %s in the hash table\n", __func__, - yak_realtime(), yak_cputime() / yak_realtime(), (long)p->batch_offset, + yak_realtime(), yak_cputime() / yak_realtime(), (long)s->n_seq0 + s->n_seq, (long)(p->pt? p->pt->tot_pos : p->ct->tot), p->pt? "positions" : "distinct k-mers"); free(s); } return 0; } -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) +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, All_reads *rs) { pl_data_t pl; gzFile fp; @@ -620,6 +639,11 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt pl.ks = kseq_init(fp); pl.flt_tab = flt_tab; pl.opt = opt; + pl.flag = flag; + if (flag & (HAF_RS_WRITE_LEN|HAF_RS_WRITE_SEQ)) + pl.rs_out = rs; + else if (flag & HAF_RS_READ) + pl.rs_in = rs; if (p0) { pl.pt = p0, pl.create_new = 0; assert(p0->k == opt->k && p0->pre == opt->pre); @@ -636,28 +660,26 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt return pl.ct; } -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) +ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const void *flt_tab, All_reads *rs) { int i; - ha_ct_t *h = 0; - for (i = 0; i < n_fn; ++i) - 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, int flag, ha_pt_t *p0, const void *flt_tab) -{ yak_copt_t opt; - ha_ct_t *h; + ha_ct_t *h = 0; + assert(!(flag & HAF_RS_WRITE_LEN) || !(flag & HAF_RS_WRITE_SEQ)); // not both + if (flag & HAF_RS_WRITE_LEN) + init_All_reads(rs); + else if (flag & HAF_RS_WRITE_SEQ) + malloc_All_reads(rs); yak_copt_init(&opt); opt.k = asm_opt->k_mer_length; opt.is_HPC = !asm_opt->no_HPC; 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, flag, p0, asm_opt->num_reads, asm_opt->read_file_names, flt_tab); + for (i = 0; i < asm_opt->num_reads; ++i) + h = yak_count(&opt, asm_opt->read_file_names[i], flag, p0, h, flt_tab, rs); + if (h && opt.bf_shift > 0) + ha_ct_destroy_bf(h); return h; } @@ -704,13 +726,14 @@ void ha_ft_destroy(void *h) * High-level interfaces * *************************/ -void *ha_gen_flt_tab(const hifiasm_opt_t *asm_opt) +void *ha_gen_flt_tab(const hifiasm_opt_t *asm_opt, All_reads *rs) { yak_ft_t *flt_tab; int64_t cnt[YAK_N_COUNTS]; int peak_hom, peak_het, cutoff; ha_ct_t *h; - h = ha_count(asm_opt, HAF_COUNT_ALL, NULL, NULL); + h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN, NULL, NULL, rs); + if (rs) fprintf(stderr, "%ld,%ld,%ld\n", (long)rs->total_reads, (long)rs->total_reads_bases, (long)rs->total_name_length); 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); @@ -724,13 +747,13 @@ void *ha_gen_flt_tab(const hifiasm_opt_t *asm_opt) return (void*)flt_tab; } -void *ha_gen_mzidx(const hifiasm_opt_t *asm_opt, const void *flt_tab) +void *ha_gen_mzidx(const hifiasm_opt_t *asm_opt, const void *flt_tab, All_reads *rs) { int64_t cnt[YAK_N_COUNTS], tot_cnt; int peak_hom, peak_het, i; ha_ct_t *ct; ha_pt_t *pt; - ct = ha_count(asm_opt, HAF_COUNT_EXACT, NULL, flt_tab); + ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_RS_WRITE_SEQ, NULL, flt_tab, rs); 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); @@ -740,7 +763,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, HAF_COUNT_EXACT, pt, flt_tab); + ha_count(asm_opt, HAF_COUNT_EXACT, pt, flt_tab, rs); 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__, diff --git a/htab.h b/htab.h index bb07a98..eac0c98 100644 --- a/htab.h +++ b/htab.h @@ -2,6 +2,8 @@ #define __HA_HTAB_H__ #define __STDC_LIMIT_MACROS #include +#include "Process_Read.h" +#include "CommandLines.h" typedef struct { uint64_t x; @@ -16,12 +18,14 @@ typedef struct { uint32_t n, m; ha_mz1_t *a; } ha_mz1_v; extern const unsigned char seq_nt4_table[256]; +void *ha_gen_flt_tab(const hifiasm_opt_t *asm_opt, All_reads *rs); +void *ha_gen_mzidx(const hifiasm_opt_t *asm_opt, const void *flt_tab, All_reads *rs); +void trio_partition(void); + int ha_ft_isflt(const void *hh, uint64_t y); void ha_ft_destroy(void *h); void ha_idx_destroy(void *h); -void trio_partition(void); - double yak_cputime(void); void yak_reset_realtime(void); double yak_realtime(void); diff --git a/main.cpp b/main.cpp index e2a5cdd..d1acb21 100644 --- a/main.cpp +++ b/main.cpp @@ -17,13 +17,14 @@ int main(int argc, char *argv[]) if (!CommandLine_process(argc, argv, &asm_opt)) return 1; yak_reset_realtime(); - flt_tab = ha_gen_flt_tab(&asm_opt); - idx = ha_gen_mzidx(&asm_opt, flt_tab); + flt_tab = ha_gen_flt_tab(&asm_opt, &R_INF); + idx = ha_gen_mzidx(&asm_opt, flt_tab, &R_INF); ha_idx_destroy(idx); ha_ft_destroy(flt_tab); if (0) { Correct_Reads(asm_opt.number_of_round); } + destory_All_reads(&R_INF); destory_opt(&asm_opt); fprintf(stderr, "[M::%s] Version: %s\n", __func__, "dummy"); fprintf(stderr, "[M::%s] CMD:", __func__);