diff --git a/bwamem.c b/bwamem.c index b32c166..99a6d88 100644 --- a/bwamem.c +++ b/bwamem.c @@ -114,7 +114,7 @@ static void smem_aux_destroy(smem_aux_t *a) static void mem_collect_intv(const mem_opt_t *opt, const bwt_t *bwt, int len, const uint8_t *seq, smem_aux_t *a) { int i, k, x = 0, old_n; - int start_width = (opt->flag & MEM_F_SELF_OVLP)? 2 : 1; + int start_width = 1; int split_len = (int)(opt->min_seed_len * opt->split_factor + .499); a->mem.n = 0; // first pass: find all SMEMs @@ -488,13 +488,6 @@ int mem_sort_dedup_patch(const mem_opt_t *opt, const bntseq_t *bns, const uint8_ return m; } -int mem_test_and_remove_exact(const mem_opt_t *opt, int n, mem_alnreg_t *a, int qlen) -{ - if (!(opt->flag & MEM_F_SELF_OVLP) || n == 0 || a->truesc != qlen * opt->a) return n; - memmove(a, a + 1, (n - 1) * sizeof(mem_alnreg_t)); - return n - 1; -} - typedef kvec_t(int) int_v; static void mem_mark_primary_se_core(const mem_opt_t *opt, int n, mem_alnreg_t *a, int_v *z) @@ -1046,8 +1039,6 @@ mem_alnreg_v mem_align1_core(const mem_opt_t *opt, const bwt_t *bwt, const bntse } free(chn.a); regs.n = mem_sort_dedup_patch(opt, bns, pac, (uint8_t*)seq, regs.n, regs.a); - if (opt->flag & MEM_F_SELF_OVLP) - regs.n = mem_test_and_remove_exact(opt, regs.n, regs.a, l_seq); if (bwa_verbose >= 4) { err_printf("* %ld chains remain after removing duplicated chains\n", regs.n); for (i = 0; i < regs.n; ++i) { @@ -1168,12 +1159,8 @@ static void worker2(void *data, int i, int tid) worker_t *w = (worker_t*)data; if (!(w->opt->flag&MEM_F_PE)) { if (bwa_verbose >= 4) printf("=====> Finalizing read '%s' <=====\n", w->seqs[i].name); - if (w->opt->flag & MEM_F_ALN_REG) { - mem_reg2ovlp(w->opt, w->bns, w->pac, &w->seqs[i], &w->regs[i]); - } else { - mem_mark_primary_se(w->opt, w->regs[i].n, w->regs[i].a, w->n_processed + i); - mem_reg2sam(w->opt, w->bns, w->pac, &w->seqs[i], &w->regs[i], 0, 0); - } + mem_mark_primary_se(w->opt, w->regs[i].n, w->regs[i].a, w->n_processed + i); + mem_reg2sam(w->opt, w->bns, w->pac, &w->seqs[i], &w->regs[i], 0, 0); free(w->regs[i].a); } else { if (bwa_verbose >= 4) printf("=====> Finalizing read pair '%s' <=====\n", w->seqs[i<<1|0].name); diff --git a/bwamem.h b/bwamem.h index 8ffe729..730b603 100644 --- a/bwamem.h +++ b/bwamem.h @@ -16,8 +16,6 @@ typedef struct __smem_i smem_i; #define MEM_F_ALL 0x8 #define MEM_F_NO_MULTI 0x10 #define MEM_F_NO_RESCUE 0x20 -#define MEM_F_SELF_OVLP 0x40 -#define MEM_F_ALN_REG 0x80 #define MEM_F_REF_HDR 0x100 #define MEM_F_SOFTCLIP 0x200 #define MEM_F_SMARTPE 0x400 diff --git a/bwamem_extra.c b/bwamem_extra.c index 02e817a..d2e37f4 100644 --- a/bwamem_extra.c +++ b/bwamem_extra.c @@ -87,31 +87,6 @@ mem_alnreg_v mem_align1(const mem_opt_t *opt, const bwt_t *bwt, const bntseq_t * return ar; } -void mem_reg2ovlp(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, bseq1_t *s, mem_alnreg_v *a) -{ - int i; - kstring_t str = {0,0,0}; - for (i = 0; i < a->n; ++i) { - const mem_alnreg_t *p = &a->a[i]; - int is_rev, rid, qb = p->qb, qe = p->qe; - int64_t pos, rb = p->rb, re = p->re; - pos = bns_depos(bns, rb < bns->l_pac? rb : re - 1, &is_rev); - rid = bns_pos2rid(bns, pos); - assert(rid == p->rid); - pos -= bns->anns[rid].offset; - kputs(s->name, &str); kputc('\t', &str); - kputw(s->l_seq, &str); kputc('\t', &str); - if (is_rev) qb ^= qe, qe ^= qb, qb ^= qe; // swap - kputw(qb, &str); kputc('\t', &str); kputw(qe, &str); kputc('\t', &str); - kputs(bns->anns[rid].name, &str); kputc('\t', &str); - kputw(bns->anns[rid].len, &str); kputc('\t', &str); - kputw(pos, &str); kputc('\t', &str); kputw(pos + (re - rb), &str); kputc('\t', &str); - ksprintf(&str, "%.3f", (double)p->truesc / opt->a / (qe - qb > re - rb? qe - qb : re - rb)); - kputc('\n', &str); - } - s->sam = str.s; -} - static inline int get_pri_idx(double XA_drop_ratio, const mem_alnreg_t *a, int i) { int k = a[i].secondary_all; diff --git a/fastmap.c b/fastmap.c index 3cca4de..912e31d 100644 --- a/fastmap.c +++ b/fastmap.c @@ -130,7 +130,7 @@ int main_mem(int argc, char *argv[]) aux.opt = opt = mem_opt_init(); memset(&opt0, 0, sizeof(mem_opt_t)); - while ((c = getopt(argc, argv, "1epaFMCSPVYjk:c:v:s:r:t:R:A:B:O:E:U:w:L:d:T:Q:D:m:I:N:W:x:G:h:y:K:X:H:")) >= 0) { + while ((c = getopt(argc, argv, "1paMCSPVYjk:c:v:s:r:t:R:A:B:O:E:U:w:L:d:T:Q:D:m:I:N:W:x:G:h:y:K:X:H:")) >= 0) { if (c == 'k') opt->min_seed_len = atoi(optarg), opt0.min_seed_len = 1; else if (c == '1') no_mt_io = 1; else if (c == 'x') mode = optarg; @@ -145,8 +145,6 @@ int main_mem(int argc, char *argv[]) else if (c == 'p') opt->flag |= MEM_F_PE | MEM_F_SMARTPE; else if (c == 'M') opt->flag |= MEM_F_NO_MULTI; else if (c == 'S') opt->flag |= MEM_F_NO_RESCUE; - else if (c == 'e') opt->flag |= MEM_F_SELF_OVLP; - else if (c == 'F') opt->flag |= MEM_F_ALN_REG; else if (c == 'Y') opt->flag |= MEM_F_SOFTCLIP; else if (c == 'V') opt->flag |= MEM_F_REF_HDR; else if (c == 'c') opt->max_occ = atoi(optarg), opt0.max_occ = 1; @@ -251,7 +249,6 @@ int main_mem(int argc, char *argv[]) fprintf(stderr, " -m INT perform at most INT rounds of mate rescues for each read [%d]\n", opt->max_matesw); fprintf(stderr, " -S skip mate rescue\n"); fprintf(stderr, " -P skip pairing; mate rescue performed unless -S also in use\n"); - fprintf(stderr, " -e discard full-length exact matches\n"); fprintf(stderr, "\nScoring options:\n\n"); fprintf(stderr, " -A INT score for a sequence match, which scales options -TdBOELU unless overridden [%d]\n", opt->a); fprintf(stderr, " -B INT penalty for a mismatch [%d]\n", opt->b); @@ -263,7 +260,6 @@ int main_mem(int argc, char *argv[]) fprintf(stderr, " pacbio: -k17 -W40 -r10 -A1 -B1 -O1 -E1 -L0 (PacBio reads to ref)\n"); fprintf(stderr, " ont2d: -k14 -W20 -r10 -A1 -B1 -O1 -E1 -L0 (Oxford Nanopore 2D-reads to ref)\n"); fprintf(stderr, " intractg: -B9 -O16 -L5 (intra-species contigs to ref)\n"); -// fprintf(stderr, " pbread: -k13 -W40 -c1000 -r10 -A1 -B1 -O1 -E1 -N25 -FeaD.001\n"); fprintf(stderr, "\nInput/output options:\n\n"); fprintf(stderr, " -p smart pairing (ignoring in2.fq)\n"); fprintf(stderr, " -R STR read group header line such as '@RG\\tID:foo\\tSM:bar' [null]\n"); @@ -296,21 +292,14 @@ int main_mem(int argc, char *argv[]) if (!opt0.b) opt->b = 9; if (!opt0.pen_clip5) opt->pen_clip5 = 5; if (!opt0.pen_clip3) opt->pen_clip3 = 5; - } else if (strcmp(mode, "pacbio") == 0 || strcmp(mode, "pbref") == 0 || strcmp(mode, "pbread") == 0 || strcmp(mode, "ont2d") == 0) { + } else if (strcmp(mode, "pacbio") == 0 || strcmp(mode, "pbref") == 0 || strcmp(mode, "ont2d") == 0) { if (!opt0.o_del) opt->o_del = 1; if (!opt0.e_del) opt->e_del = 1; if (!opt0.o_ins) opt->o_ins = 1; if (!opt0.e_ins) opt->e_ins = 1; if (!opt0.b) opt->b = 1; if (opt0.split_factor == 0.) opt->split_factor = 10.; - if (strcmp(mode, "pbread") == 0) { // pacbio read-to-read setting; NOT working well! - opt->flag |= MEM_F_ALL | MEM_F_SELF_OVLP | MEM_F_ALN_REG; - if (!opt0.min_chain_weight) opt->min_chain_weight = 40; - if (!opt0.max_occ) opt->max_occ = 1000; - if (!opt0.min_seed_len) opt->min_seed_len = 13; - if (!opt0.max_chain_extend) opt->max_chain_extend = 25; - if (opt0.drop_ratio == 0.) opt->drop_ratio = .001; - } else if (strcmp(mode, "ont2d") == 0) { + if (strcmp(mode, "ont2d") == 0) { if (!opt0.min_chain_weight) opt->min_chain_weight = 20; if (!opt0.min_seed_len) opt->min_seed_len = 14; if (!opt0.pen_clip5) opt->pen_clip5 = 0; @@ -359,8 +348,7 @@ int main_mem(int argc, char *argv[]) opt->flag |= MEM_F_PE; } } - if (!(opt->flag & MEM_F_ALN_REG)) - bwa_print_sam_hdr(aux.idx->bns, hdr_line); + bwa_print_sam_hdr(aux.idx->bns, hdr_line); aux.actual_chunk_size = fixed_chunk_size > 0? fixed_chunk_size : opt->chunk_size * opt->n_threads; kt_pipeline(no_mt_io? 1 : 2, process, &aux, 3); free(hdr_line);