Fast way to count number of reads and number of bases in a fastq file? I am looking for a tool, preferably written in C or C++, that can quickly and efficiently count the number of reads and the number of bases in a compressed fastq file. I am currently doing this using zgrep and awk:

zgrep . foo.fasq.gz |
     awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c; 
                print "Number of bases in reads: "l
              }'


The zgrep . will print non-blank lines from the input file and the awk 'NR%4==2 will process every 4th line starting with the second (the sequence). 
This works fine, but can take a very long time when dealing with large files such as WGS data. Is there a tool I can use (on Linux) that will give me these values? Or, if not, I'm also open to suggestions for speeding up the above command.
It's difficult to get this to go massively quicker I think - as with this question working with large gzipped FASTQ files is mostly IO-bound. We could instead focus on making sure we are getting the right answer.

People deride them too often, but this is where a well-written parser is worth it's weight in gold. Heng Li gives us this FASTQ Parser in C. 

I downloaded the example tarball and modified the example code (excuse my C...):

#include &lt;zlib.h&gt;
#include &lt;stdio.h&gt;
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
    gzFile fp;
    kseq_t *seq;
    int l;
    if (argc == 1) {
        fprintf(stderr, "Usage: %s &lt;in.seq&gt;\n", argv[0]);
        return 1;
    }
    fp = gzopen(argv[1], "r");
    seq = kseq_init(fp);
    int seqcount = 0;
    long seqlen = 0;
    while ((l = kseq_read(seq)) &gt;= 0) {
        seqcount = seqcount + 1;
        seqlen = seqlen + (long)strlen(seq-&gt;seq.s);
    }
    kseq_destroy(seq);
    gzclose(fp);
    printf("Number of sequences: %d\n", seqcount);
    printf("Number of bases in sequences: %ld\n", seqlen);
    return 0;
}


Then make and kseq_test foo.fastq.gz.

For my example file (~35m reads of ~75bp) this took:

real    0m49.670s
user    0m49.364s
sys     0m0.304s


Compared with your example:

real    0m43.616s
user    1m35.060s
sys     0m5.240s


Konrad's solution (in my hands):

real    0m39.682s
user    1m11.900s
sys     0m5.112s


(By the way, just zcat-ing the data file to /dev/null):

real    0m38.736s
user    0m38.356s
sys     0m0.308s


So, I get pretty close in speed, but am likely to be more standards compliant. Also this solution gives you more flexibility with what you can do with the data.

And my horrible C can almost certainly be optimised.



Same test, with kseq.h from Github, as suggested in the comments:

My machine is under different load this morning, so I've retested. Wall clock times:

OP: 0m44.813s

Konrad: 0m40.061s

zcat &gt; /dev/null: 0m34.508s

kseq.h (Github): 0m32.909s

So most recent version of kseq.h is faster than simply zcat-ing the file (consistently in my tests...).