Project 3 - Scoring Protein Structures

import pyrosetta
pyrosetta.init("-mute all")
print("PyRosetta ready!")

from pyrosetta.toolbox import pose_from_rcsb

# Download one structure and score it
scorefxn = pyrosetta.get_fa_scorefxn()

pdb_id = "1LYZ"
pose = pose_from_rcsb(pdb_id)
score = scorefxn(pose)

print(f"PDB ID: {pdb_id}")
print(f"Number of residues: {pose.total_residue()}")
print(f"Rosetta score: {score:.2f} REU")

import pandas as pd

# The five Lysozyme structures we will compare
pdb_ids = ["1LYZ", "2LYZ", "4LYZ", "1LSE", "2LZT"]

results = []

for pdb_id in pdb_ids:
    try:
        print(f"Processing {pdb_id}...", end=" ")
        pose = pose_from_rcsb(pdb_id)
        score = scorefxn(pose)
        n_res = pose.total_residue()

        results.append({
            "PDB ID": pdb_id,
            "Residues": n_res,
            "Rosetta Score (REU)": round(score, 2),
            "Score per Residue": round(score / n_res, 2)
        })
        print(f"done — {score:.2f} REU")

    except Exception as e:
        print(f"failed — {e}")
        results.append({
            "PDB ID": pdb_id,
            "Residues": "N/A",
            "Rosetta Score (REU)": "N/A",
            "Score per Residue": "N/A"
        })

df = pd.DataFrame(results)
print("\n" + df.to_string(index=False))

# Filter out any failed entries and sort
df_clean = df[df["Score per Residue"] != "N/A"].copy()
df_clean["Score per Residue"] = df_clean["Score per Residue"].astype(float)
df_clean = df_clean.sort_values("Score per Residue", ascending=True)
df_clean["Rank"] = range(1, len(df_clean) + 1)

print("=" * 60)
print("STRUCTURE RANKING (most negative per-residue score = best)")
print("=" * 60)
print(df_clean[["Rank", "PDB ID", "Residues", "Score per Residue"]].to_string(index=False))
print("=" * 60)

best = df_clean.iloc[0]
print(f"\nBest structure: {best['PDB ID']} ({best['Score per Residue']:.2f} REU/residue)")

import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np

df_plot = df_clean.sort_values("Score per Residue", ascending=True)

# Colour gradient: darker = better score
colors = cm.YlOrRd_r(np.linspace(0.2, 0.8, len(df_plot)))

fig, ax = plt.subplots(figsize=(10, 5))
bars = ax.barh(df_plot["PDB ID"], df_plot["Score per Residue"],
               color=colors, height=0.5)

ax.set_xlabel("Score per Residue (REU)", fontsize=12)
ax.set_title("Lysozyme Structure Quality Comparison\n(More negative = better Rosetta score)",
             fontsize=13, fontweight="bold")
ax.axvline(x=0, color="black", linestyle="--", linewidth=0.8)

# Annotate each bar with the score value
for bar, val in zip(bars, df_plot["Score per Residue"]):
    ax.text(val - 0.1, bar.get_y() + bar.get_height() / 2,
            f"{val:.2f}", va="center", ha="right",
            fontsize=10, fontweight="bold", color="white")

plt.tight_layout()
plt.savefig("structure_comparison.png", dpi=150, bbox_inches="tight")
plt.show()
print("Chart saved as structure_comparison.png")

df_plot = df_clean.sort_values("Score per Residue")

best  = df_plot.iloc[0]
worst = df_plot.iloc[-1]
spread = abs(best["Score per Residue"]) - abs(worst["Score per Residue"])

print("INTERPRETATION")
print("=" * 55)
print(f"\nBest structure:  {best['PDB ID']} ({best['Score per Residue']:.2f} REU/residue)")
print(f"Worst structure: {worst['PDB ID']} ({worst['Score per Residue']:.2f} REU/residue)")
print(f"Spread:          {spread:.2f} REU/residue")
print()
print("What the spread tells us:")
print("  • A larger spread means some structures are significantly")
print("    better modelled than others.")
print("  • If the spread is small (<2 REU/residue), all structures")
print("    are roughly equivalent quality.")
print("  • If the spread is large (>5 REU/residue), you should")
print("    prefer the top-ranked structure for drug design work.")
print()
print(f"Recommendation: Use {best['PDB ID']} as your starting structure.")