MDAnalysis：MD 轨迹分析工具

MDAnalysis 是分析分子动力学轨迹的主力 Python 库。它能读几乎所有主流格式（DCD、XTC、TRR、NetCDF、PDB、GRO…），提供强大的原子选择语法和一整套分析模块。跑完 MD 之后的所有定量分析，基本都从它开始。

安装与基本对象

pip install MDAnalysis MDAnalysisTests

import MDAnalysis as mda

u = mda.Universe("topology.pdb", "traj.dcd")   # 拓扑 + 轨迹
print(u.atoms.n_atoms, u.trajectory.n_frames)
print(u.trajectory.dt, "ps/frame")

protein = u.select_atoms("protein")
ligand  = u.select_atoms("resname MOL")

for ts in u.trajectory[::10]:      # 每 10 帧取一帧
    print(ts.frame, protein.center_of_mass())

核心概念只有三个：Universe（体系）、AtomGroup（原子子集）、trajectory（可迭代的帧序列）。轨迹是惰性读取的——每次迭代只把当前帧载入内存，所以能处理远大于内存的轨迹文件。

选择语法：最值得先学的部分

表达式	选中什么
protein and name CA	蛋白的所有 Cα
resname MOL	配体（按残基名）
around 5 resname MOL	配体周围 5 Å 内的原子
byres around 4 protein	靠近蛋白的完整残基
resid 45:60 and backbone	45–60 号残基的主链
not (name H* or resname SOL)	去掉氢和水
sphlayer 3 8 protein	球壳层选择

around 这类几何选择默认是静态的（按当前帧计算一次）。要让选择随轨迹变化，需加 updating=True：u.select_atoms("around 5 resname MOL", updating=True)。这是分析结合位点水分子时最容易踩的坑。

四个最常用的分析

from MDAnalysis.analysis import rms, align, distances, hydrogenbonds
import numpy as np

# 1) RMSD：结构相对参考的偏离
R = rms.RMSD(u, u, select="backbone",
             groupselections=["resname MOL"]).run()
time, rmsd_bb, rmsd_lig = R.rmsd[:, 1], R.rmsd[:, 2], R.rmsd[:, 3]

# 2) RMSF：逐残基柔性（需先对齐）
align.AlignTraj(u, u, select="protein and name CA", in_memory=True).run()
ca = u.select_atoms("protein and name CA")
rmsf = rms.RMSF(ca).run().rmsf

# 3) 关键距离随时间（如配体与催化残基）
d = []
for ts in u.trajectory:
    d.append(distances.dist(ligand, u.select_atoms("resid 145 and name OG"))[2][0])

# 4) 氢键分析
hb = hydrogenbonds.HydrogenBondAnalysis(
        u, between=["resname MOL", "protein"]).run()
print(hb.count_by_ids()[:10])       # 各氢键出现的帧数

怎么读这些数字

• 主链 RMSD：稳定体系通常在 1～3 Å 后进入平台。持续上升说明还没平衡，或蛋白正在发生构象变化。
• 配体 RMSD：对齐蛋白后看配体的 RMSD。若持续 < 2 Å，说明结合姿势稳定；若跳到 5 Å 以上并不再回落，往往意味着对接姿势不正确——这是用 MD 验证对接结果的核心判据。
• RMSF：识别柔性区域。loop 与末端高是正常的，若关键结合残基 RMSF 很高，说明相互作用不稳定。
• 氢键占有率：某个氢键在轨迹中出现的帧数比例。>70% 算稳定相互作用，<30% 基本是偶发接触，不该写进结合模式结论。

性能要点

• 大轨迹优先用 XTC 等压缩格式；分析前用 u.atoms.select_atoms(...).write() 抽出感兴趣的子集，能大幅提速。
• in_memory=True 会把整条轨迹载入内存，只对小体系可行。
• 逐帧 Python 循环慢，尽量用内置分析类（底层有向量化实现）；确需自定义时考虑 MDAnalysis.analysis.base.AnalysisBase 配合并行。
• 周期性边界会让分子被「切开」，计算距离前需要正确处理 PBC（用 transformations 模块做 unwrap）。

延伸资源

• 模拟引擎：200《OpenMM》；轻量替代：206《MDTraj》；
• 相互作用指纹：207《ProLIF》、208《PLIP》；
• MD 概念与采样问题见「结构与模拟」模块。