Counting specific sequence stretches in fastq files

[exo]

Hey all,

For a while ive been trying to find a way how to count specific sequence stretches directly from fastq files. For example i would like to see how many reads contain the following sequence: "TGACCCGATAGA". It seems to work with grep -c when i use fasta files but not fastq files.

Is there anything im missing?

Best,
Exo

[dpryan]

You can still use grep -c, but pipe the file through awk first:

awk '{if(NR%4 == 2) print $0}' sample.fastq | grep -c ...

[Brian Bushnell]

This is a perfect application for BBDuk. I would not recommend grep for this approach with multiline fasta as it does not distinguish between lines and reads, and cannot handle newlines in the middle of the target sequence. And I have no idea why grep would fail with fastq files; that should not be a problem. But anyway -

bbduk.sh in=reads.fq literal=TGACCCGATAGA k=12 mm=f rcomp=f int=f

That will tell you exactly how many reads contain that sequence. If you are also interested in reads containing the reverse-complement of the sequence, remove the flag "rcomp=f". The "int=f" flag tells it to treat reads individually rather than as pairs, if your reads are interleaved. k is the kmer length, which in this case is the length of your sequence, since you want an exact match. Incidentally, BBDuk can also tell you how many reads contain a sequence that is 1 substitution away from the target with the flag "hdist=1", and so forth. It can also handle degenerate symbols like 'N' if you add the flag "copyundefined", which is occasionally useful when dealing with amplification primers.

[rhinoceros]

You can still use grep -c, but pipe the file through awk first:

awk '{if(NR%4 == 2) print $0}' sample.fastq | grep -c ...

Returns at max 1 hit per sequence though..

[dpryan]

Returns at max 1 hit per sequence though..

Yes, though that's typically what's desired.

[dariober]

Returns at max 1 hit per sequence though..

Hi- Is that a problem? The OP asks for the number of reads containing pattern. So if a read contains the pattern twice it should still be counted as 1, which is what grep -c does, isn't it?