diff --git a/CommandLines.cpp b/CommandLines.cpp index 739e44a..1da7d60 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -21,7 +21,7 @@ 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 }, + { "hom-cov", ko_required_argument, 309 }, { "pri-range", ko_required_argument, 310 }, { "lowQ", ko_required_argument, 312 }, { "min-hist-cnt", ko_required_argument, 313 }, @@ -39,6 +39,7 @@ static ko_longopt_t long_options[] = { { "n-hap", ko_required_argument, 325 }, { "n-weight", ko_required_argument, 326 }, { "l-msjoin", ko_required_argument, 327 }, + { "purge-max", ko_required_argument, 328 }, { 0, 0, 0 } }; @@ -55,9 +56,8 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, "Options:\n"); fprintf(stderr, " Input/Output:\n"); fprintf(stderr, " -o STR prefix of output files [%s]\n", asm_opt->output_file_name); - fprintf(stderr, " -i ignore saved read correction and overlaps\n"); fprintf(stderr, " -t INT number of threads [%d]\n", asm_opt->thread_num); - fprintf(stderr, " -z INT length of adapters that should be removed [%d]\n", asm_opt->adapterLen); + fprintf(stderr, " -h show help information\n"); fprintf(stderr, " --version show version number\n"); fprintf(stderr, " Overlap/Error correction:\n"); fprintf(stderr, " -k INT k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length); @@ -68,12 +68,16 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " -r INT round of correction [%d]\n", asm_opt->number_of_round); fprintf(stderr, " Assembly:\n"); fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round); + fprintf(stderr, " -z INT length of adapters that should be removed [%d]\n", asm_opt->adapterLen); fprintf(stderr, " -m INT pop bubbles of large_pop_bubble_size); fprintf(stderr, " -p INT pop bubbles of small_pop_bubble_size); fprintf(stderr, " -n INT remove tip unitigs composed of <=INT reads [%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, " -i ignore saved read correction and overlaps\n"); fprintf(stderr, " -u disable post join contigs step which may improve N50\n"); + fprintf(stderr, " --hom-cov INT\n"); + fprintf(stderr, " homozygous read coverage [auto]\n"); fprintf(stderr, " --lowQ INT\n"); fprintf(stderr, " output contig regions with >=INT%% inconsistency in BED format; 0 to disable [%d]\n", asm_opt->bed_inconsist_rate); fprintf(stderr, " --b-cov INT\n"); @@ -105,7 +109,7 @@ void Print_H(hifiasm_opt_t* asm_opt) asm_opt->purge_simi_rate_l2, asm_opt->purge_simi_rate_l3); fprintf(stderr, " -O INT min number of overlapped reads for duplicate haplotigs [%d]\n", asm_opt->purge_overlap_len); - fprintf(stderr, " --purge-cov INT\n"); + fprintf(stderr, " --purge-max INT\n"); fprintf(stderr, " coverage upper bound of Purge-dups [auto]\n"); fprintf(stderr, " --n-hap INT\n"); fprintf(stderr, " number of haplotypes [%d]\n", asm_opt->polyploidy); @@ -183,6 +187,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->recover_atg_cov_max = INT_MAX; asm_opt->hom_global_coverage = -1; asm_opt->hom_global_coverage_set = 0; + asm_opt->pur_global_coverage = -1; asm_opt->bed_inconsist_rate = 70; asm_opt->hic_inconsist_rate = 30; ///asm_opt->bub_mer_length = 3; @@ -410,7 +415,13 @@ int check_option(hifiasm_opt_t* asm_opt) 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"); + fprintf(stderr, "[ERROR] homozygous read coverage should be >= 0 (--hom-cov)\n"); + return 0; + } + + if(asm_opt->pur_global_coverage < 0 && asm_opt->pur_global_coverage != -1) + { + fprintf(stderr, "[ERROR] purge duplication coverage threshold should be >= 0 (--purge-max)\n"); return 0; } @@ -663,6 +674,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 325) asm_opt->polyploidy = atoi(opt.arg); else if (c == 326) asm_opt->n_weight = atoi(opt.arg); else if (c == 327) asm_opt->misjoin_len = atol(opt.arg); + else if (c == 328) asm_opt->pur_global_coverage = atoi(opt.arg); 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); diff --git a/CommandLines.h b/CommandLines.h index 0da1f99..6c285c2 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.15.4-r348" +#define HA_VERSION "0.15.4-r349" #define VERBOSE 0 @@ -77,6 +77,7 @@ typedef struct { int recover_atg_cov_max; int hom_global_coverage; int hom_global_coverage_set; + int pur_global_coverage; int bed_inconsist_rate; int hic_inconsist_rate; diff --git a/Overlaps.cpp b/Overlaps.cpp index c224a11..560510e 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -16995,7 +16995,8 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) if(asm_opt.recover_atg_cov_min == -1024) { - asm_opt.recover_atg_cov_max = asm_opt.hom_global_coverage/HOM_PEAK_RATE; + asm_opt.recover_atg_cov_max = (asm_opt.hom_global_coverage_set? + (asm_opt.hom_global_coverage):(((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE))); asm_opt.recover_atg_cov_min = asm_opt.recover_atg_cov_max * 0.85; asm_opt.recover_atg_cov_max = INT32_MAX; } @@ -24806,7 +24807,8 @@ uint32_t collect_p_trans, uint32_t collect_p_trans_f) if(asm_opt.recover_atg_cov_min == -1024) { - asm_opt.recover_atg_cov_max = asm_opt.hom_global_coverage/HOM_PEAK_RATE; + asm_opt.recover_atg_cov_max = (asm_opt.hom_global_coverage_set? + (asm_opt.hom_global_coverage):(((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE))); asm_opt.recover_atg_cov_min = asm_opt.recover_atg_cov_max * 0.85; asm_opt.recover_atg_cov_max = INT32_MAX; } diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index ba6d81d..269cc9c 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -5267,9 +5267,14 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle hap_alignment_struct_pip hap_buf; long long k_mer_only, coverage_only; - if(asm_opt.hom_global_coverage != -1) + if(asm_opt.pur_global_coverage != -1) { - hap_buf.cov_threshold = asm_opt.hom_global_coverage; + hap_buf.cov_threshold = asm_opt.pur_global_coverage; + } + else if(asm_opt.hom_global_coverage != -1) + { + hap_buf.cov_threshold = (asm_opt.hom_global_coverage_set? + (((double)asm_opt.hom_global_coverage)*((double)HOM_PEAK_RATE)):(asm_opt.hom_global_coverage)); } else { diff --git a/hic.cpp b/hic.cpp index 242a9a1..ab44a93 100644 --- a/hic.cpp +++ b/hic.cpp @@ -43,6 +43,7 @@ KRADIX_SORT_INIT(u_trans_m, u_trans_t, u_trans_m_key, 8) KRADIX_SORT_INIT(u_trans_occ, u_trans_t, u_trans_occ_key, member_size(u_trans_t, occ)) #define is_hom_hit(a) ((a).id == (uint64_t)-1) +#define HC_PT_MA 65 typedef struct { kv_gg_status sg; @@ -555,7 +556,9 @@ int write_hc_pt_index(ha_ug_index* idx, char* file_name) free(gfa_name); return 0; } + uint64_t i = HC_PT_MA; + fwrite(&i, sizeof(i), 1, fp); fwrite(&idx->uID_bits, sizeof(idx->uID_bits), 1, fp); fwrite(&idx->uID_mode, sizeof(idx->uID_mode), 1, fp); fwrite(&idx->pos_bits, sizeof(idx->pos_bits), 1, fp); @@ -566,7 +569,6 @@ int write_hc_pt_index(ha_ug_index* idx, char* file_name) fwrite(&idx->tot, sizeof(idx->tot), 1, fp); fwrite(&idx->tot_pos, sizeof(idx->tot_pos), 1, fp); - uint64_t i = 0; for (i = 0; i < idx->tot; i++) { fwrite(&idx->idx_buf[i].n, sizeof(idx->idx_buf[i].n), 1, fp); @@ -596,7 +598,19 @@ int load_hc_pt_index(ha_ug_index** r_idx, ma_ug_t *ug, char* file_name) return 0; } ha_ug_index* idx = NULL; CALLOC(idx, 1); + uint64_t i; + flag += fread(&i, sizeof(i), 1, fp); + if(i != HC_PT_MA) + { + free(gfa_name); + destory_hc_pt_index(idx); + free(idx); + (*r_idx) = NULL; + fclose(fp); + fprintf(stderr, "[M::%s::] ==> Renew Hi-C index\n", __func__); + return 0; + } flag += fread(&idx->uID_bits, sizeof(idx->uID_bits), 1, fp); flag += fread(&idx->uID_mode, sizeof(idx->uID_mode), 1, fp); flag += fread(&idx->pos_bits, sizeof(idx->pos_bits), 1, fp); @@ -607,7 +621,7 @@ int load_hc_pt_index(ha_ug_index** r_idx, ma_ug_t *ug, char* file_name) flag += fread(&idx->tot, sizeof(idx->tot), 1, fp); flag += fread(&idx->tot_pos, sizeof(idx->tot_pos), 1, fp); MALLOC(idx->idx_buf, idx->tot); - uint64_t i = 0; + for (i = 0; i < idx->tot; i++) { flag += fread(&idx->idx_buf[i].n, sizeof(idx->idx_buf[i].n), 1, fp); @@ -876,6 +890,14 @@ ha_ug_index* build_unitig_index(ma_ug_t *ug, int k, uint64_t up_occ, uint64_t lo fprintf(stderr, "[M::%s::%.3f] ==> Sorting pos\n", __func__, yak_realtime()-beg_time); fprintf(stderr, "[M::%s::%.3f] ==> HiC index has been built\n", __func__, yak_realtime()-index_time); + + uint64_t i; + for (i = 0; i < idx->tot; i++) + { + kv_destroy(pl.cnt[i].a); + kv_destroy(pl.buf[i].a); + } + free(pl.cnt); free(pl.buf); return idx; } @@ -4043,6 +4065,7 @@ void update_containment_distance(asg_t *sg, kv_u_trans_t *ta, hc_links* link) } kv_destroy(buf.a); + free(uc_idx); } @@ -4877,7 +4900,12 @@ int load_hc_hits(kvec_pe_hit* hits, ma_ug_t* ug, const char *fn) FILE* fp = NULL; fp = fopen(buf, "r"); - if(!fp) return 0; + if(!fp) + { + free(buf); + return 0; + } + kv_init(hits->a); flag += fread(&hits->a.n, sizeof(hits->a.n), 1, fp); @@ -9506,11 +9534,20 @@ H_partition* hap, int8_t *s, mc_gg_status *sa, trans_idx* dis) if((buf.a[k]&1) == 1) f_idx = k; } buf.n = MIN(r_idx, f_idx); - + trans_p_t* p = NULL; dis->n = 0; uint64_t bin_size = MIN(2250, buf.n>>8), m; + if(bin_size == 0) + { + for (i = 8; i > 0; i--) + { + bin_size = buf.n>>i; + if(bin_size > 0) break; + } + if(bin_size == 0) bin_size = buf.n; + } i = 0; while (i < buf.n) { @@ -9603,7 +9640,7 @@ H_partition* hap, int8_t *s, mc_gg_status *sa, trans_idx* dis) dis->a[i].end += ((dis->a[i+1].beg - dis->a[i].end)/2); dis->a[i+1].beg = dis->a[i].end; } - + // for (i = 0; i < dis->n; i++) // { // if(i > 0 && dis->a[i].beg != dis->a[i-1].end) fprintf(stderr, "ERROR: dis->a[i].beg: %lu, dis->a[i-1].end: %lu\n", dis->a[i].beg, dis->a[i-1].end); @@ -9617,7 +9654,7 @@ H_partition* hap, int8_t *s, mc_gg_status *sa, trans_idx* dis) dis->max = dis->a[dis->n-1].end; - + if(idx->a < 0) idx->a = 0; if(idx->a == 0) @@ -16166,7 +16203,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o idx->hap_cnt = asm_opt.hap_occ; kv_init(sl.hits.a); kv_init(sl.hits.idx); kv_init(sl.hits.occ); - + if(!load_hc_hits(&sl.hits, idx->ug, asm_opt.output_file_name)) { alignment_worker_pipeline(&sl, fn1, fn2); @@ -16448,7 +16485,7 @@ void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, else hic_short_align_poy(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt); - destory_hc_pt_index(ug_index); + destory_hc_pt_index(ug_index);free(ug_index); } diff --git a/khashl.h b/khashl.h index 4ecf294..6fddd39 100644 --- a/khashl.h +++ b/khashl.h @@ -136,25 +136,33 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26 #define __KHASHL_IMPL_S_L(SCOPE, HType, prefix, khkey_t) \ SCOPE khint_t prefix##_save(HType *h, FILE* fp) { \ if (!h) return 0; \ + uint8_t ff; \ khint_t n_buckets = (h->keys? 1U<bits : 0U); \ fwrite(&n_buckets, sizeof(n_buckets), 1, fp); \ fwrite(&h->bits, sizeof(h->bits), 1, fp); \ fwrite(&h->count, sizeof(h->count), 1, fp); \ - fwrite(h->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \ - fwrite(h->keys, sizeof(khkey_t), n_buckets, fp); \ + ff = h->used? 1:0; fwrite(&ff, sizeof(ff), 1, fp); \ + if(ff) fwrite(h->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \ + ff = h->keys? 1:0; fwrite(&ff, sizeof(ff), 1, fp); \ + if(ff) fwrite(h->keys, sizeof(khkey_t), n_buckets, fp); \ return 1; \ } \ SCOPE khint_t prefix##_load(HType **h, FILE* fp) { \ (*h) = prefix##_init(); \ khint_t n_buckets; \ uint64_t flag = 0;\ + uint8_t ff; \ flag += fread(&n_buckets, sizeof(n_buckets), 1, fp); \ flag += fread(&(*h)->bits, sizeof((*h)->bits), 1, fp); \ flag += fread(&(*h)->count, sizeof((*h)->count), 1, fp); \ + flag += fread(&ff, sizeof(ff), 1, fp); \ + if(ff) {\ (*h)->used = (khint32_t*)kmalloc(__kh_fsize(n_buckets) * sizeof(khint32_t)); \ + flag += fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); }\ + flag += fread(&ff, sizeof(ff), 1, fp); \ + if(ff) {\ (*h)->keys = (khkey_t*)kmalloc(n_buckets * sizeof(khkey_t)); \ - flag += fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \ - flag += fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); \ + flag += fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); }\ return 1; \ } \