pythonforbio.
[WASM idle]
Biopython14/15

SeqRecord

Starting sandbox…
Beginnerlesson

SeqRecord

In one line: a Seq plus everything else — identifier, description, annotations, features and per-letter data.

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:

rec = next(SeqIO.parse("reads.fastq", "fastq"))
q = rec.letter_annotations["phred_quality"]     # list of ints, one per base
min(q); sum(q) / len(q)

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:

trimmed = rec[10:90]        # both seq AND quality are sliced ✓

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

sub = rec[100:200]
Attribute Slicing behaviour
seq sliced
letter_annotations sliced
features only those fully within the slice, coordinates shifted
annotations dropped (they describe the whole record)
dbxrefs dropped
id, name, description kept

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

rc = rec.reverse_complement(
    id=rec.id + "_rc",
    description=rec.description,
    annotations=rec.annotations,       # NOT carried over by default
    features=True,                      # features are flipped correctly
)

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

for f in rec.features:
    if f.type == "CDS":
        print(f.location, f.qualifiers.get("gene"))
        cds = f.extract(rec.seq)              # handles joins and strand!
        protein = cds.translate(cds=True)

cds_feats = [f for f in rec.features if f.type == "CDS"]
gene_names = {f.qualifiers["gene"][0] for f in rec.features if "gene" in f.qualifiers}

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

lesson example

No output yet — run the code to populate this drawer.