diff --git a/CommandLines.cpp b/CommandLines.cpp index 04821b5..81d7f4a 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -158,8 +158,11 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->thread_num = 1; asm_opt->k_mer_length = 51; asm_opt->hic_mer_length = 31; + asm_opt->ul_mer_length = 19; asm_opt->mz_win = 51; + asm_opt->ul_mz_win = 19; asm_opt->mz_rewin = 1000; + asm_opt->ul_mz_rewin = 360; asm_opt->mz_sample_dist = 500; asm_opt->bf_shift = 37; asm_opt->max_kmer_cnt = 2000; diff --git a/CommandLines.h b/CommandLines.h index 2f64e03..8a15dc9 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -44,9 +44,12 @@ typedef struct { int thread_num; int k_mer_length; int hic_mer_length; + int ul_mer_length; int bub_mer_length; int mz_win; + int ul_mz_win; int mz_rewin; + int ul_mz_rewin; int mz_sample_dist; int bf_shift; int max_kmer_cnt; diff --git a/Makefile b/Makefile index 030ab95..3d85a74 100644 --- a/Makefile +++ b/Makefile @@ -7,7 +7,7 @@ INCLUDES= OBJS= CommandLines.o Process_Read.o Assembly.o Hash_Table.o \ POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \ htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o rcut.o horder.o \ - tovlp.o + tovlp.o inter.o EXE= hifiasm LIBS= -lz -lpthread -lm @@ -76,3 +76,4 @@ hic.o: hic.h rcut.o: rcut.h horder.o: horder.h tovlp.o: tovlp.h +inter.o: inter.h diff --git a/anchor.cpp b/anchor.cpp index 08b9c31..83e03d2 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -76,7 +76,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg rlen = Get_READ_LENGTH(R_INF, rid); // read length // get the list of anchors - ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin); + ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0); // minimizer of queried read if (ab->mz.m > ab->old_mz_m) { ab->old_mz_m = ab->mz.m; diff --git a/htab.cpp b/htab.cpp index 0a01a65..516c971 100644 --- a/htab.cpp +++ b/htab.cpp @@ -38,6 +38,75 @@ void *ha_flt_tab_hp; ha_pt_t *ha_idx_hp; void *ha_ct_table; +#define MZ_FUNC_INIT(sf, HType) \ +static inline void sf##_init_kuf(pl_data_t *p, st_data_t *s){\ + int i, n_pre = 1<opt->pre, m;\ + /**allocate the k-mer buffer**/\ + CALLOC(s->buf, n_pre);\ + m = (int)(s->nk * 1.2 / n_pre) + 1;\ + /**pre-allocate memory for each of 4096 buffer**/\ + for (i = 0; i < n_pre; ++i) {\ + s->buf[i].m = m;\ + /**for 0-th counting, p->pt = NULL**/\ + if (p->pt && !(p->flag&HAF_COUNT_REFINE)) MALLOC(s->buf[i].b_##sf, m);\ + else MALLOC(s->buf[i].a, m);\ + }\ +}\ +static inline void sf##_destory_kuf(pl_data_t *p, st_data_t *s, int n){\ + int i;\ + uint64_t n_ins = 0;\ + /**n_ins is number of distinct k-mers**/\ + for (i = 0; i < n; ++i) {\ + n_ins += s->buf[i].n_ins;\ + if (p->pt && !(p->flag&HAF_COUNT_REFINE)) free(s->buf[i].b_##sf);\ + else free(s->buf[i].a);\ + }\ + if (p->ct) p->ct->tot += n_ins, p->ct->bs += s->sum_len;\ + if (p->pt) p->pt->tot_pos += n_ins;\ + free(s->buf);\ + /**#if 0\ + fprintf(stderr, "[M::%s::%.3f*%.2f] processed %ld sequences; %ld %s in the hash table\n", __func__,\ + yak_realtime(), yak_cpu_usage(), (long)s->n_seq0 + s->n_seq,\ + (long)(p->pt? p->pt->tot_pos : p->ct->tot), p->pt? "positions" : "distinct k-mers");\ + #endif**/\ + free(s);\ +}\ +static inline void sf##_pt_insert_buf(ch_buf_t *buf, int p, const HType *y){\ + /**assign minimizer to one of 4096 bins by low 12 bits**/\ + int pre = y->x & ((1<n == b->m) {\ + b->m = b->m < 8? 8 : b->m + (b->m>>1);\ + REALLOC(b->b_##sf, b->m);\ + }\ + b->b_##sf[b->n++] = *y;\ +}\ +static inline void sf##_mselect(pl_data_t *p, st_data_t *s){\ + int i; uint32_t j;\ + /**s->n_seq is how many reads at this buffer**/\ + /**s->mz && s->mz_buf are lists of minimzer vectors**/\ + CALLOC(s->sf, s->n_seq), CALLOC(s->sf##_buf, p->opt->n_thread), CALLOC(s->mt, p->opt->n_thread);\ + /**calculate minimzers for each read, each read corresponds to one thread**/\ + kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq);\ + for (i = 0; i < p->opt->n_thread; ++i) free(s->mt[i].a), free(s->sf##_buf[i].a);\ + free(s->mt), free(s->sf##_buf);\ + /**insert minimizers**/\ + if (p->pt && !(p->flag&HAF_COUNT_REFINE)) {/**insert whole minimizer**/\ + for (i = 0; i < s->n_seq; ++i)\ + for (j = 0; j < s->sf[i].n; ++j)\ + sf##_pt_insert_buf(s->buf, p->opt->pre, &s->sf[i].a[j]);\ + } else {/**just insert the hash key of minimizer**/\ + for (i = 0; i < s->n_seq; ++i)\ + for (j = 0; j < s->sf[i].n; ++j)\ + ct_insert_buf(s->buf, p->opt->pre, s->sf[i].a[j].x);\ + }\ + for (i = 0; i < s->n_seq; ++i) {\ + p->n_mz += s->sf[i].n;\ + free(s->sf[i].a);\ + if (!p->is_store) free(s->seq[i]);\ + }\ + free(s->sf);} + /*************************** * Yak specific parameters * ***************************/ @@ -303,6 +372,7 @@ typedef struct { yak_pt_t *h; uint64_t n; ha_idxpos_t *a; + ha_idxposl_t *al; } ha_pt1_t; struct ha_pt_s { @@ -314,6 +384,7 @@ struct ha_pt_s { typedef struct { const ha_ct_t *ct; ha_pt_t *pt; + int is_l; } pt_gen_aux_t; @@ -370,10 +441,11 @@ static void worker_pt_gen(void *data, long i, int tid) // callback for kt_for() } yak_ct_destroy(g); a->ct->h[i].h = 0; - CALLOC(b->a, b->n); + if(a->is_l) CALLOC(b->al, b->n); + else CALLOC(b->a, b->n); } -ha_pt_t *ha_pt_gen(ha_ct_t *ct, int n_thread) +ha_pt_t *ha_pt_gen(ha_ct_t *ct, int n_thread, int is_l) { pt_gen_aux_t a; int i; @@ -386,7 +458,7 @@ ha_pt_t *ha_pt_gen(ha_ct_t *ct, int n_thread) pt->h[i].h = yak_pt_init(); yak_pt_resize(pt->h[i].h, kh_size(ct->h[i].h)); } - a.ct = ct, a.pt = pt; + a.ct = ct, a.pt = pt, a.is_l = is_l; kt_for(n_thread, worker_pt_gen, &a, 1<pre); free(ct->h); free(ct); return pt; @@ -534,7 +606,8 @@ typedef struct { int n, m; uint64_t n_ins; uint64_t *a; - ha_mz1_t *b; + ha_mz1_t *b_mz; + ha_mzl_t *b_mzl; } ch_buf_t; ///p = 12 @@ -552,18 +625,6 @@ static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k b->a[b->n++] = y; } -static inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y) -{ - ///assign minimizer to one of 4096 bins by low 12 bits - int pre = y->x & ((1<n == b->m) { - b->m = b->m < 8? 8 : b->m + (b->m>>1); - REALLOC(b->b, b->m); - } - b->b[b->n++] = *y; -} - ///buf is the read block, k is the k-mer length, p = 12, len is the read length, seq is the read static void count_seq_buf(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf { @@ -639,11 +700,13 @@ typedef struct { // data structure for each step in kt_pipeline() pl_data_t *p; uint64_t n_seq0; ///the start index of current buffer block at R_INF ///sum_len = total bases, nk = number of k-mers - int n_seq, m_seq, sum_len, nk; + int n_seq, m_seq, sum_len, nk, uq; int *len; char **seq; ha_mz1_v *mz_buf; ha_mz1_v *mz; + ha_mzl_v *mzl_buf; + ha_mzl_v *mzl; ch_buf_t *buf; st_mt_t *mt; } st_data_t; @@ -654,14 +717,8 @@ static void worker_for_insert(void *data, long i, int tid) // callback for kt_fo ch_buf_t *b = &s->buf[i]; if (s->p->pt) { - if(s->p->flag&HAF_COUNT_REFINE) - { - b->n_ins += ha_pt_cnt_insert_list(s->p->pt, b->n, b->a); - } - else - { - b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b); - } + if(s->p->flag&HAF_COUNT_REFINE) b->n_ins += ha_pt_cnt_insert_list(s->p->pt, b->n, b->a); + else b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b_mz); } else///for 0-th count, go into here { @@ -676,12 +733,14 @@ static void worker_for_mz(void *data, long i, int tid) ha_mz1_v *b = &s->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->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, asm_opt.mz_sample_dist, 0, 0, - (s->p->pt&&(s->p->flag&HAF_COUNT_REFINE))?s->p->pt:NULL, s->p->opt->min_rcnt, asm_opt.dp_min_len, asm_opt.dp_e, &(s->mt[tid]), asm_opt.mz_rewin); + (s->p->pt&&(s->p->flag&HAF_COUNT_REFINE))?s->p->pt:NULL, s->p->opt->min_rcnt, asm_opt.dp_min_len, asm_opt.dp_e, &(s->mt[tid]), asm_opt.mz_rewin, s->uq); 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)); } +MZ_FUNC_INIT(mz, ha_mz1_t) +MZ_FUNC_INIT(mzl, ha_mzl_t) static void *worker_count(void *data, int step, void *in) // callback for kt_pipeline() { pl_data_t *p = (pl_data_t*)data; @@ -715,7 +774,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip break; } } else if(p->us_in) { - ma_utg_t *u; + ma_utg_t *u; s->uq = 1; while (p->n_seq < p->us_in->n) { u = &(p->us_in->a[p->n_seq]); if (s->n_seq == s->m_seq) { @@ -778,22 +837,13 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip } else if (step == 1) { // step 2: extract k-mers ///s is the block of reads st_data_t *s = (st_data_t*)in; - ///for 0-th counting, n_pre = 4096 - int i, n_pre = 1<opt->pre, m; - // allocate the k-mer buffer - CALLOC(s->buf, n_pre); - m = (int)(s->nk * 1.2 / n_pre) + 1; - //pre-allocate memory for each of 4096 buffer - for (i = 0; i < n_pre; ++i) { - s->buf[i].m = m; - ///for 0-th counting, p->pt = NULL - if (p->pt && !(p->flag&HAF_COUNT_REFINE)) MALLOC(s->buf[i].b, m); - else MALLOC(s->buf[i].a, m); - } + if(p->us_in) mzl_init_kuf(p, s); + else mz_init_kuf(p, s); // fill the buffer ///for 0-th counting, p->opt->w == 1 if (p->opt->w == 1) { // enumerate all k-mers ///scan all reads + int i; for (i = 0; i < s->n_seq; ++i) { if (p->opt->is_HPC) count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]); @@ -802,39 +852,8 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip if (!p->is_store) free(s->seq[i]); } } else { // minimizers only - uint32_t j; - // compute minimizers - // s->n_seq is how many reads at this buffer - // s->mz && s->mz_buf are lists of minimzer vectors - CALLOC(s->mz, s->n_seq); - CALLOC(s->mz_buf, p->opt->n_thread); - CALLOC(s->mt, p->opt->n_thread); - ///calculate minimzers for each read, each read corresponds to one thread - kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq); - - for (i = 0; i < p->opt->n_thread; ++i) - { - free(s->mt[i].a); - free(s->mz_buf[i].a); - } - free(s->mt); - free(s->mz_buf); - // insert minimizers - if (p->pt && !(p->flag&HAF_COUNT_REFINE)) {///insert whole minimizer - for (i = 0; i < s->n_seq; ++i) - for (j = 0; j < s->mz[i].n; ++j) - pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]); - } else {///just insert the hash key of minimizer - for (i = 0; i < s->n_seq; ++i) - for (j = 0; j < s->mz[i].n; ++j) - ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x); - } - for (i = 0; i < s->n_seq; ++i) { - p->n_mz += s->mz[i].n; - free(s->mz[i].a); - if (!p->is_store) free(s->seq[i]); - } - free(s->mz); + if(p->us_in) mzl_mselect(p, s); + else mz_mselect(p, s); } ///just clean seq free(s->seq); free(s->len); @@ -842,25 +861,10 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip return s; } else if (step == 2) { // step 3: insert k-mers to hash table st_data_t *s = (st_data_t*)in; - int i, n = 1<opt->pre; - uint64_t n_ins = 0; ///for 0-th counting, p->pt = NULL - kt_for(p->opt->n_thread, worker_for_insert, s, n); - ///n_ins is number of distinct k-mers - for (i = 0; i < n; ++i) { - n_ins += s->buf[i].n_ins; - if (p->pt && !(p->flag&HAF_COUNT_REFINE)) free(s->buf[i].b); - else free(s->buf[i].a); - } - if (p->ct) p->ct->tot += n_ins, p->ct->bs += s->sum_len; - if (p->pt) p->pt->tot_pos += n_ins; - free(s->buf); - #if 0 - fprintf(stderr, "[M::%s::%.3f*%.2f] processed %ld sequences; %ld %s in the hash table\n", __func__, - yak_realtime(), yak_cpu_usage(), (long)s->n_seq0 + s->n_seq, - (long)(p->pt? p->pt->tot_pos : p->ct->tot), p->pt? "positions" : "distinct k-mers"); - #endif - free(s); + kt_for(p->opt->n_thread, worker_for_insert, s, 1<opt->pre); + if(p->us_in) mzl_destory_kuf(p, s, 1<opt->pre); + else mz_destory_kuf(p, s, 1<opt->pre); } return 0; } @@ -975,11 +979,11 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi malloc_All_reads(rs); } yak_copt_init(&opt); - opt.k = asm_opt->k_mer_length; + opt.k = us? asm_opt->ul_mer_length:asm_opt->k_mer_length; ///always 0 opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC); ///for ft-counting, shoud be 1 - opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win; + opt.w = flag & HAF_COUNT_ALL? 1 : (us? asm_opt->ul_mz_win:asm_opt->mz_win); ///for ft-counting, shoud be 37 ///for ha_pt_gen, shoud be 0 opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift; @@ -987,7 +991,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi opt.adaLen = (keep_adapter? asm_opt->adapterLen : 0); opt.min_rcnt = (low_freq?*low_freq:-1); ///asm_opt->num_reads is the number of fastq files - for (i = n_bs = 0; i < asm_opt->num_reads; ++i){ + for (i = n_bs = 0; i < (us?1:asm_opt->num_reads); ++i){ h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq); if(h) n_bs += h->bs; } @@ -1076,50 +1080,29 @@ void debug_ct_index(void* q_ct_idx, void* r_ct_idx) void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int hap_n) { yak_ft_t *flt_tab; - int64_t cnt[YAK_N_COUNTS]; - int cutoff = hap_n + 1; ha_ct_t *h; + ///HAF_COUNT_EXACT ---> no bf; HAF_COUNT_ALL ---> no minimizer h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, NULL, NULL, NULL, us, 0, NULL); - - ha_ct_hist(h, cnt, asm_opt->thread_num); - print_hist_lines(YAK_N_COUNTS, 1, cnt); - - ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num); + ha_ct_shrink(h, 1, YAK_MAX_COUNT-1, asm_opt->thread_num); flt_tab = gen_hh(h, asm_opt->max_kmer_cnt); ha_ct_destroy(h); - fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> filtered out %ld k-mers occurring %d or more times\n", __func__, - yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb(), (long)kh_size(flt_tab), cutoff); return (void*)flt_tab; } ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int hap_n) { - int64_t cnt[YAK_N_COUNTS], tot_cnt; - int i; ha_ct_t *ct; ha_pt_t *pt; ///HAF_COUNT_EXACT: no bf ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, NULL, flt_tab, NULL, us, 0, NULL); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__, yak_realtime(), yak_cpu_usage(), (long)ct->tot); + ///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen + ha_ct_shrink(ct, 1, YAK_MAX_COUNT - 1, asm_opt->thread_num); - ha_ct_hist(ct, cnt, asm_opt->thread_num); - print_hist_lines(YAK_N_COUNTS, 1, cnt); - - ///here ha_ct_shrink is mostly used to remove k-mer appearing only 1 time - if (flt_tab == 0) { - ha_ct_shrink(ct, 2, hap_n, asm_opt->thread_num); - for (i = 2, tot_cnt = 0; i <= hap_n; ++i) tot_cnt += cnt[i] * i; - } else { - ///Note: here is just to remove minimizer appearing YAK_MAX_COUNT times - ///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen - 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); + pt = ha_pt_gen(ct, asm_opt->thread_num, 1); ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, pt, flt_tab, NULL, us, 0, NULL); - 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__, yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos); @@ -1219,7 +1202,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f if(!(asm_opt->flag & HA_F_FAST)) { fprintf(stderr, "[M::%s::] counting in normal mode\n", __func__); - pt = ha_pt_gen(ct, asm_opt->thread_num); + pt = ha_pt_gen(ct, asm_opt->thread_num, 0); ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, pt, flt_tab, rs, NULL, 1, NULL); assert((uint64_t)tot_cnt == pt->tot_pos); } diff --git a/htab.h b/htab.h index cd46dea..28be901 100644 --- a/htab.h +++ b/htab.h @@ -25,16 +25,16 @@ typedef struct { uint32_t n, m; ha_mz1_t *a; } ha_mz1_v; typedef struct { uint64_t x; ///x is the hash key - ///rid is the read id, pos is the end pos of this minimizer, rev is the direction - ///span is the length of this k-mer. For non-HPC k-mer, span may not be equal to k - uint64_t rid:30, pos:34; - uint16_t rev:1, span:15; + uint64_t rid:31, rev:1; + uint32_t pos; + uint8_t span; } ha_mzl_t; typedef struct { - uint64_t rid:30, pos:34; - uint16_t rev:1, span:15; -} ha_mzl_idxpos_t; + uint64_t rid:31, rev:1; + uint32_t pos; + uint8_t span; +} ha_idxposl_t; typedef struct { uint32_t n, m; ha_mzl_t *a; } ha_mzl_v; @@ -101,7 +101,7 @@ double yak_cpu_usage(void); void ha_triobin(const hifiasm_opt_t *opt); -void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws); +void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique); int ha_analyze_count(int n_cnt, int start_cnt, int m_peak_hom, const int64_t *cnt, int *peak_het); int adj_m_peak_hom(int m_peak_hom, int max_i, int max2_i, int max3_i, int *peak_het); void print_hist_lines(int n_cnt, int start_cnt, const int64_t *cnt); diff --git a/inter.cpp b/inter.cpp new file mode 100644 index 0000000..b3c685a --- /dev/null +++ b/inter.cpp @@ -0,0 +1,25 @@ +#include +#include +#include +#include +#include "inter.h" +#include "Overlaps.h" +#include "CommandLines.h" +#include "htab.h" + +void uidx_build(ma_ug_t *ug, int hap_n) +{ + int flag = asm_opt.flag; + asm_opt.flag |= HA_F_NO_HPC; + ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), hap_n); + ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), hap_n); + + + + + + + asm_opt.flag = flag; + ha_ft_destroy(ha_flt_tab); + ha_pt_destroy(ha_idx); +} diff --git a/inter.h b/inter.h new file mode 100644 index 0000000..b8efb52 --- /dev/null +++ b/inter.h @@ -0,0 +1,5 @@ +#ifndef __INTER__ +#define __INTER__ + + +#endif diff --git a/sketch.cpp b/sketch.cpp index c9a92ab..a7b08b7 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -582,7 +582,7 @@ void debug_pl(const char *str, int len, int w, int k, int is_hpc, ha_mz1_v *p, c * @param is_hpc homopolymer-compressed or not * @param p minimizers */ -void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws) +void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique) { ///in default, w = 51, k = 51, is_hpc = 1 /** uint64_t x; @@ -648,6 +648,10 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]); cnt = hf? ha_ft_cnt(hf, y) : 0; filtered = (cnt >= 1<<28); + if(is_unique){ + filtered = (cnt == 0); + cnt = (cnt == 1? 0:cnt); + } if (dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i)); if (!filtered) info.x = y, info.rid = cnt, info.pos = i, info.rev = z, info.span = kmer_span; // initially ha_mz1_t::rid keeps the k-mer count if (k_flag != NULL) k_flag->a.a[i]++; @@ -711,128 +715,9 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, kv_push(ha_mz1_t, *p, min); kv_push(uint64_t, *mt, min_s); } // debug_pl(str, len, w, k, is_hpc, p, hf, mt); - // if (sample_dist > w) select_mz(p, len, MAX_HIGH_OCC, dp_min_len); - select_mz_h(p, mt, len, sample_dist, ws, k, tl); + if (sample_dist > w) select_mz_h(p, mt, len, sample_dist, ws, k, tl); if (dp_min_len > 0 && pt && mt) refine_sketch(p, pt, len, dp_min_len, dp_e, min_freq, mt); - for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts - p->a[i].rid = rid; -} - -void ha_sketch_worse(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt) -{ ///in default, w = 51, k = 51, is_hpc = 1 - /** - uint64_t x; - uint64_t rid:28, pos:27, rev:1, span:8; - **/ - extern void *ha_ct_table; - static const ha_mz1_t dummy = { UINT64_MAX, (1<<28) - 1, 0, 0 }; - uint64_t shift1 = k - 1, mask = (1ULL< 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63)); - if (dbg_ct != NULL) dbg_ct->a.n = 0; - if (k_flag != NULL) { - kv_resize(uint8_t, k_flag->a, (uint64_t)len); - k_flag->a.n = len; - memset(k_flag->a.a, 0, k_flag->a.n); - } - - memset(buf, 0xff, w * sizeof(ha_mz1_t)); - memset(&tq, 0, sizeof(tiny_queue_t)); - ///len/w is the evaluated minimizer numbers - kv_resize(ha_mz1_t, *p, p->n + len/w); - - for (i = l = buf_pos = min_pos = 0; i < len; ++i) { - int c = seq_nt4_table[(uint8_t)str[i]]; - ha_mz1_t info = dummy; - if (c < 4) { // not an ambiguous base - int z; - if (is_hpc) { - int skip_len = 1; - if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) { - for (skip_len = 2; i + skip_len < len; ++skip_len) - if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c) - break; - i += skip_len - 1; // put $i at the end of the current homopolymer run - } - tq_push(&tq, skip_len); - kmer_span += skip_len; - ///how many bases that are covered by this HPC k-mer - ///kmer_span includes at most k HPC elements - if (tq.count > k) kmer_span -= tq_shift(&tq); - } else kmer_span = l + 1 < k? l + 1 : k; - ///kmer_span should be used for HPC k-mer - ///non-HPC k-mer, kmer_span should be k - ///kmer_span is used to calculate anchor pos on reverse complementary strand - - if (k_flag != NULL) k_flag->a.a[i] = 1;///lable all useful base, which are not ignored by HPC - - kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer - kmer[1] = (kmer[1] << 1 | (c>>1)) & mask; - kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer - kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1; - if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand - z = kmer[1] < kmer[3]? 0 : 1; // strand - ++l; - if (l >= k && kmer_span < 256) { - uint64_t y; - int32_t cnt, filtered; - y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]); - cnt = hf? ha_ft_cnt(hf, y) : 0; - filtered = (cnt >= 1<<28); - if (dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i)); - if (!filtered) info.x = y, info.rid = pt? ha_pt_cnt(pt, y):cnt, info.pos = i, info.rev = z, info.span = kmer_span; // initially ha_mz1_t::rid keeps the k-mer count - if (k_flag != NULL) k_flag->a.a[i]++; - if (k_flag != NULL && filtered > 0) k_flag->a.a[i]++; - } - } else l = 0, tq.count = tq.front = 0, kmer_span = 0; - - //for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i - //i is the real base iterator, while l is the HPC base iterator - //only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX) - //but even if l < k, infor is still stored into buf - buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below - if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet - for (j = buf_pos + 1; j < w; ++j) - if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]); - for (j = 0; j < buf_pos; ++j) - if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]); - } - /** - * There are three cases: - * 1. info.x <= min.x, means info is a new minimizer - * 2. info.x > min.x, info is not a new minimizer - * (1) buf_pos != min_pos, do nothing - * (2) buf_pos == min_pos, means current minimizer has moved outside the window - * **/ - ///three cases: 1. - if (info.x <= min.x) { // a new minimum; then write the old min - if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); - min = info, min_pos = buf_pos; - } else if (buf_pos == min_pos) { // old min has moved outside the window - if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); - ///buf_pos == min_pos, means current minimizer has moved outside the window - ///so for now we need to find a new minimizer at the current window (w k-mers) - for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers - if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer - for (j = 0; j <= buf_pos; ++j) - if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j; - - if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers - for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted - if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]); - for (j = 0; j <= buf_pos; ++j) - if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]); - } - } - if (++buf_pos == w) buf_pos = 0; - } - if (min.x != UINT64_MAX) - kv_push(ha_mz1_t, *p, min); - if (dp_min_len > 0 && pt && mt) refine_sketch(p, pt, len, dp_min_len, dp_e, min_freq, mt); - // if (sample_dist > w) select_mz(p, len, MAX_HIGH_OCC); for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts p->a[i].rid = rid; } +