diff --git a/Assembly.cpp b/Assembly.cpp index 673c2df..9479b1a 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1203,8 +1203,8 @@ int ha_assemble(void) ha_extract_print_list(&R_INF, asm_opt.extract_iter, asm_opt.extract_list); exit(0); } - if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN) && !(asm_opt.flag & HA_F_VERBOSE_GFA)) ha_triobin(&asm_opt); - ///if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt); + if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN) && !(asm_opt.flag & HA_F_VERBOSE_GFA)) ha_triobin(&asm_opt), ovlp_loaded = 2; + ///if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt), ovlp_loaded = 2; if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads(); if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF(); if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2; diff --git a/CommandLines.cpp b/CommandLines.cpp index 853ea98..f0265a6 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -117,6 +117,9 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->purge_level_trio = 0; asm_opt->purge_simi_rate = 0.75; asm_opt->purge_overlap_len = 1; + asm_opt->recover_atg_cov_min = -1; + asm_opt->recover_atg_cov_max = -1; + asm_opt->hom_global_coverage = -1; } void destory_opt(hifiasm_opt_t* asm_opt) diff --git a/CommandLines.h b/CommandLines.h index b0e4571..60ee2b1 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -58,6 +58,9 @@ typedef struct { int purge_level_primary; int purge_level_trio; int purge_overlap_len; + int recover_atg_cov_min; + int recover_atg_cov_max; + int hom_global_coverage; float max_hang_rate; float min_drop_rate; @@ -71,6 +74,7 @@ typedef struct { long long num_recorrected_bases; long long mem_buf; long long coverage; + } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Correct.h b/Correct.h index 944583a..ee8912d 100644 --- a/Correct.h +++ b/Correct.h @@ -5,6 +5,7 @@ #include "Levenshtein_distance.h" #include "POA.h" #include "Process_Read.h" +#include "Correct.h" //#define CORRECT_THRESHOLD 0.70 #define CORRECT_THRESHOLD 0.60 diff --git a/Overlaps.cpp b/Overlaps.cpp index 032e3d0..c8d3103 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -21231,7 +21231,7 @@ void unroll_simple_case_advance(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* reve while (n_reduce > 0) { n_reduce = 0; - + ///break nearly circle, forget why... n_reduce += asg_arc_del_simple_circle_untig(NULL, NULL, nsg, 100, 0); /** for (v = 0; v < n_vtx; ++v) @@ -22039,6 +22039,205 @@ R_to_U* ruIndex) } + +void append_utg(ma_ug_t* ptg, ma_ug_t* atg) +{ + uint64_t num_nodes = 0; + asg_t* nsg = atg->g; + uint32_t v, n_vtx = nsg->n_seq; + ma_utg_t *p; + for (v = 0; v < n_vtx; ++v) + { + if(nsg->seq[v].del || atg->u.a[v].m == 0) continue; + num_nodes++; + } + + if(num_nodes == 0) return; + + ptg->u.n = ptg->u.n + num_nodes; + if(ptg->u.n > ptg->u.m) + { + ptg->u.m = ptg->u.n; + ptg->u.a = (ma_utg_t*)realloc(ptg->u.a, ptg->u.m*sizeof(ma_utg_t)); + } + ptg->u.n = ptg->u.n - num_nodes; + + for (v = 0; v < atg->g->n_seq; ++v) + { + if(atg->g->seq[v].del || atg->u.a[v].m == 0) continue; + + p = &(ptg->u.a[ptg->u.n]); + p->len = atg->u.a[v].len; + p->circ = atg->u.a[v].circ; + p->start = atg->u.a[v].start; + p->end = atg->u.a[v].end; + p->m = atg->u.a[v].m; atg->u.a[v].m = 0; + p->n = atg->u.a[v].n; atg->u.a[v].n = 0; + p->a = atg->u.a[v].a; atg->u.a[v].a = 0; + p->s = atg->u.a[v].s; atg->u.a[v].s = 0; + asg_seq_set(ptg->g, ptg->u.n, p->len, 0); + ptg->u.n++; + } + + asg_cleanup(ptg->g); +} + + +void print_utg_coverage(ma_ug_t *ug, ma_sub_t* coverage_cut, uint32_t v, ma_hit_t_alloc* sources) +{ + asg_t* nsg = ug->g; + uint32_t rId, k, j; + ma_utg_t* u = NULL; + ma_hit_t *h; + + if(nsg->seq[v].del) return; + u = &(ug->u.a[v]); + if(u->m == 0) return; + long long R_bases = 0, C_bases = 0; + long long U_R_bases = 0, U_C_bases = 0; + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + C_bases = 0; + R_bases = coverage_cut[rId].e - coverage_cut[rId].s; + for (j = 0; j < (uint64_t)(sources[rId].length); j++) + { + h = &(sources[rId].buffer[j]); + if(h->el != 1) continue; + C_bases += Get_qe((*h)) - Get_qs((*h)); + } + U_R_bases += R_bases; + U_C_bases += C_bases; + C_bases = C_bases/R_bases; + + fprintf(stderr, "%.*s\t%lld\n", (int)Get_NAME_LENGTH(R_INF, rId), Get_NAME(R_INF, rId), C_bases); + } + + fprintf(stderr, "v: %u, coverage: %lld\n\n", v, U_C_bases/U_R_bases); +} + +void recover_utg_by_coverage(ma_ug_t **ptg, asg_t* read_g, ma_sub_t* coverage_cut, +ma_hit_t_alloc* sources, R_to_U* ruIndex) +{ + if(asm_opt.recover_atg_cov_min == -1) return; + if(asm_opt.recover_atg_cov_max == -1) return; + if(asm_opt.recover_atg_cov_min > asm_opt.recover_atg_cov_max) return; + ma_ug_t *atg = NULL; + atg = ma_ug_gen_primary(read_g, ALTER_LABLE); + asg_t* nsg = atg->g; + uint32_t v, n_vtx = nsg->n_seq, k, j, rId, available_reads = 0, keep_atg = 0, tn, is_Unitig; + ma_utg_t* u = NULL; + ma_hit_t *h; + long long R_bases = 0, C_bases = 0, C_bases_primary = 0, C_bases_alter = 0; + for (v = 0; v < n_vtx; ++v) + { + if(nsg->seq[v].del) continue; + u = &(atg->u.a[v]); + if(u->m == 0) continue; + available_reads = 0; + + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + C_bases = C_bases_primary = C_bases_alter = 0; + R_bases = coverage_cut[rId].e - coverage_cut[rId].s; + for (j = 0; j < (uint64_t)(sources[rId].length); j++) + { + h = &(sources[rId].buffer[j]); + if(h->el != 1) continue; + tn = Get_tn((*h)); + if(read_g->seq[tn].del == 1) + { + ///get the id of read that contains it + get_R_to_U(ruIndex, tn, &tn, &is_Unitig); + if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue; + } + if(read_g->seq[tn].del == 1) continue; + if(read_g->seq[tn].c == ALTER_LABLE) + { + C_bases_alter += Get_qe((*h)) - Get_qs((*h)); + } + else + { + C_bases_primary += Get_qe((*h)) - Get_qs((*h)); + } + } + + C_bases = C_bases_primary + C_bases_alter; + if(C_bases_alter < C_bases * ALTER_COV_THRES) continue; + + C_bases = C_bases/R_bases; + if(C_bases >= asm_opt.recover_atg_cov_min && C_bases <= asm_opt.recover_atg_cov_max) + { + available_reads++; + } + } + + if(available_reads < (u->n * 0.8) || available_reads == 0) + { + asg_seq_del(nsg, v); + + if(u->m!=0) + { + u->m = u->n = 0; + free(u->a); + u->a = NULL; + } + } + else + { + ///print_utg_coverage(atg, coverage_cut, v, sources); + ///fprintf(stderr, "rId: %u\n", rId); + ///fprintf(stderr, "%.*s\t%lld\n", (int)Get_NAME_LENGTH(R_INF, rId), Get_NAME(R_INF, rId), C_bases); + keep_atg++; + } + } + + if(keep_atg > 0) + { + asg_cleanup(nsg); + asg_symm(nsg); + append_utg(*ptg, atg); + + + n_vtx = read_g->n_seq; + for (v = 0; v < n_vtx; v++) + { + read_g->seq[v].c = ALTER_LABLE; + } + + nsg = (*ptg)->g; + n_vtx = nsg->n_seq; + for (v = 0; v < n_vtx; ++v) + { + if(nsg->seq[v].del) continue; + u = &((*ptg)->u.a[v]); + if(u->m == 0) continue; + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + read_g->seq[rId].c = nsg->seq[v].c; + } + } + + + n_vtx = read_g->n_seq; + for (v = 0; v < n_vtx; v++) + { + if(read_g->seq[v].c == ALTER_LABLE) + { + asg_seq_drop(read_g, v); + } + } + } + + fprintf(stderr, "keep_atg: %u\n", keep_atg); + ma_ug_destroy(atg); +} + + + + void adjust_utg_by_primary(ma_ug_t **ug, asg_t* read_g, float drop_rate, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, @@ -22049,7 +22248,8 @@ kvec_asg_arc_t_warp* new_rtg_edges) uint32_t v, n_vtx = nsg->n_seq, k, rId, just_contain; ma_utg_t* u = NULL; - + ///print_utg_coverage(*ug, coverage_cut, 440, sources); + ///exit(0); /** kvec_t_u32_warp new_rtg_nodes; kv_init(new_rtg_nodes.a); @@ -22110,6 +22310,7 @@ kvec_asg_arc_t_warp* new_rtg_edges) renew_utg(ug, read_g, new_rtg_edges); } } + n_vtx = read_g->n_seq; @@ -22134,7 +22335,6 @@ kvec_asg_arc_t_warp* new_rtg_edges) } } - n_vtx = read_g->n_seq; for (v = 0; v < n_vtx; v++) { @@ -22143,7 +22343,17 @@ kvec_asg_arc_t_warp* new_rtg_edges) asg_seq_drop(read_g, v); } } - + + if(asm_opt.recover_atg_cov_min == -1 || asm_opt.recover_atg_cov_max == -1) + { + asm_opt.recover_atg_cov_max = asm_opt.hom_global_coverage/HOM_PEAK_RATE; + asm_opt.recover_atg_cov_min = asm_opt.recover_atg_cov_max * 0.8; + asm_opt.recover_atg_cov_max = asm_opt.recover_atg_cov_max * 1.2; + } + fprintf(stderr, "asm_opt.recover_atg_cov_min: %d\n", asm_opt.recover_atg_cov_min); + fprintf(stderr, "asm_opt.recover_atg_cov_max: %d\n", asm_opt.recover_atg_cov_max); + + recover_utg_by_coverage(ug, read_g, coverage_cut, sources, ruIndex); /** kv_destroy(new_rtg_nodes.a); diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index 843d449..8df0701 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -6,6 +6,9 @@ #include "Overlaps.h" #include "Correct.h" #include "kthread.h" +#include "kdq.h" + +KDQ_INIT(uint64_t) #define Cal_Off(OFF) ((long long)((uint32_t)((OFF)>>32)) - (long long)((uint32_t)((OFF)))) #define Get_xOff(OFF) ((long long)((uint32_t)((OFF)>>32))) @@ -149,7 +152,8 @@ void print_peak(long long* cov_buf, long long cov_buf_length, long long max_i) } -void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo_peak_cov, long long* hom_peak, long long* het_peak) +void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo_peak_cov, +long long* hom_peak, long long* het_peak, long long* k_mer_only, long long* coverage_only) { long long i, start, err_i, max_i, max2_i, max3_i, topo_peak_i, max, max2, max3, topo_peak, min; @@ -240,36 +244,71 @@ void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo { topo_peak_i = (*topo_peak_cov); topo_peak = cov_buf[topo_peak_i]; - if(topo_peak_i <= max_i * 1.2 && topo_peak_i >= max_i * 0.8 && topo_peak > max * 0.05) - { - (*het_peak) = max_i; - } - - ///fprintf(stderr, "topo_peak: %lld, topo_peak_i: %lld\n", topo_peak, topo_peak_i); + if (topo_peak <= max * 0.05) topo_peak_i = topo_peak = -1; } + long long k_mer_het, k_mer_hom, coverage_het, coverage_hom, alter_peak; + k_mer_het = k_mer_hom = coverage_het = coverage_hom = alter_peak = -1; - ///if we really want to use rev_sources for double checking, we should use peak instead of mean - ///if we found a small peak at the right hand of the largest peak + alter_peak = topo_peak_i; + k_mer_het = asm_opt.het_cov; + k_mer_hom = asm_opt.hom_cov; if(max3_i > 0) { - (*het_peak) = max_i; - (*hom_peak) = max3_i; + coverage_het = max_i; + coverage_hom = max3_i; } - else if((*het_peak) == -1) + else { - (*het_peak) = max2_i; - (*hom_peak) = max_i; + coverage_het = max2_i; + coverage_hom = max_i; } + if(k_mer_het != -1) + { + (*het_peak) = k_mer_het; + (*hom_peak) = k_mer_hom; + return; + } + else if(coverage_het != -1) + { + (*het_peak) = coverage_het; + (*hom_peak) = coverage_hom; + return; + } + else if(k_mer_hom > coverage_hom*1.5) + { + (*het_peak) = coverage_hom; + (*hom_peak) = k_mer_hom; + return; + } + else if(alter_peak != -1) + { + ///if peak is het, coverage peak is more reliable + if(coverage_hom >= alter_peak*0.8 && coverage_hom <= alter_peak*1.2) + { + (*het_peak) = coverage_hom; + return; + }///if peak is homo, k-mer peak is more reliable + else if(k_mer_hom >= alter_peak*0.8*2 && k_mer_hom <= alter_peak*1.2*2) + { + (*hom_peak) = k_mer_hom; + return; + } + } + + (*k_mer_only) = k_mer_hom; + (*coverage_only) = coverage_hom; + + // fprintf(stderr, "max: %lld, max_i: %lld\n", max, max_i); // fprintf(stderr, "max2: %lld, max2_i: %lld\n", max2, max2_i); // fprintf(stderr, "max3: %lld, max3_i: %lld\n", max3, max3_i); - fprintf(stderr, "[M::%s] Heterozygous k-mer peak: %d\n", __func__, asm_opt.het_cov); - fprintf(stderr, "[M::%s] Homozygous k-mer peak: %d\n", __func__, asm_opt.hom_cov); - fprintf(stderr, "[M::%s] Heterozygous coverage peak: %lld\n", __func__, (*het_peak)); - fprintf(stderr, "[M::%s] Homozygous coverage peak: %lld\n", __func__, (*hom_peak)); - fprintf(stderr, "[M::%s] Alter coverage peak: %lld\n", __func__, topo_peak_i); + // fprintf(stderr, "[M::%s] Heterozygous k-mer peak: %d\n", __func__, asm_opt.het_cov); + // fprintf(stderr, "[M::%s] Homozygous k-mer peak: %d\n", __func__, asm_opt.hom_cov); + // fprintf(stderr, "[M::%s] Heterozygous coverage peak: %lld\n", __func__, (*het_peak)); + // fprintf(stderr, "[M::%s] Homozygous coverage peak: %lld\n", __func__, (*hom_peak)); + // fprintf(stderr, "[M::%s] Alter coverage peak: %lld\n", __func__, topo_peak_i); } @@ -278,9 +317,7 @@ void get_read_peak(long long* cov_buf, long long cov_buf_length, long long* topo long long get_alter_peak(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, uint64_t* position_index, ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length) { - #define ALTER_COV_THRES 0.9 - #define REAL_ALTER_THRES 0.1 - + ma_utg_t* u = NULL; asg_t* nsg = ug->g; uint64_t v, j, k, qn, n_vtx = nsg->n_seq, primary_bases = 0, alter_bases = 0; @@ -372,7 +409,9 @@ ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, long long cov_buf_length) return max_i; } -long long get_read_coverage_thres(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, uint64_t* position_index, ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, uint64_t n_read, long long cov_buf_length) +long long get_read_coverage_thres(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, uint64_t* position_index, +ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, uint64_t n_read, long long cov_buf_length, +long long* k_mer_only, long long* coverage_only) { uint64_t i, j; long long* cov_buf = NULL; @@ -401,13 +440,14 @@ long long get_read_coverage_thres(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, u cov_buf_length); } - get_read_peak(cov_buf, cov_buf_length, alter_peak == -1? NULL: &alter_peak, &hom_peak, &het_peak); + get_read_peak(cov_buf, cov_buf_length, alter_peak == -1? NULL: &alter_peak, &hom_peak, &het_peak, + k_mer_only, coverage_only); free(cov_buf); - if(hom_peak != -1) return hom_peak*1.25; - if(het_peak != -1) return het_peak*2.50; - return cov_buf_length; + if(hom_peak != -1) return hom_peak*HOM_PEAK_RATE; + if(het_peak != -1) return het_peak*HET_PEAK_RATE; + return -1; } @@ -1901,7 +1941,7 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon sources, coverage_cut); ///fprintf(stderr, "ploid_coverage: %lu, cov_threshold: %lu\n", ploid_coverage, cov_threshold); - if(ploid_coverage >= cov_threshold) return NON_PLOID; + if(cov_threshold > 0 && ploid_coverage >= cov_threshold) return NON_PLOID; return PLOID; } @@ -2783,7 +2823,8 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) kvec_t_i32_warp* prevIndex_vec = &(hap_buf->buf[tid].u_buffer_prevIndex); kvec_t_i32_warp* begIndex_vec = &(hap_buf->buf[tid].u_buffer_beg); kvec_t_u8_warp* flag_vec = &(hap_buf->buf[tid].u_buffer_flag); - long long cov_threshold = hap_buf->cov_threshold; + uint64_t cov_threshold = hap_buf->cov_threshold; + if(hap_buf->cov_threshold < 0) cov_threshold = (uint64_t)-1; ma_utg_t *xReads = NULL, *yReads = NULL; ma_hit_t_alloc *xR = NULL; @@ -3927,6 +3968,207 @@ void print_all_purge_ovlp(ma_ug_t *ug, hap_overlaps_list* all_ovlp) } } + +inline int get_available_cnt(asg_t *g, uint32_t v, uint8_t* del, asg_arc_t* v_s) +{ + //v has direction + if(del && del[v>>1]) return 0; + uint32_t i, kv = 0; + asg_arc_t *av = asg_arc_a(g, v); + uint32_t nv = asg_arc_n(g, v); + + for (i = 0, kv = 0; i < nv; i++) + { + if(!av[i].del) + { + if(del && del[av[i].v>>1]) continue; + if(v_s) v_s[kv] = av[i]; + kv++; + } + } + + return kv; +} + +long long get_specific_contig_length(asg_t *g, uint8_t *del) +{ + asg_cleanup(g); + uint32_t v, n_vtx = g->n_seq * 2, q_occ; + uint8_t *mark = NULL; + ///is a queue + //kdq_t(uint64_t) *q; + ///each node has two directions + //q = kdq_init(uint64_t); + + + mark = (uint8_t*)calloc(n_vtx, 1); + + long long totalLen = 0; + for (v = 0; v < n_vtx; ++v) + { + uint32_t w, x, l, start, end, len; + asg_arc_t arc; + if (g->seq[v>>1].del || mark[v]) continue; + if (get_available_cnt(g, v, del, NULL) == 0 && get_available_cnt(g, (v^1), del, NULL) != 0) continue; + if (del[v>>1]) continue; + + mark[v] = 1; + //q->count = 0, start = v, end = v^1, len = 0; + q_occ =0, start = v, end = v^1, len = 0; + // forward + w = v; + + + while (1) + { + /** + * w----->x + * w<-----x + * that means the only suffix of w is x, and the only prefix of x is w + **/ + if (get_available_cnt(g, w, del, NULL) != 1) break; + get_available_cnt(g, w, del, &arc); + x = arc.v; // w->x + if (get_available_cnt(g, x^1, del, NULL) != 1) break; + + /** + * another direction of w would be marked as used (since w has been used) + **/ + mark[x] = mark[w^1] = 1; + ///l is the edge length, instead of overlap length + ///note: edge length is different with overlap length + ///l = asg_arc_len(arc_first(g, w)); + get_available_cnt(g, w, del, &arc); + l = ((uint32_t)((arc).ul)); + //kdq_push(uint64_t, q, (uint64_t)w<<32 | l); + q_occ++; + end = x^1, len += l; + w = x; + if (x == v) break; + } + + + //if (start != (end^1) || kdq_size(q) == 0) { // linear unitig + if (start != (end^1) || q_occ == 0) { // linear unitig + ///length of seq, instead of edge + l = g->seq[end>>1].len; + //kdq_push(uint64_t, q, (uint64_t)(end^1)<<32 | l); + q_occ++; + len += l; + } else { // circular unitig + start = end = UINT32_MAX; + goto add_unitig; // then it is not necessary to do the backward + } + + // backward + x = v; + while (1) { // similar to forward but not the same + if (get_available_cnt(g, x^1, del, NULL) != 1) break; + get_available_cnt(g, x^1, del, &arc); + w = arc.v ^ 1; + if (get_available_cnt(g, w, del, NULL) != 1) break; + mark[x] = mark[w^1] = 1; + ///l = asg_arc_len(arc_first(g, w)); + get_available_cnt(g, w, del, &arc); + l = ((uint32_t)((arc).ul)); + ///w is the seq id + direction, l is the length of edge + ///push element to the front of a queue + //kdq_unshift(uint64_t, q, (uint64_t)w<<32 | l); + q_occ++; + + start = w, len += l; + x = w; + } + + + add_unitig: + if (start != UINT32_MAX) mark[start] = mark[end] = 1; + totalLen += len; + } + //kdq_destroy(uint64_t, q); + return totalLen; +} + + +void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen) +{ + uint8_t *del = (uint8_t *)malloc(sizeof(uint8_t)*g->n_seq); + uint32_t v, k; + ma_utg_t* u = NULL; + memset(del, 1, g->n_seq); + (*primaryLen) = (*alterLen) = 0; + + for (v = 0; v < ug->g->n_seq; ++v) + { + if(ug->g->seq[v].del) continue; + if(ug->g->seq[v].c == ALTER_LABLE) continue; + u = &(ug->u.a[v]); + if(u->m == 0) continue; + for (k = 0; k < u->n; k++) + { + del[u->a[k]>>33] = 0; + } + } + (*primaryLen) = get_specific_contig_length(g, del); + + + for (v = 0; v < g->n_seq; ++v) + { + del[v] = 1 - del[v]; + } + + (*alterLen) = get_specific_contig_length(g, del); + + free(del); +} + + +int if_ploid_sample(ma_ug_t *ug, asg_t *read_g, R_to_U* ruIndex, +ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, +hap_alignment_struct_pip* hap_buf, hap_overlaps_list* all_ovlp, hap_overlaps_list* back_all_ovlp, +uint32_t minLen, double purge_threshold) +{ + asg_t* nsg = ug->g; + uint64_t v, k, total_bases = 0, alter_bases = 0, primary_bases = 0, purge_bases = 0; + kt_for(asm_opt.thread_num, hap_alignment_advance_worker, hap_buf, nsg->n_seq); + filter_hap_overlaps_by_length(all_ovlp, minLen); + normalize_hap_overlaps_advance(all_ovlp, back_all_ovlp, ug, read_g, reverse_sources, ruIndex); + + get_contig_length(ug, read_g, &primary_bases, &alter_bases); + total_bases = primary_bases + alter_bases; + fprintf(stderr, "primary_bases: %lu\n", primary_bases); + fprintf(stderr, "alter_bases: %lu\n", alter_bases); + fprintf(stderr, "total_bases: %lu\n", total_bases); + + + for (v = 0; v < all_ovlp->num; v++) + { + for (k = 0; k < all_ovlp->x[v].a.n; k++) + { + purge_bases += all_ovlp->x[v].a.a[k].x_end_pos - all_ovlp->x[v].a.a[k].x_beg_pos; + } + } + purge_bases = purge_bases/2; + fprintf(stderr, "purge_bases: %lu\n", purge_bases); + alter_bases = alter_bases + purge_bases; + fprintf(stderr, "new alter_bases: %lu\n", alter_bases); + + + for (v = 0; v < all_ovlp->num; v++) + { + all_ovlp->x[v].a.n = 0; + } + + for (v = 0; v < back_all_ovlp->num; v++) + { + back_all_ovlp->x[v].a.n = 0; + } + + if(alter_bases > total_bases * purge_threshold) return 1; + return 0; +} + + void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, @@ -3961,9 +4203,11 @@ uint32_t just_contain) asg_arc_t* p = NULL; int r; hap_alignment_struct_pip hap_buf; + long long k_mer_only, coverage_only; - hap_buf.cov_threshold = get_read_coverage_thres(ug, read_g, ruIndex, position_index, sources, coverage_cut, read_g->n_seq, COV_COUNT); - fprintf(stderr, "cov_threshold: %lld\n", hap_buf.cov_threshold); + hap_buf.cov_threshold = get_read_coverage_thres(ug, read_g, ruIndex, position_index, + sources, coverage_cut, read_g->n_seq, COV_COUNT, &k_mer_only, &coverage_only); + ///fprintf(stderr, "cov_threshold: %lld\n", hap_buf.cov_threshold); @@ -3998,6 +4242,24 @@ uint32_t just_contain) sources, reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp, 0.05, &all_ovlp); + if(hap_buf.cov_threshold < 0) + { + if(if_ploid_sample(ug, read_g, ruIndex, sources, reverse_sources, coverage_cut, + &hap_buf, &all_ovlp, &back_all_ovlp, purege_minLen, 0.333)) + { + ///if peak is het, coverage peak is more reliable + hap_buf.cov_threshold = coverage_only * HET_PEAK_RATE; + } + else + { + ///if peak is homo, k-mer peak is more reliable + hap_buf.cov_threshold = k_mer_only * HOM_PEAK_RATE; + } + } + asm_opt.hom_global_coverage = hap_buf.cov_threshold; + fprintf(stderr, "cov_threshold: %lld\n", hap_buf.cov_threshold); + + ///kt_for(asm_opt.thread_num, hap_alignment_worker, &hap_buf, nsg->n_seq); kt_for(asm_opt.thread_num, hap_alignment_advance_worker, &hap_buf, nsg->n_seq); diff --git a/Purge_Dups.h b/Purge_Dups.h index 382c181..4ea05a7 100644 --- a/Purge_Dups.h +++ b/Purge_Dups.h @@ -7,6 +7,10 @@ #include "Overlaps.h" #include "Hash_Table.h" #define COV_COUNT 1024 +#define HOM_PEAK_RATE 1.25 +#define HET_PEAK_RATE (HOM_PEAK_RATE*2) +#define ALTER_COV_THRES 0.9 +#define REAL_ALTER_THRES 0.1 void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, @@ -14,7 +18,7 @@ uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float uint32_t just_contain); void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge, uint32_t is_circle, uint64_t* rLen); - +void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen); void enable_debug_mode(uint32_t mode); #endif \ No newline at end of file