diff --git a/Makefile b/Makefile index 017a6ea..430ff4a 100644 --- a/Makefile +++ b/Makefile @@ -5,7 +5,7 @@ LDFLAGS= INCLUDES= LOBJS= kommon.o kalloc.o bwt.o l2bit.o options.o seed.o map-algo.o lchain.o align.o pe.o cs.o format.o \ ksw2_extz2_sse.o ksw2_extd2_sse.o ksw2_ll_sse.o -AOBJS= kthread.o libsais.o libsais64.o index.o bseq.o map-main.o fastmap.o +AOBJS= kthread.o libsais.o libsais64.o index.o bseq.o fr_fastq.o map-main.o fastmap.o MALLOC_O= mimalloc.o PROG= minibwa LIBS= -lpthread -lz -lm @@ -65,7 +65,8 @@ depend: QSufSort.o: QSufSort.h align.o: mbpriv.h minibwa.h l2bit.h bwt.h kommon.h bseq.h kalloc.h ksw2.h -bseq.o: bseq.h kommon.h kseq.h +bseq.o: bseq.h fr_fastq.h kommon.h +fr_fastq.o: fr_fastq.h bwt.o: kommon.h kalloc.h bwt.h bwtgen.o: QSufSort.h cs.o: mbpriv.h minibwa.h l2bit.h bwt.h kommon.h bseq.h kalloc.h diff --git a/bseq.c b/bseq.c index c3630d7..0511585 100644 --- a/bseq.c +++ b/bseq.c @@ -1,11 +1,16 @@ -#include #include #include +#include #include #define __STDC_LIMIT_MACROS #include "bseq.h" -#include "kseq.h" -KSEQ_INIT2(, gzFile, gzread) +#include "fr_fastq.h" +#include + +/* Byte source for the fr_fastq parser: stock zlib gzread. gzread transparently + * passes plaintext through when the input is not gzip, so the same path serves + * an in-process .gz file and a decompressed plaintext pipe on stdin. */ +static int mb_src_read(void *ctx, unsigned char *buf, int len) { return gzread((gzFile)ctx, buf, len); } #define kvec_t(type) struct { size_t n, m; type *a; } @@ -36,7 +41,7 @@ KSEQ_INIT2(, gzFile, gzread) struct mb_bseq_file_s { gzFile fp; - kseq_t *ks; + fr_fastq_t *parser; mb_bseq1_t s; }; @@ -48,46 +53,54 @@ mb_bseq_file_t *mb_bseq_open(const char *fn) if (f == 0) return 0; fp = (mb_bseq_file_t*)calloc(1, sizeof(mb_bseq_file_t)); fp->fp = f; - fp->ks = kseq_init(fp->fp); + fp->parser = fr_fastq_init(mb_src_read, fp->fp); return fp; } void mb_bseq_close(mb_bseq_file_t *fp) { - kseq_destroy(fp->ks); + fr_fastq_destroy(fp->parser); gzclose(fp->fp); + free(fp->s.name); free(fp->s.seq); free(fp->s.qual); free(fp->s.comment); free(fp); } -static inline char *kstrdup(const kstring_t *s) +/* Length-aware dup: malloc len+1 and NUL-terminate (fr_fastq slices are not + * NUL-terminated). Mirrors the storage kstrdup produced from a kstring_t. */ +static inline char *fr_dup(const char *s, size_t l) { - char *t; - t = (char*)malloc(s->l + 1); - memcpy(t, s->s, s->l + 1); + char *t = (char*)malloc(l + 1); + memcpy(t, s, l); + t[l] = 0; return t; } -static inline void kseq2bseq(kseq_t *ks, mb_bseq1_t *s, int with_qual, int with_comment) +/* Copy one parsed record into mb_bseq1_t. Byte-for-byte the same as the old + * kseq2bseq: warn on empty name, own all strings, convert U->T via the same + * `--c` decrement ('U'->'T', 'u'->'t'), and store comment/qual only when + * present and requested. */ +static inline void fr_rec_to_bseq1(const fr_fastq_rec_t *r, mb_bseq1_t *s, int with_qual, int with_comment) { int i; - if (ks->name.l == 0) + if (r->name_l == 0) fprintf(stderr, "[WARNING]\033[1;31m empty sequence name in the input.\033[0m\n"); - s->name = kstrdup(&ks->name); - s->seq = kstrdup(&ks->seq); - for (i = 0; i < (int)ks->seq.l; ++i) // convert U to T + s->name = fr_dup(r->name, r->name_l); + s->seq = fr_dup(r->seq, r->seq_l); + for (i = 0; i < (int)r->seq_l; ++i) // convert U to T if (s->seq[i] == 'u' || s->seq[i] == 'U') --s->seq[i]; - s->qual = with_qual && ks->qual.l? kstrdup(&ks->qual) : 0; - s->comment = with_comment && ks->comment.l? kstrdup(&ks->comment) : 0; - s->l_seq = ks->seq.l; + s->qual = with_qual && r->qual_l? fr_dup(r->qual, r->qual_l) : 0; + s->comment = with_comment && r->comment_l? fr_dup(r->comment, r->comment_l) : 0; + s->l_seq = r->seq_l; } mb_bseq1_t *mb_bseq_read(mb_bseq_file_t *fp, int64_t chunk_size, int with_qual, int with_comment, int frag_mode, int min_cnt, int64_t max_chunk_size, int *n_) { int64_t size = 0; int ret; + fr_fastq_rec_t rec; kvec_t(mb_bseq1_t) a = {0,0,0}; - kseq_t *ks = fp->ks; + fr_fastq_t *parser = fp->parser; *n_ = 0; if (fp->s.seq) { kv_resize(mb_bseq1_t, a, 256); @@ -97,22 +110,24 @@ mb_bseq1_t *mb_bseq_read(mb_bseq_file_t *fp, int64_t chunk_size, int with_qual, } if (max_chunk_size < chunk_size) max_chunk_size = chunk_size; - while ((ret = kseq_read(ks)) >= 0) { + for (;;) { + ret = fr_fastq_next(parser, &rec); + if (ret <= 0) break; int32_t to_stop = 0; mb_bseq1_t *s; - assert(ks->seq.l <= INT32_MAX); + assert(rec.seq_l <= INT32_MAX); if (a.m == 0) kv_resize(mb_bseq1_t, a, 256); kv_pushp(mb_bseq1_t, a, &s); - kseq2bseq(ks, s, with_qual, with_comment); - size += ks->seq.l; + fr_rec_to_bseq1(&rec, s, with_qual, with_comment); + size += rec.seq_l; if (chunk_size <= 0 || max_chunk_size <= 0) to_stop = 1; else if (size >= max_chunk_size) to_stop = 1; else if (size >= chunk_size && a.n >= min_cnt) to_stop = 1; if (to_stop) { if (frag_mode && a.a[a.n-1].l_seq < CHECK_PAIR_THRES) { - while ((ret = kseq_read(ks)) >= 0) { - kseq2bseq(ks, &fp->s, with_qual, with_comment); - size += ks->seq.l; + while ((ret = fr_fastq_next(parser, &rec)) == 1) { + fr_rec_to_bseq1(&rec, &fp->s, with_qual, with_comment); + size += rec.seq_l; if (mb_qname_same(fp->s.name, a.a[a.n-1].name)) { kv_push(mb_bseq1_t, a, fp->s); memset(&fp->s, 0, sizeof(mb_bseq1_t)); @@ -122,7 +137,7 @@ mb_bseq1_t *mb_bseq_read(mb_bseq_file_t *fp, int64_t chunk_size, int with_qual, break; } } - if (ret < -1) { + if (ret == -2) { if (a.n) fprintf(stderr, "[WARNING]\033[1;31m failed to parse the FASTA/FASTQ record next to '%s'. Continue anyway.\033[0m\n", a.a[a.n-1].name); else fprintf(stderr, "[WARNING]\033[1;31m failed to parse the first FASTA/FASTQ record. Continue anyway.\033[0m\n"); } @@ -134,13 +149,17 @@ mb_bseq1_t *mb_bseq_read_frag(int n_fp, mb_bseq_file_t **fp, int64_t chunk_size, { int i; int64_t size = 0; + fr_fastq_rec_t *rec; kvec_t(mb_bseq1_t) a = {0,0,0}; *n_ = 0; if (n_fp < 1) return 0; + /* One record slot per file: each file has its own parser, so the slices + * stay valid across the per-file reads until we copy them out below. */ + rec = (fr_fastq_rec_t*)calloc(n_fp, sizeof(fr_fastq_rec_t)); while (1) { int n_read = 0; for (i = 0; i < n_fp; ++i) - if (kseq_read(fp[i]->ks) >= 0) + if (fr_fastq_next(fp[i]->parser, &rec[i]) == 1) ++n_read; if (n_read < n_fp) { if (n_read > 0) @@ -151,16 +170,17 @@ mb_bseq1_t *mb_bseq_read_frag(int n_fp, mb_bseq_file_t **fp, int64_t chunk_size, for (i = 0; i < n_fp; ++i) { mb_bseq1_t *s; kv_pushp(mb_bseq1_t, a, &s); - kseq2bseq(fp[i]->ks, s, with_qual, with_comment); + fr_rec_to_bseq1(&rec[i], s, with_qual, with_comment); size += s->l_seq; } if (size >= chunk_size) break; } + free(rec); *n_ = a.n; return a.a; } int mb_bseq_eof(mb_bseq_file_t *fp) { - return (ks_eof(fp->ks->f) && fp->s.seq == 0); + return (fr_fastq_eof(fp->parser) && fp->s.seq == 0); } diff --git a/fr_fastq.c b/fr_fastq.c new file mode 100644 index 0000000..2ccc4d6 --- /dev/null +++ b/fr_fastq.c @@ -0,0 +1,255 @@ +/* fr_fastq: single-copy FASTQ/FASTA parser over a pluggable byte source. + * + * The scan_record() grammar mirrors kseq's exactly — that is what makes the + * output byte-identical to the kseq path. fr_refill() pulls bytes through a + * caller-supplied callback (p->read_fn(p->ctx, ...)), so the parser is not tied + * to any one byte source. + * + * Strategy: keep a large refillable buffer and guarantee the *whole* current + * record is resident before handing back field slices. For the dominant case + * (single-line FASTQ/FASTA) every field is then a contiguous run in that buffer, + * so the caller copies each field exactly once (no kstring_t intermediate). + * Multi-line seq/qual are the only case that needs concatenation; those are + * compacted into a small scratch buffer, matching what kseq did anyway. + */ +#include "fr_fastq.h" + +#include +#include +#include + +#define FR_FASTQ_BUF_INIT (1u << 20) /* 1 MiB; a record almost always fits */ +#define FR_FASTQ_READ_MAX (1u << 20) /* cap a single byte-source read */ + +/* scan_record return sentinel: the record is not yet fully buffered. */ +#define FR_UNDERRUN (-3) + +struct fr_fastq { + fr_read_fn read_fn; /* borrowed byte source */ + void *ctx; + unsigned char *buf; /* refillable read buffer */ + size_t cap, beg, end; /* capacity, parse cursor, valid byte count */ + int eof; /* byte source returned 0 */ + int err; /* byte source returned -1 (decode error) */ + /* scratch for assembling multi-line seq/qual (allocation persists across + * records; the used length is recomputed per record). */ + unsigned char *seq_scr; size_t seq_scr_cap; + unsigned char *qual_scr; size_t qual_scr_cap; +}; + +static void *xrealloc(void *p, size_t n) +{ + void *q = realloc(p, n); + if (q == NULL) { abort(); } /* OOM on the read path is unrecoverable */ + return q; +} + +fr_fastq_t *fr_fastq_init(fr_read_fn read_fn, void *ctx) +{ + fr_fastq_t *p = (fr_fastq_t *)calloc(1, sizeof(*p)); + if (p == NULL) abort(); + p->read_fn = read_fn; + p->ctx = ctx; + p->cap = FR_FASTQ_BUF_INIT; + p->buf = (unsigned char *)xrealloc(NULL, p->cap); + return p; +} + +void fr_fastq_destroy(fr_fastq_t *p) +{ + if (!p) return; + free(p->buf); free(p->seq_scr); free(p->qual_scr); + free(p); +} + +int fr_fastq_eof(const fr_fastq_t *p) +{ + return p && p->eof && p->beg >= p->end; +} + +/* Pull more bytes in: drop the consumed prefix [0, beg), then read into the + * tail, growing the buffer if a single record fills it. Returns bytes read + * (0 at EOF, -1 on decode error). */ +static int fr_refill(fr_fastq_t *p) +{ + if (p->beg > 0) { + memmove(p->buf, p->buf + p->beg, p->end - p->beg); + p->end -= p->beg; + p->beg = 0; + } + if (p->end == p->cap) { /* record larger than the buffer: grow */ + p->cap *= 2; + p->buf = (unsigned char *)xrealloc(p->buf, p->cap); + } + size_t want = p->cap - p->end; + if (want > FR_FASTQ_READ_MAX) want = FR_FASTQ_READ_MAX; + int got = p->read_fn(p->ctx, p->buf + p->end, (int)want); + if (got < 0) { p->err = 1; return -1; } + if (got == 0) p->eof = 1; + else p->end += (size_t)got; + return got; +} + +/* Append n bytes to a scratch buffer, growing as needed. */ +static void scr_put(unsigned char **bufp, size_t *capp, size_t *lenp, + const unsigned char *src, size_t n) +{ + if (*lenp + n > *capp) { + size_t ncap = *capp ? *capp : 256; + while (ncap < *lenp + n) ncap *= 2; + *bufp = (unsigned char *)xrealloc(*bufp, ncap); + *capp = ncap; + } + memcpy(*bufp + *lenp, src, n); + *lenp += n; +} + +/* Length of a line's content with a single trailing '\r' dropped, matching + * kseq's KS_SEP_LINE rule: strip only when the copied run (which includes the + * '\r') is longer than one byte. */ +static size_t line_content_len(const unsigned char *s, size_t raw_len) +{ + if (raw_len > 1 && s[raw_len - 1] == '\r') return raw_len - 1; + return raw_len; +} + +/* Index of the next '\n' in b[from..end), or `end` if none. memchr is the + * libc-vectorised newline scan (NEON/SSE/AVX) — the sequence and quality lines + * are the dominant scan cost, so this is where the speedup lives. */ +static inline size_t find_nl(const unsigned char *b, size_t from, size_t end) +{ + const unsigned char *nl = (const unsigned char *)memchr(b + from, '\n', end - from); + return nl ? (size_t)(nl - b) : end; +} + +/* Attempt to parse one record entirely out of the resident buffer. + * 1 -> *rec filled, *new_beg = byte after the record + * 0 -> clean EOF + * -2 -> malformed (matches kseq_read's -2) + * FR_UNDERRUN -> the record is not fully buffered yet (caller refills) */ +static int scan_record(fr_fastq_t *p, fr_fastq_rec_t *rec, size_t *new_beg) +{ + const unsigned char *b = p->buf; + const size_t end = p->end; + + /* Phase 1: skip to a header line. Junk/blank lines before it are consumed + * for good, so advance beg past them. */ + while (p->beg < end && b[p->beg] != '>' && b[p->beg] != '@') p->beg++; + if (p->beg >= end) return p->eof ? 0 : FR_UNDERRUN; + + const int is_fastq_marker = (b[p->beg] == '@'); + size_t i = p->beg + 1; /* first byte after the header char */ + + /* name: up to the first whitespace (kseq KS_SEP_SPACE). */ + size_t name_s = i; + while (i < end && !isspace(b[i])) i++; + if (i >= end && !p->eof) return FR_UNDERRUN; + size_t name_e = i; + + /* comment: present iff the name stopped at a non-newline whitespace. */ + size_t com_s = name_e, com_e = name_e; + size_t after_header; /* first byte of the line below the header */ + if (i < end && b[i] != '\n') { + size_t cs = i + 1, j = find_nl(b, cs, end); + if (j >= end && !p->eof) return FR_UNDERRUN; + com_s = cs; + com_e = cs + line_content_len(b + cs, j - cs); + after_header = (j < end) ? j + 1 : end; + } else { + after_header = (i < end) ? i + 1 : end; + } + + /* seq: sequence lines (newlines stripped) until the next header or '+'. */ + size_t k = after_header; + size_t seq_s = after_header, seq_e = after_header; /* single-line fast path */ + size_t seq_len = 0; int seq_lines = 0; + unsigned char term = 0; /* 0 = EOF, else '>' '@' '+' */ + for (;;) { + if (k >= end) { if (p->eof) break; else return FR_UNDERRUN; } + unsigned char c = b[k]; + if (c == '>' || c == '@' || c == '+') { term = c; break; } + if (c == '\n') { k++; continue; } /* skip empty lines */ + size_t ls = k; + k = find_nl(b, k, end); + if (k >= end && !p->eof) return FR_UNDERRUN; + size_t clen = line_content_len(b + ls, k - ls); + if (seq_lines == 0) { seq_s = ls; seq_e = ls + clen; } + else { + if (seq_lines == 1) { /* promote first line to scratch */ + size_t l = 0; + scr_put(&p->seq_scr, &p->seq_scr_cap, &l, b + seq_s, seq_e - seq_s); + seq_len = l; + } + scr_put(&p->seq_scr, &p->seq_scr_cap, &seq_len, b + ls, clen); + } + if (seq_lines == 0) seq_len = clen; + seq_lines++; + if (k < end) k++; /* step past '\n' */ + } + + rec->name = (const char *)(b + name_s); + rec->name_l = name_e - name_s; + rec->comment = (const char *)(b + com_s); + rec->comment_l = com_e - com_s; + if (seq_lines <= 1) { rec->seq = (const char *)(b + seq_s); } + else { rec->seq = (const char *)p->seq_scr; } + rec->seq_l = seq_len; + + if (term != '+') { /* FASTA (or EOF): no quality */ + (void)is_fastq_marker; + rec->qual = NULL; rec->qual_l = 0; + *new_beg = k; /* leave the next header unconsumed */ + return 1; + } + + /* skip the '+' separator line. */ + size_t pl = find_nl(b, k, end); + if (pl >= end && !p->eof) return FR_UNDERRUN; + if (pl >= end) return -2; /* '+' line with no newline and no quality */ + size_t m = pl + 1; + + /* qual: quality lines accumulated until they reach the sequence length + * (kseq reads at least one line, then stops once qual.l >= seq.l). */ + size_t qual_s = m, qual_e = m, qual_len = 0; int qual_lines = 0; + for (;;) { + if (m >= end) { if (p->eof) break; else return FR_UNDERRUN; } + size_t ls = m; + m = find_nl(b, m, end); + if (m >= end && !p->eof) return FR_UNDERRUN; + size_t clen = line_content_len(b + ls, m - ls); + if (qual_lines == 0) { qual_s = ls; qual_e = ls + clen; } + else { + if (qual_lines == 1) { + size_t l = 0; + scr_put(&p->qual_scr, &p->qual_scr_cap, &l, b + qual_s, qual_e - qual_s); + qual_len = l; + } + scr_put(&p->qual_scr, &p->qual_scr_cap, &qual_len, b + ls, clen); + } + if (qual_lines == 0) qual_len = clen; + qual_lines++; + if (m < end) m++; /* step past '\n' */ + if (qual_len >= seq_len) break; + } + if (qual_len != seq_len) return -2; /* kseq: seq.l != qual.l */ + + if (qual_lines <= 1) { rec->qual = (const char *)(b + qual_s); } + else { rec->qual = (const char *)p->qual_scr; } + rec->qual_l = qual_len; + *new_beg = m; + return 1; +} + +int fr_fastq_next(fr_fastq_t *p, fr_fastq_rec_t *rec) +{ + for (;;) { + size_t new_beg = p->beg; + int r = scan_record(p, rec, &new_beg); + if (r != FR_UNDERRUN) { + if (r == 1) p->beg = new_beg; + return r; + } + if (fr_refill(p) < 0) return -2; /* decode error: stop the batch */ + /* on EOF, scan_record will now resolve (no more underruns possible). */ + } +} diff --git a/fr_fastq.h b/fr_fastq.h new file mode 100644 index 0000000..316f6d7 --- /dev/null +++ b/fr_fastq.h @@ -0,0 +1,75 @@ +/* fr_fastq: a single-copy FASTQ/FASTA parser over a pluggable byte source. + * + * Replaces the kseq layer on the read hot path. kseq copies every field twice + * (byte source -> kstring_t, then kstring_t -> mb_bseq1_t); this parser slices + * records directly out of its own refillable buffer so the caller copies each + * field exactly once. + * + * The parser is not tied to any one decompressor: it pulls bytes through a + * caller-supplied callback (fr_read_fn) + opaque context. bseq.c feeds it + * gzread over a gzFile, which transparently handles both gzip and plaintext, so + * the same parser serves an in-process .gz file and a decompressed plaintext + * pipe on stdin. + * + * The grammar matches kseq's exactly so output is byte-identical: + * - record starts at '>' (FASTA) or '@' (FASTQ); leading junk before the + * first header is skipped, as are blank lines between records; + * - name = bytes up to the first whitespace; comment = the rest of that line + * (empty if the name ran to end-of-line), with a single trailing '\r' + * dropped when the line has more than one byte; + * - seq = the sequence line(s), newlines removed, up to the next header or + * the '+' separator; CRLF endings are normalised the same way kseq does; + * - qual (FASTQ only) = quality line(s) accumulated until they match the + * sequence length. + * + * Slice pointers returned by fr_fastq_next() reference parser-internal storage + * and are valid only until the next fr_fastq_next() / fr_fastq_destroy() call. + */ +#ifndef FR_FASTQ_H +#define FR_FASTQ_H + +#include + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct fr_fastq fr_fastq_t; + +/* Byte-source callback: read up to `len` bytes into `buf` from `ctx`. Returns + * the number of bytes read (0 at clean EOF, -1 on a decode/read error). This is + * the exact contract of zlib's gzread. */ +typedef int (*fr_read_fn)(void *ctx, unsigned char *buf, int len); + +/* One parsed record. Lengths are exact; pointers are NOT guaranteed to be + * NUL-terminated. comment_l == 0 means "no comment" (mirrors kseq's + * comment.l == 0). qual_l == 0 means no quality (FASTA, or an empty record); + * downstream should treat comment/qual as absent when their length is 0. */ +typedef struct { + const char *name; size_t name_l; + const char *comment; size_t comment_l; + const char *seq; size_t seq_l; + const char *qual; size_t qual_l; +} fr_fastq_rec_t; + +/* Bind a parser to a byte source. `ctx` is borrowed: it must outlive the parser + * and is NOT closed by fr_fastq_destroy. Aborts on OOM. */ +fr_fastq_t *fr_fastq_init(fr_read_fn read_fn, void *ctx); + +/* Free the parser's buffers. `p` may be NULL. Does not touch the byte source. */ +void fr_fastq_destroy(fr_fastq_t *p); + +/* Parse the next record into *rec. + * returns 1 -> a record was parsed; slices in *rec are valid until the next call + * returns 0 -> clean end of input + * returns -2 -> malformed record (matches kseq_read's -2). */ +int fr_fastq_next(fr_fastq_t *p, fr_fastq_rec_t *rec); + +/* 1 once the byte source has signalled EOF and every buffered byte is consumed. */ +int fr_fastq_eof(const fr_fastq_t *p); + +#ifdef __cplusplus +} +#endif + +#endif /* FR_FASTQ_H */