diff --git a/Assembly.cpp b/Assembly.cpp index 588dcd6..3603a27 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -425,6 +425,7 @@ typedef struct { overlap_region_alloc olist; overlap_region_alloc olist_hp; ha_abuf_t *ab; + ha_abufl_t *abl; // error correction related buffers int64_t num_read_base, num_correct_base, num_recorrect_base; Cigar_record cigar1; @@ -441,7 +442,7 @@ typedef struct { st_mt_t sp; } ha_ovec_buf_t; -ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov) +ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov, int is_ug) { ha_ovec_buf_t *b; CALLOC(b, 1); @@ -456,7 +457,8 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov) kv_init(b->r_buf.a); kv_init(b->k_flag.a); kv_init(b->sp); - b->ab = ha_abuf_init(); + if(!is_ug) b->ab = ha_abuf_init(); + else b->abl = ha_abufl_init(); if (!b->is_final) { init_Cigar_record(&b->cigar1); init_Graph(&b->POA_Graph); @@ -476,6 +478,7 @@ void ha_ovec_destroy(ha_ovec_buf_t *b) destory_overlap_region_alloc(&b->olist); destory_overlap_region_alloc(&b->olist_hp); ha_abuf_destroy(b->ab); + ha_abufl_destroy(b->abl); destory_fake_cigar(&(b->tmp_region.f_cigar)); kv_destroy(b->b_buf.a); kv_destroy(b->r_buf.a); @@ -526,7 +529,8 @@ int64_t ha_ovec_mem(const ha_ovec_buf_t *b) mem_olist += r->boundary_cigars.size * sizeof(window_list); } - mem = ha_abuf_mem(b->ab) + mem_clist + mem_olist; + if(b->ab) mem = ha_abuf_mem(b->ab) + mem_clist + mem_olist; + if(b->abl) mem = ha_abufl_mem(b->abl) + mem_clist + mem_olist; if (!b->is_final) { mem += sizeof(Cigar_record) + b->cigar1.lost_base_size + b->cigar1.size * 4; mem += sizeof(Correct_dumy) + b->correct.size * 8; @@ -852,7 +856,7 @@ void ha_overlap_and_correct(int round) // overlap and correct reads CALLOC(b, asm_opt.thread_num); for (i = 0; i < asm_opt.thread_num; ++i) - b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1)); + b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1),0); if(ha_idx) hom_cov = asm_opt.hom_cov; if(ha_idx == NULL) ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, 0, &R_INF, &hom_cov, &het_cov); // build the index ///debug_adapter(&asm_opt, &R_INF); @@ -1588,7 +1592,7 @@ void ha_overlap_final(void) CALLOC(b, asm_opt.thread_num); for (i = 0; i < asm_opt.thread_num; ++i) - b[i] = ha_ovec_init(asm_opt.flag & HA_F_HIGH_HET, 1);///b[i] = ha_ovec_init(1, 1); + b[i] = ha_ovec_init(asm_opt.flag & HA_F_HIGH_HET, 1,0);///b[i] = ha_ovec_init(1, 1); ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, 0, &R_INF, &hom_cov, &het_cov); // build the index if(asm_opt.flag & HA_F_HIGH_HET) { @@ -1644,7 +1648,7 @@ void ug_idx_build(ma_ug_t *ug, int hap_n) CALLOC(b, asm_opt.thread_num); for (i = 0; i < asm_opt.thread_num; ++i) { - b[i] = ha_ovec_init(1, 1); + b[i] = ha_ovec_init(1, 1, 0); b[i]->ua = &(ug->u); } diff --git a/CommandLines.h b/CommandLines.h index fd403ba..5490e8f 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.16.1-r374" +#define HA_VERSION "0.16.1-r375" #define VERBOSE 0 diff --git a/Overlaps.h b/Overlaps.h index 8f215c5..9430e15 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -190,7 +190,7 @@ void print_gfa(asg_t *g); typedef struct { size_t n, m; uint64_t *a; } asg64_v; -typedef struct { size_t n, m; ma_utg_t *a; int h;} ma_utg_v; +typedef struct { size_t n, m; ma_utg_t *a;} ma_utg_v; typedef struct { ma_utg_v u; diff --git a/anchor.cpp b/anchor.cpp index f10d607..3a07f95 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -29,6 +29,11 @@ typedef struct { const ha_idxpos_t *a; } seed1_t; +typedef struct { + int n, cnt; + const ha_idxposl_t *a; +} seedl_t; + struct ha_abuf_s { uint64_t n_a, m_a;///number of anchors (seed positions) uint32_t old_mz_m;///number of seeds @@ -37,20 +42,23 @@ struct ha_abuf_s { anchor1_t *a; }; -ha_abuf_t *ha_abuf_init(void) -{ - return (ha_abuf_t*)calloc(1, sizeof(ha_abuf_t)); +struct ha_abufl_s { + uint64_t n_a, m_a;///number of anchors (seed positions) + uint32_t old_mz_m;///number of seeds + ha_mzl_v mz; + seedl_t *seed; + anchor1_t *a; +}; + +#define HA_ABUF_INIT(HType, MZType, SDType, sf) \ +HType *sf##_init(void){return (HType*)calloc(1, sizeof(HType));}\ +void sf##_destroy(HType *ab){if(ab){free(ab->seed); free(ab->a); free(ab->mz.a); free(ab);}}\ +uint64_t sf##_mem(const HType *ab){\ + return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(MZType) + sizeof(SDType)) + sizeof(HType);\ } -void ha_abuf_destroy(ha_abuf_t *ab) -{ - free(ab->seed); free(ab->a); free(ab->mz.a); free(ab); -} - -uint64_t ha_abuf_mem(const ha_abuf_t *ab) -{ - return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(ha_mz1_t) + sizeof(seed1_t)) + sizeof(ha_abuf_t); -} +HA_ABUF_INIT(ha_abuf_s, ha_mz1_t, seed1_t, ha_abuf) +HA_ABUF_INIT(ha_abufl_s, ha_mzl_t, seedl_t, ha_abufl) int ha_ov_type(const overlap_region *r, uint32_t len) { @@ -257,122 +265,248 @@ void calculate_ug_chaining(Candidates_list* candidates, overlap_region_alloc* ov } } +/** +void ha_get_inter_candidates(ha_abufl_t *ab, int64_t uid, ma_utg_v *ua, overlap_region_alloc *ovlp, +Candidates_list *cl, double bw_thres, int max_n_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, +void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, +double chain_match_rate, uint32_t uk, uint32_t uw, uint32_t is_hpc, uint32_t h_occ) +{ + + uint32_t i, rlen; + uint64_t k, l; + ma_utg_t *u = &(ua->a[uid]); + // prepare + clear_Candidates_list(cl); + clear_overlap_region_alloc(ovlp); + ab->mz.n = 0, ab->n_a = 0; + + mz2_ha_sketch(u->s, u->len, uw, uk, uid, is_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, 1); + // minimizer of queried read + if (ab->mz.m > ab->old_mz_m) { + ab->old_mz_m = ab->mz.m; + REALLOC(ab->seed, ab->old_mz_m); + } + for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) { + int n; + ab->seed[i].a = ha_ptl_get(ha_idx, ab->mz.a[i].x, &n); + ab->seed[i].n = n; + ab->seed[i].cnt = ha_ft_cnt(ha_flt_tab, ab->mz.a[i].x); + ab->n_a += n; + } + if (ab->n_a > ab->m_a) { + ab->m_a = ab->n_a; + kroundup64(ab->m_a); + REALLOC(ab->a, ab->m_a); + } + for (i = 0, k = 0; i < ab->mz.n; ++i) { + int j; + ///z is one of the minimizer + ha_mzl_t *z = &ab->mz.a[i]; + seedl_t *s = &ab->seed[i]; + for (j = 0; j < s->n; ++j) { + const ha_idxposl_t *y = &s->a[j]; + anchor1_t *an = &ab->a[k++]; + uint8_t rev = z->rev == y->rev? 0 : 1; + an->other_off = y->pos; + an->self_off = rev? u->len - 1 - (z->pos + 1 - z->span) : z->pos; + an->cnt = s->cnt; + an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off; + } + } + + // sort anchors + radix_sort_ha_an1(ab->a, ab->a + ab->n_a); + for (k = 1, l = 0; k <= ab->n_a; ++k) { + if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) { + if (k - l > 1) + radix_sort_ha_an2(ab->a + l, ab->a + k); + l = k; + } + } + + + // copy over to _cl_ + if (ab->m_a >= (uint64_t)cl->size) { + cl->size = ab->m_a; + REALLOC(cl->list, cl->size); + } + for (k = 0; k < ab->n_a; ++k) { + k_mer_hit *p = &cl->list[k]; + p->readID = ab->a[k].srt >> 33; + p->strand = ab->a[k].srt >> 32 & 1; + p->offset = ab->a[k].other_off; + p->self_offset = ab->a[k].self_off; + p->cnt = (ab->a[k].cnt == 1? 1 : 16); + } + cl->length = ab->n_a; + // calculate_overlap_region_by_chaining(cl, ovlp, chain_idx, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain, f_cigar); + // calculate_ug_chaining(cl, overlap_list, chain_idx, rid, ua, bw_thres, keep_whole_chain, f_cigar, ab->mz.n, chain_match_rate); + + + #if 0 + if (overlap_list->length > 0) { + fprintf(stderr, "B\t%ld\t%ld\t%d\n", (long)rid, (long)overlap_list->length, rlen); + for (int i = 0; i < (int)overlap_list->length; ++i) { + overlap_region *r = &overlap_list->list[i]; + fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand], + (int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed, ha_ov_type(r, rlen)); + } + } + #endif + + if ((int)ovlp->length > max_n_chain) { + int32_t w, n[4], s[4]; + n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0; + ks_introsort_or_ss(ovlp->length, ovlp->list); + for (i = 0; i < (uint32_t)ovlp->length; ++i) { + const overlap_region *r = &ovlp->list[i]; + w = ha_ov_type(r, rlen); + ++n[w]; + if ((int)n[w] == max_n_chain) s[w] = r->shared_seed; + } + if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) { + // n[0] = n[1] = n[2] = n[3] = 0; + for (i = 0, k = 0; i < (uint32_t)ovlp->length; ++i) { + overlap_region *r = &ovlp->list[i]; + w = ha_ov_type(r, rlen); + // ++n[w]; + // if (((int)n[w] <= max_n_chain) || (r->shared_seed >= s[w] && s[w] >= (asm_opt.k_mer_length<<1))) { + if (r->shared_seed >= s[w]) { + if ((uint32_t)k != i) { + overlap_region t; + t = ovlp->list[k]; + ovlp->list[k] = ovlp->list[i]; + ovlp->list[i] = t; + } + ++k; + } + } + ovlp->length = k; + } + } + ///ks_introsort_or_xs(overlap_list->length, overlap_list->list); +} +**/ + + void ha_get_ug_candidates(ha_abuf_t *ab, int64_t rid, ma_utg_t *u, ma_utg_v *ua, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, double chain_match_rate) { - uint32_t i; - uint64_t k, l; + uint32_t i; + uint64_t k, l; - // prepare + // prepare clear_Candidates_list(cl); clear_overlap_region_alloc(overlap_list); - ab->mz.n = 0, ab->n_a = 0; + ab->mz.n = 0, ab->n_a = 0; - // get the list of anchors - //should use the new version... - ///ha_sketch_query(u->s, u->len, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag, dbg_ct); - // minimizer of queried read - if (ab->mz.m > ab->old_mz_m) { - ab->old_mz_m = ab->mz.m; - REALLOC(ab->seed, ab->old_mz_m); - } - for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) { - int n; - ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &n); - ab->seed[i].n = n; - ab->n_a += n; - } - if (ab->n_a > ab->m_a) { - ab->m_a = ab->n_a; - kroundup64(ab->m_a); - REALLOC(ab->a, ab->m_a); - } - for (i = 0, k = 0; i < ab->mz.n; ++i) { - int j; - ///z is one of the minimizer - ha_mz1_t *z = &ab->mz.a[i]; - seed1_t *s = &ab->seed[i]; - for (j = 0; j < s->n; ++j) { - const ha_idxpos_t *y = &s->a[j]; - anchor1_t *an = &ab->a[k++]; - uint8_t rev = z->rev == y->rev? 0 : 1; - an->other_off = y->pos; - an->self_off = rev? u->len - 1 - (z->pos + 1 - z->span) : z->pos; - an->cnt = 1; - an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off; - } - } + // get the list of anchors + //should use the new version... + ///ha_sketch_query(u->s, u->len, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag, dbg_ct); + // minimizer of queried read + if (ab->mz.m > ab->old_mz_m) { + ab->old_mz_m = ab->mz.m; + REALLOC(ab->seed, ab->old_mz_m); + } + for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) { + int n; + ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &n); + ab->seed[i].n = n; + ab->n_a += n; + } + if (ab->n_a > ab->m_a) { + ab->m_a = ab->n_a; + kroundup64(ab->m_a); + REALLOC(ab->a, ab->m_a); + } + for (i = 0, k = 0; i < ab->mz.n; ++i) { + int j; + ///z is one of the minimizer + ha_mz1_t *z = &ab->mz.a[i]; + seed1_t *s = &ab->seed[i]; + for (j = 0; j < s->n; ++j) { + const ha_idxpos_t *y = &s->a[j]; + anchor1_t *an = &ab->a[k++]; + uint8_t rev = z->rev == y->rev? 0 : 1; + an->other_off = y->pos; + an->self_off = rev? u->len - 1 - (z->pos + 1 - z->span) : z->pos; + an->cnt = 1; + an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off; + } + } - // sort anchors - radix_sort_ha_an1(ab->a, ab->a + ab->n_a); - for (k = 1, l = 0; k <= ab->n_a; ++k) { - if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) { - if (k - l > 1) - radix_sort_ha_an2(ab->a + l, ab->a + k); - l = k; - } - } + // sort anchors + radix_sort_ha_an1(ab->a, ab->a + ab->n_a); + for (k = 1, l = 0; k <= ab->n_a; ++k) { + if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) { + if (k - l > 1) + radix_sort_ha_an2(ab->a + l, ab->a + k); + l = k; + } + } - // copy over to _cl_ - if (ab->m_a >= (uint64_t)cl->size) { - cl->size = ab->m_a; - REALLOC(cl->list, cl->size); - } - for (k = 0; k < ab->n_a; ++k) { - k_mer_hit *p = &cl->list[k]; - p->readID = ab->a[k].srt >> 33; - p->strand = ab->a[k].srt >> 32 & 1; - p->offset = ab->a[k].other_off; - p->self_offset = ab->a[k].self_off; - p->cnt = 1; - } - cl->length = ab->n_a; + // copy over to _cl_ + if (ab->m_a >= (uint64_t)cl->size) { + cl->size = ab->m_a; + REALLOC(cl->list, cl->size); + } + for (k = 0; k < ab->n_a; ++k) { + k_mer_hit *p = &cl->list[k]; + p->readID = ab->a[k].srt >> 33; + p->strand = ab->a[k].srt >> 32 & 1; + p->offset = ab->a[k].other_off; + p->self_offset = ab->a[k].self_off; + p->cnt = 1; + } + cl->length = ab->n_a; - calculate_ug_chaining(cl, overlap_list, chain_idx, rid, ua, bw_thres, keep_whole_chain, f_cigar, ab->mz.n, chain_match_rate); + calculate_ug_chaining(cl, overlap_list, chain_idx, rid, ua, bw_thres, keep_whole_chain, f_cigar, ab->mz.n, chain_match_rate); - #if 0 - if (overlap_list->length > 0) { - fprintf(stderr, "B\t%ld\t%ld\t%d\n", (long)rid, (long)overlap_list->length, rlen); - for (int i = 0; i < (int)overlap_list->length; ++i) { - overlap_region *r = &overlap_list->list[i]; - fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand], - (int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed, ha_ov_type(r, rlen)); - } - } - #endif + #if 0 + if (overlap_list->length > 0) { + fprintf(stderr, "B\t%ld\t%ld\t%d\n", (long)rid, (long)overlap_list->length, rlen); + for (int i = 0; i < (int)overlap_list->length; ++i) { + overlap_region *r = &overlap_list->list[i]; + fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand], + (int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed, ha_ov_type(r, rlen)); + } + } + #endif - if ((int)overlap_list->length > max_n_chain) { - int32_t w, n[4], s[4]; - n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0; - ks_introsort_or_ss(overlap_list->length, overlap_list->list); - for (i = 0; i < (uint32_t)overlap_list->length; ++i) { - const overlap_region *r = &overlap_list->list[i]; - w = ha_ov_type(r, u->len); - ++n[w]; - if ((int)n[w] == max_n_chain) s[w] = r->shared_seed; - } - if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) { - for (i = 0, k = 0; i < (uint32_t)overlap_list->length; ++i) { - overlap_region *r = &overlap_list->list[i]; - w = ha_ov_type(r, u->len); - if (r->shared_seed >= s[w]) { - if ((uint32_t)k != i) { - overlap_region t; - t = overlap_list->list[k]; - overlap_list->list[k] = overlap_list->list[i]; - overlap_list->list[i] = t; - } - ++k; - } - } - overlap_list->length = k; - } - } + if ((int)overlap_list->length > max_n_chain) { + int32_t w, n[4], s[4]; + n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0; + ks_introsort_or_ss(overlap_list->length, overlap_list->list); + for (i = 0; i < (uint32_t)overlap_list->length; ++i) { + const overlap_region *r = &overlap_list->list[i]; + w = ha_ov_type(r, u->len); + ++n[w]; + if ((int)n[w] == max_n_chain) s[w] = r->shared_seed; + } + if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) { + for (i = 0, k = 0; i < (uint32_t)overlap_list->length; ++i) { + overlap_region *r = &overlap_list->list[i]; + w = ha_ov_type(r, u->len); + if (r->shared_seed >= s[w]) { + if ((uint32_t)k != i) { + overlap_region t; + t = overlap_list->list[k]; + overlap_list->list[k] = overlap_list->list[i]; + overlap_list->list[i] = t; + } + ++k; + } + } + overlap_list->length = k; + } + } - ///ks_introsort_or_xs(overlap_list->length, overlap_list->list); + ///ks_introsort_or_xs(overlap_list->length, overlap_list->list); } + void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf) { uint64_t j = 0, inner_j = 0; diff --git a/horder.cpp b/horder.cpp index 53761f5..7fe2e55 100644 --- a/horder.cpp +++ b/horder.cpp @@ -4016,7 +4016,7 @@ void update_sc_lay(sc_lay_t *sl, h_covs *b) kv_init(*p); p->n = p->m = (cidx - pidx + 1)<<1; MALLOC(p->a, p->n); - mempcpy(p->a, s->a + (pidx<<1), sizeof(*(p->a))*p->n); + memcpy(p->a, s->a + (pidx<<1), sizeof(*(p->a))*p->n); pidx = cidx + 1; } @@ -4026,7 +4026,7 @@ void update_sc_lay(sc_lay_t *sl, h_covs *b) kv_init(*p); p->n = p->m = (cidx - pidx + 1)<<1; MALLOC(p->a, p->n); - mempcpy(p->a, s->a + (pidx<<1), sizeof(*(p->a))*p->n); + memcpy(p->a, s->a + (pidx<<1), sizeof(*(p->a))*p->n); free(s->a); s->n = s->m = 0; l = k; } diff --git a/htab.cpp b/htab.cpp index d6bec2a..02a6fef 100644 --- a/htab.cpp +++ b/htab.cpp @@ -504,7 +504,13 @@ void ha_pt_destroy(ha_pt_t *h) if (h == 0) return; for (i = 0; i < 1<pre; ++i) { yak_pt_destroy(h->h[i].h); - free(h->h[i].a); + if(h->h[i].a){ + free(h->h[i].a); h->h[i].a = NULL; + } + if(h->h[i].al){ + free(h->h[i].al); h->h[i].al = NULL; + } + } free(h->h); free(h); } @@ -520,6 +526,17 @@ const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n) return &g->a[kh_val(g->h, k)]; } +const ha_idxposl_t *ha_ptl_get(const ha_pt_t *h, uint64_t hash, int *n) +{ + khint_t k; + const ha_pt1_t *g = &h->h[hash & ((1ULL<pre) - 1)]; + *n = 0; + k = yak_pt_get(g->h, hash >> h->pre << YAK_COUNTER_BITS); + if (k == kh_end(g->h)) return 0; + *n = kh_key(g->h, k) & YAK_MAX_COUNT; + return &g->al[kh_val(g->h, k)]; +} + const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash) { khint_t k; @@ -703,7 +720,7 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for break;\ }\ } else if(p->us_in) {\ - ma_utg_t *u; s->uq = p->us_in->h;\ + 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) {\ @@ -1052,7 +1069,7 @@ void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k ///HAF_COUNT_EXACT ---> no bf; HAF_COUNT_ALL ---> no minimizer h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, is_HPC, k, w, NULL, NULL, NULL, us, 0, NULL); ha_ct_shrink(h, min_freq, max_freq>YAK_MAX_COUNT-1?YAK_MAX_COUNT-1:max_freq, asm_opt->thread_num); - flt_tab = gen_hh(h, asm_opt->max_kmer_cnt); + flt_tab = gen_hh(h, YAK_MAX_COUNT); ha_ct_destroy(h); return (void*)flt_tab; } diff --git a/htab.h b/htab.h index b2f9351..259de57 100644 --- a/htab.h +++ b/htab.h @@ -62,6 +62,9 @@ typedef struct ha_pt_s ha_pt_t; struct ha_abuf_s; typedef struct ha_abuf_s ha_abuf_t; +struct ha_abufl_s; +typedef struct ha_abufl_s ha_abufl_t; + extern const unsigned char seq_nt4_table[256]; extern void *ha_flt_tab; extern ha_pt_t *ha_idx; @@ -78,6 +81,7 @@ ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov); void ha_pt_destroy(ha_pt_t *h); const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n); +const ha_idxposl_t *ha_ptl_get(const ha_pt_t *h, uint64_t hash, int *n); const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash); int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name); @@ -89,6 +93,9 @@ int query_ct_index(void* ct_idx, uint64_t hash); ha_abuf_t *ha_abuf_init(void); void ha_abuf_destroy(ha_abuf_t *ab); uint64_t ha_abuf_mem(const ha_abuf_t *ab); +ha_abufl_t *ha_abufl_init(void); +void ha_abufl_destroy(ha_abufl_t *ab); +uint64_t ha_abufl_mem(const ha_abufl_t *ab); double yak_cputime(void); void yak_reset_realtime(void); diff --git a/inter.cpp b/inter.cpp index 3b5ef63..ce6b2d5 100644 --- a/inter.cpp +++ b/inter.cpp @@ -2,18 +2,45 @@ #include #include #include +#include "kseq.h" // FASTA/Q parser +#include "kthread.h" #include "inter.h" #include "Overlaps.h" #include "CommandLines.h" #include "htab.h" +KSEQ_INIT(gzFile, gzread) -void uidx_build(ma_ug_t *ug, int is_HPC, int k, int w, int hap_n) +typedef struct { + int w, k, bw, max_gap, is_HPC, hap_n; +} mg_idxopt_t; + +typedef struct { // global data structure for kt_pipeline() + const void *ha_flt_tab; + const ha_pt_t *ha_idx; + const mg_idxopt_t *opt; + kseq_t *ks; + int64_t chunk_size; + uint64_t n_thread; + uint64_t total_base; + uint64_t total_pair; +} uldat_t; + +void init_mg_opt(mg_idxopt_t *opt, int is_HPC, int k, int w, int hap_n) +{ + opt->k = k; + opt->w = w; + opt->hap_n = hap_n; + opt->is_HPC = is_HPC; + opt->bw = 2000; + opt->max_gap = 5000; +} + +void uidx_build(ma_ug_t *ug, mg_idxopt_t *opt) { int flag = asm_opt.flag; asm_opt.flag |= HA_F_NO_HPC; - ug->u.h = hap_n; - ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), is_HPC, k, w, hap_n, hap_n*10); - ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), is_HPC, k, w, hap_n); + ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), opt->is_HPC, opt->k, opt->w, 1, opt->hap_n*10); + ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), opt->is_HPC, opt->k, opt->w, 1); asm_opt.flag = flag; } @@ -23,10 +50,110 @@ void uidx_destory() ha_pt_destroy(ha_idx); } +static void *worker_ul_pipeline(void *data, int step, void *in) // callback for kt_pipeline() +{ + /** + uldat_t *p = (uldat_t*)data; + ///uint64_t total_base = 0, total_pair = 0; + if (step == 0) { // step 1: read a block of sequences + int ret1, ret2; + uint64_t l1, l2; + stepdat_t *s; + CALLOC(s, 1); + s->idx = p->idx; s->id = p->total_pair; s->t_ch = p->t_ch; + while (((ret1 = kseq_read(p->ks1)) >= 0)&&((ret2 = kseq_read(p->ks2)) >= 0)) + { + if (p->ks1->seq.l < p->idx->k || p->ks2->seq.l < p->idx->k) continue; + if (s->n == s->m) { + s->m = s->m < 16? 16 : s->m + (s->n>>1); + REALLOC(s->len, s->m); + REALLOC(s->seq, s->m); + } + + l1 = p->ks1->seq.l; l2 = p->ks2->seq.l; + MALLOC(s->seq[s->n], l1+l2); + s->sum_len += l1+l2; + memcpy(s->seq[s->n], p->ks1->seq.s, l1); + memcpy(s->seq[s->n]+l1, p->ks2->seq.s, l2); + s->len[s->n++] = (uint64_t)(l1<<32)|(uint64_t)l2; + + if (s->sum_len >= p->chunk_size) break; + } + p->total_pair += s->n; + if (s->sum_len == 0) free(s); + else return s; + } + else if (step == 1) { // step 2: alignment + stepdat_t *s = (stepdat_t*)in; + CALLOC(s->pos_buf, p->n_thread); + CALLOC(s->pos, s->n); + int i; + kt_for(p->n_thread, worker_for_alignment, s, s->n); + for (i = 0; i < s->n; ++i) { + free(s->seq[i]); + p->total_base += (s->len[i]>>32) + (uint32_t)s->len[i]; + } + + free(s->seq); free(s->len); + for (i = 0; i < (int)p->n_thread; ++i) { + free(s->pos_buf[i].a.a); + } + free(s->pos_buf); + return s; + } + else if (step == 2) { // step 3: dump + stepdat_t *s = (stepdat_t*)in; + int i; + for (i = 0; i < s->n; ++i) { + // if(s->pos[i].a == NULL) continue; + // kv_push(pe_hit_hap, p->hits, s->pos[i]); + if(s->pos[i].s == (uint64_t)-1) continue; + kv_push(pe_hit, p->hits.a, s->pos[i]); + } + free(s->pos); + free(s); + } + **/ + return 0; +} + +int alignment_ul_pipeline(uldat_t* sl, const enzyme *fn) +{ + double index_time = yak_realtime(); + int i; + for (i = 0; i < fn->n; i++){ + gzFile fp; + if ((fp = gzopen(fn->a[i], "r")) == 0) return 0; + sl->ks = kseq_init(fp); + kt_pipeline(3, worker_ul_pipeline, sl, 3); + kseq_destroy(sl->ks); + gzclose(fp); + } + fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time); + return 1; +} + +int ul_align(mg_idxopt_t *opt, const enzyme *fn, void *ha_flt_tab, ha_pt_t *ha_idx) +{ + uldat_t sl; memset(&sl, 0, sizeof(sl)); + sl.ha_flt_tab = ha_flt_tab; + sl.ha_idx = ha_idx; + sl.opt = opt; + sl.chunk_size = 20000000; + sl.n_thread = asm_opt.thread_num; + alignment_ul_pipeline(&sl, fn); + return 1; +} void ul_resolve(ma_ug_t *ug, int hap_n) { - uidx_build(ug, 1, 63, 63, hap_n); - + mg_idxopt_t opt; + init_mg_opt(&opt, 0, 19, 10, hap_n); + uidx_build(ug, &opt); + ul_align(&opt, asm_opt.ar, ha_flt_tab, ha_idx); + + + + uidx_destory(); } \ No newline at end of file diff --git a/sketch.cpp b/sketch.cpp index 1039eea..df6f60e 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -506,10 +506,10 @@ void sf##_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is 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(is_unique){\ - filtered = (cnt < is_unique);\ - cnt = (cnt == is_unique? 0:cnt);\ + filtered = (cnt >= 1<<28);\ + if(is_unique && (!filtered)) {\ + 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**/\