Bio Structural Biology Structure Validation

Version Compatibility

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

MolProbity, phenix (phenix.molprobity, phenix.process_predicted_model), and DSSP (mkdssp) are external CLI tools invoked via subprocess, not pip packages; install them separately and confirm they are on PATH before use.

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

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

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

Structure Validation

"Is this structure good enough to build on?" -> Read the data-quality metadata, run geometry validation, and decide per-region rather than per-file.

• Python: Bio.PDB.MMCIF2Dict (resolution/R-free), Bio.PDB.calc_dihedral (phi/psi/omega), subprocess -> phenix.molprobity (clashscore, rotamer, Ramachandran).

Governing Principle

A coordinate file is an INTERPRETED MODEL fit to experimental data, not the data and not ground truth, and its reliability varies PER-ATOM within one file. The single global resolution is a DATA CEILING on what the map can resolve, never a per-region quality certificate: a 1.5 Angstrom structure can still carry a guesswork surface loop, and a 3.2 Angstrom structure can have a rigid, locally excellent active site. The number that answers "can I trust THIS region" is the LOCAL signal - the per-residue B-factor and the real-space fit (RSRZ) for the residues actually used - not the headline resolution. Validate before any geometric interpretation: clashscore, Ramachandran and rotamer outliers, cis-peptides, and bond/angle deviations tell whether the model is even self-consistent before a distance or an angle computed on it means anything.

R-free, not R-work, is the cross-validation statistic. R-work is computed on the reflections used in refinement and can be driven down by adding parameters (waters, alt-confs, loose restraints) that fit noise; R-free is the same metric on a held-out test set that never entered refinement (Brunger 1992 Nature 355:472-475). The R-free-MINUS-R-work GAP is the overfitting flag - a gap much larger than expected for the resolution signals a model fitting its own noise, and a suspiciously small gap signals test-set leakage (Kleywegt & Brunger 1996 Structure 4:897-904). B-factors conflate genuine thermal motion, static disorder, and MODEL ERROR into one number, so they are a within-structure relative signal ("which atoms are least certain here"), NOT a portable cross-structure dynamics readout; normalize before comparing across structures and exclude high-B atoms (say >60-80 Angstrom^2 at moderate resolution) from precise geometric claims.

For a PREDICTED model the confidence scores bound SELF-CONSISTENCY, NOT correctness. pLDDT is per-residue and PAE is inter-residue confidence in the model's OWN frame (Jumper et al. 2021 Nature 596:583-589); a confident model can be confidently wrong when AlphaFold modeled a monomer/apo state where the truth is a complex/holo/alternative state. Trim low-pLDDT residues and split by PAE-defined domains before docking or molecular replacement (Oeffner et al. 2022 Acta Crystallogr D 78:1303-1314). Cryo-EM has GLOBAL vs LOCAL resolution - flexible peripheries in a big map are often docked hypotheses at 6-8 Angstrom local resolution inside a "2.5 Angstrom" map - and TWO different FSC thresholds: 0.143 on the gold-standard half-map FSC for the resolution claim, 0.5 for map-vs-model agreement (Rosenthal & Henderson 2003 J Mol Biol 333:721-745). An NMR deposit is an ENSEMBLE of models; validate per-model and report the per-residue spread - never average coordinates, which produces a physically impossible structure with distorted bonds and clashes (Montelione et al. 2013 Structure 21:1563-1570).

Decision: what the resolution buys (X-ray / cryo-EM global)

Decision: R-free and the overfitting gap

Absolute R rises with resolution (weaker high-angle reflections are noisier); a gap wider than the row's threshold flags overfitting, a near-zero gap flags work/test contamination. These are heuristics - verify against contemporaneous PDB statistics for the resolution.

Decision: geometry validation targets (MolProbity)

The poor-rotamer target tightened from ~1.0% to 0.3% with the updated reference distributions (Williams 2018). Read the wwPDB report's PERCENTILE sliders alongside the raw numbers: they rank an entry against the whole archive AND against same-resolution entries (Gore et al. 2017 Structure 25:1916-1927).

Decision: how to validate this model (experimental vs predicted fork)

Read Validation Metadata From the mmCIF Header

Goal: Pull resolution, R-work, R-free, and method from a deposited entry and flag the overfitting gap - the numbers Bio.PDB's thin structure.header omits.

Approach: Read the raw mmCIF categories with MMCIF2Dict (which reaches anything in the file), cast the strings, and compare the R-free-minus-R-work gap against a resolution-scaled expectation.

from Bio.PDB.MMCIF2Dict import MMCIF2Dict

def read_refinement_metadata(cif_path):
    meta = MMCIF2Dict(cif_path)
    def first(key):
        val = meta.get(key, ['NA'])[0]
        try:
            return float(val)
        except ValueError:
            return val
    method = meta.get('_exptl.method', ['NA'])[0]
    resolution = first('_refine.ls_d_res_high')
    r_work = first('_refine.ls_R_factor_R_work')
    r_free = first('_refine.ls_R_factor_R_free')
    gap = r_free - r_work if isinstance(r_free, float) and isinstance(r_work, float) else None
    # A gap wider than ~0.05 flags overfitting at typical (~2 A) resolution (Kleywegt & Brunger 1996).
    overfit_flag = gap is not None and gap > 0.05
    return {'method': method, 'resolution': resolution, 'r_work': r_work, 'r_free': r_free, 'gap': gap, 'overfit_flag': overfit_flag}

Sanity-Check B-factors Within One Structure

Goal: Locate the least-certain atoms of THIS model so precise-distance claims avoid them - a relative, within-structure read, never a cross-structure comparison.

Approach: Collect per-residue mean B for a chain, then flag residues whose B sits far above the structure's own median as low-confidence.

import numpy as np
from Bio.PDB import PDBParser

def bfactor_outliers(structure, chain_id, z_cut=2.0):
    residue_b = {}
    for residue in structure[0][chain_id]:
        if residue.id[0] != ' ':  # skip hetero/water; validate the polymer only
            continue
        residue_b[residue.id] = np.mean([a.get_bfactor() for a in residue])
    vals = np.array(list(residue_b.values()))
    median, mad = np.median(vals), np.median(np.abs(vals - np.median(vals))) + 1e-9
    # Robust z on B: |B - median| / (1.4826*MAD); high B marks disorder/model error, not portable dynamics.
    return {rid: b for rid, b in residue_b.items() if (b - median) / (1.4826 * mad) > z_cut}

Flag Ramachandran Outliers and cis Non-Proline Peptides

Goal: Catch backbone geometry that is almost always a modeling error - residues in disallowed phi/psi and cis peptide bonds that are not proline.

Approach: Get per-residue phi/psi from PPBuilder, coarsely classify against the canonical basins, and compute omega (Ca-C-N-Ca) directly to find cis bonds; a cis assignment at a non-proline is a validation flag until proven by density.

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

# Coarse general-allowed basins (deg): alpha, beta/PPII, left-handed. Authoritative
# favored/allowed/outlier percentages need MolProbity's rama8000 contours (Williams 2018);
# this only screens gross outliers to decide whether to run MolProbity.
_BASINS = [(-160, -40, -80, 30), (-180, -40, 90, 180), (30, 90, -30, 90)]

def is_rama_allowed(phi, psi):
    d = np.degrees([phi, psi])
    return any(lo_p <= d[0] <= hi_p and lo_s <= d[1] <= hi_s for lo_p, hi_p, lo_s, hi_s in _BASINS)

def geometry_flags(structure):
    outliers, cis_nonpro = [], []
    ppb = PPBuilder()
    for pp in ppb.build_peptides(structure[0]):
        for residue, (phi, psi) in zip(pp, pp.get_phi_psi_list()):
            if phi is not None and psi is not None and not is_rama_allowed(phi, psi):
                outliers.append((residue.get_parent().id, residue.id[1], residue.resname))
        residues = list(pp)
        for prev, curr in zip(residues, residues[1:]):
            omega = calc_dihedral(prev['CA'].get_vector(), prev['C'].get_vector(), curr['N'].get_vector(), curr['CA'].get_vector())
            # omega ~180 = trans, ~0 = cis; |omega| < 30 deg is cis. cis at non-Pro is rare and usually an error.
            if abs(np.degrees(omega)) < 30 and curr.resname != 'PRO':
                cis_nonpro.append((curr.get_parent().id, curr.id[1], curr.resname))
    return {'rama_outliers': outliers, 'cis_nonproline': cis_nonpro}

Run MolProbity for Authoritative Geometry (the real answer)

Goal: Get archive-calibrated clashscore, rotamer, and Ramachandran outlier percentages rather than the coarse Python screen above.

Approach: Shell out to phenix.molprobity (or the MolProbity web service / molprobity.molprobity), which carries the reference contour and rotamer distributions the Python screen cannot reproduce.

import subprocess

def run_molprobity(model_path):
    # phenix.molprobity writes molprobity.out with clashscore, rotamer_outliers,
    # ramachandran_outliers, ramachandran_favored, molprobity_score. Parse that file.
    result = subprocess.run(['phenix.molprobity', model_path], capture_output=True, text=True)
    return result.stdout

The coarse Python screen decides WHETHER to run MolProbity; MolProbity (Chen et al. 2010 Acta Crystallogr D 66:12-21; Williams et al. 2018 Protein Sci 27:293-315) gives the numbers to report. For a deposited entry, prefer the pre-computed wwPDB validation report (https://files.rcsb.org/pub/pdb/validation_reports/<xy>/<id>/<id>_validation.xml.gz) - it already carries the percentile sliders and per-residue RSRZ.

Validate a Predicted Model Before Docking or MR

Goal: Turn an AlphaFold/ESMFold model into a trustworthy input by trimming low-confidence residues and reading inter-domain confidence - because pLDDT/PAE bound self-consistency, not correctness.

Approach: Read pLDDT from the B-factor column into bands, read the PAE matrix for domain segmentation, then hand the raw file to phenix.process_predicted_model (which converts pLDDT to a pseudo-B, trims, and splits by PAE-defined domains).

import json
import numpy as np
import subprocess
from Bio.PDB import MMCIFParser

_PLDDT_BANDS = [(90, 'very_high'), (70, 'confident'), (50, 'low'), (0, 'very_low')]

def plddt_bands(cif_path):
    # pLDDT rides in the B-factor column but is confidence (high = good), OPPOSITE polarity to a real B-factor.
    structure = MMCIFParser(QUIET=True).get_structure('pred', cif_path)
    counts = {label: 0 for _, label in _PLDDT_BANDS}
    for residue in structure[0].get_residues():
        if 'CA' not in residue:
            continue
        score = residue['CA'].get_bfactor()
        counts[next(label for cut, label in _PLDDT_BANDS if score >= cut)] += 1
    return counts  # a long very_low run usually flags an intrinsically disordered region, not an error

def pae_interdomain_confident(pae_json, block_a, block_b, cutoff=5.0):
    # Off-diagonal PAE (A) between two domain blocks; low = relative orientation trusted, high = independent guess.
    pae = np.array(json.load(open(pae_json))[0]['predicted_aligned_error'])
    return pae[np.ix_(block_a, block_b)].mean() < cutoff

def process_for_mr(model_path):
    # Trims below ~0.7 fractional pLDDT, converts pLDDT->pseudo-B, splits into PAE-defined domains (Oeffner 2022).
    return subprocess.run(['phenix.process_predicted_model', model_path], capture_output=True, text=True).stdout

Report an NMR Ensemble Spread (never average coordinates)

Goal: Turn a multi-model NMR deposit into a per-residue uncertainty map instead of over-claiming precision from model 1.

Approach: Superpose all models on a reference and report per-residue Ca RMSD across the ensemble; wide spread marks flexible or under-restrained regions.

import numpy as np
from Bio.PDB import PDBParser

def ensemble_ca_spread(structure, chain_id):
    coords = []
    for model in structure:
        coords.append(np.array([res['CA'].get_coord() for res in model[chain_id] if 'CA' in res]))
    stack = np.stack(coords)  # (n_models, n_residues, 3); assumes consistent residue set across models
    # Per-residue spread = mean distance of each model's Ca from the ensemble mean position.
    return np.linalg.norm(stack - stack.mean(axis=0), axis=2).mean(axis=0)

Secondary-structure validation (DSSP, Bio.PDB.DSSP(model, path, dssp='mkdssp')) needs backbone geometry to place the amide H and computes its own H-bond energy; it processes only the first model and its output drifts across the dssp->mkdssp v2->v4 rewrites, so name the version (Kabsch & Sander 1983 Biopolymers 22:2577-2637). See geometric-analysis for dihedral and DSSP mechanics.