3. 结构比对基础与工程实现

在蛋白结构预测中，模型输出的是一个三维坐标集合，而“预测得准不准”必须通过与实验结构（或参考结构）比对来回答。结构比对的核心任务有两层：第一，找到模型与参考结构之间合理的残基对应关系（序列对齐）；第二，在对应关系确定后，通过空间叠加使两组坐标尽可能重合，并用数值指标量化重合程度（结构叠加与打分）。不同指标对“误差”的定义不同，因此适用场景也不同。

一、四个指标的定义与计算原理

RMSD（均方根偏差）是最直观的指标。给定 N 对对应原子，先通过 Kabsch 算法求最优旋转和平移，使叠加后偏差平方和最小，再计算：

RMSD = sqrt( (1/N) * Σ ||x_i - y_i||² )

其中 x_i、y_i 是第 i 对对应原子的坐标。RMSD 的单位是 Å，数值越小越好。它的致命弱点是“全局平均”：一个局部严重错误的区域会被大量正确残基稀释，而且对结构域相对运动、柔性 loop 极其敏感。因此 RMSD 更适合比较同一蛋白的相近构象、或已裁剪到共同核心区的结构。

GDT-TS（Global Distance Test - Total Score）换了一种思路：不追求所有残基都对齐，而是统计“在某个距离阈值内能被对齐的残基比例”。常用阈值为 1、2、4、8 Å，对每个阈值分别计算在最优叠加下距离小于阈值的残基占比，再取平均：

GDT-TS = (1/4) * ( f(1Å) + f(2Å) + f(4Å) + f(8Å) )

其中 f(d) 表示在距离阈值 d 内可对齐的残基比例。GDT-TS 取值 0～1，越大越好。它对局部错误更宽容，能反映“大部分结构是否对”，因此长期作为 CASP 等结构预测评估的核心指标之一。

TM-score（Template Modeling score）在 GDT 思想基础上引入距离权重和长度归一化：

TM-score = max[ (1/L_target) * Σ 1/(1 + (d_i/d0)²) ]

其中 L_target 是参考结构长度，d_i 是第 i 对残基在最优叠加后的距离，d0 是长度相关的尺度因子（约 d0 = 1.24 * (L_target-15)^(1/3) - 1.8）。TM-score 取值 0～1，通常认为 >0.5 表示拓扑大致正确，>0.7 表示接近实验精度。由于做了长度归一化，TM-score 对长短蛋白可比性更好，且对局部大偏差不敏感。

lDDT（local Distance Difference Test）不依赖全局叠加，而是逐个残基检查“局部距离是否保持”。对参考结构中每一对距离小于包含半径（默认 15 Å）的原子对，计算模型中的对应距离与参考距离之差，统计差值小于 0.5、1、2、4 Å 的比例并取平均。lDDT 取值 0～1，越大越好。它不需要先做全局叠加，因此对结构域相对运动、多链装配、柔性区域更稳健，特别适合评估局部质量和复合物界面。

二、指标对比与适用场景

工程上常见的做法是组合使用：用 TM-score 判断整体拓扑是否可信，用 lDDT 定位局部低质量区域，用 RMSD 在裁剪后的核心区做精细比较，用 GDT-TS 与文献或 CASP 结果对齐口径。

三、序列对齐与结构叠加策略

结构比对的第一步是确定残基对应关系。若模型与参考序列完全相同，直接按残基编号一一对应即可；若存在缺失、插入或突变，则需要先做序列对齐（如 Needleman-Wunsch 全局对齐或 Smith-Waterman 局部对齐），只保留对齐上的残基对参与后续计算。对于多域蛋白，建议先按域分别对齐和叠加，再合并统计，避免域间相对运动污染全局指标。

结构叠加的数学核心是 Kabsch 算法：给定两组对应坐标，先各自去质心，计算协方差矩阵 H = XᵀY，对 H 做 SVD 分解 H = UΣVᵀ，则最优旋转 R = V·diag(1,1,det(VUᵀ))·Uᵀ，平移向量为两组质心之差。得到 R 和 t 后即可把模型坐标变换到参考坐标系，再计算各指标。

四、代码实现示例

下面给出一个基于 NumPy 的最小实现，覆盖 Kabsch 叠加、RMSD、TM-score 与 lDDT 的核心计算逻辑。示例中 P 为模型坐标，Q 为参考坐标，均为 (N,3) 数组，且已按序列对齐确定对应关系。

import numpy as np

def kabsch(P, Q):
    # P, Q: (N,3) 已对应的坐标
    Pc = P - P.mean(0)
    Qc = Q - Q.mean(0)
    H = Pc.T @ Qc
    U, S, Vt = np.linalg.svd(H)
    d = np.sign(np.linalg.det(Vt.T @ U.T))
    D = np.diag([1, 1, d])
    R = Vt.T @ D @ U.T
    t = Q.mean(0) - R @ P.mean(0)
    return R, t

def rmsd(P, Q):
    R, t = kabsch(P, Q)
    P_aligned = (R @ P.T).T + t
    return np.sqrt(((P_aligned - Q) ** 2).sum(1).mean())

def tm_score(P, Q):
    L = len(Q)
    d0 = 1.24 * (L - 15) ** (1/3) - 1.8
    R, t = kabsch(P, Q)
    P_aligned = (R @ P.T).T + t
    d = np.linalg.norm(P_aligned - Q, axis=1)
    return (1.0 / (1 + (d / d0) ** 2)).sum() / L

def lddt(P, Q, cutoff=15.0, thresholds=(0.5, 1, 2, 4)):
    # 基于参考结构 Q 的局部距离对
    Dq = np.linalg.norm(Q[:, None] - Q[None, :], axis=-1)
    Dp = np.linalg.norm(P[:, None] - P[None, :], axis=-1)
    mask = (Dq < cutoff) & (~np.eye(len(Q), dtype=bool))
    diff = np.abs(Dp[mask] - Dq[mask])
    scores = [(diff < t).mean() for t in thresholds]
    return float(np.mean(scores))

实际工程中，建议直接使用成熟库：TM-score 官方程序、US-align、OpenStructure、lDDT 的 lddt 模块、Bio.PDB 的 Superimposer 等。自实现主要用于理解原理和做轻量集成。

五、结果判读与常见陷阱

判读时不要只看单一数值。TM-score 0.5 附近往往意味着拓扑大体正确但细节偏差明显；lDDT 0.7 以上通常表示局部几何可信；RMSD 在核心区小于 2 Å 可视为高精度。若出现“RMSD 很小但 lDDT 很低”，通常说明全局叠加被少数刚性区域主导，而局部实际存在较大偏差，此时应检查是否裁剪了共同核心区、是否按域分别评估。

把上述流程固化为脚本后，每次拿到预测模型只需三步：准备参考结构与模型坐标、按序列对齐确定残基对、调用指标函数输出报告。建议在报告中同时记录参与计算的残基数量、是否裁剪核心区、是否按域分别评估，这些元信息往往比单一分数更能说明模型的实际可用性。