Extracting A Subset Of Sequences From A Fastq File (Biopython Speed)

Question

[Darked89]

Initially my problem was to extract all entries from a FASTQ file with names not present in a FASTA file. Using biopython I wrote:

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 = []
for line in open(corrected_fn):
        if line[0] == ">":
                read_name = line.split()[0][1:] 
                corrected_names.append(read_name)

input_fastq_fn = uncorrected_fn
corrected_names.sort()

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()

Problem is, it is very slow. On 2Ghz workstation starting from a local disc it can take two days per pair of files:

• 4870868 seqs in FASTQ
• 4299464 seqs in FASTA

Removing title from corrected_names speeds up things a bit (this version I used for running).

Am I doing something obviously silly or simply FastqGeneralIterator is not a best construct to use here? While I like Python best, I am open to answers in Perl/Ruby.

Slicing and dicing FASTQ files based on lists seems to be fairly common task.

Edit: Python 2.6.4, biopython 1.53, Linux Fedora 8.

Edit 2:

• corrected one line of code, see comment to giovanni
• code snippet taken from: http://news.open-bio.org/news/2009/09/biopython-fast-fastq/

[Giovanni M Dall'Olio]

input_fastq_fn is not defined (in the 4th-to-last line)

[Darked89]

@giovanni: while cleaning up script for posting I cut one line to many (input_fastq_fn = uncorrected_fn) see also Edit 2 in the post

[Giovanni M Dall'Olio]

well, you don't need the input_fastq_fn variable anyway, you can use uncorrected_fn directly.

Answer 1

[Eric Normandeau]

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

[Will]

A set() might be even more pythonic. Replacing corrected_names[read_name] = 0 with correct_names.add(read_name])

[Peter]

As Will suggests, using a set rather than a dict would be more elegant and will waste less memory. Both the set and the dict will uses hashes to make the membership test MUCH faster than a list.

[Darked89]

Yes, it works ant it is a dramatic speed improvement. Resulting files are identical to ones produced by my slow-as-hell variant :-). Thank you a lot!

[Darked89]

One more thing: the old version was writing ca 100 sequences / minute (after creating the list, which was fairy fast). The new one writes almost 11 MILLIONS in less than 3 mins.

[Eric Normandeau]

Nice :) Python dictionaries are often a better choice than lists, mostly if execution time is limiting and you have to do multiple searches in the data structure.

[Eric Normandeau]

This is one useful thing I learned about Python thanks to you!

Answer 2

[Paulo Nuin]

A smallish change would be from

if line[0] == ">":

to

if line.startswith('>'):

[Eric Normandeau]

Sure is more elegant!

Answer 3

[brentp]

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.

[Peter]

That code ended up being called the Bio.SeqIO.index_db function, and will be in Biopython 1.57 due this month.

[brentp]

changed, thanks. and glad that's going in, it'll be quite useful.

[Nick Crawford]

Still not changed. The first line should be

fi = SeqIO.index_db(corrected_fn, corrected_fn + ".idx", "fastq")

[Nick Crawford]

The first line should be:

fi = SeqIO.index_db(corrected_fn + ".idx", corrected_fn, "fastq")

[brentp]

@Nick, indeed. updated.

Answer 4

[Giovanni M Dall'Olio]

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.

[Peter]

WARNING: You can get false positives using this:

grep '^@' my_input_fastq.ftq

This is because quality lines can start with an @ character.

[brentp]

cool. i didn't know about comm

[Darked89]

grep '^>' my_input_fasta.fas

will be even faster.

For FASTQ I will have to use:

grep '^@' my_input_fastq.ftq

Thanks for comm tip. I run away from (f)grep after getting a bit wrong results (names vary in length and can be substrings of other names)

[Darked89]

Also originally I wanted to get sequences for some reason omitted from the output by Correction for SOAPdenovo. So getting names alone will not do. BTW, Eric's solution is more than fast enough.

[Darked89]

@Peter: you are right that Illumina fastq standart contains @. Somehow quality strings in my files do not have @ or A, but even for a sequence with a string of Ns have B-quality instead. In general, yes, grep is way more dangerous than i.e. proper fastq format parsers.

[Giovanni M Dall'Olio]

it is not grep that gives you false positives, it is the regular expression you write. If that is wrong, it doesn't matter whether you match it with grep or with python's re. Anyway, sorry for the error in the post, I never used fastq files..

Answer 5

[Yifangt]

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

[Brad Chapman]

Welcome to BioStar. You should post this as a new question, instead of as an answer to an existing question.