Bio Structural Biology Interface Analysis

Version Compatibility

Reference examples tested with: biopython 1.83+, numpy 1.26+, freesasa 2.2+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Interface Analysis

"Which residues contact the ligand / the partner chain?" -> Threshold interatomic distances and collect residues within a cutoff.

• Python: Bio.PDB.NeighborSearch(atoms).search_all(cutoff, level='R') or .search(center, cutoff, level='R')

"Compute the buried surface area of this interface" -> Subtract complex SASA from the summed SASA of the isolated partners.

• Python: Bio.PDB.SASA.ShrakeRupley().compute(entity), then BSA = SASA_A + SASA_B - SASA_complex

"Is this interface biological or a crystal-packing artifact?" -> Score interface size and chemistry on the biological assembly, and treat the assignment as a hypothesis.

• Python: Bio.PDB for BSA / H-bonds; PDBePISA for the assembly call (external service)

Governing Principle

A "contact" is not a physical fact, it is a thresholded distance, and the residue count changes with the cutoff. Heavy-atom pairs within 4-5A capture direct van der Waals contact; CA-CA within 8A captures topological proximity (contact maps, coevolution features) but says nothing about side-chain interaction; 3.5-4.0A heavy-atom is the H-bond / salt-bridge regime. Changing 4A to 5A can shift the contact count substantially, so a contact result is meaningless without its atom set and cutoff stated explicitly (Chakrabarti & Janin 2002 Proteins 47:334-343). The trap is reporting "N interface residues" or "N contacts" as if the number were intrinsic.

A contact list is not an interface. The physical interface measure is BURIED SURFACE AREA (BSA, also dSASA): BSA = SASA(part A alone) + SASA(part B alone) - SASA(complex), conventionally halved to report the area buried per partner. All three SASA terms must be computed with identical parameters (same probe radius, radii set, algorithm) or the subtraction is garbage (see geometric-analysis for SASA fundamentals). SASA itself depends on the probe radius (1.4A water default) and the algorithm, so an absolute BSA is only comparable to another BSA computed the same way. Because heavy-atom contacts and BSA are both defined on non-hydrogen atoms, hydrogens are NOT required for either - add them (structure-preparation) only for H-bond/salt-bridge angle geometry, and if H are present keep them consistent across all three SASA terms.

The deepest trap: an interface seen in the deposited ASYMMETRIC UNIT may be a CRYSTAL-PACKING ARTIFACT, not biology. The asymmetric unit is a crystallographic bookkeeping object; the functional molecule is the BIOLOGICAL ASSEMBLY, which may be a subset of the ASU or built from several ASUs by symmetry. Compute interfaces on the biological assembly, not blindly on the ASU (see structure-io for downloading the assembly). PDBePISA (Krissinel & Henrick 2007 J Mol Biol 372:774-797) predicts the biological assembly and scores interface stability, but it recovers the correct assembly only ~80-90% of the time and has known false positives, so "biological interface" is a HYPOTHESIS. Larger BSA, more H-bonds and salt bridges, shape complementarity, and evolutionary conservation of interface residues each raise confidence, but each is probabilistic, not proof (Levy 2010 J Mol Biol 403:660-670). Corroborate anything load-bearing with solution data (SEC-MALS, SAXS, native MS).

Decision: contact / interface definition

The one-line rule: heavy-atom 4-5A answers "who touches"; CA-CA 8A answers "who is near"; BSA answers "how big is the interface". State the atom set and cutoff every time.

Decision: biological interface vs crystal contact

Every row is probabilistic. Interface size (BSA) is the single most-used signal, but small biological interfaces (transient/weak complexes) and large crystal contacts both exist, so no one number settles it.

Contact residues between two chains

Goal: List the residues of chain A and chain B that form the interface, under an explicit cutoff.

Approach: Build one KD-tree over the interface atoms, query all close pairs at residue level, and keep pairs whose two residues belong to different chains. Heavy-atom cutoff 4.5A (midpoint of the 4-5A vdW-contact regime; excludes H so it is robust to whether H atoms were modeled).

from Bio.PDB import PDBParser, NeighborSearch, Selection

parser = PDBParser(QUIET=True)
structure = parser.get_structure('complex', 'complex.pdb')
model = structure[0]

cutoff = 4.5  # heavy-atom contact; 4-5A captures direct vdW contact, state it always
atoms = [a for a in model.get_atoms() if a.element != 'H']
ns = NeighborSearch(atoms)

interface_a, interface_b = set(), set()
for res1, res2 in ns.search_all(cutoff, level='R'):
    c1, c2 = res1.get_parent().id, res2.get_parent().id
    if c1 == 'A' and c2 == 'B':
        interface_a.add(res1); interface_b.add(res2)
    elif c1 == 'B' and c2 == 'A':
        interface_b.add(res1); interface_a.add(res2)

print(f'Chain A interface residues ({cutoff}A): {len(interface_a)}')
print(f'Chain B interface residues ({cutoff}A): {len(interface_b)}')

Ligand-contact (binding-site / epitope) residues

Goal: Identify the protein residues lining a bound ligand or the residues an antibody contacts (structural epitope).

Approach: Select the ligand atoms (a HETATM group, hetflag starts with 'H_'), search protein atoms within the cutoff of each, collect unique parent residues. The same pattern with two protein chains yields a structural epitope.

from Bio.PDB import PDBParser, NeighborSearch

parser = PDBParser(QUIET=True)
structure = parser.get_structure('complex', 'complex.pdb')
model = structure[0]

ligand_resname = 'ATP'  # target HETATM group
cutoff = 4.5

ligand_atoms = [a for r in model.get_residues() if r.resname == ligand_resname
                for a in r if a.element != 'H']
protein_atoms = [a for a in model.get_atoms()
                 if a.element != 'H' and a.get_parent().id[0] == ' ']
ns = NeighborSearch(protein_atoms)

pocket = set()
for a in ligand_atoms:
    for res in ns.search(a.coord, cutoff, level='R'):
        pocket.add((res.get_parent().id, res.id[1], res.resname))

for chain, num, name in sorted(pocket):
    print(f'{chain} {name}{num}')

Buried surface area (BSA / dSASA)

Goal: Quantify the physical size of a two-chain interface as area buried on complex formation.

Approach: Compute SASA on the intact complex, then on each chain in isolation (same ShrakeRupley settings), and take BSA = SASA_A + SASA_B - SASA_complex. Halve for per-partner area. Probe radius 1.4A models a water molecule; keep it identical across all three computations or the subtraction is meaningless.

from Bio.PDB import PDBParser
from Bio.PDB.SASA import ShrakeRupley

parser = PDBParser(QUIET=True)
sr = ShrakeRupley(probe_radius=1.4)  # 1.4A ~ water; MUST match across all three terms

def chain_sasa(path, keep_chains):
    structure = parser.get_structure('s', path)
    model = structure[0]
    for chain in list(model):
        if chain.id not in keep_chains:
            model.detach_child(chain.id)
    sr.compute(model, level='C')
    return sum(chain.sasa for chain in model)

sasa_complex = chain_sasa('complex.pdb', {'A', 'B'})
sasa_a = chain_sasa('complex.pdb', {'A'})
sasa_b = chain_sasa('complex.pdb', {'B'})

bsa_total = sasa_a + sasa_b - sasa_complex
print(f'Total buried surface area: {bsa_total:.0f} A^2')
print(f'Per partner: {bsa_total / 2:.0f} A^2')  # convention: split half to each side

For Lee-Richards SASA or full control of the radii set and probe, use freesasa instead of ShrakeRupley (Mitternacht 2016 F1000Research 5:189); Bio.PDB provides only Shrake-Rupley. Compute all three terms in the same tool.