diff --git a/CommandLines.cpp b/CommandLines.cpp index 2af9327..256a35a 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -79,6 +79,8 @@ static ko_longopt_t long_options[] = { { "chem-c", ko_required_argument, 361}, { "chem-f", ko_required_argument, 362}, { "ul-m", ko_required_argument, 363}, + { "rl-cut", ko_required_argument, 364}, + { "sc-cut", ko_required_argument, 365}, // { "path-round", ko_required_argument, 348}, { 0, 0, 0 } }; @@ -228,6 +230,10 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " detect chimeric 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 chimeric read detection [%lu]\n", asm_opt->chemical_flank); + fprintf(stderr, " --rl-cut INT\n"); + fprintf(stderr, " filter out ONT simplex reads shorter than for assembly [%ld]\n", asm_opt->rl_cut); + fprintf(stderr, " --sc-cut INT\n"); + fprintf(stderr, " filter out ONT simplex reads with a mean base quality score below [%ld]\n", asm_opt->sc_cut); fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n"); @@ -364,6 +370,9 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->chemical_cov = 1; asm_opt->chemical_flank = 256; asm_opt->ul_mod = 0; + + asm_opt->rl_cut = 1000; + asm_opt->sc_cut = 10; } void destory_enzyme(enzyme* f) @@ -990,6 +999,10 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) } else if (c == 363) { asm_opt->ul_mod = atol(opt.arg); **/ + } else if (c == 364) { + asm_opt->rl_cut = atol(opt.arg); + } else if (c == 365) { + asm_opt->sc_cut = 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); } @@ -1027,6 +1040,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) // fprintf(stderr, "[M::%s::] post join::%u\n", __func__, (uint32_t)(!(asm_opt->flag & HA_F_BAN_POST_JOIN))); // exit(1); + if(!(asm_opt->is_ont)) { + asm_opt->rl_cut = -1; asm_opt->sc_cut = 1; + } return check_option(asm_opt); } diff --git a/CommandLines.h b/CommandLines.h index 511a5f5..3bebc87 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.25.0-r721" +#define HA_VERSION "0.25.0-r726" #define VERBOSE 0 @@ -166,6 +166,10 @@ typedef struct { uint64_t is_sc; uint64_t chemical_cov; uint64_t chemical_flank; + + int64_t rl_cut; + int64_t sc_cut; + } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/htab.cpp b/htab.cpp index d268df4..bd16ded 100644 --- a/htab.cpp +++ b/htab.cpp @@ -546,6 +546,21 @@ const int ha_pt_cnt(const ha_pt_t *h, uint64_t hash) return kh_key(g->h, k) & YAK_MAX_COUNT; } +inline uint64_t flt_quals(char *sc_a, uint64_t sc_l, uint64_t sc_off, int64_t sc_cut) +{ + int64_t sc_min = sc_l * sc_cut, sc_tot; uint64_t k; + for (k = sc_tot = 0; (k < sc_l) && (sc_tot < sc_min); k++) { + sc_tot += (((uint8_t)sc_a[k]) - sc_off); + } + + // if(sc_tot < sc_min) { + // fprintf(stderr, "[M::%s] sc_tot::%ld, sc_min::%ld, sc_l::%lu\n", __func__, sc_tot, sc_min, sc_l); + // } + + if(sc_tot < sc_min) return 0; + return 1; +} + /********************************** * Buffer for counting all k-mers * **********************************/ @@ -745,7 +760,8 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for } else {\ while ((ret = kseq_read(p->ks)) >= 0) {\ int l = (int)(p->ks->seq.l) - (int)(p->opt->adaLen) - (int)(p->opt->adaLen);\ - if(l <= 0) continue;\ + if((l <= 0) || (l < asm_opt.rl_cut)) continue;\ + if((asm_opt.is_sc) && (asm_opt.sc_cut > 0) && (!flt_quals(p->ks->qual.s+p->opt->adaLen, l, 33, asm_opt.sc_cut))) continue;\ if (p->n_seq >= 1<<28) {\ fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);\ exit(1);\