Biopython Motif To Find Patterns In Sequences

I have a sequence file where for each read, I want to check for instances of a pattern(s). I will demultiplex the file into several files based on which pattern is found in each read. I am trying out the Biopython Motif module as -

motifs=Motif.Motif(alphabet=IUPAC.unambiguous_dna)  
motif.add_instance(Seq('ATAGCATAG',motifs.alphabet))
motif.add_instance(Seq('AAGCATAAG',motifs.alphabet))
motif.add_instance(Seq('CTAGCAGGG',motifs.alphabet))

etc. Then I loop over the file, and print out to the files as -

for read in SeqIO.parse(r_file,"fastq"):       
    for motif in motifs.search_instances(read.seq):                                                 
        motif_seq= str(motif[1])
        out_handle=dict_out_file_handles[motif_seq]
        out_handle.write(read.format("fastq"))

This is turning out to be rather slow. Is there are faster way to do this? Should I use plain regex?

Thanks very much.

biopython motif next-gen

• 10,279 views • link

written 13.9 years ago by dami.gupta •

You're doing a lot of IO during the loops! Perhaps you could read the FASTQ file into memory first, push any identified sequences with the motifs to a list and then output to disk?

• 0 views • link

written 13.9 years ago by Steve Moss

5 answers

Hi, have you checked

MotifMetrics -- Kitchen sink of routines for evaluating motifs on genomes This module is used both as a command-line script and as a repository for motif metrics. For a summary of the command-line usage, just type
python MotifMetrics.py

It may be able to help.

• 0 views • link

written 13.9 years ago by Whetting

The first rule about optimisation is to profile the code to find out where it is slow. However, based on experience I would guess that most of the time taken is overhead is reading FASTQ records into SeqRecord objects (including decoding the quality scores), and doing the reverse when writing them. You could stick with parsing the FASTQ records as strings instead, see http://news.open-bio.org/news/2009/09/biopython-fast-fastq/

• 1 views • link

written 13.9 years ago by Peter

Thanks - TAMO.MotifMetrics - not free :(

• 0 views • link

written 13.9 years ago by dami.gupta •

From: http://bioinformatics.oxfordjournals.org/content/21/14/3164.full.pdf

"Availability: TAMO is a Python/C++ package and requires Python 2.3 or higher. Source code and documentation are available at http://web.wi.mit.edu/fraenkel/TAMO/ maybe contact the corresponding author?

EDIT: http://fraenkel.mit.edu/TAMO/

• 0 views • link

updated 13.9 years ago by Istvan Albert

• written 13.9 years ago by Whetting

Further to my comment, something like this would remove the disk IO during the for loops:

reads = list(SeqIO.parse(r_file, "fastq"))

for i in range(len(reads)):
    for motif in motifs.search_instances(reads[i].seq):
        motif_seq = str(motif[1])
        outseqs.append([motif_seq, i])

for seq in outseqs:
    out_handle = dict_out_file_handles[seq[0]]
    out_handle.write(reads[seq[1]].format("fastq"))

• 1 views • link

written 13.9 years ago by Steve Moss

Not very memory friendly however :S

• 0 views • link

written 13.9 years ago by Steve Moss

My read file is huge. Creating a list would be too expensive

• 0 views • link

written 13.9 years ago by dami.gupta •

How huge is huge?

• 0 views • link

written 13.9 years ago by Steve Moss

I think that what is taking time is all_motifs.search_instances(read.seq) - which is the searching on the motif. I have taken that off, and am using

if bcode in str(rseq):
... and it is much faster

but, let me look at the FastqGeneralIterator now. Thanks everyone.

• 1 views • link

written 13.9 years ago by dami.gupta •