From 80fa5ed4360bed6aa0fbd8ac1418c798b6069ced Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Tue, 12 Nov 2024 23:53:29 -0500 Subject: [PATCH] ONT EC --- Assembly.cpp | 38 +- CommandLines.cpp | 8 +- CommandLines.h | 3 +- Correct.cpp | 987 +++++++++++++++++++++++++++++++++++++++++------ Correct.h | 8 +- Hash_Table.h | 1 + Process_Read.cpp | 169 ++++++++ Process_Read.h | 18 + ecovlp.cpp | 102 ++++- ecovlp.h | 1 + htab.cpp | 11 + 11 files changed, 1211 insertions(+), 135 deletions(-) diff --git a/Assembly.cpp b/Assembly.cpp index aa500f7..7b4bfdd 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -904,6 +904,29 @@ void Output_corrected_reads() fclose(output_file); } +void Output_corrected_fastq() +{ + long long i; uint64_t k; + UC_Read g_read; asg8_v dv; + init_UC_Read(&g_read); kv_init(dv); + char* gfa_name = (char*)malloc(strlen(asm_opt.output_file_name)+35); + sprintf(gfa_name, "%s.ec.fq", asm_opt.output_file_name); + FILE* fp = fopen(gfa_name, "w"); + free(gfa_name); + + for (i = 0; i < (long long)R_INF.total_reads; i++) { + recover_UC_Read(&g_read, &R_INF, i); + fprintf(fp, "@%.*s\n", (int32_t)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + fprintf(fp, "%.*s\n", (int32_t)g_read.length, g_read.seq); + fprintf(fp, "+\n"); + retrive_bqual(&dv, NULL, i, -1, -1, 0, sc_bn); + for (k = 0; k < dv.n; k++) fprintf(fp, "%c", (char)(sc_tb[dv.a[k]] + 33 - 1)); + fprintf(fp, "\n"); + } + destory_UC_Read(&g_read); kv_destroy(dv); + fclose(fp); +} + void debug_print_pob_regions() { uint64_t i, total = 0; @@ -992,9 +1015,12 @@ void ha_ec(int64_t round, int num_pround, int des_idx, uint64_t *tot_b, uint64_t if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name); + // Output_corrected_fastq(); + + cal_ec_r(asm_opt.thread_num, round, num_pround, R_INF.total_reads, (round == (asm_opt.number_of_round-1))?1:0, tot_b, tot_e); - // exit(1); + // exit(1); // if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name); if(des_idx) { @@ -2014,7 +2040,10 @@ int ha_assemble(void) ha_extract_print_list(&R_INF, asm_opt.extract_iter, asm_opt.extract_list); exit(0); } - if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads(); + if (asm_opt.flag & HA_F_WRITE_EC) { + if(asm_opt.is_sc) Output_corrected_fastq(); + else Output_corrected_reads(); + } if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF(); if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2; } @@ -2042,7 +2071,10 @@ int ha_assemble(void) // asm_opt.num_bases, asm_opt.num_corrected_bases, asm_opt.num_recorrected_bases); // fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0); } - if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads(); + if (asm_opt.flag & HA_F_WRITE_EC) { + if(asm_opt.is_sc) Output_corrected_fastq(); + else Output_corrected_reads(); + } // overlap between corrected reads ha_opt_reset_to_round(&asm_opt, asm_opt.number_of_round); // ha_overlap_final(); diff --git a/CommandLines.cpp b/CommandLines.cpp index bb1dcb0..a53f865 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -72,6 +72,7 @@ static ko_longopt_t long_options[] = { { "telo-s", ko_required_argument, 357}, { "ctg-n", ko_required_argument, 358}, { "ont", ko_no_argument, 359}, + { "sc-n", ko_no_argument, 360}, // { "path-round", ko_required_argument, 348}, { 0, 0, 0 } }; @@ -340,7 +341,8 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->telo_mic_sc = 500; asm_opt->is_ont = 0; -} + asm_opt->is_sc = 0; +} void destory_enzyme(enzyme* f) { @@ -912,7 +914,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 357) asm_opt->telo_mic_sc = atol(opt.arg); else if (c == 358) asm_opt->max_contig_tip = atol(opt.arg); else if (c == 359) { - asm_opt->is_ont = 1; asm_opt->max_ov_diff_ec = 0.07; + asm_opt->is_ont = 1; asm_opt->max_ov_diff_ec = 0.07; ///asm_opt->mz_win = 37; asm_opt->k_mer_length = 37; + } else if (c == 360) { + asm_opt->is_sc = 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); } diff --git a/CommandLines.h b/CommandLines.h index a3bbbcf..bd483d6 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.20.0-r656" +#define HA_VERSION "0.21.0-r666" #define VERBOSE 0 @@ -162,6 +162,7 @@ typedef struct { int64_t telo_mic_sc; uint64_t is_ont; + uint64_t is_sc; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Correct.cpp b/Correct.cpp index e7b445d..fa0c3e0 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include "Correct.h" #include "Levenshtein_distance.h" #include "Assembly.h" @@ -15,6 +16,7 @@ #include "htab.h" #include "Overlaps.h" #include "inter.h" +#include "Process_Read.h" #define A_L 16 #define ext_w 6 #define r_simi_w 0.05 @@ -30,6 +32,9 @@ KRADIX_SORT_INIT(haplotype_evdience_srt, haplotype_evdience, haplotype_evdience_ #define haplotype_evdience_id_key(x) ((x).overlapID) KRADIX_SORT_INIT(haplotype_evdience_id_srt, haplotype_evdience, haplotype_evdience_id_key, member_size(haplotype_evdience, overlapID)) +#define haplotype_evdience_os_key(x) ((x).overlapSite) +KRADIX_SORT_INIT(haplotype_evdience_os_srt, haplotype_evdience, haplotype_evdience_os_key, member_size(haplotype_evdience, overlapSite)) + #define overlap_region_dp_key(x) ((x).x_pos_e) KRADIX_SORT_INIT(overlap_region_dp_srt, overlap_region, overlap_region_dp_key, member_size(overlap_region, x_pos_e)) @@ -8883,7 +8888,7 @@ 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) { for (i = l, o = 0; i < k; i++) { - if(hap->list[i].type!=1) continue;///mismatch + if(hh_tp(hap->list[i])!=1) continue;///mismatch s = &(hap->snp_stat.a[hap->list[i].overlapSite]); assert(s->site == hap->list[i].site); if(s->occ_0 < 2 || s->occ_1 < 2) continue; @@ -8894,7 +8899,7 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio // if(overlap_list->list[hap->list[l].overlapID].y_id == 4290) { // fprintf(stderr, "[M::%s-id::%u] o->%lu(%c)\n", __func__, overlap_list->list[hap->list[l].overlapID].y_id, o, "+-"[overlap_list->list[hap->list[l].overlapID].y_pos_strand]); // for (i = l, o = 0; i < k; i++) { - // if(hap->list[i].type!=1) continue;///mismatch + // if(hh_tp(hap->list[i])!=1) continue;///mismatch // s = &(hap->snp_stat.a[hap->list[i].overlapSite]); // assert(s->site == hap->list[i].site); // if(s->occ_0 < 2 || s->occ_1 < 2) continue; @@ -8907,7 +8912,7 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio // } // else { // for (i = l, o = 0; i < k; i++) { - // if(hap->list[i].type!=0) continue;///mismatch + // if(hh_tp(hap->list[i])!=0) continue;///mismatch // s = &(hap->snp_stat.a[hap->list[i].overlapSite]); // assert(s->site == hap->list[i].site); // if(s->occ_0 < 2 || s->occ_1 < 2) continue; @@ -8925,7 +8930,7 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio // } // } - // if(overlap_list->list[hap->list[l].overlapID].y_id == 317 || overlap_list->list[hap->list[l].overlapID].y_id == 287) { + // if(overlap_list->list[hap->list[l].overlapID].y_id == 20835) { // fprintf(stderr, "***0***[M::%s-id::%u] o->%lu\n", __func__, overlap_list->list[hap->list[l].overlapID].y_id, o); // } @@ -8943,13 +8948,13 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio for (k = 0; k < hap->snp_srt.n; k++) { o = 0; l = (uint32_t)hap->snp_srt.a[k]; for (i = l; i < hap->length && hap->list[i].overlapID == hap->list[l].overlapID; i++) { - if(hap->list[i].type!=1) continue; + if(hh_tp(hap->list[i])!=1) continue; s = &(hap->snp_stat.a[hap->list[i].overlapSite]); if(s->occ_0 < 2 || s->occ_1 < 2) continue; if(is_st_bs((*s), st_rate, st_max)) continue; if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) { o++; - // if(overlap_list->list[hap->list[l].overlapID].y_id == 3276) { + // if(overlap_list->list[hap->list[l].overlapID].y_id == 20835) { // fprintf(stderr, "[M::%s-id::%u] occ_0->%u, occ_1->%u, occ_2->%u, site->%u\n", __func__, overlap_list->list[hap->list[l].overlapID].y_id, s->occ_0, s->occ_1, s->occ_2, s->site); // } } @@ -8962,16 +8967,19 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio ii = hap->list[l].overlapID; if(overlap_list->list[ii].is_match == 1) overlap_list->list[ii].is_match = 2; for (i = l; i < hap->length && hap->list[i].overlapID == hap->list[l].overlapID; i++) { - if(hap->list[i].type==1){ + if(hh_tp(hap->list[i])==1){ s = &(hap->snp_stat.a[hap->list[i].overlapSite]); s->score = 1; - } ///else if((hap->list[i].type==0) && (o>=(overlap_list->list[ii].align_length*up))) { - else if(hap->list[i].type==0) { + } ///else if((hh_tp(hap->list[i])==0) && (o>=(overlap_list->list[ii].align_length*up))) { + else if(hh_tp(hap->list[i])==0) { ///not real allels z = hap->list[i].overlapSite; s = &(hap->snp_stat.a[z]); for (z = hap->list[i].overlapSite; z >= 0; z--) { t = &(hap->snp_stat.a[z]); if(s->site!=t->site) break; + // if(t->site == 14217) { + // fprintf(stderr, "[M::%s]\tsite::%u\tn0::%u\tn1::%u\to::%lu\t%.*s\n", __func__, t->site, t->occ_0, t->occ_1, o, (int)Get_NAME_LENGTH(R_INF, overlap_list->list[ii].y_id), Get_NAME(R_INF, overlap_list->list[ii].y_id)); + // } t->occ_0 -= hap->list[i].cov; assert(t->occ_0 >= 1); if((st_max != ((uint64_t)-1)) && (overlap_list->list[ii].y_pos_strand == 0)) { @@ -8986,7 +8994,7 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio for (k = 0; k < hap->snp_srt.n; k++) {///sorted by how many allels in each overlap; more -> less o = 0; l = (uint32_t)hap->snp_srt.a[k]; for (i = l; i < hap->length && hap->list[i].overlapID == hap->list[l].overlapID; i++) { - if(hap->list[i].type!=1) continue; + if(hh_tp(hap->list[i])!=1) continue; s = &(hap->snp_stat.a[hap->list[i].overlapSite]); if(s->occ_0 < 2 || s->occ_1 < 2) continue; if(is_st_bs((*s), st_rate, st_max)) continue; @@ -9008,7 +9016,7 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio ii = hap->list[l].overlapID; if(overlap_list->list[ii].is_match==1) { for (i = l; i < k; i++) { - if(hap->list[i].type==1) { + if(hh_tp(hap->list[i])==1) { hap->snp_stat.a[hap->list[i].overlapSite].score = -1; } } @@ -9028,7 +9036,7 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio continue; } for (i = l, o = 0; i < k; i++) { - if(hap->list[i].type!=1) continue; + if(hh_tp(hap->list[i])!=1) continue; s = &(hap->snp_stat.a[hap->list[i].overlapSite]); if(s->occ_0 < 2 || s->occ_1 < 2) continue; if(is_st_bs((*s), st_rate, st_max)) continue; @@ -9078,11 +9086,11 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio 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(hap->list[i].type==1 || hap->list[i].type==0) { + 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(hap->list[i].type==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; @@ -9098,6 +9106,673 @@ void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_regio } } +inline int64_t comput_sc_rphase(SnpStats *ai, uint64_t id, SnpStats *aj, uint64_t jd, haplotype_evdience *za, uint64_t occ0_cut) +{ + if(ai->site == aj->site) return INT64_MIN; + // if(ai->occ_0 < occ0_cut || aj->occ_0 < occ0_cut) return INT64_MIN; + haplotype_evdience *iz = NULL, *jz = NULL; int64_t in, jn, ik, jk, nn[2]; uint8_t fi, fj; + iz = za + ai->non_homopolymer_num; in = ai->homopolymer_num - ai->non_homopolymer_num; + jz = za + aj->non_homopolymer_num; jn = aj->homopolymer_num - aj->non_homopolymer_num; + + for (ik = jk = nn[0] = nn[1] = 0; (ik < in) && (jk < jn); ik++) { + for (; (jk < jn) && (jz[jk].overlapID < iz[ik].overlapID); jk++); + if((jk < jn) && (jz[jk].overlapID == iz[ik].overlapID)) { + + fi = 2; + if(hh_tp(iz[ik]) == 0) { + fi = 0; + } else if(iz[ik].overlapSite == id){ + fi = 1; + } + + fj = 2; + if(hh_tp(jz[jk]) == 0) { + fj = 0; + } else if(jz[jk].overlapSite == jd){ + fj = 1; + } + + if((fi == 2) && (fj == 2) && (iz[ik].overlapSite != ((uint32_t)-1)) && (jz[jk].overlapSite != ((uint32_t)-1))) {///for rare cases + fi = fj = 0; + } + + if(fi == 2 || fj == 2) return INT64_MIN; + if(fi != fj) return INT64_MIN; + nn[fi]++; + } + } + + if(nn[0] > 0 && nn[1] > 0) return 1; + return INT64_MIN; +} + +///idx->a:: [0, ch_n) -> tree; +void gen_rphase_path(asg64_v *idx, int64_t *cn, int32_t *f, int64_t *p, uint64_t v0, asg64_v *res) +{ + uint64_t n0 = idx->n, pn = 0, v, w, ct, *ca, k, gi; + kv_push(uint64_t, *idx, v0); + while (idx->n > n0) { + v = idx->a[--idx->n]; + pn = v>>32; v = (uint32_t)v; + f[v] = 1; p[pn++] = v; + ct = ((v>0)?(cn[v]-cn[v-1]):(cn[v])); + ca = ((v>0)?(idx->a+cn[v-1]):(idx->a)); + if(ct) { + for (k = 0; k < ct; k++) { + w = ((pn<<32)|ca[k]); + kv_push(uint64_t, *idx, w); + } + } else { + gi = ((res->n)?((res->a[res->n-1]>>32)+1):(0)); + for (k = 0; k < pn; k++) { + w = (gi<<32)|((uint32_t)p[pn-k-1]); + kv_push(uint64_t, *res, w); + } + } + } +} + +void dbg_prt_rphase_vec(SnpStats *ai, uint64_t id, haplotype_evdience *za) +{ + haplotype_evdience *iz = NULL; int64_t in, ik; uint8_t fi; + iz = za + ai->non_homopolymer_num; in = ai->homopolymer_num - ai->non_homopolymer_num; + + for (ik = 0; ik < in; ik++) { + fi = 2; + if(hh_tp(iz[ik]) == 0) { + fi = 0; + } else if(iz[ik].overlapSite == id){ + fi = 1; + } + fprintf(stderr, "oid::%u(f::%u)\n", iz[ik].overlapID, fi); + } +} + +void gen_rphase_dp0_multiple_path(SnpStats *a, int64_t an, haplotype_evdience *za, Chain_Data *dp, asg64_v *idx, asg64_v *res) +{ + if(an <= 0) return; + int64_t *p, *t, i, k, j, max_f, sc, ch_n, plus = 0, rn; int32_t *f, *ii; uint64_t m; + resize_Chain_Data(dp, an, NULL); idx->n = res->n = 0; + t = dp->tmp; f = dp->score; p = dp->pre; ii = dp->occ; + + for (i = 0; i < an; ++i) { + max_f = 1; ch_n = idx->n; + for (j = i - 1; j >= 0; --j) { + sc = comput_sc_rphase(&a[i], i, &a[j], j, za, 0); + if (sc == INT64_MIN) continue; + sc += f[j]; + if (sc > max_f) { + max_f = sc; + idx->n = ch_n; kv_push(uint64_t, *idx, j); + } else if (sc == max_f) { + kv_push(uint64_t, *idx, j); + } + } + + f[i] = max_f; p[i] = idx->n; ii[i] = 0; + if(f[i] < plus) plus = f[i]; + } + ch_n = idx->n; + + for (i = 0; i < an; i++) { + f[i] -= plus; + m = f[i]; m = ((uint32_t)-1) - m; m <<= 32; m |= i; + kv_push(uint64_t, *idx, m); + } + radix_sort_bc64(idx->a + ch_n, idx->a + idx->n); + + for (i = 0; i < an; i++) {///idx->a:: [0, ch_n) -> tree; [ch_n, ch_n + an) -> sort; [ch_n + an, -1) -> path + k = (uint32_t)idx->a[ch_n + i]; + if(ii[k]) continue; + gen_rphase_path(idx, p, ii, t, k, res); + } + + rn = res->n; + for (k = 1, i = 0; k <= rn; k++) { + if(k == rn || (res->a[i]>>32) != (res->a[k]>>32)) { + if(k - i > 0) { + fprintf(stderr, "\n[M::%s]\tcnt::%lu\n", __func__, k - i); + for (j = i; j < k; j++) { + fprintf(stderr, "pos::%u\n", a[(uint32_t)res->a[j]].site); + // if(k - i == 2) { + // dbg_prt_rphase_vec(&(a[(uint32_t)res->a[j]]), (uint32_t)res->a[j], za); + // } + } + } + i = k; + } + } + +} + +int64_t is_hpc_vec(SnpStats *ai, uint64_t id, haplotype_evdience *za) +{ + haplotype_evdience *iz = NULL; int64_t in, ik, n0 = ai->occ_0, n1 = ai->occ_1, f = 0; + iz = za + ai->non_homopolymer_num; in = ai->homopolymer_num - ai->non_homopolymer_num; + + for (ik = 0; ik < in; ik++) { + if(!hh_hp(iz[ik])) continue; + if(hh_tp(iz[ik]) == 0) { + ai->occ_0 -= iz[ik].cov; + } else if(iz[ik].overlapSite == id){ + ai->occ_1 -= iz[ik].cov; + } + } + if((ai->occ_0 < 2 || ai->occ_1 < 2) || (!(ai->occ_0 >= asm_opt.s_hap_cov && ai->occ_1 >= asm_opt.infor_cov))) f = 1; + ai->occ_0 = n0; ai->occ_1 = n1; + // if((sec_check) && (((n0)<=((n0+n1)*0.333333)) || ((n1)<=((n0+n1)*0.333333)))) f = 1; + return f; +} + + +void prt_dbg_vec_comp() +{ + ; +} + +void get_hq_value(SnpStats *ai, uint64_t id, haplotype_evdience *za, int64_t *lq0, int64_t *hq0, int64_t *lq1, int64_t *hq1, uint8_t *ref_a) +{ + haplotype_evdience *z; int64_t n, k, occ0[2], occ1[2]; uint8_t fi; + z = za + ai->non_homopolymer_num; n = ai->homopolymer_num - ai->non_homopolymer_num; + occ0[0] = occ0[1] = occ1[0] = occ1[1] = 0; + for (k = 0; k < n; k++) { + if((hh_tp(z[k]) == 1) && (z[k].overlapSite != id)) continue; + fi = 0; + if((hh_bq(z[k]) > 1)) fi = 1; + + if(hh_tp(z[k]) == 0) occ0[fi]++; + else occ1[fi]++; + } + if(ref_a) { + if(ref_a[ai->site] > 1) occ0[1]++; + else occ0[0]++; + } + (*lq0) = occ0[0]; (*hq0) = occ0[1]; + (*lq1) = occ1[0]; (*hq1) = occ1[1]; +} + +void gen_rphase_dp0_single_path(SnpStats *a, int64_t an, haplotype_evdience *za, Chain_Data *dp, asg64_v *idx, int64_t het_cov, int64_t hom_cov, int64_t n_hap, double cut_rate, uint64_t cut_bd, asg64_v *res, uint8_t *qual_a) +{ + if(an <= 0) return; + int64_t *p, i, k, j, st, max_f, max_j, sc, ch_n, plus = 0, rn, rn0; int32_t *f, *ii; uint64_t m, cc = 0; + resize_Chain_Data(dp, an, NULL); idx->n = res->n = 0; + f = dp->score; p = dp->pre; ii = dp->occ; + cc = ((het_cov > 0)?(het_cov):(hom_cov/n_hap)); cc *= cut_rate; if(cc < cut_bd) cc = cut_bd; + + // fprintf(stderr, "\n\n\n[M::%s]\tcc::%lu\n\n\n", __func__, cc); + + for (i = 0; i < an; ++i) { + max_f = 1; max_j = -1; + st = 0; ///if(a[i].occ_0 < cc) st = i; + for (j = i - 1; j >= st; --j) { + sc = comput_sc_rphase(&a[i], i, &a[j], j, za, 0/**cc**/); + if (sc == INT64_MIN) continue; + sc += f[j]; + if (sc > max_f) { + max_f = sc; max_j = j; + } + } + + f[i] = max_f; p[i] = max_j; ii[i] = 0; + if(f[i] < plus) plus = f[i]; + } + + for (i = 0; i < an; i++) { + f[i] -= plus; + m = f[i]; m = ((uint32_t)-1) - m; m <<= 32; m |= i; + kv_push(uint64_t, *idx, m); + } + radix_sort_bc64(idx->a, idx->a + idx->n); + + kv_resize(uint64_t, *res, ((uint64_t)an)); + for (i = rn = 0; i < an; i++) {///idx->a:: [0, ch_n) -> tree; [ch_n, ch_n + an) -> sort; [ch_n + an, -1) -> path + rn0 = rn; + for(k = (uint32_t)idx->a[i]; (k >= 0) && (!ii[k]);) { + res->a[rn++] = k; ii[k] = 1; k = p[k]; + } + if(rn0 == rn) continue; + m = rn - rn0; m = ((uint32_t)-1) - m; m <<= 32; m |= rn0; + kv_push(uint64_t, *idx, m); + } + + radix_sort_bc64(idx->a + an, idx->a + idx->n); + ch_n = idx->n; + if(!qual_a) { + for (k = an; k < ch_n; k++) { + rn0 = ((uint32_t)idx->a[k]); + rn = ((uint32_t)-1) - (idx->a[k]>>32); + // break; + plus = -1; + if(rn > 1) { + // for (i = 0; (i < rn) && (is_hpc_vec(&(a[res->a[rn0+rn-i-1]]), res->a[rn0+rn-i-1], za)); i++); + // plus = ((i < rn)?(1):(-1)); + plus = 1; + + // for (i = 0; i < rn; i++) { + // a[res->a[rn0+rn-i-1]].score = plus; + // // fprintf(stderr, "pos::%u\tsc::%ld\tn0::%u\tn1::%u\tn2::%u\tk::%lu\n", a[res->a[rn0+rn-i-1]].site, plus, a[res->a[rn0+rn-i-1]].occ_0, a[res->a[rn0+rn-i-1]].occ_1, + // // 1 + a[res->a[rn0+rn-i-1]].homopolymer_num - a[res->a[rn0+rn-i-1]].non_homopolymer_num - a[res->a[rn0+rn-i-1]].occ_0 - a[res->a[rn0+rn-i-1]].occ_1, res->a[rn0+rn-i-1]); + // } + + } else { + if((!is_hpc_vec(&(a[res->a[rn0]]), res->a[rn0], za)) && (a[res->a[rn0]].occ_0 >= cc)) plus = 1; + ///a[res->a[rn0]].score = plus; + // fprintf(stderr, "pos::%u\tsc::%ld\tn0::%u\tn1::%u\tn2::%u\tk::%lu\n", a[res->a[rn0]].site, plus, a[res->a[rn0]].occ_0, a[res->a[rn0]].occ_1, + // 1 + a[res->a[rn0]].homopolymer_num - a[res->a[rn0]].non_homopolymer_num - a[res->a[rn0]].occ_0 - a[res->a[rn0]].occ_1, res->a[rn0]); + } + + // fprintf(stderr, "\n[M::%s]\tcnt::%lu\tplus::%ld\n", __func__, rn, plus); + // if(plus > 0) { + // for (i = 0; (i < rn) && (a[res->a[rn0 + i]].occ_0 < cc); i++); + // if(i < rn) { + // for (i = 0; i < rn; i++) a[res->a[rn0 + i]].score = plus; + // } + // } + for (i = 0; i < rn; i++) { + if(a[res->a[rn0 + i]].occ_0 >= cc) { + a[res->a[rn0 + i]].score = plus; + } else { + a[res->a[rn0 + i]].score = -1; + } + // fprintf(stderr, "pos::%u\tsc::%d\tn0::%u\tn1::%u\tn2::%u\tk::%lu\n", a[res->a[rn0 + i]].site, a[res->a[rn0 + i]].score, a[res->a[rn0 + i]].occ_0, a[res->a[rn0 + i]].occ_1, + // 1 + a[res->a[rn0 + i]].homopolymer_num - a[res->a[rn0 + i]].non_homopolymer_num - a[res->a[rn0 + i]].occ_0 - a[res->a[rn0 + i]].occ_1, res->a[rn0 + i]); + } + } + } else { + int64_t b0l, b0h, b1l, b1h, krn; + /** + for (k = an; k < ch_n; k++) { + rn0 = ((uint32_t)idx->a[k]); + rn = ((uint32_t)-1) - (idx->a[k]>>32); + // break; + fprintf(stderr, "\n[M::%s]\tcnt::%lu\n", __func__, rn); + if(rn > 1) { + for (i = 0; i < rn; i++) { + get_hq_value(&(a[res->a[rn0 + i]]), res->a[rn0 + i], za, &b0l, &b0h, &b1l, &b1h, qual_a); + if(b0h > b0l && b0h >= asm_opt.s_hap_cov && b1h > b1l && b1h >= asm_opt.infor_cov) { + a[res->a[rn0 + i]].score = 1; + } else { + a[res->a[rn0 + i]].score = -1; + } + fprintf(stderr, "pos::%u\tsc::%d\tn0::%u\tn1::%u\tn2::%u\tk::%lu\tb0l::%ld\tb0h::%ld\tb1l::%ld\tb1h::%ld\n", a[res->a[rn0 + i]].site, a[res->a[rn0 + i]].score, a[res->a[rn0 + i]].occ_0, a[res->a[rn0 + i]].occ_1, + 1 + a[res->a[rn0 + i]].homopolymer_num - a[res->a[rn0 + i]].non_homopolymer_num - a[res->a[rn0 + i]].occ_0 - a[res->a[rn0 + i]].occ_1, res->a[rn0 + i], b0l, b0h, b1l, b1h); + } + } else { + if((!is_hpc_vec(&(a[res->a[rn0]]), res->a[rn0], za))) { + get_hq_value(&(a[res->a[rn0]]), res->a[rn0], za, &b0l, &b0h, &b1l, &b1h, qual_a); + if(((b0h > b0l) && (b0h > ((b0h + b0l)*0.7)) && (b0h >= asm_opt.s_hap_cov) && (b0h >= (int64_t)cc)) && + ((b1h > b1l) && (b1h > ((b1h + b1l)*0.7)) && (b1h >= asm_opt.infor_cov) && (b1h >= (int64_t)cc))) { + a[res->a[rn0]].score = 1; + } else { + a[res->a[rn0]].score = -1; + } + } + } + } + **/ + for (k = an; k < ch_n; k++) { + rn0 = ((uint32_t)idx->a[k]); + rn = ((uint32_t)-1) - (idx->a[k]>>32); + krn = rn; + // break; + + for (i = 1; (i < rn) && ((a[res->a[rn0 + i]].site + 8) >= a[res->a[rn0 + i - 1]].site); i++); + if(i >= rn) krn = 1; + // fprintf(stderr, "\n[M::%s]\tcnt::%lu\tkrn::%ld\n", __func__, rn, krn); + + for (i = 0; i < rn; i++) { + get_hq_value(&(a[res->a[rn0 + i]]), res->a[rn0 + i], za, &b0l, &b0h, &b1l, &b1h, qual_a); + if(krn > 1) { + if(b0h > b0l && b0h >= asm_opt.s_hap_cov && b1h > b1l && b1h >= asm_opt.infor_cov) { + a[res->a[rn0 + i]].score = 1; + } else { + a[res->a[rn0 + i]].score = -1; + } + } else { + if(((b0h > b0l) && (b0h > ((b0h + b0l)*0.7)) && (b0h >= asm_opt.s_hap_cov) && (b0h >= (int64_t)cc)) && + ((b1h > b1l) && (b1h > ((b1h + b1l)*0.7)) && (b1h >= asm_opt.infor_cov) && (b1h >= (int64_t)cc))) { + a[res->a[rn0 + i]].score = 1; + } else { + a[res->a[rn0 + i]].score = -1; + } + } + + // fprintf(stderr, "pos::%u\tsc::%d\tn0::%u\tn1::%u\tn2::%u\tk::%lu\tb0l::%ld\tb0h::%ld\tb1l::%ld\tb1h::%ld\n", a[res->a[rn0 + i]].site, a[res->a[rn0 + i]].score, a[res->a[rn0 + i]].occ_0, a[res->a[rn0 + i]].occ_1, + // 1 + a[res->a[rn0 + i]].homopolymer_num - a[res->a[rn0 + i]].non_homopolymer_num - a[res->a[rn0 + i]].occ_0 - a[res->a[rn0 + i]].occ_1, res->a[rn0 + i], b0l, b0h, b1l, b1h); + } + + for (i = j = 0; i < rn; i = j) { + plus = a[res->a[rn0 + i]].score; + for (j = i + 1; (j < rn) && ((a[res->a[rn0 + i]].site-a[res->a[rn0 + j]].site) == (j-i)); j++) { + if(a[res->a[rn0 + j]].score == -1) plus = -1; + } + // fprintf(stderr, "[M::%s]\ti::%ld\tj::%ld\tplus::%ld\n", __func__, i, j, plus); + if(plus == -1) { + for (; i < j; i++) { + a[res->a[rn0 + i]].score = plus; + } + } + } + + } + } + +} + +inline void fill_incom(asg64_v *om, uint64_t oid, uint64_t pe, uint64_t *idx_a, int64_t idx_n, uint64_t qid) +{ + if(om->a[oid] == ((uint64_t)-1)) { + om->a[oid] = pe; + return; + } + + uint64_t ps = om->a[oid]; int64_t k; + for (k = idx_n - 1; idx_a[k] != ps; k--); + // if(!(k >= 0)) { + // fprintf(stderr, "qid::%lu, oid::%lu, pe::%lu, ps::%lu, idx_n::%ld\n", qid, oid, pe, ps, idx_n); + // } + assert(k >= 0); + for (k++; k < idx_n; k++) { + kv_push(uint64_t, *om, ((oid<<32)|(idx_a[k]))); + } + om->a[oid] = pe; +} + +void get_wqual(uint64_t zid, uint64_t zpos, uint64_t zrev, asg8_v *v, uint8_t *va, uint64_t scw, uint64_t *tqual, uint64_t *wqual) +{ + (*tqual) = (*wqual) = 0; + uint64_t s, e, k, wk, wsc, msc, l = Get_READ_LENGTH((R_INF), (zid)); uint8_t *a; + s = ((zpos + 1 >= scw)?(zpos + 1 - scw):(0)); + e = zpos + scw; if(e > l) e = l; + + if(v) { + retrive_bqual(v, NULL, zid, s, e, zrev, sc_bn); a = v->a; + } else { + a = va + s; + } + + (*tqual) = a[zpos - s]; + for (k = s, wk = wsc = 0; k < e && wk < scw; k++, wk++) { + wsc += a[k - s]; + } + msc = wsc; + + if(k < e) { + assert(wk == scw); + for (; k < e; k++) { + wsc -= a[k - scw - s]; + wsc += a[k - s]; + if(wsc < msc) msc = wsc; + } + } + + (*wqual) = round(((double)msc)/((double)wk)); + // fprintf(stderr, "[M::%s] s::%lu, e::%lu, wk::%lu, scw::%lu, msc/wk::%lu, wqual::%lu\n", __func__, s, e, wk, scw, msc/wk, (*wqual)); +} + +void gen_rphase_dp(haplotype_evdience_alloc* hl, overlap_region_alloc* ol, UC_Read* g_read, double st_rate, uint64_t st_max, Chain_Data *dp, asg64_v *idx, asg64_v *res, uint64_t rid, asg8_v *qv) +{ + if(hl->length <= 0) return; + + uint64_t k, l, i, i0, o, ii, m_snp_stat, m_snp_stat0, m_list /**m_off**/, m1, c0, c1, rev_n, tqual, wqual, hq_cut = 2; uint32_t hq[5], hp[4], is_st; SnpStats *s; haplotype_evdience ev; char mc; + retrive_bqual(qv, NULL, rid, -1, -1, 0, sc_bn); + + idx->n = 0; + kv_resize(uint64_t, *res, ol->length); memset(res->a, -1, sizeof((*(res->a)))*ol->length); res->n = ol->length; + for (k = 1, l = 0, i = m_snp_stat = m_list = 0; k <= hl->snp_stat.n; ++k) {///filter snps + if(k == hl->snp_stat.n || hl->snp_stat.a[k].site != hl->snp_stat.a[l].site) { + for (; i < hl->length && hl->list[i].site != hl->snp_stat.a[l].site; i++); + assert(i < hl->length && hl->list[i].site == hl->snp_stat.a[l].site); + for (o = l, m1 = 0; o < k; o++) { + s = &(hl->snp_stat.a[o]); + // fprintf(stderr, "+[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, s->site, s->occ_0, s->occ_1); + if((s->occ_0 < 2 || s->occ_1 < 2) || (is_st_bs((*s), st_rate, st_max)) || (!(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov))) { + continue; + } + // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, s->site, s->occ_0, s->occ_1); + m1++; + } + + if(m1) { + hq[0] = hq[1] = hq[2] = hq[3] = hq[4] = 0; hp[0] = hp[1] = hp[2] = hp[3] = ((uint32_t)-1); + for (i0 = i, c0 = rev_n = 0; i0 < hl->length && hl->list[i0].site == hl->snp_stat.a[l].site; i0++) { + hl->list[i0].site = ((uint32_t)-1); + if(hh_tp(hl->list[i0]) == 0) { + c0 += 1/**hl->list[i0].cov**/; + if((ol->list[hl->list[i0].overlapID].y_pos_strand == 0)) rev_n += 1/**hl->list[i0].cov**/; + } + if(qv) { + if(hh_tp(hl->list[i0]) == 0) { + if((hh_bq(hl->list[i0]) > 1) || (hh_wq(hl->list[i0]) > 1)) hq[4]++; + } else { + hp[seq_nt6_table[(uint8_t)(hl->list[i0].misBase)]] = hl->list[i0].overlapSite; + if((hh_bq(hl->list[i0]) > 1) || (hh_wq(hl->list[i0]) > 1)) { + hq[seq_nt6_table[(uint8_t)(hl->list[i0].misBase)]]++; + } + } + } + } + assert(i0 > i); + if(qv) { + if(hp[0] != ((uint32_t)-1) && hq[0] < hq_cut) { + // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, hl->snp_stat.a[hp[0]].site, hl->snp_stat.a[hp[0]].occ_0, hl->snp_stat.a[hp[0]].occ_1); + hl->snp_stat.a[hp[0]].occ_1 = 0; + } + + if(hp[1] != ((uint32_t)-1) && hq[1] < hq_cut) { + // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, hl->snp_stat.a[hp[1]].site, hl->snp_stat.a[hp[1]].occ_0, hl->snp_stat.a[hp[1]].occ_1); + hl->snp_stat.a[hp[1]].occ_1 = 0; + } + + if(hp[2] != ((uint32_t)-1) && hq[2] < hq_cut) { + // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, hl->snp_stat.a[hp[2]].site, hl->snp_stat.a[hp[2]].occ_0, hl->snp_stat.a[hp[2]].occ_1); + hl->snp_stat.a[hp[2]].occ_1 = 0; + } + + if(hp[3] != ((uint32_t)-1) && hq[3] < hq_cut) { + // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, hl->snp_stat.a[hp[3]].site, hl->snp_stat.a[hp[3]].occ_0, hl->snp_stat.a[hp[3]].occ_1); + hl->snp_stat.a[hp[3]].occ_1 = 0; + } + + get_wqual(rid, hl->snp_stat.a[l].site, 0, NULL, qv->a, sc_wn, &tqual, &wqual); + if(tqual > 1 || wqual > 1) hq[4]++; + // if(hq[4] == 0) { + // fprintf(stderr, "+[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, hl->snp_stat.a[l].site, hl->snp_stat.a[l].occ_0, hl->snp_stat.a[l].occ_1); + // } + } else { + hq[0] = hq[1] = hq[2] = hq[3] = hq[4] = ((uint32_t)-1); + } + + if((hq[4] >= hq_cut) && (hq[0] >= hq_cut || hq[1] >= hq_cut || hq[2] >= hq_cut || hq[3] >= hq_cut)) { + for (o = l, m_snp_stat0 = m_snp_stat; o < k; o++) { + s = &(hl->snp_stat.a[o]); + if((s->occ_0 < 2 || s->occ_1 < 2) || (is_st_bs((*s), st_rate, st_max)) || (!(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov))) { + continue; + } + + hl->snp_stat.a[m_snp_stat] = *s; + for (ii = i, c1 = 0; ii < i0; ii++) { + if(((hl->list[ii].overlapSite == o) && (hl->list[ii].site == ((uint32_t)-1))) || (hh_tp(hl->list[ii]) == 0)) { + hl->list[ii].overlapSite = m_snp_stat; + hl->list[ii].site = hl->snp_stat.a[o].site; + if(hh_tp(hl->list[ii]) == 1) c1 += 1/**hl->list[ii].cov**/;///not correct if we don't do read correction + } + } + // fprintf(stderr, "-0-[M::%s]\tsite::%u\tc0::%lu\tc1::%lu\n", __func__, s->site, c0, c1); + hl->snp_stat.a[m_snp_stat].occ_0 = c0 + 1; + hl->snp_stat.a[m_snp_stat].occ_1 = c1; + hl->snp_stat.a[m_snp_stat].overlap_num = rev_n + 1; + hl->snp_stat.a[m_snp_stat].non_homopolymer_num = m_list; + hl->snp_stat.a[m_snp_stat++].score = -1; + } + + for (ii = i, m1 = m_list; ii < i0; ii++) { + if((hl->list[ii].site == ((uint32_t)-1))) { + hl->list[ii].overlapSite = ((uint32_t)-1); + hl->list[ii].site = hl->snp_stat.a[l].site; + } + fill_incom(res, hl->list[ii].overlapID, hl->list[ii].site, idx->a, idx->n, ol->list[0].x_id); + hl->list[m_list++] = hl->list[ii]; + } + radix_sort_haplotype_evdience_id_srt(hl->list + m1, hl->list + m_list); + for (ii = m_snp_stat0; ii < m_snp_stat; ii++) { + hl->snp_stat.a[ii].homopolymer_num = m_list; + } + + i = i0; + + kv_push(uint64_t, *idx, hl->snp_stat.a[l].site); + } + } + + l = k; + } + } + hl->snp_stat.n = m_snp_stat; hl->length = m_list; + if(hl->snp_stat.n == 0 || hl->length == 0) return; + + + ///fill incomplete values + if(res->n > ol->length) { + // fprintf(stderr, "\n[M::%s]\tin_comp::%lu\n", __func__, ((uint64_t)res->n - ol->length)); + for (k = ol->length; k < res->n; k++) { + ev.misBase = 0; + ev.overlapID = res->a[k]>>32; + ev.site = ((uint32_t)res->a[k]); + ev.overlapSite = ((uint32_t)-1); + ev.type = 0; + ev.cov = 1; + addHaplotypeEvdience(hl, &ev, NULL); + // fprintf(stderr, "[M::%s]\toid::%u\tsite::%u\n", __func__, ev.overlapID, ev.site); + } + radix_sort_haplotype_evdience_srt(hl->list, hl->list + hl->length); + + for (k = 1, l = 0, i = 0; k <= hl->snp_stat.n; ++k) { + if(k == hl->snp_stat.n || hl->snp_stat.a[k].site != hl->snp_stat.a[l].site) { + for (; i < hl->length && hl->list[i].site != hl->snp_stat.a[l].site; i++); + assert(i < hl->length && hl->list[i].site == hl->snp_stat.a[l].site); + + m1 = ((uint32_t)-1); mc = 0; + for (i0 = i, c0 = rev_n = 0; i0 < hl->length && hl->list[i0].site == hl->snp_stat.a[l].site; i0++) { + if(hh_tp(hl->list[i0]) == 0) { + c0+=1/**hl->list[i0].cov**/; + if((ol->list[hl->list[i0].overlapID].y_pos_strand == 0)) rev_n += 1/**hl->list[i0].cov**/; + if(hl->list[i0].overlapSite != ((uint32_t)-1)) { + m1 = hl->list[i0].overlapSite; mc = hl->list[i0].misBase; + } + } + // if(hl->snp_stat.a[l].site == 6995) { + // fprintf(stderr, "[M::%s]\tsite::%u\tosite::%u\toid::%u\ttp::%u\n", __func__, hl->list[i0].site, hl->list[i0].overlapSite, hl->list[i0].overlapID, hh_tp(hl->list[i0])); + // } + } + assert(i0 > i); + assert(m1 != ((uint32_t)-1)); + assert(mc != 0); + + + for (ii = i, is_st = 1; ii < i0; ii++) { + if(hh_tp(hl->list[ii]) == 0) { + assert(hl->list[ii].overlapSite == ((uint32_t)-1) || hl->list[ii].overlapSite == m1); + hl->list[ii].overlapSite = m1; hl->list[ii].misBase = mc; + } + if((ii > i) && (hl->list[ii].overlapID <= hl->list[ii-1].overlapID)) { + is_st = 0; + // fprintf(stderr, "[M::%s]\tid::[%u,%u]\tsite::[%u,%u]\n", __func__, hl->list[ii].overlapID, hl->list[ii-1].overlapID, hl->list[ii].overlapSite, hl->list[ii-1].overlapSite); + } + } + if(is_st == 0) radix_sort_haplotype_evdience_id_srt(hl->list + i, hl->list + i0); + // assert(is_st == 1); + + + + for (o = l; o < k; o++) { + s = &(hl->snp_stat.a[o]); + for (ii = i, c1 = 0; ii < i0; ii++) { + if((hl->list[ii].overlapSite == o) && (hh_tp(hl->list[ii]) == 1)) c1 += 1/**hl->list[ii].cov**/;///not correct if we don't do read correction + } + // fprintf(stderr, "-1-[M::%s]\tsite::%u\tc0::%lu\tc1::%lu\n", __func__, s->site, c0, c1); + s->occ_0 = c0 + 1; + s->occ_1 = c1; + s->overlap_num = rev_n + 1; + s->non_homopolymer_num = i; + s->homopolymer_num = i0; + } + + i = i0; + l = k; + } + } + + } + + + + + + // gen_rphase_dp0_multiple_path(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, res); + gen_rphase_dp0_single_path(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, asm_opt.het_cov, asm_opt.hom_cov, asm_opt.polyploidy, 0.7, 6, res, qv->a); + + for (k = 1, l = 0, i = m_snp_stat = m_list = 0; k <= hl->snp_stat.n; ++k) {///filter snps + if(k == hl->snp_stat.n || hl->snp_stat.a[k].site != hl->snp_stat.a[l].site) { + for (; i < hl->length && hl->list[i].site != hl->snp_stat.a[l].site; i++); + assert(i < hl->length && hl->list[i].site == hl->snp_stat.a[l].site); + for (o = l, m1 = 0; o < k; o++) { + s = &(hl->snp_stat.a[o]); + if(s->score == -1) continue; + m1++; + } + + if(m1) { + for (i0 = i, c0 = rev_n = 0; i0 < hl->length && hl->list[i0].site == hl->snp_stat.a[l].site; i0++) { + hl->list[i0].site = ((uint32_t)-1); + if(hh_tp(hl->list[i0]) == 0) { + c0 += 1/**hl->list[i0].cov**/; + if((ol->list[hl->list[i0].overlapID].y_pos_strand == 0)) rev_n += 1/**hl->list[i0].cov**/; + } + } + assert(i0 > i); + + for (o = l, m_snp_stat0 = m_snp_stat; o < k; o++) { + s = &(hl->snp_stat.a[o]); + if(s->score == -1) continue; + + hl->snp_stat.a[m_snp_stat] = *s; + for (ii = i, c1 = 0; ii < i0; ii++) { + if(((hl->list[ii].overlapSite == o) && (hl->list[ii].site == ((uint32_t)-1))) || (hh_tp(hl->list[ii]) == 0)) { + hl->list[ii].overlapSite = m_snp_stat; + hl->list[ii].site = hl->snp_stat.a[o].site; + if(hh_tp(hl->list[ii]) == 1) c1 += 1/**hl->list[ii].cov**/;///not correct if we don't do read correction + } + } + + // if(hl->snp_stat.a[m_snp_stat].site == 4) { + // fprintf(stderr, "[M::%s] occ_0->%lu, occ_1->%lu, site->%u\n", __func__, c0 + 1, c1, hl->snp_stat.a[m_snp_stat].site); + // } + + hl->snp_stat.a[m_snp_stat].occ_0 = c0 + 1; + hl->snp_stat.a[m_snp_stat].occ_1 = c1; + hl->snp_stat.a[m_snp_stat].overlap_num = rev_n + 1; + hl->snp_stat.a[m_snp_stat].non_homopolymer_num = m_list; + hl->snp_stat.a[m_snp_stat++].score = -1; + } + + for (ii = i, m1 = m_list; ii < i0; ii++) { + if((hl->list[ii].site == ((uint32_t)-1))) continue; + hl->list[m_list++] = hl->list[ii]; + } + for (ii = m_snp_stat0; ii < m_snp_stat; ii++) { + hl->snp_stat.a[ii].homopolymer_num = m_list; + } + + i = i0; + } + + l = k; + } + } + hl->snp_stat.n = m_snp_stat; hl->length = m_list; + +} void generate_haplotypes_naive_UL(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, void *km) { @@ -9585,23 +10260,102 @@ int insert_snp_ee(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a return m; } -int push_info(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, haplotype_evdience* u_a, overlap_region *oa) + +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) { - uint64_t i, m, occ_0, occ_1[6], occ_2, diff, rev_n; + int64_t s = ((p>=hpc_flk)?(p-hpc_flk):0), e = (((p+hpc_flk)<=sn)?(p+hpc_flk):(sn)), k, r, rc, zs, ze; + + for (r = 1; r <= hpc_rr; r++) { + rc = r * HPC_CC; + + ///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) > r) && ((ze - zs) >= rc)) { + // fprintf(stderr, "-0-[M::%s] p::%ld, hh::[%ld,%ld), f::%u, %.*s\n", __func__, p, zs, ze, f?1:0, (int32_t)(ze - zs), sa + zs); + if(f) { + for (k = MAX(p, zs); (k < ze) && (k - fsift < fn); k++) f[k - fsift] = 0; f[p - fsift] = 0; + } + return 1; + } + + ///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) > r) && ((ze - zs) >= rc)) { + // fprintf(stderr, "-1-[M::%s] p::%ld, hh::[%ld,%ld), f::%u, %.*s\n", __func__, p, zs, ze, f?1:0, (int32_t)(ze - zs), sa + zs); + if(f) { + for (k = MAX(p, zs); (k < ze) && (k - fsift < fn); k++) f[k - fsift] = 0; f[p - fsift] = 0; + } + return 1; + } + + ///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) > r) && ((ze - zs) >= rc)) { + // fprintf(stderr, "-2-[M::%s] p::%ld, hh::[%ld,%ld), f::%u, %.*s\n", __func__, p, zs, ze, f?1:0, (int32_t)(ze - zs), sa + zs); + if(f) { + for (k = MAX(p, zs); (k < ze) && (k - fsift < fn); k++) f[k - fsift] = 0; f[p - fsift] = 0; + } + return 1; + } + + ///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) > r) && ((ze - zs) >= rc)) { + // fprintf(stderr, "-3-[M::%s] p::%ld, hh::[%ld,%ld), f::%u, %.*s\n", __func__, p, zs, ze, f?1:0, (int32_t)(ze - zs), sa + zs); + if(f) { + for (k = MAX(p, zs); (k < ze) && (k - fsift < fn); k++) f[k - fsift] = 0; f[p - fsift] = 0; + } + return 1; + } + } + // fprintf(stderr, "-6-[M::%s] p::%ld, hh::[,), %.*s\n", __func__, p, (int32_t)(e - s), sa + s); + return 0; +} + +int push_info(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, haplotype_evdience* u_a, overlap_region *oa, asg8_v *v8, uint64_t scw) +{ + uint64_t i, k, m, occ_0, occ_1[6], occ_2, diff, rev_n; uint8_t ihpc = 0; occ_0 = occ_2 = diff = rev_n = 0; memset(occ_1, 0, sizeof(uint64_t)*6); - for (i = 0; i < a_n; i++) { - if(a[i].type == 0){ - occ_0 += a[i].cov; - if((oa) && (oa[a[i].overlapID].y_pos_strand == 0)) { - rev_n += a[i].cov; + radix_sort_haplotype_evdience_id_srt(a, a + a_n); + for (k = 1, m = i = 0; k <= a_n; k++) { + if((k == a_n) || (a[k].overlapID != a[m].overlapID)) { + a[i] = a[m]; + if(hh_tp(a[i]) == 0){ + occ_0 += a[i].cov; + if((oa) && (oa[a[i].overlapID].y_pos_strand == 0)) { + rev_n += a[i].cov; + } + if((!ihpc) && (hh_hp(a[i]))) ihpc = 1; + }else if(hh_tp(a[i]) == 1){ + occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]] += a[i].cov; + diff += a[i].cov; } - }else if(a[i].type == 1){ - occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]] += a[i].cov; - diff += a[i].cov; + occ_2 += a[i].cov; + + i++; m = k; } - occ_2 += a[i].cov; } + a_n = i; + + // for (i = 0; i < a_n; i++) { + // if(hh_tp(a[i]) == 0){ + // occ_0 += a[i].cov; + // if((oa) && (oa[a[i].overlapID].y_pos_strand == 0)) { + // rev_n += a[i].cov; + // } + // if((!ihpc) && (hh_hp(a[i]))) ihpc = 1; + // }else if(hh_tp(a[i]) == 1){ + // occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]] += a[i].cov; + // diff += a[i].cov; + // } + // occ_2 += a[i].cov; + // } + /** 1. if occ_0 = 0, that means all overlaps are different with this read at this site @@ -9617,6 +10371,7 @@ int push_info(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, } else { rev_n++; } + for (i = m = 0; i < 4; i++) { if(occ_1[i] >= 2){ kv_pushp(SnpStats, h->snp_stat, &p); @@ -9627,9 +10382,10 @@ int push_info(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, p->site = a[0].site; p->score = -1; p->overlap_num = rev_n; - p->is_homopolymer = 0; + p->is_homopolymer = ihpc; occ_1[i] = p->id; m++; + // fprintf(stderr, "[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, p->site, p->occ_0, p->occ_1); } else { occ_1[i] = (uint64_t)-1; } @@ -9637,18 +10393,39 @@ int push_info(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, occ_1[4] = occ_1[5] = (uint64_t)-1; if(m == 0) return 0; - for (i = m = 0; i < a_n; i++) { - // fprintf(stderr, "[M::%s] a[%lu].misBase->%c\n", __func__, i, a[i].misBase); - if(a[i].type == 0) { - a[i].overlapSite = h->snp_stat.n-1; - } else if(occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]]!=(uint64_t)-1){ - a[i].cov = a[i].overlapSite;///note: only renew cov here!!! - a[i].overlapSite = occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]]; - } else { - continue; + if(!v8) { + for (i = m = 0; i < a_n; i++) { + // fprintf(stderr, "[M::%s] a[%lu].misBase->%c\n", __func__, i, a[i].misBase); + if(hh_tp(a[i]) == 0) { + a[i].overlapSite = h->snp_stat.n-1; + } else if(occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]]!=(uint64_t)-1){ + a[i].cov = a[i].overlapSite;///note: only renew cov here!!! + a[i].overlapSite = occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]]; + } else { + continue; + } + u_a[m++] = a[i]; + } + } else { + uint64_t tpos, tqual, wqual; + for (i = m = 0; i < a_n; i++) { + tpos = a[i].overlapSite; + if(hh_tp(a[i]) == 0) { + a[i].overlapSite = h->snp_stat.n-1; + } else if(occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]]!=(uint64_t)-1){ + a[i].cov = a[i].overlapSite;///note: only renew cov here!!! + a[i].overlapSite = occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]]; + } else { + continue; + } + get_wqual(oa[a[i].overlapID].y_id, tpos, oa[a[i].overlapID].y_pos_strand, v8, NULL, scw, &tqual, &wqual); + a[i].type |= (tqual<<2); a[i].type |= (wqual<<(sc_bn+2)); + // assert(hh_bq(a[i]) == tqual); assert(hh_wq(a[i]) == wqual); + // fprintf(stderr, "[M::%s] site::%u, type::%u, occ0::%lu, occ2::%lu, tqual::%lu, wqual::%lu\n", __func__, a[i].site, hh_tp(a[i]), occ_0, occ_2, tqual, wqual); + u_a[m++] = a[i]; } - u_a[m++] = a[i]; } + return m; } @@ -17526,61 +18303,6 @@ int64_t extract_sub_cigar_err_rr(overlap_region *z, int64_t s, int64_t e, ul_ov_ return err; } -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 s = ((p>=hpc_flk)?(p-hpc_flk):0), e = (((p+hpc_flk)<=sn)?(p+hpc_flk):(sn)), k, r, rc, zs, ze; - - for (r = 1; r <= hpc_rr; r++) { - rc = r * HPC_CC; - - ///inlcuding p - for (k = p + r; (k < e) && (sa[k] == sa[k-r]); k++); ze = k; if(ze > e) ze = e; - for (k = p - 1; (k >= s) && (sa[k] == sa[k+r]); k--); zs = k + 1; if(zs < s) zs = s; - if(((ze - zs) > r) && ((ze - zs) >= rc)) { - // fprintf(stderr, "-0-[M::%s] p::%ld, hh::[%ld,%ld), f::%u, %.*s\n", __func__, p, zs, ze, f?1:0, (int32_t)(ze - zs), sa + zs); - if(f) { - for (k = MAX(p, zs); (k < ze) && (k - fsift < fn); k++) f[k - fsift] = 0; f[p - fsift] = 0; - } - return 1; - } - - ///do not inlcude p - for (k = p + r + 1; (k < e) && (sa[k] == sa[k-r]); k++); - zs = p + 1; if(zs < s) zs = s; ze = k; if(ze > e) ze = e; - if(((ze - zs) > r) && ((ze - zs) >= rc)) { - // fprintf(stderr, "-1-[M::%s] p::%ld, hh::[%ld,%ld), f::%u, %.*s\n", __func__, p, zs, ze, f?1:0, (int32_t)(ze - zs), sa + zs); - if(f) { - for (k = MAX(p, zs); (k < ze) && (k - fsift < fn); k++) f[k - fsift] = 0; f[p - fsift] = 0; - } - return 1; - } - - ///inlcuding p - for (k = p - r; (k >= s) && (sa[k] == sa[k+r]); k--); zs = k + 1; if(zs < s) zs = s; - for (k = p + 1; (k < e) && (sa[k] == sa[k-r]); k++); ze = k; if(ze > e) ze = e; - if(((ze - zs) > r) && ((ze - zs) >= rc)) { - // fprintf(stderr, "-2-[M::%s] p::%ld, hh::[%ld,%ld), f::%u, %.*s\n", __func__, p, zs, ze, f?1:0, (int32_t)(ze - zs), sa + zs); - if(f) { - for (k = MAX(p, zs); (k < ze) && (k - fsift < fn); k++) f[k - fsift] = 0; f[p - fsift] = 0; - } - return 1; - } - - ///do not inlcude p - for (k = p - r - 1; (k >= s) && (sa[k] == sa[k+r]); k--); - zs = k + 1; if(zs < s) zs = s; ze = p; if(ze > e) ze = e; - if(((ze - zs) > r) && ((ze - zs) >= rc)) { - // fprintf(stderr, "-3-[M::%s] p::%ld, hh::[%ld,%ld), f::%u, %.*s\n", __func__, p, zs, ze, f?1:0, (int32_t)(ze - zs), sa + zs); - if(f) { - for (k = MAX(p, zs); (k < ze) && (k - fsift < fn); k++) f[k - fsift] = 0; f[p - fsift] = 0; - } - return 1; - } - } - // fprintf(stderr, "-6-[M::%s] p::%ld, hh::[,), %.*s\n", __func__, p, (int32_t)(e - s), sa + s); - return 0; -} - ///[s, e) inline int64_t detect_near_cc_tlen(bit_extz_t *ez, int64_t ck0, int64_t xk0, int64_t yk0, uint8_t rev) { @@ -17611,6 +18333,7 @@ inline int64_t detect_near_cc_tlen(bit_extz_t *ez, int64_t ck0, int64_t xk0, int ///[s, e) int64_t extract_sub_cigar_hc(overlap_region *z, All_reads *rref, haplotype_evdience_alloc* hp, char *qstr, uint64_t ql, UC_Read* tu, int64_t s, int64_t e, ul_ov_t *p, int64_t set_f, uint8_t *f, uint8_t occ_thres, uint64_t hpc_len/**, uint8_t is_dbg**/) { + // fprintf(stderr, "\n[M::%s]\ts::%ld\te::%ld\tset_f::%ld\tovlp_id::%u\twid::%u\n", __func__, s, e, set_f, ovlp_id(*p), ovlp_cur_wid(*p)); // if((!set_f) && (!ovlp_cur_ylen(*p))) return 1;///no potential informative site int64_t wk = ovlp_cur_wid(*p), xk = ovlp_cur_xoff(*p), yk = ovlp_cur_yoff(*p), ck = ovlp_cur_coff(*p), os, oe, t; bit_extz_t ez; int64_t bd = ovlp_bd(*p), s0, e0; char *ystr = NULL; @@ -17623,7 +18346,7 @@ int64_t extract_sub_cigar_hc(overlap_region *z, All_reads *rref, haplotype_evdie set_bit_extz_t(ez, (*z), wk); if(!ez.cigar.n) return -1; - int64_t cn = ez.cigar.n, op; int64_t ws, we, ovlp, xk0, yk0, ck0, yk1 = -1, /**xk1 = -1, ck1 = -1,**/ yl; haplotype_evdience ev; + int64_t cn = ez.cigar.n, op; int64_t ws, we, ovlp, xk0, yk0, ck0, yk1 = -1, /**xk1 = -1, ck1 = -1,**/ yl; haplotype_evdience ev; uint8_t om; xk0 = xk; yk0 = yk; ck0 = ck; ///for assertion if((ck < 0) || (ck > cn)) {//(*ck) == cn is allowed ck = 0; xk = ez.ts; yk = ez.ps; @@ -17677,19 +18400,20 @@ int64_t extract_sub_cigar_hc(overlap_region *z, All_reads *rref, haplotype_evdie } else { if(op == 0) { for (t = os; t < oe; t++) { - if(f[t-s] > occ_thres) { + if(f[t-s]) { + om = ((f[t-s]==3)?1:0); ev.misBase = qstr[t]; ev.overlapID = ovlp_id(*p); ev.site = t; ev.overlapSite = t-xk+yk; - ev.type = op; + ev.type = (om<<1); ev.cov = 1; addHaplotypeEvdience(hp, &ev, NULL); } } } else if(op == 1) { for (t = os; t < oe; t++) { - if(f[t-s] > occ_thres) { + if(f[t-s]) { if(!ystr) { // yk0 = ((hpc_len)?(detect_near_cc_tlen(&ez, ck, xk, yk, 1)):(t-xk+yk)); yk0 = t-xk+yk; @@ -17698,15 +18422,19 @@ int64_t extract_sub_cigar_hc(overlap_region *z, All_reads *rref, haplotype_evdie recover_UC_Read_sub_region(ystr, yk0, (yk1 - yk0), z->y_pos_strand, rref, z->y_id); } - if((!hpc_len) || (!hpc_mask_ff(ystr, yk1 - yk0, t-xk+yk-yk0, hpc_len, HPC_RR, NULL, -1, -1))) { - ev.misBase = ystr[t-xk+yk-yk0]; - ev.overlapID = ovlp_id(*p); - ev.site = t; - ev.overlapSite = t-xk+yk; - ev.type = op; - ev.cov = 1; - addHaplotypeEvdience(hp, &ev, NULL); - } + om = ((f[t-s]==3)?1:0); + if((!om) && (hpc_len) && (hpc_mask_ff(ystr, yk1 - yk0, t-xk+yk-yk0, hpc_len, HPC_RR, NULL, -1, -1))) om = 1; + + ev.misBase = ystr[t-xk+yk-yk0]; + ev.overlapID = ovlp_id(*p); + ev.site = t; + ev.overlapSite = t-xk+yk; + ev.type = (om<<1) + 1; + ev.cov = 1; + addHaplotypeEvdience(hp, &ev, NULL); + // if(ev.overlapID == 21 && ev.site == 17129) { + // fprintf(stderr, "[M::%s]\ts::%ld\te::%ld\tck::%ld\txk::%ld\tyk::%ld\tset_f::%ld\n", __func__, s, e, ck, xk, yk, set_f); + // } } } @@ -18453,6 +19181,7 @@ void rphase_hc_back(overlap_region_alloc* ol, All_reads *rref, haplotype_evdienc l0 = hp->length; // if(is_dbg) fprintf(stderr, "-1-[M::%s]\ts::%ld\te::%ld\n", __func__, s, e); rr = hc_phase_robust_rr(ol->list, rref, hp, qu->seq, qu->length, tu, idx->a + srt_n, idx->n - srt_n, s, e, c_idx->a, 1, occ_thres, hpc_len/**, is_dbg**/); + /** for (wi = fi = ei = 0, si = ((uint64_t)-1), wl0 = e - s; wi < wl0; wi++) { if(hp->flag[wi] > 0) { if((hp->flag[wi] > occ_thres) && ((!hpc_len) || (!hpc_mask_ff(qu->seq, qu->length, wi + s, hpc_len, HPC_RR, hp->flag, e - s, s)))) { @@ -18461,6 +19190,18 @@ void rphase_hc_back(overlap_region_alloc* ol, All_reads *rref, haplotype_evdienc ei = wi + 1; if(si == ((uint64_t)-1)) si = wi; } } + **/ + for (wi = fi = ei = 0, si = ((uint64_t)-1), wl0 = e - s; wi < wl0; wi++) { + if(hp->flag[wi] > 0) { + if(hp->flag[wi] > occ_thres) { + fi = 1; hp->nn_snp++; hp->flag[wi] = 1; + if((hpc_len) && (hpc_mask_ff(qu->seq, qu->length, wi + s, hpc_len, HPC_RR, NULL, 0, 0))) hp->flag[wi] = 3; + ei = wi + 1; if(si == ((uint64_t)-1)) si = wi; + } else { + hp->flag[wi] = 0; + } + } + } if(fi) { // if(is_dbg) fprintf(stderr, "-2-[M::%s]\ts::%ld\te::%ld\n", __func__, s, e); @@ -18481,7 +19222,7 @@ void rphase_hc_back(overlap_region_alloc* ol, All_reads *rref, haplotype_evdienc srt_n = hp->length; z = ((std_bs)?(ol->list):(NULL)); for (k = 1, i = t = 0; k <= srt_n; ++k) { if (k == srt_n || hp->list[k].site != hp->list[i].site) { - t += push_info(hp, hp->list+i, k-i, hp->list+t, z); + t += push_info(hp, hp->list+i, k-i, hp->list+t, z, NULL, sc_wn); i = k; } } @@ -18514,7 +19255,7 @@ void rphase_hc_back(overlap_region_alloc* ol, All_reads *rref, haplotype_evdienc } -void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, 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) +void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, 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) { int64_t on = ol->length, k, i, zwn, q[2]; uint64_t m, l0, wi, wl0, si, ei, fi; overlap_region *z; ul_ov_t *cp; @@ -18597,19 +19338,22 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all // debug_inter(ol, c_idx, idx->a, srt_n, idx->a + srt_n, idx->n - srt_n, s, e); l0 = hp->length; // if(is_dbg) fprintf(stderr, "-1-[M::%s]\ts::%ld\te::%ld\n", __func__, s, e); - rr = hc_phase_robust_rr(ol->list, rref, hp, qu->seq, qu->length, tu, idx->a + srt_n, idx->n - srt_n, s, e, c_idx->a, 1, occ_thres, 0/**hpc_len**//**, is_dbg**/); + rr = hc_phase_robust_rr(ol->list, rref, hp, qu->seq, qu->length, tu, idx->a + srt_n, idx->n - srt_n, s, e, c_idx->a, 1, occ_thres, hpc_len/**, is_dbg**/); for (wi = fi = ei = 0, si = ((uint64_t)-1), wl0 = e - s; wi < wl0; wi++) { if(hp->flag[wi] > 0) { - if((hp->flag[wi] > occ_thres)/** && ((!hpc_len) || (!hpc_mask_ff(qu->seq, qu->length, wi + s, hpc_len, HPC_RR, hp->flag, e - s, s)))**/) { - fi = 1; hp->nn_snp++; + if(hp->flag[wi] > occ_thres) { + fi = 1; hp->nn_snp++; hp->flag[wi] = 1; + if((hpc_len) && (hpc_mask_ff(qu->seq, qu->length, wi + s, hpc_len, HPC_RR, NULL, 0, 0))) hp->flag[wi] = 3; + ei = wi + 1; if(si == ((uint64_t)-1)) si = wi; + } else { + hp->flag[wi] = 0; } - ei = wi + 1; if(si == ((uint64_t)-1)) si = wi; - } + } } if(fi) { // if(is_dbg) fprintf(stderr, "-2-[M::%s]\ts::%ld\te::%ld\n", __func__, s, e); - rr = hc_phase_robust_rr(ol->list, rref, hp, qu->seq, qu->length, tu, idx->a + srt_n, idx->n - srt_n, s, e, c_idx->a, 0, occ_thres, 0/**hpc_len**//**, is_dbg**/); + rr = hc_phase_robust_rr(ol->list, rref, hp, qu->seq, qu->length, tu, idx->a + srt_n, idx->n - srt_n, s, e, c_idx->a, 0, occ_thres, hpc_len/**, is_dbg**/); if(hp->length > l0) radix_sort_haplotype_evdience_srt(hp->list + l0, hp->list + hp->length); } @@ -18623,20 +19367,25 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all SetSnpMatrix(hp, &(hp->nn_snp), &(ol->length), 0, NULL); - srt_n = hp->length; z = ((std_bs)?(ol->list):(NULL)); + srt_n = hp->length; z = (((std_bs)||(t8))?(ol->list):(NULL)); for (k = 1, i = t = 0; k <= srt_n; ++k) { if (k == srt_n || hp->list[k].site != hp->list[i].site) { - t += push_info(hp, hp->list+i, k-i, hp->list+t, z); + t += push_info(hp, hp->list+i, k-i, hp->list+t, z, t8, sc_wn); i = k; } } hp->length = t; - // generate_haplotypes_naive_advance(hap, overlap_list, NULL); - generate_haplotypes_naive_HiFi(hp, ol, 0.04, qu, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1))); - // 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); - + if(!dp) { + // generate_haplotypes_naive_advance(hap, overlap_list, NULL); + generate_haplotypes_naive_HiFi(hp, ol, 0.04, qu, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1))); + // 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 { + gen_rphase_dp(hp, ol, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), dp, idx, buf, rid, q8); + generate_haplotypes_naive_HiFi(hp, ol, 0.04, qu, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1))); + } + // lable_large_indels(overlap_list, g_read->length, dumy, asm_opt.max_ov_diff_ec); diff --git a/Correct.h b/Correct.h index 7dfea9f..33ceacf 100644 --- a/Correct.h +++ b/Correct.h @@ -144,7 +144,10 @@ typedef struct char misBase; }haplotype_evdience; - +#define hh_tp(z) (((z).type&1)) +#define hh_hp(z) ((((z).type>>1)&1)) +#define hh_bq(z) ((((z).type>>2))&sc_bm) +#define hh_wq(z) ((((z).type>>(sc_bn+2)))&sc_bm) typedef struct { @@ -1391,10 +1394,11 @@ bit_extz_t *exz, double e_rate, int64_t qs); void gen_hc_r_alin(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); void gen_hc_r_alin_nec(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); 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); -void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, 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); +void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, 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); 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); +void get_wqual(uint64_t zid, uint64_t zpos, uint64_t zrev, asg8_v *v, uint8_t *va, uint64_t scw, uint64_t *tqual, uint64_t *wqual); #define ovlp_id(x) ((x).tn) diff --git a/Hash_Table.h b/Hash_Table.h index c40b672..8fc79ac 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -10,6 +10,7 @@ #define WINDOW_BOUNDARY 375 #define WINDOW_HC 775 #define WINDOW_OHC 475 +// #define WINDOW_OHC 375 ///ONT high error #define WINDOW_HC_FAST 512 ///for one side, the first or last WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY bases should not be corrected #define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25 diff --git a/Process_Read.cpp b/Process_Read.cpp index 8c17fa5..c95b7cc 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -51,6 +51,7 @@ void destory_All_reads(All_reads* r) if (r->read_sperate[i]) free(r->read_sperate[i]); if (r->paf && r->paf[i].buffer) free(r->paf[i].buffer); if (r->reverse_paf && r->reverse_paf[i].buffer) free(r->reverse_paf[i].buffer); + if(r->rsc && r->rsc[i]) free(r->rsc[i]); ///if (r->pb_regions) kv_destroy(r->pb_regions[i].a); } free(r->paf); @@ -61,6 +62,7 @@ void destory_All_reads(All_reads* r) free(r->name_index); free(r->read_length); free(r->trio_flag); + free(r->rsc); ///if (r->pb_regions) free(r->pb_regions); } @@ -108,6 +110,14 @@ void write_All_reads(All_reads* r, char* read_file_name) fwrite(&(asm_opt.hom_cov), sizeof(asm_opt.hom_cov), 1, fp); fwrite(&(asm_opt.het_cov), sizeof(asm_opt.het_cov), 1, fp); + uint64_t mm = 1; + if(asm_opt.is_sc) { + fwrite(&mm, sizeof(mm), 1, fp); + for (i = 0; i < r->total_reads; i++) { + fwrite(r->rsc[i], sizeof(uint8_t), ((r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0)), fp); + } + } + free(index_name); fflush(fp); fclose(fp); @@ -201,6 +211,19 @@ int load_All_reads(All_reads* r, char* read_file_name) } ///r->pb_regions = NULL; + uint64_t mm = 0; + if (!feof(fp)) { + if((fread(&mm, sizeof(mm), 1, fp)) && (mm == 1)) { + MALLOC(r->rsc, r->total_reads); + for (i = 0; i < r->total_reads; i++) { + MALLOC(r->rsc[i], (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0)); + f_flag += fread(r->rsc[i], sizeof(uint8_t), (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0), fp); + } + } + } + + + free(index_name); fclose(fp); fprintf(stderr, "Reads has been loaded.\n"); @@ -412,10 +435,13 @@ void malloc_All_reads(All_reads* r) memcpy(r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads); r->read_sperate = (uint8_t**)malloc(sizeof(uint8_t*)*r->total_reads); + if(asm_opt.is_sc) MALLOC(r->rsc, r->total_reads); + long long i = 0; for (i = 0; i < (long long)r->total_reads; i++) { r->read_sperate[i] = (uint8_t*)malloc(sizeof(uint8_t)*(r->read_length[i]/4+1)); + if(r->rsc) MALLOC(r->rsc[i], (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0)); } r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads); @@ -823,6 +849,131 @@ void ha_compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_sit } } +void convert_qual(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitu, uint64_t rev, uint64_t sc_off) +{ + uint64_t i = 0; uint8_t c = 0, sc; + // fprintf(stderr, "\n[M::%s]\n", __func__); + for (i = 0; i < src_l; i++) { + for (c = 0, sc = ((uint8_t)src[i]) - sc_off; (c < bitu) && (sc_tb[c] < sc); c++); + if(c >= bitu) c = bitu - 1; + dest[(rev?(src_l-i-1):(i))] = c; + // fprintf(stderr, "%u->%u\n", sc, c); + } +} + +void ha_compress_qual_bit(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitn) +{ + + uint64_t i = 0, k, bit_r = 8/bitn, dest_i = 0; + uint8_t tmp = 0, c = 0; + + for (i = 0; i + bit_r <= src_l;) { + for (k = tmp = 0; k < bit_r; k++) { + c = ((uint8_t)src[i]); + tmp <<= bitn; tmp |= c; i++; + } + dest[dest_i++] = tmp; + } + + if(i < src_l) { + for (k = tmp = 0; i < src_l; k++) { + c = ((uint8_t)src[i]); + tmp <<= bitn; tmp |= c; i++; + } + + dest[dest_i++] = (tmp<<(8-(bitn*k))); + } +} + +void ha_compress_qual(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitn, uint64_t sc_off) +{ + + uint64_t i = 0, k, bit_r = 8/bitn, dest_i = 0, bitu = (1<= bitu) c = bitu - 1; + tmp <<= bitn; tmp |= c; i++; + } + dest[dest_i++] = tmp; + } + + if(i < src_l) { + for (k = tmp = 0; i < src_l; k++) { + for (c = 0, sc = ((uint8_t)src[i]) - sc_off; (c < bitu) && (sc_tb[c] < sc); c++); + if(c >= bitu) c = bitu - 1; + tmp <<= bitn; tmp |= c; i++; + } + + dest[dest_i++] = (tmp<<(8-(bitn*k))); + } +} + +///[s, e) +int64_t retrive_bqual(asg8_v *dv, uint8_t *ds, uint64_t id, int64_t s, int64_t e, uint8_t rev, int64_t bitn) +{ + int64_t rl = Get_READ_LENGTH(R_INF, id), l; + if(s < 0) s = 0; if(e < 0) e = rl; + if(s >= e || e > rl) return -1; + + uint8_t *da = NULL, *src = Get_QUAL(R_INF, id), mm = (((uint8_t)1)<a; + } else { + da = ds; + } + + if(!rev) { + dk = 0; sk = s; + + mrf = ((s%bitr)*bitn); + // if(s == 21519 && e == 22332) { + // fprintf(stderr, "+[M::%s] id::%lu, in::[%ld, %ld), rev::%u, bitn::%ld, bitr::%ld, mrf::%u\n", __func__, id, s, e, rev, bitn, bitr, mrf); + // } + if(mrf) { + for (swk = sk/bitr; mrf < 8 && sk < e; mrf += bitn, sk++) da[dk++] = ((src[swk]<>mlf)&mm; + } + + for (swk = sk/bitr; (sk + bitr) <= e; sk += bitr, swk++) { + for (mrf = 0; mrf < 8; mrf += bitn) da[dk++] = ((src[swk]<>mlf)&mm; + } + + if(sk < e) { + for (mrf = 0; sk < e; mrf += bitn, sk++) da[dk++] = ((src[swk]<>mlf)&mm; + } + // if(dk != l) { + // fprintf(stderr, "+[M::%s] id::%lu, in::[%ld, %ld), rev::%u, bitn::%ld, bitr::%ld\n", __func__, id, s, e, rev, bitn, bitr); + // } + assert(dk == l); + } else { + sk = s; s = e; e = sk; + s = rl - s; e = rl - e; + dk = l; sk = s; + + mrf = ((s%bitr)*bitn); + if(mrf) { + for (swk = sk/bitr; mrf < 8 && sk < e; mrf += bitn, sk++) da[--dk] = ((src[swk]<>mlf)&mm; + } + + for (swk = sk/bitr; (sk + bitr) <= e; sk += bitr, swk++) { + for (mrf = 0; mrf < 8; mrf += bitn) da[--dk] = ((src[swk]<>mlf)&mm; + } + + if(sk < e) { + for (mrf = 0; sk < e; mrf += bitn, sk++) da[--dk] = ((src[swk]<>mlf)&mm; + } + + assert(dk == 0); + } + dv->n = l; + + return l; +} + void reverse_complement(char* pattern, uint64_t length) { uint64_t i = 0; @@ -845,6 +996,24 @@ void reverse_complement(char* pattern, uint64_t length) } } +void print_fastq(FILE *fp, char *id, char *bs, char *qual, uint64_t bitu, uint64_t sc_off) +{ + uint64_t i = 0, ql = strlen(qual); uint8_t c = 0, sc; + + if(fp) fprintf(fp, "@%s\n%s\n+\n", id, bs); + else fprintf(stdout, "@%s\n%s\n+\n", id, bs); + + for (i = 0; i < ql; i++) { + for (c = 0, sc = ((uint8_t)qual[i]) - sc_off; (c < bitu) && (sc_tb[c] < sc); c++); + if(c >= bitu) c = bitu - 1; + if(fp) fprintf(fp, "%u", c); + else fprintf(stdout, "%u", c); + } + + if(fp) fprintf(fp, "\n"); + else fprintf(stdout, "\n"); +} + void init_Debug_reads(Debug_reads* x, const char* file) { diff --git a/Process_Read.h b/Process_Read.h index 4711b24..99d7e30 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -22,6 +22,7 @@ #define Get_NAME_LENGTH(R_INF, ID) ((R_INF).name_index[(ID)+1] - (R_INF).name_index[(ID)]) ///#define Get_READ(R_INF, ID) R_INF.read + (R_INF.index[ID]>>2) + ID #define Get_READ(R_INF, ID) (R_INF).read_sperate[(ID)] +#define Get_QUAL(R_INF, ID) (R_INF).rsc[(ID)] #define Get_NAME(R_INF, ID) ((R_INF).name + (R_INF).name_index[(ID)]) #define CHECK_BY_NAME(R_INF, NAME, ID) (Get_NAME_LENGTH((R_INF),(ID))==strlen((NAME)) && \ memcmp((NAME), Get_NAME((R_INF), (ID)), Get_NAME_LENGTH((R_INF),(ID))) == 0) @@ -38,6 +39,8 @@ extern char rc_Table[6]; void init_aux_table(); +typedef struct { size_t n, m; uint8_t *a; } asg8_v; + typedef struct { uint64_t x_id; @@ -119,6 +122,7 @@ typedef struct uint64_t* read_length; uint64_t* read_size; uint8_t* trio_flag; + uint8_t** rsc; ///seq start pos in uint8_t* read ///do not need it @@ -223,6 +227,7 @@ void init_All_reads(All_reads* r); void malloc_All_reads(All_reads* r); void ha_insert_read_len(All_reads *r, int read_len, int name_len); void ha_compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ); +void ha_compress_qual(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitn, uint64_t sc_off); void init_UC_Read(UC_Read* r); void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID); void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID); @@ -254,4 +259,17 @@ scaf_res_t *init_scaf_res_t(uint32_t n); void destroy_scaf_res_t(scaf_res_t *p); void read_ma(ma_hit_t* x, FILE* fp); +const uint64_t sc_tb[8] = { + 10, 20, 30, 40, 50, 60, 70, 80 +}; + +#define sc_bn 2 +#define sc_bm ((((uint64_t)1)<v64); kv_init(z->v32); kv_init(z->v16); + kv_init(z->v8q); + kv_init(z->v8t); init_bit_extz_t(&(z->exz), 31); z->ab = ha_abuf_init(); @@ -117,6 +119,8 @@ void destroy_ec_ovec_buf_t(ec_ovec_buf_t *p) kv_destroy(z->v64); kv_destroy(z->v32); kv_destroy(z->v16); + kv_destroy(z->v8q); + kv_destroy(z->v8t); destroy_bit_extz_t(&(z->exz)); ha_abuf_destroy(z->ab); @@ -148,6 +152,8 @@ inline void refresh_ec_ovec_buf_t0(ec_ovec_buf_t0 *z, uint64_t n) kv_destroy(z->v64); kv_init(z->v64); kv_destroy(z->v32); kv_init(z->v32); kv_destroy(z->v16); kv_init(z->v16); + kv_destroy(z->v8q); kv_init(z->v8q); + kv_destroy(z->v8t); kv_init(z->v8t); destroy_bit_extz_t(&(z->exz)); init_bit_extz_t(&(z->exz), 31); @@ -2930,6 +2936,56 @@ void dedup_chains(overlap_region_alloc* ol) } } +void debug_retrive_bqual(asg8_v *vq, asg8_v *vt, uint64_t id, uint64_t rn) +{ + uint64_t k, n, z[2], s, e, rev; + retrive_bqual(vq, NULL, id, -1, -1, 0, sc_bn); n = vq->n; + retrive_bqual(vt, NULL, id, -1, -1, 1, sc_bn); + assert(vq->n == vt->n); + for (k = 0; k < vq->n && vq->a[k] == vt->a[vt->n - k - 1]; k++); + // if((k == vq->n)) { + // fprintf(stderr, "[M::%s] id::%lu, k::%lu, n::%lu\n", __func__, id, k, ((uint64_t)vq->n)); + // } + assert(k == vq->n); + + // if(id == 0) { + // s = 21519; e = 22332; rev = 0; + // retrive_bqual(vt, NULL, id, s, e, rev, sc_bn); + // if(memcmp(vq->a + s, vt->a, e - s)) { + // fprintf(stderr, "+[M::%s] id::%lu, t::[%lu, %lu), rev::%lu\n", __func__, id, s, e, rev); + // exit(1); + // } + + // } + // return; + + for (k = 0, rev = 0; k < rn; k++) { + z[0] = rand()%(n + 1); + z[1] = rand()%(n + 1); + if(z[0] == z[1]) continue; + s = MIN(z[0], z[1]); e = MAX(z[0], z[1]); + retrive_bqual(vt, NULL, id, s, e, rev, sc_bn); + if(memcmp(vq->a + s, vt->a, e - s)) { + fprintf(stderr, "[M::%s] id::%lu, t::[%lu, %lu), rev::%lu\n", __func__, id, s, e, rev); + exit(1); + } + } + + + retrive_bqual(vq, NULL, id, -1, -1, 1, sc_bn); + for (k = 0, rev = 1; k < rn; k++) { + z[0] = rand()%(n + 1); + z[1] = rand()%(n + 1); + if(z[0] == z[1]) continue; + s = MIN(z[0], z[1]); e = MAX(z[0], z[1]); + retrive_bqual(vt, NULL, id, s, e, rev, sc_bn); + if(memcmp(vq->a + s, vt->a, e - s)) { + fprintf(stderr, "[M::%s] id::%lu, t::[%lu, %lu), rev::%lu\n", __func__, id, s, e, rev); + exit(1); + } + } +} + static void worker_hap_ec(void *data, long i, int tid) { ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); @@ -2948,7 +3004,9 @@ static void worker_hap_ec(void *data, long i, int tid) **/ // if(i < 1100000 || i > 1400000) return; // if(i % 100000 == 0) fprintf(stderr, "-a-[M::%s-beg] rid->%ld\n", __func__, i); - // if (memcmp("4da034b0-a94d-4576-8481-c0d9a9f96d40", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { + // if (memcmp("c42804f3-0e13-43a0-8a71-b91b40accf9a", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { + // if (memcmp("b2e68ecf-381a-439c-b676-c1e6831d6acf", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { + // if (memcmp("0aec8c4f-c849-4c31-85ba-4ffb297eeb28", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { // fprintf(stderr, "-a-[M::%s-beg] rid->%ld\n", __func__, i); // } else { // return; @@ -2961,9 +3019,11 @@ static void worker_hap_ec(void *data, long i, int tid) // if(i != 2243244) return; // if(i != 19350) return; + // debug_retrive_bqual(D, &b->v8t, i, 256); return; + recover_UC_Read(&b->self_read, &R_INF, i); qlen = b->self_read.length; - h_ec_lchain(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, 0.02, asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32); + h_ec_lchain(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, /**((asm_opt.is_ont)?(0.05):(0.02))**/0.02, asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32);///ONT high error // b->num_read_base += b->olist.length; b->cnt[0] += b->self_read.length; @@ -2987,7 +3047,7 @@ static void worker_hap_ec(void *data, long i, int tid) // b->num_correct_base += b->olist.length; copy_asg_arr(buf0, b->sp); - rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &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); + rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &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)); copy_asg_arr(b->sp, buf0); // stderr_phase_ovlp(&b->olist); @@ -3495,6 +3555,12 @@ static void worker_hap_post_rev(void *data, long i, int tid) } ha_compress_base(Get_READ(R_INF, i), a, l, &R_INF.N_site[i], nn); + + if(asm_opt.is_sc) { + retrive_bqual(&(b->v8q), NULL, i, -1, -1, 0, sc_bn); + for (k = 0; k < l; k++) a[l - k - 1] = b->v8q.a[k]; + ha_compress_qual_bit(Get_QUAL(R_INF, i), a, l, sc_bn); + } } static void worker_hap_dc_ec_gen(void *data, long i, int tid) @@ -5018,7 +5084,7 @@ static void worker_hap_dc_ec0(void *data, long i, int tid) b->cnt[0] += b->self_read.length; copy_asg_arr(buf0, b->sp); - rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &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); + rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &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), ((asm_opt.is_sc)?&(b->v8t):NULL)); copy_asg_arr(b->sp, buf0); copy_asg_arr(buf0, b->sp); @@ -5165,7 +5231,7 @@ static void worker_sl_ec(void *data, long i, int tid) { // if(i != 0) return; - sl_v *p = &(((sl_v*)data)[tid]); + sl_v *p = &(((sl_v*)data)[tid]); uint8_t *oa = NULL; char *na = NULL; uint64_t tqual, wqual; uint32_t ci = 0, len, xk, yk, wx[2], wy[2], k, Nn, yn = 0, tot_e; uint16_t c, bq, bt; @@ -5213,11 +5279,13 @@ static void worker_sl_ec(void *data, long i, int tid) } // if(i == 700) fprintf(stderr, "|\n"); + if(asm_opt.is_sc) retrive_bqual(&(p->q), NULL, i, -1, -1, 0, sc_bn); if (R_INF.read_size[i] < yn) { R_INF.read_size[i] = yn; REALLOC(R_INF.read_sperate[i], R_INF.read_size[i]/4+1); + if(asm_opt.is_sc) REALLOC(R_INF.rsc[i], ((R_INF.read_size[i]/sc_bn) + ((R_INF.read_size[i]%sc_bn)?1:0))); } R_INF.read_length[i] = yn; // if(Nn > 0) fprintf(stderr, "[M::%s] Nn->%u\n", __func__, Nn); @@ -5236,7 +5304,25 @@ static void worker_sl_ec(void *data, long i, int tid) ha_compress_base(Get_READ(R_INF, i), p->a, yn, &R_INF.N_site[i], Nn); - + if(asm_opt.is_sc) { + oa = p->q.a; na = p->a; + ci = 0; xk = yk = 0; Nn = 0; + while (ci < scc.a[i].n) { + wx[0] = xk; wy[0] = yk; + ci = pop_trace_bp_f(&scc.a[i], ci, &c, &bq, &bt, &len); + if(c != 2) xk += len; + if(c != 3) yk += len; + wx[1] = xk; wy[1] = yk; + if(c == 0 || c == 1) { + memcpy(na + wy[0], oa + wx[0], (wx[1]-wx[0])*sizeof((*oa))); + } else if(c == 2) { + get_wqual(i, wx[0], 0, NULL, oa, sc_wn, &tqual, &wqual); + for (k = wy[0]; k < wy[1]; k++) na[k] = wqual; + } + } + assert(yk == yn); + ha_compress_qual_bit(Get_QUAL(R_INF, i), na, yn, sc_bn); + } } uint64_t cal_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, uint64_t *r_base) @@ -5496,7 +5582,7 @@ void sl_ec_r(uint64_t n_thre, uint64_t n_a) sl_v *b = NULL; uint64_t k; MALLOC(b, n_thre); for (k = 0; k < n_thre; k++) { b[k].a = NULL; b[k].n = b[k].m = 0; - init_UC_Read(&b[k].z); + init_UC_Read(&b[k].z); kv_init(b[k].q); } kt_for(n_thre, worker_sl_ec, b, n_a);///debug_for_fix diff --git a/ecovlp.h b/ecovlp.h index fe5cc89..06467e3 100644 --- a/ecovlp.h +++ b/ecovlp.h @@ -42,6 +42,7 @@ typedef struct { asg64_v v64; asg32_v v32; asg16_v v16; + asg8_v v8q, v8t; kvec_t_u8_warp k_flag; st_mt_t sp; diff --git a/htab.cpp b/htab.cpp index c061bf0..eb94fdd 100644 --- a/htab.cpp +++ b/htab.cpp @@ -692,6 +692,7 @@ static inline void sf##_pt_insert_buf(sf##_ch_buf_t *buf, int p, const HType *y) static void *sf##_worker_count(void *data, int step, void *in) /** callback for kt_pipeline()**/\ {\ pl_data_t *p = (pl_data_t*)data;\ + /**uint8_t src_a[1000000], des_a[1000000];**/\ if (step == 0) { /** step 1: read a block of sequences**/\ int ret;\ sf##_st_data_t *s;\ @@ -762,6 +763,16 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for ++n_N;\ ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s+p->opt->adaLen, l, &p->rs_out->N_site[p->n_seq], n_N);\ memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);\ + if(p->rs_out->rsc) {\ + ha_compress_qual(Get_QUAL(*p->rs_out, p->n_seq), p->ks->qual.s+p->opt->adaLen, l, sc_bn, 33);\ + /**print_fastq(NULL, p->ks->name.s, p->ks->seq.s, p->ks->qual.s, (1<ks->qual.s+p->opt->adaLen, l, (1<n_seq, -1, -1, 0, sc_bn);\ + if(memcmp(src_a, des_a, l)!=0) fprintf(stderr, "ERROR: incorrect qual values\n");\ + else fprintf(stderr, "Correct: correct qual values\n");\ + }**/\ + }\ }\ }\ if (s->n_seq == s->m_seq) {\