"""Recompute pDockQ and pDockQ2 for two AlphaFold Database dimer models. Reproduces the predicted scores in content/guides/dockq.mdx: python analyze.py Inputs in this directory, downloaded unchanged on October 11, 2026 from https://alphafold.ebi.ac.uk/files/: - AF-0000000211046780-*: hemoglobin subunit alpha (chain A) with subunit beta (chain B) - AF-0000000210768563-*: CDK2 (chain A) with CKS1B (chain B) Each model has a structure file (*-model_v1.cif), a PAE file (*-predicted_aligned_error_v1.json) and a pLDDT file (*-confidence_v1.json). Both scores follow the reference implementation in ipsae.py (https://github.com/DunbrackLab/IPSAE): residues are represented by their CB atom (CA for glycine), interface contacts are interchain pairs within 8 A, pDockQ uses Bryant et al. (2022) and pDockQ2 uses Zhu et al. (2023). The global metrics AlphaFold DB stores in each structure file are read back for comparison. DockQ itself needs the experimental structures and DockQ 2.1.3; the commands are in README.md and the outputs in dockq_*.json. """ import json import math from pathlib import Path from statistics import mean HERE = Path(__file__).resolve().parent CONTACT_CUTOFF = 8.0 PAE_D0 = 10.0 KEPT_DOCKQ_FIELDS = ("DockQ", "iRMSD", "LRMSD", "fnat", "nat_correct", "nat_total", "fnonnat", "F1") MODELS = { "AF-0000000211046780": "Hemoglobin alpha / beta", "AF-0000000210768563": "CDK2 / CKS1B", } def read_cif_loop(lines, prefix): """Return the rows of one mmCIF loop as dicts keyed by column name.""" columns, rows = [], [] in_loop = False for line in lines: if line.startswith(prefix + "."): columns.append(line.strip()[len(prefix) + 1 :]) in_loop = True continue if in_loop: if line.startswith("#") or line.startswith("loop_") or line.startswith("_"): if rows: break continue values = line.split() if len(values) == len(columns): rows.append(dict(zip(columns, values))) return rows def representative_atoms(lines): """One coordinate per residue: CB, or CA for glycine, in file order.""" residues = [] for atom in read_cif_loop(lines, "_atom_site"): if atom["group_PDB"] != "ATOM": continue name, residue = atom["label_atom_id"], atom["label_comp_id"] if name == "CB" or (residue == "GLY" and name == "CA"): residues.append( { "chain": atom["auth_asym_id"], "coords": tuple(float(atom[k]) for k in ("Cartn_x", "Cartn_y", "Cartn_z")), } ) return residues def stored_metrics(lines): """Global metrics that AlphaFold DB wrote into the structure file.""" names = {row["id"]: row["name"].strip("'") for row in read_cif_loop(lines, "_ma_qa_metric")} values = {} for row in read_cif_loop(lines, "_ma_qa_metric_global"): values[names[row["metric_id"]]] = float(row["metric_value"]) return values def interface(residues, chain1, chain2): pairs = [] for i, a in enumerate(residues): if a["chain"] != chain1: continue for j, b in enumerate(residues): if b["chain"] == chain2 and math.dist(a["coords"], b["coords"]) <= CONTACT_CUTOFF: pairs.append((i, j)) return pairs def pdockq(mean_plddt, n_contacts): x = mean_plddt * math.log10(n_contacts) return 0.724 / (1 + math.exp(-0.052 * (x - 152.611))) + 0.018 def pdockq2(mean_plddt, mean_pae_term): x = mean_plddt * mean_pae_term return 1.31 / (1 + math.exp(-0.075 * (x - 84.733))) + 0.005 def analyze(model_id): lines = (HERE / f"{model_id}-model_v1.cif").read_text().splitlines() residues = representative_atoms(lines) pae = json.loads((HERE / f"{model_id}-predicted_aligned_error_v1.json").read_text())[0][ "predicted_aligned_error" ] plddt = json.loads((HERE / f"{model_id}-confidence_v1.json").read_text())["confidenceScore"] assert len(residues) == len(pae) == len(plddt) pairs = interface(residues, "A", "B") interface_residues = sorted({i for i, _ in pairs} | {j for _, j in pairs}) interface_plddt = mean(plddt[k] for k in interface_residues) directions = {} for label, oriented in (("A->B", pairs), ("B->A", [(j, i) for i, j in pairs])): # Rows of the PAE matrix are aligned residues, columns are scored residues. pae_term = mean(1 / (1 + (pae[i][j] / PAE_D0) ** 2) for i, j in oriented) directions[label] = { "mean_interface_pae": round(mean(pae[i][j] for i, j in oriented), 2), "mean_pae_term": round(pae_term, 4), "pDockQ2": round(pdockq2(interface_plddt, pae_term), 4), } stored = stored_metrics(lines) return { "model": model_id, "complex": MODELS[model_id], "chain_lengths": {c: sum(r["chain"] == c for r in residues) for c in ("A", "B")}, "interface_contacts": len(pairs), "interface_residues": len(interface_residues), "mean_interface_plddt": round(interface_plddt, 2), "pDockQ": round(pdockq(interface_plddt, len(pairs)), 4), "pDockQ2": directions, "pDockQ2_max": max(d["pDockQ2"] for d in directions.values()), "alphafold_db": { "ipTM": stored.get("ipTM"), "ipSAE_AB": stored.get("ipsae_AB"), "ipSAE_BA": stored.get("ipsae_BA"), "pDockQ2_AB": stored.get("pDockQ2_AB"), "pDockQ2_BA": stored.get("pDockQ2_BA"), }, } def main(): summary = {"models": [analyze(model_id) for model_id in MODELS]} for name in ( "dockq_hemoglobin_2HHB.json", "dockq_hemoglobin_alpha1beta2_2HHB.json", "dockq_cdk2_cks1_1BUH.json", ): native = json.loads((HERE / name).read_text()) summary.setdefault("dockq", {})[name] = { "GlobalDockQ": round(native["GlobalDockQ"], 3), "interfaces": { key: {k: round(v, 3) for k, v in result.items() if k in KEPT_DOCKQ_FIELDS} for key, result in native["best_result"].items() }, } text = json.dumps(summary, indent=2) + "\n" (HERE / "summary.json").write_text(text) print(text) if __name__ == "__main__": main()