applied patches from Alec Wysoker

This commit is contained in:
Heng Li 2011-05-04 09:46:50 -04:00
parent 4d064c69ce
commit 243e735431
5 changed files with 123 additions and 38 deletions

68
bwase.c
View File

@ -437,15 +437,15 @@ void bwa_print_sam1(const bntseq_t *bns, bwa_seq_t *p, const bwa_seq_t *mate, in
if (mate->strand) flag |= SAM_FMR; if (mate->strand) flag |= SAM_FMR;
} else flag |= SAM_FMU; } else flag |= SAM_FMU;
} }
printf("%s\t%d\t%s\t", p->name, flag, bns->anns[seqid].name); err_printf("%s\t%d\t%s\t", p->name, flag, bns->anns[seqid].name);
printf("%d\t%d\t", (int)(p->pos - bns->anns[seqid].offset + 1), p->mapQ); err_printf("%d\t%d\t", (int)(p->pos - bns->anns[seqid].offset + 1), p->mapQ);
// print CIGAR // print CIGAR
if (p->cigar) { if (p->cigar) {
for (j = 0; j != p->n_cigar; ++j) for (j = 0; j != p->n_cigar; ++j)
printf("%d%c", __cigar_len(p->cigar[j]), "MIDS"[__cigar_op(p->cigar[j])]); err_printf("%d%c", __cigar_len(p->cigar[j]), "MIDS"[__cigar_op(p->cigar[j])]);
} else if (p->type == BWA_TYPE_NO_MATCH) printf("*"); } else if (p->type == BWA_TYPE_NO_MATCH) err_printf("*");
else printf("%dM", p->len); else err_printf("%dM", p->len);
// print mate coordinate // print mate coordinate
if (mate && mate->type != BWA_TYPE_NO_MATCH) { if (mate && mate->type != BWA_TYPE_NO_MATCH) {
@ -454,12 +454,12 @@ void bwa_print_sam1(const bntseq_t *bns, bwa_seq_t *p, const bwa_seq_t *mate, in
am = mate->seQ < p->seQ? mate->seQ : p->seQ; // smaller single-end mapping quality am = mate->seQ < p->seQ? mate->seQ : p->seQ; // smaller single-end mapping quality
// redundant calculation here, but should not matter too much // redundant calculation here, but should not matter too much
m_is_N = bns_coor_pac2real(bns, mate->pos, mate->len, &m_seqid); m_is_N = bns_coor_pac2real(bns, mate->pos, mate->len, &m_seqid);
printf("\t%s\t", (seqid == m_seqid)? "=" : bns->anns[m_seqid].name); err_printf("\t%s\t", (seqid == m_seqid)? "=" : bns->anns[m_seqid].name);
isize = (seqid == m_seqid)? pos_5(mate) - pos_5(p) : 0; isize = (seqid == m_seqid)? pos_5(mate) - pos_5(p) : 0;
if (p->type == BWA_TYPE_NO_MATCH) isize = 0; if (p->type == BWA_TYPE_NO_MATCH) isize = 0;
printf("%d\t%lld\t", (int)(mate->pos - bns->anns[m_seqid].offset + 1), isize); err_printf("%d\t%lld\t", (int)(mate->pos - bns->anns[m_seqid].offset + 1), isize);
} else if (mate) printf("\t=\t%d\t0\t", (int)(p->pos - bns->anns[seqid].offset + 1)); } else if (mate) err_printf("\t=\t%d\t0\t", (int)(p->pos - bns->anns[seqid].offset + 1));
else printf("\t*\t0\t0\t"); else err_printf("\t*\t0\t0\t");
// print sequence and quality // print sequence and quality
if (p->strand == 0) if (p->strand == 0)
@ -468,42 +468,42 @@ void bwa_print_sam1(const bntseq_t *bns, bwa_seq_t *p, const bwa_seq_t *mate, in
putchar('\t'); putchar('\t');
if (p->qual) { if (p->qual) {
if (p->strand) seq_reverse(p->len, p->qual, 0); // reverse quality if (p->strand) seq_reverse(p->len, p->qual, 0); // reverse quality
printf("%s", p->qual); err_printf("%s", p->qual);
} else printf("*"); } else err_printf("*");
if (bwa_rg_id) printf("\tRG:Z:%s", bwa_rg_id); if (bwa_rg_id) err_printf("\tRG:Z:%s", bwa_rg_id);
if (p->bc[0]) printf("\tBC:Z:%s", p->bc); if (p->bc[0]) err_printf("\tBC:Z:%s", p->bc);
if (p->clip_len < p->full_len) printf("\tXC:i:%d", p->clip_len); if (p->clip_len < p->full_len) err_printf("\tXC:i:%d", p->clip_len);
if (p->type != BWA_TYPE_NO_MATCH) { if (p->type != BWA_TYPE_NO_MATCH) {
int i; int i;
// calculate XT tag // calculate XT tag
XT = "NURM"[p->type]; XT = "NURM"[p->type];
if (nn > 10) XT = 'N'; if (nn > 10) XT = 'N';
// print tags // print tags
printf("\tXT:A:%c\t%s:i:%d", XT, (mode & BWA_MODE_COMPREAD)? "NM" : "CM", p->nm); err_printf("\tXT:A:%c\t%s:i:%d", XT, (mode & BWA_MODE_COMPREAD)? "NM" : "CM", p->nm);
if (nn) printf("\tXN:i:%d", nn); if (nn) err_printf("\tXN:i:%d", nn);
if (mate) printf("\tSM:i:%d\tAM:i:%d", p->seQ, am); if (mate) err_printf("\tSM:i:%d\tAM:i:%d", p->seQ, am);
if (p->type != BWA_TYPE_MATESW) { // X0 and X1 are not available for this type of alignment if (p->type != BWA_TYPE_MATESW) { // X0 and X1 are not available for this type of alignment
printf("\tX0:i:%d", p->c1); err_printf("\tX0:i:%d", p->c1);
if (p->c1 <= max_top2) printf("\tX1:i:%d", p->c2); if (p->c1 <= max_top2) err_printf("\tX1:i:%d", p->c2);
} }
printf("\tXM:i:%d\tXO:i:%d\tXG:i:%d", p->n_mm, p->n_gapo, p->n_gapo+p->n_gape); err_printf("\tXM:i:%d\tXO:i:%d\tXG:i:%d", p->n_mm, p->n_gapo, p->n_gapo+p->n_gape);
if (p->md) printf("\tMD:Z:%s", p->md); if (p->md) err_printf("\tMD:Z:%s", p->md);
// print multiple hits // print multiple hits
if (p->n_multi) { if (p->n_multi) {
printf("\tXA:Z:"); err_printf("\tXA:Z:");
for (i = 0; i < p->n_multi; ++i) { for (i = 0; i < p->n_multi; ++i) {
bwt_multi1_t *q = p->multi + i; bwt_multi1_t *q = p->multi + i;
int k; int k;
j = pos_end_multi(q, p->len) - q->pos; j = pos_end_multi(q, p->len) - q->pos;
nn = bns_coor_pac2real(bns, q->pos, j, &seqid); nn = bns_coor_pac2real(bns, q->pos, j, &seqid);
printf("%s,%c%d,", bns->anns[seqid].name, q->strand? '-' : '+', err_printf("%s,%c%d,", bns->anns[seqid].name, q->strand? '-' : '+',
(int)(q->pos - bns->anns[seqid].offset + 1)); (int)(q->pos - bns->anns[seqid].offset + 1));
if (q->cigar) { if (q->cigar) {
for (k = 0; k < q->n_cigar; ++k) for (k = 0; k < q->n_cigar; ++k)
printf("%d%c", __cigar_len(q->cigar[k]), "MIDS"[__cigar_op(q->cigar[k])]); err_printf("%d%c", __cigar_len(q->cigar[k]), "MIDS"[__cigar_op(q->cigar[k])]);
} else printf("%dM", p->len); } else err_printf("%dM", p->len);
printf(",%d;", q->gap + q->mm); err_printf(",%d;", q->gap + q->mm);
} }
} }
} }
@ -512,16 +512,16 @@ void bwa_print_sam1(const bntseq_t *bns, bwa_seq_t *p, const bwa_seq_t *mate, in
ubyte_t *s = p->strand? p->rseq : p->seq; ubyte_t *s = p->strand? p->rseq : p->seq;
int flag = p->extra_flag | SAM_FSU; int flag = p->extra_flag | SAM_FSU;
if (mate && mate->type == BWA_TYPE_NO_MATCH) flag |= SAM_FMU; if (mate && mate->type == BWA_TYPE_NO_MATCH) flag |= SAM_FMU;
printf("%s\t%d\t*\t0\t0\t*\t*\t0\t0\t", p->name, flag); err_printf("%s\t%d\t*\t0\t0\t*\t*\t0\t0\t", p->name, flag);
for (j = 0; j != p->len; ++j) putchar("ACGTN"[(int)s[j]]); for (j = 0; j != p->len; ++j) putchar("ACGTN"[(int)s[j]]);
putchar('\t'); putchar('\t');
if (p->qual) { if (p->qual) {
if (p->strand) seq_reverse(p->len, p->qual, 0); // reverse quality if (p->strand) seq_reverse(p->len, p->qual, 0); // reverse quality
printf("%s", p->qual); err_printf("%s", p->qual);
} else printf("*"); } else err_printf("*");
if (bwa_rg_id) printf("\tRG:Z:%s", bwa_rg_id); if (bwa_rg_id) err_printf("\tRG:Z:%s", bwa_rg_id);
if (p->bc[0]) printf("\tBC:Z:%s", p->bc); if (p->bc[0]) err_printf("\tBC:Z:%s", p->bc);
if (p->clip_len < p->full_len) printf("\tXC:i:%d", p->clip_len); if (p->clip_len < p->full_len) err_printf("\tXC:i:%d", p->clip_len);
putchar('\n'); putchar('\n');
} }
} }
@ -541,8 +541,8 @@ void bwa_print_sam_SQ(const bntseq_t *bns)
{ {
int i; int i;
for (i = 0; i < bns->n_seqs; ++i) for (i = 0; i < bns->n_seqs; ++i)
printf("@SQ\tSN:%s\tLN:%d\n", bns->anns[i].name, bns->anns[i].len); err_printf("@SQ\tSN:%s\tLN:%d\n", bns->anns[i].name, bns->anns[i].len);
if (bwa_rg_line) printf("%s\n", bwa_rg_line); if (bwa_rg_line) err_printf("%s\n", bwa_rg_line);
} }
void bwase_initialize() void bwase_initialize()

View File

@ -176,7 +176,7 @@ void bwa_aln_core(const char *prefix, const char *fn_fa, const gap_opt_t *opt)
} }
// core loop // core loop
fwrite(opt, sizeof(gap_opt_t), 1, stdout); err_fwrite(opt, sizeof(gap_opt_t), 1, stdout);
while ((seqs = bwa_read_seq(ks, 0x40000, &n_seqs, opt->mode, opt->trim_qual)) != 0) { while ((seqs = bwa_read_seq(ks, 0x40000, &n_seqs, opt->mode, opt->trim_qual)) != 0) {
tot_seqs += n_seqs; tot_seqs += n_seqs;
t = clock(); t = clock();
@ -213,8 +213,8 @@ void bwa_aln_core(const char *prefix, const char *fn_fa, const gap_opt_t *opt)
fprintf(stderr, "[bwa_aln_core] write to the disk... "); fprintf(stderr, "[bwa_aln_core] write to the disk... ");
for (i = 0; i < n_seqs; ++i) { for (i = 0; i < n_seqs; ++i) {
bwa_seq_t *p = seqs + i; bwa_seq_t *p = seqs + i;
fwrite(&p->n_aln, 4, 1, stdout); err_fwrite(&p->n_aln, 4, 1, stdout);
if (p->n_aln) fwrite(p->aln, sizeof(bwt_aln1_t), p->n_aln, stdout); if (p->n_aln) err_fwrite(p->aln, sizeof(bwt_aln1_t), p->n_aln, stdout);
} }
fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC); t = clock(); fprintf(stderr, "%.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC); t = clock();

5
main.c
View File

@ -1,9 +1,10 @@
#include <stdio.h> #include <stdio.h>
#include <string.h> #include <string.h>
#include "main.h" #include "main.h"
#include "utils.h"
#ifndef PACKAGE_VERSION #ifndef PACKAGE_VERSION
#define PACKAGE_VERSION "0.5.9-r16" #define PACKAGE_VERSION "0.5.9-r18-dev"
#endif #endif
static int usage() static int usage()
@ -59,5 +60,7 @@ int main(int argc, char *argv[])
fprintf(stderr, "[main] unrecognized command '%s'\n", argv[1]); fprintf(stderr, "[main] unrecognized command '%s'\n", argv[1]);
return 1; return 1;
} }
err_fflush(stdout);
err_fclose(stdout);
return 0; return 0;
} }

66
utils.c
View File

@ -30,6 +30,7 @@
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include <zlib.h> #include <zlib.h>
#include <errno.h>
#include "utils.h" #include "utils.h"
FILE *err_xopen_core(const char *func, const char *fn, const char *mode) FILE *err_xopen_core(const char *func, const char *fn, const char *mode)
@ -80,3 +81,68 @@ void err_fatal_simple_core(const char *func, const char *msg)
fprintf(stderr, "[%s] %s Abort!\n", func, msg); fprintf(stderr, "[%s] %s Abort!\n", func, msg);
abort(); abort();
} }
size_t err_fwrite(const void *ptr, size_t size, size_t nmemb, FILE *stream)
{
size_t ret = fwrite(ptr, size, nmemb, stream);
if (ret != nmemb)
{
err_fatal_simple_core("fwrite", strerror(errno));
}
return ret;
}
int err_printf(const char *format, ...)
{
va_list arg;
int done;
va_start(arg, format);
done = vfprintf(stdout, format, arg);
int saveErrno = errno;
va_end(arg);
if (done < 0)
{
err_fatal_simple_core("vfprintf(stdout)", strerror(saveErrno));
}
return done;
}
int err_fprintf(FILE *stream, const char *format, ...)
{
va_list arg;
int done;
va_start(arg, format);
done = vfprintf(stream, format, arg);
int saveErrno = errno;
va_end(arg);
if (done < 0)
{
err_fatal_simple_core("vfprintf", strerror(saveErrno));
}
return done;
}
int err_fflush(FILE *stream)
{
int ret = fflush(stream);
if (ret != 0)
{
err_fatal_simple_core("fflush", strerror(errno));
}
return ret;
}
int err_fclose(FILE *stream)
{
int ret = fclose(stream);
if (ret != 0)
{
err_fatal_simple_core("fclose", strerror(errno));
}
return ret;
}

16
utils.h
View File

@ -31,6 +31,15 @@
#include <stdio.h> #include <stdio.h>
#include <zlib.h> #include <zlib.h>
#ifdef __GNUC__
// Tell GCC to validate printf format string and args
#define ATTRIBUTE(list) __attribute__ (list)
#else
#define ATTRIBUTE(list)
#endif
#define err_fatal_simple(msg) err_fatal_simple_core(__func__, msg) #define err_fatal_simple(msg) err_fatal_simple_core(__func__, msg)
#define xopen(fn, mode) err_xopen_core(__func__, fn, mode) #define xopen(fn, mode) err_xopen_core(__func__, fn, mode)
#define xreopen(fn, mode, fp) err_xreopen_core(__func__, fn, mode, fp) #define xreopen(fn, mode, fp) err_xreopen_core(__func__, fn, mode, fp)
@ -46,6 +55,13 @@ extern "C" {
FILE *err_xopen_core(const char *func, const char *fn, const char *mode); FILE *err_xopen_core(const char *func, const char *fn, const char *mode);
FILE *err_xreopen_core(const char *func, const char *fn, const char *mode, FILE *fp); FILE *err_xreopen_core(const char *func, const char *fn, const char *mode, FILE *fp);
gzFile err_xzopen_core(const char *func, const char *fn, const char *mode); gzFile err_xzopen_core(const char *func, const char *fn, const char *mode);
size_t err_fwrite(const void *ptr, size_t size, size_t nmemb, FILE *stream);
int err_fprintf(FILE *stream, const char *format, ...)
ATTRIBUTE((format(printf, 2, 3)));
int err_printf(const char *format, ...)
ATTRIBUTE((format(printf, 1, 2)));
int err_fflush(FILE *stream);
int err_fclose(FILE *stream);
#ifdef __cplusplus #ifdef __cplusplus
} }