From 7a35bd7fcc69651040954bae19fe6b85653174f8 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 9 Apr 2020 20:40:01 -0400 Subject: [PATCH] r213: renamed trio binning output files --- CommandLines.cpp | 10 +-- CommandLines.h | 2 +- Overlaps.cpp | 13 ++-- Trio.cpp | 6 +- hifiasm.1 | 156 +++++++++++++++++++++++++---------------------- 5 files changed, 97 insertions(+), 90 deletions(-) diff --git a/CommandLines.cpp b/CommandLines.cpp index 5ff1472..fd4f849 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -50,8 +50,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\" []\n"); - fprintf(stderr, " -M FILE maternal trio index generated by \"yak count\" []\n"); + fprintf(stderr, " -1 FILE hap1/paternal k-mer dump generated by \"yak count\" []\n"); + fprintf(stderr, " -2 FILE hap2/maternal k-mer dump 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); @@ -313,7 +313,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:w: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:1:2:", long_options)) >= 0) { if (c == 'h') { Print_H(asm_opt); @@ -339,8 +339,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 'b') asm_opt->required_read_name = opt.arg; else if (c == 'c') asm_opt->min_cnt = atoi(opt.arg); else if (c == 'd') asm_opt->mid_cnt = atoi(opt.arg); - else if (c == 'P') asm_opt->pat_index = opt.arg; - else if (c == 'M') asm_opt->mat_index = opt.arg; + else if (c == '1' || c == 'P') asm_opt->pat_index = opt.arg; // -P/-M reserved for backward compatibility + else if (c == '2' || c == 'M') asm_opt->mat_index = opt.arg; else if (c == 'x') asm_opt->max_drop_rate = atof(opt.arg); else if (c == 'y') asm_opt->min_drop_rate = atof(opt.arg); else if (c == 'p') asm_opt->small_pop_bubble_size = atoll(opt.arg); diff --git a/CommandLines.h b/CommandLines.h index e5901cc..7e251df 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.3.0-dirty-r212" +#define HA_VERSION "0.3.0-dirty-r213" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index f1b90bf..a715cf6 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -9664,13 +9664,13 @@ void ma_ug_print2(const ma_ug_t *ug, All_reads *RNF, const ma_sub_t *coverage_cu uint32_t x = p->a[j]>>33; if(xtotal_reads) { - fprintf(fp, "A\t%s\t%d\t%c\t%.*s\t%d\t%d\tid:i:%d\t%c\n", name, l, "+-"[p->a[j]>>32&1], + fprintf(fp, "A\t%s\t%d\t%c\t%.*s\t%d\t%d\tid:i:%d\tHG:A:%c\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, "apmaaa"[RNF->trio_flag[x]]); } else { - fprintf(fp, "A\t%s\t%d\t%c\t%s\t%d\t%d\tid:i:%d\t%c\n", name, l, "+-"[p->a[j]>>32&1], + fprintf(fp, "A\t%s\t%d\t%c\t%s\t%d\t%d\tid:i:%d\tHG:A:%c\n", name, l, "+-"[p->a[j]>>32&1], "FAKE", coverage_cut[x].s, coverage_cut[x].e, x, '*'); } @@ -13527,7 +13527,7 @@ long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIn float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp) { char* gfa_name = (char*)malloc(strlen(output_file_name)+100); - sprintf(gfa_name, "%s.%s.r_utg.gfa", output_file_name, (flag==FATHER?"p":"m")); + sprintf(gfa_name, "%s.%s.p_ctg.gfa", output_file_name, (flag==FATHER?"hap1":"hap2")); fprintf(stderr, "Writing %s to disk... \n", gfa_name); FILE* output_file = fopen(gfa_name, "w"); @@ -13548,7 +13548,7 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp) ma_ug_print(ug, &R_INF, coverage_cut, output_file); fclose(output_file); - sprintf(gfa_name, "%s.%s.r_utg.noseq.gfa", output_file_name, (flag==FATHER?"p":"m")); + sprintf(gfa_name, "%s.%s.p_ctg.noseq.gfa", output_file_name, (flag==FATHER?"hap1":"hap2")); output_file = fopen(gfa_name, "w"); ma_ug_print_simple(ug, &R_INF, coverage_cut, output_file); fclose(output_file); @@ -26495,7 +26495,10 @@ ma_sub_t **coverage_cut_ptr, int debug_g) if ((asm_opt.pat_index && asm_opt.mat_index)) { - output_unitig_graph(sg, coverage_cut, output_file_name, sources, max_hang_length, mini_overlap_length); + char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); + sprintf(buf, "%s.dip", output_file_name); + output_unitig_graph(sg, coverage_cut, buf, sources, max_hang_length, mini_overlap_length); + free(buf); output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, reverse_sources, bubble_dist, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, diff --git a/Trio.cpp b/Trio.cpp index 7e9d974..adef200 100644 --- a/Trio.cpp +++ b/Trio.cpp @@ -255,9 +255,6 @@ void trio_partition(void) return; } - double start_time = Get_T(); - fprintf(stderr, "Start trio binning ...... \n"); - yak_ch_t *ch; int i /**, min_cnt = 2, mid_cnt = 5**/; tb_shared_t aux; @@ -286,6 +283,5 @@ void trio_partition(void) free(aux.bseq); yak_ch_destroy(ch); - fprintf(stderr, "Trio binning has been done.\n"); - fprintf(stderr, "%-30s%18.2f\n\n", "Trio binning time:", Get_T() - start_time); + fprintf(stderr, "[M::%s::%.3f*%.2f] ==> partitioned reads using yak dumps\n", __func__, yak_realtime(), yak_cpu_usage()); } diff --git a/hifiasm.1 b/hifiasm.1 index fa3a6c3..2d0c2ba 100644 --- a/hifiasm.1 +++ b/hifiasm.1 @@ -8,7 +8,7 @@ hifiasm - haplotype-resolved de novo assembler for PacBio Hifi reads. .PP hifiasm .RB [ -o -.IR outPrefix ] +.IR prefix ] .RB [ -t .IR numThres ] .RB [ -r @@ -51,8 +51,8 @@ outputs consist of multiple types of assembly graph in GFA format. .BI -o \ FILE Prefix of output files [hifiasm.asm]. The outputs of hifiasm include error corrected reads in fasta format, all-to-all overlaps in paf format, and four types of assembly -graph in GFA format. For detailed description of all assembly graphs, please see -.I 'Outputs' +graph in GFA format. For detailed description of all assembly graphs, please see the +.B OUTPUTS section of this man-page. .TP 10 @@ -83,11 +83,11 @@ Rounds of haplotype-aware error corrections [2]. This option affects all outputs .TP 10 .BI -a \ INT Rounds of assembly graph cleaning [4]. This option is used with -.I [-x maxDropRatio] +.B -x and -.I [-y minDropRatio]. +.BR -y . Note that unlike -.I [-r], +.BR -r , this option does not affect error corrected reads and all-to-all overlaps. @@ -103,11 +103,11 @@ significantly improve the assembly quality. .TP 10 .BI -m \ INT -Maximal probing distance for bubble popping when generating primary/alternate assembly +Maximal probing distance for bubble popping when generating primary/alternate contig graphs [10000000]. Bubbles longer than .I INT -bases will not be popped. For detailed description of these graphs, please see -.I 'Outputs' +bases will not be popped. For detailed description of these graphs, please see the +.B OUTPUTS section of this man-page. @@ -117,8 +117,8 @@ Maximal probing distance for bubble popping when generating haplotype-resolved p without small bubbles [100000]. Bubbles longer than .I INT bases will not be popped. Small bubbles might be caused by somatic mutations or noise in data, which -are not the real haplotype information. For detailed description of this graph, please see -.I 'Outputs' +are not the real haplotype information. For detailed description of this graph, please see the +.B OUTPUTS section of this man-page. @@ -133,7 +133,7 @@ reads [3]. Hifiasm may try to remove small unitigs at various steps. .TP 10 .BI -x \ FLOAT, -y \ FLOAT Max and min overlap drop ratio [0.8, 0.2]. This option is used with -.I [-r roundCorrection]. +.BR -r . Given a node .I N in the assembly graph, let max(N) @@ -143,15 +143,15 @@ Hifiasm iteratively drops overlaps of .I N if their length / max(N) are below a threshold controlled by -.I [-x maxDropRatio] +.B -x and -.I [-y minDropRatio]. +.BR -y . Hifiasm applies -.I [-r roundCorrection] +.B -r rounds of short overlap removal with an increasing threshold between -.I [-x maxDropRatio] +.B -x and -.I [-y minDropRatio]. +.BR -y . .TP 10 .BI -i @@ -167,38 +167,37 @@ with different parameters. .SS Trio-partition options .TP 10 -.BI -P \ FILE -Paternal trio index. This index should be generated by -.I [yak count] -with the paternal short reads. For details of yak, please see +.BI -1 \ FILE +Paternal/haplotype1 k-mer dump generated by +.B yak count +from the paternal/haplotype1 reads. For details of yak, please see .I [https://github.com/lh3/yak] - .TP 10 -.BI -M \ FILE -Maternal trio index. This index should be generated by -.I [yak count] -with the maternal short reads. For details of yak, please see +.BI -2 \ FILE +Maternal/haplotype2 k-mer dump generated by +.B yak count +from the maternal/haplotype2 reads. For details of yak, please see .I [https://github.com/lh3/yak] .TP 10 .BI -c \ INT Lower bound of the binned k-mer's frequency [2]. When doing trio binning, a k-mer is said to be differentiating if it occurs >= -.I [-d upper_bound] +.B -d times in one sample but occurs < -.I [-c lower_bound] +.B -c times in the other sample. .TP 10 .BI -d \ INT Upper bound of the binned k-mer's frequency [5]. When doing trio binning, a k-mer is said to be differentiating if it occurs >= -.I [-d upper_bound] +.B -d times in one sample but occurs < -.I [-c lower_bound] +.B -c times in the other sample. @@ -235,60 +234,69 @@ Build maternal trio index from mat.fq.gz. Build paternal trio index from pat.fq.gz. .TP -.BR ./hifiasm " " \-o " " NA12878.asm " " \-t " " 32 " " \-P " " pat.yak " " \-M " " mat.yak " " NA12878_1.fq.gz " " NA12878_2.fq.gz -In this example, hifiasm will do trio assembly with 32 CPU threads. The paternal assembly can be found at [NA12878.asm.p.r_utg.gfa], -and the maternal assembly can be found at [NA12878.asm.m.r_utg.gfa]. +.BR ./hifiasm " " \-o " " NA12878.asm " " \-t " " 32 " " \-1 " " pat.yak " " \-2 " " mat.yak " " NA12878_1.fq.gz " " NA12878_2.fq.gz +In this example, hifiasm will do trio assembly with 32 CPU threads. The paternal assembly can be found at [NA12878.asm.hap1.p_ctg.gfa], +and the maternal assembly can be found at [NA12878.asm.hap2.p_ctg.gfa]. .SH OUTPUTS +.PP +Without trio partition options +.B -1 +and +.BR -2 , +hifiasm generates the following assembly graphs in the GFA format: + +.RS 2 +.TP 2 +* +.IR prefix .r_utg.gfa: +haplotype-resolved raw unitig graph. This graph keeps all haplotype information. + +.TP +* +.IR prefix .p_utg.gfa: +haplotype-resolved processed unitig graph without small bubbles. Small bubbles +might be caused by somatic mutations or noise in data, which are not the real +haplotype information. The size of popped small bubbles should be specified by +.BR -p . + +.TP +* +.IR prefix .p_ctg.gfa: +assembly graph of primary contigs. This graph collapses different haplotypes. + +.TP +* +.IR prefix .a_ctg.gfa: +assembly graph of alternate contigs. This graph consists of all assemblies that +are discarded in primary contig graph. + +.RE .PP -Consider the prefix of output files has been specified by -.I [-o outPrefix]. -During the error correction step, hifiasm outputs the following two files: +With trio partition, hifiasm outputs the following assembly graphs: -.IP -1. Haplotype-aware error corrected reads in fasta format [outPrefix.ec.fa]. +.RS 2 +.TP 2 +* +.IR prefix .dip.r_utg.gfa: +haplotype-resolved raw unitig graph. This graph keeps all haplotype information. -2. All-to-all overlaps in paf format [outPrefix.ovlp.paf]. - -.PP -During the non-trio assembly step, hifiasm outputs the following four assembly graphs in GFA format: - - -.IP -1. Haplotype-resolved raw unitig graph [outPrefix.r_utg.gfa]. -This graph keeps all haplotype information. - - -2. Haplotype-resolved processed unitig graph without small bubbles [outPrefix.p_utg.gfa]. -Small bubbles might be caused by somatic mutations or noise in data, which are not the real haplotype information. -The size of popped small bubbles should be specified by -.I [-p maxSmallBubbles]. - - -3. Primary assembly contig graph [outPrefix.p_ctg.gfa]. -This graph collapses different haplotypes. - -4. Alternate assembly contig graph [outPrefix.a_ctg.gfa]. -This graph consists of all assemblies that are discarded in primary assembly contig graph. - - -.PP -If you have trio information, hifiasm outputs the following three assembly graphs in GFA format: - -.IP -1. Phased maternal unitig graph [outPrefix.m.r_utg.gfa]. -This graph keeps the phased maternal assembly. - -2. Phased paternal unitig graph [outPrefix.p.r_utg.gfa]. -This graph keeps the phased paternal assembly. - -3. Haplotype-resolved raw unitig graph [outPrefix.r_utg.gfa]. -This graph keeps all haplotype information. +.TP +* +.IR prefix .hap1.p_ctg.gfa: +phased paternal/haplotype1 contig graph. This graph keeps the phased +paternal/haplotype1 assembly. +.TP +* +.IR prefix .hap2.p_ctg.gfa: +phased maternal/haplotype2 contig graph. This graph keeps the phased +maternal/haplotype2 assembly. +.RE .PP For each graph, hifiasm also outputs a simplified version without sequences. These simplified