MDAnalysis主成分与构象聚类

本篇用已有泛素轨迹的76个Cα做降维和几何聚类。取前26帧与后26帧，共52行坐标，其中中间一帧重复一次；全部结构已对齐同一参考。这个重复是两段定义的结果，必须在解释聚类时保留说明。

1 运行完整示例#

下载降维、聚类与分段采样检查示例，保存为ensemble-convergence.zip：

python3 -m zipfile -e ensemble-convergence.zip .
cd mdanalysis-ensemble-convergence
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

环境为MDAnalysis 2.9.0、NumPy 1.26.4、SciPy 1.13.1、scikit-learn 1.9.1；setup.sh在本例目录读取已有依赖，不修改共享安装。数值在results/FULL-FIELDS.npz，4张图在plots/。重复分析可以用bash run.sh --output results-new指定新目录。

2 原生PCA如何计算#

method = PrincipalComponentAnalysis(
    dimension=2, svd_solver='full', whiten=False, random_state=334
)
coordinates, details = method(cartesian_coordinates.astype(float))

cartesian_coordinates为52×228数组（76×3），不是52×52距离矩阵。配套脚本在原生接口入口显式用双精度做完整SVD，坐标本身没有变化；这样重复的同一构象不会因单精度SVD舍入出现人为差别。whiten=False保留各主成分原本的幅度。

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
f = np.load('results/FULL-FIELDS.npz')
print('PCA coordinate shape:', f['PCA_2D_coordinates'].shape)
print('first two eigenvalues:', f['PCA_eigenvalues_A2'])
np.savetxt('protein-pca.dat', f['PCA_2D_coordinates'].T)
PYCODE

前两个主成分分别解释约25.19%和10.75%的坐标方差，合计35.94%。因此二维图不能代表全部运动，更不能仅凭点云的两个团块断言出现两个稳定态。主成分的正负方向可以翻转，比较不同运行时应核对特征值和子空间，而非图形朝向。

3 三种聚类方法输入不同#

kmeans = KMeans(n_clusters=3, n_init=10, max_iter=300, random_state=334)
kmeans(cartesian_coordinates.astype(float))
dbscan = DBSCAN(eps=0.35, min_samples=2)
dbscan(common_reference_rmsd_matrix)
affinity = AffinityPropagation(
    damping=0.9, preference=-0.5, max_iter=500,
    convergence_iter=50, random_state=334
)
affinity(common_reference_rmsd_matrix)

这里的后两项矩阵使用ENCORE的TriangularMatrix对象，距离单位为Å。KMeans直接用全部228维坐标；DBSCAN与AffinityPropagation用76个Cα的共同参考RMSD。本例没有对每一对构象额外重新叠合。

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
f = np.load('results/FULL-FIELDS.npz')
for key in ['KMeans_labels', 'DBSCAN_raw_labels', 'Affinity_labels']:
    label = f[key]
    print(key, dict(zip(*np.unique(label, return_counts=True))))
np.savetxt('cluster-labels.dat', np.c_[f['KMeans_labels'],
           f['DBSCAN_raw_labels'], f['Affinity_labels']], fmt='%d',
           header='KMeans DBSCAN AffinityPropagation')
PYCODE

KMeans被要求分为3组，不能用“得到3组”证明体系天然有3个状态。DBSCAN在本例eps=0.35 Å条件下给出50个噪声点；唯一的双成员组来自那两个完全相同的中间帧，不能称为新发现的亚稳态。读取DBSCAN_raw_labels时，−1表示噪声；不要把代表帧编号或噪声编号当作物理状态标签。AffinityPropagation本次得到34组，也说明小参数和短时间样本会产生很细的几何分组。

应用到自己的轨迹时，先用相同原子顺序统一参考对齐，检查是否混入重复帧，再调整距离尺度、聚类参数和采样间隔。聚类是几何整理工具，动力学状态与自由能解释还需要独立的采样和时间尺度检查。

参考：MDAnalysis 2.9.0 ENCORE聚类接口。