FASTQ sequence converter

[Eugeni]

Hi everyone
Does anybody knows any program to convert a sff file from 454 data to fastq format incorporating information of sequence quality??; i have tried fq_all2std.pl from MAQ but does not incorporate information about sequence quality.
Thanks

[maubp]

The MAQ fq_all2std.pl script does not support SFF files.

The Roche off instrument applications are an "official" way, see:

http://seqanswers.com/forums/showthread.php?t=114

The open source Python tool sff_extract will do it too, it is used in MIRA:

http://bioinf.comav.upv.es/sff_extract/index.html

And in future I expect Biopython will offer this too.

Update: Biopython 1.54 onwards (April 2010) can read and write SFF files, and this includes converting them to FASTQ, FASTA, QUAL, etc.

http://news.open-bio.org/news/2010/05/biopython-release-154/

[BaCh]

The open source Python tool sff_extract will do it too, it is used in MIRA:

http://bioinf.comav.upv.es/sff_extract/index.html

Did Jose add a FASTQ option to sff_extract? Last version of sff_extract exported as FASTA + FASTA quality file (which could then be converted into a FASTQ using convert_project from the MIRA package).

I think I'll quickly add a FASTQ option to sff_extract if it's not already in ... I already played around with the idea for quite a while.

B.

[maubp]

Did Jose add a FASTQ option to sff_extract? Last version of sff_extract exported as FASTA + FASTA quality file (which could then be converted into a FASTQ using convert_project from the MIRA package).

I think I'll quickly add a FASTQ option to sff_extract if it's not already in ... I already played around with the idea for quite a while.

B.

Oh right - I assumed it was all done by sff_extract. But yes, adding FASTQ output would be great.

[kmcarr]

Here is a perl script to convert FASTA + QUAL files to FASTQ. You would need to first generate the FASTA and QUAL files from the SFF file using a tool like sffinfo from Roche or sff_extract.

#!/usr/bin/perl

use warnings;
use strict;
use File::Basename;

my $inFasta = $ARGV[0];
my $baseName = basename($inFasta, qw/.fasta .fna/);
my $inQual = $baseName . ".qual";
my $outFastq = $baseName . ".fastq";

my %seqs;

$/ = ">";

open (FASTA, "<$inFasta");
my $junk = (<FASTA>);

while (my $frecord = <FASTA>) {
	chomp $frecord;
	my ($fdef, @seqLines) = split /\n/, $frecord;
	my $seq = join '', @seqLines;
	$seqs{$fdef} = $seq;
}

close FASTA;

open (QUAL, "<$inQual");
$junk = <QUAL>;
open (FASTQ, ">$outFastq");

while (my $qrecord = <QUAL>) {
	chomp $qrecord;
	my ($qdef, @qualLines) = split /\n/, $qrecord;
	my $qualString = join ' ', @qualLines;
	my @quals = split / /, $qualString;
	print FASTQ "@","$qdef\n";
	print FASTQ "$seqs{$qdef}\n";
	print FASTQ "+\n";
	foreach my $qual (@quals) {
		print FASTQ chr($qual + 33);
	}
	print FASTQ "\n";
}

close QUAL;
close FASTQ;

Usage notes:

- Run the program just pass it the name of the fasta sequence file, e.g.

%> fastaQual2fastq.pl foo.fasta

(assuming you saved the above code with the name 'fastaQual2fastq.pl')

- The fasta filename must end in either .fasta or .fna

- The quality filename must have the same basename as the fasta file and end with .qual. For example, if your sequence file is "foo.fna" then the quality file must be named "foo.qual".

NOTE: Look downthread. There is a small bug in this script; a patched version is posted. Thanks to drio for finding the bug.

[maubp]

Here is a perl script to convert FASTA + QUAL files to FASTQ. You would need to first generate the FASTA and QUAL files from the SFF file using a tool like sffinfo from Roche or sff_extract.

If I'm reading that right (and I am not a Perl coder), you store the entire FASTA file in memory (as sequences keyed by ID), and then loop over the QUAL file. That isn't a great idea with large files.

I can provide a Python script using Biopython 1.51+ which would be memory efficient for comparison if you like

But anyway - having SFF to FASTQ handled by sff_extract would be very welcome

[kmcarr]

It was a conscious decision to use the hash in memory for the sequences. I understand that it would be much more efficient, memory wise, to stream the FASTA and QUAL files but there is big assumption to that design -- the order of the sequences in the FASTA and QUAL files must be identical. I wrote the script the way I did because I had a situation where the FASTA and QUAL files were not in the same order. By removing the assumption of identical order this design is more flexible.

Yes the hash takes memory but by todays standards not that much. It takes ~1.3X the size of your FASTA file. I tested it on a FASTA file from a full 454 Titanium run which contained just under 900,000 sequences. The FASTA file is 296 MB; the memory usage by the script topped out at ~ 400MB. 400MB of RAM is a pittance these days.

[maubp]

It was a conscious decision to use the hash in memory for the sequences. I understand that it would be much more efficient, memory wise, to stream the FASTA and QUAL files but there is big assumption to that design -- the order of the sequences in the FASTA and QUAL files must be identical. I wrote the script the way I did because I had a situation where the FASTA and QUAL files were not in the same order. By removing the assumption of identical order this design is more flexible.

Good point - although in the case of SFF output, this is a safe assumption to make. In fact, I would have the script check this assumption and verify the FASTA and QUAL files match up (same entries in same order with same sequence lengths).

How did you get FASTA and QUAL files which were out of sync?

Yes the hash takes memory but by todays standards not that much. It takes ~1.3X the size of your FASTA file. I tested it on a FASTA file from a full 454 Titanium run which contained just under 900,000 sequences. The FASTA file is 296 MB; the memory usage by the script topped out at ~ 400MB. 400MB of RAM is a pittance these days.

Fair enough - assuming the machine isn't doing anything else memory demanding at the same time.

[kmcarr]

Good point - although in the case of SFF output, this is a safe assumption to make.

Yes, tis true that the output from sffinfo or sff_extract will have the FASTA and QUAL file entries in the same order. If you can always count on that then by all means design your script around that.

How did you get FASTA and QUAL files which were out of sync?

The sequences were run through the SeqClean cleaning & trimming pipeline first (http://compbio.dfci.harvard.edu/tgi/software/). The final, cleaned FASTA and QUAL files are not matched in terms of order.

[maubp]

Yes, tis true that the output from sffinfo or sff_extract will have the FASTA and QUAL file entries in the same order. If you can always count on that then by all means design your script around that.

This is also a safe assumption for the Roche FASTA + QUAL files.

The sequences were run through the SeqClean cleaning & trimming pipeline first (http://compbio.dfci.harvard.edu/tgi/software/). The final, cleaned FASTA and QUAL files are not matched in terms of order.

Interesting - I wonder why they do that, and if it would be easy to fix their pipeline...

[kmcarr]

Interesting - I wonder why they do that, and if it would be easy to fix their pipeline...

The pipeline script (seqclean) is written in Perl so you could download it from the link above and check it out.

[Eugeni]

Here is a perl script to convert FASTA + QUAL files to FASTQ. You would need to first generate the FASTA and QUAL files from the SFF file using a tool like sffinfo from Roche or sff_extract.

#!/usr/bin/perl

use warnings;
use strict;
use File::Basename;

my $inFasta = $ARGV[0];
my $baseName = basename($inFasta, qw/.fasta .fna/);
my $inQual = $baseName . ".qual";
my $outFastq = $baseName . ".fastq";

my %seqs;

$/ = ">";

open (FASTA, "<$inFasta");
my $junk = (<FASTA>);

while (my $frecord = <FASTA>) {
	chomp $frecord;
	my ($fdef, @seqLines) = split /\n/, $frecord;
	my $seq = join '', @seqLines;
	$seqs{$fdef} = $seq;
}

close FASTA;

open (QUAL, "<$inQual");
$junk = <QUAL>;
open (FASTQ, ">$outFastq");

while (my $qrecord = <QUAL>) {
	chomp $qrecord;
	my ($qdef, @qualLines) = split /\n/, $qrecord;
	my $qualString = join ' ', @qualLines;
	my @quals = split / /, $qualString;
	print FASTQ "@","$qdef\n";
	print FASTQ "$seqs{$qdef}\n";
	print FASTQ "+\n";
	foreach my $qual (@quals) {
		print FASTQ chr($qual + 33);
	}
	print FASTQ "\n";
}

close QUAL;
close FASTQ;

Usage notes:

- Run the program just pass it the name of the fasta sequence file, e.g.

%> fastaQual2fastq.pl foo.fasta

(assuming you saved the above code with the name 'fastaQual2fastq.pl')

- The fasta filename must end in either .fasta or .fna

- The quality filename must have the same basename as the fasta file and end with .qual. For example, if your sequence file is "foo.fna" then the quality file must be named "foo.qual".

Hi, kmcarr
Thanks for you help, the script has been worked wery well, has generated the fastq file in the sanger format, although in the stdout of the script gives this message:
Argument "" isn't numeric in addition (+) at fastaQual2fastq.pl line 41, <QUAL> chunk 380185.
Dou you know what happens, if it is important?
Thanks a lot

[kmcarr]

Hi, kmcarr
Thanks for you help, the script has been worked wery well, has generated the fastq file in the sanger format, although in the stdout of the script gives this message:
Argument "" isn't numeric in addition (+) at fastaQual2fastq.pl line 41, <QUAL> chunk 380185.
Dou you know what happens, if it is important?
Thanks a lot

Does the warning only appear once? How many entries are in your FASTA/QUAL files?

[Eugeni]

Does the warning only appear once? How many entries are in your FASTA/QUAL files?

The warning appears associated to all sequences; i have 380185 fasta/qual entries

[maubp]

Just a guess, but you could check your line endings (DOS/Windows versus Unix).