How To Write A Perl Script To Parse Fastq Files

I am undergrad learning Perl programming for bioinformatics, and having problems writing a script. I wanted to know if any one has script that counts the number of sequences in a fastq but excludes everything else such as the line that begins with @ and + and the quality score. Thanks

perl script fastq

• 10,018 views • link

updated 12.8 years ago by JacobS

• written 12.8 years ago by freddy •

First try it yourself. If it doesnt work, try these posts. They are pretty relevant to what you are trying to do Parse fastq file - pad reads with N's

Fastq Quality Read and Score Length Check

• 327 views • link

updated 6.6 years ago by Ram

• written 12.8 years ago by Ashutosh Pandey

  You can follow this approach:
  Step 1: use next if for skipping any line 
  Step 2: take a variable as string and concatenate every line to this strings then take the length of whole string.
  Always remember that perl always read file line by line

• 0 views • link

written 12.8 years ago by always_learning

4 answers

Due to the vagaries of the FastQ definition parsing it correctly is a little trickier than one might initially think. This is due to the fact that the symbol indicating the start of a new record @ is also a valid choice for quality measure. Thus one always needs to keep track of the lenght of the sequence and use that as guidance on how many quality measures to read in.

On the other hand most fastq files tend to be formatted with the entire sequence (and quality) on a single line. In these cases parsing is trivial as it turns into the problem of correctly identifying the 4 lines that form a fastq record. For example to find the sequences use a modulo division to find the remainder and print the line if the remainder is equal to 1 (assuming that the line numbering starts at 0).

All data produced by short read sequencers are in this latter format.

• 0 views • link

written 12.8 years ago by Istvan Albert

It depends what you want to do with your file, every sequence is (almost always) spread over 4 lines, so simply dividing number of rows by 4 will give you the result. If you want to access the sequences, you may want to open the file and read 4 lines (the second will be sequence of interest), then do some processing, and then again read 4 lines.. in a loop until you reach end of the file. You can find good manual how to read from file in perl here: http://www.perlfect.com/articles/perlfile.shtml

Here is a simple example how to read file:

open (MYFILE, 'file.fastq'); 
while (<MYFILE>) { 
     chomp;
   print "$_\n"; 
}
close (MYFILE);

You can modify it to count the number of sequences. However, for more complex task it is a good idea to take a look at Bioperl: http://www.bioperl.org/wiki/Main_Page

• 0 views • link

written 12.8 years ago by Biomonika (Noolean)

At the risk of sounding pedantic, I have to say that's not how you open a file in Perl (at least 3 things wrong in that first line) and that link you provided is outdated (see the perfunc page for open). The only reason you should ever use that form is to show people how it was done once upon a time, but it is very problematic and should never be used. There is really no debate on the matter (the Perl authors say they are trying to "rid the world" of this in the latest camel book) but it is still very prevalent because there are so many old books around and many webpages pop up with these guides (try to stick to the perldoc pages).

The modern perl books page has some good advice and some specific recommendations for Perl guides. Also, the Modern Perl book is free and the author (chromatic) maintains an active Perl blog that is very informative. I hope it doesn't seem like I'm nitpicking, it's just that I see this everyday on biostars and it bothers me, but for someone to say "do it this way" is not something I can ignore. Hopefully, this will save you some time. I talked to a couple of programmers at the Perl conference this summer and they told me they had just spent a month tracking down a bug in a large code base for their company, and it was related to this.

• 0 views • link

written 12.8 years ago by SES

You are right, I will delete my post. I personally don't use perl, so I got it wrong.

• 0 views • link

written 12.8 years ago by Biomonika (Noolean)

I made such a script here. You can see how I parsed it.

http://www.cbcb.umd.edu/software/PBcR/data/run_assembly_convertMultiFastqToStandard.pl

It converts multi-line fastq file entries to a four-line style.

• 0 views • link

written 12.8 years ago by Lee Katz

@freddy Here is my super simple answer for a beginner using perl to count the number of reads in a fastq file. Of course there are perl and awk one-liners that get the job done, but the following script really spells things out for you. All you are doing is reading through every 4 lines and increasing the total read count by 1. The following code can be copy and pasted to make a complete perl script:

#!/usr/bin/perl -w ## This specifies the file as a perl script

$line_position = 0;

$count = 0;

open(INPUT,$ARGV[0]) || die("Can't open file"); ## This opens the first user-privded argument as the input file.

while(<INPUT>) ## This reads each line of the file in, one at a time

{

$line_position++;

if($line_position == 4) ## If the script has already seen 4 lines, then reset the line counter and add 1 to the read counter!

{

$line_position = 0;

$count++;

}

close(INPUT);

print"Number of reads: $count\n"; ## This pring out your read counts!

An easier method would be to just do a UNIX line count on the file and divide by 4. wc -l myFile.fastq. This is all assuming you have a standard .fastq that has 4 lines per read. Good luck!

• 0 views • link

written 12.8 years ago by JacobS