xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-09-30 04:38:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

removed unused ha_sketch

Browse Source
This commit is contained in:
Heng Li
2020-10-29 11:28:43 -04:00
parent eddb8173da
commit ef8c047559
5 changed files with 43 additions and 141 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Overlaps.cpp
+2 -3
View File
@@ -9974,8 +9974,7 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source
void debug_info_of_specfic_read(const char* name, ma_hit_t_alloc* sources,
ma_hit_t_alloc* reverse_sources, int id, const char* command)
void debug_info_of_specfic_read(const char* name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int id, const char* command)
{
long long i, j, Len;
uint32_t tn;
@@ -9999,7 +9998,7 @@ ma_hit_t_alloc* reverse_sources, int id, const char* command)
fprintf(stderr, "\n\n\nafter %s\n", command);
fprintf(stderr, "****************ma_hit_t (%lld)ref_read: %.*s, len: %lu****************\n",
i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), Get_READ_LENGTH(R_INF, i));
i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), (unsigned long)Get_READ_LENGTH(R_INF, i));
fprintf(stderr, "sources Len: %d, is_fully_corrected: %d\n",
anchor.cpp
+5 -4
View File
@@ -57,8 +57,8 @@ int ha_ov_type(const overlap_region *r, uint32_t len)
else return r->x_pos_s == 0? 0 : 1;
}
void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, 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)
void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, 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)
{
uint32_t i, rlen;
uint64_t k, l;
@@ -205,7 +205,8 @@ void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf)
void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, 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, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct)
int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar,
kvec_t_u64_warp* dbg_ct)
{
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
@@ -302,4 +303,4 @@ int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp*
void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list)
{
ks_introsort_or_xs(overlap_list->length, overlap_list->list);
}
}
htab.cpp
+15 -5
View File
@@ -557,7 +557,7 @@ 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;
///s->p->opt->w = 51, s->p->opt->k
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);
ha_sketch_query(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, 0, 0);
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));
@@ -803,7 +803,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi
* High count filter table *
***************************/
KHASHL_SET_INIT(static klib_unused, yak_ft_t, yak_ft, uint64_t, kh_hash_dummy, kh_eq_generic)
KHASHL_MAP_INIT(static klib_unused, yak_ft_t, yak_ft, uint64_t, int16_t, kh_hash_dummy, kh_eq_generic)
static yak_ft_t *gen_hh(const ha_ct_t *h)
{
@@ -813,12 +813,14 @@ static yak_ft_t *gen_hh(const ha_ct_t *h)
yak_ft_resize(hh, h->tot * 2);
for (i = 0; i < 1<<h->pre; ++i) {
yak_ct_t *ht = h->h[i].h;
khint_t k;
khint_t k, l;
for (k = 0; k < kh_end(ht); ++k) {
if (kh_exist(ht, k)) {
uint64_t y = kh_key(ht, k) >> h->pre << YAK_COUNTER_BITS | i;
int absent;
yak_ft_put(hh, y, &absent);
l = yak_ft_put(hh, y, &absent);
if (absent)
kh_val(hh, l) = kh_key(ht, k)&YAK_MAX_COUNT;
}
}
}
@@ -833,6 +835,14 @@ int ha_ft_isflt(const void *hh, uint64_t y)
return k == kh_end(h)? 0 : 1;
}
int ha_ft_cnt(const void *hh, uint64_t y)
{
yak_ft_t *h = (yak_ft_t*)hh;
khint_t k;
k = yak_ft_get(h, y);
return k == kh_end(h)? 0 : kh_val(h, k);
}
void ha_ft_destroy(void *h)
{
if (h) yak_ft_destroy((yak_ft_t*)h);
@@ -1092,7 +1102,7 @@ int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_op
ha_pt_t *ha_idx = NULL;
char mode = 0;
int f_flag, absent, i;
int f_flag = 0, absent, i;
double index_time, index_s_time, pos_time, pos_s_time;
htab.h
+1 -1
View File
@@ -33,7 +33,7 @@ extern void *ha_ct_table;
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode);
int ha_ft_isflt(const void *hh, uint64_t y);
int ha_ft_cnt(const void *hh, uint64_t y);
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, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov);
sketch.cpp
+20 -128
View File
@@ -36,12 +36,13 @@ static inline int tq_shift(tiny_queue_t *q)
* @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)
void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct)
{ ///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, 0, 0, 0 };
uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0};
int i, j, l, buf_pos, min_pos, kmer_span = 0;
@@ -49,6 +50,13 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
tiny_queue_t tq;
assert(len > 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);
}
///sizeof(ha_mz1_t) = 16
memset(buf, 0xff, w * 16);
memset(&tq, 0, sizeof(tiny_queue_t));
@@ -70,12 +78,16 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
}
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
///so for non-HPC k-mer, kmer_span should be k in any case?
///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
@@ -85,13 +97,17 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
++l;
if (l >= k && kmer_span < 256) {
uint64_t y;
int filtered = 0;
y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
if (hf == 0 || ha_ft_isflt(hf, y) == 0)
if (hf != 0) filtered = (ha_ft_cnt(hf, y) > 0);
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 == 0)
info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;
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)
@@ -135,127 +151,3 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
if (min.x != UINT64_MAX)
kv_push(ha_mz1_t, *p, min);
}
void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf,
kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct)
{ ///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;
if(dbg_ct != NULL) dbg_ct->a.n = 0;
static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 };
uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0}, filtered;
int i, j, l, buf_pos, min_pos, kmer_span = 0;
ha_mz1_t buf[256], min = dummy;
tiny_queue_t tq;
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);
}
assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
///sizeof(ha_mz1_t) = 16
memset(buf, 0xff, w * 16);
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;
y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
filtered = 0;
if(hf != 0) filtered = ha_ft_isflt(hf, y);
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 (hf == 0 || ha_ft_isflt(hf, y) == 0)
if(filtered == 0)
info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;
if(k_flag != NULL) k_flag->a.a[i]++;
if(k_flag != NULL && filtered == 1) 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 (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
for (j = 0; j < buf_pos; ++j)
if (min.x == buf[j].x && 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 (min.x >= buf[j].x) 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 (min.x >= buf[j].x) 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 (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
for (j = 0; j <= buf_pos; ++j)
if (min.x == buf[j].x && 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);
}