0

I have many fastq files ending with fastq.gz.

rep1.fastq.gz
rep2.fastq.gz
rep3.fastq.gz
rep4.fastq.gz
.....

I expect my output as

rep1.fastq.gz 23516782
rep2.fastq.gz 45126780
rep3.fastq.gz 67543908
rep4.fastq.gz 76425368

Where row 1 show each of my input file and row 2 show the count of number of sequences in each file.

To achieve this I wrote a small bash script to count the number of sequences in each file with the number written after each file as an output

for sample in *.fastq.gz;do echo -en $sample "\t";(zcat $sample|wc -l)/4|bc ;done

I am getting an error: -bash: syntax error near unexpected token `/4'

1
  • Your error is because you are using a bc command out of bc: the shell has no idea that /4 is supposed to be "divide a number by 4". But are you sure this is what you want? There is no guarantee that fastq files will only have 4 lines per sequence. If you know your data well you can use this approach, but if you need to deal with arbitrary fastq files, you cannot assume only 4 lines.
    – terdon
    Jul 6, 2022 at 15:53

2 Answers 2

1

Ignoring all that can be wrong in the assumption of 4 lines per sequence... the command that you've shown above should work with the following format

for file in *.fastq.gz; do echo -en $file "\t";echo "$(zcat $file| wc -l)"/4 |bc;done
0

You are running this:

(zcat $sample|wc -l)/4|bc

The lone /4 there isn't understood by the shell, and this is why it fails. I suspect what you want is to get the result of the command zcat $sample|wc -l and then print that value along with /4 and pass that to bc. If so, you need $() and not just () and you also need to quote it:

echo "$(zcat $sample|wc -l)/4" | bc

So this would mean:

for sample in *.fastq.gz; do 
    echo -en $sample "\t"; echo "$(zcat $sample|wc -l)/4" | bc 
done

Or, a bit more portably:

for sample in *.fastq.gz; do 
    printf '%s\t%s\n' "$sample" "$(echo "$(zcat "$sample" | wc -l)/4" | bc )"
done

Alternatively, you can do the whole thing in awk:

for sample in *.fastq.gz; do 
    printf '%s\t' "$sample"
    zcat "$sample" | awk '!(NR % 4){k++}END{print k}'
done

However, do note that there is nothing in the definition of the fastq format that says that files will only have 4 lines per sequence. If you know your data well you can use this approach, but if you need to deal with arbitrary fastq files, you cannot assume only 4 lines and would be better off using a dedicated tool.

You might find this Q&A interesting: Fast way to count number of reads and number of bases in a fastq file?.

And also the FASTQ file format specification which clarified that you cannot assume only 4 lines per entry. That said, in my significant experience in working on human NGS data these past ~7 years now in a clinical setting, every single file I've seen only ever had 4 lines per sample. But I don't work with long read data and the format itself allows for more so it's something to consider.

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .