mdtraj-trajectory-analysis

mdtraj 轨迹分析

概述

mdtraj 是一个依赖较少的 Python 库，用于分析 MD 轨迹。它可将 DCD/XTC/TRR/NetCDF/H5/AMBER/GROMACS/CHARMM/OpenMM 文件读入由 NumPy 数组支持的 Trajectory 对象，然后提供几何分析、RMSD/Rg/RMSF/SASA、接触、氢键、扭转角和 8 状态 DSSP 等纯 Python 函数。

单位：mdtraj 内部使用 nm 和 ps。距离乘以 10 可转换为 Å，时间除以 1000 可转换为 ns。扭转角以 弧度 为单位——使用 np.degrees() 转换。

适用场景

• 分析 RMSD 随时间的变化、Rg，以及各残基的 RMSF，以比较不同 MD 重复模拟的稳定性
• 根据轨迹集合生成残基间接触频率图
• 计算 Ramachandran 图（常规、Gly、Pro）所需的主链 phi/psi
• 逐帧计算各残基的 8 状态 DSSP，生成二级结构时间序列
• 需要轻量级临时分析，而 MDAnalysis 的完整框架功能过于繁重时
• NMR 观测量（J 耦合、通过 SHIFTX2 计算化学位移）
• 如果需要 MDAnalysis 的选择语法、AnalysisBase 并行功能，或 LAMMPS/NAMD 专用读取器，请改用 mdanalysis-trajectory
• 使用 OpenMM/GROMACS 运行模拟——此技能仅用于模拟后的分析

前提条件

• mdtraj、numpy、pandas、matplotlib
• 拓扑文件（PDB/GRO/PRMTOP/PSF）+ 轨迹文件（DCD/XTC/TRR/NetCDF/H5）
• Python 3.9+，建议使用 conda-forge 构建版本

安装前请先检查——在 pixi/conda 环境中通常已安装 mdtraj：

python3 -c "import mdtraj" 2>/dev/null || conda install -c conda-forge mdtraj numpy pandas matplotlib

快速入门

import mdtraj as md

traj = md.load("traj.xtc", top="topology.pdb")
ca = traj.topology.select("name CA")
traj.superpose(traj, frame=0, atom_indices=ca)

rmsd_ang = md.rmsd(traj, traj, frame=0, atom_indices=ca) * 10.0  # nm -> Å
print(f"Frames: {traj.n_frames}  RMSD: {rmsd_ang.min():.2f}–{rmsd_ang.max():.2f} Å")

核心 API

模块 1：轨迹 I/O

可整体加载，也可流式加载。根据扩展名自动检测格式。

import mdtraj as md

traj = md.load("rep1.xtc", top="protein.pdb")

# Stream large trajectories — avoids OOM
for chunk in md.iterload("rep1.xtc", top="protein.pdb", chunk=500):
    rmsd_chunk = md.rmsd(chunk, chunk, frame=0)

# Save subset
ca = traj.topology.select("name CA")
traj.atom_slice(ca).save_dcd("ca_only.dcd")

选择原子并切分帧：

backbone = traj.topology.select("backbone")
chain_a  = traj.topology.select("chainid 0")

first_ns      = traj[:1000]
every_10th    = traj[::10]
last_half_bb  = traj[traj.n_frames // 2:].atom_slice(backbone)

模块 2：RMSD、Rg、RMSF

md.rmsd 会在内部进行叠合；计算 RMSF 时，需先显式叠合，再手动计算。

import mdtraj as md, numpy as np

traj = md.load("rep1.xtc", top="protein.pdb")
ca = traj.topology.select("name CA")

rmsd_ang = md.rmsd(traj, traj, frame=0, atom_indices=ca) * 10.0    # Å
rg_ang   = md.compute_rg(traj) * 10.0                              # Å
time_ns  = traj.time / 1000.0
print(f"<RMSD>={rmsd_ang.mean():.2f} Å,  <Rg>={rg_ang.mean():.2f} Å")

# Per-CA RMSF — average-structure reference
ca_traj = traj.atom_slice(ca)
ca_traj.superpose(ca_traj, frame=0)
diff     = ca_traj.xyz - ca_traj.xyz.mean(axis=0)
rmsf_ang = np.sqrt((diff ** 2).sum(axis=2).mean(axis=0)) * 10.0
res_ids  = [a.residue.resSeq for a in ca_traj.topology.atoms]
print(f"Max RMSF: residue {res_ids[np.argmax(rmsf_ang)]} = {rmsf_ang.max():.2f} Å")

模块 3：接触和距离图

通过距离阈值计算接触频率。

import mdtraj as md, numpy as np

traj = md.load("rep1.xtc", top="protein.pdb")
distances_nm, pairs = md.compute_contacts(traj, contacts="all", scheme="closest-heavy")
# distances_nm: (n_frames, n_pairs);  pairs: (n_pairs, 2) of residue indices

contact_freq = (distances_nm < 0.5).mean(axis=0)                  # 5 Å cutoff
n_res    = traj.n_residues
freq_map = np.zeros((n_res, n_res))
for (i, j), f in zip(pairs, contact_freq):
    freq_map[i, j] = freq_map[j, i] = f
print(f"Persistent contacts (>0.8): {(contact_freq > 0.8).sum()}")

模块 4：扭转角（Ramachandran）

phi/psi 以弧度返回；由于首个残基没有 phi，末个残基没有 psi，因此需取交集。

import mdtraj as md, numpy as np

traj = md.load("rep1.xtc", top="protein.pdb")
phi_ix, phi_rad = md.compute_phi(traj)
psi_ix, psi_rad = md.compute_psi(traj)

def res_of(indices): return np.array([traj.topology.atom(ix[1]).residue.index for ix in indices])
phi_res, psi_res = res_of(phi_ix), res_of(psi_ix)
common  = np.intersect1d(phi_res, psi_res)
phi_deg = np.degrees(phi_rad[:, np.isin(phi_res, common)])
psi_deg = np.degrees(psi_rad[:, np.isin(psi_res, common)])

筛选 Gly / Pro 残基：

res_names = [traj.topology.residue(r).name for r in common]
gly_cols  = [i for i, n in enumerate(res_names) if n == "GLY"]
pro_cols  = [i for i, n in enumerate(res_names) if n == "PRO"]
phi_gly, psi_gly = phi_deg[:, gly_cols].ravel(), psi_deg[:, gly_cols].ravel()
phi_pro, psi_pro = phi_deg[:, pro_cols].ravel(), psi_deg[:, pro_cols].ravel()

模块 5：DSSP 二级结构（8 状态）

md.compute_dssp(traj, simplified=False) 返回由单字符代码组成的 (n_frames, n_residues) 数组：

代码	含义
H	α 螺旋
B	β 桥
E	β 折叠
G	3_10 螺旋
I	π 螺旋
T	转角
S	弯曲
(空格) 或 C	无规卷曲

import mdtraj as md

traj = md.load("rep1.xtc", top="protein.pdb")
dssp = md.compute_dssp(traj, simplified=False)
dssp[dssp == " "] = "C"
print(f"DSSP grid: {dssp.shape}  unique codes: {set(dssp.flatten())}")

模块 6：氢键和 SASA

import mdtraj as md
import pandas as pd

traj = md.load("rep1.xtc", top="protein.pdb")

# Baker–Hubbard: returns (n_hbonds, 3) of [donor, H, acceptor] atom indices
hbonds = md.baker_hubbard(traj, freq=0.5, periodic=False)
label  = lambda a: f"{traj.topology.atom(a).residue}-{traj.topology.atom(a).name}"
hbonds_df = pd.DataFrame({"donor": [label(d) for d, _, _ in hbonds],
                         "acceptor": [label(a) for _, _, a in hbonds]})

# SASA in nm^2 per residue per frame
sasa_nm2 = md.shrake_rupley(traj, mode="residue")
total_ang2 = sasa_nm2.sum(axis=1) * 100.0                # nm^2 -> Å^2
print(f"Persistent H-bonds: {len(hbonds_df)},  <SASA>: {total_ang2.mean():.0f} Ų")

核心概念

单位（nm、ps、弧度）

traj.xyz 的单位是 nm，traj.time 的单位是 ps，扭转角的单位是弧度。绘图前请先转换：

import numpy as np
xyz_ang  = traj.xyz * 10.0
time_ns  = traj.time / 1000.0
phi_deg  = np.degrees(phi_rad)

配对索引与 resSeq

md.compute_contacts(..., contacts="all") 返回的是 从 0 开始的残基索引（跨链连续编号），而不是 PDB 的 resSeq。绘制热图前，请重塑为 (n_res, n_res)。仅在添加标签时，才通过 traj.topology.residue(i).resSeq 映射回原编号。

8 状态与简化 DSSP

simplified=False → 8 种状态。simplified=True → 3 种状态（H/E/C），会舍弃 3_10 螺旋（G）、π 螺旋（I）、β 桥（B）、弯曲（S）和转角（T）。

常见工作流

工作流 1：多次重复模拟的稳定性、柔性、接触、Ramachandran 图和 DSSP 分析

目标：为三次 MD 重复模拟生成结构分析报告——包括 RMSD/Rg 随时间的变化、各残基的 RMSF、残基接触频率图（rep1）、Ramachandran 常规图 + Gly + Pro（rep1），以及 DSSP 8 状态时间序列（rep1）。绘图使用 matplotlib 默认设置；调色板由下游使用者选择。

import mdtraj as md, numpy as np, matplotlib.pyplot as plt
from pathlib import Path

replicas = {"rep1": "rep1.xtc", "rep2": "rep2.xtc", "rep3": "rep3.xtc"}
topology = "protein.pdb"
outdir = Path("figures"); outdir.mkdir(exist_ok=True)

trajs = {name: md.load(p, top=topology) for name, p in replicas.items()}
ca    = next(iter(trajs.values())).topology.select("name CA")

# 1. RMSD vs time + Rg (Å) for all replicates
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
for name, traj in trajs.items():
    traj.superpose(traj, frame=0, atom_indices=ca)
    ax1.plot(traj.time / 1000.0, md.rmsd(traj, traj, frame=0, atom_indices=ca) * 10.0, lw=1, label=name)
    ax2.plot(traj.time / 1000.0, md.compute_rg(traj) * 10.0, lw=1, label=name)
ax1.set(xlabel="Time (ns)", ylabel="RMSD (Å)", title="Backbone RMSD")
ax2.set(xlabel="Time (ns)", ylabel="Rg (Å)",   title="Radius of gyration")
for ax in (ax1, ax2): ax.legend(frameon=False)
fig.tight_layout(); fig.savefig(outdir / "01_rmsd_rg.png"); plt.close(fig)

# 2. Per-residue RMSF (Å) across replicates
fig, ax = plt.subplots(figsize=(10, 4))
for name, traj in trajs.items():
    ct = traj.atom_slice(ca); ct.superpose(ct, frame=0)
    rmsf = np.sqrt(((ct.xyz - ct.xyz.mean(axis=0)) ** 2).sum(axis=2).mean(axis=0)) * 10.0
    ax.plot([a.residue.resSeq for a in ct.topology.atoms], rmsf, lw=1, label=name)
ax.set(xlabel="Residue", ylabel="RMSF (Å)", title="Per-residue flexibility")
ax.legend(frameon=False); fig.tight_layout()
fig.savefig(outdir / "02_rmsf.png"); plt.close(fig)

# 3. Contact frequency map (rep1, 5 Å cutoff)
rep1 = trajs["rep1"]
dist_nm, pairs = md.compute_contacts(rep1, contacts="all", scheme="closest-heavy")
freq = (dist_nm < 0.5).mean(axis=0)
freq_map = np.zeros((rep1.n_residues, rep1.n_residues))
for (i, j), f in zip(pairs, freq):
    freq_map[i, j] = freq_map[j, i] = f
np.fill_diagonal(freq_map, 1.0)
fig, ax = plt.subplots(figsize=(6, 5))
im = ax.imshow(freq_map, vmin=0, vmax=1, origin="lower")
ax.set(xlabel="Residue", ylabel="Residue", title="Contact frequency (rep1)")
fig.colorbar(im, ax=ax, label="Fraction of frames in contact")
fig.tight_layout(); fig.savefig(outdir / "03_contact_map.png"); plt.close(fig)

# 4. Ramachandran (rep1) — phi/psi 2D density
phi_ix, phi_rad = md.compute_phi(rep1)
psi_ix, psi_rad = md.compute_psi(rep1)
def res_of(ix, t): return np.array([t.topology.atom(i[1]).residue.index for i in ix])
phi_res, psi_res = res_of(phi_ix, rep1), res_of(psi_ix, rep1)
common = np.intersect1d(phi_res, psi_res)
phi_deg = np.degrees(phi_rad[:, np.isin(phi_res, common)]).ravel()
psi_deg = np.degrees(psi_rad[:, np.isin(psi_res, common)]).ravel()
H, _, _ = np.histogram2d(phi_deg, psi_deg, bins=72, range=[[-180, 180], [-180, 180]])
fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(H.T, origin="lower", extent=[-180, 180, -180, 180], aspect="equal")
ax.set(xlabel=r"$\phi$ (°)", ylabel=r"$\psi$ (°)", title="Ramachandran (rep1)")
fig.tight_layout(); fig.savefig(outdir / "04_ramachandran_all.png"); plt.close(fig)

# 5. Glycine + proline Ramachandran
res_names = [rep1.topology.residue(r).name for r in common]
for label, cols, fname in [("Glycine", [i for i, n in enumerate(res_names) if n == "GLY"], "05_rama_gly.png"),
                           ("Proline", [i for i, n in enumerate(res_names) if n == "PRO"], "06_rama_pro.png")]:
    if not cols: continue
    phi = np.degrees(phi_rad[:, cols]).ravel()
    psi = np.degrees(psi_rad[:, cols]).ravel()
    H, _, _ = np.histogram2d(phi, psi, bins=60, range=[[-180, 180], [-180, 180]])
    fig, ax = plt.subplots(figsize=(5, 5))
    ax.imshow(H.T, origin="lower", extent=[-180, 180, -180, 180], aspect="equal")
    ax.set(xlabel=r"$\phi$ (°)", ylabel=r"$\psi$ (°)", title=f"{label} Ramachandran (rep1)")
    fig.tight_layout(); fig.savefig(outdir / fname); plt.close(fig)

# 6. DSSP 8-state time series (rep1)
dssp = md.compute_dssp(rep1, simplified=False); dssp[dssp == " "] = "C"
codes_order = ["C", "E", "B", "S", "T", "H", "I", "G"]
code_to_int = {c: i for i, c in enumerate(codes_order)}
dssp_int = np.vectorize(lambda c: code_to_int.get(c, 0))(dssp)
fig, ax = plt.subplots(figsize=(10, 5))
ax.imshow(dssp_int.T, aspect="auto", origin="lower", interpolation="nearest",
          extent=[rep1.time[0] / 1000.0, rep1.time[-1] / 1000.0, 0, rep1.n_residues])
ax.set(xlabel="Time (ns)", ylabel="Residue", title="DSSP 8-state (rep1)")
fig.tight_layout(); fig.savefig(outdir / "07_dssp_8state.png"); plt.close(fig)

# Dump DSSP code mapping so downstream code can pick its own legend/palette
np.savez(outdir / "07_dssp_8state.npz", dssp=dssp, dssp_int=dssp_int, codes_order=codes_order)
print("All figures written to figures/")

工作流 2：跨重复模拟的氢键持久性和 SASA

目标：计算持久氢键（≥50% 的帧中存在）以及各重复模拟的总 SASA 时间序列。

import mdtraj as md, numpy as np, pandas as pd

replicas = {"rep1": "rep1.xtc", "rep2": "rep2.xtc", "rep3": "rep3.xtc"}
top = "protein.pdb"
rows = []
for name, p in replicas.items():
    t = md.load(p, top=top)
    hbonds = md.baker_hubbard(t, freq=0.5, periodic=False)
    sasa_ang2 = (md.shrake_rupley(t, mode="residue").sum(axis=1) * 100.0)
    rows.append({"replica": name, "n_persistent_hbonds": len(hbonds),
                 "sasa_mean_ang2": sasa_ang2.mean(), "sasa_std_ang2": sasa_ang2.std()})
print(pd.DataFrame(rows).to_string(index=False))

关键参数

参数	模块	默认值	范围 / 选项	作用
chunk (md.iterload)	I/O	100	100-5000	每个流式数据块的帧数；在速度与 RAM 之间权衡
atom_indices (md.rmsd)	稳定性	None	整数数组	限定叠合与 RMSD 的计算范围（例如仅 CA）
frame (md.rmsd)	稳定性	0	0-n_frames-1	参考帧
contacts (md.compute_contacts)	接触	"all"	"all"、[(i,j),...]	所有残基对或自定义列表
scheme (md.compute_contacts)	接触	"closest-heavy"	"ca"、"closest"、"closest-heavy"、"sidechain"、"sidechain-heavy"	残基间距离计算所用的原子子集
cutoff（接触阈值，由用户设置）	接触	0.5 nm	0.4-0.8 nm	接触距离（CA 为 5 Å，closest-heavy 为 4 Å）
simplified (md.compute_dssp)	DSSP	True	True/False	False → 8 态代码
freq (md.baker_hubbard)	氢键	0.1	0.0-1.0	报告所需的最小帧比例
mode (md.shrake_rupley)	SASA	"atom"	"atom"、"residue"	按原子计算或按残基求和
probe_radius (md.shrake_rupley)	SASA	0.14 nm	0.12-0.20 nm	溶剂探针半径（标准值为 1.4 Å）
periodic（几何函数）	几何	True	bool	存在盒子信息时采用最小镜像约定

最佳实践

1. 绘图前先转换单位 — nm → Å（×10）、ps → ns（÷1000）、rad → deg（np.degrees）。这是“我的 RMSD 小了 10 倍”的首要原因。
2. 计算 RMSF 前先叠合 — md.rmsd 会在内部叠合；RMSF 不会。若未显式调用 traj.superpose(...)，RMSF 会包含刚体运动。
3. 大轨迹使用 md.iterload — 1 µs、每帧间隔 1 ps 的轨迹包含数百万帧。流式读取可避免内存溢出。
4. 精细分析 DSSP 时使用 simplified=False — 三态模式会将 3_10 螺旋、π 螺旋和弯曲转变合并。
5. 使用配对索引，而非 resSeq — md.compute_contacts 返回从 0 开始的连续索引。仅在标注时映射到 resSeq。
6. 尽早去除水分子 — 对溶剂化体系使用 traj.atom_slice(traj.topology.select("protein"))，可将内存占用和耗时减少 5–20 倍。

常用操作示例

操作示例：转换轨迹格式

import mdtraj as md

traj = md.load("input.nc", top="topology.parm7")
traj.save_xtc("output.xtc")
traj.save_pdb("output.pdb")
print(f"Converted {traj.n_frames} frames -> output.xtc")

操作示例：计算重复模拟中逐残基 RMSF 的均值 ± 标准差

为考虑柔性的对接或同源精修提供集成柔性输入。

import mdtraj as md, numpy as np, pandas as pd

paths = {"rep1": "rep1.xtc", "rep2": "rep2.xtc", "rep3": "rep3.xtc"}
top = "protein.pdb"
rmsf_by_rep = {}
for name, p in paths.items():
    t  = md.load(p, top=top)
    ct = t.atom_slice(t.topology.select("name CA")); ct.superpose(ct, frame=0)
    rmsf_by_rep[name] = np.sqrt(((ct.xyz - ct.xyz.mean(axis=0)) ** 2).sum(axis=2).mean(axis=0)) * 10.0

df = pd.DataFrame(rmsf_by_rep)
df["mean"], df["std"] = df.mean(axis=1), df.std(axis=1)
df.to_csv("rmsf_replicates.csv", index_label="residue_index")

操作示例：保存 RMSD 超过阈值的帧

筛选“开放”构象，用于聚类或后续分析。

import mdtraj as md

traj = md.load("rep1.xtc", top="protein.pdb")
ca   = traj.topology.select("name CA")
traj.superpose(traj, frame=0, atom_indices=ca)
rmsd_ang = md.rmsd(traj, traj, frame=0, atom_indices=ca) * 10.0

open_state = traj[rmsd_ang > 4.0]
open_state.save_dcd("open_state.dcd")
print(f"Saved {open_state.n_frames}/{traj.n_frames} frames")

操作示例：比较重复模拟中的螺旋比例

import mdtraj as md, numpy as np

paths = ["rep1.xtc", "rep2.xtc", "rep3.xtc"]
helix_fracs = [(md.compute_dssp(md.load(p, top="protein.pdb"), simplified=False) == "H").mean() for p in paths]
print(f"α-helix per replicate: {[f'{h:.2f}' for h in helix_fracs]}")
print(f"Mean ± std: {np.mean(helix_fracs):.2f} ± {np.std(helix_fracs):.2f}")

故障排除

问题	原因	解决方案
RuntimeError: No topology	XTC/DCD/TRR 不包含拓扑信息	在 md.load 中传入 top="protein.pdb"（或 .prmtop）
RMSD 小了 10 倍	mdtraj 返回的单位是 nm	乘以 10.0
RMSF 在末端出现峰值	末端本身具有柔性，并非错误	去除末端片段，或明确报告
DSSP 仅返回 H/E/C	默认值为 simplified=True	传入 simplified=False
DSSP 单元格为空或为 ' '	空格表示卷曲，而非缺失	dssp[dssp == " "] = "C"
大轨迹触发 MemoryError	整条轨迹都载入了 RAM	改用 md.iterload(..., chunk=500)
phi/psi 数组的残基数不同	第一个残基没有 phi；最后一个残基没有 psi	配对前求残基索引集合的交集
Ramachandran 图偏离中心	角度仍以弧度表示	使用 np.degrees 转换；设置 range=[[-180, 180], [-180, 180]]
md.compute_contacts 返回 NaN	某些残基对缺少用于 closest-heavy 的重原子（例如甘氨酸侧链）	使用 scheme="closest" 或 "ca"；删除 NaN 列
接触图不对称	仅填充了上三角	镜像：freq_map[j, i] = freq_map[i, j]
traj.time 全为零	某些写入器未填写时间	重建：traj.time = np.arange(traj.n_frames) * dt_ps

相关技能

• mdanalysis-trajectory — 更丰富的原子选择语法和 AnalysisBase 框架
• autodock-vina-docking / smina-molecular-docking — 将代表性 MD 帧作为对接目标
• chembl-database-bioactivity — 在 MD 得出的柔性图谱上叠加配体生物活性

参考资料

• mdtraj 文档 — 官方文档、API 参考和教程
• McGibbon 等人（2015） — “MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories”，Biophys J
• DSSP 算法 — Kabsch 和 Sander（1983）
• Baker–Hubbard 氢键判据 — Baker 和 Hubbard（1984）
• Shrake–Rupley SASA — Shrake 和 Rupley（1973）
• mdtraj GitHub — 源代码、问题和示例