How does a small drug fit into a protein pocket?

from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
import py3Dmol
from scipy.spatial.distance import cdist
from IPython.display import display as _display
def display(value):
    if isinstance(value,pd.DataFrame):
        table=value.to_html(max_rows=12,escape=True)
        style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
        return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
    return _display(value)
RECIPE_ID='07-imatinib-pocket'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
                      "axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")

def load_snapshot(recipe_id, expected_sha256, snapshot_key):
    """Read embedded data from this notebook, including a renamed copy."""
    for path in sorted(Path.cwd().glob("*.ipynb")):
        try:
            notebook = json.loads(path.read_text(encoding="utf-8"))
            metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
            blob = metadata.get("snapshots", {}).get(snapshot_key)
            if metadata.get("recipe_id") != recipe_id or not blob:
                continue
            if blob.get("sha256") != expected_sha256:
                continue
            raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
            if hashlib.sha256(raw).hexdigest() != expected_sha256:
                raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
            return json.loads(raw)
        except (OSError, json.JSONDecodeError):
            continue
    raise FileNotFoundError(
        "Save/open the complete .ipynb in the active working folder. "
        "The embedded snapshot could not be located for " + recipe_id)

SNAPSHOT_HASHES={'imatinib': '5c917b32818fda4715c7887577a22d275c25e7a80f9ebf0e6b55e923f2c1f187'}
def snapshot(key):
    return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))

def parse_pdb(text):
    """First model; choose highest occupancy altloc, then blank/A/lexical order."""
    rows, first_model, inside = {}, None, True
    for line in text.splitlines():
        record = line[:6].strip()
        if record == "MODEL":
            number = line[10:14].strip()
            if first_model is None:
                first_model = number
            inside = number == first_model
        if record == "ENDMDL" and inside:
            break
        if record not in ("ATOM","HETATM") or not inside:
            continue
        try:
            atom = dict(record=record,serial=int(line[6:11]),atom=line[12:16].strip(),altloc=line[16:17].strip(),
                residue=line[17:20].strip(),chain=line[21:22].strip(),resnum=int(line[22:26]),icode=line[26:27].strip(),
                x=float(line[30:38]),y=float(line[38:46]),z=float(line[46:54]),
                occupancy=float(line[54:60].strip() or 0),bfactor=float(line[60:66].strip() or 0),element=line[76:78].strip())
        except ValueError as exc:
            raise ValueError("Invalid PDB numeric fields") from exc
        key = (record,atom["chain"],atom["resnum"],atom["icode"],atom["residue"],atom["atom"])
        rank = (atom["occupancy"],atom["altloc"]=="",atom["altloc"]=="A")
        old = rows.get(key)
        if old is None or rank > old[0] or (rank == old[0] and atom["altloc"] < old[1]["altloc"]):
            rows[key] = (rank,atom)
    columns = ["record","serial","atom","altloc","residue","chain","resnum","icode","x","y","z","occupancy","bfactor","element"]
    atoms = pd.DataFrame([entry[1] for entry in rows.values()],columns=columns)
    if atoms.empty:
        raise ValueError("No atoms in the first model")
    if not np.isfinite(atoms[["x","y","z","occupancy"]].to_numpy()).all():
        raise ValueError("Nonfinite coordinates or occupancy")
    return atoms

print("Method ready: parse_pdb")

def residue_minimum_distances(left, right):
    """Minimum Euclidean atom distance per residue pair, retaining chain/insertion IDs.

    Caller explicitly selects positive-occupancy heavy atoms. This measures
    geometric proximity, not a hydrogen bond, affinity, buried area or energy.
    At most one residue's atom-distance block is held in memory at a time.
    """
    keys = ['chain', 'resnum', 'icode', 'residue']
    columns = keys + ['atom_count']
    if left.empty or right.empty:
        raise ValueError('Both atom selections must be nonempty')
    if not np.isfinite(left[['x','y','z']].to_numpy()).all() or not np.isfinite(right[['x','y','z']].to_numpy()).all():
        raise ValueError('Coordinates must be finite')
    left = left.reset_index(drop=True)
    right = right.reset_index(drop=True)
    left_groups = list(left.groupby(keys, sort=False).indices.items())
    right_groups = list(right.groupby(keys, sort=False).indices.items())
    left_residues = pd.DataFrame([list(key)+[len(indices)] for key,indices in left_groups], columns=columns)
    right_residues = pd.DataFrame([list(key)+[len(indices)] for key,indices in right_groups], columns=columns)
    right_codes = np.empty(len(right), dtype=int)
    for j, (_, indices) in enumerate(right_groups):
        right_codes[indices] = j
    matrix = np.full((len(left_groups),len(right_groups)), np.inf)
    right_xyz = right[['x','y','z']].to_numpy(dtype=float)
    for i, (_, indices) in enumerate(left_groups):
        block = cdist(left.iloc[indices][['x','y','z']].to_numpy(dtype=float), right_xyz)
        np.minimum.at(matrix[i], right_codes, block.min(axis=0))
    return left_residues, right_residues, matrix

print("Method ready: residue_minimum_distances")

def residue_labels(frame):
    """Display author residue identifiers without discarding insertion codes."""
    return [f'{r.chain}:{r.residue}{int(r.resnum)}{r.icode}' for r in frame.itertuples()]

print("Method ready: residue_labels")

def pdb_atom_text(atoms):
    """Serialize an explicitly selected atom table for 3D display, without reorienting it."""
    lines = []
    for r in atoms.itertuples():
        lines.append(f'{r.record:<6}{r.serial:5d} {r.atom:>4} {r.residue:>3} {r.chain:1}{r.resnum:4d}{r.icode:1}   '
                     f'{r.x:8.3f}{r.y:8.3f}{r.z:8.3f}{r.occupancy:6.2f}{r.bfactor:6.2f}          {r.element:>2}  ')
    return '\n'.join(lines) + '\nEND\n'

print("Method ready: pdb_atom_text")

KEY = "imatinib"
data = snapshot(KEY)
atoms = parse_pdb(data["pdb"])
heavy = atoms[(atoms.occupancy > 0) & ~atoms.element.isin(["H", "D"])].copy()
assert heavy.serial.is_unique
ca = heavy[(heavy.record == "ATOM") & (heavy.atom == "CA")]
identity = []
for entity in data["entities"]:
    for chain in entity["entity_poly"]["pdbx_strand_id"].split(","):
        identity.append(dict(chain=chain, molecule=entity["rcsb_polymer_entity"]["pdbx_description"],
            deposited_sequence_length=entity["entity_poly"]["rcsb_sample_sequence_length"],
            observed_CA_residues=int((ca.chain == chain).sum())))
display(pd.DataFrame(identity))
print("PDB:", data["pdb_id"], "; resolution (Å):", data["entry"]["rcsb_entry_info"]["resolution_combined"])
print("Revision:", data["entry"]["rcsb_accession_info"]["revision_date"])
print("Selected atom rows:", len(atoms), "; usable heavy atoms:", len(heavy),
      "; zero occupancy:", int((atoms.occupancy == 0).sum()))
print("Non-polymer components:", sorted(atoms.loc[atoms.record == "HETATM", "residue"].unique()))
compound = data["compound"]["chem_comp"]
display(pd.DataFrame([dict(component=compound["id"], formula=compound["formula"],
    formula_weight_Da=compound["formula_weight"], formal_charge=compound["pdbx_formal_charge"])]))
drug = heavy[(heavy.chain == "A") & (heavy.residue == "STI")].copy()
target = heavy[(heavy.chain == "A") & (heavy.record == "ATOM")].copy()
assert len(drug) == 37 and not target.empty
print("Bound drug heavy atoms:", len(drug), "; protein heavy atoms:", len(target))

drug_view = py3Dmol.view(width="100%", height=350)
drug_view.addModel(pdb_atom_text(drug), "pdb")
drug_view.setStyle({}, {"stick":{"colorscheme":"orangeCarbon"}})
drug_view.setHoverable({}, True, "function(a,v){if(!a.label)a.label=v.addLabel(a.atom,{position:a,backgroundColor:'white',fontColor:'black'});}",
    "function(a,v){if(a.label){v.removeLabel(a.label);delete a.label;}}")
drug_view.setBackgroundColor("white"); drug_view.zoomTo(); drug_view.zoom(.7)
display(HTML(drug_view._make_html()))
complex_view = py3Dmol.view(width="100%", height=440)
complex_view.addModel(pdb_atom_text(pd.concat([target, drug])), "pdb")
complex_view.setStyle({"chain":"A"}, {"cartoon":{"color":"#18785f"}})
complex_view.setStyle({"resn":"STI"}, {"stick":{"colorscheme":"orangeCarbon"}})
complex_view.setBackgroundColor("white"); complex_view.zoomTo(); complex_view.zoom(.7)
display(HTML(complex_view._make_html()))

CUTOFF = 4.0
left_residues, right_residues, distances = residue_minimum_distances(drug, target)
target_profile = right_residues.copy()
target_profile["label"] = residue_labels(right_residues)
target_profile["min_distance_A"] = distances.min(axis=0)
pocket = target_profile[target_profile.min_distance_A <= CUTOFF].sort_values("min_distance_A", kind="stable")
print("Pocket residues:", len(pocket))
display(pocket.head(12))
fig, ax = plt.subplots(figsize=(8,5))
if len(pocket):
    ax.barh(pocket.label.iloc[::-1], pocket.min_distance_A.iloc[::-1], color="#18785f")
else:
    ax.text(.5,.5,"No residues within cutoff",ha="center",transform=ax.transAxes)
ax.axvline(CUTOFF, color="#bc6541", linestyle="--", label=f"Cutoff {CUTOFF:g} Å")
ax.set(xlabel="Minimum heavy-atom distance (Å)", title="1IEP: ABL residues near imatinib")
ax.legend(); plt.tight_layout(); plt.show()

atom_to_target = cdist(drug[["x","y","z"]].to_numpy(), target[["x","y","z"]].to_numpy())
groups = list(target.reset_index(drop=True).groupby(["chain","resnum","icode","residue"], sort=False).indices.values())
atom_residue = np.column_stack([atom_to_target[:,indices].min(axis=1) for indices in groups])
assert np.allclose(atom_residue.min(axis=0), target_profile.min_distance_A)
selected = np.flatnonzero(target_profile.min_distance_A.to_numpy() <= CUTOFF)
fig, ax = plt.subplots(figsize=(10,7))
if len(selected):
    image = ax.imshow(atom_residue[:,selected], aspect="auto", cmap="viridis_r", vmin=2, vmax=8)
    ax.set(xticks=range(len(selected)), xticklabels=target_profile.label.iloc[selected],
           yticks=range(len(drug)), yticklabels=drug.atom)
    plt.setp(ax.get_xticklabels(), rotation=70, ha="right")
    fig.colorbar(image, ax=ax, label="Minimum atom-to-residue distance (Å; scale capped at 8)")
else:
    ax.text(.5,.5,"Empty pocket at this cutoff",ha="center",transform=ax.transAxes)
ax.set(xlabel="ABL author residue ID", ylabel="STI drug atom name", title="Observed drug–pocket geometry")
plt.tight_layout(); plt.show()

target_ca = target[target.atom == "CA"].copy()
target_ca["label"] = residue_labels(target_ca)
pocket_ca = target_ca[target_ca.label.isin(pocket.label)]
fig = plt.figure(figsize=(8,5)); ax = fig.add_subplot(111, projection="3d")
ax.scatter(target_ca.x,target_ca.y,target_ca.z,s=8,c="#18785f",alpha=.25,label="ABL Cα")
ax.scatter(drug.x,drug.y,drug.z,s=25,c="#bc6541",label="Imatinib heavy atoms")
ax.scatter(pocket_ca.x,pocket_ca.y,pocket_ca.z,s=30,c="#9974a5",label="Pocket residue Cα")
ax.set(xlabel="x (Å)",ylabel="y (Å)",zlabel="z (Å)",title="Sequence-distant residues form a drug pocket")
xyz=target[["x","y","z"]].to_numpy(); ax.set_box_aspect(np.ptp(xyz,axis=0))
ax.legend(loc="upper left"); plt.tight_layout(); plt.show()
target_labels = pd.Series(residue_labels(target), index=target.index)
pocket_atoms = target[target_labels.isin(pocket.label)]
pocket_view = py3Dmol.view(width="100%", height=440)
pocket_view.addModel(pdb_atom_text(pd.concat([pocket_atoms, drug])), "pdb")
pocket_view.setStyle({}, {"stick":{"color":"#9974a5","radius":.14}})
pocket_view.setStyle({"resn":"STI"}, {"stick":{"colorscheme":"orangeCarbon","radius":.2}})
for row in pocket.iloc[:4].itertuples():
    position = pocket_ca[pocket_ca.label == row.label]
    if len(position):
        r = position.iloc[0]
        pocket_view.addLabel(row.label, {"position":{"x":float(r.x),"y":float(r.y),"z":float(r.z)},
            "backgroundColor":"white","fontColor":"black","fontSize":12,"backgroundOpacity":.8})
pocket_view.setBackgroundColor("white"); pocket_view.zoomTo(); pocket_view.zoom(.7)
display(HTML(pocket_view._make_html()))
print("Close-up atoms:", len(pocket_atoms)+len(drug))

threshold_table = pd.DataFrame([dict(cutoff_A=t,
    target_residues=int((target_profile.min_distance_A <= t).sum())) for t in [0.01,3.,3.5,4.,4.5,5.,6.]])
assert np.all(np.diff(threshold_table.target_residues) >= 0)
display(threshold_table)
fig,ax=plt.subplots(figsize=(7,3))
ax.plot(threshold_table.cutoff_A,threshold_table.target_residues,"o-",color="#416aa6")
ax.axvline(CUTOFF,ls="--",color="#bc6541")
ax.set(xlabel="Geometric cutoff (Å)",ylabel="Unique ABL residues",title="Definition sensitivity; not dose–response")
plt.tight_layout(); plt.show()

# Independent atom-pair calculation for the closest residue pair.
i, j = np.unravel_index(np.argmin(distances), distances.shape)
def residue_atoms(frame, row):
    return frame[(frame.chain == row.chain) & (frame.resnum == row.resnum) &
                 (frame.icode == row.icode) & (frame.residue == row.residue)]
left_test = residue_atoms(drug, left_residues.iloc[i])
right_test = residue_atoms(target, right_residues.iloc[j])
manual_minimum = min(float(np.linalg.norm(a-b))
    for a in left_test[["x","y","z"]].to_numpy() for b in right_test[["x","y","z"]].to_numpy())
assert np.isclose(manual_minimum, distances[i,j], atol=1e-12)
# A common rigid translation cannot change any distance.
shifted_drug = drug.copy(); shifted_target = target.copy()
shifted_drug[["x","y","z"]] += np.array([10., -7., 3.])
shifted_target[["x","y","z"]] += np.array([10., -7., 3.])
_, _, shifted_distances = residue_minimum_distances(shifted_drug, shifted_target)
assert np.allclose(distances, shifted_distances, atol=1e-12)
# Empty selections fail explicitly; a zero-contact threshold is a valid result.
try:
    residue_minimum_distances(drug.iloc[:0], target)
except ValueError:
    print("Empty selection: reported correctly")
else:
    raise AssertionError("Empty selection was accepted")
assert not (distances <= 0.01).any()
assert np.isfinite(distances).all() and (distances >= 0).all()
print("PASS: independent atom minimum, rigid translation, finite distances, empty selections and empty contacts")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
target_profile.to_csv(OUTPUT_DIR / "target-residue-distances.csv", index=False)
threshold_table.to_csv(OUTPUT_DIR / "threshold-sensitivity.csv", index=False)
pd.DataFrame(distances, index=residue_labels(left_residues), columns=residue_labels(right_residues)).to_csv(
    OUTPUT_DIR / "residue-minimum-distances.csv")
(OUTPUT_DIR / "selected-complex.pdb").write_text(pdb_atom_text(pd.concat([target, drug])), encoding="ascii")
summary = dict(pdb_id=data["pdb_id"], cutoff_angstrom=CUTOFF,
    atom_selection="positive occupancy heavy atoms; parser selects one altloc per atom",
    target_chains=sorted(target.chain.unique()), drug_chains=sorted(drug.chain.unique()),
    target_residues_within_cutoff=int((target_profile.min_distance_A <= CUTOFF).sum()),
    contacting_residue_pairs=int((distances <= CUTOFF).sum()),
    method="minimum atom distance, not affinity, interaction energy or buried surface area")
(OUTPUT_DIR / "analysis-summary.json").write_text(json.dumps(summary, indent=2)+"\n")
print("Exported:", sorted(p.name for p in OUTPUT_DIR.iterdir()))

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://data.rcsb.org/rest/v1/core/entry/1IEP',timeout=20)
        response.raise_for_status()
        print("Source:",response.url)
        print(response.text[:700])
    except requests.RequestException as error:
        print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
    print("Online extension skipped; embedded core data unchanged")