From be8e1e3117e1f810e959d309e38c26e465955379 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Mon, 17 May 2021 00:12:41 -0400 Subject: [PATCH] sc debug --- Overlaps.cpp | 337 ++++++++++++++++++++++++++++++--------- Overlaps.h | 10 ++ hic.cpp | 2 +- horder.cpp | 436 ++++++++++++++++++++++++++++++++++++++++++++++----- horder.h | 3 +- 5 files changed, 670 insertions(+), 118 deletions(-) diff --git a/Overlaps.cpp b/Overlaps.cpp index 549b746..8b9aa7a 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -9103,6 +9103,183 @@ kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, kvec_asg_arc_t_warp *E, u } +void polish_unitig_scaffold(uint32_t uid, ma_utg_t* collection, asg_t* read_g, All_reads *RNF, +ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, +UC_Read* r_read, UC_Read* q_read, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* newE) +{ + uint32_t i, k, m, len, p_idx; + ma_utg_t qu; + for (i = p_idx = m = len = 0; i < collection->n; i++) + { + if(collection->a[i] == (uint64_t)-1) + { + qu.a = collection->a + p_idx; + qu.circ = 0; + qu.m = qu.n = i - p_idx; + for (k = qu.len = 0; k < qu.n; k++) + { + qu.len += (uint32_t)qu.a[k]; + } + + polish_unitig(&qu, read_g, sources, coverage_cut, edge, max_hang, min_ovlp, newE); + polish_unitig_advance(&qu, read_g, &R_INF, sources, coverage_cut, edge, r_read, q_read, max_hang, min_ovlp, newE); + + for (k = 0; k < qu.n; k++) + { + collection->a[m] = qu.a[k]; + m++; + } + len += qu.len; + + collection->a[m] = (uint64_t)-1; + m++; + len += GAP_LEN; + p_idx = i + 1; + } + } + + if(i - p_idx > 0) + { + qu.a = collection->a + p_idx; + qu.circ = collection->circ; + qu.m = qu.n = i - p_idx; + for (k = qu.len = 0; k < qu.n; k++) + { + qu.len += (uint32_t)qu.a[k]; + } + + polish_unitig(&qu, read_g, sources, coverage_cut, edge, max_hang, min_ovlp, newE); + polish_unitig_advance(&qu, read_g, &R_INF, sources, coverage_cut, edge, r_read, q_read, max_hang, min_ovlp, newE); + + for (k = 0; k < qu.n; k++) + { + collection->a[m] = qu.a[k]; + m++; + } + len += qu.len; + } + collection->n = m; + collection->len = len; +} + +int ma_ug_seq_scaffold(ma_ug_t *g, asg_t *read_g, ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, +kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, kvec_asg_arc_t_warp *E, uint32_t is_polish) +{ + UC_Read g_read; + init_UC_Read(&g_read); + UC_Read tmp; + init_UC_Read(&tmp); + ///utg_intv_t *tmp; + uint32_t i, j, k; + uint32_t rId, /**uId,**/ori, start, eLen, readLen; + char* readS = NULL; + if(!(g->g)) + { + g->g = asg_init(); + for (i = 0; i < g->u.n; ++i) + { + asg_seq_set(g->g, i, g->u.a[i].len, 0); + } + asg_cleanup(g->g); + g->g->r_seq = g->g->n_seq; + } + + ///why we need n_read here? it is just beacuse one read can only be included in one untig + ///but it is not true + ///tmp = (utg_intv_t*)calloc(n_read, sizeof(utg_intv_t)); + ///number of unitigs + for (i = 0; i < g->u.n; ++i) { + ma_utg_t *u = &g->u.a[i]; + if(u->m == 0) continue; + if(is_polish) polish_unitig_scaffold(i, u, read_g, &R_INF, sources, coverage_cut, edge, &g_read, &tmp, max_hang, min_ovlp, E); + g->g->seq[i].len = u->len; + + uint32_t l = 0; + u->s = (char*)calloc(1, u->len + 1); + memset(u->s, 'N', u->len); + for (j = 0; j < u->n; ++j) { + // fprintf(stderr, "j=%u, u->a[j]: %lu\n", j, u->a[j]); + if(u->a[j] == (uint64_t)-1) + { + l += GAP_LEN; + continue; + } + rId = u->a[j]>>33; + ///uId = i; + ori = u->a[j]>>32&1; + start = l; + eLen = (uint32_t)u->a[j]; + l += eLen; + + if(eLen == 0) continue; + if(rId < read_g->r_seq) + { + recover_UC_Read(&g_read, &R_INF, rId); + } + else + { + recover_fake_read(&g_read, &tmp, &(read_g->F_seq[rId-read_g->r_seq]), + &R_INF, coverage_cut); + } + + readS = g_read.seq + coverage_cut[rId].s; + readLen = coverage_cut[rId].e - coverage_cut[rId].s; + + if (!ori) // forward strand + { + for (k = 0; k < eLen; k++) + { + u->s[start + k] = readS[k]; + } + } + else + { + for (k = 0; k < eLen; k++) + { + uint8_t c = (uint8_t)readS[readLen - 1 - k]; + u->s[start + k] = c >= 128? 'N' : comp_tab[c]; + } + } + } + } + + destory_UC_Read(&g_read); + destory_UC_Read(&tmp); + + + uint32_t n_vtx = g->g->n_seq * 2, v, nv; + uint32_t vLen = 0; + asg_arc_t* av = NULL; + for (v = 0; v < n_vtx; ++v) + { + if (g->g->seq[v>>1].del) continue; + av = asg_arc_a(g->g, v); + nv = asg_arc_n(g->g, v); + for (i = 0; i < nv; i++) + { + if(av[i].del) continue; + vLen = g->g->seq[(av[i].ul>>33)].len - av[i].ol; + av[i].ul = (av[i].ul>>32)<<32; + av[i].ul = av[i].ul | vLen; + } + } + + if(E && E->a.n > 0) + { + asg_arc_t* p = NULL; + for (k = 0; k < E->a.n; k++) + { + p = asg_arc_pushp(read_g); + *p = E->a.a[k]; + } + + free(read_g->idx); + read_g->idx = 0; + read_g->is_srt = 0; + asg_cleanup(read_g); + } + return 0; +} uint32_t get_ug_coverage(ma_utg_t* u, asg_t* read_g, const ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag) { @@ -9113,12 +9290,14 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag) for (k = 0; k < u->n; k++) { + if(u->a[k] == (uint64_t)-1) continue; rId = u->a[k]>>33; r_flag[rId] = 1; } for (k = 0; k < u->n; k++) { + if(u->a[k] == (uint64_t)-1) continue; rId = u->a[k]>>33; R_bases += (coverage_cut[rId].e - coverage_cut[rId].s); for (j = 0; j < (uint64_t)(sources[rId].length); j++) @@ -9141,6 +9320,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag) for (k = 0; k < u->n; k++) { + if(u->a[k] == (uint64_t)-1) continue; rId = u->a[k]>>33; r_flag[rId] = 0; } @@ -9259,22 +9439,28 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FIL // if (print_seq) fprintf(fp, "S\t%s\t%s\tLN:i:%d\n", name, p->s? p->s : "*", p->len); // else fprintf(fp, "S\t%s\t*\tLN:i:%d\n", name, p->len); - for (j = l = 0; j < p->n; l += (uint32_t)p->a[j++]) { - uint32_t x = p->a[j]>>33; - if(xtotal_reads) + for (j = l = 0; j < p->n; j++) { + if(p->a[j] != (uint64_t)-1) { - fprintf(fp, "A\t%s\t%d\t%c\t%.*s\t%d\t%d\tid:i:%d\tHG:A:%c\n", name, l, "+-"[p->a[j]>>32&1], - (int)Get_NAME_LENGTH((*RNF), x), Get_NAME((*RNF), x), - coverage_cut[x].s, coverage_cut[x].e, x, - "apmaaa"[((RNF->trio_flag[x]!=FATHER && RNF->trio_flag[x]!=MOTHER)?AMBIGU:RNF->trio_flag[x])]); + uint32_t x = p->a[j]>>33; + if(xtotal_reads) + { + fprintf(fp, "A\t%s\t%d\t%c\t%.*s\t%d\t%d\tid:i:%d\tHG:A:%c\n", name, l, "+-"[p->a[j]>>32&1], + (int)Get_NAME_LENGTH((*RNF), x), Get_NAME((*RNF), x), + coverage_cut[x].s, coverage_cut[x].e, x, + "apmaaa"[((RNF->trio_flag[x]!=FATHER && RNF->trio_flag[x]!=MOTHER)?AMBIGU:RNF->trio_flag[x])]); + } + else + { + fprintf(fp, "A\t%s\t%d\t%c\t%s\t%d\t%d\tid:i:%d\tHG:A:%c\n", name, l, "+-"[p->a[j]>>32&1], + "FAKE", coverage_cut[x].s, coverage_cut[x].e, x, '*'); + } } else { - fprintf(fp, "A\t%s\t%d\t%c\t%s\t%d\t%d\tid:i:%d\tHG:A:%c\n", name, l, "+-"[p->a[j]>>32&1], - "FAKE", coverage_cut[x].s, coverage_cut[x].e, x, '*'); + fprintf(fp, "A\t%s\t%d\t*\t*\t*\t*\tid:i:*\tHG:A:*\n", name, l); } - - + l += (uint32_t)p->a[j]; } } // for (i = 0; i < ug->g->n_arc; ++i) { // the Link lines in GFA @@ -9283,49 +9469,52 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int print_seq, const char* prefix, FIL // prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], // prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], ug->g->arc[i].ol, asg_arc_len(ug->g->arc[i])); // } - asg_arc_t* au = NULL; - uint32_t nu, u, v; - for (i = 0; i < ug->u.n; ++i) { - if(ug->u.a[i].m == 0) continue; - if(ug->u.a[i].circ) - { - fprintf(fp, "L\t%s%.6dc\t+\t%s%.6dc\t+\t%dM\tL1:i:%d\n", - prefix, i+1, prefix, i+1, 0, ug->u.a[i].len); - fprintf(fp, "L\t%s%.6dc\t-\t%s%.6dc\t-\t%dM\tL1:i:%d\n", - prefix, i+1, prefix, i+1, 0, ug->u.a[i].len); - } - u = i<<1; - au = asg_arc_a(ug->g, u); - nu = asg_arc_n(ug->g, u); - for (j = 0; j < nu; j++) - { - if(au[j].del) continue; - v = au[j].v; - fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\n", - prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], - prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j])); - } + if(ug->g) + { + asg_arc_t* au = NULL; + uint32_t nu, u, v; + for (i = 0; i < ug->u.n; ++i) { + if(ug->u.a[i].m == 0) continue; + if(ug->u.a[i].circ) + { + fprintf(fp, "L\t%s%.6dc\t+\t%s%.6dc\t+\t%dM\tL1:i:%d\n", + prefix, i+1, prefix, i+1, 0, ug->u.a[i].len); + fprintf(fp, "L\t%s%.6dc\t-\t%s%.6dc\t-\t%dM\tL1:i:%d\n", + prefix, i+1, prefix, i+1, 0, ug->u.a[i].len); + } + u = i<<1; + au = asg_arc_a(ug->g, u); + nu = asg_arc_n(ug->g, u); + for (j = 0; j < nu; j++) + { + if(au[j].del) continue; + v = au[j].v; + fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\n", + prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], + prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j])); + } - u = (i<<1) + 1; - au = asg_arc_a(ug->g, u); - nu = asg_arc_n(ug->g, u); - for (j = 0; j < nu; j++) - { - if(au[j].del) continue; - v = au[j].v; - fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\n", - prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], - prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j])); + u = (i<<1) + 1; + au = asg_arc_a(ug->g, u); + nu = asg_arc_n(ug->g, u); + for (j = 0; j < nu; j++) + { + if(au[j].del) continue; + v = au[j].v; + fprintf(fp, "L\t%s%.6d%c\t%c\t%s%.6d%c\t%c\t%dM\tL1:i:%d\n", + prefix, (u>>1)+1, "lc"[ug->u.a[u>>1].circ], "+-"[u&1], + prefix, (v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], au[j].ol, asg_arc_len(au[j])); + } } } free(primary_flag); } -void ma_ug_print(const ma_ug_t *ug, All_reads *RNF, asg_t* read_g, const ma_sub_t *coverage_cut, +void ma_ug_print(const ma_ug_t *ug, asg_t* read_g, const ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, const char* prefix, FILE *fp) { - ma_ug_print2(ug, RNF, read_g, coverage_cut, sources, ruIndex, 1, prefix, fp); + ma_ug_print2(ug, &R_INF, read_g, coverage_cut, sources, ruIndex, 1, prefix, fp); } int asg_cut_internal(asg_t *g, int max_ext) @@ -9543,10 +9732,10 @@ void ma_sg_print(const asg_t *g, const All_reads *RNF, const ma_sub_t *sub, FILE } -void ma_ug_print_simple(const ma_ug_t *ug, All_reads *RNF, asg_t* read_g, const ma_sub_t *coverage_cut, +void ma_ug_print_simple(const ma_ug_t *ug, asg_t* read_g, const ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, const char* prefix, FILE *fp) { - ma_ug_print2(ug, RNF, read_g, coverage_cut, sources, ruIndex, 0, prefix, fp); + ma_ug_print2(ug, &R_INF, read_g, coverage_cut, sources, ruIndex, 0, prefix, fp); } void ma_ug_print_bed(const ma_ug_t *g, asg_t *read_g, All_reads *RNF, ma_sub_t *coverage_cut, @@ -12163,11 +12352,11 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) 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); + ma_ug_print(ug, 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); + ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file); fclose(output_file); if(asm_opt.bed_inconsist_rate != 0) { @@ -12504,12 +12693,12 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp, kvec_asg_a sprintf(gfa_name, "%s.p_ctg.gfa", output_file_name); fprintf(stderr, "Writing %s to disk... \n", gfa_name); FILE* output_file = fopen(gfa_name, "w"); - ma_ug_print(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "ptg", output_file); + ma_ug_print(ug, sg, coverage_cut, sources, ruIndex, "ptg", output_file); fclose(output_file); sprintf(gfa_name, "%s.p_ctg.noseq.gfa", output_file_name); output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "ptg", output_file); + ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "ptg", output_file); fclose(output_file); if(asm_opt.bed_inconsist_rate != 0) { @@ -12710,8 +12899,6 @@ trans_chain* load_hc_hits(const char *fn) void output_contig_graph_alternative(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp); -void reduce_hamming_error(asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, -int max_hang, int min_ovlp, long long gap_fuzz); void clean_u_trans_t_idx(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g); void output_hic_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, @@ -12734,6 +12921,7 @@ long long gap_fuzz, bub_label_t* b_mask_t) opt.min_ovlp = min_ovlp; opt.is_bench = 0; opt.b_mask_t = b_mask_t; + opt.gap_fuzz = gap_fuzz; kvec_asg_arc_t_warp new_rtg_edges, d_edges; @@ -12781,12 +12969,13 @@ long long gap_fuzz, bub_label_t* b_mask_t) if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_trans_chain(cov->t_ch, output_file_name); } - char* gfa_name = (char*)malloc(strlen(output_file_name)+25); - sprintf(gfa_name, "%s.pre.clean_d_utg.noseq.gfa", output_file_name); - FILE* output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); - fclose(output_file); - free(gfa_name); + // char* gfa_name = (char*)malloc(strlen(output_file_name)+50); + // FILE* output_file = NULL; + // sprintf(gfa_name, "%s.pre.clean_d_utg.noseq.gfa", output_file_name); + // FILE* output_file = fopen(gfa_name, "w"); + // ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file); + // fclose(output_file); + // free(gfa_name); hic_analysis(ug, sg, cov?cov->t_ch:t_ch, &opt); @@ -12794,12 +12983,12 @@ long long gap_fuzz, bub_label_t* b_mask_t) if(t_ch) destory_trans_chain(&t_ch); - gfa_name = (char*)malloc(strlen(output_file_name)+25); - sprintf(gfa_name, "%s.after.clean_d_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); - free(gfa_name); + // gfa_name = (char*)malloc(strlen(output_file_name)+25); + // sprintf(gfa_name, "%s.after.clean_d_utg.noseq.gfa", output_file_name); + // output_file = fopen(gfa_name, "w"); + // ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file); + // fclose(output_file); + // free(gfa_name); ma_ug_destroy(ug); kv_destroy(new_rtg_edges.a); @@ -13785,7 +13974,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) char* gfa_name = (char*)malloc(strlen(output_file_name)+25); sprintf(gfa_name, "%s.r_utg.noseq.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, read_g, coverage_cut, sources, ruIndex, "utg", output_file); + ma_ug_print_simple(ug, read_g, coverage_cut, sources, ruIndex, "utg", output_file); fclose(output_file); free(gfa_name); @@ -16197,12 +16386,12 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench, kv_destroy(new_rtg_edges.a); return ug; } - ma_ug_print(ug, &R_INF, sg, coverage_cut, sources, ruIndex, (flag==FATHER?"h1tg":"h2tg"), output_file); + ma_ug_print(ug, sg, coverage_cut, sources, ruIndex, (flag==FATHER?"h1tg":"h2tg"), output_file); fclose(output_file); sprintf(gfa_name, "%s.%s.p_ctg.noseq.gfa", output_file_name, (flag==FATHER?"hap1":"hap2")); output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, (flag==FATHER?"h1tg":"h2tg"), output_file); + ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, (flag==FATHER?"h1tg":"h2tg"), output_file); fclose(output_file); if(asm_opt.bed_inconsist_rate != 0) { @@ -24067,12 +24256,12 @@ R_to_U* ruIndex, int max_hang, int min_ovlp) char* gfa_name = (char*)malloc(strlen(output_file_name)+35); sprintf(gfa_name, "%s.p_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); + ma_ug_print(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file); fclose(output_file); sprintf(gfa_name, "%s.p_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); + ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file); fclose(output_file); if(asm_opt.bed_inconsist_rate != 0) { @@ -24119,12 +24308,12 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov char* gfa_name = (char*)malloc(strlen(output_file_name)+35); sprintf(gfa_name, "%s.p_ctg.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); - ma_ug_print(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "ptg", output_file); + ma_ug_print(ug, sg, coverage_cut, sources, ruIndex, "ptg", output_file); fclose(output_file); sprintf(gfa_name, "%s.p_ctg.noseq.gfa", output_file_name); output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "ptg", output_file); + ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "ptg", output_file); fclose(output_file); if(asm_opt.bed_inconsist_rate != 0) { @@ -24162,12 +24351,12 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) char* gfa_name = (char*)malloc(strlen(output_file_name)+35); sprintf(gfa_name, "%s.a_ctg.gfa", output_file_name); FILE* output_file = fopen(gfa_name, "w"); - ma_ug_print(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "atg", output_file); + ma_ug_print(ug, sg, coverage_cut, sources, ruIndex, "atg", output_file); fclose(output_file); sprintf(gfa_name, "%s.a_ctg.noseq.gfa", output_file_name); output_file = fopen(gfa_name, "w"); - ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "atg", output_file); + ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "atg", output_file); fclose(output_file); if(asm_opt.bed_inconsist_rate != 0) { diff --git a/Overlaps.h b/Overlaps.h index c9fb08c..4b870b9 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -27,6 +27,7 @@ #define DOUBLE_CHECK_THRES 0.1 #define FINAL_DOUBLE_CHECK_THRES 0.2 #define CHIMERIC_TRIM_THRES 4 +#define GAP_LEN 100 // #define PRIMARY_LABLE 1 // #define ALTER_LABLE 2 // #define HAP_LABLE 4 @@ -1209,6 +1210,7 @@ typedef struct{ int max_hang; int min_ovlp; int is_bench; + long long gap_fuzz; bub_label_t* b_mask_t; }ug_opt_t; @@ -1218,6 +1220,14 @@ 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, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t); uint32_t cmp_untig_graph(ma_ug_t *src, ma_ug_t *dest); +void reduce_hamming_error(asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, +int max_hang, int min_ovlp, long long gap_fuzz); +int ma_ug_seq_scaffold(ma_ug_t *g, asg_t *read_g, ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, +kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, kvec_asg_arc_t_warp *E, uint32_t is_polish); +void ma_ug_print(const ma_ug_t *ug, asg_t* read_g, const ma_sub_t *coverage_cut, +ma_hit_t_alloc* sources, R_to_U* ruIndex, const char* prefix, FILE *fp); +void ma_ug_print_simple(const ma_ug_t *ug, asg_t* read_g, const ma_sub_t *coverage_cut, +ma_hit_t_alloc* sources, R_to_U* ruIndex, const char* prefix, FILE *fp); #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/hic.cpp b/hic.cpp index d7f8840..aad0c24 100644 --- a/hic.cpp +++ b/hic.cpp @@ -15536,7 +15536,7 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o skip_flipping: verbose_het_stat(&bub); - horder_t *ho = init_horder_t(&sl.hits, idx->uID_bits, idx->pos_mode, idx->read_g, idx->ug, &bub, opt); + horder_t *ho = init_horder_t(&sl.hits, idx->uID_bits, idx->pos_mode, idx->read_g, idx->ug, &bub, opt, 3); ///print_hc_links(&link, 0, &hap); // print_kv_u_trans(&k_trans, &link, s->s); diff --git a/horder.cpp b/horder.cpp index 6b7567d..76f1486 100644 --- a/horder.cpp +++ b/horder.cpp @@ -47,7 +47,7 @@ KRADIX_SORT_INIT(osg, osg_arc_t, osg_arc_key, member_size(osg_arc_t, u)) #define BREAK_THRES 5000000 #define BREAK_CUTOFF 0.1 #define BREAK_BOUNDARY 0.015 -#define GAP_LEN 100 + typedef struct { uint64_t ruid; @@ -60,6 +60,16 @@ typedef struct { kvec_t(uint64_t) idx; } u_hits_t; +typedef struct { + uint64_t e; + double w; +} hw_aux_t; + +typedef struct { + hw_aux_t *a; + size_t n, m; +} h_w_t; + typedef struct { uint64_t s, e, dp; } h_cov_t; @@ -633,7 +643,7 @@ void get_consensus_break(h_covs *res, h_covs *tmp) } } if(tmp->n == 0) fprintf(stderr, "ERROR-break-0\n"); - if(tmp->n == 1) return; + // if(tmp->n == 1) return; for (i = m = 0; i < res->n; ++i) { @@ -872,8 +882,10 @@ uint64_t get_utg_len(ma_ug_t *ug) } void break_utg_horder(horder_t *h, h_covs *b_points) { + if(b_points->n == 0) return; + kvec_t(uint64_t) join; kv_init(join); ma_ug_t *ug = h->ug; - uint64_t k, l, i, idx, m, pidx, de_u, u_n; + uint64_t k, l, i, idx, m, pidx, de_u, u_n, oug_n = ug->u.n, dug_n = 0, puid, nuid[2], ps, pe; radix_sort_h_cov_s(b_points->a, b_points->a+b_points->n); for (k = 1, l = 0; k <= b_points->n; ++k) @@ -888,7 +900,11 @@ void break_utg_horder(horder_t *h, h_covs *b_points) idx = b_points->a[i].e + 1; if(idx > pidx && idx - pidx < u_n) { - de_u |= append_sub_utg(h, b_points->a[l].s, pidx, idx); + if(append_sub_utg(h, b_points->a[l].s, pidx, idx)) + { + de_u++; + kv_push(uint64_t, join, (b_points->a[l].s<<32)|(ug->u.n-1)); + } } pidx = idx; } @@ -896,23 +912,39 @@ void break_utg_horder(horder_t *h, h_covs *b_points) idx = u_n; if(idx > pidx && idx - pidx < u_n) { - de_u |= append_sub_utg(h, b_points->a[l].s, pidx, idx); + if(append_sub_utg(h, b_points->a[l].s, pidx, idx)) + { + de_u++; + kv_push(uint64_t, join, (b_points->a[l].s<<32)|(ug->u.n-1)); + } } if(de_u) { free(ug->u.a[b_points->a[l].s].a); free(ug->u.a[b_points->a[l].s].s); memset(&(ug->u.a[b_points->a[l].s]), 0, sizeof(ug->u.a[b_points->a[l].s])); + dug_n++; } l = k; } } + // fprintf(stderr, "oug_n-%lu, dug_n-%lu\n", oug_n, dug_n); + // for (i = 0; i < join.n; i++) + // { + // fprintf(stderr, "+puid-%lu, nuid-%u\n", join.a[i]>>32, (uint32_t)join.a[i]); + // } + + for (i = 0; i < join.n; i++) + { + join.a[i] -= dug_n; + } for (i = m = 0; i < ug->u.n; i++) { if(!ug->u.a[i].a) continue; + if(i < oug_n) kv_push(uint64_t, join, (i<<32)|(m)); ug->u.a[m] = ug->u.a[i]; m++; } @@ -925,6 +957,63 @@ void break_utg_horder(horder_t *h, h_covs *b_points) } ug->u.n = m; } + + // for (i = 0; i < join.n; i++) + // { + // fprintf(stderr, "-puid-%lu, nuid-%u\n", join.a[i]>>32, (uint32_t)join.a[i]); + // } + + oug_n = h->avoid.n; + radix_sort_ho64(join.a, join.a+join.n); + for (k = 1, l = 0; k <= join.n; ++k) + { + if (k == join.n || ((join.a[k]>>32) != (join.a[l]>>32))) + { + puid = (join.a[l]>>32); + nuid[0] = (uint32_t)join.a[l]; + nuid[0] <<= 1; + + nuid[1] = (uint32_t)join.a[k-1]; + nuid[1] <<= 1; nuid[1] += 1; + + for (i = 0; i < oug_n; i++) + { + ps = h->avoid.a[i]>>32; + pe = (uint32_t)h->avoid.a[i]; + + if((ps>>1) == puid) ps = nuid[ps&1]; + if((pe>>1) == puid) pe = nuid[pe&1]; + + h->avoid.a[i] = (ps<<32)|pe; + } + + if(k - l > 1) + { + for (i = l; i + 1 < k; i++) + { + nuid[0] = (uint32_t)join.a[i]; + nuid[0] <<=1; nuid[0] += 1; + + nuid[1] = (uint32_t)join.a[i+1]; + nuid[1] <<=1; + kv_push(uint64_t, h->avoid, (nuid[0]<<32)|(nuid[1])); + } + } + + l = k; + } + } + + radix_sort_ho64(h->avoid.a, h->avoid.a + h->avoid.n); + + // for (i = 0; i < h->avoid.n; i++) + // { + // fprintf(stderr, "break-s-%lu (dir: %lu), break-e-%u (dir: %u)\n", + // h->avoid.a[i]>>33, (h->avoid.a[i]>>32)&1, + // ((uint32_t)h->avoid.a[i])>>1, ((uint32_t)h->avoid.a[i])&1); + // } + + kv_destroy(join); } void break_contig(horder_t *h, uint64_t cutoff_s, uint64_t cutoff_e) @@ -1076,12 +1165,42 @@ uint64_t rs, uint64_t re, uint64_t limit_s, uint64_t limit_e, int unique_only) fprintf(stderr, "******cnt-%lu, cnt_no_lim-%lu, rs-%lu, re-%lu, limit_s-%lu, limit_e-%lu\n", cnt, cnt_no_lim, rs, re, limit_s, limit_e); } +uint64_t get_sub_cov(kvec_pe_hit *hits, uint64_t sidx, uint64_t eidx, uint64_t ulen, +uint64_t rs, uint64_t re, uint64_t limit_s, uint64_t limit_e, int unique_only) +{ + uint64_t i, p0s, p0e, p1s, p1e, span_s, span_e, cnt = 0; + + + for (i = sidx; i < eidx; i++)///keep all hic hits that contain interval we want + { + if(get_hit_suid(*hits, i) != get_hit_euid(*hits, i)) continue; + if(unique_only && hits->a.a[i].id == 0) continue; + p0s = get_hit_spos(*hits, i); + p0e = get_hit_spos_e(*hits, i); + p1s = get_hit_epos(*hits, i); + p1e = get_hit_epos_e(*hits, i); + + span_s = MIN(MIN(p0s, p0e), MIN(p1s, p1e)); + span_s = MIN(span_s, ulen-1); + span_e = MAX(MAX(p0s, p0e), MAX(p1s, p1e)); + span_e = MIN(span_e, ulen-1) + 1; + + //if(span_e - span_s <= ulen*BREAK_CUTOFF)//need it or not? + if(span_s <= rs && span_e >= re) + { + if(span_s >= limit_s && span_e <= limit_e) cnt++; + } + } + + return cnt; +} + void detect_lowNs(kvec_pe_hit *hit, uint64_t sHit, uint64_t eHit, kvec_t_u64_warp *b, h_cov_t *Np, uint64_t len, uint64_t cutoff_s, uint64_t cutoff_e, h_covs *res, h_covs *cov_buf, h_covs *b_points, uint64_t local_bound, int unique_only) { uint64_t cov_hic, cov_utg, cov_ava, i, p0s, p0e, p1s, p1e, span_s, span_e, cutoff, bs, be, occ = 0; - uint64_t sPos, ePos; + uint64_t sPos, ePos, min_cutoff; h_cov_t *p = NULL; b->a.n = 0; cov_hic = cov_utg = 0; sPos = (Np->s>=local_bound? Np->s-local_bound:0); @@ -1154,14 +1273,16 @@ h_covs *cov_buf, h_covs *b_points, uint64_t local_bound, int unique_only) // fprintf(stderr, "consensus_break-s: %lu, e: %lu\n", cov_buf->a[i].s, cov_buf->a[i].e); // } /*******************************for debug************************************/ - for (i = 0; i < cov_buf->n; i++) + for (i = 0, min_cutoff = (uint64_t)-1; i < cov_buf->n; i++) { if(cov_buf->a[i].s<=Np->s && cov_buf->a[i].e>=Np->e) { break; } + min_cutoff = MIN(min_cutoff, cov_buf->a[i].dp); } - if(i < cov_buf->n) + if(i < cov_buf->n || + get_sub_cov(hit, sHit, eHit, len, Np->s, Np->e, sPos, ePos, unique_only) <= min_cutoff) { kv_pushp(h_cov_t, *b_points, &p); p->s = get_hit_suid(*hit, sHit); p->e = Np->dp; p->dp = 0; @@ -1351,7 +1472,7 @@ double get_max_weight(uint32_t u, uint32_t v, osg_t *g) void update_scg(horder_t *h) { - uint64_t i, k, l, p0s, p0e, p1s, p1e, span_s, span_e, suid, euid, v, w, slen, elen, *ep = NULL; + uint64_t i, k, l, p0s, p0e, p1s, p1e, span_s, span_e, suid, euid, v, w, slen, elen, *ep = NULL, dis; uint64_t t_hits = 0, a_hits = 0; double div, max_div; kvec_t(uint64_t) e; kv_init(e); @@ -1416,23 +1537,37 @@ void update_scg(horder_t *h) { if (k == e.n || e.a[k] != e.a[l]) { - div = h->sg.g->seq[e.a[l]>>33].ez[(e.a[l]>>32)&1] + - h->sg.g->seq[((uint32_t)e.a[l])>>1].ez[e.a[l]&1]; - p = osg_arc_pushp(h->sg.g); - p->u = p->v = p->occ = p->del = p->w = p->nw = 0; - p->u = e.a[l]>>32; p->v = (uint32_t)e.a[l]; - p->occ = k - l; - if(div != 0) p->w = (double)(k - l)*(max_div/div); - p = osg_arc_pushp(h->sg.g); - p->u = p->v = p->occ = p->del = p->w = p->nw = 0; - p->u = (uint32_t)e.a[l]; p->v = e.a[l]>>32; - p->occ = k - l; - if(div != 0) p->w = (double)(k - l)*(max_div/div); + for (i = 0; i < h->avoid.n; i++) + { + v = e.a[l]; + w = e.a[l]<<32; w |= (e.a[l]>>32); + if(h->avoid.a[i] == v || h->avoid.a[i] == w) + { + break; + } + } + + if(i >= h->avoid.n) + { + div = h->sg.g->seq[e.a[l]>>33].ez[(e.a[l]>>32)&1] + + h->sg.g->seq[((uint32_t)e.a[l])>>1].ez[e.a[l]&1]; + p = osg_arc_pushp(h->sg.g); + p->u = p->v = p->occ = p->del = p->w = p->nw = 0; + p->u = e.a[l]>>32; p->v = (uint32_t)e.a[l]; + p->occ = k - l; + if(div != 0) p->w = (double)(k - l)*(max_div/div); + p = osg_arc_pushp(h->sg.g); + p->u = p->v = p->occ = p->del = p->w = p->nw = 0; + p->u = (uint32_t)e.a[l]; p->v = e.a[l]>>32; + p->occ = k - l; + if(div != 0) p->w = (double)(k - l)*(max_div/div); + + h->sg.g->seq[e.a[l]>>33].mw[(e.a[l]>>32)&1] + = MAX(h->sg.g->seq[e.a[l]>>33].mw[(e.a[l]>>32)&1], p->w); + h->sg.g->seq[((uint32_t)e.a[l])>>1].mw[e.a[l]&1] + = MAX(h->sg.g->seq[((uint32_t)e.a[l])>>1].mw[e.a[l]&1], p->w); + } - h->sg.g->seq[e.a[l]>>33].mw[(e.a[l]>>32)&1] - = MAX(h->sg.g->seq[e.a[l]>>33].mw[(e.a[l]>>32)&1], p->w); - h->sg.g->seq[((uint32_t)e.a[l])>>1].mw[e.a[l]&1] - = MAX(h->sg.g->seq[((uint32_t)e.a[l])>>1].mw[e.a[l]&1], p->w); l = k; } } @@ -1458,6 +1593,14 @@ void update_scg(horder_t *h) __func__, h->sg.g->n_seq, h->sg.g->n_arc, eg_edges, t_hits, a_hits); /*******************************for debug************************************/ + for (i = 0; i < h->sg.g->n_arc; i++) + { + p = &(h->sg.g->arc[i]); + fprintf(stderr, "u-stg%.6ul(%c)(div:%f)\tv-stg%.6ul(%c)(div:%f)\tocc:%u\tw:%f\tnw:%f\n", + (p->u>>1)+1, "+-"[p->u&1], h->sg.g->seq[p->u>>1].ez[p->u&1], + (p->v>>1)+1, "+-"[p->v&1], h->sg.g->seq[p->v>>1].ez[p->v&1], p->occ, p->w, p->nw); + } + fprintf(stderr, "sbsbsbsb\n\n\n\n\n\n"); // uint32_t u, nv, f; // osg_arc_t *av = NULL; // for (k = 0; k < h->sg.g->n_arc; k++) @@ -1993,9 +2136,6 @@ void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis) uint32_t *idx = NULL; MALLOC(idx, h->sg.g->n_seq); memset(idx, -1, sizeof(uint32_t)*h->sg.g->n_seq); - // fprintf(stderr, "***0***vis-occ: %u, sl-occ: %u\n", - // get_vis_occ(vis, h->sg.g->n_seq), get_sl_occ(sl)); - for (k = 0; k < sl->n; k++) { p = &(sl->a[k]); @@ -2005,7 +2145,7 @@ void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis) } } - + /** while (get_max_anchor(h, sl, vis, w, sgv, idx, &max_utg, &max_sc)) { @@ -2013,9 +2153,8 @@ void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis) vis[max_utg<<1] = vis[(max_utg<<1)+1] = 1; idx[max_utg] = max_sc; - - // fprintf(stderr, "max_utg-%u, max_sc-%u\n", max_utg, max_sc); } + **/ for (k = 0; k < h->sg.g->n_seq; k++) { @@ -2031,6 +2170,8 @@ void refine_layout(horder_t *h, sc_lay_t *sl, uint8_t *vis) free(w); free(idx); free(sgv); } + + void generate_scaffold(ma_utg_t *su, lay_t *ly, ma_ug_t *pug, asg_t *rg) { ma_utg_t *uu = NULL; @@ -2053,7 +2194,7 @@ void generate_scaffold(ma_utg_t *su, lay_t *ly, ma_ug_t *pug, asg_t *rg) if(i < ly->n - 2) kv_push(uint64_t, *su, (uint64_t)-1); } if(ly->n != 2) is_circle = 0; - + for (i = 0, totalLen = 0; i < su->n-1; i++) { if(su->a[i] == (uint64_t)-1) @@ -2143,12 +2284,52 @@ void generate_scaffold(ma_utg_t *su, lay_t *ly, ma_ug_t *pug, asg_t *rg) } } +uint64_t get_nuid(sc_lay_t *sl, uint64_t *p) +{ + uint64_t k, ouid[2]; + for (k = 0; k < sl->n; k++) + { + ouid[0] = sl->a[k].a[0]; + ouid[1] = sl->a[k].a[sl->a[k].n - 1]; + if((*p) == ouid[0]) + { + (*p) = (k<<1); + return 1; + } + + if((*p) == ouid[1]) + { + (*p) = (k<<1)+1; + return 1; + } + } + + return 0; +} + +void update_avoids(horder_t *h, sc_lay_t *sl) +{ + uint64_t i, m, ps, pe; + for (i = m = 0; i < h->avoid.n; i++) + { + ps = h->avoid.a[i]>>32; + pe = (uint32_t)h->avoid.a[i]; + if(get_nuid(sl, &ps) && get_nuid(sl, &pe)) + { + h->avoid.a[m] = (ps<<32)|pe; + m++; + } + } + h->avoid.n = m; +} + void update_ug_by_layout(horder_t *h, sc_lay_t *sl) { uint32_t i; lay_t *p = NULL; ma_utg_t *pu = NULL; ma_ug_t *sug = NULL; + kvec_t_u64_warp n_avoids; kv_init(n_avoids.a); sug = (ma_ug_t*)calloc(1, sizeof(ma_ug_t)); for (i = 0; i < sl->n; i++) { @@ -2158,6 +2339,74 @@ void update_ug_by_layout(horder_t *h, sc_lay_t *sl) } ma_ug_destroy(h->ug); h->ug = sug; + kv_destroy(n_avoids.a); + update_avoids(h, sl); +} + +void get_long_switch_scaffolds(horder_t *h, sc_lay_t *sl, osg_t *lg) +{ + fprintf(stderr, "\n[M::%s::]\n", __func__); + uint32_t i, k, r_i, ori, sw[3], sw_inner[3]; + uint64_t v, len; + kvec_t(uint64_t) idx; kv_init(idx); + lay_t *p = NULL; + ma_utg_t *u = NULL; + for (i = 0; i < sl->n; i++) + { + p = &(sl->a[i]); + sw[0] = sw[1] = sw[2] = 0; + for (k = 0; k < p->n; k += 2) + { + u = &(h->ug->u.a[p->a[k]>>1]); + ori = p->a[k]&1; + for (r_i = 0; r_i < u->n; r_i++) + { + v = (ori?u->a[u->n - r_i - 1]:u->a[r_i]); + if(v != (uint64_t)-1) + { + v >>= 33; + sw[R_INF.trio_flag[v]]++; + } + } + } + + v = MIN(sw[FATHER], sw[MOTHER]); + v = ((uint32_t)-1) - v; + v <<= 32; v |= i; + kv_push(uint64_t, idx, v); + } + + radix_sort_ho64(idx.a, idx.a+idx.n); + for (i = 0; i < sl->n; i++) + { + p = &(sl->a[(uint32_t)(idx.a[i])]); + fprintf(stderr, "\nscaf-%u-th, occ-%u\n", (uint32_t)(idx.a[i]), (uint32_t)(p->n>>1)); + sw[0] = sw[1] = sw[2] = len = 0; + for (k = 0; k < p->n; k += 2) + { + sw_inner[0] = sw_inner[1] = sw_inner[2] = 0; + u = &(h->ug->u.a[p->a[k]>>1]); + ori = p->a[k]&1; + for (r_i = 0; r_i < u->n; r_i++) + { + v = (ori?u->a[u->n - r_i - 1]:u->a[r_i]); + if(v != (uint64_t)-1) + { + v >>= 33; + if(R_INF.trio_flag[v] == FATHER || R_INF.trio_flag[v] == MOTHER) + { + sw[R_INF.trio_flag[v]]++; + sw_inner[R_INF.trio_flag[v]]++; + } + } + } + fprintf(stderr, "utg%.6ul (ori: %u), u->len-%u, sw_in[FATHER]-%u, sw_in[MOTHER]-%u\n", + (p->a[k]>>1)+1, p->a[k]&1, u->len, sw_inner[FATHER], sw_inner[MOTHER]); + len += u->len + ((k + 2)< p->n? GAP_LEN:0); + } + fprintf(stderr, "sw[FATHER]-%u, sw[MOTHER]-%u, len-%lu\n", sw[FATHER], sw[MOTHER], len); + } + kv_destroy(idx); } void layout_scg(horder_t *h, double nw_thres, uint32_t occ_thres) @@ -2190,16 +2439,19 @@ void layout_scg(horder_t *h, double nw_thres, uint32_t occ_thres) radix_sort_osg(h->sg.g->arc, h->sg.g->arc + h->sg.g->n_arc); get_backbone_layout(h, &sl, lg, vis); + + get_long_switch_scaffolds(h, &sl, lg); refine_layout(h, &sl, vis); - print_N50_layout(h->ug, &sl); + // print_N50_layout(h->ug, &sl); update_ug_by_layout(h, &sl); print_N50(h->ug); kv_destroy(sl); + osg_destroy(lg); free(vis); } @@ -2209,28 +2461,126 @@ void renew_scaffold(horder_t *h) while (1) { update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); - if(!break_scaffold(h, /**5**/10, /**15**/20, 2500000, 1)) break; + if(!break_scaffold(h, 5, 15, 2500000, 1)) break; print_N50(h->ug); } fprintf(stderr, "[M::%s::%.3f] \n", __func__, yak_realtime()-index_time); } -horder_t *init_horder_t(kvec_pe_hit *i_hits, uint64_t i_hits_uid_bits, uint64_t i_hits_pos_mode, -asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt) +void print_scaffold(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, +ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, +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) { + char* gfa_name = (char*)malloc(strlen(output_file_name)+100); + sprintf(gfa_name, "%s.%s.p_ctg.gfa", output_file_name, "stg"); + fprintf(stderr, "Writing %s to disk... \n", gfa_name); + FILE* output_file = NULL; + output_file = fopen(gfa_name, "w"); + + + kvec_asg_arc_t_warp new_rtg_edges; + kv_init(new_rtg_edges.a); + ma_ug_seq_scaffold(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 1); + ma_ug_print(ug, sg, coverage_cut, sources, ruIndex, "stg", output_file); + fclose(output_file); + + sprintf(gfa_name, "%s.%s.p_ctg.noseq.gfa", output_file_name, "stg"); + output_file = fopen(gfa_name, "w"); + ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "stg", output_file); + fclose(output_file); + + free(gfa_name); + kv_destroy(new_rtg_edges.a); +} + +void scaffold_hap(horder_t *h, ug_opt_t *opt, uint32_t round, char *output_file_name, uint8_t flag) +{ + uint32_t i; + kv_destroy(h->u_hits.a); + kv_destroy(h->u_hits.idx); + kv_destroy(h->u_hits.occ); + memset(&(h->u_hits), 0, sizeof(h->u_hits)); + kv_destroy(h->avoid); + h->avoid.m = h->avoid.n = 0; + h->avoid.a = NULL; + osg_destroy(h->sg.g); + h->sg.g = NULL; + ma_ug_destroy(h->ug); + h->ug = NULL; + + h->ug = get_trio_unitig_graph(h->r_g, flag, opt); + asg_destroy(h->ug->g); + h->ug->g = NULL; + + print_N50(h->ug); + + for (i = 0; i < round; i++) + { + update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); + update_scg(h); + layout_scg(h, 1.001, 19); + renew_scaffold(h); + } + + char* gfa_name = (char*)malloc(strlen(output_file_name)+100); + sprintf(gfa_name, "%s.%s", output_file_name, (flag==FATHER?"hap1":"hap2")); + + print_scaffold(h->ug, h->r_g, opt->coverage_cut, gfa_name, + opt->sources, opt->reverse_sources, opt->tipsLen, opt->tip_drop_ratio, + opt->stops_threshold, opt->ruIndex, opt->chimeric_rate, opt->drop_ratio, + opt->max_hang, opt->min_ovlp); + + free(gfa_name); +} + +void output_hic_rtg(ma_ug_t *ug, asg_t *rg, ug_opt_t *opt, char* output_file_name) +{ + char* gfa_name = (char*)malloc(strlen(output_file_name)+50); + sprintf(gfa_name, "%s.all.noseq.gfa", output_file_name); + FILE* output_file = fopen(gfa_name, "w"); + ma_ug_print_simple(ug, rg, opt->coverage_cut, opt->sources, opt->ruIndex, "utg", output_file); + fclose(output_file); + free(gfa_name); +} + +horder_t *init_horder_t(kvec_pe_hit *i_hits, uint64_t i_hits_uid_bits, uint64_t i_hits_pos_mode, +asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt, uint32_t round) +{ + uint32_t i; horder_t *h = NULL; CALLOC(h, 1); get_r_hits(i_hits, &(h->r_hits), i_rg, i_ug, bub, i_hits_uid_bits, i_hits_pos_mode); h->r_g = copy_read_graph(i_rg); horder_clean_sg_by_utg(h->r_g, i_ug); + // output_hic_rtg(i_ug, h->r_g, opt, asm_opt.output_file_name); + + // reduce_hamming_error(h->r_g, opt->sources, opt->coverage_cut, opt->max_hang, opt->min_ovlp, opt->gap_fuzz); + /** + scaffold_hap(h, opt, round, asm_opt.output_file_name, FATHER); + scaffold_hap(h, opt, round, asm_opt.output_file_name, MOTHER); + **/ + + generate_haplotypes(h, opt); print_N50(h->ug); - update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); - // break_contig(h, 10, 20); - print_N50(h->ug); - update_scg(h); - layout_scg(h, 1.001, 19); - renew_scaffold(h); + // update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); + // break_contig(h, 10, 20); + for (i = 0; i < round; i++) + { + update_u_hits(&(h->u_hits), &(h->r_hits), h->ug, h->r_g); + update_scg(h); + layout_scg(h, 1.001, 19); + renew_scaffold(h); + } + + + print_scaffold(h->ug, h->r_g, opt->coverage_cut, asm_opt.output_file_name, + opt->sources, opt->reverse_sources, opt->tipsLen, opt->tip_drop_ratio, + opt->stops_threshold, opt->ruIndex, opt->chimeric_rate, opt->drop_ratio, + opt->max_hang, opt->min_ovlp); + + exit(1); return h; } @@ -2244,6 +2594,8 @@ void destory_horder_t(horder_t **h) kv_destroy((*h)->u_hits.idx); kv_destroy((*h)->u_hits.occ); + kv_destroy((*h)->avoid); + osg_destroy((*h)->sg.g); ma_ug_destroy((*h)->ug); diff --git a/horder.h b/horder.h index 67738f6..8e10fde 100644 --- a/horder.h +++ b/horder.h @@ -29,6 +29,7 @@ typedef struct { osg_t *g; }scg_t; typedef struct { + kvec_t(uint64_t) avoid; kvec_pe_hit r_hits, u_hits; ma_ug_t *ug; asg_t *r_g; @@ -36,6 +37,6 @@ typedef struct { }horder_t; horder_t *init_horder_t(kvec_pe_hit *i_hits, uint64_t i_hits_uid_bits, uint64_t i_hits_pos_mode, -asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt); +asg_t *i_rg, ma_ug_t* i_ug, bubble_type* bub, ug_opt_t *opt, uint32_t round); void destory_horder_t(horder_t **h); #endif