diff --git a/CommandLines.cpp b/CommandLines.cpp index 7f92410..04821b5 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -43,6 +43,7 @@ static ko_longopt_t long_options[] = { { "fast", ko_no_argument, 329 }, { "dp-er", ko_required_argument, 330}, { "max-kocc", ko_required_argument, 331}, + { "hg-size", ko_required_argument, 332}, { 0, 0, 0 } }; @@ -70,6 +71,10 @@ void Print_H(hifiasm_opt_t* asm_opt) 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, " -r INT round of correction [%d]\n", asm_opt->number_of_round); fprintf(stderr, " -z INT length of adapters that should be removed [%d]\n", asm_opt->adapterLen); + fprintf(stderr, " --max-kocc INT\n"); + fprintf(stderr, " employ k-mers occurring max_kmer_cnt); + fprintf(stderr, " --hg-size INT(k, m or g)\n"); + fprintf(stderr, " estimated haploid genome size used for inferring read coverage [auto]\n"); fprintf(stderr, " Assembly:\n"); fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round); fprintf(stderr, " -m INT pop bubbles of large_pop_bubble_size); @@ -213,6 +218,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->scffold = 0; asm_opt->dp_min_len = 2000; asm_opt->dp_e = 0.0025; + asm_opt->hg_size = -1; } void destory_enzyme(enzyme* f) @@ -514,6 +520,18 @@ int check_option(hifiasm_opt_t* asm_opt) return 0; } + if(asm_opt->max_kmer_cnt < 0) + { + fprintf(stderr, "[ERROR] [--max-kocc] must >= 0\n"); + return 0; + } + + if(asm_opt->hg_size < -1) + { + fprintf(stderr, "[ERROR] [--hg-size] wrong genome size\n"); + return 0; + } + return 1; } @@ -599,6 +617,23 @@ void get_hic_enzymes(char *argv, enzyme** x, int check_name) (*x)->a[k][(*x)->l[k]] = '\0'; } +int64_t inter_gsize(char *argv) +{ + int64_t len = strlen(argv); + double s; + if(len <= 1) return -2; + char t = argv[len-1]; + if(t != 'k' && t != 'K' && t != 'm' && t != 'M' && t != 'g' && t != 'G') return -2; + char *ss=(char*)malloc(len); + memcpy(ss, argv, len-1); ss[len-1] = '\0'; + s = atof(ss); + free(ss); + if(t == 'k' || t == 'K') return s*1000; + if(t == 'm' || t == 'M') return s*1000000; + if(t == 'g' || t == 'G') return s*1000000000; + return s; +} + int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) { ketopt_t opt = KETOPT_INIT; @@ -686,6 +721,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 329) asm_opt->flag |= HA_F_FAST; else if (c == 330) asm_opt->dp_e = atof(opt.arg); else if (c == 331) asm_opt->max_kmer_cnt = atol(opt.arg); + else if (c == 332) asm_opt->hg_size = inter_gsize(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); @@ -716,6 +752,5 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) get_queries(argc, argv, &opt, asm_opt); - return check_option(asm_opt); } diff --git a/CommandLines.h b/CommandLines.h index db780f7..2f64e03 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.15.5-r366" +#define HA_VERSION "0.16.0-r369" #define VERBOSE 0 @@ -113,6 +113,7 @@ typedef struct { uint64_t scffold; int32_t dp_min_len; float dp_e; + int64_t hg_size; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index e84bf4e..bc15e08 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -122,8 +122,8 @@ void print_peak(long long* cov_buf, long long cov_buf_length, long long max_i) } -void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo_peak_cov, -long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* coverage_only) +void get_read_peak(asg_t *read_g, long long* cov_buf, long long cov_buf_length, long long* topo_peak_cov, +long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* coverage_only, long long g_size) { long long i, start, err_i, max_i, max2_i, max3_i, topo_peak_i, max, max2, max3, topo_peak, min; @@ -242,6 +242,18 @@ long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* cove coverage_hom = max_i; } + if(g_size > 0) { + long long n_bs, m_peak_hom = -1; + int p_ht = -1; + for (i = n_bs = 0; i < read_g->n_seq; i++) n_bs += read_g->seq[i].len; + m_peak_hom = n_bs/g_size; + if(m_peak_hom > 0) { + p_ht = -1; + coverage_hom = adj_m_peak_hom(m_peak_hom, max_i, max2_i, max3_i, &p_ht); + coverage_het = p_ht; + } + } + if(k_mer_het != -1) { (*het_peak) = k_mer_het; @@ -410,8 +422,8 @@ long long* k_mer_only, long long* coverage_only) cov_buf_length); } - get_read_peak(cov_buf, cov_buf_length, alter_peak == -1? NULL: &alter_peak, &hom_peak, &het_peak, - k_mer_only, coverage_only); + get_read_peak(read_g, cov_buf, cov_buf_length, alter_peak == -1? NULL: &alter_peak, &hom_peak, &het_peak, + k_mer_only, coverage_only, asm_opt.hg_size); free(cov_buf); diff --git a/docs/source/conf.py b/docs/source/conf.py index 4a4f792..21cb778 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -38,9 +38,9 @@ copyright = u'2021, Haoyu Cheng, Heng Li' # built documents. # # The short X.Y version. -version = '0.15.5-r350' +version = '0.16.0-r369' # The full version, including alpha/beta/rc tags. -release = '0.15.5' +release = '0.16.0' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/docs/source/faq.rst b/docs/source/faq.rst index d27497d..c83e5f5 100644 --- a/docs/source/faq.rst +++ b/docs/source/faq.rst @@ -44,7 +44,7 @@ Why one Hi-C integrated assembly is larger than another one? Another possibility is that hifiasm misidentifies coverage threshold for homozygous reads. For instance, hifiasm prints the following information during assembly: :: - [M::purge_dups] purge duplication coverage threshold: 36 + [M::purge_dups] homozygous read coverage threshold: 36 In this example, hifiasm identifies the coverage threshold for homozygous reads as ``36``. If it is significantly smaller than the homozygous coverage peak, hifiasm will generate two unbalanced assemblies. In this case, please set ``--hom-cov`` to homozygous coverage peak. Please note that tuning ``--hom-cov`` may affect ``*p_utg*gfa`` so that ``*hic*.bin`` should be deleted. Since v0.15.5, hifiasm can detect such changes and renew Hi-C bin files automatically. @@ -60,7 +60,7 @@ For Hi-C integrated assembly, why the assembly size of both haplotypes are much If most bases of a diploid sample are homozygous, the coverage threshold is wrongly determined by hifiasm. For instance, hifiasm prints the following information during assembly: :: - [M::purge_dups] purge duplication coverage threshold: 36 + [M::purge_dups] homozygous read coverage threshold: 36 In this example, hifiasm identifies the coverage threshold for homozygous reads as ``36``. If it is much smaller than homozygous coverage peak, hifiasm thinks most reads are homozygous and assign them to both assemblies, making both of them much larger than the estimated haploid genome size. In this case, please set ``--hom-cov`` to homozygous coverage peak. Please note that tuning ``--hom-cov`` may affect ``*p_utg*gfa`` so that ``*hic*.bin`` should be deleted. Since v0.15.5, hifiasm can detect such changes and renew Hi-C bin files automatically. diff --git a/docs/source/parameter-reference.rst b/docs/source/parameter-reference.rst index e6b0b96..cac0ed9 100644 --- a/docs/source/parameter-reference.rst +++ b/docs/source/parameter-reference.rst @@ -98,6 +98,17 @@ Error correction options **\-z ** Length of adapters that should be removed. This option remove ``INT`` bases from both ends of each read. Some old HiFi reads may consist of short adapters (e.g. 20bp adapter at one end). For such data, trimming short adapters would significantly improve the assembly quality. +.. _max-kocc-opt: + +**\-\-max-kocc ** + Employ k-mers occurring < ``INT`` times to rescue repetitive overlaps. This option may improve the resolution of repeats. + + +.. _hg-size-opt: + +**\-\-hg-size ** + Estimated haploid genome size used for inferring read coverage. This option is used to get accurate homozygous read coverage during error correction. Common suffices are required, for example, 100m or 3g. + .. _min-hist-cnt-opt: diff --git a/hifiasm.1 b/hifiasm.1 index ecda27e..6f19a24 100644 --- a/hifiasm.1 +++ b/hifiasm.1 @@ -1,4 +1,4 @@ -.TH hifiasm 1 "25 July 2021" "hifiasm-0.15.5 (r350)" "Bioinformatics tools" +.TH hifiasm 1 "22 August 2021" "hifiasm-0.16.0 (r369)" "Bioinformatics tools" .SH NAME .PP @@ -161,6 +161,19 @@ Some old Hifi reads may consist of short adapters (e.g., 20bp adapter at one end). For such data, trimming short adapters would significantly improve the assembly quality. +.TP +.BI --max-kocc \ INT +Employ k-mers occurring < +.IR INT +times to rescue repetitive overlaps [2000]. +This option may improve the resolution of repeats. + +.TP +.BI --hg-size \ INT (k/m/g) +Estimated haploid genome size used for inferring read coverage [auto]. +This option is used to get accurate homozygous read coverage during +error correction. Common suffices are required, for example, 100m or 3g. + .TP .BI --min-hist-cnt \ INT When analyzing the k-mer spectrum, ignore counts below diff --git a/hist.cpp b/hist.cpp index 98d9220..88ee752 100644 --- a/hist.cpp +++ b/hist.cpp @@ -43,7 +43,35 @@ void print_hist_lines(int n_cnt, int start_cnt, const int64_t *cnt) } } -int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het) +int adj_m_peak_hom(int m_peak_hom, int max_i, int max2_i, int max3_i, int *peak_het) +{ + int64_t mm[3], d, min_i, min_d, i; + mm[0] = max2_i; mm[1] = max_i; mm[2] = max3_i; + for (i = 0, min_i = -1, min_d = -1; i < 3; i++){ + if(mm[i] <= 0) continue; + d = (mm[i] >= m_peak_hom?mm[i]-m_peak_hom:m_peak_hom-mm[i]); + if(min_d == -1 || min_d > d || (min_d == d && i == 1)){ + min_d = d; min_i = i; + } + } + if(min_i < 0) return m_peak_hom; + if(mm[min_i] < m_peak_hom){ + d = m_peak_hom - mm[min_i]; + if(d >= mm[min_i]*0.51) { + *peak_het = mm[min_i]; + return m_peak_hom; + } + } + + for (i = min_i-1; i >= 0; i--){ + if(mm[i] <= 0) continue; + *peak_het = mm[i]; + break; + } + return mm[min_i]; +} + +int ha_analyze_count(int n_cnt, int start_cnt, int m_peak_hom, const int64_t *cnt, int *peak_het) { const int hist_max = 100; int i, start, low_i, max_i, max2_i, max3_i; @@ -117,6 +145,8 @@ int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het } if (max3 > 0) fprintf(stderr, "[M::%s] right: count[%d] = %ld\n", __func__, max3_i, (long)cnt[max3_i]); else fprintf(stderr, "[M::%s] right: none\n", __func__); + + if(m_peak_hom > 0) return adj_m_peak_hom(m_peak_hom, max_i, max2_i, max3_i, peak_het); if (max3_i > 0) { *peak_het = max_i; return max3_i; diff --git a/htab.cpp b/htab.cpp index 04c02e4..0a01a65 100644 --- a/htab.cpp +++ b/htab.cpp @@ -130,6 +130,7 @@ typedef struct { typedef struct { int k, pre, n_hash, n_shift; uint64_t tot; ///number of distinct k-mers + uint64_t bs; ha_ct1_t *h; } ha_ct_t; @@ -851,7 +852,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip if (p->pt && !(p->flag&HAF_COUNT_REFINE)) free(s->buf[i].b); else free(s->buf[i].a); } - if (p->ct) p->ct->tot += n_ins; + if (p->ct) p->ct->tot += n_ins, p->ct->bs += s->sum_len; if (p->pt) p->pt->tot_pos += n_ins; free(s->buf); #if 0 @@ -945,6 +946,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt ///building a large hash table consisting of 4096 small hash tables pl.ct = ha_ct_init(opt->k, opt->pre, opt->bf_n_hash, opt->bf_shift); } + if(pl.ct) pl.ct->bs = 0; kt_pipeline(3, worker_count, &pl, 3); if (read_rs) { destory_UC_Read(&pl.ucr); @@ -961,6 +963,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi { int i; int64_t n_seq = 0; + uint64_t n_bs = 0; yak_copt_t opt; ha_ct_t *h = 0; assert(!(flag & HAF_RS_WRITE_LEN) || !(flag & HAF_RS_WRITE_SEQ)); // not both @@ -984,8 +987,11 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi opt.adaLen = (keep_adapter? asm_opt->adapterLen : 0); opt.min_rcnt = (low_freq?*low_freq:-1); ///asm_opt->num_reads is the number of fastq files - for (i = 0; i < asm_opt->num_reads; ++i) + for (i = n_bs = 0; i < asm_opt->num_reads; ++i){ h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq); + if(h) n_bs += h->bs; + } + if(h) h->bs = n_bs; if (h && opt.bf_shift > 0) ha_ct_destroy_bf(h); return h; @@ -1140,7 +1146,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i if(!(ex_flag & HAF_SKIP_READ)) { ha_ct_hist(h, cnt, asm_opt->thread_num); - peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, cnt, &peak_het); + peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, asm_opt->hg_size>0?(h->bs/asm_opt->hg_size):(-1), cnt, &peak_het); if (hom_cov) *hom_cov = peak_hom; if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); ///in default, asm_opt->high_factor = 5.0 @@ -1155,38 +1161,6 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i return (void*)flt_tab; } -void *ha_ft_gen_worse(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode) -{ - yak_ft_t *flt_tab; - int64_t cnt[YAK_N_COUNTS]; - int peak_hom, peak_het, cutoff = YAK_MAX_COUNT - 1, ex_flag = 0; - if(is_hp_mode) ex_flag = HAF_RS_READ|HAF_SKIP_READ; - ha_ct_t *h; - h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs, NULL, 1, NULL); - if((asm_opt->flag & HA_F_VERBOSE_GFA)) - { - write_ct_index((void*)h, asm_opt->output_file_name); - // load_ct_index(&ha_ct_table, asm_opt->output_file_name); - // debug_ct_index((void*)h, ha_ct_table); - // debug_ct_index(ha_ct_table, (void*)h); - // ha_ct_destroy((ha_ct_t *)ha_ct_table); - } - - if(!(ex_flag & HAF_SKIP_READ)) - { - ha_ct_hist(h, cnt, asm_opt->thread_num); - peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, cnt, &peak_het); - if (hom_cov) *hom_cov = peak_hom; - if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); - } - ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num); - flt_tab = gen_hh(h, 1/**asm_opt->max_kmer_cnt**/); - ha_ct_destroy(h); - fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> filtered out %ld k-mers occurring %d or more times\n", __func__, - yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb(), (long)kh_size(flt_tab), cutoff); - return (void*)flt_tab; -} - ha_pt_t *ha_pt_gen_dp(const hifiasm_opt_t *asm_opt, ha_ct_t *ct, int flag, int n_thread, const void *flt_tab, All_reads *rs, int peak_hom, int peak_het) { int low_freq = mz_low_b(peak_hom, peak_het); @@ -1203,42 +1177,6 @@ ha_pt_t *ha_pt_gen_dp(const hifiasm_opt_t *asm_opt, ha_ct_t *ct, int flag, int n return pt; } -ha_pt_t *ha_pt_gen_worst(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov) -{ - int64_t cnt[YAK_N_COUNTS]; - int peak_hom, peak_het, extra_flag1, extra_flag2; - ha_ct_t *ct; - ha_pt_t *pt; - if (read_from_store) {///if reads have already been read - extra_flag1 = extra_flag2 = HAF_RS_READ; - } else if (rs->total_reads == 0) {///if reads & length have not been scanned - extra_flag1 = HAF_RS_WRITE_LEN; - extra_flag2 = HAF_RS_WRITE_SEQ; - } else {///if length has been loaded but reads have not - extra_flag1 = HAF_RS_WRITE_SEQ; - extra_flag2 = HAF_RS_READ; - } - if(is_hp_mode) extra_flag1 |= HAF_SKIP_READ, extra_flag2 |= HAF_SKIP_READ; - - ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs, NULL, 1, NULL); - fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__, - yak_realtime(), yak_cpu_usage(), (long)ct->tot); - ha_ct_hist(ct, cnt, asm_opt->thread_num); - fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]); - peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, cnt, &peak_het); - if (hom_cov) *hom_cov = peak_hom; - if (het_cov) *het_cov = peak_het; - if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); - - ///here ha_ct_shrink is mostly used to remove k-mer appearing only 1 time - ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num); - pt = ha_pt_gen_dp(asm_opt, ct, HAF_COUNT_EXACT|extra_flag2, asm_opt->thread_num, flt_tab, rs, peak_hom, peak_het); - //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); - return pt; -} - ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov) { int64_t cnt[YAK_N_COUNTS], tot_cnt; @@ -1261,7 +1199,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f yak_realtime(), yak_cpu_usage(), (long)ct->tot); ha_ct_hist(ct, cnt, asm_opt->thread_num); fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]); - peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, cnt, &peak_het); + peak_hom = ha_analyze_count(YAK_N_COUNTS, asm_opt->min_hist_kmer_cnt, asm_opt->hg_size>0?(ct->bs/asm_opt->hg_size):(-1), cnt, &peak_het); if (hom_cov) *hom_cov = peak_hom; if (het_cov) *het_cov = peak_het; if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); diff --git a/htab.h b/htab.h index 9e2caa0..cd46dea 100644 --- a/htab.h +++ b/htab.h @@ -102,7 +102,8 @@ double yak_cpu_usage(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 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); -int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het); +int ha_analyze_count(int n_cnt, int start_cnt, int m_peak_hom, const int64_t *cnt, int *peak_het); +int adj_m_peak_hom(int m_peak_hom, int max_i, int max2_i, int max3_i, int *peak_het); void print_hist_lines(int n_cnt, int start_cnt, const int64_t *cnt); void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs);