r215: prepare for list-based binning

This commit is contained in:
Heng Li
2020-04-10 00:12:19 -04:00
parent 7a35bd7fcc
commit f32bfc904a
6 changed files with 50 additions and 33 deletions
+2 -2
View File
@@ -1156,7 +1156,7 @@ int ha_assemble(void)
if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) {
ovlp_loaded = 1;
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage());
if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) trio_partition();
if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt);
if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();
if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
}
@@ -1187,7 +1187,7 @@ int ha_assemble(void)
ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads);
ha_ft_destroy(ha_flt_tab);
if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
trio_partition();
ha_triobin(&asm_opt);
}
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,
+11 -10
View File
@@ -61,14 +61,13 @@ void Print_H(hifiasm_opt_t* asm_opt)
void init_opt(hifiasm_opt_t* asm_opt)
{
memset(asm_opt, 0, sizeof(hifiasm_opt_t));
asm_opt->flag = 0;
asm_opt->coverage = -1;
asm_opt->num_reads = 0;
asm_opt->read_file_names = NULL;
asm_opt->output_file_name = (char*)(DEFAULT_OUTPUT);
asm_opt->required_read_name = NULL;
asm_opt->pat_index = NULL;
asm_opt->mat_index = NULL;
asm_opt->thread_num = 1;
asm_opt->k_mer_length = 51;
asm_opt->mz_win = 51;
@@ -121,7 +120,7 @@ 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);
}
int check_file(char* name, const char* opt)
static int check_file(char* name, const char* opt)
{
if(!name)
{
@@ -256,8 +255,10 @@ int check_option(hifiasm_opt_t* asm_opt)
}
if(asm_opt->pat_index != NULL && check_file(asm_opt->pat_index, "P") == 0) return 0;
if(asm_opt->mat_index != NULL && check_file(asm_opt->mat_index, "M") == 0) return 0;
if(asm_opt->fn_bin_yak[0] != NULL && check_file(asm_opt->fn_bin_yak[0], "YAK1") == 0) return 0;
if(asm_opt->fn_bin_yak[1] != NULL && check_file(asm_opt->fn_bin_yak[1], "YAK2") == 0) return 0;
if(asm_opt->fn_bin_list[0] != NULL && check_file(asm_opt->fn_bin_list[0], "LIST1") == 0) return 0;
if(asm_opt->fn_bin_list[1] != NULL && check_file(asm_opt->fn_bin_list[1], "LIST2") == 0) return 0;
// fprintf(stderr, "input file num: %d\n", asm_opt->num_reads);
// fprintf(stderr, "output file: %s\n", asm_opt->output_file_name);
@@ -272,8 +273,6 @@ int check_option(hifiasm_opt_t* asm_opt)
// fprintf(stderr, "size of popped large bubbles: %lld\n", asm_opt->large_pop_bubble_size);
// fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip);
// fprintf(stderr, "small removed unitig threshold: %d\n", asm_opt->max_short_tip);
// fprintf(stderr, "pat_index: %s\n", asm_opt->pat_index);
// fprintf(stderr, "mat_index: %s\n", asm_opt->mat_index);
// fprintf(stderr, "min_cnt: %d\n", asm_opt->min_cnt);
// fprintf(stderr, "mid_cnt: %d\n", asm_opt->mid_cnt);
@@ -313,7 +312,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
int c;
while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:w:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:1:2:", long_options)) >= 0) {
while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:w:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:1:2:3:4:", long_options)) >= 0) {
if (c == 'h')
{
Print_H(asm_opt);
@@ -339,8 +338,10 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 'b') asm_opt->required_read_name = opt.arg;
else if (c == 'c') asm_opt->min_cnt = atoi(opt.arg);
else if (c == 'd') asm_opt->mid_cnt = atoi(opt.arg);
else if (c == '1' || c == 'P') asm_opt->pat_index = opt.arg; // -P/-M reserved for backward compatibility
else if (c == '2' || c == 'M') asm_opt->mat_index = opt.arg;
else if (c == '1' || c == 'P') asm_opt->fn_bin_yak[0] = opt.arg; // -P/-M reserved for backward compatibility
else if (c == '2' || c == 'M') asm_opt->fn_bin_yak[1] = opt.arg;
else if (c == '3') asm_opt->fn_bin_list[0] = opt.arg;
else if (c == '4') asm_opt->fn_bin_list[1] = opt.arg;
else if (c == 'x') asm_opt->max_drop_rate = atof(opt.arg);
else if (c == 'y') asm_opt->min_drop_rate = atof(opt.arg);
else if (c == 'p') asm_opt->small_pop_bubble_size = atoll(opt.arg);
+8 -3
View File
@@ -3,7 +3,7 @@
#include <pthread.h>
#define HA_VERSION "0.3.0-dirty-r213"
#define HA_VERSION "0.3.0-dirty-r214"
#define VERBOSE 0
@@ -20,8 +20,8 @@ typedef struct {
char** read_file_names;
char* output_file_name;
char* required_read_name;
char* pat_index;
char* mat_index;
char *fn_bin_yak[2];
char *fn_bin_list[2];
int thread_num;
int k_mer_length;
int mz_win;
@@ -66,4 +66,9 @@ void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov);
int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt);
double Get_T(void);
static inline int ha_opt_triobin(const hifiasm_opt_t *opt)
{
return ((opt->fn_bin_yak[0] && opt->fn_bin_yak[1]) || (opt->fn_bin_list[0] && opt->fn_bin_list[1]));
}
#endif
+8 -8
View File
@@ -22786,7 +22786,7 @@ void pre_clean(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *sg, long
///remove isoloated single read
tri_flag += asg_arc_del_single_node_directly(sg, asm_opt.max_short_tip, sources);
if (!((asm_opt.pat_index && asm_opt.mat_index)))
if (!ha_opt_triobin(&asm_opt))
{
tri_flag += asg_arc_del_triangular_advance(sg, bubble_dist);
///remove the cross at the bubble carefully, just remove inexact cross
@@ -24105,7 +24105,7 @@ kvec_t_u32_warp* new_rtg_nodes)
int tri_flag = 0;
///remove very simple circle
tri_flag += asg_arc_del_simple_circle_untig(sources, coverage_cut, r_g, 100, 0);
if (!((asm_opt.pat_index && asm_opt.mat_index)))
if (!ha_opt_triobin(&asm_opt))
{
///remove isoloated single read
tri_flag += asg_arc_del_triangular_advance(r_g, bubble_dist);
@@ -24376,7 +24376,7 @@ kvec_asg_arc_t_warp* new_rtg_edges)
int tri_flag = 0;
///remove very simple circle
tri_flag += asg_arc_del_simple_circle_untig(sources, coverage_cut, r_g, 100, 0);
if (!((asm_opt.pat_index && asm_opt.mat_index)))
if (!ha_opt_triobin(&asm_opt))
{
///remove isoloated single read
tri_flag += asg_arc_del_triangular_advance(r_g, bubble_dist);
@@ -24754,7 +24754,7 @@ uint32_t is_bubble_check, uint32_t is_primary_check)
int tri_flag = 0;
///remove very simple circle
tri_flag += asg_arc_del_simple_circle_untig(sources, coverage_cut, r_g, 100, 0);
if (!((asm_opt.pat_index && asm_opt.mat_index)))
if (!ha_opt_triobin(&asm_opt))
{
///remove isoloated single read
tri_flag += asg_arc_del_triangular_advance(r_g, bubble_dist);
@@ -26312,7 +26312,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
// -1, "clean");
if ((asm_opt.pat_index && asm_opt.mat_index))
if (ha_opt_triobin(&asm_opt))
{
drop_edges_by_trio(sources, n_read);
}
@@ -26401,7 +26401,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
///asg_arc_identify_simple_bubbles_multi(sg, 1);
asg_arc_identify_simple_bubbles_multi(sg, 0);
///asg_arc_del_short_diploid_unclean_exact(sg, drop_ratio, sources);
if((asm_opt.pat_index && asm_opt.mat_index))
if (ha_opt_triobin(&asm_opt))
{
asg_arc_del_short_diploid_by_exact_trio(sg, asm_opt.max_short_tip, sources);
}
@@ -26414,7 +26414,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
asg_arc_identify_simple_bubbles_multi(sg, 1);
if((asm_opt.pat_index && asm_opt.mat_index))
if (ha_opt_triobin(&asm_opt))
{
asg_arc_del_short_diploid_by_length_trio(sg, drop_ratio, asm_opt.max_short_tip, reverse_sources,
asm_opt.max_short_tip, 1, 1, 0, 0, ruIndex);
@@ -26493,7 +26493,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
// rescue_no_coverage_aggressive(sg, sources, reverse_sources, &coverage_cut, ruIndex, max_hang_length,
// mini_overlap_length, bubble_dist, 10);
if ((asm_opt.pat_index && asm_opt.mat_index))
if (ha_opt_triobin(&asm_opt))
{
char *buf = (char*)calloc(strlen(output_file_name) + 25, 1);
sprintf(buf, "%s.dip", output_file_name);
+20 -9
View File
@@ -24,6 +24,9 @@
#define yak_ch_hash(a) ((a)>>YAK_COUNTER_BITS)
KHASHL_SET_INIT(static klib_unused, yak_ht_t, yak_ht, uint64_t, yak_ch_hash, yak_ch_eq)
typedef const char *ha_cstr_t;
KHASHL_MAP_INIT(static klib_unused, cstr_ht_t, cstr_ht, ha_cstr_t, int64_t, kh_hash_str, kh_eq_str)
typedef struct {
struct yak_ht_t *h;
} yak_ch1_t;
@@ -248,23 +251,18 @@ static void tb_worker(void *_data, long k, int tid)
if(type == 'm') aux->seq->trio_flag[k] = MOTHER;
}
void trio_partition(void)
static void ha_triobin_yak(const hifiasm_opt_t *opt)
{
if (asm_opt.pat_index == NULL || asm_opt.mat_index == NULL) {
memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t));
return;
}
yak_ch_t *ch;
int i /**, min_cnt = 2, mid_cnt = 5**/;
tb_shared_t aux;
memset(&aux, 0, sizeof(tb_shared_t));
aux.n_threads = asm_opt.thread_num, aux.print_diff = 0;
aux.n_threads = opt->thread_num, aux.print_diff = 0;
aux.ratio_thres = 0.33;
aux.seq = &R_INF;
ch = yak_ch_restore_core(0, asm_opt.pat_index, YAK_LOAD_TRIOBIN1, asm_opt.min_cnt, asm_opt.mid_cnt);
ch = yak_ch_restore_core(ch, asm_opt.mat_index, YAK_LOAD_TRIOBIN2, asm_opt.min_cnt, asm_opt.mid_cnt);
ch = yak_ch_restore_core(0, opt->fn_bin_yak[0], YAK_LOAD_TRIOBIN1, opt->min_cnt, opt->mid_cnt);
ch = yak_ch_restore_core(ch, opt->fn_bin_yak[1], YAK_LOAD_TRIOBIN2, opt->min_cnt, opt->mid_cnt);
aux.k = ch->k;
aux.ch = ch;
@@ -285,3 +283,16 @@ void trio_partition(void)
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> partitioned reads using yak dumps\n", __func__, yak_realtime(), yak_cpu_usage());
}
static void ha_triobin_list(const hifiasm_opt_t *opt)
{
}
void ha_triobin(const hifiasm_opt_t *opt)
{
memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads * sizeof(uint8_t));
if (opt->fn_bin_list[0] && opt->fn_bin_list[1])
ha_triobin_list(opt);
else if (opt->fn_bin_yak[0] && opt->fn_bin_yak[1])
ha_triobin_yak(opt);
}
+1 -1
View File
@@ -45,7 +45,7 @@ long yak_peakrss(void);
double yak_peakrss_in_gb(void);
double yak_cpu_usage(void);
void trio_partition(void);
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 ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het);