!pip install -q biopythonSession 3 — Biological Databases in code
This notebook builds the Session 3 page’s HBB database dossier with live code, following the page’s five questions: where is the sequence (NCBI Gene and RefSeq), what is known about the protein (UniProt, with its protein-existence level), which family it belongs to (InterPro), is there a structure (PDB entry 2HHB through two APIs, with its primary citation, and the AlphaFold model), and what it interacts with (IntAct, Complex Portal, STRING). It then counts Swiss-Prot and TrEMBL entries, checks the RefSeq prefix convention (including what a predicted XP_ record says about itself), and runs the Boolean PubMed searches from the page, checking the set arithmetic.
Everything talks to the real databases over the network. If you are offline, the cells will not run, but the saved outputs show what a live run produces.
Setup
from Bio import Entrez, SeqIO
import requests
# NCBI asks for a real contact email on every Entrez request -- not authentication,
# just so they can reach you if a script is hammering their servers by mistake.
Entrez.email = "kb8029-book@example.org"Gene → mRNA → protein: following HBB with real records
Start from the official gene symbol, HBB, and follow its cross-references down through the levels – the worked example on the Session 3 page.
# The gene record itself
handle = Entrez.esummary(db="gene", id="3043")
gene = Entrez.read(handle)["DocumentSummarySet"]["DocumentSummary"][0]
handle.close()
print("Gene symbol: ", gene["Name"])
print("Description: ", gene["Description"])
print("Chromosome: ", gene["Chromosome"], f"({gene['MapLocation']})")Gene symbol: HBB
Description: hemoglobin subunit beta
Chromosome: 11 (11p15.4)# The mRNA (RNA-level) record
handle = Entrez.efetch(db="nucleotide", id="NM_000518.5", rettype="gb", retmode="text")
mrna = SeqIO.read(handle, "genbank")
handle.close()
print("Accession: ", mrna.id)
print("Description: ", mrna.description)
print("Organism: ", mrna.annotations["organism"])
print("Length: ", len(mrna.seq), "bases")
print("First 60 bases:", str(mrna.seq[:60]))Accession: NM_000518.5
Description: Homo sapiens hemoglobin subunit beta (HBB), mRNA
Organism: Homo sapiens
Length: 628 bases
First 60 bases: ACATTTGCTTCTGACACAACTGTGTTCACTAGCAACCTCAAACAGACACCATGGTGCATC# The protein (protein-level) record -- and its cross-reference back to the mRNA above
handle = Entrez.efetch(db="protein", id="NP_000509.1", rettype="gb", retmode="text")
protein = SeqIO.read(handle, "genbank")
handle.close()
print("Accession: ", protein.id)
print("Description: ", protein.description)
print("Length: ", len(protein.seq), "residues")
print("Cross-reference (DBSOURCE):", protein.annotations["db_source"])
assert "NM_000518" in protein.annotations["db_source"], "the protein record should point back to the mRNA fetched above"Accession: NP_000509.1
Description: hemoglobin subunit beta [Homo sapiens]
Length: 147 residues
Cross-reference (DBSOURCE): REFSEQ: accession NM_000518.5The same protein, from UniProt’s side
RefSeq and UniProt are entirely separate databases, run by different organizations – so a real cross-reference between them, not just a shared name, is the only trustworthy way to know they describe the same protein.
r = requests.get("https://rest.uniprot.org/uniprotkb/P68871.json", timeout=10)
r.raise_for_status()
entry = r.json()
print("Entry type: ", entry["entryType"])
print("UniProt ID: ", entry["uniProtkbId"])
print("Protein name: ", entry["proteinDescription"]["recommendedName"]["fullName"]["value"])
print("Protein existence:", entry["proteinExistence"]) # graded 1 (protein level) to 5 (uncertain)
refseq_xrefs = [x for x in entry["uniProtKBCrossReferences"] if x["database"] == "RefSeq"]
for xref in refseq_xrefs:
mrna_ids = [p["value"] for p in xref["properties"] if p["key"] == "NucleotideSequenceId"]
print(f"RefSeq protein {xref['id']} <- cross-referenced mRNA {mrna_ids}")
all_ids = [x["id"] for x in refseq_xrefs] + [p["value"] for x in refseq_xrefs for p in x["properties"] if p["key"] == "NucleotideSequenceId"]
assert any("NM_000518" in x for x in all_ids), "UniProt's own cross-references should list the same mRNA accession"
assert any("NP_000509" in x for x in all_ids), "UniProt's own cross-references should list the same protein accession"Entry type: UniProtKB reviewed (Swiss-Prot)
UniProt ID: HBB_HUMAN
Protein name: Hemoglobin subunit beta
Protein existence: 1: Evidence at protein level
RefSeq protein NP_000509.1 <- cross-referenced mRNA ['NM_000518.5']Family and domain membership: InterPro
Which protein family/domain entries does P68871 (HBB) actually match, straight from InterPro’s own API?
r = requests.get(
"https://www.ebi.ac.uk/interpro/api/entry/interpro/protein/uniprot/P68871/",
timeout=30,
)
r.raise_for_status()
interpro_entries = r.json()["results"]
for entry in interpro_entries:
meta = entry["metadata"]
print(meta["accession"], "-", meta["name"])
assert any(e["metadata"]["accession"] == "IPR000971" for e in interpro_entries), \
"HBB should match the Globin family entry"IPR000971 - Globin
IPR002337 - Hemoglobin, beta-type
IPR009050 - Globin-like superfamily
IPR012292 - Globin/Protoglobin
IPR050056 - Hemoglobin and related oxygen transportersStructure: PDB entry 2HHB, two different APIs
PDBe’s REST API and RCSB’s GraphQL Data API, both queried live for the same real entry (2HHB, human deoxyhemoglobin) – and the AlphaFold DB prediction for the same protein, for comparison.
# PDBe's REST API (EBI)
r = requests.get("https://www.ebi.ac.uk/pdbe/api/pdb/entry/summary/2hhb", timeout=30)
r.raise_for_status()
pdbe_summary = r.json()["2hhb"][0]
print("Title: ", pdbe_summary["title"])
print("Experimental method: ", pdbe_summary["experimental_method"])
print("Release date: ", pdbe_summary["release_date"])
assert "DEOXYHAEMOGLOBIN" in pdbe_summary["title"].upper()Title: THE CRYSTAL STRUCTURE OF HUMAN DEOXYHAEMOGLOBIN AT 1.74 ANGSTROMS RESOLUTION
Experimental method: ['X-ray diffraction']
Release date: 19840718# RCSB's Data API (GraphQL) -- same entry, different API, ask for exactly one field
query = '{ entry(entry_id: "2HHB") { struct { title } rcsb_entry_info { resolution_combined } rcsb_primary_citation { pdbx_database_id_PubMed year journal_abbrev journal_volume page_first } } }'
r = requests.post("https://data.rcsb.org/graphql", json={"query": query}, timeout=30)
r.raise_for_status()
rcsb_data = r.json()["data"]["entry"]
print("Title: ", rcsb_data["struct"]["title"])
print("Resolution:", rcsb_data["rcsb_entry_info"]["resolution_combined"], "Angstrom")
cite = rcsb_data["rcsb_primary_citation"]
print(f"Primary citation: PMID {cite['pdbx_database_id_PubMed']} ({cite['journal_abbrev']} {cite['journal_volume']}:{cite['page_first']}, {cite['year']})")Title: THE CRYSTAL STRUCTURE OF HUMAN DEOXYHAEMOGLOBIN AT 1.74 ANGSTROMS RESOLUTION
Resolution: [1.74] Angstrom
Primary citation: PMID 6726807 (J.Mol.Biol. 175:159, 1984)# AlphaFold DB's predicted model for the same protein (P68871) -- experimental vs. predicted
r = requests.get("https://alphafold.ebi.ac.uk/api/prediction/P68871", timeout=30)
r.raise_for_status()
af_entry = r.json()[0]
print("Model created: ", af_entry["modelCreatedDate"])
print("Mean pLDDT confidence:", af_entry["globalMetricValue"])
print("Model version: ", af_entry["latestVersion"])
print("PDB file: ", af_entry["pdbUrl"])Model created: 2025-08-01T00:00:00Z
Mean pLDDT confidence: 97.19
Model version: 6
PDB file: https://alphafold.ebi.ac.uk/files/AF-P68871-F1-model_v6.pdbInteractions and complexes: IntAct and Complex Portal
A real physical interaction for HBB (IntAct), and the real curated tetramer it’s part of (Complex Portal) – the same real biological assembly PDB entry 2HHB is a structure of.
r = requests.get(
"https://www.ebi.ac.uk/intact/ws/interaction/findInteractions/P68871",
timeout=30,
)
r.raise_for_status()
interactions = r.json()["content"]
first = interactions[0]
print(f"{first['moleculeA']} <-> {first['moleculeB']} (IntAct AC: {first['ac']})")
print(f"{len(interactions)} interaction records found for P68871 in IntAct")HBB <-> HBZ (IntAct AC: EBI-22144391)
20 interaction records found for P68871 in IntActr = requests.get(
"https://www.ebi.ac.uk/intact/complex-ws/search/HBB?facets=species_f",
timeout=30,
)
r.raise_for_status()
complexes = r.json()["elements"]
hba_complex = next(c for c in complexes if c["complexAC"] == "CPX-2158")
print(hba_complex["complexName"], "-", hba_complex["complexAC"])
for member in hba_complex["interactors"]:
print(" ", member["stochiometry"], "x", member["name"], f"({member['identifier']})")Hemoglobin HbA complex - CPX-2158
minValue: 4, maxValue: 4 x heme (CHEBI:30413)
minValue: 2, maxValue: 2 x HBA1 (P69905)
minValue: 2, maxValue: 2 x HBB (P68871)Networks: STRING
STRING’s functional-association network, queried for this page’s other running example, CDC28 (yeast).
r = requests.get(
"https://string-db.org/api/json/get_string_ids",
params={"identifiers": "CDC28", "species": 4932},
timeout=30,
)
r.raise_for_status()
string_hits = r.json()
hit = string_hits[0]
print("STRING ID: ", hit["stringId"])
print("Preferred name:", hit["preferredName"])
print("Annotation: ", hit["annotation"][:120], "...")STRING ID: 4932.YBR160W
Preferred name: CDC28
Annotation: Cyclin-dependent kinase (CDK) catalytic subunit; master regulator of mitotic and meiotic cell cycles; alternately associ ...Swiss-Prot vs. TrEMBL, by the numbers
Reviewed (Swiss-Prot) entries are a tiny, curated fraction of UniProtKB; unreviewed (TrEMBL) entries make up almost all of it. Fetched live, so these numbers will keep growing every time this cell is re-run.
def uniprotkb_count(query):
r = requests.get(
"https://rest.uniprot.org/uniprotkb/search",
params={"query": query, "size": 0},
timeout=10,
)
r.raise_for_status()
return int(r.headers["X-Total-Results"])
swiss_prot = uniprotkb_count("reviewed:true")
trembl = uniprotkb_count("reviewed:false")
print(f"Swiss-Prot (reviewed): {swiss_prot:>12,}")
print(f"TrEMBL (unreviewed): {trembl:>12,}")
print(f"TrEMBL is {trembl / swiss_prot:.0f}x larger than Swiss-Prot")Swiss-Prot (reviewed): 575,748
TrEMBL (unreviewed): 149,430,635
TrEMBL is 260x larger than Swiss-ProtAccession prefixes predict the data level
RefSeq’s prefixes are a real, documented convention (see the Session 3 page’s Further Reading) – a tiny lookup table is enough to classify any RefSeq accession without fetching anything.
REFSEQ_PREFIXES = {
"NC": "genomic (complete chromosome)",
"NM": "mRNA",
"NP": "protein",
"NR": "non-coding RNA",
"XM": "predicted mRNA",
"XP": "predicted protein",
}
def refseq_level(accession):
prefix = accession.split("_")[0]
return REFSEQ_PREFIXES.get(prefix, "unknown prefix")
for acc in ["NC_000011.10", "NM_000518.5", "NP_000509.1", "XP_011536997.1"]:
print(f"{acc:16s} -> {refseq_level(acc)}")NC_000011.10 -> genomic (complete chromosome)
NM_000518.5 -> mRNA
NP_000509.1 -> protein
XP_011536997.1 -> predicted protein# A predicted record says so itself: the COMMENT of an XP_ record
handle = Entrez.efetch(db="protein", id="XP_011536997.1", rettype="gb", retmode="text")
record = SeqIO.read(handle, "genbank")
handle.close()
print(record.description)
print(record.annotations["comment"].split("\n")[0][:110])SWI/SNF-related matrix-associated actin-dependent regulator of chromatin subfamily D member 1 isoform X3 [Homo sapiens]
MODEL REFSEQ: This record is predicted by automated computationalBuilding a Boolean PubMed query, verified against the real index
The Session 3 page claims hemoglobin OR malaria behaves like set union, and hemoglobin NOT malaria like set difference. Check that directly against PubMed’s real, live search counts rather than taking it on faith.
import time
def pubmed_count(term):
# NCBI allows about three requests per second without an API key: pause, and retry if refused
for attempt in range(5):
time.sleep(0.5)
r = requests.get(
"https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi",
params={"db": "pubmed", "term": term, "retmode": "json"},
timeout=10,
)
if r.status_code != 429:
break
time.sleep(2 * (attempt + 1))
r.raise_for_status()
return int(r.json()["esearchresult"]["count"])
a = pubmed_count("hemoglobin")
b = pubmed_count("malaria")
a_and_b = pubmed_count("hemoglobin AND malaria")
a_or_b = pubmed_count("hemoglobin OR malaria")
a_not_b = pubmed_count("hemoglobin NOT malaria")
print(f"hemoglobin {a:>10,}")
print(f"malaria {b:>10,}")
print(f"hemoglobin AND malaria {a_and_b:>10,}")
print(f"hemoglobin OR malaria {a_or_b:>10,} (predicted: {a + b - a_and_b:,})")
print(f"hemoglobin NOT malaria {a_not_b:>10,} (predicted: {a - a_and_b:,})")hemoglobin 288,577
malaria 125,442
hemoglobin AND malaria 4,458
hemoglobin OR malaria 409,561 (predicted: 409,561)
hemoglobin NOT malaria 284,119 (predicted: 284,119)Both predictions should match PubMed’s real counts exactly – OR is set union (|A| + |B| - |A ∩ B|) and NOT is set difference (|A| - |A ∩ B|), not just a loose figure of speech.
This is also the kind of question this session’s live quiz and discussion problems build on: given two real search-result counts and their overlap, predict the third before running the query.