Merge pull request #382 from chhylp123/hifiasm_dev_debug

Hifiasm dev debug
This commit is contained in:
chhylp123
2023-01-17 04:53:24 -05:00
committed by GitHub
10 changed files with 1068 additions and 43 deletions
+1 -29
View File
@@ -1517,31 +1517,6 @@ void Output_PAF()
fprintf(stderr, "PAF has been written.\n");
}
void Output_yak_binning()
{
fprintf(stderr, "Writing binning to disk ...... \n");
char* paf_name = (char*)malloc(strlen(asm_opt.output_file_name)+50);
sprintf(paf_name, "%s.hap1.bin.log", asm_opt.output_file_name);
FILE* oh1 = fopen(paf_name, "w");
sprintf(paf_name, "%s.hap2.bin.log", asm_opt.output_file_name);
FILE* oh2 = fopen(paf_name, "w");
uint64_t i;
for (i = 0; i < R_INF.total_reads; i++) {
if(R_INF.trio_flag[i]==FATHER) {
fprintf(oh1, "%.*s\n", (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
}
if(R_INF.trio_flag[i]==MOTHER) {
fprintf(oh2, "%.*s\n", (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
}
}
free(paf_name);
fclose(oh1); fclose(oh2);
fprintf(stderr, "Binning has been written.\n");
}
int check_cluster(uint64_t* list, long long listLen, ma_hit_t_alloc* paf, float threshold)
{
long long i, k;
@@ -1785,13 +1760,9 @@ int ha_assemble(void)
ha_extract_print_list(&R_INF, asm_opt.extract_iter, asm_opt.extract_list);
exit(0);
}
// if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN) && !(asm_opt.flag & HA_F_VERBOSE_GFA)) ha_triobin(&asm_opt);
if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt);
// if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt), ovlp_loaded = 2;
if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();
if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2;
if (asm_opt.fn_bin_yak[0] && asm_opt.fn_bin_yak[1]) Output_yak_binning();
}
if (!ovlp_loaded) {
ha_flt_tab = ha_idx = NULL;
@@ -1826,6 +1797,7 @@ int ha_assemble(void)
ha_triobin(&asm_opt);
}
if(ovlp_loaded == 2) ovlp_loaded = 0;
ha_opt_update_cov_min(&asm_opt, asm_opt.hom_cov, MIN_N_CHAIN);
build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf,
R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round,
+14 -3
View File
@@ -52,6 +52,7 @@ static ko_longopt_t long_options[] = {
{ "ul-tip", ko_required_argument, 338},
{ "low-het", ko_no_argument, 339},
{ "s-base", ko_required_argument, 340},
{ "bin-only", ko_no_argument, 341},
{ 0, 0, 0 }
};
@@ -193,7 +194,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->max_ov_diff_final = 0.03;
asm_opt->hom_cov = 20;
asm_opt->het_cov = -1024;
asm_opt->max_n_chain = 100;
asm_opt->max_n_chain = MIN_N_CHAIN;
asm_opt->min_hist_kmer_cnt = 5;
asm_opt->load_index_from_disk = 1;
asm_opt->write_index_to_disk = 1;
@@ -260,6 +261,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->is_read_trans = 1;
asm_opt->is_topo_trans = 1;
asm_opt->is_bub_trans = 1;
asm_opt->bin_only = 0;
}
void destory_enzyme(enzyme* f)
@@ -304,6 +306,14 @@ void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov)
fprintf(stderr, "[M::%s] updated max_n_chain to %d\n", __func__, opt->max_n_chain);
}
void ha_opt_update_cov_min(hifiasm_opt_t *opt, int hom_cov, int min_chain)
{
int max_n_chain = (int)(hom_cov * opt->high_factor + .499);
opt->hom_cov = hom_cov; opt->max_n_chain = max_n_chain;
if(opt->max_n_chain < min_chain) opt->max_n_chain = min_chain;
fprintf(stderr, "[M::%s] updated max_n_chain to %d\n", __func__, opt->max_n_chain);
}
static int check_file(char* name, const char* opt)
{
if(!name)
@@ -782,8 +792,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 340) {
asm_opt->trans_base_rate_sec = atof(opt.arg);
if(asm_opt->trans_base_rate_sec < 0) asm_opt->is_base_trans = 0;
} else if (c == 'l')
{ ///0: disable purge_dup; 1: purge containment; 2: purge overlap
}
else if (c == 341) asm_opt->bin_only = 1;
else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap
asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);
}
else if (c == 's') asm_opt->purge_simi_rate_l2 = asm_opt->purge_simi_rate_l3 = atof(opt.arg);
+5 -2
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.18.4-r496"
#define HA_VERSION "0.18.5-r499"
#define VERBOSE 0
@@ -21,9 +21,10 @@
#define HA_F_HIGH_HET 0x400
#define HA_F_PARTITION 0x800
#define HA_F_FAST 0x1000
#define HA_F_USKEW 0x2000
#define HA_F_USKEW 0x2000
#define HA_MIN_OV_DIFF 0.02 // min sequence divergence in an overlap
#define MIN_N_CHAIN 100
typedef struct{
int *l, n;
@@ -135,6 +136,7 @@ typedef struct {
uint8_t is_read_trans;
uint8_t is_topo_trans;
uint8_t is_bub_trans;
uint8_t bin_only;
} hifiasm_opt_t;
extern hifiasm_opt_t asm_opt;
@@ -143,6 +145,7 @@ void init_opt(hifiasm_opt_t* asm_opt);
void destory_opt(hifiasm_opt_t* asm_opt);
void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round);
void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov);
void ha_opt_update_cov_min(hifiasm_opt_t *opt, int hom_cov, int min_chain);
int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt);
double Get_T(void);
+1016 -3
View File
File diff suppressed because it is too large Load Diff
+1 -2
View File
@@ -238,8 +238,7 @@ void print_gfa(asg_t *g);
typedef struct { size_t n, m; uint64_t *a; } asg64_v;
typedef struct { size_t n, m; uint32_t *a; } asg32_v;
typedef struct { size_t n, m; ma_utg_t *a;} ma_utg_v;
typedef struct {
+21
View File
@@ -349,6 +349,27 @@ static void ha_triobin_list(const hifiasm_opt_t *opt)
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> partitioned reads with external lists\n", __func__, yak_realtime(), yak_cpu_usage());
}
uint32_t test_yak_binning(char* fn, char *cmd)
{
gzFile fp; kstream_t *ks; kstring_t str = {0,0,0};
int dret, eq = 0; fp = gzopen(fn, "r");
if (fp == 0) return eq;
ks = ks_init(fp);
if(ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) {
uint64_t sl = strlen(str.s), z;
for (z = 0; (z < sl) && (str.s[z]!='\t'); z++);
if(((z+1)<sl) && (str.s[z]=='\t')) {
if((strlen(str.s+z+1)==strlen(cmd)) && (!memcmp(str.s+z+1, cmd, strlen(cmd)))) {
eq = 1;
}
}
}
free(str.s);
ks_destroy(ks);
gzclose(fp);
return eq;
}
inline void phrase_hstatus(char *s, char **rname, uint32_t *hid)
{
char *p = NULL, *id = NULL; *rname = NULL; *hid = (uint32_t)-1;
+4 -2
View File
@@ -31,7 +31,8 @@ KRADIX_SORT_INIT(osg, osg_arc_t, osg_arc_key, member_size(osg_arc_t, u))
#define BREAK_THRES 5000000
#define BREAK_CUTOFF 0.1
#define BREAK_BOUNDARY 0.015
void reduce_hamming_error_adv(ma_ug_t *iug, asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
int max_hang, int min_ovlp, long long gap_fuzz, R_to_U *ru, bubble_type* bub);
typedef struct {
uint64_t ruid;
@@ -3962,7 +3963,8 @@ asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, kv_u_trans_t *ref, ug_opt_t *opt,
t_idx = init_trans_col(i_ug, h->r_g->n_seq, ref);
// output_hic_rtg(i_ug, h->r_g, opt, asm_opt.output_file_name);
reduce_hamming_error(h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz);
// reduce_hamming_error(h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz);
reduce_hamming_error_adv(NULL, h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz, opt->ruIndex, NULL);
/**
scaffold_hap(h, t_idx, opt, round, asm_opt.output_file_name, FATHER);
scaffold_hap(h, t_idx, opt, round, asm_opt.output_file_name, MOTHER);
+1
View File
@@ -1546,6 +1546,7 @@ int uidx_load(void **r_flt_tab, ha_pt_t **r_ha_idx, char* file_name, ma_ug_t *ug
free(gfa_name); fclose(fp);
return 0;
}
// fprintf(stderr, "[M::%s]\t%s\tftell::%ld\n", __func__, file_name, ftell(fp));
ha_pt_t *ha_idx = NULL;
char mode = 0;
+1
View File
@@ -114,6 +114,7 @@ double yak_peakrss_in_gb(void);
double yak_cpu_usage(void);
void ha_triobin(const hifiasm_opt_t *opt);
uint32_t test_yak_binning(char* fn, char *cmd);
uint32_t *ha_polybin_list(const hifiasm_opt_t *opt);
void mz1_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, void *km);
+4 -2
View File
@@ -10449,7 +10449,7 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac
}
// fprintf(stderr, "[M::%s::Start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
kt_for(p->n_thread, worker_for_ul_scall_alignment, s, s->n);
fprintf(stderr, "[M::%s::Done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
// fprintf(stderr, "[M::%s::Done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
///debug
/**
uint64_t i;
@@ -10580,7 +10580,7 @@ static void *worker_ul_rescall_pipeline(void *data, int step, void *in) // callb
}
// fprintf(stderr, "[M::%s::Start] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
kt_for(p->n_thread, worker_for_ul_rescall_alignment, s, s->n);
fprintf(stderr, "[M::%s::Done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
// fprintf(stderr, "[M::%s::Done] ==> s->id: %lu, s->n:% d\n", __func__, s->id, s->n);
// get_utepdat_t_mem(s, 1);
for (i = 0; i < p->n_thread; ++i) {
@@ -15594,6 +15594,7 @@ void gen_UL_ovlps(uldat_t *sl, int32_t cutoff)
ul_idx_t *uu = dedup_HiFis(sl->uopt, 1, 0);
int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name, NULL) : 0);
if(exist == 0) uidx_l_build(uu->ug, (mg_idxopt_t *)sl->opt, cutoff);
// print_debug_gfa(sl->uopt, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len);
if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name, NULL);
sl->ha_flt_tab = ha_flt_tab; sl->ha_idx = (ha_pt_t *)ha_idx; sl->uu = uu;
ul_v_call(sl, asm_opt.ar);
@@ -15699,6 +15700,7 @@ int32_t write_all_ul_t(all_ul_t *x, char* file_name, ma_ug_t *ug)
fprintf(stderr, "[M::%s] Index has been written.\n", __func__);
fclose(fp);
if(asm_opt.bin_only) exit(1);
return 1;
}