I am new to next generation sequencing (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 that have 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 some suggestions or share any code?
Perl should be fine for this, I think the problem is trying to load the whole file in memory. It would help if you try to parse the file one line at a time, something like:
while () {
s/\//-/;
print $_;
}
This will ensure that perl is not trying to load the whole file into memory at once.
Perl should be fine for this, I think the problem is trying to load the whole file in memory. It would help if you try to parse the file one line at a time, something like:
while () {
s/\//-/;
print $_;
}
This will ensure that perl is not trying to load the whole file into memory at once.
Hello,
It depends of the change that you want to do in the file. The faster tools that you can use (without write your own C/C++ program) is to use a combination or Linux commands. For example, if you want to change the end of ID for a fastq file you can use:
sed -r 's/\/1$/-1/' myfastq.fq > out.fq
(for Illumina 1.3 or Illumina 1.5 where the pairs are @ID/1 or @ID/2)
(it changes /1 at the end of the line ($) for -1).
Keep in mind that the quality lines contains the characters '/' and '1' or '2' for Illumina 1.8+ (quality as Phred+33) but the Illumina 1.8+ uses a different way to show the pairs (@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG) where pair is ' 1'.
More information http://en.wikipedia.org/wiki/FASTQ_format
Good luck
Hi Shishir,
Just one tip, you have to be sure that the substitution is applied only to the sequence name line, because you could find those characters in the quality line.
In ruby it could be:
#!/usr/bin/env ruby
File.open("your_file.fastq").each_line do |line|
line.chomp!
line.strip!
if (line =~ /^@HWUSI-EAS/)
line.sub!(/\/1$/,'-1')
line.sub!(/\/2$/,'-2')
puts line
else
puts line
end
end
Hi,
you may do this and a lot of work more using pygr, which is a python module higly efficient and scalable to manipulate genomic data
http://code.google.com/p/pygr/
You only need a few lines of pyhon just to glue all together. Python is quite similar to Perl.
Have a look at Biostars (http://www.biostars.org/search/?q=fastq) or seqanswers. There is plenty of solution to similar problems...
for this particular case, use unix SED:
sed 's\/1$\-1\' read1.fastq > read1.corrected.fastq
sed 's\/2$\-2\' read2.fastq > read2.corrected.fastq
All of these are viable solutions; however in general I prefer tools that are designed for the task at hand. The stream editor 'sed' is designed to edit files of arbitrary length and in my experience, runs very quickly. The commands Leszek Pryszcz suggests seem to me the best solution.
Thanks Rogers.. I agree with you 'sed' works well, however it speeds down the cluster. Therefore, I have to move again for perl
Also see 'screed' in Python -- http://screed.readthedocs.org/en/latest/
Sed is very fast but could change the quality values line since you could have the characters /1 at the end of a quality line
You could give Biopython a try. Their SeqIO module works with many different bio-file formats, including every variant of FASTQ (http://biopython.org/wiki/SeqIO). PS: They released the latest version of Biopython (1.61) today (5-Feb-2013) and it has some updates for the SeqIO module.
I would recommend adding an extra line in your python code like:
if line.startswith("@"):
save.write(re.sub( r'/(\d)', r'-\1', line )
Otherwise it will run the regular expression on every line and this could really slow down the program. What happens to the new line character?
- Badges
- Science topic
- Similar topics
- Bioinformatics
- Bioinformatics and Computational Biology