How To Parse Fastq File?

Question

[SK]

I am new in NGS. I need to format some fastq files. For example I want to replace characters /1 and /2 with -1 and - 2 in the files having data like @HWUSI-EAS1671_0001:5:1:1022:10290/1. I tried to do it in perl with pattern matching but the files are too big to parse. Could anyone please give me suggestion on share any code?

[lh3]

Well, files are never too big to parse with Perl if you don't keep the file in RAM.

[Pierre Lindenbaum]

note: as far as I can see, all the examples below assume that it takes four lines for one FASTQ record .

[Obi Griffith]

This is generally the case with fastq parsers. I'd be curious if you have an alternative. One of the flaws of fastq in my opinion is that each record is spread over multiple lines with no guaranteed record separator

[Pierre Lindenbaum]

see

How common are multi-line fastq files?

[Lee Katz]

I made a tool to convert the nonstandard multi-line fastq files into 4-line entries. Fastq files should be 4 lines per entry because they are too difficult to parse otherwise. https://sourceforge.net/p/cg-pipeline/code/425/tree/cg_pipeline/branches/lkatz/scripts/run_assembly_convertMultiFastqToStandard.pl

[lh3]

While I agree multi-line fastq is inconvenient, it is not nonstandard as no official documentations invalidate multi-line fastq. Seqtk is probably the most efficient way so far to convert between multi-line and 4-line fastq.

[Alex Reynolds]

Are there cases where there are fewer lines for a FASTQ record? (I write quickie .fq parsers every now and again, so I'm genuinely curious.)

Answer 1

[Irsan]

Many times, you can use UNIX commands do file slicing and dicing tasks.

Consider the file old.fastq

@HWI-ST1019:196:D121WACXX:5:1101:1538:2300/1
CCGCGACCTCTGTTCTGCAGCCCCTTCCCTTCCCCGCCTCCTGCTCTGCCGGGACTACGCACCGGCCTGATTGGTTACCCCCGGGGTGTCCTCGGTCACCA
+
1+++4)<@<A<+2A9A2++:3C8:)1?BDBDBBDC@::6@(.8..7)777:<?@@##############################################
@HWI-ST1019:196:D121WACXX:5:1101:1538:2300/2
CCGCGACCTCTGTTCTGCAGCCCCTTCCCTTCCCCGCCTCCTGCTCTGCCGGGACTACGCACCGGCCTGATTGGTTACCCCCGGGGTGTCCTCGGTCACCA
+
1+++4)<@<A<+2A9A2++:3C8:)1?BDBDBBDC@::6@(.8..7)777:<?@@##############################################

you can use the stream editor (sed) to search and replace /1 for -1

sed 's:\/\([12]\):-\1:g' test.fastq

that will give you the output

@HWI-ST1019:196:D121WACXX:5:1101:1538:2300-1
CCGCGACCTCTGTTCTGCAGCCCCTTCCCTTCCCCGCCTCCTGCTCTGCCGGGACTACGCACCGGCCTGATTGGTTACCCCCGGGGTGTCCTCGGTCACCA
+
1+++4)<@<A<+2A9A2++:3C8:)1?BDBDBBDC@::6@(.8..7)777:<?@@##############################################
@HWI-ST1019:196:D121WACXX:5:1101:1538:2300-2
CCGCGACCTCTGTTCTGCAGCCCCTTCCCTTCCCCGCCTCCTGCTCTGCCGGGACTACGCACCGGCCTGATTGGTTACCCCCGGGGTGTCCTCGGTCACCA
+
1+++4)<@<A<+2A9A2++:3C8:)1?BDBDBBDC@::6@(.8..7)777:<?@@##############################################

[toni]

What happens here if you have '/1' in the quality strings ?

[Irsan]

Hmm, that also gets replaced: then use this to change any sequence id that ens with /1 to -1 assuming you have 1 sequence id every 4 lines

awk '{if(NR%4==1){gsub(/\/1$/,"-1");print $0}else{print $0}}' test.fastq

and for /2 to -2

awk '{if(NR%4==1){gsub(/\/2$/,"-2");print $0}else{print $0}}' test.fastq

[Biomonika (Noolean)]

I personally rather use part of the header, such as "HWI" when trying to match header than relying on 1 sequence always spread over 4 lines.

[Neilfws]

+1 for use of sed; see my blog post on handling fastq - http://nsaunders.wordpress.com/2011/12/22/sequencing-for-relics-from-the-sanger-era-part-1-getting-the-raw-data/

Answer 2

[JC]

As @lh3 pointed, you don't to read the full file in memory in Perl, just read line by line and process if required:

perl -plane 's/\/(\d)/-$1/ if (m/\@HWI/)' < file.fastq > new_file.fastq

Answer 3

[Alastair Kerr]

Some code golf for a Friday afternoon? Here is a perl one liner,

perl -pi.bak -e '/^@H/ && s#/([12])$#-\1#' filename

[Alastair Kerr]

beaten to the send button

Answer 4

[toni]

Hi, Whatever language you choose the files will remain big anyway, but I do not think it's worth to write a C code just for this. I do not know on which operating system you are working on, but supposing that you are under Linux (or equivalent) you can have a look at the sed command.

You must pay attention to the fact that the characters '/', '1' and '2' are in the range of those used for encoding quality. So when you do a substitution, you have to make sure that it is indeed a line of read name and not a quality string.

Take this example (file example.fasq):

@HWI-ST1019:196:/121WACXX:5:1101:1538:2300/1
CCGCGACCTCTGTTCTGCAGCCCCTTCCCTTCCCCGCCTCCTGCTCTGCCGGGACTACGCACCGGCCTGATTGGTTACCCCCGGGGTGTCCTCGGTCACCA
+
1+++4)<@<A<+2A9A2++:3C8:)1?BDBDBBDC@::6@(.8..7)777:<?@@######################################/1######

There are 3 '/1' appearances : 2 in the read name, one in the quality string, you only want the one at the end of the read name to be modified.

Then execute (assuming all your read names begin with HW):

sed "/^@HW/ s/\/1$/-1/g" example.fastq

And it produces :

@HWI-ST1019:196:/121WACXX:5:1101:1538:2300-1
CCGCGACCTCTGTTCTGCAGCCCCTTCCCTTCCCCGCCTCCTGCTCTGCCGGGACTACGCACCGGCCTGATTGGTTACCCCCGGGGTGTCCTCGGTCACCA
+
1+++4)<@<A<+2A9A2++:3C8:)1?BDBDBBDC@::6@(.8..7)777:<?@@######################################/1######

Is it REALLY necessary that you read such a huge file just to substitute a '/' with a '-' ?

Answer 5

[Nikolay Vyahhi]

In Python:

for line in open('file.fastq'):
  print line.replace('/1', '-1').replace('/2', '-2')

Answer 6

[Mikael Huss]

I find that modules like seqIO in BioPython are convenient for this kind of thing. Especially when you start needing to do other things than just change one character in each header. Of course there are analogues in HTSeq, BioPerl etc.