From 809ca18dda31b3789ed76d68a98b5c6308272e40 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Fri, 10 Sep 2021 17:17:44 -0400 Subject: [PATCH] r373 --- Assembly.cpp | 4 +- CommandLines.cpp | 13 + CommandLines.h | 4 +- Overlaps.cpp | 100 ++--- Overlaps.h | 8 +- Purge_Dups.cpp | 326 +++++++++++++++-- Purge_Dups.h | 1 + anchor.cpp | 2 +- hic.cpp | 235 ++++++++++-- hic.h | 3 +- horder.cpp | 248 +++++++++++-- horder.h | 2 + htab.cpp | 630 +++++++++++++++----------------- htab.h | 13 +- inter.cpp | 25 +- rcut.cpp | 50 ++- rcut.h | 2 +- sketch.cpp | 934 ++++++++++++++++++++--------------------------- 18 files changed, 1569 insertions(+), 1031 deletions(-) diff --git a/Assembly.cpp b/Assembly.cpp index ba5aa97..588dcd6 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1636,8 +1636,8 @@ void ug_idx_build(ma_ug_t *ug, int hap_n) { int flag = asm_opt.flag&HA_F_NO_HPC, i; asm_opt.flag -= flag; - ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), hap_n); - ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), hap_n); + // ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), hap_n, hap_n); + // ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), hap_n); ha_ovec_buf_t **b = NULL; // overlap and correct reads diff --git a/CommandLines.cpp b/CommandLines.cpp index 81d7f4a..064e2b7 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -44,6 +44,8 @@ static ko_longopt_t long_options[] = { { "dp-er", ko_required_argument, 330}, { "max-kocc", ko_required_argument, 331}, { "hg-size", ko_required_argument, 332}, + { "ul", ko_required_argument, 333}, + { "unskew", ko_no_argument, 334}, { 0, 0, 0 } }; @@ -155,6 +157,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->hic_enzymes = NULL; asm_opt->hic_reads[0] = NULL; asm_opt->hic_reads[1] = NULL; + asm_opt->ar = NULL; asm_opt->thread_num = 1; asm_opt->k_mer_length = 51; asm_opt->hic_mer_length = 31; @@ -245,6 +248,7 @@ void destory_opt(hifiasm_opt_t* asm_opt) if(asm_opt->hic_enzymes != NULL) destory_enzyme(asm_opt->hic_enzymes); if(asm_opt->hic_reads[0] != NULL) destory_enzyme(asm_opt->hic_reads[0]); if(asm_opt->hic_reads[1] != NULL) destory_enzyme(asm_opt->hic_reads[1]); + if(asm_opt->ar != NULL) destory_enzyme(asm_opt->ar); } void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round) @@ -493,6 +497,13 @@ int check_option(hifiasm_opt_t* asm_opt) return 0; } + if(asm_opt->ar != NULL && check_hic_reads(asm_opt->ar, "UL") == 0) return 0; + if(asm_opt->ar != NULL && asm_opt->ar->n == 0) + { + fprintf(stderr, "[ERROR] wrong UL reads (--ul)\n"); + return 0; + } + if(asm_opt->b_low_cov < 0) { fprintf(stderr, "[ERROR] must >= 0 (--b-cov)\n"); @@ -725,6 +736,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 330) asm_opt->dp_e = atof(opt.arg); else if (c == 331) asm_opt->max_kmer_cnt = atol(opt.arg); else if (c == 332) asm_opt->hg_size = inter_gsize(opt.arg); + else if (c == 333) get_hic_enzymes(opt.arg, &(asm_opt->ar), 0); + else if (c == 334) asm_opt->flag |= HA_F_USKEW; else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg); diff --git a/CommandLines.h b/CommandLines.h index 8a15dc9..911c498 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.16.0-r369" +#define HA_VERSION "0.16.1-r373" #define VERBOSE 0 @@ -21,6 +21,7 @@ #define HA_F_HIGH_HET 0x400 #define HA_F_PARTITION 0x800 #define HA_F_FAST 0x1000 +#define HA_F_USKEW 0x2000 #define HA_MIN_OV_DIFF 0.02 // min sequence divergence in an overlap @@ -40,6 +41,7 @@ typedef struct { char *extract_list; enzyme *hic_reads[2]; enzyme *hic_enzymes; + enzyme *ar; int extract_iter; int thread_num; int k_mer_length; diff --git a/Overlaps.cpp b/Overlaps.cpp index a0dbade..a1b794b 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -1475,14 +1475,6 @@ R_to_U* ruIndex, int max_hang, int min_ovlp) delete_all_edges(sources, coverage_cut, Get_qn(*h)); set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0, NULL); - - // 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) { @@ -1491,15 +1483,6 @@ R_to_U* ruIndex, int max_hang, int min_ovlp) delete_all_edges(sources, coverage_cut, Get_tn(*h)); set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0, NULL); - - // 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); } } } @@ -9673,7 +9656,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag) r_flag[rId] = 0; } - return C_bases/R_bases; + return R_bases == 0? 0:C_bases/R_bases; } uint32_t get_ug_coverage_aggressive(ma_ug_t *ug, uint32_t uID, asg_t* read_g, @@ -13158,7 +13141,7 @@ long long gap_fuzz, bub_label_t* b_mask_t) ///asm_opt.purge_simi_thres = asm_opt.purge_simi_rate_hic; adjust_utg_by_primary(©_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, - max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 1/**0**/); + max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 0); print_utg(copy_ug, copy_sg, coverage_cut, output_file_name, sources, ruIndex, max_hang, min_ovlp, &new_rtg_edges); @@ -13182,7 +13165,6 @@ long long gap_fuzz, bub_label_t* b_mask_t) } // char* gfa_name = (char*)malloc(strlen(output_file_name)+50); - // FILE* output_file = NULL; // sprintf(gfa_name, "%s.pre.clean_d_utg.noseq.gfa", output_file_name); // FILE* output_file = fopen(gfa_name, "w"); // ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file); @@ -13744,7 +13726,7 @@ void kt_u_trans_t_symm(kv_u_trans_t *ta, ma_ug_t *ug) if (i == n || a[i].tn != a[st].tn) { get_u_trans_spec(ta, a[st].tn, a[st].qn, &r_a, &r_n); - if(i - st == 1 && r_n == 0) + if(i - st == 1 && r_n == 0)///should be always here { st = i; continue; @@ -14811,9 +14793,37 @@ void label_r_set(buf_t* b, R_to_U* ruIndex, ma_ug_t *ug, uint32_t flag) } } +inline int trio_check(ma_ug_t *ug, uint32_t *a, uint32_t a_n, uint32_t flag) +{ + if(flag != FATHER && flag != MOTHER) return 0; + uint32_t flag_occ = 0, non_flag_occ = 0, ambigious = 0, u_n = 0, f, nf, ab, k; + for (k = 0; k < a_n; k++) { + get_unitig_trio_flag(&(ug->u.a[a[k]>>1]), flag, &f, &nf, &ab); + flag_occ += f; + non_flag_occ += nf; + ambigious += ab; + u_n += ug->u.a[a[k]>>1].n; + } + if((flag_occ+non_flag_occ) == 0) return 0; + if(flag_occ <= ((non_flag_occ+flag_occ)*0.75)) return 0; + if(non_flag_occ == 0 && flag_occ >= 20) return 1; + if(u_n >= 100) + { + if(flag_occ < u_n*DOUBLE_CHECK_THRES) return 0; + } + else if(u_n >= 50) + { + if(flag_occ < u_n*DOUBLE_CHECK_THRES*0.5) return 0; + } + else + { + if(flag_occ < u_n*DOUBLE_CHECK_THRES*0.25) return 0; + } + return 1; +} int asg_arc_cut_trio_long_tip_primary(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, -R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov, utg_trans_t *o) +R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t trio_flag, hap_cov_t *cov, utg_trans_t *o) { double startTime = Get_T(); ///the reason is that each read has two direction (query->target, target->query) @@ -14874,6 +14884,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov, utg } if(ll >= (v_maxLen*drop_ratio)) continue; + if(trio_check(ug, b.b.a, b.b.n, trio_flag)) continue; n_reduced++; operation = TRIM; @@ -14934,7 +14945,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hap_cov_t *cov, utg } int asg_arc_cut_trio_long_tip_primary_complex(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, -R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_threshold, hap_cov_t *cov, utg_trans_t *o) +R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_threshold, hap_cov_t *cov, utg_trans_t *o, uint32_t trio_flag) { double startTime = Get_T(); uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, in, flag, operation; @@ -14958,6 +14969,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_thre &max_stop_baseLen, 1, &b); if(return_flag != MUL_INPUT) continue; + if(trio_check(ug, b.b.a, b.b.n, trio_flag)) continue; in = convex^1; get_real_length(g, convex, &convex); convex = convex^1; @@ -15177,17 +15189,11 @@ hap_cov_t *cov, utg_trans_t *o) &max_stop_baseLen, 1, &b); if(return_flag != END_TIPS) continue; + if(trio_check(ug, b.b.a, b.b.n, trio_flag)) continue; flag = check_different_haps(g, ug, read_sg, av[base_maxLen_i].v, av[i].v, reverse_sources, &b_0, &b_1, ruIndex, cov->is_r_het, miniedgeLen, 1); - // if((av[i].v>>1) == 255 && (av[base_maxLen_i].v>>1) == 33) - // if((av[i].v>>1) == 1852 && (av[base_maxLen_i].v>>1) == 2441) - // { - // fprintf(stderr, "max-utg%.6ul (%u), p-utg%.6ul (%u)\n", - // (av[base_maxLen_i].v>>1)+1, av[base_maxLen_i].v, (av[i].v>>1)+1, av[i].v); - // } - // #define UNAVAILABLE (uint32_t)-1 // #define PLOID 0 // #define NON_PLOID 1 @@ -15445,7 +15451,7 @@ R_to_U* ruIndex, uint32_t positive_flag, float drop_rate) int asg_arc_cut_trio_long_equal_tips_assembly_complex(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_t stops_threshold, -hap_cov_t *cov, utg_trans_t *o) +hap_cov_t *cov, utg_trans_t *o, uint32_t trio_flag) { double startTime = Get_T(); uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, in, flag; @@ -15471,6 +15477,7 @@ hap_cov_t *cov, utg_trans_t *o) &max_stop_baseLen, 1, &b); if(return_flag != MUL_INPUT) continue; + if(trio_check(ug, b.b.a, b.b.n, trio_flag)) continue; in = convex^1; get_real_length(g, convex, &convex); convex = convex^1; @@ -15955,10 +15962,11 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov) if(just_bubble_pop == 0) { ///need consider tangles - asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, NULL); + asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, trio_flag, cov, NULL); + // if(trio_flag == MOTHER) print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, cov, NULL); - asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL); - asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL); + asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL, trio_flag); + asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL, trio_flag); detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL, cov->is_r_het); ///need consider tangles ///note we need both the read graph and the untig graph @@ -15972,10 +15980,11 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov) } ///print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); - + magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate); resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, cov->is_r_het, trio_flag, drop_ratio); drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het); all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, cov->is_r_het, DOUBLE_CHECK_THRES, asm_opt.trio_flag_occ_thres); + // if(trio_flag == MOTHER) print_untig_by_read(ug, "m54329U_190827_173812/30214441/ccs", (uint32_t)-1, NULL, NULL, "bf-16"); if(is_first) { is_first = 0; @@ -16031,10 +16040,10 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov) if(just_bubble_pop == 0) { ///need consider tangles - asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, NULL); + asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, (uint32_t)-1, cov, NULL); asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, cov, NULL); - asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL); - asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL); + asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL, (uint32_t)-1); + asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL, (uint32_t)-1); detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL, cov->is_r_het); if(round != T_ROUND) { @@ -16094,10 +16103,10 @@ int min_ovlp, hap_cov_t *cov) asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, o, 1); ///need consider tangles - asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, cov, o); + asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, (uint32_t)-1, cov, o); asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, cov, o); - asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, o); - asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, o); + asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, o, (uint32_t)-1); + asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, o, (uint32_t)-1); detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, o, cov->is_r_het); cur_cons = get_graph_statistic(g); @@ -23920,6 +23929,8 @@ int write_ruIndex(R_to_U* ruIndex, char* read_file_name) fwrite(ruIndex->index, sizeof(ruIndex->index[0]), ruIndex->len, fp); fwrite(R_INF.trio_flag, sizeof(R_INF.trio_flag[0]), ruIndex->len, fp); // fwrite(ruIndex->is_het, 1, ruIndex->len, fp); + fwrite(&(asm_opt.hom_global_coverage_set), sizeof(asm_opt.hom_global_coverage_set), 1, fp); + fwrite(&(asm_opt.hom_global_coverage), sizeof(asm_opt.hom_global_coverage), 1, fp); free(index_name); fflush(fp); fclose(fp); @@ -23947,6 +23958,9 @@ int load_ruIndex(R_to_U* ruIndex, char* read_file_name) // CALLOC(ruIndex->is_het, ruIndex->len); // f_flag += fread(ruIndex->is_het, 1, ruIndex->len, fp); + f_flag += fread(&(asm_opt.hom_global_coverage_set), sizeof(asm_opt.hom_global_coverage_set), 1, fp); + f_flag += fread(&(asm_opt.hom_global_coverage), sizeof(asm_opt.hom_global_coverage), 1, fp); + free(index_name); fflush(fp); fclose(fp); @@ -24707,7 +24721,7 @@ uint32_t collect_p_trans, uint32_t collect_p_trans_f) fprintf(stderr, "[M::%s] primary contig coverage range: [%d, infinity]\n", __func__, asm_opt.recover_atg_cov_min); } - + skip_purge: recover_utg_by_coverage(ug, read_g, coverage_cut, sources, ruIndex, cov->t_ch); if(i_cov) @@ -30869,12 +30883,10 @@ ma_sub_t **coverage_cut_ptr, int debug_g) output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources, asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length); - if (asm_opt.flag & HA_F_VERBOSE_GFA) { write_debug_graph(sg, sources, coverage_cut, output_file_name, n_read, reverse_sources, ruIndex); debug_gfa:; - set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length); } if (ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt)) diff --git a/Overlaps.h b/Overlaps.h index b4eec32..8f215c5 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -190,7 +190,7 @@ void print_gfa(asg_t *g); typedef struct { size_t n, m; uint64_t *a; } asg64_v; -typedef struct { size_t n, m; ma_utg_t *a; } ma_utg_v; +typedef struct { size_t n, m; ma_utg_t *a; int h;} ma_utg_v; typedef struct { ma_utg_v u; @@ -891,6 +891,12 @@ typedef struct { buf_t b0, b1; } utg_trans_t; +typedef struct { + ma_ug_t *ug; + kvec_t(uint64_t) idx; + kvec_t(uint32_t) dst; +} spg_t; + void init_hc_links(hc_links* link, uint64_t ug_num, trans_chain* t_ch); void destory_hc_links(hc_links* link); uint64_t get_bub_pop_max_dist(asg_t *g, buf_t *b); diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index bc15e08..0ef3be0 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -51,7 +51,7 @@ typedef struct { typedef struct { hap_alignment_struct* buf; uint32_t num_threads; - + uint8_t *hh; ma_ug_t *ug; asg_t *read_g; ma_hit_t_alloc* sources; @@ -471,6 +471,30 @@ void destory_hap_alignment_struct(hap_alignment_struct* x) kv_destroy(x->u_can.a); } +uint8_t *init_pip_hh(asg_t *rg, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, long long sc) +{ + uint8_t *c = NULL; CALLOC(c, rg->n_seq); + ma_hit_t *h = NULL; + uint32_t i, k; + long long R_Base, C_Base; + for (i = 0; i < rg->n_seq; i++) { + if(sc <= 0){ + c[i] = 1; + continue; + } + R_Base = (coverage_cut[i].e - coverage_cut[i].s); + for (k = 0, C_Base = 0; k < reverse_sources[i].length; k++){ + h = &(reverse_sources[i].buffer[k]); + C_Base += (Get_qe((*h)) - Get_qs((*h))); + } + C_Base = (R_Base!=0?(C_Base/R_Base):0); + C_Base /= sc; + c[i] = 1; + if(C_Base < REV_W) c[i] = REV_W - C_Base; + } + return c; +} + void init_hap_alignment_struct_pip(hap_alignment_struct_pip* x, uint32_t num_threads, uint32_t n_seq, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, ma_sub_t *coverage_cut, uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chain_rate, hap_overlaps_list* all_ovlp, hap_cov_t *cov) @@ -496,6 +520,17 @@ uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chai x->chain_rate = chain_rate; x->all_ovlp = all_ovlp; x->cov = cov; + long long sc = -1; + if(asm_opt.hom_global_coverage_set) { + sc = asm_opt.hom_global_coverage*1.75; + } + else { + if(asm_opt.hom_global_coverage > 0){ + sc = ((int)(((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE)))*1.75; + } + } + if(sc <= 0) sc = -1; + x->hh = init_pip_hh(read_g, reverse_sources, coverage_cut, sc); } @@ -508,6 +543,7 @@ void destory_hap_alignment_struct_pip(hap_alignment_struct_pip* x) } free(x->buf); + free(x->hh); } void init_hap_overlaps_list(hap_overlaps_list* x, uint32_t num) @@ -677,16 +713,10 @@ void deduplicate_edge(kvec_asg_arc_t_offset* u_buffer) long long i = u_buffer->a.n - 1, k, i_off; uint32_t v = u_buffer->a.a[i].x.ul>>33, m; - for (; i >= 0; i--) - { - if((u_buffer->a.a[i].x.ul>>33) != v) - { - break; - } + for (; i >= 0; i--) { + if((u_buffer->a.a[i].x.ul>>33) != v) break; } - ///fprintf(stderr, "u_buffer->a.n: %u, i: %lld\n", u_buffer->a.n, i); - i = i + 1; for (m = i; i < (long long)u_buffer->a.n; i++) { @@ -960,7 +990,7 @@ ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match, uint32_t i, j, qn, tn, is_Unitig, uId, min_count = 0, max_count = 0, cutoff = 0;; for (i = 0; i < Len; i++) { - if(cutoff > CUTOFF_THRES) + if(cutoff > CUTOFF_THRES && cutoff > (Len>>1)) { max_count = 0; min_count = Len; @@ -2254,6 +2284,7 @@ uint32_t* xBeg, uint32_t* xEnd, uint32_t* yBeg, uint32_t* yEnd) #define generic_key(x) (x) KRADIX_SORT_INIT(i32, int32_t, generic_key, sizeof(int32_t)) +KRADIX_SORT_INIT(ru32, uint32_t, generic_key, sizeof(uint32_t)) long long get_chain_score(ma_utg_t *xReads, asg_t *read_g, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* idx, ma_hit_t_alloc* reverse_sources, long long xBegPos, long long xEndPos) @@ -2820,7 +2851,7 @@ long long y_readLen) { u_can->a.n = 0; if(u_buffer->a.n == 0) return; - uint32_t i = 0, anchor_i = 0, m = 1, break_point = (uint32_t)-1, is_merge; + uint32_t i = 0, /**anchor_i = 0, **/m = 1, break_point = (uint32_t)-1, is_merge; qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment_chaining); @@ -2832,14 +2863,14 @@ long long y_readLen) if(u_buffer->a.a[m-1].x.el == u_buffer->a.a[i].x.el) { if(u_buffer->a.a[m-1].Off == u_buffer->a.a[i].Off) is_merge = 1; - if(is_merge == 0 && (Get_xOff(u_buffer->a.a[m-1].Off)==Get_xOff(u_buffer->a.a[i].Off))) - { - if((Get_yOff(u_buffer->a.a[i].Off)-(Get_yOff(u_buffer->a.a[m-1].Off))) == - (i-anchor_i))///not sure why, does it use for tolerate indels in overlaps? - { - is_merge = 1; - } - } + ///I think we don't need the following merging + // if(is_merge == 0 && (Get_xOff(u_buffer->a.a[m-1].Off)==Get_xOff(u_buffer->a.a[i].Off))) + // { + // if((Get_yOff(u_buffer->a.a[i].Off)-(Get_yOff(u_buffer->a.a[m-1].Off))) == (i-anchor_i))///not sure why, does it use for tolerate indels in overlaps? + // { + // is_merge = 1; + // } + // } if(is_merge) { @@ -2848,7 +2879,7 @@ long long y_readLen) } } u_buffer->a.a[m] = u_buffer->a.a[i]; - anchor_i = i; + // anchor_i = i; if(u_buffer->a.a[m].x.el != u_buffer->a.a[m-1].x.el) break_point = m; m++; } @@ -2957,6 +2988,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) kvec_t_u8_warp* flag_vec = &(hap_buf->buf[tid].u_buffer_flag); uint64_t cov_threshold = hap_buf->cov_threshold; hap_cov_t *cov = hap_buf->cov; + uint8_t *hh = hap_buf->hh, hhc; if(hap_buf->cov_threshold < 0) cov_threshold = (uint64_t)-1; ma_utg_t *xReads = NULL, *yReads = NULL; ma_hit_t_alloc *xR = NULL; @@ -3036,7 +3068,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) for (k = 0; k < xReads->n; k++) { xR = &(reverse_sources[xReads->a[k]>>33]); - + hhc = hh[xReads->a[k]>>33]; for (j = 0; j < xR->length; j++) { h = &(xR->buffer[j]); @@ -3078,7 +3110,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) if(((t_offset.Off>>32) == (uint32_t)-1) || (((uint32_t)t_offset.Off) == (uint32_t)-1)) continue; t_offset.x = t; - t_offset.weight = 1; + t_offset.weight = hhc; kv_push(asg_arc_t_offset, u_buffer->a, t_offset); } @@ -3155,7 +3187,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) } all_ovlp->x[xUid].a.n = m; **/ - + hap_align_x = NULL; for (k = m; k < all_ovlp->x[xUid].a.n; k++) { @@ -5255,6 +5287,243 @@ void collect_purge_trans_cov(ma_ug_t *ug, hap_overlaps_list* ha, hap_cov_t *cov, } } +/** +typedef struct { + uint32_t qn, qs, qe; + uint32_t tn, ts, te; + uint32_t oid; + uint8_t rev; +}scg_hits; + +typedef struct { + scg_hits *a; + size_t n,m; + kvec_t(uint64_t) idx; +}scg_hits_v; + +#define scg_key_qtn(a) ((((uint64_t)(a).qn)<<32)|((uint64_t)(a).tn)) +KRADIX_SORT_INIT(scg_qtn, scg_hits, scg_key_qtn, 8) +#define scg_key_qts(a) ((((uint64_t)(a).qs)<<32)|((uint64_t)(a).ts)) +KRADIX_SORT_INIT(scg_qts, scg_hits, scg_key_qts, 8) +#define scg_key_qte(a) ((((uint64_t)(a).qe)<<32)|((uint64_t)(a).te)) +KRADIX_SORT_INIT(scg_qte, scg_hits, scg_key_qte, 8) +#define scg_key_rev(a) ((a).rev) +KRADIX_SORT_INIT(scg_rev, scg_hits, scg_key_rev, member_size(scg_hits, rev)) + +inline void rev_scg_hits(scg_hits *p, spg_t *scg) +{ + if(p->rev){ + uint32_t t; + p->ts = scg->ug->u.a[p->tn].len - p->ts - 1; + p->te = scg->ug->u.a[p->tn].len - (p->te - 1) - 1; + t = p->ts; p->ts = p->te; p->te = t; p->te++; + } +} + +#define arc_first(g, v) ((g)->arc[(g)->idx[(v)]>>32]) +scg_hits_v *get_scg_hits_v(scg_hits_v *vp, spg_t *scg) +{ + scg_hits_v *hh = NULL; CALLOC(hh, 1); + ma_utg_v *u = &(scg->ug->u); + uint32_t i, mn, *ma = NULL; + uint64_t offset, *idx = NULL; + + for (i = 0; i < scg->idx.n; i++) { + mn = (uint32_t)scg->idx.a[i]; + ma = scg->dst.a + (scg->idx.a[i]>>32); + } + + + + + return hh; +} + +void refine_scg(spg_t *scg, ma_ug_t *lug, hap_overlaps_list *ha, hap_cov_t *cov, uint64_t* position_index) +{ + scg_hits_v vp; kv_init(vp); + uint32_t v, i, k, st, c[2]; + hap_overlaps *x = NULL; + u_trans_t *z = NULL; + scg_hits *p = NULL; + + for (v = 0; v < cov->t_ch->k_trans.n; v++){ + z = &(cov->t_ch->k_trans.a[v]); + if(z->del) continue; + kv_pushp(scg_hits, vp, &p); + p->rev = z->rev; p->oid = (uint32_t)-1; + p->qn = z->qn; p->qs = z->qs; p->qe = z->qe; + p->tn = z->tn; p->ts = z->ts; p->te = z->te; + // rev_scg_hits(p, scg); + kv_pushp(scg_hits, vp, &p); + p->rev = z->rev; p->oid = (uint32_t)-1; + p->qn = z->tn; p->qs = z->ts; p->qe = z->te; + p->tn = z->qn; p->ts = z->qs; p->te = z->qe; + // rev_scg_hits(p, scg); + } + + for (v = 0; v < ha->num; v++){ + for (i = 0; i < ha->x[v].a.n; i++){ + x = &(ha->x[v].a.a[i]); + st = cov->t_ch->k_trans.n; + chain_origin_trans_uid_by_purge(x, lug, cov, position_index); + for (k = st; k < cov->t_ch->k_trans.n; k++){ + z = &(cov->t_ch->k_trans.a[k]); + if(z->del) continue; + kv_pushp(scg_hits, vp, &p); + p->rev = z->rev; p->oid = v; + p->qn = z->qn; p->qs = z->qs; p->qe = z->qe; + p->tn = z->tn; p->ts = z->ts; p->te = z->te; + // rev_scg_hits(p, scg); + kv_pushp(scg_hits, vp, &p); + p->rev = z->rev; p->oid = v; + p->qn = z->tn; p->qs = z->ts; p->qe = z->te; + p->tn = z->qn; p->ts = z->qs; p->te = z->qe; + // rev_scg_hits(p, scg); + } + cov->t_ch->k_trans.n = st; + } + } + + ///two scg_hits might be totally equal; must remove first + radix_sort_scg_qtn(vp.a, vp.a + vp.n); + for (st = 0, i = 1; i <= vp.n; ++i){ + if (i == vp.n || vp.a[i].qn != vp.a[st].qn || vp.a[i].tn != vp.a[st].tn){ + if(i - st > 1) radix_sort_scg_rev(vp.a+st, vp.a+i); + for (v = st, c[0] = c[1] = 0; v < i; v++) c[vp.a[v].rev]++; + if(c[0]>1) radix_sort_scg_qts(vp.a+st, vp.a+st+c[0]); + if(c[1]>1) radix_sort_scg_qts(vp.a+st+c[0], vp.a+st+c[0]+c[1]); + st = i; + } + } + for (st = 0, i = 1; i <= vp.n; ++i){ + if (i == vp.n || vp.a[i].rev != vp.a[st].rev || + vp.a[i].qn != vp.a[st].qn || vp.a[i].tn != vp.a[st].tn || + vp.a[i].qs != vp.a[st].qs || vp.a[i].ts != vp.a[st].ts) + { + if(i - st > 1) radix_sort_scg_qte(vp.a+st, vp.a+i); + st = i; + } + } + for (st = 0, i = 1, k = 0; i <= vp.n; ++i){ + if (i == vp.n || vp.a[i].rev != vp.a[st].rev || + vp.a[i].qn != vp.a[st].qn || vp.a[i].tn != vp.a[st].tn || + vp.a[i].qs != vp.a[st].qs || vp.a[i].ts != vp.a[st].ts || + vp.a[i].qe != vp.a[st].qe || vp.a[i].te != vp.a[st].te) + { + vp.a[k] = vp.a[st]; + k++; + st = i; + } + } + ///build idx + vp.n = k; kv_resize(uint64_t, vp.idx, scg->ug->u.n); vp.idx.n = scg->ug->u.n; + memset(vp.idx.a, 0, vp.idx.n*sizeof(uint64_t)); + for (st = 0, i = 1; i <= vp.n; ++i) + { + if (i == vp.n || vp.a[i].qn != vp.a[st].qn) + { + vp.idx.a[vp.a[st].qn] = (uint64_t)st << 32 | (i - st); + st = i; + } + } + + + kv_destroy(vp); kv_destroy(vp.idx); +} +**/ +uint32_t seed_uid(ma_utg_t *vu, uint64_t* ps_idx, R_to_U* ruIndex, ma_ug_t *rug) +{ + int64_t v_i, v, w, w_i, wb, we, vb, ve, k; + uint32_t uid, is_u; + ma_utg_t *wu = NULL; + for (v_i = 0; v_i < vu->n; v_i++) { + v = vu->a[v_i]>>32; + get_R_to_U(ruIndex, v>>1, &uid, &is_u); + if(is_u == 0 || uid == (uint32_t)-1 || ps_idx[v>>1] == (uint64_t)-1) continue; + w_i = (uint32_t)ps_idx[v>>1]; + wu = &(rug->u.a[uid]); + w = wu->a[w_i]>>32; + if((v>>1)!=(w>>1)) continue; + vb = 0; ve = vu->n; ///[vb, ve) + if(v == w){ ///[wb, we) + wb = w_i - v_i; + we = wb + vu->n; + if(wb < 0 || we > wu->n) continue; + for (k = 0; k < vu->n; k++){ + if((vu->a[k+vb]>>32) != (wu->a[k+wb]>>32)) break; + } + if(k >= vu->n) return uid; + } else { + wb = w_i + 1 - (ve - v_i); + we = wb + vu->n; + if(wb < 0 || we > wu->n) continue; + for (k = 0; k < vu->n; k++){ + if((vu->a[k+vb]>>32) != ((wu->a[we-k-1]>>32)^1)) break; + } + if(k >= vu->n) return uid; + } + } + return (uint32_t)-1; +} + +void filter_ovlp_vecs(hap_overlaps_list* ha, uint32_t *a, uint32_t a_n) +{ + uint32_t st, k, i, m, v, w; + int idx; + radix_sort_ru32(a, a + a_n); + for (st = 0, m = 0, k = 1; k <= a_n; k++){ + if(k == a_n || a[k] != a[st]){ + a[m++] = a[st]; + st = k; + } + } + a_n = m; + if(a_n < 2) return; + for (k = 0; k < a_n; k++){ + v = a[k]; + for (i = k+1; i < a_n; i++) { + w = a[i]; + idx = get_specific_hap_overlap(&(ha->x[v]), v, w); + if(idx != -1) ha->x[v].a.a[idx].status = DELETE; + idx = get_specific_hap_overlap(&(ha->x[w]), w, v); + if(idx != -1) ha->x[w].a.a[idx].status = DELETE; + } + } +} + +void filter_ovlp_scg(hap_overlaps_list* ha, uint64_t* ps_idx, R_to_U* ruIndex, ma_ug_t *rug, spg_t *scg) +{ + uint32_t i, k, v, *ma = NULL, mn, luid; + ma_utg_v *pp = &(scg->ug->u); + kvec_t(uint32_t) vv; kv_init(vv); + for (i = 0; i < scg->idx.n; i++){ + ma = scg->dst.a + (scg->idx.a[i]>>32); + mn = (uint64_t)scg->idx.a[i]; + if(mn < 2) continue; + for (k = 0, vv.n = 0; k < mn; k++){ + luid = seed_uid(&(pp->a[ma[k]>>1]), ps_idx, ruIndex, rug); + if(luid == (uint32_t)-1) { + fprintf(stderr, "ERROR-scg\n"); + continue; + } + kv_push(uint32_t, vv, luid); + } + if(vv.n < 2) continue; + filter_ovlp_vecs(ha, vv.a, vv.n); + } + + for (v = 0; v < ha->num; v++){ + for (i = 0, k = 0; i < ha->x[v].a.n; i++){ + if(ha->x[v].a.a[i].status == DELETE) continue; + ha->x[v].a.a[k++] = ha->x[v].a.a[i]; + } + ha->x[v].a.n = k; + } + + kv_destroy(vv); +} + void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain, @@ -5343,7 +5612,9 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle ///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp); filter_hap_overlaps_by_length(&all_ovlp, purege_minLen); - normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex); + // normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex); + pg = init_p_g_t(ug, cov, read_g); + normalize_hap_overlaps_advance_by_p_g_t(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex, pg, cov, 0.8); if(collect_p_trans && collect_p_trans_f == 0) { @@ -5352,18 +5623,13 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans, uint32_t colle if(asm_opt.polyploidy <= 2) { - mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL); + mc_solve(&all_ovlp, cov->t_ch, NULL, ug, read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL, 1); } if(collect_p_trans && collect_p_trans_f == 1) { collect_purge_trans_cov(ug, &all_ovlp, cov, position_index); } - - pg = init_p_g_t(ug, cov, read_g); - - normalize_hap_overlaps_advance_by_p_g_t(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, - ruIndex, pg, cov, 0.8); ///normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex); ///debug_hap_overlaps(&all_ovlp, &back_all_ovlp); diff --git a/Purge_Dups.h b/Purge_Dups.h index 7b16bd2..018ab74 100644 --- a/Purge_Dups.h +++ b/Purge_Dups.h @@ -12,6 +12,7 @@ #define ALTER_COV_THRES 0.9 #define REAL_ALTER_THRES 0.25 #define CHAIN_FILTER_RATE 0.7 +#define REV_W 8 #define SELF_EXIST 0 #define REVE_EXIST 1 diff --git a/anchor.cpp b/anchor.cpp index 83e03d2..f10d607 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -76,7 +76,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg rlen = Get_READ_LENGTH(R_INF, rid); // read length // get the list of anchors - ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0); + mz1_ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0); // minimizer of queried read if (ab->mz.m > ab->old_mz_m) { ab->old_mz_m = ab->mz.m; diff --git a/hic.cpp b/hic.cpp index 309a741..fcfa8f8 100644 --- a/hic.cpp +++ b/hic.cpp @@ -2378,7 +2378,7 @@ void dfs_bubble(asg_t *g, kvec_t_u32_warp* stack, kvec_t_u32_warp* result, uint3 uint32_t get_unitig_het_arb(ma_ug_t* ug, uint32_t uid, uint8_t *r_het_flag, kv_u_trans_t *ref, uint32_t m_het_occ, uint32_t m_het_label, uint32_t p_het_label, uint32_t n_het_label) { - if(ref && u_trans_n(*ref, uid) > 0) return m_het_label; + if(ref && ref->idx.n > 0 && u_trans_n(*ref, uid) > 0) return m_het_label; ma_utg_t *u = &(ug->u.a[uid]); uint32_t k, rId; uint32_t het_occ, hom_occ; @@ -2411,7 +2411,6 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t uint64_t pathLen; bub->ug = ug; bub->b_bub = bub->b_end_bub = bub->tangle_bub = bub->cross_bub = bub->mess_bub = 0; - if(bub->round_id == 0) { buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); @@ -2421,10 +2420,8 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t bub->b_ug = NULL; kv_init(bub->chain_weight); bub->b_s_idx.n = ug->g->n_seq; memset(bub->b_s_idx.a, -1, bub->b_s_idx.n * sizeof(uint64_t)); - CALLOC(bub->index, n_vtx); for (i = 0; i < ug->g->n_seq; i++) ug->g->seq[i].c = 0; - for (v = 0; v < n_vtx; ++v) { if(ug->g->seq[v>>1].del) continue; @@ -2445,7 +2442,6 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t } } - kvec_t_u32_warp stack, result; kv_init(stack.a); kv_init(result.a); for (v = 0; v < n_vtx; ++v) @@ -2498,16 +2494,14 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t } } } - kv_push(uint32_t, bub->num, bub->list.n); free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); bub->f_bub = bub->num.n - 1; ///bub->s_bub = bub->num.n - 1; - for (i = 0; i < ug->g->n_seq; i++) - { - bub->index[i] = get_unitig_het_arb(ug, i, r_het_flag, ref, 20, M_het(*bub), P_het(*bub), (uint32_t)-1); - } - + for (i = 0; i < ug->g->n_seq; i++) + { + bub->index[i] = get_unitig_het_arb(ug, i, r_het_flag, ref, 20, M_het(*bub), P_het(*bub), (uint32_t)-1); + } for (i = 0; i < bub->f_bub; i++) { get_bubbles(bub, i, &beg, &sink, &a, &n, &pathLen); @@ -2555,7 +2549,6 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t bub->b_s_idx.a[v] |= i; } } - for (i = 0; i < ug->g->n_seq; i++) { if(bub->index[i] == M_het(*bub)) bub->index[i] = P_het(*bub); @@ -4923,7 +4916,7 @@ int load_hc_hits(kvec_pe_hit* hits, ma_ug_t* ug, const char *fn) } fclose(fp); - // fprintf(stderr, "[M::%s::] ==> Hi-C linkages have been loaded\n", __func__); + fprintf(stderr, "[M::%s::] ==> Hi-C linkages have been loaded\n", __func__); return 1; } @@ -14883,14 +14876,56 @@ inline uint32_t trans_checking_pass(bubble_type* bub, kv_u_trans_t *ref, uint32_ return 1; } +uint32_t get_u_trans_spec_idx(kv_u_trans_t *ta, uint32_t qn, uint32_t tn, u_trans_t **r_a, uint32_t *occ, uint32_t *idx) +{ + if(r_a) (*r_a) = NULL; + if(occ) (*occ) = 0; + if(idx) (*idx) = (uint32_t)-1; + u_trans_t *a = NULL; + uint32_t n, st, i; + a = u_trans_a(*ta, qn); + n = u_trans_n(*ta, qn); + for (st = 0, i = 1; i <= n; ++i) + { + if (i == n || a[i].tn != a[st].tn) + { + if(a[st].tn == tn) + { + if(r_a) (*r_a) = a + st; + if(occ) (*occ) = i - st; + if(idx) (*idx) = st + ((*ta).idx.a[(qn)]>>32); + return 1; + } + st = i; + } + } + return 0; +} + +void interpr_hit(ha_ug_index* idx, uint64_t x, uint32_t rLen, uint32_t *uid, uint32_t *beg, uint32_t *end); +int hic_sc_type(ha_ug_index* idx, kvec_pe_hit* hits, uint64_t k) +{ + uint32_t s_uid, s_beg, s_end, e_uid, e_beg, e_end, slen, elen, x = 0; + + interpr_hit(idx, hits->a.a[k].s, hits->a.a[k].len>>32, &s_uid, &s_beg, &s_end); + s_beg = (s_beg+s_end)>>1; slen = idx->ug->u.a[s_uid].len; + if(s_beg >= (slen>>1)) x+=1; + + interpr_hit(idx, hits->a.a[k].e, (uint32_t)hits->a.a[k].len, &e_uid, &e_beg, &e_end); + e_beg = (e_beg+e_end)>>1; elen = idx->ug->u.a[e_uid].len; + if(e_beg >= (elen>>1)) x+=2; + return x; +} + void weight_kv_u_trans(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, -kv_u_trans_t *ta, trans_idx* dis) +kv_u_trans_t *ta, trans_idx* dis, int sc_weight) { uint64_t k, i, shif = 64 - idx->uID_bits, beg, end, t_d; u_trans_t *e1 = NULL, *e2 = NULL; long double weight; + double *sw = NULL; u_trans_t *p = NULL; - uint32_t is_cc; + uint32_t is_cc, ii1, ii2; for (i = 0, ta->idx.n = ta->n = 0; i < link->a.n; i++) { @@ -14908,6 +14943,11 @@ kv_u_trans_t *ta, trans_idx* dis) } } kt_u_trans_t_idx(ta, idx->ug->g->n_seq); + if(sc_weight) { + k = ta->n*3; + MALLOC(sw, k); + for (i = 0; i < k; i++) sw[i] = 0; + } for (k = 0; k < hits->a.n; ++k) @@ -14925,16 +14965,38 @@ kv_u_trans_t *ta, trans_idx* dis) // if(t_d == (uint64_t)-1 && is_cc == 0) continue; // if(t_d == (uint64_t)-1 && !dis) continue; - get_u_trans_spec(ta, beg, end, &e1, NULL); - get_u_trans_spec(ta, end, beg, &e2, NULL); + get_u_trans_spec_idx(ta, beg, end, &e1, NULL, &ii1); + get_u_trans_spec_idx(ta, end, beg, &e2, NULL, &ii2); if(e1 == NULL || e2 == NULL) continue; weight = 1; if(dis) weight = get_trans_weight_advance(idx, t_d, dis); - e1->nw -= weight; e1->occ++; - e2->nw -= weight; e2->occ++; + if(sc_weight){ + i = hic_sc_type(idx, hits, k); + if(i == 0){ + e1->nw -= weight; e2->nw -= weight; + } + else{ + i--; + sw[(ii1*3)+i] -= weight; + sw[(ii2*3)+i] -= weight; + } + } else{ + e1->nw -= weight; e2->nw -= weight; + } + e1->occ++; e2->occ++; } + + if(sc_weight) { + for (i = 0; i < ta->n; ++i){ + if(ta->a[i].nw > sw[(i*3)]) ta->a[i].nw = sw[(i*3)]; + if(ta->a[i].nw > sw[(i*3)+1]) ta->a[i].nw = sw[(i*3)+1]; + if(ta->a[i].nw > sw[(i*3)+2]) ta->a[i].nw = sw[(i*3)+2]; + ta->a[i].nw *= 4/**2**/; + } + free(sw); + } } void interpr_hit(ha_ug_index* idx, uint64_t x, uint32_t rLen, uint32_t *uid, uint32_t *beg, uint32_t *end) @@ -15450,7 +15512,7 @@ ha_ug_index* idx, bubble_type* bub, int8_t *s, mc_gg_status *sa, uint32_t ignore if(hits->idx.n == 0) idx_hc_links(hits, idx, bub); - weight_kv_u_trans(idx, hits, lk, bub, ta, is_comples_weight == 1? &dis : NULL); + weight_kv_u_trans(idx, hits, lk, bub, ta, is_comples_weight == 1? &dis : NULL, asm_opt.flag&HA_F_USKEW?0:1); // if(bub->round_id == bub->n_round-1) // { // print_debug_hc_links(idx, bub, lk, ta, hits); @@ -15686,12 +15748,9 @@ void resolve_tangles_hic(ha_ug_index *idx, bubble_type *bub, kvec_pe_hit *hits, uint64_t shif = 64 - idx->uID_bits, qn, tn; pe_hit *h_a = NULL; u_trans_t *p = NULL; - identify_bubbles(idx->ug, bub, idx->t_ch->ir_het, &(idx->t_ch->k_trans)); - ta->idx.n = ta->n = 0; if(hits->idx.n == 0) idx_hc_links(hits, idx, NULL); - for (qn = 0; qn < hits->idx.n; qn++) { h_a = hits->a.a + (hits->idx.a[qn]>>32); @@ -15713,9 +15772,7 @@ void resolve_tangles_hic(ha_ug_index *idx, bubble_type *bub, kvec_pe_hit *hits, } } } - radix_sort_u_trans_m(ta->a, ta->a + ta->n); - for (k = 1, l = 0, m = 0; k <= ta->n; ++k) { if (k == ta->n || ta->a[l].qn != ta->a[k].qn || ta->a[l].tn != ta->a[k].tn) //same qn and tn @@ -15736,23 +15793,21 @@ void resolve_tangles_hic(ha_ug_index *idx, bubble_type *bub, kvec_pe_hit *hits, } ta->n = m; kt_u_trans_t_idx(ta, idx->ug->g->n_seq); - resolve_bubble_chain_by_hic(idx, ta, bub); - ta->idx.n = ta->n = 0; } -void print_kv_u_trans_t(kv_u_trans_t *ta) +void print_kv_u_trans_t(kv_u_trans_t *ta, ma_ug_t* ug) { uint32_t i; u_trans_t *p = NULL; for (i = 0; i < ta->n; i++) { p = &(ta->a[i]); - fprintf(stderr, "q-utg%.6ul\tqs(%u)\tqe(%u)\tt-utg%.6ul\tts(%u)\tte(%u)\trev(%u)\tw(%f)\tf(%u)\n", - p->qn+1, p->qs, p->qe, p->tn+1, p->ts, p->te, p->rev, p->nw, p->f); + fprintf(stderr, "q-utg%.6ul\tql(%u)\tqs(%u)\tqe(%u)\tt-utg%.6ul\ttl(%u)\tts(%u)\tte(%u)\trev(%u)\tw(%f)\tf(%u)\n", + p->qn+1, ug->u.a[p->qn].len, p->qs, p->qe, p->tn+1, ug->u.a[p->tn].len, p->ts, p->te, p->rev, p->nw, p->f); } fprintf(stderr, "[M::%s::] \n", __func__); } @@ -16191,6 +16246,61 @@ void renew_idx_para(ha_ug_index* idx, ma_ug_t* ug) idx->rev_mode = ((uint64_t)1) << 63; } +uint32_t get_oe_occ(uint32_t qn, uint32_t tn, kvec_pe_hit* hits, ha_ug_index* idx) +{ + uint64_t shif = 64 - idx->uID_bits, occ = 0; + pe_hit *h_a = hits->a.a + (hits->idx.a[qn]>>32); + uint32_t h_occ = (uint32_t)(hits->idx.a[qn]), i; + for (i = 0; i < h_occ; i++) { + if(((h_a[i].s<<1)>>shif)==qn && ((h_a[i].e<<1)>>shif)==tn) occ++; + } + return occ; +} + +void optimize_u_trans(kv_u_trans_t *ovlp, kvec_pe_hit* hits, ha_ug_index* idx) +{ + if(hits->idx.n == 0) idx_hc_links(hits, idx, NULL); + uint64_t i, m, occ; + u_trans_t *x = NULL, *p = NULL; + kv_u_trans_t k_trans; + kv_init(k_trans); kv_init(k_trans.idx); + for (i = 0; i < ovlp->n; i++){ + x = &(ovlp->a[i]); + if(x->qn > x->tn) continue; + if(x->f != RC_2 || x->del) continue; + occ = get_oe_occ(x->qn, x->tn, hits, idx) + get_oe_occ(x->tn, x->qn, hits, idx); + kv_pushp(u_trans_t, k_trans, &p); + (*p) = (*x); p->nw = (x->nw*(1-(((double)(occ<<1))/((double)(hits->occ.a[x->qn]+hits->occ.a[x->tn]))))); + if(p->nw < 0) fprintf(stderr, "ERROR-nw\n"); + if(p->nw == 0) p->nw = x->nw*0.005; + if(p->nw == 0) { + k_trans.n--; + } + else { + kv_pushp(u_trans_t, k_trans, &p); + (*p) = k_trans.a[k_trans.n-2]; + p->qn = k_trans.a[k_trans.n-2].tn; p->qs = k_trans.a[k_trans.n-2].ts; p->qe = k_trans.a[k_trans.n-2].te; + p->tn = k_trans.a[k_trans.n-2].qn; p->ts = k_trans.a[k_trans.n-2].qs; p->te = k_trans.a[k_trans.n-2].qe; + } + } + kt_u_trans_t_idx(&k_trans, idx->ug->g->n_seq); + mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, 1, NULL, 1, NULL, NULL, 1); + for (i = m = 0; i < ovlp->n; i++){ + x = &(ovlp->a[i]); + if(x->del) continue; + if(x->f == RC_2){ + get_u_trans_spec(&k_trans, x->qn, x->tn, &p, NULL); + if(!p) continue; + } + ovlp->a[m++] = ovlp->a[i]; + } + ovlp->n = m; + kt_u_trans_t_idx(ovlp, idx->ug->g->n_seq); + free(hits->idx.a); hits->idx.a = NULL; hits->idx.n = hits->idx.m = 0; + free(hits->occ.a); hits->occ.a = NULL; hits->occ.n = hits->occ.m = 0; + kv_destroy(k_trans); kv_destroy(k_trans.idx); +} + int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_opt_t *opt, kvec_pe_hit **rhits) { double index_time = yak_realtime(); @@ -16209,15 +16319,16 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o alignment_worker_pipeline(&sl, fn1, fn2); write_hc_hits(&sl.hits, idx->ug, asm_opt.output_file_name); } + sl.hits.uID_bits = idx->uID_bits; sl.hits.pos_mode = idx->pos_mode; + optimize_u_trans(&(idx->t_ch->k_trans), &sl.hits, idx); filter_kv_u_trans_t(&(idx->t_ch->k_trans), idx->ug, 0.5); + // print_kv_u_trans_t(&(idx->t_ch->k_trans), idx->ug); // flter_by_cov(idx, &sl.hits, 2); // update_hits(idx, &sl.hits, idx->t_ch->is_r_het); ///debug_hc_hits_v14(&sl.hits, asm_opt.output_file_name, sl.idx); ////dedup_hits(&(sl.hits), sl.idx); ///write_hc_hits_v14(&sl.hits, asm_opt.output_file_name); - - sl.hits.uID_bits = idx->uID_bits; sl.hits.pos_mode = idx->pos_mode; if(asm_opt.misjoin_len > 0) { update_switch_unitig(idx->ug, idx->read_g, &(sl.hits), &(idx->t_ch->k_trans), 10, 20, asm_opt.misjoin_len, 0.15); @@ -16255,7 +16366,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o // update_trans_g(idx, &k_trans, &bub); /*******************************for debug************************************/ mc_solve(NULL, NULL, &k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag, - (bub.round_id == 0? 1 : 0), s->s, 1, /**&bub**/NULL, &(idx->t_ch->k_trans)); + (bub.round_id == 0? 1 : 0), s->s, 1, /**&bub**/NULL, &(idx->t_ch->k_trans), 0); /*******************************for debug************************************/ label_unitigs_sm(s->s, NULL, idx->ug); @@ -16321,6 +16432,49 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o return 1; } +spg_t *hic_short_pre_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_opt_t *opt, kvec_pe_hit **rhits) +{ + // double index_time = yak_realtime(); + sldat_t sl; + sl.idx = idx; + sl.t_ch = idx->t_ch; + sl.chunk_size = 20000000; + sl.n_thread = asm_opt.thread_num; + sl.total_base = sl.total_pair = 0; + idx->hap_cnt = asm_opt.hap_occ; + kv_init(sl.hits.a); kv_init(sl.hits.idx); kv_init(sl.hits.occ); + + if(!load_hc_hits(&sl.hits, idx->ug, asm_opt.output_file_name)) + { + alignment_worker_pipeline(&sl, fn1, fn2); + fprintf(stderr, "sb0sb\n"); + write_hc_hits(&sl.hits, idx->ug, asm_opt.output_file_name); + fprintf(stderr, "sb1sb\n"); + } + sl.hits.uID_bits = idx->uID_bits; sl.hits.pos_mode = idx->pos_mode; + kv_u_trans_t k_trans; + kv_init(k_trans); kv_init(k_trans.idx); + bubble_type bub; + memset(&bub, 0, sizeof(bubble_type)); + bub.round_id = 0; bub.n_round = asm_opt.n_weight; + fprintf(stderr, "sb2sb\n"); + resolve_tangles_hic(idx, &bub, &sl.hits, &k_trans); + fprintf(stderr, "sb3sb\n"); + spg_t *scg = horder_utg(&sl.hits, idx->uID_bits, idx->pos_mode, idx->read_g, idx->ug, &bub, opt); + fprintf(stderr, "sb4sb\n"); + if(rhits){ + CALLOC(*rhits, 1); + (**rhits) = sl.hits; + sl.hits.a.a = NULL; + sl.hits.idx.a = NULL; + sl.hits.occ.a = NULL; + } + kv_destroy(sl.hits.a); kv_destroy(sl.hits.idx); kv_destroy(sl.hits.occ); + kv_destroy(k_trans); kv_destroy(k_trans.idx); + destory_bubbles(&bub); + return scg; +} + int load_psg_t(psg_t **sg, const char *fn) { uint64_t flag = 0; @@ -16433,7 +16587,7 @@ int hic_short_align_poy(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, // } /*******************************for debug************************************/ } - write_psg_t(s, asm_opt.output_file_name); + if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_psg_t(s, asm_opt.output_file_name); skip_flipping: verbose_het_stat(&bub); @@ -16488,6 +16642,19 @@ void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, destory_hc_pt_index(ug_index);free(ug_index); } +spg_t *hic_pre_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, kvec_pe_hit **rhits) +{ + ug_index = NULL; + int exist = (asm_opt.load_index_from_disk? + load_hc_pt_index(&ug_index, ug, asm_opt.output_file_name) : 0); + if(exist == 0) ug_index = build_unitig_index(ug, asm_opt.hic_mer_length, asm_opt.hap_occ, 0, asm_opt.thread_num); + if(exist == 0) write_hc_pt_index(ug_index, asm_opt.output_file_name); + ug_index->ug = ug; + ug_index->read_g = read_g; + ug_index->t_ch = t_ch; + return hic_short_pre_align(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt, rhits); +} + void init_ug_idx(ma_ug_t *ug, uint64_t k, uint64_t up_bound, uint64_t low_bound, uint64_t build_idx) { diff --git a/hic.h b/hic.h index f0ebbc9..c0f8c82 100644 --- a/hic.h +++ b/hic.h @@ -75,6 +75,7 @@ typedef struct{ }pdq; + #define P_het(B) ((B).num.n) #define M_het(B) ((B).num.n + 1) // #define IF_BUB(ID, B) ((B).index[(ID)] < (B).num.n) @@ -106,5 +107,5 @@ long long *dis); void set_utg_by_dis(uint32_t v, pdq* pq, asg_t *g, kvec_t_u32_warp *res, uint32_t dis); void dedup_hits(kvec_pe_hit* hits, uint64_t is_dup); void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, uint32_t is_poy, kvec_pe_hit **rhits); - +spg_t *hic_pre_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, kvec_pe_hit **rhits); #endif diff --git a/horder.cpp b/horder.cpp index 1313c61..53761f5 100644 --- a/horder.cpp +++ b/horder.cpp @@ -2178,9 +2178,9 @@ uint64_t rs, uint64_t re, uint64_t limit_s, uint64_t limit_e, int unique_only) return cnt; } -void detect_lowNs(kvec_pe_hit *hit, uint64_t sHit, uint64_t eHit, kvec_t_u64_warp *b, -h_cov_t *Np, uint64_t len, uint64_t cutoff_s, uint64_t cutoff_e, h_covs *res, -h_covs *cov_buf, h_covs *b_points, uint64_t local_bound, int unique_only) +int detect_lowNs(kvec_pe_hit *hit, uint64_t sHit, uint64_t eHit, kvec_t_u64_warp *b, +h_cov_t *Np, uint64_t len, uint64_t cutoff_s, uint64_t cutoff_e, uint64_t force_cutoff, uint64_t force_cutoff_cov, +h_covs *res, h_covs *cov_buf, h_covs *b_points, uint64_t local_bound, int unique_only) { uint64_t cov_hic, cov_utg, cov_ava, i, p0s, p0e, p1s, p1e, span_s, span_e, cutoff, bs, be, occ = 0; uint64_t sPos, ePos, min_cutoff; @@ -2219,6 +2219,18 @@ h_covs *cov_buf, h_covs *b_points, uint64_t local_bound, int unique_only) radix_sort_ho64(b->a.a, b->a.a+b->a.n); cov_utg = get_hic_cov_interval(b->a.a, b->a.n, 1, NULL, NULL, NULL); cov_ava = (cov_utg? cov_hic/cov_utg:0); + if(force_cutoff != (uint64_t)-1 || force_cutoff_cov != (uint64_t)-1) + { + min_cutoff = get_sub_cov(hit, sHit, eHit, len, Np->s, Np->e, sPos, ePos, unique_only); + if((force_cutoff != (uint64_t)-1 && min_cutoff <= (cov_ava/force_cutoff)) || + (force_cutoff_cov != (uint64_t)-1 && min_cutoff <= force_cutoff_cov)) + { + kv_pushp(h_cov_t, *b_points, &p); + p->s = get_hit_suid(*hit, sHit); p->e = Np->dp; p->dp = 0; + return 1; + } + } + ///if cov_ava == 0, do nothing or break? /*******************************for debug************************************/ // fprintf(stderr, "\n[M::%s::] utg%.6lul, ulen: %lu, # hic hits: %lu, map cov: %lu, utg cov: %lu, average: %lu\n", @@ -2269,24 +2281,29 @@ h_covs *cov_buf, h_covs *b_points, uint64_t local_bound, int unique_only) { kv_pushp(h_cov_t, *b_points, &p); p->s = get_hit_suid(*hit, sHit); p->e = Np->dp; p->dp = 0; + return 1; /*******************************for debug************************************/ // fprintf(stderr, "consensus_break-rid: %lu\n", p->e); /*******************************for debug************************************/ } } + return 0; } -uint64_t break_scaffold(horder_t *h, uint64_t cutoff_s, uint64_t cutoff_e, uint64_t local_bound, int unique_only) +uint64_t break_scaffold(horder_t *h, uint64_t cutoff_s, uint64_t cutoff_e, uint64_t force_cutoff, uint64_t force_cutoff_cov, +uint64_t local_bound, int unique_only, h_covs *r_b_points) { uint64_t k, l, i, ulen; kvec_t_u64_warp b; kv_init(b.a); h_covs cov_buf; kv_init(cov_buf); h_covs res; kv_init(res); - h_covs b_points; kv_init(b_points); + h_covs *b_points = NULL; + if(r_b_points) b_points = r_b_points; + else CALLOC(b_points, 1); + b_points->n = 0; h_covs Ns; kv_init(Ns); ma_ug_t *ug = h->ug; kvec_pe_hit *hits = &(h->u_hits); - b_points.n = 0; for (k = 1, l = 0; k <= hits->a.n; ++k) { if (k == hits->a.n || (get_hit_suid(*hits, k) != get_hit_suid(*hits, l))) @@ -2302,8 +2319,10 @@ uint64_t break_scaffold(horder_t *h, uint64_t cutoff_s, uint64_t cutoff_e, uint6 { for (i = 0; i < Ns.n; i++) { - detect_lowNs(hits, l, k, &b, &(Ns.a[i]), ulen, cutoff_s, cutoff_e, - &res, &cov_buf, &b_points, local_bound, unique_only); + if(detect_lowNs(hits, l, k, &b, &(Ns.a[i]), ulen, cutoff_s, cutoff_e, force_cutoff, + force_cutoff_cov, &res, &cov_buf, b_points, local_bound, unique_only) && r_b_points){ + b_points->a[b_points->n-1].dp = i; + } } } } @@ -2311,15 +2330,68 @@ uint64_t break_scaffold(horder_t *h, uint64_t cutoff_s, uint64_t cutoff_e, uint6 } } - break_utg_horder(h, &b_points); - + if(!r_b_points){ + break_utg_horder(h, b_points); + kv_destroy(*b_points); + } + kv_destroy(b.a); kv_destroy(cov_buf); kv_destroy(res); - kv_destroy(b_points); kv_destroy(Ns); - return b_points.n; + return b_points->n; +} + +uint64_t break_scaffold_mean(horder_t *h, uint64_t cutoff_s, uint64_t cutoff_e, uint64_t force_cutoff, uint64_t force_cutoff_cov, +uint64_t local_bound, int unique_only, h_covs *r_b_points) +{ + uint64_t k, l, i, ulen; + kvec_t_u64_warp b; kv_init(b.a); + h_covs cov_buf; kv_init(cov_buf); + h_covs res; kv_init(res); + h_covs *b_points = NULL; + if(r_b_points) b_points = r_b_points; + else CALLOC(b_points, 1); + b_points->n = 0; + h_covs Ns; kv_init(Ns); + ma_ug_t *ug = h->ug; + kvec_pe_hit *hits = &(h->u_hits); + for (k = 1, l = 0; k <= hits->a.n; ++k) + { + if (k == hits->a.n || (get_hit_suid(*hits, k) != get_hit_suid(*hits, l))) + { + ulen = ug->u.a[get_hit_suid(*hits, l)].len; + Ns.n = 0; + // if(ulen >= BREAK_THRES) + { + get_Ns(&(ug->u.a[get_hit_suid(*hits, l)]), &Ns); + if(Ns.n) + { + for (i = 0; i < Ns.n; i++) + { + if(detect_lowNs(hits, l, k, &b, &(Ns.a[i]), ulen, cutoff_s, cutoff_e, force_cutoff, + force_cutoff_cov, &res, &cov_buf, b_points, local_bound, unique_only) && r_b_points){ + b_points->a[b_points->n-1].dp = i; + } + } + } + } + l = k; + } + } + + if(!r_b_points){ + break_utg_horder(h, b_points); + kv_destroy(*b_points); + } + + kv_destroy(b.a); + kv_destroy(cov_buf); + kv_destroy(res); + kv_destroy(Ns); + + return b_points->n; } @@ -2729,7 +2801,7 @@ void update_scg(horder_t *h, trans_col_t *t_idx) h->sg.g->seq[i].ez[0] = ug->u.a[i].len>>1; h->sg.g->seq[i].ez[1] = ug->u.a[i].len - (ug->u.a[i].len>>1); } - idx = build_interval_idx(hits, ug); + idx = build_interval_idx(hits, ug);///idx is used to get density for (i = 0, e.n = 0; i < hits->a.n; i++) { @@ -3397,10 +3469,12 @@ uint32_t get_sl_occ(sc_lay_t *sl) void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis) { - uint32_t k, m/**, max_utg, max_sc**/; + uint32_t k/**m, max_utg, max_sc**/; lay_t *p = NULL; uint8_t *sgv = NULL; MALLOC(sgv, sl->n); double *w = NULL; MALLOC(w, sl->n); + + /** uint32_t *idx = NULL; MALLOC(idx, h->sg.g->n_seq); memset(idx, -1, sizeof(uint32_t)*h->sg.g->n_seq); @@ -3412,8 +3486,6 @@ void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis) idx[p->a[m]>>1] = k; } } - - /** while (get_max_anchor(h, sl, vis, w, sgv, idx, &max_utg, &max_sc)) { @@ -3422,6 +3494,7 @@ void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis) vis[max_utg<<1] = vis[(max_utg<<1)+1] = 1; idx[max_utg] = max_sc; } + free(idx); **/ for (k = 0; k < h->sg.g->n_seq; k++) @@ -3434,8 +3507,7 @@ void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis) vis[(k<<1)] = vis[(k<<1)+1] = 1; } - - free(w); free(idx); free(sgv); + free(w); free(sgv); } @@ -3591,7 +3663,7 @@ void update_avoids(horder_t *h, sc_lay_t *sl) h->avoid.n = m; } -void update_ug_by_layout(horder_t *h, sc_lay_t *sl) +void update_ug_by_layout(horder_t *h, sc_lay_t *sl, ma_ug_t* i_ug) { uint32_t i; lay_t *p = NULL; @@ -3603,7 +3675,7 @@ void update_ug_by_layout(horder_t *h, sc_lay_t *sl) { p = &(sl->a[i]); kv_pushp(ma_utg_t, sug->u, &pu); - generate_scaffold(pu, p, h->ug, h->r_g); + generate_scaffold(pu, p, i_ug?i_ug:h->ug, h->r_g); } ma_ug_destroy(h->ug); h->ug = sug; @@ -3677,7 +3749,14 @@ void get_long_switch_scaffolds(horder_t *h, sc_lay_t *sl, osg_t *lg) kv_destroy(idx); } -void layout_scg(horder_t *h, double nw_thres, uint32_t occ_thres) +void destory_sc_lay_t(sc_lay_t *sl) +{ + uint32_t i; + for (i = 0; i < sl->n; i++) kv_destroy(sl->a[i]); + kv_destroy(*sl); +} + +void layout_scg(horder_t *h, double nw_thres, uint32_t occ_thres, sc_lay_t *r_sl) { uint32_t k; osg_arc_t *p = NULL, *lp = NULL; @@ -3708,17 +3787,23 @@ void layout_scg(horder_t *h, double nw_thres, uint32_t occ_thres) get_backbone_layout(h, &sl, lg, vis); - get_long_switch_scaffolds(h, &sl, lg); + // get_long_switch_scaffolds(h, &sl, lg); refine_layout(h, &sl, vis); // print_N50_layout(h->ug, &sl); - update_ug_by_layout(h, &sl); + update_ug_by_layout(h, &sl, NULL); print_N50(h->ug); - kv_destroy(sl); + if(r_sl){ + r_sl->a = sl.a; sl.a = NULL; + r_sl->m = sl.m; sl.m = 0; + r_sl->n = sl.n; sl.n = 0; + + } + destory_sc_lay_t(&sl); osg_destroy(lg); free(vis); } @@ -3729,7 +3814,7 @@ void renew_scaffold(horder_t *h) while (1) { update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); - if(!break_scaffold(h, 5, 15, 2500000, 1)) break; + if(!break_scaffold(h, 5, 15, (uint64_t)-1, (uint64_t)-1, 2500000, 1, NULL)) break; print_N50(h->ug); } fprintf(stderr, "[M::%s::%.3f] \n", __func__, yak_realtime()-index_time); @@ -3789,7 +3874,7 @@ void scaffold_hap(horder_t *h, ug_opt_t *opt, trans_col_t *t_idx, uint32_t round { update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); update_scg(h, t_idx); - layout_scg(h, 1.001, 19); + layout_scg(h, 1.001, 19, NULL); renew_scaffold(h); } @@ -3832,7 +3917,7 @@ void scaffold_ug(horder_t *h, ma_ug_t *ug, ug_opt_t *opt, uint32_t round, char * fprintf(stderr, "[M::%s::]**i->%u**\n", __func__, i); update_scg(h, NULL); fprintf(stderr, "[M::%s::]***i->%u***\n", __func__, i); - layout_scg(h, 1.001, 19); + layout_scg(h, 1.001, 19, NULL); fprintf(stderr, "[M::%s::]****i->%u****\n", __func__, i); renew_scaffold(h); fprintf(stderr, "[M::%s::]*****i->%u*****\n", __func__, i); @@ -3887,7 +3972,7 @@ asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, kv_u_trans_t *ref, ug_opt_t *opt, { update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); update_scg(h, t_idx); - layout_scg(h, 1.001, 19); + layout_scg(h, 1.001, 19, NULL); renew_scaffold(h); } @@ -3902,6 +3987,113 @@ asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, kv_u_trans_t *ref, ug_opt_t *opt, return h; } +void cpy_u_hits(kvec_pe_hit *u_hits, kvec_pe_hit *i_hits, uint32_t u_n) +{ + uint64_t i; + memset(u_hits, 0, sizeof(kvec_pe_hit)); + u_hits->pos_mode = i_hits->pos_mode; + u_hits->uID_bits = i_hits->uID_bits; + u_hits->a.n = u_hits->a.m = i_hits->a.n; + MALLOC(u_hits->a.a, u_hits->a.n); + memcpy(u_hits->a.a, i_hits->a.a, u_hits->a.n*sizeof(pe_hit)); + for (i = 0; i < u_hits->a.n; i++) u_hits->a.a[i].id = 1; + idx_hits(u_hits, u_n); +} + +void update_sc_lay(sc_lay_t *sl, h_covs *b) +{ + uint64_t k, l, i, pidx, cidx; + lay_t *s = NULL; + lay_t *p = NULL; + for (k = 1, l = 0; k <= b->n; ++k) + { + if (k == b->n || (b->a[k].s != b->a[l].s)) + { + s = &(sl->a[b->a[l].s]); + for (i = l, pidx = 0; i < k; i++){ + cidx = b->a[i].dp; + kv_pushp(lay_t, *sl, &p); + kv_init(*p); + p->n = p->m = (cidx - pidx + 1)<<1; + MALLOC(p->a, p->n); + mempcpy(p->a, s->a + (pidx<<1), sizeof(*(p->a))*p->n); + pidx = cidx + 1; + } + + if(pidx >= (s->n>>1)) fprintf(stderr, "ERROR-update\n"); + cidx = (s->n>>1)-1; + kv_pushp(lay_t, *sl, &p); + kv_init(*p); + p->n = p->m = (cidx - pidx + 1)<<1; + MALLOC(p->a, p->n); + mempcpy(p->a, s->a + (pidx<<1), sizeof(*(p->a))*p->n); + free(s->a); s->n = s->m = 0; + l = k; + } + } + + for (i = k = 0; i < sl->n; i++){ + if(!sl->a[i].a) continue; + sl->a[k] = sl->a[i]; + sl->a[i].a = NULL; sl->a[i].n = sl->a[i].m = 0; + if(sl->a[k].n == 2){ + sl->a[k].a[0] >>= 1; sl->a[k].a[0] <<= 1; + sl->a[k].a[1] >>= 1; sl->a[k].a[1] <<= 1; sl->a[k].a[1]++; + } + k++; + } + sl->n = k; +} + +void renew_scaffold_utg(horder_t *h, sc_lay_t *sl, ma_ug_t* i_ug) +{ + double index_time = yak_realtime(); + h_covs b; kv_init(b); + while (1) + { + update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); + if(!break_scaffold(h, 5, 15, 15, 25, 2500000, 1, &b)) break; + update_sc_lay(sl, &b); + update_ug_by_layout(h, sl, i_ug); + print_N50(h->ug); + } + fprintf(stderr, "[M::%s::%.3f] \n", __func__, yak_realtime()-index_time); + kv_destroy(b); +} + +spg_t *scf_g(sc_lay_t *sl, ma_ug_t* ug) +{ + spg_t *scg = NULL; CALLOC(scg, 1); scg->ug = ug; + uint64_t i, k; + lay_t *p = NULL; + for (i = 0; i < sl->n; i++){ + p = &(sl->a[i]); + kv_push(uint64_t, scg->idx, (uint64_t)scg->dst.n << 32 | (p->n>>1)); + for (k = 0; k < p->n; k+=2) kv_push(uint32_t, scg->dst, p->a[k]); + } + return scg; +} + +spg_t *horder_utg(kvec_pe_hit *i_hits, uint64_t i_hits_uid_bits, uint64_t i_hits_pos_mode, +asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt) +{ + horder_t *h = NULL; CALLOC(h, 1); + sc_lay_t sl; kv_init(sl); + get_r_hits(i_hits, &(h->r_hits), i_rg, i_ug, bub, i_hits_uid_bits, i_hits_pos_mode); + h->r_g = copy_read_graph(i_rg); + horder_clean_sg_by_utg(h->r_g, i_ug); + h->ug = copy_untig_graph(i_ug); asg_destroy(h->ug->g); h->ug->g = NULL; + cpy_u_hits(&(h->u_hits), i_hits, h->ug->u.n); + + update_scg(h, NULL); + layout_scg(h, ((double)1)/((double)0.75), 19, &sl); + renew_scaffold_utg(h, &sl, i_ug); + + spg_t *scg = scf_g(&sl, i_ug); + destory_sc_lay_t(&sl); + destory_horder_t(&h); + return scg; +} void ha_aware_order(kvec_pe_hit *r_hits, asg_t *rg, ma_ug_t *ug_fa, ma_ug_t *ug_mo, kv_u_trans_t *ref, diff --git a/horder.h b/horder.h index e262075..8e136a9 100644 --- a/horder.h +++ b/horder.h @@ -67,4 +67,6 @@ kvec_pe_hit *get_r_hits_order(kvec_pe_hit *uhits, uint64_t hits_uid_bits, uint64 asg_t *rg, ma_ug_t* ug, bubble_type* bub); void ha_aware_order(kvec_pe_hit *r_hits, asg_t *rg, ma_ug_t *ug_fa, ma_ug_t *ug_mo, kv_u_trans_t *ref, ug_opt_t *opt, uint32_t round); +spg_t *horder_utg(kvec_pe_hit *i_hits, uint64_t i_hits_uid_bits, uint64_t i_hits_pos_mode, +asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt); #endif diff --git a/htab.cpp b/htab.cpp index 516c971..d6bec2a 100644 --- a/htab.cpp +++ b/htab.cpp @@ -38,75 +38,6 @@ void *ha_flt_tab_hp; ha_pt_t *ha_idx_hp; void *ha_ct_table; -#define MZ_FUNC_INIT(sf, HType) \ -static inline void sf##_init_kuf(pl_data_t *p, st_data_t *s){\ - int i, n_pre = 1<opt->pre, m;\ - /**allocate the k-mer buffer**/\ - CALLOC(s->buf, n_pre);\ - m = (int)(s->nk * 1.2 / n_pre) + 1;\ - /**pre-allocate memory for each of 4096 buffer**/\ - for (i = 0; i < n_pre; ++i) {\ - s->buf[i].m = m;\ - /**for 0-th counting, p->pt = NULL**/\ - if (p->pt && !(p->flag&HAF_COUNT_REFINE)) MALLOC(s->buf[i].b_##sf, m);\ - else MALLOC(s->buf[i].a, m);\ - }\ -}\ -static inline void sf##_destory_kuf(pl_data_t *p, st_data_t *s, int n){\ - int i;\ - uint64_t n_ins = 0;\ - /**n_ins is number of distinct k-mers**/\ - for (i = 0; i < n; ++i) {\ - n_ins += s->buf[i].n_ins;\ - if (p->pt && !(p->flag&HAF_COUNT_REFINE)) free(s->buf[i].b_##sf);\ - else free(s->buf[i].a);\ - }\ - if (p->ct) p->ct->tot += n_ins, p->ct->bs += s->sum_len;\ - if (p->pt) p->pt->tot_pos += n_ins;\ - free(s->buf);\ - /**#if 0\ - fprintf(stderr, "[M::%s::%.3f*%.2f] processed %ld sequences; %ld %s in the hash table\n", __func__,\ - yak_realtime(), yak_cpu_usage(), (long)s->n_seq0 + s->n_seq,\ - (long)(p->pt? p->pt->tot_pos : p->ct->tot), p->pt? "positions" : "distinct k-mers");\ - #endif**/\ - free(s);\ -}\ -static inline void sf##_pt_insert_buf(ch_buf_t *buf, int p, const HType *y){\ - /**assign minimizer to one of 4096 bins by low 12 bits**/\ - int pre = y->x & ((1<n == b->m) {\ - b->m = b->m < 8? 8 : b->m + (b->m>>1);\ - REALLOC(b->b_##sf, b->m);\ - }\ - b->b_##sf[b->n++] = *y;\ -}\ -static inline void sf##_mselect(pl_data_t *p, st_data_t *s){\ - int i; uint32_t j;\ - /**s->n_seq is how many reads at this buffer**/\ - /**s->mz && s->mz_buf are lists of minimzer vectors**/\ - CALLOC(s->sf, s->n_seq), CALLOC(s->sf##_buf, p->opt->n_thread), CALLOC(s->mt, p->opt->n_thread);\ - /**calculate minimzers for each read, each read corresponds to one thread**/\ - kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq);\ - for (i = 0; i < p->opt->n_thread; ++i) free(s->mt[i].a), free(s->sf##_buf[i].a);\ - free(s->mt), free(s->sf##_buf);\ - /**insert minimizers**/\ - if (p->pt && !(p->flag&HAF_COUNT_REFINE)) {/**insert whole minimizer**/\ - for (i = 0; i < s->n_seq; ++i)\ - for (j = 0; j < s->sf[i].n; ++j)\ - sf##_pt_insert_buf(s->buf, p->opt->pre, &s->sf[i].a[j]);\ - } else {/**just insert the hash key of minimizer**/\ - for (i = 0; i < s->n_seq; ++i)\ - for (j = 0; j < s->sf[i].n; ++j)\ - ct_insert_buf(s->buf, p->opt->pre, s->sf[i].a[j].x);\ - }\ - for (i = 0; i < s->n_seq; ++i) {\ - p->n_mz += s->sf[i].n;\ - free(s->sf[i].a);\ - if (!p->is_store) free(s->seq[i]);\ - }\ - free(s->sf);} - /*************************** * Yak specific parameters * ***************************/ @@ -601,77 +532,7 @@ const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash) /********************************** * Buffer for counting all k-mers * **********************************/ - -typedef struct { - int n, m; - uint64_t n_ins; - uint64_t *a; - ha_mz1_t *b_mz; - ha_mzl_t *b_mzl; -} ch_buf_t; - -///p = 12 -static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k-mer $y to a linear buffer -{ - ///assign k-mer to one of the 4096 bins - ///using low 12 bits for assigning - ///so all elements at b have the same low 12 bits - int pre = y & ((1<n == b->m) { - b->m = b->m < 8? 8 : b->m + (b->m>>1); - REALLOC(b->a, b->m); - } - b->a[b->n++] = y; -} - -///buf is the read block, k is the k-mer length, p = 12, len is the read length, seq is the read -static void count_seq_buf(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf -{ - int i, l; - uint64_t x[4], mask = (1ULL<>1)) & mask; - x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; - x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; - if (++l >= k) - ct_insert_buf(buf, p, yak_hash_long(x)); - } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart - } -} - -static void count_seq_buf_HPC(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf -{ - int i, l, last = -1; - uint64_t x[4], mask = (1ULL<>1)) & mask; - x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; - x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; - if (++l >= k) - ct_insert_buf(buf, p, yak_hash_long(x)); - last = c; - } - } else l = 0, last = -1, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart - } -} - -/****************** - * K-mer counting * - ******************/ - KSEQ_INIT(gzFile, gzread) - #define HAF_COUNT_EXACT 0x1 #define HAF_COUNT_ALL 0x2 #define HAF_RS_WRITE_LEN 0x4 @@ -696,178 +557,285 @@ typedef struct { // global data structure for kt_pipeline() const ma_utg_v *us_in; } pl_data_t; -typedef struct { // data structure for each step in kt_pipeline() - pl_data_t *p; - uint64_t n_seq0; ///the start index of current buffer block at R_INF - ///sum_len = total bases, nk = number of k-mers - int n_seq, m_seq, sum_len, nk, uq; - int *len; - char **seq; - ha_mz1_v *mz_buf; - ha_mz1_v *mz; - ha_mzl_v *mzl_buf; - ha_mzl_v *mzl; - ch_buf_t *buf; - st_mt_t *mt; -} st_data_t; - -static void worker_for_insert(void *data, long i, int tid) // callback for kt_for() -{ - st_data_t *s = (st_data_t*)data; - ch_buf_t *b = &s->buf[i]; - if (s->p->pt) - { - if(s->p->flag&HAF_COUNT_REFINE) b->n_ins += ha_pt_cnt_insert_list(s->p->pt, b->n, b->a); - else b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b_mz); - } - else///for 0-th count, go into here - { - b->n_ins += ha_ct_insert_list(s->p->ct, s->p->create_new, b->n, b->a); - } +#define MZ_TEST_INIT(sf, HType, VType, IType, Ia) \ +typedef struct {int n, m; uint64_t n_ins; uint64_t *a; HType *b;} sf##_ch_buf_t;\ +static inline void sf##_ct_insert_buf(sf##_ch_buf_t *buf, int p, uint64_t y) /** insert a k-mer $y to a linear buffer**/\ +{\ + /**assign k-mer to one of the 4096 bins**/\ + /**using low 12 bits for assigning**/\ + /**so all elements at b have the same low 12 bits**/\ + int pre = y & ((1<n == b->m) {\ + b->m = b->m < 8? 8 : b->m + (b->m>>1);\ + REALLOC(b->a, b->m);\ + }\ + b->a[b->n++] = y;\ +}\ +/**buf is the read block, k is the k-mer length, p = 12, len is the read length, seq is the read**/\ +static void sf##_count_seq_buf(sf##_ch_buf_t *buf, int k, int p, int len, const char *seq) /**insert k-mers in $seq to linear buffer $buf**/\ +{\ + int i, l;\ + uint64_t x[4], mask = (1ULL<>1)) & mask;\ + x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift;\ + x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift;\ + if (++l >= k)\ + sf##_ct_insert_buf(buf, p, yak_hash_long(x));\ + } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; /** if there is an "N", restart**/\ + }\ +}\ +static void sf##_count_seq_buf_HPC(sf##_ch_buf_t *buf, int k, int p, int len, const char *seq) /**insert k-mers in $seq to linear buffer $buf**/\ +{\ + int i, l, last = -1;\ + uint64_t x[4], mask = (1ULL<>1)) & mask;\ + x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift;\ + x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift;\ + if (++l >= k)\ + sf##_ct_insert_buf(buf, p, yak_hash_long(x));\ + last = c;\ + }\ + } else l = 0, last = -1, x[0] = x[1] = x[2] = x[3] = 0; /**if there is an "N", restart**/\ + }\ +}\ +int sf##_ha_pt_insert_list(ha_pt_t *h, int n, const HType *a)\ +{\ + int j, mask = (1<pre) - 1, n_ins = 0;\ + ha_pt1_t *g;\ + if (n == 0) return 0;\ + g = &h->h[a[0].x&mask];\ + for (j = 0; j < n; ++j) {\ + uint64_t x = a[j].x >> h->pre;\ + khint_t k;\ + int n;\ + IType *p;\ + assert((a[j].x&mask) == (a[0].x&mask));\ + k = yak_pt_get(g->h, x<h)) continue; \ + n = kh_key(g->h, k) & YAK_MAX_COUNT;\ + assert(n < YAK_MAX_COUNT);\ + p = &g->Ia[kh_val(g->h, k) + n];\ + p->rid = a[j].rid, p->rev = a[j].rev, p->pos = a[j].pos, p->span = a[j].span;\ + /**(uint64_t)a[j].rid<<36 | (uint64_t)a[j].rev<<35 | (uint64_t)a[j].pos<<8 | (uint64_t)a[j].span;**/\ + ++kh_key(g->h, k);\ + ++n_ins;\ + }\ + return n_ins;\ +}\ +/** data structure for each step in kt_pipeline()**/\ +typedef struct {pl_data_t *p;uint64_t n_seq0; int n_seq, m_seq, sum_len, nk, uq, *len; char **seq; VType *mz_buf; VType *mz;sf##_ch_buf_t *buf;st_mt_t *mt;} sf##_st_data_t;\ +static void sf##_worker_for_insert(void *data, long i, int tid) /** callback for kt_for()**/\ +{\ + sf##_st_data_t *s = (sf##_st_data_t*)data;\ + sf##_ch_buf_t *b = &s->buf[i];\ + if (s->p->pt){\ + if(s->p->flag&HAF_COUNT_REFINE) b->n_ins += ha_pt_cnt_insert_list(s->p->pt, b->n, b->a);\ + else b->n_ins += sf##_ha_pt_insert_list(s->p->pt, b->n, b->b);\ + }else{\ + b->n_ins += ha_ct_insert_list(s->p->ct, s->p->create_new, b->n, b->a);\ + }\ +}\ +static void sf##_worker_for_mz(void *data, long i, int tid)\ +{\ + sf##_st_data_t *s = (sf##_st_data_t*)data;\ + /**get the corresponding minimzer vector of this read**/\ + VType *b = &s->mz_buf[tid];\ + s->mz_buf[tid].n = 0;\ + sf##_ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, asm_opt.mz_sample_dist, 0, 0, \ + (s->p->pt&&(s->p->flag&HAF_COUNT_REFINE))?s->p->pt:NULL, s->p->opt->min_rcnt, asm_opt.dp_min_len, asm_opt.dp_e, &(s->mt[tid]), asm_opt.mz_rewin, s->uq);\ + s->mz[i].n = s->mz[i].m = b->n;\ + MALLOC(s->mz[i].a, b->n);\ + memcpy(s->mz[i].a, b->a, b->n * sizeof(VType));\ +}\ +static inline void sf##_pt_insert_buf(sf##_ch_buf_t *buf, int p, const HType *y){\ + /**assign minimizer to one of 4096 bins by low 12 bits**/\ + int pre = y->x & ((1<n == b->m) {\ + b->m = b->m < 8? 8 : b->m + (b->m>>1);\ + REALLOC(b->b, b->m);\ + }\ + b->b[b->n++] = *y;\ +}\ +static void *sf##_worker_count(void *data, int step, void *in) /** callback for kt_pipeline()**/\ +{\ + pl_data_t *p = (pl_data_t*)data;\ + if (step == 0) { /** step 1: read a block of sequences**/\ + int ret;\ + sf##_st_data_t *s;\ + CALLOC(s, 1);\ + s->p = p;\ + s->n_seq0 = p->n_seq;\ + if (p->rs_in && (p->flag & HAF_RS_READ)) {\ + while (p->n_seq < p->rs_in->total_reads) {\ + if ((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU) {\ + ++p->n_seq;\ + continue;\ + }\ + int l;\ + recover_UC_Read(&p->ucr, p->rs_in, p->n_seq);\ + l = p->ucr.length;\ + if (s->n_seq == s->m_seq) {\ + s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);\ + REALLOC(s->len, s->m_seq);\ + REALLOC(s->seq, s->m_seq);\ + }\ + MALLOC(s->seq[s->n_seq], l);\ + memcpy(s->seq[s->n_seq], p->ucr.seq, l);\ + s->len[s->n_seq++] = l;\ + ++p->n_seq;\ + s->sum_len += l;\ + s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;\ + if (s->sum_len >= p->opt->chunk_size)\ + break;\ + }\ + } else if(p->us_in) {\ + ma_utg_t *u; s->uq = p->us_in->h;\ + while (p->n_seq < p->us_in->n) {\ + u = &(p->us_in->a[p->n_seq]);\ + if (s->n_seq == s->m_seq) {\ + s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);\ + REALLOC(s->len, s->m_seq);\ + REALLOC(s->seq, s->m_seq);\ + }\ + MALLOC(s->seq[s->n_seq], u->len);\ + memcpy(s->seq[s->n_seq], u->s, u->len);\ + s->len[s->n_seq++] = u->len;\ + ++p->n_seq;\ + s->sum_len += u->len;\ + s->nk += u->len >= p->opt->k? u->len - p->opt->k + 1 : 0;\ + if (s->sum_len >= p->opt->chunk_size)\ + break;\ + }\ + } else {\ + while ((ret = kseq_read(p->ks)) >= 0) {\ + int l = (int)(p->ks->seq.l) - (int)(p->opt->adaLen) - (int)(p->opt->adaLen);\ + if(l <= 0) continue;\ + if (p->n_seq >= 1<<28) {\ + fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);\ + exit(1);\ + }\ + if (p->rs_out) {\ + /**for 0-th count, just insert read length to R_INF, instead of read**/\ + if (p->flag & HAF_RS_WRITE_LEN) {\ + assert(p->n_seq == p->rs_out->total_reads);\ + ha_insert_read_len(p->rs_out, l, p->ks->name.l);\ + } else if (p->flag & HAF_RS_WRITE_SEQ) {\ + int i, n_N;\ + assert(l == (int)p->rs_out->read_length[p->n_seq]);\ + for (i = n_N = 0; i < l; ++i) /** count number of ambiguous bases**/\ + if (seq_nt4_table[(uint8_t)p->ks->seq.s[i+p->opt->adaLen]] >= 4)\ + ++n_N;\ + ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s+p->opt->adaLen, l, &p->rs_out->N_site[p->n_seq], n_N);\ + memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);\ + }\ + }\ + if (s->n_seq == s->m_seq) {\ + s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);\ + REALLOC(s->len, s->m_seq);\ + REALLOC(s->seq, s->m_seq);\ + }\ + MALLOC(s->seq[s->n_seq], l);\ + memcpy(s->seq[s->n_seq], p->ks->seq.s+p->opt->adaLen, l);\ + s->len[s->n_seq++] = l;\ + ++p->n_seq;\ + s->sum_len += l;\ + s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;\ + /**p->opt->chunk_size is the block max size**/\ + if (s->sum_len >= p->opt->chunk_size)\ + break;\ + }\ + }\ + if (s->sum_len == 0) free(s);\ + else return s;\ + } else if (step == 1) { /** step 2: extract k-mers**/\ + /**s is the block of reads**/\ + sf##_st_data_t *s = (sf##_st_data_t*)in;\ + int i, n_pre = 1<opt->pre, m;\ + /**allocate the k-mer buffer**/\ + CALLOC(s->buf, n_pre);\ + m = (int)(s->nk * 1.2 / n_pre) + 1;\ + /**pre-allocate memory for each of 4096 buffer**/\ + for (i = 0; i < n_pre; ++i) {\ + s->buf[i].m = m;\ + /**for 0-th counting, p->pt = NULL**/\ + if (p->pt && !(p->flag&HAF_COUNT_REFINE)) MALLOC(s->buf[i].b, m);\ + else MALLOC(s->buf[i].a, m);\ + }\ + if (p->opt->w == 1) { /** enumerate all k-mers**/\ + int i;\ + for (i = 0; i < s->n_seq; ++i) {\ + if (p->opt->is_HPC)\ + sf##_count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);\ + else\ + sf##_count_seq_buf(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);\ + if (!p->is_store) free(s->seq[i]);\ + }\ + } else { /** minimizers only**/\ + uint32_t j;\ + /**s->n_seq is how many reads at this buffer**/\ + /**s->mz && s->mz_buf are lists of minimzer vectors**/\ + CALLOC(s->mz, s->n_seq), CALLOC(s->mz_buf, p->opt->n_thread), CALLOC(s->mt, p->opt->n_thread);\ + /**calculate minimzers for each read, each read corresponds to one thread**/\ + kt_for(p->opt->n_thread, sf##_worker_for_mz, s, s->n_seq);\ + for (i = 0; i < p->opt->n_thread; ++i) free(s->mt[i].a), free(s->mz_buf[i].a);\ + free(s->mt), free(s->mz_buf);\ + /**insert minimizers**/\ + if (p->pt && !(p->flag&HAF_COUNT_REFINE)) {/**insert whole minimizer**/\ + for (i = 0; i < s->n_seq; ++i)\ + for (j = 0; j < s->mz[i].n; ++j)\ + sf##_pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]);\ + } else {/**just insert the hash key of minimizer**/\ + for (i = 0; i < s->n_seq; ++i)\ + for (j = 0; j < s->mz[i].n; ++j)\ + sf##_ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);\ + }\ + for (i = 0; i < s->n_seq; ++i) {\ + p->n_mz += s->mz[i].n;\ + free(s->mz[i].a);\ + if (!p->is_store) free(s->seq[i]);\ + }\ + free(s->mz);\ + }\ + /**just clean seq**/\ + free(s->seq); free(s->len);\ + s->seq = 0, s->len = 0;\ + return s;\ + } else if (step == 2) { /** step 3: insert k-mers to hash table**/\ + sf##_st_data_t *s = (sf##_st_data_t*)in;\ + int i, n = 1<opt->pre;uint64_t n_ins = 0;\ + /**for 0-th counting, p->pt = NULL**/\ + kt_for(p->opt->n_thread, sf##_worker_for_insert, s, n);\ + /**n_ins is number of distinct k-mers**/\ + for (i = 0; i < n; ++i) {\ + n_ins += s->buf[i].n_ins;\ + if (p->pt && !(p->flag&HAF_COUNT_REFINE)) free(s->buf[i].b);\ + else free(s->buf[i].a);\ + }\ + if (p->ct) p->ct->tot += n_ins, p->ct->bs += s->sum_len;\ + if (p->pt) p->pt->tot_pos += n_ins;\ + free(s->buf);\ + free(s);\ + }\ + return 0;\ } -static void worker_for_mz(void *data, long i, int tid) -{ - st_data_t *s = (st_data_t*)data; - ///get the corresponding minimzer vector of this read - ha_mz1_v *b = &s->mz_buf[tid]; - s->mz_buf[tid].n = 0; - ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, asm_opt.mz_sample_dist, 0, 0, - (s->p->pt&&(s->p->flag&HAF_COUNT_REFINE))?s->p->pt:NULL, s->p->opt->min_rcnt, asm_opt.dp_min_len, asm_opt.dp_e, &(s->mt[tid]), asm_opt.mz_rewin, s->uq); - s->mz[i].n = s->mz[i].m = b->n; - MALLOC(s->mz[i].a, b->n); - memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t)); -} +MZ_TEST_INIT(mz1, ha_mz1_t, ha_mz1_v, ha_idxpos_t, a) +MZ_TEST_INIT(mz2, ha_mzl_t, ha_mzl_v, ha_idxposl_t, al) -MZ_FUNC_INIT(mz, ha_mz1_t) -MZ_FUNC_INIT(mzl, ha_mzl_t) -static void *worker_count(void *data, int step, void *in) // callback for kt_pipeline() -{ - pl_data_t *p = (pl_data_t*)data; - if (step == 0) { // step 1: read a block of sequences - int ret; - st_data_t *s; - CALLOC(s, 1); - s->p = p; - s->n_seq0 = p->n_seq; - if (p->rs_in && (p->flag & HAF_RS_READ)) { - while (p->n_seq < p->rs_in->total_reads) { - if ((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU) { - ++p->n_seq; - continue; - } - int l; - recover_UC_Read(&p->ucr, p->rs_in, p->n_seq); - l = p->ucr.length; - if (s->n_seq == s->m_seq) { - s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1); - REALLOC(s->len, s->m_seq); - REALLOC(s->seq, s->m_seq); - } - MALLOC(s->seq[s->n_seq], l); - memcpy(s->seq[s->n_seq], p->ucr.seq, l); - s->len[s->n_seq++] = l; - ++p->n_seq; - s->sum_len += l; - s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0; - if (s->sum_len >= p->opt->chunk_size) - break; - } - } else if(p->us_in) { - ma_utg_t *u; s->uq = 1; - while (p->n_seq < p->us_in->n) { - u = &(p->us_in->a[p->n_seq]); - if (s->n_seq == s->m_seq) { - s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1); - REALLOC(s->len, s->m_seq); - REALLOC(s->seq, s->m_seq); - } - MALLOC(s->seq[s->n_seq], u->len); - memcpy(s->seq[s->n_seq], u->s, u->len); - s->len[s->n_seq++] = u->len; - ++p->n_seq; - s->sum_len += u->len; - s->nk += u->len >= p->opt->k? u->len - p->opt->k + 1 : 0; - if (s->sum_len >= p->opt->chunk_size) - break; - } - } else { - while ((ret = kseq_read(p->ks)) >= 0) { - int l = (int)(p->ks->seq.l) - (int)(p->opt->adaLen) - (int)(p->opt->adaLen); - if(l <= 0) continue; - - if (p->n_seq >= 1<<28) { - fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28); - exit(1); - } - if (p->rs_out) { - ///for 0-th count, just insert read length to R_INF, instead of read - if (p->flag & HAF_RS_WRITE_LEN) { - assert(p->n_seq == p->rs_out->total_reads); - ha_insert_read_len(p->rs_out, l, p->ks->name.l); - } else if (p->flag & HAF_RS_WRITE_SEQ) { - int i, n_N; - assert(l == (int)p->rs_out->read_length[p->n_seq]); - for (i = n_N = 0; i < l; ++i) // count number of ambiguous bases - if (seq_nt4_table[(uint8_t)p->ks->seq.s[i+p->opt->adaLen]] >= 4) - ++n_N; - ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s+p->opt->adaLen, l, &p->rs_out->N_site[p->n_seq], n_N); - memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l); - } - } - ///for 0-th count, insert both seq and length to local block - if (s->n_seq == s->m_seq) { - s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1); - REALLOC(s->len, s->m_seq); - REALLOC(s->seq, s->m_seq); - } - MALLOC(s->seq[s->n_seq], l); - memcpy(s->seq[s->n_seq], p->ks->seq.s+p->opt->adaLen, l); - s->len[s->n_seq++] = l; - ++p->n_seq; - s->sum_len += l; - s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0; - ///p->opt->chunk_size is the block max size - if (s->sum_len >= p->opt->chunk_size) - break; - } - } - if (s->sum_len == 0) free(s); - else return s; - } else if (step == 1) { // step 2: extract k-mers - ///s is the block of reads - st_data_t *s = (st_data_t*)in; - if(p->us_in) mzl_init_kuf(p, s); - else mz_init_kuf(p, s); - // fill the buffer - ///for 0-th counting, p->opt->w == 1 - if (p->opt->w == 1) { // enumerate all k-mers - ///scan all reads - int i; - for (i = 0; i < s->n_seq; ++i) { - if (p->opt->is_HPC) - count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]); - else - count_seq_buf(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]); - if (!p->is_store) free(s->seq[i]); - } - } else { // minimizers only - if(p->us_in) mzl_mselect(p, s); - else mz_mselect(p, s); - } - ///just clean seq - free(s->seq); free(s->len); - s->seq = 0, s->len = 0; - return s; - } else if (step == 2) { // step 3: insert k-mers to hash table - st_data_t *s = (st_data_t*)in; - ///for 0-th counting, p->pt = NULL - kt_for(p->opt->n_thread, worker_for_insert, s, 1<opt->pre); - if(p->us_in) mzl_destory_kuf(p, s, 1<opt->pre); - else mz_destory_kuf(p, s, 1<opt->pre); - } - return 0; -} void debug_adapter(const hifiasm_opt_t *asm_opt, All_reads *rs) { @@ -951,7 +919,8 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt pl.ct = ha_ct_init(opt->k, opt->pre, opt->bf_n_hash, opt->bf_shift); } if(pl.ct) pl.ct->bs = 0; - kt_pipeline(3, worker_count, &pl, 3); + if(ug_rs) kt_pipeline(3, mz2_worker_count, &pl, 3); + else kt_pipeline(3, mz1_worker_count, &pl, 3); if (read_rs) { destory_UC_Read(&pl.ucr); } else if(!read_rs && !ug_rs) { @@ -963,7 +932,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt return pl.ct; } -ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const void *flt_tab, All_reads *rs, ma_utg_v *us, int keep_adapter, int *low_freq) +ha_ct_t *ha_count(const hifiasm_opt_t *asm_o, int flag, int HPC, int k, int w, ha_pt_t *p0, const void *flt_tab, All_reads *rs, ma_utg_v *us, int keep_adapter, int *low_freq) { int i; int64_t n_seq = 0; @@ -979,20 +948,19 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi malloc_All_reads(rs); } yak_copt_init(&opt); - opt.k = us? asm_opt->ul_mer_length:asm_opt->k_mer_length; - ///always 0 - opt.is_HPC = !(asm_opt->flag&HA_F_NO_HPC); + opt.k = k; + opt.is_HPC = HPC; ///for ft-counting, shoud be 1 - opt.w = flag & HAF_COUNT_ALL? 1 : (us? asm_opt->ul_mz_win:asm_opt->mz_win); + opt.w = flag & HAF_COUNT_ALL? 1 : w; ///for ft-counting, shoud be 37 ///for ha_pt_gen, shoud be 0 - opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift; - opt.n_thread = asm_opt->thread_num; - opt.adaLen = (keep_adapter? asm_opt->adapterLen : 0); + opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_o->bf_shift; + opt.n_thread = asm_o->thread_num; + opt.adaLen = (keep_adapter? asm_o->adapterLen : 0); opt.min_rcnt = (low_freq?*low_freq:-1); ///asm_opt->num_reads is the number of fastq files - for (i = n_bs = 0; i < (us?1:asm_opt->num_reads); ++i){ - h = yak_count(&opt, asm_opt->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq); + for (i = n_bs = 0; i < (us?1:asm_o->num_reads); ++i){ + h = yak_count(&opt, asm_o->read_file_names[i], flag|HAF_CREATE_NEW, p0, h, flt_tab, rs, us, &n_seq); if(h) n_bs += h->bs; } if(h) h->bs = n_bs; @@ -1077,32 +1045,32 @@ void debug_ct_index(void* q_ct_idx, void* r_ct_idx) * High-level interfaces * *************************/ -void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int hap_n) +void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq) { yak_ft_t *flt_tab; ha_ct_t *h; ///HAF_COUNT_EXACT ---> no bf; HAF_COUNT_ALL ---> no minimizer - h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, NULL, NULL, NULL, us, 0, NULL); - ha_ct_shrink(h, 1, YAK_MAX_COUNT-1, asm_opt->thread_num); + h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ|HAF_COUNT_EXACT, is_HPC, k, w, NULL, NULL, NULL, us, 0, NULL); + ha_ct_shrink(h, min_freq, max_freq>YAK_MAX_COUNT-1?YAK_MAX_COUNT-1:max_freq, asm_opt->thread_num); flt_tab = gen_hh(h, asm_opt->max_kmer_cnt); ha_ct_destroy(h); return (void*)flt_tab; } -ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int hap_n) +ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int is_HPC, int k, int w, int min_freq) { ha_ct_t *ct; ha_pt_t *pt; ///HAF_COUNT_EXACT: no bf - ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, NULL, flt_tab, NULL, us, 0, NULL); + ct = ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, NULL, flt_tab, NULL, us, 0, NULL); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__, yak_realtime(), yak_cpu_usage(), (long)ct->tot); ///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen - ha_ct_shrink(ct, 1, YAK_MAX_COUNT - 1, asm_opt->thread_num); + ha_ct_shrink(ct, min_freq, YAK_MAX_COUNT - 1, asm_opt->thread_num); pt = ha_pt_gen(ct, asm_opt->thread_num, 1); - ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, pt, flt_tab, NULL, us, 0, NULL); + ha_count(asm_opt, HAF_COUNT_EXACT|HAF_UG_READ, is_HPC, k, w, pt, flt_tab, NULL, us, 0, NULL); //ha_pt_sort(pt, asm_opt->thread_num); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__, yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos); @@ -1116,7 +1084,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i int peak_hom, peak_het, cutoff = YAK_MAX_COUNT - 1, ex_flag = 0; if(is_hp_mode) ex_flag = HAF_RS_READ|HAF_SKIP_READ; ha_ct_t *h; - h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, NULL, NULL, rs, NULL, 1, NULL); + h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN|ex_flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, NULL, rs, NULL, 1, NULL); if((asm_opt->flag & HA_F_VERBOSE_GFA)) { write_ct_index((void*)h, asm_opt->output_file_name); @@ -1148,12 +1116,12 @@ ha_pt_t *ha_pt_gen_dp(const hifiasm_opt_t *asm_opt, ha_ct_t *ct, int flag, int n { int low_freq = mz_low_b(peak_hom, peak_het); ha_pt_t *pt = ha_pt_gen_count(ct, n_thread); ///key = cnt, val = 0 - ha_count(asm_opt, HAF_COUNT_EXACT|HAF_COUNT_REFINE|flag, pt, flt_tab, rs, NULL, 1, &low_freq); + ha_count(asm_opt, HAF_COUNT_EXACT|HAF_COUNT_REFINE|flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, &low_freq); uint64_t occ = ha_pt_shrink(pt, n_thread); if(flag&HAF_RS_WRITE_LEN) flag -= HAF_RS_WRITE_LEN; if(flag&HAF_RS_WRITE_SEQ) flag -= HAF_RS_WRITE_SEQ; flag |= HAF_RS_READ; pt->tot_pos = 0; - ha_count(asm_opt, HAF_COUNT_EXACT|flag, pt, flt_tab, rs, NULL, 1, NULL); + ha_count(asm_opt, HAF_COUNT_EXACT|flag, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, NULL); // fprintf(stderr, "[M::%s::] counted %lu distinct minimizer k-mers\n", __func__, pt->tot); // fprintf(stderr, "[M::%s::] collected %lu minimizers\n\n\n", __func__, pt->tot_pos); assert(occ == pt->tot_pos); @@ -1177,7 +1145,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f } if(is_hp_mode) extra_flag1 |= HAF_SKIP_READ, extra_flag2 |= HAF_SKIP_READ; - ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs, NULL, 1, NULL); + ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, NULL, flt_tab, rs, NULL, 1, NULL); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__, yak_realtime(), yak_cpu_usage(), (long)ct->tot); ha_ct_hist(ct, cnt, asm_opt->thread_num); @@ -1203,7 +1171,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f { fprintf(stderr, "[M::%s::] counting in normal mode\n", __func__); pt = ha_pt_gen(ct, asm_opt->thread_num, 0); - ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, pt, flt_tab, rs, NULL, 1, NULL); + ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, !(asm_opt->flag&HA_F_NO_HPC), asm_opt->k_mer_length, asm_opt->mz_win, pt, flt_tab, rs, NULL, 1, NULL); assert((uint64_t)tot_cnt == pt->tot_pos); } else diff --git a/htab.h b/htab.h index 28be901..b2f9351 100644 --- a/htab.h +++ b/htab.h @@ -25,14 +25,12 @@ typedef struct { uint32_t n, m; ha_mz1_t *a; } ha_mz1_v; typedef struct { uint64_t x; ///x is the hash key - uint64_t rid:31, rev:1; - uint32_t pos; + uint64_t rid:31, rev:1, pos:32; uint8_t span; } ha_mzl_t; typedef struct { - uint64_t rid:31, rev:1; - uint32_t pos; + uint64_t rid:31, rev:1, pos:32; uint8_t span; } ha_idxposl_t; @@ -71,12 +69,12 @@ extern void *ha_flt_tab_hp; extern ha_pt_t *ha_idx_hp; extern void *ha_ct_table; -void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int hap_n); +void *ha_ft_ug_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq); void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode); int32_t ha_ft_cnt(const void *hh, uint64_t y); void ha_ft_destroy(void *h); -ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int hap_n); +ha_pt_t *ha_pt_ug_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_v *us, int is_HPC, int k, int w, int min_freq); ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov); void ha_pt_destroy(ha_pt_t *h); const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n); @@ -101,7 +99,8 @@ double yak_cpu_usage(void); void ha_triobin(const hifiasm_opt_t *opt); -void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique); +void mz1_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique); +void mz2_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mzl_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique); int ha_analyze_count(int n_cnt, int start_cnt, int m_peak_hom, const int64_t *cnt, int *peak_het); int adj_m_peak_hom(int m_peak_hom, int max_i, int max2_i, int max3_i, int *peak_het); void print_hist_lines(int n_cnt, int start_cnt, const int64_t *cnt); diff --git a/inter.cpp b/inter.cpp index b3c685a..3b5ef63 100644 --- a/inter.cpp +++ b/inter.cpp @@ -7,19 +7,26 @@ #include "CommandLines.h" #include "htab.h" -void uidx_build(ma_ug_t *ug, int hap_n) +void uidx_build(ma_ug_t *ug, int is_HPC, int k, int w, int hap_n) { int flag = asm_opt.flag; asm_opt.flag |= HA_F_NO_HPC; - ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), hap_n); - ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), hap_n); - - - - - - + ug->u.h = hap_n; + ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), is_HPC, k, w, hap_n, hap_n*10); + ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), is_HPC, k, w, hap_n); asm_opt.flag = flag; +} + +void uidx_destory() +{ ha_ft_destroy(ha_flt_tab); ha_pt_destroy(ha_idx); } + + +void ul_resolve(ma_ug_t *ug, int hap_n) +{ + uidx_build(ug, 1, 63, 63, hap_n); + + uidx_destory(); +} \ No newline at end of file diff --git a/rcut.cpp b/rcut.cpp index 22932e5..5a8019b 100644 --- a/rcut.cpp +++ b/rcut.cpp @@ -2767,6 +2767,48 @@ void clean_ovlp_by_mc(mc_g_t *mg, hap_overlaps_list* ha) } +void filter_ta_by_mc(ma_ug_t *ug, asg_t *read_g, uint32_t uID, kv_u_trans_t* ta, mc_match_t *ma, int8_t *s) +{ + mc_edge_t *o = pt_a(*ma, uID); + uint32_t n = pt_n(*ma, uID), k, qn, tn; + u_trans_t *p = NULL; + for (k = 0; k < n; ++k) + { + qn = ma_x(o[k]); tn = ma_y(o[k]); p = NULL; + if((s[qn]*s[tn])!=-1) continue; + get_u_trans_spec(ta, qn, tn, &p, NULL); + if(p && p->nw == o[k].w) { + p->del = 0; + } + else { + get_u_trans_spec(ta, tn, qn, &p, NULL); + if(p && p->nw == o[k].w) p->del = 0; + } + if(!p) fprintf(stderr, "ERROR-ta-p\n"); + } +} + + +void clean_ta_by_mc(mc_g_t *mg, kv_u_trans_t *ta) +{ + uint32_t v, i; + u_trans_t *p = NULL; + for (i = 0; i < ta->n; i++) ta->a[i].del = 1; + for (i = 0; i < mg->e->n_seq; ++i) filter_ta_by_mc(mg->ug, mg->rg, i, ta, mg->e, mg->s.a); + for (i = v = 0; i < ta->n; i++) { + if(ta->a[i].del) continue; + ta->a[v++] = ta->a[i]; + } + ta->n = v; + for (i = 0; i < v; i++) { + kv_pushp(u_trans_t, *ta, &p); + (*p) = ta->a[i]; + p->qn = ta->a[i].tn; p->qs = ta->a[i].ts; p->qe = ta->a[i].te; + p->tn = ta->a[i].qn; p->ts = ta->a[i].qs; p->te = ta->a[i].qe; + } + kt_u_trans_t_idx(ta, mg->ug->g->n_seq); +} + void p_nodes(mc_g_t *mg, trans_chain* t_ch, uint8_t* trio_flag) { uint32_t i; @@ -2844,7 +2886,7 @@ void print_hap_s(int8_t *s, uint32_t sn) } } -void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref) +void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov) { mc_opt_t opt; mc_opt_init(&opt, asm_opt.n_perturb, asm_opt.f_perturb, asm_opt.seed); @@ -2863,8 +2905,10 @@ void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_u p_nodes(mg, t_ch, trio_flag); } - if(ovlp) clean_ovlp_by_mc(mg, ovlp); - + if(clean_ov){ + if(ovlp) clean_ovlp_by_mc(mg, ovlp); + if(ta) clean_ta_by_mc(mg, ta); + } // print_hap_s(s, ug->u.n); destory_mc_g_t(&mg); diff --git a/rcut.h b/rcut.h index 606d64d..6510cb8 100644 --- a/rcut.h +++ b/rcut.h @@ -118,7 +118,7 @@ static inline double kr_drand_r(uint64_t *x) return u.d - 1.0; } -void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref); +void mc_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag, uint32_t renew_s, int8_t *s, uint32_t is_sys, bubble_type* bub, kv_u_trans_t *ref, int clean_ov); void debug_mc_g_t(const char* name); void mc_solve_general(kv_u_trans_t *ta, uint32_t un, kv_gg_status *s, uint16_t hapN, uint16_t update_ta, uint16_t write_dump); kv_gg_status *init_mc_gg_status(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, diff --git a/sketch.cpp b/sketch.cpp index a7b08b7..1039eea 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -10,14 +10,8 @@ #define MAX_HIGH_OCC 8 // TODO: don't hard code if we need to tune this parameter #define MAX_MAX_HIGH_OCC 16 #define GMC(a, x,y,xn) ((a)[(x)*(xn)+(y)]) - -static inline int mzcmp(const ha_mz1_t *a, const ha_mz1_t *b) -{ - return a->rid < b->rid? -1 : a->rid > b->rid? 1 : ((a->x > b->x) - (a->x < b->x)); -} - -#define mz_lt(a, b) (mzcmp(&(a), &(b)) < 0) -KSORT_INIT(mz, ha_mz1_t, mz_lt) +#define GL(x, i) ((int64_t)((uint32_t)((x).a[(i)]))) +#define A_M(p, i) ((i) >= 0 && (p).a[(i)].rid > 0) void debug_refine(ha_mz1_t *ma, uint64_t *mmt, int32_t sn, int32_t n, int32_t m, int32_t end) { @@ -42,287 +36,6 @@ void debug_refine(ha_mz1_t *ma, uint64_t *mmt, int32_t sn, int32_t n, int32_t m, if(nt != tot) fprintf(stderr, "ERROR-TOT, nt: %ld, tot: %ld\n", nt, tot); } -void refine_select(ha_mz1_v *mz, int32_t sidx, int32_t eidx, int32_t sn, int32_t min_freq, st_mt_t *mm, -int32_t *rsi, int32_t *rei) -{ - int32_t n = sn, m = eidx + 1 - sidx, i, k, t, mk=-1; - uint64_t ix, kx, ks; - kv_resize(uint64_t, *mm, mm->n+n*m); - ha_mz1_t *ma = mz->a + sidx; - uint64_t *mmt = mm->a + mm->n; - // fprintf(stderr, "[M::%s::] ==> +n: %d, m: %d, sn: %d, sidx: %d, eidx: %d\n", __func__, n, m, sn, sidx, eidx); - - for (i = 0; i < n; i++) ///how many selected minimizers - { - for (k = 0, mk = -1; k < m; k++) ///how many minimizers in total - { - if((int32_t)(ma[k].rid) 0) - { - for (t = k-1; t >= 0 && (ma[t].pos >= ks||(int32_t)(ma[t].rid)0?(i-1)*m+t:0xffffffff)<<32; - else if(ks == kx) ks |= (uint64_t)(mk>=0?i*m+mk:0xffffffff)<<32; - - GMC(mmt, i,k,m) = ks; - mk = k; - } - } - // fprintf(stderr, "[M::%s::] ==> ++n: %d, m: %d, sn: %d, sidx: %d, eidx: %d\n", __func__, n, m, sn, sidx, eidx); - - ks = (n-1)*m + mk; ix = (uint64_t)-1; kx = 0; - while (ks != 0xffffffff) - { - i = ks/m; k = ks%m; - ks = mmt[ks]>>32; - // fprintf(stderr, "i: %d, k: %d, ks: %lu\n", i, k, ks); - if(ks == 0xffffffff || (int32_t)(ks/m) == (i-1)) - { - mm->a[sidx+k] = 1; - ix = MIN((uint64_t)k, ix); kx = MAX((uint64_t)k, kx); - } - } - ///debug - // debug_refine(ma, mmt, sn, n, m, (n-1)*m + mk); - - if(rsi) (*rsi) = ix + sidx; - if(rei) (*rei) = kx + sidx; -} - -void refine_sketch(ha_mz1_v *p, ha_pt_t *pt, int32_t rlen, int32_t dp_min_len, float er, int32_t min_freq, st_mt_t *mt) -{ - // fprintf(stderr, "[M::%s::] ==> #########10#########, rlen: %d\n", __func__, rlen); - - int32_t i, n = p->n, bd, len = MIN(rlen, dp_min_len), sublen, cnt, ei, li, ri; - int32_t sn = len*er + 1; - kv_resize(uint64_t, *mt, (int64_t)p->n); - mt->n = p->n; memset(mt->a, 0, sizeof(uint64_t)*p->n); - for (i = 0; i < n; i++) p->a[i].rid = ha_pt_cnt(pt, p->a[i].x); - - for (i = cnt = 0, bd = -1, ei = -1; i < n; i++) - { - if((int32_t)(p->a[i].rid)a[i].pos + 1; - if(sublen > len) break; - else ei = i; - - if((int32_t)(p->a[i].pos + 1 - p->a[i].span) > bd) - { - bd = p->a[i].pos; - cnt++; - } - } - - // fprintf(stderr, "[M::%s::] ==> +cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n); - - - if(cnt >= sn) refine_select(p, 0, ei, sn, min_freq, mt, NULL, &li); - else - { - li = i-1; - for (i = 0; i <= li; i++) mt->a[i] = 1; - } - - - if(len < rlen) - { - for (i = n-1, cnt = 0, bd = rlen+1, ei = -1; i >= 0; i--) - { - if((int32_t)(p->a[i].rid)a[i].pos + 1 - p->a[i].span); - if(sublen > len) break; - else ei = i; - - if((int32_t)(p->a[i].pos) < bd) - { - bd = p->a[i].pos + 1 - p->a[i].span; - cnt++; - } - } - - // fprintf(stderr, "[M::%s::] ==> -cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n); - - - if(cnt >= sn) refine_select(p, ei, n-1, sn, min_freq, mt, &ri, NULL); - else - { - ri = i+1; - for (i = ri; i <= n-1; i++) mt->a[i] = 1; - } - - // fprintf(stderr, "[M::%s::] ==> --cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n); - - if(ri - li >= 2) - { - li++; ri--; - sn = (p->a[ri].pos - p->a[li].pos + p->a[li].span)*er + 1; - for (i = li, cnt = 0, bd = -1; i <= ri; i++) - { - if((int32_t)(p->a[i].rid)a[i].pos + 1 - p->a[i].span) > bd) - { - bd = p->a[i].pos; - cnt++; - if(cnt >= sn) break; - } - } - - if(cnt >= sn) refine_select(p, li, ri, sn, min_freq, mt, NULL, NULL); - else for (i = li; i <= ri; i++) mt->a[i] = 1; - } - } - - // fprintf(stderr, "[M::%s::] ==> #########20#########, p->n: %u, n: %d\n", __func__, p->n, n); - for (i = sn = 0; i < n; i++) - { - if(mt->a[i]) - { - p->a[sn] = p->a[i]; - sn++; - } - } - // if(p->n != sn) fprintf(stderr, "[M::%s::] ==> #########21#########, p->n: %u, sn: %d\n", __func__, p->n, sn); - p->n = sn; - -} - -inline int hf_dp(ha_mz1_v *mz, int32_t sidx, int32_t eidx, int32_t sn, int32_t min_freq, st_mt_t *mm, -int32_t *rsi, int32_t *rei) -{ - return 0; -} - -inline void hf_select(ha_mz1_v *p, int32_t si, int32_t ei, int32_t n, int32_t len, int32_t sample_dist, ha_mz1_t *b, int32_t force) -{ - if(ei - si <= 1) return; - int32_t ps = si < 0? 0 : p->a[si].pos; - int32_t pe = ei == n? len : p->a[ei].pos; - int32_t j, k, st = si + 1, en = ei; - int32_t max_high_occ = (int32_t)((double)(pe - ps) / sample_dist + .499); - if (max_high_occ > MAX_MAX_HIGH_OCC) - max_high_occ = MAX_MAX_HIGH_OCC; - for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k) - b[k] = p->a[j], b[k].pos = j; // b[].pos keeps the index in p->a[] - ks_heapmake_mz(k, b); // initialize the binomial heap - for (; j < en; ++j) { // if there are more, choose top max_high_occ - if (mz_lt(p->a[j], b[0])) { // then update the heap - b[0] = p->a[j], b[0].pos = j; - ks_heapdown_mz(0, k, b); - } - } - //ks_heapsort_mz(k, b); // sorting is not needed for now - for (j = 0; j < k; ++j) - if (b[j].rid < pe - ps || force) - p->a[b[j].pos].rid = 0; -} - -void select_mz(ha_mz1_v *p, int len, int sample_dist, int32_t dp_min_len) -{ // for high-occ minimizers, choose up to max_high_occ in each high-occ streak - int32_t i, last0 = -1, n = (int32_t)p->n, m = 0, nw[2], min_len; - ha_mz1_t b[MAX_MAX_HIGH_OCC]; // this is to avoid a heap allocation - - if (n == 0 || n == 1) return; - assert(n < 1<<27); // 27 is the number of bits for ha_mz1_t::pos; this should be safe as there are more bases than minimizers - for (i = 0; i < n; ++i) - if (p->a[i].rid != 0) ++m; - if (m == 0) return; // no high-frequency k-mers; do nothing - for (i = 0; i <= n; ++i) { - if (i == n || p->a[i].rid == 0) { - if (i - last0 > 1) { - hf_select(p, last0, i, n, len, sample_dist, b, 0); - // int32_t ps = last0 < 0? 0 : p->a[last0].pos; - // int32_t pe = i == n? len : p->a[i].pos; - // int32_t j, k, st = last0 + 1, en = i; - // int32_t max_high_occ = (int32_t)((double)(pe - ps) / sample_dist + .499); - // if (max_high_occ > MAX_MAX_HIGH_OCC) - // max_high_occ = MAX_MAX_HIGH_OCC; - // for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k) - // b[k] = p->a[j], b[k].pos = j; // b[].pos keeps the index in p->a[] - // ks_heapmake_mz(k, b); // initialize the binomial heap - // for (; j < en; ++j) { // if there are more, choose top max_high_occ - // if (mz_lt(p->a[j], b[0])) { // then update the heap - // b[0] = p->a[j], b[0].pos = j; - // ks_heapdown_mz(0, k, b); - // } - // } - // //ks_heapsort_mz(k, b); // sorting is not needed for now - // for (j = 0; j < k; ++j) - // if (b[j].rid < pe - ps) - // p->a[b[j].pos].rid = 0; - } - last0 = i; - } - } - - min_len = MAX(dp_min_len, (p->a[0].pos+1)+sample_dist); - for (i = 0, nw[0] = nw[1] = 0; i < n; i++) - { - nw[(p->a[i].rid!=0)]++; - if((p->a[i].pos + 1) > min_len) break; - } - if(nw[0]==0 && nw[1]>0) hf_select(p, -1, i, n, len, sample_dist, b, 1); - - min_len = MAX(dp_min_len, (len - (p->a[n-1].pos + 1 - p->a[n-1].span))+sample_dist); - for (i = n-1, nw[0] = nw[1] = 0; i >= 0; i--) - { - nw[(p->a[i].rid!=0)]++; - if((len - (p->a[i].pos + 1 - p->a[i].span)) > min_len) break; - } - if(nw[0]==0 && nw[1]>0) hf_select(p, i, n, n, len, sample_dist, b, 1); - - for (i = n = 0; i < (int32_t)p->n; ++i) // squeeze out filtered minimizers - if (p->a[i].rid == 0) - p->a[n++] = p->a[i]; - p->n = n; -} - -static inline int mzcmp_l(const ha_mz1_v *p, int32_t ai, int32_t bi) -{ - if(ai >= 0 && bi >= 0){ - ha_mz1_t *a = &(p->a[ai]), *b = &(p->a[bi]); - if(a->rid > 0 && b->rid > 0) return mzcmp(a, b); - return (a->rid == 0) - (b->rid == 0); - } - return (ai < 0) - (bi < 0); -} - -#define GL(x, i) ((int64_t)((uint32_t)((x).a[(i)]))) -#define A_M(p, i) ((i) >= 0 && (p).a[(i)].rid > 0) -int32_t qfw(ha_mz1_v *p, st_mt_t *mt, int32_t n, int32_t tot_l, int32_t ws, int32_t i, int32_t *mi) -{ - int32_t m, si; - for (si = i, (*mi) = -1; i < n; i++){ - if(GL(*mt, i) >= ws || (i+1 < n && GL(*mt, i) < ws && GL(*mt, i+1) > ws) || - (i+1 == n && tot_l >= ws && GL(*mt, i) < ws)){ - for (m = si; m <= i; m++){ - if(!A_M(*p, m)) continue; - if(mzcmp_l(p, *mi, m) >= 0) (*mi) = m; - } - if((*mi) >= 0 && A_M(*p, *mi)){ - for (m = si; m <= i; m++){ - if(!A_M(*p, m)) continue; - if(mzcmp_l(p, *mi, m) == 0) mt->a[m] |= 0x100000000; - } - } - break; - } - } - return i; -} - void dbg_boundary(ha_mz1_v *p, st_mt_t *mt, int32_t w, int32_t k, int32_t tot_l) { if(tot_l < w + k -1) return; @@ -401,110 +114,6 @@ void dbg_boundary(ha_mz1_v *p, st_mt_t *mt, int32_t w, int32_t k, int32_t tot_l) } } -static void select_mz_h(ha_mz1_v *p, st_mt_t *mt, int len, int sample_dist, int32_t w, int32_t k, int32_t tot_l) -{ // for high-occ minimizers, choose up to max_high_occ in each high-occ streak - int32_t i, mi = -1, si, last0 = -1, n = (int32_t)p->n, m = 0, ws = w + k - 1; - - if (n == 0) return; - assert(n < 1<<27); // 27 is the number of bits for ha_mz1_t::pos; this should be safe as there are more bases than minimizers - - for (i = m = 0, last0 = -1; i <= n; ++i) { - if (i == n || p->a[i].rid == 0) { - if (i - last0 > 1) { - int32_t ps = last0 < 0? 0 : p->a[last0].pos; - int32_t pe = i == n? len : p->a[i].pos; - if(((int32_t)((double)(pe - ps) / sample_dist + .499)) > 0){ - last0 = -2; - m++; - break; - } - } - last0 = i; - } - } - if (m == 0) return; // no high-frequency k-mers; do nothing - if(last0 >= -1) goto ff; - i = 0; - i = qfw(p, mt, n, tot_l, ws, i, &mi); - if(i == n) goto ff; - - for (si = 0, i++; i < n; i++){ - for (; si < i; si++){ - if(GL(*mt, si) + w > GL(*mt, i)) break; - } - - // a new minimum; then write the old min - if(mzcmp_l(p, i, mi) <= 0) { - if(A_M(*p, mi)) mt->a[mi] |= 0x100000000; - mi = i; - }// old min has moved outside the window - else if(si > mi){ - if(A_M(*p, mi)) mt->a[mi] |= 0x100000000; - for (m = si, mi = -1; m <= i; m++){ - if(mzcmp_l(p, mi, m) >= 0) mi = m; - } - if(A_M(*p, mi)){ - for (m = si; m <= i; m++){ - if(!A_M(*p, m)) continue; - if(mzcmp_l(p, mi, m) == 0) mt->a[m] |= 0x100000000; - } - } - } - } - - if(A_M(*p, mi)) mt->a[mi] |= 0x100000000; - - for (i = n - 1; si < n && GL(*mt, si) + w <= tot_l + 1; si++){ - if(si > mi){ - if(A_M(*p, mi)) mt->a[mi] |= 0x100000000; - for (m = si, mi = -1; m <= i; m++){ - if(mzcmp_l(p, mi, m) >= 0) mi = m; - } - if(A_M(*p, mi)){ - for (m = si; m <= i; m++){ - if(!A_M(*p, m)) continue; - if(mzcmp_l(p, mi, m) == 0) mt->a[m] |= 0x100000000; - } - } - } - } - - /** - dbg_boundary(p, mt, w, k, tot_l); - fprintf(stderr, "\n"); - for (i = 0; i < (int32_t)p->n; ++i){ - if(p->a[i].rid == 0) continue; - fprintf(stderr, "%cl: %u, pos: %lu, cnt: %lu, key: %lu, i: %d\n", "+-"[!!(mt->a[i]&0x100000000)], - (uint32_t)mt->a[i], p->a[i].pos, p->a[i].rid, p->a[i].x, i); - // if (mt->a[i]&0x100000000){ - // fprintf(stderr, "+l: %u, pos: %lu, cnt: %lu\n", (uint32_t)mt->a[i], p->a[i].pos, p->a[i].rid); - // } - } - **/ - ha_mz1_t b[MAX_MAX_HIGH_OCC]; - for (i = 0, last0 = -1; i <= n; ++i) { - if (i == n || p->a[i].rid == 0) { - if (i - last0 > 1) { - int32_t ps = last0 < 0? 0 : p->a[last0].pos; - int32_t pe = i == n? len : p->a[i].pos; - if(((int32_t)((double)(pe - ps) / sample_dist + .499)) > 0){ - for (m = last0 + 1, mi = 0; m < i; ++m){ - if(mt->a[m]&0x100000000) p->a[m].rid = 0, mi++; - } - if(mi == 0) hf_select(p, last0, i, n, len, sample_dist, b, 0); - } - } - last0 = i; - } - } - - ff: - for (i = n = 0; i < (int32_t)p->n; ++i) // squeeze out filtered minimizers - if (p->a[i].rid == 0) - p->a[n++] = p->a[i]; - p->n = n; -} - void debug_pl(const char *str, int len, int w, int k, int is_hpc, ha_mz1_v *p, const void *hf, st_mt_t *mt) { int i, l, dbi, dbcnt = 0, kmer_span = 0; @@ -571,153 +180,402 @@ void debug_pl(const char *str, int len, int w, int k, int is_hpc, ha_mz1_v *p, c } } -/** - * Find symmetric (w,k)-minimizers on a DNA sequence - * - * @param str DNA sequence - * @param len length of $str - * @param w find a minimizer for every $w consecutive k-mers - * @param k k-mer size - * @param rid reference ID; will be copied to the output $p array - * @param is_hpc homopolymer-compressed or not - * @param p minimizers - */ -void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique) -{ ///in default, w = 51, k = 51, is_hpc = 1 - /** - uint64_t x; - uint64_t rid:28, pos:27, rev:1, span:8; - **/ - extern void *ha_ct_table; - static const ha_mz1_t dummy = { UINT64_MAX, (1<<28) - 1, 0, 0, 0}; - uint64_t shift1 = k - 1, mask = (1ULL<rid < b->rid? -1 : a->rid > b->rid? 1 : ((a->x > b->x) - (a->x < b->x));} +#define mz1_mz_lt(a, b) (mz1_mzcmp(&(a), &(b)) < 0) +KSORT_INIT(mz1_mz, ha_mz1_t, mz1_mz_lt) - assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63)); - if (dbg_ct != NULL) dbg_ct->a.n = 0; - if (k_flag != NULL) { - kv_resize(uint8_t, k_flag->a, (uint64_t)len); - k_flag->a.n = len; - memset(k_flag->a.a, 0, k_flag->a.n); - } +static inline int mz2_mzcmp(const ha_mzl_t *a, const ha_mzl_t *b){return a->rid < b->rid? -1 : a->rid > b->rid? 1 : ((a->x > b->x) - (a->x < b->x));} +#define mz2_mz_lt(a, b) (mz2_mzcmp(&(a), &(b)) < 0) +KSORT_INIT(mz2_mz, ha_mzl_t, mz2_mz_lt) - memset(buf, 0xff, w * sizeof(ha_mz1_t)); - memset(&tq, 0, sizeof(tiny_queue_t)); - ///len/w is the evaluated minimizer numbers - kv_resize(ha_mz1_t, *p, p->n + len/w); - kv_resize(uint64_t, *mt, (int64_t)p->m); mt->n = p->n; - for (i = l = tl = buf_pos = min_pos = 0; i < len; ++i) { - int c = seq_nt4_table[(uint8_t)str[i]]; - ha_mz1_t info = dummy; - if (c < 4) { // not an ambiguous base - int z; - if (is_hpc) { - int skip_len = 1; - if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) { - for (skip_len = 2; i + skip_len < len; ++skip_len) - if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c) - break; - i += skip_len - 1; // put $i at the end of the current homopolymer run - } - tq_push(&tq, skip_len); - kmer_span += skip_len; - ///how many bases that are covered by this HPC k-mer - ///kmer_span includes at most k HPC elements - if (tq.count > k) kmer_span -= tq_shift(&tq); - } else kmer_span = l + 1 < k? l + 1 : k; - ///kmer_span should be used for HPC k-mer - ///non-HPC k-mer, kmer_span should be k - ///kmer_span is used to calculate anchor pos on reverse complementary strand - - if (k_flag != NULL) k_flag->a.a[i] = 1;///lable all useful base, which are not ignored by HPC - - kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer - kmer[1] = (kmer[1] << 1 | (c>>1)) & mask; - kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer - kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1; - if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand - z = kmer[1] < kmer[3]? 0 : 1; // strand - ++l; tl++; - if (l >= k && kmer_span < 256) { - uint64_t y; - int32_t cnt, filtered; - y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]); - cnt = hf? ha_ft_cnt(hf, y) : 0; - filtered = (cnt >= 1<<28); - if(is_unique){ - filtered = (cnt == 0); - cnt = (cnt == 1? 0:cnt); - } - if (dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i)); - if (!filtered) info.x = y, info.rid = cnt, info.pos = i, info.rev = z, info.span = kmer_span; // initially ha_mz1_t::rid keeps the k-mer count - if (k_flag != NULL) k_flag->a.a[i]++; - if (k_flag != NULL && filtered > 0) k_flag->a.a[i]++; - } - } else l = 0, tq.count = tq.front = 0, kmer_span = 0; - - buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below - buf_p[buf_pos] = l; - if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet - for (j = buf_pos + 1; j < w; ++j){ - if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos){ - kv_push(ha_mz1_t, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]); - } - } - for (j = 0; j < buf_pos; ++j){ - if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos){ - kv_push(ha_mz1_t, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]); - } - } - } - - /** - * There are three cases: - * 1. info.x <= min.x, means info is a new minimizer - * 2. info.x > min.x, info is not a new minimizer - * (1) buf_pos != min_pos, do nothing - * (2) buf_pos == min_pos, means current minimizer has moved outside the window - * **/ - ///three cases: 1. - if (mzcmp(&min, &info) >= 0) { // a new minimum; then write the old min - if (l >= w + k && min.x != UINT64_MAX){ - kv_push(ha_mz1_t, *p, min); kv_push(uint64_t, *mt, min_s); - } - min = info, min_pos = buf_pos, min_s = buf_p[buf_pos]; - } else if (buf_pos == min_pos) { // old min has moved outside the window - if (l >= w + k - 1 && min.x != UINT64_MAX){ - kv_push(ha_mz1_t, *p, min); kv_push(uint64_t, *mt, min_s); - } - ///buf_pos == min_pos, means current minimizer has moved outside the window - ///so for now we need to find a new minimizer at the current window (w k-mers) - for (j = buf_pos + 1, min = dummy; j < w; ++j) // the two loops are necessary when there are identical k-mers - if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j, min_s = buf_p[j]; // >= is important s.t. min is always the closest k-mer - for (j = 0; j <= buf_pos; ++j) - if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j, min_s = buf_p[j]; - - if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers - for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted - if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos){ - kv_push(ha_mz1_t, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]); - } - for (j = 0; j <= buf_pos; ++j) - if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos){ - kv_push(ha_mz1_t, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]); - } - } - } - if (++buf_pos == w) buf_pos = 0; - } - if (min.x != UINT64_MAX){ - kv_push(ha_mz1_t, *p, min); kv_push(uint64_t, *mt, min_s); - } - // debug_pl(str, len, w, k, is_hpc, p, hf, mt); - if (sample_dist > w) select_mz_h(p, mt, len, sample_dist, ws, k, tl); - if (dp_min_len > 0 && pt && mt) refine_sketch(p, pt, len, dp_min_len, dp_e, min_freq, mt); - for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts - p->a[i].rid = rid; +#define HA_SC_INIT(sf, HType, VType, RidBits, PosBits)\ +inline void sf##_hf_select(VType *p, int32_t si, int32_t ei, int32_t n, int32_t len, int32_t sample_dist, HType *b, int32_t force)\ +{\ + if(ei - si <= 1) return;\ + int32_t ps = si < 0? 0 : p->a[si].pos;\ + int32_t pe = ei == n? len : p->a[ei].pos;\ + int32_t j, k, st = si + 1, en = ei;\ + int32_t max_high_occ = (int32_t)((double)(pe - ps) / sample_dist + .499);\ + if (max_high_occ > MAX_MAX_HIGH_OCC)\ + max_high_occ = MAX_MAX_HIGH_OCC;\ + for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)\ + b[k] = p->a[j], b[k].pos = j; /** b[].pos keeps the index in p->a[]**/\ + ks_heapmake_##sf##_mz(k, b); /** initialize the binomial heap**/\ + for (; j < en; ++j) { /** if there are more, choose top max_high_occ**/\ + if (sf##_mz_lt(p->a[j], b[0])) { /** then update the heap**/\ + b[0] = p->a[j], b[0].pos = j;\ + ks_heapdown_##sf##_mz(0, k, b);\ + }\ + }\ + /**ks_heapsort_mz(k, b); // sorting is not needed for now**/\ + for (j = 0; j < k; ++j)\ + if (b[j].rid < pe - ps || force)\ + p->a[b[j].pos].rid = 0;\ +}\ +static inline int sf##_mzcmp_l(const VType *p, int32_t ai, int32_t bi)\ +{\ + if(ai >= 0 && bi >= 0){\ + HType *a = &(p->a[ai]), *b = &(p->a[bi]);\ + if(a->rid > 0 && b->rid > 0) return sf##_mzcmp(a, b);\ + return (a->rid == 0) - (b->rid == 0);\ + }\ + return (ai < 0) - (bi < 0);\ +}\ +int32_t sf##_qfw(VType *p, st_mt_t *mt, int32_t n, int32_t tot_l, int32_t ws, int32_t i, int32_t *mi)\ +{\ + int32_t m, si;\ + for (si = i, (*mi) = -1; i < n; i++){\ + if(GL(*mt, i) >= ws || (i+1 < n && GL(*mt, i) < ws && GL(*mt, i+1) > ws) || \ + (i+1 == n && tot_l >= ws && GL(*mt, i) < ws)){\ + for (m = si; m <= i; m++){\ + if(!A_M(*p, m)) continue;\ + if(sf##_mzcmp_l(p, *mi, m) >= 0) (*mi) = m;\ + }\ + if((*mi) >= 0 && A_M(*p, *mi)){\ + for (m = si; m <= i; m++){\ + if(!A_M(*p, m)) continue;\ + if(sf##_mzcmp_l(p, *mi, m) == 0) mt->a[m] |= 0x100000000;\ + }\ + }\ + break;\ + }\ + }\ + return i;\ +}\ +static void sf##_select_mz_h(VType *p, st_mt_t *mt, int len, int sample_dist, int32_t w, int32_t k, int32_t tot_l)\ +{ /**for high-occ minimizers, choose up to max_high_occ in each high-occ streak**/\ + int32_t i, mi = -1, si, last0 = -1, n = (int32_t)p->n, m = 0, ws = w + k - 1;\ + if (n == 0) return;\ + assert((int64_t)(n) < (int64_t)((((uint64_t)1)<a[i].rid == 0) {\ + if (i - last0 > 1) {\ + int32_t ps = last0 < 0? 0 : p->a[last0].pos;\ + int32_t pe = i == n? len : p->a[i].pos;\ + if(((int32_t)((double)(pe - ps) / sample_dist + .499)) > 0){\ + last0 = -2;\ + m++;\ + break;\ + }\ + }\ + last0 = i;\ + }\ + }\ + if (m == 0) return; /**no high-frequency k-mers; do nothing**/\ + if(last0 >= -1) goto sf##_ff;\ + i = 0;\ + i = sf##_qfw(p, mt, n, tot_l, ws, i, &mi);\ + if(i == n) goto sf##_ff;\ + for (si = 0, i++; i < n; i++){\ + for (; si < i; si++){\ + if(GL(*mt, si) + w > GL(*mt, i)) break;\ + }\ + /**a new minimum; then write the old min**/\ + if(sf##_mzcmp_l(p, i, mi) <= 0) {\ + if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;\ + mi = i;\ + }/**old min has moved outside the window**/\ + else if(si > mi){\ + if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;\ + for (m = si, mi = -1; m <= i; m++){\ + if(sf##_mzcmp_l(p, mi, m) >= 0) mi = m;\ + }\ + if(A_M(*p, mi)){\ + for (m = si; m <= i; m++){\ + if(!A_M(*p, m)) continue;\ + if(sf##_mzcmp_l(p, mi, m) == 0) mt->a[m] |= 0x100000000;\ + }\ + }\ + }\ + }\ + if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;\ + for (i = n - 1; si < n && GL(*mt, si) + w <= tot_l + 1; si++){\ + if(si > mi){\ + if(A_M(*p, mi)) mt->a[mi] |= 0x100000000;\ + for (m = si, mi = -1; m <= i; m++){\ + if(sf##_mzcmp_l(p, mi, m) >= 0) mi = m;\ + }\ + if(A_M(*p, mi)){\ + for (m = si; m <= i; m++){\ + if(!A_M(*p, m)) continue;\ + if(sf##_mzcmp_l(p, mi, m) == 0) mt->a[m] |= 0x100000000;\ + }\ + }\ + }\ + }\ + /**dbg_boundary(p, mt, w, k, tot_l);**/\ + HType b[MAX_MAX_HIGH_OCC];\ + for (i = 0, last0 = -1; i <= n; ++i) {\ + if (i == n || p->a[i].rid == 0) {\ + if (i - last0 > 1) {\ + int32_t ps = last0 < 0? 0 : p->a[last0].pos;\ + int32_t pe = i == n? len : p->a[i].pos;\ + if(((int32_t)((double)(pe - ps) / sample_dist + .499)) > 0){\ + for (m = last0 + 1, mi = 0; m < i; ++m){\ + if(mt->a[m]&0x100000000) p->a[m].rid = 0, mi++;\ + }\ + if(mi == 0) sf##_hf_select(p, last0, i, n, len, sample_dist, b, 0);\ + }\ + }\ + last0 = i;\ + }\ + }\ + sf##_ff:\ + for (i = n = 0; i < (int32_t)p->n; ++i) /**squeeze out filtered minimizers**/\ + if (p->a[i].rid == 0)\ + p->a[n++] = p->a[i];\ + p->n = n;\ +}\ +void sf##_refine_select(VType *mz, int32_t sidx, int32_t eidx, int32_t sn, int32_t min_freq, st_mt_t *mm, int32_t *rsi, int32_t *rei)\ +{\ + int32_t n = sn, m = eidx + 1 - sidx, i, k, t, mk=-1;\ + uint64_t ix, kx, ks;\ + kv_resize(uint64_t, *mm, mm->n+n*m);\ + HType *ma = mz->a + sidx;\ + uint64_t *mmt = mm->a + mm->n;\ + /**fprintf(stderr, "[M::%s::] ==> +n: %d, m: %d, sn: %d, sidx: %d, eidx: %d\n", __func__, n, m, sn, sidx, eidx);**/\ + for (i = 0; i < n; i++) /**how many selected minimizers**/\ + {\ + for (k = 0, mk = -1; k < m; k++) /**how many minimizers in total**/\ + {\ + if((int32_t)(ma[k].rid) 0){\ + for (t = k-1; t >= 0 && (ma[t].pos >= ks||(int32_t)(ma[t].rid)0?(i-1)*m+t:0xffffffff)<<32;\ + else if(ks == kx) ks |= (uint64_t)(mk>=0?i*m+mk:0xffffffff)<<32;\ + GMC(mmt, i,k,m) = ks;\ + mk = k;\ + }\ + }\ + /**fprintf(stderr, "[M::%s::] ==> ++n: %d, m: %d, sn: %d, sidx: %d, eidx: %d\n", __func__, n, m, sn, sidx, eidx);**/\ + ks = (n-1)*m + mk; ix = (uint64_t)-1; kx = 0;\ + while (ks != 0xffffffff)\ + {\ + i = ks/m; k = ks%m;\ + ks = mmt[ks]>>32;\ + /**fprintf(stderr, "i: %d, k: %d, ks: %lu\n", i, k, ks);**/\ + if(ks == 0xffffffff || (int32_t)(ks/m) == (i-1)){\ + mm->a[sidx+k] = 1;\ + ix = MIN((uint64_t)k, ix); kx = MAX((uint64_t)k, kx);\ + }\ + }\ + /**debug_refine(ma, mmt, sn, n, m, (n-1)*m + mk);**/\ + if(rsi) (*rsi) = ix + sidx;\ + if(rei) (*rei) = kx + sidx;\ +}\ +void sf##_refine_sketch(VType *p, ha_pt_t *pt, int32_t rlen, int32_t dp_min_len, float er, int32_t min_freq, st_mt_t *mt)\ +{\ + /**fprintf(stderr, "[M::%s::] ==> #########10#########, rlen: %d\n", __func__, rlen);**/\ + int32_t i, n = p->n, bd, len = MIN(rlen, dp_min_len), sublen, cnt, ei, li, ri;\ + int32_t sn = len*er + 1;\ + kv_resize(uint64_t, *mt, (int64_t)p->n);\ + mt->n = p->n; memset(mt->a, 0, sizeof(uint64_t)*p->n);\ + for (i = 0; i < n; i++) p->a[i].rid = ha_pt_cnt(pt, p->a[i].x);\ + for (i = cnt = 0, bd = -1, ei = -1; i < n; i++){\ + if((int32_t)(p->a[i].rid)a[i].pos + 1;\ + if(sublen > len) break;\ + else ei = i;\ + if((int32_t)(p->a[i].pos + 1 - p->a[i].span) > bd){\ + bd = p->a[i].pos;\ + cnt++;\ + }\ + }\ + /**fprintf(stderr, "[M::%s::] ==> +cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);**/\ + if(cnt >= sn) sf##_refine_select(p, 0, ei, sn, min_freq, mt, NULL, &li);\ + else{\ + li = i-1;\ + for (i = 0; i <= li; i++) mt->a[i] = 1;\ + }\ + if(len < rlen){\ + for (i = n-1, cnt = 0, bd = rlen+1, ei = -1; i >= 0; i--){\ + if((int32_t)(p->a[i].rid)a[i].pos + 1 - p->a[i].span);\ + if(sublen > len) break;\ + else ei = i;\ + if((int32_t)(p->a[i].pos) < bd){\ + bd = p->a[i].pos + 1 - p->a[i].span;\ + cnt++;\ + }\ + }\ + /**fprintf(stderr, "[M::%s::] ==> -cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);**/\ + if(cnt >= sn) sf##_refine_select(p, ei, n-1, sn, min_freq, mt, &ri, NULL);\ + else {\ + ri = i+1;\ + for (i = ri; i <= n-1; i++) mt->a[i] = 1;\ + }\ + /**fprintf(stderr, "[M::%s::] ==> --cnt: %d, sn: %d, ei: %d, n: %d\n", __func__, cnt, sn, ei, n);**/\ + if(ri - li >= 2){\ + li++; ri--;\ + sn = (p->a[ri].pos - p->a[li].pos + p->a[li].span)*er + 1;\ + for (i = li, cnt = 0, bd = -1; i <= ri; i++){\ + if((int32_t)(p->a[i].rid)a[i].pos + 1 - p->a[i].span) > bd){\ + bd = p->a[i].pos;\ + cnt++;\ + if(cnt >= sn) break;\ + }\ + }\ + if(cnt >= sn) sf##_refine_select(p, li, ri, sn, min_freq, mt, NULL, NULL);\ + else for (i = li; i <= ri; i++) mt->a[i] = 1;\ + }\ + }\ + /**fprintf(stderr, "[M::%s::] ==> #########20#########, p->n: %u, n: %d\n", __func__, p->n, n);**/\ + for (i = sn = 0; i < n; i++){\ + if(mt->a[i]){\ + p->a[sn] = p->a[i];\ + sn++;\ + }\ + }\ + /**if(p->n != sn) fprintf(stderr, "[M::%s::] ==> #########21#########, p->n: %u, sn: %d\n", __func__, p->n, sn);**/\ + p->n = sn;\ +}\ +/**\ + * Find symmetric (w,k)-minimizers on a DNA sequence\ + *\ + * @param str DNA sequence\ + * @param len length of $str\ + * @param w find a minimizer for every $w consecutive k-mers\ + * @param k k-mer size\ + * @param rid reference ID; will be copied to the output $p array\ + * @param is_hpc homopolymer-compressed or not\ + * @param p minimizers\ + */\ +void sf##_ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, VType *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, ha_pt_t *pt, int min_freq, int32_t dp_min_len, float dp_e, st_mt_t *mt, int32_t ws, int32_t is_unique)\ +{ /**in default, w = 51, k = 51, is_hpc = 1**/\ + extern void *ha_ct_table;\ + static const HType dummy = { UINT64_MAX, (((uint64_t)1)< 0 && (int64_t)(len) < (int64_t)((((uint64_t)1)< 0 && w < 256) && (k > 0 && k <= 63));\ + if (dbg_ct != NULL) dbg_ct->a.n = 0;\ + if (k_flag != NULL) {\ + kv_resize(uint8_t, k_flag->a, (uint64_t)len);\ + k_flag->a.n = len;\ + memset(k_flag->a.a, 0, k_flag->a.n);\ + }\ + memset(buf, 0xff, w * sizeof(HType));\ + memset(&tq, 0, sizeof(tiny_queue_t));\ + /**len/w is the evaluated minimizer numbers**/\ + kv_resize(HType, *p, p->n + len/w);\ + kv_resize(uint64_t, *mt, (int64_t)p->m); mt->n = p->n;\ + for (i = l = tl = buf_pos = min_pos = 0; i < len; ++i) {\ + int c = seq_nt4_table[(uint8_t)str[i]];\ + HType info = dummy;\ + if (c < 4) { /**not an ambiguous base**/\ + int z;\ + if (is_hpc) {\ + int skip_len = 1;\ + if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) {\ + for (skip_len = 2; i + skip_len < len; ++skip_len)\ + if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c)\ + break;\ + i += skip_len - 1; /**put $i at the end of the current homopolymer run**/\ + }\ + tq_push(&tq, skip_len);\ + kmer_span += skip_len;\ + /**how many bases that are covered by this HPC k-mer\ + kmer_span includes at most k HPC elements**/\ + if (tq.count > k) kmer_span -= tq_shift(&tq);\ + } else kmer_span = l + 1 < k? l + 1 : k;\ + /**kmer_span should be used for HPC k-mer\ + non-HPC k-mer, kmer_span should be k\ + kmer_span is used to calculate anchor pos on reverse complementary strand**/\ + if (k_flag != NULL) k_flag->a.a[i] = 1;/**lable all useful base, which are not ignored by HPC**/\ + kmer[0] = (kmer[0] << 1 | (c&1)) & mask;/**forward k-mer**/\ + kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;\ + kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; /**reverse k-mer**/\ + kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1;\ + if (kmer[1] == kmer[3]) continue; /** skip "symmetric k-mers" as we don't know it strand**/\ + z = kmer[1] < kmer[3]? 0 : 1; /** strand**/\ + ++l; tl++;\ + if (l >= k && kmer_span < 256) {\ + uint64_t y;\ + int32_t cnt, filtered;\ + y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);\ + cnt = hf? ha_ft_cnt(hf, y) : 0;\ + filtered = (cnt >= 1<<28);\ + if(is_unique){\ + filtered = (cnt < is_unique);\ + cnt = (cnt == is_unique? 0:cnt);\ + }\ + if (dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i));\ + if (!filtered) info.x = y, info.rid = cnt, info.pos = i, info.rev = z, info.span = kmer_span; /** initially ha_mz1_t::rid keeps the k-mer count**/\ + if (k_flag != NULL) k_flag->a.a[i]++;\ + if (k_flag != NULL && filtered > 0) k_flag->a.a[i]++;\ + }\ + } else l = 0, tq.count = tq.front = 0, kmer_span = 0;\ + buf[buf_pos] = info; /**need to do this here as appropriate buf_pos and buf[buf_pos] are needed below**/\ + buf_p[buf_pos] = l;\ + if (l == w + k - 1 && min.x != UINT64_MAX) { /**special case for the first window - because identical k-mers are not stored yet**/\ + for (j = buf_pos + 1; j < w; ++j){\ + if (sf##_mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos){\ + kv_push(HType, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);\ + }\ + }\ + for (j = 0; j < buf_pos; ++j){\ + if (sf##_mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos){\ + kv_push(HType, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);\ + }\ + }\ + }\ + /**\ + * There are three cases:\ + * 1. info.x <= min.x, means info is a new minimizer\ + * 2. info.x > min.x, info is not a new minimizer\ + * (1) buf_pos != min_pos, do nothing\ + * (2) buf_pos == min_pos, means current minimizer has moved outside the window\ + * **/\ + /**three cases: 1.**/\ + if (sf##_mzcmp(&min, &info) >= 0) { /**a new minimum; then write the old min**/\ + if (l >= w + k && min.x != UINT64_MAX){\ + kv_push(HType, *p, min); kv_push(uint64_t, *mt, min_s);\ + }\ + min = info, min_pos = buf_pos, min_s = buf_p[buf_pos];\ + } else if (buf_pos == min_pos) { /**old min has moved outside the window**/\ + if (l >= w + k - 1 && min.x != UINT64_MAX){\ + kv_push(HType, *p, min); kv_push(uint64_t, *mt, min_s);\ + }\ + /**buf_pos == min_pos, means current minimizer has moved outside the window\ + so for now we need to find a new minimizer at the current window (w k-mers)**/\ + for (j = buf_pos + 1, min = dummy; j < w; ++j) /**the two loops are necessary when there are identical k-mers**/\ + if (sf##_mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j, min_s = buf_p[j]; /** >= is important s.t. min is always the closest k-mer**/\ + for (j = 0; j <= buf_pos; ++j)\ + if (sf##_mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j, min_s = buf_p[j];\ + if (l >= w + k - 1 && min.x != UINT64_MAX) { /**write identical k-mers**/\ + for (j = buf_pos + 1; j < w; ++j) /**these two loops make sure the output is sorted**/\ + if (sf##_mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos){\ + kv_push(HType, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);\ + }\ + for (j = 0; j <= buf_pos; ++j)\ + if (sf##_mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos){\ + kv_push(HType, *p, buf[j]); kv_push(uint64_t, *mt, buf_p[j]);\ + }\ + }\ + }\ + if (++buf_pos == w) buf_pos = 0;\ + }\ + if (min.x != UINT64_MAX){\ + kv_push(HType, *p, min); kv_push(uint64_t, *mt, min_s);\ + }\ + /**debug_pl(str, len, w, k, is_hpc, p, hf, mt);**/\ + if (sample_dist > w) sf##_select_mz_h(p, mt, len, sample_dist, ws, k, tl);\ + if (dp_min_len > 0 && pt && mt) sf##_refine_sketch(p, pt, len, dp_min_len, dp_e, min_freq, mt);\ + for (i = 0; i < (int)p->n; ++i) /**populate .rid as this was keeping counts**/\ + p->a[i].rid = rid;\ } +HA_SC_INIT(mz1, ha_mz1_t, ha_mz1_v, 28, 27) +HA_SC_INIT(mz2, ha_mzl_t, ha_mzl_v, 31, 32) \ No newline at end of file