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

output multiple hits

This commit is contained in:
Heng Li 2013-02-24 13:23:43 -05:00
parent 6bdccf2a8a
commit 85775c3384
3 changed files with 12 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
bwamem.c
View File

@ -554,7 +554,9 @@ void bwa_hit2sam(kstring_t *str, const int8_t mat[25], int q, int r, int w, cons
} else kputw(0, str); } else kputw(0, str);
kputc('\t', str); kputc('\t', str);
} else kputsn("\t*\t0\t0\t", 7, str); } else kputsn("\t*\t0\t0\t", 7, str);
if (!(p->flag&0x10)) { // print SEQ and QUAL, the forward strand if (p->flag&0x100) { // for secondary alignments, don't write SEQ and QUAL
kputsn("*\t*", 3, str);
} else if (!(p->flag&0x10)) { // print SEQ and QUAL, the forward strand
int i, qb = 0, qe = s->l_seq; int i, qb = 0, qe = s->l_seq;
if (!(p->flag&4) && is_hard) qb = p->qb, qe = p->qe; if (!(p->flag&4) && is_hard) qb = p->qb, qe = p->qe;
ks_resize(str, str->l + (qe - qb) + 1); ks_resize(str, str->l + (qe - qb) + 1);
@ -610,7 +612,7 @@ void mem_alnreg2hit(const mem_alnreg_t *a, bwahit_t *h)
{ {
h->rb = a->rb; h->re = a->re; h->qb = a->qb; h->qe = a->qe; h->rb = a->rb; h->re = a->re; h->qb = a->qb; h->qe = a->qe;
h->score = a->score; h->score = a->score;
h->sub = a->sub > a->csub? a->sub : a->csub; h->sub = a->secondary >= 0? -1 : a->sub > a->csub? a->sub : a->csub;
h->qual = 0; // quality unset h->qual = 0; // quality unset
h->flag = a->secondary >= 0? 0x100 : 0; // only the "secondary" bit is set h->flag = a->secondary >= 0? 0x100 : 0; // only the "secondary" bit is set
} }
@ -623,10 +625,10 @@ void mem_sam_se(const mem_opt_t *opt, const bntseq_t *bns, const uint8_t *pac, b
if (a->n > 0) { if (a->n > 0) {
for (k = 0; k < a->n; ++k) { for (k = 0; k < a->n; ++k) {
bwahit_t h; bwahit_t h;
if (a->a[k].secondary >= 0) continue; if (a->a[k].secondary >= 0 && !(opt->flag&MEM_F_ALL)) continue;
mem_alnreg2hit(&a->a[k], &h); mem_alnreg2hit(&a->a[k], &h);
h.flag |= extra_flag; h.flag |= extra_flag;
h.qual = mem_approx_mapq_se(opt, &a->a[k]); h.qual = a->a[k].secondary >= 0? 0 : mem_approx_mapq_se(opt, &a->a[k]);
bwa_hit2sam(&str, opt->mat, opt->q, opt->r, opt->w, bns, pac, s, &h, opt->flag&MEM_F_HARDCLIP, m); bwa_hit2sam(&str, opt->mat, opt->q, opt->r, opt->w, bns, pac, s, &h, opt->flag&MEM_F_HARDCLIP, m);
} }
} else bwa_hit2sam(&str, opt->mat, opt->q, opt->r, opt->w, bns, pac, s, 0, opt->flag&MEM_F_HARDCLIP, m); } else bwa_hit2sam(&str, opt->mat, opt->q, opt->r, opt->w, bns, pac, s, 0, opt->flag&MEM_F_HARDCLIP, m);

1
bwamem.h
View File

@ -14,6 +14,7 @@ typedef struct __smem_i smem_i;
#define MEM_F_HARDCLIP 0x1 #define MEM_F_HARDCLIP 0x1
#define MEM_F_PE 0x2 #define MEM_F_PE 0x2
#define MEM_F_NOPAIRING 0x4 #define MEM_F_NOPAIRING 0x4
#define MEM_F_ALL 0x8
typedef struct { typedef struct {
int a, b, q, r, w; int a, b, q, r, w;

6
fastmap.c
View File

@ -26,11 +26,12 @@ int main_mem(int argc, char *argv[])
void *ko = 0, *ko2 = 0; void *ko = 0, *ko2 = 0;
opt = mem_opt_init(); opt = mem_opt_init();
while ((c = getopt(argc, argv, "CPHk:c:v:s:r:t:R:")) >= 0) { while ((c = getopt(argc, argv, "aCPHk:c:v:s:r:t:R:")) >= 0) {
if (c == 'k') opt->min_seed_len = atoi(optarg); 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 == '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 == 'P') opt->flag |= MEM_F_NOPAIRING;
else if (c == 'H') opt->flag |= MEM_F_HARDCLIP; else if (c == 'H') opt->flag |= MEM_F_HARDCLIP;
else if (c == 'a') opt->flag |= MEM_F_ALL;
else if (c == 'c') opt->max_occ = atoi(optarg); else if (c == 'c') opt->max_occ = atoi(optarg);
else if (c == 'v') bwa_verbose = atoi(optarg); else if (c == 'v') bwa_verbose = atoi(optarg);
else if (c == 'r') opt->split_factor = atof(optarg); else if (c == 'r') opt->split_factor = atof(optarg);
@ -49,6 +50,9 @@ int main_mem(int argc, char *argv[])
fprintf(stderr, " -r FLOAT look for internal seeds inside a seed longer than {-k} * FLOAT [%g]\n", opt->split_factor); 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, " -R STR read group header line such as '@RG\tID:foo\tSM:bar' [null]\n");
fprintf(stderr, " -v INT verbose level [%d]\n", bwa_verbose); fprintf(stderr, " -v INT verbose level [%d]\n", bwa_verbose);
fprintf(stderr, " -a output all alignments for SE or unpaired PE\n");
fprintf(stderr, " -P perform mate SW only but skip pairing\n");
fprintf(stderr, " -H hard clipping\n");
fprintf(stderr, " -C append FASTA/FASTQ comment to SAM output\n"); fprintf(stderr, " -C append FASTA/FASTQ comment to SAM output\n");
fprintf(stderr, "\n"); fprintf(stderr, "\n");
free(opt); free(opt);