SeqIO

Core Principle: Bio.SeqIO is the unified sequence I/O interface across biological formats, offering stream parsing via parse(), singular record ingestion via read(), and dictionary-based indexing via index().

The four functions

from Bio import SeqIO

SeqIO.parse(path, fmt)        # → lazy iterator of SeqRecord
SeqIO.read(path, fmt)         # → ONE SeqRecord; raises if 0 or >1
SeqIO.write(records, path, fmt)   # → number written
SeqIO.convert(in_path, in_fmt, out_path, out_fmt)
SeqIO.index(path, fmt)        # → dict-like, random access, low memory
SeqIO.index_db(idx, files, fmt)   # → SQLite-backed index across many files
SeqIO.to_dict(SeqIO.parse(...))    # → dict, ALL records in memory

parse is lazy: and one-pass

Laziness is the feature: you can stream a 50 GB FASTQ in constant memory. But it means you get one pass. If you need two, re-open the file, or materialise it (list(...)) only when you know it fits.

n = sum(1 for _ in SeqIO.parse(path, "fasta"))     # count without loading
first = next(SeqIO.parse(path, "fasta"))            # just peek at one

read raises, deliberately

rec = SeqIO.read("single.gb", "genbank")     # ValueError if 0 or 2+ records

That exception is a feature: it validates your assumption that the file has exactly one record. Use read whenever you believe it does.

The formats you will use

Format string	What
"fasta"	FASTA
"fastq"	FASTQ, Sanger/Illumina 1.8+ (Phred+33)
"fastq-solexa", "fastq-illumina"	legacy encodings
"genbank" / "gb"	GenBank, with features
"embl"	EMBL
"swiss"	UniProt/SwissProt
"fasta-2line"	FASTA with no line wrapping
"tab"	two-column id/sequence
"abi"	Sanger chromatogram
"sff"	454
"nexus", "phylip", "clustal", "stockholm"	alignments: usually via AlignIO

Not supported: BAM, CRAM, VCF, BCF, bigWig. Biopython does not read alignment files. Use pysam.

Random access without loading

idx = SeqIO.index("genome.fasta", "fasta")
idx["chr17"]                     # a SeqRecord, parsed on demand
len(idx); "chr17" in idx; idx.keys()
idx.get_raw("chr17")             # raw bytes, no parsing: fastest
idx.close()

SeqIO.index stores only file offsets, so a 3 GB genome index costs a few megabytes. Records are parsed only when you ask for them.

For many files, or to persist the index across runs:

idx = SeqIO.index_db("genome.idx", ["chr1.fa", "chr2.fa"], "fasta")

The SQLite index is built once and reused, which matters for a pipeline that restarts.

Compare with SeqIO.to_dict(SeqIO.parse(...)), which loads everything into memory: fine for a few thousand short sequences, fatal for a genome.

For heavy random access to a reference genome, pysam.FastaFile with a standard .fai index is faster still and is what the rest of the ecosystem expects.

Compressed files

SeqIO takes a handle, so compression is just a different opener:

import gzip
with gzip.open("reads.fastq.gz", "rt") as fh:      # "rt" = text mode!
    for rec in SeqIO.parse(fh, "fastq"):
        ...

The "rt" matters: in binary mode you get bytes and a confusing parse error.

Writing

SeqIO.write(records, "out.fasta", "fasta")        # accepts a list OR a generator
SeqIO.convert("in.gb", "genbank", "out.fasta", "fasta")

with open("out.fasta", "w") as fh:
    SeqIO.write(records, fh, "fasta")

Passing a generator to write keeps the whole thing streaming:

def filtered(path, min_len=200):
    for rec in SeqIO.parse(path, "fasta"):
        if len(rec) >= min_len:
            yield rec

SeqIO.write(filtered("contigs.fasta"), "long.fasta", "fasta")

Constant memory regardless of file size. This generator-in, generator-out style is the idiomatic way to build Biopython pipelines.

Writing GenBank requires record.annotations["molecule_type"].

Format conversion

SeqIO.convert("in.fastq", "fastq", "out.fasta", "fasta")     # drop quality
SeqIO.convert("in.gb", "genbank", "out.fasta", "fasta")

convert is optimised: for some format pairs it copies records without fully parsing them.

Bioinformatics examples

# stream a large FASTQ, filter, write: constant memory
def qc_filter(path, min_mean_q=25, min_len=50):
    for rec in SeqIO.parse(path, "fastq"):
        q = rec.letter_annotations["phred_quality"]
        if len(rec) >= min_len and sum(q) / len(q) >= min_mean_q:
            yield rec

with gzip.open("in.fastq.gz", "rt") as fin, gzip.open("out.fastq.gz", "wt") as fout:
    n = SeqIO.write(qc_filter(fin), fout, "fastq")
print(f"kept {n:,} reads")

# extract a genomic region using an index
idx = SeqIO.index("hg38.fa", "fasta")
region = idx["chr17"].seq[7_668_401:7_687_550]     # 0-based! see [[Coordinate Systems]]

# split a multi-FASTA into per-sequence files
for rec in SeqIO.parse("all.fasta", "fasta"):
    SeqIO.write(rec, f"seqs/{rec.id}.fasta", "fasta")

# summarise a FASTA into a DataFrame
df = pd.DataFrame([
    {"id": r.id, "length": len(r), "gc": gc_fraction(r.seq)}
    for r in SeqIO.parse("contigs.fasta", "fasta")
])

# assembly N50
lengths = sorted((len(r) for r in SeqIO.parse("asm.fasta", "fasta")), reverse=True)
total = sum(lengths)
cum = 0
for L in lengths:
    cum += L
    if cum >= total / 2:
        print(f"N50 = {L:,}"); break

# subset by an id list
wanted = set(pd.read_csv("ids.txt", header=None)[0])
SeqIO.write((r for r in SeqIO.parse("all.fasta", "fasta") if r.id in wanted),
            "subset.fasta", "fasta")

Performance

Biopython's parsers are correct and general, not maximally fast. For tens of millions of reads:

• pyfastx: much faster FASTA/FASTQ with built-in random access
• pysam.FastxFile: fast, and already a dependency in most pipelines
• seqkit: a command-line tool; often the right answer for pure filtering/subsetting

Use Biopython when you need its object model, format breadth, or feature handling. Use a specialist when you are just streaming reads.

Common mistakes

• Iterating a parse result twice.
• SeqIO.to_dict on a genome: out of memory. Use index.
• gzip.open in binary mode. Needs "rt".
• Expecting BAM/VCF support. Use pysam.
• SeqIO.read on a multi-record file → ValueError (this is intended).
• Writing GenBank without molecule_type.
• list(SeqIO.parse(...)) reflexively, discarding the memory benefit.
• Wrong FASTQ variant, giving quality scores offset by 31.
• Forgetting to close() a SeqIO.index, leaving file handles open.

See also

SeqRecord · Seq · AlignIO · Applied - FASTA and FASTQ Workflows · Applied - Genomic File Formats · Coordinate Systems

Exercise

Biopython is not installed here, so parse the FASTQ yourself. That is worth doing once, because it exposes the trap the format carries: a quality line can begin with @, since @ is a legal Phred character, so a parser that splits on lines beginning with @ corrupts its own records. Read strictly in groups of four lines instead, and ignore a trailing group that is incomplete.

Write qc(offset, min_len, min_q, text) returning the report lines. Take the record id as the header text after @ up to the first whitespace. Skip any group whose header does not begin with @ or whose third line does not begin with +, and skip any record whose quality string length differs from its sequence length. Convert each quality character with ord(c) - offset: the offset is a parameter because the legacy Solexa and Illumina encodings differ from Sanger by 31, which silently shifts every score. Keep a record when its length is at least min_len and its mean quality is at least min_q. Print one <id> len=<n> meanq=<x> line per kept record in file order with the mean to two decimals, then a final kept=<k>/<total> line counting every well-formed record as the total. Empty input prints only kept=0/0.