4. PAE矩阵可视化：热力图绘制工具与颜色解读

好，咱们进入第四章。前面聊了PAE是什么，怎么算的，现在该上手看看它长什么样了。

我个人觉得，PAE矩阵这东西，光看数字没用。你必须把它变成一张图，一眼扫过去，哪些区域粘得紧，哪些区域晃得厉害，心里就有数了。说白了，可视化是把抽象的距离预测误差，变成你能直接“看见”的蛋白质行为地图。

这一章咱们就解决两个问题：用什么工具画？画出来怎么读？

4.1 工具链：从py3Dmol到matplotlib

我最早接触AlphaFold2的时候，首选的可视化工具是py3Dmol。为什么？因为它原生支持AF2的输出格式，你下载下来的result_model_*.pkl文件里，直接就有pae这个字段。

先给个快速入门的例子，看看py3Dmol怎么把PAE叠到结构上：

import py3Dmol
import numpy as np

# 假装你已经从pkl里读出了pae矩阵和pdb文件
pae_matrix = np.load('pae_matrix.npy')  # shape (L, L)
pdb_file = 'predicted_structure.pdb'

view = py3Dmol.view(js='https://3dmol.org/build/3Dmol.js')
view.addModel(open(pdb_file).read(), 'pdb')

# 关键一步：用PAE值给残基上色
# 这里我用的是py3Dmol的setStyle接口，把pae映射到色标上
view.setStyle({'cartoon': {'colorscheme': {'prop': 'b', 
                                            'map': 'RdYlBu_r', 
                                            'min': 0, 
                                            'max': 30}}})
view.zoomTo()

等一下，上面这段其实不太严谨。py3Dmol原生的颜色映射不支持直接把PAE矩阵当属性喂进去，通常的做法是用setStyle里的colorfunc，或者先在pdb文件里把PAE存成B-factor字段。我平时更习惯第二种方法：

import pandas as pd
from Bio.PDB import PDBParser, PDBIO

# 把pae矩阵里的对角数据提取出来，比如残基i的平均PAE
mean_pae = pae_matrix.mean(axis=1)

# 修改pdb的B-factor列
parser = PDBParser(QUIET=True)
structure = parser.get_structure('protein', pdb_file)
for i, chain in enumerate(structure.get_chains()):
    for j, residue in enumerate(chain.get_residues()):
        for atom in residue:
            atom.set_bfactor(mean_pae[j])

io = PDBIO()
io.set_structure(structure)
io.save('pae_colored.pdb')

# 再在py3Dmol里加载
view = py3Dmol.view()
view.addModel(open('pae_colored.pdb').read(), 'pdb')
view.setStyle({'cartoon': {'colorscheme': {'prop': 'b', 'gradient': 'rdylbu', 'min': 0, 'max': 30}}})
view.zoomTo()
view.show()

这样渲染出来，深蓝色的区域就是PAE低、预测很稳的地方；红色区域是PAE高、预测不确定的地方。你想想看，一个蛋白的结构域内部通常都是蓝色，而结构域之间会有一道“红墙”，这就是domain boundary所在的位置。

但py3Dmol有个问题——它只适用于交互查看，你不能把它直接放到论文或报告里作为静态图。再者，想要观察PAE矩阵的整体模式（比如模块化分块），在3D结构上反而不如在2D热力图里直观。

实用建议：我通常把py3Dmol用在探索阶段。第一天拿到模型，先快速在3D结构上扫一遍，标记出高PAE区域。等到要分析domain之间的具体交互关系了，再切换到matplotlib画热力图。

4.2 用matplotlib绘制PAE热力图

再说回2D热力图。其实matplotlib自己就能画，而且自由度更高。你需要控制的东西不多：坐标轴的残基索引、颜色映射、还有最重要的——怎么标出domain boundaries。

先来一个最基础版本：

import matplotlib.pyplot as plt
import numpy as np

def plot_pae_heatmap(pae_matrix, title='PAE Heatmap'):
    plt.figure(figsize=(8, 6))
    plt.imshow(pae_matrix, cmap='viridis_r', vmin=0, vmax=30)
    plt.colorbar(label='Predicted Aligned Error (Å)')
    plt.xlabel('Residue index')
    plt.ylabel('Residue index')
    plt.title(title)
    plt.tight_layout()
    plt.show()

viridis_r这个颜色映射我经常用——深色表示低误差，浅色表示高误差。推荐几个适合PAE的色标：viridis_r、magma_r、RdYlBu_r。其实RdYlBu_r很直观，红色代表高PAE（不可信），蓝色代表低PAE（可信），和直觉一致。

但光会画不行，你得会“读”。我经常在热力图上加辅助线，把domain boundaries标出来：

# 假设domain_boundaries是已知的切割点列表，比如[50, 120, 200]
domain_boundaries = [50, 120, 200]
for b in domain_boundaries:
    plt.axvline(b, color='white', linestyle='--', linewidth=0.8)
    plt.axhline(b, color='white', linestyle='--', linewidth=0.8)

为什么要加这个？因为PAE热力图有显著的“块状结构”。如果蛋白只有两个域，你会看到一个2×2的方块：对角线两个方块是深色的（域内预测准），副对角线的两个方块是浅色的（域间预测不准）。

我曾经踩过一个坑：有一次我看到的PAE热力图副对角线区域特别亮，但不算特别高（只有15Å左右），我就以为这两个域之间有微弱相互作用。后来用实验验证才发现，根本不接触。那为什么PAE不高？因为AlphaFold2的「平均」把很多不同的构象模式抹平了，低的不一定就是对的——你可能看到的是多种构象的“平均”结果。这个坑后面会细讲。

4.3 解读颜色模式：从色标到物理含义

现在进入核心问题：PAE热力图的颜色到底代表什么？

直接给结论：颜色代表“如果我把残基i和残基j对齐到预测结构上，它们真实位置的预期偏差是多少（单位：Å）”。

换句话说，你看到矩阵中一个像素点的颜色深（比如深紫色或深蓝色），表示模型对这两个残基的相对位置很有信心。如果颜色浅（黄色、红色），表明模型拿不准这俩残基到底离多远。

我来总结一下AlphaFold2和RoseTTAFold在颜色含义上的细微差别：

特征 AlphaFold2 PAE RoseTTAFold pLDDT / PAE
单位 Å（埃） pLDDT无单位（0~100）；PAE单位为Å
取值范围 0 ~ 31.75（封顶） 不封顶，但一般也落在0~30之间
低值含义 DOM对的内部构象高度确定 同左
高值含义 不同DOM之间存在不确定的距离 同左
对角线特征 几乎全为低值（自身对齐完美） 对角线区域也可能出现噪声（RoseTTAFold稍弱）
常见色标 蓝-白-红（BlueWhiteRed） 绿-黄-红（或viridis_r）

你看到这里可能会问：PAE热力图上对角线像素的颜色和边上一样深吗？答案是：不一定。AlphaFold2会把对角线上的PAE强制设为0，所以你看到的对角线永远是深色；但RoseTTAFold不会强制，偶尔对角线上会出现一些浅色条纹。这说明RoseTTAFold在局部残基的相对位置预测上，偶尔会出现自相矛盾的情况。

核心解读法则：看PAE热力图，永远盯着副对角线区域。如果副对角线上出现两个明显的浅色方块，那就代表这个蛋白很可能是多域结构，且域间相对位置不确定。如果整个矩阵几乎全是深色（比如<5Å），那就说明这是一个紧凑的单域蛋白，或者是一个刚性非常大的多域复合物。

4.4 实战：一张图读懂domain关系

来一个实际案例。假设我们有一个蛋白，长度是500个残基。你用AlphaFold2跑了一遍，得到如下PAE矩阵。画出来之后，你看到什么？

我给你画一张SVG框架图示意一下：

PAE热力图：典型多域蛋白结构 DOM1 DOM2 DOM3 DOM1 DOM2 DOM3 图例说明 域内：低PAE（高置信度） 域间：高PAE（低置信度）

这张图里，你看得很清楚：每个domain内部（对角线方块）是蓝色的，因为预测准确；但当你跨domain去对齐（比如DOM1对DOM2的副对角线区域），颜色就变红或变成黄色，代表预测误差很大。

这就是典型的“多域蛋白”信号。我记得有一次我帮合作者分析一个激酶，他坚持认为自己的蛋白是单域。我把PAE矩阵甩给他，指着那个显眼的红色交叉说：“你看，这里130到180的残基和280到340的残基之间的PAE达到了20Å以上，这绝对是一个铰链区，两个domain是相对独立的。”

后来他用SAXS验证了，确实是一个双域构象。从那以后，他每次做AF2都会先看PAE热力图。

避坑指南（来自我自己的教训）：

千万不要只看PAE矩阵的‘平均值’或‘中位数’来判断结构质量。你可能遇到一种情况：整体平均PAE只有8Å，看着还行。但当你仔细去看热力图，发现整张图都是均匀的浅蓝色（10-15Å），没有任何明显的深色块。这意味着什么？意味着模型对你的蛋白整个都是模糊的——没有块状结构，说明它根本没预测出任何稳定的domain。

这时候你应该警惕：要么序列太短（<50个残基），要么多序列比对深度不够，要么蛋白本身在生理条件下就是高度动态的。

4.5 自定义颜色映射：让关键信号更突出

前面用了viridis_r和RdYlBu_r，但我个人觉得，在PAE热力图里，最好用的还是自己定义颜色映射。为什么？因为默认色标往往在大约10-15Å区间出现渐变平缓区，你很难肉眼分辨一个点是12Å还是18Å。

我常用的自定义方案：

import matplotlib.colors as mcolors

# 定义一个分段色标：0-10Å蓝色，10-20Å白色渐变到红色，20-30Å深红色
colors = ['#2166ac', '#4393c3', '#92c5de', '#f7f7f7', '#f4a582', '#d6604d', '#b2182b']
nodes = [0.0, 0.15, 0.35, 0.5, 0.7, 0.85, 1.0]
cmap_custom = mcolors.LinearSegmentedColormap.from_list('pae_custom', 
                                                        list(zip(nodes, colors)))

plt.imshow(pae_matrix, cmap=cmap_custom, vmin=0, vmax=30)

这个方案的妙处在哪里？它在10-20Å这个“犹疑区”用了一个白色到淡红色的渐变，视觉敏感度更高。你一眼就能看出哪些区域是稍微有点不确定但还没到完全不可信的程度。我很多次都是在白色区域发现“域间微弱接触”的信号——PAE在12-14Å之间，但这其实反映了模型对那个界面只有中等程度的确定，不是完全随机。

嗯，关于颜色映射我再说一个小技巧：不要太迷信默认的蓝色代表误差小的映射。有些团队喜欢用“冷色低误差，暖色高误差”的方案，但如果你色弱，可能红色和绿色分辨不了。我建议给组里做汇报时，用viridis_r——它的深浅变化对色弱人群相对友好。

← 上一章 📖 返回目录 下一章 →