12
$\begingroup$

How do you write a .gz (or .bgz) fastq file using Biopython?

I'd rather avoid a separate system call.

The typical way to write an ASCII .fastq is done as follows:

for record in SeqIO.parse(fasta, "fasta"):
    SeqIO.write(record, fastq, "fastq")

The record is a SeqRecord object, fastq is the file handle, and "fastq" is the requested file format. The file format may be fastq, fasta, etc., but I do not see an option for .gz.

Here is the SeqIO API.

$\endgroup$
17
$\begingroup$

I'm not sure I'm doing it the best way, but here is an example where I read a compressed gzip fastq file and write the records in block gzip fastq:

from Bio import SeqIO, bgzf
# Used to convert the fastq stream into a file handle
from io import StringIO
from gzip import open as gzopen

records = SeqIO.parse(
    # There is actually simpler (thanks @peterjc)
    # StringIO(gzopen("random_10.fastq.gz").read().decode("utf-8")),
    gzopen("random_10.fastq.gz", "rt"),
    format="fastq")

with bgzf.BgzfWriter("test.fastq.bgz", "wb") as outgz:
    SeqIO.write(sequences=records, handle=outgz, format="fastq")
$\endgroup$
  • 3
    $\begingroup$ You shouldn't need the StringIO bit for parsing, just make sure you open the gzip handle in text mode under Python3: records = SeqIO.parse(gzopen("random_10.fastq.gz", "rt"), "fastq") $\endgroup$ – peterjc Jun 26 '17 at 16:31
  • $\begingroup$ @peterjc I didn't know about the "text mode", thanks. That explains the error message I got:ValueError: Is this handle in binary mode not text mode? $\endgroup$ – bli Jun 26 '17 at 17:15
2
$\begingroup$

Here is how I was able to export to a standard gzipped FASTQ file using Biopython. Basically, instead of using SeqIO.write, I directly called the .format method of the SeqRecord object. The example code below imports a gzipped FASTQ file, removes reads that do not contain a G in positions 7, 8, and 9, and writes the results to a gzipped FASTQ file.

from Bio import SeqIO
import gzip

path_in = "path/to/in.fastq.gz"
path_out = "path/to/out.fastq.gz"
handle_in = gzip.open(path_in, "rt")
handle_out = gzip.open(path_out, "wt")

fq = SeqIO.parse(handle_in, "fastq")
for read in fq:
    # Only export reads that have a G in positions 7, 8,
    # and 9
    if read.seq[6:9] == "GGG":
        handle_out.write(read.format("fastq"))

handle_in.close()
handle_out.close()

Using the code example from the original question that converts a FASTA file to a FASTQ file, it would look something like the following:

fastq = gzip.open("path/to/out.fastq.gz", "wt")
for record in SeqIO.parse(fasta, "fasta"):
    fastq.write(record.format("fastq"))
fastq.close()

Caveat: This only works for Python 3. Python 3 differs from Python 2 in how it imports data from binary files (source). In fact, ideally you should be using Python 3.4 or greater since Biopython has compatibility issues with Python 3.3 (source).

$\endgroup$
0
$\begingroup$

As far as I know you can not write compressed fastq on the fly. You will have to write the entire fastq and then compress it.

$\endgroup$
  • $\begingroup$ Hi Aaron Berlin, thanks for your answer and welcome to Bioinformatics stack exchange. Your answer is a bit short; do you have any source you can provide about not being able to write compressed fastq with Biopython, or an alternative python hook that could be used to work around the inability for it to support writing out a compressed file? $\endgroup$ – gringer Jun 24 '17 at 12:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.