diff --git a/Assembly.cpp b/Assembly.cpp index d2b9134..bb457b5 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -10,6 +10,7 @@ #include "Correct.h" #include "htab.h" #include "kthread.h" +#include "utils.h" void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain); diff --git a/Hash_Table.cpp b/Hash_Table.cpp index b8a3318..d6c7b57 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -4,6 +4,7 @@ #include #include "Hash_Table.h" #include "ksort.h" +#include "utils.h" pthread_mutex_t output_mutex; #define overlap_region_key(a) ((a).y_id) diff --git a/Makefile b/Makefile index 250c1ae..cda0636 100644 --- a/Makefile +++ b/Makefile @@ -4,7 +4,7 @@ CPPFLAGS= INCLUDES= OBJS= CommandLines.o Process_Read.o Assembly.o Hash_Table.o \ POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \ - htab.o hist.o sketch.o anchor.o extract.o sys.o + htab.o hist.o sketch.o anchor.o extract.o overlap2.o sys.o EXE= hifiasm LIBS= -lz -lpthread -lm @@ -34,17 +34,17 @@ depend: Assembly.o: Assembly.h CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h Assembly.o: Hash_Table.h htab.h POA.h Correct.h Levenshtein_distance.h -Assembly.o: kthread.h +Assembly.o: kthread.h utils.h CommandLines.o: CommandLines.h ketopt.h Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h Correct.o: kdq.h CommandLines.h Levenshtein_distance.h POA.h Assembly.h Hash_Table.o: Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h -Hash_Table.o: CommandLines.h ksort.h +Hash_Table.o: CommandLines.h ksort.h utils.h Levenshtein_distance.o: Levenshtein_distance.h Output.o: Output.h CommandLines.h Overlaps.o: Overlaps.h kvec.h kdq.h ksort.h Process_Read.h CommandLines.h Overlaps.o: Hash_Table.h htab.h Correct.h Levenshtein_distance.h POA.h -Overlaps.o: Purge_Dups.h +Overlaps.o: Purge_Dups.h utils.h POA.o: POA.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h POA.o: CommandLines.h Correct.h Levenshtein_distance.h Process_Read.o: Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h @@ -52,16 +52,17 @@ Purge_Dups.o: ksort.h Purge_Dups.h kvec.h kdq.h Overlaps.h Hash_Table.h Purge_Dups.o: htab.h Process_Read.h CommandLines.h Correct.h Purge_Dups.o: Levenshtein_distance.h POA.h kthread.h Trio.o: khashl.h kthread.h kseq.h Process_Read.h Overlaps.h kvec.h kdq.h -Trio.o: CommandLines.h htab.h +Trio.o: CommandLines.h htab.h utils.h anchor.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h -anchor.o: ksort.h Hash_Table.h +anchor.o: ksort.h utils.h Hash_Table.h extract.o: Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h khashl.h -extract.o: kseq.h +extract.o: kseq.h utils.h hist.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h htab.o: kthread.h khashl.h kseq.h ksort.h htab.h Process_Read.h Overlaps.h -htab.o: kvec.h kdq.h CommandLines.h +htab.o: kvec.h kdq.h CommandLines.h utils.h kthread.o: kthread.h main.o: CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h Assembly.h -main.o: Levenshtein_distance.h htab.h +main.o: Levenshtein_distance.h utils.h +overlap2.o: utils.h CommandLines.h Overlaps.h kvec.h kdq.h sketch.o: kvec.h htab.h Process_Read.h Overlaps.h kdq.h CommandLines.h -sys.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h +sys.o: utils.h diff --git a/Overlaps.cpp b/Overlaps.cpp index 6ff6259..b62f08c 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -9,6 +9,7 @@ #include "Hash_Table.h" #include "Correct.h" #include "Purge_Dups.h" +#include "utils.h" uint32_t debug_purge_dup = 0; @@ -76,21 +77,17 @@ void asg_destroy(asg_t *g) { if (g == 0) return; free(g->seq); free(g->idx); free(g->arc); free(g->seq_vis); - - if(g->n_F_seq > 0 && g->F_seq) - { - uint32_t i = 0; - for (i = 0; i < g->n_F_seq; i++) - { - if(g->F_seq[i].a) free(g->F_seq[i].a); - if(g->F_seq[i].s) free(g->F_seq[i].s); - } + free(g->contain); free(g->con_idx); - free(g->F_seq); - } - - - free(g); + if (g->n_F_seq > 0 && g->F_seq) { + uint32_t i = 0; + for (i = 0; i < g->n_F_seq; ++i) { + if (g->F_seq[i].a) free(g->F_seq[i].a); + if (g->F_seq[i].s) free(g->F_seq[i].s); + } + free(g->F_seq); + } + free(g); } void asg_arc_sort(asg_t *g) @@ -132,8 +129,7 @@ void remove_overlaps(ma_hit_t_alloc* source_paf, uint64_t* source_index, long lo } -void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, -uint64_t* source_index, long long listLen) +void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, uint64_t* source_index, long long listLen) { long long i; ma_hit_t ele; @@ -186,7 +182,6 @@ uint64_t *asg_arc_index_core(size_t max_seq, size_t n, const asg_arc_t *a) uint64_t *idx; idx = (uint64_t*)calloc(max_seq * 2, 8); - /** * ul: |____________31__________|__________1___________|______________32_____________| qns direction of overlap length of this node (not overlap length) @@ -196,7 +191,6 @@ uint64_t *asg_arc_index_core(size_t max_seq, size_t n, const asg_arc_t *a) if (i == n || a[i-1].ul>>32 != a[i].ul>>32) idx[a[i-1].ul>>32] = (uint64_t)last<<32 | (i - last), last = i; - return idx; } @@ -276,16 +270,15 @@ void asg_cleanup(asg_t *g) qns direction of overlap length of this node (not overlap length) **/ asg_arc_sort(g); + asg_con_sort(g); g->is_srt = 1; } ///index the overlaps in graph with query id if (g->idx == 0) asg_arc_index(g); + if (g->con_idx == 0 && g->n_con > 0) asg_con_index(g); } - - - // delete multi-arcs /** * remove edges like: v has two out-edges to w @@ -872,51 +865,6 @@ void drop_edges_by_trio(ma_hit_t_alloc* sources, long long num_sources) } - -ma_hit_t* get_specific_overlap_with_del(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, -uint32_t qn, uint32_t tn) -{ - if(coverage_cut[qn].del || coverage_cut[tn].del) return NULL; - ma_hit_t_alloc* x = &(sources[qn]); - uint32_t i; - for (i = 0; i < x->length; i++) - { - if(x->buffer[i].del) continue; - if(coverage_cut[Get_qn(x->buffer[i])].del) continue; - if(coverage_cut[Get_tn(x->buffer[i])].del) continue; - - if(Get_tn(x->buffer[i])==tn - && - Get_qn(x->buffer[i])==qn) - { - return &(x->buffer[i]); - } - } - - return NULL; -} - - - -void delete_single_edge(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, uint32_t qn, uint32_t tn) -{ - ma_hit_t* tmp = get_specific_overlap_with_del(sources, coverage_cut, qn, tn); - if(tmp != NULL) tmp->del = 1; -} - -void delete_all_edges(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, uint32_t qn) -{ - ma_hit_t_alloc* x = &(sources[qn]); - uint32_t i; - for (i = 0; i < x->length; i++) - { - x->buffer[i].del = 1; - delete_single_edge(sources, coverage_cut, Get_tn(x->buffer[i]), Get_qn(x->buffer[i])); - } - coverage_cut[qn].del = 1; -} - - uint32_t get_real_sources_length(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, int max_hang, int min_ovlp, uint32_t query) { @@ -1027,108 +975,6 @@ int max_hang, int min_ovlp, uint32_t qn) return keep_edge; } - -void ma_hit_contained_advance(ma_hit_t_alloc* sources, long long n_read, ma_sub_t *coverage_cut, -R_to_U* ruIndex, int max_hang, int min_ovlp) -{ - ///uint32_t qn_num = 0, no_fully_qn_num = 0, tn_num = 0, no_fully_tn_num = 0; - double startTime = Get_T(); - int32_t r; - long long i, j, m; - asg_arc_t t; - ma_hit_t *h = NULL; - ma_sub_t *sq = NULL; - ma_sub_t *st = NULL; - - - for (i = 0; i < n_read; ++i) - { - if(coverage_cut[i].del) continue; - - for (j = 0; j < (long long)sources[i].length; j++) - { - h = &(sources[i].buffer[j]); - //check the corresponding two reads - sq = &(coverage_cut[Get_qn(*h)]); - st = &(coverage_cut[Get_tn(*h)]); - /****************************may have trio bugs********************************/ - if(sq->del || st->del) continue; - if(h->del) continue; - /****************************may have trio bugs********************************/ - r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); - ///r could not be MA_HT_SHORT_OVLP or MA_HT_INT - if (r == MA_HT_QCONT) - { - h->del = 1; - delete_single_edge(sources, coverage_cut, Get_tn(*h), Get_qn(*h)); - - delete_all_edges(sources, coverage_cut, Get_qn(*h)); - set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0); - - // if(delete_all_edges_carefully(sources, coverage_cut, max_hang, min_ovlp, - // Get_qn(*h))==0) - // { - // set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0); - // } - // sq->del = 1; - // set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0); - } - else if (r == MA_HT_TCONT) - { - h->del = 1; - delete_single_edge(sources, coverage_cut, Get_tn(*h), Get_qn(*h)); - - delete_all_edges(sources, coverage_cut, Get_tn(*h)); - set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0); - - // if(delete_all_edges_carefully(sources, coverage_cut, max_hang, - // min_ovlp, Get_tn(*h)) == 0) - // { - // set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0); - // no_fully_tn_num++; - // } - // st->del = 1; - // set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0); - } - } - } - - transfor_R_to_U(ruIndex); - - - for (i = 0; i < n_read; ++i) - { - m = 0; - for (j = 0; j < (long long)sources[i].length; j++) - { - ma_hit_t *h = &(sources[i].buffer[j]); - if(h->del) continue; - ///both the qn and tn have not been deleted - if(coverage_cut[Get_qn(*h)].del != 1 && coverage_cut[Get_tn(*h)].del != 1) - { - h->del = 0; - m++; - } - else - { - h->del = 1; - } - } - - ///if sources[i].length == 0, that means all overlapped reads with read i are the contained reads - if(m == 0) - { - coverage_cut[i].del = 1; - } - } - - if(VERBOSE >= 1) - { - fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime); - } -} - - void ma_hit_flt(ma_hit_t_alloc* sources, long long n_read, ma_sub_t *coverage_cut, int max_hang, int min_ovlp) { double startTime = Get_T(); @@ -1694,9 +1540,7 @@ void print_overlaps(ma_hit_t_alloc* paf, long long rLen, long long interval_s, l } - -void detect_chimeric_reads(ma_hit_t_alloc* paf, long long n_read, uint64_t* readLen, -ma_sub_t* coverage_cut, float shift_rate) +void detect_chimeric_reads(ma_hit_t_alloc* paf, long long n_read, uint64_t* readLen, ma_sub_t* coverage_cut, float shift_rate) { double startTime = Get_T(); init_aux_table(); @@ -2108,71 +1952,6 @@ static inline int asg_is_single_edge(const asg_t *g, uint32_t v, uint32_t start_ } -asg_t *ma_sg_gen(const ma_hit_t_alloc* sources, long long n_read, const ma_sub_t *coverage_cut, -int max_hang, int min_ovlp) -{ - double startTime = Get_T(); - size_t i, j; - asg_t *g; - ///just calloc - g = asg_init(); - - ///add seq to graph, seq just save the length of each read - for (i = 0; i < (uint64_t)n_read; ++i) - { - ///if a read has been deleted, should we still add them? - asg_seq_set(g, i, coverage_cut[i].e - coverage_cut[i].s, coverage_cut[i].del); - g->seq[i].c = coverage_cut[i].c; - } - - g->seq_vis = (uint8_t*)calloc(g->n_seq*2, sizeof(uint8_t)); - - for (i = 0; i < (uint64_t)n_read; ++i) - { - for (j = 0; j < sources[i].length; j++) - { - int r; - asg_arc_t t, *p; - const ma_hit_t *h = &(sources[i].buffer[j]); - if(h->del) continue; - - //high coverage region [sub[qn].e, sub[qn].s) in query - int ql = coverage_cut[Get_qn(*h)].e - coverage_cut[Get_qn(*h)].s; - //high coverage region [sub[qn].e, sub[qn].s) in target - int tl = coverage_cut[Get_tn(*h)].e - coverage_cut[Get_tn(*h)].s; - r = ma_hit2arc(h, ql, tl, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); - /** - #define MA_HT_INT (-1) - #define MA_HT_QCONT (-2) - #define MA_HT_TCONT (-3) - #define MA_HT_SHORT_OVLP (-4) - the short overlaps and the overlaps with contain reads have already been removed - here we should have overhang - so r should always >= 0 - **/ - if (r >= 0) - { - ///push node? - p = asg_arc_pushp(g); - *p = t; - } - else - { - fprintf(stderr, "error\n"); - } - } - } - - asg_cleanup(g); - g->r_seq = g->n_seq; - - if(VERBOSE >= 1) - { - fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime); - } - - return g; -} @@ -4484,143 +4263,6 @@ int asg_arc_del_cross_bubble(asg_t *g, long long max_dist) return n_reduced; } - - -// transitive reduction; see Myers, 2005 -int asg_arc_del_trans(asg_t *g, int fuzz) -{ - double startTime = Get_T(); - - uint8_t *mark; - ///n_vtx = number of seq * 2 - ///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; - ///at first, all nodes should be set to vacant - mark = (uint8_t*)calloc(n_vtx, 1); - - /**v is the id+direction of a node, - * the high 31-bit is the id, - * and the lowest 1-bit is the direction - * (0 means query-to-target, 1 means target-to-query)**/ - for (v = 0; v < n_vtx; ++v) { - ///nv is the number of overlaps with v(qn+direction) - uint32_t L, i, nv = asg_arc_n(g, v); - ///av is the array of v - asg_arc_t *av = asg_arc_a(g, v); - ///that means in this direction, read v is not overlapped with any other reads - if (nv == 0) continue; // no hits - - ///if the read itself has been removed - if (g->seq[v>>1].del) - { - for (i = 0; i < nv; ++i) av[i].del = 1, ++n_reduced; - continue; - } - - - - /** - ********************************query-to-target overlap**************************** - case 1: u = 0, rev = 0 in the view of target: direction is 1 - query: CCCCCCCCTAATTAAAAT target: TAATTAAAATGGGGGG (use ex-target as query) - |||||||||| <---> |||||||||| - target: TAATTAAAATGGGGGG query: CCCCCCCCTAATTAAAAT (use ex-query as target) - - case 2: u = 0, rev = 1 in the view of target: direction is 0 - query: CCCCCCCCTAATTAAAAT target: CCCCCCATTTTAATTA (use ex-target as query) - |||||||||| <---> |||||||||| - target: TAATTAAAATGGGGGG query: ATTTTAATTAGGGGGGGG (use ex-query as target) - ********************************query-to-target overlap**************************** - - ********************************target-to-query overlap**************************** - case 3: u = 1, rev = 0 in the view of target: direction is 0 - query: AAATAATATCCCCCCGCG target: GGGCCGGCAAATAATAT (use ex-target as query) - ||||||||| <---> ||||||||| - target: GGGCCGGCAAATAATAT query: AAATAATATCCCCCCGCG (use ex-query as target) - - case 4: u = 1, rev = 1 in the view of target: direction is 1 - query: AAATAATATCCCCCCGCG target: ATATTATTTGCCGGCCC (use ex-target as query) - ||||||||| <---> ||||||||| - target: GGGCCGGCAAATAATAT query: CGCGGGGGATATTATTT (use ex-query as target) - ********************************target-to-query overlap**************************** - - p->ul: |____________31__________|__________1___________|______________32_____________| - qns direction of overlap length of this node (not overlap length) - (in the view of query) - p->v : |___________31___________|__________1___________| - tns reverse direction of overlap - (in the view of target) - p->ol: overlap length - **/ - - - //all outnode of v should be set to "not reduce" - for (i = 0; i < nv; ++i) mark[av[i].v] = 1; - - ///length of node (not overlap length) - ///av[nv-1] is longest out-dege - /** - * v--------------- - * w1--------------- - * w2-------------- - * w3-------------- - * w4-------------- - * w5------------- - * for v, the longest out-edge is v->w5 - **/ - L = asg_arc_len(av[nv-1]) + fuzz; - - - for (i = 0; i < nv; ++i) { - //w is an out-node of v - uint32_t w = av[i].v; - - uint32_t j, nw = asg_arc_n(g, w); - asg_arc_t *aw = asg_arc_a(g, w); - ///if w has already been reduced - if (mark[av[i].v] != 1) continue; - - for (j = 0; j < nw && asg_arc_len(aw[j]) + asg_arc_len(av[i]) <= L; ++j) - if (mark[aw[j].v]) mark[aw[j].v] = 2; - } - #if 0 - for (i = 0; i < nv; ++i) { - uint32_t w = av[i].v; - uint32_t j, nw = asg_arc_n(g, w); - asg_arc_t *aw = asg_arc_a(g, w); - for (j = 0; j < nw && (j == 0 || asg_arc_len(aw[j]) < fuzz); ++j) - if (mark[aw[j].v]) mark[aw[j].v] = 2; - } - #endif - //remove edges - for (i = 0; i < nv; ++i) { - if (mark[av[i].v] == 2) av[i].del = 1, ++n_reduced; - mark[av[i].v] = 0; - } - } - free(mark); - - if(VERBOSE >= 1) - { - fprintf(stderr, "[M::%s] transitively reduced %d arcs\n", __func__, n_reduced); - } - - if (n_reduced) { - asg_cleanup(g); - asg_symm(g); - } - - if(VERBOSE >= 1) - { - fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime); - } - return n_reduced; -} - - - - - ///max_ext is 4 int asg_cut_tip(asg_t *g, int max_ext) { @@ -26086,13 +25728,12 @@ int max_hang, int min_ovlp) -void clean_graph( -int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, -long long n_read, uint64_t* readLen, long long mini_overlap_length, -long long max_hang_length, long long clean_round, long long gap_fuzz, -float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, -long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr, -ma_sub_t **coverage_cut_ptr, int debug_g) +void clean_graph(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, + long long n_read, uint64_t* readLen, long long mini_overlap_length, + long long max_hang_length, long long clean_round, long long gap_fuzz, + float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, + long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr, + ma_sub_t **coverage_cut_ptr, int debug_g) { ma_sub_t *coverage_cut = *coverage_cut_ptr; asg_t *sg = *sg_ptr; @@ -26330,52 +25971,50 @@ ma_sub_t **coverage_cut_ptr, int debug_g) *sg_ptr = sg; } -void build_string_graph_without_clean( -int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, -long long n_read, uint64_t* readLen, long long mini_overlap_length, -long long max_hang_length, long long clean_round, long long gap_fuzz, -float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, -long long bubble_dist, int read_graph, int write) +void build_string_graph_without_clean(int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, + long long n_read, uint64_t* readLen, long long mini_overlap_length, + long long max_hang_length, long long clean_round, long long gap_fuzz, + float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, + long long bubble_dist, int read_graph, int write) { - R_to_U ruIndex; - init_R_to_U(&ruIndex, n_read); - asg_t *sg = NULL; - ma_sub_t* coverage_cut = NULL; + R_to_U ruIndex; + init_R_to_U(&ruIndex, n_read); + asg_t *sg = NULL; + ma_sub_t* coverage_cut = NULL; - // debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs", sources, reverse_sources, -1, "clean"); - - ///actually min_thres = asm_opt.max_short_tip + 1 there are asm_opt.max_short_tip reads - min_thres = asm_opt.max_short_tip + 1; + // debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs", sources, reverse_sources, -1, "clean"); - if (asm_opt.flag & HA_F_VERBOSE_GFA) - { - if(load_debug_graph(&sg, &sources, &coverage_cut, output_file_name, &reverse_sources, &ruIndex)) - { - fprintf(stderr, "debug gfa has been loaded\n"); - - clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length, - max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio, - output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 1); - asg_destroy(sg); - free(coverage_cut); - destory_R_to_U(&ruIndex); - return; - } - } - - if (asm_opt.write_index_to_disk && write) - { - write_all_data_to_disk(sources, reverse_sources, - &R_INF, output_file_name); - } + ///actually min_thres = asm_opt.max_short_tip + 1 there are asm_opt.max_short_tip reads + min_thres = asm_opt.max_short_tip + 1; - try_rescue_overlaps(sources, reverse_sources, n_read, 4); + if (asm_opt.flag & HA_F_VERBOSE_GFA) + { + if (load_debug_graph(&sg, &sources, &coverage_cut, output_file_name, &reverse_sources, &ruIndex)) + { + fprintf(stderr, "debug gfa has been loaded\n"); - clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length, - max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio, - output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 0); - - asg_destroy(sg); - free(coverage_cut); - destory_R_to_U(&ruIndex); + clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length, + max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio, + output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 1); + asg_destroy(sg); + free(coverage_cut); + destory_R_to_U(&ruIndex); + return; + } + } + + if (asm_opt.write_index_to_disk && write) + { + write_all_data_to_disk(sources, reverse_sources, &R_INF, output_file_name); + } + + try_rescue_overlaps(sources, reverse_sources, n_read, 4); + + clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length, + max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio, + output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 0); + + asg_destroy(sg); + free(coverage_cut); + destory_R_to_U(&ruIndex); } diff --git a/Overlaps.h b/Overlaps.h index bebb4be..05cd11b 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -83,6 +83,7 @@ char* output_file_name); typedef struct { uint32_t s:31, del:1, e; uint8_t c; + uint8_t weak_contain; } ma_sub_t; void ma_hit_sub(int min_dp, ma_hit_t_alloc* sources, long long n_read, uint64_t* readLen, @@ -103,7 +104,6 @@ typedef struct { uint8_t no_l_indel; } asg_arc_t; - typedef struct { uint32_t len:31, circ:1; // len: length of the unitig; circ: circular if non-zero uint32_t start, end; // start: starting vertex in the string graph; end: ending vertex @@ -112,11 +112,10 @@ typedef struct { char *s; // unitig sequence is not null } ma_utg_t; - - typedef struct { uint32_t len:31, del:1; uint8_t c; + uint8_t weak_contain; } asg_seq_t; typedef struct { @@ -128,6 +127,10 @@ typedef struct { asg_seq_t *seq; uint64_t *idx; + uint32_t m_con, n_con; + uint64_t *contain; + uint64_t *con_idx; + uint8_t* seq_vis; uint32_t n_F_seq; @@ -505,6 +508,21 @@ void init_Edge_iter(asg_t* g, uint32_t v, asg_arc_t* new_edges, uint32_t new_edg int get_arc_t(Edge_iter* x, asg_arc_t* get); int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_flag, uint32_t negative_flag); +/************************************* + * Routines modified for containment * + *************************************/ + +void delete_single_edge(ma_hit_t_alloc *sources, const ma_sub_t *coverage_cut, uint32_t qn, uint32_t tn); +void delete_all_edges(ma_hit_t_alloc *sources, ma_sub_t *coverage_cut, uint32_t qn); +void ma_hit_contained_advance(ma_hit_t_alloc *sources, long long n_read, ma_sub_t *coverage_cut, R_to_U *ruIndex, int max_hang, int min_ovlp); +void asg_con_sort(asg_t *g); +void asg_con_index(asg_t *g); +asg_t *ma_sg_gen(const ma_hit_t_alloc* sources, long long n_read, const ma_sub_t *coverage_cut, int max_hang, int min_ovlp); +int asg_arc_del_trans(asg_t *g, int fuzz); + +/******************* + * Other rountines * + *******************/ inline int get_real_length(asg_t *g, uint32_t v, uint32_t* v_s) { diff --git a/Trio.cpp b/Trio.cpp index b4934a9..7473ba0 100644 --- a/Trio.cpp +++ b/Trio.cpp @@ -9,6 +9,7 @@ #include "kseq.h" #include "Process_Read.h" #include "htab.h" +#include "utils.h" #include "CommandLines.h" #define YAK_MAX_KMER 31 diff --git a/anchor.cpp b/anchor.cpp index 4f2ce43..6e14049 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -1,6 +1,7 @@ #include #include "htab.h" #include "ksort.h" +#include "utils.h" #include "Hash_Table.h" #define HA_KMER_GOOD_RATIO 0.333 diff --git a/extract.cpp b/extract.cpp index b739a30..b561e9f 100644 --- a/extract.cpp +++ b/extract.cpp @@ -3,21 +3,19 @@ #include "Process_Read.h" #include "khashl.h" #include "kseq.h" +#include "utils.h" typedef const char *cstr_t; KHASHL_CSET_INIT(KH_LOCAL, strset_t, ss, cstr_t, kh_hash_str, kh_eq_str) KHASHL_MAP_INIT(KH_LOCAL, hm64_t, h64, uint64_t, int, kh_hash_uint64, kh_eq_generic) KSTREAM_INIT(gzFile, gzread, 65536) -#define GFA_MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr)))) -#define GFA_REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr)))) - char *gfa_strdup(const char *src) { int32_t len; char *dst; len = strlen(src); - GFA_MALLOC(dst, len + 1); + MALLOC(dst, len + 1); memcpy(dst, src, len + 1); return dst; } @@ -25,7 +23,7 @@ char *gfa_strdup(const char *src) char *gfa_strndup(const char *src, size_t n) { char *dst; - GFA_MALLOC(dst, n + 1); + MALLOC(dst, n + 1); strncpy(dst, src, n); dst[n] = 0; return dst; @@ -42,7 +40,7 @@ char **gv_read_list(const char *o, int *n_) if (*p == ',' || *p == 0) { if (n == m) { m = m? m<<1 : 16; - GFA_REALLOC(s, m); + REALLOC(s, m); } s[n++] = gfa_strndup(q, p - q); if (*p == 0) break; @@ -63,7 +61,7 @@ char **gv_read_list(const char *o, int *n_) for (p = str.s; *p && !isspace(*p); ++p); if (n == m) { m = m? m<<1 : 16; - GFA_REALLOC(s, m); + REALLOC(s, m); } s[n++] = gfa_strndup(str.s, p - str.s); } @@ -94,7 +92,7 @@ void ha_extract_print(const All_reads *rs, int n_rounds, int n, char **list) for (j = 0; j < rs->total_reads; ++j) if (max_len < (int)Get_NAME_LENGTH(*rs, j)) max_len = Get_NAME_LENGTH(*rs, j); - GFA_MALLOC(s, max_len + 1); + MALLOC(s, max_len + 1); h = h64_init(); for (j = 0; j < rs->total_reads; ++j) { strncpy(s, Get_NAME(*rs, j), Get_NAME_LENGTH(*rs, j)); diff --git a/htab.cpp b/htab.cpp index 0c44400..202d768 100644 --- a/htab.cpp +++ b/htab.cpp @@ -9,6 +9,7 @@ #include "kseq.h" #include "ksort.h" #include "htab.h" +#include "utils.h" #define YAK_COUNTER_BITS 12 #define YAK_N_COUNTS (1<>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) -#endif - -#ifndef kroundup64 -#define kroundup64(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, x|=(x)>>32, ++(x)) -#endif - #ifndef klib_unused #if (defined __clang__ && __clang_major__ >= 3) || (defined __GNUC__ && __GNUC__ >= 3) #define klib_unused __attribute__ ((__unused__)) diff --git a/main.cpp b/main.cpp index ccc4ed1..5f42313 100644 --- a/main.cpp +++ b/main.cpp @@ -4,7 +4,7 @@ #include "Process_Read.h" #include "Assembly.h" #include "Levenshtein_distance.h" -#include "htab.h" +#include "utils.h" int main(int argc, char *argv[]) { diff --git a/overlap2.cpp b/overlap2.cpp new file mode 100644 index 0000000..c2a2157 --- /dev/null +++ b/overlap2.cpp @@ -0,0 +1,319 @@ +#include +#include "utils.h" +#include "CommandLines.h" +#include "Overlaps.h" + +/******************************* + * Dropping strong containment * + *******************************/ + +static ma_hit_t *get_specific_overlap_with_del(ma_hit_t_alloc *sources, const ma_sub_t *coverage_cut, uint32_t qn, uint32_t tn) +{ + if (coverage_cut[qn].del || coverage_cut[tn].del) return NULL; + ma_hit_t_alloc *x = &sources[qn]; + uint32_t i; + for (i = 0; i < x->length; i++) { + if (x->buffer[i].del) continue; + if (coverage_cut[Get_qn(x->buffer[i])].del) continue; + if (coverage_cut[Get_tn(x->buffer[i])].del) continue; + if (Get_tn(x->buffer[i]) == tn && Get_qn(x->buffer[i]) == qn) + return &x->buffer[i]; + } + return NULL; +} + +void delete_single_edge(ma_hit_t_alloc *sources, const ma_sub_t *coverage_cut, uint32_t qn, uint32_t tn) +{ + ma_hit_t *tmp = get_specific_overlap_with_del(sources, coverage_cut, qn, tn); + if (tmp != NULL) tmp->del = 1; +} + +void delete_all_edges(ma_hit_t_alloc *sources, ma_sub_t *coverage_cut, uint32_t qn) +{ + ma_hit_t_alloc* x = &sources[qn]; + uint32_t i; + for (i = 0; i < x->length; i++) { + x->buffer[i].del = 1; + delete_single_edge(sources, coverage_cut, Get_tn(x->buffer[i]), Get_qn(x->buffer[i])); + } + coverage_cut[qn].del = 1; +} + +void ma_hit_contained_advance(ma_hit_t_alloc *sources, long long n_read, ma_sub_t *coverage_cut, R_to_U *ruIndex, int max_hang, int min_ovlp) +{ + int32_t r; + long long i, j, n_strong_contain = 0, n_weak_contain = 0; + asg_arc_t t; + ma_hit_t *h = NULL; + ma_sub_t *sq = NULL; + ma_sub_t *st = NULL; + + for (i = 0; i < n_read; ++i) { + if (coverage_cut[i].del) continue; + for (j = 0; j < (long long)sources[i].length; j++) { + h = &sources[i].buffer[j]; + //check the corresponding two reads + sq = &coverage_cut[Get_qn(*h)]; + st = &coverage_cut[Get_tn(*h)]; + /****************************may have trio bugs********************************/ + if (sq->del || st->del) continue; + if (h->del) continue; + /****************************may have trio bugs********************************/ + r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + //assert(r != MA_HT_INT && r != MA_HT_SHORT_OVLP); + if (r == MA_HT_QCONT) { + if (h->ml || 1) { + h->del = 1; + ++n_strong_contain; + delete_single_edge(sources, coverage_cut, Get_tn(*h), Get_qn(*h)); + delete_all_edges(sources, coverage_cut, Get_qn(*h)); + set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0); + + // if (delete_all_edges_carefully(sources, coverage_cut, max_hang, min_ovlp, Get_qn(*h)) == 0) + // set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0); + // sq->del = 1; + // set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0); + } else { + sq->weak_contain = 1; + ++n_weak_contain; + } + } else if (r == MA_HT_TCONT) { + if (h->ml || 1) { + h->del = 1; + ++n_strong_contain; + delete_single_edge(sources, coverage_cut, Get_tn(*h), Get_qn(*h)); + delete_all_edges(sources, coverage_cut, Get_tn(*h)); + set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0); + + // if (delete_all_edges_carefully(sources, coverage_cut, max_hang, min_ovlp, Get_tn(*h)) == 0) + // set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0); + // st->del = 1; + // set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0); + } else { + st->weak_contain = 1; + ++n_weak_contain; + } + } + } + } + + transfor_R_to_U(ruIndex); + + for (i = 0; i < n_read; ++i) { + int m = 0; + for (j = 0; j < (long long)sources[i].length; j++) { + ma_hit_t *h = &(sources[i].buffer[j]); + if (h->del) continue; + /// both the qn and tn have not been deleted + if (coverage_cut[Get_qn(*h)].del != 1 && coverage_cut[Get_tn(*h)].del != 1) + h->del = 0, ++m; + else h->del = 1; + } + /// sources[i].length == 0 means all overlapped reads with read i are the contained reads + if (m == 0) coverage_cut[i].del = 1; + } + + fprintf(stderr, "[M::%s] %lld strong containments; %lld weak containments\n", __func__, + n_strong_contain, n_weak_contain); +} + +/************************************ + * Graph construction and reduction * + ************************************/ + +static inline void asg_con_push(asg_t *g, uint32_t lower, uint32_t upper, int rev) +{ + if (g->n_con == g->m_con) { + g->m_con = g->m_con? g->m_con<<1 : 16; + REALLOC(g->contain, g->m_con); + } + g->contain[g->n_con++] = (uint64_t)lower << 32 | upper << 1 | (!!rev); +} + +void asg_con_sort(asg_t *g) +{ + if (g->n_con > 1) radix_sort_ha64(g->contain, g->contain + g->n_con); +} + +void asg_con_index(asg_t *g) +{ + uint32_t i, k; + if (g->n_con == 0 || g->contain == 0) return; + if (g->con_idx) free(g->con_idx); + CALLOC(g->con_idx, g->n_seq); + for (k = 0, i = 1; i < g->n_con; ++i) + if (g->contain[k] != g->contain[i]) + g->contain[k++] = g->contain[i]; + g->n_con = k; + for (i = 1, k = 0; i <= g->n_con; ++i) + if (i == g->n_con || g->contain[i-1]>>32 != g->contain[i]>>32) + g->con_idx[g->contain[i-1]>>32] = (uint64_t)k << 32 | (i - k), k = i; +} + +asg_t *ma_sg_gen(const ma_hit_t_alloc* sources, long long n_read, const ma_sub_t *coverage_cut, int max_hang, int min_ovlp) +{ + size_t i, j; + asg_t *g; + g = asg_init(); + + // add seq to graph, seq just save the length of each read + for (i = 0; i < (uint64_t)n_read; ++i) { + ///if a read has been deleted, should we still add them? + asg_seq_set(g, i, coverage_cut[i].e - coverage_cut[i].s, coverage_cut[i].del); + g->seq[i].c = coverage_cut[i].c; + } + + g->seq_vis = (uint8_t*)calloc(g->n_seq*2, sizeof(uint8_t)); + + for (i = 0; i < (uint64_t)n_read; ++i) { + for (j = 0; j < sources[i].length; ++j) { + int r, ql, tl; + asg_arc_t t, *p; + const ma_hit_t *h = &sources[i].buffer[j]; + uint32_t qn, tn; + if (h->del) continue; + qn = Get_qn(*h); + tn = Get_tn(*h); + ql = coverage_cut[qn].e - coverage_cut[qn].s; + tl = coverage_cut[tn].e - coverage_cut[tn].s; + r = ma_hit2arc(h, ql, tl, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + assert(r >= 0); + if (r >= 0) { + p = asg_arc_pushp(g); + *p = t; + } else if (r == MA_HT_QCONT) { + assert(h->ml == 0); + asg_con_push(g, h->qns>>32, h->tn, h->rev); + } else if (r == MA_HT_TCONT) { + assert(h->ml == 0); + asg_con_push(g, h->tn, h->qns>>32, h->rev); + } + } + } + asg_cleanup(g); + g->r_seq = g->n_seq; + return g; +} + +typedef struct { + uint32_t len; + uint8_t mark; // can only be 0, 1 or 2 +} trinfo_t; + +// transitive reduction; see Myers, 2005 +int asg_arc_del_trans(asg_t *g, int fuzz) +{ + trinfo_t *info; + ///n_vtx = number of seq * 2; 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; + ///at first, all nodes should be set to vacant + CALLOC(info, n_vtx); + + /**v is the id+direction of a node, + * the high 31-bit is the id, + * and the lowest 1-bit is the direction + * (0 means query-to-target, 1 means target-to-query)**/ + for (v = 0; v < n_vtx; ++v) { + ///nv is the number of overlaps with v(qn+direction) + uint32_t L, i, nv = asg_arc_n(g, v); + ///av is the array of v + asg_arc_t *av = asg_arc_a(g, v); + ///that means in this direction, read v is not overlapped with any other reads + if (nv == 0) continue; // no hits + + // if the read itself has been removed + if (g->seq[v>>1].del) { + for (i = 0; i < nv; ++i) av[i].del = 1, ++n_reduced; + continue; + } + + /** + ********************************query-to-target overlap**************************** + case 1: u = 0, rev = 0 in the view of target: direction is 1 + query: CCCCCCCCTAATTAAAAT target: TAATTAAAATGGGGGG (use ex-target as query) + |||||||||| <---> |||||||||| + target: TAATTAAAATGGGGGG query: CCCCCCCCTAATTAAAAT (use ex-query as target) + + case 2: u = 0, rev = 1 in the view of target: direction is 0 + query: CCCCCCCCTAATTAAAAT target: CCCCCCATTTTAATTA (use ex-target as query) + |||||||||| <---> |||||||||| + target: TAATTAAAATGGGGGG query: ATTTTAATTAGGGGGGGG (use ex-query as target) + ********************************query-to-target overlap**************************** + + ********************************target-to-query overlap**************************** + case 3: u = 1, rev = 0 in the view of target: direction is 0 + query: AAATAATATCCCCCCGCG target: GGGCCGGCAAATAATAT (use ex-target as query) + ||||||||| <---> ||||||||| + target: GGGCCGGCAAATAATAT query: AAATAATATCCCCCCGCG (use ex-query as target) + + case 4: u = 1, rev = 1 in the view of target: direction is 1 + query: AAATAATATCCCCCCGCG target: ATATTATTTGCCGGCCC (use ex-target as query) + ||||||||| <---> ||||||||| + target: GGGCCGGCAAATAATAT query: CGCGGGGGATATTATTT (use ex-query as target) + ********************************target-to-query overlap**************************** + + p->ul: |____________31__________|__________1___________|______________32_____________| + qns direction of overlap length of this node (not overlap length) + (in the view of query) + p->v : |___________31___________|__________1___________| + tns reverse direction of overlap + (in the view of target) + p->ol: overlap length + **/ + + // all outnode of v should be set to "not reduce" + for (i = 0; i < nv; ++i) { + uint32_t w = av[i].v; + info[w].mark = g->seq[w>>1].del? 2 : 1; + info[w].len = asg_arc_len(av[i]); + } + + // length of node (not overlap length) + // av[nv-1] is longest out-dege + /** + * v--------------- + * w1--------------- + * w2-------------- + * w3-------------- + * w4-------------- + * w5------------- + * for v, the longest out-edge is v->w5 + **/ + L = asg_arc_len(av[nv-1]) + fuzz; + + for (i = 0; i < nv; ++i) { + uint32_t w = av[i].v; + uint32_t j, nw = asg_arc_n(g, w); + asg_arc_t *aw = asg_arc_a(g, w); + if (info[w].mark != 1) continue; + for (j = 0; j < nw; ++j) { + uint32_t x, sum = asg_arc_len(aw[j]) + asg_arc_len(av[i]); + if (sum > L) break; + x = aw[j].v; + if (info[x].mark == 1 && sum < info[x].len + fuzz && sum + fuzz > info[x].len) + info[x].mark = 2; + } + } + #if 0 + for (i = 0; i < nv; ++i) { + uint32_t w = av[i].v; + uint32_t j, nw = asg_arc_n(g, w); + asg_arc_t *aw = asg_arc_a(g, w); + for (j = 0; j < nw && (j == 0 || asg_arc_len(aw[j]) < fuzz); ++j) + if (info[aw[j].v].mark) info[aw[j].v].mark = 2; + } + #endif + // remove edges + for (i = 0; i < nv; ++i) { + if (info[av[i].v].mark == 2) av[i].del = 1, ++n_reduced; + info[av[i].v].mark = 0; + } + } + free(info); + + if (n_reduced) { + asg_cleanup(g); + asg_symm(g); + } + return n_reduced; +} diff --git a/sys.cpp b/sys.cpp index 87b5ac2..5a9ebeb 100644 --- a/sys.cpp +++ b/sys.cpp @@ -1,6 +1,6 @@ #include #include -#include "htab.h" +#include "utils.h" int yak_verbose = 3; diff --git a/utils.h b/utils.h new file mode 100644 index 0000000..e151257 --- /dev/null +++ b/utils.h @@ -0,0 +1,33 @@ +#ifndef HA_UTILS_H +#define HA_UTILS_H + +#include + +#ifndef MALLOC +#define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr)))) +#endif +#ifndef CALLOC +#define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr)))) +#endif +#ifndef REALLOC +#define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr)))) +#endif + +#ifndef kroundup32 +#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) +#endif + +#ifndef kroundup64 +#define kroundup64(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, x|=(x)>>32, ++(x)) +#endif + +void radix_sort_ha64(uint64_t *st, uint64_t *en); + +double yak_cputime(void); +void yak_reset_realtime(void); +double yak_realtime(void); +long yak_peakrss(void); +double yak_peakrss_in_gb(void); +double yak_cpu_usage(void); + +#endif