<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>

CC BY-SA 3.0
2017-06-02T07:53:08.817Z
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
Chris_Rands
CC BY-SA 3.0
2017-06-02T07:57:21.810Z
@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.
Devon Ryan
CC BY-SA 3.0
2017-06-02T08:08:24.853Z
@Chris_Rands I've updated with a python benchmark. It does as well as R, which is better than I would have guessed.
Devon Ryan
CC BY-SA 3.0
2017-06-02T08:14:11.577Z
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
Chris_Rands
CC BY-SA 3.0
2017-06-02T08:20:46.863Z
Perl: ```perl -ne 'if(!/^>/){$count +=  tr/A//} END{print "$count\n"}' dm6.fa```
gringer
CC BY-SA 3.0
2017-06-02T08:37:00.513Z
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)).
gringer
CC BY-SA 3.0
2017-06-02T08:37:18.690Z
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.
Devon Ryan
CC BY-SA 3.0
2017-06-02T08:38:13.573Z
@gringer There goes me being productive today...
Devon Ryan
CC BY-SA 3.0
2017-06-02T09:39:03.207Z
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.
bli
CC BY-SA 3.0
2017-06-02T09:39:43.203Z
@bli: Yup, I have 5 others running at the moment...
Devon Ryan
CC BY-SA 3.0
2017-06-02T10:29:14.590Z
The R code can count Ns once you write it correctly. ;-) (My bad)
Konrad Rudolph
CC BY-SA 3.0
2017-06-02T10:45:04.937Z
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".
bli
CC BY-SA 3.0
2017-06-02T11:34:17.793Z
I've added a fastq comparison with some methods (including bioawk).
Devon Ryan
CC BY-SA 3.0
2017-06-02T11:35:25.853Z
@bli There should just be a `-1` in the awk code. Among the many reasons to use `bioawk` instead!
Devon Ryan
CC BY-SA 3.0
2017-06-02T11:39:36.447Z
I think the zcat is not necessary with bioawk. Does it improve the speed?
bli
CC BY-SA 3.0
2017-06-02T11:41:12.797Z
@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!
Devon Ryan
CC BY-SA 3.0
2017-06-02T12:34:18.883Z
Why are you counting the `A`s and not the `N`s? Or am I missing something?
terdon
CC BY-SA 3.0
2017-06-02T12:35:28.177Z
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."
Devon Ryan
CC BY-SA 3.0
2017-06-02T13:04:55.853Z
Yes, sorry, I missed that when I first read the answer.
terdon
CC BY-SA 3.0
2017-06-02T21:07:06.653Z
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.
Matt Bashton
CC BY-SA 3.0
2017-06-02T21:20:41.933Z
The mean, not the median, thus the caveat. Yes, run time should be normally distributed, but occasional weird IO spikes and such do happen.
Devon Ryan
CC BY-SA 3.0
2017-06-03T02:14:44.960Z
Did you flush system cache between all these tests?
Torst
CC BY-SA 3.0
2017-06-03T06:28:14.220Z
Nah, in practice the files were in the cache for all of the runs.
Devon Ryan
CC BY-SA 3.0
2017-06-23T14:32:41.410Z
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.
Alex Reynolds
CC BY-SA 3.0
2017-06-01T22:20:21.293Z
2017-06-02T11:41:38.503Z
2017-06-02T11:41:38.503Z
Devon Ryan
<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>

CC BY-SA 3.0
2017-06-01T17:29:00.240Z
Your FASTQ answer unfortunately has the same problem, fundamentally, as Iakov’s. That is, FASTQ records have a variable number of lines.
Konrad Rudolph
CC BY-SA 3.0
2017-06-01T17:33:24.777Z
I know that the original standard admits FASTQ line splitting. However, have you ever seen such a file in practice?
Karel Břinda
CC BY-SA 3.0
2017-06-01T17:34:39.940Z
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.
Konrad Rudolph
CC BY-SA 3.0
2017-06-01T22:45:51.213Z
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.
Kamil S Jaron
CC BY-SA 3.0
2017-06-01T17:22:20.160Z
2017-06-01T17:22:20.160Z
Karel Břinda
<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>

CC BY-SA 3.0
2017-06-01T17:03:05.813Z
Now you’re counting lines containing Ns, not Ns. Also, FASTQ records can span more than four lines. :-(
Konrad Rudolph
CC BY-SA 3.0
2017-06-01T17:07:41.883Z
@KonradRudolph with -o grep outputs individual symbols, I will think how to correct fastq code
Iakov Davydov
CC BY-SA 3.0
2017-06-01T17:14:58.503Z
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
Konrad Rudolph
CC BY-SA 3.0
2017-06-01T18:11:57.120Z
@KonradRudolph, I updated the answer to indicate that this works only in the simple case
Iakov Davydov
CC BY-SA 3.0
2017-06-01T17:00:25.523Z
2017-06-02T13:41:57.900Z
2017-06-02T13:41:57.900Z
Iakov Davydov
<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>

CC BY-SA 3.0
2017-06-09T11:38:32.843Z
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'))))"``
Peter Cock
CC BY-SA 3.0
2017-06-09T11:39:05.597Z
Why do all your examples count the letter A, rather than N (and perhaps n) as per the question?
Peter Cock
CC BY-SA 3.0
2017-06-09T12:43:10.990Z
I count letter "A" because the first large file I thought about for my testing happened to have no "N".
bli
CC BY-SA 3.0
2017-06-10T21:59:36.417Z
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.
Peter Cock
CC BY-SA 4.0
2017-06-02T09:28:46.290Z
2020-08-27T15:33:15.750Z
2020-08-27T15:33:15.750Z
bli
<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>

CC BY-SA 3.0
2017-06-01T17:19:11.827Z
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.
Konrad Rudolph
CC BY-SA 3.0
2017-06-01T17:26:00.533Z
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?
BaCh
CC BY-SA 3.0
2017-06-08T11:57:15.183Z
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.
bli
CC BY-SA 3.0
2017-06-01T17:15:22.907Z
2017-06-01T17:15:22.907Z
BaCh
<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>

CC BY-SA 3.0
2017-06-01T17:10:02.993Z
2017-06-02T14:53:35.553Z
2017-06-02T14:53:35.553Z
Konrad Rudolph
<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>

CC BY-SA 3.0
2017-06-15T15:47:20.327Z
2017-06-15T15:47:20.327Z
Matt Shirley
<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>

CC BY-SA 3.0
2017-06-02T20:28:26.517Z
2017-06-02T20:57:56.203Z
2017-06-02T20:57:56.203Z
Matt Bashton
<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>

CC BY-SA 3.0
2018-03-01T09:02:34.747Z
If you could explain what each perl line do, it would improve the answer, specially commenting the differences between both approaches
llrs
CC BY-SA 3.0
2018-03-01T05:55:22.123Z
2018-03-01T23:15:30.270Z
2018-03-01T23:15:30.270Z
Edward Kirton
<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>

CC BY-SA 4.0
2021-04-21T09:35:25.373Z
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.
Konrad Rudolph
CC BY-SA 3.0
2017-06-02T08:34:23.513Z
based on Devon Ryan's benchmarking you are the fastest so far...
Kamil S Jaron
CC BY-SA 3.0
2017-06-02T08:54:52.737Z
Can you add `#include <inttypes.h>` to you first C example?
Devon Ryan
CC BY-SA 3.0
2017-06-02T10:33:44.673Z
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`.
Konrad Rudolph
CC BY-SA 3.0
2017-06-02T14:35:38.837Z
@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.
user172818
CC BY-SA 3.0
2017-06-02T15:01:26.313Z
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.
BaCh
CC BY-SA 3.0
2017-06-02T20:29:52.297Z
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`
Matt Bashton
CC BY-SA 3.0
2017-06-23T14:35:55.220Z
Are your numbers based on flushing the cache between samples? Otherwise, those numbers will almost certainly be incorrect.
Alex Reynolds
CC BY-SA 3.0
2017-06-01T21:46:22.750Z
2017-06-09T01:02:55.120Z
2017-06-09T01:02:55.120Z
user172818
<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>

CC BY-SA 3.0
2017-06-02T11:56:01.010Z
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
Chris_Rands
<fasta><fastq><benchmarking>
What is the fastest way to calculate the number of unknown nucleotides in FASTA / FASTQ files?
# 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:
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 to compile):
// to compile: gcc -O2 -o count-N this-prog.c -lz
#include 
#include 
#include 
#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 \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:
#include 

using namespace seqan;

int main(int argc, char *argv[]) {
    if (argc == 1) {
        std::cerr << "Usage: count-N " << 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): ~2.65 sec

> **Comments:**
> **username_0**: based on Devon Ryan'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**: @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.  
> **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 @user172818 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**: @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.  
> **username_2**: @Chris_Rands 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**: @gringer 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**: @bli: 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**: @bli 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**: @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!  
> **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 about how to count characters in BASH using different approaches.

> **Comments:**
> **username_5**: Your FASTQ answer unfortunately has the same problem, fundamentally, as Iakov’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**: @KonradRudolph 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**: @KonradRudolph, 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 (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 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. This bug is now fixed, so here is a pyGATB-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.)
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 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 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 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 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 Konrad'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 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  

---

