02 — 轨迹分析与可视化

本文档介绍 MD 模拟结束后如何对轨迹进行周期性矫正，以及如何计算并可视化 RMSD、回旋半径、氢键数量、RMSF 等关键指标。

---

前置条件

完成生产模拟后，工作目录应包含以下文件：

文件	说明
md100.tpr	模拟参数文件（拓扑+参数）
md100.xtc	压缩轨迹文件（坐标数据）
md100.log	模拟日志

---

1. 轨迹周期性矫正

模拟过程中分子可能跨越盒子边界，导致视觉上"断裂"。需要将轨迹平移到连续盒子中：

gmx trjconv -s md100.tpr -f md100.xtc -o md100_nojump.xtc -pbc nojump -ur compact

交互：

0    # 选择 System 组（整个体系）输出

参数	含义
-pbc nojump	去除跨边界跳跃，使分子轨迹连续
-ur compact	将所有分子紧凑排列在盒子中心

提示：后续所有分析都基于矫正后的 md100_nojump.xtc，不要使用原始 md100.xtc。

---

2. RMSD 分析

RMSD（均方根偏差）衡量模拟结构与参考结构之间的构象差异。

2.1 计算蛋白质 RMSD

gmx rms -s md100.tpr -f md100_nojump.xtc -o rmsd.xvg -tu ns

交互：

1    # 选择骨架 (Backbone) 作为拟合组
1    # 选择骨架 (Backbone) 作为计算组

2.2 计算配体 RMSD

gmx rms -s md100.tpr -f md100_nojump.xtc -o rmsd_ligand.xvg -tu ns

交互：

13   # 选择 UNK（配体）作为拟合组
13   # 选择 UNK（配体）作为计算组

参数	含义
-tu ns	时间单位以纳秒（ns）输出
-o rmsd.xvg	输出文件（XVG 格式，可用 Python 绘图）

组号（1、13 等）取决于 index.ndx 中的定义，请根据实际编号调整。

---

3. 回旋半径分析

回旋半径（Rg）衡量蛋白质的紧凑程度：

gmx gyrate -s md100.tpr -f md100_nojump.xtc -o gyrate.xvg

交互：

1    # 选择蛋白质
1    # 选择蛋白质

---

4. 氢键分析

计算蛋白质与配体之间的氢键数量随时间的变化：

gmx hbond -f md100_nojump.xtc -s md100.tpr -num hbond_num.xvg -life hbond_life.xvg -dist hbond_dist.xvg

交互：

1    # 选择蛋白质
13   # 选择配体 (UNK)

参数	含义
-num hbond_num.xvg	氢键数量随时间变化
-life hbond_life.xvg	氢键寿命
-dist hbond_dist.xvg	氢键距离分布

---

5. RMSF 分析

RMSF（均方根涨落）衡量每个残基在模拟中的柔性：

gmx rmsf -s md100.tpr -f md100_nojump.xtc -o rmsf.xvg -ox avg.pdb -res -oq bfac.pdb

参数	含义
-res	按残基输出（而非按原子）
-ox avg.pdb	输出平均结构 PDB
-oq bfac.pdb	输出 B-factor 着色的 PDB（可用 PyMOL 可视化）

---

6. 导出模拟快照

从轨迹中提取特定时间点的结构：

# 导出 0 ns（起始）结构
gmx trjconv -s md100.tpr -f md100_nojump.xtc -o start.pdb -dump 0

# 导出 50 ns 结构
gmx trjconv -s md100.tpr -f md100_nojump.xtc -o 50ns.pdb -dump 50

# 导出最终结构
gmx trjconv -s md100.tpr -f md100_nojump.xtc -o final.pdb -dump 100

交互：

0    # 选择 System 组

---

7. 使用 Python 绘制分析图表

7.1 绘制 RMSD 曲线

使用脚本 scripts/rmsd曲线图.py：

cd /path/to/your/project
python ../scripts/rmsd曲线图.py

功能： - 读取 GROMACS 输出的 rmsd.xvg（自动跳过注释行） - 绘制 RMSD vs Time 曲线 - 自动保存为 rmsd_plot.png（300 DPI）

7.2 绘制氢键数量曲线

使用脚本 scripts/氢键数量.py：

python ../scripts/氢键数量.py

功能： - 读取 hbond_num.xvg - 将时间从皮秒(ps)转换为纳秒(ns) - 绘制氢键数量随时间变化曲线 - 自动保存为 hbond_num.png（300 DPI）

---

分析流程速查

# 1. 轨迹矫正
gmx trjconv -s md100.tpr -f md100.xtc -o md100_nojump.xtc -pbc nojump -ur compact
# 选 0

# 2. RMSD
gmx rms -s md100.tpr -f md100_nojump.xtc -o rmsd.xvg -tu ns
# 选 1, 1

# 3. 配体 RMSD
gmx rms -s md100.tpr -f md100_nojump.xtc -o rmsd_ligand.xvg -tu ns
# 选 13, 13

# 4. 回旋半径
gmx gyrate -s md100.tpr -f md100_nojump.xtc -o gyrate.xvg
# 选 1, 1

# 5. 氢键
gmx hbond -f md100_nojump.xtc -s md100.tpr -num hbond_num.xvg

# 6. RMSF
gmx rmsf -s md100.tpr -f md100_nojump.xtc -o rmsf.xvg -res

# 7. Python 绘图
python ../scripts/rmsd曲线图.py
python ../scripts/氢键数量.py