prep_reads file in Tophat run shows a different number of reads

[gen2prot]

Hi All,

I have a very naive question. When I use "grep" on one of my fastq files to check the number of reads I get a certain value.

grep "@" SRR072905.fastq | wc -l
40042015

However, When I check the prep_reads.log file after a tophat run I get a different number of reads

more prep_reads.log 
prep_reads v1.1.4 (1709)
---------------------------
4060 out of 26976249 reads have been filtered out

40042015 are in the fastq file, but Tophat says the number is 26976249. Is the program not taking all the reads as input?

Please help so I know if I have to run the program again.

Thanks
Abhijit

[kmcarr]

Using grep "@" to count reads in a fastq file will not work the way you think it will.

Read this thread for an explanation.

[gen2prot]

But fastq is a set of four lines. The "@" line is followed by the sequence, while the "+" line is followed by the sequence quality. Therefore, counting the "@" or the "+" sign should give you the number of reads. Am I wrong?

[gen2prot]

Ignore my previous comment kmcarr. I figured out my mistake. I have to use the "^" at the beginning of the search while using "grep". Thnaks

[jasonwood]

The problem is that will flag every line with the @ character in it, so if your quality strings have that character you will double count. grep ^@ instead, or just count lines and divide by 4.

[kmcarr]

Originally Posted by gen2prot

Ignore my previous comment kmcarr. I figured out my mistake. I have to use the "^" at the beginning of the search while using "grep". Thnaks

Even using grep ^@ will not work perfectly. As the thread linked above notes, "@" is a valid quality character in Illumina FASTQ files which may even appear at the beginning of a quality line. grep ^@ will count these as well.

You need to use a search pattern for grep which will be absolutely unique to the header line. I typically use ^"@HW". The read IDs start with the machine names; our machine names all start with HW and "W" is not a valid quality character so the string @HW can only appear in read IDs.

[cram]

Couldn't you just do:

wc -l foo.fastq

and then divide by 4?

That should run faster too.

[nicolallias]

Originally Posted by jasonwood

just count lines and divide by 4.

Originally Posted by cram

wc -l foo.fastq

and then divide by 4?

Not if the fastq is generated with line wrapping - I've seen one of the bwa tool, qualfa2fq.pl doing this on the quality.

But you may overestimate the number of reads if the fastq is encoding quality Sanger's like (in quality line, "@" is available in Sanger's encoding ASCII+33)

You could try to export from fastq-sanger to fastq-illumina, do count with

grep -c ^@ yourfile.fastq

and then you're sure about the number of sequences.

A shorter way could be to grep some elements of the title (flowcell name, delimitation character(s)...).

Originally Posted by kmcarr

You need to use a search pattern for grep which will be absolutely unique to the header line.

Can you show us the head of your fastq ?