reverted to the older k_mer_hit definition

This commit is contained in:
Heng Li
2020-03-27 10:40:31 -04:00
parent 10ed8b36ec
commit 7271c106e4
4 changed files with 40 additions and 48 deletions
+18 -20
View File
@@ -347,7 +347,7 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp)
if(j != -1) if(j != -1)
{ {
distance_self_pos = a[current_j].self_offset - a[j].self_offset; 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; distance_gap = distance_pos > distance_self_pos? distance_pos - distance_self_pos : distance_self_pos - distance_pos;
indels += distance_gap; 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 // fill the score and backtrack arrays
for (i = 0; i < a_n; ++i) for (i = 0; i < a_n; ++i)
{ {
pos = ha_hit_get_offset(&a[i]); pos = a[i].offset;
self_pos = a[i].self_offset; self_pos = a[i].self_offset;
max_j = -1; max_j = -1;
max_score = min_score; 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 ///may have a pre-cut condition for j
for (j = i - 1; j >= 0; --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; distance_self_pos = self_pos - a[j].self_offset;
///a has been sorted by a[].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 ///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_score = dp->score[i];
max_i = i; max_i = i;
mini_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen, 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) else if(dp->score[i] == max_score)
{ {
tmp_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen, 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) 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 ///not a has been sorted by offset, that means has been sorted by query offset
i = max_i; i = max_i;
result->x_pos_e = a[i].self_offset; 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->shared_seed = max_score;
result->overlapLen = mini_xLen; result->overlapLen = mini_xLen;
distance_self_pos = result->x_pos_e - a[i].self_offset; 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; long long pre_distance_gap = distance_pos - distance_self_pos;
///record first site ///record first site
///the length of f_cigar should be at least 1 ///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) while (i >= 0)
{ {
distance_self_pos = result->x_pos_e - a[i].self_offset; 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; distance_gap = distance_pos - distance_self_pos;
if(distance_gap != pre_distance_gap) 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++; chainLen++;
result->x_pos_s = a[i].self_offset; 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]; 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) while (i >= 0)
{ {
distance_self_pos = result->x_pos_e - a[i].self_offset; 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; distance_gap = distance_pos - distance_self_pos;
if(distance_gap == pre_distance_gap) 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++; chainLen++;
result->x_pos_s = a[i].self_offset; 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]; i = dp->pre[i];
} }
} }
@@ -604,8 +604,8 @@ void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_r
i = 0; i = 0;
while (i < candidates->length) 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;
///reference read ///reference read
tmp_region.x_id = readID; tmp_region.x_id = readID;
@@ -623,9 +623,9 @@ void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_r
while (i < candidates->length 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; sub_region_end = i;
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++) for (; j < candidates->length; j++)
{ {
if ( 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 end_pos == (uint64_t)candidates->list[j].self_offset
&& &&
strand == ha_hit_get_rev(&candidates->list[j]) strand == candidates->list[j].strand
) )
{ {
break; break;
@@ -794,7 +794,6 @@ void init_Candidates_list(Candidates_list* l)
l->length = 0; l->length = 0;
l->size = 0; l->size = 0;
l->list = NULL; l->list = NULL;
l->tmp = NULL;
init_Chain_Data(&(l->chainDP)); init_Chain_Data(&(l->chainDP));
} }
@@ -807,7 +806,6 @@ void clear_Candidates_list(Candidates_list* l)
void destory_Candidates_list(Candidates_list* l) void destory_Candidates_list(Candidates_list* l)
{ {
free(l->list); free(l->list);
free(l->tmp);
destory_Chain_Data(&(l->chainDP)); destory_Chain_Data(&(l->chainDP));
} }
+2 -19
View File
@@ -72,7 +72,6 @@ typedef struct
CIGAR cigar; CIGAR cigar;
} window_list; } window_list;
typedef struct typedef struct
{ {
window_list* buffer; window_list* buffer;
@@ -128,25 +127,10 @@ typedef struct
typedef struct typedef struct
{ {
uint64_t opos; uint32_t readID:31, strand:1;
uint32_t self_offset; // offset on the target read uint32_t offset, self_offset;
} k_mer_hit; } 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 typedef struct
{ {
long long* score; long long* score;
@@ -160,7 +144,6 @@ typedef struct
typedef struct typedef struct
{ {
k_mer_hit* list; k_mer_hit* list;
k_mer_hit* tmp;
long long length; long long length;
long long size; long long size;
Chain_Data chainDP; Chain_Data chainDP;
+12 -8
View File
@@ -278,7 +278,7 @@ KRADIX_SORT_INIT(ha64, uint64_t, generic_key, 8)
typedef struct { typedef struct {
yak_pt_t *h; yak_pt_t *h;
uint64_t n; uint64_t n;
uint64_t *a; ha_idxpos_t *a;
} ha_pt1_t; } ha_pt1_t;
struct ha_pt_s { 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) { for (j = 0; j < n; ++j) {
uint64_t x = a[j].x >> h->pre; uint64_t x = a[j].x >> h->pre;
khint_t k; khint_t k;
int n;
ha_idxpos_t *p;
if ((a[j].x&mask) != (a[0].x&mask)) continue; if ((a[j].x&mask) != (a[0].x&mask)) continue;
k = yak_pt_get(g->h, x<<YAK_COUNTER_BITS); k = yak_pt_get(g->h, x<<YAK_COUNTER_BITS);
if (k == kh_end(g->h)) continue; if (k == kh_end(g->h)) continue;
assert((kh_key(g->h, k)&YAK_MAX_COUNT) < YAK_MAX_COUNT); n = kh_key(g->h, k) & YAK_MAX_COUNT;
g->a[kh_val(g->h, k) + (kh_key(g->h, k)&YAK_MAX_COUNT)] assert(n < 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; 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); ++kh_key(g->h, k);
++n_ins; ++n_ins;
} }
return n_ins; return n_ins;
} }
/*
static void worker_pt_sort(void *data, long i, int tid) static void worker_pt_sort(void *data, long i, int tid)
{ {
ha_pt_t *h = (ha_pt_t*)data; 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<<h->pre); kt_for(n_thread, worker_pt_sort, h, 1<<h->pre);
} }
*/
void ha_pt_destroy(ha_pt_t *h) void ha_pt_destroy(ha_pt_t *h)
{ {
int i; int i;
@@ -383,7 +387,7 @@ void ha_pt_destroy(ha_pt_t *h)
free(h->h); free(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; khint_t k;
const ha_pt1_t *g = &h->h[hash & ((1ULL<<h->pre) - 1)]; const ha_pt1_t *g = &h->h[hash & ((1ULL<<h->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); pt = ha_pt_gen(ct, asm_opt->thread_num);
ha_count(asm_opt, HAF_COUNT_EXACT|HAF_RS_READ, pt, flt_tab, rs); ha_count(asm_opt, HAF_COUNT_EXACT|HAF_RS_READ, pt, flt_tab, rs);
assert((uint64_t)tot_cnt == pt->tot_pos); 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__, fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__,
yak_realtime(), yak_cputime() / yak_realtime(), (long)pt->tot_pos); yak_realtime(), yak_cputime() / yak_realtime(), (long)pt->tot_pos);
return pt; return pt;
+8 -1
View File
@@ -10,11 +10,18 @@ typedef struct {
uint64_t rid:28, pos:27, rev:1, span:8; uint64_t rid:28, pos:27, rev:1, span:8;
} ha_mz1_t; } 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; typedef struct { uint32_t n, m; ha_mz1_t *a; } ha_mz1_v;
struct ha_pt_s; struct ha_pt_s;
typedef struct ha_pt_s ha_pt_t; 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]; extern const unsigned char seq_nt4_table[256];
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs); 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); 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); 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); double yak_cputime(void);
void yak_reset_realtime(void); void yak_reset_realtime(void);