📚 课程目录第3章坐标变换基础

坐标变换基础：从原子坐标到分子叠合

做3D构象比对，本质上是跟坐标打交道。你拿来两个分子，一个叫参考构象，一个叫待比对构象，怎么判断它们像不像？说白了就是：能不能通过旋转、平移，把一个分子上的原子挪到另一个分子对应原子上。

这就是刚体变换要做的事。不改变原子间的键长、键角，只做整体旋转和平移。我刚开始做分子叠合的时候，觉得这不就是中学学过的坐标加减吗？后来一上手，才发现坑很多。

平移与旋转的数学表达

先解决最简单的问题：把一个原子从点A挪到点B。

假设原始坐标是个三维向量 x = (x, y, z)^T。平移就是加一个向量 t = (t_x, t_y, t_z)^T：

x' = x + t

这个太简单，直接写代码就行。但我建议你把所有变换都用矩阵形式统一起来，后面会省很多事。

平移在齐次坐标下可以写成4×4矩阵：

[ x' ]   [ 1 0 0 tx ] [ x ]
[ y' ] = [ 0 1 0 ty ] [ y ]
[ z' ]   [ 0 0 1 tz ] [ z ]
[ 1  ]   [ 0 0 0 1  ] [ 1 ]

齐次坐标的核心价值：把平移塞进矩阵里，这样平移和旋转就可以串联成一个矩阵了。你在做迭代叠合时，只需维护一个4×4变换矩阵，而不是分开处理平移向量和旋转矩阵。

旋转稍微复杂一点。绕不同轴的旋转，矩阵不一样。

绕X轴旋转θ角：

R_x(θ) = [ 1      0       0    ]
         [ 0   cosθ   -sinθ   ]
         [ 0   sinθ    cosθ   ]

绕Y轴旋转θ角：

R_y(θ) = [  cosθ   0   sinθ ]
         [   0     1    0   ]
         [ -sinθ   0   cosθ ]

绕Z轴旋转θ角：

R_z(θ) = [ cosθ  -sinθ  0 ]
         [ sinθ   cosθ  0 ]
         [   0     0    1 ]

组合旋转就是矩阵乘法。注意顺序很重要——矩阵乘法不满足交换律。

R = R_z(γ) · R_y(β) · R_x(α)

这个顺序对应Z-Y-X欧拉角。先绕X转α，再绕新Y转β，最后绕新Z转γ。我在做蛋白质骨架叠合时，经常有人搞错顺序，叠合结果就跑偏了。

欧拉角与四元数

说到欧拉角，我想提一个著名的坑：万向锁（Gimbal Lock）。

当中间旋转角β接近90°时，第一个和第三个旋转轴会共线，你会丢失一个自由度。我在一次药物分子叠合项目中遇到过这个情况——连续优化了三天，叠合分数一直下不去。最后发现是欧拉角参数化导致优化过程中旋转矩阵退化。

一旦出现万向锁，分子就转不动了。叠合误差卡死在0.6Å，怎么都降不下来。

解决办法：用四元数。

四元数是推广的复数，形式是：

q = w + xi + yj + zk

其中w是实部，(x, y, z)是虚部。满足i² = j² = k² = ijk = -1。

一个单位四元数可以表示三维空间中的旋转。绕单位向量 u = (u_x, u_y, u_z) 转θ角，对应的四元数是：

q = cos(θ/2) + sin(θ/2)(u_x·i + u_y·j + u_z·k)

用四元数旋转点 p 的步骤：

把p写成纯四元数 p' = (0, x, y, z)
计算 p'_rot = q · p' · q⁻¹
从p'_rot的虚部取出旋转后的坐标

四元数没有万向锁，插值平滑，在分子构象搜索中超级好用。我建议你用四元数做构象叠合，用欧拉角做理解。

参数化方式	参数个数	有无万向锁	插值平滑性	适用场景
旋转矩阵	9	无（但有冗余约束）	差	理论推导、正变换
欧拉角	3	有	中等	可视化、简单旋转
四元数	4	无	好（球面插值）	优化叠合、构象搜索

用四元数时一定要归一化！浮点运算累积误差会让四元数偏离单位球面，导致旋转方向失真。我曾经在一次批处理叠合中发现中间结果异常，查了三天才发现是四元数没归一化。

归一化操作很简单：

norm = sqrt(w*w + x*x + y*y + z*z)
q = q / norm   # 每个分量除以norm

刚体变换在分子叠合中的应用

回到我们最关心的问题：给定两组对应原子坐标，怎么找到最优旋转和平移？

目标是让一组原子叠到另一组上，均方根偏差（RMSD）最小。数学表述是：

minimize_over_R,t   (1/N) * Σ || R*x_i + t - y_i ||²

其中R是旋转矩阵（3x3正交，det=1），t是平移向量（3x1）。

这个问题有解析解法，叫Kabsch算法。大致步骤：

质心对齐：两组坐标都减去各自质心
计算协方差矩阵H = X^TY
对H做SVD分解：H = UΣV^T
最优旋转矩阵R = VU^T（注意处理反射情况）
最优平移向量t = μ_y - R·μ_x

你用Python可以这样实现：

import numpy as np

def kabsch_align(X, Y):
    # X, Y: Nx3 矩阵，对应原子坐标
    # 步骤1：质心对齐
    centroid_X = np.mean(X, axis=0)
    centroid_Y = np.mean(Y, axis=0)
    X_centered = X - centroid_X
    Y_centered = Y - centroid_Y
    
    # 步骤2：协方差矩阵
    H = X_centered.T @ Y_centered
    
    # 步骤3：SVD
    U, S, Vt = np.linalg.svd(H)
    V = Vt.T
    
    # 步骤4：旋转（确保右手系）
    R = V @ U.T
    if np.linalg.det(R) < 0:
        V[:, -1] *= -1
        R = V @ U.T
    
    # 步骤5：平移
    t = centroid_Y.T - R @ centroid_X.T
    
    return R, t

这个代码我在实际项目中反复用了上百次。有一个细节——计算质心前一定要确认两个矩阵的原子顺序完全对应。

Kabsch算法的核心是SVD。协方差矩阵H = X^TY，这里的X和Y都是质心对齐后的坐标。然后SVD把H分解成旋转、缩放、再旋转。

有一个特殊情况：当det(VU^T) = -1时，你得到了一个反射——一个镜像变换，这不是合法的刚体旋转。这时候需要把V的第三列取反，再重新计算R。

我个人习惯：在计算RMSD之前，先检查一下两个分子是不是同一个参考系。有些PDB结构坐标系是文献里的，质心不归零。这时候直接跑Kabsch就会出问题——平移项应该是负质心，不是零。

用上面的代码，你从两个PDB文件里读出原子坐标，构建N×3矩阵，然后调用kabsch_align。返回的R和t就是你要的刚体变换，可以直接用于对齐。

如果你自己写分子叠合工具，我建议用四元数做优化搜索，用Kabsch做精确求解。前者在构象空间中找方向，后者在匹配点上做精确对齐。两个配合使用，效果最好。