📚 课程目录第4章 TM-score计算方法

TM-score计算方法：从结构比对到分数计算

好，咱们进入正题。TM-score这个东西，说白了就是衡量两个蛋白质结构有多像的一个分数。但它的计算方式，跟传统的RMSD还真不太一样。我刚开始接触的时候也觉得有点绕，咱们今天把它彻底理清楚。

一句话总结：TM-score先做结构对齐，再算加权距离，最后归一化成0到1之间的分数。0.5以上就算“同一个折叠”，0.2以下基本就算“没啥关系”。

第一步：结构比对——找到最优叠合

计算TM-score，首先得把两个结构“叠”在一起。这个步骤叫结构比对。我个人习惯把它拆成两部分：旋转平移和迭代优化。

你为什么需要做这一步？因为两个结构的坐标往往是不同的坐标系。一个蛋白可能在PBD文件里原点位置，另一个可能在(100, 100, 100)处。不先对齐，直接算距离根本没法比。

标准的做法是这样的：

选一组初始等价残基对（通常按序列顺序或者二级结构匹配）
做最小二乘旋转平移，让这组残基的坐标叠得最近
根据叠合后的距离重新选等价残基对
重复第2、3步，直到收敛

嗯，这里要注意。这套流程跟ABI-Net、SPRING、DeepAlign这些工具内部用的思路是一样的。迭代过程一般5-10次就够了，再多就是浪费。

第二步：坐标旋转与平移——数学细节

这一步的数学其实没那么可怕。给你看看实际代码会怎么写。

假设我们有两个结构，每个都有N个Cα原子坐标。我们用矩阵R和向量t来表示旋转和平移：

import numpy as np

def superimpose(coords1, coords2):
    """
    用Kabsch算法计算最优旋转平移
    coords1: (N, 3) 参考结构坐标
    coords2: (N, 3) 待对齐结构坐标
    返回旋转矩阵R，平移向量t
    """
    # 去质心
    c1 = np.mean(coords1, axis=0)
    c2 = np.mean(coords2, axis=0)
    coords1_centered = coords1 - c1
    coords2_centered = coords2 - c2
    
    # 协方差矩阵
    H = coords2_centered.T @ coords1_centered
    
    # SVD分解
    U, S, Vt = np.linalg.svd(H)
    R = Vt.T @ U.T
    
    # 检查行列式，确保是旋转不是反射
    if np.linalg.det(R) < 0:
        Vt[-1,:] *= -1
        R = Vt.T @ U.T
    
    t = c1 - R @ c2
    return R, t

我以前在项目里遇到过一个坑：有时候协方差矩阵是奇异的，SVD分解会出问题。解决方案很简单——加一个小正则项。不过那是另一回事了，今天先不展开。

第三步：TM-score公式——核心中的核心

对齐之后，我们得到了每个残基对的距离d_i。接下来就是TM-score公式上场了。

TMscore = (1 / N) * Σ [ 1 / (1 + (d_i / d0)^2) ]

其中 d0 = 1.24 * cube_root(N - 15) - 1.8
N 是比对上的残基总数
d_i 是第i对残基对齐后的欧氏距离

这个公式里，d0是个关键参数。它控制了距离的“宽容度”。为什么要有d0？你想想看，对于大蛋白，RMSD天然就比小蛋白大。如果不用d0做归一化，大蛋白和小蛋白的分数没法直接比。d0的设计让它随蛋白大小而调整——蛋白越大，d0也越大，对远端残基的惩罚就相对小一些。

我的经验：当N < 15时，直接用d0 = 1.5。这是结构生物学界的惯例。TM-score工具包里也是这么处理的，别自己去瞎改。

第四步：计算流程——一张图说清楚

这个流程看起来有点长，我给你画张图总结一下。

第五步：分数含义——多少算好？

分数算出来了，怎么解读？我把常见情况整理成一个表格，你对着看就行。

TM-score范围	结构相似性	实际含义
0.8 - 1.0	极高相似	基本上同一个结构，只差些侧链细节
0.5 - 0.8	中等相似	同一个折叠类型，核心区域对得不错
0.2 - 0.5	低相似	局部可能有些像，整体不是一回事
< 0.2	几乎不相似	随机对齐的水平，没有生物学意义

注意：TM-score的值跟比对长度非常敏感。如果你只比对了一小段（比如50个残基），而两个蛋白实际全长500个残基，那这个分数没有代表性。我见过很多人拿局部比对跑到0.6就高兴得不得了，结果全长一算只有0.1。切记，全长比对才可靠。

实际计算的代码示例

给你一个完整的计算示例。这段代码我在多个项目里都用过，改改路径就能跑。

from Bio.PDB import PDBParser
from Bio.PDB.Superimposer import Superimposer
import numpy as np

def compute_tm_score(pdb_file1, pdb_file2, chain_id='A'):
    """
    计算两个蛋白的TM-score
    """
    parser = PDBParser(QUIET=True)
    
    # 提取Cα原子坐标
    def get_ca_coords(pdb_file):
        structure = parser.get_structure('prot', pdb_file)
        coords = []
        for model in structure:
            for chain in model[chain_id]:
                for residue in chain:
                    if residue.get_id()[0] == ' ':  # 标准残基
                        if 'CA' in residue:
                            coords.append(residue['CA'].get_vector().get_array())
        return np.array(coords)
    
    coords1 = get_ca_coords(pdb_file1)
    coords2 = get_ca_coords(pdb_file2)
    
    N = min(len(coords1), len(coords2))
    
    # Kabsch对齐
    sup = Superimposer()
    sup.set_atoms(coords1[:N], coords2[:N])
    sup.apply(coords2)
    
    # 计算距离
    distances = np.linalg.norm(coords1[:N] - coords2[:N], axis=1)
    
    # d0公式
    d0 = 1.24 * (N - 15) ** (1/3) - 1.8
    if N < 15:
        d0 = 1.5
    
    # TM-score公式
    tm_score = np.mean(1 / (1 + (distances / d0) ** 2))
    
    return tm_score

这段代码跑起来很快，小蛋白秒出。但要注意——它假设两个蛋白的残基数量一样，而且序列已经对齐。实际场景中，你还需要额外处理不等长情况和序列对齐。TM-align这些工具做了一套很聪明的启发式方法，我们后面会讲。

避坑指南

我这些年做结构比对，踩过不少坑。挑几个最常见的跟你说说。

比对方向别搞反。TM-score不是对称的？其实不对称。拿结构A对齐到B，和B对齐到A，结果可能不一样。我建议总是用短链对齐到长链，或者取平均值。
d0别写死。有些人图省事，直接d0=3.0。对于小蛋白还行，大蛋白就偏差很大了。老老实实按公式算。
输出精度。TM-score通常保留3位小数。0.523就是0.523，别写0.52更别写0.5，丢掉分辨力了。

一个重要提醒：千万不要把你的预测结构和实验结构做RMSD比较时，直接用TM-score公式替代。它们是不同维度的评价指标。RMSD对全局刚体位移敏感，TM-score对局部偏差更宽容。两个都看，才能全面了解。

嗯，TM-score的计算方法就这些。核心公式很简单，但实现过程中的细节——迭代比对、d0计算、长度处理——才是决定分数质量的关键。下一块我们讲怎么在实际项目中批量跑TM-score，以及结果的统计分析。今天先消化这些。