From 676385cf8e422c3d2b198abd4c1db7cea5aa6635 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Wed, 27 Nov 2024 13:27:13 -0500 Subject: [PATCH] ont simplex support --- CommandLines.cpp | 49 ++++- CommandLines.h | 4 +- Correct.cpp | 110 +++++++++- Correct.h | 4 +- Hash_Table.h | 5 +- Overlaps.cpp | 67 ++++++- Overlaps.h | 1 + ecovlp.cpp | 513 +++++++++++++++++++++++++++++++++++++++++++++-- ecovlp.h | 50 +---- gfa_ut.cpp | 181 ++++++++++++++++- 10 files changed, 896 insertions(+), 88 deletions(-) diff --git a/CommandLines.cpp b/CommandLines.cpp index a53f865..25cbd8c 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -7,6 +7,9 @@ #include #include "CommandLines.h" #include "ketopt.h" +#include "kseq.h" + +KSEQ_INIT(gzFile, gzread) #define DEFAULT_OUTPUT "hifiasm.asm" @@ -72,7 +75,9 @@ static ko_longopt_t long_options[] = { { "telo-s", ko_required_argument, 357}, { "ctg-n", ko_required_argument, 358}, { "ont", ko_no_argument, 359}, - { "sc-n", ko_no_argument, 360}, + // { "sc-n", ko_no_argument, 360}, + { "chem-c", ko_required_argument, 361}, + { "chem-f", ko_required_argument, 362}, // { "path-round", ko_required_argument, 348}, { 0, 0, 0 } }; @@ -210,6 +215,14 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " --telo-s INT\n"); fprintf(stderr, " min score for telomere reads [%ld]\n", asm_opt->telo_mic_sc); + fprintf(stderr, " ONT simplex assembly (beta):\n"); + fprintf(stderr, " --ont assemble ONT simplex reads in fastq format\n"); + // fprintf(stderr, " --sc-n consider base qual value for assembly\n"); + fprintf(stderr, " --chem-c INT\n"); + fprintf(stderr, " detect chemical reads with <=INT other reads support [%lu]\n", asm_opt->chemical_cov); + fprintf(stderr, " --chem-f INT\n"); + fprintf(stderr, " length of flanking regions for chemical read detection [%lu]\n", asm_opt->chemical_flank); + fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n"); fprintf(stderr, "See `https://hifiasm.readthedocs.io/en/latest/' or `man ./hifiasm.1' for complete documentation.\n"); @@ -342,6 +355,8 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->is_ont = 0; asm_opt->is_sc = 0; + asm_opt->chemical_cov = 1; + asm_opt->chemical_flank = 256; } void destory_enzyme(enzyme* f) @@ -692,25 +707,33 @@ int check_option(hifiasm_opt_t* asm_opt) void get_queries(int argc, char *argv[], ketopt_t* opt, hifiasm_opt_t* asm_opt) { - if(opt->ind == argc) - { + if(opt->ind == argc) { return; } asm_opt->num_reads = argc - opt->ind; asm_opt->read_file_names = (char**)malloc(sizeof(char*)*asm_opt->num_reads); - long long i; - gzFile dfp; - for (i = 0; i < asm_opt->num_reads; i++) - { + long long i; int ret; + gzFile dfp; kseq_t *ks = NULL; + for (i = 0; i < asm_opt->num_reads; i++) { asm_opt->read_file_names[i] = argv[i + opt->ind]; dfp = gzopen(asm_opt->read_file_names[i], "r"); - if (dfp == 0) - { + if (dfp == 0) { fprintf(stderr, "[ERROR] Cannot find the input read file: %s\n", asm_opt->read_file_names[i]); exit(0); + } else if(asm_opt->is_sc){ + ks = kseq_init(dfp); + while (((ret = kseq_read(ks)) >= 0)) { + if((ks->qual.l == 0) || (ks->qual.s == NULL)) { + fprintf(stderr, "[ERROR] %s is in fasta format rather than fastq format\n", asm_opt->read_file_names[i]); + asm_opt->is_sc = 0; + exit(0); + } + break; + } + kseq_destroy(ks); ks = NULL; } gzclose(dfp); } @@ -914,9 +937,13 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 357) asm_opt->telo_mic_sc = atol(opt.arg); else if (c == 358) asm_opt->max_contig_tip = atol(opt.arg); else if (c == 359) { - asm_opt->is_ont = 1; asm_opt->max_ov_diff_ec = 0.07; ///asm_opt->mz_win = 37; asm_opt->k_mer_length = 37; - } else if (c == 360) { + asm_opt->is_ont = 1; asm_opt->max_ov_diff_ec = 0.07; asm_opt->is_sc = 1; ///asm_opt->mz_win = 37; asm_opt->k_mer_length = 37; + } /**else if (c == 360) { asm_opt->is_sc = 1; + }**/ else if (c == 361) { + asm_opt->chemical_cov = atol(opt.arg); + } else if (c == 362) { + asm_opt->chemical_flank = atol(opt.arg); } else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg); } diff --git a/CommandLines.h b/CommandLines.h index bd483d6..e6505bf 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.21.0-r666" +#define HA_VERSION "0.21.0-r686" #define VERBOSE 0 @@ -163,6 +163,8 @@ typedef struct { uint64_t is_ont; uint64_t is_sc; + uint64_t chemical_cov; + uint64_t chemical_flank; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Correct.cpp b/Correct.cpp index fa0c3e0..c71036b 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -24531,7 +24531,108 @@ void reassign_gaps(overlap_region *z, overlap_region *aux, char* qstr, int64_t q // } } -void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf) +uint32_t is_ovlp_debug(int64_t s, int64_t e, int64_t ws, int64_t we, int64_t op) +{ + int64_t os, oe, ovlp; + os = MAX(s, ws); oe = MIN(e, we); + ovlp = ((oe>os)? (oe-os):0); + if(op != 2) { + return (!!ovlp); + } else { + if(ws >= s && we <= e) return 1; + } + return 0; +} + +uint32_t inline ff_tend(overlap_region *z, int64_t wn, int64_t dn, double dr, double er, int64_t min_err) +{ + int64_t k, zwn = z->w_list.n, err, mm, qi, ci, cn, ql, ws, we, zs, ze, s, e, os, oe, ovlp; bit_extz_t ez; uint32_t cl; uint16_t c; + zs = z->x_pos_s; ze = z->x_pos_e + 1; ql = ze - zs; + if(ql < wn) return 0; if(dn > (ql*dr)) dn = ql*dr; if(dn < wn) return 0; if(ql < dn) return 0; + + s = zs; e = zs + dn; + // if(z->y_id == 27) fprintf(stderr, "-a-[M::%s] tid::%u\t%.*s\tz::[%ld,%ld)\ti::[%ld,%ld)\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), zs, ze, s, e); + for (k = err = mm = 0, qi = zs; (k < zwn) && (z->w_list.a[k].x_start < e); k++) { + // if(z->y_id == 27) fprintf(stderr, "-a-[M::%s] tid::%u\t%.*s\tw::[%d,%d)\terr::%d\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), z->w_list.a[k].x_start, z->w_list.a[k].x_end + 1, z->w_list.a[k].error); + if(!(is_ualn_win(z->w_list.a[k]))) { + set_bit_extz_t(ez, (*z), k); + ci = 0; cn = ez.cigar.n; qi = ez.ts; //ti = ez.ps; + while (ci < cn && qi < e) { + ws = qi; + ci = pop_trace(&(ez.cigar), ci, &c, &cl); + if(c!=2) qi += cl; + // if(c!=3) ti += cl; + we = qi; + + if(c == 0) { + os = MAX(s, ws); oe = MIN(e, we); + ovlp = ((oe>os)? (oe-os):0); mm += ovlp; + } else { + err += cl; + } + // assert(is_ovlp_debug(s, e, ws, we, c)); + + if((err > min_err) && ((mm + err) > wn) && (err > ((mm + err)*er))) { + // fprintf(stderr, "-0-[M::%s] tid::%u\t%.*s\tmm::%ld\terr::%ld\tdn::%ld\twn::%ld\tdif::%f\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), mm, err, dn, wn, er); + return 1; + } + } + } else { + ws = z->w_list.a[k].x_start; we = z->w_list.a[k].x_end + 1; + err += we - ws; + + // assert(is_ovlp_debug(s, e, ws, we, -1)); + + if((err > min_err) && ((mm + err) > wn) && (err > ((mm + err)*er))) { + // fprintf(stderr, "-1-[M::%s] tid::%u\t%.*s\tmm::%ld\terr::%ld\tdn::%ld\twn::%ld\tdif::%f\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), mm, err, dn, wn, er); + return 1; + } + } + } + + s = ze - dn; e = ze; + for (k = zwn - 1, err = mm = 0, qi = ze; (k >= 0) && ((z->w_list.a[k].x_end + 1) > s); k--) { + if(!(is_ualn_win(z->w_list.a[k]))) { + set_bit_extz_t(ez, (*z), k); + ci = ((int64_t)ez.cigar.n) - 1; qi = ez.te + 1; //ti = ez.pe + 1; + while (ci >= 0 && qi > s) { + we = qi; + ci = pop_trace_back(&(ez.cigar), ci, &c, &cl); + if(c!=2) qi -= cl; + // if(c!=3) ti += cl; + ws = qi; + + if(c == 0) { + os = MAX(s, ws); oe = MIN(e, we); + ovlp = ((oe>os)? (oe-os):0); mm += ovlp; + } else { + err += cl; + } + + // assert(is_ovlp_debug(s, e, ws, we, c)); + + if((err > min_err) && ((mm + err) > wn) && (err > ((mm + err)*er))) { + // fprintf(stderr, "-2-[M::%s] tid::%u\t%.*s\tmm::%ld\terr::%ld\tdn::%ld\twn::%ld\tdif::%f\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), mm, err, dn, wn, er); + return 1; + } + } + } else { + ws = z->w_list.a[k].x_start; we = z->w_list.a[k].x_end + 1; + err += we - ws; + + // assert(is_ovlp_debug(s, e, ws, we, -1)); + + if((err > min_err) && ((mm + err) > wn) && (err > ((mm + err)*er))) { + // fprintf(stderr, "-3-[M::%s] tid::%u\t%.*s\tmm::%ld\terr::%ld\tdn::%ld\twn::%ld\tdif::%f\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), mm, err, dn, wn, er); + return 1; + } + } + } + + return 0; +} + +void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop) { uint64_t i, bs, k, ql = qu->length; Window_Pool w; double err, e_max, rr; int64_t re; overlap_region t; overlap_region *z; //asg64_v iidx, buf, buf1; @@ -24565,6 +24666,9 @@ void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rre if(!gen_hc_fast_cigar(z, cl, rref, w.window_length, qu->seq, tu, exz, aux_o, e_rate, ql, rid, khit, &re)) continue; + if(chem_drop && ff_tend(z, 384, 2000, 0.1, (((e_rate*10)<0.36)?(e_rate*10):(0.36)), 128)) continue; + + // if(z->x_id == 3196 && z->y_id == 3199) fprintf(stderr, "-1-[M::%s] tid::%u\t%.*s\trr::%f\tre::%ld\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), rr, re); reassign_gaps(z, aux_o, qu->seq, ql, NULL, -1, rref, tu, buf); @@ -24586,7 +24690,7 @@ void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rre } -void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf) +void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop) { uint64_t i, bs, k, ql = qu->length; Window_Pool w; double err, e_max, rr; int64_t re; overlap_region t; overlap_region *z; //asg64_v iidx, buf, buf1; @@ -24628,6 +24732,8 @@ void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads if(!gen_hc_fast_cigar(z, cl, rref, w.window_length, qu->seq, tu, exz, aux_o, e_rate, ql, rid, khit, &re)) continue; + if(chem_drop && ff_tend(z, 384, 2000, 0.1, (((e_rate*10)<0.36)?(e_rate*10):(0.36)), 128)) continue; + // if(z->x_id == 3196 && z->y_id == 3199) fprintf(stderr, "-1-[M::%s] tid::%u\t%.*s\trr::%f\tre::%ld\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), rr, re); ///debug for memory // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); diff --git a/Correct.h b/Correct.h index 33ceacf..eef942a 100644 --- a/Correct.h +++ b/Correct.h @@ -1391,8 +1391,8 @@ int64_t get_rid_backward_cigar_err(rtrace_iter *it, ul_ov_t *aln, kv_rtrace_t *t const ul_idx_t *uref, char* qstr, UC_Read *tu, overlap_region_alloc *ol, overlap_region *o, bit_extz_t *exz, double e_rate, int64_t qs); -void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf); -void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf); +void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop); +void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop); uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf); void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, Chain_Data *dp, asg8_v *q8, asg8_v *t8); void set_exact_exz(bit_extz_t *exz, int64_t qs, int64_t qe, int64_t ts, int64_t te); diff --git a/Hash_Table.h b/Hash_Table.h index 8fc79ac..b7ad647 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -9,8 +9,9 @@ #define WINDOW 375 #define WINDOW_BOUNDARY 375 #define WINDOW_HC 775 -#define WINDOW_OHC 475 -// #define WINDOW_OHC 375 ///ONT high error +///ONT high error +// #define WINDOW_OHC 475 +#define WINDOW_OHC 375 #define WINDOW_HC_FAST 512 ///for one side, the first or last WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY bases should not be corrected #define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25 diff --git a/Overlaps.cpp b/Overlaps.cpp index fb5a169..d4fd9ba 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -19,6 +19,7 @@ #include "gfa_ut.h" #include "assert.h" #include "khash.h" +#include "ecovlp.h" KHASH_SET_INIT_INT64(64) @@ -63,7 +64,6 @@ KRADIX_SORT_INIT(ha_mzl_t_srt1, ha_mzl_t, ha_mzl_t_key, member_size(ha_mzl_t, x) #define Uc_beg(z) ((uint32_t)((z).a[0]>>32)) #define Uc_end(z) ((uint32_t)((z).a[(z).n-1]>>32)^1) -#define UL_COV_THRES 2 #define PHASE_SEP 64 #define PHASE_SEF 2 #define PHASE_SEP_RATE 0.04 @@ -2969,10 +2969,18 @@ int max_hang, int min_ovlp) void prt_specific_overlap(ma_hit_t_alloc *src, uint64_t qn, uint64_t tn, const char *cmd) { - int64_t idx = get_specific_overlap(&(src[qn]), qn, tn); - const ma_hit_t *h = &(src[qn].buffer[idx]); - fprintf(stderr, "%s::idx::%ld[M::%s::] qn::%u, tn::%u, del::%u, bl::%u, ml::%u\n", cmd, idx, __func__, - Get_qn(*h), Get_tn(*h), h->del, h->bl, h->ml); + ma_hit_t *h = NULL; int64_t idx = -1, k; + if(tn != ((uint64_t)-1)) { + idx = get_specific_overlap(&(src[qn]), qn, tn); h = &(src[qn].buffer[idx]); + fprintf(stderr, "%s::idx::%ld[M::%s::] qn::%u, tn::%u, del::%u, bl::%u, ml::%u\n", cmd, idx, __func__, + Get_qn(*h), Get_tn(*h), h->del, h->bl, h->ml); + } else { + for (k = 0; k < src[qn].length; k++) { + h = &(src[qn].buffer[k]); + fprintf(stderr, "%s::idx::%ld[M::%s::] qn::%u, tn::%u, del::%u, bl::%u, ml::%u\n", cmd, idx, __func__, + Get_qn(*h), Get_tn(*h), h->del, h->bl, h->ml); + } + } } asg_t *ma_sg_gen_ul(ma_hit_t_alloc* sources, int64_t n_read, const ma_sub_t *coverage_cut, @@ -5667,8 +5675,12 @@ long long weakID, uint32_t w_qs, uint32_t w_qe) { strongID = Get_tn(aim_paf->buffer[i]); index = get_specific_overlap(&(reverse_paf_list[strongID]), strongID, weakID); - if(index != -1) - { + if(index != -1) { + // if((Get_qn(aim_paf->buffer[i]) == 27087 && weakID == 27128) || (Get_qn(aim_paf->buffer[i]) == 27128 && weakID == 27087)) { + // ma_hit_t *h = &(aim_paf->buffer[i]); + // fprintf(stderr, "[M::%s]\t%.*s(qn::%u)\t%u\t%u\t%u\t%c\t%.*s(tn::%u)\t%u\t%u\t%u\t%u\t%u\t255\n", __func__, (int)Get_NAME_LENGTH(R_INF, Get_qn(*h)), Get_NAME((R_INF), Get_qn(*h)), Get_qn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_qn(*h)), Get_qs(*h), Get_qe(*h), "+-"[h->rev], + // (int)Get_NAME_LENGTH(R_INF, Get_tn(*h)), Get_NAME((R_INF), Get_tn(*h)), Get_tn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_tn(*h)), Get_ts(*h), Get_te(*h), h->ml, h->bl); + // } return 0; } } @@ -11193,6 +11205,11 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source index = get_specific_overlap(&(sources[tn]), tn, qn); // if(index < 0 || index >= sources[tn].length) fprintf(stderr, "sb, tn: %u, qn: %u, index: %ld, length: %u\n", tn, qn, index, sources[tn].length); sources[tn].buffer[index].bl |= ((uint32_t)0x40000000); + // ma_hit_t *h = &(sources[i].buffer[j]); + // if((Get_qn(*h) == 27087 && Get_tn(*h) == 27128) || (Get_tn(*h) == 27087 && Get_qn(*h) == 27128)) { + // fprintf(stderr, "[M::%s]\t%.*s(qn::%u)\t%u\t%u\t%u\t%c\t%.*s(tn::%u)\t%u\t%u\t%u\t%u\t%u\t255\n", __func__, (int)Get_NAME_LENGTH(R_INF, Get_qn(*h)), Get_NAME((R_INF), Get_qn(*h)), Get_qn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_qn(*h)), Get_qs(*h), Get_qe(*h), "+-"[h->rev], + // (int)Get_NAME_LENGTH(R_INF, Get_tn(*h)), Get_NAME((R_INF), Get_tn(*h)), Get_tn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_tn(*h)), Get_ts(*h), Get_te(*h), h->ml, h->bl); + // } } } } @@ -38528,17 +38545,43 @@ void prt_dbg_gfa(asg_t *sg, const char *suffix, ma_sub_t *cov, ma_hit_t_alloc* s free(gfa_name); free(o_file); } +void prt_dbg_rid_ovlp(ma_hit_t_alloc *ov, int64_t rid, char *rn, const char *cmd) +{ + uint64_t k; ma_hit_t *h = NULL; + if(rid < 0) { + for (k = 0; k < R_INF.total_reads; k++) { + if (memcmp(rn, Get_NAME((R_INF), k), Get_NAME_LENGTH((R_INF), k)) == 0) break; + } + if(k >= R_INF.total_reads) return; + rid = k; + } + + fprintf(stderr, "\n[M::%s::%s::id::%ld]\n", __func__, cmd, rid); + for (k = 0; k < ov[rid].length; k++) { + h = &(ov[rid].buffer[k]); + if(h->del) continue; + fprintf(stderr, "%.*s(qn::%u)\t%u\t%u\t%u\t%c\t%.*s(tn::%u)\t%u\t%u\t%u\t%u\t%u\t255\n", (int)Get_NAME_LENGTH(R_INF, Get_qn(*h)), Get_NAME((R_INF), Get_qn(*h)), Get_qn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_qn(*h)), Get_qs(*h), Get_qe(*h), "+-"[h->rev], + (int)Get_NAME_LENGTH(R_INF, Get_tn(*h)), Get_NAME((R_INF), Get_tn(*h)), Get_tn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_tn(*h)), Get_ts(*h), Get_te(*h), h->ml, h->bl); + } +} + asg_t *gen_init_sg(int32_t min_dp, uint64_t n_read, int64_t mini_overlap_length, int64_t max_hang_length, int64_t gap_fuzz, ma_hit_t_alloc* src, uint64_t* readLen, R_to_U* ruIndex, bub_label_t *b_mask_t, ma_sub_t** cov, all_ul_t *ul, telo_end_t *te) { asg_t *sg = NULL; - // prt_specific_overlap(src, 22233, 22235, "1"); + // prt_dbg_rid_ovlp(src, 3700, NULL, "1"); if(ul) rescue_src_ul(src, n_read, UL_COV_THRES); + // prt_dbg_rid_ovlp(src, 3700, NULL, "2"); ma_hit_sub(min_dp, src, n_read, readLen, mini_overlap_length, cov); + // prt_dbg_rid_ovlp(src, 3700, NULL, "3"); detect_chimeric_reads(src, n_read, readLen, *cov, asm_opt.max_ov_diff_final*2.0, ul, UL_COV_THRES); + // prt_dbg_rid_ovlp(src, 3700, NULL, "4"); ma_hit_cut(src, n_read, readLen, mini_overlap_length, cov); + // prt_dbg_rid_ovlp(src, 3700, NULL, "5"); ma_hit_flt(src, n_read, *cov, max_hang_length, mini_overlap_length); + // prt_dbg_rid_ovlp(src, 3700, NULL, "6"); ma_hit_contained_advance(src, n_read, *cov, ruIndex, max_hang_length, mini_overlap_length); + // prt_dbg_rid_ovlp(src, 3700, NULL, "7"); if(!ul) { sg = ma_sg_gen(src, n_read, *cov, max_hang_length, mini_overlap_length); @@ -38650,6 +38693,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g) } ///just for debug renew_graph_init(sources, reverse_sources, sg, coverage_cut, ruIndex, n_read); + // if(asm_opt.is_ont) handle_chemical_arc(asm_opt.thread_num, R_INF.total_reads); + // if(asm_opt.is_ont) handle_chemical_r(asm_opt.thread_num, R_INF.total_reads); ///it's hard to say which function is better ///normalize_ma_hit_t_single_side(sources, n_read); @@ -38681,7 +38726,9 @@ ma_sub_t **coverage_cut_ptr, int debug_g) } // prt_specific_overlap(sources, 22233, 22235, "0-b"); // prt_specific_overlap(sources, 22235, 22233, "0-b"); - clean_weak_ma_hit_t(sources, reverse_sources, n_read, asm_opt.ar?UL_COV_THRES:(uint32_t)-1); + // prt_dbg_rid_ovlp(sources, 27087, NULL, "0-a"); + if(!(asm_opt.is_ont)) clean_weak_ma_hit_t(sources, reverse_sources, n_read, asm_opt.ar?UL_COV_THRES:(uint32_t)-1); + // prt_dbg_rid_ovlp(sources, 27087, NULL, "0-b"); // prt_specific_overlap(sources, 22233, 22235, "0-c"); // prt_specific_overlap(sources, 22235, 22233, "0-c"); sg = gen_init_sg(min_dp, n_read, mini_overlap_length, max_hang_length, gap_fuzz, sources, readLen, ruIndex, @@ -39006,7 +39053,9 @@ long long bubble_dist, int read_graph, int write) if (!(asm_opt.flag & HA_F_BAN_ASSEMBLY)) { + if(asm_opt.is_ont) handle_chemical_r(asm_opt.thread_num, R_INF.total_reads); try_rescue_overlaps(sources, reverse_sources, n_read, 4); + clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length, max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio, output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 0); diff --git a/Overlaps.h b/Overlaps.h index 24699f2..98b5153 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -35,6 +35,7 @@ // #define ALTER_LABLE 2 // #define HAP_LABLE 4 #define ug_ext_len 75000 +#define UL_COV_THRES 2 #define Get_qn(RECORD) ((uint32_t)((RECORD).qns>>32)) #define Get_qs(RECORD) ((uint32_t)((RECORD).qns)) diff --git a/ecovlp.cpp b/ecovlp.cpp index a180263..91d7c7c 100644 --- a/ecovlp.cpp +++ b/ecovlp.cpp @@ -14,6 +14,56 @@ #define del_cns_nn(z, nn_i) ((z).a[(nn_i)].sc == CNS_DEL_V) #define REFRESH_N 128 +KDQ_INIT(uint32_t) + +typedef struct { + uint32_t v:31, f:1; + uint32_t sc; +} cns_arc; +typedef struct {size_t n, m, nou; cns_arc *a; } cns_arc_v; + +typedef struct { + // uint16_t c:2, t:2, f:1, sc:3; + uint32_t c:2, f:1, sc:29; + cns_arc_v arc; +}cns_t; + +typedef struct { + size_t n, m; + cns_t *a; + uint32_t si, ei, off, bn, bb0, bb1, cns_g_wl; + kdq_t(uint32_t) *q; +}cns_gfa; + +typedef struct { + // chaining and overlapping related buffers + UC_Read self_read, ovlp_read; + Candidates_list clist; + overlap_region_alloc olist; + ha_abuf_t *ab; + // int64_t num_read_base, num_correct_base, num_recorrect_base; + uint64_t cnt[6], rr; + haplotype_evdience_alloc hap; + bit_extz_t exz; + kv_ul_ov_t pidx; + asg64_v v64; + asg32_v v32; + asg16_v v16; + asg8_v v8q, v8t; + + kvec_t_u8_warp k_flag; + st_mt_t sp; + cns_gfa cns; +} ec_ovec_buf_t0; + +typedef struct { + ec_ovec_buf_t0 *a; + uint32_t n, rev; +} ec_ovec_buf_t; + +ec_ovec_buf_t* gen_ec_ovec_buf_t(uint32_t n); +void destroy_ec_ovec_buf_t(ec_ovec_buf_t *p); + #define generic_key(x) (x) KRADIX_SORT_INIT(ec16, uint16_t, generic_key, 2) @@ -2608,6 +2658,7 @@ void push_ff_ovlp(ma_hit_t_alloc* paf, overlap_region_alloc* ov, uint32_t flag, if(z->el == 1) cnt[4]++; if(z->no_l_indel) cnt[5]++; } + z->del = 0; } } @@ -2729,7 +2780,7 @@ inline uint64_t exact_ec_check(char *qstr, uint64_t ql, char *tstr, uint64_t tl, return 0; } -void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in) +void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in, uint8_t chem_drop) { if(ol->length <= 0) return; @@ -2744,7 +2795,7 @@ void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads * } if(!(srt->n)) { - gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf); + gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop); } else { ///debug for memory // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); @@ -2778,7 +2829,7 @@ void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads * ///debug for memory // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); - if(on > nec) gen_hc_r_alin_nec(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf); + if(on > nec) gen_hc_r_alin_nec(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop); ///debug for memory // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); @@ -2986,6 +3037,158 @@ void debug_retrive_bqual(asg8_v *vq, asg8_v *vt, uint64_t id, uint64_t rn) } } +uint32_t is_uncorrected_read(overlap_region_alloc* ov, asg64_v *idx, int64_t len, int64_t min_len) +{ + uint64_t k, s, e; int64_t dp, old_dp, st = 0, ed; + for (k = idx->n = 0; k < ov->length; k++) { + s = ov->list[k].x_pos_s; e = ov->list[k].x_pos_e + 1; + kv_push(uint64_t, (*idx), (s<<1)); + kv_push(uint64_t, (*idx), (e<<1)|1); + } + + radix_sort_ec64(idx->a, idx->a + idx->n); + for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) { + old_dp = dp; + ///if a[j] is qe + if (idx->a[k]&1) --dp; + else ++dp; + + ed = idx->a[k]>>1; + if(ed > st) { + if(old_dp == 0) { + if((ed - st) >= min_len) return 1; + } + } + st = ed; + } + + + ed = len; old_dp = dp; + if(ed > st) { + if(old_dp == 0) { + if((ed - st) >= min_len) return 1; + if((ed - st) >= len) return 1; + } + } + + return 0; +} + +uint32_t is_chemical_r_qual(overlap_region_alloc *ov, asg64_v *idx, int64_t len, int64_t cov, int64_t frank_len, asg8_v *qv, uint64_t rid) +{ + uint64_t k, s, e; int64_t dp, old_dp, st = 0, ed, s0, s1, e0, e1, rr, qk; + if((frank_len) > (len *0.01)) frank_len = len *0.01; + for (k = idx->n = 0; k < ov->length; k++) { + s = ov->list[k].x_pos_s; e = ov->list[k].x_pos_e + 1; + kv_push(uint64_t, (*idx), (s<<1)); + kv_push(uint64_t, (*idx), (e<<1)|1); + // fprintf(stderr, "[M::%s]\trid::%lu\ts::%lu\te::%lu\n", __func__, rid, s, e); + } + + radix_sort_ec64(idx->a, idx->a + idx->n); s0 = s1 = e0 = e1 = rr = -1; + for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) { + old_dp = dp; + ///if a[j] is qe + if (idx->a[k]&1) --dp; + else ++dp; + + ed = idx->a[k]>>1; + if(ed > st) { + if(old_dp <= cov) { + rr = 1; + } else { + if(s0 < 0) { + s0 = st; s1 = ed; + } + e0 = st; e1 = ed; + } + } + st = ed; + } + + + ed = len; old_dp = dp; + if(ed > st) { + if(old_dp <= cov) { + rr = 1; + } else { + if(s0 < 0) { + s0 = st; s1 = ed; + } + e0 = st; e1 = ed; + } + } + + + if((s0 != e0) && (s1 != e1) && (s0 <= frank_len) && ((len - e1) <= frank_len) && (rr > 0)) { + for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) { + old_dp = dp; + ///if a[j] is qe + if (idx->a[k]&1) --dp; + else ++dp; + + ed = idx->a[k]>>1; + if(ed > st) { + if((old_dp <= cov) && (st >= s0) && (ed <= e1)) { + retrive_bqual(qv, NULL, rid, -1, -1, 0, sc_bn); + fprintf(stderr, "[M::%s]\tlf::[%ld,%ld)\trt::[%ld,%ld)\tmd::[%ld,%ld)\tcov::%ld\n", __func__, s0, s1, e0, e1, st, ed, old_dp); + for (qk = st; qk < ed; qk++) fprintf(stderr, "%u", qv->a[qk]); + fprintf(stderr, "\n"); + return 1; + } + } + st = ed; + } + + + ed = len; old_dp = dp; + if(ed > st) { + if((old_dp <= cov) && (st >= s0) && (ed <= e1)) { + retrive_bqual(qv, NULL, rid, -1, -1, 0, sc_bn); + fprintf(stderr, "[M::%s]\tlf::[%ld,%ld)\trt::[%ld,%ld)\tmd::[%ld,%ld)\tcov::%ld\n", __func__, s0, s1, e0, e1, st, ed, old_dp); + for (qk = st; qk < ed; qk++) fprintf(stderr, "%u", qv->a[qk]); + fprintf(stderr, "\n"); + return 1; + } + } + } + + + // for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) { + // old_dp = dp; + // ///if a[j] is qe + // if (idx->a[k]&1) --dp; + // else ++dp; + + // ed = idx->a[k]>>1; + // if(ed > st) { + // if((old_dp <= cov)) { + // retrive_bqual(qv, NULL, rid, -1, -1, 0, sc_bn); + // fprintf(stderr, "[M::%s]\tlf::[%ld,%ld)\trt::[%ld,%ld)\tmd::[%ld,%ld)\tcov::%ld\n", __func__, s0, s1, e0, e1, st, ed, old_dp); + // for (qk = st; qk < ed; qk++) fprintf(stderr, "%u", qv->a[qk]); + // fprintf(stderr, "\n"); + // return 1; + // } + // } + // st = ed; + // } + + + // ed = len; old_dp = dp; + // if(ed > st) { + // if((old_dp <= cov)) { + // retrive_bqual(qv, NULL, rid, -1, -1, 0, sc_bn); + // fprintf(stderr, "[M::%s]\tlf::[%ld,%ld)\trt::[%ld,%ld)\tmd::[%ld,%ld)\tcov::%ld\n", __func__, s0, s1, e0, e1, st, ed, old_dp); + // for (qk = st; qk < ed; qk++) fprintf(stderr, "%u", qv->a[qk]); + // fprintf(stderr, "\n"); + // return 1; + // } + // } + + return 0; +} + + static void worker_hap_ec(void *data, long i, int tid) { ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); @@ -3006,7 +3209,7 @@ static void worker_hap_ec(void *data, long i, int tid) // if(i % 100000 == 0) fprintf(stderr, "-a-[M::%s-beg] rid->%ld\n", __func__, i); // if (memcmp("c42804f3-0e13-43a0-8a71-b91b40accf9a", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { // if (memcmp("b2e68ecf-381a-439c-b676-c1e6831d6acf", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { - // if (memcmp("0aec8c4f-c849-4c31-85ba-4ffb297eeb28", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { + // if (memcmp("64b2c27d-86b8-451e-9330-6ba62be2ffcc", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { // fprintf(stderr, "-a-[M::%s-beg] rid->%ld\n", __func__, i); // } else { // return; @@ -3023,7 +3226,7 @@ static void worker_hap_ec(void *data, long i, int tid) recover_UC_Read(&b->self_read, &R_INF, i); qlen = b->self_read.length; - h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, /**((asm_opt.is_ont)?(0.05):(0.02))**/0.02, asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32);///ONT high error + h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, ((asm_opt.is_ont)?(0.05):(0.02)), asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32);///ONT high error // b->num_read_base += b->olist.length; b->cnt[0] += b->self_read.length; @@ -3033,7 +3236,7 @@ static void worker_hap_ec(void *data, long i, int tid) ///debug for memory // snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n); - gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT/**asm_opt.k_mer_length**/, 1, &b->v16, &b->v64, &(R_INF.paf[i])); + gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT/**asm_opt.k_mer_length**/, 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont); // prt_ovlp_sam(&b->olist, &b->ovlp_read, b->self_read.seq, b->self_read.length); @@ -3061,7 +3264,10 @@ static void worker_hap_ec(void *data, long i, int tid) push_nec_re(aux_o, &(scc.a[i])); push_nec_re(aux_o, &(scb.a[i])); - + // if((asm_opt.is_ont) && is_chemical_r_qual(&b->olist, &b->v64, qlen, 1, 16, &(b->v8q), i)/**(is_uncorrected_read(&b->olist, &b->v64, qlen, 1600))**/) { + // // b->olist.length = 0; + // fprintf(stderr, "[M::%s] rid::%ld\t%.*s\n\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + // } push_ne_ovlp(&(R_INF.paf[i]), &b->olist, 1, &R_INF, &(scc.a[i])/**, i, &b->self_read, &b->ovlp_read**/); push_ne_ovlp(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, NULL/**, i, NULL, NULL**/); @@ -3605,9 +3811,9 @@ static void worker_hap_dc_ec_gen_new_idx(void *data, long i, int tid) ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO); - uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; + uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; uint32_t qlen = 0; - recover_UC_Read(&b->self_read, &R_INF, i); + recover_UC_Read(&b->self_read, &R_INF, i); qlen = b->self_read.length; h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, /**0.02**/0.001, asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32); @@ -3624,6 +3830,11 @@ static void worker_hap_dc_ec_gen_new_idx(void *data, long i, int tid) h_ec_lchain_fast_new(b->ab, i, &b->self_read, &b->ovlp_read, &R_INF, &b->olist, &b->clist, &b->exz, &b->v16, &b->v64, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), 0.866666); + if((asm_opt.is_ont) && (is_uncorrected_read(&b->olist, &b->v64, qlen, 1600))) { + b->olist.length = 0; + // fprintf(stderr, "[M::%s] rid::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + } + push_ff_ovlp(&(R_INF.paf[i]), &b->olist, 1, &R_INF, b->cnt); push_ff_ovlp(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, b->cnt); @@ -3641,6 +3852,236 @@ static void worker_hap_dc_ec_gen_new_idx(void *data, long i, int tid) refresh_ec_ovec_buf_t0(b, REFRESH_N); } +uint32_t is_chemical_r(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_t cov, int64_t frank_len) +{ + uint64_t k, s, e; int64_t dp, old_dp, st = 0, ed, s0, s1, e0, e1, rr; + if((frank_len) > (len *0.01)) frank_len = len *0.01; + for (k = idx->n = 0; k < ov->length; k++) { + s = (uint32_t)ov->buffer[k].qns; e = ov->buffer[k].qe; + kv_push(uint64_t, (*idx), (s<<1)); + kv_push(uint64_t, (*idx), (e<<1)|1); + } + + radix_sort_ec64(idx->a, idx->a + idx->n); s0 = s1 = e0 = e1 = rr = -1; + for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) { + old_dp = dp; + ///if a[j] is qe + if (idx->a[k]&1) --dp; + else ++dp; + + ed = idx->a[k]>>1; + if(ed > st) { + if(old_dp <= cov) { + rr = 1; + } else { + if(s0 < 0) { + s0 = st; s1 = ed; + } + e0 = st; e1 = ed; + } + } + st = ed; + } + + + ed = len; old_dp = dp; + if(ed > st) { + if(old_dp <= cov) { + rr = 1; + } else { + if(s0 < 0) { + s0 = st; s1 = ed; + } + e0 = st; e1 = ed; + } + } + + + if((s0 != e0) && (s1 != e1) && (s0 <= frank_len) && ((len - e1) <= frank_len) && (rr > 0)) { + for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) { + old_dp = dp; + ///if a[j] is qe + if (idx->a[k]&1) --dp; + else ++dp; + + ed = idx->a[k]>>1; + if(ed > st) { + if((old_dp <= cov) && (st >= s0) && (ed <= e1)) { + // if((ov->buffer[0].qns>>32) == 3364) fprintf(stderr, "[M::%s]\tlf::[%ld,%ld)\trt::[%ld,%ld)\tmd::[%ld,%ld)\tcov::%ld\n", __func__, s0, s1, e0, e1, st, ed, old_dp); + return 1; + } + } + st = ed; + } + + + ed = len; old_dp = dp; + if(ed > st) { + if((old_dp <= cov) && (st >= s0) && (ed <= e1)) { + // if((ov->buffer[0].qns>>32) == 3364) fprintf(stderr, "[M::%s]\tlf::[%ld,%ld)\trt::[%ld,%ld)\tmd::[%ld,%ld)\tcov::%ld\n", __func__, s0, s1, e0, e1, st, ed, old_dp); + return 1; + } + } + } + + return 0; +} + + +uint32_t is_chemical_r_adv(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_t cov, int64_t cut_len, double dup_rate) +{ + uint64_t k, s, e; int64_t dp, old_dp, st = 0, ed, s0, e0, rr, lt; + for (k = idx->n = 0; k < ov->length; k++) { + s0 = (uint32_t)ov->buffer[k].qns; e0 = ov->buffer[k].qe; + if(s0 > 0) s0 += cut_len; + if(e0 < len) e0 -= cut_len; + if(e0 <= s0) continue; + s = s0; e = e0; + + lt = Get_READ_LENGTH((R_INF), ov->buffer[k].tn); + rr = (lt >= len)?(lt - len):(len - lt); + if((rr <= (len*dup_rate)) && (rr <= (lt*dup_rate)) && (ov->buffer[k].rev)) { + dp = (ov->buffer[k].qe) - ((uint32_t)ov->buffer[k].qns); dp = len - dp; + old_dp = ov->buffer[k].te - ov->buffer[k].ts; old_dp = lt - old_dp; + if((dp <= (len*dup_rate)) && (old_dp <= (lt*dup_rate))) continue; + } + + kv_push(uint64_t, (*idx), (s<<1)); + kv_push(uint64_t, (*idx), (e<<1)|1); + } + + radix_sort_ec64(idx->a, idx->a + idx->n); s0 = e0 = rr = -1; + for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) { + old_dp = dp; + ///if a[j] is qe + if (idx->a[k]&1) --dp; + else ++dp; + + ed = idx->a[k]>>1; + if(ed > st) { + // if(ov->length && ((ov->buffer[0].qns>>32) == 5045637)) { + // fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\tid::%lu\n", __func__, st, ed, old_dp, len, ov->buffer[0].qns>>32); + // } + if(old_dp <= cov) { + // if(ov->length && (ov->buffer[0].qns>>32) == 22344) fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\n", __func__, st, ed, old_dp, len); + return 1; + } + } + st = ed; + } + + + ed = len; old_dp = dp; + if(ed > st) { + // if(ov->length && ((ov->buffer[0].qns>>32) == 5045637)) { + // fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\tid::%lu\n", __func__, st, ed, old_dp, len, ov->buffer[0].qns>>32); + // } + if(old_dp <= cov) { + // if(ov->length && (ov->buffer[0].qns>>32) == 22344) fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\n", __func__, st, ed, old_dp, len); + return 1; + } + } + + return 0; +} + +void prt_dbg_rid_paf(ma_hit_t_alloc *ov, UC_Read *ra, asg8_v *qa) +{ + if(!(ov->length)) return; + uint64_t k, qn = (ov->buffer[0].qns>>32), qn_n, i, m; char *nn = NULL; FILE *fp = NULL; ma_hit_t *h = NULL; + qn_n = Get_NAME_LENGTH((R_INF), qn) + 64; MALLOC(nn, qn_n); + + sprintf(nn, "%.*s.qry.fq", (int)Get_NAME_LENGTH(R_INF, qn), Get_NAME((R_INF), qn)); fp = fopen(nn, "w"); + for (k = 0; k < ov->length; k++) { + i = ov->buffer[k].tn; + recover_UC_Read(ra, &R_INF, i); + fprintf(fp, "@%.*s\n", (int32_t)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + fprintf(fp, "%.*s\n", (int32_t)ra->length, ra->seq); + fprintf(fp, "+\n"); + retrive_bqual(qa, NULL, i, -1, -1, 0, sc_bn); + for (m = 0; m < qa->n; m++) fprintf(fp, "%c", (char)(sc_tb[qa->a[m]] + 33 - 1)); + fprintf(fp, "\n"); + } + fclose(fp); + + sprintf(nn, "%.*s.ref.fq", (int)Get_NAME_LENGTH(R_INF, qn), Get_NAME((R_INF), qn)); fp = fopen(nn, "w"); + i = qn; + recover_UC_Read(ra, &R_INF, i); + fprintf(fp, "@%.*s\n", (int32_t)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + fprintf(fp, "%.*s\n", (int32_t)ra->length, ra->seq); + fprintf(fp, "+\n"); + retrive_bqual(qa, NULL, i, -1, -1, 0, sc_bn); + for (m = 0; m < qa->n; m++) fprintf(fp, "%c", (char)(sc_tb[qa->a[m]] + 33 - 1)); + fprintf(fp, "\n"); + fclose(fp); + + sprintf(nn, "%.*s.ref.fa", (int)Get_NAME_LENGTH(R_INF, qn), Get_NAME((R_INF), qn)); fp = fopen(nn, "w"); + i = qn; + recover_UC_Read(ra, &R_INF, i); + fprintf(fp, ">%.*s\n", (int32_t)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + fprintf(fp, "%.*s\n", (int32_t)ra->length, ra->seq); + // fprintf(fp, "+\n"); + // retrive_bqual(qa, NULL, i, -1, -1, 0, sc_bn); + // for (m = 0; m < qa->n; m++) fprintf(fp, "%c", (char)(sc_tb[qa->a[m]] + 33 - 1)); + // fprintf(fp, "\n"); + fclose(fp); + + sprintf(nn, "%.*s.ov.paf", (int)Get_NAME_LENGTH(R_INF, qn), Get_NAME((R_INF), qn)); fp = fopen(nn, "w"); + for (k = 0; k < ov->length; k++) { + h = &(ov->buffer[k]); + fprintf(fp, "%.*s(qn::%u)\t%u\t%u\t%u\t%c\t%.*s(tn::%u)\t%u\t%u\t%u\t%u\t%u\t255\n", (int)Get_NAME_LENGTH(R_INF, Get_qn(*h)), Get_NAME((R_INF), Get_qn(*h)), Get_qn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_qn(*h)), Get_qs(*h), Get_qe(*h), "+-"[h->rev], + (int)Get_NAME_LENGTH(R_INF, Get_tn(*h)), Get_NAME((R_INF), Get_tn(*h)), Get_tn(*h), (uint32_t)Get_READ_LENGTH(R_INF, Get_tn(*h)), Get_ts(*h), Get_te(*h), h->ml, h->bl); + } + fclose(fp); + + free(nn); +} + +static void worker_hap_dc_ec_chemical_r(void *data, long i, int tid) +{ + ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); + ma_hit_t_alloc *paf = &(R_INF.paf[i]); uint64_t k, m; + + // if (memcmp("3ed80bc4-1169-4948-a9ff-9c2463b7f7a2", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { + // fprintf(stderr, "-a-[M::%s-beg] rid->%ld, b->rr->%lu\n", __func__, i, b->rr); + // } + if(b->cnt[1] == 0) { + // if(i == 6204620) prt_dbg_rid_paf(&(R_INF.paf[i]), &(b->self_read), &(b->v8q)); + // if(is_chemical_r(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), 3, 16)) { + if(is_chemical_r_adv(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), asm_opt.chemical_cov, asm_opt.chemical_flank, 0.02)) { + // fprintf(stderr, "-um-[M::%s]\tqn::%u::%.*s\n\n", __func__, (uint32_t)(i), (int)Get_NAME_LENGTH(R_INF, i), Get_NAME((R_INF), i)); + R_INF.paf[i].length = 0; b->cnt[0]++; + } + } else if(b->cnt[1] == 1) { + for (k = 0; k < paf->length; k++) { + if(R_INF.paf[paf->buffer[k].tn].length == 0) { + paf->buffer[k].tn = (uint32_t)-1; b->cnt[0]++; + } + } + } else { + for (k = m = 0; k < paf->length; k++) { + if(paf->buffer[k].tn == ((uint32_t)-1)) continue; + paf->buffer[m++] = paf->buffer[k]; + } + paf->length = m; + } + + refresh_ec_ovec_buf_t0(b, REFRESH_N); +} + +static void worker_hap_dc_ec_chemical_arc(void *data, long i, int tid) +{ + ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); + ma_hit_t_alloc *paf = &(R_INF.paf[i]); uint64_t k; + + if(is_chemical_r_adv(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), asm_opt.chemical_cov, asm_opt.chemical_flank, 0.02)) { + // fprintf(stderr, "-um-[M::%s]\tqn::%u::%.*s\n\n", __func__, (uint32_t)(i), (int)Get_NAME_LENGTH(R_INF, i), Get_NAME((R_INF), i)); + for (k = 0; k < paf->length; k++) paf->buffer[k].del = 1; b->cnt[0]++; + } + + refresh_ec_ovec_buf_t0(b, REFRESH_N); +} + void gen_ovlst_paf(ma_hit_t_alloc *in_e, ma_hit_t_alloc *in_r, asg64_v *ou) { uint32_t n = 0, k; @@ -3843,7 +4284,7 @@ overlap_region* h_ec_lchain_re1(ha_abuf_t *ab, uint32_t rid, UC_Read *qu, UC_Rea // fprintf(stderr, "-0-[M::%s]\n", __func__); - gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf); rs = qu->seq; + gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf, 0); rs = qu->seq; // fprintf(stderr, "-1-[M::%s]\n", __func__); @@ -4507,7 +4948,7 @@ overlap_region* h_ec_lchain_re3(ha_abuf_t *ab, uint32_t rid, UC_Read *qu, UC_Rea // fprintf(stderr, "-0-[M::%s]\n", __func__); - gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf); rs = qu->seq; + gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf, 0); rs = qu->seq; // fprintf(stderr, "-1-[M::%s]\n", __func__); @@ -5588,7 +6029,55 @@ void sl_ec_r(uint64_t n_thre, uint64_t n_a) kt_for(n_thre, worker_sl_ec, b, n_a);///debug_for_fix for (k = 0; k < n_thre; k++) { - free(b[k].a); destory_UC_Read(&b[k].z); + free(b[k].a); destory_UC_Read(&b[k].z); kv_destroy(b[k].q); } free(b); +} + +void handle_chemical_r(uint64_t n_thre, uint64_t n_a) +{ + ec_ovec_buf_t *b = NULL; uint64_t k, chem_n = 0, dedup = 0; + b = gen_ec_ovec_buf_t(n_thre); + for (k = 0; k < n_thre; ++k) { + b->a[k].cnt[0] = 0; b->a[k].cnt[1] = 0; + } + + kt_for(n_thre, worker_hap_dc_ec_chemical_r, b, n_a); + + for (k = 0; k < n_thre; ++k) { + chem_n += b->a[k].cnt[0]; + b->a[k].cnt[0] = 0; b->a[k].cnt[1] = 1; + } + + kt_for(n_thre, worker_hap_dc_ec_chemical_r, b, n_a); + + for (k = 0; k < n_thre; ++k) { + dedup += b->a[k].cnt[0]; + b->a[k].cnt[1] = 2; + } + + kt_for(n_thre, worker_hap_dc_ec_chemical_r, b, n_a); + + fprintf(stderr, "[M::%s] # chemical reads: %lu, # arcs:: %lu\n", __func__, chem_n, dedup); + + destroy_ec_ovec_buf_t(b); +} + +void handle_chemical_arc(uint64_t n_thre, uint64_t n_a) +{ + ec_ovec_buf_t *b = NULL; uint64_t k, chem_n = 0; + b = gen_ec_ovec_buf_t(n_thre); + for (k = 0; k < n_thre; ++k) { + b->a[k].cnt[0] = 0; + } + + kt_for(n_thre, worker_hap_dc_ec_chemical_arc, b, n_a); + + for (k = 0; k < n_thre; ++k) { + chem_n += b->a[k].cnt[0]; + } + + fprintf(stderr, "[M::%s] # chemical reads: %lu\n", __func__, chem_n); + + destroy_ec_ovec_buf_t(b); } \ No newline at end of file diff --git a/ecovlp.h b/ecovlp.h index 06467e3..34c645d 100644 --- a/ecovlp.h +++ b/ecovlp.h @@ -7,58 +7,12 @@ #include "Process_Read.h" #include "kdq.h" -KDQ_INIT(uint32_t) -typedef struct { - uint32_t v:31, f:1; - uint32_t sc; -} cns_arc; -typedef struct {size_t n, m, nou; cns_arc *a; } cns_arc_v; - -typedef struct { - // uint16_t c:2, t:2, f:1, sc:3; - uint32_t c:2, f:1, sc:29; - cns_arc_v arc; -}cns_t; - -typedef struct { - size_t n, m; - cns_t *a; - uint32_t si, ei, off, bn, bb0, bb1, cns_g_wl; - kdq_t(uint32_t) *q; -}cns_gfa; - -typedef struct { - // chaining and overlapping related buffers - UC_Read self_read, ovlp_read; - Candidates_list clist; - overlap_region_alloc olist; - ha_abuf_t *ab; - // int64_t num_read_base, num_correct_base, num_recorrect_base; - uint64_t cnt[6], rr; - haplotype_evdience_alloc hap; - bit_extz_t exz; - kv_ul_ov_t pidx; - asg64_v v64; - asg32_v v32; - asg16_v v16; - asg8_v v8q, v8t; - - kvec_t_u8_warp k_flag; - st_mt_t sp; - cns_gfa cns; -} ec_ovec_buf_t0; - -typedef struct { - ec_ovec_buf_t0 *a; - uint32_t n, rev; -} ec_ovec_buf_t; - -ec_ovec_buf_t* gen_ec_ovec_buf_t(uint32_t n); -void destroy_ec_ovec_buf_t(ec_ovec_buf_t *p); void prt_chain(overlap_region_alloc *o); void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, uint64_t is_sv, uint64_t *tot_b, uint64_t *tot_e); void sl_ec_r(uint64_t n_thre, uint64_t n_a); void cal_ov_r(uint64_t n_thre, uint64_t n_a, uint64_t new_idx); +void handle_chemical_r(uint64_t n_thre, uint64_t n_a); +void handle_chemical_arc(uint64_t n_thre, uint64_t n_a); #endif \ No newline at end of file diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 56aefae..9702f31 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -1407,6 +1407,173 @@ uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, ma_hit_t_alloc *rev, R_to_ if (cnt > 0) asg_cleanup(g); } +uint32_t is_dedup_weak_arc(asg_arc_t *av, uint32_t an, uint32_t ak, ma_hit_t_alloc *rev) +{ + uint32_t k, w, m; ma_hit_t_alloc *z = NULL; + for (k = 0; k < an; k++) { + if((k == ak) || (av[k].del) || (!av[k].strong)) continue; + if(av[k].ol > av[ak].ol) { + z = &(rev[av[k].v>>1]); w = (av[ak].v^(av[k].v&1)); + for (m = 0; m < z->length; m++) { + if((z->buffer[m].tn<<1|(z->buffer[m].rev)) == w) break; + } + if(m < z->length) return 1; + } + } + + return 0; +} + +void asg_arc_cut_weak(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio, +uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, uint32_t test_bub, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len) +{ + asg64_v tx = {0,0,0}, *b = NULL; + uint32_t i, k, v, w, wz, n_vtx = g->n_seq<<1, nv, nw, kv, kw, trioF = (uint32_t)-1, ntrioF = (uint32_t)-1, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou, olw[2], m; + asg_arc_t *av, *aw, *ve, *we, *vl_max, *wl_max; ma_hit_t_alloc *z = NULL; + + if(in) b = in; + else b = &tx; + b->n = 0; + + for (v = 0; v < n_vtx; ++v) { + // if((v>>1)==17078) fprintf(stderr, "[M::%s::] v:%u, del:%u, seq_vis:%u\n", __func__, v, g->seq[v>>1].del, g->seq_vis[v]); + if (g->seq[v>>1].del) continue; + if((test_bub == 0) || (g->seq_vis[v] == 0)) { + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + if (nv < 2) continue; + + for (i = kv = 0; i < nv; ++i) { + if(av[i].del) continue; + kv++; + } + if(kv < 2) continue; + + for (i = olw[1] = 0, olw[0] = (uint32_t)-1; i < nv; ++i) { + if(av[i].del) continue; + if(av[i].strong) { + if(av[i].ol > olw[av[i].strong]) olw[av[i].strong] = av[i].ol; + } else { + if(av[i].ol < olw[av[i].strong]) olw[av[i].strong] = av[i].ol; + } + } + if(olw[1] <= olw[0]) continue; + for (i = 0; i < nv; ++i) { + if(av[i].del || av[i].strong) continue; + if(av[i].ol >= olw[1]) continue; + if(max_drop_len && av[i].ol >= (*max_drop_len)) continue; + if(is_dedup_weak_arc(av, nv, i, rev)) { + kv_push(uint64_t, *b, (((uint64_t)av[i].ol)<<32) | ((uint64_t)(av-g->arc+i))); + } + } + } + } + + if(rev && rI) memset(g->seq_vis, 0, g->n_seq*2*sizeof(uint8_t)); + radix_sort_srt64(b->a, b->a + b->n); + for (k = 0; k < b->n; k++) { + if(g->arc[(uint32_t)b->a[k]].del) continue; + + v = g->arc[(uint32_t)b->a[k]].ul>>32; w = g->arc[(uint32_t)b->a[k]].v^1; + if(g->seq[v>>1].del || g->seq[w>>1].del) continue; + nv = asg_arc_n(g, v); nw = asg_arc_n(g, w); + av = asg_arc_a(g, v); aw = asg_arc_a(g, w); + if(nv<=1 && nw <= 1) continue; + + if(is_trio) { + if(get_arcs(g, v, NULL, 0)<=1 && get_arcs(g, w, NULL, 0)<=1) continue;///speedup + trioF = get_tip_trio_infor(g, v^1); + ntrioF = (trioF==FATHER? MOTHER : (trioF==MOTHER? FATHER : (uint32_t)-1)); + } + + ve = &(g->arc[(uint32_t)b->a[k]]); + for (i = 0; i < nw; ++i) { + if (aw[i].v == (v^1)) { + we = &(aw[i]); + break; + } + } + ///mm_ol and mm_ou are used to make edge with long indel more easy to be cutted + mm_ol = MIN(ve->ol, we->ol); mm_ou = MIN(ve->ou, we->ou); + + for (i = kv = ol_max = ou_max = 0, vl_max = NULL; i < nv; ++i) { + if(av[i].del) continue; + kv++; + if((av[i].v == ve->v) || (!av[i].strong) || (av[i].ol <= ve->ol)) continue; + if(is_trio && get_tip_trio_infor(g, av[i].v) == ntrioF) continue; + + z = &(rev[av[i].v>>1]); wz = (ve->v^(av[i].v&1)); + for (m = 0; m < z->length; m++) { + if((z->buffer[m].tn<<1|(z->buffer[m].rev)) == wz) break; + } + if(m >= z->length) continue; + + if(ol_max < av[i].ol) ol_max = av[i].ol, vl_max = &(av[i]); + if(ou_max < av[i].ou) ou_max = av[i].ou; + } + if (kv < 1) continue; + if (kv >= 2) { + if (mm_ol > ol_max*len_rat) continue; + if (is_ou && mm_ou > ou_max*ou_rat && mm_ou > min_ou) continue; + if ((mm_ol + min_diff) > ol_max) continue; + } + + + for (i = kw = ol_max = ou_max = 0, wl_max = NULL; i < nw; ++i) { + if(aw[i].del) continue; + kw++; + if((aw[i].v == we->v) || (!aw[i].strong) || (aw[i].ol <= we->ol)) continue; + if(is_trio && get_tip_trio_infor(g, aw[i].v) == ntrioF) continue; + + z = &(rev[aw[i].v>>1]); wz = (we->v^(aw[i].v&1)); + for (m = 0; m < z->length; m++) { + if((z->buffer[m].tn<<1|(z->buffer[m].rev)) == wz) break; + } + if(m >= z->length) continue; + + if(ol_max < aw[i].ol) ol_max = aw[i].ol, wl_max = &(aw[i]); + if(ou_max < aw[i].ou) ou_max = aw[i].ou; + } + if (kw < 1) continue; + if (kw >= 2) { + if (mm_ol > ol_max*len_rat) continue; + if (is_ou && mm_ou > ou_max*ou_rat && mm_ou > min_ou) continue; + if ((mm_ol + min_diff) > ol_max) continue; + } + + if (kv <= 1 && kw <= 1) continue; + + to_del = 0; + if(is_topo) { + if (kv > 1 && kw > 1) { + to_del = 1; + } else if (kw == 1) { + if (asg_topocut_aux(g, w^1, max_ext) < max_ext) to_del = 1; + } else if (kv == 1) { + if (asg_topocut_aux(g, v^1, max_ext) < max_ext) to_del = 1; + } + } + + if(rev && rI) { + if((to_del == 0) && vl_max && (ve->v!=vl_max->v) && (trans_path_check(ve->v, vl_max->v, g, rev, rI, max_ext, b)==0)) { + to_del = 1; + } + if((to_del == 0) && wl_max && (we->v!=wl_max->v) && (trans_path_check(we->v, wl_max->v, g, rev, rI, max_ext, b)==0)) { + to_del = 1; + } + if(vl_max && wl_max) assert(ve->v!=vl_max->v||we->v!=wl_max->v); + } + + + if (to_del) { + ve->del = we->del = 1, ++cnt; + } + } + // stats_sysm(g); + if(!in) free(tx.a); + if (cnt > 0) asg_cleanup(g); +} + + void asg_arc_cut_length_adv(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio, uint32_t is_topo, uint32_t min_diff, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len) @@ -2784,7 +2951,9 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i // fprintf(stderr, "%.*s\tid::%u\tis_c::%u\n", // (int)Get_NAME_LENGTH(R_INF, 10819), Get_NAME(R_INF, 10819), 10819, is_contain_r((*rI), 10819)); // debug_info_of_specfic_node("m64011_190830_220126/47516220/ccs", sg, rI, "beg-0"); - // debug_info_of_specfic_node("m64012_190920_173625/163644465/ccs", sg, rI, "beg-0"); + // debug_info_of_specfic_node("bcb40bcc-d9cf-48e6-88ee-47ac3dde22ff", sg, rI, "beg-0"); + + if(asm_opt.is_ont) asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);///p_telo // fprintf(stderr, "[M::%s] count_edges_v_w(sg, 49778, 49847)->%ld\n", __func__, count_edges_v_w(sg, 49778, 49847)); @@ -2795,6 +2964,11 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i if(drop <= 0.500001) min_diff = step_diff>>1; else min_diff = step_diff; } + + if(asm_opt.is_ont) { + asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te); + } // fprintf(stderr, "(0):i->%ld, drop->%f\n", i, drop); // prt_specfic_sge(sg, 10531, 10519, "--0--"); @@ -2843,6 +3017,11 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i } } + if(asm_opt.is_ont) { + asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te); + } + if(is_ou) min_diff = step_diff; // print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty4.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0); // debug_info_of_specfic_node("m64012_190921_234837/111673711/ccs", sg, rI, "end");