zzh
/
fast-bwa
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

append fasta/q comment

This commit is contained in:
Heng Li 2013-02-23 16:57:34 -05:00
parent 33236de32e
commit b4c38bcc1c
4 changed files with 19 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
bwa.c
View File

@ -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);

7
bwamem.c
View File

@ -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);

2
bwamem_pair.c
View File

@ -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];

18
fastmap.c
View File

@ -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);
}