ProDy — 蛋白质动力学与结构分析 Python 库¶

一句话说明¶

ProDy 是用 Python 分析蛋白质结构动力学的工具包，可以读写 PDB/DCD 轨迹文件、做正态模分析（NMA）、主成分分析（PCA）、弹性网络模型（ANM/GNM），理解蛋白质运动模式。白话理解：用 Python 分析蛋白质"怎么抖动"的工具，不用跑 MD 就能预测蛋白质运动。

安装与配置¶

# pip 安装（推荐）
pip install prody  # 安装 ProDy（当前 v2.4+）

# conda 安装
conda install -c conda-forge prody

# 安装可视化扩展（用于 VMD/NMWiz 接口）
pip install prody[nmwiz]   # 含 NMWiz 可视化支持

# 确认安装
python -c "import prody; print(prody.__version__)"  # 查看版本

# 常用配套工具
pip install matplotlib numpy scipy  # 绘图和数值计算

核心用法¶

读取蛋白质结构¶

import prody as pd_prody  # 导入 ProDy（避免与 pandas 的 pd 混淆）

# 从 PDB 数据库下载结构
protein = pd_prody.parsePDB("1ake")   # 下载并解析 PDB 1AKE（腺苷酸激酶）
print(f"原子数: {protein.numAtoms()}")        # 总原子数
print(f"残基数: {protein.numResidues()}")     # 总残基数
print(f"链: {list(protein.getChids())}")      # 链标识符列表

# 从本地文件读取
protein_local = pd_prody.parsePDB("my_protein.pdb")  # 读取本地 PDB

# 选择原子（用 CHARMM 选择语法）
calpha = protein.select("calpha")              # 只选 Cα 原子
chain_a = protein.select("chain A")           # 只选 A 链
active_site = protein.select("resnum 45 85 120 and calpha")  # 活性位点 Cα

弹性网络模型（GNM/ANM）— 预测蛋白质运动¶

import prody as pd_prody
import matplotlib.pyplot as plt

# 加载腺苷酸激酶结构
protein = pd_prody.parsePDB("1ake")
calpha = protein.select("calpha and chain A")  # 选 A 链的 Cα 原子

# === 高斯网络模型（GNM）===
# GNM：各向同性运动，预测哪些残基"动得多"（柔性区域）
gnm = pd_prody.GNM("Adenylate Kinase")        # 创建 GNM 对象
gnm.buildKirchhoff(calpha, cutoff=7.5)         # 构建连接矩阵（7.5Å截断）
gnm.calcModes(n_modes=10)                      # 计算前10个低频运动模式

# 可视化 GNM 慢模式（最大运动方向）
pd_prody.showSqFlucts(gnm[:3])    # 显示前3个模式的均方根波动
plt.title("GNM Slow Modes")
plt.savefig("gnm_sqflucts.pdf")

# === 各向异性网络模型（ANM）===
# ANM：各向异性运动，预测运动方向（三维向量）
anm = pd_prody.ANM("Adenylate Kinase")
anm.buildHessian(calpha, cutoff=13.0)          # 构建 Hessian 矩阵（13Å截断）
anm.calcModes(n_modes=20)                      # 计算前20个模式

# 查看最慢模式（最重要的整体运动）
mode1 = anm[0]                                 # 第一个（最慢）模式
print(f"模式1 频率: {mode1.getEigval():.4f}")  # 本征值（越小=运动越慢=越重要）

参数详解¶

参数/方法	含义	常用值
`GNM.buildKirchhoff(cutoff)`	GNM 接触截断距离（Å）	7-8 Å
`ANM.buildHessian(cutoff)`	ANM 接触截断距离（Å）	10-15 Å
`calcModes(n_modes)`	计算模式数量	10-50
`showSqFlucts()`	显示均方根波动	找柔性区域
`showCrossCorr()`	显示残基交叉相关矩阵	找协同运动残基
`calcCollectivity()`	计算集体性指数（0-1）	1=整体运动，0=局部运动

实战案例¶

案例：分析 MD 轨迹做 PCA（主成分分析）¶

import prody as pd_prody
import matplotlib.pyplot as plt
import numpy as np

# 读取 MD 轨迹（DCD 格式，GROMACS 导出）
# 先把 xtc 转为 dcd：gmx trjconv -f md.xtc -o md.dcd

# 建立结构集合（ensemble）
ref_structure = pd_prody.parsePDB("initial.pdb")   # 参考结构
calpha_ref = ref_structure.select("calpha")         # 只要 Cα

# 读取轨迹
ensemble = pd_prody.PDBEnsemble("MD Trajectory")
traj = pd_prody.Trajectory("md.dcd")              # 读取 DCD 轨迹
traj.link(ref_structure)                           # 链接参考结构
traj.setAtoms(calpha_ref)                          # 只分析 Cα 原子

# 逐帧读取并叠合到参考结构
for frame in traj:
    ensemble.addCoordset(frame.getCoords())        # 添加每帧坐标

ensemble.setCoords(calpha_ref.getCoords())         # 设置参考坐标
ensemble.superpose()                               # 叠合所有帧（去掉刚体运动）

# 做 PCA（主成分分析）
pca = pd_prody.PCA("MD PCA")
pca.buildCovariance(ensemble)                      # 计算协方差矩阵
pca.calcModes(n_modes=10)                          # 计算前10个主成分

# 把轨迹投影到 PC1 和 PC2
proj = pd_prody.calcProjection(ensemble, pca[:2])  # 投影到前2个 PC

# 绘制 PCA 散点图
plt.figure(figsize=(8, 6))
plt.scatter(proj[:, 0], proj[:, 1],               # PC1 vs PC2 散点图
            c=range(len(proj)),                    # 按帧数着色（时间序列）
            cmap="RdYlBu_r",                       # 蓝→红表示时间进展
            s=20, alpha=0.7)
plt.colorbar(label="Frame Number")                 # 颜色条标签
plt.xlabel(f"PC1 ({pca[0].getVariance()*100:.1f}%)")  # PC1 方差贡献
plt.ylabel(f"PC2 ({pca[1].getVariance()*100:.1f}%)")  # PC2 方差贡献
plt.title("MD Trajectory PCA")
plt.savefig("pca_trajectory.pdf")
print(f"PC1 解释方差: {pca[0].getVariance()*100:.1f}%")
print(f"PC2 解释方差: {pca[1].getVariance()*100:.1f}%")

案例：计算蛋白质间动态相关性¶

import prody as pd_prody
import matplotlib.pyplot as plt

# 加载蛋白质
protein = pd_prody.parsePDB("1ake")
calpha = protein.select("calpha and chain A")

# 构建 ANM 模型
anm = pd_prody.ANM()
anm.buildHessian(calpha, cutoff=13.0)  # 构建 Hessian
anm.calcModes(20)                       # 计算前20个模式

# 计算交叉相关矩阵（Cross-Correlation Matrix）
# 正值=正相关（同向运动），负值=负相关（反向运动）
cross_corr = pd_prody.calcCrossCorr(anm[:20])   # 用前20个模式

# 可视化
fig, ax = plt.subplots(figsize=(8, 6))
img = ax.imshow(cross_corr,                      # 显示矩阵为热图
                cmap="RdBu_r",                   # 红=正相关，蓝=负相关
                vmin=-1, vmax=1)                 # 颜色范围 -1 到 1
plt.colorbar(img, label="Cross-Correlation")    # 颜色条
ax.set_xlabel("Residue Index")
ax.set_ylabel("Residue Index")
ax.set_title("ANM Cross-Correlation Matrix")
plt.tight_layout()
plt.savefig("cross_correlation.pdf")

常见报错与解决¶

报错	原因	解决方法
`No atoms match the selection`	选择语法错误	检查链名、残基编号、原子名拼写
`Trajectory file not found`	DCD 文件路径错误	使用绝对路径；确认文件格式为 DCD
`Hessian matrix is singular`	体系刚性不足或截断距离太小	增大 `cutoff` 值（ANM 用 13-15Å）
PCA 结果不合理	叠合未去旋转平移	确认执行了 `ensemble.superpose()`
内存不足	轨迹帧数太多	用 `stride=10` 每10帧读一帧

速查表¶

# 安装
pip install prody matplotlib

# 读取结构
import prody
protein = prody.parsePDB("1ake")   # 网络下载
protein = prody.parsePDB("x.pdb") # 本地文件

# 选择语法（CHARMM 风格）
protein.select("calpha")            # 所有 Cα
protein.select("chain A")           # A 链
protein.select("protein and backbone")  # 骨架原子
protein.select("resnum 50 to 100")  # 残基编号范围

# GNM（各向同性运动预测）
gnm = prody.GNM(); gnm.buildKirchhoff(ca); gnm.calcModes(10)

# ANM（各向异性运动预测，含方向）
anm = prody.ANM(); anm.buildHessian(ca, cutoff=13); anm.calcModes(20)

# 可视化快捷函数
prody.showSqFlucts(gnm[:3])         # 均方根波动
prody.showCrossCorr(anm[:10])       # 交叉相关矩阵
prody.showMechanicalStiffness(anm)  # 机械刚度

# 当前版本：ProDy v2.4+