0.16.0 release

This commit is contained in:
chhylp123
2021-08-22 21:58:00 -04:00
parent 37b07e4d33
commit e230b3a55f
10 changed files with 126 additions and 85 deletions
+36 -1
View File
@@ -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 <INT times to rescue repetitive overlaps [%d]\n", asm_opt->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 <INT in size in contig graphs [%lld]\n", asm_opt->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);
}
+2 -1
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#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;
+16 -4
View File
@@ -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);
+2 -2
View File
@@ -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.
+2 -2
View File
@@ -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.
+11
View File
@@ -98,6 +98,17 @@ Error correction options
**\-z <INT=0>**
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 <INT=2000>**
Employ k-mers occurring < ``INT`` times to rescue repetitive overlaps. This option may improve the resolution of repeats.
.. _hg-size-opt:
**\-\-hg-size <INT(k/m/g)>**
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:
+14 -1
View File
@@ -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
+31 -1
View File
@@ -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;
+10 -72
View File
@@ -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);
+2 -1
View File
@@ -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);