xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-10-11 06:20:56 +08:00
Code Issues Packages Projects Releases Wiki Activity

tmp macro

Browse Source
This commit is contained in:
chhylp123
2021-08-30 09:30:48 -04:00
parent e230b3a55f
commit 17b5cb21d3
9 changed files with 157 additions and 252 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.cpp
+3
View File
@@ -158,8 +158,11 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->thread_num = 1; asm_opt->thread_num = 1;
asm_opt->k_mer_length = 51; asm_opt->k_mer_length = 51;
asm_opt->hic_mer_length = 31; asm_opt->hic_mer_length = 31;
asm_opt->ul_mer_length = 19;
asm_opt->mz_win = 51; asm_opt->mz_win = 51;
asm_opt->ul_mz_win = 19;
asm_opt->mz_rewin = 1000; asm_opt->mz_rewin = 1000;
asm_opt->ul_mz_rewin = 360;
asm_opt->mz_sample_dist = 500; asm_opt->mz_sample_dist = 500;
asm_opt->bf_shift = 37; asm_opt->bf_shift = 37;
asm_opt->max_kmer_cnt = 2000; asm_opt->max_kmer_cnt = 2000;
CommandLines.h
+3
View File
@@ -44,9 +44,12 @@ typedef struct {
int thread_num; int thread_num;
int k_mer_length; int k_mer_length;
int hic_mer_length; int hic_mer_length;
int ul_mer_length;
int bub_mer_length; int bub_mer_length;
int mz_win; int mz_win;
int ul_mz_win;
int mz_rewin; int mz_rewin;
int ul_mz_rewin;
int mz_sample_dist; int mz_sample_dist;
int bf_shift; int bf_shift;
int max_kmer_cnt; int max_kmer_cnt;
Makefile
+2 -1
View File
@@ -7,7 +7,7 @@ INCLUDES=
OBJS= CommandLines.o Process_Read.o Assembly.o Hash_Table.o \ 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 \ 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 \ 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 EXE= hifiasm
LIBS= -lz -lpthread -lm LIBS= -lz -lpthread -lm
@@ -76,3 +76,4 @@ hic.o: hic.h
rcut.o: rcut.h rcut.o: rcut.h
horder.o: horder.h horder.o: horder.h
tovlp.o: tovlp.h tovlp.o: tovlp.h
inter.o: inter.h
anchor.cpp
+1 -1
View File
@@ -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 rlen = Get_READ_LENGTH(R_INF, rid); // read length
// get the list of anchors // 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 // minimizer of queried read
if (ab->mz.m > ab->old_mz_m) { if (ab->mz.m > ab->old_mz_m) {
ab->old_mz_m = ab->mz.m; ab->old_mz_m = ab->mz.m;
htab.cpp
+103 -120
View File
@@ -38,6 +38,75 @@ void *ha_flt_tab_hp;
ha_pt_t *ha_idx_hp; ha_pt_t *ha_idx_hp;
void *ha_ct_table; 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<<p->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<<p) - 1);\
ch_buf_t *b = &buf[pre];\
if (b->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 * * Yak specific parameters *
***************************/ ***************************/
@@ -303,6 +372,7 @@ typedef struct {
yak_pt_t *h; yak_pt_t *h;
uint64_t n; uint64_t n;
ha_idxpos_t *a; ha_idxpos_t *a;
ha_idxposl_t *al;
} ha_pt1_t; } ha_pt1_t;
struct ha_pt_s { struct ha_pt_s {
@@ -314,6 +384,7 @@ struct ha_pt_s {
typedef struct { typedef struct {
const ha_ct_t *ct; const ha_ct_t *ct;
ha_pt_t *pt; ha_pt_t *pt;
int is_l;
} pt_gen_aux_t; } 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); yak_ct_destroy(g);
a->ct->h[i].h = 0; 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; pt_gen_aux_t a;
int i; 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(); pt->h[i].h = yak_pt_init();
yak_pt_resize(pt->h[i].h, kh_size(ct->h[i].h)); 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<<pt->pre); kt_for(n_thread, worker_pt_gen, &a, 1<<pt->pre);
free(ct->h); free(ct); free(ct->h); free(ct);
return pt; return pt;
@@ -534,7 +606,8 @@ typedef struct {
int n, m; int n, m;
uint64_t n_ins; uint64_t n_ins;
uint64_t *a; uint64_t *a;
ha_mz1_t *b; ha_mz1_t *b_mz;
ha_mzl_t *b_mzl;
} ch_buf_t; } ch_buf_t;
///p = 12 ///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; 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<<p) - 1);
ch_buf_t *b = &buf[pre];
if (b->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 ///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 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; pl_data_t *p;
uint64_t n_seq0; ///the start index of current buffer block at R_INF uint64_t n_seq0; ///the start index of current buffer block at R_INF
///sum_len = total bases, nk = number of k-mers ///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; int *len;
char **seq; char **seq;
ha_mz1_v *mz_buf; ha_mz1_v *mz_buf;
ha_mz1_v *mz; ha_mz1_v *mz;
ha_mzl_v *mzl_buf;
ha_mzl_v *mzl;
ch_buf_t *buf; ch_buf_t *buf;
st_mt_t *mt; st_mt_t *mt;
} st_data_t; } 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]; ch_buf_t *b = &s->buf[i];
if (s->p->pt) if (s->p->pt)
{ {
if(s->p->flag&HAF_COUNT_REFINE) 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);
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);
}
} }
else///for 0-th count, go into here 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]; ha_mz1_v *b = &s->mz_buf[tid];
s->mz_buf[tid].n = 0; 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, 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; s->mz[i].n = s->mz[i].m = b->n;
MALLOC(s->mz[i].a, b->n); MALLOC(s->mz[i].a, b->n);
memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t)); 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() static void *worker_count(void *data, int step, void *in) // callback for kt_pipeline()
{ {
pl_data_t *p = (pl_data_t*)data; 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; break;
} }
} else if(p->us_in) { } else if(p->us_in) {
ma_utg_t *u; ma_utg_t *u; s->uq = 1;
while (p->n_seq < p->us_in->n) { while (p->n_seq < p->us_in->n) {
u = &(p->us_in->a[p->n_seq]); u = &(p->us_in->a[p->n_seq]);
if (s->n_seq == s->m_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 } else if (step == 1) { // step 2: extract k-mers
///s is the block of reads ///s is the block of reads
st_data_t *s = (st_data_t*)in; st_data_t *s = (st_data_t*)in;
///for 0-th counting, n_pre = 4096 if(p->us_in) mzl_init_kuf(p, s);
int i, n_pre = 1<<p->opt->pre, m; else mz_init_kuf(p, s);
// 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);
}
// fill the buffer // fill the buffer
///for 0-th counting, p->opt->w == 1 ///for 0-th counting, p->opt->w == 1
if (p->opt->w == 1) { // enumerate all k-mers if (p->opt->w == 1) { // enumerate all k-mers
///scan all reads ///scan all reads
int i;
for (i = 0; i < s->n_seq; ++i) { for (i = 0; i < s->n_seq; ++i) {
if (p->opt->is_HPC) if (p->opt->is_HPC)
count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]); 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]); if (!p->is_store) free(s->seq[i]);
} }
} else { // minimizers only } else { // minimizers only
uint32_t j; if(p->us_in) mzl_mselect(p, s);
// compute minimizers else mz_mselect(p, s);
// 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);
} }
///just clean seq ///just clean seq
free(s->seq); free(s->len); 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; return s;
} else if (step == 2) { // step 3: insert k-mers to hash table } else if (step == 2) { // step 3: insert k-mers to hash table
st_data_t *s = (st_data_t*)in; st_data_t *s = (st_data_t*)in;
int i, n = 1<<p->opt->pre;
uint64_t n_ins = 0;
///for 0-th counting, p->pt = NULL ///for 0-th counting, p->pt = NULL
kt_for(p->opt->n_thread, worker_for_insert, s, n); kt_for(p->opt->n_thread, worker_for_insert, s, 1<<p->opt->pre);
///n_ins is number of distinct k-mers if(p->us_in) mzl_destory_kuf(p, s, 1<<p->opt->pre);
for (i = 0; i < n; ++i) { else mz_destory_kuf(p, s, 1<<p->opt->pre);
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);
} }
return 0; 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); malloc_All_reads(rs);
} }
yak_copt_init(&opt); 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 ///always 0
opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC); opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC);
///for ft-counting, shoud be 1 ///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 ft-counting, shoud be 37
///for ha_pt_gen, shoud be 0 ///for ha_pt_gen, shoud be 0
opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift; 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.adaLen = (keep_adapter? asm_opt->adapterLen : 0);
opt.min_rcnt = (low_freq?*low_freq:-1); opt.min_rcnt = (low_freq?*low_freq:-1);
///asm_opt->num_reads is the number of fastq files ///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); 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; 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) void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int hap_n)
{ {
yak_ft_t *flt_tab; yak_ft_t *flt_tab;
int64_t cnt[YAK_N_COUNTS];
int cutoff = hap_n + 1;
ha_ct_t *h; 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); h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, NULL, NULL, NULL, us, 0, NULL);
ha_ct_shrink(h, 1, YAK_MAX_COUNT-1, asm_opt->thread_num);
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);
flt_tab = gen_hh(h, asm_opt->max_kmer_cnt); flt_tab = gen_hh(h, asm_opt->max_kmer_cnt);
ha_ct_destroy(h); 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; 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) 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_ct_t *ct;
ha_pt_t *pt; ha_pt_t *pt;
///HAF_COUNT_EXACT: no bf ///HAF_COUNT_EXACT: no bf
ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, NULL, flt_tab, NULL, us, 0, NULL); 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__, fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
yak_realtime(), yak_cpu_usage(), (long)ct->tot); 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); pt = ha_pt_gen(ct, asm_opt->thread_num, 1);
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);
ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, pt, flt_tab, NULL, us, 0, NULL); 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); //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_cpu_usage(), (long)pt->tot_pos); 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)) if(!(asm_opt->flag & HA_F_FAST))
{ {
fprintf(stderr, "[M::%s::] counting in normal mode\n", __func__); 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); ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, pt, flt_tab, rs, NULL, 1, NULL);
assert((uint64_t)tot_cnt == pt->tot_pos); assert((uint64_t)tot_cnt == pt->tot_pos);
} }
htab.h
+8 -8
View File
@@ -25,16 +25,16 @@ typedef struct { uint32_t n, m; ha_mz1_t *a; } ha_mz1_v;
typedef struct { typedef struct {
uint64_t x; ///x is the hash key 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 uint64_t rid:31, rev:1;
///span is the length of this k-mer. For non-HPC k-mer, span may not be equal to k uint32_t pos;
uint64_t rid:30, pos:34; uint8_t span;
uint16_t rev:1, span:15;
} ha_mzl_t; } ha_mzl_t;
typedef struct { typedef struct {
uint64_t rid:30, pos:34; uint64_t rid:31, rev:1;
uint16_t rev:1, span:15; uint32_t pos;
} ha_mzl_idxpos_t; uint8_t span;
} ha_idxposl_t;
typedef struct { uint32_t n, m; ha_mzl_t *a; } ha_mzl_v; 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_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 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); 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); void print_hist_lines(int n_cnt, int start_cnt, const int64_t *cnt);
inter.cpp
+25
View File
@@ -0,0 +1,25 @@
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <zlib.h>
#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);
}
inter.h
+5
View File
@@ -0,0 +1,5 @@
#ifndef __INTER__
#define __INTER__
#endif
sketch.cpp
+7 -122
View File
@@ -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 is_hpc homopolymer-compressed or not
* @param p minimizers * @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 { ///in default, w = 51, k = 51, is_hpc = 1
/** /**
uint64_t x; 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]); y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
cnt = hf? ha_ft_cnt(hf, y) : 0; cnt = hf? ha_ft_cnt(hf, y) : 0;
filtered = (cnt >= 1<<28); 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 (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 (!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]++; 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); 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); // 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); if (sample_dist > w) select_mz_h(p, mt, len, sample_dist, ws, k, tl);
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); 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<<k) - 1, kmer[4] = {0,0,0,0};
int i, j, l, buf_pos, min_pos, kmer_span = 0;
ha_mz1_t buf[256], min = dummy;
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);
}
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 for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts
p->a[i].rid = rid; p->a[i].rid = rid;
} }