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（各向异性网络模型）：预测运动的方向和幅度。

集合

以便于快速统计分析和 PCA 的方式存储的一组结构（例如，多个 NMR 模型或 MD 帧）。

快速参考

安装

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 对象筛选为仅包含必要的链/原子，以便在计算 Hessian 矩阵时节省内存。
• 检查特征求解器的收敛性 - 确保计算出的模式代表了大部分方差。
• 保留原子顺序 - 比较结构时，确保原子选择产生匹配的索引，使用 matchAlign()。

❌ 不应该做

• 直接对原始 PDB 执行 NMA - PDB 文件通常存在缺失环或多个占有率。分析前应清理数据或选择特定链。
• 忽略“零模式” - ANM 的前 6 个模式是刚体平移/旋转，频率为零。真实的生物运动从模式索引 6 开始。
• 为大型蛋白质计算全原子 Hessian - 全原子 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()]

性能优化

大型 Hessian 的内存管理

对于非常大的复合物（核糖体、衣壳），如果可用，请使用稀疏矩阵或分层网络模型（HNM）。

并行解析 MSA

处理海量比对（10 万条以上序列）时，使用带有特定内存高效标志的 parseMSA。

常见问题与解决方案

原子映射不匹配

比较同一种蛋白质的两个 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')

Hessian 矩阵奇异性

如果蛋白质包含不连通的部分，Hessian 矩阵将有超过 6 个零特征值。

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

最佳实践

1. 对大型蛋白质进行 NMA 时，始终选择 C-alpha 原子，以降低计算开销。
2. 使用 iterpose() 执行 PCA 前，先对齐集成中的所有结构。
3. 尽早使用选择代数筛选 PDB 结构，以减少内存使用。
4. 请记住，ANM 模式 0-5 是刚体运动；生物运动从模式 6 开始。
5. 比较原子数不同或缺少残基的结构时，使用 matchAlign()。
6. 构建 Hessian 矩阵前，检查是否存在不连通的部分。
7. 使用适当的截断距离（C-alpha 网络通常为 10-15 Å）。
8. 验证特征值求解器是否收敛，以确保结果有意义。
9. 将模式保存为 NMD 格式，以便在 VMD 或 PyMOL 中可视化。
10. 对于非常大的系统，可考虑使用稀疏矩阵表示。

ProDy 是结构生物学中“动态”这一环节不可或缺的工具包。它将蛋白质视为物理网络，架起了静态快照与分子尺度上生命振动现实之间的桥梁。