From 5a0b32bfd24e89f70460c85bdb03d7dd1d87045a Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 6 Feb 2013 14:38:40 -0500 Subject: [PATCH] updated to the latest kseq.h --- Makefile | 2 +- bntseq.c | 2 +- bseq.c | 4 ++ bseq.h | 15 ++++++ bwamem.c | 30 +++++++++++ bwamem.h | 2 +- bwaseqio.c | 2 +- bwtsw2_aux.c | 2 +- fastmap.c | 2 +- kseq.h | 137 ++++++++++++++++++++++++++++++--------------------- simple_dp.c | 2 +- 11 files changed, 138 insertions(+), 62 deletions(-) create mode 100644 bseq.c create mode 100644 bseq.h diff --git a/Makefile b/Makefile index 04fd7a0..46e0b80 100644 --- a/Makefile +++ b/Makefile @@ -5,7 +5,7 @@ CXXFLAGS= $(CFLAGS) AR= ar DFLAGS= -DHAVE_PTHREAD #-D_NO_SSE2 #-D_FILE_OFFSET_BITS=64 LOBJS= bwa.o bamlite.o utils.o bwt.o bwtio.o bwtaln.o bwtgap.o bntseq.o bwamem.o stdaln.o \ - bwaseqio.o bwase.o kstring.o + bseq.o bwaseqio.o bwase.o kstring.o AOBJS= QSufSort.o bwt_gen.o \ is.o bwtmisc.o bwtindex.o ksw.o simple_dp.o \ bwape.o cs2nt.o \ diff --git a/bntseq.c b/bntseq.c index 18abb2b..06d82a0 100644 --- a/bntseq.c +++ b/bntseq.c @@ -35,7 +35,7 @@ #include "utils.h" #include "kseq.h" -KSEQ_INIT(gzFile, gzread) +KSEQ_DECLARE(gzFile) unsigned char nst_nt4_table[256] = { 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, diff --git a/bseq.c b/bseq.c new file mode 100644 index 0000000..0ec57fa --- /dev/null +++ b/bseq.c @@ -0,0 +1,4 @@ +#include +#include "bseq.h" +#include "kseq.h" +KSEQ_INIT2(, gzFile, gzread) diff --git a/bseq.h b/bseq.h new file mode 100644 index 0000000..73afb63 --- /dev/null +++ b/bseq.h @@ -0,0 +1,15 @@ +#ifndef BATCHSEQ_H_ +#define BATCHSEQ_H_ + +typedef struct { + char *name, *comment, *seq, *qual; +} bseq1_t; + +typedef struct { + int n, m; + bseq1_t *seqs; +} bseq_t; + +int bseq_read(int chunk_size, bseq_t *bs); + +#endif diff --git a/bwamem.c b/bwamem.c index c5c6366..6d08bb8 100644 --- a/bwamem.c +++ b/bwamem.c @@ -274,6 +274,32 @@ int mem_chain_flt(const mem_opt_t *opt, int n_chn, mem_chain1_t *chains) return n; } +#define alnreg_lt(a, b) ((a).score > (b).score) +KSORT_INIT(mem_ar, mem_alnreg_t, alnreg_lt) + +int mem_choose_alnreg_se(const mem_opt_t *opt, int n, mem_alnreg_t *a) +{ // similar to the loop in mem_chain_flt() + int i, j, m; + if (n <= 1) return n; + ks_introsort(mem_ar, n, a); + for (i = 0; i < n; ++i) a[i].sub = 0; + for (i = 1, m = 1; i < n; ++i) { + for (j = 0; j < m; ++j) { + int b_max = a[j].qb > a[i].qb? a[j].qb : a[i].qb; + int e_min = a[j].qe < a[i].qe? a[j].qe : a[i].qe; + if (e_min > b_max) { // have overlap + int min_l = a[i].qe - a[i].qb < a[j].qe - a[j].qb? a[i].qe - a[i].qb : a[j].qe - a[j].qb; + if (e_min - b_max >= min_l * opt->mask_level) { // significant overlap + if (a[j].sub == 0) a[j].sub = a[i].score; + break; + } + } + } + if (j == m) a[m++] = a[i]; + } + return m; +} + /**************************************** * Construct the alignment from a chain * ****************************************/ @@ -388,3 +414,7 @@ ret_gen_cigar: free(rseq); return cigar; } + +/**************** + * Sequence I/O * + ****************/ diff --git a/bwamem.h b/bwamem.h index ef951c3..d69f15a 100644 --- a/bwamem.h +++ b/bwamem.h @@ -31,7 +31,7 @@ typedef struct { typedef struct { int64_t rb, re; - int score, qb, qe; + int score, qb, qe, sub; } mem_alnreg_t; #ifdef __cplusplus diff --git a/bwaseqio.c b/bwaseqio.c index e22d4cd..c1e9f97 100644 --- a/bwaseqio.c +++ b/bwaseqio.c @@ -5,7 +5,7 @@ #include "bamlite.h" #include "kseq.h" -KSEQ_INIT(gzFile, gzread) +KSEQ_DECLARE(gzFile) extern unsigned char nst_nt4_table[256]; static char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 }; diff --git a/bwtsw2_aux.c b/bwtsw2_aux.c index 5e8161c..619930b 100644 --- a/bwtsw2_aux.c +++ b/bwtsw2_aux.c @@ -15,7 +15,7 @@ #include "kstring.h" #include "kseq.h" -KSEQ_INIT(gzFile, gzread) +KSEQ_DECLARE(gzFile) #include "ksort.h" #define __left_lt(a, b) ((a).end > (b).end) diff --git a/fastmap.c b/fastmap.c index c5667e7..475667b 100644 --- a/fastmap.c +++ b/fastmap.c @@ -7,7 +7,7 @@ #include "bwamem.h" #include "kvec.h" #include "kseq.h" -KSEQ_INIT(gzFile, gzread) +KSEQ_DECLARE(gzFile) extern unsigned char nst_nt4_table[256]; diff --git a/kseq.h b/kseq.h index ad8937c..a5cec7c 100644 --- a/kseq.h +++ b/kseq.h @@ -1,6 +1,6 @@ /* The MIT License - Copyright (c) 2008, by Heng Li + Copyright (c) 2008, 2009, 2011 Attractive Chaos Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the @@ -23,6 +23,8 @@ SOFTWARE. */ +/* Last Modified: 05MAR2012 */ + #ifndef AC_KSEQ_H #define AC_KSEQ_H @@ -30,9 +32,14 @@ #include #include +#define KS_SEP_SPACE 0 // isspace(): \t, \n, \v, \f, \r +#define KS_SEP_TAB 1 // isspace() && !' ' +#define KS_SEP_LINE 2 // line separator: "\n" (Unix) or "\r\n" (Windows) +#define KS_SEP_MAX 2 + #define __KS_TYPE(type_t) \ typedef struct __kstream_t { \ - char *buf; \ + unsigned char *buf; \ int begin, end, is_eof; \ type_t f; \ } kstream_t; @@ -45,7 +52,7 @@ { \ kstream_t *ks = (kstream_t*)calloc(1, sizeof(kstream_t)); \ ks->f = f; \ - ks->buf = (char*)malloc(__bufsize); \ + ks->buf = (unsigned char*)malloc(__bufsize); \ return ks; \ } \ static inline void ks_destroy(kstream_t *ks) \ @@ -82,10 +89,10 @@ typedef struct __kstring_t { #endif #define __KS_GETUNTIL(__read, __bufsize) \ - static int ks_getuntil(kstream_t *ks, int delimiter, kstring_t *str, int *dret) \ + static int ks_getuntil2(kstream_t *ks, int delimiter, kstring_t *str, int *dret, int append) \ { \ if (dret) *dret = 0; \ - str->l = 0; \ + str->l = append? str->l : 0; \ if (ks->begin >= ks->end && ks->is_eof) return -1; \ for (;;) { \ int i; \ @@ -97,14 +104,20 @@ typedef struct __kstring_t { if (ks->end == 0) break; \ } else break; \ } \ - if (delimiter) { \ + if (delimiter == KS_SEP_LINE) { \ + for (i = ks->begin; i < ks->end; ++i) \ + if (ks->buf[i] == '\n') break; \ + } else if (delimiter > KS_SEP_MAX) { \ for (i = ks->begin; i < ks->end; ++i) \ if (ks->buf[i] == delimiter) break; \ - } else { \ + } else if (delimiter == KS_SEP_SPACE) { \ for (i = ks->begin; i < ks->end; ++i) \ if (isspace(ks->buf[i])) break; \ - } \ - if (str->m - str->l < i - ks->begin + 1) { \ + } else if (delimiter == KS_SEP_TAB) { \ + for (i = ks->begin; i < ks->end; ++i) \ + if (isspace(ks->buf[i]) && ks->buf[i] != ' ') break; \ + } else i = 0; /* never come to here! */ \ + if (str->m - str->l < (size_t)(i - ks->begin + 1)) { \ str->m = str->l + (i - ks->begin) + 1; \ kroundup32(str->m); \ str->s = (char*)realloc(str->s, str->m); \ @@ -117,9 +130,15 @@ typedef struct __kstring_t { break; \ } \ } \ + if (str->s == 0) { \ + str->m = 1; \ + str->s = (char*)calloc(1, 1); \ + } else if (delimiter == KS_SEP_LINE && str->l > 1 && str->s[str->l-1] == '\r') --str->l; \ str->s[str->l] = '\0'; \ return str->l; \ - } + } \ + static inline int ks_getuntil(kstream_t *ks, int delimiter, kstring_t *str, int *dret) \ + { return ks_getuntil2(ks, delimiter, str, dret, 0); } #define KSTREAM_INIT(type_t, __read, __bufsize) \ __KS_TYPE(type_t) \ @@ -127,19 +146,16 @@ typedef struct __kstring_t { __KS_GETC(__read, __bufsize) \ __KS_GETUNTIL(__read, __bufsize) -#define __KSEQ_BASIC(type_t) \ - static inline kseq_t *kseq_init(type_t fd) \ +#define kseq_rewind(ks) ((ks)->last_char = (ks)->f->is_eof = (ks)->f->begin = (ks)->f->end = 0) + +#define __KSEQ_BASIC(SCOPE, type_t) \ + SCOPE kseq_t *kseq_init(type_t fd) \ { \ kseq_t *s = (kseq_t*)calloc(1, sizeof(kseq_t)); \ s->f = ks_init(fd); \ return s; \ } \ - static inline void kseq_rewind(kseq_t *ks) \ - { \ - ks->last_char = 0; \ - ks->f->is_eof = ks->f->begin = ks->f->end = 0; \ - } \ - static inline void kseq_destroy(kseq_t *ks) \ + SCOPE void kseq_destroy(kseq_t *ks) \ { \ if (!ks) return; \ free(ks->name.s); free(ks->comment.s); free(ks->seq.s); free(ks->qual.s); \ @@ -152,44 +168,46 @@ typedef struct __kstring_t { -1 end-of-file -2 truncated quality string */ -#define __KSEQ_READ \ - static int kseq_read(kseq_t *seq) \ - { \ - int c; \ - kstream_t *ks = seq->f; \ +#define __KSEQ_READ(SCOPE) \ + SCOPE int kseq_read(kseq_t *seq) \ + { \ + int c; \ + kstream_t *ks = seq->f; \ if (seq->last_char == 0) { /* then jump to the next header line */ \ - while ((c = ks_getc(ks)) != -1 && c != '>' && c != '@'); \ - if (c == -1) return -1; /* end of file */ \ - seq->last_char = c; \ - } /* the first header char has been read */ \ - seq->comment.l = seq->seq.l = seq->qual.l = 0; \ - if (ks_getuntil(ks, 0, &seq->name, &c) < 0) return -1; \ - if (c != '\n') ks_getuntil(ks, '\n', &seq->comment, 0); \ + while ((c = ks_getc(ks)) != -1 && c != '>' && c != '@'); \ + if (c == -1) return -1; /* end of file */ \ + seq->last_char = c; \ + } /* else: the first header char has been read in the previous call */ \ + seq->comment.l = seq->seq.l = seq->qual.l = 0; /* reset all members */ \ + if (ks_getuntil(ks, 0, &seq->name, &c) < 0) return -1; /* normal exit: EOF */ \ + if (c != '\n') ks_getuntil(ks, KS_SEP_LINE, &seq->comment, 0); /* read FASTA/Q comment */ \ + if (seq->seq.s == 0) { /* we can do this in the loop below, but that is slower */ \ + seq->seq.m = 256; \ + seq->seq.s = (char*)malloc(seq->seq.m); \ + } \ while ((c = ks_getc(ks)) != -1 && c != '>' && c != '+' && c != '@') { \ - if (isgraph(c)) { /* printable non-space character */ \ - if (seq->seq.l + 1 >= seq->seq.m) { /* double the memory */ \ - seq->seq.m = seq->seq.l + 2; \ - kroundup32(seq->seq.m); /* rounded to next closest 2^k */ \ - seq->seq.s = (char*)realloc(seq->seq.s, seq->seq.m); \ - } \ - seq->seq.s[seq->seq.l++] = (char)c; \ - } \ - } \ + if (c == '\n') continue; /* skip empty lines */ \ + seq->seq.s[seq->seq.l++] = c; /* this is safe: we always have enough space for 1 char */ \ + ks_getuntil2(ks, KS_SEP_LINE, &seq->seq, 0, 1); /* read the rest of the line */ \ + } \ if (c == '>' || c == '@') seq->last_char = c; /* the first header char has been read */ \ - seq->seq.s[seq->seq.l] = 0; /* null terminated string */ \ - if (c != '+') return seq->seq.l; /* FASTA */ \ - if (seq->qual.m < seq->seq.m) { /* allocate enough memory */ \ - seq->qual.m = seq->seq.m; \ - seq->qual.s = (char*)realloc(seq->qual.s, seq->qual.m); \ - } \ + if (seq->seq.l + 1 >= seq->seq.m) { /* seq->seq.s[seq->seq.l] below may be out of boundary */ \ + seq->seq.m = seq->seq.l + 2; \ + kroundup32(seq->seq.m); /* rounded to the next closest 2^k */ \ + seq->seq.s = (char*)realloc(seq->seq.s, seq->seq.m); \ + } \ + seq->seq.s[seq->seq.l] = 0; /* null terminated string */ \ + if (c != '+') return seq->seq.l; /* FASTA */ \ + if (seq->qual.m < seq->seq.m) { /* allocate memory for qual in case insufficient */ \ + seq->qual.m = seq->seq.m; \ + seq->qual.s = (char*)realloc(seq->qual.s, seq->qual.m); \ + } \ while ((c = ks_getc(ks)) != -1 && c != '\n'); /* skip the rest of '+' line */ \ - if (c == -1) return -2; /* we should not stop here */ \ - while ((c = ks_getc(ks)) != -1 && seq->qual.l < seq->seq.l) \ - if (c >= 33 && c <= 127) seq->qual.s[seq->qual.l++] = (unsigned char)c; \ - seq->qual.s[seq->qual.l] = 0; /* null terminated string */ \ + if (c == -1) return -2; /* error: no quality string */ \ + while (ks_getuntil2(ks, KS_SEP_LINE, &seq->qual, 0, 1) >= 0 && seq->qual.l < seq->seq.l); \ seq->last_char = 0; /* we have not come to the next header line */ \ - if (seq->seq.l != seq->qual.l) return -2; /* qual string is shorter than seq string */ \ - return seq->seq.l; \ + if (seq->seq.l != seq->qual.l) return -2; /* error: qual string is of a different length */ \ + return seq->seq.l; \ } #define __KSEQ_TYPE(type_t) \ @@ -199,10 +217,19 @@ typedef struct __kstring_t { kstream_t *f; \ } kseq_t; -#define KSEQ_INIT(type_t, __read) \ - KSTREAM_INIT(type_t, __read, 4096) \ +#define KSEQ_INIT2(SCOPE, type_t, __read) \ + KSTREAM_INIT(type_t, __read, 16384) \ __KSEQ_TYPE(type_t) \ - __KSEQ_BASIC(type_t) \ - __KSEQ_READ + __KSEQ_BASIC(SCOPE, type_t) \ + __KSEQ_READ(SCOPE) + +#define KSEQ_INIT(type_t, __read) KSEQ_INIT2(static, type_t, __read) + +#define KSEQ_DECLARE(type_t) \ + __KS_TYPE(type_t) \ + __KSEQ_TYPE(type_t) \ + extern kseq_t *kseq_init(type_t fd); \ + void kseq_destroy(kseq_t *ks); \ + int kseq_read(kseq_t *seq); #endif diff --git a/simple_dp.c b/simple_dp.c index 7c078c2..d2b4b71 100644 --- a/simple_dp.c +++ b/simple_dp.c @@ -8,7 +8,7 @@ #include "utils.h" #include "kseq.h" -KSEQ_INIT(gzFile, gzread) +KSEQ_DECLARE(gzFile) typedef struct { int l;