Contacts, Interfaces, and Hydrogen Bonds

from Bio.PDB import PDBParser, NeighborSearch

parser = PDBParser(QUIET=True)
structure = parser.get_structure("3hfm", "3hfm.pdb")  # HyHEL-10 Fab + lysozyme
model = structure[0]

# Protein heavy atoms only (drop hydrogens, water, and hetero groups)
atoms = [a for a in model.get_atoms()
         if a.element != "H" and a.get_parent().id[0] == " "]
ns = NeighborSearch(atoms)                 # builds a KD-tree over the coordinates

# Every heavy-atom pair within 4.0 A, answered by radius queries on the tree
pairs = ns.search_all(4.0, level="A")
contacts = [(a, b) for a, b in pairs if a.get_parent() is not b.get_parent()]

print(f"{len(atoms)} protein heavy atoms")
print(f"{len(contacts)} inter-residue contacts <= 4.0 A")

import numpy as np
from Bio.PDB.vectors import calc_angle

DONORS = {"N", "NE", "NH1", "NH2", "ND1", "ND2", "NE1", "NE2", "NZ",
          "OG", "OG1", "OH", "SG"}
ACCEPTORS = {"O", "OD1", "OD2", "OE1", "OE2", "ND1", "NE2",
             "OG", "OG1", "OH", "SD"}

hbonds = []
for a, b in ns.search_all(3.5, level="A"):        # D...A heavy-atom cutoff 3.5 A
    if a.element not in ("N", "O", "S") or b.element not in ("N", "O", "S"):
        continue
    for donor, acc in ((a, b), (b, a)):
        dres, ares = donor.get_parent(), acc.get_parent()
        if dres is ares:
            continue
        # Peptide-bond geometry puts N(i+1) ~2.3 A from O(i) without any H-bond,
        # so skip sequentially adjacent residues within the same chain
        if dres.get_parent().id == ares.get_parent().id and \
           abs(dres.id[1] - ares.id[1]) <= 1:
            continue
        if donor.get_name() not in DONORS or acc.get_name() not in ACCEPTORS:
            continue
        dist = donor - acc                        # Atom.__sub__ returns distance (A)
        # Bonded hydrogens live in the donor residue; the KD-tree holds heavy
        # atoms only, so query the residue directly. On a protonated model this
        # enables the D-H...A angle >= 120 deg test; on 3hfm it stays n/a.
        hs = [h for h in dres.get_atoms()
              if h.element == "H" and h - donor < 1.2]
        angle = max((np.degrees(calc_angle(donor.get_vector(),
                                           h.get_vector(),
                                           acc.get_vector())) for h in hs),
                    default=None)
        if angle is not None and angle < 120:
            continue
        hbonds.append((donor, acc, dist, angle))

print(f"{len(hbonds)} candidate hydrogen bonds")
for donor, acc, dist, angle in hbonds[:4]:
    dr, ar = donor.get_parent(), acc.get_parent()
    ang = f"{angle:5.1f}" if angle is not None else "  n/a"
    print(f"{dr.resname} {dr.id[1]:>3} {donor.get_name():<4} -> "
          f"{ar.resname} {ar.id[1]:>3} {acc.get_name():<4}  d={dist:.2f}  angle={ang}")

CATIONIC = {"ARG": {"NH1", "NH2", "NE"}, "LYS": {"NZ"}, "HIS": {"ND1", "NE2"}}
ANIONIC  = {"ASP": {"OD1", "OD2"}, "GLU": {"OE1", "OE2"}}
NONPOLAR = {"ALA", "VAL", "LEU", "ILE", "MET", "PHE", "TRP", "PRO"}

salt_bridges, hydrophobic, interface = set(), set(), set()
for a, b in contacts:
    ra, rb = a.get_parent(), b.get_parent()

    # Salt bridge (Barlow & Thornton): N...O <= 4.0 A across opposite charges
    if (a - b) <= 4.0:
        if a.get_name() in CATIONIC.get(ra.resname, ()) and \
           b.get_name() in ANIONIC.get(rb.resname, ()):
            salt_bridges.add((ra.get_full_id()[2:4], rb.get_full_id()[2:4]))
        if b.get_name() in CATIONIC.get(rb.resname, ()) and \
           a.get_name() in ANIONIC.get(ra.resname, ()):
            salt_bridges.add((rb.get_full_id()[2:4], ra.get_full_id()[2:4]))

    # Hydrophobic packing: side-chain C...C between nonpolar residues
    # (exclude backbone C and CA so sequential main-chain carbons don't count)
    if a.element == "C" and b.element == "C" and \
       a.get_name() not in ("C", "CA") and b.get_name() not in ("C", "CA") and \
       ra.resname in NONPOLAR and rb.resname in NONPOLAR:
        hydrophobic.add((ra.get_full_id()[2:4], rb.get_full_id()[2:4]))

    # Interface residues: the contact spans two different chains
    if ra.get_parent().id != rb.get_parent().id:
        interface.add((ra.get_parent().id, ra.resname, ra.id[1]))
        interface.add((rb.get_parent().id, rb.resname, rb.id[1]))

print(f"{len(salt_bridges)} salt bridges, {len(hydrophobic)} hydrophobic pairs")
print(f"{len(interface)} interface residues across "
      f"{len({c for c, _, _ in interface})} chains")
for chain, resn, resi in sorted(interface)[:5]:
    print(f"  chain {chain}  {resn}{resi}")

# Cross-check the same contacts visually in PyMOL.
fetch 3hfm, async=0
remove solvent

# Polar contacts (potential h-bonds) between lysozyme (Y) and Fab (H, L)
distance epitope_hb, chain Y, (chain H or chain L), mode=2

# Interface residues within 4.0 A across the antigen-antibody boundary
select epitope,  byres (chain Y within 4.0 of (chain H or chain L))
select paratope, byres ((chain H or chain L) within 4.0 of chain Y)
show sticks, epitope or paratope
iterate epitope and name CA, print(f"{chain}/{resn}{resi}")