diff --git a/Assembly.cpp b/Assembly.cpp index df337b3..0844f85 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1717,7 +1717,7 @@ void Output_PAF() fprintf(stderr, "Writing PAF to disk ...... \n"); char* paf_name = (char*)malloc(strlen(asm_opt.output_file_name)+5); - sprintf(paf_name, "%s.paf", asm_opt.output_file_name); + sprintf(paf_name, "%s.ovlp.paf", asm_opt.output_file_name); FILE* output_file = fopen(paf_name, "w"); uint64_t i, j; ma_hit_t_alloc* sources = R_INF.paf; diff --git a/CommandLines.cpp b/CommandLines.cpp index b1ddb41..91d7ace 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -21,8 +21,7 @@ void Print_H(hifiasm_opt_t* asm_opt) { fprintf(stderr, "Usage: hifiasm [options] <...>\n"); fprintf(stderr, "Options:\n"); - ///fprintf(stderr, " -q FILE input in the fastq(.gz)/fasta(.gz) formats\n"); - fprintf(stderr, " -o FILE output assembly (in gfa format) and corrected reads [%s]\n", asm_opt->output_file_name); + fprintf(stderr, " -o FILE prefix of output files [%s]\n", asm_opt->output_file_name); fprintf(stderr, " -t INT number of threads [%d]\n", asm_opt->thread_num); fprintf(stderr, " -r INT round of correction [%d]\n", asm_opt->number_of_round); fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round); @@ -36,7 +35,7 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " -y FLOAT min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate); fprintf(stderr, " -v show version number\n"); fprintf(stderr, " -h show help information\n"); - fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -k 40 -t 32 -r 2 NA12878.fq.gz\n"); + fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n"); } @@ -247,4 +246,4 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) get_queries(argc, argv, &opt, asm_opt); return check_option(asm_opt); -} \ No newline at end of file +} diff --git a/Overlaps.cpp b/Overlaps.cpp index d16daca..bcc784d 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -7041,29 +7041,34 @@ const long long n_read) return 0; } -void ma_ug_print(const ma_ug_t *ug, All_reads *RNF, const ma_sub_t *coverage_cut, FILE *fp) +void ma_ug_print2(const ma_ug_t *ug, All_reads *RNF, const ma_sub_t *coverage_cut, int print_seq, FILE *fp) { uint32_t i, j, l; char name[32]; for (i = 0; i < ug->u.n; ++i) { // the Segment lines in GFA ma_utg_t *p = &ug->u.a[i]; sprintf(name, "utg%.6d%c", i + 1, "lc"[p->circ]); - fprintf(fp, "S\t%s\t%s\tLN:i:%d\n", name, p->s? p->s : "*", p->len); + if (print_seq) fprintf(fp, "S\t%s\t%s\tLN:i:%d\n", name, p->s? p->s : "*", p->len); + else fprintf(fp, "S\t%s\t*\tLN:i:%d\n", name, p->len); for (j = l = 0; j < p->n; l += (uint32_t)p->a[j++]) { uint32_t x = p->a[j]>>33; - fprintf(fp, "a\t%s\t%d\t%.*s(%u):%d-%d\t%c\t%d\n", name, l, - (int)Get_NAME_LENGTH((*RNF), x), Get_NAME((*RNF), x), x, - coverage_cut[x].s + 1, coverage_cut[x].e, "+-"[p->a[j]>>32&1], (uint32_t)p->a[j]); + fprintf(fp, "A\t%s\t%d\t%c\t%.*s\t%d\t%d\tid:i:%d\n", name, l, "+-"[p->a[j]>>32&1], (int)Get_NAME_LENGTH((*RNF), x), + Get_NAME((*RNF), x), coverage_cut[x].s, coverage_cut[x].e, x); } } for (i = 0; i < ug->g->n_arc; ++i) { // the Link lines in GFA uint32_t u = ug->g->arc[i].ul>>32, v = ug->g->arc[i].v; - fprintf(fp, "L\tutg%.6d%c\t%c\tutg%.6d%c\t%c\t%dM\tSD:i:%d\n", (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], + fprintf(fp, "L\tutg%.6d%c\t%c\tutg%.6d%c\t%c\t%dM\tL1:i:%d\n", (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], ug->g->arc[i].ol, asg_arc_len(ug->g->arc[i])); } } +void ma_ug_print(const ma_ug_t *ug, All_reads *RNF, const ma_sub_t *coverage_cut, FILE *fp) +{ + ma_ug_print2(ug, RNF, coverage_cut, 1, fp); +} + int asg_cut_internal(asg_t *g, int max_ext) { asg64_v a = {0,0,0}; @@ -7310,18 +7315,7 @@ void ma_sg_print(const asg_t *g, const All_reads *RNF, const ma_sub_t *sub, FILE void ma_ug_print_simple(const ma_ug_t *ug, All_reads *RNF, const ma_sub_t *coverage_cut, FILE *fp) { - uint32_t i; - char name[32]; - for (i = 0; i < ug->u.n; ++i) { // the Segment lines in GFA - ma_utg_t *p = &ug->u.a[i]; - sprintf(name, "utg%.6d%c", i + 1, "lc"[p->circ]); - fprintf(fp, "S\t%s\t%s\tLN:i:%d\n", name, "*", p->len); - } - for (i = 0; i < ug->g->n_arc; ++i) { // the Link lines in GFA - uint32_t u = ug->g->arc[i].ul>>32, v = ug->g->arc[i].v; - fprintf(fp, "L\tutg%.6d%c\t%c\tutg%.6d%c\t%c\t%dM\tSD:i:%d\n", (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], - (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], ug->g->arc[i].ol, asg_arc_len(ug->g->arc[i])); - } + ma_ug_print2(ug, RNF, coverage_cut, 0, fp); } @@ -7737,16 +7731,16 @@ void output_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_na ug = ma_ug_gen(sg); ma_ug_seq(ug, &R_INF, coverage_cut, n_read); - fprintf(stderr, "Writing unitig GFA to disk... \n"); + fprintf(stderr, "Writing raw unitig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+25); - sprintf(gfa_name, "%s.utg.gfa", output_file_name); + sprintf(gfa_name, "%s.r_utg.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); ma_ug_print(ug, &R_INF, coverage_cut, output_file); fclose(output_file); - // sprintf(gfa_name, "%s.simple.gfa", output_file_name); - // output_file = fopen(gfa_name, "w"); - // ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); - // fclose(output_file); + sprintf(gfa_name, "%s.r_utg.noseq.gfa", output_file_name); + output_file = fopen(gfa_name, "w"); + ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); + fclose(output_file); free(gfa_name); ma_ug_destroy(ug); @@ -7887,13 +7881,13 @@ void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sou All_reads *RNF, char* output_file_name) { char* gfa_name = (char*)malloc(strlen(output_file_name)+25); - sprintf(gfa_name, "%s.gfa.aux", output_file_name); + sprintf(gfa_name, "%s.ovlp", output_file_name); write_All_reads(RNF, gfa_name); - sprintf(gfa_name, "%s.gfa.aux.source", output_file_name); + sprintf(gfa_name, "%s.ovlp.source", output_file_name); write_ma_hit_ts(sources, RNF->total_reads, gfa_name); - sprintf(gfa_name, "%s.gfa.aux.reverse", output_file_name); + sprintf(gfa_name, "%s.ovlp.reverse", output_file_name); write_ma_hit_ts(reverse_sources, RNF->total_reads, gfa_name); free(gfa_name); @@ -7903,20 +7897,20 @@ int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_s char* output_file_name) { char* gfa_name = (char*)malloc(strlen(output_file_name)+25); - sprintf(gfa_name, "%s.gfa.aux", output_file_name); + sprintf(gfa_name, "%s.ovlp.aux", output_file_name); if(!load_All_reads(&R_INF, gfa_name)) { return 0; } - sprintf(gfa_name, "%s.gfa.aux.source", output_file_name); + sprintf(gfa_name, "%s.ovlp.source", output_file_name); if(!load_ma_hit_ts(sources, gfa_name)) { return 0; } - sprintf(gfa_name, "%s.gfa.aux.reverse", output_file_name); + sprintf(gfa_name, "%s.ovlp.reverse", output_file_name); if(!load_ma_hit_ts(reverse_sources, gfa_name)) { return 0; @@ -8770,17 +8764,17 @@ char* output_file_name, long long n_read, long long bubble_dist, long long tipsL ug = ma_ug_gen_primary(sg, 0); ma_ug_seq(ug, &R_INF, coverage_cut, n_read); - fprintf(stderr, "Writing unitig GFA to disk... \n"); + fprintf(stderr, "Writing processed unitig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+35); - sprintf(gfa_name, "%s.wsb.utg.gfa", output_file_name); + sprintf(gfa_name, "%s.p_utg.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); ma_ug_print(ug, &R_INF, coverage_cut, output_file); fclose(output_file); - // sprintf(gfa_name, "%s.simple.no_s_bub.gfa", output_file_name); - // output_file = fopen(gfa_name, "w"); - // ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); - // fclose(output_file); + sprintf(gfa_name, "%s.p_utg.noseq.gfa", output_file_name); + output_file = fopen(gfa_name, "w"); + ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); + fclose(output_file); free(gfa_name); ma_ug_destroy(ug); @@ -8883,21 +8877,17 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen) ug = ma_ug_gen_primary(sg, 0); ma_ug_seq(ug, &R_INF, coverage_cut, n_read); - fprintf(stderr, "Writing unitig GFA to disk... \n"); + fprintf(stderr, "Writing primary contig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+35); - sprintf(gfa_name, "%s.ctg.gfa", output_file_name); + sprintf(gfa_name, "%s.p_ctg.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); ma_ug_print(ug, &R_INF, coverage_cut, output_file); fclose(output_file); - - - - - // sprintf(gfa_name, "%s.simple.contig.gfa", output_file_name); - // output_file = fopen(gfa_name, "w"); - // ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); - // fclose(output_file); + sprintf(gfa_name, "%s.p_ctg.noseq.gfa", output_file_name); + output_file = fopen(gfa_name, "w"); + ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); + fclose(output_file); free(gfa_name); ma_ug_destroy(ug); @@ -8910,17 +8900,17 @@ void output_contig_graph_alternative(asg_t *sg, ma_sub_t* coverage_cut, char* ou ug = ma_ug_gen_primary(sg, 1); ma_ug_seq(ug, &R_INF, coverage_cut, n_read); - fprintf(stderr, "Writing unitig GFA to disk... \n"); + fprintf(stderr, "Writing alternate contig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+35); - sprintf(gfa_name, "%s.alter.ctg.gfa", output_file_name); + sprintf(gfa_name, "%s.a_ctg.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); ma_ug_print(ug, &R_INF, coverage_cut, output_file); fclose(output_file); - // sprintf(gfa_name, "%s.simple.alter.ctg.gfa", output_file_name); - // output_file = fopen(gfa_name, "w"); - // ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); - // fclose(output_file); + sprintf(gfa_name, "%s.a_ctg.noseq.gfa", output_file_name); + output_file = fopen(gfa_name, "w"); + ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); + fclose(output_file); free(gfa_name); ma_ug_destroy(ug); diff --git a/README.md b/README.md index 859c306..9849247 100644 --- a/README.md +++ b/README.md @@ -2,24 +2,46 @@ ```sh # Install hifiasm (requiring g++ and zlib) -git clone https://github.com/chhylp123/hifiasm.git +git clone https://github.com/chhylp123/hifiasm cd hifiasm && make # Assembly -./hifiasm -t 32 NA12878.fq.gz +./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz ``` ## Introduction -Hifiasm is an ultrafast haplotype-resolved de novo assembler based on PacBio Hifi reads. Unlike most existing assemblers, hifiasm starts from uncollapsed genome. Thus, it is able to keep the haplotype information as much as possible. The input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and its outputs consist of: -1. Haplotype-resolved assembly [unitig][unitig] graph in [GFA][gfa] format (hifiasm.asm.utg.gfa in dafault). -2. Haplotype-resolved assembly [unitig][unitig] graph in [GFA][gfa] format without small bubbles (hifiasm.asm.wsb.utg.gfa in dafault). Small bubbles might be caused by somatic mutations, which are useless for some applications. -3. Primary assembly [contig][unitig] graph in [GFA][gfa] format (hifiasm.asm.ctg.gfa in dafault). -4. Alternate assembly [contig][unitig] graph in [GFA][gfa] format (hifiasm.asm.alter.ctg.gfa in dafault). -5. Haplotype-aware error corrected reads in fasta format (hifiasm.asm.ec.fa in dafault). + +Hifiasm is an ultrafast haplotype-resolved de novo assembler based on PacBio +Hifi reads. Unlike most existing assemblers, hifiasm starts from uncollapsed +genome. Thus, it is able to keep the haplotype information as much as possible. +The input of hifiasm is the PacBio Hifi reads in fasta/fastq format, and its +outputs consist of: + +1. Haplotype-resolved raw [unitig][unitig] graph in [GFA][gfa] format + (hifiasm.asm.r\_utg.gfa by default). +2. Haplotype-resolved processed [unitig][unitig] graph in [GFA][gfa] format + without small bubbles (hifiasm.asm.p\_utg.gfa by default). Small bubbles + might be caused by somatic mutations, which are useless for some + applications. +3. Primary assembly [contig][unitig] graph in [GFA][gfa] format + (hifiasm.asm.p\_ctg.gfa by default). +4. Alternate assembly [contig][unitig] graph in [GFA][gfa] format + (hifiasm.asm.a\_ctg.gfa by default). +5. Haplotype-aware error corrected reads in fasta format (hifiasm.asm.ec.fa by + default). 6. All-to-all overlaps in [paf][paf] format (hifiasm.asm.paf). -So far hifiasm is still in early development stage, it will output phased chromosome-level high-quality assembly in the near future. In addition, hifiasm also outputs three binary files that save all overlap inforamtion (hifiasm.asm.gfa.aux.bin, hifiasm.asm.gfa.aux.reverse.bin, hifiasm.asm.gfa.aux.source.bin in default). With these files, hifiasm can avoid the time-consuming all-to-all overlap calculation step, and do the assembly directly and quickly. This might be helpful when you want to get an optimized assembly by multiple round of experiments with different parameters. +So far hifiasm is still in early development stage, it will output phased +chromosome-level high-quality assembly in the near future. In addition, hifiasm +also outputs three binary files that save all overlap inforamtion +(hifiasm.asm.gfa.aux.bin, hifiasm.asm.gfa.aux.reverse.bin, +hifiasm.asm.gfa.aux.source.bin in default). With these files, hifiasm can avoid +the time-consuming all-to-all overlap calculation step, and do the assembly +directly and quickly. This might be helpful when you want to get an optimized +assembly by multiple round of experiments with different parameters. -Hifiasm is a standalone and lightweight assembler, which does not need external libraries (except zlib). For large genomes, it can generate high-quality assembly in a few hours. Hifiasm has been tested on the following datasets: +Hifiasm is a standalone and lightweight assembler, which does not need external +libraries (except zlib). For large genomes, it can generate high-quality +assembly in a few hours. Hifiasm has been tested on the following datasets: |Dataset|GSize|Cov|Asm options|CPU time|Wall time|RAM|[unitig][unitig]/[contig][unitig] N50[1]| |:---------------|-----:|-----:|:---------------------|-------:|--------:|----:|----------------:| @@ -33,35 +55,48 @@ Hifiasm is a standalone and lightweight assembler, which does not need external [3] Butterfly has high heterozygous rate, so that most chromosomes have been fully separated into two haplotypes. In this case, contig N50 makes no sense. ## Usage + For Hifi reads assembly, a typical command line looks like: ```sh -./hifiasm -o NA12878.asm -k 40 -t 32 -r 2 NA12878.fq.gz +./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz ``` -where `NA12878.fq.gz` is the input reads and `-o` specifies the output files. In this example, all output files can be found at `NA12878.asm.*`. `-k`, `-t` and `-r` specify the length of k-mer, the number of CPU threads, and the number of correction rounds, respectively. Note that at first run, hifiasm will save all overlaps to disk, which can avoid the time-consuming all-to-all overlap calculation next time. For hifiasm, once the overlap information has been obtained during the previous run in advance, it is able to load all overlaps from disk and then directly do assembly. If you want to ignore the pre-computed overlap information, please specify `-i` or simply delete `*gfa.aux.bin`, `*gfa.aux.reverse.bin` and `*gfa.aux.source.bin`. +where `NA12878.fq.gz` is the input reads and `-o` specifies the output files. +In this example, all output files can be found at `NA12878.asm.*`. `-k`, `-t` +and `-r` specify the length of k-mer, the number of CPU threads, and the number +of correction rounds, respectively. Note that at first run, hifiasm will save +all overlaps to disk, which can avoid the time-consuming all-to-all overlap +calculation next time. For hifiasm, once the overlap information has been +obtained during the previous run in advance, it is able to load all overlaps +from disk and then directly do assembly. If you want to ignore the pre-computed +overlap information, please specify `-i` or simply delete `*gfa.aux.bin`, +`*gfa.aux.reverse.bin` and `*gfa.aux.source.bin`. -Please note that some old Hifi reads may consist of short adapters. To improve the assembly quality, adapters should be removed by `-z` as follow: +Please note that some old Hifi reads may consist of short adapters. To improve +the assembly quality, adapters should be removed by `-z` as follow: ```sh -./hifiasm -k 40 -t 42 -r 2 -z 20 butterfly.fq.gz +./hifiasm -o butterfly.asm -t 42 -z 20 butterfly.fq.gz ``` In this example, hifiasm will remove 20 bases from both ends of each read. - - [unitig]: http://wgs-assembler.sourceforge.net/wiki/index.php/Celera_Assembler_Terminology [gfa]: https://github.com/pmelsted/GFA-spec/blob/master/GFA-spec.md [paf]: https://github.com/lh3/miniasm/blob/master/PAF.md ## Getting Help -The `-h` option of hifiasm provides detailed description of options. If you have further questions, -please raise an issue at the issue page. -## Limitations and future works -1. For genome with low heterozygous rate, hifiasm only outputs haplotype-resolved assembly graph, instead of the phased chromosome-level assembly (will support such output in the near future). +The `-h` option of hifiasm provides detailed description of options. If you +have further questions, please raise an issue at the issue page. + +## Limitations and future works + +1. For genome with low heterozygous rate, hifiasm only outputs + haplotype-resolved assembly graph, instead of the phased chromosome-level + assembly (will support such output in the near future). 2. The running time and memory usage should be further reduced. -3. The N50 should be further improved. \ No newline at end of file +3. The N50 should be further improved.