diff --git a/bwa.c b/bwa.c index aff9aff..c8400b1 100644 --- a/bwa.c +++ b/bwa.c @@ -25,7 +25,7 @@ static inline void trim_readno(kstring_t *s) static inline void kseq2bseq1(const kseq_t *ks, bseq1_t *s) { // TODO: it would be better to allocate one chunk of memory, but probably it does not matter in practice s->name = strdup(ks->name.s); - s->comment = ks->comment.l? strdup(s->comment) : 0; + s->comment = ks->comment.l? strdup(ks->comment.s) : 0; s->seq = strdup(ks->seq.s); s->qual = ks->qual.l? strdup(ks->qual.s) : 0; s->l_seq = strlen(s->seq); diff --git a/bwamem.c b/bwamem.c index ce55cad..553fe1c 100644 --- a/bwamem.c +++ b/bwamem.c @@ -13,8 +13,6 @@ #include "kvec.h" #include "ksort.h" -int mem_verbose = 3; // 1: error only; 2: error+warning; 3: message+error+warning; >=4: debugging - void mem_fill_scmat(int a, int b, int8_t mat[25]) { int i, j, k; @@ -462,7 +460,7 @@ void mem_chain2aln(const mem_opt_t *opt, int64_t l_pac, const uint8_t *pac, int a->score = ksw_extend(l_query - qe, query + qe, rmax[1] - rmax[0] - re, rseq + re, 5, opt->mat, opt->q, opt->r, opt->w, a->score, &qle, &tle); a->qe = qe + qle; a->re = rmax[0] + re + tle; } else a->qe = l_query, a->re = s->rbeg + s->len; - if (mem_verbose >= 4) printf("[%d] score=%d\t[%d,%d) <=> [%ld,%ld)\n", k, a->score, a->qb, a->qe, (long)a->rb, (long)a->re); + if (bwa_verbose >= 4) printf("[%d] score=%d\t[%d,%d) <=> [%ld,%ld)\n", k, a->score, a->qb, a->qe, (long)a->rb, (long)a->re); // compute seedcov for (i = 0, a->seedcov = 0; i < c->n; ++i) { const mem_seed_t *t = &c->seeds[i]; @@ -572,6 +570,7 @@ void bwa_hit2sam(kstring_t *str, const int8_t mat[25], int q, int r, int w, cons if (p->score >= 0) { kputsn("\tAS:i:", 6, str); kputw(p->score, str); } if (p->sub >= 0) { kputsn("\tXS:i:", 6, str); kputw(p->sub, str); } if (bwa_rg_id) { kputsn("\tRG:i:", 6, str); kputs(bwa_rg_id, str); } + if (s->comment) { kputc('\t', str); kputs(s->comment, str); } kputc('\n', str); free(cigar); #undef is_mapped @@ -633,7 +632,7 @@ static mem_alnreg_v find_alnreg(const mem_opt_t *opt, const bwt_t *bwt, const bn s->seq[i] = nst_nt4_table[(int)s->seq[i]]; chn = mem_chain(opt, bwt, s->l_seq, (uint8_t*)s->seq); chn.n = mem_chain_flt(opt, chn.n, chn.a); - if (mem_verbose >= 4) mem_print_chain(bns, &chn); + if (bwa_verbose >= 4) mem_print_chain(bns, &chn); kv_init(regs); kv_init(tmp); for (i = 0; i < chn.n; ++i) { mem_chain2aln(opt, bns->l_pac, pac, s->l_seq, (uint8_t*)s->seq, &chn.a[i], &tmp); diff --git a/bwamem_pair.c b/bwamem_pair.c index 3dce119..57a128a 100644 --- a/bwamem_pair.c +++ b/bwamem_pair.c @@ -47,7 +47,7 @@ void mem_pestat(const mem_opt_t *opt, int64_t l_pac, int n, const mem_alnreg_v * dir = mem_infer_dir(l_pac, r[0]->a[0].rb, r[1]->a[0].rb, &is); if (is && is <= opt->max_ins) kv_push(uint64_t, isize[dir], is); } - if (mem_verbose >= 3) fprintf(stderr, "[M::%s] # candidate unique pairs for (FF, FR, RF, RR): (%ld, %ld, %ld, %ld)\n", __func__, isize[0].n, isize[1].n, isize[2].n, isize[3].n); + if (bwa_verbose >= 3) fprintf(stderr, "[M::%s] # candidate unique pairs for (FF, FR, RF, RR): (%ld, %ld, %ld, %ld)\n", __func__, isize[0].n, isize[1].n, isize[2].n, isize[3].n); for (d = 0; d < 4; ++d) { // TODO: this block is nearly identical to the one in bwtsw2_pair.c. It would be better to merge these two. mem_pestat_t *r = &pes[d]; uint64_v *q = &isize[d]; diff --git a/fastmap.c b/fastmap.c index adbe04c..437192a 100644 --- a/fastmap.c +++ b/fastmap.c @@ -14,7 +14,7 @@ extern unsigned char nst_nt4_table[256]; int main_mem(int argc, char *argv[]) { mem_opt_t *opt; - int c, n; + int i, c, n, copy_comment = 0; gzFile fp, fp2 = 0; kseq_t *ks, *ks2 = 0; bseq1_t *seqs; @@ -22,14 +22,15 @@ int main_mem(int argc, char *argv[]) char *rg_line = 0; opt = mem_opt_init(); - while ((c = getopt(argc, argv, "PHk:c:v:s:r:t:R:")) >= 0) { + while ((c = getopt(argc, argv, "CPHk:c:v:s:r:t:R:")) >= 0) { if (c == 'k') opt->min_seed_len = atoi(optarg); else if (c == 't') opt->n_threads = atoi(optarg), opt->n_threads = opt->n_threads > 1? opt->n_threads : 1; else if (c == 'P') opt->flag |= MEM_F_NOPAIRING; else if (c == 'H') opt->flag |= MEM_F_HARDCLIP; else if (c == 'c') opt->max_occ = atoi(optarg); - else if (c == 'v') mem_verbose = atoi(optarg); + else if (c == 'v') bwa_verbose = atoi(optarg); else if (c == 'r') opt->split_factor = atof(optarg); + else if (c == 'C') copy_comment = 1; else if (c == 'R') { if ((rg_line = bwa_set_rg(optarg)) == 0) return 1; // FIXME: memory leak } else if (c == 's') opt->split_width = atoi(optarg); @@ -43,7 +44,8 @@ int main_mem(int argc, char *argv[]) fprintf(stderr, " -s INT look for internal seeds inside a seed with less than INT occ [%d]\n", opt->split_width); fprintf(stderr, " -r FLOAT look for internal seeds inside a seed longer than {-k} * FLOAT [%g]\n", opt->split_factor); fprintf(stderr, " -R STR read group header line such as '@RG\tID:foo\tSM:bar' [null]\n"); - fprintf(stderr, " -v INT verbose level [%d]\n", mem_verbose); + fprintf(stderr, " -v INT verbose level [%d]\n", bwa_verbose); + fprintf(stderr, " -C append FASTA/FASTQ comment to SAM output\n"); fprintf(stderr, "\n"); free(opt); return 1; @@ -60,6 +62,14 @@ int main_mem(int argc, char *argv[]) opt->flag |= MEM_F_PE; } while ((seqs = bseq_read(opt->chunk_size, &n, ks, ks2)) != 0) { + int64_t size = 0; + if (!copy_comment) + for (i = 0; i < n; ++i) { + free(seqs[i].comment); seqs[i].comment = 0; + } + for (i = 0; i < n; ++i) size += seqs[i].l_seq; + if (bwa_verbose >= 3) + fprintf(stderr, "[M::%s] read %d sequences (%ld bp)...\n", __func__, n, (long)size); mem_process_seqs(opt, idx->bwt, idx->bns, idx->pac, n, seqs); free(seqs); }