Biopython Databases

BioPython — Biological Databases

When to Use

• Fetching sequences/annotations from NCBI (Entrez), UniProt, PDB, or Ensembl by accession or gene name
• Reading/writing FASTA, GenBank, EMBL, or FASTQ files with Bio.SeqIO
• Manipulating sequences: reverse complement, transcription, translation, GC%, ORF finding
• Batch-downloading many records programmatically instead of one-by-one browser downloads
• Converting between sequence file formats (e.g. FASTQ → FASTA, GenBank → FASTA)

Version Compatibility

• Biopython ≥ 1.81, Python ≥ 3.9
• NCBI E-utilities (Entrez) REST API — no version string, but requires Entrez.email
• UniProt REST API at rest.uniprot.org (current since 2022; the old www.uniprot.org/uniprot/ endpoint is deprecated)
• RCSB PDB Data API (data.rcsb.org) and Search API v2 (search.rcsb.org/rcsbsearch/v2)
• Ensembl REST API at rest.ensembl.org (release-independent JSON endpoints)

Prerequisites

• pip install biopython
• stdlib urllib.request / json for the UniProt, PDB, and Ensembl REST calls (no extra dependencies needed)
• A registered Entrez.email; a free NCBI API key raises the rate limit from 3 to 10 requests/sec
• Familiarity with FASTA/GenBank/FASTQ file formats

Quick Reference

NCBI Accession Prefixes

SeqIO Format Strings

Key Patterns

Setup

from Bio.Seq import Seq, MutableSeq
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation
from Bio import SeqIO, Entrez, AlignIO, Align
from Bio.SeqUtils import gc_fraction, molecular_weight

Entrez.email = "your.email@example.com"  # required before any E-utility call

Seq object

dna = Seq("ATGCGATCGATCGTAA")
dna.complement()             # 3'->5' complement
dna.reverse_complement()     # 5'->3' complement of the other strand
dna.transcribe()             # DNA -> RNA (T -> U)
dna.translate()               # returns protein, * = stop codon
dna.translate(to_stop=True)  # stops before stop codon (most common)
dna.translate(table=2)       # mitochondrial genetic code

gc_fraction(dna) * 100       # GC percentage
molecular_weight(dna)        # DNA MW in Da
molecular_weight(protein, seq_type="protein")

m = MutableSeq("ATGCGATCG")
m[3] = "T"  # point mutation

SeqRecord object

record = SeqRecord(
    Seq("ATGCGATCG"),
    id="GENE_001",
    name="GENE",
    description="Description string"
)
record.annotations["organism"] = "Homo sapiens"
record.annotations["molecule_type"] = "DNA"  # required for GenBank output

feat = SeqFeature(FeatureLocation(0, 27), type="CDS",
                  qualifiers={"gene": ["BRCA1"], "product": ["BRCA1 protein"]})
record.features.append(feat)

record.letter_annotations["phred_quality"] = [30, 28, 35]
cds_seq = feat.location.extract(record.seq)

Goal: fetch a GenBank record by accession/gene name and translate its CDS.

Approach: Entrez.esearch to resolve a gene/query to accessions, Entrez.efetch with rettype="gb" to pull the full annotated record, then walk record.features for the CDS feature and use feature.location.extract() + .translate().

from Bio import Entrez, SeqIO

Entrez.email = "your.email@example.com"

def search_and_fetch_gene(term, db="nucleotide", retmax=5):
    """Resolve a free-text query to accessions, then fetch the top GenBank record."""
    handle = Entrez.esearch(db=db, term=term, retmax=retmax)
    results = Entrez.read(handle)
    handle.close()
    ids = results["IdList"]
    if not ids:
        return None
    handle = Entrez.efetch(db=db, id=ids[0], rettype="gb", retmode="text")
    record = SeqIO.read(handle, "genbank")
    handle.close()
    return record

record = search_and_fetch_gene(
    "insulin[Gene] AND Homo sapiens[Organism] AND mRNA[Filter] AND RefSeq[Filter]"
)
for feature in record.features:
    if feature.type == "CDS":
        protein = feature.location.extract(record.seq).translate(to_stop=True)
        print(feature.qualifiers.get("product", ["?"])[0], protein)

Batch fetch by accession list

accessions = ["NM_000518.5", "NM_000207.3", "NM_000546.6"]  # HBB, INS, TP53
handle = Entrez.efetch(db="nucleotide", id=",".join(accessions), rettype="fasta", retmode="text")
records = list(SeqIO.parse(handle, "fasta"))
handle.close()
SeqIO.write(records, "batch.fasta", "fasta")

elink: cross-database navigation

handle = Entrez.elink(dbfrom="nucleotide", db="protein", id="NM_000207.3")
link_results = Entrez.read(handle)
handle.close()
protein_ids = [link["Id"] for linkset in link_results
               for linkdb in linkset["LinkSetDb"]
               for link in linkdb["Link"]]

handle = Entrez.elink(dbfrom="gene", db="pubmed", id="7157")  # TP53

Goal: query UniProt, PDB, and Ensembl over their REST APIs.

Approach: these databases don't need Biopython — plain urllib.request + json against their public REST endpoints is simpler and has no extra dependency.

import json
import urllib.parse
import urllib.request


def fetch_uniprot(accession, fmt="json"):
    """Fetch a single UniProtKB entry by accession, e.g. fetch_uniprot('P01308') for human insulin."""
    url = f"https://rest.uniprot.org/uniprotkb/{accession}.{fmt}"
    with urllib.request.urlopen(url) as r:
        return json.loads(r.read()) if fmt == "json" else r.read().decode()


def search_uniprot(query, limit=5):
    """Full-text/field search against UniProtKB, e.g. 'insulin AND reviewed:true AND organism_id:9606'."""
    q = urllib.parse.quote(query)
    url = f"https://rest.uniprot.org/uniprotkb/search?query={q}&size={limit}&format=json"
    with urllib.request.urlopen(url) as r:
        return json.loads(r.read())


def get_pdb_info(pdb_id):
    """Fetch entry-level metadata (title, resolution, method) for a PDB ID, e.g. '1ATP'."""
    url = f"https://data.rcsb.org/rest/v1/core/entry/{pdb_id}"
    with urllib.request.urlopen(url) as r:
        return json.loads(r.read())


def download_pdb(pdb_id, fmt="pdb"):
    """Download a structure file (fmt='pdb' or 'cif') to the current directory."""
    url = f"https://files.rcsb.org/download/{pdb_id}.{fmt}"
    urllib.request.urlretrieve(url, f"{pdb_id}.{fmt}")


def ensembl_lookup(ensembl_id, expand=1):
    """Look up gene/transcript metadata by Ensembl stable ID, e.g. 'ENSG00000141510' (TP53)."""
    url = f"https://rest.ensembl.org/lookup/id/{ensembl_id}?expand={expand}"
    req = urllib.request.Request(url, headers={"Content-Type": "application/json"})
    with urllib.request.urlopen(req) as r:
        return json.loads(r.read())


def ensembl_sequence(ensembl_id, seq_type="cds"):
    """Fetch sequence (seq_type='genomic'|'cds'|'cdna'|'protein') for an Ensembl stable ID."""
    url = f"https://rest.ensembl.org/sequence/id/{ensembl_id}?type={seq_type}"
    req = urllib.request.Request(url, headers={"Content-Type": "application/json"})
    with urllib.request.urlopen(req) as r:
        return json.loads(r.read())

SeqIO: reading, writing, converting

records = list(SeqIO.parse("file.fasta", "fasta"))   # iterator -> wrap in list() to reuse
record = SeqIO.read("single.gb", "genbank")           # exactly one record

SeqIO.write(records, "out.fasta", "fasta")

d = SeqIO.to_dict(SeqIO.parse("file.fasta", "fasta"))
d["INS_HUMAN"].seq

count = SeqIO.convert("reads.fastq", "fastq", "reads.fasta", "fasta")  # quality info is lost

Pairwise alignment

aligner = Align.PairwiseAligner()
aligner.mode = "global"   # Needleman-Wunsch; "local" = Smith-Waterman
alignments = aligner.align(seq1, seq2)
print(alignments[0].score)

ORF finder (all 6 frames)

def find_orfs(seq, min_length=100):
    """Return (strand, frame, start, end, orf_seq) for every ORF >= min_length nt on both strands."""
    orfs = []
    for strand, s in [('+', seq), ('-', seq.reverse_complement())]:
        for frame in range(3):
            for i in range(frame, len(s) - 2, 3):
                if s[i:i+3] == 'ATG':
                    for j in range(i, len(s) - 2, 3):
                        if str(s[j:j+3]) in ('TAA', 'TAG', 'TGA'):
                            if j - i >= min_length:
                                orfs.append((strand, frame, i, j + 3, s[i:j + 3]))
                            break
    return orfs

Pitfalls

• Always set Entrez.email before any E-utility call; NCBI blocks unidentified clients.
• SeqIO.parse() returns an iterator — iterating twice yields nothing. Wrap in list() to reuse.
• SeqIO.read() raises if the file has 0 or >1 records — use parse() for multi-record files.
• record.annotations["molecule_type"] is required for GenBank output; omitting it raises ValueError.
• Use translate(to_stop=True) for clean protein sequences (omits the trailing *).
• efetch with rettype="gb" + SeqIO: use "genbank" (not "gb") as the SeqIO format string.
• Rate-limit NCBI: max 3 requests/sec without an API key, 10/sec with one. Add time.sleep(0.4) in loops.
• UniProt REST base URL changed in 2022: use rest.uniprot.org, not the deprecated www.uniprot.org/uniprot/.

See Also

• bio-database-access-entrez-search — building complex Entrez search queries
• bio-database-access-entrez-fetch — fetching records at scale with retries/rate limiting
• pdb-database — deeper PDB structure search and download workflows
• ensembl-database — Ensembl-specific comparative genomics and homology queries