SunHongxun'S Site
计算生物学 · 药物发现
← 从侧边栏选择分类开始阅读
📖
实验手册
03 — 自由能形貌图（FEL）绘制完整指南
本文档详细说明如何从 MD 轨迹生成数据、通过 gmx sham 构建自由能形貌图（Free Energy Landscape, FEL），并分别绘制 2D 平面图 和 3D 曲面图。
目录
- 1. 原理简介
- 2. 准备输入数据
- 3. 生成 sham.xvg
- 4. 使用 gmx sham 生成 XPM
- 5. 2D 自由能形貌图绘制
- 6. 3D 自由能形貌图绘制
- 7. 方法二：使用 DuIvyTools 快速绘图
- 8. 完整流程速查
1. 原理简介
自由能形貌图（FEL）以 RMSD（均方根偏差） 为 X 轴、Rg（回旋半径） 为 Y 轴，通过 Boltzmann 反演将构象分布概率转化为自由能：
$$G_i = -k_B T \ln \frac{P_i}{P_{\text{max}}}$$
其中 $P_i$ 是第 $i$ 个构象出现的概率，$T$ 是模拟温度。
2. 准备输入数据
2.1 生成 RMSD 数据
# 生成蛋白质 RMSD 时间序列（单位：ns）
gmx rms -s md100.tpr -f md100_nojump.xtc -o rmsd.xvg -tu ns
交互：
1 # 选择骨架（Backbone）作为拟合组
1 # 选择骨架（Backbone）作为计算组
2.2 生成回旋半径数据
# 生成蛋白质 Rg 时间序列
gmx gyrate -s md100.tpr -f md100_nojump.xtc -o gyrate.xvg
交互：
1 # 选择蛋白质
1 # 选择蛋白质
3. 生成 sham.xvg
gmx sham 需要的输入文件是一个 3 列数据：时间、RMSD、Rg。
3.1 检查时间单位是否一致
常见问题：rmsd.xvg 使用 -tu ns 输出时间为纳秒，而 gyrate.xvg 默认输出时间为皮秒(ps)。
# 查看 RMSD 文件的头部
head -50 rmsd.xvg | grep -E "^[0-9]"
# 查看 Rg 文件的头部
head -50 gyrate.xvg | grep -E "^[0-9]"
如果 RMSD 的时间是 0, 1, 2, ...（ns），而 Rg 的时间是 0, 1000, 2000, ...（ps），需要将 Rg 时间除以 1000。
3.2 使用脚本合并数据
使用以下 Python 脚本生成 sham.xvg：
import numpy as np
# 读取数据（跳过注释行）
rmsd_data = np.loadtxt("rmsd.xvg", comments=["@", "#"])
gyrate_data = np.loadtxt("gyrate.xvg", comments=["@", "#"])
time_rmsd = rmsd_data[:, 0] # RMSD 时间（ns）
rmsd = rmsd_data[:, 1] # RMSD 值（nm）
time_gyr = gyrate_data[:, 0] # Rg 时间
gyrate = gyrate_data[:, 1] # Rg 值（nm）
# 检查 Rg 时间是否以 ps 为单位（数值远大于 RMSD 时间）
if time_gyr.mean() > time_rmsd.mean() * 10:
print("检测到 Rg 时间单位为 ps，正在转换为 ns ...")
time_gyr = time_gyr / 1000.0
# 确保时间范围一致
t_min = max(time_rmsd.min(), time_gyr.min())
t_max = min(time_rmsd.max(), time_gyr.max())
# 输出 sham.xvg
with open("sham.xvg", "w") as f:
f.write("# Time(ns) RMSD(nm) Rg(nm)\n")
for t_r, r, t_g, g in zip(time_rmsd, rmsd, time_gyr, gyrate):
if t_min <= t_r <= t_max and abs(t_r - t_g) < 0.01:
f.write(f"{t_r:.3f} {r:.5f} {g:.5f}\n")
print(f"已生成 sham.xvg，包含对齐的时间点数据")
将上述内容保存为 prepare_sham.py，然后执行：
python prepare_sham.py
生成的 sham.xvg 格式为：
# Time(ns) RMSD(nm) Rg(nm)
0.000 0.12345 1.23456
0.500 0.23456 1.34567
...
4. 使用 gmx sham 生成 XPM
4.1 生成标准形貌图（2D/3D 通用）
gmx sham -f sham.xvg -ls FEL_sham.xpm -b 8500 -tsham 310 -nlevels 100
| 参数 | 含义 |
|---|---|
-f sham.xvg |
输入数据文件（3 列：时间、RMSD、Rg） |
-ls FEL_sham.xpm |
输出 XPM 格式自由能形貌图 |
-b 8500 |
从第 8500 ps 开始分析（跳过前段平衡） |
-tsham 310 |
模拟温度 (K)，用于 Boltzmann 转换 |
-nlevels 100 |
能量等值线层数 |
-bin |
网格 bins 数（可选，默认自动） |
4.2 为 3D 图设置合适的等值线层数
如果使用 xpm2txt.py + 3D 绘图脚本，建议将 nlevels 设为 25（过高的 levels 会使 txt 文件过大），但 2D 图用 100 效果更好：
# 如果是为 3D 图准备，使用较低的 levels
gmx sham -f sham.xvg -ls FEL_3d.xpm -b 8500 -tsham 310 -nlevels 25
# 如果是为 2D 图准备，用较高的 levels 更平滑
gmx sham -f sham.xvg -ls FEL_2d.xpm -b 8500 -tsham 310 -nlevels 100
注意：
-b 8500表示从 8500 ps（8.5 ns）开始分析，通常建议丢弃前 10% 的轨迹作为平衡时间，具体值需根据体系收敛情况调整。
5. 2D 自由能形貌图绘制
方法一：使用 xpm2png.py（推荐）
xpm2png.py 可将 XPM 文件直接导出为高质量的 PNG 图片：
# 基本用法：直接转换
python ../scripts/xpm2png.py -f FEL_sham.xpm
# 常用用法：指定输出文件名并启用插值（推荐）
python ../scripts/xpm2png.py -f FEL_sham.xpm -o FEL_2D.png -ip yes
| 参数 | 说明 |
|---|---|
-f FEL_sham.xpm |
输入 XPM 文件 |
-o FEL_2D.png |
输出 PNG 文件 |
-ip yes |
启用插值（使形貌图更平滑），仅对 Continuous 类型有效 |
-show yes |
显示图片窗口（默认） |
xpm2png.py来自 B 站 UP 主 DuIvyTools，支持： - 自动解析 XPM 中的颜色映射 - 支持 Continuous 和 Discrete 两种类型 - 插值模式使 FEL 图像更平滑美观
方法二：使用 xpm2txt.py + matplotlib
先转换为文本，再用 Python 绘图：
步骤 1：将 XPM 转换为三列文本
python ../scripts/xpm2txt.py -f FEL_sham.xpm -o FEL.txt
步骤 2：用 Python 绘制 2D 填充等高线图
import numpy as np
import matplotlib.pyplot as plt
# 读取数据
data = np.loadtxt("FEL.txt")
x = np.unique(data[:, 0])
y = np.unique(data[:, 1])
X, Y = np.meshgrid(x, y)
Z = data[:, 2].reshape(len(y), len(x))
# 绘制
plt.figure(figsize=(8, 6))
contour = plt.contourf(X, Y, Z, levels=50, cmap='jet')
plt.colorbar(contour, label='Free Energy (kJ/mol)')
plt.xlabel('RMSD (nm)')
plt.ylabel('Rg (nm)')
plt.title('Free Energy Landscape')
plt.tight_layout()
plt.savefig('FEL_2D_matplotlib.png', dpi=300)
plt.show()
6. 3D 自由能形貌图绘制
6.1 流程概览
FEL_sham.xpm → xpm2txt.py → FEL.txt → 3D自由能形貌图.py → FEL_3D.png
6.2 步骤 1：转换 XPM 到 TXT
# 使用 nlevels=25 生成的 XPM 文件
python ../scripts/xpm2txt.py -f FEL_3d.xpm -o FEL_3d.txt
6.3 步骤 2：绘制 3D 曲面图
使用脚本 scripts/3D自由能形貌图.py：
python ../scripts/3D自由能形貌图.py
该脚本执行以下操作：
- 数据预处理：将空白分隔的列转换为逗号分隔（兼容不同格式）
- 数据插值：使用
scipy.interpolate.griddata将离散数据插值到更细的网格（4 倍密度），使曲面更平滑 - 绘图：
- 3D 曲面（
plot_surface，cmap='jet'着色） - 底部等高线投影（
contourf） - 颜色条
6.4 输出文件
gibbs_comma.txt：预处理后的逗号分隔数据gibbs_3d_final.png：最终 3D 形貌图（300 DPI）
6.5 自定义 3D 图形
根据需要修改 scripts/3D自由能形貌图.py 中的以下参数：
# 视角调整
ax.view_init(elev=15, azim=45) # elev: 俯仰角, azim: 方位角
# 颜色映射
cmap='jet' # 可改为 'hot', 'plasma', 'viridis' 等
# 插值密度
xi = np.linspace(xu.min(), xu.max(), len(xu)*4) # 4 倍密度，可改为 2 或 6
7. 方法二：使用 DuIvyTools 快速绘图
DuIvyTools 是一个第三方 GROMACS 分析工具包，提供了更简便的绘图方式：
# 安装
pip install DuIvyTools
# 绘制 3D 自由能形貌图
dit xpm_show -f FEL_sham.xpm -3d -o FEL_3D.pdf
# 绘制 2D 自由能形貌图
dit xpm_show -f FEL_sham.xpm -o FEL_2D.pdf -cmap jet
# 绘制带颜色的 3D 图
dit xpm_show -f FEL_sham.xpm -m 3d -o FEL_3D_color.pdf -cmap jet
优点：命令简洁，一步到位，支持 PDF 矢量格式。
xpm2png.py脚本也是 DuIvyTools 的一部分。
8. 完整流程速查
8.1 2D 形貌图（推荐方案）
# 1. 准备数据
gmx rms -s md100.tpr -f md100_nojump.xtc -o rmsd.xvg -tu ns # 选 1, 1
gmx gyrate -s md100.tpr -f md100_nojump.xtc -o gyrate.xvg # 选 1, 1
# 2. 生成 sham.xvg（需处理时间单位，见第 3 节）
python prepare_sham.py
# 3. 生成 XPM
gmx sham -f sham.xvg -ls FEL_sham.xpm -b 8500 -tsham 310 -nlevels 100
# 4. 绘制 2D 图
python ../scripts/xpm2png.py -f FEL_sham.xpm -o FEL_2D.png -ip yes
8.2 3D 形貌图（推荐方案）
# 1-2. 同上
# 3. 生成低 levels 的 XPM（供 3D 使用）
gmx sham -f sham.xvg -ls FEL_3d.xpm -b 8500 -tsham 310 -nlevels 25
# 4. 转换为 TXT
python ../scripts/xpm2txt.py -f FEL_3d.xpm -o FEL_3d.txt
# 5. 绘制 3D 图
python ../scripts/3D自由能形貌图.py
8.3 快速方案（使用 DuIvyTools）
pip install DuIvyTools
gmx sham -f sham.xvg -ls FEL_sham.xpm -b 8500 -tsham 310 -nlevels 100
# 2D
dit xpm_show -f FEL_sham.xpm -o FEL_2D.pdf -cmap jet
# 3D
dit xpm_show -f FEL_sham.xpm -3d -o FEL_3D.pdf
💬 留言