Skip to content

read FASTQ records with a single copy instead of two #1

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
nh13 wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

read FASTQ records with a single copy instead of two #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
80 changes: 50 additions & 30 deletions bseq.c
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,11 +1,16 @@
#include <zlib.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#define __STDC_LIMIT_MACROS
#include "bseq.h"
#include "kseq.h"
KSEQ_INIT2(, gzFile, gzread)
#include "fr_fastq.h"
#include <zlib.h>

/* 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; }

Expand Down Expand Up @@ -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;
};

Expand All @@ -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);
Expand All @@ -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));
Expand All @@ -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");
}
Expand All @@ -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)
Expand All @@ -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);
}
Loading