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

change ug extraction

Browse Source
This commit is contained in:
chhylp123
2021-12-12 20:59:48 -05:00
parent 9ec286e64e
commit 8a0ec3d575
5 changed files with 167 additions and 36 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Overlaps.cpp
+2 -1
View File
@@ -31117,7 +31117,7 @@ float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name,
long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr,
ma_sub_t **coverage_cut_ptr, int debug_g)
{
// if(asm_opt.ar) create_ul_info();
if(asm_opt.ar) create_ul_info();
@@ -31400,6 +31400,7 @@ long long bubble_dist, int read_graph, int write)
init_R_to_U(&ruIndex, n_read);
asg_t *sg = NULL;
ma_sub_t* coverage_cut = NULL;
init_aux_table();
///actually min_thres = asm_opt.max_short_tip + 1 there are asm_opt.max_short_tip reads
min_thres = asm_opt.max_short_tip + 1;
Process_Read.cpp
+14 -1
View File
@@ -881,7 +881,6 @@ void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str,
}
}
void retrieve_ul_t(UC_Read* r, all_ul_t *ref, uint64_t ID, uint8_t strand) {
ul_vec_t *p = &(ref->a[ID]);
uc_block_t *b = NULL;
@@ -931,4 +930,18 @@ void retrieve_ul_t(UC_Read* r, all_ul_t *ref, uint64_t ID, uint8_t strand) {
}
for (k = 0; k < (int64_t)p->N_site.n; k++) r->seq[r->length - p->N_site.a[k] - 1] = 'N';
}
}
void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand)
{
if(i_r) {
}
if (r->length + 4 > r->size)
{
r->size = r->length + 4;
r->seq = (char*)realloc(r->seq,sizeof(char)*(r->size));
}
}
htab.cpp
+77 -32
View File
@@ -8,6 +8,7 @@
#include "kseq.h"
#include "ksort.h"
#include "htab.h"
#include "Process_Read.h"
#define YAK_COUNTER_BITS 12
#define YAK_N_COUNTS (1<<YAK_COUNTER_BITS)
@@ -44,7 +45,7 @@ void *ha_ct_table;
typedef struct {
int32_t bf_shift, bf_n_hash;
int32_t k, w, is_HPC;
int32_t k, w, is_HPC, uq;
int32_t pre;
int32_t n_thread;
int64_t chunk_size;
@@ -562,7 +563,7 @@ KSEQ_INIT(gzFile, gzread)
typedef struct { // global data structure for kt_pipeline()
const yak_copt_t *opt;
const void *flt_tab;
int flag, create_new, is_store;
int flag, create_new, is_store, uq;
uint64_t n_mz, n_seq; ///number of total reads
kseq_t *ks;
UC_Read ucr;
@@ -697,6 +698,7 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for
CALLOC(s, 1);\
s->p = p;\
s->n_seq0 = p->n_seq;\
s->uq = p->opt->uq;\
if (p->rs_in && (p->flag & HAF_RS_READ)) {\
while (p->n_seq < p->rs_in->total_reads) {\
if ((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU) {\
@@ -721,7 +723,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 = 1;\
ma_utg_t *u;\
while (p->n_seq < p->us_in->n) {\
u = &(p->us_in->a[p->n_seq]);\
if (s->n_seq == s->m_seq) {\
@@ -730,7 +732,8 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for
REALLOC(s->seq, s->m_seq);\
}\
MALLOC(s->seq[s->n_seq], u->len);\
memcpy(s->seq[s->n_seq], u->s, u->len);\
if(u->s) memcpy(s->seq[s->n_seq], u->s, u->len);\
else memcpy(s->seq[s->n_seq], u->s, u->len);\
s->len[s->n_seq++] = u->len;\
++p->n_seq;\
s->sum_len += u->len;\
@@ -905,7 +908,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
pl_data_t pl;
gzFile fp = 0;
memset(&pl, 0, sizeof(pl_data_t));
pl.n_seq = *n_seq;
pl.n_seq = *n_seq;
if(ug_rs) {
pl.us_in = us;
} else if (read_rs) {
@@ -950,7 +953,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
return pl.ct;
}
ha_ct_t *ha_count(const hifiasm_opt_t *asm_o, int flag, int HPC, int k, int w, ha_pt_t *p0, const void *flt_tab, All_reads *rs, ma_utg_v *us, int keep_adapter, int *low_freq)
ha_ct_t *ha_count(const hifiasm_opt_t *asm_o, int flag, int HPC, int k, int w, ha_pt_t *p0, const void *flt_tab, All_reads *rs, ma_utg_v *us, int keep_adapter, int *low_freq, int unique_only)
{
int i;
int64_t n_seq = 0;
@@ -976,6 +979,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_o, int flag, int HPC, int k, int w, h
opt.n_thread = asm_o->thread_num;
opt.adaLen = (keep_adapter? asm_o->adapterLen : 0);
opt.min_rcnt = (low_freq?*low_freq:-1);
opt.uq = (unique_only?1:0);
///asm_opt->num_reads is the number of fastq files
for (i = n_bs = 0; i < (us?1:asm_o->num_reads); ++i){
h = yak_count(&opt, asm_o->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq);
@@ -1062,39 +1066,35 @@ void debug_ct_index(void* q_ct_idx, void* r_ct_idx)
/*************************
* High-level interfaces *
*************************/
void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff)
{
yak_ft_t *flt_tab;
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, !(asm_opt->flag&HA_F_NO_HPC), k, w, NULL, NULL, NULL, us, 0, NULL, 0);
// cutoff = (int)(asm_opt->hom_cov * asm_opt->high_factor);
if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num);
flt_tab = gen_hh(h, asm_opt->max_kmer_cnt);
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;
}
void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq)
{
yak_ft_t *flt_tab;
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, is_HPC, k, w, NULL, NULL, NULL, us, 0, NULL);
h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, is_HPC, k, w, NULL, NULL, NULL, us, 0, NULL, 0);
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, YAK_MAX_COUNT);
ha_ct_destroy(h);
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 is_HPC, int k, int w, int min_freq)
{
ha_ct_t *ct;
ha_pt_t *pt;
///HAF_COUNT_EXACT: no bf
ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, NULL, flt_tab, NULL, us, 0, NULL);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
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, min_freq, YAK_MAX_COUNT - 1, asm_opt->thread_num);
pt = ha_pt_gen(ct, asm_opt->thread_num, 1);
ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, pt, flt_tab, NULL, us, 0, NULL);
//ha_pt_sort(pt, asm_opt->thread_num);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__,
yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos);
return pt;
}
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode)
{
yak_ft_t *flt_tab;
@@ -1102,7 +1102,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i
int peak_hom, peak_het, cutoff = YAK_MAX_COUNT - 1, ex_flag = 0;
if(is_hp_mode) ex_flag = HAF_RS_READ|HAF_SKIP_READ;
ha_ct_t *h;
h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, NULL, rs, NULL, 1, NULL);
h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, NULL, rs, NULL, 1, NULL, 0);
if((asm_opt->flag & HA_F_VERBOSE_GFA))
{
write_ct_index((void*)h, asm_opt->output_file_name);
@@ -1130,16 +1130,61 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i
return (void*)flt_tab;
}
ha_pt_t *ha_pt_ul_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int k, int w, int cutoff)
{
ha_ct_t *ct;
ha_pt_t *pt;
///HAF_COUNT_EXACT: no bf
ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, !(asm_opt->flag&HA_F_NO_HPC), k, w, NULL, flt_tab, NULL, us, 0, NULL, 0);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
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
if (flt_tab == 0) {
if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
ha_ct_shrink(ct, 2, cutoff, asm_opt->thread_num);
} 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);
}
pt = ha_pt_gen(ct, asm_opt->thread_num, 1);
ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, !(asm_opt->flag&HA_F_NO_HPC), k, w, pt, flt_tab, NULL, us, 0, NULL, 0);
//ha_pt_sort(pt, asm_opt->thread_num);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__,
yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos);
return pt;
}
ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int is_HPC, int k, int w, int min_freq)
{
ha_ct_t *ct;
ha_pt_t *pt;
///HAF_COUNT_EXACT: no bf
ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, NULL, flt_tab, NULL, us, 0, NULL, 1);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
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, min_freq, YAK_MAX_COUNT - 1, asm_opt->thread_num);
pt = ha_pt_gen(ct, asm_opt->thread_num, 1);
ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, pt, flt_tab, NULL, us, 0, NULL, 1);
//ha_pt_sort(pt, asm_opt->thread_num);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__,
yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos);
return pt;
}
ha_pt_t *ha_pt_gen_dp(const hifiasm_opt_t *asm_opt, ha_ct_t *ct, int flag, int n_thread, const void *flt_tab, All_reads *rs, int peak_hom, int peak_het)
{
int low_freq = mz_low_b(peak_hom, peak_het);
ha_pt_t *pt = ha_pt_gen_count(ct, n_thread); ///key = cnt, val = 0
ha_count(asm_opt, HAF_COUNT_EXACT|HAF_COUNT_REFINE|flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, &low_freq);
ha_count(asm_opt, HAF_COUNT_EXACT|HAF_COUNT_REFINE|flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, &low_freq, 0);
uint64_t occ = ha_pt_shrink(pt, n_thread);
if(flag&HAF_RS_WRITE_LEN) flag -= HAF_RS_WRITE_LEN;
if(flag&HAF_RS_WRITE_SEQ) flag -= HAF_RS_WRITE_SEQ;
flag |= HAF_RS_READ; pt->tot_pos = 0;
ha_count(asm_opt, HAF_COUNT_EXACT|flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, NULL);
ha_count(asm_opt, HAF_COUNT_EXACT|flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, NULL, 0);
// fprintf(stderr, "[M::%s::] counted %lu distinct minimizer k-mers\n", __func__, pt->tot);
// fprintf(stderr, "[M::%s::] collected %lu minimizers\n\n\n", __func__, pt->tot_pos);
assert(occ == pt->tot_pos);
@@ -1163,7 +1208,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
}
if(is_hp_mode) extra_flag1 |= HAF_SKIP_READ, extra_flag2 |= HAF_SKIP_READ;
ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, flt_tab, rs, NULL, 1, NULL);
ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, flt_tab, rs, NULL, 1, NULL, 0);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
yak_realtime(), yak_cpu_usage(), (long)ct->tot);
ha_ct_hist(ct, cnt, asm_opt->thread_num);
@@ -1189,7 +1234,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
{
fprintf(stderr, "[M::%s::] counting in normal mode\n", __func__);
pt = ha_pt_gen(ct, asm_opt->thread_num, 0);
ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, NULL);
ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, NULL, 0);
assert((uint64_t)tot_cnt == pt->tot_pos);
}
else
htab.h
+2
View File
@@ -72,11 +72,13 @@ extern void *ha_flt_tab_hp;
extern ha_pt_t *ha_idx_hp;
extern void *ha_ct_table;
void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff);
void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq);
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode);
int32_t ha_ft_cnt(const void *hh, uint64_t y);
void ha_ft_destroy(void *h);
ha_pt_t *ha_pt_ul_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int k, int w, int cutoff);
ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int is_HPC, int k, int w, int min_freq);
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);
inter.cpp
+72 -2
View File
@@ -322,6 +322,13 @@ void init_mg_opt(mg_idxopt_t *opt, int is_HPC, int k, int w, int hap_n)
opt->pri_ratio = 0.8f;
}
void uidx_l_build(ma_ug_t *ug, mg_idxopt_t *opt, int cutoff)
{
ha_flt_tab = ha_ft_ul_gen(&asm_opt, &(ug->u), opt->k, opt->w, cutoff);
ha_idx = ha_pt_ul_gen(&asm_opt, ha_flt_tab, &(ug->u), opt->k, opt->w, cutoff);
fprintf(stderr, "[M::%s] Index has been built.\n", __func__);
}
void uidx_build(ma_ug_t *ug, mg_idxopt_t *opt)
{
int flag = asm_opt.flag;
@@ -2345,7 +2352,6 @@ int scall_ul_pipeline(uldat_t* sl, const enzyme *fn)
{
double index_time = yak_realtime();
int i;
init_aux_table();
init_all_ul_t(&UL_INF, &R_INF);
for (i = 0; i < fn->n; i++){
@@ -2885,10 +2891,74 @@ int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const asg_t *rg, const enz
return 1;
}
ma_ug_t *dedup_HiFis(ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, uint64_t n_read, int64_t readLen)
{
uint64_t i, k, qn, tn;
int32_t r; asg_arc_t t, *p = NULL;
asg_t *rg = asg_init();
rg->m_seq = n_read; MALLOC(rg->seq, rg->m_seq);
for (i = 0; i < n_read; ++i) rg->seq[i].len = Get_READ_LENGTH(R_INF, i);
for (i = 0; i < n_read; i++) {
if(rg->seq[i].del) continue;
for (k = 0; k < src[i].length; k++) {
if(!src[i].buffer[k].el) continue;
qn = Get_qn(src[i].buffer[k]); tn = Get_tn(src[i].buffer[k]);
if(rg->seq[qn].del || rg->seq[tn].del) continue;
r = ma_hit2arc(&(src[i].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
if (r == MA_HT_QCONT) {
rg->seq[qn].del = 1;
} else if(r == MA_HT_TCONT) {
rg->seq[tn].del = 1;
}
if(rg->seq[i].del) break;
}
}
for (i = 0; i < n_read; i++) {
if(rg->seq[i].del) continue;
for (k = 0; k < src[i].length; k++) {
if(!src[i].buffer[k].el) continue;
qn = Get_qn(src[i].buffer[k]); tn = Get_tn(src[i].buffer[k]);
if(rg->seq[qn].del || rg->seq[tn].del) continue;
r = ma_hit2arc(&(src[i].buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
if (r >= 0) {
p = asg_arc_pushp(rg);
*p = t;
}
}
}
asg_cleanup(rg); asg_symm(rg);
ma_ug_t *ug = NULL;
ug = ma_ug_gen(rg);
asg_destroy(rg); ma_ug_destroy(ug);
}
void ul_load(const ug_opt_t *uopt)
{
fprintf(stderr, "[M::%s::] ==> UL\n", __func__);
mg_idxopt_t opt;
init_mg_opt(&opt, 0, 19, 10, 4095);
ma_ug_t *ug = NULL;
int cutoff = asm_opt.hom_cov * asm_opt.high_factor;
init_aux_table();
init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff);
int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name) : 0);
if(exist == 0) uidx_l_build(ug, &opt, cutoff);
if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name);
ul_v_call(&opt, uopt, /**rg**/NULL, asm_opt.ar, ha_flt_tab, ha_idx, /**ug**/NULL);
}