MDAnalysis蛋白轨迹路径比较

两条轨迹经过相近构象，不一定按相同顺序经过。Hausdorff距离比较“到过哪些构象”，离散Fréchet距离还考虑访问顺序。本例使用集群已有泛素轨迹的51帧、每帧76个Cα原子，先统一对齐同一参考结构，再做路径比较。

1 下载、运行和看图#

下载水分子动力学与蛋白路径比较完整示例，保存为water-paths.zip。该包使用MDAnalysis 2.9.0，独立扩展waterdynamics 1.2.0和PathSimAnalysis 1.2.1；两项扩展只安装在示例目录的.deps中。

python3 -m zipfile -e water-paths.zip .
cd mdanalysis-water-paths
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

results/FULL-FIELDS.npz是完整数值，results/VALIDATION.json列出逐项核对，plots/保存4张图片。脚本发现已有同名输出目录会停止，第二次运行可以用bash run.sh --output results-new指定新目录；不要把旧结果当成新轨迹。

compact-inputs/protein-paths.npz保存共同参考坐标。比较4条路径：完整51帧、前26帧、后26帧、完整路径倒序。前后半段共享中间帧；倒序只是相同坐标的顺序对照，不是又跑了一条逆时间MD。

2 直接调用路径距离#

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
import pathsimanalysis as psa
x = np.load('compact-inputs/protein-paths.npz')['positions_A'].astype(float)
print('full vs reversed Hausdorff:', psa.hausdorff(x, x[::-1]))
print('full vs reversed Frechet:', psa.discrete_frechet(x, x[::-1]))
print('first vs second Hausdorff:', psa.hausdorff(x[:26], x[25:]))
print('first vs second Frechet:', psa.discrete_frechet(x[:26], x[25:]))
PYCODE

本例距离单位为Å，每一对构象用全部76个Cα的均方根距离，不再逐对重新最优拟合：

完整与倒序经过同一组构象，所以Hausdorff为0；Fréchet仍可区分访问顺序。数值越大只代表此坐标定义下的几何差异越大，不直接等于能垒、自由能差或动力学速率。

3 多条轨迹用PSAnalysis组织#

配套脚本写出真实Cα坐标的PDB/DCD，并实际运行以下流程：

analysis = psa.PSAnalysis(
    universes, select='all', path_select='all', labels=labels,
    targetdir='results/psa'
)
analysis.generate_paths(align=False, save=True, store=True)
analysis.paths = [p.astype(np.float64) for p in analysis.paths]
analysis.run(metric='hausdorff')
hausdorff_matrix = analysis.D.copy()
analysis.run(metric='discrete_frechet')
frechet_matrix = analysis.D.copy()

这里align=False是因为全部轨迹已用同一参考对齐；不是建议跳过结构对齐。显式转换为双精度只改变计算精度，每个坐标值保持不变，避免Fréchet动态规划继续使用默认单精度。自己的轨迹若未统一对齐，应先选相同原子、相同顺序、相同参考，不能对每一对构象各用不同参考后再宣称距离具有同一几何含义。

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
f = np.load('results/FULL-FIELDS.npz')
print('Hausdorff matrix:')
print(f['path_hausdorff_A'])
print('Frechet matrix:')
print(f['path_discrete_frechet_A'])
np.savetxt('path-hausdorff.dat', f['path_hausdorff_A'])
np.savetxt('path-frechet.dat', f['path_discrete_frechet_A'])
PYCODE

4 平均距离和最近邻帧#

包中还实际运行hausdorff_avg、hausdorff_wavg和hausdorff_neighbors。前两项是最近邻距离的平均变体，会减弱单个极端构象的作用，但不满足真正距离度量的全部性质；不要与取最大值的标准Hausdorff混用。hausdorff_neighbors给出两个方向每个构象对应的最近帧，适合找出两条路径在哪一段分开。

所有4×4矩阵均保存在FULL-FIELDS.npz，最近邻索引和距离也完整保留。这里的51帧来自短时间轨迹，路径不同不等于已发现两个稳定态，更不能据此宣称构象空间采样收敛。

参考：PathSimAnalysis原生API。