Or Kbps/Gbps. It feels like it should be conceptually very simple, but I can't seem to figure out the right combination of keywords to find it via my search engine. Help would be appreciated!
I have BBMAP, SRAtoolkit and MEGAHIT already installed, and also use bash. I'd be very happy if this can be answered with software that I already have, but if not that's perfectly fine.
I've been using this:
cat file.fastq | paste - - - - | cut -f 2 | tr -d '\n' | wc -c
Explanation :
paste - - - - : print four consecutive lines in one row (tab delimited), to merge the info for each read
cut -f2 : print only the second column, to access the sequence after the paste
wc -c : count the characters
tr -d '\n': to remove from count the eventual newline characters
(a tip for your googling: try search for "counting number of bases in fastq file")
The number of bases in a fastq file can be counted in bash with awk and wc
awk 'NR % 4 == 0' ORS="" fastqfile|wc -m
The awk code prints every fourth line (which is actually the quality scores, but that doesn't matter here). wc -m returns the total number of characters.
Edited to add ORS="" to prevent count of newlines (the double quotes are not strictly required)
Using Perl:
Print the sequence lines (line number 2, 6, 10, etc). Remove the newlines with chomp. Count the bytes (here, bases) using wc -c:
perl -ne 'if ( $. % 4 == 2 ) { chomp; print; }' | wc -c
Using a fastq parser in Python (for instance the one provided by mappy), it is rather straightforward to compute the sum of the lengths of the sequences present in a fastq file:
$ python3 -m pip install mappy # Only if mappy is not installed yet
$ python3 -c 'from mappy import fastx_read; print(sum(len(seq) for (_, seq, *_) in fastx_read("SRR077487_2.filt.fastq.gz")))'
2386161200
python3 -c '<insert some python code here>' uses the -c option of python3 that tells Python to execute the provided code (instead of starting an interactive interpreter, or running code present in a file).
Here the code has two steps (separated by a ;, in a standard script, we could use a new line instead): first making the fastx_read function available (because it is part of an optional module), second displaying the total length of the sequences.
fastx_read("path_to_a_file") generates (name, sequence, quality) triplets when a fastq file is provided, or (name, sequence) pairs when a fasta file is provided.
The (_, seq, *_) syntax is a form of "tuple unpacking" where we store the second element (the sequence) in a seq variable.
The _ in is just a way to ignore the first element in the triplet/pair. The *_ is a way to ignore whatever is after the second element in the triplet/pair.
len(seq) for (_, seq, *_) in fastx_read("SRR077487_2.filt.fastq.gz") is an expression representing the lengths (obtained with the predefined len function) of the sequences present in the file "SRR077487_2.filt.fastq.gz".
So the code uses the predefined sum function to compute the sum of values that are the lengths of the sequences extracted from the fastq file, and this is displayed using the print function.
This should also work with either fasta or fastq format, and with or without gz compression.
I tested more parsers on a slightly different problem here: https://bioinformatics.stackexchange.com/a/380/292 and you can find inspiration with the other approaches proposed in the other answers there.
seqkit stats fastq.gz. Get SeqKit from https://bioinf.shenwei.me/seqkit/usage/
With bbmap:
reformat.sh </path/to/fastq>
