mdanalysis-trajectory

MDAnalysis — 分子动力学轨迹分析

概述

MDAnalysis 提供统一的 Python 接口，用于读取和分析不同分子动力学引擎（GROMACS、AMBER、NAMD、CHARMM、LAMMPS、OpenMM）生成的轨迹。它将分子系统表示为 Universe 对象，其中包含带有位置、速度、力和拓扑数据的 AtomGroup。可以逐帧遍历轨迹，也可以使用分析模块批量分析 RMSD、RMSF、回转半径、氢键、溶剂可及表面积和 PCA。MDAnalysis 可与 NumPy、pandas 和 matplotlib 集成，是计算化学和药物发现领域进行模拟后结构分析的标准工具。

适用场景

• 分子动力学模拟后，计算蛋白质主链或特定残基组的 RMSD 和 RMSF
• 分析配体结合稳定性：口袋 RMSD、接触持久性、氢键占有率
• 对轨迹构象进行主成分分析（PCA）
• 计算溶剂可及表面积（SASA）、回转半径和端到端距离
• 从长分子动力学轨迹中提取具有代表性的聚类结构，以便可视化
• 在 HPC 工作流中进行引擎专属分析时，改用 GROMACS 或 AMBER 分析工具（gmx rms、cpptraj）
• 直接使用 OpenMM 或 GROMACS 运行分子动力学模拟；MDAnalysis 用于模拟后分析

前置条件

• Python 包：MDAnalysis、numpy、matplotlib、pandas
• 输入：拓扑文件（.psf、.prmtop、.gro、.pdb）+ 轨迹文件（.dcd、.trr、.xtc、.nc、.dms）

# Install MDAnalysis
pip install MDAnalysis

# Install with all analysis extras
pip install "MDAnalysis[analysis]"

# Verify
python -c "import MDAnalysis as mda; print(mda.__version__)"
# 2.7.0

快速开始

import MDAnalysis as mda
import numpy as np

# Load a GROMACS topology + trajectory
u = mda.Universe("protein.gro", "trajectory.xtc")

print(f"Atoms: {u.atoms.n_atoms}")
print(f"Residues: {u.residues.n_residues}")
print(f"Frames: {u.trajectory.n_frames}")
print(f"First frame positions (first 3 atoms):\n{u.atoms.positions[:3]}")

核心 API

模块 1：Universe 和 AtomGroup — 加载和选择原子

加载轨迹并选择原子子集。

import MDAnalysis as mda

# Load topology + trajectory (GROMACS xtc format)
u = mda.Universe("system.gro", "md_production.xtc")

# AtomGroup selections (CHARMM-style selection language)
protein = u.select_atoms("protein")
backbone = u.select_atoms("backbone")
ca_atoms = u.select_atoms("name CA")
ligand = u.select_atoms("resname LIG")
binding_site = u.select_atoms("protein and around 5.0 resname LIG")

print(f"Protein atoms: {protein.n_atoms}")
print(f"CA atoms: {ca_atoms.n_atoms}")
print(f"Ligand atoms: {ligand.n_atoms}")
print(f"Binding site residues: {binding_site.residues.n_residues}")

# Access atom properties at current frame
print(f"CA positions shape: {ca_atoms.positions.shape}")  # (N, 3)
print(f"Protein mass: {protein.total_mass():.1f} Da")

模块 2：轨迹遍历 — 逐帧分析

遍历轨迹帧，进行时间序列分析。

import MDAnalysis as mda
import numpy as np

u = mda.Universe("protein.gro", "trajectory.xtc")
backbone = u.select_atoms("backbone")

times = []
rg_values = []

for ts in u.trajectory:
    times.append(u.trajectory.time)
    rg_values.append(backbone.radius_of_gyration())

import pandas as pd
df = pd.DataFrame({"time_ps": times, "Rg_A": rg_values})
print(f"Frames analyzed: {len(df)}")
print(f"Mean Rg: {df['Rg_A'].mean():.2f} Å")
print(f"Rg std: {df['Rg_A'].std():.2f} Å")
df.to_csv("radius_of_gyration.csv", index=False)

模块 3：RMSD 分析 — 随时间变化的结构漂移

计算相对于参考结构的主链 RMSD。

import MDAnalysis as mda
from MDAnalysis.analysis import rms
import numpy as np
import matplotlib.pyplot as plt

u = mda.Universe("protein.gro", "trajectory.xtc")

# RMSD of Cα atoms relative to first frame
rmsd = rms.RMSD(u, select="name CA")
rmsd.run()

# Results: frame, time (ps), RMSD (Å)
results = rmsd.results.rmsd
print(f"Mean RMSD: {results[:, 2].mean():.2f} Å")
print(f"Max RMSD: {results[:, 2].max():.2f} Å")

# Plot
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(results[:, 1] / 1000, results[:, 2], color="steelblue", lw=0.8)
ax.set_xlabel("Time (ns)")
ax.set_ylabel("RMSD (Å)")
ax.set_title("Backbone RMSD")
plt.tight_layout()
plt.savefig("rmsd.png", dpi=150)
print("Saved: rmsd.png")

模块 4：RMSF 分析 — 残基柔性

计算均方根涨落，以识别柔性区域。

import MDAnalysis as mda
from MDAnalysis.analysis import rms
import numpy as np
import matplotlib.pyplot as plt

u = mda.Universe("protein.gro", "trajectory.xtc")

# RMSF per Cα atom (after aligning trajectory)
ca_atoms = u.select_atoms("name CA")
rmsf_analysis = rms.RMSF(ca_atoms)
rmsf_analysis.run()

rmsf_values = rmsf_analysis.results.rmsf
resids = ca_atoms.resids

print(f"Most flexible residue: {resids[np.argmax(rmsf_values)]} ({rmsf_values.max():.2f} Å)")
print(f"Most rigid residue:    {resids[np.argmin(rmsf_values)]} ({rmsf_values.min():.2f} Å)")

# Plot B-factor-like profile
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(resids, rmsf_values, color="coral", lw=1)
ax.fill_between(resids, 0, rmsf_values, alpha=0.3, color="coral")
ax.set_xlabel("Residue ID")
ax.set_ylabel("RMSF (Å)")
ax.set_title("Per-residue RMSF")
plt.tight_layout()
plt.savefig("rmsf.png", dpi=150)

模块 5：氢键分析

识别并统计轨迹中蛋白质与配体之间的氢键。

import MDAnalysis as mda
from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis
import pandas as pd

u = mda.Universe("complex.gro", "trajectory.xtc")

# Protein-ligand hydrogen bond analysis
hbonds = HydrogenBondAnalysis(
    universe=u,
    donors_sel="protein",
    acceptors_sel="resname LIG",
    d_h_cutoff=1.2,      # donor-hydrogen distance cutoff (Å)
    d_a_cutoff=3.0,      # donor-acceptor distance cutoff (Å)
    d_h_a_angle_cutoff=150.0,  # angle cutoff (degrees)
)
hbonds.run()

# Get occupancy: fraction of frames with each H-bond
hbonds.generate_table()
df = pd.DataFrame(hbonds.table)
print(f"Total unique H-bonds observed: {len(df['donor_resid'].unique())}")

# Count H-bonds per frame
counts = hbonds.count_by_time()
print(f"Mean H-bonds per frame: {counts[:, 1].mean():.2f}")
print(f"Max H-bonds in one frame: {counts[:, 1].max():.0f}")

模块 6：PCA 和构象聚类

提取主要运动模式并对构象进行聚类。

import MDAnalysis as mda
from MDAnalysis.analysis import pca, align
import numpy as np
import matplotlib.pyplot as plt

u = mda.Universe("protein.gro", "trajectory.xtc")

# Align trajectory to first frame
aligner = align.AlignTraj(u, u, select="backbone", in_memory=True)
aligner.run()

# PCA on backbone Cα atoms
ca = u.select_atoms("name CA")
pc = pca.PCA(u, select="backbone")
pc.run()

# Explained variance
cumvar = np.cumsum(pc.results.variance / pc.results.variance.sum())
n_for_90 = np.searchsorted(cumvar, 0.90) + 1
print(f"PCs to explain 90% variance: {n_for_90}")
print(f"PC1 variance: {pc.results.variance[0] / pc.results.variance.sum() * 100:.1f}%")

# Project trajectory onto PC1-PC2 space
transformed = pc.transform(ca, n_components=2)
plt.scatter(transformed[:, 0], transformed[:, 1], c=range(len(transformed)),
            cmap="viridis", s=2)
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.colorbar(label="Frame")
plt.title("PCA Conformational Landscape")
plt.savefig("pca.png", dpi=150)
print("Saved: pca.png")

关键参数

参数	模块	默认值	作用
select (Universe)	AtomGroup	—	MDAnalysis 选择语言字符串（例如，"protein and name CA"）
in_memory	align.AlignTraj	False	将完整轨迹载入 RAM，以加快分析速度
step	trajectory loop	1	每隔 N 帧分析一帧；step=10 可将计算量减少 10 倍
start/stop	analysis.run()	all	帧范围；run(start=100, stop=500) 分析第 100 至 500 帧
d_a_cutoff	HydrogenBondAnalysis	3.0	供体-受体距离截断值（Å）
d_h_a_angle_cutoff	HydrogenBondAnalysis	150.0	氢键角度截断值（度）
n_components	PCA.transform	all	要投影到的主成分数量
groupselection	RMSD	None	用于拟合的次级选择；RMSD 计算使用主要选择
weights	RMSD/align	None	使用 "mass" 进行质量加权 RMSD
verbose	analysis.run()	False	分析期间显示进度条

常见工作流

工作流 1：完整的蛋白质-配体结合稳定性分析

import MDAnalysis as mda
from MDAnalysis.analysis import rms, align
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

u = mda.Universe("complex.gro", "md_100ns.xtc")

# Align trajectory to initial frame
align.AlignTraj(u, u, select="protein and backbone", in_memory=True).run()

protein_bb = u.select_atoms("backbone")
ligand = u.select_atoms("resname LIG")

results = {"time_ns": [], "protein_rmsd": [], "ligand_rmsd": [], "pocket_rg": []}
ref_positions = {"protein": protein_bb.positions.copy(), "ligand": ligand.positions.copy()}

for ts in u.trajectory:
    results["time_ns"].append(ts.time / 1000)
    results["protein_rmsd"].append(
        rms.rmsd(protein_bb.positions, ref_positions["protein"], superposition=False))
    results["ligand_rmsd"].append(
        rms.rmsd(ligand.positions, ref_positions["ligand"], superposition=False))
    results["pocket_rg"].append(
        u.select_atoms("protein and around 6.0 resname LIG").radius_of_gyration())

df = pd.DataFrame(results)
df.to_csv("binding_stability.csv", index=False)
print(f"Ligand mean RMSD: {df['ligand_rmsd'].mean():.2f} Å")
print(f"Pocket stable: {'Yes' if df['pocket_rg'].std() < 0.5 else 'No'}")

工作流 2：提取 RMSD 最小的代表性结构

import MDAnalysis as mda
from MDAnalysis.analysis import rms
import numpy as np

u = mda.Universe("protein.gro", "trajectory.xtc")

# Compute RMSD for all frames
rmsd_analysis = rms.RMSD(u, select="backbone")
rmsd_analysis.run()
rmsd_values = rmsd_analysis.results.rmsd[:, 2]

# Find frame with lowest RMSD (most representative)
min_frame = np.argmin(rmsd_values)
print(f"Most representative frame: {min_frame} (RMSD: {rmsd_values[min_frame]:.2f} Å)")

# Extract and save that frame as PDB
u.trajectory[min_frame]
with mda.Writer("representative_structure.pdb", u.atoms.n_atoms) as writer:
    writer.write(u.atoms)
print("Saved: representative_structure.pdb")

常见操作示例

操作示例 1：计算两个蛋白质结构域之间的接触图

import MDAnalysis as mda
from MDAnalysis.analysis import contacts
import numpy as np

u = mda.Universe("protein.gro", "trajectory.xtc")

# Define two domains
domain_A = u.select_atoms("resid 1-100 and name CA")
domain_B = u.select_atoms("resid 200-300 and name CA")

# Count domain-domain contacts (< 8 Å) over time
contact_counts = []
for ts in u.trajectory[::10]:  # every 10th frame
    dist_matrix = np.sqrt(np.sum(
        (domain_A.positions[:, None] - domain_B.positions[None, :]) ** 2, axis=-1
    ))
    contact_counts.append((dist_matrix < 8.0).sum())

print(f"Mean contacts: {np.mean(contact_counts):.1f}")
print(f"Contact count range: {min(contact_counts)} – {max(contact_counts)}")

操作示例 2：将轨迹子集写入新文件

import MDAnalysis as mda

u = mda.Universe("system.gro", "long_trajectory.xtc")

# Write only protein atoms, every 10th frame, frames 500-2000
protein = u.select_atoms("protein")
with mda.Writer("protein_subset.xtc", protein.n_atoms) as writer:
    for ts in u.trajectory[500:2000:10]:
        writer.write(protein)

print(f"Subset trajectory written: {(2000-500)//10} frames")

故障排除

问题	原因	解决方案
ValueError: Universe has no file	缺少轨迹参数	同时提供拓扑文件和轨迹文件：mda.Universe(top, traj)
对齐后 RMSD 仍然漂移	用于拟合的选择不正确	使用 "backbone" 进行拟合；确认拓扑与轨迹匹配
KeyError: 'resname LIG'	残基名称不是标准名称	检查 u.residues.resnames；使用拓扑中的确切名称
帧迭代速度很慢	内存中的轨迹过大	使用 step=10 跳过帧；将 xtc 转换为较小的 DCD 文件
氢键计数为 0	拓扑中没有显式氢原子	使用含有显式 H 原子的拓扑；使用 psfgen 或 tleap 添加氢原子
NoDataError: xtc has no charges	轨迹中没有该属性	使用包含电荷的拓扑文件（.psf、.prmtop、带有 itp 的 .gro）
使用 in_memory=True 时发生内存错误	轨迹对于 RAM 来说过大	移除 in_memory=True；增加 RAM；或分块处理
Apple Silicon 上导入错误	二进制文件未针对 arm64 编译	通过 conda-forge 安装：conda install -c conda-forge MDAnalysis

参考资料

• MDAnalysis 文档 — 官方 API 参考和用户指南
• MDAnalysis GitHub：MDAnalysis/mdanalysis — 源代码、教程和问题跟踪器
• Michaud-Agrawal N 等（2011）“MDAnalysis：用于分子动力学模拟分析的工具包” — J Comput Chem 32:2319-2327。DOI:10.1002/jcc.21787
• Gowers RJ 等（2016）“MDAnalysis：用于快速分析分子动力学模拟的 Python 包” — Proc SciPy 2016。DOI:10.25080/majora-629e541a-00e