Answers.Body: <p>I get fairly quick results with my <a href="https://github.com/gringer/bioinfscripts/blob/master/fastx-fetch.pl" rel="noreferrer">fastx-length.pl</a> script, with the added bonus of being able to handle multi-line FASTQ files and displaying additional read-length QC statistics:</p>
  
  <pre><code>time zcat albacored_all.fastq.gz | /bioinf/scripts/fastx-length.pl &gt; /dev/null
  Total sequences: 301135
  Total length: 283.902419 Mb
  Longest sequence: 5.601 kb
  Shortest sequence: 6 b
  Mean Length: 942 b
  Median Length: 999 b
  N50: 111835 sequences; L50: 1.103 kb
  N90: 245243 sequences; L90: 608 b
  
  real    0m8,802s
  user    0m16,584s
  sys 0m0,260s
  </code></pre>
  
  <p>Versus the script you have provided:</p>
  
  <pre><code>zcat albacored_all.fastq.gz | awk 'NR%4==2{c++; l+=length($0)}
            END{
                  print "Number of reads: "c; 
                  print "Number of bases in reads: "l
                }'
  Number of reads: 301135
  Number of bases in reads: 283902419
  
  real    0m8,382s
  user    0m10,216s
  sys 0m0,332s
  </code></pre>
  
  <p>Cat to <code>/dev/null</code> for comparison:</p>
  
  <pre><code>time zcat albacored_all.fastq.gz &gt; /dev/null
  
  real    0m7,877s
  user    0m7,856s
  sys 0m0,020s
  </code></pre>
  
  <p>I suspect that something using <a href="https://github.com/lh3/bioawk" rel="noreferrer">bioawk</a> might be a bit faster (and similarly FASTQ-compliant).</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-30T13:52:05.287Z
Answers.Text: bioawk is not that efficient, at least in my tests: https://bioinformatics.stackexchange.com/a/961/292
Answers.UserDisplayName: bli
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-28T20:00:12.073Z
Answers.LastActivityDate: 2017-06-28T20:00:12.073Z
Answers.UserDisplayName: gringer
Answers.Body: <h1>pigz | awk | wc is the fastest method</h1>
  
  <p>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:</p>
  
  <p><a href="ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01815/sequence_read/ERR047740_1.filt.fastq.gz" rel="noreferrer">ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01815/sequence_read/ERR047740_1.filt.fastq.gz</a></p>
  
  <p>as my test file, the correct answers being:</p>
  
  <pre><code>Number of reads: 67051220
  Number of bases in reads: 6034609800
  </code></pre>
  
  <p>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 <code>time</code> on an otherwise unloaded system:</p>
  
  <h3>zgrep</h3>
  
  <pre><code>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
                }'
  </code></pre>
  
  <p>This is the slowest method with an average run-time of 125.35 seconds</p>
  
  <h3>gzip awk</h3>
  
  <p>Using <code>gzip</code> we gain about another 10 seconds:</p>
  
  <pre><code>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
                }'
  </code></pre>
  
  <p>Average run-time is 116.69 seconds</p>
  
  <h3>Konrad's gzip awk wc variant</h3>
  
  <pre><code>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
  </code></pre>
  
  <p>This runs slower on this test file than the <code>gzip awk</code> variant of the solution, average run-time is 122.28 seconds.</p>
  
  <h3>kseq_test using latest <code>kseq.h</code> from <a href="https://github.com/attractivechaos/klib" rel="noreferrer">klib</a></h3>
  
  <p>Code compiled with: <code>gcc -O2 -o kseq_test kseq_test.c -lz</code> where <code>kseq_test.c</code> is Simon's adaptation of Heng Li's FASTQ parser. </p>
  
  <p><code>kseq_test ERR047740_1.filt.fastq.gz</code></p>
  
  <p>Average run-time is 99.14 seconds, which is better than the <code>gzip</code> core utilities based solution so far, but we can do better! </p>
  
  <h3>piz awk</h3>
  
  <p>Using Mark Adler's <a href="https://zlib.net/pigz/" rel="noreferrer">pigz</a> as a drop-in replacement for <code>gzip</code>, note that <code>pigz</code> gives us a speed gain as on top of <code>gzip</code> as in addition to the main deflate thread it uses another 3 threads for reading, writing and checksum calculations, see the <a href="https://zlib.net/pigz/pigz.pdf" rel="noreferrer">man page</a> for details.    </p>
  
  <pre><code>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
                }'
  </code></pre>
  
  <p>Average run-time is now 93.86 seconds, <strong>this is ~5 seconds faster than the kseq based C code but we can further improve the benchmark</strong>.</p>
  
  <h3>pigz awk wc</h3>
  
  <p>Next we use <code>pigz</code> as a drop in replacment for Konrad's <code>wc</code> variant of the <code>awk</code> based solution.</p>
  
  <pre><code>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
  </code></pre>
  
  <p><strong>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 <code>zgrep</code> based solution.</strong></p>
  
  <p>Next as a baseline lets see just how much of this run-time is due to decompression of the input <code>fastq.gz</code> file.</p>
  
  <h3>gzip alone</h3>
  
  <p><code>gzip -dc ERR047740_1.filt.fastq.gz &gt; /dev/null</code></p>
  
  <p>Average run-time: 105.95 seconds, so the <code>gzip</code> based solutions (which also includes <code>zcat</code> and <code>zgrep</code> as these are provided by <code>gzip</code>) are never going to be faster than <code>kseq_test</code>.</p>
  
  <h3>pigz alone</h3>
  
  <p><code>pigz -dc ERR047740_1.filt.fastq.gz &gt; /dev/null</code></p>
  
  <p><strong>Average run-time: 77.66 seconds, so quite clearly the additional three threads for read, write and checksum calculation offer a useful advantage</strong>.  What's more this speed-up is greater when leveraging the <code>awk | wc</code> based solution, it's not clear why, but I expect this is due to the extra write thread.  </p>
  
  <p>Interestingly average CPU usage across all threads is quite revealing for the various answers, I've collated these stats using GNU time <code>/usr/bin/time --verbose</code></p>
  
  <p><code>zgrep</code> based solution 133% - must be more than one thread somehow</p>
  
  <p><code>gzip | awk</code> based solution 99% - all <code>gzip</code> based solutions run single-threaded at 99% CPU usage</p>
  
  <p><code>pigz | awk</code> 147% </p>
  
  <p><code>gzip | awk | wc</code> 99% as with <code>gzip</code></p>
  
  <p><code>pgiz | awk | wc</code> 155%</p>
  
  <p><code>kseq_test</code> 99%</p>
  
  <p><code>gzip &gt; dev/null</code> 99%</p>
  
  <p><code>pigz &gt; dev/null</code> 155%</p>
  
  <p>Whilst the main deflate thread in <code>pigz</code> 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.</p>
  
  <p>I'm using Ubuntu 16.04.2 LTS**, my <code>gzip</code>, <code>zcat</code>, <code>zgrep</code> versions are all gzip 1.6 and <code>pigz</code> is version 2.3.1.  <code>gcc</code> is version 5.4.0 </p>
  
  <p>** I think my patch level is actually 16.04.4 but I've not rebooted for 170 days :p</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-08-03T17:23:35.123Z
Answers.Text: wow excellent. thanks
Answers.UserDisplayName: Ahdee
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-03T08:08:52.987Z
Answers.Text: Nice and thorough! Thank you. I don't know why your tests showed pigz to be faster when it was slower in mine. I suspect it might be because I was testing using files stored on an NFS volume so my times could well have depended more on the network latency than anything else. I'll repeat your tests using a local drive and report back.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-03T10:28:16.943Z
Answers.Text: OK, confirmed. On your test file and on a local drive, `pigz` was indeed faster than `gzip` or `zcat | grep`.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-02T21:56:20.447Z
Answers.LastActivityDate: 2017-07-02T22:46:28.427Z
Answers.LastEditDate: 2017-07-02T22:46:28.427Z
Answers.UserDisplayName: Matt Bashton
Answers.Body: <p>The following is more than twice as fast; however, <code>wc</code> counts newline characters as well. We thus need to subtract the line count from the base count (using Bash):</p>
  <pre><code>fix_base_count() {
      local counts=(<span class="math-container">$(cat))
      echo "$</span>{counts[0]} $((counts[1] - counts[0]))&quot;
  }
  
  gunzip -c &quot;$file&quot; \
  | awk 'NR % 4 == 2' \
  | wc -cl \
  | fix_base_count
  </code></pre>
  <p>All the caveats from Simon’s comment apply: this assumes the “simple” FASTQ format, where each record consists of exactly four lines. I think this is true for all files produced by Illumina sequencers and downstream tools.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-29T09:19:24.487Z
Answers.Text: What did you test this on? I used a 2.6G exome fastq file and this approach (without the bash function, so even faster) took  1m6.754s while my version using awk to count everything took 1m14.228s. What did you measure to see a two-fold increase in speed?
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-29T10:01:03.443Z
Answers.Text: @terdon I’ve used the `time` command for benchmarking, and working on a 2.2GiB FASTQ file. The times on my machine were 9m5.628s and 3m5.825s real time, respectively. The user time in both cases was slightly elevated (11m vs 5m) to make the  difference a factor 2.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-29T10:07:46.587Z
Answers.Text: @terdon But given Simon’s answer and common sense, this is a fishy result indeed. I don’t know what my `awk` version is doing here.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2017-06-28T15:22:24.697Z
Answers.LastActivityDate: 2021-12-17T12:06:15.390Z
Answers.LastEditDate: 2021-12-17T12:06:15.390Z
Answers.UserDisplayName: Konrad Rudolph
Answers.Body: <p>I hava implemented <a href="https://github.com/jameslz/seqtk_utils/blob/master/seqtk_counts" rel="nofollow noreferrer">seqtk_counts</a> using kseq.h from  <a href="https://github.com/attractivechaos/klib" rel="nofollow noreferrer">klib</a> </p>
  
  <p>Just a few line of Codes:</p>
  
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;zlib.h&gt;
  #include "kseq.h"
  
  KSEQ_INIT(gzFile, gzread)
  
  int main(int argc, char *argv[]){
  
  gzFile fp;
  kseq_t *seq;
  
  int l = 0;
  
  int64_t total = 0;
  int64_t lines = 0;
  
  if (argc == 1) {
      fprintf(stderr, "Usage: %s &lt;fastq&gt; &lt;sample&gt;\n", argv[0]);
  
      return 1;
  }
  
  fp = strcmp(argv[1], "-")? gzopen(argv[1], "r") : gzdopen(fileno(stdin), "r");
  seq = kseq_init(fp);
  
  while ((l = kseq_read(seq)) &gt;= 0){
      total += seq-&gt;seq.l;
      lines += 1;
  }
  
  printf("%s\t%lld\t%lld\n", argv[2] ,(long long)lines, (long long)total);
  kseq_destroy(seq);
  gzclose(fp);
  return 0;
  </code></pre>
  
  <p>}</p>
  
  <p>Compile it:</p>
  
  <pre><code>gcc  -O2  seqtk_counts.c  -o  seqtk_counts  -Iklib  -lz
  </code></pre>
  
  <p>Usage:</p>
  
  <pre><code>seqtk_counts foo.fasq.gz foo 
  or
  cat foo.fasq.gz | seqtk_counts  - foo
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-29T09:45:03.563Z
Answers.Text: Thanks for the edit! :) But isn't this basically the same program as shown in  [sjcockell's answer](https://bioinformatics.stackexchange.com/a/937/298)?
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-29T09:05:43.780Z
Answers.LastActivityDate: 2017-06-29T09:33:47.703Z
Answers.LastEditDate: 2017-06-29T09:33:47.703Z
Answers.UserDisplayName: Lei ZHANG
Answers.Body: <p>if the data is in SRA, there is sra-stat utility that returns reads,bases and quality distribution. these are stored in the SRA file. use --quick to get the stored stats or --statistics to calculate additional values broken down per readgroup/barcode.
  sra-stat --quick SRR077487</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-26T17:47:55.083Z
Answers.LastActivityDate: 2017-07-26T17:47:55.083Z
Answers.UserDisplayName: user1203
Answers.Body: <p>Use assembly-stats!</p>
  <pre><code>$ assembly-stats barcode13_filtered.fastq 
  stats for barcode13_filtered.fastq
  sum = 2080834976, n = 656192, ave = 3171.08, largest = 15321
  N50 = 3598, n = 225863
  N60 = 3263, n = 286569
  N70 = 2920, n = 353904
  N80 = 2548, n = 429982
  N90 = 2081, n = 519671
  N100 = 500, n = 656192
  N_count = 0
  Gaps = 0
  
  real    0m7.418s
  user    0m5.976s
  sys     0m1.275s
  </code></pre>
  <p>Works on any fasta/fastq files, even those that aren't assemblies. And it's pretty fast!</p>
  <p>Download from here: <a href="https://anaconda.org/bioconda/assembly-stats" rel="nofollow noreferrer">https://anaconda.org/bioconda/assembly-stats</a></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-05-18T09:03:07.243Z
Answers.Text: Thanks, that looks interesting. How does the time compare to the existing answers here? Is this faster? And what are the numbers reported? Is that 656192 reads? And 2080834976 bases?
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-05-18T10:07:36.010Z
Answers.Text: Yes, yes, and yes! 656192 reads, 2080834976 bases, with the respective N50 etc values. And compared to the ones above, it is lightning fast.
  real 0m7.418s
  user 0m5.976s
Answers.UserDisplayName: winanonanona
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-05-18T10:11:13.847Z
Answers.Text: Could you show an example of that? You are running this on a single, decompressed file here, while the question is asking about compressed files and the other answers have shown that the rate limiting step is usually the decompression. So how fast is this approach when run on a compressed file and taking the decompression time into account? Also, in order for the reported times to be meaningful, please also include the time taken by one of the other methods (you can't compare times across different machines and files).
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-05-19T11:51:28.190Z
Answers.Text: This should too good to be true, is it because the `fastq` is uncompressed? @sjcockell benchmark is on a comparable ~2.6e9 bases file, but it's about 5x slower... But it's a compressed file that was used.
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-05-18T08:20:46.717Z
Answers.LastActivityDate: 2021-05-18T10:10:12.057Z
Answers.LastEditDate: 2021-05-18T10:10:12.057Z
Answers.UserDisplayName: winanonanona
Answers.Body: <h3>Using pyGATB</h3>
  
  <p>(I use the same file as in <a href="https://bioinformatics.stackexchange.com/a/400/292">https://bioinformatics.stackexchange.com/a/400/292</a>, same workstation as in <a href="https://bioinformatics.stackexchange.com/a/380/292">https://bioinformatics.stackexchange.com/a/380/292</a>)</p>
  
  <pre><code>$ time python3 -c "from gatb import Bank; seq_lens = [len(seq) for seq in Bank('SRR077487_2.filt.fastq.gz')]; print('Number of reads: %d' % len(seq_lens), 'Number of bases in reads: %d' % sum(seq_lens), sep='\n')"
  Number of reads: 23861612
  Number of bases in reads: 2386161200
  
  real    0m41.122s
  user    0m40.788s
  sys     0m0.312s
  </code></pre>
  
  <p>It is quite faster than bioawk:</p>
  
  <pre><code>$ time bioawk -c fastx '{nb_seq+=1; nb_char+=length($seq)} END {print "Number of reads: "nb_seq"\nNumber of bases in reads: "nb_char}' SRR077487_2.filt.fastq.gz
  Number of reads: 23861612
  Number of bases in reads: 2386161200
  
  real    1m3.182s
  user    1m2.916s
  sys     0m0.268s
  </code></pre>
  
  <p>But not so much than the OP example:</p>
  
  <pre><code>$ time zgrep . SRR077487_2.filt.fastq.gz | awk 'NR%4==2{c++; l+=length($0)} END{print "Number of reads: "c; print "Number of bases in reads: "l}'
  Number of reads: 23861612
  Number of bases in reads: 2386161200
  
  real    0m47.127s
  user    1m36.292s
  sys     0m6.796s
  </code></pre>
  
  <p>Or than <a href="https://bioinformatics.stackexchange.com/a/936/292">the <code>wc</code> based solution</a>:</p>
  
  <pre><code>$ fix_base_count() {
  &gt;     local counts=($(cat))
  &gt;     echo "${counts[0]} $((${counts[1]} - ${counts[0]}))"
  &gt; }
  $ time gunzip -c SRR077487_2.filt.fastq.gz | awk 'NR % 4 == 2' | wc -cl | fix_base_count
  23861612 2386161200
  
  real    0m44.915s
  user    1m12.000s
  sys     0m6.972s
  </code></pre>
  
  <p>I didn't compare with C-based solutions.</p>
  
  <p>The <a href="https://bioinformatics.stackexchange.com/a/937/292"><code>zcat</code> to <code>/dev/null</code> reference</a> is the following:</p>
  
  <pre><code>$ time zcat SRR077487_2.filt.fastq.gz &gt; /dev/null
  
  real    0m39.745s
  user    0m39.464s
  sys     0m0.252s
  </code></pre>
  
  <p>I'm still impressed by pyGATB speed</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-30T13:50:55.503Z
Answers.LastActivityDate: 2017-06-30T13:50:55.503Z
Answers.UserDisplayName: bli
Answers.Body: <p>It's difficult to get this to go massively quicker I think - as with <a href="https://bioinformatics.stackexchange.com/questions/361/what-is-the-fastest-way-to-calculate-the-number-of-unknown-nucleotides-in-fasta">this question</a> working with large gzipped FASTQ files is mostly IO-bound. We could instead focus on making sure we are getting the <em>right</em> answer.</p>
  
  <p>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 <a href="http://lh3lh3.users.sourceforge.net/parsefastq.shtml" rel="noreferrer">FASTQ Parser in C</a>. </p>
  
  <p>I downloaded the <a href="http://lh3lh3.users.sourceforge.net/download/kseq.tar.bz2" rel="noreferrer">example tarball</a> and modified the example code (excuse my C...):</p>
  
  <pre class="lang-c prettyprint-override"><code>#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;
  }
  </code></pre>
  
  <p>Then <code>make</code> and <code>kseq_test foo.fastq.gz</code>.</p>
  
  <p>For my example file (~35m reads of ~75bp) this took:</p>
  
  <pre><code>real    0m49.670s
  user    0m49.364s
  sys     0m0.304s
  </code></pre>
  
  <p>Compared with your example:</p>
  
  <pre><code>real    0m43.616s
  user    1m35.060s
  sys     0m5.240s
  </code></pre>
  
  <p>Konrad's solution (in my hands):</p>
  
  <pre><code>real    0m39.682s
  user    1m11.900s
  sys     0m5.112s
  </code></pre>
  
  <p>(By the way, just zcat-ing the data file to /dev/null):</p>
  
  <pre><code>real    0m38.736s
  user    0m38.356s
  sys     0m0.308s
  </code></pre>
  
  <p>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.</p>
  
  <p>And my horrible C can almost certainly be optimised.</p>
  
  <hr>
  
  <p>Same test, with <code>kseq.h</code> from Github, as suggested in the comments:</p>
  
  <p>My machine is under different load this morning, so I've retested. Wall clock times:</p>
  
  <p>OP: 0m44.813s</p>
  
  <p>Konrad: 0m40.061s</p>
  
  <p><code>zcat &gt; /dev/null</code>: 0m34.508s</p>
  
  <p>kseq.h (Github): 0m32.909s</p>
  
  <p>So most recent version of <code>kseq.h</code> is faster than simply zcat-ing the file (consistently in my tests...).</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-28T16:44:18.583Z
Answers.Text: Wow. So the C approach is actually marginally slower than the zgrep/awk I used and significantly slower than Konrad's trick of switching to `wc`.  That's surprising. I guess there must be some fancy parsing going on inside `kseq_read`.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-29T00:17:14.303Z
Answers.Text: The kseq.h from the tarball is outdated. The one [from github](https://github.com/attractivechaos/klib/blob/master/kseq.h) is faster. That said, decompression is the bottleneck as you have observed.
Answers.UserDisplayName: user172818
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-29T08:08:20.197Z
Answers.Text: @user172818 will that always be the case, irrespective of the size of the file? I would have thought that the decompression would be less of a bottleneck than the processing of each line when we're dealing with millions of lines. Is that wrong?
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-29T11:15:29.697Z
Answers.Text: I think with respect to speed of core utility base solutions, people often overlook just how optimised these programs are (which of course are themselves C binaries), they have been around for a long time and are not going away.
Answers.UserDisplayName: Matt Bashton
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-29T11:16:46.203Z
Answers.Text: Decompression will always be the bottleneck, it's the same situation we saw with https://bioinformatics.stackexchange.com/questions/361/what-is-the-fastest-way-to-calculate-the-number-of-unknown-nucleotides-in-fasta
Answers.UserDisplayName: Matt Bashton
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-28T16:15:06.267Z
Answers.LastActivityDate: 2017-06-29T09:28:15.997Z
Answers.LastEditDate: 2017-06-29T09:28:15.997Z
Answers.UserDisplayName: sjcockell
Body: <p>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 <code>zgrep</code> and <code>awk</code>:</p>
  <pre><code>zgrep . foo.fastq.gz |
       awk 'NR%4==2{c++; l+=length($0)}
            END{
                  print &quot;Number of reads: &quot;c; 
                  print &quot;Number of bases in reads: &quot;l
                }'
  </code></pre>
  <p>The files I need to process are regular ASCII text files (fastq) compressed using gzip (usually GNU gzip, possibly BSD gzip sometimes if clients are using macs).</p>
  <p><code>zgrep .</code> will print non-blank lines from the input file and the <code>awk 'NR%4==2</code> will process every 4th line starting with the second (the sequence).
  This works fine, but can take a <em>very</em> 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.</p>
  <p>I know that the FASTQ specification (such as it is) allows for line breaks in the sequence and qual strings, so simply taking the second of every group of 4 lines is not guaranteed to work (see <a href="https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2847217/#__sec7title" rel="nofollow noreferrer">here</a>). That's another reason why I'd rather use a more sophisticated, dedicated tool. However, I have yet to actually encounter a file with &gt;4 lines per record in the wild, so I am not too worried about that.</p>
  
Tags: <ngs><fastq><software-recommendation><benchmarking>
Title: Fast way to count number of reads and number of bases in a fastq file?
ThreadText: # Fast way to count number of reads and number of bases in a fastq file?
  
  (Asked by: username_0 on 2017-06-28)
  
  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.fastq.gz |
       awk 'NR%4==2{c++; l+=length($0)}
            END{
                  print "Number of reads: "c; 
                  print "Number of bases in reads: "l
                }'
  
  ```
  
  The files I need to process are regular ASCII text files (fastq) compressed using gzip (usually GNU gzip, possibly BSD gzip sometimes if clients are using macs).
  
  `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.
  
  I know that the FASTQ specification (such as it is) allows for line breaks in the sequence and qual strings, so simply taking the second of every group of 4 lines is not guaranteed to work (see [here](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2847217/#__sec7title)). That's another reason why I'd rather use a more sophisticated, dedicated tool. However, I have yet to actually encounter a file with >4 lines per record in the wild, so I am not too worried about that.
  
  **Tags:** ngs, fastq, software-recommendation, benchmarking
  
  
  ---
  
  ## 8 answers
  
  ### Answer by username_10 on 2017-06-28. Score: 14 | Accepted answer
  
  It's difficult to get this to go massively quicker I think - as with [this question](https://bioinformatics.stackexchange.com/questions/361/what-is-the-fastest-way-to-calculate-the-number-of-unknown-nucleotides-in-fasta) 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](http://lh3lh3.users.sourceforge.net/parsefastq.shtml).
  
  I downloaded the [example tarball](http://lh3lh3.users.sourceforge.net/download/kseq.tar.bz2) and modified the example code (excuse my C...):
  
  ```c
  #include <zlib.h>
  #include <stdio.h>
  #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 <in.seq>\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)) >= 0) {
          seqcount = seqcount + 1;
          seqlen = seqlen + (long)strlen(seq->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
  
  ```
  
  username_5'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
  
  username_5: 0m40.061s
  
  `zcat > /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...).
  
  > **Comments:**
  > **username_0**: Wow. So the C approach is actually marginally slower than the zgrep/awk I used and significantly slower than username_5's trick of switching to `wc`.  That's surprising. I guess there must be some fancy parsing going on inside `kseq_read`.  
  > **username_11**: The kseq.h from the tarball is outdated. The one [from github](https://github.com/attractivechaos/klib/blob/master/kseq.h) is faster. That said, decompression is the bottleneck as you have observed.  
  > **username_0**: @username_11 will that always be the case, irrespective of the size of the file? I would have thought that the decompression would be less of a bottleneck than the processing of each line when we're dealing with millions of lines. Is that wrong?  
  > **username_3**: I think with respect to speed of core utility base solutions, people often overlook just how optimised these programs are (which of course are themselves C binaries), they have been around for a long time and are not going away.  
  > **username_3**: Decompression will always be the bottleneck, it's the same situation we saw with https://bioinformatics.stackexchange.com/questions/361/what-is-the-fastest-way-to-calculate-the-number-of-unknown-nucleotides-in-fasta  
  
  ---
  
  ### Answer by username_3 on 2017-07-02. Score: 7
  
  # 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
  
  ### username_5'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](https://github.com/attractivechaos/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](https://zlib.net/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](https://zlib.net/pigz/pigz.pdf) 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 username_5'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
  
  > **Comments:**
  > **username_0**: Nice and thorough! Thank you. I don't know why your tests showed pigz to be faster when it was slower in mine. I suspect it might be because I was testing using files stored on an NFS volume so my times could well have depended more on the network latency than anything else. I'll repeat your tests using a local drive and report back.  
  > **username_0**: OK, confirmed. On your test file and on a local drive, `pigz` was indeed faster than `gzip` or `zcat | grep`.  
  > **username_4**: wow excellent. thanks  
  
  ---
  
  ### Answer by username_5 on 2017-06-28. Score: 7
  
  The following is more than twice as fast; however, `wc` counts newline characters as well. We thus need to subtract the line count from the base count (using Bash):
  
  ```
  fix_base_count() {
      local counts=($(cat))
      echo "${counts[0]} $((counts[1] - counts[0]))"
  }
  
  gunzip -c "$file" \
  | awk 'NR % 4 == 2' \
  | wc -cl \
  | fix_base_count
  
  ```
  
  All the caveats from Simon’s comment apply: this assumes the “simple” FASTQ format, where each record consists of exactly four lines. I think this is true for all files produced by Illumina sequencers and downstream tools.
  
  > **Comments:**
  > **username_0**: What did you test this on? I used a 2.6G exome fastq file and this approach (without the bash function, so even faster) took  1m6.754s while my version using awk to count everything took 1m14.228s. What did you measure to see a two-fold increase in speed?  
  > **username_5**: @username_0 I’ve used the `time` command for benchmarking, and working on a 2.2GiB FASTQ file. The times on my machine were 9m5.628s and 3m5.825s real time, respectively. The user time in both cases was slightly elevated (11m vs 5m) to make the  difference a factor 2.  
  > **username_5**: @username_0 But given Simon’s answer and common sense, this is a fishy result indeed. I don’t know what my `awk` version is doing here.  
  
  ---
  
  ### Answer by username_1 on 2017-06-28. Score: 5
  
  I get fairly quick results with my [fastx-length.pl](https://github.com/username_1/bioinfscripts/blob/master/fastx-fetch.pl) script, with the added bonus of being able to handle multi-line FASTQ files and displaying additional read-length QC statistics:
  
  ```
  time zcat albacored_all.fastq.gz | /bioinf/scripts/fastx-length.pl > /dev/null
  Total sequences: 301135
  Total length: 283.902419 Mb
  Longest sequence: 5.601 kb
  Shortest sequence: 6 b
  Mean Length: 942 b
  Median Length: 999 b
  N50: 111835 sequences; L50: 1.103 kb
  N90: 245243 sequences; L90: 608 b
  
  real    0m8,802s
  user    0m16,584s
  sys 0m0,260s
  
  ```
  
  Versus the script you have provided:
  
  ```
  zcat albacored_all.fastq.gz | awk 'NR%4==2{c++; l+=length($0)}
            END{
                  print "Number of reads: "c; 
                  print "Number of bases in reads: "l
                }'
  Number of reads: 301135
  Number of bases in reads: 283902419
  
  real    0m8,382s
  user    0m10,216s
  sys 0m0,332s
  
  ```
  
  Cat to `/dev/null` for comparison:
  
  ```
  time zcat albacored_all.fastq.gz > /dev/null
  
  real    0m7,877s
  user    0m7,856s
  sys 0m0,020s
  
  ```
  
  I suspect that something using [bioawk](https://github.com/lh3/bioawk) might be a bit faster (and similarly FASTQ-compliant).
  
  > **Comments:**
  > **username_2**: bioawk is not that efficient, at least in my tests: https://bioinformatics.stackexchange.com/a/961/292  
  
  ---
  
  ### Answer by username_6 on 2017-06-29. Score: 2
  
  I hava implemented [seqtk\_counts](https://github.com/jameslz/seqtk_utils/blob/master/seqtk_counts) using kseq.h from [klib](https://github.com/attractivechaos/klib)
  
  Just a few line of Codes:
  
  ```
  #include <stdio.h>
  #include <zlib.h>
  #include "kseq.h"
  
  KSEQ_INIT(gzFile, gzread)
  
  int main(int argc, char *argv[]){
  
  gzFile fp;
  kseq_t *seq;
  
  int l = 0;
  
  int64_t total = 0;
  int64_t lines = 0;
  
  if (argc == 1) {
      fprintf(stderr, "Usage: %s <fastq> <sample>\n", argv[0]);
  
      return 1;
  }
  
  fp = strcmp(argv[1], "-")? gzopen(argv[1], "r") : gzdopen(fileno(stdin), "r");
  seq = kseq_init(fp);
  
  while ((l = kseq_read(seq)) >= 0){
      total += seq->seq.l;
      lines += 1;
  }
  
  printf("%s\t%lld\t%lld\n", argv[2] ,(long long)lines, (long long)total);
  kseq_destroy(seq);
  gzclose(fp);
  return 0;
  
  ```
  
  }
  
  Compile it:
  
  ```
  gcc  -O2  seqtk_counts.c  -o  seqtk_counts  -Iklib  -lz
  
  ```
  
  Usage:
  
  ```
  seqtk_counts foo.fasq.gz foo 
  or
  cat foo.fasq.gz | seqtk_counts  - foo
  
  ```
  
  > **Comments:**
  > **username_0**: Thanks for the edit! :) But isn't this basically the same program as shown in  [username_10's answer](https://bioinformatics.stackexchange.com/a/937/298)?  
  
  ---
  
  ### Answer by username_7 on 2017-07-26. Score: 1
  
  if the data is in SRA, there is sra-stat utility that returns reads,bases and quality distribution. these are stored in the SRA file. use --quick to get the stored stats or --statistics to calculate additional values broken down per readgroup/barcode.
  sra-stat --quick SRR077487
  
  
  ---
  
  ### Answer by username_8 on 2021-05-18. Score: 1
  
  Use assembly-stats!
  
  ```
  $ assembly-stats barcode13_filtered.fastq 
  stats for barcode13_filtered.fastq
  sum = 2080834976, n = 656192, ave = 3171.08, largest = 15321
  N50 = 3598, n = 225863
  N60 = 3263, n = 286569
  N70 = 2920, n = 353904
  N80 = 2548, n = 429982
  N90 = 2081, n = 519671
  N100 = 500, n = 656192
  N_count = 0
  Gaps = 0
  
  real    0m7.418s
  user    0m5.976s
  sys     0m1.275s
  
  ```
  
  Works on any fasta/fastq files, even those that aren't assemblies. And it's pretty fast!
  
  Download from here: <https://anaconda.org/bioconda/assembly-stats>
  
  > **Comments:**
  > **username_0**: Thanks, that looks interesting. How does the time compare to the existing answers here? Is this faster? And what are the numbers reported? Is that 656192 reads? And 2080834976 bases?  
  > **username_8**: Yes, yes, and yes! 656192 reads, 2080834976 bases, with the respective N50 etc values. And compared to the ones above, it is lightning fast. real 0m7.418s user 0m5.976s  
  > **username_0**: Could you show an example of that? You are running this on a single, decompressed file here, while the question is asking about compressed files and the other answers have shown that the rate limiting step is usually the decompression. So how fast is this approach when run on a compressed file and taking the decompression time into account? Also, in order for the reported times to be meaningful, please also include the time taken by one of the other methods (you can't compare times across different machines and files).  
  > **username_9**: This should too good to be true, is it because the `fastq` is uncompressed? @username_10 benchmark is on a comparable ~2.6e9 bases file, but it's about 5x slower... But it's a compressed file that was used.  
  
  ---
  
  ### Answer by username_2 on 2017-06-30. Score: 1
  
  ### Using pyGATB
  
  (I use the same file as in <https://bioinformatics.stackexchange.com/a/400/292>, same workstation as in <https://bioinformatics.stackexchange.com/a/380/292>)
  
  ```
  $ time python3 -c "from gatb import Bank; seq_lens = [len(seq) for seq in Bank('SRR077487_2.filt.fastq.gz')]; print('Number of reads: %d' % len(seq_lens), 'Number of bases in reads: %d' % sum(seq_lens), sep='\n')"
  Number of reads: 23861612
  Number of bases in reads: 2386161200
  
  real    0m41.122s
  user    0m40.788s
  sys     0m0.312s
  
  ```
  
  It is quite faster than bioawk:
  
  ```
  $ time bioawk -c fastx '{nb_seq+=1; nb_char+=length($seq)} END {print "Number of reads: "nb_seq"\nNumber of bases in reads: "nb_char}' SRR077487_2.filt.fastq.gz
  Number of reads: 23861612
  Number of bases in reads: 2386161200
  
  real    1m3.182s
  user    1m2.916s
  sys     0m0.268s
  
  ```
  
  But not so much than the OP example:
  
  ```
  $ time zgrep . SRR077487_2.filt.fastq.gz | awk 'NR%4==2{c++; l+=length($0)} END{print "Number of reads: "c; print "Number of bases in reads: "l}'
  Number of reads: 23861612
  Number of bases in reads: 2386161200
  
  real    0m47.127s
  user    1m36.292s
  sys     0m6.796s
  
  ```
  
  Or than [the `wc` based solution](https://bioinformatics.stackexchange.com/a/936/292):
  
  ```
  $ fix_base_count() {
  >     local counts=($(cat))
  >     echo "${counts[0]} $((${counts[1]} - ${counts[0]}))"
  > }
  $ time gunzip -c SRR077487_2.filt.fastq.gz | awk 'NR % 4 == 2' | wc -cl | fix_base_count
  23861612 2386161200
  
  real    0m44.915s
  user    1m12.000s
  sys     0m6.972s
  
  ```
  
  I didn't compare with C-based solutions.
  
  The [`zcat` to `/dev/null` reference](https://bioinformatics.stackexchange.com/a/937/292) is the following:
  
  ```
  $ time zcat SRR077487_2.filt.fastq.gz > /dev/null
  
  real    0m39.745s
  user    0m39.464s
  sys     0m0.252s
  
  ```
  
  I'm still impressed by pyGATB speed
  
  
  ---
  
  