Project 5 - Alanine Scanning of a Binding Site
Drug Design with PyRosetta — Identifying Critical Residues for Drug Binding
What You Will Learn
- How alanine scanning reveals which residues are essential for drug binding
- How to compute ΔΔG (delta-delta-G) — the energetic cost of mutating a residue
- How to distinguish “hot spot” residues from dispensable ones
- How structure-activity relationships (SAR) connect computational results to drug design decisions
- How to build a prioritized list of residues for experimental validation
| Time: ~45–60 minutes | Prerequisite: Projects 1–4 completed | Difficulty: Intermediate |
Background: What Is Alanine Scanning?
You performed alanine scanning in Project 1 to test protein stability — mutating each residue to Alanine and watching the stability score change. Here we apply the same technique with a critical twist: we are not just measuring whether the protein becomes more or less stable. We are measuring how much each residue contributes to holding the drug in place.
This is called binding site alanine scanning, and it is one of the most powerful experimental and computational tools in drug design.
The logic:
Alanine is the simplest amino acid with a side chain — it has just a methyl group (-CH₃) with no special chemical properties. When you mutate a residue to Alanine, you erase whatever that residue’s side chain was doing (forming hydrogen bonds, van der Waals contacts, electrostatic interactions) while keeping the protein backbone intact.
- If erasing that side chain dramatically weakens drug binding → the residue was critical (a “hot spot”)
- If erasing it barely affects binding → the residue was dispensable
The ΔΔG Concept
We measure binding contribution using ΔΔG (pronounced “delta-delta-G”):
ΔG_binding = binding energy of wild-type protein + drug
ΔG_binding_mut = binding energy of mutant protein + drug
ΔΔG = ΔG_binding_mut − ΔG_binding
ΔΔG > 0 → mutation weakens binding (hot spot residue — critical)
ΔΔG < 0 → mutation strengthens binding (rare, but possible)
ΔΔG ≈ 0 → mutation has little effect (dispensable residue)
A ΔΔG above +1.0 REU is considered significant. Residues above +2.0 REU are considered strong hot spots — validated targets for structure-activity relationship (SAR) studies.
Setup: Using the Binding Site from Project 4
This project uses the CSV file generated in Project 4’s Cell 6. If you have not completed Project 4, you can recreate the contact residues using the code from that project.
Project Workflow
Here is the complete pipeline this project follows. Each numbered step maps to one code cell below.
Cell 1 — Initialize and Load the Contact Residues
What this cell does: Starts PyRosetta and loads the binding site CSV file produced by Project 4 Cell 6. If that file is missing, the except block prints a helpful message instead of crashing.
This cell introduces pd.read_csv() — the counterpart to pd.to_csv() used in Project 4 — and a named exception type, FileNotFoundError.
pd.read_csv("filename.csv") reads a CSV file from disk and reconstructs it as a DataFrame, with the same columns and values that were saved. It is the exact inverse of df.to_csv() used in Project 4.
This is how data flows between projects: Project 4 saved binding_site_report.csv; Project 5 loads it. This pattern — one script saves, another loads — is the foundation of reusable scientific pipelines.
In Projects 3 and 4 we used except Exception as e to catch any error. Here we use except FileNotFoundError — a named exception type that only fires when a file does not exist.
Using a specific exception type is better practice when you know exactly what can go wrong: it catches the one expected failure (missing file) without accidentally hiding other bugs. Python has many named exception types — FileNotFoundError, ValueError, KeyError, ZeroDivisionError — each describing a specific category of problem.
import pyrosetta
from pyrosetta import *
from pyrosetta.toolbox import *
import pandas as pd
pyrosetta.init(extra_options="-ignore_unrecognized_res false -load_PDB_components true -mute all")
print("PyRosetta ready!")
# Load the binding site report from Project 4
try:
df_contacts = pd.read_csv("binding_site_report.csv")
print(f"Loaded {len(df_contacts)} contact residues from Project 4.")
print(df_contacts[["AA (3-letter)", "Position (PDB)", "Min Distance (Å)", "Chemical Class"]].to_string(index=False))
except FileNotFoundError:
print("binding_site_report.csv not found.")
print("Please complete Project 4 Cell 6 first, or manually define contact residues.")
Cell 2 — Load COX-2 and Calculate Wild-Type Binding Energy
What this cell does: Downloads COX-2, cleans the structure (same steps as Project 4), then computes the wild-type binding energy as a baseline. The key new idea here is wrapping the binding energy calculation inside a reusable function using def and return.
We need to call this calculation once for the wild-type protein, then once for every mutant — potentially 20+ times. Writing it as a function means we write the logic once and call it by name each time.
def function_name(parameter1, parameter2): defines a new function. Everything indented underneath it is the function body — the code that runs when the function is called.
return value sends a result back to wherever the function was called from. Without return, a function runs its code but hands nothing back.
In this cell, compute_binding_energy(pose, ligand_res, scorefxn) takes three inputs (a pose, the ligand position, and the scoring function), runs the subtraction calculation, and returns the binding energy as a single number. Every time you call it later — wt_binding = compute_binding_energy(wt_pose, ligand_pos_clean, scorefxn) — Python jumps into the function body, runs all the code, and substitutes the returned value back in place of the call.
| Part of the function definition | What it means |
|---|---|
def compute_binding_energy( | Start defining a function called compute_binding_energy |
pose, ligand_res, scorefxn | The three inputs (parameters) the function expects when called |
): | End of the parameter list; the body follows, indented |
return score_complex - (...) | Calculate and send the result back to the caller |
score_complex.score_protein.score_ligand.score_complex − (score_protein + score_ligand). A negative result means the complex is more stable than the separated parts — the drug binds favourably.import urllib.request
from pyrosetta import pose_from_pdb
# Download and load COX-2 + Ibuprofen
url = "https://files.rcsb.org/download/4PH9.pdb"
urllib.request.urlretrieve(url, "4PH9.pdb")
pose_full = pose_from_pdb("4PH9.pdb")
scorefxn = pyrosetta.get_fa_scorefxn()
# Isolate chain A protein + IBP
ligand_pos_full = None
keep_positions = []
for i in range(1, pose_full.total_residue() + 1):
res = pose_full.residue(i)
chain = pose_full.pdb_info().chain(i)
if res.name() == "pdb_IBP":
ligand_pos_full = i
keep_positions.append(i)
elif chain == "A" and res.is_protein():
keep_positions.append(i)
wt_pose = pose_full.clone()
delete_positions = [i for i in range(1, pose_full.total_residue() + 1)
if i not in keep_positions]
for i in reversed(delete_positions):
wt_pose.delete_residue_slow(i)
# Find IBP in the clean pose
ligand_pos_clean = None
for i in range(1, wt_pose.total_residue() + 1):
if not wt_pose.residue(i).is_protein():
ligand_pos_clean = i
break
print(f"Clean structure: {wt_pose.total_residue()} residues")
print(f"Ibuprofen at clean position: {ligand_pos_clean}")
# Compute wild-type binding energy using the subtraction method (see Project 2)
def compute_binding_energy(pose, ligand_res, scorefxn):
"""
Binding energy = Complex score − (protein-only score + ligand-only score)
"""
score_complex = scorefxn(pose)
# Score the protein alone (delete ligand from a copy)
prot_pose = pose.clone()
prot_pose.delete_residue_slow(ligand_res)
score_protein = scorefxn(prot_pose)
# Score the ligand alone (delete all protein from a copy)
lig_pose = pose.clone()
protein_positions = [i for i in range(1, pose.total_residue() + 1)
if pose.residue(i).is_protein()]
for i in reversed(protein_positions):
lig_pose.delete_residue_slow(i)
score_ligand = scorefxn(lig_pose)
return score_complex - (score_protein + score_ligand)
wt_binding = compute_binding_energy(wt_pose, ligand_pos_clean, scorefxn)
print(f"\nWild-type binding energy: {wt_binding:.2f} REU")
print("This is the baseline. All mutant energies will be compared to this.")
Cell 3 — Mutate Each Binding Site Residue to Alanine
What this cell does: The core of the project — loops over every residue in the binding site, mutates it to Alanine, computes ΔΔG, and classifies the result. This cell runs the function you defined in Cell 2 once per residue, which is why writing it as a reusable function mattered.
Expected runtime: 5–15 minutes. Each iteration scores the mutant three times (complex, protein alone, ligand alone) and there may be 15–25 residues to scan.
for _, row in df.iterrows(): loops through a DataFrame one row at a time. On each iteration:
_ is the row index (we use _ as a name to signal "I don't need this value").
row is a Series object — it behaves like a dictionary keyed by column name. Access values with row["Column Name"].
| (index) | AA (1-letter) | Position (Rosetta) | Position (PDB) | Chemical Class |
|---|---|---|---|---|
| 0 ← iteration 1 | V | 102 | 349 | Hydrophobic |
| 1 ← iteration 2 | L | 117 | 364 | Hydrophobic |
| 2 ← iteration 3 | S | 121 | 368 | Polar |
row["AA (1-letter)"] → "V", row["Position (Rosetta)"] → 102, and so on.mutate_residue(pose, rosetta_position, "A") replaces the side chain of the residue at rosetta_position with an Alanine side chain (a plain methyl group), leaving the protein backbone unchanged.
We always call it on a clone() of the wild-type pose — never on wt_pose directly. This preserves the original structure so the next loop iteration starts from the same unmodified baseline.
continue immediately jumps to the next iteration of a loop, skipping all remaining code in the current one. Here it is used to skip residues that are already Alanine — there is no point in mutating Alanine to Alanine.
Without continue, you would need to wrap the entire mutation logic in an if aa_wt != "A": block, which would indent everything further. continue keeps the code flatter and more readable by handling the skip case first and early.
rosetta_pos, pdb_pos, aa_wt, and chem_class from the current DataFrame row.aa_wt == "A", append a row with ΔΔG = 0.0 and call continue to move to the next residue immediately.mutant_pose = wt_pose.clone() then mutate_residue(mutant_pose, rosetta_pos, "A"). The original wt_pose is unchanged.compute_binding_energy() on the mutant, subtract the wild-type value: ddg = mut_binding - wt_binding.scan_results.scan_results to a DataFrame and sort by ΔΔG descending — largest (most critical) residues appear first.from pyrosetta.toolbox import mutate_residue
import numpy as np
print("Running alanine scan on binding site residues...")
print("(This may take 5–15 minutes depending on your machine)\n")
scan_results = []
for _, row in df_contacts.iterrows():
rosetta_pos = int(row["Position (Rosetta)"])
pdb_pos = int(row["Position (PDB)"])
aa_wt = row["AA (1-letter)"]
chem_class = row["Chemical Class"]
# Skip residues that are already Alanine — no meaningful mutation to make
if aa_wt == "A":
print(f" Skipping {row['AA (3-letter)']} {pdb_pos} — already Alanine")
scan_results.append({
"PDB Position": pdb_pos,
"Wild-Type AA": aa_wt,
"Chemical Class": chem_class,
"ΔΔG (REU)": 0.0,
"Effect": "Already Ala"
})
continue
# Create mutant pose
mutant_pose = wt_pose.clone()
mutate_residue(mutant_pose, rosetta_pos, "A")
# Compute mutant binding energy
try:
mut_binding = compute_binding_energy(mutant_pose, ligand_pos_clean, scorefxn)
ddg = mut_binding - wt_binding # ΔΔG
if ddg >= 2.0:
effect = "Hot spot"
elif ddg >= 1.0:
effect = "Moderate"
elif ddg <= -1.0:
effect = "Stabilizing"
else:
effect = "Dispensable"
print(f" {aa_wt}{pdb_pos}A: ΔΔG = {ddg:+.2f} REU — {effect}")
scan_results.append({
"PDB Position": pdb_pos,
"Wild-Type AA": aa_wt,
"Chemical Class": chem_class,
"ΔΔG (REU)": round(ddg, 2),
"Effect": effect
})
except Exception as e:
print(f" {aa_wt}{pdb_pos}A: ERROR — {e}")
print("\nAlanine scan complete!")
df_scan = pd.DataFrame(scan_results).sort_values("ΔΔG (REU)", ascending=False)
Cell 4 — Display the Ranked Results
What this cell does: Prints the scan results as a neatly aligned table, lists the hot spots, and saves the full results to a CSV for Project 6. The new skill here is f-string column formatting — controlling the width and alignment of each printed value to produce clean, readable columns.
Inside an f-string, you can specify how a value is padded and aligned by adding a format spec after a colon: {value:format_spec}. This project uses three format specs together:
| Format spec | Meaning | Example |
|---|---|---|
:<12 | Left-align in a 12-character-wide field (pad with spaces on the right) | "V102A " |
:>10 | Right-align in a 10-character-wide field (pad with spaces on the left) | " +2.45" |
:>+10.2f | Right-align, always show sign (+/−), 2 decimal places, width 10 | " +2.45" |
By giving every value a fixed width, all rows line up into neat columns — the same technique spreadsheets use. Without fixed widths, short values like "V10A" and long values like "PHE349A" would cause columns to shift unpredictably.
print("\n" + "=" * 65)
print("ALANINE SCAN RESULTS — COX-2 Binding Site (4PH9 + Ibuprofen)")
print("=" * 65)
print(f"{'Mutation':<12} {'ΔΔG (REU)':>10} {'Chemical Class':<20} {'Effect'}")
print("-" * 65)
for _, row in df_scan.iterrows():
mutation = f"{row['Wild-Type AA']}{row['PDB Position']}A"
print(f"{mutation:<12} {row['ΔΔG (REU)']:>+10.2f} "
f"{row['Chemical Class']:<20} {row['Effect']}")
print("=" * 65)
hot_spots = df_scan[df_scan["Effect"] == "Hot spot"]
print(f"\nHot spots identified ({len(hot_spots)} residues with ΔΔG ≥ 2.0 REU):")
for _, row in hot_spots.iterrows():
print(f" → {row['Wild-Type AA']}{row['PDB Position']}A : ΔΔG = {row['ΔΔG (REU)']:+.2f} REU")
df_scan.to_csv("alanine_scan_results.csv", index=False)
print("\nFull results saved to alanine_scan_results.csv")
Cell 5 — Plot the Alanine Scan Results
What this cell does: Draws a vertical bar chart — one bar per mutation, coloured by effect category, with horizontal threshold lines at ΔΔG = 1.0 and 2.0. This cell introduces three new Matplotlib ideas: ax.bar() for vertical bars, ax.axhline() for horizontal reference lines, plt.xticks(rotation=) for angled labels, and .tolist() for extracting a column's values.
In Project 3 we used ax.barh() for horizontal bars (labels on the y-axis, values on the x-axis). Here we use ax.bar() for vertical bars (labels on the x-axis, values on the y-axis).
Vertical bars are better here because we want the ΔΔG threshold lines to be horizontal — easy reference lines crossing the chart left-to-right. With horizontal bars, the threshold lines would have to be vertical, which is visually harder to compare.
ax.axhline(y=2.0, color="red", linestyle="--") draws a horizontal dashed line across the entire chart at y = 2.0. axhline = "axis horizontal line." The counterpart axvline (used in Projects 3 and 4) draws a vertical line.
We add two threshold lines — at y = 1.0 and y = 2.0 — to visually separate dispensable residues from moderate ones and hot spots. Whenever your data has a scientifically meaningful cutoff, a reference line makes the chart self-explanatory without needing to read the axis carefully.
With 15–25 mutation labels on the x-axis (e.g. V349A, L364A, S368A…), horizontal text would overlap and become unreadable. plt.xticks(rotation=45, ha="right") rotates each label 45° clockwise. ha="right" (horizontal alignment) anchors each label's right edge to its tick mark so the angled text sits directly below the correct bar.
df["column"].tolist() extracts all values from a DataFrame column into a plain Python list. Some functions — like the list comprehension used to build colors — work more reliably with a plain list than with a pandas Series. It is a small but common conversion in data visualisation code.
df_scan by ΔΔG descending. Extract mutation labels (built with a list comprehension over iterrows()), values (.tolist()), and colours (mapped from the EFFECT_COLORS dict).ax.bar(labels, values, color=colors) — one bar per mutation, height = ΔΔG, colour = effect category.ax.axhline(y=2.0) and ax.axhline(y=1.0) draw dashed horizontal lines across the chart at the hot spot and moderate cutoffs.plt.xticks(rotation=45, ha="right") angles the mutation labels so they do not overlap.mpatches.Patch entries (seen in Project 4) build the colour legend. plt.savefig() writes to disk.import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
# Assign colors by effect
EFFECT_COLORS = {
"Hot spot": "#c0392b",
"Moderate": "#e67e22",
"Dispensable": "#95a5a6",
"Stabilizing": "#27ae60",
"Already Ala": "#bdc3c7"
}
df_plot = df_scan.sort_values("ΔΔG (REU)", ascending=False)
labels = [f"{row['Wild-Type AA']}{row['PDB Position']}A" for _, row in df_plot.iterrows()]
values = df_plot["ΔΔG (REU)"].tolist()
colors = [EFFECT_COLORS.get(e, "#95a5a6") for e in df_plot["Effect"]]
fig, ax = plt.subplots(figsize=(13, 6))
bars = ax.bar(labels, values, color=colors, edgecolor="white", linewidth=0.5)
# Reference lines
ax.axhline(y=2.0, color="#c0392b", linestyle="--", linewidth=1.0, label="Hot spot threshold (2.0)")
ax.axhline(y=1.0, color="#e67e22", linestyle="--", linewidth=0.8, label="Moderate threshold (1.0)")
ax.axhline(y=0.0, color="black", linestyle="-", linewidth=0.5)
ax.set_xlabel("Mutation", fontsize=11)
ax.set_ylabel("ΔΔG (REU)", fontsize=11)
ax.set_title("Alanine Scan — COX-2 Binding Site\n(Positive ΔΔG = residue contributes to Ibuprofen binding)",
fontsize=13, fontweight="bold")
plt.xticks(rotation=45, ha="right", fontsize=9)
# Legend
patches = [mpatches.Patch(color=v, label=k) for k, v in EFFECT_COLORS.items()]
ax.legend(handles=patches, loc="upper right", fontsize=9)
plt.tight_layout()
plt.savefig("alanine_scan.png", dpi=150, bbox_inches="tight")
plt.show()
print("Chart saved as alanine_scan.png")
Cell 6 — Structure-Activity Relationship (SAR) Summary
What this cell does: Filters the results into three groups (hot spots, moderate, dispensable), prints a count summary, then generates a specific drug design recommendation for each hot spot based on its chemical class. The new idea here is the inline conditional expression used inside an f-string.
Python allows a compact one-line conditional: value_if_true if condition else value_if_false. This is called a ternary expression or inline conditional.
In this cell it is used to print different recommendations depending on the residue's chemical class:
'add hydrophobic group to drug' if chem == 'Hydrophobic' else 'add H-bond donor/acceptor' if chem == 'Polar' else 'add complementary charge to drug'
This chains two conditionals: first check if Hydrophobic, then if Polar, otherwise fall through to the charge recommendation. It is equivalent to a three-branch if/elif/else block, written in one line so it fits neatly inside a print statement.
Use inline conditionals sparingly — they are compact but can be hard to read if they get long. A full if/elif/else block is always clearer for complex logic.
hot_spots = df_scan[df_scan["Effect"] == "Hot spot"]
moderate = df_scan[df_scan["Effect"] == "Moderate"]
dispensable = df_scan[df_scan["Effect"] == "Dispensable"]
print("STRUCTURE-ACTIVITY RELATIONSHIP (SAR) SUMMARY")
print("=" * 60)
print(f"\nTotal binding site residues scanned: {len(df_scan)}")
print(f" Hot spots (ΔΔG ≥ 2.0): {len(hot_spots)} residues")
print(f" Moderate (ΔΔG ≥ 1.0): {len(moderate)} residues")
print(f" Dispensable (ΔΔG < 1.0): {len(dispensable)} residues")
print("\nDrug design implications:")
print("-" * 60)
if len(hot_spots) > 0:
print("\n✓ HOT SPOTS — Design new drugs to maximally engage these:")
for _, row in hot_spots.iterrows():
chem = row["Chemical Class"]
pos = row["PDB Position"]
aa = row["Wild-Type AA"]
print(f" {aa}{pos}: {chem} residue — "
f"{'add hydrophobic group to drug' if chem == 'Hydrophobic' else 'add H-bond donor/acceptor' if chem == 'Polar' else 'add complementary charge to drug'}")
print("\n✗ DISPENSABLE — Drug modifications here will not improve binding:")
for _, row in dispensable.iterrows():
print(f" {row['Wild-Type AA']}{row['PDB Position']}: "
f"low contribution (ΔΔG = {row['ΔΔG (REU)']:+.2f} REU)")
print("\nNext step: Use Project 6 to simulate resistance mutations at hot spot positions.")
alanine_scan_results.csv saved in Cell 4 will be loaded by Project 6. Keep it in the same folder as your Project 6 notebook. Discussion Questions
1. How many hot spot residues did you find in the COX-2 binding site? Are they all the same chemical class, or a mix?
2. If you were designing a new drug to replace Ibuprofen, which residue(s) would you most want your new drug to interact with? Why?
3. One assumption we make is that Alanine-mutant proteins still fold correctly and maintain the same overall shape. When might this assumption break down? (Hint: think about Glycine and Proline.)
4. We used the subtraction method to calculate binding energy. What are the limitations of this approach compared to more sophisticated methods like FEP (Free Energy Perturbation)?
5. Experimental alanine scanning involves actually synthesizing and testing each mutant protein. How might the computational results from this project guide which mutations to prioritize for expensive lab testing?
Key Vocabulary
| Term | Definition |
|---|---|
| Alanine scanning | Systematically mutating each residue to Alanine to measure its contribution |
| ΔΔG | The change in binding energy caused by a mutation. Positive = weaker binding |
| Hot spot | A residue whose side chain contributes strongly to binding (ΔΔG ≥ 2.0 REU) |
| Wild type | The original, naturally occurring protein sequence before any mutations |
| Dispensable residue | A residue that contributes little to binding; mutating it has minimal effect |
| SAR | Structure-Activity Relationship — the study of how molecular structure affects activity |
| Binding affinity | How tightly a drug binds to its target. Higher affinity = tighter binding = better drug |
pd.read_csv() | Loads a saved CSV file back into a DataFrame; the inverse of .to_csv() |
FileNotFoundError | A named Python exception raised when a file does not exist |
def / return | Define a reusable function (def) and send a result back to the caller (return) |
mutate_residue() | PyRosetta function that replaces one residue’s side chain with a different amino acid |
continue | Skips the rest of the current loop iteration and jumps to the next one |
df.iterrows() | Iterates over a DataFrame row by row; each iteration yields the index and a row Series |
| f-string column alignment | Format specs like :<12 and :>+10.2f that pad values to fixed widths for neat columns |
.tolist() | Converts a DataFrame column into a plain Python list |
ax.bar() | Draws vertical bars (vs ax.barh() for horizontal bars) |
ax.axhline() | Draws a horizontal reference line across the chart at a specified y value |
plt.xticks(rotation=) | Rotates x-axis tick labels to prevent overlapping |
| Inline conditional | Compact one-line value_if_true if condition else value_if_false expression |
| Functional group | A specific chemical group on a drug molecule (e.g. -OH, -NH₂) that drives interactions |
| Medicinal chemistry | The science of designing and optimizing drug molecules for maximum efficacy and safety |