From fda858f9e78d612654ab111c1ce2dab6540452bd Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Mon, 13 Jul 2020 22:34:26 -0400 Subject: [PATCH] update coverage --- Assembly.cpp | 2 +- CommandLines.cpp | 87 ++++++---- CommandLines.h | 2 +- Overlaps.cpp | 406 ++++++++++++++++++++++++++++++++++++++++------- Overlaps.h | 4 +- Purge_Dups.cpp | 37 +++-- Purge_Dups.h | 2 +- hifiasm.1 | 13 +- 8 files changed, 449 insertions(+), 104 deletions(-) diff --git a/Assembly.cpp b/Assembly.cpp index 9479b1a..9671816 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1203,7 +1203,7 @@ 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), ovlp_loaded = 2; + 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), ovlp_loaded = 2; if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads(); if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF(); diff --git a/CommandLines.cpp b/CommandLines.cpp index f0265a6..141915e 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -19,6 +19,8 @@ static ko_longopt_t long_options[] = { { "max-od-final", ko_no_argument, 306 }, { "ex-list", ko_required_argument, 307 }, { "ex-iter", ko_required_argument, 308 }, + { "purge-cov", ko_required_argument, 309 }, + { "pri-range", ko_required_argument, 310 }, { 0, 0, 0 } }; @@ -34,41 +36,45 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, "Usage: hifiasm [options] <...>\n"); fprintf(stderr, "Options:\n"); fprintf(stderr, " Assembly:\n"); - fprintf(stderr, " -o FILE prefix of output files [%s]\n", asm_opt->output_file_name); - fprintf(stderr, " -t INT number of threads [%d]\n", asm_opt->thread_num); - fprintf(stderr, " -r INT round of correction [%d]\n", asm_opt->number_of_round); - fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round); - fprintf(stderr, " -k INT k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length); - fprintf(stderr, " -w INT minimizer window size [%d]\n", asm_opt->mz_win); - fprintf(stderr, " -f INT number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift); - fprintf(stderr, " -D FLOAT drop k-mers occuring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor); - fprintf(stderr, " -N INT consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain); - fprintf(stderr, " -i ignore saved overlaps in *.ovlp* files\n"); - fprintf(stderr, " -z INT length of adapters that should be removed [%d]\n", asm_opt->adapterLen); - fprintf(stderr, " -m INT size of popped large bubbles for contig graph [%lld]\n", asm_opt->large_pop_bubble_size); - fprintf(stderr, " -p INT size of popped small bubbles for haplotype-resolved unitig graph [%lld]\n", asm_opt->small_pop_bubble_size); - fprintf(stderr, " -n INT small removed unitig threshold [%d]\n", asm_opt->max_short_tip); - fprintf(stderr, " -x FLOAT max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate); - fprintf(stderr, " -y FLOAT min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate); - fprintf(stderr, " -u disable post join contigs step which may improve N50. Don't disable in default.\n"); - fprintf(stderr, " --version show version number\n"); - fprintf(stderr, " -h show help information\n"); + fprintf(stderr, " -o FILE prefix of output files [%s]\n", asm_opt->output_file_name); + fprintf(stderr, " -t INT number of threads [%d]\n", asm_opt->thread_num); + fprintf(stderr, " -r INT round of correction [%d]\n", asm_opt->number_of_round); + fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round); + fprintf(stderr, " -k INT k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length); + fprintf(stderr, " -w INT minimizer window size [%d]\n", asm_opt->mz_win); + fprintf(stderr, " -f INT number of bits for bloom filter; 0 to disable [%d]\n", asm_opt->bf_shift); + fprintf(stderr, " -D FLOAT drop k-mers occuring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor); + fprintf(stderr, " -N INT consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain); + fprintf(stderr, " -i ignore saved overlaps in *.ovlp* files\n"); + fprintf(stderr, " -z INT length of adapters that should be removed [%d]\n", asm_opt->adapterLen); + fprintf(stderr, " -m INT size of popped large bubbles for contig graph [%lld]\n", asm_opt->large_pop_bubble_size); + fprintf(stderr, " -p INT size of popped small bubbles for haplotype-resolved unitig graph [%lld]\n", asm_opt->small_pop_bubble_size); + fprintf(stderr, " -n INT small removed unitig threshold [%d]\n", asm_opt->max_short_tip); + fprintf(stderr, " -x FLOAT max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate); + fprintf(stderr, " -y FLOAT min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate); + fprintf(stderr, " --pri-range INT,INT keep contigs with coverage in this range at p_ctg.gfa. Inferred automatically in default. Set -1,-1 to disable\n"); + + + fprintf(stderr, " -u disable post join contigs step which may improve N50. Don't disable in default.\n"); + fprintf(stderr, " --version show version number\n"); + fprintf(stderr, " -h show help information\n"); fprintf(stderr, " Trio-partition:\n"); - fprintf(stderr, " -1 FILE hap1/paternal k-mer dump generated by \"yak count\" []\n"); - fprintf(stderr, " -2 FILE hap2/maternal k-mer dump generated by \"yak count\" []\n"); - fprintf(stderr, " -3 FILE list of hap1/paternal read names []\n"); - fprintf(stderr, " -4 FILE list of hap2/maternal read names []\n"); - fprintf(stderr, " -c INT lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt); - fprintf(stderr, " -d INT upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt); + fprintf(stderr, " -1 FILE hap1/paternal k-mer dump generated by \"yak count\" []\n"); + fprintf(stderr, " -2 FILE hap2/maternal k-mer dump generated by \"yak count\" []\n"); + fprintf(stderr, " -3 FILE list of hap1/paternal read names []\n"); + fprintf(stderr, " -4 FILE list of hap2/maternal read names []\n"); + fprintf(stderr, " -c INT lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt); + fprintf(stderr, " -d INT upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt); fprintf(stderr, " Purge-dups:\n"); - fprintf(stderr, " -l INT level of purge-dup. In default, [%d] for non-trio; [%d] for trio (see hifiasm.1 for details)\n", + fprintf(stderr, " -l INT level of purge-dup. In default, [%d] for non-trio; [%d] for trio (see hifiasm.1 for details)\n", asm_opt->purge_level_primary, asm_opt->purge_level_trio); - fprintf(stderr, " -s FLOAT similarity threshold for duplicate haplotigs [%g]\n", + fprintf(stderr, " -s FLOAT similarity threshold for duplicate haplotigs [%g]\n", asm_opt->purge_simi_rate); - fprintf(stderr, " -O INT min number of overlapped reads for duplicate haplotigs [%d]\n", + fprintf(stderr, " -O INT min number of overlapped reads for duplicate haplotigs [%d]\n", asm_opt->purge_overlap_len); + fprintf(stderr, " --purge-cov INT coverage upper bound of Purge-dups, which is inferred automatically in default.\n"); fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n"); @@ -117,8 +123,8 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->purge_level_trio = 0; asm_opt->purge_simi_rate = 0.75; asm_opt->purge_overlap_len = 1; - asm_opt->recover_atg_cov_min = -1; - asm_opt->recover_atg_cov_max = -1; + asm_opt->recover_atg_cov_min = -1024; + asm_opt->recover_atg_cov_max = -1024; asm_opt->hom_global_coverage = -1; } @@ -303,6 +309,19 @@ int check_option(hifiasm_opt_t* asm_opt) return 0; } + if(asm_opt->hom_global_coverage < 0 && asm_opt->hom_global_coverage != -1) + { + fprintf(stderr, "[ERROR] purge duplication coverage threshold should be >= 0 (--purge-cov)\n"); + return 0; + } + + if(asm_opt->recover_atg_cov_min != asm_opt->recover_atg_cov_max && + (asm_opt->recover_atg_cov_min == -1024 || asm_opt->recover_atg_cov_max == -1024)) + { + fprintf(stderr, "[ERROR] primary contig coverage range should be set correctly (--primary-range)\n"); + return 0; + } + if(asm_opt->fn_bin_yak[0] != NULL && check_file(asm_opt->fn_bin_yak[0], "YAK1") == 0) return 0; @@ -411,6 +430,13 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 306) asm_opt->max_ov_diff_final = atof(opt.arg); else if (c == 307) asm_opt->extract_list = opt.arg; else if (c == 308) asm_opt->extract_iter = atoi(opt.arg); + else if (c == 309) asm_opt->hom_global_coverage = atoi(opt.arg); + else if (c == 310) + { + char* s = NULL; + asm_opt->recover_atg_cov_min = strtol(opt.arg, &s, 10); + if (*s == ',') asm_opt->recover_atg_cov_max = strtol(s + 1, &s, 10); + } 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); @@ -429,6 +455,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) } } + if (argc == opt.ind) { Print_H(asm_opt); diff --git a/CommandLines.h b/CommandLines.h index 60ee2b1..7395673 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.7-dirty-r256" +#define HA_VERSION "0.8-dirty-r280" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index c8d3103..d8f8c9f 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -9383,17 +9383,52 @@ ma_hit_t_alloc* sources, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp) } +uint32_t get_ug_coverage(ma_utg_t* u, asg_t* read_g, const ma_sub_t* coverage_cut, +ma_hit_t_alloc* sources, R_to_U* ruIndex) +{ + uint32_t k, j, rId, tn, is_Unitig; + long long R_bases = 0, C_bases = 0; + ma_hit_t *h; + if(u->m == 0) return 0; -void ma_ug_print2(const ma_ug_t *ug, All_reads *RNF, const ma_sub_t *coverage_cut, int print_seq, FILE *fp) + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + R_bases += (coverage_cut[rId].e - coverage_cut[rId].s); + for (j = 0; j < (uint64_t)(sources[rId].length); j++) + { + h = &(sources[rId].buffer[j]); + if(h->el != 1) continue; + tn = Get_tn((*h)); + if(read_g->seq[tn].del == 1) + { + ///get the id of read that contains it + get_R_to_U(ruIndex, tn, &tn, &is_Unitig); + if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue; + } + if(read_g->seq[tn].del == 1) continue; + C_bases += (Get_qe((*h)) - Get_qs((*h))); + } + } + + return C_bases/R_bases; +} + +void ma_ug_print2(const ma_ug_t *ug, All_reads *RNF, asg_t* read_g, const ma_sub_t *coverage_cut, +ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FILE *fp) { uint32_t i, j, l; char name[32]; for (i = 0; i < ug->u.n; ++i) { // the Segment lines in GFA ma_utg_t *p = &ug->u.a[i]; if(p->m == 0) continue; - sprintf(name, "utg%.6d%c", i + 1, "lc"[p->circ]); - if (print_seq) fprintf(fp, "S\t%s\t%s\tLN:i:%d\n", name, p->s? p->s : "*", p->len); - else fprintf(fp, "S\t%s\t*\tLN:i:%d\n", name, p->len); + sprintf(name, "%s%.6d%c", prefix, i + 1, "lc"[p->circ]); + if (print_seq) fprintf(fp, "S\t%s\t%s\tLN:i:%d\trd:i:%u\n", name, p->s? p->s : "*", p->len, + get_ug_coverage(p, read_g, coverage_cut, sources, ruIndex)); + else fprintf(fp, "S\t%s\t*\tLN:i:%d\trd:i:%u\n", name, p->len, + get_ug_coverage(p, read_g, coverage_cut, sources, ruIndex)); + // if (print_seq) fprintf(fp, "S\t%s\t%s\tLN:i:%d\n", name, p->s? p->s : "*", p->len); + // else fprintf(fp, "S\t%s\t*\tLN:i:%d\n", name, p->len); for (j = l = 0; j < p->n; l += (uint32_t)p->a[j++]) { uint32_t x = p->a[j]>>33; @@ -9414,14 +9449,16 @@ void ma_ug_print2(const ma_ug_t *ug, All_reads *RNF, const ma_sub_t *coverage_cu } for (i = 0; i < ug->g->n_arc; ++i) { // the Link lines in GFA uint32_t u = ug->g->arc[i].ul>>32, v = ug->g->arc[i].v; - fprintf(fp, "L\tutg%.6d%c\t%c\tutg%.6d%c\t%c\t%dM\tL1:i:%d\n", (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], - (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], ug->g->arc[i].ol, asg_arc_len(ug->g->arc[i])); + fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\n", + prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], + prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], ug->g->arc[i].ol, asg_arc_len(ug->g->arc[i])); } } -void ma_ug_print(const ma_ug_t *ug, All_reads *RNF, const ma_sub_t *coverage_cut, FILE *fp) +void ma_ug_print(const ma_ug_t *ug, All_reads *RNF, asg_t* read_g, const ma_sub_t *coverage_cut, +ma_hit_t_alloc* sources, R_to_U* ruIndex, const char* prefix, FILE *fp) { - ma_ug_print2(ug, RNF, coverage_cut, 1, fp); + ma_ug_print2(ug, RNF, read_g, coverage_cut, sources, ruIndex, 1, prefix, fp); } int asg_cut_internal(asg_t *g, int max_ext) @@ -9666,9 +9703,10 @@ void ma_sg_print(const asg_t *g, const All_reads *RNF, const ma_sub_t *sub, FILE } -void ma_ug_print_simple(const ma_ug_t *ug, All_reads *RNF, const ma_sub_t *coverage_cut, FILE *fp) +void ma_ug_print_simple(const ma_ug_t *ug, All_reads *RNF, asg_t* read_g, const ma_sub_t *coverage_cut, +ma_hit_t_alloc* sources, R_to_U* ruIndex, const char* prefix, FILE *fp) { - ma_ug_print2(ug, RNF, coverage_cut, 0, fp); + ma_ug_print2(ug, RNF, read_g, coverage_cut, sources, ruIndex, 0, prefix, fp); } @@ -11134,7 +11172,7 @@ R_to_U* ruIndex) void output_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, -ma_hit_t_alloc* sources, int max_hang, int min_ovlp) +ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) { kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); @@ -11147,11 +11185,11 @@ ma_hit_t_alloc* sources, int max_hang, int min_ovlp) char* gfa_name = (char*)malloc(strlen(output_file_name)+25); sprintf(gfa_name, "%s.r_utg.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); - ma_ug_print(ug, &R_INF, coverage_cut, output_file); + ma_ug_print(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); fclose(output_file); sprintf(gfa_name, "%s.r_utg.noseq.gfa", output_file_name); output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); + ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); fclose(output_file); free(gfa_name); @@ -12119,15 +12157,94 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_ return n_reduced; } -int magic_trio_phasing(asg_t *g, ma_ug_t *ug, asg_t *read_sg, /**ma_hit_t_alloc* sources,**/ -ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, -float drop_rate) +uint8_t if_primary_unitig(ma_utg_t* u, asg_t* read_g, ma_sub_t* coverage_cut, +ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag) +{ + if(asm_opt.recover_atg_cov_min < 0 || asm_opt.recover_atg_cov_max < 0) + { + return 0; + } + long long R_bases = 0, C_bases = 0, C_bases_primary = 0, C_bases_alter = 0; + uint32_t available_reads = 0, k, j, rId, tn, is_Unitig; + ma_hit_t *h; + if(u->m == 0) return 0; + available_reads = 0; + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + r_flag[rId] = 1; + } + + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + C_bases = C_bases_primary = C_bases_alter = 0; + R_bases = coverage_cut[rId].e - coverage_cut[rId].s; + for (j = 0; j < (uint64_t)(sources[rId].length); j++) + { + h = &(sources[rId].buffer[j]); + if(h->el != 1) continue; + tn = Get_tn((*h)); + if(read_g->seq[tn].del == 1) + { + ///get the id of read that contains it + get_R_to_U(ruIndex, tn, &tn, &is_Unitig); + if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue; + } + if(read_g->seq[tn].del == 1) continue; + if(r_flag[tn]) + { + C_bases_primary += Get_qe((*h)) - Get_qs((*h)); + } + else + { + C_bases_alter += Get_qe((*h)) - Get_qs((*h)); + } + } + + ///fprintf(stderr, "C_bases_primary: %lld, C_bases_alter: %lld\n", C_bases_primary, C_bases_alter); + + C_bases = C_bases_primary + C_bases_alter; + if(C_bases_primary < C_bases * ALTER_COV_THRES) continue; + + C_bases = C_bases/R_bases; + if(C_bases >= asm_opt.recover_atg_cov_min && C_bases <= asm_opt.recover_atg_cov_max) + { + available_reads++; + } + } + + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + r_flag[rId] = 0; + } + + ///fprintf(stderr, "available_reads: %u, u->n: %u\n", available_reads, u->n); + + if(available_reads < (u->n * 0.8) || available_reads == 0) + { + return 0; + } + else + { + ///fprintf(stderr, "*****************\n"); + return 1; + } +} + + + +int magic_trio_phasing(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_sub_t* coverage_cut, +ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long miniedgeLen, +R_to_U* ruIndex, uint32_t positive_flag, float drop_rate) { uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0; ma_utg_t* nsu = NULL; uint32_t flag = (uint32_t)-1, flag_occ, non_flag_occ, ambigious, del_node, keep_node; if(positive_flag == FATHER) flag = MOTHER; if(positive_flag == MOTHER) flag = FATHER; + uint8_t* primary_flag = (uint8_t*)calloc(read_sg->n_seq, sizeof(uint8_t)); for (v = 0; v < n_vtx; ++v) { @@ -12167,6 +12284,11 @@ float drop_rate) continue; } + if(if_primary_unitig(nsu, read_sg, coverage_cut, sources, ruIndex, primary_flag)) + { + continue; + } + g->seq[av[i].v>>1].c = ALTER_LABLE; asg_seq_drop(g, av[i].v>>1); n_reduced++; @@ -12177,6 +12299,7 @@ float drop_rate) asg_cleanup(g); + free(primary_flag); return n_reduced; } @@ -12622,8 +12745,9 @@ ma_ug_t* copy_untig_graph(ma_ug_t *src) return ug; } -void clean_trio_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, -long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, +void clean_trio_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, +ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist, +long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen, uint32_t miniBiGraph, float chimeric_rate, int is_final_clean, int just_bubble_pop, float drop_ratio, uint32_t trio_flag, float trio_drop_rate) @@ -12635,7 +12759,7 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate) ///if(trio_flag == MOTHER) fprintf(stderr, "(0) c: %u, del: %u, n: %u\n", ug->g->seq[28141].c, ug->g->seq[28141].del, ug->u.a[28141].n); asg_pop_bubble_primary_trio(ug, bubble_dist, trio_flag, DROP); untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, trio_flag, DROP); - magic_trio_phasing(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate); + magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate); ///drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex); /**********debug**********/ @@ -12693,7 +12817,7 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate) resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, trio_flag, drop_ratio); drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex); - all_to_all_deduplicate(ug, trio_flag, trio_drop_rate, reverse_sources, ruIndex); + all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex); if(is_first) { @@ -12802,12 +12926,14 @@ void set_drop_trio_flag(ma_ug_t *ug) -void update_unitig_graph(ma_ug_t* ug, asg_t* read_g, ma_hit_t_alloc* reverse_sources, -R_to_U* ruIndex, uint8_t is_double_check, uint8_t flag, float drop_rate) +void update_unitig_graph(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, +ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, +uint8_t is_double_check, uint8_t flag, float drop_rate) { asg_t* nsg = ug->g; uint32_t v, n_vtx = nsg->n_seq, k, rId, flag_occ, non_flag_occ, hap_label_occ, n_reduce = 1; ma_utg_t *u; + uint8_t* primary_flag = (uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t)); drop_semi_circle(ug, nsg, read_g, reverse_sources, ruIndex); @@ -12839,6 +12965,10 @@ R_to_U* ruIndex, uint8_t is_double_check, uint8_t flag, float drop_rate) if(non_flag_occ > ((non_flag_occ+flag_occ)*drop_rate)) { + if(if_primary_unitig(u, read_g, coverage_cut, sources, ruIndex, primary_flag)) + { + continue; + } if(u->m != 0) { u->circ = u->end = u->len = u->m = u->n = u->start = 0; @@ -12853,6 +12983,7 @@ R_to_U* ruIndex, uint8_t is_double_check, uint8_t flag, float drop_rate) drop_semi_circle(ug, nsg, read_g, reverse_sources, ruIndex); asg_cleanup(nsg); + free(primary_flag); } @@ -12974,8 +13105,10 @@ uint32_t* non_require, uint32_t* ambigious) } } + ///note: to use this function, don't renew unitig graph!!!!!!!!! -void all_to_all_deduplicate(ma_ug_t* ug, uint8_t postive_flag, float drop_rate, +void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, +ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex) { @@ -12988,6 +13121,7 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex) nsg = ug->g; if(postive_flag == FATHER) flag = MOTHER; if(postive_flag == MOTHER) flag = FATHER; + uint8_t* primary_flag = (uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t)); for (v = 0; v < nsg->n_seq; v++) { @@ -13062,6 +13196,10 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex) if(k != u_vecs.a.n) { + if(if_primary_unitig(nsu, read_g, coverage_cut, sources, ruIndex, primary_flag)) + { + continue; + } if(nsu->m != 0) { nsu->circ = nsu->end = nsu->len = nsu->m = nsu->n = nsu->start = 0; @@ -13081,6 +13219,7 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex) } kv_destroy(u_vecs.a); + free(primary_flag); } void delete_useless_nodes(ma_ug_t **ug) @@ -13124,6 +13263,61 @@ void delete_useless_nodes(ma_ug_t **ug) asg_cleanup(nsg); } + + + +void delete_useless_trio_nodes(ma_ug_t **ug, asg_t* read_g, ma_sub_t* coverage_cut, +ma_hit_t_alloc* sources, R_to_U* ruIndex) +{ + asg_t* nsg = (*ug)->g; + uint32_t v, n_vtx = nsg->n_seq, convex; + long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen; + uint8_t* primary_flag = (uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t)); + + + for (v = 0; v < n_vtx; ++v) + { + if(nsg->seq[v].del) continue; + if(nsg->seq[v].c == ALTER_LABLE && + (if_primary_unitig(&((*ug)->u.a[v]), read_g, coverage_cut, sources, + ruIndex, primary_flag) == 0)) + { + asg_seq_del(nsg, v); + if((*ug)->u.a[v].m!=0) + { + (*ug)->u.a[v].m = (*ug)->u.a[v].n = 0; + free((*ug)->u.a[v].a); + (*ug)->u.a[v].a = NULL; + } + + continue; + } + + //note: after cleaning, some cirle might be gone, or we have some new circles + //so need to renew .circ + if(get_unitig(nsg, NULL, (v<<1), &convex, &nodeLen, &baseLen, + &max_stop_nodeLen, &max_stop_baseLen, 1, NULL)==LOOP) + { + (*ug)->u.a[v].circ = 1; + (*ug)->u.a[v].start = (*ug)->u.a[v].end = UINT32_MAX; + } + else + { + (*ug)->u.a[v].circ = 0; + + (*ug)->u.a[v].start = (*ug)->u.a[v].a[0]>>32; + (*ug)->u.a[v].end = ((*ug)->u.a[v].a[(*ug)->u.a[v].n-1]>>32)^1; + } + } + + asg_cleanup(nsg); + free(primary_flag); +} + + + + + void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex) { uint32_t v, n_vtx = nsg->n_seq*2, convex_f, convex_b, i, nv; @@ -13210,6 +13404,87 @@ void update_hap_label(ma_ug_t *ug, asg_t* read_g) } +uint8_t* get_utg_attributes(ma_ug_t *ug, asg_t* read_g, ma_sub_t* coverage_cut, +ma_hit_t_alloc* sources, R_to_U* ruIndex) +{ + asg_t* nsg = ug->g; + uint32_t v, n_vtx = nsg->n_seq, k, j, rId, available_reads = 0, tn, is_Unitig; + ma_utg_t* u = NULL; + ma_hit_t *h; + long long R_bases = 0, C_bases = 0, C_bases_primary = 0, C_bases_alter = 0; + uint8_t* r_flag = (uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t)); + uint8_t* u_flag = (uint8_t*)calloc(n_vtx, sizeof(uint8_t)); + for (v = 0; v < n_vtx; ++v) + { + if(nsg->seq[v].del) continue; + u = &(ug->u.a[v]); + if(u->m == 0) continue; + available_reads = 0; + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + r_flag[rId] = 1; + } + + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + C_bases = C_bases_primary = C_bases_alter = 0; + R_bases = coverage_cut[rId].e - coverage_cut[rId].s; + for (j = 0; j < (uint64_t)(sources[rId].length); j++) + { + h = &(sources[rId].buffer[j]); + if(h->el != 1) continue; + tn = Get_tn((*h)); + if(read_g->seq[tn].del == 1) + { + ///get the id of read that contains it + get_R_to_U(ruIndex, tn, &tn, &is_Unitig); + if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue; + } + if(read_g->seq[tn].del == 1) continue; + if(r_flag[tn]) + { + C_bases_primary += Get_qe((*h)) - Get_qs((*h)); + } + else + { + C_bases_alter += Get_qe((*h)) - Get_qs((*h)); + } + } + + C_bases = C_bases_primary + C_bases_alter; + if(C_bases_alter < C_bases * ALTER_COV_THRES) continue; + + C_bases = C_bases/R_bases; + if(C_bases >= asm_opt.recover_atg_cov_min && C_bases <= asm_opt.recover_atg_cov_max) + { + available_reads++; + } + } + + + if(available_reads < (u->n * 0.8) || available_reads == 0) + { + u_flag[v] = 0; + } + else + { + u_flag[v] = 1; + } + + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + r_flag[rId] = 0; + } + } + + free(r_flag); + return u_flag; +} + + void adjust_utg_by_trio(ma_ug_t **ug, asg_t* read_g, uint8_t flag, float drop_rate, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, @@ -13217,8 +13492,7 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov kvec_asg_arc_t_warp* new_rtg_edges) { asg_t* nsg = (*ug)->g; - uint32_t v, n_vtx = nsg->n_seq; - + uint32_t v, n_vtx = nsg->n_seq; //if(flag == MOTHER) print_untig_by_read(*ug, "m54329U_190617_231905/65340614/ccs", (uint32_t)-1, sources, reverse_sources, "beg1"); //if(flag == MOTHER) print_untig_by_read(*ug, "m54329U_190827_173812/67108993/ccs", (uint32_t)-1, sources, reverse_sources, "beg2"); @@ -13226,7 +13500,22 @@ kvec_asg_arc_t_warp* new_rtg_edges) kvec_t_u32_warp new_rtg_nodes; kv_init(new_rtg_nodes.a); **/ - update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, 0, flag, drop_rate); + purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, + asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, bubble_dist, + drop_ratio, 1, 1); + if(asm_opt.recover_atg_cov_min == -1024 || asm_opt.recover_atg_cov_max == -1024) + { + asm_opt.recover_atg_cov_max = asm_opt.hom_global_coverage/HOM_PEAK_RATE; + asm_opt.recover_atg_cov_min = asm_opt.recover_atg_cov_max * 0.8; + asm_opt.recover_atg_cov_max = asm_opt.recover_atg_cov_max * 1.2; + } + fprintf(stderr, "[M::%s] primary contig coverage range: [%d, %d]\n", + __func__, asm_opt.recover_atg_cov_min, asm_opt.recover_atg_cov_max); + + ///primary_flag = get_utg_attributes(*ug, read_g, coverage_cut, sources, ruIndex); + + + update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0, flag, drop_rate); adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex); nsg = (*ug)->g; n_vtx = nsg->n_seq; @@ -13237,17 +13526,18 @@ kvec_asg_arc_t_warp* new_rtg_edges) EvaluateLen((*ug)->u, v) = (*ug)->u.a[v].n; } - clean_trio_untig_graph(*ug, read_g, reverse_sources, bubble_dist, tipsLen, - tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, + clean_trio_untig_graph(*ug, read_g, coverage_cut, sources, reverse_sources, bubble_dist, + tipsLen, tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, chimeric_rate, 0, 0, drop_ratio, flag, drop_rate); ///if(flag == MOTHER) fprintf(stderr, "(o.1) c: %u, del: %u, n: %u\n", (*ug)->g->seq[28141].c, (*ug)->g->seq[28141].del, (*ug)->u.a[28141].n); - delete_useless_nodes(ug); + ///delete_useless_nodes(ug); + delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex); update_hap_label(*ug, read_g); - update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, 0, flag, drop_rate); + update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0, flag, drop_rate); renew_utg(ug, read_g, new_rtg_edges); @@ -13264,22 +13554,25 @@ kvec_asg_arc_t_warp* new_rtg_edges) renew_utg(ug, read_g, new_rtg_edges); } - update_unitig_graph((*ug), read_g, reverse_sources, ruIndex, 1, flag, drop_rate); + update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 1, flag, drop_rate); update_hap_label(NULL, read_g); renew_utg(ug, read_g, new_rtg_edges); - delete_useless_nodes(ug); + ///delete_useless_nodes(ug); + delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex); if(asm_opt.purge_level_trio == 1) { purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, bubble_dist, - drop_ratio, 1); - delete_useless_nodes(ug); + drop_ratio, 1, 0); + ///delete_useless_nodes(ug); + delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex); } + set_drop_trio_flag(*ug); /** @@ -13330,12 +13623,12 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp) ///debug_untig_length(ug, tipsLen, gfa_name); ///print_untig_by_read(ug, "m64011_190901_095311/125831121/ccs", 2310925, "end"); ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp); - ma_ug_print(ug, &R_INF, coverage_cut, output_file); + ma_ug_print(ug, &R_INF, sg, coverage_cut, sources, ruIndex, (flag==FATHER?"h1tg":"h2tg"), output_file); fclose(output_file); sprintf(gfa_name, "%s.%s.p_ctg.noseq.gfa", output_file_name, (flag==FATHER?"hap1":"hap2")); output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); + ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, (flag==FATHER?"h1tg":"h2tg"), output_file); fclose(output_file); free(gfa_name); @@ -22231,7 +22524,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex) } } - fprintf(stderr, "keep_atg: %u\n", keep_atg); + ///fprintf(stderr, "keep_atg: %u\n", keep_atg); ma_ug_destroy(atg); } @@ -22272,7 +22565,6 @@ kvec_asg_arc_t_warp* new_rtg_edges) tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, chimeric_rate, 0, 0, drop_ratio); - ///delete_useless_nodes(ug); delete_useless_nodes(ug); renew_utg(ug, read_g, new_rtg_edges); @@ -22283,7 +22575,7 @@ kvec_asg_arc_t_warp* new_rtg_edges) purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, bubble_dist, - drop_ratio, just_contain); + drop_ratio, just_contain, 0); delete_useless_nodes(ug); renew_utg(ug, read_g, new_rtg_edges); } @@ -22305,12 +22597,19 @@ kvec_asg_arc_t_warp* new_rtg_edges) purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, bubble_dist, - drop_ratio, just_contain); + drop_ratio, just_contain, 0); delete_useless_nodes(ug); renew_utg(ug, read_g, new_rtg_edges); } } + if(asm_opt.purge_level_primary == 0) + { + purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, + asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, bubble_dist, + drop_ratio, 0, 1); + } + n_vtx = read_g->n_seq; @@ -22344,14 +22643,15 @@ kvec_asg_arc_t_warp* new_rtg_edges) } } - if(asm_opt.recover_atg_cov_min == -1 || asm_opt.recover_atg_cov_max == -1) + if(asm_opt.recover_atg_cov_min == -1024 || asm_opt.recover_atg_cov_max == -1024) { asm_opt.recover_atg_cov_max = asm_opt.hom_global_coverage/HOM_PEAK_RATE; asm_opt.recover_atg_cov_min = asm_opt.recover_atg_cov_max * 0.8; asm_opt.recover_atg_cov_max = asm_opt.recover_atg_cov_max * 1.2; } - fprintf(stderr, "asm_opt.recover_atg_cov_min: %d\n", asm_opt.recover_atg_cov_min); - fprintf(stderr, "asm_opt.recover_atg_cov_max: %d\n", asm_opt.recover_atg_cov_max); + + fprintf(stderr, "[M::%s] primary contig coverage range: [%d, %d]\n", + __func__, asm_opt.recover_atg_cov_min, asm_opt.recover_atg_cov_max); recover_utg_by_coverage(ug, read_g, coverage_cut, sources, ruIndex); @@ -22392,12 +22692,12 @@ long long tipsLen, R_to_U* ruIndex, int max_hang, int min_ovlp) char* gfa_name = (char*)malloc(strlen(output_file_name)+35); sprintf(gfa_name, "%s.p_utg.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); - ma_ug_print(ug, &R_INF, coverage_cut, output_file); + ma_ug_print(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); fclose(output_file); sprintf(gfa_name, "%s.p_utg.noseq.gfa", output_file_name); output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); + ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); fclose(output_file); free(gfa_name); @@ -22428,12 +22728,12 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov char* gfa_name = (char*)malloc(strlen(output_file_name)+35); sprintf(gfa_name, "%s.p_ctg.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); - ma_ug_print(ug, &R_INF, coverage_cut, output_file); + ma_ug_print(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "ptg", output_file); fclose(output_file); sprintf(gfa_name, "%s.p_ctg.noseq.gfa", output_file_name); output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); + ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "ptg", output_file); fclose(output_file); free(gfa_name); @@ -22445,7 +22745,7 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov void output_contig_graph_alternative(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, -ma_hit_t_alloc* sources, int max_hang, int min_ovlp) +ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) { kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); @@ -22457,12 +22757,12 @@ ma_hit_t_alloc* sources, int max_hang, int min_ovlp) char* gfa_name = (char*)malloc(strlen(output_file_name)+35); sprintf(gfa_name, "%s.a_ctg.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); - ma_ug_print(ug, &R_INF, coverage_cut, output_file); + ma_ug_print(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "atg", output_file); fclose(output_file); sprintf(gfa_name, "%s.a_ctg.noseq.gfa", output_file_name); output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); + ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "atg", output_file); fclose(output_file); free(gfa_name); @@ -26507,7 +26807,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) { char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); sprintf(buf, "%s.dip", output_file_name); - output_unitig_graph(sg, coverage_cut, buf, sources, max_hang_length, mini_overlap_length); + output_unitig_graph(sg, coverage_cut, buf, sources, ruIndex, max_hang_length, mini_overlap_length); free(buf); output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, @@ -26520,7 +26820,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) else { - output_unitig_graph(sg, coverage_cut, output_file_name, sources, max_hang_length, mini_overlap_length); + output_unitig_graph(sg, coverage_cut, output_file_name, sources, ruIndex, max_hang_length, mini_overlap_length); if(VERBOSE >= 1) { @@ -26534,7 +26834,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) bubble_dist, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length); - output_contig_graph_alternative(sg, coverage_cut, output_file_name, sources, max_hang_length, mini_overlap_length); + output_contig_graph_alternative(sg, coverage_cut, output_file_name, sources, ruIndex, max_hang_length, mini_overlap_length); } *coverage_cut_ptr = coverage_cut; diff --git a/Overlaps.h b/Overlaps.h index d8a1b4b..5340ab9 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1047,8 +1047,8 @@ R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t is_ uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges); void deduplicate(ma_ug_t *src, asg_t *read_g, ma_hit_t_alloc* reverse_sources, long long minLongUntig, long long maxShortUntig, float l_untig_rate, float max_node_threshold, R_to_U* ruIndex, uint32_t resolve_tangle); -void all_to_all_deduplicate(ma_ug_t* ug, uint8_t postive_flag, float drop_rate, -ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex); +void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, +ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex); void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex); void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges); diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index 8df0701..f840f1f 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -4136,9 +4136,9 @@ uint32_t minLen, double purge_threshold) get_contig_length(ug, read_g, &primary_bases, &alter_bases); total_bases = primary_bases + alter_bases; - fprintf(stderr, "primary_bases: %lu\n", primary_bases); - fprintf(stderr, "alter_bases: %lu\n", alter_bases); - fprintf(stderr, "total_bases: %lu\n", total_bases); + // fprintf(stderr, "primary_bases: %lu\n", primary_bases); + // fprintf(stderr, "alter_bases: %lu\n", alter_bases); + // fprintf(stderr, "total_bases: %lu\n", total_bases); for (v = 0; v < all_ovlp->num; v++) @@ -4149,9 +4149,9 @@ uint32_t minLen, double purge_threshold) } } purge_bases = purge_bases/2; - fprintf(stderr, "purge_bases: %lu\n", purge_bases); + ///fprintf(stderr, "purge_bases: %lu\n", purge_bases); alter_bases = alter_bases + purge_bases; - fprintf(stderr, "new alter_bases: %lu\n", alter_bases); + ///fprintf(stderr, "new alter_bases: %lu\n", alter_bases); for (v = 0; v < all_ovlp->num; v++) @@ -4172,7 +4172,7 @@ uint32_t minLen, double purge_threshold) void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, -uint32_t just_contain) +uint32_t just_contain, uint32_t just_coverage) { asg_t *purge_g = NULL; purge_g = asg_init(); @@ -4205,12 +4205,17 @@ uint32_t just_contain) hap_alignment_struct_pip hap_buf; long long k_mer_only, coverage_only; - hap_buf.cov_threshold = get_read_coverage_thres(ug, read_g, ruIndex, position_index, - sources, coverage_cut, read_g->n_seq, COV_COUNT, &k_mer_only, &coverage_only); - ///fprintf(stderr, "cov_threshold: %lld\n", hap_buf.cov_threshold); - - - + if(asm_opt.hom_global_coverage != -1) + { + hap_buf.cov_threshold = asm_opt.hom_global_coverage; + } + else + { + hap_buf.cov_threshold = get_read_coverage_thres(ug, read_g, ruIndex, position_index, + sources, coverage_cut, read_g->n_seq, COV_COUNT, &k_mer_only, &coverage_only); + } + + for (v = 0; v < nsg->n_seq; v++) { uId = v; @@ -4256,9 +4261,9 @@ uint32_t just_contain) hap_buf.cov_threshold = k_mer_only * HOM_PEAK_RATE; } } - asm_opt.hom_global_coverage = hap_buf.cov_threshold; - fprintf(stderr, "cov_threshold: %lld\n", hap_buf.cov_threshold); - + if(asm_opt.hom_global_coverage == -1) asm_opt.hom_global_coverage = hap_buf.cov_threshold; + fprintf(stderr, "[M::%s] purge duplication coverage threshold: %lld\n", __func__, hap_buf.cov_threshold); + if(just_coverage) goto end_coverage; ///kt_for(asm_opt.thread_num, hap_alignment_worker, &hap_buf, nsg->n_seq); kt_for(asm_opt.thread_num, hap_alignment_advance_worker, &hap_buf, nsg->n_seq); @@ -4367,6 +4372,8 @@ uint32_t just_contain) } } + end_coverage: + uint32_t is_Unitig; for (v = 0; v < ruIndex->len; v++) { diff --git a/Purge_Dups.h b/Purge_Dups.h index 4ea05a7..5824858 100644 --- a/Purge_Dups.h +++ b/Purge_Dups.h @@ -15,7 +15,7 @@ void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, -uint32_t just_contain); +uint32_t just_contain, uint32_t just_coverage); void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge, uint32_t is_circle, uint64_t* rLen); void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen); diff --git a/hifiasm.1 b/hifiasm.1 index 9d8622e..036e80a 100644 --- a/hifiasm.1 +++ b/hifiasm.1 @@ -1,4 +1,4 @@ -.TH hifiasm 1 "12 Apr 2020" "hifiasm-0.5.0" "Bioinformatics tools" +.TH hifiasm 1 "27 June 2020" "hifiasm-0.8 (r279)" "Bioinformatics tools" .SH NAME .PP @@ -192,6 +192,12 @@ and do the assembly directly and quickly. This might be helpful when users want to get an optimized assembly by multiple rounds of experiments with different parameters. +.TP +.BI --pri-range \ INT,INT +Min and max coverage cutoff of primary contigs. +Keep contigs with coverage in this range at p_ctg.gfa. +Inferred automatically in default. +Set -1,-1 to disable .SS Trio-partition options @@ -252,6 +258,11 @@ Similarity threshold for duplicate haplotigs that should be purged [0.75]. .BI -O \ FLOAT Min number of overlapped reads for duplicate haplotigs that should be purged [1]. +.TP +.BI --purge-cov \ INT +Coverage upper bound of Purge-dups, which is inferred automatically in default. +If the coverage of a contig is higher than this bound, don't apply Purge-dups. + .SS Debugging options .TP 10