xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-09-29 05:08:13 +08:00
Code Issues Packages Projects Releases Wiki Activity

r201: automatically adjust -N

Browse Source
This commit is contained in:
Heng Li
2020-04-05 22:43:01 -04:00
parent 3773610532
commit adbc65feda
5 changed files with 43 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+15 -11
View File
@@ -382,7 +382,7 @@ typedef struct {
overlap_region_alloc olist;
ha_abuf_t *ab;
// error correction related buffers
int64_t num_read_base, num_correct_base, num_recorrect_base, mem_buf;
int64_t num_read_base, num_correct_base, num_recorrect_base;
Cigar_record cigar1;
Graph POA_Graph;
Graph DAGCon;
@@ -626,7 +626,7 @@ void Output_corrected_reads()
void ha_overlap_and_correct(int round)
{
int i;
int i, hom_cov;
ha_ovec_buf_t **b;
ha_ecsave_buf_t *e;
@@ -634,7 +634,9 @@ void ha_overlap_and_correct(int round)
CALLOC(b, asm_opt.thread_num);
for (i = 0; i < asm_opt.thread_num; ++i)
b[i] = ha_ovec_init(0);
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF); // build the index
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF, &hom_cov); // build the index
if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated
ha_opt_update_cov(&asm_opt, hom_cov);
if (asm_opt.required_read_name)
kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads);
else
@@ -647,7 +649,7 @@ void ha_overlap_and_correct(int round)
asm_opt.num_bases += b[i]->num_read_base;
asm_opt.num_corrected_bases += b[i]->num_correct_base;
asm_opt.num_recorrected_bases += b[i]->num_recorrect_base;
asm_opt.mem_buf += b[i]->mem_buf;
asm_opt.mem_buf += ha_ovec_mem(b[i]);
ha_ovec_destroy(b[i]);
}
free(b);
@@ -1135,12 +1137,12 @@ long long readNum, long long rescue_threshold, float cluster_threshold)
void ha_overlap_final(void)
{
int i;
int i, hom_cov;
ha_ovec_buf_t **b;
CALLOC(b, asm_opt.thread_num);
for (i = 0; i < asm_opt.thread_num; ++i)
b[i] = ha_ovec_init(1);
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF); // build the index
ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov); // build the index
kt_for(asm_opt.thread_num, worker_ov_final, b, R_INF.total_reads);
ha_pt_destroy(ha_idx);
ha_idx = 0;
@@ -1151,19 +1153,21 @@ void ha_overlap_final(void)
int ha_assemble(void)
{
int r, ovlp_loaded = 0;
int r, hom_cov = -1, ovlp_loaded = 0;
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 (!ovlp_loaded) {
// construct hash table for high occurrence k-mers
if (!asm_opt.no_kmer_flt)
ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF);
if (!asm_opt.no_kmer_flt) {
ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov);
ha_opt_update_cov(&asm_opt, hom_cov);
}
// error correction
assert(asm_opt.number_of_round > 0);
for (r = 0; r < asm_opt.number_of_round; ++r) {
clear_opt(&asm_opt, r); // this update asm_opt.roundID and a few other fields
ha_opt_reset_to_round(&asm_opt, r); // this update asm_opt.roundID and a few other fields
ha_overlap_and_correct(r);
fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(),
yak_cpu_usage(), yak_peakrss_in_gb(), r + 1);
@@ -1174,7 +1178,7 @@ int ha_assemble(void)
//Output_corrected_reads();
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> written corrected reads to disk\n", __func__, yak_realtime(), yak_cpu_usage());
// overlap between corrected reads
clear_opt(&asm_opt, asm_opt.number_of_round);
ha_opt_reset_to_round(&asm_opt, asm_opt.number_of_round);
ha_overlap_final();
fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(),
yak_cpu_usage(), yak_peakrss_in_gb());
CommandLines.cpp
+17 -14
View File
@@ -35,15 +35,11 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " -w INT minimizer window size [%d]\n", asm_opt->mz_win);
fprintf(stderr, " -f INT number of bits for bloom filter [%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 INT overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
///fprintf(stderr, " -l load all overlaps from disk, can avoid overlap calculation [%d]\n", asm_opt->load_index_from_disk);
///fprintf(stderr, " -i ignore saved overlaps in *.ovlp*.bin files\n");
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, " -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);
@@ -51,8 +47,8 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, " -h show help information\n");
fprintf(stderr, " Trio-partition:\n");
fprintf(stderr, " -P FILE paternal trio index generated by \"yak count\" [NULL]\n");
fprintf(stderr, " -M FILE maternal trio index generated by \"yak count\" [NULL]\n");
fprintf(stderr, " -P FILE paternal trio index generated by \"yak count\" []\n");
fprintf(stderr, " -M FILE maternal trio index generated by \"yak count\" []\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);
@@ -73,10 +69,10 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->k_mer_length = 51;
asm_opt->mz_win = 51;
asm_opt->bf_shift = 37;
asm_opt->high_factor = 5.0f;
asm_opt->no_HPC = 0;
asm_opt->no_kmer_flt = 0;
asm_opt->max_n_chain = 400;
asm_opt->high_factor = 5.0f;
asm_opt->max_n_chain = 100;
asm_opt->k_mer_min_freq = 3;
asm_opt->k_mer_max_freq = 66;
asm_opt->load_index_from_disk = 1;
@@ -107,7 +103,7 @@ void destory_opt(hifiasm_opt_t* asm_opt)
}
}
void clear_opt(hifiasm_opt_t* asm_opt, int round)
void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round)
{
asm_opt->num_bases = 0;
asm_opt->num_corrected_bases = 0;
@@ -116,6 +112,14 @@ void clear_opt(hifiasm_opt_t* asm_opt, int round)
asm_opt->roundID = round;
}
void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov)
{
int max_n_chain = (int)(hom_cov * opt->high_factor + .499);
if (opt->max_n_chain < max_n_chain)
opt->max_n_chain = max_n_chain;
fprintf(stderr, "[M::%s] updated max_n_chain to %d\n", __func__, opt->max_n_chain);
}
int check_file(char* name, const char* opt)
{
if(!name)
@@ -308,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:lw:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:", 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:", long_options)) >= 0) {
if (c == 'h')
{
Print_H(asm_opt);
@@ -325,7 +329,6 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 'r') asm_opt->number_of_round = atoi(opt.arg);
else if (c == 'k') asm_opt->k_mer_length = atoi(opt.arg);
else if (c == 'i') asm_opt->load_index_from_disk = 0;
else if (c == 'l') asm_opt->load_index_from_disk = 1;
else if (c == 'w') asm_opt->mz_win = atoi(opt.arg);
else if (c == 'D') asm_opt->high_factor = atof(opt.arg);
else if (c == 'F') asm_opt->no_kmer_flt = 1;
CommandLines.h
+5 -4
View File
@@ -3,7 +3,7 @@
#include <pthread.h>
#define HA_VERSION "0.3.0-dirty-r200"
#define HA_VERSION "0.3.0-dirty-r201"
#define VERBOSE 0
@@ -18,10 +18,10 @@ typedef struct {
int k_mer_length;
int mz_win;
int bf_shift;
float high_factor;
int no_HPC;
int no_kmer_flt;
int max_n_chain;
float high_factor; // coverage cutoff set to high_factor*hom_cov
int max_n_chain; // fall-back max number of chains to consider
int k_mer_min_freq;
int k_mer_max_freq;
int load_index_from_disk;
@@ -56,7 +56,8 @@ extern hifiasm_opt_t asm_opt;
void init_opt(hifiasm_opt_t* asm_opt);
void destory_opt(hifiasm_opt_t* asm_opt);
void clear_opt(hifiasm_opt_t* asm_opt, int round);
void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round);
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);
htab.cpp
+4 -2
View File
@@ -771,7 +771,7 @@ void ha_ft_destroy(void *h)
* High-level interfaces *
*************************/
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs)
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov)
{
yak_ft_t *flt_tab;
int64_t cnt[YAK_N_COUNTS];
@@ -780,6 +780,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs)
h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN, NULL, NULL, rs);
ha_ct_hist(h, cnt, asm_opt->thread_num);
peak_hom = ha_analyze_count(YAK_N_COUNTS, 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);
cutoff = (int)(peak_hom * asm_opt->high_factor);
if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
@@ -791,7 +792,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs)
return (void*)flt_tab;
}
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs)
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov)
{
int64_t cnt[YAK_N_COUNTS], tot_cnt;
int peak_hom, peak_het, i, extra_flag1, extra_flag2;
@@ -812,6 +813,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
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, 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);
if (flt_tab == 0) {
int cutoff = (int)(peak_hom * asm_opt->high_factor);
htab.h
+2 -2
View File
@@ -26,11 +26,11 @@ extern const unsigned char seq_nt4_table[256];
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs);
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov);
int ha_ft_isflt(const void *hh, uint64_t y);
void ha_ft_destroy(void *h);
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs);
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs, int *hom_cov);
void ha_pt_destroy(ha_pt_t *h);
const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);