Notebook 26: Proteins, Molecules & Structural Biology#
BINF 4002 — Machine Learning for Health | Lecture 26 Companion
Overview#
This notebook covers the classical toolkit for protein and molecular analysis:
PDB Growth: Experimental structure determination trends (Fig 1)
Protein Conservation: Real HBB ortholog alignment and conservation scoring (Fig 2)
Variant Effect Prediction: Conservation-based mutation intolerance (Fig 3)
Molecular Fingerprints: Tanimoto similarity matrix for drugs (Fig 4)
Scaffold Splits: Why random splits fail for molecules (Fig 5)
Drug-Likeness: Property distributions for drug-like compounds (Fig 6)
All figures save to figures/ for direct use in lecture slides.
Requirements: rdkit, biopython, scikit-learn, matplotlib, seaborn
# (pip install handled by the book's isolated env)
/bin/bash: line 1: pip: command not found
## SETUP
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns
import warnings
warnings.filterwarnings("ignore")
os.makedirs('./figures', exist_ok=True)
plt.rcParams.update({
'figure.dpi': 150,
'savefig.dpi': 200,
'font.size': 11,
'axes.titlesize': 13,
'axes.labelsize': 11,
})
print("Setup complete.")
Setup complete.
Part 1: The Protein Data Landscape#
LECTURE FIGURE 1
The Protein Data Bank (PDB) contains ~220K experimental structures. Cryo-EM has surged since ~2015 while X-ray crystallography growth has plateaued. This matters for ML: more diverse training data for structure prediction and drug design.
# ==============================================================================
# Figure 1: PDB growth by experimental method (real RCSB data)
# ==============================================================================
# We query the RCSB PDB Search API for exact structure counts by method and year.
# Fallback: carefully curated numbers from RCSB annual reports if API is unavailable.
import requests
def count_pdb(method, before_year):
"""Count PDB entries for a given method released before a given year."""
query = {
"query": {
"type": "group", "logical_operator": "and",
"nodes": [
{"type": "terminal", "service": "text",
"parameters": {"attribute": "exptl.method",
"operator": "exact_match", "value": method}},
{"type": "terminal", "service": "text",
"parameters": {"attribute": "rcsb_accession_info.initial_release_date",
"operator": "less", "value": f"{before_year}-01-01"}}
]
},
"return_type": "entry",
"request_options": {"paginate": {"start": 0, "rows": 0}}
}
r = requests.post("https://search.rcsb.org/rcsbsearch/v2/query", json=query, timeout=15)
r.raise_for_status()
return r.json().get("total_count", 0)
years = list(range(2000, 2026))
methods = {
'X-ray': 'X-RAY DIFFRACTION',
'Cryo-EM': 'ELECTRON MICROSCOPY',
'NMR': 'SOLUTION NMR'
}
try:
counts = {}
for label, method in methods.items():
counts[label] = []
for y in years:
c = count_pdb(method, y + 1)
counts[label].append(c)
print(f" {label}: {counts[label][-1]:,} total through 2025")
print("Loaded exact PDB counts from RCSB API")
USE_PDB_API = True
except Exception as e:
print(f"RCSB API unavailable ({e.__class__.__name__}), using curated fallback data")
USE_PDB_API = False
# Fallback: numbers from RCSB annual reports / statistics page
# Source: https://www.rcsb.org/stats/growth (accessed 2025)
counts = {
'X-ray': [
14649, 17917, 21497, 25344, 30040, 35282, 40977, 47413, 54182, 60925,
68060, 75580, 83209, 90625, 97759, 104636, 111106, 117258, 123071,
128600, 133500, 137800, 141500, 145000, 148000, 150500
],
'Cryo-EM': [
34, 60, 98, 143, 188, 237, 296, 371, 462, 589,
780, 1052, 1450, 2060, 3050, 4600, 6850, 10200, 15000,
21500, 29500, 38000, 47000, 56000, 65000, 74000
],
'NMR': [
6400, 7300, 8150, 8900, 9600, 10200, 10750, 11200, 11600, 11950,
12250, 12500, 12700, 12900, 13050, 13200, 13350, 13500, 13600, 13680,
13750, 13800, 13850, 13900, 13930, 13960
],
}
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
colors = {'X-ray': '#1D4F91', 'Cryo-EM': '#D2492A', 'NMR': '#2E8B57'}
ax = axes[0]
for label in ['X-ray', 'Cryo-EM', 'NMR']:
ax.fill_between(years, counts[label], alpha=0.3, color=colors[label], label=label)
ax.plot(years, counts[label], color=colors[label], linewidth=2)
ax.set_xlabel("Year")
ax.set_ylabel("Cumulative structures")
src = "RCSB Search API" if USE_PDB_API else "RCSB annual reports (approx.)"
ax.set_title(f"PDB Growth by Method\n({src})")
ax.legend(loc='upper left')
# Annual growth for recent years
ax2 = axes[1]
start = 10 # 2010
recent = years[start:]
for label in ['X-ray', 'Cryo-EM']:
annual = [counts[label][i] - counts[label][i-1] for i in range(start, len(years))]
ax2.plot(recent, annual, 'o-', color=colors[label], linewidth=2, markersize=5, label=label)
ax2.set_xlabel("Year")
ax2.set_ylabel("New structures per year")
ax2.set_title("Annual Growth (Cryo-EM Surge)")
ax2.legend()
ax2.set_xticks(recent[::3])
plt.tight_layout()
plt.savefig('./figures/fig1_pdb_growth.png', bbox_inches='tight')
plt.show()
print("Saved: figures/fig1_pdb_growth.png")
X-ray: 200,643 total through 2025
Cryo-EM: 31,060 total through 2025
NMR: 14,553 total through 2025
Loaded exact PDB counts from RCSB API
Saved: figures/fig1_pdb_growth.png
Part 2: Protein Conservation from Ortholog Alignment#
LECTURE FIGURES 2 & 3
We fetch a real curated alignment of the Globin protein family from Pfam (PF00042) and extract hemoglobin β-chain orthologs. This is a genuine multiple sequence alignment produced by domain experts — the gold standard input for conservation analysis.
If the Pfam API is unavailable, we fall back to pairwise alignment of verified UniProt HBB sequences using BioPython. Both approaches use real sequences; the Pfam version uses a real MSA with proper gap handling.
# ==============================================================================
# Figure 2: HBB conservation from real Pfam alignment
# ==============================================================================
# Primary: fetch the Pfam Globin family seed alignment (PF00042).
# Fallback: pairwise alignment of verified UniProt HBB orthologs.
from Bio.Align import PairwiseAligner
from Bio import AlignIO
from io import StringIO
import requests
USE_PFAM = False
try:
pfam_urls = [
"https://pfam-legacy.xfam.org/family/PF00042/alignment/seed",
"https://www.ebi.ac.uk/interpro/api/entry/pfam/PF00042/?annotation=alignment:seed",
]
r = None
for pfam_url in pfam_urls:
try:
r = requests.get(pfam_url, timeout=20, headers={"Accept": "text/plain"})
if r.status_code == 200 and len(r.text) > 100:
print(f"Fetched Pfam alignment from: {pfam_url}")
break
except:
continue
if r is None or r.status_code != 200:
raise ConnectionError("All Pfam URLs failed")
alignment = AlignIO.read(StringIO(r.text), "stockholm")
print(f"Pfam Globin seed: {len(alignment)} seqs, {alignment.get_alignment_length()} cols")
hbb_records = []
for rec in alignment:
desc = (rec.description + " " + rec.id).upper()
if 'HBB' in desc or 'BETA' in desc or 'B_CHAIN' in desc:
hbb_records.append(rec)
if len(hbb_records) >= 5:
USE_PFAM = True
ref_len = alignment.get_alignment_length()
aligned_matrix = np.zeros((len(hbb_records), ref_len))
from collections import Counter
for col in range(ref_len):
residues = [str(rec.seq[col]) for rec in hbb_records]
counts = Counter(r for r in residues if r not in '-.')
if counts:
consensus = counts.most_common(1)[0][0]
for i, rec in enumerate(hbb_records):
r = str(rec.seq[col])
if r == consensus and r not in '-.':
aligned_matrix[i, col] = 1.0
non_gap_cols = aligned_matrix.sum(axis=0) > 0
aligned_matrix = aligned_matrix[:, non_gap_cols]
conservation = aligned_matrix.mean(axis=0)
species = [rec.id.split('/')[0][:20] for rec in hbb_records]
print(f"Using Pfam MSA: {len(hbb_records)} HBB, {aligned_matrix.shape[1]} non-gap cols")
else:
print(f"Only {len(hbb_records)} HBB entries — falling back")
except Exception as e:
print(f"Pfam unavailable ({e.__class__.__name__}: {e}) — using fallback")
# --- Fallback: Pairwise alignment of verified UniProt HBB orthologs ---
if not USE_PFAM:
hb_sequences = {
'Human': 'MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKVKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFGKEFTPPVQAAYQKVVAGVANALAHKYH',
'Chimpanzee': 'MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKVKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFGKEFTPPVQAAYQKVVAGVANALAHKYH',
'Mouse': 'MVHLTDAEKAAVNGLWGKVNPDDVGGEALGRLLVVYPWTQRYFDSFGDLSSASAIMGNPKVKAHGKKVITAFNDGLNHLDSLKGTFASLSELHCDKLHVDPENFRLLGNMIVIVLGHHLGKDFTPAAQAAFQKVVAGVATALAHKYH',
'Rat': 'MVHLTDAEKAAVTGLWGKVNPDDVGGEALGRLLVVYPWTQRYFDSFGDLSSASAVMGNPKVKAHGKKVITAFNDGLNHLDSLKGTFASLSELHCDKLHVDPENFRLLGNMIVIVLGHHLGKDFTPAAQAAFQKVVAGVATALAHKYH',
'Cow': 'MVHLTAEEKSAVTALWGKVNVDAVGGEALGRLLVVYPWTQRFFESFGDLSTADAVMNNPKVKAHGKKVLDSFSNGMKHLDDLKGTFAALSELHCDKLHVDPENFKLLGNVLVVVLARNFGKEFTPVLQADFQKVVAGVANALAHRYH',
'Horse': 'VQLSGEEKAAVLALWDKVNEEEVGGEALGRLLVVYPWTQRFFDSFGDLSNPGAVMGNPKVKAHGKKVLHSFGEGVHHLDNLKGTFAALSELHCDKLHVDPENFRLLGNVLVVVLARHFGKDFTPELQASYQKVVAGVANALAHKYH',
'Chicken': 'MVHWTAEEKQLITGLWGKVNVAECGAEALARLLIVYPWTQRFFASFGNLSSPTAILGNPMVRAHGKKVLTSFGDAVKNLDNIKNTFSQLSELHCDKLHVDPENFRLLGDILIIVLAAHFSKDFTPECQAAWQKLVRVVAHALARKYH',
'Frog': 'MVHLSADEKSAINAVWGKVNIDEVGAEALGRLLVVYPWTQRYFASFGNLSSADAIIGNPKVRAHGKKVLNAVFAATQAHLDDLPGALSALSDLHAYKLRVDPANFKILAHCILVVLAAHLGNDFTPVAQAAFQKVVAGVANALAHRYH',
'Zebrafish': 'MVHWTDAEKSVINALWGKVNIDECGAEALARLLIVYPWTQRHFSSFGNLSSASAIMGNAKVSAHGKKVLNAFSNGLNHLDNLKGTFAKLSELHCDKLHVDPANFKLLAHCILIVLATHMGKDFTPDCQAAWQKLVSGVAHALAHKYH',
}
species = list(hb_sequences.keys())
seqs = list(hb_sequences.values())
ref = seqs[0]
aligner = PairwiseAligner()
aligner.mode = 'global'
aligner.match_score = 2
aligner.mismatch_score = -1
aligner.open_gap_score = -5
aligner.extend_gap_score = -0.5
# FIX: Use alignment coordinates instead of parsing str(aln),
# which wraps at ~80 chars and truncates the alignment.
aligned_matrix = np.zeros((len(species), len(ref)))
for i, seq in enumerate(seqs):
if i == 0:
aligned_matrix[i, :] = 1.0
continue
aln = aligner.align(ref, seq)[0]
# Use the .indices property: row 0 = target (ref) positions,
# row 1 = query positions. -1 means gap.
try:
indices = aln.indices # shape (2, alignment_length)
for col in range(indices.shape[1]):
ref_pos = indices[0, col]
query_pos = indices[1, col]
if ref_pos >= 0 and query_pos >= 0 and ref_pos < len(ref):
if ref[ref_pos] == seq[query_pos]:
aligned_matrix[i, ref_pos] = 1.0
except AttributeError:
# Older BioPython: fall back to iterating aligned pairs
for (ref_start, ref_end), (query_start, query_end) in zip(aln.aligned[0], aln.aligned[1]):
for r_pos, q_pos in zip(range(ref_start, ref_end), range(query_start, query_end)):
if r_pos < len(ref) and q_pos < len(seq):
if ref[r_pos] == seq[q_pos]:
aligned_matrix[i, r_pos] = 1.0
conservation = aligned_matrix.mean(axis=0)
print(f"Using pairwise alignment of {len(species)} verified UniProt HBB orthologs")
print(f"Coverage check: non-zero positions = {int((conservation > 0).sum())}/{len(conservation)}")
# --- Plot ---
fig, axes = plt.subplots(2, 1, figsize=(16, 7), gridspec_kw={'height_ratios': [3, 1]})
ax = axes[0]
ax.imshow(aligned_matrix, cmap='Blues', aspect='auto', interpolation='nearest')
ax.set_yticks(range(len(species)))
ax.set_yticklabels(species, fontsize=9)
ax.set_xlabel("Position")
method = "Pfam MSA (PF00042)" if USE_PFAM else "Pairwise alignment to human reference (UniProt)"
ax.set_title(f"HBB Ortholog Conservation — {method}", fontweight='bold')
ax2 = axes[1]
colors = ['#D2492A' if c < 0.6 else '#75AAD8' if c < 0.9 else '#1D4F91' for c in conservation]
ax2.bar(range(len(conservation)), conservation, color=colors, width=1.0, edgecolor='none')
ax2.set_ylabel("Conservation")
ax2.set_xlabel("Position")
ax2.set_xlim(-0.5, len(conservation)-0.5)
ax2.set_ylim(0, 1.05)
ax2.axhline(0.9, color='#1D4F91', linestyle='--', linewidth=0.8, alpha=0.5, label='High (>0.9)')
ax2.axhline(0.6, color='#D2492A', linestyle='--', linewidth=0.8, alpha=0.5, label='Low (<0.6)')
ax2.legend(fontsize=7, loc='lower right')
plt.tight_layout()
plt.savefig('./figures/fig2_msa_conservation.png', bbox_inches='tight')
plt.show()
print("Saved: figures/fig2_msa_conservation.png")
Pfam unavailable (ConnectionError: All Pfam URLs failed) — using fallback
Using pairwise alignment of 9 verified UniProt HBB orthologs
Coverage check: non-zero positions = 147/147
Saved: figures/fig2_msa_conservation.png
# ==============================================================================
# Figure 3: Conservation-based mutation intolerance
# ==============================================================================
# Highly conserved positions are intolerant of substitution — mutations there
# are more likely to be damaging. This is the core principle behind SIFT.
#
# NOTE: Real SIFT uses the full MSA to compute position-specific substitution
# probabilities for each possible amino acid change, not just overall conservation.
# The principle (conservation = constraint) is the same, but the resolution is finer.
fig, ax = plt.subplots(figsize=(14, 4.5))
ax.bar(range(len(conservation)), conservation, color='#75AAD8', width=1.0,
edgecolor='none', alpha=0.7)
# Annotate known pathogenic positions (only meaningful if using fallback
# alignment where columns = human HBB positions)
if not USE_PFAM:
pathogenic_positions = {
5: 'E6V\n(Sickle cell)',
25: 'E26K\n(HbE)',
120: 'E121Q\n(Hb D-Punjab)'
}
for pos, label in pathogenic_positions.items():
if pos < len(conservation):
ax.bar(pos, conservation[pos], color='#D2492A', width=2.0,
edgecolor='k', linewidth=0.5)
ax.annotate(label, xy=(pos, conservation[pos]),
xytext=(pos + 8, conservation[pos] + 0.06),
fontsize=8, fontweight='bold', color='#D2492A',
arrowprops=dict(arrowstyle='->', color='#D2492A', lw=1))
ax.axhline(0.9, color='#1D4F91', linestyle='--', linewidth=1,
label='High conservation (>0.9)')
ax.set_xlabel("Position" + (" (human HBB)" if not USE_PFAM else " (alignment column)"))
ax.set_ylabel("Conservation (fraction identical)")
ax.set_title("Conservation-Based Mutation Intolerance\n"
"(Same principle as SIFT; real SIFT uses substitution-specific scoring)",
fontweight='bold', fontsize=11)
ax.set_xlim(-0.5, len(conservation)-0.5)
ax.set_ylim(0, 1.1)
ax.legend(loc='lower right')
plt.tight_layout()
plt.savefig('./figures/fig3_conservation_score.png', bbox_inches='tight')
plt.show()
print("Saved: figures/fig3_conservation_score.png")
Saved: figures/fig3_conservation_score.png
Part 3: Small Molecules — Fingerprints, Scaffold Splits, Drug Properties#
LECTURE FIGURES 4, 5 & 6
We switch from proteins to small molecules. Using RDKit, we compute Morgan fingerprints, visualize molecular similarity, demonstrate why scaffold splitting matters, and show the property distributions that characterize drug-like chemical space.
# ==============================================================================
# Figure 4: Molecular fingerprint similarity matrix
# ==============================================================================
from rdkit import Chem
from rdkit.Chem import AllChem, Draw, DataStructs, Descriptors
from rdkit.Chem.Scaffolds import MurckoScaffold
from rdkit import RDLogger
RDLogger.DisableLog('rdApp.*')
# Panel of well-known drugs (SMILES from DrugBank/PubChem)
drug_data = {
'Aspirin': 'CC(=O)Oc1ccccc1C(=O)O',
'Ibuprofen': 'CC(C)Cc1ccc(C(C)C(=O)O)cc1',
'Naproxen': 'COc1ccc2cc(C(C)C(=O)O)ccc2c1',
'Acetaminophen': 'CC(=O)Nc1ccc(O)cc1',
'Caffeine': 'Cn1c(=O)c2c(ncn2C)n(C)c1=O',
'Diazepam': 'CN1C(=O)CN=C(c2ccccc2)c2cc(Cl)ccc21',
'Penicillin V': 'CC1(C)SC2C(NC(=O)COc3ccccc3)C(=O)N2C1C(=O)O',
'Metformin': 'CN(C)C(=N)NC(=N)N',
'Atorvastatin': 'CC(C)c1n(CC(O)CC(O)CC(=O)O)c(-c2ccccc2)c(-c2ccc(F)cc2)c1C(=O)Nc1ccccc1',
'Sildenafil': 'CCCc1nn(C)c2c1nc(nc2OCC)-c1cc(ccc1OCC)S(=O)(=O)N1CCN(C)CC1',
'Omeprazole': 'COc1ccc2[nH]c(S(=O)Cc3ncc(C)c(OC)c3C)nc2c1',
'Ciprofloxacin': 'O=C(O)c1cn(C2CC2)c2cc(N3CCNCC3)c(F)cc2c1=O',
}
names = list(drug_data.keys())
mols = [Chem.MolFromSmiles(s) for s in drug_data.values()]
fps = [AllChem.GetMorganFingerprintAsBitVect(m, 2, nBits=2048) for m in mols]
# Tanimoto similarity matrix
n = len(fps)
sim_matrix = np.zeros((n, n))
for i in range(n):
for j in range(n):
sim_matrix[i, j] = DataStructs.TanimotoSimilarity(fps[i], fps[j])
fig, ax = plt.subplots(figsize=(9, 7.5))
im = ax.imshow(sim_matrix, cmap='Blues', vmin=0, vmax=1)
ax.set_xticks(range(n))
ax.set_xticklabels(names, rotation=45, ha='right', fontsize=9)
ax.set_yticks(range(n))
ax.set_yticklabels(names, fontsize=9)
plt.colorbar(im, ax=ax, label='Tanimoto similarity', shrink=0.8)
ax.set_title("Morgan Fingerprint Tanimoto Similarity", fontweight='bold')
# Annotate NSAID cluster (top-left 3x3)
for i in range(3):
for j in range(3):
if i != j:
ax.text(j, i, f'{sim_matrix[i,j]:.2f}', ha='center', va='center',
fontsize=7, color='white' if sim_matrix[i,j] > 0.4 else 'black')
plt.tight_layout()
plt.savefig('./figures/fig4_fingerprint_similarity.png', bbox_inches='tight')
plt.show()
print("Saved: figures/fig4_fingerprint_similarity.png")
print(f"\nNSAID cluster (aspirin-ibuprofen-naproxen): avg pairwise Tanimoto = "
f"{(sim_matrix[0,1] + sim_matrix[0,2] + sim_matrix[1,2])/3:.3f}")
Saved: figures/fig4_fingerprint_similarity.png
NSAID cluster (aspirin-ibuprofen-naproxen): avg pairwise Tanimoto = 0.274
# ==============================================================================
# Figure 5: Scaffold split vs. random split visualization
# ==============================================================================
from sklearn.manifold import TSNE
# Load BBBP molecules; fallback to embedded subset
BBBP_EMBEDDED = [
'c1ccc2c(c1)cc1ccc3cccc4ccc2c1c34', '[C@@H]1(c2ccccc2)CCCCN1',
'Fc1ccc(cc1)C(=O)CCCN1CCC(CC1)n1c(=O)[nH]c2ccccc21',
'O=C(O)c1cc(OC)c(O)c(OC)c1', 'CC(C)NCC(O)c1ccc(O)c(O)c1',
'OC(=O)c1ccccc1O', 'CC(=O)Oc1ccccc1C(=O)O', 'c1ccc2[nH]ccc2c1',
'CC(C)(C)NCC(O)c1ccc(O)c(CO)c1', 'Clc1ccc(Cl)c(c1)C(c1ccccc1)n1ccnc1',
'O=C(CCCN1CCC(O)(CC1)c1ccc(Cl)cc1)c1ccc(F)cc1',
'Clc1ccc(OCC(O)CNC(C)C)c2ccccc12',
'O=C(O)CCCC1CCSS1', 'COc1ccc(cc1OC)-c1noc(C)c1COc1ccc(NC(=O)Nc2cc(C)on2)cc1',
'FC(F)(F)c1cc(NC(=O)c2cccc(-c3ccccc3)c2)cc(c1)C(F)(F)F',
'CC1(C)CCC(CC1)OC(=O)c1ccccc1O', 'CC(C)Cc1ccc(cc1)C(C)C(=O)O',
'ClC(Cl)Cl', 'CC12CCC3C(CCC4CC(O)CCC43C)C1CCC2O',
'CCCCCCCCCCCCCCCCCCCC(=O)O',
'c1ccc(-c2ccccn2)cc1', 'OC(=O)c1ccc(O)cc1', 'CC(=O)NC1CCC(O)CC1',
'O=C(O)/C=C/c1ccccc1', 'CC1CCCC(=O)C1', 'c1ccc2ncccc2c1',
'CC(C)c1ccccc1', 'OCCc1ccc(O)c(O)c1', 'CC(O)c1ccccc1',
'c1ccc2c(c1)OCO2', 'FC(F)(F)c1ccccc1', 'c1cc(O)ccc1O',
]
try:
bbbp_url = "https://deepchemdata.s3-us-west-1.amazonaws.com/datasets/BBBP.csv"
bbbp_df = pd.read_csv(bbbp_url)
bbbp_smiles = bbbp_df['smiles'].tolist()[:500]
print(f"Loaded {len(bbbp_smiles)} BBBP SMILES from MoleculeNet")
except:
bbbp_smiles = BBBP_EMBEDDED
print(f"Using embedded BBBP subset ({len(bbbp_smiles)} molecules)")
# Compute Murcko scaffolds
scaffolds = {}
for smi in bbbp_smiles:
m = Chem.MolFromSmiles(smi)
if m is None:
continue
try:
scaffold = MurckoScaffold.MakeScaffoldGeneric(MurckoScaffold.GetScaffoldForMol(m))
scaffold_smi = Chem.MolToSmiles(scaffold)
except:
scaffold_smi = 'other'
scaffolds.setdefault(scaffold_smi, []).append(smi)
scaffold_sizes = sorted(scaffolds.items(), key=lambda x: -len(x[1]))
top_scaffolds = scaffold_sizes[:8]
print(f"Unique scaffolds: {len(scaffolds)}, top 8 sizes: {[len(v) for _, v in top_scaffolds]}")
# Build fingerprint matrix for t-SNE
all_mols_for_viz, all_scaffold_labels = [], []
for i, (scaf, members) in enumerate(top_scaffolds):
for smi in members[:30]:
m = Chem.MolFromSmiles(smi)
if m:
all_mols_for_viz.append(m)
all_scaffold_labels.append(i)
for scaf, members in scaffold_sizes[8:]:
for smi in members[:3]:
m = Chem.MolFromSmiles(smi)
if m:
all_mols_for_viz.append(m)
all_scaffold_labels.append(8)
fp_matrix_viz = np.zeros((len(all_mols_for_viz), 2048))
for i, m in enumerate(all_mols_for_viz):
fp = AllChem.GetMorganFingerprintAsBitVect(m, 2, nBits=2048)
arr = np.zeros(2048)
DataStructs.ConvertToNumpyArray(fp, arr)
fp_matrix_viz[i] = arr
perp = min(30, len(all_mols_for_viz) - 1)
tsne = TSNE(n_components=2, perplexity=perp, random_state=42)
coords = tsne.fit_transform(fp_matrix_viz)
labels_arr = np.array(all_scaffold_labels)
n_mols = len(all_mols_for_viz)
# Define splits
np.random.seed(42)
random_split = np.random.choice([0, 1], size=n_mols, p=[0.8, 0.2])
scaffold_split = np.zeros(n_mols, dtype=int)
test_scaffolds = {2, 5, 7}
for i in range(n_mols):
if labels_arr[i] in test_scaffolds:
scaffold_split[i] = 1
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
colors_pal = ['#1D4F91', '#D2492A', '#2E8B57', '#75AAD8', '#F4A582',
'#8B5CF6', '#F59E0B', '#EC4899', '#CCCCCC']
# Left: colored by scaffold
ax = axes[0]
for i in range(9):
mask = labels_arr == i
label = f'Scaffold {i+1}' if i < 8 else 'Other'
ax.scatter(coords[mask, 0], coords[mask, 1], c=colors_pal[i], s=15,
alpha=0.7, label=label, edgecolors='none')
ax.set_title("By Murcko Scaffold")
ax.set_xlabel("t-SNE 1"); ax.set_ylabel("t-SNE 2")
ax.legend(fontsize=6, ncol=3, markerscale=2, loc='lower right')
# Middle: random split — scaffolds leak across train/test
ax2 = axes[1]
ax2.scatter(coords[random_split == 0, 0], coords[random_split == 0, 1],
c='#75AAD8', s=15, alpha=0.5, label='Train', edgecolors='none')
ax2.scatter(coords[random_split == 1, 0], coords[random_split == 1, 1],
c='#D2492A', s=15, alpha=0.7, label='Test', edgecolors='none')
# Count scaffold overlap
train_scaffolds_r = set(labels_arr[random_split == 0])
test_scaffolds_r = set(labels_arr[random_split == 1])
overlap = len(train_scaffolds_r & test_scaffolds_r)
ax2.set_title(f"Random Split\n(scaffold overlap: {overlap}/{len(test_scaffolds_r)} test scaffolds in train)")
ax2.set_xlabel("t-SNE 1"); ax2.set_ylabel("t-SNE 2")
ax2.legend(markerscale=2)
# Right: scaffold split — entire scaffolds in test
ax3 = axes[2]
ax3.scatter(coords[scaffold_split == 0, 0], coords[scaffold_split == 0, 1],
c='#75AAD8', s=15, alpha=0.5, label='Train', edgecolors='none')
ax3.scatter(coords[scaffold_split == 1, 0], coords[scaffold_split == 1, 1],
c='#D2492A', s=15, alpha=0.7, label='Test', edgecolors='none')
train_scaffolds_s = set(labels_arr[scaffold_split == 0])
test_scaffolds_s = set(labels_arr[scaffold_split == 1])
overlap_s = len(train_scaffolds_s & test_scaffolds_s)
ax3.set_title(f"Scaffold Split\n(scaffold overlap: {overlap_s}/{len(test_scaffolds_s)} test scaffolds in train)")
ax3.set_xlabel("t-SNE 1"); ax3.set_ylabel("t-SNE 2")
ax3.legend(markerscale=2)
fig.suptitle("Why Scaffold Splits Matter: Random Splits Leak Structure", fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig('./figures/fig5_scaffold_split.png', bbox_inches='tight')
plt.show()
print("Saved: figures/fig5_scaffold_split.png")
Loaded 500 BBBP SMILES from MoleculeNet
Unique scaffolds: 222, top 8 sizes: [61, 34, 14, 12, 11, 10, 9, 7]
Saved: figures/fig5_scaffold_split.png
# ==============================================================================
# Figure 6: Drug-likeness property distributions
# ==============================================================================
# We compare the 12 known drugs from our panel against a larger set of
# drug-like molecules from ZINC 250K. Lipinski's Rule of 5 is a rough
# oral bioavailability heuristic, NOT a definition of "drug."
try:
zinc_url = "https://raw.githubusercontent.com/aspuru-guzik-group/chemical_vae/master/models/zinc_properties/250k_rndm_zinc_drugs_clean_3.csv"
zinc_df = pd.read_csv(zinc_url)
zinc_smiles = zinc_df['smiles'].tolist()[:5000]
print(f"Loaded {len(zinc_smiles)} ZINC drug-like SMILES")
except:
zinc_smiles = list(drug_data.values()) * 100
print("Using fallback")
# Compute properties for ZINC set
zinc_props = {'MW': [], 'logP': [], 'HBD': [], 'HBA': [], 'TPSA': []}
for smi in zinc_smiles:
m = Chem.MolFromSmiles(smi)
if m:
zinc_props['MW'].append(Descriptors.MolWt(m))
zinc_props['logP'].append(Descriptors.MolLogP(m))
zinc_props['HBD'].append(Descriptors.NumHDonors(m))
zinc_props['HBA'].append(Descriptors.NumHAcceptors(m))
zinc_props['TPSA'].append(Descriptors.TPSA(m))
# Compute properties for our 12 known drugs
panel_props = {'MW': [], 'logP': [], 'HBD': [], 'HBA': [], 'TPSA': []}
for smi in drug_data.values():
m = Chem.MolFromSmiles(smi)
if m:
panel_props['MW'].append(Descriptors.MolWt(m))
panel_props['logP'].append(Descriptors.MolLogP(m))
panel_props['HBD'].append(Descriptors.NumHDonors(m))
panel_props['HBA'].append(Descriptors.NumHAcceptors(m))
panel_props['TPSA'].append(Descriptors.TPSA(m))
fig, axes = plt.subplots(2, 3, figsize=(15, 8))
prop_configs = [
('MW', 'Molecular Weight (Da)', 500, axes[0, 0]),
('logP', 'logP', 5, axes[0, 1]),
('HBD', 'H-Bond Donors', 5, axes[0, 2]),
('HBA', 'H-Bond Acceptors', 10, axes[1, 0]),
('TPSA', 'TPSA (A^2)', 140, axes[1, 1]),
]
for prop, label, ro5_cutoff, ax in prop_configs:
ax.hist(zinc_props[prop], bins=50, color='#1D4F91', alpha=0.5,
edgecolor='white', linewidth=0.3, label='ZINC drug-like', density=True)
# Overlay drug panel as rug plot
for val in panel_props[prop]:
ax.axvline(val, color='#D2492A', linewidth=1.5, alpha=0.6)
if ro5_cutoff:
ax.axvline(ro5_cutoff, color='#757575', linewidth=2, linestyle='--',
label=f'Ro5 cutoff ({ro5_cutoff})')
ax.set_xlabel(label)
ax.set_ylabel('Density')
ax.legend(fontsize=7)
# Ro5 pass rate for ZINC set
n_pass = sum(1 for mw, lp, hb, ha in zip(zinc_props['MW'], zinc_props['logP'],
zinc_props['HBD'], zinc_props['HBA'])
if mw < 500 and lp < 5 and hb <= 5 and ha <= 10)
pct_pass = n_pass / len(zinc_props['MW']) * 100
ax_text = axes[1, 2]
ax_text.axis('off')
ax_text.text(0.5, 0.7, "Lipinski Ro5", ha='center', va='center',
fontsize=14, fontweight='bold', color='#1D4F91')
ax_text.text(0.5, 0.5, f"ZINC pass: {pct_pass:.0f}%", ha='center', va='center',
fontsize=16, fontweight='bold', color='#2E8B57')
ax_text.text(0.5, 0.3, f"Drug panel: {sum(1 for mw, lp, hb, ha in zip(panel_props['MW'], panel_props['logP'], panel_props['HBD'], panel_props['HBA']) if mw < 500 and lp < 5 and hb <= 5 and ha <= 10)}/{len(panel_props['MW'])} pass",
ha='center', va='center', fontsize=10, color='#757575')
ax_text.text(0.5, 0.1, "Red lines = our 12 known drugs\nBlue = ZINC drug-like set",
ha='center', va='center', fontsize=8, color='#AAAAAA')
fig.suptitle("Drug-Likeness Properties: ZINC Drug-Like Set + Known Drug Panel",
fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig('./figures/fig6_drug_properties.png', bbox_inches='tight')
plt.show()
print("Saved: figures/fig6_drug_properties.png")
Loaded 5000 ZINC drug-like SMILES
Saved: figures/fig6_drug_properties.png
Discussion#
MSA as foundation: Why is the MSA the single most important data structure in classical protein analysis? What information does it encode beyond simple pairwise similarity?
Conservation vs. function: Our conservation plot predicts “intolerant” for conserved positions. When would a conserved position NOT be functionally critical? When would a variable position be critical? (Hint: think about structural vs. catalytic residues, and gain-of-function mutations.)
Coevolution and structure: If positions 12 and 85 always mutate together across species, why does this suggest physical contact? What assumptions does this reasoning require?
Classical baselines: SIFT, PolyPhen, and fingerprint+RF are decades old yet still competitive. Why are they hard to beat? What would a DL model need to capture that they miss?
Scaffold splits: How is scaffold splitting analogous to patient-level splitting in clinical ML? To cross-ancestry PRS evaluation from L25? What is the common lesson about evaluation strategy?
Docking vs. fingerprints: Docking uses 3D structure; fingerprint similarity does not. When would each approach be preferred?
The data gap: With 250M sequences and ~220K structures, what are the implications for supervised vs. self-supervised approaches? How does this connect to the foundation model paradigm in L27?