Bio Structural Biology Geometric Analysis

Version Compatibility

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

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.

Geometric Analysis

"Calculate the RMSD between two conformations" -> superimpose on a defined atom correspondence, then report deviation.

• Python: Bio.PDB.Superimposer (SVD/Kabsch); Bio.PDB.qcprot for speed in tight loops "How buried is this residue?" -> compute solvent accessible surface area and normalize to a per-residue maximum.
• Python: Bio.PDB.SASA.ShrakeRupley, then relative SASA against a max-ASA scale

Governing Principle: the metric IS the question

RMSD is NOT a property of two structures. It is a property of two structures GIVEN a superposition AND an atom selection - change either and the number changes. Report what was aligned (CA-only? a defined core? all-atom?) or the number is uninterpretable.

Global all-atom RMSD is a mean of SQUARED per-atom deviations after a least-squares rigid-body fit, so it is dominated by the worst-fitting atoms. A structure whose 150-residue core is essentially identical but whose 10-residue loop or a hinge-rotated domain swings out by 15A reports a "bad" whole-molecule RMSD (often 4-8A) that hides a near-perfect core. RMSD is also length-dependent (longer proteins accumulate larger RMSD for the same local quality) and is NOT a cross-protein similarity metric - it is only meaningful when a genuine 1:1 correspondence exists (same protein, two states; a model vs its native).

Superimposer requires an equal-length, ORDERED atom-to-atom correspondence. Feeding it mismatched or unequal atom lists is the classic error; it computes the optimal fit (Kabsch via SVD), it does NOT solve which atom maps to which. Structures with different sequences need a structure-based alignment FIRST to establish the correspondence (see alignment/structural-alignment), then superposition.

SASA depends on the PROBE RADIUS (1.4A water is a convention, not a constant of nature), the algorithm (Bio.PDB ShrakeRupley is Shrake-Rupley only; freesasa offers Lee-Richards and LCPO), and whether hydrogens are present. A SASA number without its probe radius is meaningless, and absolute SASA in A^2 is not portable across tools. Prefer RELATIVE SASA (residue SASA / max-ASA of that residue type) using the Tien et al 2013 max-ASA scale.

Backbone phi/psi and omega/cis-peptides are a VALIDATION signal, not a description: a residue in a sterically disallowed Ramachandran region usually means a modeling error, not exotic biology. This skill computes the angles; interpreting outliers as quality flags belongs to structural-biology/structure-validation.

Decision: which comparison metric

For cross-protein or fold-similarity work, do not stretch RMSD - route to alignment/structural-alignment (TM-align / Foldseek / DALI). For Ramachandran/omega as a quality gate, route to structural-biology/structure-validation.

Distance Between Atoms

from Bio.PDB import PDBParser
import numpy as np

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

chain = structure[0]['A']
atom1, atom2 = chain[100]['CA'], chain[200]['CA']

distance = atom1 - atom2                              # Atom subtraction returns the distance directly
print(f'Distance: {distance:.2f} A')
print(np.linalg.norm(atom1.coord - atom2.coord))     # Equivalent via numpy on the .coord arrays

Distance Matrix

import numpy as np
from Bio.PDB import PDBParser

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

ca_atoms = [r['CA'] for r in structure.get_residues() if r.has_id('CA') and r.id[0] == ' ']  # id[0]==' ' drops waters/hetero
n = len(ca_atoms)

dist = np.zeros((n, n))
for i in range(n):
    for j in range(i + 1, n):
        dist[i, j] = dist[j, i] = ca_atoms[i] - ca_atoms[j]
print(f'Distance matrix: {dist.shape}')

Bond Angle

import numpy as np
from Bio.PDB import PDBParser, calc_angle

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

res = structure[0]['A'][100]
angle = calc_angle(res['N'].get_vector(), res['CA'].get_vector(), res['C'].get_vector())  # calc_angle needs Vector, not .coord
print(f'N-CA-C angle: {np.degrees(angle):.1f} deg')

Backbone Dihedrals (phi / psi)

import numpy as np
from Bio.PDB import PDBParser, calc_dihedral

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')
chain = structure[0]['A']

prev, curr, nxt = chain[99], chain[100], chain[101]
phi = calc_dihedral(prev['C'].get_vector(), curr['N'].get_vector(), curr['CA'].get_vector(), curr['C'].get_vector())
psi = calc_dihedral(curr['N'].get_vector(), curr['CA'].get_vector(), curr['C'].get_vector(), nxt['N'].get_vector())
print(f'phi={np.degrees(phi):.1f}  psi={np.degrees(psi):.1f}')

Ramachandran Angles for All Residues

import numpy as np
from Bio.PDB import PDBParser, PPBuilder

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

ppb = PPBuilder()                                    # PPBuilder builds peptides from connectivity, so chain breaks split them
rama = []
for pp in ppb.build_peptides(structure):
    for res, (phi, psi) in zip(pp, pp.get_phi_psi_list()):
        if phi is not None and psi is not None:      # None at termini and chain breaks by design - skip, do not fabricate
            rama.append((res.resname, np.degrees(phi), np.degrees(psi)))
print(f'{len(rama)} residues with phi/psi')

Outliers in disallowed regions are usually refinement errors, not biology - interpret them as a quality gate in structural-biology/structure-validation.

Chi Angles (Sidechain Dihedrals)

import numpy as np
from Bio.PDB import PDBParser, calc_dihedral

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')
res = structure[0]['A'][100]

if res.has_id('CB') and res.has_id('CG'):            # Gly/Ala lack CB/CG; chi atom quartets are residue-type specific
    chi1 = calc_dihedral(res['N'].get_vector(), res['CA'].get_vector(), res['CB'].get_vector(), res['CG'].get_vector())
    print(f'Chi1: {np.degrees(chi1):.1f} deg')

Superimposing Structures and RMSD

from Bio.PDB import PDBParser, Superimposer

parser = PDBParser(QUIET=True)
ref = parser.get_structure('ref', 'reference.pdb')
mob = parser.get_structure('mobile', 'mobile.pdb')

ref_ca = [r['CA'] for r in ref.get_residues() if r.has_id('CA') and r.id[0] == ' ']
mob_ca = [r['CA'] for r in mob.get_residues() if r.has_id('CA') and r.id[0] == ' ']
n = min(len(ref_ca), len(mob_ca))                    # Superimposer needs EQUAL-LENGTH ORDERED lists; it does NOT solve correspondence
ref_ca, mob_ca = ref_ca[:n], mob_ca[:n]              # Naive truncation is only valid when residues already correspond 1:1

sup = Superimposer()
sup.set_atoms(ref_ca, mob_ca)                        # Optimal rigid-body fit via SVD/Kabsch
print(f'RMSD (CA): {sup.rms:.2f} A')
rotation, translation = sup.rotran                   # The fitted transform, for reuse on other atoms
sup.apply(mob.get_atoms())                           # Mutates mob in place - copy first if the originals are still needed

QCP alternative for speed in tight loops (MD, all-vs-all): from Bio.PDB.qcprot import QCPSuperimposer (module Bio.PDB.qcprot; historically Bio.PDB.QCPSuperimposer), same set_atoms / .rms / .rotran / apply interface and identical optimum.

Per-Residue Deviation After Superposition

Fitting minimizes the squared mean, so a global scalar hides where the structures actually differ. A per-residue deviation plot exposes the outlier domination directly.

import numpy as np
from Bio.PDB import PDBParser, Superimposer

parser = PDBParser(QUIET=True)
ref = parser.get_structure('ref', 'reference.pdb')
mob = parser.get_structure('mobile', 'mobile.pdb')

ref_ca = [r['CA'] for r in ref.get_residues() if r.has_id('CA') and r.id[0] == ' ']
mob_ca = [r['CA'] for r in mob.get_residues() if r.has_id('CA') and r.id[0] == ' ']
n = min(len(ref_ca), len(mob_ca))
ref_ca, mob_ca = ref_ca[:n], mob_ca[:n]

sup = Superimposer()
sup.set_atoms(ref_ca, mob_ca)
sup.apply([a for a in mob_ca])                        # Move only the paired CA set into the fitted frame
deviation = np.array([r - m for r, m in zip(ref_ca, mob_ca)])
print(f'core (<2A) residues: {(deviation < 2.0).sum()} / {n}')   # 2A is a common rigid-core cutoff, not a law
print(f'max deviation: {deviation.max():.2f} A at index {deviation.argmax()}')

Center of Mass and Radius of Gyration

import numpy as np
from Bio.PDB import PDBParser

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

atoms = list(structure.get_atoms())
coords = np.array([a.coord for a in atoms])
masses = np.array([{'C': 12.0, 'N': 14.0, 'O': 16.0, 'S': 32.0, 'H': 1.0}.get(a.element, 12.0) for a in atoms])

com = (masses[:, None] * coords).sum(axis=0) / masses.sum()
print(f'Center of mass: {com}')

rg = np.sqrt(np.mean(np.sum((coords - coords.mean(axis=0)) ** 2, axis=1)))  # Unweighted radius of gyration
print(f'Radius of gyration: {rg:.2f} A')

Vector Operations

from Bio.PDB import PDBParser

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

v1 = structure[0]['A'][100]['CA'].get_vector()
v2 = structure[0]['A'][101]['CA'].get_vector()

diff = v2 - v1
print(f'length: {diff.norm():.2f}  unit: {diff.normalized()}')
cross = v1 ** v2                                      # ** is cross product on Vector objects
dot = v1 * v2                                         # * is dot product on Vector objects

Solvent Accessible Surface Area (SASA)

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

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')

sr = ShrakeRupley(probe_radius=1.40, n_points=100)   # 1.4A = water radius (convention); a SASA number is meaningless without its probe radius
sr.compute(structure, level='R')                     # level R attaches .sasa on each residue; children sum to parents
print(f'total SASA: {sum(r.sasa for r in structure.get_residues() if hasattr(r, "sasa")):.1f} A^2')

Relative SASA and Burial

Absolute SASA is not portable across tools. For burial, normalize to a per-residue maximum (Tien et al 2013 theoretical Gly-X-Gly max-ASA).

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

MAX_ASA = {'ALA': 129.0, 'ARG': 274.0, 'ASN': 195.0, 'ASP': 193.0, 'CYS': 167.0, 'GLU': 223.0, 'GLN': 225.0, 'GLY': 104.0, 'HIS': 224.0, 'ILE': 197.0, 'LEU': 201.0, 'LYS': 236.0, 'MET': 224.0, 'PHE': 240.0, 'PRO': 159.0, 'SER': 155.0, 'THR': 172.0, 'TRP': 285.0, 'TYR': 263.0, 'VAL': 174.0}

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')
ShrakeRupley().compute(structure, level='R')

buried = 0
for res in structure.get_residues():
    if res.resname in MAX_ASA and hasattr(res, 'sasa'):
        rsa = res.sasa / MAX_ASA[res.resname]
        if rsa < 0.20:                               # RSA < 0.20 is the common buried heuristic (a rule of thumb, not a law)
            buried += 1
print(f'buried residues (RSA < 0.20): {buried}')

Secondary-Structure Assignment (DSSP)

Secondary structure is an INTERPRETATION, not a value stored in the file: DSSP, STRIDE, and P-SEA legitimately disagree by 1-2 residues at helix and strand termini, so name the tool and version and never mix assignments from two tools in one analysis. DSSP places the backbone amide hydrogen itself and scores an electrostatic H-bond energy, so it needs no explicit hydrogens, and it processes only the FIRST model of an ensemble. The binary was renamed dssp -> mkdssp (v4) and must be installed separately.

from Bio.PDB import PDBParser, DSSP

parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', 'protein.pdb')
model = structure[0]                                  # DSSP runs on ONE model only

dssp = DSSP(model, 'protein.pdb', dssp='mkdssp')      # pass the current binary name explicitly
codes = [dssp[k][2] for k in dssp.keys()]             # 8-state: H G I helix, E B strand, T S P - other
helix = sum(c in 'HGI' for c in codes)
strand = sum(c in 'EB' for c in codes)
print(f'helix {helix}, strand {strand}, other {len(codes) - helix - strand} of {len(codes)}')

Common Errors

Related Skills

• structure-io - Parse and write PDB/mmCIF structure files
• structure-navigation - Walk chains, residues, atoms; handle altlocs and disordered residues
• structure-modification - Transform coordinates and edit structures in place
• structural-biology/interface-analysis - Residue contacts, contact maps, and buried-surface interface analysis (NeighborSearch)
• structural-biology/structure-validation - Ramachandran and omega/cis-peptide outliers as a quality gate
• structural-biology/alphafold-predictions - overlay pLDDT/PAE when a compared structure is a predicted model
• structural-biology/modern-structure-prediction - reconcile predicted models via pLDDT/PAE/pTM before RMSD claims
• alignment/structural-alignment - Cross-protein fold comparison and correspondence (TM-align, Foldseek, DALI)

References

• Cock PJA, et al. (2009) Biopython. Bioinformatics 25(11):1422-1423.
• Kabsch W (1976) A solution for the best rotation to relate two sets of vectors. Acta Crystallogr A 32:922-923.
• Theobald DL (2005) Rapid calculation of RMSDs using a quaternion-based characteristic polynomial. Acta Crystallogr A 61(4):478-480.
• Zhang Y, Skolnick J (2004) Scoring function for automated assessment of protein structure template quality. Proteins 57(4):702-710.
• Xu J, Zhang Y (2010) How significant is a protein structure similarity with TM-score = 0.5? Bioinformatics 26(7):889-895.
• Mariani V, Biasini M, Barbato A, Schwede T (2013) lDDT: a local superposition-free score for comparing protein structures and models. Bioinformatics 29(21):2722-2728.
• Shrake A, Rupley JA (1973) Environment and exposure to solvent of protein atoms. Lysozyme and insulin. J Mol Biol 79(2):351-371.
• Tien MZ, Meyer AG, Sydykova DK, Spielman SJ, Wilke CO (2013) Maximum allowed solvent accessibilities of residues in proteins. PLoS ONE 8(11):e80635.