diff --git a/Assembly.cpp b/Assembly.cpp index 8ea37b1..df337b3 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -297,7 +297,7 @@ void Counting_multiple_thr() if (asm_opt.roundID == 0) { - init_kseq(asm_opt.read_file_name); + init_gz_files(&asm_opt); init_All_reads(&R_INF); init_R_buffer(asm_opt.thread_num); pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert); @@ -336,7 +336,8 @@ void Counting_multiple_thr() if (asm_opt.roundID == 0) { pthread_join(inputReadsHandle, NULL); - destory_kseq(); + ///destory_kseq(); + destory_gz_files(); } fprintf(stderr, "Counting has been completed.\n"); @@ -372,12 +373,9 @@ void Build_hash_table_multiple_thr() if (asm_opt.roundID == 0) { - init_kseq(asm_opt.read_file_name); - + init_gz_files(&asm_opt); clear_R_buffer(); - malloc_All_reads(&R_INF); - pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert); } @@ -422,7 +420,8 @@ void Build_hash_table_multiple_thr() if (asm_opt.roundID == 0) { - destory_kseq(); + ///destory_kseq(); + destory_gz_files(); destory_R_buffer(); } @@ -1234,7 +1233,10 @@ void Output_corrected_reads() long long i; UC_Read g_read; init_UC_Read(&g_read); - FILE* output_file = fopen(asm_opt.output_file_name, "w"); + char* gfa_name = (char*)malloc(strlen(asm_opt.output_file_name)+35); + sprintf(gfa_name, "%s.ec.fa", asm_opt.output_file_name); + FILE* output_file = fopen(gfa_name, "w"); + free(gfa_name); for (i = 0; i < (long long)R_INF.total_reads; i++) { @@ -1325,18 +1327,18 @@ void Overlap_calculate_multipe_thr() } -int load_pre_cauculated_index() -{ - if(load_Total_Pos_Table(&PCB, asm_opt.read_file_name) && load_All_reads(&R_INF, asm_opt.read_file_name)) - { - return 1; - } - else - { - return 0; - } +// int load_pre_cauculated_index() +// { +// if(load_Total_Pos_Table(&PCB, asm_opt.read_file_name) && load_All_reads(&R_INF, asm_opt.read_file_name)) +// { +// return 1; +// } +// else +// { +// return 0; +// } -} +// } void update_overlaps(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf, diff --git a/Assembly.h b/Assembly.h index 3469c7b..21f1cf8 100644 --- a/Assembly.h +++ b/Assembly.h @@ -10,7 +10,7 @@ void Counting_multiple_thr(); void Build_hash_table_multiple_thr(); -int load_pre_cauculated_index(); +///int load_pre_cauculated_index(); void Overlap_calculate_multipe_thr(); void Correct_Reads(int last_round); #endif diff --git a/CommandLines.cpp b/CommandLines.cpp index 361e2f4..1456f73 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -1,10 +1,12 @@ #include "CommandLines.h" +#include #include #include #include "ketopt.h" #include #define VERSION "0.0.0.1" +#define DEFAULT_OUTPUT "hifiasm.asm" hifiasm_opt_t asm_opt; @@ -17,38 +19,39 @@ double Get_T(void) void Print_H(hifiasm_opt_t* asm_opt) { - fprintf(stderr, "Usage: hifiasm [options] -q -o \n"); + 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, " -k FILE output assembly (in gfa format) and corrected reads (in fasta format)\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, " -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); fprintf(stderr, " -k INT k-mer length [%d] (must be < 64)\n", asm_opt->k_mer_length); - fprintf(stderr, " -w write all overlaps to disk, can accelerate assembly next time\n"); - fprintf(stderr, " -l load all overlaps from disk, can avoid overlap calculation\n"); + fprintf(stderr, " -w write all overlaps to disk, can accelerate assembly next time [%d]\n", asm_opt->write_index_to_disk); + fprintf(stderr, " -l load all overlaps from disk, can avoid overlap calculation [%d]\n", asm_opt->load_index_from_disk); fprintf(stderr, " -z INT length of adapters that should be removed [%d]\n", asm_opt->adapterLen); fprintf(stderr, " -p INT size of popped bubbles [%lld]\n", asm_opt->pop_bubble_size); 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); fprintf(stderr, " -v show version number\n"); fprintf(stderr, " -h show help information\n"); - fprintf(stderr, "Example: ./hifiasm -w -l -q NA12878.fq.gz -o NA12878.asm -k 40 -t 32 -r 2\n"); + fprintf(stderr, "Example: ./hifiasm -w -l -o NA12878.asm -k 40 -t 32 -r 2 NA12878.fq.gz\n"); } void init_opt(hifiasm_opt_t* asm_opt) { - asm_opt->read_file_name = NULL; - asm_opt->output_file_name = NULL; + asm_opt->num_reads = 0; + asm_opt->read_file_names = NULL; + asm_opt->output_file_name = (char*)(DEFAULT_OUTPUT); asm_opt->required_read_name = NULL; asm_opt->thread_num = 1; asm_opt->k_mer_length = 40; asm_opt->k_mer_min_freq = 3; asm_opt->k_mer_max_freq = 66; - asm_opt->load_index_from_disk = 0; - asm_opt->write_index_to_disk = 0; + asm_opt->load_index_from_disk = 1; + asm_opt->write_index_to_disk = 1; asm_opt->number_of_round = 2; asm_opt->adapterLen = 0; asm_opt->clean_round = 4; @@ -58,6 +61,14 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->max_drop_rate = 0.8; } +void destory_opt(hifiasm_opt_t* asm_opt) +{ + if(asm_opt->read_file_names != NULL) + { + free(asm_opt->read_file_names); + } +} + void clear_opt(hifiasm_opt_t* asm_opt, int last_round) { asm_opt->complete_threads = 0; @@ -70,7 +81,7 @@ void clear_opt(hifiasm_opt_t* asm_opt, int last_round) int check_option(hifiasm_opt_t* asm_opt) { - if(asm_opt->read_file_name == NULL) + if(asm_opt->read_file_names == NULL || asm_opt->num_reads == 0) { fprintf(stderr, "[ERROR] missing input: please specify a read file\n"); return 0; @@ -141,7 +152,7 @@ int check_option(hifiasm_opt_t* asm_opt) } - // fprintf(stderr, "input file: %s\n", asm_opt->read_file_name); + // 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); // fprintf(stderr, "number of rounds for correction: %d\n", asm_opt->number_of_round); @@ -155,6 +166,32 @@ int check_option(hifiasm_opt_t* asm_opt) return 1; } +void get_queries(int argc, char *argv[], ketopt_t* opt, hifiasm_opt_t* asm_opt) +{ + if(opt->ind == argc) + { + return; + } + + asm_opt->num_reads = argc - opt->ind; + asm_opt->read_file_names = (char**)malloc(sizeof(char*)*asm_opt->num_reads); + + long long i; + gzFile dfp; + for (i = 0; i < asm_opt->num_reads; i++) + { + asm_opt->read_file_names[i] = argv[i + opt->ind]; + dfp = gzopen(asm_opt->read_file_names[i], "r"); + if (dfp == 0) + { + fprintf(stderr, "[ERROR] Cannot find the input read file: %s\n", + asm_opt->read_file_names[i]); + exit(0); + } + gzclose(dfp); + } +} + int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) { @@ -162,7 +199,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) int c; - while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:q:k:lwm:n:r:a:b:z:x:y:p:", 0)) >= 0) { + while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:lwm:n:r:a:b:z:x:y:p:", 0)) >= 0) { if (c == 'h') { Print_H(asm_opt); @@ -175,7 +212,6 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) } else if (c == 't') asm_opt->thread_num = atoi(opt.arg); else if (c == 'o') asm_opt->output_file_name = opt.arg; - else if (c == 'q') asm_opt->read_file_name = opt.arg; else if (c == 'n') asm_opt->k_mer_min_freq = atoi(opt.arg); else if (c == 'm') asm_opt->k_mer_max_freq = atoi(opt.arg); else if (c == 'r') asm_opt->number_of_round = atoi(opt.arg); @@ -206,6 +242,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) return 0; } + get_queries(argc, argv, &opt, asm_opt); return check_option(asm_opt); } \ No newline at end of file diff --git a/CommandLines.h b/CommandLines.h index f451c7f..52c591a 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -6,7 +6,8 @@ #define VERBOSE 0 typedef struct { - char* read_file_name; + int num_reads; + char** read_file_names; char* output_file_name; char* required_read_name; int thread_num; @@ -32,6 +33,7 @@ typedef struct { 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 last_round); int CommandLine_process (int argc, char *argv[], hifiasm_opt_t* asm_opt); double Get_T(void); diff --git a/Overlaps.cpp b/Overlaps.cpp index 1e40eda..3cdac32 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -1074,15 +1074,12 @@ long long n_read, uint64_t* readLen, ma_sub_t* coverage_cut, float shift_rate) if(max_left.e > max_right.s && (max_left.e - max_right.s >= rLen * shift_rate)) { - ///coverage_cut[i].c = 0; continue; } ///simple chimeric reads if(max_left.e <= max_right.s) { - ///coverage_cut[i].c = 2; - coverage_cut[i].c = 1; coverage_cut[i].del = 1; paf[i].length = 0; n_simple_remove++; @@ -1106,7 +1103,6 @@ long long n_read, uint64_t* readLen, ma_sub_t* coverage_cut, float shift_rate) cov += intersection_check(&(rev_paf[i]), rLen, interval_s, interval_e); if(interval_e - interval_s < WINDOW && cov <= 2) { - coverage_cut[i].c = 1; coverage_cut[i].del = 1; paf[i].length = 0; n_complex_remove_real++; @@ -1119,7 +1115,6 @@ long long n_read, uint64_t* readLen, ma_sub_t* coverage_cut, float shift_rate) /**|| intersection_check_by_base(&(rev_paf[i]), rLen, interval_s, interval_e, b_q.a, b_t.a)**/) { - coverage_cut[i].c = 1; coverage_cut[i].del = 1; paf[i].length = 0; n_complex_remove_real++; @@ -4215,14 +4210,16 @@ int asg_cut_tip(asg_t *g, int max_ext) ///max_ext is 4 -int debug_asg_cut_tip(asg_t *g, int max_ext) +int asg_cut_tip_primary(asg_t *g, int max_ext) { + double startTime = Get_T(); + asg64_v a = {0,0,0}; uint32_t n_vtx = g->n_seq * 2, v, i, cnt = 0; for (v = 0; v < n_vtx; ++v) { //if this seq has been deleted - if (g->seq[v>>1].del) continue; + if (g->seq[v>>1].del || g->seq[v>>1].c) continue; ///check if the another direction of v has no overlaps ///if the self direction of v has no overlaps, we don't have the overlaps of them ///here is check if the reverse direction of v @@ -4250,21 +4247,32 @@ int debug_asg_cut_tip(asg_t *g, int max_ext) * | * ----->n(5) **/ - for (i = 0; i < a.n; ++i) + for (i = 0; i < a.n; ++i) { - asg_seq_del(g, (uint32_t)a.a[i]>>1); - fprintf(stderr, "removed node: %u\n", (uint32_t)a.a[i]>>1); + g->seq[((uint32_t)a.a[i]>>1)].c = 1; } - + + for (i = 0; i < a.n; ++i) + { + asg_seq_drop(g, (uint32_t)a.a[i]>>1); + } + ++cnt; } free(a.a); if (cnt > 0) asg_cleanup(g); - fprintf(stderr, "[M::%s] cut %d tips\n", __func__, cnt); + if(VERBOSE >= 1) + { + fprintf(stderr, "[M::%s] cut %d tips\n", __func__, cnt); + fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime); + } return cnt; } + + + // delete short arcs ///for best graph? int asg_arc_del_short(asg_t *g, float drop_ratio) @@ -4969,7 +4977,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen) int asg_arc_del_short_false_link(asg_t *g, float drop_ratio, float o_drop_ratio, int max_dist, -ma_hit_t_alloc* reverse_sources, long long miniedgeLen) +ma_hit_t_alloc* reverse_sources, long long miniedgeLen, int is_drop) { double startTime = Get_T(); @@ -4994,6 +5002,8 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen) for (v = 0; v < n_vtx; ++v) { + if(is_drop && g->seq[v>>1].c) continue; + if(g->seq_vis[v] == 0) { asg_arc_t *av = asg_arc_a(g, v); @@ -6311,6 +6321,16 @@ int asg_arc_del_false_node(asg_t *g, int max_ext) } +void check_node_lable(asg_t *g) +{ + uint32_t v, n_vtx = g->n_seq * 2; + for (v = 0; v < n_vtx; ++v) + { + if(g->seq[v>>1].c != 0) fprintf(stderr, "error\n"); + } +} + + #define arc_cnt(g, v) ((uint32_t)(g)->idx[(v)]) #define arc_first(g, v) ((g)->arc[(g)->idx[(v)]>>32]) @@ -6325,6 +6345,7 @@ ma_ug_t *ma_ug_gen(asg_t *g) ug = (ma_ug_t*)calloc(1, sizeof(ma_ug_t)); ug->g = asg_init(); + ///each node has two directions mark = (int32_t*)calloc(n_vtx, 4); q = kdq_init(uint64_t); @@ -6332,7 +6353,7 @@ ma_ug_t *ma_ug_gen(asg_t *g) uint32_t w, x, l, start, end, len; ma_utg_t *p; ///what's the usage of mark array - ///mark array is used to select another direction of node + ///mark array is used to mark if this node has already been included in a contig if (g->seq[v>>1].del || arc_cnt(g, v) == 0 || mark[v]) continue; mark[v] = 1; q->count = 0, start = v, end = v^1, len = 0; @@ -6400,7 +6421,7 @@ add_unitig: //ug saves all unitigs for (v = 0; v < n_vtx; ++v) mark[v] = -1; - + //mark all start nodes and end nodes of all unitigs for (i = 0; i < ug->u.n; ++i) { if (ug->u.a[i].circ) continue; mark[ug->u.a[i].start] = i<<1 | 0; @@ -6422,6 +6443,138 @@ add_unitig: **/ ///to connect two unitigs, we need to connect the end of unitig x to the start of unitig y ///so we need to ^1 to get the reverse direction of (x's end)? + ///>=0 means this node is a start/end node of an unitig + ///means this node is a intersaction node + if (mark[p->ul>>32^1] >= 0 && mark[p->v] >= 0) { + asg_arc_t *q; + uint32_t u = mark[p->ul>>32^1]^1; + int l = ug->u.a[u>>1].len - p->ol; + if (l < 0) l = 1; + q = asg_arc_pushp(ug->g); + q->ol = p->ol, q->del = 0; + q->ul = (uint64_t)u<<32 | l; + q->v = mark[p->v]; + } + } + for (i = 0; i < ug->u.n; ++i) + asg_seq_set(ug->g, i, ug->u.a[i].len, 0); + asg_cleanup(ug->g); + free(mark); + return ug; +} + + + +ma_ug_t *ma_ug_gen_primary(asg_t *g, uint8_t flag) +{ + int32_t *mark; + uint32_t i, v, n_vtx = g->n_seq * 2; + ///is a queue + kdq_t(uint64_t) *q; + ma_ug_t *ug; + + ug = (ma_ug_t*)calloc(1, sizeof(ma_ug_t)); + ug->g = asg_init(); + ///each node has two directions + mark = (int32_t*)calloc(n_vtx, 4); + + q = kdq_init(uint64_t); + for (v = 0; v < n_vtx; ++v) { + uint32_t w, x, l, start, end, len; + ma_utg_t *p; + ///what's the usage of mark array + ///mark array is used to mark if this node has already been included in a contig + if (g->seq[v>>1].del || arc_cnt(g, v) == 0 || mark[v] || g->seq[v>>1].c != flag) continue; + mark[v] = 1; + q->count = 0, start = v, end = v^1, len = 0; + // forward + w = v; + while (1) { + /** + * w----->x + * w<-----x + * that means the only suffix of w is x, and the only prefix of x is w + **/ + if (arc_cnt(g, w) != 1) break; + x = arc_first(g, w).v; // w->x + if (arc_cnt(g, x^1) != 1) break; + + /** + * another direction of w would be marked as used (since w has been used) + **/ + mark[x] = mark[w^1] = 1; + ///l is the edge length, instead of overlap length + ///note: edge length is different with overlap length + l = asg_arc_len(arc_first(g, w)); + kdq_push(uint64_t, q, (uint64_t)w<<32 | l); + end = x^1, len += l; + w = x; + if (x == v) break; + } + if (start != (end^1) || kdq_size(q) == 0) { // linear unitig + ///length of seq, instead of edge + l = g->seq[end>>1].len; + kdq_push(uint64_t, q, (uint64_t)(end^1)<<32 | l); + len += l; + } else { // circular unitig + start = end = UINT32_MAX; + goto add_unitig; // then it is not necessary to do the backward + } + // backward + x = v; + while (1) { // similar to forward but not the same + if (arc_cnt(g, x^1) != 1) break; + w = arc_first(g, x^1).v ^ 1; // w->x + if (arc_cnt(g, w) != 1) break; + mark[x] = mark[w^1] = 1; + l = asg_arc_len(arc_first(g, w)); + ///w is the seq id + direction, l is the length of edge + ///push element to the front of a queue + kdq_unshift(uint64_t, q, (uint64_t)w<<32 | l); + start = w, len += l; + x = w; + } +add_unitig: + if (start != UINT32_MAX) mark[start] = mark[end] = 1; + kv_pushp(ma_utg_t, ug->u, &p); + p->s = 0, p->start = start, p->end = end, p->len = len, p->n = kdq_size(q), p->circ = (start == UINT32_MAX); + p->m = p->n; + kv_roundup32(p->m); + p->a = (uint64_t*)malloc(8 * p->m); + //all elements are saved here + for (i = 0; i < kdq_size(q); ++i) + p->a[i] = kdq_at(q, i); + } + kdq_destroy(uint64_t, q); + + // add arcs between unitigs; reusing mark for a different purpose + //ug saves all unitigs + for (v = 0; v < n_vtx; ++v) mark[v] = -1; + + //mark all start nodes and end nodes of all unitigs + for (i = 0; i < ug->u.n; ++i) { + if (ug->u.a[i].circ) continue; + mark[ug->u.a[i].start] = i<<1 | 0; + mark[ug->u.a[i].end] = i<<1 | 1; + } + + //scan all edges + for (i = 0; i < g->n_arc; ++i) { + asg_arc_t *p = &g->arc[i]; + if (p->del) continue; + /** + p->ul: |____________31__________|__________1___________|______________32_____________| + qns direction of overlap length of this node (not overlap length) + (based on query) + p->v : |___________31___________|__________1___________| + tns reverse direction of overlap + (based on target) + p->ol: overlap length + **/ + ///to connect two unitigs, we need to connect the end of unitig x to the start of unitig y + ///so we need to ^1 to get the reverse direction of (x's end)? + ///>=0 means this node is a start/end node of an unitig + ///means this node is a intersaction node if (mark[p->ul>>32^1] >= 0 && mark[p->v] >= 0) { asg_arc_t *q; uint32_t u = mark[p->ul>>32^1]^1; @@ -6864,6 +7017,81 @@ int asg_arc_cut_long_tip(asg_t *g, float drop_ratio) return n_reduced; } + +int asg_arc_cut_long_tip_primary(asg_t *g, float drop_ratio) +{ + double startTime = Get_T(); + ///the reason is that each read has two direction (query->target, target->query) + uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex; + long long ll, v_maxLen; + + buf_t b; + memset(&b, 0, sizeof(buf_t)); + + for (v = 0; v < n_vtx; ++v) + { + uint32_t i, n_arc = 0, nv = asg_arc_n(g, v); + asg_arc_t *av = asg_arc_a(g, v); + ///some node could be deleted + if (nv < 2 || g->seq[v>>1].del || g->seq[v>>1].c) continue; + n_arc = get_real_length(g, v, NULL); + if (n_arc < 2) continue; + + v_maxLen = -1; + + for (i = 0, n_arc = 0; i < nv; i++) + { + if (!av[i].del) + { + detect_single_path_with_dels(g, av[i].v, &convex, &ll, NULL); + if(v_maxLen < ll) + { + v_maxLen = ll; + } + } + } + + for (i = 0, n_arc = 0; i < nv; i++) + { + if (!av[i].del) + { + b.b.n = 0; + if(detect_single_path_with_dels(g, av[i].v, &convex, &ll, &b) == END_TIPS) + { + if(v_maxLen*drop_ratio > ll) + { + n_reduced++; + uint64_t k; + + for (k = 0; k < b.b.n; k++) + { + g->seq[b.b.a[k]].c = 1; + } + + for (k = 0; k < b.b.n; k++) + { + asg_seq_drop(g, b.b.a[k]); + } + } + } + } + } + } + + + asg_cleanup(g); + asg_symm(g); + + if(VERBOSE >= 1) + { + fprintf(stderr, "[M::%s] removed %d long tips\n", + __func__, n_reduced); + fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime); + } + + return n_reduced; +} + uint32_t detect_single_path_with_dels_contigLen(asg_t *g, uint32_t begNode, uint32_t* endNode, long long* baseLen, buf_t* b) { @@ -7029,6 +7257,104 @@ int asg_arc_cut_long_equal_tips(asg_t *g, ma_hit_t_alloc* reverse_sources, long return n_reduced; } + + +int asg_arc_cut_long_equal_tips_assembly(asg_t *g, ma_hit_t_alloc* reverse_sources, long long miniedgeLen) +{ + double startTime = Get_T(); + ///the reason is that each read has two direction (query->target, target->query) + uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, flag; + long long ll, base_maxLen, base_maxLen_i; + + buf_t b; + memset(&b, 0, sizeof(buf_t)); + + for (v = 0; v < n_vtx; ++v) + { + uint32_t i, n_arc = 0, nv = asg_arc_n(g, v), n_tips; + asg_arc_t *av = asg_arc_a(g, v); + ///some node could be deleted + if (nv < 2 || g->seq[v>>1].del || g->seq[v>>1].c) continue; + n_arc = get_real_length(g, v, NULL); + if (n_arc < 2) continue; + + base_maxLen = -1; + base_maxLen_i = -1; + n_tips = 0; + + for (i = 0; i < nv; i++) + { + if (!av[i].del) + { + flag = detect_single_path_with_dels_contigLen(g, av[i].v, &convex, &ll, NULL); + + if(base_maxLen < ll) + { + base_maxLen = ll; + base_maxLen_i = i; + } + + if(flag == END_TIPS) + { + n_tips++; + } + } + } + + ///at least one tip + if(n_tips > 0) + { + for (i = 0; i < nv; i++) + { + if(i == base_maxLen_i) continue; + if (!av[i].del) + { + b.b.n = 0; + if(detect_single_path_with_dels_contigLen(g, av[i].v, &convex, &ll, &b) + != END_TIPS) + { + continue; + } + //we can only cut tips + if(check_if_diploid(av[base_maxLen_i].v, av[i].v, g, + reverse_sources, miniedgeLen)==1) + { + n_reduced++; + uint64_t k; + // for (k = 0; k < b.b.n; k++) + // { + // asg_seq_del(g, b.b.a[k]); + // } + + for (k = 0; k < b.b.n; k++) + { + g->seq[b.b.a[k]].c = 1; + } + for (k = 0; k < b.b.n; k++) + { + asg_seq_drop(g, b.b.a[k]); + } + } + } + } + } + } + + + asg_cleanup(g); + asg_symm(g); + + if(VERBOSE >= 1) + { + fprintf(stderr, "[M::%s] removed %d long tips\n", + __func__, n_reduced); + fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime); + } + + + return n_reduced; +} + void output_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, long long n_read) { ma_ug_t *ug = NULL; @@ -7037,19 +7363,14 @@ void output_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_na fprintf(stderr, "Writing unitig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+25); - sprintf(gfa_name, "%s.gfa", output_file_name); + sprintf(gfa_name, "%s.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.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); free(gfa_name); ma_ug_destroy(ug); @@ -7406,9 +7727,166 @@ int asg_pop_bubble(asg_t *g, int max_dist) } +// in a resolved bubble, mark unused vertices and arcs as "reduced" +static void asg_bub_backtrack_primary(asg_t *g, uint32_t v0, buf_t *b) +{ + uint32_t i, v; + + ///assert(b->S.n == 1); + ///first remove all nodes in this bubble + for (i = 0; i < b->b.n; ++i) + { + g->seq[b->b.a[i]>>1].c = 1; + } + + ///v is the sink of this bubble + v = b->S.a[0]; + ///recover node + do { + uint32_t u = b->a[v].p; // u->v + g->seq[v>>1].c = 0; + v = u; + } while (v != v0); + + for (i = 0; i < b->b.n; ++i) + { + v = b->b.a[i]; + ///if v is not at primary + if(g->seq[v>>1].c) + { + asg_seq_drop(g, v>>1); + } + } +} + +// pop bubbles from vertex v0; the graph MJUST BE symmetric: if u->v present, v'->u' must be present as well +static uint64_t asg_bub_pop1_primary(asg_t *g, uint32_t v0, int max_dist, buf_t *b) +{ + uint32_t i, n_pending = 0; + uint64_t n_pop = 0; + ///if this node has been deleted + if (g->seq[v0>>1].del || g->seq[v0>>1].c) return 0; // already deleted + ///asg_arc_n(n0) + if ((uint32_t)g->idx[v0] < 2) return 0; // no bubbles + ///S saves nodes with all incoming edges visited + b->S.n = b->T.n = b->b.n = b->e.n = 0; + ///for each node, b->a saves all related information + b->a[v0].c = b->a[v0].d = 0; + ///b->S is the nodes with all incoming edges visited + kv_push(uint32_t, b->S, v0); + + do { + ///v is a node that all incoming edges have been visited + ///d is the distance from v0 to v + uint32_t v = kv_pop(b->S), d = b->a[v].d, c = b->a[v].c; + uint32_t nv = asg_arc_n(g, v); + asg_arc_t *av = asg_arc_a(g, v); + ///why we have this assert? + ///assert(nv > 0); + + ///all out-edges of v + for (i = 0; i < nv; ++i) { // loop through v's neighbors + /** + p->ul: |____________31__________|__________1___________|______________32_____________| + qn direction of overlap length of this node (not overlap length) + (in the view of query) + p->v : |___________31___________|__________1___________| + tn reverse direction of overlap + (in the view of target) + p->ol: overlap length + **/ + uint32_t w = av[i].v, l = (uint32_t)av[i].ul; // v->w with length l + binfo_t *t = &b->a[w]; + ///that means there is a circle, directly terminate the whole bubble poping + ///if (w == v0) goto pop_reset; + if ((w>>1) == (v0>>1)) goto pop_reset; + ///if this edge has been deleted + if (av[i].del) continue; + + ///push the edge + ///high 32-bit of g->idx[v] is the start point of v's edges + //so here is the point of this specfic edge + kv_push(uint32_t, b->e, (g->idx[v]>>32) + i); + + ///find a too far path? directly terminate the whole bubble poping + if (d + l > (uint32_t)max_dist) break; // too far + + ///if this node + if (t->s == 0) { // this vertex has never been visited + kv_push(uint32_t, b->b, w); // save it for revert + ///t->p is the parent node of + ///t->s = 1 means w has been visited + ///d is len(v0->v), l is len(v->w), so t->d is len(v0->w) + t->p = v, t->s = 1, t->d = d + l; + ///incoming edges of w + t->r = count_out(g, w^1); + ++n_pending; + } else { // visited before + ///c is the weight (is very likely the number of node in this edge) of the parent node + ///select the longest edge (longest meams most reads/longest edge) + if (c + 1 > t->c || (c + 1 == t->c && d + l > t->d)) t->p = v; + if (c + 1 > t->c) t->c = c + 1; + ///update len(v0->w) + ///node: t->d is not the length from this node's parent + ///it is the shortest edge + if (d + l < t->d) t->d = d + l; // update dist + } + ///assert(t->r > 0); + //if all incoming edges of w have visited + //push it to b->S + if (--(t->r) == 0) { + uint32_t x = asg_arc_n(g, w); + if (x) kv_push(uint32_t, b->S, w); + ///else kv_push(uint32_t, b->T, w); // a tip + else goto pop_reset; + --n_pending; + } + } + ///if i < nv, that means (d + l > max_dist) + if (i < nv || b->S.n == 0) goto pop_reset; + } while (b->S.n > 1 || n_pending); + asg_bub_backtrack_primary(g, v0, b); + n_pop = 1; +pop_reset: + for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices + binfo_t *t = &b->a[b->b.a[i]]; + t->s = t->c = t->d = 0; + } + return n_pop; +} +// pop bubbles +int asg_pop_bubble_primary(asg_t *g, int max_dist) +{ + uint32_t v, n_vtx = g->n_seq * 2; + uint64_t n_pop = 0; + buf_t b; + if (!g->is_symm) asg_symm(g); + memset(&b, 0, sizeof(buf_t)); + ///set information for each node + b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); + //traverse all node with two directions + for (v = 0; v < n_vtx; ++v) { + uint32_t i, n_arc = 0, nv = asg_arc_n(g, v); + asg_arc_t *av = asg_arc_a(g, v); + ///some node could be deleted + if (nv < 2 || g->seq[v>>1].del || g->seq[v>>1].c) continue; + ///some edges could be deleted + for (i = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs + if (!av[i].del) ++n_arc; + if (n_arc > 1) + n_pop += asg_bub_pop1_primary(g, v, max_dist, &b); + } + free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); + if (n_pop) asg_cleanup(g); + if(VERBOSE >= 1) + { + fprintf(stderr, "[M::%s] popped %lu bubbles\n", __func__, (unsigned long)n_pop); + } + return n_pop; +} @@ -7771,7 +8249,7 @@ uint32_t detect_single_path_with_dels_by_length -long long asg_arc_del_self_circle_untig(asg_t *g, long long circleLen) +long long asg_arc_del_self_circle_untig(asg_t *g, long long circleLen, int is_drop) { uint32_t v, w, n_vtx = g->n_seq * 2, n_reduced = 0, convex, flag; long long ll; @@ -7783,6 +8261,8 @@ long long asg_arc_del_self_circle_untig(asg_t *g, long long circleLen) asg_arc_t *av = asg_arc_a(g, v); ///some node could be deleted if (g->seq[v>>1].del) continue; + if(is_drop && g->seq[v>>1].c) continue; + n_arc = get_real_length(g, v, NULL); if (n_arc != 1) continue; @@ -7814,12 +8294,6 @@ long long asg_arc_del_self_circle_untig(asg_t *g, long long circleLen) { if ((!aw[k].del) && (aw[k].v == (convex^1))) { - // fprintf(stderr, "v: %u, w: %u, aw[%d].v: %u, convex: %u, ll: %d\n", - // v>>1, w>>1, k, aw[k].v>>1, convex>>1, ll); - - // fprintf(stderr, "*aw[k].v: %u, *convex: %u\n\n", - // aw[k].v, convex); - aw[k].del = 1; asg_arc_del(g, aw[k].v^1, aw[k].ul>>32^1, 1); n_reduced++; @@ -7845,6 +8319,8 @@ long long asg_arc_del_self_circle_untig(asg_t *g, long long circleLen) } + + void output_unitig_graph_without_small_bubbles(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, long long n_read, long long bubble_dist, long long tipsLen) { @@ -7874,6 +8350,35 @@ char* output_file_name, long long n_read, long long bubble_dist, long long tipsL } +void output_unitig_graph_without_small_bubbles_primary(asg_t *sg, ma_sub_t* coverage_cut, +char* output_file_name, long long n_read, long long bubble_dist, long long tipsLen) +{ + asg_cut_tip_primary(sg, tipsLen); + asg_pop_bubble_primary(sg, bubble_dist); + asg_cut_tip_primary(sg, tipsLen); + + ma_ug_t *ug = NULL; + 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"); + char* gfa_name = (char*)malloc(strlen(output_file_name)+35); + sprintf(gfa_name, "%s.wsb.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); + + free(gfa_name); + ma_ug_destroy(ug); +} + + + void output_contig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, long long n_read, long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long circleLen, ma_hit_t_alloc* reverse_sources, long long miniedgeLen) { @@ -7889,7 +8394,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen) pre_cons = sg->n_seq + sg->n_arc; n_ac = 0; n_ac += asg_pop_bubble(sg, bubble_dist); - n_ac += asg_arc_del_self_circle_untig(sg, circleLen); + n_ac += asg_arc_del_self_circle_untig(sg, circleLen, 0); n_ac += asg_arc_cut_long_tip(sg, tip_drop_ratio); n_ac += asg_arc_cut_long_equal_tips(sg, reverse_sources, 2); @@ -7897,7 +8402,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen) } asg_arc_identify_simple_bubbles_multi(sg, 1); - asg_arc_del_short_false_link(sg, 0.6, 0.85, bubble_dist, reverse_sources, MAX_SHORT_TIPS); + asg_arc_del_short_false_link(sg, 0.6, 0.85, bubble_dist, reverse_sources, MAX_SHORT_TIPS, 0); ///asg_arc_del_self_circle_untig(sg, circleLen); @@ -7925,6 +8430,92 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen) ma_ug_destroy(ug); } +long long get_graph_statistic(asg_t *g) +{ + long long num_arc = 0; + uint32_t n_vtx = g->n_seq * 2, v; + + for (v = 0; v < n_vtx; ++v) + { + if (g->seq[v>>1].del || g->seq[v>>1].c) continue; + num_arc += asg_arc_n(g, v); + } + + return num_arc; +} + +void output_contig_graph_primary(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, long long n_read, long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long circleLen, +ma_hit_t_alloc* reverse_sources, long long miniedgeLen) +{ + asg_cut_tip_primary(sg, tipsLen); + long long n_ac = 1; + long long pre_cons = get_graph_statistic(sg); + long long cur_cons = 0; + ///while(n_ac > 0) + while(pre_cons != cur_cons) + { + pre_cons = get_graph_statistic(sg); + n_ac = 0; + n_ac += asg_pop_bubble_primary(sg, bubble_dist); + ///we don't need a special function here since it just removes edges instead of nodes + n_ac += asg_arc_del_self_circle_untig(sg, circleLen, 1); + n_ac += asg_arc_cut_long_tip_primary(sg, tip_drop_ratio); + n_ac += asg_arc_cut_long_equal_tips_assembly(sg, reverse_sources, 2); + + cur_cons = get_graph_statistic(sg); + } + + asg_arc_identify_simple_bubbles_multi(sg, 1); + ///we don't need a special function here since it just removes edges instead of nodes + asg_arc_del_short_false_link(sg, 0.6, 0.85, bubble_dist, reverse_sources, MAX_SHORT_TIPS, 1); + + + ma_ug_t *ug = NULL; + 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"); + char* gfa_name = (char*)malloc(strlen(output_file_name)+35); + sprintf(gfa_name, "%s.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); + + free(gfa_name); + ma_ug_destroy(ug); +} + + +void output_contig_graph_alternative(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, long long n_read) +{ + ma_ug_t *ug = NULL; + 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"); + char* gfa_name = (char*)malloc(strlen(output_file_name)+35); + sprintf(gfa_name, "%s.alter.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); + + free(gfa_name); + ma_ug_destroy(ug); +} int output_tips(asg_t *g, const All_reads *RNF) { @@ -8259,10 +8850,18 @@ long long bubble_dist, int read_graph, int write) asg_t *sg = NULL; sg = ma_sg_gen(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length); asg_arc_del_trans(sg, GAP_FUZZ); - char* unlean_name = (char*)malloc(strlen(output_file_name)+25); - sprintf(unlean_name, "%s.unclean", output_file_name); - output_read_graph(sg, coverage_cut, unlean_name, n_read); - free(unlean_name); + + + if(VERBOSE >= 1) + { + char* unlean_name = (char*)malloc(strlen(output_file_name)+25); + sprintf(unlean_name, "%s.unclean", output_file_name); + output_read_graph(sg, coverage_cut, unlean_name, n_read); + free(unlean_name); + } + + + asg_cut_tip(sg, MAX_SHORT_TIPS); // debug_info_of_specfic_node("m64016_190918_162737/72220752/ccs", sg, "cut_tip"); @@ -8342,7 +8941,7 @@ long long bubble_dist, int read_graph, int write) asg_cut_tip(sg, MAX_SHORT_TIPS); asg_arc_identify_simple_bubbles_multi(sg, 1); - asg_arc_del_short_false_link(sg, 0.6, 0.85, bubble_dist, reverse_sources, MAX_SHORT_TIPS); + asg_arc_del_short_false_link(sg, 0.6, 0.85, bubble_dist, reverse_sources, MAX_SHORT_TIPS, 0); asg_arc_identify_simple_bubbles_multi(sg, 1); asg_arc_del_complex_false_link(sg, 0.6, 0.85, bubble_dist, reverse_sources, MAX_SHORT_TIPS); @@ -8455,20 +9054,22 @@ long long bubble_dist, int read_graph, int write) ///out: ///output_tips(sg, &R_INF); + ///check_node_lable(sg); output_unitig_graph(sg, coverage_cut, output_file_name, n_read); - output_read_graph(sg, coverage_cut, output_file_name, n_read); - /****************************may have bugs********************************/ - output_unitig_graph_without_small_bubbles(sg, coverage_cut, output_file_name, n_read, + if(VERBOSE >= 1) + { + output_read_graph(sg, coverage_cut, output_file_name, n_read); + } + + output_unitig_graph_without_small_bubbles_primary(sg, coverage_cut, output_file_name, n_read, pop_bubble_size, MAX_SHORT_TIPS); - /****************************may have bugs********************************/ - - output_contig_graph(sg, coverage_cut, output_file_name, n_read, 10000000, MAX_SHORT_TIPS, 0.1, 20, + output_contig_graph_primary(sg, coverage_cut, output_file_name, n_read, 10000000, MAX_SHORT_TIPS, 0.1, 20, reverse_sources, MAX_SHORT_TIPS); - ///output_contig_graph(sg, coverage_cut, output_file_name, n_read, 10000000); + output_contig_graph_alternative(sg, coverage_cut, output_file_name, n_read); + - asg_destroy(sg); free(coverage_cut); } diff --git a/Overlaps.h b/Overlaps.h index 167d7e7..41c8f50 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -235,6 +235,32 @@ static inline void asg_seq_del(asg_t *g, uint32_t s) } } +static inline void asg_seq_drop(asg_t *g, uint32_t s) +{ + ///s is not at primary + if(g->seq[s].c) + { + uint32_t k; + for (k = 0; k < 2; ++k) + { + ///two directions of this node + uint32_t i, v = s<<1 | k; + uint32_t nv = asg_arc_n(g, v); + asg_arc_t *av = asg_arc_a(g, v); + for (i = 0; i < nv; ++i) + { + if(av[i].del) continue; + ///if output node is at primary + if(g->seq[(av[i].v>>1)].c == 0) + { + av[i].del = 1; + asg_arc_del(g, av[i].v^1, v^1, 1); + } + } + } + } +} + typedef struct { uint32_t len:31, circ:1; // len: length of the unitig; circ: circular if non-zero diff --git a/POA.h b/POA.h index bb7529b..6621759 100644 --- a/POA.h +++ b/POA.h @@ -410,8 +410,8 @@ inline long long get_bi_Edge(Graph* DAGCon, Node* inNode, Node* outNode, Edge** inline long long get_Edge_Weight(Graph* DAGCon, Node* inNode, Node* outNode) { - Edge* e_forward; - Edge* e_backward; + Edge* e_forward = NULL; + Edge* e_backward = NULL; get_bi_Edge(DAGCon, inNode, outNode, &e_forward, &e_backward); return e_forward->weight; } diff --git a/Process_Read.cpp b/Process_Read.cpp index 751fb70..a7aad95 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -6,8 +6,8 @@ #include -gzFile fp; -kseq_t *seq; +gz_files fps; + R_buffer RDB; static uint64_t total_reads; @@ -700,21 +700,58 @@ void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_l } -void init_kseq(char* file) + + + +void open_file(gz_files* nfps, char* name) { - fp = gzopen(file, "r"); - if (fp == 0) + nfps->fp = gzopen(name, "r"); + if(nfps->fp == 0) { - fprintf(stderr, "[ERROR] Cannot find the input file: %s\n", file); + fprintf(stderr, "[ERROR] Cannot find the input file: %s\n", name); exit(0); - } - seq = kseq_init(fp); + } + nfps->seq = kseq_init(nfps->fp); } -void destory_kseq() + +void close_file(gz_files* nfps) { - kseq_destroy(seq); - gzclose(fp); + kseq_destroy(nfps->seq); + gzclose(nfps->fp); +} + +void init_gz_files(hifiasm_opt_t* asm_opt) +{ + fps.idx = 0; + fps.num_reads = asm_opt->num_reads; + fps.reads = asm_opt->read_file_names; + fps.seq = NULL; + fps.fp = NULL; + if(fps.num_reads > 0) + { + open_file(&fps, fps.reads[fps.idx]); + fps.idx++; + } +} + +void destory_gz_files() +{ + close_file(&fps); +} + +int read_item() +{ + int l = kseq_read(fps.seq); + if(l >= 0 || (l < 0 && fps.idx >= fps.num_reads)) + { + return l; + } + ///l < 0 && fps.idx < fps.num_reads + close_file(&fps); + open_file(&fps, fps.reads[fps.idx]); + fps.idx++; + return read_item(); } @@ -729,14 +766,15 @@ inline void exchage_kstring_t(kstring_t* a, kstring_t* b) int get_read(kseq_t *s, int adapterLen) { int l; - - if ((l = kseq_read(seq)) >= 0) + + ///if ((l = kseq_read(seq)) >= 0) + if ((l = read_item()) >= 0) { - exchage_kstring_t(&seq->comment, &s->comment); - exchage_kstring_t(&seq->name, &s->name); - exchage_kstring_t(&seq->qual, &s->qual); - exchage_kstring_t(&seq->seq, &s->seq); + exchage_kstring_t(&(fps.seq->comment), &s->comment); + exchage_kstring_t(&(fps.seq->name), &s->name); + exchage_kstring_t(&(fps.seq->qual), &s->qual); + exchage_kstring_t(&(fps.seq->seq), &s->seq); if(adapterLen > 0) { diff --git a/Process_Read.h b/Process_Read.h index 5e57c5e..ced5c44 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -41,8 +41,6 @@ extern char rc_Table[5]; #define RC_CHAR(x) rc_Table[seq_nt6_table[(uint8_t)x]] void init_aux_table(); -void init_kseq(char* file); -void destory_kseq(); int get_read(kseq_t *s, int adapterLen); @@ -172,6 +170,15 @@ typedef struct long long RID; } UC_Read; +typedef struct +{ + gzFile fp; + kseq_t *seq; + char** reads; + int num_reads; + int idx; +} gz_files; + void init_R_buffer(int thread_num); void init_All_reads(All_reads* r); void* input_reads_muti_threads(void*); @@ -192,5 +199,8 @@ void destory_R_buffer_block(R_buffer_block* curr_sub_block); void destory_R_buffer(); void clear_R_buffer(); +void init_gz_files(hifiasm_opt_t* asm_opt); +void destory_gz_files(); + #endif diff --git a/README.md b/README.md index 0c2c223..e037bc2 100644 --- a/README.md +++ b/README.md @@ -4,12 +4,18 @@ # Install hifiasm (requiring g++ and zlib) git clone https://github.com/chhylp123/hifiasm.git cd hifiasm && make -# Assembly (corrected reads are at NA12878.asm.fa, assembly graph is at NA12878.asm.fa.gfa) -./hifiasm -w -l -q NA12878.fq.gz -o NA12878.asm.fa -k 40 -t 32 -r 2 +# Assembly +./hifiasm -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 there are two types of output: (1) haplotype-resolved assembly graph in [GFA][gfa] format; (2) haplotype-aware error corrected reads. So far hifiasm is still in early development stage, it will output phased chromosome-level high-quality assembly in the near future. +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 output consists of: +1. haplotype-resolved assembly [unitig][unitig] graph in [GFA][gfa] format ([unitig][unitig] graph). +2. haplotype-resolved assembly [unitig][unitig] graph in [GFA][gfa] format without small bubbles. +3. primary assembly [contig][unitig] graph in [GFA][gfa] format. +4. alternate assembly [contig][unitig] graph in [GFA][gfa] format. +5. haplotype-aware error corrected reads +So far hifiasm is still in early development stage, it will output phased chromosome-level high-quality assembly in the near future. 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: @@ -28,7 +34,7 @@ Hifiasm is a standalone and lightweight assembler, which does not need external For Hifi reads assembly, a typical command line looks like: ```sh -./hifiasm -w -l -q NA12878.fq.gz -o NA12878.asm.fa -k 40 -t 32 -r 2 +./hifiasm -t 32 NA12878.fq.gz ``` where `-q` specifies the input reads and `-o` specifies the output files. In this example, the assembly graph can be found at NA12878.asm.fa.gfa, and the corrected reads can be found at NA12878.asm.fa. `-k`, `-t` and `-r` specify the length of k-mer, the number of CPU threads, and the number of correction rounds, respectively. Note that `-w` means hifiasm will save all overlaps to disk, which can avoid the time-consuming all-to-all overlap calculation next time. For hifiasm with `-l`, if the overlap information has been obtained by `-w` in advance, it is able to load all overlaps from disk and then directly do assembly. diff --git a/main.cpp b/main.cpp index 351dfed..322e633 100644 --- a/main.cpp +++ b/main.cpp @@ -12,6 +12,8 @@ int main(int argc, char *argv[]) if (!CommandLine_process(argc, argv, &asm_opt)) return 1; Correct_Reads(asm_opt.number_of_round); + + destory_opt(&asm_opt); return 1; }