Bio Compressed Files

Version Compatibility

Reference examples tested with: BioPython 1.83+, htslib/bgzip 1.19+, samtools 1.19+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Compressed Files

Read, write, and randomly access gzip, bzip2, xz, and BGZF compressed sequence files.

"Read a compressed sequence file" -> Open a decompression handle in TEXT mode, then parse with the standard SeqIO interface.

• gzip: gzip.open(path, 'rt') (Python stdlib)
• bzip2: bz2.open(path, 'rt') (Python stdlib)
• xz/LZMA: lzma.open(path, 'rt') (Python stdlib)
• BGZF: bgzf.open(path, 'rt') (BioPython) - BGZF input ONLY

"Make a compressed file randomly accessible" -> Re-compress as BGZF, then index. Only BGZF supports SeqIO.index(), samtools faidx, and tabix on compressed data.

The Governing Principle: BGZF vs plain gzip is asymmetric

A BGZF (Blocked GNU Zip) file IS a valid gzip file - gunzip and zcat read it transparently. The reverse is FALSE: a plain .gz is NOT BGZF, so faidx/tabix/SeqIO.index reject it, and bgzf.open refuses to read it.

The reason is structural. BGZF is a series of concatenated gzip blocks, each <=64 KiB and independently decodable, so any record can be reached by seeking to its block. Plain gzip is one continuous DEFLATE stream with no block boundaries: reaching byte N means decompressing every byte before it (O(n)). Random access therefore REQUIRES BGZF; on a plain .gz, SeqIO.index() would re-decompress huge prefixes on every lookup, which is why Biopython forbids it outright (it raises rather than running slowly).

Consequences the agent must respect:

• Reading sequentially: gzip, bzip2, xz, and BGZF all work via the matching *.open(path, 'rt') handle.
• Random access / indexing: BGZF only. Convert plain gzip to BGZF first.
• bgzf.open() reads BGZF input only. Pointing it at a plain .gz raises ValueError: A BGZF block should start with b'\x1f\x8b\x08\x04'.... To read plain gzip use gzip.open().
• BAM and tabix-indexed files use BGZF natively; bzip2/xz are archive-only (no seekable index).

Required Imports

import gzip
import bz2
import lzma
from Bio import SeqIO
from Bio import bgzf

Reading Compressed Files

Goal: Parse sequence records from a compressed file without decompressing to disk.

Approach: Open a decompression handle in TEXT mode ('rt'), then pass the handle to SeqIO.parse(). The parser is format-agnostic about the underlying compression.

with gzip.open('reads.fastq.gz', 'rt') as handle:
    for record in SeqIO.parse(handle, 'fastq'):
        print(record.id, len(record.seq))

Swap gzip.open for bz2.open (.bz2), lzma.open (.xz), or bgzf.open (.bgz) - the parse loop is identical.

The 'rt' vs 'rb' trap

SeqIO.parse() in Python 3 needs a TEXT handle that yields str. A binary 'rb' handle yields bytes and raises TypeError: a bytes-like object is required (or a decode error). Always use 'rt' for reading and 'wt' for writing through SeqIO. The low-level SimpleFastaParser/FastqGeneralIterator also require text handles.

Writing Compressed Files

Goal: Save records straight to a compressed file with no intermediate plain copy.

Approach: Open a compression handle in TEXT mode ('wt'), then pass it to SeqIO.write().

with gzip.open('output.fasta.gz', 'wt') as handle:
    SeqIO.write(records, handle, 'fasta')

For an indexable result write BGZF instead:

with bgzf.open('output.fasta.bgz', 'wt') as handle:
    SeqIO.write(records, handle, 'fasta')

BgzfWriter.close() (and the with block exit) automatically appends the 28-byte empty-block EOF marker that htslib tools check for; let the context manager close the handle.

Random Access: index a BGZF file

Goal: Pull individual records by id from a large compressed file without a linear scan.

Approach: Compress as BGZF, then build a Biopython offset index. SeqIO.index() keeps virtual offsets in RAM; SeqIO.index_db() stores them in an on-disk SQLite index that persists across sessions and spans multiple files.

records = SeqIO.index('sequences.fasta.bgz', 'fasta')
target = records['gene_042'].seq
records.close()

# Persistent, multi-file, scales beyond RAM:
db = SeqIO.index_db('idx.sqlite', ['a.fasta.bgz', 'b.fasta.bgz'], 'fasta')

SeqIO.index() on a plain .gz raises ValueError: Gzipped files are not suitable for indexing, please use BGZF (blocked gzip format) instead. Convert first (below).

Convert plain gzip to BGZF

"Convert gzip to an indexable format" -> Decompress the gzip stream and re-compress it as BGZF.

CLI (fastest, htslib-native):

# Either decompress then bgzip in place...
gzip -d sequences.fasta.gz && bgzip sequences.fasta      # -> sequences.fasta.gz (now BGZF)
# ...or stream without touching disk:
zcat sequences.fasta.gz | bgzip -@ 4 > sequences.fasta.bgz

# Index a BGZF FASTA for region extraction:
samtools faidx sequences.fasta.bgz                        # writes BOTH .fai and .gzi
samtools faidx sequences.fasta.bgz gene_042:1-200

samtools faidx on a BGZF FASTA writes TWO index files: .fai (record offsets in uncompressed coordinates) AND .gzi (the compressed-to-uncompressed block map). Deleting .gzi breaks region extraction even though .fai survives. Note that bgzip keeps the .gz extension, so a .gz may be EITHER plain gzip or BGZF - check with bgzip -t file.gz (tests for a valid BGZF stream) rather than trusting the suffix.

Pure-Python equivalent (no external tools):

import gzip
from Bio import SeqIO, bgzf

with gzip.open('input.fasta.gz', 'rt') as in_handle:
    with bgzf.open('output.fasta.bgz', 'wt') as out_handle:
        SeqIO.write(SeqIO.parse(in_handle, 'fasta'), out_handle, 'fasta')