mdanalysis-trajectory

MDAnalysis — 分子動力学トラジェクトリ解析

概要

MDAnalysis は、使用する MD エンジン（GROMACS、AMBER、NAMD、CHARMM、LAMMPS、OpenMM）にかかわらず、分子動力学トラジェクトリを読み込み、解析するための統一された Python インターフェースを提供します。位置、速度、力、トポロジーデータを持つ AtomGroup を含む Universe オブジェクトとして分子系を表現します。トラジェクトリはフレームごとに反復処理することも、RMSD、RMSF、回転半径、水素結合、溶媒アクセス可能表面積、PCA の解析モジュールを使って一括解析することもできます。MDAnalysis は NumPy、pandas、matplotlib と統合されており、計算化学や創薬におけるシミュレーション後の構造解析の標準ツールです。

使用する場面

• MD シミュレーション後に、タンパク質の主鎖または特定の残基グループの RMSD と RMSF を計算する
• リガンド結合の安定性を解析する：ポケットの RMSD、接触の持続性、水素結合の占有率
• トラジェクトリのコンフォメーションに対して主成分分析（PCA）を行う
• 溶媒アクセス可能表面積（SASA）、回転半径、両端間距離を計算する
• 可視化のために、長時間の MD トラジェクトリから代表的なクラスタ構造を抽出する
• HPC パイプライン内でエンジン固有の解析を行う場合は、代わりに GROMACS または AMBER の解析ツール（gmx rms、cpptraj）を使用する
• MD シミュレーションの実行には 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 分の 1 になる
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	投影する PC の数
groupselection	RMSD	None	フィッティング用の二次選択。RMSD の計算には一次選択を使用する
weights	RMSD/align	None	質量加重 RMSD には "mass" を指定する
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：2 つのタンパク質ドメイン間の接触マップを計算する

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 et al. (2011) 「MDAnalysis: A toolkit for the analysis of molecular dynamics simulations」 — J Comput Chem 32:2319-2327. DOI:10.1002/jcc.21787
• Gowers RJ et al. (2016) 「MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations」 — Proc SciPy 2016. DOI:10.25080/majora-629e541a-00e