From 5d506e82c91eca2cebb90bb88e62f87be03a2824 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Mon, 28 Dec 2020 03:35:17 -0500 Subject: [PATCH] phasing --- Overlaps.cpp | 89 +++++--- Overlaps.h | 1 + hic.cpp | 582 ++++++++++++++++++++++++++++++++++++++++++++++----- 3 files changed, 585 insertions(+), 87 deletions(-) diff --git a/Overlaps.cpp b/Overlaps.cpp index 796f77e..c93dae7 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -11970,6 +11970,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, int max_hang, int min_ovlp) asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, 0, 0, 0, 1, NULL); dip_thres = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE)*0.75; asm_opt.hom_global_coverage = tmp_cov; + for (i = 0; i < ug->g->n_seq; i++) { @@ -12026,52 +12027,48 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov { kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); - ma_ug_t *ug = NULL, *copy = NULL; + ma_ug_t *ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); + + + + hc_links link; init_hc_links(&link, ug->g->n_seq, R_INF.total_reads); - copy = copy_untig_graph(ug); + asg_t *copy_sg = copy_read_graph(sg); + ma_ug_t *copy_ug = copy_untig_graph(ug); asm_opt.purge_overlap_len = asm_opt.purge_overlap_len_hic; asm_opt.purge_simi_rate = asm_opt.purge_simi_rate_hic; - adjust_utg_by_primary(©, sg, TRIO_THRES, sources, reverse_sources, coverage_cut, + adjust_utg_by_primary(©_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, bubble_dist, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, max_hang, min_ovlp, &new_rtg_edges, &link); - ma_ug_destroy(copy); + ma_ug_destroy(copy_ug); + asg_destroy(copy_sg); + + + + + new_rtg_edges.a.n = 0; ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp); - - /** - fprintf(stderr, "Writing raw unitig GFA to disk... \n"); - char* gfa_name = (char*)malloc(strlen(output_file_name)+25); - sprintf(gfa_name, "%s.r_utg.gfa", output_file_name); - FILE* output_file = fopen(gfa_name, "w"); - ma_ug_print(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); - fclose(output_file); - sprintf(gfa_name, "%s.r_utg.noseq.gfa", output_file_name); - output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); - fclose(output_file); - if(asm_opt.bed_inconsist_rate != 0) - { - sprintf(gfa_name, "%s.r_utg.lowQ.bed", output_file_name); - output_file = fopen(gfa_name, "w"); - ma_ug_print_bed(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, - max_hang, min_ovlp, asm_opt.bed_inconsist_rate, "utg", output_file); - fclose(output_file); - } - - free(gfa_name); - **/ classify_untigs(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, &new_rtg_edges, max_hang, min_ovlp); hic_analysis(ug, sg, &link); destory_hc_links(&link); ma_ug_destroy(ug); kv_destroy(new_rtg_edges.a); + + output_unitig_graph(sg, coverage_cut, output_file_name, sources, ruIndex, max_hang, min_ovlp); + output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, + reverse_sources, bubble_dist, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, + 0.05, 0.9, max_hang, min_ovlp, 0); + output_trio_unitig_graph(sg, coverage_cut, output_file_name, MOTHER, sources, + reverse_sources, bubble_dist, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, + 0.05, 0.9, max_hang, min_ovlp, 0); } ma_ug_t* merge_utg(ma_ug_t **dest, ma_ug_t **src) @@ -13811,6 +13808,41 @@ ma_ug_t* copy_untig_graph(ma_ug_t *src) return ug; } +asg_t* copy_read_graph(asg_t *src) +{ + asg_t *dest = NULL; + dest = asg_init(); + + dest->r_seq = src->r_seq; + dest->m_seq = dest->n_seq = src->n_seq; + dest->seq = (asg_seq_t*)malloc(dest->n_seq * sizeof(asg_seq_t)); + memcpy(dest->seq, src->seq, sizeof(asg_seq_t)*dest->n_seq); + + dest->m_arc = dest->n_arc = src->n_arc; + dest->arc = (asg_arc_t*)malloc(dest->n_arc*sizeof(asg_arc_t)); + memcpy(dest->arc, src->arc, sizeof(asg_arc_t)*dest->n_arc); + + dest->is_srt = src->is_srt; + dest->is_symm = src->is_symm; + dest->idx = (uint64_t*)malloc(dest->n_seq*2*8); + memcpy(dest->idx, src->idx, dest->n_seq*2*8); + asg_cleanup(dest); + + if(src->seq_vis) + { + dest->seq_vis = (uint8_t*)malloc(dest->n_seq*2*sizeof(uint8_t)); + memcpy(dest->seq_vis, src->seq_vis, dest->n_seq*2*sizeof(uint8_t)); + } + + if(src->n_F_seq > 0 && src->F_seq) + { + dest->n_F_seq = src->n_F_seq; + dest->F_seq = (ma_utg_t*)malloc(dest->n_F_seq*sizeof(ma_utg_t)); + memcpy(dest->F_seq, src->F_seq, dest->n_F_seq*sizeof(ma_utg_t)); + } + return dest; +} + void clean_trio_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, @@ -28063,8 +28095,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); sprintf(buf, "%s.hic", output_file_name); output_hic_graph(sg, coverage_cut, buf, sources, reverse_sources, bubble_dist, - (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, - mini_overlap_length); + (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length); free(buf); } else diff --git a/Overlaps.h b/Overlaps.h index 8419147..5f25be4 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1101,6 +1101,7 @@ ma_ug_t* output_trio_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* outpu uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench); +asg_t* copy_read_graph(asg_t *src); #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/hic.cpp b/hic.cpp index 2f7b737..b248df7 100644 --- a/hic.cpp +++ b/hic.cpp @@ -61,6 +61,8 @@ typedef struct{ kvec_t(uint32_t) a; uint32_t h[2]; uint8_t full_bub; + int status[2]; + double weight[2]; }partition_warp; typedef struct{ @@ -73,10 +75,10 @@ typedef struct{ uint64_t n; uint8_t* lock; uint32_t* hap; - double* weight; uint32_t m[3]; uint32_t label, label_add; hc_links* link; + G_partition g_p; }H_partition; typedef struct { @@ -2713,7 +2715,16 @@ int load_hc_hits(kvec_pe_hit* hits, const char *fn) return 1; } -void print_hc_links(hc_links* link, int dir) +inline int get_phase_status(H_partition* hap, uint32_t uID) +{ + int d = -2; + if(hap->hap[uID] & hap->m[0]) d = 1; + if(hap->hap[uID] & hap->m[1]) d = -1; + if(hap->hap[uID] & hap->m[2]) d = 0; + return d; +} + +void print_hc_links(hc_links* link, int dir, H_partition* hap) { uint64_t i, k; if(dir == 0) @@ -2723,9 +2734,11 @@ void print_hc_links(hc_links* link, int dir) for (k = 0; k < link->a.a[i].e.n; k++) { if(link->a.a[i].e.a[k].del) continue; - fprintf(stderr, "s-utg%.6d(%c)\td-utg%.6d(%c)\t%lu\t%c\t%f\te\n", + fprintf(stderr, "s-utg%.6d(%c)\tCLU:%u:%d\td-utg%.6d(%c)\tCLU:%u:%d\t%lu\t%c\t%f\te\n", (int)(i+1), "01"[!!(link->a.a[i].e.a[k].dis&(uint64_t)2)], + hap->hap[i]>>3, get_phase_status(hap, i), (int)(link->a.a[i].e.a[k].uID+1), "01"[!!(link->a.a[i].e.a[k].dis&(uint64_t)1)], + hap->hap[link->a.a[i].e.a[k].uID]>>3, get_phase_status(hap, link->a.a[i].e.a[k].uID), link->a.a[i].e.a[k].dis == (uint64_t)-1? (uint64_t)-1 : link->a.a[i].e.a[k].dis>>3, "fb"[!!(link->a.a[i].e.a[k].dis&(uint64_t)4)], link->a.a[i].e.a[k].weight); } @@ -3750,11 +3763,18 @@ inline void get_phased_block(G_partition* x, bubble_type* bub, uint64_t id, uint32_t* beg, uint32_t* sink, uint32_t** h0, uint32_t* h0_n, uint32_t** h1, uint32_t* h1_n, uint32_t* phased, uint32_t* bub_id) { - if(bub && beg && sink) + if(bub && beg && sink && bub_id) { - (*bub_id) = bub->index[x->a[id].a.a[0]]; - (*beg) = bub->list.a[bub->num.a[(*bub_id)]]; - (*sink) = bub->list.a[bub->num.a[(*bub_id)] + 1]; + (*bub_id) = (*beg) = (*sink) = (uint32_t)-1; + if(x->a[id].a.n > 0) + { + (*bub_id) = bub->index[x->a[id].a.a[0]]; + if((*bub_id) < bub->num.n) + { + (*beg) = bub->list.a[bub->num.a[(*bub_id)]]; + (*sink) = bub->list.a[bub->num.a[(*bub_id)] + 1]; + } + } } (*h0) = x->a[id].a.a; @@ -3762,7 +3782,9 @@ uint32_t* phased, uint32_t* bub_id) (*h1) = x->a[id].a.a + x->a[id].h[0]; (*h1_n) = x->a[id].h[1]; - (*phased) = x->a[id].full_bub; + if(phased) (*phased) = x->a[id].full_bub; + if((*h0_n) == 0) (*h0) = NULL; + if((*h1_n) == 0) (*h1) = NULL; } double get_co_weight(uint32_t *query, uint32_t query_n, uint32_t *target, uint32_t target_n, min_cut_t* m) @@ -3972,7 +3994,7 @@ void debug_hc_links(ha_ug_index* idx, hc_links* link, sldat_t* sl, bubble_type* destory_hc_links(link); init_hc_links(link, sl->idx->ug->g->n_seq, R_INF.total_reads); collect_hc_links(sl->idx, &sl->hits, link, bub); - print_hc_links(link, 0); + ///print_hc_links(link, 0); } uint64_t get_hic_distance(pe_hit* hit, hc_links* link, const ha_ug_index* idx) @@ -4035,8 +4057,8 @@ inline double get_trans(const ha_ug_index* idx, uint64_t x) inline double get_trans_weight(const ha_ug_index* idx, uint64_t x) { - #define OFFSET_RATE 0.0000001 - #define OFFSET_RATE_THRES 16.118095551 + #define OFFSET_RATE 0.000000001 + #define OFFSET_RATE_THRES 20.7232658359 long double rate = get_trans(idx, x); if(rate < 0) rate = 0; rate += OFFSET_RATE; @@ -4627,7 +4649,7 @@ uint32_t get_available_com(H_partition* hap, bubble_type* bub, ma_ug_t *ug, uint { hc_links* link = hap->link; uint32_t beg, sink, n, *a, i, j, k, uID, max_bub_i, max_non_bub_i, max_i, is_ava; - double w, max_bub_w, max_non_bub_w; + double w, cur_w0 = 0, cur_w1 = 0, max_bub_w, max_non_bub_w; max_i = (uint32_t)-1; for (i = 0, max_bub_w = -1, max_bub_i = (uint32_t)-1; i < bub->num.n-1; i++) @@ -4638,12 +4660,25 @@ uint32_t get_available_com(H_partition* hap, bubble_type* bub, ma_ug_t *ug, uint { uID = a[j]>>1; if(check_self && is_hap_set(uID, *hap)) break; - for (k = 0; k < link->a.a[uID].e.n; k++) + if(check_others) { - if(link->a.a[uID].e.a[k].del) continue; - if(check_others && (!is_hap_set(link->a.a[uID].e.a[k].uID, *hap))) continue; - w += link->a.a[uID].e.a[k].weight; - is_ava = 1; + get_related_weight(uID, hap, &cur_w0, &cur_w1); + w += (MAX(cur_w0, cur_w1) - MIN(cur_w0, cur_w1)); + for (k = 0; k < link->a.a[uID].e.n; k++) + { + if(link->a.a[uID].e.a[k].del) continue; + if(!is_hap_set(link->a.a[uID].e.a[k].uID, *hap)) continue; + is_ava = 1; + } + } + else + { + for (k = 0; k < link->a.a[uID].e.n; k++) + { + if(link->a.a[uID].e.a[k].del) continue; + w += link->a.a[uID].e.a[k].weight; + is_ava = 1; + } } } if(j != n) continue; @@ -4662,13 +4697,27 @@ uint32_t get_available_com(H_partition* hap, bubble_type* bub, ma_ug_t *ug, uint { uID = i; if(check_self && is_hap_set(uID, *hap)) continue; - for (k = 0, w = 0, is_ava = 0; k < link->a.a[uID].e.n; k++) + if(check_others) { - if(link->a.a[uID].e.a[k].del) continue; - if(check_others && (!is_hap_set(link->a.a[uID].e.a[k].uID, *hap))) continue; - w += link->a.a[uID].e.a[k].weight; - is_ava = 1; + get_related_weight(uID, hap, &cur_w0, &cur_w1); + w += (MAX(cur_w0, cur_w1) - MIN(cur_w0, cur_w1)); + for (k = 0, w = 0, is_ava = 0; k < link->a.a[uID].e.n; k++) + { + if(link->a.a[uID].e.a[k].del) continue; + if(!is_hap_set(link->a.a[uID].e.a[k].uID, *hap)) continue; + is_ava = 1; + } } + else + { + for (k = 0, w = 0, is_ava = 0; k < link->a.a[uID].e.n; k++) + { + if(link->a.a[uID].e.a[k].del) continue; + w += link->a.a[uID].e.a[k].weight; + is_ava = 1; + } + } + if(is_ava == 0) continue; if(w > max_non_bub_w) @@ -4773,15 +4822,16 @@ void phase_com(H_partition* hap, ma_ug_t *ug, bub_p_t_warp* b, bubble_type* bub, { uint32_t beg = (uint32_t)-1, sink = (uint32_t)-1, n, *a; get_bubbles(bub, bid>>1, &beg, &sink, &a, &n, NULL); - fprintf(stderr, "bubble-%uth, beg: %u, sink: %u, phasing ID: %u\n", bid>>1, beg>>1, sink>>1, hap->label>>3); + ///fprintf(stderr, "+bubble-%uth, beg: %u, sink: %u, phasing ID: %u\n", bid>>1, beg>>1, sink>>1, hap->label>>3); get_phase_path(ug, beg, sink, b, hap); get_phase_path(ug, beg, sink, b, hap); + ///fprintf(stderr, "-bubble-%uth, beg: %u, sink: %u, phasing ID: %u\n", bid>>1, beg>>1, sink>>1, hap->label>>3); } else { double cur_w0, cur_w1; get_related_weight(bid>>1, hap, &cur_w0, &cur_w1); - fprintf(stderr, "utg-%uth, phasing ID: %u\n", bid>>1, hap->label>>3); + ///fprintf(stderr, "utg-%uth, phasing ID: %u\n", bid>>1, hap->label>>3); if(cur_w0 >= cur_w1) { hap->hap[bid>>1] |= (hap->label | hap->m[0]); @@ -4793,13 +4843,169 @@ void phase_com(H_partition* hap, ma_ug_t *ug, bub_p_t_warp* b, bubble_type* bub, } } -inline uint32_t get_phase_status(H_partition* hap, uint32_t uID) + +double get_cluster_weight(H_partition* hap, hc_links* link, uint32_t *h, uint32_t h_n) { - int d = -2; - if(hap->hap[uID] & hap->m[0]) d = 1; - if(hap->hap[uID] & hap->m[1]) d = -1; - if(hap->hap[uID] & hap->m[2]) d = 0; - return d; + int o_d = 0; + double weight = 0; + uint32_t j, k, m; + for (j = 0, weight = 0; j < h_n; j++) + { + for (k = 0; k < link->a.a[h[j]].e.n; k++) + { + if(link->a.a[h[j]].e.a[k].del) continue; + for (m = 0; m < h_n; m++) + { + if(h[m] == link->a.a[h[j]].e.a[k].uID) break; + } + if(m < h_n) continue; + + o_d = get_phase_status(hap, link->a.a[h[j]].e.a[k].uID); + ///if(o_d < -1) fprintf(stderr, "ERROR\n"); + weight += (o_d*link->a.a[h[j]].e.a[k].weight); + } + } + + return weight; +} + +void update_partition_flag(H_partition* hap, hc_links* link, uint32_t id) +{ + uint32_t k, *h0, h0_n, *h1, h1_n, uID, flag = 0; + int status; + get_phased_block(&(hap->g_p), NULL, id, NULL, NULL, &h0, &h0_n, &h1, &h1_n, NULL, NULL); + + status = hap->g_p.a[id].status[0]; + if(status == 1) flag = hap->m[0]; + else if(status == -1) flag = hap->m[1]; + else if(status == 0) flag = hap->m[2]; + else if(status == -2) flag = 0; + for (k = 0; k < h0_n; k++) + { + uID = h0[k]; + hap->hap[uID] >>= 3; + hap->hap[uID] <<= 3; + hap->hap[uID] |= flag; + } + + status = hap->g_p.a[id].status[1]; + if(status == 1) flag = hap->m[0]; + else if(status == -1) flag = hap->m[1]; + else if(status == 0) flag = hap->m[2]; + else if(status == -2) flag = 0; + for (k = 0; k < h1_n; k++) + { + uID = h1[k]; + hap->hap[uID] >>= 3; + hap->hap[uID] <<= 3; + hap->hap[uID] |= flag; + } +} + +void print_contig_partition(H_partition* hap, const char* debug) +{ + uint32_t i; + int status; + for (i = 0; i < hap->n; i++) + { + status = get_phase_status(hap, i); + fprintf(stderr, "%s\tutg%.6d\tP:%u\tHG:A:%d\n", debug, (int)(i+1), hap->hap[i]>>3, status); + } +} + +void adjust_contig_partition(H_partition* hap, hc_links* link) +{ + uint32_t i, k, *h0, h0_n, *h1, h1_n; + uint32_t h0_status[4], h1_status[4], h0_status_max; + int h0_h, h1_h; + for (i = 0; i < hap->g_p.n; i++) + { + + get_phased_block(&(hap->g_p), NULL, i, NULL, NULL, &h0, &h0_n, &h1, &h1_n, NULL, NULL); + hap->g_p.a[i].status[0] = hap->g_p.a[i].status[1] = -2; + hap->g_p.a[i].weight[0] = hap->g_p.a[i].weight[1] = 0; + + + h0_status[0] = h0_status[1] = h0_status[2] = h0_status[3] = 0; + for (k = 0; k < h0_n; k++) + { + ///fprintf(stderr, "(0) utg%.6ul\n", h0[k] + 1); + ///if((hap->g_p.index[h0[k]]>>1) != i) fprintf(stderr, "ERROR\n"); + h0_status[get_phase_status(hap, h0[k])+2]++; + } + + hap->g_p.a[i].weight[0] = get_cluster_weight(hap, link, h0, h0_n); + if(h1_n == 0) + { + if(h0_status[0] == h0_n) + { + hap->g_p.a[i].status[0] = -2; + } + else + { + if(h0_status[1] > 0 || h0_status[3] > 0) + { + h0_status[0] = h0_status[2] = 0; + + h0_h = -2; + h0_status_max = 0; + for (k = 0; k < 4; k++) + { + if(h0_status[k] > h0_status_max) h0_status_max = h0_status[k], h0_h = (int)(k) - 2; + } + hap->g_p.a[i].status[0] = h0_h; + } + else if(h0_status[2] > 0) + { + hap->g_p.a[i].status[0] = 0; + } + else + { + hap->g_p.a[i].status[0] = -2; + } + } + } + else + { + hap->g_p.a[i].weight[1] = get_cluster_weight( hap, link, h1, h1_n); + h1_status[0] = h1_status[1] = h1_status[2] = h1_status[3] = 0; + for (k = 0; k < h1_n; k++) + { + ///fprintf(stderr, "(1) utg%.6ul\n", h1[k] + 1); + ///if((hap->g_p.index[h1[k]]>>1) != i) fprintf(stderr, "ERROR\n"); + h1_status[get_phase_status(hap, h1[k])+2]++; + } + + + h0_h = h1_h = 0; + for (k = 0; k < 4; k++) + { + if(h0_status[k] == h0_n) h0_h = (int)(k) - 2; + if(h1_status[k] == h1_n) h1_h = (int)(k) - 2; + } + + if(h0_h * h1_h == -1) + { + hap->g_p.a[i].status[0] = h0_h; + hap->g_p.a[i].status[1] = h1_h; + } + else + { + if(hap->g_p.a[i].weight[0] >= hap->g_p.a[i].weight[1]) + { + hap->g_p.a[i].status[0] = 1; + hap->g_p.a[i].status[1] = -1; + } + else + { + hap->g_p.a[i].status[0] = -1; + hap->g_p.a[i].status[1] = 1; + } + } + } + + update_partition_flag(hap, link, i); + } } uint32_t init_contig_partition(H_partition* hap, ha_ug_index* idx, bubble_type* bub) @@ -4809,8 +5015,7 @@ uint32_t init_contig_partition(H_partition* hap, ha_ug_index* idx, bubble_type* bub_p_t_warp b; memset(&b, 0, sizeof(bub_p_t_warp)); CALLOC(b.a, ug->g->n_seq*2); - uint32_t i, k, nv = ug->g->n_seq * 2; - int s_d, o_d; + uint32_t i, k, k_n, nv = ug->g->n_seq * 2; for (i = 0; i < nv; i++) { b.a[i].w[0] = b.a[i].w[1] = b.a[i].nh = 0; @@ -4823,11 +5028,11 @@ uint32_t init_contig_partition(H_partition* hap, ha_ug_index* idx, bubble_type* memset(hap->hap, 0, hap->n*sizeof(uint32_t)); MALLOC(hap->lock, hap->n); memset(hap->lock, 0, hap->n); - MALLOC(hap->weight, hap->n); hap->m[0] = 1; hap->m[1] = 2; hap->m[2] = 4; hap->link = link; hap->label = 0; hap->label_add = 8; + init_G_partition(&(hap->g_p), hap->n); uint32_t max_i = get_available_com(hap, bub, ug, 0, 0); @@ -4841,52 +5046,308 @@ uint32_t init_contig_partition(H_partition* hap, ha_ug_index* idx, bubble_type* phase_com(hap, ug, &b, bub, max_i); } - free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); - - s_d = 0; for (i = 0; i < hap->n; i++) { - if(bub->index[i] > bub->num.n) continue; - s_d = get_phase_status(hap, i); - ///if(s_d == -2 && link->a.a[i].e.n > 0) fprintf(stderr, "ERROR\n"); - if(s_d == -2) continue; - hap->weight[i] = 0; - - for (k = 0; k < link->a.a[i].e.n; k++) + if((hap->hap[i]&hap->m[0])&&(hap->hap[i]&hap->m[1])) { - if(link->a.a[i].e.a[k].del) continue; - o_d = get_phase_status(hap, link->a.a[i].e.a[k].uID); - ///if(o_d == -2) fprintf(stderr, "ERROR\n"); - hap->weight[i] += (s_d*o_d*link->a.a[i].e.a[k].weight); + hap->hap[i] >>= 3; + hap->hap[i] <<= 3; + hap->hap[i] |= hap->m[2]; } } - + + free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); + + partition_warp* res = NULL; + hc_edge *a = NULL; + uint32_t a_n, v, u, uv = (uint32_t)-1, k_nv, k_nu; + for (i = 0; i < hap->n; i++) + { + v = i; + a = link->a.a[v].f.a; + a_n = link->a.a[v].f.n; + for (k = k_n = 0; k < a_n; k++) + { + if(a[k].del) continue; + if(a[k].dis != 0) break; + u = a[k].uID; + k_n++; + } + if(k_n != 1) + { + u = (uint32_t)-1; + goto push_uv; + } + + a = link->a.a[u].f.a; + a_n = link->a.a[u].f.n; + for (k = k_n = 0; k < a_n; k++) + { + if(a[k].del) continue; + if(a[k].dis != 0) break; + uv = a[k].uID; + k_n++; + } + if(k_n != 1 || uv != v) + { + u = (uint32_t)-1; + goto push_uv; + } + + push_uv: + k_nv = 0;k_nu = 0; + + a = link->a.a[v].e.a; + a_n = link->a.a[v].e.n; + for (k = 0; k < a_n; k++) + { + if(a[k].del) continue; + k_nv++; + } + + if(u != (uint32_t)-1) + { + a = link->a.a[u].e.a; + a_n = link->a.a[u].e.n; + for (k = 0; k < a_n; k++) + { + if(a[k].del) continue; + k_nu++; + } + } + if(k_nv == 0) continue; + if(k_nv > 0 && k_nu > 0 && v > u) continue; + + kv_pushp(partition_warp, hap->g_p, &res); + kv_init(res->a); + res->full_bub = 0; + res->h[0] = 1; res->h[1] = 0; + kv_push(uint32_t, res->a, v); + if(u != (uint32_t)-1) + { + res->h[1] = 1; + kv_push(uint32_t, res->a, u); + } + + for (k = 0; k < res->h[0]; k++) + { + if(hap->g_p.index[res->a.a[k]] != (uint32_t)-1) fprintf(stderr, "ERROR\n"); + hap->g_p.index[res->a.a[k]] = hap->g_p.n-1; + hap->g_p.index[res->a.a[k]] = hap->g_p.index[res->a.a[k]] << 1; + } + + for (; k < res->a.n; k++) + { + if(hap->g_p.index[res->a.a[k]] != (uint32_t)-1) fprintf(stderr, "ERROR\n"); + hap->g_p.index[res->a.a[k]] = hap->g_p.n-1; + hap->g_p.index[res->a.a[k]] = (hap->g_p.index[res->a.a[k]] << 1) + 1; + } + } + + ///print_contig_partition(hap, "first"); + + adjust_contig_partition(hap, link); + + ///print_contig_partition(hap, "second"); + return 1; } + uint32_t get_max_unitig(H_partition* hap, hc_links* link, ma_ug_t *ug, bubble_type* bub) { - uint32_t i; - for (i = 0; i < hap->n; i++) + double min, weight; + uint32_t i, min_i; + + for (i = 0, min = 1, min_i = (uint32_t)-1; i < hap->g_p.n; i++) { + if(hap->lock[i]) continue; + weight = 0; + if(hap->g_p.a[i].h[0] > 0 && (hap->g_p.a[i].status[0] == 1 || hap->g_p.a[i].status[0] == -1)) + { + weight += (hap->g_p.a[i].weight[0] * hap->g_p.a[i].status[0]); + } + if(hap->g_p.a[i].h[1] > 0 && (hap->g_p.a[i].status[1] == 1 || hap->g_p.a[i].status[1] == -1)) + { + weight += (hap->g_p.a[i].weight[1] * hap->g_p.a[i].status[1]); + } + + if(weight >= 0) continue; + if(weight < min) + { + min = weight; + min_i = i; + } } - return 0; + return min_i; } -uint32_t phasing_improevment(H_partition* hap, ha_ug_index* idx, bubble_type* bub) -{ - ; +void flip_unitig(H_partition* hap, hc_links* link, ma_ug_t *ug, bubble_type* bub, uint32_t id) +{ + if(hap->g_p.a[id].h[0] > 0 && hap->g_p.a[id].status[0] != 1 && hap->g_p.a[id].status[0] != -1) return; + if(hap->g_p.a[id].h[1] > 0 && hap->g_p.a[id].status[1] != 1 && hap->g_p.a[id].status[1] != -1) return; + uint32_t k, j, m, *h0, h0_n, *h1, h1_n, uID, *h = NULL, h_n; + int status; + double weight; + get_phased_block(&(hap->g_p), NULL, id, NULL, NULL, &h0, &h0_n, &h1, &h1_n, NULL, NULL); + if(h0_n > 0) + { + status = hap->g_p.a[id].status[0]; + h = h0; h_n = h0_n; + for (j = 0; j < h_n; j++) + { + for (k = 0; k < link->a.a[h[j]].e.n; k++) + { + if(link->a.a[h[j]].e.a[k].del) continue; + for (m = 0; m < h_n; m++) + { + if(h[m] == link->a.a[h[j]].e.a[k].uID) break; + } + if(m < h_n) continue; - return 0; + uID = link->a.a[h[j]].e.a[k].uID; + weight = link->a.a[h[j]].e.a[k].weight; + hap->g_p.a[hap->g_p.index[uID]>>1].weight[hap->g_p.index[uID]&1] -= (2*status*weight); + } + } + hap->g_p.a[id].status[0] *= -1; + } + + if(h1_n > 0) + { + status = hap->g_p.a[id].status[1]; + h = h1; h_n = h1_n; + for (j = 0; j < h_n; j++) + { + for (k = 0; k < link->a.a[h[j]].e.n; k++) + { + if(link->a.a[h[j]].e.a[k].del) continue; + for (m = 0; m < h_n; m++) + { + if(h[m] == link->a.a[h[j]].e.a[k].uID) break; + } + if(m < h_n) continue; + + uID = link->a.a[h[j]].e.a[k].uID; + weight = link->a.a[h[j]].e.a[k].weight; + hap->g_p.a[hap->g_p.index[uID]>>1].weight[hap->g_p.index[uID]&1] -= (2*status*weight); + } + } + hap->g_p.a[id].status[1] *= -1; + } +} + + +uint32_t phasing_improvement(H_partition* hap, ha_ug_index* idx, bubble_type* bub) +{ + uint32_t i, occ = 0; + memset(hap->lock, 0, sizeof(uint8_t)*hap->g_p.n); + while (1) + { + i = get_max_unitig(hap, idx->link, idx->ug, bub); + if(i == (uint32_t)-1) break; + hap->lock[i] = 1; + flip_unitig(hap, idx->link, idx->ug, bub, i); + occ++; + } + + for (i = 0; i < hap->g_p.n; i++) + { + update_partition_flag(hap, idx->link, i); + } + + return !!occ; } void destory_contig_partition(H_partition* hap) { free(hap->lock); free(hap->hap); - free(hap->weight); + destory_G_partition(&(hap->g_p)); +} + +void label_unitigs(H_partition* hap, ma_ug_t* ug) +{ + memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads * sizeof(uint8_t)); + uint32_t i, k, j, *h0, h0_n, *h1, h1_n, uID, *h = NULL, h_n, flag = AMBIGU; + int status; + ma_utg_t *u = NULL; + + for (i = 0; i < hap->g_p.n; i++) + { + if(hap->g_p.a[i].h[0] > 0 && hap->g_p.a[i].status[0] != 1 && hap->g_p.a[i].status[0] != -1) continue; + if(hap->g_p.a[i].h[1] > 0 && hap->g_p.a[i].status[1] != 1 && hap->g_p.a[i].status[1] != -1) continue; + get_phased_block(&(hap->g_p), NULL, i, NULL, NULL, &h0, &h0_n, &h1, &h1_n, NULL, NULL); + + status = hap->g_p.a[i].status[0]; + h = h0; h_n = h0_n; + if(status == 1) + { + flag = FATHER; + } + else if (status == -1) + { + flag = MOTHER; + } + for (j = 0; j < h_n; j++) + { + uID = h[j]; + u = &ug->u.a[uID]; + if(u->m == 0) continue; + for (k = 0; k < u->n; k++) + { + R_INF.trio_flag[u->a[k]>>33] = flag; + } + } + + + + status = hap->g_p.a[i].status[1]; + h = h1; h_n = h1_n; + if(status == 1) + { + flag = FATHER; + } + else if (status == -1) + { + flag = MOTHER; + } + for (j = 0; j < h_n; j++) + { + uID = h[j]; + u = &ug->u.a[uID]; + if(u->m == 0) continue; + for (k = 0; k < u->n; k++) + { + R_INF.trio_flag[u->a[k]>>33] = flag; + } + } + } + + + uint64_t occ = 0; + for (i = 0; i < ug->u.n; i++) + { + occ += ug->u.a[i].n; + } + + fprintf(stderr, "# reads: %lu\n", occ); + + for (i = occ = 0; i < R_INF.total_reads; i++) + { + if(R_INF.trio_flag[i] == FATHER) occ++; + } + + fprintf(stderr, "# Father reads: %lu\n", occ); + + for (i = occ = 0; i < R_INF.total_reads; i++) + { + if(R_INF.trio_flag[i] == MOTHER) occ++; + } + + fprintf(stderr, "# Mother reads: %lu\n", occ); } int hic_short_align(const char *fn1, const char *fn2, ha_ug_index* idx) @@ -4934,11 +5395,16 @@ int hic_short_align(const char *fn1, const char *fn2, ha_ug_index* idx) ///debug_hc_links(idx, idx->link, &sl, &bub, fn1); init_hic_p((ha_ug_index*)sl.idx, &sl.hits, idx->link, &bub); - ///print_hc_links(idx->link, 0); + H_partition hap; init_contig_partition(&hap, idx, &bub); - destory_contig_partition(&hap); + phasing_improvement(&hap, idx, &bub); + label_unitigs(&hap, idx->ug); + + ///print_hc_links(idx->link, 0, &hap); + ///print_contig_partition(&hap, "final"); + destory_contig_partition(&hap); return 1;