I can suggest a modification that should improve the speed issue, though I don't know by how much. Instead of making a list of your corrected_names, you could make a set and look if the values exists. Searching through a huge set for values will be much faster than sorting a list and searching through it. This is because a set is optimized for search speed.

It is, in addition, much more pythonic :)

from Bio.SeqIO.QualityIO import FastqGeneralIterator

corrected_fn   = "my_input_fasta.fas"
uncorrected_fn = "my_input_fastq.ftq"
output_fn      = "differences_fastq.ftq"

corrected_names = set() # Use a set instead of a list
for line in open(corrected_fn):
        if line[0] == ">":
                read_name = line.split()[0][1:] 
                corrected_names.add(read_name) # Add value to set

# corrected_names.sort() # No need, a set is orderless and optimized for search

handle = open(output_fn, "w")
for title, seq, qual in FastqGeneralIterator(open(input_fastq_fn)) :
        if title not in corrected_names:
                handle.write("@%s\n%s\n+\n%s\n" % (title, seq, qual))
handle.close()

Tell me if this improves things! (Note: I didn't test this as a whole)

Cheers

A smallish change would be from

if line[0] == ">":

to

if line.startswith('>'):

it's also worth noting that there are a couple good options for random access in fasta sequences. with a branch in biopython:

#fi = SeqIO.indexed_dict(corrected_fn, corrected_fn + ".idx", "fastq")
fi = SeqIO.indexed_db(corrected_fn + ".idx", corrected_fn, "fastq")
handle = open(output_fn, "w")
for title, seq, qual in FastqGeneralIterator(open(input_fastq_fn)):
        if title not in fi: continue
        handle.write("@%s\n%s\n+\n%s\n" % (title, seq, qual))
handle.close()

and also check out [?]screed[?] which does something very similar. these will both have overhead when first indexing the file but should be quite fast after that.

You can do it with Unix command line tools instead of python: it should be faster because these tools are written in C directly. First, use grep on the first file to get the list of all the ids:

$: grep '>' my_input_fasta.fas|sed 's/>//g'|sort > my_input_fasta.ids
$: grep '>' my_input_fastq.ftq|sed 's/>//g'|sort > my_input_fastq.ids

Then, use comm to get the differences:

$: comm fas_ids ftq_ids 
       seq1
seq2
    seq3
    seq4

The first column tells you the ids belonging to fas_ids, the second the ones belonging to ftq_ids, and the third the ones in common. You can use the -1, -2, -3 options to get only the ids belonging to one file, e.g.:

$: comm fas_ids ftq_ids  -23
seq2

will give you the ids present in fas_ids but not in the other.

My question is related, but not quite the same: I need to extract the substring of each sequence from a FASTQ file a sub-string of the sequence and the corresponding quality score, the rest untouched. My Illumina reads consists of 101bp long. I want remove the first 10bp and last 25bp then only the middle 66bp left.

The reason I need do this is my DNA sequence is methylated and the first 10 and last 25bp seem not having good quality in general so that I want get rid of them. My challenge is to remove both quality score and sequence correspondingly.

Seems quite easy, but I need to learn python above perl. Any tools there for me? Thanks a lot!

Yifang