diff --git a/Assembly.cpp b/Assembly.cpp index 2821086..d9c2bdc 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1023,10 +1023,10 @@ void ha_ec(int64_t round, int num_pround, int des_idx, uint64_t *tot_b, uint64_t if(round == asm_opt.number_of_round-1 && asm_opt.is_dbg_het_cnt) CALLOC(het_cnt, R_INF.total_reads); if (r_out) { - write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name); + write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, 0); if((asm_opt.flag & HA_F_VERBOSE_GFA) && (asm_opt.bin_only == 1)) exit(1);///just for debug } - if (w_tmp) tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, round, asm_opt.number_of_round, 0); + if (w_tmp) tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, round, asm_opt.number_of_round, 0, 0); // Output_corrected_fastq(); @@ -1112,7 +1112,7 @@ void ha_overlap_and_correct(int round) // fprintf(stderr, "[M::%s::%.3f] ==> chaining\n", __func__, yak_realtime_0()-tt0); // exit(1); - if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name); + if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, 0); ha_pt_destroy(ha_idx); ha_idx = NULL; @@ -1727,6 +1727,7 @@ void Output_PAF() fclose(output_file); fprintf(stderr, "PAF has been written.\n"); + exit(1); } @@ -1952,12 +1953,12 @@ void ha_overlap_final(void) asm_opt.het_cov = het_cov; } -void ha_ec_ff(int renew_idx) +void ha_ec_ff(int renew_idx, int8_t pre_load_idx, uint64_t w_tmp) { int hom_cov, het_cov; ha_flt_tab_hp = ha_idx_hp = NULL; - if(ha_idx && renew_idx) { + if((ha_idx) && (renew_idx) && (!pre_load_idx)) { ha_pt_destroy(ha_idx); ha_idx = NULL; } @@ -1966,6 +1967,8 @@ void ha_ec_ff(int renew_idx) asm_opt.hom_cov = hom_cov; asm_opt.het_cov = het_cov; } + if (w_tmp) tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, asm_opt.number_of_round, asm_opt.number_of_round, 0, 1); + cal_ov_r(asm_opt.thread_num, R_INF.total_reads, renew_idx); if(asm_opt.write_pos_idx) { @@ -2076,7 +2079,7 @@ int ha_assemble(void) // debug_mc_gg_t(MC_NAME, 0, 0); // quick_debug_phasing(MC_NAME); extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o); - int r, r0 = -1, hom_cov = -1, ovlp_loaded = 0; uint64_t tot_b, tot_e; + int r, r0 = -1, hom_cov = -1, ovlp_loaded = 0; uint64_t tot_b, tot_e; int8_t pre_load_ff = 0; if ((asm_opt.load_index_from_disk) && (asm_opt.dbg_ec_rr < 0) && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) { ovlp_loaded = 1; fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage()); @@ -2093,13 +2096,13 @@ int ha_assemble(void) } if (!ovlp_loaded) { ha_flt_tab = ha_idx = NULL; - if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name), load_ct_index(&ha_ct_table, asm_opt.output_file_name); + if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, 0), load_ct_index(&ha_ct_table, asm_opt.output_file_name); r = ha_idx?asm_opt.number_of_round-1:0; if((!ha_idx) && (asm_opt.restart)) { - r = asm_opt.number_of_round - 1; + r = asm_opt.number_of_round/** - 1**/; if(asm_opt.dbg_ec_rr >= 0) r = asm_opt.dbg_ec_rr; for (; r >= 0; --r) { - if(tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, r, asm_opt.number_of_round, 1)) { + if(tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, r, asm_opt.number_of_round, 1, ((r==asm_opt.number_of_round)?(1):(0)))) { load_ct_index(&ha_ct_table, asm_opt.output_file_name); r0 = r; break; } @@ -2108,7 +2111,11 @@ int ha_assemble(void) fprintf(stderr, "[E::%s] no matching debug error-correction bins found\n", __func__); exit(1); } - if(r < 0) r = 0; + if(r < 0) { + r = 0; + } else if(r == asm_opt.number_of_round) { + pre_load_ff = 1; + } } // construct hash table for high occurrence k-mers @@ -2138,8 +2145,9 @@ int ha_assemble(void) // overlap between corrected reads ha_opt_reset_to_round(&asm_opt, asm_opt.number_of_round); // ha_overlap_final(); - ha_ec_ff(1/**0**/); + ha_ec_ff(1/**0**/, pre_load_ff, ((r > r0) && (asm_opt.restart))?1:0); fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb()); + if(asm_opt.dbg_ec_rr >= 0) exit(1); // fprintf(stderr, "\n[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb()); // ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads); if(!(asm_opt.write_pos_idx)) { @@ -2193,7 +2201,7 @@ int ha_assemble_pair(void) // Output_corrected_reads(); exit(0); ha_flt_tab = ha_idx = NULL; r = asm_opt.number_of_round - 1; - if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name), load_ct_index(&ha_ct_table, asm_opt.output_file_name); + if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, 0), load_ct_index(&ha_ct_table, asm_opt.output_file_name); // construct hash table for high occurrence k-mers if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL) { diff --git a/CommandLines.cpp b/CommandLines.cpp index debc942..c3edfdc 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -12,6 +12,11 @@ KSEQ_INIT(gzFile, gzread) #define DEFAULT_OUTPUT "hifiasm.asm" +#define RER_H_HiFi 0.06 +#define RER_H_ONT 0.08 +#define RER_N_HiFi 0.03 +#define RER_N_ONT 0.05 +#define MIN_R_COV 3 hifiasm_opt_t asm_opt; @@ -94,7 +99,14 @@ static ko_longopt_t long_options[] = { { "hyb-syn", ko_required_argument, 376}, { "simd-m", ko_required_argument, 377}, { "del-hf", ko_no_argument, 378}, - { "dbg-rr", ko_required_argument, 379}, + { "dbg-rr", ko_required_argument, 379},///start from 0/1/2 step of correction + { "re-aln", ko_required_argument, 380}, + { "syn", ko_required_argument, 381}, + { "rec", ko_required_argument, 382}, + { "ref", ko_required_argument, 383}, + { "h_rec", ko_required_argument, 384}, + { "h_ref", ko_required_argument, 385}, + { "ret", ko_no_argument, 386}, // { "path-round", ko_required_argument, 348}, { 0, 0, 0 } }; @@ -106,6 +118,7 @@ double Get_T(void) return t.tv_sec+t.tv_usec/1000000.0; } + void Print_H(hifiasm_opt_t* asm_opt) { fprintf(stderr, "Usage: hifiasm [options] <...>\n"); @@ -140,7 +153,34 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " discard overlaps supported by chn_occ); fprintf(stderr, " --ec-only error correction only; disable overlapping and assembly\n"); fprintf(stderr, " --simd-m use SIMD acceleration when supported: AVX-512 (2), AVX2 (1), or non-SIMD (0)\n"); + fprintf(stderr, " --re-aln keep raw reads for error correction [%d]\n", asm_opt->realn_raw); + fprintf(stderr, " --syn rescue overlaps for error correction [%d]\n", asm_opt->post_syn); + + + fprintf(stderr, " Recurrent sequencing-error filtering:\n"); + fprintf(stderr, " --ret enable recurrent sequencing-error filtering; disabled by default\n"); + fprintf(stderr, " thresholds are controlled by --rec/--ref and --h_rec/--h_ref\n"); + fprintf(stderr, " --rec/--ref and --h_rec/--h_ref are ignored unless --ret is set\n"); + fprintf(stderr, " --rec INT\n"); + fprintf(stderr, " discard candidate variants supported by recurrent_err_normal_min); + fprintf(stderr, " --ref FLOAT\n"); + fprintf(stderr, " discard candidate variants supported by recurrent_err_hpc_min); + fprintf(stderr, " --h_ref FLOAT\n"); + fprintf(stderr, " discard candidate variants near homopolymers if supported by clean_round); fprintf(stderr, " -m INT pop bubbles of large_pop_bubble_size); @@ -365,8 +405,8 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->dp_e = 0.0025; asm_opt->hg_size = -1; asm_opt->kpt_rate = -1; - asm_opt->infor_cov = 3; - asm_opt->s_hap_cov = 3; + asm_opt->infor_cov = MIN_R_COV; + asm_opt->s_hap_cov = MIN_R_COV; asm_opt->ul_error_rate = 0.2/**0.15**/; asm_opt->ul_error_rate_low = 0.1; asm_opt->ul_error_rate_hpc = 0.2; @@ -447,6 +487,22 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->del_hf = 0; asm_opt->dbg_ec_rr = -1; + + asm_opt->realn_raw = 0; + asm_opt->post_syn = 1; + + asm_opt->recurrent_err_normal_min = MIN_R_COV; + asm_opt->recurrent_err_normal_min_set = -1; + asm_opt->recurrent_err_normal_rat = RER_N_HiFi; + asm_opt->recurrent_err_normal_rat_set = -1; + + asm_opt->recurrent_err_hpc_min = MIN_R_COV; + asm_opt->recurrent_err_hpc_min_set = -1; + asm_opt->recurrent_err_hpc_rat = RER_H_HiFi; + asm_opt->recurrent_err_hpc_rat_set = -1; + + asm_opt->recurrent_err_test = 0; + } void destory_enzyme(enzyme* f) @@ -837,6 +893,26 @@ int check_option(hifiasm_opt_t* asm_opt) return 0; } + if((asm_opt->realn_raw != 0) && (asm_opt->realn_raw != 1)) { + fprintf(stderr, "[ERROR] [--re-aln] must be 0/1\n"); + return 0; + } + + if((asm_opt->post_syn != 0) && (asm_opt->post_syn != 1)) { + fprintf(stderr, "[ERROR] [--syn] must be 0/1\n"); + return 0; + } + + if(asm_opt->recurrent_err_normal_rat > 1.0) { + fprintf(stderr, "[ERROR] [--ref] must be >= 0 && <= 1.0; -1 to disable\n"); + return 0; + } + + if(asm_opt->recurrent_err_hpc_rat > 1.0) { + fprintf(stderr, "[ERROR] [--h_ref] must be >= 0 && <= 1.0; -1 to disable\n"); + return 0; + } + return 1; } @@ -1126,6 +1202,24 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) asm_opt->del_hf = 1; } else if (c == 379) { asm_opt->dbg_ec_rr = atoi(opt.arg); + } else if (c == 380) { + asm_opt->realn_raw = atoi(opt.arg); + } else if (c == 381) { + asm_opt->post_syn = atoi(opt.arg); + } else if (c == 382) { + asm_opt->recurrent_err_normal_min = atoi(opt.arg); + asm_opt->recurrent_err_normal_min_set = 1; + } else if (c == 383) { + asm_opt->recurrent_err_normal_rat = atof(opt.arg); + asm_opt->recurrent_err_normal_rat_set = 1; + } else if (c == 384) { + asm_opt->recurrent_err_hpc_min = atoi(opt.arg); + asm_opt->recurrent_err_hpc_min_set = 1; + } else if (c == 385) { + asm_opt->recurrent_err_hpc_rat = atof(opt.arg); + asm_opt->recurrent_err_hpc_rat_set = 1; + } else if (c == 386) { + asm_opt->recurrent_err_test = 1; } 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); } @@ -1167,5 +1261,35 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) asm_opt->rl_cut = -1; asm_opt->sc_cut = 1; } + if(asm_opt->recurrent_err_normal_min_set != 1) { + asm_opt->recurrent_err_normal_min = MIN_R_COV; + } + if(asm_opt->recurrent_err_normal_rat_set != 1) { + asm_opt->recurrent_err_normal_rat = ((asm_opt->is_ont)?(RER_N_ONT):(RER_N_HiFi)); + } + if(asm_opt->recurrent_err_hpc_min_set != 1) { + asm_opt->recurrent_err_hpc_min = MIN_R_COV; + } + if(asm_opt->recurrent_err_hpc_rat_set != 1) { + asm_opt->recurrent_err_hpc_rat = ((asm_opt->is_ont)?(RER_H_ONT):(RER_H_HiFi)); + } + + if(asm_opt->recurrent_err_normal_min < 0) asm_opt->recurrent_err_normal_min = -1; + if(asm_opt->recurrent_err_normal_rat < 0) asm_opt->recurrent_err_normal_rat = -1; + if(asm_opt->recurrent_err_hpc_min < 0) asm_opt->recurrent_err_hpc_min = -1; + if(asm_opt->recurrent_err_hpc_rat < 0) asm_opt->recurrent_err_hpc_rat = -1; + + if(asm_opt->recurrent_err_test == 0) { + asm_opt->recurrent_err_normal_min = -1; + asm_opt->recurrent_err_normal_min_set = -1; + asm_opt->recurrent_err_normal_rat = -1; + asm_opt->recurrent_err_normal_rat_set = -1; + + asm_opt->recurrent_err_hpc_min = -1; + asm_opt->recurrent_err_hpc_min_set = -1; + asm_opt->recurrent_err_hpc_rat = -1; + asm_opt->recurrent_err_hpc_rat_set = -1; + } + return check_option(asm_opt); } diff --git a/CommandLines.h b/CommandLines.h index 639a6ce..48b4185 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.25.1-r920" +#define HA_VERSION "0.25.1-r933" #define VERBOSE 0 @@ -206,6 +206,24 @@ typedef struct { int8_t del_hf; int64_t dbg_ec_rr; + + int8_t realn_raw; + int8_t post_syn; + + int64_t recurrent_err_normal_min; + int8_t recurrent_err_normal_min_set; + + double recurrent_err_normal_rat; + int8_t recurrent_err_normal_rat_set; + + int64_t recurrent_err_hpc_min; + int8_t recurrent_err_hpc_min_set; + + double recurrent_err_hpc_rat; + int8_t recurrent_err_hpc_rat_set; + + int8_t recurrent_err_test; + } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Correct.cpp b/Correct.cpp index 2b39ad7..19677ab 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -24,6 +24,7 @@ #define HF_W_P 4 #define OT_W_P 1 #define h0_w_p 1 +#define HPC_COMP_LEN 6 #define generic_key(x) (x) KRADIX_SORT_INIT(b32, uint32_t, generic_key, 4) @@ -8848,7 +8849,7 @@ void prt_sub_cigar(overlap_region* z, uint64_t str_l, uint64_t site, uint64_t wi #define is_st_bs(s, rr, mm) (((mm) != ((uint64_t)-1)) && (((s).overlap_num + mm) >= ((s).occ_0)) && ((((s).occ_0*(rr) + (s).overlap_num)) >= ((s).occ_0))) void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, uint64_t multi_check, double st_rate, uint64_t st_max, - int64_t hap_cov_match, int64_t hap_cov_unmatch) + int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t flag_hf_ov_cut) { // fprintf(stderr, "[M::%s::] Done\n", __func__); if(hap->length == 0) return; @@ -9021,33 +9022,70 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio } } - for (k = 1, l = 0; k <= hap->length; ++k) { - if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { - ii = hap->list[l].overlapID; - if(overlap_list->list[ii].is_match==2) { - overlap_list->list[ii].strong = 1; - overlap_list->mapped_overlaps_length -= - overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; - } else if(overlap_list->list[ii].is_match==1) { - for (i = l; i < k; i++) { - if(hh_tp(hap->list[i])==1 || hh_tp(hap->list[i])==0) { - s = &(hap->snp_stat.a[hap->list[i].overlapSite]); - if(s->score == 1 && (!(s->occ_0 < 2 || s->occ_1 < 2)) && (!(is_st_bs((*s), st_rate, st_max)))) { - overlap_list->list[ii].strong = 1; - if(hh_tp(hap->list[i])==1) { - overlap_list->list[ii].is_match = 2; - overlap_list->mapped_overlaps_length -= - overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; - break; + if(flag_hf_ov_cut >= INT64_MAX) { + for (k = 1, l = 0; k <= hap->length; ++k) { + if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { + ii = hap->list[l].overlapID; + if(overlap_list->list[ii].is_match==2) { + overlap_list->list[ii].strong = 1; + overlap_list->mapped_overlaps_length -= + overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; + } else if(overlap_list->list[ii].is_match==1) { + for (i = l; i < k; i++) { + if(hh_tp(hap->list[i])==1 || hh_tp(hap->list[i])==0) { + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->score == 1 && (!(s->occ_0 < 2 || s->occ_1 < 2)) && (!(is_st_bs((*s), st_rate, st_max)))) { + overlap_list->list[ii].strong = 1; + if(hh_tp(hap->list[i])==1) { + overlap_list->list[ii].is_match = 2; + overlap_list->mapped_overlaps_length -= + overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; + break; + } + } + } + + } + } + l = k; + } + } + } else { + int64_t n_asnp = 0, n_hits = 0; + for (k = 1, l = o = 0; k <= hap->length; ++k) { + if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { + ii = hap->list[l].overlapID; n_asnp = n_hits = 0; + for (; o <= ii; o++) overlap_list->list[o].without_large_indel = 0; + if((overlap_list->list[ii].is_match == 1) || (overlap_list->list[ii].is_match == 2)) { + for (i = l; i < k; i++) { + if(hh_tp(hap->list[i])==1 || hh_tp(hap->list[i])==0) { + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->score == 1 && (!(s->occ_0 < 2 || s->occ_1 < 2)) && (!(is_st_bs((*s), st_rate, st_max)))) { + n_hits++; + if(hh_tp(hap->list[i])==1) { + n_asnp++; + if(n_asnp > flag_hf_ov_cut) break; + } } } } - + if((n_asnp) || (overlap_list->list[ii].is_match == 2)) { + overlap_list->list[ii].strong = 1; + overlap_list->list[ii].is_match = 2; + overlap_list->mapped_overlaps_length -= + overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; + } + if((n_asnp <= flag_hf_ov_cut) && (overlap_list->list[ii].is_match == 2)) { + overlap_list->list[ii].without_large_indel = 1; + } + if(n_hits) { + overlap_list->list[ii].strong = 1; + } } + l = k; } - l = k; - } - } + } + } } @@ -9927,7 +9965,7 @@ void generate_haplotypes_naive_HiFi_adv(haplotype_evdience_alloc* hap, overlap_r -void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, double st_rate, uint64_t st_max, asg32_v *b32, uint64_t rid, int64_t hap_cov_match, int64_t hap_cov_unmatch) +void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, double st_rate, uint64_t st_max, asg32_v *b32, uint64_t rid, int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t flag_hf_ov_cut) { if(hap->length == 0) return; uint64_t k, l, i, o, obs, ii, m_snp_stat, m_list, m_off; @@ -9985,7 +10023,9 @@ void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overla s = &(hap->snp_stat.a[hap->list[i].overlapSite]); assert(s->site == hap->list[i].site); zf = is_plus_sc_obs(s, st_rate, st_max, hap_cov_match, hap_cov_unmatch); - // if(overlap_list->list[hap->list[l].overlapID].y_id == 5305) { + /**if(overlap_list->list[hap->list[l].overlapID].y_id == 3646295 || overlap_list->list[hap->list[l].overlapID].y_id == 3203512 + || overlap_list->list[hap->list[l].overlapID].y_id == 3149588) **/ + // if(overlap_list->list[hap->list[l].overlapID].y_id == 2198541) { // fprintf(stderr, "[M::%s]\t%.*s\ts->site::%u\ts->occ_0::%u\ts->occ_1::%u\ts->overlap_num::%u\n", __func__, // (int)Get_NAME_LENGTH(R_INF, overlap_list->list[hap->list[l].overlapID].y_id), Get_NAME(R_INF, overlap_list->list[hap->list[l].overlapID].y_id), // s->site, s->occ_0, s->occ_1, s->overlap_num); @@ -10110,33 +10150,70 @@ void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overla // fprintf(stderr, "sss[M::%s::]\tsite::%u\tocc0::%u\tocc1::%u\tocc2::%u\tsc::%d\n", __func__, s->site, s->occ_0, s->occ_1, s->occ_2, s->score); // } - for (k = 1, l = 0; k <= hap->length; ++k) { - if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { - ii = hap->list[l].overlapID; - if(overlap_list->list[ii].is_match==2) { - overlap_list->list[ii].strong = 1; - overlap_list->mapped_overlaps_length -= - overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; - } else if(overlap_list->list[ii].is_match==1) { - for (i = l; i < k; i++) { - if(hh_tp(hap->list[i])==1 || hh_tp(hap->list[i])==0) { - s = &(hap->snp_stat.a[hap->list[i].overlapSite]); - if(s->score > 0 && (!(s->occ_0 < 2 || s->occ_1 < 2)) && (!(is_st_bs((*s), st_rate, st_max)))) { - overlap_list->list[ii].strong = 1; - if(hh_tp(hap->list[i])==1) { - // fprintf(stderr, "-2-[M::%s]\t%.*s\n", __func__, (int)Get_NAME_LENGTH(R_INF, overlap_list->list[ii].y_id), Get_NAME(R_INF, overlap_list->list[ii].y_id)); - overlap_list->list[ii].is_match = 2; - overlap_list->mapped_overlaps_length -= - overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; - break; + if(flag_hf_ov_cut >= INT64_MAX) { + for (k = 1, l = 0; k <= hap->length; ++k) { + if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { + ii = hap->list[l].overlapID; + if(overlap_list->list[ii].is_match==2) { + overlap_list->list[ii].strong = 1; + overlap_list->mapped_overlaps_length -= + overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; + } else if(overlap_list->list[ii].is_match==1) { + for (i = l; i < k; i++) { + if(hh_tp(hap->list[i])==1 || hh_tp(hap->list[i])==0) { + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->score > 0 && (!(s->occ_0 < 2 || s->occ_1 < 2)) && (!(is_st_bs((*s), st_rate, st_max)))) { + overlap_list->list[ii].strong = 1; + if(hh_tp(hap->list[i])==1) { + // fprintf(stderr, "-2-[M::%s]\t%.*s\n", __func__, (int)Get_NAME_LENGTH(R_INF, overlap_list->list[ii].y_id), Get_NAME(R_INF, overlap_list->list[ii].y_id)); + overlap_list->list[ii].is_match = 2; + overlap_list->mapped_overlaps_length -= + overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; + break; + } + } + } + + } + } + l = k; + } + } + } else { + int64_t n_asnp = 0, n_hits = 0; + for (k = 1, l = o = 0; k <= hap->length; ++k) { + if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { + ii = hap->list[l].overlapID; n_asnp = n_hits = 0; + for (; o <= ii; o++) overlap_list->list[o].without_large_indel = 0; + if((overlap_list->list[ii].is_match == 1) || (overlap_list->list[ii].is_match == 2)) { + for (i = l; i < k; i++) { + if(hh_tp(hap->list[i])==1 || hh_tp(hap->list[i])==0) { + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->score > 0 && (!(s->occ_0 < 2 || s->occ_1 < 2)) && (!(is_st_bs((*s), st_rate, st_max)))) { + n_hits++; + if(hh_tp(hap->list[i])==1) { + n_asnp++; + if(n_asnp > flag_hf_ov_cut) break; + } } } } - + if((n_asnp) || (overlap_list->list[ii].is_match == 2)) { + overlap_list->list[ii].strong = 1; + overlap_list->list[ii].is_match = 2; + overlap_list->mapped_overlaps_length -= + overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; + } + if((n_asnp <= flag_hf_ov_cut) && (overlap_list->list[ii].is_match == 2)) { + overlap_list->list[ii].without_large_indel = 1; + } + if(n_hits) { + overlap_list->list[ii].strong = 1; + } } + l = k; } - l = k; - } + } } copy_asg_arr(hap->snp_srt, buf); @@ -10328,7 +10405,7 @@ void generate_haplotypes_naive_HiFi_adv_back(haplotype_evdience_alloc* hap, over void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, uint64_t multi_check, double st_rate, uint64_t st_max, uint64_t snp_dis, int64_t snp_cut, - int64_t hap_cov_match, int64_t hap_cov_unmatch) + int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t flag_hf_ov_cut) { // fprintf(stderr, "-0-[M::%s::] Done\n", __func__); if(hap->length == 0) return; @@ -10569,32 +10646,69 @@ void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_al } } - for (k = 1, l = 0; k <= hap->length; ++k) { - if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { - ii = hap->list[l].overlapID; - if(overlap_list->list[ii].is_match==2) { - overlap_list->list[ii].strong = 1; - overlap_list->mapped_overlaps_length -= - overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; - } else if(overlap_list->list[ii].is_match==1) { - for (i = l; i < k; i++) { - if(hh_tp(hap->list[i])==1 || hh_tp(hap->list[i])==0) { - s = &(hap->snp_stat.a[hap->list[i].overlapSite]); - if((s->id == UINT32_MAX/**s->score == 1**/) && (!(s->occ_0 < 2 || s->occ_1 < 2)) && (!(is_st_bs((*s), st_rate, st_max)))) { - overlap_list->list[ii].strong = 1; - if(hh_tp(hap->list[i])==1) { - overlap_list->list[ii].is_match = 2; - overlap_list->mapped_overlaps_length -= - overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; - break; + if(flag_hf_ov_cut >= INT64_MAX) { + for (k = 1, l = 0; k <= hap->length; ++k) { + if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { + ii = hap->list[l].overlapID; + if(overlap_list->list[ii].is_match==2) { + overlap_list->list[ii].strong = 1; + overlap_list->mapped_overlaps_length -= + overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; + } else if(overlap_list->list[ii].is_match==1) { + for (i = l; i < k; i++) { + if(hh_tp(hap->list[i])==1 || hh_tp(hap->list[i])==0) { + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if((s->id == UINT32_MAX/**s->score == 1**/) && (!(s->occ_0 < 2 || s->occ_1 < 2)) && (!(is_st_bs((*s), st_rate, st_max)))) { + overlap_list->list[ii].strong = 1; + if(hh_tp(hap->list[i])==1) { + overlap_list->list[ii].is_match = 2; + overlap_list->mapped_overlaps_length -= + overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; + break; + } + } + } + + } + } + l = k; + } + } + } else { + int64_t n_asnp = 0, n_hits = 0; + for (k = 1, l = o = 0; k <= hap->length; ++k) { + if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { + ii = hap->list[l].overlapID; n_asnp = n_hits = 0; + for (; o <= ii; o++) overlap_list->list[o].without_large_indel = 0; + if((overlap_list->list[ii].is_match == 1) || (overlap_list->list[ii].is_match == 2)) { + for (i = l; i < k; i++) { + if(hh_tp(hap->list[i])==1 || hh_tp(hap->list[i])==0) { + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if((s->id == UINT32_MAX/**s->score == 1**/) && (!(s->occ_0 < 2 || s->occ_1 < 2)) && (!(is_st_bs((*s), st_rate, st_max)))) { + n_hits++; + if(hh_tp(hap->list[i])==1) { + n_asnp++; + if(n_asnp > flag_hf_ov_cut) break; + } } } } - + if((n_asnp) || (overlap_list->list[ii].is_match == 2)) { + overlap_list->list[ii].strong = 1; + overlap_list->list[ii].is_match = 2; + overlap_list->mapped_overlaps_length -= + overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; + } + if((n_asnp <= flag_hf_ov_cut) && (overlap_list->list[ii].is_match == 2)) { + overlap_list->list[ii].without_large_indel = 1; + } + if(n_hits) { + overlap_list->list[ii].strong = 1; + } } + l = k; } - l = k; - } + } } for (k = 0; k < hap->snp_stat.n; k++) { @@ -10609,7 +10723,7 @@ void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_al } -void generate_haplotypes_sv(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, uint64_t rid, int64_t hap_cov_match, int64_t hap_cov_unmatch) +void generate_haplotypes_sv(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, uint64_t rid, int64_t hap_cov_match, int64_t hap_cov_unmatch, uint8_t flag_hf_ov) { uint64_t k, l, i, o, ii; int64_t z; SnpStats *s = NULL; @@ -10714,6 +10828,7 @@ void generate_haplotypes_sv(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list->list[ii].strong = 1; overlap_list->mapped_overlaps_length -= overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; + if(flag_hf_ov) overlap_list->list[ii].without_large_indel = 0; } else if(overlap_list->list[ii].is_match==1) { for (i = l; i < k; i++) { if(hh_tp(hap->list[i])==1 || hh_tp(hap->list[i])==0) { @@ -10724,6 +10839,7 @@ void generate_haplotypes_sv(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list->list[ii].is_match = 2; overlap_list->mapped_overlaps_length -= overlap_list->list[ii].x_pos_e + 1 - overlap_list->list[ii].x_pos_s; + if(flag_hf_ov) overlap_list->list[ii].without_large_indel = 0; // fprintf(stderr, "[M::%s] rid::%u\t%.*s\n", __func__, overlap_list->list[ii].y_id, (int)Get_NAME_LENGTH(R_INF, overlap_list->list[ii].y_id), Get_NAME(R_INF, overlap_list->list[ii].y_id)); break; } @@ -11402,6 +11518,7 @@ void gen_rphase_dp0_single_path_hybrid_0_multi(SnpStats *a, int64_t an, haplotyp rz[rn++] = rk; /**ii[rk] = 1;**/ rk = p[rk];///r826 } // assert(rn > 0); + // fprintf(stderr, "\n+[M::%s]\trn::%ld\thf_only::%lu\tcut_rate::%f\tcc0::%lu\tcut_bd::%lu\n", __func__, rn, hf_only, cut_rate, cc0, cut_bd); if(qual_a) { krn = rn; @@ -11438,7 +11555,8 @@ void gen_rphase_dp0_single_path_hybrid_0_multi(SnpStats *a, int64_t an, haplotyp } } - // fprintf(stderr, "+[M::%s]\tsite::%u\tsc::%d\tn0::%u\tn1::%u\trn::%ld\tkrn::%ld\n", __func__, a[rz[rk]].site, a[rz[rk]].score, a[rz[rk]].occ_0, a[rz[rk]].occ_1, rn, krn); + // fprintf(stderr, "+[M::%s]\tsite::%u\tsc::%d\tn0::%u\tn1::%u\trn::%ld\tkrn::%ld\tb0l::%ld\tb0h::%ld\tb1l::%ld\tb1h::%ld\tcc::%lu\ttot_occ::%u\n", __func__, a[rz[rk]].site, a[rz[rk]].score, a[rz[rk]].occ_0, a[rz[rk]].occ_1, rn, krn, + // b0l, b0h, b1l, b1h, cc, a[rz[rk]].occ_2 + a[rz[rk]].occ_1 + a[rz[rk]].occ_0); // if(hf_only) { // fprintf(stderr, "pos::%u\tsc::%d\tn0::%u\tn1::%u\tn2::%u\tk::%lu\tb0l::%ld\tb0h::%ld\tb1l::%ld\tb1h::%ld\n", a[ra[rn0 + i]].site, a[ra[rn0 + i]].score, a[ra[rn0 + i]].occ_0, a[ra[rn0 + i]].occ_1, @@ -13469,6 +13587,668 @@ int insert_snp_ee(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a } + +uint8_t hpc_mask_ff_adv(char *sa, int64_t sn, int64_t p, int64_t hpc_flk, int64_t hpc_rr, int64_t hpc_cutoff, + int64_t *hpc_s0, int64_t *hpc_e0, int64_t *hpc_r0, int64_t *hpc_s1, int64_t *hpc_e1, int64_t *hpc_r1) +{ + int64_t s = ((p>=hpc_flk)?(p-hpc_flk):0), e = (((p+hpc_flk+1)<=sn)?(p+hpc_flk+1):(sn)), k, r, rm, zs, ze; + int64_t max_rr[2], max_zs[2], max_ze[2], max_rl[2], max_st[2], ld[2], ll, rr; uint8_t fl, fr; + + // if(p == 3265) { + // fprintf(stderr, "+[M::%s]\tp::%ld\tsn::%ld\tp::%ld\thpc_flk::%ld\thpc_rr::%ld\thpc_cutoff::%ld\t\n", + // __func__, p, sn, p, hpc_flk, hpc_rr, hpc_cutoff); + // } + + max_rr[0] = max_rr[1] = -1; max_rl[0] = max_rl[1] = -1; max_st[0] = max_st[1] = -1; + max_zs[0] = max_zs[1] = -1; max_ze[0] = max_ze[1] = -1; + for (r = 1; r <= hpc_rr; r++) { + rm = r<<1; + + ///inlcuding p + for (k = p + r; (k < e) && ((k-r) >= s) && (sa[k] == sa[k-r]); k++){;} ze = k; if(ze > e) {ze = e;} + for (k = p - 1; (k >= s) && ((k+r) < e) && (sa[k] == sa[k+r]); k--){;} zs = k + 1; if(zs < s) {zs = s;} + if((ze - zs) >= rm) { + ld[0] = p - zs; ld[1] = ze - p - 1; ll = ze - zs; rr = ll/r; + if(ld[0] >= ld[1]) { + if((rr > max_rr[0]) || ((rr == max_rr[0]) && ((ll%r) > (max_rl[0]%max_st[0])))) { + max_rr[0] = rr; max_rl[0] = ll; + max_zs[0] = zs; max_ze[0] = ze; + max_st[0] = r; + } + } + + if(ld[0] <= ld[1]) { + if((rr > max_rr[1]) || ((rr == max_rr[1]) && (((ll%r) > (max_rl[1]%max_st[1]))))) { + max_rr[1] = rr; max_rl[1] = ll; + max_zs[1] = zs; max_ze[1] = ze; + max_st[1] = r; + } + } + } + + + ///do not inlcude p + for (k = p + r + 1; (k < e) && ((k-r) >= s) && (sa[k] == sa[k-r]); k++); + zs = p + 1; if(zs < s) zs = s; ze = k; if(ze > e) ze = e; + if((ze - zs) >= rm) { + ld[0] = p - zs; ld[1] = ze - p - 1; ll = ze - zs; rr = ll/r; + if(ld[0] >= ld[1]) { + if((rr > max_rr[0]) || ((rr == max_rr[0]) && ((ll%r) > (max_rl[0]%max_st[0])))) { + max_rr[0] = rr; max_rl[0] = ll; + max_zs[0] = zs; max_ze[0] = ze; + max_st[0] = r; + } + } + + if(ld[0] <= ld[1]) { + if((rr > max_rr[1]) || ((rr == max_rr[1]) && (((ll%r) > (max_rl[1]%max_st[1]))))) { + max_rr[1] = rr; max_rl[1] = ll; + max_zs[1] = zs; max_ze[1] = ze; + max_st[1] = r; + } + } + } + + ///inlcuding p + for (k = p - r; (k >= s) && ((k+r) < e) && (sa[k] == sa[k+r]); k--){;} zs = k + 1; if(zs < s) {zs = s;} + for (k = p + 1; (k < e) && ((k-r) >= s) && (sa[k] == sa[k-r]); k++){;} ze = k; if(ze > e) {ze = e;} + if((ze - zs) >= rm) { + ld[0] = p - zs; ld[1] = ze - p - 1; ll = ze - zs; rr = ll/r; + if(ld[0] >= ld[1]) { + if((rr > max_rr[0]) || ((rr == max_rr[0]) && ((ll%r) > (max_rl[0]%max_st[0])))) { + max_rr[0] = rr; max_rl[0] = ll; + max_zs[0] = zs; max_ze[0] = ze; + max_st[0] = r; + } + } + + if(ld[0] <= ld[1]) { + if((rr > max_rr[1]) || ((rr == max_rr[1]) && (((ll%r) > (max_rl[1]%max_st[1]))))) { + max_rr[1] = rr; max_rl[1] = ll; + max_zs[1] = zs; max_ze[1] = ze; + max_st[1] = r; + } + } + } + + ///do not inlcude p + for (k = p - r - 1; (k >= s) && ((k+r) < e) && (sa[k] == sa[k+r]); k--); + zs = k + 1; if(zs < s) zs = s; ze = p; if(ze > e) ze = e; + if((ze - zs) >= rm) { + ld[0] = p - zs; ld[1] = ze - p - 1; ll = ze - zs; rr = ll/r; + if(ld[0] >= ld[1]) { + if((rr > max_rr[0]) || ((rr == max_rr[0]) && ((ll%r) > (max_rl[0]%max_st[0])))) { + max_rr[0] = rr; max_rl[0] = ll; + max_zs[0] = zs; max_ze[0] = ze; + max_st[0] = r; + } + } + + if(ld[0] <= ld[1]) { + if((rr > max_rr[1]) || ((rr == max_rr[1]) && (((ll%r) > (max_rl[1]%max_st[1]))))) { + max_rr[1] = rr; max_rl[1] = ll; + max_zs[1] = zs; max_ze[1] = ze; + max_st[1] = r; + } + } + } + } + + // if(p == 3265) { + // fprintf(stderr, "-[M::%s]\tp::%ld\tmax_zs[0]::%ld\tmax_ze[0]::%ld\tmax_st[0]::%ld\tmax_zs[1]::%ld\tmax_ze[1]::%ld\tmax_st[1]::%ld\thpc_cutoff::%ld\n", + // __func__, p, max_zs[0], max_ze[0], max_st[0], max_zs[1], max_ze[1], max_st[1], hpc_cutoff); + // } + + fl = fr = 0; + if((max_zs[0] >= 0) && (max_ze[0] >= 0) && (max_ze[0] > max_zs[0]) && ((max_ze[0] - max_zs[0]) >= (max_st[0]*hpc_cutoff))) fl = 1; + if((max_zs[1] >= 0) && (max_ze[1] >= 0) && (max_ze[1] > max_zs[1]) && ((max_ze[1] - max_zs[1]) >= (max_st[1]*hpc_cutoff))) fr = 1; + + if(fl == 0 && fr == 0) { + *hpc_s0 = *hpc_e0 = *hpc_r0 = -1; + *hpc_s1 = *hpc_e1 = *hpc_r1 = -1; + return 0; + } + + if((max_zs[0] != max_zs[1]) || (max_ze[0] != max_ze[1])) { + if((fr) && (max_zs[0] >= max_zs[1]) && (max_ze[0] <= max_ze[1])) { + max_zs[0] = max_zs[1]; max_ze[0] = max_ze[1]; max_st[0] = max_st[1]; + } else if((fl) && (max_zs[1] >= max_zs[0]) && (max_ze[1] <= max_ze[0])) { + max_zs[1] = max_zs[0]; max_ze[1] = max_ze[0]; max_st[1] = max_st[0]; + } + } + + if(max_zs[0] < 0 || max_ze[0] < 0) { + if(max_zs[1] >= p && max_zs[1] > 0) { + max_ze[0] = max_zs[1]; max_zs[0] = max_zs[1] - 1; max_st[0] = 1; + } + } + + if(max_zs[1] < 0 || max_ze[1] < 0) { + if((max_ze[0] <= (p + 1)) && (max_ze[0] < sn)) { + max_zs[1] = max_ze[0]; max_ze[1] = max_zs[1] + 1; max_st[1] = 1; + } + } + + *hpc_s0 = max_zs[0]; *hpc_e0 = max_ze[0]; *hpc_r0 = max_st[0]; + *hpc_s1 = max_zs[1]; *hpc_e1 = max_ze[1]; *hpc_r1 = max_st[1]; + + int64_t os = MAX((*hpc_s0), (*hpc_s1)), oe = MIN((*hpc_e0), (*hpc_e1)); + if(oe < os || os < 0 || oe < 0) return 0;///not close + + // if(((*hpc_s0) == (*hpc_s1)) && ((*hpc_e0) == (*hpc_e1)) && ((*hpc_r0) == (*hpc_r1))) { + // *inner_block = 1; + // } + + fr = 0; + if((*hpc_s0) >= 0) fr |= 1; + if((*hpc_s1) >= 0) fr |= 2; + return fr; +} + + + +inline void update_Nlst(char *s, int64_t s_off, int64_t bs, int64_t be, int64_t *Nk, uint64_t *Nlst, uint8_t is_rc, int64_t tl) +{ + int64_t nk0, nn, rp; + // if((*Nk) <= 0) *Nk = Nlst[0]; + if ((!Nlst) || (Nlst[0] == 0)) return; + if ((*Nk) <= 0 || (*Nk) > (int64_t)Nlst[0]) *Nk = Nlst[0]; + + if(is_rc == 0) { + nk0 = *Nk; nn = Nlst[0]; + for (; (*Nk) <= nn; (*Nk)++) { + if ((((int64_t)Nlst[*Nk]) >= bs) && (((int64_t)Nlst[*Nk]) < be)) { + s[Nlst[*Nk] - s_off] = 'N'; + } else if(((int64_t)Nlst[*Nk]) >= be) { + break; + } + } + + for (*Nk = nk0-1; (*Nk) >= 1; (*Nk)--) { + if ((((int64_t)Nlst[*Nk]) >= bs) && (((int64_t)Nlst[*Nk]) < be)) { + s[Nlst[*Nk] - s_off] = 'N'; + } else if(((int64_t)Nlst[*Nk]) < bs) { + break; + } + } + } else { + nk0 = *Nk; nn = Nlst[0]; + for (; (*Nk) <= nn; (*Nk)++) { + rp = tl - Nlst[*Nk] - 1; + if ((rp >= bs) && (rp < be)) { + s[rp - s_off] = 'N'; + } else if(rp < bs) { + break; + } + } + + for (*Nk = nk0-1; (*Nk) >= 1; (*Nk)--) { + rp = tl - Nlst[*Nk] - 1; + if ((rp >= bs) && (rp < be)) { + s[rp - s_off] = 'N'; + } else if(rp >= be) { + break; + } + } + } +} + +///[ps, pe] but [s, e) +void inline iter_rr_match_adv(int64_t pid, char *ss, int64_t ps, int64_t pe, int64_t s, int64_t e, int64_t pl, int64_t r, int64_t *ms, int64_t *me, int64_t *Nk, uint8_t rev, uint8_t is_forward, uint8_t* src, uint64_t *Nlst, int64_t *rs, int64_t *re) +{ + assert(src); + int64_t bs, be, t2e, p2s, p2e, b2s, b2e, b2b, b2of, n_ms = INT64_MAX, n_me = INT64_MIN, os, oe; + t2e = pl - s; + if(is_forward == 1) { + if((ps <= pe) && (ps >= s) && (pe < e)) {///[bs, be) and [s, e) + if(rev == 0) { + for (; pe < e; ps++, pe++) { + bs = ((ps>>2)<<2); be = (((pe + 1) + 3)>>2)<<2; if(be > e) be = e; + if((bs < (*ms)) || (be > (*me))) { + be = bs + 4; if(be > e) be = e; + while ((bs <= pe) && (bs < be)) { + if(((bs < (*ms)) || (be > (*me))) && ((bs < n_ms) || (be > n_me))) { + memcpy(ss+bs-s, bit_t_seq_table[src[bs>>2]], be - bs); + // if(bs < (*ms)) *ms = bs; + // if(be > (*me)) *me = be; + if(bs < n_ms) n_ms = bs; + if(be > n_me) n_me = be; + if(Nlst) update_Nlst(ss, s, bs, be, Nk, Nlst, rev, pl); + } + bs = be; be += 4; if(be > e) be = e; + } + } + if(ss[ps-s] != ss[pe-s]) break; + } + } else { + for (; pe < e; ps++, pe++) { + p2s = pl - pe - 1; p2e = pl - ps - 1; + + b2s = ((p2s>>2)<<2); b2e = (((p2e + 1) + 3)>>2)<<2; if(b2e > t2e) b2e = t2e; + bs = pl - b2e; be = pl - b2s; + + if((bs < (*ms)) || (be > (*me))) { + b2s = b2e - ((b2e&3)?(b2e&3):(4));///align of 4 + be = pl - b2s; + while ((bs <= pe) && (bs < be)) { + if(((bs < (*ms)) || (be > (*me))) && ((bs < n_ms) || (be > n_me))) { + b2s = pl - be; b2e = pl - bs; + b2b = (b2s >> 2) << 2; b2of = 4 - (b2e - b2b); + memcpy(ss+bs-s, bit_t_seq_table_rc[src[b2s>>2]] + b2of, b2e - b2s); + // if(bs < (*ms)) *ms = bs; + // if(be > (*me)) *me = be; + if(bs < n_ms) n_ms = bs; + if(be > n_me) n_me = be; + if(Nlst) update_Nlst(ss, s, bs, be, Nk, Nlst, rev, pl); + } + bs = be; be += 4; if(be > e) be = e; + } + } + if(ss[ps-s] != ss[pe-s]) break; + } + } + } + // ze = pe; if(ze > e) {ze = e;} + } else { + if((ps <= pe) && (ps >= s) && (pe < e)) {///[bs, be) and [s, e) + if(rev == 0) { + for (; ps >= s; ps--, pe--) { + bs = ((ps>>2)<<2); be = (((pe + 1) + 3)>>2)<<2; if(be > e) be = e; + if((bs < (*ms)) || (be > (*me))) { + bs = be - ((be&3)?(be&3):(4));///align of 4 + while ((be > ps) && (bs < be)) { + if(((bs < (*ms)) || (be > (*me))) && ((bs < n_ms) || (be > n_me))) { + memcpy(ss+bs-s, bit_t_seq_table[src[bs>>2]], be - bs); + // if(bs < (*ms)) *ms = bs; + // if(be > (*me)) *me = be; + if(bs < n_ms) n_ms = bs; + if(be > n_me) n_me = be; + if(Nlst) update_Nlst(ss, s, bs, be, Nk, Nlst, rev, pl); + } + be = bs; bs -= 4; if(bs < s) bs = s; + } + } + + if(ss[ps-s] != ss[pe-s]) break; + } + } else { + for (; ps >= s; ps--, pe--) { + p2s = pl - pe - 1; p2e = pl - ps - 1; + b2s = ((p2s>>2)<<2); b2e = (((p2e + 1) + 3)>>2)<<2; if(b2e > t2e) b2e = t2e;///align of 4 + // b2s = ((p2s>>2)<<2); b2e = b2s + 4; if(b2e > t2e) b2e = t2e;///align of 4 + bs = pl - b2e; be = pl - b2s; + if((bs < (*ms)) || (be > (*me))) { + b2e = b2s + 4; if(b2e > t2e) b2e = t2e;///align of 4 + bs = pl - b2e; + while ((be > ps) && (bs < be)) { + if(((bs < (*ms)) || (be > (*me))) && ((bs < n_ms) || (be > n_me))) { + b2s = pl - be; b2e = pl - bs; + b2b = (b2s >> 2) << 2; b2of = 4 - (b2e - b2b); + memcpy(ss+bs-s, bit_t_seq_table_rc[src[b2s>>2]] + b2of, b2e - b2s); + // if(bs < (*ms)) *ms = bs; + // if(be > (*me)) *me = be; + if(bs < n_ms) n_ms = bs; + if(be > n_me) n_me = be; + if(Nlst) update_Nlst(ss, s, bs, be, Nk, Nlst, rev, pl); + } + be = bs; bs -= 4; if(bs < s) bs = s; + } + } + + if(ss[ps-s] != ss[pe-s]) break; + } + } + } + // zs = ps + 1; if(zs < s) {zs = s;} + } + + *rs = ps + 1; if((*rs) < s) (*rs) = s; + *re = pe; if((*re) > e) *re = e; + + if(n_ms >= n_me) return; + if((*ms) >= (*me)) { + (*ms) = n_ms; (*me) = n_me; + } else { + os = MAX(n_ms, (*ms)); oe = MIN(n_me, (*me)); + if(oe < os) { + recover_UC_Read_sub_region(ss+oe-s, oe, os-oe, rev, &R_INF, pid); + } + if(n_ms < (*ms)) (*ms) = n_ms; + if(n_me > (*me)) (*me) = n_me; + } +} + +uint8_t hpc_mask_ff_match0_test(UC_Read *su, int64_t rid, uint8_t rrev, int64_t p, UC_Read* tu, int64_t r, int64_t rr_max) +{ + int64_t k, s, e, cs, ce, zs, ze, ms = INT64_MAX, me = INT64_MIN, Nk = -1, zs1, ze1; int64_t rlen = Get_READ_LENGTH(R_INF, rid); + s = ((p>=rr_max)?(p-rr_max):0); + e = (((p+rr_max+1)<=rlen)?(p+rr_max+1):(rlen)); + if(rrev == 0) {///align of 4 forward + s = (s>>2)<<2; + e = ((e + 3)>>2)<<2; if(e > rlen) e = rlen; + } else { + cs = rlen - e; ce = rlen - s; + cs = (cs>>2)<<2; + ce = ((ce + 3)>>2)<<2; if(ce > rlen) ce = rlen; + e = rlen - cs; s = rlen - ce; + } + if(rrev == 0) { + recover_UC_Read(su, &R_INF, rid); + } else { + recover_UC_Read_RC(su, &R_INF, rid); + } + + resize_UC_Read(tu, e - s); + char *sa = su->seq; + + ///inlcuding p + for (k = p + r; (k < e) && ((k-r) >= s) && (sa[k] == sa[k-r]); k++){;} ze = k; if(ze > e) {ze = e;} + iter_rr_match_adv(rid, tu->seq, p, p+r, s, e, rlen, r, &ms, &me, &Nk, rrev, 1, Get_READ(R_INF, rid), R_INF.N_site[rid], &cs, &ce); ze1 = ce; + for (k = p - 1; (k >= s) && ((k+r) < e) && (sa[k] == sa[k+r]); k--){;} zs = k + 1; if(zs < s) {zs = s;} + iter_rr_match_adv(rid, tu->seq, p-1, p+r-1, s, e, rlen, r, &ms, &me, &Nk, rrev, 0, Get_READ(R_INF, rid), R_INF.N_site[rid], &cs, &ce); zs1 = cs; + assert(zs == zs1 && ze == ze1); + + + ///do not inlcude p + for (k = p + r + 1; (k < e) && ((k-r) >= s) && (sa[k] == sa[k-r]); k++); + zs = p + 1; if(zs < s) zs = s; ze = k; if(ze > e) ze = e; + iter_rr_match_adv(rid, tu->seq, p+1, p+r+1, s, e, rlen, r, &ms, &me, &Nk, rrev, 1, Get_READ(R_INF, rid), R_INF.N_site[rid], &cs, &ce); + zs1 = p + 1; if(zs1 < s) zs1 = s; ze1 = ce; + assert(zs == zs1 && ze == ze1); + + ///inlcuding p + for (k = p - r; (k >= s) && ((k+r) < e) && (sa[k] == sa[k+r]); k--){;} zs = k + 1; if(zs < s) {zs = s;} + iter_rr_match_adv(rid, tu->seq, p-r, p, s, e, rlen, r, &ms, &me, &Nk, rrev, 0, Get_READ(R_INF, rid), R_INF.N_site[rid], &cs, &ce); zs1 = cs; + for (k = p + 1; (k < e) && ((k-r) >= s) && (sa[k] == sa[k-r]); k++){;} ze = k; if(ze > e) {ze = e;} + iter_rr_match_adv(rid, tu->seq, p+1-r, p+1, s, e, rlen, r, &ms, &me, &Nk, rrev, 1, Get_READ(R_INF, rid), R_INF.N_site[rid], &cs, &ce); ze1 = ce; + assert(zs == zs1 && ze == ze1); + + ///do not inlcude p + for (k = p - r - 1; (k >= s) && ((k+r) < e) && (sa[k] == sa[k+r]); k--); + zs = k + 1; if(zs < s) zs = s; ze = p; if(ze > e) ze = e; + iter_rr_match_adv(rid, tu->seq, p-r-1, p, s, e, rlen, r, &ms, &me, &Nk, rrev, 0, Get_READ(R_INF, rid), R_INF.N_site[rid], &cs, &ce); + zs1 = cs; ze1 = p; if(ze1 > e) ze1 = e; + assert(zs == zs1 && ze == ze1); + + return 1; +} + +uint8_t same_repeat_pattern(const char *a, const char *b, int64_t n) +{ + int64_t shift, i; + + if (n <= 0) return 0; + + for (shift = 0; shift < n; shift++) { + for (i = 0; i < n; i++) { + if (a[i] != b[(i + shift) % n]) break; + } + + if (i == n) return 1; + } + + return 0; +} + +uint8_t hpc_mask_ff_match0(char *sp, int64_t rid, int64_t rlen, uint8_t rrev, int64_t p, char *ref, int64_t r, int64_t s, int64_t e, uint8_t is_forward, + int64_t *ms, int64_t *me, int64_t *Nk, int64_t *rs, int64_t *re) +{ + int64_t cs, ce, zs0, ze0, zs1, ze1, zl0 = 0, zl1 = 0, zs = -1, ze = -1; *rs = *re = -1; + if(!is_forward) { + ///inlcuding p + iter_rr_match_adv(rid, sp, p-r, p, s, e, rlen, r, ms, me, Nk, rrev, 0, Get_READ(R_INF, rid), R_INF.N_site[rid], &cs, &ce); zs0 = cs; + iter_rr_match_adv(rid, sp, p+1-r, p+1, s, e, rlen, r, ms, me, Nk, rrev, 1, Get_READ(R_INF, rid), R_INF.N_site[rid], &cs, &ce); ze0 = ce; + + ///do not inlcude p + iter_rr_match_adv(rid, sp, p-r-1, p-1, s, e, rlen, r, ms, me, Nk, rrev, 0, Get_READ(R_INF, rid), R_INF.N_site[rid], &cs, &ce); + zs1 = cs; ze1 = p; if(ze1 > e) ze1 = e; + } else { + ///inlcuding p + iter_rr_match_adv(rid, sp, p, p+r, s, e, rlen, r, ms, me, Nk, rrev, 1, Get_READ(R_INF, rid), R_INF.N_site[rid], &cs, &ce); ze1 = ce; + iter_rr_match_adv(rid, sp, p-1, p+r-1, s, e, rlen, r, ms, me, Nk, rrev, 0, Get_READ(R_INF, rid), R_INF.N_site[rid], &cs, &ce); zs1 = cs; + + + ///do not inlcude p + iter_rr_match_adv(rid, sp, p+1, p+r+1, s, e, rlen, r, ms, me, Nk, rrev, 1, Get_READ(R_INF, rid), R_INF.N_site[rid], &cs, &ce); + zs0 = p + 1; if(zs0 < s) zs0 = s; ze0 = ce; + } + + if(zs0 < ze0 && zs0 >= s && ze0 <= e) zl0 = ze0 - zs0; + if(zs1 < ze1 && zs1 >= s && ze1 <= e) zl1 = ze1 - zs1; + + if(zl0 >= zl1) { + zs = zs0; ze = ze0; + if((ze > zs) && (ze - zs >= r) && (same_repeat_pattern(ref, sp + zs - s, r))) { + *rs = zs; *re = ze; + return 1; + } + + zs = zs1; ze = ze1; + if((ze > zs) && (ze - zs >= r) && (same_repeat_pattern(ref, sp + zs - s, r))) { + *rs = zs; *re = ze; + return 1; + } + } else { + zs = zs1; ze = ze1; + if((ze > zs) && (ze - zs >= r) && (same_repeat_pattern(ref, sp + zs - s, r))) { + *rs = zs; *re = ze; + return 1; + } + + zs = zs0; ze = ze0; + if((ze > zs) && (ze - zs >= r) && (same_repeat_pattern(ref, sp + zs - s, r))) { + *rs = zs; *re = ze; + return 1; + } + } + + return 0; +} + + +uint8_t hpc_mask_ff_adv_0(char *sa, int64_t sn, int64_t p, int64_t hpc_flk, int64_t hpc_rr, uint8_t is_forward, int64_t *hpc_s, int64_t *hpc_e, int64_t *hpc_r) +{ + int64_t s = ((p>=hpc_flk)?(p-hpc_flk):0), e = (((p+hpc_flk+1)<=sn)?(p+hpc_flk+1):(sn)), k, r, zs, ze; + int64_t max_rr, max_zs, max_ze, max_rl, max_st, ll, rr; + *hpc_s = *hpc_e = *hpc_r = -1; max_rr = max_zs = max_ze = max_rl = max_st = -1; + + if(is_forward) { + for (r = 1; r <= hpc_rr; r++) { + ///inlcuding p + for (k = p + r; (k < e) && ((k-r) >= s) && (sa[k] == sa[k-r]); k++){;} ze = k; if(ze > e) {ze = e;} + zs = p; if(zs < s) {zs = s;} + if((ze - zs) >= r) { + ll = ze - zs; rr = ll/r; + if((max_st < 0) || (rr > max_rr) || ((rr == max_rr) && ((ll%r) > (max_rl%max_st)))) { + max_rr = rr; max_rl = ll; + max_zs = zs; max_ze = ze; + max_st = r; + } + } + } + } else { + for (r = 1; r <= hpc_rr; r++) { + ///inlcuding p + for (k = p - r; (k >= s) && ((k+r) < e) && (sa[k] == sa[k+r]); k--){;} zs = k + 1; if(zs < s) {zs = s;} + ze = p + 1; if(ze > e) {ze = e;} + if((ze - zs) >= r) { + ll = ze - zs; rr = ll/r; + if((max_st < 0) ||(rr > max_rr) || ((rr == max_rr) && ((ll%r) > (max_rl%max_st)))) { + max_rr = rr; max_rl = ll; + max_zs = zs; max_ze = ze; + max_st = r; + } + } + } + } + + if((max_zs == -1) || (max_ze == -1) || (max_st == -1)) return 0; + if((max_ze - max_zs < 1) || (max_st < 1)) return 0; + *hpc_s = max_zs; *hpc_e = max_ze; *hpc_r = max_st; + + return 1; +} + + +///needs to handle the HPC match in different directions, which is possible.... +uint8_t hpc_mask_ff_match(int64_t tid, int64_t tpos, uint8_t trev, UC_Read *tu, int64_t qpos, char *qstr, int64_t ql, int64_t rr_max, + int64_t hpc_s0, int64_t hpc_e0, int64_t hpc_r0, int64_t hpc_s1, int64_t hpc_e1, int64_t hpc_r1, int64_t sside_cut_r, int64_t bside_cut_r, int64_t sside_unio_r, uint8_t *is_strong) +{ + int64_t os, oe, roq, toq; (*is_strong) = 0; uint8_t rf = 0; + + int64_t s, e, s0, cs, ce, ms = INT64_MAX, me = INT64_MIN, Nk = -1, rts0, rte0, rts1, rte1, rqs0, rqe0, rqs1, rqe1, mqs, mqe, mts, mte; + int64_t tlen = Get_READ_LENGTH(R_INF, tid); uint8_t mm_f0 = 0, mm_f1 = 0; + s = ((tpos>=rr_max)?(tpos-rr_max):0); + e = (((tpos+rr_max+1)<=tlen)?(tpos+rr_max+1):(tlen)); + if(trev == 0) {///align of 4 forward + s = (s>>2)<<2; + e = ((e + 3)>>2)<<2; if(e > tlen) e = tlen; + } else { + cs = tlen - e; ce = tlen - s; + cs = (cs>>2)<<2; + ce = ((ce + 3)>>2)<<2; if(ce > tlen) ce = tlen; + e = tlen - cs; s = tlen - ce; + } + + rqs0 = rqe0 = rqs1 = rqe1 = -1; rts0 = rte0 = rts1 = rte1 = -1; + resize_UC_Read(tu, e - s); + if((hpc_r0 > 0) && (hpc_e0 > hpc_s0)) { + rqs0 = hpc_s0; rqe0 = hpc_e0; + mm_f0 = hpc_mask_ff_match0(tu->seq, tid, tlen, trev, tpos, qstr + rqe0 - hpc_r0, hpc_r0, s, e, 0, &ms, &me, &Nk, &rts0, &rte0); + // if((!mm_f0) && (!is_inner_block)) return 0; + } else { + return 0; + } + + if((hpc_r1 > 0) && (hpc_e1 > hpc_s1)) { + rqs1 = hpc_s1; rqe1 = hpc_e1; + mm_f1 = hpc_mask_ff_match0(tu->seq, tid, tlen, trev, tpos, qstr + rqe1 - hpc_r1, hpc_r1, s, e, 1, &ms, &me, &Nk, &rts1, &rte1); + // if((!mm_f1) && (!is_inner_block)) return 0; + } else { + return 0; + } + + if((!mm_f0) && (!mm_f1)) return 0; + + // if(/**(tid == 2198541) &&**/ (qpos == 140951)) { + // fprintf(stderr, "[M::%s]\trid::%ld(%.*s)\tqpos::%ld\ttpos::%ld\trq0::[%ld,%ld)\trt0::[%ld,%ld)\thpc_r0::%ld\tmm_f0::%u\trq1::[%ld,%ld)\trt1::[%ld,%ld)\thpc_r1::%ld\tmm_f1::%u\n", + // __func__, tid, (int)Get_NAME_LENGTH(R_INF, tid), Get_NAME(R_INF, tid), + // qpos, tpos, rqs0, rqe0, rts0, rte0, hpc_r0, mm_f0, rqs1, rqe1, rts1, rte1, hpc_r1, mm_f1); + // } + + if((mm_f0) && (mm_f1)) { + rf = 3; + } else { + ///q[rqs0, rqe0) is matched with t[rts0, rte0) + ///need to check q[rqe0, ql) and t[rte0, tl) + if(!mm_f1) { + roq = rqe0; toq = rte0; + if(roq < 0 || roq >= ql || toq < 0 || toq >= tlen) return 0; + if(roq < rqs1 || roq >= rqe1) { + if(!hpc_mask_ff_adv_0(qstr, ql, roq, 64, 4, 1, &rqs1, &rqe1, &hpc_r1)) return 0; + } + + if(toq != tpos) { + s0 = s; + s = ((toq>=rr_max)?(toq-rr_max):0); + e = (((toq+rr_max+1)<=tlen)?(toq+rr_max+1):(tlen)); + if(trev == 0) {///align of 4 forward + s = (s>>2)<<2; + e = ((e + 3)>>2)<<2; if(e > tlen) e = tlen; + } else { + cs = tlen - e; ce = tlen - s; + cs = (cs>>2)<<2; + ce = ((ce + 3)>>2)<<2; if(ce > tlen) ce = tlen; + e = tlen - cs; s = tlen - ce; + } + resize_UC_Read(tu, e - s); + if(s != s0) { + ms = INT64_MAX; + me = INT64_MIN; + Nk = -1; + } + } + + if((hpc_r1 > 0) && (rqe1 > rqs1)) { + mm_f1 = hpc_mask_ff_match0(tu->seq, tid, tlen, trev, toq, qstr + rqe1 - hpc_r1, hpc_r1, s, e, 1, &ms, &me, &Nk, &rts1, &rte1); + } + + if(!mm_f1) return 0; + rf = 1; + } + + ///q[rqs1, rqe1) is matched with t[rts1, rte1) + ///need to check q[0, rqs1) and t[0, rts1) + if(!mm_f0) { + roq = rqs1 - 1; toq = rts1 - 1; + if(roq < 0 || roq >= ql || toq < 0 || toq >= tlen) return 0; + if(roq < rqs0 || roq >= rqe0) { + if(!hpc_mask_ff_adv_0(qstr, ql, roq, 64, 4, 0, &rqs0, &rqe0, &hpc_r0)) return 0; + } + + if(toq != tpos) { + s0 = s; + s = ((toq>=rr_max)?(toq-rr_max):0); + e = (((toq+rr_max+1)<=tlen)?(toq+rr_max+1):(tlen)); + if(trev == 0) {///align of 4 forward + s = (s>>2)<<2; + e = ((e + 3)>>2)<<2; if(e > tlen) e = tlen; + } else { + cs = tlen - e; ce = tlen - s; + cs = (cs>>2)<<2; + ce = ((ce + 3)>>2)<<2; if(ce > tlen) ce = tlen; + e = tlen - cs; s = tlen - ce; + } + resize_UC_Read(tu, e - s); + if(s != s0) { + ms = INT64_MAX; + me = INT64_MIN; + Nk = -1; + } + } + + if((hpc_r0 > 0) && (rqe0 > rqs0)) { + mm_f0 = hpc_mask_ff_match0(tu->seq, tid, tlen, trev, toq, qstr + rqe0 - hpc_r0, hpc_r0, s, e, 0, &ms, &me, &Nk, &rts0, &rte0); + } + + if(!mm_f0) return 0; + rf = 2; + } + } + + + if(((rte0 - rts0) <= 1) && ((rte1 - rts1) <= 1)) return 0; + + os = MAX(rqs0, rqs1); oe = MIN(rqe0, rqe1); + if(oe < os || os < 0 || oe < 0) return 0;///not close + + os = MAX(rts0, rts1); oe = MIN(rte0, rte1); + if(oe < os || os < 0 || oe < 0) return 0;///not close + + mqs = MIN(rqs0, rqs1); mqe = MAX(rqe0, rqe1); + mts = MIN(rts0, rts1); mte = MAX(rte0, rte1); + if(qpos >= mqs && qpos < mqe && tpos >= mts && tpos < mte) { + // if(/**(tid == 2198541) &&**/ (qpos == 140951)) { + // fprintf(stderr, "[M::%s] rid::%ld(%.*s)\tqpos::%ld\tmq::[%ld,%ld)\ttpos::%ld\tmt::[%ld,%ld)\n", __func__, tid, (int)Get_NAME_LENGTH(R_INF, tid), Get_NAME(R_INF, tid), + // qpos, mqs, mqe, tpos, mts, mte); + // } + + if(((rte0 - rts0) >= (sside_cut_r*hpc_r0)) && ((rte0 - rts0) >= (sside_unio_r*hpc_r0))) { + (*is_strong) = 1; + } else if(((rte1 - rts1) >= (sside_cut_r*hpc_r1)) && ((rte1 - rts1) >= (sside_unio_r*hpc_r1))) { + (*is_strong) = 1; + } else { + int64_t tot = 0; + if((rte0 - rts0) >= (sside_unio_r*hpc_r0)) tot += (rte0 - rts0)/hpc_r0; + if((rte1 - rts1) >= (sside_unio_r*hpc_r1)) tot += (rte1 - rts1)/hpc_r1; + if(tot >= bside_cut_r) (*is_strong) = 1; + } + return rf; + } + return 0; +} + + uint8_t hpc_mask_ff(char *sa, int64_t sn, int64_t p, int64_t hpc_flk, int64_t hpc_rr, uint8_t *f, int64_t fn, int64_t fsift, int64_t hpc_cutoff, int64_t *hpc_s, int64_t *hpc_e) { int64_t s = ((p>=hpc_flk)?(p-hpc_flk):0), e = (((p+hpc_flk)<=sn)?(p+hpc_flk):(sn)), k, r, rc, zs, ze; @@ -13589,12 +14369,189 @@ uint8_t hpc_mask_ff_region(char *sa, int64_t sn, int64_t s0, int64_t e0, int64_t return 0; } +int64_t cal0_ew_ww(overlap_region *z, bit_extz_t *ez, int64_t wki, int64_t s, int64_t e, int64_t max_e) +{ + int64_t wn = z->w_list.n, wk = wki, ws, we, os, oe, ovlp; + if((wn <= 0) || (e - s <= 0)) return 0; + if(wk < 0 || wk >= wn) wk = 0; -char qck_cigar_chk(overlap_region *z, char *tstr, int64_t tl, int64_t qtarget, int64_t flk_len, int64_t *wk0, int64_t *ck0, int64_t *qk0, int64_t *tk0, uint8_t *r_op) + ws = z->w_list.a[wk].x_start; we = z->w_list.a[wk].x_end + 1; + + os = MAX(s, ws); oe = MIN(e, we); + ovlp = ((oe>os)? (oe-os):0); + if((ovlp == 0) && ((we != ws) || (ws < s) || (we >= e))) { + // if(s >= ws) { + // for (; wk < wn; wk++) { + // if((s >= z->w_list.a[wk].x_start) && (s <= z->w_list.a[wk].x_end)) break; + // } + // } else { + // for (; wk >= 0; wk--) { + // if((s >= z->w_list.a[wk].x_start) && (s <= z->w_list.a[wk].x_end)) break; + // } + // } + if(we <= s) { + for (; wk < wn; wk++) { + if(z->w_list.a[wk].x_end + 1 > s) break; + } + } else if(ws >= e) { + for (; wk >= 0; wk--) { + if(z->w_list.a[wk].x_start < e) break; + } + } + if(wk < 0 || wk >= wn) return 0; + + ws = z->w_list.a[wk].x_start; + we = z->w_list.a[wk].x_end + 1; + if(we <= s || ws >= e) return 0; + } + + + ws = z->w_list.a[wk].x_start; we = z->w_list.a[wk].x_end + 1; + os = MAX(s, ws); oe = MIN(e, we); + // ovlp = ((oe>os)? (oe-os):0); + set_bit_extz_t((*ez), (*z), wk); + + int64_t ck, qk, tk, cs = os, ce = oe, cn, tot_e = 0, ol; uint16_t op; + if(is_ualn_win(z->w_list.a[wk])) { + tot_e += z->w_list.a[wk].x_end + 1 - z->w_list.a[wk].x_start; + tot_e += z->w_list.a[wk].y_end + 1 - z->w_list.a[wk].y_start; + if(tot_e > max_e) return INT32_MAX; + } + if((os - ws) < (we - oe)) {///closer to the left end, then starting from the left + ck = 0; qk = z->w_list.a[wk].x_start; tk = z->w_list.a[wk].y_start; cn = ez->cigar.n; + + while (ck < cn && qk < e) {//[s, e) + ws = qk; + op = ez->cigar.a[ck]>>14; ol = ez->cigar.a[ck]&(0x3fff); + if(op!=2) qk += ol; + if(op!=3) tk += ol; + ck++; we = qk; + + os = MAX(s, ws); oe = MIN(e, we); + ovlp = ((oe>os)? (oe-os):0); + if((ovlp) || (op == 2 && ws >= s && we < e)) { + if(op != 0) { + if(op != 2) { + tot_e += ovlp; + } else { + tot_e += ol; + } + } + } + } + + } else { + ck = ez->cigar.n; qk = z->w_list.a[wk].x_end + 1; tk = z->w_list.a[wk].y_end + 1; + + while ((ck > 0) && (qk > s)) { + --ck; + we = qk; + op = ez->cigar.a[ck]>>14; ol = ez->cigar.a[ck]&(0x3fff); + if(op!=2) qk -= ol; + if(op!=3) tk -= ol; + ws = qk; + + os = MAX(s, ws); oe = MIN(e, we); + ovlp = ((oe>os)? (oe-os):0); + + if((ovlp) || (op == 2 && ws >= s && we < e)) { + if(op != 0) { + if(op != 2) { + tot_e += ovlp; + } else { + tot_e += ol; + } + } + } + } + } + + if(tot_e > max_e) return INT32_MAX; + + int64_t s0 = s, e0 = e, wk0 = wk;; + if(ce < e) { + s = ce; wk++; + for (; wk < wn; wk++) { + ws = z->w_list.a[wk].x_start; we = z->w_list.a[wk].x_end + 1; + if(ws >= e) break; + set_bit_extz_t((*ez), (*z), wk); + if(is_ualn_win(z->w_list.a[wk])) { + tot_e += z->w_list.a[wk].x_end + 1 - z->w_list.a[wk].x_start; + tot_e += z->w_list.a[wk].y_end + 1 - z->w_list.a[wk].y_start; + if(tot_e > max_e) return INT32_MAX; + } + ck = 0; qk = z->w_list.a[wk].x_start; tk = z->w_list.a[wk].y_start; cn = ez->cigar.n; + + while (ck < cn && qk < e) {//[s, e) + ws = qk; + op = ez->cigar.a[ck]>>14; ol = ez->cigar.a[ck]&(0x3fff); + if(op!=2) qk += ol; + if(op!=3) tk += ol; + ck++; we = qk; + + os = MAX(s, ws); oe = MIN(e, we); + ovlp = ((oe>os)? (oe-os):0); + if((ovlp) || (op == 2 && ws >= s && we < e)) { + if(op != 0) { + if(op != 2) { + tot_e += ovlp; + } else { + tot_e += ol; + } + } + } + } + } + } + if(tot_e > max_e) return INT32_MAX; + + s = s0; e = e0; wk = wk0; + if(cs > s) { + e = cs; wk--; + for (; wk >= 0; wk--) { + ws = z->w_list.a[wk].x_start; we = z->w_list.a[wk].x_end + 1; + if(we <= s) break; + set_bit_extz_t((*ez), (*z), wk); + if(is_ualn_win(z->w_list.a[wk])) { + tot_e += z->w_list.a[wk].x_end + 1 - z->w_list.a[wk].x_start; + tot_e += z->w_list.a[wk].y_end + 1 - z->w_list.a[wk].y_start; + if(tot_e > max_e) return INT32_MAX; + } + ck = ez->cigar.n; qk = z->w_list.a[wk].x_end + 1; tk = z->w_list.a[wk].y_end + 1; + + while ((ck > 0) && (qk > s)) { + --ck; + we = qk; + op = ez->cigar.a[ck]>>14; ol = ez->cigar.a[ck]&(0x3fff); + if(op!=2) qk -= ol; + if(op!=3) tk -= ol; + ws = qk; + + os = MAX(s, ws); oe = MIN(e, we); + ovlp = ((oe>os)? (oe-os):0); + + if((ovlp) || (op == 2 && ws >= s && we < e)) { + if(op != 0) { + if(op != 2) { + tot_e += ovlp; + } else { + tot_e += ol; + } + } + } + } + } + } + + if(tot_e > max_e) return INT32_MAX; + return tot_e; +} + +char qck_cigar_chk(overlap_region *z, char *tstr, int64_t tl, int64_t qtarget, int64_t flk_len, int64_t flk_err_add, int64_t *wk0, int64_t *ck0, int64_t *qk0, int64_t *tk0, uint8_t *r_op) { (*r_op) = ((uint8_t)-1); if((qtarget < z->x_pos_s) || (qtarget > z->x_pos_e)) return ((char)-1); - int64_t wk = *wk0, ck = *ck0, qk = *qk0, tk = *tk0, wn = z->w_list.n, cn = -1; uint16_t op; + int64_t wk = *wk0, ck = *ck0, qk = *qk0, tk = *tk0, wn = z->w_list.n, cn = -1; char rc = ((char)-1); uint16_t op; if(wn <= 0) return ((char)-1); if((wk >= 0) && (wk < wn) && (qtarget >= z->w_list.a[wk].x_start) && (qtarget <= z->w_list.a[wk].x_end)) { @@ -13636,29 +14593,40 @@ char qck_cigar_chk(overlap_region *z, char *tstr, int64_t tl, int64_t qtarget, i bit_extz_t ez; set_bit_extz_t(ez, (*z), wk); if(!ez.cigar.n) return ((char)-1); - int64_t ws, we, ttarget = -1; cn = ez.cigar.n; + int64_t ws, we, ttarget = -1, ol, qs = qtarget - flk_len, qe = qtarget + flk_len; cn = ez.cigar.n; if((ck < 0) || (ck > cn)) {//(*ck) == cn is allowed ck = 0; qk = ez.ts; tk = ez.ps; } - while ((ck > 0) && (qk > qtarget)) {///x -> t; y -> p + while ((ck > 0) && (qk > qs)) {///x -> t; y -> p --ck; op = ez.cigar.a[ck]>>14; if(op!=2) qk -= (ez.cigar.a[ck]&(0x3fff)); if(op!=3) tk -= (ez.cigar.a[ck]&(0x3fff)); } - // int64_t qs = qtarget - flk_len, qe = qtarget + flk_len, os, oe; + int64_t os, oe, ovlp, min_q = INT64_MAX, max_q = INT64_MIN, tot_e = 0; while (ck < cn && qk <= qtarget) {//[s, e) ws = qk; - op = ez.cigar.a[ck]>>14; - if(op!=2) qk += (ez.cigar.a[ck]&(0x3fff)); - if(op!=3) tk += (ez.cigar.a[ck]&(0x3fff)); + op = ez.cigar.a[ck]>>14; ol = ez.cigar.a[ck]&(0x3fff); + if(op!=2) qk += ol; + if(op!=3) tk += ol; ck++; we = qk; - // os = MAX(qs, ws); oe = MIN(qe, we); - // ovlp = ((oe>os)? (oe-os):0); + os = MAX(qs, ws); oe = MIN(qe, we); + ovlp = ((oe>os)? (oe-os):0); + if((ovlp) || (op == 2 && ws >= qs && we < qe)) { + if(os < min_q) min_q = os; + if(oe > max_q) max_q = oe; + if(op != 0) { + if(op != 2) { + tot_e += ovlp; + } else { + tot_e += ol; + } + } + } if(op != 0 && op != 1) continue;///only collect match/snp @@ -13666,33 +14634,201 @@ char qck_cigar_chk(overlap_region *z, char *tstr, int64_t tl, int64_t qtarget, i (*r_op) = op; ttarget = qtarget-qk+tk; (*wk0) = wk; (*ck0) = ck; (*qk0) = qk; (*tk0) = tk; if((op == 1) && (tl > ttarget) && (ttarget >= 0)) { - - - - - - - return tstr[ttarget]; + if((tot_e) > (flk_err_add + 1)) return ((char)-1); + rc = tstr[ttarget]; + break; } return ((char)-1); } } (*wk0) = wk; (*ck0) = ck; (*qk0) = qk; (*tk0) = tk; + if(rc != ((char)-1)) { + while (ck < cn && qk < qe) {//[s, e) + ws = qk; + op = ez.cigar.a[ck]>>14; ol = ez.cigar.a[ck]&(0x3fff); + if(op!=2) qk += ol; + if(op!=3) tk += ol; + ck++; we = qk; + + os = MAX(qs, ws); oe = MIN(qe, we); + ovlp = ((oe>os)? (oe-os):0); + if((ovlp) || (op == 2 && ws >= qs && we < qe)) { + if(os < min_q) min_q = os; + if(oe > max_q) max_q = oe; + if(op != 0) { + if(op != 2) { + tot_e += ovlp; + } else { + tot_e += ol; + } + } + } + } + + if((tot_e) > (flk_err_add + 1)) return ((char)-1); + if((min_q > qs) && (wk > 0)) { + tot_e += cal0_ew_ww(z, &ez, wk - 1, qs, min_q, flk_err_add + 1 - tot_e); + if((tot_e) > (flk_err_add + 1)) return ((char)-1); + } + if((max_q < qe) && ((wk + 1) < wn)) { + tot_e += cal0_ew_ww(z, &ez, wk + 1, max_q, qe, flk_err_add + 1 - tot_e); + if((tot_e) > (flk_err_add + 1)) return ((char)-1); + } + + return rc; + } return ((char)-1); } -int push_info_ssl(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, haplotype_evdience* u_a, overlap_region *oa, asg8_v *v8, uint64_t tot_cov, uint64_t scw, uint64_t tcut, - char r1a, overlap_region *rchn, char *tstr, int64_t tlen, int64_t flk_len, int64_t *wk, int64_t *ck, int64_t *qk, int64_t *tk) +inline uint8_t update_hpc_rr(char *qstr, int64_t ql, haplotype_evdience* a, int64_t a_n, overlap_region *oa, uint64_t oid_mm, UC_Read *tu, char r0a, int64_t hpc_cut, uint64_t tcut, + uint64_t *r_occ_0, uint64_t r_occ_1[], uint64_t r_occ_1_tqn[], uint64_t *r_occ_0_tqn, uint64_t *r_occ_2, uint64_t *r_diff, uint64_t *r_rev_n, uint8_t *r_ihpc, + double site_hpc_rate, double site_hpc_strong_rate, uint8_t *is_site_hpc) +{ + int64_t hpc_s0, hpc_e0, hpc_r0, hpc_s1, hpc_e1, hpc_r1, k; uint32_t qsite = UINT32_MAX, qsm = (((uint32_t)1)<<29)-1, qstrong = (((uint32_t)1)<<29); + if(a_n) qsite = a[0].site; + (*is_site_hpc) = 0; + + if(hpc_mask_ff_adv(qstr, ql, qsite, 64, 4, 3, &hpc_s0, &hpc_e0, &hpc_r0, &hpc_s1, &hpc_e1, &hpc_r1)) { + uint8_t is, ifl, rst[2] = {0, 0}, rcle = 0, ihp, sf; int64_t nfh[2] = {1, 1}, nss[2] = {1, 1}; uint32_t ft, fh; ///qstr itself + for (k = 0; k < a_n; k++) { + fh = hpc_mask_ff_match(oa[a[k].overlapID&oid_mm].y_id, a[k].overlapSite, oa[a[k].overlapID&oid_mm].y_pos_strand, tu, qsite, qstr, ql, 64, + hpc_s0, hpc_e0, hpc_r0, hpc_s1, hpc_e1, hpc_r1, 4, 7, 3, &is); + if(fh) { + if(hh_tp(a[k])) { + a[k].site = qsm; + if(is) a[k].site |= qstrong; + a[k].site |= (fh<<30); + rcle = 1; + + // fprintf(stderr, "sa[M::%s]\trid::%u(%.*s)\tqpos::%u\ttpos::%u\tis_strong::%u\tis_forward::%u\tis_backward::%u, hpc_s0::%ld, hpc_e0::%ld, hpc_r0::%ld, hpc_s1::%ld, hpc_e1::%ld, hpc_r1::%ld\n", + // __func__, oa[a[k].overlapID&oid_mm].y_id, + // (int)Get_NAME_LENGTH(R_INF, oa[a[k].overlapID&oid_mm].y_id), Get_NAME(R_INF, oa[a[k].overlapID&oid_mm].y_id), + // qsite, a[k].overlapSite, is, fh&1, fh&2, hpc_s0, hpc_e0, hpc_r0, hpc_s1, hpc_e1, hpc_r1); + } + if(fh&1) { + nfh[0]++; if(is) nss[0]++; + } + if(fh&2) { + nfh[1]++; if(is) nss[1]++; + } + } + } + // fprintf(stderr, "+[M::%s::site->%u]\tnss[0]::%ld\tnss[1]::%ld\tnfh[0]::%ld\tnfh[1]::%ld\n", __func__, qsite, nss[0], nss[1], nfh[0], nfh[1]); + + rst[0] = 0; + if((nss[0] > 0) && (nss[0] >= (nfh[0]*0.8))) { + nss[0] = nfh[0]; rst[0] = 2; + } + if((nss[0] > 0) && (nss[0] >= hpc_cut)) { + rst[0] |= 1; + if((nss[0] >= ((a_n+1)*site_hpc_strong_rate)) && (nfh[0] >= ((a_n+1)*site_hpc_rate))) { + (*is_site_hpc) = 1; + } + } + + rst[1] = 0; + if((nss[1] > 0) && (nss[1] >= (nfh[1]*0.8))) { + nss[1] = nfh[1]; rst[1] = 2; + } + if((nss[1] > 0) && (nss[1] >= hpc_cut)) { + rst[1] |= 1; + if((nss[1] >= ((a_n+1)*site_hpc_strong_rate)) && (nfh[1] >= ((a_n+1)*site_hpc_rate))) { + (*is_site_hpc) = 1; + } + } + // fprintf(stderr, "-[M::%s::site->%u]\tnss[0]::%ld\tnss[1]::%ld\tnfh[0]::%ld\tnfh[1]::%ld\trst_%u\n", __func__, qsite, nss[0], nss[1], nfh[0], nfh[1], ((rst[0])|(rst[1]))?1:0); + + if(rcle) { + (*r_occ_0) = (*r_occ_0_tqn) = (*r_occ_2) = (*r_diff) = (*r_ihpc) = 0; if(r_rev_n) (*r_rev_n) = 0; + memset(r_occ_1, 0, sizeof(uint64_t)*6); memset(r_occ_1_tqn, 0, sizeof(uint64_t)*6); + for (k = 0; k < a_n; k++) { + ihp = is = 0; + ihp = (a[k].site>=qsm)?(1):(0); + if(ihp) { + is = (a[k].site&qstrong)?(1):(0); sf = a[k].site>>30; ifl = 0; + + if((rst[0]&1) && (sf&1) && (!ifl)) { + if((is) || (rst[0]&2)) { + a[k].type ^= ((uint32_t)1); a[k].misBase = r0a; ifl = 1; + // fprintf(stderr, "sb[M::%s]\trid::%u(%.*s)\tqpos::%u\ttpos::%u\n", __func__, oa[a[k].overlapID&oid_mm].y_id, + // (int)Get_NAME_LENGTH(R_INF, oa[a[k].overlapID&oid_mm].y_id), Get_NAME(R_INF, oa[a[k].overlapID&oid_mm].y_id), + // qsite, a[k].overlapSite); + } + } + if((rst[1]&1) && (sf&2) && (!ifl)) { + if((is) || (rst[1]&2)) { + a[k].type ^= ((uint32_t)1); a[k].misBase = r0a; ifl = 1; + // fprintf(stderr, "sb[M::%s]\trid::%u(%.*s)\tqpos::%u\ttpos::%u\n", __func__, oa[a[k].overlapID&oid_mm].y_id, + // (int)Get_NAME_LENGTH(R_INF, oa[a[k].overlapID&oid_mm].y_id), Get_NAME(R_INF, oa[a[k].overlapID&oid_mm].y_id), + // qsite, a[k].overlapSite); + } + } + } + + a[k].site = qsite; + + if(tcut == ((uint64_t)-1)) { + if(r0a == a[k].misBase) { + (*r_occ_0) += a[k].cov; + if((r_rev_n) && (oa) && (oa[a[k].overlapID&oid_mm].y_pos_strand == 0)) { + (*r_rev_n) += a[k].cov; + } + if((!(*r_ihpc)) && (hh_hp(a[k]))) (*r_ihpc) = 1; + ft = 0; + } else { + r_occ_1[seq_nt6_table[(uint8_t)(a[k].misBase)]] += a[k].cov; + (*r_diff) += a[k].cov; + ft = 1; + } + if(hh_tp(a[k]) != ft) a[k].type ^= ((uint32_t)1); + + (*r_occ_2) += a[k].cov; + } else if(oa) { + if(r0a == a[k].misBase){ + (*r_occ_0) += a[k].cov; + if((r_rev_n) && (oa[a[k].overlapID&oid_mm].y_pos_strand == 0)) { + (*r_rev_n) += a[k].cov; + } + if((!(*r_ihpc)) && (hh_hp(a[k]))) (*r_ihpc) = 1; + if(oa[a[k].overlapID&oid_mm].y_id < tcut) { + (*r_occ_0_tqn) += a[k].cov; + } + ft = 0; + } else { + r_occ_1[seq_nt6_table[(uint8_t)(a[k].misBase)]] += a[k].cov; + if(oa[a[k].overlapID&oid_mm].y_id < tcut) { + r_occ_1_tqn[seq_nt6_table[(uint8_t)(a[k].misBase)]] += a[k].cov; + } + (*r_diff) += a[k].cov; + ft = 1; + } + if(hh_tp(a[k]) != ft) a[k].type ^= ((uint32_t)1); + + (*r_occ_2) += a[k].cov; + } + } + } + } + + // fprintf(stderr, "[M::%s::site->%u]\tr_occ_0::%lu\tr_diff::%lu\n", __func__, a[0].site, *r_occ_0, *r_diff); + if((*r_occ_0) == 0 || (*r_diff) <= 1) return 0; + if((r_rev_n) && ((*r_rev_n) == (*r_occ_0))) return 0; + return 1; +} + +int push_info_ssl(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, haplotype_evdience* u_a, overlap_region *oa, asg8_v *v8, uint64_t tot_cov, uint64_t scw, uint64_t tcut, uint64_t rid, + char *qstr, uint64_t ql, char r1a, overlap_region *rchn, char *tstr, int64_t tlen, int64_t flk_len, int64_t *wk, int64_t *ck, int64_t *qk, int64_t *tk, UC_Read *tu, uint64_t hpc_cut, overlap_region *ohp_a, + int64_t min_re_cut, double min_re_rt, int64_t min_hpc_re_cut, double min_hpc_re_rt, uint64_t hom_cov_a) { if(a_n < 2) return 0; - uint64_t i, hi = 0, k, m, occ_0, occ_1[6], occ_2, diff, rev_n; char r0a = r1a; uint8_t ihpc = 0; const uint32_t HQ_MASK = 1u << 31; const uint32_t OD_MASK = ~HQ_MASK; + uint64_t i = 0, hi = 0, k, m, occ_0, occ_1[6], occ_1_tqn[6], occ_0_tqn, occ_2, diff, rev_n; char r0a = r1a; uint8_t ihpc = 0, site_hpc = 0; const uint32_t HQ_MASK = 1u << 31; const uint32_t OD_MASK = ~HQ_MASK; haplotype_evdience at; haplotype_evdience *ra = NULL; uint8_t op; uint32_t ft; - occ_0 = occ_2 = diff = rev_n = 0; memset(occ_1, 0, sizeof(uint64_t)*6); + occ_0 = occ_0_tqn = occ_2 = diff = rev_n = 0; memset(occ_1, 0, sizeof(uint64_t)*6); memset(occ_1_tqn, 0, sizeof(uint64_t)*6); radix_sort_haplotype_evdience_id_srt(a, a + a_n); if((rchn) && (tstr) && (tlen > 0)) { - r0a = qck_cigar_chk(rchn, tstr, tlen, a[0].site, flk_len, wk, ck, qk, tk, &op); + r0a = qck_cigar_chk(rchn, tstr, tlen, a[0].site, flk_len, h0_w_p, wk, ck, qk, tk, &op); if((op != 1) || (r0a == ((char)-1))) r0a = r1a; if(r0a != r1a) { for (k = 0; k < a_n; k++) { @@ -13734,9 +14870,15 @@ int push_info_ssl(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a occ_0 += a[i].cov; if(oa[a[i].overlapID].y_pos_strand == 0) rev_n += a[i].cov; if((!ihpc) && (hh_hp(a[i]))) ihpc = 1; + if(oa[a[i].overlapID].y_id < tcut) { + occ_0_tqn += a[i].cov; + } ft = 0; } else { occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]] += a[i].cov; + if(oa[a[i].overlapID].y_id < tcut) { + occ_1_tqn[seq_nt6_table[(uint8_t)(a[i].misBase)]] += a[i].cov; + } diff += a[i].cov; ft = 1; } @@ -13781,9 +14923,46 @@ int push_info_ssl(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a 3. if occ_1 = 1, there are only one difference. It must be a sequencing error. (for repeat, it maybe a snp at repeat. but ...) **/ +// if(a[0].site == 50637) { +// fprintf(stderr, "[M::%s::site->%u]\tocc_0::%lu\tocc_0_tqn::%lu\trev_n::%lu\tdiff::%lu\tocc_1[0]::%lu\tocc_1[1]::%lu\tocc_1[2]::%lu\tocc_1[3]::%lu\t\n", +// __func__, a[0].site, occ_0 + 1, occ_0_tqn + ((rid%u]\tocc_0::%lu\tocc_1[A]::%lu\tocc_1[C]::%lu\tocc_1[G]::%lu\tocc_1[T]::%lu\n", __func__, a[0].site, + // occ_0 + 1, occ_1[0], occ_1[1], occ_1[2], occ_1[3]); + return 0; + } + } + if(site_hpc) { + ccut = min_hpc_re_cut; + if((min_hpc_re_rt >= 0) && (ccut < (mcut*min_hpc_re_rt))) { + ccut = mcut*min_hpc_re_rt; + } + } else { + ccut = min_re_cut; + if((min_re_rt >= 0) && (ccut < (mcut*min_re_rt))) { + ccut = mcut*min_re_rt; + } + } + if(ccut < 0) ccut = 0; + } + + if((ccut >= 0) && ((occ_0 + 1) < ((uint64_t)ccut))) { + // fprintf(stderr, "-0-[M::%s::site->%u]\tocc_0::%lu\tocc_1[A]::%lu\tocc_1[C]::%lu\tocc_1[G]::%lu\tocc_1[T]::%lu\tccut::%lu\n", __func__, a[0].site, + // occ_0 + 1, occ_1[0], occ_1[1], occ_1[2], occ_1[3], ccut); + return 0; + } + if(!oa) { rev_n = occ_2; } else { @@ -13791,7 +14970,15 @@ int push_info_ssl(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a } for (i = m = 0; i < 4; i++) { - if(occ_1[i] >= 2){ + // if((occ_1[i] >= 2) && (occ_1[i] < ccut)){ + // fprintf(stderr, "-1-[M::%s::site->%u]\tocc_0::%lu\tocc_1[i]::%lu\tccut::%lu\n", __func__, a[0].site, + // occ_0 + 1, occ_1[i], ccut); + // } + + if((occ_1[i] >= 2) && (occ_1[i] >= ((uint64_t)(ccut)))){ + // if(tcut != ((uint64_t)-1)) { + // fprintf(stderr, "[M::%s::site->%u]\tocc_1::%lu\tocc_1_tqn::%lu\tocc_0::%lu\tocc_0_tqn::%lu\n", __func__, a[0].site, occ_1[i], occ_1_tqn[i], occ_0 + 1, occ_0_tqn + ((ridsnp_stat, &p); p->id = h->snp_stat.n-1; p->occ_0 = 1 + occ_0; @@ -14066,7 +15253,7 @@ void partition_overlaps_advance(overlap_region_alloc* overlap_list, All_reads* R hap->length = m; // generate_haplotypes_naive_advance(hap, overlap_list, NULL); - generate_haplotypes_naive_HiFi(hap, overlap_list, 0.04, 1, 0, ((uint64_t)-1), asm_opt.s_hap_cov, asm_opt.infor_cov); + generate_haplotypes_naive_HiFi(hap, overlap_list, 0.04, 1, 0, ((uint64_t)-1), asm_opt.s_hap_cov, asm_opt.infor_cov, INT64_MAX); // generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat); // generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat); @@ -16239,6 +17426,11 @@ uint32_t push_hc_wlst_exz(const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, o // if(ol->x_id == 1 && ol->y_id == 24) { // fprintf(stderr, "-a-[M::%s]\tqid::%u\ttid::%uovl::%ld\taln::%ld\tovlp_cut::%f\n", __func__, ol->x_id, ol->y_id, ovl, aln, ovlp_cut); // } + // if(ol->y_id == 4192664 || ol->y_id == 1184850) { + // fprintf(stderr, "[M::%s::]\ttid::%u(%u)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\tovl::%ld\tcur_aln::%u\tun_aln::%ld\tinfer_aln::%ld\tovlp_cut::%f\n", __func__, + // ol->y_id, ol->y_pos_strand, (int)Get_NAME_LENGTH(R_INF, ol->y_id), Get_NAME(R_INF, ol->y_id), ol->x_pos_s, ol->x_pos_e + 1, ol->y_pos_s, ol->y_pos_e + 1, + // ovl, ol->align_length, ualn, aln, ovlp_cut); + // } if(!is_srt) { kv_push(window_list, ol->w_list, p); return 0; @@ -26234,7 +27426,8 @@ void est_rep_err_rate(overlap_region_alloc* ol, asg64_v *ix, kv_ul_ov_t *c_idx, } void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, asg16_v *qc0, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, Chain_Data *dp, asg8_v *q8, asg8_v *t8, uint8_t lindel, uint64_t tcut, uint64_t site_sc, int64_t h0_w, asg32_v *b32, - int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate, overlap_region *rchn, int64_t rref_len) + int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate, overlap_region *rchn, int64_t rref_len, int64_t flag_hf_ov_cut, + int64_t min_re_cut, double min_re_rt, int64_t min_hpc_re_cut, double min_hpc_re_rt, int8_t is_hpc_flt) { int64_t on = ol->length, k, i, zwn, q[2]; if(occ_thres + 1 < hap_cov_unmatch) occ_thres = hap_cov_unmatch-1; uint64_t m, l0, wi, wl0, si, ei, fi; overlap_region *z; ul_ov_t *cp; uint8_t *qhf = NULL, *qual = NULL; @@ -26380,7 +27573,7 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all } **/ - int64_t wk = -1, ck = -1, qk = -1, tk = -1; + int64_t wk = -1, ck = -1, qk = -1, tk = -1; overlap_region *zre = ((is_hpc_flt)?(ol->list):(NULL)); SetSnpMatrix(hp, &(hp->nn_snp), &(ol->length), 0, NULL); srt_n = hp->length; z = (((std_bs)||(t8))?(ol->list):(NULL)); for (k = 1, m = /**ms =**/ i = t = 0; k <= srt_n; ++k) { @@ -26396,8 +27589,9 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all // if(m >= idx->n) then the region will be [ms, ql), odp = 0; // dbg_cp_cov(ol, c_idx, hp->list[i].site, odp); // t += push_info(hp, hp->list+i, k-i, hp->list+t, z, t8, 0/**odp**/, sc_wn, tcut); - t += push_info_ssl(hp, hp->list+i, k-i, hp->list+t, z, t8, 0/**odp**/, sc_wn, tcut, qu->seq[hp->list[i].site], - rchn, qu->seq + qu->length, rref_len, &wk, &ck, &qk, &tk); + t += push_info_ssl(hp, hp->list+i, k-i, hp->list+t, z, t8, 0/**odp**/, sc_wn, tcut, rid, qu->seq, qu->length, qu->seq[hp->list[i].site], + rchn, qu->seq + qu->length, rref_len, MIN(5, h0_w), &wk, &ck, &qk, &tk, tu, HPC_COMP_LEN, zre, + min_re_cut, min_re_rt, min_hpc_re_cut, min_hpc_re_rt, hom_cov_a); i = k; } } @@ -26412,7 +27606,7 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all if(!dp) { // generate_haplotypes_naive_advance(hap, overlap_list, NULL); - generate_haplotypes_naive_HiFi(hp, ol, 0.04, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), hap_cov_match, hap_cov_unmatch); + generate_haplotypes_naive_HiFi(hp, ol, 0.04, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), hap_cov_match, hap_cov_unmatch, flag_hf_ov_cut); // generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat); // generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat); } else { @@ -26421,13 +27615,13 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all // return; // if(!site_sc) generate_haplotypes_naive_HiFi_adv(hp, ol, 0.04, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), b32, rid);//r835 - if(!site_sc) generate_haplotypes_naive_HiFi_adv_hc(hp, ol, 0.04, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), b32, rid, hap_cov_match, hap_cov_unmatch);///r835 - else generate_haplotypes_weight(hp, ol, 0.04, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), 32, 2, hap_cov_match, hap_cov_unmatch); + if(!site_sc) generate_haplotypes_naive_HiFi_adv_hc(hp, ol, 0.04, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), b32, rid, hap_cov_match, hap_cov_unmatch, flag_hf_ov_cut);///r835 + else generate_haplotypes_weight(hp, ol, 0.04, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), 32, 2, hap_cov_match, hap_cov_unmatch, flag_hf_ov_cut); } if(lindel) { if(rphase_lidel(ol, rref, hp, qu, tu, c_idx, idx, buf, bd, wl, ql, occ_thres, rid, hpc_len, std_bs, hom_cov_a, het_cov_a, n_hap)) { - generate_haplotypes_sv(hp, ol, rid, hap_cov_match, hap_cov_unmatch); + generate_haplotypes_sv(hp, ol, rid, hap_cov_match, hap_cov_unmatch, (flag_hf_ov_cut>=INT64_MAX)?(0):(1)); } } } @@ -37138,6 +38332,209 @@ double fcov_rat, uint64_t ch_occ, uint64_t ch_sc) return wsrt; } +#define is_exact_matched_ov(a) ((a).is_match == 1) +#define is_hwh_matched_ov(a) ((a).is_match == 2) + +uint64_t* mmp_chn_select_adv(overlap_region_alloc* ol, Candidates_list *cl, asg64_v *sp, uint64_t ocw, uint32_t *ocn, uint32_t *osc, uint64_t ql, uint64_t *rwsrt_n, uint8_t set_match, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, +double fcov_rat, uint64_t ch_occ, uint64_t ch_sc) +{ + (*rwsrt_n) = 0; + if(ol->length <= 0) return NULL; + uint8_t fc = ((ol->length > max_n_chain)?(1):(0)), ff, wf, lch = 0; overlap_region *z; + uint64_t spn0 = sp->n, *wsrt = NULL, *wcut = NULL, focv = 0, fcov_0 = 0, i, k, m, wcut_n = 0, wsrt_n = 0, wma_n = 0, bs, tot_b = 0, t_cov0, z_cov, tz; int32_t s[4]; uint32_t scn[4]; + scn[0] = scn[1] = scn[2] = scn[3] = 0; s[0] = s[1] = s[2] = s[3] = 0; + sp->n += ol->length; kv_resize(uint64_t, (*sp), sp->n); wsrt = sp->a + spn0; + + // fprintf(stderr, "[M::%s::]\tol->length::%lu\tmax_n_chain::%lu\n", __func__, ol->length, max_n_chain); + if(fc) { + permute_in_place_ss(ol->list, ol->length, wsrt, osc, ocn); + wcut = infer_chn_bar_0(ol, max_n_chain, ocn, osc, ql, ocw, s, sp, &wcut_n, &focv); + if ((s[0] <= 0) && (s[1] <= 0) && (s[2] <= 0) && (s[3] <= 0)) fc = 0; + wsrt = sp->a + spn0; + + focv = focv/ql; if(focv < 0) focv = 1; focv *= 1.05; + if(focv < ave_cov_min) focv = ave_cov_min; + if(focv > (max_n_chain<<1)) focv = max_n_chain<<1; + // if(focv > ave_cov_max) focv = ave_cov_max; + // fcov_0 = max_n_chain>>1; if(fcov_0 > focv) fcov_0 = focv; + fcov_0 = max_n_chain*fcov_rat; + } + + + s[0] = s[1] = s[2] = s[3] = 0;///reset it + for (i = wsrt_n = wma_n = 0; i < ol->length; i++) {///primary chain + z = &(ol->list[i]); + + if(set_match) z->is_match = 0; + + if(!is_exact_matched_ov((*z))) { + z->shared_seed = z->non_homopolymer_errors; + z->non_homopolymer_errors = UINT32_MAX;///for index + } + + wf = ha_ov_type(z, ql); ff = 1; + + if(fc) { + /** + if(((scn[wf] <= max_n_chain) || ((scn[wf] <= (max_n_chain<<1)) && (((int64_t)osc[i]) == s[wf])) || (z->is_match == 1))) { + ff = update_mm_wins(z, wcut, wcut_n, ocw, ql, ((scn[wf]>max_n_chain_f) && (z->is_match == 0))?1:0, 0.15, 16, focv, 0); + } else { + ff = 0; + } + **/ + if((scn[wf] <= max_n_chain) || ((scn[wf] <= (max_n_chain<<1)) && (((int64_t)osc[i]) == s[wf])) || (is_exact_matched_ov((*z))) || (is_hwh_matched_ov((*z)))) { + ///change ((scn[wf]>max_n_chain_f) && (!is_exact_matched_ov((*z)) && (is_hwh_matched_ov((*z))))?1:0 to (scn[wf]>max_n_chain_f)?1:0 + ///because we don't want to is_exact_matched_ov((*z) or is_hwh_matched_ov((*z)) consume the space + ///this will let more overlaps to be test ((!is_exact_matched_ov((*z)) && (is_hwh_matched_ov((*z)))) + ff = update_mm_wins(z, wcut, wcut_n, ocw, ql, (scn[wf]>max_n_chain_f)?1:0, 0.15, 16, focv, 0); + if((is_exact_matched_ov((*z))) || (is_hwh_matched_ov((*z)))) ff = 1; + } else { + ff = 0; + } + } + + if((!ff) && (!is_exact_matched_ov((*z))) && (!is_hwh_matched_ov((*z)))) {///fitered out due to coverage + wsrt[wsrt_n++] = (((uint64_t)osc[i])<<32)|(i); + continue; + } + + scn[wf]++; + if (scn[wf] == max_n_chain) s[wf] = osc[i]; + + if((ocn[i] < chain_cutoff) && (!is_exact_matched_ov((*z))) && (!is_hwh_matched_ov((*z)))) {///fitered out due to no enough minimizers + lch = 1; wsrt[wsrt_n++] = (((uint64_t)-1)<<32)|(i); + continue; + } + + wsrt[wsrt_n] = (((uint64_t)z->x_pos_s)<<32)|(i); + if(wsrt_n != wma_n) { + bs = wsrt[wsrt_n]; wsrt[wsrt_n] = wsrt[wma_n]; wsrt[wma_n] = bs; + } + + wsrt_n++; wma_n++; tot_b += z->x_pos_e + 1 - z->x_pos_s; + z->non_homopolymer_errors = UINT32_MAX - 1;///primary chain that needs to be verfied + if(is_hwh_matched_ov((*z))) z->is_match = 0; + } + + if(fc && wsrt_n > wma_n) { + t_cov0 = cal_mm_wins_cov(wcut, wcut_n, ocw, ql, fcov_0); + + for (tz = m = wma_n, z_cov = 0, lch = 0; (tz < wsrt_n) && (z_cov <= t_cov0); tz++) { + i = (uint32_t)wsrt[tz]; + + z = &(ol->list[i]); + assert(z->non_homopolymer_errors == UINT32_MAX); + assert((!is_exact_matched_ov((*z))) && (!is_hwh_matched_ov((*z)))); + + if((wsrt[tz]>>32) == ((uint32_t)-1)) {///fitered out due to no enough minimizers + lch = 1; continue; + } + + if(update_mm_wins(z, wcut, wcut_n, ocw, ql, 1, 0.15, 16, fcov_0, 0) == 0) { + wsrt[tz] = (((uint64_t)osc[i])<<32)|(i); continue; + } + + if(ocn[i] < chain_cutoff) {///fitered out due to no enough minimizers + lch = 1; wsrt[tz] = (((uint64_t)-1)<<32)|(i); continue; + } + + wsrt[tz] = (((uint64_t)z->x_pos_s)<<32)|(i); + if(m != tz) { + bs = wsrt[tz]; wsrt[tz] = wsrt[m]; wsrt[m] = bs; + } + m++; z_cov += z->x_pos_e + 1 - z->x_pos_s; + z->non_homopolymer_errors = UINT32_MAX - 1;///primary chain that needs to be verfied + } + wma_n = m; tot_b += z_cov; + } + + if(lch) { + assert(wsrt_n > wma_n); + uint64_t *sb = wsrt, sb_n = wma_n, *sa = wsrt + wma_n, sa_n = 0, ncut = ch_occ; ///ncut = 1*ch_occ as mini_chain_occ is 1 + for (i = tz = 0; i < sb_n; i++) { + if(ocn[(uint32_t)sb[i]] < ncut) continue; + if(tz < i) { + bs = sb[i]; sb[i] = sb[tz]; sb[tz] = bs; + } + tz++; + } + sb_n = tz; + + if(sb_n > 0) { + sa_n = wsrt_n - wma_n; + for (i = tz = 0; i < sa_n; i++) { + if((sa[i]>>32) != ((uint32_t)-1)) continue; + if(tz < i) { + bs = sa[i]; sa[i] = sa[tz]; sa[tz] = bs; + } + tz++; + } + sa_n = tz; + } + + if(sa_n > 0 && sb_n > 0) { + radix_sort_bc64(sa, sa + sa_n); radix_sort_bc64(sb, sb + sb_n); + + overlap_region *zm = NULL, *rm = NULL; uint64_t zs, ze, ob, zsc, zcn, rs, re, os, oe, oi, rr, kn, cn = cl->length, cs, ce; uint8_t f; + for (m = tz = 0; m < sa_n; m++) { + zm = &(ol->list[(uint32_t)sa[m]]); + zs = zm->x_pos_s; ze = zm->x_pos_e + 1; + ob = (ze - zs)*0.95; if(ob < 16) ob = 16; + zsc = osc[(uint32_t)sa[m]]*ch_sc; + zcn = ocn[(uint32_t)sa[m]]*ch_occ; + + for (k = f = 0; (k < sb_n) && (ze > ol->list[(uint32_t)sb[k]].x_pos_s); k++) { + rm = &(ol->list[(uint32_t)sb[k]]); + if(osc[(uint32_t)sb[k]] < zsc) continue; + if(ocn[(uint32_t)sb[k]] < zcn) continue; + rs = rm->x_pos_s; re = rm->x_pos_e + 1; + os = ((rs>=zs)?rs:zs); oe = ((re<=ze)?re:ze); + if((oe > os) && (oe - os) >= ob) { + oi = rm->shared_seed; rr = cl->list[oi].readID; kn = 0; + for (; (oi < cn) && (cl->list[oi].readID == rr) && (kn < zcn); oi++) { + ce = cl->list[oi].self_offset; cs = ce - (cl->list[oi].cnt&(0xffu)); + if((cs >= os) && (ce <= oe)) kn++; + } + if(kn >= zcn) { + f = 1; break; + } + } + } + if(f) continue; + + if(tz < m) { + bs = sa[m]; sa[m] = sa[tz]; sa[tz] = bs; + } + sa[tz] = (uint32_t)sa[tz]; sa[tz] = (((uint64_t)zm->x_pos_s)<<32)|(sa[tz]); tz++; + tot_b += zm->x_pos_e + 1 - zm->x_pos_s; + zm->non_homopolymer_errors = UINT32_MAX - 1;///primary chain that needs to be verfied + } + wma_n += tz; + } + } + + + if(wma_n > 0) { + radix_sort_bc64(wsrt, wsrt + wma_n); + for (k = 1, i = 0; k <= wma_n; k++) { + if (k == wma_n || (wsrt[k]>>32) != (wsrt[i]>>32)) { + if(k - i > 1) { + for (tz = i; tz < k; tz++) { + z = &(ol->list[(uint32_t)wsrt[tz]]); + m = z->x_pos_e + 1; + m <<= 32; m += (uint32_t)wsrt[tz]; wsrt[tz] = m; + } + radix_sort_bc64(wsrt + i, wsrt + k); + } + i = k; + } + } + } + + (*rwsrt_n) = wma_n; + return wsrt; +} + uint64_t gen_hc_r_alin_adp_mmp_0(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *sp, uint64_t ocw, uint8_t *hpf, asg32_v *v32, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, uint8_t set_match) @@ -37156,7 +38553,7 @@ uint64_t gen_hc_r_alin_adp_mmp_0(overlap_region_alloc* ol, Candidates_list *cl, bs = (w.window_length)+(THRESHOLD_MAX_SIZE<<1)+1; resize_UC_Read(tu, bs<<1); spn0 = sp->n; - wsrt = mmp_chn_select(ol, cl, sp, ocw, ocn, osc, ql, &wsrt_n, set_match, max_n_chain, max_n_chain_f, chain_cutoff, ave_cov_min, 0.333333, 16, 16); + wsrt = /**mmp_chn_select**/mmp_chn_select_adv(ol, cl, sp, ocw, ocn, osc, ql, &wsrt_n, set_match, max_n_chain, max_n_chain_f, chain_cutoff, ave_cov_min, 0.333333, 16, 16); for (i = 0; i < wsrt_n; i++) { z = &(ol->list[(uint32_t)wsrt[i]]); @@ -37218,7 +38615,7 @@ uint64_t gen_gc_r_alin_adp_mmp_0(overlap_region_alloc* ol, Candidates_list *cl, bs = (w.window_length)+(THRESHOLD_MAX_SIZE<<1)+1; resize_UC_Read(tu, bs<<1); spn0 = sp->n; - wsrt = mmp_chn_select(ol, cl, sp, ocw, ocn, osc, ql, &wsrt_n, set_match, max_n_chain, max_n_chain_f, chain_cutoff, ave_cov_min, 0.333333, 16, 16); + wsrt = /**mmp_chn_select**/mmp_chn_select_adv(ol, cl, sp, ocw, ocn, osc, ql, &wsrt_n, set_match, max_n_chain, max_n_chain_f, chain_cutoff, ave_cov_min, 0.333333, 16, 16); for (i = 0; i < wsrt_n; i++) { z = &(ol->list[(uint32_t)wsrt[i]]); @@ -37618,11 +39015,12 @@ uint64_t gen_hc_r_alin_adp_mmp_1(overlap_region_alloc* ol, Candidates_list *cl, bs = (w.window_length)+(THRESHOLD_MAX_SIZE<<1)+1; resize_UC_Read(tu, bs<<1); spn0 = sp->n; nwl = w.window_length; - wsrt = mmp_chn_select(ol, cl, sp, ocw, ocn, osc, ql, &wsrt_n, set_match, max_n_chain, max_n_chain_f, chain_cutoff, ave_cov_min, 0.333333, 16, 16); + wsrt = /**mmp_chn_select**/mmp_chn_select_adv(ol, cl, sp, ocw, ocn, osc, ql, &wsrt_n, set_match, max_n_chain, max_n_chain_f, chain_cutoff, ave_cov_min, 0.333333, 16, 16); for (i = 0; i < wsrt_n; i++) { z = &(ol->list[(uint32_t)wsrt[i]]); + // fprintf(stderr, "-s-[M::%s]\tqid::%u\ttid::%u\t\n", __func__, z->x_id, z->y_id); if(z->is_match == 0) { z->w_list.n = 0; z->align_length = 0; } @@ -38325,7 +39723,7 @@ void gen_hc_r_alin_adv_adp_smp_0(gen_hc_aln_t *ez, uint8_t set_match) bs = (MAX(ez->wl[0], ez->wl[1]))+(THRESHOLD_MAX_SIZE<<1)+1; resize_UC_Read(ez->tu, bs<<1); spn0 = sp->n; - wsrt = mmp_chn_select(ez->ol, ez->cl, sp, ez->ocw, ocn, osc, ql, &wsrt_n, set_match, ez->max_n_chain, ez->max_n_chain_f, ez->chain_cutoff, ez->ave_cov_min, 0.333333, 16, 16); + wsrt = /**mmp_chn_select**/mmp_chn_select_adv(ez->ol, ez->cl, sp, ez->ocw, ocn, osc, ql, &wsrt_n, set_match, ez->max_n_chain, ez->max_n_chain_f, ez->chain_cutoff, ez->ave_cov_min, 0.333333, 16, 16); for (i = 0; i < wsrt_n; i++) { z = &(ez->ol->list[(uint32_t)wsrt[i]]); @@ -38452,7 +39850,7 @@ void gen_hc_r_alin_adv_adp_smp_1(gen_hc_aln_t *ez, uint8_t set_match) bs = (MAX(ez->wl[0], ez->wl[1]))+(THRESHOLD_MAX_SIZE<<1)+1; resize_UC_Read(ez->tu, bs<<1); spn0 = sp->n; - wsrt = mmp_chn_select(ez->ol, ez->cl, sp, ez->ocw, ocn, osc, ql, &wsrt_n, set_match, ez->max_n_chain, ez->max_n_chain_f, ez->chain_cutoff, ez->ave_cov_min, 0.333333, 16, 16); + wsrt = /**mmp_chn_select**/mmp_chn_select_adv(ez->ol, ez->cl, sp, ez->ocw, ocn, osc, ql, &wsrt_n, set_match, ez->max_n_chain, ez->max_n_chain_f, ez->chain_cutoff, ez->ave_cov_min, 0.333333, 16, 16); for (i = nol_1 = 0; i < wsrt_n; i++) { diff --git a/Correct.h b/Correct.h index 3889a72..ee0f352 100644 --- a/Correct.h +++ b/Correct.h @@ -1458,7 +1458,7 @@ void gen_hc_r_alin_nec_adv(gen_hc_aln_t *ez); uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf); uint64_t gen_hc_r_alin_self(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, asg16_v* scc); void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, asg16_v *qc0, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, Chain_Data *dp, asg8_v *q8, asg8_v *t8, uint8_t lindel, uint64_t tcut, uint64_t site_sc, int64_t h0_w, asg32_v *b32, - int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate, overlap_region *rchn, int64_t rref_len); + int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate, overlap_region *rchn, int64_t rref_len, int64_t flag_hf_ov_cut, int64_t min_re_cut, double min_re_rt, int64_t min_hpc_re_cut, double min_hpc_re_rt, int8_t is_hpc_flt); void set_exact_exz(bit_extz_t *exz, int64_t qs, int64_t qe, int64_t ts, int64_t te); void push_alnw(overlap_region *aux_o, bit_extz_t *exz); void cal_exz_global(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez); diff --git a/Overlaps.cpp b/Overlaps.cpp index 2799503..b316883 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -178,6 +178,12 @@ typedef struct { int64_t min_sc, penalty, max_drop; } telo_end_pip_t; +typedef struct { + ma_hit_t_alloc* src; + ma_hit_t_alloc* rec; + int64_t n_thread, n_a; +} rclean_aux; + ///this value has been updated at the first line of build_string_graph_without_clean long long min_thres; @@ -188,6 +194,24 @@ kv_u_trans_t *get_utg_ovlp(ma_ug_t **ug, asg_t* read_g, ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t, uint8_t* r_het); void delete_useless_nodes(ma_ug_t **ug); +static void reset_ma_hit_t_0(void *data, long i, int tid) // callback for kt_for() +{ + rclean_aux *p = (rclean_aux *)data; uint64_t k; + for (k = 0; k < p->src[i].length; k++) { + p->src[i].buffer[k].del = 0; + } + for (k = 0; k < p->rec[i].length; k++) { + p->rec[i].buffer[k].del = 0; + } +} + +void inital_ovlap(ma_hit_t_alloc *src, ma_hit_t_alloc *rec, int64_t n_a, int64_t n_thre) +{ + rclean_aux p; memset(&p, 0, sizeof(p)); + p.n_thread = n_thre; p.src = src; p.rec = rec; p.n_a = n_a; + kt_for(p.n_thread, reset_ma_hit_t_0, &p, p.n_a); +} + static void mark_telo_ends(void *data, long i, int tid) // callback for kt_for() { telo_end_pip_t *sl = (telo_end_pip_t *)data; @@ -23758,7 +23782,7 @@ void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sou } int load_debug_graph(asg_t** sg, ma_hit_t_alloc** sources, ma_sub_t** coverage_cut, -char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_ul_t *ul_r_inf); +char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_ul_t *ul_r_inf, uint8_t **cmk); int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, char* output_file_name) { char* gfa_name = (char*)malloc(strlen(output_file_name)+25); @@ -23776,7 +23800,7 @@ int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_s } } - if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_debug_graph(NULL, NULL, NULL, output_file_name, NULL, NULL, NULL)) { + if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_debug_graph(NULL, NULL, NULL, output_file_name, NULL, NULL, NULL, NULL)) { (*sources) = NULL; (*reverse_sources) = NULL; free(gfa_name); @@ -32191,6 +32215,42 @@ int load_coverage_cut(ma_sub_t** coverage_cut, char* read_file_name) } +int write_cmk(uint8_t* cmk, char* read_file_name, uint64_t n_read) +{ + char* index_name = (char*)malloc(strlen(read_file_name)+15); + sprintf(index_name, "%s.bin", read_file_name); + FILE* fp = fopen(index_name, "w"); + fwrite(&n_read, sizeof(n_read), 1, fp); + fwrite(cmk, sizeof((*(cmk))), n_read, fp); + free(index_name); + fflush(fp); + fclose(fp); + + return 1; +} + +int load_cmk(uint8_t** cmk, char* read_file_name) +{ + char* index_name = (char*)malloc(strlen(read_file_name)+15); + sprintf(index_name, "%s.bin", read_file_name); + FILE* fp = fopen(index_name, "r"); + if(!fp) + { + return 0; + } + int f_flag = 0; + uint64_t n_read; + f_flag += fread(&n_read, sizeof(n_read), 1, fp); + (*cmk) = (uint8_t*)malloc(sizeof((*(*(cmk))))*n_read); + + fread((*cmk), sizeof((*((*cmk)))), n_read, fp); + + free(index_name); + fflush(fp); + fclose(fp); + return 1; +} + int write_coverage_cut(ma_sub_t* coverage_cut, char* read_file_name, uint64_t n_read) { char* index_name = (char*)malloc(strlen(read_file_name)+15); @@ -32441,7 +32501,7 @@ int load_asg_t(asg_t **sg, char* read_file_name) } int write_debug_graph(asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, -char* output_file_name, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, all_ul_t *ul_r_inf) +char* output_file_name, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, all_ul_t *ul_r_inf, uint8_t *cmk) { char* gfa_name = (char*)malloc(strlen(output_file_name)+55); @@ -32465,6 +32525,10 @@ char* output_file_name, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, all_ul sprintf(gfa_name, "%s.all.debug.ul.rinfor", output_file_name); write_all_ul_t(ul_r_inf, gfa_name, NULL); } + if(cmk) { + sprintf(gfa_name, "%s.all.debug.cmk", output_file_name); + write_cmk(cmk, gfa_name, R_INF.total_reads); + } free(gfa_name); fprintf(stderr, "debug_graph has been written.\n"); return 1; @@ -32472,7 +32536,7 @@ char* output_file_name, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, all_ul int load_debug_graph(asg_t** sg, ma_hit_t_alloc** sources, ma_sub_t** coverage_cut, -char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_ul_t *ul_r_inf) +char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_ul_t *ul_r_inf, uint8_t **cmk) { FILE* fp = NULL; char* gfa_name = (char*)malloc(strlen(output_file_name)+55); @@ -32494,6 +32558,11 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_u sprintf(gfa_name, "%s.all.debug.ul.rinfor.ul.ovlp.bin", output_file_name); fp = fopen(gfa_name, "r"); if(!fp) return 0; fclose(fp); } + + if(cmk) { + sprintf(gfa_name, "%s.all.debug.cmk.bin", output_file_name); + fp = fopen(gfa_name, "r"); if(!fp) return 0; fclose(fp); + } if((sources == NULL) || (reverse_sources == NULL) || (ruIndex == NULL)) { @@ -32567,6 +32636,11 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex, all_u if(!load_all_ul_t(ul_r_inf, gfa_name, &R_INF, NULL)) return 0; } + if(cmk) { + sprintf(gfa_name, "%s.all.debug.cmk", output_file_name); + if(!load_cmk(cmk, gfa_name)) return 0; + } + R_INF.paf = (*sources); R_INF.reverse_paf = (*reverse_sources); return 1; @@ -39767,19 +39841,24 @@ ma_hit_t_alloc* src, uint64_t* readLen, R_to_U* ruIndex, bub_label_t *b_mask_t, // prt_dbg_rid_ovlp(src, -1, (char*)"c85c2e91-0490-438b-977b-b7d056973996", "2"); ma_hit_sub(min_dp, src, n_read, readLen, mini_overlap_length, cov); // prt_dbg_rid_ovlp(src, -1, (char*)"c85c2e91-0490-438b-977b-b7d056973996", "3"); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "init_0", src, NULL, NULL, NULL); detect_chimeric_reads(src, n_read, readLen, *cov, asm_opt.max_ov_diff_final*2.0, ul, UL_COV_THRES); // prt_dbg_rid_ovlp(src, -1, (char*)"c85c2e91-0490-438b-977b-b7d056973996", "4"); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "init_1", src, NULL, *cov, NULL); ma_hit_cut(src, n_read, readLen, mini_overlap_length, cov); // prt_dbg_rid_ovlp(src, -1, (char*)"c85c2e91-0490-438b-977b-b7d056973996", "5"); ma_hit_flt(src, n_read, *cov, max_hang_length, mini_overlap_length); // prt_dbg_rid_ovlp(src, -1, (char*)"c85c2e91-0490-438b-977b-b7d056973996", "6"); ma_hit_contained_advance(src, n_read, *cov, ruIndex, max_hang_length, mini_overlap_length); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "init_2", src, ruIndex, *cov, NULL); // prt_dbg_rid_ovlp(src, -1, (char*)"7b70a587-f56c-48ac-adf8-b98a67365063_2", "7"); if(!ul) { sg = ma_sg_gen(src, n_read, *cov, max_hang_length, mini_overlap_length); - if(asm_opt.prt_dbg_gfa) prt_dbg_gfa(sg, "raw", *cov, src, ruIndex, max_hang_length, mini_overlap_length); + if(asm_opt.prt_dbg_gfa) prt_dbg_gfa(sg, "raw", *cov, src, ruIndex, max_hang_length, mini_overlap_length); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "init_3", NULL, ruIndex, NULL, sg); asg_arc_del_trans(sg, gap_fuzz); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "init_4", NULL, ruIndex, NULL, sg); } else { ug_opt_t uopt; sg = ma_sg_gen_ul(src, n_read, *cov, ruIndex, max_hang_length, mini_overlap_length, UL_COV_THRES); @@ -39871,6 +39950,17 @@ float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr, ma_sub_t **coverage_cut_ptr, uint8_t *cmk, int debug_g) { + // const char *dbg_ids[] = { + // "6152e53b-a2ec-4f23-898f-462c9dee8e0f", + // "c5bbf86e-3146-4757-9882-55a5baebf770", + // }; + // const uint64_t dbg_ids_n[] = { + // 3177240, + // 2198541, + // }; + + // hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "test_0", sources, NULL, NULL, NULL); + char *o_file = get_outfile_name(output_file_name); ma_sub_t *coverage_cut = *coverage_cut_ptr; asg_t *sg = *sg_ptr; @@ -39894,8 +39984,10 @@ ma_sub_t **coverage_cut_ptr, uint8_t *cmk, int debug_g) ///normalize_ma_hit_t_single_side(sources, n_read); normalize_ma_hit_t_single_side_advance(sources, n_read, asm_opt.is_ont, cmk); + // hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "test_1", sources, NULL, NULL, NULL); // normalize_ma_hit_t_single_side_advance_mult(sources, n_read, asm_opt.thread_num); normalize_ma_hit_t_single_side_advance(reverse_sources, n_read, 0, cmk); + // hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "test_2", sources, NULL, NULL, NULL); // normalize_ma_hit_t_single_side_advance_mult(reverse_sources, n_read, asm_opt.thread_num); if (ha_opt_triobin(&asm_opt)) @@ -39925,6 +40017,21 @@ ma_sub_t **coverage_cut_ptr, uint8_t *cmk, int debug_g) // prt_dbg_rid_ovlp(sources, 27087, NULL, "0-b"); // prt_specific_overlap(sources, 22233, 22235, "0-c"); // prt_specific_overlap(sources, 22235, 22233, "0-c"); + + + + ///@brief debug + if (asm_opt.flag & HA_F_VERBOSE_GFA) { + write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF, cmk); + debug_gfa:; + } + + + + // hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "ini_s", sources, NULL, NULL, NULL); + // hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "ini_r", reverse_sources, NULL, NULL, NULL); + // exit(1); + sg = gen_init_sg(min_dp, n_read, mini_overlap_length, max_hang_length, gap_fuzz, sources, readLen, ruIndex, &b_mask_t, &coverage_cut, asm_opt.ar?&UL_INF:NULL, te); // if(asm_opt.ar) exit(1); @@ -39955,24 +40062,29 @@ ma_sub_t **coverage_cut_ptr, uint8_t *cmk, int debug_g) free(unlean_name); } - - // if (asm_opt.flag & HA_F_VERBOSE_GFA) { - // write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF); - // debug_gfa:; - // } + /** + ///@brief debug + if (asm_opt.flag & HA_F_VERBOSE_GFA) { + write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF, cmk); + debug_gfa:; + } + **/ + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_z", sources, ruIndex, coverage_cut, sg); gen_ug_opt_t(&uopt, sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, &b_mask_t, te); ul_clean_gfa(&uopt, sg, sources, reverse_sources, ruIndex, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.6, asm_opt.max_short_tip, gap_fuzz, &b_mask_t, !!asm_opt.ar, ha_opt_triobin(&asm_opt), UL_COV_THRES, cmk, o_file); + /** ///@brief debug if (asm_opt.flag & HA_F_VERBOSE_GFA) { - write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF); + write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex, &UL_INF, cmk); debug_gfa:; gen_ug_opt_t(&uopt, sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, &b_mask_t, te); // set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length); } + **/ if(asm_opt.ar) { gradually_renew_g(&sources, &reverse_sources, &n_read, &readLen, &coverage_cut, ruIndex, @@ -40229,7 +40341,7 @@ long long bubble_dist, int read_graph, int write) min_thres = asm_opt.max_short_tip + 1; if (asm_opt.flag & HA_F_VERBOSE_GFA) { - if(load_debug_graph(/**NULL**/&sg, &sources, /**NULL**/&coverage_cut, output_file_name, &reverse_sources, &ruIndex, &UL_INF)) + if(load_debug_graph(NULL/**&sg**/, &sources, NULL/**&coverage_cut**/, output_file_name, &reverse_sources, &ruIndex, &UL_INF, (asm_opt.is_ont)?(&cmk):(NULL))) { fprintf(stderr, "debug gfa has been loaded\n"); @@ -40249,20 +40361,39 @@ long long bubble_dist, int read_graph, int write) } ///debug_info_of_specfic_read("m64011_190830_220126/31720629/ccs", sources, reverse_sources, -1, "beg"); + // const char *dbg_ids[] = { + // "6152e53b-a2ec-4f23-898f-462c9dee8e0f", + // "c5bbf86e-3146-4757-9882-55a5baebf770", + // }; + // const uint64_t dbg_ids_n[] = { + // 3177240, + // 2198541, + // }; + + if (!(asm_opt.flag & HA_F_BAN_ASSEMBLY)) { // if(asm_opt.is_ont) handle_chemical_r(asm_opt.thread_num, R_INF.total_reads); if(asm_opt.is_ont) { - uint64_t i = 0, j = 0; + // uint64_t i = 0, j = 0; memset(ruIndex.index, -1, sizeof(uint32_t)*(ruIndex.len)); - for (i = 0; i < n_read; i++) { - for (j = 0; j < sources[i].length; j++) sources[i].buffer[j].del = 0; - for (j = 0; j < reverse_sources[i].length; j++) reverse_sources[i].buffer[j].del = 0; - } + inital_ovlap(sources, reverse_sources, n_read, asm_opt.thread_num); + + + // hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "ttst_0", sources, NULL, NULL, NULL); + // hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "trev_0", reverse_sources, NULL, NULL, NULL); + if(asm_opt.is_ont) cmk = gen_chemical_arc_rf(asm_opt.thread_num, R_INF.total_reads); + + // hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "ttst_1", sources, NULL, NULL, NULL); + if(asm_opt.del_hf) clean_arc_rf(asm_opt.thread_num, R_INF.total_reads); + + // hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "ttst_2", sources, NULL, NULL, NULL); } try_rescue_overlaps(sources, reverse_sources, n_read, 4, asm_opt.is_ont); + + // hc_dbg_prt_ma_hit_t(dbg_ids, sizeof(dbg_ids)/sizeof(dbg_ids[0]), dbg_ids_n, sizeof(dbg_ids_n)/sizeof(dbg_ids_n[0]), "ttst_3", sources, NULL, NULL, NULL); clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length, max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio, diff --git a/ecovlp.cpp b/ecovlp.cpp index 03be9ec..4635ec7 100644 --- a/ecovlp.cpp +++ b/ecovlp.cpp @@ -218,8 +218,8 @@ cc_v scc = {0, 0, NULL, NULL, NULL, 0}; cc_v scb = {0, 0, NULL, NULL, NULL, 0}; cc_v sca = {0, 0, NULL, NULL, NULL, 0}; -typedef struct {uint64_t p, pn, pm, tov, tov_size, tqn; asg64_v *idx; ma_hit_t_alloc *pf;} tsrt_v_buf; -typedef struct {uint64_t p, pn, pm, rid, tot, chunk_size, tqn, n_thr; uint64_t n_ov, n_bl; ma_hit_t_alloc *pf;} tsrt_v_m; +typedef struct {uint64_t p, pn, pm, tov, tov_size, tqn, tot; asg64_v *idx; ma_hit_t_alloc *pf;} tsrt_v_buf; +typedef struct {uint64_t p, pn, pm, rid, tot, chunk_size, tqn, n_thr; uint64_t n_ov[2], n_bl; ma_hit_t_alloc *pf;} tsrt_v_m; @@ -2672,7 +2672,7 @@ void print_debug_ovlp_cigar(overlap_region_alloc* ol, asg64_v* idx, kv_ul_ov_t * } uint64_t wcns_gen(overlap_region_alloc* ol, All_reads *rref, uint64_t qid, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, uint64_t wl, int64_t ql, uint64_t occ_tot, double occ_exact, overlap_region *aux_o, asg32_v* b32, cns_gfa *cns, uint64_t cns_g_wl, uint32_t rid, uint64_t tcut, - uint64_t tot_ont_b, uint64_t tot_hf_b, uint64_t ont_rate_w, uint64_t hf_rate_w, uint64_t hf_rate_w_max, asg64_v *hf_idx) + uint64_t tot_ont_b, uint64_t tot_hf_b, uint64_t ont_rate_w, uint64_t hf_rate_w, uint64_t hf_rate_w_max, asg64_v *hf_idx, uint8_t reset_no_l_idel) { int64_t on = ol->length, k, i, zwn, q[2]; cns->cns_g_wl = cns_g_wl; uint64_t m, *ra, rn, nec = 0, n_id, l_nid, p[2], li; uint64_t o_rate = ((uint64_t)-1), h_rate = ((uint64_t)-1); overlap_region *z; ul_ov_t *cp; @@ -2684,9 +2684,10 @@ uint64_t wcns_gen(overlap_region_alloc* ol, All_reads *rref, uint64_t qid, UC_Re } for (k = idx->n = c_idx->n = 0; k < on; k++) { - z = &(ol->list[k]); zwn = z->w_list.n; z->without_large_indel = l_nid = 0; + z = &(ol->list[k]); zwn = z->w_list.n; l_nid = 0; + if(reset_no_l_idel) z->without_large_indel = 0; if((!zwn) || (z->is_match != 1)) continue; - hf = ((z->y_id >= tcut)?(1):(0)); + hf = ((z->y_id >= tcut)?(1):(0)); z->without_large_indel = 0; for (i = 0, li = (uint64_t)-1; i < zwn; i++) { if(is_ualn_win(z->w_list.a[i])) { n_id = z->w_list.a[i].x_end + 1 - z->w_list.a[i].x_start; @@ -3179,6 +3180,7 @@ void push_ne_ovlp(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag, // fprintf(stderr, "@%s\tSN:%.*s(id::%u)\terr::%u\n", flag==1?"SQ":"RQ", (int32_t)Get_NAME_LENGTH((*R_INF), ov->list[k].y_id), Get_NAME((*R_INF), ov->list[k].y_id), ov->list[k].y_id, ov->list[k].non_homopolymer_errors); z = &(paf->buffer[paf->length++]); + z->del = 0; z->qns = ov->list[k].x_id; z->qns = z->qns << 32; @@ -3195,6 +3197,7 @@ void push_ne_ovlp(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag, z->bl = Get_READ_LENGTH((*R_INF), ov->list[k].y_id); z->ml = ov->list[k].strong; z->no_l_indel = ov->list[k].without_large_indel; + z->cc = 0; z->el = 0; if(ec) { @@ -3291,6 +3294,7 @@ void push_ne_ovlp_syn(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t fl // fprintf(stderr, "@%s\tSN:%.*s(id::%u)\terr::%u\n", flag==1?"SQ":"RQ", (int32_t)Get_NAME_LENGTH((*R_INF), ov->list[k].y_id), Get_NAME((*R_INF), ov->list[k].y_id), ov->list[k].y_id, ov->list[k].non_homopolymer_errors); z = &(paf->buffer[paf->length++]); + z->del = 0; z->qns = ov->list[k].x_id; z->qns = z->qns << 32; @@ -3743,6 +3747,165 @@ uint64_t gen_hc_r_alin_ea_flt_mmp(overlap_region_alloc* ol, Candidates_list *cl, return tot_b; } +uint8_t is_match_ov(overlap_region *z, ma_hit_t *p, double ov_rate) +{ + int64_t pq[2], pt[2], zq[2], zt[2], os, oe, ovlp; + + // if(p->cc == 0x3fffffffu) { + // fprintf(stderr, "\n-p-[M::%s] qn::%lu, q::[%u,%u), %c, tn::%u, t::[%u,%u)\n", __func__, p->qns>>32, (uint32_t)p->qns, p->qe, "+-"[p->rev], p->tn, p->ts, p->te); + // fprintf(stderr, "-z-[M::%s] qn::%u, q::[%u,%u), %c, tn::%u, t::[%u,%u)\n", __func__, z->x_id, z->x_pos_s, z->x_pos_e + 1, "+-"[z->y_pos_strand], z->y_id, z->y_pos_s, z->y_pos_e + 1); + // } + + pq[0] = (uint32_t)p->qns; pq[1] = p->qe; + pt[0] = p->ts; pt[1] = p->te; + + zq[0] = z->x_pos_s; zq[1] = z->x_pos_e + 1; + zt[0] = z->y_pos_s; zt[1] = z->y_pos_e + 1; + + os = MAX(pq[0], zq[0]); oe = MIN(pq[1], zq[1]); + ovlp = ((oe>os)? (oe-os):0); + if(!((ovlp) && (ovlp >= ((pq[1] - pq[0])*ov_rate)) && ((ovlp >= ((zq[1] - zq[0])*ov_rate))))) return 0; + + os = MAX(pt[0], zt[0]); oe = MIN(pt[1], zt[1]); + ovlp = ((oe>os)? (oe-os):0); + if(!((ovlp) && (ovlp >= ((pt[1] - pt[0])*ov_rate)) && ((ovlp >= ((zt[1] - zt[0])*ov_rate))))) return 0; + + // if(p->cc == 0x3fffffffu) { + // fprintf(stderr, "-m-[M::%s]\n", __func__); + // } + return 1; +} + +uint64_t slash_overlap(overlap_region* za, uint64_t *ei, uint64_t en, uint64_t *oi, uint64_t on, ma_hit_t_alloc *in0, uint64_t icn, ma_hit_t_alloc *in1, bit_extz_t *exz, UC_Read* qu, UC_Read* tu, All_reads *rref, uint8_t post_syn) +{ + uint64_t ko, lo, ke, le, io, ie, nec = 0, tid, trev; uint8_t exc; overlap_region *z; ma_hit_t *p; + le = ke = 0; + for (ko = 1, lo = 0; ko <= on; ko++) { + if ((ko == on) || ((oi[ko]>>32) != (oi[lo]>>32))) { + tid = oi[lo]>>33; trev = (oi[lo]>>32)&1; + for (; (ke < en) && ((ei[ke]>>32)<((tid<<1)|trev)); ke++); + if((ke < en) && ((ei[ke]>>32)==((tid<<1)|trev))) { + le = ke; + for (; (ke < en) && ((ei[ke]>>32)==((tid<<1)|trev)); ke++); + ///oi[lo, ko) & ei[le, ke) + for (io = lo; io < ko; io++) { + z = &(za[(uint32_t)oi[io]]); z->is_match = 0; + for (ie = le; ie < ke; ie++) { + exc = ei[ie]&1; + if((((uint32_t)ei[ie])>>1) < icn) { + p = &(in0->buffer[(((uint32_t)ei[ie])>>1)]); + } else { + p = &(in1->buffer[(((uint32_t)ei[ie])>>1)-icn]); + } + if((z->x_pos_s == ((uint32_t)p->qns)) && (z->x_pos_e + 1 == p->qe) && + (z->y_pos_s == p->ts) && (z->y_pos_e + 1 == p->te)) { + if(post_syn) z->is_match = 2; + if(exc) { + resize_UC_Read(tu, p->te - p->ts); recover_UC_Read_sub_region(tu->seq, p->ts, p->te - p->ts, trev, rref, tid); + if(exact_ec_check(qu->seq, qu->length, tu->seq, p->te - p->ts, ((uint32_t)p->qns), p->qe, 0, p->te - p->ts)) { + z->is_match = 1; z->shared_seed = z->non_homopolymer_errors;///for index + z->non_homopolymer_errors = 0; z->strong = z->without_large_indel = 0; + set_exact_exz(exz, z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1); push_alnw(z, exz); + nec++; + break; + } + } + if(post_syn) break; + } + + if(post_syn && is_match_ov(z, p, 0.666666)) { + z->is_match = 2; + break; + } + } + } + } + lo = ko; + } + } + + return nec; +} + +uint64_t gen_hc_r_alin_ea_flt_mmp_adv(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in0, ma_hit_t_alloc *in1, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, + uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *bp, uint64_t ocw, asg8_v *hpz, asg32_v *v32, uint8_t post_syn) +{ + if(ol->length <= 0) return 0; + + uint64_t k, tot_b = 0, icn; v32->n = 0; uint64_t i, m, *ei, en, *oi, on, nec; overlap_region *z; + v32->n = ol->length<<1; kv_resize(uint32_t, *v32, v32->n); + for (i = 0; i < ol->length; i++) { + v32->a[i] = ol->list[i].align_length; + v32->a[i+ol->length] = ol->list[i].shared_seed; + ol->list[i].align_length = 0; + } + + srt->n = 0; + if(post_syn) { + for (k = 0; k < in0->length; k++) { + m = in0->buffer[k].tn; m <<= 1; m |= in0->buffer[k].rev; m <<= 32; + m |= (k<<1); if(in0->buffer[k].el) m |= 1; + kv_push(uint64_t, (*srt), m); + } + } else { + for (k = 0; k < in0->length; k++) { + if(!(in0->buffer[k].el)) continue; + m = in0->buffer[k].tn; m <<= 1; m |= in0->buffer[k].rev; m <<= 32; + m |= ((k<<1) + 1); + kv_push(uint64_t, (*srt), m); + } + } + + icn = in0->length; + + if(post_syn) { + for (k = 0; k < in1->length; k++) { + if(in1->buffer[k].no_l_indel == 0) continue; + m = in1->buffer[k].tn; m <<= 1; m |= in1->buffer[k].rev; m <<= 32; + m |= ((k+icn)<<1); kv_push(uint64_t, (*srt), m); + } + } + + + if(!(srt->n)) { + // gen_ff_hpc(hpz, qu->seq, qu->length, HPC_RR_Q, HPC_CC_Q); + if(asm_opt.simd_mm > 0) { + tot_b = gen_hc_r_alin_adp_mmp_1(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max, sec_aln_win, sec_aln_cov, sec_aln_err_rate, sec_aln_max, srt, ocw, NULL/**hpz->a**/, + v32, bp, max_n_chain>0?max_n_chain:1, max_n_chain_f>0?max_n_chain_f:1, chain_cutoff, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), 1); + } else { + tot_b = gen_hc_r_alin_adp_mmp_0(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max, sec_aln_win, sec_aln_cov, sec_aln_err_rate, sec_aln_max, srt, ocw, NULL/**hpz->a**/, + v32, bp, max_n_chain>0?max_n_chain:1, max_n_chain_f>0?max_n_chain_f:1, chain_cutoff, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), 1); + } + + } else { + kv_resize(uint64_t, *srt, (srt->n + ol->length)); + ei = srt->a; en = srt->n; oi = srt->a + srt->n; on = ol->length; + for (k = 0; k < on; k++) { + z = &(ol->list[k]); z->is_match = z->strong = z->without_large_indel = 0; + oi[k] = z->y_id; oi[k] <<= 1; oi[k] |= z->y_pos_strand; + oi[k] <<= 32; oi[k] |= k; + } + + radix_sort_ec64(ei, ei + en); radix_sort_ec64(oi, oi + on); + nec = slash_overlap(ol->list, ei, en, oi, on, in0, icn, in1, exz, qu, tu, rref, post_syn); + + if(on > nec) { + // gen_ff_hpc(hpz, qu->seq, qu->length, HPC_RR_Q, HPC_CC_Q); + if(asm_opt.simd_mm > 0) { + tot_b = gen_hc_r_alin_adp_mmp_1(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max, sec_aln_win, sec_aln_cov, sec_aln_err_rate, sec_aln_max, srt, ocw, NULL/**hpz->a**/, + v32, bp, max_n_chain>0?max_n_chain:1, max_n_chain_f>0?max_n_chain_f:1, chain_cutoff, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), 0); + } else { + tot_b = gen_hc_r_alin_adp_mmp_0(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max, sec_aln_win, sec_aln_cov, sec_aln_err_rate, sec_aln_max, srt, ocw, NULL/**hpz->a**/, + v32, bp, max_n_chain>0?max_n_chain:1, max_n_chain_f>0?max_n_chain_f:1, chain_cutoff, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), 0); + } + } + } + + if(ol->length) srt_olst(ol); + + return tot_b; +} + uint64_t gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *kp, asg8_v *hpz) { @@ -4088,6 +4251,79 @@ void gen_hc_r_alin_ea_adv_flt_mmp(gen_hc_aln_t *ez) } +void gen_hc_r_alin_ea_adv_flt_mmp_adv(gen_hc_aln_t *ez, ma_hit_t_alloc *in0, ma_hit_t_alloc *in1, uint8_t post_syn) +{ + if(ez->ol->length <= 0) return; + + uint64_t i, k, m, *ei, en, *oi, on, icn, nec; overlap_region *z; + ez->v32->n = ez->ol->length<<1; kv_resize(uint32_t, *(ez->v32), ez->v32->n); + for (i = 0; i < ez->ol->length; i++) { + ez->v32->a[i] = ez->ol->list[i].align_length; + ez->v32->a[i+ez->ol->length] = ez->ol->list[i].shared_seed; + ez->ol->list[i].align_length = 0; + } + + ez->srt->n = 0; + if(post_syn) { + for (k = 0; k < in0->length; k++) { + m = in0->buffer[k].tn; m <<= 1; m |= in0->buffer[k].rev; m <<= 32; + m |= (k<<1); if(in0->buffer[k].el) m |= 1; + kv_push(uint64_t, (*(ez->srt)), m); + } + } else { + for (k = 0; k < in0->length; k++) { + if(!(in0->buffer[k].el)) continue; + m = in0->buffer[k].tn; m <<= 1; m |= in0->buffer[k].rev; m <<= 32; + m |= ((k<<1) + 1); + kv_push(uint64_t, (*(ez->srt)), m); + } + } + + icn = in0->length; + + if(post_syn) { + for (k = 0; k < in1->length; k++) { + if(in1->buffer[k].no_l_indel == 0) continue; + m = in1->buffer[k].tn; m <<= 1; m |= in1->buffer[k].rev; m <<= 32; + m |= ((k+icn)<<1); kv_push(uint64_t, (*(ez->srt)), m); + } + } + + + + if(!(ez->srt->n)) { + // gen_hc_r_alin_adv_adp_smp(ez, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 1); + if(asm_opt.simd_mm > 0) gen_hc_r_alin_adv_adp_smp_1(ez, 1); + else gen_hc_r_alin_adv_adp_smp_0(ez, 1); + } else { + ///debug for memory + // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); + kv_resize(uint64_t, *(ez->srt), (ez->srt->n + ez->ol->length)); + ei = ez->srt->a; en = ez->srt->n; oi = ez->srt->a + ez->srt->n; on = ez->ol->length; + for (k = 0; k < on; k++) { + z = &(ez->ol->list[k]); z->is_match = z->strong = z->without_large_indel = 0; + oi[k] = z->y_id; oi[k] <<= 1; oi[k] |= z->y_pos_strand; + oi[k] <<= 32; oi[k] |= k; + } + + radix_sort_ec64(ei, ei + en); radix_sort_ec64(oi, oi + on); + nec = slash_overlap(ez->ol->list, ei, en, oi, on, in0, icn, in1, ez->exz, ez->qu, ez->tu, ez->rref, post_syn); + ///debug for memory + // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); + + if(on > nec) { + // gen_hc_r_alin_adv_adp_smp(ez, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 0); + if(asm_opt.simd_mm > 0) gen_hc_r_alin_adv_adp_smp_1(ez, 0); + else gen_hc_r_alin_adv_adp_smp_0(ez, 0); + } + // fprintf(stderr, "[M::%s] srt->n::%u, nec::%lu, on::%lu\n", __func__, (uint32_t)srt->n, nec, on); + ///debug for memory + // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); + } + + if(ez->ol->length) srt_olst(ez->ol); +} + void prt_ovlp_sam_0(char *cm, FILE *fp, char *ref_id, int32_t ref_id_n, char *qry_id, int32_t qry_id_n, char *qry_seq, uint64_t qry_seq_n, uint64_t rs, uint64_t re, uint64_t qs, uint64_t qe, uint64_t flag, uint64_t err0, bit_extz_t *ez) { @@ -4709,7 +4945,10 @@ static void worker_hap_ec(void *data, long i, int tid) overlap_region *aux_o = NULL/**, *rse_o = NULL**/, *rcc = NULL; asg64_v buf0; uint32_t qlen = 0, qw = 0; uint64_t tot_b = 0; double tt0 = 0, tt1 = 0; b->v8q.n = b->v8t.n = 0; set_ec_cov(asm_opt.het_cov, asm_opt.hom_cov, asm_opt.het_cov_set, asm_opt.polyploidy, het_a, hom_a); - + // if(i != 3799659) return; + // if(i !=3373007) return; + // if(i != 508213) return; + // if(i != 3177240) return; // if((i != 733166) && (i != 858708) && (i != 858732) && (i != 859819) && (i != 859899) && (i != 863486) && (i != 872165) && (i != 899887) && (i != 902298) && // (i != 906808) && (i != 946173) && (i != 952685) && (i != 983977) && (i != 1000227) && (i != 1011228) && (i != 1042858) && (i != 1045860) && (i != 1118558) && // (i != 1143886) && (i != 1155956) && (i != 1159490) && (i != 1179151) && (i != 1180199) && (i != 1230524) && (i != 1232338) && (i != 1244031) && (i != 1268467) && @@ -4738,7 +4977,7 @@ static void worker_hap_ec(void *data, long i, int tid) // if (memcmp("9edb4aa1-3a56-40dc-b687-77e4176d1053", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0 || // memcmp("25e30873-9737-4e6d-bdb6-a504a26ef3de", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0 || // memcmp("19699b82-2883-43e1-a11e-ec0c95eaccd4", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { - // fprintf(stderr, "\n+[M::%s]\trid-target::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + // fprintf(stderr, "\n+[M::%s]\trid-target::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); // } @@ -4806,15 +5045,12 @@ static void worker_hap_ec(void *data, long i, int tid) // fprintf(stderr, "\n+[M::%s]\trid::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); - ///r769: kp (gen_hc_r_alin_ea) -> NULL; site_sc (rphase_hc) -> 0 - ///r770: kp (gen_hc_r_alin_ea) -> buf0; site_sc (rphase_hc) -> 0 - ///r789: kp (gen_hc_r_alin_ea) -> NULL; site_sc (rphase_hc) -> 0 - ///r791: kp (gen_hc_r_alin_ea) -> buf0; site_sc (rphase_hc) -> 0 copy_asg_arr(buf0, b->sp); // tot_b = gen_hc_r_alin_ea_flt(b->ab, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_n_chain, asm_opt.max_n_chain*HC_MF_R, asm_opt.chn_occ, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, // 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1), (asm_opt.is_ont)?(512):(0), (asm_opt.is_ont)?(6):(0), (asm_opt.is_ont)?(1.5):(-1), (asm_opt.is_ont)?(0.1):(-1), &buf0, qw, &b->v8q, &b->v32, 1); - tot_b = gen_hc_r_alin_ea_flt_mmp(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_n_chain, asm_opt.max_n_chain*HC_MF_R, asm_opt.chn_occ, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, - 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1), (asm_opt.is_ont)?(512):(0), (asm_opt.is_ont)?(6):(0), (asm_opt.is_ont)?(1.5):(-1), (asm_opt.is_ont)?(0.1):(-1), &buf0, qw, &b->v8q, &b->v32); + ///gen_hc_r_alin_ea_flt_mmp -> r922 + tot_b = gen_hc_r_alin_ea_flt_mmp_adv(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_n_chain, asm_opt.max_n_chain*HC_MF_R, asm_opt.chn_occ, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, + 1, &b->v16, &b->v64, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1), (asm_opt.is_ont)?(512):(0), (asm_opt.is_ont)?(6):(0), (asm_opt.is_ont)?(1.5):(-1), (asm_opt.is_ont)?(0.1):(-1), &buf0, qw, &b->v8q, &b->v32, asm_opt.post_syn); copy_asg_arr(b->sp, buf0); // exit(1); @@ -4863,7 +5099,7 @@ static void worker_hap_ec(void *data, long i, int tid) gen_reseed_re(&b->olist, &b->clist, aux_o, rse_o, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, asm_opt.mz_win, 19, i, asm_opt.max_ov_diff_ec, asm_opt.max_ov_diff_ec, &b->v16, R_INF.tqn, b->v8q.a); copy_asg_arr(b->sp, buf0); **/ - if(scb.a[i].n) { + if(scb.a[i].n && asm_opt.realn_raw) { regen_scb(b->ab, &b->clist, i, &(scb.a[i]), &b->self_read, &b->ovlp_read, &b->v64, asm_opt.mz_win, asm_opt.k_mer_length, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, ((asm_opt.is_ont)?(0.05):(0.02)), 1, 1, &b->olist, &b->exz, ((asm_opt.max_ov_diff_ec>0.1)?(asm_opt.max_ov_diff_ec):(0.1)), (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), &b->v16, &aux_o, &rcc, &rl0); } @@ -4871,7 +5107,8 @@ static void worker_hap_ec(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); //site_sc: r765 -> r766: 1 -> 0 rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &(scb.a[i]), &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), ((asm_opt.is_sc)?&(b->v8t):NULL)/**&(b->v8t)**/, (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32, - asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0, rcc, rl0); + asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0, rcc, rl0, ((asm_opt.post_syn)?(5):(INT64_MAX)), + asm_opt.recurrent_err_normal_min, asm_opt.recurrent_err_normal_rat, asm_opt.recurrent_err_hpc_min, asm_opt.recurrent_err_hpc_rat, asm_opt.recurrent_err_test); copy_asg_arr(b->sp, buf0); ///for debug indel // stderr_phase_ovlp(&b->olist); @@ -4893,7 +5130,7 @@ static void worker_hap_ec(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), - R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL); + R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL, ((asm_opt.post_syn)?(0):(1))); copy_asg_arr(b->sp, buf0); if(DBG_TIME && dbg_a) { @@ -4905,7 +5142,7 @@ static void worker_hap_ec(void *data, long i, int tid) push_nec_re(aux_o, &(scc.a[i])); // push_nec_re(aux_o, &(scb.a[i])); - /**if(asm_opt.dbg_bam)**/ { + if(asm_opt.realn_raw/**asm_opt.dbg_bam**/) { update_scb(&R_INF, &(scc.a[i]), &(scb.a[i]), &(b->v16), &b->self_read, &b->ovlp_read, &b->v64, &b->exz, i); kv_resize(uint16_t, scb.a[i], b->v16.n); scb.a[i].n = b->v16.n; memcpy(scb.a[i].a, b->v16.a, b->v16.n*sizeof(*(scb.a[i].a))); } @@ -5116,10 +5353,6 @@ void worker_hap_ec_back_dbg(void *data, long i, int tid) // fprintf(stderr, "\n+[M::%s]\trid::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); - ///r769: kp (gen_hc_r_alin_ea) -> NULL; site_sc (rphase_hc) -> 0 - ///r770: kp (gen_hc_r_alin_ea) -> buf0; site_sc (rphase_hc) -> 0 - ///r789: kp (gen_hc_r_alin_ea) -> NULL; site_sc (rphase_hc) -> 0 - ///r791: kp (gen_hc_r_alin_ea) -> buf0; site_sc (rphase_hc) -> 0 /** copy_asg_arr(buf0, b->sp); // tot_b = gen_hc_r_alin_ea_flt(b->ab, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_n_chain, asm_opt.max_n_chain*HC_MF_R, asm_opt.chn_occ, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, @@ -5181,7 +5414,7 @@ void worker_hap_ec_back_dbg(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); //site_sc: r765 -> r766: 1 -> 0 rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32, - asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0, NULL, -1); + asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0, NULL, -1, ((asm_opt.post_syn)?(5):(INT64_MAX)), asm_opt.recurrent_err_normal_min, asm_opt.recurrent_err_normal_rat, asm_opt.recurrent_err_hpc_min, asm_opt.recurrent_err_hpc_rat, asm_opt.recurrent_err_test); copy_asg_arr(b->sp, buf0); ///for debug indel // stderr_phase_ovlp(&b->olist); @@ -5203,7 +5436,7 @@ void worker_hap_ec_back_dbg(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), - R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL); + R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL, ((asm_opt.post_syn)?(0):(1))); copy_asg_arr(b->sp, buf0); if(DBG_TIME && dbg_a) { @@ -5388,10 +5621,6 @@ static void worker_hap_ec_step(void *data, long i, int tid) // fprintf(stderr, "\n+[M::%s]\trid::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); - ///r769: kp (gen_hc_r_alin_ea) -> NULL; site_sc (rphase_hc) -> 0 - ///r770: kp (gen_hc_r_alin_ea) -> buf0; site_sc (rphase_hc) -> 0 - ///r789: kp (gen_hc_r_alin_ea) -> NULL; site_sc (rphase_hc) -> 0 - ///r791: kp (gen_hc_r_alin_ea) -> buf0; site_sc (rphase_hc) -> 0 copy_asg_arr(buf0, b->sp); tot_b = gen_hc_r_alin_ea_flt(b->ab, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_n_chain, asm_opt.max_n_chain*HC_MF_R, asm_opt.chn_occ, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1), (asm_opt.is_ont)?(512):(0), (asm_opt.is_ont)?(6):(0), @@ -5447,7 +5676,7 @@ static void worker_hap_ec_step(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); //site_sc: r765 -> r766: 1 -> 0 rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32, - asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0, NULL, -1); + asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0, NULL, -1, ((asm_opt.post_syn)?(5):(INT64_MAX)), asm_opt.recurrent_err_normal_min, asm_opt.recurrent_err_normal_rat, asm_opt.recurrent_err_hpc_min, asm_opt.recurrent_err_hpc_rat, asm_opt.recurrent_err_test); copy_asg_arr(b->sp, buf0); ///for debug indel // stderr_phase_ovlp(&b->olist); @@ -5466,7 +5695,7 @@ static void worker_hap_ec_step(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), - R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL); + R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL, ((asm_opt.post_syn)?(0):(1))); copy_asg_arr(b->sp, buf0); if(DBG_TIME && dbg_a) { @@ -5662,7 +5891,7 @@ static void worker_hap_ec_ss(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); //site_sc: r765 -> r766: 1 -> 0 rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), /**((asm_opt.is_sc)?&(b->v8t):NULL)**/&(b->v8t), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32, - asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0, NULL, -1); + asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0, NULL, -1, ((asm_opt.post_syn)?(5):(INT64_MAX)), asm_opt.recurrent_err_normal_min, asm_opt.recurrent_err_normal_rat, asm_opt.recurrent_err_hpc_min, asm_opt.recurrent_err_hpc_rat, asm_opt.recurrent_err_test); copy_asg_arr(b->sp, buf0); ///for debug indel // stderr_phase_ovlp(&b->olist); @@ -5674,7 +5903,7 @@ static void worker_hap_ec_ss(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), - R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL); + R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL, ((asm_opt.post_syn)?(0):(1))); copy_asg_arr(b->sp, buf0); @@ -5693,6 +5922,130 @@ static void worker_hap_ec_ss(void *data, long i, int tid) } +int64_t inline chk_sync_ovlp(ma_hit_t *p, uint64_t qn, uint64_t tn, double orate, uint8_t *rev_ov, ma_hit_t **rz) +{ + uint64_t k, tl; int64_t pq[2], pt[2], zq[2], zt[2], os, oe, ovlp; + ma_hit_t_alloc *pz = NULL; ma_hit_t *z; (*rev_ov) = (uint8_t)-1; *rz = NULL; + + pz = &(R_INF.paf[qn]); + for (k = 0; k < pz->length; k++) { + if((pz->buffer[k].tn != tn) || (pz->buffer[k].rev != p->rev)) continue; + z = &(pz->buffer[k]); + + // fprintf(stderr, "+[M::%s] qn::%lu, q::[%u,%u), %c, tn::%u, t::[%u,%u)\n", __func__, z->qns>>32, (uint32_t)z->qns, z->qe, "+-"[z->rev], z->tn, z->ts, z->te); + pq[0] = (uint32_t)p->qns; pq[1] = p->qe; + if(p->rev) { + tl = Get_READ_LENGTH(R_INF, p->tn); + pt[0] = (tl >=p->te)?(tl-p->te):(0); + pt[1] = (tl >=p->ts)?(tl-p->ts):(0); + } else { + pt[0] = p->ts; pt[1] = p->te; + } + + zt[0] = (uint32_t)z->qns; zt[1] = z->qe; + if(z->rev) { + tl = Get_READ_LENGTH(R_INF, z->tn); + zq[0] = (tl >=z->te)?(tl-z->te):(0); + zq[1] = (tl >=z->ts)?(tl-z->ts):(0); + } else { + zq[0] = z->ts; zq[1] = z->te; + } + + os = MAX(pq[0], zq[0]); oe = MIN(pq[1], zq[1]); + ovlp = ((oe>os)? (oe-os):0); + if(!((ovlp) && (ovlp >= ((pq[1] - pq[0])*orate)) && ((ovlp >= ((zq[1] - zq[0])*orate))))) continue; + + os = MAX(pt[0], zt[0]); oe = MIN(pt[1], zt[1]); + ovlp = ((oe>os)? (oe-os):0); + if(!((ovlp) && (ovlp >= ((pt[1] - pt[0])*orate)) && ((ovlp >= ((zt[1] - zt[0])*orate))))) continue; + (*rev_ov) = 0; *rz = z; + return k; + } + + pz = &(R_INF.reverse_paf[qn]); + for (k = 0; k < pz->length; k++) { + if((pz->buffer[k].tn != tn) || (pz->buffer[k].rev != p->rev)) continue; + // if(pz->buffer[k].no_l_indel == 0) continue; + z = &(pz->buffer[k]); + + // fprintf(stderr, "-[M::%s] qn::%lu, q::[%u,%u), %c, tn::%u, t::[%u,%u)\n", __func__, z->qns>>32, (uint32_t)z->qns, z->qe, "+-"[z->rev], z->tn, z->ts, z->te); + pq[0] = (uint32_t)p->qns; pq[1] = p->qe; + if(p->rev) { + tl = Get_READ_LENGTH(R_INF, p->tn); + pt[0] = (tl >=p->te)?(tl-p->te):(0); + pt[1] = (tl >=p->ts)?(tl-p->ts):(0); + } else { + pt[0] = p->ts; pt[1] = p->te; + } + + zt[0] = (uint32_t)z->qns; zt[1] = z->qe; + if(z->rev) { + tl = Get_READ_LENGTH(R_INF, z->tn); + zq[0] = (tl >=z->te)?(tl-z->te):(0); + zq[1] = (tl >=z->ts)?(tl-z->ts):(0); + } else { + zq[0] = z->ts; zq[1] = z->te; + } + + os = MAX(pq[0], zq[0]); oe = MIN(pq[1], zq[1]); + ovlp = ((oe>os)? (oe-os):0); + if(!((ovlp) && (ovlp >= ((pq[1] - pq[0])*orate)) && ((ovlp >= ((zq[1] - zq[0])*orate))))) continue; + + os = MAX(pt[0], zt[0]); oe = MIN(pt[1], zt[1]); + ovlp = ((oe>os)? (oe-os):0); + if(!((ovlp) && (ovlp >= ((pt[1] - pt[0])*orate)) && ((ovlp >= ((zt[1] - zt[0])*orate))))) continue; + (*rev_ov) = 1; *rz = z; + return k; + } + + return -1; +} + +static void select_sync_ovlp(void *data, long i, int tid) +{ + uint64_t k, qn = i; int64_t ok; uint8_t rev_ov; + ma_hit_t_alloc *pz = NULL; ma_hit_t *rz = NULL; + pz = &(R_INF.paf[i]); + for (k = 0; k < pz->length; k++) { + pz->buffer[k].bl = Get_READ_LENGTH(R_INF, pz->buffer[k].tn); + if(pz->buffer[k].tn == qn) continue; + // fprintf(stderr, "\n-a-[M::%s] qn::%lu, q::[%u,%u), %c, tn::%u, t::[%u,%u)\n", __func__, pz->buffer[k].qns>>32, (uint32_t)pz->buffer[k].qns, pz->buffer[k].qe, "+-"[pz->buffer[k].rev], + // pz->buffer[k].tn, pz->buffer[k].ts, pz->buffer[k].te); + ok = chk_sync_ovlp(&(pz->buffer[k]), pz->buffer[k].tn, qn, 0.666666, &rev_ov, &rz); + // fprintf(stderr, "[M::%s] ok::%ld, rev_ov::%u\n", __func__, ok, rev_ov); + + + if((ok == -1) || (rev_ov == 1 && rz->no_l_indel == 0)) { + pz->buffer[k].bl = (((uint64_t)1) << 30); + if((ok != -1) && (ok < 0x3fffffffu)) { + pz->buffer[k].bl |= ((uint64_t)ok); + } else { + pz->buffer[k].bl |= 0x3fffffffu; + } + } + } + + pz = &(R_INF.reverse_paf[i]); + for (k = 0; k < pz->length; k++) { + pz->buffer[k].bl = Get_READ_LENGTH(R_INF, pz->buffer[k].tn); + if(pz->buffer[k].tn == qn) continue; + if(pz->buffer[k].no_l_indel == 0) continue; + // fprintf(stderr, "\n-b-[M::%s] qn::%lu, q::[%u,%u), %c, tn::%u, t::[%u,%u)\n", __func__, pz->buffer[k].qns>>32, (uint32_t)pz->buffer[k].qns, pz->buffer[k].qe, "+-"[pz->buffer[k].rev], + // pz->buffer[k].tn, pz->buffer[k].ts, pz->buffer[k].te); + ok = chk_sync_ovlp(&(pz->buffer[k]), pz->buffer[k].tn, qn, 0.666666, &rev_ov, &rz); + // fprintf(stderr, "[M::%s] ok::%ld, rev_ov::%u\n", __func__, ok, rev_ov); + + if((ok == -1) || (rev_ov == 1 && rz->no_l_indel == 0)) { + pz->buffer[k].bl = (((uint64_t)1) << 30); + if((ok != -1) && (ok < 0x3fffffffu)) { + pz->buffer[k].bl |= ((uint64_t)ok); + } else { + pz->buffer[k].bl |= 0x3fffffffu; + } + } + } +} + static void worker_hap_ec_hybrid(void *data, long i, int tid) { ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); @@ -5710,6 +6063,13 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid) } b->v8q.n = b->v8t.n = 0; set_ec_cov(asm_opt.het_cov, asm_opt.hom_cov, asm_opt.het_cov_set, asm_opt.polyploidy, het_a, hom_a); + // if(i != 3799659) return; + // if(i != 3646295) return; + // if(i != 3373007) return; + // if(i != 508213) return; + // if(i != 3177240) return; + // e_h = e_l = 0.1;///this is for debug + // if(i != 10) return; // if((i%16) != 0) return; @@ -5723,7 +6083,7 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid) // if (memcmp("c7ecbd6b-e09d-4042-93ac-2400839feaf6", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { - // fprintf(stderr, "-a-[M::%s-beg] rid->%ld, rlen->%lu, scb.a[i].n::%u\n", __func__, i, Get_READ_LENGTH((R_INF),i), (uint32_t)scb.a[i].n); + // fprintf(stderr, "-a-[M::%s-beg] rid->%ld, rlen->%lu, scb.a[i].n::%u\n", __func__, i, Get_READ_LENGTH((R_INF),i), (uint32_t)scb.a[i].n); // } else { // return; // } @@ -5740,6 +6100,8 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid) h_ec_lchain_hybrid(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, bw_h, bw_l, /**((asm_opt.is_ont)?(0.05):(0.02)),**/ asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32, COV_W, R_INF.tqn, 0, 1);///ONT high error + // stderr_phase_ovlp(&b->olist); + // fprintf(stderr, "-b-[M::%s] rid::%ld\n", __func__, i); // b->num_read_base += b->olist.length; @@ -5756,7 +6118,7 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid) asm_opt.chn_occ, ((asm_opt.hom_cov*HC_AV_MIN)>0)?(asm_opt.hom_cov*HC_AV_MIN):(1), qw, R_INF.tqn, asm_opt.hom_cov); // gen_hc_r_alin_ea_adv(&ez); // gen_hc_r_alin_ea_adv_flt(&ez); - gen_hc_r_alin_ea_adv_flt_mmp(&ez); + gen_hc_r_alin_ea_adv_flt_mmp_adv(&ez, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), asm_opt.post_syn); copy_asg_arr(b->sp, buf0); // fprintf(stderr, "-c-[M::%s] rid::%ld\n", __func__, i); @@ -5795,21 +6157,15 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid) ///for debug indel // prt_ovlp_sam(&b->olist, &b->ovlp_read, b->self_read.seq, b->self_read.length); - if(scb.a[i].n) { + if(scb.a[i].n && asm_opt.realn_raw) { regen_scb(b->ab, &b->clist, i, &(scb.a[i]), &b->self_read, &b->ovlp_read, &b->v64, asm_opt.mz_win, asm_opt.k_mer_length, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, bw_h, 1, 1, &b->olist, &b->exz, ((e_h>0.1)?(e_h):(0.1)), (qid < R_INF.tqn)?(WINDOW_OHC):(WINDOW_HC), &b->v16, &aux_o, &rcc, &rl0); - // if(rcc) { - // fprintf(stderr, "[M::%s]\t%.*s(qid::%u)\tql::%ld\tq::[%u,\t%u)\t%c\t%.*s(tid::%u)\ttl::%ld\tt::[%u,\t%u)\terr::%u\n", __func__, - // (int32_t)Get_NAME_LENGTH(R_INF, rcc->x_id), Get_NAME(R_INF, rcc->x_id), rcc->x_id, (int64_t)b->self_read.length, rcc->x_pos_s, rcc->x_pos_e + 1, "+-"[rcc->y_pos_strand], - // (int32_t)Get_NAME_LENGTH(R_INF, rcc->y_id), Get_NAME(R_INF, rcc->y_id), rcc->y_id, rl0, rcc->y_pos_s, rcc->y_pos_e + 1, rcc->non_homopolymer_errors); - // } else { - // fprintf(stderr, "[M::%s]\tunalined\n", __func__); - // } } copy_asg_arr(buf0, b->sp); rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &(scb.a[i]), &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32, - asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, (((double)R_INF.tr[1])/((double)(R_INF.tr[0] + R_INF.tr[1]))), rcc, rl0); + asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, (((double)R_INF.tr[1])/((double)(R_INF.tr[0] + R_INF.tr[1]))), rcc, rl0, ((asm_opt.post_syn)?(5):(INT64_MAX)), + asm_opt.recurrent_err_normal_min, asm_opt.recurrent_err_normal_rat, asm_opt.recurrent_err_hpc_min, asm_opt.recurrent_err_hpc_rat, asm_opt.recurrent_err_test); copy_asg_arr(b->sp, buf0); ///for debug indel // stderr_phase_ovlp(&b->olist); @@ -5819,7 +6175,7 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); copy_asg_arr(buf1, b->hap.snp_srt); b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, R_INF.tqn, - R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, &buf1); + R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, &buf1, ((asm_opt.post_syn)?(0):(1))); copy_asg_arr(b->sp, buf0); copy_asg_arr(b->hap.snp_srt, buf1); // if(DBG_TIME && dbg_a) { @@ -5829,7 +6185,7 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid) push_nec_re(aux_o, &(scc.a[i])); // cmp_smp_ac(&b->self_read, &(scc.a[i]), &b->ovlp_read, i);///for debug // push_nec_re(aux_o, &(scb.a[i])); - /**if(asm_opt.dbg_bam)**/ { + if(asm_opt.realn_raw/**asm_opt.dbg_bam**/) { update_scb(&R_INF, &(scc.a[i]), &(scb.a[i]), &(b->v16), &b->self_read, &b->ovlp_read, &b->v64, &b->exz, i); kv_resize(uint16_t, scb.a[i], b->v16.n); memcpy(scb.a[i].a, b->v16.a, b->v16.n * sizeof((*(b->v16.a)))); scb.a[i].n = b->v16.n; // fprintf(stderr, "-b-[M::%s-beg] rid->%ld, rlen->%lu, scb.a[i].n::%u\n", __func__, i, Get_READ_LENGTH((R_INF),i), (uint32_t)scb.a[i].n); @@ -5987,7 +6343,7 @@ static void worker_hap_ec_hybrid_sync(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32, - asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, hf_rate, NULL, -1); + asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, hf_rate, NULL, -1, ((asm_opt.post_syn)?(5):(INT64_MAX)), asm_opt.recurrent_err_normal_min, asm_opt.recurrent_err_normal_rat, asm_opt.recurrent_err_hpc_min, asm_opt.recurrent_err_hpc_rat, asm_opt.recurrent_err_test); copy_asg_arr(b->sp, buf0); ///for debug indel // if(i == 23863) stderr_phase_ovlp(&b->olist); @@ -5996,7 +6352,7 @@ static void worker_hap_ec_hybrid_sync(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); copy_asg_arr(buf1, b->hap.snp_srt); b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, R_INF.tqn, - R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, hf_only?(NULL):(&buf1)); + R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, hf_only?(NULL):(&buf1), ((asm_opt.post_syn)?(0):(1))); copy_asg_arr(b->sp, buf0); copy_asg_arr(b->hap.snp_srt, buf1); push_nec_re(aux_o, &(scc.a[i])); @@ -7229,6 +7585,12 @@ uint32_t is_chemical_r_adv(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_ int64_t cal_chemical_r_adv(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_t cut_len, double dup_rate, uint64_t is_del) { + // if ((ov->length) && memcmp("573d5842-81dc-4c84-b28f-e1cd7bddb2b0", + // Get_NAME((R_INF), (ov->buffer[0].qns>>32)), Get_NAME_LENGTH((R_INF), (ov->buffer[0].qns>>32))) == 0) { + // fprintf(stderr, "-a-[M::%s-beg]\trid->%lu(%.*s)\trlen->%lu\tn_ov->%u\n", __func__, (ov->buffer[0].qns>>32), + // (int32_t)Get_NAME_LENGTH(R_INF, (ov->buffer[0].qns>>32)), Get_NAME(R_INF, (ov->buffer[0].qns>>32)), Get_READ_LENGTH((R_INF), (ov->buffer[0].qns>>32)), ov->length); + // } + uint64_t k, s, e; int64_t dp, old_dp, st = 0, ed, s0, e0, rr, lt, min_cov; for (k = idx->n = 0; k < ov->length; k++) { if(is_del && ov->buffer[k].del) continue; @@ -7240,12 +7602,20 @@ int64_t cal_chemical_r_adv(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_ lt = Get_READ_LENGTH((R_INF), ov->buffer[k].tn); rr = (lt >= len)?(lt - len):(len - lt); - if((rr <= (len*dup_rate)) && (rr <= (lt*dup_rate)) && (ov->buffer[k].rev)) { + if((rr <= (len*dup_rate)) && (rr <= (lt*dup_rate)) && (ov->buffer[k].rev)) {///this is for ONT-specific issue dp = (ov->buffer[k].qe) - ((uint32_t)ov->buffer[k].qns); dp = len - dp; old_dp = ov->buffer[k].te - ov->buffer[k].ts; old_dp = lt - old_dp; if((dp <= (len*dup_rate)) && (old_dp <= (lt*dup_rate))) continue; } + // if ((ov->length) && memcmp("573d5842-81dc-4c84-b28f-e1cd7bddb2b0", + // Get_NAME((R_INF), (ov->buffer[0].qns>>32)), Get_NAME_LENGTH((R_INF), (ov->buffer[0].qns>>32))) == 0) { + // fprintf(stderr, "[M::%s]\tqid::%lu(%.*s)\tqlen::%lu\tq::[%u,%u)\t%c\ttid::%u(%.*s)\ttlen::%lu\tt::[%u,%u)\n", __func__, + // (ov->buffer[k].qns>>32), (int32_t)Get_NAME_LENGTH(R_INF, (ov->buffer[k].qns>>32)), Get_NAME(R_INF, (ov->buffer[k].qns>>32)), Get_READ_LENGTH((R_INF), (ov->buffer[k].qns>>32)), + // (uint32_t)ov->buffer[k].qns, ov->buffer[k].qe, "+-"[ov->buffer[k].rev], + // ov->buffer[k].tn, (int32_t)Get_NAME_LENGTH(R_INF, ov->buffer[k].tn), Get_NAME(R_INF, ov->buffer[k].tn), Get_READ_LENGTH((R_INF), ov->buffer[k].tn), ov->buffer[k].ts, ov->buffer[k].te); + // } + kv_push(uint64_t, (*idx), (s<<1)); kv_push(uint64_t, (*idx), (e<<1)|1); } @@ -7259,7 +7629,7 @@ int64_t cal_chemical_r_adv(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_ ed = idx->a[k]>>1; if(ed > st) { - // if(ov->length && ((ov->buffer[0].qns>>32) == 5045637)) { + // if(ov->length && ((ov->buffer[0].qns>>32) == 3373007)) { // fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\tid::%lu\n", __func__, st, ed, old_dp, len, ov->buffer[0].qns>>32); // } if(old_dp <= min_cov) { @@ -7273,7 +7643,7 @@ int64_t cal_chemical_r_adv(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_ ed = len; old_dp = dp; if(ed > st) { - // if(ov->length && ((ov->buffer[0].qns>>32) == 5045637)) { + // if(ov->length && ((ov->buffer[0].qns>>32) == 3373007)) { // fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\tid::%lu\n", __func__, st, ed, old_dp, len, ov->buffer[0].qns>>32); // } if(old_dp <= min_cov) { @@ -7405,6 +7775,9 @@ static void worker_hap_dc_ec_chemical_arc_mark(void *data, long i, int tid) if(b->cnt[1] == 0) { msk[i] = (uint8_t)-1; cov = cal_chemical_r_adv(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), asm_opt.chemical_flank, 0.02, 1); + // if(i == 3373007/** || i == 3313150**/) { + // fprintf(stderr, "-um-[M::%s]\tqn::%u::%.*s\tcov::%ld\tpaf_n::%u\n\n", __func__, (uint32_t)(i), (int)Get_NAME_LENGTH(R_INF, i), Get_NAME((R_INF), i), cov, R_INF.paf[i].length); + // } if(cov <= msk_cut) msk[i] = cov; if(cov <= msk_cut/**FORCE_CUT**/) { // fprintf(stderr, "-um-[M::%s]\tqn::%u::%.*s\n\n", __func__, (uint32_t)(i), (int)Get_NAME_LENGTH(R_INF, i), Get_NAME((R_INF), i)); @@ -8971,12 +9344,12 @@ static void worker_hap_dc_ec0(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 1**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), /**((asm_opt.is_sc)?&(b->v8q):NULL)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32, - asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1, NULL, -1); + asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1, NULL, -1, ((asm_opt.post_syn)?(5):(INT64_MAX)), asm_opt.recurrent_err_normal_min, asm_opt.recurrent_err_normal_rat, asm_opt.recurrent_err_hpc_min, asm_opt.recurrent_err_hpc_rat, asm_opt.recurrent_err_test); copy_asg_arr(b->sp, buf0); copy_asg_arr(buf0, b->sp); b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), - R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL); + R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL, ((asm_opt.post_syn)?(0):(1))); copy_asg_arr(b->sp, buf0); push_nec_re(aux_o, &(scc.a[i])); @@ -9499,7 +9872,7 @@ static void *ff_ihyb_syn_worker_count(void *data, int step, void *in) tn = z->buffer[k].tn; if(tn < p->tqn) continue;///no ont-2-ont qn = z->buffer[k].qns>>32; - m = qn<<=32; m |= (k<<1); + m = qn<<32; m |= (k<<1); if(z->buffer[k].bl == 0x7FFFFFFF) { m |= 1; p->n_bl++; // z->buffer[k].bl = Get_READ_LENGTH(R_INF, tn); @@ -9508,7 +9881,7 @@ static void *ff_ihyb_syn_worker_count(void *data, int step, void *in) s->tov_size -= zi->m; kv_push(uint64_t, *zi, m); s->tov++; s->tov_size += zi->m; - p->n_ov++; + p->n_ov[0]++; // if(((z->buffer[k].qns>>32) == 23863 && z->buffer[k].tn == 4) || ((z->buffer[k].qns>>32) == 4 && z->buffer[k].tn == 23863)) { // fprintf(stderr, "[M::%s]\t%.*s(qid::%u)\tql::%lu\tq::[%u,%u)\t%c\t%.*s(tid::%u)\t\ttl::%lu\tt::[%u,%u)\tel::%u\n", __func__, // (int32_t)Get_NAME_LENGTH(R_INF, (z->buffer[k].qns>>32)), Get_NAME(R_INF, (z->buffer[k].qns>>32)), (uint32_t)(z->buffer[k].qns>>32), Get_READ_LENGTH(R_INF, (z->buffer[k].qns>>32)), @@ -9539,20 +9912,195 @@ static void *ff_ihyb_syn_worker_count(void *data, int step, void *in) return 0; } + +static void ff_ia_syn_worker_insert(void *data, long i, int tid) /** callback for kt_for()**/ +{ + tsrt_v_buf *s = ((tsrt_v_buf*)data); + asg64_v *za = &(s->idx[i]); int64_t n0; ma_hit_t_alloc *rr[2] = {R_INF.paf, R_INF.reverse_paf}; + uint64_t k, qn, tn, ok, ss, rp_k, ll; uint8_t is_replac; ma_hit_t *fv = NULL, *rv = NULL; ma_hit_t_alloc *rva = NULL; + if(s->tqn == 0) { + for (k = 0; k < za->n; k++) { + qn = za->a[k]>>32; ok = ((uint32_t)za->a[k])>>1; ss = za->a[k]&1; + fv = &(rr[ss][qn].buffer[ok]); tn = fv->tn; + assert((fv->bl>>30)); + if((fv->bl) < (0x7fffffffu)) continue;///replace, no increase + + rva = &(rr[ss][tn]); + if(rva->length >= rva->size) { + rva->length++; + } else { + rv = &(rva->buffer[rva->length++]); + rv->qns = (uint64_t)-1; + rv->tn = (uint32_t)-1; + } + } + } else if(s->tqn == 1) { + uint64_t m0 = i*s->pm, m1 = (i+1)*s->pm; if(m1 > s->tot) m1 = s->tot; + + for (k = m0, ss = 0; k < m1; k++) { + rva = &(rr[ss][k]); + + n0 = MIN(rva->length, rva->size); + if(rva->length > rva->size) { + rva->size = rva->length; + REALLOC(rva->buffer, rva->size); + } + + for (n0--; (n0 >= 0) && ((rva->buffer[n0].qns == ((uint64_t)-1)) || (rva->buffer[n0].tn == ((uint32_t)-1))); n0--); + rva->length = n0 + 1; + } + + for (k = m0, ss = 1; k < m1; k++) { + rva = &(rr[ss][k]); + + n0 = MIN(rva->length, rva->size); + if(rva->length > rva->size) { + rva->size = rva->length; + REALLOC(rva->buffer, rva->size); + } + + for (n0--; (n0 >= 0) && ((rva->buffer[n0].qns == ((uint64_t)-1)) || (rva->buffer[n0].tn == ((uint32_t)-1))); n0--); + rva->length = n0 + 1; + } + } else { + for (k = 0; k < za->n; k++) { + qn = za->a[k]>>32; ok = ((uint32_t)za->a[k])>>1; ss = za->a[k]&1; + fv = &(rr[ss][qn].buffer[ok]); tn = fv->tn; rp_k = (uint64_t)-1; + assert((fv->bl>>30)); + is_replac = ((fv->bl) < (0x7fffffffu))?1:0; + + if(is_replac) { + rva = &(rr[1][tn]); rp_k = (fv->bl)&(0x3fffffffu); + assert(rp_k < rva->length && rp_k < rva->size); + } else { + rva = &(rr[ss][tn]); + } + fv->bl = Get_READ_LENGTH(R_INF, fv->tn); + + if(rp_k == ((uint64_t)-1)) { + rv = &(rva->buffer[rva->length++]); + rv->qns = Get_tn(*fv); + rv->qns = rv->qns << 32; + rv->tn = Get_qn(*fv); + rv->rev = fv->rev; + rv->el = fv->el; + rv->ml = fv->ml; + + if(fv->rev == 0) { + rv->qns = rv->qns | Get_ts(*fv); + rv->qe = Get_te(*fv); + rv->ts = Get_qs(*fv); + rv->te = Get_qe(*fv); + } else { + ll = Get_READ_LENGTH(R_INF, Get_tn(*fv)); + rv->qns |= ((ll>=Get_te(*fv))?(ll-Get_te(*fv)):(0)); + rv->qe = ((ll>=Get_ts(*fv))?(ll-Get_ts(*fv)):(0)); + + ll = Get_READ_LENGTH(R_INF, Get_qn(*fv)); + rv->ts = ((ll>=Get_qe(*fv))?(ll-Get_qe(*fv)):(0)); + rv->te = ((ll>=Get_qs(*fv))?(ll-Get_qs(*fv)):(0)); + } + } else { + rv = &(rva->buffer[rp_k]); + assert(Get_qn(*rv) == Get_tn(*fv)); + assert(Get_tn(*rv) == Get_qn(*fv)); + assert(rv->no_l_indel == 0); + } + + rv->no_l_indel = 1; + rv->bl = Get_READ_LENGTH(R_INF, rv->tn); + rv->del = 0; + rv->cc = 0x3fffffffu; + } + } +} + + +static void ia_syn_worker_count(tsrt_v_m *p) +{ + ma_hit_t_alloc *z = NULL; uint64_t k, qn, tn, m, mm = (((uint64_t)1)<<30); asg64_v *zi = NULL; + tsrt_v_buf *s = NULL; CALLOC(s, 1); + s->p = p->p; s->pm = p->pm; s->pn = p->pn; s->tot = p->tot; s->tov = s->tov_size = 0; CALLOC(s->idx, s->pn); + + while (p->rid < p->tot) { + z = &(R_INF.paf[p->rid]); qn = p->rid; + for (k = 0; k < z->length; k++) { + if((z->buffer[k].bl&mm) == 0) continue; + if(k > 0x7fffffffu) { + z->buffer[k].bl = Get_READ_LENGTH(R_INF, z->buffer[k].tn); + continue; + } + tn = z->buffer[k].tn; + m = qn<<32; m |= (k<<1); + zi = &(s->idx[tn/s->pm]); + kv_push(uint64_t, *zi, m); s->tov++; s->tov_size++; + if((z->buffer[k].bl) < (0x7fffffffu)) p->n_bl++; + else p->n_ov[0]++; + } + z = &(R_INF.reverse_paf[p->rid]); + for (k = 0; k < z->length; k++) { + if(z->buffer[k].no_l_indel == 0) continue; + if((z->buffer[k].bl&mm) == 0) continue; + if(k > 0x7fffffffu) { + z->buffer[k].bl = Get_READ_LENGTH(R_INF, z->buffer[k].tn); + continue; + } + tn = z->buffer[k].tn; + m = qn<<32; m |= (k<<1) + 1; + zi = &(s->idx[tn/s->pm]); + kv_push(uint64_t, *zi, m); s->tov++; s->tov_size++; + if((z->buffer[k].bl) < (0x7fffffffu)) p->n_bl++; + else p->n_ov[1]++; + } + p->rid++; + if(s->tov_size >= p->chunk_size) { + if(s->tov) { + s->tqn = 0; + kt_for(p->n_thr, ff_ia_syn_worker_insert, s, s->pn); + s->tqn = 1; + kt_for(p->n_thr, ff_ia_syn_worker_insert, s, s->pn); + s->tqn = 2; + kt_for(p->n_thr, ff_ia_syn_worker_insert, s, s->pn); + for (k = 0; k < s->pn; k++) { + s->idx[k].n = 0; + } + } + s->tov = s->tov_size = 0; + } + } + + if(s->tov) { + s->tqn = 0; + kt_for(p->n_thr, ff_ia_syn_worker_insert, s, s->pn); + s->tqn = 1; + kt_for(p->n_thr, ff_ia_syn_worker_insert, s, s->pn); + s->tqn = 2; + kt_for(p->n_thr, ff_ia_syn_worker_insert, s, s->pn); + for (k = 0; k < s->pn; k++) { + s->idx[k].n = 0; + } + } + + for (k = 0; k < s->pn; k++) { + free(s->idx[k].a); + } + free(s->idx); free(s); +} + void ff_ihyb_syn_tid(ec_ovec_buf_t *b, uint64_t pre, uint64_t n_a, uint64_t n_thre) { tsrt_v_m sp = {0, 0, 0, 0, 0, 0}; - sp.p = pre; sp.pn = ((uint64_t)1) << pre; sp.pm = (((uint64_t)1) << pre) - 1; sp.n_ov = sp.n_bl = 0; + sp.p = pre; sp.pn = ((uint64_t)1) << pre; sp.pm = (((uint64_t)1) << pre) - 1; sp.n_ov[0] = sp.n_ov[1] = sp.n_bl = 0; sp.rid = 0; sp.tot = n_a; sp.chunk_size = 10000000; sp.tqn = R_INF.tqn; sp.n_thr = n_thre; - sp.pf = R_INF.paf; sp.n_ov = sp.n_bl = 0; sp.rid = 0; + sp.pf = R_INF.paf; sp.n_ov[0] = sp.n_ov[1] = sp.n_bl = 0; sp.rid = 0; kt_pipeline(n_thre, ff_ihyb_syn_worker_count, &sp, 2); - fprintf(stderr, "[M::%s::cis-paf] # syn overlaps::%lu, # syn informative overlaps::%lu\n", __func__, sp.n_ov, sp.n_bl); + fprintf(stderr, "[M::%s::cis-paf] # syn overlaps::%lu, # syn informative overlaps::%lu\n", __func__, sp.n_ov[0], sp.n_bl); - sp.pf = R_INF.reverse_paf; sp.n_ov = sp.n_bl = 0; sp.rid = 0; + sp.pf = R_INF.reverse_paf; sp.n_ov[0] = sp.n_ov[1] = sp.n_bl = 0; sp.rid = 0; kt_pipeline(n_thre, ff_ihyb_syn_worker_count, &sp, 2); - fprintf(stderr, "[M::%s::trans-paf] # syn overlaps::%lu, # syn informative overlaps::%lu\n", __func__, sp.n_ov, sp.n_bl); + fprintf(stderr, "[M::%s::trans-paf] # syn overlaps::%lu, # syn informative overlaps::%lu\n", __func__, sp.n_ov[0], sp.n_bl); kt_for(n_thre, worker_hap_ec_hybrid_sync, b, n_a - R_INF.tqn);///HiFi-only @@ -9572,7 +10120,7 @@ void ff_ihyb_syn_tid(ec_ovec_buf_t *b, uint64_t pre, uint64_t n_a, uint64_t n_th for (i = 0; i < z->length; i++) { tn = z->buffer[i].tn; qn = z->buffer[i].qns>>32; - m = qn<<=32; m |= (i<<1); + m = qn<<32; m |= (i<<1); if(z->buffer[i].bl == 0x7FFFFFFF) { m |= 1; z->buffer[i].bl = Get_READ_LENGTH(R_INF, tn); @@ -9605,6 +10153,24 @@ void gen_ihyb_syn(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a) } +void ff_ia_syn_tid(uint64_t pre, uint64_t n_a, uint64_t n_thre) +{ + double tt0 = yak_realtime_0(); + kt_for(n_thre, select_sync_ovlp, NULL, n_a); + // exit(1); + + tsrt_v_m sp = {0, 0, 0, 0, 0, 0}; + sp.p = pre; sp.pn = ((uint64_t)1) << pre; sp.pm = (n_a+sp.pn-1)/(sp.pn)/**(((uint64_t)1) << pre) - 1**/; sp.n_ov[0] = sp.n_ov[1] = sp.n_bl = 0; + sp.rid = 0; sp.tot = n_a; sp.chunk_size = 10000000; sp.tqn = n_a; sp.n_thr = n_thre; + + sp.pf = NULL; sp.n_ov[0] = sp.n_ov[1] = sp.n_bl = 0; sp.rid = 0; + // kt_pipeline(n_thre, ff_ia_syn_worker_count, &sp, 2); + ia_syn_worker_count(&sp); + fprintf(stderr, "[M::%s::cis-paf::%.3f] # syn overlaps(+)::%lu, # syn overlaps(-)::%lu, # replace overlaps::%lu\n", __func__, yak_realtime_0()-tt0, sp.n_ov[0], sp.n_ov[1], sp.n_bl); + + // kt_for(n_thre, worker_hap_ec_hybrid_sync, b, n_a - R_INF.tqn);///HiFi-only +} + uint64_t cal_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, uint64_t *r_base) { double tt0 = yak_realtime_0(); @@ -9957,7 +10523,7 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u cal_update_ec_multiple(b, n_thre, n_a);///update overlaps // if(is_sv) kt_for(n_thre, worker_hap_dc_ec, b, n_a);///update overlaps - fprintf(stderr, "-2-[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]); + // fprintf(stderr, "-2-[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]); // prt_nel_ovlp(R_INF.paf, n_a); // exit(1); @@ -9968,15 +10534,16 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u kt_for(n_thre, worker_hap_post_rev, b, n_a); } - fprintf(stderr, "-3-[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]); + // fprintf(stderr, "-3-[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]); // cal_sec_ec_multiple(b, n_thre, n_a, -1); // gen_sec_ec_multiple(b, n_thre, n_a); destroy_ec_ovec_buf_t(b); - fprintf(stderr, "-4-[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]); + // fprintf(stderr, "-4-[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]); + if((!is_sv) && (asm_opt.post_syn)) ff_ia_syn_tid(10, n_a, n_thre); // dbg_write_ec_reads("ec16.fa", round, &scb, 0/**!is_cr**/); // exit(1); diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 495c4fe..b663cfb 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -411,6 +411,153 @@ void init_integer_ml_t(integer_ml_t *x, ul_resolve_t *u, uint64_t n_thread) int32_t if_sup_chimeric(ma_hit_t_alloc* src, uint64_t rLen, asg64_v *b, int if_exact); +typedef struct { + const char *nn; + uint64_t nid; +} dbg_prt_t; + +static void qry_nn_id(void *data, long i, int tid) // callback for kt_for() +{ + dbg_prt_t *sl = (dbg_prt_t *)data; uint64_t ri = sl[i].nid, rl = 0; const char *nn = sl[i].nn; + rl = strlen(nn); + if((ri < R_INF.total_reads) && (Get_NAME_LENGTH((R_INF), ri) == rl) && (memcmp(nn, Get_NAME((R_INF), ri), Get_NAME_LENGTH((R_INF), ri)) == 0)) return; + + for (ri = 0; ri < R_INF.total_reads; ri++) { + if((Get_NAME_LENGTH((R_INF), ri) == rl) && (memcmp(nn, Get_NAME((R_INF), ri), Get_NAME_LENGTH((R_INF), ri)) == 0)) { + sl[i].nid = ri; + return; + } + } + sl[i].nid = (uint64_t)-1; +} + +void hc_dbg_prt_ma_hit_t(const char *ids0[], uint64_t n0, const uint64_t iids0[], uint64_t ni0, const char* cmd, ma_hit_t_alloc *paf, R_to_U* ruIndex, ma_sub_t* cov, asg_t* g) +{ + if((!paf) && (!ruIndex) && (!cov) && (!g)) return; + uint64_t k, z, n, ni, rid, uri, n_thre = asm_opt.thread_num; dbg_prt_t *ri = NULL; const char **ids = NULL; const uint64_t *iids = NULL; + if (cmd == NULL) cmd = ""; + fprintf(stderr, "\n[M::%s::%s]\t********************************\n", __func__, cmd); + const char *ids1[] = { + ///"m84039_230928_213653_s3/181277190/ccs",///tip 0, contained read issue + + ///"7e61bf83-7f09-49df-88b9-2a864587b19f", ///tip1 -> 916fced8-524e-4a65-b04d-5b5ea6937e6e(qn::199123) + ///phasing error -> ONT specific bias && HPC intoduce fake SNPs + + ///"43691917-c944-4102-8a01-d63a9c635dc3",///tip3 -> 29ee66bc-ef49-42ae-b7b4-15c37d46beac(qn::3373007) overcorrected in the first round + //when only 2 reads supported (the read itself has high sequencing error rate so hard to be aligned) + //we may also consider increase the error rate threshold + + ///"21b83ce9-9b94-445d-93b0-b9e83aaa4326",///tip3 + + "6f240062-2451-4212-a22d-abecc73ae0d9",///tip4->63a321ea-1103-4ca5-97e1-531e55f6bc84(qn::3646295) + ///phasing error -> ONT specific bias && HPC intoduce fake SNPs && the DP was to nice for the following SNPs, which we need to fix as high-sequencing depth we have more such cases + // +[M::gen_rphase_dp0_single_path_hybrid_0_multi] rn::2 hf_only::0 + // +[M::gen_rphase_dp0_single_path_hybrid_0_multi] site::175009 sc::2 n0::57 n1::3 rn::2 krn::2 b0l::16 b0h::41 b1l::0 b1h::3 + // +[M::gen_rphase_dp0_single_path_hybrid_0_multi] site::82558 sc::-1 n0::80 n1::4 rn::2 krn::2 b0l::4 b0h::76 b1l::2 b1h::2 + + "a72d3f7f-d74f-48ab-a3f7-5ecc74ad08ff",///tip5, contained read issue + }; + if(ids0) { + ids = ids0; n = n0; + } else { + ids = ids1; n = sizeof(ids1) / sizeof(ids1[0]); + } + + const uint64_t iids1[] = { + ///5038496,///tip 0, contained read issue + ///3799659,///tip 1, contained read issue + ///3313150,///tip 2, contained read issue + ///4192664,///tip3, contained read issue + 508213,///tip4 + 3203512,///tip5 + }; + if(iids0) { + iids = iids0; ni = ni0; + } else { + iids = iids1; ni = sizeof(iids1) / sizeof(iids1[0]); + } + CALLOC(ri, n); + + for (k = uri = 0; k < n; k++) { + rid = ((k= R_INF.total_reads) uri++; + ri[k].nid = rid; ri[k].nn = ids[k]; + } + if(uri) { + if(n_thre > n) n_thre = n; + kt_for(n_thre, qry_nn_id, ri, n); + } + + ma_hit_t *h = NULL; uint64_t nv, v, w; asg_arc_t *av; uint32_t cid, is_u; + for (k = 0; k < n; k++) { + if(ri[k].nid != ((uint64_t)-1)) { + fprintf(stderr, "\n[M::%s::%s]\trid::%lu\t%.*s\trlen::%lu\n", __func__, cmd, ri[k].nid, (int)Get_NAME_LENGTH(R_INF, ri[k].nid), Get_NAME(R_INF, ri[k].nid), Get_READ_LENGTH(R_INF, ri[k].nid)); + if(paf) { + fprintf(stderr, "[M::%s::%s]\tpaf::(abnm->%u\tfcor->%u)\n", __func__, cmd, paf[ri[k].nid].is_abnormal, paf[ri[k].nid].is_fully_corrected); + for (z = 0; z < paf[ri[k].nid].length; z++) { + h = &(paf[ri[k].nid].buffer[z]); + fprintf(stderr, "[%s::paf]\t%.*s(qn::%u)\t%u\t%u\t%u\t%c\t%.*s(tn::%u)\t%u\t%u\t%u\tml::%u\tbl::%u\tel::%u\tdel::%u\tnl_idl::%u", cmd, (int)Get_NAME_LENGTH(R_INF, Get_qn(*h)), Get_NAME((R_INF), Get_qn(*h)), Get_qn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_qn(*h)), Get_qs(*h), Get_qe(*h), "+-"[h->rev], + (int)Get_NAME_LENGTH(R_INF, Get_tn(*h)), Get_NAME((R_INF), Get_tn(*h)), Get_tn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_tn(*h)), Get_ts(*h), Get_te(*h), h->ml, h->bl, + h->el, h->del, h->no_l_indel); + if(ruIndex) { + get_R_to_U(ruIndex, Get_tn(*h), &cid, &is_u); + if((is_u != (uint32_t)-1) && (is_u == 0)) { + fprintf(stderr, "\tis_contain::1\tis_cr::1(crid->::%u::%.*s)", cid, (int)Get_NAME_LENGTH(R_INF, cid), Get_NAME(R_INF, cid)); + } + } + fprintf(stderr, "\n"); + } + } + if(g) { + fprintf(stderr, "[M::%s::%s]\tsg::(del->%u\tc->%u)\n", __func__, cmd, g->seq[ri[k].nid].del, g->seq[ri[k].nid].c); + v = ri[k].nid<<1; + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + for (z = 0; z < nv; ++z) { + w = av[z].v; + fprintf(stderr, "[M::%s::sg]\t%.*s(%c)\t%.*s(%c)\tol::%u\tel::%u\tdel::%u\tstrong::%u\tnl_idl::%u\n", __func__, + (int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), "+-"[v&1], + (int)Get_NAME_LENGTH(R_INF, (w>>1)), Get_NAME(R_INF, (w>>1)), "+-"[w&1], av[z].ol, av[z].el, av[z].del, av[z].strong, av[z].no_l_indel); + } + + v = (ri[k].nid<<1) + 1; + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + for (z = 0; z < nv; ++z) { + w = av[z].v; + fprintf(stderr, "[M::%s::sg]\t%.*s(%c)\t%.*s(%c)\tol::%u\tel::%u\tdel::%u\tstrong::%u\tnl_idl::%u\n", __func__, + (int)Get_NAME_LENGTH(R_INF, (v>>1)), Get_NAME(R_INF, (v>>1)), "+-"[v&1], + (int)Get_NAME_LENGTH(R_INF, (w>>1)), Get_NAME(R_INF, (w>>1)), "+-"[w&1], av[z].ol, av[z].el, av[z].del, av[z].strong, av[z].no_l_indel); + } + } + + if(cov) { + fprintf(stderr, "[M::%s::%s]\tcov::(del->%u\tc->%u)\n", __func__, cmd, cov[ri[k].nid].del, cov[ri[k].nid].c); + } + + if(ruIndex) { + fprintf(stderr, "[M::%s::%s]\t", __func__, cmd); + get_R_to_U(ruIndex, ri[k].nid, &cid, &is_u); + if((is_u != (uint32_t)-1) && (is_u == 0)) { + fprintf(stderr, "is_contain::1\tis_cr::1(rid->::%u::%.*s)\tis_cu::0\n", cid, (int)Get_NAME_LENGTH(R_INF, cid), Get_NAME(R_INF, cid)); + } else if(is_u != (uint32_t)-1) { + fprintf(stderr, "is_contain::1\tis_cr::0\tis_cu::1(uid->::%u)\n", cid); + } else { + fprintf(stderr, "is_contain::0\tis_cr::0\tis_cu::0\n"); + } + } + + } else { + fprintf(stderr, "\n[M::%s::%s]\trid::-1\t%.*s\n", __func__, cmd, (int)strlen(ri[k].nn), ri[k].nn); + } + } + + free(ri); +} + void print_edge(asg_arc_t *t, const char *cmd) { uint32_t v = t->ul>>32, w = t->v; @@ -3057,41 +3204,53 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i if(asm_opt.is_ont) asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_a", NULL, rI, NULL, sg); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);///p_telo // fprintf(stderr, "[M::%s] count_edges_v_w(sg, 49778, 49847)->%ld\n", __func__, count_edges_v_w(sg, 49778, 49847)); // if(is_ou) dedup_contain_g(uopt, sg); // debug_info_of_specfic_node("c7ecbd6b-e09d-4042-93ac-2400839feaf6", sg, rI, "beg-1"); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_b", NULL, rI, NULL, sg); for (i = 0; i < clean_round; i++, drop += step) { if(drop > max_ovlp_drop_ratio) drop = max_ovlp_drop_ratio; if(is_ou) { if(drop <= 0.500001) min_diff = step_diff>>1; else min_diff = step_diff; } + // fprintf(stderr, "\n(0):i->%ld, drop->%f\n", i, drop); if(asm_opt.is_ont) { asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_0", NULL, rI, NULL, sg); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_1", NULL, rI, NULL, sg); } - // fprintf(stderr, "(0):i->%ld, drop->%f\n", i, drop); // prt_specfic_sge(sg, 22708, 22646, "--0--"); // print_vw_edge(sg, 34156, 34090, "0"); // stats_chimeric(sg, src, &bu); if(!is_ou) asg_iterative_semi_circ(sg, src, &bu, max_tip, 1, uopt->te);///p_telo + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_3", NULL, rI, NULL, sg); asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); asg_arc_cut_chimeric(sg, src, &bu, is_ou?ou_thres:(uint32_t)-1, uopt->te);///p_telo + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_4", NULL, rI, NULL, sg); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_5", NULL, rI, NULL, sg); // prt_specfic_sge(sg, 22708, 22646, "--1--"); asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); asg_arc_cut_inexact(sg, src, &bu, max_tip, is_ou, is_trio, min_diff, ou_drop_rate/**, NULL**//**&dbg**/); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_6", NULL, rI, NULL, sg); // debug_edges(&dbg, d, 2); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_7", NULL, rI, NULL, sg); // prt_specfic_sge(sg, 22708, 22646, "--2--"); asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); asg_arc_cut_length(sg, &bu, max_tip, drop, ou_drop_rate, is_ou, is_trio, 1, min_diff, 1, NULL, NULL, NULL); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_8", NULL, rI, NULL, sg); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_9", NULL, rI, NULL, sg); // if(i == 0) { // print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty3.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0); @@ -3107,11 +3266,14 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); asg_arc_cut_bub_links(sg, &bu, HARD_OL_DROP, HARD_OL_SEC_DROP, HARD_OU_DROP, is_ou, asm_opt.large_pop_bubble_size, rev, rI, max_tip); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_10", NULL, rI, NULL, sg); // prt_specfic_sge(sg, 22708, 22646, "--4--"); asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); asg_arc_cut_complex_bub_links(sg, &bu, HARD_OL_DROP, HARD_OU_DROP, is_ou, b_mask_t); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_11", NULL, rI, NULL, sg); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_c_12", NULL, rI, NULL, sg); // prt_specfic_sge(sg, 22708, 22646, "--5--"); /** @@ -3133,7 +3295,9 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i if(asm_opt.is_ont) { asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_0", NULL, rI, NULL, sg); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_1", NULL, rI, NULL, sg); } if(is_ou) min_diff = step_diff; @@ -3147,6 +3311,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i if(clean_contain_g(uopt, sg, 1)) update_sg_uo(sg, src); } if(!is_ou) asg_iterative_semi_circ(sg, src, &bu, max_tip, 1, uopt->te); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_2", NULL, rI, NULL, sg); // prt_specfic_sge(sg, 22708, 22646, "--sb-0---"); @@ -3155,23 +3320,29 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); asg_cut_large_indel(sg, &bu, max_tip, HARD_OU_DROP, is_ou, min_diff);///shoule we ignore ou here? + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_3", NULL, rI, NULL, sg); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_4", NULL, rI, NULL, sg); // prt_specfic_sge(sg, 22708, 22646, "--sb-1---"); - // print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty6.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0); + // print_debug_gfa(sg, NULL, uopt->coverage_cut, "hybrid.dirty6.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0); if(!is_ou) { ///asg_arc_del_triangular_directly might be unnecessary asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); asg_arc_cut_length(sg, &bu, max_tip, HARD_ORTHOLOGY_DROP/**min_ovlp_drop_ratio**/, ou_drop_rate, is_ou, 0/**is_trio**/, is_ou?1:0, min_diff, 1, rev, rI, NULL); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_5", NULL, rI, NULL, sg); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_6", NULL, rI, NULL, sg); // print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty7.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0); asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); asg_arc_cut_length(sg, &bu, max_tip, min_ovlp_drop_ratio, ou_drop_rate, is_ou, 0/**is_trio**/, is_ou?1:0, min_diff, 1, rev, rI, &l_drop); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_7", NULL, rI, NULL, sg); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_8", NULL, rI, NULL, sg); } else { min_diff = step_diff; l_drop = 6000; asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); @@ -3184,6 +3355,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i // print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty8.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0); if(!is_ou) asg_cut_semi_circ(sg, LIM_LEN, 1); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_9", NULL, rI, NULL, sg); // prt_specfic_sge(sg, 22708, 22646, "--sb-3---"); /** @@ -3198,10 +3370,12 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i rescue_bubble_by_chain(sg, uopt->coverage_cut, src, rev, (asm_opt.max_short_tip*2), 0.15, 3, rI, 0.05, 0.9, uopt->max_hang, uopt->min_ovlp, 10, uopt->gap_fuzz, b_mask_t); **/ post_rescue(uopt, sg, src, rev, rI, b_mask_t, is_ou, cmk); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_10", NULL, rI, NULL, sg); // prt_specfic_sge(sg, 22708, 22646, "--sb-4---"); ug_ext_gfa(uopt, sg, ug_ext_len); + // hc_dbg_prt_ma_hit_t(NULL, 0, NULL, 0, "clen_d_11", NULL, rI, NULL, sg); // if(is_ou) dedup_contain_g(uopt, sg); diff --git a/gfa_ut.h b/gfa_ut.h index 57c6f68..e1b78f4 100644 --- a/gfa_ut.h +++ b/gfa_ut.h @@ -41,5 +41,6 @@ void ug_ext_gfa(ug_opt_t *uopt, asg_t *sg, uint32_t max_len); void update_sg_uo(asg_t *g, ma_hit_t_alloc *src); uint32_t get_arcs(asg_t *g, uint32_t v, uint32_t* idx, uint32_t idx_n); uint64_t ug_occ_w(uint64_t is, uint64_t ie, ma_utg_t *u); +void hc_dbg_prt_ma_hit_t(const char *ids0[], uint64_t n0, const uint64_t iids0[], uint64_t ni0, const char* cmd, ma_hit_t_alloc *paf, R_to_U* ruIndex, ma_sub_t* cov, asg_t* g); #endif diff --git a/htab.cpp b/htab.cpp index 6061cfc..8c1a705 100644 --- a/htab.cpp +++ b/htab.cpp @@ -1420,7 +1420,7 @@ int load_ct_index(void **i_ct_idx, char* file_name) return 1; } -int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name) +int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name, uint8_t force_rpaf_load) { char* gfa_name = (char*)malloc(strlen(file_name)+64); if(r) sprintf(gfa_name, "%s.pt_flt", file_name); @@ -1479,6 +1479,19 @@ int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* fwrite(&(r->paf[k].length), sizeof(r->paf[k].length), 1, fp); fwrite(r->paf[k].buffer, sizeof((*(r->paf[k].buffer))), r->paf[k].length, fp); } + + + int8_t rff = opt->post_syn; + if(force_rpaf_load) rff = 1; + fwrite(&rff, sizeof(rff), 1, fp); + if(rff) { + for (k = 0; k < r->total_reads; k++) { + fwrite(&(r->reverse_paf[k].is_fully_corrected), sizeof(r->reverse_paf[k].is_fully_corrected), 1, fp); + fwrite(&(r->reverse_paf[k].is_abnormal), sizeof(r->reverse_paf[k].is_abnormal), 1, fp); + fwrite(&(r->reverse_paf[k].length), sizeof(r->reverse_paf[k].length), 1, fp); + fwrite(r->reverse_paf[k].buffer, sizeof((*(r->reverse_paf[k].buffer))), r->reverse_paf[k].length, fp); + } + } } fprintf(stderr, "[M::%s] Index has been written.\n", __func__); @@ -1487,7 +1500,7 @@ int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* return 1; } -int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_opt_t* opt, char* file_name) +int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_opt_t* opt, char* file_name, uint8_t force_rpaf_load) { char* gfa_name = (char*)malloc(strlen(file_name)+64); if(r) sprintf(gfa_name, "%s.pt_flt", file_name); @@ -1604,6 +1617,27 @@ int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, hifiasm_op r->paf[k].buffer = (ma_hit_t*)malloc(sizeof(ma_hit_t)*r->paf[k].length); fread(r->paf[k].buffer, sizeof((*(r->paf[k].buffer))), r->paf[k].length, fp); } + + int8_t rff = 0; + f_flag += fread(&rff, sizeof(rff), 1, fp); + if(!force_rpaf_load) { + opt->post_syn = rff; + fprintf(stderr, "[M::%s] overwritten post_syn to::%u\n", __func__, (uint32_t)opt->post_syn); + } + if(rff) { + for (k = 0; k < r->total_reads; k++) { + f_flag += fread(&(r->reverse_paf[k].is_fully_corrected), sizeof(r->reverse_paf[k].is_fully_corrected), 1, fp); + f_flag += fread(&(r->reverse_paf[k].is_abnormal), sizeof(r->reverse_paf[k].is_abnormal), 1, fp); + f_flag += fread(&(r->reverse_paf[k].length), sizeof(r->reverse_paf[k].length), 1, fp); + r->reverse_paf[k].size = r->reverse_paf[k].length; + + r->reverse_paf[k].buffer = NULL; + if(r->reverse_paf[k].length == 0) continue; + + r->reverse_paf[k].buffer = (ma_hit_t*)malloc(sizeof(ma_hit_t)*r->reverse_paf[k].length); + fread(r->reverse_paf[k].buffer, sizeof((*(r->reverse_paf[k].buffer))), r->reverse_paf[k].length, fp); + } + } } fclose(fp); @@ -1619,16 +1653,16 @@ void refresh_pt_idx(void **flt_tab, ha_pt_t **ha_idx, All_reads *r, hifiasm_opt_ sprintf(gfa_name, "%s.ad", file_name); if(is_w) { - write_pt_index(*flt_tab, *ha_idx, NULL, opt, gfa_name); + write_pt_index(*flt_tab, *ha_idx, NULL, opt, gfa_name, 0); } else { - load_pt_index(flt_tab, ha_idx, NULL, opt, gfa_name); + load_pt_index(flt_tab, ha_idx, NULL, opt, gfa_name, 0); } free(gfa_name); } -uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rcc, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load) +uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rcc, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load, uint8_t force_rpaf_load) { char* gfa_name = (char*)malloc(strlen(file_name)+64); FILE *fp = NULL; int f_flag = 0; uint64_t rr0 = (uint64_t)-1, tot_rr0 = (uint64_t)-1; @@ -1649,7 +1683,7 @@ uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rc } sprintf(gfa_name, "%s.r%lu", file_name, rr); - if(!load_pt_index(r_flt_tab, r_ha_idx, r, opt, gfa_name)) { + if(!load_pt_index(r_flt_tab, r_ha_idx, r, opt, gfa_name, force_rpaf_load)) { destory_All_reads(r); ha_pt_destroy(*r_ha_idx); (*r_ha_idx) = NULL; ha_ft_destroy(*r_flt_tab); (*r_flt_tab) = NULL; @@ -1673,7 +1707,7 @@ uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rc write_cc_v(rcc, gfa_name); sprintf(gfa_name, "%s.r%lu", file_name, rr); - write_pt_index(*r_flt_tab, *r_ha_idx, r, opt, gfa_name); + write_pt_index(*r_flt_tab, *r_ha_idx, r, opt, gfa_name, force_rpaf_load); sprintf(gfa_name, "%s.r%lu.ht.bin", file_name, rr); fp = fopen(gfa_name, "w"); diff --git a/htab.h b/htab.h index 827df03..44d7983 100644 --- a/htab.h +++ b/htab.h @@ -86,15 +86,15 @@ const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n); const ha_idxposl_t *ha_ptl_get(const ha_pt_t *h, uint64_t hash, int *n); const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash); -int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name); -int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name); +int write_pt_index(void *flt_tab, ha_pt_t *ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name, uint8_t force_rpaf_load); +int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_opt_t* opt, char* file_name, uint8_t force_rpaf_load); void refresh_pt_idx(void **flt_tab, ha_pt_t **ha_idx, All_reads *r, hifiasm_opt_t *opt, char *file_name, uint8_t is_w); int uidx_write(void *flt_tab, ha_pt_t *ha_idx, char* file_name, ma_ug_t *ug); int uidx_load(void **r_flt_tab, ha_pt_t **r_ha_idx, char* file_name, ma_ug_t *ug); int write_ct_index(void *ct_idx, char* file_name); int load_ct_index(void **ct_idx, char* file_name); int query_ct_index(void* ct_idx, uint64_t hash); -uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rcc, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load); +uint64_t tmp_pt_pro(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads *r, cc_v* rcc, hifiasm_opt_t *opt, char *file_name, uint64_t rr, uint64_t tot_rr, uint64_t is_load, uint8_t force_rpaf_load); ha_abuf_t *ha_abuf_init_buf(void *km); ha_abufl_t *ha_abufl_init_buf(void *km); diff --git a/main.cpp b/main.cpp index 3ce4a0c..e521a3e 100644 --- a/main.cpp +++ b/main.cpp @@ -87,6 +87,14 @@ int main(int argc, char *argv[]) fprintf(stderr, "[M::%s::final] using %s mode\n", __func__, (asm_opt.simd_mm == 2) ? "AVX-512" : ((asm_opt.simd_mm == 1) ? "AVX2" : "non-SIMD")); + fprintf(stderr, "[M::%s::]\traw_aln::%d\tpost_syn::%d\n", __func__, asm_opt.realn_raw, asm_opt.post_syn); + + if(asm_opt.recurrent_err_test) { + fprintf(stderr, "[M::%s::] Enable recurrent sequencing-error filtering\n", __func__); + } else { + fprintf(stderr, "[M::%s::] Disable recurrent sequencing-error filtering\n", __func__); + } + if(asm_opt.sec_in) ret = ha_assemble_pair(); else if(asm_opt.dbg_ovec_cal) ret = ha_ec_dbg();