diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 09667b0..d387da7 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -347,7 +347,7 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp) if(j != -1) { distance_self_pos = a[current_j].self_offset - a[j].self_offset; - distance_pos = ha_hit_get_offset(&a[current_j]) - ha_hit_get_offset(&a[j]); + distance_pos = a[current_j].offset - a[j].offset; distance_gap = distance_pos > distance_self_pos? distance_pos - distance_self_pos : distance_self_pos - distance_pos; indels += distance_gap; @@ -421,7 +421,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul // fill the score and backtrack arrays for (i = 0; i < a_n; ++i) { - pos = ha_hit_get_offset(&a[i]); + pos = a[i].offset; self_pos = a[i].self_offset; max_j = -1; max_score = min_score; @@ -433,7 +433,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul ///may have a pre-cut condition for j for (j = i - 1; j >= 0; --j) { - distance_pos = pos - ha_hit_get_offset(&a[j]); + distance_pos = pos - a[j].offset; distance_self_pos = self_pos - a[j].self_offset; ///a has been sorted by a[].offset ///note for a, we do not have any two elements that have both equal offsets and self_offsets @@ -502,12 +502,12 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul max_score = dp->score[i]; max_i = i; mini_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen, - ha_hit_get_offset(&a[i]), ha_hit_get_offset(&a[i]), y_readLen); + a[i].offset, a[i].offset, y_readLen); } else if(dp->score[i] == max_score) { tmp_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen, - ha_hit_get_offset(&a[i]), ha_hit_get_offset(&a[i]), y_readLen); + a[i].offset, a[i].offset, y_readLen); if(tmp_xLen < mini_xLen) { @@ -524,12 +524,12 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul ///not a has been sorted by offset, that means has been sorted by query offset i = max_i; result->x_pos_e = a[i].self_offset; - result->y_pos_e = ha_hit_get_offset(&a[i]); + result->y_pos_e = a[i].offset; result->shared_seed = max_score; result->overlapLen = mini_xLen; distance_self_pos = result->x_pos_e - a[i].self_offset; - distance_pos = result->y_pos_e - ha_hit_get_offset(&a[i]); + distance_pos = result->y_pos_e - a[i].offset; long long pre_distance_gap = distance_pos - distance_self_pos; ///record first site ///the length of f_cigar should be at least 1 @@ -541,7 +541,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul while (i >= 0) { distance_self_pos = result->x_pos_e - a[i].self_offset; - distance_pos = result->y_pos_e - ha_hit_get_offset(&a[i]); + distance_pos = result->y_pos_e - a[i].offset; distance_gap = distance_pos - distance_self_pos; if(distance_gap != pre_distance_gap) { @@ -552,7 +552,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul chainLen++; result->x_pos_s = a[i].self_offset; - result->y_pos_s = ha_hit_get_offset(&a[i]); + result->y_pos_s = a[i].offset; i = dp->pre[i]; } } @@ -562,7 +562,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul while (i >= 0) { distance_self_pos = result->x_pos_e - a[i].self_offset; - distance_pos = result->y_pos_e - ha_hit_get_offset(&a[i]); + distance_pos = result->y_pos_e - a[i].offset; distance_gap = distance_pos - distance_self_pos; if(distance_gap == pre_distance_gap) { @@ -577,7 +577,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul chainLen++; result->x_pos_s = a[i].self_offset; - result->y_pos_s = ha_hit_get_offset(&a[i]); + result->y_pos_s = a[i].offset; i = dp->pre[i]; } } @@ -604,8 +604,8 @@ void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_r i = 0; while (i < candidates->length) { - current_ID = ha_hit_get_readID(&candidates->list[i]); - current_stand = ha_hit_get_rev(&candidates->list[i]); + current_ID = candidates->list[i].readID; + current_stand = candidates->list[i].strand; ///reference read tmp_region.x_id = readID; @@ -623,9 +623,9 @@ void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_r while (i < candidates->length && - current_ID == ha_hit_get_readID(&candidates->list[i]) + current_ID == candidates->list[i].readID && - current_stand == ha_hit_get_rev(&candidates->list[i])) + current_stand == candidates->list[i].strand) { sub_region_end = i; i++; @@ -733,13 +733,13 @@ void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n for (; j < candidates->length; j++) { if ( - n_list[i].offset == (uint64_t)ha_hit_get_offset(&candidates->list[j]) + n_list[i].offset == (uint64_t)candidates->list[j].offset && - n_list[i].readID == ha_hit_get_readID(&candidates->list[j]) + n_list[i].readID == candidates->list[j].readID && end_pos == (uint64_t)candidates->list[j].self_offset && - strand == ha_hit_get_rev(&candidates->list[j]) + strand == candidates->list[j].strand ) { break; @@ -794,7 +794,6 @@ void init_Candidates_list(Candidates_list* l) l->length = 0; l->size = 0; l->list = NULL; - l->tmp = NULL; init_Chain_Data(&(l->chainDP)); } @@ -807,7 +806,6 @@ void clear_Candidates_list(Candidates_list* l) void destory_Candidates_list(Candidates_list* l) { free(l->list); - free(l->tmp); destory_Chain_Data(&(l->chainDP)); } diff --git a/Hash_Table.h b/Hash_Table.h index e5c5fed..1dd5ccd 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -72,7 +72,6 @@ typedef struct CIGAR cigar; } window_list; - typedef struct { window_list* buffer; @@ -128,25 +127,10 @@ typedef struct typedef struct { - uint64_t opos; - uint32_t self_offset; // offset on the target read + uint32_t readID:31, strand:1; + uint32_t offset, self_offset; } k_mer_hit; -static inline uint32_t ha_hit_get_readID(const k_mer_hit *h) -{ - return h->opos >> 33; -} - -static inline uint32_t ha_hit_get_rev(const k_mer_hit *h) -{ - return h->opos >> 32 & 1; -} - -static inline uint32_t ha_hit_get_offset(const k_mer_hit *h) -{ - return (uint32_t)h->opos; -} - typedef struct { long long* score; @@ -160,7 +144,6 @@ typedef struct typedef struct { k_mer_hit* list; - k_mer_hit* tmp; long long length; long long size; Chain_Data chainDP; diff --git a/htab.cpp b/htab.cpp index 6b696b9..ba139fb 100644 --- a/htab.cpp +++ b/htab.cpp @@ -278,7 +278,7 @@ KRADIX_SORT_INIT(ha64, uint64_t, generic_key, 8) typedef struct { yak_pt_t *h; uint64_t n; - uint64_t *a; + ha_idxpos_t *a; } ha_pt1_t; struct ha_pt_s { @@ -340,18 +340,22 @@ int ha_pt_insert_list(ha_pt_t *h, int n, const ha_mz1_t *a) for (j = 0; j < n; ++j) { uint64_t x = a[j].x >> h->pre; khint_t k; + int n; + ha_idxpos_t *p; if ((a[j].x&mask) != (a[0].x&mask)) continue; k = yak_pt_get(g->h, x<h)) continue; - assert((kh_key(g->h, k)&YAK_MAX_COUNT) < YAK_MAX_COUNT); - g->a[kh_val(g->h, k) + (kh_key(g->h, k)&YAK_MAX_COUNT)] - = (uint64_t)a[j].rid<<36 | (uint64_t)a[j].rev<<35 | (uint64_t)a[j].pos<<8 | (uint64_t)a[j].span; + n = kh_key(g->h, k) & YAK_MAX_COUNT; + assert(n < YAK_MAX_COUNT); + p = &g->a[kh_val(g->h, k) + n]; + p->rid = a[j].rid, p->rev = a[j].rev, p->pos = a[j].pos, p->span = a[j].span; + //(uint64_t)a[j].rid<<36 | (uint64_t)a[j].rev<<35 | (uint64_t)a[j].pos<<8 | (uint64_t)a[j].span; ++kh_key(g->h, k); ++n_ins; } return n_ins; } - +/* static void worker_pt_sort(void *data, long i, int tid) { ha_pt_t *h = (ha_pt_t*)data; @@ -371,7 +375,7 @@ void ha_pt_sort(ha_pt_t *h, int n_thread) { kt_for(n_thread, worker_pt_sort, h, 1<pre); } - +*/ void ha_pt_destroy(ha_pt_t *h) { int i; @@ -383,7 +387,7 @@ void ha_pt_destroy(ha_pt_t *h) free(h->h); free(h); } -const uint64_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n) +const ha_idxpos_t *ha_pt_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)]; @@ -799,7 +803,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f pt = ha_pt_gen(ct, asm_opt->thread_num); ha_count(asm_opt, HAF_COUNT_EXACT|HAF_RS_READ, pt, flt_tab, rs); assert((uint64_t)tot_cnt == pt->tot_pos); - ha_pt_sort(pt, asm_opt->thread_num); + //ha_pt_sort(pt, asm_opt->thread_num); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__, yak_realtime(), yak_cputime() / yak_realtime(), (long)pt->tot_pos); return pt; diff --git a/htab.h b/htab.h index d7232cf..18567f8 100644 --- a/htab.h +++ b/htab.h @@ -10,11 +10,18 @@ typedef struct { uint64_t rid:28, pos:27, rev:1, span:8; } ha_mz1_t; +typedef struct { + uint64_t rid:28, pos:27, rev:1, span:8; +} ha_idxpos_t; + typedef struct { uint32_t n, m; ha_mz1_t *a; } ha_mz1_v; struct ha_pt_s; typedef struct ha_pt_s ha_pt_t; +struct ha_abuf_s; +typedef struct ha_abuf_s ha_abuf_t; + extern const unsigned char seq_nt4_table[256]; void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs); @@ -23,7 +30,7 @@ void ha_ft_destroy(void *h); ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs); void ha_pt_destroy(ha_pt_t *h); -const uint64_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n); +const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n); double yak_cputime(void); void yak_reset_realtime(void);