[Sunil Bhavsar]

I am a beginner for using Perl command. Wen I am trying FastaQual2fastq.pl script for making my fastq file but they like this error - readline() on closed filehandle.
So please help me to give a right solution for putting fasta seq.

[coolFlame]

Nice catch drio, thanks. One of those really subtle things you don't catch until you work with a different set of files.

Eugeni, sorry I didn't get back to you on this; got really crushed at work. I have uploaded a modified version of the script incorporating drio's fix.

@kmcarr: I found your script very useful and I am currently as a MSc Bioinformatics students working on an assignment which involves developing a web interface to a little mapping pipeline. This is purely for educational purposes. Would I be allowed to use your script to prepare the fastq file for the pipeline?
I really would appreciated it.

[SES]

Hi,
Is there a simple way to reduce the fasta name (e.g /
"> HH42GP401CAJLD length=118 xy=0823_0287 region=1 run=R_2012_01_27_13_59_03_ "

to ">HH42GP401CAJLD"?

I don't think you need a script for that. If your file is "454reads.fas" then just do:

sed 's/\s.*//' 454reads.fas > 454reads_trimmedheader.fas

[maubp]

Try something like this, untested:

from Bio import SeqIO

in_file = "example.fasta"
out_file = "new.fasta"
file_format = "fasta"

def remove_descr(record):
    record.description=""
    return record

#This is a generator expression - not all in memory at once!
wanted = (remove_descr(r) for r in SeqIO.parse(in_file, file_format))
count = SeqIO.write(wanted, out_file, file_format)

print "Saved %i records" % count

[lplough81]

how to trim FASTA name

Hi,
Is there a simple way to reduce the fasta name (e.g /
"> HH42GP401CAJLD length=118 xy=0823_0287 region=1 run=R_2012_01_27_13_59_03_ "

to ">HH42GP401CAJLD"?

Similar to trimming an SFF file to FASTA with biopython SeqIOconvert(), but taking a fasta file as the input and then outputting another fasta file?

Thanks,

Louis

[maubp]

Got it. Fairly new work for me, so I appreciate the patient replies. Can I specify the quality cutoff for trimming? Or what is the default that the biopython fastq trimmer uses?

There are two things to consider - getting rid of the adapter sequences and quality trimming. Roche does a good job of this as part of the base calling and production of the SFF file. When reading SFF files, Biopython (and other tools like sff_extract and Roche's own tools) will just apply the trimming information recorded in the SFF file. Using the Roche trimming is usually fine.

You may need to further trim off PCR primers or other library specific adapters if the Roche software wasn't told about them.

You may decide to further apply some quality cutoff trimming as well. This may be a good idea for some downstream analysis, not for others.

It is possible to do this kind of trimming in Biopython, but not in one line. There are some examples in the tutorial. I've written some SFF trimming tools using Biopython available within the Galaxy Tool Shed (if your institute runs its own Galaxy instance that may be interesting).

There are also other tools which will do it for you - especially if you want to work with the FASTQ file (or FASTA+QUAL) instead of the SFF file.

[lplough81]

OK!

Got it. Fairly new work for me, so I appreciate the patient replies. Can I specify the quality cutoff for trimming? Or what is the default that the biopython fastq trimmer uses?

Thanks again.

LP

[maubp]

Hi,
I was actually able to get it to run today.. Not sure what the problem was yesterday. But i got some funny results anyhow. Some of the nt's are uppercase and some are lowercase.

You'll see the same from Roche's own tools. The lower case are the bits which would be trimmed off as adapters or low quality bases.

This caused problems for some of the Galaxy fastx tools that summarize quality data.

Any thoughts?

That could be an oversight in fastx - ask them about it.

Or, what you probably want to do is ask for the trimmed sequences (which will be all upper case):

SeqIO.convert("454Reads.JA11255_155_RL13.sff", "sff-trim", "trimmed.fastq", "fastq")

[lplough81]

Hi,
I was actually able to get it to run today.. Not sure what the problem was yesterday. But i got some funny results anyhow. Some of the nt's are uppercase and some are lowercase. This caused problems for some of the Galaxy fastx tools that summarize quality data.

Any thoughts?

@HH42GP401CAJLD
gactagactcgacgtGTACTCAGGCTCGCACCGTGGCATGTCGCACTGTACTCAAGGCTCGCACCGTGGCATGTCGCACTGTACTTAAGGCTCACACCGTGGCATGTCGCACTGTACTCAAGGCACACAGGGGntaggnn
+
IIIIIIIIIIIIIIIIIIIGD666IIIIIIIIGDDDIIIIIIIIIIIIIIIGB;;;;IIIGGGGGCC>>>CIHID@@@C==:99==GGIIIIHIIIIIIIGGGCCCHIDDDC@777@C>1111AA@>;84445!;:44!!
@HH42GP401B4BC5
gactagactcgacgtGCAGTAGCTGCAATGGCGCAGAAGGCGTGCTTCtctctcncacgcacacacgagagagagngnnn
+
FFFFFFFFFFFFFFFIIIIIIIIIFFFFDDAAAB?<4444<>>9422323663/!//5///59=///2222////!2!!!

The code that I ran is here, (117,221 is the right number of reads for this file)
>>> SeqIO.convert("454Reads.JA11255_155_RL13.sff", "sff", "untrimmed.fastq", "fastq")
117221

[westerman]

Might help us if you demonstrated that the file is indeed not empty. How about a 'ls -l' on the file. Or an 'od -c yourfile.sff | head --lines 4' or the actual command you sent to SeqIO.convert so that we can be sure that you did send your file to it.

[lplough81]

Error on Fastq convert

HI,
I tried the fastq convert module in Biopython;

from Bio import SeqIO
SeqIO.convert("example.sff", "sff", "untrimmed.fastq", "fastq")

(I used my sff file though)

and I recieved this error:

File "/usr/lib/pymodules/python2.7/Bio/SeqIO/SffIO.py", line 258, in _sff_file_header
raise ValueError("Empty file.")
ValueError: Empty file.

Does this mean that there is an open line in the sff file? Any thoughts?

Thanks,
Louis

[maasha]

Using Biopieces you can do:

read_sff -i data.sff | write_fastq -o data.fq -x

or

read_sff -i data.sff | write_454 -o data.fna -q data.fna.qual -x

or both in one go:

read_sff -i data.sff | write_fastq -o data.fq | write_454 -o data.fna -q data.fna.qual -x

[prisnirath]

thanks...the script worked for me with a little alterations (minor ones).

[ketil]

Thanks for the benchmarks! What machine was used for this? I've written a program (flower - http://blog.malde.org/index.php/flower) to extract various information from SFF files, including Fasta and (Illumina or Sanger style) FastQ. It takes about 20 seconds to convert at 2.1G SFF to FastQ, but this is on a beefy server (Xeon 3.4GHz), so it's probably not directly comparable. Nice to see that we're in the same league, at least.

[nt2010]

Thanks BaCh and idas for your answers. All clear.

I'm not sure if i should continue here or start another thread. My questions would be that some of trimmed reads output by the converter(s) can still be very long with low quality at the end (Phred ~ 10). Should i trim then further, or it's acceptable to keep them as 454 works differently from illumina?