AlphaFold 结果怎么看：pLDDT、PAE、pTM 与 ipTM 的含义和判断阈值

#!/usr/bin/env python3
"""读取 AlphaFold 3 / AlphaFold Server / ColabFold 的置信度 JSON，输出链对指标并画 PAE 热图。

AlphaFold 3 本地版:  python af_pae.py job_confidences.json --summary job_summary_confidences.json
AlphaFold Server:    python af_pae.py fold_job_full_data_0.json --summary fold_job_summary_confidences_0.json
ColabFold:           python af_pae.py job_scores_rank_001_xxx.json --lengths 110,89
"""
import argparse
import json

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np


def segments(ids):
    """把逐 token 的链 ID 切成 [(链, 起点, 终点)]，终点不含。"""
    segs, start = [], 0
    for i in range(1, len(ids) + 1):
        if i == len(ids) or ids[i] != ids[start]:
            segs.append((ids[start], start, i))
            start = i
    return segs


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("conf", help="含 pae 矩阵的 JSON")
    ap.add_argument("--summary", help="AF3 的 summary_confidences JSON")
    ap.add_argument("--lengths", help="ColabFold 用：各链长度，逗号分隔")
    ap.add_argument("--cutoff", type=float, default=10.0, help="统计界面用的 PAE 阈值 (Å)")
    ap.add_argument("--out", default="pae.png")
    args = ap.parse_args()

    with open(args.conf) as f:
        d = json.load(f)
    pae = np.asarray(d.get("pae", d.get("predicted_aligned_error")), dtype=float)
    n = pae.shape[0]

    if "token_chain_ids" in d:  # AF3 本地版与 AlphaFold Server
        chain_of_token = d["token_chain_ids"]
    elif args.lengths:  # ColabFold：scores JSON 不含链 ID
        lens = [int(x) for x in args.lengths.split(",")]
        if sum(lens) != n:
            raise SystemExit(f"--lengths 之和 {sum(lens)} 与 PAE 维度 {n} 不一致")
        chain_of_token = [chr(65 + k) for k, m in enumerate(lens) for _ in range(m)]
    else:
        chain_of_token = ["A"] * n
    segs = segments(chain_of_token)

    # 每条链的平均 pLDDT
    if "atom_plddts" in d and "atom_chain_ids" in d:  # AF3：逐原子
        pl = np.asarray(d["atom_plddts"])
        ids = np.asarray(d["atom_chain_ids"])
        for c in dict.fromkeys(d["atom_chain_ids"]):
            print(f"chain {c}: mean pLDDT (atoms) = {pl[ids == c].mean():.1f}")
    elif "plddt" in d:  # ColabFold：逐残基
        pl = np.asarray(d["plddt"])
        for c, s, e in segs:
            print(f"chain {c}: mean pLDDT (residues) = {pl[s:e].mean():.1f}")

    # 链对 PAE 统计：行 = 对齐所用的链，列 = 被评估的链
    print(f"\naligned scored  min_PAE  mean_PAE  pairs<{args.cutoff:g}A")
    for ci, si, ei in segs:
        for cj, sj, ej in segs:
            if ci == cj:
                continue
            b = pae[si:ei, sj:ej]
            print(f"{ci:>7} {cj:>6} {b.min():8.2f} {b.mean():9.2f} {int((b < args.cutoff).sum()):10d}")

    # 汇总指标
    for key in ("ptm", "iptm", "ipsae", "pdockq", "pdockq2", "actifptm"):
        if key in d:
            print(f"{key}: {d[key]}")
    if args.summary:
        with open(args.summary) as f:
            s = json.load(f)
        for key in ("ptm", "iptm", "ranking_score", "fraction_disordered", "has_clash"):
            print(f"{key}: {s.get(key)}")
        # AF3 本地版的 chain_ids 是逐 token 的列表，Server 版没有该字段；去重后得到链顺序
        cids = list(dict.fromkeys(s.get("chain_ids") or [])) or [c for c, _, _ in segs]
        for key in ("chain_pair_iptm", "chain_pair_pae_min"):
            if key in s:
                print(f"\n{key} (rows/cols = {cids})")
                for c, row in zip(cids, s[key]):
                    print(c, " ".join("  nan" if v is None else f"{v:5.2f}" for v in row))

    # PAE 热图
    fig, ax = plt.subplots(figsize=(6, 5), dpi=200)
    im = ax.imshow(pae, cmap="Greens_r", vmin=0, vmax=30)
    for _, s, _ in segs[1:]:
        ax.axhline(s - 0.5, color="black", lw=0.8)
        ax.axvline(s - 0.5, color="black", lw=0.8)
    ax.set_xlabel("Scored token")
    ax.set_ylabel("Aligned token")
    fig.colorbar(im, ax=ax, label="Expected position error (A)")
    fig.tight_layout()
    fig.savefig(args.out)
    print(f"\nsaved {args.out}")


if __name__ == "__main__":
    main()

# 只需要 numpy 和 matplotlib
pip install numpy matplotlib

# AlphaFold 3 本地版（排名第一的模型）
python af_pae.py myjob/myjob_confidences.json --summary myjob/myjob_summary_confidences.json --out myjob_pae.png

# AlphaFold Server（解压后，N=0 为排名第一）
python af_pae.py fold_myjob_full_data_0.json --summary fold_myjob_summary_confidences_0.json --out pae_0.png

# ColabFold：scores JSON 不含链信息，用 --lengths 按顺序给出各链长度
python af_pae.py out/barnase_barstar_scores_rank_001_*.json --lengths 110,89 --out pae_rank1.png

load fold_myjob_model_0.cif, m
hide everything, m
show cartoon, m
# B-factor 列即 pLDDT。按 AlphaFold DB 的四段配色：<=50 橙、50-70 黄、70-90 浅蓝、>90 深蓝
color 0xFF7D45, m
color 0xFFDB13, m and b > 50
color 0x65CBF3, m and b > 70
color 0x0053D6, m and b > 90
# 连续配色（可选）：50 以下全红，90 以上全蓝
# spectrum b, red_yellow_green_cyan_blue, m, minimum=50, maximum=90
# 界面残基：与对方链任意原子距离 5 A 以内
select ifA, byres ((m and chain A) within 5 of (m and chain B))
select ifB, byres ((m and chain B) within 5 of (m and chain A))
select iface, ifA or ifB
show sticks, iface and not name N+C+O
# 打印界面残基的 pLDDT（取 CA 原子）
iterate iface and name CA, print(chain, resi, resn, round(b, 1))
# 另存为 PDB
save model_0.pdb, m

open fold_myjob_model_0.cif
# ChimeraX 默认的 pLDDT 连续配色
color bfactor #1 palette alphafold
# 读入 PAE：AF3 的 full_data / confidences JSON、ColabFold 的 scores JSON 都可直接读
alphafold pae #1 file fold_myjob_full_data_0.json plot true
# 把 A、B 链之间距离 <=4 A 且 PAE <=5 A 的残基对画成伪键
alphafold contacts #1/A toAtoms #1/B distance 4 maxPae 5
# 按 PAE 聚类出刚性结构域并着色（默认 connectMaxPae 5 A）
alphafold pae #1 colorDomains true

pip install gemmi
# --shorten 把链名缩到 1-2 个字符；配体用了 5 字符单体名时再加 --shorten-tlc
gemmi convert --shorten fold_myjob_model_0.cif model_0.pdb

cat > barnase_barstar.fasta <<'EOF'
>barnase_barstar
AQVINTFDGVADYLQTYHKLPDNYITKSEAQALGWVASKGNLADVAPGKSIGGDIFSNREGKLPGKSGRTWREADINYTSGFRNSDRILYSSDWLIYKTTDHYQTFTKIR:KKAVINGEQIRSISDLHQTLKKELALPEYYGENLDALWDALTGWVEYPLVLEWRQFEQSKQLTENGAESVLQVFREAKAEGADITIILS
EOF

# 5 个模型 x 5 个 seed = 25 个预测；复合物默认按 80*ipTM + 20*pTM 排序
colabfold_batch --num-models 5 --num-seeds 5 --calc-extra-ptm \
  --amber --num-relax 1 --use-gpu-relax \
  barnase_barstar.fasta out/

# 汇总 25 个模型的 ipTM、pTM、ipSAE、pDockQ2（1.6.3 起 scores JSON 自带后两项）
# ipsae、pdockq2 按链对存成字典（A-B、B-A），这里取较大值，与 log.txt 一致
python - <<'PY'
import glob, json
for f in sorted(glob.glob("out/barnase_barstar_scores_rank_*.json")):
    s = json.load(open(f))
    rank = f.split("_scores_")[1].split("_alphafold2")[0]
    mx = lambda v: max(v.values()) if isinstance(v, dict) and v else v
    print(rank, s.get("iptm"), s.get("ptm"), mx(s.get("ipsae")), mx(s.get("pdockq2")))
PY

# 与实验结构 1BRS 对比：模型 A、B 链对应 1BRS 的 A（barnase）、D（barstar）链
wget https://files.rcsb.org/download/1BRS.pdb
pip install DockQ
for f in out/barnase_barstar_unrelaxed_rank_*.pdb; do
  DockQ "$f" 1BRS.pdb --mapping AB:AD --short | tail -1
done