Filtering FASTQ file based on read names from other file (how to increase performance) Python

Question

I have some code here that basically takes in a list of IDs (sequencing read names) using STDIN, and then filters fin (a gzipped FASTQ file) to just those reads with a read name which matches one of the IDs, and then writes the name and sequence to fout (a gzipped FASTA file). Fairly standard data processing problem in bioinformatics.

The code below can run through about 100k/sec lines of fin, however fin is about 500M lines, so it still takes over an hour.

Can this be sped up significantly? I'm hoping for a 10x or better improvement in performance. I may be taking a performance hit running line.decode('utf-8'), however it's necessary for matching since gzip files are read as bytestring. I'm willing to drop down and rewrite this in C if this will substantially increase speed.

import sys
import gzip
import re

id_array = []
for line in sys.stdin:
    id_array.append(line.rstrip())

ids = set(id_array)

with gzip.open(sys.argv[1], "r") as fin:
    with gzip.open(sys.argv[2], "w") as fout:

        take_next_line = False
        for line in fin.readlines():
            seq_search = re.search('[@+]([^\s]+)', line.decode('utf-8'))

            if take_next_line:
                fout.write(line)
                take_next_line = False

            elif seq_search and seq_search.group(1) in ids:
                take_next_line = True
                fout.write('>' + line[1:])

Answer 1

Do not reinvent the wheel. There are bioinformatics tools that accomplish this task.

To extract reads from fastq files by IDs, use seqtk subseq.

Extract sequences with names in file name.lst, one sequence name per line:
seqtk subseq in.fq name.lst > out.fq

It works with fasta files as well.

Use conda install seqtk or conda create --name seqtk seqtk to install the seqtk package, which has other useful functionalities as well, and is very fast.

Alternatively, use BBMap filterbyname.sh, see docs:

By default, "filterbyname" discards reads with names in your name list, and keeps the rest. To include them and discard the others, do this:
filterbyname.sh in=003.fastq out=filter003.fq names=names003.txt include=t

SEE ALSO:
Extracting specific sequences from FASTQ using Seqtk
How To Extracting Fastq Sequence For Given Fastq Ids And Fastq File
How To Extract A Subset Of Reads In Fastq Using An Id List?

Answer 2

Opened Python file objects are directly iterable, so you almost never need to use readlines(), which reads the entire file and creates a list storing the whole thing. Usually one avoids readlines() to prevent out-of-memory errors, but the change might help with speed, so it's worth trying. Just iterate directly:

for line in fin:
    ...

A couple regex issues, not related to speed: (1) it's a good habit to create all regexes with r'...' strings to avoid unwanted interpretation of escape characters; and (2) you don't need [^\s], which does the same thing as \S

rgx = re.compile(r'[@+](\S+)')

I haven't used the gzip library very many times, but the documentation implies that it will perform the decoding for you. The library might be faster at decoding than your code, because it performs fewer decoding operations (entire blocks at a time, rather than line-by-line). Anyway, all of this speculative, but it might be worth reversing things as an experiment: let gzip decode, and you encode whatever you need to write.

That said, your regular expression isn't very complex and Python does support regex searches on bytes, provided that the regex and the target string are both bytes. That might eliminate the need for large-scale decoding; maybe you just need to decode seq_search.group(1) when you get a match. Or just pre-convert ids to be bytes as well, and live in the land of bytes!

rgx = re.compile(rb'[@+](\S+)')
for line in fin:
    seq_search = rgx.search(line)
    ...

Beyond that, you could try parallelism: partitioning the work into blocks within the file and handing those partitions off to different Python processes. These heavy file reading/writing scripts sometimes benefit from that kind of approach, but not always.