pigz | awk | wc is the fastest method
First off for benchmarks with FASTQ it's best to use a specific real-world example with a known answer. I've chosen this file:
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01815/sequence_read/ERR047740_1.filt.fastq.gz
as my test file, the correct answers being:
Number of reads: 67051220
Number of bases in reads: 6034609800
Next we want to find the fastest way possible to count these, all timings are the average wall-clock time (real) of 10 runs collected with the bash time on an otherwise unloaded system:
zgrep
zgrep . ERR047740_1.filt.fastq.gz |
awk 'NR%4==2{c++; l+=length($0)}
END{
print "Number of reads: "c;
print "Number of bases in reads: "l
}'
This is the slowest method with an average run-time of 125.35 seconds
gzip awk
Using gzip we gain about another 10 seconds:
gzip -dc ERR047740_1.filt.fastq.gz |
awk 'NR%4==2{c++; l+=length($0)}
END{
print "Number of reads: "c;
print "Number of bases in reads: "l
}'
Average run-time is 116.69 seconds
Konrad's gzip awk wc variant
fix_base_count() {
local counts=($(cat))
echo "${counts[0]} $((${counts[1]} - ${counts[0]}))"
}
gzip -dc ERR047740_1.filt.fastq.gz \
| awk 'NR % 4 == 2' \
| wc -cl \
| fix_base_count
This runs slower on this test file than the gzip awk variant of the solution, average run-time is 122.28 seconds.
kseq_test using latest kseq.h from klib
Code compiled with: gcc -O2 -o kseq_test kseq_test.c -lz where kseq_test.c is Simon's adaptation of Heng Li's FASTQ parser.
kseq_test ERR047740_1.filt.fastq.gz
Average run-time is 99.14 seconds, which is better than the gzip core utilities based solution so far, but we can do better!
piz awk
Using Mark Adler's pigz as a drop-in replacement for gzip, note that pigz gives us a speed gain as on top of gzip as in addition to the main deflate thread it uses another 3 threads for reading, writing and checksum calculations, see the man page for details.
pigz -dc ERR047740_1.filt.fastq.gz |
awk 'NR%4==2{c++; l+=length($0)}
END{
print "Number of reads: "c;
print "Number of bases in reads: "l
}'
Average run-time is now 93.86 seconds, this is ~5 seconds faster than the kseq based C code but we can further improve the benchmark.
pigz awk wc
Next we use pigz as a drop in replacment for Konrad's wc variant of the awk based solution.
fix_base_count() {
local counts=($(cat))
echo "${counts[0]} $((${counts[1]} - ${counts[0]}))"
}
gzip -dc ERR047740_1.filt.fastq.gz \
| awk 'NR % 4 == 2' \
| wc -cl \
| fix_base_count
Average run-time is now down to 83.03 seconds, this is ~16 seconds faster than the kseq based solution and ~42 seconds faster than the OPs zgrep based solution.
Next as a baseline lets see just how much of this run-time is due to decompression of the input fastq.gz file.
gzip alone
gzip -dc ERR047740_1.filt.fastq.gz > /dev/null
Average run-time: 105.95 seconds, so the gzip based solutions (which also includes zcat and zgrep as these are provided by gzip) are never going to be faster than kseq_test.
pigz alone
pigz -dc ERR047740_1.filt.fastq.gz > /dev/null
Average run-time: 77.66 seconds, so quite clearly the additional three threads for read, write and checksum calculation offer a useful advantage. What's more this speed-up is greater when leveraging the awk | wc based solution, it's not clear why, but I expect this is due to the extra write thread.
Interestingly average CPU usage across all threads is quite revealing for the various answers, I've collated these stats using GNU time /usr/bin/time --verbose
zgrep based solution 133% - must be more than one thread somehow
gzip | awk based solution 99% - all gzip based solutions run single-threaded at 99% CPU usage
pigz | awk 147%
gzip | awk | wc 99% as with gzip
pgiz | awk | wc 155%
kseq_test 99%
gzip > dev/null 99%
pigz > dev/null 155%
Whilst the main deflate thread in pigz will run at 100% CPU load the extra 3 don't quite fully occupy additional cores to 100% (as is evidenced by average CPU usage of ~150%) they do however clearly result in reduced run-time.
I'm using Ubuntu 16.04.2 LTS**, my gzip, zcat, zgrep versions are all gzip 1.6 and pigz is version 2.3.1. gcc is version 5.4.0
** I think my patch level is actually 16.04.4 but I've not rebooted for 170 days :p
zgrep .fulfils any tangible purpose. You should be able to leave it off entirely (replaced withzcat). $\endgroup$zgrep .will only print non-blank lines (granted, that will still count lines that have nothing but whitespace since those aren't technically empty, but it's better than nothing). $\endgroup$