3. 结构比对中的坐标系与变换

结构比对的核心问题可以一句话概括：在不让分子发生形变的前提下，把两个结构“摆”到最合适的位置，再衡量它们差多少。这句话里藏着两个独立的技术点——一是“摆”的过程，即刚体变换；二是“差多少”的度量，即RMSD。二者必须分开理解，否则很容易在实现时把优化目标和评价指标混为一谈。

刚体变换：只允许旋转和平移

刚体变换（rigid-body transformation）指只包含旋转与平移、不改变任意两点间距离的变换。对蛋白结构而言，这意味着键长、键角、二面角全部保持不变，我们只是整体搬动和转动一个结构去贴合另一个。数学上，把结构A的原子坐标 $x_i$ 变换到新位置 $x_i'$ 的公式为：

x_i' = R · x_i + t

其中 $R$ 是 3×3 旋转矩阵，$t$ 是 3×1 平移向量。整个比对问题就变成：找到一组 $(R, t)$，使变换后的A与B尽可能重合。注意这里只有6个自由度——3个旋转角加3个平移分量，这也是为什么刚体比对是“低维优化”，不需要像分子动力学那样处理成千上万个坐标变量。

旋转矩阵：直观但有约束

旋转矩阵 $R$ 必须满足两个条件：正交性 $R^T R = I$，以及行列式 $\det(R) = +1$。第二个条件排除了镜像反射——如果允许 $\det(R) = -1$，你会得到一个“镜像分子”，这在化学上毫无意义，因为天然蛋白质都是L-氨基酸构成的手性分子。这是结构比对中第一个常见陷阱：用最小二乘直接解出的矩阵可能带反射，必须检查行列式符号。

旋转矩阵的优点是直观：它的三个列向量分别代表变换后坐标系的三个基轴方向。缺点是参数冗余——9个元素只有3个独立自由度，且优化时必须维持正交约束，数值上容易漂移。因此在实际算法中，旋转矩阵常作为最终输出形式，而优化过程往往借助其他参数化方式。

四元数：无奇异的旋转参数化

四元数用一个标量 $q_0$ 和三个分量 $(q_1, q_2, q_3)$ 表示旋转，满足单位约束 $q_0^2 + q_1^2 + q_2^2 + q_3^2 = 1$。相比欧拉角，四元数没有万向节死锁；相比旋转矩阵，它只有4个参数、约束更简单，且插值和求导都更平滑。在Kabsch算法的四元数变体中，问题被转化为求一个4×4对称矩阵的最大特征值对应的特征向量，该特征向量就是最优旋转四元数。

需要强调的是：四元数 $q$ 与 $-q$ 表示同一个旋转。这不是bug，而是双重覆盖性质。如果你在比对结果中看到四元数符号翻转，不要惊慌，它对应的物理旋转完全相同。但如果你在做四元数插值或聚类，就必须先统一符号，否则会把同一个姿态当成两个不同结果。

RMSD：定义与计算

RMSD（Root Mean Square Deviation，均方根偏差）的定义为：

RMSD = sqrt( (1/N) · Σ_i |x_i' - y_i|² )

其中 $x_i'$ 是结构A经过最优刚体变换后的第 $i$ 个原子坐标，$y_i$ 是结构B对应原子的坐标，$N$ 是参与计算的原子对数。注意：RMSD是在最优叠合之后计算的，不是随便摆两个结构直接算。先叠合、后算RMSD，这个顺序不能颠倒。

计算流程通常分三步：第一，确定原子对应关系（哪些原子和哪些原子配对）；第二，用Kabsch或四元数方法求最优 $(R, t)$；第三，代入公式计算RMSD。其中第一步往往比后两步更棘手——序列比对、结构对齐、原子匹配策略都会直接影响结果。

计算陷阱与工程实践

陷阱一：反射矩阵。用SVD求解时，若 $\det(R) = -1$，说明最优解包含反射。标准修正方法是把SVD中最小奇异值对应的向量取反，重新构造旋转矩阵。很多开源实现默认不检查这一点，导致手性翻转的隐蔽错误。

陷阱二：原子对应决定上限。RMSD再小，如果配对错了也没有意义。例如比对两个同源蛋白时，若把活性位点残基和表面loop残基错误配对，即使RMSD数值很低，生物学结论也是错的。实践中应先做序列比对确定等价残基，再只用这些残基的CA原子参与叠合。

陷阱三：RMSD对离群点敏感。由于是平方和，一个偏差5 Å的原子对RMSD的贡献是偏差1 Å原子的25倍。如果结构中有柔性loop或构象变化区域，全局RMSD会被少数大偏差原子主导。此时应改用截断RMSD或只计算核心区域。

陷阱四：单位与原子选择。RMSD单位是Å，但计算时若坐标单位是nm，结果会差10倍。另外，全原子RMSD和CA-only RMSD数值差异很大，报告时必须注明原子选择范围。

总结来说，刚体变换是“摆”的数学工具，旋转矩阵和四元数是它的两种参数化方式，RMSD是“摆好之后”的评分标准。理解三者的分工与边界，比记住公式更重要。工程实现中，反射检查、原子对应、离群点处理和单位一致性，是四个必须主动防守的环节。