第3章：AlphaFold2中的PAE：输出格式与解读方法

PAE这玩意儿，说白了就是AlphaFold2告诉你：它对两个残基之间的相对位置，到底有多没底。很多新手拿到PDB文件就以为万事大吉了，实际上真正的宝藏藏在PAE文件里。我个人习惯是把PAE作为结构可信度的核心判据，比pLDDT还要先看。

3.1 PAE的原始输出长什么样？

AlphaFold2跑完后，会在输出目录里生成一堆文件。PAE并不是一个单独的文件名，而是藏在其中的数据。具体来说，AF2会输出一个result_model_*.pkl文件，PAE就打包在里面。

我见过不少结构生物学的朋友，跑去输出目录里翻来翻去找PAE，结果发现根本找不到单独的PAE文件，甚至开始怀疑模型没跑对。实际上它就在pickle文件里，需要你用Python提取出来。

3.2 如何把PAE从pickle里拽出来？

我们来看看实际操作的代码。假设你跑完一个目标蛋白，输出目录在./output/下。

import pickle
import numpy as np

# 加载pickle文件
with open('./output/result_model_1_ptm.pkl', 'rb') as f:
    data = pickle.load(f)

# 查看所有key，找找PAE在哪
print(data.keys())

# 输出中你会看到类似：
# dict_keys(['plddt', 'pae', 'ptm', 'iptm', 'pae_args'])

# 提取PAE矩阵
pae_matrix = data['pae']  # 形状: [L, L]
print(pae_matrix.shape)
# 输出例子: (256, 256)

嗯，这里要注意一下。PAE矩阵的形状是[L, L]，L是蛋白的序列长度。比如一个256个残基的蛋白，PAE矩阵就是256×256。

3.3 PAE矩阵到底在说什么？

PAE矩阵中，位置(i, j)的值代表：当参考残基i对齐时，残基j的预测位置的预期误差，单位是埃（Å）。你想想看，这个值越大，就说明AlphaFold2认为残基j相对于残基i的位置越不确定。

搞明白了这个，再来看一个完整的解读流程。我习惯把PAE矩阵可视化出来，用热图的形式看。

import matplotlib.pyplot as plt
import seaborn as sns

plt.figure(figsize=(8, 6))
sns.heatmap(pae_matrix, cmap='viridis_r', vmin=0, vmax=31.75)
plt.xlabel('残基索引 (j)')
plt.ylabel('残基索引 (i)')
plt.title('PAE矩阵热图')
plt.show()

3.4 从PAE矩阵能看出什么门道？

我总结了几条实用经验，帮你快速读图。别一上来就盯着数字看，要先看颜色的全局分布。

第一种情况：结构域内部稳定。如果矩阵中某个对角块的区域都是深色（低PAE值），说明这个区域内部残基之间的相对位置预测得很准。这通常就是一个结构域。

第二种情况：结构域之间松散。不同对角块之间的区域如果呈现亮色（高PAE值），说明这些区域之间的相对位置不确定。典型例子：两个结构域之间是柔性连接，或者蛋白本身就是多结构域的。AlphaFold2对于它们之间的摆位没有把握。

第三种情况：全局高PAE。整张图都是亮的，几乎没有深色块。这说明模型对整体结构都没谱——要么是天然无序蛋白（IDP），要么是预测失败了。这时候pLDDT应该也很低。

3.5 PAE-pairwise图：更精细的解读工具

除了整个矩阵的热图，我还经常用PAE-pairwise图。说白了，就是只看某个残基对所有其他残基的PAE曲线。

举个例子，你想知道第50号残基对整条链的误差估计，可以画这样一条曲线：

import matplotlib.pyplot as plt

residue_idx = 49  # Python下标从0开始，所以第50号残基是下标49
pae_row = pae_matrix[residue_idx, :]

plt.figure(figsize=(10, 4))
plt.plot(range(pae_matrix.shape[0]), pae_row, linewidth=1.5)
plt.xlabel('目标残基索引')
plt.ylabel('相对于残基50的PAE (Å)')
plt.title('第50号残基对每条残基的预期位置误差')
plt.axhline(y=15, color='red', linestyle='--', label='15 Å 警告线')
plt.legend()
plt.show()

3.6 我在项目中踩过的坑

我曾参与一个多结构域蛋白的项目，膜蛋白加胞外域。初始预测结果看着pLDDT还挺正常，膜区70多，胞外域80多。我就信心满满地拿去做分子对接。结果怎么也对接不上。

回头仔细看PAE矩阵，发现胞外域内部的PAE是深色，膜区内部的PAE也是深色，但这两个域之间的PAE全部超过30 Å。说白了，AlphaFold2把这两个域各自的位置预测得不错，但它不知道它们之间怎么摆。整个模型里，两个域的相对取向可能是随意摆放的。

从此我养成了一个习惯：先看PAE矩阵的域间区，再决定是否用全链做下游计算。

3.7 多聚体与PAE：一个容易被忽视的细节

如果你跑的是多聚体（AlphaFold-Multimer），PAE矩阵的形状就不是[L, L]了，而是[L_total, L_total]，L_total是所有链的总长度。位置(i, j)依然表示残基i和残基j之间的预期误差，无论它们是否在同一条链上。

这就有意思了。链间的PAE区域如果亮，说明这两个链之间的接触界面预测得不可靠。如果链间的PAE区域呈现出深色的对角块，说明模型对链间界面很有把握。我在做抗体-抗原复合物预测时，全靠PAE矩阵的链间区域来判断结合界面的可信度。

# 从AlphaFold-Multimer的pickle中提取链间PAE
# 假设chain_A长100个残基，chain_B长80个残基
chain_A_len = 100
chain_B_len = 80
total_len = chain_A_len + chain_B_len

# 链间区域在矩阵中的位置
inter_pae = pae_matrix[chain_A_len:total_len, 0:chain_A_len]
inter_pae_avg = np.mean(inter_pae)
print(f"链间平均PAE: {inter_pae_avg:.2f} Å")

3.8 核心流程图：PAE解读的决策路径

我把它整理成一个简单的流程图。如果你拿到一个AlphaFold2输出的PAE结果，按照这个路径走就行。

3.9 几个常用的PAE处理函数

最后我分享几个自己常用的子函数，你可以直接复制进脚本里用。这些函数帮我省了不少时间。

def calculate_domain_pae(pae_matrix, domain_indices):
    """
    计算指定结构域内部的平均PAE。
    参数:
        pae_matrix: PAE矩阵，形状[L, L]
        domain_indices: 结构域残基索引的列表
    返回:
        该结构域内部平均PAE
    """
    sub_matrix = pae_matrix[np.ix_(domain_indices, domain_indices)]
    return np.mean(sub_matrix)

def calculate_inter_domain_pae(pae_matrix, domain1, domain2):
    """
    计算两个结构域之间的平均PAE。
    参数:
        pae_matrix: PAE矩阵，形状[L, L]
        domain1, domain2: 结构域残基索引列表
    返回:
        两个结构域之间的平均PAE（值越大，相对位置越不可靠）
    """
    sub_matrix = pae_matrix[np.ix_(domain1, domain2)]
    return np.mean(sub_matrix)

def pae_confidence_region(pae_matrix, threshold=15.0):
    """
    返回PAE低于阈值的残基对索引。
    用于判断哪些区域是模型比较有把握的。
    """
    confident_pairs = np.where(pae_matrix < threshold)
    return list(zip(confident_pairs[0], confident_pairs[1]))

3.10 别忘了AlphaFold-Multimer的PAE变体

跑AlphaFold-Multimer时，输出里可能同时有pae和pae_plddt两个字段。前者是纯粹的预期位置误差矩阵，后者是对PAE按pLDDT做了加权平均后得到的标量值。我个人习惯直接用原始的pae，因为它能提供更细粒度的空间信息。

# 从AlphaFold-Multimer输出中检查PAE相关字段
data = pickle.load(open('result_multimer_v2.pkl', 'rb'))
print('pae' in data)  # True
print('pae_plddt' in data)  # True，注意这是标量，不是矩阵
print('ptm' in data)  # True，pTM评分
print('iptm' in data)  # True，界面pTM