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

debug ul

Browse Source
This commit is contained in:
chhylp123
2021-10-10 13:37:32 -04:00
parent 7dd4a848b9
commit cb320b88f3
7 changed files with 1212 additions and 325 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.h
+1 -1
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.16.1-r377"
#define HA_VERSION "0.16.2-r378"
#define VERBOSE 0
Overlaps.cpp
+152
View File
@@ -15,6 +15,7 @@
#include "Assembly.h"
#include "rcut.h"
#include "horder.h"
#include "inter.h"
uint32_t debug_purge_dup = 0;
@@ -13280,6 +13281,151 @@ long long gap_fuzz, bub_label_t* b_mask_t)
}
}
void print_debug_gfa(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name,
ma_hit_t_alloc* sources, R_to_U* ruIndex)
{
char* gfa_name = (char*)malloc(strlen(output_file_name)+50);
sprintf(gfa_name, "%s.after.clean_d_utg.gfa", output_file_name);
FILE* output_file = fopen(gfa_name, "w");
ma_ug_print(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file);
fclose(output_file);
free(gfa_name);
}
void output_ul_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name,
ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
long long tipsLen, float tip_drop_ratio, long long stops_threshold,
R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp,
long long gap_fuzz, bub_label_t* b_mask_t)
{
hic_clean(sg);
ug_opt_t opt; memset(&opt, 0, sizeof(opt));
// kvec_pe_hit *rhits = NULL;
// ma_ug_t *ug_fa = NULL, *ug_mo = NULL;
opt.coverage_cut = coverage_cut;
opt.sources = sources;
opt.reverse_sources = reverse_sources;
opt.tipsLen = (asm_opt.max_short_tip*2);
opt.tip_drop_ratio = 0.15;
opt.stops_threshold = 3;
opt.ruIndex = ruIndex;
opt.chimeric_rate = 0.05;
opt.drop_ratio = 0.9;
opt.max_hang = max_hang;
opt.min_ovlp = min_ovlp;
opt.is_bench = 0;
opt.b_mask_t = b_mask_t;
opt.gap_fuzz = gap_fuzz;
kvec_asg_arc_t_warp new_rtg_edges, d_edges;
kv_init(new_rtg_edges.a); kv_init(d_edges.a);
ma_ug_t *ug = NULL;
ug = ma_ug_gen_primary(sg, PRIMARY_LABLE);
new_rtg_edges.a.n = 0;
ma_ug_seq(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, &d_edges, 1);///polish
print_debug_gfa(ug, sg, coverage_cut, output_file_name, sources, ruIndex);
ul_resolve(ug, asm_opt.polyploidy);
/**
hap_cov_t *cov = NULL;
trans_chain* t_ch = NULL;
if((asm_opt.flag & HA_F_VERBOSE_GFA)) t_ch = load_hc_trans(output_file_name);
if(!t_ch)
{
new_rtg_edges.a.n = 0;
asg_t *copy_sg = copy_read_graph(sg);
ma_ug_t *copy_ug = copy_untig_graph(ug);
///asm_opt.purge_overlap_len = asm_opt.purge_overlap_len_hic;
///asm_opt.purge_simi_thres = asm_opt.purge_simi_rate_hic;
adjust_utg_by_primary(&copy_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut,
tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio,
max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 0);
print_utg(copy_ug, copy_sg, coverage_cut, output_file_name, sources, ruIndex, max_hang,
min_ovlp, &new_rtg_edges);
if(asm_opt.is_alt)
{
output_contig_graph_alternative(copy_sg, coverage_cut, output_file_name, sources, ruIndex, max_hang,
min_ovlp);
}
ma_ug_destroy(copy_ug);
asg_destroy(copy_sg);
clean_u_trans_t_idx(&(cov->t_ch->k_trans), ug, sg);
new_rtg_edges.a.n = 0;
ma_ug_print_bed(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges,
max_hang, min_ovlp, asm_opt.hic_inconsist_rate, NULL, NULL, cov->t_ch);
if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_trans_chain(cov->t_ch, output_file_name);
}
// char* gfa_name = (char*)malloc(strlen(output_file_name)+50);
// sprintf(gfa_name, "%s.pre.clean_d_utg.noseq.gfa", output_file_name);
// FILE* output_file = fopen(gfa_name, "w");
// ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file);
// fclose(output_file);
// free(gfa_name);
hic_analysis(ug, sg, cov?cov->t_ch:t_ch, &opt, 0, asm_opt.scffold?&rhits:NULL);
if(!rhits && cov) destory_hap_cov_t(&cov);
if(!rhits && t_ch) destory_trans_chain(&t_ch);
// char* gfa_name = (char*)malloc(strlen(output_file_name)+50);
// sprintf(gfa_name, "%s.after.clean_d_utg.noseq.gfa", output_file_name);
// FILE* output_file = fopen(gfa_name, "w");
// ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file);
// fclose(output_file);
// free(gfa_name);
ma_ug_destroy(ug);
kv_destroy(new_rtg_edges.a);
asg_arc_t* av = NULL;
uint32_t v, w, k, i, nv;
for (i = 0; i < d_edges.a.n; i++)
{
v = d_edges.a.a[i].ul>>32;
w = d_edges.a.a[i].v;
av = asg_arc_a(sg, v);
nv = asg_arc_n(sg, v);
for (k = 0; k < nv; k++)
{
if(av[k].del) continue;
if(av[k].v == w)
{
av[k].del = 1;
break;
}
}
}
kv_destroy(d_edges.a);
asg_cleanup(sg);
reduce_hamming_error(sg, sources, coverage_cut, max_hang, min_ovlp, gap_fuzz);
ug_fa = output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex,
0.05, 0.9, max_hang, min_ovlp, rhits?1:0, b_mask_t, NULL);
ug_mo = output_trio_unitig_graph(sg, coverage_cut, output_file_name, MOTHER, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex,
0.05, 0.9, max_hang, min_ovlp, rhits?1:0, b_mask_t, NULL);
if(rhits)
{
ha_aware_order(rhits, sg, ug_fa, ug_mo, cov?&(cov->t_ch->k_trans):&(t_ch->k_trans), &opt, 3);
kv_destroy(rhits->a); kv_destroy(rhits->idx); kv_destroy(rhits->occ); free(rhits);
if(cov) destory_hap_cov_t(&cov);
if(t_ch) destory_trans_chain(&t_ch);
ma_ug_destroy(ug_fa); ma_ug_destroy(ug_mo);
}
**/
}
ma_ug_t *get_poly_ug(asg_t *sg, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
int max_hang, int min_ovlp, R_to_U* ruIndex, bub_label_t* b_mask_t)
{
@@ -30974,6 +31120,12 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
output_poly_trio(sg, coverage_cut, o_file, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex,
0.05, 0.9, max_hang_length, mini_overlap_length, 0, &b_mask_t, asm_opt.polyploidy);
}
else if(asm_opt.ar)
{
if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION;
output_ul_graph(sg, coverage_cut, o_file, sources, reverse_sources, (asm_opt.max_short_tip*2),
0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, gap_fuzz, &b_mask_t);
}
else if (ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt))
{
if(asm_opt.flag & HA_F_PARTITION) asm_opt.flag -= HA_F_PARTITION;
htab.cpp
+129 -2
View File
@@ -509,8 +509,7 @@ void ha_pt_destroy(ha_pt_t *h)
}
if(h->h[i].al){
free(h->h[i].al); h->h[i].al = NULL;
}
}
}
free(h->h); free(h);
}
@@ -634,6 +633,7 @@ int sf##_ha_pt_insert_list(ha_pt_t *h, int n, const HType *a)\
ha_pt1_t *g;\
if (n == 0) return 0;\
g = &h->h[a[0].x&mask];\
/**fprintf(stderr, "a[0].x&mask: %lu, a[0].x&mask: %lu, n: %d\n", a[0].x&mask, a[0].x, n);**/\
for (j = 0; j < n; ++j) {\
uint64_t x = a[j].x >> h->pre;\
khint_t k;\
@@ -643,6 +643,7 @@ int sf##_ha_pt_insert_list(ha_pt_t *h, int n, const HType *a)\
k = yak_pt_get(g->h, x<<YAK_COUNTER_BITS);\
if (k == kh_end(g->h)) continue; \
n = kh_key(g->h, k) & YAK_MAX_COUNT;\
/**fprintf(stderr, "j: %d, n: %d\n", j, n);**/\
assert(n < YAK_MAX_COUNT);\
p = &g->Ia[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;\
@@ -973,6 +974,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_o, int flag, int HPC, int k, int w, h
///for ha_pt_gen, shoud be 0
opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_o->bf_shift;
opt.n_thread = asm_o->thread_num;
// opt.n_thread = p0? 1: asm_o->thread_num;
opt.adaLen = (keep_adapter? asm_o->adapterLen : 0);
opt.min_rcnt = (low_freq?*low_freq:-1);
///asm_opt->num_reads is the number of fastq files
@@ -1083,11 +1085,15 @@ ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_
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);
exit(1);
///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);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> sb0\n", __func__, yak_realtime(), yak_cpu_usage());
pt = ha_pt_gen(ct, asm_opt->thread_num, 1);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> sb1\n", __func__, yak_realtime(), yak_cpu_usage());
ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, pt, flt_tab, NULL, us, 0, NULL);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> sb2\n", __func__, yak_realtime(), yak_cpu_usage());
//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);
@@ -1436,3 +1442,124 @@ int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_op
free(gfa_name);
return 1;
}
int uidx_write(void *flt_tab, ha_pt_t *ha_idx, char* file_name)
{
char* gfa_name = (char*)malloc(strlen(file_name)+25);
sprintf(gfa_name, "%s.uidx.bin", file_name);
FILE* fp = fopen(gfa_name, "w");
if (!fp) {
free(gfa_name);
return 0;
}
yak_ft_t *ha_flt_tab = (yak_ft_t*)flt_tab;
if(ha_flt_tab)
{
fwrite("f", 1, 1, fp);
yak_ft_save(ha_flt_tab, fp);
}
if(ha_idx)
{
int i;
ha_pt1_t *g;
fwrite("h", 1, 1, fp);
fwrite(&ha_idx->k, sizeof(ha_idx->k), 1, fp);
fwrite(&ha_idx->pre, sizeof(ha_idx->pre), 1, fp);
fwrite(&ha_idx->tot, sizeof(ha_idx->tot), 1, fp);
fwrite(&ha_idx->tot_pos, sizeof(ha_idx->tot_pos), 1, fp);
for (i = 0; i < 1<<ha_idx->pre; ++i)
{
g = &(ha_idx->h[i]);
yak_pt_save(g->h, fp);
fwrite(&g->n, sizeof(g->n), 1, fp);
fwrite(g->al, sizeof(ha_idxposl_t), g->n, fp);
}
}
fprintf(stderr, "[M::%s] Index has been written.\n", __func__);
free(gfa_name);
fclose(fp);
return 1;
}
int uidx_load(void **r_flt_tab, ha_pt_t **r_ha_idx, char* file_name)
{
char* gfa_name = (char*)malloc(strlen(file_name)+25);
sprintf(gfa_name, "%s.uidx.bin", file_name);
FILE* fp = fopen(gfa_name, "r");
if (!fp) {
free(gfa_name);
return 0;
}
ha_pt_t *ha_idx = NULL;
char mode = 0;
int f_flag = 0, i;
double index_time, index_s_time, pos_time, pos_s_time;
f_flag += fread(&mode, 1, 1, fp);
if(mode == 'f')
{
index_time = yak_realtime();
yak_ft_load((yak_ft_t **)r_flt_tab, fp);
f_flag += fread(&mode, 1, 1, fp);
fprintf(stderr, "[M::%s::%.3f] ==> Loaded flt table\n", __func__, yak_realtime()-index_time);
}
///insert using multiple threads???
if(mode == 'h')
{
pos_time = index_time = 0;
CALLOC(ha_idx, 1);
ha_pt1_t *g;
f_flag += fread(&ha_idx->k, sizeof(ha_idx->k), 1, fp);
f_flag += fread(&ha_idx->pre, sizeof(ha_idx->pre), 1, fp);
f_flag += fread(&ha_idx->tot, sizeof(ha_idx->tot), 1, fp);
f_flag += fread(&ha_idx->tot_pos, sizeof(ha_idx->tot_pos), 1, fp);
CALLOC(ha_idx->h, 1<<ha_idx->pre);
for (i = 0; i < 1<<ha_idx->pre; ++i)
{
index_s_time = yak_realtime();
g = &(ha_idx->h[i]);
yak_pt_load(&(g->h), fp);
index_time += yak_realtime() - index_s_time;
pos_s_time = yak_realtime();
f_flag += fread(&g->n, sizeof(g->n), 1, fp);
MALLOC(g->al, g->n);
f_flag += fread(g->al, sizeof(ha_idxposl_t), g->n, fp);
pos_time += yak_realtime() - pos_s_time;
}
(*r_ha_idx) = ha_idx;
fprintf(stderr, "[M::%s::%.3f(index)/%.3f(pos)] ==> Loaded pos table\n", __func__, index_time, pos_time);
}
if(mode != 'h' && mode != 'f')
{
free(gfa_name);
fclose(fp);
return 0;
}
fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__);
fclose(fp);
free(gfa_name);
return 1;
}
htab.h
+2
View File
@@ -86,6 +86,8 @@ const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash);
int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name);
int uidx_write(void *flt_tab, ha_pt_t *ha_idx, char* file_name);
int uidx_load(void **r_flt_tab, ha_pt_t **r_ha_idx, char* file_name);
int write_ct_index(void *ct_idx, char* file_name);
int load_ct_index(void **ct_idx, char* file_name);
int query_ct_index(void* ct_idx, uint64_t hash);
inter.cpp
+915 -322
View File
File diff suppressed because it is too large Load Diff
inter.h
+2
View File
@@ -1,5 +1,7 @@
#ifndef __INTER__
#define __INTER__
#include "Overlaps.h"
void ul_resolve(ma_ug_t *ug, int hap_n);
#endif
ksort.h
+11
View File
@@ -52,6 +52,17 @@ typedef struct {
} \
l[i] = tmp; \
} \
void ks_heapup_##name(size_t n, type_t l[]) \
{ \
size_t i, k = n - 1; \
type_t tmp = l[k]; \
while (k) { \
i = (k - 1) >> 1; \
if (__sort_lt(tmp, l[i])) break; \
l[k] = l[i]; k = i; \
} \
l[k] = tmp; \
} \
void ks_heapmake_##name(size_t lsize, type_t l[]) \
{ \
size_t i; \