MDAnalysis 分析轨迹：从 RMSD 到相互作用统计

MDAnalysis 能读取几乎所有主流 MD 软件的轨迹格式，提供统一的 Python 接口。它的核心是原子选择语法与可组合的分析模块——掌握这两点，就能应对绝大多数轨迹分析需求。

基本概念

pip install MDAnalysis MDAnalysisTests

import MDAnalysis as mda
import numpy as np

# Universe = 拓扑 + 轨迹
u = mda.Universe("topology.pdb", "trajectory.dcd")

print(f"原子数: {len(u.atoms)}")
print(f"残基数: {len(u.residues)}")
print(f"帧数: {len(u.trajectory)}")
print(f"时间步: {u.trajectory.dt} ps")

# 支持的格式：
#   拓扑: PDB, GRO, PSF, PRMTOP, TPR, ...
#   轨迹: DCD, XTC, TRR, NetCDF, ...

# 【重要】：轨迹是按需读取的（惰性加载）
#   u.trajectory[10]  会跳到第 10 帧
#   遍历时每次只有一帧在内存中
#   → 这让处理超大轨迹成为可能
#   → 但也意味着不能同时访问多帧的坐标

选择语法：最需要掌握的部分

# 基础选择
protein = u.select_atoms("protein")
backbone = u.select_atoms("backbone")          # N, CA, C, O
ca = u.select_atoms("name CA")
ligand = u.select_atoms("resname LIG")
water = u.select_atoms("resname HOH SOL TIP3")

# 按残基/链
res = u.select_atoms("resid 100-150")
chain_a = u.select_atoms("segid A")
specific = u.select_atoms("resname TYR and resid 123")

# 【几何选择】—— 最常用也最强大
around = u.select_atoms("around 5.0 resname LIG")
#   距配体 5 Å 内的原子（不含配体本身）

pocket = u.select_atoms("byres (around 5.0 resname LIG)")
#   byres 把选择扩展到完整残基 —— 【几乎总是需要的】

sphere = u.select_atoms("sphlayer 0 8 resname LIG")
#   以配体为中心的球层

# 组合与逻辑
sel = u.select_atoms("protein and not name H*")
sel = u.select_atoms("(resname ARG LYS) and name NZ NH1 NH2")
sel = u.select_atoms("protein and around 5 resname LIG", updating=True)
#   updating=True 让选择【每帧动态更新】
#   → 分析「随时间变化的邻近残基」时必需
#   → 但会明显变慢

# 【常见错误】：
#   忘记 byres → 只选到部分原子，导致结果错误
#   忘记 updating → 用第一帧的选择分析整条轨迹
#
# 【调试技巧】：选择后先打印检查
print(f"选中 {len(pocket)} 个原子, "
      f"{len(pocket.residues)} 个残基")
print(sorted(set(pocket.resids)))

常用分析

from MDAnalysis.analysis import rms, distances, hydrogenbonds, contacts

# ---- 距离与接触 ----
lig = u.select_atoms("resname LIG")
res = u.select_atoms("resid 123 and name CA")

dists = []
for ts in u.trajectory:
    d = np.linalg.norm(lig.center_of_mass() - res.positions[0])
    dists.append(d)
dists = np.array(dists)
print(f"平均距离: {dists.mean():.2f} ± {dists.std():.2f} Å")

# ---- 氢键分析 ----
hbonds = hydrogenbonds.HydrogenBondAnalysis(
    universe=u,
    donors_sel="resname LIG",
    hydrogens_sel="resname LIG and name H*",
    acceptors_sel="protein and name O* N*",
    d_a_cutoff=3.5,
    d_h_a_angle_cutoff=150,
)
hbonds.run()

# 计算每个氢键的占有率
counts = hbonds.count_by_ids()
n_frames = len(u.trajectory)
for donor, hydrogen, acceptor, count in counts:
    print(f"占有率 {count / n_frames:.1%}: "
          f"{u.atoms[donor].name} → {u.atoms[acceptor].resname}"
          f"{u.atoms[acceptor].resid}")

# ---- 原生接触（判断结构是否保持）----
ref = u.copy()
q = contacts.Contacts(u,
                      select=("protein and name CA", "protein and name CA"),
                      refgroup=(ref.select_atoms("protein and name CA"),
                                ref.select_atoms("protein and name CA")),
                      radius=8.0).run()
# q.results.timeseries[:, 1] 是保持的接触比例

# ---- 回旋半径（紧密程度）----
rg = []
for ts in u.trajectory:
    rg.append(protein.radius_of_gyration())

# ---- 二级结构（需要 MDTraj）----
import mdtraj as md
traj = md.load("trajectory.dcd", top="topology.pdb")
dssp = md.compute_dssp(traj, simplified=True)
helix_fraction = (dssp == "H").mean(axis=1)
print(f"平均螺旋含量: {helix_fraction.mean():.1%}")

结合位点水分子分析

# 这是一个实用但少见于教程的分析
# 目标：找出结合位点中稳定存在的水分子（见 091）

from collections import Counter

water_counts = Counter()
n_frames = 0

for ts in u.trajectory[::5]:
    # updating 选择：每帧重新找配体附近的水
    near_water = u.select_atoms(
        "resname HOH and around 4.0 resname LIG")
    for res in near_water.residues:
        water_counts[res.resid] += 1
    n_frames += 1

print("稳定的结合位点水（占有率 > 50%）:")
for resid, count in water_counts.most_common(20):
    occ = count / n_frames
    if occ > 0.5:
        print(f"  HOH {resid}: {occ:.1%}")

# 判读：
#   高占有率的水 = 结构水，设计时应保留或谨慎置换
#   低占有率的水 = 交换频繁，可能是可置换的
#
# 注意：
#   水分子会互相交换位置，
#   按 resid 统计可能低估某个「位点」的占有率
#   → 更严格的做法是按空间位置聚类而非按分子 ID

性能优化

# 问题：大轨迹的分析可能很慢

# 1) 抽帧
for ts in u.trajectory[::10]:       # 每 10 帧取一次
    ...
# 相邻帧高度相关，抽帧几乎不损失信息

# 2) 只在需要时用 updating 选择
#    updating=True 每帧重算选择，很慢
#    如果口袋残基基本不变，用静态选择即可

# 3) 预先切出感兴趣的部分，写成小轨迹
sel = u.select_atoms("protein or resname LIG")
with mda.Writer("reduced.dcd", sel.n_atoms) as W:
    for ts in u.trajectory:
        W.write(sel)
# 去掉水分子后，文件可能小 10 倍，后续分析快得多

# 4) 用 in_memory 加载（小轨迹）
u = mda.Universe("top.pdb", "traj.dcd", in_memory=True)
# 全部载入内存，随机访问快
# 但大轨迹会 OOM

# 5) 并行分析（新版支持）
from MDAnalysis.analysis import rms
R = rms.RMSD(u, select="backbone")
R.run(n_workers=4, backend="multiprocessing")

# 6) 用 MDTraj 做某些计算
#    MDTraj 的某些分析（如 DSSP、SASA）实现更快
#    两个库可以配合使用

与 ProLIF 的配合

MDAnalysis 提供轨迹读取与原子选择，ProLIF（见 109《ProLIF》）在其之上做相互作用指纹分析——这是分析配体结合稳定性的标准组合。ProLIF 直接接受 MDAnalysis 的 AtomGroup 作为输入，两者无缝衔接。

延伸资源

• RMSD/RMSF：124《RMSD 与 RMSF》；ProLIF：109《ProLIF》；OpenMM：128《OpenMM 跑 10–100 ns MD》；
• MD 入门：123《分子动力学 MD 入门》；口袋与水：091《蛋白口袋 Protein Pocket》。