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

update 0.14-r310

Browse Source
This commit is contained in:
chhylp123
2021-02-13 23:21:45 -05:00
parent 5bb22f701c
commit 18086e2c2e
6 changed files with 1329 additions and 175 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.cpp
+38 -3
View File
@@ -30,6 +30,8 @@ static ko_longopt_t long_options[] = {
{ "h2", ko_required_argument, 315 },
{ "enzyme", ko_required_argument, 316 },
{ "b-cov", ko_required_argument, 317 },
{ "h-cov", ko_required_argument, 318 },
{ "m-rate", ko_required_argument, 319 },
{ 0, 0, 0 }
};
@@ -68,7 +70,12 @@ void Print_H(hifiasm_opt_t* asm_opt)
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");
fprintf(stderr, " break contigs at breakpoints with coverage drop at <INT-fold coverage [%d]\n", asm_opt->break_cov);
fprintf(stderr, " break contigs at positions with <INT-fold coverage; work with '--m-rate'; 0 to disable [%d]\n", asm_opt->b_low_cov);
fprintf(stderr, " --h-cov INT\n");
fprintf(stderr, " break contigs at positions with >INT-fold coverage; work with '--m-rate'; -1 to disable [%d]\n", asm_opt->b_high_cov);
fprintf(stderr, " --m-rate FLOAT\n");
fprintf(stderr, " break contigs at positions with <=FLOAT*coverage exact overlaps;\n");
fprintf(stderr, " only work with '--b-cov' or '--h-cov'[%.2f]\n", asm_opt->m_rate);
// fprintf(stderr, " --pri-range INT1[,INT2]\n");
// fprintf(stderr, " keep contigs with coverage in this range in p_ctg.gfa; -1 to disable [auto,inf]\n");
@@ -153,7 +160,9 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->hic_inconsist_rate = 30;
///asm_opt->bub_mer_length = 3;
asm_opt->bub_mer_length = 1000000;
asm_opt->break_cov = 0;
asm_opt->b_low_cov = 0;
asm_opt->b_high_cov = -1;
asm_opt->m_rate = 0.75;
}
void destory_enzyme(enzyme* f)
@@ -419,6 +428,30 @@ int check_option(hifiasm_opt_t* asm_opt)
return 0;
}
if(asm_opt->b_low_cov < 0)
{
fprintf(stderr, "[ERROR] must >= 0 (--b-cov)\n");
return 0;
}
if(asm_opt->b_high_cov != -1 && asm_opt->b_high_cov < 0)
{
fprintf(stderr, "[ERROR] must >= 0 (--h-cov)\n");
return 0;
}
if(asm_opt->m_rate < 0)
{
fprintf(stderr, "[ERROR] must >= 0 (--m-rate)\n");
return 0;
}
if(asm_opt->b_high_cov != -1 && asm_opt->b_high_cov <= asm_opt->b_low_cov)
{
fprintf(stderr, "[ERROR] [--h-cov] must >= [--b-cov]\n");
return 0;
}
// fprintf(stderr, "input file num: %d\n", asm_opt->num_reads);
// fprintf(stderr, "output file: %s\n", asm_opt->output_file_name);
// fprintf(stderr, "number of threads: %d\n", asm_opt->thread_num);
@@ -592,7 +625,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 314) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[0]), 0);
else if (c == 315) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[1]), 0);
else if (c == 316) get_hic_enzymes(opt.arg, &(asm_opt->hic_enzymes), 1);
else if (c == 317) asm_opt->break_cov = atoi(opt.arg);
else if (c == 317) asm_opt->b_low_cov = atoi(opt.arg);
else if (c == 318) asm_opt->b_high_cov = atoi(opt.arg);
else if (c == 319) asm_opt->m_rate = atof(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);
CommandLines.h
+4 -2
View File
@@ -3,7 +3,7 @@
#include <pthread.h>
#define HA_VERSION "0.14-r309"
#define HA_VERSION "0.14-r310"
#define VERBOSE 0
@@ -49,7 +49,9 @@ typedef struct {
double max_ov_diff_final;
int hom_cov;
int het_cov;
int break_cov;
int b_low_cov;
int b_high_cov;
double m_rate;
int max_n_chain; // fall-back max number of chains to consider
int min_hist_kmer_cnt;
int load_index_from_disk;
Overlaps.cpp
+1243 -41
View File
File diff suppressed because it is too large Load Diff
Overlaps.h
+1 -1
View File
@@ -477,7 +477,7 @@ void destory_R_to_U(R_to_U* x);
void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig, uint8_t* flag);
void get_R_to_U(R_to_U* x, uint32_t rID, uint32_t* uID, uint32_t* is_Unitig);
void transfor_R_to_U(R_to_U* x);
void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, int require_equal_nv, int test_tangle);
void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, kvec_asg_arc_t_warp* edge, int require_equal_nv, int test_tangle);
int asg_pop_bubble_primary(asg_t *g, int max_dist);
long long asg_arc_del_simple_circle_untig(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *g, long long circleLen, int is_drop);
hic.cpp
+5 -127
View File
@@ -7369,7 +7369,7 @@ void init_chain_hic_warp(ma_ug_t* ug, hc_links* link, bubble_type* bub, chain_hi
}
}
fprintf(stderr, "# chain: %u, # pre chain: %u\n", m, (uint32_t)(*c_w).n);
///fprintf(stderr, "# chain: %u, # pre chain: %u\n", m, (uint32_t)(*c_w).n);
(*c_w).n = m;
for (i = 0; i < (*c_w).n; i++)
{
@@ -11272,7 +11272,7 @@ int alignment_worker_pipeline(sldat_t* sl, const enzyme *fn1, const enzyme *fn2)
kt_pipeline(3, worker_pipeline, sl, 3);
fprintf(stderr, "fn1->a[i]: %s, fn2->a[i]: %s, sl->hits.a.n: %u\n", fn1->a[i], fn2->a[i], (uint32_t)sl->hits.a.n);
///fprintf(stderr, "fn1->a[i]: %s, fn2->a[i]: %s, sl->hits.a.n: %u\n", fn1->a[i], fn2->a[i], (uint32_t)sl->hits.a.n);
kseq_destroy(sl->ks1);
@@ -11281,11 +11281,11 @@ int alignment_worker_pipeline(sldat_t* sl, const enzyme *fn1, const enzyme *fn2)
gzclose(fp2);
}
fprintf(stderr, "+sl->hits.a.n: %u\n", (uint32_t)sl->hits.a.n);
///fprintf(stderr, "+sl->hits.a.n: %u\n", (uint32_t)sl->hits.a.n);
dedup_hits(&(sl->hits));
fprintf(stderr, "-sl->hits.a.n: %u\n", (uint32_t)sl->hits.a.n);
///fprintf(stderr, "-sl->hits.a.n: %u\n", (uint32_t)sl->hits.a.n);
return 1;
}
@@ -11401,127 +11401,6 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx)
return 1;
}
int hic_short_align_back(const char *fn1, const char *fn2, ha_ug_index* idx)
{
double index_time = yak_realtime();
sldat_t sl;
gzFile fp1, fp2;
kvec_hc_edge back_hc_edge;
kv_init(back_hc_edge.a);
if ((fp1 = gzopen(fn1, "r")) == 0) return 0;
if ((fp2 = gzopen(fn2, "r")) == 0) return 0;
sl.ks1 = kseq_init(fp1);
sl.ks2 = kseq_init(fp2);
sl.idx = idx;
sl.link = idx->link;
sl.chunk_size = 20000000;
sl.n_thread = asm_opt.thread_num;
sl.total_base = sl.total_pair = 0;
idx->max_cnt = 5;
kv_init(sl.hits.a);
if(!load_hc_hits(&sl.hits, asm_opt.output_file_name))
{
/*******************************for debug************************************/
// load_reads(&R1, fn1);
// test_reads(&R1, fn1);
// load_reads(&R2, fn2);
// test_reads(&R1, fn1);
/*******************************for debug************************************/
fprintf(stderr, "+sl->hits.a.n: %u\n", (uint32_t)sl.hits.a.n);
kt_pipeline(3, worker_pipeline, &sl, 3);
fprintf(stderr, "-sl->hits.a.n: %u\n", (uint32_t)sl.hits.a.n);
fprintf(stderr, "fn1: %s, fn2: %s\n", fn1, fn2);
/*******************************for debug************************************/
// sort_hits(&sl.hits);
// print_hits(idx, &sl.hits, fn1);
/*******************************for debug************************************/
dedup_hits(&sl.hits);
write_hc_hits(&sl.hits, asm_opt.output_file_name);
}
fprintf(stderr, "u.n: %d, uID_bits: %lu, pos_bits: %lu, sl.hits.a.n: %u\n", (uint32_t)idx->ug->u.n, idx->uID_bits, idx->pos_bits, (uint32_t)sl.hits.a.n);
H_partition hap;
MT M;
init_MT(&M, idx->ug->g->n_seq<<1);
bubble_type bub;
memset(&bub, 0, sizeof(bubble_type));
bub.round_id = 0; bub.n_round = 2;
for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++)
{
identify_bubbles(idx->ug, &bub, idx->link);
if(bub.round_id == 0)
{
collect_hc_links(sl.idx, &sl.hits, idx->link, &bub, &M);
collect_hc_reverse_links(idx->link, idx->ug, &bub);
}
init_hic_p((ha_ug_index*)sl.idx, &sl.hits, idx->link, &bub, &back_hc_edge, &M, &hap, 0);
///init_hic_p_new((ha_ug_index*)sl.idx, &sl.hits, idx->link, &bub, &back_hc_edge, &M);
reset_H_partition(&hap, (bub.round_id == 0? 1 : 0));
init_contig_partition(&hap, idx, &bub);
phasing_improvement(&hap, &(hap.g_p), idx, &bub);
label_unitigs(&(hap.g_p), idx->ug);
///print_hc_links(idx->link, 0, &hap);
}
cluster_contigs(&bub, idx, &sl.hits, &M, &hap);
destory_MT(&M);
///print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, idx->link, idx);
///print_hits(idx, &sl.hits, fn1);
///print_debug_bubble_graph(&bub, idx->ug, asm_opt.output_file_name);
// print_bubble_chain(&bub);
// print_hc_links(idx->link, 0, &hap);
///print_contig_partition(&hap, "final");
// uint32_t i;
// for (i = 0; i < idx->ug->g->n_seq; i++)
// {
// fprintf(stderr, "utg%.6ul, index: %u\n", (int)(i+1), bub.index[i]);
// }
destory_contig_partition(&hap);
kv_destroy(back_hc_edge.a);
return 1;
/*******************************for debug************************************/
// destory_reads(&R1);
// destory_reads(&R2);
/*******************************for debug************************************/
print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, idx->link, idx);
collect_hc_reverse_links(idx->link, idx->ug, &bub);
normalize_hc_links(idx->link);
/*******************************for debug************************************/
///print_hc_links(&link);
/*******************************for debug************************************/
min_cut_t* cut = clean_hap(idx->link, &bub, idx->ug);
///print_bubbles(idx->ug, &bub, NULL, &link, idx);
G_partition* gp = clean_bubbles(idx->link, &bub, cut, idx->ug);
///print_hc_links(&link);
destory_min_cut_t(cut); free(cut);
destory_G_partition(gp); free(gp);
kv_destroy(sl.hits.a);
destory_bubbles(&bub);
kseq_destroy(sl.ks1);
kseq_destroy(sl.ks2);
gzclose(fp1);
gzclose(fp2);
fprintf(stderr, "[M::%s::%.3f] processed %lu pairs; %lu bases\n", __func__, yak_realtime()-index_time, sl.total_pair, sl.total_base);
return 1;
}
void hic_analysis(ma_ug_t *ug, asg_t* read_g, hc_links* link)
{
@@ -11533,8 +11412,7 @@ void hic_analysis(ma_ug_t *ug, asg_t* read_g, hc_links* link)
ug_index->read_g = read_g;
ug_index->link = link;
///test_unitig_index(ug_index, ug);
///hic_short_align(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index);
hic_short_align_back(asm_opt.hic_reads[0]->a[0], asm_opt.hic_reads[1]->a[0], ug_index);
hic_short_align(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index);
destory_hc_pt_index(ug_index);
}
hifiasm.1
+38 -1
View File
@@ -1,4 +1,4 @@
.TH hifiasm 1 "19 July 2020" "hifiasm-0.9 (r289)" "Bioinformatics tools"
.TH hifiasm 1 "13 Feb 2021" "hifiasm-0.14 (r310)" "Bioinformatics tools"
.SH NAME
.PP
@@ -212,6 +212,32 @@ with suffix
.B lowQ.bed
[70]. Set 0 to disable.
.TP
.BI --b-cov \ INT
Break contigs at potential misassemblies with <INT-fold coverage [0].
Work with
.B --m-rate.
Set 0 to disable.
.TP
.BI --h-cov \ INT
Break contigs at potential misassemblies with >INT-fold coverage [-1].
Work with
.B --m-rate.
Set -1 to disable.
.TP
.BI --m-rate \ FLOAT
Break contigs with <=FLOAT*coverage exact overlaps [0.75].
Only work with
.B --b-cov
and
.B --h-cov.
.SS Trio-partition options
.TP 10
@@ -289,6 +315,17 @@ For ordinary samples, no need to enable this mode [experimental, not stable].
Write additional files to speed up the debugging of graph cleaning.
.SS Hi-C-partition options [experimental, not stable]
.TP
.BI --h1 \ FILEs
File names of input Hi-C R1 [r1_1.fq,r1_2.fq,...]
.TP
.BI --h2 \ FILEs
File names of input Hi-C R2 [r2_1.fq,r2_2.fq,...]
.SH OUTPUTS
.PP