"""Reproduce the worked example in the ProteinIQ structure prediction confidence guide. Reads the public ProteinIQ AlphaFold2 example job for the GCN4 leucine zipper homodimer (https://proteiniq.io/app/alphafold-2?jobId=leiy7ew8yacrd7k2srsqqkdk), summarizes its pLDDT and PAE, and collects the ipSAE, DockQ, US-align and MolProbity results stored next to this script into summary.json. Python 3 standard library only. Run from this directory: python -I analyze.py """ import json import math import re from pathlib import Path from statistics import mean HERE = Path(__file__).resolve().parent CHAIN_LENGTH = 33 # residues per chain in the model; chain A is 0-32, chain B is 33-65 def chain_of(index): return "A" if index < CHAIN_LENGTH else "B" def read_cb_coordinates(path): """One CB atom per residue (CA for glycine), in file order.""" coords = [] for line in path.read_text().splitlines(): if not line.startswith("ATOM"): continue atom = line[12:16].strip() residue = line[17:20] if atom == "CB" or (atom == "CA" and residue == "GLY"): coords.append((float(line[30:38]), float(line[38:46]), float(line[46:54]))) return coords def block_stats(pae, rows, cols): values = [pae[i][j] for i in rows for j in cols] return { "mean_pae": round(mean(values), 2), "fraction_below_5": round(sum(v < 5 for v in values) / len(values), 3), } def main(): scores = json.loads((HERE / "gcn4_af2_scores.json").read_text()) plddt = scores["plddt"] pae = scores["pae"] n = len(plddt) a = range(0, CHAIN_LENGTH) b = range(CHAIN_LENGTH, n) plddt_summary = { "mean": round(mean(plddt), 1), "chain_A_mean": round(mean(plddt[i] for i in a), 1), "chain_B_mean": round(mean(plddt[i] for i in b), 1), "residues_above_90": sum(v > 90 for v in plddt), "residues": n, "min": min(plddt), "chain_A_first_two": plddt[0:2], "chain_A_last_three": plddt[30:33], "chain_A_residues_3_to_30_mean": round(mean(plddt[2:30]), 1), } # PAE blocks: rows are the aligned residue, columns the scored residue. blocks = { "A_to_A": block_stats(pae, a, a), "B_to_B": block_stats(pae, b, b), "A_to_B": block_stats(pae, a, b), "B_to_A": block_stats(pae, b, a), } # Alignment frames: mean PAE of each row (aligned residue) over the whole dimer. # Core residues 3-30 of each chain versus the frayed ends (residue 1 and 31-33). row_mean = [mean(row) for row in pae] core = [i for c in (0, CHAIN_LENGTH) for i in range(c + 2, c + 30)] ends = [i for c in (0, CHAIN_LENGTH) for i in (c, c + 30, c + 31, c + 32)] core_set = set(core) near = [pae[i][j] for i in a for j in a if 0 < abs(i - j) <= 4] far_core = [pae[i][j] for i in core_set & set(a) for j in core_set & set(a) if abs(i - j) >= 20] along_chain = { "within_chain_A_separation_1_to_4_mean": round(mean(near), 2), "core_rows_mean_pae": round(mean(row_mean[i] for i in core), 2), "end_rows_mean_pae_range": [round(min(row_mean[i] for i in ends), 1), round(max(row_mean[i] for i in ends), 1)], "chain_A_row_means": [round(row_mean[i], 1) for i in a], "chain_A_core_pairs_20_plus_apart_mean": round(mean(far_core), 2), } # Interchain contacts: CB (CA for Gly) within 8 Angstrom, as in pDockQ. cb = read_cb_coordinates(HERE / "gcn4_af2_model.pdb") contacts = [ (i, j) for i in a for j in b if math.dist(cb[i], cb[j]) <= 8.0 ] interface = sorted({i for i, _ in contacts} | {j for _, j in contacts}) contact_pae = [(pae[i][j] + pae[j][i]) / 2 for i, j in contacts] noncontact_pae = [ (pae[i][j] + pae[j][i]) / 2 for i in a for j in b if (i, j) not in set(contacts) ] interface_summary = { "contacts_cb_8A": len(contacts), "interface_residues": len(interface), "interface_plddt_mean": round(mean(plddt[i] for i in interface), 1), "contact_pair_pae_mean": round(mean(contact_pae), 2), "noncontact_interchain_pae_mean": round(mean(noncontact_pae), 2), } # ipSAE report (ipsae.py at commit 6174cf9, PAE cutoff 10, distance cutoff 15). lines = (HERE / "ipsae_10_15.txt").read_text().split("\n") header = lines[1].split() if not lines[0].strip() else lines[0].split() ipsae_rows = [dict(zip(header, line.split())) for line in lines if line.startswith(("A ", "B "))] dockq = json.loads((HERE / "dockq_2ZTA.json").read_text())["best_result"]["AB"] def usalign(name): text = (HERE / name).read_text() rmsd = re.search(r"Aligned length=\s*(\d+), RMSD=\s*([\d.]+)", text) tm = re.findall(r"TM-score= ([\d.]+) \(normalized by length of Structure_(\d)", text) return { "aligned_length": int(rmsd.group(1)), "ca_rmsd": float(rmsd.group(2)), "tm_score_by_model": float(tm[0][0]), "tm_score_by_reference": float(tm[1][0]), } molprobity_text = (HERE / "molprobity_model.txt").read_text() summary_block = molprobity_text.split("= Summary =")[1] def mp(label): return float(re.search(label + r"\s*=\s*([\d.]+)", summary_block).group(1)) summary = { "job": "https://proteiniq.io/app/alphafold-2?jobId=leiy7ew8yacrd7k2srsqqkdk", "model_reported": { "ptm": scores["ptm"], "iptm": scores["iptm"], "actifptm": scores["actifptm"], "per_chain_ptm": scores["per_chain_ptm"], }, "plddt": plddt_summary, "pae_blocks": blocks, "pae_along_chain": along_chain, "interface": interface_summary, "ipsae": ipsae_rows, "dockq_vs_2ZTA": {k: dockq[k] for k in ("DockQ", "fnat", "nat_correct", "nat_total", "iRMSD", "LRMSD")}, "usalign_complex_vs_2ZTA": usalign("usalign_complex.txt"), "usalign_chainA_vs_2ZTA_chainA": usalign("usalign_chainA.txt"), "molprobity_model": { "ramachandran_outliers_pct": mp("Ramachandran outliers"), "ramachandran_favored_pct": mp("favored"), "rotamer_outliers_pct": mp("Rotamer outliers"), "clashscore": mp("Clashscore"), "molprobity_score": mp("MolProbity score"), }, } (HERE / "summary.json").write_text(json.dumps(summary, indent=2) + "\n") print(json.dumps(summary, indent=2)) if __name__ == "__main__": main()