SeqRecord

Core Principle: The SeqRecord object combines a primary sequence (Seq) with accession identifiers, descriptions, structured metadata dictionaries, and per-base feature annotations.

The anatomy

from Bio import SeqIO
rec = SeqIO.read("NM_000546.gb", "genbank")

rec.id                  # 'NM_000546.6'  : primary accession
rec.name                # 'NM_000546'    : locus name
rec.description         # 'Homo sapiens tumor protein p53 (TP53), mRNA'
rec.seq                 # Seq(...)
rec.annotations         # dict: molecule_type, organism, taxonomy, references, ...
rec.features            # list of SeqFeature: see [[SeqFeature and Location]]
rec.dbxrefs             # cross-references
rec.letter_annotations  # per-position data (quality scores!)
len(rec)                # length of the sequence

id vs name vs description

For a FASTA header >NM_000546.6 Homo sapiens tumor protein p53 (TP53), mRNA:

• rec.id = "NM_000546.6": everything up to the first whitespace
• rec.description = the whole line including the id
• rec.name = the id, for FASTA

description including the id catches people out when writing FASTA:

rec.description = ""        # ← otherwise the id appears twice in the output header
SeqIO.write(rec, "out.fasta", "fasta")

letter_annotations: per-base data

The FASTQ quality scores live here:

letter_annotations is a restricted dict: every value must be a sequence exactly as long as rec.seq. Biopython enforces this, so slicing keeps the correspondence:

That automatic co-slicing is the whole reason to work with SeqRecord rather than parallel lists.

Assignment order matters:

rec.seq = new_shorter_seq     # ValueError while old quality is still attached
# do this instead:
rec.letter_annotations = {}
rec.seq = new_seq
rec.letter_annotations["phred_quality"] = new_quals

Slicing behaviour

Features that straddle the boundary are silently discarded. If you slice out a region and find a gene has vanished, that is why.

Also note annotations being dropped means molecule_type disappears, which then breaks GenBank writing:

sub = rec[100:200]
sub.annotations["molecule_type"] = "DNA"      # restore it before writing
SeqIO.write(sub, "sub.gb", "genbank")

Building one

from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq

rec = SeqRecord(
    Seq("ATGGCC"),
    id="my_seq_1",
    name="my_seq_1",
    description="synthetic construct",
    annotations={"molecule_type": "DNA"},
)
SeqIO.write([rec], "out.fasta", "fasta")

molecule_type is required for GenBank/EMBL output and ignored for FASTA.

Reverse complementing a record

By default reverse_complement() drops the id, description and annotations, on the reasoning that they may no longer be accurate. You must opt in to keeping them. Feature locations are correctly reversed, including strand.

Working with features

f.extract(rec.seq) is the method that matters: it stitches together multi-exon join(...) locations and reverse-complements minus-strand features automatically. Doing this by hand with slicing is the classic source of off-by-one and strand bugs. See SeqFeature and Location.

Note f.qualifiers values are always lists, even when there is one value: hence f.qualifiers["gene"][0].

Bioinformatics examples

# quality-filter and trim reads
def process(path, min_q=20, trim=10):
    for rec in SeqIO.parse(path, "fastq"):
        q = rec.letter_annotations["phred_quality"]
        if sum(q) / len(q) < min_q:
            continue
        yield rec[trim:len(rec) - trim]        # seq and quality both trimmed

SeqIO.write(process("in.fastq"), "clean.fastq", "fastq")

# extract all CDS features as a protein FASTA
proteins = []
for rec in SeqIO.parse("genome.gb", "genbank"):
    for f in rec.features:
        if f.type != "CDS":
            continue
        gene = f.qualifiers.get("gene", ["unknown"])[0]
        proteins.append(SeqRecord(
            f.extract(rec.seq).translate(to_stop=True),
            id=gene, description=f"from {rec.id}",
        ))
SeqIO.write(proteins, "proteins.faa", "fasta")

# rename sequences from a mapping table
def rename(records, mapping):
    for rec in records:
        rec.id = mapping.get(rec.id, rec.id)
        rec.description = ""
        yield rec

# a summary table of a multi-FASTA
df = pd.DataFrame([
    {"id": r.id, "length": len(r), "gc": gc_fraction(r.seq),
     "n_count": str(r.seq).upper().count("N")}
    for r in SeqIO.parse("contigs.fasta", "fasta")
])

That last pattern: comprehension over SeqIO.parse into a DataFrame: is the standard bridge from Biopython into the pandas half of your analysis.

Common mistakes

• Duplicated id in FASTA output because description still contains it.
• Losing straddling features when slicing.
• annotations dropped by slicing, then GenBank writing fails on missing molecule_type.
• Setting seq while old letter_annotations are attached.
• Forgetting qualifiers values are lists.
• Manual slicing instead of feature.extract(), getting strand or joins wrong.
• Assuming reverse_complement() keeps the metadata. It does not.
• Holding every record in memory. SeqIO.parse is lazy; keep it that way. See SeqIO.

See also

Seq · SeqIO · SeqFeature and Location · Coordinate Systems · Applied - FASTA and FASTQ Workflows · Applied - Quality Control Plots

Exercise

Slicing a SeqRecord co-slices letter_annotations, shifts the features that fit and silently discards the ones that straddle the boundary. Biopython is not installed here, so reproduce that behaviour on plain Python objects, which is the clearest way to see which parts of a record survive a slice and which do not.

Write slice_record(start, end, seq, qual, features), where qual is a Phred+33 string or None and features is a list of (type, start, end) triples already in 0-based half-open coordinates. Return the report lines. Print seq=<sliced>, then meanq=<x> rounded to two decimals over the correspondingly sliced quality, or meanq=NA when there is no quality, when the slice is empty, or when the quality string length does not match the sequence length: a mismatched pair is exactly what Biopython refuses to construct. Then print one <type> <start> <end> line for each feature that lies entirely within the slice, in the original order and with its coordinates shifted to be relative to the slice, and finally kept=<k> dropped=<d>. A feature that overlaps the boundary in either direction is dropped. The bounds arrive as start end on the first line, the sequence on the second, the quality string or - on the third, and one type start end feature per line after that.