grooming fastq dots to N's

[Giles]

I am trying to use various datasets from hudson alpha and the fastq files have dots in them. I don't know how to remove the dots. I'm not really a programmer, but I did try to remove them using grep from the command line. I thought it worked, because it got rid of some lines. But I didn't get rid of all of them. I think the problem has something to do with grooming them. Does anyone know how to get rid of these dots?

[arvid]

I am trying to use various datasets from hudson alpha and the fastq files have dots in them. I don't know how to remove the dots. I'm not really a programmer, but I did try to remove them using grep from the command line. I thought it worked, because it got rid of some lines. But I didn't get rid of all of them. I think the problem has something to do with grooming them. Does anyone know how to get rid of these dots?

What kind of reads are those? Are they short reads, and only have a single line per read, or do they span several lines per read? If they span several lines per read, I'd parse them with BioPython's SeqIO and do the replacement in there:

import sys
from Bio import Seq, SeqRecord, SeqIO

for record in SeqIO.parse(sys.stdin,"fastq-sanger"):
	SeqIO.write(SeqRecord.SeqRecord(seq=Seq.Seq(str(record.seq).replace(".","N")),letter_annotations=record.letter_annotations,id=record.id,name=record.name,description=record.description),sys.stdout,"fastq-sanger")

The above code can be used as follows (save the code in fastq_replace_dots_with_Ns_in_seq.py):
cat fastq_with_dots.fq | python fastq_replace_dots_with_Ns_in_seq.py > fastq_with_Ns.fq

If they are single line reads, an awk solution should be quite easily possible as well. A simple search/replace of "." will not work if you have a Sanger Phred encoded file, however, since "." is a valid quality score in a +33 encoding. Also, just simply dropping lines with dots will lead to a corrupted file - don't do that.
Anyways, beware that a IUPAC dot means an gap of unknown size, which is distinct from an N (any *one* base), so be careful with how the downstream alignment is done, especially if the reads come from a technology that delivers long gaps in the reads.

[maubp]

What kind of reads are those? Are they short reads, and only have a single line per read, or do they span several lines per read? If they span several lines per read, I'd parse them with BioPython's SeqIO and do the replacement in there:

If that turns out to be a bit slow, you could drop down to plain Python strings instead, e.g.

import sys
from Bio.SeqIO.QualityIO import FastqGeneralIterator
for title, seq, qual in FastqGeneralIterator(sys.stdin):
    print "@%s\n%s\n+\n%s" % (title, seq.replace(".", "N"), qual)

See also http://news.open-bio.org/news/2009/0...on-fast-fastq/