Skip to content

Commit

Permalink
removed the pbread mode
Browse files Browse the repository at this point in the history 
It is not working well. Not at all.
  • Loading branch information
@lh3
lh3 committed Jan 27, 2016
1 parent 616a534 commit 1247dc2
Show file tree
Hide file tree
Showing 4 changed files with 7 additions and 59 deletions.
19 changes: 3 additions & 16 deletions bwamem.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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);
Expand Down
2 changes: 0 additions & 2 deletions bwamem.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
25 changes: 0 additions & 25 deletions bwamem_extra.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
20 changes: 4 additions & 16 deletions fastmap.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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);
Expand All @@ -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");
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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);
Expand Down

0 comments on commit 1247dc2

Please sign in to comment.