MDAnalysis蛋白接触、GNM与PCA

想知道蛋白哪些区域一起运动，可以比较参考接触、Cα网络和坐标主成分。三个工具回答的问题不同：接触是距离规则，GNM是几何网络，PCA是所采样轨迹的坐标方差。

1 下载并运行真实示例#

下载泛素蛋白八类分析输入、轨迹、结果与绘图包，保存为 mdanalysis-protein-wave326-v3.zip：

module load pymatgen/2026.9.24
python -m zipfile -e mdanalysis-protein-wave326-v3.zip .
cd mdanalysis-protein
bash run.sh
source env.sh

run.sh在新目录中完成八类分析并生成results/和plots/，已有同名输出时先换一个解压目录。环境使用集群已安装的Python 3.11 / MDAnalysis 2.9.0；env.sh选择相应Python并限制数值库为单线程，没有调用qvasp或VASP。其他机器可把ANALYSIS_PYTHON设为含相同依赖的Python路径用于run.sh。

本包使用鸢算科技集群已完成的GROMACS泛素轨迹：1231个原子、76个残基、1237条键。原始轨迹有1001帧，间隔0.1 ps；本例每20帧取一帧，得到0、2、…、100 ps共51帧。这里重用真实轨迹进行分析，没有重新做分子动力学；100 ps的变化不能说明蛋白折叠、长时稳定性或构象自由能已经收敛。

2 先处理拓扑和周期性边界#

本例GROMACS 2026.4生成的TPR格式为138，MDAnalysis 2.9.0不能直接读取。包中已提供由原生GROMACS导出的protein-conect.pdb及质量、电荷、元素元数据，不靠PDB猜测来代替拓扑。需要重做导出时：

bash export-topology.sh
source env.sh
python - <<'PYCODE'
from load import load
u = load()
print(len(u.atoms), len(u.residues), len(u.bonds))
print(len(u.trajectory), u.trajectory.dt)
print(len(u.atoms.fragments), u.atoms.charges.sum())
PYCODE

前两行分别得到1231 76 1237及1001 0.1；只有一个连通蛋白片段，总电荷约为数值零。export-topology.sh实际使用：

module load gromacs/2026.4-p100-threadmpi
# export-topology.sh已建立topology-export目录
gmx editconf -f inputs/protein.tpr -o topology-export/protein-conect.pdb -conect
gmx dump -s inputs/protein.tpr > topology-export/native-topology.dump
python topology_metadata.py

重新导出的文件在topology-export/，不会覆盖已有inputs/。topology_metadata.py从原生TPR转储提取质量、电荷和原子序号；本例脚本限定为这个单蛋白拓扑，不要直接当成任意多分子溶剂体系的通用转换器。load()随后加载PDB和XTC、补齐明确元数据，并用键连接执行unwrap让蛋白完整。自己的轨迹必须先确认原子顺序、残基命名、键和盒信息匹配；配体、断链或多链体系要相应修改选择。

3 Contacts：相对于第一帧的接触比例#

python - <<'PYCODE'
from load import load
from MDAnalysis.analysis.contacts import Contacts
u = load()
u.trajectory[0]
ca = u.select_atoms('name CA')
a = Contacts(u, select=('name CA','name CA'), refgroup=(ca,ca),
             method='radius_cut', radius=8.0, pbc=False)
a.run(stop=1001, step=20)
print(a.results.timeseries.shape)
print(a.results.timeseries[:,1].min(), a.results.timeseries[:,1].max())
PYCODE

本例使用完整展开的蛋白，8 Å内的第一帧Cα对作为参考，输出(51,2)，两列为帧索引和接触比例。实际范围为 0.9695291–1.0000000。原生同组选取包含自身及序列局部对，因此它不是排除了局部接触的折叠Q值，也不是“约97%蛋白已折叠”。研究折叠时应另行定义非局部残基对和距离函数。

4 GNM：从7 Å邻接网络读柔性模式#

python - <<'PYCODE'
from load import load
from MDAnalysis.analysis.gnm import GNMAnalysis
u = load()
a = GNMAnalysis(u, select='name CA', cutoff=7.0).run(stop=1001, step=20)
print(a.results.eigenvalues.shape)
print(a.results.eigenvalues.min(), a.results.eigenvalues.max())
print(a.results.eigenvectors.shape)
PYCODE

每帧76个Cα构建Kirchhoff网络，原生输出最低非零本征值及对应模式。51帧本征值范围约 0.157966–0.378247。这是几何网络本征值，没有赋予真实质量、力常数或时间单位，不能直接标成THz声子频率。

本征矢可以整体翻号；简并模式也可以在子空间内旋转。比较两个结果时应比较本征值、子空间或接触网络，不能仅因一条模式曲线正负颠倒就认为计算错误。

5 PCA：先对齐，再看方差#

python - <<'PYCODE'
from load import load
from MDAnalysis.analysis.pca import PCA
u = load()
a = PCA(u, select='name CA', align=True).run(stop=1001, step=20)
print(a.results.variance.shape)
print(a.results.cumulated_variance[2])
PYCODE

76个Cα有228个坐标维度，本例输出228项方差。前三个主成分合计解释 44.7501% 的坐标方差，剩余约55.25%仍在其他方向，不能只看二维散点就认为捕捉了全部运动。PCA官方手册介绍了协方差和累计方差。

本包analyze.py对每帧做一致的Kabsch对齐、减去该采样的平均坐标，再投影到原生主轴。新轨迹投影也必须使用一致的原子顺序、对齐参照及中心化方式；不要对未对齐的原始轨迹直接调用投影并把整体转动当内部运动。

python - <<'PYCODE'
import numpy as np
z = np.load('results/FULL-FIELDS.npz')
print(z['PCA_projection'].shape)
print(z['PCA_variance'][:3])
print(z['Contacts_Q'][:5])
PYCODE

PCA_projection形状为(51,3)。主轴正负本身任意，不同软件图形镜像仍可能表示同一子空间。这些是100 ps样本的运动描述，不能直接给出状态间自由能或转变速率。二级结构与二面角见上一篇，内坐标与螺旋见下一篇。

6 逐帧结构距离矩阵#

需要直接比较任意两个时刻的结构差异时，可以使用原生DistanceMatrix。下载水、聚合物和蛋白距离矩阵扩展包，按前述方式在新目录解压、运行bash run.sh，再读取实际51帧Cα共同对齐坐标：

source env.sh
python - <<'PYCODE'
import numpy as np
import MDAnalysis as mda
from MDAnalysis.analysis.diffusionmap import DistanceMatrix
x = np.load('results/FULL-FIELDS.npz')['protein_aligned_CA_A']
u = mda.Universe.empty(x.shape[1], trajectory=True)
u.load_new(x, order='fac', dt=2.0)
a = DistanceMatrix(u, select='all', cutoff=0).run()
print(a.results.dist_matrix.shape, a.results.dist_matrix.max())
np.savetxt('my-protein-distance.csv', a.results.dist_matrix, delimiter=',')
PYCODE

这里得到 51×51 矩阵，最大距离约 1.199687 Å。每一帧已按同一第一帧参照对齐；矩阵随后比较这些共同对齐的坐标，不再为每一对额外做独立最优叠合。两种距离定义可能不同，自己的计算要明确选择。

图的横纵轴均为0–100 ps，主对角线为零，矩阵对称。颜色显示这份短轨迹采样到的结构差异，不代表构象间自由能、转变速率或扩散系数。DistanceMatrix位于diffusionmap模块里，但只生成距离矩阵并不等于已经完成扩散映射或状态动力学研究。