diff --git a/CommandLines.h b/CommandLines.h index 8b2830e..dee0ca4 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.16.7-r418" +#define HA_VERSION "0.16.8-r420" #define VERBOSE 0 diff --git a/Correct.cpp b/Correct.cpp index 1ea4c1f..c312f07 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -3282,7 +3282,11 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea y_id = z->y_id; y_strand = z->y_pos_strand; y_readLen = (rref?(Get_READ_LENGTH((*rref), y_id)):(uref->ug->u.a[y_id].len)); overlap_length = z->x_pos_e + 1 - z->x_pos_s; //z->is_match = 0; - + // if(y_id == 0 || y_id == 1) { + // fprintf(stderr, "[M::%s::j->%lld] utg%.6dl(%c), align_length::%u, overlap_length::%lld\n", __func__, + // j, (int32_t)z->y_id + 1, "+-"[z->y_pos_strand], z->align_length, overlap_length); + // } + ///debug_scan_cigar(&(overlap_list->list[j])); ///only calculate cigar for high quality overlaps if ((rref && (overlap_length*OVERLAP_THRESHOLD_FILTER <= z->align_length)) || @@ -7502,11 +7506,33 @@ void generate_haplotypes_naive_advance(haplotype_evdience_alloc* hap, overlap_re 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].overlapSite == 55) { + // fprintf(stderr, "5555555555[M::%s::] utg%.6dl(%c), x_site::%u, x::[%u, %u), misBase::%c, occ_0::%u, occ_1::%u, occ_2::%u, snp_idx::%u\n", __func__, + // ((int32_t)(overlap_list->list[hap->list[i].overlapID].y_id)) + 1, + // "+-"[overlap_list->list[hap->list[i].overlapID].y_pos_strand], + // hap->list[i].site, overlap_list->list[hap->list[i].overlapID].x_pos_s, + // overlap_list->list[hap->list[i].overlapID].x_pos_e + 1, hap->list[i].misBase, + // s->occ_0, s->occ_1, s->occ_2, hap->list[i].overlapSite); + // } if(hap->list[i].type!=1) continue;///mismatch s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + // if(hap->list[i].overlapID == 125 || hap->list[i].overlapID == 127) { + // fprintf(stderr, "[M::%s::] utg%.6dl(%c), x_site::%u, x::[%u, %u), misBase::%c, occ_0::%u, occ_1::%u, occ_2::%u, snp_idx::%u\n", __func__, + // ((int32_t)(overlap_list->list[hap->list[i].overlapID].y_id)) + 1, + // "+-"[overlap_list->list[hap->list[i].overlapID].y_pos_strand], + // hap->list[i].site, overlap_list->list[hap->list[i].overlapID].x_pos_s, + // overlap_list->list[hap->list[i].overlapID].x_pos_e + 1, hap->list[i].misBase, + // s->occ_0, s->occ_1, s->occ_2, hap->list[i].overlapSite); + // } if(s->occ_0 < 2 || s->occ_1 < 2) continue; if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) o++;///allels must be real } + + // if(hap->list[l].overlapID == 125 || hap->list[l].overlapID == 127) { + // fprintf(stderr, "+++[M::%s::] utg%.6dl(%c), o::%lu\n", __func__, + // ((int32_t)(overlap_list->list[hap->list[l].overlapID].y_id)) + 1, + // "+-"[overlap_list->list[hap->list[l].overlapID].y_pos_strand], o); + // } if(o > 0) { o = ((uint32_t)-1) - o; o <<= 32; o += l; @@ -7528,11 +7554,11 @@ void generate_haplotypes_naive_advance(haplotype_evdience_alloc* hap, overlap_re if(s->occ_0 < 2 || s->occ_1 < 2) continue; if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) o++; } - // fprintf(stderr, "o: %lu, op: %lu, %c, ovlp_id: %u, x_s: %u, x_e: %u\n", - // o, ((uint32_t)-1) - (hap->snp_srt.a[k]>>32), - // "ze"[o==(((uint32_t)-1) - (hap->snp_srt.a[k]>>32))], hap->list[l].overlapID, - // overlap_list->list[hap->list[l].overlapID].x_pos_s, - // overlap_list->list[hap->list[l].overlapID].x_pos_e); + // if(hap->list[l].overlapID == 125 || hap->list[l].overlapID == 127) { + // fprintf(stderr, "sbsbsb[M::%s::] utg%.6dl(%c), o::%lu\n", __func__, + // ((int32_t)(overlap_list->list[hap->list[l].overlapID].y_id)) + 1, + // "+-"[overlap_list->list[hap->list[l].overlapID].y_pos_strand], o); + // } if(o == 0) continue; ii = hap->list[l].overlapID; @@ -7546,9 +7572,8 @@ void generate_haplotypes_naive_advance(haplotype_evdience_alloc* hap, overlap_re for (z = hap->list[i].overlapSite; z >= 0; z--) { t = &(hap->snp_stat.a[z]); if(s->site!=t->site) break; - // t->occ_0--; t->occ_0 -= hap->list[i].cov; - if(t->occ_0 < 1) fprintf(stderr, "WRONG-CORRECTION\n"); + assert(t->occ_0 >= 1);// if(t->occ_0 < 1) fprintf(stderr, "WRONG-CORRECTION\n"); } } } @@ -7560,9 +7585,24 @@ void generate_haplotypes_naive_advance(haplotype_evdience_alloc* hap, overlap_re if(hap->list[i].type!=1) continue; s = &(hap->snp_stat.a[hap->list[i].overlapSite]); if(s->occ_0 < 2 || s->occ_1 < 2) continue; - if(s->score == 1) o++; + if(s->score == 1) { + o++; + // if(hap->list[i].overlapID == 125 || hap->list[i].overlapID == 127) { + // fprintf(stderr, "[M::%s::] utg%.6dl(%c), x_site::%u, x::[%u, %u), misBase::%c, occ_0::%u, occ_1::%u, occ_2::%u, snp_idx::%u\n", __func__, + // ((int32_t)(overlap_list->list[hap->list[i].overlapID].y_id)) + 1, + // "+-"[overlap_list->list[hap->list[i].overlapID].y_pos_strand], + // hap->list[i].site, overlap_list->list[hap->list[i].overlapID].x_pos_s, + // overlap_list->list[hap->list[i].overlapID].x_pos_e + 1, hap->list[i].misBase, + // s->occ_0, s->occ_1, s->occ_2, hap->list[i].overlapSite); + // } + } } ii = hap->list[l].overlapID; + // if(hap->list[l].overlapID == 125 || hap->list[l].overlapID == 127) { + // fprintf(stderr, ">>>[M::%s::] utg%.6dl(%c), o::%lu\n", __func__, + // ((int32_t)(overlap_list->list[hap->list[l].overlapID].y_id)) + 1, + // "+-"[overlap_list->list[hap->list[l].overlapID].y_pos_strand], o); + // } if(overlap_list->list[ii].is_match == 2 && o == 0) { overlap_list->list[ii].is_match = 1; } @@ -7588,6 +7628,11 @@ void generate_haplotypes_naive_advance(haplotype_evdience_alloc* hap, overlap_re if(!km) kv_push(uint64_t, hap->snp_srt, hap->list[i].overlapSite); else kv_push_km(km, uint64_t, hap->snp_srt, hap->list[i].overlapSite); } + // if(hap->list[l].overlapID == 125 || hap->list[l].overlapID == 127) { + // fprintf(stderr, "---[M::%s::] utg%.6dl(%c), o::%lu\n", __func__, + // ((int32_t)(overlap_list->list[hap->list[l].overlapID].y_id)) + 1, + // "+-"[overlap_list->list[hap->list[l].overlapID].y_pos_strand], o); + // } hap->snp_srt.n -= o; if(o >= 2) {///there are at least two variants at one read radix_sort_bc64(hap->snp_srt.a + hap->snp_srt.n, hap->snp_srt.a + hap->snp_srt.n + o); @@ -7651,6 +7696,397 @@ void generate_haplotypes_naive_advance(haplotype_evdience_alloc* hap, overlap_re } + +void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, void *km) +{ + // fprintf(stderr, "[M::%s::] Done\n", __func__); + if(hap->length == 0) return; + uint64_t k, l, i, o, *a, ii, m_snp_stat, m_list, m_off; + int64_t z; + SnpStats *s = NULL, *t = NULL; + + for (k = 1, l = 0, i = m_snp_stat = m_list = 0; k <= hap->snp_stat.n; ++k) {///filter snps + if(k == hap->snp_stat.n || hap->snp_stat.a[k].site != hap->snp_stat.a[l].site) { + if((l > 0) && (hap->snp_stat.a[l].site == (hap->snp_stat.a[l-1].site + 1))) { + l = k; continue; + } + if((k < hap->snp_stat.n) && ((hap->snp_stat.a[l].site+1) == hap->snp_stat.a[k].site)) { + l = k; continue; + } + + for (; i < hap->length && hap->list[i].site != hap->snp_stat.a[l].site; i++); + assert(i < hap->length && hap->list[i].site == hap->snp_stat.a[l].site); + m_off = l - m_snp_stat; + for (; i < hap->length && hap->list[i].site == hap->snp_stat.a[l].site; i++) { + assert(hap->list[i].overlapSite>=l && hap->list[i].overlapSite < k); + // assert(hap->snp_stat.a[hap->list[i].overlapSite].site==hap->list[i].site); + hap->list[m_list] = hap->list[i]; hap->list[m_list++].overlapSite -= m_off; + } + + for (; l < k; l++) hap->snp_stat.a[m_snp_stat++] = hap->snp_stat.a[l]; + } + } + hap->snp_stat.n = m_snp_stat; hap->length = m_list; + if(hap->snp_stat.n == 0 || hap->length == 0) return; + + hap->snp_srt.n = 0; + radix_sort_haplotype_evdience_id_srt(hap->list, hap->list + hap->length); + 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 + 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; + if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) o++;///allels must be real + } + if(o > 0) { + o = ((uint32_t)-1) - o; + o <<= 32; o += l; + if(!km) kv_push(uint64_t, hap->snp_srt, o); + else kv_push_km(km, uint64_t, hap->snp_srt, o); + } + l = k; + } + } + + if (hap->snp_srt.n > 0) { + radix_sort_bc64(hap->snp_srt.a, hap->snp_srt.a + hap->snp_srt.n);///sort by how many snps in one overlap + 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; + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->occ_0 < 2 || s->occ_1 < 2) continue; + if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) o++; + } + if(o == 0) continue; + + 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){ + 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) { + ///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; + t->occ_0 -= hap->list[i].cov; + assert(t->occ_0 >= 1); + } + } + } + } + + 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; + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->occ_0 < 2 || s->occ_1 < 2) continue; + if(s->score == 1) o++; + } + ii = hap->list[l].overlapID; + ///for HiFi, do not flip trans to cis + // if(overlap_list->list[ii].is_match == 2 && o == 0) { + // overlap_list->list[ii].is_match = 1; + // } + if(overlap_list->list[ii].is_match == 1 && o > 0) { + overlap_list->list[ii].is_match = 2; + } + } + + + for (k = 1, l = 0; k <= hap->length; ++k) { ///reset snp_stat + if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { + 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) { + hap->snp_stat.a[hap->list[i].overlapSite].score = -1; + } + } + } + l = k; + } + } + } + + + + hap->snp_srt.n = 0; + for (k = 1, l = 0; k <= hap->length; ++k) { + if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { + if(overlap_list->list[hap->list[l].overlapID].is_match == 2) { + l = k; + continue; + } + for (i = l, o = 0; i < k; i++) { + if(hap->list[i].type!=1) continue; + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->occ_0 < 2 || s->occ_1 < 2) continue; + if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) continue; + if(s->score == 1) continue; + o++; + if(!km) kv_push(uint64_t, hap->snp_srt, hap->list[i].overlapSite); + else kv_push_km(km, uint64_t, hap->snp_srt, hap->list[i].overlapSite); + } + hap->snp_srt.n -= o; + ///there are at least two variants at one read + if(o>=(overlap_list->list[hap->list[l].overlapID].align_length*up)) { + radix_sort_bc64(hap->snp_srt.a + hap->snp_srt.n, hap->snp_srt.a + hap->snp_srt.n + o); + a = hap->snp_srt.a + hap->snp_srt.n; + for (i = z = 0; i < o; i++) { + if(i > 0) s = &(hap->snp_stat.a[a[i-1]]); + if(i + 1 < o) t = &(hap->snp_stat.a[a[i+1]]); + if(s && s->site + 32 > hap->snp_stat.a[a[i]].site) continue; + if(t && hap->snp_stat.a[a[i]].site + 32 > t->site) continue; + a[z] = a[i]; + z++; + } + if(z >= 2) hap->snp_srt.n += z; + } + l = k; + } + } + if (hap->snp_srt.n > 0) { + radix_sort_bc64(hap->snp_srt.a, hap->snp_srt.a + hap->snp_srt.n); + for (k = 1, l = 0; k <= hap->snp_srt.n; ++k) { + if(k == hap->snp_srt.n || hap->snp_srt.a[k] != hap->snp_srt.a[l]) { + if(k - l >= 2) hap->snp_stat.a[hap->snp_srt.a[l]].score = 1; + } + l = k; + } + } + + 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(hap->list[i].type==1 || hap->list[i].type==0) { + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->score == 1 && (!(s->occ_0 < 2 || s->occ_1 < 2))) { + overlap_list->list[ii].strong = 1; + if(hap->list[i].type==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; + } + } +} + + +void generate_haplotypes_naive_UL(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, void *km) +{ + if(hap->length == 0) return; + uint64_t k, l, i, o, *a, ii, m_snp_stat, m_list, m_off; + int64_t z; + SnpStats *s = NULL, *t = NULL; + + for (k = 1, l = 0, i = m_snp_stat = m_list = 0; k <= hap->snp_stat.n; ++k) {///filter snps + if(k == hap->snp_stat.n || hap->snp_stat.a[k].site != hap->snp_stat.a[l].site) { + if((l > 0) && (hap->snp_stat.a[l].site == (hap->snp_stat.a[l-1].site + 1))) { + l = k; continue; + } + if((k < hap->snp_stat.n) && ((hap->snp_stat.a[l].site+1) == hap->snp_stat.a[k].site)) { + l = k; continue; + } + + for (; i < hap->length && hap->list[i].site != hap->snp_stat.a[l].site; i++); + assert(i < hap->length && hap->list[i].site == hap->snp_stat.a[l].site); + m_off = l - m_snp_stat; + for (; i < hap->length && hap->list[i].site == hap->snp_stat.a[l].site; i++) { + assert(hap->list[i].overlapSite>=l && hap->list[i].overlapSite < k); + // assert(hap->snp_stat.a[hap->list[i].overlapSite].site==hap->list[i].site); + hap->list[m_list] = hap->list[i]; hap->list[m_list++].overlapSite -= m_off; + } + + for (; l < k; l++) hap->snp_stat.a[m_snp_stat++] = hap->snp_stat.a[l]; + } + } + hap->snp_stat.n = m_snp_stat; hap->length = m_list; + if(hap->snp_stat.n == 0 || hap->length == 0) return; + + hap->snp_srt.n = 0; + radix_sort_haplotype_evdience_id_srt(hap->list, hap->list + hap->length); + 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 + 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; + if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) o++;///allels must be real + } + if(o > 0) { + o = ((uint32_t)-1) - o; + o <<= 32; o += l; + if(!km) kv_push(uint64_t, hap->snp_srt, o); + else kv_push_km(km, uint64_t, hap->snp_srt, o); + } + l = k; + } + } + + if (hap->snp_srt.n > 0) { + radix_sort_bc64(hap->snp_srt.a, hap->snp_srt.a + hap->snp_srt.n);///sort by how many snps in one overlap + 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; + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->occ_0 < 2 || s->occ_1 < 2) continue; + if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) o++; + } + if(o == 0) continue; + + 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){ + 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) { + ///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; + t->occ_0 -= hap->list[i].cov; + assert(t->occ_0 >= 1); + } + } + } + } + + 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; + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->occ_0 < 2 || s->occ_1 < 2) continue; + if(s->score == 1) o++; + } + ii = hap->list[l].overlapID; + ///for HiFi, do not flip trans to cis + // if(overlap_list->list[ii].is_match == 2 && o == 0) { + // overlap_list->list[ii].is_match = 1; + // } + if(overlap_list->list[ii].is_match == 1 && o > 0) { + overlap_list->list[ii].is_match = 2; + } + } + + + for (k = 1, l = 0; k <= hap->length; ++k) { ///reset snp_stat + if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { + 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) { + hap->snp_stat.a[hap->list[i].overlapSite].score = -1; + } + } + } + l = k; + } + } + } + + + + hap->snp_srt.n = 0; + for (k = 1, l = 0; k <= hap->length; ++k) { + if (k == hap->length || hap->list[k].overlapID != hap->list[l].overlapID) { + if(overlap_list->list[hap->list[l].overlapID].is_match == 2) { + l = k; + continue; + } + for (i = l, o = 0; i < k; i++) { + if(hap->list[i].type!=1) continue; + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->occ_0 < 2 || s->occ_1 < 2) continue; + if(s->occ_0 >= asm_opt.s_hap_cov && s->occ_1 >= asm_opt.infor_cov) continue; + if(s->score == 1) continue; + o++; + if(!km) kv_push(uint64_t, hap->snp_srt, hap->list[i].overlapSite); + else kv_push_km(km, uint64_t, hap->snp_srt, hap->list[i].overlapSite); + } + hap->snp_srt.n -= o; + ///there are at least two variants at one read + if(o>=(overlap_list->list[hap->list[l].overlapID].align_length*up)) { + radix_sort_bc64(hap->snp_srt.a + hap->snp_srt.n, hap->snp_srt.a + hap->snp_srt.n + o); + a = hap->snp_srt.a + hap->snp_srt.n; + for (i = z = 0; i < o; i++) { + if(i > 0) s = &(hap->snp_stat.a[a[i-1]]); + if(i + 1 < o) t = &(hap->snp_stat.a[a[i+1]]); + if(s && s->site + 32 > hap->snp_stat.a[a[i]].site) continue; + if(t && hap->snp_stat.a[a[i]].site + 32 > t->site) continue; + a[z] = a[i]; + z++; + } + if(z >= 2) hap->snp_srt.n += z; + } + l = k; + } + } + if (hap->snp_srt.n > 0) { + radix_sort_bc64(hap->snp_srt.a, hap->snp_srt.a + hap->snp_srt.n); + for (k = 1, l = 0; k <= hap->snp_srt.n; ++k) { + if(k == hap->snp_srt.n || hap->snp_srt.a[k] != hap->snp_srt.a[l]) { + if(k - l >= 2) hap->snp_stat.a[hap->snp_srt.a[l]].score = 1; + } + l = k; + } + } + + 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(hap->list[i].type==1 || hap->list[i].type==0) { + s = &(hap->snp_stat.a[hap->list[i].overlapSite]); + if(s->score == 1 && (!(s->occ_0 < 2 || s->occ_1 < 2))) { + overlap_list->list[ii].strong = 1; + if(hap->list[i].type==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; + } + } +} + /** void partition_overlaps(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy, haplotype_evdience_alloc* hap, @@ -7992,7 +8428,8 @@ 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_advance(hap, overlap_list, NULL); + generate_haplotypes_naive_HiFi(hap, overlap_list, 0.04, NULL); // 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); @@ -8362,6 +8799,18 @@ haplotype_evdience_alloc* hap, UC_Read* g_read, int64_t flanking, uint64_t yid) } **/ +void print_ovlp_occ_stat(overlap_region_alloc* overlap_list, uint32_t xlen, uint8_t is_match) +{ + uint64_t k; + for (k = 0; k < overlap_list->length; k++) { + if(overlap_list->list[k].is_match != is_match) continue; + fprintf(stderr, "[M::%s::xlen::%u] utg%.6dl(%c), is_match::%u, x::[%u, %u)\n", __func__, xlen, + (int32_t)overlap_list->list[k].y_id + 1, + "+-"[overlap_list->list[k].y_pos_strand], overlap_list->list[k].is_match, + overlap_list->list[k].x_pos_s, overlap_list->list[k].x_pos_e+1); + } +} + void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref, UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g, Graph* DAGCon, Cigar_record* current_cigar, @@ -8411,6 +8860,8 @@ void correct_ul_overlap(overlap_region_alloc* overlap_list, const ul_idx_t *uref // ovlp_occ(overlap_list, 0), ovlp_occ(overlap_list, 1), ovlp_occ(overlap_list, 2), ovlp_occ(overlap_list, 3)); ///after this function, overlap_list is sorted by x_pos_e; used for g_chain partition_ul_overlaps_advance(overlap_list, uref, g_read, overlap_read, dumy, hap, force_repeat, max_ov_diff_ec, w_inf.window_length, km); + // print_ovlp_occ_stat(overlap_list, g_read->length, 1); + // print_ovlp_occ_stat(overlap_list, g_read->length, 2); // fprintf(stderr, "[M::%s-end] occ[0]->%lu, occ[1]->%lu, occ[2]->%lu, occ[3]->%lu\n", __func__, // ovlp_occ(overlap_list, 0), ovlp_occ(overlap_list, 1), ovlp_occ(overlap_list, 2), ovlp_occ(overlap_list, 3)); // debug_phasing_status(overlap_list, uref->ug, 0, hap, g_read, 20, 1176); diff --git a/gfa_ut.cpp b/gfa_ut.cpp index cb7c5ca..914e6e3 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -6734,12 +6734,12 @@ void print_raw_uls_seq(ul_resolve_t *uidx, const char *nn) ma_ug_t *ug = uidx->init_ug; all_ul_t *aln = uidx->idx; uint64_t id; uc_block_t *a = NULL; int64_t k, a_n; for (id = 0; id < aln->n; id++) { - a = aln->a[id].bb.a; a_n = aln->a[id].bb.n; + a = aln->a[id].bb.a; a_n = aln->a[id].bb.n; k = 0; if(a_n == 0) continue; fprintf(fp,"%.*s\tid::%lu\t", (int32_t)aln->nid.a[id].n, aln->nid.a[id].a, id); - for (k = 0; k < a_n && ug_occ_w(a[k].ts, a[k].te, &(ug->u.a[a[k].hid])) == 0; k++); + // for (k = 0; k < a_n && ug_occ_w(a[k].ts, a[k].te, &(ug->u.a[a[k].hid])) == 0; k++); for (; k < a_n; k++) { - if(ug_occ_w(a[k].ts, a[k].te, &(ug->u.a[a[k].hid])) == 0) break; + // if(ug_occ_w(a[k].ts, a[k].te, &(ug->u.a[a[k].hid])) == 0) break; fprintf(fp, "utg%.6d%c(%c)\t", a[k].hid + 1, "lc"[ug->u.a[a[k].hid].circ], "+-"[a[k].rev]); } fprintf(fp,"\n"); @@ -11914,7 +11914,7 @@ double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_ ma_ug_t *init_ug = ul_realignment(uopt, sg, 0); // exit(1); filter_sg_by_ug(sg, init_ug, uopt); - // print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1); + // print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-0"); bub = gen_bubble_chain(sg, init_ug, uopt, &r_het); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-1"); diff --git a/inter.cpp b/inter.cpp index 7c3242c..f112bcf 100644 --- a/inter.cpp +++ b/inter.cpp @@ -3801,12 +3801,15 @@ uint64_t mode, All_reads *ridx, ma_ug_t *ug, int64_t debug_i, void *km) // if(!mode) { // fprintf(stderr, "[M::utg%.6d%c] qs->%u; qe->%u\n", li->tn+1, "lc"[uref->ug->u.a[li->tn].circ], li->qs, li->qe); // } + // fprintf(stderr, "[M::utg%.6d%c] qs->%u; qe->%u; sc->%ld; pre->%ld; el->%u;\n", + // li->tn+1, "lc"[uref->ug->u.a[li->tn].circ], li->qs, li->qe, pop_sc(track[i]), pop_pre(track[i]), li->el); } int64_t n_v, n_u, n_v0, le, lnv; radix_sort_gfa64(srt, srt+res->n); for (k = (int64_t)res->n-1, n_v = n_u = 0; k >= 0; --k) { n_v0 = n_v; i = (uint32_t)srt[k]; + // fprintf(stderr, "+[M::utg%.6d%c] n_v0->%ld;\n", res->a[i].tn+1, "lc"[uref->ug->u.a[res->a[i].tn].circ], n_v0); if(res->a[i].el) { ///chain must start from cis alignments for (le = -1; i >= 0 && (track[i]&((uint64_t)0x80000000)) == 0;) { if(res->a[i].el) { @@ -3819,6 +3822,7 @@ uint64_t mode, All_reads *ridx, ma_ug_t *ug, int64_t debug_i, void *km) i = pop_pre(track[i]); } } + // fprintf(stderr, "-[M::] n_v->%ld;\n", n_v); if(n_v0 == n_v) continue; if(le >= 0) { i = le; n_v = lnv; @@ -4548,6 +4552,7 @@ void l2g_chain(const ul_idx_t *uref, kv_ul_ov_t *lidx, vec_mg_lchain_t *res) int64_t l2g_res_chain(ma_ug_t *ug, ul_ov_t *a, uint64_t a_n, vec_mg_lchain_t *gchains, double diff_rate) { + // fprintf(stderr, "[M::%s::] a_n::%lu\n", __func__, a_n); if(a_n <= 0) return 0; uint64_t k, m; int64_t l, rs, re, qs, qe, dq, dr, dif, mm; a_n++; asg_t *g = ug->g; gchains->n = 0; kv_resize(mg_lchain_t, *gchains, a_n); gchains->n = a_n; @@ -4564,6 +4569,7 @@ int64_t l2g_res_chain(ma_ug_t *ug, ul_ov_t *a, uint64_t a_n, vec_mg_lchain_t *gc assert(gchains->a[k-1].dist_pre >= 0); l += g->seq[gchains->a[k-1].v>>1].len + gchains->a[k-1].dist_pre; } + // fprintf(stderr, "[M::%s::k->%lu] utg%.6dl(%c)\n", __func__, k, (int32_t)(gchains->a[k].v>>1)+1, "+-"[gchains->a[k].v&1]); } if(diff_rate < 0) return 1; @@ -4972,6 +4978,8 @@ void update_ul_vec_t_ug(const ul_idx_t *uref, ul_vec_t *rch, vec_mg_lchain_t *uc z->qs = a[m].qs; z->qe = a[m].qe; z->te = a[m].re; z->ts = a[m].rs; z->pidx = k + 1 + m; z->pdis = z->aidx = (uint32_t)-1; + // fprintf(stderr, "[M::%s::k->%ld] m->%ld, utg%.6dl(%c)\n", + // __func__, k, m, (int32_t)z->hid+1, "+-"[z->rev]); } } @@ -5017,6 +5025,9 @@ void update_ul_vec_t_ug(const ul_idx_t *uref, ul_vec_t *rch, vec_mg_lchain_t *uc } else if(l < ((int64_t)rch->rlen)*0.001) { rch->dd = 2; } + + // fprintf(stderr, "[M::%s::] rch->dd::%u, rch->bb.n::%u\n", + // __func__, rch->dd, (uint32_t)rch->bb.n); } void print_raw_chains(vec_mg_lchain_t *uc, int64_t ulid) @@ -5298,9 +5309,10 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // if(s->id+i!=41927 && s->id+i!=47072 && s->id+i!=67641 && s->id+i!=90305 && s->id+i!=698342 && s->id+i!=329421) { // return; // } - // if(s->id+i!=41927) return; + // if(s->id+i!=43) return; - // fprintf(stderr, "\n[M::%s] rid:%ld, len:%lu\n", __func__, s->id+i, s->len[i]); + // fprintf(stderr, "\n[M::%s] rid::%ld, len::%lu, name::%.*s\n", __func__, s->id+i, s->len[i], + // (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a); // if (memcmp(UL_INF.nid.a[s->id+i].a, "d0aab024-b3a7-40fb-83cc-22c3d6d951f8", UL_INF.nid.a[s->id+i].n-1)) return; // fprintf(stderr, "[M::%s::] ==> len: %lu\n", __func__, s->len[i]); ha_get_ul_candidates_interface(b->abl, i, s->seq[i], s->len[i], s->opt->w, s->opt->k, s->uu, &b->olist, &b->olist_hp, &b->clist, s->opt->bw_thres, @@ -10673,8 +10685,9 @@ ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_c // detect_outlier_len("ul_realignment"); clear_all_ul_t(&UL_INF); ///for debug interval - if(!load_all_ul_t(&UL_INF, gfa_name, &R_INF, ug)) { + if(!load_all_ul_t(&UL_INF, gfa_name, &R_INF, ug)/**1**/) { gen_UL_reovlps(&sl, ug, sg, gfa_name, cutoff); + // exit(1); write_all_ul_t(&UL_INF, gfa_name, ug); } else if(double_check_cache){ if(drenew_UL_reovlps(&sl, ug, sg, gfa_name, cutoff)) {