How does an antibody meet a protein surface?

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='08-trastuzumab-interface'
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={'trastuzumab': '147b0a3329c44a3ac92e5eda4c4e0e5b18511bf95c9938e3e0936567f9c252b5'}
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 = "trastuzumab"
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()))
drug = heavy[(heavy.record == "ATOM") & heavy.chain.isin(["A","B"])].copy()
target = heavy[(heavy.record == "ATOM") & (heavy.chain == "C")].copy()
assert set(drug.chain) == {"A","B"} and not target.empty
print("Fab heavy atoms:", len(drug), "; HER2 heavy atoms:", len(target))
print("Residues carrying insertion codes:",
      len(drug.loc[drug.icode != "", ["chain","resnum","icode"]].drop_duplicates()))

viewer = py3Dmol.view(width="100%",height=480)
viewer.addModel(pdb_atom_text(pd.concat([target,drug])),"pdb")
for chain,color in [("C","#18785f"),("A","#bc6541"),("B","#416aa6")]:
    viewer.setStyle({"chain":chain},{"cartoon":{"color":color}})
viewer.setBackgroundColor("white");viewer.zoomTo();viewer.zoom(.7)
display(HTML(viewer._make_html()))
fig=plt.figure(figsize=(8,5));ax=fig.add_subplot(111,projection="3d")
for chain,color,label in [("C","#18785f","HER2 C"),("A","#bc6541","Fab light A"),("B","#416aa6","Fab heavy B")]:
    frame=ca[ca.chain==chain]
    ax.scatter(frame.x,frame.y,frame.z,s=10,c=color,alpha=.6,label=label)
ax.set(xlabel="x (Å)",ylabel="y (Å)",zlabel="z (Å)",title="1N8Z: an antibody meets a target surface")
ax.set_box_aspect(np.ptp(pd.concat([target,drug])[["x","y","z"]].to_numpy(),axis=0))
ax.legend(loc="upper left");plt.tight_layout();plt.show()

CUTOFF=4.0
left_residues,right_residues,distances=residue_minimum_distances(drug,target)
left_profile=left_residues.copy();left_profile["label"]=residue_labels(left_residues)
left_profile["min_distance_A"]=distances.min(axis=1)
target_profile=right_residues.copy();target_profile["label"]=residue_labels(right_residues)
target_profile["min_distance_A"]=distances.min(axis=0)
contact_mask=distances<=CUTOFF
contributions=pd.DataFrame([dict(Fab_chain=chain,
    unique_Fab_residues=int(((left_profile.chain==chain)&(left_profile.min_distance_A<=CUTOFF)).sum()),
    residue_pair_edges=int(contact_mask[left_profile.chain.to_numpy()==chain].sum())) for chain in ["A","B"]])
display(contributions)
print("Unique HER2 interface residues:",int((target_profile.min_distance_A<=CUTOFF).sum()),
      "; residue-pair edges:",int(contact_mask.sum()))
display(target_profile[target_profile.min_distance_A<=CUTOFF].sort_values("min_distance_A").head(12))
fig,ax=plt.subplots(figsize=(6,3))
ax.bar(["Fab light A","Fab heavy B"],contributions.unique_Fab_residues,color=["#bc6541","#416aa6"])
ax.set(ylabel="Unique Fab residues within cutoff",title="Geometric contribution of each Fab chain (≤4 Å)")
plt.tight_layout();plt.show()

left_indices=np.flatnonzero(contact_mask.any(axis=1))
right_indices=np.flatnonzero(contact_mask.any(axis=0))
fig,ax=plt.subplots(figsize=(10,7))
if len(left_indices) and len(right_indices):
    image=ax.imshow(distances[np.ix_(left_indices,right_indices)],aspect="auto",cmap="viridis_r",vmin=2,vmax=8)
    ax.set(xticks=range(len(right_indices)),xticklabels=target_profile.label.iloc[right_indices],
        yticks=range(len(left_indices)),yticklabels=left_profile.label.iloc[left_indices])
    plt.setp(ax.get_xticklabels(),rotation=70,ha="right")
    fig.colorbar(image,ax=ax,label="Minimum heavy-atom distance (Å; scale capped at 8)")
else:
    ax.text(.5,.5,"No interface at this cutoff",ha="center",transform=ax.transAxes)
ax.set(xlabel="HER2 author residue ID",ylabel="Fab author residue ID",title="Antibody–target geometry; not interaction energy")
plt.tight_layout();plt.show()

fig,ax=plt.subplots(figsize=(9,3))
ax.scatter(target_profile.resnum,target_profile.min_distance_A,s=9,c="#18785f",alpha=.6,label="Observed HER2 residues")
interface=target_profile[target_profile.min_distance_A<=CUTOFF]
ax.scatter(interface.resnum,interface.min_distance_A,s=28,c="#9974a5",label="Within cutoff")
ax.axhline(CUTOFF,ls="--",c="#bc6541")
ax.set(xlabel="HER2 deposited author residue number",ylabel="Minimum distance to Fab (Å)",title="A localized protein surface is recognized")
ax.legend();plt.tight_layout();plt.show()

threshold_table=pd.DataFrame([dict(cutoff_A=t,Fab_residues=int((distances.min(axis=1)<=t).sum()),
    HER2_residues=int((distances.min(axis=0)<=t).sum()),residue_pair_edges=int((distances<=t).sum()))
    for t in [0.01,3.,3.5,4.,4.5,5.,6.]])
assert np.all(np.diff(threshold_table[["Fab_residues","HER2_residues","residue_pair_edges"]],axis=0)>=0)
display(threshold_table)
fig,ax=plt.subplots(figsize=(7,3))
for column in ["Fab_residues","HER2_residues","residue_pair_edges"]:
    ax.plot(threshold_table.cutoff_A,threshold_table[column],"o-",label=column)
ax.axvline(CUTOFF,ls="--",c="#bc6541")
ax.set(xlabel="Geometric cutoff (Å)",ylabel="Count (denominator shown in legend)",title="Interface definition sensitivity")
ax.legend();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()))
left_profile.to_csv(OUTPUT_DIR / "Fab-residue-distances.csv",index=False)
rows,cols=np.nonzero(contact_mask)
contact_edges=pd.DataFrame(dict(Fab_residue=np.array(residue_labels(left_residues))[rows],
    HER2_residue=np.array(residue_labels(right_residues))[cols],distance_A=distances[rows,cols]))
contact_edges.to_csv(OUTPUT_DIR / "interface-residue-pairs.csv",index=False)
print("Exported Fab profile and",len(contact_edges),"geometric interface edges")

RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
    import requests
    try:
        response=requests.get('https://data.rcsb.org/rest/v1/core/entry/1N8Z',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")