Answers.Body: <p>5 hours and no benchmarks posted? I am sorely disappointed.</p>
  
  <p>I'll restrict the comparison to just be fasta files, since fastq will end up being the same. So far, the contenders are:</p>
  
  <ol>
  <li>R with the <code>ShortRead</code> package (even if not the fastest, certainly a super convenient method).</li>
  <li>A pipeline of <code>grep -v "^&gt;" | tr -cd A | wc -c</code></li>
  <li>A pipeline of <code>grep -v "^&gt;" | grep -io A | wc -l</code></li>
  <li>A pipeline of <code>seqtk comp | awk</code></li>
  <li>Something custom in C/C++</li>
  </ol>
  
  <p>The R code can't actually count N records for some reason when I try it, so I counted <code>A</code> for everything.</p>
  
  <p>I'll exclude all python solutions, because there's no conceivable universe in which python is fast. For the C/C++ solution I just used what @user172818 posted (everything I tried writing took the same amount of time).</p>
  
  <p>Repeating 100 times and taking the average (yes, one should take the median):</p>
  
  <ol>
  <li>1.7 seconds</li>
  <li>0.65 seconds</li>
  <li>15 seconds</li>
  <li>1.2 seconds</li>
  <li>0.48 seconds</li>
  </ol>
  
  <p>Unsurprisingly, anything in straight C or C++ is going to win. <code>grep</code> with <code>tr</code> is surprisingly good, which surprises me since even though <code>grep</code> is very very optimized I still expected it to have more overhead. Piping to <code>grep -io</code> is a performance killer. The R solution is surprisingly good given the typical R overhead.</p>
  
  <p><strong>Update1:</strong> As suggested in the comments</p>
  
  <ol start="6">
  <li><code>sum(x.count('A') for x in open('seqs.fa','r') if x[0] != '&gt;')</code> in python</li>
  </ol>
  
  <p>This yields a time of</p>
  
  <ol start="6">
  <li>1.6 seconds (I'm at work now, so I've adjusted the time as best I can to account for the different computers)</li>
  </ol>
  
  <p><strong>Update2:</strong> More benchmarks from the comments</p>
  
  <ol start="7">
  <li><code>awk -F A '!/^&gt;"/ {cnt+=NF-1}END{print cnt}' foo.fa</code></li>
  <li><code>perl -ne 'if(!/^&gt;/){$count += tr/A//} END{print "$count\n"}' foo.fa</code></li>
  </ol>
  
  <p>These yield times of:</p>
  
  <ol start="7">
  <li>5.2 seconds</li>
  <li>1.18 seconds</li>
  </ol>
  
  <p>The awk version is the best hack I could come up with, there are likely better ways.</p>
  
  <p><strong>Update 3</strong>: Right, so regarding fastq files, I'm running this on a 3.3GB gzipped file with 10 repetitions (this takes a bit of time to run), so I'm not going to initially limit tests to commands that I can trivially modify to handle compressed files (after all, uncompressed fastq files are an abomination).</p>
  
  <ol>
  <li><code>seqtk fqchk</code></li>
  <li><code>sed -n '2~4p' &lt;(zcat foo.fastq.gz) | grep -io A | wc -l</code></li>
  <li><code>sed -n '1d;N;N;N;P;d' &lt;(zcat Undetermined_S0_R1_001.fastq.gz) | tr -cd A | wc -c</code></li>
  <li>The C example from user1728181 (the comparison to C++ is already provided there, so no need to include it).</li>
  <li><code>bioawk -c fastx '{n+=gsub(/A/, "", $seq)} END {print n}' foo.fastq.gz</code></li>
  </ol>
  
  <p>The average times are:</p>
  
  <ol>
  <li>1 minute 44 seconds</li>
  <li>2 minutes 6 seconds</li>
  <li>1 minute 52 seconds</li>
  <li>1 minute 15 seconds</li>
  <li>3 minutes 8 seconds</li>
  </ol>
  
  <p>N.B., these will often not handle fastq files with entries spanning more than 4 lines. However, such files are the dodo birds of bioinformatics, so that's not practically important. Also, I couldn't get <code>ShortRead</code> to work with the compressed fastq file. I suspect there's a way to fix that.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T07:53:08.817Z
Answers.Text: Nice benchmark, what is the size of your test fasta file? Also, why not add Python if you have time, yes it will not be fastest, but interesting to see for comparison, I would try: `sum(x.count('A') for x in open('seqs.fa','r') if x[0] != '>')`, probably `BioPython` will be slower for this
Answers.UserDisplayName: Chris_Rands
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T07:57:21.810Z
Answers.Text: @Chris_Rands I used dm6 (uncompressed), so it's not terribly large. I figured that python would be so self-evidently slower that I didn't want to bother, but I suppose it wouldn't hurt to add.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T08:08:24.853Z
Answers.Text: @Chris_Rands I've updated with a python benchmark. It does as well as R, which is better than I would have guessed.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T08:14:11.577Z
Answers.Text: Thanks! I'm not so surprised since python's `str.count` is implemented in C. I also feel there should be a fast `awk` only solution, but I'm not able to figure it out
Answers.UserDisplayName: Chris_Rands
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T08:20:46.863Z
Answers.Text: Perl: ```perl -ne 'if(!/^>/){$count +=  tr/A//} END{print "$count\n"}' dm6.fa```
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T08:37:00.513Z
Answers.Text: Okay, now that the question has been updated, you can modify these to work for the general FASTQ case (multi-line sequence + quality string, with '@' possibly appearing in the quality line(s)).
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T08:37:18.690Z
Answers.Text: My awk-fu is sufficient to come up with a working example, but not good enough to come up with one that's not painfully slow.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T08:38:13.573Z
Answers.Text: @gringer There goes me being productive today...
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T09:39:03.207Z
Answers.Text: Can you try my `bioawk -c fastx '{n+=gsub(/A/, "", $seq)} END {print n}'` ? I'm curious to see how it compares with the other solutions on your data.
Answers.UserDisplayName: bli
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T09:39:43.203Z
Answers.Text: @bli: Yup, I have 5 others running at the moment...
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T10:29:14.590Z
Answers.Text: The R code can count Ns once you write it correctly. ;-) (My bad)
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T10:45:04.937Z
Answers.Text: The awk and perl versions do not give the same results on my test fasta file: The awk version find one more "A" in C. elegans genome than the perl version. I guess it depends whether one of the sequence starts or ends with "A".
Answers.UserDisplayName: bli
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T11:34:17.793Z
Answers.Text: I've added a fastq comparison with some methods (including bioawk).
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T11:35:25.853Z
Answers.Text: @bli There should just be a `-1` in the awk code. Among the many reasons to use `bioawk` instead!
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T11:39:36.447Z
Answers.Text: I think the zcat is not necessary with bioawk. Does it improve the speed?
Answers.UserDisplayName: bli
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T11:41:12.797Z
Answers.Text: @bli It's not needed and that's a mistake on my part (I didn't actually use zcat in the benchmark). Thanks for catching that!
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T12:34:18.883Z
Answers.Text: Why are you counting the `A`s and not the `N`s? Or am I missing something?
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T12:35:28.177Z
Answers.Text: Because the original ShortRead variant couldn't handle N and the performance comparisons are the same regardless. "The R code can't actually count N records for some reason when I try it, so I counted A for everything."
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T13:04:55.853Z
Answers.Text: Yes, sorry, I missed that when I first read the answer.
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T21:07:06.653Z
Answers.Text: You say: "Repeating 100 times and taking the average (yes, one should take the median):" But the average is the mean not median - so which one did you use? I appreciate the median is less susceptible to skew but in this case surely run time is normally distributed (in which case mean is fine) unless you're running on a system with very high instantaneous CPU or IO load from other users which is creating contention.
Answers.UserDisplayName: Matt Bashton
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T21:20:41.933Z
Answers.Text: The mean, not the median, thus the caveat. Yes, run time should be normally distributed, but occasional weird IO spikes and such do happen.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-03T02:14:44.960Z
Answers.Text: Did you flush system cache between all these tests?
Answers.UserDisplayName: Torst
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-03T06:28:14.220Z
Answers.Text: Nah, in practice the files were in the cache for all of the runs.
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-23T14:32:41.410Z
Answers.Text: You should flush cache between tests. Your numbers will be misleading, otherwise, (well, wrong, really) regardless of whether you draw a mean or a median from test results.
Answers.UserDisplayName: Alex Reynolds
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T22:20:21.293Z
Answers.LastActivityDate: 2017-06-02T11:41:38.503Z
Answers.LastEditDate: 2017-06-02T11:41:38.503Z
Answers.UserDisplayName: Devon Ryan
Answers.Body: <p>I think that this should be pretty fast:</p>
  
  <p>FASTA:</p>
  
  <pre><code>grep -v "^&gt;" seqs.fa | tr -cd N | wc -c
  </code></pre>
  
  <p>FASTQ:</p>
  
  <pre><code>sed -n '1d;N;N;N;P;d' seqs.fq | tr -cd N | wc -c
  </code></pre>
  
  <p>See this <a href="https://stackoverflow.com/a/16679640/4641846">answer on SO</a> about how to count characters in BASH using different approaches.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T17:29:00.240Z
Answers.Text: Your FASTQ answer unfortunately has the same problem, fundamentally, as Iakov’s. That is, FASTQ records have a variable number of lines.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T17:33:24.777Z
Answers.Text: I know that the original standard admits FASTQ line splitting. However, have you ever seen such a file in practice?
Answers.UserDisplayName: Karel Břinda
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T17:34:39.940Z
Answers.Text: I haven’t, but I’ve never myself worked with Sanger sequencing/PacBio data. I think you might encounter them there. For short read data I guess the “record = 4 lines” assumption is safe.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T22:45:51.213Z
Answers.Text: Long read data do not use fastq. Currently PacBio uses bam (yes, really bam) and ONT fast5 (h5 archive). I think it is very reasonable to assume 4 line fastq.
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T17:22:20.160Z
Answers.LastActivityDate: 2017-06-01T17:22:20.160Z
Answers.UserDisplayName: Karel Břinda
Answers.Body: <p><strong>FASTQ</strong></p>
  
  <p>As it was pointed out, fastq can be complicated. But in a simple case when you have four lines per record, one possible solution in bash is:</p>
  
  <pre><code>sed -n '2~4p' seqs.fastq | grep -io N | wc -l
  </code></pre>
  
  <ul>
  <li><code>sed -n '2~4p'</code> will print every fourth line</li>
  <li><code>grep -o N</code> will output a line with <code>N</code> for every matching symbol</li>
  <li><code>wc -l</code> will count the lines</li>
  </ul>
  
  <p>I suspect this python approach will work faster:</p>
  
  <pre><code>cat seqs.fastq | python3 -c "import sys; print(sum(line.upper().count('N') for line in sys.stdin))"
  </code></pre>
  
  <p><strong>FASTA</strong></p>
  
  <p>Coreutils:</p>
  
  <pre><code>grep -v "^&gt;" seqs.fasta | grep -io N | wc -l
  </code></pre>
  
  <p>Python:</p>
  
  <pre><code>cat seqs.fasta | python3 -c "import sys; print(sum(line.upper().count('N') for line in sys.stdin if not line.startswith('&gt;')))"
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T17:03:05.813Z
Answers.Text: Now you’re counting lines containing Ns, not Ns. Also, FASTQ records can span more than four lines. :-(
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T17:07:41.883Z
Answers.Text: @KonradRudolph with -o grep outputs individual symbols, I will think how to correct fastq code
Answers.UserDisplayName: Iakov Davydov
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T17:14:58.503Z
Answers.Text: FASTQ is tricky because `@` can also start the quality string line; see http://chat.stackexchange.com/transcript/message/37809200#37809200 and following lines. For reference, here’s a Perl script that parses FASTQ and trims the mate information off of the identifier. It can probably serve as a starting point; I think this is very close to the minimal code: https://gist.github.com/klmr/81a2b86cd93c706d28611f40bdfe2702
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T18:11:57.120Z
Answers.Text: @KonradRudolph, I updated the answer to indicate that this works only in the simple case
Answers.UserDisplayName: Iakov Davydov
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T17:00:25.523Z
Answers.LastActivityDate: 2017-06-02T13:41:57.900Z
Answers.LastEditDate: 2017-06-02T13:41:57.900Z
Answers.UserDisplayName: Iakov Davydov
Answers.Body: <h2>Using bioawk</h2>
  <p>With <a href="https://github.com/lh3/bioawk" rel="nofollow noreferrer">bioawk</a> (counting &quot;A&quot; in the <em>C. elegans</em> genome, because there seem to be no &quot;N&quot; in this file), on my computer:</p>
  <pre><code><span class="math-container">$ time bioawk -c fastx '{n+=gsub(/A/, "", $</span>seq)} END {print n}' genome.fa 
  32371810
  
  real    0m1.645s
  user    0m1.548s
  sys     0m0.088s
  </code></pre>
  <p>bioawk is an extension of awk with convenient parsing options. For instance <code>-c fastx</code> makes the sequence available as <code>$seq</code> in fasta and fastq format (<strong>including gzipped versions</strong>), so the above command <strong>should also handle fastq format robustly</strong>.</p>
  <p>The <code>gsub</code> command is a normal awk function. It returns the number of substituted characters (got it from <a href="https://unix.stackexchange.com/a/169575/55127">https://unix.stackexchange.com/a/169575/55127</a>). I welcome comments about how to count inside awk in a less hackish way.</p>
  <p>On this fasta example it is almost 3 times slower than the <code>grep</code>, <code>tr</code>, <code>wc</code> pipeline:</p>
  <pre><code>$ time grep -v &quot;^&gt;&quot; genome.fa | tr -cd &quot;A&quot; | wc -c
  32371810
  
  real    0m0.609s
  user    0m0.744s
  sys     0m0.184s
  </code></pre>
  <p>(I repeated the timings, and the above values seem representative)</p>
  <h2>Using python</h2>
  <h3>readfq</h3>
  <p>I tried the python readfq retrieved from the repository mentioned in this answer: <a href="https://bioinformatics.stackexchange.com/a/365/292">https://bioinformatics.stackexchange.com/a/365/292</a></p>
  <pre><code>$ time python -c &quot;from sys import stdin; from readfq import readfq; print sum((seq.count('A') for _, seq, _ in readfq(stdin)))&quot; &lt; genome.fa
  32371810
  
  real    0m0.976s
  user    0m0.876s
  sys     0m0.100s
  </code></pre>
  <h3>&quot;pure python&quot;</h3>
  <p>Comparing with something adapted from this answer: <a href="https://bioinformatics.stackexchange.com/a/363/292">https://bioinformatics.stackexchange.com/a/363/292</a></p>
  <pre><code>$ time python -c &quot;from sys import stdin; print sum((line.count('A') for line in stdin if not line.startswith('&gt;')))&quot; &lt; genome.fa
  32371810
  
  real    0m1.143s
  user    0m1.100s
  sys     0m0.040s
  </code></pre>
  <p>(It is slower with python 3.6 than with python 2.7 on my computer)</p>
  <h3>pyGATB</h3>
  <p>I just learned (08/06/2017) that <a href="http://gatb.inria.fr/" rel="nofollow noreferrer">GATB</a> includes a fasta/fastq parser and has recently released a python API. I tried to use it yesterday to test another answer to the present question and <a href="https://github.com/GATB/pyGATB/issues/2" rel="nofollow noreferrer">found a bug</a>. This bug is now fixed, so here is a <a href="https://pypi.python.org/pypi/pyGATB/0.1.2" rel="nofollow noreferrer">pyGATB</a>-based answer:</p>
  <pre><code>$ time python3 -c &quot;from gatb import Bank; print(sum((seq.sequence.count(b\&quot;A\&quot;) for seq in Bank(\&quot;genome.fa\&quot;))))&quot;
  32371810
  
  real    0m0.663s
  user    0m0.568s
  sys     0m0.092s
  </code></pre>
  <p>(You can also do <code>sequence.decode(&quot;utf-8&quot;).count(&quot;A&quot;)</code> but this seems <a href="https://github.com/GATB/pyGATB/issues/2#issuecomment-307131176" rel="nofollow noreferrer">a little slower</a>.)</p>
  <p>Although I used python3.6 here (pyGATB seems python3-only), this is faster than the other two python approaches (for which the reported timings are obtained with python 2.7). This is even almost as fast as the <code>grep</code>, <code>tr</code>, <code>wc</code> pipeline.</p>
  <h3>Biopython</h3>
  <p>And, to have even more comparisons, here is a solution using <code>SeqIO.parse</code> from Biopython (with python2.7):</p>
  <pre><code>$ time python -c &quot;from Bio import SeqIO; print(sum((rec.seq.count(\&quot;A\&quot;) for rec in SeqIO.parse(\&quot;genome.fa\&quot;, \&quot;fasta\&quot;))))&quot;
  32371810
  
  real    0m1.632s
  user    0m1.532s
  sys     0m0.096s
  </code></pre>
  <p>This is a bit slower than the &quot;pure python&quot; solution, but perhaps more robust.</p>
  <p>There seems to be a slight improvement with @peterjc's suggestion to use the lower level <code>SimpleFastaParser</code>:</p>
  <pre><code>$ time python -c &quot;from Bio.SeqIO.FastaIO import SimpleFastaParser; print(sum(seq.count('A') for title, seq in SimpleFastaParser(open('genome.fa'))))&quot;
  32371810
  
  real    0m1.618s
  user    0m1.500s
  sys     0m0.116s
  </code></pre>
  <p>(I did a series of timings and tried to take one that seemed representative, but there's a lot of overlap with the higher-level parser's timings.)</p>
  <h3>scikit-bio</h3>
  <p>I read about <a href="http://scikit-bio.org" rel="nofollow noreferrer">scikit-bio</a> today (23/06/2017), which has a sequence reader.</p>
  <pre><code>$ time python3 -c  &quot;import skbio; print(sum(seq.count('A') for seq in skbio.io.read('genome.fa', format='fasta', verify=False)))&quot;
  32371810
  
  real    0m3.643s
  user    0m3.440s
  sys     0m1.228s
  </code></pre>
  <h3>pyfastx</h3>
  <p>I read about <a href="https://pyfastx.readthedocs.io/en/latest/index.html" rel="nofollow noreferrer">pyfastx</a> today (27/08/2020), which has a sequence reader.</p>
  <pre><code>$ time python3 -c &quot;from pyfastx import Fasta; print(sum(seq.count('A') for (_, seq) in Fasta('genome.fa', build_index=False)))&quot;
  32371810
  
  real    0m0.781s
  user    0m0.740s
  sys     0m0.040s
  </code></pre>
  <p>(The timing is done using python 3.6, on the same machine as previous timings from 3 years ago, so it should be comparable.)</p>
  <h3><code>fastx_read</code> from mappy</h3>
  <p>I realize (27/08/2020) that I hadn't added a benchmark for <code>fastx_read</code> from <a href="https://pypi.org/project/mappy/" rel="nofollow noreferrer">mappy</a> which I regularly use since quite some time already:</p>
  <pre><code>$ time python3 -c &quot;from mappy import fastx_read; print(sum(seq.count('A') for (_, seq, _) in fastx_read('genome.fa')))&quot;
  32371810
  
  real    0m0.731s
  user    0m0.648s
  sys     0m0.080s
  </code></pre>
  <h2>Benchmarks on a fastq.gz file</h2>
  <p>I tested some of the above solutions (or adaptations thereof), counting &quot;N&quot; in the same file that was used here: <a href="https://bioinformatics.stackexchange.com/a/400/292">https://bioinformatics.stackexchange.com/a/400/292</a></p>
  <h3>bioawk</h3>
  <pre><code><span class="math-container">$ time bioawk -c fastx '{n+=gsub(/N/, "", $</span>seq)} END {print n}' SRR077487_2.filt.fastq.gz
  306072
  
  real    1m9.686s
  user    1m9.376s
  sys     0m0.304s
  </code></pre>
  <h3>pigz + readfq python module</h3>
  <p>readfq doesn't complain and is very fast when I pass directly the compressed fastq, but returns something wrong, so don't forget to manually take care of the decompression.</p>
  <p>Here I tried with <code>pigz</code>:</p>
  <pre><code>$ time pigz -dc SRR077487_2.filt.fastq.gz | python -c &quot;from sys import stdin; from readfq import readfq; print sum((seq.count('N') for _, seq, _ in readfq(stdin)))&quot;
  306072
  
  real    0m52.347s
  user    1m40.716s
  sys     0m8.604s
  </code></pre>
  <p>The computer has 16 cores, but I suspect the limiting factor for <code>pigz</code> is reading from the disk: the processors are very far from running full speed.</p>
  <p>And with <code>gzip</code>:</p>
  <pre><code>$ time gzip -dc SRR077487_2.filt.fastq.gz | python -c &quot;from sys import stdin; from readfq import readfq; print sum((seq.count('N') for _, seq, _ in readfq(stdin)))&quot;
  306072
  
  real    0m49.448s
  user    1m31.984s
  sys     0m2.312s
  </code></pre>
  <p>Here gzip and python both used a full processor. Resource usage balance was slightly better. I work on a desktop computer. I suppose a computing server would take advantage of <code>pigz</code> better.</p>
  <h3>pyGATB</h3>
  <pre><code>$ time python3 -c &quot;from gatb import Bank; print(sum((seq.sequence.count(b\&quot;N\&quot;) for seq in Bank(\&quot;SRR077487_2.filt.fastq.gz\&quot;))))&quot;
  306072
  
  real    0m46.784s
  user    0m46.404s
  sys     0m0.296s
  </code></pre>
  <h3>biopython</h3>
  <pre><code>$ time python -c &quot;from gzip import open as gzopen; from Bio.SeqIO.QualityIO import FastqGeneralIterator; print(sum(seq.count('N') for (_, seq, _) in FastqGeneralIterator(gzopen('SRR077487_2.filt.fastq.gz'))))&quot;
  306072
  
  real    3m18.103s
  user    3m17.676s
  sys     0m0.428s
  </code></pre>
  <h3>pyfastx (added 27/08/2020)</h3>
  <pre><code>$ time python3 -c &quot;from pyfastx import Fastq; print(sum(seq.count('N') for (_, seq, _) in Fastq('SRR077487_2.filt.fastq.gz', build_index=False)))&quot;
  306072
  
  real    0m51.140s
  user    0m50.784s
  sys     0m0.352s
  </code></pre>
  <h3><code>fastx_read</code> from mappy (added 27/08/2020)</h3>
  <pre><code>$ time python3 -c &quot;from mappy import fastx_read; print(sum(seq.count('N') for (_, seq, _) in fastx_read('SRR077487_2.filt.fastq.gz')))&quot;
  306072
  
  real    0m52.336s
  user    0m52.024s
  sys     0m0.300s
  </code></pre>
  <h2>Conclusion</h2>
  <p>pyGATB and bioawk both handle transparently compression (gzipped or not) and format differences (fasta or fastq). pyGATB is quite a new tool, but seems more efficient compared to the other python modules I tested.</p>
  <p>(Update 27/08/2020) pyfastx and mappy seem quite convenient too, and are only slightly slower than pyGATB.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-09T11:38:32.843Z
Answers.Text: Could you add the Biopython low-level string FASTA parser numbers? ``time python -c "from Bio.SeqIO.FastaIO import SimpleFastaParser; print(sum(seq.count('A') for title, seq in SimpleFastaParser(open('genome.fa'))))"``
Answers.UserDisplayName: Peter Cock
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-09T11:39:05.597Z
Answers.Text: Why do all your examples count the letter A, rather than N (and perhaps n) as per the question?
Answers.UserDisplayName: Peter Cock
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-09T12:43:10.990Z
Answers.Text: I count letter "A" because the first large file I thought about for my testing happened to have no "N".
Answers.UserDisplayName: bli
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-10T21:59:36.417Z
Answers.Text: I would think that all the Python examples will score badly on this small example file due to the overhead of loading Python and importing libraries. But if the user only wants to count one genome, that is a fair test all the same.
Answers.UserDisplayName: Peter Cock
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2017-06-02T09:28:46.290Z
Answers.LastActivityDate: 2020-08-27T15:33:15.750Z
Answers.LastEditDate: 2020-08-27T15:33:15.750Z
Answers.UserDisplayName: bli
Answers.Body: <p>If it is raw speed you're after, then writing an own little C/C++ program is probably what you need to do.</p>
  
  <p>Fortunately, the worst part (a fast and reliable parser) has already been tackled: the <a href="https://github.com/lh3/readfq" rel="noreferrer">readfq</a> from Heng Li is probably the fastest FASTA/FASTQ parser around.</p>
  
  <p>And it's easy to use, the example on GitHub can easily be expanded to do what you need. Just add in a parameter parser (for the filename you want to analyse), program a simple 'N'-counter loop and have the results printed to stdout. Done.</p>
  
  <p>Else, if it is simplicity you are after, see Konrad's answer for an R based solution. Or a very simple script using Biopython.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T17:19:11.827Z
Answers.Text: God that’s absolutely atrocious code. ☹️ I think that for raw speed, [SeqAn](http://docs.seqan.de/seqan/master/?p=SeqFileIn) (C++) is probably on par, and a lot higher quality.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T17:26:00.533Z
Answers.Text: It's not the prettiest around. Then again, it does what it needs to do and doesn't come with a whole library dependency. Do you have a benchmark with SeqAn & readfq?
Answers.UserDisplayName: BaCh
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-08T11:57:15.183Z
Answers.Text: I'm curious to see how readfq ans SeqAn compare with the fasta/fastq parser from GATB: http://gatb-core.gforge.inria.fr/doc/api/classgatb_1_1core_1_1bank_1_1impl_1_1BankFasta.html#details
  I'm not really able to code in C/C++ though.
Answers.UserDisplayName: bli
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T17:15:22.907Z
Answers.LastActivityDate: 2017-06-01T17:15:22.907Z
Answers.UserDisplayName: BaCh
Answers.Body: <p>Honestly, the easiest way (especially for FASTQ) is probably to use a dedicated parsing library, such as R/Bioconductor:</p>
  
  
  
  <pre><code>suppressMessages(library(ShortRead))
  
  seq = readFasta(commandArgs(TRUE)[1]) # or readFastq
  cat(colSums(alphabetFrequency(sread(seq))[, 'N', drop = FALSE]), '\n')
  </code></pre>
  
  <p>This may not be the fastest, but it’s <em>pretty fast</em>, since the relevant methods are highly optimised. The biggest overhead is actually probably from loading <code>ShortRead</code>, which is an unjustifiable slog.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T17:10:02.993Z
Answers.LastActivityDate: 2017-06-02T14:53:35.553Z
Answers.LastEditDate: 2017-06-02T14:53:35.553Z
Answers.UserDisplayName: Konrad Rudolph
Answers.Body: <p>For FASTA files, I've implemented a relatively efficient method in <a href="https://github.com/mdshw5/pyfaidx" rel="nofollow noreferrer">pyfaidx</a> v0.4.9.1. This post made me realize that my previous code was quite slow and easy to replace:</p>
  
  <pre><code>$ pip install pyfaidx
  $ time faidx -i nucleotide ~/Downloads/hg38.fa
  name    start   end A   T   C   G   N   others
  chr1    1   248956422   67070277    67244164    48055043    48111528    18475410    
  chr10   1   133797422   38875926    39027555    27639505    27719976    534460  
  chr11   1   135086622   39286730    39361954    27903257    27981801    552880  
  chr11_KI270721v1_random 1   100316  18375   19887   31042   31012   0   
  ...
  chrX    1   156040895   46754807    46916701    30523780    30697741    1147866 
  chrY    1   57227415    7886192 7956168 5285789 5286894 30812372    
  chrY_KI270740v1_random  1   37240   8274    12978   7232    8756    0   
  
  real    2m28.251s
  user    2m15.548s
  sys 0m9.951s
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-15T15:47:20.327Z
Answers.LastActivityDate: 2017-06-15T15:47:20.327Z
Answers.UserDisplayName: Matt Shirley
Answers.Body: <h3>Decompression of gzipped FASTQ is the main issue</h3>
  
  <p>If we take real world gzipped FASTQ (which as the OP suggested would be beneficial) rather than trivial FASTA as the starting point then the real issue is actually decompressing the file not counting the Ns and in this case the C program <code>count-N</code> is no longer the fastest solution.</p>
  
  <p>Additionally it would be good to use a specific file for benchmarking which actually has Ns, because you'll get some quite interesting execution time differences with some methods counting the more frequently occurring As rather than Ns.</p>
  
  <p>One such file is:</p>
  
  <p><a href="http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/sequence_read/SRR077487_2.filt.fastq.gz" rel="noreferrer">http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/sequence_read/SRR077487_2.filt.fastq.gz</a></p>
  
  <p>It's also worth checking the various solutions return the correct answer there should be 306072 Ns in the above file.</p>
  
  <p>Next note that decompression of this file redirected to <code>/dev/null</code> is slower with <code>zcat</code> and <code>gzip</code> (which are both gzip 1.6 on my system) than say a parallel implementation of gzip like <a href="http://zlib.net/pigz/" rel="noreferrer">Mark Adler's</a> <code>pigz</code>, which appears to use 4 threads for decompression.  All timings represent an average of 10 runs reporting real (wall clock time).</p>
  
  <pre><code>time pigz -dc SRR077487_2.filt.fastq.gz &gt; /dev/null
  
  real    0m29.0132s
  
  time gzip -dc SRR077487_2.filt.fastq.gz &gt; /dev/null
  
  real    0m40.6996s
  </code></pre>
  
  <p>There is an ~11.7 second difference between the two.  Next if I then try to benchmark a one-liner which performs on FASTQ and gives the correct answer (Note I've yet to encounter FASTQ which is not 4 line, and seriously who generates these files!)</p>
  
  <pre><code>time pigz -dc SRR077487_2.filt.fastq.gz | awk 'NR%4==2{print $1}' | tr -cd N | wc -c
  306072
  
  real    0m34.793s
  </code></pre>
  
  <p>As you can see the counting adds ~5.8 second to the total run time versus <code>pigz</code> based decompression.  Additionally this time delta is higher when using <code>gzip</code> ~6.7 seconds above gzip decompression alone.   </p>
  
  <pre><code>time gzip -dc SRR077487_2.filt.fastq.gz | awk 'NR%4==2{print $1}' | tr -cd N | wc -c                                                                
  306072
  
  real    0m44.399s
  </code></pre>
  
  <p>The <code>pigz awk tr wc</code> based solution is however ~4.5 seconds faster than the <code>count-N</code> based C code solution:</p>
  
  <pre><code>time count-N SRR077487_2.filt.fastq.gz 
  2385855128      306072  0
  
  real    0m39.266s
  </code></pre>
  
  <p>This difference appears to be robust to re-running as many times as I like.  I expect if you could use pthread in the C based solution or alter it to take the standard out from pigz it would also show an increase in performance.</p>
  
  <p>Benchmarking another alternative <code>pigz grep</code> variant appears to take more or less the same time as the <code>tr</code> based variant:</p>
  
  <pre><code>time pigz -dc SRR077487_2.filt.fastq.gz | awk 'NR%4==2{print $1}' | grep -o N | wc -l
  306072
  
  real    0m34.869s
  </code></pre>
  
  <p>Note that the <code>seqtk</code> based solution discussed above is noticeably slower,</p>
  
  <pre><code>time seqtk comp  SRR077487_2.filt.fastq.gz | awk '{x+=$9}END{print x}'
  306072
  
  real    1m42.062s 
  </code></pre>
  
  <p>However it's worth noting <code>seqtk comp</code> is doing a bit more than the other solutions.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T20:28:26.517Z
Answers.LastActivityDate: 2017-06-02T20:57:56.203Z
Answers.LastEditDate: 2017-06-02T20:57:56.203Z
Answers.UserDisplayName: Matt Bashton
Answers.Body: <p>Here are perl 1-liners for fasta and fastq files that are fast enough if you just want totals:</p>
  
  <p><code>$ perl -ne 'BEGIN { my $n=0 }; /^&gt;/ or $n += s/N//gi; END { print "N=$n\n" }' contigs.fasta</code></p>
  
  <p>For each line of a FASTA file, if not header (/^>/) then it's a sequence line. 
   Substitute Ns, which returns the number of substitutions made, which is added to the total.</p>
  
  <p><code>$ zcat reads.fastq | perl -ne 'BEGIN { my $n=0 }; /^[ATCGN]+$/i and $n += s/N//gi; END { print "N=$n\n" }'</code></p>
  
  <p>For each sequence line of a FASTQ file, matched by /^[ATCGN]+$/i (i.e. all NTs), then count substitutions made in the same way.</p>
  
  <p>You may use pigz instead of zcat (or gunzip -c) if you have pigz installed.</p>
  
  <p>An uncompressed fastq file of 14M 250bp Illumina reads completed in 40sec.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-01T09:02:34.747Z
Answers.Text: If you could explain what each perl line do, it would improve the answer, specially commenting the differences between both approaches
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-01T05:55:22.123Z
Answers.LastActivityDate: 2018-03-01T23:15:30.270Z
Answers.LastEditDate: 2018-03-01T23:15:30.270Z
Answers.UserDisplayName: Edward Kirton
Answers.Body: <p>For FASTQ:</p>
  
  <pre><code>seqtk fqchk in.fq | head -2
  </code></pre>
  
  <p>It gives you percentage of "N" bases, not the exact count, though.</p>
  
  <p>For FASTA:</p>
  
  <pre class="lang-sh prettyprint-override"><code>seqtk comp in.fa | awk '{x+=$9}END{print x}'
  </code></pre>
  
  <p>This command line also works with FASTQ, but it will be slower as awk is slow.</p>
  
  <p>EDIT: ok, based on @BaCH's reminder, here we go (you need <a href="https://github.com/lh3/readfq/blob/master/kseq.h" rel="noreferrer">kseq.h</a> to compile):</p>
  
  <pre class="lang-c prettyprint-override"><code>// to compile: gcc -O2 -o count-N this-prog.c -lz
  #include &lt;zlib.h&gt;
  #include &lt;stdio.h&gt;
  #include &lt;stdint.h&gt;
  #include "kseq.h"
  KSEQ_INIT(gzFile, gzread)
  
  unsigned char dna5tbl[256] = {
      0, 1, 2, 3,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 5, 4, 4,
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4
  };
  
  int main(int argc, char *argv[]) {
      long i, n_n = 0, n_acgt = 0, n_gap = 0;
      gzFile fp;
      kseq_t *seq;
      if (argc == 1) {
          fprintf(stderr, "Usage: count-N &lt;in.fa&gt;\n");
          return 1;
      }
      if ((fp = gzopen(argv[1], "r")) == 0) {
          fprintf(stderr, "ERROR: fail to open the input file\n");
          return 1;
      }
      seq = kseq_init(fp);
      while (kseq_read(seq) &gt;= 0) {
          for (i = 0; i &lt; seq-&gt;seq.l; ++i) {
              int c = dna5tbl[(unsigned char)seq-&gt;seq.s[i]];
              if (c &lt; 4) ++n_acgt;
              else if (c == 4) ++n_n;
              else ++n_gap;
          }
      }
      kseq_destroy(seq);
      gzclose(fp);
      printf("%ld\t%ld\t%ld\n", n_acgt, n_n, n_gap);
      return 0;
  }
  </code></pre>
  
  <p>It works for both FASTA/Q and gzip'ed FASTA/Q. The following uses SeqAn:</p>
  
  <pre class="lang-cpp prettyprint-override"><code>#include &lt;seqan/seq_io.h&gt;
  
  using namespace seqan;
  
  int main(int argc, char *argv[]) {
      if (argc == 1) {
          std::cerr &lt;&lt; "Usage: count-N &lt;in.fastq&gt;" &lt;&lt; std::endl;
          return 1;
      }
      std::ios::sync_with_stdio(false);
      CharString id;
      Dna5String seq;
      SeqFileIn seqFileIn(argv[1]);
      long i, n_n = 0, n_acgt = 0;
      while (!atEnd(seqFileIn)) {
          readRecord(id, seq, seqFileIn);
          for (i = beginPosition(seq); i &lt; endPosition(seq); ++i)
              if (seq[i] &lt; 4) ++n_acgt;
              else ++n_n;
      }
      std::cout &lt;&lt; n_acgt &lt;&lt; '\t' &lt;&lt; n_n &lt;&lt; std::endl;
      return 0;
  }
  </code></pre>
  
  <p>On a FASTQ with 4-million 150bp reads:</p>
  
  <ul>
  <li>The C version: ~0.74 sec</li>
  <li>The C++ version: ~2.15 sec</li>
  <li>An older C version without a lookup table (see <a href="https://bioinformatics.stackexchange.com/posts/369/revisions">the previous edit</a>): ~2.65 sec</li>
  </ul>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-04-21T09:35:25.373Z
Answers.Text: Out of curiosity I just checked the performance with a straightforward implementation in seqan3 and to my great disappointment it’s still a factor 3 slower than kseq.h. It’s of course entirely possible that I did something wrong but I don’t think so.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T08:34:23.513Z
Answers.Text: based on Devon Ryan's benchmarking you are the fastest so far...
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T08:54:52.737Z
Answers.Text: Can you add `#include <inttypes.h>` to you first C example?
Answers.UserDisplayName: Devon Ryan
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T10:33:44.673Z
Answers.Text: I can’t not upvote SeqAn, but I’m seriously concerned that your C++ code is so much slower than your C code. They should be virtually identical. That said, you didn’t specify `std::ios::sync_with_stdio(false)`, which might account for some of the difference; and the loop could be simplified by using iterators, though I’d be surprised if this impacted the speed on any setting above `-O1`.
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T14:35:38.837Z
Answers.Text: @KonradRudolph `sync_with_stdio` helps little when reading from a normal file. I don't know the internal of seqan, so I can't explain why it is slower. That said, I wouldn't worry too much about seqan performance. Gzip decompression is much slower than parsing fastq. @DevonRyan added.
Answers.UserDisplayName: user172818
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T15:01:26.313Z
Answers.Text: The difference is in readfq. That beast is *fast*. I know, I tried to replace Heng's code by something less ugly in C++ and couldn't get anywhere near it in simple read performance. Best I could come up with in clean C++ was a speed factor of 2x slower.
Answers.UserDisplayName: BaCh
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-02T20:29:52.297Z
Answers.Text: you need to add `-lz` to your compile flags (at least for the current Ubuntu LTS) `gcc -O2 -o count-N this-prog.c -lz`
Answers.UserDisplayName: Matt Bashton
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-23T14:35:55.220Z
Answers.Text: Are your numbers based on flushing the cache between samples? Otherwise, those numbers will almost certainly be incorrect.
Answers.UserDisplayName: Alex Reynolds
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-01T21:46:22.750Z
Answers.LastActivityDate: 2017-06-09T01:02:55.120Z
Answers.LastEditDate: 2017-06-09T01:02:55.120Z
Answers.UserDisplayName: user172818
Body: <p>I used to work with publicly available genomic references, where basic statistics are usually available and if they are not, you have to compute them only once so there is no reason to worry about performance.</p>
  
  <p>Recently I started sequencing project of couple of different species with mid-sized genomes (~Gbp) and during testing of different assembly pipelines I had compute number of unknown nucleotides many times in both raw reads (in fastq) and assembly scaffolds (in fasta), therefore I thought that I would like to optimize the computation.</p>
  
  <ul>
  <li>For me it is reasonable to expect 4-line formatted fastq files, but general solution is still prefered</li>
  <li>It would be nice if solution would work on gzipped files as well</li>
  </ul>
  
  <p>Q : What is the fastest way (performance-wise) to compute the number of unknown nucleotides (Ns) in fasta and fastq files?</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-06-02T11:56:01.010Z
Comments.Text: Presumably, you might want the total number of nucleotides as well as the number of Ns (to give proportion of Ns). This is a fairly fast implementation that provides this plus a few other metrics for fasta: https://github.com/ENCODE-DCC/kentUtils/blob/master/src/utils/faSize/faSize.c
Comments.UserDisplayName: Chris_Rands
Tags: <fasta><fastq><benchmarking>
Title: What is the fastest way to calculate the number of unknown nucleotides in FASTA / FASTQ files?
ThreadText: # What is the fastest way to calculate the number of unknown nucleotides in FASTA / FASTQ files?
  
  (Asked by: username_0 on 2017-06-01)
  
  I used to work with publicly available genomic references, where basic statistics are usually available and if they are not, you have to compute them only once so there is no reason to worry about performance.
  
  Recently I started sequencing project of couple of different species with mid-sized genomes (~Gbp) and during testing of different assembly pipelines I had compute number of unknown nucleotides many times in both raw reads (in fastq) and assembly scaffolds (in fasta), therefore I thought that I would like to optimize the computation.
  
  * For me it is reasonable to expect 4-line formatted fastq files, but general solution is still prefered
  * It would be nice if solution would work on gzipped files as well
  
  Q : What is the fastest way (performance-wise) to compute the number of unknown nucleotides (Ns) in fasta and fastq files?
  
  **Tags:** fasta, fastq, benchmarking
  
  > **Comments:**
  > **username_1**: Presumably, you might want the total number of nucleotides as well as the number of Ns (to give proportion of Ns). This is a fairly fast implementation that provides this plus a few other metrics for fasta: https://github.com/ENCODE-DCC/kentUtils/blob/master/src/utils/faSize/faSize.c  
  
  ---
  
  ## 10 answers
  
  ### Answer by username_17 on 2017-06-01. Score: 26 | Accepted answer
  
  For FASTQ:
  
  ```
  seqtk fqchk in.fq | head -2
  
  ```
  
  It gives you percentage of "N" bases, not the exact count, though.
  
  For FASTA:
  
  ```sh
  seqtk comp in.fa | awk '{x+=$9}END{print x}'
  
  ```
  
  This command line also works with FASTQ, but it will be slower as awk is slow.
  
  EDIT: ok, based on @BaCH's reminder, here we go (you need [kseq.h](https://github.com/lh3/readfq/blob/master/kseq.h) to compile):
  
  ```c
  // to compile: gcc -O2 -o count-N this-prog.c -lz
  #include <zlib.h>
  #include <stdio.h>
  #include <stdint.h>
  #include "kseq.h"
  KSEQ_INIT(gzFile, gzread)
  
  unsigned char dna5tbl[256] = {
      0, 1, 2, 3,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 5, 4, 4,
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4, 
      4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4
  };
  
  int main(int argc, char *argv[]) {
      long i, n_n = 0, n_acgt = 0, n_gap = 0;
      gzFile fp;
      kseq_t *seq;
      if (argc == 1) {
          fprintf(stderr, "Usage: count-N <in.fa>\n");
          return 1;
      }
      if ((fp = gzopen(argv[1], "r")) == 0) {
          fprintf(stderr, "ERROR: fail to open the input file\n");
          return 1;
      }
      seq = kseq_init(fp);
      while (kseq_read(seq) >= 0) {
          for (i = 0; i < seq->seq.l; ++i) {
              int c = dna5tbl[(unsigned char)seq->seq.s[i]];
              if (c < 4) ++n_acgt;
              else if (c == 4) ++n_n;
              else ++n_gap;
          }
      }
      kseq_destroy(seq);
      gzclose(fp);
      printf("%ld\t%ld\t%ld\n", n_acgt, n_n, n_gap);
      return 0;
  }
  
  ```
  
  It works for both FASTA/Q and gzip'ed FASTA/Q. The following uses SeqAn:
  
  ```cpp
  #include <seqan/seq_io.h>
  
  using namespace seqan;
  
  int main(int argc, char *argv[]) {
      if (argc == 1) {
          std::cerr << "Usage: count-N <in.fastq>" << std::endl;
          return 1;
      }
      std::ios::sync_with_stdio(false);
      CharString id;
      Dna5String seq;
      SeqFileIn seqFileIn(argv[1]);
      long i, n_n = 0, n_acgt = 0;
      while (!atEnd(seqFileIn)) {
          readRecord(id, seq, seqFileIn);
          for (i = beginPosition(seq); i < endPosition(seq); ++i)
              if (seq[i] < 4) ++n_acgt;
              else ++n_n;
      }
      std::cout << n_acgt << '\t' << n_n << std::endl;
      return 0;
  }
  
  ```
  
  On a FASTQ with 4-million 150bp reads:
  
  * The C version: ~0.74 sec
  * The C++ version: ~2.15 sec
  * An older C version without a lookup table (see [the previous edit](https://bioinformatics.stackexchange.com/posts/369/revisions)): ~2.65 sec
  
  > **Comments:**
  > **username_0**: based on username_2 username_2's benchmarking you are the fastest so far...  
  > **username_2**: Can you add `#include <inttypes.h>` to you first C example?  
  > **username_5**: I can’t not upvote SeqAn, but I’m seriously concerned that your C++ code is so much slower than your C code. They should be virtually identical. That said, you didn’t specify `std::ios::sync_with_stdio(false)`, which might account for some of the difference; and the loop could be simplified by using iterators, though I’d be surprised if this impacted the speed on any setting above `-O1`.  
  > **username_17**: @username_5 `sync_with_stdio` helps little when reading from a normal file. I don't know the internal of seqan, so I can't explain why it is slower. That said, I wouldn't worry too much about seqan performance. Gzip decompression is much slower than parsing fastq. @username_2 added.  
  > **username_13**: The difference is in readfq. That beast is *fast*. I know, I tried to replace Heng's code by something less ugly in C++ and couldn't get anywhere near it in simple read performance. Best I could come up with in clean C++ was a speed factor of 2x slower.  
  > **username_7**: you need to add `-lz` to your compile flags (at least for the current Ubuntu LTS) `gcc -O2 -o count-N this-prog.c -lz`  
  > **username_9**: Are your numbers based on flushing the cache between samples? Otherwise, those numbers will almost certainly be incorrect.  
  > **username_5**: Out of curiosity I just checked the performance with a straightforward implementation in seqan3 and to my great disappointment it’s still a factor 3 slower than kseq.h. It’s of course entirely possible that I did something wrong but I don’t think so.  
  
  ---
  
  ### Answer by username_2 on 2017-06-01. Score: 18
  
  5 hours and no benchmarks posted? I am sorely disappointed.
  
  I'll restrict the comparison to just be fasta files, since fastq will end up being the same. So far, the contenders are:
  
  1. R with the `ShortRead` package (even if not the fastest, certainly a super convenient method).
  2. A pipeline of `grep -v "^>" | tr -cd A | wc -c`
  3. A pipeline of `grep -v "^>" | grep -io A | wc -l`
  4. A pipeline of `seqtk comp | awk`
  5. Something custom in C/C++
  
  The R code can't actually count N records for some reason when I try it, so I counted `A` for everything.
  
  I'll exclude all python solutions, because there's no conceivable universe in which python is fast. For the C/C++ solution I just used what @username_17 posted (everything I tried writing took the same amount of time).
  
  Repeating 100 times and taking the average (yes, one should take the median):
  
  1. 1.7 seconds
  2. 0.65 seconds
  3. 15 seconds
  4. 1.2 seconds
  5. 0.48 seconds
  
  Unsurprisingly, anything in straight C or C++ is going to win. `grep` with `tr` is surprisingly good, which surprises me since even though `grep` is very very optimized I still expected it to have more overhead. Piping to `grep -io` is a performance killer. The R solution is surprisingly good given the typical R overhead.
  
  **Update1:** As suggested in the comments
  
  6. `sum(x.count('A') for x in open('seqs.fa','r') if x[0] != '>')` in python
  
  This yields a time of
  
  6. 1.6 seconds (I'm at work now, so I've adjusted the time as best I can to account for the different computers)
  
  **Update2:** More benchmarks from the comments
  
  7. `awk -F A '!/^>"/ {cnt+=NF-1}END{print cnt}' foo.fa`
  8. `perl -ne 'if(!/^>/){$count += tr/A//} END{print "$count\n"}' foo.fa`
  
  These yield times of:
  
  7. 5.2 seconds
  8. 1.18 seconds
  
  The awk version is the best hack I could come up with, there are likely better ways.
  
  **Update 3**: Right, so regarding fastq files, I'm running this on a 3.3GB gzipped file with 10 repetitions (this takes a bit of time to run), so I'm not going to initially limit tests to commands that I can trivially modify to handle compressed files (after all, uncompressed fastq files are an abomination).
  
  1. `seqtk fqchk`
  2. `sed -n '2~4p' <(zcat foo.fastq.gz) | grep -io A | wc -l`
  3. `sed -n '1d;N;N;N;P;d' <(zcat Undetermined_S0_R1_001.fastq.gz) | tr -cd A | wc -c`
  4. The C example from user1728181 (the comparison to C++ is already provided there, so no need to include it).
  5. `bioawk -c fastx '{n+=gsub(/A/, "", $seq)} END {print n}' foo.fastq.gz`
  
  The average times are:
  
  1. 1 minute 44 seconds
  2. 2 minutes 6 seconds
  3. 1 minute 52 seconds
  4. 1 minute 15 seconds
  5. 3 minutes 8 seconds
  
  N.B., these will often not handle fastq files with entries spanning more than 4 lines. However, such files are the dodo birds of bioinformatics, so that's not practically important. Also, I couldn't get `ShortRead` to work with the compressed fastq file. I suspect there's a way to fix that.
  
  > **Comments:**
  > **username_1**: Nice benchmark, what is the size of your test fasta file? Also, why not add Python if you have time, yes it will not be fastest, but interesting to see for comparison, I would try: `sum(x.count('A') for x in open('seqs.fa','r') if x[0] != '>')`, probably `BioPython` will be slower for this  
  > **username_2**: @username_1 I used dm6 (uncompressed), so it's not terribly large. I figured that python would be so self-evidently slower that I didn't want to bother, but I suppose it wouldn't hurt to add.  
  > **username_2**: @username_1 I've updated with a python benchmark. It does as well as R, which is better than I would have guessed.  
  > **username_1**: Thanks! I'm not so surprised since python's `str.count` is implemented in C. I also feel there should be a fast `awk` only solution, but I'm not able to figure it out  
  > **username_3**: Perl: ```perl -ne 'if(!/^>/){$count +=  tr/A//} END{print "$count\n"}' dm6.fa```  
  > **username_3**: Okay, now that the question has been updated, you can modify these to work for the general FASTQ case (multi-line sequence + quality string, with '@' possibly appearing in the quality line(s)).  
  > **username_2**: My awk-fu is sufficient to come up with a working example, but not good enough to come up with one that's not painfully slow.  
  > **username_2**: @username_3 There goes me being productive today...  
  > **username_4**: Can you try my `bioawk -c fastx '{n+=gsub(/A/, "", $seq)} END {print n}'` ? I'm curious to see how it compares with the other solutions on your data.  
  > **username_2**: @username_4: Yup, I have 5 others running at the moment...  
  > **username_5**: The R code can count Ns once you write it correctly. ;-) (My bad)  
  > **username_4**: The awk and perl versions do not give the same results on my test fasta file: The awk version find one more "A" in C. elegans genome than the perl version. I guess it depends whether one of the sequence starts or ends with "A".  
  > **username_2**: I've added a fastq comparison with some methods (including bioawk).  
  > **username_2**: @username_4 There should just be a `-1` in the awk code. Among the many reasons to use `bioawk` instead!  
  > **username_4**: I think the zcat is not necessary with bioawk. Does it improve the speed?  
  > **username_2**: @username_4 It's not needed and that's a mistake on my part (I didn't actually use zcat in the benchmark). Thanks for catching that!  
  > **username_6**: Why are you counting the `A`s and not the `N`s? Or am I missing something?  
  > **username_2**: Because the original ShortRead variant couldn't handle N and the performance comparisons are the same regardless. "The R code can't actually count N records for some reason when I try it, so I counted A for everything."  
  > **username_6**: Yes, sorry, I missed that when I first read the answer.  
  > **username_7**: You say: "Repeating 100 times and taking the average (yes, one should take the median):" But the average is the mean not median - so which one did you use? I appreciate the median is less susceptible to skew but in this case surely run time is normally distributed (in which case mean is fine) unless you're running on a system with very high instantaneous CPU or IO load from other users which is creating contention.  
  > **username_2**: The mean, not the median, thus the caveat. Yes, run time should be normally distributed, but occasional weird IO spikes and such do happen.  
  > **username_8**: Did you flush system cache between all these tests?  
  > **username_2**: Nah, in practice the files were in the cache for all of the runs.  
  > **username_9**: You should flush cache between tests. Your numbers will be misleading, otherwise, (well, wrong, really) regardless of whether you draw a mean or a median from test results.  
  
  ---
  
  ### Answer by username_10 on 2017-06-01. Score: 10
  
  I think that this should be pretty fast:
  
  FASTA:
  
  ```
  grep -v "^>" seqs.fa | tr -cd N | wc -c
  
  ```
  
  FASTQ:
  
  ```
  sed -n '1d;N;N;N;P;d' seqs.fq | tr -cd N | wc -c
  
  ```
  
  See this [answer on SO](https://stackoverflow.com/a/16679640/4641846) about how to count characters in BASH using different approaches.
  
  > **Comments:**
  > **username_5**: Your FASTQ answer unfortunately has the same problem, fundamentally, as username_11’s. That is, FASTQ records have a variable number of lines.  
  > **username_10**: I know that the original standard admits FASTQ line splitting. However, have you ever seen such a file in practice?  
  > **username_5**: I haven’t, but I’ve never myself worked with Sanger sequencing/PacBio data. I think you might encounter them there. For short read data I guess the “record = 4 lines” assumption is safe.  
  > **username_0**: Long read data do not use fastq. Currently PacBio uses bam (yes, really bam) and ONT fast5 (h5 archive). I think it is very reasonable to assume 4 line fastq.  
  
  ---
  
  ### Answer by username_11 on 2017-06-01. Score: 8
  
  **FASTQ**
  
  As it was pointed out, fastq can be complicated. But in a simple case when you have four lines per record, one possible solution in bash is:
  
  ```
  sed -n '2~4p' seqs.fastq | grep -io N | wc -l
  
  ```
  
  * `sed -n '2~4p'` will print every fourth line
  * `grep -o N` will output a line with `N` for every matching symbol
  * `wc -l` will count the lines
  
  I suspect this python approach will work faster:
  
  ```
  cat seqs.fastq | python3 -c "import sys; print(sum(line.upper().count('N') for line in sys.stdin))"
  
  ```
  
  **FASTA**
  
  Coreutils:
  
  ```
  grep -v "^>" seqs.fasta | grep -io N | wc -l
  
  ```
  
  Python:
  
  ```
  cat seqs.fasta | python3 -c "import sys; print(sum(line.upper().count('N') for line in sys.stdin if not line.startswith('>')))"
  
  ```
  
  > **Comments:**
  > **username_5**: Now you’re counting lines containing Ns, not Ns. Also, FASTQ records can span more than four lines. :-(  
  > **username_11**: @username_5 with -o grep outputs individual symbols, I will think how to correct fastq code  
  > **username_5**: FASTQ is tricky because `@` can also start the quality string line; see http://chat.stackexchange.com/transcript/message/37809200#37809200 and following lines. For reference, here’s a Perl script that parses FASTQ and trims the mate information off of the identifier. It can probably serve as a starting point; I think this is very close to the minimal code: https://gist.github.com/klmr/81a2b86cd93c706d28611f40bdfe2702  
  > **username_11**: @username_5, I updated the answer to indicate that this works only in the simple case  
  
  ---
  
  ### Answer by username_4 on 2017-06-02. Score: 6
  
  ## Using bioawk
  
  With [bioawk](https://github.com/lh3/bioawk) (counting "A" in the *C. elegans* genome, because there seem to be no "N" in this file), on my computer:
  
  ```
  $ time bioawk -c fastx '{n+=gsub(/A/, "", $seq)} END {print n}' genome.fa 
  32371810
  
  real    0m1.645s
  user    0m1.548s
  sys     0m0.088s
  
  ```
  
  bioawk is an extension of awk with convenient parsing options. For instance `-c fastx` makes the sequence available as `$seq` in fasta and fastq format (**including gzipped versions**), so the above command **should also handle fastq format robustly**.
  
  The `gsub` command is a normal awk function. It returns the number of substituted characters (got it from <https://unix.stackexchange.com/a/169575/55127>). I welcome comments about how to count inside awk in a less hackish way.
  
  On this fasta example it is almost 3 times slower than the `grep`, `tr`, `wc` pipeline:
  
  ```
  $ time grep -v "^>" genome.fa | tr -cd "A" | wc -c
  32371810
  
  real    0m0.609s
  user    0m0.744s
  sys     0m0.184s
  
  ```
  
  (I repeated the timings, and the above values seem representative)
  
  ## Using python
  
  ### readfq
  
  I tried the python readfq retrieved from the repository mentioned in this answer: <https://bioinformatics.stackexchange.com/a/365/292>
  
  ```
  $ time python -c "from sys import stdin; from readfq import readfq; print sum((seq.count('A') for _, seq, _ in readfq(stdin)))" < genome.fa
  32371810
  
  real    0m0.976s
  user    0m0.876s
  sys     0m0.100s
  
  ```
  
  ### "pure python"
  
  Comparing with something adapted from this answer: <https://bioinformatics.stackexchange.com/a/363/292>
  
  ```
  $ time python -c "from sys import stdin; print sum((line.count('A') for line in stdin if not line.startswith('>')))" < genome.fa
  32371810
  
  real    0m1.143s
  user    0m1.100s
  sys     0m0.040s
  
  ```
  
  (It is slower with python 3.6 than with python 2.7 on my computer)
  
  ### pyGATB
  
  I just learned (08/06/2017) that [GATB](http://gatb.inria.fr/) includes a fasta/fastq parser and has recently released a python API. I tried to use it yesterday to test another answer to the present question and [found a bug](https://github.com/GATB/pyGATB/issues/2). This bug is now fixed, so here is a [pyGATB](https://pypi.python.org/pypi/pyGATB/0.1.2)-based answer:
  
  ```
  $ time python3 -c "from gatb import Bank; print(sum((seq.sequence.count(b\"A\") for seq in Bank(\"genome.fa\"))))"
  32371810
  
  real    0m0.663s
  user    0m0.568s
  sys     0m0.092s
  
  ```
  
  (You can also do `sequence.decode("utf-8").count("A")` but this seems [a little slower](https://github.com/GATB/pyGATB/issues/2#issuecomment-307131176).)
  
  Although I used python3.6 here (pyGATB seems python3-only), this is faster than the other two python approaches (for which the reported timings are obtained with python 2.7). This is even almost as fast as the `grep`, `tr`, `wc` pipeline.
  
  ### Biopython
  
  And, to have even more comparisons, here is a solution using `SeqIO.parse` from Biopython (with python2.7):
  
  ```
  $ time python -c "from Bio import SeqIO; print(sum((rec.seq.count(\"A\") for rec in SeqIO.parse(\"genome.fa\", \"fasta\"))))"
  32371810
  
  real    0m1.632s
  user    0m1.532s
  sys     0m0.096s
  
  ```
  
  This is a bit slower than the "pure python" solution, but perhaps more robust.
  
  There seems to be a slight improvement with @peterjc's suggestion to use the lower level `SimpleFastaParser`:
  
  ```
  $ time python -c "from Bio.SeqIO.FastaIO import SimpleFastaParser; print(sum(seq.count('A') for title, seq in SimpleFastaParser(open('genome.fa'))))"
  32371810
  
  real    0m1.618s
  user    0m1.500s
  sys     0m0.116s
  
  ```
  
  (I did a series of timings and tried to take one that seemed representative, but there's a lot of overlap with the higher-level parser's timings.)
  
  ### scikit-bio
  
  I read about [scikit-bio](http://scikit-bio.org) today (23/06/2017), which has a sequence reader.
  
  ```
  $ time python3 -c  "import skbio; print(sum(seq.count('A') for seq in skbio.io.read('genome.fa', format='fasta', verify=False)))"
  32371810
  
  real    0m3.643s
  user    0m3.440s
  sys     0m1.228s
  
  ```
  
  ### pyfastx
  
  I read about [pyfastx](https://pyfastx.readthedocs.io/en/latest/index.html) today (27/08/2020), which has a sequence reader.
  
  ```
  $ time python3 -c "from pyfastx import Fasta; print(sum(seq.count('A') for (_, seq) in Fasta('genome.fa', build_index=False)))"
  32371810
  
  real    0m0.781s
  user    0m0.740s
  sys     0m0.040s
  
  ```
  
  (The timing is done using python 3.6, on the same machine as previous timings from 3 years ago, so it should be comparable.)
  
  ### `fastx_read` from mappy
  
  I realize (27/08/2020) that I hadn't added a benchmark for `fastx_read` from [mappy](https://pypi.org/project/mappy/) which I regularly use since quite some time already:
  
  ```
  $ time python3 -c "from mappy import fastx_read; print(sum(seq.count('A') for (_, seq, _) in fastx_read('genome.fa')))"
  32371810
  
  real    0m0.731s
  user    0m0.648s
  sys     0m0.080s
  
  ```
  
  ## Benchmarks on a fastq.gz file
  
  I tested some of the above solutions (or adaptations thereof), counting "N" in the same file that was used here: <https://bioinformatics.stackexchange.com/a/400/292>
  
  ### bioawk
  
  ```
  $ time bioawk -c fastx '{n+=gsub(/N/, "", $seq)} END {print n}' SRR077487_2.filt.fastq.gz
  306072
  
  real    1m9.686s
  user    1m9.376s
  sys     0m0.304s
  
  ```
  
  ### pigz + readfq python module
  
  readfq doesn't complain and is very fast when I pass directly the compressed fastq, but returns something wrong, so don't forget to manually take care of the decompression.
  
  Here I tried with `pigz`:
  
  ```
  $ time pigz -dc SRR077487_2.filt.fastq.gz | python -c "from sys import stdin; from readfq import readfq; print sum((seq.count('N') for _, seq, _ in readfq(stdin)))"
  306072
  
  real    0m52.347s
  user    1m40.716s
  sys     0m8.604s
  
  ```
  
  The computer has 16 cores, but I suspect the limiting factor for `pigz` is reading from the disk: the processors are very far from running full speed.
  
  And with `gzip`:
  
  ```
  $ time gzip -dc SRR077487_2.filt.fastq.gz | python -c "from sys import stdin; from readfq import readfq; print sum((seq.count('N') for _, seq, _ in readfq(stdin)))"
  306072
  
  real    0m49.448s
  user    1m31.984s
  sys     0m2.312s
  
  ```
  
  Here gzip and python both used a full processor. Resource usage balance was slightly better. I work on a desktop computer. I suppose a computing server would take advantage of `pigz` better.
  
  ### pyGATB
  
  ```
  $ time python3 -c "from gatb import Bank; print(sum((seq.sequence.count(b\"N\") for seq in Bank(\"SRR077487_2.filt.fastq.gz\"))))"
  306072
  
  real    0m46.784s
  user    0m46.404s
  sys     0m0.296s
  
  ```
  
  ### biopython
  
  ```
  $ time python -c "from gzip import open as gzopen; from Bio.SeqIO.QualityIO import FastqGeneralIterator; print(sum(seq.count('N') for (_, seq, _) in FastqGeneralIterator(gzopen('SRR077487_2.filt.fastq.gz'))))"
  306072
  
  real    3m18.103s
  user    3m17.676s
  sys     0m0.428s
  
  ```
  
  ### pyfastx (added 27/08/2020)
  
  ```
  $ time python3 -c "from pyfastx import Fastq; print(sum(seq.count('N') for (_, seq, _) in Fastq('SRR077487_2.filt.fastq.gz', build_index=False)))"
  306072
  
  real    0m51.140s
  user    0m50.784s
  sys     0m0.352s
  
  ```
  
  ### `fastx_read` from mappy (added 27/08/2020)
  
  ```
  $ time python3 -c "from mappy import fastx_read; print(sum(seq.count('N') for (_, seq, _) in fastx_read('SRR077487_2.filt.fastq.gz')))"
  306072
  
  real    0m52.336s
  user    0m52.024s
  sys     0m0.300s
  
  ```
  
  ## Conclusion
  
  pyGATB and bioawk both handle transparently compression (gzipped or not) and format differences (fasta or fastq). pyGATB is quite a new tool, but seems more efficient compared to the other python modules I tested.
  
  (Update 27/08/2020) pyfastx and mappy seem quite convenient too, and are only slightly slower than pyGATB.
  
  > **Comments:**
  > **username_12**: Could you add the Biopython low-level string FASTA parser numbers? ``time python -c "from Bio.SeqIO.FastaIO import SimpleFastaParser; print(sum(seq.count('A') for title, seq in SimpleFastaParser(open('genome.fa'))))"``  
  > **username_12**: Why do all your examples count the letter A, rather than N (and perhaps n) as per the question?  
  > **username_4**: I count letter "A" because the first large file I thought about for my testing happened to have no "N".  
  > **username_12**: I would think that all the Python examples will score badly on this small example file due to the overhead of loading Python and importing libraries. But if the user only wants to count one genome, that is a fair test all the same.  
  
  ---
  
  ### Answer by username_7 on 2017-06-02. Score: 6
  
  ### Decompression of gzipped FASTQ is the main issue
  
  If we take real world gzipped FASTQ (which as the OP suggested would be beneficial) rather than trivial FASTA as the starting point then the real issue is actually decompressing the file not counting the Ns and in this case the C program `count-N` is no longer the fastest solution.
  
  Additionally it would be good to use a specific file for benchmarking which actually has Ns, because you'll get some quite interesting execution time differences with some methods counting the more frequently occurring As rather than Ns.
  
  One such file is:
  
  <http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/sequence_read/SRR077487_2.filt.fastq.gz>
  
  It's also worth checking the various solutions return the correct answer there should be 306072 Ns in the above file.
  
  Next note that decompression of this file redirected to `/dev/null` is slower with `zcat` and `gzip` (which are both gzip 1.6 on my system) than say a parallel implementation of gzip like [Mark Adler's](http://zlib.net/pigz/) `pigz`, which appears to use 4 threads for decompression. All timings represent an average of 10 runs reporting real (wall clock time).
  
  ```
  time pigz -dc SRR077487_2.filt.fastq.gz > /dev/null
  
  real    0m29.0132s
  
  time gzip -dc SRR077487_2.filt.fastq.gz > /dev/null
  
  real    0m40.6996s
  
  ```
  
  There is an ~11.7 second difference between the two. Next if I then try to benchmark a one-liner which performs on FASTQ and gives the correct answer (Note I've yet to encounter FASTQ which is not 4 line, and seriously who generates these files!)
  
  ```
  time pigz -dc SRR077487_2.filt.fastq.gz | awk 'NR%4==2{print $1}' | tr -cd N | wc -c
  306072
  
  real    0m34.793s
  
  ```
  
  As you can see the counting adds ~5.8 second to the total run time versus `pigz` based decompression. Additionally this time delta is higher when using `gzip` ~6.7 seconds above gzip decompression alone.
  
  ```
  time gzip -dc SRR077487_2.filt.fastq.gz | awk 'NR%4==2{print $1}' | tr -cd N | wc -c                                                                
  306072
  
  real    0m44.399s
  
  ```
  
  The `pigz awk tr wc` based solution is however ~4.5 seconds faster than the `count-N` based C code solution:
  
  ```
  time count-N SRR077487_2.filt.fastq.gz 
  2385855128      306072  0
  
  real    0m39.266s
  
  ```
  
  This difference appears to be robust to re-running as many times as I like. I expect if you could use pthread in the C based solution or alter it to take the standard out from pigz it would also show an increase in performance.
  
  Benchmarking another alternative `pigz grep` variant appears to take more or less the same time as the `tr` based variant:
  
  ```
  time pigz -dc SRR077487_2.filt.fastq.gz | awk 'NR%4==2{print $1}' | grep -o N | wc -l
  306072
  
  real    0m34.869s
  
  ```
  
  Note that the `seqtk` based solution discussed above is noticeably slower,
  
  ```
  time seqtk comp  SRR077487_2.filt.fastq.gz | awk '{x+=$9}END{print x}'
  306072
  
  real    1m42.062s 
  
  ```
  
  However it's worth noting `seqtk comp` is doing a bit more than the other solutions.
  
  
  ---
  
  ### Answer by username_13 on 2017-06-01. Score: 5
  
  If it is raw speed you're after, then writing an own little C/C++ program is probably what you need to do.
  
  Fortunately, the worst part (a fast and reliable parser) has already been tackled: the [readfq](https://github.com/lh3/readfq) from Heng Li is probably the fastest FASTA/FASTQ parser around.
  
  And it's easy to use, the example on GitHub can easily be expanded to do what you need. Just add in a parameter parser (for the filename you want to analyse), program a simple 'N'-counter loop and have the results printed to stdout. Done.
  
  Else, if it is simplicity you are after, see username_5's answer for an R based solution. Or a very simple script using Biopython.
  
  > **Comments:**
  > **username_5**: God that’s absolutely atrocious code. ☹️ I think that for raw speed, [SeqAn](http://docs.seqan.de/seqan/master/?p=SeqFileIn) (C++) is probably on par, and a lot higher quality.  
  > **username_13**: It's not the prettiest around. Then again, it does what it needs to do and doesn't come with a whole library dependency. Do you have a benchmark with SeqAn & readfq?  
  > **username_4**: I'm curious to see how readfq ans SeqAn compare with the fasta/fastq parser from GATB: http://gatb-core.gforge.inria.fr/doc/api/classgatb_1_1core_1_1bank_1_1impl_1_1BankFasta.html#details I'm not really able to code in C/C++ though.  
  
  ---
  
  ### Answer by username_5 on 2017-06-01. Score: 5
  
  Honestly, the easiest way (especially for FASTQ) is probably to use a dedicated parsing library, such as R/Bioconductor:
  
  ```
  suppressMessages(library(ShortRead))
  
  seq = readFasta(commandArgs(TRUE)[1]) # or readFastq
  cat(colSums(alphabetFrequency(sread(seq))[, 'N', drop = FALSE]), '\n')
  
  ```
  
  This may not be the fastest, but it’s *pretty fast*, since the relevant methods are highly optimised. The biggest overhead is actually probably from loading `ShortRead`, which is an unjustifiable slog.
  
  
  ---
  
  ### Answer by username_14 on 2017-06-15. Score: 3
  
  For FASTA files, I've implemented a relatively efficient method in [pyfaidx](https://github.com/mdshw5/pyfaidx) v0.4.9.1. This post made me realize that my previous code was quite slow and easy to replace:
  
  ```
  $ pip install pyfaidx
  $ time faidx -i nucleotide ~/Downloads/hg38.fa
  name    start   end A   T   C   G   N   others
  chr1    1   248956422   67070277    67244164    48055043    48111528    18475410    
  chr10   1   133797422   38875926    39027555    27639505    27719976    534460  
  chr11   1   135086622   39286730    39361954    27903257    27981801    552880  
  chr11_KI270721v1_random 1   100316  18375   19887   31042   31012   0   
  ...
  chrX    1   156040895   46754807    46916701    30523780    30697741    1147866 
  chrY    1   57227415    7886192 7956168 5285789 5286894 30812372    
  chrY_KI270740v1_random  1   37240   8274    12978   7232    8756    0   
  
  real    2m28.251s
  user    2m15.548s
  sys 0m9.951s
  
  ```
  
  
  ---
  
  ### Answer by username_15 on 2018-03-01. Score: 0
  
  Here are perl 1-liners for fasta and fastq files that are fast enough if you just want totals:
  
  `$ perl -ne 'BEGIN { my $n=0 }; /^>/ or $n += s/N//gi; END { print "N=$n\n" }' contigs.fasta`
  
  For each line of a FASTA file, if not header (/^>/) then it's a sequence line.
  Substitute Ns, which returns the number of substitutions made, which is added to the total.
  
  `$ zcat reads.fastq | perl -ne 'BEGIN { my $n=0 }; /^[ATCGN]+$/i and $n += s/N//gi; END { print "N=$n\n" }'`
  
  For each sequence line of a FASTQ file, matched by /^[ATCGN]+$/i (i.e. all NTs), then count substitutions made in the same way.
  
  You may use pigz instead of zcat (or gunzip -c) if you have pigz installed.
  
  An uncompressed fastq file of 14M 250bp Illumina reads completed in 40sec.
  
  > **Comments:**
  > **username_16**: If you could explain what each perl line do, it would improve the answer, specially commenting the differences between both approaches  
  
  ---
  
  