prody

ProDy - タンパク質ダイナミクスと構造生物学

ProDy は、タンパク質の集団運動をモデル化するために設計されています。タンパク質を弾性ネットワークとして扱うことで、単一の PDB ファイルまたは構造アンサンブルから、機能的な運動や構造の柔軟性を予測できます。

使用する場面

• タンパク質の柔軟性と集団運動の予測（ANM/GNM）。
• 構造アンサンブルまたは MD 軌跡に対する主成分分析（PCA）の実施。
• 構造保存性と共進化の解析（Evol）。
• 複数のタンパク質構造の比較（アンサンブル解析）。
• タンパク質のヒンジ領域と剛体ドメインの特定。
• ドッキングの準備と結合部位の解析（ドラッガビリティ）。
• 集団モードに基づく MD 軌跡のフィルタリング。

リファレンスドキュメント

公式ドキュメント: http://prody.csb.pitt.edu/
マニュアル: http://prody.csb.pitt.edu/manual/
検索パターン: prody.parsePDB、prody.ANM、prody.GNM、prody.select、prody.Ensemble

基本原則

原子選択の代数

ProDy は、VMD や PyMOL に似た強力な選択言語を備えています。鎖、残基、属性、または近接性によって原子を選択できます（例：'protein and resname TRP and within 5 of resname HEM'）。

弾性ネットワークモデル（ENM）

• GNM（ガウスネットワークモデル）：揺らぎの大きさ（B 因子）を予測します。
• ANM（異方性ネットワークモデル）：運動の方向と大きさを予測します。

アンサンブル

複数の NMR モデルや MD フレームなどの構造を、統計解析や PCA を迅速に実行できる形式で格納した集合です。

クイックリファレンス

インストール

pip install prody

標準的なインポート

import numpy as np
from prody import *
# Optional: for plotting
# confProDy(auto_show=False)

基本的なパターン - 基準振動解析

from prody import *

# 1. Parse structure
atoms = parsePDB('1p38')
calphas = atoms.select('protein and calpha')

# 2. Build and solve ANM
anm = ANM('p38_anm')
anm.buildHessian(calphas)
anm.calcModes(n_modes=20)

# 3. Analyze results
for mode in anm[:3]:
    print(f"Mode {mode.getIndex()}: Variance = {mode.getVariance():.2f}")

# 4. Save for visualization (NMD format for VMD/PyMOL)
writeNMD('p38_modes.nmd', anm, calphas)

重要なルール

✅ 推奨事項

• NMA では Cα を選択する - 大規模なシステムでは、Cα 原子のみに適用することで、ENM（ANM/GNM）の効果が最も高くなり、計算効率も向上します。
• アンサンブルは必ずアラインメントする - 構造アンサンブルに対して PCA を実行する前に、ensemble.iterpose() を使ってすべてのフレームが参照構造にアラインメントされていることを確認してください。
• 早い段階で select() を使う - ヘッセ行列の計算時にメモリを節約するため、PDB オブジェクトから必要な鎖や原子のみを選択してください。
• 固有値ソルバーの収束を確認する - 計算されたモードが分散の大部分を表していることを確認してください。
• 原子の順序を保持する - 構造を比較する際は、matchAlign() を使い、原子選択のインデックスが一致するようにしてください。

❌ 禁止事項

• 未処理の PDB に NMA を実行しない - PDB にはループの欠損や複数の占有率が含まれていることがよくあります。解析前にデータを整理するか、特定の鎖を選択してください。
• 「ゼロモード」を無視しない - ANM の最初の 6 つのモードは剛体の並進・回転であり、周波数はゼロです。生物学的な実際の運動はモードインデックス 6 から始まります。
• 大きなタンパク質の全原子ヘッセ行列を計算しない - 全原子 ENM は 3N×3N の行列を生成します。残基数が 1000 のタンパク質では 30,000×30,000 の行列となり、大量のメモリを消費します。

避けるべきパターン（絶対に行わない）

from prody import *

# ❌ BAD: Iterating over atoms to find distance
# for a1 in atoms:
#     for a2 in atoms: ... # O(N^2) Python loop

# ✅ GOOD: Use selection algebra
nearby = atoms.select('within 5 of resname LIG')

# ❌ BAD: PCA on unaligned frames
# pca = PCA('test'); pca.buildCovariance(coord_array) # Wrong!

# ✅ GOOD: Create Ensemble and interpose
ens = Ensemble(atoms)
ens.addCoordset(trajectory)
ens.iterpose() # Crucial step
pca = PCA('test')
pca.buildCovariance(ens)

# ❌ BAD: Using NMA modes 0-5 for biology
# slow_mode = anm[0] # This is just a translation/rotation

原子の選択と操作

強力なクエリ

atoms = parsePDB('3hhr')

# Chain and residue range
heavy_chain = atoms.select('chain H and resnum 1 to 120')

# Chemical properties
backbone = atoms.select('backbone')
hydrophobic = atoms.select('resname ALA VAL ILE LEU MET PHE TYR TRP')

# Proximity (Binding site)
site = atoms.select('protein and within 10 of resname ATP')

# Geometric center
center = calcCenter(site)

弾性ネットワークモデル（ENM）

GNM（揺らぎ）

gnm = GNM('1p38_gnm')
gnm.buildKirchhoff(calphas, cutoff=10.0)
gnm.calcModes()

# Cross-correlations (How atoms move together)
cross_corr = calcCrossCorr(gnm)

# Square fluctuations (Theoretical B-factors)
sq_flucts = calcSqFlucts(gnm)

ANM（運動の方向）

anm = ANM('1p38_anm')
anm.buildHessian(calphas, cutoff=15.0)
anm.calcModes()

# Getting the hinge regions (where motion changes direction)
hinges = findHinges(anm[0]) # From the slowest mode

アンサンブル解析と PCA

構造比較

# Parse multiple structures
pdb_ids = ['1p38', '1zz2', '1ywr']
structures = [parsePDB(pid) for pid in pdb_ids]

# Align and match
ensemble = Ensemble('p38_set')
for s in structures:
    # Match calphas of s to the reference first structure
    mappings = matchAlign(s, structures[0])
    ensemble.addCoordset(mappings[0][0]) # Add matched coords

# PCA
pca = PCA('p38_pca')
pca.buildCovariance(ensemble)
pca.calcModes()

# Project structures onto PCs
projection = ensemble.getProjection(pca[:2])

進化解析（Evol）

配列保存性と共進化

# Load Multiple Sequence Alignment (MSA)
msa = parseMSA('p38_alignment.fasta')

# Calculate conservation (Shannon Entropy)
entropy = calcShannonEntropy(msa)

# Mutual Information (Co-evolution)
mi = calcMutualInformation(msa)

# Direct Coupling Analysis (DCA) - requires external tools or specific plugins
# dca = calcDirectCovariance(msa)

実践的なワークフロー

1. 機能的な「ヒンジ」の特定

def get_protein_hinges(pdb_id):
    atoms = parsePDB(pdb_id)
    calphas = atoms.select('protein and calpha')
    
    anm = ANM(pdb_id)
    anm.buildHessian(calphas)
    anm.calcModes()
    
    # Hinge residues for the first two functional modes
    hinges_m1 = findHinges(anm[0])
    hinges_m2 = findHinges(anm[1])
    
    return list(set(hinges_m1) | set(hinges_m2))

2. MD 軌跡と ANM モードの比較

def compare_md_to_anm(md_traj, p_pdb):
    # 1. ANM from static structure
    atoms = parsePDB(p_pdb)
    anm = ANM('static')
    anm.buildHessian(atoms.select('calpha'))
    anm.calcModes()
    
    # 2. PCA from MD
    ens = Ensemble('md')
    ens.setCoords(atoms)
    ens.addCoordset(md_traj)
    ens.iterpose()
    pca = PCA('md_pca')
    pca.buildCovariance(ens)
    pca.calcModes()
    
    # 3. Overlap (Inner product of modes)
    overlap = calcOverlap(anm[0], pca[0])
    return overlap

3. ドラッガビリティ解析（TRAWLER/ProDy の統合）

# Note: Full druggability analysis usually involves 'hotspot' calculations
def binding_site_flexibility(atoms, lig_resname):
    site = atoms.select(f'protein and within 8 of resname {lig_resname}')
    # Calculate GNM just for the site context
    gnm = GNM('site')
    gnm.buildKirchhoff(atoms.select('calpha'))
    gnm.calcModes()
    
    # High fluctuations = likely flexible binding site
    flucts = calcSqFlucts(gnm)
    return flucts[site.getIndices()]

パフォーマンスの最適化

大きなヘッセ行列のメモリ管理

非常に大きな複合体（リボソーム、カプシド）では、利用可能な場合、疎行列または階層ネットワークモデル（HNM）を使用してください。

並列 MSA パース

大規模なアラインメント（100k+ 配列）を扱う場合は、メモリ効率に優れた特定のフラグを指定して parseMSA を使用してください。

よくある落とし穴と解決策

原子マッピングの不一致

同じタンパク質の 2 つの PDB を比較する際、一方に残基が欠けている場合があります。

# ❌ Problem: ensemble.addCoordset(pdb2) fails due to different atom counts
# ✅ Solution: Use matchAlign
matches = matchAlign(pdb2, pdb1)
if matches:
    ensemble.addCoordset(matches[0][0])

標準外の残基

ProDy は、通常とは異なるリガンドを「ヘテロ」や「タンパク質」として認識しない場合があります。

# ✅ Solution: Use specific residue names or 'all'
ligand = atoms.select('resname MYL')

ヘッセ行列の特異性

タンパク質に分離した部分がある場合、ヘッセ行列には 6 個を超えるゼロ固有値が存在します。

# ✅ Solution: Check connectivity
if not atoms.select('protein').connected:
    print("Warning: Disconnected components found!")

ベストプラクティス

1. 大きなタンパク質で NMA を行う場合は、計算コストを削減するため、常に Cα 原子を選択してください。
2. PCA を行う前に、iterpose() を使ってアンサンブル内のすべての構造をアラインメントしてください。
3. 選択代数を使って PDB 構造を早い段階でフィルタリングし、メモリ使用量を削減してください。
4. ANM のモード 0～5 は剛体運動であり、生物学的な運動はモード 6 から始まることを忘れないでください。
5. 原子数が異なる構造や残基が欠落している構造を比較する場合は、matchAlign() を使用してください。
6. ヘッセ行列を構築する前に、分離した成分がないか確認してください。
7. 適切なカットオフ距離（Cα ネットワークでは通常 10～15 Å）を使用してください。
8. 意味のある結果を得るため、固有値ソルバーの収束を検証してください。
9. VMD または PyMOL で可視化できるよう、モードを NMD 形式で保存してください。
10. 非常に大きなシステムでは、疎行列による表現の使用を検討してください。

ProDy は、構造生物学における「ダイナミクス」を扱うために欠かせないツールキットです。タンパク質を物理ネットワークとして捉えることで、静的なスナップショットと、分子スケールで絶えず振動する生命の現実をつなぎます。