xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-09-29 23:28:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

revised HiFi phasing

Browse Source
This commit is contained in:
chhylp123
2022-08-23 10:36:35 -04:00
parent 54465dd9ac
commit 8b7ca5f0c3
4 changed files with 482 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.h
+1 -1
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.16.7-r418"
#define HA_VERSION "0.16.8-r420"
#define VERBOSE 0
Correct.cpp
+461 -10
View File
@@ -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);
gfa_ut.cpp
+4 -4
View File
@@ -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");
inter.cpp
+16 -3
View File
@@ -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)) {