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 whitespacerec.description= the whole line including the idrec.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
descriptionstill contains it. - Losing straddling features when slicing.
annotationsdropped by slicing, then GenBank writing fails on missingmolecule_type.- Setting
seqwhile oldletter_annotationsare attached. - Forgetting
qualifiersvalues 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.parseis 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